# A Virtual Motion Camouflage Approach for Cooperative Trajectory Planning of Multiple UCAVs

Authors: Xueqiang Gu; Lincheng Shen; Jing Chen; Yu Zhang; Wanpeng Zhang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper investigates cooperative trajectory planning of multiple unmanned combat aerial vehicles (multi-UCAV) in performing autonomous cooperative air-to-ground target attack missions. By integrating an approximate allowable attack region model, several constraint models, and a multicriterion objective function, the problem is formulated as a cooperative trajectory optimal control problem (CTOCP). Then, a virtual motion camouflage (VMC) for cooperative trajectory planning of multi-UCAV, combining with the differential flatness theory, Gauss pseudospectral method (GPM), and nonlinear programming, is designed to solve the CTOCP. In particular, the notion of the virtual time is introduced to the VMC problem formulation to handle the temporal cooperative constraints. The simulation experiments validate that the CTOCP can be effectively solved by the cooperative trajectory planning algorithm based on VMC which integrates the spatial and temporal constraints on the trajectory level, and the comparative experiments illustrate that VMC based algorithm is more efficient than GPM based direct collocation method in tackling the CTOCP.


---

## Body

## 1. Introduction

Nowadays, it is an active research area to perform autonomous cooperative air-to-ground target attack (CA/GTA) missions using multiple unmanned combat aerial vehicles (multi-UCAV) [1]. However, compared with single UCAV planning and coordinated formation control problems [2], new technical challenges in the CA/GTA missions are emerging. The cooperative trajectory planning is 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.Recently, various algorithms have been developed to solve this cooperative trajectory planning problem [3, 4], including artificial neural network methods [5], sample-based planning methods [6], maneuver automation (MA) [7], and optimal control methods. There is no doubt that the optimal control theory is the most natural framework for this type of problem with dynamic constraints [8]. However, the rapid solution to optimal control problems (OCPs) for complicated nonlinear systems, such as UCAVs, is a challenging task [9]. Analytical solutions are seldom available or even possible. As a result, one usually resorts to numerical techniques [3]. The techniques can be classified into two general types, namely, indirect and direct methods. Indirect methods [10] solve the OCPs by formulating the first-order optimality conditions, applying Pontryagin’s Maximum Principle and nonlinear programming (NLP) to tackle the resulting two-point boundary value problem numerically. Direct methods [11], on the other hand, are devoted to reduce the OCPs to finite-dimensional NLP problems by discretization and parameterization of subsets of the state and control vectors and then are solved by developed optimizers. As one of direct methods, Gauss pseudospectral method (GPM) is applied widely and efficiently in the trajectory optimization field [12–14], due to its advantages of fewer parameters and higher precision in the calculation procedure. However, the method to achieve the trajectory planning for a single vehicle or cooperative multivehicle required a high computational load. To reduce the computational complexity, the inverse dynamics methods [15, 16] and differential flatness theory (DFT) based methods [17–20] are introduced. Compared to pseudospectral methods, these methods can use any models and any performance indexes and transform the OCP into a low dimensional NLP problem. However, they are still time-consuming. Particularly, when the number of the vehicles is too large, the computational cost will be unacceptable. For this purpose, the virtual motion camouflage (VMC) algorithm [21–23] is introduced into this work. Inspired by the biological motion known as the motion camouflage, the VMC algorithm can dramatically reduce the problem dimension and the computational cost. Some numerical simulations [22–24] suggest that its computational speed could be much faster than the pseudospectral methods. Motivated by these advantages, this paper employs virtual motion camouflage approach to develop cooperative trajectory planning algorithms.For the cooperative trajectory planning of multi-UCAV, the objective is to generate dynamically feasible trajectories for UCAVs, which can guide them to the goals in the shortest time or distance, without collisions with obstacles or each other. To date, a number of theories and techniques have been developed to accomplish the cooperative tasks, including several collision avoidance techniques and time adjustment strategies [17, 25–30]. Bollino and Lewis [25] addressed the optimal path planning of UAVs in obstacle-rich environments and proposed a collision-free path planning for multi-UAV using optimal control techniques and pseudospectral methods. Kuwata and How [26] presented a cooperative distributed robust trajectory optimization approach, using RH-MILP with independent dynamics but coupled objectives and hard constraints. To improve the convergence rate, the virtual motion camouflage method was applied to the cooperative electronic combat air vehicles (ECAVs) [27] and unmanned aerial vehicle (UAV) formation control [28]. Reference [27] described how an interesting bioinspired motion strategy can be used to design real-time trajectories for cooperative electronic combat air vehicles. The constant speed motion camouflage law was developed to derive the feasible condition of a constant speed ECAVs and PT coherent mission, and its feasibility conditions were found. Reference [28] proposed a divide and conquer hierarchical approach in three levels to solve the UAV formation flight trajectory plan problem considering dynamics, state, and control variable inequality and equality constraints. These approaches mentioned above could generate collision-free trajectories, but the temporal cooperation was ignored. In contrast, Lian [17] introduced a differential flatness theory (DFT) based method to optimally formulate the cooperative path planning for multiagent dynamical systems considering spatial and temporal constraints. McLain and Beard  [29] proposed the coordination variables (CV) and coordination function (CF) based strategy to achieve cooperative timing among teams of vehicles by coordinating the velocity and path length of each vehicle. But these approaches failed to deal with the dynamics constraints of vehicles, and the generated paths were not always flyable and smooth. Kaminer et al. [30] presented a general framework for coordinated control problem of multiple autonomous agents. On the basis of decoupling space and time in the problem formulation, it reduced the number of optimization parameters and made it easy to implement optimization in real time. However, the time coordination lays in the path following and the design of control laws in path following algorithms was much more difficult. Although the previous investigations have described several valuable strategies in the cooperative path planning, these methods cannot tackle the point-to-region cooperative trajectory planning for CA/GTA missions directly, which needs to integrate both the spatial and temporal constraints on the trajectory level.To address the problems mentioned previously, a novel cooperative trajectory planning method for multi-UCAV in performing the CA/GTA missions is presented in this paper. Firstly, some constraints including individual and cooperative constraints are modeled. Particularly, an approximate allowable attack region is built for the critical terminal constraints. Then, after the multicriterion objective function is constructed, the cooperative trajectory planning problem is formulated as a cooperative trajectory optimal control problem (CTOCP). Owing to the temporal constraints, a notion of the virtual time is introduced and the VMC problem is reformulated in the virtual time domain. Inspired by VMC, DFT, and GPM algorithms, a new cooperative trajectory planning algorithm in an optimal control framework is proposed. The proposed approach is demonstrated by two typical CA/GTA examples, and the comparative experiments between GPM and VMC are carried out. The results show that the proposed approach is feasible, effective, and efficient.The rest of the paper is organized as follows. In Section2, the integrated model is set up. In Section 3, the problem is formulated as a CTOCP followed by a multicriteria objective function. Section 4 develops a novel cooperative trajectory planning algorithm based on virtual motion camouflage in virtual time domain, which can solve the problem efficiently. Section 5 presents several numerical examples, and finally the conclusions and future works are outlined in Section 6.


## 2. Modeling

### 2.1. Aircraft Model

The kinematic and dynamic model of vehicles is needed in the cooperative trajectory planning of multi-UCAV. In this work, a team ofN
v homogeneous UCAVs is considered, each of which is described by the kinematic and dynamic model according to a full-blown three-degree of freedom (3-DOF) model as follows [31]:
(1)
x
˙
=
V
cos
⁡
⁡
γ
cos
⁡
⁡
ψ
,
y
˙
=
V
cos
⁡
⁡
γ
sin
⁡
ψ
,
z
˙
=
V
sin
⁡
γ
,
V
˙
=
g
(
n
x
-
sin
⁡
γ
)
,
γ
˙
=
g
V
(
n
z
cos
⁡
⁡
μ
-
cos
⁡
⁡
γ
)
,
h
h
h
h
h
h
h
h
h
h
h
ψ
˙
=
g
V
cos
⁡
⁡
γ
n
z
sin
⁡
μ
,
where x, y
, and z are the aircraft coordinates, that is, longitude, latitude, and height, V is the aircraft velocity, γ is the flight-path angle, ψ is the heading angle, μ is the roll angle, g is the gravity acceleration, and n
x and n
z are the longitudinal and normal components of the load factor, respectively.For the optimal control problem, the state vectorx and control vector u can be defined as
(2)
x
=
{
(
x
,
y
,
z
,
V
,
γ
,
ψ
)
∣
g
z
min
⁡
≤
z
≤
z
max
⁡
g
V
min
⁡
≤
V
≤
V
max
⁡
g
γ
min
⁡
≤
γ
≤
γ
max
⁡
ψ
min
⁡
≤
ψ
≤
ψ
max
⁡
}
∈
X
⊂
ℝ
6
,
u
=
{
(
μ
,
n
x
,
n
z
)
∣
g
g
i
μ
min
⁡
≤
μ
≤
μ
max
⁡
n
x
,
min
⁡
≤
n
x
≤
n
x
,
max
⁡
n
z
,
min
⁡
≤
n
z
≤
n
z
,
max
⁡
}
∈
U
⊂
ℝ
3
.


### 2.2. Allowable Attack Region Model

For the point-to-region trajectory planning problem, the allowable attack regions (AARs) of targets are defined as the areas where weapon delivery operations can be effectively performed by UCAVs. Therefore, in order to plan accurate and optimal attack trajectories, the AARs and the delivery parameters need to be obtained. The AARs of thei
th target TA
R
i, denoted as R
(
TA
R
i
), are such sets of all feasible release states that TA
R
i can be effectively attacked whenever the aircraft belongs to those state sets. When n
v (n
v
≥
2) UCAVs are assigned to attack the same target TA
R
i cooperatively, there will be n
v AARs in R
(
TA
R
i
); that is, R
(
TA
R
i
)
=
⋃
j
n
v
AA
R
j. Each AA
R
j can be formulated as an abstract 6-dimensional space [32]:
(3)
AA
R
j
=
{
x
,
y
,
z
,
V
,
γ
,
ψ
}
⊂
ℝ
6
.Obviously,AA
R
j is a high-dimensional nonlinear space, which is difficult to handle. By presetting an appropriate weapon release speed V
⌢
r and the flight-path angle γ
⌢
r based on estimating the weapon impact effects and destruction requirements to the target and predetermining release heading ψ
⌢
r, AA
R
j can be reduced to a 3-dimensional space as(4)
LA
R
sp
j
(
TA
R
i
,
[
V
⌢
r
,
γ
⌢
r
]
)
=
{
(
x
,
y
,
z
,
V
⌢
r
,
γ
⌢
r
,
ψ
⌢
r
)
∣
(
x
,
y
,
z
)
∈
LA
R
sp
j
(
TA
R
i
,
[
V
⌢
r
,
γ
⌢
r
,
ψ
⌢
r
]
)
ψ
⌢
r
=
azimuth
of
TA
R
i
relative
to
position
(
x
,
y
,
z
)
}
.In the point-to-region trajectory planning, weapon delivery points (WDPt) as terminals of trajectories need to be included in the AARs, that is, meeting terminal constraints [33]. The formula can be denoted as WDP
t
j
∈
AA
R
j; that is,
(5)
|
x
A
j
-
x
f
j
|
≤
Δ
x
,
|
y
A
j
-
y
f
j
|
≤
Δ
y
,
h
h
h
h
h
h
h
h
|
z
A
j
-
z
f
j
|
≤
Δ
z
,
where (
x
A
j
,
y
A
j
,
z
A
j
) are the coordinates of the center of the AA
R
j, (
x
f
j
,
y
f
j
,
z
f
j
) are the coordinates of WDP
t
j, and (
Δ
x
,
Δ
y
,
Δ
z
) are the thresholds of errors.


### 2.3. No-Fly Zone Model

In the battlefield environment, the no-fly zone (NFZ) is an area where vehicles are not permitted to fly over, due to the presence of military restrictions (e.g., armed enemies and the missile killing range), the physical obstacles (e.g., mountains and buildings within the natural and urban environment, resp.), and civil restrictions mainly due to safety reasons (e.g., densely populated areas and severe weather condition zones).It is complex and unnecessary to describe a NFZ’s exact shape and size. In this work, thep-norm method [34] is used to mathematically model the shapes of the NFZs, given by
(6)
h
(
x
,
y
,
z
)
=
∥
(
x
-
x
c
k
1
+
b
s
,
y
-
y
c
k
2
+
b
s
,
z
-
z
c
k
3
+
b
s
)
∥
p
p
-
k
4
p
,
(7)
∥
(
x
,
y
,
z
)
∥
p
=
(
|
x
|
p
+
|
y
|
p
+
|
z
|
p
)
1
/
p
,
p
∈
N
,
where (x
c
,
y
c
,
z
c) indicates the location of the geometric center of the NFZ, h
(
x
,
y
,
z
) is the distance between a point (x
,
y
,
z) and the boundary of the NFZ, while k
1
~
k
4 are the constant parameters chosen to define the size of the NFZ, and b
s represents the width of a safe buffer which can effectively expand the NFZ boundary by an amount that accounts for the size of the vehicle and improve the robustness of plans. By varying parameters (p
,
k
1
~
k
4), one can easily model a number of generic shapes (see Figure 1). And any NFZ can be approximately modeled by fitting one of those shapes around it.Figure 1
Several shapes generated from (6).Accordingly, the above NFZs can be expressed by the following path constraints in the OCP formulation as(8)
h
n
(
x
,
y
,
z
)
≥
0
,
n
=
1,2
,
…
,
N
NFZ
,
where h
n represents the n
th NFZ and N
NFZ is the total of NFZs. However, (8) is not well scaled. If a large p is chosen or the vehicle is far from the NFZ, h
n
(
x
,
y
,
z
) can produce very large numbers; hence, for computational efficiency, the path constraint is scaled by the natural logarithm function as follows:
(9)
h
n
(
x
,
y
,
z
)
=
ln
⁡
(
1
(
k
4
n
)
p
n
(
(
x
-
x
c
n
k
1
n
+
b
s
)
p
n
+
(
y
-
y
c
n
k
2
n
+
b
s
)
p
n
h
h
h
h
h
h
h
h
h
h
+
(
z
-
z
c
n
k
3
n
+
b
s
)
p
n
)
)
≥
0
.


### 2.4. Cooperative Constraint Model

In order to make the UCAVs arrive in the AARs of targets safely in an expected time sequence, cooperative trajectories should satisfy the following constraints.


#### 2.4.1. Spatial Constraints

During the mission, multi-UCAV should maintain a safe distance to guarantee them not to collide with each other. The model can be denoted as(10)
∥
ρ
j
(
t
i
)
-
ρ
k
(
t
i
)
∥
2
≥
max
⁡
(
d
safe
j
,
d
safe
k
)
,
h
h
h
h
h
h
h
i
h
h
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
,
where ρ
j
(
t
i
)
=
{
x
j
(
t
i
)
,
y
j
(
t
i
)
,
z
j
(
t
i
)
} is the spatial position of each UCAV
j at the time t
i, d
safe
j is the predetermined minimum safety radius of UCAV
j(based on the vehicle size, maneuver agility, etc.), and N
v is the total of UCAVs.


#### 2.4.2. Temporal Constraints

The temporal constraints include the simultaneous arrival constraint and sequencing constraint, which are described as(11)
t
f
j
+
Δ
j
k
-
t
f
k
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
,
where t
f
j is the terminal time of each UCAV
j and Δ
j
k is the arrival interval between UCAV
j and UCAV
k. For the simultaneous arrival constraint, Δ
=
0, and for the sequencing constraint, Δ is the plus fixed quantity. Equation (11) is of the general form f
T
(
t
f
1
,
…
,
t
f
N
v
)
≤
0.


## 2.1. Aircraft Model

The kinematic and dynamic model of vehicles is needed in the cooperative trajectory planning of multi-UCAV. In this work, a team ofN
v homogeneous UCAVs is considered, each of which is described by the kinematic and dynamic model according to a full-blown three-degree of freedom (3-DOF) model as follows [31]:
(1)
x
˙
=
V
cos
⁡
⁡
γ
cos
⁡
⁡
ψ
,
y
˙
=
V
cos
⁡
⁡
γ
sin
⁡
ψ
,
z
˙
=
V
sin
⁡
γ
,
V
˙
=
g
(
n
x
-
sin
⁡
γ
)
,
γ
˙
=
g
V
(
n
z
cos
⁡
⁡
μ
-
cos
⁡
⁡
γ
)
,
h
h
h
h
h
h
h
h
h
h
h
ψ
˙
=
g
V
cos
⁡
⁡
γ
n
z
sin
⁡
μ
,
where x, y
, and z are the aircraft coordinates, that is, longitude, latitude, and height, V is the aircraft velocity, γ is the flight-path angle, ψ is the heading angle, μ is the roll angle, g is the gravity acceleration, and n
x and n
z are the longitudinal and normal components of the load factor, respectively.For the optimal control problem, the state vectorx and control vector u can be defined as
(2)
x
=
{
(
x
,
y
,
z
,
V
,
γ
,
ψ
)
∣
g
z
min
⁡
≤
z
≤
z
max
⁡
g
V
min
⁡
≤
V
≤
V
max
⁡
g
γ
min
⁡
≤
γ
≤
γ
max
⁡
ψ
min
⁡
≤
ψ
≤
ψ
max
⁡
}
∈
X
⊂
ℝ
6
,
u
=
{
(
μ
,
n
x
,
n
z
)
∣
g
g
i
μ
min
⁡
≤
μ
≤
μ
max
⁡
n
x
,
min
⁡
≤
n
x
≤
n
x
,
max
⁡
n
z
,
min
⁡
≤
n
z
≤
n
z
,
max
⁡
}
∈
U
⊂
ℝ
3
.


## 2.2. Allowable Attack Region Model

For the point-to-region trajectory planning problem, the allowable attack regions (AARs) of targets are defined as the areas where weapon delivery operations can be effectively performed by UCAVs. Therefore, in order to plan accurate and optimal attack trajectories, the AARs and the delivery parameters need to be obtained. The AARs of thei
th target TA
R
i, denoted as R
(
TA
R
i
), are such sets of all feasible release states that TA
R
i can be effectively attacked whenever the aircraft belongs to those state sets. When n
v (n
v
≥
2) UCAVs are assigned to attack the same target TA
R
i cooperatively, there will be n
v AARs in R
(
TA
R
i
); that is, R
(
TA
R
i
)
=
⋃
j
n
v
AA
R
j. Each AA
R
j can be formulated as an abstract 6-dimensional space [32]:
(3)
AA
R
j
=
{
x
,
y
,
z
,
V
,
γ
,
ψ
}
⊂
ℝ
6
.Obviously,AA
R
j is a high-dimensional nonlinear space, which is difficult to handle. By presetting an appropriate weapon release speed V
⌢
r and the flight-path angle γ
⌢
r based on estimating the weapon impact effects and destruction requirements to the target and predetermining release heading ψ
⌢
r, AA
R
j can be reduced to a 3-dimensional space as(4)
LA
R
sp
j
(
TA
R
i
,
[
V
⌢
r
,
γ
⌢
r
]
)
=
{
(
x
,
y
,
z
,
V
⌢
r
,
γ
⌢
r
,
ψ
⌢
r
)
∣
(
x
,
y
,
z
)
∈
LA
R
sp
j
(
TA
R
i
,
[
V
⌢
r
,
γ
⌢
r
,
ψ
⌢
r
]
)
ψ
⌢
r
=
azimuth
of
TA
R
i
relative
to
position
(
x
,
y
,
z
)
}
.In the point-to-region trajectory planning, weapon delivery points (WDPt) as terminals of trajectories need to be included in the AARs, that is, meeting terminal constraints [33]. The formula can be denoted as WDP
t
j
∈
AA
R
j; that is,
(5)
|
x
A
j
-
x
f
j
|
≤
Δ
x
,
|
y
A
j
-
y
f
j
|
≤
Δ
y
,
h
h
h
h
h
h
h
h
|
z
A
j
-
z
f
j
|
≤
Δ
z
,
where (
x
A
j
,
y
A
j
,
z
A
j
) are the coordinates of the center of the AA
R
j, (
x
f
j
,
y
f
j
,
z
f
j
) are the coordinates of WDP
t
j, and (
Δ
x
,
Δ
y
,
Δ
z
) are the thresholds of errors.


## 2.3. No-Fly Zone Model

In the battlefield environment, the no-fly zone (NFZ) is an area where vehicles are not permitted to fly over, due to the presence of military restrictions (e.g., armed enemies and the missile killing range), the physical obstacles (e.g., mountains and buildings within the natural and urban environment, resp.), and civil restrictions mainly due to safety reasons (e.g., densely populated areas and severe weather condition zones).It is complex and unnecessary to describe a NFZ’s exact shape and size. In this work, thep-norm method [34] is used to mathematically model the shapes of the NFZs, given by
(6)
h
(
x
,
y
,
z
)
=
∥
(
x
-
x
c
k
1
+
b
s
,
y
-
y
c
k
2
+
b
s
,
z
-
z
c
k
3
+
b
s
)
∥
p
p
-
k
4
p
,
(7)
∥
(
x
,
y
,
z
)
∥
p
=
(
|
x
|
p
+
|
y
|
p
+
|
z
|
p
)
1
/
p
,
p
∈
N
,
where (x
c
,
y
c
,
z
c) indicates the location of the geometric center of the NFZ, h
(
x
,
y
,
z
) is the distance between a point (x
,
y
,
z) and the boundary of the NFZ, while k
1
~
k
4 are the constant parameters chosen to define the size of the NFZ, and b
s represents the width of a safe buffer which can effectively expand the NFZ boundary by an amount that accounts for the size of the vehicle and improve the robustness of plans. By varying parameters (p
,
k
1
~
k
4), one can easily model a number of generic shapes (see Figure 1). And any NFZ can be approximately modeled by fitting one of those shapes around it.Figure 1
Several shapes generated from (6).Accordingly, the above NFZs can be expressed by the following path constraints in the OCP formulation as(8)
h
n
(
x
,
y
,
z
)
≥
0
,
n
=
1,2
,
…
,
N
NFZ
,
where h
n represents the n
th NFZ and N
NFZ is the total of NFZs. However, (8) is not well scaled. If a large p is chosen or the vehicle is far from the NFZ, h
n
(
x
,
y
,
z
) can produce very large numbers; hence, for computational efficiency, the path constraint is scaled by the natural logarithm function as follows:
(9)
h
n
(
x
,
y
,
z
)
=
ln
⁡
(
1
(
k
4
n
)
p
n
(
(
x
-
x
c
n
k
1
n
+
b
s
)
p
n
+
(
y
-
y
c
n
k
2
n
+
b
s
)
p
n
h
h
h
h
h
h
h
h
h
h
+
(
z
-
z
c
n
k
3
n
+
b
s
)
p
n
)
)
≥
0
.


## 2.4. Cooperative Constraint Model

In order to make the UCAVs arrive in the AARs of targets safely in an expected time sequence, cooperative trajectories should satisfy the following constraints.


### 2.4.1. Spatial Constraints

During the mission, multi-UCAV should maintain a safe distance to guarantee them not to collide with each other. The model can be denoted as(10)
∥
ρ
j
(
t
i
)
-
ρ
k
(
t
i
)
∥
2
≥
max
⁡
(
d
safe
j
,
d
safe
k
)
,
h
h
h
h
h
h
h
i
h
h
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
,
where ρ
j
(
t
i
)
=
{
x
j
(
t
i
)
,
y
j
(
t
i
)
,
z
j
(
t
i
)
} is the spatial position of each UCAV
j at the time t
i, d
safe
j is the predetermined minimum safety radius of UCAV
j(based on the vehicle size, maneuver agility, etc.), and N
v is the total of UCAVs.


### 2.4.2. Temporal Constraints

The temporal constraints include the simultaneous arrival constraint and sequencing constraint, which are described as(11)
t
f
j
+
Δ
j
k
-
t
f
k
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
,
where t
f
j is the terminal time of each UCAV
j and Δ
j
k is the arrival interval between UCAV
j and UCAV
k. For the simultaneous arrival constraint, Δ
=
0, and for the sequencing constraint, Δ is the plus fixed quantity. Equation (11) is of the general form f
T
(
t
f
1
,
…
,
t
f
N
v
)
≤
0.


## 2.4.1. Spatial Constraints

During the mission, multi-UCAV should maintain a safe distance to guarantee them not to collide with each other. The model can be denoted as(10)
∥
ρ
j
(
t
i
)
-
ρ
k
(
t
i
)
∥
2
≥
max
⁡
(
d
safe
j
,
d
safe
k
)
,
h
h
h
h
h
h
h
i
h
h
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
,
where ρ
j
(
t
i
)
=
{
x
j
(
t
i
)
,
y
j
(
t
i
)
,
z
j
(
t
i
)
} is the spatial position of each UCAV
j at the time t
i, d
safe
j is the predetermined minimum safety radius of UCAV
j(based on the vehicle size, maneuver agility, etc.), and N
v is the total of UCAVs.


## 2.4.2. Temporal Constraints

The temporal constraints include the simultaneous arrival constraint and sequencing constraint, which are described as(11)
t
f
j
+
Δ
j
k
-
t
f
k
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
,
where t
f
j is the terminal time of each UCAV
j and Δ
j
k is the arrival interval between UCAV
j and UCAV
k. For the simultaneous arrival constraint, Δ
=
0, and for the sequencing constraint, Δ is the plus fixed quantity. Equation (11) is of the general form f
T
(
t
f
1
,
…
,
t
f
N
v
)
≤
0.


## 3. Problem Formulation

### 3.1. Objective Function

The objective function of the entire team is a sum of each individual cost, which is a multicriterion objective function, and each criterion could compete with each other as follows:(12)
J
=
min
⁡
⁡
∑
j
=
1
N
v
J
j
=
min
⁡
⁡
∑
j
=
1
N
v
(
w
t
j
J
t
j
+
w
p
j
J
prd
j
+
w
r
j
J
rob
j
)
,
and the criterions can be, respectively, defined as
(13)
J
t
j
=
∫
t
d
t
=
t
f
j
-
t
0
j
,
(14)
J
prd
j
=
∫
t
PR
D
j
(
t
)
=
∫
t
(
1
-
∏
r
=
1
n
r
(
1
-
P
d
j
(
t
,
r
)
)
)
,
(15)
J
rob
j
=
∫
t
R
j
(
t
)
=
∫
t
∑
n
=
1
N
NFZ
w
r
n
,
j
(
e
e
-
h
n
-
1
)
,
where t
0
j and t
f
j denote the initial and terminal time of UCAV
j. The first term, J
t
j, denotes the total fight time of UCAV
j. The second term, J
prd
j, denotes the detection-probability of the n
r-radar system to UCAV
j at the time t, which describes the threat risk criterion of UCAV
j as follows:
(16)
PR
D
j
(
t
)
=
1
-
∏
r
=
1
n
r
(
1
-
P
d
j
(
t
,
r
)
)
,
where P
d
j
(
t
,
r
) is the model of the radar detection probability between the trajectory point of UCAV
j at the time t and the r
th radar [35]. The final term, J
rob
j, guarantees the robustness of the vehicle collision avoidance against the model imprecision, disturbing in the battlefield and operating errors that could otherwise cause vehicles to collide with obstacles or other vehicles. The robustness cost is defined as
(17)
R
j
(
t
)
=
∑
n
=
1
N
NFZ
w
r
n
,
j
(
e
e
-
h
n
-
1
)
,
where w
r
n
,
j are weights that are incremented to increase the maneuver robustness when the trajectory of a vehicle passes in close proximity to an obstacle or another vehicle, which can be chosen by the method presented by Hurni et al. [36].As can be seen, the individual objective functionJ
j is defined by the weighted sum of the three separate running cost terms with appropriate weighting factors w
t
j, w
p
j, and w
r
j. The three criteria (13)–(15) represent different physical meanings, and they are difficult to be compared directly. Hence the selection of these weighting factors is a skilled technique. In the experimentations, the multi-objective fuzzy optimization method proposed by Wang et al. [37] is employed. The method includes two main steps: (a) normalizing each single objective function and (b) solving the membership function about each criterion by using fuzzy distribution function.


### 3.2. Cooperative Trajectory Planning Problem Formulation

After establishing the above models, the cooperative trajectory planning problem is formulated in this section. The problem under consideration is a cooperative scheme, consisting ofN
v homogeneous UCAVs. According to the system equations (1), the general form of the system for UCAV
j is given by
(18)
x
˙
j
(
t
)
=
f
j
[
x
j
(
t
)
,
u
j
(
t
)
,
t
]
,
j
=
1
,
…
,
N
v
.As stated above, to obtain optimal or suboptimal cooperative trajectories, the cooperative trajectory planning for the CA/GTA missions can be formulated as a cooperative trajectory optimal control problem (CTOCP).Problem 1 (CTOCP).
Find the trajectories, which drive the system from given initial conditions to desired final conditions over time horizons[
t
0
,
t
f
], while the cooperative objective function is minimized as
(19)
min
⁡
J
=
∑
j
=
1
N
v
J
j
(
x
,
u
)
=
∑
j
=
1
N
v
[
∑
n
=
1
N
NFZ
w
r
n
,
j
(
e
e
-
h
n
-
1
)
w
t
j
(
t
f
j
-
t
0
j
)
+
w
p
j
∫
t
0
j
t
f
j
(
1
-
P
d
j
(
t
,
r
)
)
d
t
h
h
h
h
h
+
w
r
j
∫
t
0
j
t
f
j
∑
n
=
1
N
NFZ
w
r
n
,
j
(
e
e
-
h
n
-
1
)
d
t
]
,
subject to the system equation (18) and the boundary constraints (i.e., the initial and terminal states (5))
(20)
Φ
[
x
j
(
t
0
j
)
,
u
j
(
t
0
j
)
,
t
0
j
,
x
j
(
t
f
j
)
,
u
j
(
t
f
j
)
,
t
f
j
]
=
0
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
∀
j
=
1,2
,
…
,
N
v
,
and several inequality and equality constraints, the individual and cooperative constraints, including the state and control vectors (2) and (9)–(11), are denoted as
(21)
δ
j
[
x
j
(
t
)
,
u
j
(
t
)
,
t
]
≤
0
,
φ
j
[
x
j
(
t
)
,
u
j
(
t
)
,
x
k
(
t
)
,
u
k
(
t
)
,
t
]
≤
0
,
f
T
j
(
t
f
j
,
t
f
k
)
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
.


## 3.1. Objective Function

The objective function of the entire team is a sum of each individual cost, which is a multicriterion objective function, and each criterion could compete with each other as follows:(12)
J
=
min
⁡
⁡
∑
j
=
1
N
v
J
j
=
min
⁡
⁡
∑
j
=
1
N
v
(
w
t
j
J
t
j
+
w
p
j
J
prd
j
+
w
r
j
J
rob
j
)
,
and the criterions can be, respectively, defined as
(13)
J
t
j
=
∫
t
d
t
=
t
f
j
-
t
0
j
,
(14)
J
prd
j
=
∫
t
PR
D
j
(
t
)
=
∫
t
(
1
-
∏
r
=
1
n
r
(
1
-
P
d
j
(
t
,
r
)
)
)
,
(15)
J
rob
j
=
∫
t
R
j
(
t
)
=
∫
t
∑
n
=
1
N
NFZ
w
r
n
,
j
(
e
e
-
h
n
-
1
)
,
where t
0
j and t
f
j denote the initial and terminal time of UCAV
j. The first term, J
t
j, denotes the total fight time of UCAV
j. The second term, J
prd
j, denotes the detection-probability of the n
r-radar system to UCAV
j at the time t, which describes the threat risk criterion of UCAV
j as follows:
(16)
PR
D
j
(
t
)
=
1
-
∏
r
=
1
n
r
(
1
-
P
d
j
(
t
,
r
)
)
,
where P
d
j
(
t
,
r
) is the model of the radar detection probability between the trajectory point of UCAV
j at the time t and the r
th radar [35]. The final term, J
rob
j, guarantees the robustness of the vehicle collision avoidance against the model imprecision, disturbing in the battlefield and operating errors that could otherwise cause vehicles to collide with obstacles or other vehicles. The robustness cost is defined as
(17)
R
j
(
t
)
=
∑
n
=
1
N
NFZ
w
r
n
,
j
(
e
e
-
h
n
-
1
)
,
where w
r
n
,
j are weights that are incremented to increase the maneuver robustness when the trajectory of a vehicle passes in close proximity to an obstacle or another vehicle, which can be chosen by the method presented by Hurni et al. [36].As can be seen, the individual objective functionJ
j is defined by the weighted sum of the three separate running cost terms with appropriate weighting factors w
t
j, w
p
j, and w
r
j. The three criteria (13)–(15) represent different physical meanings, and they are difficult to be compared directly. Hence the selection of these weighting factors is a skilled technique. In the experimentations, the multi-objective fuzzy optimization method proposed by Wang et al. [37] is employed. The method includes two main steps: (a) normalizing each single objective function and (b) solving the membership function about each criterion by using fuzzy distribution function.


## 3.2. Cooperative Trajectory Planning Problem Formulation

After establishing the above models, the cooperative trajectory planning problem is formulated in this section. The problem under consideration is a cooperative scheme, consisting ofN
v homogeneous UCAVs. According to the system equations (1), the general form of the system for UCAV
j is given by
(18)
x
˙
j
(
t
)
=
f
j
[
x
j
(
t
)
,
u
j
(
t
)
,
t
]
,
j
=
1
,
…
,
N
v
.As stated above, to obtain optimal or suboptimal cooperative trajectories, the cooperative trajectory planning for the CA/GTA missions can be formulated as a cooperative trajectory optimal control problem (CTOCP).Problem 1 (CTOCP).
Find the trajectories, which drive the system from given initial conditions to desired final conditions over time horizons[
t
0
,
t
f
], while the cooperative objective function is minimized as
(19)
min
⁡
J
=
∑
j
=
1
N
v
J
j
(
x
,
u
)
=
∑
j
=
1
N
v
[
∑
n
=
1
N
NFZ
w
r
n
,
j
(
e
e
-
h
n
-
1
)
w
t
j
(
t
f
j
-
t
0
j
)
+
w
p
j
∫
t
0
j
t
f
j
(
1
-
P
d
j
(
t
,
r
)
)
d
t
h
h
h
h
h
+
w
r
j
∫
t
0
j
t
f
j
∑
n
=
1
N
NFZ
w
r
n
,
j
(
e
e
-
h
n
-
1
)
d
t
]
,
subject to the system equation (18) and the boundary constraints (i.e., the initial and terminal states (5))
(20)
Φ
[
x
j
(
t
0
j
)
,
u
j
(
t
0
j
)
,
t
0
j
,
x
j
(
t
f
j
)
,
u
j
(
t
f
j
)
,
t
f
j
]
=
0
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
∀
j
=
1,2
,
…
,
N
v
,
and several inequality and equality constraints, the individual and cooperative constraints, including the state and control vectors (2) and (9)–(11), are denoted as
(21)
δ
j
[
x
j
(
t
)
,
u
j
(
t
)
,
t
]
≤
0
,
φ
j
[
x
j
(
t
)
,
u
j
(
t
)
,
x
k
(
t
)
,
u
k
(
t
)
,
t
]
≤
0
,
f
T
j
(
t
f
j
,
t
f
k
)
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
.


## 4. Virtual Motion Camouflage Based Cooperative Trajectory Planning

To tackle the cooperative trajectory planning for CA/GTA missions, an efficient virtual motion camouflage based planning method is proposed for numerically solving the CTOCP, which combines with the benefits of several other classical trajectory generation techniques, including DFT, GPM, and NLP. Correspondingly, this task contains three primary steps. The first is to determine outputs by VMC and DFT in the virtual time domain, such that the dynamics system can be mapped to a lower-dimensional output space. Then, GPM is used to discretize the outputs. Finally, the system is transformed into a NLP problem and the sequential quadratic programming (SQP) is used to solve the NLP to minimize the cost and subject to constraints in output space.


### 4.1. Virtual Motion Camouflage and Problem Formulation in the Output Space

#### 4.1.1. Virtual Motion Camouflage

Inspired by mating hoverflies, Srinivasan and Davey [38] described a new form of the stealth strategy which is the so-called motion camouflage (MC). MC involves two moving objects, a prey and an aggressor, and occurs when the aggressor conceals its motion (i.e., maintains a constant bearing in the prey’s coordinates), except the inevitable size change, as viewed by the moving prey. Here the MC concept is introduced into the optimal trajectory planning and is called virtual motion camouflage (VMC) [21–23] because the “prey” is a virtual one. In the VMC framework, the aggressor trajectory (i.e., the aggressor position vector) ρ
a
(
t
)
∈
ℝ
n
a is confined by the virtual prey trajectory (VPT) ρ
p
(
t
)
∈
ℝ
n
a, the selected reference point ρ
ref
(
t
)
∈
ℝ
n
a, and the path control parameter (PCP) v
(
t
)
∈
ℝ as follows:
(22)
ρ
a
(
t
)
=
ρ
ref
+
v
(
t
)
(
ρ
p
(
t
)
-
ρ
ref
)
.The reference point is considered fixed over time, so the derivatives of the trajectory can be described as(23)
ρ
˙
a
(
t
)
=
v
˙
(
t
)
(
ρ
p
(
t
)
-
ρ
ref
)
+
v
(
t
)
ρ
˙
p
(
t
)
,
ρ
¨
a
(
t
)
=
v
¨
(
t
)
(
ρ
p
(
t
)
-
ρ
ref
)
+
v
(
t
)
ρ
¨
p
(
t
)
+
2
v
˙
(
t
)
ρ
˙
p
(
t
)
.From [23], it can be known that when the proposed VMC suboptimal trajectory design method is used, the complexity of the problem can be reduced by (
n
s
v
+
m
c
v
) times from the collocation method or n
s
v times from the differential inclusion, where n
s
v and m
c
v are the numbers of the state variables and the control variables. In VMC framework, if the dimension of the position vector x
a
(
t
) is one (i.e., n
a
=
1), the dimension of the problem is the same as that of the differential inclusion method. Otherwise, when n
a
>
1, the reference point can be fixed or regarded as a parameter to be optimized, whereas VPT can be predefined or predetermined by the initial and terminal conditions of the aggressor trajectory. Thus, the dimension of the problem is lowered from n
a to one, such that the solving speed can increase.


#### 4.1.2. VMC Problem Formulation in Virtual Time Domain

Obviously, the time factor of trajectories is an argument of the state and control. To deal with temporal constraints and achieve the time cooperation of multi-UCAV, the time along the trajectories should be considered separately. Therefore, the independent intermediate variable (called virtual time here,t
-
∈
[
0,1
]) is introduced and described as
(24)
t
-
=
(
t
-
t
0
)
(
t
f
-
t
0
)
,
such that the trajectories can be generated in the virtual time domain from t
-
0
=
0 to t
-
f
=
1.For UCAVs, which are required to take off at the same time, it can be assumed that the initial time of all UCAVs is zero (t
0
=
0). Thus, the terminal time t
f can be written as t
f
=
t
/
t
-. That is, t
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 to be optimized in the dynamical system, designated as T
f. Then, the following relationship between the virtual time domain and true time domain can be obtained for an arbitrary variable χ as follows:
(25)
χ
′
(
t
-
)
=
d
χ
(
t
)
d
t
-
=
T
f
d
χ
(
t
)
d
t
=
T
f
χ
˙
(
t
)
,
χ
′
′
(
t
-
)
=
d
χ
′
(
t
-
)
d
t
-
=
T
f
d
χ
′
(
t
-
)
d
t
=
T
f
d
(
T
f
χ
˙
(
t
)
)
d
t
=
T
f
2
χ
¨
(
t
)
.Particularly, the derivative of the speed variable can be denoted as(26)
V
(
t
-
)
=
(
x
′
)
2
+
(
y
′
)
2
+
(
z
′
)
2
=
T
f
x
˙
2
+
y
˙
2
+
z
˙
2
=
T
f
V
(
t
)
,
where the superscript “′” represents the derivative with respect to the virtual time.According to (18) and (24), the system equations can be rearranged as
(27)
x
′
j
(
t
-
)
=
T
f
x
˙
j
(
t
)
=
T
f
f
j
[
x
j
(
t
)
,
u
j
(
t
)
,
t
]
=
g
j
[
x
j
(
t
-
)
,
u
j
(
t
-
)
,
T
f
j
]
,
T
f
′
j
=
0
,
j
=
1
,
…
,
N
v
.Hence, the Problem1 can be reformulated in virtual time domain (VTD) as follows.Problem 2 (CTOCP-VTD).
Minimize the cooperative objective function (19) of all UCAVs represented with respect to the new independent variable t
- as
(28)
min
⁡
J
=
∑
j
=
1
N
v
J
j
(
x
,
u
,
T
f
j
)
=
∑
j
=
1
N
v
[
w
t
j
T
f
j
+
w
p
j
∫
0
1
PR
D
j
(
t
-
)
d
t
-
h
h
h
h
+
w
r
j
∫
0
1
R
j
(
t
-
)
d
t
-
]
,
subject to the system equations (27), and the boundary constraints (20), written as
(29)
ϕ
[
x
j
(
0
)
,
u
j
(
0
)
,
x
j
(
1
)
,
u
j
(
1
)
,
T
f
j
]
=
0
,
Meanwhile satisfying the inequality and equality constraints (21) and additional temporal constraints as follows:
(30)
δ
j
[
x
j
(
t
-
)
,
u
j
(
t
-
)
,
T
f
j
]
≤
0
,
φ
j
[
x
j
(
t
-
)
,
u
j
(
t
-
)
,
x
k
(
t
-
)
,
u
k
(
t
-
)
,
T
f
j
,
T
f
k
]
≤
0
,
f
T
j
(
t
f
j
,
t
f
k
)
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
.From (1), it can be clearly known that the dynamics system under consideration is a higher-dimensional space system. Due to the complexity in solving this type of the problem, the VMC based method is used to lower the system dimension and the state and control vectors can be related by the differential flatness theory based method [39]. In this system, the state vector x
∈
ℝ
6 can be separated into two parts: the trajectory position (e.g., spatial trajectory) ρ
a
=
[
x
,
y
,
z
]
∈
ℝ
3 and the corresponding state rate ρ
r
=
[
V
,
ψ
,
γ
]
∈
ℝ
3, by which the trajectory position in virtual time domain and the terminal time T
f can be defined as the flat output vector
(31)
ϑ
=
[
ρ
a
(
t
-
)
,
T
f
]
=
[
x
(
t
-
)
,
y
(
t
-
)
,
z
(
t
-
)
,
T
f
]
T
.Thus, the original state vectorx and control vector u can be recovered from the flat outputs and their derivatives as
(32)
x
=
α
(
ϑ
,
ϑ
′
,
ϑ
′
′
)
=
α
(
ρ
a
,
ρ
a
′
,
ρ
a
′
′
,
T
f
)
,
u
=
β
(
ϑ
,
ϑ
′
,
ϑ
′
′
)
=
β
(
ρ
a
,
ρ
a
′
,
ρ
a
′
′
,
T
f
)
.According to (1), the relation of the state vector and control vector in the VTD is described as follows:
(33)
V
=
(
x
′
)
2
+
(
y
′
)
2
+
(
z
′
)
2
,
ψ
=
arctan
⁡
(
y
′
x
′
)
,
γ
=
arctan
⁡
(
z
′
(
x
′
)
2
+
(
y
′
)
2
)
,
n
x
=
V
′
(
T
f
2
g
)
+
sin
⁡
γ
,
n
z
=
(
V
γ
′
+
T
f
2
g
cos
⁡
⁡
γ
)
2
+
(
V
ψ
′
cos
⁡
⁡
γ
)
2
T
f
g
,
μ
=
arctan
⁡
(
V
ψ
′
cos
⁡
γ
(
V
γ
′
+
T
f
2
g
cos
⁡
γ
)
)
.In the VMC framework, the trajectory positionρ
a can be defined as the aggressor position vector. As described above, the state and the control variables of the dynamics system can be recovered from the flat outputs. Hence, according to (22), these variables are also functions of the PCP, the predefined virtual prey motion, the reference point, and their corresponding derivatives. Once the predefined virtual prey motion and the reference point are selected, the system of cooperative trajectory planning of UCAVs, including the objective function and the constraints, is mapped to a lower-dimensional output space (here, the dimension is one). Obviously, the dynamic constraints of this system (27) can be automatically satisfied. Therefore, the Problem 2 can be redefined as follows.Problem 3 (CTOCP-VTD in the output space).
Consider(34)
min
⁡
J
(
v
)
=
∑
j
=
1
N
v
J
j
(
v
j
)
=
∑
j
=
1
N
v
[
∫
(
v
j
,
v
′
j
,
…
,
T
f
j
)
h
h
h
h
h
h
h
+
∫
0
1
f
j
(
v
j
,
v
′
j
,
…
,
t
-
)
d
t
-
]
s
.
t
.
i
ϕ
j
[
v
j
(
0
)
,
v
j
(
1
)
,
T
f
j
]
=
0
,
h
h
h
h
h
δ
j
(
v
j
,
v
′
j
,
…
,
)
≤
0
,
h
h
h
h
h
φ
j
(
v
j
,
v
′
j
,
…
,
v
k
,
v
′
k
,
…
)
≤
0
,
h
h
h
h
f
T
j
(
t
f
j
,
t
f
k
)
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
.


### 4.2. Collocation Using Gauss Pseudospectral Method

In order to solve the Problem3 through a NLP algorithm, the PCP history v
(
t
-
) should be discretized into c
=
0,1
,
…
,
N
c
,
N
c
+
1 nodes, with v
0
=
v
(
t
-
0
) and v
N
c
+
1
=
v
(
t
-
N
c
+
1
). In this paper, GPM is selected as the discretization method [12], which is an orthogonal collection method where the collocation points are the Legendre-Gauss (LG) points. Similar to other pseudospectral methods, a finite basis of global interpolating polynomials is used to approximate the state and control space at a set of discretization nodes in the GPM and then the optimal control equations can be transformed into nonlinear algebra equations, such that the OCP can be solved by a NLP algorithm.The standard interval considered here is denoted asτ
∈
[
-
1,1
]. By using a linear transformation, the virtual time t
- can be expressed as a function of τ via
(35)
t
-
=
[
(
t
-
f
-
t
-
0
)
τ
+
(
t
-
f
+
t
-
0
)
]
2
=
(
τ
+
1
)
2
.Letτ
1
<
⋯
<
τ
N
c be collocation points in (
-
1,1
) and τ
0
=
-
1, τ
f
=
1. The time history of the approximated PCP at the N
c LG points is given by
(36)
v
(
τ
)
≈
v
-
(
τ
)
=
∑
c
=
0
N
c
v
(
τ
c
)
L
c
(
τ
)
,
where v
(
τ
c
) is an approximant at the node τ
c and the base function L
c
(
τ
) is the Lagrange interpolating polynomial of degree N
c, expressed as
(37)
L
c
(
τ
)
=
∏
l
=
0
,
l
≠
c
N
c
τ
-
τ
l
τ
c
-
τ
l
,
c
=
0,1
,
…
,
N
c
,
satisfying the isolation property
(38)
L
c
(
τ
l
)
=
δ
c
l
=
{
1
,
c
=
l
,
0
,
c
≠
l
.And the endpoint can be approximated according to the following formula:(39)
v
(
τ
f
)
≈
v
-
(
τ
f
)
=
v
-
(
τ
0
)
+
t
-
f
-
t
-
0
2
×
∑
c
=
1
N
c
w
c
f
[
v
-
(
τ
c
)
,
τ
c
]
,
where τ
c
(
c
=
1
,
…
,
N
c
) are the LG points and the variables ω
c are the LG weights given by
(40)
w
c
=
∫
-
1
1
L
c
(
τ
)
d
τ
.The derivative of the state can be approximated as the exact derivative of the interpolating polynomial. Evaluating the derivative at the LG points results in(41)
v
′
(
τ
l
)
≈
v
-
′
(
τ
l
)
=
∑
c
=
0
N
c
v
(
τ
c
)
L
c
′
(
τ
l
)
=
∑
c
=
0
N
c
v
(
τ
c
)
D
l
c
,
l
=
1
,
…
,
N
c
,
where D
=
[
D
l
c
]
∈
ℝ
N
c
×
(
N
c
+
1
) is a differential matrix and can be offline determined as
(42)
D
l
c
=
∑
k
1
=
0
N
c
∏
k
2
=
0
,
k
2
≠
c
,
k
1
N
c
(
τ
l
-
τ
k
2
)
∏
k
2
=
0
,
k
2
≠
c
N
c
(
τ
c
-
τ
k
2
)
,
h
h
h
h
h
h
h
l
=
1
,
…
N
c
,
c
=
0
,
…
,
N
c
.The discretized states and controls satisfying the vehicle dynamics can be computed by ensuring that equation(43)
∑
c
=
0
N
c
v
(
τ
c
)
D
l
c
-
t
-
f
-
t
-
0
2
f
[
v
(
τ
l
)
,
τ
l
]
=
0
,
l
=
1
,
…
,
N
c
.As described above, the continuous OCP can be transformed to the discretized NLP problem. The parameters to be optimized are the PCP nodesv
=
[
v
c
]
c
=
0,1
,
…
N
c
+
1. When the prey motion ρ
p
,
c is equal to the aggressor motion ρ
a
,
c for c
=
0
,
N
c
+
1, one can assume that v
0
=
v
N
c
=
1, so the boundary conditions are no longer considered. And the discretized NLP problem can be denoted as follows.Problem 4 (CTP-NLP).
Its standard form is denoted as(44)
min
⁡
J
(
v
)
=
∑
j
=
1
N
v
J
j
(
v
j
)
=
∑
j
=
1
N
v
[
∑
j
=
1
N
v
Φ
j
[
v
j
,
T
f
j
]
h
h
h
h
h
h
h
h
+
t
-
f
-
t
-
0
2
∑
c
=
1
N
c
f
j
(
v
j
,
T
f
j
)
ω
c
]
s
.
t
.
i
δ
l
j
[
v
j
]
≤
0
,
l
=
1
,
…
,
N
c
,
h
h
h
h
h
φ
j
[
v
j
,
v
k
]
≤
0
,
h
h
h
h
h
f
T
j
(
t
f
j
,
t
f
k
)
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
,
where v
j
=
[
v
c
j
]
c
=
1
,
…
N
c is the vector form of the discretized PCP.Then the resulting Problem4 can be solved through well-developed algorithms, such as the SQP algorithm. In this paper, TOMLAB/SNOPT software toolbox is chosen due to its advantages on solution effectiveness for the large-scale NLP problems [40].


### 4.3. Convergence and Initial Guess

As analyzed in [41], the GPM exhibits global convergence properties in many applications. In addition, the convergence of the virtual motion camouflage has been proved by Xu et al. [21, 23], in which [23] showed that if the problem was only one-dimensional, the VMC based trajectory design method could converge to the optimal solution, but when the problem was multi-dimensional and highly nonlinear, the result might be suboptimal in the full solution space, and [21] showed that the VMC framework together with the pseudospectral discretization method could converge to the optimal solution.By Introducing elastic programming concept and some advances in both mathematical programming techniques and pseudospectral methods, Ross and Gong [42] designed a guess-free trajectory optimization algorithm. These advances have the effect of eliminating the guessing problem in the trajectory optimization, which has been used in the current GPM algorithm. In the proposed VMC method, the motion of the virtual prey and a reference point are guessed. However, both of them have direct physical relation to the problem. Typically the virtual prey motion is selected according to the boundary conditions and a proper guess is not a difficult task [23]. In this paper, the trivial guess for the reference point can be the initial point (IP), and an initial guess of v
c
=
1 is chosen for c
=
1
,
…
,
N
c. The virtual prey motion is defined as a 2-order curve, determined by the initial and terminal conditions.


## 4.1. Virtual Motion Camouflage and Problem Formulation in the Output Space

### 4.1.1. Virtual Motion Camouflage

Inspired by mating hoverflies, Srinivasan and Davey [38] described a new form of the stealth strategy which is the so-called motion camouflage (MC). MC involves two moving objects, a prey and an aggressor, and occurs when the aggressor conceals its motion (i.e., maintains a constant bearing in the prey’s coordinates), except the inevitable size change, as viewed by the moving prey. Here the MC concept is introduced into the optimal trajectory planning and is called virtual motion camouflage (VMC) [21–23] because the “prey” is a virtual one. In the VMC framework, the aggressor trajectory (i.e., the aggressor position vector) ρ
a
(
t
)
∈
ℝ
n
a is confined by the virtual prey trajectory (VPT) ρ
p
(
t
)
∈
ℝ
n
a, the selected reference point ρ
ref
(
t
)
∈
ℝ
n
a, and the path control parameter (PCP) v
(
t
)
∈
ℝ as follows:
(22)
ρ
a
(
t
)
=
ρ
ref
+
v
(
t
)
(
ρ
p
(
t
)
-
ρ
ref
)
.The reference point is considered fixed over time, so the derivatives of the trajectory can be described as(23)
ρ
˙
a
(
t
)
=
v
˙
(
t
)
(
ρ
p
(
t
)
-
ρ
ref
)
+
v
(
t
)
ρ
˙
p
(
t
)
,
ρ
¨
a
(
t
)
=
v
¨
(
t
)
(
ρ
p
(
t
)
-
ρ
ref
)
+
v
(
t
)
ρ
¨
p
(
t
)
+
2
v
˙
(
t
)
ρ
˙
p
(
t
)
.From [23], it can be known that when the proposed VMC suboptimal trajectory design method is used, the complexity of the problem can be reduced by (
n
s
v
+
m
c
v
) times from the collocation method or n
s
v times from the differential inclusion, where n
s
v and m
c
v are the numbers of the state variables and the control variables. In VMC framework, if the dimension of the position vector x
a
(
t
) is one (i.e., n
a
=
1), the dimension of the problem is the same as that of the differential inclusion method. Otherwise, when n
a
>
1, the reference point can be fixed or regarded as a parameter to be optimized, whereas VPT can be predefined or predetermined by the initial and terminal conditions of the aggressor trajectory. Thus, the dimension of the problem is lowered from n
a to one, such that the solving speed can increase.


### 4.1.2. VMC Problem Formulation in Virtual Time Domain

Obviously, the time factor of trajectories is an argument of the state and control. To deal with temporal constraints and achieve the time cooperation of multi-UCAV, the time along the trajectories should be considered separately. Therefore, the independent intermediate variable (called virtual time here,t
-
∈
[
0,1
]) is introduced and described as
(24)
t
-
=
(
t
-
t
0
)
(
t
f
-
t
0
)
,
such that the trajectories can be generated in the virtual time domain from t
-
0
=
0 to t
-
f
=
1.For UCAVs, which are required to take off at the same time, it can be assumed that the initial time of all UCAVs is zero (t
0
=
0). Thus, the terminal time t
f can be written as t
f
=
t
/
t
-. That is, t
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 to be optimized in the dynamical system, designated as T
f. Then, the following relationship between the virtual time domain and true time domain can be obtained for an arbitrary variable χ as follows:
(25)
χ
′
(
t
-
)
=
d
χ
(
t
)
d
t
-
=
T
f
d
χ
(
t
)
d
t
=
T
f
χ
˙
(
t
)
,
χ
′
′
(
t
-
)
=
d
χ
′
(
t
-
)
d
t
-
=
T
f
d
χ
′
(
t
-
)
d
t
=
T
f
d
(
T
f
χ
˙
(
t
)
)
d
t
=
T
f
2
χ
¨
(
t
)
.Particularly, the derivative of the speed variable can be denoted as(26)
V
(
t
-
)
=
(
x
′
)
2
+
(
y
′
)
2
+
(
z
′
)
2
=
T
f
x
˙
2
+
y
˙
2
+
z
˙
2
=
T
f
V
(
t
)
,
where the superscript “′” represents the derivative with respect to the virtual time.According to (18) and (24), the system equations can be rearranged as
(27)
x
′
j
(
t
-
)
=
T
f
x
˙
j
(
t
)
=
T
f
f
j
[
x
j
(
t
)
,
u
j
(
t
)
,
t
]
=
g
j
[
x
j
(
t
-
)
,
u
j
(
t
-
)
,
T
f
j
]
,
T
f
′
j
=
0
,
j
=
1
,
…
,
N
v
.Hence, the Problem1 can be reformulated in virtual time domain (VTD) as follows.Problem 2 (CTOCP-VTD).
Minimize the cooperative objective function (19) of all UCAVs represented with respect to the new independent variable t
- as
(28)
min
⁡
J
=
∑
j
=
1
N
v
J
j
(
x
,
u
,
T
f
j
)
=
∑
j
=
1
N
v
[
w
t
j
T
f
j
+
w
p
j
∫
0
1
PR
D
j
(
t
-
)
d
t
-
h
h
h
h
+
w
r
j
∫
0
1
R
j
(
t
-
)
d
t
-
]
,
subject to the system equations (27), and the boundary constraints (20), written as
(29)
ϕ
[
x
j
(
0
)
,
u
j
(
0
)
,
x
j
(
1
)
,
u
j
(
1
)
,
T
f
j
]
=
0
,
Meanwhile satisfying the inequality and equality constraints (21) and additional temporal constraints as follows:
(30)
δ
j
[
x
j
(
t
-
)
,
u
j
(
t
-
)
,
T
f
j
]
≤
0
,
φ
j
[
x
j
(
t
-
)
,
u
j
(
t
-
)
,
x
k
(
t
-
)
,
u
k
(
t
-
)
,
T
f
j
,
T
f
k
]
≤
0
,
f
T
j
(
t
f
j
,
t
f
k
)
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
.From (1), it can be clearly known that the dynamics system under consideration is a higher-dimensional space system. Due to the complexity in solving this type of the problem, the VMC based method is used to lower the system dimension and the state and control vectors can be related by the differential flatness theory based method [39]. In this system, the state vector x
∈
ℝ
6 can be separated into two parts: the trajectory position (e.g., spatial trajectory) ρ
a
=
[
x
,
y
,
z
]
∈
ℝ
3 and the corresponding state rate ρ
r
=
[
V
,
ψ
,
γ
]
∈
ℝ
3, by which the trajectory position in virtual time domain and the terminal time T
f can be defined as the flat output vector
(31)
ϑ
=
[
ρ
a
(
t
-
)
,
T
f
]
=
[
x
(
t
-
)
,
y
(
t
-
)
,
z
(
t
-
)
,
T
f
]
T
.Thus, the original state vectorx and control vector u can be recovered from the flat outputs and their derivatives as
(32)
x
=
α
(
ϑ
,
ϑ
′
,
ϑ
′
′
)
=
α
(
ρ
a
,
ρ
a
′
,
ρ
a
′
′
,
T
f
)
,
u
=
β
(
ϑ
,
ϑ
′
,
ϑ
′
′
)
=
β
(
ρ
a
,
ρ
a
′
,
ρ
a
′
′
,
T
f
)
.According to (1), the relation of the state vector and control vector in the VTD is described as follows:
(33)
V
=
(
x
′
)
2
+
(
y
′
)
2
+
(
z
′
)
2
,
ψ
=
arctan
⁡
(
y
′
x
′
)
,
γ
=
arctan
⁡
(
z
′
(
x
′
)
2
+
(
y
′
)
2
)
,
n
x
=
V
′
(
T
f
2
g
)
+
sin
⁡
γ
,
n
z
=
(
V
γ
′
+
T
f
2
g
cos
⁡
⁡
γ
)
2
+
(
V
ψ
′
cos
⁡
⁡
γ
)
2
T
f
g
,
μ
=
arctan
⁡
(
V
ψ
′
cos
⁡
γ
(
V
γ
′
+
T
f
2
g
cos
⁡
γ
)
)
.In the VMC framework, the trajectory positionρ
a can be defined as the aggressor position vector. As described above, the state and the control variables of the dynamics system can be recovered from the flat outputs. Hence, according to (22), these variables are also functions of the PCP, the predefined virtual prey motion, the reference point, and their corresponding derivatives. Once the predefined virtual prey motion and the reference point are selected, the system of cooperative trajectory planning of UCAVs, including the objective function and the constraints, is mapped to a lower-dimensional output space (here, the dimension is one). Obviously, the dynamic constraints of this system (27) can be automatically satisfied. Therefore, the Problem 2 can be redefined as follows.Problem 3 (CTOCP-VTD in the output space).
Consider(34)
min
⁡
J
(
v
)
=
∑
j
=
1
N
v
J
j
(
v
j
)
=
∑
j
=
1
N
v
[
∫
(
v
j
,
v
′
j
,
…
,
T
f
j
)
h
h
h
h
h
h
h
+
∫
0
1
f
j
(
v
j
,
v
′
j
,
…
,
t
-
)
d
t
-
]
s
.
t
.
i
ϕ
j
[
v
j
(
0
)
,
v
j
(
1
)
,
T
f
j
]
=
0
,
h
h
h
h
h
δ
j
(
v
j
,
v
′
j
,
…
,
)
≤
0
,
h
h
h
h
h
φ
j
(
v
j
,
v
′
j
,
…
,
v
k
,
v
′
k
,
…
)
≤
0
,
h
h
h
h
f
T
j
(
t
f
j
,
t
f
k
)
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
.


## 4.1.1. Virtual Motion Camouflage

Inspired by mating hoverflies, Srinivasan and Davey [38] described a new form of the stealth strategy which is the so-called motion camouflage (MC). MC involves two moving objects, a prey and an aggressor, and occurs when the aggressor conceals its motion (i.e., maintains a constant bearing in the prey’s coordinates), except the inevitable size change, as viewed by the moving prey. Here the MC concept is introduced into the optimal trajectory planning and is called virtual motion camouflage (VMC) [21–23] because the “prey” is a virtual one. In the VMC framework, the aggressor trajectory (i.e., the aggressor position vector) ρ
a
(
t
)
∈
ℝ
n
a is confined by the virtual prey trajectory (VPT) ρ
p
(
t
)
∈
ℝ
n
a, the selected reference point ρ
ref
(
t
)
∈
ℝ
n
a, and the path control parameter (PCP) v
(
t
)
∈
ℝ as follows:
(22)
ρ
a
(
t
)
=
ρ
ref
+
v
(
t
)
(
ρ
p
(
t
)
-
ρ
ref
)
.The reference point is considered fixed over time, so the derivatives of the trajectory can be described as(23)
ρ
˙
a
(
t
)
=
v
˙
(
t
)
(
ρ
p
(
t
)
-
ρ
ref
)
+
v
(
t
)
ρ
˙
p
(
t
)
,
ρ
¨
a
(
t
)
=
v
¨
(
t
)
(
ρ
p
(
t
)
-
ρ
ref
)
+
v
(
t
)
ρ
¨
p
(
t
)
+
2
v
˙
(
t
)
ρ
˙
p
(
t
)
.From [23], it can be known that when the proposed VMC suboptimal trajectory design method is used, the complexity of the problem can be reduced by (
n
s
v
+
m
c
v
) times from the collocation method or n
s
v times from the differential inclusion, where n
s
v and m
c
v are the numbers of the state variables and the control variables. In VMC framework, if the dimension of the position vector x
a
(
t
) is one (i.e., n
a
=
1), the dimension of the problem is the same as that of the differential inclusion method. Otherwise, when n
a
>
1, the reference point can be fixed or regarded as a parameter to be optimized, whereas VPT can be predefined or predetermined by the initial and terminal conditions of the aggressor trajectory. Thus, the dimension of the problem is lowered from n
a to one, such that the solving speed can increase.


## 4.1.2. VMC Problem Formulation in Virtual Time Domain

Obviously, the time factor of trajectories is an argument of the state and control. To deal with temporal constraints and achieve the time cooperation of multi-UCAV, the time along the trajectories should be considered separately. Therefore, the independent intermediate variable (called virtual time here,t
-
∈
[
0,1
]) is introduced and described as
(24)
t
-
=
(
t
-
t
0
)
(
t
f
-
t
0
)
,
such that the trajectories can be generated in the virtual time domain from t
-
0
=
0 to t
-
f
=
1.For UCAVs, which are required to take off at the same time, it can be assumed that the initial time of all UCAVs is zero (t
0
=
0). Thus, the terminal time t
f can be written as t
f
=
t
/
t
-. That is, t
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 to be optimized in the dynamical system, designated as T
f. Then, the following relationship between the virtual time domain and true time domain can be obtained for an arbitrary variable χ as follows:
(25)
χ
′
(
t
-
)
=
d
χ
(
t
)
d
t
-
=
T
f
d
χ
(
t
)
d
t
=
T
f
χ
˙
(
t
)
,
χ
′
′
(
t
-
)
=
d
χ
′
(
t
-
)
d
t
-
=
T
f
d
χ
′
(
t
-
)
d
t
=
T
f
d
(
T
f
χ
˙
(
t
)
)
d
t
=
T
f
2
χ
¨
(
t
)
.Particularly, the derivative of the speed variable can be denoted as(26)
V
(
t
-
)
=
(
x
′
)
2
+
(
y
′
)
2
+
(
z
′
)
2
=
T
f
x
˙
2
+
y
˙
2
+
z
˙
2
=
T
f
V
(
t
)
,
where the superscript “′” represents the derivative with respect to the virtual time.According to (18) and (24), the system equations can be rearranged as
(27)
x
′
j
(
t
-
)
=
T
f
x
˙
j
(
t
)
=
T
f
f
j
[
x
j
(
t
)
,
u
j
(
t
)
,
t
]
=
g
j
[
x
j
(
t
-
)
,
u
j
(
t
-
)
,
T
f
j
]
,
T
f
′
j
=
0
,
j
=
1
,
…
,
N
v
.Hence, the Problem1 can be reformulated in virtual time domain (VTD) as follows.Problem 2 (CTOCP-VTD).
Minimize the cooperative objective function (19) of all UCAVs represented with respect to the new independent variable t
- as
(28)
min
⁡
J
=
∑
j
=
1
N
v
J
j
(
x
,
u
,
T
f
j
)
=
∑
j
=
1
N
v
[
w
t
j
T
f
j
+
w
p
j
∫
0
1
PR
D
j
(
t
-
)
d
t
-
h
h
h
h
+
w
r
j
∫
0
1
R
j
(
t
-
)
d
t
-
]
,
subject to the system equations (27), and the boundary constraints (20), written as
(29)
ϕ
[
x
j
(
0
)
,
u
j
(
0
)
,
x
j
(
1
)
,
u
j
(
1
)
,
T
f
j
]
=
0
,
Meanwhile satisfying the inequality and equality constraints (21) and additional temporal constraints as follows:
(30)
δ
j
[
x
j
(
t
-
)
,
u
j
(
t
-
)
,
T
f
j
]
≤
0
,
φ
j
[
x
j
(
t
-
)
,
u
j
(
t
-
)
,
x
k
(
t
-
)
,
u
k
(
t
-
)
,
T
f
j
,
T
f
k
]
≤
0
,
f
T
j
(
t
f
j
,
t
f
k
)
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
.From (1), it can be clearly known that the dynamics system under consideration is a higher-dimensional space system. Due to the complexity in solving this type of the problem, the VMC based method is used to lower the system dimension and the state and control vectors can be related by the differential flatness theory based method [39]. In this system, the state vector x
∈
ℝ
6 can be separated into two parts: the trajectory position (e.g., spatial trajectory) ρ
a
=
[
x
,
y
,
z
]
∈
ℝ
3 and the corresponding state rate ρ
r
=
[
V
,
ψ
,
γ
]
∈
ℝ
3, by which the trajectory position in virtual time domain and the terminal time T
f can be defined as the flat output vector
(31)
ϑ
=
[
ρ
a
(
t
-
)
,
T
f
]
=
[
x
(
t
-
)
,
y
(
t
-
)
,
z
(
t
-
)
,
T
f
]
T
.Thus, the original state vectorx and control vector u can be recovered from the flat outputs and their derivatives as
(32)
x
=
α
(
ϑ
,
ϑ
′
,
ϑ
′
′
)
=
α
(
ρ
a
,
ρ
a
′
,
ρ
a
′
′
,
T
f
)
,
u
=
β
(
ϑ
,
ϑ
′
,
ϑ
′
′
)
=
β
(
ρ
a
,
ρ
a
′
,
ρ
a
′
′
,
T
f
)
.According to (1), the relation of the state vector and control vector in the VTD is described as follows:
(33)
V
=
(
x
′
)
2
+
(
y
′
)
2
+
(
z
′
)
2
,
ψ
=
arctan
⁡
(
y
′
x
′
)
,
γ
=
arctan
⁡
(
z
′
(
x
′
)
2
+
(
y
′
)
2
)
,
n
x
=
V
′
(
T
f
2
g
)
+
sin
⁡
γ
,
n
z
=
(
V
γ
′
+
T
f
2
g
cos
⁡
⁡
γ
)
2
+
(
V
ψ
′
cos
⁡
⁡
γ
)
2
T
f
g
,
μ
=
arctan
⁡
(
V
ψ
′
cos
⁡
γ
(
V
γ
′
+
T
f
2
g
cos
⁡
γ
)
)
.In the VMC framework, the trajectory positionρ
a can be defined as the aggressor position vector. As described above, the state and the control variables of the dynamics system can be recovered from the flat outputs. Hence, according to (22), these variables are also functions of the PCP, the predefined virtual prey motion, the reference point, and their corresponding derivatives. Once the predefined virtual prey motion and the reference point are selected, the system of cooperative trajectory planning of UCAVs, including the objective function and the constraints, is mapped to a lower-dimensional output space (here, the dimension is one). Obviously, the dynamic constraints of this system (27) can be automatically satisfied. Therefore, the Problem 2 can be redefined as follows.Problem 3 (CTOCP-VTD in the output space).
Consider(34)
min
⁡
J
(
v
)
=
∑
j
=
1
N
v
J
j
(
v
j
)
=
∑
j
=
1
N
v
[
∫
(
v
j
,
v
′
j
,
…
,
T
f
j
)
h
h
h
h
h
h
h
+
∫
0
1
f
j
(
v
j
,
v
′
j
,
…
,
t
-
)
d
t
-
]
s
.
t
.
i
ϕ
j
[
v
j
(
0
)
,
v
j
(
1
)
,
T
f
j
]
=
0
,
h
h
h
h
h
δ
j
(
v
j
,
v
′
j
,
…
,
)
≤
0
,
h
h
h
h
h
φ
j
(
v
j
,
v
′
j
,
…
,
v
k
,
v
′
k
,
…
)
≤
0
,
h
h
h
h
f
T
j
(
t
f
j
,
t
f
k
)
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
.


## 4.2. Collocation Using Gauss Pseudospectral Method

In order to solve the Problem3 through a NLP algorithm, the PCP history v
(
t
-
) should be discretized into c
=
0,1
,
…
,
N
c
,
N
c
+
1 nodes, with v
0
=
v
(
t
-
0
) and v
N
c
+
1
=
v
(
t
-
N
c
+
1
). In this paper, GPM is selected as the discretization method [12], which is an orthogonal collection method where the collocation points are the Legendre-Gauss (LG) points. Similar to other pseudospectral methods, a finite basis of global interpolating polynomials is used to approximate the state and control space at a set of discretization nodes in the GPM and then the optimal control equations can be transformed into nonlinear algebra equations, such that the OCP can be solved by a NLP algorithm.The standard interval considered here is denoted asτ
∈
[
-
1,1
]. By using a linear transformation, the virtual time t
- can be expressed as a function of τ via
(35)
t
-
=
[
(
t
-
f
-
t
-
0
)
τ
+
(
t
-
f
+
t
-
0
)
]
2
=
(
τ
+
1
)
2
.Letτ
1
<
⋯
<
τ
N
c be collocation points in (
-
1,1
) and τ
0
=
-
1, τ
f
=
1. The time history of the approximated PCP at the N
c LG points is given by
(36)
v
(
τ
)
≈
v
-
(
τ
)
=
∑
c
=
0
N
c
v
(
τ
c
)
L
c
(
τ
)
,
where v
(
τ
c
) is an approximant at the node τ
c and the base function L
c
(
τ
) is the Lagrange interpolating polynomial of degree N
c, expressed as
(37)
L
c
(
τ
)
=
∏
l
=
0
,
l
≠
c
N
c
τ
-
τ
l
τ
c
-
τ
l
,
c
=
0,1
,
…
,
N
c
,
satisfying the isolation property
(38)
L
c
(
τ
l
)
=
δ
c
l
=
{
1
,
c
=
l
,
0
,
c
≠
l
.And the endpoint can be approximated according to the following formula:(39)
v
(
τ
f
)
≈
v
-
(
τ
f
)
=
v
-
(
τ
0
)
+
t
-
f
-
t
-
0
2
×
∑
c
=
1
N
c
w
c
f
[
v
-
(
τ
c
)
,
τ
c
]
,
where τ
c
(
c
=
1
,
…
,
N
c
) are the LG points and the variables ω
c are the LG weights given by
(40)
w
c
=
∫
-
1
1
L
c
(
τ
)
d
τ
.The derivative of the state can be approximated as the exact derivative of the interpolating polynomial. Evaluating the derivative at the LG points results in(41)
v
′
(
τ
l
)
≈
v
-
′
(
τ
l
)
=
∑
c
=
0
N
c
v
(
τ
c
)
L
c
′
(
τ
l
)
=
∑
c
=
0
N
c
v
(
τ
c
)
D
l
c
,
l
=
1
,
…
,
N
c
,
where D
=
[
D
l
c
]
∈
ℝ
N
c
×
(
N
c
+
1
) is a differential matrix and can be offline determined as
(42)
D
l
c
=
∑
k
1
=
0
N
c
∏
k
2
=
0
,
k
2
≠
c
,
k
1
N
c
(
τ
l
-
τ
k
2
)
∏
k
2
=
0
,
k
2
≠
c
N
c
(
τ
c
-
τ
k
2
)
,
h
h
h
h
h
h
h
l
=
1
,
…
N
c
,
c
=
0
,
…
,
N
c
.The discretized states and controls satisfying the vehicle dynamics can be computed by ensuring that equation(43)
∑
c
=
0
N
c
v
(
τ
c
)
D
l
c
-
t
-
f
-
t
-
0
2
f
[
v
(
τ
l
)
,
τ
l
]
=
0
,
l
=
1
,
…
,
N
c
.As described above, the continuous OCP can be transformed to the discretized NLP problem. The parameters to be optimized are the PCP nodesv
=
[
v
c
]
c
=
0,1
,
…
N
c
+
1. When the prey motion ρ
p
,
c is equal to the aggressor motion ρ
a
,
c for c
=
0
,
N
c
+
1, one can assume that v
0
=
v
N
c
=
1, so the boundary conditions are no longer considered. And the discretized NLP problem can be denoted as follows.Problem 4 (CTP-NLP).
Its standard form is denoted as(44)
min
⁡
J
(
v
)
=
∑
j
=
1
N
v
J
j
(
v
j
)
=
∑
j
=
1
N
v
[
∑
j
=
1
N
v
Φ
j
[
v
j
,
T
f
j
]
h
h
h
h
h
h
h
h
+
t
-
f
-
t
-
0
2
∑
c
=
1
N
c
f
j
(
v
j
,
T
f
j
)
ω
c
]
s
.
t
.
i
δ
l
j
[
v
j
]
≤
0
,
l
=
1
,
…
,
N
c
,
h
h
h
h
h
φ
j
[
v
j
,
v
k
]
≤
0
,
h
h
h
h
h
f
T
j
(
t
f
j
,
t
f
k
)
≤
0
,
∀
j
≠
k
,
j
,
k
=
1,2
,
…
,
N
v
,
where v
j
=
[
v
c
j
]
c
=
1
,
…
N
c is the vector form of the discretized PCP.Then the resulting Problem4 can be solved through well-developed algorithms, such as the SQP algorithm. In this paper, TOMLAB/SNOPT software toolbox is chosen due to its advantages on solution effectiveness for the large-scale NLP problems [40].


## 4.3. Convergence and Initial Guess

As analyzed in [41], the GPM exhibits global convergence properties in many applications. In addition, the convergence of the virtual motion camouflage has been proved by Xu et al. [21, 23], in which [23] showed that if the problem was only one-dimensional, the VMC based trajectory design method could converge to the optimal solution, but when the problem was multi-dimensional and highly nonlinear, the result might be suboptimal in the full solution space, and [21] showed that the VMC framework together with the pseudospectral discretization method could converge to the optimal solution.By Introducing elastic programming concept and some advances in both mathematical programming techniques and pseudospectral methods, Ross and Gong [42] designed a guess-free trajectory optimization algorithm. These advances have the effect of eliminating the guessing problem in the trajectory optimization, which has been used in the current GPM algorithm. In the proposed VMC method, the motion of the virtual prey and a reference point are guessed. However, both of them have direct physical relation to the problem. Typically the virtual prey motion is selected according to the boundary conditions and a proper guess is not a difficult task [23]. In this paper, the trivial guess for the reference point can be the initial point (IP), and an initial guess of v
c
=
1 is chosen for c
=
1
,
…
,
N
c. The virtual prey motion is defined as a 2-order curve, determined by the initial and terminal conditions.


## 5. Numerical Examples

The basic ideas presented in this paper are illustrated in the following three examples. The specific vehicle platform used in simulations is the Storm shadow UCAV. Its relevant parameters are all taken from [43], summarized in Table 1.Table 1
State and control constraints of UCAVs.
Item
Minimum value
Maximum value
Flight altitude (m)
z
min
⁡
=
200
z
max
⁡
=
8
000
Airspeed (m/s)
V
min
⁡
=
60
V
max
⁡
=
300
Flight-path angle (deg)
γ
min
⁡
=
-
89
γ
max
⁡
=
89
Roll angle (deg)
μ
min
⁡
=
-
80
μ
max
⁡
=
80
Tangential load factor (g)
n
x
,
min
⁡
=
-
0.725
n
x
,
max
⁡
=
0.91
Normal load factor (g)
n
z
,
min
⁡
=
-
3.2
n
z
,
max
⁡
=
8
Roll angle rate (deg/s)
μ
˙
min
⁡
=
-
30
μ
˙
max
⁡
=
30
Rate of change of tangential load factor (g/s)
n
˙
x
,
min
⁡
=
-
0.2
n
˙
x
,
max
⁡
=
0.2
Rate of change of normal load factor (g/s)
n
˙
z
,
min
⁡
=
-
2
n
˙
z
,
max
⁡
=
2The experimental test environment is a rectangle area of 30 × 40 km2, as shown in Figure 2, where the NFZ1 caused by the severe weather condition is modeled as a cuboid with infinite altitude, and the NFZs caused by enemy threats are modeled as two hemispheres, denoted as THT1 and THT2. According to Figure 1, the shape parameters can be set as follows: for the NFZ1, k
1
=
k
2
=
k
3
=
k
4
=
1, p
=
2, and the length of a side is set as 8 km; for THT1 and THT2, k
1
=
k
2
=
k
4
=
1, k
3
=
2, p
=
6, and the radius of the threats can be set as 7 km and 9 km, respectively. In order to avoid contact, a safe buffer b
s
(
b
s
=
200
m
) is added around the outside edge of each obstacle/threat boundary. All the results presented below are generated using TOMLAB/SNOPT software toolbox on a 2.4 GHz CPU and 2 G RAM computer running Windows XP and MATLAB R2009b. Table 2 summarizes the parameters used in the algorithm. To simplify the problem, the same weight coefficients for the objective function of each UCAV are set, and it is assumed that the target assignment is already completed before.Table 2
Parameter values for the algorithm.
Item
Parameter
Value
Minimum safety radius (m)
d
safe
500
Weighting factors of the objective function
w
t
0.3
w
p
0.4
w
r
0.3
Weighting factors of the robust subobjective
w
r
1
0.5
w
r
2
0.3
w
r
3
0.2
Safe buffer of NFZs (m)
b
s
200Figure 2
Two collision-free UCAV trajectories with arriving simultaneously.


### 5.1. Example  1: Cooperative Trajectories of Two UCAVs Arriving Simultaneously

In this example, two UCAVs cooperatively attack two stationary ground targets while avoiding a series of static obstacles/threats detected and collision with each other en route and meeting aircraft dynamical constraints, especially simultaneous arrival constraint. UCAV1 and UCAV2 start at each initial point, that is, IP1 (10, 2, 2) km and IP2 (17, 2, 2) km. Then they fly into the AARs of two targets, TAR1 (4, 34, 0) km and TAR2 (13, 40, 0) km, respectively. The initial remaining three states and control inputs of each UCAV are preset asρ
r
,
0
1
=
(
V
0
1
,
γ
0
1
,
ψ
0
1
)
=
(
220
m
/
s
,
0
∘
,
9
0
∘
), ρ
r
,
0
2
=
(
V
0
2
,
γ
0
2
,
ψ
0
2
)
=
(
220
m
/
s
,
0
∘
,
14
0
∘
), and u
0
1
=
u
0
2
=
(
μ
0
,
n
x
,
0
,
n
z
,
0
)
=
(
0
∘
,
0
g
,
1
g
), and the terminal remaining three states of each UCAV are predetermined as ρ
r
,
f
1
=
(
V
⌢
f
1
,
γ
⌢
f
1
,
ψ
⌢
f
1
)
=
(
250
m
/
s
,
0
∘
,
11
0
∘
) and ρ
r
,
f
2
=
(
V
⌢
f
2
,
γ
⌢
f
2
,
ψ
⌢
f
2
)
=
(
250
m
/
s
,
0
∘
,
16
0
∘. In VMC framework, the reference point ρ
ref for each UCAV is presumed as its IP, and the virtual prey path is predefined as a 2-order curve, determined by the initial and terminal conditions.Figure2 shows the overall collision-free attack trajectories of two UCAVs, whose total flight time is equal, equaling to 124 s. It means that two UCAVs arrive simultaneously. As can be seen from Figure 2, along the resulting smooth trajectories, the UCAVs can avoid collision with all obstacles, threats, and the other en route and then successfully fly into the AARs (fuchsin areas) to perform weapon delivery missions. In addition, the approximate weapon trajectories are drawn to illustrate the attack process. The PCP histories of two UCAVs are shown in Figure 3, which are smooth and lie in the range from 0 to 1. And the time histories of the UCAVs’ states (
V
,
γ
,
ψ
) and control inputs (
μ
,
n
x
,
n
z
) are shown in Figure 4. It is clear that the constraints on these variables, especially the cooperative constraints, are all satisfied (see Tables 1 and 2), which means that the resulting trajectories are feasible and safe.Figure 3
PCP histories of two UCAVs.Figure 4
Example
1—state and control time histories of two UCAVs.


### 5.2. Example  2: Cooperative Trajectories of Multi-UCAV Arriving in Sequence

In this example, 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 set to 20 s and 30 s, respectively, denoted asΔ
12
=
20
sand Δ
23
=
30
s. The coordinates of the initial points IP1, IP2 and two targets TAR1, TAR2 are the same with those in Example  1, and the third initial point is IP3 (4, 2, 2) km. The result of the previously finished target assignment is that the UCAV1 and UCAV2 attack TAR1 and UCAV3 attacks the other one. The initial remaining three states and control inputs of the UCAV3 are preset as ρ
r
,
0
3
=
(
V
0
3
,
γ
0
3
,
ψ
0
3
)
=
(
220
m
/
s
,
0
∘
,
4
5
∘
) and u
0
3
=
(
μ
0
,
n
x
,
0
,
n
z
,
0
)
=
(
0
∘
,
0
g
,
1
g
), and the terminal remaining three states of the UCAV3 are predetermined as ρ
r
,
f
3
=
(
V
⌢
f
3
,
γ
⌢
f
3
,
ψ
⌢
f
3
)
=
(
220
m
/
s
,
0
∘
,
16
0
∘
). In VMC framework, the reference point and the virtual prey path are predetermined, using the same strategy as the one in Example  1.The overall collision-free attack trajectories of multi-UCAV are shown in Figure5, and the total flight time of three UCAVs is 100 s, 120 s, and 150 s, respectively. It can be clearly demonstrated that UCAVs can avoid all obstacles or threats and successfully fly into the AARs in sequence to perform the weapon delivery. The PCP histories of three UCAVs are shown in Figure 6. And Figure 7 shows the distance between each pair of UCAVs. From it, one can find that the minimum distance between each UCAV is more than the minimal safety radius of d
safe
(
d
safe
=
500
m
); that is, the UCAVs can avoid collision with each other en route. Figure 8 shows the time histories of the UCAVs’ states (
V
,
γ
,
ψ
) and control inputs (
μ
,
n
x
,
n
z
). Obviously, the resulting trajectories are feasible and safe, since all the constraints listed in Tables 1 and 2 are satisfied.Figure 5
Three collision-free UCAV trajectories with arriving in sequence.Figure 6
Path control parameter histories of three UCAVs.Figure 7
Distance between each pair of  UCAVs.Figure 8
Example
2—state and control time histories of three UCAVs.In order to validate the convergence rate of VMC based algorithm in cooperative trajectory planning for multi-UCAV further, the comparative experiments are carried out with different numbers of UCAVs (N
v is set as 1, 2, 3, 5, 8, and 10, resp.) and the same number of nodes (N
c
=
20). The performance comparison of multi-UCAV is summarized in Table 3. As can be seen clearly, along with the increase of the number of UCAVs, the average planning time is rising at a slow rate, which is acceptable for the preplanning applications. Particularly, for one UCAV, the average planning time is 2.21 s. When the number of UCAVs increases to 10, the average planning time increases by only more than quadruple of that for one UCAV, which is 9.62 s.Table 3
Performance comparison with the increase of the number of UCAVs.
Number of UCAVs (N
v)
Number of optimization variables
Average planning time (s)
Ratio of the planning time (%)
1
(
20
×
1
+
1
)
×
1
=
21
2.21
—
2
(
20
×
1
+
1
)
×
2
=
42
2.83
128.05
3
(
20
×
1
+
1
)
×
3
=
63
3.32
150.23
5
(
20
×
1
+
1
)
×
5
=
105
4.93
223.08
8
(
20
×
1
+
1
)
×
8
=
168
7.70
348.42
10
(
20
×
1
+
1
)
×
10
=
210
9.62
435.29


### 5.3. Comparison between VMC Based Method and GPM Based Direct Collocation Method

To further evaluate the performance of the proposed algorithm, a recently developed optimal trajectory generation method, GPM based direct collocation method, is used for comparison, which is basically converting the original optimal control problem to a NLP problem via the GPM directly. In the following comparative experiments, Example  1 is simulated again using the two algorithms, respectively, performed on the same desktop with different numbers of nodes (N
c is set as 10, 15, 20, and 30, resp.). All the results are presented below.Figures9 and 10 show the comparative solutions attained by using the GPM and VMC based methods, respectively. It can be seen that the trajectories and the state and control time histories of these two algorithms are similar. The detailed performance comparison is summarized in Table 4. With comparison, when the number of nodes increases, the number of optimization variables in GPM is much more than VMC. As a result, the computational speed of VMC is more than an order of magnitude faster than GPM, while the optimization performance has only small loss.Table 4
Comparison between GPM and VMC for Example  1.
Method
N
c
Number of optimization variables
Objective function value
Average planning time (s)
Ratio of the planning time (%)
GPM
10
(
10
×
9
+
1
)
×
2
=
182
1.776
12.31
—
15
(
15
×
9
+
1
)
×
2
=
272
1.731
16.28
—
20
(
20
×
9
+
1
)
×
2
=
362
1.698
20.31
—
30
(
30
×
9
+
1
)
×
2
=
542
1.705
28.67
—
VMC
10
(
10
×
1
+
1
)
×
2
=
22
1.851
2.23
18.12
15
(
15
×
1
+
1
)
×
2
=
32
1.840
2.46
15.11
20
(
20
×
1
+
1
)
×
2
=
42
1.769
2.83
13.93
30
(
30
×
1
+
1
)
×
2
=
62
1.731
3.18
11.09Figure 9
Comparison of two collision-free UCAV trajectories with arriving simultaneously.Figure 10
Comparison of state and control time histories.In order to compare the convergence rate and analyze the impacts of the number of UCAVs on the computation time between VMC and GPM further, the comparative experiment is run multiple times with different numbers of UCAVs (N
v is set as an integer from 1 to 10, resp.) and the same number of nodes (N
c
=
20). Figure 11 shows the change of the average planning time required to solve the optimization problem as the number of UCAVs increases. It can be seen that the increase of the number of UCAVs results in an exponential increase of the required computation time for GPM. In contrast, the computation time for VMC based algorithm increases slowly. When the number of UCAVs increases to 10, the number of the optimization variables for GPM is up to 1810. However, for VMC, the number is only 210. Since the required computation time strongly depends on the number of optimization variables, the computational time required by VMC is more than an order of magnitude less than that required by GPM. Maybe it can be concluded that the larger the number of UCAVs is, the more obvious the advantage of VMC to GPM is.Figure 11
Comparison between GPM and VMC with the increase of the number of UCAVs.


## 5.1. Example  1: Cooperative Trajectories of Two UCAVs Arriving Simultaneously

In this example, two UCAVs cooperatively attack two stationary ground targets while avoiding a series of static obstacles/threats detected and collision with each other en route and meeting aircraft dynamical constraints, especially simultaneous arrival constraint. UCAV1 and UCAV2 start at each initial point, that is, IP1 (10, 2, 2) km and IP2 (17, 2, 2) km. Then they fly into the AARs of two targets, TAR1 (4, 34, 0) km and TAR2 (13, 40, 0) km, respectively. The initial remaining three states and control inputs of each UCAV are preset asρ
r
,
0
1
=
(
V
0
1
,
γ
0
1
,
ψ
0
1
)
=
(
220
m
/
s
,
0
∘
,
9
0
∘
), ρ
r
,
0
2
=
(
V
0
2
,
γ
0
2
,
ψ
0
2
)
=
(
220
m
/
s
,
0
∘
,
14
0
∘
), and u
0
1
=
u
0
2
=
(
μ
0
,
n
x
,
0
,
n
z
,
0
)
=
(
0
∘
,
0
g
,
1
g
), and the terminal remaining three states of each UCAV are predetermined as ρ
r
,
f
1
=
(
V
⌢
f
1
,
γ
⌢
f
1
,
ψ
⌢
f
1
)
=
(
250
m
/
s
,
0
∘
,
11
0
∘
) and ρ
r
,
f
2
=
(
V
⌢
f
2
,
γ
⌢
f
2
,
ψ
⌢
f
2
)
=
(
250
m
/
s
,
0
∘
,
16
0
∘. In VMC framework, the reference point ρ
ref for each UCAV is presumed as its IP, and the virtual prey path is predefined as a 2-order curve, determined by the initial and terminal conditions.Figure2 shows the overall collision-free attack trajectories of two UCAVs, whose total flight time is equal, equaling to 124 s. It means that two UCAVs arrive simultaneously. As can be seen from Figure 2, along the resulting smooth trajectories, the UCAVs can avoid collision with all obstacles, threats, and the other en route and then successfully fly into the AARs (fuchsin areas) to perform weapon delivery missions. In addition, the approximate weapon trajectories are drawn to illustrate the attack process. The PCP histories of two UCAVs are shown in Figure 3, which are smooth and lie in the range from 0 to 1. And the time histories of the UCAVs’ states (
V
,
γ
,
ψ
) and control inputs (
μ
,
n
x
,
n
z
) are shown in Figure 4. It is clear that the constraints on these variables, especially the cooperative constraints, are all satisfied (see Tables 1 and 2), which means that the resulting trajectories are feasible and safe.Figure 3
PCP histories of two UCAVs.Figure 4
Example
1—state and control time histories of two UCAVs.


## 5.2. Example  2: Cooperative Trajectories of Multi-UCAV Arriving in Sequence

In this example, 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 set to 20 s and 30 s, respectively, denoted asΔ
12
=
20
sand Δ
23
=
30
s. The coordinates of the initial points IP1, IP2 and two targets TAR1, TAR2 are the same with those in Example  1, and the third initial point is IP3 (4, 2, 2) km. The result of the previously finished target assignment is that the UCAV1 and UCAV2 attack TAR1 and UCAV3 attacks the other one. The initial remaining three states and control inputs of the UCAV3 are preset as ρ
r
,
0
3
=
(
V
0
3
,
γ
0
3
,
ψ
0
3
)
=
(
220
m
/
s
,
0
∘
,
4
5
∘
) and u
0
3
=
(
μ
0
,
n
x
,
0
,
n
z
,
0
)
=
(
0
∘
,
0
g
,
1
g
), and the terminal remaining three states of the UCAV3 are predetermined as ρ
r
,
f
3
=
(
V
⌢
f
3
,
γ
⌢
f
3
,
ψ
⌢
f
3
)
=
(
220
m
/
s
,
0
∘
,
16
0
∘
). In VMC framework, the reference point and the virtual prey path are predetermined, using the same strategy as the one in Example  1.The overall collision-free attack trajectories of multi-UCAV are shown in Figure5, and the total flight time of three UCAVs is 100 s, 120 s, and 150 s, respectively. It can be clearly demonstrated that UCAVs can avoid all obstacles or threats and successfully fly into the AARs in sequence to perform the weapon delivery. The PCP histories of three UCAVs are shown in Figure 6. And Figure 7 shows the distance between each pair of UCAVs. From it, one can find that the minimum distance between each UCAV is more than the minimal safety radius of d
safe
(
d
safe
=
500
m
); that is, the UCAVs can avoid collision with each other en route. Figure 8 shows the time histories of the UCAVs’ states (
V
,
γ
,
ψ
) and control inputs (
μ
,
n
x
,
n
z
). Obviously, the resulting trajectories are feasible and safe, since all the constraints listed in Tables 1 and 2 are satisfied.Figure 5
Three collision-free UCAV trajectories with arriving in sequence.Figure 6
Path control parameter histories of three UCAVs.Figure 7
Distance between each pair of  UCAVs.Figure 8
Example
2—state and control time histories of three UCAVs.In order to validate the convergence rate of VMC based algorithm in cooperative trajectory planning for multi-UCAV further, the comparative experiments are carried out with different numbers of UCAVs (N
v is set as 1, 2, 3, 5, 8, and 10, resp.) and the same number of nodes (N
c
=
20). The performance comparison of multi-UCAV is summarized in Table 3. As can be seen clearly, along with the increase of the number of UCAVs, the average planning time is rising at a slow rate, which is acceptable for the preplanning applications. Particularly, for one UCAV, the average planning time is 2.21 s. When the number of UCAVs increases to 10, the average planning time increases by only more than quadruple of that for one UCAV, which is 9.62 s.Table 3
Performance comparison with the increase of the number of UCAVs.
Number of UCAVs (N
v)
Number of optimization variables
Average planning time (s)
Ratio of the planning time (%)
1
(
20
×
1
+
1
)
×
1
=
21
2.21
—
2
(
20
×
1
+
1
)
×
2
=
42
2.83
128.05
3
(
20
×
1
+
1
)
×
3
=
63
3.32
150.23
5
(
20
×
1
+
1
)
×
5
=
105
4.93
223.08
8
(
20
×
1
+
1
)
×
8
=
168
7.70
348.42
10
(
20
×
1
+
1
)
×
10
=
210
9.62
435.29


## 5.3. Comparison between VMC Based Method and GPM Based Direct Collocation Method

To further evaluate the performance of the proposed algorithm, a recently developed optimal trajectory generation method, GPM based direct collocation method, is used for comparison, which is basically converting the original optimal control problem to a NLP problem via the GPM directly. In the following comparative experiments, Example  1 is simulated again using the two algorithms, respectively, performed on the same desktop with different numbers of nodes (N
c is set as 10, 15, 20, and 30, resp.). All the results are presented below.Figures9 and 10 show the comparative solutions attained by using the GPM and VMC based methods, respectively. It can be seen that the trajectories and the state and control time histories of these two algorithms are similar. The detailed performance comparison is summarized in Table 4. With comparison, when the number of nodes increases, the number of optimization variables in GPM is much more than VMC. As a result, the computational speed of VMC is more than an order of magnitude faster than GPM, while the optimization performance has only small loss.Table 4
Comparison between GPM and VMC for Example  1.
Method
N
c
Number of optimization variables
Objective function value
Average planning time (s)
Ratio of the planning time (%)
GPM
10
(
10
×
9
+
1
)
×
2
=
182
1.776
12.31
—
15
(
15
×
9
+
1
)
×
2
=
272
1.731
16.28
—
20
(
20
×
9
+
1
)
×
2
=
362
1.698
20.31
—
30
(
30
×
9
+
1
)
×
2
=
542
1.705
28.67
—
VMC
10
(
10
×
1
+
1
)
×
2
=
22
1.851
2.23
18.12
15
(
15
×
1
+
1
)
×
2
=
32
1.840
2.46
15.11
20
(
20
×
1
+
1
)
×
2
=
42
1.769
2.83
13.93
30
(
30
×
1
+
1
)
×
2
=
62
1.731
3.18
11.09Figure 9
Comparison of two collision-free UCAV trajectories with arriving simultaneously.Figure 10
Comparison of state and control time histories.In order to compare the convergence rate and analyze the impacts of the number of UCAVs on the computation time between VMC and GPM further, the comparative experiment is run multiple times with different numbers of UCAVs (N
v is set as an integer from 1 to 10, resp.) and the same number of nodes (N
c
=
20). Figure 11 shows the change of the average planning time required to solve the optimization problem as the number of UCAVs increases. It can be seen that the increase of the number of UCAVs results in an exponential increase of the required computation time for GPM. In contrast, the computation time for VMC based algorithm increases slowly. When the number of UCAVs increases to 10, the number of the optimization variables for GPM is up to 1810. However, for VMC, the number is only 210. Since the required computation time strongly depends on the number of optimization variables, the computational time required by VMC is more than an order of magnitude less than that required by GPM. Maybe it can be concluded that the larger the number of UCAVs is, the more obvious the advantage of VMC to GPM is.Figure 11
Comparison between GPM and VMC with the increase of the number of UCAVs.


## 6. Conclusions

This paper is devoted to explore the cooperative collision-free trajectory planning for multiple UCAVs in performing the CA/GTA missions. The main contributions of the paper are as follows. Firstly, the cooperative trajectory planning problem under consideration is mathematically formulated as a cooperative trajectory optimal control problem (CTOCP). In order to consider the weapon delivery constraints, an approximate allowable attack region model is proposed and integrated into the problem formulation. Secondly, a particular virtual motion camouflage approach combining with differential flatness theory, Gauss pseudospectral method, and nonlinear programming is used to develop the trajectory planner for solving CTOCP. The advantages of this method are that it can automatically satisfy all boundary conditions and use any models and any performance indexes, does not require numerical integration of differentiation of the state equations, and transforms the OCP into a NLP problem of very low dimension. Finally, in order to handle the temporal constraints, the notion of the virtual time is introduced into the virtual motion camouflage approach to realize the spatial-temporal cooperation.The proposed approach is validated by numerical examples. The results show that the approach is able to generate both feasible and near-optimal attack trajectories with meeting the spatial and temporal constraints very efficiently. Moreover, the convergence rate and average planning time of the method and optimality of the generated trajectories are evaluated via a detailed comparison with GPM based direct collocation method. The results show that the computational speed of virtual motion camouflage approach is more than an order of magnitude faster than GPM, at small loss of optimality.For the future work, we will analyze some uncertain factors in the true battlefield environment and carry out the research on the real-time cooperative trajectory planning. For the presence of a larger number of UCAVs, we will make further efforts to exploit more efficient trajectory planning algorithm and improve the time cooperative strategy. Moreover, we will try to study another important aspect about how to plan in the opposability battlefield environment.



---

*Source: 748974-2014-01-22.xml*