Collision-free Multiple Unmanned Combat Aerial Vehicles Cooperative Trajectory Planning for Time-critical Missions using Differential Flatness Approach

This paper investigates the cooperative trajectory planning for multiple unmanned combat aerial vehicles in performing autonomous cooperative air-to-ground target attack missions. Firstly, the collision-free cooperative trajectory planning problem for time-critical missions is formulated as a cooperative trajectory optimal control problem (CTP-OCP), which is based on an approximate allowable attack region model, several constraints model, and a multi-criteria objective function. Next, a planning algorithm based on the differential flatness, B-spline curves and nonlinear programming is designed to solve the CTP-OCP. In particular, the notion of the virtual time is introduced to deal with the temporal constraints. Finally, the proposed approach is validated by two typical scenarios and the simulation results show the feasibility and effectiveness of the proposed planning approach.

Keywords:    Cooperative trajectory planning,   unmanned combat aerial vehicle,   differential flatness,  B-spline curve,  allowable attack region

Nowadays, it is an active research area to perform autonomous cooperative air-to-ground target attack (CA/GTA) missions by multiple unmanned combat aerial vehicles (multi-UCAV)1 . However, it is quite difficult to coordinate the strike operation of multi-UCAV, especially for time-critical missions that require precise timing and sequencing of tasks and operations. In order to accomplish the mission successfully, UCAVs need to generate detailed cooperative execution plans to lead themselves well to penetrate through enemy threat areas, avoid the collisions with obstacles or each other, fly into the allowable attack regions (AARs) simultaneously or in sequence, and then perform weapon delivery operations. The cooperative trajectory planning for CA/GTA is vital to achieve the mission goals. It is really one of the key challenging technologies, due to its high dimensionality, severe equality and inequality constraints involved, and the requirement of spatial-temporal cooperation of multi-UCAV, and has recently received extensive attentions2 .

To date, various algorithms3 - 7 have been developed to solve this cooperative planning problem, including several collision-avoidance techniques and time adjustment strategies. McLain3 , et al. used coordination variables and coordination functions based strategies to achieve cooperative timing among teams of vehicles, by coordinating the velocity and path length of each vehicle. Kaminer4 , et al. proposed a solution to the problem of coordinated control of multiple unmanned aerial vehicles (multi-UAV) to ensure collision-free maneuvers under strict spatial and temporal constraints. Bollino5 , et al. addressed the optimal path planning of UAVs in obstacle-rich environments and proposed the collision-free path planning for multi-UAV using optimal control techniques and pseudospectral methods. Lian6 introduced a differential flatness based approach to optimally formulate the cooperative path planning for multi-agent dynamical systems considering spatial and temporal constraints, and parameterized the curves by B-spline representations. Kuwata7 , et al. presented a cooperative distributed robust trajectory optimization approach, using RH-MILP with independent dynamics but coupled objectives and hard constraints. The above investigations have given several valuable strategies in the cooperative planning, but they failed to tackle the point-to-region cooperative trajectory planning for CA/GTA missions under consideration directly, which needs to integrate both the spatial and temporal constraints on the level of the trajectory planning. To address the trajectory planning, a novel cooperative trajectory planning algorithm for multi-UCAV in performing the CA/GTA missions is presented

Given a set of stationary ground targets in a terrain region, the mission objective is to plan several cooperative optimal or suboptimal, dynamically feasible flight trajectories from the initial points (IPs) into the AARs, such that multi-UCAV can effectively attack the targets with minimum time, maximum survivability and target destruction effectiveness. Therefore, many factors need to be considered synthetically in the cooperative trajectory planning. In this section, the cooperative trajectory planning problem is formulated after the modelling of AAR, constraints and the objective function.

2.1  Allowable Attack Region Model

For the point-to-region trajectory planning problem, the AARs of targets are defined as the areas where UCAVs can effectively perform weapon delivery operations. So, in order to plan the accurate and optimal attack trajectories for weapon delivery, the AARs and the delivery parameters need to be obtained. The AAR of the ith target TARi, denoted as R(TARi), is such a set of all feasible release states that TARi can be effectively attacked whenever the aircraft is in that states. R(TARi) can be formulated as an abstract 6-dimensional space8

$R\left(TA{R}_{i}\right)=\left\{x,y,h,v,\gamma ,\psi \right\}\subset {ℝ}^{6}\text{.}$      (1)

Obviously, R(TARi) is high-dimensional nonlinear space, which will make the problem solving much difficult. By presetting an appropriate weapon release speed, vr and the flight-path angle, ${\gamma }_{\text{r}}$ based on estimating the weapon impact effects and destruction requirements to the target, and predetermining release heading ${{\psi }^{\prime }}_{\text{r}}$ ,it can be reduced to a 3-dimensional space

(2)

2.2  Constraint Model

2.2.1 Maneuverability Constraints of UCAVs

The maneuverability constraints impact every phase during target attack missions, so they should be met for an executable plan. The constraints can be denoted as

(3)

where v is the aircraft velocity;$\gamma$ denotes the flight-path angle; denotes the heading angle$\psi$ ; denotes the roll angle $\mu$ ; nx and nh are the longitudinal and normal components of the load factor. and ${\left(\cdot \right)}_{\mathrm{min}}$ and ${\left(\cdot \right)}_{\mathrm{max}}$ represent the minimum and maximum boundary values of the value $\left(\cdot \right)$ .

2.2.2 Battlefield Environment Constraints

(a) Threat avoidance: without loss of generality, the threat impact range can be assumed as hemispheres approximately, thus threat avoidance constraint can be defined as

(4)

where the norm ${‖\cdot ‖}_{2}$ denotes the Euclidean distance between two points and$\left({x}_{i}^{\text{o}},{y}_{i}^{\text{o}},{h}_{i}^{\text{o}}\right)$ and ${R}_{i}^{\text{o}}$ denote the origin coordinates and the detection radius of the ith threat, respectively.
(b) No-fly zones (NFZs): herein, an infinite-length cylinder is used to describe the NFZs:

(5)

where $\left({x}_{i}^{\text{w}},{y}_{i}^{\text{w}}\right)$ and ${R}_{i}^{\text{w}}$ denote the origin coordinates and the radius of the ith NFZ.

2.2.3 Terminal Constraints

In the point-to-region trajectory planning, weapon delivery point (WDPt) as the trajectory terminal of CA/GTA needs to be in the AAR, i.e. satisfying the terminal constraints. The formula can be denoted as ${\text{WDPt}}_{i}\in R\text{(}TA{R}_{i}\text{)}$, that is

(6)

where tf is the terminal time, (xa, ya, ha) and (x(tf), y(tf), h(tf)) are the coordinates of AAR’s center point and WDPt, are the thresholds of errors.

2.2.4 Cooperative Constraints

During the mission, multi-UCAV should maintain a safe distance to avoid collision with each other, i.e. satisfying spatial constraint. The model can be denoted as

(7)

where ${\rho }^{j}\left({t}_{i}\right)=\left\{{x}^{j}\left({t}_{i}\right),{y}^{j}\left({t}_{i}\right),{h}^{j}\left({t}_{i}\right)\right\}$ is the spatial position of UCAVj at the time ti, is the minimum safety radius of UCAVj, and Nv is the total number of UCAVs.

The temporal constraints considered here include simultaneous arrival constraint and tight-sequencing constraint, which can be described as

(8)

where ${t}_{\text{f}}^{j}$ is the terminal time of UCAVj, ${\Delta }_{jk}$ is the arrival interval between UCAVs.

2.3   Objective Function

The objective function of each UCAVj can be defined by the weighted sum of the three separate running cost terms with appropriate weighting factors ${w}_{\text{t}}^{j}$ ${w}_{\text{r}}^{j}$ and ${w}_{\text{d}}^{j}$ , where ${w}_{\text{t}}^{j}+{w}_{\text{r}}^{j}+{w}_{\text{d}}^{j}=1$ . And the objective function of the entire team can defined as:

$J=\mathrm{min}\sum _{j=1}^{{N}_{v}}{J}^{j}=\mathrm{min}\sum _{j=1}^{{N}_{v}}\left({w}_{\text{t}}^{j}{J}_{\text{t}}^{j}+{w}_{\text{r}}^{j}{J}_{\text{prd}}^{j}+{w}_{\text{d}}^{j}{J}_{\text{dest}}^{j}\right),$          (9)

and three separate running cost terms can be respectively defined as

${J}_{\text{t}}^{j}=\frac{1}{{\Gamma }^{j}}{\int }_{{t}_{0}^{j}}^{{t}_{\text{f}}^{j}}dt=\frac{1}{{\Gamma }^{j}}\left({t}_{\text{f}}^{j}-{t}_{0}^{j}\right),$         (10)

(11)

${J}_{\text{dest}}^{j}={f}_{dest}\left({x}^{j}\left({t}_{\text{f}}^{j}\right)\right)=1-\sum _{i}^{{n}_{s}}{D}_{i}^{j}\left({x}^{j}\left({t}_{\text{f}}^{j}\right)\right)/{n}_{s},$          (12)

where ${t}_{0}^{j}$ and ${t}_{\text{f}}^{j}$ denote the initial and terminal time of UCAVj. The first term, ${J}_{\text{t}}^{j}$ denotes the total fight time of UCAVj $\Gamma =L/{v}_{\mathrm{min}}$ . is maximum flight time along the straight line (L denotes its length) from the IP to the center of AARs. The value is predetermined, and can be used to insure the first cost term consistent with others. The second term ${J}_{\text{prd}}^{j}$ , describes the average detection-probability of the nr-radar system to UCAVj, where ${P}_{\text{d}}^{j}\left(t,r\right)$ is the radar detection probability model between the trajectory point of UCAVj at the time t and the rth radar9. And the third term ${J}_{dest}^{j}$ , is the minimal form of the target damage probability, where ${D}_{i}^{j}\left({x}^{j}\left({t}_{\text{f}}^{j}\right)\right)$ denotes the target damage probability of UCAVj from the ith Monte Carlo simulation, and ns is the total simulation times.

2.4   Cooperative Trajectory Planning Problem Formulation

After establishing the above three models, the cooperative trajectory planning problem is formulated in the subsection. The problem under consideration is a cooperative scenario, consisting of Nv similar UCAVs, and the dynamics of UCAVj is given by

(13)

where $x={\left[x,y,h,v,\gamma ,\psi \right]}^{\text{T}}\in {ℝ}^{6}$ and $u={\left[\mu ,{n}_{\text{x}},{n}_{\text{h}}\right]}^{\text{T}}\in {ℝ}^{3}$ denote the state and control vectors, which are in accordance with the following UCAV kinematic and dynamics equations10

(14)

(15)

where x, y, h are the east, north, and up components of the earth-fixed reference frame, and denote longitude, latitude and altitude respectively; g denotes the gravity acceleration.Fig 1. shows the spatial relation of the states.

As mentioned above, to generate optimal or suboptimal cooperative trajectories, the trajectory planning problem for the CA/GTA missions can be formulated as a cooperative trajectory optimal control problem (CTP-OCP).

Figure 1. Relation of aircraft position, velocity, flight-path and heading angels.

2.4.1 Problem 1 (CTP-OCP)

Find the trajectories, which drive the system from given initial conditions to desired final conditions over time horizons [t0, tf], while the cooperative objective function is minimized as follows

$\begin{array}{l}\mathrm{min}J=\sum _{j=1}^{{N}_{\text{v}}}{J}^{j}\left(x,u\right)=\sum _{j=1}^{{N}_{\text{v}}}\left[{w}_{\text{t}}^{j}\left({t}_{\text{f}}^{j}-{t}_{0}^{j}\right)/\Gamma +{w}_{\text{r}}^{j}{\int }_{{t}_{0}^{j}}^{{t}_{\text{f}}^{j}}{\text{PRD}}^{j}\left(t\right)dt/\left({t}_{\text{f}}^{j}-{t}_{0}^{j}\right)+\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{w}_{\text{d}}^{j}{f}_{dest}\left({x}^{j}\left({t}_{\text{f}}^{j}\right)\right)\right]\end{array}$     (16)

subject to the dynamics equation (i.e., Eqn (13)) and the boundary constraints (i.e., the initial and terminal states (i.e., Eqn (6))

$\begin{array}{l}{\Phi }_{0}\left[{x}^{j}\left({t}_{0}^{j}\right),{u}^{j}\left({t}_{0}^{j}\right),{t}_{0}^{j}\right]=0,\\ {\Phi }_{1}\left[{x}^{j}\left({t}_{\text{f}}^{j}\right),{u}^{j}\left({t}_{\text{f}}^{j}\right),{t}_{\text{f}}^{j}\right]\le \text{0,}\end{array}$     (17)

and several inequality and equality constraints, the individual and cooperative constraints, including the state and control (i.e., Eqns (3)-(5) and Eqns (7)-(8)), are denoted as:

(18)

Note that, the terminal time tf is free. Thus the problem is a free terminal time problem. In addition, it is important in the problem formulation to scale the variables. The choice of scaling will balance the equations for numerical analysis, thus improving the accuracy of the solution and the computation time11.

To solve the previous CTP-OCP as quickly as possible, an efficient numerically cooperative trajectory planning algorithm is introduced, which combines several classical techniques, including differential flatness theory, B-spline curves, and nonlinear programming. In addition, for the time-critical cooperative missions, a novel strategy is proposed to handle the temporal constraints.

3.1  Time Cooperative Strategy

The time factor of trajectories is an argument of the state and control. To deal with temporal constraints, the time along the trajectories should be considered separately. In the work, the independent intermediate variable (called virtual time here $\tau \in \left[0,1\right]$ ) is introduced and described as

$\tau \triangleq \left(t-{t}_{0}\right)/\left({t}_{\text{f}}-{t}_{0}\right).$      (19)

Such that the trajectories can be generated in the virtual time domain from ${\tau }_{0}=0$ to ${\tau }_{\text{f}}=1$ .

For UCAVs are required to takeoff at the same time, it is assumed that the initial time of all UCAVs is zero (i.e.${t}_{\text{0}}=0$ ). Thus, the terminal time ${t}_{\text{f}}$ can be written as ${t}_{\text{f}}=t/\tau$ . That is ${t}_{\text{f}}$ , denotes the ratio between the true time variable and the newly defined virtual time variable. To coordinate the arrival time of all UCAVs, the terminal time can be defined as an argument in the dynamics to be optimized, designated as ${T}_{\text{f}}$ . Then, the following relationship between the virtual time and true time domain can be obtained for an arbitrary variable $\chi$

$\stackrel{˙}{\chi }\left(t\right)=\frac{d\chi \left(t\right)}{dt}=\frac{d\chi \left(t\right)}{{T}_{\text{f}}d\tau }={\chi }^{\prime }\left(\tau \right)/{T}_{\text{f}}.$       (20)

Especially, the derivative of the speed variable can be denoted as

$v\left(t\right)=\sqrt{\stackrel{˙}{x}{\left(t\right)}^{2}+\stackrel{˙}{y}{\left(t\right)}^{2}+\stackrel{˙}{h}{\left(t\right)}^{2}}=\sqrt{{\left({x}^{\prime }\right)}^{2}+{\left({y}^{\prime }\right)}^{2}+{\left({h}^{\prime }\right)}^{2}}/{T}_{\text{f}}=v\left(\tau \right)/{T}_{\text{f}}$       (21)

where the superscript () represents the derivative with respect to the virtual time.
According to Eqns (13) and (19), the dynamics can be rewritten as

(22)

Hence, the Problem 1 can be reformulated in virtual time domain (VTD) as Problem 2.

3.1.1 Problem 2 (CTP-OCP-VTD)

Minimize the cooperative cost function (Eqn (16)) of all UCAVs represented with respect to the new independent variable $\tau$ as:

$J=\sum _{j=1}^{{N}_{\text{v}}}{J}^{j}\left(x,u,{T}_{\text{f}}\right)=\sum _{j=1}^{{N}_{\text{v}}}\left[{w}_{\text{t}}^{j}{T}_{\text{f}}^{j}/\Gamma +{w}_{\text{r}}^{j}{\int }_{0}^{1}{\text{PRD}}^{j}\left(\tau \right)d\tau /{T}_{\text{f}}^{j}+{w}_{\text{d}}^{j}{f}_{dest}\left({x}^{j}\left(1\right)\right)\right]$       (23)

subject to the dynamics in Eqn (22), and the boundary constraints in Eqn (17), written as:

$\begin{array}{l}{\stackrel{˜}{\Phi }}_{0}\left[{x}^{j}\left(0\right),{u}^{j}\left(0\right)\right]=0\\ {\stackrel{˜}{\Phi }}_{1}\left[{x}^{j}\left(1\right),{u}^{j}\left(1\right),{T}_{\text{f}}^{j}\right]\le \text{0}\end{array}$       (24)

and the inequality and equality constraints (Eqn (18)), and additional temporal constraints

(25)

3.2  Problem Formulation in the Output Space

Due to the complexity in solving this higher dimensional space system, the differential flatness approach is introduced to transform the system dynamics to a lower dimensional space12. To confirm that the dynamics system is differentially flat, the spatial trajectory $\rho \left(\tau \right)$ in virtual time domain and the terminal time ${T}_{\text{f}}$ as the flat output vector can be defined as

$z=\left[\rho \left(\tau \right),{T}_{\text{f}}\right]={\left[x\left(\tau \right),y\left(\tau \right),h\left(\tau \right),{T}_{\text{f}}\right]}^{T}$       (26)

The original state x and control u can be recovered from the flat outputs and their derivatives as follows13

(27)

According to the kinematics equations in Eqn (14), the remaining three states using the flatness outputs in virtual time domain can be easily described as below

(28)

Thus, according to Eqn (15), the control variables can be determined as

(29)

Obviously, the dynamics constraints of this system (i.e. Eqns (14-15)) can be automatically satisfied. And the system of cooperative planning of UCAVs, including the objective function and the constraints, is mapped to a lower dimensional output space. Hence, Problem 2 can be modified as problem 3.

3.2.1 Problem 3 (CTP-OCP-VTD in Output Space)
Solving the problem

$\begin{array}{l}\mathrm{min}J=\sum _{j=1}^{{N}_{\text{v}}}{J}^{j}\left(\xi \left(\tau \right)\right)=\sum _{j=1}^{{N}_{\text{v}}}\left[{\varphi }^{j}\left(\phi \left({\xi }^{j}\left(0\right)\right),\alpha \left({\xi }^{j}\left(0\right)\right),\phi \left({\xi }^{j}\left(1\right)\right),\alpha \left({\xi }^{j}\left(1\right)\right),{T}_{\text{f}}^{j}\right)+\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\int }_{0}^{1}{L}^{j}\left(\phi \left({\xi }^{j}\left(\tau \right)\right),\alpha \left({\xi }^{j}\left(\tau \right)\right),{T}_{\text{f}}^{j}\right)d\tau \right)\right]\end{array}$       (30)

subject to

(31)

3.3  Parameterization of the Spatial Trajectory

To generate feasible trajectories in a finite parameter space, flatness outputs are parameterized in terms of B-spline curves, which have been used to represent the trajectory of UCAV to minimize computation loads and successfully applied to path/trajectory planning14. Without loss of generality, to decrease the computational time, the 3-degree B-spline curves are chosen to represent the trajectory, whose order equals to 4. Therefore, choosing every four neighboring points from the point set as control points, a 3-degree B-spline curve can be defined, which is a trajectory segment. In addition, it is clearly seen from Eqn (19) that the non-decreasing variable $\tau$ can be the knot sequence of B-spline. Consequently, the ith trajectory segment can be described as follows

${B}_{i}\left(\tau \right)={b}_{1,4}^{i}\left(\tau \right){C}_{1}^{i}+{b}_{2,4}^{i}\left(\tau \right){C}_{2}^{i}+{b}_{3,4}^{i}\left(\tau \right){C}_{3}^{i}+{b}_{4,4}^{i}\left(\tau \right){C}_{4}^{i}$       (32)

Obviously, as is known, once the control point serial ${C}_{m}^{i}$ is determined, the ith trajectory segment can be generated by Eqn (32). Thus, the trajectory of UCAV is mapped into the control point serial of B-spline curves.

As is previously stated, the spatial trajectory segment ${B}_{i}\left(\tau \right)$ can be replaced by flatness output variables, ${\left[{x}_{i}\left(\tau \right),{y}_{i}\left(\tau \right),{h}_{i}\left(\tau \right)\right]}^{\text{T}}$ and simplified as

(33)

where ${a}_{i}^{m-1}$, ${b}_{i}^{m-1}$, ${c}_{i}^{m-1}$ are the evaluation coefficients determined by the coefficients of the B-spline basis functions.

Once all control points are obtained, the spatial trajectory, i.e., the flatness outputs will be generated. Such that, the original control variables and remaining state variables, as well as the objective function and the constraints, can be parameterized (i.e., discretized) by the control point serial. Thus, the trajectory will be mapped to control point serial, and the complex trajectory planning is transformed into a parameter optimization problem.

3.4  Transformation into a Nonlinear Programming Problem

After the flatness outputs have been parameterized in terms of B-spline curves, the coefficients of the B-spline basis functions can be found by nonlinear programming.

First of all, the virtual time horizon of UCAVj denoted by $\tau \in \left[0,1\right]$ need be uniformity discretized as the $\text{(}{N}_{\text{c}}^{j}\text{+1)}$ collocation points below

(34)

And the coordinates of collocation points can be written as , which can be calculated by the Eqn (33). Thus, the continuous-spline curves can be discretized in $\text{(}{N}_{\text{c}}^{j}\text{+1)}$ points, and the constraints on the B-spline curves will be enforced at collocation points.
For the 3-degree B-spline curve under consideration, B-spline coefficients of UCAVj for all position outputs can be denoted as

${\zeta }_{j}:=\left[\begin{array}{l}\left({C}_{1}^{1j},{C}_{2}^{1j},\dots ,{C}_{{p}_{1j}}^{1j},{C}_{1}^{2j},{C}_{2}^{2j},\\ \text{\hspace{0.17em}}\dots ,{C}_{{p}_{2j}}^{2j},\dots ,{C}_{1}^{qj},{C}_{2}^{qj},\dots ,{C}_{{p}_{qj}}^{qj}\right)\end{array}\right]\in {ℝ}^{{M}_{j}}$     (35)

where ${M}_{j}=\sum _{i=1}^{q}{p}_{ij}$ Pij is the number of the control points for the ith position output of UCAVj , and these coefficients are used as the decision variables in the NLP problem.
Thus, the multi-UCAV trajectory planning problem in optimal control framework is transformed into a NLP problem, i.e. CTP-NLP, given by Problem 4.

3.4.1 Problem 4 (CTP-NLP)
Solving the problem

(36)

where $\zeta$ is a decision vector of cooperative planning, described as $\zeta :=\left({\zeta }_{1},{\zeta }_{2},\dots ,{\zeta }_{{N}_{\text{v}}}\right)$, ${r}_{n}$ and ${s}_{n}$ are the total number of the inequality and equality constraints of all UCAVs, respectively.

Then the resulting CTP-NLP problem can be solved through well developed algorithms, such as the sequential quadratic programming. In the paper, the SNOPT software toolbox is used due to its advantages in solution effectiveness for the large-scale NLP problems15.

The basic ideas presented in this paper are illustrated in the following two scenarios. The common parameters of models in our simulations are listed in Table 1. .

Table 1. State and control constraints of UCAVs

The experimental test environment is a rectangle area of 30 × 40 km2, as shown in Fig 2. and Fig. 4 All the results presented below are generated using TOMLAB/SNOPT software toolbox on a 2.4-GHz Core 2 CPU, 2G RAM computer running with MATLAB R2009b. The weighting factors, ${w}_{\text{t}}^{j}$ , ${w}_{\text{r}}^{j}$ , and ${w}_{\text{d}}^{j}$ are set as 0.4, 0.3, and 0.3, by using trail and error. And the minimum safety radius of each UCAV dsafe is set as 500 m. The first initial guesses of the optimization parameters can be generated by using the B-spline interpolation of the lines from the IPs to the AARs of the targets. Then, the multiresolution-based iterative strategy8 is used to generate some ‘good’ initial guesses to the solver, which can reduce the number of iterations required to solve the NLP problem. The specified maximum number of collocation points is set as. ${N}_{\text{c,max}}=100$ Moreover, it is assumed that the target assignment is already completed before.

Figure 2. Two collision-free UCAV trajectories with arriving simultaneously.

4.1 Scenario 1: Cooperative Trajectories of Two UCAVs Arriving Simultaneously

In this scenario, two UCAVs cooperatively attack two stationary ground targets while avoiding a series of static obstacles/threats detected and collision en route, and satisfying aircraft dynamics constraint, especially simultaneous arrival constraint. UCAV1 and UCAV2 start at each IP, i.e., IP1 (10 km, 2 km, 2 km) and IP2 (17 km, 2 km, 2 km). Then they fly into the AARs of two targets: Target 1 (4.2 km, 34 km, 0 km) and Target 2 (14 km, 40 km, 0 km), respectively. The other three initial states and control inputs are given by and. Thus and . The collision-free trajectories of the UCAVs and their arrival time, i.e. the total flight time, are shown in Fig 2. In addition, the approximate weapon trajectories are drawn to simulate the attack process. The time histories of the UCAVs’ states ( ) and control inputs are shown in Fig. 3 The average detection-probabilities of UCAVs are 0.12 and 0.17, and the target damage probabilities of UCAVs are 0.851 and 0.762, respectively. It is clear that the resulting trajectories are smooth, and the constraints on these variables, especially the cooperative constraints are all satisfied (Table 1.), which means the resulting trajectories are feasible and safe.

Figure 3. State and control time histories of the two UCAVs.

4.2 Scenario 2: Cooperative Trajectories of Multi-UCAV Arriving in Sequence

In this scenario, three UCAVs attack two stationary ground targets cooperatively. The only additional requirement is that the UCAVs arrive at their AARs in sequence rather than simultaneously, and the intervals between UCAVs are equal, denoted as . The coordinates of IP1, IP2 and two targets are the same as Scenario 1, and the third IP is IP3 (4 km, 2 km, 2 km). The result of the previous finished target assignment is that the UCAV1 and UCAV2 attack Target 1 and UCAV3 attacks the other one. The initial remaining three states and control inputs of the UCAV3 are preset as and . Thus , and Fig. 4 shows the overall collision-free attack trajectories of multi-UCAV and their arrival time. It can be clearly demonstrated that the UCAVs can avoid all obstacles or threats and successfully fly into the AARs in sequence to perform weapon delivery. Fig. 5 shows the distance between each pair UCAVs. From it, one can find that the minimum distance is more than the minimum safety radius of dsafe = 500 m. Fig. 6 shows the time histories of the UCAVs’ states ( ) and control inputs . The average detection-probabilities of UCAVs are 0.133, 0.092 and 0.12, and the target damage probabilities of UCAVs are 0.813, 0.901 and 0.724, respectively. Obviously, the resulting trajectories are feasible and safe for all the constraints listed in Table 1. are satisfied.

Figure 4. Three collision-free UCAV trajectories with arriving in sequence.

Figure 5. Distance between each pair UCAVs.

Figure 6. State and control time histories of the three UCAVs.

This paper is devoted to exploring the cooperative collision-free trajectory planning for multi-UCAV performing the CA/GTA missions. A novel cooperative planning approach is proposed in optimal control framework, based on differential flatness, B-spline curves and nonlinear programming. To integrate weapon delivery constraints in the problem formulation, an approximate AAR model is established. Moreover, the notion of the virtual time is introduced to handle the temporal constraint. Instead of solving an OCP over a high-dimensional continuous space, the NLP problem of very low dimension has been solved successfully, over a much smaller space. The proposed approach can efficiently solve the point-to-region trajectory planning problem with integrating the spatial and temporal constraints on the trajectory level, whose validity is illustrated with some simulation results finally. Further efforts may be made to analyse some uncertain factors in the true battlefield environment and carry out the research on the real-time cooperative trajectory planning.

1.     Murray, R.M. Recent research in cooperative control of multi-vehicle systems. J. Dyn. Syst.-T. ASME, 2007, 129(5), 571-583.

2.     Hurni, M. A.; Sekhavat, P.; Karpenko, M. & Ross, I. M. A pseudospectral optimal motion planner for autonomous unmanned vehicles. In American Control Conference. Baltimore, MD, USA. 2010.[Full text via CrossRef]

3.     McLain, T.W. & Beard, R.W. Coordination variables, coordination functions, and cooperative-timing missions. J. Guid. Control Dynam., 2005, 28(1), 150-161.[Full text via CrossRef]

4.     Kaminer, I.; Yakimenko, O.; Dobrokhodov, V.; A. Pascoal; Hovakimyan, N.; Cao, C.; Young, A. & Patel, V. Coordinated path pollowing for time-critical missions of multiple UAVs via L1 adaptive output feedback controllers. In AIAA Guidance, Navigation and Control Conference and Exhibit. Hilton Head, South Carolina, USA. 2007, pp.1-23.[Full text via CrossRef]

5.     Bollino, K.P. & Lewis, L.R. Collision-free multi-UAV optimal path planning and cooperative control for tactical applications. In AIAA Guidance, Navigation and Control Conference and Exhibit. Honolulu, Hawaii, USA. 2008, pp.1-18.[Full text via CrossRef]

6.     Lian, F.-L. Cooperative path planning of dynamical multi-agent systems using differential flatness approach. Int. J. Control Autom., 2008, 6(3), 401-412.

7.     Kuwata, Y. & How, J.P. Cooperative distributed robust trajectory optimization using receding horizon MILP. IEEE T. Contr. Syst. T., 2011, 19(2), 423-431.[Full text via CrossRef]

8.     Zhang, Y.; Chen, J. & Shen, L. Hybrid hierarchical trajectory planning for a fixed-wing UCAV performing air-to-surface multi-target attack. J. Syst. Eng. Electron., 2012, 23(4), 256-264.[Full text via CrossRef]

9.     Besada-Portas, E.; Torre, L.d.l.; Cruz, J.M.d.l. & Andres-Toro, B.d. Evolutionary trajectory planner for multiple UAVs in realistic scenarios. IEEE Trans.. Robot., 2010, 26(4), 619-634.[Full text via CrossRef]

10.   Stevens, B.L. & Lewis, F.L. Aircraft control and simulation. John Wiley & Sons, Inc. Ed 2, New Jersey, 2003. 680 p.

11.   Lewis, L.-P.R. Rapid motion planning and autonomous obstacle avoidance for unmanned vehicles. Naval Postgraduate School, Monterey, California, USA, 2006. MA Thesis.

12.   Michel, F. Flatness and defect of nonlinear systems: Introductory theory and examples. Int. J. Control, 1995, 61(6), 1327-1361.[Full text via CrossRef]

13.   Zhang, Y.; Chen, J. & Shen, L. Real-time trajectory planning for UCAV air-to-surface attack using inverse dynamics optimization method and receding horizon control. Chinese J. Aeronaut., 2013, 26(4), 1038-1056.[Full text via CrossRef]

14.   Foo, J.L.; Knutzon, J.; Kalivarapu, V.; Oliver, J. & Winer, E. Path planning of unmanned aerial vehicles using B-splines and particle swarm optimization. J. Aeros. Comp. Inf. Com., 2009, 6, 271-290.[Full text via CrossRef]

15.   E. Rutquist, P. & M. Edvall, M. PROPT - Matlab optimal control software. TOMLABoptimization. 2010.

 Dr Xueqiang Gu is a PhD scholar in College of Mechatronic Engineering and Automation, National University of Defense Technology (NUDT), Changsha, Hunan, China. His current research interests included: combat vehicle mission planning and intelligence control. Dr Yu Zhang received his PhD in control science and engineering from NUDT, Changsha, Hunan, China in 2012. He is a Docent of College of Mechatronic Engineering and Automation in NUDT. His current research interests included: combat vehicle mission planning and artificial intelligence Dr Jing Chen obtained his PhD in control science and engineering from NUDT in 1999. He is a Professor of College of Mechatronic Engineering and Automation in NUDT, Changsha, Hunan, China. His current research interests included: artificial intelligence and mission planning of aircraft Mr Lincheng Shen is the Dean and Professor of College of Mechatronic Engineering and Automation in NUDT, Changsha, Hunan, China. He has published over 100 technical papers in refereed international journals and academic conferences. His current research interests included: mission planning, SAR image processing, biomimetic robotics, automation and control engineering.