**1. Introduction**

As technology has developed, USVs have been extensively used in various fields. However, due to disturbances from the sea wind, waves, and current, trajectory tracking control is of widespread concern. There have many studies undertaken on the control of the USV trajectory tracking technologies, including the PID controller [1–4], sliding mode control [5–8], backstepping control [9–12], MPC [13–16], adaptive control [17–20] and intelligent control [21–24].

MPC has been developed as an advanced optimization control algorithm based on the superiorities of feedback correction and rolling optimization. MPC has low requirements for the model, and it can solve constrained and multivariable problems. Although the computational load increases with fractional order, presently the development of microprocessors has made it possible for such computation. Predictive compensator-based event-triggered MPC with an NDO strategy has been proposed, and the trajectory tracking control problem of USV subject to input constraints, external disturbances, and cyberattacks has been addressed [14]. However, the output constraints are considered in this paper, and the NDO has been discretized with an approximate discretization method. The adaptive line-of-sight algorithm was developed to obtain an expected heading angle. In addition, MPC was applied to reduce the lateral error, where the sideslip angle compensation was considered [25]. In addition, to obtain accurate state variables in real time, a linear extended-state observer was designed to overcome the influence of environmental disturbances and the nonlinearity of the model. However, the linearization still caused certain deviations in model disturbance estimation. To adjust the controller parameters, the MPC controller has been used to carry out both control allocation and trajectory tracking

**Citation:** Fu, H.; Yao, W.; Cajo, R.; Zhao, S. Trajectory Tracking Predictive Control for Unmanned Surface Vehicles with Improved Nonlinear Disturbance Observer. *J. Mar. Sci. Eng.* **2023**, *11*, 1874. https:// doi.org/10.3390/jmse11101874

Academic Editor: Alessandro Ridolfi

Received: 2 September 2023 Revised: 22 September 2023 Accepted: 23 September 2023 Published: 26 September 2023

**Copyright:** © 2023 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 real time [26]. Furthermore, it has concurrently optimized closed-loop performance with reinforcement learning-based and system-identification methods. To convert chance constraints into deterministic convex constraints, a convex conditional value of risk approximation has been introduced [27]. The converted constraints were further transformed into second-order cone constraints. Then, to account for the external disturbances and fulfill physical constraints, a stochastic model predictive control (SMPC) scheme was used to design the controller. The authors in [28] designed the path planning and control of the USVs simultaneously to overcome the disadvantage of the "first-planning-then-tracking" structure, and the artificial potential field and MPC were combined to solve the planning and tracking problem. The authors in [29] proposed the finite control set model predictive control (FCS-MPC). The more practical control commands formed a limited set of control, namely the thruster propulsion angle and speed of the USV. In addition, a fast and safe collision-avoidance system was designed according to the basics of FCS-MPC, which was applicable to varying environments. The authors in [30], to guarantee precise and stable trajectory tracking performance for AUVs, proposed a novel control architecture based on model-free control principles. The combination of intelligent PID and PD feedforward control had good performance for trajectory tracking accuracy, disturbance rejections, and initial tracking error compensations. The authors in [31], to solve the surge-motion tracking control problem of an autonomous undersea vehicle (AUV) with system constraints, proposed an adaptive backstepping control scheme. Both a state feedback control scheme and an output feedback control scheme were developed for AUVs with deferred output constraints. The authors in [32], to pay more attention to the characteristics of flexibility and accessibility, proposed a fusion framework of field theoretical planning and an MPC algorithm. In addition, the trajectory smoothness and collision-avoidance constraints under a complex environment were considered. The authors in [33], to solve the fault-tolerant trajectory tracking control problem of twin-propeller non-rudder USVs subject to propeller faults, proposed an adaptive fault-tolerant trajectory tracking control scheme by utilizing the excellent nonlinearity approximation performance of neural networks. The authors in [34], to achieve autonomous cooperative formation control of underactuated USVs in a complex ocean environment, adopted a dual MPC approach based on a virtual trajectory. The authors in [35] showed that time-varying external disturbances affected the accuracy of trajectory tracking. To ensure trajectory tracking accuracy, a reduced-order extended-state observer and the super-twisting second-order sliding mode controller were proposed. The authors in [36], to solve the lumped uncertainties caused by input quantization, actuator faults, and dead zones, proposed an adaptive sliding mode tracking controller for USVs with predefined time performance. In the control design, adaptive control gains were established based on barrier functions. In addition, a predefined time-adaptive SMC scheme was adopted by introducing an auxiliary function. The authors in [37], to set the velocity of the USV converge to zero at the berth, adopted an interpolation approach to densify the waypoints at the end of the berthing trajectory. In addition, to improve computational performance during the USV berthing, an event-triggered adaptive horizon MPC approach was adopted.

Controllers are usually poorly tuned to USV motion systems, and the disturbancerejection performance of controllers is not satisfactory due to disturbances from ocean waves, wind, and currents. Therefore, an improved NDO-based MPC method for trajectory tracking control of USVs is proposed in this paper. First, the MPC is designed for USV trajectory tracking. Then, an NDO is designed to estimate the disturbance of ocean wind, waves, and currents, which has fewer parameters. In addition, Lyapunov stability is analyzed for the overall system. Finally, the proposed method is verified by simulation experiments. The main contributions of this work can be summarized as follows: the output constraints are considered in this paper, and the NDO has been discretized with an approximate discretization method; the disturbances are compensated for by an improved NDO, and better trajectory tracking performances are obtained for USVs; also, with the improved NDO, calculation time is saved compared with the traditional NDO.

The rest of the paper is structured as follows: Section 2 introduces the USV kinematics model and dynamic model. In Section 3, the improved NDO-based MPC is designed for the USV, and Lyapunov stability is analyzed for the overall system. Comparison simulations and discussion of the results is performed in Section 4. Section 5 presents the conclusions.

#### **2. State-Space Model of Unmanned Surface Vehicles**

In general, the kinematics modeling of USVs is represented by the North–East–Down coordinate system, while the dynamics model is built in the ship coordinate system. The North–East–Down coordinate system is also called the geodetic coordinate system or the inertial coordinate system. It is usually used as the reference system, and any point on the sea can be used as the origin of the system. The ship coordinate system changes with the motion of the ship and can describe the force, moment, linear velocity, and angular velocity of a USV in various degrees of freedom.

To simplify the model, a USV model is utilized with three degrees of freedom for trajectory tracking control. The motions of yaw, surge, and sway are the most important for the trajectory tracking of the USVs, so the roll, pitch, and heave of the USVs are ignored. Thus, the USV model can be represented in the ship coordinate system and the geodetic coordinate system. This is shown in Figure 1.

**Figure 1.** Three degrees of freedom of motion model for USVs.

In Figure 1, it can be seen the North–East–Down coordinate system is represented by *NOEE*, and the coordinate system for the ship is described by *XOY*. The kinematics and dynamics model of the ship can be represented as:

$$
\dot{\boldsymbol{\eta}} = \mathbf{R}(\boldsymbol{\psi})\boldsymbol{\upsilon} \tag{1}
$$

$$
\mathbf{M}\dot{\boldsymbol{\nu}} + \mathbf{C}(\boldsymbol{\nu})\boldsymbol{\nu} + \mathbf{D}(\boldsymbol{\nu})\boldsymbol{\nu} = \boldsymbol{\pi} + \boldsymbol{\pi}\_{\mathbf{d}} \tag{2}
$$

where η = [*x*, *y*, *ψ*] *<sup>T</sup>* represents the *x*, *y* position and heading angle vector of the USV in the inertial coordinate system; *x* and *y* represent the position of the ship with regard to the North–East–Down coordinate system, *ψ* represents the yaw angle information of the USV; υ = [*u*, *v*,*r*] *<sup>T</sup>* is the vector of the velocity and angular velocity for the USV in the ship coordinate system; *u*, *v*, and *r* are the surge velocity, sway velocity, and yaw angular velocity of USV, respectively; *τu*, *τv*, and *τ<sup>r</sup>* represent the surge thrust, sway thrust, and yaw moment, which are the control inputs of the system; τ**<sup>d</sup>** = [*τud*, *τvd*, *τrd*] *<sup>T</sup>* is the corresponding thrust and moment caused by the time-varying external disturbances. **R**(*ψ*) is the rotation matrix with the relationship of **R**−1(*ψ*) = **R***T*(*ψ*); **M** represents the inertial

matrix of USV, where **M** = **M***<sup>T</sup>* > 0; **C**(υ) represents the Coriolis centripetal force matrix, and **C**(υ) = −**C**(υ) *<sup>T</sup>*; and **D**(υ) is the nonlinear hydrodynamic damping matrix. The detailed information for the matrixes is shown as follows:

$$\mathbf{R}(\boldsymbol{\psi}) = \begin{bmatrix} \cos \psi & -\sin \psi & 0 \\ \sin \psi & \cos \psi & 0 \\ 0 & 0 & 1 \end{bmatrix}, \mathbf{M} = \begin{bmatrix} m\_{11} & 0 & 0 \\ 0 & m\_{22} & m\_{23} \\ 0 & m\_{32} & m\_{33} \end{bmatrix} \tag{3}$$

$$\mathbf{C}(\boldsymbol{\upsilon}) = \begin{bmatrix} 0 & 0 & -m\_{22}\boldsymbol{\upsilon} \\ 0 & 0 & m\_{11}\boldsymbol{\mu} \\ m\_{22}\boldsymbol{\upsilon} & -m\_{11}\boldsymbol{\mu} & 0 \end{bmatrix}, \ \mathbf{D}(\boldsymbol{\upsilon}) = -\begin{bmatrix} d\_{11} & 0 & 0 \\ 0 & d\_{22} & d\_{23} \\ 0 & d\_{32} & d\_{33} \end{bmatrix} \tag{4}$$

According to Equations (3) and (4), the reduced kinematics and dynamics equations can be represented as:

$$\begin{cases} \dot{x} = u\cos\psi - v\sin\psi\\ \dot{y} = u\sin\psi + v\cos\psi\\ \dot{\psi} = r \end{cases} \tag{5}$$

$$\begin{cases} m\_{11}\dot{u} - m\_{22}vr + d\_{11}u = \mathfrak{r}\_{\mathfrak{u}} + \mathfrak{r}\_{\mathfrak{u}d} \\ m\_{22}\dot{v} + m\_{23}\dot{r} + m\_{11}ur + d\_{22}v + d\_{23}r = \mathfrak{r}\_{\mathfrak{v}} + \mathfrak{r}\_{\mathfrak{v}d} \\ m\_{32}\dot{v} + m\_{33}\dot{r} + (m\_{22} - m\_{11})uv + d\_{32}v + d\_{33}r = \mathfrak{r}\_{\mathfrak{r}} + \mathfrak{r}\_{\mathfrak{r}d} \end{cases} \tag{6}$$

For the USVs, there are constraints for the actuators and the outputs. In addition, they are described as follows.

The increment constraints for inputs can be represented as:

$$\begin{array}{c} \Delta \mathbf{u}\_{\min}(t+k) \le \Delta \mathbf{u}(t+k) \le \Delta \mathbf{u}\_{\max}(t+k) \\\ k = 0, 1, \cdots, N\_c - 1 \end{array} \tag{7}$$

where *Nc* denotes the control horizon.

The upper and lower-limit constraints for inputs can be represented as:

$$\begin{array}{c} \mathbf{u}\_{\min}(t+k) \le \mathbf{u}(t+k) \le \mathbf{u}\_{\max}(t+k) \\\ k = 0, 1, \cdots, N\_c - 1 \end{array} \tag{8}$$

The outputs of speed increment constraints can be represented as:

$$\begin{array}{c} \Delta \mathfrak{v}\_{\min}(t+k) \le \Delta \mathfrak{v}(t+k) \le \Delta \mathfrak{v}\_{\max}(t+k) \\\ k = 0, 1, \dots, N\_c - 1 \end{array} \tag{9}$$

The outputs of speed upper and lower-limit constraints can be represented as:

$$\begin{array}{c} \mathfrak{w}\_{\min}(t+k) \le \mathfrak{w}(t+k) \le \mathfrak{w}\_{\max}(t+k) \\\ k = 0, 1, \cdots, N\_c - 1 \end{array} \tag{10}$$

The terminal equality constraint can be represented as:

$$\left\|\mathfrak{n}(k+N|t) - \mathfrak{n}\_{\mathfrak{r}}(k+N|t)\right\|^2 = 0\tag{11}$$

where η denotes the vector to be controlled, and η**<sup>r</sup>** denotes the reference trajectory.

#### **3. Nonlinear Disturbance Observer-Based Model Predictive Control**

In this section, NDO-based MPC is designed for the three-degrees-of-freedom kinematics and dynamics of a USV with the state-space model. The design schematic diagram of NDO-based MPC is shown in Figure 2.

**Figure 2.** Schematic diagram of MPC based on a nonlinear observer.

According to Equations (1) and (2), the state-space equation for the ship can be rewritten as:

$$\begin{cases} \dot{\boldsymbol{\eta}} = \mathbf{R}(\boldsymbol{\varrho}) \boldsymbol{\upsilon} \\ \dot{\boldsymbol{\upsilon}} = \mathbf{M}^{-1} (\boldsymbol{\tau} + \boldsymbol{\tau\_d} - \mathbf{C}(\boldsymbol{\upsilon}) \boldsymbol{\upsilon} - \mathbf{D}(\boldsymbol{\upsilon}) \boldsymbol{\upsilon}) \end{cases} \tag{12}$$

From Equation (12), the matrixes **R**, **C**, and **D** are nonlinear variables, which gives the model strong nonlinearity.

#### *3.1. Model Predictive Control Design of an Unmanned Surface Vehicle*

3.1.1. Discrete Linearization of an Unmanned Surface Vehicle Model

USV is a complex system with large inertia and time-delay characteristics. During navigation, it is subjected to various forces such as thrust, hydrodynamic force, hydrostatic force, and external disturbances. Therefore, the USV model has obvious nonlinear characteristics. In this paper, a simplified model of the USV is applied, and the uncertainty is dealt with by the NDO together with the disturbances from the ocean environment. In addition, the model is linearized.

In this paper, the NDO is designed to compensate for the disturbance. Therefore, a linear model of the USV is sufficient for MPC design, where the uncertainty from the linearization of the USV model can be solved with an NDO. Therefore, the model is linearized first and then discretized. Finally, the optimal control sequence is obtained according to the linear model predictive control.

The linearization of a nonlinear system can be divided into approximate linearization and exact linearization. Among them, the approximate linearization method is relatively simple with high applicability but low accuracy. The precision of the accurate linearization method is high, but it should be noted that a special case analysis of a single system is difficult and has poor universality.

The reference trajectory is represented in Equation (13) with environmental disturbances. The first-order Taylor expansion can be obtained at any point (**xr**, **ur**); then approximate linearization is achieved by leaving the higher-order terms. This can be seen in Equation (14). **.**

$$
\dot{\mathbf{x}}\_{\mathbf{t}} = f(\mathbf{x}\_{\mathbf{t}}, \mathbf{u}\_{\mathbf{t}}) \tag{13}
$$

$$\begin{aligned} \dot{\mathbf{x}} &= f(\mathbf{x}\mathbf{r}, \mathbf{u}\mathbf{r}) + \left. \frac{\partial f}{\partial \mathbf{x}} \right| \mathbf{x} = \mathbf{x}\mathbf{r} \left( \mathbf{x} - \mathbf{x}\mathbf{r} \right) + \left. \frac{\partial f}{\partial \mathbf{u}} \right| \mathbf{x} = \mathbf{x}\mathbf{r} \left( \mathbf{u} - \mathbf{u}\mathbf{r} \right) \\ \mathbf{u} &= \mathbf{u}\mathbf{r} \end{aligned} \tag{14}$$

Subtract Equation (13) from (14), and the following equation can be achieved.

$$\dot{\tilde{\mathbf{x}}} = \mathbf{A}\tilde{\mathbf{x}} + \mathbf{B}\tilde{\mathbf{u}} \tag{15}$$
  $\text{and } \tilde{\mathbf{x}} = \mathbf{x} - \mathbf{x}\_t, \tilde{\mathbf{u}} = \mathbf{u} - \mathbf{u}\_t, \mathbf{A} = \frac{\partial f}{\partial \mathbf{x}} \Big|\_{\mathbf{x}} = \mathbf{x} \mathbf{r}' \quad \mathbf{B} = \frac{\partial f}{\partial \mathbf{u}} \Big|\_{\mathbf{x}} \mathbf{x} = \mathbf{x} \mathbf{r}' $  
$$\mathbf{u} = \mathbf{u} \mathbf{r} \qquad \mathbf{u} = \mathbf{u} \mathbf{r}$$

The discrete form for Equation (15) is shown as follows:

$$\begin{cases} \tilde{\mathbf{x}}(k+1) = \mathbf{A\_d}\tilde{\mathbf{x}}(k) + \mathbf{B\_d}\tilde{\mathbf{u}}(k) \\ \tilde{\mathbf{y}}(k) = \mathbf{C}\tilde{\mathbf{x}}(k) \end{cases} \tag{16}$$

where **Ad** <sup>=</sup> *<sup>T</sup>***<sup>A</sup>** <sup>+</sup> **<sup>I</sup>**, **Bd** <sup>=</sup> *<sup>T</sup>***B**, **<sup>y</sup>**1(*k*) = **<sup>y</sup>**(*k*) <sup>−</sup> **yr**(*k*), *<sup>T</sup>* is the discretization step.

#### 3.1.2. Objective Function Design

⎣

**R**

⎦ *Nc*×*Nc*

To ensure that the USV can track the trajectory smoothly and quickly, the cost function shown in Equation (17) is applied, which concerns the increments of the control variables and the errors of the system states.

$$\min \quad J = \min \left\{ \sum\_{i=1}^{N\_p} \left\| \eta(k+i) - \eta\_\mathbf{r}(k+i) \right\|\_\mathbf{Q}^2 + \sum\_{i=1}^{N\_c-1} \left\| \Delta \mathbf{u}(k+i) \right\|\_\mathbf{R}^2 \right\} \tag{17}$$

where **Q** and **R** are the weight matrixes for tracking errors and increments of the control variables, respectively; and *Np* is the prediction horizon. Δ**u**(*k* + *i*) is the variable of the increments of the control variables, so it can be obtained as follows:

$$\mathfrak{E}(k) = \begin{bmatrix} \tilde{\mathbf{x}}(k) \\ \tilde{\mathbf{u}}(k-1) \end{bmatrix} \tag{18}$$

The new state-space equation is represented as:

$$\begin{cases} \pounds(k+1) = \dot{\mathbf{A}}\pounds(k) + \dot{\mathbf{B}}\Delta\mathbf{u}(k) \\ \vounds(k) = \overline{\mathbf{C}}\pounds(k) \end{cases} \tag{19}$$

where **<sup>A</sup>**<sup>1</sup> <sup>=</sup> **Ad Bd 0 INu** , **<sup>B</sup>**<sup>1</sup> <sup>=</sup> **Bd INu** , **<sup>C</sup>**<sup>1</sup> = [**INx <sup>0</sup>**], *Nu* denotes the number of control variables, *Nx* denotes the number of state variables.

Therefore, the system prediction outputs can be calculated as follows:

$$\begin{array}{l} \mathbf{Y} = \mathbf{Y}\boldsymbol{\xi}(k) + \mathbf{H}\boldsymbol{\tilde{u}}(k) \\ J = \frac{1}{2}\Delta\mathbf{u}(k)^{T}\mathbf{H}\_{I}\Delta\mathbf{u}(k) + \mathbf{f}\_{I}\Delta\mathbf{u}(k) \end{array} \tag{20}$$

$$\begin{array}{rclcrcl}\hline\\\mathbf{0}&\mathbf{H}\_{\parallel}&=&\mathbf{2}\begin{pmatrix}\mathbf{H}^{T}\mathbf{Q}\_{\mathbf{c}}\mathbf{H}+\mathbf{R}\_{\mathbf{c}}\end{pmatrix},&\mathbf{f}\_{\parallel}&=&\mathbf{2}\mathbf{V}\mathbf{\hat{x}}(k)\mathbf{Q}\_{\mathbf{c}}\mathbf{H},&\mathbf{Q}\_{\mathbf{c}}&=&\begin{bmatrix}\mathbf{Q} & &\\&\ddots& &\\&\mathbf{Q}\end{bmatrix}\_{N\_{\mathbf{p}}\times N\_{\mathbf{p}}},\\\\ \mathbf{R}\_{\mathbf{c}}&=&\begin{bmatrix}\mathbf{R} & &\\&\ddots& \\\\&&\ddots& \end{bmatrix}\_{N\_{\mathbf{p}}\times N\_{\mathbf{p}}},\\\\ \end{array}$$

#### *3.2. Nonlinear Disturbance Observer Design of Unmanned Surface Vehicles*

To make MPC for a USV more applicable, a nonlinear disturbance observer is designed, which can estimate and compensate for the external environmental disturbance received by the USV. Therefore, the stability and anti-disturbance performance is improved for the USV, while capsizing and unnecessary navigation accidents are avoided for the USV.

According to the mathematical model of the USV, the state equation is designed as follows:

$$\dot{\mathbf{\dot{r}\_d}} = \mathbf{K\_0}(\tau\_\mathbf{d} - \dot{\tau}\_\mathbf{d}) = -\mathbf{K\_0}\dot{\tau}\_\mathbf{d} + \mathbf{K\_0}(\mathbf{M}\dot{\mathbf{v}} + \mathbf{C}(\mathbf{v})\mathbf{v} + \mathbf{D}(\mathbf{v})\mathbf{v} - \mathbf{r}) \tag{21}$$

where **K0** is a three-dimensional positive definite matrix; the estimated disturbance values τˆ **<sup>d</sup>** can be specifically written as τˆ **<sup>d</sup>** = [*τ*ˆ*ud*, *τ*ˆ*vd*, *τ*ˆ*rd*] *<sup>T</sup>*, which are the estimated values of surge disturbance, sway disturbance, and yaw direction disturbance.

It can be seen from Equation (12) that υ of a USV can be obtained directly, but the derivative term of the speed state variable **.** υ cannot be obtained directly. Therefore, it is necessary to improve the disturbance observer. The variable β can be selected as the intermediate assignment variable of the observer, which is expressed as follows:

$$
\beta = \dot{\mathbf{\tau}}\_{\mathbf{d}} - \mathbf{K}\_{\mathbf{0}} \mathbf{M} \mathbf{v} \tag{22}
$$

Then, **.**

$$\begin{array}{lcl}\dot{\mathfrak{F}} &= \dot{\mathfrak{r}}\_{\mathsf{d}} - \mathsf{K}\_{\mathsf{0}} \mathsf{M} \dot{\mathfrak{v}} \\ &= \mathsf{K}\_{\mathsf{0}} (\mathsf{M} \dot{\mathfrak{v}} + \mathsf{C}(\mathsf{v}) \mathfrak{v} + \mathsf{D}(\mathsf{v}) \mathfrak{v} - \mathsf{r}) - \mathsf{K}\_{\mathsf{0}} \dot{\mathfrak{r}}\_{\mathsf{d}} - \mathsf{K}\_{\mathsf{0}} \mathsf{M} \dot{\mathfrak{v}} \\ &= -\mathsf{K}\_{\mathsf{0}} (\mathsf{f} + \mathsf{K}\_{\mathsf{0}} \mathsf{M} \mathsf{v}) + \mathsf{K}\_{\mathsf{0}} (\mathsf{C}(\mathsf{v}) \mathfrak{v} + \mathsf{D}(\mathsf{v}) \mathfrak{v} - \mathsf{r}) \\ &= -\mathsf{K}\_{\mathsf{0}} \mathsf{f} - \mathsf{K}\_{\mathsf{0}} (\mathsf{K}\_{\mathsf{0}} \mathsf{M} \mathsf{v} - \mathsf{C}(\mathsf{v}) \mathsf{v} - \mathsf{D}(\mathsf{v}) \mathsf{v} + \mathsf{r}) \end{array} \tag{23}$$

Therefore, the new equation of the improved NDO is:

$$\begin{aligned} \dot{\mathbf{\dot{r}}}\_{\mathbf{d}} &= \boldsymbol{\beta} + \mathbf{K}\_{\mathbf{0}} \mathbf{M} \boldsymbol{\upsilon} \\ \dot{\boldsymbol{\beta}} &= -\mathbf{K}\_{\mathbf{0}} \boldsymbol{\beta} - \mathbf{K}\_{\mathbf{0}} (\mathbf{K}\_{\mathbf{0}} \mathbf{M} \boldsymbol{\upsilon} - \mathbf{C}(\boldsymbol{\upsilon}) \boldsymbol{\upsilon} - \mathbf{D}(\boldsymbol{\upsilon}) \boldsymbol{\upsilon} + \boldsymbol{\pi}) \end{aligned} \tag{24}$$

The improved equation avoids the calculation of **.** υ and simplifies the calculation process. Therefore, it can improve calculation efficiency and save calculation time.

Equation (24) is in a continuous form, and cannot be directly used for MPC design. Therefore, it is discretized with the approximate discretization method:

$$\begin{array}{c} \mathsf{T\_{d}}(k+1) = \mathsf{\mathcal{B}\_{d}} + \mathsf{K\_{d1}}\mathsf{M}\mathsf{v} + \mathsf{\tau\_{d}}(k) \\ \mathsf{\mathcal{B}}(k+1) = -\mathsf{K\_{d2}}\mathsf{\mathcal{B}}(k) - \mathsf{K\_{d1}}(\mathsf{K\_{0}}\mathsf{M}\mathsf{v} - \mathsf{C}(\mathsf{v})\mathsf{v} - \mathsf{D}(\mathsf{v})\mathsf{v} + \mathsf{\tau}) \end{array} \tag{25}$$

and β**<sup>d</sup>** = *T*β, **Kd1** = *T***K0**, and **Kd2** = *T***K0** − **I**.

#### *3.3. Stability Analysis of Unmanned Surface Vehicles*

3.3.1. Stability Analysis of the Model Predictive Control of Unmanned Surface Vehicles

To verify the stability of the USV control system under MPC, the Lyapunov function defined as *V*0(*k*) is selected:

$$V^0(k) = \min\_{\Delta v} \sum\_{i=1}^{N\_p} \left\| \mathbf{n}(k+i|t|) - \mathbf{n}\_\mathbf{d}(k+i|t|) \right\|\_\mathbf{Q} + \sum\_{i=1}^{N\_\varepsilon - 1} \left\| \Delta \mathbf{u}(k+i|t|) \right\|\_\mathbf{R} \tag{26}$$

If the control horizon *N*c is defined to be equal to the prediction horizon *Np*, the above equation can be simplified as follows:

$$V^0(k) = \min\_{\Delta v} \sum\_{i=1}^N \left\| \eta(k + i|t) - \eta\_\mathbf{d}(k + i|t) \right\|\_\mathbf{Q}^2 + \left\| \Delta \mathbf{u}(k + i|t) \right\|\_\mathbf{R}^2 \tag{27}$$

The quadratic function is always greater than 0, so its positive definiteness is proven:

$$V^0(k) \ge 0\tag{28}$$

Then, we only need to prove that *V*0(*k*) is decreasing, then its stability is proven.

$$\begin{aligned} V^{0}(k+1) &= \underset{\Delta\mathbf{r}}{\text{min}} \left\{ \sum\_{i=1}^{N} \left\| \begin{array}{l} \left( \left\| \begin{array}{l} \left( \left\| \begin{array}{l} \left( \left\| \begin{array}{l} \left( \left| \left< i \left( \left| \left.left \left| \left.left \left. \right| \right. \right. \right. \right. \right. \right. \right. \right. \right. \end{array} \right. \end{array} \right. \right) \right\|\_{\mathbf{R}}^{2} \right\} \\ &= \underset{\Delta\mathbf{r}}{\text{min}} \left\{ \left\| \begin{array}{l} \left( \left\| \begin{array}{l} \left( \left| \left. \left$$

The terminal equation constraint is:

$$\|\|\eta(k+1+N|t) - \eta\_{\mathbf{r}}(k+1+N|t)\|\|\_{\mathbf{Q}}^2 = 0 \tag{30}$$

Furthermore,

$$\min\_{\mathbf{A}\mathbf{u}} \left\{ \left\| \boldsymbol{\eta} \left( k + 1 + N|t \right) - \boldsymbol{\eta}\_{\mathbf{d}} \left( k + 1 + N|t \right) \right\|\_{\mathbf{Q}}^2 + \left\| \Delta \mathbf{u} \left( k + 1 + N|t \right) \right\|\_{\mathbf{R}}^2 \right\} = 0 \tag{31}$$

$$\|\|\mathbf{\eta}(k+i|t) - \mathbf{\eta}\_{\mathbf{r}}(k+i|t)\|\|\_{\mathbf{Q}}^2 + \|\Delta \mathbf{u}(k+i|t)\|\_{\mathbf{R}}^2 = 0 \tag{32}$$

Therefore, *<sup>V</sup>*0(*<sup>k</sup>* + <sup>1</sup>) ≤ *<sup>V</sup>*0(*k*), and the stability of MPC is proven.

3.3.2. Stability Analysis of the Nonlinear Disturbance Observer of Unmanned Surface Vehicles

To verify the stability of the NDO and ensure that it can be applied to the trajectory tracking control system, it is necessary first to define the variable of the difference between the observer's estimated value and the actual value of the external disturbance to the USV:

$$
\widetilde{\mathbf{r}}\_{\mathbf{d}} = \mathbf{r}\_{\mathbf{d}} - \mathbf{\hat{r}}\_{\mathbf{d}} \tag{33}
$$

Considering the kinematic Equations (2), (24) and (33) of the USV, then by calculating the derivative of time on both sides of Equation (27), the formula can be represented as follows:

$$\begin{aligned} \dot{\mathsf{T}\_{\mathsf{d}}} &= \dot{\mathsf{B}} + \mathsf{K\_{0}} \mathsf{M} \dot{\mathsf{v}} \\ &= -\mathsf{K\_{0}} \mathsf{f} - \mathsf{K\_{0}} (\mathsf{K\_{0}} \mathsf{M} \mathsf{v} - \mathsf{C}(\mathsf{v}) \mathsf{v} - \mathsf{D}(\mathsf{v}) \mathsf{v} + \mathsf{r}) + \mathsf{K\_{0}} (-\mathsf{C}(\mathsf{v}) \mathsf{v} - \mathsf{D}(\mathsf{v}) \mathsf{v} + \mathsf{r} + \mathsf{r\_{d}}) \\ &= \mathsf{K\_{0}} (\mathsf{r\_{d}} - (\mathsf{g} + \mathsf{K\_{0}} \mathsf{M} \mathsf{v})) \\ &= \mathsf{K\_{0}} \tilde{\mathsf{r\_{d}}} \end{aligned} \tag{34}$$

From Equation (34), one can calculate the derivative of both sides of Equation (33) with respect to time, and simplify it to obtain:

$$
\dot{\tilde{\mathbf{r}}}\_{\mathbf{d}} = \dot{\mathbf{r}}\_{\mathbf{d}} - \dot{\tilde{\mathbf{r}}}\_{\mathbf{d}} = \dot{\mathbf{r}}\_{\mathbf{d}} - \mathbf{K}\_{\mathbf{0}} \tilde{\mathbf{r}}\_{\mathbf{d}} \tag{35}
$$

The Lyapunov method is used to verify the stability of the disturbance observer, and the appropriate Lyapunov function is selected as follows:

$$V\_d = \frac{1}{2} \widetilde{\mathbf{r}}\_\mathbf{d}^\mathrm{T} \widetilde{\mathbf{r}}\_\mathbf{d} \tag{36}$$

According to Equation (35), the derivative of both sides of Equation (36) with respect to time can be obtained as follows:

$$
\dot{V}\_d = \overline{\mathbf{r}}\_\mathbf{d}^\mathrm{T} \dot{\overline{\mathbf{r}}}\_\mathbf{d} = -\overline{\mathbf{r}}^\mathrm{T} \mathbf{K}\_0 \overline{\mathbf{r}} + \overline{\mathbf{r}}\_\mathbf{d}^\mathrm{T} \dot{\mathbf{r}}\_\mathbf{d} \tag{37}
$$

According to Young's inequality theory,

$$
\widetilde{\mathbf{r}}\_{\mathbf{d}}^{\mathrm{T}} \dot{\mathbf{r}}\_{\mathbf{d}} \le a\_1 \widetilde{\mathbf{r}}\_{\mathbf{d}}^{\mathrm{T}} \widetilde{\mathbf{r}}\_{\mathbf{d}} + \frac{\mathbf{C}\_d^2}{4a\_1} \tag{38}
$$

and *a*<sup>1</sup> > 0, *Cd* is the limit of the disturbance change rate.

From Equations (37) and (38), this can be written as inequality:

$$\dot{V}\_d \le -\tilde{\mathbf{r}}^T \mathbf{K}\_\mathbf{0} \tilde{\mathbf{r}} + a\_1 \tilde{\mathbf{r}}\_\mathbf{d}^T \tilde{\mathbf{r}}\_\mathbf{d} + \frac{\mathbf{C}\_d^2}{4a\_1} \le -2(\lambda\_{\text{min}}(\mathbf{K}\_\mathbf{0}) - a\_1)V\_d + \frac{\mathbf{C}\_d^2}{4a\_1} \tag{39}$$

Take

$$\begin{cases} \; \mu\_0 = 2(\lambda\_{\text{min}}(\mathbf{K\_0}) - a\_1) > 0\\ \; \mathcal{C}\_0 = \frac{\mathcal{C}\_d^2}{4a\_1} > 0 \end{cases} \tag{40}$$

and *λ*min(**K0**) > *a*1. Equation (39) can be abbreviated as follows:

$$
\dot{V}\_d \le -\mu\_0 V\_d + \mathcal{C}\_0 \tag{41}
$$

According to Equation (36), *Vd* is always greater than 0, and from (41), the result can be obtained as follows:

$$0 \le V\_d \le \frac{C\_0}{\mu\_0} + (V\_d(0) - \frac{C\_0}{\mu\_0})e^{-\mu\_0 t} \tag{42}$$

According to Equation (42), it shows that the Lyapunov function *Vd* stays in a closed ball of some radius whose origin is the center of the sphere. In addition, it is uniformly ultimately bounded. In addition, the radius of the sphere is *RVd* <sup>=</sup> *<sup>C</sup>*<sup>0</sup> *<sup>μ</sup>*<sup>0</sup> <sup>=</sup> *<sup>C</sup>*<sup>2</sup> *d* <sup>8</sup>*a*1(*λ*min(**K0**)−*a*1). According to Equation (36), it can be found that the disturbance estimation error variable <sup>τ</sup>1**<sup>d</sup>** also converges to the sphere radius *<sup>R</sup>*τ1**<sup>d</sup>** <sup>=</sup> 2*C*<sup>0</sup> *<sup>μ</sup>*<sup>0</sup> = √ *Cd* <sup>4</sup>*a*1(*λ*min(**K0**)−*a*1) with the origin as the center of the sphere. At the same time, it can also be known that if the external environmental disturbance value τ**<sup>d</sup>** of USV is an arbitrary unknown constant value, then the boundary of the disturbance charge rate is *Cd* = 0. According to Equation (39), the observer estimation error value <sup>τ</sup>1**<sup>d</sup>** can converge to the origin.

According to Equation (40), as long as the appropriate observer parameters *a*<sup>1</sup> and **K0** can be selected, an arbitrarily small error convergence radius *<sup>R</sup>τ*1*<sup>d</sup>* can be obtained. In other words, NDO can estimate the external environmental disturbance suffered by the USV according to an arbitrarily small error, and the estimation accuracy depends on the selected parameters.

#### **4. Results and Discussions**

To verify the influence of the improved NDO-based MPC, it is applied for trajectory tracking control of a USV called CyberShip II. The tracking errors and performance of the USV are shown in this section. In addition, the computational efficiency of the improved NDO is verified.

#### *4.1. Model Parameters of Unmanned Surface Vehicle*


In the simulation, η = [000◦] T, which is set as the initial state of USV; υ = [000] T, which is set as the initial speed of USV. In addition, the reference trajectories of USV are shown as follows:

$$\begin{cases} \mathbf{x}\_d = 2\sin(0.02t), & 0 \le t \le 500 \\ y\_d = 2 - 2\cos(0.02t), & 0 \le t \le 500 \\ \psi\_d = \arctan(\dot{\mathbf{x}}\_d/\dot{y}\_d), & 0 \le t \le 500 \end{cases} \tag{43}$$

$$\begin{cases} \mathbf{x}\_d = 8 \sin(0.02t), & 0 \le t \le 500 \\ y\_d = t, & 0 \le t \le 500 \\ \boldsymbol{\upvarphi}\_d = \arctan(\dot{\boldsymbol{x}}\_d / \dot{\boldsymbol{y}}\_d), & 0 \le t \le 500 \end{cases} \tag{44}$$

*T* = 0.1, which is set as the simulation sampling time; *Nx* = 6, which is set as the number of states; *Nu* = 3, which is set as the number of control variables; *Np* = 20, which is set as the prediction horizon; *Nc* = 10, which is set as the control horizon.

#### *4.2. Simulation Results and Analysis*

In MATLAB 2021a, NDO-based MPC for trajectory tracking control of a USV is compared with MPC without an observer. Simulations with different disturbances are performed to verify the anti-disturbance and robustness performances. The performances of the improved NDO are discussed.

The composite model of external disturbances meets the requirements of Level 3 sea conditions. The specific external wind, wave, and current disturbances settings are as follows:

$$\begin{cases} \tau\_{du} = m\_{11}h(\mathbf{s})w\_{ll}(\mathbf{s})\\ \tau\_{dv} = m\_{22}h(\mathbf{s})w\_{v}(\mathbf{s})\\ \tau\_{dr} = m\_{33}h(\mathbf{s})w\_{r}(\mathbf{s}) \end{cases} \tag{45}$$

where wave transfer function *h*(*s*) = *<sup>K</sup>ω<sup>s</sup> s*2+2*ζω*0*s*+*ω*<sup>2</sup> 0 , *K<sup>ω</sup>* = 2*ζω*0*σ*0, *ω*0, *ζ*, and *σ*<sup>0</sup> represent wave frequency, wave strength gain, and damping constant, respectively; *K<sup>ω</sup>* = 0.255, *ω*<sup>0</sup> = 0.808. *wu*(*s*), *wv*(*s*), and *wr*(*s*) represent random white-noise disturbances, then the noise power is set to 0.01, 0.005, and 0.1, respectively.

The simulation results of improved NDO-based MPC for trajectory tracking control and traditional MPC with disturbances, when the reference trajectory is a circle, are shown in Figure 3. The comparisons of improved NDO-based MPC for trajectory tracking control and traditional MPC with disturbances, when the reference trajectory is sinusoidal, are shown in Figure 4. The calculation times between the two NDO methods with disturbances are shown in Figure 5. The trajectory tracking errors of the three trajectory tracking methods with different parameters are shown in Figure 6.

**Figure 3.** Comparisons of improved NDO-based MPC for trajectory tracking control and traditional MPC with disturbances when the reference trajectory is a circle.

**Figure 4.** *Cont*.

**Figure 4.** Comparisons of improved NDO-based MPC for trajectory tracking control and traditional MPC with disturbances when the reference trajectory is sinusoidal.

**Figure 5.** Comparison of the calculation time between the two NDO methods with disturbances.

**Figure 6.** Trajectory tracking errors of the three trajectory tracking methods with different parameters.

Figure 3 shows that the improved NDO-based MPC has better disturbance-rejection performance than MPC without an observer when the reference trajectory is a circle. It also can be seen from the simulations that both methods can track the reference trajectory. However, the former can track the reference trajectory smoothly, while the latter fluctuates a lot. The surge position error range of MPC without an observer for trajectory tracking control is −0.1 to 0.1 m; the sway position error range is −0.1 to 0.1 m; and the yaw angle error range is −0.05 to 0.05. Although the surge position error range of the NDO-based MPC is −0.005 to 0.005 m, the sway position error range is −0.005 to 0.005 m and the yaw angle error range is −0.005 to 0.005. Figure 4 shows the comparisons of the two trajectory tracking methods with disturbances when the reference trajectory is sinusoidal. In addition, the improved NDO-based MPC has lower tracking errors than MPC without an observer. The surge position error range of MPC without an observer for trajectory tracking control is −0.05 to 0.05 m; the sway position error range is −0.02 to 0.02 m; and the yaw angle error range is 0.00 to 0.05. Although the surge position error range of NDO-based MPC is 0.05 to 0.05 m, the sway position error range is 0.02 to 0.06 m, and the yaw angle error range is 0.000 to 0.005. Figure 5 shows that the improved NDO has better performance than the unimproved NDO in terms of the calculation time with disturbances.

We do not have a real ship for this experiment. To overcome this shortage, the ship parameters were changed when we did some comparisons between different methods to ensure the robustness of the proposed method. The USV model has uncertainty, so three methods were used for trajectory tracking for the USV with different model parameters. The three methods are unimproved NDO-based MPC, improved NDO-based MPC, and MPC without NDO. Figure 6 shows that the improved NDO and unimproved NDO effectively reduce the roughness caused by model uncertainty. In addition, the unimproved NDO has similar performance of tracking errors compared to the improved NDO-based MPC.

In addition, a comparison of the calculation time of two NDOs with disturbances is shown in Table 1. It shows the average calculation time and maximum single calculation time of the two NDOs. *IAE* = ) *<sup>t</sup>* <sup>0</sup> <sup>|</sup>*e*(*ζ*)|*d<sup>ζ</sup>* and *RMSE* = ( <sup>1</sup> *t* ) *t* <sup>0</sup> *<sup>e</sup>*2(*ζ*)*dζ*) 1/2 are used to evaluate the tracking effect and steady-state performance. The smaller the values of IAE and RMSE, the better the control performance of the scheme applied. In addition, the comprehensive performance comparisons of the position and speed tracking errors of the two methods are shown in Table 2. Table 2 shows the *IAE* and *RMSE* of the two methods.

From Table 1, the calculation time of the improved NDO is much lower than that of the traditional NDO. The average calculation time of the improved NDO is 0.0020, and it is 0.0024 for the traditional NDO. The maximum single calculation time of the improved NDO is 0.0046, and it is 0.0064 for the traditional NDO. The results show that, compared with the traditional NDO, the average calculation time of the improved NDO is decreased by 16.67%, and the maximum individual calculation time is decreased by 28.13%.


**Table 1.** Comparison of the calculation time of the two methods.

**Table 2.** Comparison of position and velocity errors of the two methods.


From Table 2, the *IAE* and *RMSE* of NDO-based MPC are lower than the MPC without an observer. NDO-based MPC effectively enhances the anti-disturbance performance of the system. MPC without observer trajectory tracking control has the characteristics of the predictive model, rolling optimization, and feedback correction, which can resist external disturbances and model mismatches to some extent.

With the data shown in Tables 1 and 2, NDO-based MPC is superior to MPC without an observer for trajectory tracking control in terms of position and speed tracking errors.

The NDO is designed to estimate the external disturbances suffered, to improve the anti-disturbance performance of the USV. Therefore, the comparison of estimated values of the NDO and actual disturbances is shown in Figure 7.

**Figure 7.** Comparison of estimated values of the NDO and actual disturbances.

Figure 7 shows the relationship between the estimated values of the NDO and the actual values of the disturbance. The NDO has good estimation performance in terms of the disturbances, including surge disturbance force, sway disturbance force, and yaw disturbance force.

### **5. Conclusions**

In this paper, an improved NDO-based MPC for trajectory tracking is proposed to guarantee the stable motion of a USV, which suffers various disturbances from the ocean wind, waves, and currents. MPC is used to optimize the system torque based on the measured position and speed state variables. Then, the NDO is designed to estimate the disturbances, and the estimated torque is compensated for in the controller. Estimation errors can converge to zero in a finite time. The simulation results show that NDObased MPC can effectively compensate for external disturbances and obtain good tracking and disturbance-rejection performance. The proposed method has a similar tracking performance to the USVs with the MPC based on unimproved NDO, but the improved NDO-based MPC is far quicker. However, due to the linearization of the model of the USV, the method only shows good performance in a near-neighbor area around the operation point. For large-difference operation points, the parameters need to be retuned.

**Author Contributions:** Conceptualization, H.F.; methodology, H.F. and S.Z.; software, W.Y.; validation, W.Y., R.C. and S.Z.; writing—original draft preparation, W.Y. and S.Z.; writing—review and editing, H.F., R.C. and S.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 52271313; Innovative Research Foundation of Ship General Performance, grant number 21822216; Fundamental Research Funds for the Central Universities, grant number XK2040021004025.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
