Structural System Identification in the Time Domain using Evolutionary and Behaviorally Inspired Algorithms and their Hybrids

In this study, parametric identification of structural properties such as stiffness and damping is carried out using acceleration responses in the time domain. The process consists of minimizing the difference between the experimentally measured and theoretically predicted acceleration responses. The unknown parameters of certain numerical models, viz., a ten degree of freedom lumped mass system, a nine member truss and a non-uniform simply supported beam are thus identified. Evolutionary and behaviorally inspired optimization algorithms are used for minimization operations. The performance of their hybrid combinations is also investigated. Genetic Algorithm (GA) is a well known evolutionary algorithm used in system identification. Recently Particle Swarm Optimization (PSO), a behaviorally inspired algorithm, has emerged as a strong contender to GA in speed and accuracy. The discrete Ant Colony Optimization (ACO) method is yet another behaviorally inspired method studied here. The performance (speed and accuracy) of each algorithm alone and in their hybrid combinations such as GA with PSO, ACO with PSO and ACO with GA are extensively investigated using the numerical examples with effects of noise added for realism. The GA+PSO hybrid algorithm was found to give the best performance in speed and accuracy compared to all others. The next best in performance was pure PSO followed by pure GA. ACO performed poorly in all the cases.


Introduction
Inverse problems often occur in many branches of engineering and mathematics where the values of some physical model parameter(s) must be deduced from observed data.System identification comes under the category of inverse problems.It is the process of determining the parameters of a system based on the observed input and

Classical Methods
The least square method is perhaps the simplest classical methods for structural identification.It estimates the unknown parameters of structural system by minimizing the sum of squared errors between the predicted outputs and the measured outputs.Many other classical methods such as maximum likelihood method, instrument variable method and extended Kalman filter method have been developed.Koh et al. (1991) used the Kalman filter method for identification of stiffness and damping coefficients from measurement of dynamic responses in the time domain.Detailed information about classical methods and non-linear identification are presented in the review paper by (Gaetan et al. 2006).Some classical optimization methods such as gradient methods (typically variations of Cauchy, Newton and Levenburg-Marquardt methods) need a good initial guess of the parameters to converge fast and accurately.Some studies have been conducted to supply good initial values using non-traditional methods such as GA which would sample the search domain efficiently.Thus certain classical methods such as Levensburg-Marquardt were hybridized with GA to identify system parameters with less computational time and improved accuracy (Kishore et al. 2007 andFriswell et al. 1998).

Evolutionary Algorithms and Behaviourally Inspired Algorithms
Evolutionary algorithms are robust optimization algorithms based on the heuristic concept of natural selection and genetic operations.Optimization algorithms find the minima or maxima of functions in a given function domain.Over the past few decades the area of genetic algorithm (GA) has been widely developed and applied for structural identification.Unlike classical methods, there is no requirement for calculation of gradients and second derivatives which frequently lead to numerical instability.GA applications in system identification such the identification of elastic properties of composite plates from dynamic test data is presented in (Jesiel et al. 1999;Chakraborthy and Mukhopadhyaya 2002).Recently efforts have been made to alter the architecture of GA and to incorporate local search algorithms to further improve its performance, see for example see (Koh et al. 2003 andPerry et al. 2004).
Behaviorally inspired optimization algorithms have been developed out of attempts to model the natural behavior of a flock of birds or a colony of ants.Social insects have diverse foraging systems that reflect the enormous range of ecological conditions in which they operate.Ant colony optimization (ACO) and Particle swarm optimization (PSO) techniques are inspired by real ant colonies and flock of birds respectively.

Ant Colony Optimization (ACO)
Ant colonies continually adjust foraging effort to changing conditions.Individuals use local information in foraging decisions and colonies can tune foraging effort to the location, quality and abundance of food.ACO is a general purpose Combinatorial Optimization (CO) technique.CO optimization is one of the youngest and most active areas of discrete mathematics.As the name suggest CO deals with finding optimal combination of available problem components.Hence, it is required that the problem is portioned into a finite set of components and CO algorithm attempts to find the optimal combination.ACO, a meta-heuristic which is based on the Ant System introduced in the early nineties by (Dorigo 1992 andDorigo et al. 1996).It has been used to solve a variety of combinatorial optimization problems, eg. the vehicle routing problem by (Bullnheimer et al. 1999) the traveling salesman problem by (Dorigo et al. 1997) and the industrial layout problems studied by (Hami et al. 2007;Corry andKozan 2004 andRajamani andAdil 1996).Abbospour et al. (2001) proposed the ACO-IM (Inverse Modeling) method for estimating soil hydraulic parameters.They compared results of ACO-IM with Levenberg-Marquardt (LM) method and finally complimented ACO and LM to obtain final solution which was better than pure ACO itself.Some extensions and variants of ACO are presented in the review paper by (Blum 2005).
Recently ACO has been applied to structural optimization and topology optimization.The goal of truss optimization is to optimally utilize the geometry and material of the proposed design elements to result in the lightest structure satisfying all the design, manufacturing constraints and other physical constraints.Since ACO is a discrete CO problem, several studies using it to optimize steel trusses for minimum weight using discrete cross sectional areas and other parameters were conducted, see (Camp et al. 2004;Camp et al. 2005 andSerra et al. 2006).These studies mapped the length of the tour of the ant to the weight of the truss and represented the volume of the element as virtual paths in ACO model for truss design.
A biomechanical application of PSO is presented by (Jaco et al. 2005) where ankle joint model parameters were identified and compared with gradient methods such as sequential quadratic programming, quasi-Newton algorithm and GA.The PSO algorithm was found insensitive to design variable scaling and gradient methods are highly sensitive.The parameters of Lorenz chaotic system are estimated using PSO in reference (Qie et al. 2007).It was found that PSO converges to exact value with high population size and more effective than GA.A 10-dof structural dynamic model was identified using frequency response function by GA and PSO in (Mouser and Dunn 2005), without using hybrids.PSO was found more likely than GA to produce better parametric models.A PSO-GA hybrid was used to successfully identify faults in the water supply system in Japanese cities by (Furuta and Yasui 2005).However only a single case is studied, and comparisons with pure GA and PSO are not shown for that case study, but only for benchmark equations.Recently PSO has been successfully applied in many research and application areas such as pattern recognition (Peng-Yeng et al. 2006), scheduling (Binghui et al. (2007)), layout optimization (Zhang et al. (2007), and the design optimization of Unmanned Aerial Vehicle (UAV) with respect to flight loads (Hu et al. 2004).

The ACO Algorithm
The inspiration for ACO is the foraging behavior of ants, particularly how they find the shortest paths between food sources and their nest.While searching for food, ants initially explore the area surrounding their nest randomly.When ants move they leave a chemical pheromone trail on the ground.When choosing their way, ants tend to favor paths marked by strong pheromone concentrations.As soon as an ant finds a food source, it evaluates the quantity and the quality of the food and carries some of it back to the nest.During the return trip, the quantity of pheromone that an ant leaves on the ground may depend on the quantity and quality of the food.The pheromone trails will guide other ants to the food source.
ACO model consists of graph G = (V, E) (Blum et al. (2005)) where V consists of two nodes, namely v s (nest) and v d (food source).E consists of two links e 1 and e 2 between v s and v d and corresponding length are l 1 and l 2 as shown in Fig. 1.
Path e 1 represents the short path between v s and v d , and e 2 represents the long path.Real ants deposit pheromone on the paths on which they move.Thus, the chemical pheromone trails are modeled as follows.
(1) where p i is the probability of an ant choosing the i th path, τ i is artificial pheromone value for each of the two links e i , i = 1, 2. Obviously, if τ 1 > τ 2 , the probability of choosing e 1 is higher, and vice versa.For returning from v d to v s , an ant uses the same path as it chose to reach v d and it changes the artificial pheromone value associated to the used edge.After all the ants have returned to the nest, the pheromone information is updated using Eq.(2) (2) wher p ε (0,1] is evaporation rate, Q is a constant, L k is length of the path traversed k th ant and na is number ants in the colony.The aim of pheromone update is to increase the pheromone value associated with good or promising paths.Pheromone evaporation is needed to avoid too rapid convergence of the algorithm. Implementing the discrete ACO algorithm for a continuous problem requires substantial modifications.The physical problems need to be represented in a graphical from as shown in Fig. 2 and Fig 3 .The elemental stiffness value is divided in to 'np' different values randomly with in the specified range.Each of these values can be represented as local path between the two levels and each of which leads to global path as shown in Fig. 2b and 3b.As mentioned earlier, ACO algorithm requires finite set of components, i.e. resolution of the local paths, so each stiffness value is divided in 'np' different paths (stiffness) in the given range.ACO algorithm has to find shortest path out of np nv (nv-number of parameters to be identified) available combinations.In real ant colony initially ants take random paths and deposit pheromone.As the iteration proceeds, all the ants in colony take the shortest path and there by increasing the intensity of pheromone trail on the shortest path.The main goal in structural identification is minimizing the objective function, which in turn decides the amount of pheromone deposition.Generally ACO minimizes the length of the L k in Eq. 2, traversed by k th ant.In structural identification, ACO minimizes the fitness value, ε Eq. 8, which is analogous to length of the path L k in Eq. 2 and the amount of pheromone deposited depends on the value of ε.

Particle Swarm Optimization (PSO)
Particle swarm optimization (PSO) is a population based continuous optimization technique developed by (Eberhart and Kennedy 1995), inspired by social behavior of bird flocking or fish schooling.PSO shares some similarities with evolutionary computation techniques such as Genetic Algorithms (GA).The system is initialized with a population of random solutions and searches for optima by updating generations.However, unlike GA, PSO has no evolution operators such as crossover and mutation.In PSO, the potential solutions, called particles, move through the problem space by following the current optimum particles.
The basic PSO algorithm consists of the velocity and position equation of the k th generation: (3) (4) i -particle number k -iteration number v -velocity of i th particle x -position of the ith particle/ present solution pi -historically best position/solution found by i th particle G -historically best position/solution found by all particles γ 1,2 -random number in the interval (0,1) applied to i th The frequently used PSO form includes an inertia term and acceleration constants, hence the velocity equation of PSO algorithm is modified as: (5) The inertia function is commonly either taken as constant or as a linearly decreasing function from 0.9 to 0.4.The acceleration constants are most commonly set to both equal to 2 as in (Shi and Eberhart 1998).The algorithm starts with initializing the i particles randomly with initial zero velocities.PSO then searches for optima by updating the positions of particle through several generations.At every generation, the location of each particle is updated based on two "best" values defined as follows.
The first is the historically best solution (fitness) achieved so far by all particles and stored as gbest.Another "best" value is the historically best value obtained so far by the i th particle in the population and is called pbest.The position of each particle is updated by a quantity which reflects the difference between gbest and pbest (equation ( 4) and ( 5)) and eventually all the particles tend to converge to the global best (gbest) position.The superiority of PSO over other comparable heuristic algorithms such as GA could be attributed to the explicit tracking and updating of gbest and pbest over the generations.

Genetic Algorithm (GA)
Genetic Algorithms are optimization algorithms based on the mechanism of natural selection and survival of the fittest.Over the past few decades GA has been widely developed and applied for structural identification.GA combines the explorative ability of large search domains as well as reasonable guided search to the global optima (Michaelwicz 1994).GA creates an initial random sample within the specified domain of variables, called 'population'.It then ranks them in the order of fitness and conducts crossover operations from among a pool of 'parents' through the Roulette wheel selection.Parents having higher fitness have a greater probability of being selected for crossovers and their offsprings contribute to the next generation.These offspring have marginally better fitness than the parents and over many generations they attain the global maxima.GA can be programmed in the Binary or Continuous versions.Here the Continuous (decimal number) version is used to avoid the computationally intensive conversion from binary to decimal and vice-versa.It been indicated in a few studies that continuous GA is superior to binary GA in computational performance (see (Michaelwicz 1994 andHaupt 1994) ).It may be noted that unlike PSO, GA does not explicitly keep track of the globally best solution (gbest) and particle best (pbest) solutions.

Numerical Examples and Discussion
For structural identification problems it is usually assumed that the mass of the structure is known and the identification is limited to structural stiffness and damping ( 1) properties.It is assumed that the structure is excited by known forces and that the responses of the structure in terms of accelerations, are recorded at some given points.
The dynamic equation of a structural system can be written as where M, C and K are the mass matrix, damping matrix and stiffness matrix, respectively, u is the displacement vector and F is the force vector.Proportional damping is adopted as follows ( 7) where α and β are two damping constants, which can be related to the modal damping ratio.In the examples which follow, experiments are simulated numerically using the known parameters of the structure ie. the forward analysis is carried out by numerically solving the dynamic equation using the parameter values (such as in Tables 1 and  4).The acceleration responses in the time domain are obtained at desired places and they are treated as experimentally obtained values.Noise may be added to them to be more realistic.The difference between the estimated acceleration response and measured acceleration response is used to compute the fitness value (objective function) as given by the Eq. ( 8) given below, which has to be minimized (see reference Koh et al. 1991). ( The parameters used in GA, PSO and ACO are given in the Table 2.They were applied to all the numerical examples.To compare GA and PSO they were given the same population size (200) and number of generations (50).These values are taken from Koh and Shankar (2003) which dealt with identification of structural systems similar to this paper.The crossover rate (40%) and mutation rate (5%) for GA, and PSO inertia function value (0.3) and acceleration constant (2) were selected from the best parameters recommended from standard literature on these algorithms (Michaelwicz 1994;Haupt 1994;Shi and Eberhart 1998).Comparable ACO parameters were however more difficult to establish since it is a discrete approach whereas PSO and GA are continuous with their resolution decided by the precision of the smallest decimal number -which is by default set to double precision in MATLAB.It is impractical to obtain a comparable precision in ACO by discrete division of the search interval (ie.number of paths).Thus the crucial ACO parameters such as number of ants (400), number of paths (100) and number of iterations (100) were chosen after several trials as a compromise, which would converge in a reasonably comparable time as GA, with mean errors which are acceptable for system identification purposes.Each algorithm has been run ten times separately to minimize the objective function, ε, for both impulse and random excitation with noise-free and noisy data and the final results shown is the average.

Example 1: 10-DOF Lumped Mass System
In the numerical study, a 10-DOF lumped mass system studied in (Koh et al. 2003) and shown in Fig. 2(a) is considered.The structural parameters are tabulated in the Table 1.Impulse and random excitation is applied at the 5 th level.The impulse was given as an initial velocity of 10 m/s to the first mass, the displacements and velocities of all other DOF set to zero.The impulse excitation was simulated by imposing equivalent initial velocity conditions obtained from impulse momentum relation using the method followed in (Hanagud 1985).The Gaussian random force was applied as with a maximum amplitude of 10 N. Accelerations at alternate levels, ie.levels 2, 4, 6, 8, 10 are assumed to be available for the purpose of structural identification.Referring to Eq. 8, here m=5, and n the time steps are 500 in the range from 0 to 2 seconds.Here the objective is to find out the unknown spring stiffnesses k 1 to k 10 .In all the problems hereafter the search limits for unknown parameters are taken as -50% (lower limit) to +100% (upper limit) of the exact value.In GA and PSO initial population/particles are generated within this specified range.In ACO a matrix of size nv × np is generated within the specified parameter range.ACO has to choose best combination of paths which minimize the objective function.
The results for this system are shown in the Fig. 4. GA and PSO have identified the parameters more accurately as shown in Fig. 4 (a) and (b).The accuracy of identified values by PSO is more than GA which can be noticed in the almost constant plateau in Fig. 4 (b) as compared to est u i j are, respectively the measured and estimated responses of i th measurement location at j th time step, m is the number of measurement location and n is the number of time steps.The "measured" responses are obtained from numerical simulation.The "estimated" responses are obtained from the mathematical model (Eq.6 and 7) where the optimization variables are t he now unknown stiffness properties and damping ratio.The values of the variables which minimize Eq. 8 give the required structural properties.To excite higher modes for better identification the input forces are impulse and random excitations (the latte r by means of Gaussian white noise).The dynamic equations are solved using Newmark's constant acceleration method.Two seconds of acceleration response is computed with constant time step of 0.004s in 500 steps.To investigate the effect of noisy data on identification, the I/O time signals are artificially contaminated by zero -mean Gaussian white noise.The noise level is defined as the ratio of standard deviation of the noise to the root -mean-square value of the uncontaminated time history; here it is taken as 10%.

GA (Fig 4(a)
) for stiffness levels 6-10.From Fig. 5 and 6 it is clear that GA and PSO performs better with presence of noise for both types of excitations (impulse and random).But ACO fails to identify the parameters with noisy data.Table 3 gives maximum and average errors (averaged over 10 runs) in identification of the stiffness parameters using impulse and random excitation.The maximum errors observed for the noisy case are 15.9% for ACO, 9.12% for GA and 4.17% for PSO (using random excitation).The mean errors are 12.34%, 3.47% and 2.88% respectively (using random excitation).The Table shows that using impulse excitation gives lesser errors than the random type.In such as case PSO produces a maximum error of only 3.02% with a mean error of 1.13%.Thus it is clear from the Fig. 5, Fig. 6 and Table 3, that the PSO gives superior results as compared to the other two algorithms for both types of excitations and with noisy data from the point of view of maximum and mean errors as well as computational time which is similar for both load cases.

Example 2: 11-Member Planar Truss
In this example an 11-member planar truss is considered for identification of the axial rigidity (EA) of all the 11 (ie.number of variables, nv = 11) members.Configuration and structural properties of the truss are given in Fig. 3(a) and Table 4, respectively.To implement ACO, the physical truss has to be represented as a graphical model with axial rigidity (stiffness) of each member discretized into 100 units ie.this would correspond to 100 paths per member.Fig. 3(b) shows schematically this division of each member into several paths.For clarity only 3 paths are shown per member.

Damping
Critical damping ratio for first 2 modes 5%   The truss is excited with impulse and random input at node 2 in Y-direction and acceleration responses (using the same excitation and times as the previous case) are assumed to be available at nodes 1, 2, 3, 4 and 5 in Ydirection.Acceleration responses are simulated using properties given in Table 4, and 10% noise added to the simulated responses.The truss is modeled using finite element method.The same set of algorithm parameters (Table 2) are used in this example also for minimizing the objective function ε (Eq.8).In all the algorithms finite element model is updated till the error between simulated acceleration responses and estimated by the algorithms is minimized.
Figure 7 shows the typical percentage error in identification of axial rigidity of the truss members for impulse excitation with and without 10% noise.From Fig. 7 it can be noticed that the error in identification of all the members is small in the case of PSO and GA.But ACO appears to oscillate about zero with significant variation giving maximum errors of about 20% -30%.The maximum and mean errors in identification by all algorithms for both excitations are tabulated in Table 5.It is seen that when noise is added to the responses the impulse excitation gives better identification than random excitation.The maximum and mean errors and computational time of PSO is much better than the others in the presence of noise.The computational time is roughly same for both load cases hence only a typical average value of time is shown to compare between the various algorithms.

Example 3: Simply Supported Beam
In this example a simply supported beam, Fig. 8, with linearly varying stiffness from support to the middle is considered for element wise stiffness identification (Modulus of Elasticity) in the time domain.Cross section of the beam is square with an area of 9x10 -4 m 2 and density of 7800 kg/m 3 .The beam is modeled with 2 noded, 11 Euler beam elements with two degrees of freedom per node, translation and rotation.
The stiffness variation is 200 GPa at the support to 50 GPa at the center of the beam.Impulse excitation is applied at node 3 in transverse direction and acceleration responses are assumed to be available at nodes 1, 4, 7 and 10 for identification.In identification, modulus of elasticity of all eleven elements are considered unknown and damping ratio of 5% is assumed to be known.Here also the finite element model is updated as mentioned in the previous example to minimize the objective function given in Eq. ( 8).Results are represented in Fig. 9 for impulse excitation only.Again it is noted that the performance of PSO is better than other algorithms.This is seen in the results where the mean error in identification is 5.8%, 3.1% and 15% for GA, PSO and ACO with noisy data, respectively.In all the examples thus presented it is found that ACO cannot identify the parameters with the same degree of accuracy as PSO and GA.So in order to combine the advantages of each of these algorithms, a hybrid approach is next investigated.

The Hybrid Approach
The results of hybrid approach are tabulated in the Tables 7 to 10. Table 7a shows the identification results with known damping for three cases A, B and C which represent the lumped mass, truss and beam models.The extensive data in Table 7a is summarized in Table 7b by averaging over the three cases.Table 7b data can be studied from the point of view of a) average time of solution : from this point of view it is seen that GA+PSO hybrid is clearly superior with only 39.7s, followed by pure PSO with 66.7s, then ACO+PSO hybrid with 110s and pure GA at 143s.The worst time performer is pure ACO which takes 219s.b) accuracy of solution: we examine the values for mean and maximum error in identification .Here also GA+PSO hybrid gives the smallest mean and maximum errors (0.72% and 2.47%) followed by pure PSO     (1.4% and 3.71%) and then pure GA (2.73% and 6.29%).
The ACO+PSO hybrid has marginally less accuracy than pure GA although speed wise it is better.The worst accuracy is shown by ACO.The discrete nature of ACO and its lack of resolution compared to GA and PSO are obvious handicaps.
Next, the same structural cases are studied, but damping is considered as an additional unknown.Tables 8, 9 and 10 deal with these cases.In each table the values of all unknown parameters are presented in detail, in addition to mean error, maximum error and time taken.The accuracy of identification of the unknown damping ratio by the algorithms is also an important parameter considered here.
In Table .8 it is seen that the GA+PSO hybrid is the fastest (ie.with least time of 23s) and next is PSO (34s) followed closely by ACO+PSO (78s) and GA (80s) and the worst is pure ACO(145s).The time trend is the same as in Table .7. Damping parameter has been identified most accurately by GA+PSO (-1.21% error) and pure PSO (-1.05%) followed by ACO+PSO and GA.Pure ACO gives the worst value here(16%).Regarding the overall accuracy (including damping as well as stiffness) it is seen that GA+PSO has identified fairly accurately with a small mean and maximum error of 1.47% and 4.67% followed by ACO+PSO (1.7%, 4.89%) , PSO (1.85%, 7.2%) and GA (2.01% and 5.41%).Thus there are some small deviations in the trend of accuracy when compared to Table 7.
Summary of observations: The general trend shown in Tables 7 to 10 is as follows.From the point of view of speed (minimum computation time) GA+PSO is the best, followed by pure PSO, then ACO+PSO and pure GA.In some cases GA is nearly comparable to ACO+PSO in speed.Pure ACO is the worst performer.From the point of Table 8.Percentage error in identified value of 10-DOF lumped mass system with unknown damping and noise view of accuracy, again GA+PSO is the best performer followed by pure PSO.Then follows GA and ACO+PSO -their relative superiority in accuracy varies from case to case.Even in these circumstances pure GA would still be preferred as it is simpler to program than the hybrid.The ACO+GA hybrid and pure ACO was not found promising in speed or accuracy.The observations from this comprehensive numerical study would be useful in future work, such as damage detection of structural members using these algorithms.

Conclusions
The paper studies the application of behaviorally inspired algorithms and their hybrids for structural parametric identification in the time domain.The stiffness and Colony Optimization (discrete) are used here.Three different hybridizations of GA, PSO and ACO have been investigated using three numerical models i.e., a lumped mass model, a truss and a non-uniform beam.Comparing the performance between pure GA and PSO it is found that the latter consistently outperforms the former in computational time and mean error, especially in the presence of noise.The computational implementation of PSO is also simpler than GA.Unlike GA, PSO keeps track of the global best and particle best solution from iteration to iteration which explains its superiority over GA and fast convergence.Regarding the performance of hybrid algorithms, the GA+PSO has shown its superiority over all the other algorithms in both speed and accuracy.The next best performer is pure PSO followed by pure GA.The performance of pure GA is however comparable to ACO+PSO hybrid in a few cases.However pure GA would be preferred as it is simpler to program than the hybrid.Pure ACO performed worst in all cases as expected because of its discrete nature.

Figure 4 .Figure 5 .
Figure 4. Identified values of the level stiffness v/s no. of runs with noise-free data for impulse excitation (a) GA (b) PSO (c) ACO

Figure 6 .
Figure 6.Identified values of the level stiffness v/s no. of runs with 10% noise data for random excitation (a) GA (b) PSO (c) ACO Figure 7. Percentage error in identification of axial rigidities of the members with impulse excivation (a) with out noise (b) with 10% noise

Figure 9 .
Figure 9. Identified stiffness along the length of the beam for impulse excitation (a) without noise (b) with 10% noise

Table 10 . Prcentage error in identified value of simply supported beam with unknown damping and noise damping
parameters are predicted from the time domain acceleration responses.Impulse and random excitations cases are studied.Particle Swarm Optimization (continuous), Genetic Algorithm (continuous) and classical Ant