A sectorial fuzzy consensus algorithm for the formation flight of multiple quadrotor unmanned aerial vehicles

Abstract

This paper presents an algorithm based on fuzzy theory for the formation flight of the multi-quadrotors. For this purpose, the mathematical model of N-quadrotor unmanned aerial vehicles is presented using the Newton-Euler formulation. The strategy of the formation flight is based on a structure composed by a sectorial fuzzy controller and the linear systems whose state variables are the position and velocity of the i^th quadrotor. The stability analysis is described as a generalized form for N-quadrotor unmanned aerial vehicles and it is based on the Lyapunov theory. This analysis demonstrates that the closed-loop system is globally asymptotically stable so that the quadrotors unmanned aerial vehicles reach the consensus. Numerical simulation demonstrates the robustness of the proposed scheme for the formation flight even in the presence of disturbances. Finally, experimental results show the feasibility of the proposed algorithm for the formation flight of multiple unmanned aerial vehicles.

Introduction

Formation flight for multiple unmanned aerial vehicles (UAVs) has been studied in the past years due to its civilian and military applications. These applications include mapping, precision agriculture, environmental monitoring, and search and rescue. In effect, some applications require robust flights where two or more UAVs are used to execute group tasks. Some research works have used the formation flight as a way to improve the performance of the aircraft in their missions, getting benefits such as the time reduction of the task, coverage in iterative applications, and even reductions in the induced drag, allowing maximum endurance and range.^1,2 In Cao and Ren ³ a distributed and coordinated tracking based on variable structure approach for group of autonomous vehicles was presented, and the consensus tracking problem was solved in which a group of autonomous vehicles can track a virtual leader.

Gu and Hu ⁴ presented the classical Reynolds rules based on a fuzzy logic controller; specifically, with the separation component, which is considered as a gradient-like function in order to stabilize the overall formation with constrained magnitude control inputs and to represent the separation as a repulsive. In Yang et al. ⁵ a simulation of a second order dynamic model for a multi-agent system (MAS) was presented. The flocking behavior is achieved by combining the rate consensus algorithm, based on optimal control theory, plus a fuzzy logical controller, which is built as an attraction and repulsion function. The whole consensus is achieved by using network topology communication, with directed graphs, in the overall swarm. Similarly, Chang et al. ⁶ worked on a first order formation protocol with a leader-follower configuration. A locally distributed connectivity of the MAS is addressed as well as the collision-avoidance phenomena. The previous work is solved with a fuzzy separation controller, and a neural fuzzy formation controller, with the capability of online learning, with the gradient descent method. The simulation results provide better responses in comparison with conventional consensus formation and potential-based collision-avoidance algorithms.

Geometric flight protocols were studied in literature. ⁷ In this work, a leader-follower configuration is used for a group of quadrotors. To achieve a desired formation in the XY plane, a fuzzy logic controller is used as coordination protocol. The main development of this research involves the online tuning of the parameters by means of a genetic algorithm. The effectiveness of the application is only presented in simulation. Nonlinear MASs were studied in literature ⁸ where an adaptive fuzzy control structure describes the formation control of nonlinear agents by using the position errors in an artificial potential function. The adaptive controller is designed using a fuzzy logic scheme and the $H_{\infty }$ robust optimal control concept. Simulations results are used to test the performance of the application. Fuzzy logic controllers were proposed for formation control, in literature ⁹ where the attitude and coordination of quadrotors are addressed. The orientation controller is based on the Takagi-Sugeno inference fuzzy methodology, while the coordination among vehicles is performed by using a quasi e-lattice approach. Numerical simulations are presented to validate the theoretical scheme. Muslimov and Munasypov ¹⁰ presented a new three-dimensional formation flight control algorithm for fixed-wing UAVs based on non-uniform guidance vector field. Consensus-based model of living organisms motor neural network was used for generalized law of UAVs interaction, and the adaptive loop stability is provided with fuzzy Lyapunov synthesis. Simulation results are obtained in order to show the performance of the full dynamic of the low-level nonlinear models. In Kumbasar and Tekinalp ¹¹ a fuzzy logic control algorithm was proposed in order to control the formation flight of fixed wing aircraft in a leader-follower form. The algorithm consists of a control with two loop structure in which the inner loop of the aircraft is controlled using the state dependent Riccati equation method and the outer loop. The fuzzy logic control approach is used to control the formation guidance, and the performance of the proposed algorithm is demonstrated through simulation results. Li et al. ¹² proposed a fuzzy algorithm to track and maintain close relative planar spacing during close formation maneuver and to account for the aerodynamic impacts of leading UAV on the trailing UAV. Simulation results showed the effectiveness of the proposed method.

Muñoz et al. ¹³ proposed a control strategy based on distributed adaptive leader-follower consensus algorithms, for MASs affected by switching network events. Liu et al. ¹⁴ presented the consensus tracking of the leader-follower MASs via second-order Super Twisting sliding mode control and only the simulation results were presented in order to demonstrate the effectiveness of the algorithm. Yang and Fan ¹⁵ presented the distributed formation control for second-order MASs under the leader-follower control structure. The formation control problem with constant and time-varying reference states is proposed in the presence external disturbances. In this sense, a non-saturated proportional-integral controller is proposed to attenuate the constant disturbances, and a saturated distributed control algorithm is proposed to address the formation control, which is subject to bounded control inputs. To demonstrate the effectiveness of the proposed formation controllers, the numerical simulation results are presented. Wang et al. ¹⁶ proposed a new predictor-based consensus disturbance rejection method for high order MASs with Lipschitz nonlinearity and input delay. To estimate the disturbance under the delay constraint, an observer is developed for consensus control for each agent. An unconventional predictor control scheme is constructed for each agent by utilizing the estimate of the disturbance and the prediction of the relative state information. Simulation results are presented to demonstrate the performance of the proposed controller. Qin et al. ¹⁷ investigated the formation consensus controller problem for nonlinear MASs in which a virtual leader and a distributed formation control strategy based on linear extended state observer is proposed without velocity measurement of the neighboring agents. Simulation results with different scenarios are presented to demonstrate the performance of formation consensus controller. Yu et al. ¹⁸ investigated practical time-varying formation tracking problems for second order nonlinear MASs with multiple leaders using adaptive neural networks. A practical protocol using adaptive is proposed, which is constructed using only local neighboring information. The proposed control protocol processed the matched-mismatched heterogeneous nonlinearities and disturbances, the unknown control inputs of the leaders, and the stability of MAS was presented by using the Lyapunov theory. A simulation result was shown to illustrate the effectiveness of the obtained theoretical approach. In Rojo-Rodriguez et al. ¹⁹ two schemes for formation flight were proposed based on a PID and a second order sliding mode. In this paper, the two algorithms were studied as a comparison, and mainly these algorithms control the orientation and position dynamics of each agent computing its own algorithm based on local information and information from its neighbors in the MASs. The algorithms are validated in the real-time experiments in order to demonstrate the performance of the algorithms. Compared with those research works, our proposed algorithm provides a robust consensus based on a sectorial fuzzy control and distributed structure. This control operates as position and velocity sectors emulating smooth saturation functions and allows the MAS to reach the global asymptotic convergence about the origin so that the consensus is reached. The stability analysis based on Lyapunov theory is presented for the UAVs as a generalized form. In this sense, the tuning strategy is based on heuristic methods which allows the quadrotor UAVs to perform a robust formation flight. Most references about fuzzy consensus only present simulation results; nevertheless, in this paper the proposed algorithm is implemented and run in real-time experiments in order to demonstrate the effectiveness of our approach.

The main contribution of this paper focuses on a fuzzy consensus algorithm for the formation flight of quadrotor UAVs. The quadrotors are modeled using the Newton-Euler formulation which describes the translational and rotational dynamics. For the formation flight strategy, a sectorial fuzzy consensus algorithm for formation flight of multiple quadrotor UAVs is proposed and the stability analysis is presented to prove the global asymptotic convergence about the origin. Thus, the main contribution of this paper is mentioned as follows:

The fuzzy consensus strategy, based on a sectorial fuzzy control, is proposed for multiple quadrotor UAVs.

The stability analysis based on Lyapunov theory is presented for the N-quadrotor UAVs as a generalized form.

The proposed algorithms of fuzzy consensus for multiple quadrotor UAVs are implemented in real-time experiments.

The organization of the paper is as follows: section “Equations of motion for N-quadrotor UAVs” presents the equations of motion for the N-quadrotor UAVs using the Newton-Euler formulation. The section “Sectorial fuzzy consensus for formation flight” describes the sectorial fuzzy consensus protocol, as well as consensus considerations and fuzzy logic mathematical tools. The section “Numerical simulation” shows the simulation results of the proposed scheme. The section “Experimental results” describes the real-time results showing the validation of the proposed algorithms of the formation flight, and the conclusions are given in Section “Conclusions”.

Equations of motion for N-quadrotor UAVs

The mathematical model of the UAVs agents considers three reference frames. A ground-fixed inertial framed defined as ${\mathcal{I}}_i$  =  $\{x_{{\mathcal{I}}_i},y_{{\mathcal{I}}_i},z_{{\mathcal{I}}_i}\}$ , a body frame attached to the center of gravity of the aircraft as ${\mathcal{B}}_i$  =  $\{x_{{\mathcal{B}}_i},y_{{\mathcal{B}}_i},z_{{\mathcal{B}}_i}\}$ , and a wind frame that considers the aerodynamics forces of the vehicle structure, ${\mathcal{W}}_i$  =  $\{x_{{\mathcal{W}}_i},y_{{\mathcal{W}}_i},z_{{\mathcal{W}}_i}\}$ , ²⁰ see Figure 1.

Figure 1.

The quadrotor UAV.

In order to describe the equations of motion, a Newton-Euler formulation is used, considering $i=1\dots N$ agents. The model is defined as

$${\dot{\xi }}_i=V_i$$

(1)

$$m_i{\dot{V}}_i=(-T_{T_i})R_i{e}_3+m_{i}g{e}_3+D_{{\xi }_i}$$

(2)

$${\dot{R}}_i=R_i{\widehat{\Omega }}_i$$

(3)

$$J_i{\dot{\Omega }}_i=-{\Omega }_i\times J_i{\Omega }_i+{\tau }_{a_i}+D_{{\eta }_i}$$

(4)

where ${\xi }_i={(x_i,y_i,z_i)}^{\mathrm{\top }}\in {ℝ}^3$ are the position coordinates relative to the inertial frame and ${\eta }_i={({\varphi }_i,{\theta }_i,{\psi }_i)}^{\mathrm{\top }}\in {ℝ}^3$ describes the rotation coordinates for the ith UAV; this given by an orthogonal rotation matrix $R_i\in SO(3):{\mathcal{B}}_i\to {\mathcal{I}}_i$ parametrized by the Euler angles ${\varphi }_i$ roll, θ_i pitch and ψ_i yaw or heading.

In the same way, ${\Omega }_i={\left(p_i,q_i,r_i\right)}^{\mathrm{\top }}\in {ℝ}^3$ is the angular velocity vector in ${\mathcal{B}}_i,V_i={({\dot{x}}_i,{\dot{y}}_i,{\dot{z}}_i)}^{\mathrm{\top }}\in {ℝ}^3$ is the translational velocity vector in ${\mathcal{I}}_i,T_{T_i}\in {ℝ}_{>0}$ is the total thrust and ${\tau }_{a_i}\in {ℝ}^3$ are the moments due the actuators. In addition, vectors of canonical basis of ${ℝ}^3$ are considered, these are represented by e₁, e₂, and e₃. The term $m_i\in ℝ$ denotes the mass of the ith UAV, while $J_i\in {ℝ}^{3\times 3}$ contains the moments of inertia of the ith UAV. ${\widehat{\Omega }}_i$ is the skew-symmetric matrix associated with the cross product (i.e. $\widehat{a}b=a\times b\forall a,b\in {ℝ}^3$ ). $D_{{\xi }_i}={(d_{{\xi }_{i_1}},d_{{\xi }_{i_2}},d_{{\xi }_{i_3}})}^{\mathrm{\top }}\in {ℝ}^3$ and $D_{{\eta }_i}={(d_{{\eta }_{i_1}},d_{{\eta }_{i_2}},d_{{\eta }_{i_3}})}^{\mathrm{\top }}\in {ℝ}^3$ are disturbances vectors.^21–24

The total thrust $T_{T_i}$ and the actuator moments ${\tau }_{a_i}$ generated by the four rotors are describe as follows

$$\begin{array}{l}{T}_{T_i}=\sum _{k=1}^4{f}_k,\\ {\tau }_{a_i}=\left(\begin{array}{c}{\tau }_{{\varphi }_i}\\ {\tau }_{{\theta }_i}\\ {\tau }_{{\psi }_i}\end{array}\right)=\left(\begin{array}{c}d\left[\right(f_2+f_4)-(f_1+f_3\left)\right]\\ d\left[\right(f_3+f_4)-(f_1+f_2\left)\right]\\ Q_2+Q_3-Q_1-Q_4\end{array}\right)\end{array}$$

where d denotes the distance from the center of mass to the rotor axis, $f_k={\rho }_a{A}_{d_k}{\omega }_k^2{r}_k^2{c}_{T_k}$ is the lift force or thrust force of the propeller with ρ_a that is the air density, $A_{d_k}$ denotes the disk area, ω_k denotes the angular velocity of the rotor, r_k is the rotor radius and $c_{T_k}$ is the thrust coefficient. $Q_k={\rho }_a{A}_{d_k}{\omega }_k^2{r}_k^3{c}_{Q_k}$ is the reaction moment for the rotors with $c_{Q_k}$ that denotes the rotor shaft moment coefficient for $k=1,2,3,4$ .

$$\left(\begin{array}{c}{T}_{T_i}\\ {\tau }_{{\varphi }_i}\\ {\tau }_{{\theta }_i}\\ {\tau }_{{\psi }_i}\end{array}\right)=\left(\begin{array}{cccc}{\kappa }_0& {\kappa }_0& {\kappa }_0& {\kappa }_0\\ -d{\kappa }_0& d{\kappa }_0& -d{\kappa }_0& d{\kappa }_0\\ -d{\kappa }_0& -d{\kappa }_0& d{\kappa }_0& d{\kappa }_0\\ -{\kappa }_1& {\kappa }_1& {\kappa }_1& -{\kappa }_1\end{array}\right)\left(\begin{array}{c}{\omega }_1^2\\ {\omega }_2^2\\ {\omega }_3^2\\ {\omega }_4^2\end{array}\right)$$

where ${\kappa }_0={\rho }_a{A}_{d_k}{r}_k^2{c}_{T_k}$ and ${\kappa }_1={\rho }_a{A}_{d_k}{r}_k^3{c}_{Q_k}$ .

In hover flight, the main purposes of the rotors are to provide vertical lifting force in opposition to the weight of the quadrotor UAV; however, in forward flight the rotors must also provide a propulsive force P_k to overcome the drag of the vehicle for $k=1,2,3,4$ , see Figure 2. Due to the quadrotor frame, aerodynamics forces are presented in forward flight in a 3D environment. ²⁵

Figure 2.

Schematic of the quadrotor UAV in forward flight.

Gyroscopic moment. The gyroscopic moment generated by the rotation of the airframe and the four propellers is described as

$$d_{{\eta }_{\mathit{gyr}{o}_i}}={\displaystyle \sum _{k=1}^4{(-1)}^{k+1}}{I}_{r_k}[\Omega \times ({\omega }_k{e}_3)]$$

where $I_{r_k}$ is the moment of inertia of the rotor k and ω_k denotes the angular velocity of the rotor k, with k =  $1,2,3,4$ ^22,26–28

Disturbances due to wind

The aerodynamic forces produced during the flight are written as

$$D_{{\xi }_i}=\left(\begin{array}{c}{d}_{{\xi }_1}\\ d_{{\xi }_2}\\ d_{{\xi }_3}\end{array}\right)=R{W}^T\left(\begin{array}{c}D\\ Y\\ 0\end{array}\right)$$

with the rotation aerodynamic matrix $W:\mathcal{B}\to \mathcal{A}$ that transforms a force from the body frame to wind frame is described as

$$W=\left(\begin{array}{ccc}{c}_{\alpha }{c}_{\beta }& s_{\beta }& s_{\alpha }{c}_{\beta }\\ -c_{\alpha }{s}_{\beta }& c_{\beta }& -s_{\alpha }{s}_{\beta }\\ -s_{\alpha }& 0& c_{\alpha }\end{array}\right)$$

where α is the angle of attack and β are the sideslip angle. Y and D are the aerodynamic side force, and drag force, respectively. ²¹

The aerodynamic moments. The aerodynamic moments generated during the flight are written as described as

$$d_{{\eta }_{\mathit{aer}{o}_i}}=\left(\begin{array}{c}\mathcal{L}\\ \mathcal{M}\\ \mathcal{N}\end{array}\right)$$

(5)

where $\mathcal{L},\mathcal{M}$ and $\mathcal{N}$ are the aerodynamic rolling, pitching and yawing moments respectively.^20,22

Rewriting the disturbances in the time-varying vectors $D_{{\xi }_i}$ and $D_{{\eta }_i}$ results in

$$D_{{\xi }_i}=\left(\begin{array}{c}{d}_{{\xi }_1}\\ d_{{\xi }_2}\\ d_{{\xi }_3}\end{array}\right),D_{{\eta }_i}=d_{{\eta }_{\mathit{gyr}{o}_i}}+d_{{\eta }_{\mathit{aer}{o}_i}}=\left(\begin{array}{c}{d}_{{\eta }_1}\\ d_{{\eta }_2}\\ d_{{\eta }_3}\end{array}\right)$$

which are continuously differentiable.

Guidance, navigation and control for the i aerial vehicle

In order to propose a guidance navigation and control algorithm for the i aerial vehicle as a distributed strategy, equations (1) to (4) are rewritten as follows

$${\dot{\xi }}_i=V_i$$

(6)

$${\dot{V}}_i=u_{p_i}+d_{{\xi }_i}$$

(7)

$${\dot{R}}_i=R_i{\widehat{\Omega }}_i$$

(8)

$${\dot{\Omega }}_i=u_{a_i}+d_{R_i}$$

(9)

with $m_i{u}_{p_i}=m_{i}g{e}_3-T_{T_i}(R_i{e}_3),{\tau }_{a_i}=J{u}_{a_i},d_{{\xi }_i}=\frac{D_{{\xi }_i}}{m_i},d_{R_i}=J_i^{-1}\left[-{\Omega }_i\times J_i{\Omega }_i+D_{{\eta }_i}\right]$ . In this sense, $u_{p_i}\in {ℝ}^3$ and $u_{a_i}\in {ℝ}^3$ are virtual control inputs for the position and orientation dynamics of the i agent. The control inputs $u_{p_i}$ and $u_{a_i}$ are second order sliding mode controllers in order to obtain a guidance and navigation for the i aerial vehicle of the MAS as a networked control system.^19,21,29

Remark 1: As the rotation of the four propellers on the quadrotor is balanced, gyroscopic moment $d_{{\eta }_{\mathit{gyr}{o}_i}}$ will essentially be zero. The only case in which gyroscopic moments will not be zero is if there is a significant difference in the RPM of the four motors and in the presence of a strong sideways cross-wind.

Remark 2: The real-time flight scenario consists of quadrotors operating at a determined position and attitude. The ground effect and sideways cross-wind may cause undesired changes during the flight of the aerial vehicles. In this work, the aerodynamic forces as D and Y and moments are considered as time-varying disturbances.

Sectorial fuzzy consensus for formation flight

For this work, a formation flight is defined in which a group of aircraft is accommodated in a predetermined or desired position by means of a mathematically defined geometry, which is a predefined geometric shape for formation flight in three-dimensional space. Thus, the aerial vehicles perform the trajectory tracking within the time-varying geometric shape of the desired formation flight. While the geometry of the formation flight is specified by the desired displacements or movements with respect to an inertial coordinate system, under the assumption that each aircraft is capable of detecting the relative positions and velocities of its neighboring aircraft, with respect to the inertial coordinate system. This implies that the aerial vehicle needs to know its positions and velocities in the inertial coordinate system. ³⁰

Regarding the geometric shape that is defined for the formation flight, it is necessary to assign a position displacement in the three-dimensional plane XYZ for each vehicle considering a desired position reference and respecting the position and velocity among vehicles in order to maintain the desired geometric formation flight. This procedure is mentioned later in the section of simulation results and in the section of experimental results, where several formation flights of quadrotor UAVs were implemented.

In order to describe information exchange among agents of the MAS, the concept of graph is used. A graph is denoted as $G=(\mathcal{V},\mathcal{E})$ , where $\mathcal{V}=\{v_1,\dots ,v_N\}$ represents the set of N nodes and $\mathcal{E}$ represents the set of edges. An edge is denoted as $\mathcal{E}(i,j)$ and is graphically represented by an arrow with tail node i and head node j with information from agent i to agent j. An agent i is called a neighbor of agent j if $(i,j)\in \mathcal{E}$ and the set of neighbors of agent i is denoted as $N_i=\left\{j\right|(i,j)\in \mathcal{E}\}$ . the adjacency o connectivity matrix is defined as $\mathcal{A}=[a_{ij}]\in {ℝ}^{N\times N}$ with weights a_ij > 0 if $(j,i)\in \mathcal{E}$ and a_ij = 0 otherwise. A weighted in-degree of node i is defined as $d_i={\displaystyle {\sum }_{i=1}^N{a}_{ij}}$ and in-degree matrix as $\mathcal{D}=\text{diag}\left\{d_i\right\}\in {ℝ}^{N\times N}$ . Then, the Laplacian matrix is defined as $\mathcal{L}=D-A$ . A node is balanced if its in-degree is equal to its out-degree, i.e., $d_i={\displaystyle {\sum }_{i=1}^N{a}_{ij}}=d_i={\displaystyle {\sum }_{i=1}^N{a}_{ji}}$ . In this work, graphs are considered to be invariant over time; it means, $\mathcal{A}$ is formed by constant terms.^31–34

The distributed consensus formation protocol of the MASs allows an agent to take distributed decisions on the local information resulting in a collective motion of all the agents in the group.

Definition 1: In order to perform consensus, a linear system is defined as follows

$$\dot{\lambda }=A\lambda +B{u}_{\zeta }$$

(10)

$$\zeta =C\lambda$$

(11)

where $\lambda ={[x_i,{\dot{x}}_i,y_i,{\dot{y}}_i,z_i,{\dot{z}}_i,{\varphi }_i,{\dot{\varphi }}_i,{\theta }_i,{\dot{\theta }}_i,{\psi }_i,{\dot{\psi }}_i]}^{\mathrm{\top }}\in {ℝ}^{2DN\times 1}$ is the position and velocity vector of the agent $i=1,\dots ,N,u_{\zeta }\in {ℝ}^{DN\times 1}$ is the corresponding control input vector, $\zeta \in {ℝ}^{DN\times 1}$ is the output vector of the agent i, and $A\in {ℝ}^{2DN\times 2DN},B\in {ℝ}^{2DN\times DN}$ and $C\in {ℝ}^{DN\times 2ND}$ are constant matrices with compatible dimensions. The parameters D and N in the dimension of the system are the degrees of freedom of the agent i to perform consensus and the number of agents, respectively.

$$A=\left(\begin{array}{ccccc}0& 1& 0& \dots & 0\\ 0& 0& 0& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & 0\\ 0& 0& 0& \dots & 1\end{array}\right)B=\left(\begin{array}{ccccc}0& 0& 0& \dots & 0\\ 1& 0& 0& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & 0\\ 0& 0& 0& \dots & 1\end{array}\right)$$

$$C=\left(\begin{array}{ccccc}1& 0& 0& \dots & 0\\ 0& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & 0\\ 0& 0& 0& \dots & 0\end{array}\right)$$

Definition 2: A multi-agent system is said to achieve the consensus if for each agent $i=\{1,\dots ,N\}$ there exists a local state feedback control input $u_{\zeta }$ such that the closed-loop system satisfies

$$\underset{t\to \infty }{\lim }\left|\right|{\zeta }_i(t)-{\zeta }_j(t)\left|\right|=\left|\right|{\zeta }_{d_i}(t)-{\zeta }_{d_j}(t)\left|\right|$$

(12)

$$\underset{t\to \infty }{\lim }\left|\right|{\dot{\zeta }}_i(t)-{\dot{\zeta }}_j(t)\left|\right|=\left|\right|{\dot{\zeta }}_{d_i}(t)-{\dot{\zeta }}_{d_j}(t)\left|\right|=0$$

(13)

and adding the aircraft trajectory tracking control, the following conditions must be reached

$$\underset{t\to \infty }{\lim }\left|\right|{\zeta }_i(t)-{\zeta }_{d_i}(t)\left|\right|=0$$

(14)

$$\underset{t\to \infty }{\lim }\left|\right|{\dot{\zeta }}_i(t)-{\dot{\zeta }}_{d_i}(t)\left|\right|=0$$

(15)

$$\underset{t\to \infty }{\lim }\left|\right|{\zeta }_j(t)-{\zeta }_{d_j}(t)\left|\right|=0$$

(16)

$$\underset{t\to \infty }{\lim }\left|\right|{\dot{\zeta }}_j(t)-{\dot{\zeta }}_{d_j}(t)\left|\right|=0$$

(17)

for any initial condition $\lambda (0)$ , with $i=1,\dots ,N$ , where $\zeta ={[{\zeta }_1,{\zeta }_2,{\zeta }_3,\dots ,{\zeta }_{DN}]}^{\mathrm{\top }}\in {ℝ}^{DN\times 1}$ represents a formation position vector, and ${\zeta }_d={[{\zeta }_{d_1},{\zeta }_{d_2},{\zeta }_{d_3},\dots ,{\zeta }_{d_{DN}}]}^{\mathrm{\top }}\in {ℝ}^{DN\times 1}$ represents a desired formation position vector, and where $\dot{\zeta }={[{\dot{\zeta }}_1,{\dot{\zeta }}_2,{\dot{\zeta }}_3,\dots ,{\dot{\zeta }}_{DN}]}^{\mathrm{\top }}\in {ℝ}^{DN\times 1}$ represents a formation velocity vector, and ${\dot{\zeta }}_d={[{\dot{\zeta }}_{d_1},{\dot{\zeta }}_{d_2},{\dot{\zeta }}_{d_3},\dots ,{\dot{\zeta }}_{d_{DN}}]}^{\mathrm{\top }}\in {ℝ}^{DN\times 1}$ represents a desired formation velocity vector. For the consensus formation flight, and for the tracking control of aircraft trajectories, it must be fulfilled that $\underset{t\to \infty }{\lim }{\dot{\zeta }}_i(t)={\dot{\zeta }}_j(t)$ , as well as ${\dot{\zeta }}_{d_i}(t)={\dot{\zeta }}_{d_j}(t)$ .

Fuzzy logic formation controller

For the fuzzy logic controller formation, $\chi (t)={[{\tilde{e}}_{{\zeta }_i}(t){\dot{\tilde{e}}}_{{\zeta }_i}(t)]}^{\mathrm{\top }}$ is considered and represents the error vector of the position and velocity respectively and $\Phi ({\tilde{e}}_{{\zeta }_i}(t),{\dot{\tilde{e}}}_{{\zeta }_i}(t))$ represents the fuzzy logic consensus protocol that is an additional part of the formation flight controller. The proposed control input, defined as a consensus protocol, has the similar structure to this presented in Oh et al.. ³⁰ Then, the consensus protocol for the formation flight is described as

$$u_{{\zeta }_i}(t)=-K_g\chi (t)-\Phi ({\tilde{e}}_{{\zeta }_i}(t),{\dot{\tilde{e}}}_{{\zeta }_i}(t))$$

(18)

with

$$K_g=\left[\begin{array}{cc}{k}_{g_p}& 0\\ 0& k_{g_d}\end{array}\right]$$

where K_g is positive definite matrix with compatible dimensions. If for any $i,j\in \{1,2,\dots ,N\},\underset{t\to \infty }{\lim }\left|\right|{\zeta }_i(t)-{\zeta }_j(t)\left|\right|=\left|\right|{\zeta }_{d_i}(t)-{\zeta }_{d_j}(t)\left|\right|$ , and $\underset{t\to \infty }{\lim }\left|\right|{\dot{\zeta }}_i(t)-{\dot{\zeta }}_j(t)\left|\right|=\left|\right|{\dot{\zeta }}_{d_i}(t)-{\dot{\zeta }}_{d_j}(t)\left|\right|$ , then MAS (10) is said to achieve the global consensus regime for the state output variables. Moreover, if for any $i,j\in \{1,2,\dots ,N\},\underset{t\to \infty }{\lim }\left|\right|{\zeta }_i(t)-{\zeta }_{d_i}(t)\left|\right|=0,\underset{t\to \infty }{\lim }\left|\right|{\dot{\zeta }}_i(t)-{\dot{\zeta }}_{d_i}(t)\left|\right|=0,\underset{t\to \infty }{\lim }\left|\right|{\zeta }_j(t)-{\zeta }_{d_j}(t)\left|\right|=0$ and $\underset{t\to \infty }{\lim }\left|\right|{\dot{\zeta }}_j(t)-{\dot{\zeta }}_{d_j}(t)\left|\right|=0$ , then MAS (10) is said to achieve trajectory tracking control for the state output variables. Consequently, the consensus error variables as

$${\tilde{e}}_{{\zeta }_i}(t)={\displaystyle \sum _{j\in N_i}{a}_{ij}}(e_{{\zeta }_i}(t)-e_{{\zeta }_j}(t))$$

(19)

$${\dot{\tilde{e}}}_{{\zeta }_i}(t)={\displaystyle \sum _{j\in N_i}{a}_{ij}}({\dot{e}}_{{\zeta }_i}(t)-{\dot{e}}_{{\zeta }_j}(t))$$

(20)

with the position errors of the agents i and j respectively

$$e_{{\zeta }_i}(t)={\zeta }_{d_i}(t)-{\zeta }_i(t)$$

(21)

$$e_{{\zeta }_j}(t)={\zeta }_{d_j}(t)-{\zeta }_j(t)$$

(22)

and, with the velocity errors of the agents i and j respectively

$${\dot{e}}_{{\zeta }_i}(t)={\dot{\zeta }}_{d_i}(t)-{\dot{\zeta }}_i(t)$$

(23)

$${\dot{e}}_{{\zeta }_j}(t)={\dot{\zeta }}_{d_j}(t)-{\dot{\zeta }}_j(t)$$

(24)

where N_i represents the neighbor set of agent i with $i=1,2,\dots ,N$ . The vector, ${\zeta }_i(t)={[x_i,y_i,z_i,{\varphi }_i,{\theta }_i,{\psi }_i]}^{\mathrm{\top }}$ is the position of the ith vehicle, and the vector ${\dot{\zeta }}_i(t)={[{\dot{x}}_i,{\dot{y}}_i,{\dot{z}}_i,{\dot{\varphi }}_i,{\dot{\theta }}_i,{\dot{\psi }}_i]}^{\mathrm{\top }}$ is the velocity of the ith vehicle, $\Phi ({\tilde{e}}_{{\zeta }_i}(t),{\dot{\tilde{e}}}_{{\zeta }_i}(t))$ is the function vector corresponding to fuzzy logic consensus term that guarantees that the position errors $e_{{\zeta }_i},e_{{\zeta }_j}$ and that velocity errors ${\dot{e}}_{{\zeta }_i},{\dot{e}}_{{\zeta }_j}$ , all reach the consensus to maintain formation flight.

To summarize the whole coordination scheme, in Figure 3(a) block diagram is presented, where the consensus protocol is now referred as the Sectorial Fuzzy Consensus Manager.

Figure 3.

Blocks diagram: The consensus protocol is depicted as the Sectorial Fuzzy Consensus Manager.

Each quadrotor UAV is based on a guidance, navigation and control approach in SE(3) (special Euclidean group) and NED system (North-East-Down) as a local algorithm. Then, a sliding mode control (SMC) is used for each quadrotor UAV considering the complete mathematical model, for more details see literature.^19,21 In this sense, the trajectory tracking block is directly connected to the Sectorial Consensus Manager block which runs the consensus algorithm of the N quadrotor UAVs, see Figure 3. Thus, the sectorial consensus manager operates as a distributed protocol that is used to perform formation flight for the MASs or multiple quadrotor UAVs in which each agent or vehicle can take distributed decisions on the local information resulting in a collective motion of all the agents in the group. The proposed algorithm guarantees consensus of the MASs so that the states of all the systems reach the same value.

Theorem 1: Consider the linear systems, (10)-(11), and the consensus protocol (18) in a closed-loop system. Then, the consensus protocol guarantees an asymptotic convergence of the position consensus errors ${\tilde{e}}_{{\zeta }_i}(t),{\tilde{e}}_{{\zeta }_j}(t)$ , and also guarantees an asymptotic convergence of the velocity consensus errors ${\dot{\tilde{e}}}_{{\zeta }_i}(t),{\dot{\tilde{e}}}_{{\zeta }_j}(t)$ , and ensures that the consensus of a multi-agent system is achieved, $\zeta (t)\to {\zeta }_d(t)$ as $t\to \infty$ , and $\dot{\zeta }(t)\to {\dot{\zeta }}_d(t)$ as $t\to \infty$ .

The stability analysis is presented in Appendix 1 of this paper.

Rule base and the inference mechanism

The set of membership functions $A_1^{l_1}({\chi }_1)$ and $A_2^{l_2}({\chi }_2)$ , associated with each of the input variables χ₁ and χ₂, is defined as

$${\mathcal{H}}_1=\left\{A_1^{l_1}({\chi }_1),\text{ for }{l}_1=-\frac{M_1-1}{2},\dots ,0,\dots ,\frac{M_1-1}{2}\right\}$$

(25)

$${\mathcal{H}}_2=\left\{A_2^{l_2}({\chi }_2),\text{ for }{l}_2=-\frac{M_2-1}{2},\dots ,0,\dots ,\frac{M_2-1}{2}\right\}$$

(26)

where M₁ and M₂ are odd constants that represent the number of scalar membership functions associated with each input variable, thus forming a fuzzy partition of the input variables, which in this case satisfy the properties mentioned in literature ³⁵ and are expressed below:

• For all ${\chi }_r\in U_r\subset ℝ$ , the sum of all membership functions associated with each input variables ${\chi }_1(t)={\tilde{e}}_{{\zeta }_i}(t)$ and ${\chi }_2(t)={\dot{\tilde{e}}}_{{\zeta }_i}(t)$ must have a unit value (condition of orthogonality), i.e.

  $${\displaystyle \sum _{l_r=-\frac{M_r-1}{2}}^{\frac{M_r-1}{2}}{A}_r^{l_r}}({\chi }_r)=1\text{para}r=1,2$$

  (27)

• For each $k_r\in \{-\frac{M_r-1}{2},\dots ,0,\dots ,\frac{M_r-1}{2}-1\}$ and $r=1,2,$ the membership functions

  $$A_r^{k_r}({\chi }_r)\text{ and }{A}_r^{k_r+1}({\chi }_r)$$

  (28)

which are adjacent.

• for all ${\chi }_r\in U_r\subset ℝ$ , exists

• $k_r\in \{-\frac{M_r-1}{2},\dots ,0,\dots ,\frac{M_r-1}{2}-1\}$ such that

  $${\displaystyle \sum _{l_r=-\frac{M_r-1}{2}}^{\frac{M_r-1}{2}}{A}_r^{l_r}}({\chi }_r)-{\displaystyle \sum _{l_r=k_r}^{k_r+1}{A}_r^{l_r}}({\chi }_r)=0$$
  $$\text{ for }r=1,2.$$

  (29)

This means that, for any possible value of ${\chi }_r\in U_r\subset ℝ$ , at most, only two membership functions produce values (membership grade) strictly different from zero.

The input membership functions $A_1^{l_1}({\chi }_1)$ and $A_2^{l_2}({\chi }_2)$ are symmetrical with respect to the origin, such that

$$A_r^{l_r}({\chi }_r)=A_r^{-l_r}(-{\chi }_r)$$

(30)

for all ${\chi }_r\in U_r$ , with r = 1, 2.

Figures 4 and 5 illustrate the input membership functions used for the consensus formation protocol.

Figure 4.

Input membership function of consensus position.

Figure 5.

Input membership function of consensus velocity.

The set of parameters ${\overline{\chi }}_r^{l_r}$ and $p_{A_r},$ where $p_{A_r}=\{p_{-2r},p_{-1r},p_{0r,}{p}_{1r},p_{2r}\}$ and ${\overline{\chi }}_i=\{{\overline{\chi }}_r^{-2},{\overline{\chi }}_r^{-1},{\overline{\chi }}_r^0,{\overline{\chi }}_r^1,{\overline{\chi }}_r^2\}$ with $p_{A_r}={\overline{\chi }}_r$ that describe each input membership function, associated with χ₁ and χ₂.

The membership function associated with the output variable v(t) is defined by $B^l(v)$ for $v(t)\in V\subset ℝ$ , with $V\subset ℝ$ as the universe supported for the output variable, and with $l=-\frac{J-1}{2},\dots ,\frac{J-1}{2}$ , where J is the number of scalar functions of the output membership function. In this case, the function $B^l(v)$ of sigmoidal form is defined as shown in Figure 6, with $p_B=\left\{{\overline{v}}^l\right\}$ as the set of its centers.

Figure 6.

Output membership function of consensus FLC.

The parameters of $p_B=\{{\overline{v}}^{-2},{\overline{v}}^{-1},{\overline{v}}^0,{\overline{v}}^1,{\overline{v}}^2\}$ , that describe the output membership function, associated with the output variable of FLC controller v(t), see Figure 6.

Lemma 1: The rule base for the fuzzy consensus protocol term with the input variables χ₁ and χ₂, and the output variable v, contemplates $M=M_1{M}_2$ possible combinations of the membership functions inputs such that

$$v=\Phi (\chi )$$

(31)

$$\text{Rule}(l_1{l}_2):$$

$$\text{ if }{\chi }_1\in A_1^{l_1}\text{ and }{\chi }_2\in A_2^{l_2}\text{THEN}v\in B^{l_1{l}_2}$$

(32)

where ${[{\chi }_1{\chi }_2]}^T=\chi \in U_1\times U_2\subset {ℝ}^2,v\in V\subset ℝ$ , and $A_r^{l_r}(\chi )\in {\mathcal{H}}_r$ , defined in equations (25) and (26), with r = 1, 2.

Each membership function $B^{l_1{l}_2}(v)$ is associated with an i-th membership function $B^l(v)$ ; this is

$$B^{l_1{l}_2}(v)\in \{B^{-\frac{J-1}{2}}(v),\dots ,B^{\frac{J-1}{2}}(v)\}$$

(33)

where $l_1=-\frac{M_1-1}{2},\dots ,\frac{M_1-1}{2}$ and $l_2=-\frac{M_2-1}{2},\dots ,\frac{M_2-1}{2}$ . In our case, in Table 1, the input membership functions describe $A_k^{-2}$  = NB (Negative Big), $A_k^{-1}$  = NS (Negative Small), $A_k^0$  = ZO (approximately zero), $A_k^1$  = PS (Positive Small), $A_k^1$  = PB (Positive Big), with k = 1, 2. Therefore for k = 1 it results $M_1=5$ input membership functions for the position consensus error, and for k = 2 it results $M_2=5$ input membership functions for the velocity consensus error. Moreover, in Table 1 the output membership functions describe $B^{-2}$  = NB (Negative Big), $B^{-1}$  = NS (Negative Small), B⁰ = ZO (approximately zero), B¹ = PS (Positive Small), B² = PB (Positive Big), so it results J = 5 output membership functions for the fuzzy consensus controller. Then, it results that $M_1=M_2=J=5$ , such that the number of total rules for the fuzzy inference algorithm is obtained as $M=M_1{M}_2=25$ . The IF-THEN rule base is defined from the base of fuzzy rules shown in Table 1.

Table 1.

Tuning of fuzzy rules.

${\chi }_1\backslash {\chi }_2$	$A_1^{-2},NB$	$A_1^{-1},NS$	$A_1^0,ZO$	$A_1^1,PS$	$A_1^2,PB$
$A_2^{-2},NB$	$B^{-2},NB$	$B^{-2},NB$	$B^{-2},NB$	$B^{-1},NS$	$B^0,ZO$
$A_2^{-1},NS$	$B^{-2},NB$	$B^{-2},NB$	$B^{-1},NS$	$B^0,ZO$	$B^1,PS$
$A_2^0,ZO$	$B^{-2},NB$	$B^{-1},NS$	$B^0,ZO$	$B^1,PS$	$B^2,PB$
$A_2^1,PS$	$B^{-1},NS$	$B^0,ZO$	$B^1,PS$	$B^2,PB$	$B^2,PB$
$A_2^2,PB$	$B^0,ZO$	$B^1,PS$	$B^2,PB$	$B^2,PB$	$B^2,PB$

The inference mechanism uses the base of fuzzy rules and defines the operation that will be used to calculate the function implied or consequential. In our case, the minimal inference process is used.^36,37 If the membership function of the fuzzifier is singleton, the membership function involved is expressed as:

$$B_{l_1{l}_2}^{\prime }({[\chi v]}^T)=\text{min}({\omega }_{l_1{l}_2},B^{l_1{l}_2}(v))$$

(34)

where ${\omega }_{l_1{l}_2}=\text{min}(A_1^{l_1}({\chi }_1),A_2^{l_2}({\chi }_2))$ .

This inference mechanism is also called a clipping method. ³⁷

The inference mechanism produces $M=M_1{M}_2$ membership functions defined by $B_{l_1{l}_2}^{\prime }({[\chi v]}^T)$ , that must be combined in order to obtain the set involved defined by $B^{\prime }({[\chi \upsilon ]}^T)$ , as expressed in equation (35), considering the definition of the norm- s maximum represented by the union operator $\dot{+}$ ³⁸

$$B^{\prime }({[\chi v]}^T)=B_{\frac{-(M_1-1)}{2}\frac{-(M_2-1)}{2}}^{\prime }({[\chi v]}^T)\dot{+}\dots \dot{+}{B}_{\frac{M_1-1}{2}\frac{M_2-1}{2}}^{\prime }({[\chi v]}^T)$$

(35)

Defuzzification module

The process of defuzzification is the step where the membership function $B^{\prime }({[xv]}^T)$ calculated in the inference mechanism is used to obtain the value of the variable $v\subset ℝ$

To calculate the output variable of the consensus fuzzy protocol, the average center defuzzification method is used and calculated by means of the weighted arithmetic mean of the set implicated, as expressed in equations (36) and (37).

$$\Phi (\chi )=\frac{{\displaystyle {\sum }_{l_1=-\frac{M_1-1}{2}}^{\frac{M_1-1}{2}}}{\displaystyle {\sum }_{l_2=-\frac{M_2-1}{2}}^{\frac{M_2-1}{2}}{\overline{v}}^{l_1{l}_2}}\left[{\displaystyle \underset{r=1}{\overset{2}{\cap }}{A}_r^{l_r}}({\chi }_r)\right]}{{\displaystyle {\sum }_{l_1=-\frac{M_1-1}{2}}^{\frac{M_1-1}{2}}}{\displaystyle {\sum }_{l_2=-\frac{M_2-1}{2}}^{\frac{M_2-1}{2}}\left[{\displaystyle \underset{r=1}{\overset{2}{\cap }}{A}_r^{l_r}}({\chi }_r)\right]}}$$

(36)

$$v=\Phi (\chi )$$

(37)

where $\cap$ denotes the intersection operator which is chosen in this work to represent the minimum operator, and ${\overline{v}}^{l_1{l}_2}$ represents the centers of the membership functions $B_{l_1{l}_2}^{\prime }({[\chi v]}^T)$ .

Numerical simulation

As a first validation tool, numerical simulations are developed. For this case, a geometrical formation, consisting of six agents arranged in a delta-like array, is simulated.

A directed spanning tree graph is selected as the communication arrangement between agents, this is shown in Figure 7. Note that the first agent is not receiving information from any of its neighbors, meaning it is acting as the leader of the formation; thus, what affects it, affects the other agents through its neighbors. In this sense, the connectivity matrix corresponding to this connection graph is defined as

$$\mathcal{A}=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\end{array}\right]$$

(38)

which leads to its corresponding input matrix as

$$\mathcal{D}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]$$

(39)

Figure 7.

Directed connection graph.

obtaining the following Laplacian matrix

$$\mathcal{L}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ -1& 1& 0& 0& 0& 0\\ -1& 0& 1& 0& 0& 0\\ 0& 0& -1& 1& 0& 0\\ -1& 0& 0& 0& 1& 0\\ 0& -1& 0& 0& 0& 1\end{array}\right]$$

(40)

It is necessary to remark that in order to achieve the geometrical formation, independently of the leader trajectory, a physical separation is defined from the actual position states of the agent. In fact, to maintain the agents always pointing at the heading of the leader, they are commanded to always stay behind it; in other words, the geometrical offsets are driven by the body frame of the leader.

For all the following simulations, a geometrical formation with separation of 1000 mm among the agents is selected. In addition, the dynamics of the X and Y-axis movements of an aircraft is simulated by means of the following transfer function

$$G(s)=\frac{170.5}{s^2+0.7395s+1.403}\frac{1}{s}$$

(41)

which is experimentally obtained from a Parrot AR. Drone 2.0, by a step-test procedure.

To test the overall algorithm scheme, four scenarios are selected to be proven; i.e., a simple L-shaped path, to test the sudden change of direction of the formation, the same L-shaped path with the addition of a simulated wind disturbance, and a Lemniscate trajectory with and without the simulated wind disturbance.

L-shaped path

For the case of this first scenario, an L-shaped path is commanded to the leader of the formation. The previous means that all the followers will remain behind it, even in the turning part of the path, which in this case leads the neighbor agents to increment their velocities; however, due the fact that the Consensus Manager is coordinating the states, the formation is capable of performing the turn as seen in Figure 8. Note that the agents did not start within their positions in the formation, meaning they also needed to be coordinated in the arrangement process.

Figure 8.

2D path described by the six agents of the formation.

The individual axis trajectories are also splitted in the X and Y-axis and are graphed in Figure 9.

Figure 9.

Positions in (a) X and (b) Y-axis of the six agents.

As previously stated, the Consensus Manager is the layer which deals with the determination of the actual desired states that the aircraft must follow; this means that the ideal input trajectory is then modified and adjusted by the Manager, according to the overall behavior of the formation. In this case, the Fuzzy-part of the Consensus Manager output is extracted and showed in Figure 10 for the five followers and separated into individual axis components. It is worth to mention that at the beginning of the trajectory, there is a sudden rise of the signals, that can be explained by the fact that the agents started with initial conditions different from those on the corresponding formation, so the coordination algorithm responded with fast coupling signals.

Figure 10. (a)

X and (b) Y fuzzy components for the six agents.

Consensus Manager complete signals, as input to the follower agents, can be seen in Figure 11.

Figure 11.

(a) X and (b) Y-axis Consensus Manager output used as reference by the six agents.

L-shaped path with disturbances

One of the main purposes to design a coordination algorithm, is to give the formation the faculty of being disturbance-tolerant. For instance, in aircraft the most common disturbance is produced by wind gusts. For this reason, this second test also simulates a wind disturbance located at coordinates X = –4000 mm and Y = –100 mm, with an effective angle of 40 degrees, and an action range of straight 4000 mm. The wind gust is simulated using the following equation

$$Ac{c}_{\mathit{wind}}=4000\sin (2t)+500\cos (15t)+900\sin (10t)+4000$$

(42)

In addition, this simulated acceleration is splitted into local X and Y components, depending on the relative angle between the agent and the wind gust origin; meaning it affected both movements of the aircrafts. This behavior is graphed in Figure 12, as axis-component signals.

Figure 12.

(a) X and (b) Y-axis simulated wind disturbances for the six agents.

The 2D trajectories performed by the formation can be seen in Figure 13, while the individual axis components, in Figure 14. Note that despite being in a disturbance zone, the formation is maintained, meaning the coordination algorithm is capable of adjusting its response, even in the presence of intermittent disturbances.

Figure 13.

2D path described by the six agents of the formation.

Figure 14.

Positions in (a) X and (b) Y-axis of the six agents.

Figure 15 graphs the fuzzy-part of the Consensus Manager output. For this disturbed trajectory, the response of the coordinator against the simulated wind gust can be seen about the 13 s, which is the time where the agents started to cross the disturbance zone.

Figure 15.

(a) X and (b) Y fuzzy components for the six agents.

The complete components of the Consensus Manager outputs are displayed in Figure 16.

Figure 16.

(a) X and (b) Y-axis Consensus Manager output used as reference by the six agents.

Figure 17.

2D path described by the six agents of the formation.

Figure 18.

Positions in (a) X and (b) Y-axis of the six agents.

Lemniscate path

As an accuracy test, a recursive lemniscate-like path is also tested. As the previous simulations, the leader follows the original lemniscate, while the five followers always remain behind it, according with the leader heading. This trajectory is designed in the form of parametric equations, and only the leader take them as input, while the followers, as previously stated, will only receive the states of the leader and the relative orientation of the complete formation. These parametric equations are described as follows

$${\zeta }_{d_i}(t)=\left[\begin{array}{c}{\zeta }_{d_{i,X}}(t)\\ {\zeta }_{d_{i,Y}}(t)\end{array}\right]=\left[\begin{array}{c}r\text{sin}\left(\omega t\right)\\ r\text{cos}\left(2\omega t\right)\end{array}\right]$$

(43)

$${\dot{\zeta }}_{d_i}(t)=\left[\begin{array}{c}{\dot{\zeta }}_{d_{i,X}}(t)\\ {\dot{\zeta }}_{d_{i,Y}}(t)\end{array}\right]=\left[\begin{array}{c}\omega r\text{cos}\left(\omega t\right)\\ -2\omega r\text{sin}\left(2\omega t\right)\end{array}\right]$$

(44)

where r represents the magnitude of the curve, ω the angular frequency, and t the time. For the case of this simulation, the parameters are selected as $r=5000\text{mm}$ and $\omega =0.2\text{rad}/\text{s}$ .

The 2D trajectory performed by the formation is displayed in Figure 17, while the individual axis components in Figure 18.

Fuzzy components of the Consensus Manager are graphed in Figure 19. As previous simulations, the fact that the initial positions of the agents are different from those on the formation, allow seeing the arrangement process, and because of it, a sudden rise in the components at the beginning of the trajectory.

Figure 19.

(a) X and (b) Y fuzzy components for the six agents.

Finally, the complete Consensus Manager outputs can be seen in Figure 20 for the five follower agents.

Figure 20.

(a) X and (b) Y-axis Consensus Manager output used as reference by the six agents.

Lemniscate path with disturbances

As a final simulation, the same lemniscate path is used, but this time the same simulated wind disturbance, is also considered. This wind gust is previously described by equation (42), and it is now positioned in X =–2500 mm and Y = 0 mm, with an effective angle of 40 degrees, and a maximum coverage distance of 4000 mm, causing the simulated accelerations as depicted in Figure 21.

Figure 21.

(a) X and (b) Y-axis simulated wind disturbances for the six agents.

The achieved trajectories are shown in Figure 22 for the six agents. This same behavior can be seen in the individual components in Figure 23. Note that when the formation passes through the center intersection of the path; this is, the disturbance zone, the algorithm is still capable of maintaining the geometrical formation.

Figure 22.

2D path described by the six agents of the formation.

Figure 23.

Positions in (a) X and (b) Y-axis of the six agents.

Consensus Manager fuzzy components are graphed in Figure 24. In these components, multiple signal rises can be spotted; it is due to the wind disturbance affecting the formation each time the last passes by the disturbance zone.

Figure 24.

(a) X and (b) Y fuzzy components for the six agents.

The complete Consensus Manager output signals are displayed in Figure 25.

Figure 25.

(a) X and (b) Y-axis Consensus Manager output used as reference by the six agents.

Experimental results

The Navigation Laboratory in the Center Aerospace Engineering Research and Innovation Center of the Faculty of Mechanical and Electrical Engineering at the Autonomous University of Nuevo Leon has a flight area of 6 m by 10 m, with a maximum effective flight height of 8 m. The motion capture system is able to track retroreflective markers placed in each UAV. The system allows us to obtain high accuracy measurements in the translational and rotational movements of an aerial vehicle. The 16 Vicon T40 cameras possess a sensor resolution of 4 megapixels and are capable of capturing frames at velocities of up to 340 fps. The experimental platform consists of two Parrot AR.Drones 2.0 and the motion capture system, allowing us to track the position and orientation of an object with high accuracy up to 0.1 mm of translation and 0.1 degrees of rotation. The technical characteristics of the employed experimental platform and aerial vehicles can be found in detail in literature.^39,40

The technical configuration of processing units is divided into a ground station computer and individual onboard aircraft computers. The ground station gets the translational data from the Vicon motion capture system, and runs the coordination algorithm, as well as each individual position controllers of the aircraft. On the other hand, each independent onboard computer runs the corresponding orientation controller, and determines the orientation states from the embedded inertial unit. This implementation can be considered distributed due to the fact that every algorithm runs independently in a separated thread of the corresponding computer. In addition, there are individual and independent position controllers for each aircraft that gets data from the main coordination algorithm. These individual position controllers then communicate independently to their corresponding onboard aircraft computer. In fact, this allows the aircraft to take distributed decisions based on the local information obtained from the coordination algorithm, which in this sense, is a result of the collective motion. A configuration scheme can be seen in Figure 26.

An elliptical path was designed in order to test the capability of the algorithm to maintain a straight-line formation when using a bi-directional connected graph of communication by two agents, see Figure 27; and for testing all the movement axes, a variable height lemniscate-like path was followed by the agents using the same proposed graph.

Figure 26.

Experimental platform.

For both experiments, the connectivity matrix $\mathcal{A}$ is defined as

$$\mathcal{A}=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]$$

(45)

while the input degree matrix is calculated as

$$\mathcal{D}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$

(46)

obtaining the following expression for the Laplacian matrix

$$\mathcal{L}=\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$$

(47)

The Laplacian matrix is obtained from the connections between the linear dynamic models of the MAS aircraft through this matrix; it is possible to establish the connectivity of the global closed-loop linear system to maintain a feasible formation flight. In equation (47) it can be observed that when using an undirected graph for the connectivity of a MAS, the Laplacian matrix is symmetric. On the other hand, note that the Laplacian matrix is diagonally dominant and has no negative entries in its diagonal, and also the Laplacian matrix is positive semidefinite, i.e., $\mathcal{L}\ge 0$ . These properties of the Laplacian matrix confirm that the connectivity of the MAS allows formation flight to be achievable by the group of quadrotor UAVs.

Experimental tuning parameters

In order to tune the parameters of the K_g matrix of the consensus protocol presented in equation (18), it is necessary to consider the physical limitations and geometry of the formation.

As stated previously, the protocol is formed by two terms, the fuzzy logic part, and the static formation term containing K_g matrix. In this sense, this last term guarantees that UAVs reach consensus maintaining the same behavior. The previous statements mean that the physical distance and maximum bounded disturbance determine the corresponding values for the $k_{g_p}$ term of the matrix. Thus, this value is adjusted depending on the distance between the selected vehicle and its neighbors, the number of interconnections, and the selected maximum distance that the vehicle is allowed to move even in the presence of disturbances; for the case of these experiments, a value of 1500 mm is the greater allowed and 1 connection. After a series of test and calibration using these considerations, a $k_{g_p}=1$ value was found to be the adequate to achieve the desired behavior.

The subterm $k_{g_d}$ deals with the maximum velocity at which the vehicle can react against the disturbance. For the case of these experiments, a velocity of 2000 mm/s is considered as the maximum value the vehicle can reach.

The tuned values of the parameters of the two inputs and one output membership functions were assigned from the physical limits of the consensus variables of position and velocity respectively, and their numerical values are indicated in equations (48), (49) and (51) respectively with five membership functions for each input variable ( $M_1=5,M_2=5$ ) and five membership functions for the output variable J = 5. The generated surface, using these parameters, is shown in Figure 28.

$$p_{A_1}=\{(-1000),(-500),(0),(500),(1000)\}[\text{mm}]$$

(48)

$$p_{A_2}=\{(-2000),(-1000),(0),(1000),(2000)\}\left[\frac{\text{mm}}{\text{s}}\right]$$

(49)

$$p_B=\{(-1000),(-500),(0),(500),(1000)\}[\text{mm}].$$

(50)

Figure 27.

Connection graph for the two agents.

Figure 28.

Surface generated for the fuzzy term in the Sectorial Fuzzy Consensus Manager.

Elliptical path

In the first case, the desired trajectory in the airspace is assigned in the plane formed by the axes X–Y, considering that the height is constant on the Z axis, while for guidance and navigation, the quadrotor maintains its nose in the forward front following the planned trajectory. The desired flight path assigned to each vehicle belongs to the formation flight, which consists of a closed path described by an elliptic function expressed by the equations (51) and (52), adding the corresponding displacements of the positions of each quadrotor UAVs to perform a pursuit formation flight.

$${\zeta }_{d_i}(t)=\left[\begin{array}{c}{\zeta }_{d_{i,X}}(t)\\ {\zeta }_{d_{i,Y}}(t)\end{array}\right]=\left[\begin{array}{c}{r}_1\text{sin}\left(\omega t-\frac{2\pi (i-1)}{l}\right)\\ r_2\text{cos}\left(\omega t-\frac{2\pi (i-1)}{l}\right)\end{array}\right]$$

(51)

$${\dot{\zeta }}_{d_i}(t)=\left[\begin{array}{c}{\dot{\zeta }}_{d_{i,X}}(t)\\ {\dot{\zeta }}_{d_{i,Y}}(t)\end{array}\right]=\left[\begin{array}{c}\omega r_1\text{cos}\left(\omega t-\frac{2\pi (i-1)}{l}\right)\\ -\omega r_2\text{sin}\left(\omega t-\frac{2\pi (i-1)}{l}\right)\end{array}\right]$$

(52)

where r₁ is the radius of the ellipse on the X axis, while r₂ is the radius of the ellipse on the Y axis, such that for the shape of the elliptic path $r_2>r_1$ must be fulfilled. The constant l represents the displacement of the assigned position on the X and Y axes in the path described by the elliptical function for each aircraft of the pursuit formation flight, and ω is the angular frequency of the formation tracking.

For the first experiment, an elliptical path, of $r_1=2000\text{mm},r_2=3000\text{mm}$ and $\omega =0.2\text{rad}/\text{s}$ , is tracked by the agents maintaining a pursuit formation.

In Figure 29 is shown the 3D response of the navigation of both agents.

Figure 29.

3D path described by the two agents.

For a clearer view of the movement, Figure 30 shows the separate X, Y and Z-axes movements for the two agents.

Figure 30.

Positions in (a) X, (b) Y and (c) Z-axis of both agents.

As it was mentioned before, the Consensus Manager has a Fuzzy part which allows the agents to reach consensus. In this sense, the fuzzy part that composes the manager can be analyzed separately in order to study the contribution to the manager. The fuzzy components of the manager are seen in Figure 31.

Figure 31.

X and Y fuzzy components of both agents. (a) Fuzzy components for agent 1. (b) Fuzzy components for agent 2.

The Consensus Manager output, which is used as a desired reference by the position controller, is shown in Figure 32. Note that these signals already contain the fuzzy part described in Figure 31.

Figure 32.

X and Y-axis Consensus Manager output used as references by both agents. (a) Consensus Manager output for agent 1. (b) Consensus Manager output for agent 2.

Figure 33 shows the control signals for the position movement X,Y and Z, and for the heading angle ψ.

Figure 33.

X, Y, Z-axis and -angle control signals of both agents. (a) Control signals of agent 1. (b) Control signals of agent 2.

Figure 34.

3D path described by the two Agents.

Height variable lemniscate path

In this case, the desired trajectory in the airspace is assigned in the plane formed by the axes X–Y, considering that the height varies in the Z axis, while the guidance angle ψ of the vehicle maintain its nose ahead of the preset flight path. The desired flight path assigned to each quadrotor UAV belongs to the formation flight, which consists of a closed path described by a lemniscate function expressed by equations (53) and (54), adding the corresponding displacements of the positions of each quadrotor aircraft to perform a pursuit formation flight.

$${\zeta }_{d_i}(t)=\left[\begin{array}{c}{\zeta }_{d_{i,X}}(t)\\ {\zeta }_{d_{i,Y}}(t)\end{array}\right]=\left[\begin{array}{c}r\text{sin}\left(\omega t-\frac{2\pi (i-1)}{l}\right)\\ r\text{cos}\left(2\omega t-\frac{2\pi (i-1)}{l}\right)\end{array}\right]$$

(53)

$${\dot{\zeta }}_{d_i}(t)=\left[\begin{array}{c}{\dot{\zeta }}_{d_{i,X}}(t)\\ {\dot{\zeta }}_{d_{i,Y}}(t)\end{array}\right]=\left[\begin{array}{c}\omega r\text{cos}\left(\omega t-\frac{2\pi (i-1)}{l}\right)\\ -2\omega r\text{sin}\left(2\omega t-\frac{2\pi (i-1)}{l}\right)\end{array}\right]$$

(54)

where r is the radius of the curve, and where l is the displacement of the assigned position on the X and Y axes in the path described by the lemniscate function for each UAV of the pursuit formation flight, and ω is the angular frequency of the formation tracking.

A height variable path allows the controller to be tested for the two vehicles. For this reason, a lemniscate-like path, with $r=2250\text{mm},\omega =0.2\text{rad}/\text{s}$ and variable height, is performed by the same formation of two vehicles. In Figure 34, the 3D movement is shown. Note that two trajectories are depicted in the figure; this was done for the purpose of showing in a clearer way the movement of both agents meaning that agent 2 has a retardation in the change of height with respect to the agent 1.

Figure 35.

Positions in (a) X, (b) Y and (c) Z-axis of both agents.

In Figure 35, the separated axis movements are shown in order to see in a clearer way the response of the algorithm.

The fuzzy part of the Consensus Manager is seen in Figure 36, meaning that these are the components that allow the formation to reach consensus.

Figure 36.

X and Y fuzzy components of both agents. (a) Fuzzy components for the agent 1. (b) Fuzzy components for the agent 2.

The Consensus Manager output used as desired reference is shown in Figure 37, while Figure 38 shows the control signals corresponding to the movement in the X, Y, and Z axes, as well as the heading angle ψ.

Figure 37.

X and Y-axis Consensus Manager output used as references by both agents. (a) Consensus Manager output for the agent 1. (b) Consensus Manager output for the agent 2.

Figure 38.

X, Y, Z-axis and -angle control signals of both agents. (a) Control signals of the agent 1. (b) Control signals of the agent 2.

The performance of the proposed fuzzy consensus, compared with the PID consensus presented in literature ¹⁹ leads to the following results: the PID consensus tracks the desired trajectory generating some under damped responses in the displacements of the agents and the responses show a slow convergence to the references even in the presence of disturbances due to the vortices produced by the rotors. On the other hand, the sectorial fuzzy operates as position and velocity sectors emulating smooth saturation functions and allows the MAS to reach the global asymptotic convergence about the origin so that the consensus is reached.

Both experiments were captured in video and are presented in the following web link: https://youtu.be/3Tx_F2vqjlM

Conclusions

In this paper, a sectorial fuzzy consensus algorithm has been proposed in order to perform the formation flight. The mathematical model for N-quadrotor UAVs was described using the Newton-Euler formulation. A structure of formation flight was developed considering a sectorial fuzzy controller along with its properties, and a linear system in which position and velocity states of the ith quadrotor are considered in the consensus protocol. The stability analysis was obtained as a generalized form for N-quadrotor UAVs and it has shown that the closed-loop system is globally asymptotically stable reaching the consensus of the UAVs. Simulation results have demonstrated the feasibility of the proposed scheme of the formation flight for multiple quadrotors UAVs even in the presence of disturbances. Finally, in order to validate the proposed sectorial fuzzy algorithms, a series of real-time experiments were executed obtaining the behavior responses of the Multiple UAVs.