**1. Introduction**

Underactuated AUV is a kind of AUV which has fewer independent control inputs than the DOF to be controlled. Compared with the fully actuated AUV, it has more advantages in saving costs, reducing consumption, and improving system reliability, so it has a wide range of applications [1]. However, due to the challenges such as nonlinearity, model uncertainties, time-varying external disturbance, actuator saturation, etc., the control of the underactuated AUV becomes difficult. Currently, the precise motion control of the underactuated AUV is one of the research hotspots.

As a powerful system design and analysis tool, CT [2] has been applied in many fields; however, to our knowledge, there is currently no literature on the application of CT to underactuated AUVs. Based on this, this paper focuses on the horizontal trajectory tracking control of an underactuated AUV in the presence of model uncertainties, time-varying environmental disturbances, and actuator saturation. CT and its application in SPS [3] are used to design the controller and give quantitative analyses of various errors. In general, the contributions of this paper are as follows:


Compared with existing methods, the controller designed in this paper has a simple form and is easy to apply. In solving input saturation problems especially, there is no need to design auxiliary systems or run complex algorithms such as NN or FL.

(3) CT is applied to analyse the convergence properties of the DO and the SPS, and the upper bounds of estimation error, tracking error, and the error between the ideal controller and actual controller are given.

The structure of this paper will be arranged as follows. Section 1 introduces the research background, purpose, contributions, and structure of this paper. Section 2 reviews the literature related to the work of this paper. Section 3 establishes the mathematical model for an underactuated 3-DOF AUV in the presence of model uncertainties, external disturbances, and actuator saturation. A coordinate transformation is introduced to deal with the problem of underactuation. The DO and the saturated controller are designed in Section 4. Section 5 gives quantitative analysis of the control variable error and the tracking error. Section 6 conducts numerical simulations and analyses the performance of the controller and DO. Finally, Section 7 gives conclusions. In addition, Appendix A summarizes the abbreviations used in this paper and Appendix B briefly introduces the theoretical basis of this paper.

#### **2. Literature Review**

Underactuation is the first problem to be solved in the motion control of a USV or underactuated AUV. One of the common methods to tackle this problem is to derive the tracking error model in SF frame [4–7], and then design the control law on each controllable DOF. In [4], a 3D-path-following error model for an underactuated 5-DOF AUV was established based on virtual guidance method in SF frame, and then backstepping and sliding mode control were applied to design controllers in surge, pitch, and yaw directions. In [6], the path-following error dynamics were derived and several reduced-order ESOs were designed to estimate various disturbances. Coordinate transformation is also a common method [8–12]. In [8,9], the model of a USV was converted into a cascade system and the control problem of the USV was transformed into the stabilization analysis of two small subsystems. In [10,11], the output variables of the USV were redefined via a coordinate transformation and then the other control methods were applied to design the controller such that the new output variables could track their desired values. In addition, the system order of the underactuated AUV can be reduced by constructing an SPS [13–15]. In [13], SPT was used to decompose the full system of a 4-DOF underactuated AUV into two time-scale subsystems and then the independent controller design was carried out on each subsystem.

Disturbances are a ubiquitous challenge for the motion control of an AUV. It may come from model uncertainties, unknown system parameters, or external environment disturbances, or more likely the superposition of these factors. The estimation and compensation of disturbances are very important for AUVs' motion control. For mechanical systems including AUVs, a common practice is to combine various disturbances into total disturbances, and then design a DO to estimate it. The observer based on auxiliary variables and the ESO [16] are two common types of DO. Readers can find a detailed overview of the first type of DO in [17,18]. In [19], an auxiliary variable was introduced to design a nonlinear DO and then a backstepping finite-time sliding mode controller was constructed for the trajectory tracking control of a 5-DOF AUV. In [20], a reduced-order observer was proposed to estimate the total disturbances. ESO estimates the total disturbances as another state of the system. In [14], a high gain ESO based on SPT was designed to estimate the total disturbances. In [21], a proportional-integral velocity variable based third-order fast finite-time ESO was designed to estimate the lumped uncertainties and their first derivatives. In addition, with the continuous maturity of artificial intelligence, more and more researchers use NN [22] and FL [23] to approximate disturbances and uncertainties.

In practical applications, the problem of actuator saturation is almost unavoidable since the force/torque provided by the actuators is limited. Ignoring this problem may reduce the performance of the controller or even cause instability. At present, the auxiliary variable system [24–27], adaptive method [28], FL [29], NN [30], MPC [31,32], and DAI [33,34] are common means used to address actuator saturation. In [27], the saturation effects of rudder angle in diving control of AUVs were compensated by a modified auxiliary system with time-varying nonlinear gains. In [28], a smooth dead zone based model was designed to linearize the actuator model, so that the adaptive law could be applied to eliminate input saturation and actuator failure. In [31,32], a 3D trajectory tracking controller based on MPC was proposed such that the state and control input of the AUV did not exceed their respective physical constraints. In [33], DAI was recently applied to guarantee saturation avoidance and the techniques were applied to DC motor control for UUVs [34]. However, most of the methods used to solve input saturation are relatively complex, and some require additional design of new systems. In addition, methods such as NN and FL require strong computing power, which is a significant challenge for AUVs.

In the above literature, the Lyapunov method occupies an absolute position in the design and stability analysis of the system. However, for a complex nonlinear system, it is not easy to find a suitable Lyapunov function and prove that its first derivative is UND. In recent years, the continuous development of CT [2] provides researchers with new ideas. It studies the convergence properties of system trajectory, which is very applicable to tracking control problems. At present, CT has been widely used in many fields, such as controller and observer design [35–39], cooperative control [40,41], SPS [3,42], iterative learning control [43,44], convex optimization [45,46], and so on.

Compared with the extensive application in other fields, the application of CT in AUVs is still rare, and the model of AUVs in some rare applications is relatively simple. In [35], CT was used to design the position and velocity observer for an AUV with 1-DOF and analyse the convergence properties. In [47], after omitting the Coriolis and centripetal terms in the dynamic equation, CT was applied to design an UGES observer for a 6-DOF AUV. Unfortunately, the simulation results of the observer were not shown. In [48,49], combining the CT and the backstepping method, the author designed a speed stabilization controller and a trajectory tracking controller for a simplified AUV, and discussed the incremental stability of the AUV system. In [37], CT was applied to construct a trajectory tracking controller for an openframe AUV under the pH framework.

#### **3. AUV Modelling**

### *3.1. AUV Modelling in the Horizontal Plane*

This section establishes a horizontal 3-DOF kinematics and dynamics model for an underactuated AUV. To describe the motion of the AUV, we first define two reference frames: earth-fixed frame and body-fixed frame, as shown in Figure 1. Here, (*x*, *y*, *z*)*<sup>T</sup>* and (*ϕ*, *θ*, *ψ*)*<sup>T</sup>* denote the position vector and the attitude vector relative to the earth-fixed frame, respectively. The terms (*u*, *v*, *w*)*<sup>T</sup>* and (*p*, *q*,*r*)*<sup>T</sup>* represent the linear velocity vector and angular velocity vector with respect to the body-fixed frame, respectively.

**Figure 1.** Schematic diagram of coordinate system.

Since the motion considered here is the motion of an underactuated 3-DOF AUV in the horizontal plane, the dynamics in heave, roll, and pitch directions are all neglected and the AUV is supposed to be neutrally buoyant. The kinematics model is such that [50]

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

and the dynamics model

$$\begin{cases} & m\_{11}\dot{u} - m\_{22}\upsilon r - X\_{\text{u}}u - X\_{|u|u}|u|u &= \tau\_{u} + \tau\_{d1} \\ & m\_{22}\dot{\upsilon} + m\_{11}\iota r - Y\_{\upsilon}\upsilon - Y\_{|\upsilon s|,\upsilon}|\upsilon s.|\upsilon &= \tau\_{d2} \\ & m\_{33}\dot{\tau} + (m\_{22} - m\_{11})\iota\upsilon - N\_{\tau}r - N\_{|\upsilon|r}|r|r &= \tau\_{r} + \tau\_{d3} \end{cases} \tag{2}$$

where *m*<sup>11</sup> = *m* − *Xu*˙ , *m*<sup>22</sup> = *m* − *Yv*˙ and *m*<sup>33</sup> = *Iz* − *Nr*˙. *m* is the AUV's mass; *Iz* is the moment of inertia about the yaw rotation; *Xu*˙ , *Yv*˙ , and *Nr*˙ are the hydrodynamic added mass terms in the surge, the sway, and the yaw directions, respectively; *Xu*, *Yv*, and *Nr* are the linear damping terms, and *<sup>X</sup>*|*u*|*u*, *<sup>Y</sup>*|*vs*.|*v*, and *<sup>N</sup>*|*r*|*<sup>r</sup>* are the second-order damping terms. *<sup>τ</sup><sup>u</sup>* and *τ<sup>r</sup>* are control inputs for the surge force and the yaw torque. *τd*1, *τd*2, and *τd*<sup>3</sup> represent the total disturbances in three directions, they are composed of model uncertainties and time-varying environmental disturbances.

To facilitate subsequent design, (2) is simplified as follows:

$$\begin{cases} \begin{aligned} \dot{u} &= f\_u + m\_{11}^{-1} \tau\_u + m\_{11}^{-1} \tau\_{d1} \\ \dot{v} &= f\_v + m\_{22}^{-1} \tau\_{d2} \\ \dot{r} &= f\_r + m\_{33}^{-1} \tau\_r + m\_{33}^{-1} \tau\_{d3} \end{aligned} \tag{3}$$

with *fu* = *m*−<sup>1</sup> <sup>11</sup> (*m*22*vr* <sup>+</sup> *Xuu* <sup>+</sup> *<sup>X</sup>*|*u*|*u*|*u*|*u*), *fv* <sup>=</sup> *<sup>m</sup>*−<sup>1</sup> <sup>22</sup> (−*m*11*ur* + *Yvv* + *<sup>Y</sup>*|*vs*.|*v*|*vs*.|*v*), and *fr* = *m*−<sup>1</sup> <sup>33</sup> ((*m*<sup>11</sup> − *<sup>m</sup>*22)*uv* + *Nrr* + *<sup>N</sup>*|*r*|*r*|*r*|*r*).

This paper also considers the problem of actuator saturation. In order to deal with this problem, we design a control law *τ<sup>i</sup>* = *σ*(*τci*, *τmax*), *i* = *u*,*r*, where *τci* is the force/torque calculated by the controller, *τmax* represents the maximum force/torque that the actuator can provide, *τ<sup>i</sup>* is the actual output force/torque of the actuator, and *σ*(·) is a bounded smooth function, which satisfies

$$\begin{cases} & \sigma(0, \tau\_{\text{max}}) = 0 \\ & s\sigma(s, \tau\_{\text{max}}) > 0, \forall s \neq 0 \\ & \lim\_{s \to +\infty} \sigma(s, \tau\_{\text{max}}) = \tau\_{\text{max}}, \lim\_{s \to -\infty} \sigma(s, \tau\_{\text{max}}) = -\tau\_{\text{max}} \\ & \frac{\partial \sigma(s, \tau\_{\text{max}})}{\partial s} > 0, \forall s \in \mathbb{D}\_{s} \in \mathbb{R} \end{cases} \tag{4}$$

As we know, many functions including the Gaussian error function and hyperbolic tangent function satisfy the properties in (4). Here we choose the hyperbolic tangent function (Figure 2), so the form of the controller is *τ* = *σ*(*τc*, *τmax*) = *τmaxtanh*(*τc*/*τmax*), where *τ* = [*τ<sup>u</sup> τr*] *T*, *τ<sup>c</sup>* = [*τcu τcr*] *T*.

**Figure 2.** Schematic diagram of *tanh*(*x*).

The dynamics model considering actuator saturation is obtained:

$$\begin{cases} \begin{aligned} \dot{u} &= f\_u + m\_{11}^{-1} \sigma(\mathbf{r}\_{c\nu}, \mathbf{r}\_{\text{max}}) + m\_{11}^{-1} \mathbf{r}\_{d1} \\ \dot{v} &= f\_v + m\_{22}^{-1} \mathbf{r}\_{d2} \\ \dot{r} &= f\_r + m\_{33}^{-1} \sigma(\mathbf{r}\_{c\nu}, \mathbf{r}\_{\text{max}}) + m\_{33}^{-1} \mathbf{r}\_{d3} \end{aligned} \end{cases} \tag{5}$$

The following Assumptions are made to facilitate the subsequent design and analyses:

**Assumption 1.** *τd*1, *τd*2*, and τd*<sup>3</sup> *are all unknown and time-varying, and their first and second time derivatives are bounded.*

**Remark 1.** *The disturbances are bounded since the energy of the environmental disturbance and the state of AUVs are finite. Therefore, Assumption 1 is reasonable.*

#### *3.2. Coordinate Transformation*

To cope with the problem of underactuation, the following coordinate transformation is introduced [10,12]:

$$\begin{cases} \begin{array}{l} \chi\_{\!\!r} = \text{x} + l \cos \psi\\ y\_{\!\!r} = \text{y} + l \sin \psi \end{array} \end{cases} \tag{6}$$

where *l* > 0 is a constant, *xr* → *x*, *y* → *y* as *l* → 0. The first order derivatives of *xr* and *yr* with respect to time are given:

$$\begin{cases}
\dot{\mathbf{x}}\_r = \iota c \cos \psi - v \sin \psi - l r \sin \psi \\
\dot{\mathbf{y}}\_r = \iota s \sin \psi + v \cos \psi + l r \cos \psi
\end{cases} \tag{7}$$

Subsequently, the second derivatives of *xr* and *yr* are obtained:

$$\begin{cases}
\ddot{x}\_r &= \dot{u}\cos\psi - (\dot{v} + l\dot{r})\sin\psi - f\_x(\mathbf{z}) \\
\ddot{y}\_r &= \dot{u}\sin\psi + (\dot{v} + l\dot{r})\cos\psi + f\_y(\mathbf{z})
\end{cases} \tag{8}$$

where **z** = [*uvr ψ*] *T*,

$$\begin{cases} f\_x(\mathbf{z}) = \iota r \sin \psi + (\upsilon r + lr^2) \cos \psi\\ f\_y(\mathbf{z}) = \iota r \cos \psi - (\upsilon r + lr^2) \sin \psi \end{cases} \tag{9}$$

In the light of (5) and (8), we have

$$\begin{cases}
\dddot{\mathbf{x}}\_{r} = F\_{\mathbf{x}}(\mathbf{z}) - f\_{\mathbf{x}}(\mathbf{z}) + \Delta \mathbf{x} + m\_{11}^{-1} \cos \mathfrak{p} \sigma(\tau\_{\mathbf{c}u}, \tau\_{\mathbf{max}}) - l m\_{33}^{-1} \sin \mathfrak{p} \sigma(\tau\_{\mathbf{c}r}, \tau\_{\mathbf{max}}) \\
\ddot{\mathbf{y}}\_{r} = F\_{\mathbf{y}}(\mathbf{z}) + f\_{\mathbf{y}}(\mathbf{z}) + \Delta \mathbf{y} + m\_{11}^{-1} \sin \mathfrak{p} \sigma(\tau\_{\mathbf{c}u}, \tau\_{\mathbf{max}}) + l m\_{33}^{-1} \cos \mathfrak{p} \sigma(\tau\_{\mathbf{c}r}, \tau\_{\mathbf{max}})
\end{cases} \tag{10}$$

where

$$\begin{cases} F\_{\mathbf{x}}(\mathbf{z}) &= f\_{\mathbf{l}} \cos \psi - (f\_{\mathbf{v}} + lf\_{\mathbf{r}}) \sin \psi \\ \Delta \mathbf{x} &= m\_{11}^{-1} \pi\_{d1} \cos \psi - (m\_{22}^{-1} \pi\_{d2} + lm\_{33}^{-1} \pi\_{d3}) \sin \psi \end{cases} \tag{11}$$

$$\begin{cases} \begin{array}{rcl} \mathbf{F}\_{\mathbf{Y}}(\mathbf{z}) &= f\_{\mathbf{u}} \sin \psi + (f\_{\mathbf{v}} + lf\_{\mathbf{f}}) \cos \psi\\ \Delta y &= m\_{11}^{-1} \mathbf{r}\_{d1} \sin \psi + (m\_{22}^{-1} \mathbf{r}\_{d2} + lm\_{33}^{-1} \mathbf{r}\_{d3}) \cos \psi \end{array} \tag{12}$$

Let *η*<sup>1</sup> = [*xr yr*] *<sup>T</sup>* be the reference trajectory and *<sup>η</sup>*<sup>2</sup> = [*x*˙*<sup>r</sup> <sup>y</sup>*˙*r*] *<sup>T</sup>*, (10) can be rewritten as follows:

$$\begin{cases}
\dot{\eta}\_1 = \eta\_2 \\
\dot{\eta}\_2 = \Phi(\mathbf{z}) + \Delta + \mathbf{g}(\boldsymbol{\psi})\sigma(\boldsymbol{\tau}\_{c\prime}, \boldsymbol{\tau}\_{\text{max}})
\end{cases} \tag{13}$$

where **<sup>Φ</sup>**(**z**)=[*Fx*(**z**) − *fx*(**z**) *Fy*(**z**) + *fy*(**z**)]*T*, **<sup>Δ</sup>** = [Δ*<sup>x</sup>* <sup>Δ</sup>*y*] *<sup>T</sup>* is the combined disturbance vector and

$$\mathbf{g}(\boldsymbol{\psi}) = \begin{bmatrix} m\_{11}^{-1}\cos\psi & -lm\_{33}^{-1}\sin\psi\\ m\_{11}^{-1}\sin\psi & lm\_{33}^{-1}\cos\psi \end{bmatrix} \tag{14}$$

It is obvious that the matrix *g*(*ψ*) is nonsingular for any *ψ*. The control objective is to design a controller *τ<sup>c</sup>* such that the reference trajectory *η*<sup>1</sup> = [*xr yr*] *<sup>T</sup>* converges to the desired trajectory *η<sup>d</sup>* = [*xd yd*] *<sup>T</sup>* asymptotically, and the error between the actual trajectory *η* = [*x y*] and the desired trajectory is small enough. In addition, considering that the disturbance **Δ** is unknown, a DO is needed.

**Assumption 2.** *According to Assumption 1, (11) and (12), the disturbance vector* **Δ** *and its first derivative are also bounded. For the convenience of subsequent analysis, we assume that* **Δ**˙ ≤ <sup>Δ</sup>*.*

**Assumption 3.** *The desired trajectory η<sup>d</sup>* = [*xd yd*] *<sup>T</sup> and its first two time derivatives <sup>η</sup>*˙ *<sup>d</sup>, <sup>η</sup>*¨ *<sup>d</sup> are bounded. Furthermore, the desired yaw angle <sup>ψ</sup><sup>d</sup>* <sup>=</sup> *arctan <sup>y</sup>*˙*d*(*t*) *<sup>x</sup>*˙*d*(*t*) *is also bounded.*

#### **4. DO and Saturated Controller Design**

*4.1. DO Design*

Here, we introduce an auxiliary variable to design a DO to estimate **Δ**, and its expression is as follows:

$$\begin{cases} \dot{\Delta} &= \mathfrak{F} + K\_1 \mathfrak{v}\_2 \\ \dot{\mathfrak{F}} &= K\_1(-\Phi(\mathbf{z}) - \mathbf{g}\psi r(\mathbf{r}\_{c\prime}, \mathbf{r}\_{\text{max}})) - K\_1(\mathfrak{F} + K\_1 \mathfrak{v}\_2) \end{cases} \tag{15}$$

where **Δ**ˆ is the estimated value of **Δ**, *ξ* is the auxiliary variable and *K*<sup>1</sup> > 0 is the observer gain. Define the estimation error *<sup>e</sup>*<sup>1</sup> <sup>=</sup> **<sup>Δ</sup>**<sup>ˆ</sup> <sup>−</sup> **<sup>Δ</sup>**. According to (13) and (15), the dynamic of *<sup>e</sup>*<sup>1</sup> is obtained as follows:

$$\begin{split} \dot{\mathbf{e}}\_{1} &= \dot{\Delta} - \dot{\Delta} \\ &= K\_{1}(-\Phi(\mathbf{z}) - f(\mathbf{z}) - \mathbf{g}\,\psi r(\pi\_{\mathbf{c}}, \pi\_{\max})) - K\_{1}\hat{\Delta} + K\_{1}\dot{\eta}\_{2} - \dot{\Delta} \\ &= -K\_{1}\mathbf{e}\_{1} - \dot{\Delta} \end{split} \tag{16}$$

according to Assumption 2, the dynamic of *e*<sup>1</sup> can be regarded as the combination of the nominal dynamic *<sup>e</sup>*˙ <sup>0</sup> <sup>=</sup> <sup>−</sup>*K*1*e*<sup>0</sup> and the bounded disturbance **<sup>Δ</sup>**˙ . Obviously, the nominal dynamic is contracting with respect to *e*<sup>0</sup> with a contracting rate *λ*<sup>1</sup> = *K*<sup>1</sup> and a transformation matrix **Θ**1.

According to triangular inequality, we obtain *e*1≤*e*<sup>1</sup> − *e*0 + *e*0. Since *e*˙ <sup>0</sup> = −*K*1*e*<sup>0</sup> is contracting, there is *e*0(*t*)≤*e*0(0) and *e*0(0) is the initial value of *e*0. Then, applying the robust properties of contracting system, we can obtain

$$\begin{aligned} \|\mathbf{e}\_1\| &\le \|\mathbf{e}\_1 - \mathbf{e}\_0\| + \|\mathbf{e}\_0\|\\ &\le \chi\_1 \|\mathbf{e}\_1(0) - \mathbf{e}\_0(0)\| \varepsilon x p^{-\lambda\_1 t} + \frac{\chi\_1 \overline{\Delta}}{\lambda\_1} + \|\mathbf{e}\_0(0)\| \\ &= v\_1 \end{aligned} \tag{17}$$

where *χ*<sup>1</sup> is the upper bound of the condition number of **Θ**1.

We can see that *<sup>v</sup>*<sup>1</sup> consists of three parts. The first part *<sup>χ</sup>*1*e*1(0) <sup>−</sup> *<sup>e</sup>*0(0)*exp*−*λ*1*<sup>t</sup>* will converge to zero exponentially, so the size of *v*<sup>1</sup> will ultimately depend on the last two items. Obviously, *e*1 can be infinitely close to zero by selecting the appropriate observer gain and the initial estimation values of the disturbances.

#### *4.2. Saturated Controller Design*

In this subsection, we apply CT and SPT to design a saturated controller. Firstly, we define a tracking error *e*<sup>2</sup> = *η*1(*t*) − *ηd*(*t*) and then construct a new variable based on *e*2:

$$\mathbf{S} = \dot{\mathbf{e}}\mathbf{2} + K\mathbf{e}\mathbf{e}\_2 \tag{18}$$

here *K*<sup>2</sup> > 0 is the gain. From the definition of *S*, we know that

$$S \to 0 \Rightarrow \mathfrak{e}\_2 \to 0 \tag{19}$$

In the light of (13) and (18), we have

$$\mathbf{S} = \mathbf{F}(\mathbf{z}, \boldsymbol{\tau}\_c, \boldsymbol{\Delta}, \mathbf{e}\_2) = \boldsymbol{\Phi}(\mathbf{z}) + \boldsymbol{\Delta} + \mathbf{g}(\boldsymbol{\psi})\boldsymbol{\sigma}(\boldsymbol{\tau}\_c, \boldsymbol{\tau}\_{\max}) - \ddot{\boldsymbol{\eta}}\_d + K\_2 \dot{\mathbf{e}}\_2 \tag{20}$$

The present objective is to design a controller *τ<sup>c</sup>* to make *S* converge to 0. Applying SPT, we construct a dynamic process for *τc*:

$$\mu \dot{\boldsymbol{\pi}}\_{\varepsilon} = \mathbf{H}(\mathbf{z}, \boldsymbol{\pi}\_{\varepsilon}, \mathbf{A}, \mathbf{e}\_{2}, \boldsymbol{\mu}) = -\mathbf{g}^{T}(\boldsymbol{\psi})(\mathbf{K}\_{3}\mathbf{S} + \boldsymbol{\Phi}(\mathbf{z}) + \boldsymbol{\Delta} + \mathbf{g}(\boldsymbol{\psi})\boldsymbol{\sigma}(\boldsymbol{\pi}\_{\varepsilon}, \boldsymbol{\pi}\_{\max}) - \ddot{\boldsymbol{\eta}}\_{d} + \mathbf{K}\_{2}\dot{\mathbf{e}}\_{2}) \tag{21}$$

where 0 < *μ* 1 is the singular perturbation parameter and *K*<sup>3</sup> > 0 is the controller gain. Equations (20) and (21) constitute a standard SPS [51], where (20) is the slow subsystem,

and (21) is the fast subsystem. Since *μ* 1, the dynamic process of (21) is much faster than (20).

Now, we apply CT to analyse the properties of the controller. By solving the partial derivative of *τc*, the Jacobian of (21) is

$$J\_b = -\frac{1}{\mu} \mathbf{g}^T(\boldsymbol{\psi}) \mathbf{g}(\boldsymbol{\psi}) \begin{bmatrix} \frac{\partial \boldsymbol{r}(\cdot)}{\partial \boldsymbol{\tau}\_\mu} \\\\ \frac{\partial \boldsymbol{r}(\cdot)}{\partial \boldsymbol{\tau}\_\tau} \end{bmatrix} \tag{22}$$

According to the properties of *σ*(·) and *g*(*ψ*), it is obvious that *J<sup>b</sup>* is UND. Therefore, for any **z**, **Δ**, *e*2, and *μ*, *H*(**z**, *τc*, **Δ**, *e*2, *μ*) is partially contracting with respect to *τc*, we assume that the contracting rate is *<sup>λ</sup>*<sup>2</sup> *<sup>μ</sup>* and the transformation matrix is **Θ**2. Then, according to the results in [3], the algebraic equation *H*(**z**, *τc*, **Δ**, *e*2, 0) = 0 can be equivalently written as *τ<sup>d</sup>* = *ϑ*(**z**, **Δ**, *e*2), i.e., there is a unique, global mapping between *τ<sup>c</sup>* and **z**, **Δ**, *e*, where *τ<sup>d</sup>* is the root of the above algebraic equation, which is also called the quasi-steady state of the fast subsystem (21). According to the properties of contracting system, *τ<sup>c</sup>* converges to its quasisteady state *τ<sup>d</sup>* exponentially. Now, solving the algebraic equation *H*(**z**, *τd*, **Δ**, *e*2, 0) = 0, we obtain

$$\sigma(\mathbf{r}\_{d\prime}, \mathbf{r}\_{\text{max}}) = \mathbf{g}^{-1}(\boldsymbol{\psi})(-K\_{\rm 3}\mathbf{S} - \boldsymbol{\Phi}(\mathbf{z}) - \boldsymbol{\Delta} + \ddot{\eta}\_{d} - K\_{\rm 2}\dot{\mathbf{e}}\_{2})\tag{23}$$

according to SPT, the slow subsystem can be simplified via the quasi-steady state of the fast subsystem, therefore, bring (23) into (20), and the simplified slow subsystem is obtained:

$$\mathbf{S} = -\mathsf{K}\_3 \mathbf{S} \tag{24}$$

since *K*<sup>3</sup> > 0, the dynamic of *S* is contracting with a transformation metric *I* and a contracting rate *λ<sup>s</sup>* = *K*3. Therefore, *S* converges to 0 exponentially. Based on the above analysis, we know that

$$
\mu \to 0 \Rightarrow \pi\_c \to \pi\_d \Rightarrow \mathcal{S} \to 0 \tag{25}
$$

Therefore, when *μ* is small enough and *K*<sup>3</sup> is large enough, *τ<sup>c</sup>* → *τ<sup>d</sup>* and *S* → 0 quickly. At the same time, due to the existence of *σ*(·), the control input will not exceed the limit of the actuator. Compared with the literature in the introduction, the method proposed in this paper is simple in form and convenient in application when dealing with actuator saturation.

Since it has been proven in Section 3.1 that **Δ**ˆ approaches **Δ** infinitely, we replace **Δ** in (21) with its estimated value, and the practical controller can be obtained:

$$\begin{array}{rcl} \mu \dot{\boldsymbol{\tau}}\_{c} &=& -\mathbf{g}^{T}(\boldsymbol{\psi}) (\mathbf{K}\_{3}\mathbf{S} + \boldsymbol{\Phi}(\mathbf{z}) + \boldsymbol{\hat{\Lambda}} + \mathbf{g}(\boldsymbol{\psi}) \boldsymbol{\sigma}(\boldsymbol{\tau}\_{c}, \boldsymbol{\tau}\_{\max}) - \ddot{\boldsymbol{\eta}}\_{d} + \mathbf{K}\_{2} \dot{\boldsymbol{\varepsilon}}\_{2} \\ \boldsymbol{\tau} &=& \boldsymbol{\tau}\_{\max} \tanh(\boldsymbol{\tau}\_{c}/\boldsymbol{\tau}\_{\max}) \end{array} \tag{26}$$

In order to ensure the fast convergence of *τ<sup>c</sup>* to *τ<sup>d</sup>* and *e*<sup>2</sup> to 0, *K*2, *K*<sup>3</sup> and *μ* should be reasonably selected. At the same time, *K*<sup>1</sup> should be large enough to ensure the rapid convergence of the estimation error. Of course, due to the complexity of the system, each parameter has an impact on the performance of the controller. After many simulations, we find that three parameters *l*, *μ*, and *K*<sup>2</sup> have a significant impact on the performance of the controller. For *l*, when *l* is too small, the controller may diverge, while too large *l* will

reduce the tracking accuracy. Reducing *μ* can improve the convergence speed, but it also makes the changes of the control input less smooth. For *μ* , the smaller the *μ* , the faster the convergence speed will be, but it will lead to drastic changes in the control input, while an excessively large *μ* will reduce the convergence speed. At the same time, the larger the *K*<sup>2</sup> , the faster the convergence speed, but this will cause oscillation. In addition, *K*<sup>1</sup> mainly affects the performance of DO, and *K*<sup>3</sup> has a limited impact on the results. Based on the above analysis, it is necessary to make a balance between tracking accuracy, convergence speed, and smooth changes in control input when selecting parameters.

#### **5. Error Analysis**

#### *5.1. Error Analysis of Control Variables*

As analysed above, when *τ<sup>c</sup>* converges exponentially to its quasi-steady state *τd*, the contracting reduced-order slow subsystem will be obtained. However, there is always an error between *τ<sup>c</sup>* and *τd*, because *μ* can only be be as small as possible and cannot be equal to zero. Here, we make a quantitative analysis of the error between *τ<sup>c</sup>* and *τd*.

Defining the control variable error *τ<sup>e</sup>* = *τ<sup>c</sup>* − *τd*. According to SPT, the fast boundary layer dynamics in a new time scale *t*<sup>1</sup> = *<sup>t</sup> <sup>μ</sup>* can be derived as

$$\frac{d\mathbf{\dot{\tau}\_{\varepsilon}}}{dt\_{1}} = H(\mathbf{z}, \mathbf{\dot{\tau}\_{\varepsilon}} + \mathbf{\tau}\_{d\prime}\hat{\Delta} - \mathbf{e}\_{1}, \mathbf{e}\_{2\prime}\boldsymbol{\mu}) - \mu\mathbf{\dot{\tau}\_{d}}\tag{27}$$

The unperturbed error dynamic of (27) is

$$\frac{d\mathbf{\tau}\_{\varepsilon}}{dt\_{1}} = H(\mathbf{z}, \mathbf{\tau}\_{\varepsilon} + \mathbf{\tau}\_{d}, \mathbf{\hat{A}} - \mathbf{e}\_{1}, \mathbf{e}\_{2}, \mu) \tag{28}$$

Because the Jacobian *J<sup>b</sup>* in (22) is UND, it implies that (28) is also partially contracting with respect to *τe*, and the contracting rate is *<sup>λ</sup>*<sup>2</sup> *<sup>μ</sup>* and the transformation matrix is denoted as **Θ**<sup>2</sup> with a supermum condition number *χ*2. Assume that *τ*˙ *<sup>d</sup>* is Lipschitz continuous in *τ<sup>e</sup>* and *e*<sup>1</sup> [42], i.e.,

$$\|\|\pi\_d(\cdot)\|\| \le c\_1 \|\|\pi\_c\|\| + c\_2 \|\|\mathfrak{e}\_1\|\| + c\_3 \tag{29}$$

where *c*1, *c*2, *c*<sup>3</sup> are all positive constants.

According to the robust properties of contracting system and ignoring the initial value of *τe*,

$$\|\boldsymbol{\pi}\_{\varepsilon}(t)\| = \|\boldsymbol{\pi}\_{\varepsilon} - \boldsymbol{\pi}\_{d}\| \leq \chi\_{2} \|\boldsymbol{\pi}\_{\varepsilon}(0) - \boldsymbol{\pi}\_{d}(0)\| \exp^{-\frac{\lambda\_{2}}{\mu}t} + \frac{\mu \chi\_{2} (c\_{1} \|\boldsymbol{\pi}\_{\varepsilon}\| + c\_{2} \|\boldsymbol{\varepsilon}\_{1}\| + c\_{3})}{\lambda\_{2}} \tag{30}$$

then, we can obtain

$$\|\|\pi\_{\varepsilon}(t)\|\| \leq \frac{\lambda\_{2}\chi\_{2}}{\lambda\_{2} - \mu\chi\_{2}c\_{1}} \|\|\pi\_{\varepsilon}(0) - \pi\_{d}(0)\|\|\exp^{-\frac{\lambda\_{1}}{\mu^{2}}t} + \frac{\mu\chi\_{2}(c\_{2}\upsilon\_{1} + c\_{3})}{\lambda\_{2} - \mu\chi\_{2}c\_{1}}\tag{31}$$

it can be observed that *τ<sup>e</sup>* depends not only on *μ*, but also on *λ*<sup>2</sup> and *χ*2, which are related to the values of observer gain and controller gain. Therefore, *τ<sup>e</sup>* can be very small by reasonably selecting parameters.

#### *5.2. Analysis of Tracking Error*

The existence of *e*<sup>1</sup> and *τ<sup>e</sup>* causes *S* to fail to converge to zero. In this subsection, a quantitative analysis of *S* is conducted. By performing some simple transformations on (20), we obtain a new dynamic form of *S*

$$\begin{aligned} \mathbf{S} &= \Phi(\mathbf{z}) + \mathbf{g}(\boldsymbol{\psi})\sigma(\boldsymbol{\tau}\_{d}, \boldsymbol{\tau}\_{\max}) + \boldsymbol{\Delta} - \ddot{\boldsymbol{\eta}}\_{d} + K\_{2}\dot{\mathbf{e}}\_{2} + \\ &\Phi(\mathbf{z}) + \mathbf{g}(\boldsymbol{\psi})\sigma(\boldsymbol{\tau}\_{c}, \boldsymbol{\tau}\_{\max}) + \boldsymbol{\Delta} - \ddot{\boldsymbol{\eta}}\_{d} + K\_{2}\dot{\mathbf{e}}\_{2} - \\ &(\Phi(\mathbf{z}) + \mathbf{g}(\boldsymbol{\psi})\sigma(\boldsymbol{\tau}\_{d}, \boldsymbol{\tau}\_{\max}) + \boldsymbol{\Delta} - \ddot{\boldsymbol{\eta}}\_{d} + K\_{2}\dot{\mathbf{e}}\_{2}) \end{aligned} \tag{32}$$

By substituting (23) in (32), we obtain

$$\mathcal{S} = -K\_3 \mathcal{S} + \mathbf{g}(\psi)\sigma(\mathbf{r}\_{d\prime}, \mathbf{r}\_{\text{max}})(\sigma(\mathbf{r}\_{c\prime}, \mathbf{r}\_{\text{max}}) - \sigma(\mathbf{r}\_{d\prime}, \mathbf{r}\_{\text{max}})) \tag{33}$$

Similarly, the dynamic of *S* can be regarded as the coupling of the contracting nominal dynamic *<sup>S</sup>*˙ <sup>0</sup> <sup>=</sup> <sup>−</sup>*K*3*S*<sup>0</sup> and the bounded disturbance *<sup>g</sup>*(*ψ*)(*σ*(*τc*, *<sup>τ</sup>max*) <sup>−</sup> *<sup>σ</sup>*(*τd*, *<sup>τ</sup>max*)). Suppose the contracting rate is *λ*<sup>3</sup> and the transformation matrix is **Θ**3. For the convenience of analysis, it is further assumed that *g*(*ψ*)(*σ*(*τc*, *τmax*) − *σ*(*τd*, *τmax*)) is Lipschitz continuous in *τ<sup>e</sup>* with constant *L*, i.e.,

$$\left\|\left\|\mathbf{g}\left(\psi\right)\left(\sigma\left(\boldsymbol{\tau}\_{\mathcal{E}},\boldsymbol{\tau}\_{\max}\right)-\sigma\left(\boldsymbol{\tau}\_{\mathcal{d}},\boldsymbol{\tau}\_{\max}\right)\right)\right\|\leq L\left\|\left\|\boldsymbol{\tau}\_{\mathcal{E}}-\boldsymbol{\tau}\_{\mathcal{d}}\right\|\right\|=L\left\|\left\|\boldsymbol{\tau}\_{\mathcal{E}}\right\|\right\|\tag{34}$$

Continuing to apply the triangular inequality and the robustness property of the contracting system, and we can obtain:

$$\begin{aligned} \|\mathbf{S}\| &\le \|\mathbf{S} - \mathbf{S}\_0\| + \|\mathbf{S}\_0(0)\|\\ &\le \chi\_3 \|\mathbf{S}(0) - \mathbf{S}\_0(0)\| \|\exp^{-\lambda\_3 t} + \frac{\chi\_3 L \|\mathbf{r}\_\varepsilon\|}{\lambda\_3} + \|\mathbf{S}\_0(0)\| \end{aligned} \tag{35}$$

where *χ*<sup>3</sup> is the upper bound of the condition number of **Θ**<sup>3</sup> and *S*0(0) is the initial value of *S*0. If we replace *τ<sup>e</sup>* in (35) with (31), a more detailed expression of *S* can be obtained.

From the above analysis, we can see that the estimation error and tracking error can finally converge to a small range by reasonably selecting the controller gain, the observer gain, and the singular perturbation parameters.

#### **6. Numerical Simulations**

In this section, based on Matlab Simulink, numerical simulations on the Remus AUV [52] are conducted to illustrate the effectiveness of the proposed method; the parameters in (2) are given as follows: *m* = 30.58, *Iz* = 3.45, *Xu*˙ = −0.93, *Yv*˙ = −35.5, *Nr*˙ = −4.88, *Xu* = −13.5,*Yv* = −66.6, *Nr* = −6.87, *<sup>X</sup>*|*u*|*<sup>u</sup>* = −1.62, *<sup>Y</sup>*|*v*|*<sup>v</sup>* = −131 and *<sup>N</sup>*|*r*|*<sup>r</sup>* = −188.

In order to highlight the advantages of this method, the simulation results obtained based on the methods in [53] are compared with the results in this paper. In [53], the sliding mode control method was applied to deal with the trajectory tracking problem of an underactuated AUV by introducing a first-order sliding surface in terms of surge tracking errors and a second-order surface in terms of lateral motion tracking errors. The sliding surfaces are defined as

$$\begin{cases} \ S\_1 = \mathfrak{e}\_{\mathfrak{u}} + \lambda\_1 \int\_0^t \mathfrak{e}\_{\mathfrak{u}}(\tau) d\tau\\ \ S\_2 = \dot{\mathfrak{e}}\_{\overline{\mathfrak{v}}} + \lambda\_3 \mathfrak{e}\_{\overline{\mathfrak{v}}} + \lambda\_2 \int\_0^t \mathfrak{e}\_{\overline{\mathfrak{v}}}(\tau) d\tau \end{cases} \tag{36}$$

where *λ*1, *λ*2, *λ*<sup>3</sup> > 0. *eu* = *u* − *ud* and *ev* = *v* − *vd* are tracking errors for surge and sway velocity, respectively. *ud* and *vd* are desired surge velocity and sway velocity, and they are defined as follows:

$$
\begin{bmatrix} u\_d \\ v\_d \end{bmatrix} = \begin{bmatrix} \cos\psi & \sin\psi \\ -\sin\psi & \cos\psi \end{bmatrix} \begin{bmatrix} \dot{\mathbf{x}}\_d + l\_x \tanh\left(-\frac{k\_x}{l\_x}\mathbf{x}\_\varepsilon\right) \\ \dot{y}\_d + l\_y \tanh\left(-\frac{k\_y}{l\_y}y\_\varepsilon\right) \end{bmatrix} \tag{37}
$$

where *kx*, *ky* > 0 are controller gains and *lx*, *ly* > 0 are saturation constants. *xe* = *x* − *xd* and *ye* = *y* − *yd* are position tracking errors.

Then, the control law is given by

$$\begin{cases} \ \mathfrak{r}\_{\mathfrak{u}} = \mathfrak{r}\_{\mathfrak{u},\mathfrak{c}\eta} + \mathfrak{r}\_{\mathfrak{u},\mathfrak{sw}}\\ \ \mathfrak{r}\_{\mathfrak{u}} = \mathfrak{r}\_{\mathfrak{r},\mathfrak{c}\eta} + \mathfrak{r}\_{\mathfrak{r},\mathfrak{sw}} \end{cases} \tag{38}$$

where

$$\begin{aligned} \tau\_{xy,w} &= -m\_{22}\tau r - X\_{u}u - X\_{|u|}|u| |u| + m\_{11}\bar{u}\_d - m\_{11}\lambda\_1 \varepsilon\_u \\ \tau\_{x,w} &= m\_{11}(-k\_1 S\_1 - W\_l \text{sign}(S\_1) \\ \tau\_{r,w} &= -(m\_{11} - m\_{22})uv - N\_r r - N\_{|r|}|r|r + (m\_{11}m\_{22}^{-1}\bar{u}r - m\_{22}^{-1}(Y\_v + 2\varepsilon \text{sign}(v)Y\_{[v]p}v)v)/b \\ &+ (\Gamma - \lambda\_3 \varepsilon\_v - \lambda\_2 \varepsilon\_v)/b \\ \tau\_{r,w} &= (-k\_2 S\_2 - W\_l \text{sign}(S\_2))/b \\ \tau\_{\bar{x},w} &= (-k\_2 S\_2 - W\_l \text{sign}(S\_2))/b \\ \bar{b} &= \bar{m}\_{33}^{-1}(\bar{u}\_d - m\_{11}m\_{22}^{-1}) \\ \Gamma &= -\bar{\mathcal{X}}\_d \bar{\sigma} \sin \psi + \bar{\mathcal{Y}}\_d \cos \psi - \bar{\mathcal{X}}\_d \cos \psi - \bar{y}\_d r \sin \psi - \bar{u}\_d r + \gamma\_1 r \cos \psi + \gamma\_2 r \sin \psi \\ &+ \gamma\_1 \sin \psi \theta - \gamma\_2 \varepsilon\_\theta \omega \psi \\ \gamma\_1 &= k\_x x\_c \mathrm{sech}^2(-\frac{k\_x}{l\_x} x\_c) \\ \gamma\_2 &= k\_y \bar{y}\_c \mathrm{sech}^2(-\frac{k\_y}{l\_y} y\_c) \end{aligned}$$

The first derivatives of the desired velocity is given by

$$
\begin{bmatrix} \dot{u}\_d \\ \dot{v}\_d \end{bmatrix} = r \begin{bmatrix} -\sin\psi & \cos\psi \\ -\cos\psi & -\sin\psi \end{bmatrix} \begin{bmatrix} \dot{x}\_d + l\_x \tanh\left(-\frac{k\_x}{l\_x}x\_\ell\right) \\ \dot{y}\_d + l\_y \tanh\left(-\frac{k\_y}{l\_y}y\_\ell\right) \end{bmatrix} + \begin{bmatrix} \cos\psi & \sin\psi \\ -\sin\psi & \cos\psi \end{bmatrix} \begin{bmatrix} \dot{x}\_d - \gamma\_1 \\ \dot{y}\_d - \gamma\_2 \end{bmatrix} \tag{39}
$$

There are two types of trajectories that need to be tracked by the AUV, one is a sinusoidal curve, and the other is a combination of straight lines and circles.

Trajectory 1: Sinudoidal curve

$$\begin{cases} \begin{array}{l} \mathfrak{x}\_d(t) = 0.5t \\ y\_d(t) = 20 \cos(0.02t) \end{array} \end{cases}$$

The total simulation time is 200*π s*. To facilitate the distinction, we mark the simulation based on the method in this paper as Case 1, while the simulation based on the literature [53] is marked as Case 2. The various settings for Case 1 and Case 2 simulations have been given in Table 1.



Trajectory 2: A combination of straight lines and circles

$$\begin{cases} \begin{aligned} x\_d(t) &= 0.4t, y\_d(t) = 20, t < 100 \\ x\_d(t) &= 40 + 20\cos(0.02t + 1.5\pi - 2), y\_d(t) = 40 + 20\sin(0.02t + 1.5\pi / - 2), t \ge 100 \end{aligned} \end{cases}$$

The total simulation time is 100*π* + 100 *s*. Similarly, two simulations are called Case 3 and Case 4, respectively. The various settings required for them are shown in Table 2.



The following disturbances are considered in these simulations:

$$\begin{cases} \tau\_{d1} = 0.1f\_u - 3 + 2\sin(0.2t) + 2\varepsilon(5) \\ \tau\_{d2} = 0.2f\_v + 5 + 0.2\sin(0.2t) + 5\varepsilon(5) \\ \tau\_{d3} = 0.15f\_r + 2 + \sin(0.1t) + 3\varepsilon(5) \end{cases} \tag{40}$$

here, *ε*(5) is the zero-mean white noise with power intensity of 5%. Specifically, the first term denotes degrees of model uncertainties. The second, third, and forth terms, respectively, account for the constant, periodic unknown disturbances, and Gaussian white noise.

Simulation results for Case 1 and Case 2 are shown in Figures 3–6. Figure 3 illustrates the trajectory of the AUV. It can be seen that both the method in this paper and the sliding mode control method in the literature have high tracking accuracy. The reference trajectory in Case 1 and the actual trajectory in Case 2 all converge to the desired trajectory. Of course, the error between the actual trajectory and the desired trajectory in Case 1 is relatively large, which can be more clearly seen from Figure 4.

**Figure 3.** The comparison of trajectories in Case 1 and Case 2.

**Figure 4.** The tracking errors in Case 1 and Case 2: (**a**) position error, (**b**) heading angle error.

**Figure 5.** The estimation results of DO in Case 1.

**Figure 6.** The control inputs in Case 1 and Case 2.

Figure 4 shows the trajectory tracking error and heading angle error, where *<sup>η</sup>e*<sup>1</sup> = /(*<sup>x</sup>* − *xd*)<sup>2</sup> + (*<sup>y</sup>* − *yd*)<sup>2</sup> and *<sup>η</sup>e*<sup>2</sup> = /(*xr* − *xd*)<sup>2</sup> + (*yr* − *yd*)<sup>2</sup> denote the trajectory tracking error, *ψ<sup>e</sup>* = *ψ* − *ψ<sup>d</sup>* represents the heading angle error. It can be seen from Figure 4a that the trajectory error *ηe*<sup>2</sup> in Case 1 converges after about 200s, and it is almost 0, indicating that the reference state defined via coordinate transformation can accurately be traced to the desired state. However, *ηe*<sup>1</sup> in Case 1 does not converge to 0 but remains around 0.35, which is just the value of *l*. Reviewing Equation (6), we know that there is an inherent error between the reference state and the actual state. The results here just confirm this phenomenon. A very natural idea is that the value of *l* should be as small as possible to reduce *ηe*1. However, we find in the simulation that too small *l* will lead to slower convergence and even divergence. Therefore, the value of *l* should be balanced between tracking accuracy and convergence speed. In contrast, the sliding mode control methods in the literature have high tracking accuracy and convergence speed regardless of trajectory tracking error shown in Figure 4a or heading angle error shown in Figure 4b, which is a disadvantage of the method in this paper.

Figure 5 provides the results of the DO in Case 1. It reveals that the unknown disturbance including model uncertainties and time-varying environmental disturbance can be accurately estimated by the DO designed in this paper.

Figure 6 shows the control inputs of surge force and yaw torque in Case 1 and Case 2. It can be seen that the method in this paper considers the problem of input saturation, so neither *τ<sup>u</sup>* nor *τ<sup>r</sup>* exceed the limitation of the actuator, and the control input changes smoothly, which is conducive to the stable operation of the actuator. However, the methods in the literature do not solve the problem of input saturation. We can see that at the initial stage of simulation, both *τ<sup>u</sup>* and *τ<sup>r</sup>* are large, which can easily exceed the limitation of the actuator. Moreover, due to the use of symbolic functions in the controller, there is significant chattering in the control input curve, which is detrimental to the stable operation of the actuator.

In general, the sliding mode control methods in the literature excel in tracking accuracy and convergence speed. The advantage of the method in this paper lies in the simple form of the controller, and it easily solves the problem of input saturation, ensuring stable changes in control inputs and stable operation of the actuator. Of course, improving tracking accuracy and convergence speed is the work we must do in the next stage.

Simulation results for Case 3 and Case 4 are shown in Figures 7–10. It can be seen that the conclusions obtained from the analysis of Case 1 and Case 2 are still applicable here. Although the desired trajectory becomes more complex, the AUV still tracks the desired trajectory accurately. At the same time, even if the input saturation is more complex, the saturated controller proposed in this paper enables the control input in each direction to not exceed the limit of the actuator. Similarly, it can be seen that the performance of the sliding mode control method in the literature in Case 4 is the same as in Case 2.

**Figure 7.** The comparison of trajectories in Case 3 and Case 4.

**Figure 8.** The tracking errors in Case 3 and Case 4: (**a**) position error, (**b**) heading angle error.

**Figure 9.** The estimation results of DO in Case 3.

**Figure 10.** The control inputs in Case 3 and Case 4.

In all of above simulation results, the proposed controller works well for the horizontal trajectory tracking of underactuated AUVs in the presence of unknown internal and external disturbances. Therefore, we reach a conclusion that the validity and efficacy of the proposed method and proposed control scheme are sufficiently demonstrated.

#### **7. Conclusions**

This paper focusses on the horizontal trajectory tracking control of a 3-DOF underactuated AUV in the face of model uncertainties, time-varying external disturbance, and input saturation. A coordinate transformation is introduced to tackle the problem of underactuation and a DO is designed to estimate the total unknown disturbance. Applying CT and SPT, we design a saturation controller and perform quantitative analysis for the estimation error and the tracking error. Simulation results show that the controller proposed in this paper makes the AUV track the desired trajectory well and avoid the problem of input saturation. Of course, compared to the methods in the literature, the method in this paper still needs to be improved in terms of tracking accuracy and convergence speed. Therefore, the next research direction is to improve the tracking accuracy and convergence speed on the basis of this paper, while expanding the research results to the 3D trajectory tracking of underactuated AUVs.

**Author Contributions:** C.M.: Conceptualization, Methodology, Software, Validation, Investigation, Writing—original draft; J.J.: Methodology, Validation, Writing—Review and Editing; T.Z.: Writing— Review and Editing, Supervision; S.W.: Writing—Review and Editing; D.J.: Conceptualization, Resources, Writing—Review and Editing, Supervision, Project administration. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

**Informed Consent Statement:** No 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; in the collection, analysis, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

Due to the use of many abbreviations in this paper, for the sake of standardization, they are summarized in the Table A1 according to the order in which they appear in the paper.



#### **Appendix B**

*Appendix B.1. Contraction Theory*

For a nonlinear system

$$\dot{\mathbf{x}} = f(\mathbf{x}, t), \ x(t\_0) = \mathbf{x}\_0, \forall t \ge t\_0 \ge 0 \tag{A1}$$

where *<sup>x</sup>*(*t*) <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* is the state vector, *<sup>t</sup>* is the time and *<sup>f</sup>* : <sup>R</sup> <sup>×</sup> <sup>R</sup>*t*≥*t*<sup>0</sup> <sup>→</sup> <sup>R</sup>*<sup>n</sup>* is a nonlinear smooth function, meaning that all required derivatives and partial derivatives exist and are continuous. If there is a positive scalar *λ* and a uniformly positive matrix Θ, such that

$$(\mathbf{F} + \mathbf{F}^T)/2 \le -\lambda I\_n \tag{A2}$$

then (A1) is said to be contracting, where *I<sup>n</sup>* denotes the identity matrix with dimension *n*, *F* = (Θ˙ + Θ*<sup>T</sup> <sup>∂</sup> <sup>f</sup>*(*x*,*t*) *<sup>∂</sup><sup>x</sup>* )Θ−<sup>1</sup> is the generalized Jacobian. For a contracting system, the trajectories starting from any initial condition will converge together exponentially. If *λ* = 0, (A1) is called semi-contracting and all its trajectories converge together asymptotically [2].

Partial contraction is a very important concept in CT [40]. Let the auxiliary system, called virtual system

$$
\dot{\varsigma} = f(\varsigma, \varkappa, t) \tag{A3}
$$

associated with (A1) through *f*(*x*, *x*, *t*) = *f*(*x*, *t*). Assume that (A3) is contracting with respect to *ς*, i.e., the Jacobian *<sup>∂</sup> <sup>f</sup>*(*ς*,*x*,*t*) *∂ζ* is UND for any *ς* and *x*.

If a particular solution of the virtual system verifies a smooth specific property, then all trajectories of the original *x*-system verify the same property exponentially. The original system is called partial contracting.

When a contracting system is subject to bounded disturbance, the error between the trajectory of the system after the disturbance and the original system trajectory is very small, that is, the contracting system is robust.

Consider the perturbed system:

$$\dot{\mathbf{x}}\_p = f(\mathbf{x}\_p, t) + d(\mathbf{x}\_p, t) \tag{A4}$$

where *d*(*xp*, *t*) is bounded, i.e., ∃ *d*<sup>0</sup> ≥ 0, ∀*xp*, ∀*t* ≥ 0, *d*(*xp*, *t*) ≤ *d*0. Then the error between *x*(*t*) and *xp*(*t*) satisfies

$$\|\mathbf{x}\_{\mathcal{P}}(t) - \mathbf{x}(t)\| \le \chi\_{\mathcal{P}} \|\mathbf{x}\_{\mathcal{P}}(0) - \mathbf{x}(0)\| \varepsilon \mathbf{x}^{-\lambda t} + \frac{\chi\_{\mathcal{P}} d\_0}{\lambda} \tag{A5}$$

where *χ<sup>p</sup>* is the upper bound of the condition number of Θ.

*Appendix B.2. Contraction Analysis of Singular Perturbation System*

For a standard singular perturbation system (SPS) [51]:

$$\begin{cases} \text{ } \pounds &= f(\mathfrak{x}, \mathfrak{z}) \\ \text{ } \epsilon \underline{z} &= \mathfrak{g}(\mathfrak{x}, \mathfrak{z}, \mathfrak{e}) \end{cases} \tag{A6}$$

0 <  1 is the singular perturbation parameter. For any *z*, if the virtual system *y*˙*<sup>x</sup>* = *f*(*yx*, *z*) is contracting with respect to *yx*, then the system (A6) is called to be partially contracting with respect to *x*. Similarly, it is partially contracting in *z* when the virtual system  *y*˙*<sup>z</sup>* = *g*(*x*, *yz*, ) is contracting for any *x* [3,40].

If the system (A6) is partially contracting in *z*, there exist a unique, global mapping between *x*, *z* and  [3], i.e., the algebraic equation *g*(*x*, *z*, ) = 0 can be equivalently written as *z* = *h*(*x*, ), here *h*(*x*, ) is called slow manifold or the quasi-steady state of the *z*subsystem. According to SPT, the *x*-subsystem can be simplified by introducing the slow manifold into the *x*-subsystem:

$$\dot{\mathbf{x}}\_{rc} = f(\mathbf{x}\_{rc}, h(\mathbf{x}\_{rc}, \boldsymbol{\epsilon})) \tag{A7}$$

To analyse the convergence behavior between *z* and the slow manifold *h*(*x*, ), we define a error variable *y* = *z* − *h*(*x*, ), and the dynamic for *y* can be expressed as:

$$\frac{dy}{d\tau} = \mathcal{g}(\mathbf{x}, y + h(\mathbf{x}, \boldsymbol{\varepsilon}), \boldsymbol{\varepsilon}) - \boldsymbol{\varepsilon}\frac{h(\cdot)}{dt} \tag{A8}$$

where *τ* = *<sup>t</sup>*  is a new time scale. Then, the dynamic behavior of the whole SPT (A6) can be determined by analysing the behavior of the reduced order state variables *xre* and *z*.

#### **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.
