# Recurrent Neural Network-Based Model Predictive Control for Multiple Unmanned Quadrotor Formation Flight

Authors: Boyang Zhang; Xiuxia Sun; Shuguang Liu; Xiongfeng Deng
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

This paper presents a dynamical recurrent neural network- (RNN-) based model predictive control (MPC) structure for the formation flight of multiple unmanned quadrotors. A distributed hierarchical control system with the translation subsystem and rotational subsystem is proposed to handle the formation-tracking problem for each quadrotor. The RNN-based MPC is proposed for each subsystem, where the RNN is introduced as the predictive model in MPC. And to improve the modeling accuracy, an adaptive updating law is developed to tune weights online for the RNN. Besides, the adaptive differential evolution (DE) algorithm is utilized to solve the optimization problem for MPC. Furthermore, the closed-loop stability is analyzed; meanwhile, the convergence of the DE algorithm is discussed as well. Finally, some simulation examples are provided to illustrate the validity of the proposed control structure.


---

## Body

## 1. Introduction

The unmanned quadrotor is a class of pilot-less aerial vehicles which enable to hover, take off, and land vertically with aerodynamic and propulsion characteristics. Due to the capabilities of simplicity, maneuverability, and payload [1, 2], the unmanned quadrotor has been gaining extensive attention and plays an important role in the military and civilian fields, such as real-time reconnaissance, surveillance, search and rescue missions, bush fire monitoring, agricultural crop dusting, and different airborne operations [3–5]. These flight missions not only impose stringent requirements over the modeling and control of quadrotor dynamics but also increasingly require multiple quadrotors to coordinate in a formation. Here, the basis of these flight missions is the trajectory tracking. However, the dynamics of the quadrotor are highly nonlinear and subject to uncertain disturbances, so it is a challenge to design a desirable control system for the unmanned quadrotor.To solve the above problems, several control approaches have been investigated in recent years, such as proportional-integral-derivative (PID) control and linear quadratic regulator (LQR) control [6], sliding mode control [7], and backstepping control [8]. However, these methods only consider the control effects and often neglect the actual flight constraints. Besides, the optimality of the path cannot be ensured, as the selection of reactive parameters is aimless or blind.To solve the problems further, a new developing method, the model predictive control (MPC), has been proposed for the unmanned quadrotor system [9, 10], where the physical constraints of the system inputs and outputs could be easily solved by using the proper penalty terms in cost function. In [11], the MPC method was proposed to solve the path-tracking problem for the unmanned quadrotor, where a future control sequence was computed in a defined horizon, such that the prediction of the quadrotor output was driven close to the reference. In [12], a novel control strategy which combined the active disturbance rejection control and model predictive control was presented for the unmanned quadrotor. The proposed method improved the robustness for the modeling error and disturbances while performing a smooth trajectory tracking.The performance of MPC highly depends on the accuracy of the described predictive model. Conventional techniques of system identification such as the maximum likelihood estimation method, modified maximum likelihood estimation method, and Kalman filtering method [13–15] have been widely used in system identification. However, these methods require the known structure of the mathematical model, which is difficult for the unmanned quadrotor to know its accurate mathematical model, considering uncertain disturbances in the complex environment.As a nonconventional method, the neural network (NN) does not explicitly need an accurate mathematical model [16–19]. Based on that, some works on NN-based MPC have improved modeling performance for many different nonlinear systems. In [20], Nikdel et al. presented a model predictive controller with the feedforward neural network (FNN) model for a single degree of freedom-rotary manipulator achieving well closed-loop performance. In [21], a model predictive controller based on the FNN model was utilized to improve the transient stability of a power system. However, the main drawback of FNNs considered in [20, 21] is that they essentially own static input-to-output maps and their capability of representing complex and time-varying nonlinear systems is limited. Besides, they can only provide predetermined finite predictions [22, 23]. To overcome the above drawbacks, a recurrent neural network (RNN) is introduced for MPC, which can capture the system behavior dynamically and provide long-range predictions even in the presence of disturbances [24–26]. In a RNN, the current control signal and state are given as the inputs to obtain the outputs which are fed back as new states. Therefore, RNN is more suitable to model the quadrotor system for MPC in this paper.On the other hand, the performance of MPC also relies on the local optimization problem in MPC, which can be solved by either analytical or numerical solutions. The solution depends on the structure of the cost function. For the tracking problem of the multiquadrotor formation in this paper, the corresponding cost function is complex and nondifferentiable. Traditional search techniques using the characteristics of the problem to determine the next sampling point (gradients, Hessians, and linearity) are not feasible, because the characteristics of tracking performance are affected by mismatched disturbances. Nevertheless, stochastic search techniques, such as the genetic algorithm (GA), particle swarm optimization (PSO) algorithm, and chaotic predator prey biogeography-based optimization (CPPBBO) algorithm [27–29], make no such solutions. Instead, the next point is determined by the stochastic decision rules. Therefore, this paper proposes a stochastic search technique based on the differential evolution (DE) algorithm [30] to solve the optimization problem in MPC, which has been successfully applied in many fields recently [31, 32].Although the concept of NN-based MPC is not new [33, 34], the application of RNN-based MPC in multiquadrotor formation flight is rather scarce. Motivated by this fact, this paper presents the dynamical RNN-based MPC strategy for the multiquadrotor formation flight. The main contributions are summarized as follows:
(i)
A distributed hierarchical control system with the translation subsystem and rotational subsystem is proposed to simplify the control process for the unmanned quadrotor
(ii)
To overcome the drawbacks of conventional MPC, RNN is introduced as the predictive model in MPC to deal with the mismatched modeling for both subsystems
(iii)
The composite control structure based on RNN and MPC is novel for multiquadrotor formation flight, which shows a good control performance in simulation results as well as in theoretical analysisThe remainder of this paper is organized as follows. In Section2, the problem formulation is introduced briefly. Section 3 models the MPC-based tracking problem. In Section 4, the RNN-based prediction model is presented. The control structure is provided in Section 5. In Section 6, the main results concerning stability of the proposed control scheme are analyzed. Simulation results are shown in Section 7. Section 8 includes the conclusion.


## 2. Problem Formulation

Suppose that there areN (N>1) unmanned quadrotors with the same dynamic characteristics moving in R3, which together compose a multiquadrotor formation system. In what follows, the unmanned quadrotor is a miniature four-rotor quadrotor [35]. Each quadrotor can be viewed as a 6 Degree-of-Freedom (DOF) rigid body with generalized coordinates xi=xi,yi,zi,ϕi,θi,ψiT∈R6, where xi,yi,ziT denotes the ith quadrotor’s positions with reference to the earth-fixed frame (E-frame) and ϕi,θi,ψiT represents the attitudes with reference to the body-fixed frame (B-Frame) (i=1,⋯,N). Furthermore, the flight motion is controlled by a total thrust Ui1 and three control torques Ui2, Ui3, and Ui4. And the dynamical model of the ith quadrotor can be defined by the nonlinear state-space function x¨i=fxixi,ẋi,Ui1,Ui2,Ui3,Ui4. To simplify the control process, the quadrotor system can be decomposed into two subsystems in Figure 1: the translational subsystem and the rotational subsystem [36].Figure 1
The distributed control scheme for each quadrotor.Define the system state of the translational subsystem asxiζ=xi,ẋi,yi,ẏi,zi,żiT∈R6 and its control input as uiζ=Ui1,uix,uiyT∈R3 [37], where uix and uiy are two intermediate control inputs, given by
(1)uix=cosϕisinθicosψi+sinϕisinψi,uiy=cosϕisinθisinψi−sinϕicosψi.Then, the dynamics of translational subsystem can be described as(2)ẋiζ=fiζxiζ,uiζ+diζ,yiζ=Cζxiζ,where fiζ: R6×R3→R6 is a nonlinear function, diζ∈R6 is the bounded disturbance, yiζ=xi,yi,ziT is the translational subsystem output, and
(3)Cζ=1,0,0,0,0,00,0,1,0,0,00,0,0,0,1,0.Similarly, the system state and control input of rotational subsystem are represented byxiη=ϕi,ϕ̇i,θi,θ̇i,ψi,ψ̇iT∈R6 and uiη=Ui2,Ui3,Ui4T∈R3, respectively [37]. Then, the dynamic equation is described as
(4)ẋiη=fiηxiη,uiη+diη,yiη=Cηxiη,where fiη: R6×R3→R6 is a nonlinear function, diη∈R6 is the bounded disturbance, yiη=ϕi,θi,ψiT is the rotational subsystem output, and
(5)Cη=1,0,0,0,0,00,0,1,0,0,00,0,0,0,1,0.Based on (2) and (4), both subsystems can be expressed in the discrete-time domain as follows:
(6)xiζk+1=fiζxiζk,uiζk,T+diζk,yiζk=Cζxiζk,xiηk+1=fiηxiηk,uiηk,T+diηk,yiηk=Cηxiηk,where fiζ and fiη are discrete nonlinear functions and T is the sampling time.As shown in Figure2, the leader-follower structure is introduced for the multiquadrotor system. Define a quadrotor as the leader in the formation, and its state vector is represented by xL. Each quadrotor is equipped with a complete set of communication equipment and controllers. Assuming that the information exchange topology is undirected and connected, each follower receives the control commands from the leader directly or indirectly, then executes the corresponding maneuvers to keep the relative state to the leader, where Δxi,L=xi−xL denotes the relative state of the ith follower.Figure 2
Multiquadrotor system.The distance between theith quadrotor and the rth quadrotor can be described as follows:
(7)dirt=xit−xrt2+yit−yrt2+zit−zrt2.In order to avoid collision among quadrotors,dirt should be greater than the safe anti-collision distance Dsafe, i.e.,
(8)dirt≥Dsafe,∀t∈0,T,∀i≠ri,r∈1,⋯,N.Also, to ensure the real-time direct communication of adjacent quadrotors,dirt would be less than the communication distance Dcom, i.e.,
(9)dirt≤Dcom,∀t∈0,T,∃i≠ri,r∈1,⋯,N.


## 3. Model Predictive Control

Generally, MPC is a control process based on predicting future outputs and obtaining nonlinear model outputs at each time horizon [38]. In this paper, a distributed MPC feedback system is designed for each quadrotor in Figure 3. It is mainly composed of the predictive model, reference planning, feedback compensation, and rolling optimization, where uik, yik, rik, yimk, yipk, and dik denote the control input, output state, reference, predictive output, revised predictive output, and disturbance input, respectively.Figure 3
Flowchart of typical MPC.At each time horizon, the predictive outputyimk is obtained by the predictive model, to determine the revised predictive output yipk by feedback compensation. Then, the control input uik is optimally determined by comparing the revised predictive output yipk and the reference rik. Based on that, the system output yik is obtained. Meanwhile, the next time horizon is starting and the new predictive output yimk+1 is predicted by the predictive model. The whole process will be online iterative until the output of system reaches the reference ri.Due to the nonlinearity and time-variant characteristic, the predictive outputyimk is hard to match the real output yik. The predictive errors between the actual output and predictive output are used to revise the predictive output. The formulation can be described as
(10)yipk=yimk+Eyik−1−yimk−1,where E=I3×3 is the weight coefficient matrix.Considering that the quadrotor system decomposed into two subsystems, MPC is introduced separately for both subsystems. In the outer loop, the translational subsystem is controlled to follow a sequence of referencexri,yri,zriT. The optimal control inputs Ui1, uix, and uiy are obtained by minimizing the trajectory tracking errors. With ψri=π/6, the desired attitudes θri and ϕri are calculated using (1). Then, the rotational subsystem’s attitudes ϕi,θi,ψiT are adjusted to achieve the reference ϕri,θri,ψriT in the inner loop. By minimizing attitude errors, the optimal controls Ui2, Ui3, and Ui4 are calculated. Finally, the optimal control inputs Ui1, Ui2, Ui3, and Ui4 are applied to the ith quadrotor.Consider the MPC formulation based on the constrained finite-horizon optimization of tracking problem for translational subsystem [39]:
(11a)minJiζk=∑j=1Hpyipζk+j−riζk+jQζ2+ρ∑j=1HuΔuiζk+j−1Rζ2,subject to
(11b)yminζ≤yipζk+j≤ymaxζ,(11c)uminζ≤uiζk+j−1≤umaxζ,(11d)Δuminζ≤Δuiζk+j−1≤Δumaxζ,where yipζk+j is the revised predictive output and uiζk+j−1 and Δuiζk+j−1 are control inputs and incremental control moves, respectively. riζk+j=xrik+j,yrik+j,zrik+jT is the reference position. yminζ≥ymaxζ and umaxζ≥uminζ are the upper and lower bounds of the system states and control inputs, respectively. ⋅Qζ2 and ⋅Rζ2 are two 2-norms weighted by positive definite matrices Qζ and Rζ, respectively. Hp is the prediction horizon, Hu is the control prediction horizon (Hp≥Hu), and ρ represents the penalty factor.Remark 1.
The tracking problem of the quadrotor formation is decomposed into two phases. Firstly, the leader tracks a given trajectory in a certain state. Then, the followers regard the leader as the trace point to track and further achieve keeping formation. Therefore, if theith quadrotor is the leader, riζk+j denotes reference trajectory; otherwise, it represents the position of the leader for the followers.Similarly, for the rotational subsystem, the MPC optimization problem is described as(12a)minJiηk=∑j=1Hpyipηk+j−riηk+jQζ2+ρ∑j=1HuΔuiηk+j−1Rζ2,subject to
(12b)yminη≤yipηk+j≤ymaxη,(12c)uminη≤uiηk+j−1≤umaxη,(12d)Δuminη≤Δuiηk+j−1≤Δumaxη,where riηk+j=ϕrik+j,θrik+j,ψrik+jT represents the reference attitude.In addition, the nonlinear constraints of (11a) and (12a) should be considered further. A feasible approach is to introduce a penalty cost, transforming the constrained problem to an alternative unconstrained form [39]. Taking the translational subsystem for example, we substitute l1ζk+j=ymaxζ−yipζk+j and l2ζk+j=yipζk+j−yminζ (inequality constraint (11b)), m1ζk+j−1=umaxζ−uiζk+j−1 and m2ζk+j−1=uiζk+j−1−uminζ (inequality constraint (11c)), and h1ζk+j−1=Δumaxζ−Δuiζk+j−1 and h2ζk+j−1=Δuiζk+j−1−Δuminζ (inequality constraint (11d)), then the penalty cost function can be defined as follows:
(13)J¯iζk=Jiζk+μ∑j=1Hpl1ζk+jSl1ζk+j2+l2ζk+jSl2ζk+j2+λ∑j=1Hum1ζk+j−1Sm1ζk+j−12+m2ζk+j−1Sm2ζk+j−12+α∑j=1Huh1ζk+j−1Sh1ζk+j−12+h2ζk+j−1Sh2ζk+j−12,where
(14)Sa=s1000s2000s3,with a=a1,a2,a3T and si=105, if ai<0; otherwise, si=0 and i=1,2,3. The introduction of Sa makes it possible to consider the inequality constraints during the optimization.Similarly, the penalty cost function of the rotational subsystem can be given by(15)J¯iηk=Jiηk+μ∑j=1Hpl1ηk+jSl1ηk+j2+l2ηk+jSl2ηk+j2+λ∑j=1Hum1ηk+j−1Sm1ηk+j−12+m2ηk+j−1Sm2ηk+j−12+α∑j=1Huh1ηk+j−1Sh1ηk+j−12+h2ηk+j−1Sh2ηk+j−12.Based on (8) and (9), the constraints of Dsafe and Dcom are also taken into account. The corresponding auxiliary optimization term is described as
(16)Lk=σ⋅∑j=1Hp∑i≠rmax0,Dsafe−dirk+j+max0,dirk+j−Dcom2,where σ=105 and Lk as a penalty term is added in (13) according to the flight conditions of the multiquadrotor.Remark 2.
The cost functions of (13) and (15) are complex and nondifferentiable. So the optimization problem using traditional search techniques is not feasible. In this study, the DE algorithm [30] as a probabilistic search technique is proposed to solve the optimization problem for MPC. The corresponding details will be presented in Section 5.


## 4. RNN-Based Prediction Models

Although the feedback mechanism of MPC tolerates some model mismatches for multiquadrotor system, MPC still demands that the prediction model be sufficiently accurate. In this section, RNN is considered as the prediction model for MPC, due to its superiority in comparison with other conventional modeling methods [40].


### 4.1. Dynamical RNN

Generally, in a RNN, there is usually a feedback from the hidden layer output to the hidden layer input. This recurrent connection allows the NN to detect and generate time-varying patterns [41]. In this paper, the RNN predictive model is separately developed for the translational subsystem and rotational subsystem. Both subsystems belong to the class of multi-input/multioutput square nonlinear dynamical systems (that is, the systems with as many inputs as outputs), which have the general form of
(17)ẋ=fx+∑i=1mgixui=fx+Gxu,(18)y=Cx,where x∈Rn is the state, u∈Rm is the control input, y∈Rm is the system output, f:Rn→Rn is a smooth function, C is a m×n constant matrix, and Gx is a matrix with columns gi; note that f and gi contain parametric uncertainties.This paper deals with the control of both subsystems described by (17), i.e., m=3 and n=6. However, due to the complex environment and disturbances, the system (17) is difficult to use if the functions f and gi cannot be accurately known. Therefore, in order to provide a solution to this problem, a dynamical RNN is chosen to obtain a more accurate model of (17) in Figure 4.Figure 4
The block diagram of the dynamical RNN.This RNN containsn=6 dynamical hidden neurons, i.e., a number of neurons equal to the dimension of the state vector of subsystem dynamics. And it can be described by the following system [42, 43]:
(19)x^̇=Ax^+BW1Φx^+BW2Ξx^u,(20)y^=Cx^,with the state x^∈R6, the input u∈R3, the output y^∈R3, W1 is a 6×6 matrix of adjustable synaptic weights, A is a 6×6 diagonal matrix with negative eigenvalues ai, B is a 6×6 diagonal matrix with scalar elements bi, and W2 is a 6×6 diagonal matrix with adjustable synaptic weights. Φx^ is a 6-dimensional vector with elements of φx^i, and Ξx^ is a 6×3 matrix with elements of ξx^i. φx^i and ξx^i are represented by sigmoids of the form
(21)φx^i=m11+e−δ1x∧i,ξx^i=m21+e−δ2x∧i+κ,i=1,⋯,6,where m1, m2, δ1, and δ2 are constants representing the bound and slope of the sigmoid’s curvature and κ>0 is a constant that shifts the sigmoid, such that ξx^i>0 for all i=1,⋯,6. In this paper, the designed parameters are chosen as m1=0.5, m2=1.3, δ1=δ2=0.1, and κ=0.4.


### 4.2. Online Updating for RNN

Based on the above analysis, suppose that both subsystems of the quadrotor can be accurately described by the dynamical RNN (19) plus a modelling error term ωx,u as follows:
(22)ẋ=Ax+BW1∗Φx+BW2∗Ξxu+ωx,u,(23)y=Cxwhere W1∗ and W2∗ are the corresponding weight values.Note that the optimal control input is applied to both the real system and the predictive model. Define the error between the predictive model (19) and the real system (22) as e=y^−y=Cx^−x. Assuming ωx,u is zero, we obtain
(24)ė=CAC−1e+CBW~1Φx+CBW~2Ξxu,where W~1=W1−W1∗ and W~2=W2−W2∗, which are substituted online through the appropriate updating laws.To obtain the stable updating laws, the Lyapunov synthesis method is introduced. Consider a Lyapunov function candidate as(25)V=12⋅eTPe+trW~1TW~1+trW~2TW~2,where P>0 satisfies the Lyapunov equation PCA+ATCTP=−C.Differentiating (25) along the solution of (24), there exists
(26)V̇=−12⋅eTe+ΦTxW~1TBCTPe+uTΞTxW~2BCTPe+trẆ1TW~1+trẆ2TW~2.By choosing(27)trẆ1TW~1=−ΦTxW~1TBCTPe,trẆ2TW~2=−uTΞTxW~2BCTPe,then
(28)V̇=−12⋅eTe=−12⋅e2≤0.Consequently, based on (27), the weight matrices can be updated online over time horizons.Theorem 1.
The updating laws (27) can guarantee the following properties:
(i)
limt→∞et=0,
(ii)
limt→∞W~̇1t=0,
(iii)
limt→∞W~̇2t=0.Proof 1.
SinceV̇ in (28) is negative semidefinite, we have V∈L∞, which implies e∈L∞, W~1∈L∞, and W~2∈L∞. Furthermore, x^=x+C−1e is bounded. Since V is a nonincreasing function of time and bounded from below, there exists limt→∞Vt=V∞. By integrating V̇ from 0 to ∞, we obtain
(29)∫0∞12e2dt=V∞−V0<∞,which implies e∈L2. By definition, φxi and ξxi are bounded for all xi and all inputs to RNN are also bounded by assumption. From (24), we conclude ė∈L∞. Since e∈L2∩L∞ and ė∈L∞, based on Barbalat’s lemma, there exists limt→∞et=0. Due to the boundedness of u, Φx, Ξx, and the convergence of et to zero, we have that W~̇1 and W~̇2 also convergence to zero.Remark 3.
Combined with the MPC, updating online for RNN is executed over time horizons. With the actual outputsyiζk and yiηk obtained at the kth time horizon, the weight matrices of the RNNs are updated online to predict the subsystem output for the next time horizon.


## 4.1. Dynamical RNN

Generally, in a RNN, there is usually a feedback from the hidden layer output to the hidden layer input. This recurrent connection allows the NN to detect and generate time-varying patterns [41]. In this paper, the RNN predictive model is separately developed for the translational subsystem and rotational subsystem. Both subsystems belong to the class of multi-input/multioutput square nonlinear dynamical systems (that is, the systems with as many inputs as outputs), which have the general form of
(17)ẋ=fx+∑i=1mgixui=fx+Gxu,(18)y=Cx,where x∈Rn is the state, u∈Rm is the control input, y∈Rm is the system output, f:Rn→Rn is a smooth function, C is a m×n constant matrix, and Gx is a matrix with columns gi; note that f and gi contain parametric uncertainties.This paper deals with the control of both subsystems described by (17), i.e., m=3 and n=6. However, due to the complex environment and disturbances, the system (17) is difficult to use if the functions f and gi cannot be accurately known. Therefore, in order to provide a solution to this problem, a dynamical RNN is chosen to obtain a more accurate model of (17) in Figure 4.Figure 4
The block diagram of the dynamical RNN.This RNN containsn=6 dynamical hidden neurons, i.e., a number of neurons equal to the dimension of the state vector of subsystem dynamics. And it can be described by the following system [42, 43]:
(19)x^̇=Ax^+BW1Φx^+BW2Ξx^u,(20)y^=Cx^,with the state x^∈R6, the input u∈R3, the output y^∈R3, W1 is a 6×6 matrix of adjustable synaptic weights, A is a 6×6 diagonal matrix with negative eigenvalues ai, B is a 6×6 diagonal matrix with scalar elements bi, and W2 is a 6×6 diagonal matrix with adjustable synaptic weights. Φx^ is a 6-dimensional vector with elements of φx^i, and Ξx^ is a 6×3 matrix with elements of ξx^i. φx^i and ξx^i are represented by sigmoids of the form
(21)φx^i=m11+e−δ1x∧i,ξx^i=m21+e−δ2x∧i+κ,i=1,⋯,6,where m1, m2, δ1, and δ2 are constants representing the bound and slope of the sigmoid’s curvature and κ>0 is a constant that shifts the sigmoid, such that ξx^i>0 for all i=1,⋯,6. In this paper, the designed parameters are chosen as m1=0.5, m2=1.3, δ1=δ2=0.1, and κ=0.4.


## 4.2. Online Updating for RNN

Based on the above analysis, suppose that both subsystems of the quadrotor can be accurately described by the dynamical RNN (19) plus a modelling error term ωx,u as follows:
(22)ẋ=Ax+BW1∗Φx+BW2∗Ξxu+ωx,u,(23)y=Cxwhere W1∗ and W2∗ are the corresponding weight values.Note that the optimal control input is applied to both the real system and the predictive model. Define the error between the predictive model (19) and the real system (22) as e=y^−y=Cx^−x. Assuming ωx,u is zero, we obtain
(24)ė=CAC−1e+CBW~1Φx+CBW~2Ξxu,where W~1=W1−W1∗ and W~2=W2−W2∗, which are substituted online through the appropriate updating laws.To obtain the stable updating laws, the Lyapunov synthesis method is introduced. Consider a Lyapunov function candidate as(25)V=12⋅eTPe+trW~1TW~1+trW~2TW~2,where P>0 satisfies the Lyapunov equation PCA+ATCTP=−C.Differentiating (25) along the solution of (24), there exists
(26)V̇=−12⋅eTe+ΦTxW~1TBCTPe+uTΞTxW~2BCTPe+trẆ1TW~1+trẆ2TW~2.By choosing(27)trẆ1TW~1=−ΦTxW~1TBCTPe,trẆ2TW~2=−uTΞTxW~2BCTPe,then
(28)V̇=−12⋅eTe=−12⋅e2≤0.Consequently, based on (27), the weight matrices can be updated online over time horizons.Theorem 1.
The updating laws (27) can guarantee the following properties:
(i)
limt→∞et=0,
(ii)
limt→∞W~̇1t=0,
(iii)
limt→∞W~̇2t=0.Proof 1.
SinceV̇ in (28) is negative semidefinite, we have V∈L∞, which implies e∈L∞, W~1∈L∞, and W~2∈L∞. Furthermore, x^=x+C−1e is bounded. Since V is a nonincreasing function of time and bounded from below, there exists limt→∞Vt=V∞. By integrating V̇ from 0 to ∞, we obtain
(29)∫0∞12e2dt=V∞−V0<∞,which implies e∈L2. By definition, φxi and ξxi are bounded for all xi and all inputs to RNN are also bounded by assumption. From (24), we conclude ė∈L∞. Since e∈L2∩L∞ and ė∈L∞, based on Barbalat’s lemma, there exists limt→∞et=0. Due to the boundedness of u, Φx, Ξx, and the convergence of et to zero, we have that W~̇1 and W~̇2 also convergence to zero.Remark 3.
Combined with the MPC, updating online for RNN is executed over time horizons. With the actual outputsyiζk and yiηk obtained at the kth time horizon, the weight matrices of the RNNs are updated online to predict the subsystem output for the next time horizon.


## 5. Control Structure

The composite control structure based on RNN and MPC is proposed in Figure5 for each quadrotor, which consists of three main parts: the predictive model based on RNN, the designed cost function, and the evolutionary algorithm. In this figure, the RNN1 and RNN2 blocks represent the predictive models of the translational subsystem and the rotational subsystem, respectively. U~iζk≜argminJ¯iζk and U~iηk≜argminJ¯iηk are the optimal control sequence of the translational subsystem and the rotational subsystem, respectively, where the first elements, i.e., u~ξk and u~ηk, are taken as the control inputs for both subsystems.Figure 5
MPC structure based on RNN for each unmanned quadrotor.The evolutionary algorithm to introduce in this paper is the DE algorithm. It is an optimization algorithm which searches randomly in continuous space with three evolutionary operations, namely, mutation, recombination, and selection. Suppose that the initial population hasNP individuals inside the given D-dimensional search space. The solution vector in continuous space can be represented by viG=vi1,vi2,⋯,viD, i=1,⋯,NP, G=1,⋯,NC [30]. The optimization goal of the DE algorithm is to find the proper control variable that minimizes the cost (13) and (15).


### 5.1. Mutation Operation

As shown in Figure6, for each viG in the parent population, the mutation individual v^iG+1 is produced by
(30)v^iG+1=vr1G+F×vr2G−vr3G,where r1, r2, and r3 are randomly chosen and different from the running index i and F is a real constant scaling factor within 0,1, which controls the amplification of the differential variation vr2G−vr3G.Figure 6
Mutation process of the DE algorithm.


### 5.2. Recombination Operation

For eachviG, a trial individual v¯iG+1 is generated by coping components from the mutation vector v^iG+1 and the target vector viG in dependence, as shown in Figure 7. It can be formulated as follows:
(31)v¯ijG+1=v^ijG+1,ifrandb≤CR,vijG,otherwise,i=1,⋯,NP,j=1,⋯,D,where randb∈0,1 is a uniform random number and the recombination control parameter CR is a constant in the interval 0,1.Figure 7
Recombination process of the DE algorithm.


### 5.3. Selection Operation

Offspringv¯iG+1 competes one to one with its parent viG. The selection operation is expressed by
(32)viG+1=v¯iG+1,iffv¯iG+1≤fviG,viG,otherwise.This selection scheme ensures that the best control variable with the minimum cost can be retained when moving from one generation to the next, which results in a fast convergence behavior.Generally, the DE algorithm has been successfully applied in some fields due to its high search accuracy, robustness, and good convergence speed. However, when the optimization problem of multi-UAV tracking is complex and time-varying, selecting suitable mutation strategies and control parameters at different evolution stages is difficult. Therefore, in this paper, we adopt the ZEPDE algorithm [44], where the appropriate mutation strategy can be gradually adjusted by a roulette wheel, and suitable combinations of F and CR can be generated by zoning evolution. Each individual has its own control parameter combination and mutation strategy. And the corresponding details can be found in [44].The control algorithm can be summarized as follows:Step 1.
Load the state of theith quadrotor and its corresponding reference. Initialize the solution and generate NP individuals randomly.Step 2.
The decoded values of each individualUiζk and Uiηkin the population carried out the RNN1 and RNN2 models to predict the values of xik+j,yik+j,zik+jT and ϕik+j,θik+j,ψik+jT, j=1,⋯,Hp. Whereas, consecutive inputs to the RNN models can change their internal state in a detrimental way. Then, the DE algorithm is executed at fault due to this situation, leading the optimal control variables to not be found. To avoid this problem, the RNN models should return their initial condition before each individual carries out the RNN models.Step 3.
Calculate the fitness value of each individual using the predicted states ofxik+jyik+j,zik+jT and ϕik+j,θik+j,ψik+jT, j=1,⋯,Hp.Step 4.
The individual which has the best fitness value is recorded. Perform the selection method based on the fitness of individuals, and in this way, the population remains better individuals.Step 5.
Execute mutation, recombination to create new individuals in the population.Step 6.
Produce the prediction values for each new individual in the population and calculate their fitness values based on procedures in Steps 2 and 3. The best individual that has been recorded previously is located in the new generation instead of the individual with the worst fitness value.Step 7.
Repeat Steps 4-7 until convergence criteria of the evaluation function are met.Step 8.
The first elementsu~iζk and u~iηk of the best individuals U~iζk and U~iηk apply both the RNN models and the real quadrotor system. The corresponding outputs are obtained.Step 9.
Perform the feedback compensation in MPC and online updating for RNN. Then, go to Step 2 for the next time horizon.


## 5.1. Mutation Operation

As shown in Figure6, for each viG in the parent population, the mutation individual v^iG+1 is produced by
(30)v^iG+1=vr1G+F×vr2G−vr3G,where r1, r2, and r3 are randomly chosen and different from the running index i and F is a real constant scaling factor within 0,1, which controls the amplification of the differential variation vr2G−vr3G.Figure 6
Mutation process of the DE algorithm.


## 5.2. Recombination Operation

For eachviG, a trial individual v¯iG+1 is generated by coping components from the mutation vector v^iG+1 and the target vector viG in dependence, as shown in Figure 7. It can be formulated as follows:
(31)v¯ijG+1=v^ijG+1,ifrandb≤CR,vijG,otherwise,i=1,⋯,NP,j=1,⋯,D,where randb∈0,1 is a uniform random number and the recombination control parameter CR is a constant in the interval 0,1.Figure 7
Recombination process of the DE algorithm.


## 5.3. Selection Operation

Offspringv¯iG+1 competes one to one with its parent viG. The selection operation is expressed by
(32)viG+1=v¯iG+1,iffv¯iG+1≤fviG,viG,otherwise.This selection scheme ensures that the best control variable with the minimum cost can be retained when moving from one generation to the next, which results in a fast convergence behavior.Generally, the DE algorithm has been successfully applied in some fields due to its high search accuracy, robustness, and good convergence speed. However, when the optimization problem of multi-UAV tracking is complex and time-varying, selecting suitable mutation strategies and control parameters at different evolution stages is difficult. Therefore, in this paper, we adopt the ZEPDE algorithm [44], where the appropriate mutation strategy can be gradually adjusted by a roulette wheel, and suitable combinations of F and CR can be generated by zoning evolution. Each individual has its own control parameter combination and mutation strategy. And the corresponding details can be found in [44].The control algorithm can be summarized as follows:Step 1.
Load the state of theith quadrotor and its corresponding reference. Initialize the solution and generate NP individuals randomly.Step 2.
The decoded values of each individualUiζk and Uiηkin the population carried out the RNN1 and RNN2 models to predict the values of xik+j,yik+j,zik+jT and ϕik+j,θik+j,ψik+jT, j=1,⋯,Hp. Whereas, consecutive inputs to the RNN models can change their internal state in a detrimental way. Then, the DE algorithm is executed at fault due to this situation, leading the optimal control variables to not be found. To avoid this problem, the RNN models should return their initial condition before each individual carries out the RNN models.Step 3.
Calculate the fitness value of each individual using the predicted states ofxik+jyik+j,zik+jT and ϕik+j,θik+j,ψik+jT, j=1,⋯,Hp.Step 4.
The individual which has the best fitness value is recorded. Perform the selection method based on the fitness of individuals, and in this way, the population remains better individuals.Step 5.
Execute mutation, recombination to create new individuals in the population.Step 6.
Produce the prediction values for each new individual in the population and calculate their fitness values based on procedures in Steps 2 and 3. The best individual that has been recorded previously is located in the new generation instead of the individual with the worst fitness value.Step 7.
Repeat Steps 4-7 until convergence criteria of the evaluation function are met.Step 8.
The first elementsu~iζk and u~iηk of the best individuals U~iζk and U~iηk apply both the RNN models and the real quadrotor system. The corresponding outputs are obtained.Step 9.
Perform the feedback compensation in MPC and online updating for RNN. Then, go to Step 2 for the next time horizon.


## 6. Stability Analysis

In this section, we first give the stability of the proposed RNN-based MPC. Meanwhile, the convergence analysis of the DE algorithm is also discussed.


### 6.1. Stability Analysis of RNN-Based MPC

The quadrotor system in this paper is represented in state space, and the cost functions (11a) and (12a) are also based on states. Therefore, the stability of the proposed predictive scheme is investigated by checking the monotonicity of the cost function with respect to time [39].Problems (11a) and (12a) can be given in the following general form:
(33a)Jik=∑j=1Hpyipk+j−rik+jQ2+ρ∑j=1HuΔuik+j−1R2=∑j=1Hpeipk+jQ2+ρ∑j=1HuΔuik+j−1R2,(33b)s.t.ymin≤yipk+i≤ymax,(33c)umin≤uik+j−1≤umax,(33d)Δumin≤Δuik+j−1≤Δumax,(33e)eipk+Hp+j=0,∀j≥1,(33f)Δuik+Hu+j=0,∀j≥0.Assume thatU~ik=u~ik,⋯,u~ik+Hu−1T is the optimal control sequence at time k. Then, the corresponding cost function is
(34)J~ik=∑j=1Hpy~ipk+j−rik+jQ2+ρ∑j=1HuΔu~ik+j−1R2=∑j=1Hpe~ipk+jQ2+ρ∑j=1HuΔu~ik+j−1R2.Then, at timek+1, based on the optimal control sequence of time k, a suboptimal control sequence U~i∗k+1 is postulated as
(35)U~i∗k+1=u~i∗k+1⋮u~i∗k+Hu−1u~i∗k+Hu=u~ik⋮u~ik+Hu−1u~ik+Hu−1.For the suboptimal controlU~i∗k+1, the cost function can be described as
(36)J~i∗k+1=∑j=2Hp+1y~ip∗k+j−rik+jQ2+ρ∑j=1HuΔu~i∗k+jR2=∑j=2Hp+1e~ip∗k+jQ2+ρ∑j=1HuΔu~i∗k+jR2.Assuming that with the same control input, the predictionsyipk+j at time k+1 are the same as ones derived at time k, the difference between the J~ik and J~i∗k+1 can be defined as
(37)J~i∗k+1−J~ik=e~ip∗k+Hp+1Q2−e~ipk+1Q2−ρΔu~ikR2.Taking into account the terminal equality constraints (33e), it is obvious that e~ip∗k+Hp+1=0, then
(38)J~i∗k+1−J~ik=−e~ipk+1Q2−ρΔu~ikR2≤0.Furthermore, introduceU~ik+1 as the optimal solution further. It is obvious that J~i∗k+1<J~ik+1, so
(39)J~ik+1−J~ik<−eipk+1Q2−ρΔuikR2.It is clear that forρ≠0, the cost is monotonically decreased with respect to time while the proposed control system is stable.


### 6.2. Convergence Analysis of DE Algorithm

It is known that the MPC method divides the global control problem into a series of local optimizations and each local optimization problem is solved by the adaptive DE algorithm. In this section, the asymptotic convergence of the local optimization is first analyzed; on this basis, the global convergence analysis is presented as well.


#### 6.2.1. Local Convergence Analysis [45]

Suppose that at thekth time horizon, the optimal solution set of the cost function J¯k can be defined as follows:
(40)S∗k=v∗:J¯v∗∣k=minJ¯v∣k:v∈Sk,where Sk is a measurable search space and v∗ represents the optimal solution individual.However, condition (40) is hard to satisfy in the actual process. Thus, an expanded optimal solution set is considered as follows:
(41)Sε∗k=v~∈Sk:J¯v~∣k−J¯v∗∣k<ε,where ε>0 is a small positive real number. Denoting μ as the Lebesgue measure, we suppose that the μSε∗k>0 for each ε. And the expanded optimal solution (41) can be approximately regarded as the optimal one.Before analyzing the convergence of the local optimization, we need to introduce a concept of the convergence in probability as follows:DefineVG,G=1,2,⋯ as a population sequence generated by the adaptive DE algorithm for the local optimization problem, where VG=viG,i=1,2,⋯NP. It is required that the algorithm converges to the optimal solution set in probability if and only if
(42)limG→∞PVG∩Sε∗k≠∅=1,where P denotes the probability of an event.Assumption 1.
Based on the above definition, it is assumed that for each populationVG of the adaptive DE algorithm, there exists at least one individual vG such that
(43)PvG∈Sε∗k≥α>0,where α is a small positive value.Proof 2.
Suppose that an individualvrandG is generated randomly in each population VG. The probability that vrandG belongs to the optimal solution set is given by
(44)PvrandG∈Sε∗k=μSε∗kμSk=α>0.It means thatPvrandG∉Sε∗k=1−α>0. And the probability that the first d populations do not include an individual vG∈Sε∗k can be estimated by
(45)∏G=1dPVG∩Sε∗k=∅≤1−αd.Based on basic ideas of the DE algorithm, the best individual in the populationVG has the same or better fitness value than the best one from all the previous populations, implying
(46)limG→∞PVG∩Sε∗k=∅=limd→∞∏G=1dPVG∩Sε∗k=∅≤limd→∞1−αd=0,in other words,
(47)limG→∞PVG∩Sε∗k≠∅=1−limd→∞∏G=1dPVG∩Sε∗k=∅=1−0=1,which proves (42).


#### 6.2.2. Global Convergence Analysis

Based on the basic ideas of the MPC approach, the local optimizations follow the same control algorithm by using the adaptive DE algorithm. Therefore, the local convergence analysis of thekth time horizon can be applicable for the k+jth time horizon, j=1,2,⋯. Furthermore, the local optimization of the k+jth time horizon is executed on the basis of the optimization results of the k+j−1th time horizon. Finally, this local convergence of each time horizon leads to the global convergence of the entire control process.


## 6.1. Stability Analysis of RNN-Based MPC

The quadrotor system in this paper is represented in state space, and the cost functions (11a) and (12a) are also based on states. Therefore, the stability of the proposed predictive scheme is investigated by checking the monotonicity of the cost function with respect to time [39].Problems (11a) and (12a) can be given in the following general form:
(33a)Jik=∑j=1Hpyipk+j−rik+jQ2+ρ∑j=1HuΔuik+j−1R2=∑j=1Hpeipk+jQ2+ρ∑j=1HuΔuik+j−1R2,(33b)s.t.ymin≤yipk+i≤ymax,(33c)umin≤uik+j−1≤umax,(33d)Δumin≤Δuik+j−1≤Δumax,(33e)eipk+Hp+j=0,∀j≥1,(33f)Δuik+Hu+j=0,∀j≥0.Assume thatU~ik=u~ik,⋯,u~ik+Hu−1T is the optimal control sequence at time k. Then, the corresponding cost function is
(34)J~ik=∑j=1Hpy~ipk+j−rik+jQ2+ρ∑j=1HuΔu~ik+j−1R2=∑j=1Hpe~ipk+jQ2+ρ∑j=1HuΔu~ik+j−1R2.Then, at timek+1, based on the optimal control sequence of time k, a suboptimal control sequence U~i∗k+1 is postulated as
(35)U~i∗k+1=u~i∗k+1⋮u~i∗k+Hu−1u~i∗k+Hu=u~ik⋮u~ik+Hu−1u~ik+Hu−1.For the suboptimal controlU~i∗k+1, the cost function can be described as
(36)J~i∗k+1=∑j=2Hp+1y~ip∗k+j−rik+jQ2+ρ∑j=1HuΔu~i∗k+jR2=∑j=2Hp+1e~ip∗k+jQ2+ρ∑j=1HuΔu~i∗k+jR2.Assuming that with the same control input, the predictionsyipk+j at time k+1 are the same as ones derived at time k, the difference between the J~ik and J~i∗k+1 can be defined as
(37)J~i∗k+1−J~ik=e~ip∗k+Hp+1Q2−e~ipk+1Q2−ρΔu~ikR2.Taking into account the terminal equality constraints (33e), it is obvious that e~ip∗k+Hp+1=0, then
(38)J~i∗k+1−J~ik=−e~ipk+1Q2−ρΔu~ikR2≤0.Furthermore, introduceU~ik+1 as the optimal solution further. It is obvious that J~i∗k+1<J~ik+1, so
(39)J~ik+1−J~ik<−eipk+1Q2−ρΔuikR2.It is clear that forρ≠0, the cost is monotonically decreased with respect to time while the proposed control system is stable.


## 6.2. Convergence Analysis of DE Algorithm

It is known that the MPC method divides the global control problem into a series of local optimizations and each local optimization problem is solved by the adaptive DE algorithm. In this section, the asymptotic convergence of the local optimization is first analyzed; on this basis, the global convergence analysis is presented as well.


### 6.2.1. Local Convergence Analysis [45]

Suppose that at thekth time horizon, the optimal solution set of the cost function J¯k can be defined as follows:
(40)S∗k=v∗:J¯v∗∣k=minJ¯v∣k:v∈Sk,where Sk is a measurable search space and v∗ represents the optimal solution individual.However, condition (40) is hard to satisfy in the actual process. Thus, an expanded optimal solution set is considered as follows:
(41)Sε∗k=v~∈Sk:J¯v~∣k−J¯v∗∣k<ε,where ε>0 is a small positive real number. Denoting μ as the Lebesgue measure, we suppose that the μSε∗k>0 for each ε. And the expanded optimal solution (41) can be approximately regarded as the optimal one.Before analyzing the convergence of the local optimization, we need to introduce a concept of the convergence in probability as follows:DefineVG,G=1,2,⋯ as a population sequence generated by the adaptive DE algorithm for the local optimization problem, where VG=viG,i=1,2,⋯NP. It is required that the algorithm converges to the optimal solution set in probability if and only if
(42)limG→∞PVG∩Sε∗k≠∅=1,where P denotes the probability of an event.Assumption 1.
Based on the above definition, it is assumed that for each populationVG of the adaptive DE algorithm, there exists at least one individual vG such that
(43)PvG∈Sε∗k≥α>0,where α is a small positive value.Proof 2.
Suppose that an individualvrandG is generated randomly in each population VG. The probability that vrandG belongs to the optimal solution set is given by
(44)PvrandG∈Sε∗k=μSε∗kμSk=α>0.It means thatPvrandG∉Sε∗k=1−α>0. And the probability that the first d populations do not include an individual vG∈Sε∗k can be estimated by
(45)∏G=1dPVG∩Sε∗k=∅≤1−αd.Based on basic ideas of the DE algorithm, the best individual in the populationVG has the same or better fitness value than the best one from all the previous populations, implying
(46)limG→∞PVG∩Sε∗k=∅=limd→∞∏G=1dPVG∩Sε∗k=∅≤limd→∞1−αd=0,in other words,
(47)limG→∞PVG∩Sε∗k≠∅=1−limd→∞∏G=1dPVG∩Sε∗k=∅=1−0=1,which proves (42).


### 6.2.2. Global Convergence Analysis

Based on the basic ideas of the MPC approach, the local optimizations follow the same control algorithm by using the adaptive DE algorithm. Therefore, the local convergence analysis of thekth time horizon can be applicable for the k+jth time horizon, j=1,2,⋯. Furthermore, the local optimization of the k+jth time horizon is executed on the basis of the optimization results of the k+j−1th time horizon. Finally, this local convergence of each time horizon leads to the global convergence of the entire control process.


## 6.2.1. Local Convergence Analysis [45]

Suppose that at thekth time horizon, the optimal solution set of the cost function J¯k can be defined as follows:
(40)S∗k=v∗:J¯v∗∣k=minJ¯v∣k:v∈Sk,where Sk is a measurable search space and v∗ represents the optimal solution individual.However, condition (40) is hard to satisfy in the actual process. Thus, an expanded optimal solution set is considered as follows:
(41)Sε∗k=v~∈Sk:J¯v~∣k−J¯v∗∣k<ε,where ε>0 is a small positive real number. Denoting μ as the Lebesgue measure, we suppose that the μSε∗k>0 for each ε. And the expanded optimal solution (41) can be approximately regarded as the optimal one.Before analyzing the convergence of the local optimization, we need to introduce a concept of the convergence in probability as follows:DefineVG,G=1,2,⋯ as a population sequence generated by the adaptive DE algorithm for the local optimization problem, where VG=viG,i=1,2,⋯NP. It is required that the algorithm converges to the optimal solution set in probability if and only if
(42)limG→∞PVG∩Sε∗k≠∅=1,where P denotes the probability of an event.Assumption 1.
Based on the above definition, it is assumed that for each populationVG of the adaptive DE algorithm, there exists at least one individual vG such that
(43)PvG∈Sε∗k≥α>0,where α is a small positive value.Proof 2.
Suppose that an individualvrandG is generated randomly in each population VG. The probability that vrandG belongs to the optimal solution set is given by
(44)PvrandG∈Sε∗k=μSε∗kμSk=α>0.It means thatPvrandG∉Sε∗k=1−α>0. And the probability that the first d populations do not include an individual vG∈Sε∗k can be estimated by
(45)∏G=1dPVG∩Sε∗k=∅≤1−αd.Based on basic ideas of the DE algorithm, the best individual in the populationVG has the same or better fitness value than the best one from all the previous populations, implying
(46)limG→∞PVG∩Sε∗k=∅=limd→∞∏G=1dPVG∩Sε∗k=∅≤limd→∞1−αd=0,in other words,
(47)limG→∞PVG∩Sε∗k≠∅=1−limd→∞∏G=1dPVG∩Sε∗k=∅=1−0=1,which proves (42).


## 6.2.2. Global Convergence Analysis

Based on the basic ideas of the MPC approach, the local optimizations follow the same control algorithm by using the adaptive DE algorithm. Therefore, the local convergence analysis of thekth time horizon can be applicable for the k+jth time horizon, j=1,2,⋯. Furthermore, the local optimization of the k+jth time horizon is executed on the basis of the optimization results of the k+j−1th time horizon. Finally, this local convergence of each time horizon leads to the global convergence of the entire control process.


## 7. Simulation Analysis

The performance of RNN-based MPC for quadrotor formation flight is evaluated by computer simulations. All the simulations are programmed with MATLAB, version 2015, and are run on a PC with a clock speed of 2.8 GHz and 4 GB RAM, in a Microsoft Windows 8.1 environment.


### 7.1. Offline Training for RNN

Before the experiment, the weights of RNN should be determined preliminarily. The input output data for RNN training is obtained by simulating a quadrotor model [35] with a PID controller. A large training data set spanning over 1000 seconds is generated by using the random step references ranging from 5 seconds to 15 seconds. All system states and control inputs are restricted within the corresponding ranges. To enable the RNN model to learn the noisy dynamics in the presence of sensor noises, wind gusts and aerodynamic coupling effects, additive disturbances are included for data generation during the simulation. The disturbances diζ and diη are in the form of white noise with various intensity. And the network is adopted in the form of (19). Figure 8 shows one sample of the training data sets. Based on the test data, the training performance of RNN is tested in Figure 9, compared with the FNN model [20] and nonlinear partial-differential-equation-based model [35]. From the figures, the apparent advantage of the RNN model is the prediction accuracy. In other words, these results show the success of the dynamic RNN structure.Figure 8
One sample of the training data sets: (a) control inputs and system outputs of translation subsystem; (b) control inputs and system outputs of the rotation subsystem.
(a)
(b)Figure 9
Test performance of RNN.


### 7.2. Results and Discussion

A scenario is constructed to simulate a certain flight condition for the multiquadrotor formation. Assuming that a square formation is composed of five quadrotors, the relative position of each quadrotor is described in Figure10, and their initial positions are 0,5,0m, 10,5,0m, 5,0,0m, 0,−5,0m, and 10,−5,0m. Besides, the initial attitudes of all quadrotors are supposed to be ϕi,θi,ψi=0,0,0. Define the 3rd quadrotor as the leader flying along a reference trajectory called “elliptical”: x=5⋅cosπt/4, y=5⋅sinπt/4, z=0.5⋅t, and t∈0,50. And the other quadrotors as the followers fly with the following relative states: Δx1,3=−5,5,0,0,0,0, Δx2,3=5,5,0,0,0,0, Δx3,4=−5,−5,0,0,0,0, and Δx3,5=5,−5,0,0,0,0. To simulate the real flight environment, the bounded disturbances are introduced in rotation subsystem (white noise of intensity 0.3 rad/sec) and transition subsystem (white noise of intensity 0.5 m/sec). Moreover, the quadrotor’s safe anticollision and communication regions are Dsafe=1m, Dcom=50m, respectively.Figure 10
Initial formation of five quadrotors.The parameter settings of the proposed system condition are as follows: the dynamic parameters of unmanned quadrotor system in this paper are represented in Table1. The limits on the system states including the range of ϕ, θ, and ψ; the corresponding turn rates ϕ̇, θ̇, and ψ̇; and the control inputs U1, U2, U3, and U4 are shown in Table 2. Table 3 shows the DE algorithm-run parameters and the associated weight values of the fitness function.Table 1
Quadrotor dynamic parameters.
Parameters
Value
Massms (kg)
2.14
Arm lengthla (m)
0.24
Inertia aroundx-axis Ixx (kg⋅m2)
1.25
Inertia aroundy-axis Iyy (kg⋅m2)
1.25
Inertia aroundz-axis Izz (kg⋅m2)
2.5
Thrust factorb (N⋅s2)
0.3
Drag factord (Nm⋅s2)
0.014Table 2
Constraints of the system states and control inputs.
Parameters
Value
Range of roll angleϕ (rad)
−π/2,π/2
Range of pitch angleθ (rad)
−π/2,π/2
Range of yaw angleψ (rad)
−π,π
Range of turn ratesϕ̇, θ̇, and ψ̇ (rad/s)
−1,1
Range of the main thrustU1 (N)
0,60
Range of the torqueU2 (M)
−50,50
Range of the torqueU3 (M)
−50,50
Range of the torqueU4 (M)
−10,10Table 3
DE-run parameters and weights of the cost function.
Parameters
NP
NC
μ
λ
α
ρ
Q
R
Value
50
50
1
1
1
0.5
I3×3
I3×3In such a case, the MPC with the RNN online model is introduced for the formation-tracking problem over the prediction horizons, compared with the conventional MPC and traditional PID, where the parameters of the PID controller arekP=106.3, kI=530.9, and kD=414.5. For both conventional MPC and RNN-based MPC, the prediction horizon and control horizon are Hp=3 and Hu=2 and the sampling frequency fs is 4 Hz. In Figure 11, taking the tracking performance of the leader as an example, the results of three different control methods are presented. It is obvious that the performances of the RNN-based MPC approach with smallest tracking errors are better than the other two approaches. In particular, the control inputs of the PID controller are out of constraints at some time horizons.Figure 11
Comparison results: (a) the variation of attitudes; (b) the variation of positions.
(a)
(b)In order to reflect the control results intuitively, the tracking performance of the quadrotor formation can be measured using the integral of absolute error (IAE) [46].
(48)IAE=1N⋅HT∑i=1N∑k=1HTyiζk−riζkI3×32+yiηk−riηkI3×32,where N is the total number of quadrotors and HT is the total number of time horizons.The IAEs of the conventional MPC and RNN-based MPC are presented in Table4. Based on that, the tracking performance controlled by RNN-based MPC is better with the smaller errors. This is because the RNN is trained based on the test data and has learned the influences of unknown disturbances before experiments. On the other hand, due to the updating online for RNN, the proposed MPC method could estimate in a timely manner the disturbances improving the control performance.Table 4
IAEs of the proposed RNN-based MPC and conventional MPC.
Method
RNN-based MPC
Conventional MPC
IAE
1.3983
2.3453The 3-D performance of formation tracking with RNN-based MPC are depicted graphically in Figure12. The overall variations of attitude and position are shown in Figure 13. The corresponding errors of relative attitude (Figure 14(a)) and position (Figure 14(b)) between the leader quadrotor and followers are within the anticipated range. Generally, the proposed control system manages to accurately achieve formation tracking with fast and precise responses and effectively attenuate the mismatched disturbances within the acceptable error range.Figure 12
3-D performance of formation flight.Figure 13
Flight process: (a) the variation of attitudes; (b) the variation of positions.
(a)
(b)Figure 14
Tracking errors: (a) relative attitudes; (b) relative positions.
(a)
(b)


## 7.1. Offline Training for RNN

Before the experiment, the weights of RNN should be determined preliminarily. The input output data for RNN training is obtained by simulating a quadrotor model [35] with a PID controller. A large training data set spanning over 1000 seconds is generated by using the random step references ranging from 5 seconds to 15 seconds. All system states and control inputs are restricted within the corresponding ranges. To enable the RNN model to learn the noisy dynamics in the presence of sensor noises, wind gusts and aerodynamic coupling effects, additive disturbances are included for data generation during the simulation. The disturbances diζ and diη are in the form of white noise with various intensity. And the network is adopted in the form of (19). Figure 8 shows one sample of the training data sets. Based on the test data, the training performance of RNN is tested in Figure 9, compared with the FNN model [20] and nonlinear partial-differential-equation-based model [35]. From the figures, the apparent advantage of the RNN model is the prediction accuracy. In other words, these results show the success of the dynamic RNN structure.Figure 8
One sample of the training data sets: (a) control inputs and system outputs of translation subsystem; (b) control inputs and system outputs of the rotation subsystem.
(a)
(b)Figure 9
Test performance of RNN.


## 7.2. Results and Discussion

A scenario is constructed to simulate a certain flight condition for the multiquadrotor formation. Assuming that a square formation is composed of five quadrotors, the relative position of each quadrotor is described in Figure10, and their initial positions are 0,5,0m, 10,5,0m, 5,0,0m, 0,−5,0m, and 10,−5,0m. Besides, the initial attitudes of all quadrotors are supposed to be ϕi,θi,ψi=0,0,0. Define the 3rd quadrotor as the leader flying along a reference trajectory called “elliptical”: x=5⋅cosπt/4, y=5⋅sinπt/4, z=0.5⋅t, and t∈0,50. And the other quadrotors as the followers fly with the following relative states: Δx1,3=−5,5,0,0,0,0, Δx2,3=5,5,0,0,0,0, Δx3,4=−5,−5,0,0,0,0, and Δx3,5=5,−5,0,0,0,0. To simulate the real flight environment, the bounded disturbances are introduced in rotation subsystem (white noise of intensity 0.3 rad/sec) and transition subsystem (white noise of intensity 0.5 m/sec). Moreover, the quadrotor’s safe anticollision and communication regions are Dsafe=1m, Dcom=50m, respectively.Figure 10
Initial formation of five quadrotors.The parameter settings of the proposed system condition are as follows: the dynamic parameters of unmanned quadrotor system in this paper are represented in Table1. The limits on the system states including the range of ϕ, θ, and ψ; the corresponding turn rates ϕ̇, θ̇, and ψ̇; and the control inputs U1, U2, U3, and U4 are shown in Table 2. Table 3 shows the DE algorithm-run parameters and the associated weight values of the fitness function.Table 1
Quadrotor dynamic parameters.
Parameters
Value
Massms (kg)
2.14
Arm lengthla (m)
0.24
Inertia aroundx-axis Ixx (kg⋅m2)
1.25
Inertia aroundy-axis Iyy (kg⋅m2)
1.25
Inertia aroundz-axis Izz (kg⋅m2)
2.5
Thrust factorb (N⋅s2)
0.3
Drag factord (Nm⋅s2)
0.014Table 2
Constraints of the system states and control inputs.
Parameters
Value
Range of roll angleϕ (rad)
−π/2,π/2
Range of pitch angleθ (rad)
−π/2,π/2
Range of yaw angleψ (rad)
−π,π
Range of turn ratesϕ̇, θ̇, and ψ̇ (rad/s)
−1,1
Range of the main thrustU1 (N)
0,60
Range of the torqueU2 (M)
−50,50
Range of the torqueU3 (M)
−50,50
Range of the torqueU4 (M)
−10,10Table 3
DE-run parameters and weights of the cost function.
Parameters
NP
NC
μ
λ
α
ρ
Q
R
Value
50
50
1
1
1
0.5
I3×3
I3×3In such a case, the MPC with the RNN online model is introduced for the formation-tracking problem over the prediction horizons, compared with the conventional MPC and traditional PID, where the parameters of the PID controller arekP=106.3, kI=530.9, and kD=414.5. For both conventional MPC and RNN-based MPC, the prediction horizon and control horizon are Hp=3 and Hu=2 and the sampling frequency fs is 4 Hz. In Figure 11, taking the tracking performance of the leader as an example, the results of three different control methods are presented. It is obvious that the performances of the RNN-based MPC approach with smallest tracking errors are better than the other two approaches. In particular, the control inputs of the PID controller are out of constraints at some time horizons.Figure 11
Comparison results: (a) the variation of attitudes; (b) the variation of positions.
(a)
(b)In order to reflect the control results intuitively, the tracking performance of the quadrotor formation can be measured using the integral of absolute error (IAE) [46].
(48)IAE=1N⋅HT∑i=1N∑k=1HTyiζk−riζkI3×32+yiηk−riηkI3×32,where N is the total number of quadrotors and HT is the total number of time horizons.The IAEs of the conventional MPC and RNN-based MPC are presented in Table4. Based on that, the tracking performance controlled by RNN-based MPC is better with the smaller errors. This is because the RNN is trained based on the test data and has learned the influences of unknown disturbances before experiments. On the other hand, due to the updating online for RNN, the proposed MPC method could estimate in a timely manner the disturbances improving the control performance.Table 4
IAEs of the proposed RNN-based MPC and conventional MPC.
Method
RNN-based MPC
Conventional MPC
IAE
1.3983
2.3453The 3-D performance of formation tracking with RNN-based MPC are depicted graphically in Figure12. The overall variations of attitude and position are shown in Figure 13. The corresponding errors of relative attitude (Figure 14(a)) and position (Figure 14(b)) between the leader quadrotor and followers are within the anticipated range. Generally, the proposed control system manages to accurately achieve formation tracking with fast and precise responses and effectively attenuate the mismatched disturbances within the acceptable error range.Figure 12
3-D performance of formation flight.Figure 13
Flight process: (a) the variation of attitudes; (b) the variation of positions.
(a)
(b)Figure 14
Tracking errors: (a) relative attitudes; (b) relative positions.
(a)
(b)


## 8. Conclusion

In this paper, a novel predictive control scheme based on RNNs was proposed for the multiquadrotor system to achieve trajectory tracking and formation keeping simultaneously. Different from the conventional MPC, the proposed control approach designed the RNNs as the predictive model and used the adaptive DE algorithm to solve the local optimization problem for both subsystems of each quadrotor. Moreover, the closed-loop stability was analyzed; meanwhile, the convergence of the DE algorithm was discussed as well. Finally, the effectiveness of the designed control strategy was illustrated through comparative simulation results. However, the applicability of the proposed control structure is limited to the high computation burden of receding optimization. Therefore, an appealing extension of this work could focus on the time efficiency of whole control process.



---

*Source: 7272387-2019-06-09.xml*