Abstract
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.
1. 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 to . 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 to in step 1 with new values obtained at step 2, and then repeat step 1 and 2 until the identified input excitation at to 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 noise-polluted case and has good identification accuracy in parameters and inputs.
2. 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–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.
2.1. Ant System Model [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 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 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.
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 .
2.2. ACO for Continuous Domain: [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 () based on Gaussian probability density functions. Application of is briefly summarized here for the purpose of completion. More technical details of can be found in Socha [22] and Socha and Dorigo [23].
A Gaussian kernel PDF is used in to account for the multiple-peaks (multiple optimization solutions) domain. Suppose the optimization problem on hand is dimensions; for the th dimension, the Gaussian kernel PDF is defined as a weighted sum of several one-dimensional Gaussian functions as follows: where, is the number of dimensions of the problem; 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 , 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, . Each solution has an associated weight that is proportional to the solution quality, . The PDF is constructed using only the th coordinates of all solutions from the archive.
For each dimension 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 , 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 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 .
3. Integration of with CIM
3.1. Introduction of CIM
The equation of motion of a -DOFs system can be expressed as follows: where , , and represent, respectively, mass, damping, and stiffness matrices of the structures, and , while represents, respectively, the displacements, velocities, and accelerometers response vector of the structure; 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 consists of the measured structural responses; and vector is the system input. To assemble (9) at all the sampling time instants , together, will give
The components of matrix and depend on the type of structure. For a shear-type structure, the expressions for , , and are
where , and 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 () 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 and .
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, .
3.2. 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: where is the set of and the identified force is given by where is th row of the matrix , and is the result of identified parameters based on with the details as where , 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 can be easily determined by (7).
4. 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 and , respectively. Therefore the set in (18) is 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 algorithm, and the convergence threshold in the objective function is set as for all the cases. All the computation parameters used for are summarized in Table 2.
4.1. Noise-Free Measurements
The proposed method is first applied to noise-free measurements. The parameter 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 method in the noise-free case.
4.2. 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, is a uniform distribution vector with interval , and 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.
(a) Noise level of 1%
(b) Noise level of 5%
(c) Noise level of 10%
(d) Noise level of 15%
5. 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 model is shown in Figure 2 and dimensions of the model are given in Figure 3. The strength of the steel of the model was tested as 435 Mpa, and the elastic module was 200 Gpa. The overall dimension of the structure is mm (Figure 3). The structure’s column was made of a 9.5 × 75 mm steel plate and the floor was made by an 850 × 500 × 25 mm steel plate. Additional mass of 130 kg was installed on each floor leading a total floor mass to 230.3 Kg. The structural model can be treated as a shear-type building, and the shear stiffness of each floor was computed as 563651 N/m. Hammer tests, hand-shaking tests, and shaking table tests were conducted on the model. Only hammer tests and hand-shaking tests data are adopted in this paper. The hammer tests can be broadly divided into two categories. The first is an impact test where a certain floor of the building is hit once by a force hammer and then the building vibrates freely. The second is a continuous test wherein a certain floor of the building is hit continuously by the hammer.
(a) Elevation view
(b) Plan view
5.1. 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 where the lower cut-off frequency = 0.16 rad/sec.
(a) The third floor
(b) The second floor
(c) The first floor
The known information for this case is and , and the structural parameters are to be identified by the proposed 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.53 s and 11.0~11.33 s 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.
| (a) Stiffness |
| (b) Stiffness |
| (c) Stiffness |
| (a) Stiffness |
| (b) Stiffness |
| (c) Stiffness |
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.
(a) Real input on the third floor
(b) Identified input on the third floor
(c) Identified input of the second floor
(d) Identified input of the first floor
5.2. 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 method, and the results are summarized in Table 6.
(a) Acceleration response of the third floor
(b) Acceleration response of the second floor
(c) Acceleration response of the first floor
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.
(a) Real input on the third floor
(b) Identified input on the third floor
(c) Identified input of the second floor
(d) Identified input of the first floor
5.3. 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 and , where the objective function can then be established as follows:
(a) Acceleration response of the third floor
(b) Acceleration response of the second floor
(c) Acceleration response of the first floor
Table 7 shows the identification results for this case for seven different response segments. The maximum identification error is less than 1%. It is interesting to compare the identified results of input with the real input. Figures 12(a) to 12(c) show the identified inputs of 1st to 3rd floors, respectively. It is interesting to find that the input of the 2nd floor (Figure 12(b)) has a relatively large excitation force for the duration of 12 to 22 seconds. That is consistent with the test at that duration where the model was hand-shaked. The amplitudes of the inputs from the 1st and 3rd floor are relatively small compared to that of the 2nd floor. The above three cases demonstrate the applicability and feasibility of the proposed hybrid identification method using response-only measurements.
(a) Identified input on the third floor
(b) Identified input on the second floor
(c) Identified input on the third floor
6. Conclusions
This paper combines the complete inverse method with the ant colony optimization algorithm to identify the input excitation and the structural parameters from the output-only measurements. The core idea of the complete inverse method is to convert the physical features of the excitation 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 noise-free 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
Professor Jun Chen gratefully acknowledges financial support for the research represented by this paper, which came from The Project of National Key Technology R&D Program in the 12th Five Year Plan of China (2012BAJ11B02).