2.1. On-line particle filter

Particle filter (PF) methods, also referred to as sequential Monte Carlo (SMC) methods, are importance sampling-based techniques that use a set of weighted samples $ \left\{\left({\boldsymbol{\theta}}_n^{(i)},{w}_n^{(i)}\right):i=1,\dots, {N}_{\mathrm{par}}\right\} $ , called particles, to represent the posterior distribution of interest at estimation time step $ n $ , $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n}\right) $ . PFs form the following approximation to the posterior distribution of interest:

(4) $$ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n}\right)\approx \sum \limits_{i=1}^{N_{\mathrm{par}}}{w}_n^{(i)}\delta \left(\boldsymbol{\theta} -{\boldsymbol{\theta}}_n^{(i)}\right), $$

where $ \delta $ denotes the Dirac delta function.

When a new measurement $ {y}_n $ becomes available, PFs shift from $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n-1}\right) $ to $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n}\right) $ by importance sampling using an appropriately chosen importance distribution, which results in a reweighting procedure (updating of the weights). An important issue that arises from this weight-updating procedure is the sample degeneracy problem (Särkkä, Reference Särkkä2013). This relates to the fact that the importance weights $ {w}_n^{(i)} $ become more unevenly distributed with each updating step. In most cases, after a certain number of updating steps, the weights of almost all the particles assume values close to zero (see Figure 1). This problem is alleviated by the use of adaptive resampling procedures based on the effective sample size $ {N}_{\mathrm{eff}}=1/{\sum}_{i=1}^{N_{\mathrm{par}}}{\left({w}_n^{(i)}\right)}^2 $ (Liu and Chen, Reference Liu and Chen1998). Most commonly, resampling is performed with replacement according to the particle weights whenever $ {N}_{\mathrm{eff}} $ drops below a user-defined threshold $ {N}_{\mathrm{T}}={cN}_{\mathrm{par}},c\in \left[0,1\right] $ . Resampling introduces additional variance to the parameter estimates (Särkkä, Reference Särkkä2013). In the version of the PF algorithm presented in Algorithm 1, the dynamic model of equation (3) is used as the importance distribution, as originally proposed in the bootstrap filter by Gordon et al. (Reference Gordon, Salmond and Smith1993). Theoretical analysis of PF algorithms can be found in numerous seminal sources (e.g., Chopin, Reference Chopin2004; Del Moral et al., Reference Del Moral, Doucet and Jasra2006; Hu et al., Reference Hu, Schon and Ljung2008).

Figure 1. Sample degeneracy and impoverishment.

When using PFs to estimate time-invariant parameters, for which the process noise in the dynamic equation is zero, one runs into the issue of sample impoverishment (Särkkä, Reference Särkkä2013). The origin of this issue is the resampling process. More specifically, after a few resampling steps, most (or in extreme cases all) of the particles in the sample set end up assuming the exact same value, that is, the particle set consists of only few (or one) distinct particles (see Figure 1). The sample impoverishment issue poses the greatest obstacle for time-invariant parameter estimation with PFs. A multitude of techniques has been suggested in literature to alleviate the sample impoverishment issue in joint state-parameter estimation setups (see, e.g., Gilks and Berzuini, Reference Gilks and Berzuini2001; Liu and West, Reference Liu and West2001; Musso et al., Reference Musso, Oudjane and Le Gland2001; Storvik, Reference Storvik2002; Andrieu et al., Reference Andrieu, Doucet, Singh and Tadic2004; Andrieu et al., Reference Andrieu, Doucet and Holenstein2010; Carvalho et al., Reference Carvalho, Johannes, Lopes and Polson2010; Chatzi and Smyth, Reference Chatzi and Smyth2013; Chopin et al., Reference Chopin, Jacob and Papaspiliopoulos2013). Fewer works have proposed solutions for resolving this issue in parameter estimation setups (see, e.g., Chopin, Reference Chopin2002; Del Moral et al., Reference Del Moral, Doucet and Jasra2006). One of the simplest and most commonly used approaches consists in introducing artificial dynamics in the dynamic model of the parameter vector, that is, the dynamic model $ {\boldsymbol{\theta}}_n={\boldsymbol{\theta}}_{n-1}+{\boldsymbol{\varepsilon}}_{n-1} $ is employed, where $ {\boldsymbol{\varepsilon}}_{n-1} $ is a small artificial process noise (Kitagawa, Reference Kitagawa1998). In this way, the time-invariant parameter vector is transformed into a time-variant one, therefore, the parameter estimation problem deviates from the original one (Särkkä, Reference Särkkä2013; Kantas et al., Reference Kantas, Doucet, Singh and Chopin2015). This approach can introduce a bias and an artificial variance inflation in the estimates (Kantas et al., Reference Kantas, Doucet, Singh and Chopin2015). For these reasons, this approach is not considered in this article.

To resolve the sample impoverishment issue encountered when using the PF Algorithm 1 for parameter estimation, this work employs the particle filter with Gaussian mixture resampling (PFGM), described in Algorithm 2. The PFGM algorithm relates to preexisting concepts (van der Merwe and Wan, Reference van der Merwe and Wan2003; Veettil and Chakravorty, Reference Veettil and Chakravorty2016), and is here specifically tailored for the parameter estimation task, with its main goal being, in contrast to previous works, the quantification of the full posterior parameter uncertainty. A comparison between Algorithms 1 and 2 shows that the only difference lies in the way that the resampling step is performed. PFGM replaces the standard resampling process of PF by first approximating the posterior distribution at estimation step $ n $ by a Gaussian mixture model (GMM), which is fitted via the Expectation–Maximization (EM) algorithm (McLachlan and Krishnan, Reference McLachlan and Krishnan2007; Chen et al., Reference Chen, Gupta, Chen and Gupta2010) on the weighted particle set. The degenerating particle set is then rejuvenated by sampling $ {N}_{\mathrm{par}} $ new particles from the GMM of equation (5),

(5) $$ p(\boldsymbol{\theta} |{\mathbf{y}}_{1:n})\approx \sum \limits_{i=1}^{N_{\mathrm{GM}}}{\phi}_i\mathcal{N}(\boldsymbol{\theta}; {\boldsymbol{\mu}}_i,{\boldsymbol{\Sigma}}_i), $$

where $ {\phi}_i $ represents the weight of the Gaussian component $ i $ , while $ {\boldsymbol{\mu}}_i $ and $ {\boldsymbol{\Sigma}}_i $ are the respective mean vector and covariance matrix. The number of Gaussians in the mixture $ {N}_{\mathrm{GM}} $ , has to be chosen in advance, or can be estimated by use of appropriate algorithms (Schubert et al., Reference Schubert, Sander, Ester, Kriegel and Xu2017; Celeux et al., Reference Celeux, Fruewirth-Schnatter and Robert2019; Geyer et al., Reference Geyer, Papaioannou and Straub2019). In the numerical investigations of Section 3, we set $ {N}_{\mathrm{GM}} $  = 8. We point out that the efficacy of PFGM strongly depends on the quality of the GMM posterior approximation. The reason for applying a GMM (and not a single Gaussian) is that the posterior distribution can deviate from the normal distribution, and can even be multimodal or heavy-tailed.

The simple reweighting procedure used in the on-line PFs is based on the premise that $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n-1}\right) $ and $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n}\right) $ are likely to be similar, that is, that the new measurement $ {y}_n $ will not cause a very large change in the posterior. However, when that is not the case, this simple reweighting procedure is bound to perform poorly, leading to very fast degeneration of the particle set. In cases where already the first measurement set $ {y}_1 $ is strongly informative relative to the prior, the PF is bound to strongly degenerate already in the first weight updating step (e.g., we observe this in the second case study of subsec:RF in the case of 10 sensors). To counteract this issue, in this article, we incorporate the idea of simulated annealing (enforcing tempering of the likelihood function) (Neal, Reference Neal2001) when needed within the on-line PFGM algorithm, which we term the tPFGM Algorithm 3. The tPFGM algorithm draws inspiration from previous works (Deutscher et al., Reference Deutscher, Blake and Reid2000; Gall et al., Reference Gall, Potthoff, Schnörr, Rosenhahn and Seidel2007), but is here tailored for the parameter estimation task, opting for the quantification of the full posterior parameter uncertainty. The algorithm operates as follows: At estimation time step $ n $ , before performing the reweighting operation, the algorithm first checks the updated effective sample size for indication of sample degeneracy. If no degeneracy is detected, tPFGM operates exactly like PFGM. When sample degeneracy occurs, tPFGM employs adaptive tempering of the likelihood $ L\left({y}_n|\boldsymbol{\theta} \right) $ of the new measurement $ {y}_n $ in order to “sequentially” sample from $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n-1}\right) $ to $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n}\right) $ by visiting a sequence of artificial intermediate posteriors, as defined by the tempered likelihood function $ {L}^q\left({y}_n|\boldsymbol{\theta} \right) $ . The tempering factor $ q $ takes values between 0 and 1. When $ q=0 $ , the new measurement $ {y}_n $ is neglected, while $ q=1 $ entails considering the whole likelihood function of $ {y}_n $ , thus reaching to $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n}\right) $ . The intermediate values of $ q $ are adaptively selected via solution of the optimization problem in line 11 of Algorithm 3, which ensures that the effective sample size does not drop below the threshold $ {N}_{\mathrm{T}} $ for the chosen $ q $ value. Naturally, use of tPFGM can trigger more resampling events than PFGM, as resampling can occur more than once within a time step $ n $ .

The PFGM and tPFGM filters rely entirely on the posterior approximation via a GMM for sampling $ {N}_{\mathrm{par}} $ new particles during the resampling process. However, there is no guarantee that these new particles follow the true posterior distribution of interest. The IBIS filter of the following Section 2.2 aims at addressing this issue.

2.2. Iterated batch importance sampling

Implementing MCMC steps within PF methods to move the particles after a resampling step was originally proposed by Gilks and Berzuini (Reference Gilks and Berzuini2001), in the so-called resample-move algorithm. Chopin (Reference Chopin2002) introduced a special case of the resample-move algorithm, specifically tailored for application to static parameter estimation purposes, namely the iterated batch importance sampling (IBIS) filter. IBIS was originally introduced as an iterative method for solving off-line estimation tasks by incorporating the sequence of measurements one at a time. In doing this, the algorithm visits the sequence of intermediate posteriors within its process, and can therefore also be used to perform on-line estimation tasks. An on-line version of the IBIS filter is presented in Algorithm 5, used in conjuction with the MCMC routine of Algorithm 4.

The core idea of the IBIS filter is the following: At estimation step $ n $ , if sample degeneracy is identified, first the particles are resampled with replacement, and subsequently the resampled particles are moved with a Markov chain transition kernel whose stationary distribution is $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n}\right) $ . More specifically, each of the $ {N}_{\mathrm{par}} $ resampled particles is used as the seed to perform a single MCMC step. This approach is inherently different to standard applications of MCMC, where a transition kernel is applied multiple times on one particle.

A question that arises is how to choose the Markov chain transition kernel. Chopin (Reference Chopin2002) argues for choosing a transition kernel that ensures that the proposed particle only weakly depends on the seed particle value. It is therefore recommended to use an independent Metropolis-Hastings (IMH) kernel, wherein the proposed particle is sampled from a proposal distribution $ g $ , which has to be as close as possible to the target distribution $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n}\right) $ . In obtaining such a proposal distribution, along the lines of what is described in Section 2.1, in this work we employ a GMM approximation (see equation (5)) of the target distribution as the proposal density $ {g}_{\mathrm{GM}}\left(\boldsymbol{\theta} \right) $ within the IMH kernel (Papaioannou et al., Reference Papaioannou, Papadimitriou and Straub2016; South et al., Reference South, Pettitt and Drovandi2019). The IMH kernel with a GMM proposal distribution is denoted IMH-GM herein. The acceptance probability (line 6 of Algorithm 4) of the IMH-GM kernel is a function of both the initial seed particle and the GMM proposed particle. The acceptance rate can indicate how efficient the IMH-GM kernel is in performing the MCMC move step within the IBIS algorithm. It is important to note that when computing the acceptance probability, a call of the full likelihood function is invoked, which requires the whole set of measurements $ {\mathbf{y}}_{1:n} $ to be processed; this leads to a significant additional computational demand, which pure on-line methods are not supposed to accommodate (Doucet et al., Reference Doucet, de Freitas and Gordon2001).

The performance of the IBIS sampler depends highly on the mixing properties of the IMH-GM kernel. If the kernel leads to slowly decreasing chain auto-correlation, the moved particles are bound to remain in regions close to the particles obtained by the resampling step. This can lead to an underrepresentation of the parameter space of the intermediate posterior distribution. It might therefore be beneficial to add a burn-in period within the IMH-GM kernel (Del Moral et al., Reference Del Moral, Doucet and Jasra2006). Implementing that is straightforward and is shown in Algorithm 4, where $ {n}_{\mathrm{B}} $ is the user-defined number of burn-in steps. Naturally, the computational cost of the IMH-GM routine increases linearly with the number of burn-in steps.

Algorithm 5 details the workings of the IMH-GM-based IBIS filter used in this work. In line 11 of this algorithm, the IMH-GM routine of Algorithm 4 is called, which implements the IMH-GM kernel for the MCMC move step. Comparing Algorithms 2 and 5, it is clear that both filters can be used for on-line inference within a single run, but the IBIS filter has significantly larger computational cost, as will also be demonstrated in the numerical investigations of Section 3. In the same spirit as the proposed tPFGM Algorithm 3, which enforces simulated annealing (tempering of the likelihood function) in cases when $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n-1}\right) $ and $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:n}\right) $ are likely to be quite different, the same idea can be implemented also within the IBIS algorithm. That leads to what we refer to as the tIBIS algorithm in this article.

2.3. Off-line sequential Monte Carlo sampler

In Section 4 of Del Moral et al. (Reference Del Moral, Doucet and Jasra2006), the authors presented a generic approach to convert an off-line MCMC sampler into a sequential Monte Carlo (SMC) sampler tailored for performing off-line estimation tasks, that is, for estimating the single posterior density of interest $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:N}\right) $ . The off-line SMC sampler used in this work is presented in Algorithm 6 based on Del Moral et al. (Reference Del Moral, Doucet and Jasra2006) and Jasra et al. (Reference Jasra, Stephens, Doucet and Tsagaris2011). The key idea of this sampler is to adaptively construct the following artificial sequence of densities,

(6) $$ {\pi}_j\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:N}\right)\propto {L}^{q_j}\left({\mathbf{y}}_{1:N}|\boldsymbol{\theta} \right){\pi}_{\mathrm{pr}}\left(\boldsymbol{\theta} \right), $$

where $ {q}_j $ is a tempering parameter which obtains values between 0 and 1, in order to “sequentially” sample in a smooth manner from the prior to the final single posterior density of interest. Once $ {q}_j=1 $ , $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:N}\right) $ is reached. Similar to what was described in tPFGM, the intermediate values of $ {q}_j $ are adaptively found via solution of the optimization problem in line 5 of Algorithm 6. The GMM approximation of the intermediate posteriors and the IMH-GM kernel of Algorithm 4 in order to move the particles after resampling are also key ingredients of this SMC sampler. Unlike PFGM and IBIS, this SMC algorithm cannot provide the on-line solution within a single run, and has to be rerun from scratch for every new target posterior of interest. In this regard, use of Algorithm 6 for on-line inference is impractical. We choose to include this algorithm for the purpose of the comparative assessment in this article, as off-line algorithms are generally considered to be better suited for the full posterior uncertainty quantification of time-invariant parameters.

2.4. Computational remarks

The algebraic operations in all presented algorithms are implemented in the logarithmic scale, which employs evaluations of the logarithm of the likelihood function and, hence, ensures computational stability. Furthermore, the EM step for fitting the GMM is performed after initially transforming the prior joint probability density function of $ \boldsymbol{\theta} $ to an underlying vector $ \mathbf{u} $ of independent standard normal random variables (Kiureghian and Liu, Reference Kiureghian and Liu1986). In standard normal space, the parameters are decorrelated, which enhances the performance of the EM algorithm.

3.1. Low-dimensional case study: Paris–Erdogan fatigue crack growth model

A fracture mechanics-based model serves as the first case study. This describes the fatigue crack growth evolution under increasing stress cycles (Paris and Erdogan, Reference Paris and Erdogan1963; Ditlevsen and Madsen, Reference Ditlevsen and Madsen1996). The crack growth follows the following first-order differential equation (7), known as Paris–Erdogan law,

(7) $$ \frac{da\left(n,\boldsymbol{\theta} \right)}{dn}=\exp \left({C}_{\mathrm{ln}}\right){\left[\Delta S\sqrt{\pi a(n)}\right]}^m, $$

where $ \boldsymbol{\theta} =\left[a,\Delta S,{C}_{\mathrm{ln}},m\right] $ is a vector containing the uncertain time-invariant model parameters. Specifically, $ a $ is the crack length, $ n $ is the number of stress cycles, $ \Delta S $ is the stress range per cycle when assuming constant stress amplitudes, $ C $ and $ m $ represent empirically determined model parameters; $ {C}_{\mathrm{ln}} $ corresponds to the natural logarithm of $ C $ .

The solution to this differential equation, with boundary condition $ a\left(n=0\right)={a}_0 $ , can be written as a function of the number of stress cycles $ n $ and the vector $ \boldsymbol{\theta} =\left[{a}_0,\Delta S,{C}_{\mathrm{ln}},m\right] $ as (for the derivation see, e.g., Ditlevsen and Madsen (Reference Ditlevsen and Madsen1996)):

(8) $$ a(n,\boldsymbol{\theta} )={[(1-\frac{m}{2})\mathrm{exp}({C}_{\mathrm{ln}})\Delta {S}^m{\pi}^{m/2}n+{a}_0^{(1-m/2)}]}^{(1-m/2)^{-1}}. $$

We assume that noisy measurements of the crack $ {y}_n $ are obtained sequentially at different values of $ n $ . The measurement equation (9) assumes a multiplicative lognormal measurement error, $ \exp \left({\varepsilon}_n\right) $ .

(9) $$ {y}_n=a\left(n,\boldsymbol{\theta} \right)\exp \left({\varepsilon}_n\right). $$

In numerical investigations that follow, the measurement equation (9) is used for generating synthetic measurements of the deterioration state. In this context, a multiplicative lognormal measurement error ensures that nonnegative generated measurements of the deterioration state are not feasible.

Under this assumption in equation (9), the likelihood function for a measurement at a given $ n $ is shown in equation (10).

(10) $$ L\left({y}_n;a\left(n,\boldsymbol{\theta} \right)\right)=\frac{1}{\sigma_{\varepsilon_n}\sqrt{2\pi }}\exp \left[-\frac{1}{2}{\left(\frac{\ln \left({y}_n\right)-{\mu}_{\varepsilon_n}-\ln \left(a\left(n,\boldsymbol{\theta} \right)\right)}{\sigma_{\varepsilon_n}}\right)}^2\right]. $$

Table 1 shows the prior probability distribution model for each random variable in the vector $ \boldsymbol{\theta} $ (Ditlevsen and Madsen, Reference Ditlevsen and Madsen1996; Straub, Reference Straub2009), as well as the assumed probabilistic model of the measurement error. In this case study we are dealing with a nonlinear model and a parameter vector with non-Gaussian prior distribution.

Table 1. Prior distribution model for the fatigue crack growth model parameters and the measurement error.

3.1.1. Markovian state-space representation for application of on-line filters

A Markovian state-space representation of the deterioration process is required for application of on-line filters. The number of stress cycles is discretized as $ n=k\Delta n $ , with $ k=1,\dots, 100 $ denoting the estimation time step and $ \Delta n=1\times {10}^5 $ the number of stress cycles per time step. The dynamic and measurement equations of the discrete-time state-space representation of the fatigue crack growth model with unknown time-invariant parameters $ \boldsymbol{\theta} =\left[{a}_0,\Delta S,{C}_{\mathrm{ln}},m\right] $ are shown below:

(11) $$ {\displaystyle \begin{array}{ll}{\boldsymbol{\theta}}_k& \hskip-1em ={\boldsymbol{\theta}}_{k-1}\\ {}{y}_k& \hskip-1em =a\left(k\Delta n,{\boldsymbol{\theta}}_k\right)\exp \left({\varepsilon}_k\right)={\left[\left(1-\frac{m_k}{2}\right)\exp \left({C}_{\ln_k}\right)\Delta {S}_k^{m_k}{\pi}^{m_k/2}k\Delta n+{a}_0^{\left(1-{m}_k/2\right)}\right]}^{{\left(1-{m}_k/2\right)}^{-1}}\exp \left({\varepsilon}_k\right).\end{array}} $$

The state-space model of equation (11) is nonlinear and the prior is non-Gaussian. For reasons explained in Section 2, the subscript $ k $ in $ {\boldsymbol{\theta}}_k $ is dropped in the remainder of this section.

3.1.2. Reference posterior solution

For the purpose of performing a comparative assessment of the different filters, an underlying “true” realization of the fatigue crack growth process $ {a}^{\ast}\left(n,\boldsymbol{\theta} \right) $ is generated for $ n=k\Delta n $ , with $ k=1,\dots, 100 $ and $ \Delta n=1\times {10}^5 $ . This realization corresponds to the randomly generated “true” vector of time-invariant parameters $ {\boldsymbol{\theta}}^{\ast }=\left[{a}_0^{\ast }=2.0,\Delta {S}^{\ast }=50.0,{C}_{\mathrm{ln}}^{\ast }=-33.5,{m}^{\ast }=3.7\right] $ . Sequential synthetic crack monitoring measurements $ {y}_k $ are sampled from the measurement equation (9) for $ a\left(k\Delta n,{\boldsymbol{\theta}}^{\ast}\right) $ , and for randomly generated measurement noise samples $ \exp \left({\varepsilon}_k\right) $ . These measurements are scattered in green in Figure 3.

Based on the generated measurements, the sequence of reference posterior distributions $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:k}\right) $ is obtained using the prior distribution as an envelope distribution for rejection sampling (Smith and Gelfand, Reference Smith and Gelfand1992; Rubinstein and Kroese, Reference Rubinstein and Kroese2016). More specifically, for each of the 100 posterior distributions of interest $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:k}\right) $ , $ {10}^5 $ independent samples are generated. The results of this reference posterior estimation of the four time-invariant model parameters are plotted in Figure 2. With posterior samples, the reference filtered estimate of the crack length $ a\left(k\Delta n,\boldsymbol{\theta} \right) $ at each estimation step is also obtained via the model of equation (8) and plotted in Figure 3. In the left panel of this figure, the filtered state is plotted in logarithmic scale. In an off-line estimation, a single posterior density is of interest. One such reference posterior estimation result for the last estimation step, $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:100}\right) $ , is plotted for illustration in Figure 4.

Figure 2. Reference posterior solution: mean and credible intervals for the sequence of posterior distributions $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:k}\right) $ .

Figure 3. Reference mean and credible intervals for the filtered crack growth state $ a\left(k\Delta n,\boldsymbol{\theta} \right) $ .

Figure 4. Reference final posterior: prior and single posterior distribution of interest $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:100}\right) $ .

3.1.3. Comparative assessment of the investigated on-line and off-line filters

We apply the PFGM filter with 5,000 and 50,000 particles, the IBIS filter with 5,000 particles, and the SMC filter with 5,000 particles for performing on-line and off-line time-invariant parameter estimation tasks. We evaluate the performance of each filter by taking the relative error of the estimated mean and standard deviation of each of the four parameters with respect to the reference posterior solution. For example, the relative error in the estimation of the mean of parameter $ {a}_0 $ at a certain estimation step $ k $ is computed as $ |\frac{\mu_{a_0,k}-{\hat{\mu}}_{a_0,k}}{\mu_{a_0,k}}| $ , where $ {\mu}_{a_0,k} $ is the reference posterior mean from rejection sampling (Section 3.1.2), and $ {\hat{\mu}}_{a_0,k} $ is the posterior mean estimated with each filter. Each filter is run 50 times, and the mean relative error of the mean and the standard deviation of each parameter, together with the 90% credible intervals (CI), are obtained. These are plotted in Figure 5.

Figure 5. Comparison of the relative error of the mean and standard deviation of the parameters evaluated for each filter. The solid lines show the mean and the shaded areas the 90% credible intervals inferred from 50 repeated runs of each filter. In the horizontal axis, $ n $ is the number of stress cycles.

Figure 6 plots the $ {L}^2 $ relative error norm of the mean and the standard deviation of all four parameters, that is, the quantity of equation (12) (here formulated for the mean at estimation step $ k $ ).

(12) $$ \sqrt{\frac{\sum \limits_{i=1}^d{\left({\mu}_{i,k}-{\hat{\mu}}_{i,k}\right)}^2}{\sum \limits_{i=1}^d{\left({\mu}_{i,k}\right)}^2}}, $$

where $ d $ is the dimensionality of the time-invariant parameter vector $ \boldsymbol{\theta} $ (in this example $ d=4 $ ). More specifically, Figure 6 plots the mean and credible intervals of the $ {L}^2 $ relative error norm of the estimated mean and standard deviation, as obtained from 50 runs of each filter.

Figure 6. Comparison of the $ {L}^2 $ relative error norm of the mean and the standard deviation of the parameters evaluated for each filter. The solid lines show the mean and the shaded areas the 90% credible intervals inferred from 50 repeated runs of each filter. In the horizontal axis, $ n $ is the number of stress cycles.

Figures 5 and 6 reveal that, when all three filters are run with the same number of particles, the IBIS and SMC filters yield superior performance over PFGM. When the number of particles in the PFGM filter is increased to 50,000, the PFGM filter performance is comparable to the one of the IBIS and SMC filters. In estimating the mean, the mean $ {L}^2 $ relative error norm obtained from the PFGM filter with 50,000 particles is slightly larger than the corresponding error obtained from IBIS and SMC with 5,000 particles, while the 90% credible intervals of the PFGM filter estimation are still wider. In estimating the standard deviation, the PFGM filter with 50,000 particles proves competitive.

Figures 5 and 6 show the estimation accuracy of each filter when used for on-line inference, that is, for estimating the whole sequence of 100 posterior distributions $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:k}\right) $ , $ k=1,\dots, 100 $ . The PFGM and IBIS filters, being intrinsically on-line filters, provide the whole posterior sequence with one run. On the other hand, the off-line SMC filter is run anew for each of the 100 required posterior estimations. Hence, Figures 5 and 6 enclose the results of both the on-line and the off-line inference. If one is interested in the off-line estimation accuracy at a specific stress cycle $ n $ , one can simply consider a vertical “cut” at $ n $ .

Table 2 documents the computational cost associated with each filter, expressed in the form of required model evaluations induced by calls of the likelihood function. When using the term model, we refer to the model defined in equation (8), which corresponds to an analytical expression with negligible associated runtime. However, unlike the simple measurement equation that we have assumed in this example, in many realistic deterioration monitoring settings, the deterioration state cannot be measured directly (e.g., in vibration-based structural health monitoring Kamariotis et al., Reference Kamariotis, Chatzi and Straub2022). In such cases, each deterioration model evaluation often entails evaluation of a finite element (FE) model, which has substantial runtime. It, therefore, appears appropriate to evaluate the filters’ computational cost in terms of required model evaluations. The on-line PFGM filter with 5,000 particles requires $ 5\times {10}^5 $ model evaluations, and yields by far the smallest computational cost, while at the same time providing the solution to both on-line and off-line estimation tasks. However, it also yields the worst performance in terms of accuracy of the posterior estimates. Running the IBIS filter with 5,000 particles, which performs MCMC move steps, leads to $ 3.4\times {10}^6 $ model evaluations. Comparing this value against the $ 5\times {10}^5 $ model evaluations required by the PFGM filter with 5,000 particles for performing the same task distinctly shows the computational burden associated with MCMC move steps, which require a complete browsing of the whole measurement data set in estimating the acceptance probability. However, the IBIS filter also leads to enhanced estimation accuracy, which might prove significant when the subsequent tasks entail prognosis of the deterioration evolution, the structural reliability or the remaining useful lifetime, and eventually the predictive maintenance planning. Using 50,000 particles, the PFGM filter performance increases significantly with a computational cost that is comparable to the IBIS filter with 5,000 particles. For the off-line SMC algorithm, $ 4.5\times {10}^6 $ model evaluations are required only for the task of estimating the final posterior density. The $ 1.9\times {10}^8 $ model evaluations required by the SMC for obtaining the whole sequence of posteriors $ {\pi}_{\mathrm{pos}}\left(\boldsymbol{\theta} |{\mathbf{y}}_{1:k}\right) $ , $ k=1,\dots, 100 $ , clearly demonstrate that off-line MCMC techniques are unsuited to on-line estimation tasks.

Table 2. Average number of model evaluations (equation (8)) for the fatigue crack growth model parameter estimation.

3.2. High-dimensional case study: Corrosion deterioration spatially distributed across beam

As a second case study, we employ the deterioration model of equation (13), which describes the spatially and temporally varying corrosion deterioration across the structural beam shown in Figure 7.

(13) $$ D\left(t,x\right)=A(x){t}^{B(x)},\hskip2em t=0,\dots, 50. $$

Figure 7. Structural beam subjected to spatially and temporally varying corrosion deterioration. The deterioration process is monitored from sensors deployed at specific sensor locations (in green).

$ A(x),B(x),x\in \Omega $ are random fields defined on $ \Omega =\left[0,L\right] $ , with $ L $ denoting the length of the beam taken as $ L=4 $ m. $ A(x) $ models the deterioration rate, while $ B(x) $ models the nonlinearity effect of the deterioration process in terms of a power law in time. The corrosion deterioration $ D\left(t,x\right) $ is therefore also a spatial random field.

A random field, by definition, contains an infinite number of random variables, and must therefore be discretized (Vanmarcke, Reference Vanmarcke2010). One of the most common methods for discretization of random fields is the midpoint method (Der Kiureghian and Ke, Reference Der Kiureghian and Ke1988). Thereby the spatial domain $ \Omega $ is discretized into $ m $ elements, and the random field is approximated within each element through the random variable that corresponds to midpoint of the element. In that case, the uncertain time-invariant deterioration model parameter vector is $ \boldsymbol{\theta} =\left[{A}_1,\dots, {A}_m,{B}_1,\dots, {B}_m\right] $ , where are the random variables corresponding to the element midpoints. With the midpoint discretization, the spatial deterioration $ D\left(t,x\right) $ is parametrized by $ \boldsymbol{\theta} $ .

We assume that noisy measurements of the corrosion deterioration state $ {D}_{t,l}\left(\boldsymbol{\theta} \right):= D\left(t,{x}_l,\boldsymbol{\theta} \right) $ at time $ t $ and at certain locations $ l $ of the beam are obtained sequentially (summarized in one measurement per year) from $ {n}_l $ sensors deployed at these locations ( $ {n}_l=10 $ sensor locations are shown in Figure 7). The measurement equation (14), describing the corrosion measurement at time $ t $ and sensor location $ l $ , assumes a multiplicative measurement error, $ \exp \left({\varepsilon}_{t,l}\right) $ .

(14) $$ {y}_{t,l}={D}_{t,l}\left(\boldsymbol{\theta} \right)\exp \left({\varepsilon}_{t,l}\right)={A}_{i_l}{t}^{B_{i_l}}\exp \left({\varepsilon}_{t,l}\right), $$

where $ {i}_l $ returns the discrete element number of the midpoint discretization within which the measurement location $ l $ lies. Table 3 shows the prior distribution model for the two random fields of the deterioration model of equation (13) and the assumed probabilistic model of the multiplicative measurement error. Since $ A(x) $ models a lognormal random field, $ \ln \left(A(x)\right) $ follows the normal distribution. For both random fields $ \ln \left(A(x)\right) $ and $ B(x) $ , the exponential correlation model with correlation length of 2 m is applied (Sudret and Der Kiureghian, Reference Sudret and Der Kiureghian2000).

Table 3. Prior distribution model for the corrosion deterioration model parameters and the measurement error.

The goal is to update the time-invariant deterioration model parameters $ \boldsymbol{\theta} =\left[{A}_1,\dots, {A}_m,{B}_1,\dots, {B}_m\right] $ given sequential noisy corrosion measurements $ {y}_{t,l} $ from $ {n}_l $ deployed sensors. The dimensionality of the problem is $ d=2\times m $ . Hence, the more elements in the midpoint discretization, the higher the dimensionality of the parameter vector.

The main goal of this second case study is to investigate the effect of the problem dimensionality and the amount of sensor information on the posterior results obtained with each filter. We choose the following three midpoint discretization schemes:

1. 1. $ m=25 $ elements: $ d=50 $ time-invariant parameters to estimate.

2. 2. $ m=50 $ elements: $ d=100 $ time-invariant parameters to estimate.

3. 3. $ m=100 $ elements: $ d=200 $ time-invariant parameters to estimate.

Furthermore, we choose the following three potential sensor arrangements:

1. 1. $ {n}_l=2 $ sensors (the $ 4\mathrm{th} $ and $ 7\mathrm{th} $ sensors of Figure 7).

2. 2. $ {n}_l=4 $ sensors (the $ 1\mathrm{st},4\mathrm{th},7\mathrm{th} $ , and $ 10\mathrm{th} $ sensors of Figure 7).

3. 3. $ {n}_l=10 $ sensors of Figure 7.

We, therefore, study nine different cases of varying problem dimensionality and the number of sensors.

3.2.1. Markovian state-space representation for application of on-line filters

A Markovian state-space representation of the deterioration process is required for application of on-line filters. The dynamic and measurement equations are shown in equation (15). The measurement equation is written in the logarithmic scale. Time $ t $ is discretized in yearly estimation time steps $ k $ , that is, $ k=1,\dots, 50 $ , and the subscript $ l=1,\dots, {n}_l $ corresponds to the sensor location.

(15) $$ {\displaystyle \begin{array}{ll}\hskip3em {\boldsymbol{\theta}}_k& \hskip-1em ={\boldsymbol{\theta}}_{k-1}\\ {}\ln \left({y}_{k,l}\right)& \hskip-1em =\ln \left({D}_{k,l}\left(\;{\boldsymbol{\theta}}_k\right)\right)+{\varepsilon}_{k,l}\Rightarrow \ln \left({y}_{k,l}\right)=\ln \left({A}_{k,{i}_l}\right)+{B}_{k,{i}_l}\ln (k)+{\varepsilon}_{k,l},\end{array}} $$

where $ {\boldsymbol{\theta}}_k $ denotes the time-invariant parameter vector at time step $ k $ . In the logarithmic scale, both the dynamic and measurement equations are linear functions of Gaussian random variables. For reasons explained in Section 2, the subscript $ k $ in $ {\boldsymbol{\theta}}_k $ is dropped in the following.

3.2.2. Underlying “true” realization

To generate a high-resolution underlying “true” realization of the two random fields $ A(x) $ and $ B(x) $ , and the corresponding synthetic monitoring data set, we employ the Karhunen–Loève (KL) expansion using the first 400 KL modes. The KL expansion is an alternative random field discretization scheme to the midpoint method that represents the random field in terms of the eigenfunctions of its autocovariance function (Ghanem and Spanos, Reference Ghanem and Spanos1991; Sudret and Der Kiureghian, Reference Sudret and Der Kiureghian2000). The implementation of the KL expansion can be found in the Matlab codes accompanying this article.Footnote ² We remark that a fine resolution KL expansion is chosen to represent the “true” realization in order to avoid the inverse crime (Wirgin, Reference Wirgin2004). These realizations are shown in the left panel of Figure 8. Given these $ A(x) $ and $ B(x) $ realizations, the underlying “true” realizations of the deterioration process at 10 specific beam locations are generated, which correspond to the 10 potential sensor placement locations shown in Figure 7. Subsequently, a synthetic corrosion sensor measurement data set (one measurement per year) at these 10 locations is generated from the measurement equation (14). These are shown in the right panel of Figure 8.

Figure 8. Left: the blue solid line plots the underlying “true” realization of $ \ln \left(A(x)\right) $ and $ B(x) $ created using the KL expansion. Right: the blue solid line plots the underlying “true” realization of $ \hskip0.1em \ln \left(D\left(t,x\right)\right) $ at 10 specific sensor locations and the corresponding synthetic sensor monitoring data are scattered in black. In both figures, the black dashed lines plot the prior mean and the black solid lines the prior 90% credible intervals.

3.2.3. Reference posterior solution

For the investigated linear Gaussian state space representation of equation (15), we create reference on-line posterior solutions for each of the nine considered cases by applying the Kalman filter (KF) (Kalman, Reference Kalman1960), which is the closed-form solution to the Bayesian filtering equations. The process noise covariance matrix in the KF equations is set equal to zero. The linear Gaussian nature of the chosen problem ensures existence of an analytical reference posterior solution obtained with the KF. One such reference on-line posterior solution for the case described by $ m=25 $ elements ( $ d=50 $ ) and $ {n}_l=4 $ sensors is shown in Figure 9.

Figure 9. Case with $ m=25,{n}_l=4 $ : reference on-line posterior solution at 10 locations across the beam obtained by applying the Kalman filter for solving equation (15). The solid blue horizontal line represents the underlying “true” values of $ \ln \left(A(x)\right) $ and $ B(x) $ at these locations. The black dashed lines plot the posterior mean and the black solid lines the posterior 90% credible intervals. Locations 1, 4, 7, and 10 correspond to the four assumed sensor placement locations.

3.2.4. Comparative assessment of the investigated on-line and off-line filters

The goal of this section is to offer a comparative assessment among the three Bayesian filtering algorithms presented in this article when applied on a high-dimensional problem. To be able to derive a reference solution, as described above, a linear Gaussian state space representation of a structural deterioration problem has been defined (equation (15)). We apply the tPFGM filter, the tIBIS filter, and the SMC filter, all with $ {N}_{\mathrm{par}} $  = 2000 particles, for estimating the time-invariant parameter vector $ \boldsymbol{\theta} $ . For each of the nine cases of varying problem dimensionality and number of sensors described above, we compute the $ {L}^2 $ relative error norm of the estimated means, correlation coefficients, and standard deviations of the parameters with respect to the corresponding KF reference posterior solution, that is, we estimate a quantity as in equation (12) for all estimation steps $ k=1,\dots, 50 $ . In Figures 10–12, we plot the mean and credible intervals of these relative errors as obtained from 50 different runs. The off-line SMC filter, which does not provide the on-line solution within a single run, is run anew for estimating the single posterior density of interest at years 10, 20, 30, 40, and 50, and in between, the relative error is linearly interpolated. Although each of the nine panels in the figures corresponds to a different case with a different underlying KF reference solution, their y-axes have the same scaling. Table 4 documents the computational cost of each filter in each considered case, measured by average number of evaluations of the model of equation (13).

Figure 10. Comparison of the $ {L}^2 $ relative error norm of the means of the parameters evaluated for each filter. The solid lines show the mean and the shaded areas the 90% credible intervals inferred from 50 repeated runs of each filter.

Table 4. Average number of model evaluations (equation (13)) for the high-dimensional case study.

Note. For the SMC, the required model evaluations for obtaining the single final posterior density are reported.

Figures 10 and 11 show that the off-line IMH-GM-based SMC filter yields the best performance in estimating the KF reference posterior mean and correlation, for all nine considered cases, while at the same time producing the narrowest credible intervals. Comparison of the relative errors obtained with the SMC and tIBIS filters reveals that, although they are both reliant on the IMH-GM MCMC move step, the on-line tIBIS filter leads to larger estimation errors. The on-line tPFGM and tIBIS filters generate quite similar results in estimating the reference posterior mean and correlation, thus rendering the benefit of the MCMC move step in tIBIS unclear, except in cases with more sensors and lower parameter dimension. Figures 10 and 11 reveal a slight trend, indicating that for fixed dimensionality, the availability of more sensors (i.e., stronger information content in the likelihood function) leads to a slight decrease in the relative errors when using the SMC and tIBIS filters, whereas the opposite trend can be identified for the tPFGM filter. Increasing problem dimensionality (for fixed number of sensors) does not appear to have strong influence on the posterior results in any of the columns of Figures 10–12, a result that initially appears puzzling.

Figure 11. Comparison of the $ {L}^2 $ relative error norm of the correlation coefficients of the parameters evaluated for each filter. The solid lines show the mean and the shaded areas the 90% credible intervals inferred from 50 repeated runs of each filter.

Figure 12. Comparison of the $ {L}^2 $ relative error norm of the standard deviations of the parameters evaluated for each filter. The solid lines show the mean and the shaded areas the 90% credible intervals inferred from 50 repeated runs of each filter.

Figure 12 conveys that the tPFGM filter, which entirely depends on the GMM posterior approximation, induces the smallest relative errors for the estimation of the standard deviation of the parameters in all considered cases. This result reveals a potential inadequacy of the single application of the IHM-GM kernel for the move step within the tIBIS and SMC filters in properly exploring the space of $ \boldsymbol{\theta} $ . In all 50 runs of the tIBIS and SMC filters, the standard deviation of the parameters is consistently underestimated compared to the reference, unlike when applying the tPFGM filter.

Based on the discussion of Section 2.2, we introduce a burn-in period of $ {n}_B $  = 5 in the IMH-GM kernel of Algorithm 4 and perform 50 new runs of the tIBIS and SMC filters. One can expect that inclusion of a burn-in is more likely to ensure sufficient exploration of the intermediate posterior distributions. However, at the same time, the computational cost of tIBIS and SMC increases significantly, with a much larger number of required model evaluations than in Table 4. In Figures 13 and 14, we plot the mean and credible intervals for the relative errors in the estimation of the mean and standard deviation of the parameters. Comparing Figures 10 and 13, inclusion of burn-in is shown to lead to an improved performance of tIBIS and SMC in estimating the mean of the parameters in all cases. This improvement is more evident in the lower-dimensional case with 25 elements, and lessens as the problem dimension increases. Hence, it is only after the inclusion of burn-in, which leads to an enhanced posterior solution, that one starts observing the anticipated deterioration of the tIBIS and SMC filters’ performance with increasing dimensionality. This effect was masked in the results of Figure 10 without burn-in. This point becomes more evident when looking at the relative errors of the estimated standard deviation in Figure 14. With burn-in, the tIBIS and SMC filters provide better results than the tPFGM filter in estimating the standard deviation in the case of 25 elements, but perform progressively worse as the dimensionality increases, where they underestimate the KF reference standard deviation. This underestimation is clearly illustrated in Figure 15. The reason for this behavior is the poor performance of the IMH-GM algorithm in high dimensions, which is numerically demonstrated in Papaioannou et al. (Reference Papaioannou, Papadimitriou and Straub2016). We suspect that this behavior is related to the degeneracy of the acceptance probability of MH samplers in high dimensions, which has been extensively discussed in the literature for random walk samplers, for example, in Gelman et al. (Reference Gelman, Gilks and Roberts1997), Au and Beck (Reference Au and Beck2001), Katafygiotis and Zuev (Reference Katafygiotis and Zuev2008), Beskos and Stuart (Reference Beskos, Stuart, L’ Ecuyer and Owen2009), Cotter et al. (Reference Cotter, Roberts, Stuart and White2013), and Papaioannou et al. (Reference Papaioannou, Betz, Zwirglmaier and Straub2015). Single application of the IHM-GM kernel without burn-in yielded acceptance rates of around 50% for all cases. With inclusion of burn-in, in higher dimensions, the acceptance rate in IMH-GM drops significantly in the later burn-in steps, leading to rejection of most proposed particles. To alleviate this issue, one could consider using the preconditioned Crank Nicolson (pCN) sampler to perform the move step within the IBIS and SMC filters, whose performance is shown to be independent of the dimension of the parameter space when the prior is Gaussian (Cotter et al., Reference Cotter, Roberts, Stuart and White2013).

Figure 13. Comparison of the $ {L}^2 $ relative error norm of the mean of the parameters evaluated for each filter. The solid lines show the mean and the shaded areas the 90% credible intervals inferred from 50 repeated runs of each filter. Burn-in $ {n}_B $  = 5.

Figure 14. Comparison of the $ {L}^2 $ relative error norm of the standard deviation of the parameters evaluated for each filter. The solid lines show the mean and the shaded areas the 90% credible intervals inferred from 50 repeated runs of each filter. Burn-in $ {n}_B $  = 5.

Figure 15. Updating of the random field $ \ln \left(D\left(t=50,x\right)\right) $ in three different cases of varying problem dimensionality. The solid lines show the mean and the shaded areas the 90% credible intervals inferred from 10 repeated runs of each filter. The black dashed line represented the posterior mean obtained via the KF, and the black solid lines the KF 90% credible intervals.

Increase of dimensionality does not seem to have any influence on the results obtained with the tPFGM filter. The illustrated efficacy of the tPFGM filter in estimating the time-invariant parameters in all considered cases of increasing dimensionality is related to the nature of the studied problem. The tPFGM filter relies entirely on the GMM approximation of the posterior distribution within its resampling process, in that it simply “accepts” all the $ {N}_{\mathrm{par}} $ GMM-proposed particles, unlike the tIBIS and SMC filters, which contain the degenerating acceptance-rejection step within the IMH-GM move step. Clearly, the worse the GMM fit, the worse the expected performance of the tPFGM filter. The particular case investigated here has a Gaussian reference posterior solution, hence the GMM fitted by EM proves effective in approximating the posterior with a relatively small number of particles, even when going up to $ d $ = 200 dimensions, thus leading to a good proposal distribution for sampling $ {N}_{\mathrm{par}} $ new particles in tPFGM. As reported in Table 4, the tPFGM filter is associated with a significantly lower computational cost than its MCMC-based counterparts.