A Realistic Model to Support Rescue Operations After an Earthquake via UAVs

In this paper, we consider the problem of completely flying over an area just hit by an earthquake with a fleet of Unmanned Aerial Vehicles (UAVs) to opportunely direct rescue teams. The cooperation between UAVs ensures that the search for possible survivors can be faster and more effective than the solutions currently implemented by civil protection. To study this scenario, we introduce the Cover by Multitrips with Priorities (CMP) problem, which tries to keep into account all the main real-life issues connected to the flight and coordination of the UAVs. We conduct a theoretical study to estimate the best number of UAVs and additional batteries, to give indications to the organization that leads the rescue teams to be able to guarantee rapid and effective rescue. Finally, based on some theoretical considerations, we propose some heuristics that tackle the problem of flying over the whole area with a fleet of UAVs in the shortest possible time. Simulations show that they work efficiently in both the proposed scenarios and provide better performance than previous solutions once they are arranged to work in our scenarios. The main advantages of our approach w.r.t. the current drone-based solutions used by the civil defense are that UAVs do not need drivers so the time of all available rescue workers can be invested in doing something else. In our model, we take into account that some sites (e.g. buildings with a high fire risk or schools and hospitals) have a higher priority and must be inspected first, and the possibility that UAVs can make a decision based on what they detect. Finally, our approach allows UAVs to collaborate so that the same sites will be flown over exactly once in order to speed up the rescue mission.


I. INTRODUCTION
In case of natural disasters, such as earthquakes, rescue teams must complete their operations within a few hours of the event to avoid increased loss of life. Furthermore, in the aftermath of these disasters, we cannot count on the presence of infrastructure (e.g., there is no guarantee that internet will be available), and there may not be time to re-establish it to undertake coordinated activities. In this context, Unmanned Aerial Vehicles (UAVs) can make a difference to ensure rapid and effective aid. We recognize that the current legislation of several countries is not ready to immediately accept the implementation of solutions exploiting UAVs [1]. Nevertheless, civil defense, policies, and lawmakers are positively responding to the tremendous advantages given by innovative solutions based on the cooperation capabilities of UAVs [2].
The associate editor coordinating the review of this manuscript and approving it for publication was Halil Ersin Soken .
Researchers strongly believe that, as we have already seen for artificial intelligence and robotic systems in many fields, UAVs will gain the trust of common people and lawmakers. This is also the reason why, with recent advances in UAV technologies, many papers have been published proposing different theoretical models for handling many algorithmic optimization problems related to UAV rescue operations. Nevertheless, since the problems related to these real-life issues are very complex, most of such models appear too simplified and assume some hypotheses that are not reasonable in practice. For example, in the majority of the papers dealing with Search and Rescue (e.g., [3], [4]), UAVs are supposed not to have battery constraints at all, while these are usually rather pressing; in other cases, UAVs are assumed to be able to fly back to the base in a very short time (e.g. in [5] UAVs go to recharge their battery when they are at less than 5 minutes out of 28 of battery life left) but there is no certainty that this time will be enough; in [3], [6], [7] mobile phone signals are exploited to localize a person, but it is not certain that people keep their cell phones close to them when they are at home, especially at night; in [8], [9] UAVs send stream videos to the base station through a cellular network in mountain environments, where usually there is no guarantee of coverage; often UAVs have to fly over an entire area and not focus on certain buildings, making it impossible to assign priorities to certain places [10], [11]; finally, the costs that a UAV will face during its trip are usually considered as fully known [12]- [14], while in real-life scenarios, a UAV could spend a time even very different from the expected one to explore or traverse an area.
This paper addresses the problem of completely flying over an area just hit by an earthquake with a fleet of UAVs to opportunely direct rescue teams. The main advantages of our approach w.r.t. the current drone-based solutions used by the civil defense are that UAVs do not need drivers, so the time of all available rescue workers can be better invested; moreover, we take into account that some sites (e.g. buildings with a high fire risk or schools and hospitals) have a higher priority and must be inspected first. Finally, in order to speed up the rescue mission, our approach allows UAVs to collaborate so that the same sites will be flown over exactly once and consider the possibility that UAVs make a decision based on what they detect. This is important in order to have fast and efficient rescue operations.
More in detail, our main contributions are the following: • we define a new problem, called Cover by Multitrips with Priorities (CMP), that tries to keep into account all the main real-life issues connected to the flight and coordination of the UAVs; • we provide a complete graph model for it in two different realistic scenarios; • we make a theoretical study to provide an a priori estimate on the completion time of the rescue operations and to estimate the best number of UAVs and additional batteries that guarantee a fast and effective rescue to the organization leading the rescue teams; • we propose a mixed-integer linear programming formulation of our graph model; • based on some theoretical considerations, we propose some heuristics that address the problem of flying over the whole area with a fleet of UAVs in the shortest possible time; • we extend the algorithm proposed by Kim et al. [4], [15], which is one of the most complete solutions known in the literature, to make it work under the same operative scenarios addressed by our heuristics; we call this extension A TSPN ; • we perform extensive simulations that show that our heuristics work efficiently in both the proposed scenarios; furthermore, we experimentally compare our heuristics with A TSPN and we show that our heuristics significantly outperform it in all the key performance metrics.
This paper is organized as follows. In the next section we describe our problem and compare it with some routing problems already studied in the literature that can appear somehow similar; we detail two possible real-life scenarios; for each, we propose the associated graph model, the adopted performance metrics, and a mixed-integer linear programming formulation. Section III is devoted to the description of a meta-algorithm that we exploit to design new heuristics.
In Section IV, we consider the possibility of replacing exhausted batteries with already loaded ones: we associate each possible value of the number of additional batteries with the corresponding completion time so that the civil defense organizers can have an immediate idea of the time needed to complete the operations. Knowing all the parameters, the minimum number of batteries necessary to guarantee the minimum completion time is also obtained (i.e. we impose to have no time in which the UAVs remain inactive waiting for their battery to be recharged).
In Section V we study the possibility of concluding the overflight of the area within a single flight (that is, without the need to recharge batteries) and we propose an algorithm to exactly determine the number of necessary UAVs.
Section VI is devoted to the design of five heuristics that we experimentally compare in Section VII. Section VIII concludes the paper with a list of open problems.

II. THE PROBLEM: FROM REAL LIFE SITUATION TO MODELS
In this section, we first describe the real-life situation that we want to model in two possible scenarios, then we detail our corresponding graph problem, including a discussion on some useful optimization functions, and finally, we provide a mixed-integer linear programming formulation of it.

A. REAL LIFE SITUATION
We consider the very critical situation after a natural disaster (e.g. an earthquake) and the main problem of flying over the entire area using a fleet of UAVs equipped with sensors and/or cameras, acquiring information as soon as possible to have a clear idea of where rescue teams are most needed.
It is reasonable to assume that civil defense has an updated map of each area considered at risk and, for every relevant building, has an estimate of the information on the time needed to fly over it entirely and the level of importance. With the latter, we want to prioritize the most critical buildings, such as hospitals and schools or buildings with a high fire risk.
We consider to have a fleet of q homogeneous UAVs u 1 , . . . , u q that can work independently, taking off from the operations center v 0 . Each UAV has a battery corresponding to B ∈ R + units of flying time; when the battery is discharged, it takes R ∈ R + time units to be fully recharged, and this can only be done at v 0 . Technology currently sets 1.5B ≤ R ≤ 2B, see for example [16], so it is suitable to assume R ≥ B.
It is worth noting that, instead of a fleet of homogeneous UAVs, we could consider also heterogeneous devices, each one with a different battery constraint B i and different cruise speeds and endurance. In this case, handling parameters becomes more complicated (e.g., not necessarily battery units would correspond to units of flying time). For the sake of simplicity, we avoid this case but we discuss it in Subsection II-E how to change the MILP formulation in order to adapt it to this more general situation.
Since our scenario is a time-critical situation, we allow having a certain number b ≥ 0 of batteries that can be recharged separately while the UAVs are in flight and can replace the discharged ones in a few seconds, without causing a delay in the takeoff for the next flight. Indeed, already in 2015, Lee et al. [17] used a small robot arm to do the battery swap, the whole process taking about 60 seconds; more recently, the authors of [18], [19] propose a low-cost system only taking 10 seconds to finish the battery swapping process, from landing to take off.
We explicitly avoid handling the collision problem (which is not the focus of this document and for which there is a large literature [20]). Connected to this, we assume that the UAVs are not equipped with infrared cameras, in fact, thermal imaging cameras have proven to be reliable in detecting humans over short distances (hence low and fixed altitude is required, making it more difficult to avoid collisions), but it results in an ambiguous blob formation when used from long distances or in the presence of numerous obstacles [21].
We consider two scenarios, depending on whether the UAV fleet consists of a set of small-scale UAVs or some highcost UAVs that can be equipped with additional sensors. These represent the two extremes among the many realistic possibilities: on the one hand we assume that the rescue teams have many cheap UAVs that cannot take decisions, on the other hand, they may have a lower number of more powerful UAVs, that can both communicate with the base and take decisions.
In the first scenario, each UAV is equipped with a commercial RGB camera and has very small computing power and no communication devices. In this case, the q UAVs can only follow an assigned path through some sensible targets, take a video, and bring it back to the operations center where image detection/recognition tools will be able to detect, in parallel, possible survivors needing help in real-time. The positive side of this scenario is that the fleet can easily consist of a large number of UAVs due to their low cost.
In the second scenario, we consider high-cost UAVs, which have a processing unit with good computing power and can be equipped with a device that allows them to communicate with the base at any time (e.g. a radio); therefore they can produce and scan the videos in real-time(see e.g. [22]- [24]) in order to immediately detect and recognize if there is an emergency; if people needing help are detected, the UAVs will promptly get closer to possible survivors, to acquire more information and send a message to the base for a rescue team to be sent.
It should be noted that, due to legal constraints, both proposed scenarios are more advanced than those currently adopted by the civil defense, which either uses human-guided drones and/or allows UAVs to fly only by sight from the base.

B. RELATED PROBLEMS
Before detailing the graph model, the performance metrics, and the MILP formulation for our problem, we briefly discuss here similarities and differences between our scenarios and some routing problems already studied in the literature.
Although the problem addressed in this paper deals with limited autonomy of the vehicles due to battery capacity, it is deeply different from the routing problems involving e.g. electric vehicles, which have been broadly analysed in the literature [25]. In fact, those problems allow in-route battery recharging (or swapping) at specific locations (recharging stations) spread across the territory. Instead, in our case, batteries can be swapped only at the depot and not during the trip's execution.
A problem that could seem close to ours is the multi-trip vehicle routing problem (MTVRP) [26], in which vehicles, whose trips are not limited by the driving range but by the load capacity, may perform several trips in sequence. At the end of each trip, the vehicle returns to the depot in order to load other goods and perform the next trip. Generally, a maximum cumulative working time is imposed for each vehicle. The main difference between MTVRP and our problem concerns the objective function. In fact, while in MTVRP the goal is to minimize covered distances, in our case, we look for weighted completion time minimization. Moreover, MTVRP applications are mainly related to freight delivery (or pick-up). In such a context, the minimization of costs, which are strictly dependent on the traveled distance, is the goal of the company. In our setting, consisting in monitoring an area hit by a natural disaster (such as an earthquake), the primary objective is to serve the points of interest in the shortest possible time, considering that some sites must be flown over before than others. Indeed, most densely populated areas and some crucial buildings are marked with a higher priority because the probability of finding injured people is larger in these zones. Although completion time minimization has been already considered in the literature for a single vehicle problem [27], we are the first researchers to apply it in a multitrip and multi-vehicle context.

C. THE GRAPH MODEL
In order to model the described context, we consider a complete node-and edge-weighted graph G = (V , E, dist, σ, p), where: • V = {v 0 , v 1 , . . . , v n }, and v 1 , . . . , v n are the sites representing the points of interest inside the area, i.e., buildings (or areas) that need to be explored, while v 0 represents the operations center.
• Functions p : V → {p min , p med , p max }, with p min < p med < p max , and σ : V → R + associate to each node two values, representing the priority and the time needed to VOLUME 10, 2022 explore the site, respectively. The value of σ (v) is computed based on the shape and dimension of the target and is given as part of the input (We do not further detail how σ is computed because it is out of the scope of this paper, and simply assume that it is a known value.) In the first scenario, it represents the only (fixed) time to fly over each site. In the second scenario, σ (v) includes both the time necessary to fly over the site and to process the images; moreover, it may be necessary to add a second contribution to σ (v), i.e. the time σ (v) needed for the UAV to decrease its altitude, take more precise information concerning detected people possibly in difficulty and send to the base a radio message. We will denote by σ e (v) the effective time used by a UAV to completely fly over a site, i.e. σ e (v) = σ (v)+σ (v) and σ e (v) becomes known only after that a UAV has finished to fly over v. As far as the effective overflight time is conceived, it always holds σ (v) ≤ σ e (v); indeed, if no people needing help are detected, σ e (v) = σ (v); on the contrary, if a rescue team is deemed necessary, In both scenarios, it is not restrictive to define σ (v 0 ) = 0; indeed, even if the base had to have a not negligible area, and each UAV had to take a time s to fly over it and reach the destination, this would be equivalent to set σ (v 0 ) = 0 and battery constraint equal to B − s.
is associated a distance function dist that represents the traveling distance between two nodes.
Moreover, we assume that for dist the triangle inequality holds. For the feasibility of the problem, we need that each site is reachable and completely flown over by a UAV, i.e. for If there is a very large site v that is reachable but takes a too long time to fly over (so that 2 dist(v 0 , v)+σ (v)+σ (v) > B), we replace v with a set v 1 , . . . , v m v of m v nodes, so that the overflying time is split among the duplicates: m v for each i = 1, . . . , m v ; these m v sites are considered as different sites at null mutual distance and they keep the same set of adjacent nodes as v. The value of m v must be enough to guarantee that every v i is reachable and completely flown over by a UAV; in other words, m v is such that, for every . . , m v and it is not necessary to highlight which nodes in the input graph are of this kind. Hence, in the rest of this paper, we implicitly assume that all the edges respect the feasibility requirement.
Since we can translate the value of dist in terms of flight time, as in [28], [29], we assign to each node an edge weight function l : 2 for each i, j = 0. The triangle inequality still holds for l. In this way, from now on we will handle the only edge weight function l instead of dist and σ .
The goal of our problem is to inspect all points of interest in the shortest possible time, trying to explore sites with higher priority first; to this aim, we partition the sites into disjoint sets and assign each UAV one of such sets to fly over, so that the entire fleet can cooperate to the complete exploration, and avoid that some sites are flown over by more than one UAV.
Since UAVs are prone to their battery constraint, they may not be able to visit all assigned sites in one flight but may need to perform multiple flights, in order to either recharge or substitute their batteries. So, we have to assign to each UAV an ordered set of cycles, whose weights are bounded by B, and all nodes of the graph are flown over, giving precedence to those with higher priority, in such a way that an appropriate optimization function is minimized. Of course, in the two scenarios, the situation is different and hence also the solution is. In particular, it is worth noting that, while in the first scenario the time it takes to explore each site is exact (i.e. σ (v) = 0 ∀v), in the second one it is not (and hence σ (v) can be either positive or null). The optimization problem is not run on-board but on a dedicated computer at the operations center in order not to require the usage of wireless connectivity among UAVs.
To consider this difference and, at the same time, to give definitions suitable for both scenarios, we introduce the concepts of presumed and effective cycles.
Informally, a presumed cycle is a cycle that can be flown over in a time-bounded by B assuming that, for each node involved in the cycle, the value of σ is null; after this cycle is assigned to a UAV and it begins to fly over it, if some values of σ are not null, it may become impossible to complete the exploration of the cycle within time B; so only a proper subset of sites involved in the presumed cycle will be really flown over, and the UAV will have to go back to the base: the flown over part of the presumed cycle is the effective cycle. Trivially, in the first scenario presumed and effective cycles coincide, while in the second scenario they can be different. We now formalize these definitions and state the properties of a feasible solution and the objective functions in the two scenarios.
sites that a UAV can explore in this order in a single flight, assuming σ (v) = 0 ∀v ∈ C.
We fix time 0 at the beginning of operations for the whole system; denoting by t s (C) and t f (C) the starting and finishing times of each cycle C, and setting i 0 = i m+1 = 0, in both scenarios it holds that This inequality is not guaranteed to be true, when passing to the effective times, i.e. for We hence introduce the following: We define t s (C e ), t f (C e ) and t(C e ) analogously to is defined as an ordered set of presumed (effective) cycles assigned to UAV u i .
A set SOL (SOL e ), defined as a collection of presumed (effective) sequences, one for each UAV, is called The final result we would like to obtain is an effective solution SOL e . In particular, in the first scenario, we are able to compute a priori a set of (presumed = effective) sequences that will take part in the (presumed = effective) solution. On the contrary, in the second scenario, we will not compute all cycles of a sequence at the same time: as we will detail in Section III, we will compute a presumed cycle of a sequence only after that the corresponding UAV is back from the effective previous cycle.

D. PEFORMANCE METRICS
We want the completion time of a site to represent the time necessary to know whether there are people needing help on it. In the two scenarios, we consider a site as ''served'' by a solution SOL if different conditions are verified. Namely, in the first scenario, a node is ''served'' only after that the video of the cycle including it has been delivered to the base and the portion corresponding to it has been analyzed, while in the second scenario a node is ''served'' as soon as it has been completely flown over since UAVs can immediately detect people needing help. We formalize these facts in the following: Definition 4: Given a solution SOL e and assuming that site v k is the k-th in the effective cycle C e (belonging to a sequence of SOL e ), the completion time of v k in a solution SOL e in the first scenario is: where while in the second scenario is: where Note that the first term of cost (I ) is the time necessary to completely fly over the cycle containing v k (C e ) and the second term is the time necessary to analyze the recorded video from the beginning to the moment when v k has been VOLUME 10, 2022 flown over. On the contrary, the first term of cost (II ) is the time spent before starting flying over C e and the second term is the time necessary to fly over the portion of C e preceding v k and v k itself (including both the time necessary to fly over v k and to process the images).
We now define two functions that exploit the concept of completion time of a site and are useful to evaluate the goodness of a solution: in the first one we generalize a classical parameter (i.e. latency) in order to take into account the priority, while the second one is the classical definition of a very natural optimization function, that is completion time.
Definition 5: The weighted latency of a solution SOL e , wL(SOL e ), is the mean of the completion times of all sites (either in the first or in the second scenario), taking into account their priorities: Changing parameters p min , p med , p max influences how much it is pressing to serve higher priority sites first. Namely, sites with the highest priority should be served first, in order to associate a higher multiplicative factor (corresponding to the priority) to a lower sum of costs. Hence, we can claim the following: Problem Requirement HPSF (Higher Priority Sites First): The sites with higher priority should be flown over before than the ones with lower priority.
Definition 6: We denote as the completion time of a solution SOL as: Clearly, we consider a solution better than another one if it has smaller values of both ct and wSOL. Observe that wL (i) (SOL e ) and ct (i) (SOL e ) are not independent indeed, calling wA(p) = 1 n v∈V p(v) the weighted average of the priorities over all nodes, it holds that wL (i) Note that ct (i) does not take into account the priority of nodes; nevertheless, it is a primary parameter to evaluate the goodness of the solution.
We conclude this subsection with the formal definition of our problem: Definition 7 (CMP): Given a complete node-and edgeweighted graph G = (V , E, dist, σ, p), defined as in Subsection II-C, and a set of q homogeneous UAVs u 1 , . . . , u q , the problem Cover by Multitrips with Priorities (CMP) consists of finding a set SOL e of q effective sequences such that ct (i) (SOL e ), where i = I , II , is minimized and Requirement HPSF is addressed.
It is easy to see that CMP is a generalization of the well known multiple Travelling Salesperson Problem, m-TSP [30] (achievable with B = ∞, p(v) = 1 and σ (v) = 0 ∀v) so, in its turn, it is strongly NP-hard. Nevertheless, the wide literature devoted to m-TSP and all its variants cannot be exploited in this case; indeed, this class of problems misses an important parameter, that is the battery constraint B. Modifying algorithms for m-TSP to address the battery constraint does not seem an easy issue as pointed out in [28] for a different problem, and it seems neither possible to inherit from m-TSP any approximability result without relaxing the battery constraint (and accepting, for instances, cycles with completion time B + for some > 0), that in our context is inadmissible.

E. A MIXED-INTEGER LINEAR PROGRAMMING FORMULATION
We exploit the nomenclature defined in Subsection II-C to formulate CMP as a mixed-integer linear programming problem in the first scenario.
The following sets and parameters are known from the problem: (i) V : set of n + 1 sites; (ii) Q: set of q UAVs; (iii) B: UAVs' battery capacity; (iv) p(v i ): priority of site v i .
Moreover, we define the following: (i) K : set of feasible (effective) cycles and K = K ∪ {k 0 } where k 0 is an empty cycle introduced for modeling reasons; for all C k ∈ K , t k (v i ) (completion time for site v i on cycle C k ) and t(C k ) (flying time of cycle C k ) are considered as known; (ii) Cov i : set of cycles passing through site i; (iii) : a very small constant; (iv) M : a very large constant.
We introduce the following decision variables: is a binary variable assuming value equal to 1 if cycle C k is executed by UAV j; (iii) Z j k 1 k 2 ∈ {0, 1}, k 1 , k 2 = 0, . . . , |K | and j = 1, . . . , q, is a binary variable assuming value equal to 1 if cycle C k 1 and cycle C k 2 are executed in sequence by UAV j and 0 otherwise; (iv) t s (C k ), k = 0, . . . , |K |, is a nonnegative variable representing starting time of cycle C k ; (v) τ is a non-negative variable representing the total completion time; (vi) cost(v i ) is a non-negative variable representing latency of site i.
Note that t s (C k ) and cost(v i ) are not independent variables, on which we can directly act, as they are expressed as a function of other decision variables. Nevertheless, in a mathematical programming model, they must be defined as decision variables as well. Indeed, in an integer programming model, two types of entities are considered: input data, which are fixed and cannot be modified by the model, and decision variables, whose value must be decided by the model. The latter group includes both independent variables (such as those related to the selection and assignment of cycles) and variables that are expressed as a function of them (such as the total cost or the starting time of a cycle).
We claim that, in the first scenario, CMP can be completely described by the following set of constraints: The hierarchic objective function is reported in Equation (of). This function firstly aims at minimizing the total completion time and secondly at minimizing targets' weighted latency. The very small constant is defined such that a solution with a lower total completion time would be always preferred w.r.t. one with a higher one, whichever is the value of the secondary objective. This is a standard technique commonly used in multi-objective optimization dealing with hierarchical objective functions [31]- [33].
Constraints (C1) guarantee that each target is covered by exactly one cycle while Constraints (C2) impose that, if a cycle is selected, it must be assigned to exactly one UAV. Similarly, a cycle must occupy a position in the schedule list for a given UAV j if and only if it has been assigned to that UAV (Constraints (C3)). Constraints (C4) delineate the sequence of cycles for each UAV. If variable Z j 0k is equal to 1 it means that k is the first cycle executed by UAV j. Similarly, if variable Z j k0 is equal to 1, k is the last cycle executed by j. Constraints (C5) imply that each UAV is used and explicates that, for each UAV, there is a cycle that is flown over as the first one and a cycle that is flown over as the last one. These two cycles can coincide in the case of single-cycle UAV routes.
Starting times of cycles are ruled by Constraints (C6) and (C7): a cycle can start only after the previous cycles assigned to the same UAV has been completely flown over. Note that the order relationship between the starting times of the different cycles k 1 and k 2 is given by the term M (1 − j∈Q Z j k 1 k 2 ), that guarantees the correctness of this constraint. Indeed, if C k 1 and C k 2 are executed in sequence by the same UAV, then Z j k 1 k 2 = 1 and thus we have t s (C k 2 ) ≥ t s (C k 1 ) + t(C k 1 ), thus guaranteeing that cycle C k 2 can only starts after the end of the previous cycle. On the other hand, if C k 1 and C k 2 are either executed by different UAVs or are not consecutive, we have that Z j k 1 k 2 = 0 and thus t s (C k 2 ) ≥ t s (C k 1 )+t(C k 1 )−M , making this constraint not binding.
The total completion time is computed through Constraint (C8). Finally, Constraints (C9) allow one to identify the latency of each site.
Note that M is a big constant, large enough to make Constraints (C6) actually binding only if cycle k 2 is performed by UAV j just after cycle k 1 . A similar role is played by M also in Constraints (C8), which are actually binding only if cycle k is performed, i.e. only if X k = 1.
The model involves |K |+|K |q+|K | 2 q binary variables and |K |+n+1 continuous variables. The number of constraints is n+|K |+2|K |q+q+|K | 2 +n|K |+1. If we want to model the second scenario, we can run the following recovery procedure: we solve the model as if we were in the first scenario and send all the UAVs following the optimal solution one cycle at the time; whenever a UAV is not able to fly over all the sites assigned to it, we construct a reduced graph starting from the input one and removing all the sites effectively visited together with their incident edges; we solve again the model on this reduced graph and send again all the UAVs following the new solution.
Finally, the model can be adapted to the case of a heterogeneous fleet of UAVs, characterized by a different battery duration and consumption rate. Indeed, feasible trips can be generated considering the less binding type of UAV, i.e. the one performing the longest duration trips. Then a set of constraints can be added in order to allow us to assign a cycle j to a UAV k only if the battery of k is enough to completely fly over j. We can define the compatibility between cycle j and UAV k as an input constant φ kj , equal to 1 if UAV k can perform cycle j and equal to 0 otherwise. Then, to avoid forbidden assignments, we add the following set of constraints to the model:

III. A META-ALGORITHM FOR CMP
In this paper, we will propose several heuristics, each one based on a different consideration. Nevertheless, all of them are built upon on the same meta-algorithm that is run at the operations center; it consists of a number of iterations going on until all sites have been flown over; at each iteration, a presumed cycle for each UAV is produced as output.
Here we describe an outline of this meta-algorithm, through some facts on which it is based: • At the beginning, v 0 has complete knowledge of the structure of graph G, including the distances between the sites and the presumed values σ (v), but (only in the second scenario) not the effective values σ e (v).
• The solution is computed in any case at v 0 , in order not to overburden UAVs with computation and intercommunication issues. This is possible because UAVs require to come back to v 0 to either recharge or substitute their battery and there they can get new information, e.g. a new cycle to perform.
• The meta-algorithm consists of several iterations to be executed until all sites have been flown over. At each VOLUME 10, 2022 iteration, a certain algorithm A is run; A computes a set of (presumed) cycles, one for each UAV, and for each such cycle C, t(C) ≤ B. Changing A, either in one or in all the iterations, produces a different heuristic, as we will detail later in Section VI.
• In the second scenario, each UAV is equipped with a device (e.g. a radio) to directly communicate with the base. Such a device is primarily exploited to call rescue teams when it realizes that there are survivors to be rescued. Moreover, as soon as a UAV decides to go back to v 0 , for example, because it has not had enough battery to continue, it sends a message indicating its latter overflown site, and v 0 will immediately update the set of sites to be still flown over. In this way, v 0 has in real-time a precise knowledge about the (im)possibility to successfully complete the flown over of the assigned cycle by each UAV.
• In the first scenario, presumed and effective cycles coincide, so the iterations of the meta-algorithm can be run all together, outputting a final solution.
In the second scenario, presumed and effective cycles can be even very different, so the iterations are performed one by one, and the next iteration is run only after that UAVs inform the operations center v 0 of the effective cycles. Namely, at each iteration, a cycle for each UAV is computed at v 0 and it is assigned as a presumed cycle; it is assumed that all the sites in these assigned cycles will be flown over; therefore sites currently assigned to a UAV are considered committed and cannot be assigned to any other UAV. Based on the effective value of σ on the already flown over sites in the current cycle and of the remaining battery, every UAV can periodically infer whether some sites assigned to it will remain not flown over. Specifically, before visiting the next node, a UAV is able to calculate if the remaining energy is not enough to continue the trip and thus it has to return to v 0 . All the remaining nodes in the presumed cycle will be automatically discarded from the current effective cycle. The UAV communicates this information to v 0 , where the not flown over sites will be released and included again in the current graph; then, a new set of presumed cycles is computed. In this way, as soon as each UAV comes back to the base, it will have ready a new presumed cycle to be assigned to it. This iteration continues until all the sites have been flown over.
In the following sections, we will show that this metaalgorithm can also be used as a preprocessing algorithm to make some preliminary estimates of the required completion time and the best number of additional batteries and UAVs required.

IV. CHANGING VS RECHARGING BATTERIES
Suppose that b ≥ 0 additional batteries are available for the q UAVs. Assume that the meta-algorithm is run immediately after an earthquake as a pre-processing and a (presumed) solution SOL = {S u 1 , · · · , S u q } is output; we call N (SOL) = q i=1 |S u i | the total number of cycles in SOL. As SOL is clear from the context, we call it simply N .
In all this section we assume that max u i |S u i | > 1, otherwise there were no need of additional batteries, and S u ct = {C 1 (u ct ), . . . , C y ct (u ct )} is the sequence flown over by UAV u ct whose cost equals the completion time of SOL, and its last site before coming back to v 0 is v ct ; in other words, C y ct (u ct ) is the last cycle to be oversight in the sequence assigned to u ct and v ct is the last site lying on it. (Note that sequence S u ct does not necessarily coincide with the one with the maximum number of cycles.) We now analyze each possible value of b and associate to it the corresponding value of the completion time, so that civil defense organizers can have an immediate idea about the necessary time to conclude operations, knowing all the parameters. We derive also the minimum number of batteries b necessary to guarantee minimum completion time (i.e. no time in which UAVs remain inactive waiting for battery recharging).
We focus first on the first scenario and, for the sake of simplicity in calculations, we assume that every cycle of each sequence of every UAV is flown over in time exactly B. This is an upper bound on the completion time of cycles, but this approximation is, in fact, negligible especially for some heuristics, as shown in Table 1, so we use it to simplify our calculations. It is not difficult to see that more precise values give anyway similar results. We distinguish different cases according to the relation between b and q. case b = 0: If there are no additional batteries, the completion time is trivially given by Informally, we can read this result as follows: we consider time as marked in intervals B + R long, and an immediate upper bound on the maximum time necessary to complete operations by solution SOL is given by the product between (B + R) and y ct ; in view of the definition of completion time, a more precise approximation substitutes the last R long interval with t (I ) (v ct ) that is bounded by B. As a side effect, this calculation proves also the following finitetime result: Theorem 1 (Termination 1 st sc): In the first scenario, there always exists a solution to CMP and an upper bound to its completion time is given by (y ct − 1)(B + R) + 2B (reached in absence of additional batteries and when t(C) = B for every C in SOL).
If there are few additional batteries (where ''few'' means that they are not enough to cover all the operations without UAVs' inactivity time due to recharging batteries), b is necessarily less than: • N − q, that is the total number of flights in SOL after the subtraction of the first one for each UAV, performed with the supplied battery; • R B q, because more batteries would imply no UAVs' inactivity due to recharge operations. Observe that the value of ct (I ) (SOL) as function of b is not decreasing (in other words, for any k > 0, the completion time using b + k additional batteries is not greater than the completion time using b additional batteries). We will use the monotonicity of this function to simplify calculations and consider only the cases when b is a multiple of q. In this special case, the whole fleet of UAVs will be able to go on and back exactly b q + 1 times without waiting for recharging batteries. Every time the q UAVs come back to the base, exactly q batteries are put to be recharged. After b q +1 flights, all UAVs (and among them u ct , so it is not restrictive to focus on it) necessarily have to wait for the first batteries to be completely recharged; in this case, the waiting time of u ct due to battery recharge after cycle C b q +1 (u ct ) represents the waiting time after the first group of flights and is wtb (1 i=2 t(C i (u ct )). After this time, u ct (and analogously all the other UAVs) will be able to perform another group of b q + 1 flights without any waiting time. More in general, the waiting time after the r-th group of flights is wtb(r) = The number of groups of flights that need to be followed by a waiting time are y ct b q +1 − 1, so the completion time is It is not difficult to see that, with so many additive batteries, it is possible to let UAVs fly without pauses, because there are always enough charged batteries. In this case, the completion time is upper bounded by ct (I ) This proves the following result: Theorem 2 (Completion Time Minimization): In the first scenario, a guarantee for minimizing the completion time is to have at least R B · q additional batteries and to minimize the value of y ct . Figure 2 depicts the comparison between the function describing completion time in the first scenario when b grows up from 0 to R B · q on an example instance and the upper bounds provided above. The gap between the two plots in the interval between 0 and R B q is due to the approximation we did upper bounding ct (I ) (SOL) with its values reached when q divides b; this is also the reason why the corresponding function is a step function. Nevertheless, our computation is quite good when we need to estimate the number of needed additional batteries to avoid waiting times due to recharging operations.
All the previous considerations cannot be applied to the second scenario as they are, because the pre-processing phase outputs only presumed cycles. In the following, we estimate the average of an additional contribution that cannot be ignored when passing from a presumed to an effective solution.
Assume that, immediately after the earthquake, we are able to approximately quantify a uniform probability p to find people needing help during the flight and that, for simplicity, σ (v) is a constant σ and σ (v) is either 0 or a constant σ , for every v ∈ V . Given a presumed solution SOL, with N presumed cycles output by a pre-processing computation, the average of the number of sites in each cycle is E[m] = n N and the E[m] + 1 sites lie on every cycle at an average distance of +1 . For every presumed cycle C of SOL, on average, in E[m] · p sites a UAV will detect a need to spend a further time σ long.
When passing from C to the corresponding effective cycle C e (containing m ≤ m sites), m − m sites will not be in fact flown over, and the average value for m , E[m ], is the maximum value for which it holds: Executing some algebraic calculations, we get: Now, calling N (SOL, p, σ, σ ) = pσ (n+N ) B+σ , the total number of estimated cycles is at least N + N (SOL, p, σ, σ ). If we substitute this value instead of N in the computations related to the first scenario and modify them according to the definition of completion time of a site in the second scenario, we have an estimate on ct (II ) (SOL) w.r.t. the number of batteries b. We explicitly write only the following, that is the analogous of Theorem 1 and provides also a finite-time result because Theorem 2 can be re-written identical for the second scenario: Theorem 3 (Termination 2 nd sc): In the second scenario, there always exists a solution to CMP and an estimate to its completion time is given by (N+ N (SOL, p, σ, σ ))/q (B + R) − R (reached in absence of additional batteries and when t(C)=B for every C ∈ SOL).
In our experimental evaluation we observed that, although the derived formulas would provide lower bounds to the estimates of ct (II ) (SOL), in practice they are quite close to the exact values, especially for small values of p, as shown in Table 2.

V. NUMBER OF UAVs TO GUARANTEE A SINGLE CYCLE
In this section we give a simple method to immediately estimate the number of UAVs that civil defense should own to ensure that all the operations can be concluded within a single UAV flight, i.e. without any changing or recharging batteries and within time B. To do this, we recall the definition of a well-known problem in the literature.
The Team Orienteering Problem (TOP) [34] comes from an outdoor sport played on mountains where a team of hikers, with the help of a map and a compass, must visit as many nodes as possible within a specified time limit, getting a profit from each visited node. Since the goal is to maximize the total gain that the hikers all together collect, if the gains associated with nodes are equal, equivalently the problem is to maximize the number of nodes visited. TOP is also known as selective m-TSP since not all nodes are visited in a solution.
More in detail, using our nomenclature, the Team Orienteering Problem (TOP) [34] can be defined as follows: let G = (V , E, l, p) be a complete graph, where V = {v 0 , v 1 , . . . , v n } is the set of nodes, v 0 is the base and the remaining nodes are objectives. Each node v i , i = 0, 1, . . . , n, is associated with a profit p(v i ) (p(v 0 ) = 0), corresponding to our priority. There are q hikers and each one gains a profit p(v i ) if visits a still unvisited v i . A travel time l(v i , v j ) is associated to each edge (v i , v j ), that is assumed to satisfy the triangular inequality. The hikers must complete their tour within a predetermined time B. It follows that some nodes may not be visited due to the B constraint, so TOP consists in determining a set of q cycles, each passing through v 0 and respecting the time constraint such that each node is visited at most once and the total profit collected is maximized. If the profit of all sites is the same, TOP maximizes the number of visited sites.
It is evident that, if after solving TOP on our instance, all the sites are visited, we are guaranteed that a single flight (i.e. a time-bounded by B) will be enough to complete all the rescue operations. We want to give an estimate on the number q of UAVs needed to make this possible.
Let us focus on the first scenario, and first put forward the hypothesis to have n UAVs. Because of the feasibility of the problem, l(v 0 , v) + l(v, v 0 ) ≤ B, so we are guaranteed that each UAV will fly over a single site, concluding the whole operations with a single flight and within time B. Now, we can insert the computation of an algorithm solving TOP on our instance inside a loop that mimics a binary search on the value of q; so, we assume now that q = n/2, run the algorithm and output a solution for TOP; if all the sites are visited by this solution, then we eliminate all the values larger than current q, otherwise, we eliminate all the values smaller than current q, and then we will proceed again recursively with a new q equal to the median of the remaining values. This algorithm, together with the fact that presumed and effective cycles coincide in the first scenario, guarantees that, in a number of iterations logarithmic in n, we have the exact value of q allowing us to get a solution in which all UAVs fly for a single cycle.
Concerning the second scenario, we have no knowledge in advance about the exact number of cycles in an effective solution, so we can only estimate this number similarly to how we did in the previous section adding the contribution of N (SOL, p, σ, σ ) and mime again a binary search. In this case, we will get an expected value of q.
We formalize the result of this section as follows: Unfortunately, TOP is NP-hard and APX -hard [35]. So, there is a wide literature solving the problem either optimally or with good approximations (see e.g. the survey [36]). Nevertheless, even these latter algorithms require very long times, though theoretically polynomial. Thus, we use a very simple and fast heuristic for TOP [37], which guarantees an immediate calculation of the cycles, at cost of a slightly worse solution; indeed the obtained results of this heuristic have a similar quality as the results of the best-known heuristics while the computational time is reduced significantly. Of course, the better the guaranteed approximation of the algorithm used to solve TOP, the more accurate the number of UAVs produced by the algorithm described above will be.

VI. HEURISTICS
In this section, we present five heuristics to solve our problem both in the first and in the second scenario. Every time each of them is executed, it produces as output q node-disjoint cycles (except v 0 , from which all of them pass), C(u 1 ), . . . C(u q ) with the properties that t(C(u i )) ≤ B and Problem Requirement HPSF is somehow addressed.
Once more we observe that the meta-algorithm presented in Section III produces a feasible and definitive solution in the first scenario, after which the algorithm A is defined. As for the second scenario, the meta-algorithm is run in iterations so that the cycles of the next iteration are computed only after that it is known which portion of the cycles of the previous iteration has been in fact flown over and the input graph of A at each iteration is the subgraph induced by the sites that are still uncovered.

A. TOP BASED ALGORITHM A TOP
Before describing algorithm A TOP , we present a case study from which it is possible to derive a general rule that we would like to follow to get better solutions.
Assume to have an instance (either of the first or of the second scenario) in which many nodes v 1 , v 2 , . . . v n−1 lie very close to v 0 , while a single node v n is far from it; moreover, due to the distances between v 0 and all the other nodes, it is possible to group the sites in only two cycles: For simplicity, we assume that q = 1, that the n nodes have all the same priority (w.l.o.g. equal to 1) and that σ (v i ) = σ e (v i ) for each i = 1, . . . , n. There are only two possibilities for the solutions: C 1 is executed either as first (SOL 1 ) or as a second cycle (SOL 2 ). It is easy to see that ct(SOL 1 ) = ct(SOL 2 ) while wL(SOL 1 ) is much lower than wL(SOL 2 ), indeed in SOL 2 , the time contribution of flying over the first cycle followed by the possible recharge time turns out to be multiplied for the number of sites flown over in the successive cycle. More precisely, in the first scenario: n · wL (I ) (SOL 1 )= n · t(C 1 ) + R + t(C 2 ) + n i=1 t(v i ) and n · wL (I ) (SOL 2 ) = n · t(C 2 )+(n−1)(R + t(C 1 ))+ n i=1 t(v i ) while in the second scenario: n · wL (II ) (SOL 1 This is a very special example with many simplifications, but the result can be generalized as follows: The cycles with a greater number of sites should be ante-posed to the other cycles in a solution. We observe that TOP (defined in Section V) aims at maximizing the profit, which results either in choosing high priority sites or in choosing a larger number of low priority sites; both cases go in favor of the minimization of wL.
We select as A TOP the heuristic described in [37], having time complexity O(|E| log |V | + q|V | 2 ) = O(n 2 (q + log n)): indeed, to the best of our knowledge, every approximation algorithm for TOP has a very high computational time in practice and so cannot be used; instead the chosen heuristic turns out to have a similar quality as the results of the bestknown approximation algorithms while the computational time is reduced significantly.

B. GREEDY BASED ALGORITHM A Greedy
This algorithm finds all the q cycles at the same time according to a greedy approach. Namely, it first fixes an a priori order on the UAVs, then, at each step it considers one by one all the q current cycles and, for each of them, starting from the site selected last v last (at the beginning v 0 ), selects the next site v next as the one maximizing the quotient p(v next )/l(v last , v next ) still guaranteeing that the whole cycle with the addition of v next can be flown over within time B.
The reason why we decided to maximize p(v next )/ l(v last , v next ) is that we would like to choose sites with high priority (hence addressing rule HPSF) without going too far, whenever possible, to keep low the necessary number of cycles.
Calling s the number of sites involved overall in the q computed cycles, the computational complexity of this heuristic is (ns) = O(n 2 ).

C. IP FORMULATION BASED ALGORITHM A IP
We observe that, although elegant, neither for very small values of the number of sites the model described in Subsection II-E can be used in practice, because |K | is exponential in the number of sites, and the number of Constraints (C3) and (C4) is quadratic in |K |. This is the reason why we propose a heuristic based on a modification of this model. VOLUME 10, 2022 Namely, we divide the problem into two consecutive subproblems: first, we compute an optimal cycle cover of all the sites (where all cycles pass through v 0 ): the optimization function minimizes the number of cycles in the solution; secondly, we perform an optimal scheduling of the only cycles in this solution; in this case, we minimize the completion time.
Clearly, the solution to this variation is possibly worse than the optimal one. Since computing an optimal cycle cover is still a time-consuming issue, to speed up the computational time, even by deteriorating the quality of the solution, we do not generate the whole set of possible cycles K but a subset of it. Namely, when we are generating a cycle, we impose that every site can be followed in the cycle by one among the k sites closest to it, where k is a parameter to be manually tuned. Clearly, the largest is k, the larger will be the number of generated cycles, the closer the solution will be to the optimal one, but the larger will be the required computational time.
To reach a good balance between efficiency and effectiveness, in our experiments we choose k = 5.
The algorithm consisting first of (not optimally) solving the cycle cover problem and then in (optimally) solving the scheduling problem on the solution of the first sub-problem represents algorithm A IP .
Solving the cycle cover phase in this simplified model corresponds to the minimization of the cumulative length of the selected cycles, subject uniquely to Constraints (C1). This model is much simpler to solve with respect to the complete model since it involves only |K | binary variables, no continuous variables, and n constraints. In the scheduling phase, we apply the model proposed in Section II-E but on a restricted set of cycles, K containing only those selected by the simplified model adopted in the first phase. The cardinality of K can be several orders of magnitude lower than the cardinality of K making this subproblem much faster to be solved with respect to the original one.

D. ALGORITHM A TSPN
As discussed in the Introduction, no prior works are dealing with the scenarios we consider in this article. So, in order to compare the performance of our heuristics, we modified one of the best-performing algorithms in the literature, that is, one of the algorithms described in [4], [15]. In their work, Kim et al. designed many algorithms for solving different problems, all related to ours, but no one of them is the same. Among them, we selected the closest to ours, that is what they call q-TSPN (q-Travelling Salesman Problem with Neighborhood). The algorithm designed by the authors for this problem is one of the best solutions currently available. Indeed, unlike several previous proposals, it is not a heuristic but it is an approximation algorithm with a provable constant approximation ratio. Nevertheless, the considered problem requires weaker assumptions than ours, i.e. no battery constraints and no priorities for sites, so we need to modify their algorithm to compare it with our heuristics.
Note that the modifications we made to the original algorithm do not compromise its performance in any aspect, but  are only meant to extend its applicability to our two scenarios: indeed, running our modified algorithm with a very high value of B and with all equal priorities produces exactly the same output as the original algorithm. The original algorithm gives as part of the input q sites and uses them as children of the root of a minimum spanning tree rooted at v 0 ; the q sub-trees of this minimum spanning tree are then transformed into q cycles covering all sites and intersecting only in the root using the Christofides's approximation algorithm for TSP [38]; finally, some operations are executed in order to equalize the weight. We do not have in input the q sites, so we select them as a solution of metric q-center problem [39]: given set V with n sites and distance function l, the goal is to find a subset C ⊆ V with |C| = q such that the maximum distance of any point in V to the closest point in C is minimized. Unfortunately, the problem is NP-hard, and the algorithm we use guarantees an approximation ratio of 2 with a time complexity of O(nq).
Moreover, to force this algorithm to take into account the site priorities, instead of working on the entire input graph, we consider it as the union of subgraphs, each one induced by the nodes with the same priority plus v 0 : G = G min ∪ G med ∪G max , where G i is the subgraph induced by {v 0 }∪{v ∈ V | p(v) = p i } and i = min, med, max. The aim is to choose first as many nodes with high priority as possible and to choose nodes with lower priority only later while considering the consumed budget.
Recalling that G is the graph induced by the set of sites that are still uncovered, this can be done by selecting the nodes  to be added to the sub-trees from G max first; after choosing as many nodes as possible in this graph, if the computed sub-trees have not been completed yet, the nodes will be selected from G med and finally from G min ; of course, after that, the meta-algorithm has executed some times A TSPN , G max will remain empty, hence, only nodes from the other graphs will be selected. When passing from the sub-trees to the cycles, the priorities are again kept into consideration, and the sites with high priority are placed along with the cycles before the sites with a medium priority which, in turn, are placed before the sites with low priority. This is necessary especially in the second scenario, in which there is no guarantee that all the sites in a cycle will be served in that cycle and, in this way, sites with lower priorities are possibly postponed to the next cycle. We call this algorithm A TSPN .
Observing that the update of the considered subgraph does not affect the worst-case time complexity (for which the dimension of the input is given by O(n) nodes and O(n 2 ) edges), the total cost of this algorithm is O(nq) to find the centers, plus the cost of finding the approximate TSP, O(n 2 log n).
Before concluding this subsection, it is worth noting that, although the resulting algorithm is strongly inspired by [4], [15], we cannot guarantee anymore that the computed spanning tree is minimum; the reason is that we run this algorithm first on G max , then on G med and finally on G min . It follows that its nice approximation ratio cannot be inherited in our scenarios.

E. HEURISTICS
We exploit algorithms A TSPN , A TOP , A Greedy and A IP to design five different heuristics, called H TSPN , H TOP , H Greedy , H mixed and H IP : consider the meta-algorithm described in Section III; if A is equal to A TSPN at each iteration, we get H TSPN ; if A is equal to A TOP at each iteration, we get H TOP and if A is equal to A Greedy at each iteration, we get H Greedy . Then, observe that TOP is not designed to be iterated, and so it gives its best during the first iteration, choosing first sites that are as close as possible to v 0 to maximize their number; this means that, while in the first iteration we select very attractive cycles, along with the successive iterations the solution of TOP could degrade in terms of goodness. This is the reason why we introduce a fourth heuristic, H mixed , running A TOP in the first iteration of the meta-algorithm and A Greedy in the successive ones.
Finally, in the first scenario, H IP coincides with A IP and we do not need to exploit the meta-algorithm. Vice-versa, in the second scenario, H IP is obtained by running the metaalgorithm with A IP if not at each iteration, at least each time we find a cycle for which presumed and effective versions do not coincide.
The difference of this algorithm w.r.t. the other ones is that, in this case, at each iteration of the meta-algorithm, a whole solution is output. This seems to be unavoidable for an approach based on the mixed-integer linear programming model and further contributes to increasing the computational time of this heuristic.
In Figure 1, a list of symbols used along the paper is reported. In the next section, we will show the result of several experiments geared to compare the performance of these five heuristics.

VII. EXPERIMENTS
All our experiments have been performed on a computer equipped with an Intel(R) Core(TM) i5-7200U CPU (4 cores clocked at 2.5GHz) and 8GB RAM; our programs have been implemented in C++ (g++ compiler v9.2.1 with optimization level O3).
We evaluate the performance of the proposed heuristics on two types of randomly generated graphs: we position n ∈ [50, 200] points, placed: (i) either uniformly at random;  (ii) or according to a Poisson point distributed method in a unit square modified as follows: we divide the unit square into 4 × 4 sub-areas and position the points in each sub-area using a Poisson point process. This distribution is intended to simulate building collapses concentrated in some specific zones. Figure 3 depicts a graphical representation of n = 100 sites placed in a squared area when the two probability distributions are applied.
In all the experiments, we use the following parameter setting: B = 50 minutes (in agreement with the current technology, see e.g. [16]) and -to ease the reading of the results-the speed is chosen as 1000 meters per minute, corresponding to about 17 meters per second (in line with other recent papers in the literature, such as [40], [41]); σ (v) is a real value randomly set (with uniform distribution) in the interval (0, 3]; when simulating the second scenario, we assign to sites a not null value of σ randomly chosen with uniform distribution either in the interval (0, 3] or in the interval (7,10] with probability p, p varying among three possible values: 0.25, 0.5 and 0.75. Only when σ ∈ (7, 10], we need to set B = 60, instead of 50, to guarantee that there exists a feasible solution. We choose a squared area of interest with dimension 15km × 15km, which always guarantees the feasibility of the problem, positioning v 0 at its bottom-left corner. Finally, we consider different values of q, coherently with the two scenarios: in the first one, we set q = 5, 10, 15, 20 while in the second one q = 2, 4, 5, 10.
In all the experiments dealing with the second scenario, the results are very similar with different values of p and q so we show the plots of only one setting, namely the one with q = 10 and p = 0.25.
All simulation runs have been repeated 20 times with different seeds, the plots show mean values with 95% confidence intervals.
First, we studied the running time of all the heuristics; as shown in Fig. 4 they do not behave all in the same way (in particular, we had to plot the running time of H IP separately for the sake of clarity because it is significantly higher than the other ones); nevertheless, the times are all low enough to consider all the heuristics equally suitable to be run in an emergency situation such as the one proposed in this paper. We remark that an approximated algorithm for TOP instead of A T OP would have required a time several orders of magnitude larger, and this would have been realistically inapplicable.
Then we studied the two performance metrics completion time and weighted latency.
As shown in Fig. 5 and 6, for all the heuristics the completion time grows almost linearly with the number of sites, but the growth rates of H IP , H Greedy and H mixed are lower than the other ones. Note that changing the values of parameters q and σ 6 does not affect the trend of the completion times. The behavior similar to a step function of the completion time in some cases can be perhaps explained with the generation method of the graph instances (any graph with n nodes is constructed by an (n − 1) node instance and so contains it as a sub-graph).
On the contrary, the weighted latency, see Fig. 7, and 8, does not always show a linear growth, due to its definition and to Problem Requirement HPSF, for which certain algorithms keep sites with high priority in the first cycles. Also for this function, H IP and H Greedy behave better than the other heuristics, and also, in this case, there is no valuable difference changing the value of q and of σ . Although, note that when increasing the value of σ , H IP tends to behave slightly worse w.r.t. H Greedy , especially when using the poisson distribution.
Finally, to highlight which heuristics best address HPSF, and to understand the different behavior of the heuristics, we considered the time in which all sites with a certain priority have been served (see Fig. 9). Observe that H TSPN behaves well in completely flying over the nodes of G max but is then the last one to complete the operations (i.e. to completely fly over G min ); this is due to the adjustment of the algorithm from [4], [15] which forces it to process the high priority nodes first. H IP and H Greedy have the lowest completion times on all three subgraphs because both heuristics explicitly require site priorities to be taken into account. Note that the completion time of G min does not coincide with the total completion time, since the time required to return to the base is not considered in Fig. 9(c).
In general, H TOP is one of the least performing heuristics, for two different reasons: first, we recall that we are using the solution output by a heuristic since algorithms guaranteeing provably good solutions are too slow and cannot be used; second, TOP aims at maximizing the profit, and no importance is given to cycle equalization. This is, instead, an issue addressed by H Greedy , whose philosophy consists of filling up q cycles together. This is the reason why this algorithm performs better than all the others, besides being very simple and particularly fast.
We conclude this section by highlighting that H IP heuristic performs better than the others from the point of view of performance metrics even if it is not far from H Greedy . On the other hand, H IP is by far the slowest heuristic (3 orders of magnitude more than the others, which have similar behavior). So, we conclude that H Greedy is the right compromise that balances all the metrics considered: its simplicity, low execution time, and good performance make it the best heuristic among those considered.

VIII. CONCLUSION AND OPEN PROBLEMS
In this paper we introduced a new problem, called Cover by Multitrips with Priorities, trying to model most of the parameters that emerge from the real-life situation of detecting people needing help immediately after an earthquake thanks to a fleet of UAVs.
We studied the problem theoretically, providing a priori estimates of the number of additional batteries to have waitings due to recharging batteries and the number of UAVs Then, we have proposed some heuristics on which we have carried out extensive experimentation, comparing them with each other and with one of the best-known algorithms that solve a similar problem; the results show that our heuristics significantly outperform this algorithm, although we added to it some parts to let it work well in our settings.
Many issues and possible future improvements naturally arise from this paper, and we already started to tackle some of them. Here is a non-exhaustive list: • We considered a fleet of homogeneous UAVs; if they had different batteries and different cruise speeds and endurance, an easy approach would be to substitute in every computation the value of B with either min i B i (but we would get a worse performance of all the heuristics) or max i B i (but not all UAVs could terminate every cycle). An interesting open problem consists in modifying all the heuristics in order to let them work efficiently in the case of heterogeneous UAVs.
• We considered static distances; instead, distances could change in time in order to keep into account the weather and wind conditions. In this case, we can consider the graph with dynamic edge weights.
• As stated in Subsection VI-C, when we run heuristic H IP , we allow each site to be followed by k = 5 other sites and so generate cycles. We wondered how much the VOLUME 10, 2022 goodness of the results is affected by this parameter. To this aim, we performed the experiments depicted in Fig. 10, showing that increasing the value of k significantly improves the performance metrics. Nevertheless, we expect that, for a certain value of k on, there will no longer be an appreciable improvement in performance metrics, that these values are very tight upper bounds for optimal values, and that we could hence consider them as benchmarks for testing heuristics and approximation algorithms.
Unfortunately, these limit values are not easy to achieve because both the execution time and the required RAM grow exponentially with k and this makes it unfeasible for any experiment with larger values of k.
• When we compute cycles of a solution we go on until there are sites that can be added and that guarantee to remain within the battery constraint B; so a UAV can return to base with a not-negligible amount of remaining battery that goes lost. We wondered if it was possible to accept that a site is partially flown over during a certain cycle to completely drain the battery of the corresponding UAV and then its exploration is completed at a later time (possibly by a different UAV). We decided not to implement this possibility since we have some case studies showing that this approach can be either favorable or not, depending on the relation between the main parameters of the problem. We are trying to give conditions on the parameters to delimit when partial overflights are convenient.
• From a methodological point of view, future research could address the design and development of other effective heuristic methods, such as Greedy Randomized Adaptive Search (GRASP) [42] or Multistart Random Constructive Heuristic (MRCH) [43] to generate potentially good cycles to be exploited by the IP model in heuristic H IP .
• Complete exploitation of UAVs' potentialities requires to introduce cooperation, both among them and between people and UAVs; these interactions would introduce new possibilities of producing faster and more flexible solutions.