Efficient Estimation of Temporal Exceeding Probability for Ship Responses in Broadband Wave Fields

TL;DR

This paper introduces an efficient method to estimate the probability that ship motion exceeds a threshold in broadband wave fields, significantly reducing computational costs compared to traditional simulations.

Contribution

It develops a new wave group parameterization and incorporates parameter variability into the estimation, enhancing applicability to general broadband wave fields.

Findings

Achieves about 15% error with only 210 samples

Reduces computational cost by over 2300 times compared to full simulations

Sequential sampling improves accuracy over random or Latin hypercube sampling

Abstract

In this paper, we develop an efficient method to evaluate the temporal exceeding probability of ship motion (percentage of time for the ship motion to be above a given high threshold) in an irregular wave field. Our method builds on our previous work of Gong et al. (Applied Ocean Research, 2022), which converts the calculation into a sampling problem in the space of wave group parameters, within which an acquisition-based sequential sampling method is developed to reduce the number of required samples. Two critical advancements are achieved in this paper relative to that work. (1) We develop a new wave group parameterization method, which allows the framework to be applied to general broadband wave fields. (2) We incorporate the variability regarding each parameterized wave group (e.g. varying wave form and initial condition of the ship encountering the group) into the final estimation of a single value of the temporal exceeding probability. Our complete framework is tested for a ship subject to a wave field with a JONSWAP spectrum, for different ship motion dynamical models, spectral bandwidths, and exceeding thresholds. The results show that for most cases our method provides a result with O(15%) error or below within 210 samples, with the ground truth obtained from a continuous simulation that is more than 2300 times more expensive than our method. We also demonstrate the benefits of sequential sampling (with an acquisition function updated due to (2)) compared to standard random or Latin hypercube (LH) samplings, in terms of the mean error of the results.

Tables

Table 1. Table 1 : Parameter values explored in the case studies for validating the proposed method.

Parameter	Values
Peak enhancement factor $\gamma$ (–)	1, 3
Ship-motion dynamical model $({\beta }_2,{\epsilon }_1)$ ($s^{-2}{\mathrm{rad}}^{-2},m^{-1}{s}^{-2}$)*	(–0.1, 0.008), (–0.2, 0.008), (–0.1, 0.016)
Exceeding threshold $r_s$ ($\mathrm{rad}$)	0.30, 0.35
Wave-parametrization threshold $\Delta \eta$ ($m$)	5, 6, 7

Equations

$$P_{\text{temp}}=\frac{1}{T_{\text{end}}}{\int }_0^{T_{\text{end}}}1(|r(t)|-r_s)\mathrm{d}t,$$

$$1(x)=\{\begin{array}{ll}1,& x>0\\ 0,& x\le 0\end{array}.$$

$$P_{\mathrm{temp}}^{a}\;=\;\frac{1}{\,T_{\mathrm{end}}\,}\sum_{i=1}^{m}S\bigl{(}l_{i},a_{i}\bigr{)}\;=\;\frac{m}{\,T_{\mathrm{end}}\,}\,\int\;S(l,a)\;\,P_{LA}(l,a)\,\mathrm{d}l\,\mathrm{d}a,$$

$$S\bigl{(}l_{i},a_{i}\bigr{)}\;=\;\int\;\mathbf{1}\!\bigl{(}\bigl{|}\,r\bigl{(}t;l_{i},a_{i}\bigr{)}\bigr{|}\;-\;r_{s}\bigr{)}\,\mathrm{d}t,$$

$$S(l,a)\sim \text{GP}\left({\mu }_S(l,a),{\sigma }_S^2(l,a)\right),$$

$$P_{\mathrm{temp}}\;=\;\frac{m}{T_{\mathrm{end}}}\;\int\;\mathbb{E}_{\omega}\bigl{[}S(l,a,\omega)\bigr{]}\;\,P_{LA}(l,a)\,\mathrm{d}l\,\mathrm{d}a,$$

$$\bar{S}(l,a)|D\sim \text{GP}\left({\mu }_{\bar{S}}(l,a|D),{\sigma }_{\bar{S}}^2(l,a|D)\right),$$

$$\delta (\omega )\sim N(0,{\sigma }_0^2).$$

$$h(l,a,\omega )=\{\begin{array}{ll}\frac{S(l,a,\omega )}{l},& \text{if }S(l,a,\omega )>0\\ \frac{r_{\max }(l,a,\omega )-r_s}{r_s},& \text{if }S(l,a,\omega )=0\end{array},$$

$$\bar{h}(l,a)|D\sim \text{GP}\left({\mu }_{\bar{h}}(l,a|D),{\sigma }_{\bar{h}}^2(l,a|D)\right),$$

$$\delta (\omega )\sim N(0,{\sigma }_0^2).$$

$$U|D=\frac{m}{T_{\text{end}}}\int \left(E_{\omega }[S^{+}(l,a,\omega )|D]-E_{\omega }[S^{-}(l,a,\omega )|D]\right)P_{LA}(l,a)\mathrm{d}l\mathrm{d}a,$$

$${\mathbb{E}}_{\omega }[S^{\pm }(l,a,\omega )\mid D]$$

$$\mu_{\text{tr}}^{\pm}(l,a\mid D)=\bar{h}^{\pm}(l,a)-\sigma_{0}\frac{\varphi\!\Bigl{(}\frac{1-\bar{h}^{\pm}(l,a)}{\sigma_{0}}\Bigr{)}-\varphi\!\Bigl{(}-\frac{\bar{h}^{\pm}(l,a)}{\sigma_{0}}\Bigr{)}}{\Phi\!\Bigl{(}\frac{1-\bar{h}^{\pm}(l,a)}{\sigma_{0}}\Bigr{)}-\Phi\!\Bigl{(}-\frac{\bar{h}^{\pm}(l,a)}{\sigma_{0}}\Bigr{)}}.$$

$$l^{\ast },a^{\ast }=\arg \underset{l,a}{\max }\left(E_{\omega }[S^{+}(l,a,\omega )|D]-E_{\omega }[S^{-}(l,a,\omega )|D]\right)P_{LA}(l,a).$$

$$\ddot{r}+{\alpha }_1\dot{r}+{\alpha }_2\dot{r}|\dot{r}|+\left({\beta }_1+{\epsilon }_1\cos (\theta )\eta (t)\right)r+{\beta }_2{r}^3={\epsilon }_2\sin (\theta )\eta (t),$$

$$\text{NMAE}=\frac{1}{N}\sum _{i=1}^N\left|\frac{{\hat{P}}_{\text{temp}}^i-P_{\text{temp}}^e}{P_{\text{temp}}^e}\right|,$$

$$y=f(X)+R,R\sim N(0,{\gamma }_0^2).$$

$$p(f(X)|D)=\frac{p(f(X),\mathbf{y})}{p(\mathbf{y})}=\text{GP}({\mu }_f(X),{\mathrm{cov}}_f(X,X^{\prime })),$$

$${\mu }_f(X)=k_f(X,\mathbf{X}{)}^{\top }{\left(K_f(\mathbf{X},\mathbf{X})+{\gamma }_0^{2}I\right)}^{-1}\mathbf{y},$$

$${\mathrm{cov}}_f(X,X^{\prime })=k_f(X,X^{\prime })-k_f(X,\mathbf{X}{)}^{\top }{\left(K_f(\mathbf{X},\mathbf{X})+{\gamma }_0^{2}I\right)}^{-1}{k}_f(\mathbf{X},X^{\prime }),$$

$$k_f(X,X^{\prime })={\tau }^2\frac{1}{\Gamma (\nu )2^{\nu -1}}{\left(\sqrt{2\nu }\text{dist}(X,X^{\prime })\right)}^{\nu }{K}_{\nu }\left(\sqrt{2\nu }\text{dist}(X,X^{\prime })\right).$$

$$\text{dist}(X,X^{\prime })={\left((X-X^{\prime }{)}^{\top }{\Lambda }^{-1}(X-X^{\prime })\right)}^{1/2},$$

Full text

Keywords: Extreme ship responses, Broadband wave fields, Bayesian experimental design, Gaussian process regression

1 Introduction

The assessment of rare events in complex dynamical systems is crucial for reliability and safety in engineering design Mohamad and Sapsis (2018). In the marine context, one of such problems is the extreme ship motion induced by adverse sea conditions, such as parametric rolls, which poses challenge in the design and operation of ships Mohamad (2017); Spyrou et al. (2008); Umeda et al. (2004). These rare events, although occurring infrequently, can have catastrophic consequences, necessitating accurate estimation of their probabilities for informed design and decision-making processes.

Estimation of the statistics of wave-induced extreme motions is a challenging problem due to the stochastic nature of the wave field and the low probability of these extreme events. A brute-force calculation using a high-fidelity method (such as CFD) is usually computationally prohibitive, since a very long simulation of ship-wave interaction is needed to obtain sufficient exposure of the extreme events Belenky et al. (2012). Over the years, many reduced‐order methods have been developed to enable the calculation, such as the Envelope POT approach McTaggart (2000); Campbell et al. (2016); Weems et al. (2019), the split‐time approach Belenky et al. (2012, 2016, 2018), and the critical‐wave‐group approach Themelis and Spyrou (2008); Anastopoulos and Spyrou (2016, 2019); Anastopoulos et al. (2016). In spite of their success in reducing the computational cost, heuristic components are more or less involved which affects the robustness of their probability estimation in realistic situations.

More recently, a new group of methods have been proposed for estimation of the extreme ship motion probability Gong et al. (2022); Mohamad and Sapsis (2018); Cousins and Sapsis (2016); Farazmand and Sapsis (2017); Gong et al. (2021); Gong and Pan (2022b); Tang and Adcock (2022); Gong and Pan (2024, 2022a). These methods convert the probability estimation problem into a sampling problem in a low-dimensional parameter space, which describes the encountering wave condition associated with input probability. A sequential sampling algorithm based on a Gaussian process regression (GPR) surrogate model is then developed to reduce the number of necessary samples (each associated with a ship motion simulation) for estimation of extreme ship motion probability. Following the above framework, such methods have been first developed for the group-based motion probability Mohamad and Sapsis (2018); Gong et al. (2021); Tang and Adcock (2022) (i.e., distribution of largest motion within a group), and then for temporal exceeding probability (i.e., percentage of time for the ship motion to be above a given high threshold) Gong et al. (2022). As discussed in Gong et al. (2022), the temporal exceeding probability, which we will continue to explore in this paper, provides a more robust and general probabilistic description of extreme motions than the group-based counterpart. While this group of methods are developed on a more rigorous foundation than previous ones, their applications are mostly restricted to narrow-band wave fields. This limitation stems from the wave field parameterization: as shown in Figure 1, a narrow-band wave field consists of consecutive Gaussian-like wave groups that can be naturally parameterized by the group length $l$ and height $a$. However, a broadband wave field lacks such structures since more frequency components (other than simply a modulation of the main frequency) are involved (see Figure 2).

One work that deals with the broadband wave field is Guth and Sapsis (2022), which parameterizes the wave field using Proper Orthogonal Decomposition (POD) instead of the Gaussian wave group approach. Although POD naturally captures the structure of a broadband wave field, the dimension of the parameter space increases up to $20$-$30$ compared to only $2$ ($l$ and $a$) in the narrow-band cases. For such a high-dimensional parameter space, sequential sampling is more difficult (or at least significantly more expensive) and is thus not implemented in Guth and Sapsis (2022). In addition, the number of dimensions needed for the input parameter space depends on the spectral bandwidth, with the relation not clearly known a priori.

In this paper, we extend our methodology in Gong et al. (2022) to enable the evaluation of temporal exceeding probability in broadband wave fields. A key feature of the developed method is that we are able to retain a two-dimensional parameterization of the wave field (in contrast to up to 20-30 dimensions as in Guth and Sapsis (2022)). This unique feature has been made possible due to two critical developments: (a) we apply a new wave parameterization method inspired by Bassler et al. (2010), which parameterizes the wave field by length $l$ and height $a$ of wave groups exceeding a pre-selected threshold. (b) We approximate the randomness of the response due to variability associated with each $(l,a)$ group by a Gaussian distribution, and incorporate its effect into the final estimation of the temporal exceeding probability. Here the “variability” includes detailed wave form and initial condition of a ship encountering this wave group, the former associated more with the broadband wave field, and the latter general for wave fields of any spectral width but not handled in Gong et al. (2022). A sequential sampling method is correspondingly developed considering (b), in particular in terms of tracking the randomness of the response function as a part of the posterior of the GPR. We note that in the whole process, the group threshold in (a) serves as a free parameter whose optimal value is not known a priori. Our extensive tests show that a choice of it close to half of the significant wave height $H_{s}$ produces satisfactory results for all cases.

We test our method for broadband wave fields with JONSWAP spectra (characterized by small peak enhancement parameters), with the ship response calculated from a nonlinear roll model. Different spectral bandwidths, exceeding thresholds, and parameters in the roll model are chosen, resulting in a total of 36 test cases. We show that, in most cases, our method achieves a relative error of O($15\text{\,}\mathrm{\char 37\relax}$) or below within $210$ samples, compared to the ground truth obtained from one continuous simulation that is typically over $2300$ times more expensive. Moreover, we illustrate the advantages of the developed sequential sampling method over traditional random and Latin hypercube (LH) sampling methods, in terms of faster convergence of the result with the number of samples. We finally remark that the coupling of the current method with high-fidelity ship motion simulations (such as CFD) is straightforward, as has been demonstrated in Gong et al. (2022). This will not be demonstrated again in the current paper.

The remainder of this paper is organized as follows. In Section 2, we provide the problem setup on temporal exceeding probability. In Section 3, we present the methodology, especially the new wave group parameterization and the sequential sampling framework. The validation of the method is demonstrated in Section 4, with Section 5 concluding the paper.

2 Problem setup

Considering a ship navigating through a wave field with wave-induced responses, say roll motions, we are interested in evaluating the temporal exceeding probability given by:

[TABLE]

where $r(t)$ is the time series of roll angle resulted from wave signal $\eta(t)$, $r_{s}$ is a high threshold above which the roll is considered extreme, $t$ is time and $T_{\text{end}}$ is the total exposure time, and $\mathbf{1}(\cdot)$ is the Heaviside (indicator) function defined as:

[TABLE]

Therefore, Eq. (1) describes the percentage of time that the ship roll response exceeds a high threshold given sufficient exposure time.

A direct evaluation of Eq. (1) involves a computation of roll motion $r(t)$ from wave signal $\eta(t)$ in time $[0,T_{\text{end}}]$. This can become computationally prohibitive if high-fidelity simulations (such as CFD) are employed, since $T_{\text{end}}$ needs to be very large to cover enough exceeding events associated with a given wave spectrum. Our goal is to develop an efficient method (free from the long computation of $r(t)$) to enable the evaluation of Eq. (1) in general wave fields including those with broadband spectra.

3 Method

3.1 Review of the previous method in Gong et al. (2022)

In Gong et al. (2022), we have developed a method to evaluate Eq. (1) for narrow-band wave fields with some other simplifications. The method starts by describing the wave elevation $\eta(t)$ by consecutive Gaussian groups parameterized by length $l$ and height $a$ (see Figure 1(a)), a procedure only appropriate for waves with a narrow-band spectrum. By collecting $(l,a)$ from all $m$ groups in the wave time series $\eta(t)$, a probability distribution $P_{LA}(l,a)$ can be generated (see Figure 1(b)). The temporal exceeding probability can then be (approximately) computed by:

[TABLE]

where

[TABLE]

is the exceeding time of the $(l_{i},a_{i})$ group, so that the integral in Eq. (4) gives the expected exceeding time in one wave group.

Clearly, in the computation of Eq. (3), the key is to obtain the function $S(l,a)$. Since a direct computation of $S(l,a)$ for all groups can be prohibitively expensive, a GPR (the details of the algorithm and formulae are summarized in Appendix A), in particular the mean value, is used as a surrogate of

[TABLE]

trained through data from limited number of samples in $(l,a)$ space with associated function evaluations. In order to reduce the number of required samples, a sequential sampling algorithm based on Bayesian experimental design (BED) method is developed, which samples at optimal locations to accelerate the convergence of the results to the true $P_{\text{temp}}^{a}$.

We note that $P_{\text{temp}}^{a}$ in Eq. (3) serves only as an approximation to $P_{\text{temp}}$ in Eq. (1) (thus having the superscript $a$ in the former). This is because when using Eq. (3), one has to assume that the group exceeding time is deterministic given $(l,a)$, i.e., $S(l,a)$ is a deterministic function. However, in practice given $(l,a)$, the group exceeding time can vary because both ship initial condition encountering a wave group and finer details of the wave elevation in a group can be different. The former is a problem for any wave field, and the latter becomes more severe for those with a broadband spectrum (of course after a proper group parameterization is developed as will be discussed). A new approach is needed to address this variability problem regarding $S(l,a)$.

3.2 The proposed method

We now introduce our new approach, which is applicable to wave fields with a broadband spectrum and addresses the variability problem regarding $S$.

Inspired by Bassler et al. (2010), a new group parameterization method is developed which targets the critical segments in $\eta(t)$ that may lead to extreme ship responses. Specifically, we first define a single wave as a segment starting from $\eta=0$ with a positive slope and experiencing two subsequent zero-crossings. A wave group is then defined as a series of consecutive single waves with amplitudes above a given threshold $\Delta\eta$ (see Figure 2(a)). Each wave group is then parameterized by its length $l$ and maximum amplitude $a$ within the group, with the collection of all groups providing a two-dimensional parameter space with known probability density $P_{LA}(l,a)$ (Figure 2(b)).

We note that the threshold $\Delta\eta$ serves as a free parameter in the above group detection method. Clearly, if $\Delta\eta$ is too high, some critical wave groups leading to extreme ship motions may be missed. If $\Delta\eta$ is too small, $P_{LA}(l,a)$ may become more complex together with a higher variability of detailed wave elevations for each $(l,a)$, which both negatively impact the convergence of the sequential sampling method. Therefore, a proper value of $\Delta\eta$ is important for the success of the proposed method, but this can only be determined empirically. In section 4, we test values of $\Delta\eta$ around $H_{s}/2$ and show that such choices provide satisfactory results.

While the above parameterization method allows description of the wave field in a low-dimensional space, groups with the same $(l,a)$ parameter may have significantly different detailed elevation shapes (see Figure 3). These different detailed elevations, together with the random initial condition of ship encountering each wave group, lead to a random exceeding time for groups with the same parameters. To account for this randomness, the temporal exceeding probability needs to be formulated as (now dropping the superscript a as this is exact):

[TABLE]

where $S(l,a,\omega)$ with $\omega$ is a random seed accounting for the randomness mentioned above.

We further note that \mathbb{E}_{\omega}\bigl{[}S(l,a,\omega)\bigr{]} in Eq. (6) is quantitatively different from $S(l,a)$ in Eq. (3). This is because when $S(l,a)$ is approximated by Eq. (5) with training data $D$, there is no guarantee that each function evaluation provides the mean value over all randomness (which is especially the case for varying ship initial conditions, see Gong and Pan (2022a)). Our next task is to describe the procedure to approximate $S(l,a,\omega)$ through a surrogate GPR model and a sequential sampling method to reduce the number of required samples.

3.2.1 Surrogate GPR model

Construction of a surrogate model for $S(l,a,\omega)$ is a common procedure in BED. Given a dataset $D=\{(l^{i},a^{i},\omega),S^{i}\}_{i=1}^{n}$ consisting of $n$ inputs of $(l,a)$ and their exceeding time $S$ as responses, one may naively think that a GPR can be placed on the function $S(l,a,\omega)\equiv\bar{S}(l,a)+\delta(\omega)$, i.e.,

[TABLE]

[TABLE]

However, the above GPR approximation is not appropriate due to two reasons. (1) $\delta(\omega)$ does not follow a Gaussian distribution as desired in Eq. (8). To understand this point, simply consider a situation that $\bar{S}(l,a)$ is positive and close to zero, so that value of $\delta(\omega)$ is limited (since $S(l,a,\omega)$ has to be non-negative) and cannot be approximated by a Gaussian distribution. (2) $\bar{S}(l,a)$ is not a smooth function (zero for most cases and positive in a small region with a non-smooth transition) which is difficult to be represented as a Gaussian process. We note that reason (2) is also present in the deterministic case, which has been discussed properly in Gong et al. (2022). To address both issues, we define an auxiliary function (i.e., similar procedure as in Gong et al. (2022) but here we also include the randomness parameter $\omega$):

[TABLE]

where $r_{\text{max}}(l,a,\omega)=\text{max}_{t}|r(t;l,a,\omega)|$ is the group-maximum response in group $(l,a)$. When $S=0$, function $h$ takes negative values with $r_{\text{max}}-r_{s}$ serving as a “negative penalty”, quantifying how far the response is from the threshold. We note that by definition $h(l,a,\omega)$ is a much smoother function than $S(l,a,\omega)$ and is rigorously bounded in $[-1,1]$. Therefore, both issues associated with GPR on $S(l,a,\omega)$ are resolved, and we are ready to place a GPR surrogate model on $h(l,a,\omega)$.

Given the dataset $D=\{(l^{i},a^{i},\omega),h^{i}\}_{i=1}^{n}$ and the decomposition $h(l,a,\omega)=\bar{h}(l,a)+\delta(\omega)$, we can place the GPR approximation on $\bar{h}$:

[TABLE]

with the randomness term:

[TABLE]

where $\mu_{\bar{h}}$ is the posterior mean and $\sigma_{\bar{h}}^{2}$ is the posterior variance. The randomness $\delta(\omega)$ is modeled by a zero-mean Gaussian distribution with constant variance $\sigma_{0}^{2}$, implying a case of homoscedastic randomness. We note that the whole algorithm can be extended to the case of heteroscedastic randomness $\delta(l,a,\omega)\sim N(0,\sigma_{0}^{2}(l,a))$ with techniques developed in Gong and Pan (2022b). However, we find that using homoscedastic randomness is sufficient for the current purpose and we will keep it for the current paper. Once $h$ is obtained through GPR, we can reconstruct $S(l,a,\omega)$ by $S(l,a,\omega)\equiv\mathbf{1}\left(h(l,a,\omega)\right)\text{min}(1,h(l,a,\omega))l$. The min operator is because $S$ cannot physically exceed $l$, meaning that the ship cannot spend more than the entire group length above the threshold.

3.2.2 Acquisition function

Given the function $h$ (and $S$) learned from the GPR based on dataset $D$, $P_{\text{temp}}$ can be considered as a random variable with its randomness resulting from the uncertainty in $h$. The uncertainty in $P_{\text{temp}}$ can therefore be estimated as:

[TABLE]

where $S^{\pm}(l,a,\omega)|D=\mathbf{1}\{h^{\pm}(l,a,\omega)|D\}(\text{min}(1,h^{\pm}(l,a,\omega)|D))l$ with $h^{\pm}(l,a,\omega)|D=\bar{h}^{\pm}(l,a)|D+\delta(\omega)$ being the upper and lower bound of the defined function $h$ in Eq. (9), and $\bar{h}^{\pm}(l,a)|D=\mu_{\bar{h}}(l,a|D)\pm\sigma_{\bar{h}}(l,a|D)$ being the upper and lower bounds of $\bar{h}$ prediction. To analytically compute the expectation $\mathbb{E}_{\omega}[S^{\pm}(l,a,\omega)\mid D]$, we use the following expression:

[TABLE]

where $\Phi(\cdot)$ and $\varphi(\cdot)$ are the standard normal cumulative distribution function (CDF), and standard normal probability density function (PDF), respectively, and $\mu^{\pm}_{\text{tr}}$ is the truncated Gaussian mean given by:

[TABLE]

The expression in Eq.(13) has two terms representing contributions from two different intervals of $h$, where the first term accounts for the contribution when $h$ lies between $[0,1]$, and the second term accounts for the contribution when $h$ exceeds 1. We note that when $h$ is negative, the corresponding $S$ is 0, contributing nothing to the expectation.

Given the uncertainty expressed in Eq. (12), we then select the next sample which is expected to reduce the uncertainty in $P_{\text{temp}}$ significantly, i.e., a sample corresponding to the maximum of the integrand in Eq.(12):

[TABLE]

At each iteration, the acquisition function is optimized to identify the most informative point $(l^{*},a^{*})$ within the parameter space. We then randomly locate a wave group in the original wave field with parameters $(l^{*},a^{*})$ (or sufficiently close to it) to capture the randomness in detailed wave shape. In addition, we run the next ship motion simulation from 1$~T_{p}$ ahead of the $(l^{*},a^{*})$ group to naturally capture the randomness in ship initial condition encountering the $(l^{*},a^{*})$ group (see Gong and Pan ([2022a](https://arxiv.org/html/2508.20271v1#bib.bib11)) for an extensive validation of this approach). Additionally, the simulation is continued for approximately 1$~T_{p}$ after the group to account for the time-delayed responses. This process is then continued such that for each iteration, the sequential sample is placed at a location that contributes most significantly to reducing the uncertainty in $P_{\text{temp}}$. The process ends with the computation of $P_{\text{temp}}$ using Eq. (6). We summarize the full algorithm in Algorithm 1, with the code available on Github111https://github.com/umbrellagong/gpbroad.

4 Validation of the method

In this section, we present the results of the proposed method applied to a wide range of cases with varying wave spectral bandwidth, nonlinear ship motion model parameters, and exceeding thresholds. The results of $P_{\text{temp}}$ are compared to their exact values $P^{e}_{\text{temp}}$ computed using Eq. (1) from simulating the ship motion model with sufficiently long time series of wave excitation. In addition, we compare the performance of the developed sequential sampling approach with other commonly used methods such as random and LH sampling to show the superiority of the developed approach.

4.1 Case setup

To set up our test cases, we consider the nonlinear roll motion of a ship subject to a wave field with a broadband spectrum. Specifically, we use a JONSWAP spectrum with spectral peak period $T_{p}=$15\text{,}\mathrm{s}, significant wave height $H_{s}=$12\text{\,}\mathrm{m}, and a peak enhancement factor $\gamma$ with two chosen values of $\gamma$ =1 and 3 (both realistically broadband). The ship motion is assumed to be modeled by a nonlinear roll equation, which models the ship roll response $r(t)$ due to nonlinear resonance and parametric roll in oblique irregular waves with time series of elevation $\eta(t)$:

[TABLE]

where $\alpha_{1}$ and $\alpha_{2}$ are linear and nonlinear damping coefficients, respectively; $\beta_{1}$ and $\beta_{2}$ are linear and nonlinear restoring coefficients; $\varepsilon_{1}$ and $\varepsilon_{2}$ are parametric and direct wave excitation coefficients; and $\theta$ is the encounter angle between the ship and the wave. These empirical coefficients are set as $\alpha_{1}=0.35$, $\alpha_{2}=0.06$, $\beta_{1}=0.04$, $\theta=\pi/6$, and $\varepsilon_{2}=0.012$, as in Gong et al. (2022). We vary the other two parameters, namely $(\beta_{2},\varepsilon_{1})\in\{(-0.1,0.008),\,(-0.2,0.008),\,(-0.1,0.016)\}$ to represent different levels of nonlinear restoring and parametric excitations in the ship roll dynamics. We note that the simple model in Eq. (16) is used since we need the exact result $P^{e}_{\text{temp}}$ for validation, but our method is not restricted to a particular choice of ship model. For example, one can easily replace Eq. (16) by a CFD calculation as demonstrated earlier in Gong et al. (2022).

In summary, the parameter values in all tested cases are listed in Table 1, which also includes different values of $\Delta\eta$, chosen as values close to $H_{s}/2$ (in particular, with three choices: $H_{s}$ and $H_{s}/2\pm 1$) to benchmark the sensitivity of the method to its free parameter as discussed in section 3.2.

4.2 Results

We start by showing a typical result for the case with $\varepsilon_{1}=0.008$, $\beta_{2}=-0.2$, $\gamma=3$, $r_{s}=$0.35\text{,}\mathrm{rad}, and $\Delta\eta=$5\text{\,}\mathrm{m}. Figure 4 plots the computed $P_{\text{temp}}$ as a function of sampling numbers using sequential sampling, random sampling, and LH sampling, together with its exact value $P^{e}_{\text{temp}}$. All the sampling results are obtained as an ensemble of $100$ repeated experiments, i.e., we repeat the respective sampling calculations $100$ times, each time starting from $10$ different initial LH samples. The results are therefore presented as the mean plus one standard deviation shown on the upper side of the mean curve, to improve the fairness of the comparison. We see in Figure 4 that results from sequential sampling successfully approaches $P_{\text{temp}}^{e}$ in less than $200$ samples, while the results from random and LH samplings fail to capture $P_{\text{temp}}^{e}$ even at the end of $210$ samples. In addition, the results from sequential and LH samplings are associated with similar uncertainty (or standard deviation), both smaller compared to that of random sampling. This is a point which will be further discussed later.

We further characterize the computational cost saving by our method. In the result shown in Figure 4, it takes about $182$ samples for our sequential sampling method to reach within $1\text{\,}\mathrm{\char 37\relax}$ error bar around $P^{e}_{\text{temp}}$. Considering each sequential sample contain a simulation of time series of $68.8\text{\,}\mathrm{s}$ on average, the total simulation time with sequential sampling method is about $12\,520.4\text{\,}\mathrm{s}$. In contrast, the exact calculation of $P_{\text{temp}}^{e}$ requires $28\,876\,800\text{\,}\mathrm{s}$ simulation for the obtained value to reach $1\text{\,}\mathrm{\char 37\relax}$ accuracy, as shown in Figure 5. Therefore, by employing our approach equipped by the new wave group parameterization and sequential sampling, we are expected to save approximately $2300$ times computational resources.

Our next goal is to benchmark the performance of our developed method with all $36$ cases listed in Table 1. For this purpose, we first introduce two metrics to characterize the performance of a method. The first one is the Normalized Mean Absolute Error (NMAE) to describe the relative error of an estimation normalized by the true value, defined as:

[TABLE]

where $\hat{P}_{\text{temp}}^{i}$ is the estimated $P_{\text{temp}}$ using a particular sampling method (from random, LH, and sequential samplings) at the $i$-th trial, and $N$ is the total number of trials. To be consistent with results in Figure 4, we will keep using $\hat{P}_{\text{temp}}^{i}$ obtained at 210 samples and $N=100$. The second metric is the standard deviation of the estimated $P_{\text{temp}}$ normalized by the true $P_{\text{temp}}^{e}$ value (NSTD), also evaluated at $210$ samples with respect to the $N=100$ ensembles.

The values of NMAE and NSTD are plotted in Figures 6 and 7, respectively, for the 36 test cases. In Figure 6 we see that NMAE resulted from our sequential sampling method is consistently the smallest among all three tested sampling methods. For most cases, NMAEs from sequential sampling are below or around $15\text{\,}\mathrm{\char 37\relax}$, except that for cases with $\gamma=1$ and $r_{s}=$0.35\text{,}\mathrm{rad}$$, NMAE can rise up to $20\text{\,}\mathrm{\char 37\relax}$ in Figures 6d and 6e (considering the $P_{\text{temp}}^{e}$ for these cases is O(${10}^{-4}$), this error should still be considered small). This is natural because these cases are associated with the larger spectral bandwidth and higher exceeding threshold, where the former factor leads to more irregular wave patterns (thus groups) and the latter leads to rarer exceedance. For all cases (including those of Figures 6d and 6e), we have verified that if we further increase the sequential sampling numbers, the result may further improve, but the final converged results may be limited by the slight inaccuracies inherent in the group parameterization method. In addition, from all tested cases in Figure 6, we find no consistent preference on the value of $\Delta\eta$. Instead, from the results we argue that all choices of $\Delta\eta$ around $H_{s}/2$ provide reasonable and close estimation for $P_{\text{temp}}$.

The values of NSTD plotted in Figure 7 show a somewhat different trend. In general, the results of LH sampling achieves the smallest value of NSTD among the three sampling methods. This is because LH sampling ensures samples to be uniformly distributed in the $(L,A)$ space with less variability over trials. However, practically the lower NSTD of LH sampling method is not a favorable situation, since the method is meanwhile associated with larger NMAE. Such a situation of lower NSTD and higher NMAE represents a bias in the estimation of $P_{\text{temp}}$ by LH sampling. The sequential sampling, on the other hand, is associated with slightly larger NSTD but significantly smaller NMAE, therefore serving as a much more balanced method for computing $P_{\text{temp}}$.

We finally explore the reason underlying the high efficiency of sequential sampling in computing $P_{\text{temp}}$. Figure 8 shows the sampling locations of the three sampling methods in the space of $(L,A)$. We also plot the contour of the function $E_{\omega}[S(l,a,\omega)]P_{LA}(l,a)$ in the figure, which illustrates the importance of sampling location to the result of $P_{\text{temp}}$ according to Eq. (6). Preferably, we would like regions with higher $E_{\omega}[S(l,a,\omega)]P_{LA}(l,a)$ to be sampled more so that the regions contributing dominantly to $P_{\text{temp}}$ is approximated more accurately in GPR. This is indeed what happens in the sequential sampling, where we see most sample points are concentrated near the peaks of $E_{\omega}[S(l,a,\omega)]P_{LA}(l,a)$ that vary in different cases. In contrast, random sampling and LH sampling lack this feature, with the former concentrating on regions with higher $P_{LA}(l,a)$ and the latter uniformly distributed in the $(L,A)$ space. Therefore, the favorable efficiency of sequential sampling to compute $P_{\text{temp}}$ can be expected from this visualization of sampling locations.

5 Conclusion

In this paper, we continue to develop our efficient method in evaluating the temporal exceeding probability $P_{\text{temp}}$ of ship responses subject to an incoming wave field, i.e., the fraction of time that the response exceeds a specified threshold. The current work improves upon our previous paper Gong et al. (2022) in two aspects. (1) It extends the calculation into wave fields of general spectra, including the realistically broadband cases. (2) It captures the variability of response for given group parameters $l$ and $a$, which originates from the variable wave form and ship encountering initial conditions. Both improved features are integrated into our general framework including wave group parameterization and sequential sampling enabled by GPR and a newly-developed acquisition function. We test the performance of our new framework in many cases with different ship motion models, wave spectral bandwidth and exceeding thresholds. We show that in all cases our method achieves favorable results, with an error of O($15\text{\,}\mathrm{\char 37\relax}$) or below in most of the cases and a saving of computational cost of $2300$ times compared to the brute-force computation.

ACKNOWLEDGEMENT

This research is supported by the Office of Naval Research, USA grant N00014-24-1-2266.

Appendix A Gaussian process regression

We consider the task of inferring the input-to-response (ItR) function from a dataset $D=\{X^{i},y^{i}\}_{i=1}^{n}$ consisting of $n$ inputs $\boldsymbol{X}=\{X^{i}\in\mathbb{R}^{d}\}_{i=1}^{n}$ (e.g., $X=(l,a)$ in section 3.2.1) and the corresponding outputs $\boldsymbol{y}=\{y^{i}\in\mathbb{R}\}_{i=1}^{n}$ (e.g., $y=h$ in section 3.2.1).

GPR assumes the observed output equals the latent function $f(X)$ plus i.i.d. Gaussian noise with constant variance $\gamma_{0}^{2}$ at all $X$:

[TABLE]

A prior, representing our beliefs over all possible functions we expect to observe, is placed on $f$ as a Gaussian process $f(X)\sim\text{GP}(0,k_{f}(X,X^{\prime}))$ with zero mean and covariance function $k_{f}$. Following the Bayes’ theorem, the posterior prediction for $f$ given the dataset $D$ can be derived to be another Gaussian:

[TABLE]

with analytically tractable mean $\mu_{f}(X)$ and covariance $\operatorname{cov}_{f}(X,X^{\prime})$:

[TABLE]

[TABLE]

where $K_{f}(\boldsymbol{X},\boldsymbol{X})$ is the covariance matrix with elements $[K_{f}(\boldsymbol{X},\boldsymbol{X})]_{ij}=k_{f}(X_{i},X_{j})$, and $I$ is the identity matrix.

For the covariance function $k_{f}$, we use the Matérn kernel with smoothness parameter $\nu=1.5$ defined as:

[TABLE]

Here, $K_{\nu}$ is the modified Bessel function of the second kind, $\Gamma(\nu)$ is the gamma function, and the distance function $\text{dist}(X,X^{\prime})$ is computed by:

[TABLE]

where $\tau$ and diagonal matrix $\Lambda$ are the kernel hyperparameters representing the characteristic amplitude and length scales. The GPR hyperparameters $\theta=\{\tau,\Lambda,\gamma_{0}\}$ are determined by maximizing the likelihood function $p(D|\theta)\equiv p(\boldsymbol{y}|\theta)=N(0,K_{f}(\boldsymbol{X},\boldsymbol{X})+\gamma_{0}^{2}I)$.