Adaptive Auto-Berthing Control of Underactuated Vessel Based on Barrier Lyapunov Function

Abstract

This paper investigates the automatic berthing problem of underactuated surface vessels in the case of uncertain dynamics and yaw rate limitation, given the importance of yaw rate control and the unmeasurable hydrodynamic parameters of the vessel at low speeds. First, we use the differential homeomorphism coordinate transformation to solve the problem of underactuation. Second, a radial basis function network (RBF) is introduced to approximate unknown nonlinear functions. Third, we apply the barrier Lyapunov function (BLF) approach to limit the yaw rate within a safe range. Fourth, we use dynamic surface control (DSC) technology and minimum learning parameters (MLP) to tackle the differential explosion problems in backstepping and computational complexity. Finally, Lyapunov stability theory proves that signals produced by the designed control scheme are bounded and effective. The simulation results show that, compared with the control scheme without BLF, the proposed method can effectively limit the yaw rate within a specific range and effectively solves the influence of the model uncertainly.

1. Introduction

Generally, vessels berth with the assistance of thrusters or tugs. However, most ships sailing at sea are only equipped with main propellers and rudder devices, while some ships are equipped with either two independent aft thrusters or one main aft thruster and a rudder without any bow or side thrusters. The insufficient force on the sway direction makes it difficult to control the vessel, especially when berthing at a low speed. Therefore, solving this problem will help save investment and construction costs, reduce the weight of the system, and increase its flexibility and operating reliability during the design process of the ship [1,2]. Researchers usually regard the underactuated vessel mentioned above as a multi-input multi-output underactuated system [3].

Solving the problem of underactuated vessel auto-berthing means stabilizing the ship at the pier by controlling its rudder and propeller. During this process, precisely controlling the vessel’s position and speed is the key to realizing safe berthing, especially the speed control. Due to the massive inertia of the ship, excessive speed margin will cause huge damage to the wharf and hull. Therefore, limiting the yaw rate within a certain range during berthing will be meaningful. Meanwhile, the underactuated vessel system will be affected by dynamic uncertainty due to low-speed ships, diving effects, and quay wall effects. Accurately estimating the dynamic uncertainty of the vessel motion system will bring the model closer to reality and make the test results more referential.

Intelligent algorithms, classic control, and a combination of the two are the most common way to tackle the problem of underactuated vessels auto-berthing. In terms of the intelligent algorithms, Artificial Neural Network (ANN) and Proximal Policy Optimization (PPO) [4] are already used in auto-berthing research. As for the former, the training data are usually obtained manually by controlling the ship to a specific pier in advance and training the weights and biases of the ANN, which could help realize auto-berthing control in a specific wharf [5]. However, this method cannot achieve vessel auto-berthing control in different ports without training data. To solve this problem, the works in [6,7] use coordinate conversion controller switching technology and a distance measurement system and controller with ship sub-routes to tackle the problems of geographic coordinate limitation, measurement accuracy, and repeated network training and improve the adaptability of ANN berthing controllers.

In terms of the classical control methods, the auto-berthing problem tends to be classified as a stabilization control problem, i.e., to design a control scheme that makes the vessel stable in a certain position and state. Direct Liapunov method [8,9] and sliding mode control [10] are the most common methods applied to realize auto berthing. In [8], an auto-berthing controller based on concise backstepping was proposed. Bu [10] designed a dynamic feedback controller and applied the recursive decomposition iterative method to tackle the underactuated vessel control problem. However, classical control theories are all based on the precise model of the controlled object. Thus, when the vessel berths at a low speed, the underactuated vessel system is affected by dynamic uncertainty.

In summary, whether intelligent algorithms are applied to implement ship berthing controllers or simple classical control methods are limited by the need for accurate mathematical models of ships, offline data training, or data consistency issues will occur, thus causing certain limitations to the engineering application of the controller. Therefore, researchers have been trying to combine the advantages of the two methods, which mainly focus on the control itself and berthing path planning. From the perspective of ship berthing control, in [11], proportional derivative (PD) control was introduced into the ANN berthing controller designation to tackle the problem of training data consistency during auto berthing. In [12], controllers are proposed that realized ship auto-berthing in a turn-around way and solved the dynamic uncertainty with RBF. From the perspective of berthing path planning, a controller based on Immune Memory-Particle Swarm Optimization (IM-PSO) is proposed to optimize the proportional integral derivative control parameters of berthing path tracking [13]. An extended dynamic window approach for the automatic berthing of underactuated surface vessels [14] could realize the automatic berthing under the influence of wind loads and obstacles. Han [15] introduced a layered artificial potential field method to the berthing trajectory planning to tackle the problem of excessive turning of the berthing trajectory. The research mentioned above mainly focuses on the steady-state characteristics of the underactuated vessel motion control system but pays less attention to the transient performance, which might cause a collision with nearby vessels as discussed in [16]. Dai et al. [17] pointed out that when a ship is operating in a narrow waterway, the force will change. The system state will be restricted, i.e., the output or state of the ship is not allowed to exceed a certain constraint distance of the reference trajectory path, which is similar to a finite time control method [18] proposed based on barrier Lyapunov function realized course keeping to improve control performance significantly. Tee, K.P., et al. [19] used barrier Lyapunov to tackle the problem of restricted location during berthing. With the full-drive system with state constraints or output constraints as the research object [20], the barrier Lyapunov function is used to reach infinity when its parameters are close to the constraint boundary, and stable control of the constraint system is realized. However, the above research does not involve the attempt to limit the speed in automatic berthing, which is very important for safe berthing.

Given the above, this paper focuses on solving the auto-berthing of underactuated surface vessels in the presence of uncertain dynamics and yaw rate limitation. The differential homeomorphism transformation approach is applied to convert the vessel to a cascade form. The RBF network estimates the unknown nonlinear functions. Furthermore, a filter based on dynamic surface control (DSC) is constructed. An auto-berthing control scheme is derived by combining minimum learning parameters technique (MLP), RBF, DSC, and applying the BLF method. The Lyapunov theory proves the stability. The significance of this paper can be summarized as follows.

(1)

A novel auto-berthing control scheme considering BLF is proposed to successfully restrict the yaw rate within a smaller range during ship berthing, which has essential safety significance for real ship berthing.

(2)

The vessel model dynamic uncertainty caused by underactuation has been fully solved using the RBF. This is compared with the research in [8] that, without considering vessel model dynamic uncertainty, showed that an adaptive control scheme based on RBF can effectively approximate the unknown term and keep the output bounded.

The rest of this paper is organized as follows. The Section 2 is the preliminary introduction. The Section 3 introduces the process of the underactuated problem of the ship motion model with uncertain dynamics. The Section 4 presents the design of the automatic berthing controller using the BLF, RBF, DSC, and MLP. The Section 5 shows the simulation results and discussion. The last part presents conclusions and prospects for further research.

2. Preliminaries

Definition 1

([21]). Barrier Lyapunov function $V\left(X\right)$ is a scalar function about system $\dot{x}=f(x,t)$ defined in open set D with origin has the following characteristics:

(1) $V\left(X\right)$ is a continuous and positive function. (2) The derivative of $V\left(X\right)$ is available on open set D. (3) When x is close to the boundary of D, the $V\left(x\right)\to \infty$ holds. (4) If, $x\left(0\right)\in D$, for $\forall t>0$ and the $V\left(x\right(t\left)\right)\le b$ holds.

For any $V\left(x\right)$, if $e\to -{\mathrm{Y}}_a$ or $e\to {\mathrm{Y}}_b,V\left(e\right)\to \infty$ holds. Further, if ${\mathrm{Y}}_a={\mathrm{Y}}_b,V\left(x\right)$ is symmetrical barrier Lyapunov function, if ${\mathrm{Y}}_a\ne {\mathrm{Y}}_b$ the $V\left(e\right)$ is asymmetrical barrier Lyapunov function, in which, ${\mathrm{Y}}_a$ and ${\mathrm{Y}}_b$ are constraint value of output error $e,r_c$ and ${\overline{\mathrm{Y}}}_c$ are the constraint value of the control systems’ output satisfies ${\mathrm{Y}}_a=y_d\left(t\right)-{\mathrm{Y}}_c,{\mathrm{Y}}_b={\overline{\mathrm{Y}}}_c-y_d\left(t\right)$.

Lemma 1

([14]). For any error variable $\left|e\right|<{\mathrm{Y}}_b$, ${\mathrm{Y}}_b\in \mathbb{R}$ is a positive constant, and the following inequality (1) holds.

$$ln\frac{{\mathrm{Y}}_b^2}{{\mathrm{Y}}_b^2-e^2}\le \frac{e^2}{{\mathrm{Y}}_b^2-e^2}$$

(1)

Lemma 2

([18]). For any constant ${\mathrm{Y}}_{b1}$, let $S_i:=\left\{e_i\in \mathbb{R}:\left|e_i\right|<{\mathrm{Y}}_{bi}\right\}\subset \mathbb{R},N:={\mathbb{R}}^I\times S_1\subset {\mathbb{R}}^{I+1}$ is open set, consider the following system (2):

$$\dot{\eta }=h(t,\eta )$$

(2)

where $\eta ={[w,e_1]}^{\tau }\in N$ is state value of the system, $h:{\mathbb{R}}_{+}\times N\to {\mathbb{R}}^{l+1}$ is piecewise continuous function about t and satisfies local Lipschitz condition, l is a positive integer. For $\left|e_i\right|\to {\mathrm{Y}}_{b1}$, the following inequality (3) holds:

$$V_i\left(e_i\right)\to \infty ,{\gamma }_1(\parallel w\parallel )\le U\left(w\right)\le {\gamma }_2(\parallel w\parallel )$$

(3)

where ${\gamma }_1$ and ${\gamma }_2$ are $K_{\infty }$ functions. $U:{\mathbb{R}}^l\to {\mathbb{R}}_{+}$ and $V_i:e_i\to R_{+},i=1,\cdots ,n$ are continuously derivable and positive trigonometric functions. Let $V\left(\eta \right)={\sum }_{i=1}^n{V}_i\left(e_i\right)+U\left(w\right)$ and $e_i\left(0\right)\in S_i$. If the following inequality holds:

$$\dot{V}=\frac{\partial V}{\partial \eta }\le -\mu V+\lambda$$

(4)

where $\eta \in N,\mu$ and λ are positive constants. Then, for $\forall t\in [0,\infty ),e_i\left(t\right)\in S_i$.

Lemma 3

([21]). For compact set $X\subset {\mathsf{\Omega }}_x\in R^n\to R$, there must be an RBF that makes the following equation true.

$$F\left(X\right)={\omega }^{*T}S\left(X\right)+\epsilon$$

(5)

where $F\left(X\right)$ is a continuous nonlinear function. ${\omega }^{*}={\left[{\omega }_1,{\omega }_2,\cdots {\omega }_l\right]}^T$ is an ideal matrix that represent the weights of the RBF, $S\left(X\right)={\left[S_1\left(x\right),S_2\left(x\right)\cdots S_l\left(x\right)\right]}^T$ is the matrix of activation function. $S_i\left(x\right)$ is the Gaussian kernel function, and $S_i\left(x\right)=exp\left[\frac{-{\left(x-c_i\right)}^T\left(x-c_i\right)}{{\omega }_i^2}\right](i=1,2\cdots l)$, $c_i=\left[c_1,c_2\cdots c_q\right]\in R^q$ is the central vector value of $S_i\left(x\right)$. ${\omega }_i\in R$ is the width of $S_i\left(x\right)$. ε is the matrix of the approximation error.

Lemma 4

([22]). For a system

$$\dot{x}=f\left(x\right),t\ge 0$$

(6)

where $f\left(0\right)=0$. If there is a scalar function $V\left(x\right)$ with continuous first derivative, $V\left(0\right)=0$, and for all non-zero points x in the state space X the following conditions are satisfied:

(1) 

$V\left(x\right)$ is positive definite;

(2) 

$\dot{V}\left(x\right)$ is negative definite;

(3) 

when $\parallel x\parallel \to \infty$, $V\left(x\right)\to \infty$.

Then, the original equilibrium state of the system is asymptotically stable in an extensive range. In the process of the controller design, the following assumptions are made for the desired value of the system.

Assumption 1

([20]). For $Y_c\left(t\right)>0$, there exist ${\underline{Y}}_0,Y_0,A_0,Y_1,\cdots ,Y_n$ satisfy $max\left\{{\underline{Y}}_0,{\overline{Y}}_0\right\}\le A_0<{\mathrm{Y}}_c$, make the desired trajectory $x_d\left(t\right)$ and its derivative satisfy $-Y_0\le x_d\left(t\right)\le {\overline{Y}}_0,\left|{\dot{x}}_d\left(t\right)\right|\le Y_1,\left|{\ddot{x}}_d\left(t\right)\right|\le Y_2$, $\left|x_d^{\left(n\right)}\left(t\right)\right|\le Y,\forall t>0$

Assumption 2

([23]). For $Y_a\left(t\right)$ and ${\Gamma }_a\left(t\right),i=1,\dots ,n$, there exists constants ${\mathrm{Y}}_a,I_a,i=1,\dots ,n$ the derivatives of $Y_a\left(t\right),\forall t>0$.

Assumption 3.

The extremal time-varying disturbance satisfies ${\tau }_{\mathrm{wi}}\le {\tau }_{wi}^{*}(i=u,v,r)$, where ${\tau }_{\mathrm{xi}}^{*}(i=u,v,r)$ is the disturbance upper bound, and all of them are normal numbers.

3. Problem Formulation

In general, the three degrees of freedom (3-DOF) underactuated marine surface vessel motion model can be described as in [24], as shown in Figure 1.

$$\left\{\begin{array}{c}\dot{x}=ucos\left(\psi \right)-vsin\left(\psi \right)\hfill \\ \dot{y}=usin\left(\psi \right)-vcos\left(\psi \right)\hfill \\ \dot{\psi }=r\hfill \\ \dot{u}=\frac{m_{22}}{m_{11}}vr-\frac{d_u}{m_{11}}v-\sum _{i=2}^3\frac{d_{ui}}{m_{11}}{\left|u\right|}^{i-1}u+\frac{{\tau }_u}{m_{11}}+\frac{{\tau }_{wu}}{m_{11}}\hfill \\ \dot{v}=-\frac{m_{11}}{m_{22}}ur-\frac{d_v}{m_{22}}v-\sum _{i=2}^3\frac{d_{vi}}{m_{22}}{\left|v\right|}^{i-1}v+\frac{{\tau }_{wv}}{m_{22}}\hfill \\ \dot{r}=\frac{m_{11}-m_{22}}{m_{11}}uv-\frac{d_r}{m_{33}}r-\sum _{i=2}^3\frac{d_{ri}}{m_{33}}{\left|r\right|}^{i-1}r+\frac{{\tau }_r}{m_{11}}+\frac{{\tau }_{wr}}{m_{11}}\hfill \end{array}\right.$$

(7)

where $x,y,\psi$-denotes the vehicles’ position and heading. $u,v,r$ denotes the surge, sway, and yaw velocities in the body-fixed frame. ${\tau }_u,{\tau }_r$ denotes the actual input vector. $d_X(X=u,v,r)$, $d_{Xi}(i=2,3)$ are hydrodynamic damping. ${\tau }_{wX}(X=u,v,r)$ are unknown time-varying disturbance.

To solve the underactuated problem of the marine surface vessel model, the following differential homeomorphism transformation is introduced [25]:

$$h=J^{\mathrm{T}}\left(\psi \right)\eta$$

(8)

where $h={\left[h_1,h_2,h_3\right]}^{\mathrm{T}},\eta ={[x,y,\psi ]}^{\mathrm{T}},J\left(\psi \right)=\left[\begin{array}{ccc}cos\left(\psi \right)& -sin\left(\psi \right)& 0\\ sin\left(\psi \right)& cos\left(\psi \right)& 0\\ 0& 0& 1\end{array}\right]$.

Taking the time derivative of Equation (8), one can get

$${\dot{h}}_1={\dot{J}}^{\mathrm{T}}\left(\psi \right)\eta +J^{\mathrm{T}}\left(\psi \right)\dot{\eta }$$

(9)

$J\left(\psi \right)$ satisfies the following properties:

$${\dot{J}}^{\mathrm{T}}\left(\psi \right)=J\left(\psi \right)S\left(r\right)$$

(10)

where $S\left(r\right)=\left[\begin{array}{ccc}0& -r& 0\\ r& 0& 0\\ 0& 0& 0\end{array}\right]$.

Furthermore, the mathematical model of underactuated marine surface vessel after differential homeomorphism transformation can be obtained as

$$\left\{\begin{array}{c}{\dot{h}}_1=u+h_{2}r\hfill \\ {\dot{h}}_2=v-h_{1}r\hfill \\ {\dot{h}}_3=r\hfill \\ \dot{u}=\frac{m_{22}}{m_{11}}vr-\frac{d_{11}}{m_{11}}u+\frac{{\tau }_u}{m_{11}}-\frac{\Delta f_u}{m_{11}}+\frac{{\tau }_{wu}}{m_{11}}\hfill \\ \dot{v}=-\frac{m_{11}}{m_{22}}ur-\frac{d_{22}}{m_{22}}v-\frac{\Delta f_v}{m_{22}}+\frac{{\tau }_{wv}}{m_{22}}\hfill \\ \dot{r}=\frac{m_{11}-m_{22}}{m_{33}}uv-\frac{d_{33}}{m_{33}}r+\frac{{\tau }_r}{m_{33}}-\frac{\Delta f_r}{m_{33}}+\frac{{\tau }_{wr}}{m_{33}}\hfill \end{array}\right.$$

(11)

where ${\sum }_{i=2}^3\frac{d_{ui}}{m_{11}}{\left|u\right|}^{i-1}u=\frac{\Delta f_u}{m_{11}}$, ${\sum }_{i=2}^3\frac{d_{vi}}{m_{22}}{\left|v\right|}^{i-1}v=\frac{\Delta f_v}{m_{22}}$, ${\sum }_{i=2}^3\frac{d_{ri}}{m_{33}}{\left|r\right|}^{i-1}r=\frac{\Delta f_r}{m_{33}}$ are unmodeled dynamics.

Let ${\tau }_1=\dot{u},{\tau }_3=\dot{r},\frac{m_{22}}{m_{11}}=\frac{1}{A},\frac{d_{22}}{m_{22}}=B$, then Equation (11) can be written as

$$\left\{\begin{array}{c}{\dot{h}}_1=u+h_{2}r\hfill \\ {\dot{h}}_2=v-h_{1}r\hfill \\ {\dot{h}}_3=r\hfill \\ \dot{u}={\tau }_1\hfill \\ \dot{v}=-Aur-Bv-\frac{\Delta f_v}{m_{22}}+\frac{{\tau }_{wv}}{m_{22}}\hfill \\ \dot{r}={\tau }_3\hfill \end{array}\right.$$

(12)

Further, convert system Equation (12) into a chain structure system by introducing the following variable substitution and feedback transformation [26]:

$$H_2=h_2+\frac{v}{B}$$

(13)

Next, ${\dot{H}}_2=-\left(h_1+\frac{1}{B}Au\right)r-\frac{\Delta f_v}{d_{22}}+\frac{{\tau }_{wv}}{d_{22}}$, let $\mu =-\left(h_1+\frac{1}{B}Au\right)$, we have

$$\left\{\begin{array}{c}{\dot{h}}_1=-\frac{B}{A}\left(\mu +h_1\right)+H_{2}r-\frac{v}{B}r\hfill \\ {\dot{H}}_2=\mu r-\frac{\Delta f_v}{d_{22}}+\frac{{\tau }_{wv}}{d_{22}}\hfill \\ {\dot{h}}_3=r\hfill \\ \dot{v}=B\left(\mu +h_1\right)r-Bv-\frac{\Delta f_v}{m_{22}}+\frac{{\tau }_{wv}}{m_{22}}\hfill \\ \dot{\mu }={\tau }_{\mu }\hfill \\ \dot{r}={\tau }_3\hfill \end{array}\right.$$

(14)

Let $z_1=h_1,z_2=H_2,z_3=h_3,z_4=v,z_5=\mu ,z_6=r,u_1={\tau }_{\mu },u_2=\dot{r}={\tau }_3$, the system Equation (14) can be described as

$$\left\{\begin{array}{c}{\dot{z}}_1=-\frac{B}{A}\left(z_1+z_5\right)+z_2{z}_6-\frac{1}{B}{z}_4{z}_6\hfill \\ {\dot{z}}_2=z_5{z}_6+\frac{1}{d_{22}}\left({\tau }_{vv}-\Delta f_v\right)\hfill \\ {\dot{z}}_3=z_6\hfill \\ {\dot{z}}_4=B\left(z_1+z_5\right)z_6-B{z}_4+\frac{1}{m_{22}}\left({\tau }_{wv}-\Delta f_v\right)\hfill \\ {\dot{z}}_5=u_1\hfill \\ {\dot{z}}_6=u_2\hfill \end{array}\right.$$

(15)

As Equation (15) is converted from Equation (11), so if the stability of the system Equation (15) and system Equation (11) are same. In other words, when $t\to \infty$, if ${lim}_{t\to \infty }{z}_i=0(1\le i\le 6)$, $(x,y,\psi ,u,v,r)$, will gradually stable at the origin.

Theorem 1

([27]). For Equation (15), when $t\to \infty$, if $z_i(i=2,5,6)$, converges to zero, we have ${lim}_{t\to \infty }{z}_i=0(i=1,4)$ Then, Equation (14) can be simplified as

$$\left\{\begin{array}{c}{\dot{z}}_2=z_5{z}_6+\frac{1}{d_{22}}\left({\tau }_{wv}-\Delta f_v\right)\hfill \\ {\dot{z}}_3=z_6\hfill \\ {\dot{z}}_5=u_1\hfill \\ {\dot{z}}_6=u_2\hfill \end{array}\right.$$

(16)

Let $x_1=z_3$, $x_2=z_6$, $x_3=z_2$, $x_4=z_5$, we have

$$\left\{\begin{array}{c}{\dot{x}}_1=x_2\hfill \\ {\dot{x}}_2=u_2\hfill \\ {\dot{x}}_3=x_4{x}_2+\frac{1}{d_{22}}\left({\tau }_{wv}-\Delta f_v\right)\hfill \\ {\dot{x}}_4=u_1\hfill \end{array}\right.$$

(17)

According to Theorem 1, if system Equation (17) gradually stabilizes, the system Equation (15) will gradually stabilizes as well. Therefore, the remain work is make system Equation (17) stable.

4. Controller Design and Stability Analysis

The control design will be divided into two subsystems: subsystem 1 and subsystem 2.

Subsystem 1:

$$\left\{\begin{array}{c}{\dot{x}}_1=x_2\hfill \\ {\dot{x}}_2=u_2\hfill \end{array}\right.$$

(18)

Subsystem 2:

$$\left\{\begin{array}{c}{\dot{x}}_3=x_4{x}_2+\frac{1}{d_{22}}\left({\tau }_{wv}-\Delta f_v\right)\hfill \\ {\dot{x}}_4=u_1\hfill \end{array}\right.$$

(19)

4.1. Controller Design

4.1.1. Control Scheme Design of Subsystem 1

Definition of error variable

$$e_1=x_1-x_{1exp},x_{1exp}=0$$

(20)

$$e_2=x_2-{\alpha }_{op}$$

(21)

where ${\alpha }_{op}$ is the virtual filtering function DSC is introduced to the control design to eliminate the differential explosion.

$$T_1{\dot{\alpha }}_{op}+{\alpha }_{op}={\alpha }_1,{\rho }_1\left(0\right)={\alpha }_1\left(0\right)$$

(22)

where ${\alpha }_1$ and ${\alpha }_{op}$ are the input and output of the filter, respectively.

The filtering error is defined as

$${\epsilon }_1={\alpha }_{op}-{\alpha }_1$$

(23)

Taking the derivative of both sides of Equation (23), one can get

$${\dot{\epsilon }}_1=-\frac{{\dot{\epsilon }}_1}{T_1}+{\Gamma }_1\left(x_1,x_{1exp},{\dot{x}}_{1exp},{\ddot{x}}_{1exp}\right)$$

(24)

where ${\Gamma }_1\left(x_1,x_{1d},{\dot{x}}_{1d},{\ddot{x}}_{1d}\right)=-{\dot{\alpha }}_1$ is continuous, bounded, which will be abbreviated as ${\Gamma }_1$ in subsequent articles and $\left|{\Gamma }_1\right|\le L_1$.

Define the Lyapunov equation as

$$V_1=\frac{1}{2}ln\frac{{\mathrm{Y}}_{b1}^2}{{\mathrm{Y}}_{b1}^2-e_1^2}+\frac{1}{2}{\epsilon }_1^2$$

(25)

where ${\mathrm{Y}}_{b1}$ is positive scalar, ${\mathrm{Y}}_{b1}={\mathrm{Y}}_{c1}-A_0$, $A_0$, is constrained boundary, $\left|x_{1exp}\right|\le A_0$.

Taking the derivative of both sides of Equation (25), one can obtain

$${\dot{V}}_1=\frac{e_1}{{\mathrm{Y}}_{b1}^2-e_1^2}\left(e_2+{\alpha }_1+{\epsilon }_1-{\dot{x}}_{1exp}\right)-\frac{{\epsilon }_1^2}{T_1}+{\epsilon }_1{\Gamma }_1$$

(26)

Design the virtual control scheme as ${\alpha }_1=-k_1{e}_1+{\dot{x}}_{1exp}$.

Define the second Lyapunov equation as

$$V_2=V_1+\frac{1}{2}ln\frac{{\mathrm{Y}}_{b2}^2}{{\mathrm{Y}}_{b2}^2-e_2^2}$$

(27)

Then,

$${\dot{V}}_2=-k_1\frac{e_1^2}{{\mathrm{Y}}_{b1}^2-e_1^2}+\frac{e_1{e}_2}{{\mathrm{Y}}_{b1}^2-e_1^2}+\frac{e_1{\epsilon }_1}{{\mathrm{Y}}_{b1}^2-e_1^2}-\frac{{\epsilon }_1^2}{T_1}+{\epsilon }_1{\Gamma }_1+\frac{e_2^2}{{\mathrm{Y}}_{b2}^2-e_2^2}\left(u_2-{\dot{\alpha }}_{op}\right)$$

(28)

According to Lemma 3, using RBF to estimate unknown nonlinear functions $-\frac{\Delta f_r}{m_{33}}$ as follows:

$$-\frac{\Delta f_r}{m_{33}}={\omega }_r^{*T}{S}_r\left(v\right)+{\epsilon }_r$$

(29)

where ${\omega }_r^{*T}$ is the ideal weight matrix, $S_r\left(v\right)$ is the activation function of RBF, then

$$\begin{array}{cc}\hfill & {\dot{V}}_2=-k_1\frac{e_1^2}{{\mathrm{Y}}_{b1}^2-e_1^2}+\frac{e_1{e}_2}{{\mathrm{Y}}_{b1}^2-e_1^2}+\frac{e_1{\epsilon }_1}{{\mathrm{Y}}_{b1}^2-e_1^2}-\frac{{\epsilon }_1^2}{T_1}+{\epsilon }_1{\Gamma }_1\hfill \\ \hfill & +\frac{e_2^2}{{\mathrm{Y}}_{b2}^2-e_2^2}\left(g+\frac{{\tau }_r}{m_{33}}+{\omega }_r^{*T}{S}_r\left(v\right)+{\epsilon }_r+\frac{{\tau }_{wr}}{m_{33}}-{\dot{\alpha }}_{op}\right)\hfill \end{array}$$

(30)

The control scheme of subsystem 1 is designed as

$$\left\{\begin{array}{c}{\tau }_r=m_{33}\left(-g+{\dot{\alpha }}_{op}-c_1{e}_2{\widehat{\delta }}_r\xi {\vartheta }_r^2\left(v\right)-\frac{{\mu }_1}{{\mu }_2}{e}_1\right)\\ \dot{\widehat{\theta }}={\mu }_2{c}_1{\vartheta }_r^2\left(v\right)e_2^2-{\beta }_r{\widehat{\delta }}_r\end{array}\right.$$

(31)

4.1.2. Control Scheme Design of Subsystem 2

Definition of error variable

$$e_3=x_3-x_{3exp},x_{3exp}=0$$

(32)

$$e_4=x_4-{\alpha }_{2op}$$

(33)

where ${\alpha }_{2op}$ is the virtual filtering function Here, we introduce DSC to the control design to eliminate the differential explosion.

$$T_2{\dot{\alpha }}_{2op}+{\alpha }_{2op}={\alpha }_2,{\alpha }_{2op}\left(0\right)={\alpha }_2\left(0\right)$$

(34)

where ${\alpha }_2$ and ${\alpha }_{2op}$ is the input and output of the filter, respectively. The filtering error is defined as

$${\epsilon }_2={\alpha }_{2op}-{\alpha }_2$$

(35)

Take the derivative both side of of Equation (35), one can get

$${\dot{\epsilon }}_2=-\frac{{\dot{\epsilon }}_2}{T_2}+{\Gamma }_2\left(x_3,x_{3exp},{\dot{x}}_{3exp},{\ddot{x}}_{3exp}\right)$$

(36)

where ${\Gamma }_2\left(x_3,x_{3exp},{\dot{x}}_{3exp},{\ddot{x}}_{3exp}\right)=-{\dot{\alpha }}_2$ is continuous, bounded, which will be abbreviated as ${\Gamma }_2$ in subsequent articles, $\left|{\Gamma }_2\right|\le L_2$.

Define the third Lyapunov equation as

$$V_3=\frac{1}{2}{e}_3^2+\frac{1}{2}{\epsilon }_2^2$$

(37)

Then, take the derivative both side of Equation (37), one can obtain

$${\dot{V}}_3=e_3\left[\left(e_4+{\epsilon }_2+{\alpha }_2\right)x_2-{\dot{x}}_{3exp}\right]+e_3\left(\frac{{\tau }_{wv}}{d_{22}}-\frac{\Delta f_v}{d_{22}}\right)-\frac{{\epsilon }_2^2}{T_2}+{\Gamma }_2{\epsilon }_2$$

(38)

Design the virtual control scheme ${\alpha }_2=-k_3{e}_3{x}_2+\frac{{\dot{x}}_{3exp}}{x_2}$ then

$${\dot{V}}_3=-k_3{e}_3^2{x}_2^2+e_3{e}_4{x}_2+e_3{x}_2{\epsilon }_2+e_3\left(\frac{{\tau }_{wv}}{d_{22}}-\frac{\Delta f_v}{d_{22}}\right)-\frac{{\epsilon }_2^2}{T_2}+{\Gamma }_2{\epsilon }_2$$

(39)

Define the forth Lyapunov equation as

$$V_4=V_3+\frac{1}{2}{e}_4^2$$

(40)

Take the derivative both side of Equation (40), one can get

$$\begin{array}{cc}\hfill & {\dot{V}}_4=k_3{e}_3^2{x}_2^2+e_3\left(\frac{{\tau }_{wv}}{d_{22}}-\frac{\Delta f_v}{d_{22}}\right)-\frac{{\epsilon }_2^2}{T_2}+{\Gamma }_2{\epsilon }_2+e_3{x}_2{\epsilon }_2\hfill \\ \hfill & +e_4\left[\left(\frac{d_{11}}{d_{22}}-1\right)u-\left(\frac{{\tau }_u}{d_{22}}-\frac{\Delta f_u}{d_{22}}+\frac{{\tau }_{wu}}{d_{22}}\right)-{\dot{\alpha }}_{2op}\right]\hfill \end{array}$$

(41)

Let $\Theta =\left(\frac{d_{11}}{d_{22}}-1\right)u-{\dot{\alpha }}_{2op}$, then

$$\begin{array}{c}\hfill {\dot{V}}_4=k_3{e}_3^2{x}_2^2+e_3\left(\frac{{\tau }_{wv}}{d_{22}}-\frac{\Delta f_v}{d_{22}}\right)-\frac{{\epsilon }_2^2}{T_2}+{\Gamma }_2{\epsilon }_2+e_3{x}_2{\epsilon }_2+e_4\left[\Theta -\frac{{\tau }_u}{d_{22}}+\frac{\Delta f_u}{d_{22}}-\frac{{\tau }_{wu}}{d_{22}}\right]\end{array}$$

(42)

According to the principle of RBF, $-\frac{\Delta f_v}{d_{22}}$ and $\frac{\Delta f_u}{d_{22}}$ can be written as

$$-\frac{\Delta f_v}{d_{22}}={\omega }_v^{*T}{S}_v\left(v\right)+{\epsilon }_v$$

(43)

$$\frac{\Delta f_u}{d_{22}}={\omega }_u^{*T}{S}_u\left(v\right)+{\epsilon }_u$$

(44)

Substituting Equations (43) and (44) into Equation (42), one can get

$$\begin{array}{cc}\hfill & {\dot{V}}_4=k_3{e}_3^2{x}_2^2+e_3\left(\frac{{\tau }_{wv}}{d_{22}}+{\omega }_v^{*T}{S}_v\left(v\right)+{\epsilon }_v\right)-\frac{{\epsilon }_2^2}{T_2}+{\Gamma }_2{\epsilon }_2\hfill \\ \hfill & +e_3{x}_2{\epsilon }_2+e_4\left[\Theta -\frac{{\tau }_u}{d_{22}}+{\omega }_u^{*T}{S}_u\left(v\right)+{\epsilon }_u-\frac{{\tau }_{wu}}{d_{22}}\right]\hfill \end{array}$$

(45)

The control scheme of subsystem 2 is designed as

$${\tau }_u=d_{22}f+d_{22}{k}_4{e}_4+d_{22}{c}_3{\widehat{\delta }}_u{\vartheta }_u^2\left(v\right)e_4$$

(46)

$${\dot{\widehat{\delta }}}_v=c_2{\vartheta }_v\left(v\right)e_3^2+{\beta }_v{\widehat{\delta }}_v$$

(47)

$${\dot{\widehat{\delta }}}_u=c_3{\vartheta }_u\left(v\right)e_4^2+{\beta }_u{\widehat{\delta }}_u$$

(48)

Overall, the following underactuated marine surface vessel automatic berthing control scheme is designed as

$$\left\{\begin{array}{c}{\tau }_u=d_{22}f+d_{22}{k}_4{e}_4+d_{22}{c}_3{\widehat{\delta }}_u{\vartheta }_u^2\left(v\right)e_4\\ {\tau }_r=m_{33}\left(-g+{\dot{\alpha }}_{op}-c_1{e}_2{\widehat{\delta }}_r{\vartheta }_r^2\left(v\right)-\frac{{\mu }_1}{{\mu }_2}{e}_1\right)\end{array}\right.$$

(49)

4.2. Stability Analysis

For the underactuated vessels steering motion model Equation (7), under Assumptions 1–3, by design the control scheme, and choose the design parameters, the vessel can finish berthing task and the proposed control protocol can sure that

• all signals of the closed-loop system are bounded,
• the constrained state variables do not exceed the constraint bounds, and
• the tracking error signal can converge to an arbitrarily small neighborhood of the origin when appropriate design parameters are chosen.

Proof. 

The Lyapunov function for the subsystem is designed as below:

$$\begin{array}{cc}\hfill & V^{*}=\frac{1}{2}ln\left(\frac{{\mathrm{Y}}_{b1}^2}{{\mathrm{Y}}_{b1}^2-e_1^2}\right)+\frac{1}{2}ln\left(\frac{{\mathrm{Y}}_{b2}^2}{{\mathrm{Y}}_{b2}^2-e_1^2}\right)+\frac{1}{2}{\epsilon }_1^2+\frac{1}{2}{\tilde{\delta }}_r^2\hfill \\ \hfill & +\frac{1}{2}{e}_3^2+\frac{1}{2}{e}_4^2+\frac{1}{2}{\tilde{\delta }}_u^2+\frac{1}{2}{\tilde{\delta }}_v^2+\frac{1}{2}{\epsilon }_2^2\hfill \end{array}$$

(50)

Taking the derivative of both sides of Equation (50), one can get

$$\begin{array}{cc}\hfill & {\dot{V}}^{*}=\frac{e_2}{{\mathrm{Y}}_{b2}^2-e_2^2}\left(g+\frac{{\tau }_r}{m_{33}}-\frac{\Delta f_r}{m_{33}}+\frac{{\tau }_{wr}}{m_{33}}-{\dot{\alpha }}_{op}\right)-k_3{e}_3^2{x}_2^2\hfill \\ \hfill & +\frac{-k_1{e}_1^2+e_1{e}_2+e_1{\epsilon }_1}{{\mathrm{Y}}_{b1}^2-e_1^2}-\frac{{\epsilon }_1^2}{T_1}+{\epsilon }_1{\Gamma }_1+e_3\left(\frac{{\tau }_{wv}}{d_{22}}-\frac{\Delta f_v}{m_{33}}+{\epsilon }_v\right)\hfill \\ \hfill & +e_4\left(\Theta -\frac{{\tau }_u}{d_{22}}+\frac{\Delta f_u}{m_{33}}+{\epsilon }_v-\frac{{\tau }_{wu}}{d_{22}}\right)-\frac{{\epsilon }_2^2}{T_2}+{\epsilon }_2{\Gamma }_2-{\tilde{\delta }}_r{\dot{\widehat{\delta }}}_r-{\tilde{\delta }}_u{\dot{\widehat{\delta }}}_u-{\tilde{\delta }}_v{\dot{\widehat{\delta }}}_v\hfill \end{array}$$

(51)

Let ${\mu }_1=\frac{1}{{\mathrm{Y}}_{b1}^2-e_1^2}$ and ${\mu }_2=\frac{1}{{\mathrm{Y}}_{b2}^2-e_2^2}$ yield

$$\begin{array}{cc}\hfill & {\dot{V}}^{*}=\frac{e_2}{{\mathrm{Y}}_{b2}^2-e_2^2}\left(g+\frac{{\tau }_r}{m_{33}}+{\omega }_r^{*T}{S}_r\left(v\right)+{\epsilon }_r+\frac{{\tau }_{wr}}{m_{33}}-{\dot{\alpha }}_{op}\right)\hfill \\ \hfill & +\frac{-k_1{e}_1^2+e_1{e}_2+e_1{\epsilon }_1}{{\mathrm{Y}}_{b1}^2-e_1^2}-\frac{{\epsilon }_1^2}{T_1}+{\epsilon }_1{\Gamma }_1-k_3{e}_3^2{x}_2^2+e_3\left(\frac{{\tau }_{wv}}{d_{22}}+{\omega }_v{}^{*T}{S}_v\left(v\right)+{\epsilon }_v\right)\hfill \\ \hfill & +e_4\left(\Theta -\frac{{\tau }_u}{d_{22}}+{\omega }_v{}^{*T}{S}_v\left(v\right)+{\epsilon }_v-\frac{{\tau }_{wu}}{d_{22}}\right)-\frac{{\epsilon }_2^2}{T_2}+{\epsilon }_2{\Gamma }_2-{\tilde{\delta }}_r{\dot{\widehat{\delta }}}_r-{\tilde{\delta }}_u{\dot{\widehat{\delta }}}_u-{\tilde{\delta }}_v{\dot{\widehat{\delta }}}_v\hfill \end{array}$$

(52)

Use MLP [28] technology to reconstruct model uncertainties and external disturbances as follows:

$$∥{\omega }_r^{*T}{S}_r\left(v\right)+{\epsilon }_r+\frac{{\tau }_{wr}}{m_{33}}∥\le ∥{\omega }_r^{*T}{S}_r\left(v\right)∥+∥{\epsilon }_r+\frac{{\tau }_{wr}}{m_{33}}∥\le {\delta }_r{\vartheta }_r\left(v\right)$$

(53)

$$∥{\omega }_v^{*T}{S}_v\left(v\right)+{\epsilon }_v+\frac{{\tau }_{wv}}{d_{22}}∥\le ∥{\omega }_v^{*T}{S}_v\left(v\right)∥+∥{\epsilon }_v+\frac{{\tau }_{wv}}{d_{22}}∥\le {\delta }_v{\vartheta }_v\left(v\right)$$

(54)

$$∥{\omega }_u^{*T}{S}_u\left(v\right)+{\epsilon }_u-\frac{{\tau }_{wu}}{m_{33}}∥\le ∥{\omega }_u^{*T}{S}_u\left(v\right)∥+∥{\epsilon }_u-\frac{{\tau }_{wu}}{m_{33}}∥\le {\delta }_u{\vartheta }_u\left(v\right)$$

(55)

where ${\delta }_r=max\left\{∥{\omega }_r^{*T}∥+∥{\epsilon }_r+\frac{{\tau }_{wr}}{m_{33}}∥\right\}$, ${\vartheta }_v\left(v\right)=∥S_v\left(v\right)∥+1$, ${\vartheta }_r\left(v\right)=∥S_r\left(v\right)∥+1$, ${\delta }_u=max\left\{∥{\omega }_u^{*T}∥+∥{\epsilon }_u+\frac{{\tau }_{wu}}{m_{33}}∥\right\}$, ${\vartheta }_u\left(v\right)=∥S_u\left(v\right)∥+1$, ${\delta }_v=max\left\{∥{\omega }_v^{*T}∥+∥{\epsilon }_v+\frac{{\tau }_{wv}}{m_{33}}∥\right\}$

Substituting Equations (55)–(57) into Equation (54), one can get

$$\begin{array}{cc}\hfill & {\dot{V}}^{*}\le -k_1{\mu }_1{e}_1^2+{\mu }_1{e}_1{e}_2+{\mu }_1{e}_1{\epsilon }_1-\frac{{\epsilon }_1^2}{T_1}+{\epsilon }_1{\Gamma }_1+{\mu }_2{e}_2\left(g+\frac{{\tau }_r}{m_{33}}-{\dot{\alpha }}_{op}\right)\hfill \\ \hfill & -{\tilde{\delta }}_r{\dot{\widehat{\delta }}}_r+{\mu }_2\left|e_2\right|{\delta }_r{\vartheta }_r\left(v\right)-k_3{e}_3^2{x}_2^2+e_3\left(\frac{{\tau }_{wv}}{d_{22}}+{\omega }_v{}^{*T}{S}_v\left(v\right)+{\epsilon }_v\right)\hfill \\ \hfill & +e_4\left(\Theta -\frac{{\tau }_u}{d_{22}}+{\omega }_u{}^{*T}{S}_u\left(v\right)+{\epsilon }_u-\frac{{\tau }_{wu}}{d_{22}}\right)-\frac{{\epsilon }_2^2}{T_2}+{\epsilon }_2{\Gamma }_2-{\tilde{\delta }}_u{\dot{\widehat{\delta }}}_u-{\tilde{\delta }}_v{\dot{\widehat{\delta }}}_v\hfill \end{array}$$

(56)

According to Young’s inequality, ${\mu }_1{e}_1{\epsilon }_1\le {\mu }_1\left(e_1^2+\frac{1}{4}{\epsilon }_1^2\right)$, $\left|{\epsilon }_1{\Gamma }_1\right|\le \frac{1}{2{\lambda }_1}{\epsilon }_1^2{\Gamma }_1^2+\frac{{\lambda }_1}{2}$, ${\mu }_2{\delta }_r{\vartheta }_r\left(v\right)\left|e_2\right|\le {\mu }_2\left(c_1{\delta }_r{\vartheta }_r^2\left(v\right)e_2^2+\frac{{\delta }_r}{4c_1}\right)$, ${\tilde{\delta }}_r{\widehat{\delta }}_r={\widehat{\delta }}_r\left({\delta }_r-{\tilde{\delta }}_r\right)\le \frac{1}{2}{\delta }^2-\frac{1}{2}{\tilde{\delta }}_r^2$, $e_3{x}_2{\epsilon }_2\le x_2\left(e_3^2+\frac{1}{4}{\epsilon }_2^2\right)$, $\left|{\epsilon }_2{\Gamma }_2\right|\le \frac{1}{2{\lambda }_2}{\epsilon }_2^2{\Gamma }_2^2+\frac{{\lambda }_2}{2}$, ${\delta }_v{\vartheta }_v\left(v\right)\left|e_3\right|\le c_2{\delta }_v{\vartheta }_v^2\left(v\right)e_3^2+\frac{{\delta }_v}{4c_2}$, ${\delta }_u{\vartheta }_u\left(v\right)\left|e_4\right|\le c_3{\delta }_u{\vartheta }_u^2\left(v\right)e_4^2+\frac{{\delta }_u}{4c_3}$, ${\tilde{\delta }}_u{\widehat{\delta }}_u={\widehat{\delta }}_u\left({\delta }_u-{\tilde{\delta }}_u\right)\le \frac{1}{2}{\delta }_u^2-\frac{1}{2}{\tilde{\delta }}_u^2$, ${\tilde{\delta }}_v{\widehat{\delta }}_v={\widehat{\delta }}_v\left({\delta }_v-{\tilde{\delta }}_v\right)\le \frac{1}{2}{\delta }_v^2-\frac{1}{2}{\tilde{\delta }}_v^2$, where ${\lambda }_1$, ${\lambda }_2$, $c_1$, $c_2$, and $c_3$ are the parameters to be designed; they are all positive constants. Then combination of Equations (31), (46), (47) and (56) gives

$$\begin{array}{cc}\hfill & V^{*}\le -{\mu }_1\left(k_1-1\right)e_1^2-k_2{\mu }_2{e}_2^2-\left(\frac{1}{T_1}-\frac{L_1^2}{2{\lambda }_1}\right){\epsilon }_1^2-\frac{1}{2}{\beta }_r{\tilde{\delta }}_r^2\hfill \\ \hfill & +\frac{1}{2}{\beta }_r{\delta }_r^2+\frac{{\lambda }_1}{2}+{\mu }_2\frac{{\delta }_r}{4c_1}-\left(k_3{x}_2^2-x_2\right)e_3^2-k_4{e}^2-\frac{{\beta }_u}{2}{\tilde{\delta }}_u^2\hfill \\ \hfill & -\frac{{\beta }_v}{2}{\tilde{\delta }}_v^2-\left(\frac{1}{T_2}-\frac{1}{2{\lambda }_2}{L}_2^2-\frac{1}{4}\right){\epsilon }_2^2+c_2{\widehat{\delta }}_v{\vartheta }_v^2\left(v\right)e_3^2+\frac{{\lambda }_2}{2}+\frac{{\delta }_v}{4c_2}+\frac{{\delta }_u}{4c_3}\hfill \end{array}$$

(57)

where ${\beta }_u$ and ${\beta }_r$ are positive constants parameters to be designed. Let $P_1=k_1-1$, $P_2=k_2$, $P_3=\frac{1}{T_1}-\frac{L_1^2}{2{\lambda }_1}$, ${\sigma }_1=\frac{1}{2}{\beta }_r{\delta }_r^2+\frac{{\lambda }_1}{2}+{\mu }_2\frac{{\delta }_r}{4c_1}$, $P_4=k_3{x}_2^2-x_2$, $P_5=k_4,P_6=\frac{1}{T_2}-\frac{1}{2{\lambda }_2}{L}_2^2-\frac{1}{4}$, ${\sigma }_2=c_2{\widehat{\delta }}_v{\vartheta }_v^2\left(v\right)e_3^2+\frac{{\lambda }_2}{2}+\frac{{\delta }_v}{4c_2}+\frac{{\delta }_u}{4c_3}$ yield

$$V^{*}\le -{\mu }_1{P}_1{e}_1^2-P_2{\mu }_2{e}_2^2-P_3{\epsilon }_1^2-\frac{1}{2}{\beta }_r{\tilde{\delta }}_r^2-P_4{e}_1^2-P_5{e}_2^2-P_6{\epsilon }_1^2-\frac{{\beta }_u}{2}{\tilde{\delta }}_u^2-\frac{{\beta }_v}{2}{\tilde{\delta }}_v^2+{\sigma }_1+{\sigma }_2$$

(58)

Let $\Phi =min\left\{2P_1,2P_2,2P_3,{\beta }_r,2P_4,2P_5,2P_6,{\delta }_u,{\delta }_v\right\},\sigma ={\sigma }_1+{\sigma }_2$, then one can obtain

$${\dot{V}}^{*}\le -\Phi V^{*}+\sigma$$

(59)

Multiply both sides of Equation (59) by $e^{\Phi t}$, one can get

$$\frac{d\left(e^{\Phi t}{V}^{*}\right)}{dt}\le \sigma e^{\Phi t}$$

(60)

Integrate both sides of Equation (60), one can obtain

$$0\le V^{*}\left(t\right)\le \frac{\sigma }{\Phi }+\left(V^{*}\left(0\right)-\frac{\sigma }{\Phi }\right)e^{-\Phi t}$$

(61)

It can be seen from Equation (61), we know that $V^{*}\left(t\right)$ is uniformly bounded, so the whole system is stable. □

5. Simulation

In the simulation studies, model vessel named “Cyber ship I” is used to verify the effectiveness of the proposed control scheme. The parameters of “Cyber ship I” can be found in [29]: $\mathrm{m}=17.6\mathrm{kg}$, ship length $\left(L\right)=1.19$ m, ${\mathrm{m}}_{11}=200\mathrm{kg}$, ${\mathrm{m}}_{22}=250\mathrm{kg}$, ${\mathrm{m}}_{33}=80\mathrm{kg}·{\mathrm{m}}^2$, $d_{11}=70\mathrm{kg}/\mathrm{s}$, $d_{22}=100\mathrm{kg}/\mathrm{s}$, $d_{33}=50\mathrm{kg}·{\mathrm{m}}^2/\mathrm{s}$. In the simulation, The design parameters are $k_1=0.13$, $k_2=0.2$, $k_3=1$, $k_4=0.66$, $c_1=0.35$, $c_3=0.0001$, ${\mathrm{Y}}_{b1}={60}^{\circ }$, ${\mathrm{Y}}_{b2}=0.15\mathrm{rad}/\mathrm{s}$, ${\mathrm{T}}_1$ = 0.01, ${\mathrm{T}}_2$ = 0.01, ${\beta }_u$$=0.0001$, ${\beta }_r=0.5.$ The initial states of the vessel is designed as $x\left(0\right)=-10\mathrm{L}$, $y\left(0\right)=-10\mathrm{L}$, $\psi \left(0\right)=\frac{\pi }{3}$, $u\left(0\right)=0.2\mathrm{m}/\mathrm{s},v\left(0\right)=0$, $r\left(0\right)=0$.

In addition, to verify the effectiveness of the control scheme designed in this paper, a comparison is made between the adaptive control considering BLF and the adaptive control without BLF. The adaptive control scheme without BLF is listed in Equation (62). The initial states of the vessel model system setting are taken to be the same as the control scheme designed in this paper, and the parameters are $k_1=2.8$, $k_2=0.08,k_3=7.5,$$k_4=1.5,c_1=0.018,c_3=0.001$. The result of the simulation is shown in Figure 2 and Figure 3.

$$\left\{\begin{array}{c}{\tau }_u=d_{22}\left[\left(\frac{d_{11}}{d_{22}}-1\right)u-+k_4{e}_4+c_3{\widehat{\delta }}_u{\vartheta }_u^2\left(v\right)e_4\right]\hfill \\ {\tau }_r=m_{33}\left(-k_2{e}_2-k_1{x}_2-\frac{m_{11}-m_{33}}{m_{33}}uv+\frac{d_{33}}{m_{33}}r-e_1-c_1{\widehat{\delta }}_r{\vartheta }_r^2\left(v\right)e_2\right.\hfill \\ {\dot{\widehat{\delta }}}_u=c_3{\vartheta }_u^2\left(v\right)e_4^2-{\beta }_u{\widehat{\delta }}_u\hfill \\ {\dot{\widehat{\delta }}}_r=c_1{\vartheta }_r^2\left(v\right)e_2^2-{\beta }_r{\widehat{\delta }}_r\hfill \end{array}\right.$$

(62)

Figure 2 shows the trajectory comparison of the vessel with control scheme Equation (49) and control scheme Equation (62). It can be seen from this graph that both schemes helped the vessel realize berthing successfully.

Figure 3a shows curves of vessel position and heading angle. The vessel under the control scheme (49) reaches the desired lateral position x = 0 at 17.27 s, and the longitudinal position reaches the desired position y = 0 at 19.49 s, while the vessel under the control scheme Equation (62) reaches desired position of lateral and longitudinal at 14.28 s and 16.09 s, respectively. In terms of heading, the vessel under control scheme Equation (62) sharply decreases to 3.9° in 5.7 s, then gradually decreases and stabilizes at 0°. In contrast, the value of the vessel heading under control scheme Equation (49) has been used for 43.59 s from 60° to 0°. It can be seen that the curve changes more smoothly.

Figure 3b displays curves change of surge velocity (u), sway velocity (v), and yaw rate (r). For the vessel with a control scheme (62), the acceleration phase of u is from 0 to 0.6 s, and the maximum speed is 3.48 m/s, then it is reduced to 0 m/s in about 15 s, while the vessel under control scheme Equation (49) reaches the maximum velocity of 2.15 m/s at approximately 2.15 s. Compared to u, v does not need to be too big in the process of vessel berthing for the vessel with control scheme Equation (62), sway speed v reaches the peak value of 0.94 m/s at around 2.4 s, decreases to 0 m/s, and stays stable at 18.87 s, while the sway speed v of the vessel is under control scheme Equation (49) reaches the maximum velocity of 0.29 m/s and stays stable at around 3.19 s. Through the comparison of the yaw rate (r), it can be seen that the vessel controlled by the control scheme (49) reaches a stable position in 26 s, and the value of r is within (−0.15 rad/s − 0.15 rad/s) during the whole process, while the yaw rate of vessel controlled by the control scheme Equation (62) exceeds 0.15 rad/s at around 0.16 s and reaches its peak value of 0.35 rad/s at around 0.6 s. Obviously, from the comparison result of the yaw rate, it can be seen that the control scheme with BLF can better limit the yaw rate within a certain range which verifies the effectiveness of our proposed control scheme.

Figure 3c illustrates the curves change of force and control force moment. It can be seen from the simulation image that the surge control force with control scheme Equation (49) changes sharply increase from −959.8 N to 869.2 N and then gradually stabilize at 0 N in 10.9 s. The yaw control force moment changes from −1089 Nm to 22.97 Nm within approximately 2.4 s and gradually decreases to 0 Nm at 11.83 s. In contrast, the simulation image that the surge control force with control scheme Equation (62) changes from 4653 N to −10 N in 6.5 s and gradually stabilizes at 0 N. The yaw control force moment reaches a peak value of 127 Nm at around 1.78 s and gradually decreases to 0 Nm at 9.7 s. The change of control force and torque is reasonable and bounded, which shows the effectiveness of the control rate designed in this paper. Figure 3d shows the change of the adaptive parameters. It can be seen from the figure that the adaptive parameters are bounded and stable.

In short, it can be seen from Figure 2 and Figure 3 that the control scheme designed in this paper can effectively limit the yaw rate in a particular range, which meets the project’s needs.

In order to further verify the robustness of the control scheme, under the condition that the initial state of the ship and the design parameters of the control scheme remain unchanged, the bounded disturbance is selected for simulation test, and the disturbance vector is selected as $0.05\ast \left[sin\right(1.5\ast \pi \ast t+\pi /6);cos(0.2\ast \pi \ast t-\pi /4);sin(0.3\ast \pi \ast t+\pi /3\left)\right]$. Figure 4 and Figure 5 present the simulation results, which are used to analyze the control performance and system robustness.

Figure 4 compares the results of trajectory with control scheme Equation (49) under disturbance (indicated by a blue line) and without disturbance (indicated by a red line). It can be found that both berthing tasks are well finished. Figure 5a shows curves of vessel position and heading angle. There is no overshoot in the two curves in the figure. Compared to the red line, the vessel’s position showed a slight deviation, and the curve of heading angle appeared fluctuations within a small bounded range between 0.13 rad/s to 0.02 rad/s after approximately 40 s.

Figure 5b demonstrates the comparison of the curves of surge velocity (u), sway velocity (v), and yaw rate (r) of the vessel controlled by control scheme Equation (49) and control scheme Equation (49) with disturbance. For control scheme Equation (49) with disturbance, the acceleration phase of u is from 0 to 1.4 s, the maximum speed is 2.14 m/s at 1.41 s, then reduced to a stable scale between −0.00013 m/s to 0.00013 m/s in 26 s, while the sway velocity (v) stabilized in a small bounded range between from −0.00025 m/s to 0.00025 m/s after approximately 26 s. By comparing the yaw rate (r), the value of the yaw rate is within [−0.15 rad/s −0.15 rad/s] during the whole process, even under the influence of disturbance. It can be seen that in the case of interference, the berthing speed curve can still be stabilized within a small range. It is very close to the curve without interference, indicating that the system has certain robustness.

Figure 5c illustrates the curves change of force and control force moment of the vessel controlled by control scheme Equation (49) without disturbance and control scheme Equation (49) with disturbance. It can be seen from the image that the surge control force, which is indicated by a blue line, sharply increases from −2789 N to 876.7 N and then gradually stabilizes in a small bounded range between −0.01 N to 0.01 N after approximately 12 s. The yaw control force moment changes from −1089 Nm to 22.73 Nm within around 2.94 s and gradually decreases the small bounded range between −0.014 Nm to 0.014 Nm. Figure 5d shows the change of the adaptive parameters. It can be seen that the (with disturbance indicated by a blue line) reaches 0.367 and keep it steady after 9.93 s, which is a little higher as compared to that without disturbance. Under the influence of disturbance, it begins to fluctuate and stabilizes in a small bounded range between 1.58 × 10${}^{-5}$ to 4.29 × 10${}^{-6}$. In short, it can be seen from Figure 4 and Figure 5 that under the disturbance, The control scheme still shows good performance when vessel berthing.

Furthermore, simulation is done between different original statuses with the same control parameters. In this study, different initial conditions of ship position and heading angle. In simulation, the initial states are taken as $x\left(0\right)=-8\mathrm{L}$, $y\left(0\right)=-5.78\mathrm{L}$, $\psi \left(0\right)=\pi /4$, $u\left(0\right)=0.2\mathrm{m}/\mathrm{s}$, $v\left(0\right)=0$, $r\left(0\right)=0$. Figure 5 and Figure 6 show the simulation results. It can be seen that the berthing task is successful with the new initial position.

Figure 7a shows that the lateral position reaches the desired position x = 0 at 15.58 s, and the longitudinal position reaches the desired position y = 0 at 14.59 s. Figure 7b illustrations that the yaw rate does not exceed the limitation. Figure 7c shows the curves change of control force and control force moment of the vessel with new initial states. It can be seen from the portrait that the surge control force with new initial states increases from −1.731 N to 609.3 N and then gradually stabilizes at 0 N in 10.1 s. The yaw control force moment changes from −816.8 Nm to 9.614 Nm within around 2.7 s and gradually decreases to 0 Nm at 9.4 s. Figure 7d shows that the adaptive control schemes are all bounded and stable, which proves the effectiveness and practicability of the control scheme designed in this paper.

We clearly illustrate the effectiveness of the proposed control rate in

Table 1. Control scheme performance comparison.

Index	Simulation Description	Result (Exceed the Limitation or Not)
Control scheme Equation (49)	No disturbance	No
Control scheme Equation (49)	Disturbed	No
Control scheme Equation (49)	Different initial position	No
Control scheme Equation (62)	Without RBF and No disturbance	Yes

.

6. Conclusions

The BLF-based auto-berthing control scheme has been proposed for an underactuated vessel system with 3-DOF considering model uncertain dynamics and yaw rate limitation. The control design applies a differential homeomorphism transformation approach to convert the vessel to a cascade form to solve the underactuated problem. An auto-berthing control scheme is derived based on the backstepping framework, combining BLF, RBF, MLP, and DSC methods. By applying RBF, the uncertain factors affected by the low speed, shallow water, and quay–wall effects are efficiently approximated. Furthermore, a DSC filter is constructed to avoid the differential explosion, and MLP is adopted to improve calculation efficiency. Compared with the auto-berthing control scheme without considering BLF, this method successfully stints the yaw rate in a relatively small range. The Lyapunov-based theoretical analysis indicates that all signals under the proposed auto berthing control scheme are bounded.

However, the proposed method does not solve the problem of the general applicability of the port, and it can only realize automatic berthing from a specific location. In addition, this study does not consider input nonlinearities that may lead to input saturation, hysteresis, dead zones, the effect of communication loads on the system, sway-yaw coupling of vessel model, and tests on actual vessels. The above aspects will be considered in the automatic berthing control design in the future.