Decentralized H PID Team Formation Tracking Control of Large-Scale Quadrotor UAVs Under External Disturbance and Vortex Coupling

For more practical applications, a simple decentralized <inline-formula> <tex-math notation="LaTeX">$H_{\infty }$ </tex-math></inline-formula> proportional-integral-derivative (PID) team formation tracking strategy is proposed in this study for large-scale stochastic quadrotor unmanned aerial vehicles (UAVs) under external disturbance, intrinsic stochastic fluctuation and trailing vortex coupling. By adopting the virtual leader concept, the reference trajectory of each UAV is generated by the combination of virtual leader trajectory with a specific time-varying formation offset to form a desired flight formation shape. Then, by using the proposed optimal <inline-formula> <tex-math notation="LaTeX">$\text{H}_{\infty }$ </tex-math></inline-formula> decentralized PID controller, each UAV can efficiently attenuate the effect of external disturbance and trailing vortex coupling from the neighboring quadrotor UAVs on the team formation reference tracking performance simultaneously. To avoid solving a complex nonlinear Hamiltion Jacobi inequality (HJI) for PID control of each UAV, the Takagi-Sugeno (T-S) fuzzy method is employed to interpolate several local linearized UAVs to approximate the nonlinear stochastic quadrotor UAV system so that the HJI can be transformed to a set of bilinear matrix inequalities (BMIs). By using novel variable transformation method, the optimal decentralized <inline-formula> <tex-math notation="LaTeX">$\text{H}_{\infty }$ </tex-math></inline-formula> PID team formation tracking control problem of large-scale UAVs can be independently designed in terms of a set of independent LMIs-constrained optimization problems for each quadrotor UAV. Further, to simplify the design procedure in a single run, the large amount of decentralized <inline-formula> <tex-math notation="LaTeX">$\text{H}_{\infty }$ </tex-math></inline-formula> PID controller parameters for each quadrotor UAV in the team formation can be efficiently obtained via the current convex optimization technique without any conventional parameter tuning procedure for PID design. Finally, a simulation example of team formation reference tracking control for 25 quadrotor UAVs with external disturbance and vortex coupling is given to validate the team formation tracking performance of the proposed decentralized <inline-formula> <tex-math notation="LaTeX">$\text{H}_{\infty }$ </tex-math></inline-formula> PID tracking control method in comparison with conventional decentralized <inline-formula> <tex-math notation="LaTeX">$H_{\infty }$ </tex-math></inline-formula> T-S fuzzy tracking control scheme.


I. INTRODUCTION
Recently, unmanned aerial vehicles (UAVs) have attracted a growing interest in academic and industrial field due to wide and high maneuverability [1], [2]. Since the quadrotor UAV is a highly nonlinear dynamic system with 6 degrees of freedom The associate editor coordinating the review of this manuscript and approving it for publication was Zhiguang Feng . and 4 actuators (i.e., the underactuated characteristics), it is not easy to control quadrotor UAV for a desired trajectory expecially by a conventional PID control. Nonlinear adaptive control [3], [4] and nonlinear robust control [5], [6], [7], [8] have been developed for UAV. However, a single quadrotor UAV has many limitations on applications, for example, battery power and loading capacity. Unlike a single quadrotor UAV, a team formation of multiple UAVs can achieve more VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ challenging and complex tasks [9]. Consequently, a team formation of multiple quadrotor UAVs is becoming practical and popular in many application areas, i.e., environmental monitoring, exploration of unsecured regions, detection, and interception, target searching and rescuing, battlefield surveillance, etc [9]. Since the missions of quadrotor UAV are more challenging and complex, several formation control strategies have been developed, i.e., leader-follower method [10], virtual leader method [11], etc. Leader-follower method is that one of quadrotor UAVs in the team plays the role of leader to represent the behavior trajectory of the team and is controlled to track a desired reference path, while other quadrotor UAVs are asked to track the trajectory of quadrotor UAV with a formation shape [10]. Even leader-follower method was developed more earlier and commonly used in the team formation of UAVs, there exist some drawbacks. One is that the team formation is always dependent on the leader UAV. Thus, if the leader UAV is failed or disintegrated, the whole team formation can not be maintained or will collapse. When compared with leader-follower method, the virtual leader method is more robust and can avoid the above problems, i.e., virtual leader is a reference model at the central of the desired shape to generate a reference trajectory with a desired formation shape [11], [12].
Since the virtual leader is a reference model and will not suffer from perturbation, it can efficiently reduce the probability of team formation failure in reality. Recently, the formation control strategy design of UAV has been widely investigated. The robust H ∞ team formation network control is designed in [11] in consideration of internal fluctuation in UAV system. Also, to save the communication resource in UAV network, the robust event-triggered team formation tracking is developed in [12]. Further, to reduce the interconnected effect in UAV network, the decentralized H ∞ tracking control of a large-scale team formation UAV network system is proposed to passively attenuate the interconnceted effect [13]. Therein, to relieve the design difficulty of solving a set of nonlinear Hamilton Jacobi inequalities (HJIs), the Takagi-sugeno (T-S) fuzzy system is employed to interpolate several local linearized systems to approximate nonlinear quadrotor UAV. In this case, the nonlinear HJIs for the H ∞ tracking control design can be reformulated as a set of coupled linear matrix inequalities (LMIs), which can be easily solved by the current convex optimization techniques [11], [12], [13]. On the other hand, based on a specified potential function associated with network topology of UAV network, a decentralized multi-subgroup formation algorithm is developed to achieve the desired formation control and collision avoidance [14].
In the aforementioned articles, it is worth to point out that the computational complexity of proposed strategies in [11], [12], [13], and [14] is still large. More specifically, a large computational time is needed to calculate the fuzzy control law u(t) = L j=1 h j (z(t))K j (x(t) − r(t)) for the centralized fuzzy control in [11], [12] and u i (t) = L j=1 h i,j (z i,j (t))K i,j (x i (t) − r i (t)), for i = 1, . . . , N for the decentralized tracking control of team formation tracking of N UAVs [13], [15], respectively. More efforts are still needed for these fuzzy team formation tracking controls of large-scale UAVs for practical applications. On the other hand, to achieve an accurate formation control, the potential function in [14] for the controller synthesis is with high nonlinearity. As a result, to remedy the complex T-S fuzzy team formation tracking control design or other nonlinear control strategy design for large-scale UAVs, the PID tracking control is probably the most simple method to overcome the complex and highly nonlinear teamformation tracking control design problem.
In this study, a robust decentralized H ∞ PID team formation tracking control design is proposed for large-scale UAVs with a desired time-varying formation shape under wind disturbance and coupling of trailing vortex from neighboring UAVs. The proposed decentralized H ∞ PID team formation tracking control design for each UAV can efficiently attenuate the external disturbance and decouple the coupling effect of trailing vortex from the neighboring UAVs. Therefore, it can significantly reduce the control design complexity and reduce the computational complexity of practical team formation tracking of large-scale quadrotor UAVs. At the beginning, the proposed robust decentralized H ∞ PID team formation tracking control design of large-scale UAVs can be transformed to a set of independent nonlinear partial differential HJIs for each UAV in the team formation. Then, PID controller in each quadrotor UAV needs to be specified to satisfy the corresponding HJI. Further, the optimal decentralized H ∞ PID team formation tracking control design needs to solve an HJI-constrained optimization problem for PID control parameters of each quadrotor UAV in the team formation. However, at present, it is not easy to specify PID controller parameters to solve the corresponding HJIconstrained optimization problem for each quadrotor UAV to achieve the optimal decentralized H ∞ PID team formation tracking control design of large-scale UAVs. To overcome the difficulty of solving nonlinear partial differential HJI, the T-S Fuzzy interpolation scheme is employed to interpolate multiple local linearized stochastic systems to approximate the nonlinear stochastic system of UAV. Therefore, each HJIconstrained optimization problem of the optimal decentralized H ∞ PID team formation tracking control design of each quadrotor UAV in the team formation can be transformed to a corresponding Riccati-like inequalities (RLIs)-constrained optimization problem. Therefore, the HJI-constrained optimization problem of the optimal decentralized H ∞ robust PID team formation tracking control design of large-scale quadrotor UAVs can be transformed to an independent LMIs-constrained optimization problem for each quadrotor UAV in the team formation. Since the decentralized H ∞ robust PID team formation control can be independently designed, the number of UAVs in the teamformation can be increased to a very large scale. A simulation example of time-varying team formation task for 25 UAV system is provided to validate the performance of proposed decentralized H ∞ PID control strategy in comparison with the conventional decentralized H ∞ T-S fuzzy team formation tracking control strategy in [15].
The main contributions of this work are describe as follows: 1. The previous PID control designs of UAV are always dependent on the linearized model of UAV and there are little PID control designs for the team formation tracking control of UAVs. In this study, we consider the robust decentralized H ∞ PID team formation tracking control design problem of large-scale nonlinear stochastic UAVs with external disturbance and coupling of trailing vortex from neighboring UAVs. The proposed robust decentralized H ∞ PID team formation tracking control of large-scale UAVs can significantly simplify the design complexity and save much computation time of PID controller than the conventional H ∞ fuzzy team formation control of large-scale UAVs in [13], [14], [15], and [16]. Therefore, the proposed control scheme has great potential to practical application for large-scale time-varying team formation control of UAVs.
2. Based on T-S fuzzy interpolation scheme to interpolate multiple local linear stochastic systems to approximate nonlinear stochastic UAV, the proposed optimal decentralized H ∞ PID team formation tracking control design problem of UAVs can be transformed to an independent LMIsconstrained optimization problem for each UAV in the team formation. Therefore, the number of UAVs in the team formation can be easily extended to a very large scale without increasing the design difficulty.
3. When compared with other conventional PID control designs for the linearized quadrotor UAV system, the uncertain wind disturbance, intrinsic random Wiener fluctuation, and wind coupling of trailing vortex from neighborhood UAVs can be efficiently attenuated by the proposed optimal decentralized H ∞ PID team formation tracking control strategy of large-scale nonlinear stochastic UAVs. Further, unlike the conventional complex tuning process of PID control parameters, the proposed optimal H ∞ decentralized control parameters can be found in a single round of convex optimization algorithm with help of LMI toolbox in MATLAB.
The remainder of the study is organized as follows. In Section II, the large-scale quadrotor UAV systems are introduced and the problem formulation of the decentralized robust stochastic H ∞ PID control for the team formation quadrotor UAVs is presented. In Section III, the design condition is transformed into an equivalent HJI problem for each UAV. In Section IV, with the help of T-S fuzzy method, the HJI problem of each UAV is transformed into a set of bilinear matrix inequalities (BMIs). Furthermore, with the help of the novel variable transformation method, a set of independent BMIs are transformed into a set of independent LMIs for each UAV and are solved by using the LMI toolbox in MATLAB. In Section V, the decentralized robust stochastic H ∞ PID team formation tracking control method for 25 quadrotor UAVs in comparison with the traditional decentralized H ∞ T-S fuzzy team formation tracking control method is given to illustrate the effectiveness of the proposed method. Finally, the conclusion is made in Section VI.
Notation: A T denotes the transpose of matrix A; A ≥ 0 (A > 0) denotes symmetric semi-positive definite matrix (symmetric positive definite matrix); I n denotes the ndimensional identity matrix; E{·} is the expectation operator; x 2 denotes the Euclidean norm for the given vector x(t) is the finite mean energy; C 2 denotes the class of function f (x) which is twice continuously differentiable with respect to x; ∂f (x) ∂x denotes the gradient column vector of differentiable function f (x) ∈ R; ∂ 2 f (x) ∂x 2 is the Hessian matrix of twicedifferentiable function f (x) ∈ R.

II. PRELIMINARY AND PROBLEM FORMULATION A. PHYSICAL PLANT OF QUADROTOR UAVs
In this study, we simultaneously consider the three position variables with the corresponding three velocities of the quadrotor UAV and three attitude states with the corresponding three angular velocities of the quadrotor UAV as the state variables in the state space dynamic model of the quadrotor UAV. To model the quadrotor UAV dynamic, two reference frames including inertial frame and body frame are introduced, and the free body diagram of the quadrotror UAV is shown in Fig. 1. The position of the mass center of the quadrotor UAV is represented as coordinates (x, y, z) w.r.t. an inertial frame associated with the unit vector basis (e x , e y , e z ). On the other hand, the attitude of the quadrotor UAV is denoted by the three Euler angles (φ, θ, ψ) which define the orientation vector of the quadrotor UAV w.r.t. a body frame associated with the unit vector basis (e bx , e by , e bz ). Besides, φ is a roll angle (− π 2 < φ < π 2 ), θ is a pitch angle (− π 2 < θ < π 2 ), and ψ is a yaw angle (−π < ψ < π ). With the Newton Euler method, the ith quadrotor UAV dynamic model of the large-scale UAVs networked system can be described as the following state space model [12]: for i = 1, 2, . . . , N where x i,1 (t), y i,1 (t), z i,1 (t) ∈ R denote the three position states of the ith quadrotor UAV in the inertial frame, x i,2 (t), y i,2 (t), z i,2 (t) ∈ R denote the three velocity states of the ith quadrotor UAV in the inertial frame, φ i,1 (t), θ i,1 (t), ψ i,1 (t) ∈ R denote the three attitude states of the ith quadrotor UAV in the body frame, φ i,2 (t), θ i,2 (t), ψ i,2 (t) ∈ R denote the angular velocity states of the ith quadrotor UAV in the body frame, F i (t) ∈ R denotes the total thrust of the ith quadrotor UAV, , v i,ψ (t) represent the external disturbances of the quadrotor UAV, m i ∈ R + denotes the total mass of the quadrotor UAV, g ∈ R + denotes the gravitational acceleration, J i,φ , J i,θ , J i,ψ ∈ R + denote the moments of the inertia of φ i (t), θ i (t), and ψ i (t), respectively, and d i, represent aerodynamic damping coefficients of the quadrotor UAV.
After introducing a single dynamic system of quadrotor UAV, the large-scale team formation for the dynamic system of quadrotor UAVs is denoted as a set S = {S 1 , S 2 , . . . , S N }, which is composed by N quadrotor UAVs. From the practical point of view, each quadrotor UAV is affected by the external disturbances, coupling effects, internal continuous random fluctuations, and internal discontinuous random fluctuations. To make the dynamic system of quadrotor UAVs more practical, the Wiener process and Poisson counting process are considered to formulate the internal continuous and discontinuous random fluctuations due to parameters fluctuations in each quadrotor UAV system. Therefore, the ith nonlinear stochastic quadrotor UAV dynamic system S i in (1) can be reformulated as follows: with ∈ R 12 represents the ith dynamic system state of the quadrotor UAV, u i (t) ∈ R 4 represents the control input of the ith quadrotor UAV, X j (t − τ i,j (t)) ∈ R 12 represents the trailing vortex coupling term from the jth UAV to the ith UAV with a time-varying represent the interconnected couplings as the trailing vortex effect caused by the neighborhood UAVs of the ith UAV, N i represents a neighborhood set of quadrotor UAVs nearing the ith UAV, v i (t) ∈ L 2 [0; ∞) denotes the external disturbance of the ith UAV, w i (t) ∈ R denotes the Wiener process of the ith quadrotor UAV which is continuous but non-differentiable, σ i (X (t))dw i (t) ∈ R 12 denotes the effect of the continuous stochastic intrinsic fluctuation caused in the ith quadrotor UAV, p i (t) ∈ R denotes the Poisson counting process with jump intensity λ i , and E i (X i (t))dp i (t) ∈ R 12 is the discontinuous intrinsic fluctuation caused in the ith quadrotor UAV.
Remark 1: The Wiener processes {w i (t)} and Poisson counting processes {p i (t)} have the following properties [8], [11], [17], [18]: In order to let N UAVs maintain a specific time-varying formation shape during the flighting, the virtual leader method is used in this study. The team formation shape of N UAVs is constructed by one virtual leader and N UAVs and then each quadrotor UAV will track the virtual leader's trajectory with a specific time-varying formation shape to achieve the formation task. The reference model of the desired trajectory for N quadrotor UAVs of the team formation can be defined as follows: where denotes the desired trajectory of the ith quadrotor UAV to be tracked. X i,r1 (t) ∈ R 6 denotes the three positions and attitudes of the ith quadrotor UAV and X i,r2 (t) ∈ R 6 denotes the three velocities and angular velocities of the ith quadrotor UAV. r(t) ∈ R 12 is the trajectory of the virtual leader to be generated as the desired flight trajectory of the team formation of N quadrotor UAVs and d i (t) ∈ R 12 denotes the time-varying offset from the virtual leader to the ith quadrotor UAV. A i,r ∈ R 12×12 is the asymptotically stable matrix and nonsingular. C i,r ∈ R 12×12 is the reference input matrix.
Remark 2: (i) The ith reference model in (3) is to generate the reference flight trajectory X i,r (t) of the ith quadrotror UAV. The system matrix A i,r is to be specified as an asymptotically matrix. If we set C i,r = −A i,r . Then, at the steady state, the trajectory of the ith reference model in (3) . Therefore, if we can control N UAVs to track the corresponding N trajectories of reference models X i,r (t) i = 1, . . . , N in (3), then the N UAVs will form a timevarying formation shape (d 1 (t), . . . , d i (t), . . . , d N (t)) around the trajectory of the virtual leader r(t) at the steady state. This phenomenon will be confirmed by the simulation example in the sequel. (ii) In the conventional formation control for UAVs in [11] and [12], the formation control is achieved by specifying in the numerotor of the H ∞ team formation tracking strategy in (5) to simplify the team formation control design. However, it is difficult to apply to a time-varying formation Because the Proportional-Integral-Derivative (PID) controllers have been used in the most automatic process control applications in industry, the PID controller is employed to control each quadrotor UAV for the team formation in this study. The conventional PID control for the ith UAV is shown as follows: where ) ∈ R 6 denotes the tracking error of the three positions and three attitudes of the ith quadrotor UAV, In the realistic situation, the effect of external disturbances v i (t) and trailing vortex is inevitable during the team formation tracking process of quadrotor UAVs and it may deteriorate the team formation tracking performance. Besides, the reference input r(t) and the team formation offset d i (t) in (3) of the ith quadrotor UAV are unpredictable for the stochastic quadrotor UAV system in (2). Thus, the decentralized H ∞ PID team formation tracking strategy is proposed for each UAV in the team formation to efficiently attenuate the effect of the external disturbances, trailing vortex and reference input on the team formation tracking performance. Then, the decentralized robust stochastic H ∞ PID team formation tracking control strategy of the ith quadrotor UAV is proposed as follows: , respectively, and V i (X i (0), X i,r (0)) denotes the effect of the initial condition, which is to be deducted in the design procedure. By using the decentralized robust stochastic H ∞ PID team formation tracking control strategy to boost the team formation performance of in the ith quadrotor UAV, the aim is to design the specific PID tracking controller u i (t) to attenuate the worst-case effect of external disturbance, reference input, and trailing vortex on the team formation tracking performance below a prescribed level ρ i in (5) for each quadrotor UAV in the team. However, the conventional PID control methods are limited to the linear dynamic systems or simple nonlinear dynamic systems so the methods should be improved to solve the complex decentralized H ∞ team formation control problem of highly nonlinear dynamic UAV system. In this study, we employ the method of coordinate transformation so that the conventional PID control problems can be systematically analyzed and conveniently designed for the decentralized H ∞ team formation tracking control of large-scale UAVs. The PID controller u i (t) in (4) is reformulated as follows: with

III. DECENTRALIZED ROBUST STOCHASTIC H ∞ PID TRACKING CONTROL FOR THE NONLINEAR LARGE-SCALE TEAM FORMATION QUADROTOR UAVs
In this section, the decentralized robust stochastic H ∞ tracking control design for the team formation of the ith quadrotor UAV is addressed. Let us augment the ith stochastic quadrotor UAV system in (2), the reference system in (3) and the PID controller in (6), then we get the augmented system of ith quadrotor UAV as follows: whereX The system matrices are given as: ] T Based on the augmented system in (7), The decentralized robust stochastic H ∞ PID team formation tracking control strategy of the ith quadrotor UAV can be reformulated as: (0)). Clearly, with the help of the augmented system of the ith UAV in (7), the robust decentralized H ∞ PID team formation tracking control design problem in (5) is reformulated as the decentralized H ∞ PID stabilization control problem of the augmented system for each quadrotor UAV in (8) to simplify the design procedure. However, the conventional method can not be directly used in the stochastic nonlinear augmented system in (7) because the property of the Wiener process and Poisson counting process in the nonlinear ith quadrotor system is indifferentiable at every time. Fortunately, by using the Itô-Lévy formula in [22], the decentralized robust stochastic H ∞ PID control problem for the team formation of quadrotor UAVs can be transformed into a set of independent HJIs. Before we deduce the HJIs for the decentralized H ∞ PID control in Theorem 1, the Itô-Lévy formula and technical inequality are described as follows: Lemma 1: Define the Lyapunov functionV i (X i (t)) ∈ C 2 withV i (0) = 0 andV i (X i (t)) > 0, then the Itô-Lévy formula ofV i (X i (t)) for the stochastic nonlinear augmented system of the ith UAV in (7) is given as follows [22]: Lemma 2: For any two matrices G and H with appropriate dimensions, the following matrix inequality holds [23]: where Z is any positive definite symmetric matrix.
With the help of above lemmas, the decentralized robust stochastic H ∞ PID reference tracking control design problem for the team formation of quadrotor UAVs is transformed to an equivalent set of indenpendent HJIs in the following theorem.
Theorem 1: In the decentralized stochastic nonlinear H ∞ PID team formation tracking control of the ith quadrotor UAV in (7), if there exists a PID control gain (6) such that the following independent HJI i has a positive solutionV i (X i (t)) > 0 withV i (0) = 0 for each UAV then (i) the decentralized robust stochastic H ∞ PID team formation tracking performance of quadrotor UAVs in (8) is guaranteed with a prescribed disturbance level ρ i for all possiblev i (t) ∈ L 2 [0; ∞) from the perspective of the energy.
(ii) the mean square asymptotical team formation tracking ability is also achieved i.e., Proof: Please refer to Appendix A. Since the HJI i in (11) for the decentralized H ∞ PID team formation tracking control of the each quadrotor UAV does not contain the information of the other quadrotor UAVs in the formation team, the (11 ) can be solved independently to obtain PID control parameters of each quadrotor UAV and therefore, the number N of UAVs in the team formation can be increased to a very large scale. Furthermore, the optimal decentralized H ∞ PID team formation reference tracking control design of N quadrotor UAVs can also be formulated as the following N independent HJI i constrained optimization problems.
Due to nonlinear termsF i (X i (t),Ḡ(X i (t), K i ),D i (X i (t)), at present, there exists no analytical or numerical method to efficiently solve N HJI i in (11) or N HJI i -constrained optimization problems in (12) for the optimal decentralized H ∞ PID team formation tracking control of N UAVs.

IV. DECENTRALIZED H ∞ ROBUST PID CONTROL DESIGN FOR THE TEAM FORMATION OF LAREGE-SCALE QUADROTOR UAVs VIA T-S FUZZY APPROACH
In the previous section, the decentralized H ∞ PID team formation tracking control strategy of N stochastic nonlinear quadrotor UAVs in (8) is transformed into the N equivalent HJI i problems in (11). However, it is still very difficult to find an analytical or numerical method to solve N nonlinear partial differential HJIs in (11) or N HJI i -constrained optimization problems in (12) for the decentralized H ∞ PID team formation tracking control design of N stochastic nonlinear quadrotor UAVs. Hence, the T-S fuzzy interpolation method is employed to interpolate several local linear stochastic systems to approximate the ith nonlinear stochastic quadrotor UAV system in (2). First, the kth rule of the T-S fuzzy model for the ith stochastic nonlinear quadrotor UAV dynamic system in (2) can be described as follows [19], [20]: Plant rule k: where L is the number of the fuzzy rules, are the local matrices of f i (X i (t)), g i (X i (t)), σ i (X i (t)), and E(X i (t)) respectively. (F i,1,k , . . . , F i,i−1,k , F i,i+1,k , . . . , F i,N ,k ) are the local matrices of the effectively trailing vortex from other UAVs. p i,l (t) for l = 1, . . . , g are the premise variables about the ith UAV state, andF k,1 , . . . ,F k,g are the fuzzy sets and g is the number of the premise variables of each UAV. The kth membership function ν k (p i (t)) is defined as follows: whereF k,l (p i,l (t)) is the membership grade ofp i,l (t) inF k,l . p i (t) = [p i,1 (t), . . . ,p i,g (t)] T , and L k=1 ν k (p i (t)) > 0. The interpolation functions of the ith quadrotor UAV system can be inferred as follows: ≥ 0 for k = 1, 2, . . . , L, and it is clear that By the defuzzification process in [20], the nonlinear system of the ith UAV can be represented: VOLUME 10, 2022 i.e., the nonlinear stochastic system of the ith UAV in (2) can be represented by interpolating L local linearized stochastic systems through L fuzzy interpolation functions {h k (p i (t))} L k=1 . By substituting the T-S fuzzy model (16) of the ith UAV in (2) into the augmented system in (7), we have: With the help of the T-S fuzzy interpolation, the complex HJI i -constrained optimization in (12) for the decentralized robust H ∞ PID control of N quadrotor UAV systems will be transformed into an equivalent LMIs-constraint optimization problem for each quadrotor UAV, which can be easily solved by using the LMI toolbox in MATLAB. Before the discussion of the main theorem, a quadratic Lyapunov function is selected for HJI i in (11) to simplify the design procedure of the decentralized H ∞ PID control of team formation of N UAVs.
To begin with, the Lyapunov function of the ith augmented quadrotor UAV system in (17) is defined as follows: whereP i ∈ R 30×30 is a symmetric positive-definite matrix. Theorem 2: For the decentralized H ∞ PID team formation tracking control of the ith quadrotor UAV stochastic system in (17), if we can specify the PID control gains K i = [K i,P , K i,D , K i,I ] in (6) and symmetric positive-definite matrix P i such that the following independent RLIs hold: (8) is achieved with a prescribed disturbance level ρ i > 0 for all possiblev i (t) ∈ L 2 [0; ∞) from the perspective of the energy. (ii) the mean square asymptotical team formation tracking ability is also achieved i.e., X i (t) → X i,r (t) as t f → ∞ ifv i (t) ∈ L 2 [0; ∞).

then (i) the decentralized robust H ∞ PID team formation tracking control strategy of N UAVs in
Proof: Please refer to Appendix B. Remark 3: If PID controller of the ith UAV in (6) can satisfy with RLIs, i = 1, . . . , N , k = 1, . . . , L, in (19 ), i,e., N local fuzy linearized systems are took care by robust H ∞ PID controller, then the PID controller in (6) can guarantee the H ∞ team formation tracking performance in (5) for L local fuzzy linearized systems of ith nonlinear stochastic UAV. Further, N RLIS are decoupled for N UAVs, i,e., they can be solved independently without the information of other UAVs.
If the sufficient condition in (19) is satisfied, then the decentralized robust stochastic H ∞ PID tracking control performance for the team formation of N quadrotor UAVs in (8) is guaranteed with a prescribed disturbance level ρ i . However, the design parametersP i and K i are coupled and bilinear matrix inequalities in (19) can not be solved by the current convex optimization method. To deal with the problem in (19), the following transformation method is adopted.
By premultiplying and postmultiplyingP −1 i =W i , the RLIs in (19) are equivalent tō 3 ∈ R 6×6 , then we have By using the variable transformation technique Therefore, the RLI in (20) can be formulated as follows: By using Schur complement [21] to (24), we have the following LMIs: Then, the above analyses can be concluded as the following main result.
Theorem 3: If we can solveW i > 0 andȲ i > 0 from LMIs in (25) for a prescribed ρ i , i = 1, . . . , N , and k = 1, . . . , L, then the decentralized H ∞ PID team formation tracking controllers of N UAVs in (6) are designed as follows: By solving LMIs in (25) with the help of the LMI toolbox in Matlab and based on Theorem 3, the PID controllers for the decentralized H ∞ PID team formation tracking control of N stochastic quadrotor UAVs are independently designed as (26). However, in practical applications, the control input is always restricted by physical mechanism and bounded by the saturation of actuator. To meet the real situation, the saturation control method in [24] will be used as follows: Suppose thatX i (t) is restricted to stay in an invariant ellip- iȲ i holds, then u i (t) 2 ≤ ν i for all t. Therefore, by using Schur complement [21], if the following LMIs hold, then we have the result, Thus, the design condition of the decentralized H ∞ PID control for the team formation of stochastic quadrotor UAVs is also satisfied with the saturation property of actuator if the conditions in (25), (27) are achieved simultaneously. Therefore, the optimal decentralized H ∞ saturation PID control design for the team formation of stochastic quadrotor UAVs can be designed by the following LMIs-constrained optimization problem for each UAV The LMIs-constrained optimization problem in (29 ) can be solved by decreasing ρ i until there exists no feasible solution W i > 0 in LMIs.

Remark 4: The design complexity of proposed H ∞ PID team formation tracking controller lies mainly in solving
LMIs in (25), which are solved based on Newton's method. Therefore, based on [21], the computational complexity is about O(NLn(n + 1)/2), where N is the number of quadrotor UAVs in the team formation, L is the number of fuzzy rules, and n is the dimension ofP i orW i , n = 3 × 12 = 36 in the quadrotor UAV case.
Remark 5: Since quadrotor UAV is a highly nonlinear system as shown in (2), the H ∞ stabilization or tracking control design problem of quadrotor UAV needs to solve a corresponding HJI. In general, there still exists no analytic or numerical method to efficiently solve the HJI for PID control, except a very special system. To avoid this difficulty in solving HJI, the conventional PID control design only considers a local linearized system of quadrotor UAV of some operation point, so that the HJI in (11) becomes a Riccatilike inequality (RLI) in (19), which can be transformed to an equivalent LMI in (25). However, the quadrotor UAV havdto shift to many operation points during flying, for example, in taking off, landing and hovering. Therefore, the conventional PID control designs need to switch to an adequate PID controller as the operation point of quadrotor UAV shifts during flying. In the proposed design method, the T-S fuzzy interpolates L local linearized systems at L operation points by fuzzy interpolation method in (17) (in the simulation example L = 72). Therefore, the proposed PID controller can achieve the decentralized robust H ∞ tracking control performance at L local linearized systems (i.e., L operation points) of UAVs and has the ability to reject the worst-case effect of external disturbance, vortex coupling and saturation of actuator. However, the proposed decentralized H ∞ PID team formation control design needs to solve L LMIs in (25) or the LMIs-constrained optimization problem in (29) for each UAV in the team formation tracking control design.

V. SIMULATION
To illustrate the design procedure and validate the effectiveness of the proposed decentralized H ∞ PID team formation stochastic tracking control problem of N stochastic nonlinear quadrotor UAVs, a team formation tracking control of 25 UAVs is given as a design example in this section. To illustrate the implementation simplicity and validate the effectiveness of the decentralized H ∞ PID team formation tracking control, the comparison with conventional decentralized H ∞ T-S fuzzy control for team formation tracking of 25 UAVs in [15] will be showed. Since the decentralized controller is equipped with the individual information of each UAV and does not need the information of other UAVs, the number of the UAVs for the team formation can be easily increased to very large scale.

A. PARAMETERS SETTING OF LARGE-SCALE UAVs FOR TEAM FORMATION
The parameters of each stochastic dynamic of the UAV system in (2) for the 25 UAVs are identical and can be given as follows. m i = 1kg is the mass of each UAV. J i,X = J i,Y = 1.25Ns 2 /rad and J i,Z = 2.2Ns 2 /rad are the moments of the inertia of φ i , θ i , and ψ i , respectively. d i,x = d i,y = d i,z = 0.01Ns/m and d i,φ = d i,θ = d i,ψ = 0.013Ns/m are the aerodynamic coefficients of each quadrotor UAV. g = 9.8m/s 2 is the acceleration of gravity. f i,j (X i (t))X j (t − τ i,j (t)) = [0, 0, 0, 0, 0, 0, x i,2 (t)×0.01x j,2 (t −0.1), 0, z i,2 (t)× 0.01z i,2 (t − 0.1), 0, 0, 0] T , j ∈ N i denote that each UAV is suffered the trailing vortex coupling from neighboring UAVs in N i . External disturbance v i (t) for each UAV is set as zero-mean Gaussian noise with unit variance. 25 independent Wiener processes {w i (t)} 25 i=1 in Fig.14 and Poisson counting process{p i (t)} 25 i=1 in Fig.13 with Poisson counting intensity λ i = 1 are set in 25 UAVs of the team formation.
To construct the T-S fuzzy model of each UAV in (17 ) Based on the above setting, each UAV in the 25 UAVs has 32 fuzzy rules, respectively. The kth T-S fuzzy model of the ith UAV dynamic system in (16) is described as follows.
For the team formation tracking control problem of 25 quadrotor UAVs, in the design, the position and attitude tracking performance as well as velocity and angular velocity tracking performance are supposed to be the same importance. Further, in order to achieve a small optimal H ∞ PID team formation tracking performance ρ * i in (29), consequently, the weighting matrices are specified for each UAV subsystem as follows: In order to make 25 UAVs maintain a diamond time-varying formation shape, the trajectory of the ith UAV is constructed by a virtual leader r(t) with a specific timevarying distance d i (t) in (3) for i = 1, . . . , 25 and A i,r = −4I 12 , C i,r = 4I 12 in (3). However, the diamond formation shape is difficult to directly design so an indirect method is proposed. First, we design the square formation with time-varying distance d i (t) and then use the coordinate transformation to rotate d i (t) counterclockwise π The d i (t) is obtained by rotating d i (t) counterclockwise π 4 from the virtual leader r(t) in the y-z plane and the diamond formation shape is shown in Fig. 2 at t = 0. In the example, the trajectory of the virtual leader is specified with a radius 6 m in the x-y plane and is rising at a constant velocity at 2 m/s on the z-axis as shown in Fig. 3. The reference input r(t) in (3) with three positions [x r1 (t), y r1 (t), z r 1 (t)] T and three attitudes [φ r1 (t), θ r1 (t), ψ r1 (t)] T is given as follows [3]: On the other hand, the velocities of three positions and three attitudes are the differentiations of their three positions and three attitudes with time.

B. SIMULATION RESULTS
Based on the above parameters setting, the simulation results for the proposed robust decentralized H ∞ PID team formation tracking of 25 UAVs are shown in this section. Furthermore, the saturation constant ν i in (28) for each UAV is set as 100 (i.e., ν i = 100, i = 1, . . . , 25) [25]. By solving 25 LMIs constrained optimization problems in (29), the decentralized PID tracking controllers for the 25 UAVs in the team formation can be designed to achieve the optimal H ∞ decentralized PID team formation tracking performance under the external disturbance and interconnected coupling in the 25 UAVs. The average optimal attenuation level ρ * i of 25 UAVs is 0.92 when the ρ i is iteratively decreased from 1.5 until there exists no solutionW * i > 0 for the ith UAV in the LMIsconstrained optimization problem in (29). The corresponding optimal PID control parameters 1,2 are obtained by solving the LMIs-constrained optimization problem in (29).
The 3D flight trajectories of 25 UAVs in the team formation are shown in Fig. 3 with the time-varying formation shape within 50 seconds. It can be seen that each UAV can track the desired formation trajectory r(t) + d i (t) where the virtual trajectory r(t) is with radius 6m circle in the x-y plane and 2m/s ascent rate in the z-axis. To highlight the time-varying formation shape, the snapshot of 3D flight trajectories of 25 UAVs during 20-30 seconds is shown in Fig. 4. In Fig. 4, the formation shape for 25 UAVs team formation is a diamond shape and becomes bigger during 20-25 seconds and smaller during 25-30 seconds.
The three positions and the three velocities of 25 UAVs are shown in Figs. 5-7, respectively. Due to the time-varying formation shape, it can be seen that the trajectories of y-axis and z-axis of any two UAVs are derivated with a time-varying distance in Figs. 6-7. For example, for 20s-25s in Fig. 6, the distance between any two UAVs in y-axis at 20s is larger than the distance between any two UAVs at 25s. On the other hand, since the time-varying formation of 25 UAVs only affects the trajectories of y-position and z-position, 25 UAVs track the identical trajectory of x-position as shown in Fig. 5. From the Figs. 5-7, by the proposed robust PID formation tracking controller, the position tracking and the velocity tracking can be effectively achieved for 25 UAVs only with some fluctuation mainly caused by the stochastic effect due to the intrinsic Poisson counting process P i (t) in Fig.13 and Wiener process in Fig.14. Since the scales of velocities is x-axis, y-axis and z-axis are large than the scales of positions, the fluctuations are also amplified too. It is worth to point out  that the stochastic effects on UAV system are also efficiently reduced during the robust H ∞ PID tracking control process. For example, the velocity on y-axis can be controlled to reference velocity when the discontinuous Poisson counting jump occurs in Fig. 6.
The three attitude variables of 25 UAVs with the corresponding angular velocities are shown in Figs. 8-10, respectively. Caused by the different time-varying formation shapes for 25 UAVs, the desired attitudes to be tracked for each UAV are also different. Especially, due to the time-varying offset d i (t) of each UAV, it makes the desired attitudes   The control input u 1 (t) = [F 1 (t), τ 1,φ (t), τ 1,θ (t), τ 1,ψ (t)] T of the 1st UAV is shown in Fig. 11. Under the effect of intrinsic continuous Wiener fluctuation, and discontinuous Poisson fluctuation on UAV system, it can be seen that the control input has continuous fluctuation and some suddenly jumps. Indeed, once a Poisson jump occurs, the controller will use a relatively large control effort to reduce its effect. On the other hand, since the saturation effect is considered during the design, it is clear that the control input of 1st UAV is below the saturation constant, i.e., u 1 (t) ≤ ν 1 = 100, ∀t ∈ [0, 50].
In Fig. 12, in comparison with the 3D flight trajectory for 25 UAVs based on the proposed decentralized H ∞ PID team formation tracking control strategy, the team formation tracking performance of the decentralized robust H ∞ T-S fuzzy tracking control design [15] is almost equal. Based on the above simulation result in Fig.5-11, the average decentralized H ∞ PID team formation tracking performance in (5) is calculated as 0.118 and the average decentralized H ∞ fuzzy team formation tracking performance of 25 UAVs is calculated as 0.107. However, the design and computation complexity of PID control and T-S fuzzy control [15]    for 25 UAVs are very different. The PID control for 25UAVs needs 72 fuzzy rules and the T-S fuzzy control [15] for 25UAVs needs 1944 fuzzy rules, i.e., we need to solve 72 LMI s in (25) for PID control parameters of each UAV but 1944 LMIs for fuzzy control parameters of each UAV in [15]. Furthermore, we only need to compute u i (t) = K i TX i (t) in (6) at every instant time for PID control signal   of each UAV but u i (t) = 1944 k=1 h k (p i (t))K i,k X i (t) needs to be computed at every instant time for fuzzy control, where h k (p i (t)) is the complex ith Fuzzy interpolation function. Based on the above discussion, the decentralized H ∞ PID team formation tracking control for 25 UAVs can have less design and computation complexity than the decentralized H ∞ T-S fuzzy team formation tracking control for 25 UAVs but with almost the same team formation tracking performance.
Remark 6: The biggest barrier to apply the proposed algorithm to a real plant is the control input transformation. In the manuscript, we considered four control inputs, i.e., total thrust force of ith UAV,and three rotation torques of UAV. However, the UAV is controlled by the four rotors of UAV in the real situation. Therefore, the transformation between the theoretical control inputs and real control inputs should be tested by some practical experiments.

VI. CONCLUSION
In this study, the decentralized H ∞ PID team formation tracking strategy is proposed for large-scale stochastic quadrotor UAVs under external disturbance, intrinsic stochastic fluctuation and trailing vortex coupling. By using the virtual leader concept, the reference trajectory of each UAV is generated by the combination of the virtual leader trajectory with an adequate time-varying offset for each UAV to form a desired flight formation shape. Then, by the proposed optimal H ∞ decentralized PID team formation controller for each UAV, it can efficiently reject the effect of external disturbance and trailing vortex coupling from the neighboring quadrotor UAVs on the team formation reference tracking performance simultaneously. In order to avoid solving a complex nonlinear HJI for PID control of each UAV, the T-S fuzzy interpolation method is employed to approximate the stochastic quadrotor UAV so that the HJI can be transformed to a set of BMIs. Furthermore, by using a novel variable transform method, the optimal decentralized H ∞ PID team formation tracking control design problem of large-scale UAVs can be derived in terms of a set of independent LMIs-constrained optimization problems for each quadrotor UAV in the team formation, which can be efficiently solved in a single run via the current convex optimization technique of LMI toolbox in MATLAB without any tuning procedure of conventional PID design to significantly simplify the design procedure. A simulation example of 25 team formation quadrotor UAVs is given to validate the reference tracking performance of the proposed method in comparison with the conventional decentralized H ∞ fuzzy tracking control scheme. From the simulation results, the decentralized H ∞ PID team formation control for 25 UAVs in this study has a less design and computation complexity than the fuzzy control for 25 UAVs but with almost the same team formation tracking performance. Future works will focus on the decentralized H ∞ adaptive or attack-tolerant PID team formation tracking control of largescale quadrotor UAVs under event truggered control mechanism and Fos attacks [25], [26], [27] or Dos attacks [27].