Research on Multi-Alternatives Problem of Finite Element Model Updating Based on IAFSA and Kriging Model

Due to insufficient test data, insufficient constraint equations and uncertain objective function, the local optimal solution and the global optimal solution of the objective function in finite element model updating may represent the actual parameters of the structure. Based on this, this paper proposes an improved artificial fish school algorithm. By combining the niche technology with the artificial fish school algorithm, the improved algorithm can systematically find multiple global optimal solutions and local optimal solutions of the objective function. Aiming at the difficulty of determining the niche radius, an adaptive niche radius mechanism is proposed. The improved algorithm is used to study the multi-alternatives problem of finite element model updating after verifying its feasibility through numerical simulation analysis. In the case of benchmark framework model updating, it is confirmed that multi-alternative problems exist and the global optimal solution of the objective function does not necessarily represent the true parameters of the structure. In case 2, the improved algorithm combined with the Kriging model is applied to the model updating of a cable-stayed footbridge, and 15 sets of solutions are obtained, in which the error objective function values of the measured and theoretical values of the bridge modes are close but the solutions are completely different. Combining with the actual bridge condition and reanalysis technology, the author takes the suboptimal solution 2 as the most representative solution of the bridge parameters, which reduces the possibility of misjudgment of structural parameters.


Introduction
The finite element model is usually established according to the design drawings. In the process of modeling, there are geometric parameters, physical parameters, boundary conditions, and other assumptions, which make the theoretical values of the model have errors with the actual responses of the structure. The accurate finite element model is very important for the evaluation of the seismic and wind resistance performance, structural damage identification, and health condition monitoring of the existing structure [1][2][3]. In order to obtain the benchmark model, the finite element model updating (FEMU) aimed at reducing the structural theoretical responses and the actual responses error has been extensively studied by scholars at home and abroad. The research review and results of FEMU can be found in the literature [4][5][6][7][8][9][10][11][12][13][14]. In summing up the research progress of FEMU, Li et al. [7] reviewed the main methods of traditional FEMU and pointed out that the FEMU is developing from linear to nonlinear. Teughels et al. [12] carried out the FEMU research earlier; he proposed a method named coupled local minimizers, which was applicable to global optimization of functions In order to obtain the real parameters of the actual FEM of the structure, the error between the measured value of the structure and the theoretical value of the FEM is usually reduced to the greatest extent. The expression of the objective function is shown in Equation (2).
where r i (x) means the residual function that expresses the relative difference between the analytically predicted response and the experimentally measured response, α i represents the weight value of each response. The response can be modal frequency, nodal displacement, nodal acceleration, etc., represented only by reducing the error between the measured modal frequency of the structure and the theoretical frequency of FEM; the expression is shown in Equation ( where x is the value of the parameter to be modified, VLB and VUB are the upper and lower limits of parameter values, respectively, f j m (x) represents the j-th frequency calculated by FEM, f j e means the j-th frequency measured by the experiment, and m is the modal order.
The diversity of the measured response information and the uncertainty of updating parameters determine that the objective function constructed in the FEMU is uncertain, and the optimal parameters obtained by different objective functions may not be the same. The schematic diagram of multi-alternatives caused by different response information is shown in Figure 1. ( ) m f x represents the j-th frequency calculated by FEM, e f means the j-th frequency measured by the experiment, and m is the modal order.
The diversity of the measured response information and the uncertainty of updating parameters determine that the objective function constructed in the FEMU is uncertain, and the optimal parameters obtained by different objective functions may not be the same. The schematic diagram of multi-alternatives caused by different response information is shown in Figure 1. The following reasons may also lead to the multiple solutions problem of finite element model updating: The following reasons may also lead to the multiple solutions problem of finite element model updating: (1) In the process of structural test, an insufficient number of sensors leads to incomplete test data, and the objective function in finite element model updating may correspond to multiple solutions due to insufficient constraint equations; (2) There may be errors in the experimental data. When the objective function is constructed by the experimental data with errors, the global optimal solution may be a set of deceptive solutions. On the contrary, the local optimal solution of the objective function may best represent the real parameters of the structure; (3) The finite element model updating is an inverse problem in mathematics. Due to the complexity of the finite element model, one response may correspond to multiple solutions.
In fact, even if the objective function and updating parameters are determined, the same objective function value may correspond to multiple solutions due to the complexity of the FEM. Taking the two updating parameters as an example, the objective function f simultaneously has two global optimal solutions (x 1 , y 1 ) and (x 2 , y 2 ), and simultaneously satisfies Equation (4), as shown in Figure 2. For a local optimal solution (x 3 , y 3 ) that exists, if Equation (5) is satisfied, it will be used as an alternative. In the formula, ε is the threshold value, which controls whether to accept the suboptimal solution.
f (x 1 , y 1 ) = f (x 2 , y 2 ) = min( f (x, y)) (4) f (x 3 , y 3 ) − min( f (x, y)) = ε Sensors 2020, 20, x FOR PEER REVIEW 4 of 24 (1) In the process of structural test, an insufficient number of sensors leads to incomplete test data, and the objective function in finite element model updating may correspond to multiple solutions due to insufficient constraint equations; (2) There may be errors in the experimental data. When the objective function is constructed by the experimental data with errors, the global optimal solution may be a set of deceptive solutions. On the contrary, the local optimal solution of the objective function may best represent the real parameters of the structure; (3) The finite element model updating is an inverse problem in mathematics. Due to the complexity of the finite element model, one response may correspond to multiple solutions.
In fact, even if the objective function and updating parameters are determined, the same objective function value may correspond to multiple solutions due to the complexity of the FEM. Taking the two updating parameters as an example, the objective function f simultaneously has two global optimal solutions (x1, y1) and (x2, y2), and simultaneously satisfies Equation (4), as shown in Figure 2. For a local optimal solution (x3, y3) that exists, if Equation (5) is satisfied, it will be used as an alternative. In the formula, ε is the threshold value, which controls whether to accept the suboptimal solution.  Therefore, it is not reasonable for the traditional FEMU technology to take a group of global optimal solutions as the most representative structural actual parameters, so it is urgent to find a multi-alternatives technology in the FEMU. Through this technique, we can find multiple global optimal solutions or local optimal solutions which can reduce the error between the measured value and the theoretical value. The corresponding objective function values of multiple solutions are similar, but the solution space is quite different. The decision-makers can select one or more groups based on the comprehensive analysis of engineering experience, which can not only effectively avoid the misjudgment of structural parameters, but also carry out the risk assessment and reliability calculation of the structure by combining multi-alternatives.  Therefore, it is not reasonable for the traditional FEMU technology to take a group of global optimal solutions as the most representative structural actual parameters, so it is urgent to find a multi-alternatives technology in the FEMU. Through this technique, we can find multiple global optimal solutions or local optimal solutions which can reduce the error between the measured value Sensors 2020, 20, 4274 5 of 25 and the theoretical value. The corresponding objective function values of multiple solutions are similar, but the solution space is quite different. The decision-makers can select one or more groups based on the comprehensive analysis of engineering experience, which can not only effectively avoid the misjudgment of structural parameters, but also carry out the risk assessment and reliability calculation of the structure by combining multi-alternatives.

Kriging Model
Kriging model is composed of two parts, one is the deterministic polynomial regression model, which can well fit the linear and low-order nonlinear model problems, the other is the stochastic refinement simulation process using Gaussian function to fit the high-order nonlinear model. The expression of the Kriging model is shown in Equation (6).
where x and y(x) represent input and output responses respectively, f T (x)β is polynomial regression model, β indicates the coefficient of polynomial vector regression, and z(x) refines the simulation part randomly with Gaussian function.
To construct the Kriging model, we need to obtain test samples, usually adopting the design methods of central composite design, Box-Behnken design, and Latin hypercube design (LHD) [31]. After the Kriging model is built, its fitting accuracy needs to be tested. Generally, the determination coefficient R 2 , root mean square error (RMSE), mean square error (MSE), and other methods are used to test, and it can only be used after meeting the accuracy requirements.

Procedure of Multi-Alternatives of FEMU
This section summarizes the detailed procedure of Kriging model assisted multi-alternatives FEMU. The whole approach can be divided into the following six steps: Step 1: An initial FEM using commercial FE software is established based on design documents; Step 2: Determine the updating parameters by combining engineering experience and parameter sensitivity analysis or parameter significance analysis; Step 3: Use Latin hypercube design to generate samples, and the samples are substituted into the finite element model to extract the corresponding responses; Step 4: The Kriging model to describe the relationship between structural responses and updating parameters is established and test the fitting accuracy of the proxy model. If the accuracy of the proxy model meets the requirements, go to the next step, otherwise return to step 3; Step 5: Construct the error objective function between the responses of the Kriging model and the measured values of the structure; Step 6: Use the improved artificial fish swarm algorithm proposed in this paper to optimize the parameters and gets the objective function's multiple solutions; Step 7: Evaluate the results of FEMU and comprehensively analyze the multiple solutions, trying to select one or more solutions that best represent the structural damage parameters.
The flow chart is shown in Figure 3.
the measured values of the structure; Step 6: Use the improved artificial fish swarm algorithm proposed in this paper to optimize the parameters and gets the objective function's multiple solutions; Step7: Evaluate the results of FEMU and comprehensively analyze the multiple solutions, trying to select one or more solutions that best represent the structural damage parameters.
The flow chart is shown in Figure 3.

Standard Artificial Fish Swarm Algorithm
Artificial fish swarm algorithm (AFSA) [32] is one of many swarm intelligence algorithms, which simulates the behavior of a fish swarm always moving in the direction of high food concentration and low crowding degree and has the advantages of strong robustness and global search ability. The principles of AFSA are as follows: (1) Foraging behavior. Fish search other fish in the field of vision within the maximum number of attempts until they find the fish that is better than the current situation and move in this direction. The mathematical expression is shown in Equation (7): (2) Swarming behavior. The fish find the central coordinate in the field of vision. If the food concentration is better than the current situation and not crowded, they will swim to the central coordinate. The mathematical expression is shown in Equation (8): (3) Following behavior. The fish move towards the fish with the highest food concentration in the field of vision. The mathematical expression is shown in Equation (9): Sensors 2020, 20, 4274 7 of 25 (4) Random behavior. When the artificial fish does not satisfy the above three behavioral conditions, it will perform a random swimming behavior, the expression of which is shown in Equation (10): where x i (next) is the next position of the fish, x i is the current position, x j is the fish better than the current position, x c is the coordinate center of the fish in the field of vision, and x best is the optimal coordinate in the field of vision.
Step means the moving step, and visual indicates the field of vision.

Improved Artificial Fish Swarm Algorithm
AFSA has been successfully applied in aerospace, civil engineering, communication engineering, and other fields [33][34][35][36]. For example, Li et al. [33] applied AFSA to FEMU and structural damage detection. The results show that the structural damage detection method based on AFSA can effectively modify FEM of the structure and accurately identify the damage degree under different noise levels and damage conditions. Zhang et al. [34] introduced crossover operator and Gaussian mutation operator in genetic algorithm into AFSA, which accelerated the convergence speed of AFSA and improved its optimization accuracy. Then, the improved AFSA was applied to FEMU of an aircraft model, and good results were obtained. Zhu et al. [35] found that the convergence of AFSA will be slow or even stagnant in the process of optimization. Therefore, an improved method using U-shaped function transformation to deal with the random swimming step of artificial fish is proposed, which also improves the optimization accuracy and convergence speed of AFSA. Then, the improved AFSA is applied to the model updating of a cable-stayed bridge, which makes the theoretical responses of the structure match the measured responses.
The above improvement measures can either accelerate the convergence speed or improve the global optimization ability of AFSA. However, the improved algorithm can only find a global optimal solution of multi-peak functions, which is not applicable to the multi alternatives problem of FEMU. Therefore, the standard AFSA still has the following defects: (1) difficult in finding multiple global optimal solutions and suboptimal solutions for multi-peak functions; (2) when AFSA runs to later stage, the artificial fish oscillates near the global optimal solution, and the optimization accuracy is low. In view of the above two problems, the authors propose an improved artificial fish swarm algorithm (IAFSA). The improvement measures are as follows: (1) In view of the difficulty of finding multiple peaks in the multi-peak function of AFSA, the sequential niche technology is combined with AFSA, and the niche radius R is introduced to control the distance between the fish. The specific steps are as follows: Step 1: Execute AFSA once, find a global optimal solution, and liberate it into the niche core; Step 2: Repeat AFSA, calculate the distance between the current fish and the existing solution in the niche core, and compare the distance with the radius of each niche, if within the radius, initialize the artificial fish position; Step 3: Calculate the distance between the optimal solution obtained by each AFSA and the solutions in the niche core. If it is outside the niche radius, put the new solution set in the niche core, otherwise give it up; Step 4: Continue the next cycle until all local and global optimal solutions are found or the algorithm termination conditions are satisfied.
(2) To solve the problem of difficulty in determining the niche radius, an adaptive niche radius mechanism is proposed. The steps to calculate the niche radius of the multi-peaks function f (x 1 x 2 x 3 . . . x d ) are as follows: Step 1: When there is a solution X (x 1 x 2 x 3 . . . x d ) in the niche core, look for the niche radius in the increasing direction of the i-th dimension. Increase ∆R i along the x i direction, and calculate f (x 1 x 2 Sensors 2020, 20, 4274 8 of 25 x i + ∆R i . . . x d ), ∆R i is set according to the definition domain of each dimension parameter; ∆R i is recommended as shown in Equation (11): where x up i and x lb i are the upper and lower limits of definition field of the i-th dimension parameter.
The niche radius in the increasing direction of the i-th dimension parameter is n∆R i ; Step 3: Look for the niche radius in the decreasing direction of the i-th dimension. Decrease ∆R i along the x i direction and calculate f ( The niche radius in the decreasing direction of the i-th dimension is m∆R i ; Step 4: Available from step 1 to step 3, the niche radius of the i-th dimension is R i = min(n∆R i , m∆R i ); Step 5: Determine the niche radius of other dimensions, then the radius of each dimension of X ( (3) In order to improve the accuracy of AFSA optimization, the optimal solution obtained from each execution of AFSA is combined with simulated annealing algorithm to carry out local refinement and optimization. The specific steps are as follows: Step 1: Set the initial temperature T 0 and the end temperature Tend of the simulated annealing algorithm to determine the number of iterations L at each temperature; Step 2: Execute the steps of AFSA to find the best solution x best ; Step 3: Let x best have a random disturbance and get a new solution x best '; Step 4: Judge whether f (x best ') is better than f (x best ), if so, accept it; otherwise, judge whether exp(−df /T) > rand is satisfied, if yes, still accept x best ', otherwise give up it. Repeat step 3 and step 4 until the termination temperature condition is met where df = f (x best ') − f (x best ); rand means a random number between 0 and 1.

Numerical Simulation Analysis
Three classic multi-peak functions are used to test IAFSA. The function expression is shown in Equations (12)- (14).
Features: f 1 (x, y) contains 100 peaks, and the peaks gradually increase from the coordinates (0,0) to (1,1). It is used to test the performance of IAFSA in finding multiple peaks of multi-peak function.
Features: f 2 (x, y) contains 36 peaks with unequal spacing. The performance of adaptive niche radius mechanism of IAFSA is tested by this function.
Features: f 3 (x, y) contains four global optimal solutions (GOS) with theoretical values of 200. The accuracy of IAFSA is tested by this function.
Before testing, the basic parameters of IAFSA need to be set, the basic parameters are set according to the following principles: (1) Fish population size. The size of fish population ranges from 50 to 200, if the value is too large, it will increase the time consumption; if the value is too small, the optimization accuracy will be low; (2) Number of iterations. Similar to the fish population size, the recommended range is 20 to 100; (3) Visual and step. These two parameters are not sensitive to the optimization effect. Their values are suggested in Equations (15)- (16).
Step = ( 1 100 (4) Crowding factor. This parameter can range from 0.1 to 0.3; (5) The recommended value of ∆R is shown in Equation (11); (6) Number of AFSA executions. This value indicates the number of times AFSA is executed when the IAFSA is executed once; the value is greater than the peak numbers of the test function.
According to the above principles, the basic parameter settings of IAFSA are shown in Table 1. Perform IAFSA 50 times to find the peaks of f 1 . In order to compare the performance of SSGA proposed previously [17] and IAFSA, SSGA is also executed 50 times. Limited to space, extract the first five search results. Meanwhile, five groups of AFSA will be executed, and each group will execute AFSA 120 times (ensure that the number of AFSA executions in each group is the same as the number of IAFSA executions each time). The results are shown in Table 2.  Record the peak seeking process performed by IAFSA once, and when IAFSA performs AFSA 30 times, 60 times, 90 times, and 120 times. The results are shown in Figures 4-7.    Record the peak seeking process performed by IAFSA once, and when IAFSA performs AFSA 30 times, 60 times, 90 times, and 120 times. The results are shown in Figures 4-7.    From Table 2 and Figures 4-7, we can see that the first five times of IAFSA's execution   From Table 2 and Figures 4-7, we can see that the first five times of IAFSA's execution successfully found one global optimal solution and 99 local optimal solutions of f1. In fact, the author counted the results of running IAFSA 50 times and found 100 peaks of f1 successfully each time. However, SSGA can only find one global optimal solution and a few local optimal solutions of f1, and the number of solutions found is not stable. Meanwhile, the result of running AFSA for each group shows that AFSA can find only one global optimal solution of f1, and 99 local optimal solutions are missed.

LOS Found
From Table 2 and Figures 4-7, we can see that the first five times of IAFSA's execution successfully found one global optimal solution and 99 local optimal solutions of f 1 . In fact, the author counted the results of running IAFSA 50 times and found 100 peaks of f 1 successfully each time. However, SSGA can only find one global optimal solution and a few local optimal solutions of f 1 , and the number of solutions found is not stable. Meanwhile, the result of running AFSA for each group shows that AFSA can find only one global optimal solution of f 1 , and 99 local optimal solutions are missed.
In order to test the adaptive niche radius mechanism of IAFSA, IAFSA and fixed niche radius R = 0.1, R = 0.3, and R = 0.5 were executed 50 times to find f 2 peak respectively, and the results of the first 10 times are shown in Table 3. Perform IAFSA once and find the peak of f 2 as shown in Figure 8. Table 3. Results of optimizing f 2 with different niche radius. 1  1  35  1  6  1  31  1  23  2  1  35  1  6  1  31  1  21  3  1  35  1  10  1  30  1  21  4  1  35  1  9  1  31  1  23  5  1  35  1  12  1  29  1  22  6  1  35  1  11  1  31  1  24  7  1  35  1  10  1  31  1  21  8  1  35  1  9  1  31  1  21  9  1  35  1  12  1  29  1  23  10  1  35  1  7  1  31  1  24 Sensors 2020, 20, x FOR PEER REVIEW 11 of 24 It can be seen from Table 3 that IAFSA using adaptive niche radius mechanism successfully finds one global optimal solution and 35 local optimal solutions of function f2. In fact, running IAFSA 50 times successfully found 36 peaks of function f2 each time. When the niche radius is fixed, the niche artificial fish swarm algorithm can find the global optimal solution of function f2, but it can only find It can be seen from Table 3 that IAFSA using adaptive niche radius mechanism successfully finds one global optimal solution and 35 local optimal solutions of function f 2 . In fact, running IAFSA 50 times successfully found 36 peaks of function f 2 each time. When the niche radius is fixed, the niche artificial fish swarm algorithm can find the global optimal solution of function f 2 , but it can only find a small number of local optimal solutions, and the number of solutions is not stable.

LOS Found
In order to test the optimization accuracy of IAFSA, IAFSA was run five times, AFSA was run in five groups, and each group was run five times. The statistical results are shown in Table 4, one of the results of running IAFSA is shown in Figure 9.  It can be seen from Table 3 that IAFSA using adaptive niche radius mechanism successfully finds one global optimal solution and 35 local optimal solutions of function f2. In fact, running IAFSA 50 times successfully found 36 peaks of function f2 each time. When the niche radius is fixed, the niche artificial fish swarm algorithm can find the global optimal solution of function f2, but it can only find a small number of local optimal solutions, and the number of solutions is not stable.
In order to test the optimization accuracy of IAFSA, IAFSA was run five times, AFSA was run in five groups, and each group was run five times. The statistical results are shown in Table 4, one of the results of running IAFSA is shown in Figure 9.  It can be seen from Table 3  -1000 Figure 9. Results of performing IAFSA once to find the peaks of f 3.
It can be seen from Table 3 and Figure 9 that IAFSA runs five times and finds four sets of global optimal solutions for f 3  The results show that running AFSA for many times cannot guarantee the convergence to all the global optimal solutions of the function, and there is the problem of low accuracy in the optimization.

Introduction of ASCE Benchmark Model
ASCE benchmark model was proposed by Black and Ventura and built in 1998. The purpose of model construction is to be used for structural health monitoring and damage identification, providing a unified test platform for different damage identification methods. The model is a 4-story × 2-span × 2-span frame structure composed of beams, columns, braces, and other components. The size and material characteristics of each component are detailed elsewhere [17], and the model is shown in Figure 10.
ASCE benchmark model was proposed by Black and Ventura and built in 1998. The purpose of model construction is to be used for structural health monitoring and damage identification, providing a unified test platform for different damage identification methods. The model is a 4-story × 2-span × 2-span frame structure composed of beams, columns, braces, and other components. The size and material characteristics of each component are detailed elsewhere [17], and the model is shown in Figure 10. ASCE benchmark model can select 60 degrees of freedom or 120 degrees of freedom mode by modifying the MATLAB code, and also can customize the damage of each member. In this paper, 120 degrees of freedom models are selected. We set the damage to the Young's modulus of the first layer diagonal brace of the model, and the damaged part of the member is shown in Figure 11. The Young's modulus values of the damaged members are E1 = 0.75E0 (north side), E2 = 0.5E0 (east side), E3 = 0.45E0 (south side), E4 = 0.95E0 (west side). Calculate the first six frequencies of the model before and after the damage, as shown in Table 5. ASCE benchmark model can select 60 degrees of freedom or 120 degrees of freedom mode by modifying the MATLAB code, and also can customize the damage of each member. In this paper, 120 degrees of freedom models are selected. We set the damage to the Young's modulus of the first layer diagonal brace of the model, and the damaged part of the member is shown in Figure 11. The Young's modulus values of the damaged members are E 1 = 0.75E 0 (north side), E 2 = 0.5E 0 (east side), E 3 = 0.45E 0 (south side), E 4 = 0.95E 0 (west side). Calculate the first six frequencies of the model before and after the damage, as shown in Table 5.  It can be seen from Table 5 that the maximum relative error of the first six frequencies of the initial model and the damaged model reach 7.12%. The modal assurance criterion (MAC) values are all above 0.947, which show that the first six frequencies of the damage model and the initial model are in the same order mode with high confidence.

Establishment of Kriging Model
Next, 44 sets of sample data were generated using the Latin hypercube design method, and Kriging models with E 1 , E 2 , E 3 , E 4 as input parameters and f 1 -f 6 as output were established. To ensure the physical significance of the four updating parameters, when generating sample data, set the parameter range to [40% E 0 , E 0 ], where E 0 represents the initial Young's modulus of the member. To facilitate data processing, input parameters are normalized. The Kriging model is established based on the sampled data. Figure 12a shows the first frequency in the function of E 1 and E 2 , and Figure 12b further shows the MSE value of the Kriging model. It can be seen from Table 5 that the maximum relative error of the first six frequencies of the initial model and the damaged model reach 7.12%. The modal assurance criterion (MAC) values are all above 0.947, which show that the first six frequencies of the damage model and the initial model are in the same order mode with high confidence.

Establishment of Kriging Model
Next, 44 sets of sample data were generated using the Latin hypercube design method, and Kriging models with E1, E2, E3, E4 as input parameters and f1-f6 as output were established. To ensure the physical significance of the four updating parameters, when generating sample data, set the parameter range to [40% E0, E0], where E0 represents the initial Young's modulus of the member. To facilitate data processing, input parameters are normalized. The Kriging model is established based on the sampled data. Figure 12a shows the first frequency in the function of E1 and E2, and Figure 12b further shows the MSE value of the Kriging model. It can be seen from Figure 12 that almost all the sample points fall on the response surface, and the magnitude of the mean square error of the parameters E1 and E2 corresponding to f1 is 10 −11 , which indicates that the Kriging model constructed has high precision. In order to further test the fitting accuracy of f1-f6, 50 groups of sample data are regenerated for Kriging model prediction, and the accuracy is tested by using the determination coefficient (R 2 ) and root mean square error (RMSE). The test results are shown in Table 6   Table 6. Determination coefficient (R 2 ) and root mean square error (RMSE) values corresponding to the first six frequencies.  It can be seen from Figure 12 that almost all the sample points fall on the response surface, and the magnitude of the mean square error of the parameters E 1 and E 2 corresponding to f 1 is 10 −11 , which indicates that the Kriging model constructed has high precision. In order to further test the fitting accuracy of f 1 -f 6 , 50 groups of sample data are regenerated for Kriging model prediction, and the accuracy is tested by using the determination coefficient (R 2 ) and root mean square error (RMSE). The test results are shown in Table 6. It can be seen from Table 6 that R 2 values are above 0.9988, and RMSE values are all less than 3.99 × 10 −5 , which further indicates that the constructed Kriging model has a high fitting accuracy.

Model Updating and Results Analysis
The updating parameters are E 1 -E 4 , and the error objective function of fitting value of Kriging model and theoretical value of the damage model is shown in Equation (17).
where f i Kri (x) represents the i-th frequency predicted by Kriging model, f i t (x) represents the i-th frequency of damage model, and x means the updating parameters.
Running IAFSA, one group of global optimal solutions and three groups of local optimal solutions of the objective function are obtained. The ratio of each parameter after updating to before updating is shown in Table 7. Table 7. Comparison of parameters before and after updating. Four sets of solutions are substituted into ASCE benchmark model to calculate f 1 -f 6 , and the relative error between the theoretical response corresponding to each set of solutions and the damage model response is calculated as shown in Table 8. It is known from Table 8 that the four sets of solutions all effectively reduce the error between the updated model and the damaged model. From the perspective of the degree of error reduction, it is clear that the global optimal solution is the best recommended solution. However, the damage of the ASCE benchmark model is preset by the authors. The values of E 1 , E 2 , E 3 , and E 4 are 0.750 E 0 , 0.50 E 0 , 0.45 E 0 , and 0.95 E 0 , respectively. The global optimal solutions obtained by IAFSA are 0.750 E 0 , 0.950 E 0 , 0.450 E 0 , and 0.499 E 0 , respectively. Obviously, the global optimal solution is not the damage value preset by the authors. On the contrary, the closest solution to the preset damage value is local optimal solution 1, and the corresponding parameters are 0.751 E 0 , 0.497 E 0 , 0.451 E 0 , and 0.954 E 0 . In order to clearly compare the relationship between each group solution and the preset damage value, the ratio of the parameters corresponding to each group of solutions to the initial value is shown in Figure 13. The ordinate in the graph represents the ratio of the corresponding parameters of each group of solutions to the initial value, and the abscissa represents the corresponding parameters of each group of solutions. Sensors 2020, 20, x FOR PEER REVIEW 15 of 24 Figure 13. Ratio of parameters corresponding to each group of solutions to initial values.
Through the case of ASCE benchmark, the authors not only prove that IAFSA can find multiple solutions of objective function in FEMU, but also prove that the global optimal solution in FEMU does not necessarily represent the real parameters of the structure.

Engineering Background
The pedestrian bridge in this case is a steel box girder cable-stayed bridge with a total length of 128.8 m, of which the main span is 45 m and the side span is 21 m. The main tower adopts a rectangular steel box girder with a variable cross-section and a height of 24.5 m. The main tower is provided with four parallel cables on both sides along the longitudinal bridge direction, the cable diameter is 115 mm, the design tension of the cables is: fN4, fB1 − fB4 = 430 kN, fN1 − fN3 = 355 kN. The elevation of the pedestrian overpass is shown in Figure 14. We used Ansys15.0 software to establish the FEM of this bridge, as shown in Figure 15. Through the case of ASCE benchmark, the authors not only prove that IAFSA can find multiple solutions of objective function in FEMU, but also prove that the global optimal solution in FEMU does not necessarily represent the real parameters of the structure.

Engineering Background
The pedestrian bridge in this case is a steel box girder cable-stayed bridge with a total length of 128.8 m, of which the main span is 45 m and the side span is 21 m. The main tower adopts a rectangular steel box girder with a variable cross-section and a height of 24.5 m. The main tower is provided with four parallel cables on both sides along the longitudinal bridge direction, the cable diameter is 115 mm, the design tension of the cables is: f N4 , f B1 − f B4 = 430 kN, f N1 − f N3 = 355 kN. The elevation of the pedestrian overpass is shown in Figure 14. In order to obtain the dynamic response, the environmental vibration of the bridge was tested. A total of 53 measuring points are set in the whole bridge, which are divided into 17 measurement batches. In order to ensure the comparability and consistency of the measured data of each batch, each batch includes two reference points, namely point 8 and point 15. The data acquisition time of each batch is 15 min, and the acquisition frequency is set to 200 Hz. The layout of dynamic test points of the footbridge is shown in Figure 16. The modal information of the bridge was identified by the stochastic subspace method, and the order was determined by the stability diagram. Take the identified frequency as shown in Figure 17, at the same time, extract the corresponding modal response under the initial modeling parameters of the finite element model as shown in Table 9. In order to obtain the dynamic response, the environmental vibration of the bridge was tested. A total of 53 measuring points are set in the whole bridge, which are divided into 17 measurement batches. In order to ensure the comparability and consistency of the measured data of each batch, each batch includes two reference points, namely point 8 and point 15. The data acquisition time of each batch is 15 min, and the acquisition frequency is set to 200 Hz. The layout of dynamic test points of the footbridge is shown in Figure 16. The modal information of the bridge was identified by the stochastic subspace method, and the order was determined by the stability diagram. Take the identified frequency as shown in Figure 17, at the same time, extract the corresponding modal response under the initial modeling parameters of the finite element model as shown in Table 9. Sensors 2020, 20, x FOR PEER REVIEW 17 of 24 Figure 17. Stability graph of stochastic subspace identification.

Updating Parameters Selection and Establishment of Kriging Model
When selecting the updating parameters, the following aspects are mainly considered: (1) The steel box girder is a prefabricated component, so it is difficult to change the section area or the moment of inertia due to the size deviation of the section; (2) The tension of cable may relax under frequent loads; (3) There are different degrees of corrosion in straight-line segment and U-shape segment of steel box girder, which greatly affects the Young's modulus and density of beam; (4) The actual stiffness of the main tower should be greater than the design value when the section thickness of the main tower at the connection of the cable and the main tower increases; (5) T-shaped pier as a member supporting the main beam; its Young's modulus and density may change with the increase of service life.
Therefore, the following parameters are preliminarily selected as the updating parameters: density and Young's modulus of main tower (denoted by D1 and E1), density and Young's modulus of steel box girder in straight-line segment (denoted by D2 and E2), density and Young's modulus of steel box girder in U-shape segment (denoted by D3 and E3), density and Young's modulus of the cables (denoted by D4 and E4), density and Young's modulus of T-shaped pier (denoted by D5 and E5), and cable tension (denoted by CN and CB). For the preselected parameters, set 10% change based on the initial values, and substitute them into the finite element model to calculate the frequency change percentage. The sensitivity of each parameter to the extracted modal frequency is shown in Figure  18.

Updating Parameters Selection and Establishment of Kriging Model
When selecting the updating parameters, the following aspects are mainly considered: (1) The steel box girder is a prefabricated component, so it is difficult to change the section area or the moment of inertia due to the size deviation of the section; (2) The tension of cable may relax under frequent loads; (3) There are different degrees of corrosion in straight-line segment and U-shape segment of steel box girder, which greatly affects the Young's modulus and density of beam; (4) The actual stiffness of the main tower should be greater than the design value when the section thickness of the main tower at the connection of the cable and the main tower increases; (5) T-shaped pier as a member supporting the main beam; its Young's modulus and density may change with the increase of service life.
Therefore, the following parameters are preliminarily selected as the updating parameters: density and Young's modulus of main tower (denoted by D 1 and E 1 ), density and Young's modulus of steel box girder in straight-line segment (denoted by D 2 and E 2 ), density and Young's modulus of steel box girder in U-shape segment (denoted by D 3 and E 3 ), density and Young's modulus of the cables (denoted by D 4 and E 4 ), density and Young's modulus of T-shaped pier (denoted by D 5 and E 5 ), and cable tension (denoted by CN and CB). For the preselected parameters, set 10% change based on the initial values, and substitute them into the finite element model to calculate the frequency change percentage. The sensitivity of each parameter to the extracted modal frequency is shown in Figure 18.

Model Updating and Result Analysis
The error objective function of fitting value of Kriging model and measured response is shown in Equation (18).
where x means the updating parameters, and α1 and α2 are the weight coefficients; in this case, α1 = 1, α2 = 1.
( ) ( ) i MAC x means the confidence degree that the theoretical mode and the measured mode are of the same order. It should be noted that the measured value of f5 is not a constraint condition of the objective function but is used to verify the results of the finite element model updating.
Running IAFSA, one group of global optimal solutions and fourteen groups of local optimal solutions of the objective function are obtained. Limited to space, a representative group of global optimal solutions and three groups of local optimal solutions are selected, as shown in Table 10. Table  11 show the difference of parameters before and after updating.

Model Updating and Result Analysis
The error objective function of fitting value of Kriging model and measured response is shown in Equation (18).
where x means the updating parameters, and α 1 and α 2 are the weight coefficients; in this case, α 1 = 1, α 2 = 1. f i Kri (x) represents the frequency calculated by Kriging model, f i t (x) represents the measured values of frequencies. MAC i (x) means the confidence degree that the theoretical mode and the measured mode are of the same order. It should be noted that the measured value of f 5 is not a constraint condition of the objective function but is used to verify the results of the finite element model updating.
Running IAFSA, one group of global optimal solutions and fourteen groups of local optimal solutions of the objective function are obtained. Limited to space, a representative group of global optimal solutions and three groups of local optimal solutions are selected, as shown in Table 10. Table 11 show the difference of parameters before and after updating. The corresponding parameters of each group of solutions are respectively substituted into the finite element model to extract the 1st, 2nd, 3rd, and 6th vertical bending modes and 1st torsional mode, and the relative error between the theoretical response corresponding to each group of solutions and the measured value is calculated as shown in Table 12. The MAC values under each group of solutions are shown in Table 13. From Table 12, it can be seen that the four groups of solutions can effectively reduce the error between the theoretical value and the measured value of the frequencies, and the global optimal solution has the most significant reduction effect. However, it is not comprehensive that the only criterion for parameters selection is to reduce the error between the theoretical values and the measured values of the model. Table 11 shows the change degree of parameters corresponding to each group of solutions before and after update. The parameters of the global optimal solution are as follows: the density of the straight-line segment of beam is reduced by 36.94%, the Young's modulus is reduced by 39.81%, while the density and Young's modulus of U-shaped segment are increased to 27.26% and 39.81%, respectively. Considering the actual situation of bridge diseases and combined with engineering experience, the corrosion of straight-line segment of the main beam is more serious than that of U-shaped segment. However, from the practical engineering point of view, the reduction of the Young's modulus and density of the main beam to nearly 40% is a small probability event. The reasons are as follows: (1) There is no obvious defect in the section of steel box girder, and the density and inertia moment of the girder will not be greatly reduced; (2) There are auxiliary facilities such as pavement and steel guardrail on the bridge deck system, which will buffer the decrease of the density and stiffness of the main beam caused by the corrosion of the steel structure to some extent.
Compared with the global optimal solution, the parameter value corresponding to the sub optimal solution 2 is more consistent with the actual bridge condition. The analysis is as follows: Compared with the initial value, the elastic modulus of the main tower increases by 28.16%, the density and elastic modulus of the straight-line segment of beam decrease by 12.87% and 7.77%, and the density and elastic modulus of U-shaped segment increase by 12.48% and 17.48%, respectively. Due to the increase of the section thickness of the connection between the cable and the main tower, and the existence of lifting lugs, stiffeners and other components, it is possible that the elastic modulus of the main tower is larger than the initial value. It should be noted that the measured value of f 5 (the 6-th vertical bending mode) is used to verify the results of the finite element model updating. Four groups of solutions are substituted into the FEM to calculate f 5 and MAC 5 . For the global optimal solution, the error between the theoretical value and the measured value of the 6-th vertical bending frequency is reduced from 3.41% to 2.88%, and the MAC 5 value is 0.950, for suboptimal solution 2, the error is reduced to 0.07%, and the corresponding MAC 5 value is 0.955. The graphs of the theoretical and experimental mode shapes of the vertical bending modes corresponding to the suboptimal solution 2 are shown in Figures 21-24. Therefore, in this case, the suboptimal solution 2 is selected as the solution that best represents the bridge parameters. (2) There are auxiliary facilities such as pavement and steel guardrail on the bridge deck system, which will buffer the decrease of the density and stiffness of the main beam caused by the corrosion of the steel structure to some extent.
Compared with the global optimal solution, the parameter value corresponding to the sub optimal solution 2 is more consistent with the actual bridge condition. The analysis is as follows: Compared with the initial value, the elastic modulus of the main tower increases by 28.16%, the density and elastic modulus of the straight-line segment of beam decrease by 12.87% and 7.77%, and the density and elastic modulus of U-shaped segment increase by 12.48% and 17.48%, respectively. Due to the increase of the section thickness of the connection between the cable and the main tower, and the existence of lifting lugs, stiffeners and other components, it is possible that the elastic modulus of the main tower is larger than the initial value. It should be noted that the measured value of f5 (the 6-th vertical bending mode) is used to verify the results of the finite element model updating. Four groups of solutions are substituted into the FEM to calculate f5 and MAC5. For the global optimal solution, the error between the theoretical value and the measured value of the 6-th vertical bending frequency is reduced from 3.41% to 2.88%, and the MAC5 value is 0.950, for suboptimal solution 2, the error is reduced to 0.07%, and the corresponding MAC5 value is 0.955. The graphs of the theoretical and experimental mode shapes of the vertical bending modes corresponding to the suboptimal solution 2 are shown in Figures 21-24. Therefore, in this case, the suboptimal solution 2 is selected as the solution that best represents the bridge parameters.   (2) There are auxiliary facilities such as pavement and steel guardrail on the bridge deck system, which will buffer the decrease of the density and stiffness of the main beam caused by the corrosion of the steel structure to some extent.
Compared with the global optimal solution, the parameter value corresponding to the sub optimal solution 2 is more consistent with the actual bridge condition. The analysis is as follows: Compared with the initial value, the elastic modulus of the main tower increases by 28.16%, the density and elastic modulus of the straight-line segment of beam decrease by 12.87% and 7.77%, and the density and elastic modulus of U-shaped segment increase by 12.48% and 17.48%, respectively. Due to the increase of the section thickness of the connection between the cable and the main tower, and the existence of lifting lugs, stiffeners and other components, it is possible that the elastic modulus of the main tower is larger than the initial value. It should be noted that the measured value of f5 (the 6-th vertical bending mode) is used to verify the results of the finite element model updating. Four groups of solutions are substituted into the FEM to calculate f5 and MAC5. For the global optimal solution, the error between the theoretical value and the measured value of the 6-th vertical bending frequency is reduced from 3.41% to 2.88%, and the MAC5 value is 0.950, for suboptimal solution 2, the error is reduced to 0.07%, and the corresponding MAC5 value is 0.955. The graphs of the theoretical and experimental mode shapes of the vertical bending modes corresponding to the suboptimal solution 2 are shown in Figures 21-24. Therefore, in this case, the suboptimal solution 2 is selected as the solution that best represents the bridge parameters.

Conclusions
(1) An improved artificial fish swarm algorithm is proposed, which introduces niche technology, adaptive niche radius mechanism, and uses simulated annealing algorithm to locally refine and optimize in the later stage of the algorithm. Through three typical functions, IAFSA is verified to have high efficiency, high precision, and intelligent search for multiple solutions of multiple peak functions. (2) In the case of ASCE benchmark finite element model updating, one set of global optimal solutions and three sets of local optimal solutions of the objective function are obtained. By comparing the preset damage value with the four sets of solutions, it is found that the parameters corresponding to the global optimal solution are different from the preset values, while the suboptimal solution 1 can match well with the preset values, which proves that the global optimal solution does not necessarily represent the actual parameters value of the structure.

Conclusions
(1) An improved artificial fish swarm algorithm is proposed, which introduces niche technology, adaptive niche radius mechanism, and uses simulated annealing algorithm to locally refine and optimize in the later stage of the algorithm. Through three typical functions, IAFSA is verified to have high efficiency, high precision, and intelligent search for multiple solutions of multiple peak functions. (2) In the case of ASCE benchmark finite element model updating, one set of global optimal solutions and three sets of local optimal solutions of the objective function are obtained. By comparing the preset damage value with the four sets of solutions, it is found that the parameters corresponding to the global optimal solution are different from the preset values, while the suboptimal solution 1 can match well with the preset values, which proves that the global optimal solution does not necessarily represent the actual parameters value of the structure.

Conclusions
(1) An improved artificial fish swarm algorithm is proposed, which introduces niche technology, adaptive niche radius mechanism, and uses simulated annealing algorithm to locally refine and optimize in the later stage of the algorithm. Through three typical functions, IAFSA is verified to have high efficiency, high precision, and intelligent search for multiple solutions of multiple peak functions. (2) In the case of ASCE benchmark finite element model updating, one set of global optimal solutions and three sets of local optimal solutions of the objective function are obtained. By comparing the preset damage value with the four sets of solutions, it is found that the parameters corresponding to the global optimal solution are different from the preset values, while the suboptimal solution 1 can match well with the preset values, which proves that the global optimal solution does not necessarily represent the actual parameters value of the structure.