Collaborative Hybrid Aerial and Ground Vehicle Routing for Post-Disaster Assessment

Abstract

Hybrid aerial and ground vehicles are seen as a promising option for deployment in a post-disaster assessment due to the risk of infrastructure damage that may hinder the assessment operation. The efficient operation of the hybrid aerial and ground vehicle, particularly routings, remains a challenge. The present study proposed a collaborative hybrid aerial and ground vehicle to support the operation of post-disaster assessment. The study developed two models, i.e., the Two-Echelon Vehicle Routing Problem combined with Assignment (2EVRPA) and the Two-Echelon Collaborative Vehicle Routing Problem (2ECoVRP) to evaluate optimal routings for both aerial and ground vehicles. The difference lies in the second echelon in which the 2EVRPA uses a single point-to-point assignment, whereas the 2ECoVRP considers the collaborative routings between the ground vehicle and the aerial vehicle. To demonstrate its applicability, the developed models were applied to solve the post-disaster assessment for the Mount Merapi eruption in Yogyakarta, Indonesia. Sets of numerical experiments based on the empirical case were conducted. The findings indicate that the 2ECoVRP performs better than 2EVRPA in terms of the total operation time. The tabu search algorithm was found to be a promising method to solve the models due to its good quality solution and computational efficiency. The deployment of eight drones appears to be optimum for the given network configuration of the studied case. Flight altitude and battery capacity were found to be influential to the operation time, hence requiring further exploration. Other potential avenues for future research are also discussed.

1. Introduction

In humanitarian operations, a faster response is critical to save more human lives and reduce casualties [1]. Due to their quick, flexible, and cheap operation, unmanned aerial vehicles (UAVs), commonly known as drones, are seen as a promising technology to support humanitarian operations, particularly with respect to post-disaster assessment and last-mile delivery for aid. The drone has been widely implemented in diverse applications such as military, agriculture, forestry, geographical surveillance, sports, and entertainment [2,3]. The application of the drone in humanitarian operations is, however, very rare. Nevertheless, there is growing interest to deploy drones in humanitarian operations [4].

Humanitarian operations are characterized by high uncertainty, dynamicity, and the goal of alleviating human suffering/saving human life, implying that time is a critical factor [1]. Typical disasters such as earthquakes, tsunamis, volcano eruptions, and pandemics are unpredictable. They pose challenges on humanitarian operations in which affected areas, infrastructure damage (e.g., roads, bridges), and the number and location of victims are unknown. That uncertainty is increased as the conditions change over time. Although most disasters are barely predictable, several efforts have been conducted to minimize the impacts. For instance, during the preparedness phase, prepositioning inventory, dedicated shelters, and evacuation routes have been pre-arranged [5]. During emergency response (post-disaster), humanitarian operations such as search and rescue missions and aid delivery focus on providing a fast response in an efficient way (due to limited resources). To perform efficient and immediate emergency response, information on the number and location of victims, incurred damage, and logistics capability should be on-hand. Post-disaster assessment to provide accurate and timely required information is therefore necessary. The assessment should be conducted as quickly as possible: the faster the rescue operations can be performed, the more likely the victims are to be rescued. This is also the case for aid delivery. The faster the demands for aids are known, the faster the required aids can be delivered, and consequently, the more likely people will be saved [4]. Another challenge encountered by the post-disaster assessment is the widespread areas needing to be evaluated. Some areas may be inaccessible due to challenging terrain, debris blockage, damaged infrastructure, or other hazards. Resources available for the post-disaster assessment are limited, particularly in developing countries where vehicle resources are quite limited and costly and there are few means of transport. In contrast, drones that are equipped with a high-resolution camera and global positioning system (GPS) modules can provide spatial imaging in a relatively short time and at a low cost. Due to its flexibility and its speed, the drone is particularly useful for post-disaster assessment to collect critical data that are unavailable to the response team. In many cases, disaster-affected areas are inaccessible and unsafe for people. The drone, which is operated remotely, opens new possibilities for inaccessible, hazardous, or cut-off regions. McRae, et al., [6] have demonstrated that drone real-time imagery and GPS coordinates have provided rescuers with valuable information on victim status and location as well as terrain, which allows the rescuers to judiciously plan routes and allocate appropriate resources and personnel. In pandemic situations, drones are useful as transportation modes for delivering aid to reach inaccessible areas [7,8] and reducing the risk of disease spread [9,10]. Furthermore, drones are more efficient in terms of energy consumption/emissions per unit distance than conventional diesel vans [11,12], thus supporting sustainability in humanitarian operations. In light of Sopha, et al., [13], who have highlighted the importance of incorporating sustainability in humanitarian operations to enable long-term solutions, drones contribute toward a more sustainable operation in post-disaster assessment.

Despite the significant role of the drone in humanitarian operations, existing studies addressing this issue are few. Current literature related to the application of drones for humanitarian operations focuses on the technical aspects of drones, such as wireless sensor networks [14], technical evaluation of the effectiveness of the drones [15] and developing image datasets for human detection for search and rescue [4]. In addition to the drones’ technicalities, the operation of the drones such as location, allocation, and routing should also be optimized to support efficient and effective post-disaster assessment. Location addresses where the drones are launched, e.g., at a single depot or multiple depots. Owing to the limited number of drones to perform the post-disaster assessment, the allocated number of drones to each depot is important. The routing of the drones determining the optimal sequence of the area visited by the drones should also be optimized. The present paper, therefore, explores potential models for route planning of the hybrid unmanned aerial and ground vehicles during the post-disaster assessment. Two routing models, i.e., Two-Echelon Vehicle Routing Problem combined with Assignment (2EVRPA) and Two-Echelon Collaborative Vehicle Routing Problem (2ECoVRP) are proposed and contrasted. Post-disaster assessment during the Mount Merapi eruption in Yogyakarta, Indonesia, was used as a studied case to demonstrate the models’ applicability. Tabu search heuristics was developed to solve the medium- and large-scale routing problems.

Because post-disaster assessment should be conducted in widespread areas as quickly as possible, the hybrid aerial and ground vehicle approach should be deployed. Ground vehicles have more capacity and longer operation time than drones. However, ground vehicles are constrained by road network and congestion, consequently having limited accessibility to disaster-affected areas, and they also move more slowly. In contrast, drones can move faster than ground vehicles because the road network is not a constraint. Drones do not experience congestion, and they consume less energy because they are lightweight. However, the drone has a limited operation time due to its battery life. The idea of the hybrid system is basically to take advantage of both vehicles. The typical combination of ground vehicles and drones has been used in commercial last-mile delivery [16,17,18]. Hence, instead of optimizing the drone’s routings, the present study simultaneously optimizes the routings of both the ground vehicle and the drone to globally optimize the assessment operation.

Motivated by the requirement toward practical application of the hybrid ground vehicle and drone operation, the present study developed two models, 2EVRPA and 2ECoVRP, to evaluate optimal routings of both ground vehicle and drone. The first echelon addresses the routing of the ground vehicle, whereas the second echelon addresses the routing of the drone. The 2EVRPA model uses the routing problem in the first echelon and the assignment problem in the second echelon. Unlike the 2EVRP, the 2ECoVRP includes the collaborative operation between the ground vehicle and drone operations (implementing routing problems for both echelons). The 2ECoVRP allows the drone to visit multiple targets in one go in the second echelon.

Some studies using drones in humanitarian operations have existed, such as Mishra, et al., [4], Chowdhury, et al., [19], Cannioto, et al., [20], Oruc and Kara [21], and Luo, et al., [22]. For the hybrid of ground vehicle and drone for post-disaster assessment, Luo, et al., [22] developed the model and tested it through experiments using hypothetical datasets. Otto, et al., [23] indicated that the literature using the combined system of ground vehicle and drone is dominated by experimental and hypothetical studies, in which only 10% were based on the empirical case. In contrast, Asih, et al., [24] demonstrated that the metaheuristics performance applied to hypothetical data is slightly different from that using empirical data. Complementing the previous studies, the present study focuses on the quantitative evaluation of the 2EVRP and its variants. The two proposed models are then implemented in a real/empirical case of post-disaster assessment of the Mount Merapi eruption.

Hence, the present study contributes in two ways. First, the study develops and compares the two models capturing the operational characteristics of post-disaster assessment by introducing the hybrid ground vehicle and drone system. Second, it provides an empirical contribution to evaluating efficient routing strategy for the hybrid ground vehicle and drone in post-disaster assessment for the Mount Merapi eruption.

The paper is structured as follows. This section has highlighted the motivation and the contribution of the paper. Section 2 provides brief reviews of the literature on drone operations, particularly related to the routing problem. Section 3 describes the development of the mathematical models and the solution method, which is followed by the description of the empirical case, the numerical results and analysis, and the practical implications in Section 4. Section 5 concludes and discusses future research.

2. Literature Review

Literature on drones has increased rapidly in recent years and that trend seems likely to continue in the future. The first drone literature was published in 2001 for civil operations dealing with routing problems [23] and was followed by the second publication in 2005 [3]. Since then, the numbers of drone studies were stagnant until 2013, when the body of drone literature began growing dramatically. The development was triggered by Amazon, who had started to develop delivery drones in 2013 [25]. Further, Otto, et al., [23] indicated that more than 75% of the articles have been published in the last seven years. To date, four literature reviews on drone operations exist. Otto, et al., [23] surveyed the optimization approaches for drone operations in civil applications, Khoufi, et al., [26] and Viloria, et al., [3] narrowed the review focusing on the routing aspect including types of problems and resolution methods, and Kellerman, et al., [27] focused on the potential barriers, proposed solutions, and expected benefits of drones for parcel and passenger transportation.

Based on those reviews, drones have been used in diverse applications such as transportation (delivery), communication, logistics processes, surveillance, and monitoring [3]. It is worth noting that drone applications in disaster operations are still very limited, as most drone applications are for commercial delivery. The articles for commercial transportation/delivery applications account for 54% and 41% of the total surveyed articles based on Khoufi, et al., [26] and Viloria, et al., [3], respectively. In contrast, only 3% and 4% of the surveyed articles support disaster operations according to Otto, et al., [23] and Khoufi, et al., [26], respectively. The most recent survey by Viloria, et al., [3] indicated a slight increase of 6% of the articles addressing humanitarian operations, mostly addressing the delivery in humanitarian logistics. Only a few have dealt with post-disaster assessment.

Drones can support various operations such as mapping, monitoring/surveillance, communication, transporting/delivering, and data gathering. Otto, et al., [23] categorized the types of drone operations, which include area coverage (covering location problem), search operations, routing, data gathering and recharging in a wireless sensor network, allocating communication links, and flying wireless ad-hoc networks. Among the surveyed articles, 48% of the articles address routing operation, followed by allocating communication links and area coverage. Only a few papers dealt with routings in disaster operations [23]. A recent study on vehicle routing for post-disaster assessment has been conducted by Bruni, et al., [28], who developed a selective routing problem to minimize the sum of arrival times (total latency) at all nodes.

Concerning routing problems in humanitarian operations, two issues have to be addressed simultaneously: reaching the target and minimizing the operation time. For post-disaster assessment operations when the gathered information depends on the arrival time at the node, the total latency should also be considered. Two classes of routing problems, i.e., the traveling salesman problem (TSP) and the vehicle routing problem, (VRP) have commonly been used in the literature. TSP computes the shortest possible route that visits all targets and returns to the starting position, whereas VRP assigns a set of vehicles that visits all targets in such a way that the routings with the minimum completion time or total operation cost are achieved. According to Viloria, et al., [3], 69% of the surveyed articles consider the routing problem for drones only, whereas 31% introduced the combined system of drones and a ground vehicle. The combined system between the drone and the ground vehicle is motivated by the drone limitations such as limited battery capacity (which constraints operating time) and limited load capacity.

The deployment of the combined system of drone and ground vehicle has accordingly given rise to new variants of TSP and VRP [29]. For a collaborative operation between drone and vehicle, the flying sidekick traveling salesman problem (FSTSP) which determines optimal assignment for a drone working in tandem with a ground vehicle, was developed by Murray and Chu [30]. FS-TSP models the drone that executes an operation to serve targets from the ground vehicle, which traverses at a node, and the drone is retrieved at a rendezvous node. Other TSP variants such as traveling salesman problem with a drone (TSP-D) [31], TSP with multiple drones (TSP-mD) [32], and multiple TSPD (mTSPD) [33] have also been developed. Luo, et al., [34] have further expanded mTSPD to the multi-visit traveling salesman problem with multi-drones (mTSP-mD).

Similarly, the VRP also has many extended variants such as VRP-D, which models a set of ground vehicles equipped with drones to serve the targets [35]. The drone can be launched and picked up at the depot or any of the target nodes. VRPTW-D (vehicle routing problem with time window with a drone) was developed from VRPD by addressing the time window and was even further developed to CVRPTW-D (capacitated vehicle routing problem with time window with a drone) to include payload capacity of drones so that each drone cannot deliver goods exceeding its capacity. Motivated by reducing the number of drones, green vehicle routing problem with drone (GVRP-D) was also proposed by Coelho, et al., [36] and extended to multi-trip vehicle routing problem with drone (MTVRP-D) to embrace the load capacity of drone [37].

The general routing problem (GRP) has also been implemented by Oruc and Kara [21] to deal with routing for post-disaster assessment. Oruc and Kara [21] implemented a bi-objective approach to maximize total value added to assess both road segments (arcs) and population centers (nodes). The GRP differs from VRP in that VRP considers node routing problems, whereas GRP includes both node routing problems (in which assessment targets are located at nodes) and arc routing problems (in which assessment targets are located at arcs on a directed network). As a result, the developed model by Oruc and Kara [21] did not require all nodes or edges to be traversed. Some nodes that lie within the drone’s angular point of view are not visited, and some nodes can be visited more than once.

Another approach, i.e., the two-echelon vehicle routing problem (2EVRP) was developed from the supply chain/logistics field [38]. The 2EVRP is a special case of the two-echelon location routing problem in which the location of all nodes is known. The 2EVRP can be deployed to solve the hybrid aerial and ground vehicle routing by segmenting the system into two echelons. The first echelon refers to the route from starting point (depot) to the stopover point or between stopover points, whereas the second echelon refers to the route from a stopover point to target points. The 2EVRP aims to find a set of first and second echelon routes so that the targets are visited with minimum time/cost. The approach simultaneously optimizes the routes of both a ground vehicle (first echelon) and drones (second echelon). The 2EVRP has been widely used in city logistics/commercial last-mile delivery such as Perboli, et al., [39]. Within the context of humanitarian operations, Luo, et al., [22] have developed the two-echelon ground vehicle and unmanned aerial vehicle routing problem (2E-GU-RP). The difference between the 2E-GU-RP and the developed model in the present study is that the 2E-GU-RP models the drone to be launched at a depot or stopover point and picked up at rendezvous node, whereas 2EVRPA and 2ECoVRP models the drone to be launched and picked up at a stopover point where the drone is launched.

Table 1. Modeling approaches for the hybrid aerial and ground vehicle routing.

	Modelling Approach	Flying Sidekick Travelling Salesman Problem (FS-TSP) [30]	Traveling Salesman Problem with Drone (TSP-D) [31]	Traveling Salesman Problem with Multiple Drones (TSP-mD) [32]	Multiple TSPD (mTSPD) [33]	Multiple TSP with Multiple Drones (mTSP-mD) [34]	Vehicle Routing Problem with Drone (VRP-D) [35]	General Vehicle Routing Problem (GVRP) [21]	Two-Echelon Ground Vehicle and Unmanned Aerial Vehicle Routing Problem (2E-GUCRP) [22]	2EVRPA and 2ECoVRP (This Paper)
Features
Routing problem		TSP	TSP	TSP	TSP	TSP	VRP	VRP	VRP	VRP
Number of ground vehicles:		×	×	×				×
One vehicle
Multi vehicles					×	×	×		×	×
Number of drones:		×	×					×
One drone
Multi drones				×	×	×	×		×	×
Non-collaborative operation		×	×				×	×	×	×
Collaborative operation:
One drone works in tandem with one truck
A fleet of drones working in tandem with one truck				×
A fleet of drones working in tandem with many trucks					×	×
Number of echelons:		×	×	×	×	×	×	×
Non-echelon
Two-echelon									×	×
Considering service time at each target		×	×				×	×		×
Retrieving nodes		target nodes	target nodes	target nodes	target nodes	target nodes	target nodes	target nodes	stopover nodes	stopover nodes
Application		parcel delivery	parcel delivery	parcel delivery	parcel delivery	parcel delivery	parcel delivery	post-disaster assessment	intelligence, surveillance, reconnaissance (ISR)	post-disaster assessment

compares various modeling approaches for the hybrid aerial and ground vehicle.

3. Mathematical Models of the Hybrid Aerial and Ground Vehicle Routing

The 2EVRP is considered a suitable approach to solve the routing problem for post-disaster assessment because the approach allows for more efficient and flexible operations. The approach considers a circumstance when a ground vehicle collaborates with a drone. The 2EVRP can be presented as a two-echelon distribution network. The first echelon considers the ground-space network traveled from a depot to stopover points. The ground vehicle travels on a road network near the reconnaissance area and keeps launching and recycling the drone (first echelon). The second echelon considers the aerial-space network where the drone is launched from the stopover point to the target point(s). The target point is an area where the drone conducts mapping operations and collects the target information and then returns to the ground vehicle before the battery runs out (second echelon).

The first variant, the 2EVRPA, applies the routing problem in the first echelon and assignment problem as a sub-problem of the routing problem. The second variant, the 2ECoVRP, includes the collaborative operation between the ground vehicle and drone (implementing routing problems for both echelons). The models allow the drone to visit multiple targets. The schematic diagrams of both the 2EVRPA and the 2ECoVRP are presented in Figure 1.

3.1. Problem Description and Assumptions

The 2EVRP considers a set of mapping points as targets, each of which must be visited exactly once by a drone. All the mapping points are inaccessible by a ground vehicle. The ground vehicle mounted with a drone departs from a depot and visits stopover points. The drone is launched at the stopover point, visits its mapping point(s), and returns to the launching point to swap its battery for the next flight. The 2EVRPA approach models that once the drone finishes the mapping at the respective target/mapping points, it will return to the launching point although it has still enough battery capacity to perform another mapping process. In contrast, the 2ECoVRP approach allows the drone to visit more than one target/mapping point and will return to the launching point once the battery capacity is insufficient for the next mapping operation. Both models aim to minimize the total operation time for post-disaster assessment.

The assumption used in the models are the following:

• The locations and service times for all mapping points are known;
• All mapping points are visited only once by the drone;
• The locations of the stopover points are known;
• The ground vehicle can only traverse on the road network;
• Each road arc traversed by the ground vehicle corresponds to a flight route of the drone;
• The drone has limited operation time, which is known;
• The time required for battery swap is negligible.

Based on the above description, the notations used in the mathematical models are shown below.

Sets:
$S$	Set of all nodes in the first echelon (includes depot)
$S_0$	Set of all stopover points (without depot)
$C$	Set of all mapping points
$V$	Set of all vehicles in the first echelon
Parameters:
$C{s}_{ij}$	Travel time of vehicle from gathering point or stopover point $i$ to stopover point $j$ (includes depot) in the first echelon
$C{c}_{ij}$	Travel time of drone from stopover point $i$ to mapping point $j$ in the second echelon
$D_i$	Required mapping time for capturing the area of a mapping point $i$
$Q{d}_v$	Maximum flight time of the drone
$L$	Sufficiently large number
Decision variables
$x_{ijv}$	A binary variable that equals one if vehicle v is traveled from node i to j where the edge $\left(i,j\right)\in E1$; otherwise zero
$z_i$	A binary variable that equals one if the stopover point i is visited where $i\in S$; otherwise zero
$y_{ij}$	A binary variable that equals one if the drone is traveled from node i to j where the edge $\left(i,j\right)\in E2$; otherwise zero
$Q_j$	The total amount of time needed for drone taking pictures in a stopover point j where $j\in S_0$
$T_j$	The total amount of time needed for drone traveling in the second echelon where $j\in S$
$a_{iv}$	The total accumulated amount of information gathering time has been used up to stopover point i using vehicle v where $i\in S$ and $v\in V$
$b_{iv}$	The total accumulated amount of transportation time has been used up to mapping point i using vehicle v where $i\in C$ and $v\in V$

3.2. Two-Echelon Vehicle Routing Problem Combined with Assignment (2EVRPA)

The mathematical model of the 2EVRPA was developed based on the basic model of the two-echelon vehicle routing problem (2EVRP) developed by Perboli, et al., [39], which was then modified to facilitate single point-to-point assignment in the second echelon for the present study. A comparable model has been developed for the city logistics context [40]. The 2EVRPA is formulated as the following.

Objective
$Minimize Z={\displaystyle {\displaystyle \sum }_{i\in S}}{\displaystyle {\displaystyle \sum }_{j\in S}}{\displaystyle {\displaystyle \sum }_{v\in V}}C{s}_{ij}{x}_{ijv}+{\displaystyle {\displaystyle \sum }_{i\in C}}{\displaystyle {\displaystyle \sum }_{j\in S_0}}C{c}_{ij}{y}_{ij}+{\displaystyle {\displaystyle \sum }_{i\in S_0}}{\displaystyle {\displaystyle \sum }_{j\in C}}{D}_i{y}_{ij}$		(1)
Subject to
${\displaystyle {\displaystyle \sum }_{i\in S}}{\displaystyle {\displaystyle \sum }_{v\in V}}{x}_{ijv}=z_j$	$\forall j\in S_0$	(2)
${\displaystyle {\displaystyle \sum }_{j\in S}}{\displaystyle {\displaystyle \sum }_{v\in V}}{x}_{ijv}=z_i$	$\forall i\in S_0$	(3)
${\displaystyle {\displaystyle \sum }_{i\in S}}{x}_{ilv}-{\displaystyle {\displaystyle \sum }_{j\in S}}{x}_{ljv} =0$	$\forall l\in S_0$ $\forall v \in V$	(4)
${\displaystyle {\displaystyle \sum }_{j\in S_0}}{x}_{1jv}\le 1$	$\forall v \in V$	(5)
${\displaystyle {\displaystyle \sum }_{i\in S_0}}{x}_{i1v}\le 1$	$\forall v \in V$	(6)
${\displaystyle {\displaystyle \sum }_{i\in C}}{D}_i{y}_{ij}\le Q_j$	$\forall j\in S_0$	(7)
$Q_j\le L·z_j$	$\forall j\in S_0$	(8)
${\displaystyle {\displaystyle \sum }_{i\in C}}C{c}_{ij}· y_{ij}\le T_j$	$\forall j\in S_0$	(9)
$T_j\le L·z_j$	$\forall j\in S_0$	(10)
${\displaystyle {\displaystyle \sum }_{i\in S_0}}{y}_{ij}=1$	$\forall i\in C$	(11)
$a_{iv}+Q_j-L·\left(1-x_{ijv}\right)\le a_{jv}$	$\forall i\in S$ $\forall j\in S_0$ $\forall v\in V$	(12)
$b_{iv}+T_j-L·\left(1-x_{ijv}\right)\le b_{jv}$	$\forall i\in S$ $\forall j\in S_0$ $\forall v\in V$	(13)
${\displaystyle {\displaystyle \sum }_{i\in S}}{\displaystyle {\displaystyle \sum }_{j\in S_0}}{x}_{ijv}· Q_j+{\displaystyle {\displaystyle \sum }_{i\in S}}{\displaystyle {\displaystyle \sum }_{j\in S_0}}{x}_{ijv}· T_j \le Q{d}_v$	$\forall v \in V$	(14)

The objective function (1) is to minimize the total operation time, comprising the traveling time of the ground vehicle to carry the drone, the flight time of the drone, and the mapping time of the drone. Constraints (2) and (3) limit that the nodes can only be visited once. Constraint (4) ensures that the vehicle moves consecutively. Constraints (5) and (6) ensure that the ground vehicle returns to the depot once the operation is completed. Constraints (7) and (8) evaluate the total mapping time at each stopover point. Constraints (9) and (10) are used to calculate the traveling time of the drone in the second echelon. Constraint (11) ensures that the mapping point is only visited once by the drone. Constraints (12) and (13) are used for validating the continuity of both ground vehicle’s routes and the drone assignments. Constraint (14) ensures that the drone’s flight time does not exceed the maximum flight time determined by the drone’s battery capacity.

3.3. Two Echelon Collaborative Vehicle Routing Problem (2ECoVRP)

The mathematical model of the 2ECoVRP was constructed by modifying the 2EVRPA model to include the collaborative operation between the ground vehicle and drone operations. The model allows the drone to visit multiple targets in one go (routing problem). The 2ECoVRP model is formulated as below.

Sets:
$K$	Number of the drone routes in the second echelon
Decision variables
$y_{ijk}$	A binary variable that equals one if drone route k is traveled from node i to j where the edge $\left(i,j\right)\in E2$; otherwise zero
$Q{K}_{jk}$	The total amount of time needed for drone mapping in a stopover point j and every routes k where $j\in S_0$ and $k\in K$
$T_{jk}$	The total amount of time needed for drone traveling in the second echelon where $j\in S_0$ and $k\in K$
$U_j$	The total amount of time needed for drone traveling in the second echelon where $j\in S_0$
$b_{jk}$	The total accumulated amount of transportation time has been used up to mapping point j in route k where $j\in S_0{\displaystyle \cup }C$ and $k\in K$
$c_{jk}$	The total accumulated amount of information gathering time has been used up at stopover point j in route k where $j\in S_0{\displaystyle \cup }C$ and $k\in K$
Objective
$Min Z={\displaystyle {\displaystyle \sum }_{i\in S}}{\displaystyle {\displaystyle \sum }_{j\in S}}{\displaystyle {\displaystyle \sum }_{v\in V}}C{s}_{ij} x_{ijv}+{\displaystyle {\displaystyle \sum }_{i\in S_0{\displaystyle \cup }C}}{\displaystyle {\displaystyle \sum }_{j\in S_0{\displaystyle \cup }C}}{\displaystyle {\displaystyle \sum }_{k\in K}}C{c}_{ij} y_{ijk}+{\displaystyle {\displaystyle \sum }_{j\in S_0{\displaystyle \cup }C}}{\displaystyle {\displaystyle \sum }_{i\in C}}{\displaystyle {\displaystyle \sum }_{k\in K}}{D}_i y_{ijk}$		(15)
Subject to
${\displaystyle {\displaystyle \sum }_{i\in S}}{\displaystyle {\displaystyle \sum }_{v\in V}}{x}_{ijv}=z_j$	$\forall j\in S_0$	(16)
${\displaystyle {\displaystyle \sum }_{j\in S}}{\displaystyle {\displaystyle \sum }_{v\in V}}{x}_{ijv}=z_i$	$\forall i\in S_0$	(17)
${\displaystyle {\displaystyle \sum }_{i\in S}}{x}_{ilv}-{\displaystyle {\displaystyle \sum }_{j\in S}}{x}_{ljv}=0$	$\forall l\in S_0$ $\forall v \in V$	(18)
${\displaystyle {\displaystyle \sum }_{j\in S_0}}{x}_{1jv}\le 1$	$\forall v \in V$	(19)
${\displaystyle {\displaystyle \sum }_{i\in S_0}}{x}_{i1v}\le 1$	$\forall v \in V$	(20)
${\displaystyle {\displaystyle \sum }_{j\in S_0{\displaystyle \cup }C}}{\displaystyle {\displaystyle \sum }_{k\in K}}{y}_{ijk}=1$	$\forall j\in C$	(21)
${\displaystyle {\displaystyle \sum }_{j\in S_0{\displaystyle \cup }C}}{y}_{ijk}\le 1$	$\forall i\in S_0{\displaystyle \cup }C$, $\forall k\in K$	(22)
${\displaystyle {\displaystyle \sum }_{j\in S_0{\displaystyle \cup }C}}{y}_{jik}={\displaystyle {\displaystyle \sum }_{j\in S_0{\displaystyle \cup }C}}{y}_{ijk}$	$\forall i\in S_0{\displaystyle \cup }C$ $\forall k\in K$	(23)
${\displaystyle {\displaystyle \sum }_{i\in S_0{\displaystyle \cup }C}}{y}_{ilk}-{\displaystyle {\displaystyle \sum }_{j\in S_0{\displaystyle \cup }C}}{y}_{ljk}=0$	$\forall l\in C$ $\forall k \in K$	(24)
$Q{K}_{jk}\ge c_{ik}· y_{ijk}$	$\forall j\in S_0$ $\forall i\in C$ $\forall k \in K$	(25)
$Q_j\ge {\displaystyle {\displaystyle \sum }_{k\in K}}Q{K}_{jk}$	$\forall j\in S_0$	(26)
${\displaystyle {\displaystyle \sum }_{j\in S_0}}{Q}_j={\displaystyle {\displaystyle \sum }_{i\in C}}{D}_i$		(27)
$Q_j\le L·z_j$	$\forall j\in S_0$	(28)
$b_{ik}=0$	$\forall i\in S_0$ $\forall k \in K$	(29)
$c_{jk}=0$	$\forall j\in S_0$ $\forall k \in K$	(30)
$c_{jk}\ge 0$	$\forall j\in C$ $\forall k \in K$	(31)
$T_{jk}\ge b_{ik}+C{c}_{ij}$	$\forall j\in S_0$ $\forall i \in C$ $\forall k \in K$	(32)
$U_j\ge {\displaystyle {\displaystyle \sum }_{k\in K}}{T}_{jk}{z}_j$	$\forall j\in S_0$	(33)
$a_{iv}+Q_j-L\left(1-x_{ijv}\right)\le a_{jv}$	$\forall i\in S$ $\forall j\in S_0$ $\forall v\in V$	(34)
$b_{ik}+C{c}_{ij}-L\left(1-y_{ijk}\right)\le b_{jk}$	$\forall i\in S_0{\displaystyle \cup }C$ $\forall j\in C$ $\forall k\in K$	(35)
$c_{ik}+D_j-L\left(1-y_{ijk}\right)\le c_{jk}$	$\forall i\in S_0{\displaystyle \cup }C$ $\forall j\in C$ $\forall k\in K$	(36)
${\displaystyle {\displaystyle \sum }_{i\in S}}{\displaystyle {\displaystyle \sum }_{j\in S_0}}{x}_{ijv}\left(Q_j+U_j\right)\le Q{d}_v$	$\forall v \in V$	(37)

The objective function (15) is to minimize the total operation time consisting of the traveling time of the ground vehicle to carry the drone, the flight time of the drone, and the mapping time of the drone. Constraints (16) and (17) ensure that the nodes can only be visited once. Constraint (18) guarantees that the vehicle moves sequentially. Constraints (19) and (20) ensure that the ground vehicle should return to the depot once the operation is completed. Constraint (21) makes sure all mapping points are visited once, and constraint (22) ensures that the drone visits the mapping points only once. Constraint (23) validates the second echelon vehicle routes and constraint (24) ensures that the drone moves consecutively. Constraint (25) denotes the duration spent at mapping points in the visited route k. Constraint (26) denotes the total accumulated mapping time of each route k. Constraint (27) denotes the accumulated total of mapping time at each stopover point as the accumulated time required to record in that area. Constraint (28) calculates the accumulated mapping time of each mapping point at the stopover point and indicates the stopover point has unlimited capacity. Constraint (29) ensures the total travel time of the drone starts from zero. Constraint (30) ensures the calculation of mapping time starts from zero, and Constraint (31) denotes the mapping duration of the mapping points on each route k, which can be more than zero due to the mapping operation. Constraints (32) and (33) calculate the accumulated travel time of the drone and the total travel time of all drones from the stopover point, respectively. Constraints (34) and (35) are used to validate the continuity of both the ground vehicle’s and the drone’s routing. Constraint (36) validates the mapping process, whereas Constraint (37) limits the flight time of the drone according to its battery capacity, ensuring that it operates below the maximum flight time.

3.4. Tabu Search Algorithm

An optimization problem can be solved by either exact methods or heuristics methods. Exact methods are suitable for solving small-scale problems. However, they may have difficulties in solving a large-scale problem in a reasonable computational time. Heuristics methods, specifically metaheuristics, in contrast, are tailored to find a good solution to large-size and complex problems. Tabu search (TS) is one of the metaheuristics algorithms that was discovered by Glover and formulated in 1989 [41]. Tabu search is an algorithm that incorporates “memory” of the searching history or moves recently applied known as tabu list, which makes it able to search the space economically and effectively. The TS algorithm can avoid the previously visited solution using a tabu list. The tabu list is a short-term set of the solutions, which are changed by the process of moving from one solution to another determined by a set of rules. Despite its adaptive memory, the TS algorithm is also able to accept a worse solution to evade the local optimal trap and can be applied on both discrete and continuous problems. The TS algorithm has been widely used to solve complex problems such as scheduling, quadratic assignment, and routing problems. Several studies, such as Semet and Taillard [42], Chao [43], Scheuerer [44], Boccia, et al., [45], and Venkatachalam, et al., [46], have implemented the tabu search algorithm to solve the two-echelon location-routing problem and the two-echelon routing problem, i.e., truck and trailer routing problem. Furthermore, Venkatachalam, et al., [46] applied the tabu search algorithm to solve the multiple-drone routing problem. It was found that the tabu search algorithm was not only providing a good-quality solution but was also computationally efficient. Similar evidence was also demonstrated by Luo, et al., [34] who solved the multi-visit traveling salesman problem with multi-drones using the tabu search algorithm. The present study thus selects the tabu search algorithm. The tabu search algorithm deployed for the present study is explained in pseudocode as shown in Algorithm 1.

Algorithm 1. Pseudocode of the proposed tabu search algorithm.

Procedure: Tabu Search Input: a set of n neighborhood structure, Nnon-improving, Ncandidate$\mathrm{and} \theta$ Output: Xbest 1.  Construct an initial solution by using random initial solution 2.  Initialize tabu list 3.  Set Iteration = 1, NonImproveCount = 0; 4.  Let Xcurrent←Xbest←Xinit 5.  repeat 6.   repeat 7.    Generate r = random (0,1) to select the neighborhood n 8.    Let Xnew be the candidate solution X obtained by applying the selected neighborhood n(i,j) to Xcurrent 9.    Calculate the objective value of candidate solution X 10.   until the number of candidate solution in the candidate list is equal to Ncandidate 11.   Sort the candidate solutions based on its objective value ascendingly 12.   index = 1; 13.   repeat 14.    Xnew = candidate solutionsp[index] 15.    if Xnew  < Xcurrent || tabu list [n][i][j] = 0 then 16.     Xcurrent←Xnew 17.     Update the tabu list 18.     Go to 21; 19.    index++ 20.   until no more solution in the candidate solutions 21.   if F(Xcurrent) < F(Xbest) then 22.    Xbest←Xcurrent 23.    NonImproveCount = 0 24.   else 25.    NonImproveCount++ 26.   Iteration ← Iteration+1 27.  until the termination criteria is satisfied

The algorithm is initiated by constructing an initial solution (X_init) using a random initial solution and constructing an empty tabu list. Other variables such as current iteration (Iteration = 1) and the number of iterations when the solution is not being improved from the best solution (NonImproveCount = 0) were also initiated. The initial solution is set as the current best solution (X_current) and the best solution found so far (X_best←X_current←X_init). At each iteration, a list of candidate solutions (a list of X_new) is generated by applying a randomly chosen neighborhood n to the current solution (X_current). This procedure is repeated until the number of solutions in the candidate list is equal to N_candidate. The candidate list is then sorted based on its objective value in ascending order. Subsequently, the algorithm to replace the current solution (X_current) with the new solution from the candidate list (X_new) is implemented. The algorithm sets rules on the conditions allowing the replacement, i.e., when the new solution is not generated from a tabu move or it is not in the tabu list, and when the new solution passes the aspiration criteria. The aspiration criteria are set based on the condition that the new tabu solution is better than the current solution. Once the new obtained tabu solution is better than the current solution, the new tabu solution will be adopted. Otherwise, it will be taken out from the neighborhood set. The next step is to repeat the checking process of the neighborhood set. The new and the current solution are exchanged and the tabu movement is updated on the tabu list. The previous best-found solution (X_best) is replaced by the current solution (X_current) once the current solution is better. The count of non-improvement is hence reset (NonImproveCount = 0). However, if the current solution is not better than the best solution, the NonImproveCount is updated. The process is repeated until the termination criteria are met. Otherwise, another iteration for the neighborhood set will be conducted. The proposed TS uses NonImproveCount to decide when the algorithm should be terminated.

4. Computational Results and Analysis

4.1. Case Description: Post-Disaster Assessment for Mount Merapi Eruption

Mount Merapi, the most active volcano in the world, has erupted every four to six years. Within the last 100 years, the largest eruption occurred in 2010, which resulted in 61,154 people being evacuated and 341 casualties [5]. Due to its high occurrence, the regional disaster agency (BPBD–Badan Penanggulangan Bencana Daerah) has formulated a disaster contingency plan, which includes a post-disaster assessment. Due to difficult terrain, the post-disaster assessment becomes challenging. Steep and narrow roads have made some areas inaccessible. However, post-disaster assessment is a crucial operation to support accurate and timely responses.

Therefore, it is necessary to explore an efficient and effective way to support post-disaster assessment. The hybrid ground and aerial vehicle approach seems to be a favorable approach for several reasons. First, the pyroclastic flow of Mount Merapi is barely predictable, making the affected area highly unknown. Second, the pyroclastic flow may emit hazardous and poisonous gas to the surroundings so that it is hazardous for humans. Third, the infrastructure in the affected areas, such as bridges and roads, may be severely damaged, making the areas inaccessible. Using a drone that is mounted on a ground vehicle opens the opportunity to access the disaster-affected areas. Accurate and fast post-disaster assessment is urgently required as the aftermath conditions, such as the number and location of victims and incurred damage of roads and bridges should be known to plan the allocation of resources and personnel as well as to deliver aid efficiently.

Due to the narrow roads and difficult terrain, it appears that motorcycles rather than trucks are the suitable ground vehicles to be deployed for post-disaster assessment for the Mount Merapi eruption. The motorcycle has less carrying capacity than a truck, but the motorcycle has more flexibility to access areas that might be unreachable by truck. The motorcycle brings the drone to map post-disaster conditions. The drone is advantageous for reaching and maneuvering in very narrow terrain such as collapsed buildings. Due to the limited flight time of the drone, the motorcycle increases the mapping coverage. Therefore, the present study uses a motorcycle that is combined with a drone to conduct the post-disaster assessment. The motorcycle is assumed to operate at a constant speed of 45 km/h. The specification for the drones used for the present study is shown in

Table 2. Drone specification.

Attribute	Specification
Mapping Rate	8.125 × 10−5 min/m2
Speed	16 m/s
Ground Speed	45 km/h
Battery Capacity	120 min

.

Figure 2 visualizes the network configuration for assessment operation represented by nodes including depots, stopover points, and target/mapping points for post-disaster assessment, which was implemented using Google MyMaps. It includes 49 nodes comprising nine existing depots, 31 target/mapping points, and 9 stopover points. The detailed location in terms of latitude and altitude as well as areas required to be mapped in each target/mapping point is provided in the Appendix A. QGIS was used to measure the time required to travel from one node to another. The area of the target/mapping points was converted into time units, i.e., minutes, by multiplying the area and the mapping rate of 0.00008125 min/m², as shown in Table 2. The mapping rate was derived from the technical data, indicating the time required by the drones operated at the altitude of 70 m to map the area of 400 m × 400 m is 13 min.

The depots are the starting points of the post-disaster assessment operation. Both motorcycles with mounted drones depart from the depots. The stopover points are the locations where the motorcycles stop and launch the drones for mapping operations. The target/mapping points represent the areas to be mapped by the drones. Each target/mapping point has different areas to be mapped, which then determines the mapping time. The post-disaster assessment operation starts from the depot where the motorcycle that is mounted by a drone is ridden to the stopover points where the drone is launched for mapping. Once the drone nearly reaches the total flight time (running out of its energy), the drone returns to the stopover points. When the mapping operation is accomplished, both the motorcycle and the drone return to the depots. The required total time for those operations is defined as the total operation time.

4.2. Test Instances and Experimental Parameters

To evaluate the consistency of the computational results for various scales, instance benchmarking was conducted. Three groups of instances, based on the procedure developed by Liperda, et al., [47], were conducted.

Table 3. Test instances.

Problem Scale	Range of Nodes	Number of Experiments
Small instance	10–13	12
Medium instance	19–28	27
Large instance	41–49	9

presents the three scale categories, i.e., small-, medium-, and large-size instances with various ranges of nodes to have experimented. A set of experiments were conducted for each scale as shown in Table 3. The total conducted computational experiments was 48. Both the results based on individual experiments and average experiments were then reported.

To compare the performance of the E2VRPA and the 2ECoVRP, the values of the objective function (Obj) of both models were contrasted and quantitatively evaluated using Equation (38) as shown below:

$$\mathrm{Diff}=\frac{\left(Ob{j}_{2EVRPA}-Ob{j}_{2ECoVRP}\right) }{Ob{j}_{2EVRPA}}\times 100\%$$

(38)

To compare the performance of the tabu search algorithm with the exact algorithm, Equations (39) and (40) were used for 2EVRPA and 2ECoVRPA, respectively.

$$\mathrm{Diff}-2\mathrm{EVRPA}=\frac{\left(BestOb{j}_{2EVRPA\left(col. 4\right)}-Ob{j}_{2EVRPA\left(col.2\right)}\right)}{BestOb{j}_{2EVRPA\left(col.2\right)}}\times 100\%$$

(39)

$$\mathrm{Diff}-2\mathrm{ECoVRP}=\frac{\left(BestOb{j}_{2ECoVRP\left(col.6\right)}-Ob{j}_{2ECoVRP\left(col.3\right)}\right)}{Ob{j}_{2ECoVRP\left(col.3\right)}}\times 100\%$$

(40)

The experiments were conducted using the operational parameters that represent the current conditions (such as the number and location of the depots, the stopover points, mapping/target points, motorcycle speed, and drone specifications). Eight drones with the specification shown in Table 2 were used, each of which was mounted on the motorcycle.

The computational experiments were carried out on a computer with an Intel Core i7-7700 3.6 GHz CPU and 16GB RAM. The model was solved using A Mathematical Programming Language (AMPL) with Gurobi solver using branch and bound method and branch and cut method. A time limit of 6 h (21,600 s) was imposed. The tabu search algorithm was implemented using the C# programming language. The results of the tabu search were reported in terms of the best solution and average solution from ten computational runs for each instance.

The proposed tabu search algorithm used four input parameters: a set of n neighborhood structures, N_{non-improving}, N_candidate, and $\theta$. A set of n neighborhood structures was used by the algorithm to explore the solution space. The neighborhood operator was selected randomly among neighborhood $n\in N$ and N = {Swap, Insert, Reverse}. The proposed tabu search used all of the neighborhoods. The parameter N_{non-improving} determined the number of iterations where the best objective function value was not improved consecutively. It was used to terminate the algorithm after it converged to some value in several iterations. The parameter N_candidate was the number of candidate solutions that were generated at each iteration. The parameter $\theta$ (tabu tenure) was the number of iteration where a move generated by a neighborhood operator was prohibited or considered as a tabu move.

The parameter setting was conducted to help the algorithm provide its best performance in terms of solution quality and computational time. The parameters were set following the procedure known as one factor at a time (OFAT). First, a set of candidate parameters were selected based on previous literature. Second, a preliminary run on the given parameter was conducted to reduce the number of candidates. Finally, the remaining parameters were analyzed iteratively by setting one parameter while the rest of the parameters were fixed. The list of candidate parameters for the final selection was as follows:

N_{non-improving}: 25 *|C|*|S|, 50 *|C|*|S|, 100 *|C|*|S|;

N_candidate: 5, 25, 50;

$\theta$ (tabu tenure): |C|*|S|, 2 *|C|*|S|, 5*|C|*|S|.

The notation |C| and |S| represent the number of nodes in the first echelon (depot and stopover points) and the number of nodes (mapping points) in the second echelon, respectively. The experiment was conducted by selecting three instances, each of which represented small-, medium-, and large-size instances for both the 2EVRPA and the 2ECoVRP. Each instance was solved ten times and the average performance was reported. Figure 3 shows the parameter calibration of the tabu search algorithm. It indicates that the higher value of N_{non-improving} led the algorithm to take a longer time for termination. However, there is no significant improvement toward the objective value when changing the value of N_{non-improving} from 50 *|C|*|S| to 100 *|C|*|S|. Using a similar approach, the rest of the parameters, i.e., N_candidate and $\theta$ were analyzed and selected. The selected parameters for the tabu search algorithm were as follows: N_{non-improving} = 50 *|C|*|S|, N_candidate = 25, and $\theta$ = |C|*|S|.

4.3. Performance Comparison of the 2EVRPA and the 2ECoVRP

The performance of the 2EVRPA and the 2ECoVRP in solving the routings of both ground vehicle (i.e., motorcycle) and the drone in the post-disaster assessment of the Mount Merapi eruption was measured using two performance indicators, i.e., the total operation time (which is the value of the objective function) and the computational time to obtain the objective value.

Table 4. Computational results for small-size instances.

	AMPL		TS				2EVRPA vs. 2ECoVRP	AMPL vs. TS		AMPL		TS
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
No.	Obj 2EVRPA	Obj 2ECoVRP	Best Obj 2EVRPA	Avg Obj 2EVRPA	Best Obj 2ECoVRP	Avg Obj 2ECoVRP	Diff (%)	Diff-2EVRPA (%)	Diff-2ECoVRP (%)	Time 2EVRPA (s)	Time 2ECoVRP (s)	Time 2EVRPA (s)	Time 2ECoVRP (s)
1	86.90	76.37	86.90	86.99	76.37	76.39	12.12	0.00	0.00	0.06	0.22	2.22	2.37
2	89.65	81.45	91.18	91.76	81.45	81.70	9.15	1.71	0.00	0.13	5.23	3.36	3.35
3	163.69	141.60	163.69	164.08	141.60	155.73	13.49	0.00	0.00	0.16	18.61	3.68	3.51
4	172.72	138.90	173.94	174.25	139.85	139.99	19.58	0.71	0.68	0.25	1.99	2.62	2.57
5	166.12	138.36	166.12	211.57	138.36	181.40	16.71	0.00	0.00	0.13	0.77	2.60	2.97
6	171.25	138.53	171.25	184.72	142.62	142.75	19.11	0.00	2.95	0.08	0.74	2.24	2.27
7	89.00	77.84	89.00	89.02	77.84	77.89	12.54	0.00	0.00	0.09	0.39	1.94	2.31
8	94.60	83.99	96.13	97.02	83.99	84.49	11.22	1.62	0.00	0.06	16.34	3.08	3.26
9	160.55	140.71	160.55	161.12	140.71	152.12	12.36	0.00	0.00	0.16	19.69	3.41	3.55
10	162.96	132.99	162.96	163.08	132.99	132.99	18.39	0.00	0.00	0.19	0.84	2.56	2.53
11	155.14	134.24	155.14	200.56	134.24	170.41	13.47	0.00	0.00	0.09	0.53	2.49	2.57
12	160.27	131.67	160.27	174.50	131.67	138.69	17.84	0.00	0.00	0.08	0.55	2.21	2.20
Avg	139.40	118.05	139.76	149.89	118.47	127.88	14.66	0.34	0.30	0.12	5.49	2.70	2.79

,

Table 5. Computational results for medium-size instances.

	AMPL		TS				2EVRPA vs. 2ECoVRP	AMPL vs. TS		AMPL		TS
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
No.	Obj 2EVRPA	Obj 2ECoVRP	Best Obj 2EVRPA	Avg Obj 2EVRPA	Best Obj 2ECoVRP	Avg Obj 2ECoVRP	Diff (%)	Diff−2EVRPA (%)	Diff−2ECoVRP (%)	Time 2EVRPA (s)	Time 2ECoVRP (s)	Time 2EVRPA (s)	Time 2ECoVRP (s)
1	180.11	164.52	175.41	183.25	158.91	159.44	9.41	−2.61	−3.41	21,601.31	21,601.95	8.05	6.72
2	192.11	163.76	174.94	190.07	155.78	157.99	10.95	−8.94	−4.87	21,602.23	21,600.94	7.73	6.32
3	193.29	175.49	179.40	189.20	161.83	164.96	9.79	−7.19	−7.78	21,603.22	21,601.95	7.66	6.21
4	207.47	180.73	191.42	208.88	173.94	177.03	9.13	−7.74	−3.76	21,601.44	21,601.44	7.38	5.94
5	198.63	165.74	187.98	198.98	164.23	167.48	12.63	−5.36	−0.91	21,601.31	21,601.31	6.86	6.34
6	232.79	205.85	225.03	231.91	205.85	210.69	8.52	−3.33	0.00	21,601.44	21,601.44	7.34	7.06
7	185.61	158.38	167.77	180.09	153.51	158.84	8.50	−9.61	−3.07	21,601.31	21,601.31	7.61	6.31
8	178.69	155.59	173.63	182.41	151.38	154.36	12.81	−2.83	−2.71	21,602.62	21,602.23	7.79	6.07
9	175.79	151.25	168.09	182.79	150.60	152.85	10.41	−4.38	−0.43	21,603.22	21,603.22	7.69	6.23
10	391.32	367.87	379.31	388.81	350.73	372.09	7.53	−3.07	−4.66	21,601.31	21,602.01	20.31	15.23
11	380.40	347.79	386.08	386.75	347.79	361.78	9.92	1.49	0.00	21,601.31	21,602.23	20.81	17.60
12	383.32	354.80	378.20	393.33	347.81	355.83	8.04	−1.34	−1.97	21,269.94	21,602.65	23.32	15.69
13	451.51	380.86	409.38	446.93	377.22	379.86	7.86	−9.33	−0.96	21,602.01	21,603.38	21.59	15.94
14	411.13	374.16	404.27	430.90	368.33	372.34	8.89	−1.67	−1.56	21,604.54	21,602.62	18.58	16.35
15	476.11	408.08	471.97	487.89	408.08	422.95	13.54	−0.87	0.00	21,600.94	21,603.22	22.69	17.70
16	402.78	348.21	386.72	399.73	345.08	354.20	10.77	−3.99	−0.90	21,601.95	21,601.31	21.05	17.85
17	391.01	350.00	387.70	409.21	349.16	357.24	9.94	−0.85	−0.24	21,604.32	21,269.94	26.22	17.69
18	383.83	333.55	371.88	386.17	333.55	355.47	10.31	−3.11	0.00	21,604.13	21,602.01	20.58	17.58
19	631.50	456.04	629.44	630.57	456.56	462.58	27.47	−0.33	0.11	21,604.54	21,604.54	22.37	15.24
20	623.24	446.97	613.98	657.09	444.37	454.62	27.62	−1.49	−0.58	21,603.22	21,604.32	21.58	15.14
21	640.50	453.18	629.21	637.69	453.18	466.01	27.98	−1.76	0.00	21,603.86	21,604.13	24.30	15.27
22	677.78	479.61	671.08	709.99	479.42	488.94	28.56	−0.99	−0.04	21,604.13	21,604.54	19.11	14.88
23	685.06	490.96	662.80	714.08	490.96	498.56	25.93	−3.25	0.00	21,604.54	21,603.22	18.24	15.73
24	750.24	547.36	750.24	775.73	538.75	558.38	28.19	0.00	−1.57	21,603.22	21,603.86	20.02	16.04
25	618.28	451.32	618.02	628.00	451.32	453.83	26.97	−0.04	0.00	21,603.86	21,603.83	21.29	15.29
26	631.88	461.93	624.06	633.28	461.93	473.39	25.98	−1.24	0.00	21,603.54	21,603.54	18.69	14.46
27	602.32	467.87	602.32	605.79	438.85	461.51	27.14	0.00	−6.20	21,602.76	21,602.76	19.10	15.81
Avg	417.66	334.88	408.16	424.80	330.34	339.01	15.73	−3.10	−1.69	21,590.45	21,590.37	16.59	12.84

and

Table 6. Computational results for large-size instances.

	AMPL		TS				2EVRPA vs. 2ECoVRP	AMPL vs. TS		AMPL		TS
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
No.	Obj 2EVRPA	Obj 2ECoVRP	Best Obj 2EVRPA	Avg Obj 2EVRPA	Best Obj 2ECoVRP	Avg Obj 2ECoVRP	Diff (%)	Diff−2EVRPA (%)	Diff−2ECoVRP (%)	Time 2EVRPA (s)	Time 2ECoVRP (s)	Time 2EVRPA (s)	Time 2ECoVRP (s)
1	683.68	634.86	678.35	702.25	590.33	596.37	12.98	−0.78	−7.01	21,602.23	21,603.22	30.82	60.35
2	653.22	635.50	625.83	672.99	582.26	590.90	6.96	−4.19	−8.38	21,602.65	21,601.31	30.41	60.04
3	652.74	597.97	622.32	675.80	564.57	597.07	9.28	−4.66	−5.59	21,603.38	21,269.94	30.49	59.48
4	741.25	637.65	676.10	795.59	601.53	632.87	11.03	−8.79	−5.66	21,602.62	21,602.01	30.58	60.48
5	696.21	662.33	682.82	765.80	624.68	632.74	8.51	−1.92	−5.68	21,601.09	21,604.54	30.91	62.28
6	781.77	699.89	757.63	820.17	689.84	699.30	8.95	−3.09	−1.44	21,601.06	21,604.32	30.40	61.99
7	670.33	601.13	666.95	802.21	585.85	595.20	12.16	−0.50	−2.54	21,600.95	21,603.38	29.95	60.02
8	725.41	609.23	648.33	689.92	609.07	623.49	6.06	−10.63	−0.03	21,601.05	21,602.62	29.57	60.68
9	772.83	601.16	618.51	688.58	578.43	585.56	6.48	−19.97	−3.78	21,602.62	21,601.09	29.31	59.60
Avg	708.60	631.08	664.09	734.81	602.95	617.05	9.16	−6.06	−4.46	21,601.96	21,565.83	30.27	60.55

show the computational results for small-size, medium-size, and large-size instances, respectively. The first column of the tables indicates the identity of the instance. The values of the objective function based on the exact algorithm using AMPL for the 2EVRPA and the 2ECoVRP are reported in the second and the third columns, respectively. Concerning the tabu search algorithm, the values of the best and average objective function for the 2EVRPA are shown in the fourth and fifth columns, whereas those for the 2ECoVRP are shown in the sixth and seventh columns. It is worth noting that for the small-size instances, the comparison evaluation between the 2EVRPA and the 2ECoVRP is based on column 2 (2EVRPA) and column 3 (2ECoVRP). However, for the medium-size and large-size instances, the evaluation is based on column 4 (2EVRPA) and column 6 (2CoVRP) because the exact algorithm cannot obtain optimality. The higher the difference, the better performance of 2ECoVRP compared with that of 2EVRPA. In addition to the objective value, the computational time was also evaluated to measure the efficiency of the two models. The computational times for both the exact algorithm and the tabu search algorithm were also recorded in columns 11–14.

The findings for small-size instances indicate that the 2ECoVRP performs 14.66% better than the 2EVRPA in terms of the total operation time. However, the 2ECoVRP requires a slightly longer computational time of about five seconds, which is considered insignificant. It is also interesting to highlight that the tabu search algorithm produces results as good as the exact algorithm. The average difference of the total operation time resulting from the tabu search compared to the exact algorithm is 0.34% and 0.30% for the 2EVRPA and the 2ECoVRP, respectively. The average computational time using the tabu search algorithm is about 2.7 s, which is insignificant. The findings imply that the 2ECoVRP approach solved by the tabu search performs reasonably well in small-size instances.

The findings in Table 5 and Table 6 indicate that the exact algorithm was not able to reach the optimal results within 6 h of computation time. When the operation time obtained by the tabu search algorithm is compared with the best feasible solution by the exact algorithm, it is found that the tabu search algorithm performs better, as indicated by the 3.1% and 1.7% lower operation time for the 2EVRPA and the 2ECoVRP, respectively. It becomes obvious for large-size instances in which the operation time is lower by as much as 6.06% and 4.46% for the 2EVRPA and the 2ECoVRP, respectively. When all target points are considered, it would take about 10 h for the assessment.

Furthermore, the findings demonstrate that the 2ECoVRP approach is better than the 2EVRPA as shown by the lower average of the total operation time by 15.73% for medium-size instances and 9.16% for large-size instances. Regarding the computational time, the 2ECoVRP generally requires a longer computational time than the 2EVRPA. When all nodes are being evaluated, the computational time for the 2ECoVRP is twice as long as that for the 2EVRPA. Nevertheless, the average computational time of the 2ECoVRP for large-size instances, which is about 60 s or one minute, is relatively short and acceptable to be used in practice/field.

Figure 4 shows the comparative analysis for the two models, i.e., the 2EVRPA and the 2ECoVRP, and the solution methods, i.e., exact algorithm and tabu search algorithm. It is observed that the 2ECoVRP approach performs better than 2EVRPA for all scales of instances. It is important to note that a trade-off, however, exists between the achieved the best solution, i.e., the average total operation time, and the average computational time. To achieve lower total operation time, it generally requires longer computational time. An exception is, however, observed for medium-size instances in which the computational time for the best total operation time can be achieved with a shorter computational time. The result is of course desirable, but it is important to note that the medium-size instances were generated using 19–28 nodes out of 49 nodes, so the results are highly influenced by the selected node configurations. Detail analysis was further conducted to measure the standard deviation of the best solution for the 2ECoVRP approach. It was found that the standard deviation for medium-size instances was higher, at 126 min, than that of small-size instances. at 26 min, and that of large-size instances, at 33 min, indicating a wider variance of total operation time for medium-size instances.

With respect to the solution method, it seems that the tabu search algorithm is acceptable because it performs almost as well as the exact algorithm, as indicated by a slight difference when compared to the exact algorithm in the small-size instances. Because the tabu search algorithm has been able to solve the medium-size and large-size instances with a good-quality solution in sensible computational time, the tabu search is therefore suggested to be used as the solution method for the 2ECoVRP model.

4.4. Sensitivity Analysis

Sensitivity analysis is conducted to analyze the sensitivity of the total operation time due to changes in the number of drones. Figure 5 demonstrates that a lower number of drones corresponds to a higher total operation time of post-disaster assessment. In contrast, the more drones, the faster the assessment operation of post-disaster assessment. For instance, the addition of three drones, from five drones to eight drones, results in a shorter operation time of 23 min. Nevertheless, it is important to note that more drones are not always coupled with the faster operation. It appears that additional drones after the deployment of eight drones do not lead to a significant reduction of the total operation time. It can be argued that the operation time of post-disaster assessment is also constrained by the given network configuration, particularly the number and location of stopover points as shown in Figure 2, the total areas to be mapped, and the assumption that each motorcycle can only be mounted by one drone.

4.5. The Effect of Drone Configuration on the Operation Time

In addition to network configuration as previously discussed, the operation of hybrid aerial and ground vehicles is also influenced by drone configuration. This subsection, therefore, evaluates the effect of various drone configurations on the total operation time.

Table 7. Drone configurations.

	Configuration	1	2	3	4
Attribute
Flight altitude (meters)		70	150	70	150
Battery capacity (minutes)		120	100	100	90

presents four drone configurations concerning flight altitude and battery capacity because these two parameters are considered to be influential factors in drone operation. The drone operated at a higher flight altitude has a higher capability to capture wider areas; however, the flight time is lower than that when operated at a lower altitude due to faster depletion of battery capacity. Similarly, the drone with a bigger battery capacity has a longer flight time. The large-size instances were used for the analysis.

Figure 6 demonstrates that the 2ECoVRP performs better than 2EVRPA in all drone configuration settings, as indicated by the lower operation time. It is not surprising because the previous findings have demonstrated the superiority of the 2ECoVRP model. It can be argued that multi-visits in the 2ECoVRP helps to reduce the operation time significantly.

The results indicate that the best drone configuration is when the drone with a battery capacity of 120 min and is operated at a flight altitude of 70 m. It is worth noting that when the battery capacity of the drone is reduced by 20 min, the total operation time is longer by 8.65 min and 51.5 min for 2ECoVRP and 2EVRPA, respectively (see Configuration 3). Similarly, when the flight altitude was set at about twice as high, the operation took an additional time of 28.88 min for 2ECoVRP and a lower time by 21.73 min for 2EVRPA (see Configuration 2). When comparing the result of Configuration 2 with that of Configuration 4, it was found that the reduction of battery capacity by 10 min results in a longer operation time of 29.77 min for 2ECoVRP and a lower operation time by 38.2 min for 2EVRPA. It appears that some inconclusive results were found for the effect of the flight altitude on the operation time and the effect of battery capacity on the operation time for the EVRPA.

Further investigation was conducted to explore the effect of a single operation parameter, i.e., battery capacity, on the best operation time. The investigation was based on Configuration 1 as the best configuration, which was modified by varying its battery capacities. Given the flight altitude of 70 m, Figure 7 indicates that higher battery capacity results in faster operation time. In contrast, using a lower battery capacity results in a significantly longer operation time. The improvement of the total operation time, however, does not seem to be linear. The increased battery capacity of more than 110 min does not lead to a significant reduction of the operation time. This implies that it is important to find the optimal battery capacity that fits a given post-disaster assessment operation so that the assessment operation is efficient.

The aforementioned findings indicate that both flight altitude and battery capacity also contribute to the efficient assessment operation. It implies that further study exploring the optimal setting of drone operation for flight altitude and battery capacity is required.

4.6. Practical Implications

Although sustainability is not a new concept, its implementation within humanitarian operations has only recently been realized. A sustainable system generally rests on the continuous awareness that the objectives of sustainability, i.e., economic (e.g., cost), social (e.g., fairness), and ecological (e.g., emissions) goals are met. Therefore, the present study suggests that planning for post-disaster assessment such as the implementation of the hybrid aerial and ground vehicle system should be conducted not only to facilitate more quick and accurate assessment efficiently but also to reduce resources used.

Concerning the ecological goal, the deployment of the technology, i.e., the drone, for post-disaster assessment helps to reduce fuel consumption, consequently less GHG emissions and lower cost. Studies on life-cycle assessment for drones by Figliozzi [11] and Koiwanit [12] have indicated that drone is one of the most environmentally friendly transportation options. However, the deployment of the drone should also be accompanied by the efficient operation of the drone. The results imply that the efficiency of the hybrid aerial and ground vehicle operation depends on the network configuration, the routings of both vehicles, and the operation parameters such as flight altitude and battery capacity of the drone. To facilitate efficient routing of both the ground vehicle and the drone, the results imply that the 2ECoVRP approach solved by the tabu search algorithm seems to be promising as it only requires one minute to obtain the optimum routings of both vehicles. The findings also imply that operation settings should be explored further for a more efficient operation. For the studied case, it is suggested to deploy a maximum of eight drones with a flight altitude of 70 m and battery capacity of 120 min. The findings can be used by the search and rescue team of the Regional Disaster Agency to improve the efficiency and effectiveness of the post-disaster assessment.

Concerning the social aspect, the hybrid aerial and ground vehicle system allows accessing wider areas, as the drone can access areas that are inaccessible by the ground vehicle, thus facilitating fairness for the remote areas. Hence, the deployment of the hybrid aerial and ground vehicle at its optimal operation not only provides immediate and accurate assessment, which is crucial for emergency response operation, but also facilitates the sustainability aspect of the operations in terms of environmental, economic, and social goals.

5. Conclusions and Future Research

Due to its critical role of post-disaster assessment toward timely and effective response, an innovative and efficient approach to carry out the assessment is required. The present study proposes and evaluates the Two-Echelon Vehicle Routing Problem combined with Assignment (2EVRPA) and Two-Echelon Collaborative Vehicle Routing Problem (2ECoVRP) to evaluate optimal routings of both aerial and ground vehicles, which gives the minimum operation time for post-disaster assessment. The developed models were applied to solve the post-disaster assessment for the Mount Merapi eruption in Yogyakarta, Indonesia, to demonstrate its applicability to a real problem. A set of numerical experiments based on the empirical case were conducted. Both the exact algorithm and the tabu search algorithm were implemented to solve both models.

The findings indicate that the 2ECoVRP performs better than 2EVRPA in terms of the total operation time. It can be argued that the 2ECoVRP facilitates the more efficient operation of the drone by allowing multiple visits in one go. The average operation time to assess all mapping points is about 10 h, which is relatively faster than the current practice. Hence, the 2ECoVRP seems to be a promising approach to evaluate the routing of the hybrid aerial-ground vehicle system for post-disaster assessment. It is worth noting that a trade-off exists between the best solution and computational time. The computational time of the 2ECoVRP is about twice as long as that of the 2EVRPA for large instances. However, as the computational time for the 2ECoVRP is about one minute for large instances, it is considered reasonable and practical. It is also important to note that the tabu search algorithm can obtain optimal solutions that are comparable to those of exact algorithms, i.e., branch and bound as well as branch and cut methods, in a reasonable time. It is therefore suggested to be deployed to solve the 2ECoVRP.

Despite its contribution, some limitations of the present study need to be highlighted. First, the study has not compared the performance of the tabu search algorithm with other heuristics methods. The exploration for a more effective and efficient method to solve for the 2ECoVRP is therefore suggested as a potential avenue for future research. Second, the study has limited numerical analysis for different operational parameters such as flight altitude, drone specifications (e.g., battery capacity, mapping rate, speed), and the network configuration. However, the findings imply that the drone configuration settings (i.e., flight altitude and the drone’s battery capacity) influence the total operation time. Future studies could therefore extend the present study to explore the best solution by considering both network configuration and drone configuration. Last but not least, the present study deals with the deterministic problem. Due to high uncertainty corresponding to humanitarian operations, it is, therefore, necessary to evaluate the effectiveness of the obtained best solution in a dynamic environment setting. A hybrid optimization-simulation method such as that developed by Sopha, et al., [48] is thus suggested for other potential future research to test the effectiveness of the solution to support more effective and efficient post-disaster assessment practices.