Flatness-Based Aggressive Formation Tracking Control of Quadrotors with Finite-Time Estimated Feedforwards

: In this paper we address a decentralized neighbor-based formation tracking control of multiple quadrotors with leader–follower structure. Different from most of the existing work, the formation tracking controller is given in one loop without distinguishing the motion control and attitude control by means of the theory of ﬂatness. In order to achieve an aggressive formation tracking, the high-order states of the neighbors motion are estimated by using a proposed extended ﬁnite-time observer for each quadrotor. Then the estimated motion states are used as feedforwards in the formation controller design. Simulation and experimental results show that the proposed formation controller improves the formation performance, i.e., the formation pattern of the quadrotors is better maintained than that using the formation controller without high-order feedforwards, when tracking an aggressive reference formation trajectory.


Introduction
The cooperative control of multi-quadrotor systems has progressively attracted attention of researchers in both civil and military areas [1]. Especially, in the application of object transportation, the formation of multiple quadrotors has potential application significance [2][3][4][5][6].
The formation control problem can be categorized into two types. Namely, the formation producing problem and the formation tracking problem [7]. The former problem refers to the algorithm design for a group of vehicles to reach a predefined geometric pattern without a group reference [8][9][10], while the latter problem refers to the same task and meanwhile following a predefined trajectory [11][12][13]. The object transformation is a typical formation tracking problem, which is considered in this paper. An external reference formation trajectory (RFT) is given to the leading quadrotor(s) to guide the formation. We note that in the case of multiple leaders, they share the same RFT. According to [14,15], the formation tracking problem can be attributed as a consensus problem.
In most aggressive maneuver cases, the quadrotor is controlled in tracking a desired aggressive moving trajectory that is prescribed a prior, for instance, by polynomials w.r.t time. A widely used method is the control of geometry on Lie group (SE(3) for example). The aggressive formation problem of multiple quadrotors is investigated in [16], where the quadrotors are permitted to move quickly in 3-D environment with a tight formation. The formation of quadrotors is based on proper trajectory generation process, thereby, the full information of the trajectory (including each order derivatives) are known. However, in a multi-agent system with neighbor-based decentralized control, it is impossible to enforce the motion of quadrotor by generating polynomials-based trajectories.
The flatness-based control is frequently used in the trajectory planing for high-order dynamical system. For instance, the time-optimal trajectory generation problem is in-λ max (·) and λ min (·) represent the maximum and minimum eigenvalue of matrix inside the parenthesis, respectively.

Preliminaries
Some basic knowledge in graph theory and 'interaction matrix' is given in this section. Some propositions and corollaries about the interaction matrix are introduced, which are useful in the analysis of the consensus of quadrotors in Section 4.

Graph Theory
Graph theory is helpful to represent the interaction relations of multi-agent systems. A graph G = (V, E ) is normally with the sets of vertices V and edges E . The set of vertices V = {1, 2, . . . , n} is composed of the indices of agents. |V | represents the cardinality of the set V, which satisfies |V | = n. The set of edges is represented by E ⊆ V × V. If an edge exists between two vertices, the two vertices are called adjacent. In this work, simple and undirected graphs are considered.
The adjacency matrix of G is denoted by G A = [a ij ] ∈ R n×n , where a ij represents the connection of nodes i and j on a graph. Since the simple graph is considered, we have a ii = 0. Since the graph is undirected, we have a ij = a ji and a ij = 1 if (i, j) ∈ E , otherwise, a ij = 0. The degree matrix of G is denoted by G D = diag{∑ n j=1 a 1j , . . . , ∑ n j=1 a nj }. The neighbor set N i = {j ∈ V : (i, j) ∈ E } of agent i, is composed of the indices of the agents j, which has interaction with the agent i. In other words, if a ij > 0, then, agent j is a neighbor of agent i. The number of neighbors of agent i is equal to |N i |.
We define a diagonal matrix G L = diag{l 1 , . . . , l n }, which indicates the knowledge of reference formation trajectory (RFT) of each agent. If l i = 1, agent i can obtain the RFT (by sensing/detecting), i.e., a leader. Otherwise, if l i = 0, agent i is a follower, for i ∈ V. Then, the leader set is defined as V L = {i ∈ V : l i > 0}. The leader set V L ⊂ V is a subset of V, which contains the indices of the leaders. Particularly, all the quadrotors are leaders, when V L = V. The indices of the followers are contained in the complementary set of V L , namely, V − V L .

Interaction Matrix
We propose an interaction matrix G for representing the interconnecting relations of quadrotors with leader-follower (L-F) configuration as follows.
Remark 1. Let us note that the first two terms G D − G A are normally called the Laplacian in graph theory. Proposition 1. Let G be an undirected simple graph, then the interconnection matrix G in Equation (1), is positive-definite, if (i) G is connected; (ii) V L = Φ, where Φ represents a null set.
Considering the definition of G L , we have G L ≥ 0. We prove this proposition by contradiction. Firstly, we suppose that there exists a nonzero vector x ∈ R n , which renders Therefore, we must have x T Lx = 0 and x T G L x = 0. Since x T Lx = 0, we obtain that x = α1 n , where α is a nonzero scalar. According to the fact that V L = Φ, then, G L = 0. As a result, x T G L x = α 2 1 T n G L 1 n > 0, which contradicts x T G L x = 0. Therefore, such a nonzero vector x does not exist. Thus, for any nonzero vector x, x T Gx > 0, namely, G is positive-definite.
When condition (i) is not satisfied such that G is not connected, the graph can be divided into s connected sub-graphs G 1 , G 2 , . . . , G s . In that case, according to Proposition 1, we have the following corollary. Corollary 1. The interaction matrix G is positive-definite, if the leaders set satisfies V L ≥ s, and Proof. If a multi-agent system contains s connected sub-groups of agents, then, the interaction matrix is block diagonal, which has s blocks on the diagonal. Since each sub-group G i , i = {1, 2, . . . , s} is connected and V iL = Φ, then, the block in the interaction matrix is positive definite. Hence, the interaction matrix is positive definite.
The normalized interaction matrix is defined bȳ The normalized interaction matrixḠ is not symmetric, but it has the following property.

Proposition 2.
The normalized interaction matrixḠ has real eigenvalues, if G is undirected.
Proof. SinceḠ is a normalized interaction matrix, then,Ḡ = (G D + G L ) −1 · G. We recall that G is symmetric. Using (G D + G L ) − 1 2 , we apply similarity transformation toḠ, then, we obtain that , the eigenvalues ofḠ are real, moreover, they are equal to the eigenvalues of (G Corollary 2. The normalized interaction matrixḠ is positive-definite, if conditions in Proposition 1 or Corollary 1 are satisfied. Proof. According to the conditions in Proposition 1 or Corollary 1, we know that G is invertible. SinceḠ is similar to (G D + G L ) − 1 2 · G · (G D + G L ) − 1 2 according to Proposition 2, we obtain thatḠ is also positive-definite.

Quadrotor Model
The dynamics of a quadrotor is modelled as the motion of a rigid body in 3-D space under a thrust force and three moments, which are generated by the thrust forces of the four rotors. The notations used in the quadrotor's model are shown in Table 1. The orientation of the quadrotor with respect to the inertial frame is represented by the rotation matrix R i ∈ SO(3), where SO(3) represents a special orthogonal group, and whose determinant is one. Note that R T i R i = I, where I represents an identical matrix. The dynamics of a quadrotor i is shown as follows where X i = [X i , Y i , Z i ] T represents the coordinates of the center of mass of a quadrotor in the fixed inertial frame x e y e z e . The Euler angles (pitch, roll, and yaw) are represented by the vector Θ i = [φ i , θ i , ψ i ] T . The first equation of (3) represents the translational dynamics in the inertial frame, where m represents the mass of a quadrotor and g represents the gravity. The rotation matrix R i in the second equation of (3) can transform the coordinates of a point from the body-fixed frame to the inertial frame. The third equation of (3) represents the rotational dynamics of quadrotor i. The inertia matrix J of a quadrotor i is represented in the body-fixed frame, according to the assumptions on the physical structure of a quadrotor, we can obtain that J = diag{I x b , I y b , I z b } is diagonal, where scalars I x b , I y b and I z b represent the moments of inertia with respect to x b , y b and z b respectively. Additionally, we assume that all the quadrotors in the flock have the same physical coefficients such as the mass and the inertia matrix. The angular velocity of the quadrotor i in the body-fixed frame is represented by Ω i ∈ R 3 . The function S(·) : R 3 → R 3×3 represents an operation that transforms a vector in R 3 to a skew-symmetric matrix R 3×3 . Given two arbitrary vectors v 1 , v 2 ∈ R 3 , the function S(·) satisfies the property S(v 1 ) · v 2 = v 1 × v 2 . Then, according to the definition of cross product and the operation S(·), we conclude that S(Ω i ) ∈ R 3×3 is skew-symmetric matrix in terms of the members of Ω i = [p i , q i , r i ] T as follows where p i , q i and r i represent the angular velocities with respect to the body-fixed frame of the quadrotor i.

Symbol Description
x e y e z e Fixed inertial frame e 1 , e 2 , e 3 Basis vector of inertial frame Coordinates of the center of mass of a quadrotor in x e y e z e Euler angles (pitch, roll and yaw) Rotation matrices in yaw, pitch and roll Angular velocity of the quadrotor i in the body-fixed frame S(·) : R 3 → R 3×3 Operation from vector in R 3 to skew-symmetric matrix R 3×3 The moments of roll, pitch and yaw Control input of the quadrotor i We refer to φ i , θ i , ψ i as the roll, pitch, and yaw angles, represents the roll, pitch, and yaw moments. The thrust force is represented by F T i . We note that e 3 = [0, 0, 1] T is a constant vector. The attitude dynamics of each quadrotor is generated by moments τ φ i , τ θ i and τ ψ i and thrust force F T i .
The rotation matrix R i from the body-fixed frame to the inertial frame is the sequence represents the rotation matrices of yaw, pitch, and roll.
Thus, according to Equation (3), we obtain the translational dynamics as follows According to the second equation in (3), we can obtain where According to (6) and usingT i , we can obtain the rotational dynamics as follows The quadrotor is known as an under-actuated system. When the quadrotor tilts, a component of the total thrust is directed sideways and the aircraft translates horizontally. Therefore, the translational motion of quadrotor is generated by implicitly controlling its attitude (pitch and roll). This fact can be observed from (5), where the translation is controlled by the rotational angles.
In general, the control of the quadrotor is composed of inner and outer loop controllers, in which the outer loop is the translation control and the inner loop is the rotation control. Within this control framework, the control output of the outer loop is treated as an input of the rotational control loop. Nevertheless, in this paper, the idea is to get rid of the structure of inner-outer-loop control, we use the theory of flatness to directly design the torques of the quadrotor.
We assume that θ = ± π 2 . According to (3), (5), and (7) the quadrotor dynamics in inertial frame is represented in state space as Equation (8). We denote the state of a quadrotor i by κ i ∈ R 12 , which is defined by κ We denote the control input of the quadrotor i by The state space quadrotor i dynamics is given byκ where f : R 12 → R 12×12 and g : The output vector y i of the quadrotor i is composed by the positions, linear velocities, angles, and the angular velocities, which will be used in the controller design. Then, y i is given as follows where C i = C t C r i . Matrix C i is divided into two parts as follows It is important to note that the translational dynamics are modeled in the inertial frame while the rotational dynamics are modeled in the body-fixed frame, although the states (X i ,Ẋ i , Θ i andΘ i ) of the quadrotor model (8) are represented in the inertial frame. Indeed, in the body-fixed frame, the inertia matrix J is diagonal and then the rotational dynamics are easier to calculate than in the inertial frame, where J is not diagonal and time-variant. Note that in this work, the dynamics of the motors are omitted for the sake of simplicity.
The identification of the quadrotor parameters is not detailed here. Therefore, in the sequel, the parameters of the quadrotor, such as mass, inertias, and physical size are supposed to be known. For a quadrotor system, the flat output can be selected according to the following proposition.

Nonlinear-Flatness Based Decoupling Control
We firstly define the aggressive formation of quadrotors as follows.
Assumption 1 (Aggressive formation assumption). The attitude angles pitch θ i and roll φ i of each quadrotor can be greater than 15 • (but smaller than 30 • ) in aggressive formation modes.
According to Proposition 3, we conclude that the desired motion of a quadrotor can be uniquely specified by the 3D coordinates X i and the yaw angle ψ i . Therefore, the decentralized formation tracking problem is to find the guidance vector T for each quadrotor. The guidance vector is composed by the neighbors states, which will be designed in Section 4.
The control input where X d i and ψ d i represent the desired 3D position and the heading direction (yaw angle). According to Equation (A1), the controller for the altitude is designed as follows where Z d i represents the desired value of the altitude. Let us denote by e Z i = Z i − Z d i the tracking error of the altitude and substitute (10) into the third equation in (5). Then, we haveë Obviously, the altitude can exponentially track the given desired altitude trajectory Z d i for some selected positive scalars k 1z and k 2z . The desired altitude trajectory should satisfy Z d i ∈ C 3 (R, R). Then, the altitude dynamics are decoupled with the other states of the quadrotor.
We rewrite the dynamic of the attitude angles in (7) as follows Substituting (12) into (11), we obtainΘ it is trivial to observe that the attitude angles can exponentially track the desired attitude angles that are in terms of i and its derivatives. However, the controller (12) is open-loop on the translational dynamics, because the position feedback is not used in the controller design. If the dynamics of the quadrotor is precisely modeled and no external disturbances and sensors noise exist, this controller can perfectly perform. However, in practice, the dynamics of the system can be never precisely modeled, additionally, the noise and disturbance always exist. Therefore, the closed-loop control should be used, such that the feedback of the positions and translational velocities are required.
We observe from the equations in (A4) that the yaw angle is a component of the flat output i , then, the desired values of the yaw angle (ψ d i ,ψ d i andψ d i ) are trivial to obtain. However, the desired pitch and roll angles are not explicitly given by the components of the flat output d i and its derivatives. We will discuss how to obtain the terms As analyzed before, the design of τ ψ i is simple. Normally, in order to represent the control input by using the flatness output, the high-order derivatives of i are calculated. Specifically, in our case, in order to represent the torques τ φ i and τ θ i , two supplementary derivatives ofẌ i andŸ i in Equation (A2) are necessary and we obtain The function ζ represents the terms of the second-order derivative of the right side of Equation (A2) except the term that containsφ i andθ i .
The vector v i = [v i1 v i2 ] T is composed by the new control inputs v i1 and v i2 , which are decoupled on axes X and Y. Then, the new trajectory tracking controllers in closed loop on the plane X − Y are given as follows Observe that det A = (Z i + g)sec 2 φ i sec 3 θ i . Therefore, whenZ i = −g and |φ i | < π 2 , |θ i | < π 2 , det A is finite and positive such that A is invertible. In fact, if |φ i | < π 2 , |θ i | < π 2 , we haveZ i > −g.Z i will be never equal to −g, unless that all the rotors of the quadrotor stop rotating. This case is not considered in this work.
Then, according to Equation (13), we can design the torqueτ Then, substituting (15) into (12), we obtain the three torques. In practice, the velocitiesẊ i ,Ẏ i and its higher-order derivativesẌ i , ... X i ,Ÿ i and ... Yi are not directly measurable (not contained in the output y i (9)), they are obtained by using Equation (A2) and its derivatives. We will present in the following subsection how to obtain the these derivatives.
Under Assumption 1, we can assume that tan The desired yaw angle is zero. Using controller (11), ψ i ≈ 0. Then, cos ψ i ≈ 1 and sin ψ i ≈ 0. We also assume that the altitude is stabilized and keep constant, such that Z i ≈ Z d i and Z i ≈ 0. Then we have,Z i ≈ 0. Therefore, according to the analysis in Appendix A.2, we considerẌ Then, the controllers in (14) becomes Then, we obtain The control strategy in (18) is based on the knowledge of high-order derivatives of the planar positions X d i and Y d i in guidance vector, which are normally not measurable. Note that, for each quadrotor, the guidance vector depends on neighbors states (and the reference formation trajectory (RFT), if the quadrotor is a leader). If the quadrotors move slowly and the RFT varies slowly, we can assume that the higher-order derivatives are null. However, the performance will not be satisfactory for aggressive formations. The fast convergence of the observer is crucial in aggressive formation. Therefore, we will propose in the next subsection an extended finite-time observer to estimate these high-order derivatives.

Extended Finite-Time Observer
An observer of the high-order derivatives guidance vector for each quadrotor is proposed to get the feedforwards. We take X d i for example, according to v i1 in (17) Before designing the finite-time observer, we introduce the following lemma about finite-time stability. Lemma 1. [24] Consider the systemẋ = f (x, u). Suppose that there exist continuous differentiable functions V(x) : R n → R + , real numbers c > 0, and 0 < α < 1, such that V(0) = 0, V(x) > 0 for any x = 0, andV(x) ≤ −cV α (x). Then, for all x(t 0 ) = x 0 in the neighborhood µ of the origin, there exists settling time function T(x 0 ) : Denoting that where Then, we propose the extended observing model as followṡx where L t represents the gain matrix which satisfies then, according to (19) and (20) Denoting Suppose that J is the Jordan standard form of matrix A o by implementing a Jordan transformation using matrix M.
Proof. (i). The matrix A o can be rewritten in the following form.
(ii). The control gains in L t are selected such that the nearest pole of A o to imaginary axis is placed to −γ, where γ > 1. The Equation (21) can be rewritten asẋ d The quadratic term T J o inV( ) yields, where the intersection term ∑ According to (23) and (24), we have Then, the derivative of the Lyapunov function candidate yields, Therefore, . In addition,V( ) = 0, if and only if V( ) = 0. Note that V( ) = 1 2 T = 0, only if = 0 which is equivalent to that the observer error x d i = 0. Therefore, the observer error is asymptotic stable. According to Lemma 1, will converge into the compact sphere within the finitetime T converge .
The proof is completed.
According to (26), the observing error will converge to the invariant set ≤λ ρ γ−1 . Therefore, by increasing observer gains, which renders a greater γ, the invariant set will become smaller. The observing error will smaller as well. Nevertheless, the gains cannot be chosen arbitrarily great in practice, since the sensor noise will become significant and lead to instability.

Remark 3.
We note that since a decoupling quadrotor dynamics control law is derived in (16), the observer for the derivatives of Y d i can be designed independent with X d i .
The high-order derivatives are replaced by the estimated values using the observer (20). Then, the decoupled control law (17) becomes where the scalars k iX and k iY , i ∈ {1, 2, 3, 4} represent some selected positive gains. In the following section, we will show how to find the X d i and Y d i in the guidance vector (since the control of the altitude and yaw is decoupled in double integrators, it is omitted in the following section).

Formation Controller Design
Since the control of altitude (10) and rotation (12) is decoupled with the planar motion control, they are assumed to be controlled in constant, without loss of generality.

Available Guidance Vector
The formation controller of each quadrotor uses the relative states (positions, velocities, and the high-order derivatives with respect to its neighbors). If the quadrotor is a leader, besides the foregoing measurements, its formation controller also uses the relative states with respect to the RFT r(t) = [r X (t), r Y (t), r Z (t), r ψ (t)] T , where r Z and r ψ are assumed to be constant. Then, for the sake of simplicity, the formation control in (X, Y) space is considered.
The objective of the formation tracking control is to guarantee that all the quadrotors track the RFT with some constant biases d i0 = [d X i0 , d Y i0 ] ∈ R 2 , such that the quadrotors keep a formation pattern. In detail, the formation objective is to achieve the following equations.
for some initial condition x i (0), i = 1, . . . , n. An example of the formation tracking task of four quadrotors is shown in Figure 1 (left), where the solid red circle represents the RFT at time t a . The dashed red circles represent the desired positions (X d i (t a ),Ȳ d i (t a )) for the quadrotors i ∈ V at t a , where we denote bȳ The solid black circles represent the quadrotors' real positions at t a , i.e., (X i (t a ), Y i (t a )), i ∈ V. However, in decentralized formation control, the RFT r(t) is not available for the followers, therefore, it is not possible to design the formation tracking controller using the errorsX d i − X i andȲ d i − Y i for quadrotor i. For the followers, only the neighbors states are available for the formation controller design, as shown in Figure 1 (right). Now, the formation problem becomes: how to find the available guidance vector for quadrotors i ∈ V in order to attain the formation task.
Let us make a sum of the relative position state vectors. Note that we drop the explicit expression of time in the expressions for the sake of simplicity. Similar to the foregoing analysis, we still take X for example.
The inter-distance satisfies d X ij = d X i0 − d X j0 . Then, Equation (29) can be rewritten as follows Then, the guidance vector for each quadrotor is given as follows We then observe that X d i is available for quadrotor i, since X d i is composed only by neighbors states.
Then, according to definition of normalized interaction matrix in (2), we can rewrite Equation (31) in matrix form for all the quadrotors as follows whereḠ represents the normalized interaction matrix. According to Corollary 2 and the definition in Equation (2), we know thatḠ is invertible if each connected sub-graph of the multi-quadrotor system has at least one leader. Therefore, the formation task is achieved, i.e., X 1 −X d 1 , . . . , X n −X d n T → 0, as long as each quadrotor can precisely track the guidance vector X d i (t) and Y d i (t), i.e., vector

Convergence Analysis
T for each quadrotor is obtained (we note that Z d i and ψ d i are set constant, so omitted here). Their higher-order derivatives are estimated by using finite-time observer (20) in an aggressive formation. Then, these estimated derivatives are used as feedforwards in Equation (17), composing controller (12).
In practice, the estimation errors are usually high at t = 0 s, when the initial states of the observer may greatly deviate from the real values. In order to avoid the instability caused by the big initial deviation of the observer, we implement saturation functions to the estimated high-order derivatives to prevent great observation errors at the beginning.
To avoid the instability of quadrotors caused by the great initial observation error, we propose a modified controller in Equation (33) using the saturation function as follows The analysis of the convergence is given by the following theorem.

Theorem 2.
Consider the quadrotors with dynamics (8) are controlled by the formation controller (12) composed by (33). The bounds are chosen as follows, where j = {1, 2, 3, 4} andh j ∈ R + . Then, the formation error for each quadrotor is ultimately bounded and converge to the invariant set as follows.
Take X i for example and denote by Then, we can write the dynamics of e X i in state-space form as followṡ Therefore, the perturbed term δ in (36) is bounded and satisfying δ i ≤ δ max = max{d,κβ}, for all t ≥ t 0 .
The control gains are selected such that A c is Hurwitz. Since A c is Hurwitz, there exists a positive-definite matrix P, which solves the Lyapunov equation A T c P + PA c = −Q, for a positive-definite matrix Q. A Lyapunov function is selected as V = e T X i Pe X i . Then, its time derivative yieldṡ . When t > T convergence , according to LaSalle's invariance principle, for any initial condition x(t 0 ), the state x will finally stay in the set The error e X i in (36) is ultimately bounded. The same result on Y can be obtained similarly.

Remark 4.
According to Theorem 2, the tracking error of quadrotor is proportional to the bound β of the observer error after finite time T convergence . If the observer error asymptotically converges to zero, then the formation tracking error will also converge to zero asymptotically.

Remark 5.
In practice, before the time instant T convergence , the bounds b j are chosen small enough in order to prevent the unstable dynamics, since the observer error probably enormous, which is generated by great initial deviation.

Simulation and Experimental Results
In this subsection, we first use numerical simulation by means of MATLAB to qualitatively illustrate that the formation controller with estimated feedforwards has better performance than that without them. By doing this, we are able to show that our proposed formation controller makes the multi-quadrotor system performing a high bandwidth, such that the aggressive formation is achieved.
Then in the second subsection, the proposed control laws are implemented on the simulator-experiment platform comparing with that using the formation controller without high-order feedforwards.

Simulation Results
We reconsider the formation of four quadrotors in Figure 1 (left). The reference formation trajectory RFT is given by r(t) = [3 sin(0.1t), 3 cos(0.1t), 1, 0] T . The four quadrotors should keep a desired biases w.r.t. RFT defined by d i0 , i ∈ {1, 2, 3, 4}. As shown in Figure 1  The tuning process of the controller gains in Table 2 is given as follows.
• Controller gains tuning on simulator: .... Y d i be null. Give some small X d i and Y d i , for example a small step signal, then, tune the gains k jX and k jY , j = {1, 2, 3, 4}. Observe the oscillation of the rotation angles during translational motion, adjust the gains k 3X , k 4X and k 3Y , k 4Y to reduce the oscillation.
• Implement the controller gains on real quadrotor: Test on single quadrotor, slightly tune the parameters k 1X , k 2X if necessary. Then, if the performance is satisfied, implement them on the formation tracking control.  The observer gain matrix is selected according to the pole assignment as follows • Pole assignment: The gains of the observer are dominated by the finite time T converge that we want to have. Once it is fixed, the gains of the observer can be calculated by the technique of pole assignment, such that the nearest pole of A o to imaginary axis is placed to −γ. In the following simulation and experiments, we assume that the position and velocity in the guidance vector are measurable, i.e., We respectively show in Figures 2 and 3 the simulation results (i) without estimation of high-order derivatives and (ii) with estimation of high-order derivatives. We observe from Figure 2 that the quadrotors track the RFT with delays and the inter-distances are not well maintained. On the contrary, in Figure 3, the delays are small and the quadrotors keep around the desired inter-distances.
In the following subsection, we will implement the formation control protocol on the real-time simulator-experiment platform.

Experimental Setup
Heudiasyc laboratory (In the Université de Technologie de Compiègne) has developed a PC-based simulator-experiment framework (This framework is primarily designed and developed by engineer Guillaume Sanahuja under the help of other colleagues and PhD students in Heudiasyc laboratory) for controlling quadrotors in real-time. The programs (written in C++) running in the quadrotors are the same, both in the simulator and in the embedded processors of real quadrotors. The goal of the framework is to firstly validate a program on the simulator before implementing on real quadrotors, in order to avoid destroying the real drones. It is important to note that within this framework, the control algorithms are implemented on the quadrotors rather than on a PC. There does not exist a central controller that sends control signals to the quadrotors.
In the simulator, the complete quadrotor dynamics are used. The flight of the quadrotors are animated in a virtual 3D environment, which permits us to observe the behavior of the quadrotors under some formation control laws. This framework permits the simulation to reflect better the real-time experiment. Both the simulator and the real-time experiment share the ground station interface on the PC, which is responsible for displaying and sending instructions such as taking off and landing.
The experimental setup is shown in Figure 4. In the experiments, the motion capture system Optitrack is used to localize the quadrotors in the formation. The ArDrone2 quadrotors manufactured by Parrot are used for real-time experiments. The procedure of using the simulator-experiment framework is stated as follows: • Program the algorithms by using C++. • Compile the program into executable files for the quadrotors in the simulator and for the real quadrotors. • Test the algorithm in the simulator, adjust the parameters of the controller. • Send the executable file to each quadrotor. • Carry out the real-time experiment.
Remark 6. Although the system Optitrack is a global localization system, we just use it to obtain the coordinates or inter-distances of the quadrotors. For instance, a quadrotor can use its onboard cameras to detect other quadrotors and calculate their inter-distances. Additionally, the self-localization can be realized by using GPS or other sensors (laser, etc.). The system Optitrack is used here to simplify the study.

Real-Time Experiments on Simulator-Experiment Platform
In these experiments are the formation task of four quadrotors, where quadrotor 1 is a leader. The objective is to track the trajectory r(t) with constant biases d 10 = [0, , 0] T . Two tests are carried out on the simulator for the sake of comparison. In both tests, the objective is to track a circular trajectory (radius: 2 m; maximum linear velocity: 1.3 m/s; linear acceleration 0.1 m/s 2 ).
In the formation control without high-order feedforwards, the gains are chosen by k p = 1.5 and k d = 0.1. We note that these gains are not arbitrarily chosen. For the fairness of comparison, we tune the parameters on the simulator to obtain the gains with best performance.
We proposed a nonlinear-flatness-based formation controller in (12), where the highorder derivatives are estimated by using (20). In the simulator and the experiment, the control frequency is 50 Hz. Then, the states of the observer (20) are updated as followŝ where T = 1 50 = 0.02 s. The observer gains are given by (38). In Figure 5, we observe from Figure 5 that even at the beginning (when t < 20 s), the quadrotors are able to follow the given trajectory, after 20 s, the performance becomes not satisfactory, although the formation errors keep stable. On the contrary, by using the proposed formation control with the same controller gains, the performance is greatly improved and the quadrotors can always follow the given trajectory (see Figure 6).
The desired altitude is 1 m for all the quadrotors. The curves of the altitudes of the quadrotors are given in Figure 7. We can observe from Figure 7 (left) that the altitude of the quadrotor 2 oscillates after t = 25 s. It is important to note that this phenomena is caused by the actuator saturation of the quadrotors. Considering this problem, we have proposed several bounded formation controllers such as in [13], which will not be detailed in the current paper.   Figure 8, the formation control without high-order feedforwards is used, where the quadrotors fail to preserve the rectangle formation pattern after several second of flight. On the contrary, the nonlinearflatness-based formation control can preserve the rectangle formation and track the RFT that is only given to the leader.
A real-time experiment is also given in Figure 9, where four quadrotors are used in the tests. The reference is a circular trajectory. Two tests are given for comparison. We can observe that by using our proposed controller, the formation trajectory is more precisely followed. Figure 10 shows a photo of flight in a real-time experiment.

Conclusions
In this paper, a nonlinear-flatness-based formation control with feedforwards of estimated high-order guidance vector (neighbors' states, and RFT, if the quadrotor is a leader) derivatives is developed. The formation controller is decentralized, which is designed based on neighbors' states. In order to attain an aggressive formation, an extended finite-time observer of the guidance vector is investigated. In order to deal with the great initial deviation of the estimated states, the saturated feedforwards are developed. Then, the convergence of the formation tracking error is analyzed. Furthermore, the proposed formation control is illustrated by simulation and real-time experiment by means of MATLAB and the simulator-experiment platform, comparing with the formation control without high-order feedforwards. Both simulation and real-time experimental results show that the aggressive formation performance is improved by using our proposed formation controller.
One of the directions of the future work is to implement the proposed control with advanced nonlinear control methods, in order to deal with the effects of nonlinearities (uncertainties), to achieve asymptotic finite-time convergence, to control explicitly the formation tracking speed, etc.  Proof of Proposition 3. According to (5), we obtain that Replacing F T i by (A1) in the first two equations in (5), we obtain We assume thatZ i = −g. Then, according to Equation (A2), we obtain that tan φ i cos θ i =Ẍ i sin ψ i −Ÿ i cos ψ ï Z i + g tan θ i =Ẍ i cos ψ i +Ÿ i sin ψ ï Z i + g Then, we have Since i = [X T i , ψ i ] T , we denote i1 = X i , i2 = Y i , i3 = Z i and i4 = ψ i . Then, the state of the quadrotor κ i yields Equation (A4) imply that the state κ i can be parameterized by variable i and its derivatives. According to (9), the output y i is in terms of X i ,Ẋ i , Θ i andΘ i , therefore y i can also be parameterized by i . Additionally, according to (8), the control input u i is in terms of Z i ,Ż i , Θ i ,Θ i andΘ i . Therefore, the variable i is the flat output of system (8).

Appendix A.2. Model Simplification
According to Equation (A2), we can obtain the following equations θ i = arctan Ẍ i cos ψ i +Ÿ i sin ψ ï Z i + g φ i = arctan Ẍ i sin ψ i −Ÿ i cos ψ ï Z i + g cos arctan Ẍ i cos ψ i +Ÿ i sin ψ ï Z i + g Without loss of generality, we assume that the yaw angle is controlled at origin, and the altitude is stabilized at some constant value, then, θ i = arctan Ẍ i /g and φ i = arctan −(Ÿ i /g) cos θ i (A6) Furthermore, considering Assumption 1, we obtain thaẗ where the nonlinearities ∆ X i and ∆ Y i are bounded since the absolute values of θ i and φ i are bounded within 30 • by assumption 1. Then, we know that (A7) is Lyapunov stable, if the simplified system (16) is asymptotically stable.