Robust Neural Network Consensus for Multiagent UASs Based on Weights’ Estimation Error

Abstract

We propose a neural network consensus strategy to solve the leader–follower problem for multiple-rotorcraft unmanned aircraft systems (UASs), where the goal of this work was to improve the learning based on a set of auxiliary variables and first-order filters to obtain the estimation error of the neural weights and to introduce this error information in the update laws. The stability proof was conducted based on Lyapunov’s theory, where we concluded that the formation errors and neural weights’ estimation error were uniformly ultimately bounded. A set of simulation results were conducted in the Gazebo environment to show the efficacy of the novel update laws for the altitude and translational dynamics of a group of UASs. The results showed the benefits and insights into the coordinated control for multiagent systems that considered the weights’ error information compared with the consensus strategy based on classical $\sigma$-modification. A comparative study with the performance index ITAE and ITSE showed that the tracking error was reduced by around 45%.

1. Introduction

In the last couple of decades, exhaustive research and development began in the field of unmanned aerial vehicles in search and rescue, surveillance, monitoring and mapping the environment, and delivery and transportation tasks [1,2]. Moreover, it is widely recognized both in nature and in robotic systems that the development of a task or mission in a cooperative way offers different advantages, among which are a reduction in time, a robustness or tolerance to failures since a member that presents problems can easily be replaced by some other agent [3,4]. In this sense, consensus methods for multiagent systems (MAS) have a extensive use in the field of intelligent autonomous vehicles, such as unmanned aircraft systems vehicles (UASs), ground vehicles, and sea applications. In [5], for UAV multiagent systems with external disturbances, a finite-time distributed-formation-tracking strategy was used. The authors in [6] applied in sliding mode a consensus algorithm to track the leader in multi-UAS systems and maintain the desired formation at the same time.

With the aim of solving the problem of formation control, more attention was given to the consensus control of multiagent systems, because of their extensive application in several areas [7,8]. Mainly, consensus control is divided into two approaches [9,10], the leaderless and the leader–follower consensus. Firstly, the consensus protocols were designed for linear systems, ranging from first-order integrator systems to higher-order integrator systems [11,12]. Because of the applicability of the consensus algorithm, further developments were conducted. The consensus strategy has been developed for nonlinear systems [10,13], with time-delay [14,15], with the input and actuator saturated [11,16], with sliding-mode techniques [17,18], as well as regarding algorithms dealing with external uncertainties within systems with unknown dynamics [19]. With the development of computational resources, strategies focused on machine learning have been proposed, where, based on the observations of both the states in multiagent systems and the environment, they have to improve both the performance and communication within the consensus algorithm [20,21].

An extended survey of consensus algorithms can be founded in [22,23].

In recent years, in order to solve the problem of lumped disturbances, i.e., unmodelled dynamics and external disturbances, several neural network nonlinear consensus approaches have been developed. Moreover, these protocols are efficient when the model parameters change while the system is working, as is the case with UAVs performing logistics tasks [24]. In this sense, further research on consensus protocols based on neural networks has been conducted in [25,26]. For nonlinear multiagent systems with uncertainties, a leader–follower consensus with state and output feedback was presented in [27], where the uncertain dynamics were approximated through neural networks. An adaptive neural consensus strategy for nonaffine nonlinear multiagent systems with uncertainties, unknown control directions, and subject to switching topologies was developed in [28], where a radial basis neural network was used to approximate the unknown dynamics. A neuronal dynamic surface control based on a predictor was developed for nonlinear systems subject to uncertain dynamics in [29]; in contrast to existing methodologies where the neural weights are updated with information from the tracking errors, a predictor was used and the prediction errors were employed in the adaptation law for the neural weights. In [30], a robust consensus strategy was developed for high-order nonlinear multiagent systems subject to unknown dynamics, affected by unknown actuator failures and unknown control gains. In this approach, a backstepping was combined with a neural network to ensure that the tracking error was restricted to a small region. In [31], an output-feedback formation tracking control strategy with modeling uncertainties subject to communication constraints was developed, where a radial basis neural network was used to identify the unknown dynamic.

In the aforementioned works, a major limitation is that the only signal available that reflects the difference between estimated parameters and the real parameters is the tracking or prediction error signal. To increase the robustness of the adaptive laws, the neural network consensus protocols were improved by introducing variants in the adaptive laws for the neural weights, such as the modification $-\sigma$ [32] and the modification $-e$ [33]. However, these algorithms induced a longer convergence time and affected the response of the tracking error. As stated in [34,35] for adaptive control systems, a fast and accurate neural network weight estimation convergence is valuable for the stability and robustness properties in closed-loop systems.

Since the optimal weights are unknown, the difference between the estimated weights and the optimal weights is not available for measurement. Based on the previous observations, an adaptive estimation strategy for unknown parameters was first introduced in [36], where, unlike conventional adaptation mechanisms, the adaptation of parameters included information on the estimation error of the parameter to increase the convergence of the tracking errors.

Motivated by the use of the parameter estimation error in the development of adaptive laws, further research was conducted by the scientific community. The development of robust adaptive controllers for the online estimation of the mass parameter for robotic systems was reported in [37]. In [38,39], an adaptive parameter estimation strategy was developed to estimate parameters in vehicle systems, mainly the road gradient and the mass of the vehicle. Moreover, an adaptive control strategy was developed for vehicles’ active suspensions subject to unknown nonlinearities in [40]; in that approach, the authors proposed an adaptive law based on the parameter estimation error. This novel parameter estimation strategy was used in aerial systems as well. In [41], the inertia and mass parameters of a quadrotor aerial vehicle were estimated by introducing the parameter estimation error. Furthermore, following the same methodology, the inertia and mass parameters of a six-DOF spacecraft were estimated in [42]. Recently, in addition to being used to estimate unknown parameters, this parameter estimation methodology was employed to estimate the optimal weights for systems based on a neural network [43]. In [44], an adaptive control strategy for robot manipulators was developed by combining a neural controller and a robust controller. The proposed scheme guaranteed that the estimated weights converged in finite time to the optimal weights. Furthermore, for robotic systems, in [45], a sliding mode controller enhanced with a neural network was developed. In that approach, an adaptive law for the neural weights was proposed, which was based on the estimation error. For a dual-arm robot system with dynamic uncertainties, an adaptive command-filtered control strategy was designed in [46]. Moreover, a radial basis neural network was used to approximate the uncertainties of the system, where the information of the estimated error was used in the adaptive law for the neural weights.

Main Contributions

We propose a robust, adaptive-consensus cooperative strategy for a group of quadrotor aerial vehicles with lumped disturbances (i.e., partially unknown nonlinearities and external disturbances). In order to compensate for the lumped disturbances, a radial basis neural network is introduced in the consensus protocol. To calculate the information of the weights’ estimation error, a set of first-order intermediate variables are computed. Finally, the adaptive laws for the neural network weights incorporate the information of the estimation error. The main contributions of this research work are listed below:

• The development of a leader–follower consensus algorithm enhanced with a radial basis neural network for multiagent quadrotor unmanned aircraft systems subject to lumped disturbances is presented.
• The algorithm used to update the weights of the neural network is based on a first-order filter and auxiliary matrices to obtain the weight error information and use it in the adaptation law. To the best our knowledge, the consensus algorithms with neural network compensation reported until now have not taken into account the weight error information in the adaptation law.
• A stability proof based on Lyapunov’s theory is developed for the consensus algorithm and neural network compensation developed, where the uniformly ultimate boundedness of the (i) formation errors and (ii) neural weights’ estimation errors, is guaranteed.
• Simulation experiments for multiagent UAV systems are performed over the ROS platform [47] and ROS-based Gazebo environment.

The document is organized as follows: Some preliminaries and the dynamic model for the UAV are reviewed in Section 2. The consensus control strategy is introduced in Section 3. Section 4 provides the simulation results with examples in Gazebo, and Section 5 concludes this paper.

2. Preliminaries and Dynamic Model of the UAV

In this section, we introduce the graph theory necessary for the development of the consensus control strategy. After that, the general mathematical model for the quadrotor vehicle is presented. However, the quadrotor aerial vehicle simulated in Gazebo has an internal controller for the rotational dynamics. Simulations are conducted to identify the dynamics of the internal controller and finally a mathematical model that includes the internal controller dynamics is presented.

2.1. Graph Theory

Formally, the agents and a leader in a multiagent system are simply nodes, denoted by N, and the connectivity/structure among the agents is define by a communication graph $\mathcal{G}=\{⪰,ϝ\}$ where $⪰=\{n_1,\cdots ,n_N\}$ is a node set and $ϝ=\{(n_i,n_j)\in ⪰\times ⪰$ is an edge set composed by the elements $(n_i,n_j)$ that represent the communication link with i as the origin node and j as the final node. The so-called adjacency matrix is defined as $\mathcal{A}=\left(a_{ij}\right)\in {\mathbb{R}}^{N\times N}$, where the element $a_{ij}$ of $\mathcal{A}$ is selected as $a_{ij}=1$ if node $n_i$ is connected to node $n_j$ and $a_{ij}=0$ if is disconnected. A very important and practical feature is the degree of each node, which is simply the number of nodes that it is connected to; the in-degree matrix is $\mathcal{D}=diag\left(d_i\right)\in {\mathbb{R}}^{N\times N}$, it provides a single value of each node, with $d_i={\sum }_{j\in {\mathcal{N}}_i}{a}_{ij}$. It is also used for the computation of the most important graph operator, the graph’s Laplacian $\mathcal{L}=\left(l_{ij}\right)\in {\mathbb{R}}^{N\times N}$, defined as $\mathcal{L}=\mathcal{D}-\mathcal{A}$. If $a_{ij}=a_{ji}$, for $i,j=1,\cdots ,N$, then the communication graph $\mathcal{G}$ is defined as an undirected graph. If there is a path connecting every pair of nodes, then an undirected graph is connected. Furthermore, the leader’s adjacency matrix is defined as ${\mathcal{A}}_0=\mathrm{diag}\left(a_{i0}\right)\in {\mathbb{R}}^{N\times N}$, where $a_{i0}>0$ implies that agent ith has access to the state of the leader agent and $a_{i0}=0$ in another case. If the undirected graph is connected and at least one of diagonal entries of the leader’s adjacency matrix is not zero, then the matrix $\mathcal{H}$ defined as $\mathcal{H}=\mathcal{L}+{\mathcal{A}}_0\in {\mathbb{R}}^{N\times N}$ is symmetric and positive definite [48].

2.2. Dynamic Model of the Quadrotor UAV

A brief description of the mathematical model is shown below [49,50]. In Figure 1, we observe a three-dimensional UAV, where ${\Gamma }_B=\left[X_B,Y_B,Z_B\right]$ represents the body-fixed frame located in the center of gravity of the vehicle. The inertial frame is represented by ${\Gamma }_I=\left[X_I,Y_I,Z_I\right]$, and the position vector for the aerial vehicle in the inertial frame is denoted as $\xi =(x,y,z)\in {\Gamma }_I$. In this work, the rotation matrix used is given by

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

(1)

where $\mathbf{R}\left(\eta \right):{\Gamma }_I\to {\Gamma }_B$ and $\mathbf{R}\left(\eta \right)\in SO\left(3\right)$. $\psi$, $\theta$, and $\varphi$ are referred to as the Euler angles representing yaw, pitch, and roll, respectively. Moreover $s\left(·\right)$ and $c\left(·\right)$ stand for $sin\left(·\right)$ and $cos\left(·\right)$, respectively.

Let $\dot{\xi }$ and $\dot{\eta }$ be defined as the linear and angular velocities of the three-dimensional vehicle expressed in the inertial frame. Moreover, the linear and angular velocities represented in the body-fixed frame are denoted by ${\dot{\xi }}_b$ and ${\dot{\eta }}_b$, respectively. Then, $\eta ={\left[\varphi \theta \psi \right]}^T\in {\Gamma }_I$ represents the vector of Euler angles, and ${\dot{\eta }}_b={\left[pqr\right]}^T$, as the body angular velocities are related to $\dot{\eta }=\mathbf{W}{\dot{\eta }}_b$, where the matrix W is defined as

$$\mathbf{W}=\left[\begin{array}{ccc}1& s\varphi tan\theta & c\varphi tan\theta \\ 0& c\varphi & -s\varphi \\ 0& \frac{s\varphi }{c\theta }& \frac{c\varphi }{c\theta }\end{array}\right]$$

(2)

Using the Newton–Euler equations, we can represent the equations of motion of the UAV in the body frame as follows

$$\begin{array}{cc}\hfill m\ddot{\xi }& =\mathbf{R}\left(\eta \right)v_{z}F-mg{v}_z\hfill \end{array}$$

(3a)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{J}\left(\eta \right)\ddot{\eta }& =\tau -\mathbf{C}(\eta ,\dot{\eta })\dot{\eta }\hfill \end{array}\end{array}$$

(3b)

Here, $\mathbf{J}\left(\eta \right)={\mathbf{W}}^{\top }\mathbf{I}\mathbf{W}$ acts as the inertia matrix, where $\mathbf{I}=diag(I_x,I_y,I_z)$ and $\mathbf{C}(\eta ,\dot{\eta })\dot{\eta }=\left(\dot{\mathbf{J}}-\frac{\partial \left({\dot{\eta }}_b^{\top }\mathbf{J}\right)}{2\partial \eta }\right){\dot{\eta }}_b$ represents the Coriolis term, $v_z={\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^{\top }$ is a unitary vector along the z-axis in the inertial reference frame, F is the thrust force applied to the aerial vehicle, and $\tau \in {\mathbb{R}}^3$ denotes the roll, pitch, and yaw moments. The multirotor UAV model given previously is a highly coupled nonlinear system. Under some assumptions, it is possible to redefine the above model. Now, to simplify the moments $\tau$, the auxiliary vector $\tilde{\tau }$ is defined as follows

$$\tilde{\tau }=\left[\begin{array}{c}{u}_2/I_x\\ u_3/I_y\\ u_4/I_z\end{array}\right]={\mathbf{I}}^{-1}{\mathbf{W}}^{-1}\left(-\mathbf{I}\dot{\mathbf{W}}\dot{\eta }-\mathbf{W}\dot{\eta }\times \mathbf{I}\mathbf{W}\dot{\eta }+\tau \right)$$

(4)

Using (3) and (4), the UAV dynamical model can be represented by

$$\begin{array}{cc}\hfill \ddot{x}& =\frac{u_1}{m}(cos\psi sin\theta cos\varphi +sin\psi sin\varphi )\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill \ddot{y}& =\frac{u_1}{m}(sin\psi sin\theta cos\varphi -cos\psi sin\varphi )\hfill \end{array}$$

(5b)

$$\begin{array}{cc}\hfill \ddot{z}& =\frac{u_1}{m}(cos\theta cos\varphi )-g\hfill \end{array}$$

(5c)

$$\begin{array}{cc}\hfill \ddot{\theta }& =\frac{u_2}{I_x}\hfill \end{array}$$

(5d)

$$\begin{array}{cc}\hfill \ddot{\varphi }& =\frac{u_3}{I_y}\hfill \end{array}$$

(5e)

$$\begin{array}{cc}\hfill \ddot{\psi }& =\frac{u_4}{I_z}\hfill \end{array}$$

(5f)

where $m\in \mathbb{R}$ represents the mass of the aerial vehicle and $g\in \mathbb{R}$ is the force due to the gravity, and $u_1={\sum }_{i=1}^4{\mathbf{T}}_i$ is the thrust, which is generated by the rotors of the vehicle. For simplicity, the generalized torques are obtained as

$$\left[\begin{array}{c}{u}_2\\ u_3\\ u_4\end{array}\right]=\left[\begin{array}{c}-l{\mathbf{T}}_1-l{\mathbf{T}}_2+l{\mathbf{T}}_3+l{\mathbf{T}}_4\\ -l{\mathbf{T}}_1+l{\mathbf{T}}_2+l{\mathbf{T}}_3-l{\mathbf{T}}_4\\ -C_M{\mathbf{T}}_1+C_M{\mathbf{T}}_2-C_M{\mathbf{T}}_3+C_M{\mathbf{T}}_4\end{array}\right]$$

(6)

where $C_M$ is a constant defined from the rotor characteristics and l is the distance from the center of mass of the vehicle to the center of the rotor.

2.3. Modeling UAV Multirotor in Gazebo

The quadrotor available in Gazebo considers the presence of a low-level attitude controller provided by the PX4 autopilot. The internal controller stabilizes the angular dynamics for the roll, pitch, and yaw angles as well as the altitude in the z-axis direction. In a effort to consider the dynamics of the internal controller, we assumed that those dynamics were represented by the following differential equations

$$\begin{array}{cc}\hfill \ddot{z}=& -a_1\dot{z}-a_{2}z+a_3{u}_z\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill \ddot{\theta }=& -b_1\dot{\theta }-b_2\theta +b_3{u}_{\theta }\hfill \end{array}$$

(7b)

$$\begin{array}{cc}\hfill \ddot{\varphi }=& -c_1\dot{\varphi }-c_2\varphi +c_3{u}_{\varphi }\hfill \end{array}$$

(7c)

$$\begin{array}{cc}\hfill \ddot{\psi }=& -d_1\dot{\psi }-d_2\psi +d_3{u}_{\psi }\hfill \end{array}$$

(7d)

where $u_z$, $u_{\theta }$, $u_{\varphi }$, and $u_{\psi }$ denotes the control inputs to the quadrotor aerial vehicle. In order to identify the parameters $a_i$, $b_i$, $c_i$, and $d_i$ of Equation (7), a set of simulations was conducted in Gazebo by using step functions as control inputs and recording the input–output data. With the collected data, we used the System Identification Toolbox of MATLAB, which estimated the parameters of the dynamic model (7) by using the least squares method.

Table 1. Parameters identified for the dynamics of the internal controller.

$a_1=0.1862$	$b_1=25.353$	$c_1=22.747$	$d_1=9.769$
$a_2=0.1093$	$b_2=179.523$	$c_2=170.635$	$d_2=29.289$
$a_3=1.2737$	$b_3=177.246$	$c_3=167.685$	$d_3=29.283$

summarizes the parameters obtained.

The model was further validated by comparing the results obtained through simulations in Matlab with the data obtained from Gazebo. We can observe in Figure 2a the output from the experimental test for the z dynamics using as control input step function and the result for the transfer function. This procedure was similar within the attitude dynamics, where the control inputs of the Iris drone were related with the Euler angles as shown in Figure 2b and Figure 3a,b.

From Equations (5c) and (7a), we could obtain for the control signal $u_1$ the following

$$u_1=\frac{m\left(-a_1\dot{z}-a_{2}z+a_3{u}_z+g\right)}{cos\theta cos\varphi }$$

(8)

Substituting (8) in (5a) and (5b) and combining with the dynamics identified, we obtained a complete model that took into account the dynamics of the internal controller, and was given as

$$\begin{array}{cc}\hfill \ddot{x}& =f(z,u_z)(cos\psi tan\theta +\frac{sin\psi }{cos\theta }tan\varphi )\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill \ddot{y}& =f(z,u_z)(sin\psi tan\theta -\frac{cos\psi }{cos\theta }tan\varphi )\hfill \end{array}$$

(9b)

$$\begin{array}{cc}\hfill \ddot{z}& =-a_1\dot{z}-a_{2}z+a_3{u}_z\hfill \end{array}$$

(9c)

$$\begin{array}{cc}\hfill \ddot{\theta }& =-b_1\dot{\theta }-b_2\theta +b_3{u}_{\theta }\hfill \end{array}$$

(9d)

$$\begin{array}{cc}\hfill \ddot{\varphi }& =-c_1\dot{\varphi }-c_2\varphi +c_3{u}_{\varphi }\hfill \end{array}$$

(9e)

$$\begin{array}{cc}\hfill \ddot{\psi }& =-d_1\dot{\psi }-d_2\psi +d_3{u}_{\psi }\hfill \end{array}$$

(9f)

where $f(z,u_z)=(-a_1\dot{z}-a_{2}z+a_3{u}_z+g)$.

3. Control Strategy Development

In this section, we introduce the development of the consensus cooperative strategy for the multiagent system. In view of the underactuated property of the quadrotor vehicle, we rewrite the mathematical model of the quadrotor vehicle presented in Equation (9). After that, the consensus strategy for the translational dynamics of the team of quadrotor aerial vehicles is developed.

3.1. Mathematical Model for the Consensus Development

Before the development of the consensus control strategy for the multiagent quadrotor system, we rewrite Equation (9) to obtain a general equation to be used in the development of the consensus control strategy. From (9a) and (9b), we can define the virtual control inputs $u_x$ and $u_y$ as follows

$$\left[\begin{array}{c}{u}_x\\ u_y\end{array}\right]=(-a_1\dot{z}-a_{2}z+a_3{u}_z+g)\left[\begin{array}{cc}tan\theta & \frac{tan\varphi }{cos\theta }\\ -\frac{tan\varphi }{cos\theta }& tan\theta \end{array}\right]\left[\begin{array}{c}cos\psi \\ sin\psi \end{array}\right]$$

(10)

where $u_x$ and $u_y$ can be translated to the desired roll and pitch angles. Therefore, in this matrix representation, we can obtain the desired references angles ${\varphi }_d$ and ${\theta }_d$ that enable the system to follow the desired positions $x_d$ and $y_d$; the desired angles are given as

$$\begin{array}{cc}\hfill {\varphi }_d& =arctan\left[\frac{(u_{x}cos\psi +u_{y}sin\psi )}{f(z,u_z)}\right]\hfill \end{array}$$

(11a)

$$\begin{array}{cc}\hfill {\theta }_d& =arctan\left[-\frac{(u_{x}sin\psi +u_{y}cos\psi )cos{\varphi }_d}{f(z,u_z)}\right]\hfill \end{array}$$

(11b)

We assume that the Euler angles (roll, pitch, and yaw) are bounded in the following ranges: $(-\pi /2<\theta <\pi /2)$, $(-\pi /2<\varphi <\pi /2)$ and $(-\pi <\psi <\pi )$, respectively.

With the definitions of $u_x$ and $u_y$, the dynamics of translation for each quadrotor vehicle are given by

$$\begin{array}{cc}\hfill \ddot{x}& =u_x\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill \ddot{y}& =u_y\hfill \end{array}$$

(12b)

$$\begin{array}{cc}\hfill \ddot{z}& =-a_1\dot{z}-a_{2}z+a_3{u}_z\hfill \end{array}$$

(12c)

The control inputs $u_x$, $u_y$, and $u_z$ are computed for each vehicle through the consensus control strategy developed in Section 3.2. Once $u_x$ and $u_y$ are calculated, the desired angles ${\varphi }_d$ and ${\theta }_d$ are obtained by using (11). Moreover, we assume that the yaw reference angle is set to zero $\left({\psi }_d=0\right)$. Finally, the desired angles ${\varphi }_d$, ${\theta }_d$ and ${\psi }_d$ are the inputs for the internal controller of each aerial vehicle in the ROS environment.

3.2. Consensus Control Strategy

Now, the consensus cooperative strategy enhanced by a radial basis neural network is developed. First, it is developed for the altitude subsystem given in Equation (12c). For the horizontal plane, given by Equations (12a) and (12b), it is worth mentioning that the development of the cooperative strategy follows the same structure as that for the altitude subsystem.

3.2.1. Consensus Control Strategy for Altitude Dynamics

The control strategy designed assumes that the aforementioned low-level altitude controller has as its input a velocity command given by $u_z$. Therefore, we need to define the altitude controller such that it converges to the altitude desired and it has the capability to counteract the unknown disturbances. For this reason, we need to begin with the vertical dynamics given in (12c), which can be rewritten for each agent as

$${\dot{z}}_i=A_z{z}_i+B_z(u_{z_i}+d_{z_i}),i=1,\cdots ,N,$$

(13)

where $z_i={\left[z_{1_i}{z}_{2_i}\right]}^{\top }$, $A_z=\left[\begin{array}{cc}0& 1\\ -a_2& -a_1\end{array}\right]$, $B_z=\left[\begin{array}{c}0\\ a_3\end{array}\right]$, and $d_{z_i}$ denotes external uncertainties. Moreover, we introduce the leader’s dynamics as follows

$${\dot{z}}_0=A_z{z}_0+B_{z}r(z_0,t)$$

(14)

where $r(z_0,t)$ is a bounded input for the leader agent. The following assumption is stated in order to deal with the external disturbances $d_{z_i}$

Assumption 1.

The matched uncertainty $d_{iz}$ for the altitude dynamics defined in (13) is approximated by a radial basis neural network as

$$d_{z_i}=W_{z_i}^{*\top }{\Psi }_{z_i}\left(z_i\right)+{ϵ}_{z_i}\forall z_i\in \mathbf{D}$$

(15)

where $W_{z_i}^{*}\in {\mathbb{R}}^s$ is an unknown constant optimal weight vector. ${\Psi }_{z_i}(·):{\mathbb{R}}^n\to {\mathbb{R}}^s$ represents a vector function as ${\Psi }_{z_i}\left(z_i\right)={[{\Psi }_{1_i}\left(z_i\right),{\Psi }_{2_i}\left(z_i\right),\cdots ,{\Psi }_{s_i}\left(z_i\right)]}^T$ where ${\Psi }_i\left(x\right)$ is a radial basis activation function, which can be represented as: if ${\Psi }_i\left(x\right)=exp(-\frac{\left|\right|x-{c\left|\right|}^2}{2b^2})$, $i=1,2,\dots ,N$, where the center and width of the Gaussian functions are represented by $c\in {\mathbf{R}}^r$ and b, respectively. Moreover, the activation function satisfies $\left|\right|{\Psi }_{z_i}\left|\right|\le {\Psi }_{M_i}$ with ${\Psi }_{M_i}>0$, the approximation error ${ϵ}_{z_i}$ fulfills $|{ϵ}_{z_i}|\le {ϵ}_{z_i}^{+}$ where ${ϵ}_{z_i}^{+}>0$, and $\mathbf{D}\subset {\mathbb{R}}^n$ is a sufficiently large domain [29].

Now, we define the tracking error for each follower agent as ${\delta }_{z_i}=z_i-z_0$, whose time derivative is along (13), (14), and using the approximation (15), we obtain

$$\begin{array}{cc}\hfill {\dot{\delta }}_{z_i}& =A_z{z}_i+B_z(u_{z_i}+W_{z_i}^{*\top }{\Psi }_i\left(z_i\right)+{ϵ}_{z_i})-A_z{z}_0-B_{z}r\hfill \\ \hfill & =A_z{\delta }_{z_i}+B_z(u_{z_i}+W_{z_i}^{*\top }{\Psi }_i\left(z_i\right)+{ϵ}_{z_i}-r)\hfill \end{array}$$

(16)

With the aim to develop the adaptation procedure for the neural weights ${\widehat{W}}_{z_i}$, inspired by the approach presented in for single systems, a set of first-order filters is introduced. It is worth mentioning that with this approach, we can use information about the weight error ${\tilde{W}}_{z_i}$, in contrast to classical adaptation mechanisms. To conduct this strategy, the tracking error dynamics (16) are parameterized as

$$\begin{array}{cc}\hfill {\dot{\delta }}_{z_i}& =A_z{\delta }_{z_i}+B_z{u}_{z_i}+B_z{\Psi }_i{\left(z_i\right)}^T{W}_{z_i}^{*}+B_z({ϵ}_{z_i}-r)\hfill \\ \hfill & ={\phi }_{z_i}({\delta }_{z_i},u_{z_i})+{\Phi }_{z_i}{W}_{z_i}^{*}+{\epsilon }_{z_i}\hfill \end{array}$$

(17)

where

$${\phi }_{z_i}({\delta }_{z_i},u_{z_i})=A_z{\delta }_{z_i}+B_z{u}_{z_i},{\Phi }_{z_i}=B_z{\Psi }_i{\left(z_i\right)}^T,{\epsilon }_{z_i}=B_z({ϵ}_{z_i}-r),$$

(18)

with ${\phi }_{z_i}({\delta }_{z_i},u_{z_i})\in {\mathbb{R}}^2$ and ${\Phi }_{z_i}\in {\mathbb{R}}^{2\times s}$ is defined as the “regressor” matrix.

In this sense, the filtered variables ${\delta }_{z_{i}f}$,${\phi }_{z_{i}f}$ and ${\Phi }_{z_{i}f}$ are defined as

$$\begin{array}{cc}\hfill k{\dot{\delta }}_{z_{i}f}+{\delta }_{z_{i}f}& ={\delta }_{z_i},{\delta }_{z_{i}f}\left(0\right)=0\hfill \\ \hfill k{\dot{\Phi }}_{z_{i}f}+{\Phi }_{z_{i}f}& ={\Phi }_{z_i},{\Phi }_{z_{i}f}\left(0\right)=0\hfill \\ \hfill k{\dot{\phi }}_{z_{i}f}+{\phi }_{z_{i}f}& ={\phi }_{z_i},{\phi }_{z_{i}f}\left(0\right)=0\hfill \\ \hfill k{\dot{\epsilon }}_{z_{i}f}+{\epsilon }_{z_{i}f}& ={\epsilon }_{z_i},{\epsilon }_{z_{i}f}\left(0\right)=0\hfill \end{array}$$

(19)

where $k>0$ is a filter parameter. In order to obtain from (17) and (19) that

$${\dot{\delta }}_{z_{i}f}=\frac{{\delta }_{z_i}-{\delta }_{z_{i}f}}{k}-{\phi }_{z_{i}f}={\Phi }_{z_{i}f}{W}_{z_i}^{*}+{\epsilon }_{z_{i}f}.$$

(20)

let us define two auxiliary matrices $P_{z_i}\in {\mathbb{R}}^{s\times s}$ and $Q_{z_i}\in {\mathbb{R}}^s$ as

$$\begin{array}{cc}\hfill {\dot{P}}_{z_i}& =-\mathit{l}{P}_{z_i}+{\Phi }_{z_{i}f}^T{\Phi }_{z_{i}f},P_{z_i}\left(0\right)=0\hfill \\ \hfill {\dot{Q}}_{z_i}& =-\mathit{l}{Q}_{z_i}+{\Phi }_{zf}^T[\frac{{\delta }_{z_i}-{\delta }_{z_{i}f}}{k}-{\phi }_{z_{i}f}],Q_{z_i}\left(0\right)=0\hfill \end{array}$$

(21)

where the positive scalar l is a parameter to be designed. The solution of (21) is derived as

$$\begin{array}{cc}\hfill P_{z_i}& ={\int }_0^t{e}^{-\mathit{l}(t-r)}{\Phi }_{z_{i}f}^T\left(r\right){\Phi }_{z_{i}f}\left(r\right)dr\hfill \\ \hfill Q_{z_i}& ={\int }_0^t{e}^{-\mathit{l}(t-r)}{\Phi }_{z_{i}f}^T\left(r\right)\left[\frac{{\delta }_{z_i}\left(r\right)-{\delta }_{z_{i}f}\left(r\right)}{k}-{\phi }_{z_{i}f}\left(r\right)\right]dr\hfill \end{array}$$

(22)

Lemma 1.

The auxiliary matrix $P_{z_i}$ defined in (22) is positive definite, ${\lambda }_{min}\left(P_{z_i}\right)>{\sigma }_{\Phi }>0$, if the regressor ${\Phi }_{z_i}$ presented in (18) fulfills the PE condition [45]. A function ${\Phi }_{z_i}$ satisfies the PE condition if there exist positive constants T and γ satisfying the condition ${\int }_t^{t+T}{\Phi }_{z_i}{\left(r\right)}^T{\Phi }_{z_i}\left(r\right)dr\ge \gamma I$, $\forall t\ge 0$.

Notice that from (20) and (22), we can rewrite $Q_{z_i}$ as

$$Q_{z_i}=P_{z_i}{W}_{z_i}^{*}-{\varrho }_{z_i}$$

(23)

where ${\varrho }_{z_i}=-{\int }_0^t{e}^{-\mathit{l}(t-r)}{\Phi }_{z_{i}f}\left(r\right){\epsilon }_{z_{i}f}\left(r\right)dr\in {\mathbb{R}}^s$. As ${\epsilon }_{z_i}$ and the regressor ${\Phi }_{z_i}$ are bounded, we obtain that the filtered variables ${\Phi }_{z_{i}f}$ and ${\epsilon }_{z_{i}f}$ are also bounded. Therefore, ${\varrho }_{z_i}$ is bounded by ${\epsilon }_{z_{i}f}^{+}>0$.

Based on the solutions of $P_{z_i}$ and $Q_{z_i}$, we define an auxiliary vector $H_{z_i}$ as:

$$H_{z_i}=P_{z_i}{\widehat{W}}_{z_i}-Q_{z_i}$$

(24)

Then, based on (22) and the parameter estimation error ${\tilde{W}}_{z_i}=W_{z_i}^{*}-{\widehat{W}}_{z_i}$, we can verify that

$$\begin{array}{cc}\hfill H_{z_i}& =P_{z_i}{\widehat{W}}_{z_i}-Q_{z_i}\hfill \\ \hfill & =P_{z_i}{\widehat{W}}_{z_i}-P_{z_i}{W}_{z_i}^{*}+{\varrho }_{z_i}\hfill \\ \hfill & =-P_{z_i}{\tilde{W}}_{z_i}+{\varrho }_{z_i}.\hfill \end{array}$$

(25)

To guarantee the stability of the tracking error ${\delta }_{z_i}$, we design the following distributed control law with augmented adaptive control term

$$u_{z_i}=c_z{K}_z{e}_{z_i}-{\widehat{W}}_{z_i}^T{\Psi }_{z_i}\left(z_i\right),$$

(26)

where ${\widehat{W}}_{z_i}$ is an estimate of the unknown weights $W_{z_i}^{*}$, the scalar $c_z$ is a coupling gain, and $e_{z_i}$ is a neighborhood error represented as

$$e_{z_i}=\sum _{j\in N}{a}_{ij}(z_i-z_j)+a_{i0}(z_i-z_0),$$

(27)

Define $K_z\in {\mathbb{R}}^{1\times 2}$ as $K_z=-B_z^T{M}_z$, where $M_z$ is the solution of the Riccati inequality defined as

$$A_z^T{M}_z+M_z{A}_z+N_z-M_z{B}_z{B}_z^T{M}_z\le 0,$$

(28)

Notice that $M_z$ and $N_z\in {\mathbb{R}}^{2\times 2}$ are positive definite matrices. If the consensus control law (26) is substituted into (16), we arrive at

$$\begin{array}{cc}\hfill {\dot{\delta }}_{z_i}& =A_z{\delta }_{z_i}+B_z(c_z{K}_z{e}_{z_i}-{\widehat{W}}_{z_i}^T{\Psi }_{z_i}\left(z_i\right)+W_{z_i}^{*\top }{\Psi }_{z_i}\left(z_i\right)+{ϵ}_{z_i}-r)\hfill \\ \hfill & =A_z{\delta }_{z_i}+B_z(c_z{K}_z{e}_{z_i}+{\tilde{W}}_{z_i}^{\top }{\Psi }_{z_i}\left(z_i\right)+{ϵ}_{z_i}-r).\hfill \end{array}$$

(29)

It should be noted from (25) that the information about the weights’ estimation error ${\tilde{W}}_{z_i}$ is contained in $H_{z_i}$. Then, the adaptation law, for the weights of the ith follower agent in the altitude dynamics is proposed as

$${\dot{\widehat{W}}}_{z_i}={\Gamma }_{z_i}({\Psi }_{z_i}\left(z_i\right)e_{z_i}^T{M}_z{B}_z-k_{z_i}{H}_{z_i}).$$

(30)

where ${\Gamma }_{z_i}$ and $k_{z_i}$ are positive constants. Then, we state the theorem for the consensus cooperative control for the altitude dynamics of the quadrotor multiagent system as follows.

Theorem 1.

Consider the altitude dynamics for the ith follower agent given in Equation (13). If the regressor ${\Phi }_{z_i}$ defined in (18) fulfills the PE condition, the cooperative control strategy $u_{z_i}$ is proposed as in Equation (26), the gain $K_z$ is obtained from solving the Riccati Equation (28), the coupling gain $c_z$ fulfills condition (58), and the neural weights are updated by the adaptive law given by Equation (30), where the auxiliary matrix $H_{z_i}$ is defined in Equation (24). Then, the altitude tracking error ${\delta }_{z_i}$ and the estimation error of the neural weights ${\tilde{W}}_{z_i}$ for the ith agent converge to a bounded region defined by (66).

Proof. 

The proof of this Theorem is presented in Section 6.    □

3.2.2. Consensus Control Strategy for Position Control in the Horizontal Plane

To obtain the consensus control strategy for the X–Y plane, it is necessary to carry out the development of a consensus strategy for each axis. In this sense, we present the development for the X axis and notice that for the Y axis the procedure is similar. Consider the dynamics for the X axis presented in Equation (12a), which is rewritten as

$$\begin{array}{cc}\hfill {\dot{x}}_i& =A_x{x}_i+B_x(u_{x_i}+d_{x_i})\hfill \end{array}$$

(31)

where $x_i={\left[x_{1_i}{x}_{2_i}\right]}^{\top }$, $A_x=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]$, $B_x=\left[\begin{array}{c}0\\ 1\end{array}\right]$, and $d_{x_i}$ denotes the external uncertainties. Moreover, the leader dynamics are given as

$${\dot{x}}_0=A_x{x}_0+B_x{r}_x\left(x_0,t\right).$$

(32)

Let us define the tracking errors for the ith agent in the X axis as ${\delta }_{x_i}=x_i-x_0$. The distributed control strategy is proposed as

$$\begin{array}{cc}\hfill u_{x_i}& =c_x{K}_x{e}_{x_i}-{\widehat{W}}_{x_i}^T{\Psi }_{x_i}\left(x_i\right),\hfill \end{array}$$

(33)

where $c_x$ is a scalar coupling gain, ${\widehat{W}}_{x_i}$ is an estimate of the unknown weights $W_{x_i}^{*}$, and the neighborhood synchronization error in the X axis is defined as

$$\begin{array}{cc}\hfill e_{x_i}& =\sum _{j\in N}{a}_{ij}(x_i-x_j)+a_{i0}(x_i-x_0)\hfill \end{array}$$

(34)

Similar to Equation (30), the updated law for the weights of the neural network is given by

$${\dot{\widehat{W}}}_{x_i}={\Gamma }_{x_i}({\Psi }_{x_i}\left(x_i\right)e_{x_i}^T{M}_x{B}_x-k_x{H}_{x_i}),$$

(35)

where we can observe that the knowledge of the estimation error of the neural weights is incorporated in the term $H_{x_i}$, which is defined as

$$H_{x_i}=P_{x_i}{\widehat{W}}_{x_i}-Q_{x_i}$$

(36)

The auxiliary filtered matrices $P_{x_i}$ and $Q_{x_i}$ are defined as in (21). The filtered variables ${\delta }_{x_{i}f}$, ${\phi }_{x_{i}f}$, and ${\Phi }_{x_{i}f}$ are defined as in Equations (18) and (19). For the purpose of guaranteeing the convergence of the tracking errors ${\delta }_{x_i}$ and the weights’ estimation error ${\tilde{W}}_{x_i}$, Theorem 1 can be applied for the consensus control strategy in the X and Y axes.

4. Simulation Results

The simulation results for the robust cooperative control strategy developed are presented in this section. First, we introduce the simulation environment, where the required steps to carry out the simulations trials in ROS and Gazebo are summarized in Algorithm 1. After that, we present two simulations scenarios. The first one is for the altitude control of a team of four quadrotor aerial vehicles, with a comparative study with a classical distributed control and adaptive distributed control. In the second one, the simulation results for the translational cooperative control in the plane X–Y are presented.

4.1. Simulation Environment

In this section, we present the implementation of the algorithm within the Gazebo robotics simulator, it allows the testing of the behavior of multirotor vehicles without the need to use physical systems, due to the complexity of the tests and the fact it requires a greater number of components due to the number of vehicles that we defined as follower agents in this simulator. The environment can be constituted from the tools and packages developed by PX4 [51], which supports the connection with an embedded computer through a standardized protocol called Mavlink. We can observe that in the following Algorithm 1, we have to set the windy environment with the gazebo plugin that introduces external disturbances, followed by the number of agents to be used, then the topics and services are created for each vehicle. The Algorithm 1 summarizes the program developed.

Algorithm 1 Algorithm to set the multiagent system tests.

4.2. Simulation Scenarios

In these simulation results, we considered the configuration of the multiagent system with $N=4$ multirotor vehicles. Each aircraft had a mass m = 1.5 kg and the communication graphs utilized in the simulations are represented in Figure 4.

For the communication graph of Figure 4, the adjacency and in-degree matrices were defined as

$$\mathcal{A}=\left[\begin{array}{cccc}0& 1& 1& 0\\ 1& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right],\mathcal{D}=\left[\begin{array}{cccc}2& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(37)

Moreover, the Laplacian and leader’s adjacency matrices were defined as

$${\mathcal{A}}_0=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],\mathcal{L}=\left[\begin{array}{cccc}2& -1& -1& 0\\ -1& 1& 0& 0\\ -1& 0& 1& 0\\ 0& 0& -1& 1\end{array}\right]$$

(38)

In order to show the performance of our proposed approach, we present a comparative analysis with two other consensus control strategies presented in the literature.

• Classical nominal distributed control. In this approach, a static control strategy is used, given by

  $$u_{z_i}=c_z{K}_z{e}_{z_i},$$

  (39)

  which uses only the coupling gain $c_z$, the feedback gain $K_z$, and the neighborhood error $e_{z_i}$ [52].
• Adaptive neural network consensus control. In this strategy, a radial basis neural network is employed to compensate the lumped disturbances. The consensus cooperative strategy and the update law for the neural weights are presented below, where the weights are updated using the conventional method [53]

  $$\begin{array}{cc}\hfill u_{z_i}& =c_z{K}_z{e}_{z_i}-{\widehat{W}}_{z_i}^T{\Psi }_{z_i}\left(z_i\right)\hfill \\ \hfill {\dot{\widehat{W}}}_{z_i}& ={\Gamma }_{z_i}({\Psi }_{z_i}\left(z_i\right)e_{z_i}^T{M}_z{B}_z-k_{z_i}{\widehat{W}}_{z_i})\hfill \end{array}$$

  (40)
  In comparison with our approach, notice that for the neural weights, the adaptation law (40) includes the estimation of the weights given by the term ${\widehat{W}}_{z_i}$.

4.2.1. Altitude Control

The mathematical model for the altitude dynamic of the ith agent can be written as follows

$${\dot{z}}_i=\left[\begin{array}{cc}0& 1\\ -0.1093& -0.1862\end{array}\right]z_i+\left[\begin{array}{c}0\\ 1.2737\end{array}\right](u_{z_i}+d_{z_i})$$

(41)

For the leader agent, we selected a bounded input given as $r\left(t\right)=6-3cos\left(\frac{t}{3.5}\right)$.

Table 2. Parameters and initial conditions used in the simulation results.

Parameter	Value	Parameter	Value
$c_z$	2.7	$k_{wz}$	0.17
$K_z$	$[-0.75-0.35]$	$z_1\left(0\right)$	${\left[3.00.0\right]}^{\top }$
$\lambda$	13	$z_2\left(0\right)$	${\left[2.70.0\right]}^{\top }$
$\Gamma$	0.13	$z_3\left(0\right)$	${\left[2.50.0\right]}^{\top }$
$k_{fz}$	0.3	$z_4\left(0\right)$	${\left[3.70.0\right]}^{\top }$

presents the parameters used for this simulation scenario, including the control gains, the gains of the robust estimator, and the gains for updating the neural weights.

A radial-based neural network was used, where the output of this neural network is given in Equation (15). The neural network proposed had 25 neurons, where the centers of the neural weights were chosen as

$$ce{n}_z=⌈\begin{array}{ccccccccccccccccccccccccc}2& 2& 2& 2& 2& 3.5& 3.5& 3.5& 3.5& 3.5& 5& 5& 5& 5& 5& 7.5& 7.5& 7.5& 7.5& 7.5& 9& 9& 9& 9& 9\\ -0.6& -0.3& 0& 0.3& 0.6& -0.6& -0.3& 0& 0.3& 0.6& -0.6& -0.3& 0& 0.3& 0.6& -0.6& -0.3& 0& 0.3& 0.6& -0.6& -0.3& 0& 0.3& 0.6\end{array}⌉$$

(42)

with $b=0.7$ as the width of the hidden neuron. Furthermore, the initial weights of the network were selected randomly.

As we can see in Figure 5, a classic controller was used that was stationary so the tracking error performed by the agents was constant at 20 and 60 s. Furthermore, in the zoomed rectangle, we can see that the vehicles oscillated more in comparison to the reference virtual agent. In the second scenario in Figure 6, we can see the effects of the adaptive term as the tracking error began to decrease. However, the update law to estimate the parameters generated oscillations, specifically in the lower part, at 40 and 80 s. The robust estimator was introduced in the third scenario and we can observe in Figure 7 how the tracking error converged and the oscillations were reduced. In this work, to verify the effectiveness of the strategy proposed, we used two performance indices based on the tracking error between each one of the follower agents and the leader agent.

The indices used were the integral-time absolute error (ITAE) and integral-time square error (ITSE). The performance values are shown in

Table 3. Table of comparison of the error-integral performance indexes.

Performance Index List
Controller	UAV	ITAE	ITSE
Classical nominal distributed control	Agent 1	0.568	0.141
	Agent 2	0.567	0.139
	Agent 3	0.575	0.143
	Agent 4	0.598	0.157
Classical nominal distributed control with adaptive term	Agent 1	0.465	0.083
	Agent 2	0.463	0.081
	Agent 3	0.467	0.085
	Agent 4	0.494	0.091
Our proposed algorithm	Agent 1	0.259	0.025
	Agent 2	0.251	0.023
	Agent 3	0.263	0.027
	Agent 4	0.267	0.029

where as can be seen, the values in each index decrease and improve in our approach compared with the other methods.

4.2.2. Translational Cooperative Control

The dynamic model in the horizontal plane formed by the X and Y axes can be written as follows:

$$\begin{array}{cc}\hfill {\dot{x}}_i& =\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]x_i+\left[\begin{array}{c}0\\ 1\end{array}\right](u_{x_i}+d_{x_i})\hfill \\ \hfill {\dot{y}}_i& =\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]y_i+\left[\begin{array}{c}0\\ 1\end{array}\right](u_{y_i}+d_{y_i})\hfill \end{array}$$

(43)

It is worth mentioning that from the consensus control signals $u_{x_i}$ and $u_{y_i}$, the desired reference angles ${\varphi }_{d_i}$ and ${\theta }_{d_i}$ were obtained and used as control inputs for the quadrotor vehicles.

In order to follow a circular trajectory, the reference for the leader agent in the X and Y axes were given as

$$\begin{array}{cc}\hfill r_x\left(t\right)& =(-3cos\left(\frac{t\pi }{20}\right)+3)-3\hfill \\ \hfill r_y\left(t\right)& =-3sin\left(\frac{\pi }{20}t\right)\hfill \end{array}$$

(44)

The number of neurons was similar to the neural network for the height dynamics and the width $b=0.7$, and the centers of the neural weights were

$$ce{n}_x=⌈\begin{array}{ccccccccccccccccccccccccc}-6& -6& -6& -6& -6& -3& -3& -3& -3& -3& 0& 0& 0& 0& 0& 3& 3& 3& 3& 3& 6& 6& 6& 6& 6\\ -0.6& -0.3& 0& 0.3& 0.6& -0.6& -0.3& 0& 0.3& 0.6& -0.6& -0.3& 0& 0.3& 0.6& -0.6& -0.3& 0& 0.3& 0.6& -0.6& -0.3& 0& 0.3& 0.6\end{array}⌉$$

(45)

The parameters and the initial conditions changed the effect, when analyzing the performance of the proposed algorithm as we can observe in

Table 4. Parameters and initial conditions used in the horizontal example.

Parameter	Value	Parameter	Value
$c_x$	30	${\Gamma }_x$	0.3
$K_x$	$[-0.2-0.01]$	${\lambda }_x$	13
$k_{f\left(x\right)}$	0.3	$k_{w\left(x\right)}$	2.9
$x_1\left(0\right)$	${\left[0.00.0\right]}^{\top }$	$x_3\left(0\right)$	${\left[0.00.0\right]}^{\top }$
$x_2\left(0\right)$	${[-9.00.0]}^{\top }$	$x_4\left(0\right)$	${[-9.00.0]}^{\top }$
$y_1\left(0\right)$	${[-3.00.0]}^{\top }$	$y_3\left(0\right)$	${\left[3.00.0\right]}^{\top }$
$y_2\left(0\right)$	${[-3.00.0]}^{\top }$	$y_4\left(0\right)$	${\left[3.00.0\right]}^{\top }$

. Note that the centroids and parameters described in Table 4 for the X-axis dynamics were the same for the Y-axis dynamics of the system.

In Figure 8, we present the results of the simulations with the conditions mentioned above. The reference of the virtual leader agent was a circular path in the X–Y plane centered on $[0,0]$ and the circular references for each vehicle maintained an offset of their centers of $\pm 6$ m. Last, in Figure 8a,b, we can observe the performance of the position of the UAVs in the X and Y axes, respectively, where we can see two vehicles following the circular trajectory in their respective offsets despite having a wind disturbance $d_{x,y}=1\frac{m}{s}$ with a positive direction on the horizontal plane X-Y. In Figure 9, we show the point of view of the trajectory of the follower agents, where we can observe how the tracking error converges; a video showing the simulations in the Gazebo environment was shared at https://youtu.be/i6Qp5MMB2Ws (accessed on 22 June 2022).

5. Conclusions

A consensus protocol with a robust neural network estimator for the unknown elements and flock motion was proposed for multi-UAV systems to reduce the tracking error subject to unknown external perturbations and unmodeled dynamics. The closed-loop stability for the whole system was verified through Lyapunov’s theory, arriving at a result of a uniformly ultimately bounded stability, guaranteeing that the estimation and tracking errors were bounded. The consensus strategy and the update law for the neuronal weights for each vehicle were obtained through the stability proof as well. In addition, the proposed consensus protocol combined with the neural network allowed us to estimate and compensate the unmodeled dynamics and external perturbations. To demonstrate the performance of the proposed method, simulations were carried out in a virtual environment, where the method was also compared against a nominal distributed controller as well as against a consensus protocol in combination with an adaptive strategy where a traditional estimation method was used and was verified with two performance indexes to measure the tracking error of the multiagent system. The simulation results showed that the update of the neuronal weights with the first-order filters and auxiliary matrices was more robust against unmodeled dynamics, parameter uncertainties, and external disturbances. Future lines of research will consider the implementation of the cooperative strategy developed with multiple vehicles to observe the performance in real time and under external disturbances. Moreover, the strategy proposed can be extended to deal with transportation tasks with multiple vehicles.

6. Stability Proof

Proof. 

In order to conduct the stability proof, we define the following variables, ${\delta }_z={[{\delta }_{z_1}^T,\cdots ,{\delta }_{z_N}^T]}^T\in {\mathbb{R}}^{2N}$, $e_z={[e_{z_1}^T,\cdots ,e_{z_N}^T]}^T\in {\mathbb{R}}^{2N}$, ${ϵ}_z={[{ϵ}_{z_1},\cdots ,{ϵ}_{z_N}]}^T\in {\mathbb{R}}^N$, ${\Psi }_z={[{\Psi }_{z_1}^T\left(x_1\right),\cdots ,{\Psi }_{z_N}^T\left(x_N\right)]}^T\in {\mathbb{R}}^{Ns}$, ${\tilde{W}}_z=diag\left({\tilde{W}}_{z_i}\right)\in {\mathbb{R}}^{Ns\times N}$, and $\underline{r_z}={\left[r_z,\cdots ,r_z\right]}^T\in {\mathbb{R}}^N$. From the previous definitions and the tracking error dynamics for the ith agent given in Equation (16), the dynamics of the global tracking error ${\delta }_z$ can be written as

$${\dot{\delta }}_z=\left(I_N\otimes A_z+c_z\mathcal{H}\otimes B_z{K}_z\right){\delta }_z+\left(I_N\otimes B_z\right)\left({\tilde{W}}_{z_i}^T{\Psi }_z\left(z\right)+{ϵ}_z-\underline{r_z}\right).$$

(46)

Now, consider the candidate Lyapunov function defined as

$$V_z=\frac{1}{2}{\delta }_z^T\left(\mathcal{H}\otimes M_z\right){\delta }_z+\frac{1}{2}\mathrm{tr}\left\{{\tilde{W}}_z^T{\Gamma }_z^{-1}{\tilde{W}}_z\right\}$$

(47)

where ${\Gamma }_z=\mathrm{diag}\left({\overline{\Gamma }}_{z_i}\right)\in {\mathbb{R}}^{Ns\times Ns}$, ${\overline{\Gamma }}_{z_i}={\Gamma }_{z_i}{I}_{x\times s}\in {\mathbb{R}}^{s\times s}$, and ${\Gamma }_{z_i}>0,i=1,\dots ,N$. The time derivative of $V_z$ along trajectories of (46) is obtained as

$$\begin{array}{cc}\hfill {\dot{V}}_z=& {\delta }_z^T\left(\mathcal{H}\otimes M_z\right)\left[(I_N\otimes A_z+c_z\mathcal{H}\otimes B_z{K}_z){\delta }_z+(I_N\otimes B_z)({\tilde{W}}_z^T{\Psi }_z\left(z\right)+{ϵ}_z-\underline{r_z})\right]\hfill \\ \hfill & +\mathrm{tr}\left\{{\tilde{W}}_z^T{\Gamma }_z^{-1}{\dot{\tilde{W}}}_z\right\}\hfill \end{array}$$

(48)

which can be rewritten as

$$\begin{array}{cc}\hfill {\dot{V}}_z=& \frac{1}{2}{\delta }_z^T\left[\mathcal{H}\otimes \left(M_z{A}_z+A_z^T{M}_z\right)-2c_z{\mathcal{H}}^2\otimes \left(M_z{B}_z{B}_z^T{M}_z\right)\right]{\delta }_z\hfill \\ & +{\delta }_z^T\left(\mathcal{H}\otimes M_z{B}_z\right)\left({\tilde{W}}_z^T{\Psi }_z\left(z\right)+{ϵ}_z-\underline{r_z}\right)+tr\left\{{\tilde{W}}^T{\Gamma }^{-1}\dot{\tilde{W}}\right\}\hfill \end{array}$$

(49)

From the trace properties, we can obtain the following

$$\begin{array}{cc}\hfill {\dot{V}}_z=& \frac{1}{2}{\delta }_z^T\left[\mathcal{H}\otimes \left(M_z{A}_z+A_z^T{M}_z\right)-2c_z{\mathcal{H}}^2\otimes \left(M_z{B}_z{B}_z^T{M}_z\right)\right]{\delta }_z+{\delta }_z^T\left(\mathcal{H}\otimes M_z{B}_z\right)\left({ϵ}_z-\underline{r_z}\right)\hfill \\ \hfill & +\mathrm{tr}\left\{{\tilde{W}}_z^T{\Psi }_z\left(z\right){\delta }_z^T(\mathcal{H}\otimes M_z{B}_z)-{\widetilde{W_z}}^T{\Gamma }_z^{-1}{\dot{\widehat{W}}}_z\right\}\hfill \end{array}$$

(50)

Notice that the adaptation law for the weights of the ith agent presented in Equation (30), can be rewritten in a global form as

$$\begin{array}{cc}\hfill {\dot{\widehat{W}}}_z=& {\Gamma }_z\left[{\Psi }_z\left(z\right){\delta }_z^T(\mathcal{H}\otimes M_z{B}_z)-k_z{H}_z\right]\hfill \end{array}$$

(51)

where $k_z=\mathrm{diag}\left({\overline{k}}_{z_i}\right)\in {\mathbb{R}}^{Ns\times Ns}$, ${\overline{k}}_{z_i}=k_{z_i}{I}_{s\times s}\in {\mathbb{R}}^{s\times s}$, and $k_{z_i}$ are positive constants, and $H_z=\mathrm{diag}\left(H_{z_i}\right)\in {\mathbb{R}}^{Ns\times N}$ with $H_{z_i}$ defined in Equation (25). Substituting the global adaptation law (51) in the derivative of the Lyapunov function (52), we obtain

$$\begin{array}{cc}\hfill {\dot{V}}_z=& \frac{1}{2}{\delta }_z^T\left[\mathcal{H}\otimes \left(M_z{A}_z+A_z^T{M}_z\right)-2c_z{\mathcal{H}}^2\otimes \left(M_z{B}_z{B}_z^T{M}_z\right)\right]{\delta }_z+{\delta }_z^T\left(\mathcal{H}\otimes M_z{B}_z\right)\left({ϵ}_z-\underline{r_z}\right)\hfill \\ \hfill & +\mathrm{tr}\left\{{\tilde{W}}_z^T{k}_z{H}_z\right\}\hfill \end{array}$$

(52)

Defining $P_z=\mathrm{diag}\left(P_{z_i}\right)\in {\mathbb{R}}^{Ns\times Ns}$, $Q_z=\mathrm{diag}\left(Q_{z_i}\right)\in {\mathbb{R}}^{Ns\times N}$, and ${\varrho }_z=\mathrm{diag}\left({\varrho }_{z_i}\right)\in {\mathbb{R}}^{Ns\times N}$, we can express $H_z$ as

$$\begin{array}{cc}\hfill H_z=& P_z{\widehat{W}}_z-Q_z\hfill \\ \hfill =& -P_z{\tilde{W}}_z+{\varrho }_z.\hfill \end{array}$$

(53)

Then, substituting $H_z$ in Equation (52), we obtain the following

$$\begin{array}{cc}\hfill {\dot{V}}_z=& \frac{1}{2}{\delta }_z^T\left[\mathcal{H}\otimes \left(M_z{A}_z+A_z^T{M}_z\right)-2c_z{\mathcal{H}}^2\otimes \left(M_z{B}_z{B}_z^T{M}_z\right)\right]{\delta }_z+{\delta }_z^T\left(\mathcal{H}\otimes M_z{B}_z\right)\left({ϵ}_z-\underline{r_z}\right)\hfill \\ \hfill & -\mathrm{tr}\left\{{\tilde{W}}_z^T{k}_z{P}_z\tilde{W}+{\tilde{W}}_z^T{k}_z{\varrho }_z\right\}\hfill \end{array}$$

(54)

From Young’s inequality, for the second term in (54), we have

$$\begin{array}{cc}\hfill {\dot{V}}_z=& \frac{1}{2}{\delta }_z^T\left[\mathcal{H}\otimes \left(M_z{A}_z+A_z^T{M}_z\right)-2c_z{\mathcal{H}}^2\otimes \left(M_z{B}_z{B}_z^T{M}_z\right)\right]{\delta }_z\hfill \\ & +\frac{1}{2}{\delta }_z^T\left(\mathcal{H}\otimes M_z{B}_z\right)\left({\mathcal{H}}^{\top }\otimes B_z^{\top }{M}_z\right){\delta }_z-\mathrm{tr}\left\{{\tilde{W}}_z^T{k}_z{P}_z\tilde{W}+{\tilde{W}}_z^T{k}_z{\varrho }_z\right\}\hfill \\ & +\frac{1}{2}{\left({ϵ}_z-\underline{r_z}\right)}^{\top }\left({ϵ}_z-\underline{r_z}\right)\hfill \\ \hfill =& \frac{1}{2}{\delta }_z^T\left[\mathcal{H}\otimes \left(M_z{A}_z+A_z^T{M}_z\right)-\left(2c_z+1\right){\mathcal{H}}^2\otimes \left(M_z{B}_z{B}_z^T{M}_z\right)\right]{\delta }_z\hfill \\ & -\mathrm{tr}\left\{{\tilde{W}}_z^T{k}_z{P}_z{\tilde{W}}_z+{\tilde{W}}_z^T{k}_z{\varrho }_z\right\}+\frac{1}{2}{\left({ϵ}_z-\underline{r_z}\right)}^{\top }\left({ϵ}_z-\underline{r_z}\right)\hfill \end{array}$$

(55)

In order to simplify, we use a unitary matrix U such that $U^T\mathcal{H}U={\Lambda }_z=diag({\lambda }_{z_1},\cdots ,$${\lambda }_{z_N})$, where ${\lambda }_{z_i}$ are the eigenvalues of matrix $\mathcal{H}$, which is positive definite. Let us introduce a state transformation given as ${\xi }_z=(U^T\otimes I_n){\delta }_z$ with ${\xi }_z={[{\xi }_{z_1}^T,\cdots ,{\xi }_{z_N}^T]}^T$. Then, introducing the state transformation in Equation (55), we obtain

$$\begin{array}{cc}\hfill {\dot{V}}_z=& \frac{1}{2}{\xi }_z^T\left[\Lambda \otimes \left(M_z{A}_z+A_z^T{M}_z\right)-\left(2c_z+1\right){\Lambda }^2\otimes \left(M_z{B}_z{B}_z^T{M}_z\right)\right]{\xi }_z\hfill \\ & -\mathrm{tr}\left\{{\tilde{W}}_z^T{k}_z{P}_z{\tilde{W}}_z+{\tilde{W}}_z^T{k}_z{\varrho }_z\right\}+\frac{1}{2}{\left({ϵ}_z-\underline{r_z}\right)}^{\top }\left({ϵ}_z-\underline{r_z}\right)\hfill \end{array}$$

(56)

which is rewritten in an expanded form as

$$\begin{array}{cc}\hfill {\dot{V}}_z=& \frac{1}{2}{\displaystyle \sum _{i=1}^N}{\lambda }_{z_i}{\xi }_{z_i}^{\top }\left[M_z{A}_z+A_z^T{M}_z-\left(2c_z+1\right){\lambda }_{z_i}{M}_z{B}_z{B}_z^T{M}_z+N_z\right]{\xi }_{z_i}-\frac{1}{2}{\displaystyle \sum _{i=1}^N}{\lambda }_{z_i}{\xi }_{z_i}^{\top }{N}_z{\xi }_{z_i}\hfill \\ \hfill & -{\displaystyle \sum _{i=1}^N}{k}_{z_i}{\tilde{W}}_{z_i}^T{P}_{z_i}{\tilde{W}}_{z_i}-{\displaystyle \sum _{i=1}^N}{k}_{z_i}{\tilde{W}}_{z_i}^T{\varrho }_{z_i}+\frac{1}{2}{\left({ϵ}_z-\underline{r_z}\right)}^{\top }\left({ϵ}_z-\underline{r_z}\right)\hfill \\ \hfill \le & \frac{1}{2}{\displaystyle \sum _{i=1}^N}{\lambda }_{z_i}{\xi }_{z_i}^{\top }\left[M_z{A}_z+A_z^T{M}_z-\left(2c_z+1\right)min\left({\lambda }_{z_i}\right)M_z{B}_z{B}_z^T{M}_z+N_z\right]{\xi }_{z_i}\hfill \\ & -\frac{1}{2}{\displaystyle \underset{i=1,\cdots ,N}{min}}\left({\lambda }_{z_i}\right){\displaystyle \sum _{i=1}^N}{\xi }_{z_i}^{\top }{N}_z{\xi }_{z_i}-{\displaystyle \sum _{i=1}^N}{k}_{z_i}{\tilde{W}}_{z_i}^T{P}_{z_i}{\tilde{W}}_{z_i}-{\displaystyle \sum _{i=1}^N}{k}_{z_i}{\tilde{W}}_{z_i}^T{\varrho }_{z_i}+{ϵ}_z^{+}\hfill \end{array}$$

(57)

where $\frac{1}{2}\Vert {{ϵ}_z-\underline{r_z}\Vert }^2\le {ϵ}_z^{+}$ with ${ϵ}_z^{+}$ a positive scalar. Let us define the coupling gain $c_z$ as

$$c_z\ge \frac{1}{2{\displaystyle \underset{i=1,\cdots ,N}{min}}\left({\lambda }_{z_i}\right)}-\frac{1}{2},$$

(58)

such that $\left(2c_z+1\right)min\left({\lambda }_{z_i}\right)\ge 1$. Since the part $\left(A_z,B_z\right)$ is stabilizable, then the existence of the solution matrix $M_z$ is guaranteed, such that

$$M_z{A}_z+A_z^T{M}_z+N_z-\left(2c_z+1\right)min\left({\lambda }_{z_i}\right)M_z{B}_z{B}_z^T{M}_z\le M_z{A}_z+A_z^T{M}_z+N_z-M_z{B}_z{B}_z^T{M}_z\le 0.$$

(59)

Then, from Equation (59), we can rewrite (57) as

$${\dot{V}}_z\le -\frac{1}{2}{\displaystyle \underset{i=1,\cdots ,N}{min}}\left({\lambda }_{z_i}\right){\displaystyle \sum _{i=1}^N}{\xi }_{z_i}^{\top }{N}_z{\xi }_{z_i}-{\displaystyle \sum _{i=1}^N}{k}_{z_i}{\tilde{W}}_{z_i}^T{P}_{z_i}{\tilde{W}}_{z_i}-{\displaystyle \sum _{i=1}^N}{k}_{z_i}{\tilde{W}}_{z_i}^T{\varrho }_{z_i}+{ϵ}_z^{+}$$

(60)

From Lemma 1, the matrices $P_{z_i}$ are positive definite. Then, defining ${\displaystyle \underset{min}{\lambda }}\left(P_{z_i}\right)$ as the minimum eigenvalue of $P_{z_i},i=1,\cdots ,N$, we obtain

$${\dot{V}}_z\le -\frac{1}{2}{\displaystyle \underset{i=1,\cdots ,N}{min}}\left({\lambda }_{z_i}\right){\displaystyle \sum _{i=1}^N}{\xi }_{z_i}^{\top }{N}_z{\xi }_{z_i}-k_z^{+}{\displaystyle \underset{min}{\lambda }}\left(P_{z_i}\right){\displaystyle \sum _{i=1}^N}{\tilde{W}}_{z_i}^T{\tilde{W}}_{z_i}-{\displaystyle \sum _{i=1}^N}{k}_{z_i}{\tilde{W}}_{z_i}^T{\varrho }_{z_i}+{ϵ}_z^{+}$$

(61)

Removing the summations and grouping the terms, we obtain

$${\dot{V}}_z\le -\frac{1}{2}{\displaystyle \underset{i=1,\cdots ,N}{min}}\left({\lambda }_{z_i}\right){\xi }_z^{\top }\left(I_N\otimes N_z\right){\xi }_z-k_z^{+}{\displaystyle \underset{min}{\lambda }}\left(P_w\right)\mathrm{tr}\left\{{\tilde{W}}_z^{\top }{\tilde{W}}_z\right\}+\left|\mathrm{tr}\left\{{\tilde{W}}_z^T{k}_z{\varrho }_z\right\}\right|+{ϵ}_z^{+},$$

(62)

From the definition of ${\xi }_z$, we obtain by the Cauchy–Schwartz inequality and (62) the following

$${\dot{V}}_z\le -\frac{1}{2}{\displaystyle \underset{i=1,\cdots ,N}{min}}\left({\lambda }_{z_i}\right){\delta }_z^{\top }\left(I_N\otimes N_z\right){\delta }_z-k_z^{+}{\displaystyle \underset{min}{\lambda }}\left(P_{z_i}\right){∥\tilde{W}∥}_F^2+k_{z_F}^{+}{\varrho }_z^{+}{∥\tilde{W}∥}_F+{ϵ}_z^{+}$$

(63)

where $k_{z_F}^{+}={\Vert k_z\Vert }_F$, ${\Vert {\varrho }_z\Vert }_F\le {\varrho }_{z_F}^{+}$ with ${\varrho }_{z_F}^{+}$ a positive constant. By Young’s inequality $ab\le \frac{a^2}{2\eta }+\frac{\eta b^2}{2}$ with $\eta >0$, we obtain

$$\begin{array}{cc}{\dot{V}}_z& \le -\frac{1}{2}{\displaystyle \underset{i=1,\cdots ,N}{min}}\left({\lambda }_{z_i}\right){\lambda }_{min}\left(N_z\right){∥{\delta }_z∥}^2-\left(k_z^{+}{\displaystyle \underset{min}{\lambda }}\left(P_{z_i}\right)-\frac{1}{2{\eta }_z}\right){∥\tilde{W}∥}_F^2+\frac{{\eta }_z{\left(k_{z_F}^{+}{\varrho }_z\right)}^2}{2}+{ϵ}_z^{+}\hfill \\ & \le -{\alpha }_z{V}_z+{\beta }_z\hfill \end{array}$$

(64)

where

$${\alpha }_z=min\left(\frac{{\lambda }_{min}\left(\mathcal{H}\right){\lambda }_{min}\left(N_z\right)}{{\lambda }_{max}\left(\mathcal{H}\right){\lambda }_{max}\left(M_z\right)},\frac{2\left(k_z^{+}{\displaystyle \underset{min}{\lambda }}\left(P_{z_i}\right)-\frac{1}{2{\eta }_z}\right)}{{\sigma }_{max}\left({\Gamma }_z^{-1}\right)}\right),{\beta }_z=\frac{{\eta }_z{\left(k_{z_F}^{+}{\varrho }_z\right)}^2}{2}+{ϵ}_z^{+}$$

(65)

Considering the extended Lyapunov theory, it is possible to conclude that ${\delta }_z$ and $\tilde{W}$ are uniformly ultimately bounded, and converge to a small compact set given by

$${\Omega }_z:=\left\{{\delta }_z,\tilde{W}|∥{\delta }_z∥\le \sqrt{\frac{2}{{\lambda }_{min}\left(\mathcal{H}\right){\lambda }_{min}\left(M_z\right)}\frac{{\beta }_z}{{\alpha }_z}},∥\tilde{W}∥\le \sqrt{\frac{2}{{\sigma }_{min}\left({\Gamma }_z^{-1}\right)}\frac{{\beta }_z}{{\alpha }_z}}\right\}$$

(66)

Therefore, the tracking errors ${\delta }_{z_i}$ for the ith agent are uniformly ultimately bounded. □