Dynamic Model Identification of Ships and Wave Energy Converters Based on Semi-Conjugate Linear Regression and Noisy Input Gaussian Process

Abstract

Reducing the carbon emissions of ships and increasing the utilization of marine renewable energy are the important ways to achieve the goal of carbon neutrality in ocean engineering. Establishing an accurate mathematical model is the foundation of simulating the motion of marine vehicles and structures, and it is the basis of operation energy efficiency optimization and prediction of power generation. System identification from observed input–output data is a practical and powerful method. However, for modeling objects with different characteristics and known information, a single modeling framework can hardly meet the requirements of model establishment. Moreover, there are some challenges in system identification, such as parameter drift and overfitting. In this work, three robust methods are proposed for generating ocean hydrodynamic models based on Bayesian regression. Two Bayesian techniques, semi-conjugate linear regression and noisy input Gaussian process regression are used for parametric and nonparametric gray-box modeling and black-box modeling. The experimental free-running tests of the KRISO very large crude oil carrier (KVLCC2) ship model and a multi-freedom wave energy converter (WEC) are used to validate the proposed Bayesian models. The results demonstrate that the proposed schemes for system identification of the ship and WEC have good generalization ability and robustness. Finally, the developed modeling methods are evaluated considering the aspects required conditions, operating characteristics, and prediction accuracy.

1. Introduction

The demand for mitigating anthropogenic CO₂ emissions increasingly focuses on the transportation system and energy system. Ship emission reduction and ocean renewable energy power generation are the core methods to achieve carbon neutrality in ocean engineering. The International Maritime Organization has set the first target for limiting carbon emissions for international shipping: to limit emissions by at least 50% by 2050 compared to 2008 [1]. Among various types of renewable energy, wave energy resources are enormous [2] and its intensity is high [3]. Moreover, the carbon emissions of wave energy converter (WEC) are smaller than solar energy generation. For ships, the high precision of ship maneuvering systems plays a crucial role in ship controller design [4] and operation energy efficiency [5], which can significantly reduce the carbon emissions of the shipping industry. A wave energy converter needs an accurate hydrodynamic model to predict power generation efficiency under different wave conditions to be used for operational planning [6]. Various methods have been proposed to construct the hydrodynamic model in naval architecture. Depending on whether prior knowledge and physical laws are used in modeling, the modeling methods can be categorized as white-box modeling, gray-box modeling and black-box modeling methods [7].

White-box modeling is the case in which a model is perfectly known. It needs to predefine the mathematical structure entirely from prior knowledge and physical insight. However, due to the strong nonlinearity of water resistance and the randomness of turbulence [8], it is extremely difficult to establish an accurate white-box model of a marine vehicle or structure. The practical way is to first select the model through certain criteria, and then estimate the parameters in the selected model from observation data with system identification. This modeling method is called gray-box modeling. Specific to marine equipment, the most commonly used approach is to establish the equation according to Newton’s second law and then substitute the fitted regression hydrodynamic force in it.

The traditional way to fit the hydrodynamic force in gray-box model is to expand it into a linear function of velocity. For ship modeling, different parametric model structures, such as Abkowtiz model [9,10], Maneuvering Model-ing Group (MMG) model [11] and Nomoto model [12], have been proposed and validated over the years. The hydrodynamic parameters can be obtained by a captive model test with planar motion mechanism (PMM), computational fluid dynamics (CFD) and free-running tests with system identification [13]. Among the above approaches, the system identification with free-running test has been proven to be a powerful and practical method with lower experiment cost [14]. System identification is a general term for estimating parameters from observed input and output data, which provides a reliable mathematical surrogate model in multiple engineering areas [15]. The least square (LS) [16], extended Kalman filter (EKF) [17] and maximum likelihood (ML) [18] algorithms are introduced to identify the hydrodynamic derivatives and proved the effectiveness. Over the last decade, some new methods, with stronger generalization ability and robustness, based on machine learning have also been applied to the estimation of hydrodynamic parameters. Minimizing the Hausdorff metric with the genetic algorithm (GA) can alleviate the impact of noise-induced problems [14]. Mei et al. [19] introduced model reference and random forest (RM-RF) to model ship dynamic model and validated the proposed scheme with free-running test data. Wang et al. [20] presented nu-Support Vector Machine (v-SVM) to improve the robustness of the algorithm.

In the gray-box modeling of wave energy community, Cummins equation [21] is used to define the hydrodynamic model. Generally, there are two ways to determine the equation. Typically, the hydrodynamic model is predefined as the linear model and solved by the potential flow theory [22], whereby the problem is simplified and linearized through assumptions of small amplitude oscillations. However, the simplified linearizing assumptions are invalid when the WECs have large amplitude motions resulting from energetic waves or sustained wave resonance [23]. An alternative method is to use system identification. The training data can be obtained from CFD simulation or scale experiments in a towing tank [24,25]. A popular method is to estimate the real hydrodynamic force using an EKF observer. EKF assumes that the excitation force can be expressed as the sum of a finite set of harmonic components [26,27].

The gray-box modeling methods mentioned above are all parameterized. Recently, a nonparametric gray-box model has been put forward in some studies, and encouraging results have been obtained [28,29,30]. The model still follows the framework of Newton’s law, and the force element, which is difficult to determine, is directly replaced by a machine learning model of related variables. Wang et al. [28] used SVM to replace the Taylor expansion in Abkowtiz model, and they compared the accuracy and computation speed with parametric gray-box and black-box modeling. Xu and Guedes Soares [29] proposed a nonlinear implicit model with nonlinear kernel-based Least Square SVM for a maneuvering simulation of a container ship in shallow water. The forces and moments in [26] are obtained by a PMM test and then trained as outputs for an SVM model related to speed and water depth. In the study of WEC [30], an observer-based unknown input estimator is proposed to identify the wave excitation force, then a Gaussian Process (GP) model is used to predict the wave excitation force. On the one hand, the nonparametric gray-box model directly substitutes the information of the object itself. On the other hand, compared with linear expansion, it can better fit the hydrodynamic force. Therefore, this method is worth studying and comparing with the experimental data of more devices.

Black-box identification directly maps the input data into a high dimensional space without physical insight, but the chosen model structure is selected from a known family [7]. The data-driven model of time-series avoids the inaccuracy caused by the unmodeled part of the parametric gray-box model, such as low external perturbations and coupling effect among the different degrees of freedom (DOF). Machine learning techniques have recently gained popularity in black-box modeling of maneuverability. Rajesh and Bhattacharyya [31] used neural network model to identify the nonlinear maneuvering of large tankers. Fuzzy artificial neural network (ANN) using generalized ellipsoidal basis function [32] are proposed for black-box modeling of ship motion. The assessment of heave displacement for non-buoyant type WEC is investigated by means of ANN [33]. ANN is utilized to predict the wave surface elevation at the WEC location using measurements of wave elevation at ahead located sensor [34]. However, these NN methods have some drawbacks, such as the need for a large amount of training data and the sensitivity to noise. The kernel methods overcome these problems based on statistical learning theory [35]. The kernel methods, such as SVM [28], the Gaussian process (GP) [36], locally weighted learning (LWL) [37] and kernel ridge regression confidence machine [38], are used for identifying the marine dynamic model. Among them, the GP have stronger robustness and generalization with a priori introduction from Bayesian perspective than other methods.

Recently, Bayesian regression has been successful applied in multiple fields for parameter estimation and black-box modeling. Bayesian method possesses an apparent advantage over modeling with good statistical properties, predictions for missing data and forecasting [39,40,41,42]. Moreover, Bayes’ rule provides a reasonable way to update beliefs in light of training data, and the hyperparameters in the Bayesian scheme have an intuitive meaning [43]. Bayesian regression models can work well in dynamic system modeling with a relatively small number of training data points and noisy output [44]. With regards to parametric gray-box modeling, ship dynamic models based on conjugate and semi-conjugate Bayesian regression (ScBR) are adopted to estimate the hydrodynamic parameters in our previous work [44]. For the black-box modeling, Ariza Ramirez et al. [36] used multioutput GPs to identify the ship dynamic system and showed that the GP scheme has better generalization than recurrent neural network (RNN). A series of Bayesian methods is used to quantify the extremal responses of a floating production storage and offloading (FPSO) vessel in [45]. Shi et al. [46] used GP to predict short-term wave for optimal control of wave energy. A complexity penalty and an automatic regularization are introduced in GP based on Bayesian theory. The complexity penalty makes Gaussian process regression have a far smaller risk of overfitting than neural networks.

Based on the above discussion, it can be known that there are various effective system identification methods in ocean engineering. The parameter estimation methods used for gray box modeling are summarized in [14], and nonparametric methods used for black box-modeling are explored and compared in [38]. However, to the best knowledge of the authors, there are very few studies on modeling and comparison of an object using gray-box and black-box models at the same time. When modeling different marine equipment, how to choose a modeling method and design a regression framework based on the existing prior knowledge is still a difficult problem. To answer this question, semi-conjugate regression (ScBR) used for gray-box modeling [44] and noisy input Gaussian Process (NIGP) [47] are proposed in our previous work. However, these two Bayesian algorithms have only been verified by simulation data of different ships. To further verify the applicability of the algorithms and explore the possibility of applying the Bayesian methods to other marine equipment, experimental data including ships and wave energy devices are used for verification and comparison considering the aspects of prerequisite conditions, accuracy and robustness.

This article contributes to the use of Bayesian regression to identify the nonlinear dynamic model of a container ship and an oscillating buoy WEC with gray-box modeling and black-box modeling. First, the Bayesian regression algorithms, including ScBR and NIGP, are introduced. Then, parametric and nonparametric gray-box modeling and black-box modeling schemes based on Bayesian algorithms are proposed for the ship and WEC, respectively. These proposed schemes are validated and compared using experimental data. Finally, the capabilities and challenges of the proposed models are further discussed.

This paper is organized as follows. Section 2 describes the marine dynamic model. The algorithms of ScBR and NIGP are depicted in Section 3. In Section 4 and Section 5, the identification schemes of the ship and WEC and experimental examples are presented to demonstrate the distinction and effectiveness of the proposed two methods. Section 6 and Section 7 present the main conclusions and a further discussion.

2. Kinematic Model

The classical kinematic model in naval architecture is motivated by Newton’s second law, and the rigid-body kinemics equations can be expressed in vector form as [48]

$$\begin{array}{c}{\mathrm{M}}_{RB}\dot{\mathcal{V}}={\tau }_{RB}-C_{RB}\left(\mathcal{V}\right)\\ {\tau }_{RB}={\tau }_h+{\tau }_{env}+{\tau }_{control}\end{array}$$

(1)

where ${\mathrm{M}}_{RB}$ is the rigid-body inertia matrix; $C_{RB}\left(\mathcal{V}\right)$ is a matrix of rigid-body Coriolis and centripetal terms; and ${\tau }_{RB}$ is a vector of generalized forces containing hydrodynamic water resistance, ${\tau }_h$, environmental forces, ${\tau }_{env}$, and control forces, ${\tau }_{control}$. $\mathcal{V}$ denotes the generalized velocity in 6 DOF. The notation of motion variables is shown in

Table 1. Notation of motion variables.

DOF	Motions	Forces	Linear Velocity	Positions
1	Surge	$F_1$	$u$	$x$
2	Sway	$F_2$	$v$	$y$
3	Heave	$F_3$	$w$	$z$
	Rotations	Moments	Angular Velocity	Rotation Angles
4	Roll	$M_1$	$p$	$\phi$
5	Pitch	$M_2$	$q$	$\theta$
6	Yaw	$M_3$	$r$	$\psi$

.

The marine dynamic model is essentially a nonlinear autoregressive model with an exogenous input (NARX) system, and the predictions are based on the previous measurements of the input signals and output signals [43]. Figure 1 shows the first-order NARX configuration for dynamic system, where $c_k$ is the excitation signals such as propeller speed and rudder angle of the ship and the wave force for the WEC; $y_k$ is the original output, which is polluted by noise, $ϵ$; $z$ denotes the z-transformation; and subscript $k$ stands for time step.

3. Bayesian Regression Framework

3.1. Semi-Conjugate Bayesian Regression

Linear regression is to estimate the hydrodynamic parameters $\beta$ in damping matrix ${\tau }_h$, the form of multiple linear regression as

$$y_t=x_t\beta +{\epsilon }_t$$

(2)

where $t$ = 1, n is time step; $y_t$ denotes the measured response; $x_t$ is a $1\times 𝒸$ row vector of the measured values of $𝒸$ predictors; $\beta$ is a $𝒸\times$ 1 column vector of regression parameters corresponding to the argument; and ${\epsilon }_t$ denotes the random disturbance with a mean of zero and common variance of ${\sigma }^2$.

Bayes theorem treats $\beta$ and ${\sigma }^2$ as random variables belonging to some probability distributions. In the Bayesian analysis process, the probability density functions (PDF) of the parameters are updated by incorporating information from the training data. The posterior PDF from Bayes’ theorem can be given as

$$P\left(parameters|data\right)=\frac{P(data|para)P\left(para\right)}{P\left(data\right)}$$

(3)

Generally, the measured value distribution is the normal-inverse-gamma conjugate model [48], in which $\beta$ obeys the multivariate normal distribution ($\mathcal{N}$) and ${\sigma }^2$ is the inverse gamma ($IG$) distribution. Equation (3) can be abbreviated as follows:

$$\pi \left(\beta ,{\sigma }^2|y,x\right)\propto \mathcal{N}\left(\beta \right)\mathcal{N}\left({\sigma }^2\right){{\displaystyle \prod }}_{t=1}^n\varphi \left(y_t;x_t\beta ,{\sigma }^2\right)$$

(4)

where $\varphi \left(y_t;x_t\beta ,{\sigma }^2\right)$ is the Gaussian probability density with mean $x_t\beta$ and variance ${\sigma }^2$ on $y_t$. The regression model contains two types, conjugate and semi-conjugate Bayesian regression, depending on whether the parameters and noise are independent [48].

Parameters and noise are usually not independent of each other in engineering [49]. The prior distributions of $\beta$ and ${\sigma }^2$ are as follows when $\beta$ and ${\sigma }^2$ are dependent:

$$\begin{array}{c}\beta |{\sigma }^2~N_c\left(\mu ,V\right)\\ {\sigma }^2~IG\left(A,B\right)\end{array}$$

(5)

where $\mu$ denotes the mean value ($𝒸\times 1$ vector), $V$ is the $𝒸\times 𝒸$ diagonal matrix in which each element equals the prior variance factor of ${\beta }_j$ and $IG\left(A,B\right)$ stands for the inverse gamma distribution with shape $A$ and scale $B$.

The conditional posterior distribution of $\beta$ and ${\sigma }^2$ can be obtained:

$$\beta |{\sigma }^2,y,x~N_c({\left(V^{-1}+{\sigma }^{-2}{X}^{T}X\right)}^{-1}\left[{\sigma }^{-2}\left(X^{T}X\right)\widehat{\beta }+V^{-1}\mu \right],{\left(V^{-1}+X^{T}X\right)}^{-1})$$

(6)

$${\sigma }^2|\beta ,y,x~IG(A+\frac{n}{2},{(B^{-1}+\frac{1}{2}SSR\left(\beta \right))}^{-1})$$

(7)

where $X$ is an $n\times c$ matrix of training data and $SSR\left(\beta \right)$ is given by

$$SSR\left(\beta \right)={\displaystyle \sum }_{i=1}^n{\left(y_i-{\beta }^T{x}_i\right)}^2={\left(y-X\beta \right)}^T\left(y-X\beta \right)$$

(8)

The posterior distributions of $\beta$ and ${\sigma }^2$ are not analytically tractable because they are mutually influential. In the present work, the Gibbs sampler [50] is applied to approximate the posterior of $\beta$ and ${\sigma }^2$. The Gibbs sampler is an iterative algorithm based on the Markov chain Monte Carlo method. It constructs a dependent sequence of parameter values whose distribution converges to the target joint posterior distribution. The values of parameters are the mean of the posterior of $\beta$.

In multivariate linear regression, introducing the ${\mathrm{L}}_2$—norm into the algorithm to overcome the problems of multicollinearity and overfitting is a generally accepted and effective method, such as ridge regression. ScBR naturally introduces the norm through prior parameters. These types of parameters in the algorithm framework are called hyperparameters in machine learning. Compared to other algorithms, the hyperparameters of the prior distribution, such as the mean and variance, in the Bayesian approach have a clear and intuitive meaning: The value of the prior mean $\mu$ represents the parameter to be identified, which we subjectively set before the regression is performed. When there is no other prior information about the parameter to be estimated, the mean $\mu$ is usually set to zero. The prior variance is obtained by Bayesian optimization algorithm (BOA). BOA is a powerful global optimization algorithm, which is suitable for scenarios with fewer hyperparameters and slower operations of the objective model for hyperparameter optimization [51]. More details about the ScBR with BOA can be found in our previous work [44].

3.2. Noisy Input Gaussian Process

A GP is a distribution over functions and a generalization of the Gaussian distribution to an infinite-dimensional function space. GP is fully specified by a mean function, $m\left(x\right)$, and a covariance function, $k\left(x,x^{\prime }\right)$, as

$$m\left(x\right)=E\left[f\left(x\right)\right]$$

(9)

$$k\left(x,x^{\prime }\right)=E\left[\left(f\left(x\right)-m\left(x\right)\right)\left(f\left(x^{\prime }\right)-m\left(x^{\prime }\right)\right)\right]$$

(10)

where $E$ stands for the expectation operator.

The objective of system identification using GP is to approximate some function $f\left(x\right)$, which maps a D—dimensional input to a scalar output value, $f$. Some training points $n$, which include c—dimensional inputs, ${\left\{x_t\right\}}_{t=1}^n$ and noisy observations ${\left\{y_t\right\}}_{t=1}^n$, are given. The training data are denoted as the $n\times c$ input, $X$, and the $n\times 1$ output vector, $y$. The training outputs is assumed to be corrupted by noise in the regular Gaussian process (RGP).

$$y=f\left(x\right)+{ϵ}_y$$

(11)

where ${ϵ}_y$ stands for Gaussian white noise, which follows a Gaussian distribution with zero mean and variance ${\mathsf{\sigma }}_y{}^2$. RGP defines a GP prior on the function values,

$$p(f|X)=\mathcal{N}\left(m\left(X\right),k\left(X,X\right)\right)$$

(12)

The likelihood function can be obtained through the above assumptions in place,

$$p\left(y|f,X\right)={\displaystyle \prod }_{t=1}^n\mathcal{N}\left(y_t;f_t,{\sigma }_y{}^2\right)$$

(13)

Then, combining the prior Equation (12) and the likelihood function Equation (13), the posterior probability distribution and predict the function values $f^{\ast }$ can be calculated at a set of test points $X^{\ast }$.

$$\left[\begin{array}{c}{f}^{\ast }\\ y\end{array}\right]~\mathcal{N}\left(\begin{array}{cc}\left[\begin{array}{c}m\left(X^{\ast }\right)\\ m\left(X\right)\end{array}\right],& \left[\begin{array}{cc}K\left(X^{\ast },X^{\ast }\right)& K\left(X^{\ast },X\right)\\ K(X,X^{\ast })& K+{\sigma }_y{}^{2}I\end{array}\right]\end{array}\right)$$

(14)

RGP predictive equations are given as

$$p\left(f^{\ast }|X^{\ast },X,y\right)=\mathcal{N}\left(𝓂,𝓈\right)$$

(15)

$$𝓂=m\left(X^{\ast }\right)+K\left(X^{\ast },X\right){[K\left(X,X\right)+{\sigma }_y{}^{2}I]}^{-1}\left(y-m\left(X\right)\right)$$

(16)

$$𝓈=k\left(X^{\ast },X^{\ast }\right)-K\left(X^{\ast },X\right){[K\left(X,X\right)+{\sigma }_y{}^{2}I]}^{-1}K\left(X,X^{\ast }\right)$$

(17)

The RGP only considers output noise ${\sigma }_y{}^2$, while the input data are assumed to be noise-free. However, the output noise will be passed to the input in the NARX models, as shown in Figure 1. Incoming input noise will reduce the accuracy of prediction of GP. Noisy input GP, proposed to solve this problem in [52], takes into account the input noise and posterior data. NIGP further assumes that the inputs are also noisy, and the actual inputs and outputs are labeled $\tilde{x}$ and $\tilde{y}$, respectively.

$$x=\tilde{x}+{ϵ}_x$$

(18)

where ${ϵ}_x$ stands for Gaussian white noise with zero mean and variance ${{\displaystyle \sum }}_x$. Because each input dimension is independently corrupted by noise in this model, ${{\displaystyle \sum }}_x$ is diagonal. Similar to Equation (11), the output function can be written as:

$$y=f\left(\tilde{x}+{ϵ}_x\right)+{ϵ}_y$$

(19)

A first-order Taylor series expansion of the GP latent function, $f$, is used to approximate Equation (19) as,

$$y=f\left(x\right)+{ϵ}_x^T\frac{\partial f\left(\tilde{x}\right)}{\partial \tilde{x}}+{ϵ}_y$$

(20)

The derivative of one GP mean function in Equation (20) is denoted as ${\partial }_{\overline{f}}$, a c-dimensional vector. ${△}_{\overline{f}}$, an $n\times c$ matrix, stands for the value of the derivative for the $n$ functions.

The prior of NIGP is the same as RGP, $p(f|X)=\mathcal{N}\left(0,K\left(X,X\right)\right)$, where $K\left(X,X\right)$ is the $n\times n$ training data covariance matrix. The predictive posterior mean and variance can be obtained, as

$$\mathbb{E}\left[f^{\ast }|X,y,X^{\ast }\right]=K\left(X^{\ast },X\right){[K\left(X,X\right)+{\sigma }_y{}^{2}I+diag\left\{{△}_{\overline{f}}{{\displaystyle \sum }}_x{△}_{\overline{f}}^T\right\}]}^{-1}y$$

(21)

$$\begin{array}{cc}\mathbb{V}\left[f^{\ast }|X,y,X^{\ast }\right]=& k\left(X^{\ast },X^{\ast }\right)-K\left(X^{\ast },X\right)[K\left(X,X\right)+{\sigma }_y{}^{2}I\\ & +diag\left\{{△}_{\overline{f}}{{\displaystyle \sum }}_x{△}_{\overline{f}}^T\right\}{]}^{-1}K\left(X,X^{\ast }\right)\end{array}$$

(22)

where the notation “$diag$” results in a diagonal matrix.

In this way, the input is regarded as deterministic and a correction term, $diag\left\{{△}_{\overline{f}}{{\displaystyle \sum }}_x{△}_{\overline{f}}^T\right\}$, is added to the output noise. In other words, the influence of the input noise depends on the slope of the function we are approximating. Compared to the RGP, the NIGP introduces extra hyperparameters ${ϵ}_x$ per input dimension. A major advantage of this method is that these hyperparameters can be trained alongside any others by ML. The marginal likelihood function of the NIGP is given as,

$$-log p_{NIGP}\left(y|X,\theta \right)=\frac{1}{2}\log \left|K_n\right|+\frac{1}{2}{\left(m\left(X\right)-y\right)}^Tℬ+\frac{N}{2}log2\pi$$

(23)

where

$$K_n=K\left(X,X\right)+diag\left\{{△}_{\overline{f}}{{\displaystyle \sum }}_x{△}_{\overline{f}}^T\right\}+{\sigma }_y{}^{2}I$$

(24)

$$ℬ=K_n{}^{-1}\left(m\left(X\right)-y\right)$$

(25)

A two-step iteration approach is used to estimate all the hyperparameters. In Step 1, a regular GP without input noise is trained by gradient descent algorithm, and the hyperparameters except input noise can be obtained. In Step 2, the derivatives are calculated and used to approximate the posterior distribution. The marginal likelihood of the GP with corrected variance can also be obtained, and the process continues until the convergence. The computation involves chaining the derivatives of the marginal likelihood back through the slope calculation. A complete explanation of NIGP can be found in [47,48,49,50,51,52,53,54] and some supplementary notes are in Bijl’s study [55].

The squared exponential (SE) covariance function is used in NIGP, expressed as

$$k\left(x_i,x_j\right)={\sigma }_f{}^2\exp (-\frac{1}{2}{\left(x_i-x_j\right)}^T\wedge \left(x_i-x_j\right)$$

(26)

where ${\sigma }_f$ is the amplitude and $\wedge$ denotes a diagonal matrix of the squared length-scale parameter.

4. Identification of Marine Craft

4.1. Parametric Gray-Box Modeling

The essence of the parametric gray-box modeling is to construct a simplified parameterized equation to replace Equation (1). Figure 2 shows the 6 DOF reference frames for the ship. The nondimensional rigid-body kinetics using the Prime system of surface ship 3 DOF maneuvering motion is given as follows:

$$\begin{array}{c}\left[\begin{array}{ccc}{m}^{\prime }-X_{\dot{u}}^{\prime }& 0& 0\\ 0& m^{\prime }-Y_{\dot{v}}^{\prime }& m^{\prime }{x}^{\prime }{}_G-Y_{\dot{r}}^{\prime }\\ 0& m^{\prime }{x}_G^{\prime }-N_{\dot{v}}^{\prime }& I^{\prime }{}_{zz}-N_{\dot{r}}^{\prime }\end{array}\right]\left[\begin{array}{c}\dot{u\prime }\\ \dot{v}\prime \\ \dot{r}\prime \end{array}\right]=\left[\begin{array}{c}{F}_1^{\prime }\\ F_2^{\prime }\\ M_3^{\prime }\end{array}\right]\end{array}$$

(27)

where $m$ stands for the ship mass; $x_G$ stands for the longitudinal coordinate of the ship’s center of gravity in the body-fixed coordinate frame; $I_z$ stands for the moments of inertia of the ship about the ${\mathrm{z}}_0$-axes; $X_{\dot{u}}^{ }$, $Y_{\dot{v}}^{ }$, $Y_{\dot{r}}^{ }$, $N_{\dot{v}}^{ }$ and $N_{\dot{r}}^{ }$ stand for acceleration derivatives; and $F_1$ $F_2$ and $M_3$ stand for forces and moment disturbing quantity at ${\mathrm{x}}_0$-axis,${ \mathrm{y}}_0$-axis and ${\mathrm{z}}_0$-axis, respectively. The superscript “$\prime$” indicates that the corresponding variable is normalized using the Prime-system.

Model complexity and model capacity should be balanced when selecting the mathematical model for identification. The most widely used model is the Abkowitz model, a Taylor-series expansion model. It has good generalization performance but includes many coefficients. Some of the coefficients have no physical meaning. A simplified Abkowitz 3-DOF model [56] is employed to construct the white-box model as it contains fewer hydrodynamic parameters while ensuring high accuracy, which can suppress parameter drift caused by too many variables [57]. The nonlinear forces and moment are defined as:

$$\begin{array}{c}{F}_1{}^{\prime }=X_{hydro}·A\left(i\right)\\ F_2{}^{\prime }=Y_{hydro}·B\left(i\right)\\ M_3{}^{\prime }=N_{hydro}·C\left(i\right)\end{array}$$

(28)

where the hydrodynamic derivatives and speed state variables in Equation (28) are as follows:

• $X_{hydro}={[X_u^{\prime },X_{vv}^{\prime },X_{rr}^{\prime },X_{\delta \delta }^{\prime }, X_{vr}^{\prime },X_{v\delta }^{\prime },X_{r\delta }^{\prime },X_0^{\prime },]}_{1\times 8}$
• $Y_{hydro}= {[Y_v^{\prime },Y_r^{\prime },Y_{\delta }^{\prime },Y_{v\left|v\right|}^{\prime },Y_{v\left|r\right|}^{\prime },Y_{\left|r\right|r}^{\prime },Y_{r\left|v\right|}^{\prime },Y_{\delta \delta \delta }^{\prime },Y_{vv\delta }^{\prime },Y_{v\delta \delta }^{\prime },Y_{r\delta \delta }^{\prime },Y_{rr\delta }^{\prime },Y_{rv\delta }^{\prime },Y_0^{\prime }]}_{1\times 14}$
• $N_{hydro}={[N_v^{\prime },N_r^{\prime },N_{\delta }^{\prime },N_{v\left|v\right|}^{\prime },N_{v\left|r\right|}^{\prime },N_{\left|r\right|r}^{\prime },N_{r\left|v\right|}^{\prime },N_{\delta \delta \delta }^{\prime },N_{vv\delta }^{\prime },N_{v\delta \delta }^{\prime },N_{r\delta \delta }^{\prime },N_{rr\delta }^{\prime },N_{rv\delta }^{\prime },N_0^{\prime }]}_{1\times 14}$
• $A\left(i\right)=[u_a^{\prime }\left(i\right),{v^{\prime }}^2\left(i\right),{r^{\prime }}^2\left(i\right),\dots , r^{\prime }\left(i\right){\delta }^{\prime }\left(i\right) ,1 ]{ }_{1\times 8}^T$
• $B\left(i\right)={[v^{\prime }\left(i\right),{r^{\prime }}^2\left(i\right),{\delta }^{\prime }\left(i\right),v^{\prime }\left(i\right)\left|v^{\prime }\left(i\right)\right|,\dots , r^{\prime }\left(i\right)v^{\prime }\left(i\right){\delta }^{\prime }\left(i\right) ,1]}_{1\times 14}^T$
• $C\left(i\right)={[v^{\prime }\left(i\right),{r^{\prime }}^2\left(i\right),{\delta }^{\prime }\left(i\right),v^{\prime }\left(i\right)\left|v^{\prime }\left(i\right)\right|,\dots , r^{\prime }\left(i\right)v^{\prime }\left(i\right){\delta }^{\prime }\left(i\right) ,1]}_{1\times 14}^T$

with the relative speed $u_a=u-u_{nom}.$ As seen, there are 36 hydrodynamic parameters in the simplified Abkowitz model. Euler’s stepping method is utilized to discretize the equation of motions. The constructor of samples for hydrodynamic parameter estimation can be obtained as follows:

Input variables:$\left[A\left(i\right),B\left(i\right),C\left(i\right)\right]$

Output response:

$$\left[\begin{array}{c}\left(m^{\prime }-X_{\dot{u}}^{\text{'}}\right)L\frac{u_a{}^{\text{'}}\left(i+1\right)-u_a{}^{\text{'}}\left(i\right)}{U\left(i\right)△t}\\ \left(m^{\prime }-Y_{\dot{v}}^{\text{'}}\right)L\frac{v^{\prime }\left(i+1\right)-v^{\prime }\left(i\right)}{U\left(i\right)△t}+\left(m^{\prime }{x}_G^{\text{'}}-Y_{\dot{r}}^{\text{'}}\right)L\frac{r^{\prime }\left(i+1\right)-r^{\prime }\left(i\right)}{U\left(i\right)△t}\\ \left(m^{\prime }{x}_G^{\text{'}}-N_{\dot{v}}\right)L\frac{v^{\prime }\left(i+1\right)-v^{\prime }\left(i\right)}{U\left(i\right)△t}+\left(I_{zz}^{\text{'}}-N_{\dot{r}}\right)L\frac{r^{\prime }\left(i+1\right)-r^{\prime }\left(i\right)}{U\left(i\right)△t}\end{array}\right]$$

(29)

where $U=\sqrt{u^2+v^2}$ is the resultant speed in the horizontal plane and $△t$ is the time sample.

The procedure of the parametric gray-box modeling and motion prediction using ScBR is briefly depicted in Figure 3. A Bayesian optimization algorithm (BOA) is employed to tune the value of prior variance, $V$, in semi-conjugate regression. For more details regarding the use of semi-conjugate regression with BOA to identify the parameters, please refer to our previous work [44].

4.2. Black-Box Modeling

A continuous-time black-box model directly describes the relationship between the input variables and out response without any constrains. The principal parameters and the mathematical model are not required in the black-box modeling. The structure of the training data follows the form

Input variables: $\left[u\left(i-1\right),v\left(i-1\right),r\left(i-1\right),\delta \left(i-1\right)\right]$

$$\mathrm{Output} \mathrm{response}: \left[u\left(i\right),v\left(i\right),r\left(i\right)\right]$$

(30)

Figure 4 shows the process of black-box modeling and motion prediction using NIGP. The SVM is also used with the same training data for comparison with Bayesian regression. The RBF kernel function in Equation (31), with an automatic kernel scale ${\sigma }_{SVM}$, is used to train the SVM.

$$k\left(x_i,x_j\right)=\exp (-\frac{{\left|x_i-x_j\right|}^2}{2{\sigma }_{SVM}})$$

(31)

The hyperparameters in SVM are trained by BOA with the “Bayesopt” function in Matlab. In theory, this scheme can overcome the drawbacks of parametric gray-box models, such as a failure to represent the actual behavior of the system due to unmodeled components.

4.3. A Case Study of a Large Container Ship

KVLCC2 is a scale model of large tankers. It is one of the benchmark ships for verification and validation of ship maneuvering simulation methods recommended by Simulation Workshop for Ship Maneuvering (SIMMAN) [58]. Maneuvering and course maintaining tests with the KVLCC2 models were performed at the Hamburg Ship Model Basin (HSVA). The main dimensions of the scale ship are detailed in

Table 2. Particulars of KVLCC2.

Elements	Full-Scale	Model
$L_{pp}$ (m)	320.0	7.0
$B$ (m)	58.0	1.1688
$D$ ($\mathrm{m}$)	30.0	0.6563
Displacement (${\mathrm{m}}^3$)	312,622	3.2724
Draught (m)	20.8	0.4550
Beam coefficient	0.8098	0.8098
Nominal speed (m/s)	7.97	1.18
Rudder speed ($\dot{\delta }$)	2.3 deg/s	15.8 deg/s
Nondim mass ($m^{\prime })$	$1908\times {10}^{-5}$
Nondim x coordinate of CG ($x^{\prime }{}_G$)	3486$\times {10}^{-5}$
Nondim inertia in yaw ($I^{\prime }{}_z$)	119$\times {10}^{-5}$

.

Here, 35°/5° zigzag maneuver data with a cumulative time of 180 s are used for training the parametric gray-box model using ScBR, and the sample time is 0.5 s. The hyperparameters of ScBR, prior variance $V$, are tuned by BOA. The posterior hydrodynamic parameters estimated by ScBR are listed in

Table 3. The nondimensional hydrodynamic parameters for ScBR and SVM (1 × 10⁻⁵).

X-Coef.	ScBR	SVM	Y-Coef.	ScBR	SVM	N-Coef.	ScBR	SVM
$X_u^{\prime }$	−140.1	−128	$Y_v^{\prime }$	350.4	−94	$N_v^{\prime }$	−44.8	−54.9
$X_{vv}^{\prime }$	152.6	175	$Y_r^{\prime }$	1936.0	2066	$N_r^{\prime }$	−125.5	−82.9
$X_{rr}^{\prime }$	−180.0	−118	$Y_{\delta }^{\prime }$	568.3	486	$N_{\delta }^{\prime }$	−180.9	−146.8
$X_{\delta \delta }^{\prime }$	125.2	−116	$Y_{v\left|v\right|}^{\prime }$	68.7	63	$N_{v\left|v\right|}^{\prime }$	5.4	−5.4
$X_{vr}^{\prime }$	−328.2	−303	$Y_{v\left|r\right|}^{\prime }$	128.5	67	$N_{v\left|r\right|}^{\prime }$	−3.3	2.6
$X_{v\delta }^{\prime }$	245.2	196	$Y_{\left|r\right|r}^{\prime }$	932.6	737	$N_{\left|r\right|r}^{\prime }$	−48.5	−30.8
$X_{r\delta }^{\prime }$	−584.2	−455	$Y_{r\left|v\right|}^{\prime }$	30.8	177	$N_{r\left|v\right|}^{\prime }$	−14.3	−5.5
$X_0^{\prime }$	−144.0	−85	$Y_{\delta \delta \delta }^{\prime }$	216.3	−58	$N_{\delta \delta \delta }^{\prime }$	−65.2	−52.9
			$Y_{vv\delta }^{\prime }$	98.7	29	$N_{vv\delta }^{\prime }$	−9.6	−8.0
			$Y_{v\delta \delta }^{\prime }$	41.0	17	$N_{v\delta \delta }^{\prime }$	1.6	−6.3
			$Y_{r\delta \delta }^{\prime }$	306.5	−50	$N_{r\delta \delta }^{\prime }$	9.6	6.1
			$Y_{rr\delta }^{\prime }$	314.2	99	$N_{rr\delta }^{\prime }$	12.4	12.9
			$Y_{rv\delta }^{\prime }$	350.4	−40	$N_{rv\delta }^{\prime }$	−44.8	2.1
			$Y_0^{\prime }$	1936.0	−56	$N_0^{\prime }$	−125.5	1.4
Added mass not identified	$X_{\dot{u}}^{\prime }$	$-95.4$	$Y_{\dot{v}}^{\prime }$	$-1283$			$N_{\dot{v}}^{\prime }$	0
			$Y_{\dot{r}}^{\prime }$	0			$N_{\dot{r}}^{\prime }$	−107

. The added mass, including $X_{\dot{u}}^{\prime }$ and $Y_{\dot{v}}^{\prime }$, is calculated by slender-body instead of SI [59]. For comparison of the ScBR, Luo and Li’s results of SVM under the same parameterization gray-box modeling are also listed in the table. It should be noted that the mainstream algorithms for marine equipment identification are offline algorithms, which are usually trained after the data are obtained and then deployed in the controller or simulation system [14]. Therefore, the time spent on tuning the hyperparameters is not mentioned in the article.

Here, 35°/5° and 20°/5° zigzag maneuvers every 5 s are used for training the black-box modeling driven by NIGP and SVM. It is of no application value to predict the training movement of ship by using the model obtained from the training data. To verify the generalization ability of the models identified by gray-box modeling and black-box modeling driven by SVM and Bayesian regression, the 30°/5° and 15°/5° zigzag tests are predicted. Figure 5 and Figure 6 show the prediction results of each method, and the root mean square error (RMSE) is adopted to analyze the prediction performance of these methods, which is shown in

Table 4. Estimation of forecast accuracy by RMSE and computation time for the validation test.

	Parametric Gray-Box Model		Black-Box Model
	SVM	ScBR	SVM	NIGP
30°/5°u	0.053	0.040	0.115	0.094
v	0.182	0.092	0.186	0.121
r	2.530	1.226	2.213	1.834
15°/5°u	0.155	0.240	0.262	0.021
v	0.126	0.163	0.238	0.062
r	0.605	1.294	2.140	0.443
time (s/step)	0.0009		0.004	0.014

. In addition, the computation time of each step of these methods for prediction is also listed in the table. From the validation results, it can be concluded that the trends of all the predictions before 70 s are basically consistent with the experiment. After 70 s, the difference between the predictions results of various methods gradually increased. On the whole, the parametric gray-box and black-box modeling based on Bayesian regression results are in acceptable agreement with the validation samples and show a stronger ability to predict than SVM. From the perspective of the modeling framework, the prediction results of gray-box modeling are better than those of black-box modeling in 30°/5° zigzag maneuvers, but worse in 15°/5° zigzag tests. The main reason is that the training data of gray-box modeling only contain 35°/5° movement, which is closer to the 30°/5° zigzag validation test. For the prediction time, parametric gray-box modeling is significantly faster than black-box modeling because the calculation process of parametric gray-box modeling is entirely linear. Because it considers the input noise and variance in the calculation process, NIGP spends more time on the prediction than SVM. Note that the black-box modeling usually requires more training data to enhance generalization ability than parametric gray-box modeling, because the specified framework of the parametric gray-box model already contains some information about the system. In a similar study [28], four groups of ship maneuver datasets are used for training black-box models while one group dataset is used for parameter estimation.

5. Identification of WEC

5.1. Nonparametric Gray-Box Modeling

Similar to the ship model in Equation (27), the time domain 3 DOF model of the WEC buoy is given as,

$$\left[\begin{array}{ccc}m-X_{\dot{u}}^{\prime }& 0& 0\\ 0& m^{\prime }-Z_{\dot{w}}^{\prime }& m^{\prime }{y}_G^{\prime }-Z_{\dot{q}}^{\prime }\\ m{z}_G& -m{x}_G& I_{zz}\end{array}\right]\left[\begin{array}{c}\dot{u}\\ \dot{w}\\ \dot{q}\end{array}\right]=\left[\begin{array}{c}{F}_1^{\prime }\\ F_3^{\prime }\\ M_2^{\prime }\end{array}\right]$$

(32)

Different from the parametric model in Equation (28), the force and moment on the right side of the equation are not fitted by the method of multiplying the hydrodynamic coefficient and the speed. In this case, NIGP is adopted to perform nonlinear regression among forces, speed and other variables. The training sample that couples hydrodynamic forces and moment nonlinear regression for training NIGP is

Input variables: $\left[u\left(i\right),w\left(i\right),q\left(i\right),h\left(i\right)\right]$

Output response:

$$\left[\begin{array}{c}\left(m-X_{\dot{u}}^{ }\right)\frac{u\left(i+1\right)-u\left(i\right)}{△t}\\ \left(m-Z_{\dot{w}}^{ }\right)\frac{w\left(i+1\right)-w\left(i\right)}{△t}+\left(m{y}_G-Z_{\dot{q}}^{ }\right)\frac{q\left(i+1\right)-q\left(i\right)}{△t}\\ m{z}_G\frac{u\left(i+1\right)-u\left(i\right)}{△t}-m{x}_G\frac{w\left(i+1\right)-w\left(i\right)}{△t}+I_{zz}\frac{q\left(i+1\right)-q\left(i\right)}{△t}\end{array}\right]$$

(33)

The process of nonparametric gray-box modeling and motion prediction using NIGP is depicted in Figure 7.

5.2. Black-Box Modeling

In the same way as the black-box modeling of the ship, only the time series of motion state variables and wave height are used to train the NIGP model. The structure of the training data follows the form

Input variables: $\left[u\left(i-1\right),w\left(i-1\right),q\left(i-1\right),h\left(i-1\right)\right]$

$$\mathrm{Output} \mathrm{response}: \left[u\left(i\right),w\left(i\right),q\left(i\right)\right]$$

(34)

The detailed process of the black-box modeling of the WEC using NIGP is shown in Figure 8.

5.3. A Case Study of a Multi-Freedom Buoy WEC

The experiment was carried out in the wave tank of Ocean University of China [60], as shown in Figure 9. The model had three independent DOF: surge, heave and pitch. The main dimensions and parameters of the model are listed in

Table 5. Particulars and test conditions of the WEC.

Elements	Value
$\mathrm{Water} \mathrm{depth}$(m)	1.0
$\mathrm{Wave} \mathrm{heihgt}$ (m)	0.2
$\mathrm{Radius}$ ($\mathrm{m}$)	0.4
$\mathrm{Draft}$ ($\mathrm{m}$)	0.4
Height (m)	0.12
Spring stiffness coefficient	85 N/m
Mass (kg)	58
Inertia in yaw ($I_{zz}$)	2.2

. The NDI Optotrak Certus 3D investigator was used to collect the buy’s motion data. The mass of the sliding frame was 58 kg. A spring was used to provide the restoring force for the motion of surge.

The added mass can be calculated as

$$\begin{array}{c}{m}_{\infty }=A\left(\omega \right)+\frac{1}{\omega }{{\displaystyle \int }}_0^{\infty }K\left(t\right)\sin \left(\omega t\right)dt\\ K\left(t\right)=\frac{2}{\pi }{{\displaystyle \int }}_0^{\infty }B\left(\omega \right)\cos \left(\omega t\right)dt\end{array}$$

(35)

where $\omega$ is the wave frequency and $B\left(\omega \right)$ is the radiation damping matrix. The values of $A\left(\omega \right)$ and $B\left(\omega \right)$ were calculated by the ANSYS AQWA software package (AQWA-LINE suite), which implements a boundary element method algorithm.

The wave period of the experimental data is 1.6 s. The first 30 s of the experimental data were used to train each model, and the last 15 s of the data were used as the test set to verify the accuracy of the identified models. The sampling interval of training data for nonparametric gray-box modeling is 0.05 s, while the black-box modeling is 0.1 s. Figure 10 shows that, except for the nonparametric gray-box model driven by SVM, the predicted results of the other three methods are almost the same as the experimental values. The rest motion data with a wave period of 1.8 s were used to further verify the identified models, as presented in Figure 11. It should be noted that the WEC buoy has different added mass in a different wave frequency, so we recalculated the added mass with the wave period 1.8 s and substituted it into Equation (31) for gray-box modeling prediction. Figure 11 shows that the trend of the experimental and Bayesian gray-box and black-box modeling prediction fit well in the motion of surge and heave. However, it can be observed that there is some discrepancy between the prediction and the experiment in pitch. This may be mainly due to fact that the frequency of wave in the training data is higher than that of the test, and the motion is very regular, which means that the training data do not fully reflect the dynamic characteristics of the device. In the 1.8 s wave period, the prediction of the gray-box modeling based on SVM failed, and its motion state was significantly slower than the experiment. The RMSE of u, w and q and computation time of the models are listed in

Table 6. Estimation of forecast accuracy by RMSE and computation time for the WEC motion with the wave period 1.6 s and 1.8 s.

	Nonparametric Gray-Box Model		Black-Box Model
	SVM	NIGP	SVM	NIGP
T = 1.6 su	/	0.866	1.883	2.184
w	/	0.580	0.610	1.143
q	/	0.151	0.346	0.404
T = 1.8 su	/	1.503	3.872	1.142
w	/	0.979	1.935	0.620
q	/	0.396	0.664	0.211
time (s/step)	0.0012	0.086	0.0013	0.025

. Table 6 demonstrates that the black-box model based on NIGP is the most accurate identification method for WEC buoy.

6. Discussions

In this work, three different identification frameworks, parametric gray-box modeling, nonparametric gray-box modeling and black-box modeling based on Bayesian regression, are developed. The main objective is to propose a robust and widely used identification methodology for hydrodynamic models of marine equipment for carbon neutralization using experimental data. The Bayesian regression approach was compared with SVM on a KVLCC2 and WEC buoy model and showed good generalization ability.

The relative strengths and weakness of each method are summarized in

Table 7. Capabilities and challenges of Bayesian gray-box modeling and black-box modeling.

Property	Parametric Gray-Box Model Driven by ScBR	Nonparametric Gray-Box Driven by NIGP	Black-Box Model Driven by NIGP
Modeling framework	Newton’s second law equation with Taylor expansion forces	Newton’s second law equation with nonparametric forces	High-dimensional mapping of time series
Required prior knowledge	weak	fair	strong
Nonlinearities	fair	strong	strong
Training with limited data	strong	fair	weak
Noise robustness	weak	fair	strong
Execution time	strong	weak	fair

. For different modeling objects and characteristics, the corresponding modeling method should be selected according to their capabilities. For conversional ship, choosing traditional parametric modeling can produce good results under the limited data conditions. For new types of vehicles such as unmanned surface vehicle (USV) and remotely operated vehicle (ROV), as well as other irregularly shaped marine structures, nonparametric modeling could be a better choice. However, when very few training data exist, the parametric gray-box modeling method can provide a useful model with the help of prior knowledge such as the added mass of the marine equipment. The obtained experimental data is usually the velocity obtained by MRU (motion reference unit) or displacement data measured with a camera. If the force of the device can be obtained by CFD simulation or directly measured by a PMM test, the nonparametric gray-box modeling with nonlinear fluid dynamics would be a very effective method. In terms of the practicality of the algorithm, compared with SVM and ANN, Bayesian regression introduces a prior into the loss function, which has stronger generalization ability. Moreover, NIGP shows stronger predictive ability because of its additional processing of input noise. However, it needs to be acknowledged that it requires a longer execution time due to the complicated calculations in nonparametric modeling.

7. Conclusions

In this paper, we explore three modeling methods using Bayesian methods for high-dimensional marine systems. Taking the ship maneuvering model and the WEC device as examples, this work shows how to use Bayesian regression to design and train gray-box and black-box models of marine dynamic systems for prediction. The results show that the proposed scheme is more robust than ordinary SVM and has the potential to be further applied to other marine equipment. Finally, the characteristics and advantages of these methods are summarized to facilitate international managers and scientists to choose suitable modeling methods for different modeling objects.

Although the preliminary application of the proposed Bayesian methods seems encouraging thus far, the work needs further extension and investigations. (1) For a model calculation to be used in the practical application of control design, the training dataset should be richer and obtained from more abundant excitation signal, to make the identification model more accurate. The experimental data used in this article are not from experiments specially designed for system identification, so the excitation signal of the training data is not enough. Especially for the wave energy device, compared with regular waves, the motion data under irregular waves (such as Jonswap spectral waves) can better reflect its dynamic feature. (2) The experiments of ship and WEC presented in this article were all carried out in water tank, but the equipment in the ocean will be affected by various factors such as wind, water depth and current. Further study is required to introduce these factors as inputs into nonparametric modeling. (3) The WEC used in the article only uses waves as the excitation signal. To obtain better power generation effects, the WEC used in the article only uses waves as the excitation signal for motion prediction. It is better to model PTO damping as the control signal of the device and apply it to power generation control. (4) Model predictive control (MPC) based on GP allows the direct assessment of residual model uncertainly to enable cautious control. It is very interesting to integrate NIGP-based nonlinear nonparametric modeling into MPC for marine systems.