We characterize the complex, heavy-tailed probability density functions (pdfs) describing the response and its local extrema for structural systems subject to random forcing that includes extreme events. Our approach is based on recent probabilistic decomposition-synthesis (PDS) technique (Mohamad, M. A., Cousins, W., and Sapsis, T. P., 2016, “A Probabilistic Decomposition-Synthesis Method for the Quantification of Rare Events Due to Internal Instabilities,” J. Comput. Phys., 322, pp. 288–308), where we decouple rare event regimes from background fluctuations. The result of the analysis has the form of a semi-analytical approximation formula for the pdf of the response (displacement, velocity, and acceleration) and the pdf of the local extrema. For special limiting cases (lightly damped or heavily damped systems), our analysis provides fully analytical approximations. We also demonstrate how the method can be applied to high dimensional structural systems through a two-degrees-of-freedom (TDOF) example system undergoing extreme events due to intermittent forcing. The derived formulas can be evaluated with very small computational cost and are shown to accurately capture the complicated heavy-tailed and asymmetrical features in the probability distribution many standard deviations away from the mean, through comparisons with expensive Monte Carlo simulations.

## Introduction

A large class of physical systems in engineering and science can be modeled by stochastic differential equations. For many of these systems, the dominant source of uncertainty is due to the forcing which can be modeled by a stochastic process. Applications include ocean engineering systems excited by water waves (such as ship motions in large waves [1–4] or high speed crafts subjected to rough seas [5,6]) and rare events in structural systems (such as beam buckling [7,8], vibrations due to earthquakes [9,10], and wind loads [11,12]). For all of these cases, it is common that hidden in the otherwise predictable magnitude of the fluctuations are extreme events, i.e., abnormally large magnitude forces which lead to rare responses in the dynamics of the system (Fig. 1). Clearly, these events must be adequately taken into account for the effective quantification of the reliability properties of the system. In this work, we develop an efficient method to fully describe the probabilistic response of linear structural systems under general time-correlated random excitations containing rare and extreme events.

Systems with forcing having these characteristics pose significant challenges for traditional uncertainty quantification schemes. While there is a large class of methods that can accurately resolve the statistics associated with random excitations (e.g., the Fokker–Planck (FP) equation [13,14] for systems excited by white-noise and the joint response-excitation method [15–18] for arbitrary stochastic excitation), these have important limitations for high dimensional systems. In addition, even for low-dimensional systems determining the part of the probability density function (pdf) associated with extreme events poses important numerical challenges. On the other hand, Gaussian closure schemes and moment equation or cumulant closure methods [19,20] either cannot “see” rare events completely or they are very expensive and require the solution of an inverse moment problem in order to determine the pdf of interest [21]. Similarly, approaches relying on polynomial-chaos expansions [22,23] have been shown to have important limitations for systems with intermittent responses [24].

For system reliability, Monte Carlo and parameter estimation based methods are another class of strategies. It is now well established that anomalous diffusion leads to Lévy motion when particle jumps have power law tail probabilities (compared to Brownian motion, where the particle jumps have finite first two moments) for which fractional calculus is the appropriate tool for such problems [25]. Consequently, Lévy stable, Pareto, and Mittag–Leffler distributions have all been utilized in methods to understand and quantify non-Gaussian behavior in fractional-order mechanics [26–29], including fatigue modeling in solid mechanics [30], where it has been noted that the relatively simple description of Brownian diffusion is not adequate to capture observed non-Gaussian behavior.

Another popular approach for the study of rare event statistics in systems under intermittent forcing is to represent extreme events in the forcing as identically distributed independent impulses arriving at random times. The generalized FP equation or Kolmogorov–Feller equation is the governing equation that solves for the evolution of the response pdf under Poisson noise [14]. However, exact analytical solutions are available only for a limited number of special cases [31]. Although alternative methods such as the path integral method [32–34] and the stochastic averaging method [35,36] may be applied, solving the FP or Kolmogorov–Feller equations is often very expensive even for very low dimensional systems [37,38].

Here, we consider the quantification of the probability distribution of the response of linear systems subjected to stochastic forcing containing extreme events based on recently formulated probabilistic-decomposition synthesis (PDS) method [39,40]. The approach relies on the decomposition of the statistics into a “nonextreme core,” typically Gaussian, and a heavy-tailed component. This decomposition is in full correspondence with a partition of the phase space into a “stable” region where we do not have extreme events and a region where nonlinear instabilities or external forcing lead to rare transitions with high probability. We quantify the statistics in the stable region using a Gaussian approximation approach, while the non-Gaussian distribution associated with the intermittently unstable regions of phase space is performed taking into account the nontrivial character of the dynamics (either due to an instability or external forcing). The probabilistic information in the two domains is analytically synthesized through a total probability argument.

### Central Contributions.

We begin with the simplest case of a linear, single-degree-of-freedom (SDOF) system and then formulate the method for multidegree-of-freedom systems. *The main result of our work is the derivation of analytic/semi-analytic approximation formulas for the response pdf and the pdf of the local extrema of intermittently forced systems that can accurately characterize the statistics many standard deviations away from the mean*. In brief, the principal contributions of this paper are as follows:

Analytical (under certain conditions) and semi-analytical (under no restrictions) pdf expressions for the response displacement, velocity, and acceleration for single-degree-of-freedom systems under intermittent forcing.

Semi-analytical pdf expressions for the value and the local extrema of the displacement, velocity, and acceleration for multidegree-of-freedom systems under intermittent forcing.

The systems considered in this work are linear; however, the general method is directly applicable to nonlinear structural systems (see Ref. [41] for nonlinear systems).

The proposed approach circumvents the challenges that rare and extreme events pose for traditional uncertainty quantification schemes, in particular the computational burden associated when dealing with rare events in systems. We emphasize the statistical accuracy and the computational efficiency of the presented approach, which we rigorously demonstrate through extensive comparisons with direct Monte Carlo simulations.

The paper is structured as follows: In Sec. 2, we provide a general formulation of the probabilistic decomposition-synthesis method for structural systems under intermittent forcing. Next, in Sec. 3, we apply the developed method analytically, which is possible for two limiting cases: underdamped systems with $\zeta \u226a1$ or overdamped with $\zeta \u226b1$, where *ζ* is the damping ratio. The system we consider is excited by a forcing term consisting of a background time-correlated stochastic process superimposed with a random impulse train (describing the rare and extreme component). We give a detailed derivation of the response pdf of the system (displacement, velocity, and acceleration) and compare the results with expensive Monte Carlo simulations. In Sec. 4, we slightly modify the developed formulation to derive a semi-analytical scheme for the same linear system without restriction on the damping ratio *ζ*, demonstrating global applicability of the proposed approach. In Sec. 5, we consider multiple-degree-of-freedom systems, and in Sec. 6, include the response pdf results for the local extrema. Concluding remarks are made in Sec. 7.

## The Probabilistic Decomposition-Synthesis Method for Intermittently Forced Structural Systems

*σ*and a

_{b}*rare and extreme*component $Fr$ with magnitude $\sigma r\u226b\sigma b.$ The components $Fb$ and $Fr$ may both be (weakly) stationary stochastic processes, while the sum of the two processes is in general nonstationary. This can be seen if we directly consider the sum of two (weakly) stationary processes

*x*

_{1}and

*x*

_{2}, with time correlation functions $Corrx1(\tau )$ and $Corrx2(\tau )$, respectively. Then for the sum $z=x1+x2$, we have

Therefore, the process *z* is stationary if and only if the cross-covariance terms $E[x1(t)x2(t+\tau )]$ and $E[x1(t+\tau )x2(t)]$ are functions of *τ* only or they are zero (i.e., *x*_{1} and *x*_{2} are not correlated).

For the case where the excitation is given in terms of realizations, i.e., time-series, one can first separate the extreme events from the stationary background by applying time-frequency analysis methods (e.g., wavelets [42]). Then, the stationary background can be approximated with a Gaussian stationary stochastic process (with properly tuned covariance function), while the rare event component can be represented with a Poisson process with properly calibrated parameters that represent the characteristics of the extreme forcing events (frequency and magnitude).

where $xb$ accounts for the background state (nonextreme) and $xr$ captures extreme responses (due to the intermittent forcing)—see Fig. 2. More precisely, $xr$ is the system response under two conditions: (1) the forcing is given by $F=Fr$ (i.e., we have an impulse) and (2) the norm of the response is greater than the background response fluctuations according to a given criterion, e.g., $\Vert x\Vert >\gamma $. However, as we will see in Sec. 3, other criteria may be used. These rare transitions occur when we have a forcing impulse and include a period where the system relaxes to the background state $xb$. The background component $xb$ corresponds to the response without rare events $xb=x\u2212xr$, and in this regime, the dynamics are primarily governed by the background forcing term $Fb$.

We require that rare events are statistically independent from each other. In the general formulation of the PDS method, we also need to assume that rare events have negligible effects on the background state $xb$, but in the current setup this assumption is not necessary due to the linear character of the examples we consider. However, in order to apply the method for general nonlinear structural systems, we need to satisfy this condition. Also, we require ergodicity since the statistics we approximate correspond to long time averages of trajectories [40].

*t*=

*t*

_{0}, we have an arbitrary background state $u(t0)=ub$ as the initial condition. The problem in this regime can be formulated as

Under the assumption of independent rare events, we can use Eq. (5) as the starting point for analytical and numerical approximations for the statistical response during rare events.

Because of the nonintermittent character of the response 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 $Fb(t)$ is a Gaussian process or can be adequately approximated by a Gaussian process, in which case the response statistics is readily obtained from the associated Fokker–Planck equation. Consequently, this step provides us with the statistical steady-state probability distribution for the mode of interest under the condition that the dynamics “live” in the stochastic background.

*q*may be any function of interest involving the response. In the last equation, $Pr$ denotes the overall rare event probability. This is defined as the probability of the response exceeding a threshold

*γ*due to extreme event excitations

where $1(\u2009\xb7\u2009)$ is the indicator function. The rare event probability measures the total rare event duration.

The utility of the proposed decomposition is its flexibility and effectiveness in cheaply capturing rare responses, since we account for rare event dynamics directly by connecting their statistical properties with the original system response.

### Problem Formulation for Linear Single-Degree-of-Freedom Systems.

where *k* is the stiffness, *λ* is the damping, and $\zeta =\lambda /2k$ is the damping ratio. For what follows, we adopt the standard definitions: $\omega n=k$, $\omega o=\omega n\zeta 2\u22121$, and $\omega d=\omega n1\u2212\zeta 2$.

Here, *F _{b}* is the background forcing component that has a characteristic magnitude

*σ*, and

_{b}*F*is a rare and large amplitude forcing component that has a characteristic magnitude $\sigma r$, which is much larger than the magnitude of the background forcing $\sigma r\u226b\sigma b$. Despite the simplicity of the system, its response may feature a significantly complicated statistical structure with heavy-tailed characteristics (large kurtosis) and asymmetry (skewness).

_{r}where *q* controls the magnitude of the forcing.

The second forcing term in Eq. (11) describes rare and extreme events. In particular, we assume that this component is a random impulse train ($\delta (\u2009\xb7\u2009)$ is a unit impulse), where *N*(*t*) is a Poisson counting process that represents the number of impulses that arrive in the time interval $0<t\u2264T$, *α* is the impulse mean magnitude (characterizing the rare event magnitude *σ _{r}*), which we assume is normally distributed with mean $\mu \alpha $, variance $\sigma \alpha 2$ and independent from the state of the system. In addition, the arrival rate is constant and given by $\nu \alpha $ (or by the mean arrival time $T\alpha =1/\nu \alpha $ so that impulse arrival times are exponentially distributed $\tau \u223ceT\alpha $).

where $\sigma h\u0307$ is the standard deviation of $h\u0307(t)$. The factor *m* sets the severity of the extreme impulse events, with respect to the background amplitude level, so that $\sigma r\u226b\sigma b$.

The setup above is applicable to numerous applications in other domains (with appropriate modifications on the forcing, e.g., the background excitation spectrum shape), including structures under wind excitations, systems under seismic excitations, and vibrations of high-speed crafts and road vehicles [5,13,14].

## Analytical Probability Density Functions of Single-Degree-of-Freedom Systems for Limiting Cases

In this section, we apply the probabilistic decomposition-synthesis method to derive analytical approximations for the pdf of the displacement, velocity, and acceleration for two special limiting cases: underdamped motion with $\zeta \u226a1$ and overdamped motion with $\zeta \u226b1$. We first perform the analysis for the response displacement and then derive the response pdf for the velocity and acceleration.

### Background Response Probability Density Functions.

where $R(\sigma )$ denotes a Rayleigh distribution with parameter *σ*.

### Analytical Probability Density Functions for the Underdamped Case ζ≪1.

In the strongly underdamped case with $\zeta \u226a1$, due to the highly oscillatory nature of the response, we focus on deriving the statistics of the local extrema. Consequently, we present results in terms of the *envelope* of the response statistics.

#### Rare Events Response.

*α*. The actual value of the response

*x*is considered negligible. For the present case, $\zeta \u226a1$, we have that the envelopes of the response (displacement, velocity, and acceleration) during the rare event are given by (see Appendix A for details)

_{b}*α*in the term $x\u0307b+\alpha $ are both Gaussian distributed and independent, and therefore, their sum is also Gaussian distributed

where $\sigma |\eta |=\sigma x\u0307b2+\sigma \alpha 2$.

#### Rare Event Probability.

Next, we compute the rare event probability: the total rare event period over a fixed time interval, as defined in Eq. (8). This is done by employing an appropriate definition of an extreme event in terms of a threshold value. One possible option is to set an absolute threshold *γ*. However, in the current context, it is far more convenient to set this threshold relative to the value of the local maximum of the extreme event response. Specifically, the time duration *τ _{e}* a rare response takes to return back to the background state will be given by the duration starting from the initial impulse event time ($t0$) to the point where the response has decayed back to $\rho c$ (or $100\rho c%$) of its absolute maximum; here and throughout this paper, we set $\rho c=0.1$. This is a value that we considered without tuning. We emphasize that the derived approximation is not sensitive to the exact value of

*ρ*as long as this has been chosen within reasonable values.

_{c}*the typical duration τ*. Using this rare event duration

_{e}is independent of the background state or the impulse impact intensity*τ*, we can compute the total rare event probability from the arrival frequency $\nu \alpha $ (recall $\nu \alpha =1/T\alpha $) of the Poisson process that describes the frequency of impulse events

_{e}Note, since we naturally assume that extreme events are rare enough to be statistical independent, the above probability is smaller than 0.5.

#### Conditional Probability Density Functions for Rare Events.

*τ*(from Eq. (22)) when the response has relaxed back to the background dynamics

_{e}where we have already derived the pdf for $f|\eta |$ in Eq. (20). What remains is the derivation of the conditional pdf for $fur\u2223|\eta |$.

where $s(\u2009\xb7\u2009)$ denotes the step function, which is equal to 1 when the argument is greater than or equal to 0 and 0 otherwise. Refer to Appendix B for a detailed derivation of the above.

#### Summary of Results for the Underdamped Case

##### Displacement envelope.

where $\tau e=\u2212log(\rho c)/(\zeta \omega n)\u2009$ and $s(\u2009\xb7\u2009)$ denotes the step function.

##### Velocity envelope.

##### Acceleration envelope.

#### Comparison With Monte Carlo Simulations.

For Monte Carlo simulations, the excitation time series is generated by superimposing the background and rare event components together. The background excitation, described by a stationary stochastic process with a Pierson–Moskowitz spectrum (Eq. (12)), is simulated through a superposition of cosines over a range of frequencies with corresponding amplitudes and uniformly distributed random phases. The intermittent component is the random impulse train, and each impact is introduced as a velocity jump with a given magnitude at the point of the impulse impact. For each of the comparisons, we generate ten realizations of the excitation time series, each with a train of 100 impulses. Once each time series for the excitation is computed, the governing ordinary differential equations are solved using a 4th/5th-order Runge–Kutta method (we carefully account for the modifications in the momentum that an impulse imparts by integrating up to each impulse time and modifying the initial conditions that the impulse imparts before integrating the system to the next impulse time). We ensured that the time series length and ensemble size lead to converged response statistics (displacement, velocity, and acceleration).

Furthermore, we utilize a shifted Pierson–Moskowitz spectrum $Shh(\omega \u22121)$ in order to avoid resonance conditions. The other parameters and resulting system statistical properties are given in Table 1. As seen in Fig. 3, the analytical approximations compare well with the Monte Carlo simulations many standard deviations away from the zero mean. The results are robust to different parameters as long as the assumption of independent (nonoverlapping) extreme events is satisfied.

Some of the discrepancies observed between the Monte Carlo simulations and analytical results can be attributed to the envelope approximation used to derive the conditionally rare event pdf. Indeed, these discrepancies are significantly reduced (results not shown) if we utilize the semi-analytical method presented in Sec. 4, where no simplifications are made for the form of the rare event response.

### Analytical Probability Density Functions for the Overdamped Case ζ≫1.

In Sec. 3.2, we derived the analytical response pdf under the underdamped assumption $\zeta \u226a1$. Here, we briefly summarize the results for the response pdf in the overdamped case with $\zeta \u226b1$. The derivation follows the exact same steps as in Secs. 3.1 and 3.2, but instead uses the corresponding formulas for the rare event response in overdamped conditions $\zeta \u226b1$ (see Appendix A). *An important difference in the overdamped scenario is that the system does not exhibit oscillatory response, hence we directly operate on the response instead of the envelope of the response*.

#### Displacement.

where $\tau e,dis=(\pi /2\omega o)\u2212log\u2009\rho c/(\zeta \omega n+\omega o)$.

#### Velocity.

where $\tau e,vel=\u2212(1/\zeta \omega n+\omega o)log\u2009\rho c$.

#### Acceleration.

where $\tau e,acc=\u2212log\u2009\rho c/(\zeta \omega +\omega o)$.

Note that here, we do not have the simple scaling as in the underdamped case for the conditionally rare component of the pdf.

#### Comparison With Monte Carlo Simulations.

We confirm the accuracy of the analytical results given in Eqs. (38), (40), and (42) for the strongly overdamped case through comparison with direct Monte Carlo simulations. The parameters and resulting statistical quantities of the system are given in Table 3. The analytical expressions for the response show excellent agreement compared to the Monte Carlo simulations, see Fig. 4 and Table 4.

## Semi-Analytical Probability Density Functions of the Response of Single-Degree-of-Freedom Systems

We now formulate a semi-analytical approach to quantify the response pdf for arbitrary system parameters (including the severely underdamped or overdamped cases considered previously). The approach here adapts the numerical scheme described in Ref. [1] for systems undergoing internal instabilities.

Whereas for the limiting cases that were previously studied, where (time series) knowledge of the system trajectory ($xr(t)$ or $ur(t)$) could be translated to information about the corresponding pdf ($fxr$ or $fur$), this is not always possible. In addition, *nonlinear* structural systems, in general, have analytical expressions for the rare event transitions. For these cases, we can compute the rare event statistics by numerically approximating the corresponding histogram, using either analytical or numerically generated trajectories in the rare event regime.

### Numerical Computation of Rare Events Statistics.

In Sec. 3, the derivation consisted of estimating all three unknown quantities: the background density $fxb$, the rare event density $fxr$, and the rare event probability $Pr$ analytically. However, in the semi-analytical scheme, we will obtain the rare event density $fxr$ and the rare event probability $Pr$ by taking a histogram of the numerically simulated analytical form of the rare response. The background density $fxb$ will still be obtained analytically as in Sec. 3.1.

*t*= 0 (the beginning of the rare event) until the end of the rare event at $t=\tau e$. The conditional density of rare event response for the velocity and acceleration are, similarly, given by

### Numerical Estimation of the Rare Events Probability.

*τ*, the absolute maximum of the response is obtained from the numerical simulations of the trajectories. Once the rare event duration has been determined, we compute the total rare event probability, just as before, from

_{e}### Semi-Analytical Probability Density Functions.

The corresponding pdf for the velocity $fx\u0307$ and acceleration $fx\u0308$ can be computed with the same formula but with the appropriate variance for the Gaussian core ($\sigma x\u0307b$ or $\sigma x\u0308b$), rare event duration ($\tau e,vel$ or $\tau e,acc$), and histograms ($x\u0307r\u2223\eta (t\u2223n)$ or $x\u0308r\u2223\eta (t\u2223n)$).

#### Comparison With Monte Carlo Simulations.

For illustration, a SDOF configuration is considered with critical damping ratio, *ζ* = 1. The detailed parameters and relevant statistical quantities of the system are given in Table 5. This is a regime where the analytical results derived in Sec. 3 are not applicable. Even for this *ζ* value, the semi-analytical pdf for the response shows excellent agreement with direct simulations (Fig. 5 and Table 6). We emphasize that the computational cost of the semi-analytical scheme is comparable with that of the analytical approximations (order of seconds) and both are significantly lower than the cost of Monte Carlo simulation (order of hours).

## Semi-Analytical Probability Density Functions for the Response of Multidegree-of-Freedom Systems

An important advantage of the semi-analytical scheme is the straightforward applicability of the algorithm to multidegree-of-freedom systems. In this section, we demonstrate how the extension can be made for a two-degrees-of-freedom (TDOF) linear system.

where the stochastic forcing $F(t)=Fb(t)+Fr(t)$ is applied to the first mass (mass *m*) and *x*, *y* are displacements of the two masses. As before, $Fb(t)$ is the background component and $Fr(t)=\u2211i=1N(t)\alpha i\u2009\delta (t\u2212\tau i)$ is the rare event component.

The background statistics are obtained by analyzing the response spectrum of the TDOF system subject to the background excitation component. Details of this derivation are provided in Appendix C. The histograms for the rare event transitions can be computed through standard analytical expressions that can be derived for linear systems, under the set of initial conditions: $x0=y0=y\u03070=0$ and $x\u03070=n.$

where *z* can be either of the degrees-of-freedom (*x* or *y*) or the corresponding velocities or accelerations, while $\tau ez$ is the typical duration of the rare events, which is estimated numerically. For each desired quantity *z*, the corresponding variance under background forcing, temporal duration of the rare events, and the rare event histogram need to be utilized.

Results are presented for the pdf of the displacement, velocity, and acceleration of each degree-of-freedom, see Fig. 7. These compare favorably with the direct Monte Carlo simulations. The parameters and resulting statistical quantities of the system are provided in Table 7 with error measures in Table 8. Further numerical simulations (not presented) demonstrated strong robustness of the approach to different system parameters.

## Semi-Analytical Probability Density Functions of Local Maxima

It is straightforward to extend the semi-analytical framework for other quantities of interest, such as the local extrema (maxima and minima) of the response, which is of interest in reliability applications. In this case, the numerically simulated analytical trajectories of rare responses can be used to compute the numerical histogram for the local extrema instead of the full response. For the conditionally nonextreme component, we can use known results from the theory of stationary Gaussian stochastic processes to describe the corresponding background pdf analytically.

### Distribution of Local Maxima Under Background Excitation.

*ϵ*, the probability density function of positive extrema (maxima) is given by [44,45]

*x*is the magnitude of the maxima, spectral bandwidth $\u03f5=(1\u2212\mu 22/(\mu 0\mu 4))$, and $\Phi (x)=(1/2\pi )\u222b\u2212\u221exe\u2212u2/2du$ is the standard normal cumulative distribution function. The spectral moments for the background response displacement

*x*are also defined as

_{b}We note that that in the limit of an infinitesimally narrowbanded signal $(\u03f5=0)$, the pdf converges to a Rayleigh distribution. On the other hand, for an infinitely broadbanded signal $(\u03f5=1)$, the distribution converges to the Gaussian pdf. For a signal with in between spectral bandwidth $(0\u2264\u03f5\u22641)$, the pdf has a blended structure with the form in Eq. (55).

where the $x\u0302$ notation denotes the local extrema of *x*. Similar expressions can be obtained for the velocity and acceleration extrema.

### Statistics of Local Extrema During Rare Transitions.

where $M(\u2009\xb7\u2009)$ is the operator that provides all the positive/negative extrema. The positive/negative extrema are defined as the points where the derivative of the signal is zero.

### Semi-Analytical Probability Density Function for Local Extrema.

The same decomposition can be used for the velocity and acceleration local maxima. We compare the semi-analytical decomposition with Monte Carlo simulations. In Fig. 8, we present results for the two-degrees-of-freedom system. The pdf is shown for local extrema of the displacement, velocity, and acceleration for each degree-of-freedom. We emphasize the nontrivial structure of the pdf and especially the tails. We observe from these comparisons that the semi-analytical scheme accurately estimates both the tails of the response and the non-Gaussian structure of core of the distribution.

## Summary and Conclusions

We have formulated a robust approximation method to quantify the probabilistic response of structural systems subjected to stochastic excitation containing intermittent components. The foundation of our approach is the recently developed probabilistic decomposition-synthesis method for the quantification of rare events due to internal instabilities to the problem where extreme responses are triggered by external forcing. The intermittent forcing is represented as a background component, modeled by a colored stochastic processes with energy distributed across a range of frequencies, and additionally a rare and extreme component that can be represented by impulses that are Poisson distributed with large interarrival time. Owing to the nature of the forcing, even the probabilistic response of a linear system can be highly complex with asymmetry and complex non-Gaussian tail behavior, which would be expected if the forcing did not consist of intermittently extreme events.

The main result of this work is the derivation of analytical and semi-analytical expressions for the pdf of the response and its local extrema for structural systems (including the displacement, velocity, and acceleration pdf). These expressions decompose the pdf into a statistical core, capturing the response under background forcing, as well as a heavy-tailed component associated with extreme transitions due to rare impulse impacts. We performed a thorough analysis for linear SDOF systems under various system parameters and also derived analytical formulas for two special parameter ranges (lightly damped or heavily damped systems). The general semi-analytical decomposition is applicable for arbitrary system parameters and we have demonstrated its validity through comprehensive comparisons with Monte Carlo simulations. The general framework is also directly applicable to multidegree-of-freedom systems, as well as systems with nonlinearities. We demonstrated applicability to a 2DOF linear system consisting of two coupled masses, where the first mass is excited. Modifications of the method to compute statistics of local extrema have also been presented.

The developed approach allows for computation of the response pdf of structural systems many orders of magnitude faster than a direct Monte Carlo simulation, which is currently the only reliable tool for such computations. The rapid evaluation of response pdfs for systems excited by extreme forcing by the proposed method paves the way for enabling robust design of structural systems subjected to extreme events of a stochastic nature. In such cases, it is usually not feasible to run Monte Carlo simulations for various parameter sets during the design process owing to the computational costs associated with low probability rare events, let alone perform parametric optimization. In different work [41], we demonstrate application of the method to design and optimization in ocean engineering applications and show that the approach is well suited to predict optimal extreme-event mitigating designs and for system reliability.

## Acknowledgment

H. K. J. and M. A. M. have been supported through the first and third grants as graduate students. We are also grateful to the Samsung Scholarship Program for support of H. K. J. as well as the MIT Energy Initiative for support under the grant Nonlinear Energy Harvesting From Broad-Band Vibrational Sources By Mimicking Turbulent Energy Transfer Mechanisms.

## Funding Data

Air Force Office of Scientific Research (Grant No. FA9550-16-1-0231).

MITEI (Nonlinear Energy Harvesting From Broad-Band Vibrational Sources By Mimicking Turbulent Energy Transfer Mechanisms).

Office of Naval Research (Grant Nos. N00014-14-1-0520 and N00014-15-1-2381).

### Appendix A: Impulse Response of Single-Degree-of-Freedom Systems

*α*at an arbitrary time

*t*

_{0}, say $t0=0$, and with zero initial value but nonzero initial velocity ($(xr,x\u0307r)=(0,x\u0307r0)$ at $t=0\u2212$) are given by the following equations under the two limiting cases of damping: