A momentequationcopulaclosure method for nonlinear vibrational systems subjected to correlated noise
Abstract
We develop a momentequationcopulaclosure method for the inexpensive approximation of the steady state statistical structure of strongly nonlinear systems which are subjected to correlated excitations. Our approach relies on the derivation of moment equations that describe the dynamics governing the twotime statistics. These are combined with a nonGaussian pdf representation for the joint responseexcitation statistics, based on copula functions that has i) single time statistical structure consistent with the analytical solutions of the FokkerPlanck equation, and ii) twotime statistical structure with Gaussian characteristics. Through the adopted pdf representation, we derive a closure scheme which we formulate in terms of a consistency condition involving the second order statistics of the response, the closure constraint. A similar condition, the dynamics constraint, is also derived directly through the moment equations. These two constraints are formulated as a lowdimensional minimization problem with respect to the unknown parameters of the representation, the minimization of which imposes an interplay between the dynamics and the adopted closure. The new method allows for the semianalytical representation of the twotime, nonGaussian structure of the solution as well as the joint statistical structure of the responseexcitation over different time instants. We demonstrate its effectiveness through the application on bistable nonlinear singledegreeoffreedom energy harvesters with mechanical and electromagnetic damping, and we show that the results compare favorably with direct MonteCarlo simulations.
1 Introduction
In numerous systems in engineering, uncertainty in the dynamics is as important as the known conservation laws. Such an uncertainty can be introduced by external stochastic excitations, e.g. energy harvesters or structural systems subjected to ocean waves, wind excitations, earthquakes, and impact loads [1, 2, 3, 4, 5, 6]. For these cases, deterministic models cannot capture or even describe the essential features of the response and to this end, understanding of the system dynamics and optimization of its parameters for the desired performance is a challenging task. On the other hand, a probabilistic perspective can, in principle, provide such information but then the challenge is the numerical treatment of the resulted descriptive equations, which are normally associated with prohibitive computational cost.
The focal point of this work is the development of a semianalytical method for the inexpensive probabilistic description of nonlinear vibrational systems of low to moderate dimensionality subjected to correlated inputs. Depending on the system dimensionality and its dynamical characteristics, numerous techniques have been developed to quantify the response statistics, i.e. the probability density function (pdf) for the system state. For systems subjected to white noise, FokkerPlanckKolmogorov (FPK) equation provides a complete statistical description of the response statistics [7, 8, 9]. However, exact analytical solutions of the FPK equation are available only for a small class of systems. An alternative computational approach, the path integral solution (PIS) method, has been developed to provide the response pdf for general nonlinear systems at a specific time instant given the pdf of an earlier time instant. Many studies have been focused on the application of stepbystep PIS method numerically [10, 11, 12] and analytically [13, 14, 15] reporting its effectiveness on capturing the response statistics. On the other hand, for nonMarkovian systems subjected to correlated excitations the joint responseexcitation pdf method provides a computational framework for the full statistical solution [16, 17, 18]. However, such methodologies rely on the solution of transport equations for the pdf and they are associated with very high computational cost especially when it comes to the optimization of system parameters.
To avoid solving the transport equations for the pdf, semianalytical approximative approaches with significantly reduced computational cost have been developed. Among them the most popular method in the context of structural systems is the statistical linearization method [19, 20, 21, 22, 23], which can also handle correlated excitations. The basic concept of this approach is to replace the original nonlinear equation of motion with a linear equation, which can be treated analytically, by minimizing the statistical difference between those two equations. Statistical linearization performs very well for systems with unimodal statistics, i.e. close to Gaussian. However, when the response is essentially nonlinear, e.g. as it is the case for a doublewell oscillator, the application of statistical linearization is less straightforward and involves the adhoc selection of shape parameters for the response statistics [24].
An alternative class of methods relies on the derivation of moment equations, which describes the evolution of the the joint responseexcitation statistical moments or (depending on the nature of the stochastic excitation) the response statistical moments [25, 26, 27]. The challenge with moment equations arises if the equation of motion of the system contains nonlinear terms in which case we have the well known closure problem. This requires the adoption of closure schemes, which essentially truncate the infinite system of moment equations to a finite one. The simplest approach along this line is the Gaussian closure [28] but nonlinear closure schemes have also been developed (see e.g. [29, 30, 31, 32, 33, 34, 35, 36, 37]). In most cases, these nonlinear approaches may offer some improvement compared with the stochastic linearization approach applied to nonlinear systems but the associated computational cost is considerably larger [38]. For strongly nonlinear systems, such as bistable systems, these improvements can be very small. Bistable systems, whose potential functions have bimodal shapes, have become very popular in energy harvesting applications [39, 40, 41, 42, 43, 44, 45, 46], where there is a need for fast and reliable calculations that will be able to resolve the underlying nonlinear dynamics in order to provide with optimal parameters of operation (see e.g. [47, 48]).
The goal of this work is the development of a closure methodology that can overcome the limitations of traditional closure schemes and can approximate the steady state statistical structure of bistable systems excited by correlated noise. We first formulate the moment equations for the joint pdf of the response and the excitation at two arbitrary time instants [49]. To close the resulted system of moment equations, we formulate a twotime representation of the joint responseexcitation pdf using copula functions. We choose the representation so that the single time statistics are consistent in form with the FokkerPlanckKolmogorov solution in steady state, while the joint statistical structure between two different time instants is represented with a Gaussian copula density. Based on these two ingredients (dynamical information expressed as moment equations and assumed form of the response statistics), we formulate a minimization problem with respect to the unknown parameters of the pdf representation so that both the moment equations and the closure induced by the representation are optimally satisfied. For the case of unimodal systems, the described approach reproduces the statistical linearization method while for bimodal systems it still provides meaningful and accurate results with very low computational cost.
The developed approach allows for the inexpensive and accurate approximation of the second order statistics of the system even for oscillators associated with doublewell potentials. In addition, it allows for the semianalytical approximation of the full nonGaussian joint responseexcitation pdf in a postprocessing manner. We illustrate the developed approach through nonlinear singledegreeoffreedom energy harvesters with doublewell potentials subjected to correlated noise with PiersonMoskowitz power spectral density. We also consider the case of bistable oscillators coupled with electromechanical energy harvesters (one and a half degreesoffreedom systems), and we demonstrate how the proposed probabilistic framework can be used for performance optimization and parameters selection.
2 Description of the Method
In this section, we give a detailed description of the proposed method for the inexpensive computation of the response statistics for dynamical systems subjected to colored noise excitation. The computational approach relies on two basic ingredients:

Twotime statistical moment equations. These equations will be derived directly from the system equation and they will express the dynamics that govern the twotime statistics. For systems excited by whitenoise, single time statistics are sufficient to describe the response but for correlated excitation, this is not the case and it is essential to consider higher order moments. Note that higher (than two) order statistical moment equations may be used but in the context of this work twotime statistics would be sufficient.

Probability density function (pdf) representation for the joint responseexcitation statistics. This will be a family of probability density functions with embedded statistical properties such as multimodality, tail decay properties, correlation structure between response and excitation, or others. The joint statistical structure will be represented using copula functions. We will use representations inspired by the analytical solutions of the dynamical system when this is excited by white noise. These representations will reflect features of the Hamiltonian structure of the system and will be used to derive appropriate closure schemes that will be combined with the moment equations.
Based on these two ingredients, we will formulate a minimization problem with respect to the unknown parameters of the pdf representation so that both the moment equations and the closure induced by the representation are optimally satisfied. We will see that for the case of unimodal systems the described approach reproduces the statistical linearization method while for bimodal systems it still provides meaningful and accurate results with very low computational cost.
For the sake of simplicity, we will present our method through a specific system involving a nonlinear SDOF oscillator with a double well potential. This system has been studied extensively in the context of energy harvesting especially for the case of white noise excitation [41, 50, 51, 52]. However, for realistic setups it is important to be able to optimize/predict its statistical properties under general (colored) excitation. More specifically we consider a nonlinear harvester of the form
(1) 
where is the relative displacement between the harvester mass and the base, is the base excitation representing a stationary stochastic process, is normalized (with respect to mass) damping coefficient, and and are normalized stiffness coefficients.
2.1 Twotime Moment System
We consider two generic time instants, and . The twotime moment equations have been considered previously in [49] for the determination of the solution of a ‘half’ degreeoffreedom nonlinear oscillator by utilizing a Gaussian closure. We multiply the equation of motion at time with the response displacement and apply the mean value operator (ensemble average). This will give us an equation which contains an unknown term on the right hand side. To determine this term we repeat the same step but we multiply the equation of motion with . This gives us the following twotime moment equations:
(2)  
(3) 
Here the excitation is assumed to be a stationary stochastic process with zero mean and a given power spectral density; this can have an arbitrary form, e.g. monochromatic, colored, or white noise. Since the system is characterized by an odd restoring force, we expect that its response also has zero mean. Moreover, we assume that after an initial transient the system will be reaching a statistical steady state given the stationary character of the excitation. Based on properties of mean square calculus [3, 27], we interchange the differentiation and the mean value operators. Then the moment equations will take the form:
(4)  
(5) 
Expressing everything in terms of the covariance functions, above equations will result in:
(6)  
(7) 
where the covariance function is defined as
(8) 
Taking into account the assumption for a stationary response (after the system has gone through an initial transient phase), the above moment equations can be rewritten in terms of the time difference :
(9)  
(10) 
Note that all the linear terms in the original equation of motion are expressed in terms of covariance functions, while the nonlinear (cubic) terms show up in the form of fourth order moments. To compute the latter we will need to adopt an appropriate closure scheme.
2.2 Twotime PDF representations and induced closures
In the absence of higherthantwo order moments, the response statistics can be analytically obtained in a straightforward manner. However, for higher order terms it is necessary to adopt an appropriate closure scheme that closes the infinite system of moment equations. A standard approach in this case, which performs very well for unimodal systems, is the application of Gaussian closure which utilizes Isserlis’ Theorem [53] to connect the higher order moments with the second order statistical quantities. Despite its success for unimodal systems, Gaussian closure does not provide accurate results for bistable systems. This is because in this case (i.e. bistable oscillators) the closure induced by the Gaussian assumption does not reflect the properties of the system attractor in the statistical steady state.
Here we aim to solve this problem by proposing a nonGaussian representation for the joint responseresponse pdf at two different time instants and for the joint responseexcitation pdf at two different time instants. These representations will:

incorporate specific properties or information about the response pdf (single time statistics) in the statistical steady state,

capture the correlation structure between the statistics of the response and/or excitation at different time instants by employing Gaussian copula density functions,

have a consistent marginal with the excitation pdf (for the case of the joint responseexcitation pdf).
2.2.1 Representation Properties for Single Time Statistics
We begin by introducing the pdf properties for the single time statistics. The selected representation will be based on the analytical solutions of the FokkerPlanck equation which are available for the case of white noise excitation [54, 3], and for vibrational systems that has an underlying Hamiltonian structure. Here we will leave the energy level of the system as a free parameter  this will be determined later. In particular, we will consider the following family of pdf solutions (Figure 2a):
(11) 
where is the potential energy of the oscillator, is a free parameter connected with the energy level of the system, and is the normalization constant expressed as follows:
(12) 
2.2.2 Correlation Structure between Twotime Statistics
Representing the single time statistics is not sufficient since for nonMarkovian systems (i.e. correlated excitation) the system dynamics can be effectively expressed only through (at least) twotime statistics. To represent the correlation between two different time instants we introduce Gaussian copula densities [55, 56]. A copula is a multivariate probability distribution with uniform marginals. It has emerged as an useful tool for modeling stochastic dependencies allowing the separation of dependence modeling from the given marginals [57]. Based on this formulation we obtain pdf representations for the joint responseresponse and responseexcitation at different time instants.
Joint responseexcitation pdf. We first formulate the joint responseexcitation pdf at two different (arbitrary) time instants. In order to design the joint pdf based on the given marginals of response and excitation, we utilize a bivariate Gaussian copula whose density can be written as follows [56]:
(13) 
where and indicate cumulative distribution functions and the standard cumulative distribution function is given as the following form:
(14) 
Denoting with the argument that corresponds to the response at time , with the argument for the excitation at time , and with the (zeromean) marginal pdf for the excitation, we have the expression for the joint responseexcitation pdf.
(15) 
where defines the correlation between the response and the excitation and has values and and are the cumulative distribution functions obtained through the response marginal pdf, , and the excitation marginal pdf, , respectively. Note that the coefficient depends on the time difference of the response and excitation. This dependence will be recovered through the resolved secondorder moments (over time) between the response and excitation.
Joint responseresponse pdf. The joint pdf for two different time instants of the response, denoted as , is a special case of what has been presented. In order to avoid confusion, a different notation is used to represent the response at a different time instant . We have:
(16) 
where is a correlation constant (that depends on the timedifference ). Note that the response at the second time instant follows the same nonGaussian pdf corresponding to the single time statistics of the response.
In Figure 2b, we present the above joint pdf (15) with the marginal (response) having a bimodal structure and the marginal (excitation) having a Gaussian structure. For we have independence, which essentially expresses the case of very distant twotime statistics, while as we increase the correlation between the two variables increases referring to the case of small values of .
2.2.3 Induced NonGaussian Closures
Using these nonGaussian pdf representations, we will approximate the fourth order moment terms that show up in the moment equations. We numerically observe that in the context of the pdf representations given above, the relation between and is essentially linear (see Figure 3). To this end, we choose a closure of the following form for both the responseresponse and the responseexcitation terms:
(17) 
where is the closure coefficient for the joint responseresponse statistics. The value of is obtained by expanding both and with respect to keeping up to the first order terms:
(18)  
(19) 
where the error function is given by:
(20) 
Thus, we observe that the assumed copula function in combination with the marginal densities prescribe an explicit dependence between fourth and secondorder moments, expressed through the coefficient:
(21) 
We emphasize that this closure coefficient does not depend on the timedifference but only on the single time statistics and in particular the energy level of the system, defined by . To this end, for any given marginal pdf , we can analytically find what would be the closure coefficient under the assumptions of the adopted copula function.
The corresponding coefficient for the joint responseexcitation statistics can be similarly obtained through a first order expansion of the moments:
(22) 
The closure coefficient has exactly the same form with the closure coefficient and it does not depend on the statistical properties of the excitation nor on the timedifference but only on the energy level . We will refer to equations (21) and (22) as the closure constraints. This will be one of the two sets of constraints that we will include in the minimization procedure for the determination of the solution.
2.2.4 Closed Moment Equations
The next step involves the application of above closure scheme on the derived twotime moment equations. By directly applying the induced closure schemes on equations (9) and (10), we have the linear set of moment equations for the secondorder statistics:
(23)  
(24) 
Using the WienerKhinchin theorem, we transform the above equations to the corresponding power spectral density equations:
(25)  
(26) 
These equations allow us to obtain an expression for the power spectral density of the response displacement in terms of the excitation spectrum:
(27) 
Integration of the above equation will give us the variance of the response:
(28) 
The last equation is the second constraint, the dynamics constraint, which expresses the second order dynamics of the system. Our goal is to optimally satisfy it together with the closure constraints defined by equations (21) and (22).
2.2.5 Moment Equation Copula Closure (MECC) Method
The last step is the minimization of the two set of constraints, the closure constraints and the dynamics constraint, which have been expressed in terms of the system response variance . The minimization will be done in terms of the unknown energy level and the closure coefficients and .
More specifically, we define the following cost function which incorporates our constraints:
(29) 
Note that in the context of statistical linearization only the first constraint is minimized while the closure coefficient is the one that follows exactly from a Gaussian representation for the pdf. In this context there is no attempt to incorporate in an equal manner the mismatch in the dynamics and the pdf representation. The minimization of this cost function essentially allows mismatch for the equation (expressed through the dynamic constraint) but also for the pdf representation (expressed through the closure constraints). For linear systems and an adopted Gaussian pdf for the response the above cost function vanishes identically.
3 Applications to Bistable Energy Harvesters under Correlated Excitation
We apply the presented Moment Equation Copula Closure (MECC) method to two nonlinear vibrational systems. One is a singledegreeoffreedom (SDOF) bistable oscillator with linear damping that simulates energy harvesting while the other is a similar bistable oscillator coupled with an electromechanical energy harvester. For both applications, it is assumed that the stationary stochastic excitation has a power spectral density given by the PiersonMoskowitz spectrum, which is typical for excitation created by random water waves:
(30) 
where controls the intensity of the excitation.
3.1 SDOF Bistable Oscillator Excited by Colored Noise
For the colored noise excitation that we just described, we apply the MECC method. We consider a set of system parameters that correspond to a double well potential. Depending on the intensity of the excitation (which is adjusted by the factor ), the response of the bistable system ‘lives’ in three possible regimes. If is very low, the bistable system is trapped in either of the two wells while if is very high the energy level is above the homoclinic orbit and the system performs crosswell oscillations. Between these two extreme regimes, the stochastic response exhibits combined features and characteristics of both energy levels and it has a highly nonlinear, multifrequency character [58, 59].
Despite these challenges, the presented MECC method can inexpensively provide with a very good approximation of the system’s statistical characteristics as it is shown in Figure 4. In particular in Figure 4, we present the response variance as the intensity of the excitation varies for two sets of the system parameters. We also compare our results with direct MonteCarlo simulations and with a standard Gaussian closure method [3, 4, 1].
For the MonteCarlo simulations the time series for the excitation has been generated as the sum of cosines over a range of frequencies. The amplitudes and the range of frequencies are determined through the power spectrum while the phases are assumed to be random variables which follow a uniform distribution. In the presented examples, the excitation has power spectral density that follows the PiersonMoskowitz spectrum. Once each ensemble time series for the excitation has been computed, the governing ordinary differential equation is solved using a 4th/5th order RungeKutta method. For each realization the system is integrated for a sufficiently long time interval in order to guarantee that the response statistics have converged. For each problem, we generate 100 realizations in order to compute the secondorder statistics. However, for the computation of the full joint pdf, a significantly larger number of samples is needed reaching the order of .
We observe that for very large values of the computed approximation closely follows the MonteCarlo simulation. On the other hand, the Gaussian closure method systematically underestimates the variance of the response. For lower intensities of the excitation, the exact (MonteCarlo) variance presents a nonmonotonic behavior with respect to due to the coexistence of the cross and intrawell oscillations. While the Gaussian closure has very poor performance on capturing this trend, the MECC method can still provide a satisfactory approximation of the dynamics. Note that the nonsmooth transition observed in the MECC curve is due to the fact that for very low values of the minimization of the cost function (29) does not reach a zero value while this is the case for larger values of . In other words, in the strongly nonlinear regime neither the dynamics constraint nor the closure constraint is satisfied exactly, yet this optimal solution provides with a good approximation of the system dynamics.
After we have obtained the unknown parameters , and by minimizing the cost function for each given , we can then compute the covariance functions and the joint pdf in a postprocess manner. More specifically, since a known corresponds to a specific (equation (22)) we can immediately determine by taking the inverse Fourier transform of found through equation (25). The next step is the numerical integration of the closed moment equation (24) utilizing the determined value with initial conditions given by
(31) 
where the second condition follows from the symmetry properties of . Note that we integrate equation (24) instead of using the inverse Fourier transform as we did for so that we can impose the variance found in the last equation by integrating the resulted density for the determined . Using the correlation functions and we can also determine, for each case, the correlation coefficient of the copula function for each timedifference . The detailed steps are given at the end of this subsection.
The results as well as a comparison with the Gaussian closure method and a direct MonteCarlo simulation are presented in Figure 5. We can observe that through the proposed approach we are able to satisfactorily approximate the correlation function even close to the nonlinear regime , where the Gaussian closure method presents important discrepancies.
Finally, using the computed parameters and closure coefficients, and , we can also construct the threedimensional nonGaussian joint pdf for the responseresponseexcitation at different time instants. This will be derived based on the threedimensional Gaussian copula density of the following form:
(32) 
The threedimensional nonGaussian joint pdf for the responseresponseexcitation at different time instants can be expressed as follows:
(33) 
where represents the correlation matrix with all diagonal elements equal to 1:
(34) 
The time dependent parameters , , of the copula function can be found through the resolved moments, by expanding the latter as:
(35)  
(36)  
(37) 
where,
If necessary higher order terms may be retained in the Taylor expansion although for the present problem a linear approximation was sufficient. The computed approximation is presented in Figure 6 through two dimensional marginals as well as through isosurfaces of the full threedimensional joint pdf. We compare with direct MonteCarlo simulations and as we are able to observe, the computed pdf compares favorably with the expensive MonteCarlo simulation. The joint statistics using the MonteCarlo approach were computed using number of samples while the computational cost of the MECC method involved the minimization of a three dimensional function.
3.2 SDOF bistable oscillator coupled to an electromechanical harvester
In practical configurations, energy harvesting occurs through a linear electromechanical transducer coupled to the nonlinear oscillator [42, 60, 61]. In this section, we assess how our method performs for a bistable nonlinear SDOF oscillator coupled to a linear electromechanical transducer. The equations of motion in this case take the form:
(38)  
(39) 
where is the response displacement, is a stationary stochastic excitation, is the voltage across the load, is the normalized damping coefficient, and are the normalized stiffness coefficients, and are the normalized coupling coefficients, and is the normalized time coefficient for the electrical system. All the coefficients except are positive. Based on the linearity of the second equation, we express the voltage in an integral form:
(40) 
where indicates convolution and represents the Heaviside step function. We then formulate the secondorder moment equations following a similar approach with the previous section.
(41)  
(42)  
(43) 
In this case, we estimate two additional covariance functions, and before applying MECC method:
(44)  
(45)  
(46)  
(47) 
where is the time difference of two generic time instants and . Considering the power spectrum, the Fourier transform of the above gives:
(48) 
Similarly, we also obtain for :
(49) 
By applying the previously described closure scheme on equations (41) and (42), we have a linear set of moment equations for the secondorder statistics:
(50)  
(51)  
(52) 
Using the WienerKhinchin theorem, we transform the above equations to the corresponding power spectral density equations:
(53)  
(54)  
(55) 
These equations allow us to obtain an expression for the power spectral density of the response displacement and response voltage in terms of the excitation spectrum:
(56)  
(57) 
Integration of the above equation will give us the variance of the response displacement and voltage:
(58)  
(59) 
Equation (58) expresses the second order dynamics of the SDOF bistable oscillator coupled with an electromechanical harvester, and is the dynamics constraint for this system. We will minimize it together with the closure constraints defined by equations (21) and (22):
(60) 
In Figure 7, we illustrate the variance of the response displacement and the voltage as the intensity of the excitation varies for two sets of the system parameters. For both sets of system parameters, we observe that for large intensity of the excitation, the MECC method computes the response variances (displacement and voltage) very accurately, while the Gaussian closure method systematically underestimates them. For lower intensities of the excitation, the response displacement variance computed by the MonteCarlo simulation presents a nonmonotonic behavior with respect to . While the Gaussian closure has very poor performance on capturing this trend, the MECC method can still provide a satisfactory approximation of the dynamics.
Following similar steps with the previous section, we obtain the covariance functions of the response displacement and voltage and the joint pdf in a postprocess manner. The results as well as a comparison with the Gaussian closure method and the MonteCarlo simulation are illustrated in Figure 8. We can observe that through the proposed approach we are able to satisfactorily approximate the correlation function even close to the nonlinear regime , where the Gaussian closure method presents important discrepancies. In Figure 9, we illustrate two dimensional marginal pdfs as well as isosurfaces of the full threedimensional joint pdf. We compare with direct MonteCarlo simulations and as we are able to observe, the computed pdf closely approximates the expensive MonteCarlo simulation in statistical regimes which are far from Gaussian.
Finally in Figure 10, we demonstrate how the proposed MECC method can be used to study robustness over variations of the excitation parameters. In particular, we present the mean square response displacement and response voltage estimated for various amplification factors and frequencyvaried excitation spectra:
(61) 
where is the perturbation frequency. The comparison with direct MonteCarlo simulation indicates the effectiveness of the presented method to capture accurately the response characteristics over a wide range of input parameters.
4 Conclusions
We have considered the problem of determining the nonGaussian steady state statistical structure of bistable nonlinear vibrational systems subjected to colored noise excitation. We first derived moment equations that describe the dynamics governing the twotime statistics. We then combined those with a nonGaussian pdf representation for the joint responseresponse and joint responseexcitation statistics. This representation has i) single time statistical structure consistent with the analytical solutions of the FokkerPlanck equation, and ii) twotime statistical structure that follows from the adoption of a Gaussian copula function. The pdf representation takes the form of closure constraints while the moment equations have the form of a dynamics constraint. We formulated the two sets of constraints as a lowdimensional minimization problem with respect to the unknown parameters of the representation. The minimization of both the dynamics constraint and the closure constraints imposes an interplay between these two factors.
We then applied the presented method to two nonlinear oscillators in the context of vibration energy harvesting. One is a single degree of freedom (SDOF) bistable oscillator with linear damping while the other is a same SDOF bistable oscillator coupled with an electromechanical energy harvester. For both applications, it was assumed that the stationary stochastic excitation has a power spectral density given by the PiersonMoskowitz spectrum. We have shown that the presented method can provide a very good approximation of second order statistics of the system, when compared with direct MonteCarlo simulations, even in essentially nonlinear regimes, where Gaussian closure techniques fail completely to capture the dynamics. In addition, we can compute the full (nonGaussian) probabilistic structure of the solution in a postprocess manner. We emphasize that the computational cost associated with the new method is considerably smaller compared with methods that evolve the pdf of the solution since MECC method relies on the minimization of a function with a few unknown variables.
These results indicate that the new method can be a very good candidate when it comes to the calculation of the stochastic response for vibrational system with complex potentials as it is required in parameter optimization or selection. Future endeavors include the application of the presented approach in higher dimensional contexts involving nonlinear energy harvesters and passive protection of structures as well as on the development/optimization of structural configurations able to operate effectively under intermittent loads [62].
Acknowledgments
We would like to thank the anonymous referees for helpful suggestions and in particular one of them for suggesting the use of copula functions, a change that significantly improved the original version of the manuscript. We would also like to acknowledge the support from the Samsung Scholarship Program as well as the MIT Energy Initiative for support under the grant ‘Nonlinear Energy Harvesting From BroadBand Vibrational Sources By Mimicking Turbulent Energy Transfer Mechanisms’. T.P.S. is also grateful to the American Bureau of Shipping for support under a Career Development Chair.
References
 [1] M. Grigoriu, Stochastic calculus: applications in science and engineering. Springer, 2002.
 [2] R. Stratonovich, Topics in the theory of random noise, vol. 2. CRC Press, 1967.
 [3] K. Sobczyk, Stochastic differential equations: with applications to physics and engineering, vol. 40. Springer, 2001.
 [4] T. Soong and M. Grigoriu, “Random vibration of mechanical and structural systems,” NASA STI/Recon Technical Report A, vol. 93, p. 14690, 1993.
 [5] A. Naess and T. Moan, Stochastic dynamics of marine structures. Cambridge University Press, 2012.
 [6] C. To, Nonlinear random vibration: Analytical techniques and applications. CRC Press, 2011.
 [7] S. Wojtkiewicz, E. Johnson, L. Bergman, M. Grigoriu, and B. Spencer Jr, “Response of stochastic dynamical systems driven by additive gaussian and poisson white noise: Solution of a forward generalized kolmogorov equation by a spectral finite difference method,” Computer methods in applied mechanics and engineering, vol. 168, no. 1, pp. 73–89, 1999.
 [8] J. Dunne and M. Ghanbari, “Extremevalue prediction for nonlinear stochastic oscillators via numerical solutions of the stationary fpk equation,” Journal of Sound and Vibration, vol. 206, no. 5, pp. 697–724, 1997.
 [9] M. Di Paola and A. Sofi, “Approximate solution of the fokker–planck–kolmogorov equation,” Probabilistic Engineering Mechanics, vol. 17, no. 4, pp. 369–384, 2002.
 [10] M. F. Wehner and W. Wolfer, “Numerical evaluation of pathintegral solutions to fokkerplanck equations,” Physical Review A, vol. 27, no. 5, p. 2663, 1983.
 [11] A. Naess and J. Johnsen, “Response statistics of nonlinear, compliant offshore structures by the path integral solution method,” Probabilistic Engineering Mechanics, vol. 8, no. 2, pp. 91–106, 1993.
 [12] M. Di Paola and R. Santoro, “Path integral solution for nonlinear system enforced by poisson white noise,” Probabilistic Engineering Mechanics, vol. 23, no. 2, pp. 164–169, 2008.
 [13] I. Kougioumtzoglou and P. Spanos, “An analytical wiener path integral technique for nonstationary response determination of nonlinear oscillators,” Probabilistic Engineering Mechanics, vol. 28, pp. 125–131, 2012.
 [14] I. A. Kougioumtzoglou and P. D. Spanos, “Nonstationary stochastic response determination of nonlinear systems: A wiener path integral formalism,” Journal of Engineering Mechanics, vol. 140, no. 9, p. 04014064, 2014.
 [15] A. Di Matteo, I. A. Kougioumtzoglou, A. Pirrotta, P. D. Spanos, and M. Di Paola, “Stochastic response determination of nonlinear oscillators with fractional derivatives elements via the wiener path integral,” Probabilistic Engineering Mechanics, vol. 38, pp. 127–135, 2014.
 [16] T. P. Sapsis and G. A. Athanassoulis, “New partial differential equations governing the joint, response–excitation, probability distributions of nonlinear systems, under general stochastic excitation,” Probabilistic Engineering Mechanics, vol. 23, no. 2, pp. 289–306, 2008.
 [17] D. Venturi, T. P. Sapsis, H. Cho, and G. E. Karniadakis, “A computable evolution equation for the joint responseexcitation probability density function of stochastic dynamical systems,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, vol. 468, no. 2139, pp. 759–783, 2012.
 [18] H. Cho, D. Venturi, and G. E. Karniadakis, “Adaptive discontinuous galerkin method for responseexcitation pdf equations,” SIAM Journal on Scientific Computing, vol. 35, no. 4, pp. B890–B911, 2013.
 [19] T. Caughey, “Response of a nonlinear string to random loading,” Journal of Applied Mechanics, vol. 26, no. 3, pp. 341–344, 1959.
 [20] T. Caughey, “Equivalent linearization techniques,” The Journal of the Acoustical Society of America, vol. 35, no. 11, pp. 1706–1711, 1963.
 [21] I. Kazakov, “An approximate method for the statistical investigation of nonlinear systems,” Trudy VVIA im Prof. NE Zhukovskogo, vol. 394, pp. 1–52, 1954.
 [22] J. Roberts and P. Spanos, Random vibration and statistical linearization. Courier Dover Publications, 2003.
 [23] L. Socha, Linearization methods for stochastic dynamic systems, vol. 730. Springer, 2008.
 [24] S. H. Crandall, “On using nongaussian distributions to perform statistical linearization,” vol. 39, p. 1395, 2004.
 [25] N. Sancho, “Technique for finding the moment equations of a nonlinear stochastic system,” Journal of Mathematical Physics, vol. 11, no. 3, pp. 771–774, 1970.
 [26] D. Bover, “Moment equation methods for nonlinear stochastic systems,” Journal of Mathematical Analysis and Applications, vol. 65, no. 2, pp. 306–320, 1978.
 [27] J. Beran, Statistics for longmemory processes, vol. 61. CRC Press, 1994.
 [28] R. Iyengar and P. Dash, “Study of the random vibration of nonlinear systems by the gaussian closure technique,” Journal of Applied Mechanics, vol. 45, no. 2, pp. 393–399, 1978.
 [29] S. Crandall, “Nongaussian closure for random vibration of nonlinear oscillators,” International Journal of NonLinear Mechanics, vol. 15, no. 4, pp. 303–313, 1980.
 [30] S. Crandall, “Nongaussian closure techniques for stationary random vibration,” International journal of nonlinear mechanics, vol. 20, no. 1, pp. 1–8, 1985.
 [31] Q. Liu and H. Davies, “Application of nongaussian closure to the nonstationary response of a duffing oscillator,” International journal of nonlinear mechanics, vol. 23, no. 3, pp. 241–250, 1988.
 [32] W. Wu and Y. Lin, “Cumulantneglect closure for nonlinear oscillators under random parametric and external excitations,” International Journal of NonLinear Mechanics, vol. 19, no. 4, pp. 349–362, 1984.
 [33] R. Ibrahim, A. Soundararajan, and H. Heo, “Stochastic response of nonlinear dynamic systems based on a nongaussian closure,” Journal of applied mechanics, vol. 52, no. 4, pp. 965–970, 1985.
 [34] M. Grigoriu, “A consistent closure method for nonlinear random vibration,” International journal of nonlinear mechanics, vol. 26, no. 6, pp. 857–866, 1991.
 [35] A. Hasofer and M. Grigoriu, “A new perspective on the moment closure method,” Journal of applied mechanics, vol. 62, no. 2, pp. 527–532, 1995.
 [36] S. Wojtkiewicz, B. Spencer Jr, and L. Bergman, “On the cumulantneglect closure method in stochastic dynamics,” International journal of nonlinear mechanics, vol. 31, no. 5, pp. 657–684, 1996.
 [37] M. Grigoriu, “Moment closure by monte carlo simulation and moment sensitivity factors,” International journal of nonlinear mechanics, vol. 34, no. 4, pp. 739–748, 1999.
 [38] M. Noori, A. Saffar, and H. Davoodi, “A comparison between nongaussian closure and statistical linearization techniques for random vibration of a nonlinear oscillator,” Computers & structures, vol. 26, no. 6, pp. 925–931, 1987.
 [39] P. Green, K. Worden, K. Atallah, and N. Sims, “The benefits of duffingtype nonlinearities and electrical optimisation of a monostable energy harvester under white gaussian excitations,” Journal of Sound and Vibration, vol. 331, no. 20, pp. 4504–4517, 2012.
 [40] R. Harne and K. Wang, “A review of the recent research on vibration energy harvesting via bistable systems,” Smart Materials and Structures, vol. 22, no. 2, p. 023001, 2013.
 [41] M. Daqaq, “Transduction of a bistable inductive generator driven by white and exponentially correlated gaussian noise,” Journal of Sound and Vibration, vol. 330, no. 11, pp. 2554–2564, 2011.
 [42] E. Halvorsen, “Fundamental issues in nonlinear widebandvibration energy harvesting,” Physical Review E, vol. 87, no. 4, p. 042129, 2013.
 [43] P. Green, E. Papatheou, and N. Sims, “Energy harvesting from human motion and bridge vibrations: An evaluation of current nonlinear energy harvesting solutions,” Journal of Intelligent Material Systems and Structures, 2013.
 [44] Q. He and M. Daqaq, “New insights into utilizing bistability for energy harvesting under white noise,” Journal of Vibration and Acoustics, 2014.
 [45] B. Mann and N. Sims, “Energy harvesting from the nonlinear oscillations of magnetic levitation,” Journal of Sound and Vibration, vol. 319, no. 1, pp. 515–530, 2009.
 [46] D. Barton, S. Burrow, and L. Clare, “Energy harvesting from vibrations with a nonlinear oscillator,” Journal of Vibration and Acoustics, vol. 132, no. 2, p. 021009, 2010.
 [47] H. K. Joo and T. P. Sapsis, “Performance measures for singledegreeoffreedom energy harvesters under stochastic excitation,” Journal of Sound and Vibration, vol. 333, no. 19, pp. 4695–4710, 2014.
 [48] J. M. Kluger, T. P. Sapsis, and A. H. Slocum, “Robust energy harvesting from walking vibrations by means of nonlinear cantilever beams,” Journal of Sound and Vibration, vol. 341, pp. 174–194, 2015.
 [49] G. Athanassoulis, I. Tsantili, and Z. Kapelonis, “Twotime, responseexcitation moment equations for a cubic halfoscillator under gaussian and cubicgaussian colored excitation. part 1: The monostable case,” arXiv preprint arXiv:1304.2195, 2013.
 [50] M. Daqaq, “On intentional introduction of stiffness nonlinearities for energy harvesting under white gaussian excitations,” Nonlinear Dynamics, vol. 69, no. 3, pp. 1063–1079, 2012.
 [51] L. Gammaitoni, I. Neri, and H. Vocca, “Nonlinear oscillators for vibration energy harvesting,” Applied Physics Letters, vol. 94, no. 16, p. 164102, 2009.
 [52] M. Ferrari, V. Ferrari, M. Guizzetti, B. Andò, S. Baglio, and C. Trigona, “Improved energy harvesting from wideband vibrations by nonlinear piezoelectric converters,” Sensors and Actuators A: Physical, vol. 162, no. 2, pp. 425–431, 2010.
 [53] L. Isserlis, “On a formula for the productmoment coefficient of any order of a normal frequency distribution in any number of variables,” Biometrika, pp. 134–139, 1918.
 [54] C. Soize, The FokkerPlanck equation for stochastic dynamical systems and its explicit steady state solutions, vol. 17. World Scientific, 1994.
 [55] R. B. Nelsen, An introduction to copulas. Springer Science & Business Media, 2007.
 [56] C. Meyer, “The bivariate normal copula,” Communications in StatisticsTheory and Methods, vol. 42, no. 13, pp. 2402–2422, 2013.
 [57] L. Qu and W. Yin, “Copula density estimation by total variation penalized likelihood with linear equality constraints,” Computational Statistics & Data Analysis, vol. 56, no. 2, pp. 384–398, 2012.
 [58] M. Dykman, S. Soskin, and M. Krivoglaz, “Spectral distribution of a nonlinear oscillator performing brownian motion in a doublewell potential,” Physica A: Statistical Mechanics and its Applications, vol. 133, no. 1, pp. 53–73, 1985.
 [59] M. Dykman, R. Mannella, R. McClintock, F. Moss, and S. Soskin, “Spectral density of fluctuations of a doublewell duffing oscillator driven by white noise,” Physical Review A, vol. 37, no. 4, p. 1303, 1988.
 [60] M. Karami and D. Inman, “Powering pacemakers from heartbeat vibrations using linear and nonlinear energy harvesters,” Applied Physics Letters, vol. 100, no. 4, p. 042901, 2012.
 [61] R. Masana and M. Daqaq, “Response of duffingtype harvesters to bandlimited noise,” Journal of Sound and Vibration, vol. 332, no. 25, pp. 6755–6767, 2013.
 [62] M. Mohamad and T. Sapsis, “Probabilistic description of extreme events in intermittently unstable dynamical systems excited by correlated stochastic processes,” SIAM/ASA Journal on Uncertainty Quantification, vol. 3, no. 1, pp. 709–736, 2015.