Abstract
This paper proposes an efficient short-time probability approximation with Lévy excitation to capture the transient probability distribution and its evolving path. Using principal component analysis (PCA), the method constructs a probability core to exclude outliers beyond it. The statistics of samples that fall inside the core are treated, with a prescribed fiducial probability, as an easy-to-estimate Gaussian type. The idea is verified numerically by compared with Monte-Carlo results. Then, it is integrated into the path integral (PI) method, combined with evolving probabilistic vector (EPV) techniques, to efficiently obtain probability distributions in each time step of PI. This scheme is semianalytical, only dependent on a relatively small amount of response samples to form the probability core; thus, it can have very computational advantages over full Monte-Carlo simulation to capture transient responses and probability distributions. The application to investigating response transitions of a nonsmooth system driven by Lévy shock and jump has revealed the performance of the proposed method. Also, the exit times of stochastic response are characterized quantitatively from the perspective of global dynamic transition. These investigations will be helpful to achieve the efficient probability estimation for nonlinear system with non-Gaussian inputs and quantify the reliability of the mechanical system.
1. Introduction
Noise, which is ubiquitous in engineering, inevitably disturbs the dynamics of mechanical systems. Non-Gaussian Lévy noise with different relevant parameters is a more general form to describe various types of specific disturbances such as intermittent shock and jump, which has been observed in sea wave motion [1], turbulent fluid flows [2], plasmas [3], and heartbeat dynamics [4]. The uncertainty, applied to a strong nonlinear system like a nonsmooth dynamical structure, may arouse some potential dynamic behaviors underlying intrinsic global structure to induce a switch or even a transition between multiple responses. For a nonsmooth system, of course, the change of dynamics might be intentional design (such as a stopper in an energy harvesting device [5, 6]) or undesired trouble (such as clearances in a gear system [7, 8]). Anyway, plenty of studies and their results have uncovered the fact that the noise, under a basic frame of global structure, dominates the final forms of responses. Therefore, it is of significance to develop an efficient method for investigating the disturbance effects of noise on dynamical responses, particularly to determine when the response of the system likely escapes from current motion and where it will go from an engineering standpoint, which can gain insight into predictions of fatigue life and reliability, as well as structural design of a complex nonlinear system.
Noise-induced behaviors could trigger a large-range transition in the transient process [9–11]. Generally, the probabilistic description of the stochastic response is governed by Fokker–Planck–Kolmogorov (FPK) equation, but it is difficult to obtain the exact solution for the second-order parabolic partial differential equation. The most direct methods to estimate the random characteristics of the response are the statistical approaches based on the Monte-Carlo simulation (MCS), but it exhibits the widely known drawback that the accuracy of the estimates depends on the number of samples, increasing the related computational effort. Some alternative approximation methods consider the solving input–output relations of probability density function (PDF), for example, cumulant-neglect closure schemes [12] based on evolutive moments, but it is efficient only when the responses process is Gaussian. The stochastic perturbation methods [13] based on Taylor expansion can be applied to formulate the statistics’ governing equations with uncertainties, but restricted to local dynamics. The probability transformation method [14] and its variants [15, 16] have yielded much in describing the response probabilistic characterizations of linear structures with non-Gaussian input process, but it still faces a computational challenge with respect to the evolution of PDF, above all for large nonlinear systems. It is well known that the path integral method in terms of the Markov rule can work nicely to numerically estimate the probabilistic evolution of responses. Nevertheless, the computation of a one-step transition probability matrix is quite time-consuming in the traditional method. A short-time Gaussian approximation, which was proposed by Riske in the 1980s [17], has been reported in many published references [18]. It suggests that the response PDF of a nonlinear system subject to Gaussian noise approximately follows a Gaussian distribution with the low-order moments determined by moment equations [19, 20] or statistical simulations [21] when the time step is sufficiently small. The error of the short-time Gaussian approximation is within an error of order . However, for non-Gaussian cases such as Lévy process with shock and jump in presence of sea wave and Hurricane, the outliers may lead to the deviation of estimated mean value and variance from the core of real PDF, resulting in a failure of the short-time Gaussian approximation in probability evolution. Recently, Xu et al. [22] derived theoretically an approximated short-time solution aimed at a class of smooth fractional FPK equation with one-dimension, which provides a theoretical feasibility of achieving the non-Gaussian short-time approximation. The method in the cited work requires the participation of the analytic expressions of stochastic differential equations; however, the efficiency could suffer when the specific expressions are complex.
In this paper, we present an efficient numerical estimation method for approximating the distribution characteristics of short-time responses of a nonlinear system subject to Lévy noise, which is then used to compute the one-step transition probability matrix of the PI formulation. In this manner, the PCA is utilized to identify a probability core, and the short-time PDF in the core can be analytically estimated using just a small number of samples. The overall process is thus rather efficient, compared with full MCS. When a transient-state large transition is dealt with, a strategy of evolving probabilistic vector [23] is also employed in this paper to enhance the performance in storage and efficiency of computation.
This paper is devoted to promoting the idea of short-time probability approximation to the case of non-Gaussian input, by which to quantitatively uncover how the Lévy disturbance with shock and jump takes effect on response transitions in a nonsmooth system with bilateral clearances and piecewise elastic restoring forces. The interest of the paper is to detect the occurrence of response transition and uncover its path of PDF in the complex nonsmooth system, taking advantage of the proposed method. Also, the idea of the proposed method is applicable to the investigation of short-time probability approximation for other types of non-Gaussian inputs.
This paper is organized as follows: In Section 2, the idea of the proposed short-time probability approximation for non-Gaussian inputs is introduced in the frame of PI. In Section 3, the stochastic response transitions of a nonsmooth system with additive Lévy noise are investigated based on the deterministic global structure that is treated as a basic frame of responses. The performance of the proposed method is verified by comparing it with MCS. In Section 4, the conclusions are summarized.
2. Methodology
2.1. Formulation of the Short-Time Path Integral
Recall that the path integral method is used for obtaining the probability evolution of the stochastic system from a given initial distribution on a -dimensional state space , following the diffusion Markov process:where denotes probability at th instant time step, which is a probability vector distributed on the space with discrete grids. is a one-step transition probability matrix, and its th element of can be determined by the following formula:where is th grid on , denotes the total number of samples whose initial points are uniformly selected within domain, and is the number of samples that fall in domain.
An intractable difficulty is how can estimate the condition probability of efficiently. Traditionally, it can be statistically obtained by direct MCS . Though this way is very attractive due to its versatility and simplicity, a vast number of samples, up to hundreds for each initial state and even more for higher dimensions, are required to achieve an acceptable accuracy level of , which will result in a huge time consumption in the computation of whole evolutions of . The idea of short-time probability distribution approximation is a great way to evaluate efficiently a one-step transition probability matrix by an analytic probability density function (PDF). In reference [17, 18], when the evolving time step is sufficiently small, the response PDF of a nonlinear system under excitations of Gaussian noise is considered approximately as Gaussian within an error of order O , and then, the distribution is specified by the short-time mean value vector and the short-time variance matrix determined by moment evolution equations at time , which can be denoted as follows:
The success of the previous short-time Gaussian approximation is based on the cumulant-neglect closure of low-order moments. Yet, it is far from obtaining the expression of moment evolution equations for a nonsmooth system. What is more, due to the Lévy jump, low-order moments directly evaluated by statistical simulations may fail to depict the profile and core of the response distribution, resulting in an unacceptable deviation from its true PDF. This deviation will be illustrated in Figure 4(a) hereafter.
2.2. Idea of Short-Time Probability Estimation Based on PCA
To address the aforementioned problem with the deviation in estimating probability distribution due to isolated outliers, the principal component analysis is employed here to find a core region covering the majority of sample realizations. To characterize the short-time PDF evolving from th grid, the data set whose row denotes samples with dimensions is initially centered as bywhere defines a matrix centering for raw samples evolving from th grid, is a identity matrix with dimensions, and represents a -dimensional column vector whose elements are all one and is its transpose. Then, the eigenvalues and eigenvectors of can be obtained by singular value decompositionwhere and are orthogonal matrices. Columns of are eigenvectors representing the principal components of the centered data, and the diagonal matrix consists of corresponding eigenvalues representing the magnitudes of principal components. Thus, a hyperellipsoid whose center is at the is formed in a -dimensional state space. The coverage of its core is determined by axial directions denoted by the eigenvectors and axial lengths denoted by the corresponding scaled eigenvalues ( is a chosen scale factor). Of course, the dimensions of the hyperellipsoid may be reduced to be less than if there are parts of eigenvalues close to zero. The number of samples covered by the hyperellipsoid’s core is predominant in statistics, and its distribution on the core can be approximated to be Gaussian. In this manner, therefore, the short-time probability distribution can be characterized semianalytically by equation (3), where the mean value vector and variance matrix are estimated by only a few samples that fall into the hyperellipsoid in the short-time interval, avoiding lots of whole long-time integrals for governing equations to reduce computation cost.
To further efficiently capture the transient-state responses and large probability transition, the strategy of evolving probabilistic vector is introduced into the above PI method based on the short-time approximation. In this manner, by a prescribed fiducial probability , the grids are classified into the active region covered by and the inactive region outside of . Usually, the number of grids that reside in the inactive region is considerable, yet their joint probabilities are tiny [23]. For this reason, the inactive region will dose not be involved in the computation and storage to reduce computational effort. The iterative probability evolution rule is as follows:where , denotes the index of grids occupied inactive region. The fiducial probability in each short-time step is expressed by
Here, in the evaluating process, the one-step transition probability matrix is no longer a fix-sized matrix corresponding to grid resolutions, but rather a dynamical probability vector whose size is far less than the whole, depending on the scope of the active region. The performance of the proposed short-time probability approximation method can thus be enhanced significantly to capture a large-range probability transition.
3. Stochastic Dynamics of a Nonsmooth System
3.1. Stochastic Modeling
In this paper, we focus on a nonsmooth dynamical system with base motion and displacement constraint, which is a simplified model derived from a vibration isolation system with an amplitude limiter or an energy harvesting device with a designed stopper. The model is for qualitative discussion on stochastic dynamics. As shown in Figure 1, an object with lumped mass is supported vertically by a linear spring and damping. Also, the bilateral springs with stiffness coefficients are installed symmetrically along its direction of motion to exert a displacement restriction by an additional restoring force with cubic nonlinearity when the motion of mass exceeds the amplitude of clearances. The bilateral clearances are denoted by and . Then, the restoring force is modeled by a piecewise formula:
With base motion , the clearances between object and bilateral springs both are time-varying; that is, and , where is a constant denoting the initial value of clearances. In some cases, such as disturbances from sea waves, hurricanes, and environmental noise might not be typically Gaussian. The intermittent shocks in the persistent disturbances could be considered as an additive Lévy process with jump. Then, the mathematical model of the dynamic system can be derived in the Cartesian coordinates:where and are the linearized damping and stiffness coefficients. Here, the base motion is treated as a periodic excitation with amplitude and frequency ; that is, . Meanwhile, defines a random process obeying Lévy probability distribution to simulate an uncertain disturbance with shock and jump, which satisfies the characteristic function :where is a stability index, indicating the thickness of the tail of the Lévy distribution. denotes the skewness parameter, indicating the direction of the random jump. and are the mean and the generalized diffusion coefficient of random walks. The is an operator of the sign function. Actually, the Lévy process can be viewed as a compound of Wiener motion and compensated Poisson process. In particular, the Lévy process degenerates to be Gaussian for the case of parameters and .
For the simplification and wide applicability of the following analysis, the nondimensional parameters are introduced as follows:where denotes the characteristic length used for dimensionless. Equation (9) becomeswhere denotes nondimensional linearized damping, be nondimensional stiffness coefficient, and is nonlinear stiffness item, and then,
Although the expression in equation (12) should be regarded primarily as a model one for illustrating the analysis of basic effects involved in the case of nonsmooth characteristics and impact excitation, it may have potential application to the dynamics of suppression or utilization of a vibration device equipped on the moored body under ocean waves.
3.2. Short-Time Probability Approximation for the Nonsmooth System
To validate the effectiveness of the proposed method for the nonsmooth system and determine the appropriate estimation parameters, namely the scale factor and the time step , the coverage of the semianalytical probability core is here examined by overlapping 5000 stochastic samples of short-time responses obtained by MCS as a benchmark. In the following, the basic properties of the system in equation (12) are set as , , , , , unless otherwise indicated. Figure 2 illustrates that when in equation (12) with parameters of noise and , 89.6%, 92.6%, 95.7%, 97.6% (orange points) of 5000 samples generated by MCS from the identical initial state fall inside the cores of 2D-ellipses determined by two principal components with and 2.0, respectively. Clearly, the smaller the scale factor, the more the number of samples that are excluded, resulting in a drastic truncation of the probability distribution tail. However, if the scale factor is set to a too large value, the isolated outliers it contained may distort the Gaussian distribution within the core. Thus, both of the two cases will bring in an unacceptable loss of accuracy. In the stochastic analysis of the nonsmooth system, is selected to surround as many dispersed samples as there are in the ellipse denoting Gaussian core.
Based on this, Figure 3 reveals an optimal short-time interval for mapping. It can be seen that the majority of response realizations congregate in the given Gaussian core when , in accordance with the expectancy of the Gaussian core, yet the escape of responses from the estimated ellipse has occurred to the samples when , indicating the onset of dynamic transition that leads to the failure of short-time approximation. In other words, being is appropriate and suitable for the system.
(a)
(b)
(c)
(d)
Take to be the number of samples used in our method to estimate a probability core, which is much less than that of full MCS. Figure 4 shows a comparison of the short-time PDFs estimated by the proposed method, the direct low-order moments method, and the full MCS using samples as a benchmark. As can be seen, the PDF described by the first two order moments (see Figure 4(a)) based on a very small amount of samples shows a considerable deviation especially at the peak of probability, while the proposed short-time approximation technique combined with the PCA strategy (see Figure 4(b)) has a great accuracy of estimation, as compared with full MCS (see Figure 4(c)).
(a)
(b)
(c)
3.3. Probability Transition Induced by Lévy Noise
Below, the dynamics of the nonsmooth system is investigated in an interesting state space , which is discretized into grids with resolution . The full MCS will inevitably encounter a huge computation cost for such a high resolution.
When noise is absent, that is, , the deterministic global structure of the nonsmooth system is first examined using the generalized cell mapping method [24]. The purpose of the investigation is to understand the underlying basic frame that the stochastic response will obey, which involves what types of responses (attractors) exist and how attractive they are (basins of attraction). As seen in Figure 5(a), for a given designed initial clearance , there are six attractors coexisting in the state space, including two period-1 attractors ( and ) labeled by symbols and , two period-2 attractors by and as well as two period-3 attractors by and . Corresponding orbits of steady-state responses after 400 evolution periods are obtained by numerical integral from initial states locating in the different basins of attraction (see Figures 5(b)–5(d)). The response corresponding to the period-1 attractor () does not reach the minimum amplitude of dynamical clearance , indicating that the system runs without collision when staying in the basin of the attractor, while other responses mean bilateral elastic collision happened. Also, it is observed in Figure 5(a) that one chaotic saddle (red points) and isolated saddle points (red “”) are embedded in a wide range of fractal boundaries separating the corresponding basins of attraction. Obviously, the basin boundary is wrinkled and interweaved, and every grid discretized on the boundary has the nature that the points inside will go into at least three different basins. In other words, the basin boundary is Wada, which is drawn by colorful points.
Taking the Lévy excitation into account, the stochastic probability evolutions are then simulated numerically using the proposed method in Section 2. Consider that the system is initially working in the state of collision-free. It is easy to image that when a small noise is present, the response PDF can still stay around the initial basin of attraction for a long time. For a relatively large uncertainty of impact and jump in Lévy noise, as shown in Figure 6 with the parameters and , the probability distribution around the initial state is explosive increasingly, and rapidly, most response realizations are pushed out of the basin of the period-1 attractor. This suggests that a significant qualitative change in dynamics occurs since the stochastic period-1 attractor under fiducial probability touches the saddle on its basin boundary (isolated red point in Figure 6), leading to the blue-sky disappearance of the stochastic responses around the period-1 attractor. Eventually, the PDF of response escapes and flows along the unstable manifold of the saddle to where another distant period-1 attractor resides, after about 400 periods of responses.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
The transient results obtained by the proposed method are compared to those by full MCS. Here, the absolute error (AE) and accuracy (Acc) are used to quantify the validity of the proposed method, that is,where and denote the grid probabilities estimated by the proposed method and full MCS with samples, respectively. With a given fiducial probability = 95% in EPV, the AE contour maps shown in Figure 6 demonstrate that the maximum AE throughout the whole transient process is under the magnitude of , and overall Acc reaches 99.96%, indicating the proposed method performs well.
In engineering, the exit time in evolution is regarded as a measurement for evaluating the stochastic stability of the system. In this paper, it is defined as a time period during which the response PDF continues to remain in an initial basin of attraction under a certain probability threshold , which can be written as follows:where denotes probabilities of grids that cover a basin of attraction. Taking the proposed method, the exit times of transient response PDF escaped from initial collision-free motion are investigated with different noise parameters and = 90%. Figure 7(a) reveals that the response transition from one state to another in response slows down as the parameter increases. When is greater than 1.10, the response transition apparently gets put off. On the contrary, a small might hasten the transition. From the Figure 7(b), the parameter has negligible influence on the response transition, though. These investigations evidence that the proposed method performs well in capturing the large-range response transitions of a nonsmooth nonlinear system subject to the Lévy input process.
(a)
(b)
4. Conclusions
This paper develops a short-time probability approximation strategy in the PI frame for efficiently estimating one-step transition probability for a stochastic system with input of Lévy noise. In this manner, the PCA method based on data is employed to establish the principal directions of distribution and their magnitudes so that forming a probability core of hyperellipsoid dominates the majority of the probability of samples, and the distribution in the core is considered as Gaussian. The proposed method presents two advantages: by PCA, just a small amount of response realizations are necessary, and by EPV, a prescribed fiducial probability significantly reduces the region to be examined in state space. Therefore, the method takes much less computing cost than full simulations when capturing a large probability transition to predict the fatigue life and reliability of the mechanical system. Meanwhile, the proposed method has the potential for further application to the efficient investigation of transient-state probability evolution for other cases of nonsmooth and non-Gaussian excitation.
By the proposed method, the large-range probability transition of a Lévy-induced nonsmooth system is efficiently captured. Throughout this process, using just 30 samples to estimate a probability core, as opposed to hundreds of samples in MCS, the computation cost is relatively cheap, and the overall accuracy can approach 99.96%. Meanwhile, the exit time of response PDF in transition from a basin of attraction is denoted as an evaluation of the stochastic stability of the system. It is discovered that the index of Lévy noise can influence the transition speed significantly. From the perspective of engineering, the noise-induced phenomena cause a large-range transition of responses between different attractors, which may be beneficial to the structure design of vibration utilization devices such as energy harvesting system.
Data Availability
No underlying data were collected or produced in this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The first author would like to cordially thank Professor Li Ming for guiding his academic study selflessly over the past thirteen years, and the first author benefits from the outstanding contributions of the professor in the field of nonlinear dynamics. The authors gratefully acknowledge the support of the National Natural Science Foundation of China through Grant nos. 12172279, 11972282, and 12002267.