Fast Cooperative Trajectory Generation of Unmanned Aerial Vehicles Using a Bezier Curve-Based Shaping Method

This study uses a Bezier curve-based shaping (BCBS) approach based on the dynamic model to quickly generate 3D cooperative trajectories for unmanned aerial vehicles (UAVs). It can facilitate the solution of the flight trajectory problem that occurs after the reallocation of mission points. Trajectory coordinates are expanded based on the Bezier curves that naturally satisfy boundary conditions (BCs); trajectories with the minimum flight time and those satisfying the dynamic constraints and collision-free conditions can be obtained by adjusting the Bezier coefficients. A three-UAV flight simulation is carried out to verify the effectiveness of the proposed method. Further, the performance of the proposed method is compared with that of the Gauss pseudospectral method (GPM). The simulation results of the BCBS method are used as the initial values for GPM optimization. It is demonstrated that the BCBS method requires a lower computation time than the direct solver, which is only 0.07% of the latter, and obtains similar optimization results (3.34% difference). This is considerably important for the rapid generation of optimized trajectories with the limited computing power of onboard computers. Furthermore, this method is expected to achieve online collaborative trajectory generation for multiple UAVs in view of its high computational efficiency.


I. INTRODUCTION
Unmanned aerial vehicles (UAVs) have been developed rapidly and used widely in recent decades [1]. Compared to manned combat units, UAVs provide the advantages of low cost and zero casualties. Therefore, UAVs are widely used in the military and civilian fields [2] such as for reconnaissance, surveillance, interference, relay communications, forest fire detection, and meteorological observation. With the increasing complexities of application scenarios, it is difficult for a single UAV to perform missions that can satisfy demands. Multiple UAVs have greater flexibility and adaptability and can perform better in some mission scenarios. Multi-UAV collaboration has become an important direction for UAV The associate editor coordinating the review of this manuscript and approving it for publication was Abderrahmane Lakas . technology, and in the future, it is expected that UAVs will be developed with a focus on clustering, autonomy, and high intelligence.
However, complex and variable multi-vehicle missions pose challenges for UAV technology. In an unexpected situation, the mission of UAVs may be disturbed or even changed. UAVs need to quickly generate flyable cooperative trajectories to the new mission points in emergencies. Therefore the completion of a complex multi-UAV cooperative mission requires a feasible UAV cooperative trajectory planning algorithm, which can generate safe and effective multi-UAV cooperative trajectories and provide UAVs a certain degree of autonomy. Further, because of the limited computing power of the onboard computer, the calculation amount of the UAV cooperative trajectory planning algorithm needs to be as small as possible. And the difficulty with multi-UAV cooperation is the need to improve the coordination of the UAV group such that it can avoid collisions between individuals under complex circumstances and ensure cooperative operation.
Trajectory planning algorithms generally use an optimal control theory. Popular numerical solutions of the optimal control theory include the direct (''discretise, then optimise'') and indirect methods (''optimise, then discretise'') [1], both of which need to be based on dynamics and involve a large amount of calculation (integral calculation). The pseudospectral optimal control theory is a popular strategy in the direct method, and it has been used to solve many aviation and aerospace problems [3]- [5]. In more recent developments, many improvements have been suggested for the pseudospectral method. Zhang et al. [6] developed an algorithm that integrates the Legendre pseudospectral method and the artificial potential field method. Their algorithm effectively solved the trajectory planning problem of the UAV in an obstacle-rich environment. Rogowski and Maroński [7] used the Chebyshev pseudospectral method to generate a trajectory for a glider in the vertical plane. Yang et al. [8] used Gauss pseudospectral method (GPM) in the trajectory optimization of a ramjet-powered vehicle. Sun et al. [9] optimized the 6-degree freedom trajectory of the parafoil delivery system with GPM. Other direct methods include the convex optimization method [10]. Pepy and Hérissé [11] proposed an algorithm based on an indirect shooting method.
The geometric method corresponds to the method using the optimal control theory. Babel [12] used shortest path algorithms for network optimization. They assumed that there are many line elements in space; connecting the directed line elements can act as the trajectory of the UAV. Other studies used Dubins paths, considering constraints such as obstacle avoidance, collision avoidance, and curvature limitation [13]- [16]. However, all UAVs move at a constant speed, and therefore, the optimization target is the path length. The Pythagorean hodograph (PH) method is another effective geometric method used for cooperative path planning [17]- [19]; dynamical constraints are described using geometric differential characteristics such as curvature and torsion. Although these purely geometric methods have high computational efficiency, they do not consider the dynamics of the vehicles. In the task planning for designing a large number of UAV trajectories, the most effective task allocation strategy can be quickly filtered using a purely geometric method; however, it may cause some trajectories to be less flyable.
The Bezier curve-based shaping (BCBS) approach presented in this paper employs a combination of geometry and dynamics, and it can, therefore, be called a geometric-dynamics method [20]. Inspired by the introduction of the concept of the shape function in Petropoulos and Longuski's article [21], the Bezier approach was proposed in our previous paper [22]. Several works of literature [23], [24] employ Bezier curves for the path planning of robots. Bezier curves have been used to connect straight lines or circles to make the curvature of the UAV trajectories smoother and more continuous [25]- [27]; other similar methods include B-spline [28]. Yu et al. [29] used the temporal-spatial Bezier curve to optimize the trajectory of multiple UAVs directly, and they obtained the suboptimal time for multiple vehicles to arrive simultaneously. Under the same conditions, the arrival time calculated using the temporal-spatial Bezier curve is shorter, and the curvature of the trajectory is more continuous compared with [30]. Although the temporal-spatial Bezier curve considers some dynamic constraints and collision avoidance, this method is a cubic Bezier curve. Under the condition that the boundary conditions (BCs) are determined, a cubic Bezier function can determine a unique curve without optimizing the Bezier coefficients.
The contributions of this study are as follows. First, this paper considers the dynamics model and constraints of the UAV with the collision avoidance of cooperative trajectories, instead of using the kinematics model in other papers [13]- [16] where the constraints on the turning radius and path length cannot consider the dynamic characteristics of the UAV. The method proposed in this paper takes thrust, angle of attack, and roll angle as control variables, which can be expressed by the curve coordinates and their derivatives obtained by the Bezier method. Second, the BCBS approach is adopted to design a multi-UAV cooperative trajectory quickly and efficiently. Assuming that the trajectory of the UAV is in the form of the Bezier function, the optimized trajectory can be obtained by optimizing the Bezier coefficients. The characteristics of the Bezier curve can satisfy the dynamic constraints and BCs of the UAV. A simulation scenario, wherein three UAVs depart and arrive simultaneously and may collide with each other, is designed to verify the feasibility of the proposed method. In the simulation, the Bezier method can quickly obtain the trajectory of the minimum flight time.
The remainder of this paper is organized as follows. Section 2 describes the dynamic model of the UAV, dynamic constraints, and collision-free conditions. Section 3 introduces a 3D cooperative trajectories generation framework using the Bezier approach, which solves the trajectory optimization problem using nonlinear programming problem (NLP) for Bezier coefficients. Section 4 describes a simulation that is performed in the set scenario to verify the feasibility of the method; the results are compared with those of GPM. Section 5 presents more simulated cases to test and verify the applicability and effectiveness of the method. Finally, some concluding remarks are described.

II. PROBLEM DESCRIPTION A. MISSION FRAMEWORK OF MULTI-VEHICLE
The multi-UAV mission framework described in this paper is a centralized architecture. For a general multi-UAV mission framework, the ground station generates all the trajectories and sends them to each vehicle before and during the mission (Fig. 1). However, the mission environment is complex and prone to changes. Once an emergency occurs, the ground station may not have time to change trajectories as soon as possible due to the long distance between the ground VOLUME 10, 2022 station and the UAV. UAVs need to have the ability to quickly generate trajectories autonomously. In the framework without a ground station (Fig. 2), the vehicles exchange and transmit information such as their current position and speed in real time. In the UAV group, any UAV can be used as the central vehicle to command and dispatch other UAVs. The central vehicle can be selected among the UAVs in the middle of the UAV group, which can facilitate the connection between the central unit and other vehicles. When there is an emergency that require the trajectories of the UAVs to be temporarily changed, the central vehicle quickly calculates a set of flyable flight trajectories based on the flight status of all the UAVs and transmits them to the extensions.

B. DYNAMIC MODEL OF UNMANNED AERIAL VEHICLE (UAV)
A three-dimensional coordinate system needs to be determined to describe the position and motion of the vehicles conveniently. During missions, the rotation and curvature of the earth can be ignored because the flight distance and altitude are within a small range.
First, a three-dimensional inertial coordinate system is established in space, where the OX axis and OY axis are on the horizontal plane, and the OZ axis is perpendicular to the horizontal plane and points in the opposite direction of gravity. Then a speed coordinate system is established; the UAV speed vector coincides with the OX v axis. The transformation relationship between the two coordinate systems can be expressed by the flight path angle γ and the heading angle ψ, which are the vertical and horizontal projection angles of OX and OX v , respectively.
It is assumed that the mass of the vehicles does not change during flight, there is no sideslip, and the angle of attack is relatively small. The 6-DoF dynamic equations arė where x, y, and z are the position coordinates of the vehicle in the inertial frame, and v is the speed. L and D are the lift and drag forces of the UAV, respectively.
S is the effective area, and C L is the lift coefficient. In the range where the angle of attack is small, C L can be regarded as proportional to the angle of attack α. Drag coefficient C D includes the zero-lift drag coefficient and the part related to lift. K is a constant. In addition, F T , α, and φ can be passed to the autopilot of the vehicle as control variables.(x, y, z, v, γ , ψ) are regarded as state variables.
The flight path and heading angles can be expressed as Further, the derivatives of γ and ψ arė whereẍ,ÿ, andz can be obtained using the variables of the Bezier method. The higher-order part of α can be ignored when the angle of attack is small. So it can be assumed that sin α ≈ α and cos α ≈ 1. Then Eq.(1) can be solved by the coordinates of the trajectories and their derivative.
wherev =ẋẍ +ẏÿ +żz ẋ 2 +ẏ 2 +ż 2 After obtaining the flight trajectories of the UAVs using the Bezier method, the control variables of the trajectories can be solved. The flyable trajectories must satisfy the dynamic constraints of the vehicles, such as speed, acceleration, and flight angles, moreover, the control input of the UAV also has a certain range of restrictions. These constraints can be described as where a = ẍ 2 +ÿ 2 +z 2 . For multiple UAVs, each vehicle has corresponding constraints based on its capability.

C. COOPERATION CONSTRAINTS
As the number of vehicles in the mission increases, the cooperation among the vehicles need to be maintained to complete the mission, which includes time and space cooperation [31]. Some complex tasks require vehicles to reach the target points simultaneously or sequentially [32], [33]. This paper focuses on the study of multiple vehicles reaching the target points at the same time. Therefore, there are the following constraints on the flight time T i of N UAVs.
Space cooperation requires to have a safe distance between the UAVs. If the vehicles cannot maintain a safe distance, it can affect flight control and cause collisions. The distance constraint can be satisfied by maintaining the distance between the vehicles to be greater than the minimum safe distance. The interval is calculated for the time where the mission times of two UAVs overlap. The vehicle position vector is p, and the safe distance constraint between the i-th and j-th vehicles can be written as where t i,i , and t i,j indicate the mission start time of the i-th and j-th vehicles, and correspondingly, t f ,i and t f ,j indicate the mission end times. Further, E d denotes the minimum safety distance, which is far greater than the required minimum safety distance. It can avoid the insufficient safety distance caused by errors resulting from control, discrete curve, air flow, etc.

III. BEZIER CURVE-BASED SHAPING APPROACH A. STATES APPROXIMATION
A vehicle describing the flight trajectory generated by the BCBS method is considered as an example. First, the dimensionless parameter τ = t/T ∈ [0, 1] of the flight time t is introduced, where T is the total flight time of the vehicle. Therefore, the 12 BCs of the trajectory generation problem can be written as where the normal subscripts ''i'' and ''f '' represent the initial and final states; the superscripts˙and represent the derivative of the dimensionless time parameter τ and the flight time t, respectively. The three-dimensional position coordinates of the vehicle flight trajectory can be expanded into the form of the Bezier curve function using the BCBS method [22].
where n x , n y , and n z are the degrees of the Bezier function; P x,j , P y,j , and P z,j are the control coefficients of the Bezier curve; and B x,j , B y,j , and B z,j are Bezier basis functions [34], which can be written as follows (Examples are given in the x dimension to avoid repetition): From Eq.(10), the first derivative of x with respect to the time of flight t can be obtained as where B x,j (τ ) can be obtained from Eq. (11).
Combining Eq.(9), Eq.(10), and Eq. (14), we get Hence, the four Bezier coefficients can be parameterized as The Bezier curve can be adjusted by changing the remaining n x − 3 Bezier coefficients in the optimization process. The roots of the mth-degree Legendre polynomial are used to discretize dimensionless time τ .
The coordinates [x, y, z] of the vehicle flight trajectory are represented in the form of matrix products. For example, x is written as are the Bezier coefficients. P x,0 , P x,1 , P x,n x −1 , and P x,n x are the known coefficients, as described in Eq.(16), and [X x ] (n x −3)×1 = P x,2 · · · P x,n x −2 T are the unknown coefficients. Matrices [B x ] m×(n x +1) , B x m×(n x +1) and B x m×(n x +1) can be obtained by substituting τ j (j ∈ [1, m]).
Variables n x , n y , n z , and m are required to be determined before beginning the trajectory generation calculation. The coefficient matrices ([B x ] m×(n x +1) , B x m×(n x +1) , and B x m×(n x +1) ) are constant and are calculated once before starting the loop optimization calculation.
If n x = n y = n z = 3, the shape of the flight trajectory is fixed by the BCs, and all the coefficients are determined. If n x > 3, n y > 3, andn z > 3, there are unknown coefficients ([X x ] (n x −3)×1 , X y (n y −3)×1 , [X z ] (n z −3)×1 ) which need be optimized to satisfy the dynamic constraints.
The above-mentioned descriptions of vehicle trajectory generation are based on a single vehicle. In this study, the time limitation problem of cooperative trajectory optimization is that N UAVs begin missions at the same time and reach their respective flight end points simultaneously. The optimization target is to generate the minimum flight time trajectories under the condition that all UAVs satisfy the dynamic constraints (Eq.(6)), and no collisions occur (Eq. (8)). Hence, the trajectory generation NLP problem of N vehicles can be described as where [X xi ] , X yi , and [X zi ] are the unknown Bezier coefficients of the ith vehicle. The constraints come from Eq. (6)(7)(8).
For the free-time problem of trajectory generation, the number of variables that need be optimized is N (n x + n y + n z − 9) + N , which includes the total flight time T i . The Bezier coefficients can represent the entire trajectory, which is optimized by the Bezier method. In the process of information transmission between the vehicles, only the Bezier coefficients (n x + n y + n z ) need to be transferred, instead of the points where the entire curve is discretized in a certain time scale, which greatly reduces the amount of data for information interaction between the vehicles. After the other vehicles receive the Bezier coefficients from the central vehicle, they can quickly calculate the point information at the next moment or the discrete point information of the entire flight curve according to the Bezier function.

B. INITIALIZATION OF UNKNOWN COEFFICIENTS
Before the optimization, initialize the variables ought to be optimized within a reasonable range. This paper discusses the free-time optimization problem. The approximated flight time T APP can be estimated from the median flight speed of a randomly selected vehicle, and the maximum acceleration limit is required to be considered.
where S is the straight line length between the starting and final points. The approximate coordinates (x APP , y APP , andz APP ) at the m Legendre-Gauss discretization points of each vehicle can be estimated by a third-order Bezier curve (n x = n y = n z = 3).
Considering the BCs (Eq. (9)), the coefficients are as follow.  (22). Then, the initial guess for the unknown coefficients in the Bezier curves can be obtained as

A. CASE STUDY
The proposed method is a framework for the optimization of the cooperative trajectory of the UAVs. A 3D simulation flight mission is set up where three fixed-wing UAVs reach the patrol mission target points from their respective start points. The three UAVs are required to start the mission at the same time and reach the end points simultaneously as well. The goal is to determine the trajectory with the minimum flight time that satisfies the dynamics constraints. The simulation scenarios and some parameters refer to the data in Choe's paper [20]. Choe et al. used the PH Bezier curve to describe the UAV trajectories, and they expressed the timing law using the quadratic Bezier polynomials. Through the above method, all variables of the UAV are described in the form of the PH Bezier polynomials. However, Choe et al. only consider the kinematics model of the UAVs, and the adjustable range of the curves is relatively small. Our method uses the geometrical characteristics of the Bezier curves combined with dynamic equations to generate trajectories under dynamic constraints and boundary conditions. The control variables  of the UAVs can be directly obtained through BCBS. The two methods have different applications to Bezier curves. The BCs and dynamic constraints are listed in Table 1 and 2. Some aerodynamic parameters are quoted from [10]. From the trajectory in Fig. 3, collisions will occur at the midpoints of the flight trajectories if the collision avoidance constraint is not considered in this simulation scenario. In the following simulation, the orders of the Bezier curve function are n x = n y = n z = 6, and the time parameter τ is discretized into 100 points. The simulation results of the BCBS method are used as the initial values of the GPM. In the trajectories generated by using GPM, the Legendre-Gauss points are set to 50. The NLP of both methods are solved using the sequence quadratic program method. All simulations in the paper are executed on a computer with the following specifications: 16.0 GB RAM, i7-10710U Intel processor at 1.10 GHz.
As shown in Fig. 4, in order to avoid the particularity of a single simulation, both methods have performed 100 Monte Carlo simulations. Fig. 5 shows the flight trajectories and space between vehicles calculated using the BCBS method.
Further, Fig. 6 shows the results calculated using the GPM. The figure shows that the trajectories of both the methods have a collision avoidance flight tendency, and the distances between the vehicles are greater than the minimum safety distance. The resulting curve of the Bezier method is more detailed because the Bezier method has more discrete points   than GPM. If GPM uses more discrete points, the calculation time will increase significantly. The flight time obtained by using the Bezier method is 197.8359 s, and that obtained using GPM is 191.444 s. The calculation time spent by the Bezier method and GPM is 0.588 s and 754.7936 s, respectively. The difference between the results calculated by the two methods is about 3.34%; however, the calculation time of the Bezier method is only 0.078% of that used to optimize and generate the trajectories using GPM further.
Figs. 7 and 8 show the dynamics variables of the vehicles with the Bezier method and the GPM, respectively. The figures indicate that the results of using the Bezier method and GPM can satisfy dynamics constraints. However, the results using the Bezier method are smoother, and there are no mutations, which is more favorable for the UAV's flight control.
To avoid accidental results in the simulation, some parameters of the simulation cases in Section 4 and performed    Table. 3. Case 1 is the original case in Section 4. Case 2-4 changed multiple coordinates and the total flight distance; Case 5 changed the maximum thrust (F T max ); Case 6-9 changed the initial speed (v i ); and Cases 10-13 changed the maximum speed (v max ). In the simulation of the above cases, the Bezier method required less calculation time (on an average, 0.08% of that GPM takes) and received relatively satisfactory results. The average difference of the calculated total flight time was only 3.09%. Further, the GPM uses the calculation result of the Bezier method as the initial value for further optimization.  The calculated results using both methods can satisfy the dynamic constraints. Moreover, the curves calculated by the Bezier method were continuous and smooth because more discrete points were used. Furthermore, a more continuous trajectory was found to be beneficial for avoiding collisions and helpful for flight control. If the GPM has more discrete points, the calculation time will be much longer. The small amount of calculation for the Bezier method was attributed to the coefficient matrix requiring one single calculation before optimization; the result was then substituted into the optimization iteration. Another reason was the geometry of the Bezier method, which does not require a large number of integration operations. In the optimization process, adjusting the Bezier coefficient can make the trajectory meet the requirements. When calculating the coordinated trajectories of more UAVs, the difference in the amount of calculation will increase more times.
In particular, the simulation case 5 changes the maximum thrust limit, which directly reflects the dynamic constraints. The trajectories are shown in Fig. 9, where it can be seen that the UAVs can only reach the target points by circling far because of insufficient thrust. The trajectory simulation of the kinematics model cannot realize the limitation of the thrust and angle of attack, which can only be indirectly restricted by the maximum flight path angle and the turning radius reflected by the dynamics model.

B. INFLUENCE OF BEZIER ORDER
Simulations on the influence of the order of the Bezier function on the calculation results and calculation time have been carried out. In the simulation, only the order is changed, and the other conditions are not changed. The results are shown in Table. 4. Higher-order Bezier functions can obtain better results, however, the calculation time is also significantly increased higher. The difference between the simulation results of the 6th-order and 20th-order functions is only 3.07%. Higher-order Bezier functions have higher degrees of freedom and can fine-tune the flight trajectory curve. But there are also more coefficients involved in the optimization process, which increases the amount of calculation. In order to improve calculation efficiency, a smaller order needs to be selected. However, if the order is too low, the error of the optimization result will be larger. The selection of the order of the Bezier function requires a balance between the accuracy of the calculation result and the calculation efficiency.

V. CONCLUSION
The rapid generation of 3D cooperative trajectories of multiple UAVs using the proposed BCBS method was presented in this paper. The BCBS method had both the authenticity of dynamics and the quickness of calculation. The trajectory coordinates were expanded into the form of the Bezier function, and the desired trajectories were optimized by adjusting the Bezier coefficient. A simulation was performed, and the results were compared with those of the GPM to verify the feasibility of the method. The Bezier method required less calculation time (0.07%) to generate the cooperative trajectories with the minimum flight time that satisfies the dynamics constraints and avoids collisions. Furthermore, the optimization results were found to be satisfactory, with a difference was only 3.34%. The trajectories obtained by the Bezier method were more continuous and smoother than those obtained using GPM, which is beneficial for flight control and collision avoidance. The relatively small amount of computations required by the Bezier method is of considerable importance for the rapid generation of optimized trajectories with onboard computers that have limited computing power. Therefore, it is expected that the proposed method will be employed for real-time online trajectory generation.