Probabilistic Structural Model Updating with Modal Flexibility Using a Modified Firefly Algorithm

Structural model updating is one of the most important steps in structural health monitoring, which can achieve high-precision matching between finite element models and actual engineering structures. In this study, a Bayesian model updating method with modal flexibility was presented, where a modified heuristic optimization algorithm named modified Nelder–Mead firefly algorithm (m-NMFA) was proposed to find the most probable values (MPV) of model parameters for the maximum a posteriori probability (MAP) estimate. The proposed m-NMFA was compared to the original firefly algorithm (FA), the genetic algorithm (GA), and the particle swarm algorithm (PSO) through the numerical illustrative examples of 18 benchmark functions and a twelve-story shear frame model. Then, a six-story shear frame model test was performed to identify the inter-story stiffness of the structure in the original and the damage states, respectively. By comparing the two, the position and extent of damage were accurately found and quantified in a probabilistic manner. In terms of optimization, the proposed m-NMFA was powerful to find the MPVs much faster and more accurately. In the incomplete measurement case, only the m-NMFA achieved target damage identification results. The proposed Bayesian model updating method has the advantages of high precision, fast convergence, and strong robustness in MPV finding and the ability of parameter uncertainty quantification.


Introduction
Structural model updating is an important part of structural health monitoring. Using monitoring data, structural models can be updated to more realistically represent the actual state of the structures and to obtain mechanical analytical results which are closest to the real state of the structures. In fact, structural model updating is an inverse problem of structural analysis, which can be regarded as an optimization problem. Since the 1960s, the field of structural model updating has developed a variety of methods. However, errors exist widely in experimental testing and modeling which are caused by various factors [1,2]. Shinozuka [3] and Collins et al. [4] introduced an uncertainty concept and used a probability statistical method to evaluate the uncertainty of parameters while updating the model parameters. As a result, structural model updating has had a milestone leap.
Probabilistic model updating methods can be roughly divided into stochastic finite element corrections and finite element corrections based on Bayesian methods. Among them, stochastic finite element correction is a large category, which includes model correction based on the Monte Carlo method and model correction based on the perturbation method. Beck and Katafygiotis et al. [5,6] first introduced the Bayesian method into model updating in 1998 and established a preliminary theoretical framework; they then used the asymptotic approximation method to estimate a posterior probability in the parameter space. Researchers including Beck, Katafygiotis, Yuen, and Feng have conducted a series of in-depth studies on the problems of incomplete measurement, complex likelihood function, and posterior probability that are difficult to solve, and proposed "system modal shape", modal flexibility method, local observation concepts and methods, such as data fusion, to further enrich the relevant theories [7][8][9][10][11][12].
Based on the Bayesian model updating framework, model parameters are embedded in the posteriori probability function. Generally, there are two kinds of methods to solve the maximum a posteriori probability (MAP) estimate, one is the Laplace approximation and the other is random sampling of the posteriori probability distribution function. For a globally identifiable case with a large number of data, we can approximate the posterior probability distribution function as a Gaussian distribution including mean and variance, which mean is the most probable value (MPV), and the obtained MPV is the MAP estimate for the model parameters. Because MPV takes the maximum function value for the posterior probability distribution function, the solution of MPV eventually becomes an optimization problem. After completing the construction of a mathematical model for the optimization problem, the optimization algorithm can be used to find the optimal value of the mathematical model. Due to the limitations of analytical methods, numerical methods have become the most widely used methods for solving optimization problems. They can be mainly divided into two categories. The first type includes traditional optimization algorithms focusing on single-objective optimization problems, such as the Newton Iteration Method (NIM), the Conjugate Direction Method (CDM), etc. The other type includes Heuristic Algorithms (HA) for multi-objective optimization problems. Inspired by bionics, scholars have developed Metaheuristic Algorithms (MA) [13,14]. A meta-heuristic optimization algorithm allows randomness in the optimization process of the algorithm by introducing a mutation factor, thereby reducing the possibility of falling into a local optimum. At present, scholars have proposed many meta-heuristic optimization algorithms, such as the Genetic Algorithm (GA) [15], the Ant Colony Optimization (ACO) [16], the Firefly Algorithm (FA) [17], etc.
The firefly algorithm was first proposed by Yang in 2007. It is a global optimization algorithm inspired by the glowing behavior of fireflies in nature. Since the firefly algorithm was proposed, due to its simple parameters, fast convergence speed, high precision, as well as its own automatic population segmentation ability and nonlinear individual attraction mechanism, it has received extensive attention from the academic community. However, like other swarm intelligence optimization algorithms, this algorithm is prone to problems such as premature maturity, falling into local optimum, and insufficient ability to solve multi-dimensional optimization problems. Therefore, it needs to be appropriately improved when solving actual complex optimization problems.
Based on the original FA algorithm, this paper makes several modifications and introduces the Nelder-Mead Algorithm (NMA) for local search and proposes a modified Nelder-Mead firefly algorithm. Three improvements are proposed: modification in the control parameters, introduction of a boundary constrain mechanism, and introduction of a diversity threshold. Eighteen benchmark functions and a 12-story shear frame model are used to test the performance of the m-NMFA and to provide support for the comparison of the original FA, the GA, and the PSO to the m-NMFA. A six-story steel shear frame was designed to perform a vibration experiment for verifying the effectiveness of the m-NMFA in Bayesian structural model updating.

Modal Flexibility
The stiffness matrix K and flexibility matrix C of a linear dynamical system with N degrees of freedom can be written as follows [18]: where M and K are the mass matrix and the stiffness matrix, respectively; Φ = [φ 1 , φ 2 , . . . , φ n ] is the mass-normalized mode shape matrix; and Λ is the diagonal eigenvalue with squared natural circular modal frequencies.
In practice, it is always the case that the spatial arrangement of sensors attached to the structure is generally limited and only a limited number of DOFs can be observed, for which only incomplete spatial mode shapes can be gained. It can be seen from Equation (2) that the greater the modal frequency is, the smaller the contribution of one mode to the flexibility matrix C will be, which is on the contrary to that of Equation (1) and is more practical in solving vibration-related problems. Therefore, the flexibility matrix C can be established approximately with only low-frequency modes of interest. Therefore, the Equation (2) yield as follows [19]: where C m denotes the flexibility submatrix with respect to the measured DOFs; C m t denotes the truncated flexibility submatrix; N m denotes the number of selected lowfrequency modes (frequency from low to high and no interval in this research); and φ m i ∈ R N o ×1 (i = 1, 2, . . . , N m denotes the incomplete mode shapes. Needless to say, C m t is supposed to be as close as possible to C m ; if all the DOFs of a structure are well observed and all the mode shapes are gained, C m t will be exactly equal to C m . Mass-normalization is a key issue to build up the flexibility matrix. If the measurement is incomplete, the dimension of modal parameters, such as mode shapes, would be less than that of the mass matrix M. To handle this problem, the SEREP method [20] is introduced to convert M to reduced mass matrix M m . Then, the incomplete mode shape matrixes are normalized as follows:

Model Reduction
System equivalent reduction and expansion process (SEREP) is an eigenvalue-based reduction technique used to approximate the real structural state vectors with mode shape measured from the DOFs of the user's interest [20].
In the method of SEREP, the first step is to separate the mass and stiffness matrices as follows: where m and u denote the measured and unmeasured DOFs. Then, the SEREP transformation matrix can be written as follows: where Φ m denotes the retained mode shape matrix with N o (the number of observed DOFs) rows and N m columns, and (Φ m ) * represents the generalized inverse of matrix Φ m . The process of converting Φ m to (Φ m ) * is shown in [11,20,21]. Ultimately, the reduced mass and stiffness matrices can be written as follows:

Bayesian Model Updating Based on Measured Modal Flexibility Data
Bayesian method for finite element model updating is one of the most commonly used probabilistic methods for structural model updating, and it has been confirmed to be efficient and useful [22,23]. In this section, modal flexibility data are used to establish a posterior-PDF-based objective function.
As the modal flexibility matrix is symmetrical, we can use the lower triangular matrix of the flexibility matrix C l to establish the objective function. Let f denote the modal flexibility data, and the matrix C l is vectorized based on the following equation: where c denotes the element in the matrix C l . It is supposed that the test is repeated N t times and D denotes the experimental data. Then, the matrix D can be written as follows: As the modal can be described by the probability distributions of both the unknown parameters and the prediction error [5], the measured data can be written as follows: where f (θ) denotes the prediction of the analytical model and is defined as the vector for the given values of the model parameters θ. The term e t denotes the prediction error and is assumed to follow the maximum entropy principle [17]. That is to say, e t yields a multi-dimensional Gaussian distribution with zero mean and covariance matrix Σ [5]. The covariance matrix is computed based on the modal flexibility vectors as follows: where f contains the average value of N t experimental modal flexibility vectors. It can be expressed using the following equation: Then, the likelihood function of the t-th set of modal flexibility vector data can be expressed as follows: where C 1 is a constant. Assuming that the N t sets of the modal flexibility vector datas are mutually independent, the likelihood function of D can be written as follows: where C 2 is also a constant like C 1 . As the prior PDF p(θ) follows a non-informative uniform distribution, the posterior PDF p(θ|D, M) can be written in a similar form as the likelihood function in Equation (14) Given that the observed DOFs are limited, the SEREP method will be utilized here and then the measured mode shapes are mass-normalized with respect to the reduced model mass M m . It is easy to see that M m depends on θ only and so does the posterior PDF p(θ|D, M). On this basis, the most probable stiffness parameter vectorθ can be gained by maximizing the posterior PDF value or minimizing its negative value. Here, the objective function is defined as a truncated part of Equation (14) and is simplified with the following logarithm: where J(θ) denotes the objective function. Substitute Equation (14) into Equation (15), the function is expressed as follows: Till now, the objective function has been established and the model updating problem has been transformed as an optimization problem. When the most probable stiffness parameter vectorθ is found during the process of minimizing J(θ), the actual stiffness matrix of a structure is ensured naturally.
The posterior PDF can be well approximated using a set of multi-dimensional Gaussian distribution withθ and Σ θ as its mean value and covariance matrix.
where H(θ) is the Hessian matrix of objective function.
In order to search for the best-fit vectorθ exactly and reliably, a modified FA algorithm is proposed, which is introduced in the next section and utilized in the numerical and experimental examples in Sections 4 and 5.

Original Firefly Algorithm
Firefly Algorithm (FA) was first proposed by Xin-She Yang at Cambridge University in 2007 [24] and has become a widely used global optimization algorithm. In the original FA, it is not gender difference but the intensity of light I, which decreases with the increase in the square of Cartesian distance between two fireflies r 2 , that distinguishes fireflies. Fireflies that have greater I values are regarded as fitting the solution to the objective function better, and they are supposed to attract the ones with smaller values to them. Given that light will gradually weaken in a medium as distance increases, the attractiveness parameter β, which denotes the attraction between any two fireflies, is introduced in this algorithm. Similar to the intensity of light, the greater the distance is, the smaller the attractiveness parameter β will be, and the attraction between two fireflies will become weaker. Furthermore, the light absorption parameter γ is introduced to denote the ability of a medium to absorb light, which is proposed to simulate the natural environment that fireflies are in.
The behavior of firefly is idealized for algorithmic use, and three rules are introduced as follows [24]:

•
All the fireflies are unisex so they can attract each other without the influence of gender.

•
The attractiveness is proportional to the brightness, and they both decrease as their distance increases.

•
The brightness of a firefly is determined by the landscape of the objective function.
In a given medium, light intensity I decreases in a Gaussian form as follows: where I 0 represents the intensity at r = 0, and r denotes the distance between any two fireflies. It can be seen that light intensity decreases exponentially and the influence of the absorption coefficient γ on light intensity is defined in the exponential term of the function mentioned above. However, light intensity is flexible in both definition and form since it is allowed to be the fitness value of a certain objective function to a practical problem [17].
The attractiveness β of a firefly is proportional to the light intensity I and is also defined in an exponential form as follows: where β 0 is the attractiveness value at r = 0. The distance r between any adjacent fireflies x i and x j is expressed as follows: where x i,k is the k-th component of the spatial coordinate x i of i-th firefly, and n denotes the dimension of the objective function or problems to be optimized. As mentioned above, fireflies with a greater value of light intensity I is supposed to attract the ones with a smaller value of I. Assume that the i-th firefly is weaker than the j-th, then the movement of firefly i to firefly j is described as follows: where the second term is the attraction part, and the third term is the randomization part with step parameter α. i is drawn from different distributions, such as Gaussian distribution, uniform distribution, and levy flight. In many cases, a uniform distribution is utilized as the simplest form and, thus, i becomes a vector of random numbers uniformly generated in [17,[25][26][27]. An example of a levy flight distribution is shown in [24]. The firefly algorithm is an efficient algorithm, and Algorithm 1's pseudocode is shown as follows:

Modified Nelder-Mead Firefly Algorithm
The modified Nelder-Mead firefly algorithm(m-NMFA) is generally a hybrid version combining the advantages of the local optimization algorithm, the Nelder-Mead algorithm (NMA), and the global optimization algorithm. In order to strengthen the superiority of the hybrid algorithm, three improvements are proposed: modification in the parameters α and β, introduction of a boundary constrain mechanism, and introduction of a diversity threshold.
Yang in 2008 [17] used a fixed-step parameter α = 0.2 to solve a four-peak function and then pointed out that a decreasing α could improve the convergence of the FA. Two examples are given as follows: and where α 0 is the initial step parameter and α ∞ is the final value. θ ∈[0,1] is the randomness reduction constant. α iterates as the iteration counter t increases. It can be seen that the step parameter α is limited between α 0 and α ∞ in Formula (22) and gradually approaches α ∞ . Similarly, α gradually approaches zero in Equation (23). Moreover, in Equation (23), the θ and α 0 values have great influence in the decay speed of step parameter α in the process of iteration, which will produce a unknown effect-positive or negative-on the solution result. In previous studies, different θ values are taken to cope with different problems. Yang and He [34] pointed out that θ could take values from 0.95 to 0.97, Manoharan [30] and Rizk-Allah [35] took θ as 0.9, and Wang [31] took θ as 0.7 and α 0 as 0.25. Moreover, a noticeable issue has been highlighted that the scaling parameter S d should be introduced in consideration of the actual scales of the problem of interest [36].
where d denotes the d-th dimension unknown variable, and u d and l d are the upper and lower limit values. Therefore, i in Equation (21) should be replaced by i S, where S={S d , d=1,2, . . . n}. Beside the strategies mentioned above, various versions of the step parameter iteration formula have been proposed in the last decade. Brajević [27] proposed a novel strategy that iterates α with the maximum number of interaction and the initial step parameter α 0 and it is written as follows: where κ is the identifier of each firefly and takes the dimension of the problem as its maximum value; MaxGen denotes the maximum number of interaction; and t denotes the current number of iterations. The step parameter decreases more quickly and will be much smaller than that in Equation (23). Ignoring the distinction among all the fireflies, Equation (25) can be expressed by another form as follows: where Θ = (1/90, 000) 2 ≈ 1.234 × 10 −10 , and ρ is defined as the reference coefficient and its value is 1. Previous tests have shown that the results drawn from an α that is of a slightly smaller value during the previous iterations (take 50 iterations for example) tend to outperform the results drawn from a greater α value in convergence speed. Previous tests have also shown that the convergence speed of the algorithm increases but the accuracy decreases gradually if ρ is taken to be an integer constant value ranging from 1 to 3. That is to say, the decrease rate of α in a later course of iteration should be well controlled to have a good balance between exploration and exploitation to the maximum. So, let ρ be expressed in a nonlinear form using the following equation: Figure 1 shows the decreasing process of the reference coefficient. Apart from the modification on the step parameter, extensive efforts have been made to modify the attraction parameter β and the light absorption parameter γ. If β0 and γ are taken to be constants, the distance between any two adjacent fireflies will approach zero gradually. As the iteration number increases, it will ideally follow the equation below: This ensures that the attraction within a firefly colony will always be a contributing factor. In practice, however, it would be a tough problem in high-dimensional problems if the distance between any two firefly individuals is too large at the very beginning of the iteration progress and the attraction is too weak to get the individuals closer. As a consequence, the algorithm stagnates. Moreover, the higher the dimension of a project is, the more easily this phenomenon will occur.
Selvarasu [33] proposed a modified β formula: where βmin and βmax are the lower and upper boundaries of the attraction parameter β. This strategy prevents the attraction from being excessively large or small. Herein, this paper utilizes Equation (29) to improve the performance of the FA in solving high-dimensional problems.

Boundary Constrain Mechanism
Similar to other swarm intelligence methods, such as the Bee Colony Algorithm (BCA), the Genetic Algorithm (GA), and the Particle Swarm Optimization (PSO), the FA will be ill-conditioned if the variables of firefly individuals exceed the feasible region. Yang introduced a boundary constrain mechanism that forces the outer variables to the border, and the formula is expressed as follows: where x i,j denotes the i-th firefly individual and the j-th variable, and l j and u j are the lower and upper boundaries of the j-th variable. On the one hand, this mechanism ensures that all the variables of one firefly are within the feasible range and, moreover, prevents the algorithm from falling into the local optimum or leading to a completely wrong Apart from the modification on the step parameter, extensive efforts have been made to modify the attraction parameter β and the light absorption parameter γ. If β 0 and γ are taken to be constants, the distance between any two adjacent fireflies will approach zero gradually. As the iteration number increases, it will ideally follow the equation below: This ensures that the attraction within a firefly colony will always be a contributing factor. In practice, however, it would be a tough problem in high-dimensional problems if the distance between any two firefly individuals is too large at the very beginning of the iteration progress and the attraction is too weak to get the individuals closer. As a consequence, the algorithm stagnates. Moreover, the higher the dimension of a project is, the more easily this phenomenon will occur.
Selvarasu [33] proposed a modified β formula: where β min and β max are the lower and upper boundaries of the attraction parameter β. This strategy prevents the attraction from being excessively large or small. Herein, this paper utilizes Equation (29) to improve the performance of the FA in solving highdimensional problems.

Boundary Constrain Mechanism
Similar to other swarm intelligence methods, such as the Bee Colony Algorithm (BCA), the Genetic Algorithm (GA), and the Particle Swarm Optimization (PSO), the FA will be ill-conditioned if the variables of firefly individuals exceed the feasible region. Yang introduced a boundary constrain mechanism that forces the outer variables to the border, and the formula is expressed as follows: where x i,j denotes the i-th firefly individual and the j-th variable, and l j and u j are the lower and upper boundaries of the j-th variable. On the one hand, this mechanism ensures that all the variables of one firefly are within the feasible range and, moreover, prevents the algorithm from falling into the local optimum or leading to a completely wrong solution.
On the other hand, it reduces the diversity of the firefly population to some extent. This mechanism has been utilized by Ivona Brajević [37] and Rajan [38]. The optimal solution of a multi-dimensional problem is determined by all the variables. It is possible that a newly generated firefly individual would deviate from the optimal solution with only several variables being updated. Therefore, a new boundary constraint mechanism is proposed. If any one of the variables overflows, a new firefly individual will be regenerated randomly, Algorithm 2 is: where rand ∈ [0, 1] denotes a random factor that is used to generate a new variable randomly, and F ∈ [0, 0.5] is another random factor that leads the new firefly towards the current global best solution x best . Both rand and F follow a uniform distribution.

Hybrid of Nelder-Mead Algorithm and the Diverse Threshold
The Nelder-Mead Algorithm is a local optimal algorithm that adapts itself to the local landscape and contracts on to the final minimum. This method for finding the local minimum is widely used and it can be realized expediently by using the built-in function "fminsearch" in MATLAB software. The way of combining the NMA and the FA is to start the NMA in the iterative process of the FA in order to dig out a better solution nearby the current global best solution found up to t-th (t ≤ MaxGen) iteration.
However, the start point of the NMA in the iterative process needs to be well weighted because when to enable the NMA will greatly affect the solution accuracy and convergence speed. If the NMA starts too early, the current best solution will be far from reliable and the NMA will run meaninglessly, which will result in inefficiency. Meanwhile, if it starts too late, the contribution of the solution searched by the NMA to the algorithm will be weakened. Therefore, the concept of diversity threshold is introduced. When the diversity of a firefly population reaches a certain threshold, the NMA is enabled to balance the contradiction mentioned above.
Let x k,best and x k,worst denote the best and worst firefly individuals of the single k-th iteration, and x g,best and x g,worst denote the current global best and worst individuals found up till the k-th iteration. The diversity threshold is calculated as follows: Then, the NMA is enabled when . The FA would be capable of finding a relatively stable solution close to the target when e ξ −1 is from 0.0001 to 0.01, that is to say, it is feasible to take T 0 from 10 −3 to 10 −5 , and then the step parameter α would be small enough to correspond with the NMA.

Numerical Illustrative Examples
In order to demonstrate the superiority of the newly modified algorithm, 18 benchmark functions and a 12-story shear frame model are used to test the performance of the m-NMFA and to provide support for the comparison of the original FA, the GA [39] and the PSO [40] to the modified algorithm. Table 1 presents 18 well-known benchmark functions [41,42], and the order number (column 1), the function name (column 2), the formula (column 3), the feasible region (column 4), the minimum value (column 5), and the dimension (column 6) of each function are given successively. These functions can be generally divided into two categories: multimodal functions (f 1 -f 9 ) and unimodal functions (f 10 -f 18 ). The multimodal functions are used to prove the ability of the proposed algorithm to escape the local optimal solution and seek the global optimal solution, and the unimodal functions are used to illustrate the precision of the results searched out by the algorithm.

Benchmark Work of m-NMFA
In this test, all the population size and the number of chromosomes are taken to be 30. The maximum number of iterations is taken to be MaxGen = 1000. For the FA algorithms, the initial value of the step parameter is taken to be 0.5 in the m-NMFA and 0.25 in the original FA; the attractiveness parameters β min and β max take the values of 0.2 and 1.0, respectively; and the absorption coefficient γ is taken to be 1.0. For the GA, the test uses the built-in GA optimization algorithm of MATLAB R2017a Win64, and most of the input parameters take the default values of the built-in toolbox, except for the number of chromosomes, the maximum number of iterations, and FunctionTolerance which is defined as a stop criterion. When the average relative change in the best fitness function value over MaxGen generations is less than or equal to FunctionTolerance, the algorithm stops. Given that it demands the same iteration number for all the algorithms, the criterion value should be significantly small and, therefore, takes the value of 10 −100 in the test. For the PSO, the learning factors take the value of c 1 = c 2 = 2, the maximum velocity is taken to be 1.0, and the inertia weighting factor is in a linear form: where ω i is the inertia weighting factor of the i-th iteration, and ω 1 and ω 2 are the lower and upper limit values of ω. If the velocity or spatial vector of a particle in one dimension reaches or exceeds the boundary, all dimensions of the particle will be reset to the boundary value. All the functions are run independently using each of the optimization method for 100 times to avoid stochastic discrepancy in MATLAB R2017a Win64, Windows 10 system. Table 1 shows the results of the m-NMFA under different diversity threshold values, and additionally, the second to fifth columns show the results of the improved algorithm without the NMA. In each case, the mean value, the standard deviation (SD), and the maximum (Max) and minimum (Min) values of the results of each 100 runs are listed. It can be seen from the results that, compared to the case without the NMA, the other three cases are generally better, which proves that the introduction of the NMA can effectively improve the accuracy and stability of the firefly algorithm. In addition, although taking T 0 = 10 −1 can improve the optimization performance of the algorithm to a certain extent, the results are not good enough. When T = 10 −3 -10 −5 , the accuracy and stability of the optimization algorithm are better than the results of T = 10 −1 , which shows that the start point of the local optimization algorithm in the global optimization algorithm has a significant impact on the optimization result. This test also proves that a diversity threshold value in the range of T = 10 −3 ∼ 10 −5 is relatively appropriate. In subsequent tests and studies, the m-NMFA will adopt the value of T = 10 −5 . Table 2 shows the mean values, the standard deviations, and the maximum and minimum values of every 100 runs. It can be seen that the modified algorithm performs generally the same as the other three algorithms in f 1 and f 2 , and it is slightly inferior to the GA in f 7 and f 8. As for the minimum values of f 10 and f 14 , m-NMFA performs slightly worse than PSO and GA. Apart from the cases mentioned above, the m-NMFA significantly outperforms the FA, the GA, and the PSO in both precision and dispersion in most cases. It is worth mentioning that the ability of the m-NMFA to search for the minimum values of multimodal functions is well-illustrated in the tests of f 7 and f 8 for its distinguished minimum values out of the other values. where u(x i , 5, 100, 4),  Figure 2 shows the convergence diagrams of the 18 functions operated by the FA, the GA, the PSO, and the m-NMFA. The value of the ordinate is the mean value after 100 operations. It can be seen that the m-NMFA converges faster in most cases and has a strong high-precision solution searching ability.

Numerical Simulation of Shear Frame
A 12-story steel shear frame shown in Figure 3 is used to test the performance of the m-NMFA. The frame is assumed to be fixed on a solid ground and is of 1.8 m in height. The masses are uniformly distributed on each floor and are totally supported by six full-length rectangular columns with an identical cross-section of 6 mm in thickness and 25 mm in width. The mass of each floor is taken to be m e = 75 kg, and the inter-story stiffness is calculated to be k e = 1.512 × 10 5 N/m. Hence, the nominal undamaged structure mass and stiffness matrices are calculated as follows: Materials 2022, 15, x FOR PEER REVIEW 12 listed. It can be seen from the results that, compared to the case without the NMA other three cases are generally better, which proves that the introduction of the NMA effectively improve the accuracy and stability of the firefly algorithm. In addition hough taking T0 = 10 −1 can improve the optimization performance of the algorithm certain extent, the results are not good enough. When T=10 −3 -10 −5 , the accuracy and bility of the optimization algorithm are better than the results of T=10 −1 , which shows the start point of the local optimization algorithm in the global optimization algor has a significant impact on the optimization result. This test also proves that a dive threshold value in the range of T=10 -3~1 0 -5 is relatively appropriate. In subsequent and studies, the m-NMFA will adopt the value of T=10 −5 . Table A2 shows the mean values, the standard deviations, and the maximum minimum values of every 100 runs. It can be seen that the modified algorithm perfo generally the same as the other three algorithms in f1 and f2, and it is slightly inferi the GA in f7 and f8. As for the minimum values of f10 and f14, m-NMFA performs slig worse than PSO and GA. Apart from the cases mentioned above, the m-NMFA sig cantly outperforms the FA, the GA, and the PSO in both precision and dispersion in m cases. It is worth mentioning that the ability of the m-NMFA to search for the minim values of multimodal functions is well-illustrated in the tests of f7 and f8 for its di guished minimum values out of the other values.

Numerical Simulation of Shear Frame
A 12-story steel shear frame shown in Figure 3 is used to test the performance o m-NMFA. The frame is assumed to be fixed on a solid ground and is of 1.8m in he The masses are uniformly distributed on each floor and are totally supported by full-length rectangular columns with an identical cross-section of 6 mm in thickness 25 mm in width. The mass of each floor is taken to be me = 75 kg, and the inter-s stiffness is calculated to be ke = 1.512 × 10 5 N/m. Hence, the nominal undamaged struc mass and stiffness matrices are calculated as follows: Figure 2. Convergence diagram of four optimization algorithms.
It is assumed to have a stiffness reduction of overall 20%, 40%, and 20% on the 5th, 6th, and 7th inter-story columns, which is denoted by the artificially settled damage vector θ ith (i = 5, 6, 7). Gaussian white noises with mean values of zero are considered for the modal data, and the degrees of which are taken to be 1% for modal frequencies and 3% for modal shapes. Complete and incomplete measurements are considered. In the complete measurement case, all the first 12 modal parameters can be well measured. In the incomplete case, it is assumed that only the first eight modal parameters and the special displacement information of the 1st, 2nd, 4th, 6th, 8th, 10th, 11th, and 12th stories can be obtained. In order to eliminate accidental errors, the tests are executed independently for N t = 100 times and 100 sets of modal parameters are used. The simulation progress is operated using the FA, the GA, the PSO, and the m-NMFA separately and the algorithm parameters are kept the same as those in Section 4.1.
In this example, the standard deviation (SD) and the coefficient of variances (CV) are calculated for each parameter, which are listed in Tables 2 and 3. Figure 4 shows the results of the damage vector θ for both cases. It can be seen that, in the complete measurement case, all the algorithms are able to identify the position and the degree of the damage ideally. However, in the incomplete measurement case, only the m-NMFA can achieve the goal of identifying the damage of a structure.  It is assumed to have a stiffness reduction of overall 20%, 40%, and 20% on the 5th, 6th, and 7th inter-story columns, which is denoted by the artificially settled damage vector θ ith (i=5,6,7). Gaussian white noises with mean values of zero are considered for the modal data, and the degrees of which are taken to be 1% for modal frequencies and 3% for modal shapes. Complete and incomplete measurements are considered. In the complete measurement case, all the first 12 modal parameters can be well measured. In the incomplete case, it is assumed that only the first eight modal parameters and the special displacement information of the 1st, 2nd, 4th, 6th, 8th, 10th, 11th, and 12th stories can be obtained. In order to eliminate accidental errors, the tests are executed independently for = 100 times and 100 sets of modal parameters are used. The simulation progress is operated using the FA, the GA, the PSO, and the m-NMFA separately and the algorithm parameters are kept the same as those in Section 4.1.
In this example, the standard deviation (SD) and the coefficient of variances (CV) are calculated for each parameter, which are listed in Tables 2 and 3. Figure 4 shows the results of the damage vector θ for both cases. It can be seen that, in the complete measurement case, all the algorithms are able to identify the position and the degree of the damage ideally. However, in the incomplete measurement case, only the m-NMFA can achieve the goal of identifying the damage of a structure.

Experimental Illustrative Example
A six-story steel shear frame is designed to perform a vibration experiment, in which the modal parameters are identified by FDD to form a modal flexibility matrix. Then, the m-NMFA is utilized to search out the minimum value of the modal-flexibility-based objective function. Damage identification problem is introduced to illustrate the effectiveness of the proposed method.
The structure is fixed on the ground by foundation bolts. The length, width, and thickness of the interlayer plates are 450 mm, 450 mm, and 10 mm, respectively. The thickness center spacing of each plate is 250 mm. The interlayer stiffness is provided by the main columns fixed at the four vertices of the plates and additional columns at the midpoint of the side edges. An angle steel is welded with a steel plate in the center of it. An external excitation of random force and random position is applied to the structure by hammering excitation. The displacement response of each layer is measured using an IL-300 laser displacement meter. The measuring point is set at the center of the angle steel baffle. The sampling frequency of the meter is 200 Hz. The masses from m1 to m6 are 17.151 kg, 17.123 kg, 17.215 kg, 17.169 kg, 17.193 kg, and 16.731 kg, respectively, and the nominal interlaminar stiffness is 81920n/m. Before the test, a lumped mass reference finite element model is established according to the measured mass of each layer and the calculated nominal stiffness.
The displacement responses of the steel shear frames in the damaged and undamaged states are tested. In the damaged state, the additional columns of the fifth floor are missing, of which the nominal damage coefficient is −20%. In both tests, the responses of the 1st, 2nd, 4th, 5th, and 6th stories are measured, and the measurement repeats 50 times for each test and lasts for 5 min each. Figures 5 and 6 show the time domain information of the undamaged and damaged structures.

Experimental Illustrative Example
A six-story steel shear frame is designed to perform a vibration experiment, in which the modal parameters are identified by FDD to form a modal flexibility matrix. Then, the m-NMFA is utilized to search out the minimum value of the modal-flexibility-based objective function. Damage identification problem is introduced to illustrate the effectiveness of the proposed method.
The structure is fixed on the ground by foundation bolts. The length, width, and thickness of the interlayer plates are 450 mm, 450 mm, and 10 mm, respectively. The thickness center spacing of each plate is 250 mm. The interlayer stiffness is provided by the main columns fixed at the four vertices of the plates and additional columns at the midpoint of the side edges. An angle steel is welded with a steel plate in the center of it. An external excitation of random force and random position is applied to the structure by hammering excitation. The displacement response of each layer is measured using an IL-300 laser displacement meter. The measuring point is set at the center of the angle steel baffle. The sampling frequency of the meter is 200 Hz. The masses from m1 to m6 are 17.151 kg, 17.123 kg, 17.215 kg, 17.169 kg, 17.193 kg, and 16.731 kg, respectively, and the nominal interlaminar stiffness is 81920n/m. Before the test, a lumped mass reference finite element model is established according to the measured mass of each layer and the calculated nominal stiffness.
The displacement responses of the steel shear frames in the damaged and undamaged states are tested. In the damaged state, the additional columns of the fifth floor are missing, of which the nominal damage coefficient is −20%. In both tests, the responses of the 1st, 2nd, 4th, 5th, and 6th stories are measured, and the measurement repeats 50 times for each test and lasts for 5 min each. Figures 5 and 6 show the time domain information of the undamaged and damaged structures.
The frequencies of the first six modes of the frame structure are shown in Table 4, which take the average values of 50 measurements. The mode shapes shown in Figure 7 are normalized with respect to the reduced mass matrix M m .  The frequencies of the first six modes of the frame structure are shown in Table 4, which take the average values of 50 measurements. The mode shapes shown in Figure 7 are normalized with respect to the reduced mass matrix .  The frequencies of the first six modes of the frame structure are shown in Table 4, which take the average values of 50 measurements. The mode shapes shown in Figure 7 are normalized with respect to the reduced mass matrix .  In the Bayesian probabilistic approach, the probability distributions of the model updating parameters of the damaged and undamaged structures are needed to calculate the cumulative probability distribution of the damage parameters in the damaged state. The distributions of the damage parameters can be obtained using the following formula: In the Bayesian probabilistic approach, the probability distributions of the model updating parameters of the damaged and undamaged structures are needed to calculate the cumulative probability distribution of the damage parameters in the damaged state. The distributions of the damage parameters can be obtained using the following formula: where the damage extent d is used to describe the value of the stiffness reduction parameter; θ d n and θ ud n are the model updating parameters of the damaged and undamaged structures; andθ d n ,θ ud n ,σ d n andσ ud n are the mean values and standard deviations of the model updating parameters in Gaussian probability distribution. The value of θ n can be obtained by choosing the mean value of each distribution in the probability density function diagram. The standard deviations of the n-th damage parameters can be calculated using the following formula: where ∆d = d n,v+1 −d n,v , v ∈ {i, j|i ∈ N * , j ∈ N * }, and y(d) is the Gaussian probability density function of d.
The model updating parameters of the finite element model of the undamaged structure are identified using the m-NMFA, which are shown in Table 5. The model updating parameters of the damaged structure are shown in Table 6. Table 5. Model updating parameters of the undamaged structure using the Bayesian probability method.  Table 6. Model updating parameters of the damaged structure using the Bayesian probability method.  Tables 5 and 6, the cumulative probability distribution of the structural damage parameter of each floor is calculated using Formula 35. The cumulative probability distribution is plotted in Figure 8. The probability density function is plotted in Figure 9. Figure 10 shows the mean value of each distribution in the probability density function directly. Formula 35. The cumulative probability distribution is plotted in Figure 8. The probability density function is plotted in Figure 9. Figure 10 shows the mean value of each distribution in the probability density function directly.

Conclusions
Based on the modal flexibility matrix, this paper uses the Bayesian probability method to establish the objective function for the model updating problem, which is solved using the proposed m-NMFA. The results show that this method can accurately fix the damage position and determine the damage extent.
Compared to the other meta-heuristic optimization algorithms, the proposed m-NMFA has a stronger optimization ability in solving benchmark functions and the finite element model updating objective functions, which shows a higher applicability in the practical problems of model updating. The reason why the m-NMFA has more superior optimization performance is mainly due to the following three points: 1. In the research on the random step size formula, it is found that the smaller the random step size is at the beginning of the iteration, the faster the convergence speed of the algorithm will be. However, as the random step size gradually approaches zero at the end of the iteration, its value should be large enough to keep the algorithm with sufficient exploration ability. 2. In the research on the value of diversity threshold, a fraction is used to quantify the diversity of the population to enable the NMA algorithm to start at an appropriate time. The diversity threshold is quantized into a value between (0, e -1 ] through exponential form. When the diversity threshold is taken to be between 10 −3~1 0 −5 , the algorithm can obtain a more accurate optimal solution. 3. The selection of the iterative formula of the attraction parameter has great impact on the solving ability of the FA in solving the multi-dimensional optimization problems, The selection of the iterative formula of the attraction parameter has a great impact

Conclusions
Based on the modal flexibility matrix, this paper uses the Bayesian probability method to establish the objective function for the model updating problem, which is solved using the proposed m-NMFA. The results show that this method can accurately fix the damage position and determine the damage extent.
Compared to the other meta-heuristic optimization algorithms, the proposed m-NMFA has a stronger optimization ability in solving benchmark functions and the finite element model updating objective functions, which shows a higher applicability in the practical problems of model updating. The reason why the m-NMFA has more superior optimization performance is mainly due to the following three points:

1.
In the research on the random step size formula, it is found that the smaller the random step size is at the beginning of the iteration, the faster the convergence speed of the algorithm will be. However, as the random step size gradually approaches zero at the end of the iteration, its value should be large enough to keep the algorithm with sufficient exploration ability.

2.
In the research on the value of diversity threshold, a fraction is used to quantify the diversity of the population to enable the NMA algorithm to start at an appropriate time. The diversity threshold is quantized into a value between (0, e −1 ] through exponential form. When the diversity threshold is taken to be between 10 −3~1 0 −5 , the algorithm can obtain a more accurate optimal solution. 3.
The selection of the iterative formula of the attraction parameter has great impact on the solving ability of the FA in solving the multi-dimensional optimization problems, The selection of the iterative formula of the attraction parameter has a great impact on the ability of the FA in solving multi-dimensional optimization problems, which would lead to stagnation and non-convergence if an improper formula were selected. This paper avoids such problem by selecting an appropriate formula for the coefficient of attraction.
In the research of model updating based on modal flexibility and Bayesian method, the modal flexibility theory, the system equivalent reduction and expansion process, and the basic theory of Bayesian method are introduced. The objective function established based on these theories provides a relatively complete theoretical support for solving the model updating problems of limited measuring points and number of modes, so that people can obtain reliable results from limited structural information. Bayesian method allows people to empirically select the prior probability distribution of the modified parameters so as to exert influence on the identification results, which provides theoretical and practical convenience for solving the problems of model updating. Through the numerical simulation of a 12-story shear frame and the experimental test of a 6-story steel shear frame, the reliability of the Bayesian probability method in solving model updating and structural damage identification is verified.