# Cooperative Search by Multiple Unmanned Aerial Vehicles in a Nonconvex Environment

Authors: Xiaoting Ji; Xiangke Wang; Yifeng Niu; Lincheng Shen
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper presents a distributed cooperative search algorithm for multiple unmanned aerial vehicles (UAVs) with limited sensing and communication capabilities in a nonconvex environment. The objective is to control multiple UAVs to find several unknown targets deployed in a given region, while minimizing the expected search time and avoiding obstacles. First, an asynchronous distributed cooperative search framework is proposed by integrating the information update into the coverage control scheme. And an adaptive density function is designed based on the real-time updated probability map and uncertainty map, which can balance target detection and environment exploration. Second, in order to handle nonconvex environment with arbitrary obstacles, a new transformation method is proposed to transform the cooperative search problem in the nonconvex region into an equivalent one in the convex region. Furthermore, a control strategy for cooperative search is proposed to plan feasible trajectories for UAVs under the kinematic constraints, and the convergence is proved by LaSalle’s invariance principle. Finally, by simulation results, it can be seen that our proposed algorithm is effective to handle the search problem in the nonconvex environment and efficient to find targets in shorter time compared with other algorithms.


---

## Body

## 1. Introduction

Over the past decade, unmanned air vehicles (UAVs) with functional diversity and low cost have been extensively employed in many civil and military applications, such as environment surveillance, battle reconnaissance, and search and rescue in the hazardous environment [1–3]. With the development of advanced sensing and information processing technology, cooperative search has been one of the most popular utilizations of UAVs equipped with sensors (such as camera, Lidar, and sonar) [4]. The goal of cooperative search is to control multiple UAVs to find several unknown targets deployed in a given region, while maximizing the detection probability and minimizing the expected search time.The problem of cooperative search with multiple UAVs has been studied extensively due to its critical importance for a myriad of applications [3, 5, 6]. Existing methods can be classified into two categories: predefined flight paths based search and dynamic path planning based search. The former is to generate the flight paths (e.g., parallel lines or outward spirals) in advance and follow these paths during search execution. The typical method of this category is sweep-line based search [7], in which the agents sweep all the points in the given area to find the targets. This method is effective so that no search areas are missed but not efficient due to the predefined paths and cannot be used for searching the dynamic targets. The latter method is to convert the cooperative search problem into a multiagent path planning issue. Typically, the path planning problem is formulated as the optimization of a team objective function subject to a set of constraints [8]. Therefore, dynamic programming [9], artificial intelligence [5, 6], and model predictive control (MPC) [10] can be used for solving these problems. In this paper, inspired by the coverage control, a distributed cooperative search method is proposed by integrating information update into the coverage control scheme, which also belongs to path planning based search methods. In recent years, with the development of related theory in Mobile Wireless Sensor Network (MWSN), coverage control methods [11] have been widely used in multiple robots system. The purpose of coverage control is to optimally deploy the mobile sensors to maximize the coverage of environment. Due to the limited Field of View (FoV), the cooperative search problem with multiple UAVs can be treated as an optimal coverage problem with a bounded sensing size of sensors [12]. And the distributed control strategy can be easily obtained by optimizing the objective function through the Lloyd algorithm [11]. Besides, by using the coverage control scheme in the search problem, it will not only promote the detection performance of target search but also take account of coverage performance. However, few of the existing coverage control methods have considered the information update about target existence probability, which is crucial to the search problem in practice and essentially affects the movements of UAVs.Therefore, the first issue that should be addressed is how to integrate the information update into the coverage control scheme in order to solve the search problem. To our knowledge, there are only a few works that utilize the coverage control method with consideration of information update. Zhong and Cassandras [13] use coverage control and data collection to maximize the joint detection probability of random events in a given mission space. However, their goal is to find the optimal deployments for detection, which is essentially different from the search problem. And the density function (importance level of each cell in the region) in their objective formulation depends on the distance to the critical points (i.e., the detected sources). Mirzaei et al. [14] use two different types of vehicles for search and coverage tasks. Search vehicles use a limited look-ahead dynamic programming algorithm to maximize the amount of information gathered by the whole team, while service vehicles use the coverage control method to spread out over the environment to optimally cover the terrain. However, the search and coverage control are executed on different platforms by different schemes, and only the probability distribution of critical points is utilized in the objective function. Besides, the information fusion update is not considered. Hu et al. [15] formulate the path planning as a coverage control problem to find an optimal configuration of all agents that minimizes a given coverage performance cost function. They use the uncertainty information as the density function to cover the region uniformly, which cannot directly facilitate agents to search the region with high target existence. In this paper, we integrate the information update into the coverage control scheme by introducing the probability map of target existence together with the uncertainty map into the density function of coverage optimization. In most coverage control problems, the density function is fixed [16, 17] or only depends on the distance to some critical points [18, 19]. However, our new density function can be adaptively changed, depending on the real-time updated target existence probability and uncertainty map. The higher the probability map is, the larger the probability of target existence is. So utilizing the probability map as the density function can help UAV converge to the regions with high target existence and find targets in a short time. If all the regions with high target existence probability have been searched, or the UAV gets stuck in the local optimum, the uncertainty map can help the UAV escape from the local optimum and explore more regions.Due to practical requirements, realistic search areas may be arbitrarily shaped with arbitrary obstacles. Thus, the other issue in our search problem is how to handle the nonconvex environment during search execution. Pimenta et al. [12] use the geodesic distance based algorithm to solve the coverage problem with a network of heterogeneous robots in the nonconvex environment. However, it may cause the robot to get stuck in a saddle point, or even worse, to drive into an obstacle, as the projections of the gradients into target direction of the geodesic distance do not add up to zero when reaching the boundary. Teraoka et al. [19] and Adibi et al. [20] present the potential field-based approaches representing control under the influence of virtual force generated by goals (attraction) and obstacles (repulsion), which have advantages in less computation and disadvantages in being easily trapped into a local minimum of the potential field. Breitenmoser et al. [18] present an algorithm which combines classical Voronoi coverage with Lloyd algorithm and the local planning algorithm TangentBug to compute the motion of the robots around obstacles and corners. It assumes that the range of sensor is infinite, meaning that the sensors can “see” through obstacles, which is unrealistic in practice. Caicedo et al. [16, 17] map a class of connected regions with holes to an almost convex region through a diffeomorphism, such that Voronoi partition and Lloyd algorithm [11] can be used to solve the coverage problem. However, the obstacles in their work must be simple convex ones, and the kinematic model of agents is not considered. As mentioned above, cooperative search in the nonconvex environment with arbitrary obstacles is still an open problem, especially when taking more realistic detection and the kinematic model into consideration. Motivated by this, we propose a generalized method to construct a transformation that can transform the nonconvex region with arbitrary obstacles into an almost convex region. Combining with the asynchronous distributed cooperative search framework, the cooperative search problem in the complex nonconvex environment can be solved.The main contributions of this work are as follows. First, an asynchronous distributed cooperative search framework is developed that integrates information update into coverage control scheme. And an adaptive density function is formulated depending on real-time updated probability map and uncertainty map, which can balance target search and environment exploration. Second, by extending the diffeomorphism [16], we propose a new transformation method to handle the nonconvex environment with arbitrary obstacles, not limited to convex ones. Based on the transformation, the cooperative search algorithm is provided that can address search problem in the nonconvex environment with arbitrary obstacles. Finally, a control strategy is designed considering the kinematic constraints of UAV, and the convergence is proved using LaSalle’s invariance principle.The remainder of this paper is organized as follows. In Section2, some useful preliminaries are provided. The cooperative search framework and the objective formulation, as well as the explicit information update procedure, are presented in Section 3. In Section 4, the generalized method for constructing a transformation is presented to handle the nonconvex environment with arbitrary obstacles. Then the cooperative search algorithm in the nonconvex environment is proposed, and the stability and convergence are proved. Simulation results with analysis are shown in Section 5, and the conclusions are drawn in Section 6.


## 2. Preliminaries

### 2.1. Basic Notions and Definitions

LetR be the set of real numbers and let R
n denote the n-dimensional Euclidean space. For a given set Q
⊆
R
n, Q
i
n
t and ∂
Q denote the interior and the boundary of Q, respectively. One has Q
i
n
t
=
Q
/
∂
Q. d
(
x
,
y
)
=
x
-
y represents the Euclidean distance between points x
,
y
∈
R
2, and B
(
q
,
ε
) is a closed ball centered at q with radius ε. In the following, some definitions about the environment are presented.Definition 1 (connected partition [21]).
LetQ
⊂
R
2 be a finite set and let P
=
P
1
,
…
,
P
N be a subset of Q; if P
i
∩
P
j
=
∅, ∀
i
≠
j, and ⋃
i
=
1
N
P
i
=
Q, P
=
P
1
,
…
,
P
N is said to be a connected partition of Q.Definition 2 (Voronoi partition).
A connected partitionV
=
V
1
,
…
,
V
N of Q is said to be a Voronoi partition of Q, generated by a vector of distinct points c
∈
Q
N, if for each V
i there exists c
i
∈
V
i that satisfies d
(
c
i
,
q
)
≤
d
(
c
j
,
q
), ∀
q
∈
V
i, ∀
j
≠
i.Definition 3 (nonconvex [22]).
Given any finite setQ
⊂
R
2, point q of ∂
Q is nonconvex if, for all ε
>
0, there exist q
1 and q
2 in B
(
q
,
ε
)
∩
∂
Q so that the open interval (
q
1
,
q
2
) is outside Q.Definition 4 (nonconvex allowable environment [22]).
LetQ
⊂
R
2 be a finite set that is allowable if (i)
Qis compact and connected;
(ii)
∂
Q is continuously differentiable except on a finite number of points;
(iii)
∂
Q has a finite number of nonconvex points which are either isolated points or arcs continuously differentiable everywhere except, possibly, at their end points.In this paper, the search environment is denoted byQ, which is a polygon in R
2. It is assumed that there are several polygonal obstacles H
l
⊂
Q, l
=
1
,
…
,
M within Q. Thus, the overall feasible region is F
=
Q
∖
⋃
i
=
1
M
H
l. Obviously, the search region is nonconvex since it contains several obstacles.


### 2.2. Detection and Communication Model

ConsideringN UAVs in the search region, each UAV is modeled as a nonholonomic point mass moving at a limited speed with a minimum turning radius. The position of UAV i is denoted by p
i
=
[
x
i
,
y
i
]
T, and the heading angle is θ
i. So the state of UAV i can be represented as [
x
i
,
y
i
,
θ
i
].Each UAV is equipped with a camera to take measurements within its Field of View (FoV). Supposing UAVi is aiming at point q at time k (Figure 1), the formulation of FoV is(1)
S
i
,
k
=
x
,
y
∣
q
x
-
d
x
1
≤
x
≤
q
x
+
d
x
2
,
q
y
-
d
y
≤
y
≤
q
y
+
d
y
,
q
x
=
x
i
+
h
cos
⁡
θ
i
tan
⁡
α
s
,
q
y
=
y
i
+
h
sin
⁡
θ
i
tan
⁡
α
s
,
d
x
1
=
h
1
tan
⁡
α
s
-
1
tan
⁡
α
s
+
β
s
,
d
x
2
=
h
1
tan
⁡
α
s
-
β
s
-
1
tan
⁡
α
s
,
d
y
=
h
tan
⁡
γ
s
sin
⁡
α
s
,where h
≤
h
m
a
x represents the flight height of UAV and α
s, β
s, and γ
s are the installed angle and view angles of the camera, related to the physical performance.Figure 1
FoV of the UAV.The search regionQ is then discretized into several cells, defined as o
i. Z
o
i
,
k represents that a target exists in cell o
i at time k and Z
¯
o
i
,
k represents that no target exists in cell o
i at time k. Similarly, D
o
i
,
k represents that the UAV detects a target in cell o
i at time k and D
¯
o
i
,
k represents that the UAV detects no target in cell o
i at time k. Here, we only consider the influence of the Euclidean distance p
i
-
o
i between UAV i and the detected cell o
i. The detection probability for UAV i can be represented by a monotonically decreasing function as(2)
P
D
o
i
,
k
∣
Z
o
i
,
k
=
P
0
i
e
-
λ
P
i
-
o
i
2
,
o
i
∈
FoV
0
,
o
i
∉
FoV
,where P
0
i
∈
(
0,1
] is a constant of the sensing performance for UAV i and λ is a positive constant, set to be 1. All the UAVs are provided with the identical kind of camera sensors, so P
0
i
=
P
0
j
=
P
0 (i
,
j
=
1
,
…
,
N
,
i
≠
j).In this paper, the communication range is assumed to be limited. LetC
N
=
(
ξ
,
ν
) be an undirected graph, where ν
=
1,2
,
…
,
N and ξ
=
(
i
,
j
)
:
i
,
j
∈
ν
} represent the vertex set and the edge set, respectively. The communication network for UAV i at time k is modeled as C
N
i
,
k
=
(
ξ
k
,
ν
). The edge set is ξ
k
=
(
i
,
j
)
:
i
,
j
∈
ν
∣
d
(
p
i
,
p
j
)
≤
R
C, where R
C is the communication range. Therefore, the set of UAV i’s neighbors can be defined as ζ
i
,
k
=
j
∈
ν
∣
(
i
,
j
)
∈
ξ
k. “Neighbor communication strategy” is considered; that is, only when two UAVs reach into each other’s communication range, information exchange may happen. DeLima and Pack [23] have proved that this strategy is more efficient in solving the search problem than that of maintaining a global communication network by restricting UAVs’ mobility.


## 2.1. Basic Notions and Definitions

LetR be the set of real numbers and let R
n denote the n-dimensional Euclidean space. For a given set Q
⊆
R
n, Q
i
n
t and ∂
Q denote the interior and the boundary of Q, respectively. One has Q
i
n
t
=
Q
/
∂
Q. d
(
x
,
y
)
=
x
-
y represents the Euclidean distance between points x
,
y
∈
R
2, and B
(
q
,
ε
) is a closed ball centered at q with radius ε. In the following, some definitions about the environment are presented.Definition 1 (connected partition [21]).
LetQ
⊂
R
2 be a finite set and let P
=
P
1
,
…
,
P
N be a subset of Q; if P
i
∩
P
j
=
∅, ∀
i
≠
j, and ⋃
i
=
1
N
P
i
=
Q, P
=
P
1
,
…
,
P
N is said to be a connected partition of Q.Definition 2 (Voronoi partition).
A connected partitionV
=
V
1
,
…
,
V
N of Q is said to be a Voronoi partition of Q, generated by a vector of distinct points c
∈
Q
N, if for each V
i there exists c
i
∈
V
i that satisfies d
(
c
i
,
q
)
≤
d
(
c
j
,
q
), ∀
q
∈
V
i, ∀
j
≠
i.Definition 3 (nonconvex [22]).
Given any finite setQ
⊂
R
2, point q of ∂
Q is nonconvex if, for all ε
>
0, there exist q
1 and q
2 in B
(
q
,
ε
)
∩
∂
Q so that the open interval (
q
1
,
q
2
) is outside Q.Definition 4 (nonconvex allowable environment [22]).
LetQ
⊂
R
2 be a finite set that is allowable if (i)
Qis compact and connected;
(ii)
∂
Q is continuously differentiable except on a finite number of points;
(iii)
∂
Q has a finite number of nonconvex points which are either isolated points or arcs continuously differentiable everywhere except, possibly, at their end points.In this paper, the search environment is denoted byQ, which is a polygon in R
2. It is assumed that there are several polygonal obstacles H
l
⊂
Q, l
=
1
,
…
,
M within Q. Thus, the overall feasible region is F
=
Q
∖
⋃
i
=
1
M
H
l. Obviously, the search region is nonconvex since it contains several obstacles.


## 2.2. Detection and Communication Model

ConsideringN UAVs in the search region, each UAV is modeled as a nonholonomic point mass moving at a limited speed with a minimum turning radius. The position of UAV i is denoted by p
i
=
[
x
i
,
y
i
]
T, and the heading angle is θ
i. So the state of UAV i can be represented as [
x
i
,
y
i
,
θ
i
].Each UAV is equipped with a camera to take measurements within its Field of View (FoV). Supposing UAVi is aiming at point q at time k (Figure 1), the formulation of FoV is(1)
S
i
,
k
=
x
,
y
∣
q
x
-
d
x
1
≤
x
≤
q
x
+
d
x
2
,
q
y
-
d
y
≤
y
≤
q
y
+
d
y
,
q
x
=
x
i
+
h
cos
⁡
θ
i
tan
⁡
α
s
,
q
y
=
y
i
+
h
sin
⁡
θ
i
tan
⁡
α
s
,
d
x
1
=
h
1
tan
⁡
α
s
-
1
tan
⁡
α
s
+
β
s
,
d
x
2
=
h
1
tan
⁡
α
s
-
β
s
-
1
tan
⁡
α
s
,
d
y
=
h
tan
⁡
γ
s
sin
⁡
α
s
,where h
≤
h
m
a
x represents the flight height of UAV and α
s, β
s, and γ
s are the installed angle and view angles of the camera, related to the physical performance.Figure 1
FoV of the UAV.The search regionQ is then discretized into several cells, defined as o
i. Z
o
i
,
k represents that a target exists in cell o
i at time k and Z
¯
o
i
,
k represents that no target exists in cell o
i at time k. Similarly, D
o
i
,
k represents that the UAV detects a target in cell o
i at time k and D
¯
o
i
,
k represents that the UAV detects no target in cell o
i at time k. Here, we only consider the influence of the Euclidean distance p
i
-
o
i between UAV i and the detected cell o
i. The detection probability for UAV i can be represented by a monotonically decreasing function as(2)
P
D
o
i
,
k
∣
Z
o
i
,
k
=
P
0
i
e
-
λ
P
i
-
o
i
2
,
o
i
∈
FoV
0
,
o
i
∉
FoV
,where P
0
i
∈
(
0,1
] is a constant of the sensing performance for UAV i and λ is a positive constant, set to be 1. All the UAVs are provided with the identical kind of camera sensors, so P
0
i
=
P
0
j
=
P
0 (i
,
j
=
1
,
…
,
N
,
i
≠
j).In this paper, the communication range is assumed to be limited. LetC
N
=
(
ξ
,
ν
) be an undirected graph, where ν
=
1,2
,
…
,
N and ξ
=
(
i
,
j
)
:
i
,
j
∈
ν
} represent the vertex set and the edge set, respectively. The communication network for UAV i at time k is modeled as C
N
i
,
k
=
(
ξ
k
,
ν
). The edge set is ξ
k
=
(
i
,
j
)
:
i
,
j
∈
ν
∣
d
(
p
i
,
p
j
)
≤
R
C, where R
C is the communication range. Therefore, the set of UAV i’s neighbors can be defined as ζ
i
,
k
=
j
∈
ν
∣
(
i
,
j
)
∈
ξ
k. “Neighbor communication strategy” is considered; that is, only when two UAVs reach into each other’s communication range, information exchange may happen. DeLima and Pack [23] have proved that this strategy is more efficient in solving the search problem than that of maintaining a global communication network by restricting UAVs’ mobility.


## 3. Cooperative Search Framework

The goal of this paper is to develop an efficient cooperative search method for multiple UAVs to find several unknown targets in a nonconvex environment with arbitrary obstacles. Inspired by the coverage control, an asynchronous distributed cooperative search framework is proposed by integrating the information update into the coverage control scheme, as shown in Figure2. The proposed framework will promote the detection performance of target search while taking account of coverage performance when the distribution of targets is uniform. Specifically, there are two interdependent tasks. One is online information update containing local update and fusion update. The other is control strategy generation for cooperative search by optimizing the objective function. These two tasks are closely related through the probability map T
P
i and uncertainty map I
i. At first, each UAV takes a measurement s
i and updates its local probability map T
P
i through Bayesian update rule, where the inaccuracies and uncertainties in the measurements are taken into consideration. Then, UAVs in the same communication network C
N
i will exchange the information of the probability map and position, as denoted by the blue lines in Figure 2. The shared information of the probability map is fused by individual UAV to maintain a common view of the search region, and the position information is used to update the Voronoi partition V
i. The uncertainty map I
i will be obtained by using both position and probability map. Finally, each UAV generates its control law by optimizing the objective function H
i. Note that this framework is first designed for the convex environment and can also be used for the nonconvex environment, which will be described in the next section. In the following, the explicit description of the objective function and information update will be presented.Figure 2
Asynchronous distributed cooperative search framework.


### 3.1. Objective Function Formulation

#### 3.1.1. Objective Function and Control Law

First, the cooperative search method in the convex region is presented, which is fundamental to that in the nonconvex region. At each decision, the UAV should determine the optimal control input based on the current state and environment information. The objective is to find several unknown targets, while maximizing the detection probability and minimizing the expected search time. To this end, the following objective function is proposed to be minimized, which is associated with the Voronoi partition:(3)
H
P
=
∫
Q
f
p
i
-
r
γ
r
d
r
=
∑
i
=
1
N
∫
V
i
f
p
i
-
r
γ
r
d
r
,where r
∈
Q is the cell of environment, V
i represents the Voronoi region generated by p
i, and f
(
p
i
-
r
) is a nondecreasing function. The density function γ
(
r
)
:
Q
→
R
+ represents the (relative) importance level of each cell in the environment Q. The more important the cell r is, the larger the value γ
(
r
) would be. Generally speaking, more attention should be paid to the region with a larger density function.In order to optimize the objective functionH
(
P
), the partial derivative of the objective function should be calculated. Supposing f
(
p
i
-
r
)
=
p
i
-
r
2 in this paper, it can be obtained that(4)
∂
H
∂
p
i
=
∂
H
i
p
i
,
V
i
∂
p
i
=
2
∫
V
i
p
i
-
r
γ
r
d
r
=
2
M
V
i
p
i
-
C
V
i
,where M
V
i and C
V
i are the mass and the centroid with respect to the density function γ
(
r
), defined as(5)
M
V
i
=
∫
V
i
γ
r
d
r
,
C
V
i
=
1
M
V
i
∫
V
i
r
γ
r
d
r
.For simplification, the dynamics of each UAV is chosen asp
i
˙
=
u
i. The gradient descent approach is used to minimize the objective function, expressed as(6)
u
i
=
p
i
˙
=
-
k
i
p
i
-
C
V
i
,where k
i is a positive gain.However, in the nonconvex environment, the above control law may fail for two reasons. The centroid of Voronoi region may lie in the unfeasible area, making it impossible to converge to the centroid. The other is that the trajectory generated by the control law (6) may travel through the unfeasible area. Motivated by this, in the next section, a generalized transform will be presented to transform the nonconvex environment into a convex one such that the above control law can be used.


#### 3.1.2. Adaptive Density Function

The density function represents the importance level of each cell in the region, which essentially affects the control strategy of UAVs in the search process. In this paper, we propose an adaptive density functionγ
(
r
) that can utilize different types of the information in different search stages. This density function depends on the real-time updated probability map and uncertainty map.Here, the concept of the uncertainty mapI
i is introduced, defined as the degree of uncertainty about target existence in the environment. The uncertainty map can be derived from the probability map and UAV’s trajectories as follows:(7)
I
i
,
k
r
=
-
TP
i
,
k
r
ln
⁡
TP
i
,
k
r
,
if
TP
i
,
k
r
≠
0
0
,
if
TP
i
,
k
r
=
0
and
r
has  been  detected
0.5
,
if
TP
i
,
k
r
=
0
and
r
has  not  been  detected
.From the above equations, it can be seen that if the target existence probability is close to the middle section of[
0,1
], it will be difficult to determine the existence of the target. Thus, the uncertainty map will have high value in those regions. Besides, the regions that have not been detected by UAVs also have high values, in order to facilitate covering the unvisited regions. Hence, we define the adaptive density function as follows:(8)
γ
r
=
TP
i
,
k
r
,
Stage
T
1
I
i
,
k
r
,
Stage
T
2
,where S
t
a
g
e
T
1 is from the beginning of search to the time when UAV converges to the centroid of its Voronoi region and gets stuck. From that time on, the second S
t
a
g
e
T
2 is started.This adaptive density function can promote the detection performance of target search while taking account of coverage performance. The probability map is used to lead UAV to move towards the region with high target existence probability, in order to find the targets in a shorter time. However, if the UAV has detected all the regions with high target probability and the centroid of Voronoi region no longer updates, the control law in (6) always leads UAV to the centroid. So the UAV may get stuck in the local optimum. Therefore, the uncertainty map instead of the probability map is used to help UAV escape from the local optimum and explore the regions with large uncertainty.


### 3.2. Information Update Method

In the cooperative search problem, each UAV keeps an individual probability mapT
P
i. The information update is in the form of probability map update, and the uncertainty map I
i is derived from the probability map. The source of the UAV’s information consists of three parts: initial information, detection information, and communication information, so the update of individual probability map has two stages: probability map local update and probability map fusion update.


#### 3.2.1. Probability Map Local Update

Considering the uncertainty of sensor measurements, the probability map local update is related with whether a target is detected or not, provided that the target existence probabilityP
(
Z
o
j
,
k
) is known in advance. B
F
O
V
i represents that a target is detected by UAV i at time k with respect to the sensor measurement s
i
,
k
=
1. B
¯
F
O
V
i is the opposite of B
F
O
V
i with respect to measurement s
i
,
k
=
0. Considering the different sensor measurements, the update processes are different.(1) If a target is detected by UAVi, the update of target probability at time k based on Bayesian theory can be expressed as(9)
P
Z
o
j
,
k
+
1
∣
B
FOV
i
=
P
B
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∈
Q
,where P
(
B
FOV
i
∣
Z
o
j
,
k
) is the probability of detecting the target by UAV i, when a target exists in cell o
j, and P
(
Z
o
j
,
0
) is the initialized probability map. Since the FoV of UAV i contains M
i cells, P
(
B
FOV
i
∣
Z
o
j
,
k
) can be described by the probability of failing to detect the target in these M
i cells as(10)
P
B
FOV
i
∣
Z
o
j
,
k
=
1
-
P
B
¯
FOV
i
∣
Z
o
j
,
k
=
1
-
∏
i
=
1
M
i
P
D
¯
o
i
,
k
∣
Z
o
j
,
k
,where P
(
D
¯
o
i
,
k
∣
Z
o
j
,
k
) represents the probability that UAV i fails to detect the target in cell o
i when the target is in cell o
j. Assuming p
a is false alarm probability, P
(
D
¯
o
i
,
k
∣
Z
o
j
,
k
) is described as follows:(11)
P
D
¯
o
i
,
k
∣
Z
o
j
,
k
=
p
¯
a
,
o
i
≠
o
j
P
¯
0
i
,
o
i
=
o
j
,where P
¯
0
i
=
1
-
P
0
i
e
-
λ
i
p
i
-
o
i and p
¯
a
=
1
-
p
a.If the target is in the FoV of UAV, the detection fault of missing the target occurs in one cell with the probabilityP
¯
0
i. If not, the detection result is correct, and only false alarm probability is eliminated. Therefore, the probability for UAV i of failing to detect the target is(12)
∏
i
=
1
M
i
P
D
¯
o
i
,
k
∣
Z
o
j
,
k
=
P
¯
0
i
×
p
¯
a
M
i
-
1
,
o
j
∈
S
i
,
k
1
-
p
¯
a
M
i
,
o
j
∉
S
i
,
k
.Combining (10) with (12), it can be obtained that(13)
P
B
FOV
i
∣
Z
o
j
,
k
=
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
,
o
j
∈
S
i
,
k
1
-
p
¯
a
M
i
,
o
j
∉
S
i
,
k
.Based on the Bayesian complete probability formula,P
(
B
F
O
V
i
) can be expressed as(14)
P
B
FOV
i
=
∑
o
j
∈
S
i
,
k
P
B
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
+
∑
o
j
∉
S
i
,
k
P
B
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
=
∑
o
j
∈
S
i
,
k
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
+
∑
o
j
∉
S
i
,
k
1
-
p
¯
a
M
i
P
Z
o
j
,
k
.Therefore, when the target is detected by UAVi in its FoV, the update expression of the target existence probability can be obtained as follows:(15)
P
Z
o
j
,
k
∣
B
FOV
i
=
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∈
S
i
,
k
1
-
p
¯
a
M
i
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∉
S
i
,
k
.(2) If the target is not detected by UAVi in its FoV, the update expression of the target existence probability can be obtained as follows:(16)
P
Z
o
j
,
k
+
1
∣
B
¯
FOV
i
=
P
B
¯
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∈
Q
.With the same evolving process as mentioned above, it can be obtained that(17)
P
Z
o
j
,
k
+
1
∣
B
¯
FOV
i
=
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∈
S
i
,
k
p
¯
a
M
i
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∉
S
i
,
k
,where P
B
¯
FOV
i
=
∑
o
j
∈
S
i
,
k
‍
P
¯
0
i
×
p
¯
a
M
i
-
1
P
(
Z
o
j
,
k
)
+
∑
o
j
∉
S
i
,
k
‍
p
¯
a
M
i
P
(
Z
o
j
,
k
).As described in (15) and (17), the update of target probability with different sensor measurements at time k for UAV i can be summarized as(18)
TP
i
,
k
+
1
=
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∈
S
i
,
k
and
s
i
,
k
=
1
1
-
p
¯
a
M
i
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∉
S
i
,
k
and
s
i
,
k
=
1
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∈
S
i
,
k
and
s
i
,
k
=
0
p
¯
a
M
i
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∉
S
i
,
k
and
s
i
,
k
=
0
.


#### 3.2.2. Probability Map Fusion Update

Due to the “neighbor communication strategy,” only when UAVi and UAV j fly into each other’s communication range, the information exchange occurs. Each UAV updates its own probability map by fusing information from its neighbors, using the consensus protocol [15](19)
TP
i
,
k
=
∑
j
=
1
N
ω
i
,
j
,
k
TP
i
,
k
,where ω
i
,
j
,
k represents the weight of the target probability map for each UAV, which can be expressed as(20)
ω
i
,
j
,
k
=
1
-
ζ
i
,
k
N
,
if
i
=
j
1
N
,
if
i
≠
j
,
j
∈
ζ
i
,
k
0
,
if
i
≠
j
,
j
∉
ζ
i
,
k
,where ζ
i
,
k is the number of UAV i’s neighbors at time k.


## 3.1. Objective Function Formulation

### 3.1.1. Objective Function and Control Law

First, the cooperative search method in the convex region is presented, which is fundamental to that in the nonconvex region. At each decision, the UAV should determine the optimal control input based on the current state and environment information. The objective is to find several unknown targets, while maximizing the detection probability and minimizing the expected search time. To this end, the following objective function is proposed to be minimized, which is associated with the Voronoi partition:(3)
H
P
=
∫
Q
f
p
i
-
r
γ
r
d
r
=
∑
i
=
1
N
∫
V
i
f
p
i
-
r
γ
r
d
r
,where r
∈
Q is the cell of environment, V
i represents the Voronoi region generated by p
i, and f
(
p
i
-
r
) is a nondecreasing function. The density function γ
(
r
)
:
Q
→
R
+ represents the (relative) importance level of each cell in the environment Q. The more important the cell r is, the larger the value γ
(
r
) would be. Generally speaking, more attention should be paid to the region with a larger density function.In order to optimize the objective functionH
(
P
), the partial derivative of the objective function should be calculated. Supposing f
(
p
i
-
r
)
=
p
i
-
r
2 in this paper, it can be obtained that(4)
∂
H
∂
p
i
=
∂
H
i
p
i
,
V
i
∂
p
i
=
2
∫
V
i
p
i
-
r
γ
r
d
r
=
2
M
V
i
p
i
-
C
V
i
,where M
V
i and C
V
i are the mass and the centroid with respect to the density function γ
(
r
), defined as(5)
M
V
i
=
∫
V
i
γ
r
d
r
,
C
V
i
=
1
M
V
i
∫
V
i
r
γ
r
d
r
.For simplification, the dynamics of each UAV is chosen asp
i
˙
=
u
i. The gradient descent approach is used to minimize the objective function, expressed as(6)
u
i
=
p
i
˙
=
-
k
i
p
i
-
C
V
i
,where k
i is a positive gain.However, in the nonconvex environment, the above control law may fail for two reasons. The centroid of Voronoi region may lie in the unfeasible area, making it impossible to converge to the centroid. The other is that the trajectory generated by the control law (6) may travel through the unfeasible area. Motivated by this, in the next section, a generalized transform will be presented to transform the nonconvex environment into a convex one such that the above control law can be used.


### 3.1.2. Adaptive Density Function

The density function represents the importance level of each cell in the region, which essentially affects the control strategy of UAVs in the search process. In this paper, we propose an adaptive density functionγ
(
r
) that can utilize different types of the information in different search stages. This density function depends on the real-time updated probability map and uncertainty map.Here, the concept of the uncertainty mapI
i is introduced, defined as the degree of uncertainty about target existence in the environment. The uncertainty map can be derived from the probability map and UAV’s trajectories as follows:(7)
I
i
,
k
r
=
-
TP
i
,
k
r
ln
⁡
TP
i
,
k
r
,
if
TP
i
,
k
r
≠
0
0
,
if
TP
i
,
k
r
=
0
and
r
has  been  detected
0.5
,
if
TP
i
,
k
r
=
0
and
r
has  not  been  detected
.From the above equations, it can be seen that if the target existence probability is close to the middle section of[
0,1
], it will be difficult to determine the existence of the target. Thus, the uncertainty map will have high value in those regions. Besides, the regions that have not been detected by UAVs also have high values, in order to facilitate covering the unvisited regions. Hence, we define the adaptive density function as follows:(8)
γ
r
=
TP
i
,
k
r
,
Stage
T
1
I
i
,
k
r
,
Stage
T
2
,where S
t
a
g
e
T
1 is from the beginning of search to the time when UAV converges to the centroid of its Voronoi region and gets stuck. From that time on, the second S
t
a
g
e
T
2 is started.This adaptive density function can promote the detection performance of target search while taking account of coverage performance. The probability map is used to lead UAV to move towards the region with high target existence probability, in order to find the targets in a shorter time. However, if the UAV has detected all the regions with high target probability and the centroid of Voronoi region no longer updates, the control law in (6) always leads UAV to the centroid. So the UAV may get stuck in the local optimum. Therefore, the uncertainty map instead of the probability map is used to help UAV escape from the local optimum and explore the regions with large uncertainty.


## 3.1.1. Objective Function and Control Law

First, the cooperative search method in the convex region is presented, which is fundamental to that in the nonconvex region. At each decision, the UAV should determine the optimal control input based on the current state and environment information. The objective is to find several unknown targets, while maximizing the detection probability and minimizing the expected search time. To this end, the following objective function is proposed to be minimized, which is associated with the Voronoi partition:(3)
H
P
=
∫
Q
f
p
i
-
r
γ
r
d
r
=
∑
i
=
1
N
∫
V
i
f
p
i
-
r
γ
r
d
r
,where r
∈
Q is the cell of environment, V
i represents the Voronoi region generated by p
i, and f
(
p
i
-
r
) is a nondecreasing function. The density function γ
(
r
)
:
Q
→
R
+ represents the (relative) importance level of each cell in the environment Q. The more important the cell r is, the larger the value γ
(
r
) would be. Generally speaking, more attention should be paid to the region with a larger density function.In order to optimize the objective functionH
(
P
), the partial derivative of the objective function should be calculated. Supposing f
(
p
i
-
r
)
=
p
i
-
r
2 in this paper, it can be obtained that(4)
∂
H
∂
p
i
=
∂
H
i
p
i
,
V
i
∂
p
i
=
2
∫
V
i
p
i
-
r
γ
r
d
r
=
2
M
V
i
p
i
-
C
V
i
,where M
V
i and C
V
i are the mass and the centroid with respect to the density function γ
(
r
), defined as(5)
M
V
i
=
∫
V
i
γ
r
d
r
,
C
V
i
=
1
M
V
i
∫
V
i
r
γ
r
d
r
.For simplification, the dynamics of each UAV is chosen asp
i
˙
=
u
i. The gradient descent approach is used to minimize the objective function, expressed as(6)
u
i
=
p
i
˙
=
-
k
i
p
i
-
C
V
i
,where k
i is a positive gain.However, in the nonconvex environment, the above control law may fail for two reasons. The centroid of Voronoi region may lie in the unfeasible area, making it impossible to converge to the centroid. The other is that the trajectory generated by the control law (6) may travel through the unfeasible area. Motivated by this, in the next section, a generalized transform will be presented to transform the nonconvex environment into a convex one such that the above control law can be used.


## 3.1.2. Adaptive Density Function

The density function represents the importance level of each cell in the region, which essentially affects the control strategy of UAVs in the search process. In this paper, we propose an adaptive density functionγ
(
r
) that can utilize different types of the information in different search stages. This density function depends on the real-time updated probability map and uncertainty map.Here, the concept of the uncertainty mapI
i is introduced, defined as the degree of uncertainty about target existence in the environment. The uncertainty map can be derived from the probability map and UAV’s trajectories as follows:(7)
I
i
,
k
r
=
-
TP
i
,
k
r
ln
⁡
TP
i
,
k
r
,
if
TP
i
,
k
r
≠
0
0
,
if
TP
i
,
k
r
=
0
and
r
has  been  detected
0.5
,
if
TP
i
,
k
r
=
0
and
r
has  not  been  detected
.From the above equations, it can be seen that if the target existence probability is close to the middle section of[
0,1
], it will be difficult to determine the existence of the target. Thus, the uncertainty map will have high value in those regions. Besides, the regions that have not been detected by UAVs also have high values, in order to facilitate covering the unvisited regions. Hence, we define the adaptive density function as follows:(8)
γ
r
=
TP
i
,
k
r
,
Stage
T
1
I
i
,
k
r
,
Stage
T
2
,where S
t
a
g
e
T
1 is from the beginning of search to the time when UAV converges to the centroid of its Voronoi region and gets stuck. From that time on, the second S
t
a
g
e
T
2 is started.This adaptive density function can promote the detection performance of target search while taking account of coverage performance. The probability map is used to lead UAV to move towards the region with high target existence probability, in order to find the targets in a shorter time. However, if the UAV has detected all the regions with high target probability and the centroid of Voronoi region no longer updates, the control law in (6) always leads UAV to the centroid. So the UAV may get stuck in the local optimum. Therefore, the uncertainty map instead of the probability map is used to help UAV escape from the local optimum and explore the regions with large uncertainty.


## 3.2. Information Update Method

In the cooperative search problem, each UAV keeps an individual probability mapT
P
i. The information update is in the form of probability map update, and the uncertainty map I
i is derived from the probability map. The source of the UAV’s information consists of three parts: initial information, detection information, and communication information, so the update of individual probability map has two stages: probability map local update and probability map fusion update.


### 3.2.1. Probability Map Local Update

Considering the uncertainty of sensor measurements, the probability map local update is related with whether a target is detected or not, provided that the target existence probabilityP
(
Z
o
j
,
k
) is known in advance. B
F
O
V
i represents that a target is detected by UAV i at time k with respect to the sensor measurement s
i
,
k
=
1. B
¯
F
O
V
i is the opposite of B
F
O
V
i with respect to measurement s
i
,
k
=
0. Considering the different sensor measurements, the update processes are different.(1) If a target is detected by UAVi, the update of target probability at time k based on Bayesian theory can be expressed as(9)
P
Z
o
j
,
k
+
1
∣
B
FOV
i
=
P
B
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∈
Q
,where P
(
B
FOV
i
∣
Z
o
j
,
k
) is the probability of detecting the target by UAV i, when a target exists in cell o
j, and P
(
Z
o
j
,
0
) is the initialized probability map. Since the FoV of UAV i contains M
i cells, P
(
B
FOV
i
∣
Z
o
j
,
k
) can be described by the probability of failing to detect the target in these M
i cells as(10)
P
B
FOV
i
∣
Z
o
j
,
k
=
1
-
P
B
¯
FOV
i
∣
Z
o
j
,
k
=
1
-
∏
i
=
1
M
i
P
D
¯
o
i
,
k
∣
Z
o
j
,
k
,where P
(
D
¯
o
i
,
k
∣
Z
o
j
,
k
) represents the probability that UAV i fails to detect the target in cell o
i when the target is in cell o
j. Assuming p
a is false alarm probability, P
(
D
¯
o
i
,
k
∣
Z
o
j
,
k
) is described as follows:(11)
P
D
¯
o
i
,
k
∣
Z
o
j
,
k
=
p
¯
a
,
o
i
≠
o
j
P
¯
0
i
,
o
i
=
o
j
,where P
¯
0
i
=
1
-
P
0
i
e
-
λ
i
p
i
-
o
i and p
¯
a
=
1
-
p
a.If the target is in the FoV of UAV, the detection fault of missing the target occurs in one cell with the probabilityP
¯
0
i. If not, the detection result is correct, and only false alarm probability is eliminated. Therefore, the probability for UAV i of failing to detect the target is(12)
∏
i
=
1
M
i
P
D
¯
o
i
,
k
∣
Z
o
j
,
k
=
P
¯
0
i
×
p
¯
a
M
i
-
1
,
o
j
∈
S
i
,
k
1
-
p
¯
a
M
i
,
o
j
∉
S
i
,
k
.Combining (10) with (12), it can be obtained that(13)
P
B
FOV
i
∣
Z
o
j
,
k
=
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
,
o
j
∈
S
i
,
k
1
-
p
¯
a
M
i
,
o
j
∉
S
i
,
k
.Based on the Bayesian complete probability formula,P
(
B
F
O
V
i
) can be expressed as(14)
P
B
FOV
i
=
∑
o
j
∈
S
i
,
k
P
B
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
+
∑
o
j
∉
S
i
,
k
P
B
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
=
∑
o
j
∈
S
i
,
k
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
+
∑
o
j
∉
S
i
,
k
1
-
p
¯
a
M
i
P
Z
o
j
,
k
.Therefore, when the target is detected by UAVi in its FoV, the update expression of the target existence probability can be obtained as follows:(15)
P
Z
o
j
,
k
∣
B
FOV
i
=
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∈
S
i
,
k
1
-
p
¯
a
M
i
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∉
S
i
,
k
.(2) If the target is not detected by UAVi in its FoV, the update expression of the target existence probability can be obtained as follows:(16)
P
Z
o
j
,
k
+
1
∣
B
¯
FOV
i
=
P
B
¯
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∈
Q
.With the same evolving process as mentioned above, it can be obtained that(17)
P
Z
o
j
,
k
+
1
∣
B
¯
FOV
i
=
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∈
S
i
,
k
p
¯
a
M
i
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∉
S
i
,
k
,where P
B
¯
FOV
i
=
∑
o
j
∈
S
i
,
k
‍
P
¯
0
i
×
p
¯
a
M
i
-
1
P
(
Z
o
j
,
k
)
+
∑
o
j
∉
S
i
,
k
‍
p
¯
a
M
i
P
(
Z
o
j
,
k
).As described in (15) and (17), the update of target probability with different sensor measurements at time k for UAV i can be summarized as(18)
TP
i
,
k
+
1
=
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∈
S
i
,
k
and
s
i
,
k
=
1
1
-
p
¯
a
M
i
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∉
S
i
,
k
and
s
i
,
k
=
1
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∈
S
i
,
k
and
s
i
,
k
=
0
p
¯
a
M
i
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∉
S
i
,
k
and
s
i
,
k
=
0
.


### 3.2.2. Probability Map Fusion Update

Due to the “neighbor communication strategy,” only when UAVi and UAV j fly into each other’s communication range, the information exchange occurs. Each UAV updates its own probability map by fusing information from its neighbors, using the consensus protocol [15](19)
TP
i
,
k
=
∑
j
=
1
N
ω
i
,
j
,
k
TP
i
,
k
,where ω
i
,
j
,
k represents the weight of the target probability map for each UAV, which can be expressed as(20)
ω
i
,
j
,
k
=
1
-
ζ
i
,
k
N
,
if
i
=
j
1
N
,
if
i
≠
j
,
j
∈
ζ
i
,
k
0
,
if
i
≠
j
,
j
∉
ζ
i
,
k
,where ζ
i
,
k is the number of UAV i’s neighbors at time k.


## 3.2.1. Probability Map Local Update

Considering the uncertainty of sensor measurements, the probability map local update is related with whether a target is detected or not, provided that the target existence probabilityP
(
Z
o
j
,
k
) is known in advance. B
F
O
V
i represents that a target is detected by UAV i at time k with respect to the sensor measurement s
i
,
k
=
1. B
¯
F
O
V
i is the opposite of B
F
O
V
i with respect to measurement s
i
,
k
=
0. Considering the different sensor measurements, the update processes are different.(1) If a target is detected by UAVi, the update of target probability at time k based on Bayesian theory can be expressed as(9)
P
Z
o
j
,
k
+
1
∣
B
FOV
i
=
P
B
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∈
Q
,where P
(
B
FOV
i
∣
Z
o
j
,
k
) is the probability of detecting the target by UAV i, when a target exists in cell o
j, and P
(
Z
o
j
,
0
) is the initialized probability map. Since the FoV of UAV i contains M
i cells, P
(
B
FOV
i
∣
Z
o
j
,
k
) can be described by the probability of failing to detect the target in these M
i cells as(10)
P
B
FOV
i
∣
Z
o
j
,
k
=
1
-
P
B
¯
FOV
i
∣
Z
o
j
,
k
=
1
-
∏
i
=
1
M
i
P
D
¯
o
i
,
k
∣
Z
o
j
,
k
,where P
(
D
¯
o
i
,
k
∣
Z
o
j
,
k
) represents the probability that UAV i fails to detect the target in cell o
i when the target is in cell o
j. Assuming p
a is false alarm probability, P
(
D
¯
o
i
,
k
∣
Z
o
j
,
k
) is described as follows:(11)
P
D
¯
o
i
,
k
∣
Z
o
j
,
k
=
p
¯
a
,
o
i
≠
o
j
P
¯
0
i
,
o
i
=
o
j
,where P
¯
0
i
=
1
-
P
0
i
e
-
λ
i
p
i
-
o
i and p
¯
a
=
1
-
p
a.If the target is in the FoV of UAV, the detection fault of missing the target occurs in one cell with the probabilityP
¯
0
i. If not, the detection result is correct, and only false alarm probability is eliminated. Therefore, the probability for UAV i of failing to detect the target is(12)
∏
i
=
1
M
i
P
D
¯
o
i
,
k
∣
Z
o
j
,
k
=
P
¯
0
i
×
p
¯
a
M
i
-
1
,
o
j
∈
S
i
,
k
1
-
p
¯
a
M
i
,
o
j
∉
S
i
,
k
.Combining (10) with (12), it can be obtained that(13)
P
B
FOV
i
∣
Z
o
j
,
k
=
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
,
o
j
∈
S
i
,
k
1
-
p
¯
a
M
i
,
o
j
∉
S
i
,
k
.Based on the Bayesian complete probability formula,P
(
B
F
O
V
i
) can be expressed as(14)
P
B
FOV
i
=
∑
o
j
∈
S
i
,
k
P
B
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
+
∑
o
j
∉
S
i
,
k
P
B
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
=
∑
o
j
∈
S
i
,
k
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
+
∑
o
j
∉
S
i
,
k
1
-
p
¯
a
M
i
P
Z
o
j
,
k
.Therefore, when the target is detected by UAVi in its FoV, the update expression of the target existence probability can be obtained as follows:(15)
P
Z
o
j
,
k
∣
B
FOV
i
=
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∈
S
i
,
k
1
-
p
¯
a
M
i
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∉
S
i
,
k
.(2) If the target is not detected by UAVi in its FoV, the update expression of the target existence probability can be obtained as follows:(16)
P
Z
o
j
,
k
+
1
∣
B
¯
FOV
i
=
P
B
¯
FOV
i
∣
Z
o
j
,
k
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∈
Q
.With the same evolving process as mentioned above, it can be obtained that(17)
P
Z
o
j
,
k
+
1
∣
B
¯
FOV
i
=
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∈
S
i
,
k
p
¯
a
M
i
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∉
S
i
,
k
,where P
B
¯
FOV
i
=
∑
o
j
∈
S
i
,
k
‍
P
¯
0
i
×
p
¯
a
M
i
-
1
P
(
Z
o
j
,
k
)
+
∑
o
j
∉
S
i
,
k
‍
p
¯
a
M
i
P
(
Z
o
j
,
k
).As described in (15) and (17), the update of target probability with different sensor measurements at time k for UAV i can be summarized as(18)
TP
i
,
k
+
1
=
1
-
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∈
S
i
,
k
and
s
i
,
k
=
1
1
-
p
¯
a
M
i
P
Z
o
j
,
k
P
B
FOV
i
,
o
j
∉
S
i
,
k
and
s
i
,
k
=
1
P
¯
0
i
×
p
¯
a
M
i
-
1
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∈
S
i
,
k
and
s
i
,
k
=
0
p
¯
a
M
i
P
Z
o
j
,
k
P
B
¯
FOV
i
,
o
j
∉
S
i
,
k
and
s
i
,
k
=
0
.


## 3.2.2. Probability Map Fusion Update

Due to the “neighbor communication strategy,” only when UAVi and UAV j fly into each other’s communication range, the information exchange occurs. Each UAV updates its own probability map by fusing information from its neighbors, using the consensus protocol [15](19)
TP
i
,
k
=
∑
j
=
1
N
ω
i
,
j
,
k
TP
i
,
k
,where ω
i
,
j
,
k represents the weight of the target probability map for each UAV, which can be expressed as(20)
ω
i
,
j
,
k
=
1
-
ζ
i
,
k
N
,
if
i
=
j
1
N
,
if
i
≠
j
,
j
∈
ζ
i
,
k
0
,
if
i
≠
j
,
j
∉
ζ
i
,
k
,where ζ
i
,
k is the number of UAV i’s neighbors at time k.


## 4. Nonconvex Environment with Arbitrary Obstacles

Considering the search environmentQ
⊆
R
2 is nonconvex with arbitrary obstacles H, we will propose a cooperative search method to address the search problem in the nonconvex environment in this section. Generally, there are two ways to address this problem. One is to find a different objective function that works for nonconvex environments by taking path planning into account, such as the potential field-based method [19, 20]. The second is to find a transformation for the environment [16, 17]. Here, we will choose the second one in order to use the cooperative search framework and formulations in the above section.Our transformation is inspired by the diffeomorphism idea [16]. However, the diffeomorphism in [16] can only be applied to the nonconvex environment with the convex obstacles. Considering more realistic application, we propose a generalized transformation method that can transform the nonconvex region with arbitrary obstacles (including convex and concave obstacles) into an almost convex region. First, the definition of diffeomorphism is given, as well as the lemma in [16].Definition 5 (diffeomorphsim [24]).
Given two manifoldsM and N, a differentiable map f
:
M
→
N is called a diffeomorphism, if it is a bijection and its inverse f
-
1
:
N
→
M is differentiable as well.Lemma 6.
For a given setQ
⊆
R
2, let R
~
⊂
R
2 be a compact, convex, and simply connected partition of Q and let H
⊆
R
~ be a convex open set. Let ∂
H
∩
∂
R
~
=
∅ and R
=
R
~
/
H. Given p
′
∈
H
i
n
t, there is a.e. C
1 diffeomorphism φ between R
i
n
t and R
~
i
n
t
∖
{
p
′
}.The above lemma follows thatH is convex or star-shaped. However, the real obstacle H may be arbitrary; the diffeomorphism method in Lemma 6 cannot be used in this case. Therefore, we propose a generalized method to construct the transformation based on a global convex decomposition and some extensions of diffeomorphism in the following.


### 4.1. Diffeomorphism Based Transformation

#### 4.1.1. Global Convex Decomposition

As the main preprocessing procedure for constructing the transformation, the global convex decomposition method decomposes the original nonconvex region and the concave obstacles, making it possible to find the diffeomorphism based transformation in the nonconvex environment with concave obstacles.First, the convex decomposition of concave obstacles is given. Some works have been conducted on the convex decomposition of concave polygons. Schachter [25] utilizes the Delaunay method to decompose polygons into convex sets. Chazelle and Dobkin [26] construct the concave point sets of the x
k type to realize the convex decomposition. Keil [27] uses improved dynamic programming algorithms to achieve decomposing a simple polygon. In this paper, our purpose for convex decomposition of concave polygons is to find the diffeomorphism, which is related to the visibility. So a point visibility based convex decomposition method is presented, which can improve the shape quality of the decomposed polygons.(1) Find the Concave Points. Since convex decomposition of concave polygon is to eliminate the concave points, the decomposing line is always drawn from concave points. So the concave points of the concave polygon should be found. The concave point is defined as in Definition 3. The most common method for determining whether a point is concave is to use the vector cross product. Choose the successive points p
i
-
1, p
i, and p
i
+
1 from the point sequence of a polygon in counterclockwise direction. Obtain the vectors p
i
-
1
p
i
→ and p
i
p
i
+
1
→, and then compute the cross product(21)
M
=
x
i
-
x
i
-
1
×
y
i
+
1
-
y
i
-
x
i
+
1
-
x
i
×
y
i
-
y
i
-
1
,where x
i and y
i are the coordinates of point p
i. If M is negative, then p
i is a concave point; otherwise, p
i is a convex point.(2) Find the Visible Point Set. Point q
∈
Q is visible from p
∈
Q if [
p
,
q
]
⊂
Q, where [
p
,
q
] represents the closed path segment [
p
,
q
]
=
p
+
λ
(
q
-
p
)
∣
λ
∈
[
0,1
]. The visibility set is the set of points in Q visible from p. We connect P
i with other vertices of the polygon in counterclockwise direction. For a given vertex P
j, if all the rest vertices behind P
j are on the left side of line P
i
P
j
¯, then P
j is visible from point P
i, as shown in Figure 3.Figure 3
Visible points fromP
1.(3) Find the Decomposing Line Based on the Weighting Function. The decomposing line is obtained from the concave point to the visible point which is carefully searched from the visible point set of this concave point by using a weight function. In [28], a weight function is defined as f
(
α
,
β
)
=
α
-
β, where α and β are the angles between the candidate decomposing line and the neighbor lines. Since α
,
β lie within the interval [
0
,
π
] and the cosine function is monotonous in this interval, we use the following weight function for computational convenience:(22)
f
α
,
β
=
cos
⁡
α
-
cos
⁡
β
.The visible point, from which the decomposing line has the minimum weight value, is chosen. For each polygon, repeat the above steps until all the points of the decomposed polygons are convex.Further, according to the decomposing lines and the edges of the decomposed obstacles, the whole search region is divided into a connected partitionR
~
=
R
~
1
,
…
,
R
~
K of Q. Each partition contains no more than one decomposed convex obstacle. The global decomposition method of the whole region is as follows.First, extend the decomposing lines from two endpoints to intersect with region boundaries, defined as the dividing lines (for the whole region decomposition), while giving up extending those endpoints that will lead to intersecting with obstacles or other extending lines. Then, from each of those abandoned endpoints, draw a radial line along the two edges of the original concave obstacle to intersect with the region boundary or the obstacle boundary, and obtain the dividing lines. By these dividing lines, the whole region is divided into several subregions, each of which containing no more than one decomposed obstacle. Note that the feasible region of a new subregion may be separated into two parts by the decomposed obstacle. In this case, the subregion is divided again to generate a new small subregion without obstacles and the larger one with the obstacle.


#### 4.1.2. Extension of Diffeomorphism

Through convex decomposition, the whole region is divided into several subregions with no more than one convex obstacle. The method for constructing the diffeomorphism in Lemma6 requires that ∂
H
∩
∂
R
~
=
∅. However, the new decomposed obstacles may have the common edges with the subregions such that ∂
H
i
∩
∂
R
~
i
≠
∅. So the method in Lemma 6 is not appropriate. The expression of the “common edge” means that the edge of the subregion contains the edge of the obstacle. In the following, we will present a generalized method to construct the diffeomorphism for the new decomposed subregions.LetR
~
=
R
~
1
,
…
,
R
~
K be a compact, convex, and simply connected partition of Q, representing the subregions. For each R
~
i, let H
i
⊂
R
~
i be a convex open set, representing the decomposed obstacle. Thus, the feasible region is R
i
=
R
~
i
∖
H
i. For this problem, the edges generated by the decomposing lines must be considered. Let ∂
R
i
i
n
t
=
∂
R
i
∖
(
∂
R
i
∩
∂
Q
) be the boundary ∂
R
i but not belong to ∂
Q. Similarly, let ∂
R
~
i
i
n
t
=
∂
R
~
i
∖
(
∂
R
~
i
∩
∂
Q
) be the boundary ∂
R
~
i but not belong to ∂
Q.Proposition 7.
Givenp
i
′
∈
H
i, there is a.e. C
1 diffeomorphism φ
i between R
i
i
n and R
~
i
i
n
∖
{
p
′
i
}, where R
i
i
n and R
~
i
i
n are the simplified representations of R
i
i
n
t
∪
∂
R
i
i
n
t and R
~
i
i
n
t
∪
∂
R
~
i
i
n
t, respectively.Proof.
The situation is categorized according to the number of the common edges of the subregion and the obstacle.
Situation  1 (no common edge). Since ∂
H
i
∩
∂
R
~
i
=
∅, the diffeomorphism can be constructed using the method in the proof of Lemma  1 in [16].
Situation  2 (one common edge in Figure4(a)). Draw a radial line l from p
i
′ to cross H
i and intersect with the boundaries of R
~
i and H
i at the points p
1 and p
0, respectively. Note that if the radial line lies on R
~
i
i
n, p
0 and p
1 are the vertices of H
i and R
~
i on the corresponding boundary, that is, points a
1(a
2) and b
1(b
2). Obviously, it can be obtained that p
1
≠
p
0. Then a simple diffeomorphism φ
i can be constructed as(23)
φ
i
q
=
p
1
-
p
i
′
p
1
-
p
i
′
-
p
0
-
p
i
′
q
-
p
i
′
-
p
0
-
p
i
′
.So any q
∈
R
i
i
n on the radial line l can be mapped to φ
i
(
q
)
∈
R
~
i
i
n
∖
{
p
i
′
}. Due to our hypothesis on R
~
i and R
i, it can be seen that φ
i is continuous and differentiable everywhere in R
i
i
n, and it has a unique inverse φ
i
-
1, which is also continuous and differentiable everywhere in R
~
i
i
n
∖
{
p
i
′
}. Therefore, φ
i is a C
1 diffeomorphism. Note that the intersection points p
1 and p
0 are the functions of q. So even if q
1 and q
2 are equidistant to p
i
′, we still have φ
i
(
q
1
)
≠
φ
i
(
q
2
). Hence, the inverse exists and φ
i is a diffeomorphism.
Situation  3 (two common edges in Figure4(b)). The diffeomorphism is constructed as (23) in Situation  2.
Situation  4 (three common edges). There are two cases.
(i) For the first case in Figure4(d), we first define the critical vertex as the common vertex of the subregion and obstacle, while one of its adjacent vertices does not belong to the obstacle, such as p
i
′ and p
2. Choose one of the critical vertices as p
i
′, then p
′
′ is on the same edge with the other critical vertex, and it does not belong to the obstacle. Thus the convex region is divided into two subregions by the blue dividing line connecting p
i
′ with p
′
′. The first diffeomorphism φ
i
a
(
q
) can be constructed by drawing a radial line from p
i
′ using (23) as in Situation  2. For the other part, the diffeomorphism should be modified a little by drawing a radial line from p
′
′, as follows:(24)
φ
i
b
q
=
p
1
-
p
′
′
p
0
-
p
′
′
q
-
p
′
′
,where p
1 and p
0 are the two intersection points of the radial line with the boundaries of the obstacle, and if q lies on the edge ∂
R
~
i
i
n, p
0 and p
1 are the corresponding vertices of the obstacle.
Note that, through the diffeomorphismφ
i
b
(
q
), all the points (except p
i
′ and p
′
′) on the dividing line have the same mapping with the first diffeomorphism φ
i
a
(
q
). And p
′
′ has a mapping by the diffeomorphism φ
i
a. So we obtain φ
i
=
φ
i
a
∪
φ
i
b that can map R
i
i
n to R
~
i
i
n
∖
p
i
′.
(ii) For the second case in Figure4(e), two edges of H
i are parallel which is a typical situation in structural environment. First, the convex polygon is divided into two partitions such that only one partition contains the obstacle. For the partition without any obstacle, the diffeomorphism is constructed as φ
i
a
(
q
)
=
q. For the other partition, a parallel line is drawn from q and intersects with the boundaries of R
~
i at the left point p
0 and right point p
3. Given the vertices of the polygon, p
0 lying on the line b
1
c
1
¯ can be mapped to p
1 on the line p
i
′
c
1
¯, described as p
1
=
p
i
′
-
(
p
i
′
-
c
1
b
1
-
p
0
)
/
(
b
1
-
c
1
). Draw another parallel line from p
1, and it intersects with the boundary of R
~
i at point p
2. Thus the diffeomorphism can be constructed as φ
i
b
(
q
)
=
p
1
+
(
q
-
p
0
p
2
-
p
1
)
/
(
p
3
-
p
0
). Combining with the representation of p
1, the diffeomorphism can be described as(25)
φ
i
b
q
=
p
i
′
-
p
i
′
-
c
1
b
1
-
p
0
b
1
-
c
1
+
q
-
p
0
p
2
-
p
1
p
3
-
p
0
.For this case, φ
i
=
φ
i
a
∪
φ
i
b.
Situation  5 (more than three common edges in Figure4(e)). The diffeomorphism is constructed as in Situation  4.
Situation  6 (a special case). All edges of the obstacle are common edges in Figure 4(f). The diffeomorphism is constructed as (24) by choosing the only vertex of the subregion which does not belong to obstacle as p
′
′.
As claimed in the above six situations, it is obvious that the diffeomorphismφ
i can be constructed for each subregion.Figure 4
Diffeomorphism for the new decomposed regions.
(a)
(b)
(c)
(d)
(e)
(f)Remark 8.
The subregion containing no obstacle will be transformed by an identical transformation, which is also a diffeomorphism.Remark 9.
Let⋃
i
=
1
M
p
i
′
=
p
′. For each decomposed subregion, the diffeomorphism φ
i is constructed. So the transformation φ
:
F
i
n
t
→
Q
i
n
t
/
p
′ is defined by φ
(
q
)
=
φ
i
(
q
), where q
∈
R
i
i
n.Remark 10.
Given the search region and obstacles, the diffeomorphism can be computed off-line and stored in a table. Therefore, in the search procedure, the diffeomorphism can be directly used by looking up the table, which will significantly reduce the computation time.


#### 4.1.3. Numerical Examples

Two examples of applying the convex decomposition method to divide the nonconvex region with arbitrary obstacles are shown in Figure5. The green lines are the decomposing lines for convex decomposition of the concave obstacle, and the black lines are the dividing lines for the whole environment. In Figure 5(a), the concave obstacle is divided into three convex polygons. And the whole region is divided into four subregions, three of which containing one obstacle and one containing none. In Figure 5(b), a more complex nonconvex region with two obstacles is given. First, the whole region is divided into two partitions. Then the partition with a concave obstacle is divided into six subregions by using the decomposition method.Figure 5
Convex decomposition.
(a)
Nonconvex regionQ
1 with a typical concave obstacle
(b)
Nonconvex regionQ
2 with a convex obstacle and a complex concave obstacleFurther, the diffeomorphism is constructed for each example region by using the method in Proposition7. We sample some points in the original regions and then map them to the transformed regions by the diffeomorphism, as shown in Figure 6.Figure 6
Diffeomorphism of the nonconvex region. (a) and (c) are the two examples of the original nonconvex region with 10 randomly sampled points. (b) and (d) are the transformed convex regions with the transformed image of the sampled points by diffeomorphism.
(a)
Sampled points in originalQ
1
(b)
Point images in transformedφ
(
Q
1
)
(c)
Sampled points in originalQ
2
(d)
Point images in transformedφ
(
Q
2
)


### 4.2. Cooperative Search Method in Nonconvex Environment with Arbitrary Obstacles

The diffeomorphism that transforms a complex nonconvex region with arbitrary obstacles to a convex region has been constructed. Now we will propose the cooperative search method in the nonconvex environment based on the transformation.


#### 4.2.1. Objective Function Transformation

The objective function (3) should be modified by shifting the allowable search domain from the feasible region F
i
n
t to the whole region Q
i
n
t
∖
{
p
′
}. Based on the modification of Variables Theorem [29], it claims that if Ω
⊂
R
n is open, g
1
:
Ω
→
R is integrable over Ω, and g
2
:
Ω
~
→
Ω is a C
1 diffeomorphism, then g
1
∘
g
2
det
⁡
J
g
2 is integrable over Ω
~, and(26)
∫
Ω
g
1
ω
d
ω
=
∫
Ω
~
g
1
∘
g
2
ω
~
det
⁡
J
g
2
ω
~
d
ω
~
,where det
⁡
J
g
2
ω
~ is the determinant of the Jacobian of g
2. In this paper, we assume that g
1
=
f
(
p
i
-
r
)
γ
(
r
)
:
F
i
n
t
→
R and g
2
=
φ
-
1
:
Q
i
n
t
/
(
p
′
)
→
R
i
n
t. Since φ
-
1
(
ω
)
=
r, the transformed objective function in the new convex region can be expressed as(27)
∫
F
i
n
t
f
p
i
-
r
γ
r
d
r
=
∫
Q
i
n
t
∖
p
′
f
φ
-
1
q
i
-
φ
-
1
ω
γ
φ
-
1
ω
det
⁡
J
φ
-
1
ω
d
ω
.It can be seen that the allowable search domain is transformed from the feasible regionF
i
n
t to the whole region Q
i
n
t
∖
p
′, while keeping the function value unchanged. However, the control law for the above objective function (the right side of (27)) is quite hard to be obtained. In order to use the Lloyd algorithm, the new objective function is defined as(28)
H
′
P
=
∫
F
i
n
t
f
φ
p
i
-
φ
r
γ
r
d
r
.Sinceφ is the diffeomorphism, F
i
n
t and φ
(
F
i
n
t
) are topologically equivalent, and the problem remains essentially unchanged. Further, the corresponding function in the transformed domain is obtained:(29)
H
W
,
V
~
=
∫
Q
i
n
t
∖
p
′
f
q
i
-
ω
γ
φ
-
1
ω
det
⁡
J
φ
-
1
ω
d
ω
.So the gradient descent approach as (6) can be used to minimize the objective function(30)
q
˙
i
=
-
k
i
1
q
i
-
C
~
V
i
,where k
i
2 is a positive gain, q
i
=
φ
(
p
i
) is the transformed position of p
i by the diffeomorphism in the transformed region Q
i
n
t
/
(
p
′
), and C
~
V
i is the centroid of the Voronoi region generated by q
i.


#### 4.2.2. Control Strategy for Cooperative Search in the Nonconvex Environment

Since the diffeomorphismφ is a bijection, the trajectory of q
i in Q
i
n
t
/
p
′ has the unique preimage in original region F
i
n
t, which is the trajectory of p
i. So the control law for the new objective function (28) in the original region can be obtained by the inverse mapping of (30):(31)
p
i
˙
=
-
k
i
2
J
φ
-
1
φ
p
i
p
i
-
φ
-
1
C
~
V
i
,where k
i
2 is a positive gain and J
φ is the Jacobian of φ.The following is an analysis of the convergence of the control law (31) using LaSalle’s invariance principle [30], which is basically an extension of Lyapunov’s theorem [31].Theorem 11.
In the nonconvex region, the trajectories of UAVs governed by the control law (31), starting from any initial condition P
(
0
)
=
(
p
1,0
,
…
,
p
N
,
0
)
∈
F
i
n
t, will asymptotically converge to the desired location.Proof.
Since the diffeomorphismφ is a bijection, the trajectory q
i in Q
i
n
t
∖
p
′ has the unique preimage p
i in original region F
i
n
t. Therefore, to prove that the trajectory p
i governed by the control law (31), starting from any initial condition P
(
0
), can converge to φ
-
1
(
C
~
V
i
), we only need to prove that the trajectory q
i governed by the control law (30), starting from initial condition W
(
0
)
=
φ
(
P
(
0
)
)
=
(
φ
(
p
1,0
)
,
…
,
φ
(
p
N
,
0
)
), can converge to C
~
V
i.
AssumeV
(
Q
)
=
H
(
W
,
V
~
)
=
∑
i
=
1
N
‍
∫
V
~
i
‍
q
i
-
ω
2
γ
φ
-
1
ω
det
⁡
J
φ
-
1
ω
d
ω, where Q
=
(
q
1
,
q
2
,
…
,
q
N
) represent the transformed positions of UAVs in Q
i
n
t
∖
p
′ through the diffeomorphism φ. Consider(32)
V
˙
Q
=
d
H
W
,
V
~
d
t
=
∑
i
=
1
N
∂
H
W
,
V
~
∂
q
i
q
˙
i
=
∑
i
=
1
N
2
λ
M
~
V
i
q
i
-
C
~
V
i
-
k
i
1
q
i
-
C
~
V
i
=
-
2
λ
1
k
i
1
∑
i
=
1
N
M
~
V
i
q
i
-
C
~
V
i
2
⩽
0
,where λ
1
=
det
⁡
J
φ
-
1
ω and M
~
V
i is the mass of the Voronoi region generated by q
i in the convex region Q
i
n
t
∖
p
′.
The following can be observed:(i)
V
:
Q
i
n
t
∖
p
′
→
R is continuously differentiable in Q
i
n
t
∖
p
′.
(ii)
Q
i
n
t
∖
p
′ is a compact invariant set.
(iii)
V
˙ is negative semidefinite in Q
i
n
t
∖
p
′.
(iv)
E
=
V
˙
-
1
(
0
)
=
C
~
V
i.
(v)
E itself is the largest invariant subset of E by the control law (30).
Thus, by LaSalle’s invariance principle, the trajectoriesq
i governed by (30), starting from initial configuration W
(
0
), will asymptotically converge to the critical points of H
(
W
,
V
~
); that is, the trajectories of UAVs governed by the control law (31), starting from any initial condition P
(
0
)
=
(
p
1,0
,
…
,
p
N
,
0
)
∈
F
i
n
t, will asymptotically converge to the desired location.Therefore, the discrete-time motion model and the real control law for UAVi in the original nonconvex region can be written as(33)
p
i
,
m
=
p
i
,
m
-
1
+
u
i
,
m
,
u
i
,
m
=
-
k
i
2
J
φ
-
1
φ
p
i
,
m
-
1
p
i
,
m
-
1
-
φ
-
1
C
~
V
i
.Further, considering the mobility of the UAV, the control law must meet some constraints such that the generated trajectories are kinematically feasible:(1) Maximum velocity: the maximum velocityv
m
a
x will limit the control law by the following saturation rule:(34)
u
i
,
m
=
u
i
,
m
,
u
i
,
m
≤
v
max
⁡
v
max
⁡
u
i
,
m
u
i
,
m
,
Otherwise
.(2) Maximum turning angle: the maximum turning angle during a unit time periodΔ
t is derived from the turning rate ω
m
a
x, described as θ
m
a
x
=
ω
m
a
x
Δ
t. Thus the control law is limited by the following equations:(35)
u
i
,
m
=
u
i
,
m
,
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
≤
θ
max
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
-
1
A
,
Otherwise
,where A is the rotating matrix for the maximum turning angle from the last control input as (36)
A
=
cos
⁡
α
-
sin
⁡
α
sin
⁡
α
cos
⁡
α
,
α
=
θ
max
⁡
,
u
i
,
m
·
u
i
,
m
-
1
≥
0
-
θ
max
⁡
,
u
i
,
m
·
u
i
,
m
-
1
<
0represents the rotating direction.Combining (34) with (35), the control law can be expressed as(37)
u
i
,
m
=
u
i
,
m
,
u
i
,
m
≤
v
max
⁡
and
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
≤
θ
max
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
-
1
A
,
u
i
,
m
≤
v
max
⁡
and
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
>
θ
max
⁡
v
max
⁡
u
i
,
m
u
i
,
m
,
u
i
,
m
>
v
max
⁡
and
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
≤
θ
max
⁡
v
max
⁡
u
i
,
m
-
1
u
i
,
m
-
1
A
,
Otherwise
,where u
i
,
m
=
-
k
i
2
J
φ
-
1
(
φ
(
p
i
,
m
-
1
)
)
(
p
i
,
m
-
1
-
φ
-
1
(
C
~
V
i
)
), as shown in (33).


#### 4.2.3. Cooperative Search Algorithm in the Nonconvex Environment with Arbitrary Obstacles

The cooperative search algorithm for the nonconvex region with arbitrary obstacles is presented as Algorithm1.Algorithm 1.
Cooperative search algorithm in the nonconvex environment with arbitrary obstacles is as follows:(1)
Initially, all the UAVs are randomly deployed in the feasible region. The probability map and uncertainty map of each UAV are initialized.
(2)
The nonconvex region with arbitrary obstacles is divided by the global decomposition method, and the diffeomorphismφ is constructed as described in Proposition 7.
(3)
Through the diffeomorphismφ, the original nonconvex region F
i
n
t is transformed into the new convex one, Q
i
n
t
∖
p
′. And the positions of UAVs are also transformed into new virtual ones.
(4)
The Voronoi partition of transformed regionQ
i
n
t
∖
p
′ is generated by the virtual positions of UAVs. Compute the centroid C
~
V
i of each Voronoi region using probability map as the density.
(5)
At each time step, UAVi moves towards the preimages of its generalized Voronoi centroid φ
-
1
(
C
~
V
i
) in the original region, governed by the control law (31).
(6)
Update the probability mapT
P
i of UAV i. If the target probability of a cell is greater than P
m
a
x, a target is found.
(7)
Update the Voronoi partitions based on the new virtual positions of UAVs. Compute the generalized Voronoi centroidsC
~
V
i.
(8)
If UAVi reaches the preimage of its Voronoi centroid φ
-
1
(
C
~
V
i
) and does not move, then use its uncertainty map instead of probability map to compute the new Voronoi centroid φ
-
1
(
C
~
V
i
) from this time step.
(9)
If the terminal condition is satisfied, the algorithm ends; otherwise, go to(
5
).


## 4.1. Diffeomorphism Based Transformation

### 4.1.1. Global Convex Decomposition

As the main preprocessing procedure for constructing the transformation, the global convex decomposition method decomposes the original nonconvex region and the concave obstacles, making it possible to find the diffeomorphism based transformation in the nonconvex environment with concave obstacles.First, the convex decomposition of concave obstacles is given. Some works have been conducted on the convex decomposition of concave polygons. Schachter [25] utilizes the Delaunay method to decompose polygons into convex sets. Chazelle and Dobkin [26] construct the concave point sets of the x
k type to realize the convex decomposition. Keil [27] uses improved dynamic programming algorithms to achieve decomposing a simple polygon. In this paper, our purpose for convex decomposition of concave polygons is to find the diffeomorphism, which is related to the visibility. So a point visibility based convex decomposition method is presented, which can improve the shape quality of the decomposed polygons.(1) Find the Concave Points. Since convex decomposition of concave polygon is to eliminate the concave points, the decomposing line is always drawn from concave points. So the concave points of the concave polygon should be found. The concave point is defined as in Definition 3. The most common method for determining whether a point is concave is to use the vector cross product. Choose the successive points p
i
-
1, p
i, and p
i
+
1 from the point sequence of a polygon in counterclockwise direction. Obtain the vectors p
i
-
1
p
i
→ and p
i
p
i
+
1
→, and then compute the cross product(21)
M
=
x
i
-
x
i
-
1
×
y
i
+
1
-
y
i
-
x
i
+
1
-
x
i
×
y
i
-
y
i
-
1
,where x
i and y
i are the coordinates of point p
i. If M is negative, then p
i is a concave point; otherwise, p
i is a convex point.(2) Find the Visible Point Set. Point q
∈
Q is visible from p
∈
Q if [
p
,
q
]
⊂
Q, where [
p
,
q
] represents the closed path segment [
p
,
q
]
=
p
+
λ
(
q
-
p
)
∣
λ
∈
[
0,1
]. The visibility set is the set of points in Q visible from p. We connect P
i with other vertices of the polygon in counterclockwise direction. For a given vertex P
j, if all the rest vertices behind P
j are on the left side of line P
i
P
j
¯, then P
j is visible from point P
i, as shown in Figure 3.Figure 3
Visible points fromP
1.(3) Find the Decomposing Line Based on the Weighting Function. The decomposing line is obtained from the concave point to the visible point which is carefully searched from the visible point set of this concave point by using a weight function. In [28], a weight function is defined as f
(
α
,
β
)
=
α
-
β, where α and β are the angles between the candidate decomposing line and the neighbor lines. Since α
,
β lie within the interval [
0
,
π
] and the cosine function is monotonous in this interval, we use the following weight function for computational convenience:(22)
f
α
,
β
=
cos
⁡
α
-
cos
⁡
β
.The visible point, from which the decomposing line has the minimum weight value, is chosen. For each polygon, repeat the above steps until all the points of the decomposed polygons are convex.Further, according to the decomposing lines and the edges of the decomposed obstacles, the whole search region is divided into a connected partitionR
~
=
R
~
1
,
…
,
R
~
K of Q. Each partition contains no more than one decomposed convex obstacle. The global decomposition method of the whole region is as follows.First, extend the decomposing lines from two endpoints to intersect with region boundaries, defined as the dividing lines (for the whole region decomposition), while giving up extending those endpoints that will lead to intersecting with obstacles or other extending lines. Then, from each of those abandoned endpoints, draw a radial line along the two edges of the original concave obstacle to intersect with the region boundary or the obstacle boundary, and obtain the dividing lines. By these dividing lines, the whole region is divided into several subregions, each of which containing no more than one decomposed obstacle. Note that the feasible region of a new subregion may be separated into two parts by the decomposed obstacle. In this case, the subregion is divided again to generate a new small subregion without obstacles and the larger one with the obstacle.


### 4.1.2. Extension of Diffeomorphism

Through convex decomposition, the whole region is divided into several subregions with no more than one convex obstacle. The method for constructing the diffeomorphism in Lemma6 requires that ∂
H
∩
∂
R
~
=
∅. However, the new decomposed obstacles may have the common edges with the subregions such that ∂
H
i
∩
∂
R
~
i
≠
∅. So the method in Lemma 6 is not appropriate. The expression of the “common edge” means that the edge of the subregion contains the edge of the obstacle. In the following, we will present a generalized method to construct the diffeomorphism for the new decomposed subregions.LetR
~
=
R
~
1
,
…
,
R
~
K be a compact, convex, and simply connected partition of Q, representing the subregions. For each R
~
i, let H
i
⊂
R
~
i be a convex open set, representing the decomposed obstacle. Thus, the feasible region is R
i
=
R
~
i
∖
H
i. For this problem, the edges generated by the decomposing lines must be considered. Let ∂
R
i
i
n
t
=
∂
R
i
∖
(
∂
R
i
∩
∂
Q
) be the boundary ∂
R
i but not belong to ∂
Q. Similarly, let ∂
R
~
i
i
n
t
=
∂
R
~
i
∖
(
∂
R
~
i
∩
∂
Q
) be the boundary ∂
R
~
i but not belong to ∂
Q.Proposition 7.
Givenp
i
′
∈
H
i, there is a.e. C
1 diffeomorphism φ
i between R
i
i
n and R
~
i
i
n
∖
{
p
′
i
}, where R
i
i
n and R
~
i
i
n are the simplified representations of R
i
i
n
t
∪
∂
R
i
i
n
t and R
~
i
i
n
t
∪
∂
R
~
i
i
n
t, respectively.Proof.
The situation is categorized according to the number of the common edges of the subregion and the obstacle.
Situation  1 (no common edge). Since ∂
H
i
∩
∂
R
~
i
=
∅, the diffeomorphism can be constructed using the method in the proof of Lemma  1 in [16].
Situation  2 (one common edge in Figure4(a)). Draw a radial line l from p
i
′ to cross H
i and intersect with the boundaries of R
~
i and H
i at the points p
1 and p
0, respectively. Note that if the radial line lies on R
~
i
i
n, p
0 and p
1 are the vertices of H
i and R
~
i on the corresponding boundary, that is, points a
1(a
2) and b
1(b
2). Obviously, it can be obtained that p
1
≠
p
0. Then a simple diffeomorphism φ
i can be constructed as(23)
φ
i
q
=
p
1
-
p
i
′
p
1
-
p
i
′
-
p
0
-
p
i
′
q
-
p
i
′
-
p
0
-
p
i
′
.So any q
∈
R
i
i
n on the radial line l can be mapped to φ
i
(
q
)
∈
R
~
i
i
n
∖
{
p
i
′
}. Due to our hypothesis on R
~
i and R
i, it can be seen that φ
i is continuous and differentiable everywhere in R
i
i
n, and it has a unique inverse φ
i
-
1, which is also continuous and differentiable everywhere in R
~
i
i
n
∖
{
p
i
′
}. Therefore, φ
i is a C
1 diffeomorphism. Note that the intersection points p
1 and p
0 are the functions of q. So even if q
1 and q
2 are equidistant to p
i
′, we still have φ
i
(
q
1
)
≠
φ
i
(
q
2
). Hence, the inverse exists and φ
i is a diffeomorphism.
Situation  3 (two common edges in Figure4(b)). The diffeomorphism is constructed as (23) in Situation  2.
Situation  4 (three common edges). There are two cases.
(i) For the first case in Figure4(d), we first define the critical vertex as the common vertex of the subregion and obstacle, while one of its adjacent vertices does not belong to the obstacle, such as p
i
′ and p
2. Choose one of the critical vertices as p
i
′, then p
′
′ is on the same edge with the other critical vertex, and it does not belong to the obstacle. Thus the convex region is divided into two subregions by the blue dividing line connecting p
i
′ with p
′
′. The first diffeomorphism φ
i
a
(
q
) can be constructed by drawing a radial line from p
i
′ using (23) as in Situation  2. For the other part, the diffeomorphism should be modified a little by drawing a radial line from p
′
′, as follows:(24)
φ
i
b
q
=
p
1
-
p
′
′
p
0
-
p
′
′
q
-
p
′
′
,where p
1 and p
0 are the two intersection points of the radial line with the boundaries of the obstacle, and if q lies on the edge ∂
R
~
i
i
n, p
0 and p
1 are the corresponding vertices of the obstacle.
Note that, through the diffeomorphismφ
i
b
(
q
), all the points (except p
i
′ and p
′
′) on the dividing line have the same mapping with the first diffeomorphism φ
i
a
(
q
). And p
′
′ has a mapping by the diffeomorphism φ
i
a. So we obtain φ
i
=
φ
i
a
∪
φ
i
b that can map R
i
i
n to R
~
i
i
n
∖
p
i
′.
(ii) For the second case in Figure4(e), two edges of H
i are parallel which is a typical situation in structural environment. First, the convex polygon is divided into two partitions such that only one partition contains the obstacle. For the partition without any obstacle, the diffeomorphism is constructed as φ
i
a
(
q
)
=
q. For the other partition, a parallel line is drawn from q and intersects with the boundaries of R
~
i at the left point p
0 and right point p
3. Given the vertices of the polygon, p
0 lying on the line b
1
c
1
¯ can be mapped to p
1 on the line p
i
′
c
1
¯, described as p
1
=
p
i
′
-
(
p
i
′
-
c
1
b
1
-
p
0
)
/
(
b
1
-
c
1
). Draw another parallel line from p
1, and it intersects with the boundary of R
~
i at point p
2. Thus the diffeomorphism can be constructed as φ
i
b
(
q
)
=
p
1
+
(
q
-
p
0
p
2
-
p
1
)
/
(
p
3
-
p
0
). Combining with the representation of p
1, the diffeomorphism can be described as(25)
φ
i
b
q
=
p
i
′
-
p
i
′
-
c
1
b
1
-
p
0
b
1
-
c
1
+
q
-
p
0
p
2
-
p
1
p
3
-
p
0
.For this case, φ
i
=
φ
i
a
∪
φ
i
b.
Situation  5 (more than three common edges in Figure4(e)). The diffeomorphism is constructed as in Situation  4.
Situation  6 (a special case). All edges of the obstacle are common edges in Figure 4(f). The diffeomorphism is constructed as (24) by choosing the only vertex of the subregion which does not belong to obstacle as p
′
′.
As claimed in the above six situations, it is obvious that the diffeomorphismφ
i can be constructed for each subregion.Figure 4
Diffeomorphism for the new decomposed regions.
(a)
(b)
(c)
(d)
(e)
(f)Remark 8.
The subregion containing no obstacle will be transformed by an identical transformation, which is also a diffeomorphism.Remark 9.
Let⋃
i
=
1
M
p
i
′
=
p
′. For each decomposed subregion, the diffeomorphism φ
i is constructed. So the transformation φ
:
F
i
n
t
→
Q
i
n
t
/
p
′ is defined by φ
(
q
)
=
φ
i
(
q
), where q
∈
R
i
i
n.Remark 10.
Given the search region and obstacles, the diffeomorphism can be computed off-line and stored in a table. Therefore, in the search procedure, the diffeomorphism can be directly used by looking up the table, which will significantly reduce the computation time.


### 4.1.3. Numerical Examples

Two examples of applying the convex decomposition method to divide the nonconvex region with arbitrary obstacles are shown in Figure5. The green lines are the decomposing lines for convex decomposition of the concave obstacle, and the black lines are the dividing lines for the whole environment. In Figure 5(a), the concave obstacle is divided into three convex polygons. And the whole region is divided into four subregions, three of which containing one obstacle and one containing none. In Figure 5(b), a more complex nonconvex region with two obstacles is given. First, the whole region is divided into two partitions. Then the partition with a concave obstacle is divided into six subregions by using the decomposition method.Figure 5
Convex decomposition.
(a)
Nonconvex regionQ
1 with a typical concave obstacle
(b)
Nonconvex regionQ
2 with a convex obstacle and a complex concave obstacleFurther, the diffeomorphism is constructed for each example region by using the method in Proposition7. We sample some points in the original regions and then map them to the transformed regions by the diffeomorphism, as shown in Figure 6.Figure 6
Diffeomorphism of the nonconvex region. (a) and (c) are the two examples of the original nonconvex region with 10 randomly sampled points. (b) and (d) are the transformed convex regions with the transformed image of the sampled points by diffeomorphism.
(a)
Sampled points in originalQ
1
(b)
Point images in transformedφ
(
Q
1
)
(c)
Sampled points in originalQ
2
(d)
Point images in transformedφ
(
Q
2
)


## 4.1.1. Global Convex Decomposition

As the main preprocessing procedure for constructing the transformation, the global convex decomposition method decomposes the original nonconvex region and the concave obstacles, making it possible to find the diffeomorphism based transformation in the nonconvex environment with concave obstacles.First, the convex decomposition of concave obstacles is given. Some works have been conducted on the convex decomposition of concave polygons. Schachter [25] utilizes the Delaunay method to decompose polygons into convex sets. Chazelle and Dobkin [26] construct the concave point sets of the x
k type to realize the convex decomposition. Keil [27] uses improved dynamic programming algorithms to achieve decomposing a simple polygon. In this paper, our purpose for convex decomposition of concave polygons is to find the diffeomorphism, which is related to the visibility. So a point visibility based convex decomposition method is presented, which can improve the shape quality of the decomposed polygons.(1) Find the Concave Points. Since convex decomposition of concave polygon is to eliminate the concave points, the decomposing line is always drawn from concave points. So the concave points of the concave polygon should be found. The concave point is defined as in Definition 3. The most common method for determining whether a point is concave is to use the vector cross product. Choose the successive points p
i
-
1, p
i, and p
i
+
1 from the point sequence of a polygon in counterclockwise direction. Obtain the vectors p
i
-
1
p
i
→ and p
i
p
i
+
1
→, and then compute the cross product(21)
M
=
x
i
-
x
i
-
1
×
y
i
+
1
-
y
i
-
x
i
+
1
-
x
i
×
y
i
-
y
i
-
1
,where x
i and y
i are the coordinates of point p
i. If M is negative, then p
i is a concave point; otherwise, p
i is a convex point.(2) Find the Visible Point Set. Point q
∈
Q is visible from p
∈
Q if [
p
,
q
]
⊂
Q, where [
p
,
q
] represents the closed path segment [
p
,
q
]
=
p
+
λ
(
q
-
p
)
∣
λ
∈
[
0,1
]. The visibility set is the set of points in Q visible from p. We connect P
i with other vertices of the polygon in counterclockwise direction. For a given vertex P
j, if all the rest vertices behind P
j are on the left side of line P
i
P
j
¯, then P
j is visible from point P
i, as shown in Figure 3.Figure 3
Visible points fromP
1.(3) Find the Decomposing Line Based on the Weighting Function. The decomposing line is obtained from the concave point to the visible point which is carefully searched from the visible point set of this concave point by using a weight function. In [28], a weight function is defined as f
(
α
,
β
)
=
α
-
β, where α and β are the angles between the candidate decomposing line and the neighbor lines. Since α
,
β lie within the interval [
0
,
π
] and the cosine function is monotonous in this interval, we use the following weight function for computational convenience:(22)
f
α
,
β
=
cos
⁡
α
-
cos
⁡
β
.The visible point, from which the decomposing line has the minimum weight value, is chosen. For each polygon, repeat the above steps until all the points of the decomposed polygons are convex.Further, according to the decomposing lines and the edges of the decomposed obstacles, the whole search region is divided into a connected partitionR
~
=
R
~
1
,
…
,
R
~
K of Q. Each partition contains no more than one decomposed convex obstacle. The global decomposition method of the whole region is as follows.First, extend the decomposing lines from two endpoints to intersect with region boundaries, defined as the dividing lines (for the whole region decomposition), while giving up extending those endpoints that will lead to intersecting with obstacles or other extending lines. Then, from each of those abandoned endpoints, draw a radial line along the two edges of the original concave obstacle to intersect with the region boundary or the obstacle boundary, and obtain the dividing lines. By these dividing lines, the whole region is divided into several subregions, each of which containing no more than one decomposed obstacle. Note that the feasible region of a new subregion may be separated into two parts by the decomposed obstacle. In this case, the subregion is divided again to generate a new small subregion without obstacles and the larger one with the obstacle.


## 4.1.2. Extension of Diffeomorphism

Through convex decomposition, the whole region is divided into several subregions with no more than one convex obstacle. The method for constructing the diffeomorphism in Lemma6 requires that ∂
H
∩
∂
R
~
=
∅. However, the new decomposed obstacles may have the common edges with the subregions such that ∂
H
i
∩
∂
R
~
i
≠
∅. So the method in Lemma 6 is not appropriate. The expression of the “common edge” means that the edge of the subregion contains the edge of the obstacle. In the following, we will present a generalized method to construct the diffeomorphism for the new decomposed subregions.LetR
~
=
R
~
1
,
…
,
R
~
K be a compact, convex, and simply connected partition of Q, representing the subregions. For each R
~
i, let H
i
⊂
R
~
i be a convex open set, representing the decomposed obstacle. Thus, the feasible region is R
i
=
R
~
i
∖
H
i. For this problem, the edges generated by the decomposing lines must be considered. Let ∂
R
i
i
n
t
=
∂
R
i
∖
(
∂
R
i
∩
∂
Q
) be the boundary ∂
R
i but not belong to ∂
Q. Similarly, let ∂
R
~
i
i
n
t
=
∂
R
~
i
∖
(
∂
R
~
i
∩
∂
Q
) be the boundary ∂
R
~
i but not belong to ∂
Q.Proposition 7.
Givenp
i
′
∈
H
i, there is a.e. C
1 diffeomorphism φ
i between R
i
i
n and R
~
i
i
n
∖
{
p
′
i
}, where R
i
i
n and R
~
i
i
n are the simplified representations of R
i
i
n
t
∪
∂
R
i
i
n
t and R
~
i
i
n
t
∪
∂
R
~
i
i
n
t, respectively.Proof.
The situation is categorized according to the number of the common edges of the subregion and the obstacle.
Situation  1 (no common edge). Since ∂
H
i
∩
∂
R
~
i
=
∅, the diffeomorphism can be constructed using the method in the proof of Lemma  1 in [16].
Situation  2 (one common edge in Figure4(a)). Draw a radial line l from p
i
′ to cross H
i and intersect with the boundaries of R
~
i and H
i at the points p
1 and p
0, respectively. Note that if the radial line lies on R
~
i
i
n, p
0 and p
1 are the vertices of H
i and R
~
i on the corresponding boundary, that is, points a
1(a
2) and b
1(b
2). Obviously, it can be obtained that p
1
≠
p
0. Then a simple diffeomorphism φ
i can be constructed as(23)
φ
i
q
=
p
1
-
p
i
′
p
1
-
p
i
′
-
p
0
-
p
i
′
q
-
p
i
′
-
p
0
-
p
i
′
.So any q
∈
R
i
i
n on the radial line l can be mapped to φ
i
(
q
)
∈
R
~
i
i
n
∖
{
p
i
′
}. Due to our hypothesis on R
~
i and R
i, it can be seen that φ
i is continuous and differentiable everywhere in R
i
i
n, and it has a unique inverse φ
i
-
1, which is also continuous and differentiable everywhere in R
~
i
i
n
∖
{
p
i
′
}. Therefore, φ
i is a C
1 diffeomorphism. Note that the intersection points p
1 and p
0 are the functions of q. So even if q
1 and q
2 are equidistant to p
i
′, we still have φ
i
(
q
1
)
≠
φ
i
(
q
2
). Hence, the inverse exists and φ
i is a diffeomorphism.
Situation  3 (two common edges in Figure4(b)). The diffeomorphism is constructed as (23) in Situation  2.
Situation  4 (three common edges). There are two cases.
(i) For the first case in Figure4(d), we first define the critical vertex as the common vertex of the subregion and obstacle, while one of its adjacent vertices does not belong to the obstacle, such as p
i
′ and p
2. Choose one of the critical vertices as p
i
′, then p
′
′ is on the same edge with the other critical vertex, and it does not belong to the obstacle. Thus the convex region is divided into two subregions by the blue dividing line connecting p
i
′ with p
′
′. The first diffeomorphism φ
i
a
(
q
) can be constructed by drawing a radial line from p
i
′ using (23) as in Situation  2. For the other part, the diffeomorphism should be modified a little by drawing a radial line from p
′
′, as follows:(24)
φ
i
b
q
=
p
1
-
p
′
′
p
0
-
p
′
′
q
-
p
′
′
,where p
1 and p
0 are the two intersection points of the radial line with the boundaries of the obstacle, and if q lies on the edge ∂
R
~
i
i
n, p
0 and p
1 are the corresponding vertices of the obstacle.
Note that, through the diffeomorphismφ
i
b
(
q
), all the points (except p
i
′ and p
′
′) on the dividing line have the same mapping with the first diffeomorphism φ
i
a
(
q
). And p
′
′ has a mapping by the diffeomorphism φ
i
a. So we obtain φ
i
=
φ
i
a
∪
φ
i
b that can map R
i
i
n to R
~
i
i
n
∖
p
i
′.
(ii) For the second case in Figure4(e), two edges of H
i are parallel which is a typical situation in structural environment. First, the convex polygon is divided into two partitions such that only one partition contains the obstacle. For the partition without any obstacle, the diffeomorphism is constructed as φ
i
a
(
q
)
=
q. For the other partition, a parallel line is drawn from q and intersects with the boundaries of R
~
i at the left point p
0 and right point p
3. Given the vertices of the polygon, p
0 lying on the line b
1
c
1
¯ can be mapped to p
1 on the line p
i
′
c
1
¯, described as p
1
=
p
i
′
-
(
p
i
′
-
c
1
b
1
-
p
0
)
/
(
b
1
-
c
1
). Draw another parallel line from p
1, and it intersects with the boundary of R
~
i at point p
2. Thus the diffeomorphism can be constructed as φ
i
b
(
q
)
=
p
1
+
(
q
-
p
0
p
2
-
p
1
)
/
(
p
3
-
p
0
). Combining with the representation of p
1, the diffeomorphism can be described as(25)
φ
i
b
q
=
p
i
′
-
p
i
′
-
c
1
b
1
-
p
0
b
1
-
c
1
+
q
-
p
0
p
2
-
p
1
p
3
-
p
0
.For this case, φ
i
=
φ
i
a
∪
φ
i
b.
Situation  5 (more than three common edges in Figure4(e)). The diffeomorphism is constructed as in Situation  4.
Situation  6 (a special case). All edges of the obstacle are common edges in Figure 4(f). The diffeomorphism is constructed as (24) by choosing the only vertex of the subregion which does not belong to obstacle as p
′
′.
As claimed in the above six situations, it is obvious that the diffeomorphismφ
i can be constructed for each subregion.Figure 4
Diffeomorphism for the new decomposed regions.
(a)
(b)
(c)
(d)
(e)
(f)Remark 8.
The subregion containing no obstacle will be transformed by an identical transformation, which is also a diffeomorphism.Remark 9.
Let⋃
i
=
1
M
p
i
′
=
p
′. For each decomposed subregion, the diffeomorphism φ
i is constructed. So the transformation φ
:
F
i
n
t
→
Q
i
n
t
/
p
′ is defined by φ
(
q
)
=
φ
i
(
q
), where q
∈
R
i
i
n.Remark 10.
Given the search region and obstacles, the diffeomorphism can be computed off-line and stored in a table. Therefore, in the search procedure, the diffeomorphism can be directly used by looking up the table, which will significantly reduce the computation time.


## 4.1.3. Numerical Examples

Two examples of applying the convex decomposition method to divide the nonconvex region with arbitrary obstacles are shown in Figure5. The green lines are the decomposing lines for convex decomposition of the concave obstacle, and the black lines are the dividing lines for the whole environment. In Figure 5(a), the concave obstacle is divided into three convex polygons. And the whole region is divided into four subregions, three of which containing one obstacle and one containing none. In Figure 5(b), a more complex nonconvex region with two obstacles is given. First, the whole region is divided into two partitions. Then the partition with a concave obstacle is divided into six subregions by using the decomposition method.Figure 5
Convex decomposition.
(a)
Nonconvex regionQ
1 with a typical concave obstacle
(b)
Nonconvex regionQ
2 with a convex obstacle and a complex concave obstacleFurther, the diffeomorphism is constructed for each example region by using the method in Proposition7. We sample some points in the original regions and then map them to the transformed regions by the diffeomorphism, as shown in Figure 6.Figure 6
Diffeomorphism of the nonconvex region. (a) and (c) are the two examples of the original nonconvex region with 10 randomly sampled points. (b) and (d) are the transformed convex regions with the transformed image of the sampled points by diffeomorphism.
(a)
Sampled points in originalQ
1
(b)
Point images in transformedφ
(
Q
1
)
(c)
Sampled points in originalQ
2
(d)
Point images in transformedφ
(
Q
2
)


## 4.2. Cooperative Search Method in Nonconvex Environment with Arbitrary Obstacles

The diffeomorphism that transforms a complex nonconvex region with arbitrary obstacles to a convex region has been constructed. Now we will propose the cooperative search method in the nonconvex environment based on the transformation.


### 4.2.1. Objective Function Transformation

The objective function (3) should be modified by shifting the allowable search domain from the feasible region F
i
n
t to the whole region Q
i
n
t
∖
{
p
′
}. Based on the modification of Variables Theorem [29], it claims that if Ω
⊂
R
n is open, g
1
:
Ω
→
R is integrable over Ω, and g
2
:
Ω
~
→
Ω is a C
1 diffeomorphism, then g
1
∘
g
2
det
⁡
J
g
2 is integrable over Ω
~, and(26)
∫
Ω
g
1
ω
d
ω
=
∫
Ω
~
g
1
∘
g
2
ω
~
det
⁡
J
g
2
ω
~
d
ω
~
,where det
⁡
J
g
2
ω
~ is the determinant of the Jacobian of g
2. In this paper, we assume that g
1
=
f
(
p
i
-
r
)
γ
(
r
)
:
F
i
n
t
→
R and g
2
=
φ
-
1
:
Q
i
n
t
/
(
p
′
)
→
R
i
n
t. Since φ
-
1
(
ω
)
=
r, the transformed objective function in the new convex region can be expressed as(27)
∫
F
i
n
t
f
p
i
-
r
γ
r
d
r
=
∫
Q
i
n
t
∖
p
′
f
φ
-
1
q
i
-
φ
-
1
ω
γ
φ
-
1
ω
det
⁡
J
φ
-
1
ω
d
ω
.It can be seen that the allowable search domain is transformed from the feasible regionF
i
n
t to the whole region Q
i
n
t
∖
p
′, while keeping the function value unchanged. However, the control law for the above objective function (the right side of (27)) is quite hard to be obtained. In order to use the Lloyd algorithm, the new objective function is defined as(28)
H
′
P
=
∫
F
i
n
t
f
φ
p
i
-
φ
r
γ
r
d
r
.Sinceφ is the diffeomorphism, F
i
n
t and φ
(
F
i
n
t
) are topologically equivalent, and the problem remains essentially unchanged. Further, the corresponding function in the transformed domain is obtained:(29)
H
W
,
V
~
=
∫
Q
i
n
t
∖
p
′
f
q
i
-
ω
γ
φ
-
1
ω
det
⁡
J
φ
-
1
ω
d
ω
.So the gradient descent approach as (6) can be used to minimize the objective function(30)
q
˙
i
=
-
k
i
1
q
i
-
C
~
V
i
,where k
i
2 is a positive gain, q
i
=
φ
(
p
i
) is the transformed position of p
i by the diffeomorphism in the transformed region Q
i
n
t
/
(
p
′
), and C
~
V
i is the centroid of the Voronoi region generated by q
i.


### 4.2.2. Control Strategy for Cooperative Search in the Nonconvex Environment

Since the diffeomorphismφ is a bijection, the trajectory of q
i in Q
i
n
t
/
p
′ has the unique preimage in original region F
i
n
t, which is the trajectory of p
i. So the control law for the new objective function (28) in the original region can be obtained by the inverse mapping of (30):(31)
p
i
˙
=
-
k
i
2
J
φ
-
1
φ
p
i
p
i
-
φ
-
1
C
~
V
i
,where k
i
2 is a positive gain and J
φ is the Jacobian of φ.The following is an analysis of the convergence of the control law (31) using LaSalle’s invariance principle [30], which is basically an extension of Lyapunov’s theorem [31].Theorem 11.
In the nonconvex region, the trajectories of UAVs governed by the control law (31), starting from any initial condition P
(
0
)
=
(
p
1,0
,
…
,
p
N
,
0
)
∈
F
i
n
t, will asymptotically converge to the desired location.Proof.
Since the diffeomorphismφ is a bijection, the trajectory q
i in Q
i
n
t
∖
p
′ has the unique preimage p
i in original region F
i
n
t. Therefore, to prove that the trajectory p
i governed by the control law (31), starting from any initial condition P
(
0
), can converge to φ
-
1
(
C
~
V
i
), we only need to prove that the trajectory q
i governed by the control law (30), starting from initial condition W
(
0
)
=
φ
(
P
(
0
)
)
=
(
φ
(
p
1,0
)
,
…
,
φ
(
p
N
,
0
)
), can converge to C
~
V
i.
AssumeV
(
Q
)
=
H
(
W
,
V
~
)
=
∑
i
=
1
N
‍
∫
V
~
i
‍
q
i
-
ω
2
γ
φ
-
1
ω
det
⁡
J
φ
-
1
ω
d
ω, where Q
=
(
q
1
,
q
2
,
…
,
q
N
) represent the transformed positions of UAVs in Q
i
n
t
∖
p
′ through the diffeomorphism φ. Consider(32)
V
˙
Q
=
d
H
W
,
V
~
d
t
=
∑
i
=
1
N
∂
H
W
,
V
~
∂
q
i
q
˙
i
=
∑
i
=
1
N
2
λ
M
~
V
i
q
i
-
C
~
V
i
-
k
i
1
q
i
-
C
~
V
i
=
-
2
λ
1
k
i
1
∑
i
=
1
N
M
~
V
i
q
i
-
C
~
V
i
2
⩽
0
,where λ
1
=
det
⁡
J
φ
-
1
ω and M
~
V
i is the mass of the Voronoi region generated by q
i in the convex region Q
i
n
t
∖
p
′.
The following can be observed:(i)
V
:
Q
i
n
t
∖
p
′
→
R is continuously differentiable in Q
i
n
t
∖
p
′.
(ii)
Q
i
n
t
∖
p
′ is a compact invariant set.
(iii)
V
˙ is negative semidefinite in Q
i
n
t
∖
p
′.
(iv)
E
=
V
˙
-
1
(
0
)
=
C
~
V
i.
(v)
E itself is the largest invariant subset of E by the control law (30).
Thus, by LaSalle’s invariance principle, the trajectoriesq
i governed by (30), starting from initial configuration W
(
0
), will asymptotically converge to the critical points of H
(
W
,
V
~
); that is, the trajectories of UAVs governed by the control law (31), starting from any initial condition P
(
0
)
=
(
p
1,0
,
…
,
p
N
,
0
)
∈
F
i
n
t, will asymptotically converge to the desired location.Therefore, the discrete-time motion model and the real control law for UAVi in the original nonconvex region can be written as(33)
p
i
,
m
=
p
i
,
m
-
1
+
u
i
,
m
,
u
i
,
m
=
-
k
i
2
J
φ
-
1
φ
p
i
,
m
-
1
p
i
,
m
-
1
-
φ
-
1
C
~
V
i
.Further, considering the mobility of the UAV, the control law must meet some constraints such that the generated trajectories are kinematically feasible:(1) Maximum velocity: the maximum velocityv
m
a
x will limit the control law by the following saturation rule:(34)
u
i
,
m
=
u
i
,
m
,
u
i
,
m
≤
v
max
⁡
v
max
⁡
u
i
,
m
u
i
,
m
,
Otherwise
.(2) Maximum turning angle: the maximum turning angle during a unit time periodΔ
t is derived from the turning rate ω
m
a
x, described as θ
m
a
x
=
ω
m
a
x
Δ
t. Thus the control law is limited by the following equations:(35)
u
i
,
m
=
u
i
,
m
,
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
≤
θ
max
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
-
1
A
,
Otherwise
,where A is the rotating matrix for the maximum turning angle from the last control input as (36)
A
=
cos
⁡
α
-
sin
⁡
α
sin
⁡
α
cos
⁡
α
,
α
=
θ
max
⁡
,
u
i
,
m
·
u
i
,
m
-
1
≥
0
-
θ
max
⁡
,
u
i
,
m
·
u
i
,
m
-
1
<
0represents the rotating direction.Combining (34) with (35), the control law can be expressed as(37)
u
i
,
m
=
u
i
,
m
,
u
i
,
m
≤
v
max
⁡
and
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
≤
θ
max
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
-
1
A
,
u
i
,
m
≤
v
max
⁡
and
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
>
θ
max
⁡
v
max
⁡
u
i
,
m
u
i
,
m
,
u
i
,
m
>
v
max
⁡
and
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
≤
θ
max
⁡
v
max
⁡
u
i
,
m
-
1
u
i
,
m
-
1
A
,
Otherwise
,where u
i
,
m
=
-
k
i
2
J
φ
-
1
(
φ
(
p
i
,
m
-
1
)
)
(
p
i
,
m
-
1
-
φ
-
1
(
C
~
V
i
)
), as shown in (33).


### 4.2.3. Cooperative Search Algorithm in the Nonconvex Environment with Arbitrary Obstacles

The cooperative search algorithm for the nonconvex region with arbitrary obstacles is presented as Algorithm1.Algorithm 1.
Cooperative search algorithm in the nonconvex environment with arbitrary obstacles is as follows:(1)
Initially, all the UAVs are randomly deployed in the feasible region. The probability map and uncertainty map of each UAV are initialized.
(2)
The nonconvex region with arbitrary obstacles is divided by the global decomposition method, and the diffeomorphismφ is constructed as described in Proposition 7.
(3)
Through the diffeomorphismφ, the original nonconvex region F
i
n
t is transformed into the new convex one, Q
i
n
t
∖
p
′. And the positions of UAVs are also transformed into new virtual ones.
(4)
The Voronoi partition of transformed regionQ
i
n
t
∖
p
′ is generated by the virtual positions of UAVs. Compute the centroid C
~
V
i of each Voronoi region using probability map as the density.
(5)
At each time step, UAVi moves towards the preimages of its generalized Voronoi centroid φ
-
1
(
C
~
V
i
) in the original region, governed by the control law (31).
(6)
Update the probability mapT
P
i of UAV i. If the target probability of a cell is greater than P
m
a
x, a target is found.
(7)
Update the Voronoi partitions based on the new virtual positions of UAVs. Compute the generalized Voronoi centroidsC
~
V
i.
(8)
If UAVi reaches the preimage of its Voronoi centroid φ
-
1
(
C
~
V
i
) and does not move, then use its uncertainty map instead of probability map to compute the new Voronoi centroid φ
-
1
(
C
~
V
i
) from this time step.
(9)
If the terminal condition is satisfied, the algorithm ends; otherwise, go to(
5
).


## 4.2.1. Objective Function Transformation

The objective function (3) should be modified by shifting the allowable search domain from the feasible region F
i
n
t to the whole region Q
i
n
t
∖
{
p
′
}. Based on the modification of Variables Theorem [29], it claims that if Ω
⊂
R
n is open, g
1
:
Ω
→
R is integrable over Ω, and g
2
:
Ω
~
→
Ω is a C
1 diffeomorphism, then g
1
∘
g
2
det
⁡
J
g
2 is integrable over Ω
~, and(26)
∫
Ω
g
1
ω
d
ω
=
∫
Ω
~
g
1
∘
g
2
ω
~
det
⁡
J
g
2
ω
~
d
ω
~
,where det
⁡
J
g
2
ω
~ is the determinant of the Jacobian of g
2. In this paper, we assume that g
1
=
f
(
p
i
-
r
)
γ
(
r
)
:
F
i
n
t
→
R and g
2
=
φ
-
1
:
Q
i
n
t
/
(
p
′
)
→
R
i
n
t. Since φ
-
1
(
ω
)
=
r, the transformed objective function in the new convex region can be expressed as(27)
∫
F
i
n
t
f
p
i
-
r
γ
r
d
r
=
∫
Q
i
n
t
∖
p
′
f
φ
-
1
q
i
-
φ
-
1
ω
γ
φ
-
1
ω
det
⁡
J
φ
-
1
ω
d
ω
.It can be seen that the allowable search domain is transformed from the feasible regionF
i
n
t to the whole region Q
i
n
t
∖
p
′, while keeping the function value unchanged. However, the control law for the above objective function (the right side of (27)) is quite hard to be obtained. In order to use the Lloyd algorithm, the new objective function is defined as(28)
H
′
P
=
∫
F
i
n
t
f
φ
p
i
-
φ
r
γ
r
d
r
.Sinceφ is the diffeomorphism, F
i
n
t and φ
(
F
i
n
t
) are topologically equivalent, and the problem remains essentially unchanged. Further, the corresponding function in the transformed domain is obtained:(29)
H
W
,
V
~
=
∫
Q
i
n
t
∖
p
′
f
q
i
-
ω
γ
φ
-
1
ω
det
⁡
J
φ
-
1
ω
d
ω
.So the gradient descent approach as (6) can be used to minimize the objective function(30)
q
˙
i
=
-
k
i
1
q
i
-
C
~
V
i
,where k
i
2 is a positive gain, q
i
=
φ
(
p
i
) is the transformed position of p
i by the diffeomorphism in the transformed region Q
i
n
t
/
(
p
′
), and C
~
V
i is the centroid of the Voronoi region generated by q
i.


## 4.2.2. Control Strategy for Cooperative Search in the Nonconvex Environment

Since the diffeomorphismφ is a bijection, the trajectory of q
i in Q
i
n
t
/
p
′ has the unique preimage in original region F
i
n
t, which is the trajectory of p
i. So the control law for the new objective function (28) in the original region can be obtained by the inverse mapping of (30):(31)
p
i
˙
=
-
k
i
2
J
φ
-
1
φ
p
i
p
i
-
φ
-
1
C
~
V
i
,where k
i
2 is a positive gain and J
φ is the Jacobian of φ.The following is an analysis of the convergence of the control law (31) using LaSalle’s invariance principle [30], which is basically an extension of Lyapunov’s theorem [31].Theorem 11.
In the nonconvex region, the trajectories of UAVs governed by the control law (31), starting from any initial condition P
(
0
)
=
(
p
1,0
,
…
,
p
N
,
0
)
∈
F
i
n
t, will asymptotically converge to the desired location.Proof.
Since the diffeomorphismφ is a bijection, the trajectory q
i in Q
i
n
t
∖
p
′ has the unique preimage p
i in original region F
i
n
t. Therefore, to prove that the trajectory p
i governed by the control law (31), starting from any initial condition P
(
0
), can converge to φ
-
1
(
C
~
V
i
), we only need to prove that the trajectory q
i governed by the control law (30), starting from initial condition W
(
0
)
=
φ
(
P
(
0
)
)
=
(
φ
(
p
1,0
)
,
…
,
φ
(
p
N
,
0
)
), can converge to C
~
V
i.
AssumeV
(
Q
)
=
H
(
W
,
V
~
)
=
∑
i
=
1
N
‍
∫
V
~
i
‍
q
i
-
ω
2
γ
φ
-
1
ω
det
⁡
J
φ
-
1
ω
d
ω, where Q
=
(
q
1
,
q
2
,
…
,
q
N
) represent the transformed positions of UAVs in Q
i
n
t
∖
p
′ through the diffeomorphism φ. Consider(32)
V
˙
Q
=
d
H
W
,
V
~
d
t
=
∑
i
=
1
N
∂
H
W
,
V
~
∂
q
i
q
˙
i
=
∑
i
=
1
N
2
λ
M
~
V
i
q
i
-
C
~
V
i
-
k
i
1
q
i
-
C
~
V
i
=
-
2
λ
1
k
i
1
∑
i
=
1
N
M
~
V
i
q
i
-
C
~
V
i
2
⩽
0
,where λ
1
=
det
⁡
J
φ
-
1
ω and M
~
V
i is the mass of the Voronoi region generated by q
i in the convex region Q
i
n
t
∖
p
′.
The following can be observed:(i)
V
:
Q
i
n
t
∖
p
′
→
R is continuously differentiable in Q
i
n
t
∖
p
′.
(ii)
Q
i
n
t
∖
p
′ is a compact invariant set.
(iii)
V
˙ is negative semidefinite in Q
i
n
t
∖
p
′.
(iv)
E
=
V
˙
-
1
(
0
)
=
C
~
V
i.
(v)
E itself is the largest invariant subset of E by the control law (30).
Thus, by LaSalle’s invariance principle, the trajectoriesq
i governed by (30), starting from initial configuration W
(
0
), will asymptotically converge to the critical points of H
(
W
,
V
~
); that is, the trajectories of UAVs governed by the control law (31), starting from any initial condition P
(
0
)
=
(
p
1,0
,
…
,
p
N
,
0
)
∈
F
i
n
t, will asymptotically converge to the desired location.Therefore, the discrete-time motion model and the real control law for UAVi in the original nonconvex region can be written as(33)
p
i
,
m
=
p
i
,
m
-
1
+
u
i
,
m
,
u
i
,
m
=
-
k
i
2
J
φ
-
1
φ
p
i
,
m
-
1
p
i
,
m
-
1
-
φ
-
1
C
~
V
i
.Further, considering the mobility of the UAV, the control law must meet some constraints such that the generated trajectories are kinematically feasible:(1) Maximum velocity: the maximum velocityv
m
a
x will limit the control law by the following saturation rule:(34)
u
i
,
m
=
u
i
,
m
,
u
i
,
m
≤
v
max
⁡
v
max
⁡
u
i
,
m
u
i
,
m
,
Otherwise
.(2) Maximum turning angle: the maximum turning angle during a unit time periodΔ
t is derived from the turning rate ω
m
a
x, described as θ
m
a
x
=
ω
m
a
x
Δ
t. Thus the control law is limited by the following equations:(35)
u
i
,
m
=
u
i
,
m
,
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
≤
θ
max
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
-
1
A
,
Otherwise
,where A is the rotating matrix for the maximum turning angle from the last control input as (36)
A
=
cos
⁡
α
-
sin
⁡
α
sin
⁡
α
cos
⁡
α
,
α
=
θ
max
⁡
,
u
i
,
m
·
u
i
,
m
-
1
≥
0
-
θ
max
⁡
,
u
i
,
m
·
u
i
,
m
-
1
<
0represents the rotating direction.Combining (34) with (35), the control law can be expressed as(37)
u
i
,
m
=
u
i
,
m
,
u
i
,
m
≤
v
max
⁡
and
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
≤
θ
max
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
-
1
A
,
u
i
,
m
≤
v
max
⁡
and
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
>
θ
max
⁡
v
max
⁡
u
i
,
m
u
i
,
m
,
u
i
,
m
>
v
max
⁡
and
arccos
⁡
u
i
,
m
u
i
,
m
-
1
u
i
,
m
u
i
,
m
-
1
≤
θ
max
⁡
v
max
⁡
u
i
,
m
-
1
u
i
,
m
-
1
A
,
Otherwise
,where u
i
,
m
=
-
k
i
2
J
φ
-
1
(
φ
(
p
i
,
m
-
1
)
)
(
p
i
,
m
-
1
-
φ
-
1
(
C
~
V
i
)
), as shown in (33).


## 4.2.3. Cooperative Search Algorithm in the Nonconvex Environment with Arbitrary Obstacles

The cooperative search algorithm for the nonconvex region with arbitrary obstacles is presented as Algorithm1.Algorithm 1.
Cooperative search algorithm in the nonconvex environment with arbitrary obstacles is as follows:(1)
Initially, all the UAVs are randomly deployed in the feasible region. The probability map and uncertainty map of each UAV are initialized.
(2)
The nonconvex region with arbitrary obstacles is divided by the global decomposition method, and the diffeomorphismφ is constructed as described in Proposition 7.
(3)
Through the diffeomorphismφ, the original nonconvex region F
i
n
t is transformed into the new convex one, Q
i
n
t
∖
p
′. And the positions of UAVs are also transformed into new virtual ones.
(4)
The Voronoi partition of transformed regionQ
i
n
t
∖
p
′ is generated by the virtual positions of UAVs. Compute the centroid C
~
V
i of each Voronoi region using probability map as the density.
(5)
At each time step, UAVi moves towards the preimages of its generalized Voronoi centroid φ
-
1
(
C
~
V
i
) in the original region, governed by the control law (31).
(6)
Update the probability mapT
P
i of UAV i. If the target probability of a cell is greater than P
m
a
x, a target is found.
(7)
Update the Voronoi partitions based on the new virtual positions of UAVs. Compute the generalized Voronoi centroidsC
~
V
i.
(8)
If UAVi reaches the preimage of its Voronoi centroid φ
-
1
(
C
~
V
i
) and does not move, then use its uncertainty map instead of probability map to compute the new Voronoi centroid φ
-
1
(
C
~
V
i
) from this time step.
(9)
If the terminal condition is satisfied, the algorithm ends; otherwise, go to(
5
).


## 5. Simulation

In this section, some experiments are performed to verify the effectiveness of the cooperative search algorithm in different scenarios, and then the comparison with other algorithms is conducted to demonstrate the advantage of our algorithm in addressing the cooperative search problem in a nonconvex environment.


### 5.1. Simulation Environment

The algorithms have been programmed in Matlab. The simulation scenario is set as that five UAVs search for targets in two typical nonconvex regions with differently shaped obstacles. The whole search regionQ is defined as a square region of [
0,10
]
×
[
0,10
] km2. And the distribution of each target follows the Gaussian distribution function(38)
ϕ
r
=
1
σ
2
π
e
-
r
-
μ
2
/
2
σ
2
,where the mean value μ is the expected position of the target and σ is set to be 1 km, representing the predefined probability distribution range. The communication range R
C is 5 km. Initially, all the UAVs are randomly deployed in the feasible region. Some other parameters for simulation are in Table 1.Table 1
Simulation parameters.
UAV
Sensor
Detection model
Maximum speed (m/s)
100
α
s
45°
Positive sensing probabilityP
0
i
0.9
Maximum turning angle
60°
β
s
30°
False alarm probabilityp
a
0.1
Flight height (m)
100
γ
s
30°
Target existence criteriaP
m
a
x
0.9


### 5.2. Simulation Results and Analysis

#### 5.2.1. Scenario  1: Nonconvex Environment with Concave Obstacles

The first obstacle area is a typical concave region, represented as a red square. The convex decomposition and the diffeomorphism have been implemented as shown in Figures5(a), 6(a), and 6(b). Then the cooperative search algorithm is performed. The trajectories of UAVs at different time steps in the transformed region without obstacles and in the original region are shown in Figure 7, where the UAVs and targets are denoted by “o” and “+”, respectively. The probability map with time evolving is shown in Figure 8. From the figures, it can be seen that the probability map has converged. UAVs in the original region, starting from any initial positions, can find all the 6 targets at time 45 s while avoiding the obstacle.Figure 7
Search flight trajectories of UAVs with time: (a–c) in the transformed region and (d–f) in the original region.
(a)
0 s
(b)
30 s
(c)
60 s
(d)
0 s
(e)
30 s
(f)
60 sFigure 8
Probability map update for Scenario  1.
(a)
0 s
(b)
30 s
(c)
60 s


#### 5.2.2. Scenario  2: Nonconvex Environment with Several Complex Obstacles

The nonconvex region in Scenario  2 contains two obstacles. The results of cooperative search are shown in Figures9 and 10. It can be seen that the trajectories of UAVs in the original region governed by the control law (31), starting from any initial positions, can move towards the targets while avoiding the obstacle. It also can be seen that the probability map has converged. After finding the targets, the UAV uses the uncertainty map instead of the probability map to escape from the local optimum (the positions of targets with high target existence probability) and keep searching. From the simulation results, all the targets are found at time 54 s.Figure 9
Search flight trajectories of UAVs with time: (a–c) in the transformed region and (d–f) in the original region.
(a)
0 s
(b)
30 s
(c)
60 s
(d)
0 s
(e)
30 s
(f)
60 sFigure 10
Probability map update for Scenario  2.
(a)
0 s
(b)
30 s
(c)
60 s


#### 5.2.3. Comparisons with Other Algorithms

In order to verify the performance, the proposed algorithm is compared with a sweep-line search algorithm [7] and a random search algorithm [32]. In the sweep-line search algorithm, the search mood is fixed, and each UAV sweeps the region parallel to the boundary and decides the width between two search lines based on the detection range. In the random search algorithm, each UAV randomly chooses a cell as the next desired flight point without considering the target information. If the choice conflicts with other UAVs or the next flight point conflicts with the obstacle, choose another one.We use the same simulation environment in Scenario  1. The above three algorithms are implemented and compared from the following aspects: average time required to find all the 6 targets, success rate of finding all the targets in a given time period of 500 s, and repeated detection rate. All the data are averaged by50 runs. The results are shown in Table 2.Table 2
Results of different search algorithms.
Sweep-line
Random
Our cooperative
search algorithm
search algorithm
search algorithm
Average time required to search out all the targets (s)
461
N/A
56
Success rate in 500 s
64%
2%
94%
Repeated detection rate in 500 s
4.64
17.74%
6.45%From Table2, it can been seen that the average time for our proposed cooperative search algorithm to find all the targets is 56 s, which is more efficient and effective than both the sweep-line (461 s) and the random search algorithm (N/A, too many failures). For the second item, the success rate in 500 s of the cooperative search algorithm is 94
%, which is much larger than that of the sweep-line and random search algorithm. The success rate of the sweep-line search algorithm is 64
% due to the limited task time of 500 s, while the random search algorithm has an unsatisfying success rate of 2
%. For the third item, our proposed cooperative search algorithm performs the repeated detection rate as 6.45
% when escaping from the local optimum to find new targets. And the repeated search of the sweep-line algorithm only occurs when turning as 4.64
%. However, the random search algorithm may search for the same cell frequently with a repeated detection rate of 17.74
%.Further, 20 targets are randomly deployed in the nonconvex region. The number of targets found with time evolving is used as the measure of performance, as shown in Figure11. According to the simulation results, the proposed cooperative search algorithm can lead UAVs to find more targets in the same time duration than others, because the control law leads UAV to converge to the regions with high target existence probability. And if time is sufficient, the sweep-line algorithm will also find all the targets. However, the random search performs poorly due to its random mode.Figure 11
Comparison with sweep-line search and random search algorithm.
(a)
Sweep-line algorithm
(b)
Random search algorithm
(c)
Our cooperative algorithm
(d)
Number of targets found with time evolvingFrom all the results, it can be seen that our proposed cooperative search algorithm is quite efficient, for the search is led by the control law to optimize the information related objective function. And it can not only lead UAVs to converge to the regions with high target existence probability but also improve the exploration of the uncertain regions. Besides, the nonconvex environment is transformed into a convex one through the diffeomorphism at the beginning, so there is no excess step to deal with obstacle avoidance during the search process, dramatically speeding up the convergence.


## 5.1. Simulation Environment

The algorithms have been programmed in Matlab. The simulation scenario is set as that five UAVs search for targets in two typical nonconvex regions with differently shaped obstacles. The whole search regionQ is defined as a square region of [
0,10
]
×
[
0,10
] km2. And the distribution of each target follows the Gaussian distribution function(38)
ϕ
r
=
1
σ
2
π
e
-
r
-
μ
2
/
2
σ
2
,where the mean value μ is the expected position of the target and σ is set to be 1 km, representing the predefined probability distribution range. The communication range R
C is 5 km. Initially, all the UAVs are randomly deployed in the feasible region. Some other parameters for simulation are in Table 1.Table 1
Simulation parameters.
UAV
Sensor
Detection model
Maximum speed (m/s)
100
α
s
45°
Positive sensing probabilityP
0
i
0.9
Maximum turning angle
60°
β
s
30°
False alarm probabilityp
a
0.1
Flight height (m)
100
γ
s
30°
Target existence criteriaP
m
a
x
0.9


## 5.2. Simulation Results and Analysis

### 5.2.1. Scenario  1: Nonconvex Environment with Concave Obstacles

The first obstacle area is a typical concave region, represented as a red square. The convex decomposition and the diffeomorphism have been implemented as shown in Figures5(a), 6(a), and 6(b). Then the cooperative search algorithm is performed. The trajectories of UAVs at different time steps in the transformed region without obstacles and in the original region are shown in Figure 7, where the UAVs and targets are denoted by “o” and “+”, respectively. The probability map with time evolving is shown in Figure 8. From the figures, it can be seen that the probability map has converged. UAVs in the original region, starting from any initial positions, can find all the 6 targets at time 45 s while avoiding the obstacle.Figure 7
Search flight trajectories of UAVs with time: (a–c) in the transformed region and (d–f) in the original region.
(a)
0 s
(b)
30 s
(c)
60 s
(d)
0 s
(e)
30 s
(f)
60 sFigure 8
Probability map update for Scenario  1.
(a)
0 s
(b)
30 s
(c)
60 s


### 5.2.2. Scenario  2: Nonconvex Environment with Several Complex Obstacles

The nonconvex region in Scenario  2 contains two obstacles. The results of cooperative search are shown in Figures9 and 10. It can be seen that the trajectories of UAVs in the original region governed by the control law (31), starting from any initial positions, can move towards the targets while avoiding the obstacle. It also can be seen that the probability map has converged. After finding the targets, the UAV uses the uncertainty map instead of the probability map to escape from the local optimum (the positions of targets with high target existence probability) and keep searching. From the simulation results, all the targets are found at time 54 s.Figure 9
Search flight trajectories of UAVs with time: (a–c) in the transformed region and (d–f) in the original region.
(a)
0 s
(b)
30 s
(c)
60 s
(d)
0 s
(e)
30 s
(f)
60 sFigure 10
Probability map update for Scenario  2.
(a)
0 s
(b)
30 s
(c)
60 s


### 5.2.3. Comparisons with Other Algorithms

In order to verify the performance, the proposed algorithm is compared with a sweep-line search algorithm [7] and a random search algorithm [32]. In the sweep-line search algorithm, the search mood is fixed, and each UAV sweeps the region parallel to the boundary and decides the width between two search lines based on the detection range. In the random search algorithm, each UAV randomly chooses a cell as the next desired flight point without considering the target information. If the choice conflicts with other UAVs or the next flight point conflicts with the obstacle, choose another one.We use the same simulation environment in Scenario  1. The above three algorithms are implemented and compared from the following aspects: average time required to find all the 6 targets, success rate of finding all the targets in a given time period of 500 s, and repeated detection rate. All the data are averaged by50 runs. The results are shown in Table 2.Table 2
Results of different search algorithms.
Sweep-line
Random
Our cooperative
search algorithm
search algorithm
search algorithm
Average time required to search out all the targets (s)
461
N/A
56
Success rate in 500 s
64%
2%
94%
Repeated detection rate in 500 s
4.64
17.74%
6.45%From Table2, it can been seen that the average time for our proposed cooperative search algorithm to find all the targets is 56 s, which is more efficient and effective than both the sweep-line (461 s) and the random search algorithm (N/A, too many failures). For the second item, the success rate in 500 s of the cooperative search algorithm is 94
%, which is much larger than that of the sweep-line and random search algorithm. The success rate of the sweep-line search algorithm is 64
% due to the limited task time of 500 s, while the random search algorithm has an unsatisfying success rate of 2
%. For the third item, our proposed cooperative search algorithm performs the repeated detection rate as 6.45
% when escaping from the local optimum to find new targets. And the repeated search of the sweep-line algorithm only occurs when turning as 4.64
%. However, the random search algorithm may search for the same cell frequently with a repeated detection rate of 17.74
%.Further, 20 targets are randomly deployed in the nonconvex region. The number of targets found with time evolving is used as the measure of performance, as shown in Figure11. According to the simulation results, the proposed cooperative search algorithm can lead UAVs to find more targets in the same time duration than others, because the control law leads UAV to converge to the regions with high target existence probability. And if time is sufficient, the sweep-line algorithm will also find all the targets. However, the random search performs poorly due to its random mode.Figure 11
Comparison with sweep-line search and random search algorithm.
(a)
Sweep-line algorithm
(b)
Random search algorithm
(c)
Our cooperative algorithm
(d)
Number of targets found with time evolvingFrom all the results, it can be seen that our proposed cooperative search algorithm is quite efficient, for the search is led by the control law to optimize the information related objective function. And it can not only lead UAVs to converge to the regions with high target existence probability but also improve the exploration of the uncertain regions. Besides, the nonconvex environment is transformed into a convex one through the diffeomorphism at the beginning, so there is no excess step to deal with obstacle avoidance during the search process, dramatically speeding up the convergence.


## 5.2.1. Scenario  1: Nonconvex Environment with Concave Obstacles

The first obstacle area is a typical concave region, represented as a red square. The convex decomposition and the diffeomorphism have been implemented as shown in Figures5(a), 6(a), and 6(b). Then the cooperative search algorithm is performed. The trajectories of UAVs at different time steps in the transformed region without obstacles and in the original region are shown in Figure 7, where the UAVs and targets are denoted by “o” and “+”, respectively. The probability map with time evolving is shown in Figure 8. From the figures, it can be seen that the probability map has converged. UAVs in the original region, starting from any initial positions, can find all the 6 targets at time 45 s while avoiding the obstacle.Figure 7
Search flight trajectories of UAVs with time: (a–c) in the transformed region and (d–f) in the original region.
(a)
0 s
(b)
30 s
(c)
60 s
(d)
0 s
(e)
30 s
(f)
60 sFigure 8
Probability map update for Scenario  1.
(a)
0 s
(b)
30 s
(c)
60 s


## 5.2.2. Scenario  2: Nonconvex Environment with Several Complex Obstacles

The nonconvex region in Scenario  2 contains two obstacles. The results of cooperative search are shown in Figures9 and 10. It can be seen that the trajectories of UAVs in the original region governed by the control law (31), starting from any initial positions, can move towards the targets while avoiding the obstacle. It also can be seen that the probability map has converged. After finding the targets, the UAV uses the uncertainty map instead of the probability map to escape from the local optimum (the positions of targets with high target existence probability) and keep searching. From the simulation results, all the targets are found at time 54 s.Figure 9
Search flight trajectories of UAVs with time: (a–c) in the transformed region and (d–f) in the original region.
(a)
0 s
(b)
30 s
(c)
60 s
(d)
0 s
(e)
30 s
(f)
60 sFigure 10
Probability map update for Scenario  2.
(a)
0 s
(b)
30 s
(c)
60 s


## 5.2.3. Comparisons with Other Algorithms

In order to verify the performance, the proposed algorithm is compared with a sweep-line search algorithm [7] and a random search algorithm [32]. In the sweep-line search algorithm, the search mood is fixed, and each UAV sweeps the region parallel to the boundary and decides the width between two search lines based on the detection range. In the random search algorithm, each UAV randomly chooses a cell as the next desired flight point without considering the target information. If the choice conflicts with other UAVs or the next flight point conflicts with the obstacle, choose another one.We use the same simulation environment in Scenario  1. The above three algorithms are implemented and compared from the following aspects: average time required to find all the 6 targets, success rate of finding all the targets in a given time period of 500 s, and repeated detection rate. All the data are averaged by50 runs. The results are shown in Table 2.Table 2
Results of different search algorithms.
Sweep-line
Random
Our cooperative
search algorithm
search algorithm
search algorithm
Average time required to search out all the targets (s)
461
N/A
56
Success rate in 500 s
64%
2%
94%
Repeated detection rate in 500 s
4.64
17.74%
6.45%From Table2, it can been seen that the average time for our proposed cooperative search algorithm to find all the targets is 56 s, which is more efficient and effective than both the sweep-line (461 s) and the random search algorithm (N/A, too many failures). For the second item, the success rate in 500 s of the cooperative search algorithm is 94
%, which is much larger than that of the sweep-line and random search algorithm. The success rate of the sweep-line search algorithm is 64
% due to the limited task time of 500 s, while the random search algorithm has an unsatisfying success rate of 2
%. For the third item, our proposed cooperative search algorithm performs the repeated detection rate as 6.45
% when escaping from the local optimum to find new targets. And the repeated search of the sweep-line algorithm only occurs when turning as 4.64
%. However, the random search algorithm may search for the same cell frequently with a repeated detection rate of 17.74
%.Further, 20 targets are randomly deployed in the nonconvex region. The number of targets found with time evolving is used as the measure of performance, as shown in Figure11. According to the simulation results, the proposed cooperative search algorithm can lead UAVs to find more targets in the same time duration than others, because the control law leads UAV to converge to the regions with high target existence probability. And if time is sufficient, the sweep-line algorithm will also find all the targets. However, the random search performs poorly due to its random mode.Figure 11
Comparison with sweep-line search and random search algorithm.
(a)
Sweep-line algorithm
(b)
Random search algorithm
(c)
Our cooperative algorithm
(d)
Number of targets found with time evolvingFrom all the results, it can be seen that our proposed cooperative search algorithm is quite efficient, for the search is led by the control law to optimize the information related objective function. And it can not only lead UAVs to converge to the regions with high target existence probability but also improve the exploration of the uncertain regions. Besides, the nonconvex environment is transformed into a convex one through the diffeomorphism at the beginning, so there is no excess step to deal with obstacle avoidance during the search process, dramatically speeding up the convergence.


## 6. Conclusion and Future Work

This paper presents a distributed cooperative search method with multiple UAVs in a complex nonconvex environment with arbitrary obstacles. The proposed method mainly focuses on two issues: how to integrate the information update into the coverage control scheme and how to address the search problem in a complex nonconvex environment. First, an asynchronous distributed cooperative search framework is proposed, integrating information update into the coverage control scheme. The information update includes the Bayesian theory based individual update and a consensus fusion protocol based fusion update, while taking the inaccuracies and uncertainties in the sensor information into consideration. And an adaptive density function of coverage optimization is formulated based on the probability and uncertainty maps. Second, a new transformation method is proposed in order to extend the diffeomorphism to deal with the nonconvex environment with concave obstacles. The control strategy is presented considering the kinematic constraints, and the convergence is proved by LaSalle’s invariance principle. Finally, the effectiveness of the proposed algorithm is demonstrated through simulations. In future work, we will extend the algorithm to a high-dimension space and take the complex dynamics of the UAV into consideration.



---

*Source: 196730-2015-08-17.xml*