Stochastic reduced-order models for stable nonlinear ordinary differential equations

Two methods based on stochastic reducedorder models (SROM) are proposed to solve stochastic stable nonlinear ordinary differential equations. One general method available for the probabilistic characterization of the response of nonlinear systems subjected to random excitation is Monte Carlo (MC), wherein the response of the nonlinear system must be calculated for a large number of samples of the input, which can be very computationally demanding. Random vibration theory is also inadequate for calculating response statistics for both linear systems under non-Gaussian inputs and nonlinear systems subjected to any kind of excitation. The two methods proposed are based on SROM, i.e., stochastic models with a finite number of optimally selected samples. The first method uses a SROMmodel for the random input. The second method is based on a surrogate model for the response of the nonlinear system defined on a Voronoi tessellation of the input samples. The newly proposed methods are applied for stable nonlinear ordinary differential equations, with deterministic coefficients and stochastic input, that are used in engineering applications: single-degree-of-freedom Duffing and Bouc– Wen systems, and a two-degree-of-freedom nonlinear energy sink system. The numerical results suggest that SROMs are able to estimate statistics of the stochastic responses for these systems efficiently and accurately, results validated by the benchmark MC results. A. Radu (B) University of Bristol, Bristol BS8 1TR, UK e-mail: alin.radu@bristol.ac.uk


Introduction
Random vibration theory and numerical methods to date provide efficient solutions for second-moment properties of the states of arbitrary linear systems subjected to Gaussian random input. Unless the input is Gaussian and the system is linear or the output can be assumed to be Gaussian, the first two moments of the response are insufficient for calculating response statistics [1] (Sect. 5.3). This is a significant limitation since realistic dynamic systems behave nonlinearly under strong random vibrations and the distribution of the response is non-Gaussian. For example, the distribution of the pressure field acting on a spacecraft during atmospheric re-entry is highly non-Gaussian [2]; the assumption that the responses of structures under extreme wind [3] or seismic [4] excitations are linear may be inadequate; the aircraft main landing-gears are nonlinear complex systems, which require time-consuming dynamic simulations to calculate steady-state solutions [5]. Thus, the reliability of such structures cannot be obtained within the framework of the random vibration theory. The conceptually simple framework of the random vibration theory cannot be used in such complex applications since some concepts cannot even be defined. The first-and second-moment equations for linear systems in random vibration theory are formal since it views white noise as a process with constant spectral density, i.e. a process with infinite variance, which is not defined in the second-moment sense [6] (Sect. 7.2).
Monte Carlo (MC) is one general and reliable method available for solving complex nonlinear dynamic systems subjected to complex stochastic input, i.e. non-Gaussian, non-stationary non-white noise. The method is asymptotically correct since it approaches the exact solution as the number of simulations increases indefinitely. However, MC simulation is computationally impractical for complex realistic nonlinear systems [5,7,8] since it involves repeated deterministic dynamic analyses for randomly selected samples of the random excitation. Efficient alternatives to the MC method have been studied before, such as the quasi-Monte Carlo simulation [7,[9][10][11], Latin hypercube [12,13], or change of measure [6] (Sect. 5.4). The limitation of the MC method has encouraged the development of approximate methods for finding statistical properties of the dynamic systems under random excitation. Some approaches are based on heuristic arguments, such as the equivalent stochastic linearisation [14][15][16], while others rely on rigorous but rather restrictive conditions, as in the stochastic averaging method [17][18][19] or perturbation [1] (Sect. 6.2), [20]. Other well-tested methodologies used for the calculation of response statistics of nonlinear systems involve the use of the Fokker Plank equation [21][22][23], moments and cumulant closures [24,25], stochastic modelling using polynomial chaos [26,27], or reduced orders of the nonlinear system by using linear or nonlinear normal modes [28][29][30]. This enumeration of methods used for solving stochastic nonlinear equations is not exhaustive, but so far using these approaches proved to make capturing statistics of extreme responses of nonlinear systems difficult or even impossible due to their built-in limitations [31], because they have solutions only for special cases [32], or are too computationally expensive [33].
The goal of this paper is to propose novel methods that can be used to calculate response statistics of engineering systems, characterised by stable, nonlinear ordinary differential equations with deterministic coefficients and general non-Gaussian, non-stationary stochastic input. Even though the methodology proposed is general for the described setup, the study of complex systems such as turbulent or chaotic systems, or infinite-dimensional systems with multiple instabilities is beyond the scope of this paper. However, considerable advances on developing efficient methods to quantify the uncertainty in complex stochastic highly nonlinear systems like Earth's climate or oceans [28,34] must be acknowledged. The two methods proposed in the current study to solve the aforementioned engineering systems are conceptually simple, accurate, non-intrusive and computationally efficient and are based on stochastic reduced-order models (SROMs) [35]. Stochastic reduced-order models are stochastic processes that have a finite number of samples selected in an optimal manner from the samples of the target process. Like MC simulation, the methods use samples of seismic load processes to characterize structural response and are not intrusive in the sense that their construction uses deterministic solutions. However, the proposed SROMs use small numbers of input samples selected in an optimal manner. In contrast, MC simulation uses a large number of samples selected at random. The use of optimally selected samples allows to reduce the number of simulations required by the MC method by one or two orders of magnitude while retaining accuracy. The method has been originally developed for dynamic response [36]. It has been shown that for static problems the SROM-based method can be improved significantly [37], and an adaptive method to further refine SROMs used in solving stochastic equations was proposed in [38]. Recently, SROMs have already been used successfully in applications for the quantification of the uncertainty in electromagnetic-signals interference in cables [39], or in the inter-granular corrosion rates [40]. SROMs have also been used in solving inverse problems, with applications in the identification of material properties in elastodynamics [41], traditionally solved using Bayesian inferences or stochastic-optimisation approaches.
The first method presented in this paper, referred to as the SROM method, constructs a SROM using directly the input samples. The second method, referred to as the extended SROM method , i.e., the ESROM method, builds upon the SROM method to constructs an SROM-based surrogate model for the response of the dynamic system, defined on a Voronoi tessellation of the input samples. Both methods proposed in this paper may be used as alternatives to MC for calculating response statistics for dynamic nonlinear equations, with just a fraction of the computational effort required by MC. The SROM method provides quick response statistics within some error of the MC results, while the superior ESROM method can provide reliable results even in the extreme tails of the response probability distributions. Thus, the ESROM would allow for a simulation-based design of highly reliable nonlinear engineering systems for probability of no-collapse of the 10 −6 order for earthquake-resistant structures [42], or probability per flight hour of catastrophic failure for aircrafts as low as 10 −9 [43], probabilities that otherwise would require billions or even trillions of MC dynamic analyses [7]. Applications of the methods proposed are shown for simple second-order, oneand two-degree-of-freedom nonlinear dynamic ordinary differential equations (ODEs) with deterministic coefficients and stochastic input. Response statistics such as probability tail distributions and stochastic moments are compared with the reference MC results.

Nonlinear dynamic system definition
Nonlinear dynamic systems subjected to random input can be represented by nonlinear stochastic ODEs with the general form N j=0 Ψ j t, X (t), Y ( j) (t) = 0, t ≥ 0, where Ψ j are linear or nonlinear differentiable functions of time t with deterministic or stochastic coefficients; Y ( j) (t) = d j Y (t)/dt j , j = 1, . . . , N are the jth derivatives with respect to time of the state Y (t) of the dynamic system; and X (t) is the stochastic input. Functions Ψ j for the example systems in Eqs. (2)(3)(4)(5) used to demonstrate how the proposed framework works, have deterministic coefficients, and are power, or exponential functions, but any other nonlinear differential function would work too. The response of the nonlinear dynamic system is described by the state Y (t), for given initial conditions Y (k) (0), k = 1, . . . , n − 1.
The stochastic process X (t) can be any stochastic process, and thus we assume a type of process that fits the general case of a zero-mean, non-Gaussian, nonstationary process: where g(ν) = (νθ 2 √ 2π) −1 exp{−0.5(ln(ν)−θ 1 ) 2 /θ 2 2 } is a one-sided power spectral density function, with parameters (θ 1 , θ 2 ), that provides the second-order moment properties and the frequency content of X (t); f (t) = θ 3 t θ 4 exp(θ 5 t) is an amplitude-modulation function with parameters (θ 3 , θ 4 , θ 5 ) that describes the non-stationary character of X (t); Φ(.) is the cumulative probability distribution function for the standard normal distribution; G(t; g(ν)) is a standard Gaussian process with second-order moment properties given by the function g(ν); and F T (.; n T ) is the Student's T distribution with n T degrees of freedom. Functions f (t), g(ν) and F T (.; n T ) are customary and define the probability law of X (t), but can be replaced by any other functions that describe the amplitude non-stationarity character, the frequency content, and the distribution of X (t), respectively. Function f (t) is called a Gamma function and has been used before to describe, for example, the non-stationarity in seismic ground motions [44]. The choice of the Student's T distribution for X (t) is justified by its heavy tails specific to heavy-tailed phenomena to which engineering systems are subjected, such as: offshore structures subjected rough waves [45], high-rise structures subjected to strong wind [46], or earthquake [4] loads. The distribution F T (.; n T ) can be replaced by other more complex, multi-modal distributions. However, for the purpose of loads with high peaks, the Student's T distribution has a kurtosis higher than 3 that can be fit to data, and thus experience tails higher than Gaussian processes, similar to realistic loads with this characteristics, such as: the wind pressure on low-rise structures tested in wind tunnels [47]; the coastal-wave elevations measured in Duck, North Carolina [48]; the seismic ground-acceleration for earthquakes recorded on rock sites [44]; or the unevenness of a railway track in India [49]. The representation of non-Gaussian processes in Eq. (1) is known as a memoryless transformation [50] and provides an intuitive way to simulate samples of non-Gaussian processes. Numerical examples are shown for the scalar parameters (θ 1,...,5 ) = (2, 0.5, 0.5, 1.5, −0.35) and n T = 3. Note that the process X (t) has variance ν≥0 g(ν)dν = 1, which can also be customised to a value σ 2 X , by multiplying g(ν) by it. Figure 1a, b shows the function g(ν) for two sets of parameters (θ 1 , θ 2 ) = {(2, 0.5), (1, 1)}, and corresponding samples of X (t), respectively.
Three simple nonlinear ODEs with various applications in engineering are used as examples. The first system is the single-degree-of-freedom (SDOF) Duffing oscillator [51,52] used, for example, in modelling 123 Fig. 1 a Power spectral density function g(ν), b samples of the input process vibration-control systems with cubic stiffness [53] in mechanical engineering: (2) The second system described by the following ODEs is the Bouc-Wen SDOF system [54,55] used, for example, to model buildings' seismic ductility demand [56] in civil engineering: where the process W (t) is known as the hysteresis response of the Bouc-Wen oscillator. Finally, the last system solved as an example for the proposed methodologies is the two-degree-of-freedom (2DOF) nonlinear energy sink (NES) system [57], used for the reduction of the response of ocean-engineering systems subjected to extreme loading [58], whose state is described by the bi-variate stochastic process In all three systems of ODEs in Eqs. (2)(3)(4)(5), 2 are the first-and second-order derivatives of the systems' states Y (t) with respect to time t, and zero initial conditions are assumed. The following systems' parameters ν = π ,

Solutions and response statistics
Monte Carlo (MC) simulation is a general method available for finding statistical properties of the state Y (t) of an arbitrary nonlinear dynamic system subjected to stochastic input. However, the method is impractical for complex systems due to its high computational cost. The objective of this paper is to adopt practical, accurate and efficient methods for calculating statistics of the response Y (t) of nonlinear dynamic systems under general stochastic excitation. We propose two novel methods for calculating response statistics of the response Y (t) of nonlinear dynamic systems subjected to stochastic input X (t). Both methods are based on stochastic reduced-order models (SROMs) for the input: one, more straightforward, uses the SROM directly to find the response of the system for a small number of optimally selected samples; the second one, referred to as the extended SROM, uses the input SROM to construct a surrogate model for the system's response based on a first-order Taylor expansion. A SROM of X (t) is a stochastic processX (t) represented by a finite number m of samplesx k (t) of X (t) with probabilities p k , k = 1, . . . , m, such that the probability laws of X (t) andX (t) are similar in some sense.
The SROM-based methods will be tested against MC results, through the comparison of statistics of the state Y (t). In this case, we will aim to compare the tail distributions and high-order moments. MC method requires the calculation of the response y k (t) for a large number n of input samples x k (t), k = 1, . . . , n. The moments of order q, and the tail distribution function can be estimated using the MC response samples, respectively, as: where 1(.) is the indicator function. The two methods proposed are presented in the following two sections. The first one, involves the construction of a SROMX (t) for the input X (t), with m << n samples. The second method, referred to as the extended SROM (ESROM), constructs a surrogate modelỸ (t) for the state Y (t) of the system, based onX (t).

SROM-based solution
Let {(x l (t), l = 1, . . . , n} be a large enough set of samples of the input process X (t) that can characterise its probability law. Any number m << n of samples . Similar to MC, the SROM also uses random samples of X (t), but unlike in the case of MC they are not equally likely, but weighed by distinct probabilities. It is shown in [36,59] that it is possible to select a small number m of samples andX (t) have similar probability laws. We are looking for a SROMX (t) with a range {x k (t), k = 1, . . . , m} of a relatively small number m of independent samples of X (t). Our objective is to find p k such that the discrepancies between the probability laws of the SROMX (t) = {(x k (t), p k ), k = 1, . . . , m} and X (t) are minimized [60]. ProcessesX (t) and X (t) are defined on the same probability space, and in order for the SROM X (t) to characterize the original model X (t), the m samples selected need to be largely spaced from each other in order to explore the entire range of X (t) samples, rather than being clustered in small subsets. Algorithms such as the integer optimisation [61], dependent thinning [62] or pattern classification [63] may be used to select the samples that defineX (t). The probability laws of stochastic processes are defined by their moments, marginal distributions and covariance functions. Target moments μ X (t; q) of order q, marginal distribution F X (x; t) and covariance function Σ X (t, s) of X (t) can be calculated directly from its probability law, if available, or from its samples. The statistics for the SROMX (t) are estimated from a selection of m samples, weighed by their probabilities. Thus, the momentsμ X (t; q) of order q, the marginal distribu-tionF X (x; t) and the covariance functionΣ X (t, s) of the SROMX (t) are calculated, respectively, as: Note that p k are the likelihood probabilities assigned to the input samplesx k (t) of the SROM and they have the following properties: (1) if m is very large andx k are independent samples of X , then optimal probabilities p k ≈ 1/m, as it approaches the case of Monte Carlo when all samples are equally likely; (2) the moments μ X (t; q) and distributionF X (x; t) ofX (t) are unbiased estimators of the corresponding statistics of X (t) with variances depending on the range ofX (t); and is the solution of a stochastic algebraic, differential or integral equation depending on X (t), the solutionỸ (t) of this equation withX (t) in place of X (t) is given by the samplesỹ k (t) = h(x k (t)) with probabilities p k . The last property is essential for the success of the SROM method, and its proof relies on the argument that if the probabilities p k partition the probability space of the input X (t), then the expected solution h(X (t)) in a given partition char-123 acterised by (x k (t), p k ) is given by the deterministic solutionỹ k (t) = h(x k (t)) with probability p k . More comprehensive arguments and detailed proves using measure theory for these properties are presented in [35]. Optimization algorithms for the construction of SROMs have been developed in [35,41]. A least mean square algorithm is used to identify the optimal probabilities p k that minimise the discrepancies between X (t) andX (t) herein, for a given cardinal m and a set of samplesx k (t). Based on the aforementioned property (1), the initial values of p k = 1/m are assumed as the starting point in the optimisation algorithm. Figure 2 shows the first two-order moments μ X (t; q), q = 1, 2 of the SROMX (t) for (a) m = 20 and (b) m = 10 3 , in comparison with the respective MC moments μ X (t; q), q = 1, 2 calculated from n = 10 4 samples. In order to already see the advantaged of SROMs, Fig. 2c, d shows 30 trials of calculating the same moments μ X (t; q), q = 1, 2 using MC, but just with a limited number of independent equally likely samples n = 20 and n = 10 3 , respectively. By comparing panels (a) with (c), and (b) with (d), respectively, we can see that unlike SROM, MC is incapable of capturing the properties of X (t) with just a few samples, the variability between the estimated of the moments being significant. Similar observations about the efficiency of the SROM vs. MC with the same small number of samples can be done regarding to the marginal distributions and the covariance functions, but due to the space limitations and given the scope of the paper, such additional preliminary results are not presented. Figure 3 shows (a) the marginal distribution of X (t) calculated by MC with n = 10 4 samples, and the marginal distributions Finally, following the property (3) of the SROM samples (x k (t), p k ), the moments of order q and the tail distribution of the system's response Y (t) shown in Eqs. (6) and (7) can be computed by the following proxies using the samplesỹ k (t), k = 1, . . . , m calculated as the response of Eqs. (2)(3)(4)(5) to the samplesx k (t) of the SROMX (t): Note that only m deterministic dynamic analyses are needed to obtain the solution of Y (t) in this case, i.e. the solution for Eqs. (2)(3)(4)(5) to the m samples of the SROMX (t). The essential benefit of the method consists in the fact that the SROMX (t) is built for the input X (t) and that the only dynamic analyses needed are the one for the samplesx k with probabilities p k since they characterise the probability law of X (t). It must also be emphasized that the pairs (x k (t), p k ) describe fully the probability law of the SROMX (t).

ESROM-based solution
The ESROM method develops a surrogate modelỸ (t) for the state Y (t) of the nonlinear ODEs employing the following four steps. First, a parametric mode for the input process X (t) is constructed: where Z = {Z i , i = 1, . . . , d} is a random vector and {ϕ i (t), , i = 1, . . . , d} are deterministic specified functions of time. Parametric models of this type can be calculated for any stochastic process. We use a truncated Karhunen-Loève model to calculate functions {ϕ i (t), i = 1, . . . , d} by calibrating it to the first two moments of X (t), that is, the mean and covariance functions [64]. Functions ϕ i (t) are the eigenvectors of the covariance matrix Σ(t, s) defined for 0 ≤ t, s ≤ τ . Unlike the Karhunen-Loève expansion, for which , t ≥ 0, where λ i are the eigenvectors and Z i are random variables with zero mean and unit variance, the random vector Z is characterized by its samples. Samples z k,i of Z i are obtained by minimizing the difference between the actual samples x k (t) of X (t) and their approximations is converging to the original process X (t) as the dimension d of the vector Z increases. Figure 5 shows one sample of X (t) and its approximation for two values (a) d = 50, and (b) d = 100.
In the second step, we construct a SROMZ = {(z k , p k ), k = 1, . . . , m} for Z using the samples of the random vector Z , following the procedure presented in the SROM-based solution, in Sect. 3.1. Note that n samples of z k of Z are available, one for each sample of x k (t) of X (t). As discussed previously, we seek to find optimum probabilities p k , k = 1, . . . , m for each samplez k defining the SROMZ of Z , with m k=1 p k = 1, such that the probability laws of the d-dimensional random vectors Z andZ are similar.
In the third step, we construct a Voronoi tessellation {Γ k , k = 1, . . . , m} on the range of Z , with centres {z k } [37,38]. Each cell Γ k contains all samples z j of Z that are closest to its centre, i.e. Γ k = {z j : |z j −z k | ≤ |z j −z l |, l = k}. Note that the probability p k defines also the probability that Z takes values in cell Γ k .
Finally, in the fourth step, a surrogate modelỸ (t; Z ) for Y (t) is constructed using local piece-wise linear approximations of Y (t; Z ) for each Voronoi cell by using a first-order Taylor expansion at the centre Z = z k of the cell: where Y (t;z k ) is the response Y (t) corresponding to the pair (z k , p k ) and ∇Y (t;z k ) = {∂Y (t; Z )/∂ Z i |Z = z k , i = 1, . . . , d} are the gradients of these response samples. These gradients can be calculated by solving the deterministic differential equations N j=0 ∂Ψ j X (t; Z ), Y ( j) (t) /∂ Z i = 0, i = 1, . . . , d, where N j=0 Ψ j t, X (t), Y ( j) (t) = 0, t ≥ 0 may be any nonlinear differentiable ODE as described previously, such as the Eqs. (2)(3)(4)(5). The gradients ∇Y (t;z k ) for the example ODEs are represented by ∂Y (t; Z )/∂ Z i , which are the solutions of Eqs. (2-5) differentiated with respect to each coordinate Z i of Z . Explicitly, the gradients ∇Y (t;z k ) for the Duffing ODE in Eq. (2) are the solutions of the following equation: where Finally, for the NES system defined in Eqs. (4)(5), the gradient ∇Y (t;z k ) is a bi-dimensional vector with components ∂Y 1 (t, Z )/∂ Z i and ∂Y 2 (t, Z ) /∂ Z i calculated by solving the following system of deterministic ODEs:  In case of more complex systems, for which explicit ODEs for the gradients cannot be written, they can be approximated numerically from the differences between the responses Y (t; Z ) calculated at Z =z k and Z =z k + Δz k , where Δz k denotes a perturbation ofz k . Also, higher-order approximationsỸ (t; Z ) can be constructed and used in the development of the surrogate model in Eq. (14), but they require the calculation of higher-order derivatives ofỸ (t; Z ) with respect to the coordinates of Z . Given that the nonlinear ODEs considered in this study are stable and are differentiable with respect to the systems' states Y (t), the gradients exist, and are solutions of nonlinear ODEs similar to the original ones. Thus, the effort to calculate ∂Y (t; Z )/∂ Z i is similar to calculating Y (t). Therefore, the number of dynamic analyses needed to construct the surrogate modelỸ (t; Z ) in Eq. (14) is comprised of m analyses to calculate Y (t;z k ) at each Voronoi cell centrez k , from Eqs. (2)(3)(4)(5) and m × d analyses to calculate ∇Y (t;z k )(Z −z k ) from Eqs. (15)(16)(17)(18)(19) for each Voronoi cell with respect to each coordinate of the random vector Z . The major advantage of ESROM is that any number of surrogate response samplesỹ k (t) can be simulated using Eq. (14) in order to increase the accuracy of the response statistics, without any additional computational effort, accounted in the fixed number m × (d + 1) of dynamic analyses needed to construct the surrogate model. Figure 6 shows the gradients ∂Y (t; Z )/∂ Z i calculated for the first three coordinates Z i , i = 1, 2, 3 of Z for the Voronoi cell centred around the first samplez 1 of the SROMZ , for each of the three nonlinear systems considered. Figure 7 shows the response samples y(t) (black, solid line) calculated by direct integration of Eqs.  Example. The following two-dimensional example is considered to illustrate graphically how the ESROM model works. Let's assume that the loading process X (t) is replaced by the bi-variate Gaussian vector X = (X 1 , X 2 ), where X 1 , X 2 ∼ N(0, 1) and correlation ρ X 1 ,X 2 = 0.3; and the state of the nonlinear system is described by the equation U (X ) = sin(X 1 ) cos(X 2 ). Figure 8a illustrates the first three steps: (i) the green dots are the samples of x k of X , (ii) the red stars are the m = 15 samplesx k of the SROMX , which are the centres of the (iii) Voronoi cells shown as the black polygons. Figure 8b, c shows the piece-wise linear surrogate modelŨ (X ) and the exact solution U (X ) (green mesh), respectively. Note that in this simple case, for which the input is a bi-variate random variable rather than a stochastic process, X and Z are identical.
The ESROM estimates of the response statistics defined in Eq. (6) and (7) are whereỹ k (t) are the samples ofỸ (t; Z ) calculated in Eq. (14). It is important to note that the response statistics in Eq. (20) and (21)  tics is not resumed to just m, but can go up to n, the total number of samples of X (t) available, since each sampleỹ k (t) is calculated directly from Eq. (14) using the sample of Z corresponding to each sample of X (t), as per the approximation in Eq. (13), without performing any additional dynamic analyses.

Numerical results and discussion
The performance of the two methods proposed in the previous sections is discussed in terms of the accuracy with which they are able to estimate statistics such as the moments or the tail probability distributions for the three nonlinear dynamic systems consid-ered, with respect to the MC estimates of the same results. Figures 9 and 10 show estimates of the firstand second-order moments μ(t, The accuracy and the efficiency of the SROM-based methods proposed depend on the SROM and ESROM models' sizes. For example, the accuracy of the SROM X (t) depends on the model size m, but also for a selected m it depends also on the optimization algorithm used to select its samplesx(t) and their probabilities p k . It was shown in [41] that the SROM solution converges to the MC solution as m → ∞. Thus, a large dimension m of the SROM would produce very accurate statistics of the process X (t), but it comes at a large computational cost, similar to Monte Carlo [35]. Even though accurate response statistics have been achieved with low values of m [40,41], the size m must be chosen on consideration of computational time. The accuracy of the response statistics for Y (t) depends on that of X (t) and the approximation used for the mapping of X (t) into Y (t). Bounds for the errors of SROM-based methods have been established in [59] and [38], but this analysis is beyond the scope of this paper. Figures 12 and 13 show the same first-and secondorder-moment estimates, similar to Figs. 9 and 10, but with an increased size of the SROMX (t) in the SROM method to m = 10 3 , which is similar to the computational effort required to build the response sur- While in the SROM method the estimates of the response statistics are based on a rather crude approximation of the mapping between the response samplesỹ k (t) and the SROM samplesx k (t), the ESROM method is based on a more accurate representation of this mapping, as shown in the previous section. The accuracy of the ESROM method depends as well on the dimension m of the SROM modelZ of Z , but also on the dimension d of the vector Z , which is used to parametrize the input process X (t). The effi-ciency of the methods proposed will be judged by the number of deterministic nonlinear dynamic analyses involved to calculate the response statistics. While MC requires n dynamic analyses, the SROM requires only m << n dynamic analyses for the samples (x k (t), p k ), k = 1, . . . , m of the SROM X (t). The ESROM method requires a larger number of dynamic analyses than SROM, that is, m × (d + 1) analyses. As a result, the efficiency of the ESROM method depends on the balance between the dimension d of Z and the dimension m of its SROMZ . In order to get a sense of how the balance between m and d plays a role in the development of the ESROM method, Fig. 14a shows the tail distributions for the response of the Bouc- Thus, the accurate representation of X (t) in Eq. (13), which is controlled by the dimension of Z (see Fig. 5), is essential for the success of the method. This conclusion is in line with expectation since the higher the dimension d is, the better the frequency content of X (t) is modelled, and the response of nonlinear dynamic systems is sensitive to the frequency content of the excitation. Figure 14b shows that an increase in the dimension m in the SROM method, provides better results also for the tail distribution of the response of the Bouc-Wen system, as already seen for the moments in Figs. 12 and 13. Further plots and discussions on the accuracy and efficiency of the two methods proposed for the tail distributions of the response of the other two systems are shown in Figs. 15 and 16. Figure 15 shows the response tail distributions in Fig. 11 in a logarithmic scale including for extreme values of the response, i.e. y > 1. It can be seen that the solutions provided by MC (n =   Figure 16 shows the same results as in Fig. 15, but for a larger number of MC samples, increased from n = 10 4 to n = 10 6 , and for a larger dimension m in the SROM method from m = 20 to m = 10 3 . A very large number of samples (n = 10 6 ) are required for the MC method to reach convergence in the extreme tails for values of y > 1. Regarding the SROM method, we notice that a large size m = 10 3 is still insufficient to improve results in the extreme tails. However, the ESROM results for the initial values m = 10, d = 100 provide excellent results, similar to the MC results calculated for n = 10 6 . Thus, with only approximately one thousandths of the computational effort, the ESROM is able to provide accurate results even in the tail of the response of the distribution of Y (t). Nevertheless, the SROM is a very efficient and reliable alternative to MC for calculating the other response statistics, at just a very small fraction of the computational effort.
Indisputable advantages of the ESROM method are shown in the analyses of extreme values of the state Y (t) of the nonlinear systems. In order to analyse the extreme values of Y (t), Fig. 15 shows the tail distributions for the Duffing and the Bouc-Wen models in a logarithmic scale. In the case of extreme-value analysis, by comparing the Figs. 15 and 16, we can conclude that the ESROM clearly outperforms MC method, managing to produce reliable results for just a small frac- tion of the deterministic analyses required by the MC method to reach convergence. The two SROM-based methods proposed are good alternatives to MC and provide accurate results efficiently, which allows for simulation-based results to be estimated accurately for nonlinear systems subjected to stochastic input. Finally, given the superior performance of the ESROM solution compared with MC throughout the distribution of the response, even in the case of extreme values, any statistics can be calculated with limited effort.
To support this idea, a more practical example of the joint distributions of the maxima of the input X (t) and the response Y (t) for the NES system is shown in Fig. 17. Since both processes X (t) and Y (t) are non-stationary, their joint distribution is impractical to show and thus the conditional probabilities of the maxima P[max t≥0 |Y j (t)| > y| max t≥0 |X (t)| > x] for each of the 2DOFs j = 1, 2 of the NES system are calculated and shown in Fig. 17a-d, calculated using the MC and the ESROM responses, respectively. The SROM results for this statistic are not shown since this method fails to perform accurately in case of large values of the response Y (t), and it does so for the joint distribution of the maxima of X (t) since it would be restricted to evaluate the maxima only through the small m number of input samples used, insufficient for such a result.

Conclusions
Two general, efficient and reliable methods are proposed for solving stochastic stable nonlinear dynamic equations. Stochastic reduced-order models (SROM), that is, random processes with finite number of samples, have been used to represent the stochastic inputs to nonlinear dynamic systems and construct approximations for the states of these systems. Two types of SROM-based solutions were presented, and numerical results were compared with the Monte Carlo (MC) estimates. The first method, that is, the SROM-based solution involves the construction of a SROM for the input process, whose samples are used to calculate the responses of the nonlinear systems. The extended version of this method, that is, the ESROM-based solu-tion, involves the construction of a piece-wise surrogate model for the states of the nonlinear systems, based on the SROM of the input. The accuracy and the efficiency of the two methods are discussed in terms of the number of dynamic analyses required for obtaining the reliable response statistics. Results show that the SROM method is very efficient producing reliable results within an acceptable error with respect to MC, with a small number of dynamic analyses, but it fails to represent statistics at the extremes. The ESROM is superior to the SROM method, and it outperforms in terms of computational time the MC method for some statistics, such as extreme-value distributions, or joint distributions of the input-response maxima.