Probabilistic Updating of Structural Models for Damage Assessment Using Approximate Bayesian Computation

Abstract

A novel probabilistic approach for model updating based on approximate Bayesian computation with subset simulation (ABC-SubSim) is proposed for damage assessment of structures using modal data. The ABC-SubSim is a likelihood-free Bayesian approach in which the explicit expression of likelihood function is avoided and the posterior samples of model parameters are obtained using the technique of subset simulation. The novel contributions of this paper are on three fronts: one is the introduction of some new stopping criteria to find an appropriate tolerance level for the metric used in the ABC-SubSim; the second one is the employment of a hybrid optimization scheme to find finer optimal values for the model parameters; and the last one is the adoption of an iterative approach to determine the optimal weighting factors related to the residuals of modal frequency and mode shape in the metric. The effectiveness of this approach is demonstrated using three illustrative examples.

1. Introduction

Structural health monitoring (SHM) refers to the damage detection and characterization of engineering structures, which is very important for guaranteeing the safety and serviceability of structures. In recent years, this research field has gained much attention and has made many ground-breaking improvements. The SHM approaches can be roughly divided into two sorts: local approaches and global approaches. The former ones refer to the methods of assessing the highly localized behavior of structures (such as strain at a single location), which includes many local forms of non-destructive evaluation (NDE) [1,2,3]. The latter ones refer to the methods of assessing the global behavior of structures (such as modal parameters), which includes many vibration-based damage assessment methods [4]. A large portion of these damage assessment methods are built using model updating: based on the experimental (modal) data, structural damages are detected, located, and quantified by calibrating the stiffness parameters of structural model [5].

Structural model updating based on modal data is a process of correcting the initial model to a better model that reflects the dynamic behavior of the structure. Besides damage assessment, structural model updating has also been widely used in model validation, response prediction, structural control, etc. [6]. Due to a variety of uncertainties such as modeling error and measurement noise, the discrepancy between the theoretical finite element (FE) model and the actual structure is inevitable, so it is often necessary to update the initial/prior FE model. This topic has become a research hotspot in recent years [7,8,9,10,11,12,13,14,15]. The methods of model updating can be broadly divided into two types: deterministic methods and probabilistic methods. In deterministic methods, one’s objective is to find a single set of optimal model parameters by using an optimization technique to minimize a single-objective function [9,16] or multiple objective functions [17,18] that measure the goodness of fit between the quantities of model prediction and measurement, while the probabilistic approaches estimate the statistical distributions of the model parameters (not only the optimal parameters but also their uncertainties) given a family of possible models [19,20]. By providing well-defined formulations, the deterministic methods can effectively obtain a set of optimal values of model parameters for a well-posed identification or updating problem; nevertheless, they usually fail to deliver satisfactory results when many uncertainties are present. In contrast, probabilistic methods have more advantages in dealing with uncertainty because they can provide the probability distribution of model parameters, as well as the plausibility of the model itself based on measured data.

In recent years, probabilistic methods have been more and more popular in model updating, among which the Bayesian method is one of them. In this method, the available information involved in the measured data is utilized for statistical inference by constructing a suitable likelihood function [7,19]. In the updating framework of the Bayesian approach, the posterior distribution of model parameters is expressed as a product of the prior distribution and the likelihood function. The likelihood function represents the probability of obtaining measured data based on a given stochastic prediction model. Due to a lack of prior information, the prior distribution is often a non-informative distribution, so the likelihood function is crucial to the posterior distribution. If different probability models are assumed in the likelihood function, the subsequent posterior distribution will be different. Usually, the likelihood function is formulated as normal distribution based on the maximum entropy principle [19]. However, challenges still exist in construction of the likelihood function. For some complex models, the explicit expression of the likelihood function may be difficult or even impossible. Another difficulty is the correlation formula of the model parameters, i.e., it is difficult to explicitly express these correlations in the joint probability density function (PDF) of the likelihood function [21]. For finite-element model updating, the likelihood is a probabilistic model for the error between predictions and observations. Most often, a zero-mean uncorrelated Gaussian prediction error is assumed, but this assumption may be questionable in some cases. For model updating in the time domain, the prediction errors of time histories may be non-Gaussian, and they may also show temporal correlation in adjacent time instants; on the other hand, the prediction errors corresponding to different measurements may be correlated as well. For model updating in the frequency domain, the components of a mode shape may show spatial correlation when densely populated sensor grids are used. The assumption of prediction error correlation highly influences the results of Bayesian model updating, and it is a challenging task for us to choose a suitable correlation structure of prediction errors in the likelihood function.

In order to solve above problems, approximate Bayesian computation (ABC) was well-conceived to steer clear of the explicit expression of likelihood function, using a simulation-based approach instead [22]. Thus, it is a likelihood-free Bayesian method. ABC does not directly use the observed data to make statistical inference according to Bayes’ theorem, but it obtains the model prediction output through sampling and simulation, and the discrepancy between the predicted output and the observed data is acceptable under certain metric. Therefore, an approximate likelihood function is defined as the probability that the simulated predicted output falls into the neighboring region of the observed data vector. Obviously, the accuracy of this approximation depends on the selected tolerance level $ϵ$. In theory, when $ϵ$ → 0, the algorithm gives the exact posterior distribution, but this cannot be done in practice. To get a good approximation of the posterior distribution in the ABC method, the tolerance parameter ϵ has to be small enough so that only the predicted output very close to the observed data vector is accepted. As a result, this becomes a simulation problem for rare events. If a Monte Carlo simulation is used, a large number of model prediction outputs must be calculated as candidate samples in order to obtain an acceptable sample size in the given approximate region of observed data. Jin and Jung presented a feasibility study on structural model calibration by using a Markov chain Monte Carlo (MCMC)-based ABC [23]. Fang et al. adopted ABC incorporated with the Metropolis–Hastings sampling (MHS) algorithm and stochastic response surface (SRS) for probabilistic damage identification [24]. Abdessalem et al. presented ABC sequential Monte Carlo (SMC) [25] and ABC ellipsoidal nested sampling (NS) [26] as efficient tools for model selection and parameter estimation in structural dynamics. Though MCMC can reduce computation cost, its efficiency is still not very high due to its slow chain initialization and convergence speed. In order to improve the efficiency of the ABC, Chiachio et al. [27] recently put forward the efficient subset simulation technique [28] into the ABC framework, naming it as ABC-SubSim. The basic idea is to take the nested descending sequence in the subset simulation as the approximate region that is increasing close to the observed data vector [29,30].

This paper proposes a new probabilistic structural model updating approach based on the ABC-SubSim for structural model updating using modal data. Obviously, the accuracy of the method depends on the selection of the tolerance level $ϵ$ on the chosen metric. First, an appropriate metric should be chosen to measure the closeness of the predicted modal data vector and the observed modal data vector. The tolerance parameter $ϵ$ should be sufficiently small to make a good approximation of the posterior distribution. However, the question of how small is sufficient enough is also questionable. To this end, we propose some new stopping criteria to cease the subset simulation, where the tolerance parameter $ϵ$ for the metric is not predetermined. On the other hand, by using subset simulation, although the approximation region of the observed data vector for a given tolerance level can be easily reached, a refined optimal value may not be easily obtained due to the stochastic exploration nature of the subset simulation. To this end, we propose a hybrid optimization scheme that combines the subset simulation with a faster local optimization technique to obtain a more refined optimal value of the model parameter vector. The last contribution of this paper is that we adopt an iterative approach to determine the optimal weighting factors relating to the residuals of modal frequencies and mode shapes in the metric.

The remainder of the paper is organized as follows. The ABC-SubSim is reviewed in Section 2, where the choice of some important parameters such as intermediate tolerances, conditional probability and variance of proposal distribution is discussed. The procedures of model updating with ABC-SubSim is presented in Section 3. Section 4 presents a damage assessment approach using model updating. Three illustrative examples are presented in Section 5, and conclusions are drawn in Section 6.

2. Recap of ABC-SubSim

2.1. Approximate Bayesian Computation

Based on a set of measured data $y\in D\subset {ℝ}^l$, the parameter of interest $\theta \in \Theta \subset {ℝ}^d$ in a parameterized model class M can be updated from prior distribution to posterior distribution by Bayesian inference. On the basis of Bayes’ theorem, the posterior PDF $p\left(\theta |y,M\right)$ of the parameters θ in model M is written as follows [19]:

$$p\left(\theta |y,M\right)=cp(y|\theta ,M)p(\theta |M)$$

(1)

where $p\left(\theta |y,M\right)$ is a posterior distribution of $\theta$ after taking measured data $y$ into account, c is a normalizing constant to make the integral of posterior PDF to be one, $p(y|\theta ,M)$ is the likelihood function which represents the distribution of the observed data $y$ conditional on its parameters θ, and $p(\theta |M)$ is a prior distribution of the parameters θ before any data are observed. For a given model, we only care about the parameters in the model, not the selection of the model, so the conditioning on model class M can be removed in the notation. Although this Bayesian framework is powerful for uncertainties treatment, construction of the likelihood function is still a challenge in some cases. Approximate Bayesian computation methods circumvent the explicit expression of the likelihood function by using simulation-based methods. Let $x\in D\subset {ℝ}^l$ denote a dataset sampled from $p(·|\theta ,M)$, the ABC algorithms take the approximate likelihood function $P_{ϵ}\left(y|\theta \right)=\int P\left(x\in N_{ϵ}\left(y\right)|x\right)p(x|\theta )dx$, where $N_{ϵ}\left(y\right)$ is an approximate neighborhood of observed data for a chosen tolerance level $ϵ$ on the given metric ρ and $N_{ϵ}\left(y\right)=\left\{x\in D:\rho \left(\eta \left(x\right),\eta \left(y\right)\right)\le ϵ\right\}$. η(·) is a low-dimensional vector of summary statistics that provides a comparison of the closeness of x and y. From the Bayes’ theorem, the approximate posterior is given by

$$p_{ϵ}\left(\theta ,x|y\right)\propto P(x\in N_{ϵ}\left(y\right)|x)p(x|\theta )p\left(\theta \right)$$

(2)

where $P(x\in N_{ϵ}\left(y\right)|x)={\mathbb{I}}_{N_{ϵ}\left(y\right)}\left(x\right)$ is an indicator function that assigns a value of one when $\rho \left(\eta \left(x\right),\eta \left(y\right)\right)\le ϵ$ and zero otherwise. The marginal approximate posterior:

$$p_{ϵ}\left(\theta |y\right)\propto P(x\in N_{ϵ}\left(y\right)|\theta )p\left(\theta \right)$$

(3)

It should be noted that the quality of the approximated posterior in Equations (2) and (3) depends on an appropriate choice of the metric $\rho$, the tolerance level $ϵ$ and the summary statistic η [31,32].

In order to obtain a good approximate posterior distribution, the tolerance $ϵ$ need to be small enough, which leads to the simulation of rare events. Direct Monte Carlo simulation is inefficient, and although MCMC sampling can improve the efficiency, it is still not very efficient.

2.2. The ABC-SubSim Algorithm

To enhance the efficiency of ABC, a new algorithm called ABC-SubSim was conceived by Chiachio et al. [27]. In this algorithm, the subset simulation technique [28] is incorporated into the ABC framework, and the idea is to select the nested descending sequence in the subset simulation as the approximate region where the simulated dataset gets closer and closer to the observed data (Figure 1).

Let’s define z as $z=\left(\theta ,x\right)$, so that $p\left(z\right)=p(x|\theta )p\left(\theta \right)$. A sequence of nested decreasing regions D_j, j = 1…, m, are defined as

$$D_j=\left\{\left(\theta ,x\right):\rho \left(\eta \left(x\right),\eta \left(y\right)\right)\le {ϵ}_j\right\}$$

(4)

where ${ϵ}_j$ is intermediate tolerance. The sequence of intermediate tolerances ${ϵ}_1,{ϵ}_2,\dots ,{ϵ}_m$, with ${ϵ}_{j+1}<{ϵ}_j$, are chosen adaptively which will be described later. The number of levels m is chosen so that ${ϵ}_m\le ϵ$, a specified tolerance. Thus, it follows that

$$p\left(D_m\right)=p\left(D_1\right){\displaystyle \prod }_{j=2}^{m}p\left(D_j|D_{j-1}\right)$$

(5)

The samples of the first level are drawn using a direct Monte Carlo simulation from the prior distribution, while the samples of higher levels are drawn using a modified Metropolis algorithm (MMA) [28]. Instead of drawing a complete parameter vector candidate from a multidimensional proposal PDF as the original algorithm would do, MMA chooses a univariate proposal PDF for each component of the parameter vector, and the generated component candidate is accepted or rejected separately. It should be noted that although ABC-SubSim and ABC-SMC do have some similarities—they both have a series of intermediate tolerances; they differ in the way of generation of new samples from the samples in the previous level. In the ABC-SMC [25,33], samples are drawn from the previous population with some weights at first, and then they perturb the sample to obtain a new sample according to a perturbation kernel; while in the ABC-SubSim, the new samples are generated using the modified Metropolis algorithm from the seeds provided by the previous level.

2.2.1. Intermediate Tolerances

In the actual implementation of subset simulation, the intermediate tolerances are adaptively chosen so that $p\left(D_1\right)$ and $p\left(D_j|D_{j-1}\right)$, j = 1, …, m are equal to a fixed value P₀. In this way, the value of intermediate tolerance ${ϵ}_j$ is determined as the average of the (NP₀)th and (NP₀ + 1)th values of the set of distances $\rho (\eta (x_{j-1}^n),\eta (y))$, n = 1, …, N, arranged in increasing order, where N is the number of samples in each level. In this algorithm, the number of samples generated in D_j−1 falling in D_j is equal to NP₀, and these NP₀ samples are utilized as seeds to generate NP₀ Markov chains with a length of 1/P₀ samples, in which the new (1/P₀ − 1) samples of each chain are produced by MMA [28]. For the sake of implementation, the numbers NP₀ and 1/P₀ are both integers. Thus, the total number of samples in D_j is also N, of which NP₀ samples are produced in the (j−1)th level.

2.2.2. Conditional Probability

The parameter P₀ determines the number of intermediate domains D_j required to reach the target region D, which thereby influences the simulation efficiency. A smaller value of P₀ connotes that fewer intermediate levels are required to reach D, but it makes for larger number of samples N required at each level to accurately estimate small conditional probabilities. On the other hand, a larger value of P₀, although decreasing the number of samples required for each level, will increase the number of intermediate levels m. Therefore, we need to take a compromise value, and a theoretical study shows that the optimal choice is $P_0^{opt}\approx 0.2$ [34].

2.2.3. Variance of Proposal Distribution

It was observed in [28] that the type of proposal PDFs (Gaussian, uniform, etc.) has no significant impact the efficiency of MMA, while their spread (variance) has a greater impact. In case of Gaussian proposal PDFs, too large or too small variances ${\mathsf{\sigma }}_j^2$ tend to increase the correlation between successive samples. Larger variances may reduce the acceptance rate and increase the number of duplicate MCMC samples. Conversely, smaller variances may result in higher acceptance rates, but because they are closely related, the chain moves too slowly and produces many correlated samples. A theoretical study on the optimal variances for the proposal PDF has shown that the nearly optimal variances can be selected to make the acceptance rate between 30% and 50% at each conditional level [34]. Unfortunately, achieving this target acceptance rate usually requires the user’s prior experience or through inefficient brute force searches. Recently, Vakilzadeh et al. [29] proposed a self-regulating algorithm in which the proposal variance is adjusted adaptively in ABC-SubSim to realize effective sampling in the posterior PDF. By using this algorithm, the proposal variance can be dynamically adjusted to force the average acceptance rate to an ideal target value.

3. Model Updating Using ABC-SubSim

3.1. Choice of Metric Function

The quality of the posterior approximation depends on a suitable choice of the metric $\rho$ and the tolerance parameter $ϵ$. In this subsection, the choice of the metric $\rho$ is presented. In the following subsection, the tolerance parameter $ϵ$ will be discussed. We choose the metric $\rho$ for measuring the closeness of the model predicted modal data vector and the observed modal data vector as follows

$$\rho \left(\eta \left(f_r,{\varphi }_r\right),\eta \left({\widehat{f}}_r^i,{\widehat{\varphi }}_r^i\right)\right)={\displaystyle \sum }_{r=1}^{N_m}\left(w_r^f{\displaystyle \sum }_{i=1}^{N_s}{\left(\frac{f_r\left(\theta \right)-{\widehat{f}}_r^i}{{\widehat{f}}_r^i}\right)}^2+w_r^{\varphi }{\displaystyle \sum }_{i=1}^{N_s}\frac{\Vert a_r^i{\varphi }_r\left(\theta \right)-{\widehat{\varphi }}_r^i{\Vert }^2}{\Vert {\widehat{\varphi }}_r^i{\Vert }^2}\right)$$

(6)

where $w_r^f$ and $w_r^{\varphi }$ are the weighting factors for the residues of the rth modal frequency and the rth mode shape, respectively; $f_r\left(\theta \right)$ is the rth model predicted modal frequency, ${\widehat{f}}_r^i$ is the rth observed modal frequency in the ith data set; ${\varphi }_r\left(\theta \right)$ is the rth model predicted partial mode shape confined on the observed degrees of freedom (DOFs), ${\widehat{\varphi }}_r^i$ is the rth observed mode shape in the ith data set, $a_r^i={({\widehat{\varphi }}_r^i)}^T{\varphi }_r\left(\theta \right)/({\varphi }_r{\left(\theta \right)}^T{\varphi }_r\left(\theta \right))$ is a scaling factor that guaranties that the predicted model mode shape ${\varphi }_r\left(\theta \right)$ is closest to the measured mode shape ${\widehat{\varphi }}_r^i$ at the measured DOFs, and N_m and N_s is the number of modes and number of observed data sets. The determination of the weighting factors $w_r^f$ and $w_r^{\varphi }$ is discussed in details in Section 3.4. It should be noted that besides the metric presented above, one can also choose another function as the metric to measure the closeness of the mode prediction and the observations [9,35,36,37].

3.2. Stopping Criteria

The tolerance parameter $ϵ$ is another important factor that controls the quality of the posterior approximation. In theory, when $ϵ$ → 0, the algorithm can get an accurate posterior distribution, but this could not be achieved in practice. If one wants to get a good approximate posterior distribution in the ABC approach, $\mathsf{ϵ}$ has to be small enough so that the predicted outputs are accepted only if they fall into the nearest neighboring domain centered on the observation data vector. However, how small is sufficient is unknown in advance. By using a different metric, the tolerance level may be quite different to get a satisfactory accuracy. For example, if we use another set of weighting factors for the metric presented in Equation (6), for the same tolerance level $ϵ$, the accuracy of the posterior approximation may be quite different. In contrast, we set the stopping criteria by justifying the convergence of the model parameters through looking at two relative indexes. The new stopping criteria to cease the subset simulation are as follows:

$$R_j=\frac{{\rho }_j^{max}-{\rho }_j^{min}}{{\rho }_j^{min}}\le R^{tol}$$

(7)

$$S_j=\frac{\Vert {\theta }_j^{opt}-{\theta }_{j-1}^{opt}\Vert }{\Vert {\theta }_2^{opt}-{\theta }_1^{opt}\Vert }\le S^{tol}$$

(8)

where ${\rho }_j^{max}$ and ${\rho }_j^{min}$ are the maximum and minimum values of the metric in the jth conditional level, respectively; $R^{tol}$ is the tolerance of R, which can be chosen as a very small value (e.g., 10⁻³); ${\theta }_j^{opt}$ and ${\theta }_{j-1}^{opt}$ are the optimal model parameter vector in the jth and (j − 1)th conditional levels, respectively; and $S^{tol}$ is the tolerance of S, which can be chosen as a relatively larger value (e.g., 10⁻²). If both the stopping criteria in Equations (7) and (8) are satisfied, the subset simulation is stopped. R_j is defined to measure the diversity of the jth conditional level, and S_j is defined as the improvement of the jth conditional level. In such a way, although we have to specify two tolerance values, the meanings behind these two values are more clear and straightforward, and this approach is more operable in practice. It should be noted that this is only an expedient. To completely solve this problem, we still need to deeply study the relationship between likelihood-informed Bayesian inference and likelihood-free approximate Bayesian computation so as to provide a theoretical basis for selecting the appropriate tolerance [38].

3.3. Hybrid Optimization Scheme

Stochastic sampling techniques such as subset simulation have the advantage of efficient posterior sampling capacity, but their ability of finding an optimum or “maximum a-posteriori” is relatively low. The stochastic sampling techniques can reach the region near an optimum point relatively quickly, but it can take many levels (generations) to achieve the convergence to the exact optimum. Therefore, in order to improve the overall search ability, it is necessary to combine a powerful global search sampler with an efficient local search operator to develop a hybrid scheme. For example, Jung and Kim proposed a hybrid genetic algorithm to update their structural model by combining the genetic algorithm and the improved Nelder–Mead’s simplex method [39]. Sun and Betti proposed a hybrid optimization algorithm that combines a modified artificial bee colony (MABC) algorithm and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method for Bayesian model updating [40]. As a global stochastic sampler, the subset simulation is able to quickly reach the approximation region of the posterior for a given tolerance level, and this region stands for the approximate posterior distribution of the model parameter vector θ. However, the optimal value given by the subset simulation may not be the “exact” optimum, and this “exact” optimum can be rapidly refined by a local optimization search operator.

In this study, the BFGS method was used to fine-tune the optimal solution given by the subset simulation to obtain a better solution. Herein, in order to make sure that the hybrid scheme did not disturb the stochastic process of the subset simulation, the BFGS method was just introduced at the final level of the subset simulation. The optimal value at the final level of the subset simulation was taken as the initial guess in the BFGS method. The BFGS method is a gradient-based quasi-Newton method that uses the BFGS approximation of the hessian matrix. The details of the BFGS method can be found in the literature [41,42,43,44]

3.4. Weighting Factors Determination

In the metric function as presented in Equation (6), weight factors are imposed a relative difference between the modal frequency and the mode shape, as well as between different modes, because the measurement precision of these quantities may be different. The identification accuracy of experimental modal frequencies is usually higher than that of modal shapes, so appropriate weights are necessary. In the literature, most of them chose these weighting factors by engineering experience [9,45]. A rigorous theoretical study on the optimal weighting factors in the weighted modal residuals metric for structural model updating was conducted by Christodoulou and Papadimitriou [46]. The results showed that the optimal weight of the modal residual is asymptotically inversely proportional to the value of the modal residual at the optimal model parameters when there is a large amount of measured data. The weighted modal residuals metric in Equation (6) can be rewritten as:

$$J\left(\theta ,w\right)={\displaystyle \sum }_{i=1}^n{w}_i{J}_i\left(\theta \right)$$

(9)

where $J_i\left(\theta \right)$ is the residual of a specific modal group and $w_i$ is its weighting factor. The optimal weights in Equation (6) are given by

$${\widehat{w}}_i=\frac{{\alpha }_i}{J_i\left({\widehat{\theta }}_{opt}\right)}$$

(10)

where ${\alpha }_i$ is a scale reflecting the proportion of the data volume contained in each modal group to the total data volume in a modal data set. N_m denotes the number of modes, and N_o denotes the number of observed degree of freedom; as such, the total number $N_t=\left(N_o+1\right)N_m$ and ${\alpha }_i=1/N_t$ for modal frequency, while ${\alpha }_i=N_o/N_t$ for mode shape, satisfying ${\sum }_{i=1}^n{\alpha }_i=1$. However, the optimal weighting factors given in Equation (10) is dependent on the optimal model parameter vector ${\widehat{\theta }}_{opt}$, which is unknown in advance. Therefore, one should solve it in an iterative approach. The initial value for the weighting factors can be calculated using Equation (10) when taking the nominal value of the model parameter vector as its optimal value ${\widehat{\theta }}_{opt}$. By using ABC-SubSim and the BFGS method, a new optimal value for ${\widehat{\theta }}_{opt}$ is obtained, and then the new weighting factors can be calculated using Equation (10) again. These procedures are repeated several times until the convergence satisfied.

3.5. Summary of the Procedures

The proposed algorithm for probabilistic model updating can be schematically described below:

• Calculate the initial values for the weighting factors used in the metric as presented in Equation (6) using Equation (10) when taking the nominal value of the model parameter vector as its optimal value.
• Using direct Monte Carlo simulation to get N samples following the prior distribution of the model parameter vector.
• Evaluate the metric and sort the samples in order according to their metric values.
• Select NP₀ samples which corresponding to the first NP₀ smallest metric values. These samples will be taken as the seeds for Markov chains generation.
• Generating NP₀ Markov chains with length of 1/P₀ for each chain using a self-regulating modified Metropolis algorithm [29].
• Check whether the stopping criteria proposed in Equations (7) and (8) is satisfied. If it is not satisfied then repeat the steps iii–vi. If it is satisfied then output the posterior samples.
• Take the optimal value in the posterior samples as the initial guess for a local optimization search operator, and conduct the local optimization search to obtain a fine optimum.
• Using the fine optimum to calculate the weighting factors using Equation (10).
• Check the convergence of the weighting factors. If the convergence satisfied, then the algorithm ends. If not, take the newly calculated weighting factors for the metric and repeat the steps i–ix.

A flowchart of the proposed algorithm is given in Figure 2, and a pseudocode is provided in Appendix A.

4. Damage Assessment Based on Model Updating

4.1. Structural Model Parameterization

For model updating and damage assessment of linear structures, the structure is usually partitioned into $N_{\theta }$ substructures, and a substructure can be a group of elements or just one single element. The global stiffness matrix is parameterized as follows:

$$K\left(\theta \right)=K_0+{\displaystyle \sum }_{n=1}^{N_{\theta }}{\theta }_n{K}_n$$

(11)

where K₀ is the global stiffness matrix before updating, K_n is the nominal stiffness matrix of the nth substructure, and θ_n is a perturbation scaling factor used to change the nominal substructure stiffness K_n so that the updated model is more in line with the actual behavior of structure. When a substructure is damaged, the optimal value of the parameter θ_n decreases. The mass matrix can be parameterized in a similar substructural manner, but mass updates were not considered here, because the main interest of this paper was damage assessment, and the mass was assumed to be known accurately enough and to not change with the occurrence of damage.

4.2. Probabilistic Damage Assessment

According to a two-stage procedure for structural damage assessment using Bayesian model updating proposed by Ching and Beck [47], experimental modal parameters are identified from measured data at intact and damaged state in the first stage; in the second stage, the PDFs of structural stiffness parameters are updated using the Bayesian approach before and after damage, and the updated PDFs are used to evaluate the probability of damage of each substructure. Assuming that the stiffness parameters in the intact state and the possible damaged state are conditionally independent (and by using Gaussian asymptotic approximation), the probability of stiffness reduction of parameter θ_n which has been reduced by certain fraction d compared with the undamaged state can be obtained:

$$P_n^{dm}\left(d\right)=P[(1+{\theta }_n^{pd})<\left(1-d\right)\left(1+{\theta }_n^{ud}\right)]\approx \mathsf{\Phi }\left(\frac{\left(1-d\right)\left(1+{\widehat{\theta }}_n^{ud}\right)-\left(1+{\widehat{\theta }}_n^{pd}\right)}{\sqrt{{\left(1-d\right)}^2{\left({\widehat{\sigma }}_n^{ud}\right)}^2+{\left({\widehat{\sigma }}_n^{pd}\right)}^2}}\right)$$

(12)

where Φ( ) is the cumulative distribution function of the standard normal distribution; ${\widehat{\theta }}_n^{ud}$ and ${\widehat{\theta }}_n^{pd}$ stand for the a maximum a posterior probability (MAP) estimate of the stiffness parameters at the undamaged state and possibly damaged state, respectively; and ${\widehat{\sigma }}_n^{ud}$ and ${\widehat{\sigma }}_n^{pd}$ denote the corresponding standard deviations of the stiffness parameters.

5. Application Examples

5.1. Model Updating of a Shear Building

An eight-story shear building was considered in this example. The mass per floor was assumed to be 2 × 10⁴ kg, while the nominal inter-story stiffness for all floors was assumed to be 15 × 10⁶ N/m. The actual value of the inter-story stiffness for the first four stories had a 20% reduction of the nominal ones, while that for the remaining four stories had a 20% increase of the nominal ones. We denoted the model parameters vector θ as the increment ratio of the inter-story stiffness, compared to their nominal ones, so their true values were θ = {−0.2, −0.2, −0.2, −0.2, 0.2, 0.2, 0.2, 0.2}. Twenty sets of modal data were assumed to be available for the structural model updating, where the incomplete measurements were considered; that is, only five floors (1, 3, 5, 7, and 8) were measured, and only the first four modes were identified. For the simulated modal data, zero-mean Gaussian noise with a diagonal covariance matrix ${\mathsf{\Sigma }}_{\epsilon }$ was added to the exact calculated modal frequencies and mode shapes. Two cases of noise level were considered in this study. In noise level 1, a 1% coefficient of variation (CV) noise was added to the modal frequencies, while a 5% of CV noise was added to the mode shapes for all modes, a reasonable value for typical modal tests. Noise level 2 was a more severe noise level situation, where 2% of CV noise was added to the modal frequencies, while 10% of CV noise was added to the mode shapes for all modes. In order to study the effect of the BFGS method on the accuracy of the results in the fine-tuning stage, the results before and after the fine-tuning were compared. Moreover, two sets of different values for the stopping criteria parameters were used to study their impact on the final results. In tolerance level 1, we set the stopping criteria parameters R_tol = 10⁻³ and S_tol = 10⁻², while in tolerance level 2, stricter stopping criteria parameters were used: R_tol = 10⁻⁴ and S_tol = 10⁻³. Totally, there were four cases in the simulations: case 1 with noise level 1, fine-tuning, and tolerance level 1; case 2 with noise level 2, fine-tuning, and tolerance level 1; case 3 with noise level 1, without fine-tuning, and tolerance level 1; and case 4 with noise level 1, fine-tuning, and tolerance level 2.

We used the method proposed in Section 3 to update the stiffness parameters vector θ. In the subset simulation, we chose the number of samples in each level N = 2000 and the level conditional probability P₀ = 0.2. The variance in Gaussian proposal distribution was regulated using a self-regulating algorithm proposed by Vakilzadeh et al. [29].

The optimal weighting factors in the metric are obtained using an iterative approach. The evolution of the optimal weighting factors with the iteration steps is shown in Figure 3. It can be seen that the weighting factors converge very fast. It also can be seen that the weighting factors for the modal frequency are larger than those for the mode shape, which means that more uncertainness are contained in the mode shape. It agrees with our preset simulated modal data, in which we added lager noise to the mode shapes than the modal frequencies.

The convergence of the metric and the model parameters vector with the subset levels when using the optimal weighting factors in the metric for case 1 is shown in Figure 4, where R and S are the quantities defined in Equations (7) and (8) to measure the convergence of the metric and the model parameters vector, respectively, and log10 (ρ_min) is the logarithm of the minimum value of the metric in each level. It can be seen that the convergence parameters (R, S) nearly became flat from the 20th level, and the subset simulation stopped at the 29th level based on our preset tolerances of R_tol and S_tol.

The scatter plot of the posterior and prior samples for case 1 are shown in Figure 5, where the prior distributions of the model parameters were taken as uniform distribution within the interval [−1,1]. It can be seen that the posterior samples were confined in very small regions compared with the prior samples. This means that when using Bayesian inference from the informative measured modal data, the uncertainties of the model parameters contained in the posterior distribution were reduced much more compared to the prior distribution. The histograms of the prior and posterior samples of the model parameters are provided in Figure 6. It can be seen that the posterior samples were concentrated in a smaller range compared to the prior samples, leading to a peaked posterior PDF around the actual values of −0.2 or 0.2.

The updated results of the model parameters are tabulated in

Table 1. The identified results of the model parameters.

	θ1	θ2	θ3	θ4	θ5	θ6	θ7	θ8
True	−0.2000	−0.2000	−0.2000	−0.2000	0.2000	0.2000	0.2000	0.2000
Optimal 1	−0.2052	−0.2085	−0.1968	−0.2050	0.2002	0.1968	0.2052	0.2000
CV 1	0.0059	0.0072	0.0063	0.0084	0.0126	0.0165	0.0173	0.0063
Optimal 2	−0.1930	−0.2058	−0.2067	−0.2044	0.2184	0.2144	0.2104	0.2011
CV 2	0.0135	0.0157	0.0120	0.0165	0.0240	0.0336	0.0361	0.0120
Optimal 3	−0.2121	−0.2007	−0.2050	−0.1920	0.2167	0.2159	0.1799	0.1848
CV 3	0.0059	0.0072	0.0063	0.0084	0.0126	0.0165	0.0173	0.0063
Optimal 4	−0.1992	−0.1969	−0.2005	−0.2102	0.1941	0.1863	0.2043	0.2029
CV 4	0.0005	0.0008	0.0004	0.0005	0.0005	0.0015	0.0011	0.0010

Note: Optimal ⁱ and CV ⁱ denote the optimal values and the coefficients of variation of the model parameters identified in case i (i = 1,2,3,4).

. It can be seen that the optimal values identified in case 1 were very close to their true values. Though the accuracy of the results in case 2 was not as good as that in case 1, it was also acceptable. On the other hand, it can be seen that the coefficients of variation (CV) of the model parameters identified in case 2 were larger than those in case 1, which means that more uncertainties were involved in the identified model parameters based on the modal data in case 2. Comparing case 3 and case 1, it could be seen that the accuracy of the optimal values was improved through fine-tuning with the BFGS method. The coefficients of variation were reduced in case 4, which means that the uncertainty could be reduced by setting a stricter tolerance level for the stopping criteria.

5.2. Damage Assessment of a Truss Structure

A 21-bar planar truss model is shown in Figure 7, and the structure was simply supported with a span of 8 m. The cross sectional area of each bar was 100 mm². The mass density was 7850 kg/m³. The actual value of the elastic modulus was 190 GPa, and the nominal value was set to be 200 GPa. This means that a stiffness reduction factor of 5% was assumed for each bar element in the undamaged state, while additional stiffness reductions of 20% and 40% were assumed for bar 20 and 21, respectively. This means that the damage extents of bar #20 and #21 at the damage state were 21% (0.2/0.95) and 42% (0.4/0.95), respectively. The structure had 12 nodes, 21 elements, and 21 DOFs. It was assumed that only the first eight modes were available, and only 11 DOFs were measured (marked with arrows in Figure 7). Some nodes were only measured in the vertical direction, and some nodes were only measured in the horizontal direction. This means that the measurements were incomplete, both in modes and DOFs, which was also the case in real situation. It was also assumed that there were 20 data sets of modal frequencies and mode shapes contaminated with white noise of 1% and 5%, respectively.

Using the model updating method proposed in Section 3, the updated results of the perturbation scaling factors of stiffness θ are shown in Figure 8, both in the undamaged and damaged states. It can be seen that the optimal values of the updated parameters were very close to their true values. Their standard deviations were also obtained by the posterior samples that corresponded to the samples in the last level of subset simulation. From the figure, it can be seen that bars #6 and #12 had relatively larger standard deviations than the others, which means that more uncertainties were involved in these two bars based on the current sensor configuration. This reflects the fact that the information contained in the monitoring data, which was critical for damage detection, was dependent on the sensor layout. This was a sensor placement optimization issue, which is not the subject of this article, so readers can refer to some literature on this subject [48,49].

Using the damage assessment method proposed in Section 4, the probability curve of damage extent for each bar is shown in Figure 9. This figure clearly shows the probability of damage extent defined in Equation (12). It can be seen that there were damages with extents of 21% and 42% in bar #20 and #21, respectively, which agreed well with our preset values.

5.3. Damage Assessment of a Simply Supported Beam

A simply supported steel beam before and after damages at different locations was used for further verification. The beam and the sensor layout (marked with S1~S11) are shown in Figure 10. The total length of the beam was 200 cm, and the span length was 180 cm. The width and the thickness of the beam were 7 and 0.5 cm, respectively. For the finite element modeling of the beam, 15 nodes and 14 elements were used, where two ends, two supports, and eleven sensor locations were set as nodes. Only vertical displacements at the sensor locations were measured, while the rotations were difficult to be measured in the real case. The mass per unit volume was 7850 kg/m³. The true value of the modulus of elasticity was 190 GPa, while the nominal value was set to be 200 GPa, which means a stiffness reduction factor of 5% was assumed for each element in the undamaged state. It was also assumed that 20 data sets of modal frequencies and mode shapes of the first five modes contaminated with white noise of 1% and 5%, respectively, were available. The damages were introduced in elements 5, 8, and 11, where the extents were 10%, 20%, and 30% (compared with undamaged state not nominal model), respectively.

By using the proposed model updating algorithm in Section 3 and the damage assessment method proposed in Section 4, the probability curve of damage extent for each element is shown in Figure 11. It can be seen that there were damages with extents of 10%, 20%, and 30% in elements 5, 8, and 11, respectively, which agreed well with the predefined damaged model.

6. Conclusions

A structural model updating and damage assessment method using the ABC-SubSim algorithm was presented. An iterative procedure to find the optimal weighting factors related to different modal residuals and a new stopping criterion to find a proper tolerance level were proposed for the metric used in the ABC-SubSim. A hybrid optimization scheme was proposed to identify the structural parameters based on incomplete modal data with uncertainty. A high accuracy of the structural model updating is guaranteed by the proposed ABC-SubSim algorithm, even though the noise level of modal data was up to 10%. Three illustrative examples well-demonstrated the effectiveness and reliability of the proposed algorithm. Possible directions of future research related to the ABC method may be a relation study between likelihood-informed Bayesian inference and likelihood-free approximate Bayesian computation, thus leading to a more reasonable selection of key parameters in the ABC method.