Distributed Formation Control with Obstacle and Collision Avoidance for Hypersonic Gliding Vehicles Subject to Multiple Constraints

. Multiple hypersonic gliding vehicles ’ (HGVs ’ ) formation control problems with obstacle and collision avoidance are investigated in this paper, which is addressed in the stage of entry gliding. The originality of this paper stems from the formation control algorithm where constraints of dynamic pressure, heating rate, total aerodynamic load, control inputs, collision avoidance, obstacle avoidance, and the terminal states are considered simultaneously. The algorithm implements a control framework designed to be of two terms: distributed virtual controller and actual control input solver. The distributed virtual controller is based on distributed model predictive control with synchronous update strategy, where the virtual control signals are derived by the optimization simultaneously at each time step for each HGV under directed communication topology. Subsequently, according to the virtual control signals obtained, a coupled nonlinear equation set is solved to get actual control signals: each HGV ’ s bank angle together with the angle of attack. The actual control input solver adopts a feasible solution process to calculate the actual control signals while dealing with constraints. Finally, extensive numerical simulations are implemented to unveil the proposed algorithm ’ s performance and superiority.


Introduction
Hypersonic glide vehicles (HGVs), which exhibit high speeds, strong maneuverability, and controllable trajectories, have acquired significant attention in recent years [1,2].The flight process of HGVs consists of four stages: the adjustment stage, the powered stage, the entry gliding stage, and the terminal guidance stage [3].The longest working range and duration make the entry gliding phase of HGVs a significant focus of attention in its flight process.In [4], to safely direct the vehicle to the intended terminal point under specified conditions, an entry guidance system is set.For HGV which has high ratios of lift-to-drag, a viable entry flight trajectory with threedegree-of-freedom is created via combining real-time updated commands from the adaptive guidance approach proposed in [5].The guidance system employs analytical schemes for a 3-D hypersonic gliding trajectory, allowing for rapid onboard planning of reference trajectories.As a result, it exhibits exceptional computational efficiency, as studied in [6][7][8].Some novel hybrid algorithms are evolved for real-time trajectory planning of a single HGV during the reentry phase, which offers improved reliability for online applications [9,10].Extensive research on the entry gliding phase has demonstrated the significant controllability of the entry trajectory.Furthermore, many entry trajectory planning algorithms have emerged to address diverse mission requirements under increasingly complex environments.Consequently, exploring new operational paradigms for HGVs during the entry gliding phase holds immense significance.
Compared to a single HGV, employing a cooperative strategy for HGVs can bring significant improvements in efficiency, information sharing, and robustness.Coordinating HGVs to arrive at the same target simultaneously with varied trajectories are a commonly encountered challenge during the cooperative flight of multi-HGVs.Significant progress has been made through the application of modern approaches including the impact time control approach [11], consensusbased guidance [12], and trajectory planning approach [13], which have yielded interesting results.In [14,15], the focus is on achieving precise control over both the timing and angle of impact for multiple HGVs.Through merging the decentralized and centralized communication topologies, an integrated framework for the cooperative guidance presented in [16] addressed the challenge of cooperative attack with multimissiles on a single stationary target.In [17], the design and analysis of a distributed cooperative guidance strategy of encirclement hunting for multivehicles in a leader-follower coordination structure are investigated, which effectively engages and tracks a maneuvering target using a cooperative approach.Due to the challenges associated with handling complex dynamics and underactuated problems brought by HGVs, there is limited research dedicated to the HGVs' formation control.In [18], with normal positions as the coordination variables, the second-order consensus protocol is implemented to the formation controller.Following multiple HGVs' formation flight framework, with the hierarchical control theory, the fixed-time stability is used for a globally fixedtime formation control scheme [19].Considering the lack of thrust to counteract the complicated impact for lateral maneuvers of longitudinal motion, a feasible flight mode called rendezvous and formation flight mode is raised [20], where the corresponding guidance strategy is set as a combination of online guidance and offline planning.According to the methods presented in [18][19][20], various flight modes are utilized to address the challenge of formation control for multiple HGVs.These methods aim to generate and maintain desired formation configurations by an online algorithm and only depend on the neighbor-to-neighbor communication.It means that it is being called upon to require an online algorithm with faster computing power and a smaller communication load for multiple HGVs' formation control.Furthermore, due to the complexity of the dynamics involved, the formation controllers designed by the above methods may have limited consideration for certain constraints.Therefore, studying formation control with multiple constraints for HGVs is indeed important and intriguing.
Considering the practical applications of HGVs, it is crucial to account for multiple constraints that must be satisfied during the entry gliding phase.Trajectory constraints for HGVs arise from various hardware limitations and factors related to the vehicle's dynamics, aerodynamics, and thermal characteristics, such as state and control constraints, overload constraints, dynamic pressure constraints, and aerodynamic thermal constraints.One way to handle the complicated multiple flight constraints is the introduction of the quasi equilibrium gliding condition [21].In [22], under complicated multiconstraint, the trajectory optimization problems were transformed to the problems of nonlinear optimal control, which were handled with Gauss pseudospectral method.In [23], researchers focused on studying the analytical solution for the HGVs' gliding problem, which aimed to develop a solution that allows for real-time boundary corridors along with control variables that can be online corrected.Furthermore, the stochastic gradient PSO method [24,25] together with the pigeon-inspired optimization method [26] has been proven to be useful in acquiring optimal trajectories that satisfy multiple constraints.In [24], for hypersonic glide vehicles under significant constraints, the stochastic gradient particle swarm optimization (SGPSO) method is employed as a global optimization method to quickly produce smooth and feasible trajectories.Firstly, a velocity-dependent profile for bank angle is devised to narrow down the search space for parameters using the constrained particle swarm optimization method.Then, the reentry constraints are carried out by assigning an infinite fitness function value when particles exceed the maximum allowable values [25].By incorporating control profiles, the entry trajectory optimization question is effectively addressed using the pigeon-inspired optimization (PIO) algorithm [26].This approach allows for the determination of control profiles that satisfy equilibrium glide conditions, terminal conditions, and a range of load factor constraints while also minimizing the peak heat rate.Meanwhile, mission constraints such as specifying a particular destination or satisfying a particular flight trajectory are requirements that depend on different missions [27].The collision and obstacle avoidance problem is a fundamental research area in the domain of multiagent systems and has recently acquired significant attention in the robotics and control system communities.In [28], with the employment of uncertainty and disturbance estimator (UDE), the paper investigated a distributed formation control strategy for unmanned marine surface vehicle (MSV) system.The algorithm achieves the formation control objective, which includes satisfaction of prescribed performance constraints, asymptotic convergence for formation errors, connectivity preservation constraints, and compliance with collision avoidance constraints.Through introducing the improved repulsive and attractive potential fields, the algorithm transforms the avoidance for no-fly zone and the problems for waypoint passage into a problem of determination for reference heading angle [29].Under the obstacle environment, for avoiding obstacle and collision and retaining the formation configuration, the creation of the null-space-based behavioral control structure is realized via defining the fundamental missions' priorities and computing the vectors for corresponding velocity, which addresses the problem for the coordinated control of spacecraft formation flying which has a leader [30].To summarize, obstacle and collision avoidance play an indispensable role in the formation flight of HGVs.Additionally, combined with the second paragraph, the interaction between trajectory constraints and mission constraints presents considerable challenges when aiming to achieve successful formation control for HGVs.
The issue of multiple constraints in formation control has been tackled by various strategies, among which model predictive control (MPC) has emerged as a practical approach due to the ability to handle constraints and deliver satisfactory control performance [31][32][33].However, the computational complexity associated with MPC is significant since it requires solving optimization problems numerically in a repetitive manner.Consequently, this high computational burden poses a challenge, particularly when applying MPC to multiagent systems, which demand substantial computing power.Distributed model predictive control (DMPC) is highlighted by its reduced computation cost, structural flexibility, and lower communication burden [34,35].In [36], it is proposed that only one agent sequentially solves its optimization problem 2 International Journal of Aerospace Engineering within each sampling period, which allows each agent to independently optimize its control actions while considering the constraints and objectives of the overall system.In [37], a novel cooperative distributed method for trajectory optimization is introduced for systems that include independent dynamics but hard constraints and coupled targets.In this approach, vehicles or agents solve their respective subproblems iteratively.Both [36,37] highlight the benefits of distributing the optimization among the agents, resulting in reduced computation cost and improved scalability in multiagent systems.However, in both the iterative and sequential solutions, for a cycle of all agents, it takes more time to obtain their optimal control inputs compared to a synchronous solution.In the DMPC method, every agent is able to synchronously solve its local optimization problem via communicating its assumed information with its neighbors, where the assumed information is updated simultaneously.Many studies have adopted the synchronous DMPC approach due to its advantages in terms of computational efficiency and coordination.In [38], with collision avoidance, aiming at the tracking and forming for homogeneous multiagent systems, a synthesis method of synchronous DMPC is introduced.Introducing hypothetical input and state trajectories yields a computationally tractable distributed optimization problem, which also solves the obstacle avoidance problem [39].[40] proposed solution through incorporating a control Lyapunov function (CLF) into DMPC, and this scheme inherits CLF's sturdy stability and optimizes the formation performance.Based on the results in the field of DMPC, the advantages of its ability to handle constraints and the synchronous update strategy have inspired the application of DMPC to address multiple HGVs' formation control problem.However, adopting DMPC in the context of multiple HGVs, where the system dynamics are nonlinear, underactuated, and highly coupled, poses its challenges and requires further investigation.
Driven by the aforementioned current state of research and impending future challenges discussed above, in this study, under directed communication topology, multiple HGVs' distributed formation control algorithm is proposed.Meanwhile, obstacle and collision avoidance with multiple constraints such as trajectory constraints, control constraints, and terminal state constraints which can be referred to as the achievement of the desired formation configuration are satisfied simultaneously in the entry gliding phase.Firstly, equations of motion of multiple HGVs and formation coordinate systems are created for the formation control problem.Subsequently, the combination of distributed virtual controller and actual control inputs solver forms a control framework designed for multiple HGVs' formation control.In the distributed virtual controller, with directed communication topology, the virtual control inputs are brought by setting the distributed model predictive control laws where multiple HGVs implement optimization simultaneously at each time step for every HGV.In the actual control input solver, the actual control inputs can be acquired through introducing a novel process to calculate synchronously while satisfying the multiple constraints for each HGV.Finally, the distributed formation control issues for the HGVs with obstacle and collision avoidance are resolved under the multiple constraints.In line with this, this paper makes four aspects of contribution when compared to prior studies in the field: (1) Under the directed communication topology, the distributed formation control algorithm is investigated, enabling the attainment of the desired formation configuration and ensuring collision and obstacle avoidance (2) Consider the premise of satisfaction for mission constraints consists of obstacle avoidance constraints, collision avoidance constraints, control constraints, and terminal state constraints.The set of trajectory constraints encompasses dynamic pressure constraints, heating rate constraints, and total aerodynamic load constraints, all of which are effectively satisfied simultaneously during the phase of formation control (3) For a cycle of all HGVs acquiring their actual input bank and attack angle, based on a framework composed of distributed virtual controller and actual input solver separately, a novel solution process to calculate the actual inputs is investigated with synchronous solutions.This proposed approach results in reduced computing burden and further handling of multiple constraints (4) For multiple HGVs' formation control problems, without regard to the constraints first involved in this paper, which include collision avoidance, trajectory constraints, and obstacle avoidance, the proposed formation control algorithm offers a faster convergence rate and improved formation accuracy compared to existing methods In this paper, the remaining sections' architecture is as follows: Section 2 introduces preliminaries pertinent to this study.Section 3 presents an in-depth exploration of the proposed algorithm.Section 4 demonstrates extensive simulation outcomes and conducts comparative analyses.Finally, Section 5 supplies a summary of this work.

Problem Statement and Preliminaries
2.1.Basic Theory.Define ℕ a ≜ 1, 2, ⋯, N a as the index set of HGVs with N a homogeneous hypersonic gliding vehicles.Let G = V , E, A be a weighted directed network with order N a and the set of nodes V = v 1 , v 2 , ⋯, v N a , the set of directed edges E ⊆ V × V , and a weighted adjacency matrix A = a ij ∈ ℝ N a ×N a .A directed edge E ij in G is indicated by the ordered pair of nodes v i , v j with terminal node v i and initial node v j , suggesting that node v i is able to obtain information from node v j .Define a ij = 1 if a directed edge v i , v j is found in G and a ii = 0. Degree diagonal matrix D = diag If between any pair of distinct nodes v i and v j in G, G = V , E, A is tightly connected, and there shows a directed path from v i to v j , i, j ∈ ℕ a .Definition 2. I n denotes the identity matrices of n order.diag ⋯ stands for the diagonal block matrix.Define 0 as the zero vector or zero matrix which has appropriate dimensions.For column vector x, x (i.e., x = x T x) represents the two-norm, and x P (i.e., x P = x T Px) signifies the P-weight two-norm of x.Assumption 3. Treat any HGV in the multiple HGVs as a node with number i.The communication topology G = V , E, G be a weighted directed network is established for multi-HGVs.The directed communication topology is tightly connected.Meanwhile, the communication delay is not considered in the multiple HGVs.

Equations of Motion.
During the entry gliding phase, the formation flight mission assumes that the form generation process occurs within a range of several hundred kilometers, encompassing the transition from the initial states to the desired formation configuration.Hence, for the multiple HGVs, the earth's curvature and rotation are not taken into consideration.
At the earth, designating a point on the surface as the origin, a plane rectangular coordinate system can be constructed.Define xOy as a horizontal plane, while Ox axis points from the origin to the target location.Oy axis is established by the right-hand relationship, where Oz axis is perpendicular to the local horizontal plane and then extends vertically upward, and then next, the relative motion and the reference frame between the HGV and the target location can be defined.Consider M i i ∈ ℕ a to signify the mass point of i-th HGV included in multi-HGVs.The equations of motion of HGVs and formation coordinate system are displayed in Figure 1.
The i-th HGV is termed as [18,19] x i = V i cos θ i cos ψ i , where V i represents the velocity, ψ i is defined as the heading angle, θ i stands for the flight path angle, g = 9 8 m/s 2 indicates the gravity acceleration, σ i represents the angle of bank, and D i and L i are the accelerations for drag and lift and calculated by where S i and m i , respectively, stand for the reference area and mass of the i-th HGV; α i is the angle of attack; C L and C D , respectively, denote the coefficients for lift and drag which depend on α i and V i ; and ρ i indicates the air density, which can be counted by where h 0 = 7110 m and ρ 0 = 1 225 kg/m 3 indicates the air density around sea level.

Multiple Constraints of Hypersonic Gliding
Vehicles.This section provides a detailed explanation of both mission constraints and trajectory constraints, where the mission constraints encompass terminal constraints, control constraints, obstacle avoidance constraints, and collision avoidance constraints.
as the desired formation configuration parameters for the multiple HGVs, where x d i , y d i , and z d i are the desired positions of the i-th HGV influenced by the desired formation configuration.The control problem for forming is converted to the terminal constraints 2.3.2.Trajectory Constraints.Trajectory constraints of multiple HGVs are described as follows: where k Q = 7 97 × 10 −8 indicates the heat transfer coefficient and denote the maximum limits of the total aerodynamic load n i , dynamic pressure q i , and stagnation point heating rate Q i , respectively.These are rigid constraints and must be satisfied.

Control Constraints.
In the rest of this paper, the control variables which are called the actual control inputs are the bank angle σ i and the angle of attack α i and for i-th HGV.Define α min i , α max i as the maximum and minimum constrained values for α i and σ max i and σ min i as the maximum 4 International Journal of Aerospace Engineering and minimum constrained values for σ i .The constraints for σ i and α i are as follows: 2.3.4.Collision Avoidance Constraints.Firstly, the expected relative positions between HGVs are in consensus with the collision avoidance constraints: Define N i as the set of neighbors of i-th HGV.Secondly, define R as the safety radius for each HGV; the following constraints should be satisfied to avoid collisions with another HGV: 2.4.Control Objective.For unpowered HGV, however, the velocity is so unmanageable that it is challenging for x i , y i , and z i to reach the expected value determined according expected formation configuration.Owing to the characteristics of initial states of the entry gliding stage, to decrease the conundrum of design for formation control, control of y i is ignored.Define Δx ij = x d j − x d i as the expected relative distance along Ox axis between HGV i and HGV j, where The objective is to design α i and σ i for every HGV i, such that the formation of multiple HGVs is stabilized to its desired position: and tracking errors of the HGVs reach consensus:

Formation Control Algorithm
Based on a control framework that is organized into two separate terms.Firstly, the distributed virtual controller is designed.Then, an actual control input solver is introduced to obtain the actual control inputs.Finally, a formation control algorithm is devised to consolidate all of the aforementioned aspects.
where for agent i, p i and q i ∈ ℝ 2 stand for the position and the velocity states.Control input is represented as u i ∈ ℝ 2 .The discrete-time system (15) can be rewritten into where k denotes the time instant and N ∈ ℝ is the length of predictive horizon and are the states and inputs.Define x i k ∈ ℝ denote the velocity and position along Ox axis.v z,i k and z i k are defined equally.Denote ξ i k + l k and u i k + l k as the predicted states and inputs at time step k + l from time step k, where ξ i k k = ξ i k .T is a constant sampling period; there have Multiple HGV systems are mapped to a discrete-time system (16) with the knowledge of x 0,i , y 0,i , z 0,i , V 0,i , θ 0,i , ψ 0,i T at current time t = T 0 s.Based on the mapped system ( 16), the DMPC is implemented to obtain the virtual signals u x i and u z i .It should be noted that, at every time step, all agents should be optimized simultaneously with the synchronous update scheme.To facilitate this process, assumed states are introduced as the information transmitted among neighbors instead of the real states.The assumed control inputs for HGV i are as follows: where κ i ⋅ is the terminal controller which is set in the form of where ξ d i is compatible with d ij and can be defined as i , 0 T .Substituting Eq. ( 19) into ( 16), define as error variable; then, the error system is given by Δξ Next, the assumed states of HGV i are defined as where To minimize the errors between current states and desired states for HGV i, the tracking index is expressed as where P i ∈ ℝ 4×4 is a symmetrical positive definite weight matrix for HGV i.
information for HGV i , and Q i = δ i I 4 indicates positive definite weight matrix.The consensus index is given by The index on control inputs for HGV i can be designed by where R i = τ i I 4 is the positive definite diagonal weight matrices.Finally, the index of the terminal cost is designed by where for the designated matrix H i which is negative definite, ϕ i = A′ and S i indicates a matrix which is positive definite and meets this Lyapunov equation: The objective function to be optimized at every time instant can be constructed as International Journal of Aerospace Engineering 3.1.2.Construction of Constraints.ξ N i is defined as the assumed state vector for neighbors of i-th HGV.The inputs of the mapped discrete-time system (16) are derived from the states of multiple HGVs system at the current time t = T 0 s.Then, the initial state constraints are defined as For HGV i, defining Limiting the deviation of the actual predictive positions and the corresponding assumed positions, the constraints for achieving collision avoidance are presented by The constraints to ensure the closed-loop system's stability are as follows: The uncontrollable velocity of a single HGV brings difficulty for obstacle avoidance.Hence, the constraints of collision avoidance can be simplified as However, in a real case, collision avoidance constraints cannot be satisfied by replacing in Eq. (30).Constraints (30) are redesigned as In a real case, the constraints of collision avoidance are as below: The handling method for collision avoidance is the same as Eq.(30).Definex o,m = x o,m , z o,m T ; the constraints for obstacle avoidance can be simplified as

33
where the simplified constraints still guarantee the obstacle avoidance constraints.Then, define x o,m = η y 0,i − y o,m x o,m , where where d o i is a positive constant.The introduction of x o m gives a faster convergence rate for the formation control by selecting a reasonable d o i .Finally, the obstacle avoidance constraints are as follows: To guarantee the closed-loop stability and the recursive feasibility for system (16).ℤ f i represents the corresponding terminal set where To meet the constraints of dynamic pressure and stagnation point heating rate, a relatively conservative strategy is to replace V i with V 0,i in Eqs. ( 6) and (7).The replaced constraints will be added in the DMPC at time k as follows.Define 6) and ( 7) are satisfied successfully because V i is a monotonous decrease function of time t.Define the state constraint set for HGV i as Accordingly, the constraints for stagnation point heating rate Q i and dynamic pressure q i can be given by 7 International Journal of Aerospace Engineering Finally, in the multiple HGV systems, it is imperative to apply specific constraints on the magnitude of virtual signals within the system (16).Define and u max i = u max i,1 , u max i,2 ; the input constraint set contains the origin which is expressed as where u i,1 and u i,2 represent the first and second component of u i , respectively.The input constraints are detailed as Remark 4. The problem of obstacle avoidance in this paper considers obstacle avoidance before the achievement of the desired formation configuration.Hence, the terminal constraints of HGVs are consistent with the obstacle avoidance constraints for the mission constraints, i.e., Assumption 5.The terminal set (36) is defined in a manner that ensures , and the terminal controller ( 19) is specifically devised to fulfill ξ i ∈ ℤ f i , κ i ξ i ∈ U i , and i ∈ ℕ a under all circumstances.Problem 6.At time step k (i.e., t = T 0 s), given ξ i k + l k , ξ d i , ξ N i k + l k , i ∈ ℕ a , and l = 0, 1, ⋯, N − 1, DMPC optimization problem is as below: At every time step k, each HGV possesses the ability to optimize with a distributed manner using Algorithm 1 with a synchronous update strategy.Assuming Assumption 5 holds, if there is at least a feasible solution available for every HGV at the initial time step k (i.e., t = T 0 s), all of the subsequent optimization problems will remain practicable.Moreover, every HGV converges to the expected formation configuration asymptotically while satisfying collision and obstacle avoidance constraints, stagnation point heating rate constraints, terminal constraints, and dynamic pressure constraints.
Proof.With the DMPC implemented by Algorithm 1.The optimization Problem 6's recursive feasibility and the closed-loop system's stability where system (16) adopted Algorithm 1 are analyzed, respectively.
(1) Feasibility: firstly, it is presumed that the optimization of Problem 6 for i-HGV is practicable at time step k.Accordingly, is the optimal control inputs, and ξ * i k + l k , l = 0, 1, ⋯, N, is the corresponding optimal states.Then, the assumed control inputs ûi k + 1 + l k + 1 and the assumed states ξ i k + 1 + l k + 1 , l = 0, 1, ⋯, N − 1, can be calculated according to Eqs. ( 18) and (20).Taking the information discussed above into account, the feasible control input Therefore, the input sequence u i k + 1 + l k + 1 , l = 0, 1, ⋯, N − 1 is a feasible solution to HGV i at time step k + 1, which means the feasibility of Problem 6 at time step k implies the feasibility at time step k + 1. Problem 6 is practicable at all future time instants by recursion (2) Stability: the optimal cost summation for all the HGVs is defined as a Lyapunov function of the whole system, which is indicated by Then, the summation of the costs is considered at time step k + 1 related to the feasible solutions u i k + 1 + l k + 1 = ûi k + 1 + l k + 1 for all the HGVs, represented by For the ease of subsequent proof, the objective function to be optimized at every time instant for each HGV is built 8 International Journal of Aerospace Engineering up of two components: stage objective function The difference between the optimal objective function and the feasible objective function at two successive time steps is Then, there have the following formulas: 49 Owing to the satisfaction of where ), , and H i , determine K i,1 , K i,2 , S i , and D i .x 0,i , y 0,i , z 0,i , V 0,i , θ 0,i , ψ 0,i T is delivered into mapped system (16).
Step 2. At the next time step k + 1, k ≥ 0, for HGV i ∈ ℕ a , the current states and ν i k + 1 .Then, Problem 6 is solved to get the optimal control inputs u * i k + 1 + l k + 1 , l = 0, 1, ⋯, N − 1, and the optimal states ξ * for the passing of virtual signals between the distributed virtual controller and the actual input solver.

Distributed virtual controller
No Applying z' i and V' i to Eq. ( 69) for the results of ~i max Calculate z' i and V' i by fourth order Runge-Kutta method = 0 = 0 Initialization: for HGV i ∈ ℕ a , Step 3 of Algorithm 1 is completed.Then, the virtual control signals u x i and u z i are obtained accordingly.Denote x 0,i , y 0,i , z 0,i , V 0,i , θ 0,i , ψ 0,i , α 0,i , σ 0,i T as the states of multiple HGVs at current time t = T 0 s.σ * i and α * i are the actual inputs to be solved at current time t = T 0 s.
Step 1. Substitute D i = ρ 0,i z 0,i V 2 0,i C D α 0,i , V 0,i S i / 2m i into Eq.( 62) to get L i cos σ i d .
Step 2. L i sin σ i d is obtained by applying D i and L i cos σ i d into Eq.(63).Then, L i d can be calculated by Eq. ( 64) with the knowledge of L i sin σ i d and L i cos σ i d .
Procedure 1: Procedure of the calculation for the actual inputs.
10 International Journal of Aerospace Engineering Also, there have Furthermore, define variables as follows: Then, it follows that where by the definition of ν i k , there have 52), (53), and (56) into (47).Given the solution's optimality, the following result is derived accordingly:

57
Define 0 < υ i < 1 as the convergence speed for closed-loop system.Then, the mapped multiple HGV systems converge to expected states influenced by the desired configuration asymptotically.Moreover, the terminal constraints, stagnation point heating rate constraints, dynamic pressure constraints, and collision and obstacle avoidance constraints will be obtained according to Theorem 7.

Actual Input Solver.
In the horizontal plane, define u x i as the virtual control input.Meanwhile, in the longitudinal plane, u z i represents the virtual control input.The virtual control signals u x i and u z i should be designed according to the neighbor information gathered through i-th hypersonic gliding vehicle from its neighboring vehicles.The dynamics of the coordinated variables of x i and z i , as well as their first and second derivatives versus time, can be counted by  After the optimization of the distributed virtual controller at current time t = T 0 s, applying the virtual control signals u x i and u z i into Eq.( 58), there have By solving Eqs. ( 59) and ( 60), the actual control inputs α * i and σ * i to be solved are further obtained.However, it is too complicated and time-consuming to solve Eqs. ( 59) and (60) adopting numerical methods.Hence, to decrease the complexity and coupling of calculation for α * i and σ * i , a feasible process is proposed.Consider σ i and α i are inherent with regard to L i sin σ i and L i cos σ i , which satisfy According to Eq. ( 61), the bank angle σ i can be ignored temporarily to decrease the complexity of calculation.Then, consider α i is implicit in terms of D i .Obviously, L i and D i are the function of α i .Although V i also affects L i and D i , V i of L i and D i are the same, which means that the actual input α * i can be calculated inversely according to Eq. ( 61) by calculating L i cos σ i d and L i sin σ i d satisfying Eqs.
(59) and (60).Define x 0,i , y 0,i , z 0,i , V 0,i , θ 0,i , ψ 0,i T as the states of multiple HGVs at current t = T 0 s.According to Eq. ( 60), Meanwhile, according to Eq. ( 59), However, there exists a coupling among Eqs.( 62) and (63).To get the angle of attack α * i , the designed procedures are given as follows.Initially, to gather the actual inputs at the current moment t = T 0 s, treat the drag accelerations D i as a parameter rather with the value of current time than the variable to be solved, i.e., D i = ρ 0,i z 0,i V 2 0,i C D α 0,i , V 0,i S i / 2m i .Then, L i cos σ i d are obtained by Eq. (62).Then, according to Eq. ( 63), L i sin σ i d is obtained easily.There have  13 International Journal of Aerospace Engineering Therefore, α * i and σ * i are calculated by using the current states of multi-HGVs.In line with Eq. ( 2), the angle of attack α * i can be counted inversely.
Additionally, to satisfy the control constraints for HGV i, limit the minimum and maximum constrained values of α * i based on the calculation result of Eq. ( 65) by Then, the bank angle σ * i can be counted with Additionally, to satisfy the control constraints for HGV i, limit the minimum and maximum constrained values of σ * i based on the calculation result of Eq. (67) by In summary, the implementation procedure of the calculation for the actual inputs σ * i and α * i is given as follows.It is noted that the constraint of aerodynamic load cannot be considered in the design of virtual controller because of coupling.To constrain the violation of the aerodynamic load, the following process are proposed.
At the time step k in the distributed virtual controller, the virtual signals u x i and u z i are obtained, the temporary variables α * i and σ * i are calculated by Procedure 1.Then, apply the temporary variables α * i and σ * i into the system dynamic (1) at the current time t = T 0 s.The state variables z i ′ and V i ′ at time t = T 0 + T s are obtained directly with the fourth-order Runge-Kutta method with initial conditions x 0,i , y 0,i , z 0,i , V 0,i , θ 0,i , ψ 0,i , α * i , σ * i T .Define α i as feasible actual input.Then, Accordingly, the maximum value of α i , which is α max i , is calculated successfully by Eq. (69).The actual signal of α * i is described as and the actual signal of σ * i is described as Table 4: Initial states of multiple HGVs.

Simulation Examples
Three simulation examples are provided in this section: (1) simulation of formation accomplished with different control parameters, (2) simulation of formation control satisfying multiple constraints, and (3) simulation comparison with existing methods.Example 1 unveils the characteristic of the proposed algorithm by setting different control parameters and conservative constraint parameters.The desired formation achieved in example 1 serves as reasonable initial states for example 2 that will be studied in the middle of entry gliding phase.Then, the successful achievements of obstacle and collision avoidance for formation control while dealing with multiple constraints simultaneously are demonstrated in example 2. Examples 3 and 1 share the same initial states for multiple HGVs.After successfully demonstrating the characteristics of the proposed formation control algorithm in example 1 by setting various parameters, example 3 utilizes the same parameters from example 1 to compare with other algorithms, which exhibits superior performance compared to other existing formation control strategies for the multiple HGVs.Consider N a = 3.The simulations in this study utilize the CAV-H model [1].Multiple HGVs' communication topology is shown in  1. t ∈ 0, 60 s is the initial phase where each HGV flies in maximum lift-drag ratio with α i = 10 ∘ and σ i = 0 ∘ .Then, the proposed algorithm is to be studied when t ≥ 60 s.The desired formation configuration is depicted in case 1 of Table 2. To unveil the characteristic of the proposed algorithm, various parameters which are vital are shown in Table 3.It is noted that to better illustrate the characteristics of the control performance, the investigations of collision avoidance, obstacle avoidance, and trajectory constraints are all ignored in this simulation by setting conservative constraint parameters.The rationale for conducting the simulation prior to the simulation in Section 4.2 is to establish reasonable initial states for the simulation that will be studied in the middle of entry gliding phase.The other parameters are selected as T = 1, δ i = 0 1, τ i = 0 1, υ i = 0 9, u min i = −30,−30 , and u max i = 30, 30 .Define case 1 as the desired formation configuration with control parameter 1, and case 2 and case 3 are defined similarly.Figure 3 shows the different formation flight trajectories from the proposed algorithm with different control parameters.The desired formation configuration can be successfully acquired with the proposed algorithm.Case 2 brings a faster convergence rate together with a higher control accuracy than case 1 and case 3.
Figure 5 shows that the constrains of α i and σ i are all satisfied in three cases.In Figure 6, time histories of V i are given.From Figures 4-6, it can be seen that the increase of N and P i not only increases the convergence rate but also decreases the energy consumed.However, a bigger N will bring higher computational load, and a bigger P i may bring chattering problem for the actual control inputs.Thus, the selection of optimal parameters for N and P i should be made by considering the actual circumstances.

Simulation of Formation Control Satisfying Multiple
Constraints.Then, the cardinal contributions of our work are the achievement for obstacle and collision avoidance for formation control while dealing with multiple constraints simultaneously.Therefore, this simulation uses three comparative experiments to reveal the function of the proposed algorithm.The initial states in this section are given in Table 4 according to Figure 6 in Section 4.1.Then, the desired formation configuration in this section is depicted in case 2 of Table 2.And the positions of obstacles are set in Table 5.
The constraints of avoidance for both obstacle and collision are all included in three cases.In the presence of two obstacles, case 1 tends to achieve the desired formation configurations without considering the trajectory constraints for dynamic pressure, stagnation point heating rate, and total aerodynamic load.With the same goal as case 1, the total aerodynamic load is considered additionally in case 2. Finally, in case 3, which corresponds to case 2, both dynamic pressure constraints and stagnation point heating rate constraints are encompassed additionally.The constraints in Eqs. ( 6)-( 8) are given as Q max = 900 kW/m 2 , q max = 40 kPa, and n max = 3.The other parameters are selected as T = 1, δ i = 0 1, τ i = 0 1, υ i = 0 9, u min  From Figure 10, collision avoidance is ensured; the relative distances d e ij = x i − x j between every HGV are greater than the safety distance 2R.In Figure 11, when t ≤ 15 s and y i ≤ 36000 m for each HGV, the relative distances between each HGV and obstacle 1 d e i,o,1 = x i − x o,1 have a safety distance.In Figure 12, similarly, when t ≤ 35 s and y i ≤ 81000 m for each HGV, the relative distances between each HGV and obstacle 2 d e i,o,2 = x i − x o,2 have a safety distance as well.Meanwhile, Figure 9 shows the histories of α i and σ i , where the constraints of the maximum value are satisfied.The mission constraints are satisfied successfully by the proposed algorithm in case 1, case 2, and case 3.
Moreover, the trajectory constraints' satisfaction is analyzed with the mission constraints' satisfaction.For case 1, Figure 13 shows that the aerodynamic load constraints are violated.In case 2, the aerodynamic load has a maximum value of 3, which means the aerodynamic load constraints are satisfied by the proposed algorithm.Although the trajectories of case 1 and case 2 are almost identical as exhibited in Figure 7, the aerodynamic load constraints should be considered seriously.Case 3 and case 2 both satisfy the aerodynamic load constraints n max i = 3 as shown in Figure 13.In Figure 14, both in case 1 and case 2, the dynamic pressure constraints are violated.Comparing case 3 with case 2, Figure 7 suggests that the significant change of the trajectory for multiple HGVs is the trajectory of HGV 2. To avoid the obstacle 2, HGV 2 flies beneath the obstacle in case 2 but flies above the obstacle in case 3 to achieve the obstacle avoidance.From Figure 14, the proposed algorithm enables HGVs flying with different trajectories with obstacle avoidance, which makes the violation of dynamic pressure constraints disappear in case 3. Hence, the multiple constraints  18 International Journal of Aerospace Engineering are satisfied simultaneously in case 3.One thing worth noting is that the stagnation point heating rate is all satisfied in three cases.Although the simulation does not give a contrast example to show the suppression of stagnation point heating rate constraints violation as displayed in Figures 13 and 14, the state constraint set designed in Eq. ( 37) and the satisfaction of Assumption 5 make the suppression of dynamic press constraints equally persuasive towards the stagnation point heating rate constraints as exhibited in Figure 15.In summary, the aforementioned results convincingly unveil the proposed formation control algorithm's strong performance in handling multiple constraints.

Simulation Comparison with Existing Approaches.
Finally, to verify the proposed control algorithm's superiority despite its ability to deal with multiple constraints.The comparison simulation is performed in this section.An insightful study for multiple HGVs is shown in [19], where a globally fixed-time formation control algorithm with the hierarchical control theory is introduced named fixed-time formation control, which ensures the fixed-time stability.Meanwhile, the proposed algorithm is named DMPC formation control to make comparison.The initial states as demonstrated in Table 1 and the initial phase are the same as in Section 4.1.The desired formation configuration is given in case 1 of Table 2.The parameters for the two methods are shown in Table 6.In this section, to make better comparisons by control variates, the proposed control algorithm does not consider the constraints discussed in Section 4.2.Thus, the obstacle avoidance constraints, the total aerodynamic load constraints, the collision avoidance constraints, the stagnation point heating rate constraints, and the dynamic pressure constraints are all ignored by setting conservative constraint values for the proposed algorithm.The other parameters are selected as T = 1, δ i = 0 1, τ i = 0 1, υ i = 0 9, u min i = −30,−60 , and u max i = 30, 60 .Figure 16 shows the comparison of the formation flight trajectories, which demonstrates the achievement of the desired formation configuration using the fixed-time formation control and the DMPC formation control.In Figure 17, time histories of the formation error variables E x i are displayed; it can be seen that E x i converge to the vicinity near zero within 140 s under the DMPC formation control.And E x i converge to the vicinity near zero after 200 s with the fixed-time formation control.
It can be seen in Figure 18 that E z i converge to the vicinity near zero within 130 s under the DMPC formation control and converge to the vicinity near zero within 150 s under the condition of the fixed-time formation control.From Figures 16-18, DMPC formation control gives a faster convergence rate and smaller steady-state value for formation error  Therefore, the DMPC formation control has a relatively fast convergence speed rather than the fixed-time formation control and a rapidly faster convergence speed and improved convergence precision rather than other methods.In Figure 19, the maximum maneuvering ability is the same for the two methods, where the limitations of the α i and σ i are the same and are satisfied during the formation control.Accordingly, the simulation illustrates that in spite of the ability to cope with multiple constraints during the formation flight, the proposed formation algorithm has a superiority in convergence speed and convergence precision.The proposed formation algorithm exhibits superior performance compared to other existing formation control strategies for the multiple HGVs.

Conclusion
In our work, a distributed formation control algorithm is investigated for multi-HGVs under multiple constraints consisting of mission constraints and trajectory constraints.The first difference between previous algorithms and the proposed algorithm is the fulfilment of the obstacle avoidance constraints, the collision avoidance constraints, and the terminal state constraints for the mission constraints for multiple HGV formation flights.The second difference is that, with the satisfaction of the necessary mission constraints, the trajectory constraints were considered additionally the first time for the formation flight of multi-HGVs.A distributed virtual controller and an actual control input solver were designed for the control framework, where the virtual controller introduced an optimization problem with a synchronous update strategy for the formation control in a distributed way and the actual control input solver adopted a novel solution process to calculate the actual control inputs while dealing with multiple constraints.Extensive simulation results disclose that the proposed algorithm has the ability to guarantee multiple constraints for formation flight with obstacle and collision avoidance simultaneously.And the superiority in the control performance is demonstrated.

. 3 . 5 .
Obstacle Avoidance Constraints.Obstacle avoidance can be defined as no-fly zones of geographic constraints that should be satisfied before the desired formation configuration is achieved.Define ℕ o ≜ 1, 2, ⋯, N o as the index set of all the obstacles and x o,m = x o,m , y o,m , z o,m T as the position of obstacle m.The obstacle constraints are as follows:

Figure 1 :
Figure 1: Model diagram of HGV in the formation coordinate system.
the satisfaction of constraints (43) except ξ i k + N k + 1 ∈ ℤ f i follows directly from the assumed control inputs' definition and Assumption 5.The positive invariance of the terminal set ℤ f i further results in the satisfaction of constraint

Step 4 .Algorithm 1 :
Return to Step 1. Distributed model predictive control algorithm.International Journal of Aerospace Engineering Generate the actual inputs ~⁎i and ~⁎i by Eqs.(70) and (71) Input the desired formation configuration parameters d ij and x -d i

Figure 2 :
Figure 2: Flowchart of the formation control algorithm.

Figure 6 :
Figure 6: Time histories of V i .

Figure 8 :
Figure 8: Historical simulation curves: E x i and E z i .

Figure 2 .
Then, each HGV's maximum maneuvering ability is limited by the limitation of the actual inputs, where α min i = 0 ∘ , α max i = 25 ∘ , σ min i = −90 ∘ , and σ max i = 90 ∘ .Define the formation error variables E x i = ∑ N a j=1,j≠i a ij x j − x i − Δx ij for x i and E z i = ∑ N a j=1,j≠i a ij z j − z i − Δz ij for z i .4.1.Simulation of Formation Accomplished with Different Control Parameters.The initial states for the multi-HGVs are exhibited in Table

Figure 18 :
Figure 18: Time histories of E z i .

3 Figure 19 :
Figure 19: Time histories of α i and σ i in example 3.

Table 1 :
Initial states of multiple HGVs.

Table 5 :
States of obstacles.
16 Journal of Aerospace Engineering demonstrates the histories of value and convergence rate for formation error variables, where E x i and E z i converge to a small constant with different convergence rate.Comparing case 2 with case 1, a bigger prediction horizon N brings a faster convergence rate.Comparing case 2 with case 3, a bigger weight matrix P i also brings a faster convergence rate.