2.1. Dynamic Bayesian networks

PGMs represent joint probability distributions through efficient factorization, leveraging the conditional independence among random variables induced by the graph structure. Bayesian networks are a specific type of PGMs represented by directed acyclic graphs. A directed graph $ \mathcal{G}=\left\langle \mathcal{V},\mathrm{\mathcal{E}}\right\rangle $ is a graph composed of a set of vertices $ \mathcal{V} $ connected by directed edges, in $ \mathrm{\mathcal{E}}\subseteq \left\{\mathcal{V}\times \mathcal{V}\right\} $ which is a set of ordered pairs of vertices.

In a Bayesian network, $ \mathcal{V} $ represents a set of random variables, and the edges connecting the nodes represent causal relationships through conditional probability distributions (CPDs). Furthermore, given a set $ \mathcal{X}=\left\{{X}_1,\dots, {X}_N\right\} $ of $ N $ random variables (whether observed or hidden), we can factorize their joint probability distribution using the chain rule of probability, effectively leveraging the conditional independence structure induced by the graph topology. To this end, each node $ {X}_i $ , $ i=1,\dots, N $ , has an associated factor $ {\phi}^i=p\left({X}_i|{\mathrm{Pa}}_{X_i}\right) $ , which quantifies the influence of its parent nodes $ {\mathrm{Pa}}_{X_i} $ , reflecting the CPDs implied by the joint distribution. Nodes with no parents represent independent variables, specified by a prior (unconditional) distribution.

To extend Bayesian networks to settings where variables evolve over time, we denote the template variable $ {X}_i\in \mathcal{X} $ , $ i=1,\dots, N $ , at a specific point in time $ t $ as $ {X}_i^{(t)} $ . The trajectories of the (dynamic) Bayesian network refer to the complete set of possible realizations over a time horizon $ \left[0,T\right] $ , with $ t\in \left[0,T\right] $ . The joint distribution for all possible trajectories reads as follows:

(2.1) $$ p\left({\mathbf{X}}^{\left({t}_1:{t}_2\right)}\right)=p\left({\mathbf{X}}^{\left({t}_1\right)},\dots, {\mathbf{X}}^{\left({t}_2\right)}\right), $$

where $ \mathbf{X}={\left({X}_1,\dots, {X}_N\right)}^{\top } $ , and $ {\mathbf{X}}^{\left({t}_1:{t}_2\right)} $ denotes the set of variables $ \left\{{\mathbf{X}}^{(t)}:t\in \left[{t}_1,{t}_2\right]\right\} $ , $ {t}_1<{t}_2 $ .

To alleviate the computational burden associated with representing the joint distribution in Equation (2.1), a DBN typically relies on a set of simplifying assumptions. First, we consider a time discretization, which involves partitioning the time horizon $ \left[0,T\right] $ into time slices $ t=0,\dots, T $ , with an arbitrary step $ \Delta t $ . In this work, we assume regularly spaced time steps with granularity equal to the time between two consecutive observations of the system. Second, we assume that the conditional dependencies between consecutive random variables are represented by a Markov process. Accordingly, each variable $ {X}_i^{(t)} $ can have parents only in the current or immediately preceding time slice. This conditional independence structure allows us to compactly write $ p\left({\mathbf{X}}^{\left(0:T\right)}\right) $ as

(2.2) $$ p\left({\mathbf{X}}^{\left(0:T\right)}\right)=p\left({\mathbf{X}}^{(0)}\right)\prod \limits_{t=0}^{T-1}p\left({\mathbf{X}}^{\left(t+1\right)}|{\mathbf{X}}^{(t)}\right). $$

Hence, the probability distribution over different trajectories is fully characterized by the initial state distribution $ p\left({\mathbf{X}}^{(0)}\right) $ and the transition model $ p\left({\mathbf{X}}^{\left(t+1\right)}|{\mathbf{X}}^{(t)}\right) $ . This implies that probabilities at time $ t+1 $ are conditionally independent of all states prior to time $ t $ , given the state at time $ t $ .

Typically, a third assumption concerns the stationarity of the underlying Markov process, which allows the conditional dependencies between consecutive random variables to be represented by an invariant transition model:

(2.3) $$ p\left({\mathbf{X}}^{\left(t+1\right)}|{\mathbf{X}}^{(t)}\right)=p\left({\mathbf{X}}^{\prime }|\mathbf{X}\right)\hskip1em \forall t,\hskip1em \mathbf{X},{\mathbf{X}}^{\prime}\subseteq \mathcal{X}. $$

As detailed in the following, in this work, we relax this assumption, resulting in adaptive DTs capable of tracking systems that evolve according to a non-stationary process.

The transition model in Equation (2.3) can be represented by a two-time-slice Bayesian network (2-TBN) for a process over $ \mathcal{X} $ , which encodes the following CPD:

(2.4) $$ p\left({\mathbf{X}}^{\prime }|\mathbf{X}\right)=\prod \limits_{i=1}^Np\left({X}_i^{\prime }|{\mathrm{Pa}}_{X_x^{\prime }}\right), $$

where, in addition to $ {\mathbf{X}}^{\prime } $ , the 2-TBN includes only those variables in $ \mathbf{X} $ with a direct effect on $ {\mathbf{X}}^{\prime } $ .

Given a Bayesian network $ {\mathrm{B}}_0 $ representing the initial distribution $ p\left({\mathbf{X}}^{(0)}\right) $ , a DBN over the set of template variables $ \mathcal{X} $ is defined by a pair $ \left\langle {\mathrm{B}}_0,{\mathrm{B}}_{\to}\right\rangle $ , where $ {\mathrm{B}}_{\to } $ is a 2-TBN for the process over $ \mathcal{X} $ . The graph topology and conditional dependencies for $ t=0 $ follow those in $ {\mathrm{B}}_0 $ , while for $ t>0 $ they follow those in $ {\mathrm{B}}_{\to } $ . The resulting DBN provides a single compact model that captures the system dynamics and yields distributions over different trajectories. It can be interpreted as a sequence of Bayesian networks, one for each time slice.

For any desired time horizon $ T>0 $ , the distribution over $ {\mathbf{X}}^{\left(0:T\right)} $ is obtained by unrolling the DBN until $ T $ , repeating the core structure of the graph at each time step. Given a sequence of observations, a DBN can be unrolled until it fully encompasses the observed data. Additional time slices beyond the last observation do not affect inference within the observation period. Once unrolled, the DBN becomes a standard Bayesian network, and beliefs can be updated and propagated as for static Bayesian networks.

2.2. Probabilistic graphical model for predictive structural digital twins

The PGM used to represent the dynamic interaction between the physical and virtual domains is depicted in Figure 2 as a dynamic decision network, that is, a DBN with decision nodes. In the graphical representation, circle nodes denote random variables, square nodes indicate actionable inputs, and diamond-shaped nodes symbolize the objective function, all at discrete time steps. Each time the DT is updated by assimilating observational data, the PGM advances of one time step, with $ t\in \left\{0,\dots, T\right\} $ . Here, $ t=0 $ marks the time at which the DT enters operation, and $ t=T $ represents its lifetime horizon. Nodes with bold outlines represent observed quantities, while nodes with thin outlines represent unobserved quantities to be inferred. The PGM is sparsely connected, with edges encoding known or assumed conditional dependencies between variables.

Figure 2.

Dynamic Bayesian network encoding the asset–twin dynamical system, unrolled over two time slices. Circle nodes denote random variables, square nodes denote actions, and diamond nodes denote the objective function. Nodes with bold outlines denote observed quantities, while nodes with thin outlines denote unobserved quantities to be inferred. Directed edges represent the variables’ dependencies encoded via conditional probability distributions.

Figure 2. Long description

The flowchart is organized into three vertical columns corresponding to time steps t equals 0, t equals 1, and t equals 2.

At t equals 0:

* Top: Blue circle S sub 0 (thin outline) points to grey circle O sub 0 and blue circle S sub 1.

* Middle: Grey circle O sub 0 (bold outline) points to green circle D sub 0.

* Middle-Right: Red square U sub 0 (bold outline) points to S sub 1, D sub 1, and yellow diamond R sub 0.

* Bottom: Green circle D sub 0 (bold outline) points to U sub 0, R sub 0, and D sub 1.

At t equals 1:

* Top: Blue circle S sub 1 (thin outline) points to grey circle O sub 1 and blue circle S sub 2.

* Middle: Grey circle O sub 1 (bold outline) points to green circle D sub 1.

* Middle-Right: Red square U sub 1 (bold outline) points to S sub 2, D sub 2, and yellow diamond R sub 1.

* Bottom: Green circle D sub 1 (bold outline) points to U sub 1, R sub 1, and D sub 2.

At t equals 2:

* Top: Blue circle S sub 2 (thin outline) points to grey circle O sub 2 and continues right.

* Middle: Grey circle O sub 2 (bold outline) points to green circle D sub 2.

* Middle-Right: Red square U sub 2 (bold outline) points to yellow diamond R sub 2 and continues right.

* Bottom: Green circle D sub 2 (bold outline) points to U sub 2 and R sub 2.

A horizontal timeline at the bottom marks the intervals for t equals 0, t equals 1, and t equals 2.

We use capital letters to denote random variables representing the quantities involved in our abstraction, with subscripts indicating the time index. Corresponding lowercase letters represent realizations of these random variables. Calligraphic script variables denote the space of possible values that these quantities can take. For instance, the unknown physical state, which reflects the ground truth variability of interest, is denoted as $ {S}_t\sim p\left({s}_t\right) $ , where $ {s}_t $ represents the realization of the random variable $ {S}_t $ at time $ t $ , and $ p\left({s}_t\right) $ represents the probability that $ {S}_t={s}_t $ .

The digital state $ {D}_t\sim p\left({d}_t\right) $ aims to capture the variability in the physical state. It may include various types of information, such as initial conditions, boundary conditions, material properties, and other relevant quantities describing the asset (or process) of interest. However, a DT does not need to be an exact virtual replica to provide maximum value; rather, it should capture the essential features and variations relevant for diagnosis, prediction, and decision-making (Ferrari and Willcox, Reference Ferrari and Willcox2024). Although there are no strict constraints in this regard, we define the digital state as a categorical random variable that can assume values among a set of predefined configurations in the digital state space $ \mathcal{D} $ .

The physical-to-digital information flow from $ {S}_t $ to $ {D}_t $ is governed by the continuous collection and assimilation of observational data $ {O}_t\sim p\left({o}_t\right) $ . These may include sensor recordings, inspection data, or diagnostic information. Since the physical state $ {S}_t $ is only partially and indirectly observable, the digital state $ {D}_t $ cannot be uniquely determined. Instead, it provides evidence for digital state configurations that could be consistent with the given observations.

The updated digital state $ {D}_t $ informs the digital-to-physical information flow, that is, the control actions. In Figure 2, $ {U}_t\sim p\left({u}_t\right) $ is a decision variable representing the action to be taken. The space of actions $ \mathcal{U} $ may include interventions affecting the physical state, adjustments to the operational conditions, or modifications to the observational process. Each action is characterized by its own transition dynamics across the system’s state space.

Finally, the reward node $ {R}_t\sim p\left({r}_t\right) $ measures the performance of the asset–twin system. Specifically, the reward function provides values within a reward space $ \mathcal{R} $ , measuring the expected utility of DBN trajectories that we maximize when solving the planning problem induced by the PGM.

Since we assume discrete spaces for digital states and control actions, the related causal relationships are modeled using conditional probability tables (CPTs). In a CPT, each column contains the CPD of a node’s value given the values of its parent nodes. Each column sums to one, as the entries represent an exhaustive set of possible outcomes for the variable. A node with no parents has only one column, indicating the prior probabilities of each possible value. CPTs can be viewed as a set of parameters representing prior knowledge about random variables. Thanks to the discrete DBN setting, we can propagate and update beliefs exactly with a single pass of the sum-product message-passing algorithm (Pearl, Reference Pearl, Geffner, Dechter and Halpern1982).

The Bayesian network $ {\mathrm{B}}_0 $ and the 2-TBN $ {\mathrm{B}}_{\to } $ , jointly forming our PGM, are shown in Figure 3. $ {\mathrm{B}}_0 $ encodes the initial distribution of the random variables, which can be derived from expert knowledge to assign a prior distribution to the initial state of the asset–twin system. In contrast, $ {\mathrm{B}}_{\to } $ encodes the inter-slice conditional dependencies, modeling the temporal relationships between the random variables. The PGM enables the extraction of the joint distribution of the network’s trajectories, conditioned on observed data $ {O}_t={o}_t $ and enacted actions $ {U}_t={u}_t $ . Specifically, the joint belief state can be factorized from the initial time step $ t=0 $ up to the current time $ {t}_c $ , using the following sequential Bayesian inference formulation:

(2.5) $$ p\left({D}_{0:{t}_c},{U}_{0:{t}_c},{R}_{0:{t}_c}|{o}_{0:{t}_c},{u}_{0:{t}_c}\right)\propto \prod \limits_{t=0}^{t_c}\left[{\phi}_t^{\mathrm{update}}{\phi}_t^{\mathrm{control}}{\phi}_t^{\mathrm{reward}}\right], $$

Figure 3.

Dynamic Bayesian network encoding the asset–twin dynamical system: Bayesian network representing the initial distribution $ {\mathcal{B}}_0 $ ; two-time-slice Bayesian network for the process $ {\mathcal{B}}_{\to } $ . Circle nodes denote random variables, square nodes denote actions, and diamond nodes denote the objective function. Nodes with bold outlines denote observed quantities, while nodes with thin outlines denote unobserved quantities to be inferred. Directed edges represent the variables’ dependencies encoded via conditional probability distributions.

Figure 3. Long description

The diagram consists of two panels.

Left panel: Labeled B sub 0, representing the initial distribution.

* At the top left is a blue circle S sub 0 with a thin outline, pointing down to a grey circle O sub 0 with a bold outline.

* O sub 0 points down to a green circle D sub 0 with a thin outline.

* D sub 0 points up to a red square U sub 0 with a bold outline and right to a yellow diamond R sub 0.

* U sub 0 also points down to R sub 0.

Right panel: Labeled B sub arrow, representing the process across two time slices, t minus 1 and t, enclosed in dashed grey boxes.

* Time slice t minus 1 contains blue circle S sub t minus 1 and green circle D sub t minus 1. Between the two time slices is a red square U sub t minus 1 with a bold outline.

* Time slice t contains blue circle S sub t, grey circle O sub t with a bold outline, green circle D sub t, red square U sub t with a bold outline, and yellow diamond R sub t.

* Directed edges show dependencies: S sub t minus 1 points to S sub t. D sub t minus 1 points to D sub t. U sub t minus 1 points to both S sub t and D sub t.

* Within time slice t: S sub t points to O sub t. O sub t points to D sub t. D sub t points to both U sub t and R sub t. U sub t points to R sub t.

with factors:

(2.6) $$ {\phi}_t^{\mathrm{update}}=p\left({D}_t|{D}_{t-1},{U}_{t-1}={u}_{t-1},{O}_t={o}_t\right), $$

(2.7) $$ {\phi}_t^{\mathrm{control}}=p\left({U}_t|{D}_t\right), $$

(2.8) $$ {\phi}_t^{\mathrm{reward}}=p\left({R}_t|{D}_t,{U}_t={u}_t\right). $$

The factor $ {\phi}_t^{\mathrm{update}} $ updates the digital state belief $ {D}_t $ to reflect the (hidden) physical state $ {S}_t $ . First, it accounts for a control-dependent transition model, describing the expected evolution of the digital state given the previous belief $ {D}_{t-1} $ and the most recent action taken $ {U}_{t-1}={u}_{t-1} $ . Additionally, it integrates this prior with the data assimilation outcome by conditioning on $ {O}_t={o}_t $ . We can explicitly split the contributions of transition dynamics and data assimilation by further factorizing $ {\phi}_t^{\mathrm{update}} $ into a predictor–corrector form update, as follows:

(2.9) $$ {\phi}_t^{\mathrm{update}}\propto {\phi}_t^{\mathrm{data}}{\phi}_t^{\mathrm{history}}, $$

where:

(2.10) $$ {\phi}_t^{\mathrm{history}}=p\left({D}_t|{D}_{t-1},{U}_{t-1}={u}_{t-1}\right), $$

(2.11) $$ {\phi}_t^{\mathrm{data}}=p\left({O}_t={o}_t|{D}_t\right). $$

The causal structure required for this factorization is not included in the set of conditional dependencies induced by the graph depicted in Figure 2. Graphically, this would require an arrow from $ {D}_t $ to $ {O}_t $ , but in our model, the arrow is depicted in the opposite direction to reflect the information flow within the DT, consistent with MG Kapteyn et al., (Reference Kapteyn, Pretorius and Willcox2021) and Torzoni et al., (Reference Torzoni, Tezzele, Mariani, Manzoni and Willcox2024). For further details, the reader may refer to Koller and Friedman, (Reference Koller and Friedman2009). The transition model $ {\phi}_t^{\mathrm{history}}:\mathcal{D}\times \mathcal{D}\times \mathcal{U}\mapsto \left[0,1\right] $ acts as a forward-time predictor. It can be visualized as a CPT describing the expected transition dynamics from $ {D}_{t-1} $ to $ {D}_t $ for each possible action $ {U}_{t-1}={u}_{t-1} $ . This CPT should incorporate any prior knowledge the DT designer has regarding the asset and relevant operational conditions. It can be derived offline from historical data and accumulated experience (Thompson and Johnson, Reference Thompson and Johnson2005; Morcous, Reference Morcous2006; Bianchi and Biondini, Reference Bianchi and Biondini2022), as well as from physics-based or empirical evolution models (Biondini and Frangopol, Reference Biondini and Frangopol2016, Reference Biondini and Frangopol2017; Garayalde et al., Reference Garayalde, Rosafalco, Torzoni and Corigliano2025). In this paper, we introduce a stochastic parametrization of $ {\phi}_t^{\mathrm{history}} $ , similar to Tezzele et al., (Reference Tezzele, Carr, Topcu and Willcox2026), enabling online learning of state transition beliefs via Bayesian updates, as described in Section 3.1. In contrast, the factor $ {\phi}_t^{\mathrm{data}} $ is associated with the assimilation of observational data $ {O}_t={o}_t $ , as further described at the end of this section.

The factor $ {\phi}_t^{\mathrm{control}}:\mathcal{U}\times \mathcal{D}\mapsto \left[0,1\right] $ evaluates the belief about the next optimal control action. This factor corresponds to a CPT encoding a control policy $ \pi :\mathcal{D}\mapsto \mathcal{U} $ . The actual action to be taken is selected as the best point estimate $ {u}_t=\arg \hskip0.1em {\max}_{u\in \mathcal{U}}\pi \left({D}_t\right) $ , that is, the conversion from probabilistic to deterministic control. The computation of $ \pi $ is described in Section 3.2.

Finally, $ {\phi}^{\mathrm{reward}}:\mathcal{D}\times \mathcal{U}\to \mathrm{R} $ encapsulates the evaluation of the objective function. The two edges leading to $ {R}_t $ represent a reward function that quantifies the costs associated with the digital state and control actions, respectively.

Given the current digital state $ {D}_{t_c} $ , probabilistic estimates of future digital states and control actions are made by unrolling a simplified DBN up to a prediction time $ {t}_p>{t}_c $ . This prediction DBN, depicted in Figure 4, retains a structure similar to the one discussed above but is stripped of the nodes related to observations, data assimilation, and action execution. At the prediction stage, information from sensing and performed actions is no longer available. Accordingly, the actions $ {U}_t $ are depicted as random variables rather than decision nodes, as they represent predictions of what-if scenarios beyond data assimilation, made through uncertainty propagation via DBN unrolling. As a result, the factorization over the trajectories in Equation (2.5) can be extended over the prediction time steps $ t={t}_c+1,\dots, {t}_p $ by augmenting the belief state as follows:

(2.12) $$ p\left({D}_{t_c+1:{t}_p},{U}_{t_c+1:{t}_p},{R}_{t_c+1:{t}_p}|{D}_{t_c},{U}_{t_c}\right)=\prod \limits_{t={t}_c+1}^{t_p}\left[{\phi}_t^{\mathrm{history}}{\phi}_t^{\mathrm{control}}{\phi}_t^{\mathrm{reward}}\right]. $$

In this case, the factors $ {\phi}_t^{\mathrm{history}}=p\left({D}_t|{D}_{t-1},{U}_{t-1}\right) $ and $ {\phi}_t^{\mathrm{reward}}=p\left({R}_t|{D}_t,{U}_t\right) $ enable future state predictions and probabilistic reward estimates, respectively, conditioned on predicted control actions that have not yet been enacted. As clarified in Section 3.2, factors $ {\phi}_t^{\mathrm{history}} $ and $ {\phi}_t^{\mathrm{control}} $ actually encode dynamically updated transition models and control policies, respectively. Indeed, after updating the transition model $ {\phi}_t^{\mathrm{history}} $ with the evidence collected up to $ {t}_c $ , we also recompute an optimal policy $ {\pi}_{t_c} $ , such that $ {U}_{t_c:{t}_p}\sim {\pi}_{t_c}\left({D}_{t_c:{t}_p}\right) $ .

Figure 4.

Dynamic Bayesian network used to predict the future evolution of digital states and control actions. Circle nodes denote random variables, and diamond nodes denote the objective function. Directed edges represent the variables’ dependencies encoded via conditional probability distributions.

Figure 4. Long description

A horizontal timeline at the bottom marks three points: t equals t sub c plus 1, t equals t sub c plus 2, and t equals t sub p.

At time t equals t sub c plus 1, a green circle node labeled D sub t sub c plus 1 receives two incoming arrows from the left. From this node, arrows point right to the next state node, upward to a red circle node labeled U sub t sub c plus 1, and diagonally to a yellow diamond node labeled R sub t sub c plus 1. The red node U sub t sub c plus 1 also has a downward arrow pointing to the yellow diamond R sub t sub c plus 1 and a diagonal arrow pointing to the next state node.

At time t equals t sub c plus 2, the green circle node D sub t sub c plus 2 receives arrows from the previous state and control nodes. It repeats the internal structure: arrows point right to the final state, upward to red circle U sub t sub c plus 2, and diagonally to yellow diamond R sub t sub c plus 2. Red circle U sub t sub c plus 2 points down to the diamond and diagonally to the final state.

At the final time step t equals t sub p, the green circle node D sub t sub p receives the incoming arrows. It points upward to red circle U sub t sub p and diagonally to yellow diamond R sub t sub p. The red circle U sub t sub p points down to the yellow diamond R sub t sub p. An arrow from the final green node continues to the right toward an arrowhead on the timeline.

In this work, the asset of interest is a deteriorating structural system, although the framework outlined above can be adapted to other types of systems or applications beyond SHM. The asset is dynamic, meaning its state evolves over time according to laws that depend on both the system itself and external factors. These external factors can either improve or degrade the condition of the asset and influence its deterioration rate. Such factors may include long-term degradation caused by chemical, physical, or mechanical processes, as well as sudden changes due to discrete damage events or interventions by external agents (Zakić et al., Reference Zakić, Ryzynski, Guo-Hong and Jokela1991).

We simulate the monitored asset using a physics-based computational model. Specifically, the structure is modeled as a linear-elastic continuum, and its dynamic response to applied loads is described through elasto-dynamics equations, assuming linearized kinematics. The model is spatially discretized using finite elements, and its solution is advanced over time to produce synthetic observational data. This computational model serves as a data generator, controlled by a vector of $ {N}_{\mathrm{par}} $ parameters $ \boldsymbol{\mu} \in {\mathrm{\mathbb{R}}}^{N_{\mathrm{par}}} $ , which describe the operational, damage, and (possibly) environmental conditions. Further details on the modeling choices for the Hörnefors railway bridge case study considered in this work will be provided in Section 4.

The observational data consist of vibration recordings shaped as multivariate time series $ \mathbf{U}\left(\boldsymbol{\mu} \right)=\left[{\mathbf{u}}_1\left(\boldsymbol{\mu} \right),\dots, {\mathbf{u}}_{N_s}\left(\boldsymbol{\mu} \right)\right]\in {\mathrm{\mathbb{R}}}^{L\times {N}_s} $ , made up of $ {N}_s $ time series, each containing $ L $ sensor measurements equally spaced in time, such as accelerations or displacements. Accordingly, the observation space $ \mathcal{O} $ is described through the sensitivity of $ \mathbf{U} $ with respect to $ \boldsymbol{\mu} $ . For the problem settings we consider, each observation process refers to a time interval that is short enough for the operational, environmental, and damage conditions to be assumed time-invariant, yet long enough to allow accurate identification of the structural behavior.

As the unknown physical state reflects the ground truth variability in the structural health of the asset, the categorical random variable modeling the digital state describes the presence, location, and magnitude of damage among a set of predefined configurations in $ \mathcal{D} $ . We describe changes in the structural dynamic response due to damage through a localized reduction of effective stiffness. If the time duration of the observation process is small compared to the time required for damage progression, this approach allows the structural behavior to be treated as linear within an observation interval, thus enabling a timescale separation between damage evolution and health assessment (Eftekhar Azam and Mariani, Reference Eftekhar Azam and Mariani2018). Examples of damage patterns that can be represented in this way include the failure of bolted connections and corrosion in steel structures (Pandey and Biswas, Reference Pandey and Biswas1994; Mariniello et al., Reference Mariniello, Pastore, Menna, Festa and Asprone2021), cracking and delamination in composite structures (Kapteyn et al., Reference Kapteyn, Knezevic, Huynh, Tran and Willcox2020), cracking in reinforced concrete structures (Teughels et al., Reference Teughels, Maeck and De Roeck2002; Mariniello et al., Reference Mariniello, Pastore, Menna, Festa and Asprone2021), cracking and crushing-induced stiffness losses in masonry brick structures (Hernández-González and García-Macías, Reference Hernández-González and García-Macías2024), cuts in steel beams (Seventekidis et al., Reference Seventekidis, Giagopoulos, Arailopoulos and Markogiannaki2020), and the erosion of piles and abrasion of deck surfaces (Fernandez-Navamuel et al., Reference Fernandez-Navamuel, Zamora-Sánchez, Omella, Pardo, Garcia-Sanchez and Magalhães2022). Here, the local stiffness reduction shaping $ \mathcal{D} $ is achieved by means of two variables, $ y\in \mathrm{\mathbb{N}} $ and $ \delta \in \mathrm{\mathbb{R}} $ , both included in $ \boldsymbol{\mu} $ . Specifically, $ y\in \mathrm{\mathbb{N}} $ labels the specific damage region among a set of $ {N}_{\Omega}\in \mathrm{\mathbb{N}} $ predefined damageable regions, where $ y=0 $ denotes the undamaged baseline, while $ \delta \in \mathrm{\mathbb{R}} $ represents the magnitude of the stiffness reduction occurring in the region associated with $ y $ .

Data assimilation, represented by the $ {\phi}_t^{\mathrm{data}} $ factor, is performed using DL models to approximate the inverse mapping from measurements $ \mathbf{U} $ to the target parameters $ y $ and $ \delta $ . The deterministic output of the DL models is then converted into evidence for configurations in the digital state space that are consistent with the given observation, thereby serving as the likelihood mechanism informing the posterior. This conversion is carried out through a confusion matrix that measures the offline (expected) performance of the DL models in correctly identifying the digital state among the predefined configurations in $ \mathcal{D} $ . The confusion matrix can be interpreted as a CPT, where rows represent the ground truth digital state, and columns represent the predicted digital state. For a more detailed description of the evaluation of the $ {\phi}_t^{\mathrm{data}} $ , the reader is referred to (Torzoni et al., Reference Torzoni, Tezzele, Mariani, Manzoni and Willcox2024).

3.1. Parametric Markov decision processes

In a decision-making setting, an agent must choose from a set of possible actions, each of which may have uncertain effects on the system state. The goal of the decision-making process is to maximize the expected utilities associated with the outcomes of these actions. Sequential decision problems expand on this concept by requiring the agent to choose a sequence of actions that optimizes the long-term accumulation of rewards over the planning horizon.

MDPs provide a framework for modeling sequential decision problems in stochastic environments, using Markov transition models and additive rewards (Bellman, Reference Bellman1957b; Howard, Reference Howard1960). A discrete MDP is defined as a $ 4 $ -tuple $ \left\langle \mathcal{S},\mathcal{U},\mathcal{P},\mathrm{R}\right\rangle $ , consisting of a set of states $ \mathcal{S} $ , a set of actions $ \mathcal{U} $ , a CPT $ \mathcal{P}\in {\mathrm{\mathbb{R}}}^{\mathcal{D}\times \mathcal{D}\times \mathcal{U}} $ describing the Markov transition model, and a set of rewards $ \mathrm{R} $ .

The control action is selected as the best point estimate, thereby mapping the probabilistic state representation to a deterministic decision. This is consistent with the standard MDP assumption of full observability. While alternative approaches account for partial observability (e.g., (Arcieri et al., Reference Arcieri, Hoelzl, Schwery, Straub, Papakonstantinou and Chatzi2023, Reference Arcieri, Hoelzl, Schwery, Straub, Papakonstantinou and Chatzi2024; Papakonstantinou et al., Reference Papakonstantinou, Andriotis, Shinozuka, Bakker, Frangopol and Breugel2016), they require more complex and computationally demanding solution strategies. In the present work, we rely on the high accuracy of the deep learning models used for digital state classification, which provide a near-complete characterization of the system state. On this basis, we approximate the true state with its most likely estimate, effectively decoupling state estimation from control. We note, however, that this simplification neglects residual uncertainty and may lead to suboptimal decisions when classification uncertainty is non-negligible.

In our DT framework, the MDP associated with optimal control planning is encoded by the prediction DBN shown in Figure 4. While in (Torzoni et al., Reference Torzoni, Tezzele, Mariani, Manzoni and Willcox2024) the transition model $ {\phi}_t^{\mathrm{history}} $ was kept fixed, we relax the assumption of a stationary Markov process for the conditional dependencies between random variables in consecutive time slices. Specifically, we treat the transition probabilities in $ \mathcal{P} $ as random variables with conjugate distributions, similar to (Tezzele et al., Reference Tezzele, Carr, Topcu and Willcox2026). Starting with an initial prior, the transition dynamics is learned online using Bayesian updates that incorporate new evidence about the system response to actions. This approach is scalable and allows ( $ i $ ) adapting DTs starting from (possibly, inaccurate) initial priors or ( $ ii $ ) tracking non-stationary, evolving contexts.

Actions with uncertain effects produce non-deterministic (and possibly unexpected) outcomes. Within an SHM context, inter-state transitions can be described using (stochastic) transition matrices $ {\mathcal{P}}_u\in {\mathrm{\mathbb{R}}}^{\mathcal{D}\times \mathcal{D}} $ , one for each action $ u\in \mathcal{U} $ :

(3.1) $$ {\mathcal{P}}_u=\left(\begin{array}{cccc}p\left({d}_1|{d}_1,u\right)& p\left({d}_1|{d}_2,u\right)& \dots & p\left({d}_1|{d}_{\mid \mathcal{D}\mid },u\right)\\ {}p\left({d}_2|{d}_1,u\right)& p\left({d}_2|{d}_2,u\right)& \dots & p\left({d}_2|{d}_{\mid \mathcal{D}\mid },u\right)\\ {}\vdots & \vdots & \ddots & \vdots \\ {}p\left({d}_{\mid \mathcal{D}\mid }|{d}_1,u\right)& p\left({d}_{\mid \mathcal{D}\mid }|{d}_2,u\right)& \dots & p\left({d}_{\mid \mathcal{D}\mid }|{d}_{\mid \mathcal{D}\mid },u\right)\end{array}\right), $$

where $ \mid \mathcal{A}\mid \in \mathrm{\mathbb{N}} $ denotes the cardinality of a generic set $ \mathcal{A} $ . Matrix $ {\mathcal{P}}_u $ collects the probabilities $ p\left({d}^{\prime }|d,u\right) $ of transitioning from state $ d\in \mathcal{D} $ to state $ {d}^{\prime}\in \mathcal{D} $ , under a given action $ u\in \mathcal{U} $ . Here, columns represent the starting state $ d $ , and rows represent the arrival state $ {d}^{\prime } $ . Each column of $ {\mathcal{P}}_u $ sums to one, reflecting the exhaustive set of transition probabilities from a given state $ d $ under the action $ u $ . By partitioning $ \mathcal{D} $ in lexicographic order, with states ordered first by damage location and then by damage level, the diagonal entries in $ {\mathcal{P}}_u $ denote the probability that the system has of remaining in the same state, while the lower-left and upper-right triangles represent the probabilities of system deterioration and improvement, respectively. The CPT describing the complete transition model, tensorial with respect to all available actions, is given by $ \mathcal{P}=\left\{{\mathcal{P}}_1\dots, {\mathcal{P}}_{\mid \mathcal{U}\mid}\right\} $ .

The specific content of the transition matrices should account for the structural system under consideration, its operational conditions, and the surrounding environmental variability (Torzoni et al., Reference Torzoni, Rosafalco, Manzoni, Mariani and Corigliano2022). To address these sources of uncertainty, we adopt a hierarchical Dirichlet–Multinomial parametrization. For each state–action pair, the categorical transition distribution is endowed with a Dirichlet prior, whose hyperparameters encode prior beliefs over transition outcomes. This defines a hierarchical Bayesian structure in which transition probabilities are treated as latent random variables, and uncertainty is propagated through the higher-level Dirichlet prior, capturing both epistemic uncertainty and data-driven adaptation (Gelman et al., Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin2013). These parameters are then subject to temporal updating, that is, learned online as new data arrives, enabling adaptive Bayesian inference.

We model the transition probabilities $ p\left({d}^{\prime }|d,u\right) $ with categorical distributions parameterized by $ {\boldsymbol{\theta}}_{d,u}={\left({\theta}_{d,u}^1,\dots, {\theta}_{d,u}^{\mid \mathcal{D}\mid}\right)}^{\top}\in {\mathrm{\mathbb{R}}}^{\mid \mathcal{D}\mid } $ , where $ {\sum}_{k=1}^{\mid \mathcal{D}\mid }{\theta}_{d,u}^k=1 $ and $ k=1,\dots, \mid \mathcal{D}\mid $ indexes the arrival states $ {d}^{\prime } $ . The uncertainty over $ {\boldsymbol{\theta}}_{d,u} $ is represented by a Dirichlet prior:

(3.2) $$ {\boldsymbol{\theta}}_{d,u}\sim \mathrm{Dirichlet}\left({\boldsymbol{\alpha}}_{d,u}\right), $$

where $ {\boldsymbol{\alpha}}_{d,u}={\left({\alpha}_{d,u}^1,\dots, {\alpha}_{d,u}^{\mid \mathcal{D}\mid}\right)}^{\top}\in {\mathrm{\mathbb{R}}}^{\mid \mathcal{D}\mid } $ are (positive) concentration parameters encoding prior beliefs about the transition probabilities. We use the Dirichlet distribution because it is conjugate to the Multinomial likelihood, ensuring that the posterior distribution remains within the Dirichlet family after incorporating observations. This property allows efficient refinement of the transition probabilities and beliefs by updating the Dirichlet hyperpriors in closed form.

Managing $ \mid \mathcal{D}\mid \times \mid \mathcal{U}\mid $ Dirichlet distributions, each with $ \mid \mathcal{D}\mid $ parameters, can become computationally intractable for large state spaces. While we provide the formulation and update scheme for such a state-dependent transition model in A, a more computationally viable approach to handling multi-step transitions is to assign equal transition probabilities between consecutive elements of the partition, rather than a separate prior for each individual state transition (Ross et al., Reference Ross, Pineau, Chaib-draa and Kreitmann2011; Papakonstantinou et al., Reference Papakonstantinou, Andriotis, Shinozuka, Bakker, Frangopol and Breugel2016; Andriotis and Papakonstantinou, Reference Andriotis and Papakonstantinou2019; Tezzele et al., Reference Tezzele, Carr, Topcu and Willcox2026). In other words, we assign probabilities of remaining in the same state, of transitioning by one step, two steps, and so on, independently of the starting state $ d $ . This state-independent transition model is thus parametrized through $ {\mathbf{p}}_u={\left({p}_u^0,{p}_u^1,\dots, {p}_u^I\right)}^{\top}\in {\mathrm{\mathbb{R}}}^{I+1} $ , representing transition probabilities of $ 0,1,\dots, I $ steps forward under action $ u $ . The corresponding transition matrix $ {\mathcal{P}}_u $ is then constructed from $ {\mathbf{p}}_u $ , ensuring the same column structure for all starting states, except for absorbing states—that is, states from which no further transitions occur, meaning they retain probability $ 1 $ of remaining in the same state.

As an example, for a state space of cardinality $ \mid \mathcal{D}\mid =5 $ with a single damage location, the transition matrix for an action under which the structural health state can deteriorate by up to $ I=2 $ steps is:

(3.3) $$ {\mathcal{P}}_u\left({\mathbf{p}}_u\right)=\left(\begin{array}{ccccc}{p}_u^0& 0& 0& 0& 0\\ {}{p}_u^1& {p}_u^0& 0& 0& 0\\ {}{p}_u^2& {p}_u^1& {p}_u^0& 0& 0\\ {}0& {p}_u^2& {p}_u^1& {p}_u^0& 0\\ {}0& 0& {p}_u^2& {p}_u^1+{p}_u^2& 1\end{array}\right). $$

As per Equation (3.2), the probabilities $ {\mathbf{p}}_u $ are assigned a Dirichlet prior:

(3.4) $$ {\mathbf{p}}_u\sim \mathrm{Dirichlet}\left({\boldsymbol{\alpha}}_u\right), $$

where $ {\boldsymbol{\alpha}}_u={\left({\alpha}_u^0,{\alpha}_u^1,\dots, {\alpha}_u^I\right)}^{\top}\in {\mathrm{\mathbb{R}}}^{I+1} $ are the concentration parameters for the step transitions under action $ u $ .

Let $ {d}_{0:{t}_c} $ and $ {u}_{0:{t}_c-1} $ denote the history of system states and actions observed up to the current time $ {t}_c $ , respectively. Additionally, let $ {\mathbf{y}}_u={\left({N}_u^0,{N}_u^1,\dots, {N}_u^I\right)}^{\top}\in {\mathrm{\mathbb{N}}}^{I+1} $ represent the counts of observed step transitions under $ \left\{{d}_{0:{t}_c},{u}_{0:{t}_c-1}\right\} $ . The likelihood of $ {\mathbf{y}}_u $ , given the transition probabilities $ {\mathbf{p}}_u $ , is modeled as a Multinomial distribution:

(3.5) $$ {\mathbf{y}}_u\mid {\mathbf{p}}_u\sim \mathrm{Multinomial}\left({N}_u,{\mathbf{p}}_u\right), $$

where $ {N}_u={\sum}_{i=0}^I{N}_u^i $ is the total transition count observed under action $ u $ . Using Bayes’ theorem, the posterior distribution of $ {\mathbf{p}}_u $ after observing $ {\mathbf{y}}_u $ remains in the Dirichlet family due to conjugacy:

(3.6) $$ {\mathbf{p}}_u\mid {\mathbf{y}}_u\sim \mathrm{Dirichlet}\left({\overline{\boldsymbol{\alpha}}}_u\right), $$

with updated concentration parameters $ {\overline{\boldsymbol{\alpha}}}_u\in {\mathrm{\mathbb{R}}}^{I+1} $ , computed as:

(3.7) $$ {\overline{\boldsymbol{\alpha}}}_u={\left({\alpha}_u^0+{N}_u^0,{\alpha}_u^1+{N}_u^1,\dots, {\alpha}_u^I+{N}_u^I\right)}^{\top }. $$

The parametrized Dirichlet–Multinomial transition model efficiently accounts for the possibility of multi-step transitions. However, it does not address the possibility of inconsistencies arising from evidence that contradicts the transition model. This limitation can be particularly critical, as assigning a zero probability to potential state transitions could compromise the inference process (Lindley, Reference Lindley2013). To resolve this potential inconsistency, a perturbation matrix $ \Xi \in {\mathrm{\mathbb{R}}}^{\mid \mathcal{D}\mid \times \mid \mathcal{D}\mid } $ , with positive entries, arbitrarily small relative to the non-zero terms of $ {\mathcal{P}}_u $ , is introduced. This matrix is used to obtain a perturbed transition matrix as follows:

(3.8) $$ {\mathcal{P}}_u^{\Xi}\left({\mathbf{p}}_u\right)={\mathcal{P}}_u\left({\mathbf{p}}_u\right)+\Xi, $$

which is then properly normalized.

A special case of the state-independent Dirichlet–Multinomial transition model is the Beta-Binomial model, which is parameterized by a single one-step transition probability. The formal derivation follows an approach analogous to the one presented in this section and is provided in B for completeness. The numerical experiments in Section 4 use Dirichlet–Multinomial state-independent transition models.

3.2. Dynamic control policies and precision updates

The $ {\phi}_t^{\mathrm{control}} $ factor in Equation (2.5) represents a CPT encoding a health-aware control policy $ {\pi}_t\left({D}_t\right) $ . In (Torzoni et al., Reference Torzoni, Tezzele, Mariani, Manzoni and Willcox2024), a stationary policy was computed offline for an infinite planning horizon, based on the simplifying assumption of known transition dynamics. Here, the sequence of actions $ {u}_{t_c:T}^{\ast } $ defining the optimal control policy $ {\pi}^{\ast } $ is systematically recomputed using MBRL with a dynamically updated transition model. The policy is recomputed over the finite planning horizon $ \left[{t}_c,T\right] $ , representing the design lifetime of the asset. The corresponding optimization problem is expressed as:

(3.9) $$ {\pi}_t^{\ast}\left({D}_t\right)=\underset{\pi_t}{\arg \max}\hskip0.1em \unicode{x1D53C}\left[\sum \limits_{t={t}_c}^T{R}_t\right]. $$

The subscript $ t $ in $ {\pi}_t^{\ast } $ denotes a non-stationary control policy, as the optimal action depends on the remaining time before the end of the planning horizon, as well as the time $ {t}_c $ during which the transition dynamics is updated, and the control policy is recomputed. The additive reward function to be maximized is defined as follows:

(3.10) $$ R\left(d,u\right)={R}^{\mathrm{health}}(d)+\xi {R}^{\mathrm{control}}(u),\hskip2em \forall d\in \mathcal{D},\forall u\in \mathcal{U}, $$

where $ {R}^{\mathrm{health}}(d) $ and $ {R}^{\mathrm{control}}(u) $ quantify the rewards associated with the health state of the system and the action taken, respectively. In general, these reward values may represent real-world costs or abstract values, tuned using the weighting factor $ \xi \in \mathrm{\mathbb{R}} $ to achieve the desired system performance.

One way to characterize an MDP is to consider the expected utility associated with a policy $ \pi $ when starting in a certain state and following $ \pi $ thereafter. For a finite horizon $ T $ , the state-value function $ {V}^{\pi }:\mathcal{D}\to \mathrm{\mathbb{R}} $ quantifies, for every state and starting time $ t $ , the expected cumulative reward when starting at $ {d}_t $ and following policy $ \pi $ , as follows (Sutton and Barto, Reference Sutton and Barto2018):

(3.11) $$ {V}^{\pi}\left({d}_t\right)={\unicode{x1D53C}}_{\pi}\left[\left.\sum \limits_{\tau =t}^T{R}_{\tau}\right|\hskip0.2em {d}_t\right], $$

where $ {\unicode{x1D53C}}_{\pi}\left[\cdot \right] $ denotes the expected value, given that the agent follows policy $ \pi $ . For any control policy $ \pi $ and for each possible state $ {d}_t $ , with $ t<T $ , the state-value function satisfies the following recursive relationship, known as the Bellman equation for $ {V}^{\pi } $ :

(3.12) $$ {V}^{\pi}\left({d}_t\right)=\sum \limits_{u\in \mathcal{U}}\pi \left(u|{d}_t\right)\sum \limits_{d_{t+1}\in \mathcal{D}}p\left({d}_{t+1}|{d}_t,u\right)\left({r}_t+{V}^{\pi}\left({d}_{t+1}\right)\right), $$

where $ {r}_t $ is the immediate reward at time $ t $ , and $ p\left({d}_{t+1}|{d}_t,{u}_t\right)\in \mathcal{P} $ are the transition probabilities. According to the principle of maximum expected utility, for each time $ t<T $ , there exists a time-dependent optimal policy $ {\pi}_t^{\ast } $ that maximizes expected utility. The time-dependent optimal action at state $ {d}_t $ is given as:

(3.13) $$ {\pi}_t^{\ast}\left({d}_t\right)=\underset{u\in \mathcal{U}}{\arg \max}\left(\sum \limits_{d_{t+1}\in \mathcal{D}}p\left({d}_{t+1}|{d}_t,u\right)\left({r}_t+{V}^{\ast}\left({d}_{t+1}\right)\right)\right), $$

where $ {V}^{\ast}\left({d}_{t+1}\right) $ measures the utility associated with $ {d}_{t+1} $ when following $ {\pi}^{\ast } $ . The Bellman optimality equation for $ {\pi}^{\ast } $ can thus be written as:

(3.14) $$ {V}^{\pi \ast}\left({d}_t\right)=\underset{u\in \mathcal{U}}{\max}\left(\sum \limits_{d_{t+1}\in \mathcal{D}}p\left({d}_{t+1}|{d}_t,u\right)\left({r}_t+{V}^{\ast}\left({d}_{t+1}\right)\right)\right), $$

yielding a system of $ \mid \mathcal{D}\mid $ recursions in $ \mid \mathcal{D}\mid $ unknowns. The optimal state-value functions are the unique solutions to these equations, from which the corresponding optimal policy can be directly derived. However, the optimal action $ {u}_t^{\ast } $ is selected through the $ \max $ operator, which introduces nonlinearity in the system of Bellman equations.

An approximation of the unique solution to Equation (3.14) is computed using finite-horizon dynamic programming. Dynamic programming is an algorithmic paradigm introduced by Bellman (Reference Bellman1957a) that simplifies the optimization by decomposing it into smaller sub-problems that are solved recursively. Among these methods, we employ the value iteration algorithm, also known as backward induction, as it proceeds by solving the problem backward in time, from $ t=T $ down to $ t={t}_c $ .

An algorithmic description of the offline and online phases of the proposed computational procedure is reported in Algorithm 1 and Algorithm 2, respectively.

Remark 3.1. For each action, the transition matrix $ {\mathcal{P}}_u $ is parametrized by a statistic of the inferred transition posteriors leading to different control policies. The parametrizing statistic is directly connected to the risk aversion of the decision-maker, which involves balancing a trade-off between maximizing profits and ensuring structural safety. In the numerical experiments presented in Section 4, the transition matrices are parameterized using the maximum a posteriori estimate of the inferred Dirichlet posterior. Alternative choices could be adopted to induce different levels of risk aversion. For univariate probability distributions it is possible to incorporate different risk measures as proposed in (Tezzele et al., Reference Tezzele, Carr, Topcu and Willcox2026). For the Dirichlet distribution, due to its multivariate nature and its support on simplexes, it is not possible to extend commonly used risk measures such as the value at risk (Morgan, Reference Morgan1996) or the conditional value at risk (Rockafellar and Uryasev, Reference Rockafellar and Uryasev2000) without using an intermediate scalar function.

4.1. Physical asset

The Hörnefors railway bridge, shown in Figure 5, is an integral concrete bridge located along the Swedish Bothnia line. The bridge features a span of $ 15.7\mathrm{m} $ , a free height of $ 4.7\mathrm{m} $ , and a width of $ 5.9\mathrm{m} $ (edge beams excluded). The structural elements have a thickness of $ 0.5\mathrm{m} $ for the deck, $ 0.7\mathrm{m} $ for the frame walls, and $ 0.8\mathrm{m} $ for the wing walls. The bridge is founded on two plates, connected by two stay beams and supported by pile groups. The concrete is of class C35/45, with mechanical properties: $ E=34\mathrm{GPa} $ , $ \nu =0.2 $ , $ \rho =2500\mathrm{kg}/{\mathrm{m}}^3 $ . The superstructure consists of a single track with sleepers spaced $ 0.65\mathrm{m} $ apart, resting on a $ 0.6\mathrm{m} $ deep, $ 4.3\mathrm{m} $ wide ballast layer with a density of $ {\rho}_B=1800\mathrm{kg}/{\mathrm{m}}^3 $ . The bridge is subjected to the transit of Gröna Tåget trains traveling at speeds ranging from $ v\in \left[\mathrm{160,215}\right]\mathrm{km}/\mathrm{h} $ . We specifically consider trains composed of two wagons, totaling $ 8 $ axles, with each axle carrying a mass of $ \psi \in \left[16,22\right]\mathrm{ton} $ . The geometrical and mechanical data, as well as the traveling load model, are adapted from (Rosafalco et al., Reference Rosafalco, Torzoni, Manzoni, Mariani and Corigliano2021; Torzoni et al., Reference Torzoni, Tezzele, Mariani, Manzoni and Willcox2024).

Figure 5.

Physical asset and its digital twin—The physical space corresponds to the Hörnefors railway bridge. The digital space is a schematization for structural health monitoring, including details of synthetic recordings related to displacements $ {u}_1(t),\dots, {u}_{10}(t) $ , and predefined damage regions $ {\Omega}_1,\dots, {\Omega}_6 $ .

Figure 5. Long description

The top panel, labeled Physical space, is a photo of a concrete railway bridge with a passenger train crossing over it. Red dots are superimposed on the bridge’s concrete abutments and the underside of the span to indicate sensor placement.

The bottom panel, labeled Digital space, is a three-dimensional isometric schematic of the same bridge. The model is divided into six purple-shaded damage regions labeled Omega sub 1 through Omega sub 6.

* Omega sub 1 is at the base of the left support.

* Omega sub 2 is at the top corner of the left support.

* Omega sub 3, Omega sub 4, and Omega sub 5 are spaced along the horizontal span.

* Omega sub 6 is at the base of the right support.

Ten red spherical markers represent synthetic recording points for displacements, each with a red arrow indicating the direction of measurement.

* u sub 1 (t) and u sub 5 (t) are on the left support.

* u sub 2 (t), u sub 6 (t), u sub 7 (t), u sub 8 (t), and u sub 9 (t) are distributed along the horizontal span.

* u sub 3 (t), u sub 4 (t), and u sub 10 (t) are on the right support.

A coordinate system at the base of the left support defines the x, y, and z axes.

The integral design of the Hörnefors bridge minimizes the ingress of moisture, debris, and corrosive agents, thereby enhancing its long-term durability. The types of damage that can affect such bridges depend on material and structural factors, environmental conditions, and traffic loads. Some damage patterns that can be described through a localized stiffness reduction include: concrete cracking caused by thermal gradients, freeze–thaw cycles, or overloading; gradual concrete degradation due to the chemical reaction between reactive silica aggregates and alkalis in the cement paste, potentially resulting in cracking and spalling; concrete cracking and structural damage induced by stress concentrations from differential foundational settlements; erosion of concrete surfaces due to vehicular collisions or prolonged exposure to environmental factors.

4.2. Dataset assembly and training options

The monitoring system is deployed to provide displacement recordings $ \mathbf{U} $ , corresponding to the $ {N}_s=10 $ degrees of freedom highlighted in Figure 5. Each recording spans a time interval $ \left[0,T=1.5\mathrm{s}\right] $ , with a sampling frequency of $ 400\mathrm{Hz} $ . To simulate measurement errors caused by sensor self-noise, the recordings are corrupted with an additive Gaussian noise to yield a signal-to-noise ratio of $ 120 $ . An exemplary set of displacement time histories is shown in Figure 6.

Figure 6.

Exemplary multivariate time series of structural vibrations (standardized displacements $ {u}_1(t),\dots, {u}_{10}(t) $ ) induced by the transit of a pilot train across the railway bridge.

Figure 6. Long description

The graph displays ten individual time series plots arranged in two columns. The x-axis for all plots represents Time in seconds, ranging from 0 to 1.5. Each plot features a horizontal zero-baseline and a blue signal line representing standardized displacement.

Left Column (Bottom to Top):

* u sub 1: Shows three distinct downward-dipping wave cycles between 0 and 1.1 seconds, followed by low-amplitude noise.

* u sub 2: Displays three oscillatory cycles with sharp peaks and troughs, stabilizing after 1.1 seconds.

* u sub 3: Similar to u sub 2, showing three high-amplitude oscillations that dampen into noise at the 1.1-second mark.

* u sub 4: Features three broad upward-trending peaks before returning to the baseline.

Right Column (Bottom to Top):

* u sub 5: Shows three irregular downward spikes with jagged recovery phases.

* u sub 6: Displays three sawtooth-like downward oscillations with a gradual return to baseline.

* u sub 7: Shows three smooth, symmetrical downward-curving waves.

* u sub 8: Similar to u sub 7, featuring three smooth downward oscillations.

* u sub 9: Displays three broader, slightly more jagged downward wave cycles.

* u sub 10: Shows three sharp, narrow downward spikes with rapid returns to the baseline.

In all ten plots, the significant vibration activity occurs between 0 and approximately 1.1 seconds, after which the signals return to a steady-state noise level near the zero-baseline.

The parameter $ y\in \left\{0,\dots, 6\right\} $ labels the specific damage region. In addition to the damage-free baseline $ y=0 $ , the presence of damage in the structure is simulated by considering a set of $ {N}_{\Omega}=6 $ predefined damageable subdomains $ {\Omega}_m $ , with $ m=1,\dots, 6 $ , as depicted in Figure 5 and corresponding to $ y=1,\dots, 6 $ . Within these regions, a localized reduction of the material stiffness can occur with a magnitude $ \delta \in \left[30\%,80\%\right] $ , and is kept fixed while a train travels across the bridge.

The full-order model is the same as that described in Torzoni et al., (Reference Torzoni, Tezzele, Mariani, Manzoni and Willcox2024). It features $ \mathrm{17,292} $ degrees of freedom, resulting from a space discretization using linear tetrahedral finite elements with an element size of $ 0.80\mathrm{m} $ , reduced to $ 0.15\mathrm{m} $ for the deck. The ballast layer is modeled through an increased density for the deck and the edge beams. The embankments are modeled using distributed springs over the surfaces in contact with the ground, implemented as a Robin mixed boundary condition (with elastic coefficient $ {10}^8\mathrm{N}/{\mathrm{m}}^3 $ ). The structural dissipation is modeled using a Rayleigh damping matrix, calibrated to achieve a $ 5\% $ damping ratio for the first two structural modes.

The solution time interval is uniformly partitioned into $ 600 $ time steps, and the time integration is performed using an implicit Newmark scheme (Hughes, Reference Hughes2000). We use the (unconditionally stable) constant average acceleration rule, so that the integration step can be arbitrarily selected based on the structural frequencies and the sensor sampling rate.

To reduce the computational cost of solving the full-order model for any sampled $ \boldsymbol{\mu} $ , a projection-based reduced-order model is used to speed up the dataset generation for training the DL models. We use a proper orthogonal decomposition Galerkin reduced basis method (Quarteroni et al., Reference Quarteroni, Manzoni and Negri2015; Chinesta et al., Reference Chinesta, Huerta, Rozza, Willcox, Stein, de and Hughes2017). The reduced-order model is derived from a snapshot matrix assembled from $ 400 $ full-order solutions for varying input parameters $ \boldsymbol{\mu} ={\left(v,\psi, y,\delta \right)}^{\top } $ , taken as uniformly distributed and sampled via the Latin hypercube rule. The order of the reduced-order expansion is determined using an energy criterion. By prescribing a tolerance of $ {10}^{-3} $ for the fraction of energy content to be disregarded, the number of degrees of freedom is reduced to $ 133 $ . Using this reduced-order model, the training dataset is populated with $ \mathrm{10,000} $ instances. These data are subsequently used to train DL models for detecting, locating, and quantifying structural damage, as described in (Torzoni et al., Reference Torzoni, Tezzele, Mariani, Manzoni and Willcox2024).

The application of DL in SHM follows a pattern recognition paradigm (Bishop, Reference Bishop2006), where damage is assessed by comparing incoming measurements with data previously collected under known structural conditions. By capturing temporal correlations within and across time series data, DL optimizes the extraction of damage-sensitive features for modeling the relationship between inputs and damage states. Moreover, as DL models are trained offline, they provide real-time predictions during operation.

Within our framework, each time new observational data are acquired—when a train crosses the bridge, marking a new DT time slice—they are first processed by a classification model to address damage detection and localization. The target class corresponds to one of the predefined damage scenarios described by the $ y $ parameter. The classifier maps the vibration recordings $ \mathbf{U} $ onto confidence levels representing the likelihood of $ \mathbf{U} $ belonging to each damage class. The most likely class is selected as the best point estimate for categorizing the measurements. Such a classifier is trained by minimizing the categorical cross-entropy loss function between the one-hot encoding of $ y $ and its predicted counterpart.

Whenever damage is detected and localized in the $ m $ -th region $ {\Omega}_m $ , where $ m=1,\dots, 6 $ , the observational data are further processed with a regression model—one for each damageable region—to quantify the extent of damage $ \delta $ . Each regression model maps the vibration recordings associated with the $ m $ -th damage region onto the estimated magnitude of the stiffness reduction within that region. These regressors are trained by minimizing a mean-squared error loss function between the target stiffness reduction and its predicted value. For a detailed and formal description of the DL models, architectures, and their implementation, the reader is referred to (Torzoni et al., Reference Torzoni, Tezzele, Mariani, Manzoni and Willcox2024).

4.3. Digital twin framework

Since the space of digital states is discrete, the outcomes of the DL models are accommodated into the PGM by discretizing the value range of $ \delta $ into $ {N}_{\delta}\in \mathrm{\mathbb{N}} $ intervals. The number and width of these intervals are chosen arbitrarily, without specific restrictions. The $ \delta $ range is discretized into $ {N}_{\delta }=6 $ intervals $ \Big\{\left[30\%,35\%\right], $ $ \left[35\%,45\%\right], $ $ \left[45\%,55\%\right], $ $ \left[55\%,65\%\right], $ $ \left[65\%,75\%\right], $ $ \left[75\%,80\%\right]\Big\} $ . This results in $ \mid \mathcal{D}\mid =1+{N}_{\delta }{N}_{\Omega}=37 $ digital states, ordered first by damage location and then by damage level.

In the absence of experimental data, the testing phase of the DL models is carried out using $ 4000 $ noisy full-order solutions. The ground truth digital state is classified with an accuracy of $ 91.39\% $ .

We consider four different control actions, each associated with a CPT that models the transition probability $ p\left({D}_{t+1}|{D}_t,{U}_t={u}_t\right) $ . Together, these CPTs encode the $ {\phi}_t^{\mathrm{history}} $ factor. Note that these internal models predict the evolution of structural health but do not replicate the ground truth progression, which remains unknown to the DT. The considered control actions are as follows:

1. 1. Do nothing (DN) action: no maintenance is scheduled, and the structural health state evolves according to a stochastic deterioration process, while regular revenue is maintained.

2. 2. Restrict operations (RO) action: traffic flow on the bridge is restricted to lightweight trains carrying less than $ 18\mathrm{ton} $ per axle. This reduces the structure degradation rate but also decreases the revenue generated by the infrastructure.

3. 3. Minor repair (MR) action: a moderate-cost maintenance intervention is carried out to address specific damages, though the structure does not return to a damage-free state.

4. 4. Perfect repair (PR) action: a high-cost maintenance intervention is performed, almost certainly restoring the bridge to a damage-free state.

The (unknown) ground truth stochastic model varies depending on the most recent control actions. We assume that damage can develop in any of the predefined regions, without propagating between different damageable subdomains. Based on this assumption, we define below the prescribed degradation (or improvement) evolution model associated with each action.

4.3.1. Do nothing action

In the absence of maintenance actions, a (simulated) stochastic deterioration process monotonically degrades the structural health. If the (unknown) state of the system is undamaged, the presence/location of damage $ y $ is sampled from a categorical distribution. This distribution assigns half of its probability mass to the undamaged configuration $ y=0 $ , indicating that the bridge remains intact. The remaining half is evenly distributed among the remaining damage classes $ y\ne 0 $ . Damage may develop with magnitude $ \delta $ uniformly sampled within the range of the first damage interval $ \left[30\%,35\%\right] $ . In contrast, if the system is already damaged, further deterioration propagates through $ \delta $ increments sampled from a categorical distribution with Dirichlet hyperparameters $ {\boldsymbol{\alpha}}_{\mathrm{DN}}^{\mathrm{true}}={\left({\alpha}_{\mathrm{DN}}^{0,\mathrm{true}},{\alpha}_{\mathrm{DN}}^{1,\mathrm{true}},{\alpha}_{\mathrm{DN}}^{2,\mathrm{true}}\right)}^{\top } $ , scaled by the $ \delta $ discretization step $ 10\% $ , and finally perturbed by Gaussian noise with mean $ 2\% $ and standard deviation $ 0.5\% $ (negative increments are rounded to zero). The DN model can be formalized as follows:

(4.1) $$ y\sim \mathrm{Categorical}\left(\frac{1}{2},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12}\right), $$

(4.2) $$ {\delta}_0\mid y\ne 0\sim \mathrm{Uniform}\left(\mathrm{0.3,0.35}\right), $$

(4.3) $$ {\delta}_t-{\delta}_{t-1}\mid y\ne 0\sim 0.1\cdot \mathrm{Categorical}\left({\mathbf{p}}_{\mathrm{DN}}^{\mathrm{true}}\right)+{\mathrm{Normal}}_{\ge 0}\left(\mathrm{0.02,0.005}\right), $$

(4.4) $$ {\mathbf{p}}_{\mathrm{DN}}^{\mathrm{true}}\sim \mathrm{Dirichlet}\left({\boldsymbol{\alpha}}_{\mathrm{DN}}^{\mathrm{true}}\right),\hskip1em {\boldsymbol{\alpha}}_{\mathrm{DN}}^{\mathrm{true}}={\left(\mathrm{500,400,200}\right)}^{\top }. $$

The prescribed trajectory of the structural health parameters is arbitrarily chosen to fully showcase the capabilities of the DT. However, the DT would be equally capable of tracking the structural health evolution, whether considering more or less aggressive degradation processes.

4.3.2. Restrict operation action

If only lightweight trains are allowed to cross the bridge, the model resembles the one for the DN action but is characterized by a lower probability of damage initiation and a reduced deterioration rate. In this scenario, damage may occur with a probability of $ 0.25 $ , evenly distributed among the damage classes. Once damage develops, it propagates with $ \delta $ increments sampled from a categorical distribution with Dirichlet hyperparameters $ {\boldsymbol{\alpha}}_{\mathrm{RO}}^{\mathrm{true}}={\left({\alpha}_{\mathrm{RO}}^{0,\mathrm{true}},{\alpha}_{\mathrm{RO}}^{1,\mathrm{true}},{\alpha}_{\mathrm{RO}}^{2,\mathrm{true}}\right)}^{\top } $ , scaled by the $ \delta $ discretization step $ 10\% $ , and finally perturbed by Gaussian noise with mean $ 2\% $ and standard deviation $ 0.5\% $ (negative increments are rounded to zero). The RO model is formalized as follows:

(4.5) $$ y\sim \mathrm{Categorical}\left(\frac{3}{4},\frac{1}{24},\frac{1}{24},\frac{1}{24},\frac{1}{24},\frac{1}{24},\frac{1}{24}\right), $$

(4.6) $$ {\delta}_0\mid y\ne 0\sim \mathrm{Uniform}\left(\mathrm{0.3,0.35}\right), $$

(4.7) $$ {\delta}_t-{\delta}_{t-1}\mid y\ne 0\sim 0.1\cdot \mathrm{Categorical}\left({\mathbf{p}}_{\mathrm{RO}}^{\mathrm{true}}\right)+{\mathrm{Normal}}_{\ge 0}\left(\mathrm{0.02,0.005}\right), $$

(4.8) $$ {\mathbf{p}}_{\mathrm{RO}}^{\mathrm{true}}\sim \mathrm{Dirichlet}\left({\boldsymbol{\alpha}}_{\mathrm{RO}}^{\mathrm{true}}\right),\hskip1em {\boldsymbol{\alpha}}_{\mathrm{RO}}^{\mathrm{true}}={\left(\mathrm{600,300,100}\right)}^{\top }. $$

4.3.3. Minor repair action

If the structural health state of the bridge undergoes partial improvement and $ y=0 $ , the system remains undamaged. If $ y\ne 0 $ , the damage magnitude $ \delta $ decreases. Specifically, the $ \delta $ decrement is sampled from a categorical distribution with Dirichlet hyperparameters $ {\boldsymbol{\alpha}}_{\mathrm{MR}}^{\mathrm{true}}={\left({\alpha}_{\mathrm{MR}}^{0,\mathrm{true}},{\alpha}_{\mathrm{MR}}^{0,\mathrm{true}},{\alpha}_{\mathrm{MR}}^{0,\mathrm{true}}\right)}^{\top } $ , scaled by the $ \delta $ discretization step $ 10\% $ , and finally perturbed by Gaussian noise with mean $ -2\% $ and standard deviation $ 0.5\% $ (positive increments are rounded to zero). The state is updated to reflect a partial repair, ensuring that the damage magnitude does not fall below the undamaged state threshold of $ \delta =30\% $ . The MR model is formalized as follows:

(4.9) $$ {\delta}_t-{\delta}_{t-1}\mid y\ne 0\sim -0.1\cdot \mathrm{Categorical}\left({\mathbf{p}}_{\mathrm{MR}}^{\mathrm{true}}\right)+{\mathrm{Normal}}_{\le 0}\left(-\mathrm{0.02,0.005}\right), $$

(4.10) $$ {\mathbf{p}}_{\mathrm{MR}}^{\mathrm{true}}\sim \mathrm{Dirichlet}\left({\boldsymbol{\alpha}}_{\mathrm{MR}}^{\mathrm{true}}\right),\hskip1em {\boldsymbol{\alpha}}_{\mathrm{MR}}^{\mathrm{true}}={\left(\mathrm{1,100,300}\right)}^{\top }. $$

4.3.5. Transition model

All numerical results are obtained using a Dirichlet–Multinomial state-independent transition model. Transition probabilities thus depend only on the selected action and not on the current state, and are learned online through action-specific Dirichlet priors updated via observed state transitions.

The assumed transition model $ p\left({D}_{t+1}|{D}_t,{U}_t={u}_t\right) $ satisfies the assumption that only one damage scenario can be activated at a time. The transition dynamics associated with the different actions are formulated as described in Section 3.1. Starting from an initial belief over the transition probabilities, these are subsequently refined by integrating evidence about the system response to applied actions. Specifically, we instantiate two Dirichlet–Multinomial state-independent models to represent degradations of zero, one, or two $ \delta $ intervals under the DN and RO actions, respectively, and one Dirichlet–Multinomial state-independent model to represent improvements of zero, one, or two $ \delta $ intervals under the MR action. The initial Dirichlet hyperpriors are assigned as:

(4.11) $$ {\boldsymbol{\alpha}}_{\mathrm{DN}}={\left({\alpha}_{\mathrm{DN}}^0,{\alpha}_{\mathrm{DN}}^1,{\alpha}_{\mathrm{DN}}^2\right)}^{\top }={\left(\mathrm{10,7,4}\right)}^{\top }, $$

(4.12) $$ {\boldsymbol{\alpha}}_{\mathrm{RO}}={\left({\alpha}_{\mathrm{RO}}^0,{\alpha}_{\mathrm{RO}}^1,{\alpha}_{\mathrm{RO}}^2\right)}^{\top }={\left(\mathrm{10,5,2}\right)}^{\top }, $$

(4.13) $$ {\boldsymbol{\alpha}}_{\mathrm{MR}}={\left({\alpha}_{\mathrm{MR}}^0,{\alpha}_{\mathrm{MR}}^1,{\alpha}_{\mathrm{MR}}^2\right)}^{\top }={\left(\mathrm{3,10,9}\right)}^{\top }, $$

where $ {\boldsymbol{\alpha}}_{\mathrm{DN}} $ and $ {\boldsymbol{\alpha}}_{\mathrm{RO}} $ prioritize remaining in the same state over one-step and two-step degrading transitions, while $ {\boldsymbol{\alpha}}_{\mathrm{MR}} $ prioritizes one-step and two-step recovery over remaining in the same state. Finally, for the PM action, we prescribe a transition matrix that assigns unit probability to transitioning to the undamaged state, independently of the current condition.

Transitions are defined in terms of discrete step variations, corresponding to changes of zero, one, or two $ \delta $ intervals. The state space is bounded by construction, so that any transition exceeding the admissible range is mapped onto the nearest valid state, that is, the undamaged or most-damaged condition. Moreover, it is worth highlighting that our state-independent assumption is appropriate when the deterioration and repair dynamics are primarily governed by the applied action rather than the current damage level, or when available data are insufficient to reliably estimate fully state-dependent transitions. Under these conditions, pooling observations across states improves statistical efficiency and yields more stable learning, while still capturing the dominant degradation and recovery trends.

4.3.6. Reward model

At each time step, the DT selects a control action $ {u}_t\in \left\{\mathrm{DN},\mathrm{RO},\mathrm{MR},\mathrm{PR}\right\} $ . Choosing a DN or RO action yields a positive reward but also entails the risk of structural deterioration. Conversely, the MR and PR actions mitigate structural health deterioration but incur a negative reward. The costs associated with the structural health state, encoded in $ {D}_t $ , and the costs related to the control actions $ {U}_t $ , appearing in Equation (3.10), are defined as follows:

(4.14) $$ {R}_t^{\mathrm{health}}(d)=\left\{\begin{array}{ll}0& \mathrm{if}\;y=0,\delta =0,\\ {}-\exp \left(5\delta \right)+4& \mathrm{if}\;y\ne 0,0<\delta <80\%,\\ {}-250& \mathrm{if}\;\delta \ge 80\%,\end{array}\right.\hskip2em {R}_t^{\mathrm{control}}(u)=\left\{\begin{array}{ll}+30& \mathrm{if}\;u=\mathrm{DN},\\ {}+18& \mathrm{if}\;u=\mathrm{RO},\\ {}-75& \mathrm{if}\;u=\mathrm{MR},\\ {}-250& \mathrm{if}\;u=\mathrm{PM}.\end{array}\right. $$

These non-dimensional rewards represent indicative costs charged to the decision-maker. Actual values can be derived from service and cost lists provided by state agencies and companies (KG Papakonstantinou and Shinozuka, Reference Papakonstantinou and Shinozuka2014). For instance, $ {R}^{\mathrm{health}} $ is designed to penalize the progressive deterioration of the structural condition as a function of $ \delta $ . These costs may reflect various impacts, such as a reduction in service level, workplace accidents, or increased structural failure probability (KG Papakonstantinou and Shinozuka, Reference Papakonstantinou and Shinozuka2014), with the last entry imposing a significantly negative reward for excessively compromised structural states.

4.4. Results

In this section, we present the results of DT simulations carried out over an expected lifespan of $ T=60 $ time steps. At each time step, new observational data are simulated based on the (unknown) ground truth structural health and the most recently applied control action. The DT assimilates these data, estimates the digital state, and recommends the next control action.

By running a cluster of $ 100 $ simulations, each spanning $ 60 $ time steps and initialized with a different random seed for the ground-truth physical process, both the DT using a finite-horizon dynamic programming policy without precision updates—computed based on the initial prior transition model—and the adaptive DT using a finite-horizon MBRL policy with precision updates, where control policies are iteratively recomputed using the updated transition model, exhibit zero failures. Figure 7 compares the evolution of cumulative rewards, showing that the adaptive DT yields significantly higher cumulative rewards (+17%). This improvement stems from the online refinement of the transition model, which enables more informed planning and decision-making. Moreover, our experiments indicate that non-adaptive policies can lead to undesirable outcomes; depending on the interplay between the transition model, the reward model, and the ground-truth dynamics, such policies may cause the DT to fail, that is, the structural state reaches a critically compromised condition. For instance, digital failure (corresponding to $ \delta \ge 80\% $ and symbolizing structural collapse) can occur if the initial prior over the transition dynamics is poorly aligned with the ground-truth process. In contrast, the robustness of MBRL ensures that the DT consistently maintains the structural health above critical thresholds.

Figure 7.

Digital twin online phase—Comparison of cumulative rewards under control policies computed using finite-horizon model-based reinforcement learning (MBRL) with precision updates and finite-horizon dynamic programming (DP) without precision updates. Solid lines denote the mean cumulative rewards, while the shaded areas indicate the one-standard-deviation credibility intervals. Results are shown for two clusters of simulations, each comprising $ 100 $ runs initialized with different random seeds.

Figure 7. Long description

The X-axis is labeled Time step and ranges from 0 to 60. The Y-axis is labeled Cumulative reward and ranges from negative 1000 to 250.

Two main data series are plotted:

* M B R L with update mean: A blue solid line starting near 0, peaking slightly at time step 5, and then following a steady linear decrease to approximately negative 650 at time step 60. It is surrounded by a light blue shaded area representing the M B R L ensemble s t d.

* D P no update mean: A red solid line that follows a similar trajectory but remains consistently lower than the blue line after time step 5, ending at approximately negative 750 at time step 60. It is surrounded by a light orange shaded area representing the D P ensemble s t d.

At the far right of the graph, an upward-pointing arrow between the red and blue lines is annotated with plus 17 percent. A legend at the bottom left identifies the four visual elements: blue line for M B R L with update mean, red line for D P no update mean, light blue block for M B R L ensemble s t d, and light orange block for D P ensemble s t d.

In-simulation performance is further assessed through a tracking accuracy metric based on the maximum a posteriori estimate of the digital state relative to the ground-truth physical process. To account for the $ \delta $ discretization, accuracy is computed using a $ 10\% $ tolerance, corresponding to the discretization step. Under these conditions, the adaptive DT achieves a mean accuracy of $ 93.6\% $ with a standard deviation of $ 4\% $ , while the non-adaptive counterpart attains a mean accuracy of $ 93.15\% $ with a standard deviation of $ 4.3\% $ . To avoid biasing the metric, true undamaged cases are excluded from this computation, as they are easier to identify and would favor conservative policies that promote occupancy of the undamaged state. The similar tracking accuracies observed for the two configurations are expected, since the observation model itself is not modified by the adaptive updates.

Figure 8 illustrates a sample DT simulation using finite-horizon MBRL with precision updates. The DT accurately tracks the digital state evolution with relatively low uncertainty and a limited delay of at most two time steps with respect to the ground truth. This delay is attributed both to the need to update the prior belief based on data from previous time steps and to occasional misclassification by the DL model between digital states associated with the same damage location but adjacent $ \delta $ intervals. The figure also demonstrates how accurate inference of transition dynamics and optimization of finite-horizon control policies allow the PR action to be avoided as the simulation approaches the end of the planning horizon (i.e., the decommissioning time). Specifically, the belief of selecting a PR action for the same damage magnitude progressively decreases as the simulation advances. Additionally, Figure 9 shows the predicted evolution of the digital state and control actions starting from $ {t}_c=60 $ and spanning a prediction horizon of $ {t}_p={t}_c+5 $ . The DT first forecasts the expected improvement in structural health after taking the MR action, based on the learned transition model. This is followed by a high-probability prediction of DN and MR actions, yielding a pseudo-stable belief propagation pattern that maintains most of the probability mass of the digital state within the first three $ \delta $ intervals.

Figure 8.

Digital twin online phase—Simulation under finite-horizon model-based reinforcement learning with precision updates. Probabilistic and best point estimates of: (top four panels) digital state evolution against the ground truth physical state for regions $ {\Omega}_1 $ , $ {\Omega}_2 $ , $ {\Omega}_3 $ , and $ {\Omega}_6 $ ; (bottom panel) control actions informed by the digital twin against the optimal action under the ground truth. The background colors represent the digital state beliefs in the top panels, and the control action beliefs in the bottom panel.

Figure 8. Long description

A vertical stack of five panels with a shared X-axis representing Time step from 0 to 60. A vertical color bar on the right indicates probability from 0.0 in light blue to 1.0 in dark blue.

* The top four panels show Damage probability on the Y-axis from 0% to 80% for regions Omega sub 1, Omega sub 2, Omega sub 3, and Omega sub 6. Each features a dashed black line for Ground truth and a solid red line for Digital twin. Background blue shading represents belief density.

* Region Omega sub 1: Damage remains at 0% until time step 20, peaks at 60% around step 25, then returns to 0%.

* Region Omega sub 2: Similar to Omega sub 1, but the peak occurs slightly later, between steps 25 and 30.

* Region Omega sub 3: Damage is 0% until step 28, then fluctuates rapidly between 20% and 60% for the remainder of the timeline. A label Multiple steps points to these oscillations.

* Region Omega sub 6: Shows two distinct peaks of approximately 60% damage between steps 5 and 15, remaining at 0% otherwise.

* The bottom panel, Action probability, has a Y-axis with categorical labels: D N, P R, R O, and M R. The digital twin red line closely tracks the ground truth dashed line, showing frequent switching between actions, particularly after time step 30, corresponding to the damage fluctuations in Omega sub 3.

Figure 9.

Digital twin online phase—Future predictions starting at $ {t}_c=60 $ under finite-horizon model-based reinforcement learning with precision updates. The background colors represent the digital state beliefs in the top panels, and the control action beliefs in the bottom panel.

Figure 9. Long description

A multi-panel visualization with a shared color scale on the right ranging from 0.0 (light blue) to 1.0 (dark blue).

Top Panel: Damage probability in region Omega sub 3.

* The Y-axis represents damage levels from 0% to 80% with a separate bottom row for No Damage.

* At time step 60, there is a high probability (dark blue) concentrated at the 60% damage level.

* From time steps 61 to 65, the probability shifts downward and becomes more distributed, with a moderate concentration (medium blue) appearing around the 40% damage level.

Bottom Panel: Action probability.

* The Y-axis lists four actions: M R, R O, P R, and D N.

* At time step 60, the action M R has the highest probability (dark blue).

* From time steps 61 to 65, the probability for M R decreases significantly (light blue), while the probability for action D N increases and remains dominant (medium blue) through the end of the sequence.

* Actions R O and P R maintain near-zero probability (lightest blue) throughout the entire period.

Figure 10 illustrates the initial Dirichlet priors on the transition probabilities for the DN, RO, and MR actions, together with the corresponding posteriors at the end of the simulation. In all cases, the posteriors exhibit clear convergence toward the (unknown) ground-truth transition probabilities while significantly reducing uncertainty. The corresponding posterior parameters read:

(4.15) $$ {\overline{\boldsymbol{\alpha}}}_{\mathrm{DN}}={\left({\alpha}_{\mathrm{DN}}^0+{N}_{\mathrm{DN}}^0,{\alpha}_{\mathrm{DN}}^1+{N}_{\mathrm{DN}}^1,{\alpha}_{\mathrm{DN}}^2+{N}_{\mathrm{DN}}^2\right)}^{\top }, $$

(4.16) $$ {\overline{\boldsymbol{\alpha}}}_{\mathrm{RO}}={\left({\alpha}_{\mathrm{RO}}^0+{N}_{\mathrm{RO}}^0,{\alpha}_{\mathrm{RO}}^1+{N}_{\mathrm{RO}}^1,{\alpha}_{\mathrm{RO}}^2+{N}_{\mathrm{RO}}^2\right)}^{\top }, $$

(4.17) $$ {\overline{\boldsymbol{\alpha}}}_{\mathrm{MR}}={\left({\alpha}_{\mathrm{MR}}^0+{N}_{\mathrm{MR}}^0,{\alpha}_{\mathrm{MR}}^1+{N}_{\mathrm{MR}}^1,{\alpha}_{\mathrm{MR}}^2+{N}_{\mathrm{MR}}^2\right)}^{\top }, $$

whose probability masses are centered around the true modes, indicating that the updated beliefs tend to reflect the actual transition probabilities while reducing uncertainty.

Figure 10.

Prior and posterior distributions of the transition probabilities associated with the do nothing action (first row), the restrict operations action (second row), and the minor repair action (third row). Initial priors correspond to $ t=0 $ , while the posterior distributions correspond to $ t=60 $ .

Figure 10. Long description

The figure consists of three rows of paired ternary plots. Each plot has vertices labeled 0 steps (bottom-left), 1 step (bottom-right), and 2 steps (top). A color scale on the right of each row indicates density from 0 (white) to Max (dark blue).

* Row 1: Transition probability evolution for the D N action. The Prior plot shows a broad, medium-blue density concentrated near the bottom-center, slightly toward 0 steps. The Posterior plot shows a much tighter, darker blue concentrated peak shifted further toward the bottom-center, indicating increased certainty.

* Row 2: Transition probability evolution for the R O action. The Prior plot shows a wide density area in the bottom-left quadrant. The Posterior plot shows the density has narrowed and intensified in the same bottom-left region, closer to the 0 steps vertex.

* Row 3: Transition probability evolution for the M R action. The Prior plot shows an elongated density peak in the bottom-right quadrant, leaning toward the 1 step vertex. The Posterior plot shows this density has become more concentrated and moved higher along the right edge toward the 2 steps vertex.

A dedicated sensitivity analysis to observation uncertainty is not repeated in the present work; instead, we refer the reader to (Torzoni et al., Reference Torzoni, Maisto, Manzoni, Donnarumma, Pezzulo and Corigliano2026), where this aspect has been systematically investigated. The analysis considers progressively degraded observation conditions, including incomplete data streams during train passages, varying observation failure rates, and uncertainty in damage extent. The results show a gradual performance degradation: reduced data availability increases state uncertainty and may delay optimal actions, while the transition model helps preserve acceptable performance even under severe data loss, although with increased uncertainty. On the other hand, uncertainty in damage extent mainly affects estimation accuracy, with negligible impact on failure rates.