Cooperative Search by Multiple Unmanned Aerial Vehicles in a Nonconvex Environment

. 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.


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][2][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 Section 2, 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.

Preliminaries
2.1.Basic Notions and Definitions.Let R be the set of real numbers and let R  denote the -dimensional Euclidean space.For a given set  ⊆ R  ,  int and  denote the interior and the boundary of , respectively.One has  int = /.(, ) = ‖ − ‖ represents the Euclidean distance between points ,  ∈ R 2 , and (, ) is a closed ball centered at  with radius .In the following, some definitions about the environment are presented.
Definition 4 (nonconvex allowable environment [22]).Let  ⊂ R 2 be a finite set that is allowable if (i)  is compact and connected; (ii)  is continuously differentiable except on a finite number of points; (iii)  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 by , which is a polygon in R 2 .It is assumed that there are several polygonal obstacles   ⊂ ,  = 1, . . .,  within .Thus, the overall feasible region is  =  \ ⋃  =1   .Obviously, the search region is nonconvex since it contains several obstacles.

Detection and Communication Model.
Considering  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  is denoted by   = [  ,   ]  , and the heading angle is   .So the state of UAV  can be represented as [  ,   ,   ].
Each UAV is equipped with a camera to take measurements within its Field of View (FoV).Supposing UAV  is aiming at point  at time  (Figure 1), the formulation of FoV is where ℎ ≤ ℎ max represents the flight height of UAV and   ,   , and   are the installed angle and view angles of the camera, related to the physical performance.
The search region  is then discretized into several cells, defined as   .   , represents that a target exists in cell   at time  and    , represents that no target exists in cell   at time .Similarly,    , represents that the UAV detects a target in cell   at time  and    , represents that the UAV detects no target in cell   at time .Here, we only consider the influence of the Euclidean distance ‖  −   ‖ between UAV  and the detected cell   .The detection probability for UAV  can be represented by a monotonically decreasing function as where  0 ∈ (0, 1] is a constant of the sensing performance for UAV  and  is a positive constant, set to be 1.All the UAVs are provided with the identical kind of camera sensors, so  0 =  0 =  0 (,  = 1, . . ., ,  ̸ = ).In this paper, the communication range is assumed to be limited.Let CN = (, ]) be an undirected graph, where ] = {1, 2, . . ., } and  = {(, ) : ,  ∈ ]} represent the vertex set and the edge set, respectively.The communication network for UAV  at time  is modeled as CN , = (  , ]).The edge set is   = {(, ) : ,  ∈ ] | (  ,   ) ≤   }, where   is the communication range.Therefore, the set of UAV 's neighbors can be defined as  , = { ∈ ] | (, ) ∈   }. "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.

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 Figure 2. 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 TP  and uncertainty map   .At first, each UAV takes a measurement   and updates its local probability map TP  through Bayesian update rule, where the inaccuracies and uncertainties in the measurements are taken into consideration.Then, UAVs in the same communication network CN  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   .The uncertainty map   will be obtained by using both position and probability map.Finally, each UAV generates its control law by optimizing the objective function   .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.

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: where  ∈  is the cell of environment,   represents the Voronoi region generated by   , and (‖  − ‖) is a nondecreasing function.The density function () :  → R + represents the (relative) importance level of each cell in the environment .The more important the cell  is, the larger the value () would be.Generally speaking, more attention should be paid to the region with a larger density function.In order to optimize the objective function (), the partial derivative of the objective function should be calculated.Supposing (‖  − ‖) = ‖  − ‖ 2 in this paper, it can be obtained that where    and    are the mass and the centroid with respect to the density function (), defined as For simplification, the dynamics of each UAV is chosen as ṗ  =   .The gradient descent approach is used to minimize the objective function, expressed as where   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.

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 () 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 map   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: ) = 0 and  has been detected 0.5, if TP , () = 0 and  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: where Stage  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 Stage  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.

Information Update Method.
In the cooperative search problem, each UAV keeps an individual probability map TP  .The information update is in the form of probability map update, and the uncertainty map   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.

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 probability (   , ) is known in advance. FOV  represents that a target is detected by UAV  at time  with respect to the sensor measurement  , = 1. FOV  is the opposite of  FOV  with respect to measurement  , = 0. Considering the different sensor measurements, the update processes are different.
(1) If a target is detected by UAV , the update of target probability at time  based on Bayesian theory can be expressed as where ( FOV  |    , ) is the probability of detecting the target by UAV , when a target exists in cell   , and (   ,0 ) is the initialized probability map.Since the FoV of UAV  contains   cells, ( FOV  |    , ) can be described by the probability of failing to detect the target in these   cells as where (   , |    , ) represents the probability that UAV  fails to detect the target in cell   when the target is in cell   .Assuming   is false alarm probability, (   , |    , ) is described as follows: where If the target is in the FoV of UAV, the detection fault of missing the target occurs in one cell with the probability  0 .If not, the detection result is correct, and only false alarm probability is eliminated.Therefore, the probability for UAV  of failing to detect the target is Combining ( 10) with ( 12), it can be obtained that Based on the Bayesian complete probability formula, ( FOV  ) can be expressed as Therefore, when the target is detected by UAV  in its FoV, the update expression of the target existence probability can be obtained as follows: Mathematical Problems in Engineering 7 (2) If the target is not detected by UAV  in its FoV, the update expression of the target existence probability can be obtained as follows: With the same evolving process as mentioned above, it can be obtained that where (   , ).As described in ( 15) and ( 17), the update of target probability with different sensor measurements at time  for UAV  can be summarized as where  ,, represents the weight of the target probability map for each UAV, which can be expressed as where | , | is the number of UAV 's neighbors at time .

Nonconvex Environment with Arbitrary Obstacles
Considering the search environment  ⊆ R 2 is nonconvex with arbitrary obstacles , 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 manifolds  and , a differentiable map  :  →  is called a diffeomorphism, if it is a bijection and its inverse  −1 :  →  is differentiable as well.Lemma 6.For a given set  ⊆ R 2 , let R ⊂ R 2 be a compact, convex, and simply connected partition of  and let  ⊆ R be a convex open set.Let  ∩  R = 0 and  = R/.Given   ∈  int , there is a.e. 1 diffeomorphism  between  int and Rint \ {  }.
The above lemma follows that  is convex or starshaped.However, the real obstacle  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.

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  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  −1 ,   , and  +1 from the point sequence of a polygon in counterclockwise direction.Obtain the vectors   →  −1   and   →    +1 , and then compute the cross product where   and   are the coordinates of point   .If  is negative, then   is a concave point; otherwise,   is a convex point.
(  3. (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 (, ) = | − |, 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: 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 partition R = R1 , . . ., R of .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.

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 Lemma 6 requires that  ∩  R = 0.However, the new decomposed obstacles may have the common edges with the subregions such that   ∩  R ̸ = 0.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.Let R = R1 , . . ., R be a compact, convex, and simply connected partition of , representing the subregions.For each R , let   ⊂ R be a convex open set, representing the decomposed obstacle.Thus, the feasible region is   = R \  .For this problem, the edges generated by the decomposing lines must be considered.Let  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   ∩  R = 0, the diffeomorphism can be constructed using the method in the proof of Lemma 1 in [16].

Situation 4 (three common edges).
There are two cases.
(i) For the first case in Figure 4(d), we first define the critical vertex as the common vertex of the subregion and Mathematical Problems in Engineering obstacle, while one of its adjacent vertices does not belong to the obstacle, such as    and  2 .Choose one of the critical vertices as    , then   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    with   .The first diffeomorphism   () can be constructed by drawing a radial line from    using (23) as in Situation 2. For the other part, the diffeomorphism should be modified a little by drawing a radial line from   , as follows: where  1 and  0 are the two intersection points of the radial line with the boundaries of the obstacle, and if  lies on the edge  R in ,  0 and  1 are the corresponding vertices of the obstacle.Note that, through the diffeomorphism   (), all the points (except    and   ) on the dividing line have the same mapping with the first diffeomorphism   ().And   has a mapping by the diffeomorphism   .So we obtain   =   ∪   that can map  in  to R in \ {   }. (ii) For the second case in Figure 4(e), two edges of   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   () = .For the other partition, a parallel line is drawn from  and intersects with the boundaries of R at the left point  0 and right point  3 .Given the vertices of the polygon,  0 lying on the line  1  1 can be mapped to  1 on the line     1 , described as  1 =    − (‖   −  1 ‖‖ 1 −  0 ‖)/(‖ 1 −  1 ‖ ).Draw another parallel line from  1 , and it intersects with the boundary of R at point  2 .Thus the diffeomorphism can be constructed as   () =  1 + (‖ −  0 ‖‖ 2 −  1 ‖)/(‖ 3 −  0 ‖).Combining with the representation of  1 , the diffeomorphism can be described as For this case,   =   ∪   .
Situation 5 (more than three common edges in Figure 4(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   .As claimed in the above six situations, it is obvious that the diffeomorphism   can be constructed for each subregion.
Remark 8.The subregion containing no obstacle will be transformed by an identical transformation, which is also a diffeomorphism.
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.

Numerical Examples.
Two examples of applying the convex decomposition method to divide the nonconvex region with arbitrary obstacles are shown in Figure 5.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.
Further, the diffeomorphism is constructed for each example region by using the method in Proposition 7. We sample some points in the original regions and then map them to the transformed regions by the diffeomorphism, as shown in Figure 6.

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.

Objective Function Transformation.
The objective function (3) should be modified by shifting the allowable search domain from the feasible region  int to the whole region  int \ {  }.Based on the modification of Variables Theorem [29] where | det   2 (ω)| is the determinant of the Jacobian of  2 .In this paper, we assume that It can be seen that the allowable search domain is transformed from the feasible region  int to the whole region  int \ {  }, 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 Since  is the diffeomorphism,  int and ( int ) are topologically equivalent, and the problem remains essentially unchanged.Further, the corresponding function in the transformed domain is obtained: So the gradient descent approach as (6) can be used to minimize the objective function where  2  is a positive gain,   = (  ) is the transformed position of   by the diffeomorphism in the transformed region  int /(  ), and C  is the centroid of the Voronoi region generated by   .

Control Strategy for Cooperative Search in the Nonconvex
Environment.Since the diffeomorphism  is a bijection, the trajectory of   in  int /{  } has the unique preimage in original region  int , which is the trajectory of   .So the control law for the new objective function (28) in the original region can be obtained by the inverse mapping of (30): where  2  is a positive gain and   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].
Proof.Since the diffeomorphism  is a bijection, the trajectory   in  int \ {  } has the unique preimage   in original region  int .Therefore, to prove that the trajectory   governed by the control law (31), starting from any initial condition (0), can converge to  −1 ( C  ), we only need to prove that the trajectory   governed by the control law (30), starting from initial condition (0) = ((0)) = (( 1,0 ), . . ., ( ,0 )), can converge to C  .
where  1 = | det   −1 ()| and M  is the mass of the Voronoi region generated by   in the convex region  int \ {  }.
(v)  itself is the largest invariant subset of  by the control law (30).
Thus, by LaSalle's invariance principle, the trajectories   governed by (30), starting from initial configuration (0), will asymptotically converge to the critical points of (, Ṽ); that is, the trajectories of UAVs governed by the control law (31), starting from any initial condition (0) = ( 1,0 , . . .,  ,0 ) ∈  int , will asymptotically converge to the desired location.Therefore, the discrete-time motion model and the real control law for UAV  in the original nonconvex region can be written as 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 velocity V max will limit the control law by the following saturation rule: (2) Maximum turning angle: the maximum turning angle during a unit time period Δ is derived from the turning rate  max , described as  max =  max Δ.Thus the control law is limited by the following equations: where  is the rotating matrix for the maximum turning angle from the last control input as represents the rotating direction.Combining (34) with (35), the control law can be expressed as where  , = − 2   −1  (( ,−1 ))( ,−1 −  −1 ( C  )), as shown in (33).

Cooperative Search Algorithm in the Nonconvex Environment with Arbitrary
Obstacles.The cooperative search algorithm for the nonconvex region with arbitrary obstacles is presented as Algorithm 1.
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  int is transformed into the new convex one,  int \ {  }.And the positions of UAVs are also transformed into new virtual ones.
(4) The Voronoi partition of transformed region  int \ {  } is generated by the virtual positions of UAVs.
Compute the centroid C  of each Voronoi region using probability map as the density.
(5) At each time step, UAV  moves towards the preimages of its generalized Voronoi centroid  −1 ( C  ) in the original region, governed by the control law (31).

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.

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 region  is defined as a square region of [0, 10] × [0, 10] km 2 .And the distribution of each target follows the Gaussian distribution function 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   is 5 km.Initially, all the UAVs are randomly deployed in the feasible region.Some other parameters for simulation are in Table 1.The probability map with time evolving is shown in Figure 8.

Simulation Results and Analysis
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.(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.

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

Figure 4 :
Figure 4: Diffeomorphism for the new decomposed regions.

Figure 6 :𝜔‖ 2
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.

Figure 7 :
Figure 7: Search flight trajectories of UAVs with time: (a-c) in the transformed region and (d-f) in the original region.

( 6 )
Update the probability map TP  of UAV .If the target probability of a cell is greater than  max , a target is found.(7) Update the Voronoi partitions based on the new virtual positions of UAVs.Compute the generalized Voronoi centroids C  .(8) If UAV  reaches the preimage of its Voronoi centroid  −1 ( C  ) and does not move, then use its uncertainty map instead of probability map to compute the new Voronoi centroid  −1 ( C  ) from this time step.(9) If the terminal condition is satisfied, the algorithm ends; otherwise, go to (5).

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 Figures 5(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 Figure7, where the UAVs and targets are denoted by "o" and "+", respectively.
Number of targets found with time evolving

Figure 11 :
Figure 11: Comparison with sweep-line search and random search algorithm.
The visibility set is the set of points in  visible from .We connect   with other vertices of the polygon in counterclockwise direction.For a given vertex   , if all the rest vertices behind   are on the left side of line     , then   is visible from point   , as shown in Figure 2) Find the Visible Point Set.Point  ∈  is visible from  ∈  if [, ] ⊂ , where [, ] represents the closed path segment [, ] = { + ( − ) |  ∈ [0, 1]}.
int  =   \ (  ∩ ) be the boundary   but not belong to .Similarly, let  R Given    ∈   , there is a.e. 1 diffeomorphism   between  in  and R in \ {   }, where  in  and R int =  R \ ( R ∩ ) be the boundary  R but not belong to .Proposition 7.
Obviously, it can be obtained that  1 ̸ =  0 .Then a simple diffeomorphism   can be constructed as So any  ∈  in  on the radial line  can be mapped to   () ∈ R in \ {   }.Due to our hypothesis on R and   , it can be seen that   is continuous and differentiable everywhere in  in  , and it has a unique inverse  −1  , which is also continuous and differentiable everywhere in R in \ {   }.Therefore,   is a  1 diffeomorphism.Note that the intersection points  1 and  0 are the functions of .So even if  1 and  2 are equidistant to    , we still have   ( 1 ) ̸ =   ( 2 ).Hence, the inverse exists and   is a diffeomorphism.
(a)).Draw a radial line  from    to cross   and intersect with the boundaries of R and   at the points  1 and  0 , respectively.Note that if the radial line lies on  R in ,  0 and  1 are the vertices of   and R on the corresponding boundary, that is, points 1(2) and 1(2).
∘Target existence criteria  max 0.9 5.2.2.Scenario 2: Nonconvex Environment withSeveral Complex Obstacles.The nonconvex region in Scenario 2 contains two obstacles.The results of cooperative search are shown in Figures 9 and 10.It can be seen that the trajectories of UAVs in the original region governed by the control law