Comparison of stochastic stability boundaries for parametrically forced systems with application to ship rolling motion

Numerous accidents caused by parametric rolling have been reported on container ships and pure car carriers (PCCs). A number of theoretical studies have been performed to estimate the occurrence condition of parametric rolling in both regular and irregular seas. Some studies in random wave conditions have been the approximate extension of the occurrence conditions for regular waves (e.g. Maki et al). Furthermore, several researches have been based on the stochastic process in ocean engineering (Roberts and Dostal). This study tackled the parametric rolling in irregular seas from the stability of the system's origin. It provided a novel theoretical explanation of the instability mechanism for two cases: white noise parametric excitation and colored noise parametric excitation. The authors then confirmed the usefulness of the previously provided formulae by Roberts and Dostal through numerical examples.


Introduction
Parametric rolling, one of the threats to oceangoing vessels, is triggered by the change in the restoring moment caused by waves and from the viewpoint of nonlinear dynamical system theory, corresponds to parametric resonance.In the late 1990s, several accidents were reported on container vessels because of parametric rolling [4], and since the 2000s, there have been also reports of it occurring at PCTCs [5].
In our paper survey, we have found Watanabe's work [6] conducted in the 1930s.There is a long research history on parametric rolling, and theoretical studies have often been aimed at estimating the conditions of occurrence and motion amplitude of parametric rolling in regular waves.We present a review of representative studies to determine the conditions and amplitudes of parametric rolling in regular waves below.Kerwin [7] theoretically analyzed parametric rolling using the harmonic balance method with success, thus leading to subsequent works for predicting parametric rolling in regular seas by Zavodney et al. [8], Francescutto [9], Bulian [10], Spyrou [11], Umeda et al. [12], Maki et al. [1], and Sakai et al. [13].
Research on parametric rolling in irregular waves has been very active in recent years, with a history that can be traced back to the 1980s [1,14,15].Many studies have been conducted on this more generalized problem, taking a stochastic process approach in the fields of control engineering and mechanical engineering.For example, Samuels performed numerical experiments in which the damping term of the second-order differential equation was parametrically oscillated by white noise [16].Caughey discussed this topic [17,18].Khasminskii [19] discussed the stability of systems with parametric excitation terms, and Kozin [20] further developed Khasminskii's ideas, deriving stability conditions and comparing them with numerical calculations.For an overview of the progress made in these studies over time, references [21] and [22] can be consulted.
However, the problem we are addressing is not a system with a white noise parametric excitation term but one with a colored noise parametric excitation.For colored noise, obtaining results on the stability of the system's origin becomes difficult, as the solution method used in the white noise case cannot be directly applied.To address this, Infante [23] obtained the condition of stability using the Lyapunov stability theory, while Arnold [24] used perturbation expansions and calculated Lyapunov exponents to obtain analytical relations for stability.
For systems with parametric excitation terms of the colored noise type, the problem can sometimes be solved using stochastic averaging methods.Stratonovich [25] and Khasminskii [26] developed the Stratonovich-Khasminskii limit theorem, which Roberts applied to ship problems [2].More recently, Dostal proposed an energy-based stochastic averaging method [3], and Maruyama developed enhancements to this approach [27].Liu provides an extensive survey of the applicability of this method [28].Maruyama [29] also established an estimation method using the higher-order moment technique and successfully estimated the lateral angular acceleration using the results of this method [30].Most of these studies focused on obtaining the PDF of the response, but some, such as the studies by Ariaratnam and Tam [31] and [3], refer to the stability of the origin of the system.
This study aims to demonstrate the physics and mechanism of parametric rolling in irregular seas from the stability of the upright condition (the system's origin).In regular head sea cases, it is easy to understand that parametric rolling is attributed to the loss of asymptotic stability of the upright condition.For instance, as Maki et al. [1] demonstrated, its destabilization because of parametric excitation can be theoretically estimated by the averaging method.However, the occurrence of parametric rolling in irregular seas can be explained by the loss of stability of the system's origin.To explain this physics, the authors first present the system's stability with a stochastic coefficient modeled by the Wiener process.Secondly, they demonstrate the stability of the ship equation of motion with parametric excitation.In both considerations, several formulae to predict the stability threshold are examined.
Initial results from the investigation presented in this study were initially described by Maki et al. [32].In this paper, the results are presented more extensively, with more details, and with some revisions.

Notations
In this study, the n-dimensional Euclidean space is denoted by R n , and the set of real numbers for n = 1 is denoted by R. The expectation operation is denoted by E, and t represents time.The overdot of time-dependent variables indicates the derivative with respect to time t.

Equation of motion
In this study, the authors deal with the single-DoF (Degreeof-Freedom) roll equation of motion: In this equation, φ (t) is the roll angle, φ (t) is the roll angular velocity, and φ (t) is the roll angular acceleration.2ζ is the damping coefficient, c 1 is the restoring moment coefficient, f (t) is the parametric excitation term based on the restoring moment variation, and h(t) is the roll moment caused by waves, each a function of time with a stochastic variation.
In the colored noise case, the term f (t)φ is estimated using Grim's effective wave concept [33,34].This equation has linear damping and restoring components; thus, the nonlinear components are not essential for assessing the stability of the origin.

Parametric oscillation for white noise
In this section, the authors briefly review the existing results conducted for the system with white noise parametric excitation.This is useful for understanding the physics of random parametric excitations, and therefore they consider the following parametric excitation term: Here, Γ is a intensity strength of a white noise, and W (t) is a 1D standard Wiener process that satisfies the relation or and dW (t) is the increment of this Wiener process.Here, δ (t) means the Dirac's delta function.

Some remarks on the stability in stochastic systems
In this subsection, we analyze the stability of the system origin with a parametric excitation term represented by white noise.We briefly explain the mathematical concepts used in this analysis.
The problem addressed here is almost identical to that described in a recent paper by the authors, where the problem of ship maneuvering motion under white noise multiplicative noise from the aspect of stochastic disturbances is tackled [35].The governing equations, in that case, are almost identical to those in the present paper.However, attention must be paid to the presence of the Wong-Zakai's correction term [36], which appears when a real system is reduced to an Itô-type system of stochastic differential equations.If the Wong-Zakai's correction term is present, it can strongly affect the stability of the system's origin, as discussed in [21].However, as mentioned in our previous work [35], the Wong-Zakai's correction term does not exist in the system dealt with here.
In the following subsections, the authors will demonstrate the method to identify the system's stability driven by white noise.First, the authors must carefully examine the relationship between the real system and the corresponding SDE to analyze the random system.The present system is the same as the one dealt with in our previous paper [35]; hence, the Wong-Zakai's correction term [36] is zero; this topic is discussed in the paper of Kozin [21].
The roll angle and roll angular velocity are now represented as a state variable x(t) ∈ R 2 , hereafter.
Consider the following stochastic differential equation (SDE: Stochastic Differential Equation ).
where µ(x(t)) : R 2 → R 2 , and σ (x(t)) : R 2 → R 2 .Since there does not exist the Wong-Zakai's correction term [36], the drift and diffusion terms for the present system are given by µ(x(t),t) and σ (x(t),t), respectively Finally, the SDE that the authors tackled is as follows: We define the infinitesimal operator For the system that the authors tackle, the infinitesimal operator for a scalar function F(x(t)) ∈ R becomes: The final term in this equation, , is a characteristic term that appears in the framework of stochastic differential equations and can serve as a correction term for the dynamics of F(x) because of the addition of noise.This term can either destabilize or stabilize the system.
We consider the system which does not have the term of parametric excitation, that is Γ x 1 (t)dW (t): For the above system, the infinitesimal operator for a scalar function F(x) becomes: Therefore, the additive Wiener process does not stabilize or destabilize the present system without parametric excitation.This is an important point to note.

Stability of second moment
Here, the authors apply the following relation between the expectation and differentiation operators: Then, the following moment equation can be derived: This equation contains no information beyond the secondorder moments and has a closed structure; however, as Bogdanoff and Kozin [37] state, if a parametric excitation term of the colored noise is added, then the closed structure cannot be established anymore.For the above system, the characteristic equation becomes: Here, by applying the Routh-Hurwitz stability criteria, we can obtain an inequality that represents the stability boundary of the 2nd moment: 4.3 Infante's approach [23] As the authors explained in section 5.2, Infante's result applies to the white noise case.Therefore, for our system, that is Eq. ( 10), since [ f (t) 2 ] = σ 2 , the boundary of stability is: Kozin's [41] results apply to Infante's approach since f (t) is not limited to white noise, allowing us to extend the results to the colored noise problem.
4.4 Kozin's approach [41] Kozin [41] proposed a methodology to estimate a system's stability with multiplicative noise based on Khasminskii's work [19].This methodology allows us to calculate J 1 to study a system's stability with colored noise, not just with white noise, as in Infante's approach.
Here, C is a positive constant, and η ζ is defined as follows: The integrand has singularities at ±π/2, making its numerical computation challenging.Detailed descriptions can be found in [20] and our previous literature [35].Once J 1 has been successfully calculated, the system's origin is considered stable if the following equation is satisfied.
4.5 Arnold's approach for white-noise case [24] As shown in sec.5.3, Arnold obtained the stability boundary for the colored noise case, which the authors apply to the present white-noise problem.Now, they calculate eq. ( 25).
In the present case, C f defined in eq. ( 26) becomes: where δ (t) denotes the Dirac's delta function.Here, we define the spectral density of C f (t) as S f t f t .Then, Eq. ( 38) can be calculated.
4.6 Numerical results for white-noise system Next, the authors show comparisons of the theoretical formulae with numerical results.A Monte Carlo simulation is employed using Euler Maruyama's scheme [42].For each coefficient combination of ζ and Γ , 20 sample paths are generated, where the initial condition is x 1 (0) = 0.1 and x 2 (0) = 0.1.After a time of 160 sec, it is judged whether the path has converged to the system origin or not.Fig. 1 and Fig. 2 show the comparative results for c 1 = 1 and c 1 = 2, respectively.It can be seen that Kozin's method accurately predicts the boundary of the stability.Additionally, Arnold's method also predicts the boundary well in the vicinity of Γ 2 = 0, since it is based on the assumption of Γ 2 → 0. However, the methods based on the 2nd moment stability and Infante's method do not quantitatively correlate with MCS results, although the boundary trend can be qualitatively represented.

Parametric oscillation for colored noise
The stability of the following type of differential equation is analyzed in this section.
which includes the colored noise processes f (t) and h(t).

Modeling of parametric excitation term
The following expressions of the parametric excitation terms are adopted in section 5. Let f (t) be denoted as P(t).The parametric excitation terms adopted in this section are calculated by first obtaining the GM variation, which is calculated from the computed restoring arm GZ for each regular wave height with a ratio of the ship length and wavelength of 1 (λ /L = 1) based on the Froude-Krylov assumption [43].This GM variation is then combined with the time series data of Grim's effective waves to obtain the parametric excitation term P(t).
where ω 0 = √ c 1 or c 1 = ω 2 0 denotes the natural roll frequency and A w denotes the effective wave amplitude.The parametric excitation terms in the section are expressed using the polynomial approximation of the relation between the change in the restoring force ∆ GM and wave amplitude amidships.
Thereby, the processes P(t) and h(t) are stationary Gaussian processes with spectral density S f t f t , and S h t h t as: where It is worth noting that there exists a relation G(ω) = 2S(ω) between the two-sided spectrum S(ω) and the single-sided spectrum G(ω).
5.2 Infante's approach [23] In this subsection, the author slightly generalizes the method of Infante [23] to the present problem.It is now assumed that the system is driven by physical noise (colored noise).
Here, f (t) is a zero-mean, stationary, ergodic physical noise.We rewrote the equation as follows: A positive-defined matrix B is introduced as follows: According to Infante's paper [23], if the following relation is satisfied for some positive ε, then the main theory holds: let λ max [X] be the maximum eigenvalue of matrix X.
then the system is almost asymptotically stable, yielding: Here, we define the elements of the B matrix as follows: For the above, by using Schwarz's inequality, the following result can be finally obtained: 5.3 Arnold-Dostal's approach for colored-noise case [24,3] Arnold et al. [24] have obtained asymptotic results for the stability of the equation where f (t) denotes a stationary Gaussian process with spectral density S f t f t .The top Lyapunov exponent was determined for the system of y(t) = x 1 (t) exp(ζ t): where y(t) satisfies the following differential equation: The top Lyapunov exponent λ x was determined for the system of x(t), which can be represented as: A negative Lyapunov exponent yields the stability of the corresponding system; thus, if λ x < 0 in Eq. ( 34), the SDE in Eq. ( 34) is stable The condition of λ x = 0 for negligibly small P(t) becomes: This is an implicit function for ζ or c 1 because S f t f t implicitly includes H 1/3 .Therefore, iterative methods, such as Newton's method, should be applied to obtain the stability boundary.
5.4 Ariaratnam and Tam's approach [31] Ariaratnam and Tam [31] obtained the stability boundary with the use of the stochastic averaging method and utilized the outcome to Eq. ( 27).
Here, they analyze with the use of the stochastic averaging theorem proposed by Stratonovich [25] and Khasminskii [26], the derived averaged equation for roll amplitude A being as follows: where Ariaratnam and Tam [31] obtained the moment stability conditions for the equation they averaged.

First moment stability
5.4.2 Second moment stability

Condition on the PDF
By solving the Fokker-Planck equation, the stationary probability density function (PDF) was obtained as follows: For the above PDF, the condition of stability is considered to be ν > 0, yielding: The two results obtained here are identical to those obtained by Roberts [2] in later years.These results also demonstrate that the external moment moment term, h(t), does not affect the system's stability.Note that Eq. 43 exactly matches the Arnold-Dostal result (Eq.38) if ζ is sufficiently small.5.5 Results of the energy-based averaging method [3,27] Here, we present an approach using an energy-based stochastic average method.The target system is as follows: Here, the Hamiltonian H of the system is: Then, the following equation can be obtained: The 1D SDE with respect to H can be obtained as: In the above equation, drift terms, i.e. m 1 and m 2 , and diffusion term σ 2 can be represented by: where FPK equation can be obtained from Eq.(47, and then the PDF for H is obtained as: Here, C and C denote the normalization constant of the PDF.Using the transformation formula for the probability density, the probability density function for roll amplitude is: Here, C also denotes the normalization constant of the PDF.
From this, the asymptotic behavior of the probability density function can be categorized into three cases for the relation between k and the linear roll damping coefficient 2ζ .
Now calculate the equation of the boundary from (ii) in Eq.( 52).Let k denote the spectrum S P (ω) of P(t).Therefore, we have: By utilizing the above relation, 8ζ = k becomes as follows.
However, even if P(A) → +∞, it cannot be said that P(A) = 0 when A > 0, so it may be a non-conservative side estimation as shown in Fig. 4.This is because the condition in question represents the boundary between the behavior of two PDFs, Type A and Type B, as shown in Fig. 3.This is the boundary of the bifurcation phenomenon of the probability density function; see, for example, page 506 of Arnold [44].Thus, this condition does not directly indicate the stability of the system's origin.
Roberts [2] presented a conditional expression for the probability density function from the FPK equation with λ > 0, where λ is expressed as: Here, M means: Therefore, the boundary proposed by Roberts is: The above expression is obtained for the system with an external force.Furthermore, consider the following limit: By considering 1 + 2λ < 0, the following condition can be obtained: The boundary equation obtained from this method agrees with the boundary equation obtained from the energy-based stochastic averaging method.

Results and discussions
This subsection compares the analytical conditions presented so far with numerical calculations.Fig. 4 shows the final results obtained for the ITTC spectrum.The subject ship is C11, and the graph shows a comparison of the results of theoretical calculations and numerical simulations of the stochastic averaging method, Infante's approach, and Arnold's approach.

Infante's method
The results based on Infante's method show that the estimation is overly safe.This can be attributed to two factors.Firstly, Infante's method uses Eq. 29 to set the matrix B related to the Lyapunov function.Even though the Lyapunov function is used to obtain the stability condition, it is a sufficient condition and may not be the optimal choice.This likely leads to an overly safe estimation.Secondly, Infante transforms the equation using Schwartz's inequality, another factor contributing to the overly safe estimation.

Arnold-Dostal's method
From the figure, the results based on Arnold-Dostal's method explain the numerical results relatively well.However, since this theory is based on the assumption that σ → 0, the estimation accuracy is expected to deteriorate as σ increases, as seen in the results for white noise.Nonetheless, for the faced problem of vessel motion, the assumptions of Arnold-Dostal's method are considered to be valid to a great extent.

Averaging method
The analytical condition based on the averaging method appears to be essentially an estimation on the non-conservative side.This is because, as mentioned above, even if P(A) → +0 and P(A) → +∞, P(A) = 0 is not necessarily guaranteed with A > 0. Therefore, it is important to note that this is a condition that indicates the behavior of the probability density function near the origin, not directly indicating the origin's stability.Consequently, this method has the potential to be an estimation on the non-conservative side.

Mitigation of parametric rolling due to rudder control
As Söder et al. [45] have conducted, there is the potential to reduce the risk of parametric rolling using rudder control.The rudder is an inherent piece of equipment for the ship, making it an attractive option as it does not require any new active actuators or additional anti-rolling tanks.Furthermore, the control aims to restore the asymptotic stability of the system's origin, which is the major difference from regular roll mitigation controllers by rudder or fin-stabilizers.Thus, once the stability of the system's origin is restored, it is expected that roll motion will be completely eliminated in random head seas.ẍ1 (t) + 2ζ ẋ1 (t) + (c 1 + f (t))x 1 (t) = f R δ R (60) Fig. 4: Comparison of the stability diagram between 'Arnold-Dostal's approach', 'Infante's approach', 'stochastic averaging method', and numerical simulation for C11 container ship parameters and forcing due to a ITTC spectrum with significant wave height H 1/3 and mean period T 1 = 10s.
Here, f R denotes the hydrodynamic derivative of the roll moment with respect to the rudder angle δ R .Suppose that δ R has a feedback form as follows: In this research, the delay of rudder action is ignored for the sake of brevity; thereby, the system becomes: ẍ1 (t) + 2ζ ẋ1 (t) + (c 1 + f (t))x 1 (t) = 0 where 2ζ If 2ζ and c 1 , i.e., k 1 and k 2 , are selected below the thresholds for stability, then parametric rolling can be completely prevented.The authors have demonstrated the control policy in the previous sections, and it can be said that control to increase the "apparent damping force" of 2ζ ≡ 2ζ − k 2 will lead to complete prevention of parametric rolling.

Conclusion
In this study, the stability of the system with stochastically varying parametric excitation terms is discussed.Estimation formulae to predict the occurrence of instability because of the inclusion of multiplicative noise are introduced.The equations presented here are primarily based on those proposed by previous researchers, with only minor modifications to make them applicable to the parametric rolling of ships.The results of Arnold, for example, show that the stability boundaries can be captured with a relatively high accuracy, which is promising for practical use in the near future.