2.1. Sway-yaw dynamics

The mathematical model for ship dynamics is introduced for a naval vessel. The dynamic equations of motion and corresponding forces in the body fixed frame can be represented as:

(1)$$\eqalign{( m-Y_{\dot {v}} ) \dot {v}+( mx_G -Y_{\dot {r}} ) \dot {r} &= Y_{\vert U \vert v} \vert U \vert v+Y_{Ur} Ur+Y_{v\vert v \vert } v\vert v \vert \cr & \quad +Y_{v\vert r \vert } v\vert r \vert +Y_{r\vert v \vert } r\vert v \vert -mUr+Y_c} $$

(2)$$\eqalign{( mx_G -N_{\dot {v}} ) \dot {v}+( I_{zz} -N_{\dot {r}} ) \dot {r} &= N_{\vert U \vert v} \vert U \vert v+N_{\vert U \vert r} \vert U \vert r+N_{r\vert r \vert } r\vert r \vert \cr & \quad +N_{r\vert v \vert } r\vert v \vert -mx_G Ur\hbox{+}N_c} $$

(3)$${\dot {\psi }} = r $$

where m is the ship mass and I _zz is inertia. x _G are the coordinates of the centre of gravity with respect to the body fixed frame. Sway velocity is represented by v. Yaw and its angular velocity are denoted by ψ and r, respectively. Y and N are external forces with respect to sway and yaw. Subscript c denotes the force or moment produced by the control surface.

The rudder is the device used for heading control here. The rudder induced forces and moments can be expressed in the following:

(4)$$Y_c =\displaystyle{{1}\over{2}}\rho A_R C_L U^2 $$

(5)$$N_c =-\displaystyle{{1}\over{2}}\rho A_R C_L U^2LCG $$

where ρ is the water density, A _R is the area of the rudder, C _L is the lift coefficient which varies with the effective angle of attack and LCG is the distance from the centre of gravity to the rudder stock.

2.2. Steering model

The ship steering system in this study is an underactuated system where the sway motion cannot be directly controlled. The study aims to achieve ship course keeping, so only yaw motion is considered. Combining the rudder action N _c = N _δδ and the wave disturbance, the steering dynamics can be treated as:

(6)$$( I_{zz} -N_{\dot {r}} ) \dot {r}=N_{\vert U \vert r} \vert U \vert r+N_{r\vert r \vert } r\vert r \vert -mx_G Ur+N_{\delta} \delta +d_w $$

where δ is rudder angle and d _w is the yaw moment caused by waves.

In practice, the dynamics of the ship and actuator are changed by the time varying speed which causes parameter uncertainties. The upper bounds of disturbance and parameter uncertainties are often difficult to find. The uncertainties in a system are assumed to meet the matching conditions. In this section, a modified control model is presented, considering the time varying speed U(t) = U ₀ + Δ U due to the speed loss, where U ₀ is the nominal (designed) speed and Δ U is the time varying speed loss which is unknown but bounded. The steering model is written as:

(7)$$\dot {r}=a_1 ( U) r+a_2 r\vert r \vert +a_3 \delta +a_4 d_w $$

where $a_{1} ( U) =\displaystyle{{N_{\vert U \vert r} \vert U \vert -mx_{G} Ur}\over{I_{zz} -N_{\dot {r}}}}$, $a_{2} =\displaystyle{{N_{r\vert r \vert }}\over{I_{zz} -N_{\dot {r}}}}$, $a_{3} =\displaystyle{{N_{\delta}}\over{I_{zz} -N_{\dot {r}}}}$, $a_{4} =\displaystyle{{1}\over{I_{zz}-N_{\dot {r}}}}$

Considering the dynamic uncertainties caused by the speed loss, the steering model is revised as:

(8)$$\dot {r}=a_1 ( U_0 ) r+a_2 r\vert r \vert +a_3 \delta +d $$

where d = d ₁ + d ₂ = a ₁ (Δ U) r + a ₄ d _w is the lumped uncertainties, d ₁ is the internal disturbance and d ₂ is the external disturbance.

Assumption 1: The position and rate measurements of the yaw motion of the vehicle are available for feedback. The gyroscopic compass measures ψ and the yaw rate gyro measures r.

Assumption 2: The disturbance signal satisfies |d(t)| ≤ d _max and d _max is an unknown positive constant.

Remark 1: The sway motion cannot be directly controlled, but it is influenced by the yaw motion control commanded rudder angle.

2.3. Wave model

Complex sea states are considered to be the superposition of an infinite number of monochromatic waves, distributed in all directions. In order to simulate random waves, an International Towing Tank Conference (ITTC) long-crest wave spectrum is adopted to recreate a fully developed sea environment. The equation of wave Power Spectral Density (PSD) is given as follows:

(9)$$S( {\omega_i } ) =\displaystyle{{173H_{1/3}}\over{T^4\omega_i^5}}\exp \left({-\displaystyle{{691}\over{T^4\omega_i^4}}} \right)$$

where H _1/3 is the significant wave height, T is the wave period and ω _i is the wave frequency of the i-th regular wave component.

In this paper, 60 regular wave components are used to form an irregular wave. The amplitude of each regular wave component ζ _i and the resultant wave ζ can be obtained by the following equations:

(10)$$\zeta _i = \sqrt {2S( {\omega_i}) \Delta \omega} $$

(11)$$\zeta = \sum\limits_{i=1}^{60} {\zeta_i \cos ( {\omega_{i}t+\varepsilon_i}) } $$

where ε _i is the random phase angle of the i − th regular wave, which ranges from 0 to 2π. In this calculation, the resultant wave ζ is used to calculate the external forces. This calculation is performed using code prepared in MATLAB, along with the calculation of the wave excitation forces and moments.

Generally, when a vessel is sailing with a constant speed, the wave frequency is modified as the encounter frequency and it is expressed as follows:

(12)$$\omega _{ei} =\omega _i -\displaystyle{{\omega _i^2 }\over{g}}\cos \beta $$

where β is the encounter angle, g is the gravitational acceleration.

In order to simulate the time series of external disturbances accurately, the time domain model of the yaw moment caused by waves is derived with the strip theory (Bhattacharyya, Reference Bhattacharyya1978) and it is given as follows:

(13)$$d_w =\sum\limits_{i=1}^n {[ d_{i1} \zeta _i \cos ( \omega _{ei} t+\varepsilon_i ) +d_{i2} \zeta _i \sin ( \omega _{ei} t+\varepsilon_i ) ] } $$

where $d_{i1} =2\rho g\sum\limits_{j=1}^{N} {\lcub \exp ( -{k_{i} T_{j}}/2) T_{j} x_{j} \sin [ k_{i} B_{j} \sin ( \beta/2) ] \cos ( k_{i} x_{j} \cos \beta ) \Delta x\rcub } $, ρ is the water density, and $d_{i2} =2\rho g\sum\limits_{j=1}^{N} {\lcub \exp ( {-k_{i} T_{j}} /2) T_{j} x_{i} \sin [ k_{i} B_{j} \sin ( \beta/2) ] \sin ( k_{i} x_{i} \cos \beta ) \Delta x\rcub }$. ζ _i and k _i are the wave amplitude and the wave number of the wave component i. N and Δ x are the section number of the ship and the length of each section. T _j, B _j and x _j are draft, breadth and the coordinate point of the section j, respectively.

2.4. Added resistance

Generally, the design of propulsive power of a ship is based on still water resistance, which will be a constant when sailing on the open sea. The forward speed will be reduced due to the added resistance which is independent of the calm water resistance. In the following we consider the added resistance caused by ship motions both in waves and in still water.

In oblique waves, the force and moment which should be applied externally can be separated according to energy considerations. The force along the x-axis (added resistance) expends energy, whereas the one along the y-axis (drift force) does not expend energy. These forces are calculated by the method proposed by Loukakis and Sclavounos (Reference Loukakis and Sclavounos1978). As shown in Figure 1, the ship is moving forward along a fixed direction. The waves arrive at an encounter angle β with a speed of propagation c. The horizontal force R _T can be resolved into two components, the added resistance R _x and the drift force R _y.

Figure 1.

Definition of added resistance and drift force.

According to the extended radiated energy theory (Loukakis and Sclavounos, Reference Loukakis and Sclavounos1978), the work of the added resistance is assumed to equate to the energy contained in the dumping waves radiated away from the ship during each wave period. The work of the horizontal force to the energy per encounter period radiated away from from the ship which moves with five degrees of freedom can be equated with strip theory approximations. Therefore, the work transformed to the fluid can therefore be expressed as:

(14)$$P=( -R_T ) ( -c-U\cos \beta ) T_e $$

where T _e is the encounter period.

Since the purpose of this study is the control of underactuated vessels, the general Six Degree of Freedom (6-DOF) motions can be reduced to the motion in sway and yaw (that is, horizontal motions) under the assumption that the only available control input is rudder angle. For the energy radiated from the ship sideways, we obtain the energy radiated by the horizontal motions:

(15)$$P_{26} =\displaystyle{{\pi }\over{\omega _e }}\vint_L {b_{26} \vert {U_{RY} } \vert ^2\hbox{d}x} $$

where ω _e is the encounter frequency, L is the ship length and b ₂₆ is the motion damping coefficient. U _RY is the transverse relative velocity of each ship section. Hence, the formulae to calculate the added resistance and drift force in oblique waves are:

(16)$$( -R_T ) ( -c-U\cos \beta ) T_e =P_{26} $$

(17)$$\vert {R_x } \vert =\vert {R_T \cos \beta } \vert $$

(18)$$\vert {R_y } \vert =\vert {R_T \sin \beta } \vert $$

where P ₂₆ is the energy radiated by sway-yaw motion during one encounter period.

In the calm water case, when the ship travels with a yaw angle, it will create a pressure difference between the port and starboard of the ship and lead to an extra resistance, which gives a contribution to the longitudinal response of the ship. The longitudinal component of Newton's second law in the body fixed coordinate system is:

(19)$$M( {\dot u}-vr) =X_{\dot u} {\dot u}-R_T ( u) +( 1-\tau) T( u\comma \; n) +X_{vv} v^2+X_{vr} vr+X_{rr} r^2+X_{\delta \delta } \delta ^2 $$

where $X_{\dot u} $, X _vv,… are used to express longitudinal hydrodynamic forces on the hull and the rudder. R _T (u) is the ship calm resistance, τ is the thrust deduction coefficient, T(u, n) is the propeller revolutions per second and δ denotes the rudder angle.

The main cause of resistance increase in a turning motion is due to the terms Mvr and X _vr vr. Both R _T (u) and T(u, n) are independent of the steering motion. The terms X _vv and X _rr are defined to be zero if the ship is symmetrical according to potential flow theory. The term X _{δ δ}δ ² is the extra drag due to rudder angle and can be neglected because of its small value.

So, the added resistance due to yawing can be simplified to the following:

(20)$$R_{yaw} =( M+X_{vr} ) vr $$

where M is the ship mass and this resistance will be decided by the value of sway rate and yaw angle velocity.

Assumption 3: The measurement of sway can be obtained either by using a compact 6-DOF sensor or Fibre Optic Gyroscopes (FOGs) to provide continuous and accurate measurement of 6-DOF motions (MARIN, 2014).

3.1. Adaptive sliding mode controller

Considering Equation (8), the nonlinear steering system contains a nonlinear term; the feedback linearization method is adopted to convert the nonlinear part. By defining a new control signal u, the rudder angle order can be treated as follows:

(21)$$\delta _c =\displaystyle{{1}\over{a_3}}( u-a_2 r\vert r \vert ) $$

Then the nonlinear steering system is expressed in a linear equation:

(22)$$\dot {r}=a_1 ( U_0 ) r+u+d $$

The state space format is:

(23)$$\dot{{\bf x}}={\bf Ax}+{\bf B}u+d $$

where ${\bf x}=\left[\matrix{ r \cr \psi}\right]$, ${\bf A}=\left[\matrix{ {a_{1} } &0 \cr 1 &0}\right]$, ${\bf B}=\left[\matrix{ 1 \cr 0}\right]$.

Let the reference state be x_d = [0, ψ _d] and $\dot {{\bf x}}_{{\rm d}} =0$. A sliding manifold is used to obtain the control law and is defined as:

(24)$$s={\bf h}^{\rm T}{\bf x}_{{\rm e}} ={\bf h}^{\rm T}( {\bf x}-{\bf x}_{{\rm d}}) $$

where h = [h ₁, h ₂] ^T is a right eigenvector of A_c (i.e. ${\bf A}_{{\rm c}}^{\rm T} {\bf h}=\lambda {\bf h}) $. The weighting vector h is selected by computing the equation ${\bf A}_{{\rm c}}^{\rm T} {\bf h}=0$ for λ = 0 (Healey and Lienard, Reference Healey and Lienard1993).

In the SMC system, a feedback control law is written as:

(25)$$u=-{\bf kx}+u_0 $$

where the first item of the controller is a state feedback control law (that is, an equivalent controller), the second term is a nonlinear switching control law.

Substituting Equation (25) into Equation (23), we obtain:

(26)$$\dot {{\bf x}}={\bf A}_{{\rm c}} {\bf x}\hbox{+}{\bf B}u_0 +d $$

where A_c = A − Bk^T is the combined state matrix, k = [k, 0] ^T is the feedback gain vector and the zero gain in k represents the integration in the yaw angle channel.

The nonlinear switching control law to reject the disturbance is chosen as follows:

(27)$$u_0 =-( {\bf h}^{\rm T}{\bf B}) ^{-1}[ {\bf h}^{\rm T}\hat {d}+\eta \hbox{sgn}( s) ] $$

where ${\hat d}$ is the estimate of d.

Differentiating the sliding surface function, then:

(28)$$\dot {s}={\bf h}^{\rm T}{\bf A}_{{\rm c}} {\bf x}+{\bf h}^{\rm T}{\bf B}u_0 +{\bf h}^{\rm T}d-{\bf h}^{\rm T}\dot {{\bf x}}_{{\rm d}} =\lambda {\bf x}^{\rm T}{\bf h}-\eta \hbox{sgn}( s) +{\bf h}^{\rm T}\Delta d=-\eta \hbox{sgn}( s) +{\bf h}^{\rm T}\Delta d $$

where λ x^Th = 0 if h is a right eigenvector, and ${\bf h}^{\rm T}\dot {{\bf x}}_{{\rm d}} =0$ because the reference signal x_d is a constant vector. The parameter $\Delta d=d-\hat {d}$ is the estimation error.

Since the disturbance is unknown, a better guess for it is ${\hat d}=0$. Hence, the nonlinear switching control law becomes:

(29)$$u_0 =-( {\bf h}^{\rm T}{\bf B}) ^{-1}\eta \hbox{sgn}( s) $$

where η > d _max‖h‖.

Equation (29) leads to the sliding mode controller:

(30)$$\delta _c =\displaystyle{{1}\over{a_3}}[ -kr-\displaystyle{{1}\over{h_1}}\eta \hbox{sgn}( s/\varphi) -a_2 r\vert r \vert ] $$

In this controller, a larger switching gain η corresponds to a shorter time to reach s = 0 and the system robustness against the environmental disturbance is proven. However, the upper bound of disturbance is often difficult to find in practice. So, the adaptive method is adopted to tune the controller gain without knowledge about the disturbance. Still considering the estimate of the disturbance is zero, the robust switching control law of the total controller is modified as:

(31)$$u_0 =-( {\bf h}^{\rm T}{\bf B}) ^{-1}\hat {\eta }\hbox{sgn}( s) $$

where ${\hat \eta}$ is the estimate of the adjustable gain and is a positive value.

The adaptation law is written as:

(32)$$\dot {\hat {\eta }}=\displaystyle{{1}\over{\alpha}}\vert s \vert $$

where α > 0 is the adaptation gain.

Then the differentiation of the sliding surface is:

(33)$$\dot {s}={\bf h}^{\rm T}( \dot {{\bf x}}\hbox{-}\dot {{\bf x}}_{{\rm d}} ) ={\bf h}^{\rm T}( {\bf A}_{{\rm c}} {\bf x}+{\bf B}u_0 +d) ={\bf h}^{\rm T}{\bf B}u_0 +{\bf h}^{\rm T}d=-\hat {\eta }\hbox{sgn}( s) +{\bf h}^{\rm T}d $$

Selecting the Lyapunov function,

(34)$$V=\displaystyle{{1}\over{2}}s^2+\displaystyle{{1}\over{2}}\alpha \tilde {\eta }^2 $$

where $\tilde {\eta }=\hat {\eta }-\eta $ is the estimation error.

(35)$$\eqalign{\dot {V}&=\displaystyle{{1}\over{2}}s\dot {s}+\displaystyle{{1}\over{2}}\alpha \tilde {\eta }^2=s[ {\bf h}^{\rm T}d-\hat {\eta }\hbox{sgn}( s) ] +\alpha ( \hat {\eta }-\eta ) \dot {\hat {\eta }} \cr & = s{\bf h}^{\rm T}d-\hat {\eta }\vert s \vert +( \hat {\eta }-\eta ) \vert s \vert =s{\bf h}^{\rm T}d-\eta \vert s \vert \le 0} $$

Equation (35) holds for all time for the external disturbance. The parameter ${\hat \eta}$ chosen in Equation (31) implies that the system trajectory will move and reach the sliding surface in a finite time and the sliding surface declines to zero, so the control law given by Equation (25) guarantees the sliding mode sustained (Huang et al., Reference Huang, Kuo and Chang2013). The tanh function has replaced the signum function to attenuate the chattering effect. Hence, the Adaptive Sliding Mode Controller (A-SMC) is:

(36)$$\delta _c =\displaystyle{{1}\over{a_3}}[ -kr-\displaystyle{{1}\over{h_1 }}\hat {\eta }\tanh ( s/\varphi ) -a_2 r\vert r \vert ] $$

where φ is the boundary layer thickness.

3.2. Nonlinear disturbance observer

In the last section, the estimate of disturbance is treated as zero because no knowledge of the disturbance can be available. An alternative method to process this problem is called a Disturbance Observer (DO). Since the yaw acceleration is not easy to obtain, it is also difficult to construct the acceleration signal from yaw rate by differentiation. So, a modified observer called a Nonlinear Disturbance Observer (NDO) is adopted here (Chen, Reference Chen2004; Yang et al., Reference Yang, Li and Yu2013).

Define a variable:

(37)$$z=\hat {d}-p( r\comma \; \psi ) $$

Let the function p(r, ψ) be given by the following equation:

(38)$$\displaystyle{{\hbox{d}p}\over{\hbox{d}t}}=\displaystyle{{L( r\comma \; \psi ) }\over{a_4 }\dot {r}} $$

Define the observer error signal:

(39)$$\tilde {d}=d-\hat {d} $$

According to the linear disturbance observer, a DO is proposed as:

(40)$$\eqalign{\dot {\hat {d}} & =L( r\comma \; \psi ) ( d-\hat {d}) =L( r\comma \; \psi ) ( \dot {r}-a_1 ( U_0 ) r-a_2 r\vert r \vert -a_3 \delta ) -L( r\comma \; \psi ) \hat {d} \cr & = L( r\comma \; \psi ) ( \dot {r}-a_1 ( U_0 ) r-a_2 r\vert r \vert -a_3 \delta ) -L( r\comma \; \psi ) [ z+p( r\comma \; \psi ) ] \cr &=-L( r\comma \; \psi ) z+L( r\comma \; \psi ) [ \dot {r}-a_1 ( U_0 ) r-a_2 r\vert r \vert -a_3 \delta -p( r\comma \; \psi ) ] } $$

In general, prior information about the derivative of the disturbance is unavailable, and the disturbance varies slowly relative to the observer dynamics. Then it reasonable to suppose that:

(41)$$\dot {d}=0 $$

Hence, the derivative of the observer error is:

(42)$$\dot {\tilde {d}}=\dot {d}-\dot {\hat {d}}=-\dot {\hat {d}}=-L( r\comma \; \psi) \tilde {d} $$

Let the functions in Equation (36) be chosen as:

(43)$$L( r\comma \; \psi ) =a $$

where a is a positive constant.

Then the following equations are obtained:

(44)$$p( r\comma \; \psi ) = ar $$

(45)$$\displaystyle{{\hbox{d}p}\over{\hbox{d}t}} = a\dot {r} $$

Combining the above equations, the update law can be written as:

(46)$$\eqalign{\dot {z} & =\dot {\hat {d}}-\displaystyle{{\hbox{d}p}\over{\hbox{d}t}}=\dot {\hat {d}}-L( r\comma \; \psi ) \dot {r} \cr &=-L( r\comma \; \psi ) z+L( r\comma \; \psi ) [ -a_1 ( U_0 ) r-a_2 r\vert r \vert -a_3 \delta -p( r\comma \; \psi ) ] \cr &=-az+a\lcub -[ a+a_1 ( U_0 ) ] r-a_2 r\vert r \vert -a_3 \delta \rcub } $$

Hence, the NDO is given by:

(47)$$\hat {d}=z+ar $$

The nonlinear switching control law in this SMC is defined as:

(48)$$u_0 =-( {\bf h}^{\rm T}{\bf B}) ^{-1}[ {\bf h}^{\rm T}\hat {d}+\eta _0 \hbox{ sgn}( s) ] $$

The differentiation of the sliding surface is:

(49)$$\dot {s}={\bf h}^{\rm T}{\bf A}_{{\rm c}} {\bf x}-\eta _0 \hbox{ sgn}( s) +{\bf h}^{\rm T}d-{\bf h}^{\rm T}\hat {d}={\bf h}^{\rm T}{\bf A}_{{\rm c}} {\bf x}-\eta _0 \hbox{ sgn}( s) +{\bf h}^{\rm T}\tilde {d}=-\eta _0 \hbox{ sgn}( s) +{\bf h}^{\rm T}\tilde {d} $$

Let the Lyapunov function be chosen as:

(50)$$V=\displaystyle{{1}\over{2}}s^2+\displaystyle{{1}\over{2a}}\tilde {d}^2 $$

Differentiation of the Lyapunov function with respect to the trajectory of states gives:

(51)$$\dot {V}=s\dot {s}+\displaystyle{{1}\over{a}}\tilde {d}\dot {\tilde {d}}=-\eta _0 \hbox{ sgn}( s) +{\bf h}^{\rm T}\tilde {d}s-a\tilde {d}^2\le 0 $$

where the gain must satisfy $\eta _{0} \gt \Vert {\bf h} \Vert \cdot \Vert {\tilde {d}} \Vert $.

Hence, the final output of the sliding mode controller with nonlinear disturbance observer (NDO-SMC) is written as:

(52)$$\eqalign{\delta _c & =\displaystyle{{1}\over{a_3}}[ -kr-\hat {d}-\displaystyle{{1}\over{h_1}}\eta_0 \tanh ( s/\varphi ) -a_2 r\vert r \vert ] \cr & = \displaystyle{{1}\over{a_3}}[ -kr-\displaystyle{{1}\over{h_1}}\eta _0 \tanh ( s/\varphi ) -a_2 r\vert r \vert ] -\displaystyle{{1}\over{a_3}}\hat {d}} $$

where the first term in Equation (52) can be treated as a sliding mode controller that rejects the disturbance ${\tilde d}$. A diagram of this control system is shown in Figure 4.

Figure 4.

Modified autopilot control system.