Path following of underactuated surface ships based on model predictive control with neural network

Abstract

A model predictive control approach is proposed for path following of underactuated surface ships with input saturation, parameters uncertainties, and environmental disturbances. An Euler iterative algorithm is used to reduce the calculation amount of model predictive control. The matter of input saturation is addressed naturally and flexibly by taking advantage of model predictive control. The mathematical model group (MMG) model as the internal model improves the control accuracy. A radial basis function neural network is also applied to compensate the total unknowns including parameters uncertainties and environmental disturbances. The numerical simulation results show that the designed controller can force an underactuated ship to follow the desired path accurately in the case of input saturation and time-varying environmental disturbances including wind, current, and wave.

Keywords

Model predictive control path following underactuated surface ships neural network

Introduction

Path following of surface ships has been a long-standing control problem that has attracted attention from the control community for many years.¹ As a matter of fact, the path following control possesses the significant position in some marine applications, such as seawater monitoring for the fixed points, laying and maintaining the underwater pipelines, even tracking control of the ordinary merchant ships. The underactuated nature of these problems, namely with more variables to be controlled than the number of control actuators, coupled with the matter of uncertain parameters, external disturbances as well as input constraint, renders the control problem both challenging and interesting.^2,3

The path following control problem of underactuated ships has been tackled in a number of research studies in the last few years. For example, a position prediction model was presented by Nagai and Watanabe,⁴ which estimated positions by the integral in continuous time and predicted future positions by considering only constant velocities. In the study of Shen and Dai,⁵ the iterative sliding mode was designed by a hyperbolic tangent function, and both neural network (NN) and reinforced learning were used to inhibit the chattering of the control input.

Methods to deal with the uncertain parameters or external disturbances could also be found in previous studies. In the study of Wang et al.,⁶ a robust and adaptive control algorithm was proposed to cope with the random disturbances and uncertain parameters, in which the Serret–Frenet frame and output-redefinition were introduced to simplify the controller design. In the study of Herman and Adamski,⁷ the robustness was promoted by designing a non-adaptive velocity tracking algorithm for fully actuated vehicles. Li et al.⁸ introduced an active disturbance rejection control (ADRC) controller with sliding mode which can solve the path following problem of uncertainties including internal dynamic and external disturbances. In the study of Lekkas and Fossen,⁹ a sliding mode controller was employed to stabilize the heading angle, and the constant external disturbances could be handled. It is noted that the sliding mode and incremental feedback methods were used by Liu et al.¹⁰ to avoid effectively the trajectory tracking with unknown disturbances such as wind or current. In the studies of Mu et al.¹¹ and Liang et al.,¹² both surge control and yaw control were considered by designing the sliding mode trajectory tracking controller, and the unknown dynamics or disturbances were packaged for compensation via NN together with minimum learning parameter method (MLP). The ADRC was proposed by Huang and Fan¹³ to deal with the path following problem, and the extended state observer (ESO) and adaptive laws were used to estimate the external disturbances and sideslip angle, respectively. In the study of Vo et al.,¹⁴ an adaptive autopilot based on backstepping methodology was designed to control ship heading with external disturbances. In the study of Zhang et al.,¹⁵ the nonlinear disturbance observer was designed to estimate disturbances, and a predictive controller was offered based on Serret–Frenet frame, but the advantage of predictive control in saturated suboptimal solutions was not utilized. To solve the problems of both the uncertain parameters and environmental disturbances, the robust and adaptive control algorithm was also used by Tong and Li.^16,17 An adaptive takagi and sugeno (T-S) fuzzy NN control method was proposed by Dong et al.¹⁸ to deal with the external disturbance, but the simulation model was very different from a real ship. Model uncertainties were approximated and compensated by virtue of the concise radial basis function (RBF) adaptive NN in the study of Zhang et al.¹⁹

Some investigations to address the problem of input saturation could be found in previous research studies. Zheng et al.²⁰ proposed a linear model predictive control (MPC) and a nonlinear MPC to demonstrate the advantage of MPC on dealing with the matter of saturation. A novel robust MPC was presented by Li et al.²¹ to solve the vessel’s roll constraints. In the study of Shen and Jing,²² both neuron adaptive and iterative sliding mode strategy were developed for the path tracking of underactuated ships. In the study of Yu et al.,²³ the input saturation was considered via an auxiliary system and the MLP was provided to approximate the unknown dynamic without external disturbances. A backstepping method was designed by Zheng and Feroskhan,²⁴ and the disturbance observer and auxiliary system were used to estimate external disturbances and cope with the input saturation, respectively. In the study of Qiu et al.,²⁵ the external disturbances, dynamical uncertainties, and input saturation were settled effectively by the adaptive technology, MLP, and auxiliary dynamic system, respectively.

In the aforementioned studies, the problem of unknown dynamical or external disturbances could be addressed effectively. However, few have taken into account the input saturation, whereas the poor control performances or even actuator damages would happen in real implementation if these constraints were neglected in the controller design.²⁰ Some research studies have pointed out the solutions of the saturation such as introducing an iterative sliding mode theory^5,10,22 or an auxiliary system.^23
–25 Comparing to the schemes above, the MPC methodology proposed in this article is more natural and flexible^20,26 as just the value of the constrained rudder angle is considered directly in case of solving a quadratic programming (QP)²⁷ problem. The designed MPC is simpler than those proposed in the literature,^6,15,20,21 because the Euler iterative algorithm is introduced to predict the future states. The MMG model is used to improve the prediction accuracy as the internal model. Moreover, with reference to the literature,^19,23,25,28 a RBF NN is used to compensate the uncertainties for the improvement of the internal model accuracy and the controller robustness.

The remainder of this article is organized as follows. The kinematic and dynamic model of the surface ship is presented in the next section. The MPC with RBF NN method is developed in the third section. The fourth section provides the stability analysis. Simulation results are presented and analyzed in the fifth section, followed by conclusions.

An underactuated surface ship model

The ship position in the horizontal plane and the motion parameters with ocean current disturbance are shown in Figure 1. The x-axis and y-axis point to the true North and East, respectively. ϕ is the heading angle from the x-axis. u and v are the surge and sway velocity over the ground, and V = (u ² + v ²)^1/2 is the ship speed over ground. r is the yaw rate. V_c and ϕ_c are the ocean current speed and direction, respectively. u_r and v_r are the longitudinal and lateral speed through water, respectively. $Β = arctan (v / u)$ is the sideslip angle. δ denotes the rudder angle, where $|δ| \leq 35 °$ is the input saturation.

Figure 1.

Ship’s position and kinematics under current disturbance.

The ship’s MMG model was proposed by Ship Manoeuvring Mathematical Model Group⁸ in the late 1970s, which took the rudder angle as the control input. Considering the wind, ocean current and wave disturbances, as well as servo response characteristics, a ship’s MMG model could be expressed as

\{\begin{cases} \dot{x} = u_{r} cos ϕ - v_{r} sin ϕ + V_{c} cos ϕ_{c} = u cos ϕ - v sin ϕ \\ \dot{y} = u_{r} sin ϕ + v_{r} cos ϕ + V_{c} sin ϕ_{c} = u sin ϕ + v cos ϕ \\ \dot{ϕ} = r \\ \dot{u} = [(m + m_{y}) v r + X_{H} + X_{P} + X_{R} + X_{W} + X_{Wave} \\ + (m_{x} - m_{y}) V_{c} sin (ϕ_{c} - ϕ) r] / (m + m_{x}) \\ \dot{v} = [(m + m_{x}) u r + Y_{H} + Y_{P} + Y_{R} + Y_{W} + Y_{Wave} \\ - (m_{x} - m_{y}) V_{c} cos (ϕ_{c} - ϕ) r] / (m + m_{y}) \\ \dot{r} = (N_{H} + N_{P} + N_{R} + N_{W} + N_{Wave}) / (I_{z z} + J_{z z}) \\ \dot{δ} = K_{E} (δ_{r} - δ) / T_{E} \end{cases}

where m is the ship mass. m_x , m_y , $I_{z z}$ , and $J_{z z}$ are the corresponding added masses and added moment of inertia, respectively. X and Y denote the forces in x and y frames, N denotes the moment, and W denotes the wind. It is noted that only the slowly varying forces are counteracted by the steering system. The oscillatory motion caused by the first-order wave-induced forces should be prevented from entering the feedback loop. H, P, and R denote the low-frequency hydrodynamic forces on the hull, propeller, and rudder, respectively. $δ_{r}$ is the command rudder angle, K_E is the steering gear control gain, and T_E is the servo time constant. And X_R , Y_R , and N_R are the rudder forces (moments) given by

\{\begin{cases} X_{R} = - (1 - t_{R}) F_{N} sin δ \\ Y_{R} = - (1+ a_{H}) F_{N} cos δ \\ N_{R} = - (x_{R} + a_{H} x_{H}) F_{N} cos δ \end{cases}

where t_R is the factor of rudder resistance deductions, a_H is the ratio of the hull’s additional lateral force to the rudder lateral force caused by steering, x_H is the distance from the steering induced hull lateral force center to the ship’s center of gravity, and F_N is the rudder positive pressure.

The winds (moments) X_W , Y_W , and N_W in the dynamics (1) are given²⁹ by

\{\begin{cases} X_{W} = \frac{1}{2} ρ_{a} A_{f} U_{R}^{2} C_{w x} (α_{R}) \\ Y_{W} = \frac{1}{2} ρ_{a} A_{s} U_{R}^{2} C_{w y} (α_{R}) \\ N_{W} = \frac{1}{2} ρ_{a} A_{s} L_{o a} U_{R}^{2} C_{w n} (α_{R}) \end{cases}

where $ρ_{a}$ is the air density, $α_{R}$ is the wind angle of attack, U_R is the relative wind speed, and A_f and A_s are the frontal and lateral projected areas, respectively. L_oa is the length overall, and C_wx (α_R ), C_wy (α_R ), and C_wn (α_R ) are the nondimensional wind coefficients, respectively.

The waves (moments) X _Wave, Y _Wave, and N _Wave are given³⁰ by

\{\begin{cases} X_{Wave} = \frac{1}{2} ρ L α^{2} cos χ C_{X w} (λ) \\ Y_{Wave} = \frac{1}{2} ρ L α^{2} sin χ C_{Y w} (λ) \\ N_{Wave} = \frac{1}{2} ρ L α^{2} sin χ C_{N w} (λ) \end{cases}

where λ is the wavelength, χ is the encounter angle, ρ is the seawater density, α is the wave amplitude, L is the ship length, and C_Xw (λ), C_Yw (λ), and C_Nw (λ) are coefficients related to wave length λ, respectively.

Control design and stability analysis

The structure of the control

In this article, the path following controller is designed by MPC that consists of internal model, feedback correction, and dynamic optimization. Its main advantage is that the input constraint could be handled explicitly^26,31 because MPC considers a standard QP problem with saturation only. First of all, the future states over the prediction horizon N_p can be predicted according to the internal model and current states of the system. Then, the cost function including prediction errors would be built via the future states. Finally, the optimal control input will be obtained by solving the cost function with constraint. At the next sampling time, the same procedure is repeated on the basis of the feedback of the new measurements from the system, so that is also called rolling optimization.

Path following pertains to following a desired path independent of time, without imposing any restrictions on the temporal propagation along the path. Thus, the revolution of propeller has been set a constant in this article. Therefore, the control objective becomes that a suitable rudder angle δ is designed to force the ship to follow the reference path at the certain propeller revolutions. Moreover, the MMG model is employed directly as the internal model for the prediction of states on the basis of Euler iterative algorithm. The structure of the control is shown in Figure 2.

Figure 2.

The structural diagram of the control.

As shown in Figure 2, the output of ship motion system is the current states at the current sampling time, and the unknowns including parameters and external disturbances are approximated by RBF NN on the basis of the previous states. Then the internal model would predict the future states of the system over the prediction horizon N_p based on internal model (MMG) and current states. After that, the cost function will be established with respect to the prediction error e. Finally, the optimal rudder angle with saturation will be obtained by solving the cost function in the optimal solution module. At the next sampling time, the same procedure will be repeated upon the feedback of the new measurements from the ship motion system.

Assumptions

The assumptions in this article are as follows:

The states of ship such as x, y, ϕ, u, v, and r can be measured.

The derivative with respect to time of historical states such as $\dot{u}$ , $\dot{v}$ , and $\dot{r}$ could be obtained.

The unknowns d_u , d_v , and d_r are bounded, namely, $|d_{u}| \leq d_{u} _{max}$ , $|d_{v}| \leq d_{v} _{max}$ , and $| d_{r} | \leq d_{r} _{max}$ .

Remark

In assumption Ⅰ, the states could be measured through the sensors or observers in the navigation system.²⁰ Then, assumption Ⅱ can be satisfied by collecting the historical data. Assumption Ⅲ is reasonable and appeared in previous study.¹⁰

MPC controller

First of all, the prediction equations of the states may be employed according to the Euler iteration algorithm³² and equation (1)

\{\begin{cases} \bar{x} (k + 1) = \bar{x} (k) + T_{c} * {\dot{f}}_{x} \\ \bar{y} (k + 1) = \bar{y} (k) + T_{c} * {\dot{f}}_{y} \\ \bar{ϕ} (k + 1) = \bar{ϕ} (k) + T_{c} * {\dot{f}}_{ϕ} \\ \bar{u} (k + 1) = \bar{u} (k) + T_{c} * {\dot{f}}_{u} \\ \bar{v} (k + 1) = \bar{v} (k) + T_{c} * {\dot{f}}_{v} \\ \bar{r} (k + 1) = \bar{r} (k) + T_{c} * {\dot{f}}_{r} \\ \bar{δ} (k + 1) = \bar{δ} (k) + T_{c} * {\dot{f}}_{δ} \end{cases}

where T_c is the sampling time. ${\dot{f}}_{x}$ , ${\dot{f}}_{y}$ , ${\dot{f}}_{ϕ}$ , ${\dot{f}}_{u}$ , ${\dot{f}}_{v}$ , ${\dot{f}}_{r}$ , and ${\dot{f}}_{δ}$ are the discrete values of equation (1), respectively

\{\begin{cases} {\dot{f}}_{x} = u_{r} (k) cos (ϕ (k)) - v_{r} (k) sin (ϕ (k)) \\ {\dot{f}}_{y} = u_{r} (k) sin (ϕ (k)) + v_{r} (k) cos (ϕ (k)) \\ {\dot{f}}_{ϕ} = r (k) \\ {\dot{f}}_{u} = [(m + m_{y}) v (k) r (k) + X_{H} (k) + X_{P} (k) + X_{R} (k) + {\hat{d}}_{u} (k)] / (m + m_{x}) \\ {\dot{f}}_{v} = [(m + m_{x}) u (k) r (k) + Y_{H} (k) + Y_{P} (k) + Y_{R} (k) + {\hat{d}}_{v} (k)] / (m + m_{y}) \\ {\dot{f}}_{r} = (N_{H} (k) + N_{P} (k) + N_{R} (k) + {\hat{d}}_{r} (k)) / (I_{z z} + J_{z z}) \\ {\dot{f}}_{δ} = K_{E} (δ_{r} (k) - δ (k)) / T_{E} \end{cases}

where ${\hat{d}}_{u} (k)$ , ${\hat{d}}_{v} (k)$ , and ${\hat{d}}_{r} (k)$ are the estimated values of the total unknowns including parameters uncertainties and environmental disturbances. At the current sampling time k, we can predict the future lateral displacements over the prediction horizon N_p according to equation (5) as follows

\{\begin{cases} \bar{y} (k + 1) = y (k) + \sum_{j = 1}^{1} (T_{c} * {\dot{f}}_{y}) \\ \bar{y} (k + 2) = y (k) + \sum_{j = 1}^{2} (T_{c} * {\dot{f}}_{y}) \\ \dots \\ \bar{y} (k + N_{p}) = y (k) + \sum_{j = 1}^{N_{p}} (T_{c} * {\dot{f}}_{y}) \end{cases}

Then, the prediction errors ${\bar{e}}_{y} (k + j), j = 1, 2, \dots, N_{P}$ could be calculated via equation (7) and the reference path y _d

\{\begin{cases} {\bar{e}}_{y} (k + 1) = \bar{y} (k + 1) - {\bar{y}}_{d} (k + 1) \\ {\bar{e}}_{y} (k + 2) = \bar{y} (k + 2) - {\bar{y}}_{d} (k + 2) \\ \dots \\ {\bar{e}}_{y} (k + N_{p}) = \bar{y} (k + N_{p}) - {\bar{y}}_{d} (k + N_{p}) \end{cases}

At the sampling time k, the cost function can be established by the prediction errors

min_{δ_{min} \leq δ \leq δ_{max}} f^{*} ({\bar{e}}_{y}, δ) = \sum_{j = 1}^{N_{p}} \bar{e} {(k + j)}^{T} Q \bar{e} (k + j)

where Q is the weighting matrix. The optimal control law $δ^{*} (k + 1)$ will be calculated by solving the QP problem and the partial differential equation $\partial f^{*} ({\bar{e}}_{y}, δ) / \partial δ = 0$ with the constraint $δ_{min} \leq δ \leq δ_{max}$ . At each time over the prediction horizon, in most of the previous MPC applications, all the corresponding future inputs $δ_{k} _{+ 1}, δ_{k} _{+ 2}, \dots, δ_{k} _{+ N_{P}}$ are considered during predicting future states, whereas the first optimal input $δ_{k + 1}^{*}$ only is chosen as the final control law to the system. However, we only consider a common future input over the prediction horizon, namely $δ_{k} _{+ 1} = δ_{k} _{+ 2} = \dots = δ_{k} _{+ N_{P}} = δ$ here. And the control law is adjusted online without the analytical expression because MPC depends on the modern computer.²⁶ The RBF NN is used to compensate the unknowns such as ${\hat{d}}_{u} (k)$ , ${\hat{d}}_{v} (k)$ , and ${\hat{d}}_{r} (k)$ that they are included in $f^{*} ({\bar{e}}_{y}, δ)$ . N_p , N_c , and T_c , namely prediction horizon, control horizon, and prediction sampling time, respectively, are control parameters to be tuned.

Radial basis function neural network

The gradient descent method of RBF NN is used to approximate the unknowns including ${\hat{d}}_{u} (k)$ , ${\hat{d}}_{v} (k)$ , and ${\hat{d}}_{r} (k)$ on the basis of the previous information. RBF NN can approximate both the linear and nonlinear continuous functions over a compact set.²⁸ Taking ${\hat{d}}_{r} (k)$ as an example, the RBF NN method is given by

y_{n} = w_{1} h_{1} + w_{2} h_{2} + \dots + w_{n} h_{n}

where y_n is output vector, and here is ${\hat{d}}_{r} (k)$ . n is the NN node number, n > 1, and w = [w ₁, w ₂,…, w _n] is weight vector. The vector of RBF h is given by

h_{j} = exp (- \frac{{‖x - c_{j}‖}^{2}}{2 b_{j}^{2}}), j = 1, 2, \dots, n

where x is the input vector, namely $[δ (k - 1), \dot{r} (k - 1)]$ , c_j is the center of the receptive field, and b_j is the width of Gaussian function

\{\begin{cases} w_{j} (t) = w_{j} (t - 1) + Δ w_{j} (t) + α (w_{j} (t - 1) - w_{j} (t - 2)) \\ b_{j} (t) = b_{j} (t - 1) + Δ b_{j} (t) + α (b_{j} (t - 1) - b_{j} (t - 2)) \\ c_{j} (t) = c_{j} (t - 1) + Δ c_{j} (t) + α (c_{j} (t - 1) - c_{j} (t - 2)) \end{cases}

where $α \in (0, 1)$ is momentum factor, the method to determine w, b, and c is given by

\{\begin{cases} Δ w_{j} = - η \frac{\partial E}{\partial w_{j}} = η (y_{h} - y_{n}) h_{j} \\ Δ b_{j} = - η \frac{\partial E}{\partial b_{j}} = η (y_{h} - y_{n}) w_{j} h_{j} \frac{{‖x - c_{j}‖}^{2}}{b_{j}^{3}} \\ Δ c_{j i} = - η \frac{\partial E}{\partial c_{j i}} = η (y_{h} - y_{n}) w_{j} h_{j} \frac{x_{j} - c_{j i}}{b_{j}^{2}} \end{cases}

where $η \in (0, 1)$ is learning rate and y_h is previous information. w, b, and c could be adjusted online with gradient descent method. The acquisition processes of ${\hat{d}}_{u} (k)$ , ${\hat{d}}_{v} (k)$ are the same as that of ${\hat{d}}_{r} (k)$ .

Stability analysis

Lemma 1

The continuous function will obtain minimum value if it has a point ξ where the derivative equals zero (g′(ξ) = 0) and the left and right derivative symbols of this point are opposite (g′(ξ − Δ) < 0 and g′(ξ + Δ) > 0).³³

From equations (1) to (6), we can predict the value of $\bar{r} (k + h)$

\begin{matrix} \bar{r} (k + h) = r (k) + T_{c} \sum_{n = 1}^{h - 1} [\frac{{\bar{N}}_{R} (k + n)}{I_{z z} + J_{z z}} + {\bar{f}}_{r}] \\ = r (k) + T_{c} \sum_{n = 1}^{h - 1} [- F_{r} sin (δ - γ β_{R}) cos δ + {\bar{f}}_{r}] \end{matrix}

where F_r is a positive constant without δ, ${\bar{f}}_{r}$ is the sum of other terms without δ, and γ is the rectification coefficient. γ and β_R ³⁴ could be given by

\{\begin{cases} γ = - 22.2 {(C_{b} B / L)}^{2} + 0.02 (C_{b} B / L) + 0.68 \\ β_{R} = β + r L / U \end{cases}

where C_b is the square coefficient and B is the ship breadth. It is reasonable to suppose that the sideslip angle β is not very large and velocity is almost constant, so let $|β| \leq 25^{\circ}$ , $U \leq 7.5 m/s$ , and $|r| \leq 0.02 rad/s$ here, then the derivative of equation (14) with respect to δ yields

\frac{\partial \bar{r} (k + h)}{\partial δ} = T_{c} \sum_{n = 1}^{h - 1} - F_{r} cos (2 δ \pm {17.8}^{\circ})

From equations (5) and (6), the value of $\bar{ϕ} (k + j)$ could be predicted

\bar{ϕ} (k + j) = ϕ (k) + T_{c} \sum_{h = 1}^{j - 1} \bar{r} (k + h)

The derivative of equation (17) with respect to δ, and substituting equation (16) into (17) yields

\frac{\partial \bar{ϕ} (k + j)}{\partial δ} = T_{c} \sum_{h = 1}^{j - 1} \frac{\partial \bar{r} (k + h)}{\partial δ} = T_{c} \sum_{h = 1}^{j - 1} T_{c} \sum_{n = 1}^{h - 1} - F_{r} cos (2 δ \pm 17.8)

The second equation of (5) can be reformulated as follows

\bar{y} (k + i) = y (k) + T_{c} \sum_{j = 1}^{i - 1} \sqrt{{\bar{u}}^{2} (k + j) + {\bar{v}}^{2} (k + j)} sin (\bar{ϕ} (k + j) + \bar{β} (k + j))

where $\bar{β} (k + j) = arctan (\bar{v} (k + j) / \bar{u} (k + j))$ .

Combining equations (8) and (9), the cost function will be built as follows

\begin{matrix} f = Q \sum_{i = 1}^{N_{p}} {[\bar{y} (k + i) - {\bar{y}}_{d} (k + i)]}^{2} \\ = Q \sum_{i = 1}^{N_{p}} {[y (k) + T_{c} \sum_{j = 1}^{i - 1} \sqrt{{\bar{u}}^{2} (k + j) + {\bar{v}}^{2} (k + j)} sin (\bar{ϕ} (k + j) + \bar{β} (k + j)) - {\bar{y}}_{d} (k + i)]}^{2} \end{matrix}

${\bar{y}}_{d}$ does not contain δ, and the component of rudder force in u or v is smaller compared to angular rate r,⁴ so the partial derivative of equation (20) with respect to δ yields

\frac{\partial f}{\partial δ} = 2 Q \sum_{i = 1}^{N_{p}} (z_{1} \cdot z_{2} \cdot z_{3})

where z ₁, z ₂, and z ₃ are as follows

\{\begin{cases} z_{1} = y (k) + T_{c} \sum_{j = 1}^{i - 1} \sqrt{{\bar{u}}^{2} (k + j) + {\bar{v}}^{2} (k + j)} sin θ_{e} - {\bar{y}}_{d} (k + j) \\ z_{2} = \sqrt{{\bar{u}}^{2} (k + j) + {\bar{v}}^{2} (k + j)} cos θ_{e} \\ z_{3} = \frac{\partial \bar{ϕ} (k + j)}{\partial δ} = T_{c} \sum_{h = 1}^{j - 1} T_{c} \sum_{n = 1}^{h - 1} - F_{r} cos (2 δ \pm 17.8) \\ θ_{e} = \bar{ϕ} (k + j) + \bar{β} (k + j) \end{cases}

z ₃ is always negative because $|δ| \leq 35^{\circ}$ . Then, substituting equations (14) and (17) into equation (22) yields

\{\begin{cases} θ_{e} = ϕ (k) + T_{c} \sum_{h = 1}^{j - 1} r (k) + T_{c} \sum_{h = 1}^{j - 1} T_{c} \sum_{n = 1}^{h - 1} P cos δ + \bar{β} (k + j) \\ P = - F_{r} sin (δ \pm {17.8}^{\circ}) \end{cases}

where cos(δ) is positive due to $|δ| \leq 35^{\circ}$ . The sign of θ_e is only determined by P, N_p , or T_c because the sideslip angle β is small and the heading is constant at current time. The sign of P will be plus as $δ = - 35^{\circ}$ , and θ_e will increase with the increasement of N_p or T_c . If N_p and T_c are selected suitably, then $z_{1} > 0$ , $z_{2} > 0$ , and $\partial f / \partial δ < 0$ . The sign of P would be minus as $δ = 35^{\circ}$ , and θ_e will decrease with the increasement of N_p or T_c . If N_p and T_c are selected suitably, z ₁ will be negative, whereas z ₂ is positive, and $\partial f / \partial δ > 0$ . So there are suitable parameters N_p and T_c to satisfy equation (21) as follows

\{\begin{cases} \partial f / \partial δ < 0, δ = - 35^{\circ} \\ \partial f / \partial δ > 0, δ = 35^{\circ} \end{cases}

The process of ship motion is continuous, so f is also a continuous function. According to equation (24) and the zero point theorem, there is an optimal rudder angle δ* to satisfy ∂f /∂δ* = 0. And it can fulfill ∂f /∂δ < 0 as δ < δ*, or ∂f/∂δ > 0 as δ > δ*. Thus, f will be the minimum at δ* from lemma 1. Finally, f and ${\bar{e}}_{y}$ will converge to zero gradually by obtaining the minimum of the f at each time during rolling optimization, then y(k) will also converge to y _d(k) gradually.

Numerical simulation

To illustrate the effectiveness of the designed MPC scheme, the ship “Yulong” is taken for an example in the simulation.⁸ The ship’s dynamic is given by equation (1) and the main ship parameters are as follows: daft is 8.0 m, full load displacement is 14,635 tons, ship length is 126 m, ship breadth is 20.8 m, block coefficient is 0.681, front wind projection area is 369.9 m², side wind projection area is 1031.94 m², propeller pitch is 3.66 m, rudder area is 18.8 m², servo time constant T_E is 2.5 s, and propeller revolutions per minute is set 100.

Straight line path following: The initial states u = 7.2 m/s, v = 0, r = 0, ϕ = 0, (x ₀, y ₀) = (0, 200 m), and reference path y_d = 0. The external disturbances are as follows: wind speed and direction are 10 m/s and 45° + 30° sin(0.02t), respectively. Ocean current speed and direction are 1.0 m/s and 45° + 10° sin(0.005t), respectively. Wave length, wave encounter angle, and wave amplitude are 83 m, ϕ + 135° − 30° sin(0.02t), and 3 m, respectively. The parameters of the controller are as follows: N_p = 18, N_c = 1, and T_c = 2.0 s are the parameters of MPC, initial $c_{0} = [- 1 - 0.5 0 0.5 1; - 1 - 0.5 0 0.5 1]$ , b ₀ = 5, w ₀ = 0, η = 0.15, and α = 0.05 are the parameters of RBF NN.

Curve path following: The initial states u = 7.2 m/s, v = 0, r = 0, ϕ = 0, and (x ₀, y ₀) = (0, 0). Reference path y_d = 200 sin(0.00035πx). The external disturbances conditions are the same as case 1. The parameters of the controller are as follows: N_p = 15, N_c = 1, and T_c = 2.5 s are the parameters of MPC, and the parameters of RBF NN are the same as the previous simulation.

The simulation results are shown in Figures 3 to 6, where y, ϕ, and δ denote the path, heading, and rudder angle, respectively. The heading ϕ and rudder angle δ have the time-varying deviation for the reason of the influence of the time-varying disturbances such as wave, wind, and ocean current. The rudder angle δ is smooth with more stringent limitation on the rudder amplitude constraint within 10°. This illustrates that the ability of the MPC could handle the input saturation effectively. In Figure 3, the path y has the small overshoot less than the ship breadth. The maximum path tracking fluctuation is less than 5 m or less than 20% of the ship width on passage. In short, the control scheme could force the underactuated ship to follow the desired path accurately with uncertain parameters and environmental disturbances.

Figure 3.

Straight line path following in the fixed-coordinate system.

Figure 4.

Cross tracking error, rudder and ship motion variables as straight line path following.

Figure 5.

Curve path following in the fixed-coordinate system.

Figure 6.

Tracking error, rudder and ship motion variables as curve path following.

Conclusions

This article proposed an MPC design method, aiming at a path following control for underactuated surface vessels with input saturation, parameters uncertainties, and environmental disturbances. The choice of MPC for this application is primarily motivated by the need to explicitly handle input saturation. The Euler iterative algorithm reduced the amount of calculation of MPC. The total unknowns, including parameters uncertainties and environmental disturbances, were approximated by RBF NN to enhance the system robustness. The MMG model was employed as a predictive model to improve the control accuracy. The numerical simulation illustrated the effectiveness of designed controller.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported in part by the special projects of key fields (Artificial Intelligence) of Universities in Guangdong Province (2019KZDZX1035), Innovation and Entrepreneurship Team Guidance and Cultivation “Pilot Plan” of Zhanjiang City, the National Natural Science Foundation of China (nos 51979045, 51939001, and 61976033), and program for scientific research start-up funds of Guangdong Ocean University.

ORCID iDs

Ronghui Li

Qionglei Hu

References

Liang

Wang

, et al. Swarm control with collision avoidance for multiple underactuated surface vehicles. Ocean Eng 2019; 191(106516): 1–10.

Liang

Wang

, et al. A novel distributed and self-organized swarm control framework for underactuated unmanned marine vehicles. IEEE Access 2019; 7: 703–712.

Liang

Hou

, et al. Distributed coordinated tracking control of multiple unmanned surface vehicles under complex marine environments. Ocean Eng 2020; 205(107328): 1–9.

Nagai

Watanabe

. Applying position prediction model for path following of ship on curved path. In: Proceedings of the IEEE region 10 conference, Singapore, 22–25 November 2016, pp. 22–25. IEEE.

Shen

Dai

. Iterative sliding mode control based on reinforced learning and used for path tracking of under-actuated ship. J Harbin Eng Univ 2017; 38(5): 697–704.

Wang

Zou

, et al. Path following control of under-actuated ships based on nonswitch analytic model predictive control. J Control Theory Appl 2010; 8(4): 429–434.

Herman

Adamski

. Non-adaptive velocity tracking controller for a class of vehicles. Bull Pol Acad Sci Tech Sci 2017; 65(4): 459–468.

, et al. Active disturbance rejection with sliding mode control based course and path following for under-actuated ships. Math Probl Eng 2013; 2013(1): 1–9.

Lekkas

Fossen

. Integral LOS path following for curved paths based on a monotone cubic Hermite spline parametrization. IEEE Trans Control Syst Technol 2014; 22(6): 2287–2301.

10.

Liu

Gao

. Ship trajectory tracking control system design based on sliding model control algorithm. Pol Marit Res 2018; 25(3): 26–34.

11.

Wang

Fan

, et al. Adaptive trajectory tracking control for under-actuated unmanned surface vehicle subject to unknown dynamics and time-varing disturbances. Appl Sci (Basel) 2018; 8(4): 1–16.

12.

Liang

Wan

, et al. Three-dimensional path following of an underactuated AUV based on fuzzy backstepping sliding mode control. Int J Fuzzy Syst 2018; 20(2): 640–649.

13.

Huang

Fan

. Path following control for under-actuated surface vessel with disturbance. In: Proceedings of the 30th Chinese control and decision conference, Shenyang, China, 9–11 June 2018, pp. 3265–3269. IEEE.

14.

Pham

Nguyen

. Design an adaptive autopilot for an unmanned surface vessel. In: Proceedings of the 4th international conference on green technology and sustainable development, Ho Chi Minh City, Vietnam, 23–24 November 2018, pp. 323–328. IEEE.

15.

Zhang

. Disturbance-observer-based path following predictive controller for underactuated surface vessels. In: Proceedings of the 36th Chinese control conference, Dalian, China, 26–28 July 2017, pp. 668–673. IEEE.

16.

Tong

. Robust adaptive fuzzy backstepping output feedback tracking control for nonlinear system with dynamic uncertainties. Sci China Inform Sci 2010; 53(2): 307–324.

17.

Tong

. Observer-based adaptive fuzzy backstepping control of uncertain nonlinear pure-feedback systems. Sci China Inform Sci 2014; 57(1): 1–14.

18.

Dong

Liu

Wan

, et al. Straight-path tracking control of underactuated USV based on Takagi-Sugeno fuzzy neural network. China J Sci Instrum 2015; 36(04): 863–870.

19.

Zhang

Guan

. Concise robust adaptive path following control for underactuated ships. J Harbin Eng Univ 2014; 35(9): 1053–1059.

20.

Zheng

Negenborn

Lodewijks

. Trajectory tracking of autonomous vessels using model predictive control. IFAC Proc Vol 2014; 47(3): 24–29.

21.

Zhang

Liu

, et al. Predictive control for straight path following of under-actuated surface vessels with roll constraints. In: Proceedings of the 28th Chinese control and decision conference, Yinchuan, China, 28–30 May 2016, pp. 583–588. IEEE.

22.

Shen

Jing

. Neuron adaptive iterative sliding-mode control for path tracking of underactuated ship. J Harbin Eng Univ 2019; 40(3): 489–500.

23.

Guo

. Finite-time predictor line-of-sight–based adaptive neural network path following for unmanned surface vessels with unknown dynamics and input saturation. Int J Adv Robot Syst 2018; 15(6): 1–14.

24.

Zheng

Feroskhan

. Path following of a surface vessel with prescribed performance in the presence of input saturation and external disturbances. IEEE/ASME Trans Mechatron 2017; 22(6): 2564–2575.

25.

Qiu

Wang

Fan

, et al. Adaptive sliding mode trajectory tracking control for unmanned surface vehicle with modeling uncertainties and input saturation. Appl Sci (Basel) 2019; 9(6): 1–18.

26.

Lin

. Model predictive control—status and challenges. Acta Autom Sinica 2013; 39(3): 222–236.

27.

Young

Willems

. An approach to the linear multivariable servomechanism problem. Int J Control 1972; 15(5): 961–979.

28.

Bai

, et al. An adaptive neural network approach for ship roll stabilization via fin control. Neurocomputing 2016; 173(2016): 953–957.

29.

Wang

Zou

Wang

. Path following of a surface ship sailing in restricted waters under wind effect using robust H_∞ guaranteed cost control. Int J Nav Archit Ocean Eng 2019; 11(1): 606–623.

30.

Daidola

Graham

Chandrash

. A simulation program for vessel’s maneuvering at slow speeds. In: Proceedings of the 11th ship technology and research symposium (STAR), Portland, OR, USA, 21–23 May 1986, pp. 156–161.

31.

Al Seyab

. Differential recurrent neural network based predictive control. Comput Chem Eng 2007; 32(7): 1533–1545.

32.

Amat

Busquier

Gutierrez

, et al. On the global convergence of Chebyshev’s iterative method. J Comput Appl Math 2007; 220(1–2): 17–21.

33.

Craig

. MVT: a most valuable theorem. Berlin: Springer, 2017, pp. 151–444.

34.

Zhang

Nam-kyun

. Ship nonlinear-feedback course keeping algorithm based on MMG model driven by bipolar sigmoid function for berthing. Int J Nav Archit Ocean Eng 2017; 9(5): 525–536.