# An Improved Finite Element Model Updating Method Based on the Singular Values of Frequency Response Functions

Authors: Hong Yin; Zenghui Wang; Mingming Cao; Zhenrui Peng; Kangli Dong
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Aiming at the problems that Markov chain Monte Carlo algorithm is not easy to converge, has high rejection rate, and is easy to be disturbed by the noise when the parameter dimension is high, an improved model updating method combining the singular values of frequency response functions and the beetle antennae search algorithm is proposed. Firstly, the Latin hypercube sampling is used to extract the training samples. The Hankel matrix is reconstructed using the calculated frequency response functions and is decomposed by singular value decomposition. The effective singular values are retained to represent the frequency response functions. Secondly, according to the training samples and the corresponding singular values, the support vector machine surrogate model is fitted and its accuracy is tested. Then, the posterior probability distribution of parameters is estimated by introducing the beetle antennae search algorithm on the basis of standard Metropolis–Hastings algorithm to improve the performance of Markov chains and the ergodicity of samples. The results of examples show that the Markov chains have better overall performance and the acceptance rate of candidate samples is increased after updating. Even if the Gaussian white noise is introduced into the test frequency response functions under the single and multiple working damage conditions, satisfactory updating results can also be obtained.


---

## Body

## 1. Introduction

In recent years, the deterministic finite element model updating (FEMU) techniques have been widely used in structural parameters identification and health monitoring fields [1, 2]. However, due to the widespread uncertainty in engineering problems, the practicability of the deterministic methods is greatly reduced, and the deterministic methods can not get the confidence of the updated prediction values of the model [3]. Therefore, in order to improve the reliability of the parameters identification and health monitoring results, it is necessary to combine the probability and statistics methods to analyze the uncertain factors in the procedure of FEMU [4, 5].All test data is analyzed in the uncertain FEMU methods, and according to the theorem of large numbers in statistical theory, the statistical characteristics of parameters are estimated indirectly by using the statistical characteristics of dynamic responses. The influence of uncertain factors can be fully considered through the probability and statistics, and the deficiency of deterministic model updating methods can be overcome effectively [6]. The FEMU method based on Bayesian theory is a research hotspot of uncertainty methods. Beck [7] firstly put forward the idea of FEMU based on Bayesian theory and proposed an adaptive Markov chain Monte Carlo (MCMC) algorithm based on standard Metropolis–Hastings (MH) algorithm. The most important feature of this method is that it fundamentally avoids the complex inverse problem by solving the randomly selected parameter samples [8]. However, when the parameter dimension is high, the standard MH sampling algorithm will cause the phenomenon of “sampling stagnation,” and the acceptance rate of samples will decrease drastically [9]. In view of the above shortcomings, Wan et al. [10] proposed the delayed rejection adaptive MCMC algorithm combining Bayesian theory to solve the complex gradient calculation, ill-posedness, and nonunique solution. Peng et al. [11] introduced the maximum entropy method into Bayesian theory and fused the cuckoo algorithm with the standard MH sampling algorithm to improve the acceptance rate of candidate samples. Cheung et al. [12] proposed a hybrid MCMC (HMCMC) algorithm to solve the FEMU problems with high-dimensional uncertain parameters.Compared with the FEMU methods based on modal parameters, the FEMU methods based on FRFs have gradually become the mainstream methods [13]. The FRFs contain a lot of frequency point information, so the selection of frequency point is very important, but there is no unified solution to this problem at present. In [14], the FRFs are transformed into the impulse response functions in the time domain by inverse Fourier transform, and then the Hankel matrix is decomposed by singular value decomposition (SVD), and a limited number of singular values are retained to represent the FRFs for model updating. This method not only retains the advantage of FRFs, but also avoids the difficulty of frequency point selection. In [15, 16], the Hankel matrix of the FRFs is decomposed by SVD, and the appropriate number of larger singular values is retained by the principle of sudden change in the number of extreme value points. The results show that this method has a good ability to reduce the noise of the signal.Under the background of the above theory, in this paper, the FRFs are decomposed by SVD, and the effective singular values are extracted as the characteristic quantity of responses for model updating, which can avoid the complex derivation of likelihood functions of FRFs and the difficulty of frequency point selection in Bayesian model updating methods. At the same time, in order to improve the updating efficiency, the support vector machine (SVM) surrogate model is used to replace the complex finite element model for iterative solution. Then, the method of updating the centroid position of the beetle individual in BAS algorithm is introduced to generate new candidate samples to estimate the posterior probability distribution of the parameters. Finally, the feasibility and good antinoise performance of the proposed method are verified by the plane truss model and the three-dimensional truss model.The rest of this paper is structured as follows: Section2.1 starts out with the basics of the FEMU methods based on Bayesian theory and presents the standard MH sampling algorithm. Section 2.2 introduces the SVD technique for extracting the effective singular values of FRFs. In Section 2.3, the improved MCMC algorithm is introduced in detail, which is also the core of this paper. And in Section 2.4, the basic principle of Latin Hypercube sampling and SVM surrogate model is given. Section 2.5 shows the specific flowchart of the proposed FEMU method. Section 3 presents two examples to verify the performance of the novel FEMU method. The main conclusions are given in Section 4.


## 2. Basic Theory and Improvement

### 2.1. Bayesian Finite Element Model Updating

#### 2.1.1. Model Updating Method Based on Bayesian Theory

The FEMU method based on Bayesian theory combines a priori information (subjective information) and test data (objective information), and the MCMC algorithm is used to infer the posterior probability distribution of the updated parameters [17]. The procedure can be expressed by Bayesian formula as(1)pθ|x=px|θπθ∫θpx|θπθdθ=c⋅px|θπθ,where x is the observation information; π(θ) is the prior distribution of the uncertain parameter vector θ to be updated, which is usually equal to 1; c is the constant factor, called the regularization constant; and p(x|θ) is the conditional distribution under the given parameters, which is usually called the likelihood function. Assuming that N independent tests are carried out and the above likelihood function formula can be expressed as(2)px|θ=12πcovyNexp−Jθ,Jθ=∑i=1Nyi−yθiTcovy−1yi−yθi,where y is the test response vector; y(θ) is the calculated response vector of the model; covy is the covariance matrix of test response, substitute (2) into (1), and the posterior probability distribution of the parameter vector can be further written as(3)pθ|x=c′⋅exp−∑i=1Nyi−yθiTcovy−1yi−yθi,where c′ is a constant independent of θ. In engineering problems, the parameter vector is usually a complex high-dimensional vector and the response usually has no explicit expression. Therefore, the stationary Markov chains are obtained with sampling by using MCMC algorithm to estimate the posterior probability distribution of parameters, so as to avoid the problem of high-dimensional integral operation which may be involved in the above Bayesian formula.


#### 2.1.2. Markov Chain Monte Carlo (MCMC) Algorithm

At present, in the field of Bayesian statistics, the MCMC algorithm has become an effective tool to deal with complex high-dimensional integral problems [18]. Firstly, the definition of Markov chain is given: there is a random sequence (θ(0), θ(1), …, θ(n)) at any time j(j ≥ 0), and the value θ(j+1) of the next time j + 1 in the sequence only depends on the current state of time j and has nothing to do with the historical state (θ(0), θ(1), …, θ(j−1)) before time j, which is called the Markov chain.The basic idea of MCMC algorithm is to obtain a Markov chain with stable distribution by sampling to get samples of the target probability distribution and then make statistical inference based on the valid samples in the convergence phase. In Bayesian model updating, the target probability distribution is the posterior probability distributionp(θ|x) of the parameter vector θ. When the samples of the target probability distribution are obtained by MCMC sampling, the statistical characteristics of the updated parameters can be estimated according to the large number theorem; for example, the mean and variance can be approximately calculated as(4)μθ=∫θpθ|xdθ≈1n∑j=1nθj,σθ2=∫θ−μθ2pθ|xdθ≈1n∑j=1nθj−μθ2,where n is the length of the effective Markov chain, that is, the number of valid samples. Markov chain is produced by sampling according to a certain approach in MCMC method, in which Metropolis–Hastings (MH) sampling is the most widely used MCMC algorithm, and many MCMC algorithms are the extension of standard MH sampling algorithm [19].


#### 2.1.3. Standard Metropolis–Hastings (MH) Sampling Algorithm

The most prominent feature of the standard MH sampling algorithm is to obtain a posteriori samples by sampling and constructing appropriate Markov chains to estimate the statistical characteristics of parameters. The main steps of standard MH sampling algorithm are as follows:(1)
Select the appropriate initialθ0 and make the likelihood function p(θ0|D) > 0, where D is the response matrix.(2)
Select the appropriate proposal distributionq(θ∗, θt). Use the current sample θt to generate a new candidate sample θ∗ obeying the proposal distribution.(3)
Calculate the acceptance probabilityα by the current sample θt and the candidate sample θ∗; the formula of acceptance probability is α = min[1, p(θ∗|D)/p (θt|D)].(4)
Generate random numberμ from the uniform distribution U(0,1). If α>μ, accept the candidate sample θ∗, θt+1 = θ∗; otherwise, reject the candidate sample θ∗, θt+1 = θt.(5)
Repeat steps (2) to (4) until the set sampling times are reached or the variance of samples is less than the set value. Terminate the iteration to generate a convergence sequence, the Markov chain.


### 2.2. Feature Extraction of FRFs Based on SVD

#### 2.2.1. Theory of Frequency Response Functions

For ann-degree-of-freedom damping system, the basic dynamic equation can be expressed as Mx¨+Cx˙+Kx=F, where M, C, and K are the mass, damping, and stiffness matrix, respectively; F is system excitation force vector; and x is response displacement vector. Under the action of simple harmonic excitation, the steady-state response X(ω) in frequency domain can be obtained by Fourier transform: X(ω) = H(ω)F(ω), where H(ω) is the FRFs; ω is the excitation frequency. The FRFs matrix can be expressed as(5)Hω=M−iCω−Kω2−1.


#### 2.2.2. Theory of Singular Value Decomposition

Reconstruct the Hankel matrix from the original signalX = [x1, x2, …, xN] with noise, the Hankel matrix A can be expressed as(6)A=x1x2⋯xnx1×τ+1x1×τ+2⋯x1×τ+n⋮⋮⋮⋮xm−1τ+1xm−1τ+2⋯xm−1τ+n,where A is the Hankel matrix of order m × n; τ is the delay step and generally takes 1; m is the embedded dimension, when N is even, m = N/2; or m takes the median [20].First of all, the Hankel matrixA is decomposed by SVD, and then the appropriate threshold value λ is selected, and the first r larger singular values are selected to represent the signal characteristics according to the threshold, while the latter k-r singular values are set to 0. This procedure can be expressed as A=Ur∑rVrT, where Ur and VrT are orthogonal matrices; Σr = diag(σ1, σ2, …, σk) is a diagonal matrix; σi (i = 1,2, …, k) are singular values of matrix A, σ1 ≥ σ2 ≥ ··· ≥σk ≥ 0, k = min[(m − 1) × τ + 1, n].As a nonlinear filtering method, the SVD based on Hankel matrix is used to eliminate random noise in the signal. In this paper, the principle of sudden change in the number of extreme value points proposed in [15, 16] is used to test the rationality of the number of effective singular values, and the retained effective singular values are used to represent the FRFs to update the model.


### 2.3. Proposed MCMC Algorithm

#### 2.3.1. Beetle Antennae Search (BAS) Algorithm

The BAS algorithm proposed by Jiang et al. [21] is an intelligent optimization algorithm inspired by the beetle’s foraging behavior. This algorithm is divided into two stages: search and detection.The search stage is based on the beetle’s food-searching behavior and can be expressed asxl = xi + di·b, xr = xi − di·b, where xl and xr represent the position coordinate of the left and right antennae of the beetle individual, respectively; di is the distance between the two antennae of the beetle individual during the ith search, also called the perceptual distance, which determines the search space of the algorithm. The perceptual distance gradually attenuates with the increase of search time, which makes the beetle individual approach to the optimal value adaptively; b is the unitized random direction vector.The detection stage updates the centroid position of the beetle individual, which can be expressed asxi = xi−1 + δi·b·sign(F(xl, X)−F(xr, X)), where F is the objective function, that is, the food odor concentration perceived by the left and right antennae of the beetle individual; X is the location coordinate; δi is the step length at the ith search, which is very important for the performance of the algorithm. An exponential attenuation model is introduced to constantly update the values of δ and d, and the attenuation model can be expressed as(7)di=d0⋅adi/Hd,δi=δ0⋅aδi/Hδ,where d0 and δ0 are the initial values of perceptual distance and search step, respectively; ad and aδ are attenuation coefficients; i is the number of search times; Hd and Hδ are attenuation constants.


#### 2.3.2. Improved MCMC Algorithm Combining BAS Algorithm

In the procedure of Bayesian model updating, the phenomenon of “sampling stagnation” is easy to occur in standard MH sampling algorithm when the parameter dimension is high. At the same time, because the BAS algorithm does not need complex gradient calculation, the iterative speed is fast and the beetle individual can adaptively approach the optimal solution according to the objective function. In view of the above factors, the standard MH sampling algorithm can be improved by the BAS algorithm. From the expression of the likelihood function of (2) and the properties of the exponential function, the solution corresponding to the maximum value of the likelihood function is the solution corresponding to the maximum value of the exponential term [−J(θ)]. In this paper, the exponential term [−J(θ)] of likelihood function is taken as the objective function of BAS algorithm to guide the search of beetle individual. When the candidate sample point is rejected, the method of updating the centroid position of beetle individual is used to carry out local search on this point to generate a new candidate sample point. The specific sampling steps are as follows:Step 1: initialize the Markov chainsθ(t) (t = 1, 2, …, n), construct an appropriate proposal distribution q(θ∗, θt), where t is the number of iterations and n is the total number of iterations set by usersStep 2: assume the sample at the currentt iteration is θt, and a new candidate sample θ∗ is generated by standard MH samplingStep 3: calculate the acceptance probabilityα and generate random number μ from the uniform distribution U(0, 1). If α > μ, the candidate sample θ∗ is accepted, that is, θt+1 = θ∗; or θ∗ is not rejected immediately, and proceed to the next stepStep 4: fuse the idea of updating the centroid position of the beetle individual(a)
Initialize the parameters of BAS, and take the candidate sampleθ∗ as the initial centroid position of the beetle individual.(b)
Take the exponential term of the likelihood function expression in standard MH sampling algorithm as the search objective function to guide the iterative optimization of the beetle individual.(c)
Stop the search and calculate the acceptance rateαζ according to the centroid position θζ of the beetle individual when the set number of searches is reached (same as the standard MH algorithm to judge the acceptance status: if α > μ, θt+1 = θζ; otherwise, θt+1 = θt).Step 5: Repeat step 2 to step 4 until a convergence sequence is generated.Because the Gaussian distribution is widely used in engineering and its calculation is simple, the Gaussian distribution is selected as the proposal distribution in this paper.


### 2.4. Construct Support Vector Machine (SVM) Surrogate Model

#### 2.4.1. Latin Hypercube Sampling (LHS)

In order to ensure the randomness and uniformity of the training set and test set samples, the LHS is introduced to obtain the samples to construct the surrogate model. LHS is a multidimensional stratified random sampling method, and the sampling points have good distribution uniformity and representativeness [22]. It is assumed that q-group samples need to be taken within the prior range of m-dimensional model parameters. The specific sampling steps are as follows:(1)
Them prior ranges of m-dimensional model parameters [Ai, Bi] (i = 1, 2, …, m) are divided into q cells, which can be recorded as [Aij, Bij] (i = 1, 2, …, m) (j = 1, 2, …, q), resulting in a total of m × q cells.(2)
A sample point is randomly selected in any intercell [Aij, Bij] and marked as αij sampling, and a total of m × q sample points are generated to form a matrix Γ:(8)Γ=α11α12⋯α1qα21α22⋯α2q⋮⋮⋮⋮αm1αm2⋯αmq.(3)
By randomly sorting each row vector [αi1, αi2,·…, αiq] (i = 1, 2, …, m) in the matrix, the resulting vector is [βi1, βi2, …, βiq] (i = 1, 2, …, m), and the matrix Ф is obtained:(9)Φ=β11β12⋯β1qβ21β22⋯β2q⋮⋮⋮⋮βm1βm2⋯βmq.Each column vector in the matrixФ is a group of samples, and a total of q groups of samples are obtained.


#### 2.4.2. SVM Surrogate Model

The surrogate model fits the complex relationship between the FE model parameters and the response by an explicit function. The surrogate model can replace the complex FE model for iterative calculation so as to greatly improve the FEMU efficiency [23]. Clarke et al. [24] compared the performance of SVM with Radial Basis Function (RBF), Neural Network (NN), Kriging model, and the Response Surface surrogate model by engineering examples. The results show that the robustness and accuracy of SVM surrogate model is better than the other three models. So, the SVM is selected as the surrogate model in this paper. When fitting the SVM, the LHS is used to obtain the training set samples and test set samples, and then the corresponding singular values of FRFs are obtained. These data are normalized to avoid error caused by large orders of magnitude.The SVM surrogate model [25] can be expressed as f(s) = ω·ϕ(s) + b, where ϕ(·) is a nonlinear function, which maps the data S to a high-dimensional linear space; ω is a weight vector; b is a deviation.The optimization problem can be expressed as(10)minw,b,ξ12w2+C∑i=1lξi+ξi∗,s.t.yi−ω⋅ϕsi−b≤ε+ξi∗,ω⋅ϕsi+b−yi≤ε+ξi∗,ξi≥0,ξi∗≥0,i=1,…,l,where C is the penalty coefficient; ξi and ξi∗∗ are the relaxation variables, which ensure that the constraint conditions can be satisfied. The Lagrange multiplier method is introduced to transform the above quadratic programming problem to a dual problem, and the regression estimation function can be expressed as(11)fs=∑i=1nαi−αi∗Ksi,si∗+b,where αi and αi∗ are the Lagrange multipliers; K(si, si∗) is the kernel function. Here, the widely used RBF is selected as the kernel function.


### 2.5. Specific Procedure of Model Updating

The training set samples and test set samples are generated by the LHS. The modal participation variation coefficient criterion [26] is introduced to select the appropriate location of excitation points. The modal kinetic energy method [27] is introduced to select the appropriate locations of measuring points. The calculated FRFs are decomposed by SVD to extract the effective singular values according to the principle of sudden change in the number of extreme value points. According to the input parameters and the corresponding singular values, the SVM surrogate model is constructed and its accuracy is tested. Then the SVM surrogate model which meets the accuracy requirement is introduced into the proposed algorithm to iteratively solve the model updating problem, and the mean values of parameters are estimated according to the statistical characteristics of the convergent Markov chains. The specific model updating procedure is shown in Figure 1.Figure 1
Flow chart of FEMU through the proposed MCMC algorithm.


## 2.1. Bayesian Finite Element Model Updating

### 2.1.1. Model Updating Method Based on Bayesian Theory

The FEMU method based on Bayesian theory combines a priori information (subjective information) and test data (objective information), and the MCMC algorithm is used to infer the posterior probability distribution of the updated parameters [17]. The procedure can be expressed by Bayesian formula as(1)pθ|x=px|θπθ∫θpx|θπθdθ=c⋅px|θπθ,where x is the observation information; π(θ) is the prior distribution of the uncertain parameter vector θ to be updated, which is usually equal to 1; c is the constant factor, called the regularization constant; and p(x|θ) is the conditional distribution under the given parameters, which is usually called the likelihood function. Assuming that N independent tests are carried out and the above likelihood function formula can be expressed as(2)px|θ=12πcovyNexp−Jθ,Jθ=∑i=1Nyi−yθiTcovy−1yi−yθi,where y is the test response vector; y(θ) is the calculated response vector of the model; covy is the covariance matrix of test response, substitute (2) into (1), and the posterior probability distribution of the parameter vector can be further written as(3)pθ|x=c′⋅exp−∑i=1Nyi−yθiTcovy−1yi−yθi,where c′ is a constant independent of θ. In engineering problems, the parameter vector is usually a complex high-dimensional vector and the response usually has no explicit expression. Therefore, the stationary Markov chains are obtained with sampling by using MCMC algorithm to estimate the posterior probability distribution of parameters, so as to avoid the problem of high-dimensional integral operation which may be involved in the above Bayesian formula.


### 2.1.2. Markov Chain Monte Carlo (MCMC) Algorithm

At present, in the field of Bayesian statistics, the MCMC algorithm has become an effective tool to deal with complex high-dimensional integral problems [18]. Firstly, the definition of Markov chain is given: there is a random sequence (θ(0), θ(1), …, θ(n)) at any time j(j ≥ 0), and the value θ(j+1) of the next time j + 1 in the sequence only depends on the current state of time j and has nothing to do with the historical state (θ(0), θ(1), …, θ(j−1)) before time j, which is called the Markov chain.The basic idea of MCMC algorithm is to obtain a Markov chain with stable distribution by sampling to get samples of the target probability distribution and then make statistical inference based on the valid samples in the convergence phase. In Bayesian model updating, the target probability distribution is the posterior probability distributionp(θ|x) of the parameter vector θ. When the samples of the target probability distribution are obtained by MCMC sampling, the statistical characteristics of the updated parameters can be estimated according to the large number theorem; for example, the mean and variance can be approximately calculated as(4)μθ=∫θpθ|xdθ≈1n∑j=1nθj,σθ2=∫θ−μθ2pθ|xdθ≈1n∑j=1nθj−μθ2,where n is the length of the effective Markov chain, that is, the number of valid samples. Markov chain is produced by sampling according to a certain approach in MCMC method, in which Metropolis–Hastings (MH) sampling is the most widely used MCMC algorithm, and many MCMC algorithms are the extension of standard MH sampling algorithm [19].


### 2.1.3. Standard Metropolis–Hastings (MH) Sampling Algorithm

The most prominent feature of the standard MH sampling algorithm is to obtain a posteriori samples by sampling and constructing appropriate Markov chains to estimate the statistical characteristics of parameters. The main steps of standard MH sampling algorithm are as follows:(1)
Select the appropriate initialθ0 and make the likelihood function p(θ0|D) > 0, where D is the response matrix.(2)
Select the appropriate proposal distributionq(θ∗, θt). Use the current sample θt to generate a new candidate sample θ∗ obeying the proposal distribution.(3)
Calculate the acceptance probabilityα by the current sample θt and the candidate sample θ∗; the formula of acceptance probability is α = min[1, p(θ∗|D)/p (θt|D)].(4)
Generate random numberμ from the uniform distribution U(0,1). If α>μ, accept the candidate sample θ∗, θt+1 = θ∗; otherwise, reject the candidate sample θ∗, θt+1 = θt.(5)
Repeat steps (2) to (4) until the set sampling times are reached or the variance of samples is less than the set value. Terminate the iteration to generate a convergence sequence, the Markov chain.


## 2.1.1. Model Updating Method Based on Bayesian Theory

The FEMU method based on Bayesian theory combines a priori information (subjective information) and test data (objective information), and the MCMC algorithm is used to infer the posterior probability distribution of the updated parameters [17]. The procedure can be expressed by Bayesian formula as(1)pθ|x=px|θπθ∫θpx|θπθdθ=c⋅px|θπθ,where x is the observation information; π(θ) is the prior distribution of the uncertain parameter vector θ to be updated, which is usually equal to 1; c is the constant factor, called the regularization constant; and p(x|θ) is the conditional distribution under the given parameters, which is usually called the likelihood function. Assuming that N independent tests are carried out and the above likelihood function formula can be expressed as(2)px|θ=12πcovyNexp−Jθ,Jθ=∑i=1Nyi−yθiTcovy−1yi−yθi,where y is the test response vector; y(θ) is the calculated response vector of the model; covy is the covariance matrix of test response, substitute (2) into (1), and the posterior probability distribution of the parameter vector can be further written as(3)pθ|x=c′⋅exp−∑i=1Nyi−yθiTcovy−1yi−yθi,where c′ is a constant independent of θ. In engineering problems, the parameter vector is usually a complex high-dimensional vector and the response usually has no explicit expression. Therefore, the stationary Markov chains are obtained with sampling by using MCMC algorithm to estimate the posterior probability distribution of parameters, so as to avoid the problem of high-dimensional integral operation which may be involved in the above Bayesian formula.


## 2.1.2. Markov Chain Monte Carlo (MCMC) Algorithm

At present, in the field of Bayesian statistics, the MCMC algorithm has become an effective tool to deal with complex high-dimensional integral problems [18]. Firstly, the definition of Markov chain is given: there is a random sequence (θ(0), θ(1), …, θ(n)) at any time j(j ≥ 0), and the value θ(j+1) of the next time j + 1 in the sequence only depends on the current state of time j and has nothing to do with the historical state (θ(0), θ(1), …, θ(j−1)) before time j, which is called the Markov chain.The basic idea of MCMC algorithm is to obtain a Markov chain with stable distribution by sampling to get samples of the target probability distribution and then make statistical inference based on the valid samples in the convergence phase. In Bayesian model updating, the target probability distribution is the posterior probability distributionp(θ|x) of the parameter vector θ. When the samples of the target probability distribution are obtained by MCMC sampling, the statistical characteristics of the updated parameters can be estimated according to the large number theorem; for example, the mean and variance can be approximately calculated as(4)μθ=∫θpθ|xdθ≈1n∑j=1nθj,σθ2=∫θ−μθ2pθ|xdθ≈1n∑j=1nθj−μθ2,where n is the length of the effective Markov chain, that is, the number of valid samples. Markov chain is produced by sampling according to a certain approach in MCMC method, in which Metropolis–Hastings (MH) sampling is the most widely used MCMC algorithm, and many MCMC algorithms are the extension of standard MH sampling algorithm [19].


## 2.1.3. Standard Metropolis–Hastings (MH) Sampling Algorithm

The most prominent feature of the standard MH sampling algorithm is to obtain a posteriori samples by sampling and constructing appropriate Markov chains to estimate the statistical characteristics of parameters. The main steps of standard MH sampling algorithm are as follows:(1)
Select the appropriate initialθ0 and make the likelihood function p(θ0|D) > 0, where D is the response matrix.(2)
Select the appropriate proposal distributionq(θ∗, θt). Use the current sample θt to generate a new candidate sample θ∗ obeying the proposal distribution.(3)
Calculate the acceptance probabilityα by the current sample θt and the candidate sample θ∗; the formula of acceptance probability is α = min[1, p(θ∗|D)/p (θt|D)].(4)
Generate random numberμ from the uniform distribution U(0,1). If α>μ, accept the candidate sample θ∗, θt+1 = θ∗; otherwise, reject the candidate sample θ∗, θt+1 = θt.(5)
Repeat steps (2) to (4) until the set sampling times are reached or the variance of samples is less than the set value. Terminate the iteration to generate a convergence sequence, the Markov chain.


## 2.2. Feature Extraction of FRFs Based on SVD

### 2.2.1. Theory of Frequency Response Functions

For ann-degree-of-freedom damping system, the basic dynamic equation can be expressed as Mx¨+Cx˙+Kx=F, where M, C, and K are the mass, damping, and stiffness matrix, respectively; F is system excitation force vector; and x is response displacement vector. Under the action of simple harmonic excitation, the steady-state response X(ω) in frequency domain can be obtained by Fourier transform: X(ω) = H(ω)F(ω), where H(ω) is the FRFs; ω is the excitation frequency. The FRFs matrix can be expressed as(5)Hω=M−iCω−Kω2−1.


### 2.2.2. Theory of Singular Value Decomposition

Reconstruct the Hankel matrix from the original signalX = [x1, x2, …, xN] with noise, the Hankel matrix A can be expressed as(6)A=x1x2⋯xnx1×τ+1x1×τ+2⋯x1×τ+n⋮⋮⋮⋮xm−1τ+1xm−1τ+2⋯xm−1τ+n,where A is the Hankel matrix of order m × n; τ is the delay step and generally takes 1; m is the embedded dimension, when N is even, m = N/2; or m takes the median [20].First of all, the Hankel matrixA is decomposed by SVD, and then the appropriate threshold value λ is selected, and the first r larger singular values are selected to represent the signal characteristics according to the threshold, while the latter k-r singular values are set to 0. This procedure can be expressed as A=Ur∑rVrT, where Ur and VrT are orthogonal matrices; Σr = diag(σ1, σ2, …, σk) is a diagonal matrix; σi (i = 1,2, …, k) are singular values of matrix A, σ1 ≥ σ2 ≥ ··· ≥σk ≥ 0, k = min[(m − 1) × τ + 1, n].As a nonlinear filtering method, the SVD based on Hankel matrix is used to eliminate random noise in the signal. In this paper, the principle of sudden change in the number of extreme value points proposed in [15, 16] is used to test the rationality of the number of effective singular values, and the retained effective singular values are used to represent the FRFs to update the model.


## 2.2.1. Theory of Frequency Response Functions

For ann-degree-of-freedom damping system, the basic dynamic equation can be expressed as Mx¨+Cx˙+Kx=F, where M, C, and K are the mass, damping, and stiffness matrix, respectively; F is system excitation force vector; and x is response displacement vector. Under the action of simple harmonic excitation, the steady-state response X(ω) in frequency domain can be obtained by Fourier transform: X(ω) = H(ω)F(ω), where H(ω) is the FRFs; ω is the excitation frequency. The FRFs matrix can be expressed as(5)Hω=M−iCω−Kω2−1.


## 2.2.2. Theory of Singular Value Decomposition

Reconstruct the Hankel matrix from the original signalX = [x1, x2, …, xN] with noise, the Hankel matrix A can be expressed as(6)A=x1x2⋯xnx1×τ+1x1×τ+2⋯x1×τ+n⋮⋮⋮⋮xm−1τ+1xm−1τ+2⋯xm−1τ+n,where A is the Hankel matrix of order m × n; τ is the delay step and generally takes 1; m is the embedded dimension, when N is even, m = N/2; or m takes the median [20].First of all, the Hankel matrixA is decomposed by SVD, and then the appropriate threshold value λ is selected, and the first r larger singular values are selected to represent the signal characteristics according to the threshold, while the latter k-r singular values are set to 0. This procedure can be expressed as A=Ur∑rVrT, where Ur and VrT are orthogonal matrices; Σr = diag(σ1, σ2, …, σk) is a diagonal matrix; σi (i = 1,2, …, k) are singular values of matrix A, σ1 ≥ σ2 ≥ ··· ≥σk ≥ 0, k = min[(m − 1) × τ + 1, n].As a nonlinear filtering method, the SVD based on Hankel matrix is used to eliminate random noise in the signal. In this paper, the principle of sudden change in the number of extreme value points proposed in [15, 16] is used to test the rationality of the number of effective singular values, and the retained effective singular values are used to represent the FRFs to update the model.


## 2.3. Proposed MCMC Algorithm

### 2.3.1. Beetle Antennae Search (BAS) Algorithm

The BAS algorithm proposed by Jiang et al. [21] is an intelligent optimization algorithm inspired by the beetle’s foraging behavior. This algorithm is divided into two stages: search and detection.The search stage is based on the beetle’s food-searching behavior and can be expressed asxl = xi + di·b, xr = xi − di·b, where xl and xr represent the position coordinate of the left and right antennae of the beetle individual, respectively; di is the distance between the two antennae of the beetle individual during the ith search, also called the perceptual distance, which determines the search space of the algorithm. The perceptual distance gradually attenuates with the increase of search time, which makes the beetle individual approach to the optimal value adaptively; b is the unitized random direction vector.The detection stage updates the centroid position of the beetle individual, which can be expressed asxi = xi−1 + δi·b·sign(F(xl, X)−F(xr, X)), where F is the objective function, that is, the food odor concentration perceived by the left and right antennae of the beetle individual; X is the location coordinate; δi is the step length at the ith search, which is very important for the performance of the algorithm. An exponential attenuation model is introduced to constantly update the values of δ and d, and the attenuation model can be expressed as(7)di=d0⋅adi/Hd,δi=δ0⋅aδi/Hδ,where d0 and δ0 are the initial values of perceptual distance and search step, respectively; ad and aδ are attenuation coefficients; i is the number of search times; Hd and Hδ are attenuation constants.


### 2.3.2. Improved MCMC Algorithm Combining BAS Algorithm

In the procedure of Bayesian model updating, the phenomenon of “sampling stagnation” is easy to occur in standard MH sampling algorithm when the parameter dimension is high. At the same time, because the BAS algorithm does not need complex gradient calculation, the iterative speed is fast and the beetle individual can adaptively approach the optimal solution according to the objective function. In view of the above factors, the standard MH sampling algorithm can be improved by the BAS algorithm. From the expression of the likelihood function of (2) and the properties of the exponential function, the solution corresponding to the maximum value of the likelihood function is the solution corresponding to the maximum value of the exponential term [−J(θ)]. In this paper, the exponential term [−J(θ)] of likelihood function is taken as the objective function of BAS algorithm to guide the search of beetle individual. When the candidate sample point is rejected, the method of updating the centroid position of beetle individual is used to carry out local search on this point to generate a new candidate sample point. The specific sampling steps are as follows:Step 1: initialize the Markov chainsθ(t) (t = 1, 2, …, n), construct an appropriate proposal distribution q(θ∗, θt), where t is the number of iterations and n is the total number of iterations set by usersStep 2: assume the sample at the currentt iteration is θt, and a new candidate sample θ∗ is generated by standard MH samplingStep 3: calculate the acceptance probabilityα and generate random number μ from the uniform distribution U(0, 1). If α > μ, the candidate sample θ∗ is accepted, that is, θt+1 = θ∗; or θ∗ is not rejected immediately, and proceed to the next stepStep 4: fuse the idea of updating the centroid position of the beetle individual(a)
Initialize the parameters of BAS, and take the candidate sampleθ∗ as the initial centroid position of the beetle individual.(b)
Take the exponential term of the likelihood function expression in standard MH sampling algorithm as the search objective function to guide the iterative optimization of the beetle individual.(c)
Stop the search and calculate the acceptance rateαζ according to the centroid position θζ of the beetle individual when the set number of searches is reached (same as the standard MH algorithm to judge the acceptance status: if α > μ, θt+1 = θζ; otherwise, θt+1 = θt).Step 5: Repeat step 2 to step 4 until a convergence sequence is generated.Because the Gaussian distribution is widely used in engineering and its calculation is simple, the Gaussian distribution is selected as the proposal distribution in this paper.


## 2.3.1. Beetle Antennae Search (BAS) Algorithm

The BAS algorithm proposed by Jiang et al. [21] is an intelligent optimization algorithm inspired by the beetle’s foraging behavior. This algorithm is divided into two stages: search and detection.The search stage is based on the beetle’s food-searching behavior and can be expressed asxl = xi + di·b, xr = xi − di·b, where xl and xr represent the position coordinate of the left and right antennae of the beetle individual, respectively; di is the distance between the two antennae of the beetle individual during the ith search, also called the perceptual distance, which determines the search space of the algorithm. The perceptual distance gradually attenuates with the increase of search time, which makes the beetle individual approach to the optimal value adaptively; b is the unitized random direction vector.The detection stage updates the centroid position of the beetle individual, which can be expressed asxi = xi−1 + δi·b·sign(F(xl, X)−F(xr, X)), where F is the objective function, that is, the food odor concentration perceived by the left and right antennae of the beetle individual; X is the location coordinate; δi is the step length at the ith search, which is very important for the performance of the algorithm. An exponential attenuation model is introduced to constantly update the values of δ and d, and the attenuation model can be expressed as(7)di=d0⋅adi/Hd,δi=δ0⋅aδi/Hδ,where d0 and δ0 are the initial values of perceptual distance and search step, respectively; ad and aδ are attenuation coefficients; i is the number of search times; Hd and Hδ are attenuation constants.


## 2.3.2. Improved MCMC Algorithm Combining BAS Algorithm

In the procedure of Bayesian model updating, the phenomenon of “sampling stagnation” is easy to occur in standard MH sampling algorithm when the parameter dimension is high. At the same time, because the BAS algorithm does not need complex gradient calculation, the iterative speed is fast and the beetle individual can adaptively approach the optimal solution according to the objective function. In view of the above factors, the standard MH sampling algorithm can be improved by the BAS algorithm. From the expression of the likelihood function of (2) and the properties of the exponential function, the solution corresponding to the maximum value of the likelihood function is the solution corresponding to the maximum value of the exponential term [−J(θ)]. In this paper, the exponential term [−J(θ)] of likelihood function is taken as the objective function of BAS algorithm to guide the search of beetle individual. When the candidate sample point is rejected, the method of updating the centroid position of beetle individual is used to carry out local search on this point to generate a new candidate sample point. The specific sampling steps are as follows:Step 1: initialize the Markov chainsθ(t) (t = 1, 2, …, n), construct an appropriate proposal distribution q(θ∗, θt), where t is the number of iterations and n is the total number of iterations set by usersStep 2: assume the sample at the currentt iteration is θt, and a new candidate sample θ∗ is generated by standard MH samplingStep 3: calculate the acceptance probabilityα and generate random number μ from the uniform distribution U(0, 1). If α > μ, the candidate sample θ∗ is accepted, that is, θt+1 = θ∗; or θ∗ is not rejected immediately, and proceed to the next stepStep 4: fuse the idea of updating the centroid position of the beetle individual(a)
Initialize the parameters of BAS, and take the candidate sampleθ∗ as the initial centroid position of the beetle individual.(b)
Take the exponential term of the likelihood function expression in standard MH sampling algorithm as the search objective function to guide the iterative optimization of the beetle individual.(c)
Stop the search and calculate the acceptance rateαζ according to the centroid position θζ of the beetle individual when the set number of searches is reached (same as the standard MH algorithm to judge the acceptance status: if α > μ, θt+1 = θζ; otherwise, θt+1 = θt).Step 5: Repeat step 2 to step 4 until a convergence sequence is generated.Because the Gaussian distribution is widely used in engineering and its calculation is simple, the Gaussian distribution is selected as the proposal distribution in this paper.


## 2.4. Construct Support Vector Machine (SVM) Surrogate Model

### 2.4.1. Latin Hypercube Sampling (LHS)

In order to ensure the randomness and uniformity of the training set and test set samples, the LHS is introduced to obtain the samples to construct the surrogate model. LHS is a multidimensional stratified random sampling method, and the sampling points have good distribution uniformity and representativeness [22]. It is assumed that q-group samples need to be taken within the prior range of m-dimensional model parameters. The specific sampling steps are as follows:(1)
Them prior ranges of m-dimensional model parameters [Ai, Bi] (i = 1, 2, …, m) are divided into q cells, which can be recorded as [Aij, Bij] (i = 1, 2, …, m) (j = 1, 2, …, q), resulting in a total of m × q cells.(2)
A sample point is randomly selected in any intercell [Aij, Bij] and marked as αij sampling, and a total of m × q sample points are generated to form a matrix Γ:(8)Γ=α11α12⋯α1qα21α22⋯α2q⋮⋮⋮⋮αm1αm2⋯αmq.(3)
By randomly sorting each row vector [αi1, αi2,·…, αiq] (i = 1, 2, …, m) in the matrix, the resulting vector is [βi1, βi2, …, βiq] (i = 1, 2, …, m), and the matrix Ф is obtained:(9)Φ=β11β12⋯β1qβ21β22⋯β2q⋮⋮⋮⋮βm1βm2⋯βmq.Each column vector in the matrixФ is a group of samples, and a total of q groups of samples are obtained.


### 2.4.2. SVM Surrogate Model

The surrogate model fits the complex relationship between the FE model parameters and the response by an explicit function. The surrogate model can replace the complex FE model for iterative calculation so as to greatly improve the FEMU efficiency [23]. Clarke et al. [24] compared the performance of SVM with Radial Basis Function (RBF), Neural Network (NN), Kriging model, and the Response Surface surrogate model by engineering examples. The results show that the robustness and accuracy of SVM surrogate model is better than the other three models. So, the SVM is selected as the surrogate model in this paper. When fitting the SVM, the LHS is used to obtain the training set samples and test set samples, and then the corresponding singular values of FRFs are obtained. These data are normalized to avoid error caused by large orders of magnitude.The SVM surrogate model [25] can be expressed as f(s) = ω·ϕ(s) + b, where ϕ(·) is a nonlinear function, which maps the data S to a high-dimensional linear space; ω is a weight vector; b is a deviation.The optimization problem can be expressed as(10)minw,b,ξ12w2+C∑i=1lξi+ξi∗,s.t.yi−ω⋅ϕsi−b≤ε+ξi∗,ω⋅ϕsi+b−yi≤ε+ξi∗,ξi≥0,ξi∗≥0,i=1,…,l,where C is the penalty coefficient; ξi and ξi∗∗ are the relaxation variables, which ensure that the constraint conditions can be satisfied. The Lagrange multiplier method is introduced to transform the above quadratic programming problem to a dual problem, and the regression estimation function can be expressed as(11)fs=∑i=1nαi−αi∗Ksi,si∗+b,where αi and αi∗ are the Lagrange multipliers; K(si, si∗) is the kernel function. Here, the widely used RBF is selected as the kernel function.


## 2.4.1. Latin Hypercube Sampling (LHS)

In order to ensure the randomness and uniformity of the training set and test set samples, the LHS is introduced to obtain the samples to construct the surrogate model. LHS is a multidimensional stratified random sampling method, and the sampling points have good distribution uniformity and representativeness [22]. It is assumed that q-group samples need to be taken within the prior range of m-dimensional model parameters. The specific sampling steps are as follows:(1)
Them prior ranges of m-dimensional model parameters [Ai, Bi] (i = 1, 2, …, m) are divided into q cells, which can be recorded as [Aij, Bij] (i = 1, 2, …, m) (j = 1, 2, …, q), resulting in a total of m × q cells.(2)
A sample point is randomly selected in any intercell [Aij, Bij] and marked as αij sampling, and a total of m × q sample points are generated to form a matrix Γ:(8)Γ=α11α12⋯α1qα21α22⋯α2q⋮⋮⋮⋮αm1αm2⋯αmq.(3)
By randomly sorting each row vector [αi1, αi2,·…, αiq] (i = 1, 2, …, m) in the matrix, the resulting vector is [βi1, βi2, …, βiq] (i = 1, 2, …, m), and the matrix Ф is obtained:(9)Φ=β11β12⋯β1qβ21β22⋯β2q⋮⋮⋮⋮βm1βm2⋯βmq.Each column vector in the matrixФ is a group of samples, and a total of q groups of samples are obtained.


## 2.4.2. SVM Surrogate Model

The surrogate model fits the complex relationship between the FE model parameters and the response by an explicit function. The surrogate model can replace the complex FE model for iterative calculation so as to greatly improve the FEMU efficiency [23]. Clarke et al. [24] compared the performance of SVM with Radial Basis Function (RBF), Neural Network (NN), Kriging model, and the Response Surface surrogate model by engineering examples. The results show that the robustness and accuracy of SVM surrogate model is better than the other three models. So, the SVM is selected as the surrogate model in this paper. When fitting the SVM, the LHS is used to obtain the training set samples and test set samples, and then the corresponding singular values of FRFs are obtained. These data are normalized to avoid error caused by large orders of magnitude.The SVM surrogate model [25] can be expressed as f(s) = ω·ϕ(s) + b, where ϕ(·) is a nonlinear function, which maps the data S to a high-dimensional linear space; ω is a weight vector; b is a deviation.The optimization problem can be expressed as(10)minw,b,ξ12w2+C∑i=1lξi+ξi∗,s.t.yi−ω⋅ϕsi−b≤ε+ξi∗,ω⋅ϕsi+b−yi≤ε+ξi∗,ξi≥0,ξi∗≥0,i=1,…,l,where C is the penalty coefficient; ξi and ξi∗∗ are the relaxation variables, which ensure that the constraint conditions can be satisfied. The Lagrange multiplier method is introduced to transform the above quadratic programming problem to a dual problem, and the regression estimation function can be expressed as(11)fs=∑i=1nαi−αi∗Ksi,si∗+b,where αi and αi∗ are the Lagrange multipliers; K(si, si∗) is the kernel function. Here, the widely used RBF is selected as the kernel function.


## 2.5. Specific Procedure of Model Updating

The training set samples and test set samples are generated by the LHS. The modal participation variation coefficient criterion [26] is introduced to select the appropriate location of excitation points. The modal kinetic energy method [27] is introduced to select the appropriate locations of measuring points. The calculated FRFs are decomposed by SVD to extract the effective singular values according to the principle of sudden change in the number of extreme value points. According to the input parameters and the corresponding singular values, the SVM surrogate model is constructed and its accuracy is tested. Then the SVM surrogate model which meets the accuracy requirement is introduced into the proposed algorithm to iteratively solve the model updating problem, and the mean values of parameters are estimated according to the statistical characteristics of the convergent Markov chains. The specific model updating procedure is shown in Figure 1.Figure 1
Flow chart of FEMU through the proposed MCMC algorithm.


## 3. Numerical Examples

### 3.1. Example 1

As shown in Figure2, the plane truss structure consists of 14 nodes, 25 rod elements, and 25 degrees of freedom. The constraint condition is the fixed hinge connection of support 1 (node 1) and the movable hinge connection of support 2 (node 13). The rod elements are divided into two groups, numbered as1, 2∼7, 8 and7, 8∼13, 14, respectively. The elasticity modulus and density of two groups of rod elements are divided by their initial values as the parameters to be updated, which are θ1 = E1/E0, θ2 = ρ1/ρ0, θ3 = E2/E0, and θ4 = ρ2/ρ0, respectively, where the initial elasticity modulus E0 is 210 GPa and the initial density ρ0 is 7850 kg·m−3. Suppose that the test mean value of each parameter to be updated obeys the normal distribution matrix [1.5, 2, 0.8, 3]; the standard deviations are all 0.1; and the ranges of parameters are set to [1, 2], [1.5, 2.5], [0.3, 1,3], and [2.5, 3.5], respectively.Figure 2
Plane truss model.3000 groups of samples are selected by the LHS. The first 2500 groups are used as training set while the last 500 groups are as test set. The degree of freedom is analyzed by modal variation coefficient criterion. The 25th, the 22nd, and the 6th degrees of freedom are selected as excitation points. The truss nodes in the diagram are numbered as 14-nodey-direction, 12-node y-direction, and 4-node y direction. The corresponding 4th, 6th, and 18th degrees of freedom are selected as the measurement points by introducing the modal kinetic energy method, and the truss nodes in the corresponding diagram are numbered as 3-node y-direction, 4-node y-direction, and 10-node y-direction. The FRFs are calculated and SVD is carried out. According to the principle of sudden change in the number of extreme value points, the first three large singular values are extracted as the characteristics of the FRFs, and then the training set and corresponding singular values are used to construct the SVM surrogate model. The error values of the first three singular values of the first FRF are selected to test the accuracy of the SVM, as listed in Table 1.Table 1
Error of the SVM surrogate model.
OrderError (%)MinimumMaximumAverage16.70 × 10−21.42 × 10−11.27 × 10−121.30 × 10−37.67 × 10−21.43 × 10−231.10 × 10−41.64 × 10−14.20 × 10−3In order to consider the uncertainty in the structure, 50 groups of samples are randomly generated to simulate experiment, and the samples are substituted into the model to calculate the first three large singular values of the three FRFs. The standard MH sampling algorithm and the proposed algorithm sample are 6000 times, respectively. The cloud map of the sample points distribution containing 95% confidence ellipse is illustrated in Figure3. As can be seen, most of the sample points extracted by the proposed algorithm can be located within 95% of the high probability areas.Figure 3
Updated characteristic quantity and confidence ellipse.
(a)(b)The Markov chains obtained by two methods are shown in Figure4. Both of them can converge to the preset mean values of parameters. But there are many smooth segments in the Markov chains obtained by standard MH sampling algorithm, indicating the phenomenon of “sampling stagnation” in the sampling procedure. While the Markov chains obtained by the proposed algorithm can converge quickly and keep fluctuating all the time, the ergodicity of samples is obviously better than that of the former. Compared with the acceptance rate of 21% of standard MH sampling algorithm, the acceptance rate of the proposed algorithm is about 62.3%, so the acceptance rate is significantly improved.Figure 4
Comparison of four Markov chains obtained by two algorithms. (a) Standard MH sampling algorithm. (b) Proposed algorithm.
(a)(b)The posterior samples of parameters in Figure4(b) are analyzed according to the knowledge of normal probability test. Considering that the samples during the nonstationary period will interfere with the statistical characteristics of the parameters, the first 10% of the samples are removed as the combustion stage, and then the remaining samples are tested for normal distribution, and the test results are shown in Figure 5. As can be seen, most of the samples obtained by the proposed algorithm can fall on the hypothetical normal distribution line, except that some of the combustion period samples still exist and deviate sparsely at the boundary of the probability graph, so the posterior samples conform to the normal distribution hypothesis.Figure 5
Normal test probability map of four parameters.
(a)(b)(c)(d)The mean values of updated parameters are estimated indirectly from the statistical characteristics of the posterior samples. The updated results are obtained as listed in Table2. The relative errors of parameters after updating are less than 0.3%, indicating that the proposed algorithm has a good effect for the uncertainty FEMU.Table 2
Comparison of updating results of the plane truss model.
ParameterTest mean valueInitial mean valuePreupdating error (%)Updated mean valueUpdated error (%)θ11.5001.95330.1871.4980.107θ22.0002.39919.9301.9960.211θ30.8001.06933.5880.7990.106θ43.0002.7747.5372.9990.009In order to further verify the updating effect of the proposed algorithm, the acceleration FRFs corresponding to the test mean model, the initial finite element mean model, and the updated mean model are calculated according to the data in Table2, and the acceleration FRFs curves are shown in Figure 6. The updated FRF curve almost coincides with the FRF curve of test mean model.Figure 6
Comparison of FRF curve before and after updating.In order to obtain the updating effect more clearly, the real and imaginary part curves of the acceleration FRFs of the mean model before and after updating are compared, and the results are shown in Figures7 and 8. Compared with the real and imaginary part curves of the preupdating mean model FRF, the real and imaginary part curves of the updated mean model FRF are closer to and basically coincide with the real and imaginary part curves of the test mean model FRF.Figure 7
Comparison of real parts of FRF before and after updating.Figure 8
Comparison of imaginary parts of FRF before and after updating.To sum up, it is shown that, compared with the standard MH sampling algorithm, the proposed algorithm can greatly improve the acceptance rate of samples and make the overall performance of the Markov chains better while ensuring the better updating accuracy.


### 3.2. Example 2

In order to further verify the reliability of the proposed algorithm, single damage and multiple damage are added to the structure and then the proposed algorithm is used to identify the parameters. A three-dimensional truss structure shown in Figure9 is selected for verification.Figure 9
Three-dimensional truss model.The structure consists of 28 nodes, 66 rod elements, and 48 degrees of freedom. The nodes are hinged. The constraint conditions are fixed by 4 supports (node numbers 1, 8, 9, and 16). For each node, only translational degrees of freedom inY and Z directions are considered. The selected three excitation points are the 4th, 8th, and 10th degrees of freedom, and the corresponding measuring points are the 38th, 20th, and 24th degrees of freedom. The corresponding truss nodes are shown in Figure 9 and are all in Y direction.The rod elements are divided into three parts: the lower rod elements, the middle rod elements (straight rods and oblique rods), and the upper rod elements. Since the damage will cause significant change in the structural stiffness in most engineering cases, three parametersθ1 = E1/E0, θ2 = E2/E0, and θ3 = E3/E0 related to the elasticity modulus of three parts are selected as the parameters to be updated, and the initial elasticity modulus E0 is 190 GPa. It is assumed that all the test mean values of parameters obey the normal mean matrix [1.5, 1, 2], and the standard deviations are all 0.1. The ranges of parameters to be updated are set to [1, 2], [0.5, 1.5], and [1.5, 2.5], respectively. The damage of the rod elements is listed in Table 3.Table 3
Damage state of truss model.
Damage stateSpecific damage situationS0No damageS1Elasticity modulus of the lower rod elements is damaged by 20%S2Elasticity modulus of the lower rod elements is damaged by 20%; elasticity modulus of the middle rod elements is damaged by 30%; elasticity modulus of the upper rod elements is damaged by 25%The parameters are identified in three damage states. Markov chains obtained by the standard MH sampling algorithm and the proposed algorithm in the S0 state are shown in Figure 10. Through the comparison of Figures 10(a) and 10(b), it is obvious that there are many smooth segments in the Markov chains obtained by the standard MH sampling algorithm, that is, the phenomenon of "sampling stagnation", while the Markov chains obtained by the proposed algorithm can always remain in a fluctuating state, and the overall performance is better.Figure 10
Markov chains obtained by different algorithms in the S0 damage state. (a) Standard MH sampling. (b) Proposed algorithm.
(a)(b)Normal distribution test is carried out on the posterior samples of Markov chains obtained by the proposed algorithm. After excluding the first 10% nonstationary period samples, the stationary period samples (samples in the convergence phase of Markov chain) are tested, and the results are shown in Figure11. Most of the samples obtained by the proposed algorithm can fall on the hypothetical normal distribution line, so the normal distribution hypothesis of posterior samples is valid. Then the mean values of the updated parameters are estimated by the Markov chains, and the updated results are listed in Table 4. The updating errors of the proposed algorithm in the S0 state are basically within 1.5%.Figure 11
Normal probability test diagram of parameters in the S0 state.
(a)(b)(c)Table 4
Comparison of updating results in the S0 damage state.
ParameterTest mean valueInitial mean valuePreupdating error (%)Updated mean valueUpdated error (%)θ11.501.86224.151.5130.90θ21.001.31931.890.9881.18θ32.002.21510.751.9890.53The Markov chains and the posterior distribution histograms obtained by the proposed algorithm in the S1 and S2 damage state are shown in Figures 12 and 13, respectively. The updating results of the proposed algorithm in the S1 and S2 damage state are listed in Tables 5 and 6, respectively. The Markov chains obtained in the two damage states can quickly converge to the preset mean values. Compared with the sample acceptance rate of about 21% of the standard MH sampling algorithm, the sample acceptance rate of the proposed algorithm can be increased to about 67%, and the overall performance of Markov chains is improved. This reveals that the proposed algorithm can identify the changes of parameters in the case of structural damage.Figure 12
Markov chains and histograms of posterior distribution in the S1 state.
(a)(b)(c)Figure 13
Markov chains and histograms of posterior distribution in the S2 state.
(a)(b)(c)Table 5
Comparison of updating results in the S1 damage state.
ParameterTest mean valueInitial mean valuePreupdating error (%)Updated mean valueUpdated error (%)θ11.2001.562230.181.21150.96θ21.0001.318931.890.98921.08θ32.0002.215010.751.99290.35Table 6
Comparison of updating results in the S2 damage state.
ParameterTest mean valueInitial mean valuePreupdating error (%)Updated mean valueUpdated error (%)θ11.2001.563130.261.18850.95θ20.7001.018945.560.71432.04θ31.5001.715014.331.49270.48In order to verify the antinoise performance of the proposed algorithm, the white Gaussian noises with different signal-to-noise ratio (SNR) are added to the test FRFs in the S0 state, and the SNR are 50 dB, 30 dB, 20 dB, 10 dB, and 5 dB, respectively. For the length limitation, only FRF curve under 10 dB SNR is shown in Figure 14.Figure 14
FRF curve under 10 dB SNR.For the FRFs disturbed by noise, it is necessary to retain the larger singular values according to the threshold, so the selection of threshold is particularly important. According to the principle of sudden change in the number of extreme value points mentioned above, a reasonable threshold is selected. The sudden change point of the extreme points of singular values under 10 dB SNR is shown in Figure15. It is worth noting that the better characterization of the FRFs is obtained with the increase of the number of singular values within the sudden change point of the extreme points, but at the same time, it will also increase the computational cost of model updating. The number of extreme points of singular values is 13, so the first 3 larger singular values selected can meet the requirement of retaining 13 singular values at most.Figure 15
Sudden change point of the number of extreme value points under 10 dB SNR.Under each SNR noise, the average values after 10 running times are taken as the updating results listed in Table7. The proposed algorithm can still obtain satisfactory updating results under different SNR noises. It shows that the proposed algorithm has good antinoise performance.Table 7
Comparison of updating results after adding noise.
ParameterSNR (50 dB)SNR (30 dB)SNR (20 dB)SNR (10 dB)SNR (5 dB)Updated valueUpdated error (%)Updated valueUpdated error (%)Updated valueUpdated error (%)Updated valueUpdated error (%)Updated valueUpdated error (%)θ11.4900.6751.4880.8091.4801.3281.5271.8251.5473.151θ20.9950.4650.9990.0580.9811.8670.9653.5140.9623.839θ31.9950.2362.0170.8431.9741.2951.8975.1282.1256.244


## 3.1. Example 1

As shown in Figure2, the plane truss structure consists of 14 nodes, 25 rod elements, and 25 degrees of freedom. The constraint condition is the fixed hinge connection of support 1 (node 1) and the movable hinge connection of support 2 (node 13). The rod elements are divided into two groups, numbered as1, 2∼7, 8 and7, 8∼13, 14, respectively. The elasticity modulus and density of two groups of rod elements are divided by their initial values as the parameters to be updated, which are θ1 = E1/E0, θ2 = ρ1/ρ0, θ3 = E2/E0, and θ4 = ρ2/ρ0, respectively, where the initial elasticity modulus E0 is 210 GPa and the initial density ρ0 is 7850 kg·m−3. Suppose that the test mean value of each parameter to be updated obeys the normal distribution matrix [1.5, 2, 0.8, 3]; the standard deviations are all 0.1; and the ranges of parameters are set to [1, 2], [1.5, 2.5], [0.3, 1,3], and [2.5, 3.5], respectively.Figure 2
Plane truss model.3000 groups of samples are selected by the LHS. The first 2500 groups are used as training set while the last 500 groups are as test set. The degree of freedom is analyzed by modal variation coefficient criterion. The 25th, the 22nd, and the 6th degrees of freedom are selected as excitation points. The truss nodes in the diagram are numbered as 14-nodey-direction, 12-node y-direction, and 4-node y direction. The corresponding 4th, 6th, and 18th degrees of freedom are selected as the measurement points by introducing the modal kinetic energy method, and the truss nodes in the corresponding diagram are numbered as 3-node y-direction, 4-node y-direction, and 10-node y-direction. The FRFs are calculated and SVD is carried out. According to the principle of sudden change in the number of extreme value points, the first three large singular values are extracted as the characteristics of the FRFs, and then the training set and corresponding singular values are used to construct the SVM surrogate model. The error values of the first three singular values of the first FRF are selected to test the accuracy of the SVM, as listed in Table 1.Table 1
Error of the SVM surrogate model.
OrderError (%)MinimumMaximumAverage16.70 × 10−21.42 × 10−11.27 × 10−121.30 × 10−37.67 × 10−21.43 × 10−231.10 × 10−41.64 × 10−14.20 × 10−3In order to consider the uncertainty in the structure, 50 groups of samples are randomly generated to simulate experiment, and the samples are substituted into the model to calculate the first three large singular values of the three FRFs. The standard MH sampling algorithm and the proposed algorithm sample are 6000 times, respectively. The cloud map of the sample points distribution containing 95% confidence ellipse is illustrated in Figure3. As can be seen, most of the sample points extracted by the proposed algorithm can be located within 95% of the high probability areas.Figure 3
Updated characteristic quantity and confidence ellipse.
(a)(b)The Markov chains obtained by two methods are shown in Figure4. Both of them can converge to the preset mean values of parameters. But there are many smooth segments in the Markov chains obtained by standard MH sampling algorithm, indicating the phenomenon of “sampling stagnation” in the sampling procedure. While the Markov chains obtained by the proposed algorithm can converge quickly and keep fluctuating all the time, the ergodicity of samples is obviously better than that of the former. Compared with the acceptance rate of 21% of standard MH sampling algorithm, the acceptance rate of the proposed algorithm is about 62.3%, so the acceptance rate is significantly improved.Figure 4
Comparison of four Markov chains obtained by two algorithms. (a) Standard MH sampling algorithm. (b) Proposed algorithm.
(a)(b)The posterior samples of parameters in Figure4(b) are analyzed according to the knowledge of normal probability test. Considering that the samples during the nonstationary period will interfere with the statistical characteristics of the parameters, the first 10% of the samples are removed as the combustion stage, and then the remaining samples are tested for normal distribution, and the test results are shown in Figure 5. As can be seen, most of the samples obtained by the proposed algorithm can fall on the hypothetical normal distribution line, except that some of the combustion period samples still exist and deviate sparsely at the boundary of the probability graph, so the posterior samples conform to the normal distribution hypothesis.Figure 5
Normal test probability map of four parameters.
(a)(b)(c)(d)The mean values of updated parameters are estimated indirectly from the statistical characteristics of the posterior samples. The updated results are obtained as listed in Table2. The relative errors of parameters after updating are less than 0.3%, indicating that the proposed algorithm has a good effect for the uncertainty FEMU.Table 2
Comparison of updating results of the plane truss model.
ParameterTest mean valueInitial mean valuePreupdating error (%)Updated mean valueUpdated error (%)θ11.5001.95330.1871.4980.107θ22.0002.39919.9301.9960.211θ30.8001.06933.5880.7990.106θ43.0002.7747.5372.9990.009In order to further verify the updating effect of the proposed algorithm, the acceleration FRFs corresponding to the test mean model, the initial finite element mean model, and the updated mean model are calculated according to the data in Table2, and the acceleration FRFs curves are shown in Figure 6. The updated FRF curve almost coincides with the FRF curve of test mean model.Figure 6
Comparison of FRF curve before and after updating.In order to obtain the updating effect more clearly, the real and imaginary part curves of the acceleration FRFs of the mean model before and after updating are compared, and the results are shown in Figures7 and 8. Compared with the real and imaginary part curves of the preupdating mean model FRF, the real and imaginary part curves of the updated mean model FRF are closer to and basically coincide with the real and imaginary part curves of the test mean model FRF.Figure 7
Comparison of real parts of FRF before and after updating.Figure 8
Comparison of imaginary parts of FRF before and after updating.To sum up, it is shown that, compared with the standard MH sampling algorithm, the proposed algorithm can greatly improve the acceptance rate of samples and make the overall performance of the Markov chains better while ensuring the better updating accuracy.


## 3.2. Example 2

In order to further verify the reliability of the proposed algorithm, single damage and multiple damage are added to the structure and then the proposed algorithm is used to identify the parameters. A three-dimensional truss structure shown in Figure9 is selected for verification.Figure 9
Three-dimensional truss model.The structure consists of 28 nodes, 66 rod elements, and 48 degrees of freedom. The nodes are hinged. The constraint conditions are fixed by 4 supports (node numbers 1, 8, 9, and 16). For each node, only translational degrees of freedom inY and Z directions are considered. The selected three excitation points are the 4th, 8th, and 10th degrees of freedom, and the corresponding measuring points are the 38th, 20th, and 24th degrees of freedom. The corresponding truss nodes are shown in Figure 9 and are all in Y direction.The rod elements are divided into three parts: the lower rod elements, the middle rod elements (straight rods and oblique rods), and the upper rod elements. Since the damage will cause significant change in the structural stiffness in most engineering cases, three parametersθ1 = E1/E0, θ2 = E2/E0, and θ3 = E3/E0 related to the elasticity modulus of three parts are selected as the parameters to be updated, and the initial elasticity modulus E0 is 190 GPa. It is assumed that all the test mean values of parameters obey the normal mean matrix [1.5, 1, 2], and the standard deviations are all 0.1. The ranges of parameters to be updated are set to [1, 2], [0.5, 1.5], and [1.5, 2.5], respectively. The damage of the rod elements is listed in Table 3.Table 3
Damage state of truss model.
Damage stateSpecific damage situationS0No damageS1Elasticity modulus of the lower rod elements is damaged by 20%S2Elasticity modulus of the lower rod elements is damaged by 20%; elasticity modulus of the middle rod elements is damaged by 30%; elasticity modulus of the upper rod elements is damaged by 25%The parameters are identified in three damage states. Markov chains obtained by the standard MH sampling algorithm and the proposed algorithm in the S0 state are shown in Figure 10. Through the comparison of Figures 10(a) and 10(b), it is obvious that there are many smooth segments in the Markov chains obtained by the standard MH sampling algorithm, that is, the phenomenon of "sampling stagnation", while the Markov chains obtained by the proposed algorithm can always remain in a fluctuating state, and the overall performance is better.Figure 10
Markov chains obtained by different algorithms in the S0 damage state. (a) Standard MH sampling. (b) Proposed algorithm.
(a)(b)Normal distribution test is carried out on the posterior samples of Markov chains obtained by the proposed algorithm. After excluding the first 10% nonstationary period samples, the stationary period samples (samples in the convergence phase of Markov chain) are tested, and the results are shown in Figure11. Most of the samples obtained by the proposed algorithm can fall on the hypothetical normal distribution line, so the normal distribution hypothesis of posterior samples is valid. Then the mean values of the updated parameters are estimated by the Markov chains, and the updated results are listed in Table 4. The updating errors of the proposed algorithm in the S0 state are basically within 1.5%.Figure 11
Normal probability test diagram of parameters in the S0 state.
(a)(b)(c)Table 4
Comparison of updating results in the S0 damage state.
ParameterTest mean valueInitial mean valuePreupdating error (%)Updated mean valueUpdated error (%)θ11.501.86224.151.5130.90θ21.001.31931.890.9881.18θ32.002.21510.751.9890.53The Markov chains and the posterior distribution histograms obtained by the proposed algorithm in the S1 and S2 damage state are shown in Figures 12 and 13, respectively. The updating results of the proposed algorithm in the S1 and S2 damage state are listed in Tables 5 and 6, respectively. The Markov chains obtained in the two damage states can quickly converge to the preset mean values. Compared with the sample acceptance rate of about 21% of the standard MH sampling algorithm, the sample acceptance rate of the proposed algorithm can be increased to about 67%, and the overall performance of Markov chains is improved. This reveals that the proposed algorithm can identify the changes of parameters in the case of structural damage.Figure 12
Markov chains and histograms of posterior distribution in the S1 state.
(a)(b)(c)Figure 13
Markov chains and histograms of posterior distribution in the S2 state.
(a)(b)(c)Table 5
Comparison of updating results in the S1 damage state.
ParameterTest mean valueInitial mean valuePreupdating error (%)Updated mean valueUpdated error (%)θ11.2001.562230.181.21150.96θ21.0001.318931.890.98921.08θ32.0002.215010.751.99290.35Table 6
Comparison of updating results in the S2 damage state.
ParameterTest mean valueInitial mean valuePreupdating error (%)Updated mean valueUpdated error (%)θ11.2001.563130.261.18850.95θ20.7001.018945.560.71432.04θ31.5001.715014.331.49270.48In order to verify the antinoise performance of the proposed algorithm, the white Gaussian noises with different signal-to-noise ratio (SNR) are added to the test FRFs in the S0 state, and the SNR are 50 dB, 30 dB, 20 dB, 10 dB, and 5 dB, respectively. For the length limitation, only FRF curve under 10 dB SNR is shown in Figure 14.Figure 14
FRF curve under 10 dB SNR.For the FRFs disturbed by noise, it is necessary to retain the larger singular values according to the threshold, so the selection of threshold is particularly important. According to the principle of sudden change in the number of extreme value points mentioned above, a reasonable threshold is selected. The sudden change point of the extreme points of singular values under 10 dB SNR is shown in Figure15. It is worth noting that the better characterization of the FRFs is obtained with the increase of the number of singular values within the sudden change point of the extreme points, but at the same time, it will also increase the computational cost of model updating. The number of extreme points of singular values is 13, so the first 3 larger singular values selected can meet the requirement of retaining 13 singular values at most.Figure 15
Sudden change point of the number of extreme value points under 10 dB SNR.Under each SNR noise, the average values after 10 running times are taken as the updating results listed in Table7. The proposed algorithm can still obtain satisfactory updating results under different SNR noises. It shows that the proposed algorithm has good antinoise performance.Table 7
Comparison of updating results after adding noise.
ParameterSNR (50 dB)SNR (30 dB)SNR (20 dB)SNR (10 dB)SNR (5 dB)Updated valueUpdated error (%)Updated valueUpdated error (%)Updated valueUpdated error (%)Updated valueUpdated error (%)Updated valueUpdated error (%)θ11.4900.6751.4880.8091.4801.3281.5271.8251.5473.151θ20.9950.4650.9990.0580.9811.8670.9653.5140.9623.839θ31.9950.2362.0170.8431.9741.2951.8975.1282.1256.244


## 4. Conclusion

In this paper, the method of updating the centroid position of beetle individual in BAS algorithm is introduced into the standard MH sampling algorithm, and a novel uncertain FEMU method based on improved MCMC algorithm is proposed. It is verified by a plane truss structure and three-dimensional truss structure with single and multidamage conditions, and the following conclusions are obtained.(1)
The uncertainty FEMU method proposed in this paper can not only ensure the better updating accuracy, but also greatly improve the sample acceptance rate of the Markov chains and make the overall performance of the Markov chains better.(2)
Using the effective singular values as the response to construct the SVM surrogate model can not only improve the efficiency of model updating, but also indirectly introduce the FRFs into the Bayesian model updating, avoiding the complexity of the derivation of a posteriori probability density function and the difficulty of frequency point selection in traditional methods.(3)
Using the effective singular values to represent the FRFs to update the model can effectively reduce the interference of noise to the model updating. Even if white Gaussian noise is added to the test signal, the proposed method can still obtain satisfactory updating results, which proves that the proposed method is robust to noise.



---

*Source: 5543549-2021-04-30.xml*