Autonomous Underwater Vehicle (AUV) Motion Design: Integrated Path Planning and Trajectory Tracking Based on Model Predictive Control (MPC)

Abstract

This paper attempts to develop a unified model predictive control (MPC) method for integrated path planning and trajectory tracking of autonomous underwater vehicles (AUVs). To deal with the computational burden of online path planning, an event-triggered model predictive control (EMPC) method is introduced by using the environmental change as a triggering mechanism. A collision hazard function utilizing the changing rate of hazard as a triggering threshold is proposed to guarantee safety. We further give an illustration of how to calculate this threshold. Then, a Lyapunov-based model predictive control (LMPC) framework is developed for the AUV to solve the trajectory tracking problem. Leveraging a nonlinear integral sliding mode control strategy, we construct the contraction constraint within the formulated LMPC framework, thereby theoretically ensuring closed-loop stability. We derive the necessary and sufficient conditions for recursive feasibility, which are subsequently used to prove the closed-loop stability of the system. In the simulations, the proposed path planning and tracking control are verified separately and integrated and combined with static and dynamic obstacles.

1. Introduction

Autonomous underwater vehicles (AUVs), which are small crafts capable of autonomously executing underwater missions, boasting excellent concealment and high safety, have been widely applied in fields such as ocean exploration, environmental monitoring, and underwater rescue [1]. Tasks such as localization [2], path planning [1], trajectory tracking [3], and detection [4,5] are essential for AUVs. Path planning and trajectory tracking are the indispensable components of AUV control systems to perform tasks effectively [6].

Path planning requires algorithms to possess the ability to predict the future and select points within the environment [6]. Therefore, compared to commonly used methods such as the A* algorithm [7,8], the dynamic window approach [9], and theartificial potential field method [10,11], model predictive control (MPC), with its receding horizon feature and predictive characteristics, exhibits excellent compatibility with path planning [12].

An advanced approach to MPC was formulated in [3], emphasizing the synergy between localized planning and path tracking for autonomous vehicles. This framework introduced an enhanced MPC method, leveraging refined particle swarm optimization techniques, to tackle both planning and tracking challenges in an integrated manner. The work reported in [13] devised a path planning strategy for intelligent vehicles to avoid collisions with dynamic pedestrians, incorporating an attention-guided long short-term memory network, a modified social force model, and a systematic exploration of MPC. This methodology specifically addresses the scenario of pedestrian–vehicle interactions at unsignaled intersections, emphasizing dynamic conflict resolution. The research in [14] focused on addressing the challenge of robust collision-free formation navigation for a set of multiple unmanned surface vehicles, considering these vehicles as underactuated nonlinear systems operating under both state and input constraints. The objective was to devise a solution for these constrained unmanned surface vehicles to navigate safely while maintaining their formation.

It can be seen that MPC has a good ability to deal with dynamic obstacles. Moreover, the multi-horizon prediction into the future enables MPC to generate more executable paths [15]. However, MPC encounters the issue of high computational cost, which is unfavorable for online solutions [15]. To address the issue of computational burden, an event-triggered mechanism is often combined with MPC to solve problems. In [16], an EMPC strategy incorporating an adaptive APF was devised for autonomous electric vehicles, enabling both obstacle evasion and precise trajectory tracking. In [17], an EMPC approach was formulated to accomplish trajectory following and obstacle circumvention for a wheeled mobile robot, taking into account input constraints and external disturbances.

The complete AUV motion planning generally consists of path planning and trajectory tracking [18]. Common methods used to handle trajectory tracking include PID control [19,20], sliding mode control [21,22], backstepping control [23,24], and so on. In contrast to the majority of conventional control approaches, the distinctive strength of MPC resides in its inherent ability to systematically incorporate and manage system constraints directly within the controller’s design phase, which renders MPC an exceptionally robust platform for tackling diverse control engineering challenges [25].

In [25], an innovative LMPC framework was introduced for AUVs, leveraging computational resources through online optimization to enhance trajectory tracking performance. By incorporating a nonlinear backstepping tracking control law, the study formulated contraction constraints within the LMPC problem, theoretically ensuring closed-loop stability and optimizing the AUV’s trajectory execution. The work reported in [26] formulated a robust MPC approach tailored for AUV trajectory tracking amidst bounded disturbances. This tube-based MPC approach includes two distinct optimization tasks: the first tackles a standard MPC problem for the nominal system, establishing a reference state trajectory, while the second strives to confine the perturbed system’s state within a tube encircling this reference trajectory, thereby ensuring the robust stability of the closed-loop systems. Research work in [27] designed a trajectory tracking scheme using the sliding mode control technique in order to be robust against bounded disturbances. In order to align the planned position, [27] made strict assumptions on the reference state, i.e., that position tracking error converged while velocity error converged.

In summary, the existing research on path planning and trajectory tracking for AUVs has the following limitations. For path planning, EMPC is often only used for local obstacle avoidance during trajectory tracking (e.g., [16,17]) rather than the entire path planning process, which we believe provides a limited reduction in computational burden during planning. Furthermore, path planning with safety as a consideration focuses excessively on the degree of the collision hazard itself (e.g., [28,29]) while neglecting the changing rate of the collision hazard, which is often more closely related to environmental changes. For trajectory tracking, the use of sliding mode control methods (e.g., [27]) requires constraints on the reference trajectory that needs position tracking error to converge while velocity error converges, and due to the inherent chattering phenomenon in sliding mode control, the actual tracking performance is not ideal. Moreover, the chattering phenomenon in sliding mode control has multifaceted impacts on the system, including degrading system stability, affecting dynamic performance, damaging system hardware, reducing control accuracy, and limiting the scope of system application [30]. Existing LMPC methods typically use backstepping or linear sliding mode control (e.g., [25]) as auxiliary control laws and are mostly applied to other models rather than AUVs (e.g., [31]).

To address these issues, for path planning, this paper proposes an EMPC-based path planning algorithm operating during the entire process while leveraging the environmental changes as a triggering mechanism to reduce the computational burden. Additionally, we innovatively incorporate a collision hazard evaluation function with a threshold based on the changing rate of collision hazard, enabling the AUV to make selections between reachability and safety. For trajectory tracking, to avoid the chattering phenomenon caused by sliding mode control while ensuring the stability of MPC, this paper constructs an LMPC method for trajectory tracking using integral sliding mode control as an auxiliary control law. This approach not only guarantees stability but also ensures tracking performance, and there is no need to make additional assumptions on the reference state.

The main contributions of the paper are as follows:

(1)

An EMPC leveraging environmental changes as a triggering mechanism is developed for path planning. Compared to [1,12], our approach, by accepting feasible solutions according to environmental change not only ensures reachability but also imposes a lower computational burden.

(2)

We innovatively introduce the changing rate of the collision hazard as the condition for safety penalty. Considering the possibility that the environment itself is at a high risk, traditional methods such as [28,29] that solely rely on hazard evaluation functions to determine collision avoidance actions may cause AUsV to constantly attempt to flee in such an environment. The proposed method takes safety into account based on environmental changes, which allows a larger feasible area while balancing safety and reachability.

(3)

The LMPC augmented with integral sliding mode control as an auxiliary control law has been employed for the trajectory tracking of AUV. Compared to [27], the proposed tracking approach not only ensures tracking performance but also eliminates the chattering phenomenon without making any assumptions on the reference path (position tracking error is converged while velocity error converged).

The remaining part of this paper is organized as follows: In Section 2, we introduce our path planning method. In Section 3, we first design the auxiliary control law and then incorporate it into our LMPC-based trajectory tracking method. Feasibility and stability are proven rigorously. In Section 4, we introduce our integrated path planning and trajectory tracking method. In Section 5, we explain our simulation results. In Section 6, we draw a conclusion of our work.

2. Path Planning Method Based on EMPC and Changing Rate of Collision Hazard

In this section, we will first explain the environment and then establish a kinematic model that describes the movement of AUV, and subsequently, we will present an approach to path planning utilizing EMPC as the foundation.

2.1. Environment Description

To better illustrate our path planning methods, we will first explain the working environment of AUVs in this part. Thus, the following two assumptions are needed.

Assumption 1.

The environmental information obtained from the sensor of AUV is comprehensive.

Remark 1.

Assumption 1 differs from the limit that “Global information is known”. Here, we require that all environmental information within the sensor’s range is known, rather than requiring that the entire global information is previously known.

Assumption 1 imposes a prerequisite for path planning, The credibility of sensor information typically serves as the foundation for path planning and even the motion design of AUVs, as exemplified in [1,32]. Therefore, to pursue more precise path planning, the formulation of Assumption 1 is reasonable. Meanwhile, to facilitate the description of obstacle constraints, we make the following assumption.

Assumption 2.

All obstacles can be completely enclosed by circular areas.

Assumption 2 requires that circular areas be able to enclose obstacles. In fact, existing methods often use polygons or circles to enclose irregular obstacles such as [33]. In the real situation, the enclosed areas of these rules will simplify the description of obstacle constraints and simplify the calculation of MPC [33]. Here, for the sake of convenience, we uniformly use circular areas to enclose obstacles.

Assumptions 1 and 2 limit the path planning problem to an environment where sensor detection is possible and obstacle information can be formulated, laying a foundation for us to subsequently develop a path planning method based on MPC.

2.2. AUV Kinematics

Given the AUV’s sway velocity being substantially lower than its surge velocity, the nonholonomic kinematic equations effectively model the vehicle’s motion dynamics at the velocity level with sufficient accuracy, as demonstrated in [1]. Thus, the AUV motion in the horizontal plane is

$$\begin{array}{cc}\hfill & \dot{x}=u_{1}cos\psi -{\nu }_{1}sin\psi \approx vcos\psi \hfill \\ \hfill & \dot{y}=u_{1}sin\psi +{\nu }_{1}cos\psi \approx vsin\psi \hfill \\ \hfill & \dot{\psi }=r_1=\omega \hfill \end{array}$$

(1)

where $x,y,\psi$ is the desired position and orientation; $u_1,{\nu }_1,r_1$ are surge, sway, and yaw velocities; v is the resultant tangential velocity of the vehicle; and $\omega$ is the angular velocity. This can be written as

$$\dot{\mathit{x}}=\mathit{h}(\mathit{x},\mathit{u})$$

(2)

where $\mathit{x}={\left[\begin{array}{c}x,y,\psi \end{array}\right]}^{\mathrm{T}}$ is the state of the system, $\mathit{u}={\left[\begin{array}{c}v,\omega \end{array}\right]}^{\mathrm{T}}$ is the input.

MPC can be used for real-time path planning, but due to its computationally intensive nature, it is difficult to ensure the efficiency of online planning when using MPC. Additionally, existing path planning methods often prioritize reachability as the primary goal, but in hazardous situations, safety should be the primary consideration.

Therefore, the main objective in path planning of this paper is as follows: Plan a collision-free path for AUV (2) while ensuring computational efficiency, requiring the capability of trading off the priorities of reachability and safety.

2.3. Path Planning Method Based on EMPC

In this part, we aim to outline our proposed path-planning method. Crucial components of our method will be introduced next.

Problem 1.

Define $t_k$ as the event-triggered time. Let ${\mathit{x}}_{\mathit{tar}}$ = ${\left[\begin{array}{c}{x}_{tar},y_{tar},{\psi }_{tar}\end{array}\right]}^{\mathrm{T}}$ denote the target point; $g_i\left(t_k\right)$ denote the ground penalty which is the basic penalty of the environment; $A_s\left(t_k\right)$ denote the collision hazard evaluated at $t_k$, which is the risk assessment at $t_k$ [29], $M_{sum}$ denote the total number of obstacles; and $A_{\delta }$ denote the threshold of changing rate of collision hazard. For $\xi \in \left[\begin{array}{c}{t}_k,t_k+N_p\end{array}\right]$, the optimal problem can be formulated as

$$\begin{array}{cc}\hfill & \underset{\mathit{u},\mathit{x}}{min}{J}_p(\mathit{x}\left(t_k\right),\mathit{u}\left(t_k\right))\hfill \\ \hfill s.t.& \dot{\mathit{x}}\left(\xi \right|t_k)=\mathit{h}(\mathit{x}\left(\xi \right|t_k),\mathit{u}\left(\xi \right|t_k))\hfill \\ \hfill & \mathit{x}\left(\xi \right|t_k)\in {\mathit{X}}_{\mathit{p}}\hfill \\ \hfill & \mathit{u}\left(\xi \right|t_k)\in {\mathit{U}}_{\mathit{p}}\hfill \\ \hfill & \left|\right|x\left(\xi \right|t_k)-{\underline{x}}^i{\left|\right|}^2+\left|\right|y\left(\xi \right|t_k)-{\underline{y}}^i{\left|\right|}^2\ge \left|\right|{\underline{r}}^i{\left|\right|}^2\hfill \\ \hfill & (i\in [1,M_{sum}])\hfill \end{array}$$

where the cost function is

$$\begin{array}{cc}\hfill & J_p(\mathit{x}\left(t_k\right),\mathit{u}\left(t_k\right))=J_a(\mathit{x}\left(t_k\right),\mathit{u}\left(t_k\right)+{\int }_{t_k}^{t_k+N_p}\kappa {∥A_{\delta }\left(\xi \right|t_k)∥}_S)d\xi ,\hfill \\ \hfill & J_a(\mathit{x}\left(t_k\right),\mathit{u}\left(t_k\right))={\int }_{t_k}^{t_k+N_p}(||\mathit{x}\left(\xi \right|t_k)-{\mathit{x}}_{\mathit{tar}}{\left|\right|}_Q^2+||\mathit{u}\left(\xi \right|t_k){\left|\right|}_R^2+\left|\right|\sum _{j=1}^{M_{sum}}{g}_j\left(\xi \right|t_k){\left|\right|}_P)d\xi ,\hfill \end{array}$$

$N_p$ is the prediction horizon; Q, R, P and S are weighting matrices; ${\mathit{X}}_{\mathit{p}}$ and ${\mathit{U}}_{\mathit{p}}$ are, respectively, the region of state and input of AUV system (2); $({\underline{x}}^i,{\underline{y}}^i)$ is the position; and ${\underline{r}}^i$ is the radius of the ith obstacle. κ is a constant relevant to $A_{\delta }$, let $D_{A_s}\left(t\right)=\Delta A_s\left(t\right)-\Delta A_s(t-1)$,

$$\kappa =\left\{\begin{array}{cc}1\hfill & \Delta D_{A_s}\ge 0,A_s\ge A_{\delta }or{D}_{A_s}<0,A_s\ge \beta \hfill \\ 0\hfill & \Delta D_{A_s}\ge 0,A_s<A_{\delta }or{D}_{A_s}<0,A_s<\beta \hfill \end{array}\right.$$

(3)

where $\Delta A_s\left(t\right)={\displaystyle \frac{1}{T}}(A_s\left(t\right)-A_s(t-1))$, β is a positive constant to ensure that κ does not switch continuously around $A_{\delta }$.

The weights Q, R, P, and S are assigned based on the current objectives that need to be achieved, with a higher value assigned to the weight corresponding to the primary goal. Typically, to ensure reachability, Q is set to a relatively large value. Furthermore, to prioritize safety in critical situations, S, as the weight for the collision hazard function, needs to be set to a significantly high value. The ground penalty $g\left(t\right)$ is used to prevent the AUV from driving closely along the boundary of obstacles. $A_s\left(t\right)$, on the other hand, is a penalty based on safety, ensuring that the AUV prioritizes safety before reachability when it is in a dangerous situation. The setting of $A_{\delta }$ provides a boundary for switching between reachability and safety.

2.4. Event-Triggered Mechanism

Due to the significant computational load of MPC, it is challenging to ensure the efficiency of online planning. Therefore, we decide to introduce an event-triggered mechanism to enhance computational efficiency.

Since MPC can generate a sequence of optimal moves in a single computation, we propose to only re-plan when the environment changes (i.e., re-compute MPC). If the environment remains unchanged, we will utilize the previously computed feasible solution. Thus, we first need to define what a change in the environment is.

2.4.1. Definition of Environment Change

We categorize the following main cases as belonging to the realm of environmental changes.

(1)

New obstacles are added or existing obstacles are reduced.

(2)

The obstacle changes from stationary to dynamic.

(3)

The velocity or orientation of a dynamic obstacle changes (i.e., it does not match the previously predicted motion state).

Remark 2.

As a supplementary note to the third point, the presence of obstacles with constant speed and heading in the environment is not considered a change in the environment, as their subsequent positions and impact on our AUV can be predicted.

In particular, there are some other special cases in which we need to recalculate the MPC.

(1)

A new penalty term is added to the cost function (for instance, $A_s\left(t\right)$ is added when $\Delta A_s\left(t\right)\ge A_{\delta }$).

(2)

All feasible solutions obtained from a single MPC computation have been exhausted.

At these special moments, a trigger is also required.

Moreover, since it is difficult to mathematically explain how to avoid the Zeno phenomenon based on environmental changes as a triggering mechanism, we make the following assumption while considering practical situations.

Assumption 3.

There exists a minimal time $\varsigma >0$ during two changing actions of the environment.

Remark 3.

From a practical perspective, this assumption is reasonable. Firstly, the finite obstacles present in the environment are unlikely to be constantly altering their navigation postures. Additionally, in real-world scenarios, there exists a sampling time, where $\varsigma >0$ is equal to this sampling time, such as [16,17].

2.4.2. Determination of Environmental Change

Based on the conclusions of [28,29], the collision hazard evaluation function typically takes into account factors such as the relative speed, angle, and distance between ships. Therefore, the collision hazard evaluation function can effectively describe environmental information in a mathematical manner.

According to [28,29], the collision hazard evaluation function can be formulated as

$$A_s\left(t\right)=\sum _{j=1}^{M_{sum}}{K}_t{\left({\displaystyle \frac{d_{min}}{d_j\left(t\right)}}\right)}^q{∥v_R\left(t\right)∥}^2$$

(4)

where $d_{min}$ represents the minimum safe distance, $d_j\left(t\right)$ represents the distance between AUV and the jth obstacle, $v_R\left(t\right)$ is the relative speed, $K_t$ is a positive constant, and $q\ge 1$.

Therefore, we can predict the value of $A_s\left(t\right)$ at a fixed point regardless of our AUV to determine if the environment has changed, as shown in Figure 1. A is a fixed point; $O_1$, $O_2$, $O_3$ are stationary obstacles; and $T{S}_1$, $T{S}_2$ are dynamic obstacles whose velocities are $v_{TS1}$ and $v_{TS2}$, respectively.

At time t, we can predict $A_s\left(t\right)$ in the future N steps of fixed point A, written as $\left\{\begin{array}{c}{A}_{s1}^A,A_{s2}^A,A_{s3}^A,\dots ,A_{sN}^A\end{array}\right\}$, which is relevant to $\left\{\begin{array}{c}t,t+1,t+2,\dots ,t+N-1\end{array}\right\}$.

Compare the predicted collision hazard with the actual calculated collision hazard. If they are not equal, it indicates that the environment has changed, and the MPC needs to be recalculated. Therefore, the event-triggered mechanism based on environment changes can be given by

$$t_{k+1}=min\left\{\begin{array}{c}{\tilde{t}}_{k+1},t_k+N_p\end{array}\right\}$$

(5)

where

$${\tilde{t}}_{k+1}=\left\{\begin{array}{c}\xi :A_s\left(t\right)\ne A_{s.}^A\end{array}\right\}$$

The upper and lower bounds for the inter-event time ${\Delta }_{k+1}$ are specified as ${\Delta }_{max}=N_p$${\Delta }_{min}=\varsigma$, respectively. Thus, Zeno behavior can be avoided because of the designed event-triggered mechanism (5).

2.5. Design of Changing Rate-Based Collision Hazard Function

The existing literature, such as [28,29], often advocates for determining collision avoidance based on whether a ship is in a high-risk situation (i.e., checking if $A_s\left(t\right)$ exceeds a certain collision hazard threshold). However, ships may inherently operate in high-risk environments, resulting in $A_s\left(t\right)$ continuously being treated as a penalty in the cost function. This defeats the original purpose of setting a threshold. The threshold is intended to enable the system to switch between reachability and safety. If $A_s\left(t\right)$ is constantly included, the AUV will first escape to a safe area and then re-ensure reachability in a safer environment, significantly increasing the path length and consequently energy consumption. Moreover, paths in high-risk environments are not necessarily unusable. If the AUV always prioritizes safety, it may become overly conservative, losing its ability to operate in complex scenarios.

Therefore, this paper proposes using the rate of change in risk $A_{\delta }$ as a condition to determine whether to include $A_s$ as a penalty, which is analogous to how an airbag in a car deploys based on acceleration.

Calculation of $A_{\delta }$

As shown in Figure 2, if we consider the most extreme scenario, then the AUV and the dynamic obstacle are at the minimum safe distance from each other, and they are moving towards each other. The following equation can be satisfied.

$$\left|\right|v_{sum}{\left|\right|}_{max}=\left|\right|v_{OS}{\left|\right|}_{max}+\left|\right|v_{TS}{\left|\right|}_{max}$$

(6)

$$d_j=d_{min\_TS}+d_{min\_OS}=d_{min}$$

(7)

where $v_{sum}$ is the resultant velocity.

At this time,

$$A_s^{lim}=K_t{\left({\displaystyle \frac{d_{min}}{d_{min\_TS}+d_{min\_OS}}}\right)}^q{\left(v_{sum}\right)}^2=K_t{\left(v_{sum}\right)}^2$$

(8)

where $A_s^{lim}$ is the extreme collision hazard.

The extreme changing rate of the collision hazard is when it changes from 0 to $A_s^{lim}$ within a sampling step. And if we consider all the obstacles, it can be formulated as

$$A_{\delta }={\displaystyle \frac{1}{T}}\sum _{j=0}^{M_{sum}}{K}_t{\left(v_{sum}^j\right)}^2$$

(9)

where $v_{sum}^j$ is the resultant velocity with the jth obstacle.

Specifically, when the obstacle is stationary, the resultant velocity is the velocity of our AUV.

2.6. Algorithm of the Proposed Path Planning Method

The path planning algorithm is summarized as Algorithm 1.

Algorithm 1 Path Planning Method based on EMPC and changing rate of collision hazard

Input:  v (resultant tangential velocity), $\omega$ (yaw velocity) Output:  $\mathit{x}$ (position)  1: Given the initial state ${\mathit{x}}_{\mathbf{0}}$ of AUV and total duration $t_{total}$  2: Calculate the changing rate threshold $A_{\delta }$ by (8)  3: Set a fixed point A to judge if the environment changes  4: $t=0$  5: while $t\le t_{total}$ do  6:     Evaluate the collision hazard of environment $A_s\left(t\right)$ by (4)  7:     Calculate the changing rate of collision hazard $\Delta A_s\left(t\right)$  8:     Predict the collision hazard of A in the future steps, written as $\left\{\begin{array}{c}{A}_{s1}^A,A_{s2}^A,A_{s3}^A,\dots ,A_{sN}^A\end{array}\right\}$  9:     if $\Delta A_s\left(t\right)\ge A_{\delta }$ then 10:         Use $J_p$ as cost function 11:     else 12:         Use $J_a$ as cost function 13:     Solve Problem 1 to obtain the input sequence $\left\{\begin{array}{c}{u}_0,u_1,u_2,\dots ,u_{N-1}\end{array}\right\}$ 14:     Use $u_0$ as input at time t 15:     $t=t+1$, $k=1$ 16:     for $k\le N-1$ do 17:         if $A_s(t+k+1)==A_{sk}^A$ then 18:            Use $u_k$ as input at time $t+k$ 19:            $k=k+1$ 20:         else 21:            $t=t+k-1$ 22:            break 23: Record the position $\mathit{x}$ and input $\mathit{u}$

Remark 4.

In contrast to [29,32,34], which solely emphasizes the magnitude of collision hazard, the path planning method presented in this paper innovatively incorporates the changing rate of collision hazard. This additional consideration enables a more comprehensive assessment of potential hazards and allows a larger feasible area for planning since our method can work in a high collision hazard environment but [29,32,34] attempts to flee from such environment. Furthermore, by employing EMPC to balance between feasible and optimal solutions, our approach enhances planning efficiency. Consequently, the paths generated by our method are not only safer but also more rational, as they account for both the current risk level and its trend over time.

3. Trajectory Tracking Method Based on LMPC

In this section, we will introduce a trajectory-tracking method based on LMPC. Integral sliding mode control is used as an auxiliary controller.

3.1. AUV Dynamics

Based on the study of underwater vehicles in [35], we can model the AUV as follows:

$$\begin{array}{cc}\hfill & \dot{\mathit{\eta }}=R\left(\psi \right)\mathit{v}\hfill \\ \hfill & M\dot{\mathit{v}}+C\left(\mathit{v}\right)\mathit{v}+D\left(\mathit{v}\right)\mathit{v}+g\left(\mathit{\eta }\right)=\mathit{\tau }\hfill \end{array}$$

(10)

where $\mathit{\eta }={\left[\begin{array}{c}x,y,\psi \end{array}\right]}^T$ is the position and heading of the AUV; $\mathit{v}={\left[\begin{array}{c}u,\nu ,r\end{array}\right]}^T$ is the velocity of the vehicle; and $R\left(\psi \right)$ is the rotation matrix depending on the heading $\psi$.

$$R\left(\psi \right)=\left[\begin{array}{ccc}cos\left(\psi \right)& -sin\left(\psi \right)& 0\\ sin\left(\psi \right)& cos\left(\psi \right)& 0\\ 0& 0& 1\end{array}\right]$$

(11)

M is the inertia matrix including added mass, $C\left(\mathit{v}\right)$ is the state-dependent matrix of Coriolis and centripetal terms, $D\left(\mathit{v}\right)$ is the hydrodynamic damping and lift matrix, $g\left(\mathit{\eta }\right)$ is the vector of gravitational forces and moments, and $\mathit{\tau }={\left[\begin{array}{c}{\tau }_u,{\tau }_v,{\tau }_r\end{array}\right]}^{\mathrm{T}}$ is the vector of thrust forces and moments and the generalized control input is expressed as $\mathit{\tau }=\mathbf{B}\mathbf{u}$, where $\mathbf{B}$ is the control matrix and the control input is defined as $\mathit{u}={\left[\begin{array}{c}{u}_1,u_2,u_3,u_4\end{array}\right]}^{\mathrm{T}}$.

Additionally, by selecting the center of gravity as the origin of the body-fixed reference frame, the system matrices can be simplified as $M=\mathrm{diag}$$(M_{\dot{u}},M_{\dot{v}},M_{\dot{r}})$ and $D\left(\mathit{v}\right)=\mathrm{diag}$$(X_u,Y_v,N_r)+\mathrm{diag}$$(D_u\left|u\right|,D_v\left|\nu \right|,D_r\left|r\right|)$, $g\left(\mathit{\eta }\right)=0$ and

$$C\left(\mathit{v}\right)=\left[\begin{array}{ccc}0& 0& -M_{\dot{v}}\\ 0& 0& M_{\dot{u}}\\ M_{\dot{v}}\nu & -M_{\dot{u}}& 0\end{array}\right]$$

(12)

$M_x=m-X_{\dot{u}}$, $M_y=m-Y_{\dot{v}}$, and $M_{\psi }=I_z-N_{\dot{r}}$ are inertia terms including added mass; $X_u$, $Y_v$ and $N_r$ are linear drag terms; and $D_u$, $D_v$, and $D_r$ are corresponding quadratic drag terms. The overall control system model is built as follows:

$$\dot{\mathit{x}}=\left[\begin{array}{c}R\left(\psi \right)\mathit{v}\\ M^{-1}(\mathit{\tau }-C\left(\mathit{v}\right)\mathit{v}-D\left(\mathit{v}\right)\mathit{v})\end{array}\right]=\mathbf{f}(\mathit{x},\mathit{\tau }).$$

(13)

where the state is defined as $\mathit{x}={\left[\begin{array}{c}x,y,\psi ,u,\nu ,r\end{array}\right]}^{\mathrm{T}}$.

The main objective of trajectory tracking can be summarized as follows: Using the planned path ${\mathit{x}}^{\mathit{D}}$ from Algorithm 1 as a reference, control AUV (13) converges to a planned path ${\mathit{x}}^{\mathit{D}}$ by LMPC while using integral sliding mode control as an auxiliary control law.

3.2. Trajectory Tracking Method Based on LMPC

Reference states and reference inputs can be constructed based on the information obtained through path planning. The reference state can be obtained as follows.

$${\mathit{\eta }}^{\mathit{D}}=\mathit{x}=\left[\begin{array}{c}{x}^D,y^D,{\psi }^D\end{array}\right],\mathit{u}={\mathit{v}}^{\mathit{D}},u^D=v,{\nu }^D=0,r^D=\omega$$

(14)

where ${\mathit{\eta }}^{\mathit{D}}$ is the reference position and ${\mathit{v}}^{\mathit{D}}$ is the reference velocity. Therefore, the reference state ${\mathit{x}}^{\mathit{D}}$ is

$${\mathit{x}}^{\mathit{D}}=\left[\begin{array}{c}{x}^D,y^D,{\psi }^D,v,0,\omega \end{array}\right]$$

(15)

Moreover, the reference input ${\mathit{\tau }}^{\mathit{D}}$ can be obtained from (10). Thus, the problem of trajectory tracking based on LMPC can be formulated as follows.

Problem 2.

Let ${\mathit{\tau }}^{\mathit{h}}(.)$ denote the auxiliary Lyapunov-based nonlinear tracking control law and let $V(.)$ be the corresponding Lyapunov function; the optimal problem can be formulated as

$$\begin{array}{cc}\hfill & \underset{\mathit{\tau },\mathit{x}}{min}{J}_t(\mathit{x}\left(t\right),\mathit{\tau }\left(t\right))={\int }_0^{N-1}\ell (\mathit{x}\left(s\right),\mathit{\tau }\left(s\right))ds\hfill \end{array}$$

(16a)

$$\begin{array}{cc}\hfill s.t.& \dot{\mathit{x}}\left(s\right)=\mathbf{f}(\mathit{x}\left(s\right),\mathit{\tau }\left(s\right))\hfill \end{array}$$

(16b)

$$\begin{array}{cc}\hfill & \left|\mathit{\tau }\left(s\right)\right|\le {\mathit{\tau }}_{\mathbf{max}}\hfill \end{array}$$

(16c)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial V}{\partial \mathit{x}}}\mathbf{f}(\mathit{x}\left(0\right),\mathit{\tau }\left(0\right))\le {\displaystyle \frac{\partial V}{\partial \mathit{x}}}\mathbf{f}(\mathit{x}\left(0\right),{\mathit{\tau }}^{\mathit{h}}\left(0\right))\hfill \end{array}$$

(16d)

where $\ell (\mathit{x},\mathit{\tau })={\left(\mathit{x}-{\mathit{x}}^{\mathit{D}}\right)}^{\mathrm{T}}{Q}_t\left(\mathit{x}-{\mathit{x}}^{\mathit{D}}\right)+{(\mathit{\tau }-{\mathit{\tau }}^{\mathit{D}})}^{\mathrm{T}}{R}_t(\mathit{\tau }-{\mathit{\tau }}^{\mathit{D}})$, N is the prediction horizon, $Q_t$ and $R_t$ are weighting matrices.

The introduction of the contraction constraint (16d) facilitates our demonstration that the LMPC controller inherits the inherent stability traits of ${\mathit{\tau }}^{\mathit{h}}$, enabling a precise outline of a guaranteed region of attraction. Additionally, owing to the real-time optimization process, the LMPC controller dynamically executes optimal tracking control, adhering to the system’s physical constraints, thereby optimizing performance while ensuring system stability [25].

3.3. Design of Auxiliary Tracking Controller

In this part, we will construct the auxiliary tracking controller using integral sliding mode control.

Firstly, we define the following tracking position errors:

$$\begin{array}{cc}\hfill & e_x=x-x^D\hfill \\ \hfill & e_y=y-y^D\hfill \\ \hfill & e_{\psi }=\psi -{\psi }^D\hfill \end{array}$$

(17)

Taking the derivative of (17) and substituting in (13), we have

$$\begin{array}{cc}\hfill & {\dot{e}}_x=\dot{x}-{\dot{x}}^D=ucos\left(\psi \right)-\nu sin\left(\psi \right)-{\dot{x}}^D\hfill \\ \hfill & {\dot{e}}_y=\dot{y}-{\dot{y}}^D=usin\left(\psi \right)+\nu cos\left(\psi \right)-{\dot{y}}^D\hfill \\ \hfill & {\dot{e}}_{\psi }=\dot{\psi }-{\dot{\psi }}^D=r-r^D\hfill \end{array}$$

(18)

Taking the derivative of (18), we have

$$\begin{array}{cc}\hfill & {\ddot{e}}_x=\dot{u}cos\left(\psi \right)-ursin\left(\psi \right)-\dot{\nu }sin\left(\psi \right)-\nu rcos\left(\psi \right)-{\ddot{x}}^D\hfill \\ \hfill & {\ddot{e}}_y=\dot{u}sin\left(\psi \right)+urcos\left(\psi \right)+\dot{\nu }cos\left(\psi \right)-\nu rsin\left(\psi \right)-{\ddot{y}}^D\hfill \\ \hfill & {\ddot{e}}_{\psi }=\dot{r}-{\dot{r}}^D\hfill \end{array}$$

(19)

To devise the trajectory tracking auxiliary integral sliding mode controller, we shall select sliding surfaces with the objective of

$$\mathit{s}=\mathit{c}\mathit{e}+\dot{\mathit{e}}+\mathit{\lambda }{\int }_0^t\mathit{e}\left(\mathit{\tau }\right)d\mathit{\tau }$$

(20)

where $\mathit{s}={\left[\begin{array}{c}{s}_1,s_2,s_3\end{array}\right]}^{\mathrm{T}}$, $\mathit{c}={\left[\begin{array}{c}{c}_1,c_2,c_3\end{array}\right]}^{\mathrm{T}}$, $\mathit{e}={\left[\begin{array}{c}{e}_1,e_2,e_3\end{array}\right]}^{\mathrm{T}}$, $\mathit{\lambda }={\left[\begin{array}{c}{\lambda }_1,{\lambda }_2,{\lambda }_3\end{array}\right]}^{\mathrm{T}}$, $c_1$, $c_2$, $c_3$ and ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$ are positive constants. Taking the derivative of these sliding mode surfaces and substitute in (13), we have

$$\begin{array}{cc}\hfill {\dot{s}}_1& =c_1{\dot{e}}_x+{\ddot{e}}_x+{\lambda }_1{e}_x\hfill \\ \hfill & =c_1(ucos\left(\psi \right)-\nu sin\left(\psi \right)-{\dot{x}}^D)+{\lambda }_1(x-x^D)\hfill \\ \hfill & +({\displaystyle \frac{M_y}{M_x}}\nu r-{\displaystyle \frac{X_u}{M_x}}u-{\displaystyle \frac{D_u}{M_x}}u|u|+{\displaystyle \frac{{\tau }_u^h}{M_x}})cos\left(\psi \right)-ursin\left(\psi \right)\hfill \\ \hfill & -(-{\displaystyle \frac{M_x}{M_y}}ur-{\displaystyle \frac{Y_v}{M_y}}\nu -{\displaystyle \frac{D_v}{M_y}}\nu |\nu |+{\displaystyle \frac{{\tau }_v^h}{M_y}})sin\left(\psi \right)\hfill \\ \hfill & -\nu rcos\left(\psi \right)-{\ddot{x}}^D\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\dot{s}}_2& =c_2{\dot{e}}_y+{\ddot{e}}_y+{\lambda }_2{e}_y\hfill \\ \hfill & =c_1(usin\left(\psi \right)-\nu cos\left(\psi \right)-{\dot{y}}^D)+{\lambda }_2(y-y^D)\hfill \\ \hfill & +({\displaystyle \frac{M_y}{M_x}}\nu r-{\displaystyle \frac{X_u}{M_x}}u-{\displaystyle \frac{D_u}{M_x}}u|u|+{\displaystyle \frac{{\tau }_u^h}{M_x}})sin\left(\psi \right)+urcos\left(\psi \right)\hfill \\ \hfill & +(-{\displaystyle \frac{M_x}{M_y}}ur-{\displaystyle \frac{Y_v}{M_y}}\nu -{\displaystyle \frac{D_v}{M_y}}\nu |\nu |+{\displaystyle \frac{{\tau }_v^h}{M_y}})cos\left(\psi \right)\hfill \\ \hfill & -\nu rsin\left(\psi \right)-{\ddot{y}}^D\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\dot{s}}_3& =c_3{\dot{e}}_{\psi }+{\ddot{e}}_{\psi }+{\lambda }_3{e}_{\psi }\hfill \\ \hfill & =c_3(r-r^D)+{\displaystyle \frac{M_x-M_y}{M_{\psi }}}u\nu -{\displaystyle \frac{N_r}{M_{\psi }}}r-{\displaystyle \frac{D_r}{M_{\psi }}}r\left|r\right|+{\displaystyle \frac{{\tau }_r^h}{M_{\psi }}}\hfill \end{array}$$

(23)

To guarantee the convergence of these sliding surfaces to zero within a finite time frame, the subsequent dynamics will be applied to the sliding surfaces.

$$\begin{array}{cc}\hfill & {\dot{s}}_1=-k_1{s}_1-W_{1}sign\left(s_1\right)\hfill \\ \hfill & {\dot{s}}_2=-k_2{s}_2-W_{2}sign\left(s_2\right)\hfill \\ \hfill & {\dot{s}}_3=-k_3{s}_3-W_{3}sign\left(s_3\right)\hfill \end{array}$$

(24)

where $k_1$, $k_2$, $k_3$ and $W_1$, $W_2$, $W_3$ are positive constants.

Let

$$\begin{array}{cc}\hfill A_{s_1}& =-k_1{s}_1-W_{1}sign\left(s_1\right)+{\ddot{x}}^D+\nu rcos\left(\psi \right)+ursin\left(\psi \right)\hfill \\ \hfill & -c_1(ucos\left(\psi \right)-\nu sin\left(\psi \right)-{\dot{x}}^D)-{\lambda }_1(x-x^D)\hfill \\ \hfill & -{\displaystyle \frac{M_y}{M_x}}\nu rcos\left(\psi \right)+{\displaystyle \frac{X_u}{M_x}}ucos\left(\psi \right)+{\displaystyle \frac{D_u}{M_x}}u\left|u\right|cos\left(\psi \right)\hfill \\ \hfill & -{\displaystyle \frac{M_x}{M_y}}ursin\left(\psi \right)-{\displaystyle \frac{Y_v}{M_y}}\nu sin\left(\psi \right)-{\displaystyle \frac{D_v}{M_y}}\nu \left|\nu \right|sin\left(\psi \right)\hfill \end{array}$$

(25)

and

$$\begin{array}{cc}\hfill B_{s_2}& =-k_2{s}_2-W_{2}sign\left(s_2\right)+{\ddot{y}}^D+\nu rsin\left(\psi \right)-urcos\left(\psi \right)\hfill \\ \hfill & -c_2(usin\left(\psi \right)+\nu cos\left(\psi \right)-{\dot{y}}^D)-{\lambda }_2(y-y^D)\hfill \\ \hfill & -{\displaystyle \frac{M_y}{M_x}}\nu rsin\left(\psi \right)+{\displaystyle \frac{X_u}{M_x}}\nu sin\left(\psi \right)+{\displaystyle \frac{D_u}{M_x}}u\left|u\right|sin\left(\psi \right)\hfill \\ \hfill & +{\displaystyle \frac{M_x}{M_y}}urcos\left(\psi \right)+{\displaystyle \frac{Y_v}{M_y}}\nu cos\left(\psi \right)+{\displaystyle \frac{D_v}{M_y}}\nu \left|\nu \right|cos\left(\psi \right)\hfill \end{array}$$

(26)

Thus, we have

$$\begin{array}{cc}\hfill & A_{s_1}={\displaystyle \frac{{\tau }_u^h}{M_x}}cos\left(\psi \right)-{\displaystyle \frac{{\tau }_v^h}{M_y}}sin\left(\psi \right)\hfill \\ \hfill & B_{s_2}={\displaystyle \frac{{\tau }_u^h}{M_x}}sin\left(\psi \right)+{\displaystyle \frac{{\tau }_v^h}{M_y}}cos\left(\psi \right)\hfill \end{array}$$

(27)

From (27), we can obtain

$$\begin{array}{cc}\hfill & {\tau }_u^h=M_x(cos\left(\psi \right)A_{s_1}+sin\left(\psi \right)B_{s_2})\hfill \\ \hfill & {\tau }_v^h=M_y(-sin\left(\psi \right)A_{s_1}+cos\left(\psi \right)B_{s_2})\hfill \end{array}$$

(28)

Similarly, we can obtain

$$\begin{array}{cc}\hfill {\tau }_r^h& =M_{\psi }{\dot{s}}_3-M_{\psi }{c}_3(r-r^D)\hfill \\ \hfill & -M_{\psi }{\lambda }_3(\psi -{\psi }^D)-(M_x-M_y)u\nu +N_{r}r+D_{r}r\left|r\right|\hfill \end{array}$$

(29)

And the Lyapunov function candidate can be formulated as

$$V={\displaystyle \frac{1}{2}}{s}_1^2+{\displaystyle \frac{1}{2}}{s}_2^2+{\displaystyle \frac{1}{2}}{s}_3^2$$

(30)

Taking the derivative of (30), we have

$$\begin{array}{cc}\hfill \dot{V}& =s_1{\dot{s}}_1+s_2{\dot{s}}_2+s_3{\dot{s}}_3\hfill \\ \hfill & =-k_1{s}_1^2-W_1|s_1|-k_2{s}_2^2-W_2|s_2|-k_3{s}_3^2-W_3\left|s_3\right|\hfill \\ \hfill & \le 0\hfill \end{array}$$

(31)

Based on the standard Lyapunov proof, the closed-loop system under conditions (28) and (29) is globally asymptotically stable.

3.4. Recursive Feasibility and Stability of LMPC

Now, we will first analyze the recursive feasibility of Problem 2 and then explain the stability.

3.4.1. Feasibility Analysis

Theorem 1.

Suppose Assumptions 1–3 hold. Suppose the reference state ${\mathit{x}}^{\mathit{D}}$ and reference input from Algorithm 1 are obtained before trajectory tracking. Let $({\overline{\tau }}_u^h,{\overline{\tau }}_v^h,{\overline{\tau }}_r^h)$ denote the upper bound of $({\tau }_u^h,{\tau }_v^h,{\tau }_r^h)$, respectively. For Problem 2, if ${\mathit{\tau }}_{\mathbf{max}}\ge max\left\{\begin{array}{c}{\overline{\tau }}_u^h,{\overline{\tau }}_v^h,{\overline{\tau }}_r^h\end{array}\right\}$ is satisfied, then ${\mathit{\tau }}^{\mathit{h}}\le {\mathit{\tau }}_{\mathbf{max}}$ can hold. LMPC (16) admits recursive feasibility.

Proof. 

We observe that, given the current system state $\mathit{x}\left(t\right)$, ${\mathit{\tau }}^{\mathit{h}}\left(t\right)$ is always feasible for the LMPC problem (16), if the condition ${\mathit{\tau }}^{\mathit{h}}\le {\mathit{\tau }}_{\mathbf{max}}$ is satisfied. To prove that ${\mathit{\tau }}^{\mathit{h}}\le {\mathit{\tau }}_{\mathbf{max}}$ holds, we will transform this problem into a specific one.

We observe the auxiliary control laws derived from (28) and (29). Since $M_x$, $M_y$, $M_{\psi }$, $X_u$, $Y_v$, $N_r$ and $D_u$, $D_v$, $D_r$ are fundamental parameters of the vessel, they are constants, and all the reference path parameters ${\mathit{\eta }}^{\mathit{D}}$, ${\dot{\mathit{\eta }}}^{\mathit{D}}$ and ${\ddot{\mathit{\eta }}}^{\mathit{D}}$ are also determined once the path planning is completed. Therefore, the parameter that can influence the inequality ${\mathit{\tau }}^{\mathit{h}}\le {\mathit{\tau }}_{\mathbf{max}}$ is the velocity $\mathit{v}$. Hence, to ensure that the auxiliary control laws satisfy ${\mathit{\tau }}^{\mathit{h}}\le {\mathit{\tau }}_{\mathbf{max}}$, we need to guarantee that $\mathit{v}$ has an upper bound.

From (10), we can obtain

$$\mathit{v}=R^T\left(\psi \right)\dot{\mathit{\eta }}$$

(32)

Since $\parallel {\mathbf{R}}^{\mathrm{T}}{\left(\psi \right)\parallel }_{\infty }=max\left\{\right|cos\psi |+|sin\psi |,1\}\le \sqrt{2}$, we further obtain

$${\parallel \mathit{v}\parallel }_{\infty }\le \sqrt{2}{\parallel \dot{\mathit{\eta }}\parallel }_{\infty }$$

(33)

In order to obtain $\parallel \dot{\mathit{\eta }}{\parallel }_{\infty }$, we will perform the following transformation.

$$\begin{array}{cc}\hfill \parallel \dot{\mathit{\eta }}{\parallel }_{\infty }& =\parallel \dot{\mathit{\eta }}-{\dot{\mathit{\eta }}}^{\mathit{D}}+{\dot{\mathit{\eta }}}^{\mathit{D}}{\parallel }_{\infty }\hfill \\ \hfill & \le \parallel \dot{\mathit{\eta }}-{\dot{\mathit{\eta }}}^{\mathit{D}}{\parallel }_{\infty }+{\parallel {\dot{\mathit{\eta }}}^{\mathit{D}}\parallel }_{\infty }\hfill \\ \hfill & =\parallel \dot{\mathit{e}}{\parallel }_{\infty }+{\parallel {\dot{\mathit{\eta }}}^{\mathit{D}}\parallel }_{\infty }\hfill \end{array}$$

(34)

${\dot{\mathit{\eta }}}^{\mathit{D}}$ is a certain number after path planning is finished; let

$$\parallel {\dot{\mathit{\eta }}}^{\mathit{D}}{\parallel }_{\infty }={\overline{\dot{\mathit{\eta }}}}^{\mathit{D}}=\left[\begin{array}{c}{\overline{\dot{x}}}^D,{\overline{\dot{y}}}^D,{\overline{\dot{\psi }}}^D\end{array}\right]$$

(35)

From (30), we have

$$\parallel \mathit{\eta }\left(t\right)-{\mathit{\eta }}^{\mathit{D}}{\left(t\right)\parallel }_2\le {\parallel \mathit{\eta }\left(0\right)-{\mathit{\eta }}^{\mathit{D}}\left(0\right)\parallel }_2$$

(36)

In other words, we can reformulate it as

$$\begin{array}{cc}\hfill {\parallel \mathit{e}\left(t\right)\parallel }_2& =\parallel \mathit{\eta }\left(t\right)-{\mathit{\eta }}^{\mathit{D}}{\left(t\right)\parallel }_2\hfill \\ \hfill & \le \parallel \mathit{\eta }\left(0\right)-{\mathit{\eta }}^{\mathit{D}}{\left(0\right)\parallel }_2\hfill \\ \hfill & ={\parallel \mathit{e}\left(0\right)\parallel }_2\hfill \end{array}$$

(37)

Further, according to (20), we have

$$\dot{\mathit{e}}=\mathit{s}-\mathit{c}\mathit{e}-\mathit{\lambda }{\int }_0^t\mathit{e}\left(\mathit{\tau }\right)d\mathit{\tau }$$

(38)

Therefore, we have

$$\begin{array}{cc}\hfill \parallel \dot{\mathit{e}}{\parallel }_{\infty }& =\parallel \mathit{s}-\mathit{c}\mathit{e}-\mathit{\lambda }{\int }_0^t{\mathit{e}\left(\mathit{\tau }\right)d\mathit{\tau }\parallel }_{\infty }\hfill \\ \hfill & \le {\parallel \mathit{s}\parallel }_{\infty }+{\parallel \mathit{ce}\parallel }_{\infty }+{\parallel \mathit{\lambda }{\int }_0^t\mathit{e}\left(\mathit{\tau }\right)d\mathit{\tau }\parallel }_{\infty }\hfill \end{array}$$

(39)

Due to the sign-preserving property of integration and the fact that ${\parallel \mathit{e}\left(t\right)\parallel }_{\infty }\le {\parallel \mathit{e}\left(t\right)\parallel }_2\le {\parallel \mathit{e}\left(0\right)\parallel }_2$, we have

$$\begin{array}{cc}\hfill \left|\right|\dot{\mathit{e}}{\left|\right|}_{\infty }& \le 0+c_{max}{\parallel \mathit{e}\left(0\right)\parallel }_2+t{\lambda }_{max}{\parallel \mathit{e}\left(0\right)\parallel }_2\hfill \\ \hfill & =(c_{max}+t{\lambda }_{max}){\parallel \mathit{e}\left(0\right)\parallel }_2\hfill \end{array}$$

(40)

where $c_{max}=max\left\{\begin{array}{c}{c}_1,c_2,c_3\end{array}\right\}$ and ${\lambda }_{max}=\left\{\begin{array}{c}{\lambda }_1,{\lambda }_2,{\lambda }_3\end{array}\right\}$.

Substituting in (34), we have

$$\parallel \dot{\mathit{\eta }}{\parallel }_{\infty }\le (c_{max}+t{\lambda }_{max}){\parallel \mathit{e}\left(0\right)\parallel }_2+{\overline{\dot{\mathit{\eta }}}}^{\mathit{D}}$$

(41)

Thus,

$${\parallel \mathit{v}\parallel }_{\infty }\le \overline{v}=\sqrt{2}(c_{max}+t{\lambda }_{max}){\parallel \mathit{e}\left(0\right)\parallel }_2+\sqrt{2}{\overline{\dot{\mathit{\eta }}}}^{\mathit{D}}$$

(42)

From the AUV model, we have the fact that the initial matrix is symmetric positive definite and upper bounded: $\infty >\overline{M}\mathit{I}\ge M=M^T>0$ [25].

Therefore, substituting in $C\left(\mathit{v}\right)$ and $D\left(\mathit{v}\right)$, we have

$${\left|\right|C\left(\mathit{v}\right)\left|\right|}_{\infty }\le \overline{C}=\sqrt{2}(c_{max}+t{\lambda }_{max}){\parallel \mathit{e}\left(0\right)\parallel }_2\overline{M}+\sqrt{2}{\overline{\dot{\mathit{\eta }}}}^{\mathit{D}}\overline{M}$$

(43)

$${\left|\right|D\left(\mathit{v}\right)\left|\right|}_{\infty }\le \overline{D}=\sqrt{2}{\overline{D}}_1+\sqrt{2}{\overline{D}}_2\left[\begin{array}{c}(c_{max}+t{\lambda }_{max}){\parallel \mathit{e}\left(0\right)\parallel }_2\overline{M}+{\overline{\dot{\mathit{\eta }}}}^{\mathit{D}}\overline{M}\end{array}\right]$$

(44)

where ${\overline{D}}_1=max\left\{\begin{array}{c}|X_u|,|Y_v|,|N_r|\end{array}\right\}$ and ${\overline{D}}_2=max\left\{\begin{array}{c}{D}_u,D_v,D_r\end{array}\right\}$.

Substituting (37), (41), (43), (44) in (25), (26), we have

$$\begin{array}{cc}\hfill A_{s_1}& \le {\overline{A}}_{s1}\hfill \\ \hfill & ={\ddot{x}}^D+{\overline{v}}^{2}cos\left(\psi \right)+{\overline{v}}^{2}sin\left(\psi \right)+c_1\left|\right|\dot{\mathit{\eta }}-{\dot{\mathit{\eta }}}^{\mathit{D}}{\left|\right|}_{\infty }+{\lambda }_1\left|\right|\mathit{\eta }-{\mathit{\eta }}^{\mathit{D}}{\left|\right|}_{\infty }+\overline{C}\overline{v}cos\left(\psi \right)+\overline{D}\overline{v}sin\left(\psi \right)\hfill \\ \hfill & ={\ddot{x}}^D+\sqrt{2}{\overline{v}}^2+c_1\left[\begin{array}{c}(c_{max}+t{\lambda }_{max}){\parallel \mathit{e}\left(0\right)\parallel }_2\end{array}\right]+{\lambda }_1{\parallel \mathit{e}\left(0\right)\parallel }_2+\overline{v}\sqrt{{\overline{C}}^2+{\overline{D}}^2}\hfill \end{array}$$

(45)

and

$$\begin{array}{cc}\hfill B_{s_2}& \le {\overline{B}}_{s_2}\hfill \\ \hfill & ={\ddot{y}}^D+{\overline{v}}^{2}sin\left(\psi \right)+{\overline{v}}^{2}cos\left(\psi \right)+c_2\left|\right|\dot{\mathit{\eta }}-{\dot{\mathit{\eta }}}^{\mathit{D}}{\left|\right|}_{\infty }+{\lambda }_2\left|\right|\mathit{\eta }-{\mathit{\eta }}^{\mathit{D}}{\left|\right|}_{\infty }+\overline{C}\overline{v}sin\left(\psi \right)+\overline{D}\overline{v}cos\left(\psi \right)\hfill \\ \hfill & ={\ddot{y}}^D+\sqrt{2}{\overline{v}}^2+c_2\left[\begin{array}{c}(c_{max}+t{\lambda }_{max}){\parallel \mathit{e}\left(0\right)\parallel }_2\end{array}\right]+{\lambda }_2{\parallel \mathit{e}\left(0\right)\parallel }_2+\overline{v}\sqrt{{\overline{C}}^2+{\overline{D}}^2}\hfill \end{array}$$

(46)

Therefore, substituting (45) and (46) in (28), we have

$$\begin{array}{cc}\hfill {\tau }_u^h& =\sqrt{A_{s_1}^2+B_{s_2}^2}{M}_x({\displaystyle \frac{A_{s_1}}{\sqrt{A_{s_1}^2+B_{s_2}^2}}}cos\left(\psi \right)+{\displaystyle \frac{B_{s_2}}{\sqrt{B_{s_2}^2+B_{s_2}^2}}}sin\left(\psi \right))\hfill \\ \hfill & \le \sqrt{{\overline{A}}_{s_1}^2+{\overline{B}}_{s_2}^2}{M}_x\hfill \\ \hfill & ={\overline{\tau }}_u^h\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill {\tau }_v^h& =\sqrt{A_{s_1}^2+B_{s_2}^2}{M}_y(-{\displaystyle \frac{A_{s_1}}{\sqrt{A_{s_1}^2+B_{s_2}^2}}}sin\left(\psi \right)+{\displaystyle \frac{B_{s_2}}{\sqrt{B_{s_2}^2+B_{s_2}^2}}}cos\left(\psi \right))\hfill \\ \hfill & \le \sqrt{{\overline{A}}_{s_1}^2+{\overline{B}}_{s_2}^2}{M}_y\hfill \\ \hfill & ={\overline{\tau }}_v^h\hfill \end{array}$$

(48)

Similarly, we have

$$\begin{array}{cc}\hfill {\tau }_r^h& \le M_{\psi }{c}_3\left|\right|\dot{\mathit{\eta }}-{\dot{\mathit{\eta }}}^{\mathit{D}}{\left|\right|}_{\infty }+M_{\psi }{\lambda }_3\left|\right|\mathit{\eta }-{\mathit{\eta }}^{\mathit{D}}{\left|\right|}_{\infty }+\overline{C}+\overline{D}\hfill \\ \hfill & \le M_{\psi }{c}_3\left[\begin{array}{c}(c_{max}+t{\lambda }_{max}){\parallel \mathit{e}\left(0\right)\parallel }_2\end{array}\right]+{\lambda }_3{\parallel \mathit{e}\left(0\right)\parallel }_2+\overline{C}+\overline{D}\hfill \\ \hfill & ={\overline{\tau }}_r^h\hfill \end{array}$$

(49)

Thus,

$$\begin{array}{cc}\hfill \left|\right|{\mathit{\tau }}^{\mathit{h}}{\left|\right|}_{\infty }& =max\left\{\begin{array}{c}{\overline{\tau }}_u^h,{\overline{\tau }}_v^h,{\overline{\tau }}_r^h\end{array}\right\}\hfill \end{array}$$

(50)

Thus, if ${\mathit{\tau }}_{\mathbf{max}}\ge max\left\{\begin{array}{c}{\overline{\tau }}_u^h,{\overline{\tau }}_v^h,{\overline{\tau }}_r^h\end{array}\right\}$ is satisfied, then ${\mathit{\tau }}^{\mathit{h}}\le {\mathit{\tau }}_{\mathbf{max}}$ holds. The proposed LMPC is recursively feasible.    □

It is noteworthy that the integral sliding mode controller merely serves as a preliminary estimation for the optimization problem (16) without being directly implemented on the system. Consequently, the auxiliary controller ensures feasibility, while the LMPC controller refines the control performance through optimization, achieving a more effective control strategy.

3.4.2. Asymptotic Stability Analysis

Theorem 2.

Suppose Assumption 1–3 holds. Suppose the reference state ${\mathit{x}}^{\mathit{D}}$ and reference input from Algorithm 1 are obtained before trajectory tracking. Then, the closed-loop system is asymptotically stable; i.e., the AUV will converge to the planned path ${\mathit{x}}^{\mathit{D}}$ with the control law $\mathit{\tau }$ based on LMPC.

Proof. 

Regarding the stability proof, given that all feasible solutions satisfy the constraint imposed by (16d) in Problem 2, and considering the fact that the constructed auxiliary controller satisfies $\dot{V}\le 0$ according to (31), we can consequently prove the stability by combining these two conditions. The formulaic representation of this proof can be stated as follows: According to converse Lyapunov theorems [36], there exists functions ${\alpha }_i(.)$, $i=1,2,3$ which belong to class $K_{\infty }$ such that the following inequalities hold:

$$\begin{array}{cc}\hfill & {\alpha }_1(\parallel \mathit{x}\parallel )\le V\left(\mathit{x}\right)\le {\alpha }_2(\parallel \mathit{x}\parallel )\hfill \\ \hfill & {\displaystyle \frac{\partial V}{\partial \mathit{x}}}\mathbf{f}(\mathit{x},{\mathit{\tau }}^{\mathit{h}})\le -{\alpha }_3(\parallel \mathit{x}\parallel )\hfill \end{array}$$

(51)

According to (16d), we have

$${\displaystyle \frac{\partial V}{\partial \mathit{x}}}\mathbf{f}(\mathit{x},\mathit{\tau })\le {\displaystyle \frac{\partial V}{\partial \mathit{x}}}\mathbf{f}(\mathit{x},{\mathit{\tau }}^{\mathit{h}})\le -{\alpha }_3(\parallel \mathit{x}\parallel )$$

(52)

Utilizing the Lyapunov argument, we claim that the closed-loop system is asymptotically stable, with a guaranteed region of attraction denoted by $\mathbb{R}$ (where $\mathbb{R}$ is the set of rational numbers).    □

3.5. Algorithm of Proposed Trajectory Tracking Method

Our algorithm can be summarized as Algorithm 2.

Algorithm 2 Trajectory Tracking Method based on LMPC

Input:  ${\mathit{x}}_{\mathbf{0}}$ (initial state) Output:  $\mathit{\tau }$ (current input), $\mathit{x}$ (next state) 1: Given an initial state ${\mathit{x}}_{\mathbf{0}}$ 2: $t=0$ 3: for  $t\le t_{max}$  do 4:     Solve Problem 2 to get the predictive sequence of $\mathit{\tau }$ 5:     Use the first item of the sequence of $\mathit{\tau }$ as the (sub-)optimal solution at the current time 6:     $t=t+1$

Remark 5.

In Algorithm 2, the adoption of suboptimal solutions is entirely acceptable and highly desirable for any nonlinear MPC algorithm, given the following pivotal facts. Firstly, owing to the nonlinearity of the system dynamics, the best-guaranteed solution to (16) obtained through iterative methods is inherently a local optimum. Crucially, when deploying on embedded platforms with limited computational prowess, restricting the iteration count becomes essential to sustain real-time control. This implies embracing suboptimal strategies, offering a crucial balance between numerical efficiency and control effectiveness. This enables a straightforward trade-off by specifying the maximum iteration count, ensuring that the stability of the tracking control remains intact while optimizing for computational efficiency.

Remark 6.

Compared to [27], the method proposed in this paper offers several key advantages. Firstly, it does not impose restrictions on the reference trajectory, such as requiring both position and velocity to converge simultaneously in [27]. Instead, it directly addresses position errors, thereby providing a broader range of applicability. Furthermore, by utilizing LMPC for trajectory tracking, our approach avoids the chattering phenomenon inherent in simple sliding mode control methods employed in [27]. Additionally, the convergence speed and tracking performance achieved by our LMPC-based approach surpass those of simply using sliding mode trajectory tracking, demonstrating its superiority in terms of both stability and accuracy.

4. Integrated Path Planning and Trajectory Tracking Method of AUV

The integrated EMPC-based path planning and LMPC-based trajectory tracking framework of AUV is illustrated in Figure 3. The planner calculates the environmental collision hazard $A_s\left(t\right)$ and its rate $\Delta A_s\left(t\right)$ of change based on the actual state and judges whether it is necessary to re-plan or introduce new penalties; the tracking controller tracks the reference path ${\mathit{x}}^{\mathit{D}}$ obtained from planning, using LMPC $\mathit{\tau }$ with the integral sliding mode control as an auxiliary control law ${\mathit{\tau }}^{\mathit{h}}$.

5. Simulation Results

5.1. Path Planning

First, we will validate path planning. We set the parameters as follows: the sampling time is $T=0.2$ s, the prediction horizon in path planning is $N_p=10$, and the weighting matrices are as follows: $Q=\mathrm{diag}(1,5,0.1)$, $R=0.01$, $P=1$, $S=7$, maximal velocity is $v_{max}=0.8$ m/s, minimal velocity is $v_{min}=-v_{max}$, maximal angular velocity is ${\omega }_{max}=\pi /4$ rad/s, and minimal angular velocity is ${\omega }_{min}=-{\omega }_{max}$. The changing rate threshold is set as $A_{\delta }=67$ (m/s)², where $K_t=3$.

As shown in Figure 4, Figure 4a represents the initial environment, where there are two dynamic obstacles. However, they are both in a static state for an initial period of time. When they start moving, their speeds are as follows: One of them moves solely horizontally with a speed of $0.025$ m/s, while the other moves both horizontally with a speed of $0.005$ m/s and vertically with a velocity of $0.01$ m/s. The red line segments indicate the predicted sequence of the path generated by the MPC path planner. In the period when the dynamic obstacle is initially stationary, MPC is not updated (i.e., only utilizing the feasible solution without recalculation), AUV will follow the entire predicted sequence of the path, instead of strictly utilizing the first item of the predicted sequence, as illustrated in Figure 4b. When the dynamic obstacle begins to move, MPC recalculates, and at this point, instead of utilizing a feasible solution, it employs the optimal solution as input, as depicted in Figure 4c, where the black hollow circles represent the initial position of dynamic obstacles. Figure 4d–f, illustrate the collision avoidance process of AUV with dynamic obstacles. Since the speed and direction of the obstacle remained unchanged, feasible solutions were used throughout the entire process. Figure 4g demonstrates that AUV has successfully completed the collision avoidance, and Figure 4h shows that AUV has finally reached its target point. Figure 5 shows the velocity and angular velocity changes during path planning.

Other situations are shown in Figure 6 and Figure 7. Meanwhile, in order to verify the improvement of computing efficiency by an event-triggered mechanism, we recorded the number of MPC calculations in the above three environments, as shown in

Table 1. Calculation times of the MPC.

	No Event-Triggered Mechanism	Event-Based
Figure 4	141	19
Figure 6	163	21
Figure 7	172	24

. It can be seen that under the action of an event-triggered mechanism, the calculation time of MPC decreases significantly.

5.2. Trajectory Tracking

In this section, we will validate the proposed trajectory tracking algorithm using the sine function and conduct a comparative analysis with the straightforward approach of employing integral sliding mode control. The parameter of LMPC settings are outlined as follows: $M_{\dot{u}}=283.6, M_{\dot{v}}=593.2, M_{\dot{r}}=29$, $X_u=26.9, Y_v=35.8, N_r=3.5$, $D_u=241.3, D_v=503.8, D_r=76.9$, the initial state is set as $\mathit{x}={(0,0,\pi /3,0,0,0)}^{\mathrm{T}}$, the prediction horizon is $N=10$, the sampling time is $T=0.2$ s, the total duration is set to $120\ast 0.2=24$ s, the weighting matrices are $Q_t=\mathrm{diag}({10}^4,{10}^4,{10}^3,10,10,10)$, $R_t={10}^{(-8)}\ast I$, and I is unit matrix. The parameters of integral sliding mode control are $c_1=6$, $c_2=6$, $c_3=6$, $k_1=3$, $k_2=2$, $k_3=2$, $W_1=0.1$, $W_2=0.1$, $W_3=0.1$, ${\lambda }_1=1$, ${\lambda }_2=3$, ${\lambda }_3=0.4$, which is the same as the auxiliary controller.

Figure 8 illustrates the simulation results of trajectory tracking using both integral sliding mode control and LMPC. In this figure, the blue line represents the reference trajectory, the green line denotes the tracking trajectory achieved by integral sliding mode control, and the red line indicates the tracking trajectory by LMPC. It can be seen that the chattering phenomenon occurs when the sliding mode control is used, which affects the tracking effect. However, LMPC outperforms the straightforward application of integral sliding mode control in terms of tracking performance. Moreover, LMPC exhibits a significantly faster convergence speed in the initial stage compared to integral sliding mode control, and it is capable of avoiding the chattering phenomenon that may occur in integral sliding mode control. Figure 9 depicts the state tracking performance, and it can be observed that the integral sliding mode control exhibits relatively large tracking errors in the initial stage, requiring a longer convergence time to the reference trajectory compared to LMPC. Furthermore, the chattering phenomenon significantly impacts the tracking performance across each dimension. Figure 10 shows the simulation results of the input signals.

5.3. Integrated Path Planning and Trajectory Tracking

In this section, we will conduct simulation validation for the proposed integrated path planning and trajectory tracking approach. All parameter selections are identical to those previously defined, and the environment setup for the path planning component remains unchanged from previous discussions. Figure 11 presents the simulation results, demonstrating that the proposed method can effectively accomplish the task of integrated path planning and trajectory tracking with satisfactory performance. Figure 12 illustrates the input results for this simulation.

6. Conclusions

This paper presented a comprehensive approach that integrates path planning and trajectory tracking for AUVs. During the path planning phase, we employed an event-triggered mechanism to reduce the computational load of path planning. Furthermore, a collision hazard function utilizing the changing rate of hazard as a threshold was proposed to ensure safety. We subsequently elaborated on how to calculate this threshold. In the trajectory tracking phase, we constructed an LMPC trajectory tracking method with integral sliding mode control serving as an auxiliary control law. This approach not only guaranteed stability and feasibility but also outperformed the direct application of integral sliding mode control in terms of performance. In the future, we will endeavor to optimize the computational complexity of LMPC and make fuller use of the prediction horizon in MPC for path planning.