3.1 Line-of-sight guidance law

Line-of-sight (LOS) is a conventional guidance principle, and its main idea is that if a ship is able to keep its course angle aligned with the so-called LOS angle, then the convergence to the desired position is also achieved. The LOS scheme was first applied to surface ships by Fossen et al. (Reference Fossen, Breivik and Skjetne2003). Researchers found it is useful when considering the control of underactuated ships, since it renders possible an approach of reducing the desired reference from $x_d,\,y_d,\,\psi _d$ to only one reference $\psi _d$. In this way, the path-following task $\psi \rightarrow \psi _d$ is achieved using only the one control input $\delta$.

The desired reference path is composed of a set of way-points. As shown in Figure 2, if the ship's current position is $p = [x,y]$, the LOS position $p_{{\rm los}}$ is located between the previous $p_{k-1}$ and current $p_k$ way-points. Let the ship's current horizontal position $p$ be the centre of a circle with the radius of $n$ times the ship length $L_{{\rm length}}$. This circle then intersects the current straight-line segment at two points where $p_{{\rm los}}$ is selected as the point closest to the next way-point.

Figure 2. Illustration of the guidance principles

To calculate $p_{{\rm los}} = [x_{{\rm los}}, y_{{\rm los}}]$, the following two equations need to be solved:

(3.1)\begin{align} (y_{{\rm los}} - y)^2 + (x_{{\rm los}} - x)^2 = (n L_{{\rm length}})^2 \end{align}

(3.2)\begin{align} \frac{y_{{\rm los}} - y_{k-1}}{x_{{\rm los}} - x_{k-1}} = \frac{y - y_{k-1}}{x- x_{k-1}} + \tan(\alpha_{k-1}) \end{align}

Selecting way-points in the way-point table relies on a switching algorithm. A criteria for selecting the next way-point, located at $p_{k+1} = [x_{k+1},y_{k+1}]^\top$, is for the ship to be within the range of acceptance of the current $p_k$, as shown in Figure 2. If at time $t$, the ship's current position satisfies

(3.3)\begin{equation} (x_k - x(t))^2 + (y_k - y(t))^2 \le R^2_k \end{equation}

then the next way-point will be selected from the way-point set.

With the LOS position $p_{{\rm los}}$, the LOS angle can be computed:

(3.4)\begin{equation} \psi_{{\rm los}} = \arctan2 \left(\frac{y_{{\rm los}} - y}{x_{{\rm los}} - x}\right) \end{equation}

in which the four quadrant inverse tangent function $\arctan 2(y,x)$ is used to ensure that $\psi _{{\rm los}} \in [- \pi, + \pi ]$.

The ship state vector $\boldsymbol {x} = [u,v,r,x,y,\psi ]^\top \in \mathcal {R}^6$. Then ship kinematics can be discretised as follows:

(3.5)\begin{equation} \left\{\begin{aligned} \psi_{k+1}(\boldsymbol{x}) & = \psi_k + r_k \Delta t \\ x_{k+1} (\boldsymbol{x} ) & = x_k+ (u_k\cos \psi_k - v_k\sin \psi_k ) \Delta t \\ y_{k+1}(\boldsymbol{x}) & = y_k + ( u_k \sin \psi_k + v_k\cos \psi_k ) \Delta t \end{aligned}\right. \end{equation}

where $\Delta t$ refers to the discretisation time. We can rewrite it with a linear state-space form:

(3.6)\begin{align} \begin{bmatrix} \psi_{k+1} \\ x_{k+1} \\ y_{k+1} \end{bmatrix} & = \begin{bmatrix} (\nabla_{\boldsymbol{x}} \psi_{k+1}(\boldsymbol{x}))^\top \\ (\nabla_{\boldsymbol{x}} x_{k+1} (\boldsymbol{x}))^\top \\ (\nabla_{\boldsymbol{x}} y_{k+1}(\boldsymbol{x}))^\top \end{bmatrix} \begin{bmatrix} \psi_k \\ x_k \\ y_k \end{bmatrix}+\begin{bmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \end{bmatrix}\begin{bmatrix} n_k \\ \delta_k \end{bmatrix}\nonumber\\ & \quad + \begin{bmatrix} \psi_{k+1}(\boldsymbol{x}) - (\nabla_{\boldsymbol{x}} \psi_{k+1}(\boldsymbol{x}))^\top \psi_k\\ x_{k+1}(\boldsymbol{x}) - (\nabla_{\boldsymbol{x}} x_{k+1}(\boldsymbol{x}))^\top x_k \\ y_{k+1}(\boldsymbol{x}) - (\nabla_{\boldsymbol{x}} y_{k+1}(\boldsymbol{x}))^\top y_k \end{bmatrix}. \end{align}

3.2 Curvilinear reference frame

Curvature $C(s)$ is introduced to describe the curved features of the waterway segment $s \in S$, which is a known parameter that can be drawn from the map information of a waterway. As shown in Figure 2, the control task is to ensure $e_y \rightarrow 0, e_\psi \rightarrow 0$ so that the ship can follow the reference path. The evolution of ship states $e_\psi,s,e_y$ are expressed as follows:

(3.7)\begin{equation} \left\{\begin{aligned} \dot{e}_\psi & = r - \frac{u\cos(e_\psi) - v\sin(e_\psi)}{1-C(s)e_y}C(s) \\ \dot{s} & = \frac{u\cos(e_\psi)-v\sin(e_\psi)}{1-C(s)e_y} \\ \dot{e}_y & = u\sin(e_\psi) + v\cos(e_\psi) \end{aligned}\right. \end{equation}

Therefore, the ship state vector $\boldsymbol {x} = [u,v,r,e_\psi,e_y,s]^\top \in \mathcal {R}^6$. Then, the ship kinematic model in Equation (3.7) is discretised as follows:

(3.8)\begin{equation} \left\{\begin{aligned} e_{\psi_{k+1}} (\boldsymbol{x} ) & = e_{\psi_k} + \left(r_k - \frac{u_k\cos(e_{\psi_k}) - v_k\sin(e_{\psi_k})}{1-C(s_k)e_{y_k}}C(s_k)\right) \Delta t \\ s_{k+1}(\boldsymbol{x}) & = s_k + \left(\frac{u_k\cos(e_{\psi_k})-v_k\sin(e_{\psi_k})}{1-C(s_k)e_{y_k}}\right)\Delta t \\ e_{y_{k+1}}(\boldsymbol{x}) & = e_{y_k} + \left(u_k\sin(e_{\psi_k}) + v_k\cos(e_{\psi_k})\right)\Delta t \end{aligned}\right. \end{equation}

Similar to Equation (3.6), it can be reformulate as follows:

(3.9)\begin{align} \begin{bmatrix} e_{\psi_{k+1}} \\ s_{k+1} \\ e_{y_{k+1}} \end{bmatrix} & = \begin{bmatrix} (\nabla_{\boldsymbol{x}} e_{\psi_{k+1}}(\boldsymbol{x}))^\top \\ (\nabla_{\boldsymbol{x}} s_{k+1}(\boldsymbol{x}))^\top \\ (\nabla_{\boldsymbol{x}} e_{y_{k+1}}(\boldsymbol{x}))^\top \end{bmatrix}\begin{bmatrix} e_{\psi_k} \\ s_{k} \\ e_{y_k} \end{bmatrix}+\begin{bmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \end{bmatrix}\begin{bmatrix} n_k \\ \delta_k \end{bmatrix}\nonumber\\ & \quad +\begin{bmatrix} e_{\psi_{k+1}}(\boldsymbol{x}) - (\nabla_{\boldsymbol{x}} e_{\psi_{k+1}}(\boldsymbol{x}))^\top e_{\psi_k} \\ s_{k+1}(\boldsymbol{x}) - (\nabla_{\boldsymbol{x}} s_{k+1}(\boldsymbol{x}))^\top s_k\\ e_{y_{k+1}}(\boldsymbol{x}) - (\nabla_{\boldsymbol{x}} e_{y_{k+1}}(\boldsymbol{x}))^\top e_{y_k} \end{bmatrix}. \end{align}

In this paper, the control input is chosen as $\boldsymbol {u} = [n,\delta ]^\top \in \mathcal {R}^2$, in which $n$ refers to the propeller revolution rate and $\delta$ refers to the rudder angle.

According to Section 2, the dynamic equations of ship motion variables $u,\,v$ and $r$ in Equation (2.2) include many nonlinear high-order terms. To deal with the nonlinearity, this paper proposes to learn a linear model around states $(u,v,r)$, in which regression vectors $\boldsymbol {\Gamma }^l \in \mathcal {R}^5$ ($l = \{u,v,r\}$) are introduced. Based on the values of the regression vectors, the prediction model can be reformulated as follows:

(3.10)\begin{equation} \begin{bmatrix} u_{k+1} \\ v_{k+1} \\ r_{k+1} \\ \end{bmatrix} = \begin{bmatrix} \boldsymbol{\Gamma}^u_{1:3}(\boldsymbol{x}) \\ \boldsymbol{\Gamma}^v_{1:3}(\boldsymbol{x}) \\ \boldsymbol{\Gamma}^r_{1:3}(\boldsymbol{x}) \end{bmatrix} \begin{bmatrix} u_{k} \\ v_{k} \\ r_{k} \end{bmatrix}+\begin{bmatrix} \boldsymbol{\Gamma}^u_4(\boldsymbol{x}) & 0 \\ 0 & \boldsymbol{\Gamma}^v_4(\boldsymbol{x}) \\ 0 & \boldsymbol{\Gamma}^r_4(\boldsymbol{x}) \end{bmatrix} \begin{bmatrix} n_k \\ \delta_k \end{bmatrix} + \begin{bmatrix} \boldsymbol{\Gamma}^u_5(\boldsymbol{x}) \\ \boldsymbol{\Gamma}^v_5(\boldsymbol{x}) \\ \boldsymbol{\Gamma}^r_5(\boldsymbol{x}) \end{bmatrix}, \end{equation}

where $\boldsymbol {\Gamma }^l_i(\boldsymbol {x})$ denotes the $i$th element of vector $\boldsymbol {\Gamma }^l(\boldsymbol {x})$. Linear regression methods are used to identify $\boldsymbol {\Gamma }$, which will be introduced in Section 4.

Combining Equation (3.6) or (3.9) with Equation (3.10), the 6-states ship prediction model can be restructured in the form of $\boldsymbol {x}_{k+1} = A\boldsymbol {x}_k+B\boldsymbol {u}_k+C$. As can be seen, the evolution of states $(\psi,x,y)$ or $(e_\psi,s,e_y)$ are derived from ship kinematic characteristics and independent of ship parameters. The values of $\boldsymbol {\Gamma }$ are determined via regression methods based on the data collected from the ship. In other words, the formulation in Equation (3.10) does not require pre-knowledge on ship parameters and hydrodynamic coefficients. This makes it possible to construct a prediction model for unknown ship dynamics.

4.1 Kernel-based linear regression of prediction model

To facilitate a data-driven modelling of unknown ship dynamics, this paper estimates the evolution of $u,\,v$ and $r$ via linear regression instead of identifying the ship parameters first and then linearising them. It learns a linear affine time-varying ship model around each point so as to construct the prediction model in the MPC controller.

For the unknown ship dynamics, a PID controller is designed and applied to the ship, and steering the ship from the start point $\boldsymbol {x}_S$ to the finish point of the waterway for $M_1$ loops:

(4.1)\begin{equation} \delta = K_p\boldsymbol{e} + K_d\dot{\boldsymbol{e}} + K_i\int \boldsymbol{e} \,\mathrm{d}t \end{equation}

where $\delta$ represents the rudder angle and $\boldsymbol {e}=\{e_y,e_\psi \}$. After running $M_1$ loops (iterations), the ship states and control inputs are stored in dataset $X_{\mathrm {PID}}$ with vectors $\boldsymbol {x}^i = [\boldsymbol {x}^i_0,\ldots,\boldsymbol {x}^i_{T^i}]$ and $\boldsymbol {u}^i = [\boldsymbol {u}^i_0,\ldots,\boldsymbol {u}^i_{T^i}]$, in which $T^i$ denotes the time when the ship reaches the finish point, and $\boldsymbol {x}_{T^i} \in \mathcal {X}_F$. Here, $\mathcal {X}_F$ is a set of points beyond the finish point, and $\mathcal {X}_F = \{\boldsymbol {x}\in \mathcal {R}^6:[0 \, 0 \, 0 \, 0 \,1 \, 0]\cdot \boldsymbol {x} = s \ge L\}$, in which $L$ is the length of the waterway. Here, $s\ge L$ means that the ship has passed the finish point of the waterway.

Based on dataset $X_{\mathrm {PID}}$, a set time of indices $I_{i}(\boldsymbol {x})$ of $P$ nearest neighbours of point $\boldsymbol {x}$ at iteration $i$ are identified and defined as

(4.2)\begin{equation} I_{i}(\boldsymbol{x}) =[k^{i*}_1,\ldots,k^{i*}_P] = \arg \min_{k_1,\ldots,k_P} \sum_{m=1}^{P}\lVert\boldsymbol{x} -\boldsymbol{x}^i_{k_m}\rVert^2_Q \end{equation}

in which $k_m\in \{0,\ldots,T^i\},\,m \in \{1,\ldots,P\}$, and $T^i$ refers to the ship sailing time spent in iteration $i$. Additionally, $Q$ is a scaling matrix of different variables. In other words, the set $I_i$ include the associated indices of the neighbours of each state point $\boldsymbol {x}$. Based on these data with index $I$, three kernel functions are introduced as the linear regressors, including Epanechnikov, Tri-cube and Gaussian density kernel functions:

(4.3)\begin{align} & K^{\mathrm{Epanechnikov}}(z) = \left\{\begin{array}{ll} \dfrac{3}{4}(1-z^2) & \forall |z|<1 \\ \,0 & \mathrm{otherwise} \end{array}\right. \end{align}

(4.4)\begin{align} & K^{\text{Tri-cube}}(z) = \left\{\begin{array}{ll} (1-|z|^3)^3 & \forall |z|<1 \\ \hskip4pt 0 & \mathrm{otherwise} \end{array}\right. \end{align}

(4.5)\begin{align} & K^{\mathrm{Gaussian}}(z) = \left\{\begin{array}{ll} {\rm e}^{- ({1}/{2})z^2} & \forall |z|<1 \\ 0 & \mathrm{otherwise} \end{array}\right. \end{align}

where $z = {\lVert \boldsymbol {x} -\boldsymbol {x}^i_{k_m}\rVert ^2_Q }/{h}$, and $h$ is a hyperparameter that represents the bandwidth.

To find $\boldsymbol {\Gamma }^u(\boldsymbol {x}),\boldsymbol {\Gamma }^v(\boldsymbol {x}),\boldsymbol {\Gamma }^r(\boldsymbol {x}) \in \mathcal {R}^5$ in Equation (3.10) for $u,\,v$ and $r$, the following three optimisation functions are formulated:

(4.6)\begin{equation} \left\{\begin{aligned} & J_u = \min \sum_{k,i\in I(x)} K\left(\frac{\lVert\boldsymbol{x} -\boldsymbol{x}^i_{k_m}\rVert^2_Q }{h}\right)\boldsymbol{y}^{i,u} (\boldsymbol{\Gamma}^u) \\ & J_v = \min \sum_{k,i\in I(x)} K\left(\frac{\lVert\boldsymbol{x} -\boldsymbol{x}^i_{k_m}\rVert^2_Q }{h}\right)\boldsymbol{y}^{i,v}(\boldsymbol{\Gamma}^v) \\ & J_r = \min \sum_{k,i\in I(x)} K\left(\frac{\lVert\boldsymbol{x} -\boldsymbol{x}^i_{k_m}\rVert^2_Q }{h}\right) \boldsymbol{y}^{i,r}(\boldsymbol{\Gamma}^r) \end{aligned}\right. \end{equation}

in which

(4.7)\begin{equation} \left\{\begin{aligned} & \boldsymbol{y}^{i,u}(\boldsymbol{\Gamma}^u) = \lVert u^i_{k+1} - \boldsymbol{\Gamma}^u(\boldsymbol{x})\cdot[u^i_k,v^i_k,r^i_k,n^i_k,1]^\top \rVert \\ & \boldsymbol{y}^{i,v}(\boldsymbol{\Gamma}^v) = \lVert v^i_{k+1} - \boldsymbol{\Gamma}^v(\boldsymbol{x})\cdot[u^i_k,v^i_k,r^i_k,\delta^i_k,1]^\top \rVert \\ & \boldsymbol{y}^{i,r}(\boldsymbol{\Gamma}^r) = \lVert r^i_{k+1} -\boldsymbol{\Gamma}^r(\boldsymbol{x})\cdot[u^i_k,v^i_k,r^i_k,\delta^i_k,1]^\top \rVert \end{aligned}\right. \end{equation}

Problems $J_u,\,J_v$ and $J_r$ form quadratic programming problems, which can be easily solved via an optimisation solver. Here, $\boldsymbol {x}^i_k$ represents the stored ship states data in iteration $i$ at time $k$. Solutions to problems $J_u,\,J_v$ and $J_r$ are $\boldsymbol {\Gamma }^u(\boldsymbol {x}) = \arg _{\boldsymbol {\Gamma }} \min J_u,\boldsymbol {\Gamma }^v(\boldsymbol {x})=\arg _{\boldsymbol {\Gamma }}\min J_v$ and $\boldsymbol {\Gamma }^r(\boldsymbol {x})=\arg _{\boldsymbol {\Gamma }}\min J_r$, respectively.

Based on the identified model in Equation (3.10), a linear time-varying prediction model can be constructed at time $t$ of iteration $i$ as follows:

(4.8)\begin{equation} \boldsymbol{x}^i_{k+1|t} = A^i_{k|t}\boldsymbol{x}^i_{k|t} + B^i_{k|t}\boldsymbol{u}^i_{k|t} + C^i_{k|t}, k\in \{t,\ldots,t+N-1\} \end{equation}

in which it is noted that $\boldsymbol {x}^i_{k|t} = [u^i_{k|t},v^i_{k|t},r^i_{k|t},\psi ^i_{k|t},y^i_{k|t},x^i_{k|t}]$ if it is a straight-line waterway segment and that $\boldsymbol {x}^i_{k|t} = [u^i_{k|t},v^i_{k|t},r^i_{k|t},e^i_{\psi _{k|t}},e^i_{y_{k|t}},s^i_{k|t}]$ if it is a curved waterway segment. The matrices $A^i_{k|t},\,B^i_{k|t}$ and $C^i_{k|t}$ are calculated as follows:

(4.9)\begin{equation} \begin{aligned} & A^i_{k|t} =\begin{bmatrix} \boldsymbol{\Gamma}^u_{1:3}(\bar{\boldsymbol{x}}^i_{k|t}) \quad 0 \quad 0 \quad 0 \\ \boldsymbol{\Gamma}^v_{1:3}(\bar{\boldsymbol{x}}^i_{k|t}) \quad 0 \quad 0 \quad 0 \\ \boldsymbol{\Gamma}^r_{1:3}(\bar{\boldsymbol{x}}^i_{k|t}) \quad 0 \quad 0 \quad 0 \\ (\nabla_{\boldsymbol{x}} e_{\psi_{k-1}}/\psi_{k-1}(\boldsymbol{x})\rvert_{\bar{\boldsymbol{x}}^i_{k|t}})^\top \\ (\nabla_{\boldsymbol{x}} e_{y_{k-1}}/y_{k-1}(\boldsymbol{x})\rvert_{\bar{\boldsymbol{x}}^i_{k|t}})^\top \\ (\nabla_{\boldsymbol{x}} s_{k-1}/x_{k-1}(\boldsymbol{x})\rvert_{\bar{\boldsymbol{x}}^i_{k|t}})^\top \end{bmatrix} \quad B^i_{k|t} =\begin{bmatrix} \boldsymbol{\Gamma}^u_{4}(\bar{\boldsymbol{x}}^i_{k|t}) & 0 \\ 0 & \boldsymbol{\Gamma}^v_{4}(\bar{\boldsymbol{x}}^i_{k|t}) \\ 0 & \boldsymbol{\Gamma}^r_{4}(\bar{\boldsymbol{x}}^i_{k|t}) \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \end{bmatrix} \\ & C^i_{k|t} =\begin{bmatrix} \boldsymbol{\Gamma}^u_{5}(\bar{\boldsymbol{x}}^i_{k|t}) \\ \boldsymbol{\Gamma}^v_{5}(\bar{\boldsymbol{x}}^i_{k|t}) \\ \boldsymbol{\Gamma}^r_{5}(\bar{\boldsymbol{x}}^i_{k|t}) \\ e_{\psi_{k-1}}/\psi_{k-1}(\bar{\boldsymbol{x}}^i_{k|t}) - (\nabla_{\boldsymbol{x}} e_{\psi_{k-1}}/\psi_{k-1}(\boldsymbol{x})\rvert_{\bar{\boldsymbol{x}}^i_{k|t}})^\top \bar{\boldsymbol{x}}^i_{k|t} \\ e_{y_{k-1}}/y_{k-1}(\bar{\boldsymbol{x}}^i_{k|t}) - (\nabla_{\boldsymbol{x}} e_{y_{k-1}}/y_{k-1}(\boldsymbol{x})\rvert_{\bar{\boldsymbol{x}}^i_{k|t}})^\top \bar{\boldsymbol{x}}^i_{k|t} \\ s_{k-1}/x_{k-1}(\bar{\boldsymbol{x}}^i_{k|t}) - (\nabla_{\boldsymbol{x}} s_{k-1}/x_{k-1}(\boldsymbol{x})\rvert_{\bar{\boldsymbol{x}}^i_{k|t}})^\top \bar{\boldsymbol{x}}^i_{k|t} \end{bmatrix} \end{aligned} \end{equation}

where

(4.10)\begin{equation} \bar{\boldsymbol{x}}^i_{k|t} =\left\{\begin{array}{ll} \bar{\boldsymbol{x}}^{i,*}_{k|t-1}, & k\in \{t,\ldots,t+N-1\} \\ \boldsymbol{z}^i_t, & k = t+N \end{array}\right. \end{equation}

Here, $\bar {\boldsymbol {x}}^i_{k|t} \in \bar {\boldsymbol {x}}^i_t=\{\bar {x}^i_{k|t},\ldots,\bar {x}^i_{k+N|t}\}$ refers to a set of candidate solutions that are defined using the optimal solution from the previous time step $t-1$, and $z^i_t$ represents one of the candidate terminal state of the planned ship trajectory at time $t$.

Then an MPC controller can be designed using Equation (4.8) as the prediction model. The control objective is to steer the ship from initial state $\boldsymbol {x}_S$ to the terminal set $\mathcal {X}_F$. In each sampling interval $k$, the MPC controller solves an infinite time horizon optimal control problem:

(4.11)\begin{align} & J^*_{0\rightarrow\infty}(\boldsymbol{x}_S) = \min_{\boldsymbol{u}_0,\boldsymbol{u}_1,\ldots} \sum_{k=0}^{\infty} h(\boldsymbol{x}_k,\boldsymbol{u}_k) \end{align}

(4.12)\begin{align} {\rm such\ that} \quad & \boldsymbol{x}_{k+1}= A\boldsymbol{x}_k+B\boldsymbol{u}_k\quad \forall k \ge 0 \end{align}

(4.13)\begin{align} & \boldsymbol{x}_{0} = \boldsymbol{x}_S \end{align}

(4.14)\begin{align} & \boldsymbol{x}_k \in \mathcal{X} = \{\boldsymbol{x} \in \mathcal{R}^6: F_x \le b_x \} \quad \forall k \ge 0 \end{align}

(4.15)\begin{align} & \boldsymbol{u}_k \in \mathcal{U} = \{\boldsymbol{u} \in \mathcal{R}^2: F_u \le b_u \} \quad \forall k \ge 0 \end{align}

where Equation (4.12) represents the linearised ship prediction model, Equation (4.13) defines the initial ship states, Equations (4.14) and (4.15) represent input and state constraints. In Equation (4.11), it is assumed that $h(\cdot,\cdot ) = \| \boldsymbol {x_F} - \boldsymbol {x}_k\|_Q + \|\boldsymbol {u}_k\|_R$, which refers to the stage cost. It is continuous and satisfies the following conditions in all iterations:

(4.16)\begin{equation} h(\boldsymbol{x}_F,0) = 0,\quad h(\boldsymbol{x}_k,\boldsymbol{u}_k)\succ 0\quad \forall \boldsymbol{x}_k \in \mathcal{R}^6 \setminus \{\boldsymbol{x}_F\}\quad \forall \boldsymbol{u}_k \in \mathcal{R}^2 \setminus \{0\} \end{equation}

After solving Equation (4.11), an optimal control sequence $\boldsymbol {U} = \{\boldsymbol {u}_0,\ldots,\boldsymbol {u}_{N-1}\}$ with prediction horizon $N$ at each sample time $k$, such that the resulting state sequence $\boldsymbol {X}=\{\boldsymbol {x}_0,\ldots,\boldsymbol {u}_N\}$ and the control sequence $\boldsymbol {U}$ are obtained without violating Constraints (4.12)–(4.15). This controller is then applied to the ship to run for $M_2$ loops to create another dataset $X_{\mathrm {MPC}}$. Datasets $X_{\mathrm {PID}}$ and $X_ {\mathrm {MPC}}$ form the basis for LMPC controller design.

4.2 Learning-based MPC controller design

In each iteration $i$, the controller selects $N_{\mathrm {ss}}$ points from the previous $i-N_{\mathrm {s}}$ iteration to construct a sampled safe set. Here, $N_{\mathrm {ss}}$ and $N_{\mathrm {s}}$ are the control parameters to be determined. A sampled safe set $\mathcal {SS}^i$ at iteration $i$ consists of all successful trajectories performed in the previous $i-1$ number of iterations, which is defined as

(4.17)\begin{equation} \mathcal{SS}^i = \left\{ \bigcup_{j \in M^i}\bigcup_{t=0}^{\infty} \boldsymbol{x}^j_t\right\},\quad M^i \in\left\{k \in [1,i] : \lim_{t\rightarrow \infty} \boldsymbol{x}^k_t = \boldsymbol{x}_F\right\} \end{equation}

Figure 4 gives the illustration of safe set $\mathcal {SS}^i$ in iteration $i$, which is the collection of all ship states at iteration $j$ for $j \in M^i$, where $M^i$ refers to the set of indexes $k$ associated with successful iteration $k$ for $k \le i$. It is also noted that $M^j \subseteq M^i$ and $\mathcal {SS}^j \subseteq \mathcal {SS}^i,\ \forall j \le i$. For every point in $\mathcal {SS}^i$, there exists a feasible control strategy which satisfies the state constraints and steers the state towards terminal state $\boldsymbol {x}_F$.

Figure 4. Illustration of safe set $\mathcal {SS}^i$

In addition, terminal cost and constraints are updated in each time step based on the planned ship trajectory in the previous time steps, so as to reduce computational burden. The LMPC solves a finite time constrained optimal control problem at time $k$ of iteration $i$:

(4.18)\begin{align} J^{{\rm LMPC},i}_{k:k+N}(\boldsymbol{x}^i_k) & =\min_{\boldsymbol{u}_{k|k},\ldots,\boldsymbol{u}_{k+N-1|k}} \left(\sum_{t=k}^{k+N-1}h(\boldsymbol{x}_{t|k},\boldsymbol{u}_{t|k}) + Q^{i-1}(\boldsymbol{x}_{k+N|k}) \right) \end{align}

(4.19)\begin{align} {\rm such\ that} \quad & \boldsymbol{x}_{t+1|k} = A\boldsymbol{x}_{t|k}+B\boldsymbol{u}_{t|k}\quad \forall t \in [k,\ldots,k+N-1] \end{align}

(4.20)\begin{align} & \boldsymbol{x}_{k|k} = \boldsymbol{x}^i_k \end{align}

(4.21)\begin{align} & \boldsymbol{x}_{t|k} \in \mathcal{X} ,\boldsymbol{u}_{t|k} \in \mathcal{U}\quad \forall t \in [k,\ldots,k+N-1] \end{align}

(4.22)\begin{align} & \boldsymbol{x}_{k+N|k} \in \mathcal{SS}^{i-1} \end{align}

where Equation (4.19) represents the linearised ship dynamics and Equation (4.20) defines the initial condition. The state and input constraints are defined via Equation (4.21). Equation (4.22) ensures that the terminal state reaches one of the points in the sampled safe set $\mathcal {SS}^{i-1}$ of the previous iteration $i-1$. Stage cost $h(\cdot,\cdot )$ is used to quantify the controller performance, which is the same as the stage cost in MPC formulation.

Function $Q^i(\cdot )$ is defined over $\mathcal {SS}^i$, which represents the learned minimum cost from previous iterations:

(4.23)\begin{equation} \hspace{-6pt} Q^i(\boldsymbol{x}) =\left\{\begin{array}{ll} \displaystyle\min_{j,t \in F^i(\boldsymbol{x})} J^j_{t\rightarrow \infty}(\boldsymbol{x})\sum_{k=t^*}^{\infty}h(\boldsymbol{x}^{j*}_k,\boldsymbol{u}^{j*}_k) & \mathrm{if} \ \boldsymbol{x} \in \mathcal{SS}^i \\ + \infty & \mathrm{if} \ \boldsymbol{x} \notin \mathcal{SS}^i \end{array}\right. \end{equation}

where

(4.24)\begin{equation} F^i(\boldsymbol{x}) = \{(j,t):j \in [0,i], \boldsymbol{x} = \boldsymbol{x}^j_t, \boldsymbol{x}^j_t \in \mathcal{SS}^i \} \end{equation}

For every point $\boldsymbol {x} \in \mathcal {SS}^i$, the value of $Q^i$ is determined over index pairs $(j^*,t^*)$ in the $N_{\mathrm {ss}}$ points, which is the minimum cost along the ship trajectories in $\mathcal {SS}^i$, in which

(4.25)\begin{equation} J^{j*}_{t^*\rightarrow \infty}(\boldsymbol{x}) = \sum_{k=t^*}^{\infty}h(\boldsymbol{x}^{j*}_k,\boldsymbol{u}^{j*}_k) \end{equation}

After Equation (4.18) at time $k$ of iteration $i$ is solved, solutions are obtained, including $\boldsymbol {x}^{i*}_{k:k+N|k}$ and $\boldsymbol {u}^{i*}_{k:k+N|k}$:

(4.26)\begin{equation} \left\{\begin{aligned} \boldsymbol{x}^{i*}_{k:k+N|k} & = [\boldsymbol{x}^{i*}_{k|k} , \cdot, \boldsymbol{x}^{i*}_{k+N|k} ] \\ \boldsymbol{u}^{i*}_{k:k+N|k} & = [\boldsymbol{u}^{i*}_{k|k} , \cdot, \boldsymbol{u}^{i*}_{k+N-1|k} ] \end{aligned}\right. \end{equation}

Then the first element of $\boldsymbol {u}^{i*}_{k:k+N|k}$ is applied to the ship. The finite-time optimal control problem in Equation (4.18) is solved with an primal-dual interior point method based on the Nesterov–Todd scaling at time $k+1$, based on updated state $\boldsymbol {x}_{k+1|k+1} = \boldsymbol {x}^i_{k+1}$. For more details regarding the solution algorithm, we refer readers to Andersen et al. (Reference Andersen, Roos and Terlaky2003) and Sturm (Reference Sturm2002). Algorithm 1 concludes the algorithmic steps.

Algorithm 1 The algorithmic steps of the proposed LMPC.

4.3 Asymptotic stability

To guarantee the asymptotic stability of the proposed LMPC, it is desirable to use infinite prediction and control horizons. While it is not feasible to get solutions for an infinite horizon nonlinear optimisation problem, stability of LMPC can still be guaranteed by choosing suitable safe sets and setting initial conditions. This has been studied by Rosolia and Borrelli (Reference Rosolia and Borrelli2018), and the required stability conditions are summarised as follows.

1. 1. There exists a controller that keeps the ship in the waterway when the ship has passed the finish point of the waterway. Assume that $\mathcal {X}_F$ is a control invariant, $\forall \boldsymbol {x}_k \in \mathcal {X}_F$, $\exists \boldsymbol {u}_k \in \mathcal {U}: \boldsymbol {x}_{k+1} = \boldsymbol {Ax}_k + \boldsymbol {Bu}_k \in \mathcal {X}_F$.

2. 2. Let $\mathcal {SS}^i$ be the sampled safe set at iteration $i$, $\mathcal {SS}^ 0$ is non-empty, and $\boldsymbol {x}_0 \in \mathcal {SS}^i$ is feasible and convergent to control invariant set $\mathcal {X}_F$.

3. 3. At $t=0$, $J^\mathrm {LMPC,i}_{0\rightarrow N}(\boldsymbol {x}^i_0) \le J^\mathrm {LMPC,i}_{0\rightarrow N}(\boldsymbol {x}^{i-1}_0)$ holds, $\forall i \ge 1$.

If Conditions (1)–(3) hold, then the system in a closed-loop with the controller obtained by Algorithm 1 converges to a steady-state trajectory when the number of iterations goes to infinity.

5.1 Open-loop prediction performance

To evaluate the prediction performance, the changes of ship states over time are predicted with the identified linear state-space model in Equation (4.19) with a time horizon of 20 s, starting from each point on its simulated trajectory with 1-second interval. Figure 5 shows the Root Mean Square Error (RMSE) between the predicted states and true states of the ship over prediction horizon in all iterations, with the prediction models constructed with different kernel functions labelled as LMPC-$K^{\mathrm {Epanechinikov}}$, LMPC-$K^{\text {Tri-cube}}$ and LMPC-$K^{\mathrm {Gaussian}}$, respectively. Table 1 presents the maximum, minimum and average RMSEs of different prediction models. As can be seen, the RMSEs of motion variables $u,\,v$ and $r$ are much smaller than the ship's nominal speed 1${\cdot }$4 m/s and rated speed 10$^\circ /{\rm s}$. The average RMSEs of the position variable $e_\psi$ range from 9${\cdot }$396$^\circ$ to 9${\cdot }$486$^\circ$, those for variable $s$ range from 0${\cdot }$648 to 0${\cdot }$658 m and those for variable $e_y$ range from 0${\cdot }$244 to 0${\cdot }$253 m. Among these prediction models, LMPC-$K^{\mathrm {Gaussian}}$ performs best as it leads to smaller deviations from true states.

Figure 5. Open-loop prediction performance of ship states

Table 1. Comparison of Root Mean Square Errors (RMSEs) over iterations

5.2 Closed-loop path-following control performance of MPC

To evaluate the optimisation cost and controller performance improvements over iterations, $M_3=30$ loops (iterations) have been carried out with different kernel functions. Figure 6 gives the simulated ship trajectories. In the first iteration, the simulated ship trajectories of LMPC-$K^{\mathrm {Epanechinikov}}$ and LMPC-$K^{\mathrm {Gaussian}}$ deviate from the reference path in the beginning but then keep up afterwards, while the ship trajectories of LMPC-$K^{\text {Tri-cube}}$ show larger deviations over the whole path. Figure 7 gives the simulated ship trajectories in the last iteration. It can be seen that with the increase of iterations, the simulated ship trajectories converge to better path-following performance. Among these LMPC controllers, LMPC-$K^{\mathrm {Gaussian}}$ shows relatively smaller path-following errors over iterations.

Figure 6. Simulated ship trajectories over iterations

Figure 7. Simulated ship trajectories in the last iteration

Figure 8 illustrates the simulated control inputs over iterations. It can be seen that the control inputs of different LMPC controllers varies in the beginning iterations and then converges to similar values in the last iteration.

Figure 8. Simulated control inputs over iterations

5.3 Computation costs

Figure 9 presents the changes of optimisation costs of ship trajectories over iterations, which is the sum of values of the objective function in Equation (4.18) over all sampling intervals on each trajectory. It can be seen that LMPC-$K^{\mathrm {Epanechinikov}}$ converges at an earlier iteration, while LMPC-$K^{\text {Tri-cube}}$ and LMPC-$K^{\mathrm {Gaussian}}$ show larger fluctuations. Figure 10 gives the changes of path-following errors over iterations, in which the course angle error $e_\psi$ converges and stabilises at approximately 4${\cdot }$15$^\circ$, and the error on the $y$-axis $e_y$ stabilises at approximately 0${\cdot }$18 m.

Figure 9. Changes of optimisation costs over iterations

Figure 10. Changes of average path-following errors over iterations

Figure 11 gives the average computation times for solving the optimisation problem in each sampling instant. As can be seen, the computation time ranges from 0${\cdot }$16 to 0${\cdot }$26 s, with an average of 0${\cdot }$20 s. This implies that the proposed LMPC requires far less computation time than the 1-s sampling interval, which could facilitate its implementation in practice.

Figure 11. Average computation times for solving the optimisation problem in each sampling interval over iterations