- Email: [email protected]
Hydrodynamic parameter identification for ship manoeuvring mathematical models using a Bayesian approach
Ocean Engineering 195 (2020) 106612
Contents lists available at ScienceDirect
Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng
Hydrodynamic parameter identification for ship manoeuvring mathematical models using a Bayesian approach Yifan Xue a, Yanjun Liu a, b, *, Chen Ji a, Gang Xue a a b
Institute of Marine Science and Technology, Shandong University, Qingdao, 266000, China School of Mechanical Engineering, Shandong University, Jinan, 250061, China
A R T I C L E I N F O
A B S T R A C T
Keywords: Ship manoeuvrability Hydrodynamic coefficients System identification Bayesian multilinear regression White noise
A novel identification approach to nonlinear ship manoeuvring models based on Bayes’ rule is presented. The empirical Bayesian method is used to clean simulated polluted responses from a 20� /20� zigzag test. Two effi cient Bayesian models, conjugate regression and semi-conjugate regression, are adopted for hydrodynamic parameter identification. To obtain other prior parameters, the Bayesian optimization (BO) algorithm is intro duced in the Bayesian regression model. By using the identified model and an experimental model based on Blanke’s nonlinear four-degrees-of-freedom (4-DOF) model, 10� /10� zigzag motion and 35� turning circle manoeuvring are performed. A comparison between the predicted results and the test results demonstrates that both Bayesian regression models have good generalization ability for identification, and that the conjugate Bayesian model has greater prediction accuracy.
1. Introduction The highly accurate mathematical modelling of vessels plays a vital role in ship design and operation. For ships with low centre heights, such as container ships, ship manoeuvring motion is usually accompanied by a remarkable amount of roll motion. Son and Nomoto (Son and Nomoto, 1981) studied container ship motion based on experimental tests and established a four-degrees-of-freedom (4-DOF) equation of motion. Blanke.et (Blanke and Jensen, 1997; P� erez and Blanke, 2002) proposed a nonlinear 4-DOF model based on a roll planar motion mechanism (RPMM) facility and validated it via full-scale trials. There are many hydrodynamic parameters in the ship motion model. Added mass in a ship manoeuvring model can be estimated accurately (Perez and Fossen, 2006) by slender body theory. Hence, the key to the prediction accuracy is to obtain the hydrodynamic derivatives, which are regarded as the most uncertain part in the mathematical model. Many researchers have devoted themselves to solving this challenging problem and indicated that system identification (SI) theory in combi nation with full-scale trials or free-running tests (Padilla et al., 2014) is a powerful and practical method. Classic SI techniques are widely applied for ship and underwater vehicle manoeuvring, including least square (LS) estimation (Nagumo and Noda, 1967), improved least squares (Holzhüter, 1990), the
extended Kalman filter (EKF) (Abkowitz and Martin, 1980; Tiano et al., 2007), the recursive prediction error method (RPE) (Blanke et al., 1986) and ridge regression (Yoon and Rhee, 2003). However, these methods have some inherent disadvantages, such as sensitivity to noise or simultaneous drift (Zhang and Zou, 2013; Luo and Li, 2017). In order to overcome these defects, a novel algorithm has been proposed to build a purely “input-output” model on the basis of artificial neural networks (ANN) (Moreira and Guedes Soares, 2003, 2012; Woo et al., 2018). The advantage of ANN is that it is not necessary to prioritize one complex mathematical ship model. However, ANN cannot reflect any physical background. Moreover, the construction of the network structure is blind before full training. In recent years, several robust machine learning methods such as the back propagation neural network (Rajesh and Bhattacharyya, 2008) and support vector machines (SVM) (Wang et al.,2015; Wang et al.,2014) have been proposed and have achieved encouraging results in the application of hydrodynamic parameter identification. However, the hyperparameters introduced in these al gorithms, such as initial weights or penalty factors, have no intuitive physical meaning or logical explanation. The SI methods mentioned above assume that unknown hydrody namic parameters are constants and use limiting relative frequencies to define probability. Roberts.et (Roberts et al., 1994) proposed using Markov process theories to solve the linear regression problem in ship
* Corresponding author. Institute of Marine Science and Technology, shandong University, Qingdao, 266000, China. E-mail address: [email protected] (Y. Liu). https://doi.org/10.1016/j.oceaneng.2019.106612 Received 8 May 2019; Received in revised form 11 September 2019; Accepted 18 October 2019 Available online 25 November 2019 0029-8018/© 2019 Published by Elsevier Ltd.
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
2. Ship mathematical model To describe the motion characteristics of a vessel, as shown in Fig. 1, two coordinate reference frames are adopted, including the Earth-fixed coordinates E X0 Y0 Z0 and the body-fixed coordinates O x0 y0 z0 . According to Son and Nomoto’s nonlinear model (Son and Nomoto, 1981), the uniform equations, including the xG surge, sway, yaw and roll coupled motion, can be transformed as follows: 8 � m u_ vr xG r2 þ ZG rp ¼ FX > > < _ ¼ FY mðv_ þ ur xG r_ ZG pÞ (1) > IZZ r_ þ mxG ðv_ þ urÞ ¼ MN > : Ixx mZG ðv_ þ urÞ ¼ MK ρgrGZ ðφÞ where m is the ship’s mass; u; v; r; and p denote the surge velocity, sway velocity, yaw rate and roll rate, respectively; and zG are the coordinates of the ship’s centre of gravity in the body-fixed coordinate frame; Ixx and Izz are the moments of inertia of the ship about the x0 and y0 axes, respectively; FX and FY are the longitudinal and lateral force compo nents, respectively; MN is the yaw moment about the z0 axis; and MK is the roll moment about the x0 axis. In the last term of Eq. (1), r denotes the ship displacement, g is the gravity constant, ρ is the mass density of the water and Gz ðφÞ is the buoyancy function, which can be approxi mated as: � � 1 Gz ðφÞ ¼ GM þ BMtan2 ðφÞ sinðφÞ (2) 2
Fig. 1. Coordinate frames for a ship.
roll motion from the perspective of probability and statistics. However, the case of one degree of freedom cannot fully prove the applicability of these theories to the overall movement of the ship. Iseki and Terada introduced Bayesian modelling technology to estimate directional wave spectra based on ship motion data and obtained excellent results (Iseki and Terada, 2002). Bayesian methods provide an alternative viewpoint that treats the parameters as random variables and can numerically represent a set of subjective rational beliefs. Bayes’ rule offers a reasonable way to update beliefs in light of new observations. In addi tion to their formal interpretation as a means of induction, Bayesian methods have significant advantages in parameter estimation with good statistical properties, predictions for missing data, forecasting (Hoff, 2009), etc. Recently, Bayesian methods, especially in parameter esti mation, have been applied in multiple fields (Khoshravesh et al., 2015; Petra et al., 2017; Zitzmann et al., 2016) owing to the prevalence of the Bayesian framework and the applicability of modern inference methods based on the Markov chain Monte Carlo (MCMC) method. This article attempts to apply the Bayesian approach to identifying the hydrodynamic parameters in Blanke’s nonlinear 4-DOF model. First, based on the zigzag test simulated by the 4-DOF ship manoeuvring mathematical model (Blanke and Jensen, 1997; P�erez and Blanke, 2002), Gaussian white noise is added to each observed state to simulate the data obtained in the sea trial. An empirical Bayesian (EB) wavelet method (Silverman and Johnstone, 2004) is implemented for removing noise. Then, two efficient Bayesian regression models are utilized to estimate the parameters, and the Bayesian hyperparameter optimization (BO) algorithm is introduced to determine the coefficients of the prior distribution. Finally, in order to validate the proposed method, we use the posterior distribution of the obtained parameters to predict the steering motion of the ship and compare the results with RPMM experimental data.
where GM is the metacentre height, BM is the distance from the centre of buoyancy to the metacentre, and φ is the roll angle. For the convenience of comparison, expanding the force and moment in Eq. (1) with Taylor series (Blanke and Jensen, 1997; P�erez and Blanke, 2002) and rewriting them in a nondimensional matrix using the prime system of SNAME (SNAME, 1950) for substituting model test results yields: 2 3 2 0 3 0 0 2 3 F m X u_ 0 0 0 u _ 6 7 6 10 7 0 0 0 0 0 0 0 0 6 7 6 7 6 7 0 m Y m x Y m Z Y v_ r_ 6 76 v_ 7 6 F2 7 G G 6 74 5 ¼ 6 0 7 (3) 0 0 0 0 0 0 r _ 6 7 6 7 0 m x N I N N v_ r_ p_ G ZZ 4 5 4 M1 5 0 0 0 0 p_ 0 0 m ZG Kv_ Kr_ Ixx K p_ M2 (For notational convenience, the first matrix containing the mass, inertia coefficients and acceleration derivatives in Eq. (3) is labelled J .) The manoeuvring model of Eq. (3) can be discretized using Euler’s method: V
ðiþ1Þ
¼V
ðiÞ
þ
F ⋅h J
(4) T
where V ¼ ½u; v; r; p�T , F ¼ ½F1 ; F2 ; M1 ; M2 � , h is the interval time, and i and i þ 1 are the adjacent sampling time steps. 0
F’1 ¼ Xhydro ⋅AðiÞ þ m’ z’G r’ p’ F’2 ¼ Yhydro ⋅BðiÞ þ m’ u’ r’ M ’1 ¼ Nhydro ⋅CðiÞ þ m’ x’G u’ r’ M ’2 ¼ Khydro ⋅DðiÞ m’ z’G u’ r’
2
0
m’ x’G r’ 2
0
0
m’ v’ r’ (5)
m’ v’ r’
ρgr’ Gz ðφÞ
’
Ocean Engineering 195 (2020) 106612
Y. Xue et al.
where the hydrodynamic derivatives and speed state variables in Eq. (5) are as follows: # " ’ ’ X u ; X uu ; X ’uuu ; X ’vr ; X ’rr ; X ’v ; X ’vv ; X ’vφ ; X ’φ ; X ’φφ ; X ’pp ; X ’ppu ; X ’δ ; X ’δδ ; Xhydro ¼ X ’δu ; X ’δδu ; X ’vδ ; X ’vδδ 1�18
" Yhydro ¼
Nhydro ¼
�
π β; σ2 jy; x ¼ R
� �T BðiÞ ¼ v’ ðiÞ; v’ 2 ðiÞ; v’ ðiÞjv’ ðiÞj; …; δ’2 ðiÞu’a ðiÞ; δ’ 3 ðiÞu’a ðiÞ 1�28
ℓðβ; σ2 jy; xÞπðβÞπ ðσ2 Þ π ðβÞπðσ2 Þρðβ; σ 2 jy; xÞdβd σ2
(10)
where prior distributions are typically denoted πðβÞ and πðσ2 Þ. However, the ultimate goal of regression is to obtain the marginal posterior distribution of β. For different types of prior distributions selected, various Bayesian regression methods such as the hierarchical linear model (Woltman et al., 2012) and the Bayesian lasso model (Park and George, 2008) have been proposed to suit different scenarios. Considering the number of parameters, calculation efficiency and ac curacy, we choose priors so that the corresponding distribution is in the same family of distributions, called the conjugate distributions. Ac cording to the central limit theorem, most of the measured value dis tributions can be approximated by a normal distribution or a Gaussian distribution. A popular choice is the normal-inverse-gamma conjugate model (Robert and Christian, 2014), in which β obeys the multivariate normal distribution (N ) and σ 2 is the inverse gamma (IG) distribution. Eq. (10) can be abbreviated as follows: �Yn � � (11) π β; σ2 jy; x ∝N ðβÞN σ 2 ϕ yt ; xt β; σ 2 t¼1
� �T CðiÞ ¼ v’ ðiÞ; v’ 2 ðiÞ; v’ ðiÞjv’ ðiÞj; v’ ðiÞjr’ ðiÞj; …; δ2 ðiÞu’a ðiÞ; δ3 ðiÞu’a ðiÞ 1�28 � �T DðiÞ ¼ v’ ðiÞ; v’ 2 ðiÞ; v’ ðiÞjv’ ðiÞj; …; δ2 ðiÞu’a ðiÞ; δ’3 ðiÞu’a ðiÞ 1�28 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2 þ v2 ; Unom is the
3. Bayesian regression formulation The object of hydrodynamic parameter identification is to estimate β in the multiple linear regression (MLR) model, as in Eq. (6): (6)
where t ¼ 1, …n denotes time; yt is the observed response; xt is a 1� c row vector of the observed values of c predictors; β is a c� 1 column vector of regression parameters corresponding to the variables that consist of the columns of xt ; and εt is the random disturbance that has a mean of zero and common variance of σ 2 . The Bayesian approach to estimating MLR models treats β and σ 2 as random variables rather than fixed quantities. Generally, the Bayesian analysis process updates the probability distributions of the parameters by incorporating information about the parameters from the observing data, which is a kind of incremental learning. In other words, resolving the Bayesian inverse problem is a posterior probability density function (PDF). Bayes’ theorem gives the posterior PDF as
where ϕðyt ; xt β; σ 2 Þ is the Gaussian probability density with mean xt β and variance σ2 on yt . The regression model is divided into conjugate and semi-conjugate regression(Murphy, 2012) depending on whether parameters and variance are independent. 3.1. Conjugate Bayesian regression When β and σ 2 are independent, the priori distributions of β and σ2 are as follows: � � β�σ2 � Nc μ; σ 2 V (12) σ 2 � IGðA; BÞ
(7)
Eq. (7) can be simplified to posterior ∝ likelihood � prior
1�28
parameter by adjusting the variance of the prior distribution. After specifying the type and parameters of the prior distribution, the joint posterior density can be obtained by substituting the density probability into Eq. (7):
�T � AðiÞ ¼ u’a ðiÞ; u’a 2 ðiÞ; u’a 3 ðiÞ; …; v’ðiÞδ’ðiÞ ; v’ðiÞδ’ 2 ðiÞ 1�18
PðdatajparametersÞPðparametersÞ PðdataÞ
1�28
#
K ’pujpuj ; K ’vφ ; K ’vφφ ; K ’φvv ; K ’0 ; K ’0u ; K ’r ; K ’δ ; K ’δδ ; K ’δδδ ; K ’δv ; K ’δvv ; K ’δu ; K ’δδu ; K ’δδδu
PðparametersjdataÞ ¼
(9)
#
K ’v ; K ’vv ; K ’vjvj ; K ’vjrj ; K ’vrr ; K ’rjrj ; K ’rrr ; K ’rjvj ; K ’rvv ; K ’p ; K ’pp ; K ’ppp ; K ’pu ;
yt ¼ xt β þ εt
� � P yt �xt ; β; σ 2
1�28
N ’v ; N ’vv ; N ’vjvj ; N ’vjrj ; N ’vrr ; N ’r ; N ’rjrj ; N ’rrr ; N ’rjvj ; N ’rvv ; N ’p ; N ’ppp ; N ’pu ;
where the relative surge speed ua ¼ U Unom , U ¼ nominal speed, and δ is the rudder angle.
t¼1
#
Y ’v ; Y ’vv ; Y ’vjvj ; Y ’vjrj ; Y ’vrr ; Y ’r ; Y ’rr ; Y ’rrr ; Y ’rjvj ; Y ’rvv ; Y ’p ; Y ’ppp ; Y ’pu ;
N ’pujpuj ; N ’φ ; N ’vφ ; N ’vφφ ; N ’φvv ; N ’0 ; N ’0u ; N ’δ ; N ’δδ ; N ’δδδ ; N ’δv ; N ’δvv ; N ’δu ; N ’δδu ; N ’δδδu "
Yn
where prior is the parameter distribution that we assume, which allows us to incorporate information about the model before data are imported. Moreover, we can control the confidence in our knowledge about the
Y ’pujpuj ; Y ’φ ; Y ’vφ ; Y ’vφφ ; Y ’φvv ; Y ’0 ; Y ’0u ; Y ’δ ; Y ’δδ ; Y ’δδδ ; Y ’δv ; Y ’δvv ; Y ’δu ; Y ’δδu ; Y ’δδδu "
Khydro ¼
�
ℓ β; σ2 jy; x ¼
where μ is the mean value (c � 1 vector), V is the c � c diagonal matrix in which each element equals the prior variance factor of βj , and IGðA; BÞ denotes the inverse gamma distribution with shape A and scale B. Considering our limited prior information about many parameters, we cannot set different priors for each coefficient, so the diagonal elements in V are simplified to be equal to ➉ in this paper. Since both the prior and
(8)
where likelihood represents the information provided by the sample response about the parameters, and the mathematical model is
3
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
the posterior have the same family of distributions, the conditional posterior distribution of β can be calculated by Eq. (11) and Eq. (13) with σ2 held constant, as shown in Eq. (14): � � � (13) p β�σ 2 ; y; X ∝ pðβÞp yjX; β; σ 2 � � β�σ2 ; y; x �
� Nc
V
1
þ XT X
� 1�
� � XT X b β þ V 1 μ ; σ2 V
1
þ XT X
sampler is an iterative algorithm that constructs a dependent sequence of parameter values whose distribution converges to the target joint posterior distribution. The details of the Gibbs sampler are shown in Algorithm 1.
� 1� (14)
where X is an n � c matrix of training data. The marginal posterior density of σ2 from Eq. (11) is the IG distribution: � � � � 1� � nþc 1 ; B 1 þ SSRðβÞ þ ðb β μÞT V 1 ð b (15) σ 2 ��y; x � IG A þ β μÞ 2 2 where SSRðβÞ is given by n X
SSRðβÞ ¼
yi
β T xi
�2
XβÞT ðy
¼ ðy
(16)
XβÞ
i¼1
To sample from the joint posterior, first draw σ2 from its marginal posterior distribution, given in Eq. (15), and then draw β from its normal conditional posterior distribution, given in Eq. (14), using the simulated value of σ2 . Finally, from the integration of the joint posterior density with respect to σ 2 , the marginal posterior density of β can be calculated as: � � � � � 1� T � � β�y; x � tc V 1 þ X T X X X b β þ V 1 μ ; λ; 2A þ n 2B
1
þ SSRðβÞ þ ðb β
λ¼
h
μÞT V þ X T X
� 1i
1
ðb β
μÞ
Finally, through the two methods above, the forecast of ship motion can be obtained using the mean values of the c -dimensional posterior predictive distribution. 3.3. Hyperparameters in prior distribution A distinctive advantage of the Bayesian approach is that prior in formation and the observed sample are combined in the algorithm. The parameters of the prior distribution, such as the mean and variance, are called hyperparameters in machine learning. In contrast with other commonly used machine learning algorithms, the prior parameters of the Bayesian algorithm have an intuitive meaning: specifying a high prior variance implies that we know very little about the parameter, so a greater weight is given to the information in the data, whereas speci fying a low prior variance implies high confidence in prior information. Based on the feature of Bayesian regression, this paper proposes to use the BO algorithm (Bull, 2011) to determine the parameters (A;B) of the distribution for σ2 , the variance factor of parameters (V) that can reflect the confidence for the mean μ. The algorithm treats the Bayesian
(17)
2A þ n
in which tc ðl; λ; ➈Þ denotes the c-dimensional multivariate t distribution, where l is the location, λ is the scale, and ➈ represents the degrees of freedom. 3.2. Semi-conjugate Bayesian regression When β and σ2 are dependent, the priori distributions of β and σ2 are as follows: � β �σ 2 � Nc ðμ; VÞ (18) σ2 � IGðA; BÞ
Table 1 Physical parameters of the container ship.
As with conjugate regression, a conditional posterior distribution of β and σ 2 can be obtained: � � � � 1� 2 T � � 1� � β�σ2 ; y; x � Nc V 1 þ σ 2 X T X σ X X b β þ V 1 μ ; V 1 þ XT X (19) � � � � n σ 2 ��β; y; x � IG A þ ; B 2
1
1 þ SSRðβÞ 2
� 1�
(20)
Elements
Values
Length (L)
230.66 m
Beam Displacement (r)
32 m 46070 m3
Nominal speed (Unorm )
12.7 m/s
Nominal metacentric height (GM) Transverse metacentre (BM) x coordinate of CG (xG )
0.83 m 9m 0.46 m
Nondim mass (m’ )
750.81 � 10
z coordinate of CG (zG )
However, since β and σ are mutually influential, their posterior distributions are not analytically tractable. Some numerical integration techniques based on the MCMC method have been proposed to solve this problem. In the present work, the Gibbs sampler (Gelfand and Smith, 1990) is applied to approximate the posterior of β and σ2 . The Gibbs 2
Nondim inertia in roll (I’xx )
4
3.54 m 1.30 � 10
5
Nondim inertia in yaw (I’zz )
43.25 � 10
_ Rudder speed (δ)
2.3 deg/s
5
5
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
data and dimensions are given in Table 1. The choice of training manoeuvres is as important as the formulation of a regression model. As Ljung(Ljung, 1987) suggested, the input data for SI should satisfy the persistence of the excitation condition to guar antee convergence and ensure the robustness properties of the identifi cation method. According to his view, the 20� /20� zigzag test and the
regression model as a Gaussian process (GP) and uses the acquisition function (AF) of the expected improvement to find a result that can satisfy the objective function threshold. The basic pseudo-code for Bayesian optimization is shown in Algorithm 2.
Place a Gaussian process prior on Observe
at
While
points according to an initial space-filling experimental design; = do
Update the posterior probability distribution of Let
using available data
be a maximizer of the AF over , where the AF is computed using the current
posterior distribution Observe Increment n end while Return the point evaluated with the largest
Compared with other optimization methods, the BO algorithm works better in cases with fewer hyperparameters and slower operation of the objective model, which coincides with the characteristics of Bayesian regression. Therefore, the BO algorithm is applied in the present work to obtain the hyperparameters of the prior distribution. The training data are pre-processed by Z-score transformation (Altman, 2000) in order to avoid the influence of the high-impact index, and K-fold cross-validation is applied to make full use of the data in this process. K-fold cross-validation randomly divides the set of observations into K groups each time. The first fold is treated as a validation set, and other folds are used for training.
35� turning test reach the appropriate range of extrapolation, unlike the 10� /10� zigzag test, and will therefore provide more informative data for identification. However, the turning circle test is not appropriate for SI because in the motion of the circle after the rudder is turned, the observation values are very similar or even constant, which will cause the vectors in the matrix to be highly correlated. Therefore, a 20� /20� standard zigzag manoeuvre was simulated in MATLAB with 20000-point output discretization by using the hydrodynamic parameters from the RPMM facility at the Danish Maritime Institute (scale ratio 28.75), as given in Table 2 and Table 3. Note that some hydrodynamic derivatives were adjusted to match the full-scale trial results(Blanke and Jensen, 1997). The sampling interval is 0.02 s, and the total time is 400 s. The raw training data couples consist of the following:
4. Numerical example, results and discussion 4.1. Construction of the training model A container vessel was selected as the object model, and the main
input variables : ½AðiÞ;BðiÞ;CðiÞ;DðiÞ� output observed response : 2 m’
� u’ ði þ 1Þ u’ðiÞ X ’u_ L UðiÞh
2
m’ v’ ðiÞr’ ðiÞ
3
m’ x’G r’ ðiÞ þ m’ z’G r’ ðiÞp’ðiÞ
6 7 6 7 6 7 � � ’ ’ ’ ’ ’ ’ 6 7 p ði þ 1Þ p ðiÞ 6 m’ Y ’ �L v ði þ 1Þ v ðiÞ þ m’ x’ Y ’ �L r ði þ 1Þ r ðiÞ 7 ’ ’ ’ L m þ m’u’ðiÞr’ðiÞ z þ Y 6 7 v_ r_ p_ G G UðiÞh UðiÞh UðiÞh 6 7 6 7 6 7 6 7 6 7 � � ’ ’ ’ ’ ’ ’ � � 6 ’ ’ 7 v ði þ 1Þ v ðiÞ r ði þ 1Þ r ðiÞ p ði þ 1Þ p ðiÞ ’ ’ ’ ’ ’ ’ 6 m xG N v_ L 7 þ I þ N þ m L L N x u’ðiÞr’ðiÞ r_ p_ zz G 6 7 UðiÞh UðiÞh UðiÞh 6 7 6 7 6 7 6 7 � p’ ði þ 1Þ p’ ðiÞ 6 7 � v’ ði þ 1Þ v’ ðiÞ r’ ði þ 1Þ r’ ðiÞ � ’ ’ ’ ’ ’ ’ ’ 6 7 m’zG Kv_ L Kr_L þ I xx K p_ L m zG u ðiÞr ðiÞ 6 7 UðiÞh UðiÞh UðiÞh 4 5 þ ρgr0 G’z ðφÞ
5
(21)
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
Table 2 Comparison of the nondimensional hydrodynamic parameters (1 � 10 X-Coeff. X’u_
RPMM 124:4
5
) for surge and sway.
Con
ARE(%)
Semi
ARE(%)
Y-Coeff.
/
/
/
/
Y’v_ Y’r
226.2
226.4
0.05
226.3
0.1
X’uu
64.5
67.5
0.64
64.8
4.6
137.2
147.9
0.86
138.0
7.8
X’vr
24.0
23.9
3.61
24.9
0.4
X’u
X’uuu X’rr
4.4
4.4
14.38
3.8
41.0
24.0
41.51
24.0
41.5
1.0
0.9
34.22
1.4
10.9
X’vφ
108.1
108.5
0.06
108.2
0.3
X’φ
5.9
5.9
0.18
5.9
0.2
X’φφ
42.2
42.6
0.08
42.3
0.9
X’v
X’vv
X’pp
X’ppu
7.2
7.1
0.15
7.2
0.8
3.9
3.9
1.98
3.9
0.0
X’δ
1.4
1.4
0.10
1.4
1.5
X’δδ
116.8
116.8
0.01
116.8
0.0
X’δu ;
17.2
17.0
0.03
17.2
0.9
X’δδu
224.9
224.6
0.12
225.1
0.1
X’vδ
124.5
124.4
0.18
124.7
0.1
X’vδδ
341.0
341.2
0.02
341.0
878:0 48.1
Y’p_
23.3
Y’vv
98.6
725.0
Y’v
0.8
Semi
ARE(%)
/
/
/
/
/
/
/
/
/
/
/
723.3 91.0
0.2 7.7
/ 717.5
73.4
1.0 25.5
5820.0
0.3
5766.6
0.6
1192.7
1200.9
0.7
1196.8
0.3
Y’vrr
1107.9
923.9
16.6
1232.4
11.2
Y’r
Y’rr
118.2
115.4
2.4
125.6
0.0
13.1
/
11.9
6.3 /
Y’rrr
158.0
116.1
26.5
211.5
33.9
Y’rjvj
409.4
414.1
1.2
385.0
6.0
Y’rvv
994.6
787.5
20.8
1081.8
8.8
Y’p
3.4
3.2
5.0
2.1
38.9
Y’ppp
9.3
11.3
21.3
26.3
182.8
Y’pu
23.6
22.6
4.3
78.5
249.5
Y’φ
37.7
37.0
1.8
37.1
1.5
Y’vφ
144.9
150.3
3.7
164.5
13.5
Y’vφφ
2459.3
2327.1
5.4
2405.4
2.2
Y’φvv
177.2
448.9
153.3
288.8
63.0
Y’0
4.7
5.6
19.2
2.4
48.3
5.3
6.4
21.7
5.9
11.3
Y’δ
248.1
248.2
0.0
248.8
0.3
Y’δδ
13.4
13.9
3.8
13.2
1.2
52.5
Y’pujpuj
33.7 632.1
43.0 1104.1
Y’δδδ
193.0
197.7
2.4
186.4
3.4
Y’δv
100.0
98.5
1.5
97.0
3.0
189.2
Y’δvv
233.8
23.6
132.9
29.7
Y’δu
379.4
374.2
1.4
381.0
0.4
Y’δδu
55.6
47.7
14.3
48.3
13.1
232.3
182.8
21.3
240.1
Y’δδδu
It is almost meaningless to use the raw data above for identification directly, as the real experimental measurements will inevitably be affected by various disturbances and produce a great deal of highfrequency responses such as ship oscillatory motion or mechanical noise. Following Sutulo and Soares’s artificially polluted data con struction method (Sutulo and Guedes Soares, 2014), Gaussian white noise is added as follows:
ARE(%)
5801.5
Y’vjrj
0.0
Con
Y’vjvj
Y’0u
3.3
we have observations at regularly spaced points ti of the unknown function f subjected to noise: (23)
Xi ¼ f ðti Þ þ εi
where the εi are independent Nð0; σ E Þ random noises. After being suitably renormalized at a particular level, Eq. (23) can be transformed into the following form:
(22)
ζt ¼ ζ0t þ ζmax kζ ξt
RPMM
(24)
Zi ¼ ui þ εi
max
where ζ0t is the “clean” sampling data, ζ is the maximum absolute value in the original recorded data, kζ is the specific reduction factor for different variables, and ξt is a discrete random variable subject to a Gaussian distribution with a variance of 0.2. According to the ship model test data published by SIMMAN 2014 & 2020 (ITTC, 2019), the noise is mainly concentrated in the sway while the noise of the remaining responses is smaller, and the deviations of rudder angle can be ignored. This is why the value of kζ is set as 1.0 and 0.3 for the surge velocity and other records, respectively.
where Z ¼ ðZ1 ; …; Zn Þ are observations and εi are independent Nð0; 1Þ variables. The denoising objective is to recover the unknown coefficients u on the basis of the observed data. Here, we assume that the elements ui have independent prior distributions each given by the mixture fprior ðuÞ ¼ ð1
wÞδ0 ðuÞ þ w➇ðuÞ
(25)
where w is the weight used to adjust the near-zero part δ0 ðuÞ and the nonzero part ➇ðuÞ. The nonzero part of the prior, ➇; is assumed to be a fixed unimodal symmetric density. Let ➆ denote the convolution of the density ➇ with the standard normal density Φ. The marginal density of the observations Zi is
4.2. Empirical Bayes denoising EB is widely used in the field of noise filtering, and Jonhstome and Sliverman proposed a fully automatic empirical Bayes thresholding method (Silverman and Johnstone, 2004), which is applicable to various scenarios. The general model for a noisy signal has the following form, where
ð1
wÞΦð➅Þ þ w➆ð➅Þ
(26)
where w can be estimated by the marginal maximum log likelihood to 6
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
Table 3 Comparison of the nondimensional hydrodynamic parameters (1 � 10
5
) for yaw and roll.
N-Coeff.
RPMM
Con
ARE(%)
Semi
ARE(%)
K-Coeff.
RPMM
Con
ARE(%)
Semi
ARE(%)
N’v_
42.3
/
/
/
/
0.0
30
/
/
/
/
K’v_
N’p_
0.2
/
/
/
/
N’r
N’v
300.0
299.7
0.1
298.1
0.6
N’vv
0.6
2.6
336.9
0.4
35.6
N’vjvj
712.9
722.7
1.4
696.1
2.4
N’vjrj
174.7
178.0
1.9
170.0
2.7
N’vrr
36.8
N’rjrj
0.0
54.3
290.0
N’r
291.4 5.4
47.7
5.8 286.7
1.1
/
6.1
/
N’rrr
224.5
220.1
2.0
237.0
5.6
N’rjvj
778.8
780.1
0.2
765.3
1.7
N’rvv
1287.2
1291.6
0.3
1315.8
2.2
N’p
8.0
8.0
0.5
7.6
5.5
5.1
/
N’ppp
0.0
0.2
/
12.8
9.6
25.3
N’pujpuj
0.0
138.2
/
149.4
/
0.7
17.7
1.2
N’vφ
17.8
N’pu
17.9
N’φ
N’vφφ N’φvv
19.4
18.7
4.8
82.9
9309.3
28.3
3047.0
933.9
888.8
4.8
2.6
326.0
0.0
97.9
4.0
6.1
6.5
6.5
N’0u
9.0
6.2
1027.4
128.9
129.0
0.1
128.2
0.5
N’δδ
11.9
11.8
1.2
11.4
101.4
100.1
1.3
N’δv
24.6
24.3
N’δvv
349.1
324.5
N’δδδ
/
/
/
/
0.7
/
/
/
/
K’v
25.0
25.0
0.1
24.9
0.5
/
0.7
0
K’vjvj K’vrr
K’vv
0.0 99.2
98.5
0.7
100.2
1.0
K’vjrj
10.4
10.7
3.2
11.3
8.8
22.2
6.9
69.1
16.9
24.0
K’r
0.8
0.8
0.6
0.9
11.1
K’rrr
0.0
K’rjvj
20.0 41.1
0.3
20.0
0.2
20.2
0.8
3.6
/
2.1
0
40.9
0.5
41.2
0.3
K’rvv
34.6
53.7
55.3
37.5
8.3
3.0
3.0
0.5
3.0
0.6
K’ppp
1.0
1.1
10.1
1.1
11.0
K’p
K’pu
0.0 0.0
0.5
/
0.7
0
K’φ
0.0
21.3
/
22.2
0
K’pujpuj
0.7
/
0.7
0
K’vφ
14.7
14.4
2.1
15.1
103.9
95.6
8.0
85.9
17.3
K’φvv
6.2
20.8
234.9
30.9
398.4
0.8
674.6
1.1
1036.2
K’vφφ
10.0
N’δ
/
/
K’rr
14.4
0.9 0.6
N’0
18.0
14.6
/
1.0
K’p_
84.3
0.5
/
K’r
K’0
K’0u
0.1 1.1
1.1
0.9
1.1
2.5
4.2
K’δ
6.5
6.5
0.1
6.4
1.0
4.0
K’δδ
0.8
0.8
0.2
0.7
8.9
101.4
0.0
K’δδδ
4.1
4.0
3.1
3.7
10.6
1.2
24.9
1.3
K’δv
5.4
5.4
0.0
5.3
2.7
7.0
371.6
6.5
K’δvv
0.1
106.3
0.9
0.9
3.2
N’δu
196.9
198.8
1.0
197.7
0.4
K’δu
8.9
9.0
1.4
9.4
5.4
N’δδu
12.8
14.8
15.3
17.2
34.6
K’δδu
1.3
1.2
4.9
1.7
29.4
8.6
K’δδδu
125.4
N’δδδu
138.2
10.2
136.2
n X
log½ð1
wÞΦð➅i Þ þ w➆ð➅i Þ�
5.3
9.5
8.1
69.7
shown in Table 4, the values of the variance ➉ of β are quite large, indicating that the values of the parameters are scattered around 0 and that the regression results are mainly dependent on observations. The IG prior distribution of σ 2 is determined by the obtained values of A and B. Fig. 3 shows the prior and posterior distributions of σ 2 for surge using conjugate regression and semi-conjugate regression. It should be noted that the prior parameters vary within the same magnitude, which has a slight effect on the results. The mean value of the distribution of the hydrodynamic coefficients is obtained by Eq. (17) for conjugate Bayesian regression, and the semiconjugate regression results are obtained by Gibbs sampling. The sample trajectories of X’u and K’vjvj are shown in Fig. 4; note that each value of the
b: obtain the estimator w
ℓðwÞ ¼
4.8
(27)
i¼1
subject to the constraint on w that the threshold satisfies ➄ðwÞ � pffiffiffiffiffiffiffiffiffiffiffi 2logn.
Since ℓ’ ðwÞ is a monotonic function of w, the root of Eq. (27) can easily be obtained numerically if Eq. (26) is tractable. The main idea of b into the prior and then estimate this method is to insert the value for w the raw value of ui using the value of w thus obtained. We applied the above method to denoise the polluted data. Fig. 2 shows the simulated raw data, polluted data and denoised data from the 20� /20� zigzag manoeuvre. The enlarged figure shows that the EB denoising method has a significant effect. The filtered noise data that have been obtained will be used for the next step of parameter identi fication as training data.
mean of the trajectories is the value for the parameter. Results for the two models in comparison with RPMM test data are listed in Tables 2 and 3, where ARE denotes the absolute relative error, as in Eq. (28). Underscores are used to mark the values with absolute relative error higher than 50%. � � �β βRPMM �� ARE ¼ �� estimate (28) � � 100% βRPMM
4.3. Bayesian regression and validation The denoised data entered contain a total of 400 couples with a sampling frequency of 1 Hz. Since the values of the hydrodynamic co efficient are usually small in the dimensionless model, the priori mean μ is set to 0 to reflect the priori information. Then, the BO algorithm is adopted for the regression model to obtain the hyperparameters of conjugate regression and semi-conjugate regression, and the number of evaluations for BO is set at 40.
As shown in Tables 2 and 3, most results agree well with the experimental data, indicating that the application of the proposed Bayesian approach is indeed robust in the presence of noise. Neverthe less, for some high-order coefficients, such as Y ’pujpuj ; N’vφφ ; and K’φvv , the
difference is not negligible. The discrepancies may be attributed to the sensitivity of these coefficients. The noise in the data cannot be
7
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
Fig. 2. Raw simulation data, polluted data and denoising data by the EB method for the 20� /20� zigzag test.
completely filtered out, so the errors of these high-order variables accumulate by multiplication. It is difficult to judge the merits and usability of the two Bayesian methods only from the results of parameter identification in Tables 2 and 3 Therefore, 10� /10� zigzag motion and 35� turning circle manoeuvring are predicted. As shown in Fig. 5 and Fig. 6, both pre dictions have achieved good agreement with the RPMM experiments, showing the robustness and generalization ability of the proposed Bayesian method. In addition, it can be seen that the accuracy of con jugate regression is higher than that of semi-conjugate regression, especially with respect to roll motion. There are two possible reasons for this phenomenon. First, in contrast with semi-conjugate regression, conjugate regression assumes that the variance and the mean are inde pendent of each other. Considering that the artificial noise added is in dependent of the observed value, the conjugate model is a better fit in the present work. Second, Gibbs sampling, which is used in the process of semi-conjugate regression, is less stable than the analytical solution of conjugate regression. Therefore, the prior hyperparameters used in semi-conjugate regression (ultimately obtained by the BO algorithm) are
less suitable. 5. Conclusions and discussions A new method for hydrodynamic coefficient identification in ship manoeuvring mathematical models based on the Bayesian rule is pre sented and is tested and validated in a nonlinear 4-DOF model with 108 hydrodynamic derivatives. Conjugate and semi-conjugate Bayesian multivariate regression with BO arithmetic are applied for identification with 20� /20� zigzag testing polluted data after noise has been filtered by the EB denoising method. Ship motion simulations of 10� /10� zigzag testing and 35� turning manoeuvring are implemented to verify the proposed method compared to experimental data. Several conclusions can be drawn from the study: 1. After the polluted data are cleaned by the empirical Bayesian method, conjugate regression and semi-conjugate regression models for identification have good generalization performance, but some high-order hydrodynamic coefficients are hard to identify precisely because of their high sensitivity to noise. 2. The proposed Bayesian approach is suitable for nonlinear dimen sionless ship models because the parameters are small. Other prior parameters can be obtained by the BO algorithm, and these hyper parameters have intuitive meanings such as noise and confidence. Furthermore, if we have more prior information of parameters or noise in advance, the Bayesian method will exhibit superiority to the traditional method. 3. Conjugate Bayesian regression shows better prediction results than semi-conjugate Bayesian regression under artificial white noise; our future work will concentrate on full-scale trials with real environ mental disturbances.
Table 4 Prior parameters for conjugate regression and semi-conjugate regression. Con
A
Bð � 103 Þ
➉ð � 1012 Þ
Surge Sway Yaw Roll
8742 3192 5748 5267
156 532 564 741
26.5 13.4 7.8 9.8
Semi
A
Bð � 103 Þ
➉ð � 1012 Þ
Surge Sway Yaw Roll
45637 64815 59672 85241
16873 58742 81644 9438
4.9 8.2 12.7 6.1
8
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
Fig. 3. The prior and posterior distributions of σ 2 for surge.
Fig. 4. Sample trajectories of X’u ; K’vjvj drawn by the Gibbs sampler for the semi-conjugate Bayesian method.
Fig. 5. Comparison of the predicted motions for the 10� /10� zigzag test. 9
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
Fig. 6. Comparison of the predicted motions for the 35� turning circle manoeuvre.
Acknowledgement
ITTC, 2019. Simman 2014 & 2019. https://simman2014.dk/. http://www.simman2019. kr/. Khoshravesh, Mojtaba, Sefidkouhi, Mohammad Ali Gholami, Valipour, Mohammad, 2015. Estimation of reference evapotranspiration using multivariate fractional polynomial, Bayesian regression, and robust regression models in three arid environments. Appl. Water Sci. 7 (.4), 1911–1922. Print. Ljung, L., 1987. System Identification, Theory for the User, Information and System Science Series. Englewood Cliffs. Print. Luo, Weilin, Li, Xinyu, 2017. Measures to diminish the parameter drift in the modeling of ship manoeuvring using system identification. Appl. Ocean Res. 67, 9–20. Print. Moreira, L., Guedes Soares, C., 2003. Dynamic model of manoeuvrability using recursive neural networks. Ocean. Eng. 30 (13), 1669–1697. Print. Moreira, L., Guedes Soares, C., 2012. Recursive neural network model of catamaran manoeuvring. Trans. R. Ins. Naval Archit. Part A: Int. J. Mar. Eng. 154 (PART A3), 121–130. Print. Murphy, Kevin P., 2012. Machine Learning: A Probabilistic Perspective. MIT press. Print. Nagumo, Jin-Ichi, Noda, Atsuhiko, 1967. A learning method for system identification. IEEE Trans. Autom. Control 12 (3), 282–287. Print. Padilla, Arturo, Yuz, Juan I., Herzer, Benjamin, 2014. Continuous-time system identification of the steering dynamics of a ship on a river. Int. J. Control 87 (7), 1387–1405. Print. Park, Trevor, George, Casella, 2008. The Bayesian lasso. J. Am. Stat. Assoc. 103 (482), 681-86. Print. P� erez, Tristan, Blanke, Mogens, 2002. Mathematical Ship Modeling for Control Applications. Print. Perez, T., Fossen, T., 2006. Time-domain models of marine surface vessels based on seakeeping computations. In: Proc. 7th IFAC Conference on Manoeuvring and Control of Marine Craft (MCMC’06).
This study is supported in part by the National Key R&D Program of China (Grant no. SQ2018YFC030044) and in part by the National Nat ural Science Foundation of China (Grant no. U1706230). References Abkowitz, Martin, A., 1980. Measurement of Hydrodynamic Characteristics from Ship Maneuvering Trials by System Identification. Print. Altman, Edward I., 2000. Predicting Financial Distress of Companies: Revisiting the ZScore and Zeta Models. Stern School of Business, New York University, pp. 9–12. Print. Blanke, M., et al., 1986. Identification of a class of nonlinear state-space models using Rpe techniques. In: 1986 25th IEEE Conference on Decision and Control. IEEE. Print. Blanke, Mogens, Jensen, Andreas G., 1997. Dynamic properties of a container vessel with low metacentric height. Trans. Inst. Meas. Control 19 (2), 78–93. Print. Bull, Adam D., 2011. Convergence rates of efficient global optimization algorithms. J. Mach. Learn. Res. 12 (Oct), 2879–2904. Print. Gelfand, Alan E., Smith, Adrian FM., 1990. Sampling-based approaches to calculating marginal densities. J. Am. Stat. Assoc. 85 (410), 398–409. Print. Hoff, Peter D., 2009. A First Course in Bayesian Statistical Methods. Springer. Print. Holzhüter, T., 1990. Robust identification scheme in an adaptive track-controller for ships. In: Adaptive Systems in Control and Signal Processing 1989. Elsevier, pp. 461–466. Print. Iseki, Toshio, Terada, Daisuke, 2002. Bayesian estimation of directional wave spectra for ship guidance system. Int. J. Offshore Polar Eng. 12 (01). Print.
10
Y. Xue et al.
Ocean Engineering 195 (2020) 106612 Tiano, A., et al., 2007. Observer kalman filter identification of an autonomous underwater vehicle. Contr. Eng. Pract. 15 (6), 727–739. Print. Wang, Xue-Gang, et al., 2014. Sensitivity analysis and parametric identification for ship manoeuvring in 4 degrees of freedom. J. Mar. Sci. Technol. 19 (4), 394–405. Print. Wang, Xue-gang, et al., 2015. System identification modeling of ship manoeuvring motion in 4 degrees of freedom based on support vector machines. China Ocean Eng. 29 (4), 519–534. Print. Woltman, Heather, et al., 2012. An introduction to hierarchical linear modeling. Tutorials Quant. Methods Psychol. 8 (1), 52–69. Print. Woo, Joohyun, et al., 2018. Dynamic model identification of unmanned surface vehicles using deep learning network. Appl. Ocean Res. 78, 123–133. Print. Yoon, Hyeon Kyu, Rhee, Key Pyo, 2003. Identification of hydrodynamic coefficients in ship maneuvering equations of motion by estimation-before-modeling technique. Ocean. Eng. 30 (18), 2379–2404. Print. Zhang, Xin-Guang, Zou, Zao-Jian, 2013. Estimation of the hydrodynamic coefficients from captive model test results by using support vector machines. Ocean. Eng. 73, 25–31. Print. Zitzmann, Steffen, et al., 2016. A Bayesian approach for estimating multilevel latent contextual models. Struct. Equ. Model. Multidiscip. J. 23 (5), 661–679. Print.
Petra, Noemi, et al., 2017. A Bayesian approach for parameter estimation with uncertainty for dynamic power systems. IEEE Trans. Power Syst. 32 (4), 2735–2743. Print. Rajesh, G., Bhattacharyya, S.K., 2008. System identification for nonlinear maneuvering of large tankers using artificial neural network. Appl. Ocean Res. 30 (4), 256–263. Print. Robert and Christian, 2014. Machine Learning, a Probabilistic Perspective. Taylor & Francis. Print. Roberts, J Brian, Dunne, Julian F., Debonos, A., 1994. Stochastic estimation methods for non-linear ship roll motion. Probabilistic Eng. Mech. 9 (1–2), 83–93. Print. Silverman, Bernard W., Johnstone, Iain M., 2004. Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences. Ann. Stat. 32 (4), 1594–1649. Print. SNAME, 1950. Nomenclature for Treating the Motion of a Submerged Body through a Fluid. The Society of Naval Architects and Marine Engineers, Technical and Research Bulletin No, pp. 1–5. Print. Son, Kyoungho, Nomoto, Kensaku, 1981. On the coupled motion of steering and rolling of a high speed container ship (1981). J. Soc. Nav. Archit. Jpn. 150, 232–244. Print. Sutulo, Serge, Guedes Soares, C., 2014. An algorithm for offline identification of ship manoeuvring mathematical models from free-running tests. Ocean. Eng. 79, 10–25. Print.
11
Contents lists available at ScienceDirect
Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng
Hydrodynamic parameter identification for ship manoeuvring mathematical models using a Bayesian approach Yifan Xue a, Yanjun Liu a, b, *, Chen Ji a, Gang Xue a a b
Institute of Marine Science and Technology, Shandong University, Qingdao, 266000, China School of Mechanical Engineering, Shandong University, Jinan, 250061, China
A R T I C L E I N F O
A B S T R A C T
Keywords: Ship manoeuvrability Hydrodynamic coefficients System identification Bayesian multilinear regression White noise
A novel identification approach to nonlinear ship manoeuvring models based on Bayes’ rule is presented. The empirical Bayesian method is used to clean simulated polluted responses from a 20� /20� zigzag test. Two effi cient Bayesian models, conjugate regression and semi-conjugate regression, are adopted for hydrodynamic parameter identification. To obtain other prior parameters, the Bayesian optimization (BO) algorithm is intro duced in the Bayesian regression model. By using the identified model and an experimental model based on Blanke’s nonlinear four-degrees-of-freedom (4-DOF) model, 10� /10� zigzag motion and 35� turning circle manoeuvring are performed. A comparison between the predicted results and the test results demonstrates that both Bayesian regression models have good generalization ability for identification, and that the conjugate Bayesian model has greater prediction accuracy.
1. Introduction The highly accurate mathematical modelling of vessels plays a vital role in ship design and operation. For ships with low centre heights, such as container ships, ship manoeuvring motion is usually accompanied by a remarkable amount of roll motion. Son and Nomoto (Son and Nomoto, 1981) studied container ship motion based on experimental tests and established a four-degrees-of-freedom (4-DOF) equation of motion. Blanke.et (Blanke and Jensen, 1997; P� erez and Blanke, 2002) proposed a nonlinear 4-DOF model based on a roll planar motion mechanism (RPMM) facility and validated it via full-scale trials. There are many hydrodynamic parameters in the ship motion model. Added mass in a ship manoeuvring model can be estimated accurately (Perez and Fossen, 2006) by slender body theory. Hence, the key to the prediction accuracy is to obtain the hydrodynamic derivatives, which are regarded as the most uncertain part in the mathematical model. Many researchers have devoted themselves to solving this challenging problem and indicated that system identification (SI) theory in combi nation with full-scale trials or free-running tests (Padilla et al., 2014) is a powerful and practical method. Classic SI techniques are widely applied for ship and underwater vehicle manoeuvring, including least square (LS) estimation (Nagumo and Noda, 1967), improved least squares (Holzhüter, 1990), the
extended Kalman filter (EKF) (Abkowitz and Martin, 1980; Tiano et al., 2007), the recursive prediction error method (RPE) (Blanke et al., 1986) and ridge regression (Yoon and Rhee, 2003). However, these methods have some inherent disadvantages, such as sensitivity to noise or simultaneous drift (Zhang and Zou, 2013; Luo and Li, 2017). In order to overcome these defects, a novel algorithm has been proposed to build a purely “input-output” model on the basis of artificial neural networks (ANN) (Moreira and Guedes Soares, 2003, 2012; Woo et al., 2018). The advantage of ANN is that it is not necessary to prioritize one complex mathematical ship model. However, ANN cannot reflect any physical background. Moreover, the construction of the network structure is blind before full training. In recent years, several robust machine learning methods such as the back propagation neural network (Rajesh and Bhattacharyya, 2008) and support vector machines (SVM) (Wang et al.,2015; Wang et al.,2014) have been proposed and have achieved encouraging results in the application of hydrodynamic parameter identification. However, the hyperparameters introduced in these al gorithms, such as initial weights or penalty factors, have no intuitive physical meaning or logical explanation. The SI methods mentioned above assume that unknown hydrody namic parameters are constants and use limiting relative frequencies to define probability. Roberts.et (Roberts et al., 1994) proposed using Markov process theories to solve the linear regression problem in ship
* Corresponding author. Institute of Marine Science and Technology, shandong University, Qingdao, 266000, China. E-mail address: [email protected] (Y. Liu). https://doi.org/10.1016/j.oceaneng.2019.106612 Received 8 May 2019; Received in revised form 11 September 2019; Accepted 18 October 2019 Available online 25 November 2019 0029-8018/© 2019 Published by Elsevier Ltd.
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
2. Ship mathematical model To describe the motion characteristics of a vessel, as shown in Fig. 1, two coordinate reference frames are adopted, including the Earth-fixed coordinates E X0 Y0 Z0 and the body-fixed coordinates O x0 y0 z0 . According to Son and Nomoto’s nonlinear model (Son and Nomoto, 1981), the uniform equations, including the xG surge, sway, yaw and roll coupled motion, can be transformed as follows: 8 � m u_ vr xG r2 þ ZG rp ¼ FX > > < _ ¼ FY mðv_ þ ur xG r_ ZG pÞ (1) > IZZ r_ þ mxG ðv_ þ urÞ ¼ MN > : Ixx mZG ðv_ þ urÞ ¼ MK ρgrGZ ðφÞ where m is the ship’s mass; u; v; r; and p denote the surge velocity, sway velocity, yaw rate and roll rate, respectively; and zG are the coordinates of the ship’s centre of gravity in the body-fixed coordinate frame; Ixx and Izz are the moments of inertia of the ship about the x0 and y0 axes, respectively; FX and FY are the longitudinal and lateral force compo nents, respectively; MN is the yaw moment about the z0 axis; and MK is the roll moment about the x0 axis. In the last term of Eq. (1), r denotes the ship displacement, g is the gravity constant, ρ is the mass density of the water and Gz ðφÞ is the buoyancy function, which can be approxi mated as: � � 1 Gz ðφÞ ¼ GM þ BMtan2 ðφÞ sinðφÞ (2) 2
Fig. 1. Coordinate frames for a ship.
roll motion from the perspective of probability and statistics. However, the case of one degree of freedom cannot fully prove the applicability of these theories to the overall movement of the ship. Iseki and Terada introduced Bayesian modelling technology to estimate directional wave spectra based on ship motion data and obtained excellent results (Iseki and Terada, 2002). Bayesian methods provide an alternative viewpoint that treats the parameters as random variables and can numerically represent a set of subjective rational beliefs. Bayes’ rule offers a reasonable way to update beliefs in light of new observations. In addi tion to their formal interpretation as a means of induction, Bayesian methods have significant advantages in parameter estimation with good statistical properties, predictions for missing data, forecasting (Hoff, 2009), etc. Recently, Bayesian methods, especially in parameter esti mation, have been applied in multiple fields (Khoshravesh et al., 2015; Petra et al., 2017; Zitzmann et al., 2016) owing to the prevalence of the Bayesian framework and the applicability of modern inference methods based on the Markov chain Monte Carlo (MCMC) method. This article attempts to apply the Bayesian approach to identifying the hydrodynamic parameters in Blanke’s nonlinear 4-DOF model. First, based on the zigzag test simulated by the 4-DOF ship manoeuvring mathematical model (Blanke and Jensen, 1997; P�erez and Blanke, 2002), Gaussian white noise is added to each observed state to simulate the data obtained in the sea trial. An empirical Bayesian (EB) wavelet method (Silverman and Johnstone, 2004) is implemented for removing noise. Then, two efficient Bayesian regression models are utilized to estimate the parameters, and the Bayesian hyperparameter optimization (BO) algorithm is introduced to determine the coefficients of the prior distribution. Finally, in order to validate the proposed method, we use the posterior distribution of the obtained parameters to predict the steering motion of the ship and compare the results with RPMM experimental data.
where GM is the metacentre height, BM is the distance from the centre of buoyancy to the metacentre, and φ is the roll angle. For the convenience of comparison, expanding the force and moment in Eq. (1) with Taylor series (Blanke and Jensen, 1997; P�erez and Blanke, 2002) and rewriting them in a nondimensional matrix using the prime system of SNAME (SNAME, 1950) for substituting model test results yields: 2 3 2 0 3 0 0 2 3 F m X u_ 0 0 0 u _ 6 7 6 10 7 0 0 0 0 0 0 0 0 6 7 6 7 6 7 0 m Y m x Y m Z Y v_ r_ 6 76 v_ 7 6 F2 7 G G 6 74 5 ¼ 6 0 7 (3) 0 0 0 0 0 0 r _ 6 7 6 7 0 m x N I N N v_ r_ p_ G ZZ 4 5 4 M1 5 0 0 0 0 p_ 0 0 m ZG Kv_ Kr_ Ixx K p_ M2 (For notational convenience, the first matrix containing the mass, inertia coefficients and acceleration derivatives in Eq. (3) is labelled J .) The manoeuvring model of Eq. (3) can be discretized using Euler’s method: V
ðiþ1Þ
¼V
ðiÞ
þ
F ⋅h J
(4) T
where V ¼ ½u; v; r; p�T , F ¼ ½F1 ; F2 ; M1 ; M2 � , h is the interval time, and i and i þ 1 are the adjacent sampling time steps. 0
F’1 ¼ Xhydro ⋅AðiÞ þ m’ z’G r’ p’ F’2 ¼ Yhydro ⋅BðiÞ þ m’ u’ r’ M ’1 ¼ Nhydro ⋅CðiÞ þ m’ x’G u’ r’ M ’2 ¼ Khydro ⋅DðiÞ m’ z’G u’ r’
2
0
m’ x’G r’ 2
0
0
m’ v’ r’ (5)
m’ v’ r’
ρgr’ Gz ðφÞ
’
Ocean Engineering 195 (2020) 106612
Y. Xue et al.
where the hydrodynamic derivatives and speed state variables in Eq. (5) are as follows: # " ’ ’ X u ; X uu ; X ’uuu ; X ’vr ; X ’rr ; X ’v ; X ’vv ; X ’vφ ; X ’φ ; X ’φφ ; X ’pp ; X ’ppu ; X ’δ ; X ’δδ ; Xhydro ¼ X ’δu ; X ’δδu ; X ’vδ ; X ’vδδ 1�18
" Yhydro ¼
Nhydro ¼
�
π β; σ2 jy; x ¼ R
� �T BðiÞ ¼ v’ ðiÞ; v’ 2 ðiÞ; v’ ðiÞjv’ ðiÞj; …; δ’2 ðiÞu’a ðiÞ; δ’ 3 ðiÞu’a ðiÞ 1�28
ℓðβ; σ2 jy; xÞπðβÞπ ðσ2 Þ π ðβÞπðσ2 Þρðβ; σ 2 jy; xÞdβd σ2
(10)
where prior distributions are typically denoted πðβÞ and πðσ2 Þ. However, the ultimate goal of regression is to obtain the marginal posterior distribution of β. For different types of prior distributions selected, various Bayesian regression methods such as the hierarchical linear model (Woltman et al., 2012) and the Bayesian lasso model (Park and George, 2008) have been proposed to suit different scenarios. Considering the number of parameters, calculation efficiency and ac curacy, we choose priors so that the corresponding distribution is in the same family of distributions, called the conjugate distributions. Ac cording to the central limit theorem, most of the measured value dis tributions can be approximated by a normal distribution or a Gaussian distribution. A popular choice is the normal-inverse-gamma conjugate model (Robert and Christian, 2014), in which β obeys the multivariate normal distribution (N ) and σ 2 is the inverse gamma (IG) distribution. Eq. (10) can be abbreviated as follows: �Yn � � (11) π β; σ2 jy; x ∝N ðβÞN σ 2 ϕ yt ; xt β; σ 2 t¼1
� �T CðiÞ ¼ v’ ðiÞ; v’ 2 ðiÞ; v’ ðiÞjv’ ðiÞj; v’ ðiÞjr’ ðiÞj; …; δ2 ðiÞu’a ðiÞ; δ3 ðiÞu’a ðiÞ 1�28 � �T DðiÞ ¼ v’ ðiÞ; v’ 2 ðiÞ; v’ ðiÞjv’ ðiÞj; …; δ2 ðiÞu’a ðiÞ; δ’3 ðiÞu’a ðiÞ 1�28 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2 þ v2 ; Unom is the
3. Bayesian regression formulation The object of hydrodynamic parameter identification is to estimate β in the multiple linear regression (MLR) model, as in Eq. (6): (6)
where t ¼ 1, …n denotes time; yt is the observed response; xt is a 1� c row vector of the observed values of c predictors; β is a c� 1 column vector of regression parameters corresponding to the variables that consist of the columns of xt ; and εt is the random disturbance that has a mean of zero and common variance of σ 2 . The Bayesian approach to estimating MLR models treats β and σ 2 as random variables rather than fixed quantities. Generally, the Bayesian analysis process updates the probability distributions of the parameters by incorporating information about the parameters from the observing data, which is a kind of incremental learning. In other words, resolving the Bayesian inverse problem is a posterior probability density function (PDF). Bayes’ theorem gives the posterior PDF as
where ϕðyt ; xt β; σ 2 Þ is the Gaussian probability density with mean xt β and variance σ2 on yt . The regression model is divided into conjugate and semi-conjugate regression(Murphy, 2012) depending on whether parameters and variance are independent. 3.1. Conjugate Bayesian regression When β and σ 2 are independent, the priori distributions of β and σ2 are as follows: � � β�σ2 � Nc μ; σ 2 V (12) σ 2 � IGðA; BÞ
(7)
Eq. (7) can be simplified to posterior ∝ likelihood � prior
1�28
parameter by adjusting the variance of the prior distribution. After specifying the type and parameters of the prior distribution, the joint posterior density can be obtained by substituting the density probability into Eq. (7):
�T � AðiÞ ¼ u’a ðiÞ; u’a 2 ðiÞ; u’a 3 ðiÞ; …; v’ðiÞδ’ðiÞ ; v’ðiÞδ’ 2 ðiÞ 1�18
PðdatajparametersÞPðparametersÞ PðdataÞ
1�28
#
K ’pujpuj ; K ’vφ ; K ’vφφ ; K ’φvv ; K ’0 ; K ’0u ; K ’r ; K ’δ ; K ’δδ ; K ’δδδ ; K ’δv ; K ’δvv ; K ’δu ; K ’δδu ; K ’δδδu
PðparametersjdataÞ ¼
(9)
#
K ’v ; K ’vv ; K ’vjvj ; K ’vjrj ; K ’vrr ; K ’rjrj ; K ’rrr ; K ’rjvj ; K ’rvv ; K ’p ; K ’pp ; K ’ppp ; K ’pu ;
yt ¼ xt β þ εt
� � P yt �xt ; β; σ 2
1�28
N ’v ; N ’vv ; N ’vjvj ; N ’vjrj ; N ’vrr ; N ’r ; N ’rjrj ; N ’rrr ; N ’rjvj ; N ’rvv ; N ’p ; N ’ppp ; N ’pu ;
where the relative surge speed ua ¼ U Unom , U ¼ nominal speed, and δ is the rudder angle.
t¼1
#
Y ’v ; Y ’vv ; Y ’vjvj ; Y ’vjrj ; Y ’vrr ; Y ’r ; Y ’rr ; Y ’rrr ; Y ’rjvj ; Y ’rvv ; Y ’p ; Y ’ppp ; Y ’pu ;
N ’pujpuj ; N ’φ ; N ’vφ ; N ’vφφ ; N ’φvv ; N ’0 ; N ’0u ; N ’δ ; N ’δδ ; N ’δδδ ; N ’δv ; N ’δvv ; N ’δu ; N ’δδu ; N ’δδδu "
Yn
where prior is the parameter distribution that we assume, which allows us to incorporate information about the model before data are imported. Moreover, we can control the confidence in our knowledge about the
Y ’pujpuj ; Y ’φ ; Y ’vφ ; Y ’vφφ ; Y ’φvv ; Y ’0 ; Y ’0u ; Y ’δ ; Y ’δδ ; Y ’δδδ ; Y ’δv ; Y ’δvv ; Y ’δu ; Y ’δδu ; Y ’δδδu "
Khydro ¼
�
ℓ β; σ2 jy; x ¼
where μ is the mean value (c � 1 vector), V is the c � c diagonal matrix in which each element equals the prior variance factor of βj , and IGðA; BÞ denotes the inverse gamma distribution with shape A and scale B. Considering our limited prior information about many parameters, we cannot set different priors for each coefficient, so the diagonal elements in V are simplified to be equal to ➉ in this paper. Since both the prior and
(8)
where likelihood represents the information provided by the sample response about the parameters, and the mathematical model is
3
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
the posterior have the same family of distributions, the conditional posterior distribution of β can be calculated by Eq. (11) and Eq. (13) with σ2 held constant, as shown in Eq. (14): � � � (13) p β�σ 2 ; y; X ∝ pðβÞp yjX; β; σ 2 � � β�σ2 ; y; x �
� Nc
V
1
þ XT X
� 1�
� � XT X b β þ V 1 μ ; σ2 V
1
þ XT X
sampler is an iterative algorithm that constructs a dependent sequence of parameter values whose distribution converges to the target joint posterior distribution. The details of the Gibbs sampler are shown in Algorithm 1.
� 1� (14)
where X is an n � c matrix of training data. The marginal posterior density of σ2 from Eq. (11) is the IG distribution: � � � � 1� � nþc 1 ; B 1 þ SSRðβÞ þ ðb β μÞT V 1 ð b (15) σ 2 ��y; x � IG A þ β μÞ 2 2 where SSRðβÞ is given by n X
SSRðβÞ ¼
yi
β T xi
�2
XβÞT ðy
¼ ðy
(16)
XβÞ
i¼1
To sample from the joint posterior, first draw σ2 from its marginal posterior distribution, given in Eq. (15), and then draw β from its normal conditional posterior distribution, given in Eq. (14), using the simulated value of σ2 . Finally, from the integration of the joint posterior density with respect to σ 2 , the marginal posterior density of β can be calculated as: � � � � � 1� T � � β�y; x � tc V 1 þ X T X X X b β þ V 1 μ ; λ; 2A þ n 2B
1
þ SSRðβÞ þ ðb β
λ¼
h
μÞT V þ X T X
� 1i
1
ðb β
μÞ
Finally, through the two methods above, the forecast of ship motion can be obtained using the mean values of the c -dimensional posterior predictive distribution. 3.3. Hyperparameters in prior distribution A distinctive advantage of the Bayesian approach is that prior in formation and the observed sample are combined in the algorithm. The parameters of the prior distribution, such as the mean and variance, are called hyperparameters in machine learning. In contrast with other commonly used machine learning algorithms, the prior parameters of the Bayesian algorithm have an intuitive meaning: specifying a high prior variance implies that we know very little about the parameter, so a greater weight is given to the information in the data, whereas speci fying a low prior variance implies high confidence in prior information. Based on the feature of Bayesian regression, this paper proposes to use the BO algorithm (Bull, 2011) to determine the parameters (A;B) of the distribution for σ2 , the variance factor of parameters (V) that can reflect the confidence for the mean μ. The algorithm treats the Bayesian
(17)
2A þ n
in which tc ðl; λ; ➈Þ denotes the c-dimensional multivariate t distribution, where l is the location, λ is the scale, and ➈ represents the degrees of freedom. 3.2. Semi-conjugate Bayesian regression When β and σ2 are dependent, the priori distributions of β and σ2 are as follows: � β �σ 2 � Nc ðμ; VÞ (18) σ2 � IGðA; BÞ
Table 1 Physical parameters of the container ship.
As with conjugate regression, a conditional posterior distribution of β and σ 2 can be obtained: � � � � 1� 2 T � � 1� � β�σ2 ; y; x � Nc V 1 þ σ 2 X T X σ X X b β þ V 1 μ ; V 1 þ XT X (19) � � � � n σ 2 ��β; y; x � IG A þ ; B 2
1
1 þ SSRðβÞ 2
� 1�
(20)
Elements
Values
Length (L)
230.66 m
Beam Displacement (r)
32 m 46070 m3
Nominal speed (Unorm )
12.7 m/s
Nominal metacentric height (GM) Transverse metacentre (BM) x coordinate of CG (xG )
0.83 m 9m 0.46 m
Nondim mass (m’ )
750.81 � 10
z coordinate of CG (zG )
However, since β and σ are mutually influential, their posterior distributions are not analytically tractable. Some numerical integration techniques based on the MCMC method have been proposed to solve this problem. In the present work, the Gibbs sampler (Gelfand and Smith, 1990) is applied to approximate the posterior of β and σ2 . The Gibbs 2
Nondim inertia in roll (I’xx )
4
3.54 m 1.30 � 10
5
Nondim inertia in yaw (I’zz )
43.25 � 10
_ Rudder speed (δ)
2.3 deg/s
5
5
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
data and dimensions are given in Table 1. The choice of training manoeuvres is as important as the formulation of a regression model. As Ljung(Ljung, 1987) suggested, the input data for SI should satisfy the persistence of the excitation condition to guar antee convergence and ensure the robustness properties of the identifi cation method. According to his view, the 20� /20� zigzag test and the
regression model as a Gaussian process (GP) and uses the acquisition function (AF) of the expected improvement to find a result that can satisfy the objective function threshold. The basic pseudo-code for Bayesian optimization is shown in Algorithm 2.
Place a Gaussian process prior on Observe
at
While
points according to an initial space-filling experimental design; = do
Update the posterior probability distribution of Let
using available data
be a maximizer of the AF over , where the AF is computed using the current
posterior distribution Observe Increment n end while Return the point evaluated with the largest
Compared with other optimization methods, the BO algorithm works better in cases with fewer hyperparameters and slower operation of the objective model, which coincides with the characteristics of Bayesian regression. Therefore, the BO algorithm is applied in the present work to obtain the hyperparameters of the prior distribution. The training data are pre-processed by Z-score transformation (Altman, 2000) in order to avoid the influence of the high-impact index, and K-fold cross-validation is applied to make full use of the data in this process. K-fold cross-validation randomly divides the set of observations into K groups each time. The first fold is treated as a validation set, and other folds are used for training.
35� turning test reach the appropriate range of extrapolation, unlike the 10� /10� zigzag test, and will therefore provide more informative data for identification. However, the turning circle test is not appropriate for SI because in the motion of the circle after the rudder is turned, the observation values are very similar or even constant, which will cause the vectors in the matrix to be highly correlated. Therefore, a 20� /20� standard zigzag manoeuvre was simulated in MATLAB with 20000-point output discretization by using the hydrodynamic parameters from the RPMM facility at the Danish Maritime Institute (scale ratio 28.75), as given in Table 2 and Table 3. Note that some hydrodynamic derivatives were adjusted to match the full-scale trial results(Blanke and Jensen, 1997). The sampling interval is 0.02 s, and the total time is 400 s. The raw training data couples consist of the following:
4. Numerical example, results and discussion 4.1. Construction of the training model A container vessel was selected as the object model, and the main
input variables : ½AðiÞ;BðiÞ;CðiÞ;DðiÞ� output observed response : 2 m’
� u’ ði þ 1Þ u’ðiÞ X ’u_ L UðiÞh
2
m’ v’ ðiÞr’ ðiÞ
3
m’ x’G r’ ðiÞ þ m’ z’G r’ ðiÞp’ðiÞ
6 7 6 7 6 7 � � ’ ’ ’ ’ ’ ’ 6 7 p ði þ 1Þ p ðiÞ 6 m’ Y ’ �L v ði þ 1Þ v ðiÞ þ m’ x’ Y ’ �L r ði þ 1Þ r ðiÞ 7 ’ ’ ’ L m þ m’u’ðiÞr’ðiÞ z þ Y 6 7 v_ r_ p_ G G UðiÞh UðiÞh UðiÞh 6 7 6 7 6 7 6 7 6 7 � � ’ ’ ’ ’ ’ ’ � � 6 ’ ’ 7 v ði þ 1Þ v ðiÞ r ði þ 1Þ r ðiÞ p ði þ 1Þ p ðiÞ ’ ’ ’ ’ ’ ’ 6 m xG N v_ L 7 þ I þ N þ m L L N x u’ðiÞr’ðiÞ r_ p_ zz G 6 7 UðiÞh UðiÞh UðiÞh 6 7 6 7 6 7 6 7 � p’ ði þ 1Þ p’ ðiÞ 6 7 � v’ ði þ 1Þ v’ ðiÞ r’ ði þ 1Þ r’ ðiÞ � ’ ’ ’ ’ ’ ’ ’ 6 7 m’zG Kv_ L Kr_L þ I xx K p_ L m zG u ðiÞr ðiÞ 6 7 UðiÞh UðiÞh UðiÞh 4 5 þ ρgr0 G’z ðφÞ
5
(21)
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
Table 2 Comparison of the nondimensional hydrodynamic parameters (1 � 10 X-Coeff. X’u_
RPMM 124:4
5
) for surge and sway.
Con
ARE(%)
Semi
ARE(%)
Y-Coeff.
/
/
/
/
Y’v_ Y’r
226.2
226.4
0.05
226.3
0.1
X’uu
64.5
67.5
0.64
64.8
4.6
137.2
147.9
0.86
138.0
7.8
X’vr
24.0
23.9
3.61
24.9
0.4
X’u
X’uuu X’rr
4.4
4.4
14.38
3.8
41.0
24.0
41.51
24.0
41.5
1.0
0.9
34.22
1.4
10.9
X’vφ
108.1
108.5
0.06
108.2
0.3
X’φ
5.9
5.9
0.18
5.9
0.2
X’φφ
42.2
42.6
0.08
42.3
0.9
X’v
X’vv
X’pp
X’ppu
7.2
7.1
0.15
7.2
0.8
3.9
3.9
1.98
3.9
0.0
X’δ
1.4
1.4
0.10
1.4
1.5
X’δδ
116.8
116.8
0.01
116.8
0.0
X’δu ;
17.2
17.0
0.03
17.2
0.9
X’δδu
224.9
224.6
0.12
225.1
0.1
X’vδ
124.5
124.4
0.18
124.7
0.1
X’vδδ
341.0
341.2
0.02
341.0
878:0 48.1
Y’p_
23.3
Y’vv
98.6
725.0
Y’v
0.8
Semi
ARE(%)
/
/
/
/
/
/
/
/
/
/
/
723.3 91.0
0.2 7.7
/ 717.5
73.4
1.0 25.5
5820.0
0.3
5766.6
0.6
1192.7
1200.9
0.7
1196.8
0.3
Y’vrr
1107.9
923.9
16.6
1232.4
11.2
Y’r
Y’rr
118.2
115.4
2.4
125.6
0.0
13.1
/
11.9
6.3 /
Y’rrr
158.0
116.1
26.5
211.5
33.9
Y’rjvj
409.4
414.1
1.2
385.0
6.0
Y’rvv
994.6
787.5
20.8
1081.8
8.8
Y’p
3.4
3.2
5.0
2.1
38.9
Y’ppp
9.3
11.3
21.3
26.3
182.8
Y’pu
23.6
22.6
4.3
78.5
249.5
Y’φ
37.7
37.0
1.8
37.1
1.5
Y’vφ
144.9
150.3
3.7
164.5
13.5
Y’vφφ
2459.3
2327.1
5.4
2405.4
2.2
Y’φvv
177.2
448.9
153.3
288.8
63.0
Y’0
4.7
5.6
19.2
2.4
48.3
5.3
6.4
21.7
5.9
11.3
Y’δ
248.1
248.2
0.0
248.8
0.3
Y’δδ
13.4
13.9
3.8
13.2
1.2
52.5
Y’pujpuj
33.7 632.1
43.0 1104.1
Y’δδδ
193.0
197.7
2.4
186.4
3.4
Y’δv
100.0
98.5
1.5
97.0
3.0
189.2
Y’δvv
233.8
23.6
132.9
29.7
Y’δu
379.4
374.2
1.4
381.0
0.4
Y’δδu
55.6
47.7
14.3
48.3
13.1
232.3
182.8
21.3
240.1
Y’δδδu
It is almost meaningless to use the raw data above for identification directly, as the real experimental measurements will inevitably be affected by various disturbances and produce a great deal of highfrequency responses such as ship oscillatory motion or mechanical noise. Following Sutulo and Soares’s artificially polluted data con struction method (Sutulo and Guedes Soares, 2014), Gaussian white noise is added as follows:
ARE(%)
5801.5
Y’vjrj
0.0
Con
Y’vjvj
Y’0u
3.3
we have observations at regularly spaced points ti of the unknown function f subjected to noise: (23)
Xi ¼ f ðti Þ þ εi
where the εi are independent Nð0; σ E Þ random noises. After being suitably renormalized at a particular level, Eq. (23) can be transformed into the following form:
(22)
ζt ¼ ζ0t þ ζmax kζ ξt
RPMM
(24)
Zi ¼ ui þ εi
max
where ζ0t is the “clean” sampling data, ζ is the maximum absolute value in the original recorded data, kζ is the specific reduction factor for different variables, and ξt is a discrete random variable subject to a Gaussian distribution with a variance of 0.2. According to the ship model test data published by SIMMAN 2014 & 2020 (ITTC, 2019), the noise is mainly concentrated in the sway while the noise of the remaining responses is smaller, and the deviations of rudder angle can be ignored. This is why the value of kζ is set as 1.0 and 0.3 for the surge velocity and other records, respectively.
where Z ¼ ðZ1 ; …; Zn Þ are observations and εi are independent Nð0; 1Þ variables. The denoising objective is to recover the unknown coefficients u on the basis of the observed data. Here, we assume that the elements ui have independent prior distributions each given by the mixture fprior ðuÞ ¼ ð1
wÞδ0 ðuÞ þ w➇ðuÞ
(25)
where w is the weight used to adjust the near-zero part δ0 ðuÞ and the nonzero part ➇ðuÞ. The nonzero part of the prior, ➇; is assumed to be a fixed unimodal symmetric density. Let ➆ denote the convolution of the density ➇ with the standard normal density Φ. The marginal density of the observations Zi is
4.2. Empirical Bayes denoising EB is widely used in the field of noise filtering, and Jonhstome and Sliverman proposed a fully automatic empirical Bayes thresholding method (Silverman and Johnstone, 2004), which is applicable to various scenarios. The general model for a noisy signal has the following form, where
ð1
wÞΦð➅Þ þ w➆ð➅Þ
(26)
where w can be estimated by the marginal maximum log likelihood to 6
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
Table 3 Comparison of the nondimensional hydrodynamic parameters (1 � 10
5
) for yaw and roll.
N-Coeff.
RPMM
Con
ARE(%)
Semi
ARE(%)
K-Coeff.
RPMM
Con
ARE(%)
Semi
ARE(%)
N’v_
42.3
/
/
/
/
0.0
30
/
/
/
/
K’v_
N’p_
0.2
/
/
/
/
N’r
N’v
300.0
299.7
0.1
298.1
0.6
N’vv
0.6
2.6
336.9
0.4
35.6
N’vjvj
712.9
722.7
1.4
696.1
2.4
N’vjrj
174.7
178.0
1.9
170.0
2.7
N’vrr
36.8
N’rjrj
0.0
54.3
290.0
N’r
291.4 5.4
47.7
5.8 286.7
1.1
/
6.1
/
N’rrr
224.5
220.1
2.0
237.0
5.6
N’rjvj
778.8
780.1
0.2
765.3
1.7
N’rvv
1287.2
1291.6
0.3
1315.8
2.2
N’p
8.0
8.0
0.5
7.6
5.5
5.1
/
N’ppp
0.0
0.2
/
12.8
9.6
25.3
N’pujpuj
0.0
138.2
/
149.4
/
0.7
17.7
1.2
N’vφ
17.8
N’pu
17.9
N’φ
N’vφφ N’φvv
19.4
18.7
4.8
82.9
9309.3
28.3
3047.0
933.9
888.8
4.8
2.6
326.0
0.0
97.9
4.0
6.1
6.5
6.5
N’0u
9.0
6.2
1027.4
128.9
129.0
0.1
128.2
0.5
N’δδ
11.9
11.8
1.2
11.4
101.4
100.1
1.3
N’δv
24.6
24.3
N’δvv
349.1
324.5
N’δδδ
/
/
/
/
0.7
/
/
/
/
K’v
25.0
25.0
0.1
24.9
0.5
/
0.7
0
K’vjvj K’vrr
K’vv
0.0 99.2
98.5
0.7
100.2
1.0
K’vjrj
10.4
10.7
3.2
11.3
8.8
22.2
6.9
69.1
16.9
24.0
K’r
0.8
0.8
0.6
0.9
11.1
K’rrr
0.0
K’rjvj
20.0 41.1
0.3
20.0
0.2
20.2
0.8
3.6
/
2.1
0
40.9
0.5
41.2
0.3
K’rvv
34.6
53.7
55.3
37.5
8.3
3.0
3.0
0.5
3.0
0.6
K’ppp
1.0
1.1
10.1
1.1
11.0
K’p
K’pu
0.0 0.0
0.5
/
0.7
0
K’φ
0.0
21.3
/
22.2
0
K’pujpuj
0.7
/
0.7
0
K’vφ
14.7
14.4
2.1
15.1
103.9
95.6
8.0
85.9
17.3
K’φvv
6.2
20.8
234.9
30.9
398.4
0.8
674.6
1.1
1036.2
K’vφφ
10.0
N’δ
/
/
K’rr
14.4
0.9 0.6
N’0
18.0
14.6
/
1.0
K’p_
84.3
0.5
/
K’r
K’0
K’0u
0.1 1.1
1.1
0.9
1.1
2.5
4.2
K’δ
6.5
6.5
0.1
6.4
1.0
4.0
K’δδ
0.8
0.8
0.2
0.7
8.9
101.4
0.0
K’δδδ
4.1
4.0
3.1
3.7
10.6
1.2
24.9
1.3
K’δv
5.4
5.4
0.0
5.3
2.7
7.0
371.6
6.5
K’δvv
0.1
106.3
0.9
0.9
3.2
N’δu
196.9
198.8
1.0
197.7
0.4
K’δu
8.9
9.0
1.4
9.4
5.4
N’δδu
12.8
14.8
15.3
17.2
34.6
K’δδu
1.3
1.2
4.9
1.7
29.4
8.6
K’δδδu
125.4
N’δδδu
138.2
10.2
136.2
n X
log½ð1
wÞΦð➅i Þ þ w➆ð➅i Þ�
5.3
9.5
8.1
69.7
shown in Table 4, the values of the variance ➉ of β are quite large, indicating that the values of the parameters are scattered around 0 and that the regression results are mainly dependent on observations. The IG prior distribution of σ 2 is determined by the obtained values of A and B. Fig. 3 shows the prior and posterior distributions of σ 2 for surge using conjugate regression and semi-conjugate regression. It should be noted that the prior parameters vary within the same magnitude, which has a slight effect on the results. The mean value of the distribution of the hydrodynamic coefficients is obtained by Eq. (17) for conjugate Bayesian regression, and the semiconjugate regression results are obtained by Gibbs sampling. The sample trajectories of X’u and K’vjvj are shown in Fig. 4; note that each value of the
b: obtain the estimator w
ℓðwÞ ¼
4.8
(27)
i¼1
subject to the constraint on w that the threshold satisfies ➄ðwÞ � pffiffiffiffiffiffiffiffiffiffiffi 2logn.
Since ℓ’ ðwÞ is a monotonic function of w, the root of Eq. (27) can easily be obtained numerically if Eq. (26) is tractable. The main idea of b into the prior and then estimate this method is to insert the value for w the raw value of ui using the value of w thus obtained. We applied the above method to denoise the polluted data. Fig. 2 shows the simulated raw data, polluted data and denoised data from the 20� /20� zigzag manoeuvre. The enlarged figure shows that the EB denoising method has a significant effect. The filtered noise data that have been obtained will be used for the next step of parameter identi fication as training data.
mean of the trajectories is the value for the parameter. Results for the two models in comparison with RPMM test data are listed in Tables 2 and 3, where ARE denotes the absolute relative error, as in Eq. (28). Underscores are used to mark the values with absolute relative error higher than 50%. � � �β βRPMM �� ARE ¼ �� estimate (28) � � 100% βRPMM
4.3. Bayesian regression and validation The denoised data entered contain a total of 400 couples with a sampling frequency of 1 Hz. Since the values of the hydrodynamic co efficient are usually small in the dimensionless model, the priori mean μ is set to 0 to reflect the priori information. Then, the BO algorithm is adopted for the regression model to obtain the hyperparameters of conjugate regression and semi-conjugate regression, and the number of evaluations for BO is set at 40.
As shown in Tables 2 and 3, most results agree well with the experimental data, indicating that the application of the proposed Bayesian approach is indeed robust in the presence of noise. Neverthe less, for some high-order coefficients, such as Y ’pujpuj ; N’vφφ ; and K’φvv , the
difference is not negligible. The discrepancies may be attributed to the sensitivity of these coefficients. The noise in the data cannot be
7
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
Fig. 2. Raw simulation data, polluted data and denoising data by the EB method for the 20� /20� zigzag test.
completely filtered out, so the errors of these high-order variables accumulate by multiplication. It is difficult to judge the merits and usability of the two Bayesian methods only from the results of parameter identification in Tables 2 and 3 Therefore, 10� /10� zigzag motion and 35� turning circle manoeuvring are predicted. As shown in Fig. 5 and Fig. 6, both pre dictions have achieved good agreement with the RPMM experiments, showing the robustness and generalization ability of the proposed Bayesian method. In addition, it can be seen that the accuracy of con jugate regression is higher than that of semi-conjugate regression, especially with respect to roll motion. There are two possible reasons for this phenomenon. First, in contrast with semi-conjugate regression, conjugate regression assumes that the variance and the mean are inde pendent of each other. Considering that the artificial noise added is in dependent of the observed value, the conjugate model is a better fit in the present work. Second, Gibbs sampling, which is used in the process of semi-conjugate regression, is less stable than the analytical solution of conjugate regression. Therefore, the prior hyperparameters used in semi-conjugate regression (ultimately obtained by the BO algorithm) are
less suitable. 5. Conclusions and discussions A new method for hydrodynamic coefficient identification in ship manoeuvring mathematical models based on the Bayesian rule is pre sented and is tested and validated in a nonlinear 4-DOF model with 108 hydrodynamic derivatives. Conjugate and semi-conjugate Bayesian multivariate regression with BO arithmetic are applied for identification with 20� /20� zigzag testing polluted data after noise has been filtered by the EB denoising method. Ship motion simulations of 10� /10� zigzag testing and 35� turning manoeuvring are implemented to verify the proposed method compared to experimental data. Several conclusions can be drawn from the study: 1. After the polluted data are cleaned by the empirical Bayesian method, conjugate regression and semi-conjugate regression models for identification have good generalization performance, but some high-order hydrodynamic coefficients are hard to identify precisely because of their high sensitivity to noise. 2. The proposed Bayesian approach is suitable for nonlinear dimen sionless ship models because the parameters are small. Other prior parameters can be obtained by the BO algorithm, and these hyper parameters have intuitive meanings such as noise and confidence. Furthermore, if we have more prior information of parameters or noise in advance, the Bayesian method will exhibit superiority to the traditional method. 3. Conjugate Bayesian regression shows better prediction results than semi-conjugate Bayesian regression under artificial white noise; our future work will concentrate on full-scale trials with real environ mental disturbances.
Table 4 Prior parameters for conjugate regression and semi-conjugate regression. Con
A
Bð � 103 Þ
➉ð � 1012 Þ
Surge Sway Yaw Roll
8742 3192 5748 5267
156 532 564 741
26.5 13.4 7.8 9.8
Semi
A
Bð � 103 Þ
➉ð � 1012 Þ
Surge Sway Yaw Roll
45637 64815 59672 85241
16873 58742 81644 9438
4.9 8.2 12.7 6.1
8
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
Fig. 3. The prior and posterior distributions of σ 2 for surge.
Fig. 4. Sample trajectories of X’u ; K’vjvj drawn by the Gibbs sampler for the semi-conjugate Bayesian method.
Fig. 5. Comparison of the predicted motions for the 10� /10� zigzag test. 9
Y. Xue et al.
Ocean Engineering 195 (2020) 106612
Fig. 6. Comparison of the predicted motions for the 35� turning circle manoeuvre.
Acknowledgement
ITTC, 2019. Simman 2014 & 2019. https://simman2014.dk/. http://www.simman2019. kr/. Khoshravesh, Mojtaba, Sefidkouhi, Mohammad Ali Gholami, Valipour, Mohammad, 2015. Estimation of reference evapotranspiration using multivariate fractional polynomial, Bayesian regression, and robust regression models in three arid environments. Appl. Water Sci. 7 (.4), 1911–1922. Print. Ljung, L., 1987. System Identification, Theory for the User, Information and System Science Series. Englewood Cliffs. Print. Luo, Weilin, Li, Xinyu, 2017. Measures to diminish the parameter drift in the modeling of ship manoeuvring using system identification. Appl. Ocean Res. 67, 9–20. Print. Moreira, L., Guedes Soares, C., 2003. Dynamic model of manoeuvrability using recursive neural networks. Ocean. Eng. 30 (13), 1669–1697. Print. Moreira, L., Guedes Soares, C., 2012. Recursive neural network model of catamaran manoeuvring. Trans. R. Ins. Naval Archit. Part A: Int. J. Mar. Eng. 154 (PART A3), 121–130. Print. Murphy, Kevin P., 2012. Machine Learning: A Probabilistic Perspective. MIT press. Print. Nagumo, Jin-Ichi, Noda, Atsuhiko, 1967. A learning method for system identification. IEEE Trans. Autom. Control 12 (3), 282–287. Print. Padilla, Arturo, Yuz, Juan I., Herzer, Benjamin, 2014. Continuous-time system identification of the steering dynamics of a ship on a river. Int. J. Control 87 (7), 1387–1405. Print. Park, Trevor, George, Casella, 2008. The Bayesian lasso. J. Am. Stat. Assoc. 103 (482), 681-86. Print. P� erez, Tristan, Blanke, Mogens, 2002. Mathematical Ship Modeling for Control Applications. Print. Perez, T., Fossen, T., 2006. Time-domain models of marine surface vessels based on seakeeping computations. In: Proc. 7th IFAC Conference on Manoeuvring and Control of Marine Craft (MCMC’06).
This study is supported in part by the National Key R&D Program of China (Grant no. SQ2018YFC030044) and in part by the National Nat ural Science Foundation of China (Grant no. U1706230). References Abkowitz, Martin, A., 1980. Measurement of Hydrodynamic Characteristics from Ship Maneuvering Trials by System Identification. Print. Altman, Edward I., 2000. Predicting Financial Distress of Companies: Revisiting the ZScore and Zeta Models. Stern School of Business, New York University, pp. 9–12. Print. Blanke, M., et al., 1986. Identification of a class of nonlinear state-space models using Rpe techniques. In: 1986 25th IEEE Conference on Decision and Control. IEEE. Print. Blanke, Mogens, Jensen, Andreas G., 1997. Dynamic properties of a container vessel with low metacentric height. Trans. Inst. Meas. Control 19 (2), 78–93. Print. Bull, Adam D., 2011. Convergence rates of efficient global optimization algorithms. J. Mach. Learn. Res. 12 (Oct), 2879–2904. Print. Gelfand, Alan E., Smith, Adrian FM., 1990. Sampling-based approaches to calculating marginal densities. J. Am. Stat. Assoc. 85 (410), 398–409. Print. Hoff, Peter D., 2009. A First Course in Bayesian Statistical Methods. Springer. Print. Holzhüter, T., 1990. Robust identification scheme in an adaptive track-controller for ships. In: Adaptive Systems in Control and Signal Processing 1989. Elsevier, pp. 461–466. Print. Iseki, Toshio, Terada, Daisuke, 2002. Bayesian estimation of directional wave spectra for ship guidance system. Int. J. Offshore Polar Eng. 12 (01). Print.
10
Y. Xue et al.
Ocean Engineering 195 (2020) 106612 Tiano, A., et al., 2007. Observer kalman filter identification of an autonomous underwater vehicle. Contr. Eng. Pract. 15 (6), 727–739. Print. Wang, Xue-Gang, et al., 2014. Sensitivity analysis and parametric identification for ship manoeuvring in 4 degrees of freedom. J. Mar. Sci. Technol. 19 (4), 394–405. Print. Wang, Xue-gang, et al., 2015. System identification modeling of ship manoeuvring motion in 4 degrees of freedom based on support vector machines. China Ocean Eng. 29 (4), 519–534. Print. Woltman, Heather, et al., 2012. An introduction to hierarchical linear modeling. Tutorials Quant. Methods Psychol. 8 (1), 52–69. Print. Woo, Joohyun, et al., 2018. Dynamic model identification of unmanned surface vehicles using deep learning network. Appl. Ocean Res. 78, 123–133. Print. Yoon, Hyeon Kyu, Rhee, Key Pyo, 2003. Identification of hydrodynamic coefficients in ship maneuvering equations of motion by estimation-before-modeling technique. Ocean. Eng. 30 (18), 2379–2404. Print. Zhang, Xin-Guang, Zou, Zao-Jian, 2013. Estimation of the hydrodynamic coefficients from captive model test results by using support vector machines. Ocean. Eng. 73, 25–31. Print. Zitzmann, Steffen, et al., 2016. A Bayesian approach for estimating multilevel latent contextual models. Struct. Equ. Model. Multidiscip. J. 23 (5), 661–679. Print.
Petra, Noemi, et al., 2017. A Bayesian approach for parameter estimation with uncertainty for dynamic power systems. IEEE Trans. Power Syst. 32 (4), 2735–2743. Print. Rajesh, G., Bhattacharyya, S.K., 2008. System identification for nonlinear maneuvering of large tankers using artificial neural network. Appl. Ocean Res. 30 (4), 256–263. Print. Robert and Christian, 2014. Machine Learning, a Probabilistic Perspective. Taylor & Francis. Print. Roberts, J Brian, Dunne, Julian F., Debonos, A., 1994. Stochastic estimation methods for non-linear ship roll motion. Probabilistic Eng. Mech. 9 (1–2), 83–93. Print. Silverman, Bernard W., Johnstone, Iain M., 2004. Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences. Ann. Stat. 32 (4), 1594–1649. Print. SNAME, 1950. Nomenclature for Treating the Motion of a Submerged Body through a Fluid. The Society of Naval Architects and Marine Engineers, Technical and Research Bulletin No, pp. 1–5. Print. Son, Kyoungho, Nomoto, Kensaku, 1981. On the coupled motion of steering and rolling of a high speed container ship (1981). J. Soc. Nav. Archit. Jpn. 150, 232–244. Print. Sutulo, Serge, Guedes Soares, C., 2014. An algorithm for offline identification of ship manoeuvring mathematical models from free-running tests. Ocean. Eng. 79, 10–25. Print.
11