Consensus-Based Distributed Formation Control of Multi-Quadcopter Systems: Barrier Lyapunov Function Approach

The problem of formation tracking control for a group of quadcopters with nonlinear dynamics using Barrier Lyapunov Functions (BLFs) is studied in this paper where the quadcopters are following a desired predefined trajectory in a predefined formation shape. The BLFs are employed to formulate the problem of formation trajectory tracking with a predefined accuracy. For this purpose, logarithmic BLFs including both the trajectory errors and the errors between the quadcopters’ distances with the desired ones (for the formation goal) are proposed. The method is firstly developed in a centralized scheme and then extended to a distributed framework using appropriate asymptotically convergent consensus algorithms. Therefore, the asymptotic convergence of the designed distributed algorithm to the centralized one is guaranteed. Moreover, due to the under-actuated feature of a quadcopter system, a general hierarchical scheme is considered for designing the controller. To this end, firstly a formation altitude tracking control is designed and then using the generated control signal, the formation translational tracking control is developed with the assumption of virtual inputs which are then employed to generate desired trajectory signals for the attitude control subsystem. Finally, attitude controllers are designed separately for each agent using the generated desired signals through logarithmic BLFs to consider a predefined accuracy. The efficiency of the proposed method is demonstrated through simulations and comparisons with the similar approaches in MATLAB-Simulink environment.


I. INTRODUCTION
Unmanned aerial vehicles (UAVs), such as quadcopters, are extensively used in many applications such as military, search and rescue, surveillance, monitoring and navigation in highly dangerous environments and especially in applications with extreme environments where humans cannot access [1], [2], [3].However, in practice, it is required for a group of UAVs to work with each other to handle a complex mission effectively.For example, it is almost impossible for a single quadcopter to search a vast area.Missions performed by a group of UAVs (generally robots) are more reliable, The associate editor coordinating the review of this manuscript and approving it for publication was Qiang Li .time-efficient, and cost-reducing [4].When a group of UAVs works together, they should travel individually while collaborating with each other to achieve a predefined task.For this purpose, a formation control scheme generates the appropriate control signals for a UAVs' formation to complete their mission.
The problem of formation control in multi-robot systems (both aerial and terrestrial ones) is the main attention of many recent research works [5], [6].The proposed algorithms can be categorized into two general centralized and decentralized (distributed) schemes where different formation mechanisms of leader-follower [7], virtual structure [8], and behaviorbased [9] can be implemented [4].In order to design a distributed scheme for each of the aforementioned structure, consensus algorithms [10], [11], [12] can be employed to coordinate the agents for the formation control purpose.Such approaches are normally named as consensus-based formation control methods [13], [14], [15], [16], [17].In the consensus-based formation control, the main purpose is coordination of the agents through consensus algorithms to reach to a given formation.The formation control can also be designed such that agents achieve a predefined shape or formation through tracking a predefined trajectory [18].Moreover, from the perspective of the sensed and controlled variables, the formation control can be implemented in position-based, displacement-based, and distance-based control schemes [19].
The problem of formation control in multi-quadcopter systems has been studied in some research works more recently.In [20], the problem of finite-time formation control of multi-quadcopter systems in a leader-follower structure is investigated, where the finite-time stability of the proposed distributed controller is shown and the proposed approach guarantees the formation trajectory control using the Lyapunov theory.In [21], a real-time nonlinear model predictive control (NMPC) based leader-follower formation trajectory control of quadcopters are tested on a fleet of CrazyFlie 2.0 quadcopters.The main contribution of the paper is to show that a C++ implementation of MPC as a controller, and not a planner, for real-time formation control of several quadcopters, is feasible.In [22], a model reference adaptive control (MRAC) based leader-follower formation producing problem is studied for a group of quadcopters, where the relative position in the x-y plane and the relative heading angles of the agents in the presence of disturbances are considered.In [23], a consensus-based leader-follower based formation tracking control for multi-quadcopter systems is presented where the impacts of perturbations are eliminated using an adaptive technique and definition of virtual control signals.
Barrier Lyapunov Functions (BLF) are basically bounded scalar positive definite functions on a given open region defined with respect to a system and approach to infinity on the boundary of that region.They are normally used to design a controller for a system such that state or output constraints are not violated [24].BLFs in single agent systems are mostly employed to address constraints in the states, especially in the transient time operation.Using BLFs for the formation control of multi-agent systems can be considered as a relative and appropriate idea as, according to [19], multiple agents are driven to achieve prescribed constraints on their states for the formation control.One of the main advantages of using BLFs in the formation control scenario is its ability to define a predefined accuracy for the controller which can help not only for addressing constraints on the states but also the ones caused by disturbances.As a result, more recently, the BLF framework is used for the formation tracking control in multi-robot systems where both constraints of the tracking errors (between the real and the desired trajectories for each robot) and the inter-robot distances are considered for a multi-robot system in 2D space (3 degrees of freedom) [25].The proposed method is basically distributed and exponential convergence is guaranteed on the distance tracking errors.
BLFs have been employed for the trajectory control of single quadcopters under state constraints more recently [26], [27], [28], [29].For example, in [26] the problem of 3-D trajetory following for a nonlinear under-actuated quadcopter dynamics is proposed in which constraints are imposed on its states using BLFs.The constraints include singularity in Euler angles and a predefined accuracy in the trajectory following.However, in spite of many existing research works in formation control of multi-quadcopter systems, utilizing BLFs for this purpose and in multi-quadcopter systems has not been considered in the literature.
Toward this, in this paper, the BLFs are employed for the formation tracking control of multi-quadcopter systems.This goal is implemented through the steps of the backstepping control design framework [30], [31].The method is firstly developed in a centralized scheme and then extended to a distributed framework using appropriate consensus algorithms.Therefore, as the consensus algorithms asymptotically converge to the desired values, the performance of the proposed distributed method asymptotically converges to the centralized one.Logarithmic BLFs are proposed for our purpose in which both the trajectory errors and the errors between the agents' distances and the desired distances for the formation of the agents are included using the Laplacian of the connectivity graph.Since the quadcopter system is under-actuated, a hierarchical method is employed where first a formation control method is developed for the altitude subsystem and then using the generated control signal, the formation control is designed for the translational subsystem with virtual inputs which are then employed to generate desired trajectory signals for the attitude control subsystem.
The key novelty of our proposed method is by utilizing BLFs in the formation tracking control of multi-quadcopter systems, one can impose prescribed accuracy constraints on the formation control as well as predefined tracking accuracy.The contributions of this paper are summarized as follows: (1) A centralized hierarchical formation tracking control is designed for the altitude and translational subsystems of a multi-quadcopter system using a proposed logarithmic BLFs in a backstepping procedure.
(2) The designed controllers are extended to distributed ones using asymptotically convergent consensus algorithms.
(3) The attitude controllers are designed separately for each agent using the generated desired signals via BLFs to consider a predefined accuracy and the state constraints.
The rest of this paper is organized as follows.In Section II, the problem of multi-quadcopter formation control is introduced and necessary preliminaries are presented.The proposed centralized position control of multi-quadcopter system is derived in Section III.The proposed centralized controller is then transformed to a distributed one in Section IV.The attitude controllers for each agent is designed in Section V. Simulation results are presented in Section VI.Finally, Section VII provides conclusions and future work.

II. PRELIMINARIES AND PROBLEM STATEMENT
In this paper, a group of quadcopters is considered which are communicating with each other through a wireless communication network as depicted in Figure 1.This framework can be modeled using a connectivity graph considering each quadcopter (agent) as a node and each communication as an edge of the graph.

A. GRAPH THEORY AND CONSENSUS PROBLEM
Let G = (V, E) be an undirected graph with a non-negative adjacency matrix A = [a ij ] specifying the interconnection topology of a multi-agent (quadcopter) system where V = {ν 1 , . . ., ν N } and E ⊆ V × V are the set of nodes and edges, respectively.Besides, let N i = {ν j ∈ V : a ij ̸ = 0} denotes the set of neighbors of node i, and L is the corresponding Laplacian matrix [32].
Remark 1: It is assumed that the connectivity graph of the quadcopters (related to the quadcopters connection topology) is a connected undirected graph.Now, let q i (t) ∈ R m denote the value of node ν i which represents the i-agent state vector of a high-pass consensus filter.The consensus filter is a dynamical system with the input and output vectors r i (t) ∈ R m and l i (t) ∈ R m , respectively, for i-th node, with the following dynamical model: Let Then, the augmented consensus system of (1) can be expressed as follows [33]: where L = L ⊗ I m and ⊗ refers to the Kronecker product and I m is an m-dimensional identity matrix.It is worth mentioning that through computation of the improper MIMO transfer function of the consensus filter presented in (2), it can be easily shown that it is a high-pass consensus filter.Using this high-pass consensus filter in a connected undirected graph, the asymptotic convergence of l i (t) to 1

B. BARRIER LYAPUNOV THEORY
To track a desired trajectory with a predefined accuracy, Barrier Lyapunov Functions (BLF) are employed in this paper.Consider the following continuous-time nonlinear system: where x(t) ∈ R n is the system state vector and f : R n → R n is the system nonlinear vector function.Then, according to [24], a BLF is a scalar continuous positive definite function V : D → R, D ⊂ R n with respect to the system (3) which has continuous first-order partial derivatives at every point of D and V (x) → ∞ as x approaches the boundary of D.Moreover, V (x(t)) ≤ b, ∀t ≥ 0 along the solution of (3) for x(0) ∈ D with a positive constant b.Now, consider a nonlinear system as: where x 1 (t) ∈ R n 1 , x 2 (t) ∈ R n 2 are the system states, u(t) ∈ R n 2 is the system input, and f 1 , g 1 , f 2 , and g 2 are smooth vector functions.The goal is x 1 (t) to follow the desired trajectory x 1d (t) with a predefined accuracy.In other words, if one lets e 1 (t) := x 1 (t) − x 1d (t), then the goal is where k b1 is a positive predefined scalar.By extension of the idea of using BLF in a backstepping procedure as presented in [24], [31], [35] to the vector form case, the following BLF can be considered for (4): where V 2 (e 2 (t)) = 1 2 e T 2 (t)e 2 (t) with e 2 (t) := x 2 (t) − α(t) and α(t) is a stabilizing vector function to be designed.According to Lemma 1 in [24] if the inequality of V (e 1 (t), e 2 (t)

C. QUADCOPTER MODEL
A quadrotor UAV includes four rotors that generate propeller forces.Schematic of the quadcopter with body coordinate system is depicted in Figure 2. Variations on forces (f 1 to f 4 ) and moments (τ M 1 to τ M 4 ) by adjusting rotors' speeds (ω 1 to ω 4 ) produce attitude and translation change in the quadrotor [1], [36].Let a group of quadcopters (as depicted in Figure 1) consisting of N quadcopters communicating with each other constructs a connected undirected connectivity graph.Consider the index i for the i-th quadcopter parameters, i.e. f 1i to f 4i , τ Mi 1 to τ Mi 4 and ω 1i to ω 4i .
The quadrotor attitude dynamics for the i-th quadcopter, i = 1, . . ., N can be written in the following form: 142918 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where p i (t), q i (t), and r i (t) are angular velocities rotating around x-axis, y-axis, and z-axis in the body frame, respectively, the input signals u 2i (t), u 3i (t), and u 4i (t) are torques in the direction of the corresponding body frame angles.It is also worth mentioning that where k i is the lift constant and L i is the length between center of the mass and the rotor axes.Moreover, a 1i = The relationship between angular velocities in the body and inertial frames are described as follows: where the roll angle φ i (t), pitch angle θ i (t), and yaw angle ψ i (t) determine rotations around x-axis, y-axis, and z-axis, respectively.Finally, the quadcopter translational dynamics is presented in the following: where [x i (t) y i (t) z i (t)] T represents the position of quadcopter in the inertial frame, u 1i (t) defines the main thrust created by combined forces of rotors, that is u 1i (t) = 4 l=1 f li (t), g is the gravity constant, and m i is the mass of the i-th quadcopter.Now, in the rest of the paper, it is intended to solve the following problems: Problem 1: How to design a centralized controller to control the position (altitude and translation) of a group of quadcopters to follow a desired trajectory in a predefined formation when the position subsystem is under-actuated?
Problem 2: How to implement the proposed centralized formation control algorithm in a distributed manner?
Problem 3: How to control the attitude subsystems of each agent to track the generated desired trajectory from previous steps with a predefined accuracy?

III. CENTRALIZED POSITION CONTROL
In this section, to solve Problem 1, a centralized controller is designed for the position subsystem.As the position subsystem is under-actuated, at first, the altitude subsystem is controlled and then the translational subsystem is controlled using virtual inputs.It is worth mentioning that, this general hierarchical framework is similar to the one proposed in [26].However, the main difference is that our proposed method has been extended and developed for the formation tracking control of a multi quadcopter system.In other words, several agents should be able to track the desired trajectory in a formation paradigm.

A. SYSTEM DESCRIPTION AND CONTROL GOALS
In this subsection, the centralized model of altitude and translational subsystems are firstly introduced and then the control objectives for trajectory tracking and formation control are presented.

1) SYSTEM DESCRIPTION
The altitude subsystem, as presented in (8), is firstly reformulated in the form of continuous-time state space model for the i-th agent as: where ] T , then the following centralized model is obtained: where G = g1 N ×1 and It should be noted that in the singular orientation of cos(φ i (t)) cos(θ i (t)) = 0, the quadrotor altitude control is not possible.Hence, this singular orientation should be avoided in the attitude control.Now, in order to model the translational subsystem, since the position subsystem is under-actuated, two virtual inputs are considered for the translational subsystem according to (8) as follows: Reformulating the translational subsystem according to (8) and (11) in the form of the state space model for the i-th subsystem, one can conclude: where The collective state-space model is then can be formulated as: where ⊗ refers to the Kronecker product, It should be noted that u 1i (t) ̸ = 0, ∀t ≥ 0 as it is assumed that the rotors are not ever turn off during the maneuver and a minimum trust should be generated by the rotors for hovering.

2) CONTROL OBJECTIVES
Two different types of control objectives are considered in this paper.The first one is devoted to the tracking problem and the second one is related to the formation control.
Let X 1d (t) = [x 11d (t), . . ., x 1Nd (t)] T where x 1id (t), i = 1, . . ., N are continuous-time desired altitude trajectory signals with bounded first and second order derivatives and let E 1 (t) := X 1 (t) − X 1d (t) where E 1 (t) = [e 11 (t), . . ., e 1N (t)] T .The control goal for tracking the desired altitude trajectory with a predefined accuracy is to design the control input U 1 such that: where k b1 is a positive predefined scalar.Moreover, let X 3d (t) = [x T 31d (t), . . ., x T 3Nd (t)] T where x 3id (t)| N i=1 are continuous-time desired translational trajectory vector signals with bounded first and second order derivatives.Besides, let E 3 (t) := X 3 (t) − X 3d (t) where E 3 (t) = [e T 31 (t), . . ., e T 3N (t)] T .It is supposed to design U v (t) such that the following control goals are fulfilled.Firstly, the desired translational trajectory X 3d (t) should be tracked with a predefined accuracy: where k b3 is a positive predefined scalar Moreover, the following conditions are considered to keep the quadcopters in a desired formation with a predefined accuracy: i,j∈E where k e1 and k e3 are positive predefined scalars.

B. ALTITUDE CONTROL
In this section, the altitude control problem is studied considering both trajectory tracking and formation control as presented in ( 14) and ( 16), respectively.Theorem 1: Consider the altitude dynamic (10) with the following control input where with k 1 and k 2 as positive constants.Then, the state and formation constraints ( 14) and ( 16) Therefore, since x 1i (t) − x 1j (t) − x 1id (t) + x 1jd (t) = e 1i (t) − e 1j (t), the condition in ( 16) can be replaced with the following one: i,j∈E Therefore, using the Laplacian matrix properties (L) the following is obtained: Now, the following BLF candidate is chosen: For the sake of simplicity, let a b1 (t Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. where A 2 (t) is stabilizing vector parameter defined in ( 19) and E 2 (t) = [e 21 (t), . . ., e 2N (t)] T .Therefore: Since Ė1 (t) = X 2 (t) − Ẋ1d (t) = E 2 (t) + A 2 (t) − Ẋ1d (t), ( 24) can be rewritten as follow: where I N is an N × N identity matrix.The selection of A 2 (t) in ( 19) leads to: Now, a Lyapunov function candidate by augmenting V 1 (E 1 (t)) with a quadratic function is chosen in the next step to cancel out the first term in V1 (t) through an appropriate selection of the control input U 1 (t), namely: Then, it follows that: ). ( 28) Substituting ( 18) and ( 26) in (28) leads to: which proves that ( 14) and ( 22) holds, according to the barrier Lyapunov theory as presented in II-B, ∀t > 0 provided □ Remark 2: It should be noted that a e1 (t)I N + a b1 (t)L is a positive definite matrix due to the fact that L is positive semi-definite and a e1 (t) and a b1 (t) keep staying positive as e1 and this leads to V2 (0) < 0 as follows and this causes

IV. DISTRIBUTED IMPLEMENTATION OF THE PROPOSED POSITION CONTROL
In this section, Problem 2 is studied and it is shown how to implement the derived control laws in the previous section, i.e. ( 18) and (30), in a distributed way in each agent communicating with the neighboring ones.

A. DISTRIBUTED ALTITUDE CONTROL
At first, consider the designed control law for the altitude subsystem, that is U 1 (t) in (18).The i-th element of U 1 (t) is the altitude control law for i-th agent, u 1i (t), which can be written as: where α 2i (t) is the i-th element of A 2 (t).It is also worth reminding that B 2 (t) is diagonal matrix and thus only [B 2 (t)] −1 ii gets involved in u 1i (t) computation.According to (35), the i-th agent is aware of its own and its neighboring states, that is e 1i (t), e 1j (t), j ∈ N i and e 2i (t).However, an important issue is how to compute α 2i (t), a b1 (t), and a e1 (t) in a decentralized way.Toward this, it is proposed in this paper to employ consensus algorithms as introduced in II-A.Secondly, to compute E T 1 (t)LE 1 (t), the following equality is employed: where x = [x 1 , . . ., x N ] T is a N × 1 vector.Toward the goal, each agent computes firstly r ia e1 (t) = j∈N i (e 1i (t) − e 1j (t)) 2 , i = 1, . . ., N , and also receives the corresponding computed r ja e1 (t) from its neighbors and the values are substituted in (1) with the consensus state of q ia e1 (t).Again, clearly, the result of l ia e1 (t) should be multiplied by N 2 .In other words, a e1i (t) = k 2 e1 − N 2 l ia e1 (t) which asymptotically converges to a e1 (t).

2) DISTRIBUTED COMPUTATION OF α 2i (t )
The main issue in computation of α 2i (t) is how to compute the i-th element of (a e1 (t)I N + a b1 (t)L) −1 E 1 (t) according to (19) as all of the agents' states (e 1i (t), i = 1, . . ., N ) are getting involved and only neighboring states are available.For this purpose, each agent computes a be1i (t) := (a e1i (t)I N + a b1i (t)L) −1 and put r i (t) = [a be1i (t)] ii e 1i (t) and r j (t) = [a be1i (t)] ij e 1j (t), j ∈ N i in (1) with the consensus state of q iα 2i (t).Finally, the filter output computed l iα 2i (t) converges to 1 N N j=1 [a be1 ] ij (t)e 1j (t).Hence, α 2i (t) can be obtained as: Besides, it is also required for each agent i to incorporate in distributed computation of α 2l (t), l = 1, . . ., N , l ̸ = i through computation of consensus state of q iα 2l (t) for l = 1, . . ., N , l ̸ = i with inputs of r j (t) = [a be1i (t)] lj e 1j (t) and r i (t) = [a be1i (t)] li e 1i (t).In order to clarify the proposed distributed altitude control algorithm, Table 1 describes the method in details as a pseudo code.

B. DISTRIBUTED TRANSLATIONAL CONTROL
Now, consider the designed control law for the translational subsystem, that is U v (t) in (30).The i-th element of U v (t) is the altitude control law for i-th agent, u vi (t), which can be written as: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
of its own states e 3i (t) and e 4i (t) and the neighboring states e 3j (t), j ∈ N i .

1) DISTRIBUTED COMPUTATION OF a b3 (t ) AND a e3 (t )
Similar to a b1 (t), a b3 (t) can be calculated by computing E T 3 (t)E 3 (t) using the consensus filter of (1).Toward this, r i (t) = e T 3i (t)e 3i (t) and the corresponding ones in neighboring agents, r j (t) = e T 3j (t)e 3j (t), j ∈ N i , are computed and are substituted in (1) with the consensus state of q ia b3 (t).The output of the consensus system, l ia b3 (t) is then multiplied by N and is used to compute a b3i (t) = k 2 b3 − Nl ia b3 (t).Moreover, in order to compute a e3 (t), it is required to compute k 2 e3 − E T 3 (t)(L ⊗ I 2 )E 3 (t) in a distributed way.Similar to (36), consider the following equation: where y = [y T 1 , . . ., y T N ] T is a 2N × 1 matrix where y i = [y ix , y iy ] T is a 2 × 1 matrix.Accordingly, each agent computes r ia e3 (t) = j∈N i ∥e 3i (t) − e 3j (t)∥ 2 and receives the corresponding ones r ja e3 (t), j ∈ N i and uses (1) to compute 2 with the consensus state of q ia e3 (t).Then, multiplying the consensus system output, l ia e3 (t), by N  2 , is computed which asymptotically converges to a e3 (t).

V. ATTITUDE CONTROL
In this section, we are going to solve Problem 3 in which the proposed attitude control is presented for each agent.Toward this, using u vi (t) = [u xi (t), u vi (t)] T in (38), the desired trajectory of φ id (t) and θ id (t) can be generated as follows according to (11): and ψ id (t) can be freely chosen.It is worth mentioning that similar to [26], BLFs are employed to control the attitude subsystem.However, the main and most important difference is that in [26], the equation ( 7), representing the relationship between angular velocities in the body and inertial frames has not been considered.Therefore, the system structure for the attitude control is more complicated and realistic in this paper.Toward this goal, ( 6) and ( 7) are firstly rewritten as follows: where As presented in Section II-B, the desired trajectory vector x 5id (t) = [φ id (t), θ id (t), ψ id (t)] T is firstly defined and the goal is ∥e 5i (t)∥ < k b5i where e 5i (t) := x 5i (t) − x 5id (t) and k b5i is a positive constant.Moreover, e 6i (t) := x 6i (t) − A 6i (t).Now, the following BLF candidate is selected as presented in (5): Therefore: where a b5i (t) = ∥k b5i ∥ 2 −e T 5i (t)e 5i (t).Thus, A 6i (t) is designed as follows: where k 5i is a positive constant scalar.This selection gives: where the first term should be canceled out in the next step.Now, a quadratic Lyapunov function is augmented to V 5i (e 5i (t)) to cancel the additional term as follows: This selection gives: ). (48) Therefore, u i (t) is designed as follows:   where k 6i is a positive constant scalar.Therefore, V6i (t) = −k 5i e T 5i (t)e 5i (t) − k 6i e T 6i (t)e 6i (t) < 0 and it can be concluded that ∥e 5i (t)∥ < k b5i , ∀t ≥ 0, if ∥e 5i (0)∥ < k b5i .
In order to clarify the proposed method, Figure 3 depicts the general framework of the proposed method for the i-th agent where dashed lines are related to data that should be sent to or received by the neighboring agents and solid lines are the internal signals where the blue ones are related to the desired signals, the green ones are the control signals and the red ones are the states measurements from the quadcopters.

1) CENTRALIZED FORMATION CONTROL
Firstly, the results for the centralized controller is presented.In this case, the parameters are selected as: Figure 5 depicts the 3D controlled (real) versus the desired trajectories of the agents.It is obvious from the figure that the method is successful to track the desired paths.Trajectory following errors tends to zero in a short time as depicted in Figure 6.In Figure 7 the result of formation tracking control and the formation shape of the agents are depicted.The controlled and desired Euler angles are also depicted in Figure 8 which also shows the accuracy of the method in controlling the Euler angles.Figure 9 shows the input signals.It is obvious that u 1i , i = 1, 2, 3 is less than 15N  after transient time and u 2i , u 3i , u 4i , i = 1, 2, 3 converges to 0N.m.

2) DECENTRALIZED FORMATION CONTROL
Now, the results for the proposed decentralized controller is presented.In this case, the parameters are selected as: Figure 10 depicts 3D controlled versus desired trajectories of the agents for the first 30sec.However, according to Figure 11, the z axis errors do not converges to zero.However, as it can be seen in Figure 12 the errors converges asymptotically to zero as time grows to 1000sec.It is due to the asymptotic behavior of the consensus algorithms.The behavior of agent 2 is different from two other ones as it is connected to two other agents at same time and this causes the errors converges to zero significantly faster.Figure 13 depicts the 3D controlled versus desired trajectories of the agents for the times between 970sec and 1000sec which shows a more precise tracking.Moreover, the result of formation tracking control of the agents for 30sec is depicted in Figure 14.The controlled and desired Euler angles are also depicted in Figure 15 for the first 100 sec which shows a precise tracking of the Euler angles.Figure 16 shows the input signals also for the first 100sec.It is obvious that like the centralized case in the decentralized one also u 1i , i = 1, 2, 3 remains less than 15N after transient time and u 2i , u 3i , u 4i , i = 1, 2, 3 converges to 0N.m. Figure 17 depicts distances between the agents.Although it has not been proved that the proposed controller is collision free, this can be investigated through simulation as depicted in Figure 17.

3) COMPARISON STUDY
In this section, the proposed method of this paper is compared with the one presented in [20] where the finite-time formation 142926 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.tracking control is proposed in a leader-follower structure.In this paper, a special structure is considered for the topology  of the agents and the leader in which it is assumed that the position of the leader should be in the average of the positions of the agents.Accordingly, we consider agent 1 as the leader to simulate the proposed method of [20].Physical parameters (moments of inertia, weights and l i ), initial values and desired trajectories are considered similar to our paper.Figure 18 depicts 3D controlled trajectories of agent 2 and agent 3 and the leader (agent 1 in our method).The errors are shown in Figure 19.It can be concluded from Figure 18 and Figure 19 that the method is successful in the formation tracking control in finite-time.However, comparing Figure 19 with Figure 11, it is obvious that although the convergence speed in the method of [20] is less than our decentralized method (and not the centralized one), in its transient time, it suffers from severe under and overshoots which is not appropriate.Similarly, according to Figure 20, apart from its transient response the method of [20] is successful in tracking of the Euler angles.However, it can be seen from Figure 21    3D formation tracking controls and formation shapes for both centralized and decentralized methods are depicted in Figure 23 and Figure 25, respectively.

VII. CONCLUSION
The problem of formation tracking control in multi-quadcopter systems was studied in this paper.Toward this, firstly the problem was solved for a centralized case in which predefined accuracies in the trajectory tracking and formation control were satisfied through definition of logarithmic BLFs in a backstepping structure.Due to the under-actuated nature of the quadcopter system, at first, a centralized formation control was designed for the altitude subsystem and then through definition of virtual inputs, a centralized controller was designed translational subsystem.Then, the centralized controllers were implemented in a decentralized way using high pass consensus filters.Each agent, then employed the virtual inputs to generate the desired roll and pitch trajectories were the attitude control system to track the desired trajectories with a predefined accuracy were designed using logarithmic BLFs in a backstepping structure.Finally, the efficiency of the proposed method were demonstrated through simulation both centralized and decentralized cases for two different spiral and S-shape scenario with line and triangle formation shpes, respectively.As expected, it was shown in the simulation that the results of the decentralized approach asymptotically converges to the centralized one due to the use of asymptotic convergent consensus filters.The proposed methods of this paper (both centralized and decentralized ones) were compared with the proposed method of [20] from different aspects of formation tracking precision both in transient and steady state times as well as control effort aspect.
Although it was shown in the simulation that the proposed method is collision free, the extension of the work to the collision free one with mathematical proofs is suggested as the future work.It is also proposed to extend the work to case in which the agents' controllers are derived in a distributed way at the beginning to omit the embedded consensus filters and make the method faster with less computational complexities.The proposed method can also be extended in several ways of multi-quadcopter systems with directional graph topology or for the systems in the presence of communication constraints such as protocol based or event triggered ones.

FIGURE 1 .
FIGURE 1.A group of quadcopters communicating through a wireless communication network.

FIGURE 3 .
FIGURE 3. General framework of the proposed controller for i -th agent.

FIGURE 4 .
FIGURE 4. Interconnection between agents in the simulation.

FIGURE 5 .
FIGURE 5. 3D controlled versus desired trajectories of the agents corresponding to the centralized case in the first scenario.

FIGURE 6 .
FIGURE 6. Position errors of the agents corresponding to the centralized case in the first scenario.Matlab-Simulink simulation.Three interconnected quadcopters are considered where the interconnection of the agents are depicted in Figure4.The physical parameters of the quadcopters are assumed asI xxi = 1.152 × 10 −2 kgm 2 , I yyi = 1.152 × 10 −2 kgm 2 , I zzi = 2.18 × 10 −2 kgm 2 , I ri = 5.2 × 10 −5 kgm 2 , l i = 0.28m, m i = 1.47kg.Two different scenarios are considered to evaluate the proposed method.In the first one, helical paths in which the agents form a line is considered in which the reference trajectories for quadcopters are formulated as follows:

FIGURE 7 .
FIGURE 7. 3D formation tracking control and formation shape corresponding to the centralized case in the first scenario.

FIGURE 8 .
FIGURE 8. Controlled versus desired Euler angles corresponding to the centralized case in the first scenario.

FIGURE 9 .
FIGURE 9. Control inputs corresponding to the centralized case in the first scenario.

FIGURE 10 .
FIGURE 10. 3D controlled versus desired trajectories of the agents for 0 to 30 sec corresponding to the decentralized case in the first scenario.

FIGURE 11 .
FIGURE 11.Position errors of the agents for 0 to 100 sec corresponding to the decentralized case in the first scenario.

FIGURE 12 .
FIGURE 12. Position errors of the agents for 0 to 1000 sec corresponding to the decentralized case in the first scenario.

FIGURE 13 .
FIGURE 13. 3D controlled versus desired trajectories of the agenst for 970 to 1000 sec corresponding the decentralized case in the first scenario.

FIGURE 14 .
FIGURE 14. 3D formation tracking control and formation shape of the agents corresponding to the decentralized case in the first scenario.

FIGURE 15 .
FIGURE 15.Controlled versus desired euler angles for 0 to 100 sec corresponding to the decentralized case in the first scenario.

FIGURE 16 .
FIGURE 16.Control inputs for 0 to 100sec corresponding to the decentralized case in the first scenario.

FIGURE 17 .
FIGURE 17. Distances between the agents corresponding to the decentralized case in the first scenario.

FIGURE 18 .
FIGURE 18.3D controlled trajectories of the agents and the leader corresponding to the method of[20] in the first scenario.

FIGURE 23 .
FIGURE 23. 3D formation tracking control and formation shape corresponding to the centralized case in the second scenario.

FIGURE 24 .
FIGURE 24.3D controlled versus desired trajectories of the agents corresponding to the decentralized case in the second scenario.

FIGURE 25 .
FIGURE 25. 3D formation tracking control and formation shape of the agents corresponding to the decentralized case in the second scenario.
state of q ia b1 (t).Then simply, the output of the consensus system, l ia b1 (t) is multiplied by N .In other words, a b1i (t) = k 2 b1 −Nl ia b1 (t) which asymptotically converges to a b1 (t).

TABLE 1 .
Pseudo code corresponding to distributed computation of a b1 , a e1 and α 2i related to the altitude control.

TABLE 2 .
[20]numerical comparison of our proposed method (centralized and decentralized ones) with the method of[20].