Parallelized adaptive Bayesian updating with structural reliability methods for inference of large engineering models

Abstract

The reassessment of engineering structures, such as bridges, now increasingly involve the integration of models with real-world data. This integration aims to achieve accurate ‘as-is’ analysis within a digital twin framework. Bayesian model updating combines prior knowledge and data with models to enhance the modelling accuracy while consistently handling uncertainties. When updating large engineering models, numerical methods for Bayesian analysis present significant computational challenges due to the need for a substantial number of likelihood evaluations. The novelty of this contribution is to parallelize adaptive Bayesian Updating with Structural reliability methods combined with subset simulation (aBUS) to improve its computational efficiency. To demonstrate the efficiency and practical applicability of the proposed approach, we present a case study on the Maintalbrücke Gemünden, a large railway bridge. We leverage modal property data to update a linear-elastic dynamic structural model of the bridge. The parallelized aBUS approach significantly reduces computational time, making Bayesian updating of large engineering models feasible within reasonable timeframes. The improved efficiency allows for a wider implementation of Bayesian model updating in structural health monitoring and maintenance decision support systems.

Introduction

The design of engineering structures such as road and railway bridges is based on models. When re-assessing structures, models are still considered, but inspections and in-service data gathered using non-destructive testing (NDT), diagnostic load tests, or structural health monitoring (SHM) are also considered. Updating the structural models with these data provides ‘as-is’ structural analysis models that more closely represent the actual performance of an engineering structure. When performed continuously, it can be part of a ‘digital twin’ framework (Ye et al., 2022) that enables diagnostic, prognostic, and decision support tasks.

When translating research into practice, it is crucial to consider the efficiency of the updating process and the quality of the updated models (Ye et al., 2022). Bayesian methods are becoming increasingly popular (Ye et al., 2022) which is attributed to their robust handling of unavoidable uncertainties in measurement and modelling processes (Huang et al., 2019), along with their ability to incorporate and update prior knowledge with observed data for nuanced interpretations (Pai and Smith, 2022). In this contribution, we also leverage Bayesian methods for model updating. Bayesian model updating quantifies parameter uncertainty as relative plausibilities of different parameter values (Beck, 2010; Huang et al., 2019). Selecting the most suitable model class that best explains the measured data without overfitting (Beck, 2010; Beck and Yuen, 2004; Kontoroupi and Smyth, 2017) is crucial for ensuring the quality of the model updating process (Simon et al., 2024). Additionally, evaluating the relative plausibilities of competing model classes can further improve the quality of the model predictions (Beck, 2010; Beck and Yuen, 2004; Papadimitriou et al., 2001).

The requirements for Bayesian model updating methods for large engineering structures are as follows: First, they must be computationally efficient to handle large models and provide solutions within reasonable time frames. A relatively recent example of updating an engineering model of a large suspension bridge (Jang and Smyth, 2017) took 3 months to compute. Second, there should be the possibility of exploring high-dimensional parameter spaces because models of large structures may have many relevant parameters. Third, the model class selection process requires knowledge of the model evidence, also known as marginal likelihood (Beck, 2010; Huang et al., 2019).

The state-of-the-art methods for obtaining samples from the updated distribution of the model parameters are Markov Chain Monte Carlo (MCMC) methods (Lye et al., 2021). Three classes of MCMC methods provide estimates of the model evidence: (a) Transitional Markov Chain Monte Carlo (TMCMC) (Ching and Chen, 2007), (b) nested sampling (Koune et al., 2023; Skilling, 2004; Titscher et al., 2023; Xu et al., 2022), and (c) Bayesian updating with structural reliability methods (BUS) (Straub and Papaioannou, 2014), and its extension adaptive Bayesian updating with structural reliability methods using subset simulation (aBUS) (Betz et al., 2018b). The evaluation of the relative performance and accuracy of these methods is the subject of ongoing research (Smith et al., 2023). TMCMC can be parallelized (Goller et al., 2011), which greatly improves its efficiency. However, if the number of uncertain parameters is high, the potential for considerable bias arises, and the efficiency decreases (Betz et al., 2016). BUS with subset simulation promises to be applicable to high-dimensional problems. However, methods that fulfil all three requirements are missing.

We conclude that a method for Bayesian updating is needed which is efficient, works in high-dimensional parameter spaces, and provides an estimate of the evidence of a model class. A solution to closing this gap is to improve one of the three existing methods that provide an estimate of the model evidence. This could be done by addressing the performance issues of nested sampling or TMCMC when applied to problems with highly-dimensional parameter spaces. However, in our contribution we focus on improving the efficiency of aBUS through parallelization.

First, we introduce Bayesian model updating and aBUS. Then we show, how aBUS can be robustly parallelized. This is our primary contribution. Subsequently, we apply the parallelized version of aBUS to update a large finite element model of a railway bridge and evaluate the achieved speedup compared to the sequential version of aBUS. We complete the paper with a discussion and concluding remarks.

Bayesian model updating

Bayesian updating of (probabilistic) engineering models – like Bayesian inference of probabilistic models in general – can be performed on two levels: First, the probability distribution of the model parameters is updated with observed data, as briefly discussed in the following. Subsequently, the relative plausibilities of proposed candidate models are compared.

Parameter inference

To update the prior distribution p( θ ) of the uncertain parameters θ of an engineering model, the following formal procedure is applied.

As a first step, a data prediction model is defined which relates the i-th observation y _i of the system response with the system response predicted by the deterministic engineering model q ( x _i , θ ) based on the model input x _i and the model parameters θ . This is achieved by introducing an uncertain prediction error ϵ _i (Beck, 2010; Koune et al., 2023; Ntotsios et al., 2009; Simoen et al., 2013; Simon et al., 2024) which describes the discrepancy between the model-predicted and observed system response. The prediction error combines the uncertainty arising from both modelling and observation errors which includes measurement uncertainty.

${\mathit{y}}_i=\mathit{q}({\mathit{x}}_i,\mathit{\theta })+{\mathit{ϵ}}_i$

(1)

Next, a probabilistic model of the prediction error ϵ _i is defined. For the purpose of illustration, ϵ _i is here assumed to follow a multivariate normal distribution with zero mean and covariance matrix Σ _ϵ ( θ ), that is, ${\mathit{ϵ}}_i\sim \mathcal{N}(\mathbf{0},{\mathbf{\Sigma }}_{\mathbf{ϵ}}(\mathit{\theta }))$ .

Generically describing the covariance matrix Σ _ϵ in function of the model parameters θ , essentially means that the parameters defining the variance and correlation of the prediction errors – in addition to the parameters of the engineering model – included as uncertain parameters in the formulation. Consequently, they can also be updated with the observed data (Koune et al., 2023; Simoen et al., 2013). To account for bias in the prediction error, the mean of prediction error can be also modelled as uncertain parameter and updated (Simon et al., 2024).

In order to learn the probabilistic distribution of the model parameters θ based on the observed data it is necessary to formulate a likelihood function. This function probabilistically describes the likelihood of observing the data given a certain realization of the model parameters θ. This formulation is based on the probability model of the prediction error and the data prediction model. We formulate a likelihood function $L(\mathit{\theta }|{\mathcal{D}}_i)$ for a dataset ${\mathcal{D}}_i=\{{\mathit{x}}_i,{\mathit{y}}_i\}$ , which is proportional to the probability density $p({\mathcal{D}}_i|\mathit{\theta })$ of observing the data ${\mathcal{D}}_i$ given the parameters θ . For this contribution, we will assume $L(\mathit{\theta }|{\mathcal{D}}_i)=p({\mathcal{D}}_i|\mathit{\theta })=\mathcal{N}({\mathit{y}}_i-\mathit{q}({\mathit{x}}_i,\mathit{\theta });\mathbf{0},{\mathbf{\Sigma }}_{\mathbf{ϵ}}(\mathit{\theta }))$ , that is:

$L(\mathit{\theta }|{\mathcal{D}}_i)=\frac{1}{\sqrt{{(2\pi )}^{N_D}\cdot \mathrm{d}\mathrm{e}\mathrm{t}({\mathbf{\Sigma }}_{\mathbf{ϵ}}(\mathit{\theta }))}}\cdot \exp (-J(\mathit{\theta },{\mathcal{D}}_i))$

(2)with

$J(\mathit{\theta },{\mathcal{D}}_i)=\frac{1}{2}{\left({\widehat{\mathit{y}}}_i-\mathit{q}({\widehat{\mathit{x}}}_i,\mathit{\theta })\right)}^T{{\mathbf{\Sigma }}_{\mathbf{ϵ}}}^{-1}\left({\widehat{\mathit{y}}}_i-\mathit{q}({\widehat{\mathit{x}}}_i,\mathit{\theta })\right),$

(3)where $J(\mathit{\theta },{\mathcal{D}}_i)$ is the measure-of-fit function (Fang et al., 2022; Ntotsios et al., 2009; Straub and Papaioannou, 2014) and ( y _i − q ( x _i , θ )) is the realization of the prediction error corresponding to ${\mathcal{D}}_i$ , which is an N _D -dimensional column vector.

The prediction errors ϵ _i of multiple datasets ${\mathcal{D}}_i$ are commonly assumed to be statistically independent conditional on the parameters θ . Thus, the likelihood function of the entire data $\mathcal{D}$ is the product of the likelihood functions of the individual datasets:

$L(\mathit{\theta }|\mathcal{D})=\prod _{i}L(\mathit{\theta }|{\mathcal{D}}_i)$

(4)

Applying Bayes’ rule yields the posterior probability density $p\left(\mathit{\theta }|\mathcal{D}\right)$ of the parameters conditional on the data $\mathcal{D}$ based on the likelihood function $L(\mathit{\theta }|\mathcal{D})$ and the prior probability density p( θ ) of the parameters:

$p\left(\mathit{\theta }|\mathcal{D}\right)=\frac{L\left(\mathit{\theta }|\mathcal{D}\right)p\left(\mathit{\theta }\right)}{{\int }_{\mathbf{\Theta }}L\left(\mathit{\theta }|\mathcal{D}\right)p\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }}=c_E^{-1}L\left(\mathit{\theta }|\mathcal{D}\right)p\left(\mathit{\theta }\right)$

(5)

The normalizing constant c _E is also referred to as model evidence:

$c_E={\int }_{\mathbf{\Theta }}L\left(\mathit{\theta }|\mathcal{D}\right)p\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }$

(6)

Prior knowledge of the uncertain parameters θ is reflected in their prior probability density $p\left(\mathit{\theta }\right)$ . The likelihood function $L\left(\mathit{\theta }|\mathcal{D}\right)$ connects data $\mathcal{D}$ from observations with the model parameters θ . With this framework, we now have a means of integrating prior knowledge and assumptions about the system encoded in the deterministic model q ( x , θ ) and the prior probability model of its parameters with data. The result is updated knowledge of the parameters in terms of their posterior probability density.

Model class selection

There are many ways of modelling a physical system. Furthermore, there are several ways of modelling the prediction error (Koune et al., 2023; Simoen et al., 2013). Consequently, we need a way of deciding, which way of modelling is more appropriate given the data.

To this end, we formally define a stochastic model class M _i , as introduced in (Beck, 2010) which consists of: (1) a parametrized deterministic engineering model of the physical system, (2) a parametrized probabilistic model of the discrepancy between the model output and the observations, and (3) a joint prior probability density of the uncertain parameters of the deterministic engineering model and the probabilistic models, reflecting prior knowledge.

All the relations derived in the previous section depend on the modelling choices and assumptions encoded in a given model class M _i , and thus equation (5) can be written as follows to reflect this dependence:

$\begin{array}{cc}\hfill p\left(\mathit{\theta }|\mathcal{D},M_i\right)& =\frac{L\left(\mathit{\theta }|\mathcal{D},M_i\right)p\left(\mathit{\theta }|M_i\right)}{{\int }_{\mathbf{\Theta }}L\left(\mathit{\theta }|\mathcal{D},M_i\right)p\left(\mathit{\theta }|M_i\right)\mathrm{d}\mathit{\theta }}\hfill \\ \hfill & =c_{E|M_i}^{-1}L\left(\mathit{\theta }|\mathcal{D},M_i\right)p\left(\mathit{\theta }|M_i\right)\hfill \end{array}$

(7)

When the likelihood function is a proper probability distribution, that is, $L(\mathit{\theta }|\mathcal{D})=p(\mathcal{D}|\mathit{\theta })$ , the evidence of the model class is the probability of observing the data given the model class:

$c_{E|M_i}=p(\mathcal{D}|M_i)={\int }_{\mathbf{\Theta }}L\left(\mathit{\theta }|\mathcal{D},M_i\right)p\left(\mathit{\theta }|M_i\right)\mathrm{d}\mathit{\theta }$

(8)

When Bayesian model updating is performed within each model class M _i of a set $\mathcal{M}={\{M_i\}}_{i=1}^{N_{MC}}$ of N _MC competing model classes, we get the probability of observing the data for each model class $p(\mathcal{D}|M_i)$ . To each of those model classes a prior probability p(M _i ) can be assigned, so that ${\sum }_i^{N_{MC}}p(M_i)=1$ The choice of prior probabilities of the postulated model classes depends on the user and should reflect prior knowledge on the system (Beck, 2010; Beck and Yuen, 2004; Ntotsios et al., 2009). Without further information, a non-informative prior which assigns to each model class a prior probability of $N_{MC}^{-1}$ can be chosen. This way, no model class is preferred a priori.

With these conditions met, Bayes’ rule can be applied anew, yielding the posterior probability $p(M_i|\mathcal{D})$ of each model class conditional on the available data $\mathcal{D}$ :

$p(M_i|\mathcal{D})=\frac{p(D|M_i)p(M_i)}{{\sum }_i^{N_{MC}}p(D|M_i)p(M_i)}$

(9)

We now can select the model class among all postulated model classes that most appropriately explains the data. This procedure is also referred to as model class selection (Beck, 2010; Beck and Yuen, 2004). In physics-based structural health monitoring, this method can be leveraged both for selecting the most appropriate model for example (Simon et al., 2024) and for damage identification for example (Ntotsios et al., 2009; Simon et al., 2024).

Adaptive BUS combined with subset simulation

The challenge for Bayesian updating of engineering models is that there is rarely a closed-form solution for equation (5) available. Additionally, when the likelihood function contains a deterministic engineering model, it is not continuous but can only be evaluated at discrete values. To overcome this, Markov chain Monte Carlo (MCMC) sampling methods – such as the Metropolis-Hastings algorithm (Hastings, 1970; Metropolis et al., 1953) or modern derivatives (Betancourt, 2017) – are commonly employed to draw samples from the posterior distribution $p(\theta |\mathcal{D})$ .

The goal is to obtain samples representative of the posterior distribution, while at the same time minimizing the number of likelihood evaluations. An overview of various categories of currently available methods is given in (Lye et al., 2021).

For Bayesian updating of engineering models, there are two specific requirements:

(1) Large engineering models typically have many uncertain parameters, which influence the structure’s behaviour and are thus of interest in the model updating. The applied method has to be able to cope with such high-dimensional problems.

(2) Several modelling choices can be made, regarding the applied physical/mechanical theories, and the level of detail of the engineering model. To enable a comparison of competing model classes representing different modelling choices with model class selection, the model evidence has to be computed.

To meet these requirements, the adaptive variant of Bayesian Updating with Structural reliability methods combined with subset simulation (aBUS) (Betz et al., 2018b; Straub and Papaioannou, 2014) is employed in this contribution.

In the following, we describe the workings of aBUS in three steps: (1) the main idea of BUS is discussed, then (2) aBUS with subset sampling is presented, and subsequently (3) the adaptive Conditional Sampling algorithm as an MCMC sampler in aBUS is described.

The main idea of BUS

Bayesian Updating with Structural reliability methods (BUS) reframes Bayesian model updating as a structural reliability problem. BUS expands the parameter space Θ with a uniform random variable Π to an augmented space Ω = Θ × [0, 1], defining the failure event F by the corresponding domain as $F=\{\pi \le c\cdot L(\mathit{\theta }|\mathcal{D})\}$ . The scaling constant c ensures that $c\cdot L(\mathit{\theta }|\mathcal{D})\le 1$ for all θ ∈ Θ. Several methods exist for choosing this constant (Betz et al., 2018a; Straub and Papaioannou, 2014).

Straub and Papaioannou (Straub and Papaioannou, 2014) demonstrate that samples from the prior distribution p( θ ) of the parameters that lie within this failure domain follow the posterior distribution $p(\mathit{\theta }|\mathcal{D})$ , effectively expressing Bayesian inference with structural reliability principles.

In the augmented parameter space Ω with realizations of the random variables ω = [π, θ ], the limit state function $h(\mathit{\omega })=\pi -cL(\mathit{\theta }|\mathcal{D})$ describes the failure domain, where a negative value of the limit state function corresponds to failure. Structural reliability methods can be applied to compute the failure probability Pr(F) = Pr(h( ω ) ≤ 0). From this probability, the model evidence c _E can be computed as:

$c_E=\mathrm{P}\mathrm{r}(F)\cdot c^{-1}.$

(10)

BUS typically involves applying sampling-based reliability methods to generate samples in the failure domain and compute c _E . Implementing BUS with simple Monte Carlo Simulation is equivalent to rejection sampling. Subset simulation (Au and Beck, 2001) is preferred as it is efficient in high-dimensional problems and also directly produces samples within the failure domain.

aBUS algorithm

Subset simulation, proposed by Au and Beck (Au and Beck, 2001), is a sequential Monte Carlo method for solving the structural reliability problem and calculating the probability of failure. It is based on the insight, that the failure event F can be expressed as an intersection of N _L intermediate nested events, where:

$F=F_{N_L}\cap F_{N_L-1}\cap \dots F_1\cap F_0,$

(11)and F₀ is the certain event with Pr(F₀) = 1. In the context of BUS, the intermediate events are defined as:

$F_i=\{h(\mathit{\omega })\le b_i\},$

(12)where $b_0=\infty >b_1>b_2\dots >b_{N_L}$ .

The probability of failure can thus be written as a product of conditional probabilities:

$\mathrm{P}\mathrm{r}(F)=\prod _{i=1}^{N_L}\mathrm{P}\mathrm{r}(F_i|F_{i-1})$

(13)

During the algorithm, the threshold values b _i are chosen so that the conditional probability Pr(F _i |F_i−1) for all but the last intermediate event is fixed to p₀. The conditional probabilities of the intermediate events are thus much larger than Pr(F):

$\mathrm{P}\mathrm{r}(F_i|F_{i-1})=p_0,i=1\dots N_L-1$

(14)

The threshold of the last subset level is $b_{N_L}=0$ , since the failure event is defined as $F=F_{N_L}=\{h(\mathit{\omega })\le 0\}$ .

This idea of expressing the failure event as an intersection of a sequence of nested intermediate events enables a sequential sampling-based approximation of Pr(F). This approximation has three steps:

(1) Draw samples ω ₍₀₎ according to the joint prior probability density function p( ω ) which are conditional on the certain event F₀.

(2) On each subset level i corresponding to an intermediate event F _i , N _S samples ω _(i) are generated with MCMC sampling conditional on the previous intermediate event F_i−1. This conditional sampling is achieved by (a) taking samples from the previous level as MCMC seeds, and (b) simulating states of Markov chains whose stationary distribution is equal to the target distribution p( ω |F_i−1).

(3) On the last subset level N _L the number of samples ${\mathit{\omega }}_{(N_L)}$ conditional on $F_{N_L-1}$ , that fall in the failure domain F are counted as N _F . The number of samples N _F in the failure domain enable an approximation of the conditional probability $\mathrm{P}\mathrm{r}(F_{N_L}|F_{N_L-1})\approx N_F/N_S$ . Subsequently, the probability of failure can be estimated accordingly to equation (13).

From the probability Pr(F), the model evidence is estimated according to equation (10). In addition, samples of θ falling in the failure domain are samples from the posterior distribution $p(\mathit{\theta }|\mathcal{D})$ .

One of the most performance-crucial issues in the BUS approach is the choice of the scaling constant c.

To alleviate this inefficiency, the adaptive version of BUS, aBUS, adapts the scaling constant h _C = ln(c⁻¹) at every subset level i to the largest value of the likelihood function evaluated for the samples ω _(i). Thus, the need for choosing an initial value of c becomes obsolete. The change in c does not affect the distribution of the current samples, if the current threshold level is adjusted (Betz et al., 2018b). For improved efficiency and numerical stability, a logarithmic transformation of the limit state function (Betz et al., 2018b; DiazDelaO et al., 2017) is applied, which entails using the log-likelihood function $l(\mathit{\theta }|\mathcal{D})=\mathrm{l}\mathrm{n}(L(\mathit{\theta }|\mathcal{D}))$ . The limit state function becomes thus:

$g(\mathit{\omega })=\ln (\pi )+h_C-l(\mathit{\theta }|\mathcal{D})$

(15)

Algorithm 1 shows the basic aBUS approach.

Algorithm 1: aBUS.

Input:

•   p( θ ): Joint prior distribution of the parameters Θ

•   $l(\mathit{\theta }|\mathcal{D})$ : Log-likelihood function of the parameters conditional on the data $\mathcal{D}$

•   N _S : Number of samples per subset level

•   N _P : Number of desired samples of the posterior distribution; usually the same as N _S

•   p₀: Fixed conditional probability for intermediate events F _i ; Note, that the number of samples N _T = p₀ ⋅ N _S below the thresholds b _i of the subset levels must be an integer.

Output:

•  c _E : Estimate of the evidence

•   ${\mathit{\theta }}_{(N_L)}$ : Samples of the last subset level as an approximation of the posterior distribution $p(\mathit{\theta }|\mathcal{D})$

• (1) Initialization

• (ii)

  Formulate limit state function h( ω ).

• (iii)

  Draw N _S initial samples ω ₍₀₎ from the prior distribution p( ω ) and evaluate the likelihood function for each sample.

• (iv)

  Set the scaling constant h _C to the maximum of the log-likelihood function of the initial samples.

• (v)

  Set the initial threshold to b₀ = ∞ and the counter variable i = 0

(2) Subset iteration while b _i > 0:

• (i)

Increase counter by one: i = i + 1.

(ii)

Determination of the threshold b _i and seeds ${\mathit{\omega }}_{(i)}^{*}$ .

•  Set the threshold b _i so that p₀ ⋅ N _S of the samples fulfil h( ω ) ≤ b _i .

•  If b _i ≤ 0, set b _i = 0.

•  Count the number N _T of samples with the limit state function h( ω ) below the threshold b _i . Take these samples as seeds for subsequent MCMC sampling.

•  Set Pr(F _i |F_i−1) ≈ N _T /N _S , which for b _i > 0 is Pr(F _i |F_i−1) = p₀.

(iii)

Generate samples conditional on domain F_i−1 by:

•  Randomize the ordering of the seeds ${\mathit{\omega }}_{(i)}^{*}$ .

•  Take the N _T seeds ${\mathit{\omega }}_{(i)}^{*}$ and generate (N _S − N _T ) additional samples that fulfil h( ω ) ≤ b _i using an MCMC sampler. Note, that this requires the evaluation of the limit state function and thus also log-likelihood function for all candidate samples.

(iv)

Update the value of the scaling constant. This is the key step of the adaption in aBUS. For a detailed derivation, see (Betz et al., 2018b):

•  Set the (possibly) adapted scaling constant to $h_C^{*}=\max (h_C,l({\mathit{\omega }}_{(i)}|\mathcal{D}))$ .

•  Modify the threshold level to $b_i^{*}=b_i-h_C+h_C^{*}$ . Note that if the scaling constant is not adapted when $h_C^{*}=h_C$ , also the threshold level b _i doesn’t change.

•  Set $h_C=h_C^{*}$ and $b_i=b_i^{*}$ .

(v)

Decrease the dependence of the k = 1…N _S samples ω _(i). For each of the samples ω _(i,k) = ( ω _(i,k), π_(i,k)) of subset level i:

•  Draw ${\pi }_{(i,k)}^{*}$ from a uniform distribution $\mathcal{U}(1,\min (1,\exp (l({\mathit{\omega }}_{(i,k)}|\mathcal{D})-h_C+b_i)))$ .

•  Set ${\mathit{\omega }}_{(i,k)}=({\mathit{\omega }}_{(i,k)},{\pi }_{(i,k)}^{*})$ .

(3) Post-processing:

• (i)

Set N _L = i. The posterior samples are the samples ${\mathit{\theta }}_{(N_L)}$ of the last subset level.

(ii)

Estimate $\mathrm{P}\mathrm{r}(F)={\prod }_{i=1}^{N_L}\mathrm{P}\mathrm{r}(F_i|F_{i-1})$ .

(iii)

Evaluate the evidence as c _E = exp(h _C ) ⋅Pr(F).

When choosing the value of p₀, one has to consider a trade-off between efficiency, which is expressed as the total number of likelihood evaluations, and the approximation errors in the results of the algorithm. A choice of p₀ = 0.1 is “close to the optimum in many cases” (Betz et al., 2018b). Note that there are also approaches which adapt p₀ for each intermediate event (Blandfort, 2021).

In addition to selecting the conditional probability of the intermediate events p₀, the number of samples N _S per subset level must be chosen. N _S ≥1000 is recommended (Betz et al., 2018b).

Adaptive conditional sampling in standard normal space

To generate samples conditional on each intermediate event F_i−1 is to draw samples on each subset level F _i conditional on the previous level F_i−1. Specialized Markov chain Monte Carlo (MCMC) methods are employed which use the existing samples conditional on F_i−1 as seeds of the individual MCMC chains. This procedure ensures that the generated samples are distributed according to p( ω |F_i−1). In this contribution, we employ adaptive conditional sampling in standard normal space, which is simple and efficient (Papaioannou et al., 2015). The algorithm requires a suitable transformation T of samples ω _(k) from the (physical) domain Ω to samples u _(k) from the standard normal space U (Hohenbichler and Rackwitz, 1981; Liu and Kiureghian, 1986).

The main component of the algorithm is the transition from the current sample u _(k−1) to the next sample u _(k). This transition is performed as follows: The algorithm generates a candidate ${\mathit{u}}_{(k)}^{\prime }$ component-wise for each component j from a proposal distribution, which follows $\mathcal{N}(u_{(k-i)}^j\sqrt{1-s_q^{j2}},s_q^j)$ , where $s_q^j$ is the component-wise spread parameter. Candidates are accepted or rejected based on whether they fall in the intermediate failure domain F_i−1. Conditional sampling is shown in Algorithm 2.

Algorithm 2: Conditional Sampling.

Input:

• u _(k−1): Initial sample in standard normal space U .

• ω _(k−1): Initial sample in the physical domain Ω.

• $s_q^j$ : spread parameter of the proposal distributions for each parameter j.

• T: Transformation from the domain Ω to the standard normal space U .

• intermediate failure domain F_i−1, given by the limit state function h( ω ) and the threshold b_i−1.

Output:

•  u _(k): Next sample in standard normal space U .

•  ω _(k): Next sample in the physical domain Ω.

•  α_(k): indicator of acceptance of sample k, as boolean or integer variable.

• (1) Sampling of a candidate

(ii)

Transform candidate sample to physical domain: ${\mathit{\omega }}_{(k)}^{\prime }=T^{-1}(u_{(k)}^{\prime })$ .

(2) Acceptance/Rejection

• (i)

If ${\mathit{\omega }}_{(k)}^{\prime }$ is within F_i−1, that is h( ω ) ≤ b_i−1:

• Set ${\mathit{u}}_{(k)}={\mathit{u}}_{(k)}^{\prime }$ .

• Set ${\mathit{\omega }}_{(k)}={\mathit{\omega }}_{(k)}^{\prime }$ .

• Set α_(k) = 1

(ii)

Else:

• Set u _(k) = u _(k−1).

• Set ω _(k) = ω _(k−1).

• Set α_(k) = 0

The choice of the spread of the proposal distribution, which is expressed in terms of s _q , is crucial for the algorithm’s efficiency. When the spread is large, that is for small s _q , the candidate samples ${\mathit{u}}_{(k)}^{\prime }$ are likely to lie outside the transformed intermediate failure domain F _i , and thus rejected. This results in many correlated, that is, dependent, samples. When the spread is small, that is for large s _q , the samples are likely to be accepted but are also highly correlated. When samples are highly correlated, a greater number of total samples is needed to represent the target distribution appropriately. This requires more computationally costly likelihood evaluations.

An implicit measure of the dependence of the samples, and thus of the efficiency of the algorithm, is the acceptance rate α as the ratio of accepted candidate samples to the total number of candidate samples (Betz et al., 2018b). For subset simulation, α* = 0.44 is suggested as an optimal acceptance rate (Papaioannou et al., 2015). Betz et al. (Betz et al., 2018b) investigate acceptance rates for aBUS with case studies and suggest that values between 0.3 ≤ α ≤ 0.5 are efficient.

In the conditional sampling algorithm 2, the acceptance rate depends on the spread of the proposal distribution, which depends on s _q . It is impossible to select the optimal spread of the proposal distribution before evaluating the samples. Papaioannou et al. (Papaioannou et al., 2015) suggest adaptively adjusting the spread of the proposal distribution during the simulation to meet the optimal acceptance rate α*.

This is achieved as follows: A number i = 1, …, N _A of adaptions is performed while generating the MCMC chains. The adaption consists of step-wise adjusting the parameter s _q of the proposal distribution.

A starting value for the spread parameter $s_q^j$ for each component j of the proposal distribution is chosen. It can be either set to a fixed value, for example 0.8, or informed by the standard deviation of each component of the supplied seeds, see (Papaioannou et al., 2015). In each adaption i, the spread parameter $s_q^j$ is multiplied by a coefficient c _q according to

$s_q^j=\min (\exp (c_q)s_q^j,1.0)$

(16)

In each adaption step, N _C chains are generated. For each of these chains, the acceptance rate is calculated. The average acceptance rate ${\widehat{\alpha }}_i$ of these chains in adaption step i is calculated by averaging over the ratio of accepted samples to total samples of the individual chains. Based on the current adaption step i, the scaling parameter of the next adaption is:

$c_q=\frac{{\widehat{\alpha }}_i-{\alpha }^{*}}{\sqrt{i}}$

(17)

The number N _C of chains considered for adaption, should be selected such that there are 5 to 10 adaption steps, which can be calculated as N _A = p₀ ⋅ N _S /N _C (Papaioannou et al., 2015). Algorithm 3 implements the adaptive version of conditional sampling.

Algorithm 3: Adaptive Conditional Sampling (aCS).

Input:

• N _S : Number of samples to be drawn.

• N _C : Number of chains to consider for each adaptation.

• $s_{q,\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}^j$ : Initial component-wise spread parameter, optional.

• c_q,init: Initial scaling parameter, optional.

• u _(i,k): A number k = 1, …, N _T of seeds for the MCMC chains in standard normal space, with N _T = p₀ ⋅ N _S .

• α*: Target acceptance rate.

Output:

• u _(i+1): N _S samples for subset level i + 1 in standard normal space

• c_q,last: Scaling parameter of the last adaption

• (1) Initialization

• (ii)

  Calculate the number of adaptions N _A as the integer division of (N _S div N _C )

(2) Iterative Sampling for adaption steps i = 1, …, N _A :

• (i)

For each c ∈ 1, …N _C of the chains:

• Calculate seed index as k = (i − 1)N _C + c. Take the seed u _(i,k) and use algorithm 2 for conditional sampling to generate $N_D={p_0}^{-1}$ states of a Markov chain with the current spread parameters $s_q^j$ . Keep track of the acceptance indicator α_(k) of each generated sample k ∈ 1, …, N _D .

(ii)

For the states generated with the N _C chains, evaluate the average acceptance rate ${\widehat{\alpha }}_i={\sum }_{k=1}^{c\cdot N_D}{\alpha }_{(k)}/(c\cdot N_D)$ .

(iii)

Compute the scaling parameter c _q based on equation (17).

(iv)

Update the spread parameters $s_q^j$ according to equation (16).

(3) If there is a remaining number of chains (N _S mod N _C ) repeat step 2.i.

(4) Return the generated samples u _(i+1) and the last scaling parameter c_q,last.

Failure-robust parallelization of aBUS

In the application of the aBUS methodology for Bayesian updating, a substantial number of likelihood evaluations are necessitated. A key observation is that some of these evaluations are independent of one another, enabling the parallelized execution. Specifically, within the aBUS algorithm, the likelihood evaluations for the initial samples (as outlined in Step 1.iii) are independent. Similarly, in the aCS algorithm, the generation of Markov chains (referenced in Step 2.i) during the adaptation phases can also be conducted independently.

The analysis of engineering models typically involves the use of complex software systems. Running multiple instances of these systems in parallel, especially when interfaced with external control libraries, may lead to crashes due to resource over-utilization and communication-related issues between the software and the controlling scripts. We experienced crash rates α _c of up to 0.5%. This presents a challenge for the parallelization of the aBUS and aCS algorithms. To address this issue, our approach permits these crashes while keeping track of the samples of the failed analyses. This strategy allows for the subsequent re-execution of these analyses, effectively managing the instability without compromising the integrity of the overall analysis.

We integrate strategies of parallel execution and crash-handling for a failure-robust parallelization of the aBUS and aCS algorithms. First, Step 1.iii in the algorithm for aBUS is parallelized. The maximum number N _W of parallel processes is equal to the number of samples N _S . Figure 1 shows a flowchart of this algorithm, as described in Algorithm 4.

Algorithm 4: Robust parallelization for Step 1.iii in the algorithm for aBUS.

Input:

• N _S : Number of samples per subset level (this equals the number of initial samples).

• $N_R^{*}$ : A maximum number of retries. This should be set large enough so that the probability of successive crashes ${\alpha }_c^{N_R^{*}}$ is sufficiently low. Also, it should be set low enough to prevent infinite loops, in case of parameters that lead to corrupt analyses due to aspects inherent in the engineering model.

• p( ω ): Joint prior distribution of the augmented variable u in standard normal space.

• $l(\mathit{\theta }|\mathcal{D})$ : Log-likelihood function of the parameters conditional on the data $\mathcal{D}$ .

Output:

• ω ₍₀₎: Initial samples of the augmented variable u in standard normal space.

• $l({\mathit{\theta }}_{(0)}|\mathcal{D})$ : Evaluations of the log-likelihood function for the initial samples.

(1) Initialization

• (i)

Draw N _S initial samples ω ₍₀₎ from the joint prior distribution and copy them into a list L _S . This list is used to track the samples which are still to be evaluated.

(ii) Set the counter variable N _R of retries to 0.

(2) Iterate, while the list L _S of samples to be evaluated is not empty and $N_R\le N_R^{*}$ :

• (i)

Parallelly evaluate the log-likelihood function $l(\mathit{\theta }|\mathcal{D})$ for the samples of the list.

(ii)

Remove the samples from the list that were successfully evaluated.

(iii)

Increase N _R by 1.

(iv)

If $N_R>N_R^{*}$ , raise an Error.

(3) Return the initial samples and the corresponding evaluations of the log-likelihood function, or raise an exception if $N_R>N_R^{*}$ . This exception has to be handled by the user since it indicates failures related to the mathematical, physical, or computational aspects of the engineering model.

Figure 1.

Flowchart of algorithm 4: Robust parallelization for Step 1.iii in the algorithm for aBUS.

Secondly, the generation of Markov chains in the algorithm for aCS (Step 2.i) is parallelized between the adaptions of $s_q^j$ . Here, the upper limit of parallel executions is the number of chains N _C to consider for each adaptation. Figure 2 shows a flowchart of this algorithm (see Algorithm 5).

Algorithm 5: Robust parallelization for Step 2.i in the algorithm for adaptive Conditional Sampling.

Input:

• N _S : Number of samples to draw

• N _C : Number of chains to consider for each adaptation

• $s_{q,\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}^j$ : Initial component-wise spread parameter, optional.

• c_q,init: Initial scaling parameter, optional.

• u _(i,k): A number k = 1, …, N _T of seeds for the MCMC chains in standard normal space, with N _T = p₀ ⋅ N _S .

• a*: Target acceptance rate

Output:

• u _(i+1): N _S samples for subset level i + 1 in standard normal space

• c_q,last: Scaling parameter of the last adaption

(1) Initialization

• (i)

For each dimension j of the parameter space, choose an initial spread parameter $s_q^j$ either based on the standard deviation of the seeds u _(i,k) or as a fixed value $s_{q,\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}^j$ . If an initial scaling parameter c_q,init is provided, adapt the spread parameter according to equation (16).

(ii)

Store the seeds from the previous subset level u _(i,k) in a list L _S .

(2) Batch-wise parallel sampling while the list of seeds is not empty:

• (i)

Take N _C seeds from the list L _S and use algorithm 2 for conditional sampling in parallel with N _W worker processes to generate states of a Markov chain. Keep track of which seeds the generation of a Markov chain failed and note the acceptance indicator α_(k) of each generated sample.

(ii)

Remove the seeds for which the conditional sampling suceeded from the list L _S of seeds.

(iii)

Evaluate the average acceptance rate ${\widehat{\alpha }}_i$ of the N _C chains.

(iv)

Compute the scaling parameter c _q based on equation (17).

(v)

Update the spread parameters $s_q^j$ according to equation (16).

(3) Return the generated samples u _(i+1) and the last scaling parameter c_q,last.

Figure 2.

Flowchart for algorithm 5: Robust parallelization for Step 2.i in the algorithm for adaptive Conditional Sampling.

Case study: Vibration-based system identification of the Maintalbrücke Gemünden

The purpose of this case study is to demonstrate the feasibility of updating large engineering models using sampling-based Bayesian updating. To this end, a linear-elastic dynamic structural model of a large railway bridge is parametrized by several uncertain parameters. Their probabilistic distribution is learned based on modal property data.

The bridge

The Hanover-Würzburg high-speed railway line built in the 1970/80s crosses the Main valley near Gemünden am Main, approximately 70 km east of Frankfurt. At the request of the infrastructure owner Deutsche Bahn, the engineering consultancy Leonhardt, Andrä und Partner revised the initial design for the viaduct (Leonhardt, 1986) and proposed a rigid frame bridge with v-shaped piers for the main bridge instead of a tapered three-span continuous girder bridge. In the following, the approach viaducts are not considered and the term “Maintalbrücke Gemünden” is used synonymously for the main bridge across the river, which is shown in Figure 3.

Figure 3.

Picture of the Maintalbrücke Gemünden spanning over the river Main.

The superstructure acts primarily as a continuous girder with variable spans. The v-shaped piers reduce the main span from 135 m to 105 m. However, the system is semi-integral because the v-shaped piers are joined monolithically to the girder. The cross-section was designed as a box girder with variable heights ranging from 4.5 m to 6 m. A special feature of the structure are the concrete hinges at the bottom of the v-shaped piers, which are among the most heavily loaded of their kind (Schacht and Marx, 2010). Along the box girder, transverse bulkheads are located above the v-shaped piers and at both ends of the box girder.

The data

A permanent monitoring system has been installed on the bridge for research purposes. The system records data during most train passages and twice per day for longer periods. Summary statistics of the measured quantities are aggregated every 30 s. The measured quantities are strains, temperatures, inclinations, accelerations, and vibration velocities. In this contribution, we utilize vibration velocity data recorded with geophones, whose layout is shown in Figure 4. There are 13 vertical geophones distributed along the longitudinal axis of the bridge, which are placed on the floor of the box girder at the western side. Additional geophones are placed on the floor of the box girder at the eastern side to capture torsional mode shapes.

Figure 4.

Elevation of the bridge with the sensor layout of the geophones.

In this case study, a 17-min record of geophone data from October 2020 was used to identify natural frequencies and mode shapes of the bridge. To this end, covariance-driven stochastic subspace identification (SSI-COV) (Peeters and De Roeck, 2001) combined with clustering methods for automated post-processing (Reynders et al., 2012) was applied. The identified modal properties summarized in Table 1 were subsequently employed as a basis for updating the structural model of the Maintalbrücke Gemünden.

Table 1.

Identified natural frequencies and corresponding mode shapes.

	Type	Frequency	Mode shape
1	Vertical bending (1st mode main span)	1.97 Hz
2	Vertical bending (2nd mode main span; 1st mode side spans)	3.28 Hz
3	Vertical bending (2nd mode main span; 1st mode side spans)	3.38 Hz
4	Vertical bending (2nd mode main span)	4.85 Hz
5	Torsion (main span)	5.65 Hz
6	Vertical bending (3rd mode main span)	6.56 Hz
7	Torsion (side span)	6.66 Hz
8	Vertical bending (mixed)	7.68 Hz
9	Vertical bending (2nd mode side spans)	7.97 Hz
10	Vertical bending (3rd mode main span; 2nd mode side spans)	8.50 Hz
11	Vertical bending (mixed)	9.11 Hz
12	Vertical bending (mixed)	9.70 Hz

The model

The bridge was represented by a linear-elastic finite-element (FE) model, which primarily consisted of shell elements. The model was defined based on as-built drawings. Based on these drawings, the box girder was partitioned into 86 sections with different geometries. Preliminary investigations showed that the main girder alone could not explain the high torsional stiffness. Therefore, the bridge caps were modelled as additional beam elements and coupled to the outer edges of the main girders. Additional beam elements were connected to the top of the main girder with springs to represent the train tracks. In this way, a moving static load could be applied to the model. The results of this analysis were compared with deformation data recorded during load tests to validate the model (Simon et al., 2022).

Solid elements were used to model the massive blocks of concrete at the bottom of the v-shaped piers which transfer the loads to the concrete hinges. The supports on the piers at both ends of the bridge and the concrete hinges were modelled with spring elements that have a high stiffness for translational degrees of freedom (DOF) but a low stiffness for rotational DOF.

Figure 5 shows the resulting model with 87,180 DOF and highlights the described modelling details. The parameters of the FE model whose probabilistic distribution is updated (see next section) are also indicated in the figure. The rationale for selecting these parameters is as follows: The five sections of the superstructure with different elastic moduli correspond to the main sections produced during the construction of the bridge. To accommodate the aforementioned effect of the bridge caps, their density, elastic modulus, and stiffness of the shear connection were included as uncertain model parameters. The track mass was also considered as an uncertain model parameter. Thus, the masses of the various non-loadbearing members are approximately captured by both the track mass and bridge cap density. The mode shapes in Table 1 indicate that the bearings exhibit some flexibility. Therefore, they were modelled as spring elements that could be parametrized.

Figure 5.

Overview of the finite element model of the bridge, including parameters.

Three probabilistic model classes were formulated: one (referred to as M₁) which includes the stiffness K_V,T,Y of the spring elements at the bottom of the v-shaped piers as a deterministic parameter with a value of 10¹² N/mm, a second one (referred to as M₂) which includes this stiffness as an uncertain parameter, and a third one (referred to as M₃) which also includes the stiffness K_V,T,Y but models the standard deviation of the prediction error as a deterministic parameter with a value of 0.15.

These three model classes are defined to demonstrate the process of model class selection. In essence, this process can be used to determine the most plausible way of modelling the system which includes the selection of the parameters considered in the model updating. The issue of choosing the relevant model classes from the space of possible model classes remains an open research area. For a current discussion, the reader is referred to (Simon et al., 2024).

Model updating

As a basis for updating the structural model, the prior probabilistic model of the parameters is defined and the likelihood function containing the structural model is formulated. These two (the prior probabilistic model and likelihood function) form a probabilistic model class which can be updated with Bayesian updating. The results are shown subsequently.

Formulation of the likelihood function

The likelihood function is formulated based on the prediction error ϵ as introduced previously. In model updating with modal property data, the prediction error represents the discrepancies between observed and predicted natural frequencies and mode shapes (Behmanesh et al., 2015; Behmanesh and Moaveni, 2016; Hou et al., 2020; Huang et al., 2019; Ntotsios et al., 2009; Simoen et al., 2015). The realizations of the prediction errors are calculated for each observed mode, that can be assigned to a predicted mode.

In this contribution, we define a mode as a pair ( ϕ , λ) of a mode shape vector $\mathit{\varphi }\in {\mathbb{R}}^{N_{\mathit{DOF}}}$ which has N_DOF dimensions and an eigenvalue λ. An eigenvalue λ is obtained from the corresponding natural frequency f as:

$\lambda ={(2\pi f)}^2.$

(18)

There can be more model-predicted modes than observed modes. Also, the sequence of natural frequencies ${\widehat{f}}_i$ of the observed modes may differ from the sequence of natural frequencies f _j of the model predictions. This problem occurs more frequently when closely spaced modes are present. Matching predicted with observed modes is thus not trivial. For this task, methods of mode paring need to be employed based on the discrepancy between predicted and observed modes. In this contribution, an indicator δ _ij quantifying the discrepancy between each of the i observed modes and each of the j predicted modes is calculated based on the ratio of the eigenvalues λ and the discrepancy between two mode shapes via the modal assurance criterion (MAC) (Reynders, 2012). The MAC is the squared cosine between two mode shape vectors:

$\mathrm{M}\mathrm{A}\mathrm{C}({\mathit{\varphi }}_i,{\mathit{\varphi }}_j)=\frac{{\left({\mathit{\varphi }}_i^T{\mathit{\varphi }}_j\right)}^2}{\Vert {\mathit{\varphi }}_i{\Vert }_2^2\Vert {\mathit{\varphi }}_j{\Vert }_2^2},$

(19)where ‖ ⋅ ‖₂ signifies the Euclidean norm of a vector. For comparison, mode shape vectors must align, selecting only the model’s degrees of freedom (DOFs) that correspond to the sensor locations and directions in the correct order.

The discrepancy δ _ij between an observed mode $({\widehat{\mathit{\varphi }}}_i,{\widehat{\lambda }}_i)$ and a predicted mode ( ϕ _j , λ _j ) is composed of the discrepancy δ_ij,λ of the eigenvalues and the discrepancy δ_ij,m of the mode shapes as follows (Simoen et al., 2015):

${\delta }_{ij}=\stackrel{{\delta }_{ij,m}}{\overbrace{1-\mathrm{M}\mathrm{A}\mathrm{C}({\mathit{\varphi }}_j,{\widehat{\mathit{\varphi }}}_i)}}+\stackrel{{\delta }_{ij,\lambda }}{\overbrace{\left|1-\frac{{\lambda }_j}{{\widehat{\lambda }}_i}\right|}}.$

(20)

Mode pairing algorithms must seek to match the corresponding predicted and observed modes even when there is a difference in the eigenvalues and mode shapes. At the same time, modes that bear some similarity but do not correspond to each other, such as closely spaced modes, must not be matched.

Our proposed algorithm (as shown in Algorithm 6) addresses these requirements by: (1) sequentially matching the best-fitting pair and removing them from the pool of candidates, and (2) defining a threshold δ* on the discrepancy δ _ij above which no modes are paired. The first part supports correct matching even when closely spaced modes are present. The definition of the threshold δ* prevents the matching of unrelated modes.

Algorithm 6: Mode pairing.

Input:

• δ*: Threshold value for maximal discrepancy

• ${\{({\widehat{\mathit{\varphi }}}_i,{\widehat{\lambda }}_i)\}}_{i=1}^{N_I}$ : a set of N _I observed modes

• ${\{({\mathit{\varphi }}_j,{\lambda }_j)\}}_{j=1}^{N_J}$ : a set of N _J predicted modes

Output:

• L _P : list of pairs of observed and predicted modes.

(1) Initialize an empty list L _P = [] for modal pairs.

(2) Calculate the discrepancy matrix Δ, where each entry corresponds to the discrepancy δ _ij between observed mode i and predicted mode j according to equation (20). This matrix has the dimension of N _I × N _J .

(3) While the matrix Δ has dimensions $>0$ and there are entries $\le {\delta }^{*}$ do:

• Find the entry with the lowest discrepancy δ _ij . Record the indices i* and j* of the corresponding observed and predicted mode.

• Append the pair (i*, j*) to the list L _P .

• Delete row i* and column j* from the matrix Δ

(4) Return the list L _P of matched pairs of modes.

With this procedure, we can now align observed and predicted modes. This is the prerequisite for constructing prediction errors ϵ (as introduced in equation (1)) on a per-mode basis.

The usage of the mode shape discrepancy δ_ij,m assumes that even if two physical modes are close in frequency, they can be separated by their mode shapes. This implies that the mode shapes are orthogonal. Provided that the spatial resolution of the measured degrees of freedom (DOFs) is sufficient, bridges and buildings typically satisfy this condition. However, in quasi-axisymmetric structures, such as towers, special consideration is required for the modal identification. This is discussed in more detail in references (Au et al., 2021; Jonscher et al., 2023).

The available data $\mathcal{D}$ consists of N _M observed modes $({\widehat{\mathit{\varphi }}}_m,{\widehat{\lambda }}_m)$ that have to be matched with predicted modes.

The predicted modes result from an eigenvalue analysis, which is performed based on the deterministic linear-elastic structural model, which is denoted as q ( x , θ ). The eigenvalue analysis is independent of the inputs x and only depends on the model parameters θ . The dependency of the model-predicted modes on the parameters θ which govern the eigenvalue analysis is expressed through formulating the model-predicted modes as a function of θ in the following, that is ( ϕ _m ( θ ), λ _m ( θ )). The likelihood function is now formulated based on the observed and model-predicted modes.

There are different ways to model the prediction error of eigenvalues (or, likewise, frequencies), for example, additive (Behmanesh et al., 2015; Behmanesh and Moaveni, 2016; Huang et al., 2019; Simoen et al., 2015) or multiplicative (Koune et al., 2023; Ntotsios et al., 2009) formulations. The eigenvalues have different magnitudes. When the prediction error is modelled as multiplicative, the order of magnitude of the prediction error is independent of the order of magnitude of the eigenvalue. This formulation is thus deemed more suitable. It can be expressed as the ratio of predicted and observed eigenvalues. Taking the logarithm of this ratio yields the same order of magnitude of the prediction error, regardless of the order of magnitude of the numerator and denominator in the ratio (Simon et al., 2020). In this contribution, we define the prediction error for one eigenvalue of a matched mode with index m as:

${ϵ}_{\lambda ,m}=\left(\ln {\widehat{\lambda }}_m-\ln {\lambda }_m\left(\mathit{\theta }\right)\right)\sim \mathcal{N}\left(0,{\sigma }_{\lambda ,m}^2\right),$

(21)where we assume that the prediction error follows a normal distribution with a zero mean and standard deviation σ_λ,m. The prediction error depends on the parameters θ and the observed eigenvalue ${\widehat{\lambda }}_m$ . Based on this formulation of the prediction error, we can define the conditional probability density of the observed eigenvalue ${\widehat{\lambda }}_m$ given the realization of the parameters θ as (see equation (2) for reference):

$p\left({\widehat{\lambda }}_m|\mathit{\theta }\right)=\frac{1}{{\sigma }_{\lambda ,m}\sqrt{2\pi }}\cdot \exp \left(-\frac{\ln {\widehat{\lambda }}_m-\ln {\lambda }_m\left(\mathit{\theta }\right)}{{2\sigma }_{\lambda ,m}^2}\right).$

(22)

Due to the modelling choices of the current contribution, the prediction errors of different modes are statistically independent conditional on a realization of σ _λ , so that for all modes m: σ_λ,m = σ _λ . The standard deviation σ _λ of the prediction errors of the eigenvalues is modelled as uncertain parameter and are thus included in the parameter vector θ . The joint probability density of observing all the m ∈ 1, …, N _M matched eigenvalues given the realization of the parameters θ corresponds thus to the product of the marginal conditional probability densitity functions of the observed eigenvalues (Behmanesh et al., 2015; Goller et al., 2011; Hou et al., 2020; Huang et al., 2019; Jang and Smyth, 2017; Lye et al., 2021):

$\begin{array}{cc}\hfill p\left(\widehat{\lambda }|\mathit{\theta }\right)& =\prod _{m=1}^{N_M}p\left({\widehat{\lambda }}_m|\mathit{\theta }\right)\hfill \\ \hfill & =\prod _{m=1}^{N_M}\frac{1}{{\sigma }_{\lambda }\sqrt{2\pi }}\cdot \exp \left(-\frac{\ln {\widehat{\lambda }}_m-\ln {\lambda }_m\left(\mathit{\theta }\right)}{2{\sigma }_{\lambda }^2}\right)\hfill \\ \hfill & =\frac{1}{{\left({\sigma }_{\lambda }\sqrt{2\pi }\right)}^{N_M}}\hfill \\ \hfill & \cdot \prod _{m=1}^{N_M}\exp \left(-\frac{\ln {\widehat{\lambda }}_m-\ln {\lambda }_m\left(\mathit{\theta }\right)}{{2\sigma }_{\lambda }^2}\right).\hfill \end{array}$

(23)

The likelihood function $L\left(\mathit{\theta }|\widehat{\lambda }\right)$ has the same functional form as the probability density $p\left(\widehat{\lambda }|\mathit{\theta }\right)$ of observing $\widehat{\lambda }$ given θ . The difference is that $L\left(\mathit{\theta }|\widehat{\lambda }\right)$ is a function of θ and the value of $\widehat{\lambda }$ is fixed. If $L\left(\mathit{\theta }|\widehat{\lambda }\right)=p\left(\widehat{\lambda }|\mathit{\theta }\right)$ , the log-likelihood function $l\left(\mathit{\theta }|\widehat{\lambda }\right)=\ln L\left(\mathit{\theta }|\widehat{\lambda }\right)$ becomes:

$\begin{array}{cc}l\left(\mathit{\theta }|\widehat{\lambda }\right)& =\ln \left({\left({\sigma }_{\lambda }\sqrt{2\pi }\right)}^{-N_M}\right) +\sum _{m=1}^{N_M}-\frac{\ln {\widehat{\lambda }}_m-\ln {\lambda }_m\left(\mathit{\theta }\right)}{2{\sigma }_{\lambda }^2}\\ & =-N_M\ln \left({\sigma }_{\lambda }\sqrt{2\pi }\right)\\ & -\frac{1}{2{\sigma }_{\lambda }^2}\sum _{m=1}^{N_M}\ln {\widehat{\lambda }}_m-\ln {\lambda }_m\left(\mathit{\theta }\right).\end{array}$

(24)

To formulate a likelihood function based on the prediction errors of the mode shapes, we adopt a vector-valued additive prediction error model as introduced in (Beck et al., 2001):

${\mathit{ϵ}}_{\varphi ,m}={\widehat{\mathit{\varphi }}}_m-{\alpha }_m{\mathit{\varphi }}_m\left(\mathit{\theta }\right),$

(25)where α _m is a scaling factor. All components of ϵ _ϕ,m are assumed to be conditionally independent given a realization of ${\sigma }_{\varphi ,m}^2\Vert {\widehat{\mathit{\varphi }}}_m{\Vert }^2$ and follow a normal distribution $\mathcal{N}\left(0,{\sigma }_{\varphi ,m}^2\Vert {\widehat{\mathit{\varphi }}}_m{\Vert }^2\right)$ with zero mean and standard deviation ${\sigma }_{\varphi ,m}\Vert {\widehat{\mathit{\varphi }}}_m\Vert$ . This assumption leads to a multivariate normal distribution of the vector-valued prediction error (Beck et al., 2001). Additionally, we normalize matched observed and predicted mode shapes to unit length, that is $\Vert {\widehat{\mathit{\varphi }}}_m{\Vert }^2=\Vert {\mathit{\varphi }}_m\left(\mathit{\theta }\right){\Vert }^2=1$ . This normalization implies, that the standard deviation of the distribution of the prediction errors becomes σ_ϕ,m. Based on this formulation of the prediction error, we can define the conditional probability density of the observed mode shape ${\widehat{\mathit{\varphi }}}_m$ given the realization of the parameters θ as:

$\begin{array}{cc}\hfill p\left({\widehat{\mathit{\varphi }}}_m|\mathit{\theta }\right)& =\frac{1}{{\left({\sigma }_{\varphi ,m}\sqrt{2\pi }\right)}^{N_{\text{DOF}}}}\hfill \\ \hfill & \cdot \exp \left(-\frac{1-{\left({\widehat{\mathit{\varphi }}}_m^T{\mathit{\varphi }}_m\left(\mathit{\theta }\right)\right)}^2}{2{\sigma }_{\varphi ,m}^2}\right)\hfill \end{array}$

(26)(see Appendix A for a detailed derivation of this expression). Due to the modelling choices of the current contribution, the prediction errors of different mode shapes are statistically independent conditional on a realization of σ _ϕ , so that for all modes m: σ_ϕ,m = σ _ϕ . The standard deviation σ _ϕ of the prediction errors of the mode shapes is modelled as uncertain parameter and are thus included in the parameter vector θ . The joint probability density of observing all the m ∈ 1, …, N _M matched mode shapes is thus:

$\begin{array}{cc}\hfill p(\widehat{\mathit{\varphi }}|\mathit{\theta })& =\prod _{m=1}^{N_M}({\widehat{\mathit{\varphi }}}_m|\mathit{\theta })\hfill \\ \hfill & =\prod _{m=1}^{N_M}{\left({\sigma }_{\varphi }\sqrt{2\pi }\right)}^{{-N}_{\text{DOF}}}\hfill \\ \hfill & \cdot \prod _{m=1}^{N_M}\exp \left(-\frac{1-{\left({\widehat{\mathit{\varphi }}}_m^T{\mathit{\varphi }}_m\left(\mathit{\theta }\right)\right)}^2}{2{\sigma }_{\varphi }^2}\right)\hfill \\ \hfill & ={\left({\sigma }_{\varphi }\sqrt{2\pi }\right)}^{{-N}_{\text{DOF}}\cdot N_M}\hfill \\ \hfill & \cdot \exp \left(\sum _{m=1}^{N_M}-\frac{1-{\left({\widehat{\mathit{\varphi }}}_m^T{\mathit{\varphi }}_m\left(\mathit{\theta }\right)\right)}^2}{2{\sigma }_{\varphi }^2}\right).\hfill \end{array}$

(27)

When $p\left(\widehat{\mathit{\varphi }}|\mathit{\theta }\right)$ is viewed as a function of θ with $\widehat{\mathit{\varphi }}$ fixed, it becomes the likelihood function: $L\left(\mathit{\theta }|\widehat{\mathit{\varphi }}\right)=p\left(\widehat{\mathit{\varphi }}|\mathit{\theta }\right)$ , and the log-likelihood function $l\left(\mathit{\theta }|\widehat{\mathit{\varphi }}\right)=\ln L\left(\mathit{\theta }|\widehat{\mathit{\varphi }}\right)$ becomes:

$\begin{array}{cc}\hfill l(\mathit{\theta }|\widehat{\mathit{\varphi }})& =\ln \left({\left({\sigma }_{\varphi }\sqrt{2\pi }\right)}^{{-N}_{\text{DOF}}\cdot N_M}\right)\hfill \\ \hfill & +\sum _{m=1}^{N_M}\left(-\frac{1-{\left({\widehat{\mathit{\varphi }}}_m^T{\mathit{\varphi }}_m\left(\mathit{\theta }\right)\right)}^2}{2{\sigma }_{\varphi }^2}\right)\hfill \\ \hfill & ={-N}_{\text{DOF}}{N}_M\cdot \ln \left({\sigma }_{\varphi }\sqrt{2\pi }\right)\hfill \\ \hfill & -\sum _{m=1}^{N_M}\left(\frac{1-{\left({\widehat{\mathit{\varphi }}}_m^T{\mathit{\varphi }}_m\left(\mathit{\theta }\right)\right)}^2}{2{\sigma }_{\varphi }^2}\right).\hfill \end{array}$

(28)

The joint probability density of the observed data $\mathcal{D}=\{\widehat{\lambda },\widehat{\mathit{\varphi }}\}$ – conditional on the parameters θ – is equal to the product of the conditional probability density functions of the eigenvalues and mode shapes if the corresponding prediction errors are modelled to be statistically independent:

$p\left(\mathcal{D}|\mathit{\theta }\right)=p\left(\widehat{\mathit{\varphi }}|\mathit{\theta }\right)p\left(\widehat{\lambda }|\mathit{\theta }\right).$

(29)

The log-likelihood function $l\left(\mathit{\theta }|\mathcal{D}\right)$ of the parameters θ conditional on the observed data $\mathcal{D}$ is thus the sum of the log-likelihood functions derived above:

$l\left(\mathit{\theta }|\mathcal{D}\right)=l\left(\mathit{\theta }|\widehat{\mathit{\varphi }}\right)+l\left(\mathit{\theta }|\widehat{\lambda }\right).$

(30)

By constructing the likelihood function in this way, the statistical dependence among the observations of the identified modal eigenvalues $\widehat{\lambda }$ and mode shapes $\widehat{\mathit{\varphi }}$ because of the physical parameters and the standard deviations σ _ϕ and σ _λ , which are the parameters of the probabilistic model of the prediction errors, is explicitly captured.

Prior probabilistic model of the parameters

The parameters of each stochastic model class are modelled as random variables. The prior marginal distributions of the different model parameters are derived from the literature (Fédération internationale du béton, 2013; Joint Committee on Structural Safety (JCSS), 2001), design documents, and engineering judgment.

For example, considerations of upper and lower bounds were carried out for the parameters related to the stiffness of the bearings or shear connection, and these were subsequently modelled as uniformly distributed within these bounds. Since the parameters of the prediction error influence the results (Simoen et al., 2013), it is good practice to include them as uncertain parameters in the model class (Koune et al., 2023; Simon et al., 2024). They are listed in Table 2, with the relevant parameters mean μ, coefficient of variation CV, the lower limit a, and the upper limit b, where applicable. The first five parameters, which describe the elastic modulus of the box girder and v-shaped piers, are modelled with a pair-wise coefficient of correlation of ρ = 0.8.

Table 2.

Prior marginal distributions of parameters defined in terms of distribution type and parameters.

Parameter	Distribution	Parameters
Elastic modulus, Northern side span, EN	Log-normal	μ = 35000 MN/m2, CV = 0.15
Elastic modulus, v-shaped pier axis 7, EV7	Log-normal	μ = 35000 MN/m2, CV = 0.15
Elastic modulus, main span, EM	Log-normal	μ = 35000 MN/m2, CV = 0.15
Elastic modulus, v-shaped pier axis 8, EV8	Log-normal	μ = 35000 MN/m2, CV = 0.15
Elastic modulus, Southern side span, ES	Log-normal	μ = 35000 MN/m2, CV = 0.15
Elastic modulus, bridge caps, EC	Log-normal	μ = 35000 MN/m2, CV = 0.15
Shear connection between bridge caps and bridge, GC	Uniform	a = 103 N/mm, b = 109 N/mm
Track mass, μ	Log-normal	μ = 0.75 t/m2, CV = 0.10
Caps density, ρC	Log-normal	μ = 2.5 t/m3, CV = 0.5
Rotational bearing stiffness z, KB,R,Z	Uniform	a = 101Nmm/rad, b = 109Nmm/rad
Translational bearing stiffness y, KB,T,Y	Uniform	a = 104 N/mm, b = 107 N/mm
Translational bearing stiffness y, KV,T,Y a	Uniform	a = 106 N/mm, b = 1012 N/mm
Stdev. of pred. error eigenvalues, σ λ b	Log-normal	μ = 0.15 Hz, CV = 0.5
Stdev. of pred. error mode shapes, σ ϕ c	Log-normal	μ = 0.15, CV = 0.5

^aOnly included in model class M₂. In model classes M₁ and M₃ this parameter is deterministic with a value of 10¹² N/mm.

^bOnly included in model classes M₁ and M₂. In model class M₃ this parameter is deterministic with a value of 0.15.

^cOnly included in model classes M₁ and M₂. In model class M₃ this parameter is deterministic with a value of 0.15.

Alternatively, global sensitivity analyses can be performed to identify the relevant parameters to be considered in the model updating (Brownjohn et al., 2001; Jang and Smyth, 2017). This approach is, however, outside the scope of this case study.

Results

Model updating with parallelized aBUS was performed based on the three described model classes M₁, M₂, and M₃ and the identified modal property data. The updating was performed on an AMD Ryzen 7 PRO 5750G CPU with eight physical cores. A sample size of N _S = 4000 and a conditional probability for subset simulation of p₀ = 0.05 were used as settings for aBUS and a total of N _W = 10 worker processes were utilized. The computation times were 26 h, 18 h, and 23 h. Model updating provided the posterior samples of the parameters and an estimation for the evidence. The evidence is used for model class selection. Based on the posterior samples, we can update the probabilistic predictions of the model output.

Model class selection

For each postulated model class M _i , using i ∈ {1, 2, 3}, the evidence $c_{E|M_i}$ (as defined in Equation (8)) is estimated using the BUS method.

Because there is no further information to inform the prior probabilities of these model classes, we assume that they are equal to p(M _i ) = 1/3. This enables us to calculate the posterior probability $p(M_i|\mathcal{D})$ of each model class according to equation (9). Table 3 shows these results, which clearly indicate that model class M₁ (with a deterministic value for the stiffness of the springs below the v-shaped piers) explains the data better than model classes M₂ and M₃.

Table 3.

Logarithmic evidence and posterior probability used for model class selection.

Model class	$\ln (c_{E|M_i})$	$p(M_i|\mathcal{D})$
M 1	170.2	1 − (1.3 × 10−12)
M 2	102.9	5.8 × 10−30
M 3	142.8	1.3 × 10−12

Posterior parameter distributions

The prior and posterior (empirical) marginal distributions of the parameters influencing the stiffness of the main girder are shown in Figure 6. In general, the posterior distributions are more narrow and shift to higher values.

Figure 6.

Prior and (empirical) posterior marginal distributions of the stiffness-related parameters of the elastic moduli of the Northern side span, E_N, the v-shaped pier axis 7, E_V7, the main span, E_M, the v-shaped pier axis 8, E_V8, the Southern side span, E_S, and the bridge caps, E_C.

The prior and posterior (empirical) marginal distributions of the parameters governing the mass of the bridge are shown in Figure 7. The posterior distribution of the track mass shifts slightly to lower values, while the distribution of the density of the cap material shifts significantly to lower values.

Figure 7.

Prior and (empirical) posterior marginal distributions of the mass-related parameters: track mass, μ and caps density, ρ_C.

The prior and posterior (empirical) marginal distributions of the parameters related to the springs are shown in Figure 8 on logarithmic scales. No significant change is observed for the shear connection and rotational bearing stiffness. The posterior distribution of the stiffness of the translational springs shifts to higher values.

Figure 8.

Prior and (empirical) posterior marginal distributions of the parameters related to springs: The shear connection between bridge caps and bridge, G_C, the rotational bearing stiffness z, K_B,R,Z, and the translational bearing stiffness y, K_B,T,Y.

The prior and posterior (empirical) marginal distributions of the parameters related to the standard deviation of the prediction errors are shown in Figure 9. There is no change in the standard deviation of the prediction errors of the eigenvalues, but the uncertainty of the standard deviation of the prediction errors of the mode shapes significantly reduces due to the model update.

Figure 9.

Prior and (empirical) posterior marginal distributions of the parameters related to errors: The standard deviation of the prediction errors of the eigenvalues, σ _λ , and the standard deviation of the prediction errors of mode shapes, σ _ϕ .

Model predictions

One goal of model updating (e.g., in a digital twin framework) is to obtain improved model predictions. One way to evaluate the quality of model updating is to evaluate and compare the model predictions based on prior and posterior parameter distributions. For this purpose, 4000 samples of each of the prior and posterior parameter distributions were used to predict the modal properties of the bridge. These results were matched with the observations based on the same logic as described in Algorithm 6 and equation (20).

Table 4 compares the natural frequencies of the prior and posterior model predictions with the natural frequencies f of the observations, showing the average natural frequency $\overline{f}$ of the model predictions and average error ${\overline{ϵ}}_f=(\overline{f}-f)/(f)$ between the predicted and observed natural frequencies of matched modes. Using a threshold value δ* for the maximal discrepancy in the mode matching algorithm may mean that not all predicted modes can be matched to observed modes. To determine the ratio r _M of matched modes, we divide the number of matched modes by the total number of predictions. The discrepancy between the average predicted frequencies and the observed frequencies decreases after model updating. All the posterior predictions of the modes could be matched with the observed modes. For some of the observed modes, only a fraction of the prior predictions could be matched. For these modes, the average errors ${\overline{ϵ}}_f$ are only an estimate.

Table 4.

Comparison of prior and posterior predictions with observed modal properties.

	obs. f	Prior $\overline{f}$	Post. $\overline{f}$	Prior ${\overline{ϵ}}_f$	Post ${\overline{ϵ}}_f$	Prior r M	Post. r M
1	1.97 Hz	1.69 Hz	2.06 Hz	−14.2%	4.6%	100%	100%
2	3.28 Hz	2.74 Hz	3.07 Hz	−16.5%	−6.4%	47%	100%
3	3.38 Hz	2.94 Hz	3.40 Hz	−13.0%	0.6%	78%	100%
4	4.85 Hz	4.32 Hz	4.92 Hz	−10.9%	1.4%	99%	100%
5	5.65 Hz	4.21 Hz	5.44 Hz	−25.5%	−3.7%	58%	100%
6	6.56 Hz	5.44 Hz	6.12 Hz	−17.1%	−6.7%	74%	100%
7	6.66 Hz	5.25 Hz	6.25 Hz	−21.2%	−6.2%	31%	100%
8	7.68 Hz	6.31 Hz	7.82 Hz	−17.8%	1.8%	51%	100%
9	7.97 Hz	6.57 Hz	7.94 Hz	−17.6%	−0.4%	42%	100%
10	8.50 Hz	7.15 Hz	8.13 Hz	−15.9%	−4.4%	21%	100%
11	9.11 Hz	8.46 Hz	9.09 Hz	−7.1%	−0.2%	18%	100%
12	9.70 Hz	8.77 Hz	9.07 Hz	−9.6%	−6.5%	9%	100%

Figure 10(a) and (b) show the matrices of the average MAC values between predicted and observed mode shapes. The average MAC values of some observed modes with both prior and posterior predicted modes are quite low. One hypothesis is that the observed mode shapes are not well identified based on the available vibration velocity data, which result from ambient excitations. The illustrations of the mode shapes in Table 1 also supports this observation.

Figure 10.

Comparison of MAC values. (a) Matrix of the mean MAC value of the prior predicted modes versus the observed modes. (b) Matrix of the mean MAC value of the posterior predicted modes versus the observed modes. (c) Mean MAC values between the identified modes and the prior and posterior predictions and range of MAC values between the identified modes and the other observations.

To examine this hypothesis, the modes identified in the current dataset were compared to those of 11 other datasets. The MAC values of these comparisons show significant variations, further supporting this hypothesis. Figure 10(c) relates the average MAC values between the identified modes and the predicted ones to the range of MAC values between the identified modes and the other 11 observations. For some modes, there is very little congruence between observations with a lower bound of MAC values as low as 0.4.

From a system identification perspective, this result is unsatisfactory and requires further examination. For the model updating the upper bound of the range of MAC values between observations also indicates the upper bound for a plausible accordance between the posterior predictions and the identified modes. It would be implausible for the predictions to match the identified modes better than other observations. Figure 10(c) also shows improvement of the posterior predictions in comparison to the prior predictions.

Performance analysis

The motivation for this work is to provide efficient Bayesian Updating of engineering models by parallelizing aBUS. Theoretically, the reduction in execution time of the improved algorithms presented in this contribution would scale (almost) linearly with the number of processes used. However, the execution time depends on other factors, such as the synchronization of intermediate results, bottlenecks like time-consuming input-output operations, and the constraints imposed by the hardware, operating system, and engineering software used. This observation leaves us with the question of how the parallelized algorithms scale in practice.

In computer science, the performance of parallel programs can be analysed with appropriate measures of measured execution times (Rauber and Rünger, 2010: Chapter 4). To examine our case, we use the performance measure speedup S, which evaluates the reduction in execution time T of a program executed in parallel on N _W worker processes compared to the sequentially executed program as:

$S(N_W)=\frac{T(1)}{T(N_W)}$

(31)

The achievable speedup can be expressed by Amdahl’s law (Amdahl, 1967), which assumes that a fraction f _S , 0 ≤ f _S ≤ 1, of the program must be executed sequentially and the rest 1 − f _S can be executed in parallel. The execution time T(N _W ) of the parallelized program consists of the time of the sequential part and the parallelizable part: T(N _W ) = f _S ⋅ T(1) + (1 − f _S )/N _W ⋅ T(1). The achievable speedup is thus (Rauber and Rünger, 2010: 164):

$S(N_W)=\frac{T(1)}{f_S\cdot T(1)+\frac{1-f_S}{N_W}\cdot T(1)}=\frac{1}{f_S+\frac{1-f_S}{N_W}}\le \frac{1}{f_S}$

(32)

In this contribution, we examine the parallelized variant of aBUS with the case study from the previous section as a real application benchmark. This real application benchmark thus depends on the problem, the engineering software used, and the hardware. In this scenario, we try to estimate the fraction f _S that must be executed sequentially. This estimate serves as an indicator of the overall performance of the parallelized algorithm.

Bayesian Updating was performed on the model class M₂ with the measured data several times. The number N _W of worker processes was varied between 1 and 16 and the execution time was recorded. For each N _W , three tests were performed. A large variance was observed for repeated executions with the same number of worker processes. This variance can be interpreted by the fact that subset simulation and thus aBUS is a stochastic algorithm by nature, which is prone to performance variations.

To infer probable values for f _S , Bayesian updating was performed on the benchmark results. The processor used (AMD Ryzen 7 PRO 5750G CPU) has eight physical cores. The performance law for N _W is likely to be different above and below the number of physical cores. Therefore, only benchmark data with N _W ≤8 was used. For each test, the speedup was calculated with respect to the average sequential time $\overline{T}(1)$ . A uniform $\mathcal{U}(\mathrm{0,1})$ prior distribution was assumed for f _S . The likelihood was formulated based on a prediction error which models the discrepancy between the true speedup and the observed speedup . The prediction error includes thus the observation error (e.g., due to other software, e.g., anti-virus software, using the computer’s resources) as well as the model uncertainty (since Amdahl’s law is unable to capture more nuances of the algorithm’s speedup, e.g., certain thresholds regarding input-output operations). We assume the prediction error of the speedup to follow a normal distribution $\mathcal{N}(\mathrm{0,0.25})$ . The mean of the posterior distribution of f _S is 0.27, with a 90% highest density interval (HDI, (Gelman et al., 2014: 33)) of [0.24, 0.29]. Figure 11 shows the observed speedup and the mean and 90% HDI of the posterior predictions of the speedup using Amdahl’s law, including the prediction errors.

Figure 11.

Speedup over the number of worker processes. The speedup observed in the benchmark tests is shown as dots. The prediction of speedup by Amdahl’s law with the inferred parameter f _S is shown with mean and 90% HDI.

We now know that the expected value of the upper bound of the achievable speedup for the case study is 1/f _S = 3.7. This serves as an indicator of the performance of the parallelized algorithms.t

Discussion

This paper shows how sampling-based Bayesian updating of large engineering models with many uncertain parameters can be performed more efficiently, while also providing an estimate of the model evidence. For this purpose, aBUS combined with subset simulation was parallelized and applied in a case study considering a large concrete railway bridge. In the following, we discuss the algorithmic improvements and case study separately.

Algorithmic improvements

The described failure-robust parallelization of two components of aBUS is a straightforward but important improvement of the algorithm. The strengths of this approach include:

(1) Enhanced computational efficiency: Parallelization of the aBUS algorithm significantly improves computational efficiency. This is crucial for handling large engineering models with many uncertain parameters, enabling more efficient Bayesian updating.

(2) Robustness to computational failures: The design of the parallelization to accommodate and recover from computational failures enhances its robustness.

There is room for improvement in the following areas:

• Chain generation robustness: While the generation of Markov chains as a whole is robustly parallelized, a failure in the evaluation of the last candidate sample of a chain might render previous samples ineffective. This is due to the synchronous nature of the algorithm, which evaluates the acceptance rate of a batch of Markov chains to adjust the scaling parameter λ. In scenarios with a high probability of computational failure, this leads to performance degradation.

• Limit on the number of worker processes: The number of worker processes (N _W ) of Algorithm 5 is currently limited by the number of chains (N _C ) considered for adaptation, which may restrict the scalability of the algorithm in certain applications.

• Influences on scaling parameter adaptation: The adaptation of the scaling parameter λ is influenced by several factors, including the chain length (generally $p_0^{-1}$ ), the number of chains (N _C ) considered for adaptation, and the strategy chosen for the initial value of the standard deviation of the proposal distribution (see (Papaioannou et al., 2015) for a detailed discussion). Both the accuracy and performance of these influences should be examined, as the conditional subset probability p₀ and the target acceptance rate α* were examined by (Betz et al., 2018b).

• Leveraging subset simulation research: There is substantial research on subset simulation that offers critiques and improvements (e.g., Blandfort, 2021; Breitung, 2019; Chen et al., 2022; Zhao and Wang, 2022). Integrating insights from this research could further enhance the performance and accuracy of the Bayesian updating process with parallelized aBUS.

Case study

The case study considering the Maintalbrücke Gemünden bridge was conducted to demonstrate and verify the efficiency of the proposed parallelization of aBUS. By applying this approach, large engineering models can be updated within a reasonable time frame.

The approach was verified by comparing the prior and posterior predictive distributions of the modal properties with the observed data. Based on the posterior predictions, more modes could be clustered and thus matched to the observed modes. The mean natural frequency of the predictions was closer to the observed natural frequencies. The average MAC values slightly decreased.

To match the observed modes with the predicted modes, a heuristic mode pairing algorithm was used. The results differ depending on the chosen heuristic. Investigations using distance measures other than the one defined in equation (20) revealed slightly different posterior distributions for the case study presented in this contribution. In addition, the impact of the threshold value δ* on the mode pairing needs to be examined further. The optimal choice of this threshold achieves a trade-off between (a) correctly matching corresponding modes even if the value of the discrepancy δ _ij is large and (b) preventing the incorrect matching of closely spaced modes. In a related example (Brehm et al., 2010), the results of a sensitivity analysis varied depending on the measure used for mode shape similarity in mode matching. Further research is needed to compare different possible metrics (Morales, 2005; Vacher et al., 2010) and their limitations. For intrusive methods that have access to the stiffness and/or mass matrices of a finite element model, updating schemes without mode matching exist (Das and Debnath, 2020; Huang and Beck, 2015; Yuen et al., 2006).

A total of twelve well-separated modes were employed for the purpose of model updating. As is the case with other structures, closely-spaced modes require modal analysis methods capable of separating them (Au et al., 2021; Jonscher et al., 2023). The uncertainty arising from measurement uncertainty and uncertainty of the modal identification inform the prior probabilistic model of the prediction errors. An extended approach can include the uncertainties of the modal identification into the model updating process (Au and Zhang, 2016).

Different assumptions about the nature of the prediction errors and their correlation structure lead to different results (Koune et al., 2023; Simoen et al., 2013). Considering these different assumptions has not been part of this case study, but should be included in a more thorough application.

In addition, it is important to consider other modelling choices, such as parameter selection, assumptions about physical laws, and the level of modelling sophistication. Sensitivity analysis is a common practice for identifying candidate parameters for updating (see e.g., Brownjohn et al., 2001; Jang and Smyth, 2017; Moravej et al., 2019; Simoen et al., 2015). However, this practice only considers parameter selection and not other modelling choices. For further discussion on modelling choices for model-based SHM, please refer to (Simon et al., 2024). A decision-theoretic approach for model selection is presented in (Kamariotis and Chatzi, 2023).

The current case study did not consider all three sources of influencing factors in modelling and model updating: mode pairing heuristics, assumptions about prediction errors, and modelling choices. It is crucial to consider these factors in real-world scenarios. Using aBUS in a parallelized manner, it is possible to evaluate these influences based on model class selection.

Performance

The main objective of this contribution is to increase the efficiency of aBUS through parallelization. Therefore, the performance of the parallelized version was compared to the conventional unparallelized aBUS in terms of achieved speedup. The benchmark tests conducted in the case study revealed an expected upper bound of 3.7 for the speedup. Although this is an improvement, it falls short of the expected results when analysing the improved algorithms. The non-parallelizable fraction f _S of the execution time that remains is likely to be implementation-specific and/or dependent on other hardware factors besides the number of processes. Therefore, it should be examined in additional applications.

In addition to the speedup, scalability may also be of interest as a performance measure. Scalability refers to how well additional processes improve performance as the problem size increases (Rauber and Rünger, 2010: 165). The motivation for parallelizing is to address Bayesian updating of large engineering models. Therefore, scalability is an important indicator of how well the parallelized algorithms allow for the use of larger engineering models. However, Bayesian updating presents a challenge in defining the problem size, as it could refer to the execution time of the engineering model, the number of random variables, or the number and the variety of observations. Additionally, the problem size is dependent on the formulation of the problem, which is expressed in the probabilistic model class. In engineering structures, the formulation of the problem is determined by the size and type of the structure as well as the sophistication of modelling. As such, the issue of defining the problem size and investigating its effects on the efficiency of the parallelized aBUS approach is beyond the scope of this contribution and requires further examination.

In the case study, a target value of N _S = 4000 samples was used to approximate the posterior distribution of the parameters and the model evidence. Fewer samples might also suffice, but the error in the results must be considered (Betz et al., 2018b).

In addition to the parallelization of the sampling, there are three additional possible approaches to improving the efficiency of aBUS (1) optimization of the MCMC sampler used within each subset, as done by Papaioannou et al. (Papaioannou et al., 2015), and Wang et al. Wang et al. (2019), and the optimization of the sampling-related parameters of the algorithm, such as p₀ and N _S , as done by Betz et al. (Betz et al., 2018b), (2) reduction of model evaluation time by constructing surrogate models, such as artificial neural networks (Giovanis et al., 2017), or adaptive kriging (Wang and Shafieezadeh, 2020, 2023; Yuan et al., 2022, 2024, 2025), and (3) data reduction techniques (Ur Rehman et al., 2016) to remove data with redundant information and thus reduce the number of likelihood evaluations.

Surrogate models bear their own challenges in terms of evaluating trade-offs between different types of surrogate models, (Alizadeh et al., 2020; Bhosekar and Ierapetritou, 2018), robustness in the presence of failed simulations (Palar et al., 2019), validation (Bhosekar and Ierapetritou, 2018; Palar et al., 2019), and accuracy (Alizadeh et al., 2020; Palar et al., 2019).

The current contribution improves the efficiency of aBUS through parallelization and addresses the issue of failed simulations. Future research should consider the combination of parallelization and surrogate modelling.

Concluding remarks

In this contribution, we demonstrate that parallelizing aBUS enables Bayesian updating of a large engineering model within a reasonable time. This is a contribution to the field of model-based SHM and has implications for model updating of large engineering structures in general.

Further research is needed in three areas. First, the proposed method requires control parameters and employs heuristics that influence aBUS, both of which should be subject to scrutiny.

Second, the algorithm improved in this contribution fulfils the aims of efficiency, high-dimensional parameter spaces, and providing an estimate of the model evidence. However, further evaluation is necessary to assess its performance and limitations in these three areas. It is crucial to compare it with other methods, as done by (Smith et al., 2023) in their comparison of sampling-based methods. Such a comparison should include the use of surrogate models. An example of comparing surrogate models as part of structural reliability analyses is provided in (Moustapha et al., 2022).

Third, there are significant gaps between the research on model updating and its adoption by practitioners (Ye et al., 2022). Our contribution lowers barriers to adoption by making Bayesian updating of large engineering models feasible in engineering practice. However, further research is necessary to validate the quality of updated models and their modelling procedures (Simon et al., 2024).