UAV Formation Control under Communication Constraints Based on Distributed Model Predictive Control

Abstract

A distributed model predictive control method is used to transform the formation and maintenance problem of fixed-wing UAV formation during flight into an online rolling optimization problem to solve in this paper. Firstly, the state estimation model of the neighborhood UAV is established according to the relative information of the UAV. Secondly, the error state model in the three-dimensional inertial coordinate frame of the UAV is established without considering the time delay, sensor error, and external interference. Thirdly, a cost function is designed by introducing the error state of the UAV in the neighborhood. Finally, four UAVs are used to verify that under the action of the controller, the UAVs can quickly form and maintain the desired formation while tracking the reference line.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In recent years, with the improvement of UAVs performance and control technology, UAVs have been widely used in both civilian [1–3] and military fields [4, 5]. The coordinated control of UAV formations to perform reconnaissance and surveillance, search and rescue, and the realization of communication relays has become a research hotspot and a difficult problem.

Although there have been many successful studies in this field, in leader-follower formation flying, most of the formation problems are solved on the premise of the existence of active communication links [6–9]. The main problem of the above studies is that the task will fail when the communication link fails. Therefore, in [10–12], it is considered to reduce unnecessary communication or realize formation control in the case of communication link failure. However, in previous studies, vision-based formation control focused on the control of mobile robots [13–16]. In [17], they employ a line-of-sight principle to assign the desired speed and develop an extended state observer based on nonperiodic sampling to provide accurate disturbance estimation with guaranteed convergence. Then, they proposed a robust antidisturbance dynamics control protocol and verified the effectiveness of the algorithm by simulation. But With the development of UAV technology, vision-based formation flight research has become an active research topic. In [18], the formation control law is proposed by using the UAV line-of-sight angle information in the neighborhood, and the formation control of the aircraft is realized by controlling the yaw angle and speed. In [19], Wilson et al. uses a quaternion-based unscented Kalman filter to fuse UAV sensors with relative visual observers to improve the accuracy of relative state estimation. In [20], an improved transfer network is used to reduce the misclassification of images and false alarms of noise during the deep learning process which is able to control the formation of two UAVs in the case of communication refusal. In [21] an appointed-time funnel control consisting of a modified funnel variable is proposed to realize the control of quadrotors suffering from environmental disturbances and parametric uncertainties.

Since in the actual application process, it is not only necessary for a single UAV to enter the formation, but also to realize the coordinated control of multiple UAVs. Therefore, an important issue in practical design is the cooperative control of multiple UAVs. In [22], in order to solve the problem of formation control of multiple UAVs, Dehghani et al. take fixed-wing UAVs as the research object, and give a leader-follower formation kinematics formula based on seeker measurement, and adopted robust controllers and disturbance compensators are used to reduce errors in speed measurements and the effects of bounded disturbances. The effectiveness of the controller when the three UAVs fly in a straight line is verified by simulation. In [23], a double-layer formation control is proposed to solve the problem of UAV formation with relative information. Reference [22] uses sliding mode control to realize the control of UAV formation based on relative information. Reference [24] builds on paper [22] to achieve UAV formation control when one UAV cannot directly measure the relative information of the virtual leader. Also to achieve state estimation, [25] embeds a sampler based on switching thresholds to implement a state observer in the case of intermittent measurements. The paper also proposed a tracking differentiator (TD) -based prescribed performance control(TDPPC) in order to enforce system profiles evolve within the redefined performance boundaries with reduced chattering. Finally, a sigmoid function-based TD is utilized to overcome the complexity explosion.

However, when the number of UAVs increases in the above centralized control method, the amount of calculation increases and the real-time performance is poor. The distributed control algorithm has been widely concerned by researchers because it can effectively reduce the overall calculation amount and improve the autonomy of a single UAV [26, 27]. Among them, distributed model predictive control (DMPC) is more suitable for UAV formation control because of its high structural fault tolerance and strong robustness [28]. For example, in [29], a rendezvous algorithm is proposed for heterogeneous unmanned clusters based on DMPC under limited communication. Reference [30] is focused on problem of time extension and packet loss of the networked control system of UAV with DMPC and event-triggered. Moreover, DMPC is used to avoid collision [31, 32]. In [33], state keeping, obstacles avoidance and precision formation control are performed with the use of DMPC based on introducing binary variables. However, a fundamental distinction of the above methods, this paper adopts different nonlinear error models and control methods to realize formation control which can reflect the relationship between the input and the error more intuitively.

The main contribution of this paper compared with existing research in the literature devoted to three dimension UAV formation control can be expressed in two points of view:

(1) Compared with existing studies [17, 21, 22], those in two-dimensional or simple three-dimensional planes, this paper will realize the control of formation in three-dimensional complex trajectories. Unlike [22], which can only achieve formation control of UAVs with direct access to relative information, and [24], which can achieve formation control of four unmanned vehicles with one UAV cannot directly measure the relative information of the virtual leader, the algorithm in this paper can not only achieve formation control with only neighboring UAV information, but also achieve formation control in complex paths of multiple UAVs.

(2) A new cost function is designed by introducing the UAV error state model in the formation neighborhood, and the stability of the formation system is proved by using the Lyapunov function. And different from [21, 24] which requires the use of large value inputs to achieve the stability of the estimated state, this paper adopts the method of DMPC, which considers the UAV state as well as the input constraints to avoid generating large loads on the UAV and can effectively ensure the robustness of the system.

The rest of the paper is organized as follows. Section 2 establishes a stand-alone motion model of UAV, and the neighborhood UAV state estimation model. In Section 3, we derive a three-dimensional error state model through the stand-alone motion model and design the DMPC controller and give the algorithm process. In Section 4, we present the analysis process of the sufficient conditions for the convergence of the single-machine error model and the stability of the formation system. Section 5 conducts simulation and numerical analysis of three examples through MATLAB. Finally, in Section 6, we draw conclusions and give suggestions for future work.

2. Problem Modeling

In order to introduce the UAV formation movement and modeling process, we introduce three coordinate frames $I$, $V_i$, and $L_i$ as the inertial coordinate frame, the follower velocity frame of UAV $i$, and the line of sight frame. The origin of the follower velocity frame and the line of sight frame is the mass point of UAV. The relative relationship between the coordinate frames is shown in Figure 1.

[figure(s) omitted; refer to PDF]

2.1. Single UAV Motion Model

Assuming that there are $N_u$ fixed-wing UAVs in the formation, the dynamics between UAVs are decoupled, the influence of the roll of the UAVs is not considered, and the influence of the external wind effect is ignored. By considering ${\mathbf{v}}_i^{V_i}={\begin{array}{ccc}{v}_i& 0& 0\end{array}}^T$ as the follower velocity vector with respect to the frame $V_i$, we have$$\begin{array}{}\text{(1)}& {\mathbf{v}}_i^I={\dot{p}}_i={\mathbf{C}}_{V_i}^I{\mathbf{v}}_{\text{i}}^{V_i}=\begin{array}{c}{v}_i\text{cos}{\gamma }_i\text{cos}{\chi }_i\\ v_i\text{cos}{\gamma }_i\text{sin}{\chi }_i\\ -v_i\text{sin}{\gamma }_i\end{array},\end{array}$$where ${\mathbf{C}}_{V_i}^I=\begin{array}{ccc}\text{cos}{\gamma }_i\text{cos}{\chi }_i& -\text{sin}{\chi }_i& \text{sin}{\gamma }_i\text{cos}{\chi }_i\\ \text{cos}{\gamma }_i\text{sin}{\chi }_i& \text{cos}{\chi }_i& \text{sin}{\gamma }_i\text{sin}{\chi }_i\\ -\text{sin}{\gamma }_i& 0& \text{cos}{\gamma }_i\end{array}$ is the rotation matrix of $I$ with respect to $V_i$ , $v_i$ is the speed of the $i$ th follower UAV, and ${\chi }_i$ is the yaw angles of follower UAV, and ${\gamma }_i$ is the pitch angle of follower UAV. ${\mathbf{p}}_i={\begin{array}{ccc}{x}_i& y_i& z_i\end{array}}^T$ is the position vector of the UAV in the inertial coordinate frame.

By denoting ${\mathbf{D}}_I{\mathbf{v}}_i$ as the velocity differential in inertial frame, ${\omega }_i^{V_i}={\begin{array}{ccc}{\omega }_{xi}& {\omega }_{yi}& {\omega }_{zi}\end{array}}^T$ as the angular velocity vector in the follower velocity frame $V_i$, ${\mathbf{a}}_i^{V_i}={\begin{array}{ccc}{a}_{xi}^{V_i}& a_{yi}^{V_i}& a_{zi}^{V_i}\end{array}}^T$ as the acceleration vector in the follower velocity frame, we have the following Coriolis equation:$$\begin{array}{}\text{(2)}& {\mathbf{D}}_I{\mathbf{v}}_i=\begin{array}{c}{a}_{xi}^{V_i}\\ a_{yi}^{V_i}\\ a_{zi}^{V_i}\end{array}={\mathbf{D}}_{V_i}{\mathbf{v}}_i+{\mathbf{\omega }}_i^{V_i}\times {\mathbf{v}}_i=\begin{array}{c}{\dot{v}}_i\\ 0\\ 0\end{array}+\begin{array}{c}{\omega }_{xi}\\ {\omega }_{yi}\\ {\omega }_{zi}\end{array}\times \begin{array}{c}{v}_i\\ 0\\ 0\end{array}=\begin{array}{c}{\dot{v}}_i\\ {\omega }_{zi}{v}_i\\ -{\omega }_{yi}{v}_i\end{array}.\end{array}$$

By denoting ${\mathbf{x}}_i={\begin{array}{ccc}\begin{array}{cc}{x}_i& y_i\end{array}& z_i& \begin{array}{cc}\begin{array}{cc}{v}_i& {\chi }_i\end{array}& {\gamma }_i\end{array}\end{array}}^T\in {ℝ}^6$ as the states of followers, and by using equations (1) and (2), the motion equation of follower can be express as follows:$$\begin{array}{}\text{(3)}& {\dot{x}}_i=\begin{array}{c}\begin{array}{c}{\dot{x}}_i\\ {\dot{y}}_i\end{array}\\ \begin{array}{c}{\dot{z}}_i\\ {\dot{v}}_i\end{array}\\ \begin{array}{c}{\dot{\chi }}_i\\ {\dot{\gamma }}_i\end{array}\end{array}=\begin{array}{c}{v}_i\text{cos}{\gamma }_i\text{cos}{\chi }_i\\ v_i\text{cos}{\gamma }_i\text{sin}{\chi }_i\\ -v_i\text{sin}{\gamma }_i\\ a_{xi}^{V_i}\\ \frac{a_{zi}^{V_i}}{v_i}\\ -\frac{a_{yi}^{V_i}}{v_i}\end{array}.\end{array}$$

By considering ${\mathbf{u}}_i={\begin{array}{ccc}{a}_{xi}^{V_i}& a_{yi}^{V_i}& a_{zi}^{V_i}\end{array}}^T\in {ℝ}^3$ as the input of the system, we have the expression of the UAV state space as follows:$$\begin{array}{}\text{(4)}& {\dot{x}}_{i}t=f{\mathbf{x}}_{i}t,{\mathbf{u}}_{i}t,\end{array}$$where $f·$ is the corresponding nonlinear function.

2.2. Neighborhood UAV State Estimation Model

By denoting $l_{ij}$ as the relative distance between UAV $i$ and $j$, ${\chi }_{ij}$ and ${\gamma }_{ij}$ are the relative pitch and yaw angles; we have the relative kinematics of the leader-follower, as shown in Figure 2.

[figure(s) omitted; refer to PDF]

By assuming ${\mathbf{l}}_{ij}={\begin{array}{ccc}{l}_{ij}& 0& 0\end{array}}^T$ as the vector of position of UAV $j$ with respect to line of sight frame $L_i$, we have the position vector of UAV $j$ in the inertial frame$$\begin{array}{}\text{(5)}& {\mathbf{p}}_j={\mathbf{p}}_i+{\mathbf{C}}_{L_i}^{V_{\text{i}}}·{\mathbf{l}}_{ij},\end{array}$$where ${\mathbf{C}}_{L_i}^{V_i}=\begin{array}{ccc}\text{cos}{\gamma }_{ij}\text{cos}{\chi }_{ij}& \text{cos}{\gamma }_{ij}\text{cos}{\chi }_{ij}& -\text{sin}{\gamma }_{ij}\\ -\text{sin}{\chi }_{ij}& \text{cos}{\chi }_{ij}& 0\\ \text{sin}{\gamma }_{ij}\text{cos}{\chi }_{ij}& \text{sin}{\gamma }_{ij}\text{cos}{\chi }_{ij}& \text{cos}{\gamma }_{ij}\end{array}$. By introducing ${\omega }_{Ij}^{L_i}={\begin{array}{ccc}{\omega }_{xij}& {\omega }_{yij}& {\omega }_{zij}\end{array}}^T$ as the relative angular velocity in the frame $L_i$, we have$$\begin{array}{}\text{(6)}& {\mathbf{v}}_{ij}^{L_i}={\mathbf{D}}_I{\mathbf{l}}_{ij}={\mathbf{D}}_{L_i}{\mathbf{l}}_{ij}+{\mathbf{\omega }}_{Ij}^{L_i}\times {\mathbf{l}}_{ij}=\begin{array}{c}{\dot{l}}_{ij}\\ l_{ij}{\omega }_{zij}\\ -l_{ij}{\omega }_{yij}\end{array},\end{array}$$where ${\mathbf{v}}_{ij}^{L_i}={\begin{array}{ccc}{v}_{xij}& v_{yij}& v_{zij}\end{array}}^T$. Moreover, the velocity of UAV $j$ can be described as follows:$$\begin{array}{}\text{(7)}& {\mathbf{v}}_j^{V_i}={\mathbf{C}}_{L_i}^{V_i}\cdot {\mathbf{v}}_{ij}^{L_i}+{\mathbf{v}}_i^{V_i}.\end{array}$$

Then, the velocity of UAV $j$ in frame $I$ can be described as follows:$$\begin{array}{}\text{(8)}& {\text{v}}_j^{V_i}={\mathbf{C}}_{V_i}^I\cdot {\text{v}}_j^{V_i},\end{array}$$where $·$ is the Euclidean norm of vector. By taking the time derivative of equation (6), we can get$$\begin{array}{}\text{(9)}& {\mathbf{a}}_{ij}^{L_i}={\mathbf{D}}_I{\mathbf{v}}_{ij}={\mathbf{D}}_{L_i}{\mathbf{v}}_{ij}+{\mathbf{\omega }}_{Ij}^{L_i}\times {\mathbf{v}}_{ij},=\begin{array}{c}{\ddot{l}}_{ij}-l_{ij}{\omega }_{zij}^2+{\omega }_{yij}^2\\ 2{\dot{l}}_{ij}{\omega }_{zij}+l_{ij}{\dot{\omega }}_{zij}+l_{ij}{\omega }_{xij}{\omega }_{yij}\\ -2{\dot{l}}_{ij}{\omega }_{yij}-l_{ij}{\dot{\omega }}_{yij}+l_{ij}{\omega }_{xij}{\omega }_{zij}\end{array},\end{array}$$where ${\mathbf{a}}_{ij}^{L_i}={\begin{array}{ccc}{a}_{xij}& a_{yij}& a_{zij}\end{array}}^T$. Moreover, the angular velocity of $L_i$ can be given as follows:$$\begin{array}{}\text{(10)}& {\mathbf{\omega }}_{Ij}^{L_i}={\mathbf{\omega }}_{Ii}^{L_i}+{\mathbf{\omega }}_{ij}^{L_i}={\mathbf{C}}_{V_i}^{L_i}{\mathbf{\omega }}_i^{V_i}+{\mathbf{\omega }}_i^{L_i}=\begin{array}{c}{\omega }_{xi}\text{cos}{\gamma }_{ij}\text{cos}{\chi }_{ij}+{\omega }_{yi}\text{cos}{\gamma }_{ij}\text{sin}{\chi }_{ij}-{\omega }_{zi}\text{sin}{\gamma }_{ij}\\ -{\omega }_{xi}\sin {\gamma }_{ij}+{\omega }_{yi}\text{cos}{\chi }_{ij}\\ {\omega }_{xi}\text{cos}{\gamma }_{ij}\text{sin}{\chi }_{ij}+{\omega }_{yi}\text{sin}{\gamma }_{ij}\text{sin}{\chi }_{ij}+{\omega }_{zi}\text{cos}{\gamma }_{ij}\end{array}+\begin{array}{c}-{\dot{\chi }}_{ij}\text{sin}{\gamma }_{ij}\\ {\gamma }_{ij}\\ {\dot{\chi }}_{ij}\text{cos}{\gamma }_{ij}\end{array}.\end{array}$$

By assuming ${\omega }_{xi}=0$ which means that the UAV does not roll to simplify the model, we have$$\begin{array}{}\text{(11)}& {\mathbf{\omega }}_{Ij}^{L_i}=\begin{array}{c}{\omega }_{yi}\text{cos}{\gamma }_{ij}\text{sin}{\chi }_{ij}-{\omega }_{zi}\text{sin}{\gamma }_{ij}-{\dot{\chi }}_{ij}\text{sin}{\gamma }_{ij}\\ {\omega }_{yi}\text{cos}{\chi }_{ij}+{\gamma }_{ij}\\ {\omega }_{yi}\text{sin}{\gamma }_{ij}\text{sin}{\chi }_{ij}+{\omega }_{zi}\text{cos}{\gamma }_{ij}+{\dot{\chi }}_{ij}\text{cos}{\gamma }_{ij}\end{array}.\end{array}$$

By considering equations (9) and (11), we can have the acceleration of UAV $j$ in the frame $V_i$ as follows:$$\begin{array}{}\text{(12)}& {\stackrel{\wedge }{u}}_j^{V_i}={\mathbf{a}}_i^{V_i}+{\mathbf{C}}_{L_i}^{V_i}\cdot {\stackrel{\wedge }{a}}_{ij}^{L_i},\end{array}$$where ${\stackrel{\wedge }{u}}_j^{V_i}={\begin{array}{ccc}{a}_{xj}^{V_i}& a_{yj}^{V_i}& a_{zj}^{V_i}\end{array}}^T\in {ℝ}^3$. Then, the state estimator of the UAV $j$ can be described as follows:$$\begin{array}{}\text{(13)}& {\stackrel{\wedge }{\dot{x}}}_j=\begin{array}{c}{\mathbf{C}}_{V_i}^I\cdot {\mathbf{v}}_j^{V_i}\\ {\mathbf{C}}_{v_j}\cdot {\stackrel{\wedge }{u}}_j^{V_i}\end{array},\end{array}$$where ${\mathbf{C}}_{v_j}=\begin{array}{ccc}1& 0& 0\\ 0& 0& -1/{\stackrel{\wedge }{v}}_j\\ 0& 1/{\stackrel{\wedge }{v}}_j& 0\end{array}$. In the following DMPC algorithm design, the state of UAV $j$ in equation (10) will be used as the reference state of the UAV that cannot directly measure the relative state of the leader.

3. Distributed Model Predictive Controller Design

In view of the problem that there is no accurate modeling of the three-dimensional space UAV formation error model in the existing literature, this section first establishes the error model based on the state estimation model in the previous section. Then, a cost function based on the error model is given.

Without loss of generality, we give two assumptions of the algorithm as follows:

Assumption 1.

Only the leader can obtain the reference line. Follower UAVs can identify other UAV numbers and obtain relative information of UAVs in the neighborhood without delay according to the information flow transmission topology.

Assumption 2.

The algorithm does not consider the influence of external wind and sensor noise on the UAV.

3.1. State Error Model

Different from [22] and [24], in order to better realize formation control, this paper firstly establishes the UAV error state model in three-dimensional. By assuming ${\mathbf{x}}_{\text{ref}}={\begin{array}{ccc}\begin{array}{cc}{x}_{\text{ref}}& y_{\text{ref}}\end{array}& z_{\text{ref}}& \begin{array}{cc}\begin{array}{cc}{v}_{\text{ref}}& {\chi }_{\text{ref}}\end{array}& {\gamma }_{\text{ref}}\end{array}\end{array}}^T\in {ℝ}^6$ as the reference state of UAV $i$, the error between UAV $i$ and its reference state in the frame $V_i$ can be described as follows:$$\begin{array}{}\text{(14)}& {\mathbf{e}}_{ir}=\begin{array}{cc}{\mathbf{C}}_I^{V_i}& 0\\ 0& {\mathbf{I}}_3\end{array}·{\mathbf{x}}_i-{\mathbf{x}}_{\text{ref}}.\end{array}$$

By taking the time derivative of above equation, one can get$$\begin{array}{c}\text{(15)}& {\dot{e}}_{ir}={\omega }_{yi}\cdot \begin{array}{cc}{\dot{C}}_{{\gamma }_i}& 0\\ 0& {\mathbf{0}}_3\end{array}\cdot \begin{array}{cc}{\mathbf{C}}_{{\chi }_i}& 0\\ 0& {\mathbf{I}}_3\end{array}\cdot {\mathbf{x}}_i-{\mathbf{x}}_{\text{ref}}+{\omega }_{zi}\cdot \begin{array}{cc}{\mathbf{C}}_{{\gamma }_i}& 0\\ 0& {\mathbf{0}}_3\end{array}\cdot \begin{array}{cc}{\dot{C}}_{{\chi }_i}& 0\\ 0& {\mathbf{I}}_3\end{array}\cdot {\mathbf{x}}_i-{\mathbf{x}}_{\text{ref}}\\ +\begin{array}{cc}{\mathbf{C}}_{{\gamma }_i}& 0\\ 0& {\mathbf{I}}_3\end{array}\cdot \begin{array}{cc}{\mathbf{C}}_{{\chi }_i}& 0\\ 0& {\mathbf{I}}_3\end{array}\cdot {\dot{x}}_i-{\dot{x}}_{\text{ref}},\end{array}$$where ${\mathbf{C}}_{{\chi }_i}=\begin{array}{ccc}\text{cos}{\chi }_i& \text{sin}{\chi }_i& 0\\ -\text{sin}{\chi }_i& \text{cos}{\chi }_i& 0\\ 0& 0& 1\end{array}$, ${\mathbf{C}}_{{\gamma }_i}=\begin{array}{ccc}\cos {\gamma }_i& 0& -\sin {\gamma }_i\\ 0& 1& 0\\ \sin {\gamma }_i& 0& \cos {\gamma }_i\end{array}$. By using equations (3) and (11), the error state can be expressed as follows:$$\begin{array}{c}\text{(16)}& {\dot{e}}_{ir}=\begin{array}{cccccc}-\begin{array}{c}{\omega }_{yi}\text{sin}{\gamma }_i\text{cos}{\chi }_i+\\ {\omega }_{zi}\text{cos}{\gamma }_i\text{sin}{\chi }_i\end{array}& \begin{array}{c}-{\omega }_{yi}\text{sin}{\gamma }_i\text{sin}{\chi }_i+\\ {\omega }_{zi}\text{cos}{\gamma }_i\text{cos}{\chi }_i\end{array}& -{\omega }_{yi}\text{cos}{\gamma }_i& 0& 0& 0\\ -{\omega }_{zi}\text{cos}{\chi }_i& -{\omega }_{zi}\text{sin}{\chi }_i& 0& 0& 0& 0\\ \begin{array}{c}{\omega }_{yi}\text{cos}{\gamma }_i\text{cos}{\chi }_i-\\ {\omega }_{zi}\text{sin}{\gamma }_i\text{sin}{\chi }_i\end{array}& \begin{array}{c}{\omega }_{yi}\text{cos}{\gamma }_i\text{sin}{\chi }_i+\\ {\omega }_{zi}\text{sin}{\gamma }_i\text{cos}{\chi }_i\end{array}& -{\omega }_{yi}\text{sin}{\gamma }_i& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\begin{array}{c}\begin{array}{c}{e}_x\\ e_y\end{array}\\ \begin{array}{c}{e}_z\\ e_v\end{array}\\ \begin{array}{c}{e}_{\chi }\\ e_{\gamma }\end{array}\end{array}+\begin{array}{c}{v}_i-v_{\text{ref}}\text{cos}{\gamma }_{\text{ref}}\text{cos}{\gamma }_i\text{cos}{e}_{\chi }-v_{\text{ref}}\text{sin}{\gamma }_{\text{ref}}\text{sin}{\gamma }_i\\ v_{\text{ref}}\text{cos}{\gamma }_{\text{ref}}\text{sin}{e}_{\chi }\\ -v_{\text{ref}}\text{cos}{\gamma }_{\text{ref}}\text{sin}{\gamma }_i\text{cos}{e}_{\chi }+v_{\text{ref}}\text{sin}{\gamma }_{\text{ref}}\text{cos}{\gamma }_i\\ a_{xi}-a_{\text{xref}}\\ \frac{-a_{zi}}{v_i+a_{\text{zref}}/v_{\text{ref}}}\\ \frac{a_{yi}}{v_i-a_{\text{yref}}/v_{\text{ref}}}\end{array}.\end{array}$$

To simplify the calculation, define a virtual input:$$\begin{array}{}\text{(17)}& {\tilde{u}}_{i}k=\begin{array}{c}{\tilde{u}}_{ix}\\ {\tilde{u}}_{iy}\\ {\tilde{u}}_{iz}\end{array}=\begin{array}{c}{a}_{xi}-a_{\text{xref}}\\ {\omega }_{zi}-{\omega }_{\text{zref}}\\ {\omega }_{yi}-{\omega }_{\text{yref}}\end{array}=\begin{array}{c}{a}_{xi}-a_{\text{xref}}\\ \frac{-a_{zi}}{v_i+a_{\text{zref}}/v_{\text{ref}}}\\ \frac{-a_{yi}}{v_i-a_{\text{yref}}/v_{\text{ref}}}{a}_{yi}\end{array}.\end{array}$$

In actual operation, the input can be calculated according to equation (17). By linearizing equation (16) at ${\mathbf{e}}_{ir}0=0,{\dot{e}}_{ir}0=0$, we have$$\begin{array}{}\text{(18)}& {\dot{e}}_{ir}k=\mathbf{A}{\mathbf{e}}_{ir}k+\mathbf{B}{\tilde{u}}_{i}k,\end{array}$$where$$\begin{array}{c}\text{(19)}& \mathbf{A}=\begin{array}{cccccc}-\begin{array}{c}{\omega }_{\text{yref}}\text{sin}{\gamma }_{\text{ref}}\text{cos}{\chi }_{\text{ref}}+\\ {\omega }_{\text{zref}}\text{cos}{\gamma }_{\text{ref}}\text{sin}{\chi }_{\text{ref}}\end{array}& \begin{array}{c}-{\omega }_{\text{yref}}\text{sin}{\gamma }_{\text{ref}}\text{sin}{\chi }_{\text{ref}}+\\ {\omega }_{\text{zref}}\text{cos}{\gamma }_{\text{ref}}\text{cos}{\chi }_{\text{ref}}\end{array}& -{\omega }_{\text{yref}}\text{cos}{\gamma }_{\text{ref}}& 1& 0& 0\\ -{\omega }_{\text{zref}}\text{cos}{\chi }_{\text{ref}}& -{\omega }_{\text{zref}}\text{sin}{\chi }_{\text{ref}}& 0& 0& v_{\text{ref}}\text{cos}{\gamma }_{\text{ref}}& 0\\ \begin{array}{c}{\omega }_{\text{yref}}\text{cos}{\gamma }_{\text{ref}}\text{cos}{\chi }_{\text{ref}}-\\ {\omega }_{\text{zref}}\text{sin}{\gamma }_{\text{ref}}\text{sin}{\chi }_{\text{ref}}\end{array}& \begin{array}{c}{\omega }_{\text{yref}}\text{cos}{\gamma }_{\text{ref}}\text{sin}{\chi }_{\text{ref}}+\\ {\omega }_{\text{zref}}\text{sin}{\gamma }_{\text{ref}}\text{cos}{\chi }_{\text{ref}}\end{array}& -{\omega }_{\text{yref}}\text{sin}{\gamma }_{\text{ref}}& 0& 0& -v_{\text{ref}}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array},\mathbf{B}=\begin{array}{cc}{\mathbf{0}}_3& {\mathbf{0}}_3\\ {\mathbf{0}}_3& {\mathbf{I}}_3\end{array}.\end{array}$$

3.2. Formation Control Design

The DMPC-based control design can be described as, at time $k$, the UAV $i$ solves the optimal control problem with the predicted number of steps starting from the initial error $e_{ir}k$. In the UAV formation, we design different cost function for leader and followers, which can be expressed as follows:

3.2.1. Cost Function of Leader

$$\begin{array}{c}\text{(20)}& \underset{u_i^t}{\text{min}}{J}_{l}k=\underset{u_i^t}{\text{min}}{\displaystyle \sum _{l_1=0}^{N-1}{{\mathbf{e}}_{ir}k+l|k}_{{\mathbf{P}}_i}+{{\mathbf{u}}_{i}k+l|k-{\mathbf{u}}_{ref}k+l|k}_{{\mathbf{R}}_i}}+{\displaystyle \sum _{l_2=0}^{N-2}{\Delta {\mathbf{u}}_{i}k+l|k}_{{\mathbf{Q}}_i}}+{{\mathbf{e}}_{ir}k+N|k}_{{\mathbf{G}}_i}\text{s.t}.\\ {\dot{x}}_{i}k=f{\mathbf{x}}_{i}k,{\mathbf{u}}_{i}k,k=k_0,k_0+1,\dots ,k_0+N-1,\\ {\dot{e}}_{ir}t=\text{g}{\mathbf{e}}_{ir}t,{\tilde{u}}_{i}t,\\ {\mathbf{x}}_{i}k\in {\mathbf{X}}_i,\\ {\mathbf{u}}_{i}k\in {\mathbf{U}}_i.\end{array}$$

3.2.2. Cost Function of Followers

$$\begin{array}{c}\text{(21)}& \underset{u_i^t}{\text{min}}{J}_{f{i}_f}k=\underset{u_i^t}{\text{min}}{\displaystyle \sum _{l_1=0}^{N-1}{\displaystyle \sum _{j\in N_{i_f}}{{\mathbf{e}}_{ij}k+l|k}_{{\mathbf{P}}_i}}+{\displaystyle \sum _{j\in N_{i_f}}{{\mathbf{u}}_{i}k+l|k-{\mathbf{u}}_{j}k+l|k}_{{\mathbf{R}}_i}}}+{\displaystyle \sum _{l_2=0}^{N-2}{\Delta {\mathbf{u}}_{i}k+l|k}_{{\mathbf{Q}}_i}}+{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k+N|k}_{{\mathbf{G}}_i}}\text{s.t}.\\ {\dot{x}}_{i}k=f{\mathbf{x}}_{i}k,{\mathbf{u}}_{i}k,k=k_0,k_0+1,\mathrm{...},k_0+N-1\\ {\dot{e}}_{ir}t=\text{g}{\mathbf{e}}_{ir}t,{\tilde{u}}_{i}t\\ {\mathbf{x}}_{i}k\in {\mathbf{X}}_i\\ {\mathbf{u}}_{i}k\in {\mathbf{U}}_i,\end{array}$$where ${\mathbf{x}}_{\mathbf{P}}={{\mathbf{x}}^T\mathbf{P}\mathbf{x}}^{1/2}$ represents the quadratic form of the vector. $N_f$ represents the set of followers. $N_i$ represents the set of UAVs in the neighborhood of the UAV $i$. ${\mathbf{e}}_{ir}k+l|k$ indicates the error state which the UAV $i$ starts from time $k$ and applies the control sequence ${\mathbf{u}}_{i}k+l|k$ in turn to obtain within the number of predicted steps. ${\mathbf{X}}_i$ is the set of state constraints of UAV $i$ which can be expressed as follows: ${\mathbf{X}}_i=x_{i}k|v_{i}k\le v_{\text{max}},{\chi }_{i}k\le {\chi }_{\text{max}},{\gamma }_{i}k\le {\gamma }_{\text{max}}$. ${\mathbf{U}}_i$ is the input constraint set of the system, which can be described as follows:$$\begin{array}{}\text{(22)}& {\mathbf{U}}_i=\begin{array}{c}{a}_{ix}\le a_{\text{xmax}},a_{iy}\le a_{\text{ymax}},a_{iz}\le a_{\text{zmax}},\\ {\mathbf{u}}_{i}k|\hspace{1em}\hspace{1em}\Delta a_{ix}\le \Delta a_{\text{xmax}},\Delta a_{iy}\le \Delta a_{\text{ymax}},\\ \Delta a_{iz}\le \Delta a_{\text{zmax}}\end{array},\end{array}$$where $\begin{array}{ccc}{a}_{\text{xmax}}& a_{\text{ymax}}& a_{\text{zmax}}\end{array}\in {ℝ}^3$ indicates the maximum value of control input in three directions and $\begin{array}{ccc}\mathrm{\Delta }{a}_{ix}& \mathrm{\Delta }{a}_{iy}& \mathrm{\Delta }{a}_{iz}\end{array}\in {ℝ}^3$ indicates the maximum value of control input change in three directions. The weight matrix satisfies ${\mathbf{P}}_i>0,{\mathbf{Q}}_i>0,{\mathbf{R}}_i>0,{\mathbf{G}}_i>0$.

3.3. Algorithm Flow

We use leader-follower strategy to realize formation control and there is a fixed information flow topology in the formation. In order to distinguish the leader from other UAV in the formation, we use UAV1 to represent the leader. The basic principle of the algorithm is shown in Figure 3.

[figure(s) omitted; refer to PDF]

The process of DMPC can be described as follows:

  Step 1: Bind UAV number and priority. Bind the reference to UAV1.

  Step 2: At $t=k$, the UAV in formation are in a constant cruise state. Determine the UAV in the neighborhood and measure the relative state quantities between the UAV. For the follower UAV: take the current state as the initial state of the UAV, and take the measured state quantity of the UAV in the neighborhood as the reference state quantity. Calculate the current time error as the initial state of rolling time domain planning, which can be described as follows ${\mathbf{e}}_{ir}k|k={\mathbf{e}}_{ir}t=k$.

  Step 3: Calculate the optimal input sequence ${\text{U}}_i^{\ast }k$ which can be described as ${\text{U}}_i^{\ast }k={\text{u}}_i^{\ast }k|k,{\text{u}}_i^{\ast }k+1|k,\mathrm{...},{\text{u}}_i^{\ast }k+N-1|k$ according to equations (20) and (21). The first item of the optimal input sequence ${\text{u}}_i^{\ast }k|k$ is used as the input of the UAV at $t=k+1$.

  Step 4: At $t=k+1$, based on the new error state ${\mathbf{e}}_{ir}t=k+1$ of the formation, repeat to Step 2

4. System Stability Analysis

In the process of solving the UAV formation problem, the stability of the formation plays an important role, which not only reflects the performance of the controller, but also is an important guarantee for the UAV formation system to achieve the expected formation. The main idea of stability analysis is to find a suitable Lyapunov function and prove its degressivity.

According to the above ideas, in the process of cluster stability proof, the following two problems need to be considered: one is the stability of the error state of the single-machine system, and the coefficient relationship is given to ensure the stability of the single-machine system; the second is the stability of the formation system. In the following part, we first give and prove sufficient conditions for single-machine error stability. Then the cost function is used as the Lyapunov function of the system to prove the stability of the formation system.

Theorem 1.

There exists $k_1,k_2\in ℝ$, $k_1>v_{ref},k_2>v_{ref}$ such that the virtual input satisfies$$\begin{array}{}\text{(23)}& {\tilde{u}}_{i}k=-{\mathbf{A}}_{12}\cdot {\mathbf{e}}_1-k_2{\mathbf{e}}_2-{\mathbf{A}}_{12}^{-1}{\mathbf{A}}_{11}+k_1\cdot \mathbf{I}{\mathbf{A}}_{11}\cdot {\mathbf{e}}_1+{\mathbf{A}}_{12}\cdot {\mathbf{e}}_2.\end{array}$$

System equation (18) is asymptotically stable, where ${\mathbf{e}}_1={\begin{array}{ccc}{e}_x& e_y& e_z\end{array}}^T$, ${\mathbf{e}}_1={\begin{array}{ccc}{e}_v& e_{\chi }& e_{\gamma }\end{array}}^T$,$$\begin{array}{c}\text{(24)}& {\mathbf{A}}_{11}=\begin{array}{ccc}-\begin{array}{c}{\omega }_{\text{yref}}\text{sin}{\gamma }_{\text{ref}}\text{cos}{\chi }_{\text{ref}}+\\ {\omega }_{\text{zref}}\text{cos}{\gamma }_{\text{ref}}\text{sin}{\chi }_{\text{ref}}\end{array}& \begin{array}{c}-{\omega }_{\text{yref}}\text{sin}{\gamma }_{\text{ref}}\text{sin}{\chi }_{\text{ref}}+\\ {\omega }_{\text{zref}}\text{cos}{\gamma }_{\text{ref}}\text{cos}{\chi }_{\text{ref}}\end{array}& -{\omega }_{\text{yref}}\text{cos}{\gamma }_{\text{ref}}\\ -{\omega }_{\text{zref}}\text{cos}{\chi }_{\text{ref}}& -{\omega }_{\text{zref}}\text{sin}{\chi }_{\text{ref}}& 0\\ \begin{array}{c}{\omega }_{\text{yref}}\text{cos}{\gamma }_{\text{ref}}\text{cos}{\chi }_{\text{ref}}-\\ {\omega }_{\text{zref}}\text{sin}{\gamma }_{\text{ref}}\text{sin}{\chi }_{\text{ref}}\end{array}& \begin{array}{c}{\omega }_{\text{yref}}\text{cos}{\gamma }_{\text{ref}}\text{sin}{\chi }_{\text{ref}}+\\ {\omega }_{\text{zref}}\text{sin}{\gamma }_{\text{ref}}\text{cos}{\chi }_{\text{ref}}\end{array}& -{\omega }_{\text{yref}}\text{sin}{\gamma }_{\text{ref}}\end{array},A_{12}=\begin{array}{ccc}1& 0& 0\\ 0& v_{\text{ref}}\text{cos}{\gamma }_{\text{ref}}& 0\\ 0& 0& -v_{\text{ref}}\end{array}.\end{array}$$

Proof.

From equation (18), we have$$\begin{array}{}\text{(25)}& \begin{array}{c}{\dot{e}}_1\\ {\dot{e}}_2\end{array}=\begin{array}{cc}{\mathbf{A}}_{11}& {\mathbf{A}}_{12}\\ {\mathbf{0}}_3& {\mathbf{0}}_3\end{array}\begin{array}{c}{\mathbf{e}}_1\\ {\mathbf{e}}_2\end{array}+\begin{array}{cc}{\mathbf{0}}_3& {\mathbf{0}}_3\\ {\mathbf{0}}_3& {\mathbf{B}}_{22}\end{array}\begin{array}{c}{\mathbf{0}}_3\\ {\tilde{u}}_i\end{array}.\end{array}$$

We can rewrite the above equation in the following form:$$\begin{array}{}\text{(26)}& \begin{array}{c}{\dot{e}}_1\\ {\dot{e}}_2\end{array}=\begin{array}{c}{\mathbf{A}}_{11}\cdot {\mathbf{e}}_1+{\mathbf{A}}_{12}\cdot {\mathbf{e}}_2\\ {\mathbf{B}}_{22}\cdot {\tilde{u}}_i\end{array},\end{array}$$where ${\mathbf{B}}_{22}=\mathbf{I}$. By assuming ${\mathbf{e}}_{1d}$ as the expected value of ${\mathbf{e}}_1$, we can describe the error ${\delta }_1$ as follows:$$\begin{array}{}\text{(27)}& {\mathbf{\delta }}_1={\mathbf{e}}_{1d}-{\mathbf{e}}_1.\end{array}$$

By taking the time derivative of the above equation, one can get$$\begin{array}{}\text{(28)}& {\dot{\delta }}_1={\dot{e}}_{1d}-{\dot{e}}_1.\end{array}$$

Choose the Lyapunov function as$$\begin{array}{}\text{(29)}& V_1=\frac{1}{2}{\mathbf{\delta }}_1^T\cdot {\mathbf{\delta }}_1.\end{array}$$

By considering equations (26) and (27), we have$$\begin{array}{}\text{(30)}& {\dot{V}}_1={\mathbf{\delta }}_1^T\cdot {\dot{\delta }}_1={\mathbf{\delta }}_1^T\cdot {\dot{e}}_{1d}-{\mathbf{A}}_{11}{\mathbf{e}}_1-{\mathbf{A}}_{12}{\mathbf{e}}_2.\end{array}$$

Since ${\dot{V}}_1={\delta }_1\cdot {\dot{\delta }}_1<0$, we have ${\dot{\delta }}_1={\dot{e}}_{1d}-{\mathbf{A}}_{11}{\mathbf{e}}_1-{\mathbf{A}}_{12}{\mathbf{e}}_2⟶-k_1\cdot {\mathbf{e}}_1$, where $k_1>0$. By assuming ${\mathbf{e}}_{2d}$ as the expected value of ${\mathbf{e}}_2$, we have$$\begin{array}{}\text{(31)}& {\mathbf{e}}_{2d}={\mathbf{A}}_{12}^{-1}\cdot {\dot{e}}_{1d}-{\mathbf{A}}_{11}\cdot {\mathbf{e}}_1+k_1\cdot {\mathbf{\delta }}_1.\end{array}$$

The error ${\delta }_2$ for ${\mathbf{e}}_2$ can be expressed as$$\begin{array}{}\text{(32)}& {\mathbf{\delta }}_2={\mathbf{e}}_{2d}-{\mathbf{e}}_2.\end{array}$$

Substituting equations (31) and (32) into (30), we have$$\begin{array}{c}\text{(33)}& {\dot{V}}_1={\mathbf{\delta }}_1^T\cdot {\dot{e}}_{1d}-{\mathbf{A}}_{11}{\mathbf{e}}_1-{\mathbf{A}}_{12}{\mathbf{e}}_2={\mathbf{\delta }}_1^T\cdot {\dot{e}}_{1d}-{\mathbf{A}}_{11}{\mathbf{e}}_1-{\mathbf{A}}_{12}\cdot {\mathbf{A}}_{12}^{-1}\cdot {\dot{e}}_{1d}-{\mathbf{A}}_{11}\cdot {\mathbf{e}}_1+k_1\cdot {\mathbf{\delta }}_1-{\mathbf{\delta }}_2\\ =-k_1\cdot {\mathbf{\delta }}_1^T\cdot {\mathbf{\delta }}_1+{\mathbf{A}}_{12}\cdot {\mathbf{\delta }}_1^T\cdot {\mathbf{\delta }}_2.\end{array}$$

By substituting equations (26) and (31) into (32) and taking the time derivative of it, we have$$\begin{array}{}\text{(34)}& {\dot{\delta }}_2={\dot{e}}_{2d}-{\dot{e}}_2={\mathbf{A}}_{12}^{-1}{\ddot{e}}_{1d}+k_1{\dot{e}}_{1d}-{\mathbf{A}}_{11}+k_1\cdot \mathbf{I}{\mathbf{A}}_{11}\cdot {\mathbf{e}}_1+{\mathbf{A}}_{12}\cdot {\mathbf{e}}_2-{\mathbf{B}}_{22}\cdot {\tilde{u}}_i.\end{array}$$

Let the Lyapunov function with ${\delta }_2$ be$$\begin{array}{}\text{(35)}& V_2=V_1+\frac{1}{2}{\mathbf{\delta }}_2^T\cdot {\mathbf{\delta }}_2.\end{array}$$

By taking the time derivative of equation (35) and substituting equation (33) into it, we have$$\begin{array}{}\text{(36)}& {\dot{V}}_2=-k_1\cdot {\mathbf{\delta }}_1^T\cdot {\mathbf{\delta }}_1+{\mathbf{\delta }}_1^T{\mathbf{A}}_{12}\cdot {\mathbf{\delta }}_1+{\dot{\delta }}_2,\end{array}$$where $k_2>0$, By substituting equations (26), (32), and (34) into the above equation, one can get the following equation:$$\begin{array}{}\text{(37)}& {\mathbf{A}}_{12}\cdot {\mathbf{\delta }}_1+{\mathbf{A}}_{12}^{-1}{\ddot{e}}_{1d}+k_1{\dot{e}}_{1d}-{\mathbf{A}}_{11}+k_1{\mathbf{A}}_{11}\cdot {\mathbf{e}}_1+{\mathbf{A}}_{12}\cdot {\mathbf{e}}_2-{\mathbf{B}}_{22}\cdot {\tilde{u}}_i=-k_2{\mathbf{\delta }}_2.\end{array}$$

By rewriting the above equation, we can get the equation of the input:$$\begin{array}{}\text{(38)}& {\tilde{u}}_i={\mathbf{B}}_{22}^{-1}{\mathbf{A}}_{12}\cdot {\mathbf{\delta }}_1+{\mathbf{A}}_{12}^{-1}{\ddot{e}}_{1d}+k_1{\dot{e}}_{1d}-{\mathbf{A}}_{11}+k_1{\mathbf{A}}_{11}\cdot {\mathbf{e}}_1+{\mathbf{A}}_{12}\cdot {\mathbf{e}}_2+k_2{\mathbf{\delta }}_2.\end{array}$$

By substituting equations (38) into (34), we have$$\begin{array}{c}\text{(39)}& {\dot{\delta }}_2={\dot{e}}_{2d}-{\dot{e}}_2={\mathbf{A}}_{12}^{-1}{\ddot{e}}_{1d}+k_1{\dot{e}}_{1d}-{\mathbf{A}}_{11}+k_1\cdot \mathbf{I}{\mathbf{A}}_{11}\cdot {\mathbf{e}}_1+{\mathbf{A}}_{12}\cdot {\mathbf{e}}_2-{\mathbf{B}}_{22}\cdot {\mathbf{B}}_{22}^{-1}{\mathbf{\delta }}_1+{\mathbf{A}}_{12}^{-1}{\ddot{e}}_{1d}+k_1{\dot{e}}_{1d}-{\mathbf{A}}_{11}+k_1{\mathbf{A}}_{11}\cdot {\mathbf{e}}_1+{\mathbf{A}}_{12}\cdot {\mathbf{e}}_2+k_2{\mathbf{\delta }}_2\\ =-{\mathbf{A}}_{12}\cdot {\mathbf{\delta }}_1-k_2\cdot {\mathbf{\delta }}_2.\end{array}$$

By substituting equations (26) and (31) into (28), one can get$$\begin{array}{}\text{(40)}& \begin{array}{c}{\dot{\delta }}_1\\ {\dot{\delta }}_2\end{array}=\begin{array}{cc}-k_1\cdot \mathbf{I}& {\mathbf{A}}_{12}\\ -{\mathbf{A}}_{12}& -k_2\cdot \mathbf{I}\end{array}\begin{array}{c}{\mathbf{\delta }}_1\\ {\mathbf{\delta }}_2\end{array}.\end{array}$$

Let the eigenvalues of system equation (40) be ${\lambda }_i,i=\mathrm{1,2},\mathrm{...},6$, corresponding to the eigenvalues of the 1st to 6th rows of the eigenmatrix. We have$$\begin{array}{}\text{(41)}& {\displaystyle \prod _{i=1}^6{\lambda }_i}=-k_1\cdot \mathbf{I}\cdot -k_2\cdot \mathbf{I}--{\mathbf{A}}_{12}\cdot {\mathbf{A}}_{12}={k_1{k}_2}^3+{{\mathbf{A}}_{12}}^2{\displaystyle \sum _{i=1}^6{\lambda }_i}=-3\cdot k_1+k_2,\end{array}$$where $\mathbf{A}$ represents the determinant of $\mathbf{A}$. Since $k_1>0$, $k_2>0$, considering equation (41), one can get that the real part of the conjugate eigenvalue of system equation (40) has the same sign. And by considering Gershgorin’s disc theorem, we have$$\begin{array}{}\text{(42)}& \begin{array}{c}{\lambda }_1+k_1\le 1\\ {\lambda }_2+k_1\le v_{\text{ref}}\text{cos}{\gamma }_{\text{ref}}\\ {\lambda }_3+k_1\le v_{\text{ref}}\\ {\lambda }_4+k_2\le 1\\ {\lambda }_5+k_2\le v_{\text{ref}}\text{cos}{\gamma }_{\text{ref}}\\ {\lambda }_6+k_1\le v_{\text{ref}}\end{array}.\end{array}$$

By considering equation (42), when ${\lambda }_i>\max v_{ref},1,i=\mathrm{1,2},\mathrm{...},6$, the eigenvalues of system equation (40) are all in the negative half-plane and are asymptotically stable. By considering ${\mathbf{e}}_{1d}={\dot{e}}_{1d}={\ddot{e}}_{1d}=0$ and ${\mathbf{e}}_{2d}={\dot{e}}_{2d}=0$, and substituting equations (27) and (32) into (38), we have$$\begin{array}{}\text{(43)}& {\tilde{u}}_i=-{\mathbf{A}}_{12}\cdot {\mathbf{e}}_1-k_2{\mathbf{e}}_2-{\mathbf{A}}_{12}^{-1}{\mathbf{A}}_{11}+k_1\cdot \mathbf{I}{\mathbf{A}}_{11}\cdot {\mathbf{e}}_1+{\mathbf{A}}_{12}\cdot {\mathbf{e}}_2.\end{array}$$

Hence, the system is stable under the action of the input equation (43). Theorem 1 gives the sufficient conditions for the stability of the single-UAV error system. Then the cost function is taken as the Lyapunov function to prove the stability of the UAV formation system. According to the one-to-one correspondence between the state space, the error state space, the virtual input and the actual input, the solution of the actual input sequence ${\text{U}}_i^{\ast }k$ can be transformed into the solution of the virtual input sequence ${\tilde{U}}_{i}k$. By rewriting equations (20) and (21) with the virtual input ${\tilde{u}}_{i}k$ in equation (17), we have a new expression for the cost function as follows:$$\begin{array}{c}\text{(44)}& \underset{u_i^t}{\text{min}}{J}_{l}k=\underset{u_i^t}{\text{min}}{\displaystyle \sum _{l_1=0}^{N-1}{{\mathbf{e}}_{ir}k+l|k}_{{\mathbf{P}}_i}+{{\tilde{u}}_{i}k+l|k}_{{\mathbf{R}}_i}}+{\displaystyle \sum _{l_2=0}^{N-2}{\Delta {\tilde{u}}_{i}k+l|k}_{{\mathbf{Q}}_i}}+{{\mathbf{e}}_{ir}k+N|k}_{{\mathbf{G}}_i}\text{s.t}.\\ {\dot{x}}_{i}k=f{\mathbf{x}}_{i}k,{\mathbf{u}}_{i}k,k=k_0,k_0+1,\dots ,k_0+N-1\\ {\dot{e}}_{ir}t=\text{g}{\mathbf{e}}_{ir}t,{\tilde{u}}_{i}t\\ {\mathbf{x}}_{i}k\in {\mathbf{X}}_i\\ {\mathbf{u}}_{i}k\in {\mathbf{U}}_i,\\ \text{(45)}& \underset{u_i^t}{\text{min}}{J}_{f{i}_f}k=\underset{u_i^t}{\text{min}}{\displaystyle \sum _{l_1=0}^{N-1}{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k+l|k}_{{\mathbf{P}}_i}}+{\displaystyle \sum _{j\in N_i}{{\tilde{u}}_{i}k+l|k}_{{\mathbf{R}}_i}}}+{\displaystyle \sum _{l_2=0}^{N-2}{\Delta {\tilde{u}}_{i}k+l|k}_{{\mathbf{Q}}_i}}+{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k+N|k}_{{\mathbf{G}}_i}}\text{s.t}.\\ {\dot{x}}_{i}k=f{\mathbf{x}}_{i}k,{\mathbf{u}}_{i}k,k=k_0,k_0+1,\dots ,k_0+N-1\\ {\dot{e}}_{ir}t=\text{g}{\mathbf{e}}_{ir}t,{\tilde{u}}_{i}t\\ {\mathbf{x}}_{i}k\in {\mathbf{X}}_i\\ {\tilde{u}}_{i}k\in {\tilde{U}}_i.\end{array}$$

At $t=k$, the optimal solution sequence of equation (45) is ${\tilde{U}}_{i}k={\tilde{u}}_i^{\ast }k+l|k,l=0,\mathrm{..},N-1$. The optimal state sequence corresponding to the solution is ${\mathbf{X}}_{i}k={\mathbf{x}}_i^{\ast }k+l|k,l=0,\mathrm{..},N-1$. The error state sequence corresponding to the solution is ${\mathbf{e}}_i^{\ast }k={\mathbf{e}}_i^{\ast }k+l|k,l=0,\mathrm{..},N-1$. Obviously, all components of the optimal state sequence ${\mathbf{X}}_{i}k$ satisfy the constraints.

At $t=k+1$, we construct the solution of equation (45) which can be expressed as follows:$$\begin{array}{}\text{(46)}& {\tilde{u}}_{i}k+1=\begin{array}{c}{\tilde{u}}_i^{\ast }k+l|k,l=\mathrm{0,1,2},\dots ,N-2\hfill \\ \mathbf{C}{\mathbf{e}}_{1}k+N|k+\mathbf{D}{\mathbf{e}}_{2}k+N|k,l=N-1\hfill \end{array},\end{array}$$where $\mathbf{C}=-{\mathbf{A}}_{12}+{\mathbf{A}}_{12}^{-1}{\mathbf{A}}_{11}+k_1\cdot \mathbf{I}{\mathbf{A}}_{11}$, $\mathbf{D}=-k_2\cdot \mathbf{I}+{\mathbf{A}}_{12}^{-1}{\mathbf{A}}_{11}+k_1\cdot \mathbf{I}{\mathbf{A}}_{12}$. By considering equation (45), we can get the state sequence and error state sequence of the system expressed as follows:$$\begin{array}{c}\text{(47)}& {\mathbf{X}}_{i}k+1={\mathbf{x}}_{i}k+1+l|k,l=1,\dots ,N=\begin{array}{c}{\mathbf{x}}_i^{\ast }k+2|k,{\mathbf{x}}_i^{\ast }k+3|k,\dots ,\\ {\mathbf{x}}_i^{\ast }k+N|k,{\mathbf{x}}_{i}k+N|k+1\end{array}{\mathbf{e}}_{i}k+1={\mathbf{e}}_{ir}k+1+l|k,l=0,\dots ,N-1=\begin{array}{c}{\mathbf{e}}_{ir}k+2|k,{\mathbf{e}}_{ir}k+2|k,\dots ,\\ {\mathbf{e}}_{ir}k+N|k,{\mathbf{e}}_{ir}k+N|k+1\end{array}.\end{array}$$

By considering the first $N-1$ of ${\mathbf{X}}_{i}k+1$ are the components of ${\mathbf{X}}_{i}k$, the first $N-1$ of ${\mathbf{e}}_{i}k+1$ are the components of ${\mathbf{e}}_i^{\ast }k$, all of them satisfy the constraints. By substituting ${\mathbf{x}}_i^{\ast }k+N|k={\mathbf{x}}_{ref}^{\ast }k+N|k$ into (10), one can get$$\begin{array}{}\text{(48)}& {\mathbf{e}}_{ir}^{\ast }k+N|k={\mathbf{x}}_i^{\ast }k+N|k-{\mathbf{x}}_{ref}k+N|k={\mathbf{0}}_{6\times 1}.\end{array}$$

By considering equation (40) and the input which can be expressed as $\mathbf{C}{\mathbf{e}}_{ir1}^{\ast }k+N|k+\mathbf{D}{\mathbf{e}}_{ir2}^{\ast }k+N|k$, we have$$\begin{array}{}\text{(49)}& \begin{array}{c}{\mathbf{e}}_{ir1}k+N|k+1-{\mathbf{e}}_{ir1}^{\ast }k+N|k\\ {\mathbf{e}}_{ir2}k+N|k+1-{\mathbf{e}}_{ir2}^{\ast }k+N|k\end{array}=\begin{array}{c}{\dot{e}}_{ir1}k+N|k+1\\ {\dot{e}}_{ir2}k+N|k+1\end{array}=\begin{array}{cc}-k_1\cdot \mathbf{I}& {\mathbf{A}}_{12}\\ -{\mathbf{A}}_{12}& -k_2\cdot \mathbf{I}\end{array}\begin{array}{c}{\mathbf{e}}_{ir1}^{\ast }k+N|k\\ {\mathbf{e}}_{ir2}^{\ast }k+N|k\end{array}={\mathbf{0}}_{6\times 1}.\end{array}$$

It can be rewriting as ${\mathbf{e}}_{ir}k+N|k={\mathbf{x}}_{i}k+N|k+1-{\mathbf{x}}_{ref}k+N|k={\mathbf{0}}_{6\times 1}$. By considering all the above equations, we can get that ${\mathbf{X}}_{i}k+1$ meet terminal constraints and ${\tilde{u}}_{i}k+1$ is a set of feasible solutions to equation (45). Similarly, we can get a set of feasible solutions ${\tilde{u}}_{l}k+1$ to equation (44). At $t=k$, the cost functions of leader and followers can be expressed as follows:$$\begin{array}{c}\text{(50)}& J_{l}k={\displaystyle \sum _{l_1=0}^{N-1}{{\mathbf{e}}_{lr}k+l|k}_{{\mathbf{P}}_i}+{{\tilde{u}}_{l}k+l|k}_{{\mathbf{R}}_i}}+{\displaystyle \sum _{l_2=0}^{N-2}{\Delta {\mathbf{u}}_{l}k+l|k}_{{\mathbf{Q}}_i}}+{{\mathbf{e}}_{lr}k+N|k}_{{\mathbf{G}}_i},\\ J_{f{i}_f}k={\displaystyle \sum _{l_1=0}^{N-1}{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k+l|k}_{{\mathbf{P}}_i}}+{\displaystyle \sum _{j\in N_i}{{\tilde{u}}_{i}k+l|k}_{{\mathbf{R}}_i}}}\\ +{\displaystyle \sum _{l_2=0}^{N-2}{\Delta {\tilde{u}}_{i}k+l|k}_{{\mathbf{Q}}_i}}+{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k+N|k}_{{\mathbf{G}}_i}}.\end{array}$$

At $t=k+1$,$$\begin{array}{c}\text{(51)}& J_{l}k+1={\displaystyle \sum _{l_1=0}^{N-1}{{\mathbf{e}}_{lr}k+l+1|k+1}_{{\mathbf{P}}_i}+{{\tilde{u}}_{l}k+l+1|k+1}_{{\mathbf{R}}_i}}+{\displaystyle \sum _{l_2=0}^{N-2}{\Delta {\mathbf{u}}_{lr}k+l+1|k+1}_{{\mathbf{Q}}_i}}+{{\mathbf{e}}_{lr}k+l+1|k+1}_{{\mathbf{G}}_i}\\ =J_l^{\ast }k-{{\mathbf{e}}_{lr}k|k}_{{\mathbf{P}}_i}-{{\tilde{u}}_{l}k|k}_{{\mathbf{R}}_i}-{\Delta {\mathbf{u}}_{lr}k|k}_{{\mathbf{Q}}_i}\\ +{{\mathbf{e}}_{lr}k+N|k}_{{\mathbf{P}}_i}+{{\tilde{u}}_{l}k+N|k}_{{\mathbf{R}}_i}+{\Delta {\mathbf{u}}_{lr}k+N|k}_{{\mathbf{Q}}_i}\\ +{{\mathbf{e}}_{lr}k+N+1|k}_{{\mathbf{G}}_i},J_{f{i}_f}k+1\\ ={\displaystyle \sum _{l_1=0}^{N-1}{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k+l+1|k+1}_{{\mathbf{P}}_i}}+{\displaystyle \sum _{j\in N_i}{{\tilde{u}}_{i}k+l+1|k+1}_{{\mathbf{R}}_i}}}\\ +{\displaystyle \sum _{l_2=0}^{N-2}{\Delta {\tilde{u}}_{ij}k+l+1|k+1}_{{\mathbf{Q}}_i}}+{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k+l+1|k+1}_{{\mathbf{G}}_i}}\\ =J_{f{i}_f}^{\ast }k+1-{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k|k}_{{\mathbf{P}}_i}}-{\displaystyle \sum _{j\in N_i}{{\tilde{u}}_{i}k|k}_{{\mathbf{R}}_i}}\\ -{\displaystyle \sum _{j\in N_i}{\Delta {\tilde{u}}_{i}k|k}_{{\mathbf{Q}}_i}}+{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k+N|k+1}_{{\mathbf{P}}_i}}+{\displaystyle \sum _{j\in N_i}{{\tilde{u}}_{i}k+N|k}_{{\mathbf{R}}_i}}\\ +{\displaystyle \sum _{j\in N_i}{\Delta {\tilde{u}}_{ij}k+N|k}_{{\mathbf{Q}}_i}}+{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k+N+1|k}_{{\mathbf{G}}_i}}.\end{array}$$

By considering ${{\mathbf{e}}_{ir}k+N|k}_{{\mathbf{P}}_i}={{\mathbf{e}}_{ir}k+N|k}_{{\mathbf{G}}_i}=0$ and ${\mathrm{\Delta }{\tilde{u}}_{ir}k+N|k}_{{\mathbf{Q}}_i}={\mathrm{\Delta }{\tilde{u}}_{lr}k+N|k}_{{\mathbf{Q}}_i}=0$, we can rewriting the above equations as follows:$$\begin{array}{c}\text{(52)}& J_{l}k+1=J_l^{\ast }k-{{\mathbf{e}}_{lr}k|k}_{{\mathbf{P}}_i}-{{\tilde{u}}_{l}k|k}_{{\mathbf{R}}_i}-{\Delta {\mathbf{u}}_{lr}k|k}_{{\mathbf{Q}}_i}\le J_l^{\ast }k,\\ J_{f{i}_f}k+1=J_{f{i}_f}^{\ast }k+1-{\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k|k}_{{\mathbf{P}}_i}}-{\displaystyle \sum _{j\in N_i}{{\tilde{u}}_{i}k|k}_{{\mathbf{R}}_i}}\\ -{\displaystyle \sum _{j\in N_i}{\Delta {\tilde{u}}_{i}k|k}_{{\mathbf{Q}}_i}}\\ \le J_{f{i}_f}^{\ast }k.\end{array}$$

By subtracting the cost function at time $k$ from both the sides at the same time, we can get$$\begin{array}{}\text{(53)}& J_l^{\ast }k+1-J_l^{\ast }k\le 0,J_{f{i}_f}^{\ast }k+1-J_{f{i}_f}^{\ast }k\le 0.\end{array}$$

It can be rewritten as$$\begin{array}{}\text{(54)}& J_l^{\ast }k+1+{\displaystyle \sum J_{f{i}_f}^{\ast }k+1}-J_l^{\ast }k+{\displaystyle \sum J_{f{i}_f}^{\ast }k}\le 0,\end{array}$$where the equal sign holds if and only if ${{\mathbf{e}}_{lr}k|k}_{{\mathbf{P}}_i}={\displaystyle \sum _{j\in N_i}{{\mathbf{e}}_{ij}k|k}_{{\mathbf{P}}_i}}=0$. According to equation (54), the UAV DMPC system designed in this paper is asymptotically stable.

5. Simulation

In this section, four UAVs are used for simulation. The UAV formation is shown in Figure 4, and the arrows indicate the flow of information.

[figure(s) omitted; refer to PDF]

Table 1 shows the initial state of the UAV in the inertial frame and the expected relative position of each UAV and UAV1. Set the maximum value of the input in the three directions of the UAV to 12 m/s and the minimum value to -12 m/s; the maximum value of the speed in the three directions is 20 m/s; the maximum value of the acceleration change is 10 m/s². The number of simulations is 200 times, and the sampling time is 0.1s. In order to verify the effectiveness of the algorithm, the simulation will be carried out in three cases: straight flight, fixed turning radius flight and variable turning radius.

Table 1

Initial condition of UAV and desired relative position of followers.

UAV	Initial position	$v\text{m}/\text{s}$	Desired relative position
Reference line	${\begin{array}{ccc}5& 5& 0\end{array}}^T$		—
1	${\begin{array}{ccc}2& 2& 0\end{array}}^T$	8	${\begin{array}{ccc}0& 0& 0\end{array}}^T$
2	${\begin{array}{ccc}-2& 3& 0\end{array}}^T$	8	${\begin{array}{ccc}-5& -5& 0\end{array}}^T$
3	${\begin{array}{ccc}2& -1& 0\end{array}}^T$	8	${\begin{array}{ccc}5& -5& 0\end{array}}^T$
4	${\begin{array}{ccc}-2& -2& 0\end{array}}^T$	8	${\begin{array}{ccc}0& -10& 0\end{array}}^T$

Case 1.

Straight Flight

Figure 5 shows the result of the formulation of the UAV formation when they are flying in a straight line. Figure 5(a) shows the trajectories of the UAV formation. It can be found that the four UAVs can form and maintain a desired formation shape and track the reference line. Figure 5(b) shows the relative distance between UAVs. In this case the distance converges after about 4.3 s. Figure 5(c) shows the velocity of each UAV in three directions. It can be found that the velocity in each direction does not exceed the maximum value and changes smoothly. The control trajectories of UAV are shown in Figure 5(d). In Figure 5(d), the values of the first moment are not equal. This is because the value is the input calculated after the first sample period. In order to test the effect of the algorithm when the formation performs a turning flight, we also conduct experiments of Case 2.

[figure(s) omitted; refer to PDF]

Case 2.

Fixed Turning Radius Flight

Figure 6 shows the result of the formulation of the UAV formation while they are flying with a fixed turning radius. Four UAVs can form and maintain desired formation when they are flying with a fixed radius (Figure 6(a)). Compared to Case 1, the convergence rate of distance shown in Figure 6(b) is slower. It can be found in Figure 6(c) that the velocity can converge to a fixed value determined from the radius of rotation and does not exceed a maximum value. Compared to Figure 5(d), the input curve converges to a fluctuating curve. This is because the direction of the drone needs to be adjusted in time to track the trajectory. In order to further verify the effectiveness of the algorithm, we conduct simulation experiments in Case 3.

[figure(s) omitted; refer to PDF]

Case 3.

Variable Turning Radius Flight

Figure 7 shows the simulation results of the UAV formation under the condition of varying turning radius. The formation of four UAVs can ensure the tracking of the reference line when the turning radius is not large in the early stage. When the turning radius becomes larger, the formation has an error in the tracking of the reference line because the speed has reached the maximum value (Figure 7(c)). This problem will be solved in the next step by improving the performance of the UAV or the control algorithm. However, it can be found in Figure 7(b) that although the distance between UAVs is slower than when flying in a straight line, the four UAVs can form and maintain the ideal formation error. It can be seen from Figure 7(d) that the four UAVs can achieve convergence under the limit of the maximum input value. In order to further verify the control effect of the algorithm on multisortie UAV clusters, we set up Case 4.

[figure(s) omitted; refer to PDF]

Case 4.

Multi-UAVs. (Table2)

Unlike the previous cases with four UAVs, Case 4 simulates the tracking and formation control of eight UAVs in complex trajectories with changing turn radius, and the simulation results are shown in Figure 8. From Figure 8(a), we can see that the algorithm has a better control effect for large-scale clusters. At the same time, the interaircraft distance of UAVs can converge to the specified value relatively quickly (Figure 8(b)). From Figures 8(c) and 8(d), we can see that the algorithm can guarantee the control of the formation while ensuring that the UAV states and inputs satisfy the performance constraints.

Table 2

The initial condition of multi UAVs and desired relative position of followers.

UAV	Initial position	$v\text{m}/\text{s}$	Desired relative position
Reference line	${\begin{array}{ccc}5& 5& 0\end{array}}^T$		—
1	${\begin{array}{ccc}2& 2& 0\end{array}}^T$	8	${\begin{array}{ccc}0& 0& 0\end{array}}^T$
2	${\begin{array}{ccc}-2& 3& 0\end{array}}^T$	8	${\begin{array}{ccc}-5& -5& 0\end{array}}^T$
3	${\begin{array}{ccc}2& -1& 0\end{array}}^T$	8	${\begin{array}{ccc}0& -5& 0\end{array}}^T$
4	${\begin{array}{ccc}-2& -2& 0\end{array}}^T$	8	${\begin{array}{ccc}5& -5& 0\end{array}}^T$
5	${\begin{array}{ccc}2& 4& 0\end{array}}^T$	8	${\begin{array}{ccc}-5& -10& 0\end{array}}^T$
6	${\begin{array}{ccc}-4& -4& 0\end{array}}^T$	8	${\begin{array}{ccc}0& -10& 0\end{array}}^T$
7	${\begin{array}{ccc}4& -2& 0\end{array}}^T$	8	${\begin{array}{ccc}5& -10& 0\end{array}}^T$
8	${\begin{array}{ccc}-4& -4& 0\end{array}}^T$	8	${\begin{array}{ccc}0& -15& 0\end{array}}^T$

[figure(s) omitted; refer to PDF]

6. Conclusions

The DMPC-based scheme for formation control of multiple UAVs is studied in this paper. The UAV formation error state equation in three-dimensional case is proposed. By introducing the error of UAV formation, two cost functions are established according to the difference between the leader and the follower. The sufficient conditions for the stability of the UAV formation system are given and proved and the algorithm is simulated in three cases. The numerical results show that the UAV formation of four UAVs can form and maintain the desired formation under condition of straight flight and fixed turning radius. However, in the case of variable turning radius flight, when the turning radius is large, there is an error between the UAV and the reference line due to the limitation of the maximum speed of the UAV. Future work will focus on improving the performance of the DMPC algorithm in order to realize the tracking the reference line with large turning radius and resist external disturbances such as wind speed and sensor errors.

Acknowledgments

This research was funded by the Shandong Natural Science Foundation of China (ZR2020QF057).

References

[1] V. P. Bacheti, A. S. Brandao, M. Sarcinelli-Filho, "Path-Following with a UGV-UAV Formation Considering that the UAV Lands on the UGV," .

[2] L. J. R. Z. Yuanchen Zhao, "UAV formation control with obstacle avoidance using improved artificial potential fields," ,DOI: 10.23919/chicc.2017.8028347, .

[3] B. Fan, Y. Li, R. Zhang, Q. Fu, "Review on the technological development and application of UAV systems," Chinese Journal of Electronics, vol. 29 no. 2, pp. 199-207, DOI: 10.1049/cje.2019.12.006, 2020.

[4] V. Roberge, M. Labonte, G. Labonte, "Fast genetic algorithm path planner for fixed-wing military UAV using GPU," IEEE Transactions on Aerospace and Electronic Systems, vol. 54 no. 5, pp. 2105-2117, DOI: 10.1109/taes.2018.2807558, 2018.

[5] S. K. Chaturvedi, R. Sekhar, S. Banerjee, H. Kamal, "Comparative review study of military and civilian unmanned aerial vehicles (UAVs)," INCAS bulletin, vol. 11 no. 3, pp. 183-198, DOI: 10.13111/2066-8201.2019.11.3.16, 2019.

[6] O. Turnbull, J. Lawry, M. Lowenberg, A. Richards, "A cloned linguistic decision tree controller for real-time path planning in hostile environments," Fuzzy Sets and Systems, vol. 293,DOI: 10.1016/j.fss.2015.08.017, 2016.

[7] J. Zhang, J. Yan, P. Zhang, X. Kong, "Collision avoidance in fixed-wing UAV formation flight based on a consensus control algorithm," IEEE Access, vol. 6, pp. 43672-43682, DOI: 10.1109/access.2018.2864169, 2018.

[8] J. Zhang, J. Zhang, P. Zhang, "Multi-UAV formation control based on a novel back-stepping approach," IEEE Transactions on Vehicular Technology, vol. 69 no. 3, pp. 2437-2448, DOI: 10.1109/tvt.2020.2964847, 2020.

[9] C. L. Zhao, S. W. Dai, G. R. Zhao, "UAV Formation Control Based on Distributed Model Predictive," Control Control And Decision Making, 2021.

[10] Y. Kartal, K. Subbarao, N. R. Gans, A. Dogan, F. Lewis, "Distributed backstepping based control of multiple UAV formation flight subject to time delays," IET Control Theory & Applications, vol. 14 no. 12, pp. 1628-1638, DOI: 10.1049/iet-cta.2019.1151, 2020.

[11] Y. Lin, M. Wang, X. Zhou, G. Ding, S. Mao, "Dynamic spectrum interaction of UAV flight formation communication with priority: a deep reinforcement learning approach," IEEE Transactions on Cognitive Communications and Networking, vol. 6 no. 3, pp. 892-903, DOI: 10.1109/tccn.2020.2973376, 2020.

[12] Q. Hu, J. Zhang, "Relative position finite-time coordinated tracking control of spacecraft formation without velocity measurements," ISA Transactions, vol. 54, pp. 60-74, DOI: 10.1016/j.isatra.2014.08.004, 2015.

[13] C.-ming Zou, H. Shan, H.-xing Li, "Vehicle Cooperative Formation Control Based on Model Prediction," Control Engineering, 2020.

[14] K. LU, M. Wang, L. Dai L, "Vision-Based Fixed-Time Platoon Formation Control of Multi-Agent Systems," Proceeding of the 2020 39th Chinese Control Conference (CCC), .

[15] Z. Miao, H. Zhong, J. Lin, Y. Wang, "Vision-Based Formation Control of Mobile Robots with FOV Constraints and Unknown Feature Depth- ScienceDirect," IEEE Transactions on Control Systems Technology 29.5, vol. 29 no. 5, pp. 2231-2238, 2021.

[16] T. I. Fossen, "Line-of-sight path-following control utilizing an extended Kalman filter for estimation of speed and course over ground from GNSS positions," Journal of Marine Science and Technology, vol. 27 no. 1, pp. 806-813, DOI: 10.1007/s00773-022-00872-y, 2022.

[17] X. Shao, J. Zhang, W. Zhang, "Distributed cooperative surrounding control for mobile robots with uncertainties and aperiodic sampling," IEEE Transactions on Intelligent Transportation Systems,DOI: 10.1109/tits.2022.3189201, 2022.

[18] M. J. TahkTahk, C. S. Park, C. K. Ryoo, "Line-of-Sight guidance laws for formation flight," Journal of Guidance, Control, and Dynamics, vol. 28 no. 4, pp. 708-716, DOI: 10.2514/1.9605, 2005.

[19] D. B. Wilson, A. H. Goktogan, S. Sukkarieh, "A Vision Based Relative Navigation Framework for Formation Flight," .

[20] J. Dunn, R. Tron, "Temporal siamese networks for clutter mitigation applied to vision-based quadcopter formation control," IEEE Robotics and Automation Letters, vol. 6 no. 1, pp. 32-39, DOI: 10.1109/lra.2020.3028056, 2021.

[21] W. Zhang, X. Shao, W. Zhang, J. Qi, H. Li, "Unknown input observer-based appointed-time funnel control for quadrotors," Aerospace Science and Technology, vol. 126,DOI: 10.1016/j.ast.2022.107351, 2022.

[22] M. A. Dehghani, M. B. Menhaj, "Integral sliding mode formation control of fixed-wing unmanned aircraft using seeker as a relative measurement system," Aerospace Science and Technology, vol. 58, pp. 318-327, DOI: 10.1016/j.ast.2016.08.011, 2016.

[23] Z. Zhao, J. Wang, Y. Chen, S. Ju, "Iterative learning-based formation control for multiple quadrotor unmanned aerial vehicles," International Journal of Advanced Robotic Systems, vol. 17 no. 2,DOI: 10.1177/1729881420911520, 2020.

[24] D Wang, "Research on Multi UAV Distributed Formation Maintenance and Reconfiguration Control Method," Tianjin University, 2019.

[25] X. Shao, Y. Shi, W. Zhang, "Neural-network-based constrained output-feedback control for MEMS gyroscopes considering scarce transmission bandwidth," IEEE Transactions on Cybernetics. Systems,DOI: 10.1109/tcyb.2021.3070137, 2021.

[26] Y. N. Jia, X. Y Tian, Q. Li, "Review of Research Progress on UAV Cluster Acta Aeronautica Sinica," vol. 41 no. S1, 2020.

[27] X. K. Wang, "Overview and future development of small fixed wing UAV cluster," Acta Aeronautica Sinica, vol. 41 no. 04, pp. 20-45, 2020.

[28] A. Grancharova, E. I. Johansen, T. A. Johansen, "Rotary-wing UAVs trajectory planning by distributed linear MPC with reconfigurable communication network topologies," IFAC Proceedings Volumes, vol. 46 no. 27, pp. 198-205, DOI: 10.3182/20130925-2-de-4044.00020, 2013.

[29] D. Lapandic, L. Persson, D. V. Dimarogonas, B. Wahlberg, "Aperiodic Communication for MPC in Autonomous Cooperative Landing," IFAC-PapersOnLine, vol. 54 no. 6, pp. 113-118, DOI: 10.1016/j.ifacol.2021.08.532, 2021.

[30] D. Jang, C. Y. Son, J. Yoo, H. J. Kim, K. H. Johansson, "Efficient networked UAV control using event-triggered predictive control," IFAC-PapersOnLine, vol. 52 no. 15, pp. 412-417, DOI: 10.1016/j.ifacol.2019.11.710, 2019.

[31] Z. Dong, X. Xu, X. Zhang, X. Zhou, X. Liu, X. Li, "Real-time Motion Planning Based on MPC with Obstacle Constraint Convexification for Autonomous Ground Vehicles," .

[32] Y. Kuriki, T. Namerikawa, "Formation Control with Collision Avoidance for a Multi-UAV System Using Decentralized MPC and Consensus-Based Control," .

[33] A. Papen, R. Vandenhoeck, J. Bolting, F. Defay, "Collision-free rendezvous maneuvers for formations of unmanned aerial vehicles," IFAC-PapersOnLine, vol. 50 no. 1, pp. 282-289, DOI: 10.1016/j.ifacol.2017.08.047, 2017.