properties [27]:

$\displaystyle\text{MMD}^{2}(p,q)$
$\displaystyle=\mathbb{E}_{\boldsymbol{x}_{a},\boldsymbol{x}_{a}^{\prime}\sim
p}[\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{a}^{\prime})]+\mathbb{E}_{\boldsymbol{x}_{b},\boldsymbol{x}_{b}^{\prime}\sim
q}[\kappa(\boldsymbol{x}_{b},\boldsymbol{x}_{b}^{\prime})]-2\mathbb{E}_{\boldsymbol{x}_{a}\sim
p,\boldsymbol{x}_{b}\sim q}[\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{b})].$

Now, we can derive each term of this equation using our distributions $p\equiv
p_{\theta}(\boldsymbol{x}|y)$ and $q\equiv p_{\textup{data}}(\boldsymbol{x})$
for the OvA case, and $p\equiv p_{\theta}(\boldsymbol{x}|y_{a}),q\equiv
p_{\theta}(\boldsymbol{x}|y_{b})$ for the OvO case. In both scenarios, we aim
at finding an expectation over the data variable $x$ using only the
respectively known and estimable terms $p_{\theta}(y|\boldsymbol{x})$ and
$p_{\theta}(y)$.

##### OvA scenario

For the first term, we use Bayes’ theorem twice to get an expectation over two
variables $\boldsymbol{x}_{a}$ and $\boldsymbol{x}_{b}$ drawn from the data
distribution.

$\begin{split}\mathbb{E}_{\boldsymbol{x}_{a},\boldsymbol{x}_{a}^{\prime}\sim
p}&=\mathbb{E}_{\boldsymbol{x}_{a},\boldsymbol{x}_{a}^{\prime}\sim
p_{\theta}(\boldsymbol{x}|y)}\left[\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{a}^{\prime})\right],\\\
&=\mathbb{E}_{\boldsymbol{x}_{a},\boldsymbol{x}_{a}^{\prime}\sim
p_{\textup{data}}(\boldsymbol{x})}\left[\frac{p_{\theta}(y|\boldsymbol{x}_{a})p_{\theta}(y|\boldsymbol{x}_{a}^{\prime})}{p_{\theta}(y)^{2}}\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{a}^{\prime})\right].\end{split}$

For the second term, we do not need to perform anything particular as we
directy get an expectation over the data variabes $\boldsymbol{x}_{a}$ and
$\boldsymbol{x}_{b}$.

$\mathbb{E}_{\boldsymbol{x}_{b},\boldsymbol{x}_{b}^{\prime}\sim
q}=\mathbb{E}_{\boldsymbol{x}_{b},\boldsymbol{x}_{b}^{\prime}\sim
p_{\textup{data}}(\boldsymbol{x})}\left[\kappa(\boldsymbol{x}_{b},\boldsymbol{x}_{b}^{\prime})\right].$

The last term only needs Bayes theorem once, for the distribution $q$ is
directly replaced by the data distribution
$p_{\textup{data}}(\boldsymbol{x})$:

$\displaystyle\mathbb{E}_{\boldsymbol{x}_{a}\sim p,\boldsymbol{x}_{b}\sim q}$
$\displaystyle=\mathbb{E}_{\boldsymbol{x}_{a}\sim
p_{\theta}(\boldsymbol{x}|y),\boldsymbol{x}_{b}\sim
p_{\textup{data}}(\boldsymbol{x})}\left[\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{b})\right],$
$\displaystyle=\mathbb{E}_{\boldsymbol{x}_{a},\boldsymbol{x}_{b}\sim
px}\left[\frac{p_{\theta}(y|\boldsymbol{x}_{a})}{p_{\theta}(y)}\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{b})\right].$

Note that for the last term, we could replace
$p_{\theta}(y|\boldsymbol{x}_{a})$ by $p_{\theta}(y|\boldsymbol{x}_{b})$; that
would not affect the result since $\boldsymbol{x}_{a}$ and
$\boldsymbol{x}_{b}$ are independently drawn from
$p_{\textup{data}}(\boldsymbol{x})$. We thus replace all terms, and do not
forget to put a square root on the entire sum since we have computed so far
the squared MMD:

$\displaystyle\mathcal{I}_{\text{MMD}}^{\text{ova}}(x;y)$
$\displaystyle=\mathbb{E}_{p_{\theta}(y)}\left[\text{MMD}(p_{\theta}(\boldsymbol{x}|y),p_{\textup{data}}(\boldsymbol{x}))\right],$
$\displaystyle=\mathbb{E}_{p_{\theta}(y)}\left[\left(\mathbb{E}_{\boldsymbol{x}_{a},\boldsymbol{x}_{a}^{\prime}\sim
p_{\textup{data}}(\boldsymbol{x})}\left[\frac{p_{\theta}(y|\boldsymbol{x}_{a})p_{\theta}(y|\boldsymbol{x}_{a}^{\prime})}{p_{\theta}(y)^{2}}\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{a}^{\prime})\right]\right.\right.$
$\displaystyle\left.\left.\qquad+\mathbb{E}_{\boldsymbol{x}_{b},\boldsymbol{x}_{b}^{\prime}\sim
p_{\textup{data}}(\boldsymbol{x})}\left[\kappa(\boldsymbol{x}_{b},\boldsymbol{x}_{b}^{\prime})\right]-2\mathbb{E}_{\boldsymbol{x}_{a},\boldsymbol{x}_{b}\sim
p_{\textup{data}}(\boldsymbol{x})}\left[\frac{p_{\theta}(y|\boldsymbol{x}_{a})}{p_{\theta}(y)}\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{b})\right]\right)^{\frac{1}{2}}\right].$

Since all variables $\boldsymbol{x}_{a}$, $\boldsymbol{x}_{a}^{\prime}$,
$\boldsymbol{x}_{b}$ and $\boldsymbol{x}_{b}^{\prime}$ are independently drawn
from the same distribution $p_{\textup{data}}(\boldsymbol{x})$, we can replace
all of them by the variables $\boldsymbol{x}$ and $\boldsymbol{x}^{\prime}$.
We then use the linearity of the expectation and factorise by the kernel
$\kappa(\boldsymbol{x},\boldsymbol{x}^{\prime})$:

$\mathcal{I}_{\text{MMD}}^{\text{ova}}(x;y)=\mathbb{E}_{p_{\theta}(y)}\left[\mathbb{E}_{\boldsymbol{x},\boldsymbol{x}^{\prime}\sim
p_{\textup{data}}(\boldsymbol{x})}\left[\kappa(\boldsymbol{x},\boldsymbol{x}^{\prime})\left(\frac{p_{\theta}(y|\boldsymbol{x})p_{\theta}(y|\boldsymbol{x}^{\prime})}{p_{\theta}(y)^{2}}+1-2\frac{p_{\theta}(y|\boldsymbol{x})}{p_{\theta}(y)}\right)\right]^{\frac{1}{2}}\right].$

##### OvO scenario

The two first terms of the OvO MMD are the same as the first term of the OvA
setting, with a simple subscript $a$ or $b$ at the appropriate place. Only the
negative term changes. We once again use Bayes’ theorem twice:

$\displaystyle\mathbb{E}_{\boldsymbol{x}_{a}\sim p,\boldsymbol{x}_{b}\sim
q}[\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{b})]$
$\displaystyle=\mathbb{E}_{\boldsymbol{x}_{a}\sim
p_{\theta}(\boldsymbol{x}|y_{a}),\boldsymbol{x}_{b}p_{\theta}(\boldsymbol{x}|y_{b})}\left[\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{b})\right],$
$\displaystyle=\mathbb{E}_{\boldsymbol{x}_{a},\boldsymbol{x}_{b}\sim
p_{\textup{data}}(\boldsymbol{x})}\left[\frac{p_{\theta}(y_{a}|\boldsymbol{x}_{a})}{p_{\theta}(y_{a})}\frac{p_{\theta}(y_{b}|\boldsymbol{x}_{b})}{p_{\theta}(y_{b})}\kappa(\boldsymbol{x}_{a},\boldsymbol{x}_{b})\right].$

The final sum is hence similar to the OvA:

$\displaystyle\mathcal{I}_{\text{MMD}}^{\text{ovo}}(x;y)$
$\displaystyle=\mathbb{E}_{p_{\theta}(y_{a}),p_{\theta}(y_{b})}\left[\text{MMD}(p_{\theta}(\boldsymbol{x}|y_{a}),p_{\theta}(\boldsymbol{x}|y_{b}))\right],$
$\displaystyle=\mathbb{E}_{p_{\theta}(y_{a}),p_{\theta}(y_{b})}\left[\mathbb{E}_{\boldsymbol{x},\boldsymbol{x}^{\prime}\sim
p_{\textup{data}}(\boldsymbol{x})}\left[\kappa(\boldsymbol{x},\boldsymbol{x}^{\prime})\left(\frac{p_{\theta}(y_{a}|\boldsymbol{x})p_{\theta}(y_{a}|\boldsymbol{x}^{\prime})}{p_{\theta}(y_{a})^{2}}+\frac{p_{\theta}(y_{b}|\boldsymbol{x})p_{\theta}(y_{b}|\boldsymbol{x}^{\prime})}{p_{\theta}(y_{b})^{2}}\right.\right.\right.$
$\displaystyle\quad\left.\left.\left.-2\frac{p_{\theta}(y_{a}|\boldsymbol{x})p_{\theta}(y_{b}|\boldsymbol{x}^{\prime})}{p_{\theta}(y_{a})p_{\theta}(y_{b})}\right)\right]^{\frac{1}{2}}\right].$

###  F.3 Wasserstein distance (Proof of Prop. 3.2)

To compute the Wasserstein distance between the distributions
$p_{\theta}(\boldsymbol{x}|y=k)$, we estimate it using approximate
distributions. We replace $p_{\theta}(\boldsymbol{x}|y=k)$ by a weighted sum
of Dirac measures on specific samples $\boldsymbol{x}_{i}$: $p_{N}^{k}$:

$p_{\theta}(\boldsymbol{x}|y=k)\approx\sum_{i=1}^{N}m_{i}^{k}\delta_{\boldsymbol{x}_{i}}=p_{N}^{k},$

where $\\{m_{i}^{k}\\}_{i=1}^{N}$ is the set of weights. We now show that
computing the Wasserstein distance between these approximates converges to the
correct distance. We first need to show that $p_{N}^{k}$ weakly converges to
$p$. To that end, we will use the Portmanteau theorem [5]. Let $f$ be any
bounded and continuous function. Computing the expectation of such through
$p_{\theta}$ is:

$\mathbb{E}_{\boldsymbol{x}\sim
p_{\theta}(\boldsymbol{x}|y=k)}[f(\boldsymbol{x})]=\int_{\mathcal{X}}f(\boldsymbol{x})p_{\theta}(\boldsymbol{x}|y=k)d\boldsymbol{x},$

which can be estimated using self-normalised importance sampling [52, Chapter
9]. The proposal distribution we take for sampling is
$p_{\textup{data}}(\boldsymbol{x})$. Although we cannot evaluate both
$p_{\theta}(\boldsymbol{x}|y)$ and $p_{\textup{data}}(\boldsymbol{x})$ up to a
constant, we can evaluate their ratio up to a constant which is sufficient:

$\displaystyle\mathbb{E}_{\boldsymbol{x}\sim
p_{\theta}(\boldsymbol{x}|y=k)}[f(\boldsymbol{x})]$
$\displaystyle=\int_{\mathcal{X}}f(\boldsymbol{x})\frac{p_{\theta}(\boldsymbol{x}|y=k)}{p_{\textup{data}}(\boldsymbol{x})}p_{\textup{data}}(\boldsymbol{x})d\boldsymbol{x},$
$\displaystyle=\int_{\mathcal{X}}f(\boldsymbol{x})\frac{p_{\theta}(y=k|\boldsymbol{x})}{p_{\theta}(y=k)}p_{\textup{data}}(\boldsymbol{x})d\boldsymbol{x},$
$\displaystyle\approx\sum_{i=1}^{N}f(\boldsymbol{x}_{i})\frac{p_{\theta}(y=k|\boldsymbol{x}=\boldsymbol{x}_{i})}{\sum_{j=1}^{N}p_{\theta}(y=k|\boldsymbol{x}=\boldsymbol{x}_{j})}.$

Now, by noticing in the last line that the importance weights are self
normalised and add up to 1, we can identify them as the point masses of our
previous Dirac approximations:

$m_{i}^{k}=\frac{p_{\theta}(y=k|\boldsymbol{x}=\boldsymbol{x}_{i})}{\sum_{j=1}^{N}p_{\theta}(y=k|\boldsymbol{x}=\boldsymbol{x}_{j})}.$

This allows to write that the Monte Carlo estimation through importance
sampling of the expectation w.r.t $p_{\theta}(\boldsymbol{x}|y=k)$ is directly
the expectation taken on the discrete approximation $p_{N}^{k}$. We can
conclude that their is a convergence between the two expectations owing to the
law of large numbers:

$\lim_{N\rightarrow+\infty}\mathbb{E}_{\boldsymbol{x}\sim
p_{N}^{k}}[f(\boldsymbol{x})]=\mathbb{E}_{\boldsymbol{x}\sim
p_{\theta}(\boldsymbol{x}|y=k)}[f(\boldsymbol{x})].$

Since $f$ is bounded and continuous, the portmanteau theorem [5] states that
$p_{N}^{k}$ weakly converges to $p_{\theta}(\boldsymbol{x}|y=k)$ when defining
the importance weights as the normalised predictions cluster-wise.

To conclude, when two series of measures $p_{N}$ and $q_{N}$ weakly converge
respectively to $p$ and $q$, so does their Wasserstein distance [65, Corollary
6.9], hence:

$\lim_{N\rightarrow+\infty}\mathcal{W}_{c}(p_{N}^{k_{1}},p_{N}^{k_{2}})=\mathcal{W}_{c}\left(p_{\theta}(\boldsymbol{x}|y=k_{1})\|p_{\theta}(\boldsymbol{x}|y=k_{2})\right).$
(27)

For the one-vs-all Wasserstein GEMINI, we simply need to set the second
distribution to the empirical data distribution: $m_{i}=1/N$.

##  Appendix G Other maximisations of the Wasserstein distances

The Wasserstein-1 metric can be considered as an IPM defined over a set of
1-Lipschitz functions. Indeed, such writing is the dual representation of the
Wasserstein-1 metric:

$W_{c}(p\|q)=\sup_{f,\|f\|_{L}\leq 1}\mathbb{E}_{x\sim
p}[f(x)]-\mathbb{E}_{z\sim q}[f(z)].$

Yet, evaluating a supremum as an objective to maximise is hardly compatible
with the usual backpropagation in neural networks. This definition was used in
attempts to stabilise GAN training [1] by using 1-Lipschitz neural networks
[25]. However, the Lipschitz continuity was achieved at the time by weight
clipping, whereas other methods such as spectral normalisation [48] now allow
arbitrarily large weights. The restriction of 1-Lipscthiz functions to
1-Lipschitz neural networks does not equal the true Wasserstein distance, and
the term "neural net distance" is sometimes preferred [2]. Still, estimating
the Wasserstein distance using a set of Lipschitz functions derived from
neural networks architectures brings more difficulties to actually leverage
the true distance according to the energy cost $c$ of Eq. 13.

Globally, we hardly experimented the generic IPM for GEMINIs. Our efforts for
defining a set of 1-Lipschitz critics, one per cluster for OvA or one per pair
of clusters for OvO, to perform the neural net distance [2] were not fruitful.
This is mainly because such an objective requires the definition of one neural
network for the posterior distribution $p_{\theta}(y|\boldsymbol{x})$ and $K$
(resp. $K(K-1)/2$) other 1-Lipschitz neural networks for the OvA (resp. OvO)
critics, i.e. a large number of parameters. Moreover, this brings the problem
of designing not only one but many neural networks while the design of one
accurate architecture is already a sufficient problem.

##  Appendix H Choosing a GEMINI

The complexity of GEMINI increases with the distances previously mentionned
depending on the number of clusters $K$ and the number of samples per batch
$N$. It ranges from $\mathcal{O}(NK)$ for the usual MI to
$\mathcal{O}(K^{2}N^{3}\log{N})$ for the OvO Wasserstein-GEMINI. As an
example, we show in Figure 12 the average time of GEMINI as the number of
tasked clusters increases for both 10 samples per batch (Figure 12a) and 500
samples (Figure 12b). The batches consists in randomly generated prediction
and distances or kernel between randomly generated data.

(a) 10 samples per batch

(b) 500 samples per batch

Figure 12: Average performance time (in seconds) of GEMINIs as the number of
tasked clusters grows for batches of size 10 and 500 samples.

The OvO Wasserstein is the most complex, and so its usage should remain for 10
clusters or less overall. The second most time-consuming loss is the OvA
Wasserstein, however its tendency in optimisation to only find 2 clusters
makes it a inappropriate. The main difference also to notice between the
following MMD is regarding their memory complexity. The OvA MMD requires only
$\mathcal{O}(KN^{2})$ while the OvO MMD requires $\mathcal{O}(K^{2}N^{2})$.
This memory complexity should be the major guide to choosing one MMD-GEMINI or
the other. Thus, the minimal time-consuming and resource-demanding GEMINI is
the OvA MMD if we consider GEMINIs that incorporates knowledge of data through
kernels and distances. Other versions involving $f$-divergences have in fact
the same complexity as MI in our implementations, apart from the OvO TV which
reaches $\mathcal{O}(K^{2}N)$ in our implementation.

##  Appendix I Further speed-ups for the Wasserstein-OvO

The complexity of the Wasserstein metric is $\mathcal{O}(N^{3}\log{N})$ for a
batch of size $N$. Consequently, the complexity of the OvO Wasserstein reaches
$\mathcal{O}(K^{2}N^{3}\log{N})$ for $K$ clusters. This implies that this
GEMINI is hardly usable for a high number of clusters. To tackle this
complexity, we propose instead to sample $M$ independent pairs among the
$K(K-1)/2$ pairs of clusters to compare. We evaluate the OvO Wasserstein on
these pairs and scale it to the same order as if performed on all pairs. This
is stochastically equivalent to maximising the metric for all pairs:

$\hat{\mathcal{I}}^{\text{ovo}}_{\mathcal{W}}(\boldsymbol{x};y)=\frac{2M}{K(K-1)}\times\sum_{k,k^{\prime}\in
I}p_{\theta}(y=k)p_{\theta}(y=k^{\prime})\times\mathcal{W}_{\delta}\left(\\{m_{i}^{k}\\}_{i=1}^{N},\\{m_{i}^{k^{\prime}}\\}_{i=1}^{N}\right),$
(28)

with

$I=\\{(k_{n},k^{\prime}_{n})\\}_{n=1}^{M}\quad;(k_{n},k^{\prime}_{n})\sim\mathcal{U}\left(\\{1,\cdots,K\\}^{2}\right),$

a set of uniformly drawn pairs of clusters.

This optimisation requires however longer training time as a tradeoff for a
controlled complexity of $\mathcal{O}(MN^{3})$. Note that the same
optimisation can be applied to the OvA Wasserstein-GEMINI.

##  Appendix J All pair shortest paths distance

Sometimes, using distances such as the $\ell_{2}$ may not capture well the
shape of manifolds. To do so, we derive a metric using the all pair shortest
paths. Simply put, this metric consists in considering the number of closest
neighbors that separates two data samples. To compute it, we first use a sub-
metric that we note $d$, say the $\ell_{2}$ norm. This allows us to compute
all distances $d_{ij}$ between every sample $i$ and $j$. From this matrix of
sub-distances, we can build a graph adjacency matrix $W$ following the rules:

$W_{ij}=\left\\{\begin{array}[]{cr}1&d_{ij}\leq\epsilon\\\
0&d_{ij}>\epsilon\end{array}\right.,$

where $\epsilon$ is a chosen threshold such that the graph has sparse edges.
Our typical choice for $\epsilon$ is the 5% quantile of all $d_{ij}$.

We chose the graph adjacency matrix to be undirected, owing to the symmetry of
$d_{ij}$ and unweighted. Indeed, solving the all-pairs shortest paths involves
the Floyd-Warshall algorithm [66, 59] which complexity $\mathcal{O}(n^{3})$ is
not affordable when the number of samples $n$ becomes large. An undirected and
unweighted graph leverages performing $n$ times the breadth-first-search
algorithm, yielding a total complexity of $\mathcal{O}(n^{2}+ne)$ where $e$ is
the number of edges. Consequently, setting a good threshold $\epsilon$
controls the complexity of the shortest paths to finds. Our final distance
between two nodes $i$ and $j$ is eventually:

$c_{ij}=\left\\{\begin{array}[]{c r}\text{Shortest-path}^{W}(i,j)&\text{if it
exists.}\\\ n&\text{otherwise}\end{array}\right..$ (29)

This metric $c$ can then be incorporated inside the Wasserstein-GEMINI.

##  Appendix K Packages for experiments

For the implementation details, we use several packages with a python 3.8
version.

  * •

We use PyTorch [54] for all deep learning models and automatic
differentiation, as well as NumPy [30] for arrays handling.

  * •

We use Python Optimal Transport’s function `emd2`[20] to compute the
Wasserstein distances between weighted sums of Diracs.

  * •

We used the implementation of SIMCLR from PyTorch Lightning [19].

  * •

Small datasets such as isotropic Gaussian Mixture of score computations are
performed using scikit-learn[55].


∎

11institutetext: Atsuo Maki 22institutetext: Yuuki Maruyama 33institutetext:
Osaka University, 2-1 Yamadaoka, Suita, Osaka, Japan
33email<EMAIL_ADDRESS>44institutetext: Leo Dostal
55institutetext: Yaliu Liu 66institutetext: Institute of Mechanics and Ocean
Engineering, Hamburg University of Technology, 21043 Hamburg, Germany
66email<EMAIL_ADDRESS>
# Comparison of stochastic stability boundaries for parametrically forced
systems with application to ship rolling motion

Atsuo Maki Yuuki Maruyama Yaliu Liu Leo Dostal

(Received: date / Accepted: date)

###### Abstract

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. maki2011parametric ). Furthermore, several researches have been
based on the stochastic process in ocean engineering
(RobertsRoberts1982parametric and Dostaldostal2012non ). 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.

###### Keywords:

Parametric rolling irregular seas Stochastic differential equation Lyapunov
exponent

##  1 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 france2003 , and since the 2000s, there have been also
reports of it occurring at PCTCs rosen2012experience .

In our paper survey, we have found Watanabe’s work Watanabe1934 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 kerwin1955
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. zavodney1989 , Francescutto
francescutto2001 , Bulian bulian2004approximate , Spyrouspyrou2005paramet ,
Umeda et al. umeda2004nonlinear , Maki et al. maki2011parametric , and Sakai
et al. sakai2018 .

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
maki2011parametric ; blocki1980 ; themelis2008 . 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 samuels1960stability . Caughey discussed this topic caughey1960comments
; samuels1960remarks . Khasminskii Khasminskii1967 discussed the stability of
systems with parametric excitation terms, and Kozin Kozin1971 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 kozin1969 and Kozin1986 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 Infante1968 obtained the
condition of stability using the Lyapunov stability theory, while Arnold
Arnold1986Asymptotic 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 Stratonovich1963 and Khasminskii Khasminskii1966 developed the
Stratonovich–Khasminskii limit theorem, which Roberts applied to ship problems
Roberts1982parametric . More recently, Dostal proposed an energy-based
stochastic averaging method dostal2012non , and Maruyama developed
enhancements to this approach maruyama2022improved . Liu provides an extensive
survey of the applicability of this method liu2022applicability . Maruyama
maruyama2021moment also established an estimation method using the higher-
order moment technique and successfully estimated the lateral angular
acceleration using the results of this method maruyama2021moment_accel . Most
of these studies focused on obtaining the PDF of the response, but some, such
as the studies by Ariaratnam and Tam Ariaratnam1979 and dostal2012non , 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. maki2011parametric 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. Maki2022parametric_JASNAOE . In this paper, the
results are presented more extensively, with more details, and with some
revisions.

##  2 Notations

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

##  3 Equation of motion

In this study, the authors deal with the single-DoF (Degree-of-Freedom) roll
equation of motion:

$\ddot{\phi}+2\zeta\,\dot{\phi}+c_{1}\phi+f(t)\,\phi=h(\,t\,).$ (1)

In this equation, $\phi(t)$ is the roll angle, $\dot{\phi}(t)$ is the roll
angular velocity, and $\ddot{\phi}(t)$ is the roll angular acceleration.
$2\zeta$ 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)\phi$ is estimated using Grim’s effective wave concept grim1961 ;
umeda2022paramet .

This equation has linear damping and restoring components; thus, the nonlinear
components are not essential for assessing the stability of the origin.

##  4 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:

$f(t)\phi(t)\mathrm{d}t=\Gamma\phi(t)\mathrm{d}W(t)$ (2)

Here, $\Gamma$ is a intensity strength of a white noise, and $W(t)$ is a 1D
standard Wiener process that satisfies the relation

$\left\\{\begin{aligned} &\mathbb{E}[W(s)-W(t)]=0\\\
&\mathbb{E}[(W(s)-W(t))^{2}]=|t-s|\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{and}\\\ &\left\\{\begin{aligned}
&\mathbb{E}[W(s)W(t)]=\min(s,t)\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{or}\\\
&\mathbb{E}[\mathrm{d}W(s)\mathrm{d}W(t)]=\delta(t-s)\end{aligned}\right.\end{aligned}\right.$
(3)

and $\mathrm{d}W(t)$ is the increment of this Wiener process. Here,
$\delta(t)$ means the Dirac’s delta function.

###  4.1 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
tackledMaki2022stabilization . 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 Wong-Zakai1965 ,
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 kozin1969 . However, as mentioned in our previous
workMaki2022stabilization , 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 paperMaki2022stabilization ; hence, the
Wong-Zakai’s correction termWong-Zakai1965 is zero; this topic is discussed in
the paper of Kozinkozin1969 .

The roll angle and roll angular velocity are now represented as a state
variable $x(t)\in\mathbb{R}^{2}$, hereafter.

$x(t)\equiv\begin{bmatrix}x_{1}(t)\\\ x_{2}(t)\\\
\end{bmatrix}=\begin{bmatrix}\phi(t)\\\ \dot{\phi}(t)\\\ \end{bmatrix}$ (4)

Consider the following stochastic differential equation (SDE: Stochastic
Differential Equation).

$\mathrm{d}x(t)=\mu(x(t))\mathrm{d}t+\sigma(x(t))\mathrm{d}W(t)$ (5)

where $\mu(x(t)):\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$, and
$\sigma(x(t)):\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$.

Since there does not exist the Wong-Zakai’s correction termWong-Zakai1965 ,
the drift and diffusion terms for the present system are given by
$\mu(x(t),t)$ and $\sigma(x(t),t)$, respectively

$\quad\left\\{\begin{aligned} \mu(x(t),t)&=\begin{bmatrix}x_{2}(t)\\\ -\kappa
x_{2}(t)-c_{1}x_{1}(t)\end{bmatrix}\\\ \sigma(x(t),t)&=\begin{bmatrix}0\\\
\Gamma x_{1}(t)\\\ \end{bmatrix},\end{aligned}\right.$ (6)

Finally, the SDE that the authors tackled is as follows:

$\left\\{\begin{aligned} \mathrm{d}x_{1}(t)&=x_{2}(t)\mathrm{d}t\\\
\mathrm{d}x_{2}(t)&=-(2\zeta x_{2}(t)+c_{1}x_{1}(t))\mathrm{d}t+\Gamma
x_{1}(t)\mathrm{d}W(t)\\\ \end{aligned}\right.$ (7)

We define the infinitesimal operator $\mathcal{L}\left[\cdot\right]$ to
analyze the dynamics of the functional $f(x(t))\in\mathbb{R}$ of $x(t)$ as:

$\left\\{\begin{aligned}
\mathcal{L}\left[\cdot\right]\equiv&\left[\frac{\partial(\cdot)}{\partial
x}\right]\mu(x(t),t)\\\
&+\frac{1}{2}\sigma^{T}(x(t),t)\left[\frac{\partial^{2}(\cdot)}{\partial
x^{2}}\right]\sigma(x(t),t).\\\ &\text{or}\\\
\mathcal{L}\left[\cdot\right]\equiv&\sum_{i=1}^{2}\mu_{i}(x(t),t)\frac{\partial}{\partial
x_{i}}\\\
&+\frac{1}{2}\sum_{i=1}^{2}\sum_{j=1}^{2}\left[\sigma(x(t),t)\sigma^{T}(x(t),t)\right]_{ij}\frac{\partial}{\partial
x_{i}\partial x_{j}}\end{aligned}\right.$ (8)

For the system that the authors tackle, the infinitesimal operator for a
scalar function $F(x(t))\in\mathbb{R}$ becomes:

$\displaystyle\mathcal{L}F(x(t))=$ $\displaystyle x_{2}(t)\frac{\partial
F(x(t))}{\partial x_{1}}$ (9) $\displaystyle-(2\zeta
x_{2}(t)+c_{1}x_{1}(t))\frac{\partial F(x(t))}{\partial x_{2}}$
$\displaystyle+\frac{\Gamma^{2}}{2}x_{1}^{2}(t)\frac{\partial^{2}F(x(t))}{\partial
x_{2}^{2}}$

The final term in this equation, $\Gamma^{2}/2\cdot
x_{1}^{2}(t)\cdot\partial^{2}F(x)/\partial x_{2}^{2}$, 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 $\Gamma x_{1}(t)\mathrm{d}W(t)$:

$\left\\{\begin{aligned} \mathrm{d}x_{1}(t)&=x_{2}(t)\mathrm{d}t\\\
\mathrm{d}x_{2}(t)&=-(2\zeta
x_{2}(t)+c_{1}x_{1}(t))\mathrm{d}t+\Gamma\mathrm{d}W(t)\\\
\end{aligned}\right.$ (10)

For the above system, the infinitesimal operator for a scalar function $F(x)$
becomes:

$\displaystyle\mathcal{L}F(x(t))=$ $\displaystyle x_{2}(t)\frac{\partial
F(x(t))}{\partial x_{1}}$ (11) $\displaystyle-(2\zeta
x_{2}(t)+c_{1}x_{1}(t))\frac{\partial F(x(t))}{\partial x_{2}}$

Therefore, the additive Wiener process does not stabilize or destabilize the
present system without parametric excitation. This is an important point to
note.

Concerning the system in Eq. (10), several researches have been conducted in
the 1960s-1970s bogdanoff1962 ; kozin1963 ; caughey1965 ; Infante1968 ;
Kozin1971 . Here, the authors briefly introduce these theories and numerical
results.

###  4.2 Stability of second momentbogdanoff1962 ; kozin1986_2

In this subsection, the authors show the results based on the moment method
bogdanoff1962 ; kozin1986_2 . They apply the infinitesimal operator from Eq.
(8) to $x_{1}^{2}(t)$, $x_{1}(t)x_{2}(t)$, and $x_{2}^{2}(t)$, yielding:

$\left\\{\begin{aligned} \mathcal{L}x_{1}^{2}(t)&=2x_{1}(t)x_{2}(t)\\\
\mathcal{L}x_{1}(t)x_{2}(t)&=x_{2}^{2}(t)-2\zeta
x_{1}(t)x_{2}(t)-c_{1}x_{1}^{2}(t)\\\ \mathcal{L}x_{2}^{2}(t)&=-4\zeta
x_{2}^{2}(t)-2c_{1}x_{1}(t)x_{2}(t)+\Gamma^{2}x_{1}(t)\end{aligned}\right.$
(12)

Here, the authors apply the following relation between the expectation and
differentiation operators:

$\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}[g(x(t))]=\mathbb{E}[\mathcal{L}g(x(t))].$
(13)

Then, the following moment equation can be derived:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}\mathbb{E}[x_{1}^{2}(t)]\\\
\mathbb{E}[x_{1}(t)x_{2}(t)]\\\ \mathbb{E}[x_{2}^{2}(t)]\end{pmatrix}=$ (14)
$\displaystyle\begin{pmatrix}0&2&0\\\ -c_{1}&-2\zeta&1\\\
\Gamma^{2}&-2c_{1}&-4\zeta\end{pmatrix}\begin{pmatrix}\mathbb{E}[x_{1}^{2}(t)]\\\
\mathbb{E}[x_{1}(t)x_{2}(t)]\\\ \mathbb{E}[x_{2}^{2}(t)]\end{pmatrix}$

This equation contains no information beyond the second-order moments and has
a closed structure; however, as Bogdanoff and Kozin bogdanoff1962 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:

$s^{3}+6\zeta s^{2}+(8\zeta^{2}+4c_{1})s+(8c_{1}\zeta-2\Gamma^{2})=0$ (15)

Here, by applying the Routh–Hurwitz stability criteria, we can obtain an
inequality that represents the stability boundary of the 2nd moment:

$\Gamma^{2}=4c_{1}\zeta$ (16)

###  4.3 Infante’s approachInfante1968

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
$\mathbb{[}f(t)^{2}]=\sigma^{2}$, the boundary of stability is:

$\Gamma^{2}=4c_{1}\zeta^{2}$ (17)

Kozin’s kozin1973 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 approachkozin1973

Kozin kozin1973 proposed a methodology to estimate a system’s stability with
multiplicative noise based on Khasminskii’s work Khasminskii1967 . 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.

$\displaystyle\mathcal{J}_{1}=C\int_{-\pi/2}^{\pi/2}\biggl{(}$
$\displaystyle\cos 2\varphi-\frac{4\zeta}{\Gamma^{2}}\tan^{2}\varphi$ (18)
$\displaystyle+(1-c_{1})\frac{2}{\Gamma^{2}}\tan\varphi\biggr{)}\eta_{\zeta}(\varphi)\mathrm{d}\varphi$

Here, $C$ is a positive constant, and $\eta_{\zeta}$ is defined as follows:

$\displaystyle\eta_{\zeta}(\varphi)=$
$\displaystyle\exp\left[-\frac{2}{3\Gamma^{2}}\tan\varphi\left(3c_{1}+3\zeta\tan\varphi+\tan^{2}\varphi\right)\right]$
(19)
$\displaystyle\cdot\int_{-\pi/2}^{\varphi}\exp\biggl{[}\frac{2}{3\Gamma^{2}}\tan\theta\bigl{(}3c_{1}+3\zeta\tan\theta$
$\displaystyle+\tan^{2}\theta\bigr{)}\biggr{]}\sec^{2}\theta\mathrm{d}\theta$

The integrand has singularities at $\pm\pi/2$, making its numerical
computation challenging. Detailed descriptions can be found in Kozin1971 and
our previous literature Maki2022stabilization . Once $\mathcal{J}_{1}$ has
been successfully calculated, the system’s origin is considered stable if the
following equation is satisfied.

$\mathcal{J}_{1}=0$ (20)

###  4.5 Arnold’s approach for white-noise caseArnold1986Asymptotic

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:

$\displaystyle C_{f}(t)$ $\displaystyle=\mathbb{E}[f(t)f(0)]=\delta(t),$ (21)

where $\delta(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.

$\Gamma^{2}=8(1-\zeta^{2})\zeta$ (22)

###  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
Klodenbook1992 . For each coefficient combination of $\zeta$ and $\Gamma$, 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
$\Gamma^{2}=0$, since it is based on the assumption of $\Gamma^{2}\rightarrow
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.

Figure 1: Comparison of the stability diagram between numerical and analytical
results, $c_{1}=1.0$. Figure 2: Comparison of the stability diagram between
numerical and analytical results, $c_{1}=2.0$.

##  5 Parametric oscillation for colored noise

The stability of the following type of differential equation is analyzed in
this section.

$\left\\{\begin{aligned} \dot{x}_{1}(t)&=x_{2}(t)\\\ \dot{x}_{2}(t)&=-2\zeta
x_{2}(t)-(c_{1}+f(t))x_{1}(t)+h(t),\end{aligned}\right.$ (23)

which includes the colored noise processes $f(t)$ and $h(t)$.

###  5.1 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 $\mathrm{GM}$
variation, which is calculated from the computed restoring arm $\mathrm{GZ}$
for each regular wave height with a ratio of the ship length and wavelength of
$1$ ($\lambda/L=1$) based on the Froude-Krylov assumption Hamamoto1991 . This
GM variation is then combined with the time series data of Grim’s effective
waves to obtain the parametric excitation term $P(t)$.

$\displaystyle f(t)=P(t)=c_{1}P^{\prime}(t).$ (24)

where $\omega_{0}=\sqrt{c_{1}}$ or $c_{1}=\omega_{0}^{2}$ denotes the natural
roll frequency and $A_{\text{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
$\Delta\mathrm{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:

$\left\\{\begin{aligned}
S_{f_{t}f_{t}}&=\frac{1}{2\pi}\int_{-\infty}^{\infty}C_{f}(t)e^{-i\omega
t}\mathrm{d}t\\\
S_{h_{t}h_{t}}&=\frac{1}{2\pi}\int_{-\infty}^{\infty}C_{h}(t)e^{-i\omega
t}\mathrm{d}t\end{aligned}\right.$ (25)

where

$\left\\{\begin{aligned} C_{f}(t)&=\mathbb{E}[f(t)f(0)]\\\
C_{h}(t)&=\mathbb{E}[h(t)h(0)]\end{aligned}\right.$ (26)

It is worth noting that there exists a relation $G(\omega)=2S(\omega)$ between
the two-sided spectrum $S(\omega)$ and the single-sided spectrum $G(\omega)$.

###  5.2 Infante’s approachInfante1968

In this subsection, the author slightly generalizes the method of Infante
Infante1968 to the present problem. It is now assumed that the system is
driven by physical noise (colored noise).

$\ddot{x}_{1}(t)+2\zeta\dot{x}_{1}(t)+(c_{1}+f(t))x_{1}(t)=0$ (27)

Here, $f(t)$ is a zero-mean, stationary, ergodic physical noise. We rewrote
the equation as follows:

$\begin{gathered}\dot{x}(t)=[A+F(t)]x\\\ \left\\{\begin{aligned}
A=\begin{bmatrix}0&1\\\ -c_{1}&-2\zeta\\\
\end{bmatrix},F=\begin{bmatrix}0&0\\\ -f(t)&0\\\
\end{bmatrix}\end{aligned}\right.\end{gathered}$ (28)

A positive-defined matrix $B$ is introduced as follows:

$B=\begin{bmatrix}\alpha_{1}^{2}+\alpha_{2}&\alpha_{1}\\\
\alpha_{1}&1\end{bmatrix}$ (29)

According to Infante’s paper Infante1968 , if the following relation is
satisfied for some positive $\epsilon$, then the main theory holds: let
$\lambda_{\mathrm{max}}[X]$ be the maximum eigenvalue of matrix $X$.

$\mathbb{E}[\lambda_{\mathrm{max}}[A^{\top}+F^{\top}(t)+B(A+F(t))B^{-1}]]<-\epsilon,$
(30)

then the system is almost asymptotically stable, yielding:

$\displaystyle\mathbb{E}[\lambda_{\mathrm{max}}[A^{\top}+F^{\top}(t)+B(A+F(t))B^{-1}]]$
(31) $\displaystyle=-2\zeta+$
$\displaystyle\sqrt{4(\zeta-\alpha_{1})^{2}+\frac{[\alpha_{2}+\alpha_{1}^{2}-c_{1}-f(t)+2\alpha_{1}(\zeta-\alpha_{1})]^{2}}{\alpha_{2}}}$

Here, we define the elements of the $B$ matrix as follows:

$\left\\{\begin{aligned} \alpha_{1}&=\zeta\\\
\alpha_{2}&=\zeta^{2}+c_{1}\end{aligned}\right.,$ (32)

For the above, by using Schwarz’s inequality, the following result can be
finally obtained:

$\mathbb{E}[(f(t))^{2}]<4c_{1}\zeta^{2}$ (33)

###  5.3 Arnold-Dostal’s approach for colored-noise caseArnold1986Asymptotic ;
dostal2012non

Arnold et al. Arnold1986Asymptotic have obtained asymptotic results for the
stability of the equation

$\displaystyle\ddot{x}_{1}(t)+2\zeta\dot{x}_{1}(t)+(c_{1}+f(t))x_{1}(t)=0,$
(34)

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(\zeta t)$:

$\lambda_{y}=\lim_{t\rightarrow\infty}\log(|y(t)|^{2}+|\dot{y}(t)|^{2})^{\frac{1}{2}},$
(35)

where $y(t)$ satisfies the following differential equation:

$\displaystyle\ddot{y}(t)+[(c_{1}-\zeta^{2})-f(t)]y(t)=0,$ (36)

The top Lyapunov exponent $\lambda_{x}$ was determined for the system of
$x(t)$, which can be represented as:

$\displaystyle\lambda_{x}\approx-\zeta+\frac{\pi}{4(c_{1}-\zeta^{2})}S_{f_{t}f_{t}}\left(2\sqrt{c_{1}-\zeta^{2}}\right)$
(37)

A negative Lyapunov exponent yields the stability of the corresponding system;
thus, if $\lambda_{x}<0$ in Eq. (34), the SDE in Eq. (34) is stable The
condition of $\lambda_{x}=0$ for negligibly small $P(t)$ becomes:

$\displaystyle-\zeta+\frac{\pi}{4(c_{1}-\zeta^{2})}S_{f_{t}f_{t}}\left(2\sqrt{c_{1}-\zeta^{2}}\right)=0.$
(38)

This is an implicit function for $\zeta$ 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 Ariaratnam1979

Ariaratnam and Tam Ariaratnam1979 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 StratonovichStratonovich1963 and KhasminskiiKhasminskii1966 , the derived
averaged equation for roll amplitude $A$ being as follows:

$\begin{gathered}\mathrm{d}A=\left(-\alpha
A+\frac{\beta}{2A}\right)\mathrm{d}t+\left(\gamma
A^{2}+\beta\right)\mathrm{d}W(t)\\\ \text{where}\,\,\,\left\\{\begin{aligned}
\alpha&=\zeta-\frac{3\pi}{8c_{1}}S_{f_{t}f_{t}}(\omega)\\\
\beta&=\frac{\pi}{c_{1}^{2}}S_{h_{t}h_{t}}(\omega)\\\
\gamma&=\frac{1}{4c_{1}}S_{f_{t}f_{t}}(\omega)\end{aligned}\right.\end{gathered}$
(39)

Ariaratnam and TamAriaratnam1979 obtained the moment stability conditions for
the equation they averaged.

####  5.4.1 First moment stability

$\zeta>\frac{3\pi}{8c_{1}}S_{f_{t}f_{t}}(2\omega_{0})$ (40)

####  5.4.2 Second moment stability

$\zeta>\frac{\pi}{2c_{1}}S_{f_{t}f_{t}}(2\omega_{0})$ (41)

####  5.4.3 Condition on the PDF

By solving the Fokker-Planck equation, the stationary probability density
function (PDF) was obtained as follows:

$\begin{gathered}\mathcal{P}(A)=2\gamma\nu\beta^{\nu}\frac{A}{(\gamma
A^{2}+\beta)^{\nu+1}}\\\
\text{where}\,\,\,\nu=\frac{1}{2}+\frac{\alpha}{\gamma}\end{gathered}$ (42)

For the above PDF, the condition of stability is considered to be $\nu>0$,
yielding:

$\zeta>\frac{\pi}{4c_{1}}S_{f_{t}f_{t}}(2\omega_{0})$ (43)

The two results obtained here are identical to those obtained by Roberts
Roberts1982parametric 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 $\zeta$
is sufficiently small.

###  5.5 Results of the energy-based averaging method dostal2012non ;
maruyama2022improved

Here, we present an approach using an energy-based stochastic average method.
The target system is as follows:

$\ddot{x}+2\zeta\dot{x}+(c_{1}+P(t))x=0$ (44)

Here, the Hamiltonian $\mathcal{H}$ of the system is:

$\mathcal{H}(x,\dot{x})=\dfrac{\dot{x}^{2}}{2}+\dfrac{c_{1}}{2}x^{2}.$ (45)

Then, the following equation can be obtained:

$\left\\{\begin{array}[]{l}\dfrac{\mathrm{d}}{\mathrm{d}\,t}x=\dot{x}\\\ \\\
\dfrac{\mathrm{d}}{\mathrm{d}\,t}\mathcal{H}=-2\varepsilon\zeta\dot{x}^{2}-\sqrt{\varepsilon}P(t)x\dot{x}\end{array}\right.$
(46)

The 1D SDE with respect to $\mathcal{H}$ can be obtained as:

$\mathrm{d}\mathcal{H}=(m_{1}(\mathcal{H})+m_{2}(\mathcal{H}))\mathrm{d}t+\sigma(\mathcal{H})\mathrm{d}W$
(47)

In the above equation, drift terms, i.e. $m_{1}$ and $m_{2}$, and diffusion
term $\sigma^{2}$ can be represented by:

$\left\\{\begin{aligned} m_{1}(\mathcal{H})&=-2\zeta\mathcal{H}\\\
m_{2}(\mathcal{H})&=k\mathcal{H}\\\
\sigma^{2}(\mathcal{H})&=k\mathcal{H}^{2}\end{aligned}\right.$ (48)

where

$k=\dfrac{1}{c_{1}}\int_{0}^{\infty}R\left(\tau\right)\cos
2\sqrt{c_{1}}\tau\mathrm{d}\tau$ (49)

FPK equation can be obtained from Eq.(47, and then the PDF for $\mathcal{H}$
is obtained as:

$\displaystyle\mathcal{P}(\mathcal{H})$
$\displaystyle=\frac{C}{k\mathcal{H}^{2}}\exp\left(2\int_{\mathcal{H}_{\text{imt
}}}^{\mathcal{H}}\dfrac{-2\zeta\theta+k\theta}{k\theta^{2}}\mathrm{d}\theta\right)$
(50) $\displaystyle=C^{\prime}\mathcal{H}^{-\frac{4\zeta}{k}}$

Here, $C$ and $C^{\prime}$ denote the normalization constant of the PDF. Using
the transformation formula for the probability density, the probability
density function for roll amplitude is:

$\mathcal{P}(A)=\mathcal{P}(\mathcal{H})\left|\dfrac{\mathrm{d}\mathcal{H}}{\mathrm{d}A}\right|=C^{\prime\prime}A^{\left(1-\frac{8\zeta}{k}\right)}$
(51)

Here, $C^{\prime\prime}$ 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\zeta$.

$\begin{array}[]{l}\displaystyle(i)\quad\,\,8\zeta>k>0\quad\,\lim_{A\rightarrow+0}\mathcal{P}(A)\rightarrow+\infty\\\
\\\
\displaystyle(ii)\quad\,8\zeta=k\quad\quad\,\,\,\,\,\lim_{A\rightarrow+0}\mathcal{P}(A)=\mathrm{Const.}\\\
\\\ \displaystyle(iii)\quad
k>8\zeta>0\quad\lim_{A\rightarrow+0}\mathcal{P}(A)\rightarrow 0\end{array}$
(52)

Now calculate the equation of the boundary from $(ii)$ in Eq.(52). Let $k$
denote the spectrum $S_{P}(\omega)$ of $P(t)$. Therefore, we have:

$k=\dfrac{\pi}{2c_{1}}S_{P}(2\sqrt{c_{1}})$ (53)

By utilizing the above relation, $8\zeta=k$ becomes as follows.

$S_{P}(2\sqrt{c_{1}})=\dfrac{16}{\pi}c_{1}\zeta.$ (54)

However, even if $\mathcal{P}(A)\rightarrow+\infty$, it cannot be said that
$\mathcal{P}(A)\neq 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
ArnoldArnold1995random . Thus, this condition does not directly indicate the
stability of the system’s origin.

Figure 3: Schematic view of two PDFs

Roberts Roberts1982parametric presented a conditional expression for the
probability density function from the FPK equation with $\lambda>0$, where
$\lambda$ is expressed as:

$\lambda=\dfrac{\zeta-M}{M}$ (55)

Here, $M$ means:

$M=\dfrac{\pi}{8c_{1}}S_{P}(2\sqrt{c_{1}})$ (56)

Therefore, the boundary proposed by Roberts is:

$S_{P}(2\sqrt{c_{1}})=\dfrac{8}{\pi}c_{1}\zeta.$ (57)

The above expression is obtained for the system with an external force.
Furthermore, consider the following limit:

$\lim_{A\rightarrow+0}\mathcal{P}(A)\rightarrow+\infty$ (58)

By considering $1+2\lambda<0$, the following condition can be obtained:

$2\zeta<M$ (59)

The boundary equation obtained from this method agrees with the boundary
equation obtained from the energy-based stochastic averaging method.

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

Figure 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}=10$s.

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

####  5.6.2 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 $\sigma\rightarrow 0$, the estimation accuracy is expected to
deteriorate as $\sigma$ 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.

####  5.6.3 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 $\mathcal{P}(A)\rightarrow+0$ and
$\mathcal{P}(A)\rightarrow+\infty$, $\mathcal{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.

##  6 Mitigation of parametric rolling due to rudder control

As Söder et al. soder2013parametric 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.

$\ddot{x}_{1}(t)+2\zeta\dot{x}_{1}(t)+(c_{1}+f(t))x_{1}(t)=f_{\mathrm{R}}\delta_{\mathrm{R}}$
(60)

Here, $f_{R}$ denotes the hydrodynamic derivative of the roll moment with
respect to the rudder angle $\delta_{\mathrm{R}}$. Suppose that
$\delta_{\mathrm{R}}$ has a feedback form as follows:

$\delta_{\mathrm{R}}(t)=k_{1}x_{1}(t)+k_{2}x_{2}(t)$ (61)

In this research, the delay of rudder action is ignored for the sake of
brevity; thereby, the system becomes:

$\begin{gathered}\ddot{x}_{1}(t)+2\zeta^{\prime}\dot{x}_{1}(t)+(c_{1}^{\prime}+f(t))x_{1}(t)=0\\\
\text{where}\,\left\\{\begin{aligned} 2\zeta^{\prime}&\equiv 2\zeta-k_{2}\\\
c_{1}^{\prime}&\equiv c_{1}-k_{1}\end{aligned}\right.\end{gathered}$ (62)

If $2\zeta^{\prime}$ and $c_{1}^{\prime}$, 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\zeta^{\prime}\equiv 2\zeta-k_{2}$ will lead to complete
prevention of parametric rolling.

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

###### Acknowledgements.

This study was supported by a Grant-in-Aid for Scientific Research from the
Japan Society for the Promotion of Science (JSPS KAKENHI Grant #22H01701).
Further, this work was partly supported by the JASNAOE collaborative research
program / financial support. The authors are also thankful to Enago
(www.enago.jp) for reviewing the English language.

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * (1) A. Maki, N. Umeda, S. Shiotani, E. Kobayashi, Parametric rolling prediction in irregular seas using combination of deterministic ship dynamics and probabilistic wave theory, Journal of marine science and technology 16(3), 294 (2011)
  * (2) J. Roberts, Effect of parametric excitation on ship rolling motion in random waves, Journal of Ship Research 26(04), 246 (1982)
  * (3) L. Dostal, E. Kreuzer, N. Sri Namachchivaya, Non-standard stochastic averaging of large-amplitude ship rolling in random seas, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 468(2148), 4146 (2012)
  * (4) W.N. France, M. Levadou, T.W. Treakle, J.R. Paulling, R.K. Michel, C. Moore, An investigation of head-sea parametric rolling and its influence on container lashing systems, Marine Technology and SNAME news 40(01), 1 (2003)
  * (5) A. Rosén, M. Huss, M. Palmquist, Experience from parametric rolling of ships, Parametric resonance in dynamical systems pp. 147–165 (2012)
  * (6) Y. Watanabe, On the dynamical properties of the transverse instability of a ship due to pitching, Journal of Zosen Kiokai 1934(53), 51 (1934)
  * (7) J.E. Kerwin, Notes on rolling in longitudinal waves, International Shipbuilding Progress 2(16), 597 (1955)
  * (8) L.D. Zavodney, A. Nayfeh, N. Sanchez, The response of a single-degree-of-freedom system with quadratic and cubic non-linearities to a principal parametric resonance, Journal of Sound and Vibration 129(3), 417 (1989)
  * (9) A. Francescutto, An experimental investigation of parametric rolling in head waves, J. Offshore Mech. Arct. Eng. 123(2), 65 (2001)
  * (10) G. Bulian, Approximate analytical response curve for a parametrically excited highly nonlinear 1-dof system with an application to ship roll motion prediction, Nonlinear analysis: real world applications 5(4), 725 (2004)
  * (11) K. Spyrou, Design criteria for parametric rolling, Oceanic Engineering International 9(1), 11 (2005)
  * (12) N. Umeda, H. Hashimoto, D. Vassalos, S. Urano, K. Okou, Nonlinear dynamics on parametric roll resonance with realistic numerical modelling, International shipbuilding progress 51(2-3), 205 (2004)
  * (13) M. Sakai, N. Umeda, T. Yano, A. Maki, N. Yamashita, A. Matsuda, D. Terada, Averaging methods for estimating parametric roll in longitudinal and oblique waves, Journal of Marine Science and Technology 23(3), 413 (2018)
  * (14) W. Blocki, Ship safety in connection with parametric resonance of the roll, International Shipbuilding Progress 27(306), 36 (1980)
  * (15) N. Themelis, K.J. Spyrou, Probabilistic assessment of ship stability based on the concept of critical wave groups.  in _Proc. 10th Intl. Ship Stability Workshop_ (2008), pp. 18.1–18.19
  * (16) J.C. Samuels, On the stability of random systems and the stabilization of deterministic systems with random noise, The Journal of the Acoustical Society of America 32(5), 594 (1960)
  * (17) T. Caughey, Comments on “on the stability of random systems”, The Journal of the Acoustical Society of America 32(10), 1356 (1960)
  * (18) J.C. Samuels, Remarks on dr. caughey’s comments on my paper entitled “on the stability of random systems, etc.”, The Journal of the Acoustical Society of America 32(10), 1356 (1960)
  * (19) R.Z. Khas’ minskii, Necessary and sufficient conditions for the asymptotic stability of linear stochastic systems, Theory of Probability & Its Applications 12(1), 144 (1967)
  * (20) F. Kozin, S. Prodromou, Necessary and sufficient conditions for almost sure sample stability of linear ito equations, SIAM Journal on Applied Mathematics 21(3), 413 (1971)
  * (21) F. Kozin, A survey of stability of stochastic systems, Automatica 5(1), 95 (1969)
  * (22) F. Kozin, Some results on stability of stochastic dynamical systems 14, 163 (1986)
  * (23) E.F. Infante, On the Stability of Some Linear Nonautonomous Random Systems, Journal of Applied Mechanics 35(1), 7 (1968)
  * (24) L. Arnold, G. Papanicoaou, V. Wihstutz, Asymptotic analysis of the lyapunov exponent and rotation number of the random oscillator and application, Siam Journal on Applied Mathematics 46, 427 (1986)
  * (25) R. Stratonovich, _Topics in the Theory of Random Noise Vol. I: General Theory of Random Processes Nonlinear Transformations of Signals and Noise_ (Gordon & Breach, 1963)
  * (26) R.Z. Khas’ minskii, A limit theorem for the solutions of differential equations with random right-hand sides, Theory of Probability & Its Applications 11(3), 390 (1966)
  * (27) Y. Maruyama, A. Maki, L. Dostal, N. Umeda, Improved stochastic averaging method using hamiltonian for parametric rolling in irregular longitudinal waves, Journal of Marine Science and Technology 27(1), 186 (2022)
  * (28) Y. Liu, L. Liu, L. Dostal, J. Lu, The applicability of stochastic averaging method to solve the ship rolling response excited by narrow-band waves, Ocean Engineering 251, 111109 (2022)
  * (29) Y. Maruyama, A. Maki, L. Dostal, N. Umeda, Application of linear filter and moment equation for parametric rolling in irregular longitudinal waves, Journal of Marine Science and Technology (2022)
  * (30) Y. Maruyama, A. Maki, L. Dostal, N. Umeda, Stochastic assessment of acceleration probability density function for parametric rolling using moment method, Under review. Currently available in arxiv:2209.11942 [math.DS]. pp. 1–9 (2022)
  * (31) S. Ariaratnam, D. Tam, Random vibration and stability of a linear parametrically excited oscillator, Zeitschrift Angewandte Mathematik und Mechanik 59, 79 (1979)
  * (32) A. Maki, Y. Maruyama, L. Yaliu, D. Leo, Estimation of stability boundary of parametric rolling in irregular seas and its stabilization, Autumn Conference proceedings, the Japan Society of Naval Architects and Ocean Engineers 35, 419 (2022)
  * (33) O. Grim, Beitrag zu dem problem der sicherheit des schiffes im seegang, Schiff und Hafen 6, 490 (1961)
  * (34) N. Umeda, M. Sakai, H. Okamoto, Some remarks on simplified operational guidance for parametric rolling, Autumn Conference Proceedings the Japan Society of Naval Architects and Ocean Engineers (2022)
  * (35) A. Maki, K. Hoshino, L. Dostal, Y. Maruyama, F. Hane, N. Umeda, Y. Yoshimura, Stochastic stabilization and destabilization of ship maneuvering motion by multiplicative noise, Journal of Marine Science and Technology under review (2022)
  * (36) W. Eugene, Z. Moshe, On the relation between ordinary and stochastic differential equations, International Journal of Engineering Science 3(2), 213 (1965)
  * (37) J. Bogdanoff, F. Kozin, Moments of the output of linear random systems, The Journal of the Acoustical Society of America 34(8), 1063 (1962)
  * (38) F. Kozin, On almost sure stability of linear systems with random coefficients, Journal of Mathematics and Physics 42(1-4), 59 (1963)
  * (39) T. Caughey, A. Gray Jr, On the almost sure stability of linear dynamic systems with stochastic coefficients, Journal of Applied Mechanics 32, 365– (1965)
  * (40) F. Kozin, Some results on stability of stochastic dynamical systems 14, 163 (1986)
  * (41) F. Kozin, C.M. Wu, On the stability of linear stochastic differential equations, Journal of Applied Mechanics 40(1), 87 (1973)
  * (42) P.E. Kloeden, E. Platen, _Numerical solution of stochastic differential equations / Peter E. Kloeden, Eckhard Platen_ (Springer-Verlag Berlin ; New York, 1992)
  * (43) M. Hamamoto, Y. Kim, K. Uwatoko, Study on ship motions and capsizing in following seas (final report), Journal of the Society of Naval Architects of Japan 1991(170), 173 (1991)
  * (44) L. Arnold, Random dynamical systems, Dynamical systems pp. 1–43 (1995)
  * (45) C.J. Söder, A. Rosén, E. Ovegård, J. Kuttenkeuler, M. Huss, Parametric roll mitigation using rudder control, Journal of marine science and technology 18(3), 395 (2013)