A Time-Domain Structural Damage Detection Method Based on Improved Multiparticle Swarm Coevolution Optimization Algorithm

Optimization techniques have been applied to structural health monitoring and damage detection of civil infrastructures for two decades.The standard particle swarm optimization (PSO) is easy to fall into the local optimum and such deficiency also exists in the multiparticle swarm coevolution optimization (MPSCO). This paper presents an improved MPSCO algorithm (IMPSCO) firstly and then integrates it with Newmark’s algorithm to localize and quantify the structural damage by using the damage threshold proposed. To validate the proposed method, a numerical simulation and an experimental study of a seven-story steel frame were employed finally, and a comparison was made between the proposed method and the genetic algorithm (GA). The results show threefold: (1) the proposed method not only is capable of localization and quantification of damage, but also has good noisetolerance; (2) the damage location can be accurately detected using the damage threshold proposed in this paper; and (3) compared with the GA, the IMPSCO algorithm is more efficient and accurate for damage detection problems in general. This implies that the proposed method is applicable and effective in the community of damage detection and structural health monitoring.


Introduction
Structural health monitoring (SHM) systems offer the possibility to detect the damage for large civil infrastructures accurately and immediately so as to ensure structural integrity and safety.Among the techniques adopted for SHM, vibrationbased structural damage detection (SDD) has attracted considerable attention for two decades [1][2][3].However, with an increasing amount of vibration data collected from the SHM system, it is still challenging to develop more efficient and robust detection algorithms in respect of variation in structural parameters due to deterioration or damage.
In particular, time-domain analysis techniques have been studied for a relatively long time and proven to be useful for SDD.Mathematically, SDD is a highly nonlinear inverse problem and for time-domain identification it is usually solved by minimizing an objective function, which is defined in terms of the discrepancies between time-series responses measured from the laboratory or in-site and those computed from the analytical model.However, conventional optimization methods, such as the conjugate gradient, usually lead to a local optimum.Hence, a global optimization technique is needed to develop so as to obtain a more accurate and reliable solution [4].In recent years, many computational intelligence methods have been proposed and applied to SDD, such as genetic algorithm (GA) [5][6][7][8], artificial neural network (ANN) [9,10], and swarm intelligence techniques [11,12].Particle swarm optimization (PSO) proposed in 1995 [13] is a popular swarm intelligence method, which mimics the collective motion of insects and birds, known as "swarm behavior." PSO has achieved to date tremendous progress and has been successfully applied to SDD due to its simplicity, fast convergence speed, and outstanding performance.
Similar to other evolutionary algorithms, however, standard PSO also has the problem of premature convergence and can be trapped into some local optima.As a result, a number of improved PSO methods have been developed [14][15][16].In recent years, the multiparticle swarm coevolution optimization (MPSCO) by integrating the collaborative theory in ecology with the principle of automatic adjustment has become a popular hotspot [17][18][19].Nevertheless, such efforts cannot 2 Mathematical Problems in Engineering completely solve the problem of local optimum.Therefore, an improved MPSCO algorithm called IMPSCO is proposed and applied to SDD in this study.In IMPSCO the evolutionary theory is integrated with MPSCO so as to reduce the possibility of falling into the local optimum.
To develop an effective time-domain SDD method, this paper proposes an optimization-based damage strategy by integrating IMPSCO with Newmark's algorithm, and a damage threshold is also developed to accurately quantify damage.Numerical simulations and experimental test of a 7floor steel frame are conducted to validate the efficiency of the proposed approach; furthermore, the proposed method is also compared with other methods like GA.Both simulation and experimental results have demonstrated that (1) the proposed method uses any kinds of structural time-series responses and only requires very few sensors in practical applications; and (2) after the damage threshold is determined by hypothesis testing, damage location and damage extent can be identified quickly and accurately with the enormous data polluted by measurement noise which is induced by the complicated operating environment of structures.
The organization of the paper is as follows.Section 2 provides the basic theory of IMPSCO.Section 3 describes the damage detection strategy based on IMPSCO and Newmark's algorithm.Numerical simulations are carried out in Section 4, while experimental study is presented in Section 5. Section 6 gives the concluding remarks.

Improved Multiparticle Swarm Coevolution Optimization
Algorithm.With basic MPSCO, multisubpopulations are divided into two layers.All particles from the upper-layer follow the optimum of the entire population so as to obtain a faster convergence speed, while all particles from the lowerlayer follow the optimum of the subpopulation to which it belongs, so as to ensure the population diversity.Although the performance of basic MPSCO is better than standard PSO in some aspects, the subpopulations in lower-layer still perform the process of standard PSO, which makes its falling into the local optimum possible.To solve this problem, the IMPSCO algorithm is proposed and applied to locate and quantify damage of a structure in this study.
In nature, some species disappear because of environmental variation, while some new species will emerge at the same time so that the species diversity is balanced.Accordingly, the searching procedure of PSO can be regarded as a natural selection process of species in biology.If a particle is recorded as the worst particle for many times in the iteration of the algorithm, it indicates that this particle is unable to adapt to the current environment and needs to be eliminated.As a consequence, the improved strategy is presented, in which the particles to be deleted are replaced with the gravity position of the selected excellent particles in current entire population.After the replacement, the particles can quickly get out of the local optimum.All in all, the proposed IMPSCO algorithm is described as follows.
Step 1 (population initialization).Generate  subpopulations randomly, and each of them contains  particles.Then divide them into the upper-layer with only one subpopulation and the lower-layer with  − 1 subpopulations.Meanwhile, set the iteration index  to zero.
Step 2 (fitness calculation).Evaluate the fitness of each particle and save the personal best and its corresponding fitness as pbest and pbestval.Simultaneously, record the best individual of each subpopulation and the corresponding fitness as gbest j and gbestval j ( = 1, 2 . . .); meanwhile record the best individual in the entire population and its corresponding fitness as gbest and gbestval.
Step 3 (particles updating).Update the particles in the upper layer and the lower layer according to (1) and ( 2), respectively, and the worst particle in the entire population is recorded: where  is the particle's index in the swarm;   and V  represent the position and velocity of the particle, respectively;   represents the optimal position of the particle;   and   represent the optimal position of the entire population and the subpopulation to which the particle belongs, respectively;  1 and  2 are the random numbers between zero and one;  1 and  2 are the learning factors;  is the inertia weight.
Additionally, the maximum velocity of each particle cannot exceed  max which is set to be 20% of the length of the search space.
Step 4 (optimum updating).Calculate the fitness of each updated particle and compare it with the values above.If it is preferable, then update pbest, pbestval, gbest, gbestval, gbest j and gbestval j , correspondingly.Let  =  + 1.
Step 5 (worst particle replacement).Repeat Steps 3-4.When the particle is recorded as the worst for the predetermined times  w , replace it with  g as shown in (3), and its worst record returns to zero: where  and   represent the number of the selected excellent particles and their position, respectively.
Step 6. Go to Step 3, and repeat until the maximum iteration times  max is reached.

The Newmark Integration Method.
The Newmark method is a commonly used numerical integration method to solve differential equations and has been widely applied to numerical evaluation of the dynamic responses for structures.It is based on the following assumptions: where , Ẋ, and Ẍ are the displacement, velocity, and acceleration, respectively; Δ is the time step;  and  are the parameters which determine the accuracy and stability, respectively.When  = 1/2 and  = 1/4, it is Newmark's constant-average acceleration method.
The equilibrium equation at time  + 1 is also represented with where M is the mass matrix; C is the damping matrix; K is the stiffness matrix;  +1 is the external force at time  + 1.After  0 , Ẋ0 , and Ẍ0 are initialized; the structural responses , Ẋ, and Ẍ at each time step can be calculated using ( 4)-(6).

Rayleigh Damping.
In this study, the damping matrix C is obtained according to the assumption of Rayleigh damping model: where  and  represent the mass and stiffness damping coefficient, respectively, which can be calculated by where   and   are the natural frequencies of the th mode and th mode;   and   are the corresponding damping ratios.

Damage Detection Strategy
On the basis of the abovementioned, this section proposes a time-domain damage detection strategy by integrating the IMPSCO and Newmark's algorithm.The schematic diagram of the damage detection strategy is depicted in Figure 1.
From Figure 1, it is found that two steps are significant for the damage detection strategy, namely, Determination of Damage Threshold and Parameter Identification with IMP-SCO and Newmark's algorithm [20].Generally, it involves parameters encoding, establishment of fitness function, and parameters setting for the IMPSCO algorithm.

Parameters Encoding.
The first task is to encode parameters involving the IMPSCO.For a frame, generally, the floor stiffness is encoded directly as the parameter to be optimized in the intact state.However the actual stiffness values cannot be obtained accurately, so we turn to focus on detecting the stiffness reduction for damage scenarios.Therefore the floor stiffness ratio of the damaged structure to the undamaged structure, which is defined as the rigidity coefficient,is introduced and encoded in order to make the detection results simple and clear.In addition, with the Rayleigh damping taken into consideration, the mass and stiffness damping coefficients  and  should also be encoded.

Parameters Setting in IMPSCO.
To implement optimization, some parameters should be set and initialized in advance for the IMPSCO.Usually, 3∼5 subpopulations having 30 to 100 particles in total are sufficient, and the maximum iteration times  max can be set to 50∼500.When the inertia weight  is set linearly to a range from 0.9 to 0.4 in variation and learning factors  1 and  2 are set to 2 simultaneously, the PSO algorithm exhibits excellent performance.In addition, the limited times for the worst record  w and the number of the selected excellent particles  are both empirically set to 5∼10.
In this study, the parameters are set as follows: the total number of subpopulations  = 3; each population size  = 10; learning factors  1 =  2 = 2; limited times for the worst record  w = 5; maximum iteration times  max = 60; the number of the selected excellent particles  = 6; the inertia weight  is linearly decreased from 0.9 to 0.4 before the 45th iteration and afterwards it maintains at 0.4 to enhance the local search capability.

Establishment of Fitness Function.
The most important task is to determine the fitness function for the IMPSCO.Firstly, the structural mass matrix M, stiffness matrix K, and damping matrix C are constructed with the generation of particles.And then the simulated time-series responses can be obtained by using Newmark's constant-average acceleration method.Only if the analytical responses are close to the measured ones can it be determined that the structural properties represented with the particles agree well with actual damage situations.Consequently, the fitness function can be represented with where  is the parameters vector;  is the number of measuring points;  is the length of response data;  mea , and  sim , are the measured and simulated time-series responses, respectively.
When () reaches the maximum, the values of  are the optimal solution.

Determination of Damage Threshold and Results Output.
It is a difficult and significant task to determine the damage threshold because the threshold is usually used to judge damage.If the identification results are higher than the threshold, it indicates that the damage occurs; otherwise the damage is excluded because of the inevitable measurement noise.Eventually the damage detection is performed.
In this study, the damage threshold is set to 2% validated by the statistical hypothesis testing, which is a method for testing a hypothesis regarding a population parameter, such as the mean from measured sampling data.It usually involves four steps.Firstly, extract the frequency sample set of the structure before and after 2% damage, representing with X 1 and X 2 , respectively; secondly, design the null hypothesis (H 0 ) :  1 =  2 , in which  1 and  2 represent the means of X 1 and X 2 , respectively; thirdly, construct a test statistic using the obtained frequency samples and determine the rejection region on the basis of the predetermined significance level ; finally, decide whether to accept H 0 or not by comparing the value of the test statistic and the rejection region.If accept H 0 , it indicates that the mean of the frequency samplings does not change significantly when 2% damage is introduced, implying that the floor with identified damage extent below 2% is intact for the structure.More details can be seen in Section 4.2 (2)., 7) is 375 kN/m, which is simulated using the 1-D COMBIN14 spring element.The Rayleigh damping is employed with the damping ratio of 2%.Then the natural frequencies of the first two modes  1 and  2 can be obtained by modal analysis and therefore the damping coefficients ,  can be calculated by (8).The model is subjected to random excitation at the seventh floor.

Numerical Study
(2) Damage scenarios: damage is simulated by reduction of the stiffness for each floor via where  represents the damage magnitude;   and   represent the stiffness in the th floor of the undamaged structure and damaged structure, respectively.Two magnitudes of damage and three damage locations are used to validate the proposed approach.The damage magnitudes are classified quantitatively as small (4.1%) and large (16.7%), while the damage locations are denoted by their corresponding floor numbers.Thus six damage patterns are simulated (as shown in Table 1) and are discussed further later.
(3) Data Acquisition.The random force is input to the numerical model based on the vibration exciter of the experiment in Section 5.2.By setting the time step Δ = 0.0002 s and the loading time  = 0.1 s, the structural responses can be calculated, and 500 acceleration data points can be collected at each floor for the structure in the intact and damaged cases, respectively.
In order to simulate the operating environment, 10% Gaussian white noise is added to all the acceleration responses using (11): where   and    represent the contaminated and theoretical accelerations, respectively;  is a normally distributed random variable with zero mean and a derivation of 1;  is an index representing the noise level, which is set to 10% here.Because the IMPSCO algorithm is a probabilistic optimization algorithm, the evolutionary process is always accompanied with randomness.Consequently, it is difficult to judge whether the detection result is correct from a single test.In this study the detection results with bad fitness values are eliminated and ten results are kept in total.Finally the average values are regarded as the final results, as presented in Table 2.
It can be seen that the identification results agree well with the theoretical values, and then the rigidity coefficients can be calculated and introduced to identify the damaged structure by using the identification results.In addition, it is seen that there is a significant difference between the identified damping coefficients and the theoretical values.Since the damping of the structure is reasonably small and only a short time-history is required, this damping is unlikely to have a significant effect on the simulated structural responses when calculating the fitness function.Therefore, the identification results of floor stiffness are considered reliable, which indicates that the accurate assumption of structural damping model is not necessary and the assumption of Rayleigh damping model meets the requirements.
(2) The hypothesis testing of damage threshold: as presented in Section 3.4, the damage threshold is set to 2% in this paper.In order to prove its reasonability, the natural frequencies of the first 7 modes are regarded as the sensitivity indices to carry out the hypothesis testing.The specific process is as follows.
Step 1. Extract the natural frequencies of the undamaged structure using the acceleration responses of each floor by Fast Fourier Transform (FFT).Take the average values as the final natural frequencies of the first 7 modes in the intact state.
Step 2. Repeat Step 1 for 30 times to obtain 30 frequency samples of the undamaged structure.The statistical theory shows that the samples follow the Gaussian distributions approximately if the sample size is large enough.
Step 3. Introduce 2% damage to any floor of the structure (the 4th floor is chosen here) and then calculate the corresponding acceleration responses by using the developed finite element model.Similarly, repeat Step 2 to obtain 30 frequency samples of the damaged structure with 2% damage.
Step 4. Employ the natural frequencies of the second mode for hypothesis testing because their variation is the maximum for all seven modes.Therefore, the frequency sample set of the undamaged structure is represented with X 1 = {244.1,240.9, 241.Step 5. Assume that 2 ), and ), in which  1 ,  2 , and  2 1 ,  2 2 represent the mean and the variance of X 1 and X 2 , respectively.
Step 6. Assume H 0 :  = 0, H 1 :  ̸ = 0. Then the -statistic is constructed and represented with where ,  2 , and  represent the mean, the variance, and the size of the samples, respectively, so the -statistic follows the -distribution with the degree of freedom of  − 1; that is,  ∼  ( − 1).
Step 7. In this example, the confidence level (1 − ) is set to 0.99, then the critical value  /2 = 2.76, and the rejection region can be represented with  = {|| ⩾ 2.76} on account of the -distribution.Meanwhile, the -statistic value can be obtained from (12), and the value is 1.54.It can be seen that the -statistic value is out of the rejection region; that is,  = 1.54 ∉ , so accept H 0 and reject H 1 It indicates that no significant change of the frequency has occurred before and after 2% damage, implying intact structure.
As a consequence, the damage threshold of 2% is reasonable from the above demonstration.
(3) Identification of damaged structure: the identification procedure is similar to that of undamaged structure, except that the parameters encoded are replaced by  =  3.
It can be seen that although 10% Gaussian white noise is added to the accelerations, the maximum relative error is 1.6% for all damage scenarios.This implies that structural damage can be detected effectively and accurately even in the case of a noisy environment.

Comparison and Discussion.
To validate the efficiency and robustness, a comparison was made between the theoretical damage extent, identification results by the genetic algorithm (GA), and IMPSCO.It is noted that the number of species group is 80 and the number of generation is 100 in the GA employed in the literature [20].Other parameters in GA are the same as IMPSCO.The identification results are depicted in Figure 3.
According to the damage threshold determined in Section 4.2 (2), the structure whose damage extent is lower than 2% is regarded as intact structure.Consequently the following observation can be found.
(1) The damage locations of all scenarios are detected correctly using two different optimization techniques, even for multidamage patterns like D4, D5, and D6.
(2) Although the detection of damage extent is also reasonable and reliable using GA, IMPSCO is better than GA in terms of identification accuracy in general.It is noted that GA is better than IMPSCO in the detection of damage extent for D2 and D4 at the cost of longer running time.
(3) Regarding the running time, IMPSCO costs 100 CPU seconds to implement the damage detection only half of the GA.

Steel Frame.
A seven-floor steel frame of 1.4125 m in height is constructed and tested in the laboratory [21].The model is designed with flexible columns and relatively rigid beams to simulate a shear building, as shown in Figure 4.The cross-section properties of the structural members are listed in Table 4.The mass of the structure is concentrated at the floors, and the structural model is regarded as a lumped mass model.Therefore, the steel frame can be simplified as a 7-DOF spring-mass system with the masses  1 ,  2 , . . . 6 = 3.78 kg and  7 = 3.31 kg.

Dynamic Test.
(1) Damage scenarios: as the structure is constructed with six columns per floor, damage is simulated by cutting the center column partially or completely in order to keep the symmetry of the model.Small damage is produced by four partial cuts near the top and bottom of the column, whereas large damage is simulated by a complete cut at the mid-height of the column, as indicated in Figure 5.The expected reduction of floor stiffness due to the small damage is estimated by the software of ABAQUS.The finite element models of both the undamaged and damaged columns are established, and the displacements under the same nodal loads are calculated and compared in order to determine the  change in column stiffness.As there are 6 columns in each floor, the small damage results in a reduction in the floor stiffness with 4.1%.As for the case of large damage, the column damage is 100% resulting in a reduction in floor stiffness of 1/6 ≈ 16.7%.More details can be seen in the literature [21].The six damage scenarios are shown in Table 1.
(2) Experimental data acquisition: Figure 6 illustrates the dynamic testing and data acquisition system.For ease of setup, the model is mounted horizontally and excited by vibration exciter vertically at the seventh floor.The force generated is measured by an ICP (Integrated Circuit Piezoelectric) force sensor (model PCB-208C02).Meanwhile, the structural responses are measured using the ICP accelerometers mounted at each floor.The data is recorded at a sampling frequency of 5000 Hz.

Damage Detection.
(1) Identification of undamaged structure: firstly, 500 data points of the acceleration responses at each floor and the corresponding force are extracted by the installed sensors when the structure is intact.Then follow the steps presented in Section 4.2 (1).Consequently the floor stiffness is identified for the intact structure.
The detection process is repeated for 15 times using the same method in order to improve the identification precision.The average values for the 15 times are regarded as the final results, as shown in Table 5.It can be seen from Table 5 that the identification result of  1 is significantly lower, and this is due to the less rigid connection at the base of the structure or the modeling errors of the initial numerical model.In addition, the identification results of  and  are 0.75102223 and 0.00000286, respectively.Then the acceleration responses can be calculated by the identified results according to (4)-( 6), which is plotted along with the measured values in Figure 7.For the simplicity, only the comparison results of the second and fifth floors are given as shown in Figure 7.It indicates that the identified and measured acceleration responses are in good agreement, which validates the reasonability of applying the identified floor stiffness of the undamaged structure to the detection of the rigidity coefficients under the different damage scenarios.
(2) Identification of damaged structure: the identification procedure is similar to that of intact structure in Section 5.3 (1), except that the parameters encoded are replaced with  = [ 1 ,  2 , . . . 7 , , ] as Section 4.2 (3).The final identification results are listed in Table 6.
It is seen from Table 6 that the maximum relative error of identification results is 1.7%, which occurs in the third floor for the multidamage pattern D4.As the rigidity coefficient introduced in the parameter encoding reflects the structural stiffness reduction, the effect on the damage identification results owing to the modeling error of the initial numerical model can be reduced to a great extent.This implies that the proposed structural damage detection strategy is effective and reliable.8. On the basis of damage threshold determined in Section 4.2 (2), the floor whose damage extent is lower than 2% is intact.It is found from Figure 7 that damage locations canbe correctly localized using IMPSCO and GA for all damage scenarios.In addition, the detection of damage extent is also accurate and reliable using the two optimization methods.Furthermore, IMPSCO is more efficient and accurate than GA in general.Just as the numerical example, although GA is capable of detecting damage, it costs more running time than IMPSCO.All in all, the proposed damage detection strategy can not only localize the damage correctly but also quantify damage precisely.

Conclusions
This paper proposes a new method by combining the IMP-SCO algorithm and Newmark's constant-average acceleration for the structural damage detection.From the damage detection results of a seven-floor frame model, some conclusions can be drawn as follows.
(1) The proposed damage detection strategy is applicable and effective for detecting and quantifying damages using the noise polluted measured data.It is noted that the proposed strategy is implemented only by using any kinds of structural time-series responses and the excitation force.
(2) Both the numerical simulation and laboratory test of a 7-floor steel frame have proven the feasibility and efficiency of the damage threshold proposed in damage localization.
(3) On the basis of the damage threshold developed by hypothesis testing, the proposed IMPSCO is capable of accurately detecting damage and reducing the possibility of falling into the local optimum owing to the fact of integrating the evolutionary theory with MPSCO.Compared with the GA, the proposed damage detection strategy by integrating the IMPSCO and Newmark's algorithm is better than GA in accuracy and runtime.
In summary, the proposed damage detection strategy-based IMPSCO is effective and robust for the numerical simulation and laboratory test of a 7-floor steel frame.However, more experiments are to be conducted to further validate its efficiency.

Figure 1 :
Figure 1: The procedure of damage detection.

4. 1 .
Numerical Model and Analysis.(1) Numerical model: in order to validate the effectiveness of the proposed strategy, a numerical model of a 7-floor shear frame is simplified as a 7-DOF Mass-Spring-Damper system and established by the software of ANSYS, as shown in Figure 2. The mass of each floor  1 = ⋅ ⋅ ⋅ =  6 = 3.78 kg and  7 = 3.31 kg, which is simulated using mass element MASS21.The floor stiffness   ( = 1, 2, . . .

Figure 3 :
Figure 3: The detection results of all scenarios.

Figure 5 :Figure 6 :
Figure 5: Damage applied to the column.

Figure 7 :
Figure 7: The comparison results of acceleration responses.

Figure 8 :
Figure 8: Identification results of all damage scenarios.

Table 2 :
Identification results in intact state.

Table 4 :
Section properties of the frame members.

Table 5 :
Identification results of floor stiffness in intact state.

Table 6 :
Identification results of all damage scenarios.Comparison and Discussion.Just as Section 4.3, a comparison is also made among the theoretical damage extent, identification results by genetic algorithm (GA), and IMP-SCO, and the results are shown in Figure