Identification of structural models using a modified Artificial Bee Colony algorithm

Computers and Structures 116 (2013) 59–74

Contents lists available at SciVerse ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Identification of structural models using a modified Artificial Bee Colony algorithm Hao Sun a, Hilmi Lusß b, Raimondo Betti a,⇑ a b

Department of Civil Engineering and Engineering Mechanics, Columbia University, New York, NY 10027, USA Department of Civil Engineering, Bog˘aziçi University, 34342 Bebek, Istanbul, Turkey

a r t i c l e

i n f o

Article history: Received 15 June 2012 Accepted 17 October 2012 Available online 13 November 2012 Keywords: Parameter estimation Structural system identification Artificial Bee Colony algorithm Optimization

a b s t r a c t A modified version of the Artificial Bee Colony (ABC) algorithm is presented to identify structural systems. ABC is a heuristic algorithm with simple structure, ease of implementation and robustness. A nonlinear factor for convergence control is introduced in the algorithm to enhance the balance of global and local searches. To investigate the applicability of this proposed technique to system identification, three examples are studied under different conditions regarding data availability, noise pollution level, priori knowledge of parameters, etc. Simulation results show the proposed technique produces excellent parameter estimation, even with few measurements and high noise corruptions. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Structural system modeling and identification are extremely important areas in civil and mechanical engineering. Over the past few decades, they have been widely used for different purposes, with applications in structural health monitoring, nondestructive assessment, damage detection, active control, etc. From a computational point of view, system identification (SI) can be considered as an optimization or minimization process in which the objective is to determine a model of a system so that its predicted response to a given input is close enough to the measured response from the real system. Looking at the current literature, SI can be generally classified into parametric identification and nonparametric identification depending on the type of structural model used. If the SI is in terms of an assumed model defined by a set of physical parameters, such as mass and stiffness, then it is called parametric identification, while nonparametric identification is used to categorize methods that use purely mathematical representations of the system. Considerable efforts have been made in developing methods for parameter identification and state estimation of dynamical systems considering either a full or partial set of input as well as output measurements [1–10]. There are two main categories of approaches for SI, depending on whether SI is carried out in the frequency domain or in the time domain. Frequency domain methods seek to identify natural frequencies, mode shapes, damping ratios, and other modal quantities, using frequency domain based analyses of the time histories ⇑ Corresponding author. Address: Department of Civil Engineering and Engineering Mechanics, Columbia University, 610 S.W. Mudd Building, 500 West 120th Street, New York, NY 10027, USA. Tel.: +1 212 854 6388; fax: +1 212 854 6267. E-mail address: [email protected] (R. Betti). 0045-7949/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruc.2012.10.017

of the input and/or output response [1,6,11–16]. However, these tools are generally unreliable for damage assessment; in fact, frequencies are often more affected by environmental factors, such as temperature, than by the occurrence of structural damage, according to the experiments carried by Farrar and Doebling [17]. Moreover, mode shapes tend to be quite difficult to be accurately measured because of their susceptibility to noise. In addition, structural masses have to be assumed in order to draw conclusions on the stiffness from modal parameters like frequencies [18]. By contrast, time domain methods aim to estimate system physical parameters based on observed data sampled in the time domain; some of such methodologies can even track system parameters change during the duration of the records [19]. Various techniques have been developed in the time domain for both parametric and nonparametric identification, such as the least-square based estimation methods [3,4,20–23], the extended Kalman filter (EKF) [2,24–27], the unscented Kalman filter (UKF) [5,28–30], the particle filter (PF) [5,31], the H1 filter [32], the square-root unscented filter [33], and the sequential Monte Carlo method [5,34,35]. Lin et al., for example, used a variable trace approach [3] and a variable forgetting factor approach [36] to identify non-linear hysteretic structural systems on-line, while Yang and Lin [21] proposed a method for identification of parameter’s changes based on least squares estimation and adaptive tracking technique. The EKF with a weighted global iteration procedure was applied to identify structural systems by Hoshiya and Saito [24] and Jeen-Shang and Yigong [25] using earthquake records. Yang et al. [27] developed an adaptive tracking technique based on the EKF to identify the structural parameters. Such a technique is able to track the changes of system parameters on-line. The UKF, a recently developed filtering technique with many advantages

60

H. Sun et al. / Computers and Structures 116 (2013) 59–74

over the classical EKF, has been also employed in SI of both linear and nonlinear structures as evidenced by the works of Wu and Smyth [28], Mariani and Ghisi [29], and Chatzi and Smyth [5]. The study by Xie and Feng [30] shows that the iterative UKF produces better state estimation and parameter identification than the UKF for nonlinear structural SI. Sato and Qi [32] used the adaptive H1 filter for identifying multiple-floor nonlinear hysteretic structures. The implementation of sequential Monte Carlo method has proven to be successful in the identification of nonlinear systems with non-Gaussian uncertainties, but it often requires a very large number of sample points and is therefore computationally expensive [28]. However, most of the above mentioned methods require a good initial guess of the parameters and a proper function gradient. In addition, difficulties arise using these methods in the identification of large systems when few measurement data is available. Hence, methodologies based on heuristic stochastic algorithms to solve the optimization problem in structural SI have been employed in recent years: in particular, genetic algorithms (GAs), particle swarm optimization (PSO), ant colony optimization (ACO), artificial neural networks (ANN), evolutionary strategy (ES), and differential evolution (DE) algorithms have gained increasing attention. Koh et al. [37] proposed a hybrid computational strategy combining GA with a compatible local search operator for large structure parameter identification, while Wang [19] developed a hybrid GA integrated with the Gauss–Newton method, to identify both linear and nonlinear structural systems. These hybrid GA approaches have been shown to perform better than the classic GA in structural system identification. Kwok et al. [38] constructed a novel hysteretic model for magnetorheological fluid dampers and used PSO to identify the model parameters, while Charalampakis and Dimou [39] applied an enhanced PSO to the identification of Bouc–Wen hysteretic systems. The works done by Lefik and Schrefler [40] and Karimi et al. [41] show that the ANN procedure performs well in parameter identification of structural systems under dynamic loading, while Franco et al. [18] presented a parameter estimation technique based on ES for parametric identification. Tang et al. [42] adopted DE to estimate structural system parameters with and without noise pollution. All these methods have generally achieved satisfactory results in solving the multimodal optimization problem in SI. Apart from the above mentioned heuristic algorithms, a novel swarm intelligence algorithm called the Artificial Bee Colony (ABC) algorithm was introduced by Karaboga [43] in 2005 for solving complex numerical optimization problems. This algorithm is motivated by the intelligent behavior of honey bees when seeking a high quality food source. ABC is a population based stochastic algorithm with implementation simplicity since the only common control parameters are the colony size and maximum cycle number. It has the advantage of simple structure (as simple as PSO and DE), ease of use, and high stability. Karaboga and Basturk [44] used it to optimize multivariable functions and compared their results with those produced by GA, PSO, DE, and particle swarm inspired evolutionary algorithm (PS-EA). Other studies have also compared the performance of ABC with other existing heuristic algorithms [45,46], and the results reported in the literature suggest that ABC can be more effective than other methods in some optimization problems [47–49]. Recently, several researchers have paid attention to the application of ABC in civil engineering. Kang et al. [50] proposed a hybrid algorithm combining Nelder– Mead simplex method and ABC to solve inverse problems in concrete dam structures. Sonmez [51,52] applied ABC with an adaptive penalty function approach to minimize the weight of truss structures with both continuous and discrete variables. Omkar et al. [53] developed a generic model for multi-objective design optimization of laminated composite components based on Vector

Evaluated ABC. Since this algorithm has been proven to be successful in dealing with large optimization problems, it seems natural to extend this algorithm to structural parameter identification, considering the problems associated with the limitations that occur in real life applications such as incomplete set of measurements and noisy data. The paper is organized as follows. Section 2 describes the formulation of structural SI as an optimization problem. Section 3 presents a modified version of ABC algorithm which is suitable for structural parameter estimation. Section 4 discusses the numerical simulations and identification results obtained from two linear systems (5-DOF and 20-DOF) and a nonlinear system (2-DOF) identified with full and partial measurements under noise free and noise polluted cases. In section 5, we present the observations and the conclusions of this work. 2. Optimization formulation of system identification System identification can be formulated as an optimization problem in which an objective function, for example the error between the actual measured structural response and the estimated response of a model, is defined and the parameters of such models determined so to optimize (usually to maximize or to minimize) the objective function. Let’s consider that u(tk) and y(tk) are two vectors containing the input and the corresponding output, respectively, of a general system at a generic instant tk (k = 1, 2, . . . , ntime). Let us now denote ^ðtk Þ the estimated system’s response obtained from a syswith y tem’s model using the same input u(tk), expressed as:

^ðtk Þ ¼ f ðuðtk Þ; ^hÞ y

ð1Þ

where ^ h ¼ f^ h1 ; ^ h2 ; . . . ; ^ hn g 2 Rn represents a vector containing the estimated values of the model’s unknown parameters h ¼ fh1 ; h2 ; . . . ; hn g 2 Rn : in a mechanical and/or structural system, for example, the elements hj (j = 1, 2, . . . , nelement) could be the mass, damping and stiffness parameters. The objective of any identification procedure is to find the best estimates ^ h of the structural parameters h so to minimize the error between the measured response y(tk) and the predicted (or esti^ ðt k Þ. This process can be visualized in Fig. 1. In mated) response y identification problems that rely on the dynamic measurements of the response, the objective function to be minimized through the optimization process can be formulated as a sum of the mean square error (MSE) between the measured and predicted responses, such as: n

gð^hÞ ¼

nrec X time X 1 ^i ðt k Þ  yi ðtk ÞÞ2 ðy nrec ntime i¼1 k¼1

ð2Þ

where nrec is the number of recorded time histories of the response. The number of response time histories, nrec, depends on the number of sensors available in the test and can range from 1 (only one sensor available) to the total number of degrees-of-freedom of the system (full set of sensors). The overall identification problem can be summarized as:

Find ^h ¼ f^h1 ; ^h2 ; . . . ^hn g 2 C such that gð^hÞ ! min

ð3Þ

where C is the feasible n-dimensional parameter search space:

C ¼ f^h 2 Rn jhmin 6 ^hj 6 hmax ; j j

8j ¼ 1; 2; . . . ; ng

ð4Þ

where n is the number of parameters to be identified, hmin and hmax j j are the lower and upper bounds of the jth parameter. Eqs. (2)–(4) show that the problem of identification can be treated as a linearly constrained nonlinear optimization problem [18]. Hence, because of the irregularity of the multi-dimensional surface

H. Sun et al. / Computers and Structures 116 (2013) 59–74

61

Fig. 1. The process of system identification as an optimization problem.

gð^ hÞ, an efficient identification algorithm must be capable of handling problems related to multiple solutions (more than one minimum) and to local vs. global minima (many local optima with high complexities may exist in the search space). While some of these problems, such as multiple solutions, can be resolved independently of the algorithm used, for example by choosing proper sensor locations that guarantees uniqueness of the solution, others, like local vs. global minima, have to be addressed directly by the chosen algorithm. 3. Artificial Bee Colony (ABC) algorithm The Artificial Bee Colony (ABC) algorithm is inspired by the natural behavior of honey bees when foraging (seeking high quality food sources) [43–46]. In applying the ABC algorithm to system identification problems, the starting point is to accept the equivalence between ‘‘food source location’’ and ‘‘structural model’’. This can be explained by considering that a food source location is similar to a structural model in which it is characterized by a set of parameters, such as the distance from the hive, the orientation with respect to the sun, etc.; in an optimization scheme, a location represents a possible solution. Finding the ‘‘best’’ possible food source is equivalent to finding the set of parameters that corresponds to the optimized solution. The initial step of the ABC algorithm is to select an initial population of possible food source locations or, in the context of structural system identification, of possible structural models: this initial set, of size np, can be randomly generated, keeping the parameters within some specified ranges. Of all these initial np models (or source locations), those nf models that yield the smallest objective function values in Eq. (3) are selected as the initial population. nf (number of food sources) indicates the number of possible models that will be considered for upgrading in the following iterations. The generic ith possible solution (either intended as food source location or structural model) will be indicated by a vector ^ hi ¼ f^ hi1 ; ^ hi2 ; . . . ; ^ hin g, of dimension n, where n is the number of parameters to be optimized. Starting from the initial set of solutions, a search for ‘‘improved solutions’’ or ‘‘better fit solutions’’ in the neighborhood of each existing one is conducted. This is equivalent to the task, carried out by ‘‘employed bees’’ or ‘‘employed foragers’’, of trying to improve their food sources by a neighborhood search. Hence, for each original solution, a new candidate solution is randomly generated and the objective function value of such a new solution is computed by Eq. (2). If the fitness value of the candidate is better than that of the previous solution, then the previous solution will be updated to the new candidate one; otherwise, the previous solution is kept. At the generic iterative step, the new candidate solution is represented by a vector ^ hup i generated using the following updating procedure:

^hup ¼ ^hij þ u ð^hij  ^hkj Þ ij ij

ð5Þ

where ^ hup ij indicates the jth component (j = 1, 2, . . . , n) of the new ith candidate solution (i = 1, 2, . . . , nf). ^ hij represents the jth component of the current ith solution while ^ hkj is the jth component of a randomly selected kth solution (k = 1, 2, . . . , nf), with k – i. The parameter uij is a randomly distributed number in the range ½k; k (usually k ¼ 1 in a standard ABC algorithm) and is responsible for the nature of the search: if k is a small number, then the algorithm will adapt more to a local search; otherwise, the algorithm will function for more global searches. It would then be preferable to have a parameter uij that changes value during the search process, from large values at the beginning of the search (where it is important to explore the entire horizon (space)) to small values towards the end of the search (where it is important to refine the convergence to the final optimal solution). In order to balance the search preference and to control the convergence rate as the iteration progresses, a nonlinear factor viter, defined as:

  iter  dm  miter 

viter ¼ 1  

ð6Þ

is proposed to modify the factor uij, where iter denotes the current iteration number, miter is the maximum number of iterations, d is an integer in the range between [0, miter/2], and m is the power exponent. By changing the values of the parameters d and m, it is possible to control the nature of the search process. Fig. 2 depicts the behavior of the nonlinear factor viter as a function of the iteration number, iter, for different values of m and d (miter = 500). Both plots show a nonlinear, decreasing behavior of viter as the iteration number increases: at the beginning (small iteration number), the value of viter approaches 1, favoring the search over the entire search space. However, as the number of iterations increases, viter approaches, more or less rapidly, to zero, inducing small variation to the current solution. The idea behind viter is to free the method, at its initial steps, to search over the entire search space and to allow for faster convergence during the local search phase which take place towards the end of the iteration process. To avoid stagnation around local minima, a concept equivalent to the figure of a ‘‘scout bee’’ in a hive will later be described. This newly defined viter factor is now used to adapt the variability of the new candidate solution that can be expressed as:

^hup ¼ ^hij þ v u ð^hij  ^hkj Þ iter ij ij

ð7Þ

In order to select which ith food source to retain (either the current one, ^ hi , or the new candidate, ^ hup i ), a ‘‘greedy’’ selection based on roulette wheel selection strategy is carried out by comparing the fitness values corresponding to ^ hi and ^ hup i . The fitness function can be expressed as:

fitð^hi Þ ¼ ½1 þ gð^hi Þ1

ð8Þ

where gð^ hi Þ and fitð^ hi Þ are the objective function and the fitness values, respectively, for the solution ^ hi .

62

H. Sun et al. / Computers and Structures 116 (2013) 59–74

Fig. 2. The relationship between iterations and the nonlinear factor: (a) the effect of m on viter with d = 0 and (b) the effect of d on viter with m = 4.

Once all the possible candidate solutions have been identified and the values of the corresponding fitness functions have been determined, then it is possible to associate a probability of selection to each candidate. Here, a fitness-proportional selection process is used so that the probability of selection associated with the ith candidate is:

fitð^hi Þ pi ¼ Pnf ^ k fitðhk Þ

ð9Þ

In this way, a food source with high value of the fitness function will have a higher probability of being selected. At this point, to maintain a certain randomness in the process, the new set of nf solutions ^ hi ¼ f^ hi1 ; ^ hi2 ; . . . ; ^ hin gT to be used in the next iteration is selected at random by generating a uniformly distributed random number ri(i = 1, 2, . . . , nf), between 0 and 1, and comparing it with the corresponding probability pi. If the probability associated with the ith solution, pi, is greater than the associated random number ri, then the ith updated solution is selected and passed to the next iteration. If instead, pi is smaller than ri, then a new ith solution is generated following Eq. (9) and its fitness compared to the value of the current one. The solution with the largest fitness value will be chosen as representative of the ith solution for the next iteration. The new set of nf possible solution ^ hi will represent the starting point of the next iteration. The solution (i.e., food source location or system’s model) with the largest fitness value among all the nf locations is recorded and compared with that of the previous iteration. It might happen that one of the possible solutions, let’s say the lth solution, is not improved after a given number, lc, of iterations while the other solutions keep improving their fitness: this can happen if the lth solution represents a local minimum of the fitness function. At this point, it will be ‘‘abandoned’’ and a new solution, ^ hsc i , will be obtained, with the jth element expressed as:

  ^hsec ¼ ^hmin þ l hmax  hmin j j ij j ij

ð10Þ

where lij is a uniformly distributed random number in the [0, 1] range while hmax and hmin are the upper and lower bound, respecj j tively, of the jth parameter. The superscript sc is chosen to indicate the candidates ‘‘scout bee’’ nature, which is a bee with the task of searching for a new food source. The use of such a ‘‘scout bee’’ solution is used to enhance the global search capability of the algorithm, providing a way to exit stagnation points (i.e., local minima). If the newly found ith solution is ‘‘better’’ than the abandoned one, the search will continue using the new location as the new ith one. Otherwise, a new solution will be generated and tested again until a new better solution is found. To summarize, the basic steps of this modified ABC algorithm are shown in Table 1, where the phases have been named after the different types of bees: employed, onlooker and scout. It is noteworthy that the difference between the standard and the modified versions is dependent of the use of viter in ABC. The modified ABC algorithm can be reduced to the standard version by setting viter = 1. This modified ABC algorithm promises to be a robust optimization algorithm for both global and local searches. In the identification of a structural system, the unknown model’s parameters, for example the mass, stiffness and damping parameters, and/or coefficients describing nonlinear behavior, play the role of parameters ^ hj that need to be estimated and the parametric space is the solution search space in the optimization process. A schematic representation of the proposed algorithm in the context of system identification is presented in Fig. 3. 4. Numerical examples In order to investigate the effective applicability of the proposed approach to structural SI, two linear systems (5-DOF and 20-DOF)

Table 1 The basic steps of the modified ABC algorithm. Initialization: initialize the set of possible solutions (following a uniformly random distribution in a specified search space). REPEAT  EMPLOYED PHASE search: for each possible solution, create a new candidate solution by Eq. (7), evaluate its fitness by Eq. (8), and select, among the two solutions, the one with the best fitness as the new possible solution; calculate the probability of each solution using Eq. (9);  ONLOOKER PHASE search: select possible solutions from those obtained in the employed phase by using the roulette wheel selection strategy; create new candidates by Eq. (7); evaluate each of them by Eq. (8), and select the better ones as new possible solutions; if needed (The limit cycle lc is reached):  SCOUT PHASE search: investigate a new solution by Eq. (10), evaluate it by Eq. (8) and continue until a better solution is found;  Record the best solution, the values of the fitness and of the objective function; UNTIL (The maximum cycle miter is reached)

H. Sun et al. / Computers and Structures 116 (2013) 59–74

63

Fig. 3. The layout of system identification using ABC.

and a nonlinear system (2-DOF) are discussed in this work. The three systems are 2-D shear frame type structures with lumped masses at each floor. Comparisons of the proposed technique for parametric identification with other methodologies suggested in literature [18,39,42,54–56] are carried out in this part. To verify the effectiveness of ABC for SI with limited output data measurements and to learn the influence of data availability on the performance of this algorithm, different sets of time histories of the structural response have been considered: both the case when the time histories are measured at every DOF (full instrumentation setup) and the case when the time histories are available at some DOFs (partial instrumentation setup) have been analyzed. Since accelerometers are commonly used in dynamic field testing, the structural response is presented here in the form of time histories of the structural accelerations. To provide a fairly reasonable search space for the parameters, the lower bounds of the search space are taken to be half and the upper bounds twice the actual values of the parameters so that the search space can be defined as h/2 6 C 6 2h, where h is the parameter vector representing the true values. To analyze the effect of measurement noise on the identification of the parameters and so on the effectiveness of the modified ABC algorithm, various noise pollution levels in both input and output data have been taken into consideration in the identification analysis: both input and output (I/O) signals have been polluted with Gaussian zero-mean white noise sequences. The root mean square (RMS) values of the noise sequences are defined as some certain percentages (i.e., 0% (no noise), 5%, and 10%) of the original time history signals. It is noteworthy that, in addition to the white noise cases, the cases of colored input and of colored measurement noise were also considered. However, the identification results for the colored input and colored noise are almost identical or even better than those obtained for the corresponding white noise cases and so have not been included in the paper for space limitation. In some previous studies on system identification, the assumption of an a-priori known mass matrix has been used quite frequently. This assumption is justified by the assertion that a sufficiently accurate estimation of the mass can be obtained from the structural drawing. Here, to test the iden-

tification capabilities of the proposed algorithm, both the case of a priori known mass distribution and the case of the unknown mass distribution have been considered. In applying ABC to SI, the most time-consuming operation is represented by solving the second order ODEs during the fitness evaluations, whereas the time taken by the parameter updates is comparatively negligible. In the first two examples, the systems’ responses are simulated using the implicit Newmark-b integration method, while in the third example, obtained by solving (integrating) the ODEs using the fourth to fifth order embedded Runge–Kutta integration method with adaptive step-size after transforming the ODEs into state space form. Moreover, since the ABC is a stochastic algorithm, IR (i.e., 20 in this study) independent runs have been conducted so as to obtain statistical insights into the process. In these examples, the number of food sources (nf) has been selected equal to half of the number of population size, namely np/ 2. The numerical analyses are programmed in MATLABÒ on a standard Intel (R) Core (TM) 2 Q9550 2.83 GHz personal PC with 2G RAM. 4.1. A 5-DOF linear system To study the performance of the modified ABC algorithm, first let us consider a 5-DOF linear structural system, modeled as a shear type frame structure (Fig. 4). The equation of motion for such a system can be written as:

€ðtÞ þ ½CxðtÞ _ ½Mx þ ½KxðtÞ ¼ FðtÞ

ð11Þ

€ are displacement, velocity, and acceleration vecwhere x, x_ and x tors; F(t) is the external force input vector; [M], [C] and [K] are the mass, damping, and stiffness matrices. For shear type models, the mass matrix is diagonal and the stiffness matrix is tri-diagonal; here we assume the damping matrix also has a tri-diagonal structure analogous to the stiffness matrix. The properties of the system in this example are shown in Table 2. For the known mass identification case, the parameters to be identified in this example can be described by a set of ten variables:

64

H. Sun et al. / Computers and Structures 116 (2013) 59–74 Table 3 Parameter values of the modified ABC algorithm for system identification.

Fig. 4. A n-DOF shear type linear structural system.

Table 2 The structural properties of a 5-DOF linear system. Parameters

Value Mode Frequencies (rad/s) Damping ratios (%)

Mass (m1, m2, . . . , m5)

0.1

1 2

10.2625 29.9560

1.9736 5.7608

Damping (c1, c2, . . . , c5) 0.5

3

47.2227

9.0813

Stiffness (k1, k2, . . . , k5) 130

4 5

60.6637 69.1900

11.6661 13.3058

h1 ¼ fc1 ; . . . ; c5 ; k1 ; . . . ; k5 g

ð12Þ

while for the unknown mass case, the system can be fully described by a set of 15 variables:

h2 ¼ fm1 ; . . . ; m5 ; c1 ; . . . ; c5 ; k1 ; . . . ; k5 g

ð13Þ

In this numerical example, the structure is assumed to be excited by an input force vector {0, 0, f3(t), 0, 0}T, where f3(t) is supposed to be a known white noise sequence of 5 s (following Gaussian distribution with a RMS scaled to 1000) acting on the third DOF. The sampling time is DT = 0.005 s. In the full output scenario, all the floor accelerations are measured, whereas in the partial output scenario only the accelerations at the 2nd, 3rd and 4th DOFs are available. The choice of such floors was dictated by the requirements that guarantee uniqueness of the identified solution [18]. The values of the algorithm parameters used in the modified ABC are illustrated in Table 3. The settings of m and d follow the assumption that viter starts to nonlinearly decrease not early than 1/3 of miter iterations and reaches the final value not smaller than 0.5. For comparison purposes, the standard ABC algorithm was tested, in which the algorithmic parameters are the same as those used in the modified ABC algorithm (apart from the nonlinear factor viter). A classic rank-based genetic algorithm (GA) was also implemented here by using the MATLABÒ GA Toolbox™. In this GA, the population size is chosen to be equal to np, the maximum generation equals to miter, and the number of independent runs the same as IR. The rate of crossover was selected equal to 0.8 while a Gaussian function with zero mean value was adopted for mutation strategy with a scale factor of 0.5 and shrink factor of 0.75 [57]. For the unknown mass case with complete instrumentation, an algorithm presented in [7,9,10], based on the identification of a first-order model and on its convertion to second order, is used for comparison purpose.

Mass cases

np

miter

lc

m

d

IR

Known Unknown

50 80

500 1000

300 300

4 5

100 150

20 20

Simulation results for the known mass case are summarized in Tables 4 and 5 corresponding to full output scenario and partial output scenario, respectively. It can be observed that the modified ABC algorithm is able to estimate the correct values of the structural parameters in the noise free case even with partial measurements, whereas the results obtained by the standard ABC algorithm and GA have slight errors. Furthermore, it is seen from Table 4 that the modified ABC algorithm yields ‘‘finer’’ results (with errors from 0% to 0.6%) than those obtained via the standard ABC algorithm (with errors from 0.1% to 1.5%) for both noise free and 10% noise cases. However, both versions of ABC algorithm seem relatively insensitive to noise corruption. Similarly, in the partial output scenario, both the modified ABC algorithm and GA perform well and yield good estimation results with small errors, ranging from 0.1% to 1.9% for ABC, and from 2.9% to 6.2% for GA, yet with the modified ABC algorithm performing slightly better than GA. Even for the 10% noise polluted case with partial measurements, the results obtained by the modified ABC algorithm show that estimation errors pertaining to damping and stiffness parameters remain very small with maximum values of 1.9% and 0.3%, respectively. Fig. 5 presents a noise free and noise polluted comparison of identified damping and stiffness parameters for the known mass case with partial measurements: the parameter estimates are observed to converge to the true values after 100 iterations, which demonstrates the accuracy and the high efficiency of the proposed approach. A more difficult case for the 5-DOF system is when the mass matrix of the system is also unknown and needs to be identified. Tables 6 and 7 present the results of the identification for both the full output scenario and the partial output scenarios, respectively. In the full output scenario, both the modified ABC algorithm and the algorithm proposed in [7] provide good identification results, with a maximum error of 1.8%, for the noise free case. However, in the 10% RMS noise corruption in the I/O signals, the modified ABC algorithm estimated errors of mass, stiffness, and damping parameters range from 0.1% to 1.2%, from 0.3% to 2.8%, and 0.1% to 1.0%, respectively, still providing excellent identification even for large noise levels. On the contrary, while the algorithm in [7,9,10] is able to estimate mass and stiffness well, it yields very large errors in damping coefficients identification (up to 49.6%). In the partial output scenario, the methodology in [7,9,10] cannot be applied and so the results of the proposed ABC algorithm are compared with those obtained through the GA previously considered. Again, it was found that the results by the modified ABC algorithm are much better than those obtained by GA for both the noise free and 10% RMS noise contamination cases despite noise corruption and lack of measurements. The maximum errors obtained by GA jumps from 6.3% to 15.4% as noise level increases from noise free to 10% RMS noise, whereas those obtained by ABC remain basically the same (from 2.7% to 3.9%). Similar to Fig. 5 for the above known mass case, Fig. 6 shows a noise free and noise polluted comparison of the identified parameters for the unknown mass case with partial measurements. In both cases, convergences are quite fast (300 iterations). It can be seen that the maximum estimate error often exists in damping parameters and signal noise corruption has usually larger impact on damping parameter estimation. For both the known and the unknown mass cases, the modified ABC algorithm also gives

65

H. Sun et al. / Computers and Structures 116 (2013) 59–74 Table 4 Results for 5-DOF known mass system in the full output scenario. Parameters

True value

0% Noise

10% Noise

The standard ABC

c1 c2 c3 c4 c5 k1 k2 k3 k4 k5

0.5 0.5 0.5 0.5 0.5 130 130 130 130 130

Max error (%)

The modified ABC

Mean

Std.

Mean

0.4996 (0.1) 0.5000 (0.0) 0.5003 (0.1) 0.5000 (0.0) 0.5000 (0.0) 130.03 (0.02) 129.96 (0.03) 130.00 (0.0) 130.00 (0.0) 130.00 (0.0)

0.0023 0.0010 0.0009 0.0010 0.0011 0.0905 0.0914 0.0342 0.0623 0.0701

0.5000 0.5000 0.5000 0.5000 0.5000 130.00 130.00 130.00 130.00 130.00

0.1

(0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0)

The standard ABC

The modified ABC

Std.

Mean

Std.

Mean

0.0012 0.0009 0.0014 0.0010 0.0008 0.0458 0.0708 0.0322 0.0311 0.0395

0.4927 (1.5) 0.4996 (0.1) 0.5020 (0.4) 0.4994 (0.1) 0.4985 (0.3) 130.19 (0.1) 129.72 (0.2) 129.95 (0.04) 130.22 (0.2) 129.95 (0.04)

0.0260 0.0167 0.0100 0.0056 0.0140 0.6450 1.0640 0.6373 0.5367 0.6118

0.4989 0.4986 0.4994 0.4970 0.5011 129.81 129.95 130.03 130.25 129.92

0.0

1.5

Std. (0.2) (0.3) (0.1) (0.6) (0.2) (0.1) (0.04) (0.02) (0.2) (0.1)

0.0222 0.0155 0.0100 0.0086 0.0092 0.9932 0.7316 0.7476 0.4750 0.7009

0.6

Note: The values in parentheses are absolute relative errors in %. Std. represents standard deviation.

Table 5 Results for 5-DOF known mass system in the partial output scenario. Parameters

True value

0% Noise

10% Noise

GA

The modified ABC

Mean c1 c2 c3 c4 c5 k1 k2 k3 k4 k5 Max error (%)

0.5 0.5 0.5 0.5 0.5 130 130 130 130 130

0.5149 0.4885 0.5076 0.4988 0.4907 132.41 133.44 129.73 131.14 132.03 2.9

(2.9) (2.3) (1.5) (0.2) (1.9) (1.9) (2.6) (0.2) (0.9) (1.6)

Std.

Mean

0.0543 0.0487 0.0201 0.0074 0.0195 1.9753 2.0017 0.6358 1.0244 1.8628

0.5000 0.4993 0.5000 0.5000 0.5000 130.00 130.00 130.00 130.00 130.00

(0.0) (0.1) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0)

GA Std.

Mean

0.0020 0.0033 0.0006 0.0004 0.0006 0.0333 0.0415 0.0632 0.0490 0.0403

0.5238 0.4821 0.5261 0.4690 0.5281 132.59 132.03 133.22 130.84 132.57

0.1

6.2

The modified ABC

(4.8) (3.6) (5.2) (6.2) (5.6) (2.0) (1.6) (2.5) (0.6) (2.0)

Std.

Mean

Std.

0.0746 0.0397 0.0575 0.0801 0.0749 2.4832 2.0798 2.9945 0.9761 2.6185

0.4903 (1.9) 0.5030 (0.6) 0.4991 (0.2) 0.4984 (0.3) 0.5035 (0.7) 129.64 (0.3) 130.12 (0.1) 130.11 (0.1) 130.09 (0.1) 130.23 (0.2)

0.0332 0.0226 0.0126 0.0115 0.0237 1.3424 1.9780 0.9188 0.6810 1.0021

1.9

Note: The values in parentheses are absolute relative errors in %. Std. represents standard deviation.

Fig. 5. A comparison simulation results for the known mass case with partial measurements: (a) noise free and (b) 10% noise corruption in I/O signals (to be reproduced in color on the Web and in black-and-white in print).

much smaller deviation values than those by GA (Table 7), showing a good level of accuracy despite measurement noise and lack of data.

To test the robustness of the proposed technique in each independent run, the probability density functions (PDFs) of the normalized parameters have been obtained and are presented in

66

H. Sun et al. / Computers and Structures 116 (2013) 59–74

Table 6 Results for 5-DOF unknown mass system in the full output scenario. Parameters

m1 m2 m3 m4 m5 c1 c2 c3 c4 c5 k1 k2 k3 k4 k5

True value

0% Noise

0.1 0.1 0.1 0.1 0.1 0.5 0.5 0.5 0.5 0.5 130 130 130 130 130

Refs. [7,9,10]

The modified ABC

Refs. [7,9,10]

The modified ABC

Mean

Mean

Std.

Mean

Mean

Std.

0.1009 (0.9) 0.1002 (0.2) 0.1001 (0.1) 0.1009 (0.9) 0.1011(1.1) 0.5053 (1.1) 0.5092 (1.8) 0.4931 (1.4) 0.5000 (0.0) 0.5089 (1.8) 131.126 (0.9) 132.242 (1.7) 130.183 (0.1) 130.500 (0.4) 132.071 (1.6)

0.0036 0.0015 0.0002 0.0021 0.0035 0.0291 0.0217 0.0159 0.0142 0.0156 1.2925 3.0893 1.4027 1.1719 2.0012

0.0988 0.0974 0.0999 0.0973 0.0975 0.6126 0.7481 0.5014 0.4904 0.5745 129.94 120.43 129.25 130.48 123.86

0.1012 (1.2) 0.1001 (0.1) 0.1000 (0.0) 0.1012 (1.2) 0.1009 (0.9) 0.4910 (1.8) 0.5138 (2.8) 0.4926 (1.5) 0.5015 (0.3) 0.5062 (1.2) 131.189 (0.9) 131.185 (0.9) 129.901 (0.1) 130.823 (0.6) 131.333 (1.0)

0.0035 0.0014 0.0005 0.0020 0.0038 0.0241 0.0282 0.0205 0.0101 0.0174 6.6079 3.1502 1.4196 1.5459 4.0275

0.1000 0.1000 0.1000 0.1000 0.1000 0.5000 0.5000 0.5000 0.5000 0.5000 130.00 130.00 130.00 130.00 130.00

Max error (%)

10% Noise

(0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0) (0.0)

0.0

1.8

(1.2) (2.6) (0.1) (2.7) (2.5) (22.5) (49.6) (0.3) (1.9) (14.9) (0.1) (7.4) (0.6) (0.4) (4.7)

49.6

2.8

Note: The values in parentheses are absolute relative errors in %. Std. represents standard deviation.

Table 7 Results for 5-DOF unknown mass system in the partial output scenario. Parameters

True value

0% Noise

10% Noise

GA

m1 m2 m3 m4 m5 c1 c2 c3 c4 c5 k1 k2 k3 k4 k5

0.1 0.1 0.1 0.1 0.1 0.5 0.5 0.5 0.5 0.5 130 130 130 130 130

Max error (%)

The modified ABC

GA

Mean

Std.

Mean

Std.

Mean

0.1034 (3.4) 0.1032 (3.2) 0.1000 (0.0) 0.1026 (2.6) 0.1026 (2.6) 0.5309 (6.2) 0.4685 (6.3) 0.5288 (5.8) 0.4988 (0.2) 0.4907 (1.9) 137.41 (5.7) 133.44 (2.6) 129.73 (0.2) 131.14 (0.9) 134.02 (3.1)

0.0081 0.0067 0.0003 0.0072 0.0095 0.0910 0.1024 0.0782 0.0093 0.055 7.8235 3.0847 0.7119 2.3970 7.0847

0.1014 (1.4) 0.1010 (1.0) 0.1001 (0.1) 0.1014 (1.4) 0.0997 (0.3) 0.5061 (1.2) 0.5101 (2.0) 0.4979 (0.4) 0.5111 (2.2) 0.4867 (2.7) 131.23 (0.9) 131.79 (1.4) 130.75 (0.6) 130.48 (0.4) 130.73 (0.6)

0.0039 0.0024 0.0001 0.0024 0.0045 0.0399 0.0484 0.0249 0.0176 0.0222 3.3964 3.6938 0.9688 1.7643 3.8451

0.0965 0.1010 0.1023 0.1008 0.0951 0.4207 0.5616 0.4694 0.5072 0.5035 121.47 130.46 135.09 132.53 125.72

6.3

2.7

The modified ABC

(3.5) (1.0) (2.3) (0.8) (4.9) (15.9) (12.3) (6.1) (1.4) (0.7) (6.6) (0.4) (3.9) (1.9) (3.3)

Std.

Mean

Std.

0.0071 0.0041 0.0018 0.0033 0.0157 0.1215 0.1029 0.0881 0.0593 0.0192 12.5423 2.2842 3.4805 2.3086 6.1563

0.1032 (3.2) 0.1012 (1.2) 0.1002 (0.2) 0.1013 (1.3) 0.1018 (1.8) 0.5154 (3.1) 0.5198 (3.9) 0.4968 (0.6) 0.5044 (0.9) 0.5042 (0.8) 134.53 (3.5) 132.69 (2.1) 130.74 (0.6) 130.83 (0.6) 131.73 (1.3)

0.0059 0.0033 0.0005 0.0027 0.0050 0.0591 0.0575 0.0271 0.0235 0.0340 8.3124 5.4659 1.4705 1.7663 4.4569

15.4

3.9

Note: The values in parentheses are absolute relative errors in %. Std. represents standard deviation.

Fig. 7. The normalized parameter is defined as the ratio of the identified parameter value to the true value. In Fig. 7, the PDFs of the normalized parameters are obtained from IR independent runs in the unknown mass case with partial output and 10% noise corruption. It is seen that the normalized parameters approximately follow normal distributions with the mean values very close to 1 and small standard deviations, which verifies that the identified parameters are very close to the true values.

where a and b are constant coefficients, fr is modal damping ratio associated with the ith mode (set as 5% in the first and second vibration modes (r = 1, 2)), and xr is the rth natural frequency. In this example, the masses, modal damping ratios (used for determining the damping coefficient parameters a and b), and stiffness values are the unknown parameters to be identified. The parameter vectors h1, with twenty two parameters, and h2, with forty two, are respectively defined to describe the known mass and the unknown mass cases as follows:

4.2. A 20-DOF linear system To further compare the performance of the modified version of ABC with other heuristic methodologies (i.e., GAs, PSO, and DE) and to verify its applicability to SI as the number of parameters to be identified increases, a two-dimensional 20-DOF shear type linear system is analyzed, with properties given in Table 8. This system was used by Perry et al. [55] and Tang et al. [42] to test GAs and DE, respectively. The equations of motion are those shown in Eq. (11) with the Rayleigh damping matrix expressed as:

½C ¼ a½M þ b½K;

fr ¼

a 2xr

þ

bxr 2

ð14Þ

h1 ¼ fk1 ; . . . ; k20 ; f1 ; f2 g

ð15Þ

h2 ¼ fm1 ; . . . ; m20 ; k1 ; . . . ; k20 ; f1 ; f2 g

ð16Þ

Input forces are assumed to be Gaussian noise sequences with the RMS scaled to 1 kN, a duration of 10 s, and sampling time T = 0.01 s; they are applied at every five floors of the structure, in exactly the same way described in [55] and [42]. Accelerations are measured at 40% of the floors for the known mass case and at 60% of the floors for the unknown mass case as detailed in Table 9. The ABC parameters used in this example are provided in Table 10.

H. Sun et al. / Computers and Structures 116 (2013) 59–74

Fig. 6. A comparison of simulation results for the unknown mass case with partial measurements: (a) noise free and (b) 10% noise corruption in I/O signals.

Fig. 7. The PDFs of the normalized parameters obtained from IR independent runs in the unknown mass case with partial output and 10% noise corruption.

67

68

H. Sun et al. / Computers and Structures 116 (2013) 59–74

Table 8 The structural properties of a 20-DOF linear system. Mass (kg)

Levels 1–10 Levels 11–20

4000 3000

Stiffness (kN/m)

Levels 1–10 Levels 11–15 Levels 16–20

5000 4000 3500

Natural period of vibration (s)

First mode Second mode

2.123 0.797

Table 9 Location of acceleration measurements. Mass cases

Floor levels

Known mass Unknown mass

2, 4, 7, 10, 12, 14, 17, 20 1, 2, 3, 4, 6, 8, 10, 12, 14, 16, 18, 20

Table 10 Parameter values of the modified ABC algorithm for system identification. Mass cases

np

miter

lc

m

d

IR

Known Unknown

90 180

500 500

300 300

4 5

100 100

20 20

Table 11 Results for 20-DOF known mass system without noise corruption. Parameter error

GA [55] SSRM [55] PSO [42] DE [42] The modified ABC

Mean error-k (%) Max. error-k (%) Mean error-c (%) Max. error-c (%)

8.33 31.28 15.81 28.97

0.52 1.60 0.64 1.21

0.71 3.37 2.24 8.31

0.41 1.29 0.53 1.45

0.19 1.04 0.19 1.06

Note: The values in parentheses are absolute relative errors in %.

The identified results for the 20-DOF system are presented in Tables 11 and 12 in comparison with those obtained by GA, PSO, DE and SSRM [42,55]. The SSRM is a search space reduction method integrated with a modified GA. Figs. 8 and 9 show convergences of the mean results obtained via the modified ABC algorithm for the known and the unknown mass cases, respectively. For the known mass case without noise corruption, the modified ABC provides much better parameter estimation than the GA and PSO, yielding much smaller errors (a maximum value of 1.06% in damping parameter estimation), while its performance is comparable to those of DE and SSRM. From the convergence graphs shown in Fig. 8, it is observed that all the parameter estimates stabilize after 200 iterations. Moreover, the estimated values of damping ratios are quite accurate. For the unknown mass case with noise corruption (5% and 10% RMS) in I/O signals, the modified ABC algorithm performs better than the PSO and as well as the DE and the SSRM, with the maximum errors of mass and stiffness obtained by ABC being slightly

bigger than those obtained via the DE and the SSRM. On the other hand, the modified ABC algorithm seems to less sensitive to the presence of noise in the damping ratio estimation. In general, the maximum errors among all the identified parameters obtained by ABC (5.48% for 5% noise corruption and 12.26% for 10% noise corruption) are quite good when compared with those obtained by other methods. The search performance shown in Fig. 9 presents a good estimation procedure of the two damping ratios f1 and f2, converging to the correct values after 300 iterations. Compared with the results discussed by Tang et al. [42], the proposed ABC algorithm seems to be more efficient than the DE algorithm: when applied to this 20-DOFs case, the DE algorithm uses an initial population size of 400 while the proposed ABC algorithm runs with a much smaller initial population size (90 for the known mass case and 180 for the unknown mass case). The population size has a strong impact on the computational effort of the algorithm since it is linked to the number of fitness function evaluations which involve the time-consuming operation of solving the 2nd order ordinary differential equations of motion. For an initial population of 400, the DE algorithm will require about 200,400 fitness evaluations while the ABC algorithm, for initial populations of 90 (known mass case) and 180 (unknown mass case), will require about 45,045 and 95,095 evaluations, respectively, with a substantial reduction in computational time. In terms of the CPU time used for each run by the proposed algorithm coded in MATLABÒ, it took on average 18.5 min and 41.6 min for each single simulation for the known and unknown mass case, respectively. The computational cost can be significantly decreased if the scheme is set to run in parallel on a multi-processor unit coded in more efficient computer languages (i.e., C, C++, etc.). Fig. 10 shows several typical PDFs of the normalized mass and stiffness parameters obtained from IR independent runs in the unknown mass case with 10% noise corruption. It is concluded from this figure that the identified structural parameters are quite accurate. 4.3. A 2-DOF nonlinear system To test the versatility of the modified ABC algorithm in handling nonlinear problems, a known mass 2-DOF shear type hysteretic nonlinear system, considered by Xue et al. [54], is employed as the third numerical example here, as shown in Fig. 11. The equation of motion of the system can be written as:

€ g ðtÞ €ðtÞ þ ½CxðtÞ _ ½Mx þ rðtÞ ¼ ½MIu

ð17Þ

where [M] and [C] are the 2  2 mass and damping matrices; x; x_ € are relative displacement, velocity, and acceleration vectors and x to the ground, respectively; I is the ground motion influence coeffi€ g ðtÞ is the earthquake acceleration; cient vector with element 1; u and r = (r1  r2, r2)T is the vector containing the restoring forces at each floor expressed by the two following expressions [32].

r_ 1 ¼ ½k1  a1 r 1 sgnðx_ 1 Þjr 1 jn1 1  b1 jr1 jn1 x_ 1 r_ 2 ¼ ½k2  a2 sgnðx_ 2  x_ 1 Þr 2 jr 2 j

n2 1

ð18Þ n2

 b2 jr2 j ðx_ 2  x_ 1 Þ

ð19Þ

Table 12 Results for 20-DOF unknown mass system with noise corruption. Parameter error

Mean error-m (%) Max. error-m (%) Mean error-k (%) Max. error-k (%) Mean error-c (%) Max. error-c (%)

5% Noise corruption

10% Noise corruption

SSRM [55]

PSO [42]

DE [42]

The modified ABC

SSRM [55]

PSO [42]

DE [42]

The modified ABC

1.51 4.02 1.38 3.83 6.70 12.90

3.61 10.81 3.65 8.13 10.34 16.57

1.42 3.56 1.27 4.11 7.23 10.68

2.31 5.16 2.27 5.48 1.93 5.41

3.00 10.40 2.78 8.64 14.69 20.36

7.06 16.27 5.31 14.36 17.31 29.06

3.29 11.21 2.63 9.02 13.54 21.04

3.95 9.72 4.58 11.69 3.83 12.26

Note: The values in parentheses are absolute relative errors in %.

H. Sun et al. / Computers and Structures 116 (2013) 59–74

69

Fig. 8. Simulation results for the known mass case with free noise and 40% measurements.

Fig. 9. Simulation results for the unknown mass case with 60% measurements and 10% noise corruption in I/O signals.

This nonlinear model was clarified as the Bouc–Wen model [58,59]. In above two equations, x_ 1 and x_ 2 are the velocities of the first and second masses, respectively, k1 and k2 are stiffness parameters of the vertical elements, a1, a2, b1, b2, n1 and n2 are nonlinear coefficient parameters, while sgn() is signum function (i.e., sgn(#) gives 1, 0, or 1 depending on whether # is negative, zero, or positive, respectively). Hence, the system parameters to be identified can be grouped in a set of 10 variables as:

h ¼ fc1 ; c2 ; k1 ; k2 ; a1 ; a2 ; b1 ; b2 ; n1 ; n2 g

ð20Þ

The values of such parameters used in the numerical analyses are shown in Table 13. The ElCentro ARR4S40E (1940) earthquake acceleration data is used as input to the system with a sampling frequency of 100Hz (DT = 0.01 s) for the duration of 10 s. The maximum peak value of the acceleration is scaled to 3 g and the signal is filtered with

70

H. Sun et al. / Computers and Structures 116 (2013) 59–74

a high frequency cutoff of 30 Hz and a low frequency cutoff of 0.1 Hz. In order to study the effects of incomplete sets of measurements, two cases are considered as shown in Fig. 11: (1) the full output case in which a full set of measurements was available, and (2) the partial output case with only one output on the 2nd DOF. Similar to the cases in [54], two noise levels 0% and 5% were considered to emphasize the effect of measurement noise on the identification of the parameters. The ABC parameters are set as follows: np = 30, miter = 500, lc = 300, m = 4, d = 100 and IR = 20. To compare the performance of the ABC algorithm for structural SI with other heuristic methodologies, the enhanced particle swarm optimization (EPSO) method proposed by Charalampakis and Dimou [39] was used: the inertia weight factor was set as 0.6 and the two acceleration factors as 1.7 suggested in [39]; the population size, maximum iteration number and number of independent runs were the same as those used in ABC. The mean values and standard deviations are summarized in Table 14, for the full output scenario, and in Table 15, for the partial output scenario. Looking at both instrumentation setups, it can be seen that the errors obtained using ABC range from 0% to 2.4%, in the noise free case, and from 0.04% to 9.2% in the noise polluted case, while, using EPSO, the corresponding errors range from

0.1% to 5.2%, and from 0.4% to 19.9%. The results show that the modified ABC approach provides almost perfect identification in the absence of measurement noise, as the EPSO. However, in the noise polluted case, the modified ABC algorithm performs significantly better than EPSO which appears to be more sensitive to noise corruption. From further analysis, it appears that the ABC algorithm is able to yield good identification results with small errors, even at higher noise levels. For all cases, the standard deviations obtained by the modified ABC are smaller than those by EPSO. It is noteworthy that the maximum errors usually exist in the estimation of the nonlinear parameters ai, bi and ni (i = 1,2), and the signal noise corruption has larger impact on the estimation of such parameters than on the estimation of damping and stiffness parameters. Fig. 12 presents the hysteresis loops of the first floor element for the noise free and noise polluted cases in the partial output scenario. It can be seen that there is a nearly perfect agreement between the simulated and ABC estimated results in the phase-plane plots of the Bouc–Wen parameter r1 vs. the displacement x1, even with partial outputs and at a high noise corruption level. Fig. 13 shows a comparison of the identified normalized parameters for both the noise free and the noise corrupted measurements in the partial output case. It is observed that the

Fig. 10. Typical PDFs of the normalized parameters obtained from IR independent runs in the unknown mass case with 60% measurements and 10% noise corruption.

Fig. 11. A 2-DOF nonlinear system with ground motion excitation: (a) full measured acceleration outputs and (b) one measured acceleration output.

Table 13 The structural properties of a 2-DOF nonlinear system. Parameters

Values

Damping (kNs/m)

Stiffness (kN/m)

Nonlinear parameters

m1

Mass (kg) m2

c1

c2

k1

k2

a1

a2

b1

b2

n1

n2

100

80

0.55

0.50

30

24

1

2

2

1

3

2

71

H. Sun et al. / Computers and Structures 116 (2013) 59–74

estimation of parameters rapidly converges and is able to reach good results after about 200 iterations. Fig. 14 illustrates the PDFs of the normalized parameters obtained from IR independent runs in the partial output scenario with 5% noise corruption. Similar to the first two examples, the PDFs approximately follow normal distribution. The mean values of these normalized parameters are quite close to 1 accounting for good estimates of the system parameters.

4.4. Observations of the numerical examples The simulation results previously discussed lead to the following observations:  From the simulation results, it can be concluded that the identification of unknown mass system is much more challenging than that of known mass system, because the unknown mass

Table 14 Results for 2-DOF known mass system in the full output scenario. Parameters

True value

0% Noise

5% Noise

EPSO

The modified ABC

Mean c1 c2 k1 k2

a1 a2 b1 b2 n1 n2

0.55 0.50 30 24 1 2 2 1 3 2

Max error (%)

0.5529 0.4983 29.969 24.026 0.9533 2.0576 1.9451 1.0375 3.0609 2.0260

(0.5) (0.3) (0.1) (0.1) (4.7) (2.9) (2.7) (3.7) (2.0) (1.3)

EPSO

The modified ABC

Std.

Mean

Std.

Mean

0.0056 0.0043 0.1499 0.1578 0.1327 0.1420 0.2259 0.1432 0.1742 0.1341

0.5501 (0.02) 0.5002 (0.04) 30.011 (0.04) 23.997 (0.01) 1.0007 (0.1) 1.9839 (0.8) 2.0066 (0.33) 1.0035 (0.35) 2.9979 (0.1) 2.0152 (0.8)

0.0020 0.0019 0.0431 0.0162 0.0484 0.0886 0.0809 0.0213 0.0610 0.0975

0.5606 0.4771 30.111 23.744 0.9710 2.1955 2.0607 0.8470 3.0248 2.3338

4.7

0.8

(1.9) (4.6) (0.4) (1.1) (2.9) (9.8) (3.1) (15.3) (0.8) (13.3)

Std.

Mean

0.0040 0.0030 0.1437 0.0856 0.1409 0.0698 0.2853 0.1062 0.1925 0.0869

0.5488 0.4809 29.879 24.125 0.9644 2.0946 1.8516 0.9307 3.0764 2.1119

15.3

Std. (0.2) (3.8) (0.4) (0.5) (3.6) (4.7) (7.4) (6.9) (2.5) (5.6)

0.0022 0.0016 0.0476 0.0102 0.0385 0.0856 0.0738 0.0151 0.0539 0.0611

7.4

Note: The values in parentheses are absolute relative errors in %. Std. represents standard deviation.

Table 15 Results for 2-DOF known mass system in the partial output scenario. Parameters

True value

0% Noise

5% Noise

EPSO

The modified ABC

Mean c1 c2 k1 k2

a1 a2 b1 b2 n1 n2 Max error (%)

0.55 0.50 30 24 1 2 2 1 3 2

0.5546 0.4973 30.063 23.926 1.0024 2.0963 2.0749 0.9481 2.9710 2.1009 5.2

(0.8) (0.5) (0.2) (0.3) (0.2) (4.8) (3.7) (5.2) (1.0) (5.0)

EPSO

Std.

Mean

Std.

Mean

0.0036 0.0027 0.1048 0.1210 0.0853 0.1097 0.3468 0.1397 0.1014 0.1833

0.5527 (0.5) 0.5000 (0.0) 30.018 (0.1) 23.988 (0.1) 0.9755 (2.4) 1.9837 (0.8) 2.0003 (0.02) 1.0025 (0.3) 3.0264 (0.9) 2.0476 (2.4)

0.0073 0.0024 0.1709 0.1380 0.1206 0.0859 0.2437 0.1271 0.1646 0.1294

0.5383 0.4971 29.437 24.624 0.9543 1.6938 1.7163 0.9679 3.1661 2.3982

2.4

19.9

The modified ABC

(2.1) (0.6) (1.9) (2.6) (4.6) (15.3) (14.2) (3.3) (5.5) (19.9)

Std.

Mean

0.0078 0.0071 0.0879 0.0277 0.1026 0.2295 0.1802 0.1393 0.1392 0.2260

0.5608 0.5002 30.059 24.108 0.9389 2.0673 1.9143 1.0923 3.0988 1.9162 9.2

Note: The values in parentheses are absolute relative errors in %. Std. represents standard deviation.

Fig. 12. Hysteresis loops of the first DOF with partial measurements: (a) noise free and (b) 5% noise corruption in I/O signals.

Std. (2.0) (0.04) (0.2) (0.4) (6.1) (3.4) (4.3) (9.2) (3.3) (4.2)

0.0061 0.0048 0.1120 0.1525 0.1041 0.1475 0.1715 0.1788 0.1374 0.1922

72

H. Sun et al. / Computers and Structures 116 (2013) 59–74

Fig. 13. A comparison of simulation results with partial measurements: (a) noise free and (b) 5% noise corruption in I/O signals.









system identification is always a high multimodal problem with complex nonlinearities and convexities while the known mass system has a reduction of search space. The ABC algorithm has proven to be very efficient for structural parameter identification: it converges relatively quickly, it is robust, and yields small deviations of parameter estimation in each independent run. All the algorithms considered in this study (ABC, GA, EPSO, etc.) yield good estimation results with small errors for noise free cases in full output scenarios. The modified ABC algorithm, however, performs more satisfactorily than others when dealing with noisy I/O signals, unknown masses and partial output information. The parameters obtained by the modified version of ABC algorithm converge quite rapidly the correct solutions for known mass cases. For complex unknown mass systems, larger population sizes or more iteration cycles can be employed in order to find finer solutions. The modified ABC algorithm seems to be relatively insensitive to noise. It could provide accurate parameter estimations even for complex systems with partial measurements at high noise corruption levels (i.e., 10% RMS noise in I/O signals).

 The maximum error mostly occurs in the identification of damping parameters and/or nonlinear parameters. This phenomenon is enhanced if the measured signals are corrupted by noise. 5. Conclusions A modified version of a novel swarm intelligence technique called the Artificial Bee Colony (ABC) algorithm has been presented in the context of structural system identification. The ABC is a population based stochastic algorithm which requires only few common control parameters. It has a simple structure, it is easy to implement, it is robust, and it provides results with high accuracy. The application of proposed nonlinear factor for convergence control enhances the balance of global and local searches. The idea behind the nonlinear factor is to free the method, at its initial steps, to search over the entire search space and to allow for faster convergence during the local search phase which takes place towards the end of the iteration process. In this way, the search space is gradually reduced following a nonlinear trend. To avoid stagnation around local minima, a ‘‘scout bee’’ search will be carried out if needed. Thus, the modified version of ABC algorithm is supposed to be able to locate finer solutions. To investigate the applicability

H. Sun et al. / Computers and Structures 116 (2013) 59–74

73

Fig. 14. The PDFs of the normalized parameters obtained from IR independent runs in the partial output scenario with 5% noise corruption.

of this proposed technique to structural system identification, two linear systems (5-DOF and 20-DOF) and a nonlinear system (2DOF) have been studied under different conditions, addressing issues such as the number of measurements used in identification, noise in the signals, and the knowledge of the structural mass. In all cases considered, the simulation results show that the proposed ABC algorithm can produce excellent parameter estimation with small errors. The presented method is effective, robust and efficient even with reduced partial measurements at high noise pollution levels.

Acknowledgements This study was supported by the Department of Civil Engineering and Engineering Mechanics, Columbia University, whose support is greatly acknowledged. The authors also would like to thank Mr. Zifeng Yuan for his many discussions and suggestions in the development of the algorithm.

References [1] Ibáñez P. Identification of dynamic parameters of linear and non-linear structural models from experimental data. Nucl Eng Des 1973;25(1):30–41. [2] Koh CG, See LM. Identification and uncertainty estimation of structural parameters. J Eng Mech 1994;120(6):1219–36. [3] Lin J-W et al. On-line identification of non-linear hysteretic structural systems using a variable trace approach. Earthq Eng Struct Dyn 2001;30(9):1279–303. [4] Smyth AW et al. Development of adaptive modeling techniques for non-linear hysteretic systems. Int J Nonlinear Mech 2002;37(8):1435–51. [5] Chatzi EN, Smyth AW. The unscented Kalman filter and particle filter methods for nonlinear structural system identification with non-collocated heterogeneous sensing. Struct Contr Health Monit 2009;16(1):99–123. [6] Arslan Ö, Aykan M, Nevzat Özgüven H. Parametric identification of structural nonlinearities from measured frequency response data. Mech Syst Signal Pr 2011;25(4):1112–25. [7] De Angelis M et al. Extracting physical parameters of mechanical models from identified state-space representations. J Appl Mech 2002;69(5):617–25. [8] Lusß H, Betti R, Longman RW. Identification of linear structural systems using earthquake-induced vibration data. Earthq Eng Struct Dyn 1999;28(11):1449–67.

[9] Lusß et al. Constructing second-order models of mechanical systems from identified state space realizations. Part I: Theoretical discussions. J Eng Mech 2003;129(5):477–88. [10] Lusß et al. Constructing second-order models of mechanical systems from identified state space realizations. Part II: Numerical investigations. J Eng Mech 2003;129(5):489–501. [11] Loh C-H, Tsaur Y-H. Time domain estimation of structural parameters. Eng Struct 1988;10(2):95–105. [12] McVerry GH. Structural identification in the frequency domain from earthquake records. Earthq Eng Struct Dyn 1980;8(2):161–80. [13] Pintelon R, Schoukens J, Guillaume P. Parametric frequency domain modeling in modal analysis. Mech Syst Signal Pr 1989;3(4):389–403. [14] Hong K-S, Yun C-B. Improved method for frequency domain identifications of structures. Eng Struct 1993;15(3):179–88. [15] Mottershead JE, Stanway R. Identification of structural vibration parameters by using a frequency domain filter. J Sound Vib 1986;109(3):495–506. [16] Lee U, Shin J. A frequency-domain method of structural damage identification formulated from the dynamic stiffness equation of motion. J Sound Vib 2002;257(4):615–34. [17] Farrar CR, Doebling SW. Lessons learned from applications of vibration-based damage identification methods to a large bridge structure, In: Proceedings of the international workshop on structural health monitoring. Stanford, CA: Stanford University; 1997. p. 351–70. [18] Franco G et al. Identification of structural systems using an evolutionary strategy. J Eng Mech 2004;130(10):1125–39. [19] Wang GS. Application of hybrid genetic algorithm to system identification. Struct Contr Health Monit 2009;16(2):125–53. [20] Yang JNJN, Lin S. On-line identification of non-linear hysteretic structures using an adaptive tracking technique. Int J Nonlinear Mech 2004;39(9):1481–91. [21] Yang JN, Lin S. Identification of parametric variations of structures based on least squares estimation and adaptive tracking technique. J Eng Mech 2005;131(3):290–8. [22] Tang H-S et al. Online weighted LS-SVM for hysteretic structural system identification. Eng Struct 2006;28(12):1728–35. [23] Yang JN, Huang H, Lin S. Sequential non-linear least-square estimation for damage identification of structures. Int J Nonlinear Mech 2006;41(1):124–40. [24] Hoshiya M, Saito E. Structural identification by extended Kalman filter. J Eng Mech 1984;110(12):1757–70. [25] Jeen-Shang L, Yigong Z. Nonlinear structural identification using extended Kalman filter. Comput Struct 1994;52(4):757–64. [26] Wang D, Haldar A. System identification with limited observations and without input. J Eng Mech 1997;123(5):504–11. [27] Yang JN et al. An adaptive extended Kalman filter for structural damage identification. Struct Contr Health Monit 2006;13(4):849–67. [28] Wu ML, Smyth AW. Application of the unscented Kalman filter for real-time nonlinear structural system identification. Struct Contr Health Monit 2007;14(7):971–90.

74

H. Sun et al. / Computers and Structures 116 (2013) 59–74

[29] Mariani S, Ghisi A. Unscented Kalman filtering for nonlinear structural dynamics. Nonlinear Dynam 2007;49(1–2):131–50. [30] Xie Z, Feng J. Real-time nonlinear structural system identification via iterated unscented Kalman filter. Mech Syst Signal Pr 2012;28:309–22. [31] Nasrellah HA, Manohar CS. A particle filtering approach for structural system identification in vehicle-structure interaction problems. J Sound Vib 2010;329(9):1289–309. [32] Sato T, Qi K. Adaptive H1 filter: its application to structural identification. J Eng Mech 1998;124(11):1233–40. [33] Williams P. Structural damage detection from transient responses using square-root unscented filtering. Acta Astronaut 2008;63(11–12):1259–72. [34] Li SJ, Suzuki Y, Noori M. Improvement of parameter estimation for non-linear hysteretic systems with slip by a fast Bayesian bootstrap filter. Int J Nonlinear Mech 2004;39(9):1435–45. [35] Ching J, Beck JL, Porter KA. Bayesian state and parameter estimation of uncertain dynamical systems. Probab Eng Mech 2006;21(1):81–96. [36] Lin J-W, Betti R. On-line identification and damage detection in non-linear structural systems using a variable forgetting factor approach. Earthq Eng Struct Dyn 2004;33(4):419–44. [37] Koh CG, Chen YF, Liaw CY. A hybrid computational strategy for identification of structural parameters. Comput Struct 2003;81(2):107–17. [38] Kwok NM et al. A novel hysteretic model for magnetorheological fluid dampers and parameter identification using particle swarm optimization. Sensor Actuat A: Phys 2006;132(2):441–51. [39] Charalampakis AE, Dimou CK. Identification of Bouc–Wen hysteretic systems using particle swarm optimization. Comput Struct 2010;88(21–22):1197–205. [40] Lefik M, Schrefler BA. Artificial neural network for parameter identifications for an elasto-plastic model of superconducting cable under cyclic loading. Comput Struct 2002;80(22):1699–713. [41] Karimi I et al. System identification of concrete gravity dams using artificial neural networks based on a hybrid finite element-boundary element approach. Eng Struct 2010;32(11):3583–91. [42] Tang H, Xue S, Fan C. Differential evolution strategy for structural system identification. Comput Struct 2008;86(21–22):2004–12. [43] Karaboga D. An idea based on honey bee swarm for numerical optimization. In: Technical Report-TR06, Erciyes University, Engineering Faculty, Computer Engineering Department, Turkey, 2005.

[44] Karaboga D, Basturk B. A powerful and efficient algorithm for numerical function optimization: Artificial Bee Colony (ABC) algorithm. J Glob Optim 2007;39(3):459–71. [45] Karaboga D, Basturk B. On the performance of Artificial Bee Colony (ABC) algorithm. Appl Soft Comput 2008;8(1):687–97. [46] Karaboga D, Akay B. A comparative study of Artificial Bee Colony algorithm. Appl Math Comput 2009;214(1):108–32. [47] Singh A. An Artificial Bee Colony algorithm for the leaf-constrained minimum spanning tree problem. Appl Soft Comput 2009;9(2):625–31. [48] Karaboga D, Akay B. A modified Artificial Bee Colony (ABC) algorithm for constrained optimization problems. Appl Soft Comput 2011;11(3):3021–31. [49] Pan Q-K et al. A discrete Artificial Bee Colony algorithm for the lot-streaming flow shop scheduling problem. Inform Sci 2011;181(12):2455–68. [50] Kang F, Li J, Xu Q. Structural inverse analysis by hybrid simplex Artificial Bee Colony algorithms. Comput Struct 2009;87(13–14):861–70. [51] Sonmez M. Artificial Bee Colony algorithm for optimization of truss structures. Appl Soft Co 2011;11(2):2406–18. [52] Sonmez M. Discrete optimum design of truss structures using Artificial Bee Colony algorithm. Struct Multidiscip Optim 2011;43(1):85–97. [53] Omkar SN et al. Artificial Bee Colony (ABC) for multi-objective design optimization of composite structures. Appl Soft Comput 2011;11(1):489–99. [54] Xue S, Tang H, Zhou J. Identification of structural systems using particle swarm optimization. J Asian Archit Build Eng 2009;8(2):517–24. [55] Perry MJ, Koh CG, Choo YS. Modified genetic algorithm strategy for structural identification. Comput Struct 2006;84(8–9):529–40. [56] Phan M et al. Improvement of Observer/Kalman filter identification (OKID) by residual whitening. J Vib Acoust 1995;117(2):232–9. [57] Hinterding R. Gaussian mutation and self-adaption for numeric genetic algorithms. In: IEEE international conference on evolutionary computation, Perth, WA, Australia, 1995. [58] Bouc R. Forced vibration of mechanical systems with hysteresis. In: Proceedings of the fourth conference on nonlinear oscillation, Prague, Czechoslovakia, 1967. [59] Wen YK. Method for random vibration of hysteretic systems. ASCE J Eng Mech Div 1976;102(2):249–63.