**1. Introduction**

Autonomous underwater vehicles (AUVs) have assumed indispensable roles in various underwater operations, such as ocean exploration and hydrologic surveys [1]. They can autonomously perform appropriate maneuvers to achieve predefined objectives. Compared with the operational capability of a single AUV, collaborative AUVs can respond more reliably and flexibly to complex missions and extended operational ranges, thereby improving the efficiency and robustness of undersea operations. Given this backdrop, numerous application cases about AUV coordinated formation have been triggered in both civilian and industrial fields for decades [2,3]. Irrespective of the specific collaborative missions undertaken by AUVs, the core challenge lies in ensuring motion stability of AUV formations within complex underwater environments and the constraints of their own models. To tackle this problem, several mainstream methodologies have been proposed by engineers and academics. Studies by Chen et al. [4] and Zhen et al. [5] proposed AUV formation control schemes combined with the virtual structure method. However, this approach suffers from limited flexibility and applicability. Wang et al. [6] utilized the leader–follower method to address the AUV formation tracking problem, but this approach relies on the state of the leader, reducing the robustness and fault tolerance of the formation. Conversely, leaderless formations have been proposed promisingly and have received more considerable attention [7]. Munir et al. [8] proposed a new arbitrary-order distributed control strategy based on the novel sliding surfaces of error dynamics, which addresses the

**Citation:** Zhang, M.; Yan, Z.; Zhou, J.; Yue, L. Distributed Dual Closed-Loop Model Predictive Formation Control for Collision-Free Multi-AUV System Subject to Compound Disturbances. *J. Mar. Sci. Eng.* **2023**, *11*, 1897. https:// doi.org/10.3390/jmse11101897

Academic Editor: Sergei Chernyi

Received: 12 September 2023 Revised: 25 September 2023 Accepted: 28 September 2023 Published: 29 September 2023

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

cooperative tracking control of uncertain higher-order nonlinear systems. The strategies to mitigate the chattering issue caused by sliding surfaces are discussed in [9]. Despite the abundance of existing research, multiple-AUV formation tracking control remains a significantly challenging project.

One of the main challenges is the various disturbances resulting from the underwater environment and the motion model of the AUVs themselves [10]. On the one hand, unknown disturbances such as waves, tides, and currents, are inevitable in practical marine environments. On the other hand, AUVs exhibit highly nonlinear and coupled dynamics, leading to model uncertainties. These uncertainties are often induced by modeling errors and deviations in hydrodynamic coefficient measurement. According to the research by Cui et al. [11], these external disturbances and model uncertainties that degrade the system performance negatively are referred to as compound disturbances. In response to these challenges, researchers have developed diverse schemes, such as disturbance observers [12], fuzzy logic theory [13], and neural networks [14]. Among these, the extended state observer (ESO) initially proposed by Han [15] is an attractive option to estimate compound disturbances, as it does not rely on an accurate model. Lei et al. [16] designed a high-gain ESO to solve AUV horizontal trajectory tracking problems under the time-varying disturbances. Although many ESOs have been established for different platforms, most only guarantee asymptotic convergence of estimation errors, implying a potentially infinite convergence time. Some research works also lack a rigorous analysis of convergence. Considering the impact of severe underwater environments on estimation accuracy, the concept of finite-time ESO proves more beneficial for improving control performance [17]. Wang et al. [18] implemented a FTESO-based nonsingular terminal sliding mode controller to address unmanned surface vehicle (USV) trajectory tracking in disturbed environments. This approach ensured that the disturbance estimation errors converge within a finite time. However, there remains room for improvement and optimization of the design structure to further enhance observation performance.

AUV formation navigation also presents significant technical challenges due to various complex constraints. For instance, the AUV attitude has a certain desired range and navigation velocities are inherently limited. These intrinsic input and state constraints pose substantial challenges to control performance [19]. In practical applications, actuators often have input saturation constraints due to physical structure limitations. This results in a limitation of the actual active control force of the AUV. If a control signal exceeds this boundary, it may lead to system instability. However, most previous work assumes that the actuators can tolerate any level of control signals. To avoid actuator saturation, a nonlinear auxiliary system for filtering saturation errors was proposed [20]. Additionally, collisions between AUVs are undesirable during the formation configuration phase. Thus, the ability to avoid collisions is vital for AUV formation control. A wealth of solutions have been developed to this end, with Li and Wang [21] proposing a collision-free position consensus algorithm for AUVs based on potential function. Moreover, Xu et al. [22] presented an event-triggered algorithm based on deep reinforcement learning to avoid AUV collisions. However, the above studies disregard the physical constraints of AUVs. From the perspective of safe navigation, it is essential to integrate factors such as input, state restrictions, and collision avoidance into the design scheme.

Model predictive control (MPC) has garnered considerable attention due to its ability to simultaneously handle multiple composite constraints and offer superior dynamic performance. This is widely applied to MIMO systems affected by model distortions and complex constraints. Several MPC-based applications have been integrated into AUV control systems. Zhang et al. [23] proposed an MPC-based AUV trajectory tracking strategy under random disturbances. In [24], a robust model-predictive control scheme based on the active disturbance rejection control approach was developed for the AUV tracking task. The challenge of extending these systems to multi-AUV systems involves coordinating the control behavior of each subsystem and ensuring the closed-loop stability of the local MPC optimization problem under system constraints. This coordination aims to maximize the overall control performance. Hence, DMPC came into being. Zheng et al. [25] proposed a DMPC method based on local state information for MAS formation tracking. To the best of our knowledge, there are few studies that apply DMPC to multi-AUV formations. Wei et al. [26] developed a Lyapunov-based distributed predictive controller for AUV formation tracking, subject to current disturbances. The auxiliary controller was utilized to establish stability constraints to ensure the closed-loop stability of the system. However, this method only considers horizontal formations without uncertainties and state constraints. Furthermore, many works that design predictive controllers result in additional computational loads, which could impair the real-time execution capability of the controller. Shen and Shi [27] managed to reduce the MPC computational burden by decomposing the original AUV trajectory tracking optimization problem into smaller subproblems and then solving them in a distributed manner. Despite these efforts, there has been no research to address the heavy computation of DMPC applied to AUV formations. In order to improve the dynamic response and control accuracy of AUV formation tracking in three-dimensional (3-D) space, we adopt the Laguerre orthogonal function to reduce the computational load. In response to these discussions, it is imperative to develop a safe and efficient formation control scheme to solve the problems of disturbances, parameter uncertainties, and complicated constraints.

Motivated by the above observations, this paper investigates the collision-free formation tracking of multi-AUVs with compound disturbances under complicated constraints. A novel FTESO-based distributed dual closed-loop model predictive control scheme is proposed. This method satisfies the formation constraints and collision avoidance requirements while compensating for model uncertainties and external disturbances. We incorporate the Laguerre function to alleviate the computational burden of the DMPC optimization problem, also giving corresponding stability analysis. Based on the connected directed topology, comparative simulations under different schemes demonstrate the effectiveness and robustness of our proposed scheme. The main contributions of this paper are as follows:


The rest of this paper is organized as follows: Section 2 introduces some notations, lemmas, and graph theory, and describes the AUV model and control objective. Section 3 presents the methodology, including the design of the FTESO and dual closed-loop DMPC scheme, the application of the Laguerre function, and the corresponding stability analysis. Sections 4 and 5, respectively, provide simulation results and conclusions.

### **2. Preliminaries**

#### *2.1. Notations and Lemmas*

**Notation.** R*<sup>n</sup>* represents the *n*-dimensional Euclidean space, and R*m*×*<sup>n</sup>* denotes the set of (*m* × *n*) real matrix. *In*, **0***n*, and **0***p*×*<sup>q</sup>* signify (*n* × *n*) identity matrix, (*n* × *n*), and (*p* × *q*) null matrices, respectively. · refers to the Euclidean vector norm and the induced matrix norm, while the infinity norm is denoted by ·∞. *λ*min(·) represents the minimum eigenvalue of the specified matrix (·), with its maximum eigenvalue denoted as *λ*max(·). For simplicity, some notations are defined as *sigp*(*x*) <sup>=</sup> sign(*x*)|*x*<sup>|</sup> *p* , |*x*| *<sup>p</sup>* = ' |*x*1| *p* , |*x*2| *p* ,..., |*xn*| *<sup>p</sup>*(*<sup>T</sup>* , *x* = [*x*1, *x*2,..., *xn*] *<sup>T</sup>*, *<sup>p</sup>* <sup>∈</sup> <sup>R</sup>. sign(·) symbolizes the signum function with sign(0) <sup>=</sup> 0. Notably, *sig*0(*x*) <sup>=</sup> sign(*x*), *sig*0(*x*)|*x*<sup>|</sup> *<sup>p</sup>* = *sigp*(*x*).

**Lemma 1** ([31])**.** Consider the system *. <sup>x</sup>*(*t*) <sup>=</sup> *<sup>f</sup>*(*x*(*t*)), *<sup>x</sup>*(0) <sup>=</sup> *<sup>x</sup>*0, *<sup>f</sup>*(0) <sup>=</sup> 0, *<sup>x</sup>* <sup>∈</sup> <sup>R</sup>*n*, where *<sup>f</sup>* : *<sup>U</sup>* <sup>→</sup> <sup>R</sup>*<sup>n</sup>* is a continuous function. Suppose that this system has a unique solution in forward time for all initial conditions. If there exists a Lyapunov function *V*(*x*), with *V*(*x*0) denoting its initial value, the following can be assumed: (1) The trajectory of this system is finite-time uniformly ultimately bounded stable within the region of *Q*<sup>1</sup> = & *x*|*V*(*x*) *<sup>α</sup>*1−*α*<sup>2</sup> < *<sup>β</sup>*<sup>2</sup> *γ*1 ! , if . *V*(*x*) ≤ −*β*1*V*(*x*) *<sup>α</sup>*<sup>1</sup> + *<sup>β</sup>*2*V*(*x*) *<sup>α</sup>*<sup>2</sup> for *<sup>α</sup>*<sup>1</sup> > *<sup>α</sup>*2, *<sup>β</sup>*<sup>1</sup> > 0, *<sup>β</sup>*<sup>2</sup> > 0, *γ*<sup>1</sup> ∈ (0, *β*1). The settling time for the states reaching the stable residual set is subject to the constraint as *<sup>T</sup>*<sup>1</sup> <sup>≤</sup> *<sup>V</sup>*(*x*0) 1−*α*1 (*β*1−*γ*1)(1−*α*1). (2) The trajectory of this system is fast finite-time uniformly ultimately bounded stable within *Q*<sup>2</sup> = & *x*|*γ*1*V*(*x*) *<sup>α</sup>*1−*α*<sup>2</sup> <sup>+</sup> *<sup>γ</sup>*2*V*(*x*) <sup>1</sup>−*α*<sup>2</sup> < *<sup>β</sup>*<sup>3</sup> ! , if . *V*(*x*) ≤ −*β*1*V*(*x*) *<sup>α</sup>*<sup>1</sup> <sup>−</sup> *<sup>β</sup>*2*V*(*x*) <sup>+</sup> *<sup>β</sup>*3*V*(*x*) *<sup>α</sup>*<sup>2</sup> for *<sup>β</sup>*<sup>3</sup> <sup>&</sup>gt; 0, *<sup>γ</sup>*<sup>2</sup> <sup>∈</sup> (0, *<sup>β</sup>*2). The convergence time T2 is bounded as *T*<sup>2</sup> ≤ ln (*β*2−*γ*2)*V*(*x*0) <sup>1</sup>−*α*1/(*β*1−*γ*1)+<sup>1</sup> (*β*2−*γ*2)(1−*α*1) .

#### *2.2. Graph Theory*

We introduce a directed topology graph *G* = {*V*,*ε*} to describe the information interactions among the AUVs. Let the node set *V* = {*V*1, *V*2, ··· , *VN*} to represent the *N* members in the formation, and an edge set *ε* ⊆ *V* × *V* to represent the communication from the node *Vi* to the node *Vj*. *A* = ' *aij*( <sup>⊂</sup> <sup>R</sup>*N*×*<sup>N</sup>* is defined as an adjacency matrix, where *aij* represents the connection weight and *aij* = 1 if (*i*, *j*) ∈ *ε*, while *aij* = 0 if (*i*, *j*) ∈/ *ε*. It is assumed that the *i*th vehicle could receive information from the virtual leader and its neighbors N*<sup>i</sup>* = {*j* ∈ *V* : (*j*, *i*) ∈ *ε*}. The graph is termed an undirected graph if bidirectional communication links exist among all members of the formation. Otherwise, it is referred to as a directed graph. A directed graph is considered strongly connected if a directed path can connect any point in the formation to any other.

## *2.3. AUV Model*

As shown in Figure 1, it is convenient to describe the six-degree-of-freedom (DOF) AUVs with two reference frames: an earth-fixed frame {*E*} and a body-fixed frame {*B*}. This paper employs a fully actuated torpedo-type AUV, referenced from [32], based on the control objectives. In addition, the AUV uses an ultra-short baseline acoustic positioning system for underwater localization. Since this AUV can be regarded as a highly metacentric stable vehicle with self-stable roll motion, the effect of roll is ignored (roll angle *Φ<sup>i</sup>* = 0, roll angular velocity *pi* = 0). The kinematics and dynamics of the *i*th AUV are described as follows [33]: *.*

$$
\dot{\eta}\_i = J(\eta\_i)\upsilon\_i \tag{1}
$$

$$\mathbf{M}\_i \dot{\boldsymbol{\upsilon}}\_i + \mathbf{C}\_i(\boldsymbol{\upsilon}\_i) \boldsymbol{\upsilon}\_i + \mathbf{D}\_i(\boldsymbol{\upsilon}\_i) \boldsymbol{\upsilon}\_i + \mathbf{g}\_i(\boldsymbol{\eta}\_i) = \boldsymbol{\pi}\_i + \boldsymbol{\pi}\_{\rm ic} \tag{2}$$

where *i* = 1, 2, ... , *N*, *η<sup>i</sup>* = [*xi*, *yi*, *zi*, *θi*, *ψi*] *<sup>T</sup>* <sup>∈</sup> <sup>R</sup>5, and *<sup>v</sup><sup>i</sup>* <sup>=</sup> [*ui*, *vi*, *wi*, *qi*,*ri*] *<sup>T</sup>* <sup>∈</sup> <sup>R</sup><sup>5</sup> denote the states of position, orientation, and velocity of the AUV, respectively. *J*(*ηi*) is a rotation transformation matrix from the body-fixed frame to the earth-fixed frame, expressed as:


*M<sup>i</sup>* represents the inertial matrix, which includes added mass. *Ci*(*vi*) and *Di*(*vi*) denote the Coriolis and centripetal and hydrodynamic damping matrix, respectively, while *gi* (*ηi*) represents the restoring force and moment generated by gravity and buoyancy. *τ<sup>i</sup>* = ' *<sup>τ</sup>iu*, *<sup>τ</sup>iv*, *<sup>τ</sup>iw*, *<sup>τ</sup>iq*, *<sup>τ</sup>ir*(*<sup>T</sup>* represents the control input, and *<sup>τ</sup>ic* denotes the external disturbance. Detailed expressions of these matrices are available in [34].

**Figure 1.** AUV coordinate system.

In practical engineering, we may not be able to obtain accurate hydrodynamic coefficients in the model, so the matrices in (2) are typically divided into two parts: the nominal value part and the uncertainty part caused by linear shifts, i.e., *M<sup>i</sup>* = *M*<sup>∗</sup> *<sup>i</sup>* + Δ*Mi*, *Ci*(*vi*) = *C*<sup>∗</sup> *<sup>i</sup>* (*vi*) + Δ*Ci*(*vi*), *Di*(*vi*) = *D*<sup>∗</sup> *<sup>i</sup>* (*vi*) + Δ*Di*(*vi*), and *g<sup>i</sup>* (*ηi*) = *g*<sup>∗</sup> *<sup>i</sup>* (*ηi*) + Δ*g<sup>i</sup>* (*ηi*), where (·) ∗ *<sup>i</sup>* denotes the nominal value that can be obtained from the computational fluid dynamics (CFD) or experimental analysis. Δ(·)*<sup>i</sup>* symbolizes the difference between the real value and the nominal value.

Accordingly, the *i*th AUV dynamic model (2) can be reformulated as:

$$\mathbf{M}\_{i}^{\*}\dot{\boldsymbol{v}}\_{i} = -\mathbf{C}\_{i}^{\*}(\boldsymbol{v}\_{i})\boldsymbol{v}\_{i} - \mathbf{D}\_{i}^{\*}(\boldsymbol{v}\_{i})\boldsymbol{v}\_{i} - \mathbf{g}\_{i}^{\*}(\boldsymbol{\eta}\_{i}) + \boldsymbol{\tau}\_{i} + \boldsymbol{\tau}\_{id} \tag{4}$$

where *τid* = *τic* − Δ*M<sup>i</sup> . v<sup>i</sup>* − Δ*Ci*(*vi*)*v<sup>i</sup>* − Δ*Di*(*vi*)*v<sup>i</sup>* − Δ*g<sup>i</sup>* (*ηi*) is regarded as the compound disturbance, which includes uncertainties and unknown external disturbance. Typically, external disturbances are periodically varying and energy limited. The model uncertainties are related to the actual states and physical properties of the AUV. Based on the constraints of DMPC on the system state, in practice, we give the following reasonable assumption:

**Assumption 1** ([11])**.** The ocean current disturbance term *<sup>τ</sup>ic* and the first time derivative *. τic* are bounded, and the model uncertainties Δ*Mi*, Δ*Ci*, Δ*Di*, and Δ*g<sup>i</sup>* are unknown and bounded. Hence, the compound disturbance *τid* of the *i*th AUV is bounded and satisfies *τid* <sup>≤</sup> *<sup>τ</sup>id*, where *<sup>τ</sup>id* <sup>∈</sup> <sup>R</sup><sup>+</sup> represents the unknown upper bound.

It should be noted that the above assumption is untenable if there are no system state constraints [35,36].

#### *2.4. Control Objective*

In this paper, the control objective is to develop a control scheme that enables AUV formation to track a reference trajectory while maintaining a predefined configuration. Initially, a FTESO is designed to compensate for external disturbances and model uncertainties of the AUV formation, so that the estimation errors converge to the origin. Subsequently, a dual closed-loop DMPC controller is designed. In this structure, the outer-loop controller enables the *i*th AUV to track the reference trajectory *η<sup>r</sup>* by generating the desired velocity, resulting in the convergence of position tracking errors. The inner-loop controller is used to achieve the convergence of velocity tracking errors. The desired formation is implemented by setting the corresponding formation configuration vector *ri f* and the relative distance vector *rij*. The task must adhere to various constraints and ensure collision avoidance. Because the navigation trajectory has a limited range and the speed is continuous without abrupt changes, we adopt the following reasonable assumptions to avoid singularities in the reference trajectory:

**Assumption 2.** The reference trajectory *η<sup>r</sup>* = [*xr*, *yr*, *zr*, *θr*, *ψr*] *<sup>T</sup>* and its derivatives are smooth and bounded, i.e., *ηr*<sup>∞</sup> ≤ *ηr*, \$ \$ *. ηr* \$ \$ <sup>∞</sup> ≤ *<sup>η</sup>r*1, and \$ \$*.. ηr* \$ \$ <sup>∞</sup> ≤ *ηr*<sup>2</sup> with positive numbers *ηr*, *ηr*1, and *ηr*2.

#### **3. Methodology**

This section develops the FTESO-based distributed dual closed-loop model predictive control scheme for the AUV formation to perform trajectory tracking. A novel FTESO is designed to compensate the compound disturbances. Based on the model information reconstructed by FTESO, the DMPC optimization problems are formulated for the outer and inner loops under constraints such as actuator saturation and collision avoidance, respectively. The Laguerre function is applied to alleviate the computational load. The block diagram of proposed control scheme is depicted in Figure 2.

**Figure 2.** The FTESO-based DMPC dual closed-loop structure for the AUV formation.

#### *3.1. FTESO Design and Convergence Analysis*

The AUV model is fundamental to controller design, but obtaining an accurate model in practice is challenging. Considering the superiority and effectiveness of the ESO technique in estimating and compensating for synthetic uncertainty, a novel fast FTESO is designed to simultaneously reconstruct the external disturbance and model uncertainties of multiple AUVs.

First, define the auxiliary velocity variable as *ωi*(*vi*) = *M*<sup>∗</sup> *<sup>i</sup> v<sup>i</sup>* + ) *vi*, the derivative of *ωi*(*vi*) with respect to time can be obtained from (4)

$$
\omega\_i(\boldsymbol{\upsilon}\_i) = \boldsymbol{\upsilon}\_i - \mathbf{C}\_i^\*(\boldsymbol{\upsilon}\_i)\boldsymbol{\upsilon}\_i - D\_i^\*(\boldsymbol{\upsilon}\_i)\boldsymbol{\upsilon}\_i - \mathbf{g}\_i^\*(\boldsymbol{\eta}\_i) + \boldsymbol{\tau}\_i + \boldsymbol{\tau}\_{id}.\tag{5}
$$

For simplicity, denote *Gi*(*ηi*, *vi*) = *v<sup>i</sup>* − *C*<sup>∗</sup> *<sup>i</sup>* (*vi*)*v<sup>i</sup>* − *D*<sup>∗</sup> *<sup>i</sup>* (*vi*)*v<sup>i</sup>* − *g*<sup>∗</sup> *<sup>i</sup>* (*ηi*). Then, a new variable is defined as *zi*<sup>1</sup> = *ωi*(*vi*), and the order of the system is extended by additional state variables, *<sup>z</sup>i*<sup>2</sup> and *<sup>z</sup>i*3, defined as *<sup>z</sup>i*<sup>2</sup> <sup>=</sup> *<sup>τ</sup>id* and *<sup>z</sup>i*<sup>3</sup> <sup>=</sup> *. <sup>z</sup>i*<sup>2</sup> with *. zi*<sup>3</sup> = *σi*. It should be noted that the compound disturbances *zi*<sup>2</sup> are assumed to be bounded and continuously differentiable, and the components of its second derivative satisfies *σip* <sup>≤</sup> *<sup>σ</sup>i*, *<sup>p</sup>* <sup>=</sup> 1, 2, ... , 5. where *σ<sup>i</sup>* is an unknown positive constant. Afterward, the dynamic model of the *i*th AUV can be extended as follows:

$$\begin{cases}
\dot{z}\_{i1} = \mathcal{G}\_i(\eta\_i, \sigma\_i) + \mathbf{r}\_i + \mathbf{z}\_{i2} \\
\dot{z}\_{i2} = \mathbf{z}\_{i3} \\
\dot{z}\_{i3} = \sigma\_i.
\end{cases} \tag{6}$$

Denote *z*ˆ*i*1, *z*ˆ*i*2, and *z*ˆ*i*<sup>3</sup> as the observation values of states *zi*1, *zi*2, and *zi*<sup>3</sup> in the above extended system, and *ei*<sup>1</sup> = *z*ˆ*i*<sup>1</sup> − *zi*1, *ei*<sup>2</sup> = *z*ˆ*i*<sup>2</sup> − *zi*2, and *ei*<sup>3</sup> = *z*ˆ*i*<sup>3</sup> − *zi*<sup>3</sup> as the observation errors of the velocity, the compound disturbances, and its first derivatives, respectively. Then, a third-order fast FTESO is proposed as follows:

$$\begin{cases}
\dot{\hat{\mathbf{z}}}\_{i1} = \dot{\mathbf{z}}\_{i2} - \beta\_{i1}(sig^{a\_{i1}}(\mathbf{e}\_{i1}) + \mathbf{e}\_{i1}) + G\_{i}(\eta\_{i\prime}\mathbf{e}\_{i\prime}) + \mathbf{\tau}\_{i} \\
\dot{\hat{\mathbf{z}}}\_{i2} = \dot{\mathbf{z}}\_{i3} - \beta\_{i2}(sig^{a\_{i2}}(\mathbf{e}\_{i1}) + 2sig^{a\_{i1}}(\mathbf{e}\_{i1}) + \mathbf{e}\_{i1}) \\
\dot{\hat{\mathbf{z}}}\_{i3} = -\beta\_{i3}(sig^{a\_{i3}}(\mathbf{e}\_{i1}) + 2sig^{a\_{i2}}(\mathbf{e}\_{i1}) + sig^{a\_{i1}}(\mathbf{e}\_{i1}))
\end{cases} \tag{7}$$

where the observer gains satisfy *βik* > 0, *k* = 1, 2, 3, *αi*<sup>1</sup> ∈ (2/3, 1) and *αi*<sup>2</sup> = 2*αi*<sup>1</sup> − 1, and *αi*<sup>3</sup> = 3*αi*<sup>1</sup> − 2. Although the actual value of *zik* is probably unavailable, its observed value *z*ˆ*ik* can be obtained by the above FTESO. The analysis and proof that *z*ˆ*ik* tracks the actual value are described below.

According to the extended system (6) and the proposed FTESO (7), we can obtain the observation error dynamics as follows:

$$\begin{cases}
\dot{\mathbf{e}}\_{i1} = -\beta\_{i1}(sig^{a\_{i1}}(\mathbf{e}\_{i1}) + \mathbf{e}\_{i1}) + \mathbf{e}\_{i2} \\
\dot{\mathbf{e}}\_{i2} = -\beta\_{i2}(sig^{a\_{i2}}(\mathbf{e}\_{i1}) + 2sig^{a\_{i1}}(\mathbf{e}\_{i1}) + \mathbf{e}\_{i1}) + \mathbf{e}\_{i3} \\
\dot{\mathbf{e}}\_{i3} = -\beta\_{i3}(sig^{a\_{i3}}(\mathbf{e}\_{i1}) + 2sig^{a\_{i2}}(\mathbf{e}\_{i1}) + sing^{a\_{i1}}(\mathbf{e}\_{i1})) - \sigma\_i.
\end{cases} \tag{8}$$

The stability and convergence of the proposed FTESO are stated in the following theorem:

**Theorem 1.** *Consider the AUV formation control system with the dynamic model (4) under Assumption 1. If the FTESO is proposed in the form of (9), with appropriate observer gains satisfying the prescribed constraints, then the observation errors e<sup>i</sup>* = ' *eT <sup>i</sup>*1, *<sup>e</sup><sup>T</sup> <sup>i</sup>*2, *<sup>e</sup><sup>T</sup> i*3 (*<sup>T</sup> will converge to the small region* Ω*<sup>i</sup> in finite time Ti f . This implies that the error dynamics system (8) is finite-time uniformly ultimately bounded stable*.

**Proof of Theorem 1.** Consider a Lyapunov candidate function as *Vi*1(*e*) = *ε<sup>T</sup> <sup>i</sup> Piεi*, where *P<sup>i</sup>* is a positive definite symmetric matrix and *ε<sup>T</sup> <sup>i</sup>* = (*sigαi*<sup>1</sup> (*ei*1) + *ei*1) *<sup>T</sup>*, *e<sup>T</sup> <sup>i</sup>*2, *<sup>e</sup><sup>T</sup> i*3 is introduced as an auxiliary error variable. It should be noted that *ei*1, *ei*2, and *ei*<sup>3</sup> will converge to origin in finite time, if the new state *ε<sup>i</sup>* is finite-time stable. The time derivative of *εi*, invoking (8), yields:

*. ε<sup>i</sup>* = ⎡ <sup>⎣</sup> *<sup>α</sup>i*1|*ei*1<sup>|</sup> *<sup>α</sup>i*1−<sup>1</sup> *. <sup>e</sup>i*<sup>1</sup> <sup>+</sup> *. <sup>e</sup>i*<sup>1</sup> *. ei*2 *. ei*3 ⎤ <sup>⎦</sup> <sup>=</sup> ⎡ <sup>⎣</sup> *<sup>α</sup>i*1|*ei*1<sup>|</sup> *αi*1−1 (*ei*<sup>2</sup> <sup>−</sup> *<sup>β</sup>i*1(*sigαi*<sup>1</sup> (*ei*1) <sup>+</sup> *<sup>e</sup>i*1)) *<sup>e</sup>i*<sup>3</sup> <sup>2</sup> <sup>−</sup> *<sup>β</sup>i*2(*sigαi*<sup>2</sup> (*ei*1) <sup>+</sup> *sigαi*<sup>1</sup> (*ei*1)) <sup>−</sup>*βi*3(*sigαi*<sup>3</sup> (*ei*1) <sup>+</sup> *sigαi*<sup>2</sup> (*ei*1)) ⎤ ⎦ + ⎡ ⎣ *<sup>e</sup>i*<sup>2</sup> <sup>−</sup> *<sup>β</sup>i*1(*sigαi*<sup>1</sup> (*ei*1) <sup>+</sup> *<sup>e</sup>i*1) *<sup>e</sup>i*<sup>3</sup> <sup>2</sup> <sup>−</sup> *<sup>β</sup>i*2(*sigαi*<sup>1</sup> (*ei*1) <sup>+</sup> *<sup>e</sup>i*1) <sup>−</sup>*βi*3(*sigαi*<sup>2</sup> (*ei*1) <sup>+</sup> *sigαi*<sup>1</sup> (*ei*1)) ⎤ <sup>⎦</sup> <sup>+</sup> ⎡ ⎣ **0**5 **0**5 −*σ<sup>i</sup>* ⎤ <sup>⎦</sup> <sup>=</sup> *diag*|*ei*1<sup>|</sup> *αi*1−1 , |*ei*1| *αi*1−1 , |*ei*1| *αi*1−1 *Ai*1*ε<sup>i</sup>* <sup>+</sup> *<sup>A</sup>i*2*ε<sup>i</sup>* <sup>+</sup> *<sup>Φ</sup><sup>i</sup>* (9)

where *Φ<sup>i</sup>* = ' **0**<sup>5</sup> **0**<sup>5</sup> −*σ<sup>i</sup>* (*<sup>T</sup>* and the coefficient matrices *<sup>A</sup>i*<sup>1</sup> and *<sup>A</sup>i*<sup>2</sup> are expressed as:

$$A\_{i1} = \begin{pmatrix} -a\_{i1}\boldsymbol{\beta}\_{i1}\mathbf{I}\_{5} & a\_{i1}\mathbf{I}\_{5} & \mathbf{0}\_{5} \\ -\beta\_{i2}\mathbf{I}\_{5} & \mathbf{0}\_{5} & \mathbb{E}\_{i}^{-1}\mathbf{I}\_{5}/2 \\ -\beta\_{i3}\mathbf{\overline{\tau}}\_{i}\mathbf{I}\_{5} & \mathbf{0}\_{5} & \mathbf{0}\_{5} \end{pmatrix}, \quad A\_{2i} = \begin{pmatrix} -\beta\_{i1}\mathbf{I}\_{5} & \mathbf{I}\_{5} & \mathbf{0}\_{5} \\ -\beta\_{i2}\mathbf{I}\_{5} & \mathbf{0}\_{5} & \mathbf{I}\_{5}/2 \\ -\beta\_{i3}\mathbf{\overline{\tau}}\_{i}\mathbf{I}\_{5} & \mathbf{0}\_{5} & \mathbf{0}\_{5} \end{pmatrix} \tag{10}$$

with *ei* = |*ei*1| *αi*1−1 . From the characteristic polynomials of *Ai*<sup>1</sup> and *Ai*<sup>2</sup> that all their eigenvalues have negative real parts if the observer gains are set as *βik* > 0, indicating that *Ai*<sup>1</sup> and *Ai*<sup>2</sup> are Hurwitz matrices. Thus, symmetric and positive definite matrices *Qi*<sup>1</sup> and *Qi*<sup>2</sup> exist that satisfy the following Lyapunov equations:

$$\begin{cases} A\_{i1}^T P\_i + P\_i A\_{i1} = -\mathcal{Q}\_{i1} \\ A\_{i2}^T P\_i + P\_i A\_{i2} = -\mathcal{Q}\_{i2} \end{cases} \tag{11}$$

Differentiating *Vi*1(*e*) with respect to time yields the following:

$$\begin{split} \dot{V}\_{i1} &= \mathfrak{e}\_{i}^{T} \Big[ \text{diag} \left( [\overline{\boldsymbol{\varepsilon}}\_{i}, \overline{\boldsymbol{\varepsilon}}\_{i} \boldsymbol{\varepsilon}\_{i}] \right) \big( \mathbf{A}\_{i1}^{T} \mathbf{P}\_{i} + \mathbf{P}\_{i} \mathbf{A}\_{i1} \big) \big] \boldsymbol{\varepsilon}\_{i} + \mathfrak{e}\_{i}^{T} \big( \mathbf{A}\_{i2}^{T} \mathbf{P}\_{i} + \mathbf{P}\_{i} \mathbf{A}\_{i2} \big) \boldsymbol{\varepsilon}\_{i} + 2 \mathfrak{e}\_{i}^{T} \mathbf{P}\_{i} \Phi\_{i} \\ &= -\mathfrak{e}\_{i}^{T} \big[ \text{diag} \left( [\overline{\boldsymbol{\varepsilon}}\_{i}, \overline{\boldsymbol{\varepsilon}}\_{i}, \overline{\boldsymbol{\varepsilon}}\_{i}] \big) \mathbf{Q}\_{i1} \big] \boldsymbol{\varepsilon}\_{i} - \mathfrak{e}\_{i}^{T} \mathbf{Q}\_{i2} \boldsymbol{\varepsilon}\_{i} + 2 \mathfrak{e}\_{i}^{T} \mathbf{P}\_{i} \Phi\_{i} \leq -\overline{\boldsymbol{\varepsilon}}\_{i}^{\max} \mathfrak{e}\_{i}^{T} \mathbf{Q}\_{i1} \mathbf{e}\_{i} - \mathfrak{e}\_{i}^{T} \mathbf{Q}\_{i2} \boldsymbol{\varepsilon}\_{i} + 2 \| \mathbf{e}\_{i} \| \big\| \mathbf{P}\_{i} \| \| \mathbf{P}\_{i} \| \end{split} \tag{12}$$

where *e*max *<sup>i</sup>* = |*ei*1| *αi*1−1 max and |*ei*1|max = max{|*ei*11|,..., |*ei*15|}. Given the fact that |*ei*1|max ≤ *ei*1 <sup>≤</sup> *εi*1/*αi*<sup>1</sup> and *<sup>α</sup>i*<sup>1</sup> <sup>∈</sup> 2 <sup>3</sup> , 1 , we can obtain the following:

$$\begin{split} \dot{V}\_{i1} &\leq -||\boldsymbol{\varepsilon\_{i}}||^{\frac{a\_{i1}-1}{a\_{i1}}} \boldsymbol{\varepsilon\_{i}^{T}} \boldsymbol{Q}\_{i1} \boldsymbol{\varepsilon\_{i}} - \boldsymbol{\varepsilon\_{i}^{T}} \boldsymbol{Q}\_{i2} \boldsymbol{\varepsilon\_{i}} + 2||\boldsymbol{\varepsilon\_{i}}|| ||\boldsymbol{P\_{i}}|| ||\boldsymbol{\Phi\_{i}}|| \\ &\leq -\lambda\_{\text{min}}(\boldsymbol{Q}\_{i1}) ||\boldsymbol{\varepsilon\_{i}}||^{3-\frac{1}{a\_{i1}}} - \lambda\_{\text{min}}(\boldsymbol{Q}\_{i2}) ||\boldsymbol{\varepsilon\_{i}}||^{2} + 2||\boldsymbol{\varepsilon\_{i}}|| ||\boldsymbol{P\_{i}}|| ||\boldsymbol{\Phi\_{i}}||. \end{split} \tag{13}$$

Since *σ<sup>i</sup>* is assumed to be bounded reasonably by *σip* <sup>≤</sup> *<sup>σ</sup>i*, we have 2*εiPiΦi* <sup>≤</sup> 2 <sup>√</sup>5*σiεiPi* <sup>≤</sup> <sup>2</sup> <sup>√</sup>5*σiλ*min(*Pi*) − 1 <sup>2</sup> *V* 1 2 *<sup>i</sup>*<sup>1</sup> *Pi*, by using the inequality

$$\left\|\lambda\_{\min}(\mathbf{P}\_{i})\right\|\left\|\varepsilon\_{i}\right\|\right\|^{2} \leq \left|V\_{i1} \leq \lambda\_{\max}(\mathbf{P}\_{i})\right\|\left\|\varepsilon\_{i}\right\|^{2} \tag{14}$$

Then, inequality (13) becomes the following:

$$\begin{split} \dot{V}\_{l1} &\leq -\lambda\_{\text{min}}(\mathbf{Q}\_{i1})\lambda\_{\text{max}}(\mathbf{P}\_{i})^{\frac{1}{2\bar{\alpha}\_{l1}}-\frac{3}{2}}V\_{l1}^{\frac{1}{2}-\frac{1}{2\bar{\alpha}\_{l1}}} - \lambda\_{\text{min}}(\mathbf{Q}\_{i2})\lambda\_{\text{max}}(\mathbf{P}\_{i})^{-1}V\_{l1} + 2\sqrt{\mathbf{S}}\dot{r}\_{i} \|\mathbf{P}\_{i}\|\lambda\_{\text{min}}(\mathbf{P}\_{i})^{-\frac{1}{2}}V\_{i1}^{\frac{1}{2}} \\ &\leq -\lambda\_{l1}V\_{l1}^{\frac{3}{2}-\frac{1}{2\bar{\alpha}\_{l1}}} - \lambda\_{l2}V\_{l1} + \lambda\_{l3}V\_{l1}^{\frac{1}{2}} \end{split} \tag{15}$$

where *<sup>λ</sup>i*<sup>1</sup> <sup>=</sup> <sup>−</sup>*λ*min(*Qi*1)*λ*max(*Pi*) <sup>1</sup> 2*αi*1 − 3 <sup>2</sup> , *λi*<sup>2</sup> = −*λ*min(*Qi*2)*λ*max(*Pi*) −1 , and *λi*<sup>3</sup> = 2 <sup>√</sup>5*σiPiλ*min(*Pi*) − 1 2 .

It can be seen that (15) has the same form as the sufficient condition in Lemma 1 2. Thus, the error trajectories of the proposed FTESO (7) are fast finite-time uniformly ultimately bounded stable. The state observation errors *e<sup>i</sup>* will converge to a small region Ω*<sup>i</sup>* in the finite time *Ti f* . Moreover, the settling time *Ti f* is subject to the constraint:

$$T\_{if} \le \frac{\ln\left(\left(\lambda\_{i2} - \overline{\lambda}\_{i2}\right) V\_{i1} (\mathfrak{e}\_0)^{\frac{1}{2\overline{\alpha}\_{i1}} - \frac{1}{2}} / \left(\lambda\_{i1} - \overline{\lambda}\_{i1}\right) + 1\right)}{\left(\lambda\_{i2} - \overline{\lambda}\_{i2}\right) \left(\frac{1}{2\overline{\alpha}\_{i1}} - \frac{1}{2}\right)}.\tag{16}$$

And the stable region Ω*<sup>i</sup>* is denoted as

$$\Omega\_i = \left\{ \mathbf{e} \middle| \overline{\lambda}\_{i1} V\_{i1} (\mathbf{e})^{1 - \frac{1}{2\overline{\alpha}\_{i1}}} + \overline{\lambda}\_{i2} V\_{i1} (\mathbf{e})^{\frac{1}{2}} < \lambda\_{i3} \right\} \tag{17}$$

where *λi*<sup>1</sup> and *λi*<sup>2</sup> are arbitrary constants that meet the conditions *λi*<sup>1</sup> ∈ (0, *λi*1) and *λi*<sup>2</sup> ∈ (0, *λi*2). This completes the proof. -

**Remark 1.** Contrasting our proposed FTESO (7) with the FTESO in [37], our approach factors in the dynamics of disturbances and uncertainties to achieve a higher degree of estimation accuracy. Our usage of fractional powers within the FTESO allows for a quick finite-time convergence. It can be noted that the size of the attraction region Ω*<sup>i</sup>* hinges upon the selection of the observer gains *βik* and *αi*1. By increasing *βik* or decreasing *αi*1, the attraction region of the observation error system can be expanded and the convergence speed can be improved, but excessive tuning will lead to undesired overshoot and oscillation. As a result, a trade-off should be taken for *βik* and *αi*1.

#### *3.2. Outer-Loop Formation Prediction Control Law*

In this subsection, we design a DMPC-based outer-loop formation controller. This controller, which draws on the information interaction with neighbors, facilitates the positional tracking of the *i*th AUV. The controller operates under composite constraints and ensures the avoidance of collisions. Then, we formulate a constrained QP problem in accordance with the control objective to obtain the optimal driving speed.

To facilitate the recursive model prediction and the implementation of the control law, the kinematic model (1) is discretized by using the Forward-Euler method with a sampling period *Ts*, resulting in following discrete model:

$$
\eta\_i(k+1) = \eta\_i(k) + J\_i(k)\sigma\_i(k)T\_s. \tag{18}
$$

To smoothen the speed change of the AUV, the velocity increment Δ*uiv*(*k*) = *vi*(*k*) − *vi*(*k* − 1) is taken as the control input. *xiη*(*k*) = ' *ηi*(*k*) *vi*(*k* − 1) (*<sup>T</sup>* is denoted as the state variable of the prediction model. The augmented state-space model of the outer-loop subsystem can be derived as:

$$\mathbf{x}\_{i\eta}(k+1) = A\_{i\eta}\mathbf{x}\_{i\eta}(k) + B\_{i\eta}\Delta\mathbf{u}\_{i\upsilon}(k) \tag{19}$$

$$\mathbf{y}\_{i\eta}(k) = \mathbf{C}\_{i\eta} \mathbf{x}\_{i\eta}(k) \tag{20}$$

where *Ai<sup>η</sup>* = *I*<sup>5</sup> *Ji*(*k*)*Ts* **0**<sup>5</sup> *I*<sup>5</sup> <sup>∈</sup> <sup>R</sup>10×10, *<sup>B</sup>i<sup>η</sup>* <sup>=</sup> *Ji*(*k*)*Ts I*5 <sup>∈</sup> <sup>R</sup>10×5, and *<sup>C</sup>i<sup>η</sup>* <sup>=</sup> ' *I*<sup>5</sup> **0**<sup>5</sup> ( <sup>∈</sup> <sup>R</sup>5×10.

According to the state prediction model (19) and (20), we can calculate the predicted state sequence of the system when given an input sequence. Let *Np*<sup>1</sup> and *Nc*<sup>1</sup> denote the prediction and control horizon of the outer-loop controller, respectively. The predicted state sequence and the input incremental sequence are usually represented by compact vectors:

$$\mathbf{Y}\_{i\eta} = \begin{bmatrix} \mathbf{y}\_{i\eta}(k+1|k) \\ \mathbf{y}\_{i\eta}(k+2|k) \\ \vdots \\ \mathbf{y}\_{i\eta}(k+N\_{p1}|k) \end{bmatrix} \in \mathbb{R}^{5N\_{p1}}, \mathbf{x}\_{i\eta} = \begin{bmatrix} \mathbf{x}\_{i\eta}(k+1|k) \\ \mathbf{x}\_{i\eta}(k+2|k) \\ \vdots \\ \mathbf{x}\_{i\eta}\left(k+N\_{p1}|k\right) \end{bmatrix} \in \mathbb{R}^{10N\_{p1}} \tag{21}$$

$$
\Delta \mathcal{U}\_{\dot{\upsilon}\upsilon} = \begin{bmatrix}
\Delta \mu\_{\dot{\upsilon}\upsilon}(k|k) \\
\Delta \mu\_{\dot{\upsilon}\upsilon}(k+1|k) \\
\vdots \\
\Delta \mu\_{\dot{\upsilon}\upsilon}(k+N\_{\ell 1}-1|k)
\end{bmatrix} \in \mathbb{R}^{5N\_{\ell 1}} \tag{22}
$$

where *yiη*(*k* + *l*|*k*) and *xiη*(*k* + *l*|*k*) are the output vector *yiη*(*k* + *l*) and state vector *xiη*(*k* + *l*) predicted at time *k*, respectively. Δ*uiv*(*k* + *j*|*k*) denotes the input increment Δ*uiv*(*k* + *j*) predicted at the same time *k*. Then, we characterize the relationship between the predicted output vector sequence and the control increment sequence through the following prediction equation based on the recurrence relations:

$$\mathbf{Y}\_{i\eta} = H\_{i\mathbf{x}}^1 \mathbf{x}\_{i\eta}(k) + H\_{i\mathbf{u}}^1 \Delta \mathbf{U}\_{i\upsilon} \tag{23}$$

where *xiη*(*k*) is the initial state, *H*<sup>1</sup> *ix* = *CiηAiη*,*CiηA*<sup>2</sup> *iη*,...,*CiηANp*<sup>1</sup> *iη T* <sup>∈</sup> <sup>R</sup>5*Np*1×<sup>10</sup> and

$$H\_{iu}^{1} = \begin{bmatrix} \mathbf{C}\_{i\eta}\mathbf{B}\_{i\eta} & \mathbf{0}\_{5} & \cdots & \mathbf{0}\_{5} \\ \mathbf{C}\_{i\eta}\mathbf{A}\_{i\eta}\mathbf{B}\_{i\eta} & \mathbf{C}\_{i\eta}\mathbf{B}\_{i\eta} & \cdots & \mathbf{0}\_{5} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{C}\_{i\eta}\mathbf{A}\_{i\eta}^{N\_{p1}-1}\mathbf{B}\_{i\eta} & \mathbf{C}\_{i\eta}\mathbf{A}\_{i\eta}^{N\_{p1}-2}\mathbf{B}\_{i\eta} & \cdots & \mathbf{C}\_{i\eta}\mathbf{A}\_{i\eta}^{N\_{p1}-N\_{c1}}\mathbf{B}\_{i\eta} \end{bmatrix} \in \mathbb{R}^{5N\_{p1}\times 5N\_{c1}}.$$

Considering the control objective, the constraints within the outer-loop subsystem are considered. First, we set upper and lower boundaries for the amplitude of the control input *uiv*(*k*) and the input increment Δ*uiv*(*k*):

$$
\mu\_{iv}^{\min} \le \mu\_{iv}(k) \le \mu\_{iv}^{\max} \tag{24}
$$

$$
\Delta \mathfrak{u}\_{\dot{\upsilon}\upsilon}^{\min} \le \Delta \mathfrak{u}\_{\dot{\upsilon}\upsilon}(k) \le \Delta \mathfrak{u}\_{\dot{\upsilon}\upsilon}^{\max} \tag{25}
$$

where *u*min *iv* and <sup>Δ</sup>*u*min *iv* represent the predefined lower bounds, and *<sup>u</sup>*max *iv* and <sup>Δ</sup>*u*max *iv* represent the predefined upper bounds.

Next, to assure safe navigation throughout the formation construction stage, we need to consider the collision avoidance constraints between AUVs. The primitive collision avoidance constraints of the *i*th AUV can be transformed into a convex constraint, as follows:

$$\left\| \left\| \mathcal{S} \left( \mathcal{y}\_{i\eta}(k+l|k) - \mathcal{y}\_{j\eta}(k+l|k) \right) \right\| \right\| \geq r\_{s\prime} \quad j \in \Xi\_i \tag{26}$$

where *l* = 1, 2, ... , *Np*<sup>1</sup> and *rs* is the preset minimum allowable distance between the *i*th AUV and the *j*th AUV. *S* denotes a scaling matrix. Let *rd* be the radius of the safe detection zone for the *i*th AUV. *Ξ<sup>i</sup>* is the set of those AUVs that contain within *rd*. Let the nominal value *yi<sup>η</sup>* represent an initial guess of the actual value *yi<sup>η</sup>* for convexifying the collision avoidance constraint. It follows from (26) that a sufficient condition for upholding the collision avoidance constraint is the following:

$$\left\| \overline{d}\_{ij}^T (k+l|k) \mathbf{S}^T \mathbf{S} \left( \mathbf{y}\_{i\eta}(k+l|k) - \overline{\mathbf{y}}\_{j\eta}(k+l|k) \right) \right\| \geq r\_s \left\| \mathbf{S} \overline{d}\_{ij}(k+l|k) \right\| \tag{27}$$

where *dij*(*k* + *l*|*k*) = *yiη*(*k* + *l*|*k*) − *yjη*(*k* + *l*|*k*). In order to express the constraints in a compact matrix form, define *R<sup>l</sup>* = *rs* \$ \$ \$ *Sdij*(*k* + *l*|*k*) \$ \$ \$ <sup>+</sup> *<sup>d</sup> T ij*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*)*STSyjη*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*), *<sup>R</sup>ij* <sup>=</sup> *R*1, *R*2,..., *RNp*<sup>1</sup> *T* and *<sup>S</sup>ij* <sup>=</sup> *diag*& *S*1, *S*2,..., *SNp*<sup>1</sup> ! , and *S<sup>l</sup>* = *d T ij*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*)*STS*. Then, (27) can be rewritten as *SijYi<sup>η</sup>* ≥ *Rij*. Substitute (23) to derive the collision avoidance constraint as follows:

$$
\overline{\mathfrak{S}}\_{ij} H^{1}\_{iu} \Delta \mathcal{U}\_{iv} \geq \overline{\mathfrak{R}}\_{ij} - \overline{\mathfrak{S}}\_{ij} H^{1}\_{ix} \mathfrak{x}\_{i\eta}(k). \tag{28}
$$

The input amplitude constraint (24) can be converted to the input incremental constraint, associating (25) and (28), expressed in the compact linear constraint form as follows:

$$
\Gamma\_{i\eta} \Delta \mathcal{U}\_{iv} \le \gamma\_{i\eta} \tag{29}
$$

$$\begin{array}{c} \text{where } I\_{\overline{\eta}} = \begin{bmatrix} I\_{5N\_{\ell 1}} \\ -I\_{5N\_{\ell 1}} \\ I\_{\eta 1} \\ -I\_{\eta 1} \\ -\overline{S}\_{i}\overline{H}\_{i\mu}^{\perp} \end{bmatrix}, \ \gamma\_{i\eta} = \begin{bmatrix} \Delta\varPi\_{\overline{\upsilon}\nu}^{\text{max}} \\ -\Delta\varPi\_{\overline{\upsilon}\nu}^{\text{min}} \\ \Pi\_{i\overline{\upsilon}}^{\text{max}} + I\_{\eta 2}\mu\_{i}(k-1) \\ -\underline{\Pi}\_{i\overline{\upsilon}}^{\text{min}} + I\_{\eta 2}\mu\_{i}(k-1) \\ \overline{S}\_{i\overline{\upsilon}}\overline{H}\_{i\varepsilon}^{\text{x}}\mathbf{x}\_{i\eta}(k) - \overline{R}\_{i\overline{\}} \end{bmatrix}, I\_{\eta 1} = \begin{bmatrix} I\_{5} & \mathbf{0}\_{5} & \cdots & \mathbf{0}\_{5} \\ I\_{5} & I\_{5} & \cdots & \mathbf{0}\_{5} \\ \vdots & \vdots & \ddots & \vdots \\ I\_{5} & I\_{5} & \cdots & I\_{5} \end{bmatrix} \in \mathbb{R}^{5N\_{\ell} \times 5N\_{\ell} \times 10}, \text{ and } I\_{\eta 2} = [I\_{5}, I\_{5}, \dots, I\_{5}]^{T} \in \mathbb{R}^{5N\_{\ell} \times 10}. \end{array}$$

In order to achieve the control objective of formation positional tracking with low energy requirements, we define the local distributed cost function in the outer-loop subsystem of the *i*th AUV in a discretized form:

$$\begin{split} J\_{i\eta}(k) &= \sum\_{l=1}^{N\_{p1}} \left\| \left( y\_{i\eta}(k+l|k) - y\_{if}(k+l) \right) \right\|\_{\mathbf{Q}\_{if}}^{2} + \sum\_{l=0}^{N\_{\rm cl}-1} \left\| \Delta u\_{i\overline{\nu}}(k+l|k) \right\|\_{\mathbf{R}\_{i1}}^{2} \\ &+ \sum\_{l=1}^{N\_{p1}} \sum\_{j \in \mathbf{N}\_{i}} a\_{ij} \left\| \left( y\_{i\eta}(k+l|k) - y\_{ij}(k+l) \right) \right\|\_{\mathbf{Q}\_{ij}}^{2} \end{split} \tag{30}$$

where *Qi f* , *Qij*, and *Ri*<sup>1</sup> are the weight matrices. *yi f*(*k* + *l*) = *ηr*(*k* + *l*) + *ri f*(*k* + *l*) with *ri f*(*k* + *l*) represents the formation configuration. *yij*(*k* + *l*) = *yjη*(*k* + *l*) + *rij*(*k* + *l*) with *rij*(*k* + *l*) represents the predefined relative distance between the *i*th AUV and its neighbor *j*th AUV. *Np*<sup>1</sup> indicates the degree of prediction of future tracking errors. The larger it is, the better the tracking accuracy and stability. The smaller *Nc*<sup>1</sup> is, the worse the dynamic response is, and conversely the more maneuverable the control is. *Qi f* is the position tracking matrix, the larger it is, the better the tracking accuracy and dynamic response. *Qij* is the relative position matrix, the larger it is, the better the ability of the formation to maintain the preset configuration. *Ri*<sup>1</sup> is the control increment weight matrix, mainly to limit the drastic change of Δ*uiv*.

Based on the above derivations, we can formulate the optimization problem for the outer-loop subsystem of the *i*th AUV at the sampling instant *k* within the recedinghorizon framework:

$$\begin{array}{ll}\underset{\Delta \mathcal{U}\_{iv}}{\min} & I\_{i\eta}(k) \\ \text{s.t.} & \Gamma\_{i\eta} \Delta \mathcal{U}\_{iv} \leq \gamma\_{i\eta}. \end{array} \tag{31}$$

To simplify the computation of (31), it can be transformed into a convex QP problem. This problem is solved over a finite receding horizon using a QP solver. The standard convex QP form of the DMPC problem (31) can be derived:

$$\begin{array}{l} \Delta \mathbf{U}\_{iv}^{\*} = \underset{\Delta \mathbf{U}\_{iv}}{\text{argmin}} \left( \frac{1}{2} \Delta \mathbf{U}\_{iv}^{T} \mathbf{W}\_{i\eta} \Delta \mathbf{U}\_{iv} + f\_{i\eta}^{T} \Delta \mathbf{U}\_{iv} \right) \\ \text{s.t. } \Gamma\_{i\eta} \Delta \mathbf{U}\_{iv} \le \gamma\_{i\eta} \end{array} \tag{32}$$

where *Wi<sup>η</sup>* = *Ri*<sup>1</sup> + *H*1*<sup>T</sup> iu <sup>Q</sup>i f <sup>H</sup>*<sup>1</sup> *iu* <sup>+</sup> <sup>∑</sup>*j*∈N*<sup>i</sup> aijH*1*<sup>T</sup> iu <sup>Q</sup>ijH*<sup>1</sup> *iu*,

$$\begin{array}{rclcrcl}f\_{\boldsymbol{\eta}\boldsymbol{\eta}} &=& \boldsymbol{\mathsf{H}}\_{\boldsymbol{\boldsymbol{\text{in}}}}^{\mathrm{T}} \overline{\mathsf{Q}}\_{\boldsymbol{if}} \Big( \boldsymbol{\mathsf{H}}\_{\boldsymbol{\text{in}}}^{\mathrm{I}} \boldsymbol{\mathsf{x}}\_{\boldsymbol{i}\boldsymbol{\eta}} - \boldsymbol{\mathsf{Y}}\_{\boldsymbol{if}} \Big) & + \sum\_{j \in \mathcal{N}\_{i}} a\_{ij} \boldsymbol{\mathsf{H}}\_{\boldsymbol{if}}^{\mathrm{T}} \overline{\mathsf{Q}}\_{\boldsymbol{if}} \Big( \boldsymbol{\mathsf{H}}\_{\boldsymbol{\text{in}}}^{\mathrm{I}} \boldsymbol{\mathsf{x}}\_{\boldsymbol{i}\boldsymbol{\eta}} - \boldsymbol{\mathsf{Y}}\_{\boldsymbol{if}} \Big), & \text{ with } \boldsymbol{\mathsf{Y}}\_{\boldsymbol{if}} &=\\\ a\_{if}(\boldsymbol{k}+1)\_{\boldsymbol{\prime}}, \ldots, \boldsymbol{y}\_{if}(\boldsymbol{k}+N\_{\mathrm{p}1}) \Big)^{T} \boldsymbol{\mathsf{y}}\_{\boldsymbol{i}\boldsymbol{\prime}} &=& \left[ \boldsymbol{y}\_{\boldsymbol{i}i}(\boldsymbol{k}+1)\_{\boldsymbol{\prime}}, \ldots, \boldsymbol{y}\_{\boldsymbol{ii}}(\boldsymbol{k}+N\_{\mathrm{p}1}) \right]^{T} \boldsymbol{\mathsf{Q}}\_{\boldsymbol{if}} &= \end{array}$$

 *yi f*(*k* + 1),..., *yi f k* + *Np*<sup>1</sup> , *Yij* = *yij*(*k* + 1),..., *yij k* + *Np*<sup>1</sup> , *Qi f* = *diag*& *<sup>Q</sup>i f* , *<sup>Q</sup>i f* ,..., *<sup>Q</sup>i f* ! <sup>∈</sup> <sup>R</sup>5*Np*1×5*Np*<sup>1</sup> , *<sup>Q</sup>ij*, and *<sup>R</sup>i*<sup>1</sup> are similar to *<sup>Q</sup>i f* , both corresponding compact matrices.

By solving the QP optimization problem in (32) online, we obtain the optimal control input increment sequence Δ*U*∗ *iv*. Of this sequence, we only utilize the first element Δ*u*<sup>∗</sup> *iv*(*k*|*k*) for receding optimization. Once Δ*u*∗ *iv*(*k*) is determined, we obtain *vi*(*k*) which serves as the desired driving speed for the inner-loop controller of the *i*th AUV, i.e.,

$$
\sigma\_{i\dot{r}}(k) = \sigma\_i(k) = \sigma\_i(k-1) + \Delta \mathfrak{u}\_{i\dot{r}}^\*(k). \tag{33}
$$

#### *3.3. Inner-Loop Formation Prediction Control Law*

In this subsection, with the aid of the proposed FTESO, we design a DMPC-based formation controller for the inner-loop subsystem to obtain the optimal driving force and moment for the *i*th AUV to track the desired speed.

The dynamic model (4) is discretized with a sampling period *Ts*, yielding the following discretized model:

$$\boldsymbol{\sigma}\_{i}(k+1) = \left(\boldsymbol{I} - \boldsymbol{M}\_{i}^{\*-1}\boldsymbol{T}\_{s}(\mathbf{C}\_{i}^{\*} + \boldsymbol{D}\_{i}^{\*})\right)\boldsymbol{\sigma}\_{i}(k) + \boldsymbol{M}\_{i}^{\*-1}\boldsymbol{T}\_{s}\boldsymbol{\pi}\_{i}(k) + \boldsymbol{M}\_{i}^{\*-1}\boldsymbol{T}\_{s}\boldsymbol{\pi}\_{i\boldsymbol{d}}(k)\tag{34}$$

where *τ*ˆ*id* represents the compound disturbance compensated by FTESO (7), which is supposed to be invariant over a short period. It should be noted that we assume the center of gravity and buoyancy of the *i*th AUV to coincide, which allows *g<sup>i</sup>* (*ηi*) to approximate to zero. We select *xiv*(*k*) = ' *vi*(*k*) *τi*(*k* − 1) (*<sup>T</sup>* as the state variable and take the increment Δ*uiτ*(*k*) = *τi*(*k*) − *τi*(*k* − 1) as the control input. This allows us to reformulate the innerloop predictive model as follows:

$$\mathbf{x}\_{\rm iv}(k+1) = \mathbf{A}\_{\rm iv}\mathbf{x}\_{\rm iv}(k) + \mathbf{B}\_{\rm iv}\Delta\mathbf{u}\_{\rm ir}(k) + \mathbf{D}\_{\rm iv} \tag{35}$$

$$y\_{\rm iv}(k) = \mathcal{C}\_{\rm iv} \mathfrak{x}\_{\rm iv}(k) \tag{36}$$

where *Aiv* = *I*<sup>5</sup> − *M*<sup>∗</sup> *i* <sup>−</sup>1*Ts C*∗ *<sup>i</sup>* + *D*<sup>∗</sup> *i M*∗ *i* <sup>−</sup>1*Ts* **0**<sup>5</sup> *I*<sup>5</sup> <sup>∈</sup> <sup>R</sup>10×10, *<sup>B</sup>iv* <sup>=</sup> *M*∗ *i* <sup>−</sup>1*Ts I*5 <sup>∈</sup> <sup>R</sup>10×5,

*Civ* = ' *I*<sup>5</sup> **0**<sup>5</sup> ( <sup>∈</sup> <sup>R</sup>5×10, and *<sup>D</sup>iv* <sup>=</sup> *M*∗ *i* <sup>−</sup>1*Tsτ*ˆ*id* **0**5×<sup>1</sup> <sup>∈</sup> <sup>R</sup>10. Similar to our previous approach, we can characterize the relationship between the predicted output vector sequence and the

control increment sequence using the following prediction equation:

$$\mathbf{Y}\_{\rm iv} = H\_{\rm ix}^2 \mathbf{x}\_{\rm iv}(k) + H\_{\rm iul}^2 \Delta \mathbf{U}\_{\rm ir} + \overline{\mathbf{D}}\_{\rm iv} \tag{37}$$

where *Yiv* = ' *yiv*(*k* + 1|*k*), *yiv*(*k* + 2|*k*),..., *yiv k* + *Np*<sup>2</sup> *k* (*<sup>T</sup>* <sup>∈</sup> <sup>R</sup>5*Np*<sup>2</sup> , <sup>Δ</sup>*Ui<sup>τ</sup>* <sup>=</sup> [Δ*uiτ*(*k*|*k*), <sup>Δ</sup>*uiτ*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>|*k*),..., <sup>Δ</sup>*uiτ*(*<sup>k</sup>* <sup>+</sup> *Nc*<sup>2</sup> <sup>−</sup> <sup>1</sup>|*k*)]*<sup>T</sup>* <sup>∈</sup> <sup>R</sup>5*Nc*<sup>2</sup> , *H*2 *ix* = *CivAiv*,*CivA*<sup>2</sup> *iv*,...,*CivANp*<sup>2</sup> *iv <sup>T</sup>* <sup>∈</sup> <sup>R</sup>5*Np*2<sup>×</sup>10, *H*2 *iu* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ *CivBiv* **0**<sup>5</sup> ··· **0**<sup>5</sup> *CivAivBiv CivBiv* ··· **0**<sup>5</sup> . . . . . . ... . . . *<sup>C</sup>ivANp*2−<sup>1</sup> *iv <sup>B</sup>iv <sup>C</sup>ivANp*2−<sup>2</sup> *iv <sup>B</sup>iv* ··· *<sup>C</sup>ivANp*2−*Nc*<sup>2</sup> *iv Biv* ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ <sup>∈</sup> <sup>R</sup>5*Np*2×5*Nc*<sup>2</sup> , and *Div* = *CivDiv*,*CivAivDiv* + *CivDiv*,...,*Civ Np*2−1 ∑ *n*=0 *An ivDiv<sup>T</sup>* <sup>∈</sup> <sup>R</sup>5*Np*<sup>2</sup> . *Np*<sup>2</sup> and *Nc*<sup>2</sup> de-

note the prediction and control horizon of the inner-loop controller, respectively.

According to the control objective, we assess the constraints on the control input increment and the actuator saturation in the inner-loop subsystem, as follows:

$$
\Delta \mathfrak{u}\_{i\tau}^{\min} \le \Delta \mathfrak{u}\_{i\tau}(k) \le \Delta \mathfrak{u}\_{i\tau}^{\max} \tag{38}
$$

$$
\mathfrak{r}\_i^{\min} \le \mathfrak{r}\_i(k) \le \mathfrak{r}\_i^{\max} \tag{39}
$$

where *τ*min *<sup>i</sup>* and <sup>Δ</sup>*u*min *iv* represent predefined lower bounds, while *<sup>τ</sup>*max *<sup>i</sup>* and <sup>Δ</sup>*u*max *iv* denote predefined upper bounds. The actuator saturation constraint (39) can be transformed into an input incremental constraint, and we can express the above constraints in a compact linear constraint form:

$$
\Gamma\_{\rm iv} \Delta \mathcal{U}\_{i\tau} \le \gamma\_{\rm iv} \tag{40}
$$

where *Γiv* = ⎡ ⎢ ⎢ ⎣ *I*5*Nc*<sup>2</sup> −*I*5*Nc*<sup>2</sup> *Iv*<sup>1</sup> −*Iv*<sup>1</sup> ⎤ ⎥ ⎥ ⎦, *<sup>γ</sup>i<sup>η</sup>* <sup>=</sup> ⎡ ⎢ ⎢ ⎣ Δ*U*max *iτ* <sup>−</sup>Δ*U*min *iτ τ*max *<sup>i</sup>* − *Iv*2*τi*(*k* − 1) <sup>−</sup>*τ*min *iv* + *Iv*2*τi*(*k* − 1) ⎤ ⎥ ⎥ ⎦ , with *Iv*<sup>1</sup> = ⎡ ⎢ ⎢ ⎢ ⎣ *I*<sup>5</sup> **0**<sup>5</sup> ··· **0**<sup>5</sup> *I*<sup>5</sup> *I*<sup>5</sup> ··· **0**<sup>5</sup> . . . . . . ... . . . *I*<sup>5</sup> *I*<sup>5</sup> ··· *I*<sup>5</sup> ⎤ ⎥ ⎥ ⎥ ⎦ ∈

R5*Nc*2×5*Nc*<sup>2</sup> and *<sup>I</sup>v*<sup>2</sup> = [*I*5,*I*5,...,*I*5] *<sup>T</sup>* <sup>∈</sup> <sup>R</sup>5*Nc*<sup>2</sup> .

To achieve the convergence of the formation tracking velocity to the desired value, we define the local distributed cost function of the inner-loop subsystem as follows:

$$J\_{\rm iv}(k) = \sum\_{l=1}^{N\_{p2}} \left\| \left( \mathfrak{y}\_{\rm iv}(k+l|k) - \mathfrak{v}\_{\rm ir}(k+l) \right) \right\|\_{\mathbf{Q}\_{\rm iv}}^2 + \sum\_{l=0}^{N\_{\rm l2}-1} \left\| \Delta \mathfrak{u}\_{\rm ir}(k+l|k) \right\|\_{\mathbf{R}\_{\rm i2}}^2 \tag{41}$$

where *yiv*(*k* + *l*|*k*) and Δ*uiτ*(*k* + *j*|*k*) denote the predicted value of *yiv*(*k* + *l*) and Δ*uiτ*(*k* + *j*) at time *k*, respectively. *Qiv* and *Ri*<sup>2</sup> are given weight matrices.

By substituting (37) into (41), we can formulate the DMPC optimization problem for the inner-loop subsystem of the *i*th AUV at sampling instant *k* as the following QP form:

$$\begin{aligned} \Delta \mathbf{U}\_{i\tau}^{\*} &= \underset{\Delta \mathbf{U}\_{i\tau}}{\operatorname{argmin}} \left( \frac{1}{2} \Delta \mathbf{U}\_{i\tau}^{T} \mathbf{W}\_{iv} \Delta \mathbf{U}\_{i\tau} + f\_{iv}^{T} \Delta \mathbf{U}\_{i\tau} \right) \\ \text{s.t. } & \Gamma\_{iv} \Delta \mathbf{U}\_{i\tau} \le \gamma\_{iv} \end{aligned} \tag{42}$$

where *Wiv* = *Ri*<sup>2</sup> + *H*2*<sup>T</sup> iu <sup>Q</sup>ivH*<sup>2</sup> *iu* and *<sup>f</sup>iv* = *<sup>H</sup>*2*<sup>T</sup> iu Qiv H*2 *ixxiv* + *Div* − *Vir* , with *vir* = ' *vir*(*k* + 1),..., *vir k* + *Np*<sup>2</sup> (*<sup>T</sup>* <sup>∈</sup> <sup>R</sup>5*Np*<sup>2</sup> , *<sup>Q</sup>iv* <sup>=</sup> *diag*{*Qiv*, *<sup>Q</sup>iv*,..., *<sup>Q</sup>iv*} <sup>∈</sup> <sup>R</sup>5*Np*2×5*Np*<sup>2</sup> , and *<sup>R</sup>i*<sup>2</sup> <sup>=</sup> *diag*{*Ri*2, *<sup>R</sup>i*2,..., *<sup>R</sup>i*2} <sup>∈</sup> <sup>R</sup>5*Nc*2×5*Nc*<sup>2</sup> .

The solution of the QP optimization problem (42) yields the optimal control input increment sequence Δ*U*∗ *<sup>i</sup><sup>τ</sup>* at time *k*. However, only the first element Δ*u*<sup>∗</sup> *<sup>i</sup>τ*(*k*|*k*) of the sequence is used for the *i*th AUV to obtain the optimal control force and moment *τ*∗ *<sup>i</sup>* (*k*) = *τi*(*k* − 1) + Δ*u*<sup>∗</sup> *<sup>i</sup>τ*(*k*). The Δ*u*<sup>∗</sup> *<sup>i</sup>τ*(*k*) is recalculated at each sampling instant, the *i*th AUV repeatedly calculates and executes *τ*∗ *<sup>i</sup>* (*k*) to achieve receding optimization. The predicted state *xiv*(*k* + 1) and the optimal input *τ*<sup>∗</sup> *<sup>i</sup>* (*k*) are both determined solely by the current state *xiv*(*k*).

With the parallel optimization of *N* AUV subsystems, all local optimization problems are solved simultaneously at each sampling moment. One or more information interactions occur between local controllers to obtain the optimal input sequence for that moment. Thus, the proposed control law can compensate well for the compound disturbances, which consist of model uncertainties and external disturbances. This occurs throughout the iterative optimization process, while simultaneously ensuring collision avoidance and formation tracking control tasks under complex constraints.

#### *3.4. Use of Laguerre Functions in the DMPC Design*

This subsection introduces a strategy to handle the computational burden caused by a longer control horizon and dual closed-loop structure. This is the main difficulty in our theoretical analysis. The Laguerre orthogonal functions are leveraged in the DMPC design to decrease the order of the input incremental matrices. This approach permits a reduction in input variables during each control cycle, thereby reducing the computational burden in the time interval and improving real-time performance.

The Laguerre functions are a set of discrete orthogonal polynomial functions, let it be *l*1(*k*), *l*2(*k*),..., *lM*(*k*), the z-transfer of the *m*th Laguerre function is expressed as follows:

$$\chi\_{\mathfrak{m}}(z) = \frac{\sqrt{1-a^2}}{z-a} \left[ \frac{1-az}{z-a} \right]^{m-1} \tag{43}$$

where 0 ≤ *a* < 1 denotes the pole of the Laguerre function, also known as the scaling factor. It can be verified that X*m* satisfies the following orthogonality:

$$\begin{cases} \frac{1}{2\pi} \int\_{-\pi}^{\pi} \mathsf{X}\_{\pi} \left( e^{j\omega} \right) \mathsf{X}\_{\pi} \left( e^{j\omega} \right)^{\*} d\omega = 1 & m = n \\\\ \frac{1}{2\pi} \int\_{-\pi}^{\pi} \mathsf{X}\_{\pi} \left( e^{j\omega} \right) \mathsf{X}\_{\pi} \left( e^{j\omega} \right)^{\*} d\omega = 0 & m \neq n \end{cases} \tag{44}$$

where (·) <sup>∗</sup> denotes complex conjugate of (·).

The discrete Laguerre functions are defined by taking the inverse Z-transform of (43), i.e., *lm*(*k*) <sup>=</sup> *<sup>Z</sup>*−1{X*m*(*z*)}. Given the network structure of <sup>X</sup>*m*(*z*) and the recurrence relation, the set of discrete Laguerre functions satisfies the following difference equation:

$$L(k+1) = \Xi L(k) \tag{45}$$

$$\begin{array}{ccccc} \text{where} & \Xi & = & \begin{bmatrix} a & 0 & 0 & \cdots & 0 \\ \beta & a & 0 & \cdots & 0 \\ -a\beta & \beta & a & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ (-a)^{M-2}\beta & (-a)^{M-3}\beta & (-a)^{M-4}\beta & \cdots & a \end{bmatrix} \end{array} \quad \text{and} \quad L(k) &= \begin{array}{ccccc} \\ & \end{array}$$

[*l*1(*k*), *<sup>l</sup>*2(*k*),..., *lM*(*k*)]*T*, with *<sup>β</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>a</sup>*<sup>2</sup> and initial condition *<sup>L</sup>*(0) <sup>=</sup> /*β* 1, −*a*, *<sup>a</sup>*2, −*a*3,...,(−*a*) *M*−1 *T* . Note that at *a* = 0, the Laguerre functions are converted to impulse functions.

Assuming the current moment is *k*, the input increment of the single-input system at the next time *l*, represented by the Laguerre function, is:

$$
\Delta\mu(k+l) = \sum\_{m=1}^{M} \kappa\_m l\_m(l) = L(l)^T K \tag{46}
$$

where *K* = [*κ*1, *κ*2,..., *κM*] *<sup>T</sup>*. When we extend this to the multi-AUV system, each AUV has five independent control inputs, and the input increment of the *i*th AUV is as follows:

$$\left(\Delta\boldsymbol{\mu}\_{i}(\boldsymbol{k})\right)^{T} = \left[\boldsymbol{L}\_{i}^{1}(\boldsymbol{k})^{T}\boldsymbol{\mathcal{K}}\_{i}^{1}, \boldsymbol{L}\_{i}^{2}(\boldsymbol{k})^{T}\boldsymbol{\mathcal{K}}\_{i}^{2}, \dots, \boldsymbol{L}\_{i}^{5}(\boldsymbol{k})^{T}\boldsymbol{\mathcal{K}}\_{i}^{5}\right] = \boldsymbol{\mathsf{T}}\_{i}(\boldsymbol{k})^{T}\boldsymbol{\mathsf{X}}\_{i} \tag{47}$$

where *L<sup>p</sup> <sup>i</sup>* (*k*) = *l p <sup>i</sup>*1(*k*), *l p <sup>i</sup>*2(*k*),..., *l p iM*(*k*) *T* and *K<sup>p</sup> <sup>i</sup>* = *κp <sup>i</sup>*1, *<sup>κ</sup><sup>p</sup> <sup>i</sup>*2,..., *<sup>κ</sup><sup>p</sup> iM <sup>T</sup>* , with *p* = 1, 2, ... 5. *<sup>L</sup>i*(*k*) <sup>=</sup> *diag*& *L*1 *<sup>i</sup>* (*k*) *<sup>T</sup>*, *L*<sup>2</sup> *<sup>i</sup>* (*k*) *<sup>T</sup>*,..., *L*<sup>5</sup> *<sup>i</sup>* (*k*) *T* ! , and *K<sup>i</sup>* = ' *K*1*<sup>T</sup> <sup>i</sup>* , *<sup>K</sup>*2*<sup>T</sup> <sup>i</sup>* ,..., *<sup>K</sup>*5*<sup>T</sup> i* (*T* . Note that within a multi-input structure, the scaling factor *ap* and the number of polynomial terms *Mp* can be selected independently for each input signal.

For illustrative purposes, the inner-loop predictive controller of the *i*th AUV is taken as an example. If we partition the input matrix into *Biv* = ' *B*1 *iv <sup>B</sup>*<sup>2</sup> *iv* ... *<sup>B</sup>*<sup>5</sup> *iv* ( , the prediction of the system output in the next *l* steps can be derived as follows:

$$\begin{split} y\_{i\bar{\nu}}(k+l|k) &= \sum\_{j=0}^{l-1} \mathbb{C}\_{i\bar{\nu}} \mathbf{A}\_{i\bar{\nu}}^{l-j-1} \left[ \mathbf{B}\_{i\bar{\nu}}^{1} \mathbf{L}\_{i}^{1}(j)^{\mathrm{T}} \mathbf{K}\_{i\bar{\nu}}^{1} \cdot \mathbf{B}\_{i\bar{\nu}}^{2} \mathbf{L}\_{i}^{2}(j)^{\mathrm{T}} \mathbf{K}\_{i\bar{\nu}}^{2} \quad \dots \ \mathbf{B}\_{i\bar{\nu}}^{5} \mathbf{L}\_{i}^{5}(j)^{\mathrm{T}} \mathbf{K}\_{i\bar{\nu}}^{5} \right] \\ &+ \mathbb{C}\_{i\bar{\nu}} \mathbf{A}\_{i\bar{\nu}}^{l} \mathbf{x}\_{i\bar{\nu}}(k) + \sum\_{j=0}^{l-1} \mathbb{C}\_{i\bar{\nu}} \mathbf{A}\_{i\bar{\nu}}^{l-j-1} \mathbf{D}\_{i\bar{\nu}}. \end{split} \tag{48}$$

For a compact notation, we denote (48) by the following:

$$y\_{iv}(k+l|k) = \mathcal{C}\_{iv} \mathcal{A}\_{iv}^{l} \mathfrak{x}\_{iv}(k) + \mu\_{i}(l)^{T} \,\, \overline{\mathcal{K}}\_{i\tau} + \mathcal{D}\_{iv}^{l} \tag{49}$$

where *μi*(*l*) *<sup>T</sup>* <sup>=</sup> *<sup>l</sup>*−<sup>1</sup> ∑ *j*=0 *<sup>C</sup>ivAl*−*j*−<sup>1</sup> *iv B*1 *ivL*<sup>1</sup> *<sup>i</sup>* (*j*) *<sup>T</sup> B*<sup>2</sup> *ivL*<sup>2</sup> *<sup>i</sup>* (*j*) *<sup>T</sup>* ... *B*<sup>5</sup> *ivL*<sup>5</sup> *<sup>i</sup>* (*j*) *T* and *D<sup>l</sup> iv* =

*l*−1 ∑ *j*=0 *<sup>C</sup>ivAl*−*j*−<sup>1</sup> *iv Div*. *Ki<sup>τ</sup>* as the parameter vector that is to be optimized.

First, we employ the Laguerre function to optimize the constraint terms (38) and (39), leading to the following constraint form:

$$
\Delta \mu\_{i\tau}^{\min} \le \mathbf{Z}\_{i\tau}^{\top} \mathbf{Z}\_{i\tau} \le \Delta \mu\_{i\tau}^{\max} \tag{50}
$$

.

$$
\tau\_i^{\min} \le \widehat{L}\_{i\tau} \mathbb{Z}\_{i\tau} + \tau\_i(k-1) \le \tau\_i^{\max} \tag{51}
$$

$$\text{where}\\
\overline{L}\_{i\tau} = \text{diag}\left\{ L^1\_{i\tau}(l)^T, \dots, L^5\_{i\tau}(l)^T \right\} \\
\text{and}\\
\widehat{L}\_{i\tau} = \text{diag}\left\{ \sum\_{j=0}^{l-1} L^1\_{i\tau}(j)^T, \sum\_{j=0}^{l-1} L^2\_{i\tau}(j)^T, \dots, \sum\_{j=0}^{l-1} L^5\_{i\tau}(j)^T \right\}.$$

Given that the Laguerre functions are orthonormal for a sufficiently large control horizon *Nc*2. Substituting (47) into (41) and using the orthogonality (44) of the Laguerre function (i.e., the inner product of different terms is 0 and the same term is 1), the following derivation can be performed to obtain the reconstructed form of the cost function (41):

*Jiv*(*k*) = *Np*<sup>2</sup> ∑ *l*=1 (*yiv*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*) <sup>−</sup> *<sup>v</sup>ir*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*))<sup>2</sup> *<sup>Q</sup>iv* <sup>+</sup> *Nc*2−<sup>1</sup> ∑ *l*=0 Δ*uiτ*(*k* + *l*|*k*) *<sup>T</sup>Ri*2Δ*uiτ*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*) = *Np*<sup>2</sup> ∑ *l*=1 (*yiv*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*) <sup>−</sup> *<sup>v</sup>ir*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*))<sup>2</sup> *<sup>Q</sup>iv* <sup>+</sup> *Nc*2−<sup>1</sup> ∑ *l*=0 *Liτ*(*l*) *<sup>T</sup>Ki<sup>τ</sup> Ri*<sup>2</sup> *Liτ*(*l*) *<sup>T</sup>Ki<sup>τ</sup> T* = *Np*<sup>2</sup> ∑ *l*=1 [*yiv*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*) <sup>−</sup> *<sup>v</sup>ir*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*)]*TQiv*[*yiv*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*) <sup>−</sup> *<sup>v</sup>ir*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*)] + *Nc*2−1 ∑ *l*=0 *diag L*1 *<sup>i</sup>τ*(*l*), *<sup>L</sup>*<sup>2</sup> *<sup>i</sup>τ*(*l*),..., *<sup>L</sup>*<sup>5</sup> *<sup>i</sup>τ*(*l*) *Ki<sup>τ</sup> Ri*<sup>2</sup> *diag L*1 *<sup>i</sup>τ*(*l*), *<sup>L</sup>*<sup>2</sup> *<sup>i</sup>τ*(*l*),..., *<sup>L</sup>*<sup>5</sup> *<sup>i</sup>τ*(*l*) *Ki<sup>τ</sup> T* = *Np*<sup>2</sup> ∑ *l*=1 [*yiv*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*) <sup>−</sup> *<sup>v</sup>ir*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*)]*TQiv*[*yiv*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*) <sup>−</sup> *<sup>v</sup>ir*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*)] <sup>+</sup> *<sup>K</sup><sup>T</sup> <sup>i</sup>τRi*2*Kiτ*. (52)

By substituting (49) into (52), we can rewrite the DMPC optimization problem (42) for the inner-loop subsystem of the *i*th AUV as:

$$\begin{array}{ll}\underset{\begin{subarray}{c}\mathbf{K}^{\circ}\_{ir}\\\mathbf{K}^{\circ}\_{ir}\end{subarray}}{\min}J\_{iv}(k) = \underset{\begin{subarray}{c}\mathbf{K}^{\circ}\_{ir}\\\mathbf{K}^{\circ}\_{ir}\end{subarray}}{\min} \left(\ \frac{1}{2}\mathbf{K}^{T}\_{ir}\mathbf{W}\_{iL}\mathbf{Z}\_{i\tau} + f^{T}\_{iL}\mathbf{Z}\_{i\tau}\right) \\ s.t. \ (50), \end{array} \tag{53}$$

where *<sup>W</sup>iL* <sup>=</sup> <sup>∑</sup>*Np*<sup>2</sup> *<sup>l</sup>*=<sup>1</sup> *μi*(*l*)*Qivμi*(*l*) *<sup>T</sup>* <sup>+</sup> *<sup>R</sup>i*<sup>2</sup> and *<sup>f</sup>iL* <sup>=</sup> <sup>∑</sup>*Np*<sup>2</sup> *<sup>l</sup>*=<sup>1</sup> *<sup>μ</sup>i*(*l*)*Qiv CivA<sup>l</sup> ivxiv*(*k*) + *<sup>D</sup><sup>l</sup> iv* − *vir*(*k* + *l*) .

The QP optimization Equation (53), with constraints, can be solved to obtain the optimal parameter vector *K*<sup>∗</sup> *iτ*. This vector replaces the conventional DMPC method calculation of Δ*u*∗ *<sup>i</sup>τ*. Thus, the optimal input increment of the inner-loop subsystem is indirectly obtained by the rolling optimized control law, Δ*uiτ*(*k*) *<sup>T</sup>* = *<sup>L</sup>iτ*(0) *<sup>T</sup>Kiτ*, until the control variables at the next moment are calculated. This iterative process ensures the achievement of receding horizon optimization. The use of the Laguerre function in the design of the outer-loop predictive controller is not included here, as its analysis parallels that of the inner-loop controller described above.

**Remark 2.** By parameterizing the input increment sequence using the Laguerre function, the input matrix order in the prediction horizon can be lowered, thereby reducing the computational load online. This property enables its application in large-scale and realtime AUV control systems. With the employment of the Laguerre function, the coefficients *ap* and *Mp* can also be served as tuning parameters, in addition to the control and prediction horizon and weighting matrices. Larger *ap* and *Mp* lead to faster closed-loop responses [38].

#### *3.5. Stability Analysis*

A notable attribute of the MPC is the potential for establishing the stability of a closedloop system under certain conditions. Extending this to cases using Laguerre polynomials, a terminal state constraint is utilized to analyze the stability of the closed-loop system. Specifically, for the inner-loop subsystem, an additional constraint is attached to the final state of the receding optimization problem: *xiv k* + *Np*<sup>2</sup> = 0, where *xiv k* + *Np*<sup>2</sup> is the terminal state produced under the effect of the control sequence, Δ*uiτ*(*k* + *l*) *<sup>T</sup>* = *<sup>L</sup>iτ*(*l*) *<sup>T</sup>Kiτ*.

**Theorem 2.** *Consider the inner-loop subsystem (35) and (36) of the ith AUV in the formation control system, which has a local cost function (41) with constraints (38) and (39). The inner-loop predictive control subsystem is asymptotically stable if for each sampling instant k, there exists a solutionKi<sup>τ</sup> such that the performance index Jiv is minimized subject to the terminal state constraint*. **Proof of Theorem 2.** Constructing an appropriate Lyapunov function is key to ensuring the stability of the DMPC system. Select the cost function *Jiv*(*k*) as the Lyapunov function *Vi*2(*y*(*k*), *k*):

$$V\_{l2}(\boldsymbol{y}(k),k) = \sum\_{l=1}^{N\_{\mathcal{P}2}} \tilde{\boldsymbol{y}}\_{iv}(k+l|k)^{T} \mathbf{Q}\_{iv} \tilde{\boldsymbol{y}}\_{iv}(k+l|k) + \sum\_{l=0}^{N\_{\mathcal{Q}2}-1} \Delta \boldsymbol{u}\_{i\tau}(k+l)^{T} \mathbf{R}\_{l2} \Delta \boldsymbol{u}\_{i\tau}(k+l) \tag{54}$$

where *<sup>y</sup>*˜*iv*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*) <sup>=</sup> *<sup>y</sup>iv*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*) <sup>−</sup>*vir*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*) and *<sup>y</sup>iv*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup>*|*k*) <sup>=</sup> *<sup>C</sup>ivA<sup>l</sup> ivxiv*(*k*) + *μi*(*l*) *<sup>T</sup> K<sup>k</sup> <sup>i</sup><sup>τ</sup>* + *Dl iv*, *<sup>K</sup><sup>k</sup> <sup>i</sup><sup>τ</sup>* is the parameter vector solution of the cost function (41) under the original and terminal constraints at moment *k*, and input increment Δ*uiτ*(*k* + *l*) *<sup>T</sup>* = *<sup>L</sup>iτ*(*l*) *<sup>T</sup>K<sup>k</sup> <sup>i</sup>τ*. It is clear that *Vi*2(*y*(*k*), *k*) is positive definite and tends to infinity as *yiv*(*k*) tends to infinity. Similarly, the Lyapunov function at moment *k* + 1 can be derived as:

$$\begin{split} V\_{l2}(y(k+1),k+1) &= \sum\_{l=1}^{N\_{p2}} \bar{\mathbf{y}}\_{\dot{v}\dot{v}} (k+1+l|k+1)^T \mathbf{Q}\_{\dot{v}i} \bar{\mathbf{y}}\_{\dot{v}\dot{v}} (k+1+l|k+1) \\ &+ \sum\_{l=0}^{N\_{l2}-1} \Delta \mathbf{u}\_{l\dot{\tau}} (k+1+l)^T \mathbf{R}\_{l2} \Delta \mathbf{u}\_{l\dot{\tau}} (k+1+l) \end{split} \tag{55}$$

where *<sup>y</sup>iv*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup> <sup>+</sup> *<sup>l</sup>*|*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) <sup>=</sup> *<sup>C</sup>ivA<sup>l</sup> ivxiv*(*k* + 1) + *μi*(*l*) *<sup>T</sup> <sup>K</sup>k*+<sup>1</sup> *<sup>i</sup><sup>τ</sup>* + *D<sup>l</sup> iv*, *<sup>K</sup>k*+<sup>1</sup> *<sup>i</sup><sup>τ</sup>* is the parameter vector solution at time *k* + 1, and Δ*uiτ*(*k* + 1 + *l*) *<sup>T</sup>* = *<sup>L</sup>iτ*(*l*) *<sup>T</sup>Kk*+<sup>1</sup> *<sup>i</sup><sup>τ</sup>* . Given that *yiv*(*k* + 1) is the response one step ahead of *yiv*(*k*) and *yiv*(*k* + 1) = *CivAivxiv*(*k*) + *CivBiv*Δ*uiτ*(*k*) + *<sup>C</sup>ivDiv*, the feasible solution of *<sup>K</sup>k*+<sup>1</sup> *<sup>i</sup><sup>τ</sup>* corresponding to the initial output *yiv*(*k* + 1) in the receding horizon is *K<sup>k</sup> <sup>i</sup>τ*. Therefore, the feasible solution sequence at moment *k* + 1 is to move the elements in *Liτ*(0) *<sup>T</sup>K<sup>k</sup> <sup>i</sup>τ*, *Liτ*(1) *<sup>T</sup>K<sup>k</sup> <sup>i</sup>τ*, ..., *Liτ*(*Nc*<sup>2</sup> − 1) *<sup>T</sup>K<sup>k</sup> <sup>i</sup><sup>τ</sup>* one step forward and substitute the last element with *0*, i.e., *Liτ*(1) *<sup>T</sup>K<sup>k</sup> <sup>i</sup>τ*, *Liτ*(2) *<sup>T</sup>K<sup>k</sup> <sup>i</sup>τ*, ..., *Liτ*(*Nc*<sup>2</sup> − 1) *<sup>T</sup>K<sup>k</sup> <sup>i</sup>τ*, *0*5×1. Due to the optimality of the solution *<sup>K</sup>k*+<sup>1</sup> *<sup>i</sup><sup>τ</sup>* at *k* + 1, it follows that

$$\mathcal{V}\_{i2}(y(k+1), k+1) \le \mathcal{V}\_{i2}(y(k+1), k+1) \tag{56}$$

where *V*ˆ *<sup>i</sup>*2(*y*(*<sup>k</sup>* <sup>+</sup> <sup>1</sup>), *<sup>k</sup>* <sup>+</sup> <sup>1</sup>) is identical to (55) except that the parameter vector solution *<sup>K</sup>k*+<sup>1</sup> *iτ* in the control sequence is replaced by the feasible solution *K<sup>k</sup> <sup>i</sup>τ*. The difference between *Vi*2(*y*(*k*), *k*) and *Vi*2(*y*(*k* + 1), *k* + 1) is then bounded by the following:

$$V\_{i2}(\mathfrak{y}(k+1),k+1) - V\_{i2}(\mathfrak{y}(k),k) \le \mathring{V}\_{i2}(\mathfrak{y}(k+1),k+1) - V\_{i2}(\mathfrak{y}(k),k). \tag{57}$$

Eliminate the same terms in the control sequence and output sequence of *V*ˆ *<sup>i</sup>*2(*y*(*k* + 1), *k* + 1) and *Vi*2(*y*(*k*), *k*) at moments *k* + 1, *k* + 2,..., *k* + *Np*<sup>2</sup> − 1, and we can derive the following equation:

$$\begin{split} \dot{\mathcal{W}}\_{\rm l2}(\boldsymbol{y}(k+1),k+1) - \boldsymbol{V}\_{\rm l2}(\boldsymbol{y}(k),k) &= \quad \ddot{\boldsymbol{y}}\_{\rm l\bar{v}}(k+N\_{p2}|\boldsymbol{k})^{\mathsf{T}} \mathbf{Q}\_{\boldsymbol{v}\bar{v}} \ddot{\boldsymbol{y}}\_{\rm l\bar{v}}(k+N\_{p2}|\boldsymbol{k}) \\ &- \ddot{\boldsymbol{y}}\_{\rm l\bar{v}}(k+1|\boldsymbol{k})^{\mathsf{T}} \mathbf{Q}\_{\boldsymbol{v}\bar{v}} \ddot{\boldsymbol{y}}\_{\rm l\bar{v}}(k+1|\boldsymbol{k}) - \boldsymbol{\Lambda} \boldsymbol{u}\_{\rm l\bar{v}}(k)^{\mathsf{T}} \mathbf{R}\_{\bar{v}2} \boldsymbol{\Lambda} \boldsymbol{u}\_{\rm l\bar{v}}(k). \end{split} \tag{58}$$

Given the terminal constraint *xiv k* + *Np*<sup>2</sup> = 0 is applied, equivalent to *yiv k* + *Np*<sup>2</sup> = 0, we have the following:

$$\begin{split} \boldsymbol{\mathring{W}\_{l2}}(\boldsymbol{y}(k+1),k+1) - \boldsymbol{V\_{l2}}(\boldsymbol{y}(k),k) &= -\boldsymbol{\boldsymbol{v}\_{\text{ir}}}(k+N\_{p2})^{\mathsf{T}} \mathbf{Q}\_{\text{ir}} \boldsymbol{\boldsymbol{v}}\_{\text{ir}}(k+N\_{p2}) \\ &- \boldsymbol{\breve{y}}\_{\text{iv}}(k+1|k)^{\mathsf{T}} \mathbf{Q}\_{\text{iv}} \boldsymbol{\breve{y}}\_{\text{iv}}(k+1|k) - \boldsymbol{\Lambda} \boldsymbol{\boldsymbol{u}}\_{\text{i}\text{\text{\texttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttexttext$$

This allows inequality (57) to be converted into:

$$\begin{split} V\_{l2}(\boldsymbol{y}(k+1),k+1) - V\_{l2}(\boldsymbol{y}(k),k) &\leq \ -\boldsymbol{v}\_{\text{ir}}(k+N\_{p2})^{T}\boldsymbol{Q}\_{\text{ir}}\boldsymbol{v}\_{\text{ir}}(k+N\_{p2}) \\ &- \boldsymbol{\bar{y}}\_{\text{ir}}(k+1|k)^{T}\boldsymbol{\bar{Q}}\_{\text{ir}}\boldsymbol{\bar{y}}\_{\text{ir}}(k+1|k) - \boldsymbol{\Delta}\boldsymbol{u}\_{\text{ir}}(k)^{T}\boldsymbol{R}\_{l2}\boldsymbol{\Delta}\boldsymbol{u}\_{\text{ir}}(k) < 0. \end{split} \tag{60}$$

Namely, *Vi*2(*y*(*k* + 1), *k* + 1) < *Vi*2(*y*(*k*), *k*); the Lyapunov function is monotonically decreasing. This proves that the inner loop subsystem is asymptotically stable. -

Next, we analyze the stability of the entire closed-loop system. Analogous to the proof of Theorem 2, we select *Jiη*(*k*) as the Lyapunov function *Vi*3(*y*(*k*), *k*) of the outerloop subsystem:

$$\begin{split} V\_{i3}(\boldsymbol{y}(k),k) &= \sum\_{l=1}^{N\_{p1}} \left\| \left( \boldsymbol{y}\_{i\eta}(k+l|k) - \boldsymbol{y}\_{if}(k+l) \right) \right\|\_{\mathbf{Q}\_{if}}^{2} + \sum\_{l=0}^{N\_{l1}-1} \left\| \Delta \boldsymbol{u}\_{i\eta}(k+l|k) \right\|\_{\mathbf{R}\_{il}}^{2} \\ &+ \sum\_{l=1}^{N\_{p1}} \sum\_{j \in \mathcal{N}\_{l}} a\_{ij} \left\| \left( \boldsymbol{y}\_{i\eta}(k+l|k) - \boldsymbol{y}\_{ij}(k+l) \right) \right\|\_{\mathbf{Q}\_{ij}}^{2}. \end{split} \tag{61}$$

According to the idea of the proof of Theorem 2, the following inequality can be obtained:

$$\begin{split} V\_{l\overline{l}}(y(k+1),k+1) - V\_{l\overline{l}}(y(k),k) \leq & -y\_{if}(k+N\_{p1})^{\mathsf{T}} \mathbf{Q}\_{if} \mathbf{y}\_{if}(k+N\_{p1}) - \mathbf{\bar{y}}\_{i\eta}(k+1|k)^{\mathsf{T}} \mathbf{Q}\_{if} \mathbf{\bar{y}}\_{i\eta}(k+1|k) \\ & -N(N-1)\mathbf{r}\_{l\overline{l}}(k+N\_{p1})^{\mathsf{T}} \mathbf{Q}\_{if} \mathbf{r}\_{l\overline{l}}(k+N\_{p1}) - \boldsymbol{\Delta}\boldsymbol{u}\_{l\overline{v}}(k)^{\mathsf{T}} \mathbf{R}\_{l1} \boldsymbol{\Delta}\boldsymbol{u}\_{l\overline{v}}(k) < 0. \end{split} \tag{62}$$

Next, we set the Lyapunov function of the entire closed-loop system as follows:

$$V\_{i4}(\boldsymbol{y}(k),k) = V\_{i2}(\boldsymbol{y}(k),k) + V\_{i3}(\boldsymbol{y}(k),k). \tag{63}$$

From inequalities (60) and (62), we have the following:

$$\begin{aligned} V\_{i4}(y(k+1), k+1) - V\_{i4}(y(k), k) &= V\_{i2}(y(k+1), k+1) - V\_{i2}(y(k), k) \\ &+ V\_{i3}(y(k+1), k+1) - V\_{i3}(y(k), k) < 0 \end{aligned} \tag{64}$$

As a result, the entire closed-loop system is asymptotically stable.

#### **4. Simulation**

In this section, some simulation analyses are conducted to verify the effectiveness and robustness of the proposed control scheme. A formation system consisting of four AUVs (*N* = 4, *i* = 1, 2, 3, 4) with a virtual leader (AUV0) is considered. The directed communication topology for the simulation is depicted in Figure 3, the meaning of the arrows is the direction of the communication or information flow between the nodes in the formation network. Initial values for *xi*, *yi*, and *zi* are randomly distributed within the intervals [10, 40], [0, 30], and [−10, 0], respectively, while the attitude angles *θ<sup>i</sup>* and *ψ<sup>i</sup>* lie within the intervals [−*π*/18, *π*/18] and [0, *π*], respectively. The parameters related to the AUVs are based on previous research [39]. A diamond formation was predefined to facilitate omnidirectional exploration, with the desired formation configuration preset to *r*<sup>1</sup> *<sup>f</sup>* = [0, 0, 6.5, 0, 0] *<sup>T</sup>*, *<sup>r</sup>*<sup>2</sup> *<sup>f</sup>* <sup>=</sup> [0, <sup>−</sup>7.5, 0, 0, 0] *<sup>T</sup>*, *<sup>r</sup>*<sup>3</sup> *<sup>f</sup>* <sup>=</sup> [0, 0, <sup>−</sup>6.5, 0, 0] *<sup>T</sup>*, and *r*<sup>4</sup> *<sup>f</sup>* = [0, 7.5, 0, 0, 0] *<sup>T</sup>*. *<sup>r</sup>*<sup>12</sup> <sup>=</sup> <sup>−</sup>*r*<sup>21</sup> <sup>=</sup> [0, 7.5, 6.5, 0, 0] *<sup>T</sup>*, *<sup>r</sup>*<sup>13</sup> <sup>=</sup> <sup>−</sup>*r*<sup>31</sup> <sup>=</sup> [0, 0, 13, 0, 0] *<sup>T</sup>*, *r*<sup>14</sup> = −*r*<sup>41</sup> = [0, −7.5, 6.5, 0, 0] *<sup>T</sup>*, *<sup>r</sup>*<sup>23</sup> <sup>=</sup> <sup>−</sup>*r*<sup>32</sup> <sup>=</sup> [0, <sup>−</sup>7.5, 6.5, 0, 0] *<sup>T</sup>*, *<sup>r</sup>*<sup>24</sup> <sup>=</sup> <sup>−</sup>*r*<sup>42</sup> <sup>=</sup> [0, <sup>−</sup>15, 0, 0, 0] *T* and *r*<sup>34</sup> = −*r*<sup>43</sup> = [0, −7.5, −6.5, 0, 0] *<sup>T</sup>*. The safety distance during the formation construction stage is set as *rs* = 3 *m*, while the detection distance measured using sonar is set as *rd* = 6 *m*. To reflect model uncertainties, 20% of the nominal values were taken as model errors, meaning that the parameters for the AUVs in the simulation represented only 80% of the nominal system dynamics. External disturbances were applied to each AUV to evaluate the formation robustness, modeled as follows [34]:

$$\begin{cases} \tau\_{icw} = 0.1 \text{sign}(u\_i) + 0.2 \sin(0.1t) \text{ N} \\ \tau\_{icv} = 0.1 \text{sign}(v\_i) + 0.3 \sin(0.3t) \text{ N} \\ \tau\_{icw} = 0.08 \text{sign}(w\_i) + 0.2 \sin(0.5t) \text{ N} \\ \tau\_{icq} = 0.02 \text{sign}(q\_i) + 0.1 \sin(0.3t) \text{ N} \cdot m \\ \tau\_{icr} = 0.05 \text{sign}(r\_i) + 0.1 \sin(0.3t) \text{ N} \cdot m \end{cases} \tag{65}$$

**Figure 3.** Structure of communication topology.

Each control parameter has its settings guidelines: Given the low driving speed of the AUV in this paper, smaller *Np*<sup>1</sup> and *Np*<sup>2</sup> are intended to be used. During debugging, reduce it if the rapidity is not enough, and increase it if the stability is not good; the selection of *Nc*<sup>1</sup> and *Nc*<sup>2</sup> is based on a trade-off between performance and computation [40]; since we value the position tracking performance more than the velocity tracking performance, *Qi f* is set slightly larger than *Qiv*; to weaken the interaction of angles between AUVs, the orientation weight in *Qij* is set slightly smaller; when tuning *Ri*<sup>1</sup> and *Ri*2, it can be set very small first, and then increase it slightly if the system is stable and the control variable does not change too drastically [41]. By solving the Lyapunov Equation (11), the relationship between the observer gains *βik* and *αi*1, such that *Ai*<sup>1</sup> and *Ai*<sup>2</sup> are Hurwitz matrices, can be obtained, and tuned to select the appropriate values [42]; the Laguerre parameter *ap* is adjusted within the constraint interval, and a smaller *Mp* is selected to coordinate the number of constraints in the optimization problem, and to make a trade-off between response speed and control complexity [38]. Following the above guidelines, we dealt with the main difficulties in the simulation and selected the parameters that produced the optimal simulation results and listed them in Table 1.


**Table 1.** Control parameters of the proposed scheme.

Moreover, based on the actual speed limit of the thruster, we provide the state and input constraints as follows: Δ*U*max *iv* <sup>=</sup> <sup>−</sup>Δ*U*min *iv* = [0.2, 0.1, 0.2, 0.05, 0.05] *T*, *U*max *iv* <sup>=</sup> <sup>−</sup>*U*min *iv* = [1.5, 1, 1, 0.05, 0.2] *<sup>T</sup>*, and Δ*U*max *<sup>i</sup><sup>τ</sup>* <sup>=</sup> <sup>−</sup>Δ*U*min *<sup>i</sup><sup>τ</sup>* = [50, 50, 100, 5, 5] *<sup>T</sup>*. To avoid actuator saturation for each AUV, the bounds of force and moment are set as *τ*min *<sup>i</sup>* = [−200, −500, −500, −7, −10] *<sup>T</sup>* and *τ*max *<sup>i</sup>* = [300, 500, 500, 7, 10] *<sup>T</sup>*.The reference trajectory generated by the virtual leader is a 3-D spiral curve, defined as follows:

$$\begin{cases} \begin{aligned} \chi\_{\mathbb{T}}(t) &= 30 \cos(0.005 \pi t) \\ \mathcal{y}\_{\mathbb{T}}(t) &= 30 \sin(0.005 \pi t) \\ z\_{\mathbb{T}}(t) &= -0.05t - 3 \end{aligned} \end{cases} \tag{66}$$

To verify the disturbance compensation performance of the proposed FTESO (7), we conducted comparative simulations with the ESO (67) from [43] and the FTESO (68) from [18]. Figure 4 shows the norms of the compound disturbance estimation errors *ei*2 = *τ*ˆ*id* − *τid* for the four AUVs under the three observers, characterizing insights

into transient and steady-state responses. We calculated the settling time of the designed FTESO in the simulation and highlighted it on the plots. It is clear from Figure 4 that our proposed third-order fast FTESO can achieve finite-time stabilization, with the estimation errors converging to a small neighborhood of the origin. And the dynamic convergence speed and estimation accuracy of the proposed FTESO are better than ESO (67) and FTESO (68) with less chattering. This shows the advantages of our approach. Thus, each AUV can accurately compensate for model uncertainties and external disturbances of its corresponding subsystem in finite time.

$$\begin{cases}
\dot{\sharp}\_{i1} = \dot{\sharp}\_{i2} - \beta\_{i1} a\_{i1} \mathbf{e}\_{i1} + \mathbf{G}\_{i}(\eta\_{i\prime} \mathbf{e}\_{i}) + \mathbf{\tau}\_{i} \\
\dot{\sharp}\_{i2} = -\beta\_{i1}^{2} a\_{i2} \mathbf{e}\_{i1}
\end{cases}.\tag{67}$$

$$\begin{cases}
\dot{\hat{\mathbf{z}}}\_{i1} = \hat{\mathbf{z}}\_{i2} - \beta\_{i1} \text{sig}^{3/4}(\mathbf{e}\_{i1}) + \mathbf{G}\_{i}(\eta\_{i}, \mathbf{v}\_{i}) + \mathbf{\tau}\_{i} \\
\dot{\hat{\mathbf{z}}}\_{i2} = -\beta\_{i2} \text{sig}^{1/2}(\mathbf{e}\_{i1})
\end{cases}.\tag{68}$$

**Figure 4.** Compound disturbance estimation errors *ei*<sup>2</sup> of the *i*th AUV.

The collision avoidance performance of the AUV formation was tested via a set of comparison experiments with and without collision avoidance constraints based on our proposed scheme. Since the initial positions of the four AUVs are randomly distributed, the risk of collision is increased. The formation trajectory without collision avoidance constraints is shown in Figure 5 (top). Here, the four AUVs track the reference trajectory while keeping the preset shape, but AUV3 and AUV4 collide at 10 s, followed by a collision between AUV1 and AUV2 at 20 s. Specifically, as presented in Figure 6 (top), the relative distance between AUV1 and AUV2 during the formation configuration stage exceeds the safe distance, resulting in a collision. The same situation occurs with AUV3 and

AUV4. However, when collision avoidance constraints are considered, the formation trajectory (shown in Figure 5 (bottom)) indicates that the four AUVs can perform the formation tracking task while avoiding collision during the configuration stage. The collision avoidance performance is visualized in Figure 6 (bottom), where the distances among AUVs within the detection zone are always greater than the safe distance, indicating that inter-vehicle collision avoidance can be achieved. Therefore, the proposed control scheme can provide real-time collision avoidance capability for AUV formation maneuvers.

**Figure 5.** 3-D formation trajectories without (**top**) and with (**bottom**) collision avoidance constraints.

**Figure 6.** Relative distances among AUVs without (**top**) and with (**bottom**) collision avoidance constraints.

In order to assess the feasibility and superiority of the proposed scheme, we conducted three sets of comparative simulations with the same parameters, constraints, and disturbance settings: (a) the proposed FTESO-based dual closed-loop DMPC with Laguerre function; (b) a FTESO-based dual closed-loop DMPC without Laguerre function; (c) a standard DMPC without FTESO. Figures 7–16 plot the tracking performance curves of AUV positional and velocity states under the three schemes. It can be easily observed that, in all scenarios, the four AUVs are able to successfully track the desired state despite the differing tracking errors. In scheme (a), full-state stable tracking is achieved within 200 s. Meanwhile, in scheme (b), the process takes about 300 s, which suggests that the use of the Laguerre function improves both the response speed and control accuracy. Although the standard DMPC scheme (c) can also achieve formation tracking, the settling time of the state variables is longer and accompanied by oscillations due to the uncompensated compound disturbance effects. Compared with the other schemes, our proposed method delivers superior formation tracking control performance.

**Figure 7.** State *xi* of the AUV formation system. (**a**) The proposed FTESO-based dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

**Figure 8.** State *yi* of the AUV formation system. (**a**) The proposed FTESO-based dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

**Figure 9.** State *zi* of the AUV formation system. (**a**) The proposed FTESO-based dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

**Figure 10.** State *θ<sup>i</sup>* of the AUV formation system. (**a**) The proposed FTESO-based dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

**Figure 11.** State *ψ<sup>i</sup>* of the AUV formation system. (**a**) The proposed FTESO-based dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

**Figure 12.** State *ui* of the AUV formation system. (**a**) The proposed FTESO-based dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

**Figure 13.** State *vi* of the AUV formation system. (**a**) The proposed FTESO-based dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

**Figure 14.** State *wi* of the AUV formation system. (**a**) The proposed FTESO-based dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

**Figure 15.** State *qi* of the AUV formation system. (**a**) The proposed FTESO-based dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

**Figure 16.** State *ri* of the AUV formation system. (**a**) The proposed FTESO-based dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

Figure 17 intuitively presents a 3-D formation trajectory tracking. Combined with Figures 7–16, it implies that all three schemes can successfully accomplish formation spiral tracking under the specified input and state constraints. However, when the formation faces harsh compound disturbances, the tracking performance of the controller without disturbance compensation performs poorly, demonstrating a tracking error significantly larger than that of the FTESO-based controller. This is because the compound disturbances cause the AUV to deviate from the desired trajectory. By comparing the results of (a,b), it can be further observed that the proposed control scheme with Laguerre function allows the AUV to form the preset formation more quickly and converge to the desired trajectory more smoothly. This implies a faster response at the onset of the task. Thus, the dual closed-loop structure and Laguerre function enable the AUV formation to track the reference trajectory with better speed and accuracy.

**Figure 17.** 3-D trajectories of the AUV formation under different schemes. (**a**) The proposed FTESObased dual closed-loop DMPC with Laguerre function. (**b**) The FTESO-based dual closed-loop DMPC without Laguerre function. (**c**) The standard DMPC without FTESO.

Without loss of generality, Figure 18 shows the actual control forces and moments versus time for AUV1 under the three schemes. Without the benefit of FTESO to compensate for compound disturbances, the fluctuations of the control force and moment are relatively drastic and unstable (Figure 18(c1,c2)). This is attributed to the need for the AUV to significantly rectify the driving forces and moments to more rapidly approach the deviated reference trajectory. Under the proposed scheme, as shown in Figure 18(a1,a2), the AUV forces and moments vary relatively smoothly, which makes the AUV track the trajectory steadily when disturbed. Comparing Figure 18(a1,a2) and Figure 18(b1,b2), the Laguerrebased controller has the fastest control signal response with the smallest amplitude when the disturbances are accurately compensated. This confirms that our proposed scheme (a) provides superior control performance. It is worth noting that the variation of control forces and moments always remains within the prescribed limits. This reflects the ability of the DMPC to handle the actuator saturation effectively, ensuring that the control input for each DOF does not exceed the maximum force provided by the actuator, thus reducing actuator losses.

**Figure 18.** Actual control force and moment for AUV1 under different schemes. (**a1**) Control force of scheme (a). (**a2**) Control moments of scheme (a). (**b1**) Control force of scheme (b). (**b2**) Control moments of scheme (b). (**c1**) Control force of scheme (c). (**c2**) Control moments of scheme (c).

To differentiate between the computational demands among the three schemes, we recorded the emulator execution times under the same configurations. The detailed simulation times corresponding to Figure 17 are given in Table 2. It can be observed that the actual running time of the standard DMPC system is approximately 43.62 s. In contrast, the

system with a dual closed-loop DMPC requires 57.91 s, which is about 32.8% longer than the standard DMPC. This increase is due to the greater complexity of the dual closed-loop structure as opposed to the simpler DMPC structure. Although there is improvement in control efficacy, the execution of the dual closed-loop structure is sacrificed to some extent. However, the proposed system, which employs a Laguerre-based dual closed-loop DMPC, the computation time only requires 11.75 s. This suggests that, despite the inclusion of both the dual closed-loop structure and FTESO, the use of the Laguerre function makes the system solution faster. Thus, the proposed scheme can simultaneously improve the computational speed and control performance.

**Table 2.** Comparison of controller execution time.


#### **5. Conclusions**

In conclusion, this study presents a FTESO-based distributed dual closed-loop model predictive control scheme for the AUV formation subject to compound disturbances. The designed FTESO can compensate for model uncertainties and external disturbances of each AUV faster and more accurately. Control inputs are determined by solving a constrained DMPC optimization problem based on local information, while avoiding both collisions among AUVs and actuator saturation. The Laguerre orthogonal function is applied to alleviate the heavy computational burden, and the corresponding stability proof is provided. Finally, based on a connected directed topology, simulation results of different schemes are investigated under the same compound disturbances and system constraints. It is confirmed that our proposed scheme provides the best tracking effect and superior active disturbance rejection capability. Control signals show smaller oscillations and enhanced stability. In addition, the computation time of our proposed formation control system, which utilizes the Laguerre function, is reduced by 73.1% and 79.7% compared to the standard DMPC system and the dual closed-loop DMPC system, respectively. This verifies that our proposed scheme can respond quickly to minimize control costs and improve real-time execution and dynamic performance of the system.

The proposed method does not consider the impact of communication burden on AUV formation. Therefore, in future work, we will focus on the control scheme based on the event-triggered mechanisms. Considering the limitations of the optimization accuracy of discrete predictive control, we want to carry out research on continuous predictive control with faster control response. In addition, it is essential to conduct formation obstacle avoidance research and real AUV experiments.

**Author Contributions:** Conceptualization, M.Z. and Z.Y.; Data curation, J.Z. and L.Y.; Funding acquisition, Z.Y.; Investigation, J.Z. and L.Y.; Methodology, M.Z.; Resources, Z.Y.; Software, L.Y.; Validation, M.Z.; Visualization, J.Z.; Writing—original draft, M.Z.; Writing—review & editing, M.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by National Natural Science Foundation of China under grant No. 52071102, and in part by National Natural Science Foundation of China under grant No. 51679057.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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