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

Authors: Qi-jie Chen; Yu-qiang Jin; Ting-long Yan; Tao-yu Wang; Yao Wang
Journal: Mathematical Problems in Engineering
 (2022)

---

## 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.


---

## Body

## 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. Section2 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 framesI, Vi, and Li 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 1
Coordinate frames.


### 2.1. Single UAV Motion Model

Assuming that there areNu 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 viVi=vi00T as the follower velocity vector with respect to the frame Vi, we have(1)viI=p˙i=CViIviVi=vicosγicosχivicosγisinχi−visinγi,where CViI=cosγicosχi−sinχisinγicosχicosγisinχicosχisinγisinχi−sinγi0cosγi is the rotation matrix of I with respect to Vi , vi is the speed of the i th follower UAV, and χi is the yaw angles of follower UAV, and γi is the pitch angle of follower UAV. pi=xiyiziT is the position vector of the UAV in the inertial coordinate frame.By denotingDIvi as the velocity differential in inertial frame, ωiVi=ωxiωyiωziT as the angular velocity vector in the follower velocity frame Vi, aiVi=axiViayiViaziViT as the acceleration vector in the follower velocity frame, we have the following Coriolis equation:(2)DIvi=axiViayiViaziVi=DVivi+ωiVi×vi=v˙i00+ωxiωyiωzi×vi00=v˙iωzivi−ωyivi.By denotingxi=xiyiziviχiγiT∈ℝ6 as the states of followers, and by using equations (1) and (2), the motion equation of follower can be express as follows:(3)x˙i=x˙iy˙iz˙iv˙iχ˙iγ˙i=vicosγicosχivicosγisinχi−visinγiaxiViaziVivi−ayiVivi.By consideringui=axiViayiViaziViT∈ℝ3 as the input of the system, we have the expression of the UAV state space as follows:(4)x˙it=fxit,uit,where f· is the corresponding nonlinear function.


### 2.2. Neighborhood UAV State Estimation Model

By denotinglij as the relative distance between UAV i and j, χij and γij are the relative pitch and yaw angles; we have the relative kinematics of the leader-follower, as shown in Figure 2.Figure 2
The relative kinematics of the leader-follower.By assuminglij=lij00T as the vector of position of UAV j with respect to line of sight frame Li, we have the position vector of UAV j in the inertial frame(5)pj=pi+CLiVi·lij,where CLiVi=cosγijcosχijcosγijcosχij−sinγij−sinχijcosχij0sinγijcosχijsinγijcosχijcosγij. By introducing ωIjLi=ωxijωyijωzijT as the relative angular velocity in the frame Li, we have(6)vijLi=DIlij=DLilij+ωIjLi×lij=l˙ijlijωzij−lijωyij,where vijLi=vxijvyijvzijT. Moreover, the velocity of UAV j can be described as follows:(7)vjVi=CLiVi⋅vijLi+viVi.Then, the velocity of UAVj in frame I can be described as follows:(8)vjVi=CViI⋅vjVi,where · is the Euclidean norm of vector. By taking the time derivative of equation (6), we can get(9)aijLi=DIvij=DLivij+ωIjLi×vij,=l¨ij−lijωzij2+ωyij22l˙ijωzij+lijω˙zij+lijωxijωyij−2l˙ijωyij−lijω˙yij+lijωxijωzij,where aijLi=axijayijazijT. Moreover, the angular velocity of Li can be given as follows:(10)ωIjLi=ωIiLi+ωijLi=CViLiωiVi+ωiLi=ωxicosγijcosχij+ωyicosγijsinχij−ωzisinγij−ωxisinγij+ωyicosχijωxicosγijsinχij+ωyisinγijsinχij+ωzicosγij+−χ˙ijsinγijγijχ˙ijcosγij.By assumingωxi=0 which means that the UAV does not roll to simplify the model, we have(11)ωIjLi=ωyicosγijsinχij−ωzisinγij−χ˙ijsinγijωyicosχij+γijωyisinγijsinχij+ωzicosγij+χ˙ijcosγij.By considering equations (9) and (11), we can have the acceleration of UAV j in the frame Vi as follows:(12)u∧jVi=aiVi+CLiVi⋅a∧ijLi,where u∧jVi=axjViayjViazjViT∈ℝ3. Then, the state estimator of the UAV j can be described as follows:(13)x˙∧j=CViI⋅vjViCvj⋅u∧jVi,where Cvj=10000−1/v∧j01/v∧j0. 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.


## 2.1. Single UAV Motion Model

Assuming that there areNu 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 viVi=vi00T as the follower velocity vector with respect to the frame Vi, we have(1)viI=p˙i=CViIviVi=vicosγicosχivicosγisinχi−visinγi,where CViI=cosγicosχi−sinχisinγicosχicosγisinχicosχisinγisinχi−sinγi0cosγi is the rotation matrix of I with respect to Vi , vi is the speed of the i th follower UAV, and χi is the yaw angles of follower UAV, and γi is the pitch angle of follower UAV. pi=xiyiziT is the position vector of the UAV in the inertial coordinate frame.By denotingDIvi as the velocity differential in inertial frame, ωiVi=ωxiωyiωziT as the angular velocity vector in the follower velocity frame Vi, aiVi=axiViayiViaziViT as the acceleration vector in the follower velocity frame, we have the following Coriolis equation:(2)DIvi=axiViayiViaziVi=DVivi+ωiVi×vi=v˙i00+ωxiωyiωzi×vi00=v˙iωzivi−ωyivi.By denotingxi=xiyiziviχiγiT∈ℝ6 as the states of followers, and by using equations (1) and (2), the motion equation of follower can be express as follows:(3)x˙i=x˙iy˙iz˙iv˙iχ˙iγ˙i=vicosγicosχivicosγisinχi−visinγiaxiViaziVivi−ayiVivi.By consideringui=axiViayiViaziViT∈ℝ3 as the input of the system, we have the expression of the UAV state space as follows:(4)x˙it=fxit,uit,where f· is the corresponding nonlinear function.


## 2.2. Neighborhood UAV State Estimation Model

By denotinglij as the relative distance between UAV i and j, χij and γij are the relative pitch and yaw angles; we have the relative kinematics of the leader-follower, as shown in Figure 2.Figure 2
The relative kinematics of the leader-follower.By assuminglij=lij00T as the vector of position of UAV j with respect to line of sight frame Li, we have the position vector of UAV j in the inertial frame(5)pj=pi+CLiVi·lij,where CLiVi=cosγijcosχijcosγijcosχij−sinγij−sinχijcosχij0sinγijcosχijsinγijcosχijcosγij. By introducing ωIjLi=ωxijωyijωzijT as the relative angular velocity in the frame Li, we have(6)vijLi=DIlij=DLilij+ωIjLi×lij=l˙ijlijωzij−lijωyij,where vijLi=vxijvyijvzijT. Moreover, the velocity of UAV j can be described as follows:(7)vjVi=CLiVi⋅vijLi+viVi.Then, the velocity of UAVj in frame I can be described as follows:(8)vjVi=CViI⋅vjVi,where · is the Euclidean norm of vector. By taking the time derivative of equation (6), we can get(9)aijLi=DIvij=DLivij+ωIjLi×vij,=l¨ij−lijωzij2+ωyij22l˙ijωzij+lijω˙zij+lijωxijωyij−2l˙ijωyij−lijω˙yij+lijωxijωzij,where aijLi=axijayijazijT. Moreover, the angular velocity of Li can be given as follows:(10)ωIjLi=ωIiLi+ωijLi=CViLiωiVi+ωiLi=ωxicosγijcosχij+ωyicosγijsinχij−ωzisinγij−ωxisinγij+ωyicosχijωxicosγijsinχij+ωyisinγijsinχij+ωzicosγij+−χ˙ijsinγijγijχ˙ijcosγij.By assumingωxi=0 which means that the UAV does not roll to simplify the model, we have(11)ωIjLi=ωyicosγijsinχij−ωzisinγij−χ˙ijsinγijωyicosχij+γijωyisinγijsinχij+ωzicosγij+χ˙ijcosγij.By considering equations (9) and (11), we can have the acceleration of UAV j in the frame Vi as follows:(12)u∧jVi=aiVi+CLiVi⋅a∧ijLi,where u∧jVi=axjViayjViazjViT∈ℝ3. Then, the state estimator of the UAV j can be described as follows:(13)x˙∧j=CViI⋅vjViCvj⋅u∧jVi,where Cvj=10000−1/v∧j01/v∧j0. 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 xref=xrefyrefzrefvrefχrefγrefT∈ℝ6 as the reference state of UAV i, the error between UAV i and its reference state in the frame Vi can be described as follows:(14)eir=CIVi00I3·xi−xref.By taking the time derivative of above equation, one can get(15)e˙ir=ωyi⋅C˙γi0003⋅Cχi00I3⋅xi−xref+ωzi⋅Cγi0003⋅C˙χi00I3⋅xi−xref+Cγi00I3⋅Cχi00I3⋅x˙i−x˙ref,where Cχi=cosχisinχi0−sinχicosχi0001, Cγi=cosγi0−sinγi010sinγi0cosγi. By using equations (3) and (11), the error state can be expressed as follows:(16)e˙ir=−ωyisinγicosχi+ωzicosγisinχi−ωyisinγisinχi+ωzicosγicosχi−ωyicosγi000−ωzicosχi−ωzisinχi0000ωyicosγicosχi−ωzisinγisinχiωyicosγisinχi+ωzisinγicosχi−ωyisinγi000000000000000000000exeyezeveχeγ+vi−vrefcosγrefcosγicoseχ−vrefsinγrefsinγivrefcosγrefsineχ−vrefcosγrefsinγicoseχ+vrefsinγrefcosγiaxi−axref−azivi+azref/vrefayivi−ayref/vref.To simplify the calculation, define a virtual input:(17)u˜ik=u˜ixu˜iyu˜iz=axi−axrefωzi−ωzrefωyi−ωyref=axi−axref−azivi+azref/vref−ayivi−ayref/vrefayi.In actual operation, the input can be calculated according to equation (17). By linearizing equation (16) at eir0=0,e˙ir0=0, we have(18)e˙irk=Aeirk+Bu˜ik,where(19)A=−ωyrefsinγrefcosχref+ωzrefcosγrefsinχref−ωyrefsinγrefsinχref+ωzrefcosγrefcosχref−ωyrefcosγref100−ωzrefcosχref−ωzrefsinχref00vrefcosγref0ωyrefcosγrefcosχref−ωzrefsinγrefsinχrefωyrefcosγrefsinχref+ωzrefsinγrefcosχref−ωyrefsinγref00−vref000000000000000000,B=030303I3.


### 3.2. Formation Control Design

The DMPC-based control design can be described as, at timek, the UAV i solves the optimal control problem with the predicted number of steps starting from the initial error eirk. 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

(20)minuitJlk=minuit∑l1=0N−1eirk+l|kPi+uik+l|k−urefk+l|kRi+∑l2=0N−2Δuik+l|kQi+eirk+N|kGis.t.x˙ik=fxik,uik,k=k0,k0+1,…,k0+N−1,e˙irt=geirt,u˜it,xik∈Xi,uik∈Ui.


#### 3.2.2. Cost Function of Followers

(21)minuitJfifk=minuit∑l1=0N−1∑j∈Nifeijk+l|kPi+∑j∈Nifuik+l|k−ujk+l|kRi+∑l2=0N−2Δuik+l|kQi+∑j∈Nieijk+N|kGis.t.x˙ik=fxik,uik,k=k0,k0+1,...,k0+N−1e˙irt=geirt,u˜itxik∈Xiuik∈Ui,where xP=xTPx1/2 represents the quadratic form of the vector. Nf represents the set of followers. Ni represents the set of UAVs in the neighborhood of the UAV i. eirk+l|k indicates the error state which the UAV i starts from time k and applies the control sequence uik+l|k in turn to obtain within the number of predicted steps. Xi is the set of state constraints of UAV i which can be expressed as follows: Xi=xik|vik≤vmax,χik≤χmax,γik≤γmax. Ui is the input constraint set of the system, which can be described as follows:(22)Ui=aix≤axmax,aiy≤aymax,aiz≤azmax,uik|Δaix≤Δaxmax,Δaiy≤Δaymax,Δaiz≤Δazmax,where axmaxaymaxazmax∈ℝ3 indicates the maximum value of control input in three directions and ΔaixΔaiyΔaiz∈ℝ3 indicates the maximum value of control input change in three directions. The weight matrix satisfies Pi>0,Qi>0,Ri>0,Gi>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 Figure3.Figure 3
Schematic diagram of algorithm principle.The process of DMPC can be described as follows:Step 1: Bind UAV number and priority. Bind the reference to UAV1.Step 2: Att=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 eirk|k=eirt=k.Step 3: Calculate the optimal input sequenceUi∗k which can be described as Ui∗k=ui∗k|k,ui∗k+1|k,...,ui∗k+N−1|k according to equations (20) and (21). The first item of the optimal input sequence ui∗k|k is used as the input of the UAV at t=k+1.Step 4: Att=k+1, based on the new error state eirt=k+1 of the formation, repeat to Step 2


## 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 xref=xrefyrefzrefvrefχrefγrefT∈ℝ6 as the reference state of UAV i, the error between UAV i and its reference state in the frame Vi can be described as follows:(14)eir=CIVi00I3·xi−xref.By taking the time derivative of above equation, one can get(15)e˙ir=ωyi⋅C˙γi0003⋅Cχi00I3⋅xi−xref+ωzi⋅Cγi0003⋅C˙χi00I3⋅xi−xref+Cγi00I3⋅Cχi00I3⋅x˙i−x˙ref,where Cχi=cosχisinχi0−sinχicosχi0001, Cγi=cosγi0−sinγi010sinγi0cosγi. By using equations (3) and (11), the error state can be expressed as follows:(16)e˙ir=−ωyisinγicosχi+ωzicosγisinχi−ωyisinγisinχi+ωzicosγicosχi−ωyicosγi000−ωzicosχi−ωzisinχi0000ωyicosγicosχi−ωzisinγisinχiωyicosγisinχi+ωzisinγicosχi−ωyisinγi000000000000000000000exeyezeveχeγ+vi−vrefcosγrefcosγicoseχ−vrefsinγrefsinγivrefcosγrefsineχ−vrefcosγrefsinγicoseχ+vrefsinγrefcosγiaxi−axref−azivi+azref/vrefayivi−ayref/vref.To simplify the calculation, define a virtual input:(17)u˜ik=u˜ixu˜iyu˜iz=axi−axrefωzi−ωzrefωyi−ωyref=axi−axref−azivi+azref/vref−ayivi−ayref/vrefayi.In actual operation, the input can be calculated according to equation (17). By linearizing equation (16) at eir0=0,e˙ir0=0, we have(18)e˙irk=Aeirk+Bu˜ik,where(19)A=−ωyrefsinγrefcosχref+ωzrefcosγrefsinχref−ωyrefsinγrefsinχref+ωzrefcosγrefcosχref−ωyrefcosγref100−ωzrefcosχref−ωzrefsinχref00vrefcosγref0ωyrefcosγrefcosχref−ωzrefsinγrefsinχrefωyrefcosγrefsinχref+ωzrefsinγrefcosχref−ωyrefsinγref00−vref000000000000000000,B=030303I3.


## 3.2. Formation Control Design

The DMPC-based control design can be described as, at timek, the UAV i solves the optimal control problem with the predicted number of steps starting from the initial error eirk. 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

(20)minuitJlk=minuit∑l1=0N−1eirk+l|kPi+uik+l|k−urefk+l|kRi+∑l2=0N−2Δuik+l|kQi+eirk+N|kGis.t.x˙ik=fxik,uik,k=k0,k0+1,…,k0+N−1,e˙irt=geirt,u˜it,xik∈Xi,uik∈Ui.


### 3.2.2. Cost Function of Followers

(21)minuitJfifk=minuit∑l1=0N−1∑j∈Nifeijk+l|kPi+∑j∈Nifuik+l|k−ujk+l|kRi+∑l2=0N−2Δuik+l|kQi+∑j∈Nieijk+N|kGis.t.x˙ik=fxik,uik,k=k0,k0+1,...,k0+N−1e˙irt=geirt,u˜itxik∈Xiuik∈Ui,where xP=xTPx1/2 represents the quadratic form of the vector. Nf represents the set of followers. Ni represents the set of UAVs in the neighborhood of the UAV i. eirk+l|k indicates the error state which the UAV i starts from time k and applies the control sequence uik+l|k in turn to obtain within the number of predicted steps. Xi is the set of state constraints of UAV i which can be expressed as follows: Xi=xik|vik≤vmax,χik≤χmax,γik≤γmax. Ui is the input constraint set of the system, which can be described as follows:(22)Ui=aix≤axmax,aiy≤aymax,aiz≤azmax,uik|Δaix≤Δaxmax,Δaiy≤Δaymax,Δaiz≤Δazmax,where axmaxaymaxazmax∈ℝ3 indicates the maximum value of control input in three directions and ΔaixΔaiyΔaiz∈ℝ3 indicates the maximum value of control input change in three directions. The weight matrix satisfies Pi>0,Qi>0,Ri>0,Gi>0.


## 3.2.1. Cost Function of Leader

(20)minuitJlk=minuit∑l1=0N−1eirk+l|kPi+uik+l|k−urefk+l|kRi+∑l2=0N−2Δuik+l|kQi+eirk+N|kGis.t.x˙ik=fxik,uik,k=k0,k0+1,…,k0+N−1,e˙irt=geirt,u˜it,xik∈Xi,uik∈Ui.


## 3.2.2. Cost Function of Followers

(21)minuitJfifk=minuit∑l1=0N−1∑j∈Nifeijk+l|kPi+∑j∈Nifuik+l|k−ujk+l|kRi+∑l2=0N−2Δuik+l|kQi+∑j∈Nieijk+N|kGis.t.x˙ik=fxik,uik,k=k0,k0+1,...,k0+N−1e˙irt=geirt,u˜itxik∈Xiuik∈Ui,where xP=xTPx1/2 represents the quadratic form of the vector. Nf represents the set of followers. Ni represents the set of UAVs in the neighborhood of the UAV i. eirk+l|k indicates the error state which the UAV i starts from time k and applies the control sequence uik+l|k in turn to obtain within the number of predicted steps. Xi is the set of state constraints of UAV i which can be expressed as follows: Xi=xik|vik≤vmax,χik≤χmax,γik≤γmax. Ui is the input constraint set of the system, which can be described as follows:(22)Ui=aix≤axmax,aiy≤aymax,aiz≤azmax,uik|Δaix≤Δaxmax,Δaiy≤Δaymax,Δaiz≤Δazmax,where axmaxaymaxazmax∈ℝ3 indicates the maximum value of control input in three directions and ΔaixΔaiyΔaiz∈ℝ3 indicates the maximum value of control input change in three directions. The weight matrix satisfies Pi>0,Qi>0,Ri>0,Gi>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 Figure3.Figure 3
Schematic diagram of algorithm principle.The process of DMPC can be described as follows:Step 1: Bind UAV number and priority. Bind the reference to UAV1.Step 2: Att=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 eirk|k=eirt=k.Step 3: Calculate the optimal input sequenceUi∗k which can be described as Ui∗k=ui∗k|k,ui∗k+1|k,...,ui∗k+N−1|k according to equations (20) and (21). The first item of the optimal input sequence ui∗k|k is used as the input of the UAV at t=k+1.Step 4: Att=k+1, based on the new error state eirt=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 existsk1,k2∈ℝ, k1>vref,k2>vref such that the virtual input satisfies(23)u˜ik=−A12⋅e1−k2e2−A12−1A11+k1⋅IA11⋅e1+A12⋅e2.
System equation (18) is asymptotically stable, where e1=exeyezT, e1=eveχeγT,(24)A11=−ωyrefsinγrefcosχref+ωzrefcosγrefsinχref−ωyrefsinγrefsinχref+ωzrefcosγrefcosχref−ωyrefcosγref−ωzrefcosχref−ωzrefsinχref0ωyrefcosγrefcosχref−ωzrefsinγrefsinχrefωyrefcosγrefsinχref+ωzrefsinγrefcosχref−ωyrefsinγref,A12=1000vrefcosγref000−vref.Proof.
From equation (18), we have(25)e˙1e˙2=A11A120303e1e2+030303B2203u˜i.
We can rewrite the above equation in the following form:(26)e˙1e˙2=A11⋅e1+A12⋅e2B22⋅u˜i,where B22=I. By assuming e1d as the expected value of e1, we can describe the error δ1 as follows:(27)δ1=e1d−e1.
By taking the time derivative of the above equation, one can get(28)δ˙1=e˙1d−e˙1.
Choose the Lyapunov function as(29)V1=12δ1T⋅δ1.
By considering equations (26) and (27), we have(30)V˙1=δ1T⋅δ˙1=δ1T⋅e˙1d−A11e1−A12e2.
SinceV˙1=δ1⋅δ˙1<0, we have δ˙1=e˙1d−A11e1−A12e2⟶−k1⋅e1, where k1>0. By assuming e2d as the expected value of e2, we have(31)e2d=A12−1⋅e˙1d−A11⋅e1+k1⋅δ1.
The errorδ2 for e2 can be expressed as(32)δ2=e2d−e2.
Substituting equations (31) and (32) into (30), we have(33)V˙1=δ1T⋅e˙1d−A11e1−A12e2=δ1T⋅e˙1d−A11e1−A12⋅A12−1⋅e˙1d−A11⋅e1+k1⋅δ1−δ2=−k1⋅δ1T⋅δ1+A12⋅δ1T⋅δ2.
By substituting equations (26) and (31) into (32) and taking the time derivative of it, we have(34)δ˙2=e˙2d−e˙2=A12−1e¨1d+k1e˙1d−A11+k1⋅IA11⋅e1+A12⋅e2−B22⋅u˜i.
Let the Lyapunov function withδ2 be(35)V2=V1+12δ2T⋅δ2.
By taking the time derivative of equation (35) and substituting equation (33) into it, we have(36)V˙2=−k1⋅δ1T⋅δ1+δ1TA12⋅δ1+δ˙2,where k2>0, By substituting equations (26), (32), and (34) into the above equation, one can get the following equation:(37)A12⋅δ1+A12−1e¨1d+k1e˙1d−A11+k1A11⋅e1+A12⋅e2−B22⋅u˜i=−k2δ2.
By rewriting the above equation, we can get the equation of the input:(38)u˜i=B22−1A12⋅δ1+A12−1e¨1d+k1e˙1d−A11+k1A11⋅e1+A12⋅e2+k2δ2.
By substituting equations (38) into (34), we have(39)δ˙2=e˙2d−e˙2=A12−1e¨1d+k1e˙1d−A11+k1⋅IA11⋅e1+A12⋅e2−B22⋅B22−1δ1+A12−1e¨1d+k1e˙1d−A11+k1A11⋅e1+A12⋅e2+k2δ2=−A12⋅δ1−k2⋅δ2.
By substituting equations (26) and (31) into (28), one can get(40)δ˙1δ˙2=−k1⋅IA12−A12−k2⋅Iδ1δ2.
Let the eigenvalues of system equation (40) be λi,i=1,2,...,6, corresponding to the eigenvalues of the 1st to 6th rows of the eigenmatrix. We have(41)∏i=16λi=−k1⋅I⋅−k2⋅I−−A12⋅A12=k1k23+A122∑i=16λi=−3⋅k1+k2,where A represents the determinant of A. Since k1>0, k2>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(42)λ1+k1≤1λ2+k1≤vrefcosγrefλ3+k1≤vrefλ4+k2≤1λ5+k2≤vrefcosγrefλ6+k1≤vref.
By considering equation (42), when λi>maxvref,1,i=1,2,...,6, the eigenvalues of system equation (40) are all in the negative half-plane and are asymptotically stable. By considering e1d=e˙1d=e¨1d=0 and e2d=e˙2d=0, and substituting equations (27) and (32) into (38), we have(43)u˜i=−A12⋅e1−k2e2−A12−1A11+k1⋅IA11⋅e1+A12⋅e2.
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 Ui∗k can be transformed into the solution of the virtual input sequence U˜ik. By rewriting equations (20) and (21) with the virtual input u˜ik in equation (17), we have a new expression for the cost function as follows:(44)minuitJlk=minuit∑l1=0N−1eirk+l|kPi+u˜ik+l|kRi+∑l2=0N−2Δu˜ik+l|kQi+eirk+N|kGis.t.x˙ik=fxik,uik,k=k0,k0+1,…,k0+N−1e˙irt=geirt,u˜itxik∈Xiuik∈Ui,(45)minuitJfifk=minuit∑l1=0N−1∑j∈Nieijk+l|kPi+∑j∈Niu˜ik+l|kRi+∑l2=0N−2Δu˜ik+l|kQi+∑j∈Nieijk+N|kGis.t.x˙ik=fxik,uik,k=k0,k0+1,…,k0+N−1e˙irt=geirt,u˜itxik∈Xiu˜ik∈U˜i.
Att=k, the optimal solution sequence of equation (45) is U˜ik=u˜i∗k+l|k,l=0,..,N−1. The optimal state sequence corresponding to the solution is Xik=xi∗k+l|k,l=0,..,N−1. The error state sequence corresponding to the solution is ei∗k=ei∗k+l|k,l=0,..,N−1. Obviously, all components of the optimal state sequence Xik satisfy the constraints.
Att=k+1, we construct the solution of equation (45) which can be expressed as follows:(46)u˜ik+1=u˜i∗k+l|k,l=0,1,2,…,N−2Ce1k+N|k+De2k+N|k,l=N−1,where C=−A12+A12−1A11+k1⋅IA11, D=−k2⋅I+A12−1A11+k1⋅IA12. By considering equation (45), we can get the state sequence and error state sequence of the system expressed as follows:(47)Xik+1=xik+1+l|k,l=1,…,N=xi∗k+2|k,xi∗k+3|k,…,xi∗k+N|k,xik+N|k+1eik+1=eirk+1+l|k,l=0,…,N−1=eirk+2|k,eirk+2|k,…,eirk+N|k,eirk+N|k+1.
By considering the firstN−1 of Xik+1 are the components of Xik, the first N−1 of eik+1 are the components of ei∗k, all of them satisfy the constraints. By substituting xi∗k+N|k=xref∗k+N|k into (10), one can get(48)eir∗k+N|k=xi∗k+N|k−xrefk+N|k=06×1.
By considering equation (40) and the input which can be expressed as Ceir1∗k+N|k+Deir2∗k+N|k, we have(49)eir1k+N|k+1−eir1∗k+N|keir2k+N|k+1−eir2∗k+N|k=e˙ir1k+N|k+1e˙ir2k+N|k+1=−k1⋅IA12−A12−k2⋅Ieir1∗k+N|keir2∗k+N|k=06×1.
It can be rewriting aseirk+N|k=xik+N|k+1−xrefk+N|k=06×1. By considering all the above equations, we can get that Xik+1 meet terminal constraints and u˜ik+1 is a set of feasible solutions to equation (45). Similarly, we can get a set of feasible solutions u˜lk+1 to equation (44). At t=k, the cost functions of leader and followers can be expressed as follows:(50)Jlk=∑l1=0N−1elrk+l|kPi+u˜lk+l|kRi+∑l2=0N−2Δulk+l|kQi+elrk+N|kGi,Jfifk=∑l1=0N−1∑j∈Nieijk+l|kPi+∑j∈Niu˜ik+l|kRi+∑l2=0N−2Δu˜ik+l|kQi+∑j∈Nieijk+N|kGi.
Att=k+1,(51)Jlk+1=∑l1=0N−1elrk+l+1|k+1Pi+u˜lk+l+1|k+1Ri+∑l2=0N−2Δulrk+l+1|k+1Qi+elrk+l+1|k+1Gi=Jl∗k−elrk|kPi−u˜lk|kRi−Δulrk|kQi+elrk+N|kPi+u˜lk+N|kRi+Δulrk+N|kQi+elrk+N+1|kGi,Jfifk+1=∑l1=0N−1∑j∈Nieijk+l+1|k+1Pi+∑j∈Niu˜ik+l+1|k+1Ri+∑l2=0N−2Δu˜ijk+l+1|k+1Qi+∑j∈Nieijk+l+1|k+1Gi=Jfif∗k+1−∑j∈Nieijk|kPi−∑j∈Niu˜ik|kRi−∑j∈NiΔu˜ik|kQi+∑j∈Nieijk+N|k+1Pi+∑j∈Niu˜ik+N|kRi+∑j∈NiΔu˜ijk+N|kQi+∑j∈Nieijk+N+1|kGi.
By consideringeirk+N|kPi=eirk+N|kGi=0 and Δu˜irk+N|kQi=Δu˜lrk+N|kQi=0, we can rewriting the above equations as follows:(52)Jlk+1=Jl∗k−elrk|kPi−u˜lk|kRi−Δulrk|kQi≤Jl∗k,Jfifk+1=Jfif∗k+1−∑j∈Nieijk|kPi−∑j∈Niu˜ik|kRi−∑j∈NiΔu˜ik|kQi≤Jfif∗k.
By subtracting the cost function at timek from both the sides at the same time, we can get(53)Jl∗k+1−Jl∗k≤0,Jfif∗k+1−Jfif∗k≤0.
It can be rewritten as(54)Jl∗k+1+∑Jfif∗k+1−Jl∗k+∑Jfif∗k≤0,where the equal sign holds if and only if elrk|kPi=∑j∈Nieijk|kPi=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 Figure4, and the arrows indicate the flow of information.Figure 4
UAV formation.Table1 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/s2. 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.
UAVInitial positionvm/sDesired relative positionReference line550T—1220T8000T2−230T8−5−50T32−10T85−50T4−2−20T80−100TCase 1.
Straight Flight
Figure5 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 5
Result of straight flight. (a) Formation trajectories. (b) Distance between UAVs. (c) Velocity of each UAV. (d) Control trajectories of each UAV.
(a)(b)(c)(d)Case 2.
Fixed Turning Radius Flight
Figure6 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 6
Result of fixed turning radius flight. (a) Formation trajectories. (b) Distance between UAV. (c) Velocity of each UAV. (d) Control trajectories of each UAV.
(a)(b)(c)(d)Case 3.
Variable Turning Radius Flight
Figure7 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 7
Result of variable turning radius flight. (a) Formation trajectories. (b) Distance between UAV. (c) Velocity of each UAV. (d) Control trajectories of each UAV.
(a)(b)(c)(d)Case 4.
Multi-UAVs. (Table2)
Unlike the previous cases with four UAVs, Case4 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.
UAVInitial positionvm/sDesired relative positionReference line550T—1220T8000T2−230T8−5−50T32−10T80−50T4−2−20T85−50T5240T8−5−100T6−4−40T80−100T74−20T85−100T8−4−40T80−150TFigure 8
Result of multi-UAVs. (a) Formation trajectories. (b) Distance between UAV. (c) Velocity of each UAV. (d) Control trajectories of each UAV.
(a)(b)(c)(d)


## 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.



---

*Source: 7316009-2022-09-20.xml*