Bearing Damage Detection of a Reinforced Concrete Plate Based on Sensitivity Analysis and Chaotic Moth-Flame-Invasive Weed Optimization

This article proposes a novel damage detection method based on the sensitivity analysis and chaotic moth-flame-invasive weed optimization (CMF-IWO), which is utilized to simultaneously identify the damage of structural elements and bearings. First, the sensitivity coefficients of eigenvalues to the damage factors of structural elements and bearings are deduced, the regularization technology is used to solve the problem of equation undetermined, meanwhile, the modal strain energy-based index is utilized to detect the damage locations, and the regularization objective function is constructed to quantify the damage severity. Then, for the subsequent procedure of damage detection, CMF-IWO is proposed based on moth-flame optimization and invasive weed optimization as well as chaos theory, reverse learning, and evolutional strategy. The optimization effectiveness of the hybrid algorithm is verified by five benchmark functions and a damage identification numerical example of a simply supported beam; the results demonstrate it is of great global search ability and higher convergence efficiency. After that, a numerical example of an 8-span continuous beam and an experimental reinforced concrete plate are both adopted to evaluate the proposed damage identification method. The results of the numerical example indicate that the proposed method can locate and quantify the damage of structural elements and bearings with high accuracy. Furthermore, the outcomes of the experimental example show that despite the existence of some errors and uncertain factors, the method still obtains an acceptable result. Generally speaking, the proposed method is proved that it is of good feasibility.


Introduction
Civil structures suffer from traffic load, environmental temperature variation, fatigue failure, and other uncertain negative influences during the service period. Among the distinguished damage identification methods, the modal parameters-based detection method is known as the most popular one. There are some major reasons for understanding its popularity. First, the modal parameters, such as natural frequencies and modal shapes, are obtained easily by modal tests. Meanwhile, prior studies have addressed many issues, the common one of them is the sensitivity coefficient analysis.
The existing studies have provided several different sensitivity coefficient analyses, including the natural frequencies [1], modal shapes [2], modal strain energy [3], and modal kinetic energy [4,5]. For the natural frequencies, Xia and Hao [6] proposed a method of statistical damage identification based on the frequencies sensitivity analysis and finite element model updating, which was successfully applied in a numerical cantilever beam and a laboratory tested steel cantilever plate, meanwhile, where K and M denote the global stiffness and mass matrices of a structure, respectively; λ i and ϕ i indicate the i-th eigenvalue and eigenvector, respectively. Very few researches try to establish a model with structural boundary condition, and it is only treated as a simple constraint. However, this issue can be solved by dividing the structure into two systems, one is the structural system and the other is the boundary condition system [1]. Thus, Equation (1) can be rewritten as follows: where K s and K bc mean the stiffness matrices of the structural and boundary condition system, respectively; M s and M bc stand for the mass matrices of the structural system and boundary condition system, respectively. Also, the stiffness matrix of the boundary condition system can be defined as follows [1]: Ux Rx Uy Ry Uz Rz Ux Rx Uy Ry Uz Rz where the elements on the diagonal of the matrix denote the boundary stiffness in every degree of freedom (DOF); others mean the synergistic effect, which can be explained that the translation in one direction will cause rotations in the other two directions. Meanwhile, this matrix is not only available for a 3-D finite element model, but also can be applied for the 2-D model by reducing the unrelated DOFs. When the damage of structural elements and/or bearings arise in the structure, damage can be seen as the stiffness reduction with no mass change, thus, the damage of structural elements and bearings can be quantitatively measured by using stiffness reduction vectors: where θ k and α j represent the damage severity of the k-th element and the j-th bearing. Thereby, the damage in structural element and bearing can be described as follows: (1 − θ k )k s k , 0 ≤ θ k ≤ 1 (6) (1 − α j )k bc j , 0 ≤ α j ≤ 1 Sensors 2020, 20, 5488 5 of 28 where nele is the number of the structural elements; nbc is the number of bearings.

Sensitivity Analysis of Eigenvalue
The first-order sensitivity coefficients of eigenvalue corresponding to the structural damage factor θ k and the bearing damage factor α j can be obtained by derivative of Equation (2) concerning the two parameters, respectively, which can be written as follows: Then, the above equations are premultiplied by ϕ T i in both sides. At the same time, using the known relations, such as ϕ T i (M s + M bc )ϕ i = I for the unit-mass normalized mode shapes; the mass matrices M s and M bc are independent of the damage factors of structural elements and bearings; stiffness matrices K s and K bc are respectively independent of α j and θ k , namely, ∂M s /∂θ k = 0, ∂M bc /∂α j = 0, ∂K bc /∂θ k = 0 and ∂K s /∂α j = 0. Thus, the sensitivity coefficients of i-th eigenvalue corresponding to θ k and α j can be derived respectively as follows: According to Equations (6) and (7), the following equations can be derived as: where K s k and K bc j denote the k-th structural element stiffness matrix and the j-th boundary condition stiffness matrix in the global coordinate, respectively. Thus, the first-order sensitivity coefficients of eigenvalue corresponding to θ k and α j can be rewritten as follows:

Sensitivity Analysis of Eigenvector
According to the paper of Zhao and DeWolf [34], the first-order sensitivity coefficient of the eigenvector corresponding to the structural damage factor θ k can be represented as follows: Sensors 2020, 20, 5488 6 of 28 where ndof stands for the total number of structural DOFs; a in denotes the n-th undetermined coefficient in the sensitivity coefficient of the i-th eigenvector to the structural damage factor θ k . As for a in , there are two possible situations: (1) When the subscript i n, Equation (8) is premultiplied by ϕ T n on its both sides, then substituting Equation (16) into Equation (8), the following equation can be derived: Because of the orthogonality of mode shapes, the mathematical relation can be concluded, namely, if i n, ϕ T n (M s + M bc )ϕ i = 0; ϕ T n (M s + M bc )ϕ n = I; ϕ T n (K s + K bc ) = λ n ϕ T n (M s + M bc ); and ∂M s /∂θ k = 0. Thus, Equation (17) can be rewritten as follows: For a in , it can be solved as: (2) When the subscript i = n, as for the unit-mass normalized mode shapes, it can obtain: Aiming Equation (20), take derivative in terms of θ k , which can be written as: Noting the symmetry characteristic of the mass matrix: Considering Equation (16), the above equation can be rewritten as: Because of the orthogonality of mode shapes and Equation (23), it can be indicated that: Based on the above analysis, the sensitivity coefficient of a mode shape corresponding to the structural element damage factor can be summarized as follows: Sensors 2020, 20, 5488 7 of 28 Meanwhile, the same derivation procedure may be easily adapted to obtain the sensitivity coefficient of mode shapes corresponding to bearing damage factor α j , which can be defined as: Because of the eigenvalues of the structure are not only influenced by the damage of structural elements but also of the bearings, therefore, the sensitivity coefficient of eigenvalue to the two factors should be considered. Based on Equations (14) and (15), the equations that can be obtained are as follows: Expanding the above equations as: Noting ∂K bc /∂θ k = 0, ∂K s /∂α j = 0, S s λ = −ϕ T i K s k ϕ i and S bc λ = −ϕ T i K bc j ϕ i , Equations (29) and (30) are premultiplied by ϕ T i and ϕ i on their right and left sides respectively: Then the equations can be simplified by shifting terms: ∂S bc Hence, Equations (33) and (34) can be rewritten as follows: Sensors 2020, 20, 5488 , add Equations (35) and (36), then the sensitivity coefficient can be obtained as follows:

Objective Function Based on Regularization Technology
Regularization technology, as an excellent tool, has been widely applied in solving the problem of underdetermined equation in damage identification [10]. Assuming the change of modal parameters caused by damage is linear, it can be expressed as: where S denotes the sensitivity coefficient matrix, which can be calculated according to Equation (37); x = [θ, α] T is the combination of the stiffness reduction vectors of elemental damage and bearing damage; ∆λ means the change of eigenvalue. Because of only a small number of damaged elements and bearings, x will be a sparse vector. Meanwhile, to consider the limitation of sensors and incomplete modal parameters, the dimension of ∆λ is far less than x, which indicates the equation is underdetermined. Therefore, in order to obtain the expected sparsest solution, Equation (38) can be transformed into the optimization problem as follows: Noting ∆λ = λ e − λ a , the superscript e and a represent the experimental and analytical eigenvalues, respectively.
The above equation can be rewritten as follow: where λ(x) = Sx + λ a ; ε is the error tolerance. Moreover, based on Equation (40), the following l 1 regularization can be obtained: where µ > 0 represents the regularization parameter, which is determined using the L-curve approach [11] or calculated by the equation of µ = σ 2 log(p), in which p is the cardinality of S [35]. Considering the damage identification problem and obtained eigenvalues, Equation (41) is further written as: where λ a i and λ e i are the i-th analytical and experimental eigenvalues, respectively; m means the considered number of modes; n represents the length of x.

Damage Location Based on MSEBI
In a structure, the damage of structural elements and bearings often exist simultaneously, which increases the difficulty of damage detection and the cost of optimization computation. Aiming to this problem, MSEBI is utilized to locate the locations of damaged elements and bearings. The modal strain energy of the element and the bearing can be calculated as follows: (ϕ e i ) T k e k ϕ e i , k = 1, 2, · · · , nele, i = 1, 2, · · · , nm (43) where ϕ e i and ϕ bc i are the nodal displacement vectors of the k-th element and j-th bearing corresponding to the i-th mode shape, respectively; k e k and k bc j stand for the k-th elemental stiffness matrix and the j-th bearing stiffness matrix, respectively; nm represents the mode order considered; nele and nb are the total numbers of structural elements and bearings, respectively.
According to the above equations, the i-th total modal strain energy of the structure and bearings can be calculated as follows: For the convenience of calculation, the normalized i-th modal strain energy of element and bearing can be defined as follows: Then taking the mean value of first nm-order normalized modal strain energy of element and bearing, Equations (47) and (48) can be rewritten as follows: , e = 1, 2, · · · , nele (49) Thus, MSEBI can be obtained as: where max[] represents the action of taking the maximum value; the superscripts of E and A mean experimental and analytical respectively; when the analytical and experimental modal strain energy are the same, MSEBI = 0, namely, the element or bearing is intact, otherwise, the damage may occur (MSEBI > 0).

Moth-Flame Optimization
Moth-flame optimization (MFO), as a novel optimization tool, was proposed by Mirjalili [36]. The inspiration of MFO can be traced back to the transverse orientation of moths at night, but it was developed by the approach of a moth flying around the flame or candle. The basic theory of MFO can be explained that some individuals of the moth with an attribute of position are first initialized randomly in a D-dimensional solution space: where n denotes the number of moths; D represents the dimension of optimization problem. At the same time, the artificial light, namely, flame, will be marked as follows: The fitness values of each moth individual and flame are stored in two vectors, which can be shown as follows: where OM and OF represent the fitness value vectors of moth individual and flame, respectively; it means the current number of the iteration. Then, the logarithmic spiral function is utilized to update the position of each moth: where D i stands for the spatial distance from i-th moth to j-th flame, the constant b is a factor to define the spiral shape function, and t is a random number between −1 and 1. Furthermore, the adaptive decrease mechanism of flames is incorporated into optimization algorithm to ensure the powerful exploitation of optimization process, which also keeps the moth individual always flying around the optimal solution from the first iteration to iteration termination. The mathematical formula of the mechanism can be written as follows: where Fn it and Fn max are the flame number of the it-th iteration and the max number of flame; Iteration stands for the current count of iteration; Iteration max means the maximum number of iteration, and round() denotes the action of taking the integer portion.

Invasive Weed Optimization
Invasive weed optimization (IWO) is inspired by the situation of colonization of invasive weeds [37,38]. The colonizing behavior of weeds can be described in the formalization language of mathematics, namely, some weeds are randomly generated in the D-dimensional problem space with the characteristic of the position: where n denotes the initial number of weed; D represents the dimension of the problem. Then, based on the fitness of each weed, the seed of each one is calculated by the equation as follows: where f (w i ) means the fitness value of i-th weed; f min and f max are the worst and best fitness in the current population, respectively; s min and s max represent the minimum and maximum seed number that one weed can produce, respectively. Subsequently, the produced seeds are spread over the search space using a normally distributed random number, whose standard deviation at current iteration can be expressed by: where Iteration stands for the current count of iteration; Iteration max means the maximum number of iteration; σ initial and σ f inal represent the initial and final values of standard deviation, respectively; n = 3 means the nonlinear modulation coefficient. Thus, the position of weeds can be updated by the equation as follows: Furthermore, the principle of the growing competition is adopted in the algorithm, and it can be described that while evaluating the fitness of all weeds and seeds, the poor weeds and/or seeds are eliminated to reach the maximum number of preset populations.

Chaotic Moth-Flame-Invasive Weed Optimization Hybrid Algorithm
MFO is of powerful local searching ability, but it is unable to ensure the performance in global search, especially for the diversity of moth individuals in the late iterations is poor. Therefore, a hybrid algorithm, chaotic moth-flame-invasive weed optimization hybrid algorithm (CMF-IWO), is proposed to obtain a better optimization result. In the hybrid algorithm, the mechanism of seed spreading and growing competition are incorporated into MFO, meanwhile, the population initialization approach of reverse learning [39], chaos theory [40], and evolutional strategy are also adopted to enhance the diversity of the population.
The flowchart of CMF-IWO is shown in Figure 1. The basic procedures of CMF-IWO can be summarized as follows: (1) To initialize the chaotic populations of n based on the chaos theory, the related equation can be defined as follows: where lb and ub are the low and upper bounds respectively; ξ i denotes the chaotic vector, which can be generated by the logistic chaos mapping: where u is a scalar, when u = 4, the system is in chaos. (2) To apply the operation of reverse learning, the individuals of reverse learning can be produced as follows: Then merge w i and w RL i , and the fitness of each individual is evaluated. After that, according to the ranking of the fitness, the first n weeds are selected as moth individuals to be input into MFO.
Sensors 2020, 20, 5488 12 of 28 (3) In MFO, the moth individuals are spread over the search space according to Equation (61) to enhance the diversity, and the evolutional strategy is also incorporated in this stage. The operation of mutation can be described as follows: where F = (F max − F min ) * rand is the scale factor, m ij , m ik , and m ir are the j-th, k-th, and r-th element of i-th moth, respectively. The operations of crossover and selection can be defined as: where the pCR represents the probability of crossover. (4) Based on the obtained individuals of the previous step, the first n moth individuals are selected according to their fitness values. Then the remaining steps of MFO are conducted to obtain the optimization results.
Sensors 2020, 20, x FOR PEER REVIEW 12 of 27 is the scale factor, ij m , ik m , and ir m are the j-th, k-th, and r-th element of i-th moth, respectively. The operations of crossover and selection can be defined as: where the pCR represents the probability of crossover.
(4) Based on the obtained individuals of the previous step, the first n moth individuals are selected according to their fitness values. Then the remaining steps of MFO are conducted to obtain the optimization results.  (61) and (62) Merging

Evaluation Using Benchmark Functions
The optimization ability and computational accuracy of the hybrid algorithm are first evaluated with five mathematical benchmark functions (Table 1) and compared with existing optimization algorithms, such as MFO, IWO, Particle Swarm Optimization (PSO), Cuckoo Search (CS), and Differential Evolution (DE). The parameter settings of each algorithm are listed in Table 2, and each algorithm is performed 50 times with maximum iterations of 500 and maximum populations of 100, the average results and iterative curves are shown in Table 3 and Figure 2.

Function
Definition    From Table 3, it can be seen that the proposed hybrid algorithm can achieve better optimal results in the optimizations of five benchmark functions, which can be owed to the chaotic population, and the mechanism of seed space spreading can enhance the diversity of the initial population. Meanwhile, the operation of reverse learning can obtain elite populations. These two improvements can guarantee that the initial populations are of high quality and diversity. In addition, the iterative curves of Figure 2 indicate that the curves of CMF-IWO are steeper than those of other algorithms, which demonstrates that the convergence speed of CMF-IWO is good. At the same time, compared to the other algorithms, the proposed algorithm can escape from local optimal, whose main reason may be the fusion of evolutional strategy.

Evaluation Using Numerical Example of Structural Damage Identification
In order to further assess the optimization ability of the proposed hybrid algorithm, a simply supported beam with 16 elements is exploited (Figure 3). For the beam, its Young's modulus is 3.0 × 10 10 Pa, the mass density is 2450 kg/m 3 , the cross-sectional area is 0.05 m 2 , the inertia moment is 4.16 × 10 -5 m 4 , the length of each element is 0.5 m. Three damage cases, including single-point damage, double-point damage, and multiple-point damage, are introduced by the reduction of the stiffness, the details are listed in Table 4.

Case
Damage Severity @ Element Number 1 10% @ 3 2 10% @ 3, 5% @ 7 3 10% @ 3, 5% @ 7, 10% @ 13 Because of the beam with the simply supported boundary conditions, namely, it has no bearings on each side. The bearing damage factor j  is defined as zeros, only the elemental damage factor k  is considered. The parameters setting of six algorithms are the same as Table 2, each algorithm is performed seven times with maximum iteration of 500 and maximum population of 100, the average damage identification results are illustrated in Figure 4.  From Table 3, it can be seen that the proposed hybrid algorithm can achieve better optimal results in the optimizations of five benchmark functions, which can be owed to the chaotic population, and the mechanism of seed space spreading can enhance the diversity of the initial population. Meanwhile, the operation of reverse learning can obtain elite populations. These two improvements can guarantee that the initial populations are of high quality and diversity. In addition, the iterative curves of Figure 2 indicate that the curves of CMF-IWO are steeper than those of other algorithms, which demonstrates that the convergence speed of CMF-IWO is good. At the same time, compared to the other algorithms, the proposed algorithm can escape from local optimal, whose main reason may be the fusion of evolutional strategy.

Evaluation Using Numerical Example of Structural Damage Identification
In order to further assess the optimization ability of the proposed hybrid algorithm, a simply supported beam with 16 elements is exploited (Figure 3). For the beam, its Young's modulus is 3.0 × 10 10 Pa, the mass density is 2450 kg/m 3 , the cross-sectional area is 0.05 m 2 , the inertia moment is 4.16 × 10 −5 m 4 , the length of each element is 0.5 m. Three damage cases, including single-point damage, double-point damage, and multiple-point damage, are introduced by the reduction of the stiffness, the details are listed in Table 4. From Table 3, it can be seen that the proposed hybrid algorithm can achieve better optimal results in the optimizations of five benchmark functions, which can be owed to the chaotic population, and the mechanism of seed space spreading can enhance the diversity of the initial population. Meanwhile, the operation of reverse learning can obtain elite populations. These two improvements can guarantee that the initial populations are of high quality and diversity. In addition, the iterative curves of Figure 2 indicate that the curves of CMF-IWO are steeper than those of other algorithms, which demonstrates that the convergence speed of CMF-IWO is good. At the same time, compared to the other algorithms, the proposed algorithm can escape from local optimal, whose main reason may be the fusion of evolutional strategy.

Evaluation Using Numerical Example of Structural Damage Identification
In order to further assess the optimization ability of the proposed hybrid algorithm, a simply supported beam with 16 elements is exploited (Figure 3). For the beam, its Young's modulus is 3.0 × 10 10 Pa, the mass density is 2450 kg/m 3 , the cross-sectional area is 0.05 m 2 , the inertia moment is 4.16 × 10 -5 m 4 , the length of each element is 0.5 m. Three damage cases, including single-point damage, double-point damage, and multiple-point damage, are introduced by the reduction of the stiffness, the details are listed in Table 4.

Case
Damage Severity @ Element Number 1 10% @ 3 2 10% @ 3, 5% @ 7 3 10% @ 3, 5% @ 7, 10% @ 13 Because of the beam with the simply supported boundary conditions, namely, it has no bearings on each side. The bearing damage factor j α is defined as zeros, only the elemental damage factor k θ is considered. The parameters setting of six algorithms are the same as Table 2, each algorithm is performed seven times with maximum iteration of 500 and maximum population of 100, the average damage identification results are illustrated in Figure 4.   10% @ 3, 5% @ 7 3 10% @ 3, 5% @ 7, 10% @ 13 Because of the beam with the simply supported boundary conditions, namely, it has no bearings on each side. The bearing damage factor α j is defined as zeros, only the elemental damage factor θ k is considered. The parameters setting of six algorithms are the same as Table 2, each algorithm is performed seven times with maximum iteration of 500 and maximum population of 100, the average damage identification results are illustrated in Figure 4.  As shown in Figure 4, for the single-point damage, CMF-IWO can determine the damage location and quantify the damage severity with a higher accuracy. However, other algorithms cannot achieve satisfactory results. Then in Case 2, some errors occur in all the six algorithms, the main reason is that the natural frequencies do not contain the spatial information of structural damage, so it is not easy to obtain very precise results by using merely frequencies. However, comparatively, the error of the proposed algorithm is less than the other five algorithms, which indicates the optimization performance of CMF-IWO is better than other algorithms. Furthermore, the same conclusion can be obtained in Case 3.
To summarize, with limited modal characteristics, by using the proposed CMF-IWO hybrid algorithm, better damage identification results can be achieved than the other five algorithms, which has better potential in structural damage detection.

Damage Detection Methodology
Based on the proposed regularization objective function and CMF-IWO hybrid algorithm, a novel damage detection approach is put forward, which can identify the damage of bearings, as well as the damage of structural elements. The main steps of this method can be summarized as follows: (1) MSEBI of structural elements and bearings is calculated and used to detect the locations of damage;  As shown in Figure 4, for the single-point damage, CMF-IWO can determine the damage location and quantify the damage severity with a higher accuracy. However, other algorithms cannot achieve satisfactory results. Then in Case 2, some errors occur in all the six algorithms, the main reason is that the natural frequencies do not contain the spatial information of structural damage, so it is not easy to obtain very precise results by using merely frequencies. However, comparatively, the error of the proposed algorithm is less than the other five algorithms, which indicates the optimization performance of CMF-IWO is better than other algorithms. Furthermore, the same conclusion can be obtained in Case 3.
To summarize, with limited modal characteristics, by using the proposed CMF-IWO hybrid algorithm, better damage identification results can be achieved than the other five algorithms, which has better potential in structural damage detection.

Damage Detection Methodology
Based on the proposed regularization objective function and CMF-IWO hybrid algorithm, a novel damage detection approach is put forward, which can identify the damage of bearings, as well as the damage of structural elements. The main steps of this method can be summarized as follows: (1) MSEBI of structural elements and bearings is calculated and used to detect the locations of damage; (2) According to the detected damage location, the sensitivity coefficients of the suspected structural elements and bearings are calculated, which are adopted to construct the regularization objective function; (3) The regularization objective function constructed in the previous step is input into CMF-IWO hybrid algorithm to quantify the damage severities of structural elements and bearings.

Numerical Study
In this section, as shown in Figure 5, an 8-span continuous beam with 48 elements and 9 bearings is introduced to verify the proposed damage identification method. The length and cross-sectional area of each element are 0.5 m and 0.03 m 2 , respectively, its material properties, such as Young's modulus, mass density, and inertia moment are 3.45 × 10 10 Pa, 2500 kg/m 3 , and 2.5 × 10 −5 m 4 , respectively. The vertical stiffness of the bearings is 1.0 × 10 6 KN/m. The damage of structural elements and bearings are both simulated using reduction of stiffness [41][42][43][44]. Five damage cases are introduced, which are shown in Table 5. (2) According to the detected damage location, the sensitivity coefficients of the suspected structural elements and bearings are calculated, which are adopted to construct the regularization objective function; (3) The regularization objective function constructed in the previous step is input into CMF-IWO hybrid algorithm to quantify the damage severities of structural elements and bearings.

Numerical Study
In this section, as shown in Figure 5, an 8-span continuous beam with 48 elements and 9 bearings is introduced to verify the proposed damage identification method. The length and cross-sectional area of each element are 0.5 m and 0.03 m 2 , respectively, its material properties, such as Young's modulus, mass density, and inertia moment are 3.45 × 10 10 Pa, 2500 kg/m 3 , and 2.5 × 10 −5 m 4 , respectively. The vertical stiffness of the bearings is 1.0 × 10 6 KN/m. The damage of structural elements and bearings are both simulated using reduction of stiffness [41][42][43][44]. Five damage cases are introduced, which are shown in Table 5.

Numerical Study
In this section, as shown in Figure 5, an 8-span continuous beam with 48 elements and 9 bearings is introduced to verify the proposed damage identification method. The length and cross-sectional area of each element are 0.5 m and 0.03 m 2 , respectively, its material properties, such as Young's modulus, mass density, and inertia moment are 3.45 × 10 10 Pa, 2500 kg/m 3 , and 2.5 × 10 −5 m 4 , respectively. The vertical stiffness of the bearings is 1.0 × 10 6 KN/m. The damage of structural elements and bearings are both simulated using reduction of stiffness [41][42][43][44]. Five damage cases are introduced, which are shown in Table 5.

Case Damage Severity @ Element Number
Damage Severity @ Bearing Number 1 30% @ 24 25% @ 2# 2 20% @ 5, 30% @ 24 25% @ 2#, 50% @ 7# 3 20% @ 5, 30% @ 24, 40% @ 40 25% @ 2#, 75% @ 5#, 50% @ 7# 4 20% @ 1, 40% @ 48 25% @ 1#, 50% @ 9# 5 20% @ 1, 40% @ 48 99.9% @ 1#, 99.9% @ 9# The first six modal parameters are exploited to detect the damage of structural elements and bearings. Also, the maximum iterations, populations, and parameters setting are the same as Section 3.4. The results of the damage location are depicted in Figures 6 and 7.  As shown in Figure 6, the suspected damaged elements of the five damage cases are detected accurately. Also, in Cases 1 to 3, there still exist some elements whose identified damage are not zero, but compared to the damaged elements, the values are so small that they can be ignored. Hence, it is revealed that the suspected damaged elements for Cases 1 to 3 can be reduced from 48 to 1, 2, and 3 elements respectively. On the other hand, from Figure 7, it is clearly observed that MSEBI can locate the locations of damaged bearings with high accuracy. To summarize, MSEBI not only can detect the damage location of structural elements but also identify the locations of damaged bearings.  As shown in Figure 6, the suspected damaged elements of the five damage cases are detected accurately. Also, in Cases 1 to 3, there still exist some elements whose identified damage are not zero, but compared to the damaged elements, the values are so small that they can be ignored. Hence, it is revealed that the suspected damaged elements for Cases 1 to 3 can be reduced from 48 to 1, 2, and 3 elements respectively. On the other hand, from Figure 7, it is clearly observed that MSEBI can locate the locations of damaged bearings with high accuracy. To summarize, MSEBI not only can detect the damage location of structural elements but also identify the locations of damaged bearings. Because modal parameters are more sensitive to the damage of structural elements than the damage of bearings, if the damage severities of elements and bearings are simultaneously determined in the process of damage detection, the change of modal parameters caused by the damage of bearings will be masked by elemental damage. Thus, at first, assuming the bearings are intact, the damage severity of structural elements is determined. After that, the determined elemental damage condition is input to the procedure of damage detection, and the damage severity of bearing is quantified. The average results of seven times iterations are demonstrated in Figures 8 and 9. (a) (b) 6   As shown in Figure 6, the suspected damaged elements of the five damage cases are detected accurately. Also, in Cases 1 to 3, there still exist some elements whose identified damage are not zero, but compared to the damaged elements, the values are so small that they can be ignored. Hence, it is revealed that the suspected damaged elements for Cases 1 to 3 can be reduced from 48 to 1, 2, and 3 elements respectively. On the other hand, from Figure 7, it is clearly observed that MSEBI can locate the locations of damaged bearings with high accuracy. To summarize, MSEBI not only can detect the damage location of structural elements but also identify the locations of damaged bearings.
Because modal parameters are more sensitive to the damage of structural elements than the damage of bearings, if the damage severities of elements and bearings are simultaneously determined in the process of damage detection, the change of modal parameters caused by the damage of bearings will be masked by elemental damage. Thus, at first, assuming the bearings are intact, the damage severity of structural elements is determined. After that, the determined elemental damage condition is input to the procedure of damage detection, and the damage severity of bearing is quantified. The average results of seven times iterations are demonstrated in Figures 8 and 9.
determined in the process of damage detection, the change of modal parameters caused by the damage of bearings will be masked by elemental damage. Thus, at first, assuming the bearings are intact, the damage severity of structural elements is determined. After that, the determined elemental damage condition is input to the procedure of damage detection, and the damage severity of bearing is quantified. The average results of seven times iterations are demonstrated in Figures 8 and 9.   As shown in Figure 8, because of the location operation of MSEBI, the proposed method can detect the damaged elements with high accuracy. However, regarding the quantification of damage severity, the performance is not very satisfactory. Only for Case 1, the result is precise, different degrees of error occur in all the other cases, this phenomenon can be explained as follows: (1) The eigenvalues do not include the spatial information of structural damage, which cannot obtain the accurate detection results; (2) the existence of bearings damage may apply some adverse influence on the identification performance. 5 10 15 20 25   detect the damaged elements with high accuracy. However, regarding the quantification of damage severity, the performance is not very satisfactory. Only for Case 1, the result is precise, different degrees of error occur in all the other cases, this phenomenon can be explained as follows: (1) The eigenvalues do not include the spatial information of structural damage, which cannot obtain the accurate detection results; (2) the existence of bearings damage may apply some adverse influence on the identification performance.   The detection results of bearings damage are shown in Figure 9, it can be observed that the damage of Case 1, Case 4, and Case 5 can be quantified with satisfactory precision. However, there are different errors in other damage cases, such as damage severity shifting for Case 2 and inaccurate quantification results for Case 3. The reason also can be owed to the missing spatial information of eigenvalues, at the same time, the previous damage identification results of structural elements are the key factors to the damage detection of bearings. On the other hand, comparing Case 2, Case 4, and Case 5, it can be found that the more serious the damage severity of bearings, the more severe the fluctuations of the eigenvalues. Meanwhile, the comparison also indicates that the proposed method is more suitable to detect the serious damage of bearings, like bearing separation.  As shown in Figure 8, because of the location operation of MSEBI, the proposed method can detect the damaged elements with high accuracy. However, regarding the quantification of damage severity, the performance is not very satisfactory. Only for Case 1, the result is precise, different degrees of error occur in all the other cases, this phenomenon can be explained as follows: (1) The eigenvalues do not include the spatial information of structural damage, which cannot obtain the accurate detection results; (2) the existence of bearings damage may apply some adverse influence on the identification performance.

Experimental Example
The detection results of bearings damage are shown in Figure 9, it can be observed that the damage of Case 1, Case 4, and Case 5 can be quantified with satisfactory precision. However, there are different errors in other damage cases, such as damage severity shifting for Case 2 and inaccurate quantification results for Case 3. The reason also can be owed to the missing spatial information of eigenvalues, at the same time, the previous damage identification results of structural elements are the key factors to the damage detection of bearings. On the other hand, comparing Case 2, Case 4, and Case 5, it can be found that the more serious the damage severity of bearings, the more severe the fluctuations of the eigenvalues. Meanwhile, the comparison also indicates that the proposed method is more suitable to detect the serious damage of bearings, like bearing separation.

Experimental Example
In this section, an experimental example of a simply supported reinforced concrete plate is adopted to further assess the proposed damage detection approach. As shown in Figure 10a, the plate is located in the campus of Wuhan Institution of Technology; its measured size is 5.4 m × 0.6 m × 0.12 m (length × width × thickness), with the Young's modulus of 3 × 10 10 Pa, the mass density of 2410 kg/m 3 . Meanwhile, there were four rubber bearings placed on the two ends to support the plate. Additionally, the lengths of both the overhangs were 0.2 m at the two ends. In the natural environment, irregular hammer excitation was conducted to make the plate vibrate, ten accelerometers were installed on the top surface of the plate to collect the signal of acceleration. The layout of sensors can be seen as Figure 10b. The acceleration data were captured by a DH5922 vibration testing system ( Figure 11) which has the advantages of light mass and convenient use etc. The operation temperature of acceleration acquisition instrument ranges from 0 °C to 60 °C and it is also a universal dynamic signal test and analysis system which can complete the testing and analysis of stress, strain, vibration, shock, etc. The acceleration data were captured by a DH5922 vibration testing system ( Figure 11) which has the advantages of light mass and convenient use etc. The operation temperature of acceleration acquisition instrument ranges from 0 • C to 60 • C and it is also a universal dynamic signal test and analysis system which can complete the testing and analysis of stress, strain, vibration, shock, etc. The instrument has 16 24-bit IEPE input channels that are equipped with an anti-mixing filter, and supports sampling frequency up to 51.2 k Hz. The system was connected with the acceleration sensor by L5 coaxial extension wire and placed in the center of the equal dividing line to collect the acceleration signal. The operation temperature of acceleration sensor ranged from −40 • C to 80 • C and its tolerance was ±1%.
(b) Figure 10. The simply supported reinforced concrete plate: (a) view of modal test; (b) layout of sensors and mass block.
The acceleration data were captured by a DH5922 vibration testing system ( Figure 11) which has the advantages of light mass and convenient use etc. The operation temperature of acceleration acquisition instrument ranges from 0 °C to 60 °C and it is also a universal dynamic signal test and analysis system which can complete the testing and analysis of stress, strain, vibration, shock, etc. The instrument has 16 24-bit IEPE input channels that are equipped with an anti-mixing filter, and supports sampling frequency up to 51.2 k Hz. The system was connected with the acceleration sensor by L5 coaxial extension wire and placed in the center of the equal dividing line to collect the acceleration signal. The operation temperature of acceleration sensor ranged from −40 °C to 80 °C and its tolerance was ±1%. Meanwhile, in the experiment, because of its difficulties in introducing the stiffness change of concrete plate, the damage was simulated by the approach of applying additional mass block with a length of 0.3 m, width of 0.2 m, and thickness of 0.15 m ( Figure 12); moreover, the placed location was the center of the plate. Also, bearing 2# was removed to simulate the common disease, namely, bearing separation. Eight cases were set, and corresponding modal tests were conducted, the details of each case and the measured natural frequencies are listed in Tables 6 and 7, respectively. Sensors 2020, 20, x FOR PEER REVIEW 23 of 27 Meanwhile, in the experiment, because of its difficulties in introducing the stiffness change of concrete plate, the damage was simulated by the approach of applying additional mass block with a length of 0.3 m, width of 0.2 m, and thickness of 0.15 m ( Figure 12); moreover, the placed location was the center of the plate. Also, bearing 2# was removed to simulate the common disease, namely, bearing separation. Eight cases were set, and corresponding modal tests were conducted, the details of each case and the measured natural frequencies are listed in Tables 6 and 7, respectively.    Then, MATLAB is used to construct the finite element model of the plate. The plate is meshed to 66 elements, which is modeled by 20-node shell element; furthermore, the rubber bearings were simulated by the 3-D spring elements. The numbering and meshing diagram of elements and nodes are depicted in Figure 13. Because of the limitation of sensors, the modal shapes were incomplete, meanwhile, the structure is 3-dimensional, hence it is difficult to calculate the modal assurance criterion, thus, only the natural frequencies are adopted. The analytical and experimental natural frequencies are extracted and listed in Table 8.  As shown in Table 8, there exist some errors between the actual structure and the finite element model, especially for the higher-order modes, which can be owed to the environmental effects, instrument errors, and size deviations of the model. Thus, model updating is conducted using the CMF-IWO hybrid algorithm. After model updating, it can be revealed that the consistency between analytical and experimental models is really good, which means the model can be used for the baseline model of damage detection.  As shown in Table 8, there exist some errors between the actual structure and the finite element model, especially for the higher-order modes, which can be owed to the environmental effects, instrument errors, and size deviations of the model. Thus, model updating is conducted using the CMF-IWO hybrid algorithm. After model updating, it can be revealed that the consistency between analytical and experimental models is really good, which means the model can be used for the baseline model of damage detection.
Assuming the mass of the plate is uniformly distributed, that is to say, when the size is the same, the mass of each element is equal. By introducing the mass additional factors, the total mass can be calculated as: where M Mass represents the total weight of the plate and additional mass blocks; m i and β i are the i-th elemental mass and additional factor respectively. Hence, when β i is obtained, the weight of the additional mass block can be identified based on elemental mass. However, because of the limitation of sensors, the measured mode shapes are incomplete resulting in the fact that MSEBI cannot be obtained. For the purpose of bearings damage and additional mass detection, the location of mass blocks is assumed to be known, thus the search range is reduced. Additionally, the detection of mass and bearings damage are separated, namely, the mass change is first determined, after that the damage of the bearings is identified. According to Equation (42), the objective function can be defined as follows: where X = [β, α]. For the Cases 2-8, the detection procedure is carried out for seven times, the average identification results are extracted and depicted in Table 9. As listed in Table 9, the proposed method can precisely detect the separation of bearings. However, for the identification of mass change, the errors occur, namely, the inaccurate weight quantification. For Cases 3 and 6, the identified results are acceptable, but for other cases, the overestimating or underestimating has emerged, which may be attributed to the inaccurate measured data and the errors of the finite element model.

Conclusions
A novel bearings damage detection method using sensitivity analysis and chaotic moth-flameinvasive weed optimization hybrid algorithm has been put forward to determine the damage of structural elements and bearings. According to the obtained results, some conclusions and prospects can be summarized as follows: (1) The sensitivity coefficients of eigenvalues to the damage factors of structural elements and bearings provide a good evaluation approach to research the influences of damage of structural elements and bearings to eigenvalues, meanwhile, which is the basis for constructing the regularization objective function.
(2) MSEBI, as a damage location index, is able to accurately detect the damage location of structural elements as well as identify the locations of damaged bearings with high precision. At the same time, this damage location approach can greatly reduce the search range of damage detection and promote detection effectiveness. (3) Compared to PSO, CS, MFO, DE, and IWO, the proposed hybrid algorithm, CMF-IWO, is demonstrated with good convergence speed and global search performance, which is of good potential for overcoming the problem of local optimal in damage detection. (4) The proposed method is proved to be well applied in the bearings damage detection of numerical simulation, supported by the study case. Compared to the existing methods, the proposed method is easy and convenient to conduct and only the first few modal characteristics are needed. Hence complex calculation can be avoided. However, because of some uncertain factors and errors, such as the inaccuracy of instrument measurement, incomplete modal information, inevitable environmental noise, and the deviation of the finite element model, there still exist some difficulties to obtain very precise results in practical application.
Author Contributions: Conceptualization, methodology, resources, writing-review and editing, supervision and project administration must be acknowledged to M.H.; software, validation, formal analysis, investigation, data curation, writing-original draft preparation must be acknowledged to Y.L. All authors have read and agreed to the published version of the manuscript.
Funding: The paper was funded by the plan of outstanding young and middle-aged scientific and technological innovation team in universities of Hubei Province (Project No.T2020010).

Conflicts of Interest:
The authors declare no conflict of interest.