# A Short-Time Non-Gaussian Probability Approximation Based Method for Response Transitions of a Lévy-Driven Nonsmooth System

Authors: Zigang Li; Wang Yan; Jiaqi Kang; Ming Li
Journal: Mathematical Problems in Engineering
 (2023)

---

## 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.


---

## Body

## 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 Oτ2. 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 Section2, 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 aN-dimensional state space RN, following the diffusion Markov process:(1)pn+1=Pn⋅pn,where pn denotes probability at n th instant time step, which is a probability vector distributed on the space CN with Nc discrete grids. Pn is a one-step transition probability matrix, and its i,j th element of Pij can be determined by the following formula:(2)Pijn=∫CiNpxn,tnxjn−1,tn−1dxn−1=∑k=1Sipxn,tnxjn−1,tn−1,where CiN is i th grid on CN, S denotes the total number of samples whose initial points are uniformly selected within CjN domain, and Si is the number of samples that fall in CiN domain.An intractable difficulty is how can estimate the condition probabilitypx,τ|xj,0 of Pij efficiently. Traditionally, it can be statistically obtained by direct MCS Si/S. 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 Pij, which will result in a huge time consumption in the computation of whole evolutions of Pn. 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 τ2, and then, the distribution px,τ|xj,0 is specified by the short-time mean value vector mτ and the short-time variance matrix Bτ determined by moment evolution equations at time τ, which can be denoted as follows:(3)px,τ|xj,0=12πN/2Bτ1/2×exp−12x−mτTBτ−1x−mτ.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 Figure4(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 PDFpx,τ|xj,0 evolving from j th grid, the data set DS×N whose row denotes S samples with N dimensions is initially centered as D¯S×N by(4)D¯S×N=H⋅DS×N,H=Is−1S1s1sT,where H defines a matrix centering for raw samples evolving from j th grid, Is is a identity matrix with s dimensions, and 1s represents a s-dimensional column vector whose elements are all one and 1sT is its transpose. Then, the eigenvalues and eigenvectors of D¯S×N can be obtained by singular value decomposition(5)D¯S×N=UΣ2VT,where U and V are orthogonal matrices. Columns of V are eigenvectors representing the principal components of the centered data, and the diagonal matrix Σ2 consists of corresponding eigenvalues representing the magnitudes of principal components. Thus, a hyperellipsoid whose center is at the 1/sDS×NT1s is formed in a N-dimensional state space. The coverage of its core is determined by axial directions denoted by the eigenvectors V 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 N 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 probabilityPf, the grids are classified into the active region covered by Pf and the inactive region outside of Pf. 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:(6)Pij⋅pjn=0,j=r,Pij⋅pjn=pin+1,j≠r,where i,j=1,2,…,Nc, r denotes the index of grids occupied inactive region. The fiducial probability Pf in each short-time step is expressed by(7)Pf=∑i∑j≠rPijpjn.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.


## 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 aN-dimensional state space RN, following the diffusion Markov process:(1)pn+1=Pn⋅pn,where pn denotes probability at n th instant time step, which is a probability vector distributed on the space CN with Nc discrete grids. Pn is a one-step transition probability matrix, and its i,j th element of Pij can be determined by the following formula:(2)Pijn=∫CiNpxn,tnxjn−1,tn−1dxn−1=∑k=1Sipxn,tnxjn−1,tn−1,where CiN is i th grid on CN, S denotes the total number of samples whose initial points are uniformly selected within CjN domain, and Si is the number of samples that fall in CiN domain.An intractable difficulty is how can estimate the condition probabilitypx,τ|xj,0 of Pij efficiently. Traditionally, it can be statistically obtained by direct MCS Si/S. 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 Pij, which will result in a huge time consumption in the computation of whole evolutions of Pn. 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 τ2, and then, the distribution px,τ|xj,0 is specified by the short-time mean value vector mτ and the short-time variance matrix Bτ determined by moment evolution equations at time τ, which can be denoted as follows:(3)px,τ|xj,0=12πN/2Bτ1/2×exp−12x−mτTBτ−1x−mτ.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 Figure4(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 PDFpx,τ|xj,0 evolving from j th grid, the data set DS×N whose row denotes S samples with N dimensions is initially centered as D¯S×N by(4)D¯S×N=H⋅DS×N,H=Is−1S1s1sT,where H defines a matrix centering for raw samples evolving from j th grid, Is is a identity matrix with s dimensions, and 1s represents a s-dimensional column vector whose elements are all one and 1sT is its transpose. Then, the eigenvalues and eigenvectors of D¯S×N can be obtained by singular value decomposition(5)D¯S×N=UΣ2VT,where U and V are orthogonal matrices. Columns of V are eigenvectors representing the principal components of the centered data, and the diagonal matrix Σ2 consists of corresponding eigenvalues representing the magnitudes of principal components. Thus, a hyperellipsoid whose center is at the 1/sDS×NT1s is formed in a N-dimensional state space. The coverage of its core is determined by axial directions denoted by the eigenvectors V 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 N 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 probabilityPf, the grids are classified into the active region covered by Pf and the inactive region outside of Pf. 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:(6)Pij⋅pjn=0,j=r,Pij⋅pjn=pin+1,j≠r,where i,j=1,2,…,Nc, r denotes the index of grids occupied inactive region. The fiducial probability Pf in each short-time step is expressed by(7)Pf=∑i∑j≠rPijpjn.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 motionxbt 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 m is supported vertically by a linear spring and damping. Also, the bilateral springs with stiffness coefficients kn 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 xt exceeds the amplitude of clearances. The bilateral clearances are denoted by Δl and Δu. Then, the restoring force is modeled by a piecewise formula:(8)FKN=knx−Δu3,x>Δu,0,−Δl≤x≤Δu,−knx−Δl3,x<−Δl.Figure 1
The nonsmooth dynamic system with bilateral constraints.With base motionxbt, the clearances between object and bilateral springs both are time-varying; that is, Δlt=δ−xbt and Δut=δ+xbt, 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:(9)mx¨+cx˙−x˙b+kx−xb+FKN=ηt,where c and k are the linearized damping and stiffness coefficients. Here, the base motion xb is treated as a periodic excitation with amplitude A and frequency f; that is, xb=A⋅cosf⋅t. Meanwhile, ηt defines a random process obeying Lévy probability distribution to simulate an uncertain disturbance with shock and jump, which satisfies the characteristic function k:(10)ϕk=expiμk−σαkα1−iβsgnktanπα2,α≠1,expiμk−σk1+iβsgnk2πlnk,α=1,where α∈0,2 is a stability index, indicating the thickness of the tail of the Lévy distribution. β∈−1,1 denotes the skewness parameter, indicating the direction of the random jump. μ and σ are the mean and the generalized diffusion coefficient of random walks. The sgn⋅ 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 α=2 and β=0.For the simplification and wide applicability of the following analysis, the nondimensional parameters are introduced as follows:(11)X=xl,A¯=Al,Xb=xbl=A¯cosΩ¯⋅τ,ωn=km,Ω¯=fωn,δ¯=δl,τ=ωnt,Δ¯u=δ¯+Xb,Δ¯l=δ¯−Xb,where l denotes the characteristic length used for dimensionless. Equation (9) becomes(12)X¨+C¯X˙+K¯X+F¯KN=K¯Xb+C¯X˙b+ητ,where C¯=c/mωn denotes nondimensional linearized damping, K¯=k/mωn2 be nondimensional stiffness coefficient, and K¯N=kn/mωn2 is nonlinear stiffness item, and then,(13)F¯KN=K¯NX−Δ¯u3,X>Δ¯u,0,−Δl¯≤X≤Δ¯u,−K¯NX−Δl¯3,X<−Δl¯.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 C¯=0.01, K¯=0.6, K¯N=3000, A¯=0.20, Ω¯=1.0, unless otherwise indicated. Figure 2 illustrates that when τ=π/2 in equation (12) with parameters of noise σ=1.0 and α=1.5, 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 ε=1.0,1.2,1.5, 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, ε=1.5 is selected to surround as many dispersed samples as there are in the ellipse denoting Gaussian core.Figure 2
Short-time response samples obtained by MCS and cores of 2D-ellipsoid determined by PCA withε=1.0, ε=1.2, ε=1.5, and ε=2.0, respectively.Based on this, Figure3 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 τ≤π/2, in accordance with the expectancy of the Gaussian core, yet the escape of responses from the estimated ellipse has occurred to the samples when τ>π/2, indicating the onset of dynamic transition that leads to the failure of short-time approximation. In other words, τ being π/2 is appropriate and suitable for the system.Figure 3
Short-time response samples obtained by MCS and cores of 2D-ellipsoid determined by PCA with (a)τ=π/4, (b) τ=π/2, (c) τ=π, (d) τ=2π, respectively.
(a)(b)(c)(d)TakeS=30 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 106 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)).Figure 4
Comparison of short-time PDF estimated by (a) direct low-order moments by statistical simulations with 30 samples, (b) proposed probability core with 30 samples, (c) full Monte-Carlo simulations with106 samples as a benchmark.
(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−0.6≤X≤1.0,−1.0≤dX≤1.0, which is discretized into 5120×5120 grids with resolution 0.000312500×0.000390625. The full MCS will inevitably encounter a huge computation cost for such a high resolution.When noise is absent, that is,D=0, 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 δ¯=0.60, 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 Δ¯min=δ¯−A¯=0.40, 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.Figure 5
Global responses of the deterministic nonsmooth system atδ¯=0.60: (a) global dynamic structure; (b), (c), and (d) are orbits of steady-state responses for different attractors, where “+” and “⋆” stand for period-1 attractors, “♦” and “×” for period-2 attractors, “△” and “•” for period-3 attractors, and red “.” and red “⊙” for chaotic saddle and a regular saddle on Wada boundary (colorful points).Taking the Lévy excitation into account, the stochastic probability evolutions are then simulated numerically using the proposed method in Section2. 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 σ=1.0 and α=1.5, 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.Figure 6
Transient evolutions of response joint PDFs, where isolated red point is saddle point near period-1 attractor and pseudocolor maps indicate response PDFs with evolving time: (a) (c), (e) (g), (i) (k) results of the proposed method; (b) (d), (f) (h), (j) (l) AE contour map.
(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,(14)AE=Pg−P^g,Acc=1Nc∑i=1NcPgi−P^giP^gi,where Pg and P^g denote the grid probabilities estimated by the proposed method and full MCS with 106 samples, respectively. With a given fiducial probability Pf = 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 10−3, 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 thresholdPt, which can be written as follows:(15)Pt=∑i=1NBoAPBoAi,where PBoA denotes probabilities of NBoA 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 Pt = 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.Figure 7
Exit times of responses PDF for different parameters of Lévy noise with a probability thresholdPt = 90%: (a) parameter plane α-σ; (b) parameter plane α-β.
(a)(b)


## 3.1. Stochastic Modeling

In this paper, we focus on a nonsmooth dynamical system with base motionxbt 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 m is supported vertically by a linear spring and damping. Also, the bilateral springs with stiffness coefficients kn 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 xt exceeds the amplitude of clearances. The bilateral clearances are denoted by Δl and Δu. Then, the restoring force is modeled by a piecewise formula:(8)FKN=knx−Δu3,x>Δu,0,−Δl≤x≤Δu,−knx−Δl3,x<−Δl.Figure 1
The nonsmooth dynamic system with bilateral constraints.With base motionxbt, the clearances between object and bilateral springs both are time-varying; that is, Δlt=δ−xbt and Δut=δ+xbt, 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:(9)mx¨+cx˙−x˙b+kx−xb+FKN=ηt,where c and k are the linearized damping and stiffness coefficients. Here, the base motion xb is treated as a periodic excitation with amplitude A and frequency f; that is, xb=A⋅cosf⋅t. Meanwhile, ηt defines a random process obeying Lévy probability distribution to simulate an uncertain disturbance with shock and jump, which satisfies the characteristic function k:(10)ϕk=expiμk−σαkα1−iβsgnktanπα2,α≠1,expiμk−σk1+iβsgnk2πlnk,α=1,where α∈0,2 is a stability index, indicating the thickness of the tail of the Lévy distribution. β∈−1,1 denotes the skewness parameter, indicating the direction of the random jump. μ and σ are the mean and the generalized diffusion coefficient of random walks. The sgn⋅ 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 α=2 and β=0.For the simplification and wide applicability of the following analysis, the nondimensional parameters are introduced as follows:(11)X=xl,A¯=Al,Xb=xbl=A¯cosΩ¯⋅τ,ωn=km,Ω¯=fωn,δ¯=δl,τ=ωnt,Δ¯u=δ¯+Xb,Δ¯l=δ¯−Xb,where l denotes the characteristic length used for dimensionless. Equation (9) becomes(12)X¨+C¯X˙+K¯X+F¯KN=K¯Xb+C¯X˙b+ητ,where C¯=c/mωn denotes nondimensional linearized damping, K¯=k/mωn2 be nondimensional stiffness coefficient, and K¯N=kn/mωn2 is nonlinear stiffness item, and then,(13)F¯KN=K¯NX−Δ¯u3,X>Δ¯u,0,−Δl¯≤X≤Δ¯u,−K¯NX−Δl¯3,X<−Δl¯.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 C¯=0.01, K¯=0.6, K¯N=3000, A¯=0.20, Ω¯=1.0, unless otherwise indicated. Figure 2 illustrates that when τ=π/2 in equation (12) with parameters of noise σ=1.0 and α=1.5, 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 ε=1.0,1.2,1.5, 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, ε=1.5 is selected to surround as many dispersed samples as there are in the ellipse denoting Gaussian core.Figure 2
Short-time response samples obtained by MCS and cores of 2D-ellipsoid determined by PCA withε=1.0, ε=1.2, ε=1.5, and ε=2.0, respectively.Based on this, Figure3 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 τ≤π/2, in accordance with the expectancy of the Gaussian core, yet the escape of responses from the estimated ellipse has occurred to the samples when τ>π/2, indicating the onset of dynamic transition that leads to the failure of short-time approximation. In other words, τ being π/2 is appropriate and suitable for the system.Figure 3
Short-time response samples obtained by MCS and cores of 2D-ellipsoid determined by PCA with (a)τ=π/4, (b) τ=π/2, (c) τ=π, (d) τ=2π, respectively.
(a)(b)(c)(d)TakeS=30 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 106 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)).Figure 4
Comparison of short-time PDF estimated by (a) direct low-order moments by statistical simulations with 30 samples, (b) proposed probability core with 30 samples, (c) full Monte-Carlo simulations with106 samples as a benchmark.
(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−0.6≤X≤1.0,−1.0≤dX≤1.0, which is discretized into 5120×5120 grids with resolution 0.000312500×0.000390625. The full MCS will inevitably encounter a huge computation cost for such a high resolution.When noise is absent, that is,D=0, 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 δ¯=0.60, 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 Δ¯min=δ¯−A¯=0.40, 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.Figure 5
Global responses of the deterministic nonsmooth system atδ¯=0.60: (a) global dynamic structure; (b), (c), and (d) are orbits of steady-state responses for different attractors, where “+” and “⋆” stand for period-1 attractors, “♦” and “×” for period-2 attractors, “△” and “•” for period-3 attractors, and red “.” and red “⊙” for chaotic saddle and a regular saddle on Wada boundary (colorful points).Taking the Lévy excitation into account, the stochastic probability evolutions are then simulated numerically using the proposed method in Section2. 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 σ=1.0 and α=1.5, 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.Figure 6
Transient evolutions of response joint PDFs, where isolated red point is saddle point near period-1 attractor and pseudocolor maps indicate response PDFs with evolving time: (a) (c), (e) (g), (i) (k) results of the proposed method; (b) (d), (f) (h), (j) (l) AE contour map.
(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,(14)AE=Pg−P^g,Acc=1Nc∑i=1NcPgi−P^giP^gi,where Pg and P^g denote the grid probabilities estimated by the proposed method and full MCS with 106 samples, respectively. With a given fiducial probability Pf = 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 10−3, 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 thresholdPt, which can be written as follows:(15)Pt=∑i=1NBoAPBoAi,where PBoA denotes probabilities of NBoA 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 Pt = 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.Figure 7
Exit times of responses PDF for different parameters of Lévy noise with a probability thresholdPt = 90%: (a) parameter plane α-σ; (b) parameter plane α-β.
(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.



---

*Source: 4135812-2023-01-03.xml*