Complete Coverage Path Planning for a Multi-UAV Response System in Post-Earthquake Assessment

: This paper presents a post-earthquake response system for a rapid damage assessment. In this system, multiple Unmanned Aerial Vehicles (UAVs) are deployed to collect the images from the earthquake site and create a response map for extracting useful information. It is an extension of well-known coverage path problem (CPP) that is based on the grid pattern map decomposition. In addition to some linear strengthening techniques, two mathematic formulations, 4-index and 5-index models, are proposed in the approach and coded in GAMS (Cplex solver). They are tested on a number of problems and the results show that the 5-index model outperforms the 4-index model. Moreover, the proposed system could be signiﬁcantly improved by the solver-generated cuts, additional constraints, and the variable branching priority extensions.


Introduction
Rapid damage assessment in post-earthquake situations plays an important role in the early response phase activities (i.e., evacuation of injured individuals, debris collection, and relief distribution).The ground-based post-earthquake inspection is extremely time-consuming, and unhelpful in severely damaged areas; therefore, at present, aerial systems are widely used for investigations.It is highly desirable for densely populated urban areas to have a pre-planned immediate and automated post-disaster mapping and monitoring system.In the past decade, scientists have made an attempt to improve high-resolution satellite imaging and laser scanning systems to evaluate the disaster damage and loss [1][2][3][4][5][6][7][8].However, the satellite systems have many limitations for an efficient post-disaster imaging, such as weather conditions (cloudy or dust whirls), time constraints for acquiring images and uplinking the acquisition plan to the satellite, delay in satellite data delivery after collection, etc. [9,10]; therefore, much attention has been focused on utilizing small unmanned aerial vehicles (UAVs) for post-disaster mapping [11][12][13].
In recent times, UAVs have become a popular entertaining tool globally.UAVs of various sizes are available on the market that can carry film or photographic cameras for different applications.The classification of UAVs according to their classes, distinct applications, and characteristics is beyond the scope of this article, though interested readers can refer to earlier studies [13,14].Moreover, Colomina and Molina [15] have reviewed significant applications of UAVs, such as imaging and remote sensing.This study focuses on the mapping applications of small UAVs in post-earthquake situations.In a UAV-mapping system, consecutive overlapping aerial images taken by the UAVs are processed to obtain a complete map or to extract useful information.In earthquake situations, the information, such as building and bridge destructions, road blockages, and population relocation, aids the managers in organizing a more effective post-earthquake response system.Within the first 30 min, the post-earthquake survival rate is 91%, decreasing to 36.7% by the second day [11].Hence, in addition to accuracy, time is also an important factor in post-earthquake mapping.Based on this fact, it can be said that an optimal solution to the UAV coverage path-planning problem with a time-dependent objective function will have a great effect on the efficiency of the response system.
Coverage path problem (CPP), can be considered as a variant of Vehicle Routing Problem (VRP).In VRP, a fleet of vehicles start and end their tour at a single depot while visiting all nodes on the route.Unlike the classical VRP that considers a single main depot, in multi depot vehicle routing problem (MDVRP), there are several depots and the customers can be served from any of the depots.MDVRP has widely been studied in the literature and the interested readers may refer to the extensive survey of Montoya-Torres et al. [16].Based on the literature the UAV routing problem is also formulated as a VRP model where additional constraints need to be added to reflect the characteristics of the problem.UAV routing problem has attracted much attention in the past decade for different applications [17,18].CPP, tries to find an appropriate path for robots while covering the pre-defined nodes.The problem has a wide range of applications in automated harvesting, vacuum cleaning, mapping, demining, monitoring, etc.It is worth mentioning that CPP is different from the Covering Tour Problem (CTP).In CTP there are two types of nodes, those that must be visited by the tour and those that must be covered i.e., it lies within a pre-defined distance from a vertex of the tour [19].Choset [20] and Galceran and Carreras [21] conducted comprehensive surveys on CPP methods, algorithms, and recent advances.
Immediately after an earthquake, the damage and location of afflicted people can be identified by imaging and then processing the obtained images.Since time factor is very important in post-earthquake response phase, UAVs are considered to be reliable for this task.Earthquake is not the only disaster that can be monitored by UAVs, but it is the most suitable one for imaging due to absence of information about the condition of each construction over the entire area.The imaging application of UAV based on CPP necessitates further considerations.One important issue is the direction of the path while passing through a node.Another important issue that arises due to earthquake situation is that, unlike the routing problem, the maximum mission time of vehicles must be minimized.UAVs start their route from a base and end their route at the same base for refueling or at the mission completion.This study extends the CPP for an application in the post-earthquake rapid damage assessment.The CPP features on a grid-based map are adopted on a multi-depot multi-tour vehicle routing problem for an emergency situation.
Zelinsky et al. [12] were the first to investigate the CPP on a robot, starting from an origin and ending at a goal point to minimize the length, energy consumption, and travel time.Later Carvalho et al. [22] proposed an algorithm for a mobile robot in an industrial environment such that the obstacles were not pre-specified on the map.The subject of path planning has attracted a lot of attention due to the recent development of UAV systems.Li et al. [23] studied an exact cellular decomposition method for UAV path planning in a polygon region.For the purpose of precision agriculture mapping, Barrientos et al. [24] performed an experiment using an integrated tool.Initially, they subdivided the polygon area and then conducted path planning for each subarea using a multi-UAV system.Torres et al. [25] presented a path planning algorithm for a single UAV with the aim of reducing battery usage and minimizing the number of turns, coping with both convex and non-convex regions.In a study by Wang et al. [26], an algorithm was used that minimized the consumed energy by the UAV for covering 3D terrain.In addition, a distinct study [27] introduced the multi-robot boundary coverage problem with the application of inspection of blade surfaces inside a turbine.Furthermore, Galceran and Carreras [21] conducted a comprehensive survey on CPP.
Most CPP algorithms are based on boustrophedon path strategy.In this strategy, the back and forth motion in a sweep direction of a polygon tries to cover an area with a minimum number of turns.Another approach divides the decomposing area into sub-areas and then finds the visiting sequence of the sub-areas in a sweep direction base [28], describing how to find the optimal sweep direction for robots in a polygon.Huang and York claimed that the minimization of the number of turns leads to the most efficient solution (Figure 1).Li et al. [23] showed that the path with a lower number of turning motions is a more efficient coverage path for UAVs with regards to energy, distance, and travel time.
Robotics 2016, 5, 26 3 of 15 direction for robots in a polygon.Huang and York claimed that the minimization of the number of turns leads to the most efficient solution (Figure 1).Li et al. [23] showed that the path with a lower number of turning motions is a more efficient coverage path for UAVs with regards to energy, distance, and travel time.Furthermore, the map decomposition and sweep direction method were applied for a multi-UAV system [29].Due to lack of time in a post-earthquake situation, using a multi-UAV system is potentially more beneficial and practical.Avellar et al. [30] studied a multi-UAV CPP in which the nodes were pre-defined corresponding to the boustrophedon motion.Similar to earlier studies, they decomposed the region into several sub-regions and used a VRP model to find the optimal path.However, the travel time between the sub-areas and the base to a region, as well as the possibility of turning back from the middle of a region to the base for refueling, were neglected.
Although numerous studies have been conducted on the sweep direction method, path covering based on grid-based decomposition has not yet received much attention.In a study by Wang and Li [31], the desired region was divided into grids and the amount of information collected by UAV and the path-length were maximized.According to their algorithm, passing through the center of a grid is equivalent to covering the grid; however, to have a complete image of the grid, the covering is considered as entering a grid from one side and exiting from the opposite side.Therefore, the turning motions in the grids are not considered as a part of the covered area (Figure 2).The grid size can vary depending on the altitude of the UAVs and focal length of the camera; however, by increasing the number of grids the complexity of the problem will increase.In this study, it is assumed that the number of available UAVs, potential location of UAV bases, and the possible number of open bases are based on the restrictions prescribed by the decision maker.In addition, the total mission time limit should be specified earlier according to the endurance of UAVs.We have tried to address the following questions in the proposed post-earthquake mapping problem:  Furthermore, the map decomposition and sweep direction method were applied for a multi-UAV system [29].Due to lack of time in a post-earthquake situation, using a multi-UAV system is potentially more beneficial and practical.Avellar et al. [30] studied a multi-UAV CPP in which the nodes were pre-defined corresponding to the boustrophedon motion.Similar to earlier studies, they decomposed the region into several sub-regions and used a VRP model to find the optimal path.However, the travel time between the sub-areas and the base to a region, as well as the possibility of turning back from the middle of a region to the base for refueling, were neglected.
Although numerous studies have been conducted on the sweep direction method, path covering based on grid-based decomposition has not yet received much attention.In a study by Wang and Li [31], the desired region was divided into grids and the amount of information collected by UAV and the path-length were maximized.According to their algorithm, passing through the center of a grid is equivalent to covering the grid; however, to have a complete image of the grid, the covering is considered as entering a grid from one side and exiting from the opposite side.Therefore, the turning motions in the grids are not considered as a part of the covered area (Figure 2).The grid size can vary depending on the altitude of the UAVs and focal length of the camera; however, by increasing the number of grids the complexity of the problem will increase.
Robotics 2016, 5, 26 3 of 15 direction for robots in a polygon.Huang and York claimed that the minimization of the number of turns leads to the most efficient solution (Figure 1).Li et al. [23] showed that the path with a lower number of turning motions is a more efficient coverage path for UAVs with regards to energy, distance, and travel time.Furthermore, the map decomposition and sweep direction method were applied for a multi-UAV system [29].Due to lack of time in a post-earthquake situation, using a multi-UAV system is potentially more beneficial and practical.Avellar et al. [30] studied a multi-UAV CPP in which the nodes were pre-defined corresponding to the boustrophedon motion.Similar to earlier studies, they decomposed the region into several sub-regions and used a VRP model to find the optimal path.However, the travel time between the sub-areas and the base to a region, as well as the possibility of turning back from the middle of a region to the base for refueling, were neglected.
Although numerous studies have been conducted on the sweep direction method, path covering based on grid-based decomposition has not yet received much attention.In a study by Wang and Li [31], the desired region was divided into grids and the amount of information collected by UAV and the path-length were maximized.According to their algorithm, passing through the center of a grid is equivalent to covering the grid; however, to have a complete image of the grid, the covering is considered as entering a grid from one side and exiting from the opposite side.Therefore, the turning motions in the grids are not considered as a part of the covered area (Figure 2).The grid size can vary depending on the altitude of the UAVs and focal length of the camera; however, by increasing the number of grids the complexity of the problem will increase.In this study, it is assumed that the number of available UAVs, potential location of UAV bases, and the possible number of open bases are based on the restrictions prescribed by the decision maker.In addition, the total mission time limit should be specified earlier according to the endurance of UAVs.We have tried to address the following questions in the proposed post-earthquake mapping problem: ( 1  In this study, it is assumed that the number of available UAVs, potential location of UAV bases, and the possible number of open bases are based on the restrictions prescribed by the decision maker.In addition, the total mission time limit should be specified earlier according to the endurance of UAVs.We have tried to address the following questions in the proposed post-earthquake mapping problem: (1) which potential UAV base must be open; (2) how many UAVs among all the available ones should be assigned to each open base; (3) what must be the passing sequence of squares in each route; and (4) what is the minimum time required for covering all residential squares with the existing limitations.The proposed problem has a Minimax objective function that minimizes the maximum traveling time of each UAV in order to find the fastest possible mapping schedule.Each UAV belongs to a UAV base from where its mission starts and ends, and battery replacement or refueling is performed at the same base.Each UAV's mission time cannot exceed a maximum flight time (MT) due to its endurance.To the best of our knowledge, this problem has not been studied thus far, and the formulations presented here are the first mixed integer linear programming (MILP) mathematical models of the grid-based CPP problem.We call this problem multi-UAV multi-tour location coverage path planning (MUTLCPP).The application of this model is a post-earthquake mapping system due to the Minimax objective function, and a grid-based decomposition of the residential areas.
The remainder of this study is organized as follows: Section 2 represents the problem definition and two 4-index and 5-index mathematical models.Section 3 contains the extensions and techniques used to improve the performance of the branch and cut algorithm.Section 4 includes the test problems and computational results, and Section 5 contains conclusions and future prospects of CPP.

Problem Definition and Mathematical Formulations
Associated with the MUTLCPP problem, the sets, parameters, and decision variables are defined as follows (in Table 1): The map of earthquake struck region was decomposed into equal sized squares, according to the UAV's camera footprint.The camera footprint (CF) depends on the camera sensor size (CSS), lens focal length (FL), and flight altitude (FA).Figure 3 shows the camera footprint, which can be calculated as follows.The camera plays a very important role in CPPs especially for a post-earthquake situations.Despite the resolution demanded by the decision maker, the focal length of the lens and the sensor size of the camera have a great impact on the mapping schedule length.A smaller focal length with a larger sensor size ends up in a larger camera footprint (Equation ( 1)) and consequently a lower number of squares in the map and schedule length.Also, the horizontal flight disposition is preferred for covering residential areas by a sequence of overlapping images, although the new UAV systems allow changes in the viewing angle of the installed camera.
The ground sampling distance (GSD) is a scale of resolution in aerial photogrammetry.It is the distance between two successive pixel centers on the ground.The smaller value of the GSD leads to the higher resolution of the image and more differentiable details.At a fixed focal length, decreasing the altitude results in a smaller GSD, and at a fixed flight altitude, by increasing the focal length, the GSD value will decrease.In an earthquake situation, the GSDs of several meters may suffice for damage assessment of non-urban regions.However, it is quiet insufficient when analyzing highly populated earthquake-struck cities and detailed images of at least 1 m GSD are essential [32].
While the geographical mapping might not require a quick scheduling process, in post-earthquake mapping, time factor has a great importance.Since the use of UAVs for mapping is in its early days, there are potential improvements and considerable innovations related to unseen problems and technical challenges.
To ease the image processing, an overlap from each side is essential; thus, if the overlap is α% from each side, then the square edge must be equal to (1 − 2α)% of the camera footprint length.The overlap is required to discover the common details among the images, however; the minimum and maximum overlaps depend on the post-image-processing software.Considering the fact that imageprocessing discussion is beyond the scope of this study, and there is no universally accepted overlap standard [33], we only mention that the overall overlap for a specified-path mapping must be at least 20% of the image width [34].
The MUTLCPP is defined on an undirected graph (, ), where  =   ∪   is a set of squares and  = {(, )|,  ∈ } is the set of arcs.  is a set of potential depot squares and   =   ∪   , where   is the set of residential squares that must be covered and   is the set of empty squares.Non-negative travel time,   , and travel distance,   , are associated with each arc with a triangular inequality (i.e.,   +   ≥   ).Moreover, the set of homogenous UAVs is denoted by  = {1,2,3, . ., } and the set of routes by  = {1,2,3, … , }.There is a set  = {, ℎ, , , } that specifies the leaving direction form one square to another.All the squares belonging to   must be covered, at least once, such that the UAV enters and exists from opposite sides (i.e., Right-Left, Left-Right, Up-Down, and Down-Up).Among the |  | number of potential depots,  number of them can be opened, according to the restrictions prescribed by the decision maker.The flight time The camera plays a very important role in CPPs especially for a post-earthquake situations.Despite the resolution demanded by the decision maker, the focal length of the lens and the sensor size of the camera have a great impact on the mapping schedule length.A smaller focal length with a larger sensor size ends up in a larger camera footprint (Equation ( 1)) and consequently a lower number of squares in the map and schedule length.Also, the horizontal flight disposition is preferred for covering residential areas by a sequence of overlapping images, although the new UAV systems allow changes in the viewing angle of the installed camera.
The ground sampling distance (GSD) is a scale of resolution in aerial photogrammetry.It is the distance between two successive pixel centers on the ground.The smaller value of the GSD leads to the higher resolution of the image and more differentiable details.At a fixed focal length, decreasing the altitude results in a smaller GSD, and at a fixed flight altitude, by increasing the focal length, the GSD value will decrease.In an earthquake situation, the GSDs of several meters may suffice for damage assessment of non-urban regions.However, it is quiet insufficient when analyzing highly populated earthquake-struck cities and detailed images of at least 1 m GSD are essential [32].
While the geographical mapping might not require a quick scheduling process, in post-earthquake mapping, time factor has a great importance.Since the use of UAVs for mapping is in its early days, there are potential improvements and considerable innovations related to unseen problems and technical challenges.
To ease the image processing, an overlap from each side is essential; thus, if the overlap is α% from each side, then the square edge must be equal to (1 − 2α) % of the camera footprint length.The overlap is required to discover the common details among the images, however; the minimum and maximum overlaps depend on the post-image-processing software.Considering the fact that image-processing discussion is beyond the scope of this study, and there is no universally accepted overlap standard [33], we only mention that the overall overlap for a specified-path mapping must be at least 20% of the image width [34].
The MUTLCPP is defined on an undirected graph G (S, A), where S = S O ∪ S D is a set of squares and A = {(i, j)|i, j ∈ S} is the set of arcs.S D is a set of potential depot squares and S O = S C ∪ S E , where S C is the set of residential squares that must be covered and S E is the set of empty squares.Non-negative travel time, t ij , and travel distance, d ij , are associated with each arc with a triangular inequality (i.e., t ik + t kj ≥ t ij ).Moreover, the set of homogenous UAVs is denoted by U = {1, 2, 3, .., u} and the set of routes by L = {1, 2, 3, . . . ,l}.There is a set R = {Le f t, Right, U p, Down, Direct} that specifies the leaving direction form one square to another.All the squares belonging to S C must be covered, at least once, such that the UAV enters and exists from opposite sides (i.e., Right-Left, Left-Right, Up-Down, and Down-Up).Among the |S D | number of potential depots, TD number of them can be opened, according to the restrictions prescribed by the decision maker.The flight time limitation for each UAV is specified by MT.All UAVs start and end their tours at the same base.
Because time plays a vital role in the post-disaster response phase against cost, the proposed model must minimize the maximum mission time of all UAVs.
When M is a sufficiently large positive number, the 4-index and 5-index models of MUTLCPP are presented as follows:

Model:
Min T The objective Equation ( 1) is to minimize the maximum mission time of UAVs.This objective function leads to the minimum mapping schedule length for all UAVs.The MUCLTPP's constraints are as follows: Equation ( 3) ensures that the total number of open bases is TD, and Equation ( 4) imposes that each UAV has to be allocated to one depot only.
Equation ( 5) ensures that each square is passed at most once in each route.
Equation ( 6) guarantees that if a UAV enters a square then it must leave that square.In other words, the total number of entrances to a square is equal to the total number of exits from that square.
Equation (7) assures that the UAVs exit from an open base only.This constraint forces the w ul ij variable to take the value of zero if the base is close (y i = 0).
Equation ( 8) assigns zero to v ul ij variables according to the movement variable of Equation ( 9) indicates that if a UAV pass through square i ∈ S O , then the remaining number of the nodes in route l will decrease by one.
Equation (10) ensures that the duration of each route does not exceed the UAV's flight time limitation MT.
The value of T which is the maximum mission time of UAVs is determined by Equation (11).
Equation (12) indicates that UAVs can only leave their own depot and Equation ( 13) eliminates the illegal routes.
Equation ( 14) prevents the formation of arcs with the length more than dm for the horizontal and vertical directions.dm is the distance between the center of two consecutive squares; the length of horizontal or vertical arcs cannot exceed it.
Equation ( 15) specifies that the value of g lr i can be one if the entrance and exit are in the same horizontal or vertical directions. ∑ Equation ( 16) ensures that all residential squares must be covered by a UAV in either horizontal or vertical directions and Equation (17) defines the value of variable x lr ij .
x l U p ij = 0 ∀i, j ∈ S i − 1 ≤ j, l ∈ L (18) Equations ( 18)-( 21) imposes the value of zero for the wrong horizontal and vertical directions.For example Equation (18) ensures that the travel in the upward (up) direction from i to j is zero unless i − 1 is less than or equal j.The travel from one square to its right square must be impossible for the up, down, and left directions and it is shown in Equation (21).

5-Index Formulation
The 5-index model development is presented in the same structure of the 4-index formulation as follows: Min T (1) S.t (2) , (3) , and Equations ( 24)-( 41) are the same equations as Equations ( 5)-( 16) and Equations ( 18)-( 23) in the 4-index formulation.In this problem, the size of |L| is problem dependent, and an increase in the number of routes results in an exponential increase in the execution time of the problem.The minimum size of |L| equals to the size of |U|; the maximum size can be calculated according to any feasible solution that the problem solver decides.However, it is preferable to use solutions that connect the squares (horizontally or vertically) with minimum turns.

Methodology
The 4-index and 5-index models were coded in GAMS and the Cplex solver was selected.To improve the algorithm and make it more effective, some extensions and techniques were found to be appropriate for our problem.Branch-and-Bound is the main algorithm used for solving MIP problems in Cplex.Using well-known techniques based on the model's features can effectively reduce the execution time of the problem and improve the performance of MIP for effective collaborative decision-making in disaster management.The following subsections shortly describe the techniques that were effective for the models.

Solver Generated Cuts and Heuristics
In the early sixties, Gomory [35] introduced the procedure of fractional cuts.The procedure was considered impractical until a study [36] proved that Gomory's cuts can be effectively distributed over different branches of the search tree in a branch and cut framework for a mixed 0-1 problem.In the present time, most optimizers use Gomory's cuts for MILP problems.In the Cplex optimizer, the fraccuts option decides whether or not Gomory's fractional cuts should be generated for the problem; and the parameters 0, 1 and 2 control the cuts to be generated automatically, moderately, and aggressively, respectively.In the case of lack of progress in the best bound, one effective solution uses more aggressive levels of cut generation to further tighten the formulation [37].As initial observations confirm the lack of progress in the best-found-bound, parameter 2 (aggressive cut generation) was found to be more effective for the proposed models.According to the structure of the mathematical formulations, cover cuts are also useful.Therefore, Gomory's fractional cuts and cover cuts were aggressively generated for the models.
GAMS (General Algebraic Modeling System) is a high level modeling language designed for programming the linear, nonlinear, and mixed integer optimization models.GAMS can support different forms of linear and nonlinear programming problems.The system has a language compiler and an integrated high-performance commercial solver, e.g., Cplex.Cplex is a famous commercial optimizer that can solve linear, mixed integer, quadratic, and quadratically constrained programming problems.GAMS provides access to the Cplex solver with a variety of options to interact with it and improve the solution procedure.Cplex uses several alternative algorithms such as primal simplex, dual simplex, barrier, or sifting algorithms.For mixed integer programming problems, Cplex uses the branch-and-cut algorithm.
Cplex takes advantage of different heuristic techniques to avoid tree pollution, increasing the diversity of the branch-and-bound search, achieving integer feasible solution in shorter time, and fastening the final optimality.Based on the characteristics of the proposed problem, two heuristic features from the different heuristic procedures offered by Cplex were utilized, i.e., local branching heuristic (LBHeur), and node heuristic (NHeur).The LBHeur was introduced by Fischetti and Lodi [38] and has been widely used by commercial optimizers.We can infer the following facts based on the ILOG Cplex explanations:

•
LBHeur attempts to enhance new incumbents found throughout an MIP search.The LBHeur algorithm will be appealed only if a new incumbent is identified.In cases where several incumbents are identified at a single node, the last one will be considered.

•
NHeur implements the techniques to find a feasible solution from the present node in branch and bound algorithm.The parameter HeurFreq restricts this offer.A positive number specifies the frequency of invoking NHeur (in the number of nodes).

Additional Constraints
User cuts or additional constraints are the polynomial constraints that strengthen the mathematical formulation and reduce the solution space.These cuts are not supposed to affect the feasible solutions space.The following constraints eliminate unfavorable fractional solutions and expedite the execution process: Equation ( 42) imposes that the total number of squares that are to be visited by UAV u in route l right after leaving the depot is equal to the total number of squares allocated to UAV u in route l.Equation (43) ensures that the value is equal to 0 just before ending the tour at the depot.Equations ( 9), ( 42) and ( 43) are adapted from a previous study [39].
Equation ( 44) indicates that if a UAV belongs to a depot, then the depot must be open.
Equation ( 45) ensures that the total number of arcs leaving the depots is less than or equal to the number of the available routes.
According to the logic of the problem and the grid-based decomposed map, the UAVs must travel directly through the nodes of the S E set, and it is applied in Equation ( 46).As mentioned earlier, only the residential squares must be covered by horizontal or vertical movements.In other words, in favor of the objective function value, and in order to find the minimal path, the UAVs should travel in the direct direction while passing through the nonresidential squares.
Equation (47) ensures that the number of times a UAV enters a square in each route is more than or equal to the number of times the square is covered in that route.

Variable Branching Priority
The variable branching priority based on the mathematical formulation's specifications is one of the most efficient techniques that speeds up the solution procedure and tightens the model.Without using this feature, the solver automatically decides the appropriate variable to be branch on.It is imperative to select a suitable priority among the existing variables depending on the variables' dependence, relevance, and significance in the proposed mathematical model.In this study there are 4-index and 5-index mathematical formulations of a single problem.In this study, 4-index and 5-index mathematical formulations were used for a single problem.Based on the model's characteristics and variables' significance, the following priorities were suggested: I.
For the 4-index model, the variable branching priority order must be g lr i , x lr ij , w ul ij , v ul ij , p iu and y i .II.
The branching priority of the 5-index model should be g lr i , v ul ij , w lr iju , p iu , and y i .

Computational Experiments and Results
In this section, some instances are designed to test the validity of the models and features proposed in Section 3. Based on the problem's characteristics, it is assumed that the number of available UAVs, number of possible open depot squares, and maximum limit for the total traveling time of each UAV are known.Each test instance was named in the X-Y-Z-L format, where X refers to the total number of squares in the map, Y refers to the total number of residential squares, Z represents the number of potential UAV bases, and L is the number of available UAVs.The test instances are shown in Figure 4.For instance, in problem 36-11-2-3 map, where 36 is the total number of squares, 11 squares are residential areas, two squares are potential depots, and three represents the available UAVs.It is assumed that each small square side is one in terms of length unit and the speed of the UAVs is also one in terms of speed unit.The mathematical models were coded in GAMS 24.7.3 and solved using the Cplex 12.6.3.0 solver.The models were run on a laptop with the configuration of Intel ® Core™ i7-2620M CPU 2.70 GHz, RAM 8.00 GB, and windows 10 Enterprise operating system.
To evaluate the efficiency of the techniques mentioned in Section 3, each technique was considered alone in an individual run at the first step.Subsequently, the combination of efficient ones was used.All the problems were solved in a 5-h time limit, irrespective whether a feasible or optimal solution was found or not.The reason for the 5-h time limit is that no significant improvement would occur by increasing this limit.The calculation time increases exponentially by increase in the number of squares or the number of routes.Tables 2 and 3 show the results of the 4-index and 5-index models, respectively, in a combination of techniques mentioned in Section 3. It is assumed that each small square side is one in terms of length unit and the speed of the UAVs is also one in terms of speed unit.The mathematical models were coded in GAMS 24.7.3 and solved using the Cplex 12.6.3.0 solver.The models were run on a laptop with the configuration of Intel ® Core™ i7-2620M CPU 2.70 GHz, RAM 8.00 GB, and windows 10 Enterprise operating system.
To evaluate the efficiency of the techniques mentioned in Section 3, each technique was considered alone in an individual run at the first step.Subsequently, the combination of efficient ones was used.All the problems were solved in a 5-h time limit, irrespective whether a feasible or optimal solution was found or not.The reason for the 5-h time limit is that no significant improvement would occur by increasing this limit.The calculation time increases exponentially by increase in the number of squares or the number of routes.Tables 2 and 3 show the results of the 4-index and 5-index models, respectively, in a combination of techniques mentioned in Section 3.
According to the results presented in Tables 2 and 3, it can be concluded that the 5-index model shows better solutions in terms of CPU time and the absolute gap.Furthermore, as previously anticipated, the combination of the proposed LP strengthening techniques offered better solutions in a shorter amount of time.The CPU time of the 5-index problem 16-4-2-4 was 11,320 s in the original model, which decreased to 479 s by adding the techniques.In addition, for problem 36-12-4-4, a feasible solution was not obtained in the original model; however, with a combination of techniques, a solution with 6.3% relative gap was acquired.Furthermore, while no feasible solution was found for the original model of problem 39-20-3-2 in the 5-h time limit, the optimal solution was found under the influence of the techniques in 10,554 s. Figure 5   model, which decreased to 479 s by adding the techniques.In addition, for problem 36-12-4-4, a feasible solution was not obtained in the original model; however, with a combination of techniques, a solution with 6.3% relative gap was acquired.Furthermore, while no feasible solution was found for the original model of problem 39-20-3-2 in the 5-h time limit, the optimal solution was found under the influence of the techniques in 10,554 s.

Conclusions and Future Prospects
A multi-UAV and multi-tour CPP problem based on the grid map decomposition was investigated.MUTLCPP supports the application of post-earthquake CPP where high accuracy is required.Moreover, due to the speedy decrease in the survival rate in post-earthquake situations, the Minimax objective function was considered for the optimization problem.Two mathematical models of MUTLCPP were presented and coded in GAMS, and the 5-index model outperformed the 4-index formulation in terms of solution quality and execution time.The models were improved by additional constraints and other techniques described in Section 3. On the basis of the results of the test instances, it can be affirmed that the combination of the proposed techniques sped up the execution procedure and gained solutions with a lower gap in a shorter amount of time.In the small size instance of 16-7-3-2, the combination of techniques decreased the calculation time by 85% for the 5-index model but increased by 7% for the 4-index formulation.The results depicted in Tables 2 and  3 show that the linear strengthening techniques are more effective for the 5-index model.Based on the results of Tables 2 and 3, the proposed techniques decreased the CPU time by 85%, 95%, 94%, 96%, and 79% for the first five instances of the 5-index model and decreased the time by 79%, 94%, more than 91%, and 79% for the corresponding five instances of the 4-index model.For the further three instances, while the solution is not available in the original models, by the intervention of the linear strengthen techniques optimal solutions or very good feasible solutions were found.Generally, by comparing the results (Tables 2 and 3), it is clear that the 5-index model performs much better that the 4-index formulation.Another conclusion based on the obtained results is that the solver generated

Conclusions and Future Prospects
A multi-UAV and multi-tour CPP problem based on the grid map decomposition was investigated.MUTLCPP supports the application of post-earthquake CPP where high accuracy is required.Moreover, due to the speedy decrease in the survival rate in post-earthquake situations, the Minimax objective function was considered for the optimization problem.Two mathematical models of MUTLCPP were presented and coded in GAMS, and the 5-index model outperformed the 4-index formulation in terms of solution quality and execution time.The models were improved by additional constraints and other techniques described in Section 3. On the basis of the results of the test instances, it can be affirmed that the combination of the proposed techniques sped up the execution procedure and gained solutions with a lower gap in a shorter amount of time.In the small size instance of 16-7-3-2, the combination of techniques decreased the calculation time by 85% for the 5-index model but increased by 7% for the 4-index formulation.The results depicted in Tables 2 and 3 show that the linear strengthening techniques are more effective for the 5-index model.Based on the results of Tables 2 and 3, the proposed techniques decreased the CPU time by 85%, 95%, 94%, 96%, and 79% for the first five instances of the 5-index model and decreased the time by 79%, 94%, more than 91%, and 79% for the corresponding five instances of the 4-index model.For the further three instances, while the solution is not available in the original models, by the intervention of the linear strengthen techniques optimal solutions or very good feasible solutions were found.Generally, by comparing the results (Tables 2 and 3), it is clear that the 5-index model performs much better that the 4-index formulation.Another conclusion based on the obtained results is that the solver generated cuts (Section 3.1) have less effect, and additional constraints (Section 3.2) have the maximum effect on decreasing the CPU time.
There are several interesting directions for further investigations.Some new valid inequalities can be added to the model to strengthen the MILP formulations.In addition, exact methods such as the lagrangian relaxation can be developed to find the optimal or near-optimal solutions to large real size problems.The most interesting subject for further investigations is designing a heuristic or meta-heuristic algorithm to solve large-scale problems of MUTLCPP.Each randomly generated individual solution in a heuristic algorithm can be represented in such a way that each cell takes an integer <6 indicating the UAV direction (set R).Finally, one can extend the problem by reflecting the map rotation possibility to test the different covering directions.

Figure 1 .
Figure 1.The right polygon is covered by lesser waste (lesser number of turns) than the left one.(a) The vertical direction with 11 turns; (b) The horizontal direction with 4 turns.

Figure 2 .
Figure 2. The gray cells indicate the cells that have been covered by the unmanned aerial vehicle (UAV), according to entrance and exit from the opposite sides.The picture depicts a feasible solution of grid decomposition-based coverage path problem (CPP).
(1) which potential UAV base must be open; (2) how many UAVs among all the available ones should be assigned to each open base; (3) what must be the passing sequence of squares in each route; and (4) what is the minimum time required for covering all residential squares with the existing limitations.The proposed problem has a Minimax objective function that minimizes the maximum traveling time of each UAV in order to find the fastest possible mapping schedule.Each UAV belongs

Figure 1 .
Figure 1.The right polygon is covered by lesser waste (lesser number of turns) than the left one.(a) The vertical direction with 11 turns; (b) The horizontal direction with 4 turns.

Figure 1 .
Figure 1.The right polygon is covered by lesser waste (lesser number of turns) than the left one.(a) The vertical direction with 11 turns; (b) The horizontal direction with 4 turns.

Figure 2 .
Figure 2. The gray cells indicate the cells that have been covered by the unmanned aerial vehicle (UAV), according to entrance and exit from the opposite sides.The picture depicts a feasible solution of grid decomposition-based coverage path problem (CPP).
) which potential UAV base must be open; (2) how many UAVs among all the available ones should be assigned to each open base; (3) what must be the passing sequence of squares in each route; and (4) what is the minimum time required for covering all residential squares with the existing limitations.The proposed problem has a Minimax objective function that minimizes the maximum traveling time of each UAV in order to find the fastest possible mapping schedule.Each UAV belongs

Figure 2 .
Figure 2. The gray cells indicate the cells that have been covered by the unmanned aerial vehicle (UAV), according to entrance and exit from the opposite sides.The picture depicts a feasible solution of grid decomposition-based coverage path problem (CPP).
represents the number of potential UAV bases, and L is the number of available UAVs.The test instances are shown in Figure4.For instance, in problem 36-11-2-3 map, where 36 is the total number of squares, 11 squares are residential areas, two squares are potential depots, and three represents the available UAVs.

Figure 4 .
Figure 4. Test instance maps (the gray squares are assumed to be residential areas that should be covered by UAVs and squares with black triangles are potential UAV bases).

Figure 4 .
Figure 4. Test instance maps (the gray squares are assumed to be residential areas that should be covered by UAVs and squares with black triangles are potential UAV bases).

Figure 5 Figure 5 .
Figure 5.The optimal solutions of problems 39-20-3-2 and 36-11-2-3.(a) There are two available UAVs and each belongs to one base; (b) There are three available UAVs and base number 15 is selected to be open.

Figure 5 .
Figure 5.The optimal solutions of problems 39-20-3-2 and 36-11-2-3.(a) There are two available UAVs and each belongs to one base; (b) There are three available UAVs and base number 15 is selected to be open.

Table 1 .
The definition of sets, parameters, and variables.
ijuBinary decision variable equals to 1, if UAV u exists from r side of square i to square j in route l, and 0 otherwise

w ul ij
Binary variable equals to 1, if a UAV enters from r side of i and exits from the opposite side in route l glr Binary decision variable equals to 1, if UAV u travels from square i to square j in route l, and 0 otherwise

i
Binary variable equals to 1, if a UAV enters from r side of i and exits from the opposite side in route l TThe maximum mission time of UAVs

Table 2 .
Effects of techniques described in Sections 3.1-3.3on the 4-index mathematical model.

Table 3 .
Effects of techniques described in Sections 3.1-3.3on the 5-index mathematical model.