# Cooperative Solution of Multi-UAV Rendezvous Problem with Network Restrictions

Authors: Qingjie Zhang; Jianwu Tao; Fei Yu; Yingzhen Li; Hongchang Sun; Wei Xu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

By considering the complex networks, the cooperative game based optimal consensus (CGOC) algorithm is proposed to solve the multi-UAV rendezvous problem in the mission area. Firstly, the mathematical description of the rendezvous problem is established, and the solving framework is provided based on  the coordination variables and coordination function. It can decrease the transmission of the redundant information and reduce the influence of the limited network on the task. Secondly, the CGOC algorithm is presented for the UAVs in distributed cooperative manner, which can minimize the overall cost of the multi-UAV system. The CGOC control problem and  the corresponding solving protocol are given by using the cooperative game theory and sensitivity parameter method. Then, the CGOC method of multi-UAV rendezvous problem is proposed, which focuses on the trajectory control of the platform rather than the path planning. Simulation results are given to demonstrate the effectiveness of the proposed CGOC method under complex network conditions and the benefit on the overall optimality and dynamic response.


---

## Body

## 1. Introduction

In order to execute the missions such as simultaneous strike, cooperative reconnaissance, or SEAD (suppression of enemy air defense), multiple UAVs need to arrive at a selected region from different directions. Each UAV plans its path dynamically considering the restriction including enemy radars, missiles, and its own performance. For the sake of the successful task achievement, it requires that all UAVs arrive at the goal position simultaneously or sequentially, which is called rendezvous problem.In respect to the time coordination problem, the method using the coordination function and coordination variables is adopted by [1–4]. The optimal flight control sequence of the UAVs can be obtained to meet certain performance index by this method. It can not only simplify the complexity of solving the problem but also reduce the redundant communication links among the UAVs. But it is still a kind of centralized control method in the view of the command and control manner. After the UAVs flying into the mission area, the communication will be reduced or interrupted. In order to ensure the survivability and increase the probability of the successful task, a distributed control method based on multiagent average consensus algorithm is proposed by [5–9]. It can deal with the case of battlefield environment changes. Due to the adoption of the average consensus algorithm, it is a “compromise” collaborative result. Therefore, the consensus-based flight control sequence is not considered in an optimal manner. To improve the performance of consensus algorithm, literature [10] proposes a noncooperative based optimal consensus (NCOC) algorithm, which can be employed to solve this rendezvous problem. However, the overall cost of multi-UAV system is usually not the minimum, since these UAVs only minimize their own cost function rather than the overall one. This is due to the fact that these UAVs are selfish or noncooperative. Literature [11, 12] focuses on the dynamic response and optimal cost by introducing the outdated state difference and the optimized weighted matrix to the consensus algorithm. This method accelerates the speed of the consensus convergence. Similarly, the consensus algorithm with virtual leader and state predictor has been adopted in [13] to increase the convergence speed of the rendezvous problem.Different from the existing methods using average consensus algorithm [5–9], this paper focuses on the optimal solution rather than the convergence speed [11–13]. Further, the NCOC method [10] is developed, and the cooperative method of solving the rendezvous problem is proposed. The system achieves the overall optimal cost in this process. The novel method can deal with the complex network restriction such as switching topologies, communication delay.The remainder of this paper is organized as follows. Firstly, the mathematical description of rendezvous problem is given. Based on the results of [5–9], the solving framework is proposed using the coordination function and coordination variable. A novel distributed control method using the cooperative game based optimal consensus (CGOC) algorithm is proposed for the “cooperative” UAVs. Finally, numerical experiments and simulation results are illustrated to show the effectiveness and benefit of the proposed method.


## 2. Problem Description

### 2.1. Basic Assumptions and Physical Constraints

This paper focuses on the distributed control method of the rendezvous problem, rather than the flight control of the platform. Hence, we assume that(1)
all UAVs are small size;
(2)
each UAV is equipped with the autopilot, which can track the waypoint automatically;
(3)
mission control and flight control can be decoupled, respectively. Thus, the flight control problem can be simplified. In order to show the physical characteristics of the platform, the flight constraints and related performance parameters are listed as follows: (i)
the maximum and minimum speedv
m
a
x, v
m
i
n;
(ii)
the minimum turning radiusR
m
i
n;
(iii)
the maximum flight timet
m
a
x.


### 2.2.  Incomplete Kinematic Model of UAV

Similar to [1–7, 9], the 2D kinematic model of the UAV is selected to study the rendezvous problem (∀
i
=
1,2
,
…
,
n)(1)
z
˙
i
x
t
=
v
i
t
cos
⁡
φ
i
t
,
z
˙
i
y
t
=
v
i
t
sin
⁡
φ
i
t
,
φ
˙
i
t
=
w
i
t
,where z
i
(
t
)
=
[
z
i
x
(
t
)
;
z
i
y
(
t
)
]
T denotes the position vector of the ith UAV and φ
i
(
t
) is the heading angle. v
i
(
t
) is the velocity and w
i
(
t
) is the changing rate of the heading angular velocity. w
i
(
t
) and v
i
(
t
) satisfy(2)
w
i
t
∈
-
w
h
e
a
d
,
w
h
e
a
d
,
(3)
v
i
t
∈
v
m
i
n
,
v
m
a
x
.The parameters w
h
e
a
d, v
m
i
n, and v
m
a
x are determined by the physical performance.The autopilot of each UAV maintains the expected heading angle and velocity. Its mathematical model can be described by two first-order differential equations(4)
φ
˙
i
t
=
a
φ
φ
i
c
t
-
φ
i
t
,
v
˙
i
t
=
a
v
v
i
c
t
-
v
i
t
.The variable or vector with the superscript c signifies the reference instruction. Parameters a
φ and a
v are the constant coefficients of the heading and velocity channel of the autopilot.


### 2.3. Communication Model

Graph theory is used to describe the communication model among multi-UAV system; see, Figure1. Let G
=
(
V
,
E
,
A
) denote the relationship between multiple UAVs with the set of nodes V
=
{
υ
1
,
…
,
υ
n
}, the set of edges E
⊂
V
×
V, and adjacent matrix A
=
[
a
i
j
]. The node indices belong to a finite index set I
=
{
1,2
,
…
,
n
}. The edge can be depicted by (
υ
i
,
υ
j
), and the value of a
i
j corresponds to the edge of the graph; that is, (
υ
i
,
υ
j
)
∈
E
⇔
a
i
j
>
0. The neighbors set of node υ
i is defined by N
i
=
{
υ
j
∈
V
:
(
υ
i
,
υ
j
)
∈
E
}.Figure 1
The interconnection of multi-UAV system with communication delay and topology varying.Since the phenomenon such as communication delay, interruption, and limited bandwidth may appear in practical application, in this paper, two restrictions will be considered(i)
communication delayτ
i
j: it denotes the transfer delay of the message from the jth to the ith UAV;
(ii)
time varying interconnectionsG
σ
(
t
)
∈
{
G
a
,
G
b
,
…
,
G
N
s
s
}: σ
(
t
)
:
[
0
,
∞
]
→
{
a
,
b
,
…
,
N
s
s
} denotes the switching signal with successive times to describe the topology switches.


### 2.4. Mathematical Model of Rendezvous Problem

The mathematical description of the rendezvous problem will be given as follows. The initial and target position of theith UAV are defined as z
i
0 and z
i
f. Then, the flight trajectory can be defined as(5)
Z
i
T
i
≜
z
i
t
:
0
≤
t
≤
T
i
.If there exist input sequence φ
i
c
(
t
) and v
i
c
t, 0
≤
t
≤
T
i satisfying z
i
(
T
i
)
=
z
i
f, kinematic model (1), and restricts (2), (3), then Z
i
(
T
i
) is the feasible flight path of the ith UAV.For each feasible path, UAV should avoid threats and return to the base safely with adequate fuel. In fact, there is a strong coupling relationship among these restrictions. On one hand, avoiding threats means longer path and higher flight velocity. It demands less residence time of the aircraft in danger. On the other hand, fuel saving means shorter path and lower velocity. Thus, the rendezvous problem is a time coordination problem considering the exposure time in the mission area and fuel consumption.Define the cost function of the pathZ
i
(
T
i
) as(6)
J
Z
i
T
i
=
k
1
J
t
h
r
e
a
t
Z
i
T
i
+
k
2
J
f
u
e
l
Z
i
T
i
=
∑
p
=
1
N
T
∫
0
T
i
k
1
h
p
-
z
i
t
4
d
t
+
∫
0
T
i
k
2
v
i
2
t
d
t
,where k
1
,
k
2
∈
[
0,1
] and k
1
+
k
2
=
1. The threat cost J
t
h
r
e
a
t
(
·
) is determined by the exposure time under the threat radar. If the signals emitted by the radar in all directions are the same, it is proportional to the fourth power of the distance between the UAV and the threat. In (6), h
p is the position of the pth threat, and H
=
{
h
1
,
h
2
,
…
,
h
N
T
} is the threat set. Because the fuel consumption rate is determined by the air drag torque, the fuel cost J
f
u
e
l
(
·
) is proportional to the square of the velocity. c
2
∈
R is the proportion factor.With (5) and (6), the rendezvous problem is equivalent to the following global optimization equation:(7)
min
⁡
∑
i
=
1
n
J
Z
i
T
i
s.t.
Z
i
T
i
is a feasible path
T
1
=
T
2
=
⋯
=
T
n
=
T
E
T
A
.Figure2 demonstrates the rendezvous problem of three UAVs (the dashed line denotes the feasible path). For example, when number 2 UAV detects a new unknown threat, all UAVs need to negotiate with each other and determine a new ETA (estimated time of arrival). Then, a new path is replanned; see the solid line in Figure 2. The new path will ensure the UAVs arrive at their destination simultaneously. It is obvious that the costs of number 1 and number 3 UAV are not optimal for that path. But the bigger ETA negotiated will accomplish the mission. Therefore, rendezvous problem is a typical cooperative control problem of multi-UAV, including two aspects. One is path planning. Each UAV plans its path considering some constraints including radars, missile threats, and platform performance. The other one is trajectory control. Each UAV arrives at its destination simultaneously through adjusting its velocity and heading angle.Figure 2
Typical scenery of multi-UAV rendezvous problem.Obviously, there are also some challenges for solving the global optimization problem. Firstly, the expressions (1), (2), and (3) are noncomplete constraints of the motion model. So, it is difficult to generate the feasible path. Secondly, the gradient optimization technique is very sensitive to the initial path. Finally, time constraints mean that all UAVs should plan their path simultaneously. Therefore, the suboptimal or feasible solution of problem (7) is discussed.


## 2.1. Basic Assumptions and Physical Constraints

This paper focuses on the distributed control method of the rendezvous problem, rather than the flight control of the platform. Hence, we assume that(1)
all UAVs are small size;
(2)
each UAV is equipped with the autopilot, which can track the waypoint automatically;
(3)
mission control and flight control can be decoupled, respectively. Thus, the flight control problem can be simplified. In order to show the physical characteristics of the platform, the flight constraints and related performance parameters are listed as follows: (i)
the maximum and minimum speedv
m
a
x, v
m
i
n;
(ii)
the minimum turning radiusR
m
i
n;
(iii)
the maximum flight timet
m
a
x.


## 2.2.  Incomplete Kinematic Model of UAV

Similar to [1–7, 9], the 2D kinematic model of the UAV is selected to study the rendezvous problem (∀
i
=
1,2
,
…
,
n)(1)
z
˙
i
x
t
=
v
i
t
cos
⁡
φ
i
t
,
z
˙
i
y
t
=
v
i
t
sin
⁡
φ
i
t
,
φ
˙
i
t
=
w
i
t
,where z
i
(
t
)
=
[
z
i
x
(
t
)
;
z
i
y
(
t
)
]
T denotes the position vector of the ith UAV and φ
i
(
t
) is the heading angle. v
i
(
t
) is the velocity and w
i
(
t
) is the changing rate of the heading angular velocity. w
i
(
t
) and v
i
(
t
) satisfy(2)
w
i
t
∈
-
w
h
e
a
d
,
w
h
e
a
d
,
(3)
v
i
t
∈
v
m
i
n
,
v
m
a
x
.The parameters w
h
e
a
d, v
m
i
n, and v
m
a
x are determined by the physical performance.The autopilot of each UAV maintains the expected heading angle and velocity. Its mathematical model can be described by two first-order differential equations(4)
φ
˙
i
t
=
a
φ
φ
i
c
t
-
φ
i
t
,
v
˙
i
t
=
a
v
v
i
c
t
-
v
i
t
.The variable or vector with the superscript c signifies the reference instruction. Parameters a
φ and a
v are the constant coefficients of the heading and velocity channel of the autopilot.


## 2.3. Communication Model

Graph theory is used to describe the communication model among multi-UAV system; see, Figure1. Let G
=
(
V
,
E
,
A
) denote the relationship between multiple UAVs with the set of nodes V
=
{
υ
1
,
…
,
υ
n
}, the set of edges E
⊂
V
×
V, and adjacent matrix A
=
[
a
i
j
]. The node indices belong to a finite index set I
=
{
1,2
,
…
,
n
}. The edge can be depicted by (
υ
i
,
υ
j
), and the value of a
i
j corresponds to the edge of the graph; that is, (
υ
i
,
υ
j
)
∈
E
⇔
a
i
j
>
0. The neighbors set of node υ
i is defined by N
i
=
{
υ
j
∈
V
:
(
υ
i
,
υ
j
)
∈
E
}.Figure 1
The interconnection of multi-UAV system with communication delay and topology varying.Since the phenomenon such as communication delay, interruption, and limited bandwidth may appear in practical application, in this paper, two restrictions will be considered(i)
communication delayτ
i
j: it denotes the transfer delay of the message from the jth to the ith UAV;
(ii)
time varying interconnectionsG
σ
(
t
)
∈
{
G
a
,
G
b
,
…
,
G
N
s
s
}: σ
(
t
)
:
[
0
,
∞
]
→
{
a
,
b
,
…
,
N
s
s
} denotes the switching signal with successive times to describe the topology switches.


## 2.4. Mathematical Model of Rendezvous Problem

The mathematical description of the rendezvous problem will be given as follows. The initial and target position of theith UAV are defined as z
i
0 and z
i
f. Then, the flight trajectory can be defined as(5)
Z
i
T
i
≜
z
i
t
:
0
≤
t
≤
T
i
.If there exist input sequence φ
i
c
(
t
) and v
i
c
t, 0
≤
t
≤
T
i satisfying z
i
(
T
i
)
=
z
i
f, kinematic model (1), and restricts (2), (3), then Z
i
(
T
i
) is the feasible flight path of the ith UAV.For each feasible path, UAV should avoid threats and return to the base safely with adequate fuel. In fact, there is a strong coupling relationship among these restrictions. On one hand, avoiding threats means longer path and higher flight velocity. It demands less residence time of the aircraft in danger. On the other hand, fuel saving means shorter path and lower velocity. Thus, the rendezvous problem is a time coordination problem considering the exposure time in the mission area and fuel consumption.Define the cost function of the pathZ
i
(
T
i
) as(6)
J
Z
i
T
i
=
k
1
J
t
h
r
e
a
t
Z
i
T
i
+
k
2
J
f
u
e
l
Z
i
T
i
=
∑
p
=
1
N
T
∫
0
T
i
k
1
h
p
-
z
i
t
4
d
t
+
∫
0
T
i
k
2
v
i
2
t
d
t
,where k
1
,
k
2
∈
[
0,1
] and k
1
+
k
2
=
1. The threat cost J
t
h
r
e
a
t
(
·
) is determined by the exposure time under the threat radar. If the signals emitted by the radar in all directions are the same, it is proportional to the fourth power of the distance between the UAV and the threat. In (6), h
p is the position of the pth threat, and H
=
{
h
1
,
h
2
,
…
,
h
N
T
} is the threat set. Because the fuel consumption rate is determined by the air drag torque, the fuel cost J
f
u
e
l
(
·
) is proportional to the square of the velocity. c
2
∈
R is the proportion factor.With (5) and (6), the rendezvous problem is equivalent to the following global optimization equation:(7)
min
⁡
∑
i
=
1
n
J
Z
i
T
i
s.t.
Z
i
T
i
is a feasible path
T
1
=
T
2
=
⋯
=
T
n
=
T
E
T
A
.Figure2 demonstrates the rendezvous problem of three UAVs (the dashed line denotes the feasible path). For example, when number 2 UAV detects a new unknown threat, all UAVs need to negotiate with each other and determine a new ETA (estimated time of arrival). Then, a new path is replanned; see the solid line in Figure 2. The new path will ensure the UAVs arrive at their destination simultaneously. It is obvious that the costs of number 1 and number 3 UAV are not optimal for that path. But the bigger ETA negotiated will accomplish the mission. Therefore, rendezvous problem is a typical cooperative control problem of multi-UAV, including two aspects. One is path planning. Each UAV plans its path considering some constraints including radars, missile threats, and platform performance. The other one is trajectory control. Each UAV arrives at its destination simultaneously through adjusting its velocity and heading angle.Figure 2
Typical scenery of multi-UAV rendezvous problem.Obviously, there are also some challenges for solving the global optimization problem. Firstly, the expressions (1), (2), and (3) are noncomplete constraints of the motion model. So, it is difficult to generate the feasible path. Secondly, the gradient optimization technique is very sensitive to the initial path. Finally, time constraints mean that all UAVs should plan their path simultaneously. Therefore, the suboptimal or feasible solution of problem (7) is discussed.


## 3. Distributed Solution of the Rendezvous Problem

### 3.1. Distributed Solving Structure

#### 3.1.1. Coordination Function and Coordination Variable

This method was proposed to reduce the communication cost and the difficulty of the problem by [1–4, 14]. Coordination variable is the minimum amount of information to coordinate multi-UAV system. Coordination function is the performance of the system achieving effective coordination. The basic idea is listed as follows [3].For theith UAV, the state space description of the battlefield is defined as χ
i and the state as x
i
∈
χ
i. The decision of each UAV would affect the overall cost. U
i
(
x
i
) stands for the feasible decision set, and u
i
∈
U
i
(
x
i
) is the decision variable.In order to achieve effective coordination, there is a minimum amount of information, namely, the coordination variableθ
(
θ
∈
R
c
). If the variable θ and the corresponding function are known to each UAV, the system can achieve collaborative behavior.


#### 3.1.2.  Coordination Variable Selection

Distributed structure for solving multi-UAV rendezvous problem includes(i)
Waypoint Planner (WP): obtain the waypoint sequence considering minimizing the threat and fuel consumption;
(ii)
Kinematic Trajectory Smoothing (TS): generate the fine trajectory in accordance with the kinematic model (1) and the UAV platform physical constraints (2), (3);
(iii)
Distributed Coordinator (DC): receive the ETA from its adjacent UAVs and adjust its ETA by consensus algorithm. The reference velocity instructionv
i
c is generated according to the ETA. Get the reference heading instruction φ
i
c according to the trajectory point got from DS module. Send v
i
c, φ
i
c to AP module;
(iv)
autopilot (AP): ensure the aircraft fly to the destination position according tov
i
c and φ
i
c;
(v)
kinematic model (DM): describe the kinematic characteristics of the platform. Assume that it satisfies formula (1).Obviously, how to obtain the reference velocity is the key to DC module. Letd
E
T
A
,
i
(
t
) denote the remaining path length to the prespecified target location of the ith UAV. Then, d
˙
E
T
A
,
i
(
t
)
=
-
v
i
(
t
). Given d
E
T
A
,
i
(
t
) and v
i
(
t
), the ETA can be calculated as follows:(8)
T
E
T
A
,
i
t
=
d
E
T
A
,
i
t
v
i
t
.The total flight time of theith UAV is(9)
T
E
T
A
,
i
t
=
t
+
T
E
T
A
,
i
t
,
t
∈
t
0
,
t
f
.The first derivative of the T
E
T
A
,
i
(
t
) is calculated as(10)
T
˙
E
T
A
,
i
t
=
1
+
v
i
t
d
˙
l
e
n
,
i
t
-
d
l
e
n
,
i
t
v
˙
i
t
v
i
2
t
=
-
T
E
T
A
,
i
t
a
V
v
i
c
t
-
v
i
t
v
i
t
.Thus, the reference velocity can be obtained(11)
v
i
c
t
=
v
i
t
+
v
i
t
T
˙
E
T
A
,
i
t
a
V
T
E
T
A
,
i
t
.For this rendezvous problem, there are several advantages when the ETA time is selected as the coordination variable: (a) reducing the redundant information negotiating with each other, (b) lowering the difficulty of solving the optimization problem, (c) increasing the dynamic response capacity of multiple UAVs, and (d) cutting down the influence of the restricted networks on the solving structure.


## 3.1. Distributed Solving Structure

### 3.1.1. Coordination Function and Coordination Variable

This method was proposed to reduce the communication cost and the difficulty of the problem by [1–4, 14]. Coordination variable is the minimum amount of information to coordinate multi-UAV system. Coordination function is the performance of the system achieving effective coordination. The basic idea is listed as follows [3].For theith UAV, the state space description of the battlefield is defined as χ
i and the state as x
i
∈
χ
i. The decision of each UAV would affect the overall cost. U
i
(
x
i
) stands for the feasible decision set, and u
i
∈
U
i
(
x
i
) is the decision variable.In order to achieve effective coordination, there is a minimum amount of information, namely, the coordination variableθ
(
θ
∈
R
c
). If the variable θ and the corresponding function are known to each UAV, the system can achieve collaborative behavior.


### 3.1.2.  Coordination Variable Selection

Distributed structure for solving multi-UAV rendezvous problem includes(i)
Waypoint Planner (WP): obtain the waypoint sequence considering minimizing the threat and fuel consumption;
(ii)
Kinematic Trajectory Smoothing (TS): generate the fine trajectory in accordance with the kinematic model (1) and the UAV platform physical constraints (2), (3);
(iii)
Distributed Coordinator (DC): receive the ETA from its adjacent UAVs and adjust its ETA by consensus algorithm. The reference velocity instructionv
i
c is generated according to the ETA. Get the reference heading instruction φ
i
c according to the trajectory point got from DS module. Send v
i
c, φ
i
c to AP module;
(iv)
autopilot (AP): ensure the aircraft fly to the destination position according tov
i
c and φ
i
c;
(v)
kinematic model (DM): describe the kinematic characteristics of the platform. Assume that it satisfies formula (1).Obviously, how to obtain the reference velocity is the key to DC module. Letd
E
T
A
,
i
(
t
) denote the remaining path length to the prespecified target location of the ith UAV. Then, d
˙
E
T
A
,
i
(
t
)
=
-
v
i
(
t
). Given d
E
T
A
,
i
(
t
) and v
i
(
t
), the ETA can be calculated as follows:(8)
T
E
T
A
,
i
t
=
d
E
T
A
,
i
t
v
i
t
.The total flight time of theith UAV is(9)
T
E
T
A
,
i
t
=
t
+
T
E
T
A
,
i
t
,
t
∈
t
0
,
t
f
.The first derivative of the T
E
T
A
,
i
(
t
) is calculated as(10)
T
˙
E
T
A
,
i
t
=
1
+
v
i
t
d
˙
l
e
n
,
i
t
-
d
l
e
n
,
i
t
v
˙
i
t
v
i
2
t
=
-
T
E
T
A
,
i
t
a
V
v
i
c
t
-
v
i
t
v
i
t
.Thus, the reference velocity can be obtained(11)
v
i
c
t
=
v
i
t
+
v
i
t
T
˙
E
T
A
,
i
t
a
V
T
E
T
A
,
i
t
.For this rendezvous problem, there are several advantages when the ETA time is selected as the coordination variable: (a) reducing the redundant information negotiating with each other, (b) lowering the difficulty of solving the optimization problem, (c) increasing the dynamic response capacity of multiple UAVs, and (d) cutting down the influence of the restricted networks on the solving structure.


## 3.1.1. Coordination Function and Coordination Variable

This method was proposed to reduce the communication cost and the difficulty of the problem by [1–4, 14]. Coordination variable is the minimum amount of information to coordinate multi-UAV system. Coordination function is the performance of the system achieving effective coordination. The basic idea is listed as follows [3].For theith UAV, the state space description of the battlefield is defined as χ
i and the state as x
i
∈
χ
i. The decision of each UAV would affect the overall cost. U
i
(
x
i
) stands for the feasible decision set, and u
i
∈
U
i
(
x
i
) is the decision variable.In order to achieve effective coordination, there is a minimum amount of information, namely, the coordination variableθ
(
θ
∈
R
c
). If the variable θ and the corresponding function are known to each UAV, the system can achieve collaborative behavior.


## 3.1.2.  Coordination Variable Selection

Distributed structure for solving multi-UAV rendezvous problem includes(i)
Waypoint Planner (WP): obtain the waypoint sequence considering minimizing the threat and fuel consumption;
(ii)
Kinematic Trajectory Smoothing (TS): generate the fine trajectory in accordance with the kinematic model (1) and the UAV platform physical constraints (2), (3);
(iii)
Distributed Coordinator (DC): receive the ETA from its adjacent UAVs and adjust its ETA by consensus algorithm. The reference velocity instructionv
i
c is generated according to the ETA. Get the reference heading instruction φ
i
c according to the trajectory point got from DS module. Send v
i
c, φ
i
c to AP module;
(iv)
autopilot (AP): ensure the aircraft fly to the destination position according tov
i
c and φ
i
c;
(v)
kinematic model (DM): describe the kinematic characteristics of the platform. Assume that it satisfies formula (1).Obviously, how to obtain the reference velocity is the key to DC module. Letd
E
T
A
,
i
(
t
) denote the remaining path length to the prespecified target location of the ith UAV. Then, d
˙
E
T
A
,
i
(
t
)
=
-
v
i
(
t
). Given d
E
T
A
,
i
(
t
) and v
i
(
t
), the ETA can be calculated as follows:(8)
T
E
T
A
,
i
t
=
d
E
T
A
,
i
t
v
i
t
.The total flight time of theith UAV is(9)
T
E
T
A
,
i
t
=
t
+
T
E
T
A
,
i
t
,
t
∈
t
0
,
t
f
.The first derivative of the T
E
T
A
,
i
(
t
) is calculated as(10)
T
˙
E
T
A
,
i
t
=
1
+
v
i
t
d
˙
l
e
n
,
i
t
-
d
l
e
n
,
i
t
v
˙
i
t
v
i
2
t
=
-
T
E
T
A
,
i
t
a
V
v
i
c
t
-
v
i
t
v
i
t
.Thus, the reference velocity can be obtained(11)
v
i
c
t
=
v
i
t
+
v
i
t
T
˙
E
T
A
,
i
t
a
V
T
E
T
A
,
i
t
.For this rendezvous problem, there are several advantages when the ETA time is selected as the coordination variable: (a) reducing the redundant information negotiating with each other, (b) lowering the difficulty of solving the optimization problem, (c) increasing the dynamic response capacity of multiple UAVs, and (d) cutting down the influence of the restricted networks on the solving structure.


## 4. Optimal Consensus Based on Cooperative Game

### 4.1. Cooperative Game Theory

Considering a team withn players, the quadratic function (12) is constructed to describe the cost of the player in the team (12)
J
i
=
1
2
∫
0
T
x
i
T
Q
i
x
i
+
u
i
T
R
i
u
i
d
t
,where Q
i and R
i are positive definite matrices with proper dimension. Each player satisfies the following dynamical model:(13)
x
˙
=
A
x
+
B
1
u
1
+
⋯
+
B
n
u
n
,
x
0
=
x
0
,in which A and B
i
(
i
=
1
,
…
,
n
) are constant matrices with proper dimensions and u
i is the decision of the ith player. The state variable x can be reflected by the other player’s decision in the minimization process. That means that the players may have conflicting interests. If a player decides to minimize its cost in a noncooperative manner, the decision chosen by the ith player can increase the cost of the others due to the coupling relationship. However, if the players decide to cooperate, individual cost may be minimized and hence we can get a smaller team cost. This will result in the set of Pareto-efficient solutions. For the set of inequalities J
i
U
≤
J
i
U
∗
,
i
=
1,2
,
…
,
n, if there is not at least one solution U, U
∗
=
[
u
1
,
u
2
,
…
,
u
n
] is called Pareto-efficient solution and the corresponding costs J
∗
=
[
J
1
(
U
∗
)
,
J
2
(
U
∗
)
,
…
,
J
n
(
U
∗
)
] are Pareto solution.The solution of the minimization problem cannot be dominated by any other solution(14)
U
∗
α
=
arg
⁡
min
U
∈
U
⁡
∑
i
=
1
n
α
i
J
i
U
,where α
=
(
α
1
,
…
,
α
n
) and ∑
i
=
1
n
α
i
=
1. It is a set of Pareto-efficient solutions of the above problem. The solutions are the functions of the parameter α. The final solution should be selected according to an axiomatic approach as our decision for the cooperative problem. The Nash-bargaining solution is selected in this paper(15)
α
=
arg
⁡
max
α
⁡
∏
i
=
1
n
J
i
α
,
U
∗
-
J
t
p
,
i
,where J
t
p
,
i
(
i
=
1
,
…
,
n
) are the individual costs calculated by using the noncooperation solution that is obtained by minimizing the cost (12).The coefficientα can be obtained (see Theorem  6.10 in [15]) by(16)
α
j
∗
=
∏
i
≠
j
J
i
∗
α
∗
,
U
∗
-
J
t
p
,
i
∑
i
=
1
n
∏
k
≠
i
J
k
∗
α
∗
,
U
∗
-
J
t
p
,
k
.


### 4.2. CGOC Problem Description

When the multiagent systems achieve consensus, all agents get the same values. Thus, the individual cost is defined as(17)
J
i
c
=
1
2
·
∫
0
T
∑
j
∈
N
i
θ
i
-
θ
j
T
Q
i
j
θ
i
-
θ
j
+
u
i
T
R
i
u
i
d
t
,where Q
i
j and R
i are symmetric positive definite matrices.The overall cost of the multiagent systems is obtained by weighted individual cost(18)
J
c
g
=
∑
i
=
1
n
α
i
J
i
c
U
=
1
2
∫
0
∞
θ
T
Q
c
g
θ
+
U
T
R
c
g
U
d
t
,where the coefficient matrices are(19)
Q
c
g
=
Q
i
j
n
×
n
,
Q
i
i
=
∑
j
∈
N
i
α
j
Q
j
i
+
∑
k
∈
N
i
α
i
Q
i
k
,
Q
i
j
=
-
α
i
Q
i
j
-
α
j
Q
j
i
j
∈
N
i
0
j
∉
N
i
,
R
c
g
=
diag
⁡
α
1
R
1
,
…
,
α
n
R
n
.According to cooperative game theory, the smallest overall costJ
c
g can be obtained when the multiagent system achieve consensus. Taking the communication delay (τ
1
=
⋯
=
τ
m
=
τ) into account, the dynamic model of the agent system can be described as follows:(20)
θ
˙
t
=
A
c
g
θ
t
-
τ
+
B
c
g
U
t
,in which θ
=
[
θ
1
T
,
…
,
θ
n
T
]
T is the state vector. A
c
g
=
[
A
i
j
]
n
×
n, A
i
i
=
0, A
i
j
=
-
(
1
/
2
)
l
i
j
H
v
,
i
-
1
Q
i
j, and B
c
g
=
I. H
v
,
i is the solution of the Riccati equation (21), and l
i
j is the element of the Laplacian matrix, which is used to describe the interconnection of the communication(21)
-
H
v
,
i
T
R
i
-
1
H
v
,
i
+
∑
j
∈
N
i
Q
i
j
=
0
.The cooperative game based optimal consensus (CGOC) can be described as the following optimization problem:(22)
min
⁡
U
→
U
J
c
g
s.t.
θ
˙
t
=
A
c
g
θ
t
-
τ
+
B
c
g
U
t
.According to cooperative game theory, the Pareto-optimal solution set can be obtained by(23)
U
∗
α
=
arg
⁡
min
U
∈
U
⁡
∑
i
=
1
n
α
i
J
i
c
U
.It is easy to see that Pareto-efficient solution is a function of the parameterα. Further, the unique Nash bargaining solution is calculated from the above set. Consider(24)
α
=
arg
⁡
max
α
⁡
∏
i
=
1
n
J
i
c
α
,
U
∗
-
J
t
p
,
i
,in which J
t
p
,
i is the individual cost calculated using the noncooperative strategy (see [10])(25)
u
t
p
,
i
=
1
2
H
v
,
i
-
1
∑
j
∈
N
i
Q
i
j
T
θ
j
t
-
τ
i
j
-
R
i
-
1
H
v
,
i
t
θ
i
.When the multiagent systems achieve consensus, the optimal parameterα
∗ is determined. In the following, we will discuss how to get the optimal control strategy of problem (22).


### 4.3. CGOC Problem Solving

Due to the existence of the communication delay in the dynamic model (20) of the agent, it is difficult to get the exact solution or a numerical solution of such a problem. The solving method of the two-point boundary problem with delay and ahead term is adopted in this section [16].Firstly, introduce the Lagrangian operatorλ
(
t
) and define the Hamilton function(26)
H
c
g
θ
,
λ
,
U
≜
1
2
θ
T
Q
c
g
θ
+
U
T
R
c
g
U
+
λ
T
t
A
c
g
θ
t
-
τ
+
B
c
g
U
t
.By the optimal control theory, we can get(27)
θ
˙
t
=
∇
λ
H
c
g
,
t
∈
0
,
T
,
λ
˙
t
=
-
∇
θ
H
c
g
-
∇
θ
τ
H
c
g
t
∈
0
,
T
-
τ
-
∇
θ
H
c
g
t
∈
T
-
τ
,
T
,where ∇
λ
H
c
g, ∇
θ
H
c
g, and ∇
θ
τ
H
c
g are the gradients at λ, θ, and θ
τ. By formula (26), formulas (27) can be rewritten as(28)
θ
˙
t
=
A
c
g
θ
t
-
τ
+
B
c
g
U
t
,
t
∈
0
,
T
,
λ
˙
t
=
-
Q
c
g
θ
t
-
A
c
g
T
λ
t
+
τ
t
∈
0
,
T
-
τ
-
Q
c
g
θ
t
t
∈
T
-
τ
,
T
,in which the initial value is θ
t
=
ϑ
t
,
t
∈
[
-
τ
,
0
]. Thus, the optimal control strategy is(29)
U
t
=
-
R
c
g
-
1
B
c
g
T
λ
t
,
t
∈
0
,
T
.Rewriteθ
˙
(
t
) and λ
˙
(
t
) as follows:(30)
θ
˙
t
=
A
c
g
θ
t
+
B
c
g
U
t
-
A
c
g
θ
t
-
θ
t
-
τ
,
t
∈
0
,
T
,
λ
˙
t
=
-
Q
c
g
θ
t
-
A
c
g
T
λ
t
-
A
c
g
T
λ
t
+
τ
-
λ
t
,
t
∈
0
,
T
-
τ
-
Q
c
g
θ
t
-
A
c
g
T
λ
t
,
t
∈
T
-
τ
,
T
.The sensitivity parameterɛ is introduced; we have(31)
θ
˙
t
,
ɛ
=
A
c
g
θ
t
,
ɛ
+
B
c
g
U
t
,
ɛ
-
ɛ
A
c
g
θ
t
,
ɛ
-
θ
t
-
τ
,
ɛ
,
t
∈
0
,
T
,
λ
˙
t
,
ɛ
=
-
Q
c
g
θ
t
,
ɛ
-
A
c
g
T
λ
t
,
ɛ
-
ɛ
A
c
g
T
λ
t
+
τ
,
ɛ
-
λ
t
,
ɛ
,
t
∈
0
,
T
-
τ
-
Q
c
g
θ
t
,
ɛ
-
A
c
g
T
λ
t
,
ɛ
,
t
∈
T
-
τ
,
T
,
U
t
,
ɛ
=
-
R
c
g
-
1
B
c
g
T
λ
t
,
ɛ
,
k
∈
0
,
T
,and, here, ɛ is a scalar and it satisfies 0
≤
ɛ
≤
1.Assume thatθ
t
,
ɛ, λ
(
t
,
ɛ
), and U
(
t
,
ɛ
) are differentiable at ɛ
=
0, and their Maclaurin series converges at ɛ
=
1
(32)
θ
t
,
ɛ
=
∑
i
=
0
∞
1
i
!
θ
i
t
,
λ
t
,
ɛ
=
∑
i
=
0
∞
1
i
!
λ
i
t
,
U
t
,
ɛ
=
∑
i
=
0
∞
1
i
!
U
i
t
,where θ
i
t
=
lim
ɛ
→
0
⁡
(
∂
i
θ
(
t
,
ɛ
)
/
∂
ɛ
i
), λ
i
(
t
)
=
l
i
m
ɛ
→
0
(
∂
i
λ
(
t
,
ɛ
)
/
∂
ɛ
i
), and U
i
(
t
)
=
l
i
m
ɛ
→
0
(
∂
i
U
(
t
,
ɛ
)
/
∂
ɛ
i
).Thus, the suboptimal control strategy of the above problem (31) is equivalent to the sum of the solutions of the two-point boundary problem, including the 0th-order and the ith-order problems.(i)
Two-point boundary problem with the 0th-order(33)
θ
˙
0
t
=
A
c
g
θ
0
t
+
B
c
g
U
0
t
,
(34)
U
0
t
=
-
R
c
g
-
1
B
c
g
T
λ
0
t
,
(35)
λ
˙
0
t
=
-
Q
c
g
θ
0
t
-
A
c
g
T
λ
0
t
.
(ii)
Two-point boundary problem with theith-order (36)
θ
˙
i
t
=
A
c
g
θ
i
t
+
B
c
g
U
i
t
+
i
A
c
g
θ
i
-
1
t
-
θ
i
-
1
t
-
τ
,
(37)
U
i
t
=
-
R
c
g
-
1
B
c
g
T
λ
i
t
,
(38)
λ
˙
i
t
=
-
Q
c
g
θ
i
t
-
A
c
g
T
λ
i
t
-
A
c
g
T
λ
i
-
1
t
+
τ
-
λ
i
-
1
t
,
t
∈
0
,
T
-
τ
-
Q
c
g
θ
i
t
-
A
c
g
T
λ
i
t
,
t
∈
T
-
τ
,
T
.


#### 4.3.1. The Solution of the 0th-Order Problem

Assume that the Lagrangian operator is(39)
λ
0
t
≜
X
θ
0
t
,where X is a matrix with proper dimension.Substituting formula (34) with (39), we get(40)
U
0
t
=
-
R
c
g
-
1
B
c
g
T
X
θ
0
t
,
t
∈
0
,
T
,
(41)
θ
˙
0
t
=
A
c
g
θ
0
t
+
B
c
g
U
0
t
=
A
c
g
θ
0
t
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
0
t
=
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
0
t
,
t
∈
0
,
T
.Moreover, the Lagrangian operator should satisfy(42)
λ
˙
0
t
=
X
θ
˙
0
t
=
X
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
0
t
,
t
∈
0
,
T
.On the other hand, formula (35) can be rewritten as(43)
λ
˙
0
t
=
-
Q
c
g
θ
0
t
-
A
c
g
T
X
θ
0
t
=
-
Q
c
g
-
A
c
g
T
X
θ
0
t
.Comparing (42) with (43), the Riccati matrices equation can be obtained as(44)
A
c
g
T
X
+
X
A
c
g
-
X
B
c
g
R
c
g
-
1
B
c
g
T
X
+
Q
c
g
=
0
.Therefore, the optimal control strategyU
(
0
) of the 0th-order problem can be obtained through calculating X and θ
(
0
).


#### 4.3.2.  The Solution of theith-Order Problem

Similar to the0th-order problem, the Lagrangian operator is selected as(45)
λ
i
t
≜
X
θ
i
t
+
ρ
i
t
,
i
=
1,2
,
…
,where X is a matrix with proper dimension and ρ
i
(
t
) is an adjoint variable to be determined.Substitute the Lagrangian operatorλ
(
i
)
(
t
) into (36) and (38). We have(46)
ρ
˙
i
t
=
-
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
T
ρ
i
t
-
i
w
c
g
,
i
t
,
t
∈
0
,
T
-
τ
-
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
T
ρ
i
t
-
i
X
A
c
g
v
c
g
,
i
t
,
t
∈
T
-
τ
,
T
,where(47)
w
c
g
,
i
t
=
A
c
g
T
X
u
c
g
t
,
i
t
+
X
A
c
g
v
c
g
t
,
i
t
-
A
c
g
T
ρ
i
-
1
t
+
τ
-
ρ
i
-
1
t
,
v
c
g
,
i
t
=
θ
i
-
1
t
-
θ
i
-
1
t
-
τ
,
u
c
g
t
,
i
t
=
θ
i
-
1
t
+
τ
-
θ
i
-
1
t
.In addition, the state equation with the ith-order problem is(48)
θ
˙
i
t
=
A
c
g
-
B
c
g
R
c
g
B
c
g
T
X
T
θ
i
t
-
B
c
g
R
c
g
-
1
B
c
g
T
ρ
i
t
-
i
A
c
g
v
c
g
t
,
i
tin which the initial values are(49)
θ
i
0
=
0
,
i
=
1,2
,
…
.Thus, the control strategy of theith-order problem is(50)
U
i
t
=
-
R
c
g
-
1
B
c
g
T
X
θ
i
t
+
ρ
i
t
,
t
∈
0
,
T
,and the suboptimal control strategy is obtained by the former M items(51)
U
t
=
-
∑
i
=
1
M
1
i
!
R
c
g
-
1
B
c
g
T
X
θ
i
t
+
ρ
i
t
=
-
R
c
g
-
1
B
c
g
T
X
θ
t
+
∑
i
=
1
M
1
i
!
ρ
i
t
.In (51), there exists an optimal control strategy when M tends to ∞, there exists an optimal control strategy. In particular, the control strategy is U
(
t
)
=
-
R
c
g
-
1
B
c
g
T
X
θ
(
t
) when M
=
0.In the following, the optimal algorithm of the CGOC problem (22) is given with proper Maclaurin items M.Algorithm 1.
1
Step
1. Solve the Riccati equation (44) to get X.
Step
2. Get θ
(
0
)
(
t
) and λ
(
0
)
(
t
) by (41) and (39). Let i
=
1.
Step
3. Calculate ρ
i
(
t
) by (46).
Step
4. If i
<
M, solve θ
(
i
)
(
t
) by (48). Otherwise, execute Step 6.
Step
5. Let i
=
i
+
1, and then execute Step 3.
Step
6. Get the optimal control strategy U
∗
(
t
) by (51).
According to the cooperative game theory, the optimal control strategy is a function of the parameterα. Thus, a recursive algorithm can be given to get the Nash bargaining solution.Algorithm 2.
2
Step
1. Set the initial value α
0
=
[
1
/
n
,
…
,
1
/
n
].
Step
2. Calculate the suboptimal control strategy U
∗
(
α
0
) By Algorithm 1.Step
3. Check the equation J
i
c
(
U
∗
)
≤
J
t
p
,
i, ∀
i
=
1
,
…
,
n. If it does not hold, there exists an e satisfying J
e
c
(
U
∗
)
>
J
t
p
,
e. Update the parameter α
0 with the following formula, and then execute Step 2:(52)
α
e
0
=
α
e
0
+
0.001
,
α
i
0
=
α
i
0
-
0.001
n
-
1
,
i
=
1
,
…
,
n
,
i
≠
e
.Step
4. Compute α
~
j by (16).Step
5. Update parameter α by(53)
α
j
0
=
0.8
α
j
0
+
0.2
α
~
j
.If |
α
~
j
-
α
j
0
|
<
0.001, j
=
1
,
…
,
n hold, then terminate the algorithm and let α
=
α
~. Otherwise, execute Step 2.Next, we will discuss the convergence of the CGOC algorithm. WhenM
=
0, the CGOC algorithm can be described as(54)
θ
˙
t
=
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
t
+
A
c
g
θ
t
-
τ
q
.The items -
B
c
g
R
c
g
-
1
B
c
g
T
X and A
c
g can be seen as the coefficients I
n
⊗
H and L
σ
⊗
Γ in Theorem 2 (l
=
0) in [17]. Thus, the CGOC algorithm can converge with the given communication delay. When M
>
0, the θ
˙
(
t
) can also be rearranged with the items θ
(
t
) and θ
(
t
-
τ
q
). Similarly, the convergence range can be obtained.The distributed solving method of the rendezvous problem is given in Algorithm3.Algorithm 3.
2
Step
1. Get the waypoints sequence of each UAV using the route planning method [18]. Consider(55)
W
E
T
A
,
i
=
W
E
T
A
,
i
0
,
W
E
T
A
,
i
1
,
…
,
W
E
T
A
,
i
f
,
i
=
1
,
…
,
n
.
Step
2. Calculate the path length in the WGS84 coordinate system and get the ETA ranges which satisfy the restrict (3). Define the set R
E
T
A
,
i
=
{
d
E
T
A
,
i
/
v
m
a
x
,
d
E
T
A
,
i
/
v
m
i
n
}, i
=
1,2
,
…
,
n. If R
E
T
A
,
1
∩
R
E
T
A
,
2
∩
⋯
∩
R
E
T
A
,
N
≠
∅ holds, then execute Step 3. Otherwise, return to Step 1.
Step
3. Exchange the estimated arriving time θ
i
(
t
) with the other UAVs and confirm the communication delay τ
i
j
t, j
∈
N
i. Let τ
=
m
a
x
(
τ
i
j
(
t
)
).
Step
4. Get U
∗
(
α
) and α
∗ by Algorithms 1 and 2. Then, obtain the optimal control strategy U
∗
=
[
u
1
l
,
…
,
u
n
l
].
Step
5. Select the ith component u
i
l of U
∗, and substitute it into (11). The optimal reference velocity instruction v
i
c
(
t
) can be obtained(56)
v
i
c
t
=
v
i
t
+
v
i
t
·
1
/
2
H
v
,
i
-
1
∑
j
∈
N
i
Q
i
j
T
θ
j
t
-
τ
i
j
+
u
i
l
a
V
·
θ
i
t
.
Step
6. Calculate the reference heading instruction φ
i
c
(
t
) considering the waypoints information.
Step
7. Send the instructions v
i
c
(
t
) and φ
i
c
(
t
) to the autopilot system (4). If the UAVs arrive at target position, terminate the algorithm; otherwise, execute Step 3.


## 4.1. Cooperative Game Theory

Considering a team withn players, the quadratic function (12) is constructed to describe the cost of the player in the team (12)
J
i
=
1
2
∫
0
T
x
i
T
Q
i
x
i
+
u
i
T
R
i
u
i
d
t
,where Q
i and R
i are positive definite matrices with proper dimension. Each player satisfies the following dynamical model:(13)
x
˙
=
A
x
+
B
1
u
1
+
⋯
+
B
n
u
n
,
x
0
=
x
0
,in which A and B
i
(
i
=
1
,
…
,
n
) are constant matrices with proper dimensions and u
i is the decision of the ith player. The state variable x can be reflected by the other player’s decision in the minimization process. That means that the players may have conflicting interests. If a player decides to minimize its cost in a noncooperative manner, the decision chosen by the ith player can increase the cost of the others due to the coupling relationship. However, if the players decide to cooperate, individual cost may be minimized and hence we can get a smaller team cost. This will result in the set of Pareto-efficient solutions. For the set of inequalities J
i
U
≤
J
i
U
∗
,
i
=
1,2
,
…
,
n, if there is not at least one solution U, U
∗
=
[
u
1
,
u
2
,
…
,
u
n
] is called Pareto-efficient solution and the corresponding costs J
∗
=
[
J
1
(
U
∗
)
,
J
2
(
U
∗
)
,
…
,
J
n
(
U
∗
)
] are Pareto solution.The solution of the minimization problem cannot be dominated by any other solution(14)
U
∗
α
=
arg
⁡
min
U
∈
U
⁡
∑
i
=
1
n
α
i
J
i
U
,where α
=
(
α
1
,
…
,
α
n
) and ∑
i
=
1
n
α
i
=
1. It is a set of Pareto-efficient solutions of the above problem. The solutions are the functions of the parameter α. The final solution should be selected according to an axiomatic approach as our decision for the cooperative problem. The Nash-bargaining solution is selected in this paper(15)
α
=
arg
⁡
max
α
⁡
∏
i
=
1
n
J
i
α
,
U
∗
-
J
t
p
,
i
,where J
t
p
,
i
(
i
=
1
,
…
,
n
) are the individual costs calculated by using the noncooperation solution that is obtained by minimizing the cost (12).The coefficientα can be obtained (see Theorem  6.10 in [15]) by(16)
α
j
∗
=
∏
i
≠
j
J
i
∗
α
∗
,
U
∗
-
J
t
p
,
i
∑
i
=
1
n
∏
k
≠
i
J
k
∗
α
∗
,
U
∗
-
J
t
p
,
k
.


## 4.2. CGOC Problem Description

When the multiagent systems achieve consensus, all agents get the same values. Thus, the individual cost is defined as(17)
J
i
c
=
1
2
·
∫
0
T
∑
j
∈
N
i
θ
i
-
θ
j
T
Q
i
j
θ
i
-
θ
j
+
u
i
T
R
i
u
i
d
t
,where Q
i
j and R
i are symmetric positive definite matrices.The overall cost of the multiagent systems is obtained by weighted individual cost(18)
J
c
g
=
∑
i
=
1
n
α
i
J
i
c
U
=
1
2
∫
0
∞
θ
T
Q
c
g
θ
+
U
T
R
c
g
U
d
t
,where the coefficient matrices are(19)
Q
c
g
=
Q
i
j
n
×
n
,
Q
i
i
=
∑
j
∈
N
i
α
j
Q
j
i
+
∑
k
∈
N
i
α
i
Q
i
k
,
Q
i
j
=
-
α
i
Q
i
j
-
α
j
Q
j
i
j
∈
N
i
0
j
∉
N
i
,
R
c
g
=
diag
⁡
α
1
R
1
,
…
,
α
n
R
n
.According to cooperative game theory, the smallest overall costJ
c
g can be obtained when the multiagent system achieve consensus. Taking the communication delay (τ
1
=
⋯
=
τ
m
=
τ) into account, the dynamic model of the agent system can be described as follows:(20)
θ
˙
t
=
A
c
g
θ
t
-
τ
+
B
c
g
U
t
,in which θ
=
[
θ
1
T
,
…
,
θ
n
T
]
T is the state vector. A
c
g
=
[
A
i
j
]
n
×
n, A
i
i
=
0, A
i
j
=
-
(
1
/
2
)
l
i
j
H
v
,
i
-
1
Q
i
j, and B
c
g
=
I. H
v
,
i is the solution of the Riccati equation (21), and l
i
j is the element of the Laplacian matrix, which is used to describe the interconnection of the communication(21)
-
H
v
,
i
T
R
i
-
1
H
v
,
i
+
∑
j
∈
N
i
Q
i
j
=
0
.The cooperative game based optimal consensus (CGOC) can be described as the following optimization problem:(22)
min
⁡
U
→
U
J
c
g
s.t.
θ
˙
t
=
A
c
g
θ
t
-
τ
+
B
c
g
U
t
.According to cooperative game theory, the Pareto-optimal solution set can be obtained by(23)
U
∗
α
=
arg
⁡
min
U
∈
U
⁡
∑
i
=
1
n
α
i
J
i
c
U
.It is easy to see that Pareto-efficient solution is a function of the parameterα. Further, the unique Nash bargaining solution is calculated from the above set. Consider(24)
α
=
arg
⁡
max
α
⁡
∏
i
=
1
n
J
i
c
α
,
U
∗
-
J
t
p
,
i
,in which J
t
p
,
i is the individual cost calculated using the noncooperative strategy (see [10])(25)
u
t
p
,
i
=
1
2
H
v
,
i
-
1
∑
j
∈
N
i
Q
i
j
T
θ
j
t
-
τ
i
j
-
R
i
-
1
H
v
,
i
t
θ
i
.When the multiagent systems achieve consensus, the optimal parameterα
∗ is determined. In the following, we will discuss how to get the optimal control strategy of problem (22).


## 4.3. CGOC Problem Solving

Due to the existence of the communication delay in the dynamic model (20) of the agent, it is difficult to get the exact solution or a numerical solution of such a problem. The solving method of the two-point boundary problem with delay and ahead term is adopted in this section [16].Firstly, introduce the Lagrangian operatorλ
(
t
) and define the Hamilton function(26)
H
c
g
θ
,
λ
,
U
≜
1
2
θ
T
Q
c
g
θ
+
U
T
R
c
g
U
+
λ
T
t
A
c
g
θ
t
-
τ
+
B
c
g
U
t
.By the optimal control theory, we can get(27)
θ
˙
t
=
∇
λ
H
c
g
,
t
∈
0
,
T
,
λ
˙
t
=
-
∇
θ
H
c
g
-
∇
θ
τ
H
c
g
t
∈
0
,
T
-
τ
-
∇
θ
H
c
g
t
∈
T
-
τ
,
T
,where ∇
λ
H
c
g, ∇
θ
H
c
g, and ∇
θ
τ
H
c
g are the gradients at λ, θ, and θ
τ. By formula (26), formulas (27) can be rewritten as(28)
θ
˙
t
=
A
c
g
θ
t
-
τ
+
B
c
g
U
t
,
t
∈
0
,
T
,
λ
˙
t
=
-
Q
c
g
θ
t
-
A
c
g
T
λ
t
+
τ
t
∈
0
,
T
-
τ
-
Q
c
g
θ
t
t
∈
T
-
τ
,
T
,in which the initial value is θ
t
=
ϑ
t
,
t
∈
[
-
τ
,
0
]. Thus, the optimal control strategy is(29)
U
t
=
-
R
c
g
-
1
B
c
g
T
λ
t
,
t
∈
0
,
T
.Rewriteθ
˙
(
t
) and λ
˙
(
t
) as follows:(30)
θ
˙
t
=
A
c
g
θ
t
+
B
c
g
U
t
-
A
c
g
θ
t
-
θ
t
-
τ
,
t
∈
0
,
T
,
λ
˙
t
=
-
Q
c
g
θ
t
-
A
c
g
T
λ
t
-
A
c
g
T
λ
t
+
τ
-
λ
t
,
t
∈
0
,
T
-
τ
-
Q
c
g
θ
t
-
A
c
g
T
λ
t
,
t
∈
T
-
τ
,
T
.The sensitivity parameterɛ is introduced; we have(31)
θ
˙
t
,
ɛ
=
A
c
g
θ
t
,
ɛ
+
B
c
g
U
t
,
ɛ
-
ɛ
A
c
g
θ
t
,
ɛ
-
θ
t
-
τ
,
ɛ
,
t
∈
0
,
T
,
λ
˙
t
,
ɛ
=
-
Q
c
g
θ
t
,
ɛ
-
A
c
g
T
λ
t
,
ɛ
-
ɛ
A
c
g
T
λ
t
+
τ
,
ɛ
-
λ
t
,
ɛ
,
t
∈
0
,
T
-
τ
-
Q
c
g
θ
t
,
ɛ
-
A
c
g
T
λ
t
,
ɛ
,
t
∈
T
-
τ
,
T
,
U
t
,
ɛ
=
-
R
c
g
-
1
B
c
g
T
λ
t
,
ɛ
,
k
∈
0
,
T
,and, here, ɛ is a scalar and it satisfies 0
≤
ɛ
≤
1.Assume thatθ
t
,
ɛ, λ
(
t
,
ɛ
), and U
(
t
,
ɛ
) are differentiable at ɛ
=
0, and their Maclaurin series converges at ɛ
=
1
(32)
θ
t
,
ɛ
=
∑
i
=
0
∞
1
i
!
θ
i
t
,
λ
t
,
ɛ
=
∑
i
=
0
∞
1
i
!
λ
i
t
,
U
t
,
ɛ
=
∑
i
=
0
∞
1
i
!
U
i
t
,where θ
i
t
=
lim
ɛ
→
0
⁡
(
∂
i
θ
(
t
,
ɛ
)
/
∂
ɛ
i
), λ
i
(
t
)
=
l
i
m
ɛ
→
0
(
∂
i
λ
(
t
,
ɛ
)
/
∂
ɛ
i
), and U
i
(
t
)
=
l
i
m
ɛ
→
0
(
∂
i
U
(
t
,
ɛ
)
/
∂
ɛ
i
).Thus, the suboptimal control strategy of the above problem (31) is equivalent to the sum of the solutions of the two-point boundary problem, including the 0th-order and the ith-order problems.(i)
Two-point boundary problem with the 0th-order(33)
θ
˙
0
t
=
A
c
g
θ
0
t
+
B
c
g
U
0
t
,
(34)
U
0
t
=
-
R
c
g
-
1
B
c
g
T
λ
0
t
,
(35)
λ
˙
0
t
=
-
Q
c
g
θ
0
t
-
A
c
g
T
λ
0
t
.
(ii)
Two-point boundary problem with theith-order (36)
θ
˙
i
t
=
A
c
g
θ
i
t
+
B
c
g
U
i
t
+
i
A
c
g
θ
i
-
1
t
-
θ
i
-
1
t
-
τ
,
(37)
U
i
t
=
-
R
c
g
-
1
B
c
g
T
λ
i
t
,
(38)
λ
˙
i
t
=
-
Q
c
g
θ
i
t
-
A
c
g
T
λ
i
t
-
A
c
g
T
λ
i
-
1
t
+
τ
-
λ
i
-
1
t
,
t
∈
0
,
T
-
τ
-
Q
c
g
θ
i
t
-
A
c
g
T
λ
i
t
,
t
∈
T
-
τ
,
T
.


### 4.3.1. The Solution of the 0th-Order Problem

Assume that the Lagrangian operator is(39)
λ
0
t
≜
X
θ
0
t
,where X is a matrix with proper dimension.Substituting formula (34) with (39), we get(40)
U
0
t
=
-
R
c
g
-
1
B
c
g
T
X
θ
0
t
,
t
∈
0
,
T
,
(41)
θ
˙
0
t
=
A
c
g
θ
0
t
+
B
c
g
U
0
t
=
A
c
g
θ
0
t
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
0
t
=
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
0
t
,
t
∈
0
,
T
.Moreover, the Lagrangian operator should satisfy(42)
λ
˙
0
t
=
X
θ
˙
0
t
=
X
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
0
t
,
t
∈
0
,
T
.On the other hand, formula (35) can be rewritten as(43)
λ
˙
0
t
=
-
Q
c
g
θ
0
t
-
A
c
g
T
X
θ
0
t
=
-
Q
c
g
-
A
c
g
T
X
θ
0
t
.Comparing (42) with (43), the Riccati matrices equation can be obtained as(44)
A
c
g
T
X
+
X
A
c
g
-
X
B
c
g
R
c
g
-
1
B
c
g
T
X
+
Q
c
g
=
0
.Therefore, the optimal control strategyU
(
0
) of the 0th-order problem can be obtained through calculating X and θ
(
0
).


### 4.3.2.  The Solution of theith-Order Problem

Similar to the0th-order problem, the Lagrangian operator is selected as(45)
λ
i
t
≜
X
θ
i
t
+
ρ
i
t
,
i
=
1,2
,
…
,where X is a matrix with proper dimension and ρ
i
(
t
) is an adjoint variable to be determined.Substitute the Lagrangian operatorλ
(
i
)
(
t
) into (36) and (38). We have(46)
ρ
˙
i
t
=
-
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
T
ρ
i
t
-
i
w
c
g
,
i
t
,
t
∈
0
,
T
-
τ
-
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
T
ρ
i
t
-
i
X
A
c
g
v
c
g
,
i
t
,
t
∈
T
-
τ
,
T
,where(47)
w
c
g
,
i
t
=
A
c
g
T
X
u
c
g
t
,
i
t
+
X
A
c
g
v
c
g
t
,
i
t
-
A
c
g
T
ρ
i
-
1
t
+
τ
-
ρ
i
-
1
t
,
v
c
g
,
i
t
=
θ
i
-
1
t
-
θ
i
-
1
t
-
τ
,
u
c
g
t
,
i
t
=
θ
i
-
1
t
+
τ
-
θ
i
-
1
t
.In addition, the state equation with the ith-order problem is(48)
θ
˙
i
t
=
A
c
g
-
B
c
g
R
c
g
B
c
g
T
X
T
θ
i
t
-
B
c
g
R
c
g
-
1
B
c
g
T
ρ
i
t
-
i
A
c
g
v
c
g
t
,
i
tin which the initial values are(49)
θ
i
0
=
0
,
i
=
1,2
,
…
.Thus, the control strategy of theith-order problem is(50)
U
i
t
=
-
R
c
g
-
1
B
c
g
T
X
θ
i
t
+
ρ
i
t
,
t
∈
0
,
T
,and the suboptimal control strategy is obtained by the former M items(51)
U
t
=
-
∑
i
=
1
M
1
i
!
R
c
g
-
1
B
c
g
T
X
θ
i
t
+
ρ
i
t
=
-
R
c
g
-
1
B
c
g
T
X
θ
t
+
∑
i
=
1
M
1
i
!
ρ
i
t
.In (51), there exists an optimal control strategy when M tends to ∞, there exists an optimal control strategy. In particular, the control strategy is U
(
t
)
=
-
R
c
g
-
1
B
c
g
T
X
θ
(
t
) when M
=
0.In the following, the optimal algorithm of the CGOC problem (22) is given with proper Maclaurin items M.Algorithm 1.
1
Step
1. Solve the Riccati equation (44) to get X.
Step
2. Get θ
(
0
)
(
t
) and λ
(
0
)
(
t
) by (41) and (39). Let i
=
1.
Step
3. Calculate ρ
i
(
t
) by (46).
Step
4. If i
<
M, solve θ
(
i
)
(
t
) by (48). Otherwise, execute Step 6.
Step
5. Let i
=
i
+
1, and then execute Step 3.
Step
6. Get the optimal control strategy U
∗
(
t
) by (51).
According to the cooperative game theory, the optimal control strategy is a function of the parameterα. Thus, a recursive algorithm can be given to get the Nash bargaining solution.Algorithm 2.
2
Step
1. Set the initial value α
0
=
[
1
/
n
,
…
,
1
/
n
].
Step
2. Calculate the suboptimal control strategy U
∗
(
α
0
) By Algorithm 1.Step
3. Check the equation J
i
c
(
U
∗
)
≤
J
t
p
,
i, ∀
i
=
1
,
…
,
n. If it does not hold, there exists an e satisfying J
e
c
(
U
∗
)
>
J
t
p
,
e. Update the parameter α
0 with the following formula, and then execute Step 2:(52)
α
e
0
=
α
e
0
+
0.001
,
α
i
0
=
α
i
0
-
0.001
n
-
1
,
i
=
1
,
…
,
n
,
i
≠
e
.Step
4. Compute α
~
j by (16).Step
5. Update parameter α by(53)
α
j
0
=
0.8
α
j
0
+
0.2
α
~
j
.If |
α
~
j
-
α
j
0
|
<
0.001, j
=
1
,
…
,
n hold, then terminate the algorithm and let α
=
α
~. Otherwise, execute Step 2.Next, we will discuss the convergence of the CGOC algorithm. WhenM
=
0, the CGOC algorithm can be described as(54)
θ
˙
t
=
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
t
+
A
c
g
θ
t
-
τ
q
.The items -
B
c
g
R
c
g
-
1
B
c
g
T
X and A
c
g can be seen as the coefficients I
n
⊗
H and L
σ
⊗
Γ in Theorem 2 (l
=
0) in [17]. Thus, the CGOC algorithm can converge with the given communication delay. When M
>
0, the θ
˙
(
t
) can also be rearranged with the items θ
(
t
) and θ
(
t
-
τ
q
). Similarly, the convergence range can be obtained.The distributed solving method of the rendezvous problem is given in Algorithm3.Algorithm 3.
2
Step
1. Get the waypoints sequence of each UAV using the route planning method [18]. Consider(55)
W
E
T
A
,
i
=
W
E
T
A
,
i
0
,
W
E
T
A
,
i
1
,
…
,
W
E
T
A
,
i
f
,
i
=
1
,
…
,
n
.
Step
2. Calculate the path length in the WGS84 coordinate system and get the ETA ranges which satisfy the restrict (3). Define the set R
E
T
A
,
i
=
{
d
E
T
A
,
i
/
v
m
a
x
,
d
E
T
A
,
i
/
v
m
i
n
}, i
=
1,2
,
…
,
n. If R
E
T
A
,
1
∩
R
E
T
A
,
2
∩
⋯
∩
R
E
T
A
,
N
≠
∅ holds, then execute Step 3. Otherwise, return to Step 1.
Step
3. Exchange the estimated arriving time θ
i
(
t
) with the other UAVs and confirm the communication delay τ
i
j
t, j
∈
N
i. Let τ
=
m
a
x
(
τ
i
j
(
t
)
).
Step
4. Get U
∗
(
α
) and α
∗ by Algorithms 1 and 2. Then, obtain the optimal control strategy U
∗
=
[
u
1
l
,
…
,
u
n
l
].
Step
5. Select the ith component u
i
l of U
∗, and substitute it into (11). The optimal reference velocity instruction v
i
c
(
t
) can be obtained(56)
v
i
c
t
=
v
i
t
+
v
i
t
·
1
/
2
H
v
,
i
-
1
∑
j
∈
N
i
Q
i
j
T
θ
j
t
-
τ
i
j
+
u
i
l
a
V
·
θ
i
t
.
Step
6. Calculate the reference heading instruction φ
i
c
(
t
) considering the waypoints information.
Step
7. Send the instructions v
i
c
(
t
) and φ
i
c
(
t
) to the autopilot system (4). If the UAVs arrive at target position, terminate the algorithm; otherwise, execute Step 3.


## 4.3.1. The Solution of the 0th-Order Problem

Assume that the Lagrangian operator is(39)
λ
0
t
≜
X
θ
0
t
,where X is a matrix with proper dimension.Substituting formula (34) with (39), we get(40)
U
0
t
=
-
R
c
g
-
1
B
c
g
T
X
θ
0
t
,
t
∈
0
,
T
,
(41)
θ
˙
0
t
=
A
c
g
θ
0
t
+
B
c
g
U
0
t
=
A
c
g
θ
0
t
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
0
t
=
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
0
t
,
t
∈
0
,
T
.Moreover, the Lagrangian operator should satisfy(42)
λ
˙
0
t
=
X
θ
˙
0
t
=
X
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
0
t
,
t
∈
0
,
T
.On the other hand, formula (35) can be rewritten as(43)
λ
˙
0
t
=
-
Q
c
g
θ
0
t
-
A
c
g
T
X
θ
0
t
=
-
Q
c
g
-
A
c
g
T
X
θ
0
t
.Comparing (42) with (43), the Riccati matrices equation can be obtained as(44)
A
c
g
T
X
+
X
A
c
g
-
X
B
c
g
R
c
g
-
1
B
c
g
T
X
+
Q
c
g
=
0
.Therefore, the optimal control strategyU
(
0
) of the 0th-order problem can be obtained through calculating X and θ
(
0
).


## 4.3.2.  The Solution of theith-Order Problem

Similar to the0th-order problem, the Lagrangian operator is selected as(45)
λ
i
t
≜
X
θ
i
t
+
ρ
i
t
,
i
=
1,2
,
…
,where X is a matrix with proper dimension and ρ
i
(
t
) is an adjoint variable to be determined.Substitute the Lagrangian operatorλ
(
i
)
(
t
) into (36) and (38). We have(46)
ρ
˙
i
t
=
-
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
T
ρ
i
t
-
i
w
c
g
,
i
t
,
t
∈
0
,
T
-
τ
-
A
c
g
-
B
c
g
R
c
g
-
1
B
c
g
T
X
T
ρ
i
t
-
i
X
A
c
g
v
c
g
,
i
t
,
t
∈
T
-
τ
,
T
,where(47)
w
c
g
,
i
t
=
A
c
g
T
X
u
c
g
t
,
i
t
+
X
A
c
g
v
c
g
t
,
i
t
-
A
c
g
T
ρ
i
-
1
t
+
τ
-
ρ
i
-
1
t
,
v
c
g
,
i
t
=
θ
i
-
1
t
-
θ
i
-
1
t
-
τ
,
u
c
g
t
,
i
t
=
θ
i
-
1
t
+
τ
-
θ
i
-
1
t
.In addition, the state equation with the ith-order problem is(48)
θ
˙
i
t
=
A
c
g
-
B
c
g
R
c
g
B
c
g
T
X
T
θ
i
t
-
B
c
g
R
c
g
-
1
B
c
g
T
ρ
i
t
-
i
A
c
g
v
c
g
t
,
i
tin which the initial values are(49)
θ
i
0
=
0
,
i
=
1,2
,
…
.Thus, the control strategy of theith-order problem is(50)
U
i
t
=
-
R
c
g
-
1
B
c
g
T
X
θ
i
t
+
ρ
i
t
,
t
∈
0
,
T
,and the suboptimal control strategy is obtained by the former M items(51)
U
t
=
-
∑
i
=
1
M
1
i
!
R
c
g
-
1
B
c
g
T
X
θ
i
t
+
ρ
i
t
=
-
R
c
g
-
1
B
c
g
T
X
θ
t
+
∑
i
=
1
M
1
i
!
ρ
i
t
.In (51), there exists an optimal control strategy when M tends to ∞, there exists an optimal control strategy. In particular, the control strategy is U
(
t
)
=
-
R
c
g
-
1
B
c
g
T
X
θ
(
t
) when M
=
0.In the following, the optimal algorithm of the CGOC problem (22) is given with proper Maclaurin items M.Algorithm 1.
1
Step
1. Solve the Riccati equation (44) to get X.
Step
2. Get θ
(
0
)
(
t
) and λ
(
0
)
(
t
) by (41) and (39). Let i
=
1.
Step
3. Calculate ρ
i
(
t
) by (46).
Step
4. If i
<
M, solve θ
(
i
)
(
t
) by (48). Otherwise, execute Step 6.
Step
5. Let i
=
i
+
1, and then execute Step 3.
Step
6. Get the optimal control strategy U
∗
(
t
) by (51).
According to the cooperative game theory, the optimal control strategy is a function of the parameterα. Thus, a recursive algorithm can be given to get the Nash bargaining solution.Algorithm 2.
2
Step
1. Set the initial value α
0
=
[
1
/
n
,
…
,
1
/
n
].
Step
2. Calculate the suboptimal control strategy U
∗
(
α
0
) By Algorithm 1.Step
3. Check the equation J
i
c
(
U
∗
)
≤
J
t
p
,
i, ∀
i
=
1
,
…
,
n. If it does not hold, there exists an e satisfying J
e
c
(
U
∗
)
>
J
t
p
,
e. Update the parameter α
0 with the following formula, and then execute Step 2:(52)
α
e
0
=
α
e
0
+
0.001
,
α
i
0
=
α
i
0
-
0.001
n
-
1
,
i
=
1
,
…
,
n
,
i
≠
e
.Step
4. Compute α
~
j by (16).Step
5. Update parameter α by(53)
α
j
0
=
0.8
α
j
0
+
0.2
α
~
j
.If |
α
~
j
-
α
j
0
|
<
0.001, j
=
1
,
…
,
n hold, then terminate the algorithm and let α
=
α
~. Otherwise, execute Step 2.Next, we will discuss the convergence of the CGOC algorithm. WhenM
=
0, the CGOC algorithm can be described as(54)
θ
˙
t
=
-
B
c
g
R
c
g
-
1
B
c
g
T
X
θ
t
+
A
c
g
θ
t
-
τ
q
.The items -
B
c
g
R
c
g
-
1
B
c
g
T
X and A
c
g can be seen as the coefficients I
n
⊗
H and L
σ
⊗
Γ in Theorem 2 (l
=
0) in [17]. Thus, the CGOC algorithm can converge with the given communication delay. When M
>
0, the θ
˙
(
t
) can also be rearranged with the items θ
(
t
) and θ
(
t
-
τ
q
). Similarly, the convergence range can be obtained.The distributed solving method of the rendezvous problem is given in Algorithm3.Algorithm 3.
2
Step
1. Get the waypoints sequence of each UAV using the route planning method [18]. Consider(55)
W
E
T
A
,
i
=
W
E
T
A
,
i
0
,
W
E
T
A
,
i
1
,
…
,
W
E
T
A
,
i
f
,
i
=
1
,
…
,
n
.
Step
2. Calculate the path length in the WGS84 coordinate system and get the ETA ranges which satisfy the restrict (3). Define the set R
E
T
A
,
i
=
{
d
E
T
A
,
i
/
v
m
a
x
,
d
E
T
A
,
i
/
v
m
i
n
}, i
=
1,2
,
…
,
n. If R
E
T
A
,
1
∩
R
E
T
A
,
2
∩
⋯
∩
R
E
T
A
,
N
≠
∅ holds, then execute Step 3. Otherwise, return to Step 1.
Step
3. Exchange the estimated arriving time θ
i
(
t
) with the other UAVs and confirm the communication delay τ
i
j
t, j
∈
N
i. Let τ
=
m
a
x
(
τ
i
j
(
t
)
).
Step
4. Get U
∗
(
α
) and α
∗ by Algorithms 1 and 2. Then, obtain the optimal control strategy U
∗
=
[
u
1
l
,
…
,
u
n
l
].
Step
5. Select the ith component u
i
l of U
∗, and substitute it into (11). The optimal reference velocity instruction v
i
c
(
t
) can be obtained(56)
v
i
c
t
=
v
i
t
+
v
i
t
·
1
/
2
H
v
,
i
-
1
∑
j
∈
N
i
Q
i
j
T
θ
j
t
-
τ
i
j
+
u
i
l
a
V
·
θ
i
t
.
Step
6. Calculate the reference heading instruction φ
i
c
(
t
) considering the waypoints information.
Step
7. Send the instructions v
i
c
(
t
) and φ
i
c
(
t
) to the autopilot system (4). If the UAVs arrive at target position, terminate the algorithm; otherwise, execute Step 3.


## 5. Simulations and Results

### 5.1. Simulation Environment and Experimental Conditions

In order to validate the effectiveness of the proposed method, two experiments are illustrated in this section. A typical simulation environment of the multi-UAV rendezvous problem is given in Figure3. Six artillery threats are placed in the square area (100 km × 100 km). Let the initial positions of three UAV be (20 km, 20 km), (40 km, 10 km), and (60 km, 15 km), and let the target positions be (50 km, 70 km), (60 km, 60 km), and (70 km, 70 km). The results of the route planning are shown in the figure. The initial velocities of the three vehicles are 27 m/s 32 m/s and 30 m/s, respectively. The initial path lengths of the three UAVs are 64.96 km, 55.40 km, and 61.67 km using the TS module in the distributed solving structure. Their performance parameters are listed as follows: R
m
i
n
=
20
m, a
v
=
0.2, a
φ
=
0.05, w
h
e
a
d
=
1.2
r
a
d, v
m
i
n
=
10
m
/
s, and v
m
a
x
=
40
m
/
s.Figure 3
The typical scenery of multi-UAVs rendezvous problem.Two network restrictions are considered in the experiments. One is the time varying communication delayτ
(
t
); the other one is the time varying interconnection of the communication network:(1)
time varying delayτ
(
t
):
(i)
delay I: Markov random process{
0,0.092,0.2
} with the probabilities 0.3, 0.3, and 0.4;
(ii)
delay II: normal distributionτ
(
t
)
~
N
(
0.15,0.0
1
2
);
(iii)
delay III: sine functionτ
(
t
)
=
0.2
|
s
i
n
(
0.9
/
0.2
t
)
|;
(2)
jointly connected topologies{
G
a
,
G
b
,
G
c
}; see Figure 4.Figure 4
The time varying interconnection of the communication network (the interval time of signalσ
(
t
)
:
[
0
,
i
n
f
]
→
{
a
,
b
,
c
} is 1 s).


### 5.2.  Experiment Results

In the simulation, we show the response of variableθ
E
T
A and velocity V. Let M
=
0. We can get the Nash bargaining solution α
=
[
0.423,0.301,0.276
] by Algorithms 1 and 2. Further, the variable θ
E
T
A and the velocity V
c can be obtained by Algorithm 3. The results are given in Figures 5, 6, and 7.Figure 5
The response cure ofθ
E
T
A and V under delay I.
(a)
Delay I
(b)
Coordination variable
(c)
Velocity variableFigure 6
The response cure ofθ
E
T
A and V under delay II.
(a)
Delay II
(b)
Coordination variable
(c)
Velocity variableFigure 7
The response cure ofθ
E
T
A and V under delay III.
(a)
Delay III
(b)
Coordination variable
(c)
Velocity variableIt is easy to see that the coordination variableθ
E
T
A achieves consensus from the Figures 5, 6, and 7. The time used to adjust the θ
E
T
A is less than 5 s. However, there is an obvious oscillation in that process. For the case of communication delays I, II, and III, the rendezvous times are 2103 s, 2106 s, and 2102 s, respectively. The velocity instruction V
c can be computed by Algorithm 3, and their response cures are also shown in Figures 5(c), 6(c), and 7(c).Furthermore, the simulation process of the rendezvous problem is shown in Figure8. The result illustrates the effectiveness of the CGOC method.Figure 8
The simulation of the rendezvous process based on CGOC method.
(a)
242 s
(b)
576 s
(c)
913 s
(d)
1264 s
(e)
1790 s
(f)
2102 s


### 5.3. The Comparison of the Optimality and Dynamic Response

The comparison of the energy consumptions is given by using the proposed CGOC method and the NCOC method [10]. Considering the randomness and uncertainty property of the delay, 20 statistical experiments are given in Table 1 with communication delays I, II, and III.Table 1
The comparisons of the energy consumptions using NCOC method and CGOC method.
Communication restrictions
NCOC method (×1
0
3)
CGOC method (×1
0
3)
UAV1
UAV2
UAV3
UAV1
UAV2
UAV3
Delay I and topologiesG
a
,
G
b
,
G
c
1650.6
1565.0
442.2
1620.8
1461.7
258.7
Delay II and topologiesG
a
,
G
b
,
G
c
1650.6
1565.4
442.5
1611.7
1453.2
257.8
Delay III and topologiesG
a
,
G
b
,
G
c
1651.3
1566.4
443.6
1597.0
1440.3
257.0The following conclusions can be obtained from Table1.(1)
The energy consumptions of the proposed CGOC method and the NCOC method are close. The results imply that the CGOC method has better robustness adaptability than NCOC method.
(2)
The overall cost using the CGOC method is less than that using the NCOC method. The case of number 3 UAV is the biggest one, about 41.7%. The second one is number 2 UAV, about 7.2%. Number 1 UAV is the minimum one, about 2.5%. In other words, the reduction of the overall cost is about 17.1%. Thus, the proposed CGOC method is better than NCOC in terms of the optimization.Another important issue is the dynamic response when new incidents occur. For the case of the threat appearing (the length decrease 1.5 km) at 40 s and disappearing (the length increase 1.5 km) at 80 s, the dynamic responses are given in Figures9 and 10. It is easy to see that CGOC requires less adjustment time and lower oscillation.Figure 9
The response cure ofθ
E
T
A and V using NCOC method.
(a)
Coordination variable
(b)
Velocity variableFigure 10
The response cure ofθ
E
T
A and V using CGOC method.
(a)
Coordination variable
(b)
Velocity variable


## 5.1. Simulation Environment and Experimental Conditions

In order to validate the effectiveness of the proposed method, two experiments are illustrated in this section. A typical simulation environment of the multi-UAV rendezvous problem is given in Figure3. Six artillery threats are placed in the square area (100 km × 100 km). Let the initial positions of three UAV be (20 km, 20 km), (40 km, 10 km), and (60 km, 15 km), and let the target positions be (50 km, 70 km), (60 km, 60 km), and (70 km, 70 km). The results of the route planning are shown in the figure. The initial velocities of the three vehicles are 27 m/s 32 m/s and 30 m/s, respectively. The initial path lengths of the three UAVs are 64.96 km, 55.40 km, and 61.67 km using the TS module in the distributed solving structure. Their performance parameters are listed as follows: R
m
i
n
=
20
m, a
v
=
0.2, a
φ
=
0.05, w
h
e
a
d
=
1.2
r
a
d, v
m
i
n
=
10
m
/
s, and v
m
a
x
=
40
m
/
s.Figure 3
The typical scenery of multi-UAVs rendezvous problem.Two network restrictions are considered in the experiments. One is the time varying communication delayτ
(
t
); the other one is the time varying interconnection of the communication network:(1)
time varying delayτ
(
t
):
(i)
delay I: Markov random process{
0,0.092,0.2
} with the probabilities 0.3, 0.3, and 0.4;
(ii)
delay II: normal distributionτ
(
t
)
~
N
(
0.15,0.0
1
2
);
(iii)
delay III: sine functionτ
(
t
)
=
0.2
|
s
i
n
(
0.9
/
0.2
t
)
|;
(2)
jointly connected topologies{
G
a
,
G
b
,
G
c
}; see Figure 4.Figure 4
The time varying interconnection of the communication network (the interval time of signalσ
(
t
)
:
[
0
,
i
n
f
]
→
{
a
,
b
,
c
} is 1 s).


## 5.2.  Experiment Results

In the simulation, we show the response of variableθ
E
T
A and velocity V. Let M
=
0. We can get the Nash bargaining solution α
=
[
0.423,0.301,0.276
] by Algorithms 1 and 2. Further, the variable θ
E
T
A and the velocity V
c can be obtained by Algorithm 3. The results are given in Figures 5, 6, and 7.Figure 5
The response cure ofθ
E
T
A and V under delay I.
(a)
Delay I
(b)
Coordination variable
(c)
Velocity variableFigure 6
The response cure ofθ
E
T
A and V under delay II.
(a)
Delay II
(b)
Coordination variable
(c)
Velocity variableFigure 7
The response cure ofθ
E
T
A and V under delay III.
(a)
Delay III
(b)
Coordination variable
(c)
Velocity variableIt is easy to see that the coordination variableθ
E
T
A achieves consensus from the Figures 5, 6, and 7. The time used to adjust the θ
E
T
A is less than 5 s. However, there is an obvious oscillation in that process. For the case of communication delays I, II, and III, the rendezvous times are 2103 s, 2106 s, and 2102 s, respectively. The velocity instruction V
c can be computed by Algorithm 3, and their response cures are also shown in Figures 5(c), 6(c), and 7(c).Furthermore, the simulation process of the rendezvous problem is shown in Figure8. The result illustrates the effectiveness of the CGOC method.Figure 8
The simulation of the rendezvous process based on CGOC method.
(a)
242 s
(b)
576 s
(c)
913 s
(d)
1264 s
(e)
1790 s
(f)
2102 s


## 5.3. The Comparison of the Optimality and Dynamic Response

The comparison of the energy consumptions is given by using the proposed CGOC method and the NCOC method [10]. Considering the randomness and uncertainty property of the delay, 20 statistical experiments are given in Table 1 with communication delays I, II, and III.Table 1
The comparisons of the energy consumptions using NCOC method and CGOC method.
Communication restrictions
NCOC method (×1
0
3)
CGOC method (×1
0
3)
UAV1
UAV2
UAV3
UAV1
UAV2
UAV3
Delay I and topologiesG
a
,
G
b
,
G
c
1650.6
1565.0
442.2
1620.8
1461.7
258.7
Delay II and topologiesG
a
,
G
b
,
G
c
1650.6
1565.4
442.5
1611.7
1453.2
257.8
Delay III and topologiesG
a
,
G
b
,
G
c
1651.3
1566.4
443.6
1597.0
1440.3
257.0The following conclusions can be obtained from Table1.(1)
The energy consumptions of the proposed CGOC method and the NCOC method are close. The results imply that the CGOC method has better robustness adaptability than NCOC method.
(2)
The overall cost using the CGOC method is less than that using the NCOC method. The case of number 3 UAV is the biggest one, about 41.7%. The second one is number 2 UAV, about 7.2%. Number 1 UAV is the minimum one, about 2.5%. In other words, the reduction of the overall cost is about 17.1%. Thus, the proposed CGOC method is better than NCOC in terms of the optimization.Another important issue is the dynamic response when new incidents occur. For the case of the threat appearing (the length decrease 1.5 km) at 40 s and disappearing (the length increase 1.5 km) at 80 s, the dynamic responses are given in Figures9 and 10. It is easy to see that CGOC requires less adjustment time and lower oscillation.Figure 9
The response cure ofθ
E
T
A and V using NCOC method.
(a)
Coordination variable
(b)
Velocity variableFigure 10
The response cure ofθ
E
T
A and V using CGOC method.
(a)
Coordination variable
(b)
Velocity variable


## 6. Conclusions

The cooperative game based optimal consensus (CGOC) algorithm has been proposed to solve the multi-UAV rendezvous problem with complex networks. Mainly, the following contributions have been concluded in this paper: (1) the mathematical model and distributed solving framework of the rendezvous problem have been established; (2) CGOC algorithm has been presented to minimize the overall cost of the multi-UAV system. The solving strategy of CGOC has been given theoretically using the cooperative game and sensitivity parameter method. Numerical examples and simulation results have been given to demonstrate the effectiveness with different network conditions and the benefit on the overall optimality and dynamic response.



---

*Source: 878536-2015-06-09.xml*