Multi-UAV Formation Control in Complex Conditions Based on Improved Consistency Algorithm

Abstract

Formation control is a prerequisite for the formation to complete specified tasks safely and efficiently. Considering non-symmetrical communication interference and network congestion, this article aims to design a control protocol by studying the formation model with communication delay and switching topology. Based on the requirements during the flight and the features of the motion model, the three-degrees-of-freedom kinematics equation of the UAV is given by using the autopilot model of longitudinal and lateral decoupling. Acceleration, velocity, and angular velocity constraints in all directions are defined according to the requirements of flight performance and maneuverability. The control protocol is adjusted according to the constraints. The results show that the improved control protocol can quickly converge the UAV formation state to the specified value and maintain the specified formation with communication delay and switching topology.

1. Introduction

Due to their low cost, strong maneuverability, and wide application range, UAVs have good application prospects whether in the civilian or military fields [1,2]. With the complexity and diversification of mission requirements, the low efficiency of a single UAV has emerged. In order to solve this problem, in addition to improving the function and utility of a single UAV, UAV formation flight has also become a research focus [3]. Formation flight means that drones can fly in an expected formation. When the environment or tasks change, the formation can be changed according to the requirements [4]. The technology of UAV formation has broad development and application prospects, and using drones to fly in the expected formation can allow for the completion of more complex tasks and significant improvement in the efficiency of tasks [5].

The studies on multi-UAV formation focus on formation control, formation reconfiguration, real-time path planning, etc. Formation control is the basis and focus of formation flight. Commonly used formation control methods mainly include the leader−follower method, virtual structure method, and behavior control method. The virtual structure method can simplify the assignment of tasks with high accuracy. The disadvantages are that it is difficult to perform fault-tolerant processing and requires a lot of communication. The most mature traditional formation control method is the leader−follow method [6,7,8,9]. The leader−follower method simplifies the control of the multi-UAV model [10]. However, it still has certain limitations, for example, its tracking error will be propagated backward step by step and thus be amplified. Other methods are combined with the leader−follower method to solve the problems above [11,12,13]. Every aircraft receives the same information, namely the trajectory of the virtual leader in the virtual−leader method. The advantages of the virtual structure method are that it simplifies the description and assignment of tasks, and has high formation control accuracy. The disadvantages are that it is difficult to perform fault-tolerant processing and requires a large amount of communication as a centralized control method. The behavior-based method, which is based on the information obtained from the sensor to determine the responses of the UAVs, has strong robustness and flexibility, but cannot achieve accurate formation maintenance. The work of [14] studies multi-UAV formation by applying the behavior control method.

REN indicates that the above three formation control methods can be unified under the framework of the consistency theory; formation control based on the consistency theory can overcome some shortcomings of these traditional methods [15]. The formation control method based on the consistency theory is such that every agent can realize large-scale and distributed formation control through the communication between neighboring UAVs under a certain communication network without centralized coordination [16]. The impact of interaction models on the coherence of collective decision-making is discussed in [17,18,19].

Formation control methods based on the consistency theory have yielded some valuable research results in recent years [20,21,22,23]. The work of [24] studied the problem of time-varying formation control under the constraint of communication delay and designed a consistent control method that can deal with communication delay. The work of [25] studied the consensus formation control method based on time-varying communication topology. The work of [26] considered the existence of random communication noise and information packet loss constraints in the network and adopted the polygon method of information exchange based on the consistency theory to realize formation control. The work of [27] studied the cooperative formation control problem of the multi-aircraft system based on the consistency theory with a fixed connectivity of the network topology.

Control laws based on consistency are often adopted to solve the problems of multi-UAV formation, and the maneuvering performance and flight performance of UAVs will impose restrictions on the control variables and flight states in different ways. In addition, communication between aircraft is often affected by factors such as transmission speed and network congestion, resulting in communication delay; due to communication interference and complex terrain, the multi-aircraft system network topology changes. Therefore, the research on multi-UAV formation considering communication constraints and flight constraints has important value.

Aiming at the problems above, the main contributions of this paper are as follows: (1) This paper adopts the three-degrees-of-freedom kinematics model of a drone which is based on autopilot, and the lateral heading autopilot and the longitudinal autopilot are decoupled. (2) This paper proposes an improved basic consistency algorithm. During the flying process of drones, the communication constraints, such as topology switching and non-symmetrical communication delay, are considered to design the consistency algorithm. (3) In addition to communication constraints, mobility constraints are also considered to improve the consistency algorithm. Compared with other existing methods based on the consistency algorithm, the improved method considers the formation control in complex conditions. The communication constraints and flight constraints are both considered. The communication constraints include the communication delay and switching topology, and the flight performance and maneuverability constraints include the speed, acceleration, and heading angular velocity of the UAV. The improved algorithm can not only achieve multi-UAV formation control when topology switching and communication delay exist, but it also satisfies the constraints of UAV maneuverability and flight performance.

This article is organized as follows: The three-degrees-of-freedom kinematics model of a drone which is based on autopilot is adopted, and the lateral heading autopilot and the longitudinal autopilot are decoupled in Section 2. Section 3 proposes an improved consistency algorithm that is effective with topology switching and communication delay. The minimum adjustment is used to adjust for flight constraints. Then, the convergence proof of the improved consistency algorithm is given. Section 4 discusses the simulation. The results show that the consistency control protocol proposed can meet the mobility requirements with communication delay and switching topology.

2. UAV Dynamics Modeling and Consistency Algorithm

This section first establishes the coordinate system, describes the formation, and then gives the kinematics model. The consensus algorithm is presented to prepare for the subsequent proposed multi-UAV control protocol with switching topology and communication delay. Finally, the control protocol is adjusted considering the constraints of flight status and maneuverability.

2.1. UAV Formation Description

Firstly, a coordinate system is created to express the position of the UAVs. The UAV is considered a mass point. To describe the movement state of the drone, we use the ground coordinate system. On the horizontal plane, the origin $O$ can be arbitrarily selected.

There are two ways to describe the three-dimensional plane of the UAV formation, the $l-\psi$ method and the $l-l$ method in [28]. In this paper, the method $l-l$ is selected. The positional relationship between UAVs can be described by the relative positional relationship matrix $R_x,R_y,R_z$.

$$\left\{\begin{array}{r}{R}_x=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1n}\\ x_{21}& x_{22}& \cdots & x_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ x_{n1}& x_{n2}& \cdots & x_{nn}\end{array}\right]\\ R_y=\left[\begin{array}{cccc}{y}_{11}& y_{12}& \cdots & y_{1n}\\ y_{21}& y_{22}& \cdots & y_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ y_{n1}& y_{n2}& \cdots & y_{nn}\end{array}\right]\\ R_z=\left[\begin{array}{cccc}{z}_{11}& z_{12}& \cdots & z_{1n}\\ z_{21}& z_{22}& \cdots & z_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ z_{n1}& z_{n2}& \cdots & z_{nn}\end{array}\right]\end{array}\right.$$

(1)

where $\left(x_{ij},y_{ij},z_{ij}\right)\left(i,j=1,\cdots ,n\right)$ describes the difference in position between two drones. $x_{ii}=y_{ii}=z_{ii}=0$.

The conditions when multi-UAV forms a stable, desired formation are as follows:

$$\left\{\begin{array}{c}\left|x_i-x_j\right|\to x_{ij}\\ \begin{array}{l}\left|y_i-y_j\right|\to y_{ij}\\ \left|z_i-z_j\right|\to z_{ij}\end{array}\\ \left|v_i-v_j\right|\to 0\end{array}\right.$$

(2)

where $x_i,y_i,z_i$ are coordinates for UAVs. $v_i$ is velocity.

2.2. UAV Kinematics Model

In the UAV formation, the three-degrees-of-freedom kinematics model with autopilot is usually adopted. The longitudinal and lateral movements of the basic kinematic equations of UAV formation control are coupled. The work of [29] decouples the lateral heading autopilot and the longitudinal autopilot and obtains a kinematic model of lateral and longitudinal separation. The motion model of the UAV#$i$ is given by Equation (3):

$$\left\{\begin{array}{l}{\dot{x}}_i=v_i\cos {\theta }_i\\ {\dot{y}}_i=v_i\sin {\theta }_i\\ {\dot{\theta }}_i={\omega }_i\\ {\dot{v}}_i=\frac{1}{{\tau }_v}\left(v_{ci}-v_i\right)\\ {\dot{\theta }}_i=\frac{1}{{\tau }_{\theta }}\left({\theta }_{ci}-{\theta }_i\right)\\ {\ddot{z}}_i=-\frac{1}{{\tau }_{\dot{z}}}{\dot{z}}_i+\frac{1}{{\tau }_z}\left(z_{ci}-z_i\right)\end{array}\right.$$

(3)

where $v_i$ is the velocity of the aircraft on the $XOY$-plane; ${\theta }_i$ is the heading; ${\omega }_i$ is the course angular velocity; ${\dot{z}}_i$ is the climb rate; ${\ddot{z}}_i$ is the climb acceleration; ${\tau }_v$ is the speed corresponding to the flight state constant; ${\tau }_{\theta }$ is the flight state constant corresponding to the heading angle; $v_{ci}$ is the speed reference input for the UAV autopilot; ${\theta }_{ci}$ is the course angle reference input of the UAV autopilot; and $z_{ci}$ is the altitude reference input for the UAV autopilot.

In Equation (3), the relationship of ${\theta }_i$,$v_i$ and the velocity component along the OX-axis and the OY-axis is:

$$\left\{\begin{array}{c}\tan {\theta }_i=\frac{v_{yi}}{v_{xi}}\\ v_i=\sqrt{v_{xi}{}^2+v_{yi}{}^2}\end{array}\right.$$

(4)

where $v_{xi}$,$v_{yi}$ are the velocity component.

The dynamic equation with the autonomous driver can be converted into Equation (5):

$$\left\{\begin{array}{l}{\dot{x}}_i=v_{xi}\\ {\dot{y}}_i=v_{yi}\\ {\dot{z}}_i=v_{zi}\\ {\dot{v}}_{xi}=\frac{1}{{\tau }_v}(v_{xi}^c-v_{xi})\\ {\dot{v}}_{yi}=\frac{1}{{\tau }_v}(v_{yi}^c-v_{yi})\\ {\ddot{z}}_i=-\frac{1}{{\tau }_{\dot{z}}}{\dot{z}}_i+\frac{1}{{\tau }_z}\left(z_i^c-z_i\right)\end{array}\right.$$

(5)

where $v_{zi}$ is the speed of the drone along the OZ-axis.

The speed, acceleration, and heading angular velocity of the UAV must be changed within a certain range:

$$\left\{\begin{array}{c}{v}_i\in \left(v_{\min },v_{\max }\right)\\ {\dot{v}}_i\in \left(a_{\min },a_{\max }\right)\\ {\dot{\theta }}_i\in \left({\omega }_{\min },{\omega }_{\max }\right)\\ {\dot{z}}_i\in \left({\dot{z}}_{\min },{\dot{z}}_{\max }\right)\\ {\ddot{z}}_i\in \left({\ddot{z}}_{\min },{\ddot{z}}_{\max }\right)\end{array}\right.$$

(6)

2.3. The Basic Principle of Consensus Algorithm

For any vehicle, its motion states are described by differential equations:

$$\left\{\begin{array}{c}{\dot{\xi }}_i(t)={\zeta }_i(t)\\ {\dot{\zeta }}_i(t)=u_i(t)\end{array}\right.$$

(7)

where ${\xi }_i\in {\mathrm{R}}^{\mathrm{n}}$ is the coordinate vector of the drone; ${\zeta }_i\in {\mathrm{R}}^{\mathrm{n}}$ is the speed vector; $u_i(t)\in {\mathrm{R}}^{\mathrm{n}}$ is the control input vector.

In [30], the basic consensus algorithm given by Equation (7) is:

$$u_i(t)=-{\displaystyle \sum _{j=1}^n{a}_{ij}\left[\left({\xi }_i(t)-{\xi }_j(t)\right)+\alpha \left({\zeta }_i(t)-{\zeta }_j(t)\right)\right]}$$

(8)

where $\alpha >0$; $a_{ij}$ is an element of the matrix $A_n$; matrix $A_n$ is the communication topology; and $i,j$ are two different voyages. If UAV# j can send messages to UAV# $i$, then $a_{ij}=1$, else $a_{ij}=0$.

For UAVs, the information exchange topology is $G_n$, the element $l_{ij}$ in the Laplace matrix $L_n$ is defined as:

$$l_{ij}=\left\{\begin{array}{l}-a_{ij}\begin{array}{cc}\begin{array}{cc}& \end{array}& \begin{array}{cc}if& i\ne j\end{array}\end{array}\\ {\displaystyle \sum _{j=1,i\ne j,}^n{a}_{ij}\begin{array}{cc}& \begin{array}{cc}if& i=j\end{array}\end{array}}\end{array}\right.$$

(9)

We assume that the topology of ${\xi }_i,{\zeta }_i$ is consistent during the flight of the UAVs; then, the condition for the consensus algorithm to converge is:

Lemma 1 [31]

If $G_n$ has a directed spanning tree, $\alpha >\overline{\alpha }$ , the state of the UAV formation can be asymptotically consistent. If the $n-1$ non-zero eigenvalues of $-L_n$ are negative, then$\overline{\alpha }=0$ , otherwise:

$$\overline{\alpha }=\underset{\forall \mathrm{Im}\left({\eta }_i\right)>0,\mathrm{Re}\left({\eta }_i\right)<0}{\max }\sqrt{\frac{2}{\left|{\eta }_i\right|\cos \left(\arctan \frac{\mathrm{Im}\left({\eta }_i\right)}{-\mathrm{Re}\left({\eta }_i\right)}\right)}}$$

(10)

The lemma shows that, when the state is able to converge according to the consensus algorithm, then, for any initial state such as $x_i\left(0\right)$ and $v_i\left(0\right)$, when $t\to \infty$, there are $\left|x_i\left(t\right)-x_j\left(t\right)\right|\to 0$ and $\left|v_i\left(t\right)-v_j\left(t\right)\right|\to 0$.

It is necessary to set a reasonable communication topology and $\alpha$ value so that the state of the drone converges to the same level.

3. Improved Consistency Algorithm

The above Equation (8) does not consider the network communication delay and network topology switching in the formation flight of UAVs, nor does it consider the constraints of UAV maneuverability and flight performance; therefore, the consensus algorithm needs to be improved for the actual flight of UAVs.

Firstly, for the situation of non-symmetrical communication delay and topology switching in formation flight, we design the consensus control protocol for the formation flight. Then, the designed control protocol is modified to make itself and the corresponding state output meet the constraints of UAV maneuvering and flight performance.

3.1. Formation State Control

This section studies the consensus control protocol in the case of joint connectivity communication topology.

The state of the dynamic equation of UAV#$i$ is shown in Equation (7).

If the formation protocol can ensure that the states of the UAVs meet the conditions: $\left[{\xi }_i-{\xi }_j\right]\to r_{ij}$ and ${\zeta }_i\to {\zeta }_i\to {\zeta }^{\ast }$, ($r_{ij}$ is the expected difference in position between two drones, and $r_{ij}=-r_{ji}$, ${\zeta }^{\ast }$ is the desired speed vector). This shows that the control algorithm can make the multiple UAVs form our expected formation and move forward according to the expected flight speed finally.

The work of [31] gives a control protocol that can make the multi-aircraft system form the desired formation and achieve a given speed, but only for fixed communication topology, and does not consider the communication delay of the system.

This section refers to the control protocol idea of [15] for the multi-aircraft formation flight control system with non-symmetrical communication delay and a communication topology map that is jointly connected. The formation protocol for the UAV is Equation (11):

$$\begin{array}{l}{u}_i(t)={\displaystyle \sum _{j\in N_i(t)}{a}_{ij}(t)\left\{k_1\left[{\xi }_j(t-{\tau }_{ij})-{\xi }_i(t-{\tau }_{ii})-r_{ji}\right]+k_2\left[{\zeta }_j(t-{\tau }_{ij})-{\zeta }_i(t-{\tau }_{ii})\right]\right\}}\\ \begin{array}{cc}& \end{array}+{\dot{\zeta }}^{*}-k_3\left({\zeta }_i(t)-{\zeta }^{\ast }\right)\end{array}$$

(11)

$u_{xi},u_{yi},u_{zi}$ are shown in Equations (12)–(14), as follows:

$$\left\{\begin{array}{l}{v}_{xci}=v_{xi}+{\tau }_v{u}_{xi}\\ u_{xi}={\displaystyle \sum _{V_j\in N_i\left(t\right)}{a}_{ij}\left(t\right)\left\{k_1\left[x_j-x_i-{\overline{x}}_{ji}\right]+\frac{2}{k_2}\left[v_{xj}-v_{xi}\right]\right\}+{\dot{v}}^{*}{}_x-k_3^x\left(v_{xi}-v_x^{\ast }\right)}\end{array}\right.$$

(12)

$$\left\{\begin{array}{l}{v}_{yci}=v_{yi}+{\tau }_v{u}_{yi}\\ u_{yi}={\displaystyle \sum _{V_j\in N_i\left(t\right)}{a}_{ij}\left(t\right)\left\{k_1\left[y_j-y_i-{\overline{y}}_{ji}\right]+\frac{2}{k_2}\left[v_{yj}-v_{yi}\right]\right\}+{\dot{v}}^{*}{}_y-k_3\left(v_{yi}-v_y^{\ast }\right)}\end{array}\right.$$

(13)

$$\left\{\begin{array}{l}{z}_i^c=z_i+\frac{{\tau }_z}{{\tau }_{\dot{z}}}{\dot{z}}_i+{\tau }_z{u}_{zi}\\ u_{zi}={\displaystyle \sum _{V_j\in N_i\left(t\right)}{a}_{ij}\left(t\right)\left\{k_1\left[z_j-z_i-{\overline{z}}_{ji}\right]+\frac{2}{k_2}\left[v_{zj}-v_{zi}\right]\right\}+{\dot{v}}^{*}{}_z-k_3\left(v_{yi}-v_y^{\ast }\right)}\end{array}\right.$$

(14)

where ${\tau }_{ii}$ is the time change in the UAV#i itself, a type of latency which is caused by measurements or calculations; ${\tau }_{ij}$ represents the time delay for UAV#j to receive the state information from UAV#i; $N_i(t)$ is a collection of neighbors of node $i$, and $k_1,k_2,k_3>0,k_3=k_1{k}_2$.

Suppose there are M numbers of different time delay in total, it is expressed as ${\tau }_m\left(t\right)\in \left\{{\tau }_{ii}\left(t\right),{\tau }_{ij}\left(t\right),i,j\in \mathcal{l}\right\},m=1,2,3,\cdots M$, and satisfies Assumption 1.

Assumption 1. 

For specific normal values $h_m>0,d_m>0$, time-varying delay time ${\tau }_m\left(t\right),m=1,2,3,\cdots ,M$ satisfies $0\le {\tau }_m\le h_m$ and ${\dot{\tau }}_m\le d_m\le 1$.

When the network topology is switched and there is a delay in communication, this control protocol can realize the coordinated flight of multiple UAVs.

3.2. Formation Control Protocol Adjustment under Constraints

Section 3.1 does not consider the constraints of Equation (6) when designing the formation control protocol, so the generated control commands and corresponding flight states may not meet the requirements of UAV maneuverability and flight performance. This section proposes a strategy called minimum adjustment to adjust the formation control protocol in Section 2.1 so that both itself and the corresponding state output meet the constraints of UAV maneuvering and flight performance.

The control command $u_{xi},u_{yi}$ is adjusted in the $XOY$-plane, so that it satisfies the constraints of velocity $v_i$, acceleration ${\dot{v}}_i$, and heading angular velocity ${\dot{\phi }}_i$. Then, the values of $u_{xi},u_{yi}$ are fixed and the value of the control instruction $u_{zi}$ is adjusted to meet the constraints of the OZ-axis direction of the climbing speed ${\dot{z}}_i$ and the climbing acceleration ${\ddot{z}}_i$.

Then, $u_{xi},u_{yi}$ are adjusted in Equations (12) and (13), then the related constraints of speed $v_i$, acceleration ${\dot{v}}_i$, and heading angular velocity ${\dot{\phi }}_i$, $v_i\left(t+\Delta t\right)$ can be obtained through the current flight status:

$$\left\{\begin{array}{l}{\alpha }_i\left(t\right)=\sqrt{u_{xi}^2\left(t\right)+u_{yi}^2\left(t\right)}\\ v_i\left(t+\Delta t\right)=v_i\left(t\right)+{\alpha }_i\left(t\right)\Delta t\end{array}\right.$$

(15)

If $v_i\left(t+\Delta t\right)$ does not satisfy the constraint $v_i\left(t+\Delta t\right)\in \left(v_{\min },v_{\max }\right)$, the following variables can be defined as follows:

$$\left\{\begin{array}{c}{\alpha }_{\min ,i}^{’}\left(t\right)=\frac{v_{\min }-v_i(t)}{\Delta t}\\ {\alpha }_{\max ,i}^{’}\left(t\right)=\frac{v_{\max }-v_i(t)}{\Delta t}\end{array}\right.$$

(16)

where ${\alpha }_{\min ,i}^{’}\left(t\right)$ and ${\alpha }_{\max ,i}^{’}\left(t\right)$ are the accelerations of the UAV#$i$ at time $t$.

When the speeds are $v_{\min },v_{max}$ at $t+\Delta t$, ${\alpha }_i\left(t\right)\in \left[{\alpha }_{\min ,i}^{’}\left(t\right),{\alpha }_{\max ,i}^{’}\left(t\right)\right]$. $a_{\min ,i}^{’}\left(t\right)$ is compared with $a_{\min }$ and $a_{\max ,i}^{’}\left(t\right)$ is compared with $a_{\max }$, respectively, to obtain the updated constraints of acceleration:

$$\left\{\begin{array}{c}{a}_{\min ,i}^{new}\left(t\right)=\max \left(a_{\min },a_{\min ,i}^{’}\left(t\right)\right)\\ a_{\max ,i}^{new}\left(t\right)=\min \left(a_{\max },a_{\max ,i}^{’}\left(t\right)\right)\end{array}\right.$$

(17)

where $a_{\max ,i}^{new}\left(t\right)$ is the upper limit of $a_i\left(t\right)$, and $a_{\min ,i}^{new}\left(t\right)$ is the lower limit of $a_i\left(t\right)$.

Equation (17) actually includes constraints on the $v_i\left(t+\Delta t\right)$. As long as the acceleration $a_i\left(t\right)$ of the UAV# i satisfies Equation (17), the two constraints on the acceleration and velocity in the $XOY$-plane can be satisfied at the same time. If $v_i\left(t+\Delta t\right)\in \left(v_{\min },v_{\max }\right)$, then the values of $a_{\min }$ and $a_{\max }$ do not need to be updated by Equation (16). Constraint $a_i\in \left[a_{\min ,i}^{new},a_{\max ,i}^{new}\right]$ is used to adjust $u_{xi},u_{yi}$ so that it satisfies the constraints.

The adjustment of $u_{xi}\left(t\right),u_{yi}\left(t\right)$ needs to be carried out synchronously, and the influence on the original acceleration value should be as small as possible. The following adjustments can be made to $u_{xi}\left(t\right),u_{yi}\left(t\right)$:

$$\left\{\begin{array}{c}{u}_{xi}^{’}\left(t\right)=a_{\max ,i}^{new}\left(t\right)\frac{u_{xi}\left(t\right)}{a_i\left(t\right)},u_{yi}^{’}\left(t\right)=a_{\max ,i}^{new}\left(t\right)\frac{u_{yi}\left(t\right)}{a_i\left(t\right)}\\ \begin{array}{l}if\begin{array}{cc}& a_i\left(t\right)>a_{\max ,i}^{new}\left(t\right)\end{array}\\ u_{xi}^{’}\left(t\right)=a_{\min ,i}^{new}\left(t\right)\frac{u_{xi}\left(t\right)}{a_i\left(t\right)},u_{yi}^{’}\left(t\right)=a_{\min ,i}^{new}\left(t\right)\frac{u_{yi}\left(t\right)}{a_i\left(t\right)}\\ if\begin{array}{cc}& a_i\left(t\right)<a_{\min ,i}^{new}\left(t\right)\end{array}\end{array}\end{array}\right.$$

(18)

Then, the values of $u_{xi}^{’},u_{yi}^{’}$ meet the constraints of acceleration and speed after the above adjustments, and the adjustment range is the smallest.

After that, the heading angular velocity ${\dot{\theta }}_i\left(t\right)$ constraint is processed, and $u_{xi}^{’}\left(t\right)$$u_{yi}^{’}\left(t\right)$ will be adjusted in the next step.

According to the constraints of ${\dot{\theta }}_i\in \left({\omega }_{\min },{\omega }_{\max }\right)$, the allowable value range of the heading angle ${\theta }_i$ at the next sampling time can be obtained as:

$$\left\{\begin{array}{c}{\theta }_{\min ,i}\left(t+\Delta t\right)={\theta }_i\left(t\right)+{\omega }_{\min }\Delta t\\ {\theta }_{\max ,i}\left(t+\Delta t\right)={\theta }_i\left(t\right)+{\omega }_{\max }\Delta t\end{array}\right.$$

(19)

The heading angle ${\theta }_i\left(t+\Delta t\right)$ at the next sampling time is:

$${\theta }_i\left(t+\Delta t\right)=\arctan \frac{v_{yi}\left(t\right)+u_{yi}^{’}\left(t\right)\Delta t}{v_{xi}\left(t\right)+u_{xi}^{’}\left(t\right)\Delta t}$$

(20)

where $v_{xi}\left(t\right)$ and $v_{yi}\left(t\right)$ are the speeds of the drone at time $t$.

If ${\theta }_i\left(t+\Delta t\right)\notin \left[{\theta }_{\min ,i},{\theta }_{\max ,i}\right]$, then $u_{xi}^{’}\left(t\right),u_{yi}^{’}\left(t\right)$ should be adjusted by Equations (21) and (22).

$$\left\{\begin{array}{l}\frac{v_{yi}\left(t\right)+u_{yi}^{’’}\left(t\right)\Delta t}{v_{xi}\left(t\right)+u_{xi}^{’’}\left(t\right)\Delta t}=\tan \left({\theta }_{\max ,i}\left(t+\Delta t\right)\right)\\ u_{xi}^{’’2}\left(t\right)+u_{yi}^{’’2}\left(t\right)=a_i^{’2}\left(t\right)\\ {\theta }_i\left(t+\Delta t\right)>{\theta }_{\max ,i}(t+\Delta t)\end{array}\right.$$

(21)

$$\left\{\begin{array}{l}\frac{v_{yi}\left(t\right)+u_{yi}^{’’}\left(t\right)\Delta t}{v_{xi}\left(t\right)+u_{xi}^{’’}\left(t\right)\Delta t}=\tan \left({\theta }_{\min ,i}\left(t+\Delta t\right)\right)\\ u_{xi}^{’’2}\left(t\right)+u_{yi}^{’’2}\left(t\right)=a_i^{’2}\left(t\right)\\ {\theta }_i\left(t+\Delta t\right)<{\theta }_{\min ,i}(t+\Delta t)\end{array}\right.$$

(22)

where Equations (21) and (22) are binary quadratic equations; usually, there are two different sets of solutions, denoted as $u_{xi1}^{’’}\left(t\right),u_{yi1}^{’’}\left(t\right)$ and $u_{xi2}^{’’}\left(t\right),u_{yi2}^{’’}\left(t\right)$. Because both sets of solutions satisfy the constraints of the heading angular velocity ${\dot{\phi }}_i$, it is necessary to further confirm which set is finally selected as the result according to the “minimum adjustment” strategy.

Let ${\gamma }_{ai}^{’}\left(t\right)=\arctan \frac{u_{yi}^{’}\left(t\right)}{u_{xi}^{’}\left(t\right)}$ represent the direction of the acceleration $a_i^{’}\left(t\right)$; the adjusted values should not only keep the value of $a_i^{’}\left(t\right)$ unchanged, but also the direction of $a_i^{’}\left(t\right)$ should change minimally.

${\gamma }_{ai1}^{’’}\left(t\right)=\arctan \frac{u_{yi1}^{’’}\left(t\right)}{u_{xi1}^{’’}\left(t\right)}$ and ${\gamma }_{ai2}^{’’}\left(t\right)=\arctan \frac{u_{yi2}^{’’}\left(t\right)}{u_{xi2}^{’’}\left(t\right)}$ represent the directions of the UAV#$i$‘s acceleration in the $XOY$-plane after adjustment by Equations (21) or (22).

Among these two sets of solutions, the set of solutions corresponding to $\min \left(\left|{\gamma }_{ai1}^{’’}\left(t\right)-{\gamma }_{ai}^{’}\left(t\right)\right|,\left|{\gamma }_{ai2}^{’’}\left(t\right)-{\gamma }_{ai}^{’}\left(t\right)\right|\right)$ is selected as the values required.

The above procedure makes minimal adjustments to $u_{xi}\left(t\right),u_{yi}\left(t\right)$ and satisfies the constraints of the $XOY$-plane.

For the constraints in the direction of the OZ-axis, the climbing rate ${\dot{z}}_i\left(t\right)$ and climb acceleration ${\ddot{z}}_i\left(t\right)$ of the drone are limited and $u_{zi}\left(t\right)$ is adjusted. First, the climbing rate of the next sampling moment ${\dot{z}}_i\left(t+\Delta t\right)$ is calculated through the current time of the drone flight status as ${\dot{z}}_i\left(t+\Delta t\right)={\dot{z}}_i\left(t\right)+u_{zi}\left(t\right)\Delta t$.

If ${\dot{z}}_i\left(t+\Delta t\right)\notin \left[{\dot{z}}_{\min },{\dot{z}}_{\max }\right]$, then the updating constraints of the climbing acceleration ${\ddot{z}}_i\left(t\right)$ are shown in the following formula:

$$\left\{\begin{array}{c}{\ddot{z}}_{\min ,i}^{’}\left(t\right)=\frac{{\dot{z}}_{\min }-{\dot{z}}_i\left(t\right)}{\Delta t}\\ {\ddot{z}}_{\max ,i}^{’}\left(t\right)=\frac{{\dot{z}}_{\max }-{\dot{z}}_i\left(t\right)}{\Delta t}\end{array}\right.$$

(23)

where ${\ddot{z}}_{\min ,i}^{’}\left(t\right)$ and ${\ddot{z}}_{\max ,i}^{’}\left(t\right)$ are the lower limit and upper limit of the constraints after the climbing rate constraints are converted to the climbing acceleration at time $t$, respectively. If ${\dot{z}}_{\min }\le {\dot{z}}_i\left(t+\Delta t\right)\le {\dot{z}}_{\max }$, then ${\ddot{z}}_{\min ,i}^{’}\left(t\right)$ and ${\ddot{z}}_{\max ,i}^{’}\left(t\right)$ do not need to be updated.

Through Equation (24), the constraints on the climbing rate are also converted into the constraint on the climbing acceleration at time $t$. Then, ${\ddot{z}}_{\min ,i}^{’}\left(t\right)$, ${\ddot{z}}_{\min }$, ${\ddot{z}}_{\max ,i}^{’}\left(t\right)$, and ${\ddot{z}}_{\max }$ are compared, and the updated rising acceleration constraint conditions is determined as:

$$\left\{\begin{array}{l}{\ddot{z}}_{\min ,i}^{new}\left(t\right)=\max \left({\ddot{z}}_{\min },{\ddot{z}}_{\min ,i}^{’}\left(t\right)\right)\\ {\ddot{z}}_{\max ,i}^{new}\left(t\right)=\min \left({\ddot{z}}_{\max },{\ddot{z}}_{\max ,i}^{’}\left(t\right)\right)\end{array}\right.$$

(24)

where ${\ddot{z}}_{\min ,i}^{new}\left(t\right)$ and ${\ddot{z}}_{\max ,i}^{new}\left(t\right)$ are the final lower limit and upper limit values of the climbing acceleration, respectively, after the climbing rate and climbing acceleration constraints have been considered.

Finally, we limit the current climb acceleration $u_{zi}\left(t\right)$ to the allowable range:

$$u_{zi}^{’}\left(t\right)=\max \left({\ddot{z}}_{\min ,i}^{new}\left(t\right),\min \left(u_{zi}\left(t\right),{\ddot{z}}_{\max ,i}^{new}\left(t\right)\right)\right)$$

(25)

where $u_{zi}^{’}\left(t\right)$ is the adjusted climbing acceleration. From Equation (25), when $u_{zi}\left(t\right)\in \left[{\ddot{z}}_{\min ,i}^{new}\left(t\right),{\ddot{z}}_{\max ,i}^{new}\left(t\right)\right]$, we have $u_{zi}^{’}\left(t\right)=u_{zi}\left(t\right)$; when $u_{zi}\left(t\right)<{\ddot{z}}_{\min ,i}^{new}\left(t\right)$, we have $u_{zi}^{’}\left(t\right)={\ddot{z}}_{\min ,i}^{new}\left(t\right)$; when $u_{zi}\left(t\right)>{\ddot{z}}_{\max ,i}^{new}\left(t\right)$, we have $u_{zi}^{’}\left(t\right)={\ddot{z}}_{\max ,i}^{new}\left(t\right)$.

3.3. Convergence Proof of Improved Consistency Algorithm

Let ${\overline{\xi }}_i={\xi }_i-{\xi }_0-r_i,{\overline{\zeta }}_i={\zeta }_i-{\zeta }^{\ast }$, then Equation (11) can be transformed into:

$$\begin{array}{l}{u}_i\left(t\right)={\displaystyle \sum _{s_j\in N_i\left(t\right)}{a}_{ij}\left(t\right)}\left\{k_1\left[({\overline{\xi }}_j\left(t-{\tau }_{ij}\left(t\right)\right)-({\overline{\xi }}_i\left(t-{\tau }_{ii}\left(t\right)\right)\right]\right.\\ \begin{array}{cc}& +\frac{2}{k_2}\left[{\overline{\zeta }}_j\left(t-{\tau }_{ij}\left(t\right)\right)-{\overline{\zeta }}_i\left(t-{\tau }_{ii}\left(t\right)\right)\right]\end{array}+{\dot{\zeta }}^{\ast }-k_3{\overline{\zeta }}_i\left(t\right)\end{array}$$

(26)

If ${\widehat{\zeta }}_i(t)=2{\overline{\zeta }}_i(t)/k_1{k}_2+{\overline{\xi }}_i(t),\epsilon (t)={[{\overline{\xi }}_1\left(t\right),{\widehat{\zeta }}_i(t),\cdots ,{\overline{\xi }}_n\left(t\right),{\widehat{\zeta }}_n(t)]}^T$, then:

$$B=\left[\begin{array}{cc}-k_3/2& k_3/2\\ k_3/2& k_3/2\end{array}\right],Q=\left[\begin{array}{cc}0& 0\\ 0& 2/k_2\end{array}\right]$$

According to Equation (26), the closed-loop dynamic equation is Equation (27):

$$\dot{\epsilon }(t)=(I_n\otimes B)\epsilon (t)-{\displaystyle {\sum }_{m=1}^M(L_{\sigma m}\otimes Q)\epsilon (t-{\tau }_m)}$$

(27)

In fact, if we have $\underset{t\to +\infty }{\lim }\epsilon (t)=0$, then $\underset{t\to +\infty }{\lim }{\xi }_j(t)-{\xi }_i(t)=r_{ji},\underset{t\to +\infty }{\lim }\zeta (t)={\zeta }^{\ast }$.

Next, we show that the above closed-loop control system can achieve $\underset{t\to +\infty }{\lim }\epsilon \left(t\right)=0$.

Referring to the definition of switching topology, it is assumed that the time-invariant topology ${\overline{G}}_{\sigma }$ in a certain sub-interval $[t_{k_b},t_{k_{b+1}})$ has $q(q\ge 1)$ the numbers of connected parts, and its corresponding node set is denoted by ${\psi }_{k_j}^1,{\psi }_{k_j}^2,\cdots ,{\psi }_{k_j}^{d_{\sigma }}$, and $f$ is the node number in ${\psi }_{k_j}^i$. Then, a permutation matrix $P_{\sigma }\in {ℝ}^{n\times n}$ is obtained by satisfying:

$$\begin{array}{l}{P}_{\sigma }^T{L}_{\sigma }{P}_{\sigma }=diag\left\{L_{\sigma }^1,L_{\sigma }^2,\cdots ,L_{\sigma }^{d_{\sigma }}\right\}\\ P_{\sigma }^T{L}_{\sigma m}{P}_{\sigma }=diag\left\{L_{\sigma m}^1,L_{\sigma m}^2,\cdots ,L_{\sigma m}^q\right\}\end{array}$$

(28)

$${\epsilon }^T\left(t\right)\left(P_{\sigma }\otimes I_2\right)=\left[{\epsilon }_{\sigma }^{\begin{array}{cc}1& T\end{array}},{\epsilon }_{\sigma }^{\begin{array}{cc}2& T\end{array}},\cdots ,{\epsilon }_{\sigma }^{\begin{array}{cc}q& T\end{array}}\right]$$

(29)

where $L_{\sigma }^i\in {ℝ}^{f_{\sigma }^i\times f_{\sigma }^i}$ is the Laplacian matrix which corresponds to the part which is connected, and $L_{\sigma m}^i\in {ℝ}^{f\times f},L_{\sigma }^i={\displaystyle {\sum }_{m=1}^M{L}_{\sigma m}^i}$. Therefore, in $[t_{k_b},t_{k_{b+1}})$, it can be broken down into q numbers of subsystems:

$${\dot{\epsilon }}_{\sigma }^i(t)=(I_f\otimes B){\epsilon }_{\sigma }^i(t)-{\displaystyle {\sum }_{m=1}^M(L_{\sigma m}^i\otimes Q){\epsilon }_{\sigma }^i(t-{\tau }_m)},i=1,2,\cdots ,q$$

(30)

where ${\epsilon }_{\sigma }^i(t)=[{\epsilon }_{\sigma 1}^i(t),\cdots {\epsilon }_{\sigma 2f_{\sigma }^i}^i(t)]\in {ℝ}^{2f_{\sigma }^i}$.

Lemma 2. [32]

If there is $D_n=n{I}_n-11^T$, then there must be an orthogonal matrix $U_n\in {ℝ}^{n\times n}$which makes $U_n^{T}D{U}_n=diag\left\{n{I}_{n-1},0\right\}$, where the last column $U_n$ is $1/\sqrt{n}$ . We give a matrix $D\in {ℝ}^{n\times n}$ and make it satisfy $1^{T}D=0$ and $D1=0$ , then $U_n^{T}D{U}_n=diag\left\{{\overline{U}}_n^{T}D{\overline{U}}_n,0\right\}$.

Lemma 3.

[33] For any function of actual cable vector $x\left(t\right)\in R^n$, any function of cable scalar $\tau \left(t\right)\in \left[0,a\right]$, and any constant matrix $0<H=H^T\in R^{n\times n}$, there is:

$$\frac{1}{a}{\left[x\left(t\right)-x\left(t-\tau \left(t\right)\right)\right]}^{T}H\left[x\left(t\right)-x\left(t-\tau \left(t\right)\right)\right]\le {\displaystyle {\int }_{t-\tau \left(t\right)}^t{\dot{x}}^T\left(s\right)H\dot{x}\left(s\right)}ds,t\ge 0$$

(31)

where $a>0$.

Theorem 1.

Considering a multi -UAV system with non-uniform time delay and switching topology, in any sub -interval$\left[t_{r_b},t_{r_{b+1}}\right)$, if $\mathsf{\gamma }>0$, and $F_{\sigma }^i\in {ℝ}^{f\times f},i=1,2,\cdots ,q$, there is:

$$F_{\sigma }^{\begin{array}{cc}i& T\end{array}}{\mathrm{\Xi }}_{\sigma }^i{F}_{\sigma }^i<0$$

(32)

then there is $\underset{t\to \infty }{\lim }{\xi }_j\left(t\right)-{\xi }_i\left(t\right)=r_{ji},\underset{t\to \infty }{\lim }{\zeta }_i\left(t\right)={\zeta }^{\ast }$.

where $F_{\sigma }^i=diag\left\{U_{2f},I_{2Mf}\right\}$, and the definition of $U_{2f}$ is as shown in Lemma 2.

where ${\mathrm{\Xi }}_{\sigma }^i=\left[\begin{array}{cc}{\mathrm{\Xi }}_{11}& {\mathrm{\Xi }}_{12}\\ {\mathrm{\Xi }}_{12}^T& {\mathrm{\Xi }}_{22}\end{array}\right]$, and

$${\mathrm{\Xi }}_{11}=2\mathsf{\gamma }\left(I_f\otimes B\right)+{\displaystyle {\sum }_{m=1}^M{\mathrm{h}}_{\mathrm{m}}{\left(I_f\otimes B\right)}^T}\left(I_f\otimes B\right)-{\displaystyle {\sum }_{m=1}^M\frac{1-{\mathrm{d}}_{\mathrm{m}}}{{\mathrm{h}}_{\mathrm{m}}}}{I}_{2f}$$

$$\begin{array}{l}{\mathrm{\Xi }}_{12}=\left[-\mathsf{\gamma }\left(L_{\sigma 1}^i\otimes Q\right)+\frac{1-d_1}{h_1}{I}_{2f}-{\displaystyle {\sum }_{m=1}^M{h}_m{\left(I_f\otimes B\right)}^T\left(L_{\sigma 1}^i\otimes Q\right),\cdots }\right.\\ \left.\begin{array}{cc}& -\mathsf{\gamma }\left(L_{\sigma M}^i\otimes Q\right)+\frac{1-d_M}{h_M}{I}_{2f}-{\displaystyle {\sum }_{m=1}^M{h}_m{\left(I_f\otimes B\right)}^T\left(L_{\sigma M}^i\otimes Q\right)}\end{array}\right]\end{array}$$

$$\begin{array}{l}{\mathrm{\Xi }}_{22}=\left[-diag\left\{\frac{1-d_1}{h_1}{I}_{2f},\frac{1-d_2}{h_2}{I}_{2f},\cdots \frac{1-d_M}{h_M}{I}_{2f}\right\}+\right.\\ \begin{array}{cc}& {\displaystyle {\sum }_{m=1}^M{h}_m}\end{array}{\left[\left(L_{\sigma 1}^i\otimes Q\right),\cdots \left(L_{\sigma M}^i\otimes Q\right)\right]}^T\left[\left(L_{\sigma 1}^i\otimes Q\right),\cdots \left(L_{\sigma M}^i\otimes Q\right)\right]\end{array}$$

Proof of Theorem 1.

The Lyapunov−Krasovskii function for Equation (11) can be defined following:

$$V\left(t\right)=\mathsf{\gamma }{\epsilon }^T\left(t\right)\epsilon \left(t\right)+{\displaystyle {\sum }_{m=1}^M{\displaystyle {\int }_{-{\tau }_m}^0{\displaystyle {\int }_{t+a}^t{\dot{\epsilon }}^T\left(s\right)}}}\dot{\epsilon }\left(s\right)dsda,\mathsf{\gamma }>0$$

(33)

$\dot{V}\left(t\right)$ can be calculated as:

$$\begin{array}{l}\dot{V}\left(t\right)=2\mathsf{\gamma }{\epsilon }^T\left(t\right)\dot{\epsilon }\left(t\right)+{\displaystyle {\sum }_{m=1}^M{\tau }_m}{\dot{\epsilon }}^T\left(t\right)\dot{\epsilon }\left(t\right)-{\displaystyle {\sum }_{m=1}^M\left(1-{\dot{\tau }}_m\right)}{\displaystyle {\int }_{t-{\tau }_m}^t{\dot{\epsilon }}^T\left(s\right)\dot{\epsilon }\left(s\right)}ds\\ \begin{array}{cc}& \begin{array}{cc}& =\end{array}\end{array}2\mathsf{\gamma }{\epsilon }^T\left(t\right)\left[\left(I_f\otimes B\right)\epsilon \left(t\right)\right]-2\mathsf{\gamma }{\epsilon }^T\left(t\right){\displaystyle {\sum }_{m=1}^M\left[\left(L_{\sigma m}^i\otimes Q\right)\epsilon \left(t-{\tau }_m\right)\right]}+\\ \begin{array}{ccc}& & {\displaystyle {\sum }_{m=1}^M{\tau }_m}{\dot{\epsilon }}^T\left(t\right)\dot{\epsilon }\left(t\right)-\end{array}{\displaystyle {\sum }_{m=1}^M\left(1-{\dot{\tau }}_m\right)}{\displaystyle {\int }_{t-{\tau }_m}^t{\dot{\epsilon }}^T\left(s\right)\dot{\epsilon }\left(s\right)}ds\end{array}$$

(34)

According to Equation (33) and Assumption 1, $\dot{V}\left(t\right)$ is changed to the following form:

$$\begin{array}{l}\dot{V}\left(t\right)\le {\displaystyle \sum _{i=1}^q\left\{2\mathsf{\gamma }{\epsilon }_{\sigma }^{\begin{array}{cc}i& T\end{array}}\left(t\right)\left[\left(I_f\otimes B\right){\epsilon }_{\sigma }^i\left(t\right)-2\mathsf{\gamma }{\epsilon }_{\sigma }^{\begin{array}{cc}i& T\end{array}}\left(t\right){\displaystyle {\sum }_{m=1}^M\left[\left(L_{\sigma m}^i\otimes Q\right){\epsilon }_{\sigma }^i\left(t-{\tau }_m\right)\right]}\right]\right.}\\ \left.\begin{array}{cc}& +\end{array}{\displaystyle {\sum }_{m=1}^M{h}_m}{\dot{\epsilon }}_{\sigma }^{\begin{array}{cc}i& T\end{array}}\left(t\right){\dot{\epsilon }}_{\sigma }^i\left(t\right)-{\displaystyle {\sum }_{m=1}^M\left(1-d_m\right){\displaystyle {\int }_{t-{\tau }_m}^t{\dot{\epsilon }}^T\left(s\right)\dot{\epsilon }\left(s\right)}ds}\right\}\end{array}$$

(35)

According to Lemma 3, the following can be obtained:

$$\begin{array}{l}\dot{V}(t)\le \sum _{i=1}^q \left\{2\gamma {\epsilon }_{\sigma }^i{}^T(t)\left[\left(I_f\otimes B\right){\epsilon }_{\sigma }^i(t)-2\gamma {\epsilon }_{\sigma }^{i\mathrm{T}}(t)\sum _{m=1}^M \left[\left(L_{\sigma m}^i\otimes Q\right){\epsilon }_{\sigma }^i\left(t-{\tau }_m\right)\right]\right]+\sum _{m=1}^M h_m{\dot{\epsilon }}_{\sigma }^{iT}(t){\dot{\epsilon }}_{\sigma }^i(t)\right.\\ -\sum _{m=1}^M \frac{1-d_m}{h_m}\left[{\dot{\epsilon }}_{\sigma }^i{}^{\mathrm{T}}(t){\dot{\epsilon }}_{\sigma }^i(t)-{\dot{\epsilon }}_{\sigma }^{i\mathrm{T}}(t){\dot{\epsilon }}_{\sigma }^i\left(t-{\tau }_m\right)-{\epsilon }_{\sigma }^i{}^T\left(t-{\tau }_m\right){\epsilon }_{\sigma }^i(t)+{\epsilon }_{\sigma }^{i\mathrm{T}}\left(t-{\tau }_m\right){\epsilon }_{\sigma }^i\left(t-{\tau }_m\right)\right]\\ =\sum _{i=1}^q {\delta }_i^T{\mathrm{\Xi }}_{\sigma }^i{\delta }_i\end{array}$$

where ${\delta }_i=\left[{\epsilon }_{\sigma }^{\begin{array}{cc}i& T\end{array}}\left(t\right),{\epsilon }_{\sigma 1}^{\begin{array}{cc}i& T\end{array}}\left(t-{\tau }_1\right),{\epsilon }_{\sigma 2}^{\begin{array}{cc}i& T\end{array}}\left(t-{\tau }_2\right),\cdots {\epsilon }_{\sigma M}^{\begin{array}{cc}i& T\end{array}}\left(t-{\tau }_M\right)\right]$.

${\mathrm{\Xi }}_{\sigma }^i={\mathrm{\Xi }}_{\sigma }^{\begin{array}{cc}i& T\end{array}}$, and ${\mathrm{\Xi }}_{\sigma }^i{\left[1_{2f}^T,0_{2Mf}^T\right]}^T=0$. According to Lemma 2, we can conclude ${\mathrm{\Xi }}_{\sigma }^i\le 0$, when $F_{\sigma }^{\begin{array}{cc}i& T\end{array}}{\mathrm{\Xi }}_{\sigma }^i{F}_{\sigma }^i<0$, $rank({\mathrm{\Xi }}_{\sigma }^i)=2\left(M+1\right)f-1$.

Let $\eta =\left[{\epsilon }_{\sigma }^{\begin{array}{cc}i& T\end{array}}\left(t\right)-\mathrm{h}{1}^T,{\epsilon }_{\sigma 1}^{\begin{array}{cc}i& T\end{array}}\left(t\right),{\epsilon }_{\sigma 2}^{\begin{array}{cc}i& T\end{array}}\left(t\right),\cdots {\epsilon }_{\sigma M}^{\begin{array}{cc}i& T\end{array}}\left(t\right)\right]$, $\mathrm{h}>0$, then ${\mathrm{\Xi }}_{\sigma }^i\left({\delta }_i-\eta \right)=0$, and we can obtain

$${\delta }_i^T{\mathrm{\Xi }}_{\sigma }^i{\delta }_i={\eta }^{\mathrm{T}}{\mathrm{\Xi }}_{\sigma }^i\eta \le \lambda {‖\eta ‖}^2\le \lambda \left[{‖{\epsilon }_{\sigma }^i\left(t\right)-\mathrm{h}1‖}^2+{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{K=1}^f{\left({\epsilon }_{\sigma mk}^i\right)}^2\left(t\right)}}\right]$$

(36)

where $\Vert ·\Vert$ is the Standard European norm and λ < 0 represents the maximum non-zero eigenvalue of ${\mathrm{\Xi }}_{\sigma }^i$.

Therefore, $\dot{V}\left(t\right)\le \lambda {\displaystyle \sum _{i=1}^q\left[{‖{\epsilon }_{\sigma }^i\left(t\right)-\mathrm{h}1‖}^2+{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{K=1}^f{\left({\epsilon }_{\sigma mk}^i\right)}^2\left(t\right)}}\right]}\le 0$.

Through the above analysis, Equation (11) is stable, and $\underset{t\to +\infty }{\lim }V\left(t\right)=0$. Then, we can obtain $\underset{t\to +\infty }{\lim }\epsilon \left(t\right)=0$, then $\underset{t\to +\infty }{\lim }{\overline{\xi }}_i\left(t\right)=0,\underset{t\to +\infty }{\lim }{\widehat{\zeta }}_i\left(t\right)=0$, and we can have $\underset{t\to +\infty }{\lim }{\xi }_j\left(t\right)-{\xi }_i\left(t\right)=r_{ji},\underset{t\to +\infty }{\lim }{\zeta }_i\left(t\right)={\zeta }^{\ast }$. That is, under the action of the control protocol of Equation (30), the drones can eventually form the specific formation at an expected velocity. □

4. Simulation and Results

The effectiveness of the control protocol designed is verified by simulation. This section verifies the improved control protocol of the existing constraints, indicating the effectiveness of the strategies proposed. We assume that the formation has a non-symmetrical communication delay and has the jointly connected topologies in the example below.

We assume that the formation consists of six UAVs. The topology structure and the formation that we expect are shown in Figure 1 and Figure 2.

The communication topology is switched in the order of $\left(G_1,G_2,G_3,G_1\right)$, and the weight of each connected edge is 1. Assuming that there are three different time delays in the system as ${\tau }_1\left(t\right),{\tau }_2\left(t\right),{\tau }_3\left(t\right)$, for $\forall i\ne j$, then ${\tau }_{ii}\left(t\right)={\tau }_{jj}\left(t\right)={\tau }_1\left(t\right)$; ${\tau }_{12}\left(t\right)={\tau }_{23}\left(t\right)={\tau }_{34}\left(t\right)={\tau }_{45}\left(t\right)={\tau }_{56}\left(t\right)={\tau }_{61}\left(t\right)={\tau }_2\left(t\right)$; ${\tau }_{21}\left(t\right)={\tau }_{32}\left(t\right)={\tau }_{43}\left(t\right)={\tau }_{54}\left(t\right)={\tau }_{65}\left(t\right)={\tau }_{16}\left(t\right)={\tau }_3\left(t\right)$.

The time delays satisfy $0\le {\tau }_1\left(t\right)\le 0.01,0\le {\tau }_2\left(t\right)\le 0.07,0\le {\tau }_3\left(t\right)\le 0.08$. The initial state of the six UAVs and the parameters setting are listed in

Table 1. The initial state of the six UAVs.

Number	1	2	3	4	5	6
$x_i/m$	20	60	10	90	43	60
$y_i/m$	66	56	96	56	86	86
$z_i/m$	50	10	40	330	350	240
$v_i/\left(m.s^{-1}\right)$	15	35	55	75	65	90
${\theta }_i/\left(\circ \right)$	36	−36	45	−45	−20	45
$\dot{z}/(m.s^{-1})$	4	3	2	1	5	3

and

Table 2. Parameter settings.

Parameter	${\mathit{v}}_{\mathbf{min}}/\left(\mathit{m}.{\mathit{s}}^{-1}\right)$	${\mathit{v}}_{\mathbf{max}}/\left(\mathit{m}.{\mathit{s}}^{-1}\right)$	${\mathit{a}}_{\mathbf{min}}/\mathit{g}$	${\mathit{a}}_{\mathbf{max}}/\mathit{g}$	${\dot{\mathit{z}}}_{\mathbf{min}}/(\mathit{m}.{\mathit{s}}^{-1})$	${\dot{\mathit{z}}}_{\mathbf{max}}/(\mathit{m}.{\mathit{s}}^{-1})$
Value	10	600	−5	5	−30	30
Parameter	${\ddot{z}}_{\min }/\left(m.s^{-1}\right)$	${\ddot{z}}_{\max }/\left(m.s^{-1}\right)$	${\omega }_{\min }/\left(rad.s^{-1}\right)$	${\omega }_{\max }/\left(rad.s^{-1}\right)$	${\zeta }^{\ast }(m/s)$	$z^{\ast }\left(m\right)$
Value	−5	5	$-\pi /2$	$\pi /2$	50	300
Parameter	$k_1$	$k_2$	$k_3$	${\tau }_v$	${\tau }_z$	${\tau }_{\dot{z}}$
Value	0.6	1.1	0.66	10	0.3	0.3

.

Under the improved control protocol, the position curves, speed curves, course angle curves, and expected formation of the six UAVs are shown in Figure 3.

The figures show that, under the improved formation control protocol, the six drones can achieve the expected formation with the expected speed under the complex conditions of communication constraints and dynamic constraints; the composite error of the formation is 0, as shown in Figure 3f. This indicates that the formation control protocol is effective for UAV formation in the conditions of non-symmetrical communication delay and topology switching.

When drones form a stable formation, assuming that the formation needs to be changed during flight, the control protocol is still valid. The simulation results are shown in Figure 4.

The figures show that, under the formation control protocol, the six drones can achieve the expected formation with the expected speed. When the formation needs to be changed, under the control protocol, the new formation is formed. The drones can fly with the new expected speed and the designed formation control protocol is still valid. The results indicate that the improved control protocol is widely used.

5. Conclusions

This article studies the problem of formation control based on the consistency theory. This article focuses on the research of drone formation, thus ignores the gesture control of the drone. The three-degrees-of-freedom kinematics equation of the UAV is given by using the autopilot model of longitudinal and lateral decoupling. Considering the communication interference and network congestion, this paper designs the control protocol by studying the formation model with non-symmetrical communication delay and switching topology. Acceleration, velocity, and angular velocity constraints in all directions are defined according to the requirements of flight performance and maneuverability. Both communication and mobility constraints are considered in this paper. The improved control protocol is adjusted according to the constraints. The results show that the improved control protocol is effective and can quickly converge the UAV formation state to the specified value and can maintain the specified formation with communication delay and switching topology.