A Kriging Model Based Finite Element Model Updating Method for Damage Detection

Model updating is an effective means of damage identification and surrogate modeling has attracted considerable attention for saving computational cost in finite element (FE) model updating, especially for large-scale structures. In this context, a surrogate model of frequency is normally constructed for damage identification, while the frequency response function (FRF) is rarely used as it usually changes dramatically with updating parameters. This paper presents a new surrogate model based model updating method taking advantage of the measured FRFs. The Frequency Domain Assurance Criterion (FDAC) is used to build the objective function, whose nonlinear response surface is constructed by the Kriging model. Then, the efficient global optimization (EGO) algorithm is introduced to get the model updating results. The proposed method has good accuracy and robustness, which have been verified by a numerical simulation of a cantilever and experimental test data of a laboratory three-story structure.


Introduction
Damage can occur in civil, mechanical and aerospace structures after long service or under extreme environmental events.Once accumulated to a certain extent, damage can lead to structure failure, resulting in significant loss of life and property [1].Therefore, how to use effective damage identification methods to monitor and evaluate damage has been a popular topic in engineering.
As a global damage detection methodology, finite element (FE) model updating has been widely studied and applied in recent decades [2][3][4].The basic idea of it is that an FE model and the actual structure often have some deviations, and the reliability of the FE model can be improved by modifying the physical parameters of the FE model and making the response predicted by the modified FE model as close as possible to the actual structural response.The methodology was originally used as a major means of correcting a theoretical model before structural analysis and optimization, to ensure accuracy and validity of the constructed FE model.Its application in the field of damage identification is mainly to modify the parameters of the FE models through the measured response data of the initial and damaged structures, respectively.Then, the position and extent of the damage are assessed by comparing the difference between the parameters of the elements or substructures of the two updated FE models.It is a method to locate and quantify damage simultaneously.
In the early research, the sensitivity based model updating method has achieved great success [5].The principle of this method is discussed in detail by Mottershead, et al. [6] and the three most common dynamic responses used in this method are frequency, mode shape and frequency response function.Among them, the frequency and mode shape are widely studied to solve engineering problems, where usually only some degrees of freedom are measured [7].The mode shape is mainly used as a supplement to the frequency to deal with the case of more updating parameters than with the number of measured frequencies [8,9].The frequency response function itself can provide enough response information and has been used in many studies of model updating [10][11][12][13][14][15][16][17][18][19][20][21], but it also places high demand on the testing environment, which, to some extent, limits its application.
One of the main limitations of model updating methods in engineering applications is the computational efficiency.In theory, the denser the finite element mesh is, the more accurate the prediction of the structural response is-a dense FE mesh improves the reliability of the model updating results.However, a very fine finite element mesh will result in a considerable increase in the computational time of a single theoretical model response.If the model updating procedure requires a considerable number of iterations to converge, the required computation time will be long and even become unacceptable.
In this context, the surrogate model (also called the response surface method), as a quick alternative to computer experiments (FE analysis), has attracted the attention of many researchers.The core idea of the surrogate model is to replace the FE analysis model by simulating a much simplified model for the relationship between the parameters and the response.After confirming the accuracy of the surrogate model, the response of the theoretical model is calculated directly through this model in the iterative process to achieve the purpose of quick solution.
At present, the surrogate model based global optimization has been extended to nonlinear optimization problems [22].There have been many applications of surrogate models in place of the finite element models.Ren and Chen [23], and Fang and Perera [24,25] use the response surface method for model correction and damage identification.The polynomial of a particular order is used to fit the relationship between the parameters and the responses.The least-squares method is used to determine the coefficients of the polynomial based on reasonable initial sampling [23][24][25][26][27] (using D-optimal design or central composite design).A small number of sampling points can be effectively used to predict the response of the structure under different parameters.In addition, Chakraborty and Sen [28] used the moving least-squares method to solve the coefficients of the response surface polynomial in order to increase the influence of the sampling points near the prediction point on its prediction value.Vincenzi and Savoia [29] combined the response surface method with differential evolution algorithm for effective parameter identification.Zhou et al. [30] compared the effects of constructing response surfaces using different radial basis functions to update a large span cable-stayed bridge.
Furthermore, as a special response surface construction method similar to the Gaussian basis function [31], the Kriging model has received more and more attention in recent years [32].Khodaparast, et al. [33] used the Kriging model to establish the relationship between the positional parameters of two beams and the frequency of the whole structure, for fast interval model updating to define the uncertain range of the parameters.Liu et al. [34] proposed a 2-level model updating method based on Kriging predictor for large complex structures.Wang, Wang and Zhao [17] introduced the Kriging model into an acceleration frequency response function-based model updating method, which is performed successfully on a honeycomb sandwich beam.Jin and Jung [35] proposed a sequential surrogate modeling method.On the basis of a certain amount of initial sampling points, this method constructs response surfaces for different responses respectively through a specific adding point criterion.
Most of the methods mentioned above have chosen the frequency as the output response of the surrogate model.This is mainly because the response surface of a frequency has a relatively gentle shape with respect to different system parameters.Particularly when a polynomial is used as the mathematical construct of a surrogate model, the frequency is not sensitive to the cross terms of the polynomial.So an accurate response surface can be built from a small number of sample points.On the other hand, the frequency response function (FRF) can provide more structural vibration information than the frequency and its use as response in surrogate modeling is worth investigating, which has rarely been attempted.Therefore, in this paper, a Kriging model with good simulation effect of nonlinear response surface is constructed to establish the relationship between the updating parameters and the residual of the Frequency Domain Assurance Criterion (FDAC) value, and a model updating method is used on the Kriging model for damage identification.The efficient global optimization (EGO) algorithm is introduced to realize fast convergence of the updating procedure.The proposed method is applied to a simulated example of a cantilever beam with three weakened elements and a laboratory three-story frame structure to verify its effectiveness.

Kriging Based Model Updating Method
In this section, the theory of the Kriging based model updating method is presented.Generally, it requires two steps to build a precise Kriging model.The first step is to build a Kriging model through initial sample points generated by design of experiments (DOE).The second step is to add a new sample point by a specified infill criterion and rebuild the Kriging model until the precision of the prediction value meets the requirements.
In this paper, the Latin hypercube sampling (LHS) method is used to produce initial sample points, the number of which is usually 10 times the dimension.The infill criterion is described in Section 2.3.

Construction of Kriging Model
A Kriging model is a surrogate model based on a stochastic process.It was originally put forward in geostatistics and made its way into engineering design following the work of Sacks, et al. [36].For a given set of sample data (input), X = {x 1 , x 2 , • • • , x n } T , and the observed responses (output), Y = {y 1 , y 2 , • • • , y n } T , the expression of Kriging model that reflects the relationship between them is: where f(x) is a polynomial vector of the sample x, β is the vector of the linear regression coefficients to be estimated and z(x) represents errors and is assumed to be a stochastic process that follows a normal distribution of N(0, σ 2 ) with a zero mean and standard deviation σ.
It should be noted that the basic assumption of the Kriging model is that the same input will lead to an identical output.Therefore, the deviation between the output response and the polynomial regression part is only due to the modeling error itself, regardless of the measurement error and other random factors.This method does not depend on the simulated precision of the polynomial part to the response surface, but focuses on constructing the appropriate surrogate model by effective filling of the stochastic process part, which makes it more suitable for dealing with nonlinearity.Thus, the polynomial part is often taken as a constant in some other references.
To estimate the stochastic process z(x), Kriging assumes that the true response surface is continuous, any two points will tend to have the same value as the distance in between approaches zero and it is the same for z(x) of two points.Thus, the correlation between z(x) of any two sample points can be expressed as a function of their spatial distance.The most widely used Gaussian correlation model is adapted here: where x k i and x k j are the kth components of the two sample points x i and x j , m denotes the number of design variables, θ k controls the decay rate of correlation on different dimensions.And then the matrix of correlation functions between sample points is obtained as The likelihood function of the sample point can then be written as: where F is a matrix of vector f(x) for each sample point.|R| is the determinant of R which is a function of θ k .According to the maximum likelihood function method, one can get: Based on this, the logarithmic form of the maximum likelihood function can be written as: The maximum value of the function above is solved by the genetic algorithm to determine the value of the decay rate θ k on different dimensions.
At this point, a Kriging model linking the sample point and the response is constructed.The next step is to predict the value of new points.For any point x 0 , following the principle that the predicted value for the point continue to maximize the augmented likelihood function of both the sample point and the new point, the predicted response value can be obtained by: And the mean squared error (MSE) of the predictor can also be calculated to estimate the accuracy of the predicted value, which is denoted by ŝ2 (x): where r T (x 0 ) is a row vector of correlation function between the new point and each sample point: It is worth noting that when the value of the ith sample point is predicted, since r T (x i )R −1 equals the ith order unit vector (whose elements are all zero except the ith element which is one), Equation (7) which shows that the Kriging model predicts the real response value at the sample point.That is why it can be considered an interpolation technique.Interested readers can refer to references [37,38] for a more detailed presentation of the theory, and to reference [36] for a more detailed formulation, based on the idea of best linear unbiased predictor.

Kriging Model in Model Updating
Model updating technique is to make the theoretical model, usually the FE model, produce a response as closely as possible to the experimental data.The traditional sensitivity-based model updating method uses the sensitivity matrix to correlate the updating parameters with the residuals of the response, which is the first derivative of the response residual with respect to the updating parameters, and then iteratively solves the associated optimization by the (weighted) least-squares method.However, this method usually requires the knowledge of all the degrees of freedom (DoFs) of the finite element model.When some DoFs cannot be measured, the model condensation technique shall be used and if the number of measured response points is much smaller than the number of the model degrees of freedom, this method will fail to give a unique solution.
Model updating can be essentially considered an optimization problem which can be formulated in the following form: where f (x) is the objective function and usually the response residual between theoretical prediction and experimental data, design variable x j is the jth parameter of the structure to be updated, lb j and ub j are the lower and upper bound of x j , respectively.By solving this optimization problem, design variable vector x which minimizes the objective function (response residual) gives the right updating parameter values.Besides using the sensitivity method which can also be regarded as a steepest descent method, the artificial intelligence algorithm is also convenient.In this paper, the widely used genetic algorithm (GA) is selected.The advantage of this method is that it is easy to update some special parameters of which the sensitivity matrices are difficult to calculate and to specify the interval of parameters to prevent them from convergence to meaningless results.
Of course, due to the strong complementarity among the updating parameters, the optimization problem of model updating is a typical multi-peak problem.The algorithm can easily cause convergence to a local solution if the number of updating parameters is large.In addition, the core problem of this approach is that a huge number of objective functional values at different sets of design variables are required by the GA to find the global minimum.If the computation time taken to evaluate a single value of the objective functional from the theoretical model is long, the total duration of solution by the algorithm may be unacceptable.That is why a Kriging model is required here to calculate the objective functional instead of the original FE model.
There are two ways to use the Kriging model to improve the computational efficiency in the model updating procedure.Firstly, the response of the FE model is used as the output response.A precise Kriging model is constructed for every measured response to totally replace the FE analysis in the subsequent model updating procedure.However, this kind of usage is mainly for the frequency as response and for few updating parameters.When it comes to FRF, for Kriging model construction, not only many measured responses are required, but also many sample points are required, since the true response surface of the FRF is much more complicated than that of the frequency.
These are the main reasons why, in this paper, the second approach is taken in which the value of the objective function of the optimization is directly used as the output response of the Kriging model.In so doing, only one Kriging model is required to be built and the EGO algorithm is adopted to obtain the model updating results, which is described in the next subsection.
However, it is unwise to set the residual of FRF between the theoretical and experimental FRF models as the objective function, because there is easily a big difference in magnitude between the two FRFs at a specified frequency, and in this case, their residual may overwhelm the influence of the smaller residuals at other frequencies.Hence, Frequency Domain Assurance Criterion, proposed by Pascual [39,40] to compare the correlation between two FRFs, is introduced to be the objective function: where n is the number of selected frequencies, ω i is the ith frequency, and w i is the weighting.The revised FDAC given in [41] is adopted: where h is a specified column of the FRF matrix, ω is the frequency, subscripts "a" and "x" correspond to the theoretical and experimental data respectively and superscripts "R" and "I" denote the real part and imaginary part of the FRF vector h, respectively.It follows from Equation ( 13) that FDAC values are limited to the interval [-1,1], regardless of the modulus value difference between two vectors.There are 2 more main advantages of the proposed objective function: 1.
It contains more physical meaning during the model updating procedure, which can improve the reliability of the updating results.2.
The shape of response surface becomes more gentle, which reduces the fitting difficulty of the Kriging model.
Moreover, unlike the sensitivity based method, it is not necessary to use FRFs of all the DoFs to calculate the FDAC value.The consistency of several sets of FRF curves is sufficient to ensure correct updating results, which allows elimination of the measured data with low signal to noise ratio.

Efficient Global Optimization
The EGO algorithm was proposed by Jones, Schonlau and Welch [38] as an effective method using response surfaces for global optimization.In that method, a Kriging model was constructed between the design variables and the objective function of the optimization problem, and only the accuracy of the location where the minimum value might exist was of interest.By only adding sample points in these effective areas, the purpose of solving the optimization problem using only a small number of sample points is achieved.
After constructing the Kriging model by the initial sample points, it is easy to get the prediction and its MSE of any point by Equations ( 7) and (8).The EGO algorithm treats the uncertainty in the prediction as a Gaussian probability density function with the predicted value as mean value and the MSE as variance, as shown in Figure 1b.When finding a new point to be added in the sample points, the algorithm does not simply select the point with the minimum value of the current Kriging model, but rather, selects the point with the maximum value of the expected improvement (EI) relative to the minimum value of the current responses at all the sample points.Based on the Gaussian probability density function, EI is calculated by: where y min is the minimum response value of the sample points, I(x) is the improvement of the minimum value by adding a new point, erf(•) is the cumulative distribution function, and EI(x) is the expectation of that improvement.Usually, the smaller the predicted value of the unknown point is, the larger the MSE is, and the larger the value of EI is, as shown in Figure 1a, in which the right axis is only for the expected improvement.The calculation of expected improvement at x = 0.7 is illustrated in Figure 1b.Only the red region is valid for EI calculation, which is controlled by y min .Thus, the infill criterion of the EGO algorithm not only finds the minimum value of the Kriging model (local search), but also explores the region with high MSE (global search).After adding enough sample points, the Kriging model will be accurate in the region where the minimum of the objective function exits and the global optimization problem is solved efficiently.In this paper, the solution for the location of the maximum expectation is also done by genetic algorithm.And the stopping criterion is set as: where 1  is a small positive number that controls the precision and the convergence speed of the algorithm.However, a shortcoming of the EGO algorithm in this study is that its precision of the local minimum value is not particularly high since it puts a certain amount of effort into the global search.Thus it can be considered a step to find the area in which the global minimum exists and a supplementary process is implemented which keeps adding the point with the minimum value of the Kriging model as a new sample point until the following criterion is reached: where  is the Euclidean distance between the new sample point new x and min x whose output response is minimum among current sample points, and small positive number 2  denotes the tolerance for termination of iteration.
Based on all methods mentioned above in Section 2, the flow chart of the complete model updating procedure is presented in Figure 2, and can be described as follows: Step 1 Generate initial sample points of the updating parameters by DOE.
Step 2 Run the FE analysis program to calculate the objective function output vector of the sample points and construct the initial Kriging Model.
Step 3 Find the point with the maximum expected improvement value of the current Kriging model and add it into the set of sample points.
Step 4 Calculate the true response of the new sample point and reconstruct the Kriging model by the new set of sample points and their response.Thus, the infill criterion of the EGO algorithm not only finds the minimum value of the Kriging model (local search), but also explores the region with high MSE (global search).After adding enough sample points, the Kriging model will be accurate in the region where the minimum of the objective function exits and the global optimization problem is solved efficiently.In this paper, the solution for the location of the maximum expectation is also done by genetic algorithm.And the stopping criterion is set as: max(EI(x)) < ε 1 (15) where ε 1 is a small positive number that controls the precision and the convergence speed of the algorithm.However, a shortcoming of the EGO algorithm in this study is that its precision of the local minimum value is not particularly high since it puts a certain amount of effort into the global search.Thus it can be considered a step to find the area in which the global minimum exists and a supplementary process is implemented which keeps adding the point with the minimum value of the Kriging model as a new sample point until the following criterion is reached: where • is the Euclidean distance between the new sample point x new and x min whose output response is minimum among current sample points, and small positive number ε 2 denotes the tolerance for termination of iteration.
Based on all methods mentioned above in Section 2, the flow chart of the complete model updating procedure is presented in Figure 2, and can be described as follows: Step 1 Generate initial sample points of the updating parameters by DOE.
Step 2 Run the FE analysis program to calculate the objective function output vector of the sample points and construct the initial Kriging Model.Step 7 Check whether the second termination criteria (Equation ( 16)) is satisfied.If it is, then stop updating the Kriging model and the sample corresponding to the minimum value in the output vector is used as the right updating result.Otherwise, go back to Step 4 and continue adding new points.

Numerical Example
The proposed model updating method is first verified through a simulated numerical example of a cantilever beam of rectangular cross-section divided into 9 elements, as shown in Figure 3.The properties of the beam are 0.9 m in length, 200 GPa in Young's modulus, 7670 kg/m in mass density, and 6 mm and 50 mm in height and width respectively.

Numerical Example
The proposed model updating method is first verified through a simulated numerical example of a cantilever beam of rectangular cross-section divided into 9 elements, as shown in Figure 3.The properties of the beam are 0.9 m in length, 200 GPa in Young's modulus, 7670 kg/m 3 in mass density, and 6 mm and 50 mm in height and width respectively.

Numerical Example
The proposed model updating method is first verified through a simulated numerical example of a cantilever beam of rectangular cross-section divided into 9 elements, as shown in Figure 3.The properties of the beam are 0.9 m in length, 200 GPa in Young's modulus, 7670 kg/m in mass density, and 6 mm and 50 mm in height and width respectively.A total of 9 Euler-Bernoulli beam elements are used and each node has two degrees of freedom of transverse displacement and rotation.It is assumed that the damage of the structure results in stiffness reduction in element 2 by 20%, in element 4 by 25% and in element 5 by 10%.The natural frequencies of the undamaged model and the damage model are shown in Table 1.A total of 9 Euler-Bernoulli beam elements are used and each node has two degrees of freedom of transverse displacement and rotation.It is assumed that the damage of the structure results in stiffness reduction in element 2 by 20%, in element 4 by 25% and in element 5 by 10%.The natural frequencies of the undamaged model and the damage model are shown in Table 1.The stiffness coefficients of all elements are used as the updating parameters and their change interval is set to (0.7, 1).The number of frequency selected for FRF should be more than the number of the updating parameters to reduce the effect of noise and excellent frequency selection away from resonance and anti-resonance can effectively improve the reliability and convergence speed of the algorithm.However, to check the robustness of the algorithm, the frequency of the FRF is simply set from 20 Hz to 800 Hz at every 50 Hz.
A total of 90 samples are obtained by LHS as the initial sample population and the Kriging model based model updating is carried out with the stopping criterion of ε 1 = 10 −10 in Equation ( 15) and criterion of ε 2 = 10 −4 in Equation (16).At the same time, the model updating based on GA directly is also implemented as a benchmark to compare the results of the two methods.These two methods are herein after referred to as Kriging method and GA method respectively.The GA used in this paper is the ga function of MATLAB 2017 (R2017a, The Math Works Inc., Natick, MA, USA, 2017) whose Function Tolerance value is set to 10 −9 as the stopping criterion.For details, please consult Matlab Help menu.
In this example the Kriging model is constructed after 358 iterations.As the dimension of the parameter is high, only the stiffness coefficients of element 1 and 2 are taken as variables while the stiffness coefficients of the other elements are set to 1 to demonstrate the similarity between the response surface predicted by the Kriging model and the real response surface as shown in Figure 4.It can be seen that when there are many parameters to be updated, the prediction error of the objective function of the Kriging model is large in the location away from the minimum.
In this example the Kriging model is constructed after 358 iterations.As the dimension of the parameter is high, only the stiffness coefficients of element 1 and 2 are taken as variables while the stiffness coefficients of the other elements are set to 1 to demonstrate the similarity between the response surface predicted by the Kriging model and the real response surface as shown in Figure 4.It can be seen that when there are many parameters to be updated, the prediction error of the objective function of the Kriging model is large in the location away from the minimum.Furthermore, the complexity of the real response surface when all the parameters are variables can be envisaged through Figure 4b.Even building that 2-dimensional surface requires a considerable number of sample points.That is why the EGO algorithm is adopted to save the computational cost.The final updating results are shown in Figure 5.For the sake of comparison, a conventional way of Furthermore, the complexity of the real response surface when all the parameters are variables can be envisaged through Figure 4b.Even building that 2-dimensional surface requires a considerable number of sample points.That is why the EGO algorithm is adopted to save the computational cost.The final updating results are shown in Figure 5.For the sake of comparison, a conventional way of using FRF in the Kriging based model updating is also applied, in which the amplitude of FRF is used to construct the objective function, which is referred to as the AF approach in the paper.using FRF in the Kriging based model updating is also applied, in which the amplitude of FRF is used to construct the objective function, which is referred to as the AF approach in the paper.As illustrated in Figure 5, the proposed Kriging method can effectively identify the stiffness reduction caused by damage in this numerical simulation experiment.A great improvement has been made by introducing the FDAC in the objective function as the AF approach contains obvious errors.Actually, the AF approach can hardly get correct results as using the FRF to construct the objective function directly will lead to a too complicated response surface for Kriging to build and this method will not be mentioned in the next case.On the other hand, although the Kriging method contains slightly higher errors in results compared with the GA method directly, it needs far fewer number of times running the FE program than the latter, as shown in Table 2, which is the core purpose of this method.

Experimental Study
In this section, the experimental data [42] of a three-story building structure made of aluminum As illustrated in Figure 5, the proposed Kriging method can effectively identify the stiffness reduction caused by damage in this numerical simulation experiment.A great improvement has been made by introducing the FDAC in the objective function as the AF approach contains obvious errors.Actually, the AF approach can hardly get correct results as using the FRF to construct the objective function directly will lead to a too complicated response surface for Kriging to build and this method will not be mentioned in the next case.On the other hand, although the Kriging method contains slightly higher errors in results compared with the GA method directly, it needs far fewer number of times running the FE program than the latter, as shown in Table 2, which is the core purpose of this method.

Experimental Study
In this section, the experimental data [42] of a three-story building structure made of aluminum tested at Los Alamos National Laboratory (LANL), shown in Figure 6, is used to verify the effectiveness and the practicability of the proposed approach.The test-bed structure consists of 4 plates and 12 columns assembled by bolted joints and was installed on rails that allow it to slide in only one direction.Additionally, a center column suspended from the top floor and an adjustable bumper mounted on the 2nd floor was used to simulate the nonlinear behavior of damage.by adding a 1.2 kg concentrated mass to the aluminum plates, while the reduction of stiffness was introduced by replacing the corresponding column with another one with half the cross section thickness in the direction of shaking.Each state was tested separately and the detailed state conditions [43] are shown in Table 3.As recommended in [42], the structure is modeled as four lumped masses at the floors shown in Figure 7.All the mass and stiffness values are selected as the updating parameters except since the friction between the rails and the structure is negligible.Besides, the damping matrix is introduced based on the proportional damping assumption which can also be affected by the mass factors.
At first, the Baseline condition is used to update the theoretical model to get the original undamaged model.The FRF of the first floor after the updating procedure compared with the experimental data is illustrated in detail in Figure 8.It can be clearly seen that the model after updating shows good agreement with the experimental data.However, there is considerable noise An electrodynamic shaker was used to provide a band-limited random lateral excitation (20-150 Hz) to the base floor along the center line of the structure and a force transducer was connected to the tip of the stinger to measure the input force.At the same time, four accelerometers were attached at the center line of each floor on the opposite side from the excitation to measure the dynamic response of the system.The analog sensor signals were discretized into 8192 data points sampled at 3.125 m intervals corresponding to a sampling frequency of 320 Hz.
Force-and acceleration-time histories for 17 different structural state conditions were collected.However, only 9 of them are utilized as the others are concerned with the nonlinearity which is not within the scope of this article.The 9 different states can be divided into 3 categories: the baseline structure, the increase of the mass and the reduction of the stiffness.The change of mass was realized by adding a 1.2 kg concentrated mass to the aluminum plates, while the reduction of stiffness was introduced by replacing the corresponding column with another one with half the cross section thickness in the direction of shaking.Each state was tested separately and the detailed state conditions [43] are shown in Table 3.As recommended in [42], the structure is modeled as four lumped masses at the floors shown in Figure 7.All the mass and stiffness values are selected as the updating parameters except k 0 since the friction between the rails and the structure is negligible.Besides, the damping matrix is introduced based on the proportional damping assumption which can also be affected by the mass factors.At first, the Baseline condition is used to update the theoretical model to get the original undamaged model.The FRF of the first floor after the updating procedure compared with the experimental data is illustrated in detail in Figure 8.It can be clearly seen that the model after updating shows good agreement with the experimental data.However, there is considerable noise in the imaginary part of the experimental FRF especially near the anti-resonance point and modeling error in the theoretical FRF around the same point.What is more, there is also a big difference between the real parts of two FRFs near 22 Hz, which should not be selected as the frequency points of the objective function, since such a big difference cannot be eliminated by model updating.
Then the two methods are applied on the undamaged model obtained using the experimental data of the other states.The intervals of the mass and stiffness parameter ratios are set to [0.98, 1.5]  and [0.4, 1.0], respectively.This setting is based on the pre-knowledge of the structure's change and prevents all parameters from becoming larger at the same time or smaller at the same time, which would make little change in the value of the objective function.
Unlike the numerical simulation mentioned in Section 3, in this case study, a larger parameter range is beneficial in getting the right updating parameter values.However, this may lead to selected frequencies to become close to a theoretical resonant frequency, which would result in a sharply increased objective function value, as illustrated in Figure 9. Then the two methods are applied on the undamaged model obtained using the experimental data of the other states.The intervals of the mass and stiffness parameter ratios are set to 0.98, 1.5 and 0.4, 1.0 , respectively.This setting is based on the pre-knowledge of the structure's change and prevents all parameters from becoming larger at the same time or smaller at the same time, which would make little change in the value of the objective function.Unlike the numerical simulation mentioned in Section 3, in this case study, a larger parameter range is beneficial in getting the right updating parameter values.However, this may lead to selected frequencies to become close to a theoretical resonant frequency, which would result in a sharply increased objective function value, as illustrated in Figure 9. Since the experimental data has considerable noise, the number of selected frequency points for model updating should be much higher than the number of the updating parameters to reduce the effects of noise and damping.If such an increase appears in too many frequency points, there will be severe fluctuations in the response surface of the final objective function formed at these frequency points, which are undesirable and useless in identification.
It is not difficult to relate this increase to the sharp decline of the real parts of the FRFs near the natural frequency as shown in Figure 8.In contrast, the imaginary parts of the FRFs show better correlation near the natural frequency.Thus if the selected frequency is near the natural frequency of experimental data, 3 Hz on either side of a peak of the real part of the FRFs is used as a criterion in this case, the imaginary part of the FRFs only is chosen to calculate the FDAC to get a smooth response surface.Since the experimental data has considerable noise, the number of selected frequency points for model updating should be much higher than the number of the updating parameters to reduce the effects of noise and damping.If such an increase appears in too many frequency points, there will be severe fluctuations in the response surface of the final objective function formed at these frequency points, which are undesirable and useless in identification.
It is not difficult to relate this increase to the sharp decline of the real parts of the FRFs near the natural frequency as shown in Figure 8.In contrast, the imaginary parts of the FRFs show better correlation near the natural frequency.Thus if the selected frequency is near the natural frequency of experimental data, 3 Hz on either side of a peak of the real part of the FRFs is used as a criterion in this case, the imaginary part of the FRFs only is chosen to calculate the FDAC to get a smooth response surface.
It turns out that this measure reduces the complexity of the frequency response surface.This improvement is shown in Figure 10.It is clear that the response surface is much smoother than before after such a special treatment which saves a great deal of computation time in constructing the Kriging model, because the fluctuating portion of the response surface is largely reduced.It should be noted that both the Kriging model based method and the GA method benefit from the improvement, namely that their computing time decreases significantly.Since the experimental data has considerable noise, the number of selected frequency points for model updating should be much higher than the number of the updating parameters to reduce the effects of noise and damping.If such an increase appears in too many frequency points, there will be severe fluctuations in the response surface of the final objective function formed at these frequency points, which are undesirable and useless in identification.
It is not difficult to relate this increase to the sharp decline of the real parts of the FRFs near the natural frequency as shown in Figure 8.In contrast, the imaginary parts of the FRFs show better correlation near the natural frequency.Thus if the selected frequency is near the natural frequency of experimental data, 3 Hz on either side of a peak of the real part of the FRFs is used as a criterion in this case, the imaginary part of the FRFs only is chosen to calculate the FDAC to get a smooth response surface.
It turns out that this measure reduces the complexity of the frequency response surface.This improvement is shown in Figure 10.It is clear that the response surface is much smoother than before after such a special treatment which saves a great deal of computation time in constructing the Kriging model, because the fluctuating portion of the response surface is largely reduced.It should be noted that both the Kriging model based method and the GA method benefit from the improvement, namely that their computing time decreases significantly.Finally, the updating results of each state by the two methods are shown in Figure 11.
It is clear that the updating results of both methods match the real values quite well.The parameter that has experienced the main change can be clearly identified with small false alarms in other parameters, due to the measurement noise and the modeling error in the theoretical model.The successful application in these states is sufficient to demonstrate the superiority of the method using FDAC value to build the objective function.Furthermore, the Kriging method shows similar updating results to those by the GA method, which means the proposed method has the advantage of a typical good surrogate model and ensures the updating accuracy at the same time.
Additionally, the required numbers of the finite element runs of all states by the two methods are listed in Table 4. Again, the Kriging method shows a great improvement over the GA method.It should be noted that the proposed method can also converge very quickly if a point near the global minimum is selected as one of the initial sample points.
Finally, the updating results of each state by the two methods are shown in Figure 11.

Conclusions
In this study, a new surrogate model based model updating method using frequency response function (FRF) is proposed.Although FRF contains the most vibration information, its response surface is difficult to construct directly and the Frequency Domain Assurance Criterion (FDAC) is adopted to build the objective function rather than the residual of every single FRF.As a combination of three parts, the residual of FDAC for the objective function and the Kriging model for constructing the response surface and the efficient global optimization (EGO) for solving the optimization problem, this method is found, through a simulated example and a real experiment, to have the following merits: (1) it does not need to calculate the sensitivity of the parameters, ( 2) it has similar accuracy as the direct use of GA, but requires much fewer runs of the finite element analysis, and (3) it has strong robustness in the face of measurement noise and modeling errors, which are always present and cannot be ignored.

Figure 1 .
Figure 1.(a) The expected improvement of a one-variable function calculated by Kriging model constructed through five sample points; (b) A graphical interpretation of the expected improvement at x = 0.7.

Figure 1 .
Figure 1.(a) The expected improvement of a one-variable function calculated by Kriging model constructed through five sample points; (b) A graphical interpretation of the expected improvement at x = 0.7.

Appl. Sci. 2017, 7 , 1039 8 of 18 Step 3
Find the point with the maximum expected improvement value of the current Kriging model and add it into the set of sample points.Step 4 Calculate the true response of the new sample point and reconstruct the Kriging model by the new set of sample points and their response.Step 5 Check whether the first termination criterion (Equation (15)) is satisfied.If it is, go to Step 6.Otherwise, go back to Step 3 and continue adding new points.Step 6 Find the point with the minimum value of the current Kriging model and continue to add it into the set of sample points, calculate the true response output and reconstruct the Kriging model.Step 7 Check whether the second termination criteria (Equation (16)) is satisfied.If it is, then stop updating the Kriging model and the sample corresponding to the minimum value in the output vector is used as the right updating result.Otherwise, go back to Step 4 and continue adding new points.Appl.Sci.2017, 7, 1039 8 of 17

Figure 2 .
Figure 2. The flow chart of the proposed model updating method.

Figure 3 .
Figure 3.The cantilever beam used for numerical simulation.

Figure 2 .
Figure 2. The flow chart of the proposed model updating method.

Figure 3 .
Figure 3.The cantilever beam used for numerical simulation.

Figure 3 .
Figure 3.The cantilever beam used for numerical simulation.

Figure 4 .
Figure 4.The response surface of the objective function when only two parameters are variable: (a) The Kriging model prediction; (b) Real response surface.

Figure 4 .
Figure 4.The response surface of the objective function when only two parameters are variable: (a) The Kriging model prediction; (b) Real response surface.

Figure 5 .
Figure5.The updating results of the Kriging method, genetic algorithm (GA) method and AF approach compared with the 'real' damage.

Figure 5 .
Figure5.The updating results of the Kriging method, genetic algorithm (GA) method and AF approach compared with the 'real' damage.

17 Figure 7 .
Figure 7. Theoretical model for the test structure.

Figure 7 .
Figure 7. Theoretical model for the test structure.

Figure 7 .Figure 8 .
Figure 7. Theoretical model for the test structure.

Figure 8 .
Figure 8.The frequency response function (FRF) of the first floor of the real structure and the updated model: (a) Absolute value of H(1, 1); (b) Absolute value of the real part of H(1, 1); (c) Absolute value of the imaginary part of H(1, 1).

Figure 9 .
Figure 9.The objective function of 50 Hz when only m1 and m2 are variable using experimental data of state 1.

Figure 9 .
Figure 9.The objective function of 50 Hz when only m1 and m2 are variable using experimental data of state 1.

Figure 9 .
Figure 9.The objective function of 50 Hz when only m1 and m2 are variable using experimental data of state 1.

Figure 10 .
Figure 10.The response surface of the objective function when only m1 and m2 are variable: (a) Before improvement; (b) After improvement.

Figure 10 .
Figure 10.The response surface of the objective function when only m1 and m2 are variable: (a) Before improvement; (b) After improvement.

Figure 11 .
Figure 11.Updating results of the different states of the two methods compared with the real value.(a)−(h) are corresponding to the state from 2 to 9, respectively.

Figure 11 .
Figure 11.Updating results of the different states of the two methods compared with the real value.(a)−(h) are corresponding to the state from 2 to 9, respectively.

Table 1 .
Natural frequencies of the theoretical and damaged model of the cantilever (Unit: Hz).

Table 2 .
The required number of running the finite element (FE) program of the Kriging method, and the GA method (which is repeated 3 times).

Table 2 .
The required number running the finite element (FE) program of the Kriging method, and the GA method (which is repeated 3 times).

Table 3 .
Summary of structural state conditions.

3 .
Summary of structural state conditions.

Table 4 .
The required number of running the FE programs of Kriging method and the GA method for different states.