Inference and dynamic decision-making for deteriorating systems
with probabilistic dependencies through Bayesian networks
and deep reinforcement learning

The inspection and maintenance (I&M) planning decision-making problem is formulated here as a Partially Observable Markov Decision Process (POMDP), whose transition and observation models are defined by Bayesian network structures. POMDPs provide a principled mathematical framework for optimal planning and decision-making under uncertainty, formally specified by the tuple $⟨S,A,O,T,Z,R,\gamma ⟩$. A decision maker (henceforth agent) interacts with a stochastic environment, described by the state $s\in S$, taking actions $a\in A$ over a finite or infinite horizon $t_N$. The dynamics correspond to those in a MDP: At each time step $t$, an agent takes an action $a_t\in A$, and the environment evolves from state $s_t\in S$ to state $s_{t+1}\in S$, according to the transition model $T\text{coloneqq}p(s_{t+1}|s_t,a_t)$. In a MDP, the agent receives a reward based on the cost model $R\text{coloneqq}{r}_t(s_t,a_t,s_{t+1})$ discounted by factor $\gamma$, and the objective is to find the policy ${\pi }^{\ast }$ that induces the optimal value function $V^{\ast }\left(s_t\right)$:

In a POMDP, however, states $s\in S$ are not directly observed and, instead, observations $o\in O$ can be collected according to the observation model $Z\text{coloneqq}p(o_{t+1}|s_{t+1},a_t)$. Note that the observation model is the likelihood of collecting an observation $o_{t+1}\in O$ after taking an action $a_t$ and having transitioned to state $s_{t+1}$. In an I&M context, the observation model is often modeled by Probability of Detection (PoD) curves or according to inspection/monitoring measurement noise Morato et al. (2022). A POMDP policy is a mapping of the dynamically updated history of actions and observations, $a_{0:t-1},o_{0:t}$, to the current action $a_t$. This history is sufficiently encoded in belief $b$, which is the probability over system states, b(s). The optimal policy ${\pi }^{\ast }$, therefore, corresponds to the value function Kurniawati et al. (2008) that satisfies the Bellman equation:

where $p(o_{t+1}|b_t,a_t)$ is the probability of collecting an observation $o_{t+1}\in O$ given the belief $b_t$ and action $a_t\in A$. Assuming a Markovian environment is reasonable in most practical applications with the aid of state augmentation techniques (Papakonstantinou and Shinozuka, 2014a), hence, any general I&M planning decision problem can be efficiently formulated as a POMDP. The determination of the optimal I&M policy ${\pi }^{\ast }$ becomes the main objective, inducing a minimization of the expected life-cycle costs $r_{tot}$, by balancing structural failure risk against inspection and maintenance costs:

where $r_{ins}$, $r_{rep}$ and $r_F$ stand for inspection, repair and failure costs, respectively, defined as negative rewards. In terms of utilities, the failure risk $r_F$ is typically defined in a structural reliability context as the annual probability of a failure event weighted by the consequence of a structural failure, which might also include environmental and societal consequences, specified in equivalent units. The definition of the failure risk at the system level will be further elaborated in Section 3.

Existing I&M planning applications often model the deterioration evolution $d$, at the component level, conditional on a set of random variables ${\text{uptheta}}_d$ Straub (2009); Luque and Straub (2019, 2016) or as a function of the deterioration rate $\tau$ Papakonstantinou and Shinozuka (2014b); Andriotis and Papakonstantinou (2019). Both formulations are equivalent for modeling deterioration processes, as already discussed and demonstrated in Morato et al. (2022) and shown in Fig. 1. When observations are collected, through inspections or monitoring, Bayesian updating can be then conducted. Available algorithms allow exact Bayesian inference if the problem is formulated in a discrete state space Murphy (2002), as the computation of Bayes’ normalization constant is a challenging task in continuous state spaces Papakonstantinou et al. (2022). In order to utilize discrete state based algorithms, the involved continuous random variables can be discretized. The quality of the discretization has a huge impact and shall be treated carefully Straub (2009); Morato et al. (2022), especially when the problem deals with rare events, e.g., failure events. In general, an efficient discretization aims at minimizing the computational expense while preserving the required level of accuracy.

In a POMDP, the states cannot be directly observed and the decision maker reasons under partial observability, only informed by a belief $b$, which is defined as the probability over states. At each time step, the belief is dynamically updated, based on Bayes’ rule, depending on the initial belief, $b_t$, the action taken, $a_t$, and the collected observation, $o_t$, following three main steps: (i) the belief evolves according to the transition model $p(s_{t+1}|s_t,a_t)$, (ii) the belief is updated based on the collected observation with probability $p(o_{t+1}|s_{t+1},a_t)$, and (iii) the belief state is normalized. This belief update operation is denoted as forward pass within the context of hidden Markov models Murphy (2002). At the system level, the belief of each component can be updated by implementing the steps listed in Algorithm 1.

State-of-the-art POMDP solvers often require the modeling of the POMDPs in a flat structure, which can be usually encoded by augmenting the state space Papakonstantinou and Shinozuka (2014a), particularly if the process is described by multiple random variables. However, POMDPs can also be formulated in a factored fashion, exploiting the dependence structure among random variables and thus significantly alleviating the required computational effort. We specify here the transition and observation models based on conditional structures described by dynamic Bayesian networks (DBNs), and while the belief state $b$ remains the same as that for flat POMDPs, the transition and observation models are now constructed by taking advantage of the involved dependencies. For instance, the deterioration rate model can be constructed as $p(d_{t+1}|d_t,{\tau }_{t+1})p\left({\tau }_{t+1}|{\tau }_t\right)$ instead of $p(d_{t+1},{\tau }_{t+1}|d_t,{\tau }_t)$. This incorporation of conditional structures allows a reduction of the transition model dimensionality from $|S_d{|}^2|S_{\tau }{|}^2$ to $|S_d{|}^2|S_{\tau }|+|S_{\tau }{|}^2$ and can achieve significant computational benefits when multiple random variables are involved. This formulation can be seamlessly applied to simulate the deterioration environment, as will be explained in Section 4, due to the flexibility naturally offered by the proposed deep reinforcement learning approach.

Existing methods model the deterioration correlation among components either via random fields or through common influencing factors. Whereas the former are particularly useful for applications in which the dependence is attributed to the geometrical distance between components, the latter are more suitable for systems in which identical attributes of physical phenomena, e.g., similar manufacturing techniques or similar loading, lead to shared sources of model uncertainties among the components Luque and Straub (2016). In a hierarchical structure, the deterioration of each component is defined conditional on a set of common influencing variables, shared among all the components and represented at the highest level of the hierarchy. In theory, the state space of a system under deterioration dependence can be modeled directly as the joint space of all the parameters involved in the deterioration process of the system. In this case, the discretized state space would grow exponentially with the number of components $N_C$ included, into a $|S{|}^{N_C}$ dimensional space. To overcome this increase in dimensionality, we adopt the hierarchical Gaussian structure previously proposed in Luque and Straub (2016), in which the belief state of each component is encoded conditional on a hyperparameter (common influencing variable) $\alpha$ or set of hyperparameters $\alpha$. The central idea behind this hierarchical structure is that component beliefs for a given hyperparameter $b\left(s|\alpha \right)$ are independent, enabling an efficient decoupling of the components joint space. This decoupling alleviates the computational complexity from the original joint space $|S{|}^{N_C}$ to a space $|S|\cdot |S_{\alpha }|\cdot N_C$ that grows linearly with the number of considered components $N_C$. Note that the state space includes now the states of the hyperparameter(s), which should also be properly discretized. The increase of the state space due to the incorporation of the hyperparameter(s) is however less significant than when considering the joint state space.

A graphical representation of the proposed hierarchical structure is illustrated in Fig. 1, applicable to deterioration processes modeled either as a function of the deterioration rate or conditional on a set of parameters Morato et al. (2022). In either case, the deterioration process $d$ is encoded conditional on the hyperparameter(s) $\text{upalpha}$, along with the deterioration rate $\tau$ or parameters ${\text{uptheta}}_d$. Evidence collected through observations $o_{d_t}$ does not only serve for updating the damage state, but also for updating the hyperparameters. Since the hyperparameters are parent nodes for all the components, once a component is inspected, the hyperparameters are also updated, influencing all the other components, even those for which evidence was not directly available. The reliability of the system is represented in Fig. 1 by the binary node $F_{sys}$, conditional on the failure state of the components $F^{\left(l\right)}$. At a component level, the failure state is modeled by the binary variable $F^{\left(l\right)}$ and corresponds to the subset of the deterioration space classified as failure $S_F\subseteq S$.

This Gaussian hierarchical structure is a mathematically motivated model induced by the convenient formulation available for normal random variables, i.e., the conditional and joint distributions of normal random variables are also normally distributed. Let us first consider the special case in which the marginal probability of each considered component deterioration is defined as a standard normal random variable $Y_i$. Under correlation, the parameters $Y_i$ are, however, defined as normal random variables with mean ${\lambda }_i\alpha$ and standard deviation $\sqrt{1-{\lambda }_i^2}$Dunnett and Sobel (1955):

Since both $X_i$ and $\alpha$ are independent standard normal random variables, the covariance of $Y_i$ and $Y_j$ can be formulated as:

After removing all the terms associated with zero covariance, i.e., $cov(X_i,X_j)$, $cov(X_i,\alpha )$, and $cov(X_j,\alpha )$, we can define the correlation coefficient between $Y_i$ and $Y_j$ as:

If all the components are equi-correlated, then ${\lambda }_i={\lambda }_j=\sqrt{\rho (Y_i,Y_j)}$, for all $i,j$. This presented Gaussian structure is further generalized in this work for the general case of unequally correlated components, by preserving the validity of Eq. 6. For complex correlation configurations, one hyperparameter $\alpha$ might not be sufficient to satisfy Eq. 6, and in that case, one can incorporate additional hyperparameters $\alpha$, at the expense of a higher computational cost. When multiple hyperparameters $\alpha$ are included, the best fit for $\rho (Y_i,Y_j)={\sum }_{k=1}^m\left({\lambda }_{ik}{\lambda }_{jk}\right)$ can be found via optimization procedures, e.g., least squares Song and Kang (2009). Once the Gaussian correlation structure is specified through the parameters $\lambda$, the cumulative distribution of $Y_i$ conditional on the hyperparameter(s) $\alpha$ can be defined as:

For the cases in which the deterioration process is modeled by random variables other than Gaussian and considering that a Nataf transformation is applicable Luque and Straub (2016), then Eq. 7 can be also redefined as:

where $F_{D_i|\alpha }\left(d_i\right)$ stands for the cumulative distribution function of a variable $d_i$ conditional on the hyperparameter(s) $\text{upalpha}$, and $\Phi$ is the standard normal cumulative distribution function. In a discrete state space, the belief conditional on the hyperparameters is equal to the difference between the cumulative distribution function at the upper boundary and at the lower boundary of each belief interval:

We reformulate here the belief update algorithmic scheme introduced in Section 2 for a system under deterioration dependence among components. All necessary implementation steps are listed in Algorithm 2. Bayesian inference is firstly conducted for the conditional beliefs $b\left(s_{t+1}|\alpha \right)$ and hyperparameters $b\left(\alpha \right)$, propagating uncertainty according to the transition model $p(s_{t+1}|s_t,a_t)$ and observation model $p(o_{t+1}|s_{t+1},a_t)$. The likelihood of collecting an observation given the hyperparameter(s) $p\left(o_{t+1}|\alpha \right)$, later necessary to update $b\left(\alpha \right)$, can be easily computed by marginalizing out the states other than $\alpha$:

Bayesian inference is then conducted for the hyperparameters:

After updating the conditional beliefs and common influencing variables, the marginal deterioration beliefs can be computed by marginalizing out the hyperparameters $\alpha$ as:

The effect of maintenance actions on the Gaussian dependence structure has not been explored in the existing literature Luque and Straub (2019, 2016), up to the best knowledge of the authors. Whereas the defined deterioration dependence is preserved if no maintenance interventions are planned, structural interventions can potentially disrupt the underlying correlation structure. For instance, if a structural system is specified with a correlated initial crack size among fatigue hotspots, this correlation structure will be perturbed after a component is repaired, along with the correlation reduction naturally experienced by the system over time. The correlation evolution associated with the latter is intrinsically quantified through the uncertainty propagation and updating operations formulated in Eqs. 10-11, whereas the correlation disruption associated with the former can be modeled by now defining the transition model $p(s_{t+1}|s_t,a_t,\alpha )$ of the involved components also conditional to the hyperparameter(s) $\alpha$, enabling therefore the removal or modification of the deterioration dependence by redefining and implementing the relevant correlation coefficients ${\lambda }_i$ in Eqs. 7 - 9. Additional discussion and implementation of this aspect is presented in the numerical experiments section.

As input to the I&M decision-making problem (Section 2), utilities, $r_{ins}$ and $r_{rep}$, are assigned to inspection and repair actions, respectively, specified according to available options and settings in each problem. The annual risk of a system failure, $r_F$, is defined accounting for two consecutive time steps (e.g., years) and its associated system failure cost, $r_f$, as:

The system structural failure event, as illustrated in Figure 1 by the node $F_{sys}$, is specified by a binary variable $p_{F_{sys}}$, indicating the failure and survival states, conditional on the belief state $b\left(s\right)$ of the structural components, as these are determined by the performed I&M actions in time. In principle, $p_{F_{sys}}$ could be directly defined as a function of the components belief state; in practice, however, $p_{F_{sys}}$ remains only conditional to the event of component failures $p_F$, specified as:

where $S_F$ corresponds to the components state subset classified as failure $S_F\subseteq S$. Within the deep reinforcement learning approach presented in Section 4, $P_{F_{sys}}$ can be computed via closed-form procedures and/or supported by efficient matricial algorithms Song and Kang (2009); or it can be computed following a general scheme, obtaining $p_{F_{sys}}$ through a simulator Der Kiureghian (2022). By assigning utilities to the system state, the importance of each structural element to the global risk of a system failure is implicitly accounted. To illustrate the effect of defining the failure risks at the system level, I&M strategies for a redundant 2-dimensional frame structure are later explored in Section 5.

In most structural systems, from bridges to offshore platforms or wind farms, inspection and repair actions are not planned separately for each structural element. Maintenance campaigns are instead scheduled, collecting information or performing repairs on a group of structural components. The cost model can thus be adapted from Eq. 3, to include a fixed campaign cost, $r_{camp}$, incurred every time a campaign is planned, along with inspection, $r_{ins}^{\left(l\right)}$, and repair, $r_{rep}^{\left(l\right)}$, costs assigned to the individual components according to any nonlinear function of choice, $H(.)$, e.g., simple linear sum operator, as:

The components explored throughout the numerical investigations are assumed to be exposed to a similar fatigue deterioration, described according to the Markovian model, originally proposed in Ditlevsen and Madsen (2007):

where the crack depth, $d$, evolution over time, $t$, follows a linear-elastic fracture mechanics law with material parameters $C_{FM}$ and $m$, stress range $S_R$, and $n$ annual stress cycles. At the component level, failure occurs if the crack depth, $d$, exceeds a critical size, $d_c$, that corresponds to the plate thickness. In a stochastic environment, the initial crack depth, $d_0$, along with fracture mechanics model parameters are either represented by random variables or deterministic parameters as listed in Table 1.

The fatigue deterioration is encoded in a deterioration rate DBN model, and ultimately shaping a factored POMDP, as shown on the left side of Fig. 2, and presented in Section 2. The continuous crack depth, $d$, is adequately discretizated into $|S_d|=30$ states conditional on $|S_{\tau }|=31$ fully observable deterioration rates states. The intervals and state space utilized for this deterioration rate model are listed in Table 2.

In terms of observation model, the inspection quality is quantified with a Probability of Detection curve $PoD\left(d\right)\sim Exp[\mu =8]$. Further details on the fatigue deterioration or observation model, including an extensive investigation of the discretization scheme can be found in Morato et al. (2022).

The system explored in this application is composed of ten components, each of which is subjected to a non-stationary fatigue deterioration, as described earlier in this Section. The system is assumed to be functional if at least 9-out-of-10 components are operational (not failed), thus characterized with a single step change in terms of redundancy with respect to a series system, which would correspond to the case of a 10-out-of-10 system. The system failure probability $p_{F_{sys}}$ is efficiently computed here, as a function of the failure state of all components, by following the recursive method proposed in Barlow and Heidtmann (1984).

A total of eight I&M planning scenarios are investigated, exploring different deterioration, risk, and cost dependencies among components. In terms of deterioration dependence, some environments are specified with an equally correlated initial crack size, $d_0$, among components, defined by an equal correlation ${\rho }_{eq}=0$, ${\rho }_{eq}=0.4$ and ${\rho }_{eq}=0.8$, respectively. Additionally, a deterioration environment is examined with an unequally correlated ${\rho }_{uq}$ initial crack size, $d_0$, among components. The unequal deterioration dependence case is originally specified with a different correlation among components of either $\rho =0.4$, $\rho =0.6$, or $\rho =0.8$, as shown on the left side of Fig. 3. After a Gaussian hierarchical structure with two hyperparameters is optimized, by computing the $\lambda$ parameters with the objective of satisfying Eq. 6, an approximated correlation structure is obtained with relatively small errors, as shown on the right side of Fig. 3. The approximated correlation structure with two hyperparameters is deemed to be sufficiently accurate for the conducted experiments. Otherwise, a more accurate correlation structure can be achieved by adding more hyperparameters, at the expense of additional computational cost, as explained in Section 3. For each of the aforementioned environments, specified with different deterioration dependencies, two I&M cost models are further investigated, i.e., an I&M cost model that incurs inspection and repair costs individually; and an I&M cost model in which an initial campaign cost exists, if at least one component is inspected or repaired, plus a cost surplus per inspected or repaired component.

Since each component, herein denoted as fatigue hotspot, contains 930 states, defined by the joint space of 30 crack states $S_d$ and 31 deterioration rate states $S_{\tau }$, these are a total of 9,300 states in the experiments that do not consider deterioration correlation (${\rho }_{eq}=0$). For experiments under equal correlation (${\rho }_{eq}$) the total states become 744,080, rising up to $5.95\cdot {10}^7$ for deterioration environments under unequal correlation (${\rho }_{uq}$). The increase of the state space corresponds to the incorporation of the Gaussian hierarchical model, in which crack and deterioration rate states are formulated conditional on the hyperparameter(s) states. When the deterioration correlation is modeled equally for all components, only one hyperparameter is sufficient to satisfy Eq. 6, while two hyperparameters are added for the case of unequal correlation, as explained earlier. Each hyperparameter is discretized into 80 states, initially prescribed with equal probability for each state. Note here the importance of optimizing the number of hyperparameters included in the model, as the state space grows exponentially with the number of considered random variables. By formulating the POMDPs’ transition model as dynamic Bayesian networks, the dimensionality is reduced from $|S_d{|}^2|S_{\tau }{|}^2$, in a flat structure, to $|S_d{|}^2|S_{\tau }|+|S_{\tau }{|}^2$ for the uncorrelated scenario. In that case, the transition model of only one component is reduced from 864,900 to 28,861 elements. Moreover, the formulation of the environment through a hierarchical deterioration dependence model importantly enables the decoupling of the joint state space at the system level, which would grow exponentially for a flat POMDP structure, but instead grows linearly now. For instance, in the setting under unequal deterioration dependence, the joint space would be described by $\{|S_d||S_{\tau }|{\}}^{N_C}$, equaling ${930}^{10}$ states, while it is now instead defined by $\{|S_d||S_{\tau }|N_C|S_{\alpha }||S_{\beta }|+|S_{\alpha }|+|S_{\beta }|\}$, thus resulting in $930\cdot 10\cdot {80}^2+80\cdot 2\simeq 5.95\cdot {10}^7$ states in the hierarchical model, with two hyperparameters ($\alpha$ and $\beta$) discretized into 80 states.

In terms of the neural networks’ architecture, DDMAC is laid out in this application with two hidden fully-connected layers of 100 neurons for each actor network, and two hidden fully-connected layers with 200 neurons for the critic network. The learning rate is adjusted during the training of the networks from ${10}^{-4}$ to ${10}^{-5}$ for the actor, and from ${10}^{-3}$ to ${10}^{-4}$ for the critic. The exploration is set up initially with a 100% random noise, decreasing linearly over the first 20,000 episodes to a random noise of 1%, held constant for the remaining episodes. A more stable and efficient training was found when a prioritization of do-nothing actions is implemented at the beginning of the training, because this allows visitation of more states, thus better exploration.

Following typical fatigue I&M planning settings, inspection and repair decisions are combined into three available actions per component: do-nothing / no-inspection, do-nothing / inspection, and perfect repair / no-inspection. The action perfect-repair / inspection is considered a priori suboptimal, without loss of generality, as it would be unusual to plan an inspection just after a component returns to its initial state. Inspections provide binary indications, i.e., detection or no-detection of a crack according to the observation model. In terms of costs, two different scenarios are considered. In the first case, inspection and repair costs are incurred independently per component, i.e., $r_{ins}=-1$ and $r_{rep}=-20$, respectively. In the second case, a campaign cost of $r_{camp}=-5$ is incurred if at least one component is inspected or repaired, plus a surplus per inspected or repair component of $r_{ins}=-0.2$ and $r_{per}=-20$ money units, respectively. The consequence of a system structural failure is $r_F=-10,000$ money units for both cases, and the discount factor $\gamma$ is 0.95 in all the experiments.