Complete Inverse Method Using Ant Colony Optimization Algorithm for Structural Parameters and Excitation Identification from Output Only Measurements

In vibration-based structural health monitoring of existing large civil structures, it is difficult, sometimes even impossible, to measure the actual excitation applied to structures.Therefore, an identification method using output-only measurements is crucial for the practical application of structural health monitoring. This paper integrates the ant colony optimization (ACO) algorithm into the framework of the complete inverse method to simultaneously identify unknown structural parameters and input time history using output-only measurements. The complete inverse method, which was previously suggested by the authors, converts physical or spatial information of the unknown input into the objective function of an optimization problem that can be solved by the ACO algorithm. ACO is a newly developed swarm computation method that has a very good performance in solving complex global continuous optimization problems.The principles and implementation procedure of the ACO algorithm are first introduced followed by an introduction of the framework of the complete inverse method. Construction of the objective function is then described in detail with an emphasis on the common situation wherein a limited number of actuators are installed on some key locations of the structure. Applicability and feasibility of the proposed method were validated by numerical examples and experimental results from a three-story building model.


Introduction
Structural health monitoring (SHM) has remained an active research topic in structural engineering since the 1970s.SHM identifies the occurrence, location, and severity of structural damage via significant adverse changes of structural parameters or properties.Thus, the core part of an SHM system is the algorithm used to accurately identify the structural parameters.A number of methods are available nowadays for us to accomplish this identification task.Housner [1] presented an extensive summary on the state-of-the-art methods in the vibration control and health monitoring of civil engineering structures.Zou et al. [2] summarized the methods of vibration-based damage detection and health monitoring for composite structures.Very recent reviews on identification methods used in SHM can be found in Ou and Li [3] and Fan and Qiao [4] among several others.Traditional identification algorithms are generally based on the assumption that the system's input (excitations) and output (responses) are completely known (measured).However, in vibration-based structural health monitoring of existing large civil structures, it is difficult, sometimes even impossible, to measure the actual excitation applied to the structures.Therefore, output-only structural parameter identification methods are of great significance for practical application of SHM.
In contrast to a number of publications on structural parameter identification with complete output and input information, there is a paucity of publications addressing the identification method using responses only.Since the input/excitation is unknown, assumptions on the input are necessary in order to make the traditional identification method applicable.The most commonly adopted assumption is treating the input as a white-noise process, whose power spectrum is theoretically known.The classical Ibrahim method plus random decrement technique, for instance, is based on this assumption.This assumption, however, is not tenable for excitation like earthquakes or strong winds that are in fact a nonstationary random process.Moreover, the time history of the real input cannot be directly identified.The second common way is assuming that the input has some special features.For instance, Toki et al. [5] assumed that the coda of the measured structural responses during an earthquake could be treated as free vibration responses from which the structural parameters could be easily identified.Wang and Haldar [6] identified the unknown earthquake input and structural parameters using output-only measurements through a recursive identification procedure consisting mainly of three steps.Step 1: The unknown structural parameters were identified by assuming that the input excitations were zero at the beginning building up to four time instants, say  0 to  4 .Step 2: The input excitation was conjectured by the measured structural responses and the parameters estimated in the first step.Step 3: Replace excitation at  0 to  4 in step 1 with new values obtained at step 2, and then repeat step 1 and 2 until the identified input excitation at  0 to  4 converged within a preset tolerance limit.For unknown wind load, Law et al. [7] assumed the wind load model was known even though the time histories of the wind load had not been recorded.Based on this assumption, the time history of wind load and the structural parameters could be identified using the response measurements.More recently, Yang et al. [8] suggested a recursive least squares estimation with unknown inputs to identify the stiffness, damping, and other nonlinear parameters at element level.The locations of the unknown excitations are assumed known in their approach.The ASCE structural damage benchmark structure was used to show the feasibility of the method.Yang and Huang [9] further extended this method to situations where external excitations and some acceleration responses are not measured.
To deal with the unknown input identification problem, we have proposed the concept of complete inverse problem, which means identification of structural parameters and the input's time histories simultaneously from the output-only measurements as per Chen and Li [10,11], Li and Chen [12,13], and Zhao et al. [14,15].Within the framework of complete inverse problem, we have suggested a series of identification methods named as complete inverse methods (CIM) addressing different types of unknown excitations.The core idea of CIM is to convert any additional information of the excitations (whose time histories are unknown) into mathematical constraint conditions that can be further integrated into an iteration identification procedure based on least square technique.For instance, when the locations of the excitations are known, the spatial information of the external excitation can be used in the parameter identification method [10,11].For ground motion excitation like that resulting from an earthquake, its mechanical features indicate that the inertial force proportional to the mass can be introduced into the interaction procedure as a mathematical constraint to ensure a stable and unique solution [12,13].For a proportional-type excitation, like wind loads, the ratio of forces at different structural heights can be used as a mathematical condition to identify the structural parameters and inputs.For shear-type building under earthquake excitation and with limited response measurements, we proposed a hybrid identification method where the unknown structural parameters for the first floor are identified using measured modal shapes, and parameters of all the other stories are identified using measured acceleration responses [14,15].We have also provided strict mathematical proofs from Chen and Li [11] and Li and Chen [13] to demonstrate the unconditional convergence features of the proposed complete identification methods.Despite CIM's success, it does have several limitations.This paper thus tries to improve the CIM in two directions.The first arises from replacing the least square method with another efficient and more robust optimization method.The recently emerged ant colony optimization algorithm (ACO) has been adopted in this paper to identify the structural parameter.The second improvement aims to validate the proposed method by experimental examples.The feasibility and effectiveness of all the aforementioned identification methods have already been demonstrated by different kinds of numerical examples.However, few experimental investigations have been conducted to assess the practical application of these methods.The effect of measurement noise and modeling error cannot be fully investigated in a numerical simulation.However, time-domain identification methods are known to be sensitive to measurement noise.
To this end, the paper first presents the principle of ACO for solving both discrete optimization problems and continuous optimization problems as well.Then, the ACO is integrated into the CIM methods by constructing the objective function according to the type of excitation.After that, the proposed method is applied to numerical model and an experimental model.The results show that the CIM+ACO algorithm performs very well for a noise-free and For all cases, identical initial values range [1, 1𝐸7] is used for all stiffness and damping parameters;  ⬦ : identified parameter vector; Error: relative error, : number of sampling points used.noise-polluted case and has good identification accuracy in parameters and inputs.

Ant Colony Optimization
Inspired by the ants' foraging behavior, Dorigo [16] proposed the ant colony optimization algorithm (ACO).ACO emulates the behavior of a group of ants in searching food from their nest to the food sources.Every ant leaves an amount of pheromone on the path that it passes and chooses the path with more pheromone left on it from the previous ants.Then, after more and more ants pass, the path having the maximum pheromone will be the best (shortest) way from the nest to the food source.ACO algorithms have already been successfully applied for solving combinatorial optimization problems, including the traveling salesman problem (TSP) [17], the routing problem in a computer network [18], the quadratic assignment problem (QAP) [19], and structural health monitoring problems [20,21].ACO algorithms for continuous optimization have been proposed in the literature [22][23][24].All the above work has proven ACO to be an efficient and versatile tool for solving various continuous optimization problems.The ACO algorithm has already been well established so far.Principal and application procedures of ACO are briefly summarized in Section 2.1.[22].As mentioned earlier, ACO emulates the behavior of a group of ants in searching for food.Modeling of the ant system and the pheromone left by ants on the path are of great importance for application of ACO.When ACO is used for discrete optimization problems, ants construct solutions incrementally.That is, each ant starts with an empty solution  0 and a component of the solution is added at each construction step.If   denotes all the available solutions at step , we need to choose one best solution from   and add the solution to previous solution  −1 leading to the new solution as   .The definition of the available solution componentdepends on the problem tackled.For instance, in the popular travelling salesman problem (TSP), a component of the solution is a city that is added to a tour.The solution components may be defined differently for other problems.

Ant System Model
A probability-based strategy is adopted to choose the best solution components from   for constructing the current partial solution   .This decision is usually influenced by amount of pheromone  associated with available choices and by heuristic information about the problem.To avoid the loss of generality, in the following problem, no heuristic information is used.Assuming that the partial solution constructed is   so far, the probability  is usedto choose a solution component  ∈   at step  in iteration  as calculated by the following; where  is the pheromone, the subscript  denotes the solution domain, and  means the th iteration step.Hence, in case of discrete optimization problems, the ants make a probabilistic decision according to some discrete probability distribution at each construction step.For continuous optimization problems, the domain changes from discrete to continuous.The logical adaptation also would be moved from using the discrete probability distribution to a continuous one-the probability density function (PDF).Instead of choosinga component  ∈   at step , the ants would generate a random number according to a certain PDF ().

ACO for Continuous Domain:
ACO R [22,23].The original ACO algorithm applies only to discrete domain and cannot be directly introduced into continuous domain optimization problems.Structural identification in nature is an optimization problem that aims to find the best parameters for a given objective function defined in a continuous domain.The ACO for continuous domain is, therefore, necessary.We adopted the method suggested by Socha [22] and Socha and Dorigo [23], which is an ACO for continuous domain (ACO  ) based on Gaussian probability density functions.Application of ACO  is briefly summarized here for the purpose of completion.More technical details of ACO  can be found in Socha [22] and Socha and Dorigo [23].
A Gaussian kernel PDF is used in ACO  to account for the multiple-peaks (multiple optimization solutions) domain.Suppose the optimization problem on hand is  = 1, 2, . . .,  dimensions; for the th dimension, the Gaussian kernel PDF   () is defined as a weighted sum of several one-dimensional Gaussian functions as follows: where,  = 1, 2, . . .,  is the number of dimensions of the problem;  = 1, 2, . . .,  is the number of single Gaussian functions constituting the Gaussian kernel PDF;   is the weight associated with the th individual Gaussian function;    and    are the mean and standard deviation for the th function in the th dimension.They can also be expressed in vector form as {  } and {  } whose cardinality is equal to .
For an -dimension problem, an ant constructs a solution in  steps.At each step, an ant gets a value for the unknown variable   .For each solution   , ACO  will store the current results of  as unknown variables and the value of the objective function (  ) in an archive .All the solutions are ordered in  according to their qualities, that is, ( 1 ) ≤ ( 2 ) ≤ ⋅ ⋅ ⋅ ≤ (  ) ≤ ⋅ ⋅ ⋅ ≤ (  ).Each solution has an associated weight  that is proportional to the solution quality,  1 ≥  2 ≥ ⋅ ⋅ ⋅ ≥   ≥ ⋅ ⋅ ⋅ ≥   .The PDF   is constructed using only the th coordinates of all  solutions from the archive.
For each dimension  = 1, 2, . . .,  of the problem, there is a different Gaussian kernel PDF   () defined.For each   (), the values of the th variable of all the solutions in the archive  become the elements of the vector   : Before each solution is added to the archive , it must be evaluated and ranked.The solutions in the archive are sorted according to their rank-that is, solution   has rank .Better solutions will have a higher weight.The weight   of the solution   is calculated by the following: The weight is for the Gaussian function with argument , mean 1.0, and standard deviation .In the ACO  ,  is in fact a parameter to balance between diversification and intensification.When  is small, the best-ranked solutions   are strongly preferred, and when it is large, the probability becomes uniform.
In practice, generating the Gaussian kernel PDF is accomplished as follows.First, choose one of the individual Gaussian functions that compose the Gaussian kernel with probability   given by (5).Then generate the chosen Gaussian function We consider the chosen Gaussian function given in (5) with    a standard solution for other ant solutions to explore.To establish the value of the standard deviation    at step , we calculate the average distance from    to all solutions in the archive , and multiply it by the parameter : where the parameter  > 0 has a similar effect to the pheromone evaporation rate in traditional ACO.The higher the value of , the lower the convergence speed of the algorithm.
At the start of the algorithm, the solution archive  is initialized generating  solutions by uniform random sampling.Pheromone update is accomplished by adding the set of newly generated solutions to the solution archive  and then removing the same number of the worst solutions so that the total size of the archive does not change.This process ensures that only the best solutions are kept in the archive, so that they effectively guide the ants in the search process.When the algorithm is completed, the solutions are ordered in the archive according to their quality, and the best solution is s 1 = { 1  1 ,  2 1 , . . .,   1 , . . .,   1 }.

Introduction of CIM.
The equation of motion of a -DOFs system can be expressed as follows: where M, C, and K represent, respectively, mass, damping, and stiffness matrices of the structures, X() and Ẋ(), while Ẍ() represents, respectively, the displacements, velocities, and accelerometers response vector of the structure; F() is the external excitation on the structures.Equation ( 7) can be rewritten as (8), which can be further rearranged as (9) at time instant   [10]: where vector  contains all the unknown parameters to be identified; matrix H consists of the measured structural responses; and vector Y is the system input.To assemble (9) at all the sampling time instants   ,  = 1, 2, . . .,  together, will give The components of matrix H and Y depend on the type of structure.For a shear-type structure, the expressions for H, Y, and  are where  1 ,  2 , . . .,   , and  1 ,  2 , . . .,   represent, respectively, the stiffness and damping coefficients of the structure for each story, and   (  ) and   (  ) are the displacement response and external excitation force of the th DOF ( = 1, 2, . . ., ) at the time instant   .
In traditional calculations, the parameters of the structure can be identified from ( 10) by the least-squares technique as However, ( 16) cannot be easily solved to determine the stiffness and damping parameters since there are unknown quantities involved in calculating H  and Y  .
As mentioned earlier, we have suggested the complete inverse method to tackle the unknown input situation in structural identification.Following the framework of CIM, Chen and Li [11] suggested the total compensation method.The total compensation method rests on the assumption that the locations of the external forces are known even though their time histories are unknown, and the number of DOF with applied (unknown) forces is less than the number of DOF whose responses are measured.This assumption reflects the situation of a forced vibration survey of structure where a limited number of one or several actuator(s) are installed on key locations of the structure to excite it.In this case, the input excitations in (15) can be further expressed into two parts: where   (  ) denotes those DOFs with unknown external excitation and   (  ) stands for those DOFs without applied force, that is,   (  ) ≡ 0.

Objective Function.
In order to identify structural parameters and the input time history from output-only measurements, an objective function is defined as minimizing discrepancies between the force   () and the calculated force  ⬦  ().The minimization of the objective function is expressed as a bound-constrained nonlinear least squares problem: min where  is the set of  and the identified force  ⬦  (  ) is given by where H  (  ) is th row of the matrix H(  ), and  ⬦ is the result of identified parameters based on ACO  with the details as where   1 ,  = 1, 2, . . ., 2 is the value of the first row of the archive , which is the best solution for (18).
Once the iterative procedure converges, the updated parameter vectors in (20) will give the final identification result of all the structural parameters, whilst the time history of the input F() can be easily determined by (7).

Numerical Examples
The suggested algorithm has been applied to several numerical examples including a truss structure, a 6-story shear frame, and a 12-story shear frame structure [25].Since the observations for all numerical examples are similar, only the results for the 6-story shear frame structure are presented here in detail.The mass, stiffness, and damping coefficient of each story of the structure are shown in Table 1.Sinusoidal excitations are applied on the 4th and the 6th floor, which are f 4 () = 6.74 × 10 5 sin(6) and f 6 () = 8.53 × 10 5 sin(4), respectively.Therefore the set  in ( 18) is  = {1, 2, 3, 5} for this example.The dynamic responses of all six stories are first calculated in terms of displacement, velocity, and acceleration using the Newmark- method [26].Twenty ants are used in each iteration of the ACO  algorithm, and the convergence threshold in the objective function is set as 1.0 × 10 −6 for all the cases.All the computation parameters used for ACO  are summarized in Table 2  identification results of the three cases are shown in Table 3, where Cases 1 and 2 have the same initial parameters' estimate but different measurement durations, and Cases 1 and 3 have the same measurement duration but different initial parameters.It is seen from Table 3 that, for all cases, the unknown parameters can be accurately identified by the proposed method with short duration of measurements and if the method is robust to obtain the initial estimated parameters.
The input excitation can also be accurately identified by the ACO  method in the noise-free case.

Noise-Pollution Measurements.
White noise is numerically added to the calculated responses to simulate noisy measured data by the following equation: where   is the noise level expressed as a percentage,  noise is a uniform distribution vector with interval [−1, 1], and Max( ẍ ) is the maximum value of the calculated acceleration response.Four different noise levels at 1%, 5%, 10%, and 15% were considered in the calculation.The identification results are summarized in Table 4, from which we can see that the proposed method can accurately identify the stiffness parameters for a noise-level up to 10%.The maximum parameter identification error of the stiffness parameter is lower than 0.1%, 0.6%, 1.6%, and 7.5% for noise levels 1%, 5%, 10%, and 15%, respectively.The identification accuracy for the damping ratio, however, is relatively low.This is not surprise since damping is very sensitive to measurement noise.
It is interesting to compare the inverse input time history of the top floor with that of the second floor where no external force was applied.The identified input forces of the two floors for different noise levels are shown in Figure 1, where the solid blue is the actual force curve and the dotted black line is the identified time history.Visually, even for noise level of 15% the identified force of the top floor matches well with the actual force.Amplitude of the identified curve of the second floor, on the other hand, is approximately zero compared with that of the top floor for noise level of 10%.This and several other numerical examples [25] have demonstrated the applicability and accuracy of the proposed identification method.

Experimental Example
The experimental data from a hammer test on a 3-story steel frame structure was adopted to further validate the applicability of the proposed identification method.The structural   The second is a continuous test wherein a certain floor of the building is hit continuously by the hammer.

Impact Test.
In this test, the building model was hit by a hammer once on the 3rd floor and was then released for free decayed vibration in the  direction.Figure 4(a) shows the recorded accelerations for each floor.The accelerations were then integrated to obtain the corresponding velocity and displacement as shown in Figure 5.To reduce the measurement noise, the high-pass Butterworth signal filter was applied to the integration with the filter order  = 4 where the lower cut-off frequency = 0.16 rad/sec.The known information for this case is  1 () ≡ 0.0 and  2 () ≡ 0.0, and the structural parameters are to be identified by the proposed ACO  method with the following objective function: Tables 5(a), 5(b), and 5(c) show, respectively, the identified stiffness parameters from ten different response segments.Each has the same duration of 0.2, 0.33 and 1.0 second (i.e., 60, 100, and 300 data points for a sampling frequency 300 Hz).The same initial values for all the unknown parameters were used in the calculation.The maximum average identification error for all three cases is less than 2%, and the maximum identification error for a single parameter is less than 4%.Comparison between Tables 5(a), 5(b), and 5(c) also demonstrates that the identification accuracy is not sensitive to the duration of the response segment.Table 5(d) shows the identification results for response duration of 0.33 second (100 data points) and different initial values of unknown parameters.Comparison between Tables 5(b) and 5(d) shows that the proposed method is not sensitive to the initial estimates of the unknown parameters.All the above results indicate a very good identification accuracy achieved by the proposed method for a very short response duration, which has potential application for online structural health monitoring.It should be noted that the damping identification results are not satisfied in this case.The damping properties can be identified by other effective identification techniques such as the empirical mode decomposition plus Hilbert transform method suggested by the authors [27].Figures 6 and 7 show the convergence procedure of all the stiffness parameters for the computational segments as 10.2∼10.53s and 11.0∼11.33s in Table 5(b).Note that the unknown stiffness parameters can be accurately identified using a short duration of response and the results are robust to the parameters' initial guess.
Figure 8(a) shows the time history of the actual hammer force on the 3rd floor, and Figures 8(b) to 8(c) show the identified time histories of forces for the 3rd, 1st, and 2nd floor, respectively.Note that the identified input on the 3rd floor has a spike at the same time instant as the real input.The amplitudes of the identified inputs on the 1st and 2nd floor are almost zero compared to that of the 3rd floor.

Continuous Hammer Test.
In this test case, the hammer hit the 3rd floor continuously, and the resulting accelerations for each floor are shown in Figure 9.The actual test duration was 60 seconds, but only the first 20 seconds were plotted for the sake of clarity.Using the same objective function as (22), the structural parameters were identified by the ACO  method, and the results are summarized in Table 6.
Note that the average identification error for each of the seven response segments is less than 2%, and the maximum identification error for each single parameter is less than 3%.Using the global average value of all the identified parameters, the inputs of each floor can be identified.Figures 10(a) and 10(b) compare the identified input and real input on the 3rd floor.It is clear that the identified input has the continued impact spikes (peaks) at the same time as instants for the real input.The amplitudes of the identified inputs of the 2nd and 1st floor are almost zero compared to that of the 3rd floor.The result is consistent with the real test situation.

Push the Building for a While and Release It to Free Vibration.
In this test, the model was hand-shaked on the 2nd floor to vibrate for a while and was then released for free vibration.Figure 11 shows the acceleration responses of each floor for this case.The known information is  1 () ≡ 0.0 and  3 () ≡ 0.0, where the objective function can then be established as follows: ( into a mathematical confinement condition.In theory, the proposed method has no limitation on the type of excitation.This paper applies the method to the situation where spatial locations of a limited number of excitations are known.Numerical examples were carried out to evaluate the feasibility of the proposed method.For the situation of noisefree output measurements, numerical studies show that the proposed method can reliably and efficiently identify both the structural parameters and the input time history using a short duration of measurements.Moreover, the accuracy and convergence of the method are robust to the initial values selected for the unknown parameters.For  noise-pollution measurements, the stiffness parameters as well as the input excitation were identified satisfactorily even with high noise level.The identification accuracy of the damping coefficients is broadly acceptable at low noise level and become relatively poor at high noise level.The proposed method is then applied to experimental results of a three-story building model and the results also prove the feasibility of the method.

Figure 2 :
Figure 2: The 3-story building model in the test.

Figure 3 :Figure 4 :
Figure 3: Dimensions of the building model.
Velocity of the third floor (b) Velocity of the second floor (c) Velocity of the first floor (e) Displacement of the second floor (f) Displacement of the first floor (d) Displacement of the third floor

Figure 5 :
Figure 5: The integrated velocity and displacement from accelerations shown in Figure 4.

4. 1 .Figure 8 :
Figure 8: The real input on the top floor and the identified inputs of each floor.
Acceleration response of the first floor

Figure 9 :Figure 10 :
Figure 9: The recorded acceleration responses of each floor for continuous hammer test.

Figure 11 :
Figure 11: The recorded acceleration responses of each floor for hand-shaked test.

Table 2 :
Parameters used by ACO  for HIM.

Table 3 :
Identified results with noise-free measurements.

Table 4 :
Identified results with noise-pollution measurements.

Table 5 :
(a) Identified results for single hammer hit on the top floor (computation points = 60).(b) Identified results for single hammer hit on the top floor (computational points = 100).(c) Identified results for single hammer hit on the top floor (computation points = 300).(d) Identified results for single hammer hit on the top floor (computation points = 100 with different initial values for parameters).

Table 6 :
Identified results for continuous hammer hits on the top floor.