Non-Gaussian Stochastic Equivalent Linearization Method for Inelastic Nonlinear Systems with Softening Behaviour , under Seismic Ground Motions

A non-Gaussian stochastic equivalent linearization (NSEL) method for estimating the non-Gaussian response of inelastic nonlinear structural systems subjected to seismic ground motions represented as nonstationary random processes is presented. Based on a model that represents the time evolution of the joint probability density function (PDF) of the structural response, mathematical expressions of equivalent linearization coefficients are derived. The displacement and velocity are assumed jointly Gaussian and the marginal PDF of the hysteretic component of the displacement is modeled by a mixed PDF which is Gaussian when the structural behavior is linear and turns into a bimodal PDF when the structural behavior is hysteretic. The proposed NSEL method is applied to calculate the response of hysteretic single-degree-of-freedom systems with different vibration periods and different design displacement ductility values. The results corresponding to the proposed method are compared with those calculated by means of Monte Carlo simulation, as well as by a Gaussian equivalent linearization method. It is verified that the NSEL approach proposed herein leads to maximum structural response standard deviations similar to those obtained with Monte Carlo technique. In addition, a brief discussion about the extension of the method tomuti-degree-of-freedom systems is presented.


Introduction
The method of stochastic equivalent linearization (SEL) is one of the most common methods within the approximate approaches for stochastic dynamic analysis of nonlinear systems.The SEL method has proven to be a versatile and efficient technique from the computation point of view.A great impulse to the method was given by Atalik and Utku [1], who demonstrated that for Gaussian responses the calculation of linearization coefficients can be performed in a very simple way.However, an important deficiency of the method has been found when the response of systems with nonlinear inelastic behavior subjected to intense excitations has been considered as Gaussian [2][3][4][5][6].The authors [7] have found that when the hypothesis related to the response of single-degree-of-freedom (SDOF) systems with hysteretic behavior has Gaussian distribution, the standard deviation of the displacement response can be underestimated in the order of 45%; consequently, the probability of failure widely deviates from the Monte Carlo simulation results, especially for systems with high displacement ductility demands.This is mainly due to the fact that a Gaussian probability distribution is assumed for all variables; however, the restoring force should lie on a finite region, implying that its probability density function (PDF) is non-Gaussian.Several attempts have been made in the last years to overcome the problem mentioned above [8][9][10].For example, Asano and Iwan [11] obtained alternate linearization results for a hysteretic system by using equations that were identical to the usual form on all subsets having a nonzero probability density in the nonlinear 2 Mathematical Problems in Engineering system but differed in the zero probability subsets.An ample review related to applications of statistical and equivalent linearization in the analysis of structure and mechanical nonlinear stochastic dynamic systems is found in [12].
A non-Gaussian equivalent linearization criterion (NSEL) that uses mathematical expressions for the linearization coefficients is presented herein.In order to derive the non-Gaussian equivalent linearization coefficients, a model appropriately representing the time evolution of the joint probability density function of the structural response of softening systems subjected to seismic ground motions is used.The parameters of the joint PDF depend on a function which measures the non-Gaussianity of the response process.
The structural softening models are of particular interest in earthquake engineering because they are commonly used to analyze buildings made with softening materials such as steel, reinforced concrete or masonry, which present saturation of the restoring force when deformation increases [13].The nonlinearities considered in the present study are only due to material behavior.
It is noticed that the NSEL methods mentioned in the first paragraph of this section do not take into account the time evolution character of the joint PDF of the structural response in the way that is formulated here.
The mathematical expressions corresponding to the non-Gaussian equivalent linearization coefficients derived here are applied to calculate the response of hysteretic SDOF systems with different vibration periods and design displacement ductility demands, subjected to seismic ground motions that are representative of a nonstationary ergodic random process.The comparison of the results with those obtained by means of Monte Carlo simulation shows that the proposed approach predicts with acceptable accuracy the maximum standard deviation of the structural response of hysteretic SDOF systems.
An extension of the proposed method to multi-degree-offreedom (MDOF) systems is outlined at the end of the paper.

Equation of Motion
The equation of motion of hysteretic SDOF system studied here is where ẍ , ẋ , and  represent the relative acceleration, velocity, and displacement, respectively. = /2 is the fraction of viscous critical damping;  = √/ is the circular frequency of vibration of the system;  is the mass;  is the coefficient of viscous damping;  is the initial stiffness. 2 is the ratio between the postyielding stiffness and the initial stiffness; () represents the ground acceleration;  is the time and  is the hysteretic component of displacement represented here by the Bouc-Wen model [14,15] which is defined by the following nonlinear differential equation: where , , , and  are parameters that define the size and the shape of the hysteretic cycles,  controls the smoothness of the transition from elastic to plastic structural behavior in such a way that the smaller , the smoother the transition, and  → ∞ corresponds to a bilinear model. and  define the softening or the hardening of the system; the first case corresponds to  +  > 0 while the second corresponds to  +  ≤ 0. Casciati and Faravelli [16], recommend considering  =  for steel and  = −2 for reinforced concrete.In this study softening systems with smooth transition represented by  = 1 are considered.

Stochastic Model of the Seismic Excitation
The seismic motion is modeled as a nonstationary Gaussian random process generated from a filtered white noise modulated in amplitude.The transfer function proposed by Clough and Penzien [17] is used here; therefore, the power spectral density of the seismic excitation is given by )  0 . ( The square of the amplitude modulating function is [18]  2 () =     +   Exp (−) . ( In practical applications, the amplitude of the power spectral density of the white noise ( 0 ) and the filter parameters (  ,   ,   , and   ) are determined with a nonlinear fitting of the power spectral density estimated by means of the Fourier spectrum of the basic seismic ground motion representative of the excitation process.The parameters , , , , and  of the modulating function are determined with a nonlinear fitting of the Arias intensity curve of the basic accelerogram [19].In this study, the excitation process is modulated only in amplitude; to study the effect of the frequency modulation on the nonlinear systems response the reader is referred to [20][21][22][23].
In the analysis presented below, the seismic record obtained in Mexico City at the Ministry of Communications and Transportation (SCT) during the September 19, 1985 earthquake is used to obtain the parameters of the stochastic model of the seismic excitation.The amplitude of the spectral density of the white noise is  0 = 7.2776 × 10 −4 cm 2 /s 3 , the filter parameters are   = 3.1017,   = 0.0220,   = 2.2988, and   = 0.0492, and those corresponding to the modulating function are  = 5.8161 × 10 48 cm 2 /s 4 ,  = −0.3388, = 2.166×10 47 ,  = 26.461,and  = −0.1258.The bimodal power spectral densities of the seismic motion and its fitted model are shown in Figure 1.The cumulative Arias intensity curve and its fitted model ( 0 ) are in Figure 2.

Modeling the Joint Probability Density Function of the Structural Response
The NSEL approach presented here is based on an appropriate representation of the time evolution of the joint PDF of the structural seismic response of softening systems with smooth transition.Results from Monte Carlo simulation analyses have shown that the PDF of the displacement,   (), and of the velocity,  Ẋ( ẋ ), of the structural system is Gaussian, and their evolution in time is as shown in Figures 3 and 4, respectively.However, the probability density of the restoring force and of the displacement hysteretic component,   (), changes with time because it depends on the nonlinearity of the structural behavior.At the beginning and at the end of the excitation, when the structure is under low seismic intensities,   () tends to be Gaussian because the structure presents linear behavior; however, as the seismic intensity increases and the structural behavior becomes nonlinear (hysteretic), then   () becomes non-Gaussian, as shown in Figure 5.It can be seen that when the system is vibrating within the nonlinear range  equals its maximum value   and its PDF presents concentrations at +  and −  .Figure 5 also shows that as the excitation intensity grows (up to  = 59 s) and consequently the displacement ductility demand increases,   () turns into a bimodal PDF.It is observed that the support of   () is −  ≤  ≤   , where the maximum value of  is given by [13,19]   = (   +  ) Considering the concentration of values of  in the vicinity of its maximum value   and the shape that such PDF adopts (see Figure 5 for  = 59 s), the following mathematical model is used: where    is a Gaussian density function with zero mean and standard deviation   ,    is a Gaussian density function with mean   (or −  ) and standard deviation   , and  is a weighting factor.The idea is that when the nonlinear behavior becomes important, the weighting of    , functions should increase while that corresponding to    should decrease.It is noticed that (6) uses the Gaussian density functions    and    and not Dirac's pulses as has been proposed by other authors [9].The advantage of using Gaussian density functions, instead of Dirac's pulses, is that it is possible to model more appropriately the time evolution of the joint PDF and, as a consequence, the NSEL approach leads to results that better approximate the Monte Carlo results.Furthermore, it is possible to get mathematical expressions of the equivalent linearization coefficients.Using the marginal PDFs   (),  Ẋ( ẋ ), and the information about the correlation structure of the random variables, it is possible to compute the joint PDF [24,25].Here, the joint PDF is obtained using the conditional PDF of  and ẋ given , ℎ(, ẋ | ), as follows: Based on the simulation results, it is reasonable to assume that ℎ(, ẋ | ) is jointly Gaussian.The parameters of ℎ(, ẋ | ), that is, conditional means and variances, are estimated assuming that the former are linear functions of , and the latter do not depend on , as it happens in a mean square estimation of normal variables [26].Additionally, it is considered that the response process has zero mean.
Based on ( 6) and ( 7), the following expression for the joint PDF is obtained: where   ,  ẋ , and   are the standard deviations of , ẋ , and , respectively;   ,  ẋ  , and   ẋ are the corresponding correlation coefficients; and should satisfy that the area under the PDF given by ( 6) is close to unity:

Evaluating the
where Erf(⋅) is the error function given by The above mentioned condition is necessary because ( 6) is not a bounded function.If ( 9) is satisfied, the second moment (in this case, equal to the variance  2  ) can be approximated as The four parameters   ,   , , and   could be computed by means of their first four moments; however, here they are determined by means of nondimensional coefficients that depend on the system response non-Gaussianity level which is represented by [13] From (12) it is observed that  is close to zero when the response is Gaussian, that is, when the system strength is high (where   is high) and the excitation or response is low (and as a consequence   is low).On the other hand, if the response is non-Gaussian, that is, if the system strength is low (where   is low) and the excitation or response is high (and consequently   is high), then  is high.
Here, the following mathematical expressions for the normalized parameters   /  ,   /  ,   /  , and  are proposed: where   ,   ,   , and   ,  = 1, 2, 3 and 4, are constants to be determined.The general behavior of ( 13) is presented in Figures 6(a Equations (13) show that the parameters   ,   , , and   are time dependent; thus, the proposed approach considers the time evolutionary character of the joint PDF of the response, while other methods proposed in the literature do not take it into account.Figure 7 shows an example of the time evolution of   () given by (6).

Case of Narrow-Band Seismic Inputs.
Parameters   ,   ,   , and   ,  = 1, 2, 3 and 4 of (13) corresponding to the response of single-degree-of-freedom systems (SDOF) systems subjected to the action of narrow-band seismic inputs, were calculated by means of Monte Carlo simulation analysis.The critical case, where the natural vibration period of SDOF system is equal to the dominant period of seismic input (in this case  = 2.1 s, see Figure 1) and the SDOF system has a high design ductility demand  = 4, was considered.The design ductility factor  is defined here as the ratio between the expected maximum and the yield displacement of the system.In this study, the first one is found by means of a stationary Gaussian equivalent linearization analysis, and the latter is estimated by means of the initial stiffness and the yield force of the system.
The SDOF system was subjected to the action of 50,000 artificial accelerograms based on the narrow-band motion recorded in SCT on September 19,1985 (see Section 3).For this case, it was found that  lies within the interval 0 <  ≤ 0.1 and that   >   leads to good results.The parameter values obtained from the analysis are shown in Table 1.It is noticed that parameters in Table 1 are proposed for narrow-band seismic inputs and are independent of any other parameter, such as structural vibration period and design ductility factor.

Non-Gaussian Equivalent Stochastic Linearization
In the equivalent linearization method, (2) is replaced by a linear differential equation as follows [5,27]: where   ,   , and   are linearization coefficients.In order to minimize the expected value of the square of the difference of ( 2) and ( 14), the linearization coefficients must satisfy [5,28,29] where  is the response vector  = [ ẋ ]  ,   is the vector of equivalent linearization coefficients given by   = [      ]  , and ℎ is given by (2).
When the joint PDF proposed in this paper (see ( 8)) is substituted in (15), mathematical expressions for the linearization coefficients are derived (see Appendix A).In particular, for  = 1, which corresponds to SDOF systems with softening behavior, the equivalent linearization coefficients are as follows: where

Covariance of the Response of the Hysteretic System
The covariance matrix of the response of the hysteretic SDOF system, , is calculated by solving the following equation [28,29]: where  is a matrix depending on the mechanical and geometrical properties of the system, the linearization coefficients, the modulating function coefficients, and the filter parameters.In this study,  is given by   is a 7 × 7 matrix of zeros, except the element   (7, 7) which is equal to the amplitude of the white noise spectral density ( 0 ).
In order to compute Σ(), it is essential to start the integration process of ( 17) when the Clough-Penzien filter response has reached the stationary stage.Thus, the initial condition for Σ() should correspond to the covariance of the stationary response of the Clough-Penzien filter as described in Appendix B [30].

Verification of the Accuracy of the Proposed NSEL Approach
In this section the standard deviation of the response of several hysteretic SDOF systems is calculated by means of the proposed NSEL criterion using the mathematical expressions obtained above (( 13) and ( 16) using values of Table 1).The SDOF systems have vibration periods  = 2.1 s and  = 4.0 s and design ductility factors  = 1.5 and  = 4, and they are subjected to a nonstationary random process with the statistical properties of the seismic motion recorded  in Mexico City at the SCT during the September 19, 1985 earthquake (see Section 3).
The response of the SDOF hysteretic systems is calculated by means of SEL and alternatively by Monte Carlo simulation using 50,000 simulated accelerograms.Two SEL criteria are used: (a) assuming Gaussian response [1] and (b) assuming non-Gaussian response using the criterion and the mathematical expressions proposed herein.
The statistical response of the SDOF system with  = 2.1 s and  = 1.5 is shown in Figures 8, 9, and 10, and that corresponding to the system with  = 2.1 s and  = 4.0 is in Figures 11, 12, and 13.Figures 8 to 13 show that the proposed approach predicts with reasonable accuracy the maximum response; however, as well as the Gaussian method, it does not predict the permanent drift of the hysteretic system especially for high values of the design ductility (see Figures 8 and 11).The statistical response of the SDOF system with  = 4.0 s and  = 1.5 is shown in Figures 14, 15, and 16, and that corresponding to the system with  = 4.0 s and  = 4.0 is in Figures 17,18,and 19.As in the case of SDOF system with  = 2.1 s, the proposed approach predicts with reasonable accuracy the maximum response of the SDOF system with  = 4.0 s; however it does not predict with good accuracy the permanent drift.
The relative error () in the calculation of the maximum standard deviation of the response (  ,  ẋ and   ) is obtained as follows: response obtained from linearization − response obtained from simulation response obtained from simulation × 100.
The relative errors () of the maximum response, corresponding to all the cases analyzed, are shown in Tables 2 and 3 corresponding to structural vibration periods of 2.1 and 4.0 s, respectively.A negative value of  indicates that the equivalent linearization underestimates the response.Tables 2 and 3 show that the proposed NSEL method leads to results close to those obtained by means of Monte Carlo simulation.

Extension to Multiple-Degree-of-Freedom Systems
In this section an extension of the proposed method to multiple-degree-of-freedom (MDOF) systems is presented.Consider a nonlinear dynamic MDOF structural system subjected to a vector input P() with an -dimensional response q() satisfying the nonlinear differential equation: Φ (q () , q () , q ()) = P () . ( Here bold letters are related to MDOF systems.The equivalent linear system of ( 20) is defined as M e q () + C e q () + K e q () = P () , where M e , C e , and K e are defined as the equivalent mass, viscous damping, and stiffness matrices of  × .They are determined minimizing the mean square error where  is the -dimensional vector:  = Φ (q, q , q ) − (M e q + C e q + K e q) .( Denoting by Q and H e the 3-dimensional response vector and the 3 ×  equivalent structural matrix given, respectively, by it can be shown that H e can be obtained from [5,28,29] Herein, the random excitation P() is modelled as a linear combination of the responses of linear filters to white or shot noise.With the aim of calculating the statistical second-order response of the MDOF structure, it is convenient to express (21) in state space form: where Z() = [q(), q ()]  is the 2 × 1 state vector and A e the time-dependent system matrix given by and the matrix B contains the coefficients of the linear combination of the response Z F of the filters.The dynamics of the latter is governed by where the matrix D contains the coefficients of the filters and the vector W  () the white or shot noise.Equations ( 25) and ( 27) can be rewritten as where Matrix W has the same structure as W  but with higher dimension.Thus, (28) controls the response of an augmented system consisting of the linear filters and the structural system in series [5].Finally, the covariance matrix of the response of the MDOF system is calculated by solving where S F is a matrix of zeros, except the last element which is equal to the amplitude of the white noise spectral density ( 0 ).The mathematical expressions given in this section correspond to a MDOF structural system with -dimensional response, while those given in the previous sections correspond to the special case of a SDOF with 2-dimensional response (q() = [, ]  ).For example (20) took the particular form of ( 1) and ( 2); (21) took the particular form of ( 14); (24) was reduced to (15) while (30) took the particular form of (18).

Conclusions
A model appropriately representing the time evolution of the joint PDF of the structural response of softening systems   subjected to seismic ground motions was proposed.The displacement and velocity was assumed jointly Gaussian and the marginal PDF of the hysteretic component of the displacement was modeled by a mixed PDF which is Gaussian when the structural behavior is linear and turns into a bimodal PDF when the structural behavior is hysteretic.The NSEL method proposed in this study uses mathematical expressions for the equivalent linearization coefficients (16) corresponding to the non-Gaussian response of inelastic nonlinear systems.The mathematical expressions depend on the parameters of the joint PDF of the structural response, and they take into account the time evolutionary character of the PDF.It is noticed that other NSEL methods found  in the literature have proposed mathematical expressions for the equivalent linearization coefficients; however, they do not consider the time evolution of the joint The mathematical expressions derived here ( 16) can be applied to obtain the non-Gaussian response of any kind of SDOF softening system with smooth transition from elastic to plastic behavior ( = 1).It is possible to derive expressions of the linearization coefficients for other values of .
The NSEL method proposed was applied to calculate the statistical response of SDOF systems with different vibration periods and different design ductility factors.The results show that, for the SDOF systems analyzed, the proposed criterion estimates with good accuracy the time evolution of the structural response; however, it does not predict the permanent drift exhibited by hysteretic systems with high design ductility values.
It is concluded that the proposed NSEL approach is useful to estimate the maximum standard deviation of the structural response of softening systems subjected to seismic ground motions.The relative errors in the calculation of the maximum standard deviation of the displacement, velocity, The parameters   ,   ,   , and   ,  = 1, 2, 3, 4, in (13) can be found for other types of seismic excitation (e.g., wide-band ground motions).

A. Linearization Coefficients
This appendix presents the derivation of the equivalent linearization coefficients given by (16).
The linearization coefficients are obtained from where  is the response vector and   is the vector of linearization coefficients, given by  = [ ẋ ]  and The inverse of the covariance matrix is where  = 1 −  2  −  2 ẋ  −  2  ẋ + 2  ẋ    ẋ  .Substituting (A.2) in the second factor of (A.1), the following is obtained: where Considering that   Ẋ is modeled by (8),  Ẋ ,   Ẋ, and  Ẋ can be easily obtained, thus making it possible to obtain mathematical expressions for the linearization coefficients.

B. Initial Condition of the Covariance Matrix of the Response
Following the state space approach, the state equation of a linear system subjected to a nonwhite excitation modeled as a filtered white noise is [29] Z Matrix W has the same structure as W  but with a higher dimension.To avoid introducing the effect arising from the transient response of the filter, the output of the filter Z F must be allowed to reach the stationary phase before it is multiplied by the modulating function [5].Analytically, it is reached by selecting the suitable initial conditions for the covariance matrix Σ of the overall state variable vector Z, which is governed by the following differential equation: ].If the original system is at rest, at  = 0, then the response will be zero and all the elements of Σ  and Σ   will be zero; however the elements of Σ     , at  = 0, will correspond to the stationary response of the filter which is governed by DΣ      + Σ      D  + P  = 0, (B.8) where P  has the same form as P but with a lower dimension.
For the case when the Clough-Penzien filter is used, the solution of (B.

2 /s 3 )Figure 1 :
Figure 1: Bimodal power spectral density of the accelerogram and its corresponding fitted model.

Figure 3 :
Figure 3: Time evolution of the probability density function of the displacement response.

Figure 4 :
Figure 4: Time evolution of the probability density function of the velocity response.

Figure 5 :
Figure 5: Time evolution of the probability density function of the hysteretic component of displacement response.
)  = AZ () + BZ F () , (B.1) where Z() is the state vector, the matrix A is the system matrix, and the matrix B contains the coefficients of the linear combination of the response Z F of the filter.The response of the filters is controlled by Z F ()  = DZ F () + W  () , (B.2) where the matrix D contains the coefficients of the filters and W  () is a vector with zeros except the last element which correspond to white noise.(B.1) and (B.2) can be rewritten as u  = H () u + W, Σ ()  = H () Σ () + Σ () H  () + 2S F , (B.6) where S F is a matrix of zeros except the last element which is equal to the amplitude of the white noise spectral density ( 0 ).It is convenient to introduce a partition in Σ in agreement with (B.4).Thus, Σ = [ Σ  Σ   Σ    Σ     ] , (B.7) where Σ  = [ZZ  ], Σ   = [ZZ  F ], and Σ     = [Z F Z  F

Table 2 :
Relative errors of the maximum response of SDOF system with  = 2.1 s.

Table 3 :
Relative errors of the maximum response of SDOF system with  = 4.0 s.