**1. Introduction**

In recent years, with the continuous development of unmanned driving, artificial intelligence, image processing, Internet space technology, etc., the whole industry has made great progress in digitalization and intelligence, and the USV is also rapidly developing towards intelligence [1]. The research on USV has also become a research hotspot. USV has a wide application prospect [2] because of its high flexibility and cost-effective characteristics, carrying specific equipment to complete a variety of dangerous specific tasks. In the military, the USV can be used for counter-reconnaissance and other intelligence collection, while in the civilian field, it can be used for water quality detection/monitoring, maritime patrol and ocean mapping [3–6]. To accomplish such missions, USVs require good motion control ability, and trajectory tracking control is the solution. Trajectory tracking is different from path tracking as for the former, the trajectory is parameterized and depends on time [7], which requires the ship to reach specified coordinates along a specific trajectory at a specified time. The under-actuation of a USV increases the control difficulty because the number of independent control inputs is less than its degrees of freedom [8] and is thus a nonlinear system [9].

**Citation:** Xu, D.; Liu, Z.; Song, J.; Zhou, X. Finite Time Trajectory Tracking with Full-State Feedback of Underactuated Unmanned Surface Vessel Based on Nonsingular Fast Terminal Sliding Mode. *J. Mar. Sci. Eng.* **2022**, *10*, 1845. https:// doi.org/10.3390/jmse10121845

Academic Editor: Alessandro Ridolfi

Received: 22 October 2022 Accepted: 21 November 2022 Published: 1 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In practice, there are many factors in the system to deal with, such as nonlinearity, parameter uncertainty, unmeasurable noises, unmeasurable speeds and external environment disturbances, which bring great challenges to the design of the control system. A neural network can approximate any nonlinear function and can be used to solve the control problems in uncertain models. In [10], the existing model uncertainties and external disturbances were treated by radial basis function neural network and disturbance observer, respectively. Similarly, neural network techniques and adaptive techniques were used to compensate for the uncertainty of the model [11]. In [12], a wavelet neural network was used to approximate completely unknown dynamic and external disturbances. In [13], an adaptive neural network was proposed to approximate uncertain nonlinear dynamics and external environmental disturbances. Combined with the backstepping method, adaptive technology was used to approximate the disturbed boundary, and a neural network was used to approximate the uncertain function to achieve trajectory tracking [14]. The difficulty of neural network radial basis function approximating arbitrary function lies in how to properly select the center points, the number of nodes and the appropriate width of radial basis function to generate hidden layers. Based on the stability theory of immersion and invariance, the asymptotic stability was guaranteed, and the adaptive law was derived to deal with the parameter uncertainty [15]. An online constructive fuzzy approximator was designed to dynamically adjust fuzzy rules to deal with time-varying uncertainties [16]. Combined with dynamic surface technology, the desired signal is smooth and bounded, and the control output realizes trajectory tracking.

In [17], the external disturbances and model uncertainty were compensated by a disturbance observer, and trajectory tracking was realized by sliding mode control. A nonlinear disturbance observer was constructed to deal with the disturbances problem [18]. For the unmeasured velocities problem, [19] constructed a hyper-distortion observer to observe the velocity. The unknown disturbance is estimated by RBF neural network combined with adaptive method. In order to ensure sufficient safety, the state of the system is subject to certain boundary constraints to perform a specific task; Li et al. [20] used the barrier Lyapunov function to constrain all states of the system to achieve trajectory tracking control. A pre-defined performance tracking control law based on quantization state was proposed to achieve tracking control without requiring system structural parameters and function approximators for a long time [21]. In the trajectory tracking control of mobile robot [22], state feedback control and disturbance feedforward compensation control are designed to solve the problems of unmeasurable speed and disturbance. Given that the disturbances are bounded, the disturbances suppression and tracking error reduction were guaranteed by finite time convergence.

Finite time control has better convergence performance and is more desirable in practical applications. It is able to complete the trajectory tracking task within finite time and has fast response characteristics [23–28]. Sliding mode control has been widely used in nonlinear uncertain systems because of its strong robustness, simple design and easy implementation. Traditional sliding mode control does not have the characteristics of finite time control. The terminal sliding mode can control the convergence rate near or away from the equilibrium point to realize finite time control because it introduces nonlinear function into sliding hyperplane design. It is widely used to solve finite time control problems [29–31]. In [32], a finite time adaptive sliding mode controller is designed, and the adaptive law is obtained by Lyapunov function. In [33], a finite time control strategy is designed based on proportional integral-differential sliding mode control to make the error converge in finite time.

The above studies deal with the problems of external environmental disturbances by designing disturbance observers to compensate for the disturbance or neural network approximation and using state observers for velocity estimation. In this study, perturbation and velocity unmeasurable problems are considered and dealt with together. NESO and NFTSM are combined to realize finite time underactuated USV trajectory tracking control. The overall structure block diagram of this study is shown in Figure 1. The underactuated

USV mathematical model is taken as the research plant, and then NESO and NFTSM are designed and verified by numerical simulation. The main contributions are as follows:

**Figure 1.** Structure diagram of research.

A NESO is designed to achieve full-state feedback, which eliminates the need for a separate disturbance observer to compensate the disturbance, thus reducing the design difficulty. Uncertainty of model parameters and external disturbances are regarded as complex disturbances and expanded as part of the state variables of the system, and the observation values of complex disturbances feedback to the controller to simplify the design process. The observed values of velocity and position are obtained and feedback to the system, which is particularly useful in the case of sensor malfunction.

NFTSM can realize the convergence of tracking error in a shorter time, which is used to solve the problem of slow convergence of trajectory tracking error. According to Lyapunov stability theory, the corresponding surge force and yaw moment are derived to ensure the stability of the system. The numerical simulation by MATLAB verifies that the controller has fast response characteristics, and the chattering phenomenon of the controller is also reduced.

#### **2. Definition**

For first-order nonlinear systems, NTSM and NFTSM were defined [34] as follows:

$$\begin{array}{ll} \sigma(t) = \mathbf{x} + k \text{sign}^a \dot{\mathbf{x}} = \mathbf{0}, k > \mathbf{0}, \ 1 < a < 2, \\\sigma(t) = \mathbf{x} + k\_1 \text{sign}^{a\_1} \mathbf{x} + k\_2 \text{sign}^a \dot{\mathbf{x}} = \mathbf{0}, \ k\_1 > 0, k\_2 > 0, 1 < a < 2, a\_1 > a \end{array} \tag{1}$$

where *<sup>x</sup>* <sup>∈</sup> *<sup>R</sup>*, sign*ax* :<sup>=</sup> sign*<sup>x</sup>* · |*x*<sup>|</sup> *<sup>a</sup>*, *k*, *k*<sup>1</sup> and *k*<sup>2</sup> are coefficients greater than 0, *a* = *p*/*q*, *p* and *q* are positive integers *p* > 0, *q* > 0.

#### **3. Underactuated USV Modeling and Problem Formulation**

*Assumption 1: Only three degrees of freedom of the underactuated USV are considered, namely, surge, sway and yaw.*

*Assumption 2: The underactuated USV has a uniform mass distribution and a left-right symmetry.*

*Assumption 3: The (Center of Gravity) CoG and Center of Buoyancy (COB) of the underactuated USV are located on the z-axis of the body-fixed coordinate system.*

On the premise of the above assumptions, the kinematic and dynamic equations of USV can be obtained as follows:

$$\begin{aligned} \dot{\boldsymbol{\eta}} &= \mathbf{R}(\boldsymbol{\psi})\mathbf{v} \\ \mathbf{M}\dot{\mathbf{v}} &= -\mathbf{C}(\mathbf{v})\mathbf{v} - \mathbf{D}(\mathbf{v})\mathbf{v} + \boldsymbol{\pi} + \boldsymbol{\pi}\_{\mathrm{E}} \end{aligned} \tag{2}$$

where η = [*x*, *y*, *ψ*] *<sup>T</sup>* denotes the actual positions and heading angle of underactuated USV in the Earth-fixed coordinate also called North-East-Down coordinate; **v** = [*u*, *v*,*r*] *<sup>T</sup>* denotes the actual velocities and angle velocity; **R**(*ψ*) is the transformation matrix between the Earth-fixed coordinate system and the Body-fixed coordinate system, as shown in Figure 2. **M** is inertial matrix; **C**(**v**) is the Coriolis and centripetal matrix; **D**(**v**) is the hydrodynamic damping parameters matrix; τ = [τ*u*, 0, τ*r*] *<sup>T</sup>* denotes the surge force, sway force and yaw moment, where sway force is 0; τ<sup>E</sup> = [τEu, τEv, τEr] is used to represent the disturbances caused by the wind, wave and current environment of the USV. The above parameter matrices satisfy the following properties:

**Figure 2.** Earth-fixed *O* − *XoYo* and Body-fixed *A* − *XY* coordinate frames of a USV.

**Property 1:** The rotation matrix is orthogonal to satisfy **R**(*ψ*) = 1 and

$$\begin{array}{c} \mathbf{R}^T(\boldsymbol{\psi}) = \mathbf{R}^{-1}(\boldsymbol{\psi})\\ \dot{\mathbf{R}}(\boldsymbol{\psi}) = \mathbf{R}(\boldsymbol{\psi})\mathbf{S}(\boldsymbol{r})\\ \mathbf{R}(\boldsymbol{\psi})\mathbf{S}(\boldsymbol{r})\mathbf{R}^T(\boldsymbol{\psi}) = \mathbf{R}(\boldsymbol{\psi})^T\mathbf{S}(\boldsymbol{r})\mathbf{R}(\boldsymbol{\psi}) \end{array} \tag{3}$$

where

$$\mathbf{R}(\boldsymbol{\psi}) = \begin{bmatrix} \cos(\boldsymbol{\psi}) & -\sin(\boldsymbol{\psi}) & 0\\ \sin(\boldsymbol{\psi}) & \cos(\boldsymbol{\psi}) & 0\\ 0 & 0 & 1 \end{bmatrix}, \mathbf{S}(\boldsymbol{r}) = \begin{bmatrix} 0 & -\boldsymbol{r} & 0\\ \boldsymbol{r} & 0 & 0\\ 0 & 0 & 0 \end{bmatrix} \tag{4}$$

**Property 2:** The symmetric positive definite of the inertia matrix satisfies:

$$\mathbf{M} = \mathbf{M}^T > 0\tag{5}$$

**Property 3:** The Coriolis and centripetal force matrices are obliquely symmetric and satisfy:

$$\mathbf{C}(\mathbf{v}) = -\mathbf{C}^T(\mathbf{v}), \forall \mathbf{v} \in \mathbb{R}^3. \tag{6}$$

**Property 4:** Hydrodynamic damping matrix is an asymmetric and strictly positive definite real matrix, which satisfies:

$$\mathbf{D}(\mathbf{v}) > 0, \forall \mathbf{v} \in \mathbb{R}^3. \tag{7}$$

The expressions of matrix **M**, **C**(**v**), and **D**(**v**) are, respectively:

$$\mathbf{M} = \begin{bmatrix} m\_{11} & 0 & 0 \\ 0 & m\_{22} & 0 \\ 0 & 0 & m\_{33} \end{bmatrix} = \begin{bmatrix} m - X\_{\dot{u}} & 0 & 0 \\ 0 & m - Y\_{\dot{v}} & 0 \\ 0 & 0 & I\_z - N\_{\dot{r}} \end{bmatrix} \tag{8}$$

$$\mathbf{C}(\mathbf{v}) = \begin{bmatrix} 0 & 0 & -m\_{22}v \\ 0 & 0 & m\_{11}u \\ m\_{22}v & -m\_{11}u & 0 \end{bmatrix} \tag{9}$$

$$\mathbf{D}(\mathbf{v}) = \begin{bmatrix} d\_{11} & 0 & 0 \\ 0 & d\_{22} & 0 \\ 0 & 0 & d\_{33} \end{bmatrix} = \begin{bmatrix} -X\_u & 0 & 0 \\ 0 & -Y\_v & 0 \\ 0 & 0 & -N\_r \end{bmatrix} \tag{10}$$

Define auxiliary quantity ω = **R**(*ψ*)**v**, there are:

$$\begin{aligned} \dot{\boldsymbol{\omega}} &= \ddot{\boldsymbol{\eta}} \\ \dot{\boldsymbol{\eta}} &= \dot{\mathbf{R}}(\boldsymbol{\psi}) \mathbf{v} + \mathbf{R}(\boldsymbol{\psi}) \dot{\mathbf{v}} \\ \dot{\boldsymbol{\eta}} &= f(\boldsymbol{\eta}, \mathbf{v}, \boldsymbol{\tau}\_{\mathbb{E}}) + \mathbf{B}(\boldsymbol{\psi}) \boldsymbol{\tau} \end{aligned} \tag{11}$$

where *<sup>f</sup>*(η, **<sup>v</sup>**, <sup>τ</sup>E) <sup>=</sup> **<sup>R</sup>**(*ψ*)**M**−1[τ<sup>E</sup> <sup>−</sup> **<sup>C</sup>**(**v**)**<sup>v</sup>** <sup>−</sup> **<sup>D</sup>**(**v**)**v**] <sup>+</sup> . **R**(*ψ*)**v**, **B**(*ψ*) = **R**(*ψ*)**M**<sup>−</sup>1.

Therefore, the underactuated USV mathematical model can be written in the following form: .

$$\begin{aligned} \dot{\boldsymbol{\eta}} &= \boldsymbol{\omega} \\ \dot{\boldsymbol{\omega}} &= \boldsymbol{f}(\boldsymbol{\eta}, \mathbf{v}, \boldsymbol{\tau}\_{\mathrm{E}}) + \mathbf{B}(\boldsymbol{\psi})\boldsymbol{\tau} \end{aligned} \tag{12}$$

In this study, a reference trajectory is generated by providing a predefined input to the virtual underactuated USV, as shown in Figure 3.

**Figure 3.** Virtual USV generation reference estimation schematic.

Similarly for virtual underactuated USV:

$$
\begin{array}{l}
\dot{\eta}\_{\rm d} & = \mathbf{R}\_{\rm d}(\psi\_{\rm d}) \mathbf{v}\_{\rm d} \\
\mathbf{M}\_{\rm d} \dot{\mathbf{v}}\_{\rm d} & = -\mathbf{C}\_{\rm d}(\mathbf{v}\_{\rm d}) \mathbf{v}\_{\rm d} - \mathbf{D}\_{\rm d}(\mathbf{v}) \mathbf{v}\_{\rm d} + \mathbf{r}\_{\rm d}
\end{array} \tag{13}
$$

Defining an auxiliary quantity ω**<sup>d</sup>** = **Rd**(*ψ*d)**vd**, there are:

$$\begin{aligned} \dot{\boldsymbol{\omega}}\_{\rm d} &= \ddot{\boldsymbol{\eta}}\_{\rm d} \\ \dot{\boldsymbol{\tau}} &= \dot{\mathbf{R}}\_{\rm d}(\boldsymbol{\psi}\_{\rm d}) \mathbf{v}\_{\rm d} + \mathbf{R}\_{\rm d}(\boldsymbol{\psi}\_{\rm d}) \dot{\mathbf{v}}\_{\rm d} \\ &= f\_{\rm d}(\boldsymbol{\eta}\_{\rm d}, \mathbf{v}\_{\rm d}) + \mathbf{B}\_{\rm d}(\boldsymbol{\psi}\_{\rm d}) \mathbf{r}\_{\rm d} \end{aligned} \tag{14}$$

$$\begin{array}{llllll}\hline \textbf{where} & f\_{\mathrm{d}}(\boldsymbol{\eta}\_{\mathrm{d}}, \mathbf{v}\_{\mathrm{d}}) & = & \mathbf{R}\_{\mathrm{d}}(\boldsymbol{\psi}\_{\mathrm{d}}) \mathbf{M}\_{\mathrm{d}}^{-1} [-\mathbf{C}\_{\mathrm{d}}(\mathbf{v}\_{\mathrm{d}}) \mathbf{v}\_{\mathrm{d}} - \mathbf{D}\_{\mathrm{d}}(\mathbf{v}) \mathbf{v}\_{\mathrm{d}}] + & \dot{\mathbf{R}}\_{\mathrm{d}}(\boldsymbol{\psi}\_{\mathrm{d}}) \mathbf{v}\_{\mathrm{d}} & \mathbf{B}\_{\mathrm{d}}(\boldsymbol{\psi}\_{\mathrm{d}}) & = & \mathbf{R}\_{\mathrm{d}}(\boldsymbol{\psi}\_{\mathrm{d}}) \mathbf{M}\_{\mathrm{d}}^{-1} \\ \mathbf{R}\_{\mathrm{d}}(\boldsymbol{\psi}\_{\mathrm{d}}) \mathbf{M}\_{\mathrm{d}}^{-1} & & & & & & \\ \hline \end{array}$$

According to Equations (11) and (14), an error model can be obtained.

$$\begin{array}{ll} \dot{\boldsymbol{\eta}}\_{\varepsilon} & = \dot{\boldsymbol{\eta}} - \dot{\boldsymbol{\eta}}\_{\mathrm{d}} \\ \ddot{\boldsymbol{\eta}}\_{\varepsilon} & = \ddot{\boldsymbol{\eta}} - \ddot{\boldsymbol{\eta}}\_{\mathrm{d}} \\ & = f(\boldsymbol{\eta}, \mathbf{v}, \boldsymbol{\tau}\_{\mathrm{E}}) + \mathbf{B}(\boldsymbol{\psi})\boldsymbol{\tau} - f\_{\mathrm{d}}(\boldsymbol{\eta}\_{\mathrm{d}}, \mathbf{v}\_{\mathrm{d}}) - \mathbf{B}\_{\mathrm{d}}(\boldsymbol{\psi}\_{\mathrm{d}})\boldsymbol{\tau}\_{\mathrm{d}} \end{array} \tag{15}$$

#### **4. Design of NESO**

In the actual navigation process, it is inevitable that the USV is subjected to the disturbances of the external environment: wind, wave and current. At the same time, considering the parameters perturbation, that is, the system inertia matrix, Coriolis centripetal force matrix and hydrodynamic damping matrix, this part of the influence will be regarded as whole complex disturbances, and it will be expanded into part of the system states to deal with together. Further considering the situation that the system states cannot be measured, NESO will be used to observe the system states, that is, the velocities, positions and composite disturbances are feedback to the underactuated USV system.

The following NESO can be designed for the underactuated USV model (12).

$$\begin{cases} \mathbf{e}\_1 = \mathbf{z}\_1 - \mathbf{x}\_1, \mathbf{e}\_2 = \mathbf{z}\_2 - \mathbf{x}\_2, \mathbf{e}\_3 = \mathbf{z}\_3 - \mathbf{x}\_3\\ \dot{\mathbf{z}}\_1 = \mathbf{z}\_2 + \lambda\_1 \mathbf{e}\_1\\ \dot{\mathbf{z}}\_2 = \mathbf{z}\_3 + \mathbf{B} \mathbf{\tau} + \lambda\_2 \text{fal}(\mathbf{e}\_1, \beta\_1, \delta\_1)\\ \dot{\mathbf{z}}\_3 = \lambda\_3 \text{fal}(\mathbf{e}\_1, \beta\_2, \delta\_2) \end{cases} \tag{16}$$

where **<sup>x</sup>**<sup>1</sup> <sup>=</sup> <sup>η</sup>, **<sup>x</sup>**<sup>2</sup> <sup>=</sup> . η = ω, the derivative of complex disturbances **x**<sup>3</sup> = *f*(η, **v**, τE) is unknown but bounded, that is, *f*(η, **v**, τE) <= *f* > 0, **z**<sup>1</sup> is position observed value of η, **<sup>z</sup>**<sup>2</sup> is the observed value of **.** η, **z**<sup>3</sup> is complex disturbances observed value of *f*(η, **v**, τE); λ1, λ2, λ<sup>3</sup> are gain matrix of NESO, **e***<sup>i</sup>* = ' *ei*,*j* (*T* (*i* = 1, 2, 3, *j* = 1, 2, 3) are the corresponding approximation errors. The function of fal(•) expression is as follow:

$$\text{fal}(\boldsymbol{x}, \boldsymbol{\beta}, \boldsymbol{\delta}) = \begin{cases} |\boldsymbol{x}|^{\beta} \text{sign}(\boldsymbol{x}) & \boldsymbol{\varepsilon} \text{ } |\boldsymbol{x}| > \boldsymbol{\delta} \\ \frac{\boldsymbol{s}}{\delta^{1-\beta}} & \boldsymbol{s} \text{ } |\boldsymbol{x}| \le \boldsymbol{\delta} \end{cases} \tag{17}$$

where *x* is independent variable of the function. *β* ∈ (0, 1), *δ* an arbitrary small positive number.

Combining Equations (12) and (16) to calculate the derivatives of **e**1, **e**2, **e**3, the error system of NESO is as follows:

$$\begin{cases}
\dot{\mathbf{e}}\_1 = \mathbf{e}\_2 - \lambda\_1 \mathbf{e}\_1 \\
\dot{\mathbf{e}}\_2 = \mathbf{e}\_3 - \lambda\_2 \mathbf{e}\_1 \\
\dot{\mathbf{e}}\_3 = f(\boldsymbol{\eta}\_\prime \mathbf{v}\_\prime \, \mathbf{r}\_\mathrm{E}) - \lambda\_3 \mathbf{e}\_1
\end{cases} \tag{18}$$

The error state equation of NESO can be expressed as:

$$
\dot{\mathbf{e}} = \overline{\mathbf{A}} \mathbf{e} + \overline{\mathbf{B}} \dot{f}(\eta\_{\nu} \mathbf{v}\_{\nu} \mathbf{\tau}\_{\mathrm{E}}) \tag{19}
$$

where

$$
\overline{\mathbf{A}} = \begin{bmatrix} -\lambda\_1 & \mathbf{I}\_{3\times3} & \mathbf{0}\_{3\times3} \\ -\lambda\_2 & \mathbf{0}\_{3\times3} & \mathbf{I}\_{3\times3} \\ -\lambda\_3 & \mathbf{0}\_{3\times3} & \mathbf{0}\_{3\times3} \end{bmatrix}, \overline{\mathbf{B}} = \begin{bmatrix} \mathbf{0}\_{3\times3} \\ \mathbf{0}\_{3\times3} \\ \mathbf{I}\_{3\times3} \end{bmatrix} \tag{20}
$$

The eigenequation of matrix **A** is

$$\left|\mathbf{sI} - \overline{\mathbf{A}}\right| = \begin{bmatrix} \mathbf{s} + \lambda\_1 & -\mathbf{I}\_{3 \times 3} & \mathbf{0}\_{3 \times 3} \\ \lambda\_2 & \mathbf{s} & -\mathbf{I}\_{3 \times 3} \\ \lambda\_3 & \mathbf{0}\_{3 \times 3} & \mathbf{s} \end{bmatrix} \tag{21}$$

The characteristic polynomial of Equation (19) is

$$\mathbf{s}\_i^3 + \mathfrak{a}\_1 \mathbf{s}\_i^2 + \mathfrak{a}\_2 \mathbf{s}\_i^1 + \mathfrak{a}\_3 = 0 \tag{22}$$

where **s***<sup>i</sup>* = [**s**1*i*, **s**2*i*, **s**3*i*] *<sup>T</sup>*, *<sup>i</sup>* = 1, 2, 3. Appropriate choice of parameter matrixes, <sup>λ</sup>*i*, *<sup>i</sup>* = 1, 2, 3, makes **A** satisfy the Hurwitz stability condition. By Lyapunov's second method, for arbitrary given positive definite matrix **Q** and **P**, the following Lyapunov equation is satisfied.

$$
\overline{A}^T P + P\overline{A} + Q = 0 \tag{23}
$$

Define the Lyapunov function regarding NESO as

$$V = \mathfrak{e}^T \mathbf{P} \mathfrak{e} \tag{24}$$

The differential Equation (24) can be obtained.

$$\begin{split} \dot{V} &= \dot{e}^T P e + e^T P \dot{e} \\ &= e^T \overline{A}^T P e + \left[ \overline{B} f(\mathfrak{n}, \mathfrak{v}, \mathfrak{r}\_{\mathbb{E}}) \right]^T P e + e^T P \overline{A} e + e^T P \overline{B} f(\mathfrak{n}, \mathfrak{v}, \mathfrak{r}\_{\mathbb{E}}) \\ &= e^T \left( \overline{A}^T P + P \overline{A} \right) e + 2e^T P \overline{B} f(\mathfrak{n}, \mathfrak{v}, \mathfrak{r}\_{\mathbb{E}}) \\ &\leq e^T (Q) e + 2 \| e \| \cdot \| P \overline{B} \| \cdot \| f(\mathfrak{n}, \mathfrak{v}, \mathfrak{r}\_{\mathbb{E}}) \| \\ &\leq \lambda\_{\min}(Q) \| e \|^2 + 2 \overline{f} \| e \| \cdot \| P \overline{B} \| \end{split} \tag{25}$$

When . *V* ≤ 0, it can be proved that NESO is convergent, and the convergence conditions are:

$$||e|| \le \frac{2f||PB||}{\lambda\_{\text{min}}(Q)}\tag{26}$$

From Equation (16), the observed value **^** η = **z**<sup>1</sup> of position, the observed value **v**ˆ = **R**−1(*ψ*)**z**<sup>2</sup> of velocity, and the observed value ˆ *f*(η, **v**, τE) = **z**<sup>3</sup> of complex disturbance can be obtained. Using the observations obtained by NESO to replace the position states and velocity states corresponding to the actual USV, and the complex disturbances composed of unknown internal disturbances, namely model parameter disturbances and external unknown wind, wave and current environment disturbances, can simplify the design of the controller and reduce the complex operation process.

#### **5. Controller Design**

The feedback value of the underactuated USV system states obtained by NESO has been mentioned before. After redefining the error model, the sliding mode surface is designed for controller design. The control scheme structure diagram is shown in Figure 4.

**Figure 4.** Control scheme structure diagram.

According to the design idea, the error model is redefined by the error model (15) and NESO (16). **.**

$$\begin{array}{rcl} \dot{\boldsymbol{\eta}}\_{\varepsilon} & = \stackrel{\scriptstyle \hat{\mathbf{n}}}{\ddot{\boldsymbol{\eta}}} - \dot{\boldsymbol{\eta}}\_{\mathrm{d}} \\ \ddot{\boldsymbol{\eta}}\_{\varepsilon} & = \stackrel{\scriptstyle \hat{\mathbf{n}}}{\ddot{\boldsymbol{\eta}}} - \ddot{\boldsymbol{\eta}}\_{\mathrm{d}} \\ & = f\Big(\stackrel{\scriptstyle \hat{\mathbf{n}}}{\ddot{\boldsymbol{\eta}}}, \dot{\boldsymbol{\eta}}, \stackrel{\scriptstyle \hat{\mathbf{r}}}{\ddot{\boldsymbol{\tau}}}\Big) + \mathbf{B}\big(\dot{\boldsymbol{\psi}}\big)\mathbf{r} - f\_{\mathrm{d}}(\boldsymbol{\eta}\_{\mathrm{d}}, \mathbf{v}\_{\mathrm{d}}) - \mathbf{B}\_{\mathrm{d}}(\boldsymbol{\psi}\_{\mathrm{d}})\mathbf{r}\_{\mathrm{d}} \end{array} \tag{27}$$

From the definition of NFTSM, the following sliding surface can be obtained:

$$s = \mathfrak{n}\_{\varepsilon} + k\_1 \text{sign}^{a\_1} \mathfrak{n}\_{\varepsilon} + k\_2 \text{sign}^a \dot{\mathfrak{n}}\_{\varepsilon} \tag{28}$$

From Equations (27) and (28), the differential equation of the sliding mode surface can be obtained:

$$\dot{s} = \dot{\mathfrak{n}}\_{\varepsilon} + a\_1 k\_1 \text{diag}\left(\text{sign}^{a\_1 - 1} \mathfrak{r}\_{\varepsilon}\right) \dot{\mathfrak{n}}\_{\varepsilon} + a k\_2 \text{diag}\left(\text{sign}^{a - 1} \dot{\mathfrak{n}}\_{\varepsilon}\right) \ddot{\mathfrak{n}}\_{\varepsilon} \tag{29}$$

Since **..** η*<sup>e</sup>* contains the required control τ, the Lyapunov function about the sliding surface is constructed, and the desired controller output τ can be obtained by selecting the appropriate τ to satisfy the Lyapunov stability condition.

Define the following Lyapunov function,

$$V\_1 = \frac{1}{2} \mathbf{s}^T \mathbf{s} \tag{30}$$

and obtain the differential equation:

$$\begin{aligned} \dot{V}\_1 &= \mathbf{s}^T \dot{\mathbf{s}} \\ &= \mathbf{s}^T \Big[ \dot{\boldsymbol{\eta}}\_\varepsilon + a\_1 k\_1 \text{diag} \Big( \text{sign}^{a\_1 - 1} \mathbf{n}\_\varepsilon \Big) \dot{\boldsymbol{\eta}}\_\varepsilon + a k\_2 \text{diag} \Big( \text{sign}^{a - 1} \dot{\boldsymbol{\eta}}\_\varepsilon \Big) \ddot{\boldsymbol{\eta}}\_\varepsilon \Big] \end{aligned} \tag{31}$$

If the underactuated USV input is selected as:

$$\begin{split} \mathbf{r} &= \mathbf{M} \mathbf{R}^{-1} (\boldsymbol{\psi}) \left\{ -\boldsymbol{f} \left( \hat{\mathbf{\eta}}\_{\text{\prime}} \boldsymbol{\psi}, \hat{\mathbf{\tau}}\_{\text{E}} \right) + \mathbf{R}\_{\text{d}} (\boldsymbol{\psi}\_{\text{d}}) \mathbf{M} \mathbf{d}^{-1} \mathbf{r}\_{\text{d}} + f\_{\text{d}} (\mathbf{\eta}\_{\text{d}\prime} \mathbf{v}\_{\text{d}}) \\ &- \frac{1}{a k\_{2}} \left[ a k\_{1} \text{diag} \left( \text{sign}^{a\_{1}-1} \boldsymbol{\eta}\_{\text{c}} \right) + \mathbf{I} \right] \text{sign}^{2-a} \dot{\boldsymbol{\eta}}\_{\text{c}} - k\_{3} \text{s} - k\_{4} \text{sign} \mathbf{s} \right\} \end{split} \tag{32}$$

Substituting Equation (32) into Equation (31) can be obtained:

$$\begin{split} \dot{V}\_{1} &= \mathbf{s}^{T}\dot{\mathbf{s}} \\ &= \mathbf{s}^{T} \Big[\dot{\mathbf{n}}\_{\varepsilon} + a\_{1}k\_{1} \text{diag}\left(\text{sign}^{a\_{1}-1}\mathbf{n}\_{\varepsilon}\right)\dot{\mathbf{n}}\_{\varepsilon} + ak\_{2} \text{diag}\left(\text{sign}^{a}\dot{\mathbf{n}}\_{\varepsilon}\right)\ddot{\mathbf{n}}\_{\varepsilon}\Big] \\ &= ak\_{2}\mathbf{s}^{T} \text{diag}\left(\text{sign}^{a-1}\dot{\mathbf{n}}\_{\varepsilon}\right)(-k\_{3}\mathbf{s} - k\_{4}\text{sign}\mathbf{s}) \\ &\leq -ak\_{2}\sum\_{i=1}^{3}k\_{3i} \left|\text{sign}^{a-1}\dot{\mathbf{n}}\_{\varepsilon}\right|s\_{i}^{2} - ak\_{2}\sum\_{i=1}^{3}k\_{4i} \left|\text{sign}^{a-1}\dot{\mathbf{n}}\_{\varepsilon i}\right| |s\_{i}| \end{split} \tag{33}$$

The Lyapunov function *<sup>V</sup>*<sup>1</sup> <sup>≥</sup> 0 is selected, and the design <sup>τ</sup> makes . *V*<sup>1</sup> ≤ 0, so it can be concluded that the underactuated USV system is stable.

#### **6. Numerical Simulation and Analysis**

In order to verify the effectiveness of using NESO for full-state feedback with the combination of NFTSM for finite-time trajectory tracking control, the BAYCLASS long-range patrol ship in reference [35] is used for simulation, and the parameters of underactuated USV are shown in Table 1. The virtual USV, actual USV, NESO and controller in the design process are modeled using the Simulink simulation tool in MATLAB, and the solver is chosen ode45 with variable step size. The work of the simulation is divided into two stages. First of all, the simulation verifies whether the mentioned NESO has a good approximation effect on the complex disturbances, positions and velocities state of the underactuated USV system. Finally, the simulation incorporates parameters perturbation using NESO to obtain the three-part state feedback values of positions, velocities, and expansion of the underactuated USV system combined with NFTSM to verify the simulation and compare it with the NTSM in [36] in order to verify the design of the control.


**Table 1.** USV model parameters.

The virtual underactuated USV used to generate the reference trajectory uses the same model parameters as the actual underactuated USV model. The initial values of position and velocity of underactuated USV are η<sup>0</sup> = [0, 0, 0] *<sup>T</sup>* and **v**<sup>0</sup> = [0, 0, *π*/6] *<sup>T</sup>*. The initial values of positions and velocities of virtual underactuated USV are ηd0 = [5, 3.5, *π*/6] *T* and **v**d0 = [0, 0, 0] *<sup>T</sup>*. To reflect the parameter perturbation, **<sup>M</sup>** <sup>=</sup> <sup>±</sup>1.1**M**, **<sup>C</sup>**(**v**) <sup>=</sup> <sup>±</sup>1.1**C**(**v**), **D**(**v**) = ±1.1**D**(**v**) are used in the simulation. Four periods with higher probability of disturbance occurrence were chosen according to the scatter diagram of the sea states. The disturbances τ<sup>E</sup> caused by wind, wave and current are assumed to be:

$$\begin{aligned} \mathbf{r}\_{\rm E} &= \begin{bmatrix} 8 \times 10^{n} \sin(\omega \cdot t) \\ 6 \times 10^{n} \sin(\omega \cdot t) \\ \mathcal{T} \times 10^{n} \sin(\omega \cdot t) \end{bmatrix} \end{aligned} \tag{34}$$

where *n* = 3, 4, *T* = 5, 6, 7, 8*s*, four different periods can be obtained when the angular frequency as *ω* = 1.2566, 1.0472, 0.8976, 0.7854.

The reference trajectory is designed as the complex trajectory obtained by the combination of straight line and circle, and the control input τ<sup>d</sup> of the virtual underactuated USV is selected as follows:

$$\mathbf{r\_d} = \begin{cases} \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}^T \times 10^5, 0 \le t \le 100\\ \begin{bmatrix} 1 & 0 & 2 \end{bmatrix}^T \times 10^5, 100 < t \le 300 \end{cases} \tag{35}$$

The parameters selected for Equation (16) NESO are: λ1*<sup>i</sup>* = 6, λ2*<sup>i</sup>* = 9, λ3*<sup>i</sup>* = 15(*i* = 1, 2, 3), *β<sup>i</sup>* = 0.5, *δ<sup>i</sup>* = 0.001(*i* = 1, 2). Regarding the design form Equation (32) of the controller, the chosen parameters are *k*<sup>1</sup> = 3, *k*<sup>2</sup> = 1, *k*<sup>3</sup> = *k*<sup>4</sup> = diag 10 10 10 , *a* = 9/5, *a*<sup>1</sup> = 5/3. It is worth noting that the symbolic function *k*4sign*s* in Equation (32) will cause obvious chattering in the controller in the actual simulation, so the saturation function is used to replace the symbolic function in the simulation, and the form of the saturation function is shown below.

$$\text{Sat}(s\_i) = \begin{cases} \text{sign}(s\_i) & , |s\_i| > \varpi \\ |s\_i|^\xi \text{sign}(s\_i) / \varpi^\sharp & , |s\_i| \le \varpi \end{cases}, (i = 1, 2, 3) \tag{36}$$

The parameters used are = 7, *ξ* = 0.2.

For the disturbances in the form of Equation (34), there are two different amplitudes. For the first amplitude, the disturbances as in Equation (37) are chosen for simulation.

$$\begin{array}{l} \mathsf{r}\_{\mathrm{E}} = \begin{bmatrix} 8 \times 10^3 \sin(\omega \cdot t) \\ 6 \times 10^3 \sin(\omega \cdot t) \\ \mathcal{I} \times 10^3 \sin(\omega \cdot t) \end{bmatrix} \end{array} \tag{37}$$

Figure 5 shows the tracking curves of the complex disturbances *f*(η, **v**, τE) for four periods under the disturbances of Equation (37). The complex disturbances are related to the USV's own velocity, position information and external environmental disturbances, which are non-periodic disturbances. The red curve, which is the observed value, is very close to the actual value (the blue curve), as can be seen in the figures below.

**Figure 5.** Simulation results of complex disturbances under four periods. (**a**) Complex disturbances tracking curve graph at period *T* = 5*s*; (**b**) Complex disturbances tracking curve graph at period *T* = 6*s*; (**c**) Complex disturbances tracking curve graph at period *T* = 7*s*; (**d**) Complex disturbances tracking curve graph at period *T* = 8*s*.

Figure 6 shows the observation error curves of complex disturbances *f*(η, **v**, τE), positions η<sup>e</sup> and velocities **v***e*. It can be seen in the three figures that for the three different state tracking error curves converge in a minimal neighborhood of zero. The composite disturbance error is less than 10−1, and the position and velocity observation error is less than 10<sup>−</sup>3. Despite the change in the period of the disturbance, the NESO observations are still valid within the error tolerance. There is an initial error in the estimation error of the vessel positions and velocities state at the initial stage. As time goes on, the error converges very quickly. From the convergence range of the error, it can be concluded that the estimated vessel speed and position states can replace the actual measured states information and feedback to the system. Therefore, when the states of the underactuated USV system cannot be measured or is disturbed by the external environment, the structural information of the system itself can be used to estimate the states information of the underactuated USV.

**Figure 6.** Observation error results at four periods. (**a**) Graph of complex disturbances tracking error under four periods; (**b**) Graph of positions tracking error under four periods; (**c**) Graph of velocities tracking error under four periods.

Figures 7–10 show the simulation results under four periods of disturbance *T* = 5s, 6*s*, 7*s*, 8*s*. The trajectory tracking curve, position tracking error curve, velocity tracking error curve and surge force and yaw moment curve are given for each period. For the trajectory tracking curves under four periods, it can be seen that, overall, both control methods can track the reference trajectory, but from the local zoom in, it can be seen that the NFTSM control method used in this study approaches the reference trajectory earlier than the NTSM method and the approximation error is smaller. The small error indicates that the safety is higher and the risk is lower when performing the task using this method. As seen from the velocity tracking error curve, there is not much difference between the two methods because, overall, both methods can make the USV track on the reference trajectory, while achieving the position tracking needs to satisfy the velocity tracking first. For the simulation curves of surge force and yaw moment, it can be seen that the output curve of the controller is smoother than the control method used in this study, and the output curve of the NTSM control method is more strongly chattered. When the periods are *T* = 7*s* and *T* = 8*s*, it is obvious from the position tracking Figures 9b and 10b under the two sets of periods that the control method NFTSM of this paper is better than the NTSM control scheme.

**Figure 7.** Simulation results at disturbance period of *T* = 5*s*. (**a**) Trajectory tracking curves; (**b**) Position tracking error curves; (**c**) Velocity tracking error curves; (**d**) Surge force and yaw moment curves.

**Figure 8.** *Cont*.

**Figure 8.** Simulation results at disturbance period of *T* = 6*s*. (**a**) Trajectory tracking curves; (**b**) Position tracking error curves; (**c**) Velocity tracking error curves; (**d**) Surge force and yaw moment curves.

**Figure 9.** Simulation results at disturbance period of *T* = 7*s*. (**a**) Trajectory tracking curves; (**b**) Position tracking error curves; (**c**) Velocity tracking error curves; (**d**) Surge force and yaw moment curves.

**Figure 10.** Simulation results at disturbance period of *T* = 8*s*. (**a**) Trajectory tracking curves; (**b**) Position tracking error curves; (**c**) Velocity tracking error curves; (**d**) Surge force and yaw moment curves.

The second disturbances amplitude is 10 times the first disturbances amplitude, and the disturbances as in Equation (38) are chosen for simulation.

$$\begin{array}{l} \mathsf{r}\_{\mathrm{E}} = \begin{bmatrix} 8 \times 10^{4} \sin(\omega \cdot t) \\ 6 \times 10^{4} \sin(\omega \cdot t) \\ \mathcal{I} \times 10^{4} \sin(\omega \cdot t) \end{bmatrix} \end{array} \tag{38}$$

Figure 11 presents the observation error curves of complex disturbances, position and velocity under disturbances as shown in Equation (38). For disturbances with amplitude 11 times larger, the observations of the complex disturbances as shown in Figure 11a are not as good as in the case of small disturbances, but the observation errors of the position and velocity states are within 10−2. In the case of period variation, there is no significant difference between the simulation graphs under the four periods, and there is good adaptability for the period variation.

**Figure 11.** Observation error results at four periods. (**a**) Graph of complex disturbances tracking error under four periods; (**b**) Graph of positions tracking error under four periods; (**c**) Graph of velocities tracking error under four periods.

The simulation results shown in Figures 12–15 are obtained by varying the period *T* = 5, 6, 7, 8*s* of the disturbances. Since the disturbances amplitude is 10 times the disturbances in Equation (37), the errors of the observed value obtained by NESO are relatively large compared to the small amplitude disturbances, which have an impact on the position tracking results and lead to a larger error. In the presence of observation errors, both control methods can achieve trajectory tracking control as seen in the trajectory tracking result graph, but the control effects are different. From the simulation results of Figure 12 at the disturbance period of *T* = 5*s*, it is obvious from Figure 12b that the control effect of the blue curve NFTSM of the control method used in this paper has a smaller error, and Figure 12d shows that the chattering of the yaw moment is smaller. When the period becomes larger, as seen from the position error curves of the four periods, the tracking error of the control method in this paper is smaller, and its error curves are closer to the reference error curve of zero value. The simulation graph of the yaw moment at four periods shows that the control scheme proposed in this study has a smaller frequency of control output variation and less chattering.

**Figure 12.** Simulation results at disturbance period of *T* = 5*s*. (**a**) Trajectory tracking curves; (**b**) Position tracking error curves; (**c**) Velocity tracking error curves; (**d**) Surge force and yaw moment curves.

**Figure 13.** *Cont*.

**Figure 13.** Simulation results at disturbance period of *T* = 6*s*. (**a**) Trajectory tracking curves; (**b**) Position tracking error curves; (**c**) Velocity tracking error curves; (**d**) Surge force and yaw moment curves.

**Figure 14.** Simulation results at disturbance period of *T* = 7*s*. (**a**) Trajectory tracking curves; (**b**) Position tracking error curves; (**c**) Velocity tracking error curves; (**d**) Surge force and yaw moment curves.

**Figure 15.** Simulation results at disturbance period of *T* = 8*s*. (**a**) Trajectory tracking curves; (**b**) Position tracking error curves; (**c**) Velocity tracking error curves; (**d**) Surge force and yaw moment curves.

The simulation results of the above two amplitudes and four periods show that the observation errors of position and velocity are less than 10−<sup>3</sup> and the complex disturbance less than 2 × <sup>10</sup>−<sup>2</sup> in the case of small disturbance amplitude. In the case of a large disturbance of 10 times the small amplitude disturbance, the observation errors of position and velocity are less than 1 × <sup>10</sup><sup>−</sup>2, and the observation errors of the complex disturbance are less than 0.3. Under the disturbance of small amplitude and large amplitude disturbances, the position, velocity and composite disturbance observations with a small enough error can be obtained, which shows that the designed full state observer is feasible. The comparison of different periods of the two control methods highlights that the control method NFTSM in this paper has a better trajectory tracking effect, a smaller error and less controller chattering.

#### **7. Conclusions**

In this study, a method combining NESO and NFTSM is presented for trajectory control of an underactuated USV under the disturbances of wind, wave and current in external environment, in the presence of perturbation of the parameters of an underactuated USV model and in the absence of accurate information about the state of the system. Based on theoretical analysis, it is proved that the design process of NESO and NFTSM satisfies the Lyapunov stability condition. Two types of disturbances with 10 times difference in amplitude and four different periods are designed. The four period disturbances with the higher probability of generating wind and wave currents are chosen according to the scatter diagram of the sea states, and the numerical simulation verifies the effectiveness of NESO for all of the states observed. Finally, the NTSM and NFTSM are compared with the position, velocity and complex perturbation observations obtained by NESO, and the following conclusions can be drawn: (1) The stability demonstration and simulation results of NESO show that the designed observer can quickly approximate the actual complex disturbances and the underactuated USV system state, and the vessel position and vessel velocity observation errors for different periods of the two disturbances are less than 10−3, thus it is feasible and easier to implement the observer as a measurement module of the system to feedback the underactuated USV system state. (2) The simulation results show that the actual trajectory has a good approximation for the combined linear and circular trajectory, and the error curve shows the control scheme of NFTSM with a smaller error, and the error finally converges to the small neighborhood of zero, thus the probability of accidents is lower. (3) The simulation results show that under the conditions of environmental disturbances, model parameter uncertainty and state unpredictability, the control method of NFTSM can achieve rapid trajectory tracking control with less tracking error and less chattering in the process.

In future work, the focus will be on adaptive NESO, since the control accuracy is affected by the observer accuracy when using state feedback.

**Author Contributions:** Conceptualization, X.Z., Z.L. and D.X.; Data curation, J.S.; Formal analysis, X.Z., Z.L. and D.X.; Funding acquisition, X.Z. and D.X.; Methodology, D.X.; Resources, X.Z., Z.L. and D.X.; Writing—original draft, D.X.; Writing—review and editing, X.Z., D.X. and J.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Science Fund Project of Heilongjiang Province (LH2021E087), the National Natural Science Foundation of China (51779055), and Fundamental Research Funds for the Central Universities (GK2010260347).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank the fund for supporting, the Editor and Reviewers for their constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Guoqing Zhang \*, Jun Han, Jiqiang Li and Xianku Zhang**

Navigation College, Dalian Maritime University, Dalian 116026, China **\*** Correspondence: zhanggq87@dlmu.edu.cn

**Abstract:** To adapt to complex navigation conditions, this paper addresses the coordination formation of autonomous surface vehicles (ASVs) with the constraint of information interruption. For this purpose, a distributed robust fast finite-time formation control algorithm is proposed by fusion of the directed graph and neural network method. In the strategy, the graph theory is utilized for the channel of information transmission to maintain the stability of the formation system. In addition, the radial basic function (RBF) neural network is employed to approximate the structure uncertainty. Due to the merits of the robust neural damping technique, only two adaptive parameters are designed to compensate the perturbation from the model uncertainty and external environmental. Furthermore, an improved dynamic surface control (DSC) technology is developed for constituting the exponential term of the Lyapunov function. It is proven that the proposed scheme is able to achieve consensus tracking in finite time quickly, and the errors rapidly approach a small region around the origin. Finally, the feasibility and effectiveness of the algorithm are verified by two numerical simulations.

**Keywords:** fast finite-time; formation control; autonomous surface vehicles; adaptive control
