Collective Perception Using UAVs: Autonomous Aerial Reconnaissance in a Complex Urban Environment

This article examines autonomous reconnaissance in a complex urban environment using unmanned aerial vehicles (UAVs). Environments with many buildings and other types of obstacles and/or an uneven terrain are harder to be explored as occlusion of objects of interest may often occur. First, in this article, the problem of autonomous reconnaissance in a complex urban environment via a swarm of UAVs is formulated. Then, the algorithm based on the metaheuristic approach is proposed for a solution. This solution lies in deploying a number of waypoints in the area of interest to be explored, from which the monitoring is performed, and planning the routes for available UAVs among these waypoints so that the monitored area is as large as possible and the operation as short as possible. In the last part of this article, two types of main experiments based on computer simulations are designed to verify the proposed algorithms. The first type focuses on comparing the results achieved on the benchmark instances with the optimal solutions. The second one presents and discusses the results obtained from a number of scenarios, which are based on typical reconnaissance operations in real environments.


Introduction
Recently, swarm robotics has become a phenomenon in many real-world applications. It is about the coordination of multiple robots as a system to achieve a desired behavior. Collection of information about the environment is the most critical prerequisite for decision-making. For this, Unmanned Aerial Vehicles (UAVs) equipped with necessary sensors may be used in many situations and applications both in the civil and military domains.
In this article, the problem of using a swarm of UAVs in military reconnaissance operations is formulated and the solution to this problem is proposed. The principles of collective perception are used during the process in order to accelerate the operation and increase its effectiveness. The reconnaissance operations are assumed to be performed in complex urban environments and/or very uneven terrain where obstacles or terrain may cause occlusions.
The article is organized as follows: in this section, the authors' motivation is discussed along with the implementation of the results. Contributions of the research are also highlighted. Section 2 reviews

•
The problem of autonomous aerial reconnaissance in complex environments via a swarm of unmanned aerial vehicles is formulated. The monitoring is performed from a number of waypoints deployed in the area of operations; each waypoint is defined in the three-dimensional space by its coordinates and altitude.

•
The approach to effectively evaluate the coverage (monitored area) from a number of deployed waypoints is proposed. Terrain and obstacles, which may occlude the objects of interest are taken into consideration as well as the parameters of the sensors used.

•
The algorithm to estimate a minimum number of waypoints needed to explore a required portion of the area of interest is proposed.
Sensors 2020, 20,2926 3 of 27 • The metaheuristic algorithm based on the simulated annealing principles is proposed to deploy waypoints in the area of operations (in three dimensions). • A set of experiments is designed to assess the proposed algorithms. The results are compared to the optimal solutions. • A set of scenarios is designed based on the parameters and features of the real typical reconnaissance operations, and the behavior and results are discussed. The real geographic data is used in these experiments.

Literature Review
The topic concerning autonomous aerial reconnaissance using UAVs is currently discussed very frequently by the scientific as well as military personnel. From the scientific point of view, the reconnaissance problem, as is understood in this article, can be seen as the Art Gallery Problem, which is a well-studied NP-hard visibility problem in computational geometry [9]. This problem has been studied mostly in the two-dimensional space [10,11].
The three-dimensional case of the problem (referred to as 3D Art Gallery Problem) is studied less frequently. For example, Marzal [12] aimed at the determination of the number of guards required to cover the interior of a pseudo-polyhedron as well as the placement of these guards; this study models the art gallery by an orthogonal pseudo-polyhedron. Savkin and Huang [13] estimated the minimal number of drones necessary to monitor a given area of a very uneven terrain. The proposed problem may be viewed as a drone version of the 3D Art Gallery Problem. Thanou and Tzes [14] addressed the area coverage problem of a 3D-space by a group of UAVs, equipped with omnidirectional cameras; the terrain is assumed known to each UAV which has a maximum flight-altitude.
Some of the publications dealing with the problem of using UAVs for monitoring or surveillance take into consideration the fact that the area, where reconnaissance is carried out, is full of visual obstacles, particularly in an urban area. Saripalli et al. [15] presented the design and implementation of a real-time vision-based approach to detect and track features in a structured environment using an autonomous helicopter. Shaferman and Shima et al. [16] examined the problem of tracking the moving ground target in an urban environment via a set of cooperating UAVs and proposed a stochastic search method (based on a genetic algorithm) for finding in real time monotonically improving solutions. Semsch et al. [17,18] dealt with the problem of multi-UAV-based surveillance in complex urban environments with occlusions. The problem lies in controlling the flight of UAVs with on-board cameras so that the coverage and recency of the information about the designated area may be maximized. Swarming of multiple UAVs in order to increase the performance and effectiveness of the operation has been recently studied in many researches including applications such as reconnaissance or surveillance. Alfeo et al. [19] proposed a swarm coordination bio-inspired algorithm and applied it to search for a target object using a swarm of UAVs equipped with imperfect sensors. Similarly, Li et al. [20] introduced a distributed algorithm for searching moving targets via a fleet of cooperative UAVs. Silva et al. [21] examined the problem of real-time object identification and tracking through cooperative UAVs in a complex and adversarial environment involving motion, crowded scenes and varied camera angles and proposed a distributed deep learning algorithm for a solution.
Effective route planning and trajectory optimization of UAVs are critical in applications such as monitoring, reconnaissance or surveillance. Extensive surveys of planning methods are provided, for example, by Cabreira et al. [22] (methods based on coverage path planning), Zhao et al. [23] (computational-intelligence-based methods), Coutinho et al. [24] (identification of common attributes in aerial planning problems), or Geiger [25] (methods such as linear programming, dynamic programming, genetic algorithms and neural networks). In their research, Reardon and Fink [26] connected the issues of aerial reconnaissance (3D Art Gallery Problem) and path planning (Traveling Salesman Problem) and formulated the problem of searching and identification of objects using both ground and aerial autonomous robotic systems. For the solution of the problem formulated in this article, two metaheuristic approaches were adapted: the simulated annealing for the waypoints deployment and the ant colony optimization for the path planning of UAVs. The choice of these methods has been supported by the experience of authors obtained in their previous research. Metaheuristic methods are, in general, problem dependent. Further on in this paragraph, several other methods, which could possibly be used to this problem, are mentioned on examples of their recent applications which are related to some extent to the problem examined in this article. A very popular method called Particle Swarm Optimization (PSO) is often used for various types of problems. It is a population based stochastic optimization technique developed in 1995 by Eberhart and Kennedy [27]. Since that time, a thousand of application have emerged. Shao et al. [28] use this method for the problem of generating cooperative feasible paths for formation rendezvous of UAVs which was formulated as a multi-objective optimization problem with many coupled constraints; the proposed algorithm can meet the kinematic constraints of UAVs and the cooperation requirements. Another very popular search heuristic approach belongs to the family of Genetic Algorithms (GA). It is inspired by Darwin's theory of natural evolution; it reflects the process of natural selection. Cao et al. [29] adopted this algorithm for the problem of multi-base multi-UAV cooperative reconnaissance path planning; the problem is transformed into the shortest path combinatorial optimization using graph theory. However, the authors do not assume the occlusion caused by the terrain or obstacles. Tabu search formulated in 1986 by Glover [30] is another well-known local search method. Lee, Chen and Lai [31] hybridized this method with the 2-opt algorithm in the problem of the mission route planning of multiple unmanned robots in order to distribute tasks and coordinate the operation. A large family of bio-inspired algorithms has become very popular in recent years. To name a few: Artificial Bee Colony [32] (problem of UAVs used for wireless communication networks); Bat Algorithm [33] (problem of tracking a dynamic invading target by an UAV); or Grey Wolf Optimization [34] (multi-UAV multi-target urban tracking problem).

Problem Definition
In this section, the problem examined in this article is formulated. The problem is about planning routes of available UAVs for the reconnaissance mission performed in a complex urban environment. The goal is to plan the route so that the mission is carried out as fast as possible while the exploration of the area of interest is as complete as possible. The terrain and obstacles which may cause occlusions are taken into consideration. The list of all symbols and variables used in this and the following sections is included in the table at the end of the article.
Let U = {U 1 , U 2 , . . . , U M } be a finite set of UAVs where M ≥ 1 is their number available to participate in a reconnaissance operation. At the beginning of the operation, each UAV is deployed in its initial position (base) in the area of operations. The UAVs launch from their bases, and, when the operation is over, they return back. Each UAV is equipped with a sensor system capable to monitor some portion of the area directly below it. These sensors are characterized by two parameters as follows: • Angular field of view (α f ov ).
• Maximum distance from monitored objects (d max ).
Both parameters are given by the technical construction and principles of the sensors. The sensors are homogeneous, i.e., all UAVs are equipped with the same type of sensors. Each UAV is able to monitor objects of interest located within the field of view of its sensor and not farther from the UAV than allowed by the maximum distance constraint. This principle in a plane is shown in Figure 1 where the green area (called object detection range) represents the space, in which the object, if located in this range, is detected by an UAV. The aim of the reconnaissance operation is to monitor as much portion of the area of interest as possible using a fleet of UAVs . The area of interest is defined by a polygon with or without holes lying in the area of operations ( ⊆ ). During the operation, the height of flight of the UAVs is restricted. The limiting parameters are as follows: • Minimum height of flight above the ground level (ℎ ).

•
Maximum height of flight above the ground level (ℎ ).
This means that the height of flight above the ground level ℎ of all UAVs at any moment of the operation must be within the allowed limits (ℎ ≤ ℎ ≤ ℎ ). Both limits are set by the commander based on the tactical requirements for the operation and/or technical parameters of sensors. Let = , , … , be a finite set of waypoints deployed in the area of operations where ≥ 1 is their number. The monitoring of the area of interest is performed from these waypoints only. When any UAV is located at waypoint ∈ at any time of the operation, some portion of the area of interest is observed (i.e., objects of interest are detected if located in the observed area). Let = , , … be an infinite set of points lying on the ground in the area of operations. Each point ∈ is characterized by its elevation (height above the mean sea level) determined by the terrain. An UAV monitoring from waypoint ∈ is located at some altitude ℎ ; see Equation (1) where ℎ is the height above the ground level of the UAV at waypoint . ℎ = + ℎ for all = 1,2, … , Let ⊆ be a set of points lying in the area of interest, which is monitored from waypoint ; this means that an object of interest is detected by an UAV positioned at waypoint and altitude ℎ if this object is located at any point from set . All points in must satisfy the conditions as follows: 1. Points in must lie in the area of interest ( ⊆ ). 2. Points in must be within the object detection range of the sensor (see Figure 1). 3. There must be a visual line of sight (VLOS) between the sensor and all points in (see Figure  2). Let = , , … , be a finite set of not-transparent obstacles lying in the area of interest where ≥ 0 is their number. An obstacle (e.g., a building) or an uneven terrain may occlude some portion of the area in the object detection range; this happens in case that the VLOS condition in point 3 above is violated. The principle in a plane is shown in Figure 2: Figure 2a shows the situation without obstacles, whereas Figure 2b with obstacles. The green line represents points in set ⊆ (for further reference called as visible), red line remaining points in but not in (for further reference called as occluded). The aim of the reconnaissance operation is to monitor as much portion of the area of interest AoI as possible using a fleet of UAVs U. The area of interest is defined by a polygon with or without holes lying in the area of operations AoO (AoI ⊆ AoO). During the operation, the height of flight of the UAVs is restricted. The limiting parameters are as follows: • Minimum height of flight above the ground level (h min ).

•
Maximum height of flight above the ground level (h max ).
This means that the height of flight above the ground level h of all UAVs at any moment of the operation must be within the allowed limits (h min ≤ h ≤ h max ). Both limits are set by the commander based on the tactical requirements for the operation and/or technical parameters of sensors.
Let W = {W 1 , W 2 , . . . , W N } be a finite set of waypoints deployed in the area of operations where N ≥ 1 is their number. The monitoring of the area of interest is performed from these waypoints only. When any UAV is located at waypoint W i ∈ W at any time of the operation, some portion of the area of interest is observed (i.e., objects of interest are detected if located in the observed area).
Let P = {P 1 , P 2 , . . .} be an infinite set of points lying on the ground in the area of operations. Each point P i ∈ P is characterized by its elevation E(P i ) (height above the mean sea level) determined by the terrain. An UAV monitoring from waypoint W i ∈ W is located at some altitude h i alt ; see Equation (1) where h i is the height above the ground level of the UAV at waypoint W i : Let V i ⊆ P be a set of points lying in the area of interest, which is monitored from waypoint W i ; this means that an object of interest is detected by an UAV positioned at waypoint W i and altitude h i alt if this object is located at any point from set V i . All points in V i must satisfy the conditions as follows:

1.
Points in V i must lie in the area of interest (V i ⊆ AoI).

2.
Points in V i must be within the object detection range of the sensor (see Figure 1).

3.
There must be a visual line of sight (VLOS) between the sensor and all points in V i (see Figure 2).
. . , O L } be a finite set of not-transparent obstacles lying in the area of interest where L ≥ 0 is their number. An obstacle (e.g., a building) or an uneven terrain may occlude some portion of the area in the object detection range; this happens in case that the VLOS condition in point 3 above is violated. The principle in a plane is shown in Figure 2: Figure 2a shows the situation without obstacles, whereas Figure 2b with obstacles. The green line represents points in set V i ⊆ P (for further reference called as visible), red line remaining points in P but not in V i (for further reference called as occluded).  Figure 3a shows an example of the real situation from the top view. The polygon (blue line) encloses the area of interest; grey objects are obstacles (buildings in this case). Green area represents visible points in set . Figure 3b shows the same situation, but this time, five waypoints are deployed. The total visible area is given by the unification of sets for individual waypoints according to Equation (2). The total coverage of the area of interest through deployed waypoints is calculated according to Equation (3) as a ratio between the number of visible points (denoted as | |) and the number of all points in the area of interest (denoted as | |).
During the reconnaissance operation, UAVs in the fleet visit individual waypoints in the specified order. As all sensors of the UAVs are homogeneous, it is irrelevant, which UAV will visit which waypoint. Let = , , … , be a set of routes of individual UAVs. Route ∈ defines a trajectory of UAV ∈ as an order of waypoints to be visited: where ≥ 0 is the number of waypoints to be visited by UAV . Constraints in Equations (4)- (7) need to be satisfied: constraint in Equation (4) ensures that each UAV launches from its base and returns back at the end of the operation; and constraints in Equations (5)-(7) ensure that all the waypoints are visited just once.
(a) (b)    Figure 3a shows an example of the real situation from the top view. The polygon (blue line) encloses the area of interest; grey objects are obstacles (buildings in this case). Green area represents visible points in set V i . Figure 3b shows the same situation, but this time, five waypoints are deployed. The total visible area V is given by the unification of sets V i for individual waypoints according to Equation (2). The total coverage of the area of interest through deployed waypoints is calculated according to Equation (3) as a ratio between the number of visible points (denoted as |V|) and the number of all points in the area of interest (denoted as |AoI|): During the reconnaissance operation, UAVs in the fleet visit individual waypoints in the specified order. As all sensors of the UAVs are homogeneous, it is irrelevant, which UAV will visit which waypoint. Let R = {R 1 , R 2 , . . . , R M } be a set of routes of individual UAVs. Route R j ∈ R defines a trajectory of UAV U j ∈ U as an order of waypoints to be visited: where K j ≥ 0 is the number of waypoints to be visited by UAV U j . Constraints in Equations (4)-(7) need to be satisfied: constraint in Equation (4) ensures that each UAV launches from its base and returns back at the end of the operation; and constraints in Equations (5)-(7) ensure that all the waypoints are visited just once:  Figure 3a shows an example of the real situation from the top view. The polygon (blue line) encloses the area of interest; grey objects are obstacles (buildings in this case). Green area represents visible points in set . Figure 3b shows the same situation, but this time, five waypoints are deployed. The total visible area is given by the unification of sets for individual waypoints according to Equation (2). The total coverage of the area of interest through deployed waypoints is calculated according to Equation (3) as a ratio between the number of visible points (denoted as | |) and the number of all points in the area of interest (denoted as | |).
During the reconnaissance operation, UAVs in the fleet visit individual waypoints in the specified order. As all sensors of the UAVs are homogeneous, it is irrelevant, which UAV will visit which waypoint. Let = , , … , be a set of routes of individual UAVs. Route ∈ defines a trajectory of UAV ∈ as an order of waypoints to be visited: where ≥ 0 is the number of waypoints to be visited by UAV . Constraints in Equations (4)- (7) need to be satisfied: constraint in Equation (4) ensures that each UAV launches from its base and returns back at the end of the operation; and constraints in Equations (5)-(7) ensure that all the waypoints are visited just once.
(a) (b) R k j ∈ W for all j = 1, 2, . . . , M and k = 1, 2, . . . , K j (5) R k j R q p for all j, p = 1, 2, . . . , M ( j p ), k = 1, 2, . . . , K j and q = 1, 2, . . . , K p Equation (8) defines the time, in which UAV U j performs its route R j . The term R k j − R k+1 j expresses the time needed to fly from waypoint/base R k j to the next waypoint/base R k+1 j on its route; it depends on the distance between both points and the flight parameters of UAV U j . The duration of the reconnaissance operation is defined in Equation (9); all UAVs start at the same time at the beginning of the operation, and the operation ends when the last UAV returns back to its base: A simple example with two UAVs (red dots) and five waypoints (blue dots) is shown in Figure 4. The routes (violet color) are ∈ for all = 1,2, … , and = 1,2, … , ≠ for all , = 1,2, … , ( ≠ ), = 1, 2, … , and = 1, 2, … , Equation (8) defines the time, in which UAV performs its route . The term − expresses the time needed to fly from waypoint/base to the next waypoint/base on its route; it depends on the distance between both points and the flight parameters of UAV . The duration of the reconnaissance operation is defined in Equation (9); all UAVs start at the same time at the beginning of the operation, and the operation ends when the last UAV returns back to its base.
A simple example with two UAVs (red dots) and five waypoints (blue dots) is shown in Figure  4. The routes (violet color) are = , , , and = , , , , . It is desired to plan the reconnaissance operation optimally. The optimization criteria are defined in Equations (10)- (12). The first optimization criterion in Equation (10) is connected with the second optimization criterion in Equation (11). The former minimizes the number of waypoints (the goal is to find such a number so that the monitoring may be performed in the requested quality from as low number of waypoints as possible), the latter maximizes the total coverage of the area of interest (the goal is to find positions of waypoints so that the visible area may be as large as possible). The third optimization criterion in Equation (12) is connected with the planning of routes for individual UAVs (the goal is to find such routes so that the operation time may be as short as possible).  It is desired to plan the reconnaissance operation optimally. The optimization criteria are defined in Equations (10)- (12). The first optimization criterion in Equation (10) is connected with the second optimization criterion in Equation (11). The former minimizes the number of waypoints N (the goal is to find such a number so that the monitoring may be performed in the requested quality from as low number of waypoints as possible), the latter maximizes the total coverage of the area of interest C (the goal is to find positions of N waypoints so that the visible area may be as large as possible). The third optimization criterion in Equation (12) is connected with the planning of routes R for individual UAVs (the goal is to find such routes so that the operation time T may be as short as possible): Sensors 2020, 20, 2926 The first two optimization criteria go against one another (the coverage is, generally, reduced when the number of waypoints is lowered). Thus, a new parameter called minimum coverage C min is defined; it controls the requested quality of the operation as it determines the minimum necessary portion of the visible area in the area of interest-see Equation (13). The objective is to find a low number of waypoints, with which the coverage is still equal or higher than requested: In general, each reconnaissance operation is characterized by a number of constant parameters and a number of optimization parameters (decision variables). Constant parameters are as follows: • Minimum and maximum permitted height of flight above the ground level (h min , h max ).
Optimization parameters are as follows (all variables are continuous except the first one, which is discrete): •

Number of waypoints (N). •
Positions of waypoints (coordinates x i and y i for all W i ∈ W).

•
Heights above the ground level of the UAVs at waypoints (h i for all W i ∈ W).

Solution Algorithms
In this section, algorithms for solving the problem formulated in the previous section are proposed. First, the principle for calculating the coverage of the area of interest is presented. Then, algorithms for the optimization of criteria in Equations (10)-(12) follow.

Evaluation of a Solution
Let X N be a particular solution, lying in the state space, to the reconnaissance operation. This solution is characterized by N waypoints W deployed in the area of operations: , there are 3N optimization variables. For a particular problem (settings of constant parameters) and solution X N in the state space, the value of the coverage of the area of interest C N can be calculated: C N = f X N ; the value lies in interval 0, 1 .
For calculation of C N , the number of visible points |V| in the area of interest needs to be determined-see Equation (3). To get this number, any point lying inside the area of interest must be evaluated if it is visible from any waypoint or not. In general, there are an infinite number of points P inside the area of interest, which is, of course, not possible to evaluate from the practical point of view. The rasterization of the area of interest is a possible solution to this. The principle is shown in Figure 5a. The area is evenly covered by a finite number of points P using the Sukharev grid. The rasterization step d rast determines the total number of points N P ; each point lies in the middle of its square. The visibility of each point in P is evaluated independently on one another; in this way, set V of visible points of size N V is created as illustrated in Figure 5b. The approximation of C N is then calculated according to Equation (14); the precision of the approximation can be controlled by the rasterization step d rast : Sensors 2020, 20, Evaluation of solution is performed by Algorithm 1. At the beginning, the set of visible points is emptied (line 1), and set is created based on the principles presented above (line 2). Then, each point ∈ is evaluated independently (line 3 to 7); if point is evaluated as visible from at least one waypoint ∈ (line 5), then this point is included in set and the algorithm continues by evaluation of the next point. The algorithm returns the coverage of solution (lines 8 and 9). .
if is visible from then 6.

7.
continue on line 3 8.  Evaluation of solution X N is performed by Algorithm 1. At the beginning, the set of visible points V is emptied (line 1), and set P is created based on the principles presented above (line 2). Then, each point P j ∈ P is evaluated independently (line 3 to 7); if point P j is evaluated as visible from at least one waypoint W i ∈ W (line 5), then this point is included in set V and the algorithm continues by evaluation of the next point. The algorithm returns the coverage C N of solution X N (lines 8 and 9).

Algorithm 1 Evaluation of a solution in pseudocode
Evaluation of solution is performed by Algorithm 1. At the beginning, the set of visible points is emptied (line 1), and set is created based on the principles presented above (line 2). Then, each point ∈ is evaluated independently (line 3 to 7); if point is evaluated as visible from at least one waypoint ∈ (line 5), then this point is included in set and the algorithm continues by evaluation of the next point. The algorithm returns the coverage of solution (lines 8 and 9). .
if is visible from then 6.

return
Algorithm 2 elaborates the key process on line 5 of Algorithm 1. It evaluates the visibility of point P j from waypoint W i by testing four conditions: (a) distance between W i and P j must not be larger than the maximum permitted distance d max (line 1); (b) P j must be in the field of view α f ov of the camera positioned in W i (line 2); (c) there is no obstacle O k ∈ O disrupting the VLOS between W i and P j (line 3 and 4); (d) there is no point P k ∈ P of elevation, which disrupts the VLOS between W i and P j (line 5 and 6).

Algorithm 2 Visibility evaluation between a point and a waypoint in pseudocode
Sensors 2020, 20, x FOR PEER REVIEW 10 of 27

Algorithm 2 Visibility evaluation between a point and a waypoint in pseudocode
The computational complexity of Algorithm 1 is in Equation (15). The first and the second terms are related to lines 3 and 4 of Algorithm 1: the visibility needs to be evaluated between rasterized points and (up to) waypoints. The third term represents the process on lines 3 and 4 of Algorithm 2: each obstacle from set of size may disrupt the VLOS. The last term is connected with lines 5 and 6 of Algorithm 2: elevation of points in may disrupt the VLOS (the reason of square root of is that only points in lying directly between and have to be tested).
The evaluation of a solution is a key process in the optimization of the problem. Therefore, the great efforts were devoted by the authors to increase its efficiency. Two main improvements were implemented: (a) the floating horizon algorithm for computing visibility of a set of points in a line with the linear dependence was implemented; (b) the obstacles were superimposed on the terrain. Thus, the resulting computational complexity of the improved algorithm was reduced as stated in Equation (16). Moreover, the implementation was optimized from the coding point of view (only operations with integers are employed in the critical parts of the algorithm). The process is also parallelized and distributed on the cores of a multicore processor.

Optimization of Waypoint Deployment
In this section, an algorithm for optimization of a particular number of waypoints is proposed. The optimization criterion is to maximize the coverage according to Equation (11). The input of the algorithm, beside the constant parameters and algorithm settings, is the number on waypoints to be deployed. The output is a particular solution = , , … , = , , ℎ , , , ℎ , … , , , ℎ and its coverage . The proposed metaheuristic algorithm is based on the simulated annealing principles, which are inspired in annealing in metallurgy (reducing defects of material by involving heating and controlled cooling). The algorithm works in iterations where the process of solution transformation is performed. The transformed solution may replace the original-even if it is worse. This principle allows expanding the search space and prevents the algorithm to be stuck in some local optimum.
The behavior or the algorithm is controlled by 5 parameters as follows: • if disrupts VLOS between and then return false 7. return true The computational complexity of Algorithm 1 is in Equation (15). The first and the second terms are related to lines 3 and 4 of Algorithm 1: the visibility needs to be evaluated between N P rasterized points and (up to) N waypoints. The third term represents the process on lines 3 and 4 of Algorithm 2: each obstacle from set O of size L may disrupt the VLOS. The last term is connected with lines 5 and 6 of Algorithm 2: elevation of points in P may disrupt the VLOS (the reason of square root of N P is that only points in P lying directly between W i and P j have to be tested): The evaluation of a solution is a key process in the optimization of the problem. Therefore, the great efforts were devoted by the authors to increase its efficiency. Two main improvements were implemented: (a) the floating horizon algorithm for computing visibility of a set of points in a line with the linear dependence was implemented; (b) the obstacles were superimposed on the terrain. Thus, the resulting computational complexity of the improved algorithm was reduced as stated in Equation (16). Moreover, the implementation was optimized from the coding point of view (only operations with integers are employed in the critical parts of the algorithm). The process is also parallelized and distributed on the cores of a multicore processor:

Optimization of Waypoint Deployment
In this section, an algorithm for optimization of a particular number of waypoints is proposed. The optimization criterion is to maximize the coverage C N according to Equation (11). The input of the algorithm, beside the constant parameters and algorithm settings, is the number on waypoints N to be deployed. The output is a particular solution The proposed metaheuristic algorithm is based on the simulated annealing principles, which are inspired in annealing in metallurgy (reducing defects of material by involving heating and controlled cooling). The algorithm works in iterations where the process of solution transformation is performed. The transformed solution may replace the original-even if it is worse. This principle allows expanding the search space and prevents the algorithm to be stuck in some local optimum.
The behavior or the algorithm is controlled by 5 parameters as follows: • Maximum temperature T max : the initial value of temperature used in the first iteration.

•
Minimum temperature T min : the threshold value of temperature (the algorithm ends when the temperature drops below this threshold).

•
Cooling coefficient λ: it controls the speed of temperature reduction by cooling in successive iterations.
• Maximum number of transformations in iteration n 1max : it controls the higher limit of transformations performed per iteration. • Maximum number of replacements in iteration n 2max : it controls the higher limit of replacements (accepting the transformed solution and replacing the original) per iteration.
Algorithm 3 shows the principles in pseudocode. At the beginning, a random solution is generated (line 1; see Algorithm 4 for details) and evaluated (line 2; see Algorithm 1 for details). The current temperature is set to the maximum limit T max (line 3). The algorithm works in iterations (lines 4 to 18); the temperature does not change within an iteration. In each iteration, a number of transformations are performed (lines 5 to 17). Transformation of the current solution X N into the new solution X N (line 7) is the key process of the algorithm (see Algorithm 5 for details).

Algorithm 3 Optimization of waypoint deployment in pseudocode
Sensors 2020, 20, x FOR PEER REVIEW 11 of 27 • Maximum number of replacements in iteration : it controls the higher limit of replacements (accepting the transformed solution and replacing the original) per iteration.
Algorithm 3 shows the principles in pseudocode. At the beginning, a random solution is generated (line 1; see Algorithm 4 for details) and evaluated (line 2; see Algorithm 1 for details). The current temperature is set to the maximum limit (line 3). The algorithm works in iterations (lines 4 to 18); the temperature does not change within an iteration. In each iteration, a number of transformations are performed (lines 5 to 17). Transformation of the current solution into the new solution (line 7) is the key process of the algorithm (see Algorithm 5 for details).

Algorithm 3 Optimization of waypoint deployment in pseudocode
The original solution is replaced by the transformed solution with some probability (line 9) according to the Metropolis criterion in Equation (17). In case that the transformed solution is better than the original, the original is always replaced. Otherwise, the probability depends on the difference of qualities of both solutions and the current temperature : the lower the difference and the higher the temperature, the bigger the probability. This means that in the first phases of the algorithm, the state space is largely expanded by accepting often worse solutions, whereas, towards the end of the optimization, the solution is tuned in its surrounding. The acceptance of the transformed solution is decided on line 10; function RandU 0,1 is a pseudorandom number generator with a uniform distribution.
An iteration ends when either the number of transformations or the number of replacements exceed their limits and respectively (see condition on line 6). Then the temperature is lowered (line 18) by cooling coefficient (0 < < 1) and the next iteration starts. The algorithm is terminated when the temperature drops below its lower limit (see condition on line 4).

return ,
The original solution is replaced by the transformed solution with some probability (line 9) according to the Metropolis criterion in Equation (17). In case that the transformed solution is better than the original, the original is always replaced. Otherwise, the probability depends on the difference of qualities of both solutions and the current temperature T cur : the lower the difference and the higher the temperature, the bigger the probability. This means that in the first phases of the algorithm, the state space is largely expanded by accepting often worse solutions, whereas, towards the end of the optimization, the solution is tuned in its surrounding. The acceptance of the transformed solution is decided on line 10; function RandU(0, 1) is a pseudorandom number generator with a uniform distribution: Sensors 2020, 20, 2926

of 27
An iteration ends when either the number of transformations n 1 or the number of replacements n 2 exceed their limits n 1max and n 2max respectively (see condition on line 6). Then the temperature is lowered (line 18) by cooling coefficient λ (0 < λ < 1) and the next iteration starts. The algorithm is terminated when the temperature T cur drops below its lower limit T min (see condition on line 4). The best solution found during the whole optimization X N best is stored (lines 14 to 16) and returned at the end of the algorithm (line 19).
Algorithm 4 shows the process of the generation of a random solution (named on line 1 of Algorithm 3). Each variable of a solution is set in its permitted limits using a pseudorandom number generator with a uniform distribution (functions MinX, MaxX, MinY and MaxY return the minimum and maximum values of a circumscribed rectangle of the area of interest). with a normal distribution with a mean of = 0 and a standard deviation of calculated on line 6 according to Equation (18). The standard deviation depends on the current temperature and the range of the variable (calculated on lines 3, 4 and 5 for individual types of variables: coordinate, coordinate and height). The bigger the temperature, the bigger the changes. Constant in Equation (18) is a small number influencing the standard deviation in situations when the current temperature is close to its minimum limit.

Algorithm 4 Generation of a random solution in pseudocode
The computational complexity of the algorithm for optimizing waypoint deployment is given in Equation (19). It is linearly dependent on the number of iterations computed according to Equation (20) and the maximum number of transformations ( is used instead of because the number of transformations is always equal to or greater than the number of replacements). The last two terms correspond to the evaluation of the transformed solution in the loop-see Equation (16).
The algorithm was implemented in two versions: continuous and discrete. The former assumes all the variables of the solution to be continuous-waypoints can be deployed wherever in the area of operations. The latter one assumes the variables to be discreet with the distance between neighboring values set to the size of the rasterization step . It has an effect that the waypoints can be deployed only in the middle of rasterization squares. The discretization of each variable is

return
Algorithm 5 presents the process (called on line 7 of Algorithm 3), in which the current solution X N is transformed into the new solution X N . In the transformation, one randomly selected variable is changed, the remaining variables copy the original values. The variable is selected by a pseudorandom integer number generator RandI(a, b) with a uniform distribution (line 2) in the range from 1 to 3N as the solution is composed of 3N variables: X N = x 1 , y 1 , h 1 , x 2 , y 2 , h 2 , . . . , x N , y N , h N = X N 1 , X N 2 , . . . , X N 3N . The size of the change is determined by a pseudorandom number generator RandN(µ, σ) with a normal distribution with a mean of µ = 0 and a standard deviation of σ calculated on line 6 according to Equation (18). The standard deviation depends on the current temperature and the range of the variable (calculated on lines 3, 4 and 5 for individual types of variables: x coordinate, y coordinate and height). The bigger the temperature, the bigger the changes. Constant ε in Equation (18) is a small number influencing the standard deviation in situations when the current temperature is close to its minimum limit: The computational complexity of the algorithm for optimizing waypoint deployment is given in Equation (19). It is linearly dependent on the number of iterations n iter computed according to Equation (20) and the maximum number of transformations n 1max (n 1max is used instead of n 2max because the number of transformations is always equal to or greater than the number of replacements). The last two terms correspond to the evaluation of the transformed solution in the loop-see Equation (16): Sensors 2020, 20, 2926 13 of 27 The algorithm was implemented in two versions: continuous and discrete. The former assumes all the variables of the solution to be continuous-waypoints can be deployed wherever in the area of operations. The latter one assumes the variables to be discreet with the distance between neighboring values set to the size of the rasterization step d rast . It has an effect that the waypoints can be deployed only in the middle of rasterization squares. The discretization of each variable is performed by rounding to the nearest permitted value.

Algorithm 5 Transformation of a solution in pseudocode
Sensors 2020, 20, x FOR PEER REVIEW 13 of 27 Algorithm 5 Transformation of a solution in pseudocode

Optimizing the Number of Waypoints
This section presents Algorithm 6 for optimizing the number of waypoints according to the optimization criterion in Equation (10). The main idea of the algorithm consists in the first estimation of the number based on the parameters of sensors and the subsequent deployment of this number of waypoints; in the next phases, the value is updated according to the gap between the actual and required coverage. Of course, this is not the only approach to determine the number of waypoints. Another approach could be the binary search method where the value is determined between the limits by the interval halving. The algorithm proposed in this section was selected because of its low number of phases (see the results in Section 5.3).

Algorithm 6 Optimizing the number of waypoints in pseudocode
The input of the algorithm is the minimum required coverage . At first, the initial value of using Equation (21) is estimated (line 1) based on the size of the area of interest and parameters of sensors. Coefficient ≥ 1 is a constant, which can be set with regard to the parameters of scenarios (e.g., the terrain and obstacles); in most real situations, the best results were achieved with = 1.5 (the purpose of this coefficient is discussed in Section 5.3 in more detail). Then, waypoints are optimized using Algorithm 3 (lines 2 and 3).

Optimizing the Number of Waypoints
This section presents Algorithm 6 for optimizing the number of waypoints according to the optimization criterion in Equation (10). The main idea of the algorithm consists in the first estimation of the number based on the parameters of sensors and the subsequent deployment of this number of waypoints; in the next phases, the value is updated according to the gap between the actual and required coverage. Of course, this is not the only approach to determine the number of waypoints. Another approach could be the binary search method where the value is determined between the limits by the interval halving. The algorithm proposed in this section was selected because of its low number of phases (see the results in Section 5.3).

Algorithm 6 Optimizing the number of waypoints in pseudocode
Sensors 2020, 20, x FOR PEER REVIEW 13 of 27

Optimizing the Number of Waypoints
This section presents Algorithm 6 for optimizing the number of waypoints according to the optimization criterion in Equation (10). The main idea of the algorithm consists in the first estimation of the number based on the parameters of sensors and the subsequent deployment of this number of waypoints; in the next phases, the value is updated according to the gap between the actual and required coverage. Of course, this is not the only approach to determine the number of waypoints. Another approach could be the binary search method where the value is determined between the limits by the interval halving. The algorithm proposed in this section was selected because of its low number of phases (see the results in Section 5.3).

Algorithm 6 Optimizing the number of waypoints in pseudocode
The input of the algorithm is the minimum required coverage . At first, the initial value of using Equation (21) is estimated (line 1) based on the size of the area of interest and parameters of sensors. Coefficient ≥ 1 is a constant, which can be set with regard to the parameters of scenarios (e.g., the terrain and obstacles); in most real situations, the best results were achieved with = 1.5 (the purpose of this coefficient is discussed in Section 5.3 in more detail). Then, waypoints are optimized using Algorithm 3 (lines 2 and 3).

Transform_Solution(
, , ) Input: , , Output: Constant parameters: ,ℎ ,ℎ Algorithm settings: , 1.  The input of the algorithm is the minimum required coverage C min . At first, the initial value of N using Equation (21) is estimated (line 1) based on the size of the area of interest and parameters of sensors. Coefficient τ ≥ 1 is a constant, which can be set with regard to the parameters of scenarios (e.g., the terrain and obstacles); in most real situations, the best results were achieved with τ = 1.5 (the purpose of this coefficient is discussed in Section 5.3 in more detail). Then, N waypoints are optimized using Algorithm 3 (lines 2 and 3): If the coverage does not satisfy the minimum limit (line 4), the value of N is updated (line 5) according to Equation (22) and the new number of waypoints are deployed again using Algorithm 3 (lines 6 and 7). This process is repeated in a loop (lines 4 to 7) until such N is found that the minimum coverage constraint is met. The algorithm returns the number N as well as the solution X N and its quality C N (line 8):

Planning of Routes
When the number of waypoints and their deployment in the area of operations are determined, the problem of planning routes for available UAVs follows. The optimization criterion is to minimize the duration T of the reconnaissance operation as stated in Equation (12). This problem can be easily transformed into the well-known Min-Max Multi-Depot Vehicle Routing Problem (MDVRP) [35] where a set of customers should be served by a fleet of vehicles originating from multiple depots.
The authors of this article studied this problem extensively in their previous research. The metaheuristic algorithm based on the Ant Colony Optimization (ACO) theory in combination with other principles was proposed as a solution. The results were published, for example, in [36,37]. Therefore, this topic will not be examined and further pursued in this article.

Experiments and Results
The solutions proposed in the previous section were verified on a series of experiments. All experiments were carried out on a computer with configuration as follows: CPU Intel i7-7700 3.5 GHz (4 cores), 32 GB RAM.

Evaluation of a Solution
The first set of experiments is aimed at the algorithm for the evaluation of a solution, i.e., calculation of the coverage for a particular solution, with regard to its efficiency and influence of the key parameters: number of waypoints and the total number of points in the area of interest.
The constant parameters of the scenario for these experiments are as follows: number of waypoints N = 5 and N = 10, respectively, field of view of sensors α f ov = 90, maximum permitted distance between a camera and objects d max   Table 1 presents the results achieved using the original algorithm with no code optimization (see Section 4.1) and the computational complexity given by Equation (15). In Table 2, the results are achieved when using the improved (optimized) algorithm with the computational complexity from Equation (16). Table 1. Experiments with the solution evaluation via the original algorithm.

Rasterization
Step d rast (m)

Rasterization
Step d rast (m) In general, the original algorithm can be assumed more precise in evaluation than the optimized version. In order to increase its efficiency, the latter one works with the integers rather than real numbers in its critical sections; this may lead to minor approximations and inaccuracies. To calculate the coverage error, the solution achieved by the original algorithm for d r = 1 is used as benchmark (both for N = 5 and N = 10).

Coverage
The difference in the execution time (time to evaluate a solution) is immense. The main reason is repeating the VLOS tests of each point with a number of L obstacles in the original algorithm (see Equation (15)), while this is not done in the optimized version (see Equation (16)). In average (based on all experiments in Tables 1 and 2), the optimized version is almost 60 times faster than the original algorithm. For example, in case of N = 5 and d r = 1, the optimized algorithm manages to evaluate more than 30 million points per second, whereas the original version evaluates only 0.5 million points. Figure 7 compares the coverage errors. Although the original version provides more precise results than the optimized version, the error of the latter one is below 2% in all cases (with the exception of d rast = 50, but the reduction of information is too big in this case). The illustration of the influence of the rasterization step on the evaluation is shown in Figure 8 for rasterization steps d rast = 1, d rast = 10 and d rast = 50.
Sensors 2020, 20, x FOR PEER REVIEW 16 of 27 Figure 7 compares the coverage errors. Although the original version provides more precise results than the optimized version, the error of the latter one is below 2% in all cases (with the exception of = 50, but the reduction of information is too big in this case). The illustration of the influence of the rasterization step on the evaluation is shown in Figure 8 for rasterization steps = 1, = 10 and = 50.   Figure 9 shows the dependence of the execution time on the number of points to be evaluated in a solution (for the original algorithm on the left, for the optimized version on the right). The linear dependence is apparent. The linear dependence on the number of waypoints N can also be seen in the graphs. However, it is not that twice the number of waypoints means twice the execution time (it is only about 1.6 times in this case). The reason is that when a point is once evaluated as visible, it does not have to be evaluated again from other waypoints; this case is more frequent when there are more waypoints.    Figure 9 shows the dependence of the execution time on the number of points to be evaluated in a solution (for the original algorithm on the left, for the optimized version on the right). The linear dependence is apparent. The linear dependence on the number of waypoints N can also be seen in the graphs. However, it is not that twice the number of waypoints means twice the execution time (it is only about 1.6 times in this case). The reason is that when a point is once evaluated as visible, it does not have to be evaluated again from other waypoints; this case is more frequent when there are more waypoints.  Figure 9 shows the dependence of the execution time on the number of points N P to be evaluated in a solution (for the original algorithm on the left, for the optimized version on the right). The linear dependence is apparent. The linear dependence on the number of waypoints N can also be seen in the graphs. However, it is not that twice the number of waypoints means twice the execution time (it is only about 1.6 times in this case). The reason is that when a point is once evaluated as visible, it does not have to be evaluated again from other waypoints; this case is more frequent when there are more waypoints.

Optimizing the Waypoint Deployment
A set of experiments was designed to validate the algorithm proposed in Section 4.2 for the deployment of a number of waypoints in the area of operations. Conditions and parameters for the benchmark scenarios are as follows: • The area of interest is assembled by joining a number of hexagons with the circumradius 100 m.

•
The terrain is absolutely flat (the altitude does not change within the area of operations).

•
There are no obstacles in the area of operations.

•
The number of waypoints to be deployed are the same as the number of hexagons.
The principle of creating the area of interest is presented in Figure 10. The degree determines the number of hexagons in the diagonal. In the example, the degree equals 5, which means that 17 hexagons were joined to create the area of interest and the same number of waypoints = 17 are to be deployed. This principle along with the conditions mentioned above ensures that the optimal solution can be easily found. When the waypoints are deployed optimally, the whole area of interest will be covered ( % = 100%). The optimal solution is as follows:

Optimizing the Waypoint Deployment
A set of experiments was designed to validate the algorithm proposed in Section 4.2 for the deployment of a number of waypoints in the area of operations. Conditions and parameters for the benchmark scenarios are as follows: • The area of interest is assembled by joining a number of hexagons with the circumradius 100 m.

•
The terrain is absolutely flat (the altitude does not change within the area of operations).

•
There are no obstacles in the area of operations. •

Parameters of sensors are
The number of waypoints to be deployed are the same as the number of hexagons.
The principle of creating the area of interest is presented in Figure 10. The degree determines the number of hexagons in the diagonal. In the example, the degree equals 5, which means that 17 hexagons were joined to create the area of interest and the same number of waypoints N = 17 are to be deployed.

Optimizing the Waypoint Deployment
A set of experiments was designed to validate the algorithm proposed in Section 4.2 for the deployment of a number of waypoints in the area of operations. Conditions and parameters for the benchmark scenarios are as follows: • The area of interest is assembled by joining a number of hexagons with the circumradius 100 m.

•
The terrain is absolutely flat (the altitude does not change within the area of operations).

•
There are no obstacles in the area of operations.

•
The number of waypoints to be deployed are the same as the number of hexagons.
The principle of creating the area of interest is presented in Figure 10. The degree determines the number of hexagons in the diagonal. In the example, the degree equals 5, which means that 17 hexagons were joined to create the area of interest and the same number of waypoints = 17 are to be deployed. This principle along with the conditions mentioned above ensures that the optimal solution can be easily found. When the waypoints are deployed optimally, the whole area of interest will be covered ( % = 100%). The optimal solution is as follows: This principle along with the conditions mentioned above ensures that the optimal solution can be easily found. When the waypoints are deployed optimally, the whole area of interest will be covered (C N % = 100%). The optimal solution is as follows: • The waypoints (coordinates x i , y i ) lie in the centers of the hexagons (see Figure 10).

•
The monitoring height is exactly 100 m above the ground level (h i = 100 m).
Six benchmark instances were created for verification. They differ by the degree as shown in Table 3. The fifth column of Table 3 shows the number of waypoints and the last column the number of variables per a solution vector X N . For example, in the case of instance d06, the solution is composed of more than 200 independent variables. The permitted range of the flight height does not differ in the instances. Both the continuous and discrete versions of the algorithm proposed in Section 4.2 (Algorithm 3) were used to find the solutions of the benchmark problems and the results were compared with the optimal solutions. The parameters of the algorithm were set as follows: T max = 10 −2 , T min = 10 −6 , γ = 0.9, n 1max = 200, n 2max = 20. The algorithm was executed 50 times on each benchmark instance.
The results achieved with the continuous version of the algorithm are recorded in Table 4. The optimal solution was found in the case of the simplest instances (d01 and d02). In other instances, a solution was found very close to the optimum. The difference between the best solution and the optimum (referred to as error in Table 4) is below 2% in all cases. The execution time depends on the rasterization step d rast .  Table 5 shows the results achieved by the discrete version of the algorithm. The optimum was found in case of instances d01, d02, d03; the maximum error in case of the most complex problems is below 1%. The execution time is slightly higher than in case of the continuous version; this is caused by the discretization of the variables in the transformation process. The comparison of the results achieved by the continuous and discrete versions of the algorithm is shown in Figure 11. In general, the discrete version provides better solutions than the continuous version; the more complex the problem, the better the results obtained by the discrete version. Therefore, the discrete version is used further on. Figure 12 shows the best solution achieved by the discrete version of the algorithm for instance d06 (C 71 % = 99.30%). The comparison of the results achieved by the continuous and discrete versions of the algorithm is shown in Figure 11. In general, the discrete version provides better solutions than the continuous version; the more complex the problem, the better the results obtained by the discrete version. Therefore, the discrete version is used further on. Figure 12 shows the best solution achieved by the discrete version of the algorithm for instance d06 ( % = 99.30%).

Optimizing the Number of Waypoints
In this section, the algorithm proposed for optimizing the number of waypoints (Algorithm 6) is verified on the benchmark problems. The same problems as proposed in the previous section are used.
The algorithm starts with the first estimation of according to Equation (21). The part of this equation is coefficient . This coefficient takes specific parameters and characteristics of the environment into consideration (e.g., the terrain and the number of obstacles). In the ideal scenario,  The comparison of the results achieved by the continuous and discrete versions of the algorithm is shown in Figure 11. In general, the discrete version provides better solutions than the continuous version; the more complex the problem, the better the results obtained by the discrete version. Therefore, the discrete version is used further on. Figure 12 shows the best solution achieved by the discrete version of the algorithm for instance d06 ( % = 99.30%).

Optimizing the Number of Waypoints
In this section, the algorithm proposed for optimizing the number of waypoints (Algorithm 6) is verified on the benchmark problems. The same problems as proposed in the previous section are used.
The algorithm starts with the first estimation of according to Equation (21). The part of this equation is coefficient . This coefficient takes specific parameters and characteristics of the environment into consideration (e.g., the terrain and the number of obstacles Error (%)

Benchmark instance
Continuous Discreet Figure 12. The best solution found for instance d06.

Optimizing the Number of Waypoints
In this section, the algorithm proposed for optimizing the number of waypoints (Algorithm 6) is verified on the benchmark problems. The same problems as proposed in the previous section are used.
The algorithm starts with the first estimation of N according to Equation (21). The part of this equation is coefficient τ. This coefficient takes specific parameters and characteristics of the environment into consideration (e.g., the terrain and the number of obstacles). In the ideal scenario, where the waypoints could be deployed in such a way that the whole area would be covered from a number of waypoints and, at the same time, areas visible from sensors would not overlap, the coefficient should be τ = 1 and the estimation would correspond to the correct number. In the benchmark problems proposed for verification, where the terrain is flat with no obstacles, the waypoints could be deployed so that the visible areas might overlap just slightly; therefore, the coefficient was set to τ = 1.1. In the real scenarios, the best value of the coefficient was empirically found as τ = 1.5. Table 6 shows the results of the algorithm. The minimal coverage is set to C min = 99%. The reason for this is that the full coverage represents the optimal solution, which can be hardly achieved in case of the more complex problems. Therefore, the algorithm error of 1% is assumed. As can be seen, the optimal numbers of waypoints are estimated in case of all the benchmark instances. Moreover, it was done in at most 3 phases: the first estimation (line 1 of Algorithm 6) and 2 updates (line 5 of Algorithm 6). Figure 13 illustrates the situation in case of instance d04.
Sensors 2020, 20, x FOR PEER REVIEW 20 of 27 benchmark problems proposed for verification, where the terrain is flat with no obstacles, the waypoints could be deployed so that the visible areas might overlap just slightly; therefore, the coefficient was set to = 1.1. In the real scenarios, the best value of the coefficient was empirically found as = 1.5. Table 6 shows the results of the algorithm. The minimal coverage is set to = 99%. The reason for this is that the full coverage represents the optimal solution, which can be hardly achieved in case of the more complex problems. Therefore, the algorithm error of 1% is assumed. As can be seen, the optimal numbers of waypoints are estimated in case of all the benchmark instances. Moreover, it was done in at most 3 phases: the first estimation (line 1 of Algorithm 6) and 2 updates (line 5 of Algorithm 6). Figure 13 illustrates the situation in case of instance d04.

Experiments on Real Scenarios
The experiments in this section are based on the real scenarios that reflect typical reconnaissance operations. The environment of the scenarios is based on the real geographic data using two models: (a) The Digital Elevation Model (DEM), and (b) The Topographic Digital Data Model (TDDM). The former is a representation of the terrain surface in the form of a heightmap. The latter is a database of topographic and other objects; in the scenarios, buildings are used as not-transparent obstacles.
Both geographic models come from the Military Geographic and Hydrometeorologic Office of the Ministry of Defense of the Czech Republic which provides geographic data for the Army of the Czech Republic. The DEM model has, in its last version, the distance between elevations 2.5 m and elevation precision 0.3 m; it is being updated regularly by methods of digital stereophotogrammetry and airborne laser scanning. The TDDM model contains topographic objects represented by a polygon and parameters (e.g., object height); it is being regularly updated using the method of direct mapping with the support of aerial imaging.  Table 6. Results of the algorithm for the optimizing the number of waypoints.

Instance
Optimal N First Estimation First Update Second Update

Experiments on Real Scenarios
The experiments in this section are based on the real scenarios that reflect typical reconnaissance operations. The environment of the scenarios is based on the real geographic data using two models: (a) The Digital Elevation Model (DEM), and (b) The Topographic Digital Data Model (TDDM). The former is a representation of the terrain surface in the form of a heightmap. The latter is a database of topographic and other objects; in the scenarios, buildings are used as not-transparent obstacles.
Both geographic models come from the Military Geographic and Hydrometeorologic Office of the Ministry of Defense of the Czech Republic which provides geographic data for the Army of the Czech Republic. The DEM model has, in its last version, the distance between elevations 2.5 m and elevation precision 0.3 m; it is being updated regularly by methods of digital stereophotogrammetry and airborne laser scanning. The TDDM model contains topographic objects represented by a polygon and parameters (e.g., object height); it is being regularly updated using the method of direct mapping with the support of aerial imaging. Table 7 characterizes the environment of the scenarios: the size of the area of interest, the elevation difference (difference between the highest and the lowest altitude inside the area of interest) and obstacles (their number and average height). The scenarios offer various types of environment: from a small to large area to be explored, from a relatively flat to a very uneven terrain, from a low density to a high density of obstacles. For example, scenario sc04 is a typical urban environment with an irregular shape of the area of interest, a very high density of tall obstacles (buildings) with narrow streets and a flat terrain. In comparison with this configuration, scenario sc05 is a large mountain environment with a very low density of obstacles, but a very uneven terrain.  Table 8 presents the technical and tactical configurations: the number of UAVs available (at the disposal of the commander and deployed in the area of operations), the minimum requested coverage, parameters of sensors (the maximum distance to objects of interest and an angular field of view) and a minimum and maximum permitted height of flight of the UAVs. For optimization, the parameters of the algorithms were set as follows: T max = 10 −2 , T min = 10 −6 , γ = 0.9, n 1max = 10·n 2max , n 2max = 20 (for sc01, sc02, sc05) or n 2max = 50 (sc03, sc04) respectively. Table 9 presents the optimization results. The third column shows the number of waypoints estimated by Algorithm 6. The estimated number of waypoints were deployed using Algorithm 3; in total, 50 optimizations per scenario were conducted-for results, see the fourth to sixth columns. The best solution found exceeds the minimum coverage required for the scenarios (see Table 8) in all cases. The last column shows the average execution time of optimization; this depends, of course, on the size of the rasterization step. Table 9. Optimizations of the number and deployment of waypoints for the real scenarios.
Step d rast (m)  Figure 14 illustrates the deployment of waypoints and coverage for scenarios sc03 (on the left) and sc04 (on the right). In case of the former one, the area of interest is of the very irregular shape. Despite of this, the algorithm managed to deploy all the waypoints inside the area (70 waypoints means 210 variables in the solution vector). In case of the latter one, the solution vector is composed of more than 300 independent variables. For the best solutions found, the routes of UAVs available in individual scenarios were planned using the algorithm mentioned in Section 4.4. The results are recorded in Table 10. The last two columns show the parameters of the optimized routes: the total distance covered by all UAVs participating in the operation, and operation time (the time when all the UAVs are back in their base positions); the latter is the optimization criterion defined in formula (12). In all cases, UAVs with the homogeneous flight parameters were used; the average flight velocity was set to 10 m • s .  Figure 15 shows the routes for scenarios sc03 (on the left) and sc04 (on the right). The routes of individual UAVs are color coded. Moreover, the route times (times needed for UAVs to conduct their routes) are stated. The operation ends when the last UAV returns back to its base; thus, the optimization is about a good distribution of waypoints to available UAVs. The similar times of individual routes are apparent. For the best solutions found, the routes of UAVs available in individual scenarios were planned using the algorithm mentioned in Section 4.4. The results are recorded in Table 10. The last two columns show the parameters of the optimized routes: the total distance covered by all UAVs participating in the operation, and operation time T (the time when all the UAVs are back in their base positions); the latter is the optimization criterion defined in formula (12). In all cases, UAVs with the homogeneous flight parameters were used; the average flight velocity was set to 10 m·s −1 .  Figure 15 shows the routes for scenarios sc03 (on the left) and sc04 (on the right). The routes of individual UAVs are color coded. Moreover, the route times (times needed for UAVs to conduct their routes) are stated. The operation ends when the last UAV returns back to its base; thus, the optimization is about a good distribution of waypoints to available UAVs. The similar times of individual routes are apparent. The routes for UAVs are created by connecting the waypoints by straight lines (see Figure 15). This means that the routes are easily applicable for rotary-wing aircraft with the ability of vertical take-off and landing (VTOL) and abrupt changes in direction. The only requirement for them is the ability to automatically follow a set of predefined waypoints (and, of course, it needs to be equipped with an appropriate sensor). Nowadays, this ability is common not only for the state-of-the-art military UAVs but also for ordinary commercial drones. However, the model is applicable, with some minor limitations, even for fixed-wing aircraft-see [38] for more details.

Number of Waypoints
The experiments conducted so far have assumed that the monitoring can be performed only when the UAVs are located at waypoints, i.e., not during their flight between them. In many real situations, the monitoring can be performed continuously during the flight of the UAVs along their routes. This case was tested on the best solutions found (Table 9) and the routes planned for these solutions (Table 10). The results are in Table 11 where the original coverage (monitoring only from waypoints) is compared to the so-called continuous coverage (monitoring during the flight). The continuous monitoring provides better results (the worse original solution, the higher the improvement).

Conclusions
The article presented the model of aerial reconnaissance using UAVs cooperating to conduct the operation as effectively as possible. The monitoring is performed from a number of waypoints deployed in the three-dimensional space in the area of operations. The terrain and non-transparent objects, which may cause occlusion, are taken into account. For the deployment of waypoints, the metaheuristic algorithm was proposed and verified by experiments; the results were compared with the optimal solutions. The possibility of the practical use was confirmed by a set of experiments based on the typical military reconnaissance scenarios and real geographic data. The routes for UAVs are created by connecting the waypoints by straight lines (see Figure 15). This means that the routes are easily applicable for rotary-wing aircraft with the ability of vertical take-off and landing (VTOL) and abrupt changes in direction. The only requirement for them is the ability to automatically follow a set of predefined waypoints (and, of course, it needs to be equipped with an appropriate sensor). Nowadays, this ability is common not only for the state-of-the-art military UAVs but also for ordinary commercial drones. However, the model is applicable, with some minor limitations, even for fixed-wing aircraft-see [38] for more details.
The experiments conducted so far have assumed that the monitoring can be performed only when the UAVs are located at waypoints, i.e., not during their flight between them. In many real situations, the monitoring can be performed continuously during the flight of the UAVs along their routes. This case was tested on the best solutions found (Table 9) and the routes planned for these solutions ( Table 10). The results are in Table 11 where the original coverage (monitoring only from waypoints) is compared to the so-called continuous coverage (monitoring during the flight). The continuous monitoring provides better results (the worse original solution, the higher the improvement).

Conclusions
The article presented the model of aerial reconnaissance using UAVs cooperating to conduct the operation as effectively as possible. The monitoring is performed from a number of waypoints deployed in the three-dimensional space in the area of operations. The terrain and non-transparent objects, which may cause occlusion, are taken into account. For the deployment of waypoints, the metaheuristic algorithm was proposed and verified by experiments; the results were compared with the optimal solutions. The possibility of the practical use was confirmed by a set of experiments based on the typical military reconnaissance scenarios and real geographic data.
The metaheuristic algorithm proposed for waypoints deployment is based on the simulated annealing principles. It is a simple algorithm yet very efficient in positioning problems. The strengths of the algorithm are its fast convergence to the optimum, efficient mechanism preventing a solution being stuck in some local optimum and low memory demands. The most time-consuming and memory demanding part of the whole optimization process is the evaluation of a particular solution. Therefore, a lot of effort was put into increasing its efficiency. The speed of the solution evaluation, as well as the memory demands, depends, beside the configuration of a computer used (power and available memory), on the rasterization; the number of points needed to be evaluated changes quadratically with the rasterization step. However, the precision of the evaluation is also dependent on the rasterization step. This results in two contrary requirements. On the one hand, the step should be as small as possible so that the evaluation is as precise as possible; on the other hand, it should be as large as possible so that the evaluation is as fast as possible. The choice of the rasterization step is, therefore, a compromise between the two requirements.
The model proposed in this article was implemented into the decision support system for commanders of the Czech Army. Prior to any military mission there must be planning process which is called either Troops leading procedure (in case of company and below military units) or Military decision making process (in case of battalion level and higher). This process is complex, continuous and composed of many steps. As it is apparent from the results of experiments, time of optimization did not exceed 15 min in case of the most complex problems. During this time, there are other activities which are ongoing simultaneously. Moreover, ordinary planning process with no autonomous computation must be done by military personnel and lasts for significantly longer time. In conclusion, time of computation did not influence any critical military mission at all.
Although the model is intended to be used for military purposes, there are other civil real life implementations and applications such as the rapid evaluation of a disaster area and management of recovery resources, identification of threats on land, border surveillance and protection, etc.
The future work of the authors will be aimed at further optimization of the solution evaluation process. The significant speed improvement is to be achieved by implementing the algorithm on the graphics processing unit (GPU). For the further verification of the model, the experiments using the real UAVs are planned to be carried out. The current reconnaissance model will also be extended to the persistent surveillance model, in which the continuous monitoring of the area of interest by a swarm of UAVs is assumed.

K j
Number of waypoints to be visited en route R j . T j Time to perform route R j by U j . T Duration of the reconnaissance operation. T max Initial temperature (simulated annealing parameter).

T min
Minimum temperature threshold (simulated annealing parameter). T cur Current temperature (variable used in the simulated annealing algorithm). γ Cooling coefficient (simulated annealing parameter). n 1max Number of transformations (simulated annealing parameter). n 2max Number of replacements (simulated annealing parameter). τ Constant used when estimating the necessary number of waypoints. µ Mean. σ Standard deviation. ε Constant used in the solution transformation process.