A comparative study of construction methods for seismic fragility curves using numerical simulations

– A seismic fragility curve of a structure or a mechanical system, presenting its failure probability versus seismic intensity, can be established by an engineering judgment approach, an empirical approach or a numerical approach. In the numerical approach, there exist three popular methods: the scaled seismic intensity method, the maximum likelihood estimation method, and the probabilistic seismic demand/capacity models. This paper is focused on a comparison of these numerical methods. Linear/non-linear oscillators and an eight-storey frame structure were used to derive fragility curves for diﬀerent tests. The results obtained show a discrepancy of fragility curves between the methods, and their quality is commented in comparison with results of the Monte-Carlo method with a high number of simulations. The maximum likelihood estimation method is in general the recommended method, due to its good accuracy in both numerical examples.


Introduction
A seismic fragility curve expresses the probability of failure or damage states of a structure or a mechanical system due to earthquakes as a function of a ground motion index; for instance, peak ground acceleration (PGA), peak ground velocity (PGV) and so on. . .It is one of three ingredients, comprising seismic hazard, fragility curves and dominant accident sequences, which lead to core damage in a plant in the context of a seismic probabilistic assessment (PRA) in nuclear engineering applications [1].Fragility curves are applied to different structural types e.g.buildings [2], bridges [3][4][5][6][7][8] etc.They are useful for the design of a new structure, and for an existing one they are helpful for seismic retrofitting decisions, disaster response planning and quick loss estimation [3,4].Fragility curves can be obtained using one of three approaches: (i) engineering judgment; (ii) empirical approach, (iii) numerical simulation.These approaches can also be used conjunctively.This paper focuses on the ways of constructing seismic fragility curves using numerical simulations, where there are three popular a Corresponding author: thien-phu.le@ifma.frmethods (i) scaled seismic intensity (SSI) [3,4,6,9] with least square regression technique that minimizes the sum of squared errors between the probabilities predicted by the fragility function and theirs estimates using the Monte-Carlo simulation and scaled ground motion histories, (ii) maximum likelihood estimation (MLE) [10][11][12][13] in which failure event is modelled by a Bernoulli random variable and the determination of its probability distribution leads to fragility curves and (iii) probabilistic seismic demand model/probabilistic seismic capacity model (PSDM/PSCM) [7,8,14,15] in which the failure condition is defined as a "demand -capacity (resistance)" problem.The knowledge of probabilistic models for the demand and the capacity allows deducing fragility curves as the failure is reached when the seismic demand exceeds the seismic capacity.The SSI, MLE and PSDM/PSCM methods are detailed in Section 2.2 with step-by-step practical procedures.Note that the MLE and PSDM/PSCM methods can be used for both scaled and non-scaled accelerograms.However, unmodified histories were chosen for these two methods in order to avoid eventual modification of probabilistic distribution for the seismic intensity at a site due to the scaling technique.Despite the existence of these methods, the choice of a one among them for the construction of a fragility curve is not simple because (i) comparative studies between the methods are not numerous and (ii) conclusions about the effectiveness of a method can vary from one structure to another, and are sometimes contradictory.

Article published by EDP Sciences
Baker [16] studied the distribution of structural response as a function of the ground motion index using several methods.He provides in this review an useful summary of the pros and cons of the different methods and shows that in general, one must make trade-offs between the accuracy of the estimates and the amount of data required for the estimation.
In recent comparative publications, Mandal et al. [17,18] showed that the MLE method provides the most accurate [17] and very good [18] fragility estimations when they compared the MLE method with the "conventional" method and the PSDM/PSCM method (called the "regression method" by the authors) via a test case concerning the primary containment dome of a typical Indian Pressurized Heavy Water Reactor (PHWR).Note that the conventional method is a simplified method proposed by Kennedy and Ravindra [19].In this method, the parameters of fragility curve namely median capacity and logarithmic-standard deviation are evaluated by means of safety factors with respect to design study and qualification testing of components and structures.The safety factors are also modelled by lognormal distributions.The mean value and logarithmic standard deviation of each safety factor are determined based on the design study and engineering judgement.The authors used PGA and pseudo-spectral acceleration at the fundamental mode (S a ) as the seismic intensity in their fragility curves.The comparison of the three methods referred to the discrete fragility data points considered as "actual" fragility values.These points should be questioned, because they were obtained from a multi-IDA (incremental dynamic analysis) study.Note that in the IDA, a suite of ground motion histories are used to analyse responses of a structure.The ground motion histories are scaled to predefined intensity measure levels, and nonlinear response time history analysis is performed for each scaled ground motion history.The fragility defined as conditional probability at a given intensity measure is then estimated by the Monte-Carlo method.Thus, the IDA still belongs to the scaling approach.
Gelh et al. [20] try to provide a guide for the number of necessary dynamic analyses when they use the PSDM/PSCM method (called "least-squares method" by the authors), the MLE method and the sum-of-squared errors (SSE) method (a modified formulation of the MLE method) to derive fragility curves.Based on a trial inversion procedure and verified by an example of a nonlinear single-degree-of-freedom, the authors noted that the least-squares method (i.e. the PSDM/PSCM method) converges much faster than the MLE and SSE methods.The seismic intensity measure used in the example is pseudo spectral acceleration at a period of 0.5 s.
The comparison between three methods is based on the exact fragility curve in the trial inversion procedure and the fragility curves derived from a very high number of simulations.The results obtained are, however, limited to structures whose characteristic responses can be expressed by a power law of the seismic index.The assumption is thus in favor of the the PSDM/PSCM method.
The objective of the paper is thus to perform a comparative study when the three prominent methods (SSI, MLE and PSDM/PSCM) are applied to two numerical examples: oscillators (single-degree-of-freedom) and an eight-storey frame structure (multi-degree-of-freedom) with linear and non-linear behaviors.The results of these methods are compared to those of the Monte-Carlo method with a high number of simulations.For sake of completeness, useful formulas and the SSI, MLE and PSDM/PSCM methods are first summarized in Section 2. Two comparison criteria are then proposed and developed in a practical procedure.Section 3 is devoted to the results of numerical examples involving three oscillators and an 8-storey frame structure.The results of this processing lead to some discussion in Section 4.

Methodology for the comparative study 2.1 Definitions
Let A be a chosen ground motion index.A can be peak ground acceleration (PGA), peak ground velocity (PGV), peak ground displacement (PGD), pseudo spectrum acceleration (PSa) and so on.The fragility curve F r (a) is the conditioned probability of failure or a damage state given that A = a: where the failure or the specific damage state is characterized when the structural response X exceeds a critical limit x 0 .Taking into account random properties of response and critical limit due to uncertainty in loadings, material properties, dimensions and so on, the failure or the specific damage state can be considered in a more general way as the exceeding by seismic demand D of seismic capacity C. The definition of the fragility curve thus becomes Once the fragility curve F r (a) is available, using the distribution of the ground motion index A of a site, it is possible to calculate the probability of failure or damage state of the site by integral where p A (a) is the probability distribution of index A for the considered site.In general, A is positive, which is why the integral bounds are limited from zero to infinity in Equation (3).A commonly accepted assumption for a fragility curve is that F r (.) is a cumulative probability function of the log-normal distribution: where Φ(.) is the standard Gaussian cumulative distribution function while A m and β are respectively median and logarithmic standard deviation.

Analysed methods
The log-normal assumption given in Equation ( 4) greatly simplifies the original problem defined in Equation (1) or Equation (2).The analysed methods now lead to the determination of two parameters, A m and β.In the following, the SSI, MLE and PSDM/PSCM methods are summarized in step-by-step procedures.The ground motion index A used in this study is peak ground acceleration (PGA); nevertheless it is easy to apply the approach to other seismic indexes.1. Scaled Seismic Intensity (SSI) method [3,4,6,9] The SSI method consists first of determining some discrete points of a fragility curve and then performing a least square regression using the log-normal assumption.Each discrete point is obtained by ground accelerations normalized to a common PGA value, and its corresponding probability of failure is estimated by the Monte-Carlo method.The step-by-step procedure of this method is summarized as follows: -  values of pf (γ j ) using the log-normal assumption in Equation ( 4) 2. Maximum Likelihood Estimation (MLE) method [10][11][12].In this method, the basic idea is to associate the failure or damage state event with a Bernoulli random variable Y .The use of MLE to identify the distribution of Y leads to the fragility curve.The step-by-step procedure is as follows: - From data and the results of numerical simulations, the parameters of the two models are first obtained.The fragility curve is then deduced from the definition in Equation (2) as: -Step 3: Perform a linear regression from N points (a i , δ i ) to have two coefficients c 1 and c 2 following the model : ln(δ) = ln(c 1 ) + c 2 ln(a).We then -Step 4: Identify two parameters S C and β C of the seismic capacity C following a log-normal distribution.The values of S C and β C can be obtained by experimental results and/or expert opinion [21].
Due to lack of the available data, seismic capacity C can be treated as deterministic [18].Note that for the case of the deterministic limit x 0 as in the definition in Equation ( 1

Comparison criteria
In order to compare the effectiveness of the proposed method applied to numerical examples, it is necessary to propose some comparison criteria.Using results of the Monte Carlo simulation method as reference, the first criterion is based on the distance between a fragility curve and the reference results.This criterion is called the mean square error (MSE) and defined as: where F X r (a) is the curve obtained from one of three methods: SSI, MLE ou PSDM/PSCM.The higher the MSE value of a fragility curve, the farther it is from the Monte Carlo results.
A second criterion concerns the use of a fragility curve for the estimation of failure probability using Equation (3).It is the Relative ERror (RER) between the failure probability estimated by an existing method p X f and that estimated by the Monte Carlo method p MCS where p X f is obtained from F X r (a) and the probability density function p A (a) of PGA: The step-by-step procedure of the Monte Carlo simulation method to obtain F MCS r (.) and p MCS f is the following: -

Practical procedure
The comparative study considered here consists of applying the same set of ground motions to various structures.Fragility curves and failure probabilities obtained by the SSI, MLE and PSDM/PSCM methods are then compared with results of Monte Carlo simulations following the two criteria presented above.The best method should give the lowest MSE and RER values.The practical step-by-step procedure is as follows: (1) Estimate the probability density function of the ground motion index p A (a) for the considered site.This is a crucial step and not an easy task.There are following possibilities: p A (a) is already available thank to existing seismic hazard analysis or experts judgment for example.p A (a) is not available and it should elaborate from ground motion histories.A large set of ground motion histories of the site is recommended for an accurate estimation of p A (a).
• If the distribution law of p A (a) is available, for instance log-normal distribution, the parameters of p A (a) can be easily determined by a parametric method for example: maximum likelihood-based method or momentsbased method.• If the distribution law of p A (a) is not available, p A (a) can be deduced by using a nonparametric method such as the kernel method from ground motion histories.(2) For a given structure and a given failure or damage state condition, establish seismic fragility curves F X r (.) and estimate failure probabilities p X f using the three existing methods and also the Monte Carlo simulation method.
(3) Calculate MSE and RER indicators following Equations ( 12) and ( 13) respectively.Between two methods, the one with lower MSE and RER values is of better quality.

Results of numerical examples
The practical procedure is now applied for two numerical examples: linear/non-linear oscillators and an eight-storey bi-linear frame.Time ground accelerations were generated using the stochastic method proposed by Boore's studies [22,23].Fragility curves of each structure derived from a same set of time accelerations the SSI, MLE and PSDM/PSCM methods were compared via the criteria given in Section 2.3.

Ground motion histories
The stochastic method proposed by Boore [22] assumes that ground motion is distributed with random phase over a time duration related to earthquake size and propagation distance.The ground motion is characterized by its spectrum -this is where the physics of earthquake process and wave propagation are contained, usually encapsulated and put into the form of simple equations.The total spectrum of the motion at a site S(M 0 , R, f) is considered as a combination of earthquake source (E), path (P ), site (G) and type of motion (I).

S(M
where M 0 is the seismic moment that is related to the seismic magnitude M by R is the distance from source to site and f is frequency.By separating the spectrum of ground motion into source, path, and site components, the models based on the stochastic method can be easily modified to account for specific situations or to account for improved information about particular aspects of the model.Given the spectrum motion at a site, Boore [22] suggest the simulation of ground motions by 6 steps: (i) white noise (Gaussian or uniform) is generated for a duration given by the duration of the motion; (ii) this noise is then windowed; (iii) the windowed noise is transformed into frequency domain; (iv) the spectrum is normalized by the square-root of the mean square amplitud spectrum; (v) the normalized spectrum is multiplied by the ground motion spectrum S(M 0 , R, f); (vi) the resulting spectrum is transformed back to time domain to obtain a time history of ground motions.This method for simulation of time series is implemented by Boore in Fortran [23].
In order to generate time acceleration series for the numerical examples, we implemented the stochastic method in Matlab software [24].Parameters of the ground spectrum were taken as suggested values in reference [22] and summarized as bellow: -The source E(M 0 , f): is based on the source spectral shape AS00 (Atkinson and Silva (2000) for California).-The path P (R, f ): accounts the effects of geometrical spreading, attenuation and the increase of duration with distance due to wave propagation and scattering.
where c Q = 3.5, Q(f ) = 180f 0.45 and Z(R) defined by: with R 0 = 1, R 1 = 40 and p 1 = 0.5.-The site G(f ): accounts the effects due to local site geology and is separated to amplification A(f ) and attenuation D(f ) where the amplification A(f ) is taken for generic rock in reference [25] and the attenuation function D(f ) is defined as the combination of two filters and D(f ) = exp(−πκ 0 f ) (22) where f max = 100 Hz and κ 0 = 0.03.-The motion type I(f ): for acceleration is defined as Figure 1 shows an acceleration time history of the model for magnitude M = 7 and R = 9 km.Using 2 × 10 5 accelerations with the same magnitude and distance, a probability density function of PGA can be obtained by the kernel method and is indicated in Figure 1.

Oscillators
Two oscillators (linear and Bouc-Wen non-linear) used in reference [9] and a Coulomb non-linear oscillator, were considered here.
Figure 3 shows an illustration of oscillator responses under seismic excitation, obtained with the Runge-Kutta algorithm using Matlab software [24].
The failure of an oscillator was verified by comparing displacement x(t) with a displacement limit x 0 i.e. in the definition of Equation ( 2 x 0 it was not reached.Three displacement limits were tested: x 0 = 7 cm, x 0 = 10 cm and x 0 = 13 cm.These limits in increasing order can be associated to three damage states: minor, moderate and severe in real applications. Fragility curves obtained by SSI, MLE and PSDM/PSCM methods with N = 5000 ground accelerations generated from Boore's model (M = 7, R = 9 km) are given in Figures 4-6 for respectively the linear, Bouc-Wen and Coulomb oscillators.
Note that the Monte Carlo Simulation method (MCS) was also performed with a greater number of ground accelerations: N MCS = 2 × 10 5 following the procedure given in Section 2.3.These accelerations were organized in N I = 11 PGA intervals a j : [a j − 0.1; a j + 0.1] m/s 2 .These intervals have the same width of 0.2 m.s −2 .This choice allows to have enough number of ground motion time histories in each interval (for the evaluation of failure probability by MCS); and the width of each interval is sufficiently narrow, thus its central value a j can be representative for the interval.Note that the number of time histories in 11 intervals is not necessary equal.The MCS results are represented in Figures 4-6 by plus sign markers (+) together its bounds of 95% confidence interval.It can be noted from these figures that the fragility curves obtained by the MLE and PSDM/PSCM methods are very similar and closer to MCS results than those obtained using the SSI method.When fragility curves of the   MLE method and the PSDM/PSCM method are different, for instance in the case of a Coulomb oscillator with x 0 = 13 cm, the PSDM/PSCM method is better since its fragility curves are closer to the plus sign markers than those of the MLE method.
Two comparison criteria, MSE and RER, were calculated.They are given in Table 1 and also presented in Figure 7.
It can be noted that the SSI method gives the highest MSE values for every failure threshold value x 0 .It is quite the same for RER values except for x 0 = 10 cm (second case) with linear and Bouc-Wen oscillators.The MSE and RER values of the MLE and PSDM/PSCM methods are of the same order.The PSDM/PSCM method seems to be slightly better in the last case of failure threshold x 0 = 13 cm.Note that the higher the value of x 0 , the lower the number of failure events.The fragility curves of the SSI method fitted well its discrete failure probabilities that are shown in Figures 4-6 by circle markers (o).However these points obtained using the scaling technique are far from the failure probabilities given by the MCS reference results with unmodified ground motion time histories.

8-storey frame structure
An 8-storey floor shear frame inspired by references [26,27] was tested.Lumped masses are given in Table 2.The non-linear hysteretic model with K y /K 0 = 0.1 in Figure 8 was considered for the shear stiffness of one storey K.The initial inter-storey stiffness was K 0 = 2.1 × 10 8 N.m −1 and the yielding inter-storey displacement was Δ y = 0.01 m.Rayleigh damping was considered as a linear combination of mass matrix M and stiffness matrix K as : C = 0.01M + 0.005K.The height of each storey was H = 3.0 m.
The same ground motions generated with the Boore model in the previous numerical example were used for the construction of fragility curves.Structural responses were simulated using ANSYS finite element software [28] and are illustrated in Figure 9.The responses were obtained using the "full method" option that allows time non-linear analysis based on the Newmark time integration method.Time step was 0.005 sec and time length was equal to seismic duration of 50 s.
The failure of the frame depends on the inter-storey drift i.e. the relative horizontal displacement of two adjacent floors.Under the excitation of a ground acceleration a i (t), the inter-storey drift of the kth storey is defined as

Failure or damage occurs if
where δ 0 is the seismic resistance in terms of the interstorey drift limit of the frame.Three cases were considered: δ 0 = H/300, δ 0 = H/200 et δ 0 = H/100.The three inter-storey drift limits (1%, 1.5%, 3%) correspond to different damage states that can be connected respectively to minor, moderate and severe damages.They can also refer to different seismic performance levels as immediate occupancy, life safety and collapse prevention.The use of the three limits allows observing the influence of failure events on derived fragility curves.For a given set of ground motions, the more the limit increases, the less the number of failure samples is.
Fragility curves of the frame, corresponding to three critical thresholds δ 0 obtained by the SSI, MLE and PSDM/PSCM methods, are presented in Figure 10.Note that these curves are established with 5000 accelerations, while for the MCS method a larger number of N MCS = 10 5 accelerations were used following the procedure presented in Section 2.3.11 intervals of PGA were organized as the same manner as the previous example.
The fragility curves of the SSI method are one more time, very different and distant from those of the MLE and PSDM/PSCM methods.Moreover, the fragility curves of the MLE and PSDM/PSCM methods are close to the MCS results, indicating their superior quality compared with the SSI method.Note that the SSI curves are well fitted with its discrete probabilities of failure obtained using scaling technique that are however far from the MCS reference results.Using the two comparison criteria in Section 2.3, MSE and RER are calculated and their results are presented in Table 3 and in Figure 11.
MSE and RER values are highest for the SSI method in all considered cases.They are lowest for the MLE method, except for a slightly higher value in the case of δ 0 = H/100.These remarks enable us to conclude that the MLE method gives the fragility curves with the best quality in this example.

Discussion
Three methods for the construction of fragility curves using numerical simulations were presented and applied to two numerical examples: linear/non-linear oscillators and an eight-storey bilinear frame.The quality of the obtained fragility curves is commented with respect to MCS results and two comparison criteria, Mean Square Error (MSE) and Relative ERror (RER).The results obtained clearly show a discrepancy between the fragility curves of the different methods.The fragility curves of the SSI method are always far from those of the MLE and PSDM/PSCM methods and also far from the MCS results.Moreover, its MSE and RER values are the highest in most cases.This suggests that the SSI method is not a good method for the construction of fragility curves from numerical simulations.
As mentioned in numerical examples, the SSI curves fitted well its discrete failure probabilities however these points obtained using the scaling technique are far from the failure probabilities given by the MCS reference results with unmodified ground motion time histories.Moreover, the MLE and PSDM/PSCM curves based on the same number of ground motion histories used in the SSI curves but using non scaling technique, are close to the MCS reference results.It can conclude that the scaling technique of ground motion histories in the SSI method is the reason of the biased results.Grigoriou [29] also has reported a difference of probabilistic distribution between original ground motion index and its scaled version.
The PSDM/PSCM method give slightly better fragility curves in some cases for the example of the oscillators, but the MLE method give better results for the eight-storey bilinear frame.Note that the PSDM/PSCM method is based on an approximation of seismic demand: maximum displacement of mass in the first example or maximum inter-storey drift in the second example, via a linear regression with PGA.This method is adequate if the correlation between seismic demand and seismic intensity is high.
The MLE method gives fragility curves which are in general close to the MCS results for the two considered examples.Its MSE and RER values are of the same order as those of the PSDM/PSCM method in the case of oscillators, but lower in the case of the eight-storey bilinear frame.Based on this remark, the MLE method is thus the recommended method for the construction of seismic fragility curves.This conclusion is in good agreement with a recent study by Mandal et al. [17].However, it should be noted that the accuracy of the MLE method decreases when threshold δ 0 increases, i.e. failure events become more and more rare.This issue is addressed in references [30,31].
Finally, it should be insisted that in this work, only the PGA indicator was considered for numerical examples because of its popularity and commodity in the fragility assessment.Other indicators than PGA and other structures for instance real buildings with sound finite element models should be further investigated.These perspectives are under study.

) where 1 [
.] is the indicator function, equal to 1 if the failure or damage state reached and equal to 0 otherwise.-Step 5: Repeat the steps from 2 to 4 for N p PGA values γ j , i.e., PGA = γ j , j = 1 . . .N p to obtain N p probability estimates pf (γ j ), j = 1 . . .N p .-Step 6: Perform a least square regression from N p

-Step 4 : 9 ) 3 .
Identify A m and β from the optimization problem below ( Âm , β) = arg min Am,β (−L(A m , β)) (Probabilistic Seismic Demand Model/Probabilistic Seismic Capacity Model (PSDM/PSCM) method [7, 8, 14, 15].The method is based on the definition of fragility curve given in Equation (2) and an additional assumption of log-normal distribution for seismic demand (D) and seismic capacity C i.e., D ∼ p LN (δ; S D (a), β D ) and C ∼ p LN (c; S C , β C ) where the lognormal distribution p LN (x; S m , β) is defined with two parameters S m and β as: ), S C = x 0 and β C = 0. -Step 5: Estimate two parameters of the fragility curve: Âm = exp ln S C − ln c 1 c 2 and β =
): D = max t |x(t)| and C = x 0 .Thus, if D ≥ C ⇔ max t |x(t)| ≥ x 0 ,the failure or damage state was reached and if D < C ⇔ max t |x(t)| <

Fig. 7 .
Fig. 7. Oscillators: MSE and RER results in three cases of failure threshold.
Step 1: Generate or choose a set of N independent ground accelerations a i (t), i = 1 . . .N. Their PGA values are PGA i = a i = max Step 2: Simulate structural responses subjected to the ground accelerations a i (t) and verify whether the failure or damage state is reached.Consider a Bernoulli random variable Y ; if a i (t) gives rise to failure or damage, Y is allocated a value y i = 1, and t |a i (t)|.- Step 1: Generate or choose a large number N MCS of independent ground accelerations a i (t), i = 1 . . .N MCS having the corresponding PGA value PGA i = a i = max Step 2: Organize these ground accelerations with respect to their PGA interval [a j − da j , a j + da j ].There are thus N I intervals, and in an interval j there are N j ground accelerations.Note that N MCS = t |a i (t)|.f :

Table 1 .
Oscillators: Results of comparison criteria.

Table 2 .
Lumped masses of the 8-storey frame structure.