Publish and Subscribe-Based Formation and Containment Control of Heterogeneous Robotic System with Actuator Time Delay

Abstract

This paper addresses the problem of controlling a heterogeneous system composed of multiple Unmanned Aerial Vehicles (UAVs) and Autonomous Underwater Vehicles (AUVs) for formation and containment maintenance. The proposed approach considers actuator time delay and, in addition to formation and containment, considers obstacle avoidance, and offers a robust navigation algorithm and uses a reliable middleware for data transmission and exchange. The methodology followed uses both flocking technique and modified L1 adaptive control to ensure the proper navigation and coordination while avoiding obstacles. The data exchange between all the agents is provided through the data distribution services (DDS) middleware, which solves the interoperability issue when dealing with heterogeneous multiagent systems. The modified L1 controller is a local controller for stabilizing the dynamic model of each UAV and AUV, and the flocking approach is used to coordinate the followers around the leader or within the space delimited by their leaders. Potential Field (PF) allows obstacle avoidance during the agents’ movement. The performance of the proposed approach under the considerations mentioned above are verified and demonstrated using simulations.

1. Introduction

1.1. Motivation

Recently, the Unmanned Aerial Vehicle systems (UAVs) and/or Autonomous Underwater Vehicle systems (AUVs) have been involved in numerous applications, including both military and civilian applications such as defense, firefighting, pipeline inspections, exploration, photography, and monitoring. Researchers and engineers are showing great interest in UAVs and/or AUVs systems. This interest is mainly motivated by observable expansions of UAVs and/or AUVs applications, particularly now that UAVs are playing a significant role in combating the COVID-19 pandemic [1]. Due to the multiple advantages of UAV and/or AUV systems and their potential to perform complex tasks with high and reliable efficiency, multi-UAVs and AUVs are preferred to accomplish the assigned tasks more successfully and precisely compared to single-agent UAV and/or AUV systems. Numerous control algorithm designs have been reported in the literature to meet the application criteria either for a group of homogeneous or heterogeneous vehicle systems. Solutions are proposed for heterogeneous robotic agents to be allocated to corresponding tasks that impose specific execution requirements [2]. An optimal sequence of task–agent allocation to minimize the task allocation time or the consumption of the robots’ energy for a team of heterogeneous mobile robotic agents is introduced in [3]. The formation containment control designs of multi-leader follower systems have been proposed and reported in [4,5]; however, the influence of the information transmission between the agents was not considered in the design. Classical formation control designs of multi-agent systems with limited agents have been proposed by implanting controllers such as a state feedback controller, a Fuzzy Logic Control, an LQR controller, and a feedback linearization control technique [6,7]. However, the presence of parameters’ uncertainties and external disturbances may critically influence the performances of the previously mentioned controllers. DDS middleware has been proposed and designed to overcome the issue of the data loss between multiple UAVs, as well as the interoperability issue associated with heterogeneous systems. Therefore, application of the DDS middleware technique augmented with the flocking approach resulted in significant improvements in flexibility, reliability, and portability [8,9].

1.2. Related Research

Generally, multi-vehicle formation control or containment control design can be categorized based on the following three structures: virtual-based design (see [10,11]), leader-following (see [12,13]), and behavior-based [14,15]. Multi-mobile robots formation control design has been proposed in [16] where the model predictive control (MPC) for the robot’s trajectory tracking is implemented together with the adaptive terminal sliding mode control (ATSMC) for tracking the robot’s velocity. The Pythagorean theorem has been applied for cooperation between mobile robots. The proposed controller performance has been evaluated in the presence of external disturbances that may influence the mobile robot model. Relevant to this work, the formation control of multiple ground vehicle systems has been addressed in [17] and in more detail in [18]. In these studies, the authors build on the literature results to address the formation maintenance and control of several robots where the leader or virtual leader is always at the center of the formation. The control design that is used locally is a feedback linearization, and the team navigation is stabilized based on a potential field force. Multi-mecanum-wheeled mobile robots’ (MWMRs’) leader–follower formation control design has also been presented in [19], where the fuzzy sliding mode controller for trajectory tracking has been implemented considering the parameters’ uncertainties. The graph theorem has been applied for consensus between the heterogeneous systems of MWMRs. The consensus control design of a heterogeneous system of multi unmanned ground vehicles (UGVs) and multi unmanned aerial vehicles (UAVs) has been reported in [20], where the particle swarm optimization (PSO) approach was implemented for consensus between UGVs and UAVs. Furthermore, the issue of the path planning problem of multi UGVs and UAVs has been considered, while the consensus of the UGVs and the UAVs system has been investigated considering some assumptions and constraints. The leader–follower control design of multi-robot manipulator heterogeneous system has been proposed and reported in [21], where the leader–follower approach of multi-robot manipulator heterogeneous system has been studied in the presence of the parameters’ uncertainties, while the leaders’ information has been estimated by implementing the cooperative estimator technique. In addition, a graph theory approach has been utilized for consensus between the multi-robot manipulator’s heterogeneous system. The framework to integrate multi-UAV applications has been developed and reported in [22]; the new proposed framework GzUavChannel is basically a combination of the ns-3, Gazebo, and ArduPilot network simulators. The GzUavChannel middleware can guarantee the integration and communication between multiple interconnected computers with different simulator tools such as ns-3, Gazebo, and ArduPilot. The framework for the integration between multi-ground robots in real-time applications has been introduced and proposed in [23], where a Domain Specific Language (DSL) and constraints analysis of the software applications’ safety is used to provide a path between the design of architectures. The authors applied the MAUVE DSL deployment to define the model of execution; however, it is necessary to define some information about Worst-Case Execution Time (WCET) at the onset so that the Extreme Value Theory (EVT) is used to estimate the WCET execution time. The Orocos middleware is used to integrate and exchange information between the multi-ground robots. The framework of multip;e AUVs based on Mission-Oriented Operating Suite-Interval Programming (MOOSIvP) middleware has been developed in [24]. The publish–subscribe framework has been applied to exchange information between the AUVs. The authors identified the object as a sonar image for each AUV by utilizing the Sidescan Sonar Sensor (SSS). A classical PID controller has been implemented on each AUV to track the elevator, rudder, and thrust references. From the literature coverage, one can deduce that:

• the influence of the information transmission between the agents has been identified as a crucial factor in performance and reliability of the control design;
• the presence of parameters’ uncertainties and external disturbances may critically influence the performances of the prior.;
• application of the DDS middleware technique augmented with the flocking approach has not been applied for heterogeneous systems, which is closer to the real-word robotic applications.

This paper proposes an approach that can possibly address all of the above points. Indeed, navigation and control for formation and containment mainteannce of a fleet of multiagent heterogeneous systems composed of AUVs and UAVs is considered. A distributed method has been developed using an L1 adaptive control technique taking into account actuator time delay to stabilize the agents and make them more robust towards parameters change. The flocking approach has been used for the navigation, and the DDS middleware has been utilized for the data transmission between the agents. The navigation ensures formation as well as containment control in the leader–follower context. Potential field-like methodology has been adopted to avoid obstacle while maintaining the formation and containment.

The main contributions of this research are highlighted as follows:

• A robust architecture for the formation control design of heterogeneous systems has been developed.
• A robust navigation system for the containment control design of heterogeneous systems has been proposed and designed.
• The navigation algorithm avoids obstacles while still maintaining the required formation or containment.
• Reliable middleware for data transmission has been developed.

The remainder of this manuscript is organized as follows. In Section 2, the preliminaries of DDS middleware, the dynamics of the quadrotor, and the dynamics of the underwater system are presented. In Section 3, the multi UAVs formation stability has been proved using the L1 technique with an actuator time delay and Lyapunov analysis. The simulation results are presented and discussed in Section 4. Finally, the findings and the concluding remarks have been highlighted at the end.

2. Preliminaries

As mentioned above, the system that this work is considering is composed of a heterogeneous fleet of UAVs and AUVs connected using DDS middleware. This section will cover the necessary background of DDS middleware and then present the model for each type of agent (AUV and UAV). The modeling background is necessary for the development of the local agent controller, which will be presented in Section 3. In this section, the following symbols in

Table 1. List of the main symbols.

For the UAV
${\eta }_1\in {\mathbb{R}}^{3\times 1}$	the vector representing z, y, and x coordinates.
${\eta }_2\in {\mathbb{R}}^{3\times 1}$	the vector of Euler angles yaw ($\psi$), pitch ($\theta$), and roll ($\varphi$).
g	represents the gravity acceleration.
m	represents the mass.
$k_t$	is the drag coefficient.
u	is the total thrust force input created by the 4 propellers.
$\mathbf{I}\in {\mathbb{R}}^{3\times 3}=\mathrm{diag}(I_x,I_y,I_z)$	represents the quadrotor’s inertia.
${\mathbf{\nu }}_2$	is the angular velocity.
$z_e$	is equal to ${[0,0,1]}^T$.
$I_R$	represents the propeller inertia.
$k_r$	represent the rotational drag.
×	is the cross product.
For the AUV
$M\in {\mathbb{R}}^{6\times 6}$	is the inertia matrix.
$C\left(\nu \right)\in {\mathbb{R}}^{6\times 6}$	is the coriolis matrix.
$D\left(\nu \right)\in {\mathbb{R}}^{6\times 6}$	is the damping matrix.
$g\left(\eta \right)\in {\mathbb{R}}^6$	is the gravitational vector.
$\tau \in {\mathbb{R}}^{6\times 1}$	is the input control signals.
$M_{RB}$	is the rigid body’s inertia matrix.
$M_A$	represents the additional momentum and forces.
$r_g={\left[r_{g_1}{r}_{g_2}{r}_{g_3}\right]}^{\top }$	is the gravity vector.
$I_{3\times 3}$	is the 3-by-3 identity matrix.
$(I_0=I_0^{\top }>0)$	is the inertia tensor.

will be used.

2.1. Data Distribution Service (DDS)

DDS middleware has emerged as a solution for many novel distributed frameworks. It boosts the publish–subscribe technique for transmitting data, commands, and events between the agents. Agents that are producing data publish the data, and they can subscribe to the data from other agents. The middleware is responsible for providing the data to subscribers who subscribed to that topic. It has been shown that DDS middleware is a very effective approach for many-to-many communication networks as it is very successful in transmitting time-critical information. DDS offers a group of Quality-of-Service (QoS) regulations to adapt to communication behaviors. These QoS regulations might be used separately or as a set in order to influence different aspects of communications, such as reliability, durability, and data persistence. These QoS regulations are a feature of the DDS which are not available in all other middleware approaches to communication. This feature guarantees high reliability, performance, low latency integration of DDS data transmission, and DDS interoperability with the capability of DDS frameworks to communicate from various vendors. The standardized DDS Real-Time Publish–Subscribe (RTPS) protocol allows for connection interoperability between various DDS frameworks. The DDS middleware was considered for our study for the formation and containment control of heterogeneous systems. Each agent must publish its velocity and subscribe to other agents’ velocities. A selection of Quality of Service regulations must be fulfilled for interaction between the agents (publisher/subscribers), such as:

• Durability: Determines whether or not the previous data sent by the publisher has been received by a new subscriber.
• Reliability: Specifies whether or not the publisher will resend information lost over the network. Reliability involves two setups: BEST EFFORT (not re-send data lost) and RELIABLE (resend data lost).
• History: Determines a publisher’s preservation of the data received or sent to a subscriber. History involves two setups: KEEP LAST and KEEP ALL.
• Ownership: Determines whether a subscriber will simultaneously accept new samples from multiple publishers. Ownership involves two setups: SHARED and EXCLUSIVE.

2.2. Dynamics of Quadrotor

A detailed dynamic model of the quadrotor can be found in [25].

The quadrotor’s translational motion is given by

$${\ddot{\eta }}_1=-g\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]+J_1\left({\eta }_2\right)\left[\begin{array}{c}0\\ 0\\ u/m\end{array}\right]-\frac{k_t}{m}{\dot{\eta }}_1,$$

(1)

The transformation matrix $J_1\left({\eta }_2\right)$ relates to the Euler angles according to

$$J_1\left({\eta }_2\right)=\left[\begin{array}{ccc}c\psi c\theta & -s\psi c\varphi +c\psi s\theta s\varphi & s\psi s\varphi +c\psi c\varphi s\theta \\ s\psi c\theta & c\psi c\varphi +s\varphi s\theta s\psi & -c\psi s\varphi +s\psi c\varphi s\theta \\ -s\theta & c\theta s\varphi & c\varphi c\theta \end{array}\right].$$

(2)

The property of $J_1\left({\eta }_2\right)$ is

$$\begin{array}{c}\hfill J_1^{-1}\left({\eta }_2\right)=J_1^T\left({\eta }_2\right),\end{array}$$

(3)

where $s(.)=sin(.),c(.)=cos(.)$. The input u is defined as follows:

$$\begin{array}{c}\hfill u=f_1+f_2+f_3+f_4,\end{array}$$

(4)

where an upward-lifting force is $f_i=k_i{\mathsf{\Omega }}_i^2$, i = 1, 2, 3, and 4, with positive constants $k_i$ and ${\mathsf{\Omega }}_i$ represents the angular speed motors. The quadrotor’s rotational motion is given by

$$\begin{array}{c}\hfill {\dot{\nu }}_2=I^{-1}(-({\nu }_2\times I{\nu }_2)-I_R({\nu }_2\times z_e)\mathsf{\Omega }-k_r{\nu }_2+\tau ),\end{array}$$

(5)

and

$$\begin{array}{c}\hfill \mathsf{\Omega }={\mathsf{\Omega }}_1-{\mathsf{\Omega }}_2+{\mathsf{\Omega }}_3-{\mathsf{\Omega }}_4.\end{array}$$

(6)

The force u and the torques $\tau$$={[{\tau }_p,{\tau }_q,{\tau }_r]}^{\top }$ on the frame’s body. The quadrotor’s rotational and translational motion is given by

$$\left[\begin{array}{c}\tau \\ u\end{array}\right]=\left[\begin{array}{c}{\tau }_p\\ {\tau }_q\\ {\tau }_r\\ u\end{array}\right]=\left[\begin{array}{cccc}0& l& 0& -l\\ l& 0& -l& 0\\ d& -d& d& -d\\ 1& 1& 1& 1\end{array}\right]\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_4\end{array}\right],$$

(7)

where d represents the drag factor and l represents the distance between motors and the center of mass. The number actuation from the above equations is less than the degree of freedom. Therefore, the following description of underactuated structures is taken into account in our case. Consider the next system:

$$\ddot{Q}=F(Q,\dot{Q})+G\left(Q\right)U,$$

(8)

where F is a dynamic vector, Q is coordinates vector, U is control inputs, and G is input matrix. If G’s rank is less than Q’s dimensions, then Equation (8) is under-actuated. A detailed quadrotor dynamic is available in [8,9,25].

2.3. Dynamics of Underwater Vehicles

A detailed dynamic model of the underwater vehicle can be found in [26].

The underwater motion equation in the body frame is given by:

$$M\dot{\nu }+C\left(\nu \right)\nu +D\left(\nu \right)\nu +g\left(\eta \right)=\tau$$

(9)

The inertia matrix has two parts:

$$\begin{array}{c}\hfill M=M_{RB}+M_A,\end{array}$$

(10)

and

$$M_{RB}=\left[\begin{array}{cc}m{I}_{3\times 3}& -mS\left(r_g\right)\\ mS\left(r_g\right)& I_0\end{array}\right],$$

(11)

The skew symmetric matrix $S\left(r_g\right)$ is given by

$$S\left(r_g\right)=\left[\begin{array}{ccc}0& -r_{g_3}& r_{g_2}\\ r_{g_3}& 0& -r_{g_1}\\ -r_{g_2}& r_{g_1}& 0\end{array}\right].$$

(12)

In addition,

$$M_A=\left[\begin{array}{cccccc}{X}_{\dot{u}}& X_{\dot{v}}& X_{\dot{w}}& X_{\dot{p}}& X_{\dot{q}}& X_{\dot{r}}\\ Y_{\dot{u}}& Y_{\dot{v}}& Y_{\dot{w}}& Y_{\dot{p}}& Y_{\dot{q}}& Y_{\dot{r}}\\ Z_{\dot{u}}& Z_{\dot{v}}& Z_{\dot{w}}& Z_{\dot{p}}& Z_{\dot{q}}& Z_{\dot{r}}\\ K_{\dot{u}}& K_{\dot{v}}& K_{\dot{w}}& K_{\dot{p}}& K_{\dot{q}}& K_{\dot{r}}\\ M_{\dot{u}}& M_{\dot{v}}& M_{\dot{w}}& M_{\dot{p}}& M_{\dot{q}}& M_{\dot{r}}\\ N_{\dot{u}}& N_{\dot{v}}& N_{\dot{w}}& N_{\dot{p}}& N_{\dot{q}}& N_{\dot{r}}\end{array}\right],$$

(13)

where, the torques and the linear forces applied on the underwater model are $x_{.}=\frac{\partial x}{\partial (.)}$. X, Y, Z, K, M, N. For example, the hydrodynamic force in x-direction due to acceleration $\dot{u}$ is represented by $X_{\dot{u}}$, so the inertia matrix can be used to determine the Coriolis matrix by constructing it in terms of sub-matrices:

$$M=\left[\begin{array}{cc}{M}_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right].$$

(14)

The Coriolis and damping matrices are thus given respectively by:

$$C\left(\nu \right)=\left[\begin{array}{cc}{0}_{3\times 3}& -S(M_{11}{\nu }_1+M_{12}{\nu }_2)\\ -S(M_{11}{\nu }_1+M_{12}{\nu }_2)& -S(M_{21}{\nu }_1+M_{22}{\nu }_2)\end{array}\right],$$

(15)

$$D\left(\nu \right)=\left[\begin{array}{cccccc}|{\nu }_1|& 0& 0& 0& 0& 0\\ 0& |{\nu }_2|& 0& 0& 0& 0\\ 0& 0& |{\nu }_3|& 0& 0& 0\\ 0& 0& 0& |{\nu }_4|& 0& 0\\ 0& 0& 0& 0& |{\nu }_5|& 0\\ 0& 0& 0& 0& 0& |{\nu }_6|\end{array}\right].$$

(16)

The underwater equation of motion in the Earth frame is given by:

$$M_e\left(\eta \right)\ddot{\eta }+C_e(\eta ,\nu )\dot{\eta }+D_e(\nu ,\eta )\dot{\eta }+g_e\left(\eta \right)={\tau }_e,$$

(17)

where

$$\begin{array}{cc}\hfill M_e\left(\eta \right)=& R^{-T}\left(\eta \right)M\left(\eta \right)R^{-1}\left(\eta \right)\in {\mathbb{R}}^{6\times 6}\hfill \\ \hfill C_e(\eta ,\nu )=& R^{-T}\left(\eta \right)[C(\eta ,\nu )-M\left(\eta \right)R^{-1}\left(\eta \right)\dot{R}\left(\eta \right)]R^{-1}\left(\eta \right)\hfill \\ \hfill C_e(\eta ,\nu )\in & {\mathbb{R}}^{6\times 6}\hfill \\ \hfill D_e(\nu ,\eta )=& R^{-T}\left(\eta \right)D\left(\nu \right)R^{-1}\left(\eta \right)\in {\mathbb{R}}^{6\times 6}\hfill \\ \hfill g_e\left(\eta \right)=& R^{-T}\left(\eta \right)g\left(\eta \right)\in {\mathbb{R}}^6\hfill \\ \hfill {\tau }_e=& R^{-T}\left(\eta \right)\tau \in R^6,\hfill \end{array}$$

and the transformation matrix $R\left(\eta \right)$ is as follows

$$R\left(\eta \right)=\left[\begin{array}{cc}{R}_{1_{3\times 3}}\left(\eta \right)& 0_{3\times 3}\\ 0_{3\times 3}& R_{2_{3\times 3}}\left(\eta \right)\end{array}\right],$$

(18)

with

$$R_1\left(\eta \right)=\left[\begin{array}{ccc}{C}_{\psi }{C}_{\theta }& C_{\varphi }{S}_{\psi }-C_{\psi }{S}_{\varphi }{S}_{\theta }& -S_{\varphi }{S}_{\psi }-C_{\varphi }{C}_{\psi }{S}_{\theta }\\ C_{\theta }{S}_{\psi }& -C_{\varphi }{C}_{\psi }-S_{\varphi }{S}_{\psi }{S}_{\theta }& C_{\psi }{S}_{\varphi }-C_{\varphi }{S}_{\psi }{S}_{\theta }\\ -S_{\theta }& -C_{\theta }{S}_{\varphi }& -C_{\varphi }{S}_{\theta }\end{array}\right]$$

and

$$R_2\left(\eta \right)=\left[\begin{array}{ccc}0& -S_{\varphi }/C_{\theta }& C_{\varphi }/C_{\theta }\\ 0& -C_{\varphi }& S_{\varphi }\\ 1& -S_{\varphi }{S}_{\theta }/C_{\theta }& -C_{\varphi }{S}_{\theta }/C_{\theta }\end{array}\right].$$

A detailed description of the underwater dynamic can be found in [26].

3. Multiple AUVs-UAVs Control Design

3.1. L1 Adaptive Control with Actuator Time Delay

In this section, the L1 adaptive control utilized in [26] is used for a single AUV extended with an actuator time delay to stabilize multiple AUVs.

The main feature of utilizing the L1 adaptive controller with an actuator time delay is its ability to adapt rapidly and robustly, which augments the performance of the system. This capacity has been achieved by separating adaption from robustness. Utilizing a modified L1 controller, the uncertainties are predestined by using a fast method and the recompensed uncertainties proceed through a filter. Then, the overall output is included within the input signal. There are two key functions to this filter. First, it ensures that the control signal stays within an acceptable frequency range. Second, it separates the robustness from the adaptation. With this filter, we can improve robustness as well as performance. For more details, see [27,28,29].

In order to implement the modified L1 controller to the quadrotor/underwater vehicle, a representation of the state space for (5) and (17) is required.

$$\begin{array}{cc}\hfill & {\dot{x}}_1\left(m\right)=x_2,\mathrm{with}{x}_1\left(0\right)=x_{1_0}\hfill \\ \hfill & {\dot{x}}_2\left(m\right)=A_2{x}_2\left(m\right)+f_2(m,x\left(m\right))+B_{2}wu(m-\tau ),\mathrm{with}{x}_2\left(0\right)=x_{2_0}\hfill \\ \hfill & y\left(m\right)=Cx\left(m\right),\hfill \end{array}$$

(19)

where, the system’s states are denoted by $x\left(m\right)={[x_1\left(m\right),x_2\left(m\right)]}^{\top }\in {\mathbb{R}}^{2n}$, the known matrix is represented by $A_2\in {\mathbb{R}}^{n\times n}$, the constant matrix is represented by $B_2\in {\mathbb{R}}^{m\times n}$, $u(m-\tau )\in {\mathbb{R}}^m$ represents the control vector with time delay $\tau$, $w\in \mathbb{R}$ denote the input’s uncertainty, the known constant matrix is represented by $C\in {\mathbb{R}}^{m\times n}$, the system’s output is represented by $y\left(m\right)\in {\mathbb{R}}^m$, and an unknown nonlinear function is denoted by $f_2\left(m\right)$. We can rewrite (19) as follows [30]:

$$\dot{x}=Ax\left(m\right)+f(m,x)+B_{m}wu(m-\tau ),$$

(20)

with

$$A=\left[\begin{array}{cc}{0}_{n\times n}& I_{n\times n}\\ 0_{n\times n}& A_2\end{array}\right],f=\left[\begin{array}{c}{0}_{n\times 1}\\ f_2\end{array}\right]\mathrm{and}{B}_m=\left[\begin{array}{c}{0}_{n\times m}\\ B_2\end{array}\right].$$

From (17), one can have

$$\begin{array}{c}\hfill \ddot{\eta }=M_e{\left(\eta \right)}^{-1}[{\tau }_e-C_e(\eta ,\nu )\dot{\eta }-D_e(\nu ,\eta )\dot{\eta }-g_e\left(\eta \right)].\end{array}$$

(21)

Let

$$\begin{array}{cc}\hfill & {\eta }_1=\eta \hfill \\ \hfill & {\eta }_2=\dot{\eta },\mathrm{thus}\hfill \\ \hfill & \dot{{\eta }_1}=\dot{\eta }={\eta }_2\hfill \\ \hfill & \dot{{\eta }_2}=\ddot{\eta }.\hfill \end{array}$$

(22)

Therefore, the state space representation can be obtained from (21) and (22) as follows:

$$\begin{gathered} \left[\begin{array}{c}{\dot{\eta }}_1\\ {\dot{\eta }}_2\end{array}\right]=\left[\begin{array}{cc}{0}_{6\times 6}& I_{6\times 6}\\ 0_{6\times 6}& \frac{-D_e}{M_e}\end{array}\right]\left[\begin{array}{c}{\eta }_1\\ {\eta }_2\end{array}\right]-\left[\begin{array}{c}{0}_{6\times 1}\\ \frac{g_e}{M_e}\end{array}\right]+\left[\begin{array}{c}{0}_{6\times 6}\\ \frac{1}{M_e}\end{array}\right]\tau \hfill & \\ \hfill y\left(m\right)& =\left[\begin{array}{cc}{I}_{6\times 6}& 0_{6\times 6}\\ 0_{6\times 6}& I_{6\times 6}\end{array}\right]\left[\begin{array}{c}{\eta }_1\\ {\eta }_2\end{array}\right].\hfill \end{gathered}$$

(23)

Equation (23) can be rewritten as

$$\begin{array}{cc}\hfill \dot{x}\left(m\right)=& A_{m}x\left(m\right)+B_m(w\left(m\right)u_{ad}(m-\tau )+\theta \left(m\right){\parallel x\left(m\right)\parallel }_{\infty }\hfill \\ \hfill +& \sigma \left(m\right)),\mathrm{with}x\left(0\right)=x_0\hfill \\ \hfill y\left(m\right)=& c^{\top }x\left(m\right).\hfill \end{array}$$

(24)

Therefore, the dynamic of the state predictor $\hat{x}\left(m\right)\in {\mathbb{R}}^n$ can be defined as:

$$\begin{array}{cc}\hfill \dot{\widehat{x}}\left(m\right)=& A_m\widehat{x}\left(m\right)+B_m(\widehat{w}\left(m\right)u_{ad}(m-\tau )+\widehat{\theta }\left(m\right){\parallel x\left(m\right)\parallel }_{\infty }\hfill \\ \hfill +& \widehat{\sigma }\left(m\right))\mathrm{with}\widehat{x}\left(0\right)={\widehat{x}}_0\hfill \\ \hfill \widehat{y}\left(m\right)=& c^{\top }\widehat{x}\left(m\right),\hfill \end{array}$$

(25)

with the estimated parameters denoted by $\widehat{w}\left(m\right)$, $\widehat{\theta }\left(m\right)$ and $\widehat{\sigma }\left(m\right)$. Hence, the error can be defined as $\tilde{x}\left(m\right)=x\left(m\right)-\widehat{x}\left(m\right)$, $\tilde{w}\left(m\right)=w\left(m\right)-\widehat{w}\left(m\right)$, $\tilde{\theta }\left(m\right)=\theta \left(m\right)-\widehat{\theta }\left(m\right)$, and $\tilde{\sigma }\left(m\right)=\sigma \left(m\right)-\widehat{\sigma }\left(m\right)$. Therefore, the error dynamics can be found as follows:

$$\begin{array}{cc}\hfill \dot{\tilde{x}}\left(m\right)=& A_m\tilde{x}\left(m\right)+B_m(\tilde{w}\left(m\right)u_{ad}(m-\tau )+\tilde{\theta }\left(m\right){\parallel x\left(m\right)\parallel }_{\infty }\hfill \\ \hfill +& \tilde{\sigma }\left(m\right)).\hfill \end{array}$$

(26)

We build on same analysis of the L1 adaptive controller applied to a single AUV in [26] by adding some modifications based on the actuator time delay to make the proposed approach more comprehensive.

Consider the following Lyapunov function

$$\begin{array}{cc}\hfill V(\tilde{x}\left(m\right),\tilde{\theta }\left(m\right),\tilde{\sigma }\left(m\right),\tilde{\sigma }\left(m\right))& ={\tilde{x}}^{\top }\left(m\right)P\tilde{x}\left(m\right)+\frac{{\tilde{w}}^2\left(m\right)}{{\gamma }_3}\hfill \\ \hfill +& \frac{{\tilde{\theta }}^{\top }\left(m\right)\tilde{\theta }\left(m\right)}{{\gamma }_1}+\frac{{\tilde{\sigma }}^{\top }\left(m\right)\tilde{\sigma }\left(m\right)}{{\gamma }_2}.\hfill \end{array}$$

(27)

The system stability is guaranteed for any u such that $\dot{V}\left(m\right)\le 0$. Consequently,

$$\begin{array}{cc}\hfill \dot{V}\left(m\right)=& {\dot{\tilde{x}}}^{\top }\left(m\right)P\tilde{x}\left(m\right)+{\tilde{x}}^{\top }\left(m\right)P\dot{\tilde{x}}\left(m\right)\hfill \\ \hfill +& \frac{1}{{\gamma }_1}({\dot{\tilde{\theta }}}^{\top }\left(m\right)\tilde{\theta }\left(m\right)+{\tilde{\theta }}^{\top }\left(m\right)\dot{\tilde{\theta }}\left(m\right))\hfill \\ \hfill +& \frac{1}{{\gamma }_2}({\dot{\tilde{\sigma }}}^{\top }\left(m\right)\tilde{\sigma }\left(m\right)+{\tilde{\sigma }}^{\top }\left(m\right)\dot{\tilde{\sigma }}\left(m\right))\hfill \\ \hfill +& \frac{2}{{\gamma }_3}\left(\tilde{w}\left(m\right)\dot{\tilde{w}}\left(m\right)\right)\hfill \\ \hfill =& -\tilde{x}{\left(m\right)}^{\top }(A_m^{\top }P+P{A}_m)\tilde{x}\left(m\right)+B_m\tilde{\theta }\left(m\right){\parallel x\left(m\right)\parallel }_{\infty }\hfill \\ \hfill +& 2\tilde{x}{\left(m\right)}^{\top }P(B_m\tilde{w}\left(m\right)u_{ad}(m-\tau )+B_m\tilde{\sigma }\left(m\right))\hfill \\ \hfill +& \frac{2}{{\gamma }_1}\left({\tilde{\theta }}^{\top }\left(m\right)\dot{\tilde{\theta }}\left(m\right)\right)+\frac{2}{{\gamma }_2}\left({\tilde{\sigma }}^{\top }\left(m\right)\dot{\tilde{\sigma }}\left(m\right)\right)+\frac{2}{{\gamma }_3}\left(\tilde{w}\right(m)\dot{\tilde{w}}\left(m\right))\hfill \\ \hfill \mathrm{Therefore},\\ \hfill \dot{V}\left(m\right)=& -\tilde{x}{\left(m\right)}^{\top }Q\tilde{x}\left(m\right)+2\tilde{x}{\left(m\right)}^{\top }P{B}_m\tilde{w}\left(m\right)u_{ad}(m-\tau )\hfill \\ \hfill +& 2\tilde{x}{\left(m\right)}^{\top }P{B}_m\tilde{\theta }\left(m\right){\parallel x\left(m\right)\parallel }_{\infty }+2\tilde{x}{\left(m\right)}^{\top }P{B}_m\tilde{\sigma }\left(m\right)\hfill \\ \hfill +& \frac{2}{{\gamma }_1}\left({\tilde{\theta }}^{\top }\left(m\right)\dot{\tilde{\theta }}\left(m\right)\right)+\frac{2}{{\gamma }_2}\left({\tilde{\sigma }}^{\top }\left(m\right)\dot{\tilde{\sigma }}\left(m\right)\right)+\frac{2}{{\gamma }_3}\left(\tilde{w}\right(m)\dot{\tilde{w}}\left(m\right)).\hfill \end{array}$$

(28)

By utilizing the projection operators, $Proj$, the parameter derivatives can be rewritten as follows:

$$\begin{array}{cc}\hfill \dot{\tilde{\theta }}=& {\gamma }_{1}Proj(\tilde{\theta }\left(m\right),-\parallel x\left(m\right)\parallel B_m^{\top }P\tilde{x}\left(m\right)),\hfill \\ \hfill \dot{\tilde{\sigma }}=& {\gamma }_{2}Proj(\tilde{\sigma }\left(m\right),-B_m^{\top }P\tilde{x}\left(m\right)),\hfill \\ \hfill \dot{\tilde{w}}=& {\gamma }_{3}Proj(\tilde{w}\left(m\right),-x^{\top }\left(m\right)P{B}_m{u}_{ad}(m-\tau )),\hfill \end{array}$$

(29)

where the rate’s adaptation laws are denoted by (${\gamma }_1,{\gamma }_2$ and ${\gamma }_3$) $>0$, $P=P^{\top }>0$ to satisfy the Lyapunov function:

$$\begin{array}{c}\hfill A_m^{\top }P+P{A}_m=-Q\end{array}$$

(30)

and $Q=Q^{\top }>0$. Substituting Equation (29) into Equation (28), one obtains Equation (31):

$$\begin{array}{cc}\hfill \dot{V}\left(m\right)=& -\tilde{x}{\left(m\right)}^{\top }Q\tilde{x}\left(m\right)+2\tilde{x}{\left(m\right)}^{\top }P{B}_m\tilde{w}\left(m\right)u_{ad}(m-\tau )\hfill \\ \hfill +& 2\tilde{x}{\left(m\right)}^{\top }P{B}_m\tilde{\theta }\left(m\right){\parallel x\left(m\right)\parallel }_{\infty }+2\tilde{x}{\left(m\right)}^{\top }P{B}_m\tilde{\sigma }\left(m\right)\hfill \\ \hfill +& 2{\tilde{\theta }}^{\top }\left(m\right)Proj(\tilde{\theta }\left(m\right),-\parallel x\left(m\right)\parallel B_m^{\top }P\tilde{x}\left(m\right))\hfill \\ \hfill +& 2{\tilde{\sigma }}^{\top }\left(m\right)Proj(\tilde{\sigma }\left(m\right),-B_m^{\top }P\tilde{x}\left(m\right))\hfill \\ \hfill +& 2\tilde{w}\left(m\right)Proj(\tilde{w}\left(m\right),-x^{\top }\left(m\right)P{B}_m{u}_{ad}(m-\tau ))\hfill \\ \hfill =& -\tilde{x}{\left(m\right)}^{\top }Q\tilde{x}\left(m\right)+2{\tilde{\theta }}^{\top }\left(m\right){\parallel x\left(m\right)\parallel }_{\infty }{B}_m^{\top }P\tilde{x}\left(m\right)\hfill \\ \hfill +& 2{\tilde{\theta }}^{\top }\left(m\right)Proj(\tilde{\theta }\left(m\right),-{\parallel x\left(m\right)\parallel }_{\infty }{B}_m^{\top }P\tilde{x}\left(m\right))\hfill \\ \hfill +& 2{\tilde{\sigma }}^{\top }\left(m\right)B_m^{\top }P\tilde{x}\left(m\right)\hfill \\ \hfill +& 2{\tilde{\sigma }}^{\top }\left(m\right)Proj(\tilde{\sigma }\left(m\right),-B_m^{\top }P\tilde{x}\left(m\right))\hfill \\ \hfill +& 2\tilde{w}\left(m\right)x^{\top }\left(m\right)P{B}_m{u}_{ad}(m-\tau )\hfill \\ \hfill +& 2\tilde{w}\left(m\right)Proj(\tilde{w}\left(m\right),-x^{\top }\left(m\right)P{B}_m{u}_{ad}(m-\tau )).\hfill \end{array}$$

(31)

Finally, Equation (31) becomes the following:

$$\begin{array}{cc}\hfill \dot{V}\left(m\right)=& -\tilde{x}{\left(m\right)}^{\top }Q\tilde{x}\left(m\right)\hfill \\ \hfill +& 2\tilde{\theta }{\left(m\right)}^{\top }[\stackrel{-y_1}{\overbrace{{\parallel x\left(m\right)\parallel }_{\infty }{B}_m^{\top }P\tilde{x}\left(m\right)}}\hfill \\ \hfill +& Proj(\tilde{\theta }\left(m\right),\stackrel{y_1}{\overbrace{-{\parallel x\left(m\right)\parallel }_{\infty }{B}_m^{\top }P\tilde{x}\left(m\right))}}]\hfill \\ \hfill +& 2\tilde{\sigma }{\left(m\right)}^{\top }[\stackrel{-y_2}{\overbrace{B_m^{\top }P\tilde{x}\left(m\right)}}\hfill \\ \hfill +& Proj(\tilde{\sigma }\left(m\right),\stackrel{y_2}{\overbrace{-B_m^{\top }P\tilde{x}\left(m\right))}}]\hfill \\ \hfill +& 2\tilde{w}\left(m\right)[\stackrel{-y_3}{\overbrace{x^{\top }\left(m\right)P{B}_m{u}_{ad}(m-\tau )}}\hfill \\ \hfill +& Proj(\tilde{w}\left(m\right),\stackrel{y_3}{\overbrace{-x^{\top }\left(m\right)P{B}_m{u}_{ad}(m-\tau )}})],\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \Rightarrow \dot{V}\left(m\right)=& -\tilde{x}{\left(m\right)}^{\top }Q\tilde{x}\left(m\right)\le 0\hfill \end{array}$$

(33)

Therefore, the system is asymptotically stable.

The Laplace transform of the control signal is calculated as follows:

$$\begin{array}{c}\hfill u_{ad}\left(s\right)e^{-\tau s}=-\frac{H\left(s\right)}{\omega }(\widehat{\mu }\left(s\right)-k_{g}r\left(s\right)),\end{array}$$

(34)

where the Laplace of $\widehat{\mu }\left(m\right)=\widehat{\theta }\left(m\right){\parallel x\left(m\right)\parallel }_{\infty }+\widehat{\sigma }\left(m\right)$ and the reference signals $r\left(m\right)$ are $\widehat{\mu }\left(s\right)$ and $r\left(s\right)$, respectively, and the feed-forward gain is represented by $k_g=-\frac{1}{c^{\top }{A}_m^{-1}}{B}_m$.

So, the $H\left(s\right)$ filter has been chosen as follows:

$$\begin{array}{c}\hfill H\left(s\right)=\frac{\omega kC\left(s\right)}{I+\omega kC\left(s\right)}\end{array}$$

(35)

By selecting $C\left(s\right)=\frac{I}{s}$ with $H\left(0\right)=I$, Equation (35) can be written as follows:

$$\begin{array}{c}\hfill H\left(s\right)=\frac{\omega k}{sI+\omega k}\end{array}$$

(36)

The control signal can be obtained by substituting Equation (35) into Equation (34) as follows:

$$\begin{array}{c}\hfill u_{ad}\left(s\right)e^{-\tau s}=-kC\left(s\right)(\widehat{\omega }{u}_{ad}\left(s\right)e^{-\tau s}+\widehat{\mu }\left(s\right)-k_{g}r\left(s\right))\end{array}$$

(37)

Last, the following condition must be satisfied for the L1 controller:

$$\begin{array}{c}\hfill {\parallel G\left(s\right)\parallel }_{{\mathcal{L}}_1}L<1,\end{array}$$

(38)

where

$$\begin{array}{c}\hfill G\left(s\right)=D\left(s\right)(H\left(s\right)-I),D\left(s\right)={(sI-A_m)}^{-1}{B}_m\end{array}$$

(39)

and

$$\begin{array}{c}\hfill \parallel \frac{\partial f(t,x)}{\partial x}\parallel \le d_{f_x}=L.\end{array}$$

(40)

More details can be found in [26,27].

3.2. Flocking and Boids Mode

One of the most interesting sights in nature is the mutual movement of animals in sets, where each animal has its own speed and direction even though they move along with each others to form a single entity. Following social practices, ants, fish, and bird communities are examples of biological entities in their movements. The first method to simulate these movements was established in 1986 by Craig Reynolds and referred to as a bird-like object called the Boids model [31]. The Boids model explains birds’ movement and organization activities in a flock. These cooperative movements typically follow two basic rules: the need for each agent to remain in the group while avoiding collision with the adjacent agent. In order to explain these rules specifically, Reynolds split them into three groups to simplify them for computer simulation [32]. Based on the rules of Reynolds [31], the group N, which consists of n agents, is supposed to represent the whole flock.

$$\begin{array}{c}\hfill N=\{b_i,i=0,1,2\dots n-1\}.\end{array}$$

(41)

Each $b_i$ boid in the group N has p position, $\overrightarrow{u}$ unit vector, and $\overrightarrow{v}$ velocity vector in (x,y,z) plane. The visibility sphere of the $b_i$ boid is described by:

$$\begin{array}{c}\hfill v_i=\{b_j\in N,\forall b_j:|b_i-b_j|<e,j=0,1,2\dots m-1\}.\end{array}$$

(42)

The flocking model of Reynolds is based on the following three principles:

Separation: avoiding collision with the adjacent agent. The general equation of the separation rule is as follows [33]:

$$\begin{array}{c}\hfill \overrightarrow{S_i}=-\sum _{\forall b_j\in v_i}^{}(b_i-b_j).\end{array}$$

(43)

By considering all boids inside one boid’s visibility sphere, $b_i$, then link each boid vector to $b_i$. The function for separation is the sum of the corresponding vectors.

Alignment: aligning velocity and position with adjacent boids. The general equation of the alignment rule is as follows [33]:

$$\begin{array}{c}\hfill \overrightarrow{m_i}=\sum _{\forall b_j\in V_i}^{}{\displaystyle \frac{\overrightarrow{V_i}}{m}},\end{array}$$

(44)

where the boid’s motion is denoted by $m_i$.

Cohesion: staying close to the flock’s center created by its adjacent boids. Suppose that the flock centering is expressed by $K_i$ for the boid $b_i$. The important step is to obtain the mass center of all visible bodies in the sphere’s visibility, $v_i$ as follows:

$$\begin{array}{c}\hfill c_i=\sum _{\forall b_j\in v_i}^{}{\displaystyle \frac{\overrightarrow{p_j}}{m}},\end{array}$$

(45)

where $c_i$ denotes the mass center of the boids. Therefore, the vector of cohesion displacement is given as follows [33]:

$$\begin{array}{c}\hfill \overrightarrow{K_i}=c_i-p_i=\sum _{\forall b_j\in v_i}^{}{\displaystyle \frac{\overrightarrow{p_j}}{m}}-p_i.\end{array}$$

(46)

In the above equations, we clarified the rules of Reynolds through simple formulas. Therefore, we utilized the formulas in [32,34,35] as follows:

• Separation: The separation force is utilized for avoiding collision between the UAVs and AUVs, which is defined as follows:

  $$\begin{array}{c}\hfill {\stackrel{\rightharpoonup }{S}}_{sep}\left(m\right)=\frac{1}{n}\sum _{j=1}^{n}k{\xi }_j\left(m\right)\frac{{\stackrel{\rightharpoonup }{p}}_s\left(m\right)-{\overline{p}}_j\left(m\right)}{∥{\stackrel{\rightharpoonup }{p}}_s\left(m\right)-{\overline{p}}_j\left(m\right)∥}\end{array}$$

  (47)

  with

  $$\begin{array}{c}\hfill {\xi }_j\left(m\right)=\left\{\begin{array}{c}\frac{1}{\sqrt{∥{\stackrel{\rightharpoonup }{p}}_s\left(m\right)-{\stackrel{\rightharpoonup }{p}}_j\left(m\right)∥}}-\frac{1}{\sqrt{d_0}},\hfill \\ if∥{\stackrel{\rightharpoonup }{p}}_s\left(m\right)-{\stackrel{\rightharpoonup }{p}}_j\left(m\right)∥\le d_0\hfill \\ 0,if∥{\stackrel{\rightharpoonup }{p}}_s\left(m\right)-{\stackrel{\rightharpoonup }{p}}_j\left(m\right)∥>d_0,\hfill \end{array}\right.\end{array}$$

  (48)

  where the position of agent s and agent j is represented by $p_s$ and $p_j$ respectively, n represents the agent’s number, k represents a positive constant, $∥{\stackrel{\rightharpoonup }{p}}_s\left(m\right)-{\overline{p}}_j\left(m\right)∥$ represents the Euclidean distance between the agent s from the agent j, the separation gain is represented by $\xi$, and $d_0$ represents the constant distance between the agents.
• Alignment: The matching average velocity of the agents is defined as follows:

  $$\begin{array}{c}\hfill {\stackrel{\rightharpoonup }{M}}_{alig}\left(m\right)=\frac{1}{n}\sum _{j=1}^n{\stackrel{\rightharpoonup }{v}}_{s,j}\left(m\right)=\frac{1}{n}\sum _{j=1}^n({\stackrel{\rightharpoonup }{v}}_j\left(m\right)-{\stackrel{\rightharpoonup }{v}}_s\left(m\right)),\end{array}$$

  (49)

  where the velocity of agent j with respect to agent s is represented by ${\stackrel{\rightharpoonup }{v}}_{s,j}\left(m\right)$ and the velocity of agent j and agent s is represented by $v_j$ and $v_s$, respectively.
• Cohesion: The cohesion force is utilized to calculate the AUVs attraction to the desired position based on the UAVs leader’s position, which is defined as follows:

  $$\begin{array}{c}\hfill {\stackrel{\rightharpoonup }{C}}_{coh}\left(m\right)=\frac{1}{n}\sum _{j=1}^n{k}_{coh}\left({\stackrel{\rightharpoonup }{p}}_d\left(m\right)-{\stackrel{\rightharpoonup }{p}}_j\left(m\right)\right),\end{array}$$

  (50)

  where the desired agent position is represented by $p_d$ and a constant gain is represented by $k_{coh}$.
• Obstacle Avoidance: The obstacle avoidance force is utilized for avoiding collision between the obstacles and the agents (UAVs and AUVs), which is defined as follows [34]:

  $$\begin{array}{c}\hfill {\stackrel{\rightharpoonup }{S}}_{obs}\left(m\right)=\frac{1}{n}\sum _{j=1}^n{k}_o{\xi }_o\left(m\right)\frac{{\stackrel{\rightharpoonup }{p}}_j\left(m\right)-{\overline{p}}_o\left(m\right)}{∥{\stackrel{\rightharpoonup }{p}}_j\left(m\right)-{\overline{p}}_o\left(m\right)∥},\end{array}$$

  (51)

  with

  $$\begin{array}{c}\hfill {\xi }_o\left(m\right)=\left\{\begin{array}{c}{r}_0-∥{\stackrel{\rightharpoonup }{p}}_j\left(m\right)-{\stackrel{\rightharpoonup }{p}}_o\left(m\right)∥,\hfill \\ if∥{\stackrel{\rightharpoonup }{p}}_j\left(m\right)-{\stackrel{\rightharpoonup }{p}}_o\left(m\right)∥\le r_0\hfill \\ 0,if∥{\stackrel{\rightharpoonup }{p}}_j\left(m\right)-{\stackrel{\rightharpoonup }{p}}_o\left(m\right)∥>r_0,\hfill \end{array}\right.\end{array}$$

  (52)

  where the position of agent j and obstacle o are represented by $p_j$ and $p_o$, respectively, n represents the agent’s number, $k_o$ represents a positive constant, $∥{\stackrel{\rightharpoonup }{p}}_j\left(m\right)-{\overline{p}}_o\left(m\right)∥$ represents the Euclidean distance between agent j and obstacle o, the separation gain is represented by ${\xi }_o$, and $r_0$ represents the constant distance between agent j and the obstacle o.

The total control law is a combination formula of all the forces, which will be given as follows:

$$\begin{array}{c}\hfill U_r=S.S_{sep}+M.M_{alig}+C.C_{coh}+S_o.S_{obs},\end{array}$$

(53)

where the control law $U_r$ is affected by positive constants S, M, C, and $S_o$.

4. Simulation Results

4.1. Formation Control Based on L1 Adaptive Control with Actuator Time Delay

The L1 controller with actuator time delay, DDS middleware, and flocking technique built above have been introduced. The underwater and quadrotor model parameters are selected as in [25,26].

Table 2. The quadrotor model parameters.

Properties	Unit	Value
Mass	m	0.52 kg
Gravity Acceleration	g	9.8 m/s${}^2$
Drag’s Translational	$k_t$	$0.95$
Drag’s Rotational	$k_r$	$0.105$
Ratio of Drag & Thrust	d	7.5 × 10${}^{-7}$ kg × m${}^2$
Inertia of x-axis	$I_x$	0.0069 kg × m${}^2$
Inertia of y-axis	$I_y$	0.0069 kg × m${}^2$
Inertia of z-axis	$I_z$	0.0129 kg × m${}^2$
Arm Length	L	0.205 m
Propeller Inertia	$I_R$	$3.36\times {10}^{-5}$ kg × m${}^2$

list the parameters for the quadrotor.

Table 3. MARES parameters.

Properties	Value
Length	1.5 m
Diameter	20 cm
Weight in air	32 kg
Depth rating	100 m
Propulsion	2 horizontal + 2 vertical thrusters
Horizontal velocity	0–1.5 m/s, variable
Energy	Li-Ion batteries, 600 Wh
Autonomy/Range	about 10 h/40 km

,

Table 4. MARES inertia.

Properties	Value [$\mathbf{kg}\times \mathit{m}$${}^2$]
$I_{xx}$	$1.55\times {10}^{-1}$
$I_{yy}$	$4.73$
$I_{zz}$	$4.73$

and

Table 5. MARES added masses.

Properties	Value	Unit
$X_{\dot{u}}$	−1.74	kg
$Y_{\dot{v}}$	$4.28\times 10$	kg
$Z_{\dot{w}}$	−4.12 × 10	kg
$K_{\dot{p}}$	$-8.61\times {10}^{-3}$	kg × m${}^2$
$M_{\dot{q}}$	−6.07	kg × m${}^2$
$N_{\dot{r}}$	−6.40	kg × m${}^2$
$X_{\dot{q}}$	$-3.05\times {10}^{-2}$	kg × m
$Y_{\dot{p}}$	$3.05\times {10}^{-2}$	kg × m
$K_{\dot{v}}$	$3.05\times {10}^{-2}$	kg × m
$M_{\dot{u}}$	$-3.05\times {10}^{-2}$	kg × m
$Y_{\dot{r}}$	$1.13\times {10}^{-1}$	kg × m
$Z_{\dot{q}}$	$-1.23\times {10}^{-1}$	kg × m
$M_{\dot{w}}$	$-1.23\times {10}^{-1}$	kg × m
$N_{\dot{v}}$	$1.13\times {10}^{-1}$	kg × m

list the MARES parameters, inertia and added masses respectively.

With the L1 control with actuator time delay parameters chosen as ${\gamma }_1={10}^4,{\gamma }_2={10}^4,{\gamma }_3={10}^4,k={10}^3$, and $\tau =0.1$. Reynolds’ rules parameters chosen as $k=k_{coh}=10,S=1,S_o=1,M=1$, and $C=5$.

For formation control, the DDS middleware communication between the agents is shown in Figure 1 by considering some policies of QoS in

Table 6. Quality of services (QoS) policies.

QoS Policies	QoS Value
	Publisher/Subscriber
Durability	Volatile
Reliability	Reliable
History	Keep All
Ownership	Shared

that must be fulfilled before the interaction between the subscribers and the publishers. The QoS illustrated in Figure 1 will be adopted during all the tests shown in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7.

Two cases have been included for the formation control results based on the modified L1 controller.

• AUVs-UAVs formation control without disturbance and uncertain parameters. Figure 2 illustrates a group of AUVs in a triangular formation tracking their UAVs leader and coordinating themselves around it in 2D space. At the same time, Figure 3 depicts the 3D space of the set of three AUVs without disturbance and uncertainty in the inertia matrix.
• Formation of AUVs-UAVs with disturbance and uncertainty in parameters. In Figure 2 and Figure 3, the disturbance and parameter uncertainties (water density and the inertia matrix) were not taken into account. Figure 4 and Figure 5 show the 2D space and the 3D space of the group of three AUVs with disturbance and uncertainty in water density and inertia matrix. Figure 6 and Figure 7 illustrate the formation control of three AUVs in triangular form and one UAV avoiding some obstacles during their movement in 2D and 3D space. One can observe that during this maneuver, the controller still maintains the formation.

4.2. Containment Control Based on L1 Adaptive Control with Actuator Time Delay

For containment control, the communication between the agents through the DDS middleware is shown in Figure 8 by considering the same policies of QoS in Table 6. This communication structure is maintained during all the tests shown in Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14.

Two cases have been included for containment control results based on the modified L1 controller.

• AUVs–UAVs containment without disturbance and uncertain parameters: The set of five AUVs can be constituted in a certain topology inside their UAVs leaders in 2D space and 3D space without disturbance and uncertainty in the inertia matrix, which is indicated in Figure 9 and Figure 10.
• AUVs–UAVs containment with disturbance and uncertain parameters: The disturbance and the parameter uncertainty were not considered in Figure 9 and Figure 10, but now it is considered. Figure 11 and Figure 12 show the 2D space and 3D space of the set of five AUVs with disturbance and uncertainty in water density and the inertia matrix. Figure 13 and Figure 14 illustrate the containment control of five AUVs and four UAVs avoiding some obstacles during their movement in 2D and 3D space. One can see that during this maneuver the controller maintain the containment and the formation of the followers within their leaders.

5. Results Discussion and Comparative Analysis

The formation of heterogeneous AUVs–UAVs without disturbance and uncertain parameters (water density and the inertia matrix) based on the modified $L_1$ controller in 2D space and 3D space is shown in Figure 2 and Figure 3, respectively. All the AUVs are following their UAV leader and maintaining their triangle formation during tracking. Information is exchanged between all the AUVs and their UAVs leader through the DDS middleware as in Figure 1.

The formation of heterogeneous AUVs–UAVs with disturbance and uncertain parameters (water density and the inertia matrix) based on the modified L1 controller in 2D space and 3D space as shown in Figure 4 and Figure 5, respectively. In the presence of these disturbances, the fleet of AUVs maintain their triangle formation with the leader in the center. The performance of the system is not influenced by these disturbances and the large variation in the water density. The communication between all the AUVs and their leader is achieved through the DDS middleware as in Figure 1.

Figure 6 and Figure 7 show how the three AUVs and one UAV avoided obstacles during their movement. The formation of the AUVs in triangle form with the leader at the center has been maintained during obstacle avoidance in both 2D and 3D space.

The containment of heterogeneous AUVs–UAVs without disturbance and uncertain parameters (water density and the inertia matrix) based on the modified L1 controller in 2D space and 3D space is illustrated in Figure 9 and Figure 10, respectively. All the AUVs have been localized in formation inside the space delimited by the four UAVs leaders. The exchange of information between all of the AUVs and their UAVs leaders has been conducted through the DDS middleware.

The containment of heterogeneous AUVs–UAVs with disturbance and uncertain parameters (water density and the inertia matrix) based on the modified L1 controller in 2D space and 3D space is shown in Figure 11 and Figure 12, respectively. The followers are tracking the leaders and are remaining within the space delimited by these leaders. The performance of the system has not been influenced by these disturbances and the large variation in the water density. The star formation of the followers and their containment within the four leaders has been maintained. The communication between all of the AUVs themselves and their UAVs leaders has been done through the DDS middleware with QoS adapted to the leaders and followers as in Figure 8. Figure 13 and Figure 14 show how the five AUVs and four UAVs avoided the obstacles during their movement while maintaining their containment in 2D and 3D space. The obstacle has been avoided by the fleet following their leaders.

The navigation and control approach provided an excellent robustness to disturbances and also avoided obstacles while maintaining the formation or the containment. The DDS middleware provided the right environment for communication and data exchange with the appropriate QoS.

6. Conclusions

A novel design for heterogeneous formation and containment control has been addressed throughout this paper. The heterogeneous formation and containment control is proposed based on an L1 controller with an actuator time delay, the DDS middleware, and a flocking approach. The L1 controller with actuator time delay is utilized to stabilize the AUVs’ and UAVs’ dynamic models. Furthermore, the flocking technique has been utilized to manage the AUVs’ positions and keep them in their desired paths in accordance with their UAVs leaders as well as to avoid the obstacles during the agents’ movement. Moreover, the DDS middleware is used for exchanging information between all of the AUVs themselves and their UAVs leaders and for solving the interoperability issue between the UAVs and the AUVs. The modified L1 adaptive controller shows high performance. The simulation results indicate that the modified L1 adaptive controller and the DDS middleware enhanced the overall performance.