EARTHQUAKE ENGINEERING AND STRUCTURAL DYNAMICS Earthquake Engng Struct. Dyn. 2005; 34:665–685 Published online 14 January 2005 in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/eqe.449 Multi-objective framework for structural model identi cation Yiannis Haralampidis, Costas Papadimitriou∗; † and Maria Pavlidou Department of Mechanical and Industrial Engineering; University of Thessaly; Pedion Areos; 38334 Volos; Greece SUMMARY Structural identi cation based on measured dynamic data is formulated in a multi-objective context that allows the simultaneous minimization of the various objectives related to the t between measured and model predicted data. Thus, the need for using arbitrary weighting factors for weighting the relative importance of each objective is eliminated. For con icting objectives there is no longer one solution but rather a whole set of acceptable compromise solutions, known as Pareto solutions, which are optimal in the sense that they cannot be improved in any objective without causing degradation in at least one other objective. The strength Pareto evolutionary algorithm is used to estimate the set of Pareto optimal structural models and the corresponding Pareto front. The multi-objective structural identi cation framework is presented for linear models and measured data consisting of modal frequencies and modeshapes. The applicability of the framework to non-linear model identi cation is also addressed. The framework is illustrated by identifying the Pareto optimal models for a scaled laboratory building structure using experimentally obtained modal data. A large variability in the Pareto optimal structural models is observed. It is demonstrated that the structural reliability predictions computed from the identi ed Pareto optimal models may vary considerably. The proposed methodology can be used to explore the variability in such predictions and provide updated structural safety assessments, taking into consideration all Pareto structural models that are consistent with the measured data. Copyright ? 2005 John Wiley & Sons, Ltd. KEY WORDS: structural identi cation; multi-objective optimization; Pareto set; reliability 1. INTRODUCTION The problem of identifying the parameters of a structural model using dynamic data has received much attention over the years because of its importance in structural model updating, structural health monitoring and structural control. Comprehensive reviews of structural ∗ Correspondence to: Costas Papadimitriou, Department of Mechanical and Industrial Engineering, University of Thessaly, Pedion Areos, 38334 Volos, Greece. † E-mail: costasp@uth.gr Contract=grant sponsor: Greek Secretariat of Research and Technology Contract=grant sponsor: Greek Earthquake Planning and Protection Organization (OASP) Copyright ? 2005 John Wiley & Sons, Ltd. Received 9 July 2003 Revised 12 March 2004 Accepted 15 October 2004 666 Y. HARALAMPIDIS, C. PAPADIMITRIOU AND M. PAVLIDOU parameter identi cation methods can be found in References [1, 2]. The estimate of the parameter values involves uncertainties that are due to limitations of the mathematical models used to represent the behavior of the real structure, the presence of measurement error in the data, and insucient excitation and response bandwidth. Structural identi cation and nite element model updating methodologies, for example, References [3–8], are often based on modal data since these data are readily obtained from well-established experimental structural dynamics techniques based on either forced [9, 10] or ambient vibration tests [11, 12]. The optimal structural models resulting from such methods can be used for response and reliability predictions, structural health monitoring and control. Parameter identi cation problems based on measured data are often formulated as leastsquares problems in which objective functions measuring the t between measured and model predicted data are built up into a single objective using weighting factors. Standard optimization techniques are then used to nd the optimal values of the parameters that minimize the overall measure of t. The results of the optimization depend on the weighting factors assumed. The choice of the weighting factors depends on the model adequacy and the uncertainty in the available measured data, which are not known a priori. In this work, the parameter identi cation problem is formulated in a multi-objective context that allows the simultaneous minimization of the multiple objectives, eliminating the need for using arbitrary weighting factors for weighting the relative importance of each objective. The set of admissible solutions are known in multi-objective optimization terminology as Pareto optimal solutions. The characteristics of the Pareto solutions are that they cannot be improved in any objective without causing degradation in one other objective. The set of Pareto solutions can be obtained using Evolutionary Algorithms [13] well-suited to solve multi-objective optimization problems [14, 15]. A multi-objective parameter identi cation framework for linear structures based on experimentally obtained modal data is presented. For this, the measured modal properties are grouped according to their characteristics (mode type, modal frequencies and modeshapes) and their signi cance on the objective of the identi cation. For each group, an objective function is introduced measuring the mismatch between the measured and the model predicted modal properties involved in the group. The Pareto optimal structural models are identi ed by minimizing these objectives simultaneously. The multi-objective problem is solved using a recently proposed Strength Pareto Evolutionary Algorithm [16] capable of identifying and representing the Pareto solutions with a few points uniformly distributed along the Pareto front [17]. The identi ed Pareto optimal structural models constitute acceptable compromise solutions trading-o the quality of t in di erent groups of modal properties. The proposed methodology is illustrated by identifying the Pareto front and the corresponding set of Pareto optimal solutions for a model of a scaled laboratory structure using experimentally obtained modal data. These models are then used to predict the reliability of the structure to future stochastic loads modeled as stationary white noise processes. It is demonstrated that such predictions computed from the identi ed models along the Pareto optimal front may vary considerably. The presentation in this work is organized as follows. In Section 2 the structural identication problem using measured modal data is formulated as a multi-objective optimization problem. The main concepts of evolutionary strategies and a recently proposed Strength Pareto Evolutionary Algorithm [16] for solving multi-objective minimization problems are outlined in Section 3. Applications on the multi-objective structural model identi cation methodology using measured data from a laboratory scaled building structure are presented in Section 4. Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 STRUCTURAL MODEL IDENTIFICATION 667 The variability in the predictions of structural reliability obtained by the identi ed Pareto optimal structural models is also addressed. For completeness, a discussion for extending the methodology to non-linear models based on measured response time histories, instead of modal properties, is presented in Section 5. The conclusions are summarized in Section 6. 2. MULTI-OBJECTIVE IDENTIFICATION BASED ON MODAL DATA 2.1. Formulation Let !ˆ r and ˆ r ∈ RN0 (r = 1; : : : ; M ) be measured modal data from a structure, consisting of modal frequencies !ˆ r and modeshape components ˆ r at N0 measured DOFs, where M is the number of observed modes. Consider a parameterized class of linear structural models (e.g. nite element models) used to model the dynamic behavior of the structure. Let x be the set of free model parameters to be identi ed using the measured modal data. These parameters are usually associated with geometrical, material, sti ness or mass properties and boundary conditions. The modal frequencies !r (x) and modeshapes r (x) predicted from the model class for a particular value of the parameter set x are obtained by solving the following eigenvalue problem [K(x) − !r2 (x)M (x)]r (x) = 0 (1) where M (x) and K(x) are the global model mass and sti ness matrices, respectively. The objective in a modal-based structural identi cation methodology is to estimate the values of the parameter set x so that the modal data {!r (x); r (x); r = 1; : : : ; M } predicted by the linear class of models best matches, in some sense, the experimentally obtained modal data {! ˆ r ; ˆ r ; r = 1; : : : ; M }. For this, the measured modal properties are grouped into n groups. Each group contains one or more modal properties. For the i-th group, a norm Ji (x) is introduced to measure the mismatch between the measured modal properties involved in the group and the corresponding modal properties predicted from the model class for a particular value of the parameter set x. The di erence between the measured modal data and the model-based modal predictions is due to modeling and measurement errors always present in structural identi cation problems. The problem of identifying the model parameter values that give the best t in all groups of modal properties is formulated as a multi-objective optimization problem stated as follows. Find the values of the parameter set x that simultaneously minimizes the objectives y = J (x) = (J1 (x); J2 (x); : : : Jn (x)) (2) where, using multi-objective terminology, x = (x1 ; :::::::; xm ) ∈ X is the parameter vector, X is the parameter space, y = (y1 ; :::::::; yn ) ∈ Y is the objective vector, and Y is the objective space. The optimization may be constrained due to restrictions imposed on the parameter set or to other constraints involved in the formulation. The feasible parameter space is usually con ned in a hypercube by specifying lower and upper limits on each parameter. These limits depend on physical constraints, information about the physical characteristics of the system and modeling experience. Moreover, the search space can be con ned to regions in the objective space for which the t in the various modal properties are below a threshold value since solutions with high mismatch values in such properties are usually unacceptable. Constrained multi-objective Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 668 Y. HARALAMPIDIS, C. PAPADIMITRIOU AND M. PAVLIDOU optimization problems are posed by considering the parameter space X to be con ned in a region speci ed by the constraints. For con icting objectives J1 (x); : : : ; Jn (x), there is no single optimal solution, but rather a set of alternative solutions which are optimal in the sense that no other solutions in the search space are superior to them when all objectives are considered. Such alternative solutions, trading-o the ts in di erent modal group properties, are known in multi-objective optimization as Pareto optimal solutions. Next, useful de nitions related to the dominated and non-dominated vectors and Pareto solutions appearing in the multi-objective terminology [14, 15] are introduced. Speci cally, a decision vector a ∈ X is said to dominate a decision vector b ∈ X (also written as a ≺ b) if and only if Ji (a)6Ji (b) ∀i ∈ {1; : : : ; n} and ∃j ∈ {1; : : : ; n} : Jj (a)¡Jj (b) (3) Additionally, we say a covers b (a  b) if and only if a ≺ b or J (a) = J (b). Based on the above relation, we can de ne non-dominated and Pareto optimal solutions as follows. A decision vector a ∈ X is said to be non-dominated regarding a set X if and only if there is no vector in X which dominates a, that is, there is no vector a′ ∈ X such that a′ ≺ a. All non-dominated vectors constitute admissible optimal solutions known in multi-objective optimization terminology as Pareto optimal solutions. The set of objective vectors y = J (a) corresponding to the set of Pareto optimal solutions a is called Pareto optimal front. The characteristics of the Pareto solutions are that they cannot be improved in any objective without causing degradation in at least one other objective. It should be noted that the parameter identi cation problem is traditionally solved by constructing a single objective from the multiple objectives, as follows J (x) = n
wi Ji (x) (4) i=1 using some weighting factors wi , i = 1; : : : ; n. The weights depend on the adequacy of the model class and the accuracy with which the measured modal data are obtained. More uncertain modal data should be given smaller weights. The results of the identi cation depend on the weights used. However, the choice of weights is arbitrary since the modeling error and the uncertainty in the measured data are usually not known a priori. The single objective is computationally attractive since conventional minimization algorithms can be applied to solve the problem. Formulating the parameter identi cation problem as a multi-objective minimization problem, the need for using arbitrary weighting factors for weighting the relative importance of each objective is eliminated. An advantage of the multi-objective identi cation methodology is that all admissible solutions in the parameter space are obtained which constitute model tradeo s in tting the di erent modal properties. These solutions are considered optimal in the sense that the t in one modal group cannot be improved without deteriorating the t in another modal group. The optimal points along the Pareto trade-o front provide detailed information about the quality of t in the corresponding Pareto optimal models. The set of Pareto optimal solutions can be obtained using evolutionary algorithms well-suited to solve the multi-objective optimization problem. One such algorithm is the Strength Pareto Evolutionary Algorithm (SPEA) [16] described in Section 3. Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 STRUCTURAL MODEL IDENTIFICATION 669 2.2. Modal grouping and measures of t The grouping of the modal properties {!r (x); r (x); r = 1; : : : ; M } into n groups and the selection of the measures of t J1 (x); : : : ; Jn (x) are usually based on user preference. The number and type of modal properties involved in the i-th group as well as the particular form of Ji (x) may depend on the modal characteristics (mode type, modal frequencies and=or modeshapes), their expected uncertainties, and the signi cance of each modal property on the model identi cation. Among the various choices available, the following are considered for illustration purposes. A group may be selected to contain the modal frequency and all modeshape components at the measured DOFs for a particular observed mode. In this case the number of groups equals the number of observed modes (n = M ). The i-th measure of t Ji (x) accounts for the mismatch between the measured and the model predicted frequencies and modeshape components for the i-th measured mode. This grouping scheme is appropriate when the objective of the identi cation is to estimate all optimal models that trade-o the t between di erent modes. Speci cally, Ji (x) can be given in the form Ji (x) = i (x) − ai ˆ i 2 !i (x) − ! ˆ i 2 ; + ! ˆ i 2 ˆ i 2 i = 1; : : : ; M (5) where z 2 = z T z is the usual Euclidian norm, and ai (x) = Ti (x)ˆ i = ˆ i 2 is a normalization constant that accounts for the di erent scaling between the measured and the predicted modeshape for given x. The set of objectives de ned in Equation (5) is referred to as set A. A second set of objectives, referred to as set B, can be de ned by grouping the modal properties into two groups according to their type as follows. The rst group contains all modal frequencies with the measure of t selected to represent the mismatch between the measured and the model predicted frequencies for all modes, while the second group contains the modeshape components for all modes with the measure of t selected to represent the mismatch between the measured and the model predicted modeshape components for all modes. Speci cally, the two measures of t can be given by J1 (x) = and J2 (x) = M 1
!r (x) − ! ˆ r 2 = !ˆ r 2 M r=1 M 1
 (x) − ar ˆ r 2 = ˆ r 2 M r=1 r (6) (7) This selection allows one to estimate all Pareto optimal models that trade-o the overall t in modal frequencies with the overall t in the modeshapes. 3. EVOLUTIONARY ALGORITHM FOR MULTI-OBJECTIVE OPTIMIZATION Evolutionary algorithms are well-suited for performing the multi-objective optimization. They process a set of promising solutions simultaneously and therefore are capable of capturing several points along the Pareto front. These algorithms are based on an arbitrarily initialized Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 670 Y. HARALAMPIDIS, C. PAPADIMITRIOU AND M. PAVLIDOU population of search points in the parameter space, which by means of selection, mutation, and recombination evolves towards better and better regions in the search space. In this work, a recently proposed Strength Pareto Evolutionary Algorithm (SPEA) [16] based on evolution strategies is used for solving the multi-objective minimization problem. For completeness, the main concepts of evolution strategies for single-objective optimization are rst introduced and then a brief outline of the SPEA method with its key features is presented. 3.1. Evolutionary strategies (ES) Details on theoretical developments of evolution strategies (ES) can be found in Beyer [13]. The ES operates with populations of size (; ) or ( + ), where  is the size of the parent population and ¿¿1 is the size of the o spring population. The members in a population are called individuals. Each individual is characterized by a vector of object variables x = (x1 ; x2 ; : : : ; xm )T , by an additional vector of strategy parameters  = (1 ; 2 ; : : : ; m )T , and by its tness value J (x), where J (x) is the single objective function to be optimized. The population is evolved by the successive application of genetic operators on the population of individuals. These operators are selection, recombination and mutation. Each cycle of evolution is called generation. Each generation consists of a parent population of size . Using the genetic operators of recombination and mutation,  o springs are generated from the parent population, forming the population of descendants. Selection is then used to form the parent population in the next generation from the  o springs of the current population. The genetic operations of recombination and mutation are discussed next. Recombination is applied not only to the object variables x = (x1 ; x2 ; : : : ; xm )T but also the strategy parameters  = (1 ; 2 ; : : : ; m )T . Two kinds of recombination operators are used. The discrete recombination performed on the object variables and the intermediate recombination performed on the strategy parameters. In discrete recombination the o spring solution inherits its components from two parents chosen randomly with equal probability from the pool of parent population. Speci cally, in generation g, the components xig of the descendant xg are generated from the components x1;g i and x2;g i of the two randomly selected parents x1g and x2g as follows: xig = x1;g i or x2;g i with equal probability 1=2 (8) Discrete recombination generates corners of the hypercube de ned by the parents. In intermediate recombination the o spring solution inherits components, which are a weighted average of the components from two parents x1g and x2g selected randomly from the pool of  parents with equal probability 1=, according to the rule: xig = i x1;g i + (1 − i )x2;g i (9) where i is a scaling factor chosen uniformly at random over an interval [0; 1]. Intermediate recombination is capable of producing any point within a hypercube de ned by the parents. Mutations are typically implemented by adding to the parental recombination result xg a Gaussian random vector z with independent components and distribution N (0; ) as follows: x̃g = xg + z (10) The standard deviations (i.e. the square roots of the diagonal elements i2 of ) are called step sizes of the mutation and constitute the strategy parameters of the individual. A key feature Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 STRUCTURAL MODEL IDENTIFICATION 671 of ES is that they self-adapt these step sizes too. Speci cally, the step size is modi ed by taking the product of itself with a log-normally distributed random number as follows: i ← i exp(z0 + zi ) (11)  √ √ where z0 = N (0; 1)= 2m and zi = N (0; 1)= 2 m. The variable z0 is computed once and is the same for all i , while zi is computed separately for each i . As a result, general and speci c scalings can be learned at the same time. The exponential function guarantees positive values for the mutation steps. In each generation every object variable and strategy parameter of an individual is mutated with the following order: (1) mutate strategy parameter, (2) mutate object variable. 3.2. The Strength Pareto Evolutionary Algorithm (SPEA) The SPEA uses a number of features speci c to multi-objective optimization algorithms for nding the multiple Pareto optimal solutions in parallel. Speci cally, it stores in an external set the non-dominated solutions found in each generation. It uses the Pareto dominance concept in order to assign tness values to individuals since the objective vector y = J (x) de ned in Equation (2) does not qualify as a scalar tness function. The tness of an individual is determined only from the solutions stored in the external non-dominated set. The solutions in the external set participate in the selection. It accomplishes tness assignment and selection that guides the search towards the Pareto optima set. It maintains diversity in the population so that a well-distributed, wide-spread trade-o front is reached, preventing premature convergence to a part of the Pareto front. Finally, it performs clustering to reduce the number of non-dominated solutions. The algorithm as proposed by Zitzler and Thiele [16] is the following: Step Step Step Step 1: 2: 3: 4: Step 5: Step 6: Step 7: Step 8: Generate an initial population P and create the empty external set P ′ . Copy non-dominated members of P to P ′ . Remove solutions within P ′ which are covered by any other member of P ′ . If the number of externally stored non-dominated solutions exceeds a given maximum N ′ , prune P ′ by means of clustering. Calculate the tness of each individual in P as well in P ′ . Select individuals from P + P ′ (multiset union), until the mating pool is lled. Apply recombination and mutation to members of the mating pool in order to create a new population P. If the maximum number of generations is reached, then stop, else go to Step 2. In Step 1 the initial population consists of  parents chosen randomly from the feasible region in the parameter space. In Step 5, all individuals in P and P ′ are assigned a scalar tness value. This is accomplished in the following two-stage process. First, all members of the non-dominated set P ′ are ranked. Afterwards, the individuals in the population P are assigned their tness value. Stage 1: Each solution i ∈ P ′ is assigned a real value si ∈ [0; 1) called strength si which is proportional to the number of population members j ∈ P for which i  j. Let n denote the number of individuals in P that are covered by i and assume N to be the size of P. Then si is de ned as si = N n+1 . The tness fi of i is equal to its strength fi = si . Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 672 Y. HARALAMPIDIS, C. PAPADIMITRIOU AND M. PAVLIDOU Stage 2: The tness of an individual j ∈ P is calculated by summing the strengths of all external non-dominated solutions i ∈ P ′ that cover j. Add one to this sum to guarantee that members of P ′ always have better tness than members of P: fj = 1 +
(12) si i; ij The stronger a non-dominated solution, the less t are the covered individuals. The above ranking method gives a preference to individuals near the Pareto optimal front and distributes them along the trade-o surface. In Step 6 the selection of individuals is done using the tournament selection procedure. Speci cally, two individuals are selected randomly from P + P ′ and the one with better tness value is copied into the mating pool. This type of selection is elitism and ensures that good individuals do not get lost. In Step 7 recombination and mutation is taking place. Discrete recombination is applied to the control variables and intermediate recombination to the step sizes according to the formulas (8) and (9). Mutation is done using the self-adaptive mutation method where the step size in each element in the mating pool is adapted according to (11). In Step 4, the number of externally stored non-dominated solutions is limited to some number N ′ . This is necessary because otherwise P ′ would grow to in nity since there is always an in nite number of points along the Pareto-front. Moreover, one wants to be able to control the number of proposed possible solutions because, from a structural identi cation point of view, a few points along the front are often enough to provide a complete description of the wider variety of non-dominated structural models. The reason for introducing clustering is the distribution of solutions along the Pareto-front. In order to explore as much of the front as possible, the non-dominated members of P ′ should be equally distributed along the Paretofront. Without clustering, the tness assignment method would probably be biased towards a certain region of the search space, leading to an unbalanced distribution of the solutions. In this work the Single Linkage method [18] has been chosen for clustering. 4. APPLICATION The multi-objective parameter identi cation framework is demonstrated using experimental modal data from a scaled two-story aluminum building model built and tested in the Structures Laboratory of the Hong Kong University of Science and Technology. The side and front views of the shear building model are given in Figure 1. The oors of the shear building are made of aluminum plates with dimensions 300×300×8 mm. Each oor is supported by four columns and each column is made up of two identical aluminum bars with dimensions 38×4:5×274 mm. The columns are connected to the oor and base plates through angles with the help of bolts and nuts. The structure was placed on a shake table and subjected to a base excitation. Acceleration response time histories at the two oors and the base along the weakest direction were obtained and the rst two modes were identi ed using modal analysis software. The rst two modal frequencies and modeshape components at the two oors are given in Table I. Details of the model, the experiments and the modal analysis performed to identify the modal frequencies and modeshapes can be found in Reference [19]. Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 673 STRUCTURAL MODEL IDENTIFICATION Figure 1. Side and front views of the shear building model (dimensions in mm). Table I. Experimental modal frequencies and modeshapes of the two-story aluminum building model. Modal frequencies (Hz) Modeshapes 1st Floor 2nd Floor Mode 1 Mode 2 17.2 0.48 1.00 50.4 −1.23 1.00 4.1. Identi cation of a two-parameter model The multi-objective parameter identi cation methodology is applied to identify the sti ness properties of the two oors. For this, the structure is modeled by a two-DOF linear lumped mass shear building model schematically shown in Figure 2. In the modeling, the masses are treated as deterministic, while the model parameters are chosen to be the two interstory sti nesses of the two-story building. According to Reference [19], the model masses were estimated from the structural drawings to be m1 = 3:9562 kg and m2 = 4:4482 kg. The a priori best estimate of the interstory sti ness calculated from the structural drawings is the same for both stories and equal to k0 = 2:3694×105 N=m. The following parameterization of the twoDOF model shown in Figure 2 is used: ki = xi k0 ; i = 1; 2. The purpose of the identi cation is to update the values of the sti ness parameters x1 and x2 using the measured modal data reported in Table I. The set A of objectives de ned in Equation (5) is rst used to identify the Pareto optimal values of the parameter set x = (x1 ; x2 ). In order to identify and describe in detail the whole Pareto front, a high number N ′ = 500 of non-dominated solutions in the set P ′ is selected. The Pareto front, giving the Pareto solutions in the objective space, is shown in Figure 3 for Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 674 Y. HARALAMPIDIS, C. PAPADIMITRIOU AND M. PAVLIDOU m2 k2 m1 k1 Figure 2. Two DOF lumped mass model. Ngen = 10, 100, 1000 and 10000 number of generations. The corresponding Pareto optimal solutions in the parameter space are shown in Figure 4. The search for the Pareto optimal solutions was limited to the region de ned by the inequality constraints Ji (x)60:05 for i = 1; 2. Pareto solutions outside this range are considered unacceptable in structural identi cation due to the very high errors between measured and model predicted modal data involved in the measures of t. The large number of 10000 generations is used to approximate as close as possible the Pareto front and the Pareto optimal solutions. The results for 10, 100 and 1000 generations demonstrate the rate of convergence of the solutions provided by the SPEA algorithm. The results in Figure 3 suggest that as the number of generations increases the computed Pareto front converges closer to the Pareto front computed for 10000 generations. It should be noted that as the number of generations increases, the number of non-dominated solutions stored in the external set P ′ also increases. Figure 3 suggests that the allowable number of N ′ = 500 non-dominated solutions are obtained only when the number of generations is suciently high. For the small number of 10 and 100 generations, the number of non-dominated solutions in the set P ′ is smaller than N ′ = 500. In such cases clustering is not activated and as a result a non-uniform distribution of the non-dominated points along the Pareto front is obtained. The results in Figure 4 suggest that the Pareto optimal solutions are concentrated in a narrow sub-region in the parameter space that extends along a one-dimensional manifold. The size of the sub-region perpendicular to the manifold decreases as the number of generations increases suggesting that the region approaches a one-dimensional manifold in the two-dimensional parameter space. These solutions correspond to interstory sti ness values which are signi cantly lower than the nominal sti ness values x1 = x2 = 1 estimated by structural drawings. This discrepancy is mainly due to modeling errors [19]. Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 675 STRUCTURAL MODEL IDENTIFICATION 0.03 0.03 0.025 0.025 Ngen = 10 0.015 0.015 0.01 0.01 0.005 0.005 0 0 0.01 Ngen = 100 0.02 J2 J2 0.02 0.02 0.03 0.04 0 0.05 0 0.01 0.02 J1 0.03 0.05 0.03 0.04 0.05 0.025 Ngen = 1000 0.015 0.015 0.01 0.01 0.005 0.005 0 0.01 Ngen = 10000 0.02 J2 0.02 J2 0.04 0.03 0.025 0 0.03 J1 0.02 0.03 J1 0.04 0.05 0 0 0.01 0.02 J1 Figure 3. Pareto optimal solutions in objective space. In order to evaluate the e ectiveness of the proposed SPEA algorithm in adequately describing the Pareto front with a small number of points along it, the Pareto front and the corresponding Pareto solutions are also obtained for N ′ = 20 non-dominated solutions in the set P ′ using 200 generations. The results are shown in Figure 5 and are compared to the results obtained for the case of 1000 generations and N ′ = 500. The Pareto solutions for N ′ = 20 are numbered so that their correspondence in the objective and the parameter spaces becomes clear. It is seen that the N ′ = 20 points de ning the Pareto front are uniformly distributed along the more detailed Pareto front represented by a high number of N ′ = 500 non-dominated solutions. Although the convergence in the parameter space has not been achieved with high accuracy, from the practical point of view the N ′ = 20 Pareto optimal solutions seem to give a very good representation of the Pareto front and the corresponding Pareto optimal models. It is observed in Figure 5 that a wide variety of Pareto optimal solutions are obtained. These solutions are superior to all other solutions when both objectives are considered. Comparing the Pareto optimal solutions, there is no Pareto optimal solution that improves the t in both measures simultaneously. Thus, all the Pareto solutions correspond to acceptable compromise structural models trading-o the t in the modal properties of the two modes. The non-zero size Pareto front and the non-zero distance of the Pareto front from the origin are due to modeling and measurement errors. Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 676 Y. HARALAMPIDIS, C. PAPADIMITRIOU AND M. PAVLIDOU 0.9 0.8 Ngen = 100 0.8 0.7 Ngen = 10 x2 x2 0.7 0.6 0.6 0.5 0.5 0.4 0.3 0.4 0.5 0.6 0.4 0.7 0.3 0.4 x1 0.8 Ngen= 1000 0.7 0.6 0.5 Ngen = 10000 0.7 x2 x2 0.6 0.8 0.7 0.4 0.5 x1 0.6 0.5 0.3 0.4 0.5 x1 0.6 0.7 0.4 0.3 0.4 0.5 0.6 0.7 x1 Figure 4. Pareto optimal solutions in parameter space. Modeling errors, in particular, relate directly to the adequacy of the selected model class and are expected to be more signi cant than measurement errors. Therefore, the size of the Pareto front, the distance of the Pareto front from the origin, and the distance between the optimal structural models in the parameter space are mainly a ected by the adequacy of the selected model class. For model classes that can exactly represent structural behavior and for error free measured data, the multi-objective identi cation methodology will in fact yield a zero-size Pareto front that consists of an isolated point located at the origin in the objective space. This is due to the fact that a solution Ji (x) = 0 ∀i = 1; : : : ; n (n = 2) exists which dominates all other solutions in the objective space. The corresponding optimal solutions in the parameter space may consist of either isolated points or an in nite number of points along a manifold, depending on whether the parameterized class of structural models is identi able or unidenti able [20–22], respectively, for the available measured data. However, since it is unrealistic to expect that a model class can represent the underlying physical phenomena exactly, model errors will always be present and the proposed multiobjective structural identi cation will always yield a set of Pareto optimal solutions. In the present application, a simple parametric model class has been selected for updating in order to illustrate the proposed multi-objective identi cation framework. More adequate model classes such as detailed parametric nite element model classes will bring the Pareto front closer Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 677 STRUCTURAL MODEL IDENTIFICATION 0.03 Ngen = 1000, N' = 500 Ngen = 200, N' = 20 Single Objective 0.025 2 1 0.02 11 J2 10 0.015 13 0.01 18 17 0.005 19 0 0 0.005 20 5 0.01 0.015 3 4 0.02 16 0.025 8 9 0.03 15 6 0.035 7 12 0.04 14 0.045 0.05 J1 (a) 0.8 141276 9 15 Ngen = 1000, N' = 500 Ngen = 200, N' = 20 Single Objective 8 16 4 3 0.75 5 20 19 0.7 0.65 x2 17 0.6 18 0.55 13 0.5 11 1 2 0.45 0.4 10 0.3 0.35 (b) 0.4 0.45 0.5 0.55 0.6 0.65 0.7 x1 Figure 5. Pareto optimal solutions for set A in the region Ji (x)60:05: (a) objective space; and (b) parameter space. to the origin in the objective space, reducing the size of the Pareto front, and reducing the distance between the optimal models in the parameter space, thus improving the predictive accuracy of the model class. Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 678 Y. HARALAMPIDIS, C. PAPADIMITRIOU AND M. PAVLIDOU Although a Pareto optimal structural model should always be a better compromise model than any model it dominates, not all Pareto optimal structural models may constitute acceptable compromise models. Information about the size of modeling and measurement errors is important for selecting a subset of optimal structural models along the Pareto front, eliminating unacceptable Pareto optimal solutions. For example, given that the t in both measures is expected to be less than 0.01, one could disregard all Pareto solutions that fall outside the region of the objective space de ned by Ji (x)60:01 for i = 1; 2. Figure 6 shows the Pareto front and optimal solutions within the smaller region. There is still a wide variety of solutions with maximum distance between them in the parameter space equal to dmax = maxi; j∈P′ xi − xj  = 0:252. There is a signi cant amount of information that can be extracted from the Pareto front and Pareto optimal solutions. For the sake of demonstration, let us assume that Pareto solutions are only acceptable if Ji (x)6J0 ∀i = 1; : : : ; n, i.e. the corresponding measures of t are both below a pre-speci ed level J0 . The size of the Pareto front and the maximum distance dmax (J0 ) between the acceptable solutions in the parameter space depends on the level J0 selected in applications to re ect the size of modeling and measurement errors. The size of the region in the parameter space that provides alternative optimal solutions decreases as J0 decreases. Moreover, there is a bound on J0 = 3:4×10−3 below which the values of the two objectives cannot be reduced simultaneously. Such a bound is due to modeling error arising from the assumptions adopted in structural modeling, as well as measurement error present in the data. Speci cally, the Pareto solution x = [0:546 0:648]T maintains the t for both measures to values close to J1 (x) = J2 (x) = 3:4×10−3 which is the minimal value that can be achieved for both measures simultaneously. It should be noted that the aforementioned solution differs by as much as d = x − xso  = 0:0508 from the solution xso = [0:528 0:695]T obtained by minimizing the single objective (4) with weighting factors w1 = w2 = 0:5. This di erence is 20% of the maximum distance dmax between the Pareto optimal solutions in the parameter space. Returning to the Pareto optimal solutions shown in Figure 5, it is observed that there is no other Pareto optimal solution that can improve the t in both objectives to values less than 3:4×10−3 . However, there is a Pareto optimal model with x = [0:652 0:420]T that improves the t for the properties of the rst mode to values as low as J1 (x) = 3:6×10−7 at the expense of deteriorating the t in the properties of the second mode to a value as high as J2 (x) = 0:0258. This solution could be considered an acceptable one only if the objective of the identi cation is to construct a model that provides the best t in the properties of the rst mode, ignoring higher errors in the properties of the second mode. For example, this would be the case when there is information to support that measurement error is negligible in the modal properties involved in the rst mode and that the parameterized structural model can accurately predict the properties of the rst mode. Similarly, the Pareto optimal solution x = [0:297 0:798]T shown in Figure 5 ts very well the modal properties of the second mode to values as low as J2 (x) = 3×10−4 at the expense of deteriorating signi cantly the t in the modal properties of the rst mode to values as high as J1 (x) = 0:0499. Similar interpretation can be given to this solution which can be accepted or rejected depending on the information available for the size of modeling and measurement errors. Next, the Pareto optimal set is obtained for the two-parameter model using the set B of objectives de ned in Equations (6) and (7). The rst objective measures the overall mismatch of the modal frequencies, while the second one measures the overall mismatch of the Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 679 STRUCTURAL MODEL IDENTIFICATION -3 10 x 10 Ngen = 1000, N' = 500 Ngen = 200, N' = 20 Single Objective 16 9 15 7 8 6 10 7 6 J2 19 20 5 18 4 17 3 9 8 2 4 5 11 14 12 13 2 1 2 3 4 5 6 7 1 8 3 9 J1 (a) 1 2 13 12 14 5 11 -3 Ngen = 1000, N' = 500 Ngen = 200, N' = 20 Single Objective 3 0.75 10 x 10 4 0.7 8 x2 9 0.65 17 18 20 19 0.6 6 10 7 16 15 0.55 0.44 0.46 0.48 (b) 0.5 0.52 x1 0.54 0.56 0.58 0.6 Figure 6. Pareto optimal solutions for set A in the region Ji (x)60:01: (a) objective space; and (b) parameter space. Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 680 Y. HARALAMPIDIS, C. PAPADIMITRIOU AND M. PAVLIDOU modeshape components. The Pareto front is shown in Figure 7(a) for 1000 generations and N ′ = 500 as well as for 200 generations and N ′ = 20. The corresponding Pareto solutions in the parameter space are shown in Figure 7(b). It is seen that the SPEA algorithm e ectively distributes the solutions along the Pareto front so that the whole front is adequately described with a few points. In contrast to the set A, the distance dmax (J0 ) and the size of the Pareto front for the set B of objectives do not grow as J0 increases. The Pareto optimal solutions shown in Figure 7(b) are concentrated in a relatively narrow sub-region in the parameter space. The values of the rst objective, measuring the t in the modal frequencies of both modes along the Pareto front, are in the range J1 (x) ∈ [1:5×10−7 ; 0:002], while the values of the second objective, measuring the t in modeshape components of both modes, are in the range J2 (x) ∈ [0:0026; 0:0032] which is a much narrower range than that for the modal frequencies. The best t for the modeshape components cannot be less than J2 (x) = 2:6×10−3 . The t in the modal frequencies can be substantially improved to as low as J1 (x) = 1:5×10−7 at the expense of deteriorating by approximately 23% the t in the modeshape components to the value J2 (x) = 3:20×10−3 . This could suggest that if the 23% di erence in the t in the modeshape components is not of concern in the identi cation, then among all Pareto optimal models, the model x = [0:511 0:718]T corresponding to J1 (x) = 1:5×10−7 and J2 (x) = 3:20×10−3 is most representative of the structure. The distance of this solution from the one obtained by minimizing the single objective (4) with weighting factors w1 = w2 = 0:5 is d = 0:0286 which is equal to 18% of the maximum distance dmax = maxi; j∈P′ xi − xj  = 0:1571 between the Pareto optimal solutions in the parameter space. Finally, comparing the results in Figures 6 and 7 for the di erent sets A and B of objectives, it is clear that the Pareto optimal structural models also depend on the selection of modal properties involved in each modal group. 4.2. Structural reliability predictions using Pareto optimal models The purpose of identi cation is to construct faithful structural models, within a selected model class, that can be used for making improved structural reliability predictions consistent with the measured data. The Pareto optimal set provides complete information about all optimal structural models consistent with the data that can be traded-o based on the accuracy they provide in the modal properties. The alternative models along the Pareto front provide di erent predictions of structural reliability which are all acceptable based on the measured data. Next the variability in the reliability predictions from all corresponding Pareto optimal models is explored. The reliability of the structure given a model along the Pareto front can be computed using available probabilistic structural analysis tools. For simplicity, the reliability of the structure subjected to a base excitation that can be adequately modeled by white noise is considered. The system is considered to have failed under the stochastic base excitation if a response quantity y(t; x) of the structure exceeds a threshold level b over a duration T . An estimate of the failure probability of the structure is obtained using well-known approximate random vibration results for linear structures. Speci cally, for given model parameters x, the probability of failure is approximated by [23]: P(F |x) = 1 − exp[1 − 2(x)T ] Copyright ? 2005 John Wiley & Sons, Ltd. (13) Earthquake Engng Struct. Dyn. 2005; 34:665–685 681 STRUCTURAL MODEL IDENTIFICATION -3 3.3 x 10 Ngen = 1000, N' = 500 Ngen = 200, N' = 20 Single Objective 3.2 14 3.1 10 11 3 J2 6 5 2.9 12 19 2.8 20 17 2.7 16 18 15 2.6 0 0.5 13 8 7 1 (a) 2 1 3 4 9 1.5 2 2.5 J1 x 10 -3 0.72 Ngen = 1000, N' = 500 Ngen = 200, N' = 20 Single Objective 14 0.7 10 11 0.68 6 5 12 20 0.66 19 x2 17 0.64 16 18 15 0.62 13 8 7 2 0.6 1 3 4 9 0.58 0.56 0.5 0.52 (b) 0.54 0.56 x1 0.58 0.6 0.62 Figure 7. Pareto optimal solutions for set B: (a) objective space; and (b) parameter space. Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 682 Y. HARALAMPIDIS, C. PAPADIMITRIOU AND M. PAVLIDOU 0 10 Ngen = 1000, N' = 500 Ngen = 200, N' = 20 -1 Probability of Failure 10 -2 10 -3 10 -4 10 0.01 0.011 0.012 0.013 0.014 0.015 0.016 1/2 Exceedance Level b/S Figure 8. Reliability predictions based on all Pareto optimal structural models. where (x) is the rate of outcrossing the level b, given by   b2 1 ẏ (x) exp − 2 (x) = 2 y (x) 2y (x) (14) The quantities y (x) and ẏ (x) are the standard deviations of the response y(t; x) and its derivative ẏ(t; x), respectively. These standard deviations of the response are readily obtained for a linear system subjected to white excitation by using the Liapunov equation √ for the covariance matrix [23]. The reliability depends on the normalized threshold level b= S, where S is the constant spectral density of the input white noise. More realistic descriptions of the base motion could readily be incorporated in the aforementioned formulation. It should be noted that accurate estimates of the probability of failure can also be obtained using ecient Monte Carlo techniques [24]. However, this is outside the scope of the present example. In the numerical illustration, the structure is considered to have failed when the displacement response of the top oor exceeds a threshold level b. The probability of failure predicted by all Pareto optimal models √ for the set A of objectives is shown in Figure 8 for di erent normalized threshold levels b= S. It can be seen that these predictions can vary considerably. Thus, the selection of a model among all Pareto optimal models has an impact on the reliability. Di erent models along the Pareto front result in di erent predictions. In general, the variability in the predictions depends on the adequacy of the model class selected for identi cation. The results for the probability of failure in Figure 8, are computed and compared for a small number of N ′ = 20 non-dominated solutions and for a suciently large number of N ′ = 500 non-dominated solutions. It is evident that the failure probability bounds predicted by the Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 STRUCTURAL MODEL IDENTIFICATION 683 small number of N ′ = 20 non-dominated solutions match closely the more accurate failure probability bounds predicted by a very high number of N ′ = 500 non-dominated solutions. This means that from the prediction point of view, representation of the Pareto front by a small number of points, uniformly distributed along the front, is adequate. 5. EXTENSION TO NON-LINEAR STRUCTURES The multi-objective parameter identi cation framework can readily be extended to identify the free parameter set x of a class of non-linear models of structures such as, for example, hysteretic models with various degrees of modeling sophistication (e.g. elastoplastic, bilinear, degrading, etc). In this case modal analysis is no longer applicable and the identi cation must be based on the experimentally obtained response time histories (e.g. accelerations, velocities, displacements), denoted by ŷj (k
t), where the index j = {1; : : : ; M } refers to measurements at M DOFs, k = 1; : : : ; N is the time index and
t is the sampling period. For this, the measured response time histories are
grouped into n groups with the i-th group gi containing ni response time histories such that ni=1 ni = M . For the i-th group, a norm Ji (x) = M

j∈gi k=1 [yj (k
t; x) − ŷj (k
t)]2 (15) is introduced to measure the mismatch between the measured response time histories involved in the group and the corresponding response time histories yj (k
t; x) predicted from the nonlinear model class for a particular value of the parameter set x. The Pareto optimal solutions are then obtained by solving the multi-objective optimization problem stated in Equation (2). The estimation of yj (k
t; x) involved in Equation (15) requires the solution of the non-linear set of di erential equations governing the response of the non-linear structures. The selection of the response time histories that are involved in a group is based on user preference and may depend on the number, location and type (e.g. displacement, velocity or acceleration) of measurements, as well as the information contained in the measurements about the free model parameters to be identi ed. For example, responses that contain signi cant information about the local non-linear properties of the structure could be included in the same group. In this way the t of the responses in this group can be traded-o with the t in another group containing information for another set of local or global properties of the structure. Once the Pareto optimal non-linear models consistent with the measured data have been obtained, one could use these models to estimate the variability in the failure probability against various limit states which are indicative of damage, speci ed in terms of interstory drift for buildings, ultimate displacement, ductility demands, etc. For stochastic dynamic excitation, the estimation of structural reliability for each model along the Pareto front requires the integration into the formulation of ecient Monte Carlo techniques [25, 26]. 6. CONCLUSIONS A novel multi-objective framework for structural model identi cation based on modal data has been presented. Multiple objectives related to the t between measured and model predicted Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 684 Y. HARALAMPIDIS, C. PAPADIMITRIOU AND M. PAVLIDOU modal properties were simultaneously minimized, thus eliminating the need for arbitrarily weighting the relative importance of each objective in a single measure of t. In contrast to the conventional weighted least-squares t between measured and model predicted modal data, the proposed methodology provides multiple alternative optimal structural models consistent with the data in the sense that the t each model provides in a group of modal properties cannot be improved without deteriorating the t in at least one other group of modal properties. The multi-objective optimization problem can be solved using well-developed evolutionary algorithms. Speci cally, the Strength Pareto Evolutionary Algorithm was e ectively used for replacing the Pareto front by a relatively small number of representative optimal structural models, uniformly distributed along the Pareto front. These multiple Pareto optimal structural models are due to modeling and measurement errors. Information about the size of errors in the measured modal properties as well as information about the accuracy with which the selected class of structural models is expected to t the modal frequencies and modeshape components used in the identi cation are important for selecting the acceptable optimal structural models along the Pareto front, eliminating unacceptable Pareto optimal solutions. For the application considered, a wide variety of Pareto optimal structural models, uniformly distributed along the Pareto front, were obtained. The Pareto optimal solutions for the two-parameter model were concentrated in a narrow sub-region, extended along a certain direction in the parameter space. The variability of Pareto optimal models in uences the predictions of structural reliability. Using the whole set of identi ed structural models, the failure probability of structures to future uncertain loads was predicted. It was demonstrated that such predictions from the Pareto optimal models may vary considerably. The proposed methodology can provide information about the failure probability bounds that are useful for making more informed structural safety assessments consistent with measured data. Finally, the multi-objective identi cation framework is applicable to non-linear model identi cation based on response time histories. ACKNOWLEDGEMENTS This research was partially funded by the Greek Secretariat of Research and Technology within the EPAN program framework under grant DP15 and the Greek Earthquake Planning and Protection Organization. This support is gratefully acknowledged. The authors would like to thank Dr Paul H. F. Lam, Assistant Professor in the Department of Building and Construction, The City University of Hong Kong, for providing details on the experimental model and the measurements used in this work for illustrating the proposed multi-objective identi cation framework. REFERENCES 1. Mottershead JE, Friswell MI. Model updating in structural dynamics: A survey. Journal of Sound and Vibration 1993; 167(2):347–375. 2. Doebling SW, Farrar CR, Prime MB, Shevitz DW. Damage identi cation and health monitoring of structural and mechanical systems from changes in their vibration characteristics. A literature review. Los Alamos National Laboratory, Report No. LA-13070-MS, 1996. 3. Vanik MW, Beck JL, Au SK. Bayesian probabilistic approach to structural health monitoring. Journal of Engineering Mechanics (ASCE) 2000; 126(7):738 –745. 4. Bohle K, Fritzen CP. Results obtained by minimizing natural frequencies and MAC-value errors of a plate model. Mechanical Systems and Signal Processing 2003; 17(1):55 – 64. 5. Hjelmstad KD, Shin S. Crack identi cation in a cantilever beam from modal response. Journal of Sound and Vibration 1996; 198(5):527–545. Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685 STRUCTURAL MODEL IDENTIFICATION 685 6. Friswell MI, Mottershead JE. Finite Element Model Updating in Structural Dynamics. Kluwer: Dordrecht, 1995. 7. Farhat C, Hemez PM. Updating nite element dynamics models using an element-by-element sensitivity methodology. AIAA Journal 1993; 31(9):1702–1711. 8. Alvin KF. Finite element model update via Bayesian estimation and minimization of dynamic residuals. AIAA Journal 1997; 35(5):879 – 886. 9. Ewins DJ. Modal Testing—Theory, Practice and Applications, (2nd edn.), Research Studies Press: U.K., 2000. 10. McConnnel KG. Vibration Testing: Theory and Practice. Wiley: New York, 1995. 11. Katafygiotis LS, Yuen KV. Bayesian spectral density approach for modal updating using ambient data. Earthquake Engineering and Structural Dynamics 2001; 30:1103–1123. 12. Beck JL, May BS, Polidori DC. Determination of modal parameters from ambient vibration data for structural health monitoring. First World Conference on Structural Control, Los Angeles, 1994; 1395 –1402. 13. Beyer HG. The Theory of Evolution Strategies. Springer: Berlin, 2001. 14. Fonseca CM, Fleming PJ. An overview of evolutionary algorithms in multiobjective optimization. Evolutionary Computation 1995; 3(1):1–16. 15. Srinivas N, Deb K. Multiobjective optimization using nondominated sorting in genetic algorithms. Evolutionary Computation 1994; 2(3):221–248. 16. Zitzler E, Thiele L. Multiobjective evolutionary algorithms: A comparative case study and the strength Pareto approach. IEEE Transactions on Evolutionary Computation 1999; 3(4):257–271. 17. Mueller S, Walther J, Koumoutsakos P. Evolution strategies for lm cooling optimization. AIAA Journal 2001; 39(3):537–539. 18. Morse JN. Reducing the size of the nondominated set: Pruning by clustering. Computers and Operations Research 1980; 7(1–2):55 – 66. 19. Lam HF. Structural Model Updating and Health Monitoring in the Presence of Modeling Uncertainties. Ph.D. Thesis, Department of Civil Engineering, Hong Kong University of Science and Technology, Hong Kong, 1998. 20. Beck JL, Katafygiotis LS. Updating models and their uncertainties. I: Bayesian statistical framework. Journal of Engineering Mechanics (ASCE) 1998; 124(4):455 – 461. 21. Katafygiotis LS, Beck JL. Updating models and their uncertainties. II: Model identi ability. Journal of Engineering Mechanics (ASCE) 1998; 124(4):463– 467. 22. Katafygiotis LS, Lam HF. Tangential-projection algorithm for manifold representation in unidenti able model updating problems. Earthquake Engineering and Structural Dynamics 2002; 31:791– 812. 23. Lutes LD, Sarkani S. Stochastic Analysis of Structural and Mechanical Vibrations. Prentice Hall: New Jersey, 1997. 24. Au SK, Beck JL. First excursion probability for linear systems by very ecient importance sampling. Probabilistic Engineering Mechanics 2001; 16(3):193–207. 25. Pradlwalter HJ, Schueller GI. On advanced Monte Carlo simulation procedures in stochastic structural dynamics. International Journal of Non-linear Mechanics 1997; 32(4):735 –744. 26. Au SK, Beck JL. Estimation of small failure probabilities in high dimensions by subset simulation. Probabilistic Engineering Mechanics 2001; 16(4):263–277. Copyright ? 2005 John Wiley & Sons, Ltd. Earthquake Engng Struct. Dyn. 2005; 34:665–685