# Multimodel Predictive Control Approach for UAV Formation Flight

Authors: Chang-jian Ru; Rui-xuan Wei; Ying-ying Wang; Jun Che
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Formation flight problem is the most important and interesting problem of multiple UAVs (unmanned aerial vehicles) cooperative control. In this paper, a novel approach for UAV formation flight based on multimodel predictive control is designed. Firstly, the state equation of relative motion is obtained and then discretized. By the geometrical method, the characteristic points of state are determined. Afterwards, based on the linearization technique, the standard linear discrete model is obtained at each characteristic state point. Then, weighted model set is proposed using the idea of T-S (Takagi-Sugeno) fuzzy control and the predictive control is carried out based on the multimodel method. Finally, to verify the performance of the proposed method, two different simulation scenarios are performed.


---

## Body

## 1. Introduction

In recent years, as an advanced system with high autonomy, UAVs have been widely applied in the fields of both civilian and military. When a single UAV accomplishes tasks individually, it will be more likely to reduce mission success, due to its limited information accessing ability. In comparison, multiple UAVs collaborating with each other maintain a certain formation during the flight, which provides them with full access to environmental information, increases resistance to external attack capability, improves working efficiency and robustness of the system, and so forth, so it has attracted wide attention [1–3]. Formation flight is an important aspect of multiple UAVs cooperative control. When the formation shape maintains or changes according to the mission requirements, it is necessary to control the relative position, attitude, and speed between UAVs, and so forth. However, UAV’s control system is a nonlinear coupling system, coupled with complex operational environment constraints, putting forward higher design requirement for formation controller. So it is essential to propose an effective control strategy to solve those problems.Model predictive control (MPC) method is put forward by some scholars. Here, several typical researches are presented. A hierarchical approach and a set of MPC strategies for the UAV formation are proposed in [4], where obstacle and collision avoidance constraints are taken into account. A distributed collision-free formation flight control law in the framework of nonlinear model predictive control is designed in [5]. In [6], a dual mode MPC method is used for formation control. To guarantee the stability, the dual mode controller must switch from an MPC control to a terminal state controller. A simple nonlinear model predictive control (NMPC) formulation is used to adequately address the terrain avoidance problem, as presented in [7]. An online nonlinear model predictive control framework is used for the trajectory tracking of autonomous vehicles in [8], where a bicycle model is used for the prediction of future states in the NMPC framework. The validation of a formation flight control technique with obstacle avoidance capability based on nonlinear model predictive algorithms is proposed in [9].The nonlinear model predictive control method provides an effective means to solve the control problem of nonlinear systems [10, 11]. Because the close relative distance between the UAVs may lead to collision; thus, it requires higher control accuracy. However, stair-like MPC uses the way of constraining the variation of the future control quantity, which restricts the maneuverability of the vehicle and is prone to causing collision between the UAVs due to the overshoot problem [12]. So it is necessary to adopt a new predictive control method to achieve the formation flight control. For some more complex systems, the multimodel control method has stronger robustness and higher control accuracy under certain conditions [13]. Besides, multimodel control method can provide the nonlinear system with transparent model and controller, facilitating the system analysis. Compared with other nonlinear global strategies, the multimodel control method cannot greatly reduce computational complexity, but the model and structure of controller are more suitable for online adjustments and learning algorithm [14], so multiple model-based predictive control can be used to solve UAV formation control problem.This paper is organized as follows. In Section2, the discrete relative motion equations for UAV formation are established. In Section 3, a multiple models-based predictive control approach is used to design controller of the formation. Simulation results are given in Section 4. Finally, Section 5 concludes the paper.


## 2. UAV Formation Flight Control Model

### 2.1. Kinematics Model of UAV Formation Flight Control

Assume that during the formation flight, an UAV is flying horizontally and has no sideslip. In the geographic coordinate system, the relationship between the position vectors of leader UAV (leader) and follower UAV (follower) is shown in Figure1.Figure 1
The position relationship between two vehicles.From Figure1, it is easy to obtain the following equation:
(1)D→L=D→W+D→,
where, D→L, D→W are displacement vectors of two vehicles and D→ is the relative displacement vectors between two vehicles.Differentiating (1), one can obtain
(2)dD→Ldt=dD→Wdt+dD→dt.According to the relationship between the moving coordinate system, one can easily obtain(3)dD→Ldt=δD→Wδt+ωw×D→,
where ωw is the yaw angular rate.Since the vehicle is supposed to fly horizontally, the equation of motion will be(4)x˙i=vicos⁡(φi),y˙i=visin(φi),φ˙i=ωi,i=l,w,
where the subscript l and w denote leader UAV and follower UAV, respectively.Combing (3) and (4), the relative motion equation of two vehicles can be obtained as
(5)vl[cos⁡φlsinφl]=vw[cos⁡φwsinφw]+CwL([x˙dy˙d]+[-ωw·ydωw·xd]),
where xd and yd areX-axis value andY-axis value of the distance between two vehicles in the tack coordinates of leader UAV, respectively, and CwL is the coordinate transformation matrix.Then, carry on the translational process, and one can obtain(6)[x˙dy˙d]=CLW(vl[cos⁡φlsinφl]-vw[cos⁡φwsinφw])-[-ωw·ydωw·xd],
where
(7)CLWvl[cos⁡φlsinφl]=vl[cos⁡φesinφe],CLWvw[cos⁡φwsinφw]=[vw0].So (6) can be written as follows:
(8)[x˙dy˙d]=vl[cos⁡(φl-φw)sin(φl-φw)]-[vw0]-[-ωw·ydωw·xd].Additionally, there is(9)φ˙w=ωw.Combing (8) and (9), we can obtain
(10)x˙d=vlcos⁡(φl-φw)-vw+ωw·yd,y˙d=vlsin(φl-φw)-ωw·xd,φ˙w=ωw.Then, (10) can also be written as the state equations which is shown as follows
(11)[x˙dy˙dφ˙w]=[-1yd0-xd01][vwωw]+[vlcos⁡(φl-φw)vlsin(φl-φw)0].The output equation is(12)y1=xd,y2=yd.


### 2.2. Discrete Model of UAV Formation Movement and Its Predictive Control Analysis

In the previous section, the state equation of formation control is obtained. Here, discretize this equation and the following equation can be obtained(13)[xd(k+1)yd(k+1)φw(k+1)]=[xd(k)yd(k)φw(k)]+[-1yd(k)0-xd(k)01][vw(k)ωw(k)]ΔT+[vl(k)cos⁡⁡(φl(k)-φw(k))vl(k)sin⁡(φl(k)-φw(k))0]ΔT.Generally, the sampling periodic time is short during receding optimization process, so the velocity and yaw angle of leader UAV can be considered constant in sampling period [15], which means in a short sampling period, there is
(14)vl(k+i)=vl(k+i-1)=vl(k),φl(k+i)=φl(k+i-1)=φl(k).According to these two equations, the predicted value of outputs will be obtained as follows:(15)xdp(k+1),xdp(k+2),…,xdp(k+N),ydp(k+1),ydp(k+2),…,ydp(k+N).Since at timek, yd(k), vl(k), φw(k), φl(k) are known, the future output values of formation flight control are xdp(k+1),xdp(k+2),…,xdp(k+N-1),xdp(k+N), and ydp(k+1),ydp(k+2),…,ydp(k+N-1),ydp(k+N), and these outputs are only the function of the future control quantities vw(k),vw(k+1),…,vw(k+N-1) and ωw(k),ωw(k+1),…,ωw(k+N-1). Obviously, this function is a multi-input-multioutput nonlinear control problem. So, these values can be obtained using the receding optimization algorithm. Then, using vw(k) and ωw(k) as outputs of control quantity, and carry out the receding optimization algorithm in sequence, the input of control quantity in the next time can be obtained.According to Section1, it can be known that for the problem of UAV formation flight, nonlinear model predictive control and fuzzy stair-like predictive control have some limitations, but the multimodel control method has stronger robustness and higher control accuracy, the final predictive model of which is linear, so the receding optimization problem can be changed from general nonlinear optimization problem to a linear quadratic optimization problem. Since the linear quadratic optimization has faster computation than the ordinary nonlinear optimization, the multiple models can greatly improve real-time of receding optimization [16]. Therefore, multiple model-based predictive control approach is used to design the controller of UAV formation flight.


## 2.1. Kinematics Model of UAV Formation Flight Control

Assume that during the formation flight, an UAV is flying horizontally and has no sideslip. In the geographic coordinate system, the relationship between the position vectors of leader UAV (leader) and follower UAV (follower) is shown in Figure1.Figure 1
The position relationship between two vehicles.From Figure1, it is easy to obtain the following equation:
(1)D→L=D→W+D→,
where, D→L, D→W are displacement vectors of two vehicles and D→ is the relative displacement vectors between two vehicles.Differentiating (1), one can obtain
(2)dD→Ldt=dD→Wdt+dD→dt.According to the relationship between the moving coordinate system, one can easily obtain(3)dD→Ldt=δD→Wδt+ωw×D→,
where ωw is the yaw angular rate.Since the vehicle is supposed to fly horizontally, the equation of motion will be(4)x˙i=vicos⁡(φi),y˙i=visin(φi),φ˙i=ωi,i=l,w,
where the subscript l and w denote leader UAV and follower UAV, respectively.Combing (3) and (4), the relative motion equation of two vehicles can be obtained as
(5)vl[cos⁡φlsinφl]=vw[cos⁡φwsinφw]+CwL([x˙dy˙d]+[-ωw·ydωw·xd]),
where xd and yd areX-axis value andY-axis value of the distance between two vehicles in the tack coordinates of leader UAV, respectively, and CwL is the coordinate transformation matrix.Then, carry on the translational process, and one can obtain(6)[x˙dy˙d]=CLW(vl[cos⁡φlsinφl]-vw[cos⁡φwsinφw])-[-ωw·ydωw·xd],
where
(7)CLWvl[cos⁡φlsinφl]=vl[cos⁡φesinφe],CLWvw[cos⁡φwsinφw]=[vw0].So (6) can be written as follows:
(8)[x˙dy˙d]=vl[cos⁡(φl-φw)sin(φl-φw)]-[vw0]-[-ωw·ydωw·xd].Additionally, there is(9)φ˙w=ωw.Combing (8) and (9), we can obtain
(10)x˙d=vlcos⁡(φl-φw)-vw+ωw·yd,y˙d=vlsin(φl-φw)-ωw·xd,φ˙w=ωw.Then, (10) can also be written as the state equations which is shown as follows
(11)[x˙dy˙dφ˙w]=[-1yd0-xd01][vwωw]+[vlcos⁡(φl-φw)vlsin(φl-φw)0].The output equation is(12)y1=xd,y2=yd.


## 2.2. Discrete Model of UAV Formation Movement and Its Predictive Control Analysis

In the previous section, the state equation of formation control is obtained. Here, discretize this equation and the following equation can be obtained(13)[xd(k+1)yd(k+1)φw(k+1)]=[xd(k)yd(k)φw(k)]+[-1yd(k)0-xd(k)01][vw(k)ωw(k)]ΔT+[vl(k)cos⁡⁡(φl(k)-φw(k))vl(k)sin⁡(φl(k)-φw(k))0]ΔT.Generally, the sampling periodic time is short during receding optimization process, so the velocity and yaw angle of leader UAV can be considered constant in sampling period [15], which means in a short sampling period, there is
(14)vl(k+i)=vl(k+i-1)=vl(k),φl(k+i)=φl(k+i-1)=φl(k).According to these two equations, the predicted value of outputs will be obtained as follows:(15)xdp(k+1),xdp(k+2),…,xdp(k+N),ydp(k+1),ydp(k+2),…,ydp(k+N).Since at timek, yd(k), vl(k), φw(k), φl(k) are known, the future output values of formation flight control are xdp(k+1),xdp(k+2),…,xdp(k+N-1),xdp(k+N), and ydp(k+1),ydp(k+2),…,ydp(k+N-1),ydp(k+N), and these outputs are only the function of the future control quantities vw(k),vw(k+1),…,vw(k+N-1) and ωw(k),ωw(k+1),…,ωw(k+N-1). Obviously, this function is a multi-input-multioutput nonlinear control problem. So, these values can be obtained using the receding optimization algorithm. Then, using vw(k) and ωw(k) as outputs of control quantity, and carry out the receding optimization algorithm in sequence, the input of control quantity in the next time can be obtained.According to Section1, it can be known that for the problem of UAV formation flight, nonlinear model predictive control and fuzzy stair-like predictive control have some limitations, but the multimodel control method has stronger robustness and higher control accuracy, the final predictive model of which is linear, so the receding optimization problem can be changed from general nonlinear optimization problem to a linear quadratic optimization problem. Since the linear quadratic optimization has faster computation than the ordinary nonlinear optimization, the multiple models can greatly improve real-time of receding optimization [16]. Therefore, multiple model-based predictive control approach is used to design the controller of UAV formation flight.


## 3. Predictive Control for UAV Formation Based on Multimodel Approach

The basic principle of the multimodel control method is that nearby the different characteristics of nonlinear systems, different linear model is used to describe this nonlinear system, and each linear model only describes a part of nonlinear system dynamics. The multiple linearization models are used to approximate the nonlinear system in its entire operating range, and the controller is designed based on each linearization model. These controllers are combined together to constitute a multimodel controller in some way. Finally, the control of entire nonlinear system can be achieved through coordinated control between multiple linearization models. The basic steps of this method can be summarized in four steps:(1) acquisition of the model set; (2) local linearization of the model set; (3) establishment of the controller set; (4) combination of the model set.Similarly, the controller design of UAV formation flight can also include these steps, the flow chart of which is shown in Figure2. The details of controller design will be presented in the following parts using the multimodel prediction method.Figure 2
The flow chart of UAV formation flight controller design.


### 3.1. Determination of the State Characteristic Points of Formation Control Model

According to the basic principles of multimodel predictive control, the characteristic state points of the nonlinear model must be obtained first before obtaining multiple models. Characteristic state points and their region are determined using methods in [15, 17]. As is shown in Figure 3, the horizontal axis and vertical axis denote time and trajectory, respectively.Figure 3
Description of tangent error.The basic idea is as follows.Determine the first characteristic state points, then compute the error between reference trajectory tangent through the characteristic point and reference trajectory and compare the error with the maximum permissible error. If it is greater than the maximum permissible error, redetermine the characteristic point to get the next characteristic state point, and calculate repeatedly until the last characteristic state point is obtained; thus, the characteristic state points of nonlinear system and their applicable region can be determined. So the linearization model set can be obtained through the linearization process at the characteristic state point for each region.For UAV formation control, the characteristic state points are determined as follows.Assuming that the initial distance between two vehicles arexd(0) and yd(0), respectively, when UAV formation control is carried out, so the desired distances of UAV formation are xc and yc. Assuming that the expected time arrival at the desired value ofX-axis andY-axis is the same, the reference trajectory of UAV formation can be obtained as follows:
(16)y1r(t)=ξtxd(0)+(1-ξt)·xc,y2r(t)=ξtyd(0)+(1-ξt)·yc.Since there are two outputs, generating two reference trajectories, the algorithm above cannot be applied directly. But, because the expected arrival time to the desired value is the same, the reference trajectory of an output can be used to determine the characteristic point of state. Assume that the reference trajectory of relative position on theX-axis in the track coordinate system is used to determine characteristic point of the state.At the characteristic state points, there is(17)x˙d=0,y˙d=0,φ˙w=0.Assume thatith characteristic state point is (Vwi,ωwi,xdi,ydi,φdi), then there will be the following equation:
(18)[-1ydi0-xdi01][Vwiωwi]+[Vlcos⁡(φl-φwi)Vlsin(φl-φwi)0]=[000].Meanwhile, there is(19)xdi=y1r(i),ydi=y2r(i).Solve (18) and (19), then we can obtain
(20)ωwi=0,φwi=φl,vwi=vl,xdi=y1r(i),ydi=y2r(i).Assuming the initial point is as the first characteristic state point, then the state point value will be(21)ωw1=0,φw1=φl,vw1=vl,xd1=y1r(0),yd1=y2r(0).Then, the equation of the tangent is(22)y1k(t)=y˙1r(t)t+xd(0).Afterward, determine the maximum permissible errorEmax⁡, and the time t2 corresponding to the second characteristic point can be obtained using the method above. Then, the state point can be obtained as follows
(23)ωw2=0,φw2=φl,vw2=vl,xd2=xd(t2),yd2=yd(t2).So the tangent equation is(24)y1k(t)=y˙1r(t)(t-t2)+xd(t2).By repeating these procedures above, timeti corresponding to the characteristic point in the region of multimodel for formation can be obtained, and at the same time the characteristic state point corresponding to the time ti can be also obtained:
(25)ωwi=0,φwi=φl,vwi=vl,xdi=xd(ti),ydi=yd(ti).Calculate until the last characteristic state point is obtained, and then the computation will be terminated.


### 3.2. Generation of Discrete Model Sets for Formation Control

After obtaining the characteristic state points, carry on linearization at different discrete model sets of formation. Here, linearization can be realized through the following methods.Consider nonlinear systems as described in the form of discrete-time dynamic equations:(26)x(k+1)=f(x(k),u(k)),y(k)=g(x(k),u(k)).The system hasm different characteristic state points. f(x(k),u(k)) and g(x(k),u(k)) have the first continuous partial derivative. If system is linearized at each characteristic state point, the standard discrete state-space model ofm linear models of the original system is obtained as follows:
(27)x(k+1)=Aix(k)+Biu(k)-αi.y(k)=Cix(k)+Diu(k)-βi,
where
(28)Ai=∂f∂x|(xi,ui),Bi=∂f∂u|(xi,ui),Ci=∂g∂x|(xi,ui),Di=∂g∂u|(xi,ui),αi=∂f∂x|(xi,ui)xi+∂f∂u|(xi,ui)ui-xi,αi=∂g∂x|(xi,ui)xi+∂g∂u|(xi,ui)ui-yi.Here,m linearized models constitute the linearized multimodel presentation of the original system.For the UAV formation flying control, the characteristic state points are shown as(29)ωwi=0,φwi=φl,vwi=vl,xdi=xd(ti),ydi=yd(ti).
And it has the following expression:
(30)Ai=∂f∂x|(xi,ui)=[10001-vl(k)ΔT001],Bi=∂f∂u|(xi,ui)=[-1yd(ti)0-xd(ti)01]ΔT,αi=∂f∂x|(xi,ui)xi+∂f∂u|(xi,ui)ui-xi=[-vl(k)ΔT-vl(k)φl(k)ΔT0].Since outputs are linear,Ci, Di and βi will not be solved. Thus, the linearization equation at characteristic state point will be
(31)[xd(k+1)yd(k+1)φd(k+1)]=[10001-vl(k)ΔT001][xd(k)yd(k)φd(k)]+[-1yd(ti)0xd(ti)01]ΔT[vw(k)ωw(k)]-[-vl(k)ΔT-vl(k)φl(k)ΔT0].So for different characteristic state points, the linear model for UAV formation control at different horizons can be obtained, realizing the acquisition of model sets for UAV formation. And these models are denoted asM(1),M(2),…,M(S).


### 3.3. Combined Method of Model Sets for Formation Flight Control

The characteristic state points are obtained using the method above. When the error reaches the maximum value, switch to the new model, which ensures the maximum error value between the predictive trajectory and the reference trajectory. Thus, the determination of model region can be realized. It can be seen that the applicable range of desired model is divided based on the time region, so different sampling points have different models. But the predictive control is based on the future time region. So in this paper, the applicable model of predictive point is judged by the time region, and then this model is used to calculate the predictive value. The judgment rules of predictive model are described as follows.Assuming that the time corresponding to the state characteristic points ist1<t2<⋯<ts-1<ts and the predictive horizon is [t,t+N], then there will be the following.If[t,t+N]∈[ti,ti+1], then the final predictive model of all points is MT=Mi; if [t,t+N]∈[ti,ti+N], one can judge the interval [th,th+1] of all points between prediction point t+1 to t+N in sequence and denote the model of this point as MT=Mh; if [t,t+N]∈[ts,∞), the final predictive model of all the points will be MT=Ms.Based on this method, we can obtain the predictive function during the future horizon to determine the optimization index. However, the boundary point of the predictive range may be closer to the next linear model, as shown in Figure4. Sampling point p1 may be closer to the linear model at sampling point t4. But according to the method above, the calculation model used at the sampling point p1 is the linear model of the state characteristic point at sampling point t3. Sampling point p2 may be closer to the linear model of the state characteristic point at sampling point t5. But according to the method above, the calculation model used by the sampling point p2 is the linear model of the state characteristic point at sampling point t4.Figure 4
The schematic for determining model of frontier points during the predictive intervals.So the method above may decrease the performance of the approximation capability on the boundary, and each model belonging to the model set cannot switch smoothly [18]. However, T-S fuzzy model as an intelligent control method, mainly uses fuzzy reasoning to approximate the nonlinear system. Using this method, the input space can be divided into several fuzzy subspaces, where a local linear model is established, and then the local models are combined smoothly using the membership function, forming a global fuzzy model of nonlinear function, which is ultimately identified as a linear model [19]. The predictive control method based on the T-S fuzzy model belongs to multimodel predictive control with the weighted models. Compared with the common multimodel predictive controllers with weighted models, the fuzzy weighted models have more accurate nonlinear approximation performance; switch of the model is more smooth, and it is easier to understand [20, 21]. So in this section, T-S fuzzy idea is adopted for the multimodel control of UAV formation flight, as is shown in Figure 5.Figure 5
The schematic of multimodel control method for UAV formation flight.For each sampling point, use the error between the tangent of state characteristic point and the reference trajectory of this sampling point to calculate the membership degree. Assuming that the error between tangent ofjth state characteristic point and the reference trajectory at the sampling pointt is Ei(t), so for the point t+i in the predictive range, weighted function is as follows:
(32)wj(t+i)=e-(Ei(t+i)/Emax⁡)2,-Emax⁡≤Ej(t+i)≤Emax⁡,w(t+i)=0,Ej(t+i)≤-Emax⁡orEj(t+i)≥Emax⁡.This equation can ensure that the farther away from the state characteristic point the sampling point is, the lower its weighted value is. Using the weighted values, the prediction model at the sampling pointt is
(33)Mt=∑i=1m(wiMi)∑i=1m(wi).For the sampling points during the predictive horizon, there is(34)yp(t+i)=∑j=1m(wj(t+i)Mj)∑j=1m(wj(t+i)),
where
(35)wj(t+i)=e-(Ei(t+i)/Emax⁡)2,-Emax⁡≤Ej(t+i)≤Emax⁡,w(t+i)=0,Ej(t+i)≤-Emax⁡orEj(t+i)≥Emax⁡.Through this approach, the linear prediction function for UAV formation can be obtained as follows:(36)xdp(k+1),xdp(k+2),…,xdp(k+N),ydp(k+1),ydp(k+2),…,ydp(k+N).In this way, predictive outputs can change from a nonlinear function to a linear function. This nonlinear function includesvw(k),vw(k+1),…,vw(k+N-1), and ωw(k),ωw(k+1),…,ωw(k+N-1), while this linear function includes the control quantities mentioned above. Thus, the control problem will become a multi-input-multioutput linear predictive control problem.


### 3.4. Optimization Index for Formation and Receding Optimization Solution

During the predictive control process, the goal of receding optimization is to find a set ofvw(k),vw(k+1),…,vw(k+N-1) and ωw(k),ωw(k+1),…,ωw(k+N-1), making prediction outputs at entire optimization horizon as close to the reference trajectory as possible.Here, introduce the closed-loop(37)e1(t)=xd(t)-xdp(t),e2(t)=yd(t)-ydp(t).The open-loop predictive output can be directly compensated by the output feedback, and then the predictive value of the closed-loop model will be(38)xd(t+i)=xdp(t+i)+e1(t),yd(t+i)=ydp(t+i)+e2(t).In this section, there are two control objectives: the relative distances toX-axis andY-axis. Since they have equal importance and the same unit of quantity, they are set with the same weight when designing the performance index. Thus, the performance index is defined as follows:
(39)J=∑i=1N[(y1r(t+i)-xd(t+i))2+(y2r(t+i)-yd(t+i))2]+λ1∑i=0N-1(vw(t+i)-vw(t))2+λ2∑i=0N-1ωw2(t+i).Similarly, the optimization constraints of control quantity are introduced as follows:(40)vw(k+i-1)-Δv<vw(k+i)<vw(k+i-1)+Δv,ωmin⁡<ωw(k+i)<ωmax⁡,vmin⁡<vw(k+i)<vmax⁡,wherei∈0,1,2,…,N-1.After using multiple models, the performance index is linear quadratic, whose constraints are linear equality and inequality, so the optimization problem is a linear quadratic programming problem. The solution methods of linear quadratic programming problem can be used to solve the receding optimization problem. The linear quadratic programming problem is a common programming problem and has a lot of solution methods and higher speed than the ordinary nonlinear programming, which increases the speed of receding optimization solution [22].


## 3.1. Determination of the State Characteristic Points of Formation Control Model

According to the basic principles of multimodel predictive control, the characteristic state points of the nonlinear model must be obtained first before obtaining multiple models. Characteristic state points and their region are determined using methods in [15, 17]. As is shown in Figure 3, the horizontal axis and vertical axis denote time and trajectory, respectively.Figure 3
Description of tangent error.The basic idea is as follows.Determine the first characteristic state points, then compute the error between reference trajectory tangent through the characteristic point and reference trajectory and compare the error with the maximum permissible error. If it is greater than the maximum permissible error, redetermine the characteristic point to get the next characteristic state point, and calculate repeatedly until the last characteristic state point is obtained; thus, the characteristic state points of nonlinear system and their applicable region can be determined. So the linearization model set can be obtained through the linearization process at the characteristic state point for each region.For UAV formation control, the characteristic state points are determined as follows.Assuming that the initial distance between two vehicles arexd(0) and yd(0), respectively, when UAV formation control is carried out, so the desired distances of UAV formation are xc and yc. Assuming that the expected time arrival at the desired value ofX-axis andY-axis is the same, the reference trajectory of UAV formation can be obtained as follows:
(16)y1r(t)=ξtxd(0)+(1-ξt)·xc,y2r(t)=ξtyd(0)+(1-ξt)·yc.Since there are two outputs, generating two reference trajectories, the algorithm above cannot be applied directly. But, because the expected arrival time to the desired value is the same, the reference trajectory of an output can be used to determine the characteristic point of state. Assume that the reference trajectory of relative position on theX-axis in the track coordinate system is used to determine characteristic point of the state.At the characteristic state points, there is(17)x˙d=0,y˙d=0,φ˙w=0.Assume thatith characteristic state point is (Vwi,ωwi,xdi,ydi,φdi), then there will be the following equation:
(18)[-1ydi0-xdi01][Vwiωwi]+[Vlcos⁡(φl-φwi)Vlsin(φl-φwi)0]=[000].Meanwhile, there is(19)xdi=y1r(i),ydi=y2r(i).Solve (18) and (19), then we can obtain
(20)ωwi=0,φwi=φl,vwi=vl,xdi=y1r(i),ydi=y2r(i).Assuming the initial point is as the first characteristic state point, then the state point value will be(21)ωw1=0,φw1=φl,vw1=vl,xd1=y1r(0),yd1=y2r(0).Then, the equation of the tangent is(22)y1k(t)=y˙1r(t)t+xd(0).Afterward, determine the maximum permissible errorEmax⁡, and the time t2 corresponding to the second characteristic point can be obtained using the method above. Then, the state point can be obtained as follows
(23)ωw2=0,φw2=φl,vw2=vl,xd2=xd(t2),yd2=yd(t2).So the tangent equation is(24)y1k(t)=y˙1r(t)(t-t2)+xd(t2).By repeating these procedures above, timeti corresponding to the characteristic point in the region of multimodel for formation can be obtained, and at the same time the characteristic state point corresponding to the time ti can be also obtained:
(25)ωwi=0,φwi=φl,vwi=vl,xdi=xd(ti),ydi=yd(ti).Calculate until the last characteristic state point is obtained, and then the computation will be terminated.


## 3.2. Generation of Discrete Model Sets for Formation Control

After obtaining the characteristic state points, carry on linearization at different discrete model sets of formation. Here, linearization can be realized through the following methods.Consider nonlinear systems as described in the form of discrete-time dynamic equations:(26)x(k+1)=f(x(k),u(k)),y(k)=g(x(k),u(k)).The system hasm different characteristic state points. f(x(k),u(k)) and g(x(k),u(k)) have the first continuous partial derivative. If system is linearized at each characteristic state point, the standard discrete state-space model ofm linear models of the original system is obtained as follows:
(27)x(k+1)=Aix(k)+Biu(k)-αi.y(k)=Cix(k)+Diu(k)-βi,
where
(28)Ai=∂f∂x|(xi,ui),Bi=∂f∂u|(xi,ui),Ci=∂g∂x|(xi,ui),Di=∂g∂u|(xi,ui),αi=∂f∂x|(xi,ui)xi+∂f∂u|(xi,ui)ui-xi,αi=∂g∂x|(xi,ui)xi+∂g∂u|(xi,ui)ui-yi.Here,m linearized models constitute the linearized multimodel presentation of the original system.For the UAV formation flying control, the characteristic state points are shown as(29)ωwi=0,φwi=φl,vwi=vl,xdi=xd(ti),ydi=yd(ti).
And it has the following expression:
(30)Ai=∂f∂x|(xi,ui)=[10001-vl(k)ΔT001],Bi=∂f∂u|(xi,ui)=[-1yd(ti)0-xd(ti)01]ΔT,αi=∂f∂x|(xi,ui)xi+∂f∂u|(xi,ui)ui-xi=[-vl(k)ΔT-vl(k)φl(k)ΔT0].Since outputs are linear,Ci, Di and βi will not be solved. Thus, the linearization equation at characteristic state point will be
(31)[xd(k+1)yd(k+1)φd(k+1)]=[10001-vl(k)ΔT001][xd(k)yd(k)φd(k)]+[-1yd(ti)0xd(ti)01]ΔT[vw(k)ωw(k)]-[-vl(k)ΔT-vl(k)φl(k)ΔT0].So for different characteristic state points, the linear model for UAV formation control at different horizons can be obtained, realizing the acquisition of model sets for UAV formation. And these models are denoted asM(1),M(2),…,M(S).


## 3.3. Combined Method of Model Sets for Formation Flight Control

The characteristic state points are obtained using the method above. When the error reaches the maximum value, switch to the new model, which ensures the maximum error value between the predictive trajectory and the reference trajectory. Thus, the determination of model region can be realized. It can be seen that the applicable range of desired model is divided based on the time region, so different sampling points have different models. But the predictive control is based on the future time region. So in this paper, the applicable model of predictive point is judged by the time region, and then this model is used to calculate the predictive value. The judgment rules of predictive model are described as follows.Assuming that the time corresponding to the state characteristic points ist1<t2<⋯<ts-1<ts and the predictive horizon is [t,t+N], then there will be the following.If[t,t+N]∈[ti,ti+1], then the final predictive model of all points is MT=Mi; if [t,t+N]∈[ti,ti+N], one can judge the interval [th,th+1] of all points between prediction point t+1 to t+N in sequence and denote the model of this point as MT=Mh; if [t,t+N]∈[ts,∞), the final predictive model of all the points will be MT=Ms.Based on this method, we can obtain the predictive function during the future horizon to determine the optimization index. However, the boundary point of the predictive range may be closer to the next linear model, as shown in Figure4. Sampling point p1 may be closer to the linear model at sampling point t4. But according to the method above, the calculation model used at the sampling point p1 is the linear model of the state characteristic point at sampling point t3. Sampling point p2 may be closer to the linear model of the state characteristic point at sampling point t5. But according to the method above, the calculation model used by the sampling point p2 is the linear model of the state characteristic point at sampling point t4.Figure 4
The schematic for determining model of frontier points during the predictive intervals.So the method above may decrease the performance of the approximation capability on the boundary, and each model belonging to the model set cannot switch smoothly [18]. However, T-S fuzzy model as an intelligent control method, mainly uses fuzzy reasoning to approximate the nonlinear system. Using this method, the input space can be divided into several fuzzy subspaces, where a local linear model is established, and then the local models are combined smoothly using the membership function, forming a global fuzzy model of nonlinear function, which is ultimately identified as a linear model [19]. The predictive control method based on the T-S fuzzy model belongs to multimodel predictive control with the weighted models. Compared with the common multimodel predictive controllers with weighted models, the fuzzy weighted models have more accurate nonlinear approximation performance; switch of the model is more smooth, and it is easier to understand [20, 21]. So in this section, T-S fuzzy idea is adopted for the multimodel control of UAV formation flight, as is shown in Figure 5.Figure 5
The schematic of multimodel control method for UAV formation flight.For each sampling point, use the error between the tangent of state characteristic point and the reference trajectory of this sampling point to calculate the membership degree. Assuming that the error between tangent ofjth state characteristic point and the reference trajectory at the sampling pointt is Ei(t), so for the point t+i in the predictive range, weighted function is as follows:
(32)wj(t+i)=e-(Ei(t+i)/Emax⁡)2,-Emax⁡≤Ej(t+i)≤Emax⁡,w(t+i)=0,Ej(t+i)≤-Emax⁡orEj(t+i)≥Emax⁡.This equation can ensure that the farther away from the state characteristic point the sampling point is, the lower its weighted value is. Using the weighted values, the prediction model at the sampling pointt is
(33)Mt=∑i=1m(wiMi)∑i=1m(wi).For the sampling points during the predictive horizon, there is(34)yp(t+i)=∑j=1m(wj(t+i)Mj)∑j=1m(wj(t+i)),
where
(35)wj(t+i)=e-(Ei(t+i)/Emax⁡)2,-Emax⁡≤Ej(t+i)≤Emax⁡,w(t+i)=0,Ej(t+i)≤-Emax⁡orEj(t+i)≥Emax⁡.Through this approach, the linear prediction function for UAV formation can be obtained as follows:(36)xdp(k+1),xdp(k+2),…,xdp(k+N),ydp(k+1),ydp(k+2),…,ydp(k+N).In this way, predictive outputs can change from a nonlinear function to a linear function. This nonlinear function includesvw(k),vw(k+1),…,vw(k+N-1), and ωw(k),ωw(k+1),…,ωw(k+N-1), while this linear function includes the control quantities mentioned above. Thus, the control problem will become a multi-input-multioutput linear predictive control problem.


## 3.4. Optimization Index for Formation and Receding Optimization Solution

During the predictive control process, the goal of receding optimization is to find a set ofvw(k),vw(k+1),…,vw(k+N-1) and ωw(k),ωw(k+1),…,ωw(k+N-1), making prediction outputs at entire optimization horizon as close to the reference trajectory as possible.Here, introduce the closed-loop(37)e1(t)=xd(t)-xdp(t),e2(t)=yd(t)-ydp(t).The open-loop predictive output can be directly compensated by the output feedback, and then the predictive value of the closed-loop model will be(38)xd(t+i)=xdp(t+i)+e1(t),yd(t+i)=ydp(t+i)+e2(t).In this section, there are two control objectives: the relative distances toX-axis andY-axis. Since they have equal importance and the same unit of quantity, they are set with the same weight when designing the performance index. Thus, the performance index is defined as follows:
(39)J=∑i=1N[(y1r(t+i)-xd(t+i))2+(y2r(t+i)-yd(t+i))2]+λ1∑i=0N-1(vw(t+i)-vw(t))2+λ2∑i=0N-1ωw2(t+i).Similarly, the optimization constraints of control quantity are introduced as follows:(40)vw(k+i-1)-Δv<vw(k+i)<vw(k+i-1)+Δv,ωmin⁡<ωw(k+i)<ωmax⁡,vmin⁡<vw(k+i)<vmax⁡,wherei∈0,1,2,…,N-1.After using multiple models, the performance index is linear quadratic, whose constraints are linear equality and inequality, so the optimization problem is a linear quadratic programming problem. The solution methods of linear quadratic programming problem can be used to solve the receding optimization problem. The linear quadratic programming problem is a common programming problem and has a lot of solution methods and higher speed than the ordinary nonlinear programming, which increases the speed of receding optimization solution [22].


## 4. Simulation

In this section, numerical simulations are performed to demonstrate the performance of the proposed approach. Here, the formation control ability can be tested in two important scenarios. Simulation scenarios are set as follows. One scenario is the leader UAV flying straight, and the other is the leader UAV flying with turning course. Additionally, the comparison simulation between single MPC (SMPC) method and multiple MPC (MMPC) method is carried on to verify effectiveness of the method in this paper. Meanwhile, the parameters used in the simulations are set as follows. The prediction horizonN is 5 and the sampling interval is 0.2 s. The angular velocity and velocity of two vehicles are confined during the interval (-0.1,0.1) and the interval (35,45), respectively. All the computations and experiments have been on a computer with Inter Core i3 CPU, 3.30 GHz, and Windows XP operating systems. Table 1 summarizes the initial conditions of the formation.Table 1
Initial conditions of UAV formation.
Initial conditions
The role of UAV
Leader
Follower
Initial position
(
0,0
)
(
-
100
,
-
100
)
Initial angle
0
π
/
2
Initial velocity
40
40
Initial angular velocity
0
0


### 4.1. Formation Simulation of Leader UAV Flying Straight

The simulation experiment is mainly used to verify the UAV formation control capability when the leader UAV is flying straight. Here, error exists in the position measurement of leader UAV, which is ±0.5 m. There are two different control goals. One is that the relative position between follower and leader of UAV formation in the track coordinates is as follows:(41)X=-60,Y=30.The other is that the formation should be formed within 40 s.Because the leader UAV has its initial angle ofπ/4 and it flies straight, the initial relative position in the track coordinates will be obtained as follows:
(42)[xdyd]=[cos⁡π4sinπ4-sinπ4cos⁡π4][-100-100]=[-10020].Simulation is carried out by using Matlab Simulink toolbox, and the simulation results are shown from Figures6, 7, 8, 9, and 10.Figure 6
Flight trajectories of two UAVs.Figure 7
Relative position ofX-axis.Figure 8
Relative position ofY-axis.Figure 9
Change curve of follower UAV’s velocity.Figure 10
Change curve of follower UAV’s yaw rate.According to Figures6–10, it can be seen that when leader UAV is navigating in a straight line, formation control can be achieved through both SMPC and MMPC method. However, the SMPC method has a larger tracking error than the MMPC method. Meanwhile, it can also be seen that it takes a longer time for SMPC method than MMPC method to form a steady formation. The UAV formation can be realized in 40 seconds by the MMPC method, which meets the actual design demand.


### 4.2. Formation Simulation for Leader UAV with Turning Flight

The UAV formation control capability is proved in this section when the leader UAV flies with a turning flight path. The UAV flies 20 s with an initial angle of 0° between the leader UAV andX-axis, and then the UAV flies with angular velocity of π/200 for 100 seconds, and then it moves straight inY-axis direction. There are also two different control goals. One is that the relative position between follower and leader of UAV formation in the track coordinates is as follows:
(43)[xdrefydref]=[-50-50].The other is the formation should form within 40 s.From Table1, the relative position in the track coordinate system between two vehicles is obtained as follows:
(44)[xdyd]=[cos⁡0sin0-sin0cos⁡0][-100-100]=[-100-100].Simulation is carried out by using Matlab Simulink toolbox, and the results are shown in Figures11, 12, 13, 14, and 15.Figure 11
Flight trajectories of two UAVs.Figure 12
Relative position ofX-axis.Figure 13
Relative position ofY-axis.Figure 14
Change curve of follower UAV’s velocity.Figure 15
Change curve of follower UAV’s yaw rate.According to Figures10–14, when leader UAV flies with a turning flight path, using method proposed in the paper, formation control can be achieved better than the SMPC method, whenever the UAV flies straight or flies with a turning path. The SMPC method has a larger tracking error than the MMPC method. Meanwhile, it can also be seen that it takes a longer time for SMPC method than MMPC method to form a steady formation. The UAV formation can be realized in 40 seconds by the MMPC method, which meets the actual design demand.According to the Matlab simulation process of UAV formation in those two scenarios above, when the sampling interval is 0.2 s, the simulation time of the receding optimization program on the PC is less than 0.2 s each time and the time will be shorter if the simulation is done on a dedicated chip. So it meets the real-time needs. It can be seen from the relative position on theX-axis andY-axis of two vehicles in the track coordinate system that the UAV formation is realized within 40 s. All in all, the simulation shows that the control requirements and real-time requirements can be satisfied by using multimodel predictive control method for UAV formation control.


## 4.1. Formation Simulation of Leader UAV Flying Straight

The simulation experiment is mainly used to verify the UAV formation control capability when the leader UAV is flying straight. Here, error exists in the position measurement of leader UAV, which is ±0.5 m. There are two different control goals. One is that the relative position between follower and leader of UAV formation in the track coordinates is as follows:(41)X=-60,Y=30.The other is that the formation should be formed within 40 s.Because the leader UAV has its initial angle ofπ/4 and it flies straight, the initial relative position in the track coordinates will be obtained as follows:
(42)[xdyd]=[cos⁡π4sinπ4-sinπ4cos⁡π4][-100-100]=[-10020].Simulation is carried out by using Matlab Simulink toolbox, and the simulation results are shown from Figures6, 7, 8, 9, and 10.Figure 6
Flight trajectories of two UAVs.Figure 7
Relative position ofX-axis.Figure 8
Relative position ofY-axis.Figure 9
Change curve of follower UAV’s velocity.Figure 10
Change curve of follower UAV’s yaw rate.According to Figures6–10, it can be seen that when leader UAV is navigating in a straight line, formation control can be achieved through both SMPC and MMPC method. However, the SMPC method has a larger tracking error than the MMPC method. Meanwhile, it can also be seen that it takes a longer time for SMPC method than MMPC method to form a steady formation. The UAV formation can be realized in 40 seconds by the MMPC method, which meets the actual design demand.


## 4.2. Formation Simulation for Leader UAV with Turning Flight

The UAV formation control capability is proved in this section when the leader UAV flies with a turning flight path. The UAV flies 20 s with an initial angle of 0° between the leader UAV andX-axis, and then the UAV flies with angular velocity of π/200 for 100 seconds, and then it moves straight inY-axis direction. There are also two different control goals. One is that the relative position between follower and leader of UAV formation in the track coordinates is as follows:
(43)[xdrefydref]=[-50-50].The other is the formation should form within 40 s.From Table1, the relative position in the track coordinate system between two vehicles is obtained as follows:
(44)[xdyd]=[cos⁡0sin0-sin0cos⁡0][-100-100]=[-100-100].Simulation is carried out by using Matlab Simulink toolbox, and the results are shown in Figures11, 12, 13, 14, and 15.Figure 11
Flight trajectories of two UAVs.Figure 12
Relative position ofX-axis.Figure 13
Relative position ofY-axis.Figure 14
Change curve of follower UAV’s velocity.Figure 15
Change curve of follower UAV’s yaw rate.According to Figures10–14, when leader UAV flies with a turning flight path, using method proposed in the paper, formation control can be achieved better than the SMPC method, whenever the UAV flies straight or flies with a turning path. The SMPC method has a larger tracking error than the MMPC method. Meanwhile, it can also be seen that it takes a longer time for SMPC method than MMPC method to form a steady formation. The UAV formation can be realized in 40 seconds by the MMPC method, which meets the actual design demand.According to the Matlab simulation process of UAV formation in those two scenarios above, when the sampling interval is 0.2 s, the simulation time of the receding optimization program on the PC is less than 0.2 s each time and the time will be shorter if the simulation is done on a dedicated chip. So it meets the real-time needs. It can be seen from the relative position on theX-axis andY-axis of two vehicles in the track coordinate system that the UAV formation is realized within 40 s. All in all, the simulation shows that the control requirements and real-time requirements can be satisfied by using multimodel predictive control method for UAV formation control.


## 5. Conclusion

In this paper, the main work can be concluded as follows to solve the problem of UAV formation control.(1)
Discrete relative motion equations are established for UAV formation by using the leader-follower method.(2)
Multimodel sets for UAV formation are established, and the weighted model sets method is proposed.(3)
The formation controller based on multimodel predictive control is designed.(4)
Simulation in two scenarios is carried out and the effectiveness of controller designed and control strategy is verified.The multimodel predictive control method can be used for UAV formation control. This method can meet control requirements and real-time requirements well. The result of this paper is the basis of further research on formation reconfiguration control problem. In the future, we will introduce the approach proposed in this paper to the controller design of actual UAV formation flight.



---

*Source: 835301-2014-05-06.xml*