Heavy-tailed response of structural systems subjected to stochastic excitation containing extreme forcing events

Hea vy-tailed resp onse of structural systems sub jected to sto c hastic excitation con taining extreme forcing ev en ts Han Kyul Jo o, Mustafa A. Mohamad, Themistoklis P . Sapsis ∗ Departmen t of Mechanical Engineering, Massac husetts Institute of T echnology , 77 Massac husetts A v e., Cam bridge, MA 02139 Septem b er 28, 2018 Abstract W e c haracterize the complex, heavy-tailed probability distribution functions ( pdf ) de- scribing the resp onse and its local extrema for structural systems sub jected to random forcing that includes extreme even ts. Our approac h is based on the recen t probabilistic decomp osition-syn thesis technique in [ 21 ], where we decouple rare even ts regimes from the bac kground fluctuations. The result of the analysis has the form of a semi-analytical ap- pro ximation form ula for the pdf of the resp onse (displacement, velocity , and acceleration ) and the p df of the lo cal extrema. F or sp ecial limiting cases ( lightly damped or hea vily damp ed systems ) our analysis provides fully analytical approximations. W e also demon- strate how the metho d can be applied to high dimensional structural systems through a tw o-degrees-of-freedom structural system undergoing rare even ts due to intermitten t forcing. The derived formulas can b e ev aluated with very small computational cost and are shown to accurately capture the complicated hea vy-tailed and asymmetrical features in the probability distribution man y standard deviations aw ay from the mean, through comparisons with exp ensive Monte-Carlo simulations. Keyw ords Rare and extreme even ts; In termitten tly forced structural systems; Heavy-tails; Colored sto chastic excitation; Random impulse trains. 1 In tro duction A large class of physical systems in engineering and science can be mo deled by sto chastic differen tial equations. F or many of these systems, the dominan t source of uncertaint y is due to the forcing which can be described by a stochastic pro cess. Applications include ocean engineering systems excited b y w ater w av es ( such as ship motions in large w a ves [ 23 , 5 , 9 , 8 ] or high sp eed crafts sub jected to rough seas [ 26 , 25 ] ) and rare ev ents in structural systems ( suc h as b eam buc kling [ 1 , 17 ], vibrations due to earthquak es [ 16 , 7 ] and wind loads [ 18 , 31 ] ). F or all of these cases it is common that hidden in the otherwise predictable magnitude of the fluctuations are extreme even ts, i.e. abnormally large magnitude forces whic h lead to rare resp onses in the dynamics of the system ( figure 1 ). Clearly , these ev ents m ust b e adequately taken in to account ∗ Corresponding author: sapsis@mit.edu , T el: ( 617 ) 324-7508, F ax: ( 617) 253-8689 1 t 1900 2000 2100 2200 2300 2400 2500 -0.2 -0.1 0 0.1 0.2 t 1900 2000 2100 2200 2300 2400 2500 -0.02 0 0.02 0.04 0.06 F(t) time x(t) F r impulse event Ex ci ta ti on Re sp on se x r extreme response Figure 1: ( T op ) Background sto chastic excitation including impulsive loads in ( red ) upw ard arro ws. ( Bottom ) System response displacement. for the effective quantification of the reliability prop erties of the system. In this work, w e develop an efficient metho d to fully describ e the probabilistic resp onse of linear structural systems under general time-correlated random excitations containing rare and extreme even ts. Systems with forcing having these c haracteristics p ose significant c hallenges for traditional uncertain t y quan tification sc hemes. While there is a large class of metho ds that can accurately resolv e the statistics asso ciated with random excitations ( e.g. the F okker-Planc k equation [ 30 , 29 ] for systems excited by white-noise and the joint resp onse-excitation metho d [ 28 , 34 , 12 , 3 ] for arbitrary sto chastic excitation) these hav e imp ortan t limitations for high dimensional systems. In addition, even for lo w-dimensional systems determining the part of the probabilit y density function ( p df ) asso ciated with extreme ev en ts p oses imp ortant numerical challenges. On the other hand, Gaussian closure sc hemes and moment equation or cumulan t closure metho ds [ 6 , 35 ] either cannot “see” the rare even ts completely or they are very expensive and require the solution of an inv erse moment problem in order to determine the pdf of interest [ 2 ]. Similarly , approac hes relying on p olynomial-chaos expansions [ 36 , 37 ] ha v e b een sho wn to ha ve important limitations for systems with intermitten t resp onses [ 19 ]. Another p opular approach for the study of rare even t statistics in systems under in termittent forcing is to represen t extreme even ts in the forcing as iden tically distributed independent impulses arriving at random times. The generalized F okker-Planc k equation or Kolmogoro v- F eller (KF ) equation is the go verning equation that solv es for the ev olution of the response p df under Pois son noise [ 29 ]. How ever, exact analytical solutions are av ailable only for a limited n um b er of sp ecial cases [ 33 ]. Although alternative metho ds such as the path in tegral metho d [ 14 , 11 , 4 ] and the sto c hastic a veraging metho d [ 39 , 38 ] may b e applied, solving the FP or KF equations is often very expensive [ 20 , 10 ] even for v ery lo w dimensional systems. Here, we consider the problem of quan tification of the resp onse p df and the pdf associated with lo cal extrema of linear systems subjected to stochastic forcing containing extreme ev ents based on the recently form ulated probabilistic-decomposition syn thesis (PDS ) metho d [ 21 , 22 ]. The approac h relies on the decomp osition of the statistics into a ‘non-extreme core’, t ypically 2 Gaussian, and a heavy-tailed component. This decomposition is in full corresp ondence with a partition of the phase space in to a ‘stable’ region where we do not ha ve rare even ts and a region where non-linear instabilities or external forcing lead to rare transitions with high probability . W e quan tify the statistics in the stable region using a Gaussian appro ximation approac h, while the non-Gaussian distribution asso ciated with the intermitten tly unstable regions of phase space is p erformed taking into account the non-trivial character of the dynamics ( either b ecause of instabilities or external forcing ). The probabilistic information in the t w o domains is analytically syn thesized through a total probability argumen t. W e b egin with the simplest case of a linear, single-degree-of-freedom ( SDOF ) system and then formulate the method for m ulti-degree-of-freedom systems. The main r esult of our work is the derivation of analytic/semi-analytic appr oximation formulas for the r esp onse p df and the p df of the lo c al extr ema of intermittently for c e d systems that c an ac cur ately char acterize the statistics many standar d deviations away fr om the me an . Although the systems considered in this w ork are linear the metho d is directly applicable for nonlinear structural systems as w ell. This approac h circumv ents the challenges that rare even ts p ose for traditional uncertaint y quan tification schemes, in particular the computational burden asso ciated when dealing with rare even ts in systems. W e emphasize the statistical accuracy and the computational efficiency of the presen ted approac h, which we rigorously demonstrate through extensive comparisons with direct Monte-Carlo sim ulations. In brief, the principal con tributions of this pap er are: • Analytical (under certain conditions) and semi-analytical ( under no restrictions ) pdf expressions for the resp onse displacemen t, velocity and acceleration for single-degree-of- freedom systems under intermitten t forcing. • Semi-analytical p df expressions for the v alue and the lo cal extrema of the displacement, v elo cit y and acceleration for multi-degree-of-freedom systems under in termittent forcing. The pap er is structured as follo ws. In section 2 w e pro vide a general formulation of the prob- abilistic decomp osition-synthesis method for the case of structural systems under in termittent forcing. Next, in section 3 , we apply the developed metho d analytically , whic h is p ossible for t w o limiting cases: underdamp ed systems with ζ  1 or o verdamped with ζ  1 , where ζ is the damping ratio. The system w e consider is excited by a forcing term consisting of a background time-correlated sto chastic process sup erimp osed with a random impulse train ( describing the rare and extreme comp onent). W e give a detailed deriv ation of the resp onse pdf of the system ( displacement, v elo city and acceleration) and compare the results with exp ensive Mon te-Carlo sim ulations. In section 4 , w e sligh tly modify the dev elop ed form ulation to derive a semi-analytical sc heme considering the same linear system but without any restriction on the damping ratio ζ , demonstrating global applicability of our approac h. In section 5 , w e demonstrate applicability of our method for multiple-degree-of-freedom systems and in section 6 we presen t results for the lo cal extremes of the resp onse. Finally , w e offer concluding remarks in section 7 . 2 The probabilistic decomp osition-synthesis metho d for in termitten tly forced structural systems W e provide with a brief presentation of the recen tly developed probabilistic decomp osition- syn thesis ( PDS ) method adapted for the case of in termittently forced linear structural systems [ 21 ]. W e consider the following vibrational system, M ¨ x ( t ) + D ˙ x ( t ) + K x ( t ) = F ( t ) , x ( t ) ∈ R n , ( 1 ) 3 where M is a mass matrix, D is the damping matrix, and K is the stiffness matrix. W e assume F ( t ) is a sto chastic forcing with in termitten t c haracteristics that can be expressed as F ( t ) = F b ( t ) + F r ( t ) . ( 2 ) The forcing consists of a background comp onent F b of c haracteristic magnitude σ b and a r ar e and extr eme comp onen t F r with magnitude σ r  σ b . The comp onen ts F b and F r ma y b oth b e ( weakly) stationary stochastic pro cesses, while the sum of the tw o pro cesses can be, in general, non-stationary . This can b e seen if we directly consider the sum of tw o ( weakly) stationary pro cesses x 1 and x 2 , with time correlation functions Corr x 1 ( τ ) and Corr x 2 ( τ ) , resp ectively . Then for the sum z = x 1 + x 2 w e ha v e Corr z ( t, τ ) = Corr x 1 ( τ ) + E [ x 1 ( t ) x 2 ( t + τ )] + E [ x 1 ( t + τ ) x 2 ( t )] + Corr x 2 ( τ ) . Therefore, the process z is stationary if and only if the cross-co v ariance terms E [ x 1 ( t ) x 2 ( t + τ )] and E [ x 1 ( t + τ ) x 2 ( t )] are functions of τ only or they are zero ( i.e. x 1 and x 2 are not correlated ). F or the case where the excitation is given in terms of realizations, i.e. time-series, one can first separate the extreme even ts from the stationary background b y applying time-frequency analysis metho ds ( e.g. wa velets [ 32 ] ). Then the stationary background can b e approximated with a Gaussian stationary sto chastic pro cess ( with prop erly tuned co v ariance function) while the rare even t comp onent can b e represented with a P oisson process with properly c hosen parameters that represent the c haracteristics of the extreme forcing ev ents (frequency and magnitude ). T o apply the PDS metho d we also decomp ose the resp onse into t wo terms x ( t ) = x b ( t ) + x r ( t ) , ( 3 ) where x b accoun ts for the background state ( non-extreme ) and x r captures extreme resp onses ( due to the in termittent forcing ) – see figure 2 . More precisely x r is the system resp onse under t w o conditions: ( 1 ) the forcing is given by F = F r ( i.e. w e hav e an impulse ) and ( 2 ) the norm of the resp onse is greater than the bac kground response fluctuations according to a giv en criterion, e.g. k x k > γ . Ho wev er, as w e will see in the follo wing sections other criteria ma y be used. These rare transitions o ccur when we ha v e an impulse and they also include a phase that relaxes the system bac k to the bac kground state x b . The bac kground comp onen t x b corresp onds to system resp onse without rare even ts x b = x − x r , and in this regime the system is primarily gov erned b y the background forcing term F b . W e require that rare even ts are statistically indep endent from each other. In the generic form ulation of the PDS we also need to assume that rare ev ents ha ve negligible effects on the bac kground state x b but here this assumption is not necessary due to the linear c haracter of the examples considered. How ever, in order to apply the metho d for general nonlinear structural systems we need to hav e this condition satisfied. W e also need the dynamics to b e ergodic while the statistics we aim to appro ximate refer to long time av erages [ 21 ]. Next, we fo cus on the statistical characteristics of an individual mo de u ( t ) ∈ R of the original system in equation ( 1 ). The first step of the PDS method is to quan tify the conditional statistics of the rare even t regime. When the system enters the rare even t resp onse at t = t 0 w e will ha v e an arbitrary background state u b as an initial condition at t 0 and the problem will b e form ulated as: ¨ u r ( t ) + λ ˙ u r ( t ) + k u r ( t ) = F r ( t ) , with u r ( t 0 ) = u b and F = F r for t > t 0 . ( 4 ) Under the assumption of indep endent rare ev ents w e can use equation ( 4 ) as a basis to deriv e analytical or numerical estimates for the statistical resp onse during the rare even t regime. 4 F r F b x r x b Excitation Response P( x r | ∙ ) Rare component Background component P( x b | ∙ ) Figure 2: Schematic represen tation of the PDS method for an in termittently forced system. The background component, on the other hand, can b e studied through the equation, M ¨ x b ( t ) + D ˙ x b ( t ) + K x b ( t ) = F b ( t ) . ( 5 ) Because of the non-intermitten t character of the resp onse in this regime, it is sufficient to obtain the low-order statistics of this system. In the context of vibrations it is reasonable to assume that F b ( t ) follo ws a Gaussian distribution in which case the problem is straightforw ard. Consequently , this step provides us with the statistical steady state probability distribution for the mo de of in terest under the condition that the dynamics ‘liv e’ in the stochastic bac kground. Finally , after the analysis of the t wo regimes is completed we can synthesize the results through a total probability argumen t f ( q ) = f ( q | k u k > γ , F = F r ) | {z } rare events P r + f ( q | F = F b ) | {z } background (1 − P r ) , ( 6 ) where q ma y b e an y function of in terest in v olving the resp onse. In the last equation, P r denotes the o v erall rare ev ent probabilit y . This is defined as the probability of the response exceeding a threshold γ b ecause of a rare even t in the excitation: P r ≡ P ( k u k > γ , F = F r ) = 1 T Z t ∈ T 1 ( k u k > γ , F = F r ) dt, ( 7 ) where 1 ( · ) is the indicator function. The rare even t probability measures the total duration of the rare even ts taking into accoun t their frequency and duration. The utilit y of the presented decomp osition is its flexibilit y in capturing rare resp onses, since w e can account for the rare ev en t dynamics directly and connect their statistical prop erties directly to the original system resp onse. 2.1 Problem formulation for linear SDOF systems T o demonstrate the metho d we b egin with a very simple example and we consider a single- degree-of-freedom linear system ¨ x + λ ˙ x + k x = F ( t ) , ( 8 ) 5 where k is the stiffness, λ is the damping, ζ = λ/ 2 √ k is the damping ratio. F or what follows w e adopt the standard definitions: ω n = √ k , ω o = ω n p ζ 2 − 1 , and ω d = ω n p 1 − ζ 2 . F ( t ) is a sto c hastic forcing term with intermitten t characteristics, whic h can be written as F ( t ) = F b ( t ) + F r ( t ) . ( 9 ) Here F b is the background forcing comp onent that has a characteristic magnitude σ b and F r is a rare and large amplitude forcing comp onent that has a c haracteristic magnitude σ r , whic h is muc h larger than the magnitude of the background forcing, σ r  σ b . Despite the simplicity of the system, its resp onse may feature a significantly complicated statistical structure with hea vy-tailed c haracteristics. F or concreteness, w e consider a protot yp e system motiv ated from o cean engineering appli- cations, mo deling base excitation of a structural mo de: ¨ x + λ ˙ x + k x = ¨ h ( t ) + N ( t ) X i =1 α i δ ( t − τ i ) , 0 < t ≤ T . ( 10 ) Here h ( t ) denotes the zero-mean background base motion term (having opposite sign from x ) with a Pierson-Mosko witz sp ectrum: S hh ( ω ) = q 1 ω 5 exp  − 1 ω 4  , ( 11 ) where q controls the magnitude of the forcing. The second forcing term in equation ( 10 ) describ es rare and extreme ev en ts. In particular, w e assume this comp onent is a random impulse train ( δ ( · ) is a unit impulse ), where N ( t ) is a P oisson counting pro cess that represen ts the num b er of impulses that arrive in the time interv al 0 < t ≤ T , α is the impulse mean magnitude ( characterizing the rare even t magnitude σ r ), whic h we assume is normally distributed with mean µ α , v ariance σ 2 α and indep endent from the state of the system. In addition, the arriv al rate is constant and giv en by ν α ( or by the mean arriv al time T α = 1 /ν α so that impulse arriv al times are exp onen tially distributed τ ∼ e T α ). W e tak e the impulse mean magnitude as being m -times larger than the standard deviation of the excitation velocity ˙ h ( t ) : µ α = mσ ˙ h , with m > 1 , ( 12 ) where σ ˙ h is the standard deviation of ˙ h ( t ) . This protot yp e system is widely applicable to n umer- ous applications, including structures under wind excitations, systems under seismic excitations, and vibrations of high-sp eed crafts and road v ehicles [ 29 , 30 , 26 ]. 3 Analytical p df of SDOF systems f or limiting cases In this section, we apply the probabilistic decomposition-synthesis method for the special cases ζ  1 and ζ  1 to derive analytical appro ximations for the p df of the displacemen t, velocity and acceleration. W e p erform the analysis first for the resp onse displacemen t and by w ay of a minor mo dification obtain the p df for the v elo cit y and acceleration. 3.1 Bac kground resp onse p df Consider the statistical resp onse of the system to the background forcing component, ¨ x b + λ ˙ x b + k x b = ¨ h ( t ) , ( 13 ) 6 due to the Gaussian character of the statistics the resp onse is fully c haracterized by it’s spectral resp onse. The sp ectral density of the displacement, v elo city and acceleration of this system are giv en b y , S x b x b ( ω ) = ω 4 S hh ( ω ) { k − ω 2 + λ ( j ω ) } 2 , S ˙ x b ˙ x b ( ω ) = ω 2 S x b x b ( ω ) , S ¨ x b ¨ x b ( ω ) = ω 4 S x b x b ( ω ) . ( 14 ) Th us, w e can obtain the v ariance of the response displacement, v elo cit y and acceleration: σ 2 x b = ∞ Z 0 S x b x b ( ω ) dω , σ 2 ˙ x b = ∞ Z 0 S ˙ x b ˙ x b ( ω ) dω , σ 2 ¨ x b = ∞ Z 0 S ¨ x b ¨ x b ( ω ) dω . ( 15 ) Moreo v er, the en velopes are Ra yleigh distributed [ 15 ]: u b ∼ R ( σ x b ) , ˙ u b ∼ R ( σ ˙ x b ) , ¨ u b ∼ R ( σ ¨ x b ) . ( 16 ) 3.2 Analytical p df for the underdamp ed case ζ  1 Because of the underdamp ed character of the resp onse for the case of ζ  1 , we fo cus on deriving the statistics of the lo cal extrema. T o this end, we will b e presenting results for the statistics of the env elop e of the resp onse. 3.2.1 Rare ev ents resp onse T o estimate the rare even t resp onse we tak e into accoun t the non-zero background v elo city of the system ˙ x b at the moment of impact, as well as the magnitude of the impact, α . The actual v alue of the resp onse x b is considered negligible. F or this case, taking in to account ζ  1 , we ha v e the en velopes of the resp onse ( displacement, v elo city , acceleration ) during the rare ev ent giv en b y ( see app endix A for details ), u r ( t ) ' | ˙ x b + α | ω d e − ζ ω n t , ˙ u r ( t ) ' | ˙ x b + α | e − ζ ω n t , ¨ u r ( t ) ' ω d | ˙ x b + α | e − ζ ω n t . ( 17 ) In equation ( 17 ) the tw o contributions ˙ x b and α in the term ˙ x b + α are b oth Gaussian distributed and indep endent and therefore their sum is also Gaussian distributed as: η ≡ ˙ x b + α ∼ N ( µ α , σ 2 ˙ x b + σ 2 α ) . ( 18 ) Therefore, the distribution of the quantit y | η | is given by the follo wing folded normal distribution: f | η | ( n ) = 1 σ | η | √ 2 π  exp  − ( n − µ α ) 2 2 σ 2 | η |  + exp  − ( n + µ α ) 2 2 σ 2 | η |  , 0 < n < ∞ , ( 19 ) where σ | η | = q σ 2 ˙ x b + σ 2 α . 3.2.2 Rare ev ent probabilit y Next, we compute the rare ev ent probability , which is the total duration of the rare even ts o v er a time interv al, defined in equation ( 7 ). This will be done b y employing an appropriate description for extreme even ts. One p ossible option is to set an absolute threshold γ . How ev er, in the curren t context it is more con venien t to set this threshold relativ e to the local maximum of the resp onse. Sp ecifically , the time duration τ e a rare resp onse takes to return back to the 7 bac kground state will be given b y the duration starting from the initial impulse ev en t time ( t 0 ) to the p oint where the response has deca yed bac k to ρ c ( or 100 ρ c % ) of its absolute maximum; here and throughout this manuscript w e tak e ρ c = 0 . 1 . This is a v alue that we considered without any tuning. W e emphasize that the derived appro ximation is not sensitiv e to the exact v alue of ρ c as long as this has b een c hosen within reasonable v alues. This means that τ e is defined by u r ( τ e + t 0 ) = ρ c u r ( t 0 ) , ( 20 ) W e solv e the ab ov e using the derive d envelop es to obtain τ e = − 1 ζ ω n log ρ c . ( 21 ) W e note that due to the linear character of the system the typic al dur ation τ e is indep endent of the b ackgr ound state or the imp act intensity . With the obtained v alue for τ e w e compute the probabilit y of rare even ts using the frequency ν α ( equal to 1 /T α ) : P r = ν α τ e = τ e /T α . ( 22 ) Note that based on our assumption that extreme ev ents are rare enough to be statistical indep enden t, the ab ov e probability is m uch smaller than one. 3.2.3 Conditional p df for rare even ts W e now proceed with the deriv ation of the p df in the rare ev en t regime. Consider again the resp onse displacement during a rare ev ent, u r ( t ) ∼ | η | ω d e − ζ ω n t 0 , (23 ) here t 0 is a random v ariable uniformly distributed b etw een the initial impulse ev ent and the end time τ e ( equation ( 21 ) ) when the response has relaxed bac k to the bac kground dynamics: t 0 ∼ Uniform (0 , τ e ) . ( 24 ) W e condition the rare even t distribution as follows, f u r ( r ) = Z f u r || η | ( r | n ) f | η | ( n ) dn, ( 25 ) where w e hav e already deriv ed the pdf for f | η | in equation ( 19 ). What remains is the deriv ation of the conditional p df for f u r || η | . By conditioning on | η | = n , w e find the derived distribution for the conditional pdf giv en b y f u r || η | ( r | n ) = 1 r ζ ω n τ e  s  r − n ω d e − ζ ω n τ e  − s  r − n ω d  , ( 26 ) where s ( · ) denotes the step function, whic h is equal to 1 when the argumen t is greater or equal to 0 and 0 otherwise. W e refer to appendix B for a detailed deriv ation. Using the equations ( 18 ) and ( 26 ), in equation ( 25 ) we obtain the final result for the rare ev en t distribution for resp onse displacement as f u r ( r ) = Z f u r || η | ( r | n ) f | η | ( n ) dn, ( 27 ) = 1 r ζ ω n σ | η | τ e √ 2 π ∞ Z 0  exp  − ( n − µ α ) 2 2 σ 2 | η |  + exp  − ( n + µ α ) 2 2 σ 2 | η |  ×  s  r − n ω d e − ζ ω n τ e  − s  r − n ω d  dn. ( 28 ) 8 3.2.4 Summary of results for the underdamp ed case Displacemen t Env elop e Finally , combining the results of sections 3.1 , 3.2.2 , and 3.2.3 using the total probability la w, f u ( r ) = f u b ( r )(1 − P r ) + f u r ( r ) P r , ( 29 ) w e obtain the desired env elop e distribution for the displacemen t of the resp onse f u ( r ) = r σ 2 x b exp  − r 2 2 σ 2 x b  (1 − ν α τ e ) + ν α τ e r ζ ω n σ | η | τ e √ 2 π ∞ Z 0  exp  − ( n − µ α ) 2 2 σ 2 | η |  + exp  − ( n + µ α ) 2 2 σ 2 | η |  ×  s  r − n ω d e − ζ ω n τ e  − s  r − n ω d  dn, ( 30 ) where τ e = − 1 ζ ω n log ρ c and s ( · ) denotes the step function. V elo cit y Env elop e Similarly , w e obtain the en v elop e distribution of the velocity of the system. The background dynamics distribution for v elo city w as obtained in equation ( 16 ). Noting that from ( 17 ), ˙ u r = ω d u r , the rare even t p df is modified by a constant factor f ˙ u ( r ) = f ˙ u b ( r )(1 − P r ) + ω − 1 d f u r ( r /ω d ) P r . ( 31 ) The final formula for the v elo city en velope pdf is giv en b y f ˙ u ( r ) = r σ 2 ˙ x b exp  − r 2 2 σ 2 ˙ x b  (1 − ν α τ e ) + ν α r ζ ω n σ | η | √ 2 π ∞ Z 0  exp  − ( n − µ α ) 2 2 σ 2 | η |  + exp  − ( n + µ α ) 2 2 σ 2 | η |  ×  s  r − n e − ζ ω n τ e  − s ( r − n )  dn. ( 32 ) A cceleration Env elop e Lastly we also obtain the en velope distribution of the acceleration. Noting that from eq. ( 17 ), ¨ u r = ω 2 d u r , the rare even t p df for acceleration is also modified b y a constan t factor f ¨ u ( r ) = f ¨ u b ( r )(1 − P r ) + ω − 2 d f u r ( r /ω 2 d ) P r . (33 ) The final formula for the acceleration en velope pdf is then f ¨ u ( r ) = r σ 2 ¨ x b exp  − r 2 2 σ 2 ¨ x b  (1 − ν α τ e ) + ν α r ζ ω n σ | η | √ 2 π ∞ Z 0  exp  − ( n − µ α ) 2 2 σ 2 | η |  + exp  − ( n + µ α ) 2 2 σ 2 | η |  ×  s  r − nω d e − ζ ω n τ e  − s ( r − nω d )  dn. ( 34 ) 9 3.2.5 Comparison with Monte-Carlo sim ulations F or the Mon te-Carlo sim ulations the excitation time series is generated b y sup erimp osing the bac kground and rare even t comp onents. The background excitation, described by a stationary sto c hastic pro cess with a Pierson-Mosk owitz sp ectrum ( equation ( 11 ) ), is simulated through a sup erp osition of cosines ov er a range of frequencies with corresp onding amplitudes and uniformly distributed random phases. The in termittent component is the random impulse train, and each impact is introduced as a velocity jump with a given magnitude at the p oint of the impulse impact. F or each of the comparisons p erformed in this work w e generate 10 realizations of the excitation time series, eac h with a train of 100 impulses. Once eac h ensem ble time series for the excitation is computed, the go verning ordinary differential equations are solv ed using a 4th/5th order Runge-Kutta metho d ( we carefully accoun t for the modifications in the momen tum that an impulse imparts b y in tegrating up to eac h impulse time and mo difying the initial conditions that the impulse imparts before in tegrating the system to the next impulse time ). F or each realization the system is integrated for a sufficien tly long time so that w e hav e conv erged resp onse statistics for the displacement, v elo city , and acceleration. W e utilize a shifted Pierson-Mosko witz sp ectrum S hh ( ω − 1) in order to av oid resonance. The other parameters and resulted statistical quantities of the system are giv en in table 1 . As it can b e seen in figure 3 the analytical approximations compare well with the Monte-Carlo sim ulations many standard deviations aw ay from the zero mean. The results are robust to differen t parameters as far as we satisfy the assumption of indep endent ( non- ov erlapping ) random even ts. Some discrepancies shown b etw een Monte-Carlo sim ulations can b e attributed to the en- v elop e approximation used for the rare even t quantification. Indeed, these discrepancies are reduced significan tly if one utilizes the semi-analytical metho d presente d in the next section, where we do not mak e any simplifications for the form of the resp onse during extreme ev en ts. T able 1: Parameters and relev ant statistical quan tities for SDOF system 1. λ 0 . 01 k 1 T α 5000 ζ 0 . 005 ω n 1 ω d 1 µ α = 7 × σ ˙ h 0 . 1 q 1 . 582 × 10 − 4 σ α = σ ˙ h 0 . 0143 σ h 0 . 0063 σ ˙ x b 0 . 0179 σ x b 0 . 0082 σ | η | 0 . 0229 P r 0 . 0647 10 u 0 0.05 0.1 0.15 Disp Env pdf 10 -3 10 -2 10 -1 10 0 10 1 10 2 du 0 0.05 0.1 0.15 Vel Env pdf 10 -3 10 -2 10 -1 10 0 10 1 10 2 ddu 0 0.05 0.1 0.15 0.2 0.25 Acc Env pdf 10 -3 10 -2 10 -1 10 0 10 1 10 2 Numerical Analytical u u u V elocity Env . PDF Acceleration Env . PDF Displacement Env . PDF Figure 3: [Sev erely underdamp ed case] Comparison b etw een direct Mon te-Carlo sim ulation and the analytical p df for the SDOF system 1. The pdf for the env elop e of each of the sto chastic v ariables, displacement, v elo city , and acceleration, are presen ted. The dashed line indicates one standard deviation. 3.3 Analytical p df for the o v erdamp ed case ζ  1 In the previous section, w e illustrated the deriv ation of the analytical resp onse p df under the assumption ζ  1 . Here, w e briefly summarize the results for the response pdf for the case where ζ  1 . One can follow the same steps using the corresp onding formulas for the rare ev en t transitions in the presence of large damping (app endix A ). A n imp ortant differ enc e for the over damp e d c ase is that the system do es not exhibit highly oscil latory motion as opp ose d to the under damp e d c ase, and henc e we dir e ctly work on the r esp onse p df inste ad of the envelop e p df . Displacemen t The total probabilit y la w b ecomes f x ( r ) = f x b ( r )(1 − P r, dis ) + f x r ( r ) P r, dis , ( 35 ) and we obtain the follo wing p df for the displacemen t of the system f x ( r ) = 1 σ x b √ 2 π exp  − r 2 2 σ 2 x b  (1 − ν α τ e, dis ) + ν α τ e, dis r ( ζ ω n − ω o ) σ η √ 2 π ( τ e, dis − τ s ) ∞ Z −∞ exp  − ( n − µ α ) 2 2 σ 2 η  ×  s  r − n 2 ω o e − ( ζ ω n − ω o ) τ e, dis  − s  r − n 2 ω o  dn, ( 36 ) where τ e, dis = π 2 ω o − 1 ζ ω n − ω o log ρ c . V elo cit y Similarly we deriv e the total probability la w for the response velocity f ˙ x ( r ) = f ˙ x b ( r )(1 − P r, v el ) + f ˙ x r ( r ) P r, v el . ( 37 ) 11 The final result for the velocity p df is f ˙ x ( r ) = 1 σ ˙ x b √ 2 π exp  − r 2 2 σ 2 ˙ x b  (1 − ν α τ e, vel ) + ν α τ e, vel r ( ζ ω n + ω o ) σ η √ 2 π τ e, vel ∞ Z −∞ exp  − ( n − µ α ) 2 2 σ 2 η  ×  s  r − n e − ( ζ ω n + ω o ) τ e, vel  − s ( r − n )  dn, ( 38 ) where τ e, vel = − 1 ζ ω n + ω o log ρ c . A cceleration The total probabilit y la w for the response acceleration is f ¨ x ( r ) = f ¨ x b ( r )(1 − P r, acc ) + f ¨ x r ( r ) P r, acc , ( 39 ) and this gives the follo wing result for the acceleration p df f ¨ x ( r ) = 1 σ ¨ x b √ 2 π exp  − r 2 2 σ 2 ¨ x b  (1 − ν α τ e, acc ) + ν α τ e, acc r ( ζ ω n + ω o ) σ η √ 2 π τ e, acc ∞ Z −∞ exp  − ( n − µ α ) 2 2 σ 2 η  ×  s  r − n ( ζ ω n + ω o )e − ( ζ ω n + ω o ) τ e, acc  − s ( r − n ( ζ ω n + ω o ))  dn, ( 40 ) where τ e, acc = − 1 ζ ω n + ω o log ρ c . Note that in this case w e do not ha ve the simple scaling as in the underdamp ed case for the conditionally rare p df. 3.3.1 Comparison with Monte-Carlo sim ulations W e confirm the accuracy of the analytical results given in equations ( 36 ), ( 38 ), and ( 40 ) for the strongly ov erdamp ed case through comparison with direct Mon te-Carlo simulations. The parameters and resulted statistical quantities of the system are giv en in table 2 . The analytical estimates show fa vorable agreemen t with n umerical simulations for this case ( figure 4 ), just as in the previous underdamp ed case. T able 2: Parameters and relev ant statistical quan tities for SDOF system 2. λ 6 k 1 T α 1000 ζ 3 ω n 1 ω d 2 . 828 µ α = 7 × σ ˙ h 0 . 1 q 1 . 582 × 10 − 4 σ α = σ ˙ h 0 . 0143 σ h 0 . 0063 σ ˙ x b 0 . 0056 σ x b 0 . 0022 σ η 0 . 0154 P r,dis 0 . 0140 P r,vel 0 . 0004 P r,acc 0 . 0004 12 ddx -0.2 -0.1 0 0.1 0.2 Acceleration PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Numerical Analytical dx -0.1 -0.05 0 0.05 0.1 Velocity PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 x -0.03 -0.02 -0.01 0 0.01 0.02 0.03 Displacement PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 x x x V elocity PDF Acceleration PDF Displacement PDF Figure 4: [Sev erely ov erdamp ed case] Comparison b etw een direct Monte-Carlo sim ulation and the analytical p df for SDOF system 2. The p df for the v alue of each stochastic process is sho wn. The dashed line indicates one standard deviation. 4 Semi-analytical p df of the resp onse of SDOF systems W e now form ulate a semi-analytical approach to quantify the resp onse p df for an y arbitrary set of system parameters, including the sev erely underdamp ed or ov erdamp ed cases consid- ered previously . The approach here adapts the n umerical sc heme describ ed in [ 21 ] for systems undergoing internal instabilities. While for the limiting cases that we studied previously knowledge of the trajectory ( time series ) of the system ( x r ( t ) or u r ( t ) ) could analytically be translated to information ab out the corresp onding pdf ( f x r or f u r ), this is not alw ays p ossible. In addition, for nonline ar structural systems one will not, in general, hav e analytical expressions for the rare ev ent transitions. F or these cases w e can compute the rare even t statistics b y numerically approximating the corresp onding histogram, using either analytical or numerically generated tra jectories for the rare even t regime. 4.1 Numerical computation of rare ev en ts statistics Consider the same SDOF system introduced in section 3 . Recall that we hav e quantified the resp onse p df by the PDS metho d using the total probability la w f x ( r ) = f x b ( r )(1 − P r ) + f x r ( r ) P r . ( 41 ) In the previous section, the deriv ation consisted of estimating all three unknown quan tities: the bac kground distribution f x b , the rare even t distribution f x r , and the rare even t probability P r analytically . Ho wev er, in the semi-analytical sc heme w e will obtain the rare even t distribution f x r and rare even t probabilit y P r b y taking a histogram of the n umerically sim ulated analytical form of the rare resp onse. The background distribution f x b is still obtained analytically as in section 3.1 . Recall that the rare even t distribution is giv en by f x r ( r ) = Z f x r | η ( r | n ) f η ( n ) dn, ( 42 ) where f η ( n ) is kno wn analytically ( equation ( 18 ) ). It is the conditional pdf f x r | η ( r | n ) that w e estimate by a histogram: f x r | η ( r | n ) = Hist  x r | η ( t | n )  , t = [ 0 , τ e, dis ] , ( 43 ) 13 where we use the analytical solution of the oscillator with non-zero initial velocity , n : x r | η ( t | n ) = n 2 ω o  e − ( ζ ω n − ω o ) t − e − ( ζ ω n + ω o ) t  . ( 44 ) The histogram is taken from t = 0 ( the beginning of the rare ev ent) un til the end of the rare even t at t = τ e . The conditional distribution of rare ev en t resp onse for the velocity and acceleration can b e written as well: f ˙ x r | η ( r | n ) = Hist  ˙ x r | η ( t | n )  , t = [ 0 , τ e, vel ] , ( 45 ) f ¨ x r | η ( r | n ) = Hist  ¨ x r | η ( t | n )  , t = [ 0 , τ e, acc ] . ( 46 ) 4.2 Numerical estimation of the rare even ts probabilit y In order to compute the histogram of a rare impulse ev ent, the duration of a rare resp onse needs to b e obtained numerically . Recall that w e ha ve defined the duration of a rare resp onses b y x r ( τ e ) = ρ c max  | x r |  , ( 47 ) where ρ c = 0 . 1 . In the n umerical computation of τ e the absolute maxim um of the resp onse needs to b e estimated numerically as well. Once the rare ev ent duration has b een sp ecified, w e can obtain the probability of rare ev ents b y P r = ν α τ e = τ e /T α . ( 48 ) This v alue is indep endent of the conditional background magnitude. The ab ov e pro cedure is applied for the rare even t resp onse displacemen t τ e, dis , velocity τ e, vel , and acceleration τ e, acc . 4.3 Semi-analytical probability density functions W e can now compute the resp onse pdf using the describ ed semi-analytical approach. F or the displacemen t w e hav e: f x ( r ) = 1 − ν α τ e, dis σ x b √ 2 π exp  − r 2 2 σ 2 x b  + ν α τ e, dis ∞ Z 0 Hist  x r | η ( t | n )  f η ( n ) dn. ( 49 ) The corresp onding p df for the velocity f ˙ x and acceleration f ¨ x , can b e computed with the same form ula but with the appropriate v ariance for the Gaussian core ( σ ˙ x b or σ ¨ x b ), rare even t duration ( τ e, vel or τ e, acc ), and histograms ( ˙ x r | η ( t | n ) or ¨ x r | η ( t | n ) ). 4.3.1 Comparison with Mon te-Carlo simulations F or illustration, a SDOF configuration is considered with critical damping ratio, ζ = 1 . The detailed parameters and relev ant statistical quan tities of the system are given in table 3 . This is a regime where the analytical results deriv ed in section 3 are not applicable. Even for this ζ v alue, the semi-analytical p df for the resp onse sho ws excellent agreement with direct simulations ( figure 5 ). W e emphasize that the computational cost of the semi-analytical scheme is comparable with that of the analytical approximations (order of seconds) and both are significan tly low er than the cost of Monte-Carlo sim ulation ( order of hours ). 14 T able 3: Parameters and relev ant statistical quan tities for SDOF system 3. λ 2 k 1 T α 400 ζ 1 ω n 1 ω d 0 µ α = 7 × σ ˙ h 0 . 1 q 1 . 582 × 10 − 4 σ α = σ ˙ h 0 . 0143 σ h 0 . 0063 σ ˙ x b 0 . 0120 σ x b 0 . 0052 σ α = σ η 0 . 0187 P r,dis 0 . 0122 P r,vel 0 . 0075 P r,acc 0 . 0032 ddx -0.3 -0.2 -0.1 0 0.1 0.2 0.3 Acceleration PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 Numerical Semi-Analytical dx -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Velocity PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 x -0.06 -0.04 -0.02 0 0.02 0.04 0.06 Displacement PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 x x x V elocity PDF Acceleration PDF Displacement PDF Figure 5: [Critically damp ed system] Comparison b etw een direct Monte-Carlo sim ulations and the semi-analytical p df for SDOF system 3. Dashed lines indicate one standard deviation. 5 Semi-analytical p df for the resp onse of MDOF systems An imp ortant adv antage of the semi-analytical scheme is the straightforw ard applicability of the algorithm to MDOF systems. In this section w e demonstrate how the extension can b e made for a tw o-degree-of-freedom ( TDOF ) linear system ( see figure 6 ): m ¨ x + λ ˙ x + k x + λ a ( ˙ x − ˙ y ) + k a ( x − y ) = F ( t ) , ( 50 ) m a ¨ y + λ a ( ˙ y − ˙ x ) + k a ( y − x ) = 0 , ( 51 ) where the sto chastic forcing F ( t ) = F b ( t ) + F r ( t ) is applied to the first mass ( mass m ) and x, y are displacements of the tw o masses. As b efore, F b ( t ) is the background comp onent and F r ( t ) = P N ( t ) i =1 α i δ ( t − τ i ) is the rare even t comp onent. k λ m k a λ a m a F=F b +F r z 1 z 2 Figure 6: The considered TDOF system. The excitation is applied to the first mass, with mass m . 15 The background statistics are obtained by analyzing the resp onse spectrum of the TDOF system sub jected to the background excitation component. The details are giv en in App endix C . F or a nonlinear system this can be done using statistical linearization metho d [ 27 ]. The his- tograms for the rare ev ent transitions can b e computed through standard analytical expressions that one can derive for a linear system like the one we consider under the set of initial conditions: x 0 = y 0 = ˙ y 0 = 0 and ˙ x 0 = n. Once the impulse resp onse has b een obtained we n umerically quantify the rare even t distri- bution, as well as the rare ev ent duration and use the semi-analytical decomp osition. The pdf is then given b y f z ( r ) = 1 − ν α τ z e σ z b √ 2 π exp  − r 2 2 σ 2 z b  + ν α τ z e ∞ Z 0 Hist  z r | η ( t | n )  f η ( n ) dn, ( 52 ) where z can b e either of the degrees-of-freedom ( x or y ) or the corresp onding velocities or accelerations, while τ z e is the typical duration of the rare even ts and is estimated numerically . W e note that as in the previous cases the p df is comp osed of a Gaussian core ( describing the bac kground statistics ) as well as, a heavy tailed comp onent that is connected with the rare transitions. F or each case of z the corresp onding v ariance under bac kground excitation, temp oral durations of rare even ts, and histograms for rare ev ents should b e employ ed. Results are presented for the pdf of the displacement, velocity and acceleration of each degree-of-freedom ( figure 7 ). These compare fa vorably with the direct Mon te-Carlo simulations. The parameters and resulted statistical quantities of the system are given in table 4 . F urther n umerical sim ulations ( not presented) demonstrated strong robustness of the approach. T able 4: Parameters and relev ant statistical quan tities for the TDOF system. m 1 m a 1 λ 0 . 01 k 1 λ a 1 k a 0 . 1 T α 1000 σ η 0 . 0199 µ α 0 . 1 q 1 . 582 × 10 − 4 σ α 0 . 0143 σ F b 0 . 0351 P x r,dis 0 . 0177 P y r,dis 0 . 0190 P x r,vel 0 . 0098 P y r,vel 0 . 0209 P x r,acc 0 . 0066 P y r,acc 0 . 0082 16 v -0.1 -0.05 0 0.05 0.1 Diplacement PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 ddx -0.2 -0.1 0 0.1 0.2 Acceleration PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 ddv -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Acceleration PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Numerical Semi-Analytical dx -0.2 -0.1 0 0.1 0.2 Velocity PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 dv -0.1 -0.05 0 0.05 0.1 Velocity PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 x -0.15 -0.1 -0.05 0 0.05 0.1 Diplacement PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 y x V elocity PDF Displacement PDF x Displacement PDF x Acceleration PDF y V elocity PDF y Acceleration PDF Figure 7: [T wo DOF System] Comparison b etw een direct Monte-Carlo simulation and the semi-analytical appro ximation. The p df for the v alue of the time series are presented. Dashed line indicates one standard deviation. 6 Semi-analytical p df of lo cal maxima It is straigh tforward to extend the semi-analytical framework for other quan tities of interest, such as the lo cal extrema ( maxima and minima ) of the resp onse. The numerically generated histogram for the rare transitions can b e directly computed for the lo cal extrema as w ell. On the other hand, for the bac kground excitation regime, we use known results from the theory of stationary Gaussian sto chastic processes to describ e the corresponding bac kground p df analytically . 6.1 Distribution of lo cal maxima under bac kground excitation F or a stationary Gaussian pro cess with arbitrary sp ectral bandwidth  , the probability densit y function of p ositive extrema (maxima ) is given b y [ 24 , 13 ]: f m + ( ξ ) =  √ 2 π e − ξ 2 / 2  2 + p 1 −  2 ξ e − ξ 2 / 2 Φ  √ 1 −  2  ξ  , −∞ ≤ ξ ≤ ∞ , ( 53 ) where ξ = x √ µ 0 , x is the magnitude of the maxima, sp ectral bandwidth  = q 1 − µ 2 2 µ 0 µ 4 , and Φ( x ) = 1 √ 2 π R x −∞ e − u 2 / 2 du is the standard normal cumulativ e distribution function. The sp ectral momen ts for the background response displacement x b are also defined as µ n = ∞ Z 0 ω n S x b x b ( ω ) dω . ( 54 ) 17 W e note that that for the limit of an infinitesimal narrow-banded signal (  = 0) , the pdf conv erges to a Ra yleigh distribution. On the other hand for an infinitely broad-banded signal (  = 1) , the distribution con v erges to the Gaussian p df. F or a signal with in-b et ween sp ectral bandwidth (0 ≤  ≤ 1) , the pdf has a blended structure with the form in equation ( 53 ). T aking in to account the asymmetric nature of the in termittent excitation, w e need to consider b oth the p ositive and negativ e extrema. F or the background excitation the pdf can be written as: f ˆ x b ( x ) = 1 2 √ µ 0 n f m +  x µ 0  + f m +  − x µ 0 o , −∞ ≤ x ≤ ∞ , ( 55 ) where the ˆ x notation denotes the lo cal extrema of x . Similar expressions can be obtained for the velocity and acceleration extrema. 6.2 Statistics of lo cal extrema during rare transitions The conditional pdf f ˆ x r | η for the local extrema can be n umerically estimated through the histogram: f ˆ x r | η ( r | n ) = Hist  M  x r | η ( t | n )  , t = [ 0 , τ e, dis ] , ( 56 ) where M ( · ) is an op erator that giv es all the p ositive/negativ e extrema. The p ositive/negativ e extrema are defined as the p oints where the deriv ative of the signal is zero. 6.3 Semi-analytical probability density function for lo cal extrema The last step consists of applying the decomposition of the p df. This takes the form: f ˆ x ( r ) = (1 − ν α τ e, dis ) f ˆ x b ( r ) + ν α τ e, dis ∞ Z 0 Hist  M  x r | η ( t | n )  f η ( n ) dn. ( 57 ) The same decomp osition can b e used for the velocity and acceleration lo cal maxima. W e compare the semi-analytical decomp osition with Mon te-Carlo sim ulations. In figure 8 we presen t results for the tw o-degree-of-freedom system. The p df are shown for lo cal extrema of the displacemen ts, v elo cities and accelerations for each degree of freedom. W e emphasize the non-trivial structure of the p df and especially their tails. Throughout these comparisons the semi-analytical scheme demonstrates accurate estimation of b oth the heavy tails and the non-Gaussian/non-Ra yleigh structure of the background local extrema distribution. 18 ddx -0.3 -0.2 -0.1 0 0.1 0.2 0.3 A cceleration Extrema PD F 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 ddv -0.2 -0.1 0 0.1 0.2 A cceleration Extrema PD F 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Numerical Semi-Analytical dx -0.2 -0.1 0 0.1 0.2 Velocity Extrema PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 dv -0.06 -0.04 -0.02 0 0.02 0.04 0.06 Velocity Extrema PDF 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 x -0.1 -0.05 0 0.05 0.1 D isplacement Extrema PD F 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 v -0.1 -0.05 0 0.05 0.1 D isplacement Extrema PD F 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 x V elocity Extrema PDF Displacement Extrema PDF Displacement Extrema PDF x Acceleration Extrema PDF V elocity Extrema PDF Acceleration Extrema PDF y ^ y ^ y ^ x ^ ^ ^ Figure 8: [Lo cal extrema for TDOF] Comparison b etw een direct Mon te-Carlo simulation and the semi-analytical approximation. The pdf for the lo cal extrema of the resp onse are presented. Dashed line indicates one standard deviation. 7 Summary and conclusions W e hav e formulated a robust appro ximation metho d to quan tify the probabilistic response of structural systems sub jected to sto chastic excitation containing intermitten t comp onents. The foundation of our approach is the recen tly developed probabilistic decomp osition-synthesis metho d for the quantification of rare ev ents due to in ternal instabilities to the problem where extreme resp onses are triggered by external forcing. The in termittent forcing is represen ted as a background component, modeled through a colored processes with energy distributed across a range of frequencies, and additionally a rare/extreme comp onent that can b e represen ted b y impulses that are Poisson distributed with large inter-arriv al time. Owing to the nature of the forcing, even the probabilistic resp onse of a linear system can b e highly complex with asymmetry and complicated tail b ehavior that is far from Gaussian, whic h is the exp ected form of the resp onse p df if the forcing did not contain an in termittently extreme component. The main result of this work is the deriv ation of analytical/semi-analytical expressions for the p df of the resp onse and its lo cal extrema for structural systems ( including the resp onse displacemen t, velocity , and acceleration p df ). These expressions decompose the p df in to a probabilistic core, capturing the statistics under background excitation, as well as a heavy- tailed comp onent asso ciated with the extreme transitions due to the rare impacts. W e hav e p erformed a thorough analysis for linear SDOF systems under v arious system parameters and also derived analytical form ulas for tw o special cases of parameters ( lightly damped or hea vily damp ed systems). The general semi-analytical decomposition is applicable for any arbitrary set of system parameters and w e hav e demonstrated its v alidit y through comprehensive comparisons with Mon te- Carlo simulations. The general framework is also directly applicable to multi-degree- 19 of-freedom MDOF systems, as w ell as systems with nonlinearities, and w e hav e assessed its p erformance through a 2DOF linear system of tw o coupled oscillators excited through the first mass. Mo difications of the method to compute statistics of lo cal extrema ha ve also b een presen ted. The developed approac h allo ws for computation of the resp onse pdf of structural systems man y orders of magnitude faster than a direct Monte-Carlo sim ulation, which is curren tly the only reliable to ol for suc h computations. The rapid ev aluation of resp onse pdfs for systems excited b y extreme forcing by the metho d presented in this work pav es the wa y for enabling robust design of structural systems subjected to extreme even ts of a sto c hastic nature. Our future endeav ors include the application of the developed framework to the optimization of engineering systems where extreme even t mitigation is required. In such cases, it is usually not feasible to run Monte-Carlo simulations for v arious parameter sets during the design pro cess o wing to the computational costs associated with low probability rare ev ents, let alone p erform parametric optimization. W e believe our approac h is w ell suited to suc h problems and can pro ve to b e an imp ortant method for engineering design and reliability assessmen t. A c kno wledgmen ts T.P .S. has b een supp orted through the ONR grants N00014-14-1-0520 and N00014-15-1-2381 and the AF OSR gran t F A9550-16-1-0231. H.K.J. and M.A.M. ha ve b een supp orted through the first and third gran ts as graduate students. W e are also grateful to the Samsung Scholarship Program for support of H.K.J. as w ell as the MIT Energy Initiativ e for supp ort under the grant ‘Nonlinear Energy Harvesting F rom Broad-Band Vibrational Sources By Mimic king T urbulent Energy T ransfer Mechanisms’ . A Impulse resp onse of SDOF systems The resp onse of the system ¨ x r ( t ) + λ ˙ x r ( t ) + k x r ( t ) = 0 , ( 58 ) under an impulse α at an arbitrary time t 0 , say t 0 = 0 , and with zero initial v alue but nonzero initial v elo city ( ( x r , ˙ x r ) = (0 , ˙ x r 0 ) at t = 0 − ) are given by the follo wing equations under the t w o limiting cases of damping: Sev erely underdamp ed case ζ  1 With the approximation of ζ  1 (or ω d ≈ ω n ), we can simplify resp onses as x r ( t ) = α + ˙ x r 0 ω d e − ζ ω n t sin ω d t, ( 59 ) ˙ x r ( t ) = ( α + ˙ x r 0 )e − ζ ω n t cos ω d t, ( 60 ) ¨ x r ( t ) = − ( α + ˙ x r 0 ) ω d e − ζ ω n t sin ω d t. ( 61 ) Sev erely o verdamped case ζ  1 Similarly , with the appro ximation of ζ  1 ( or ω o ≈ ζ ω n ), we can simplify responses as x r ( t ) = ( α + ˙ x r 0 ) 2 ω o e − ( ζ ω n − ω o ) t , (62 ) ˙ x r ( t ) = ( α + ˙ x r 0 )e − ( ζ ω n + ω o ) t , ( 63 ) ¨ x r ( t ) = − ( ζ ω n + ω o )( α + ˙ x r 0 )e − ( ζ ω n + ω o ) t . ( 64 ) 20 B Probabilit y distribution for an arbitrarily exp onen tially deca ying function Consider an arbitrary time series in the follo wing form: x ( t ) = A e − αt , where t ∼ Uniform ( τ 1 , τ 2 ) , ( 65 ) where A , α > 0 and τ 1 < τ 2 . The cumulativ e distribution function ( cdf ) of x ( t ) is F x ( x ) = P ( A e − αt < x ) , ( 66 ) = P  t > 1 α log( A /x )  , ( 67 ) = 1 − P  t < 1 α log( A /x )  , ( 68 ) = 1 − 1 α log( A /x ) Z −∞ f T ( t ) dt, ( 69 ) where f T ( t ) is expressed using the step function s ( · ) as: f T ( t ) = 1 τ 2 − τ 1  s ( t − τ 1 ) − s ( t − τ 2 )  , τ 1 < τ 2 . ( 70 ) The p df of the resp onse x ( t ) can then b e derived b y differentiation. f x ( x ) = d dx F x ( x ) , ( 71 ) = 1 αx f T  1 α log( A x )  , ( 72 ) = 1 αx ( τ 2 − τ 1 ) n s ( x − A e − ατ 2 ) − s ( x − A e − ατ 1 ) o . ( 73 ) W e utilize the ab ov e formula for deriving analytical pdfs. Note that the step function with resp ect to x in the ab o ve has b een derived using τ 1 < t < τ 2 , ( 74 ) − ατ 2 < − αt < − ατ 1 , ( 75 ) A e − ατ 2 < x < A e − ατ 1 . (76 ) C Bac kground resp onse for TDOF system Consider the statistical resp onse of the system to the background forcing component, m ¨ x b + λ ˙ x b + k x b + λ a ( ˙ x b − ˙ y b ) + k a ( x b − y b ) = F b ( t ) , ( 77 ) m a ¨ y b + λ a ( ˙ y b − ˙ x b ) + k a ( y b − x b ) = 0 . 21 The sp ectral densities are given b y S x b x b ( ω ) = ω 4 S F b ( ω )  A ( ω ) − B ( ω ) 2 C ( ω )  A ( − ω ) − B ( − ω ) 2 C ( − ω )  , ( 78 ) S ˙ x b ˙ x b ( ω ) = ω 2 S x b x b ( ω ) , ( 79 ) S ¨ x b ¨ x b ( ω ) = ω 4 S x b x b ( ω ) , ( 80 ) S y b y b ( ω ) = ω 4 S F b ( ω )  A ( ω ) C ( ω ) B ( ω ) − B ( ω )  A ( − ω ) C ( − ω ) B ( − ω ) − B ( − ω )  , (81 ) S ˙ y b ˙ y b ( ω ) = ω 2 S y b y b ( ω ) , ( 82 ) S ¨ y b ¨ y b ( ω ) = ω 4 S y b y b ( ω ) , ( 83 ) where S F b ( ω ) is the sp ectral densit y of F b ( t ) , and A ( ω ) =( λ a + λ )( j ω ) + ( k a + k ) − mω 2 , ( 84 ) B ( ω ) = λ a ( j ω ) + k a , ( 85 ) C ( ω ) = λ a ( j ω ) + k a − m a ω 2 . ( 86 ) Th us, w e can obtain the following conditionally bac kground v ariances σ 2 x b = ∞ Z 0 S x b x b ( ω ) dω , σ 2 ˙ x b = ∞ Z 0 S ˙ x b ˙ x b ( ω ) dω , σ 2 ¨ x b = ∞ Z 0 S ¨ x b ¨ x b ( ω ) dω , ( 87 ) σ 2 y b = ∞ Z 0 S v b v b ( ω ) dω , σ 2 ˙ y b = ∞ Z 0 S ˙ y b ˙ y b ( ω ) dω , σ 2 ¨ y b = ∞ Z 0 S ¨ y b ¨ y b ( ω ) dω . ( 88 ) References [1] A. M. Ab ou-Ra y an and A. H. Nayfeh. “Sto chastic resp onse of a buc kled beam to external and parametric random excitations”. In: In: AIAA/ASME/ASCE/AHS/ASC Structur es ( 1993 ), pp. 1030–1040. [2] G. A. Athanassoulis and P . N. Gavriliadis. “The truncated Hausdorff moment problem solv ed b y using kernel density functions”. In: Pr ob abilistic Engine ering Me chanics 17.3 ( July 2002 ), pp. 273–291. [3] G. Athanassoulis, I.C. T santili, and Z.G. Kap elonis. “Bey ond the Mark ovian assumption: Resp onse-excitation probabilistic solution to random nonlinear differential equations in the long time”. In: Pr o c e e dings of the R oyal So ciety A 471 ( 2016 ), p. 20150501. [4] G. Barone, G. Na v arra, and A. Pirrotta. “Probabilistic resp onse of linear structures equipp ed with nonlinear damp er devices ( PIS metho d )”. In: Pr ob abilistic engine ering me chanics 23.2 ( 2008 ), pp. 125–133. [5] V. L. Belenky and N. B. Sev astianov. Stability and Safety of Ships: R isk of Capsizing . The So ciet y of Nav al Architects and Marine Engineers, 2007. [6] M Beran. Statistic al Continuum The ories . Interscience Publishers, 1968. [7] Linda J. Branstetter, Garrett D. Jeong, James T.P . Y ao, Y.K. W en, and Y.K. Lin. “Math- ematical mo delling of structural b eha viour during earthquakes”. In: Pr ob abilistic Engi- ne ering Me chanics 3.3 ( Sept. 1988 ), pp. 130–145. 22 [8] W. Cousins and T. P . Sapsis. “Quantification and prediction of extreme even ts in a one-dimensional nonlinear disp ersive w av e mo del”. In: Physic a D 280 ( 2014 ), pp. 48–58. [9] W. Cousins and T. P . Sapsis. “Reduced order precursors of rare even ts in unidirectional nonlinear water w av es”. In: Journal of Fluid Me chanics 790 ( 2016 ), pp. 368–388. [10] A. Di Matteo, M. Di Paola, and A. Pirrotta. “Probabilistic c haracterization of nonlin- ear systems under Poisson white noise via complex fractional moments”. In: Nonline ar Dynamics 77.3 ( 2014 ), pp. 729–738. [11] R. Iw ankiewicz and S. R. K. Nielsen. “Solution techniques for pulse problems in non-linear sto c hastic dynamics”. In: Pr ob abilistic engine ering me chanics 15.1 ( 2000 ), pp. 25–36. [12] H. K. Jo o and T. P . Sapsis. “A momen t-equation-copula-closure metho d for nonlinear vibrational systems sub jected to correlated noise”. In: Pr ob abilistic Engine ering Me chanics ( 2016 ). [13] H. Karadeniz. Sto chastic analysis of offshor e ste el structur es: an analytic al appr aisal . Springer Science & Business Media, 2012. [14] H. U. K ö ylüoglu, S. R. K. Nielsen, and R. Iwankiewicz. “Response and reliabilit y of P oisson-driv en systems b y path in tegration”. In: Journal of engine ering me chanics 121.1 ( 1995 ), pp. 117–130. [15] R S Langley. “On v arious definitions of the en velope of a random pro cess.” In: Journal of Sound and Vibr ation 105.3 ( 1986 ), pp. 503–512. [16] Y. K. Lin. “Application of nonstationary shot noise in the study of system resp onse to a class of nonstationary excitations”. In: Journal of A pplie d Me chanics 30.4 ( 1963 ), pp. 555– 558. [17] Y. K. Lin and C. Q. Cai. Pr ob abilistic Structur al Dynamics . McGraw-Hill, Inc., 1995. [18] Y.K. Lin. “Sto chastic stabilit y of wind-excited long-span bridges”. In: Pr ob abilistic Engi- ne ering Me chanics 11.4 ( Oct. 1996 ), pp. 257–261. [19] A. J. Majda and M. Branic ki. “Lessons in Uncertaint y Quan tification for T urbulent Dynamical Systems”. In: Discr ete and Continuous Dynamic al Systems 32 ( 2012 ), pp. 3133– 3221. [20] A. Masud and L. A. Bergman. “Solution of the four dimensional F okk er-Planck Equation: Still a challenge”. In: ICOSSAR 2005 (2005 ), pp. 1911–1916. [21] M. A. Mohamad, W. Cousins, and T. P . Sapsis. “A probabilistic decomp osition-synthesis metho d for the quantification of rare ev ents due to in ternal instabilities”. In: Journal of Computational Physics 322 ( 2016 ), pp. 288–308. [22] M. A. Mohamad and T. P . Sapsis. “Probabilistic description of extreme even ts in in ter- mitten tly unstable systems excited by correlated sto c hastic pro cesses”. In: SIAM ASA J. of Unc ertainty Quantific ation 3 ( 2015 ), pp. 709–736. [23] M. A. Mohamad and T. P . Sapsis. “Probabilistic resp onse and rare even ts in Mathieu’s equation under correlated parametric excitation”. In: Oc ean Engine ering Journal 120 ( 2016 ), pp. 289–297. [24] K. Ochi M. A pplie d pr ob ability and sto chastic pr o c esses: In Engine ering and Physic al Scienc es . V ol. 226. Wiley-Interscience, 1990. [25] M. R. Riley and T. W. Coats. “A Simplified Approach for Analyzing Accelerations In- duced b y W a ve- Impacts in High-Speed Planing Craft”. In: 3r d Chesap e ake Power Bo at Symp osium June ( 2012 ), pp. 14–15. 23 [26] M. R. Riley, T. Coats, K. Haupt, and D. Jacobson. “Ride Sev erity Index – A New Approach to Quantifying the Comparison of A cceleration Resp onses of High-Sp eed Craft”. In: F AST 2011 11th International Confer enc e on F ast Se a T r ansp ortation Septem b er ( 2011 ), pp. 693– 699. [27] J. Rob erts and P . Spanos. R andom V ibr ation and Statistic al Line arization . Dov er Publi- cations, 2003. [28] T. P . Sapsis and G. A. A thanassoulis. “New partial differential equations go v erning the join t, resp onse-excitation, probability distributions of nonlinear systems, under general sto c hastic excitation”. In: Pr ob ab. Eng. Me ch. 23.2-3 ( 2008 ), pp. 289–306. [29] K. Sob czyk. Sto chastic differ ential e quations: with applic ations to physics and engine ering . V ol. 40. Springer, 2001. [30] T. T. So ong and M. Grigoriu. “Random vibration of mec hanical and structural systems”. In: NASA STI/R e c on T e chnic al R ep ort A 93 ( 1993 ), p. 14690. [31] S. Spence and M. Gioffrè. “Large scale reliability-based design optimization of wind excited tall buildings”. In: Pr ob abilistic Engine ering Me chanics 28 ( 2012 ), pp. 206–215. [32] Seshan Srirangara jan, Mic hael Allen, Ami Preis, Mudasser Iqbal, Hock Beng Lim, and Andrew J. Whittle. “W av elet-based Burst Ev en t Detection and Localization in W ater Distribution Systems”. In: Journal of Signal Pr o c essing Systems 72.1 ( July 2013 ), pp. 1– 16. [33] M. V asta. “Exact stationary solution for a class of non-linear systems driven by a non- normal delta-correlated pro cess”. In: International journal of non-line ar me chanics 30.4 ( 1995 ), pp. 407–418. [34] D. V en turi, T. P . Sapsis, H. Cho, and G. E. Karniadakis. “A computable evolution equation for the joint resp onse-excitation probability densit y function of sto chastic dynamical systems”. In: Pr o c e e dings of the R oyal So ciety A 468 ( 2012 ), p. 759. [35] W. F. W u and Y. K. Lin. “Cumulan t-neglect closure for non-linear oscillators under random parametric and external excitations”. In: International Journal of Non-line ar Me chanics 19.4 ( 1984 ), pp. 349–362. [36] D. Xiu and G. Karniadakis. “The Wiener-Ask ey p olynomial chaos for sto chastic differential equations”. In: SIAM Journal on Scientific Computing 24 ( 2002 ), pp. 619–644. [37] D. Xiu and G. Karniadakis. “Mo deling uncertaint y in flow sim ulations via generalized p olynomial chaos”. In: J. Comp. Phys. 187 ( 2003 ), pp. 137–167. [38] Y. Zeng and W. Q. Zhu. “Stochastic av eraging of strongly nonlinear oscillators under P oisson white noise excitation”. In: IUT AM symp osium on nonline ar sto chastic dynamics and c ontr ol . Springer. 2011, pp. 147–155. [39] W. Q. Zhu. “Stochastic a veraging methods in random vibration”. In: A pplie d Me chanics R eviews 41.5 ( 1988 ), pp. 189–199. 24