Brief communication : Vehicle routing problem and UAV application in the post-earthquake scenario

. In this paper we simulate a Unmanned Aerial Vehicle’s (UAV) recognition after a possible case of diffuse damage after a seismic event in the town of Acireale (Sicily, Italy). Given a set of sites (84 relevant buildings), and the range of the UAV, we are able to find the number of vehicles to employ and the shortest survey path. The problem of finding the shortest survey path is an operational research problem called Vehicle Routing Problem (VRP) whose solution is known to 10 be computationally time-consuming. We used the Simulated Annealing (SA) heuristic that is able to provide stable solutions in relatively short computing time. We also examined the distribution of the cost of the solutions varying the depot on a regular grid in order to assess the best area where to execute the survey.


Introduction
In the very last decade some strong earthquakes struck central Italy (L'Aquila, 2009;Emilia, 2012;and Amatrice, 2016), causing relevant social and economical consequences (Moretti et al., 2016).In the immediate aftermath of the event it is crucial to perform the fastest possible recognition of the damaged area in order to rescue as many people as possible and to assess and map the damage scenario.
At the same time, in the last decade, the unmanned aerial vehicle (UAV) greatly diffused because of its versatility (D'Alessandro et al., 2015;Gomez and Purdie, 2016): UAVs are able to self-navigate to complete a programmed mission or can be remotely controlled and prompted to perform tasks, because they have no obstacles, unlike other ground vehicle.For all these features UAVs are becoming the safest, fastest and most reliable means of carrying out observations, surveys and mapping in various types of emergency implicating large inaccessible areas (Griffin, 2004;Chou et al., 2010;Obanawa et al., 2014;Giordan et al., 2015Giordan et al., , 2017;;D'Alessandro, 2016a;Dominici et al., 2017;Jurecka and Niedzielski, 2016;Silvagni et al., 2016).Readiness and reliability are key aspects in the post-earthquake scenario, when the first problem is to detect the hotspots: the areas where threatened people needing assistance or rescue are supposed to be and the areas where the key infrastructures (e.g.hospitals, schools) are.UAV recognitions are also essential to finding accessible and safe routes for the rescue vehicles.
The tasks that the UAV (or a group of UAVs) has to accomplish should be optimized to avoid wasting time and resources.Moreover, in emergency scenarios, some questions arise: how many UAVs should be employed?What will the autonomy be?How long will the survey take?Where is the best take-off place?What is the best route for every UAV in order to minimize the travel distance and the surveying time?To solve all these questions we will apply the concepts of the vehicle routing problem (VRP) to the UAV in a possible real scenario (Agatz et al., 2016;Golden et al., 2008;Wang et al., 2017;Yu et al., 2017).In the case study we do not have to consider the capacity of transport, since the UAVs will not carry any goods; they will just have a survey task and have to visit each site (given its geographical coordinates) as quickly as possible, taking into account their own autonomy, and return back to the depot: this VRP approach is "distance-constrained" (DVRP), which means that the main task is to minimize the total length of the routes.
Our case study is the town of Acireale (Sicily, Italy) where more than 52 000 people live.The city is located along the eastern coast of Sicily, just on the flank of the Mt Etna volcano.During historic times the town of Acireale experienced severe damages, deaths and even almost complete destruction after several strong earthquakes (Azzaro et al., 2000;Rovida et al., 2016), and according to the Italian seismic hazard map it is among the areas where the highest peak ground accelerations (up to 0.25 g) are expected (Ordinanza PCM, 2006).For these reasons, a pilot urban dense seismic network realized with Micro Electro-Mechanical Systems (MEMS) sensors is being installed in the town for early warning purposes and rapid damage assessment after strong seismic events (D'Alessandro, 2014(D'Alessandro, , 2016b;;D'Alessandro et al., 2014).MEMS accelerometers are deployed in the most vulnerable buildings such as hospitals, schools and all the facilities devoted to public security (D'Alessandro and D'Anna, 2003).In this paper we will test a routing algorithm in order to find the best solutions for a complete recognition of the 84 selected sites within the town of Acireale.

Methods and materials
The specific problem we want to solve belongs to the category of operational research problems known as vehicle routing problems (VRPs).The VRP is a combinatorial optimization and integers programming problem that generalizes the travelling salesman problem (TSP).From a computational point of view, the solution to this problem is informally identified as NP-hard (Garey and Johnson, 1990), i.e. an exhaustive solution to the problem that can be performed in polynomial time is not known.Due to such computational limitations, we have decided to provide a solution by an incomplete method for optimization called simulated annealing (SA).This method has been chosen after being tested against the genetic algorithm (GA).GA is another metaheuristic but inspired by the process of natural selection and is commonly used to compute high-quality solutions to optimization and search problems (Michalewicz, 1996).The implementation of the GA we have used adopts the same representation (permutations) and cost function of the SA, a binary tournament selection process, a mutation which simply switches two elements in the solution and a crossover operator that is usually used in the GA implementation of routing algorithms, such as the traveller salesman problem.The GA was tested using three different crossover methods (Davis, 1985;Goldberg and Linge, 1985;Oliver et al., 1987), each with a different probability and a different mutation probability.A comparison between the two methods, SA and GA, showed that the overall performances of the GA were clearly inferior.The average computational time for the GA is more than double with respect the time for the SA, and the costs of the solu-tions are systematically higher.SA is faster at exploring the solutions space, because the GA suffers from always carrying a more larger population of solutions.Therefore we chose only the SA that is described here in detail, whereas GA is not discussed here.

The simulated annealing for optimization problem solution
The simulated annealing (SA) method is based on the analogy with the so-called annealing process i.e. the process consisting of heating up and subsequently cooling down a solid that consequently freezes into a minimum energy structure.Kirkpatrick et al. (1983) introduced the method as an optimization technique for combinatorial problems, but it takes inspiration from later works by Metropolis et al. (1953) and Pincus (1970).The annealing process starts providing high temperatures to the solid.The effect is that the atoms of the solid assume high-energy states so that they are more able to arrange themselves.This atomic energy reduces while the temperature is reduced until a minimum of energy state is reached where the solid reaches the crystal structure.During the whole cooling process, it is really important to apply a slow cooling: a very quick one does not involve a minimum energy state of the solid, and some kind of irregularity and defects appear suddenly in the crystal structure.
In the case of thermal equilibrium at a temperature T , the Boltzmann distribution gives the probability P (T ,s) that the system is in a given configuration: where E(s) is the energy of the configuration, k is the Boltzmann's constant and S indicates the set of all possible configurations.A system of particles in thermal equilibrium at a given temperature T can be simulated using a technique developed by Metropolis et al. (1953).Suppose at time t, the system is in configuration q, and a new candidate configuration r is generated randomly at time t+ 1.The configuration r is accepted or rejected according to the ratio p between the probabilities of being in r and in q: In particular, if p > 1 (E(r) < E(q)), then r is accepted as the new configuration at time t+ 1.If p < = 1 (E(r) > E(q)), r is accepted as the new configuration with probability p.In this process configurations with higher energy can be achieved.Moreover, it has been shown that as time t → ∞, the probability that the system is in a given configuration s is P (T ,s), regardless of the starting configuration.
When dealing with a generic global optimization problem, the analogy is done by considering the states of the solid S as the feasible space and the energy of the states E(s) as the objective optimization function.The solution of the problem is the minimum energy state.The SA algorithm is iterative and it uses the relaxation technique by Metropolis described above as a strategy for finding the solution.The implementation of the SA algorithm for the solution of a specific optimization problem, involves the following fundamental choices: i. Representation of the solution space: by mapping S → R of the feasible space S of the optimization problem onto a new space R, a new representation of the solution is provided.
ii. Cost function: the cost function C measures the quality of a given solution.
iii.Transition mechanism: to move from one state to the next, a transition mechanism P that slightly modifies the current solution is needed.
iv. Cooling schedule: the temperature at initial state, the temperature-updating rule and the number of iterations to be carried out at each step of the cooling are fundamental for the cooling process of simulated annealing.
A stopping criterion for the termination of the search process is needed.
SA's major advantage over other incomplete methods is the ability to avoid becoming trapped in local minima.The algorithm adopts a particular random search which accepts changes that can locally improve or worsen the optimization function.It can be shown that, for any given finite problem, the probability that the simulated annealing algorithm terminates with the global optimal solution approaches 1 as the annealing schedule is extended.This theoretical result is, however, not particularly helpful, since the annealing time required to ensure a significant probability of success will usually exceed the time required for a complete search of the solution space.
In order to implement the SA solution for the VRP, we represent the solutions by permutations of I + J − 2 elements where I indicates the number of sites and J the number of UAVs.This representation is such that all the indices k of s verifying s (k) > I − 1 delimitate L (j ) j = 1, .., J subset of indices in s, which represent the routes of the j -th UAV.Taking into consideration that each UAV must take off and land from the same depot d, the final route of the j -th UAV will be {d, L (j ) , d}.This representation takes inspiration from the SA matlab implementation of VRP provided at http://www.yarpiz.comand freely available to download.The starting solution is set as a random permutation, and at each iteration of the algorithm the maximum distance MD run by the UAVs and their total distance TD are computed.Note that it can occur that the algorithm sets a solution that assigns a route to an UAV whose total distance is greater than its autonomy.To avoid this situation, the binary vector of the feasible UAVs is computed so that F (j ) = 0 means that the drone j could not complete its route.Finally the cost of a solution is defined: where µ (F ) is the mean of the values in F .The goal of the SA algorithm is to minimize the function C. The adopted transition mechanisms consist of adopting one of these operations, chosen at random: i. Swap: swap two elements in s at random positions i and j .
ii. Reversion: reverse the entire content of the vector s bounded by two random indices i and j .
iii.Move: move a random elements(i) into a random position j The transition mechanism is applied m times in order to find a set N (s) that contains m neighbours of a solution s.Among them, the best one in terms of minimum cost is selected as the new solution.The initial temperature has been set to 100.At each iteration the temperature is updated following the rule T new = 0.98•T .The maximum number of iterations have been set to 1000.All these parameters, together with the integer values in the cost function C, have been chosen after a trial and error phase.

The mission and the vehicle
The main objective of the mission is to survey the potentially damaged area as fast as possible immediately after the earthquake in order to verify the state of damage of the buildings.The adopted method is able to find the best surveying solutions for a set of sites (given their geographical coordinates), identify the best depot for the UAVs and minimize the number of vehicles.The surveying sites are the nodes of the accelerometric network installed in the town of Acireale (Sicily, Italy) and some other noteworthy buildings (D'Alessandro, 2016b), for a total of 84 sites (Fig. 1).
The damage assessment at every site is assessed at first by means of the accelerometric sensors installed at each site (D'Alessandro et al., 2017) and is then precisely verified by the UAV's survey.The pictures taken during the recognition will allow a rapid image comparison between the archived images data set and the observed scenario.
In this case study some assumptions are necessary.First, we do not intend to perform photogrammetric images for 3-D reconstruction and damage assessment of the buildings.The execution would be obviously dependent on the size of the buildings themselves and is a time consuming (tens-of-hours to days) operation.Second, the surveying time over each target site can be considered equal despite the different areal extension of the sites.At notable flying height (i.e.> 50 m, the minimum ground distance for a safe survey in an urban   area), the overflight time above a site could even be neglected considering the optical characteristics of the UAV camera.
The selected UAV is a commercial and low-cost model called Phoenix.It is able to communicate with an operative centre located up to 20 km away, accomplish preprogrammed missions by bypassing eventual obstacles and executes realtime variations of a route ordered from the operative centre.Phoenix is equipped with a gimbal camera, visual positioning system, micro-USB port, micro-SD slot, obstacle detection system, LEDs, four engines, a telemetry system and other systems able to provide information about the base-vehicle connection and the battery charge status.The UAV is able to reach speeds of 20 m s −1 and the battery ensures an energetic autonomy of about 30 min.Therefore, in ideal conditions it could travel as far as 36 km.

Results
For the selected UAVs the maximum considered travel distance is 30 km.However, when taking into account any possible unfavourable atmospheric phenomena, the relevant difference in elevation between some sites and a margin of safety, we repeated the tests with progressively increasing steps of 3 km (3, 6, 9, 12, 15, 18, 21, 24, 27 and 30 km).
No solution was found for an autonomy of 3 km.For an autonomy of 6 km the solution is found with the use of three UAVs whereas, for all the other autonomy values, the solution considers the enrollment of two UAVs.All the solutions have comparable travel distances but the estimated best depot varies.Among the performed solutions, we show the ones relative to the shortest and longest autonomies: 6 and 30 km respectively (Fig. 2).The tested method showed a good velocity of convergence reaching a stable solution within about 30 s for all the selected autonomies.Table 1 summarizes, for each UAV range autonomy, the number of UAV, the cost and computing time of the solution, the distance travelled by the fleet and the average distance travelled by every single UAV.The values of the cost, computing time and distances are almost constant and do not show any relevant dependence with a different range of autonomy, except for the autonomy of 6 km.In this case, because a third UAV is considered, the values of all the parameters are sensibly out of the range of the other considered autonomies.
Immediately after an earthquake some areas will not be accessible from the ground; therefore we should also take into account a possible set of depots for the UAVs.For this reason, after the first phase, in which we calculated the best depot, we executed the algorithm from all the nodes of a regular grid.In this second phase we thus evaluated the behaviour of the algorithm given the autonomy, varying the depot within a 10 × 10 grid (∼ 8 km 2 ) of depot superimposed onto the survey area (Fig. 1).The grid size has been chosen to have a density comparable to the survey sites.We obtained the distribution of the cost for the solution executed at every grid point for the different considered autonomy.Furthermore, the re-sults are interpolated and the maps for the shortest and the longest autonomy (i.e. 6 and 30 km) are shown in Fig. 3, where the lighter shades indicate areas where the cost is lower (i.e.solutions are better) than areas with darker shades.The blue circles indicate the depot where the algorithm did not provide any solution, while among all the feasible solutions represented by the red circles, the green circle indicates the best solution (lowest cost).Solutions for all 100 nodes of the depot grid were found with autonomies greater than 9 km and the cost pattern is somehow comparable, since the solutions with the lowest costs are located in the inner part of the grid.However, for some values of autonomy this pattern is not evident and a patchy map of the cost is obtained.

Discussion and conclusions
The SA algorithm was able to find solutions with only two UAVs (i.e.best solutions), even though solutions up to 10 UAVs were explored.Considering a more populous UAV fleet, the computing time needed to get a stable solution steadily increases.Conversely, when increasing the number of UAVs, the cost of the solution sharply decreases at the beginning, and later is almost constant (    The algorithm is able to provide the best solution in terms of number of UAVs and selection of the depot given the autonomy range of the vehicle but we also took into account the eventual impossibility of operating the survey starting from the best depot in the post-earthquake scenario.We thus examined the spatial distribution of the cost for the solution executed from a grid of depots (Fig. 3).Moreover, we computed an interpolation of the averaged cost values of all the solutions at the different autonomy values and for each grid node: the pattern of the distribution is concentric and the best solutions are around the central part of the study area (Fig. 4).In a post-earthquake scenario, the depots from the inner part of the town will be considered first, if exploitable, and later the depots from the outskirt will be considered.Even though they have greater costs, solutions from the outer depots are valid and the UAVs taking off from these depots will be able to accomplish their tasks.
The SA proved to be a good method for the routing problem in terms of computing time and reliability of the solutions.Our results also indicate that the selected UAV type is suitable to perform the survey of the area and, finally, results suggest that, under equal conditions (i.e.same UAV type an number) it would be possible to carry out a survey over a wider area.Therefore, the same procedure could also be operational during the planned expansion for the accelerometric network in the town of Acireale.In the same way it could be implemented in other areas in analogous post-earthquake scenarios or in general post-emergency scenarios.

Figure 1 .
Figure 1.Location map of the town of Acireale, Italy (a) and of the 84 survey sites (b).

Figure 2 .
Figure 2. Examples of SA solutions for 6 km autonomy (a) and 30 km autonomy (b).

Figure 3 .
Figure 3. On the left is the map of depots and sites, and on the right are the isosurfaces showing the cost of the solutions for all the depots on the grid.The lighter the colour of an isosurface, the better the solutions that start from the enclosed depots.Examples of 6 km autonomy (a) and 30 km autonomy (b).

Figure 4 .
Figure4.Interpolated average value of the cost of the solution for all the considered autonomies.The lighter the colour of isosurface, the better the solutions that start from the enclosed depots.

Table 1 .
UAV number, cost, computing time and distances for the best SA solutions.

Table 2 .
Cost-solution varying autonomy and UAV number.