Truck-drone team logistics: A heuristic approach to multi-drop route planning

Depósito de Investigación de la Universidad de Sevilla https://idus.us.es/ This is an Accepted Manuscript of an article published by Elsevier in Transportation Research Part C: Emerging Technologies, Vol. 114, on May 2020, available at: https://doi.org/10.1016/j.trc.2020.02.030 © 2020 Elsevier. En idUS Licencia Creative Commons CC BY-NC-ND Truck-Drone team logistics: a heuristic approach to multi-drop route planning Pedro L. Gonzalez-Ra*, David Cancaa, Jose L. Andrade-Pinedab, Marcos Callea, Jose M. Leon-Blancoa a Department of Industrial Engineering and Management Science, School of Engineering, University of Seville, Seville, Spain b Robotics, Vision & Control Group, School of Engineering, University of Seville, Seville, Spain Abstract Recently there have been significant developments and applications in the field of unmanned aerial vehicles (UAVs). In a few years, these applications will be fully integrated into our lives. The practical application and use of UAVs presents several problems that are of a different nature to the specific technology of the components involved. Among them, the most relevant problem deriving from the use of UAVs in logistics distribution tasks is the so-called “last mile” delivery. In the present work, we focus on the resolution of the truck-drone team logistics problem. The problems of tandem routing have a complex structure and have only been partially addressed in the scientific literature. The use of UAVs raises a series of restrictions and considerations that did not appear previously in routing problems; most notably, aspects such as the limited power-life of batteries used by the UAVs and the determination of rendezvous points where they are replaced by fully-charged new batteries. These difficulties have until now limited the mathematical formulation of truck-drone routing problems and their resolution to mainly small-size cases. To overcome these limitations we propose an iterated greedy heuristic based on the iterative process of destruction and reconstruction of solutions. This process is orchestrated by a global optimization scheme using a simulated annealing (SA) algorithm. We test our approach in a large set of instances of different sizes taken from literature. The obtained results are quite promising, even for large-size scenarios. Keywords: drones, delivery, tandem truck-drone, optimization, heuristic 1. Introduction The capability of unmanned aerial vehicles (UAVs or drones) to travel without a link to the road infrastructure is leading to new working models in a variety of logistic scenarios. For instance, in emergency and disaster management situations, UAVs equipped with cameras allow for viewing disaster scenes promptly, collecting critical data for the response teams on the ground. In the healthcare sector, UAVs can also be used to transport diagnostic samples or medical laboratory items between facilities. In all these situations, the use of drones instead of other transportation modes requires dealing with the limited life of the electrical batteries powering each UAV. Generally, the cooperative use of drones with another vehicle acting as a moving base is the most versatile solution, at the cost of the need of determining the most effective locations of rendezvous points. Beyond Preprint submitted to TRANSPORTATION RESEARCH: C. DOI 10.1007/s10479-015-1941-2 * Corresponding author. Email: [email protected] 2 the development of new and better technology solutions for UAVs, a real transference to industry still requires the resolution of certain tactical and operational issues. Currently, most of the civil uses of drones are based on the cooperative employment of a drone and a truck, as has been done for decades in the military field (Garone et al., 2011). Typically, the envisioned applications incorporate an electrical multirotor vehicle which is carried by the truck to the appropriate position, to take off and start its mission in a remotely piloted mode (see e.g., Chang and Lee 2018, Ha et al. 2018, or Yurek and Ozmutlu, 2018). Hence, a moving truck acting as a base allows the enlargement of the action radius of the drone’s operation, at the cost of complicating the determination of landing locations. Hybrid truck-drone systems have also been reported in other practical applications where the crucial feature is the drone’s ability to travel directly between two points of interest. These include: intelligence, surveillance and reconnaissance (ISR) missions to visit a set of locations (Manyam et al., 2018); first healthcare aid, especially in developing countries (Wen et al., 2016) and in emergency aid in order to reduce the worker’s exposure to danger and/or facilitate a quick and flexible access to the locations of interest (Rabta et al., 2018). The latter takes special importance in situations where the transportation network is severely compromised by natural disasters, aimed at collecting crucial information for the efficient development of the emergency mission itself (see e.g., Chowdhury et al., 2017, or Greenwood, 2015). Over the last years, such practical applications have motivated drone-related optimization approaches with objectives around to both, time and cost. The latter, mainly aimed at the degree of difficulty of operations and/or the resources consumption, oftentimes applied throughout the use of cost-drivers as the distance travelled, the number of drones or the number of batteries that a certain solution requires. However, as recently showed in Schermer et al. (2019), the more usual performance indicators in hybrid drone-truck scenarios are time-related: the minimal completion time for serving the set of locations, the minimal time in which the later vehicle retrieval occurs at the end of operations or the saving time compared to service by truck only. Our research covers some of the gaps identified in the conceptualization of new working models around the use of drones in the above-mentioned logistic scenarios. In particular, we focus on the applications where a drone trip is capable of serving several locations (every trip comprises of a launching, the visit to a typically small number of locations and a landing in the ground vehicle at the next rendezvous point). This differs from the commonly accepted hypothesis in last-mile delivery literature which considers that only one single location is visited at each drone’s trip – with the exception of only a few works, see e.g., Wang et al (2019), or Karak and Abdelghany (2019). However, notice that this limitation is not necessarily applicable; for instance, current drone technology allows for picking-up small items –e.g., vaccines, water purification tablets, or medicines- and delivering them one-at-a-time: a sequence of locations are then served by the same drone’s trip. In our approach to the problem, every location stands for a customer that can be served either by drone or by truck. We assume the truck stops at certain customer sites (which are not predefined) and hence, the drone is only allowed to merge with the ground vehicle at these locations. Once landed, the drone always gains a fully charged battery 3 and would then be ready to start a new trip or stay at the truck while the truck carries it to a new services area. The main contributions of our paper can be summarized as follows: • We present a novel truck-drone team logistic (TDTL) mathematical formulation, allowing multi-drop routes for the drone. Under the assumption of infinite travel autonomy for the truck and a limited battery life for the drone, we plan the synchronization events –namely, the sites wherein the drone will gain a fullycharged new batterywith the objective of minimizing the makespan. We assume that the truck-drone team starts from an origin (depot) and has the mission of reaching an end node, visiting the entire set of locations (customers), either by the truck or the drone. All the routes in the problem are open in the sense that are not predefined. • In order to solve real-life instances of the above challenging mixed integer model (MIP), we propose a heuristic method with an innovative coding scheme, a fast and high-quality constructive two-step procedure for obtaining the initial solution and two original local search procedures for improving solutions. We code the solutions with a tuple formed by a node-sequence vector (defining the order of visits) and a resource-type vector containing the type of resource that visits each node (truck, drone or truck-drone) in the node-sequence vector. • We compare the resolution of the instances by using the Gurobi’s Branch-and-Cut algorithm with the heuristic approach in a set of 1080 literature-based instances. The heuristic procedure exhibits a good performance in solving all the instances in a reduced time while providing good solutions. • The developed heuristic allows its transfer to real-life cases, since it is able to obtain solutions of relatively good quality for large-size problems in short computation time. 2. Literature Review The last-mile delivery is by far the most prominent research stream on hybrid truck-drone systems. Amazon, DHL, Google and Alibaba have conceptualized their logistics transportation solution around this new cooperative approach, so that a truck loaded with one or several drones gets near the delivery locations, which are ultimately served by a drone. In view of the new method for crossing the last-mile to the customer sites, the delivery service plan is to be built on the assumption that the truck acts as a mobile depot to pick the parcel up. While some authors give an active role to the trucks in the delivery process in the sense that they are able to deliver the products to the customers (see e.g., Murray and Chu, 2015, or Ha et al., 2018), others assume that the truck does not visit customers (see e.g., Karak and Abdelghany, 2019, or Ferrandez et al., 2016). The former forces the planner to face a synchronization problem, since the truck and the drone need to wait for each other. In other works (see Agatz et al., 2018) such a synchronization issue is mitigated by using a policy of “drone following the same road network as the truck”, although this is somehow contradictory with fully exploiting the last-mile drone delivery advantage. The problem of pairing drones with traditional delivery trucks was first devised by Murray and Chu (2015). These authors study two types of truck-drone delivery problems. 4 The first one, the flying sidekick traveling salesman problem (FSTSP), concerns the case in which the distribution center (DC) is far from the customer locations. In this problem the truck must visit some fixed customers due to the weight of the parcels involved, while the rest of the parcels can be delivered by the drone carried by the truck. The second one is the parallel drone scheduling traveling salesman problem (PDSTSP). This problem corresponds with the situation in which a significant proportion of customers are located within a drone’s flight range from the DC. In this model, drones depart from the depot, deliver parcels to clients and return to the depot. At the same time, the truck services some clients without carrying any drone. Both drone and truck are working simultaneously in parallel. Authors propose simple heuristics to obtain solutions to these problems. Recently, Murray and Raj (2020) have addressed the so-called multiple FSTSP, which considers an arbitrary number of heterogeneous UAVs along with the truck, with attention to real-life issues as it is the scheduling of launching and landings activities. Traditionally, the studies on the transformed last-mile delivery are focused on the routing decisions of truck and drones. The problem that arises is called Travelling Salesman Problem with Drones (TSP-D). Chang and Lee (2018) consider that a single truck carries the drones to some centers from which drones fly to serve the customers. Moshref-Javadi and Lee (2017) use a truck-multi-drone cooperative delivery system to minimize delays in a customer-oriented distribution system. At each stop site, the truck will wait until all drones come back. The truck then carries drones to the next node. Ferrandez et al. (2016), analyze a cooperative truck-drone tandem, considering that the truck follows the route generated by initially solving a traveling salesman problem (TSP). Then, at each truck stop one or more drones can be launched to deliver parcels. In this way, each truck stop is like a hub for drone deliveries. Authors first applied a 𝐾-means clustering algorithm to determine truck stops (by grouping customers) and then applied a genetic algorithm to determine the truck TSP tour. Again, concerning the parallel drone scheduling problem, Ham (2018) extends the original PDSTSP proposed by Murray and Chu (2015) by assuming that drones can perform not only delivery operations but also pick-up activities. After a drone finishes the delivery task, it can fly back to the depot and then begin a new trip, or it can fly to the next client to pick up a returned parcel, in such a way, each drone can visit only two consecutive nodes and no synchronization with the truck occurs at all. Moreover, as in the original PDSTSP, trucks and drones’ routes are decoupled. Boysen et al. (2018) provide an integer model and a complexity analysis of several variants of the PDSTSP (considering robots instead of drones) based on the number of available robots and whether or not they are launched. The paper develops scheduling procedures which determine the truck route along robot depots and drop-off points where robots are launched, such that the weighted number of late customer deliveries is minimized. Customers are only served by robots, and the truck route visits several robot depots where the truck can replenish robots. Authors propose a heuristic approach which is compared against the Gurobi Branch-and-Cut solver for small instances (up to 10 nodes). Wang et al. (2017), Poikonen et al. (2017), Schermer et al. (2019) and Sacramento et al. (2019) have considered a last-mile delivery around a fleet of homogeneous ground vehicles (trucks) equipped with multiple drones, when stating the so-called Vehicle Routing Problem with Drones (VRP-D). Wang et al. (2017) and Poikonen et al. (2017) make a comparison between the amounts of time (best case scenario) that can be saved 5 using a cooperative scenario (truck-drone) when compared with the delivery using only trucks. Sacramento et al. (2019) aim their operational costs minimization approach at the equalization of the working time of truck drivers (all the routes restricted to a maximum duration). Schermer et al. (2019) have provided a comprehensive approach to the makespan minimization VRP-D, firstly incorporating several sets of valid inequalities to its MILP formulation and secondly proposing a matheuristic procedure that exploits the structure of the VRP-D to leveraging the commercial MILP solvers. The use of multiple trucks and multiple UAVs is considered also by Kitjacharoenchai et al. (2019), which is an extension of the FSTSP, but without endurance limitations and launch/delivery time considerations. While UAVs can be launched from and retrieved at different trucks, only one UAV can be launched or retrieved at any customer location. They present a MIP formulation and propose an insertion-based heuristic to solve problems with up to 50 customers. Ulmer and Thomas (2018) study a same-day delivery problem with trucks and drones, in which customer orders come dynamically during a work shift. The involved decisions concern the acceptance or rejection of customer orders for same-day delivery, and, if accepted, the method of delivery, by selecting between a truck or a drone. Pugliese and Guerriero (2017) incorporate time windows constraints for delivery goods. In their experiments, although the use of drones is not economically positive, they conclude that a cooperative strategy can reduce negative environmental impacts and improve service quality. Campbell et al. (2017) formulate and optimize continuous approximation models of hybrid truck-drone delivery, where both trucks and drones make deliveries simultaneously and compare their performance with respect to truck-only delivery scenarios. Wang et al. (2019) propose a MIP for the VRP-D and provide a branch-and- price algorithm for its exact resolution using randomly generated instances up to 15 nodes. We can find approaches in which the role of drones is even more relevant. Dayarian et al. (2017) focused on a same-day delivery problem with drone resupply, where trucks oversee delivery orders and the role of drones is committed to resupplying the trucks. Authors present a heuristic to solve the problem with only one truck and one drone. Furthermore, Cheng et al. (2018) address a distribution problem where there are no trucks at all as a VRP which is further solved using a Branch-and-Cut algorithm. The valid cuts introduced come from the nonlinear (convex) energy function modeling drones' energy consumption (which is influenced by payload and travel distance). Dorling et al. (2017) proposed a linear approximation to the non-linear drone energy consumption, linearly depending on the drone payload and the battery weight, considered as a decision variable. In their model, each drone, departing from a given depot, can perform multiple trips and visit multiple customers per trip, in a similar way than a VRP problem, but considering battery endurance constraints. The model is approximately solved using a SA heuristic for different scenarios up to 500 customers. We refer to Dorling et al. (2017) for the variety of factors that intervene in the maximum operation time or endurance of a flight mission. Henceforth, we concentrate our revision two-fold: i) in the use of a single truck equipped with a single drone, both cooperating to serve a set of locations and ii) in the resolution methods. The tandem truck-drone has been studied for urban contexts where only the drone is responsible for delivering parcels to customers: the truck serves the aim of traveling along the street network and the drone’s purpose is to reach the customer’s doorstep (see 6 Mathew et al., 2015, and bin Othman et al., 2017). Mathew et al. (2015) consider that the truck can wait at a node until the drone returns or go to the next vertex to rendezvous with the drone. The same assumption is considered in bin Othman et al. (2017), which in addition starts from a fixed route for the truck –i.e. predefined for supply issues- and define the routing decisions for the drone. In contrast with both references, in our research we extend the FSTSP model to allow the drone to visit several customers per trip, between two consecutive rendezvous with the truck. Certain similarities arise with Luo et al. (2017), which study an ISR problem with a truck and its carried drone. However, while they force the truck to select from among a set of candidates stopping points that are adequate for rendezvous operations, in our paper the possible rendezvous nodes are not- a-priori defined. As regards the relevant resolution methods published, they mostly focus on the FSTSP. Ponza (2016) proposes a SA heuristic to solve the FSTSP. Bouman et al. (2017) use a dynamic programming (DP) approach for the FSTSP. Marinelli et al. (2017) extend the FSTSP by allowing launching and rendezvous operations along the route arcs. The authors propose a greedy randomized adaptive search procedure for solving this problem. Agatz et al. (2018) extend the FSTSP, allowing the truck the possibility of waiting for the return of the drone at the start node. They propose a route-first cluster-second heuristic method. More recently, Jeong et al. (2019) extend the FSTSP problem to include two practical considerations concerning the drone operations: the existence of non-fly zones and the effect of parcel weight on drone energy consumption. Based on the conducted review, we can state that papers addressing the truck-drone cooperative system mostly assume that during each trip a drone can only visit one customer, except for Ham (2018), Luo et al. (2017) for the FSTSP and Cheng et al. (2018) for the PDSTSP (the latter involving a pure drone-delivery scenario where no truck intervenes at all). Our research work assumes a set of locations to which we need either deliver lightweight relief items or visit for an ISR mission. Like us, Luo et al. (2017) and Wang et al. (2019) assume that both the truck and the drone can serve the customers. Luo et al. (2017) have evolved the FSTSP, when designing a truck-drone cooperative ISR mission over a set of surveillance target locations. Wang et al. (2019) have addressed the routing decisions for trucks and drones in an urban last-mile delivery problem under the following practical operational considerations: i) the existence of a predetermined set of docking hubs –i.e. large areas to carry out a controlled landing -, ii) the possibility for drones of taking-off directly from the depot, from a docking hub or from a truck and, iii) the latter is only allowed at one of the customer locations. However, both works limit the locations where the synchronization among trucks and drones can occur: Wang et al. (2019) schedule it at the docking hubs, whereas Luo et al. (2017) schedule the rendezvous points at truck stops among a preselected subset of candidate locations (a simplification that also appears in Sacramento et al. (2019) and Karak and Abdelghany (2019). On the contrary, in our approach, all the locations are open to be visited either by the drone or by the truck; therefore, all of them are potential points for both stop and rendezvous. Besides, while Sacramento et al. (2019) impose the truck cannot wait for the drone at the same location it was launched, in our approach we get the same result without enforcing this constraint: owing to the active role played by the truck, it would not be of interest that a launch and a rendezvous occur at the same location. Indeed, for us the synchronization issue is crucial, 7 since it turns the problem into finding the chain of customers to be served by each vehicle while specifying the locations where the battery swaps will take place. We next describe the formulation of our TDTL problem, a generalization of the FSTSP model proposed by Murray and Chu (2015) capable of planning multi-visit per drone’s trip in tandem with a truck that actively serves some a priori unknown locations, acting in addition as a ground base (recharging and rendezvous). 3. Problem description This problem consists of finding the routes for a drone and a truck that collaboratively have to visit a set of customers, with different starting and finishing points, as is typically required in practical ISR missions. Since a single drone cannot visit all customers (owing to a short flight range due to the limited battery life), the truck must assume the role of a battery swap station. While this fact has been typically addressed by assigning a mothership role to the truck – with the only mission of carrying the drone, see e.g., Poikonen et al. (2017), and Wang et al. (2017)–, in our study the truck itself participates also in serving customers. Moreover, we use a flexible definition of the truck routes, in contrast to what is often assumed in the literature, where the truck route must be completely pre-defined when determining the operation of the drone as in bin Othman et al. (2017), or partially defined (see e.g., Murray and Chu, 2015, Mathew et al., 2015, or Luo et al., 2017) owing to certain customers requiring a visit by truck. The latter makes some sense in delivery operations (e.g., due to heavy deliveries, as in Sacramento et al., 2019), but not in the case of small items delivery or in the case of reconnaissance missions. In short, in our study we are interested in covering the more general case in which each customer can be served or visited either by truck or by drone. In addressing this more general problem, we make the following assumptions: 1. We disregard the maximum distance the truck can operate (the truck can finish its route without the need of refueling). 2. Drones can only merge with a truck at a customer node (rendezvous locations) and not at intermediate locations. The change of the drone's battery is made at rendezvous points (see figure 1). The truck and drone routes are synchronized at the rendezvous locations; –i.e. if the truck arrives before at the customer, the truck must wait for the drone, and vice versa. Figure 1. Merging procedure at customer locations 8 3. There is no limitation to the number of batteries or swap operations made on the drone, the swap-time being negligible as compared with the total time of the mission. Moreover, since the battery is replaced by a fully charged new one each time the drone and the truck merge, the drone battery is completely full when the drone takes off from the truck. 4. As long as the drone has enough battery, it can visit multiple customers on each trip that ends by merging with the truck at a new customer location. 5. The drone, if needed, can travel a sub-route transported in the truck, in which case, there is no consumption of battery. 6. Each trip is composed of segments interconnecting the locations, measured by the Euclidean distance between them. They are crossed at constants speeds by both the truck and the drone. 7. For the sake of simplicity, the pattern of drone’s energy consumption is assumed to be linear. We do not consider the possible excess of consumption during take-off or landing operations, nor the non-linear pattern which therefore appears as an effect of the load carried. 8. Each customer must be served either by a drone or by a truck. 4. Mathematical formulation This section presents the mathematical model for the described problem. The different sets needed for modelling the problem, the parameters and variables are firstly defined. Then, the formulation is presented, grouping the constraints according to their specific function. 4.1 Notation Sets: 𝐺={𝑁,𝐴} Graph defining the set of locations or nodes to be visited and the set of directed links connecting them. 𝑁 Set of nodes of graph 𝐺. 𝑜 The origin node of the mission, 𝑜∈ 𝑁. 𝑒 The ending node of the mission, 𝑒∈ 𝑁. 𝐴 Set of directed links in 𝐺 𝛿+(𝑖) Set of nodes that can be reached from node 𝑖∈𝑁 using links in 𝐴. 𝛿−(𝑖) Set of nodes that can be used to reach node 𝑖∈𝑁 using links in 𝐴. Parameters: 𝑄 Maximum level of drone battery charge (just when it is replaced), expressed in time units. 𝑡𝑖𝑗 𝑇 Truck running time at link (𝑖,𝑗)∈ 𝐴. 𝑡𝑖𝑗 𝐷 Drone running time at link (𝑖,𝑗)∈ 𝐴. 𝑀 A big enough constant. Variables: 15 The two above comparisons point out that solving the TDTL problem, either as a permutation problem or by mathematical programming, is much more complex than solving a PFSP or TSP of similar size. Figure 4 shows the ratio between the search space using MIP and the one using the permutation approach for the TDTL, i.e. 𝑆𝑆𝑀𝐼𝑃(𝑛)/𝑆𝑆𝑃𝐸𝑅(𝑛). Figure 4. Ratio between 𝑺𝑺𝑴𝑰𝑷 and 𝑺𝑺𝑷𝑬𝑹 It can be observed that for 𝑛≥10 the search space by the MIP approach exponentially increases with respect to the permutation approach. This suggests not only the burden found for a MIP approach for instances of size greater than 10, but also that a permutation approach exhibits smaller increases on the search space beyond this limit. In practice this means that the greater the instances, the greater the comparative advantage of using a heuristic (combinatorial) method on them. 6. The algorithm 6.1. The iterated greedy approach To solve the problem and given the similarity of the proposed coded solution with permutation problems, we propose an adaptation of an iterative greedy search algorithm (IG) that was used efficiently for task completion time minimization in a permutation flowshop scheduling problem (see Ruiz and Stüzle, 2007). The IG approach was inserted into a SA framework to guide the global search, and specifically adapted to our problem as described in Algorithm 1. The process starts with an initial solution 𝜋 that feeds the block carried out by the SA. Once in the block, the procedure begins with a solution destruction and a solution construction mechanism which gives rise to a new solution 𝜋′′. After this phase, a local search is performed, which works over the two arrays 𝜋𝑠 and 𝜋𝑡 of 𝜋′′ through a local insertion mechanism over 𝜋𝑠 and a search for promising solutions on 𝜋𝑡. In the next, each of the algorithm steps is described in detail. Concerning the SA, as in previous works we select an initial temperature based on the traveling costs (distances or times) between nodes and the number of nodes (see e.g., Ruiz and Stüzle, 2007). The initial temperature 𝑇 used, for a given adjustable parameter 𝛺 is as follows: 𝑇𝑆𝑡𝑎𝑟𝑡=𝛺∑𝑡𝑖𝑗 𝑇 ∀(𝑖,𝑗)𝜖𝐴 +∑𝑡𝑖𝑗 𝐷 ∀(𝑖,𝑗)𝜖𝐴 𝑛 (30) 16 where 𝑡𝑖𝑗 𝑇 is the travel time by truck and 𝑡𝑖𝑗 𝐷 the travel time by drone for every (𝑖,𝑗)∈𝐴. The cooling factor (𝛽) considered to decrease the temperature at each iteration has been set to 0.995. Moreover, in this work 𝛺=1 and 𝑇𝑒𝑛𝑑 =0. procedure IG: 𝑇≔𝑇𝑆𝑡𝑎𝑟𝑡 𝜋≔ Initial_solution repeat (while termination criterion not satisfied): 𝜋′≔ Destruction (𝜋) 𝜋′′≔ Construction (𝜋′) 𝜋′′′≔ Combined_Local_search (𝜋′′) if 𝑓(𝜋′′′)<𝑓(𝜋): 𝜋:=𝜋′′′ If 𝑓(𝜋)<𝑓(𝜋𝑏𝑒𝑠𝑡) and 𝜋 is feasible: 𝜋𝑏𝑒𝑠𝑡:=𝜋 elseif 𝑟𝑎𝑛𝑑𝑜𝑚(0,1)<exp(−(𝑓(𝜋′′′)−𝑓(𝜋)) 𝑇): 𝜋:=𝜋′′′ update 𝑇≔𝛽·𝑇 return 𝜋𝑏 Algorithm 1. IG procedure 6.2. Computing a feasible initial solution Given the absence of direct methods to get a good initial solution to the addressed TDTL problem, we have devised a procedure for generating high quality initial solutions. It is a constructive routine that starting from an initial drone-only serving solution, evolves it by using a two-step heuristic combining the 2-opt and Dijkstra’s algorithms for identifying the rendezvous points that finally have to be part of the truck route. The idea of generating an initial solution from solving a TSP problem on all the nodes has formerly been used in Murray et al. (2015) and Luo et al. (2017). Murray et al. (2015) take as a starting route that obtained from solving a TSP assuming that initially all the nodes are visited by the truck. Later some of the nodes are selected so that they become part of the sub-routes made by the drone. In contrast, Luo et al. (2017) start from the route obtained by solving a TSP supposing all the nodes must be visited by the drone and subsequently some nodes of this route are added to the route followed by the truck. We start from this latter idea and solve the involved TSP to obtain a base route by using a 2- opt heuristic procedure which gives good quality solutions in a short computation time. Inspired by this latter idea and using a 2-opt heuristic procedure we are able to obtain good quality initial TDTL solutions in a short computation time. Despite finding a local optimum may (in the worst case) take a sub-exponential time of 2𝑂(𝑛), in practice the 2-opt heuristic converges faster. If only uncrossing moves are allowed, the 2-opt algorithm will perform at most in 𝑂(𝑛3) moves (Johnson and McGeoch, 1997). In our problem, a starting solution is denoted as a sequence 𝜋𝑇𝑆𝑃 2−𝑜𝑝𝑡 (e.g., 𝜋𝑇𝑆𝑃 2−𝑜𝑝𝑡 =(0,8,4,3,6,1,2,7,5,9) is the sequence obtained after the resolution of the TSP wherein all the nodes are initially visited by drones in the illustrative instance in 17 Figure 3). It is important to note that an optimal or very high-quality solution in this phase does not guarantee a good final solution of the problem under study and, therefore, the use of more complex and time-consuming heuristics is ruled out. In general, the 𝜋𝑇𝑆𝑃 2−𝑜𝑝𝑡 solution will not meet the battery capacity constraints for the drone, and therefore it will not be a feasible solution for the problem under study. It is therefore necessary to establish some synchronization truck-drone nodes to replace (or load) the battery. This new problem can be represented by a graph where each arc connects two possible synchronization nodes. This graph 𝐺′(𝑁,A′) is formed by the set 𝑁 containing all the nodes in the initial problem and a subset of edges, 𝐴′⊆𝐴, which is constructed respecting battery capacity constraints according to the procedure described in Algorithm 2. Let 𝑡𝑖𝑗 𝑇 be the travel time by truck from node 𝑖 to 𝑗∈𝑁, and let 𝑡𝑖𝑗 𝐷𝑟 be the travel time by drone, according to the drone route given by the 𝜋𝑇𝑆𝑃 2−𝑜𝑝𝑡. Then the arc cost, 𝑡𝑖𝑗, represents the more restrictive transportation time from node 𝑖 to 𝑗, i.e., 𝑡𝑖𝑗 = max (𝑡𝑖𝑗 𝑇,𝑡𝑖𝑗 𝐷𝑟). In this way, each feasible arc in 𝐴′ is weighted with the most restrictive travel time by comparing the travel time by truck and the drone route established in the previous phase while fulfilling the battery capacity constraint. procedure Construction of set 𝐴’: for 𝑖 = 0 to 𝑛−1: for j =𝑖+ 1 to 𝑛−1: if 𝑡𝐷𝜋𝑠 0(𝑖),𝜋𝑠 0(𝑗) 𝑟≤𝑄: Append 𝑒=(𝜋𝑠 0(𝑖),𝜋𝑠 0(𝑗)) to 𝐴′ 𝑡𝜋𝑠 0(𝑖),𝜋𝑠 0(𝑗)=max (𝑡𝜋𝑠 0(𝑖),𝜋𝑠 0(𝑗) 𝑇,𝑡𝜋𝑠 0(𝑖),𝜋𝑠 0(𝑗) 𝐷𝑟 ) exit for endif endfor endfor Algorithm 2. Set A’ construction phase Therefore, the problem of completing the construction of a feasible initial solution 𝜋0 becomes in getting the shortest path (𝑚𝑖𝑛 𝑠𝑒) between the origin and end nodes in 𝐺′, which can be efficiently solved in 𝑂(𝑛2) using Dijkstra’s algorithm (Dijkstra, 1959). Figure 5 illustrates the first of the 𝑛 steps necessary to compute all possible paths between every node and its successors using the previous 𝜋𝑇𝑆𝑃 2−𝑜𝑝𝑡 route and taking into account the battery capacity constraints. Figure 5. 1st step in A’ set construction phase 18 For the illustration used so far (shown in Figure 3), the shortest path, and therefore our initial solution 𝜋0, is shown in Figure 6, where every activated link appears labelled with the necessary travel time (either by drone in dotted lines or by truck in solid lines). Observe that when the drone is transported by the truck along links (0,8) and (5,9), the solid line is drawn and the travel time of these links is that of the truck. Figure 6. The initial solution for the previous example Proceeding as explained, we have obtained a feasible initial solution to the problem 𝜋0, with 𝑠𝑒= 70.62 , which is coded as: 𝜋𝑡=[0,0,1,0,0,1,1,1,0,0] and 𝜋𝑠= [0,8,6,4,3,1,2,7,5,9]. Once an initial good feasible solution is obtained, the heuristic must evolve by searching for improved solutions. To this end, each iteration uses a destruction-construction mechanism on the last candidate solution in order to generate a new one. Then, it proceeds with a combined local search to improve the fitness of the new solution. Thus, a new candidate solution emerges from the described two-step process, which is compared with the best stored solution, as reported in Algorithm 1. In the next subsection, we explain destruction/construction and combined local search. 6.3. The destruction and construction phases One of the most important steps in the IG algorithm consists of the destruction and construction mechanisms. Different schemes of solution’s destruction and its subsequent reconstruction are applicable, as for example the one used by Ruiz and Stüzle (2007). In that work a given number of elements in the sequence were randomly extracted obtaining a partial sequence 𝜋𝑅 (destruction phase) – some authors also order it according to a given criterion. Later, the sequence was reconstructed adding each of the extracted elements one by one in the same order than 𝜋𝑅, but following a local search procedure for the insertion. This method has proven not to be competitive enough in the preliminary tests 19 performed for our problem. It is important to highlight that while in flowshop scheduling problems the solutions are just encoded in one array, the TDTL problem requires two different arrays, what makes the direct adaptation of the previous referred method difficult. Besides, it is also important to note that the quality of the results and the computing time are both sensitive to the correct setting of the algorithm parameters, e.g., the number of elements to extract, or the sorting method. As an alternative to the previous method, another mechanism based on reverse swapping has been sought, which has proven to be robust and computationally efficient. The destruction mechanism designed for the TDTL problem is as follows: Two random numbers (𝑖,𝑗) are selected within the sequence 𝜋𝑠 so that 𝑖≠ 𝑗 and (𝑖,𝑗)≠{𝑜,𝑒}. The partial sequence 𝜋𝑅(𝑖,𝑗) is extracted from 𝑖 to 𝑗. The vector of node types 𝜋𝑡 remains the same. As a construction mechanism, the partial subsequence previously obtained 𝜋𝑅(𝑖,𝑗) is inserted in the same position as before but in reverse order. When changing the order of the sequence and not the type of node, it is possible to improve the solution by modifying 𝜋𝑡. In order to improve the construction phase, the process ends with a local search over the node type. The local search procedure which modify the array 𝜋𝑡 is reported in the Algorithm 3: procedure Local_search_in_t: improve := True repeat (while improve = True): improve := False for 𝑖=1 to 𝑛−2: Select 𝑖 at random from 𝑖=1 to 𝑛 −2 (without repetition) 𝜋′≔best solution obtained for every type of node {0,1,2} at position 𝑖 in 𝜋𝑡. if 𝑓(𝜋′)<𝑓(𝜋): 𝜋:=𝜋′ improve := True if 𝑓(𝜋)<𝑓(𝜋𝑏𝑒𝑠𝑡) and 𝜋 is feasible: 𝜋𝑏𝑒𝑠𝑡:=𝜋 endif endfor return 𝜋 Algorithm 3. Local search for ending the construction phase It should be noted that this mechanism does not need the adjustment of parameters, which makes it more robust. Figure 7 illustrates this procedure for a particular example. 20 a) Destruction phase b) Construction phase Figure 7. The destruction-construction procedure 6.4. The combined local search To carry out the local search, the mechanism of local search insertion (see e.g., Ruiz and Stüztle, 2007) is used. This technique has shown better performance than swapping or exchange techniques in flowshop scheduling problems. However, unlike the flowshop scheduling, the problem addressed in this paper also considers, apart from the sequence array 𝜋𝑠, the type of node in a second array 𝜋𝑡, so that for a certain insertion at position 𝑖, it is necessary to select the best value of 𝜋𝑡(𝑖) among the possible values in {0,1,2}. The procedure is complemented running again, for the more promising solutions, the search reported in Algorithm 2 over 𝜋𝑡. This final enhancement has improved the quality of the best solution so far. The combined local search procedure is shown below in Algorithm 4: 21 procedure Combined_local_search: improve := True repeat (while improve = True): improve := False for 𝑖=1 to 𝑛−2: Remove node 𝑘 at random from 𝜋𝑠, 𝑘≠{𝑜,𝑒}, (without repetition) 𝜋′≔ is the solution obtained by the best permutation in 𝜋𝑠 by inserting 𝑘 in all possible positions of 𝜋𝑠 and for every type of node {0,1,2}. if 𝑓(𝜋′)<𝑓(𝜋): 𝜋:=𝜋′ improve := True if 𝑓(𝜋)<𝑓(𝜋𝑏𝑒𝑠𝑡) and 𝜋 is feasible: 𝜋𝑏𝑒𝑠𝑡:=𝜋 𝜋≔ Local_search_in_t (𝜋𝑏𝑒𝑠𝑡) endif endfor return 𝜋 Algorithm 4. Local search for improving the solution The solution obtained by the adapted IG algorithm for the example followed until now is shown in Figure 8. Figure 8. Illustration of the IG final solution in a 10-nodes network 22 7. Computational study In order to test the performance of proposed heuristic, we have conducted an experimental study on a literature-based set of instances with different characteristics. The experiments have been solved on an Intel Core i7 2.8 GHz with 16 GB of RAM. As explained in Section 5.5., it is expected the Branch-and-Cut not being capable of solving medium and large size problems. This circumstance has arisen when performing a set of preliminary experiments on several instances of size greater than 20 nodes. Specifically, we have modelled and solved the TDTL problem by using the Python interface and the standard Branch-and-Cut solver algorithm from Gurobi 7.0.0 along a 24 hours computation time. Table 3 shows the obtained results. Instance size (nodes) Optimality Gap (after 24h of computation time) 20 87.1 % 50 - 75 - 100 - 175 - 250 - Table 3. Results of the Gurobi Solver after 24h of computation time. Secondly, we have addressed the same set of instances from 20 to 250 nodes but this time providing the MIP solver with an initial solution arising from considering that all the nodes are only served using the truck and stopping its execution at 1800s. The results are reported in Table 4, along with the performance provided by our IG procedure when its time limit is set at 60s and the maximum number of iterations without improvement is set at 50. Instance size (nodes) Objective function value at the initial solution (only trucks) Best solver solutions after 1800s Number of improved solutions Objective value of the heuristic after 60s % of difference between the initial solution and the 60s heuristic solution 20 329.617 182.277 2 174.012 47.207 50 564.955 564.955 0 317.814 43.745 75 564.955 564.955 0 363.221 35.707 100 794.478 794.478 0 395.899 50.168 175 1002.578 1002.578 0 506.982 49.432 250 1161.264 1161.264 0 604.862 47.913 Table 4. Results of the Gurobi Solver after feeding the problem with an initial solution. According to the conducted preliminary experiments, we can conclude that the Branch- and-Cut procedure was not expected to be able to solve medium and large problems. 23 Henceforth, we present a comprehensive test bed for testing the performance of our heuristic proposal. Coded in standard Python, it addresses these new set of experiments adjusted to a 600s time limit (IG_TL) and to a maximum number of iterations without improvement equal to 50. For the sake of comparison, we also report the Branch-and-Cut performance on the new set of experiments with time limit truncated to a value MIP_TL of 1800s. 7.1. Instances The set of instances used in our study is based in that previously reported by Agatz et al. (2018) and accessible in Bouman et al. (2018). The authors proposed a different set of randomly generated instances where every node is in a plane. The travel time between nodes is considered proportional and constant between every pair of nodes. Three different sets are considered: uniform, 1-center, and 2-center. In the uniform set, the nodes are located uniformly according an integer uniform distribution (1,100). The 1-center and 2-center set of instances try to mimic the density of customers in a city with one or two circular city centres respectively. Details about the construction of the center scenarios can be found in Agatz et al. (2018). In our study, for all the instances, we have considered a constant battery life (𝑄) proportional to the average drone travel time between each pair of nodes. Therefore 𝑄=2 𝑛∑𝑡𝐷𝑖𝑗∀(𝑖,𝑗)∈𝐴 . Our research is based on the premise of making customer visits between a source and a destination node, the former been coincident with the depot in the set of instances taken from literature. Besides, the last node in every instance specification has been taken as the destination node. For each one of the sets of instances, we apply the heuristic procedure by varying the penalization factor 𝑝=𝑣𝐷/𝑣𝑇 in the objective function for equal speed between truck and drone, twofold and threefold. The bed of instances contains 10 problems for each value of the number of clients 𝑛. We have divided this set in two different groups i) problems of small size, for 𝑛 = {5,6,7,8,9,10} and ii) problems of medium/large size with 𝑛={20,50,75,100,175,250}. In total, 1080 MIP instances have been solved (or partially solved), computing 120 cases (10 for each one of the 12 problem sizes) by varying 𝛼={1,2,3} and considering 3 scenarios (uniform, 1-center, and 2-center). Concerning to the heuristic, since by its nature it is a non-deterministic procedure, 10 replications have been made for each one of the 1,080 instances, solving a total amount of 10,800 cases. Figures 9, 10 and 11 depict three solutions for uniform instances of 75, 100 and 250 nodes respectively. 24 Figure 9. Uniform 75 nodes. Q= 54.63, Obj. function=402.21. Comp. time=163.28s Figure 10. Uniform 100 nodes. Q= 53.24, Obj. function =431.10. Comp. time=384.81s Figure 11. Uniform 250 nodes. Q=51.56, Obj. function =584.12, Comptime=605.19s 7.2. Performance measures As in many other challenging problems addressed in the scientific literature, the basic quality measure of the results is the relative percentage deviation (𝑅𝑃𝐷) over the bestknown solution. For a given solution 𝜋 and a best-known solution 𝜋𝑏 with objective values 𝑠𝑒(𝜋) and 𝑠𝑒(𝜋𝑏) respectively, the 𝑅𝑃𝐷 is obtained as follows: 𝑅𝑃𝐷=𝑠𝑒(𝜋)−𝑠𝑒(𝜋𝑏) 𝑠𝑒(𝜋𝑏)×100 (31) 31 The same conclusions can be obtained when performing the Kruskal-Wallis test with respect to the analysis of the different types of scenarios under study. See Table 9. Scenario 𝑵 Avg Range Uniform 18 29.58 1-center 18 30.86 2-center 18 22.06 Total 54 Table 9. Average Range, Uniform, 1-center, and 2-center scenarios In this case, the Kruskal-Wallis test statistic (𝐻), is calculated as 3,313. At 0.05 level of significance, the critical value obtained from the chi-square table with 2 degrees get a 𝑝value of 0.191 (>𝑙𝑒𝑣𝑒𝑙 𝑜𝑓 𝑠𝑖𝑔𝑛𝑖𝑓𝑖𝑐𝑎𝑛𝑐𝑒), and therefore we cannot reject the null Hypothesis of equal Avg. RPD(%) medians among groups. As conclusion, in both cases, at 0.05 level of significance, there is insufficient evidence to conclude that at least one of the groups differs with respect to the median Avg. RPD(%). The Figure 13 shows the results of the medians (marked with different symbols) for the average RPD (%) values grouped by 𝛼, Figure 13(a), or by type of scenario, Figure 13(b). In addition to the median, the ranges of variation of the data (minimum and maximum values) are also depicted in the graphs. (a) RPD(%) vs 𝛼 (b) RPD (%) vs Scenario Figure 13. Average RPD(%) values The previous analysis shows that the average RPD(%) obtained by the heuristic in the small instances, of which the deviation from the optimum is known in most cases, does not vary significantly according to the Kruskal-Wallis analysis. In addition, values are, in all cases, inferior to 4.1%, and the median values are less than 2%, and very close to 0. In 32 conclusion, we can say that the analysis reinforces the characteristics of robustness and quality of the proposed heuristic. 8. Conclusions In this paper, we propose a mathematical formulation for the TDTL problem, determining the best routes in the use of a single truck equipped with a single drone, with both cooperating to serve a set of locations. The problem lies in determining the best routes for both, truck and drone, between two different nodes (origin and destination), without predefining the waypoint or rendezvous locations for either the truck or the drone. Our formulation is indeed an extension of the FSTSP model, allowing the drone to visit several customers per trip between two consecutive rendezvous with the truck. We go beyond the general assumption in hybrid truck-drone systems concerning that the locations are served one at a time. Although the above general assumption can be acceptable in the last-mile delivery case, it makes no sense in other scenarios such as the delivery of lightweight relief items in first aid scenarios or in ISR missions. Aside from the multi-drop assumption other features of our novel model is the absence of a preestablished route (open) for both truck and drone and the consideration of each location as a potential synchronization point between them. We suppose infinite travel length capacity for the truck and a limited battery life for the drone, but rechargeable at certain (a priori undetermined) rendezvous points. Since the MIP model stated becomes intractable for medium or large-size instances, we have also provided an approximate IG heuristic method to solve it in an efficient way. To this aim, we have defined an original coding scheme of solutions based on a tuple formed by two vectors, a sequence vector and a vector containing the type of resource (truck and/or drone) used to visit every node. We have taken advantage of this coding scheme, firstly in developing a fast-constructive two-step procedure for generating high-quality initial solutions; and secondly, in improving the solution through a destructionconstruction mechanism followed by two original local search procedures. This process is conducted by a global optimization scheme using a SA algorithm. We have checked the developed heuristic, addressing a set of 1080 problems from literature-based scenarios. In this benchmark set of instances, we have compared the heuristic’s performance with the exact resolution using the standard Branch-and-Cut algorithm of the Gurobi solver. As the extensive computational study reveals that the heuristic shows a good performance in both its ability to solve all the instances in a short computation time, and the provision of high-quality solutions. Since our proposal is able to produce good solutions for large-sized problems, we can state it is ready for application in solving current real-life truck-drone routing problems. References Agatz, N., Bouman, P., Schmidt, M. (2018). Optimization approaches for the traveling salesman problem with drone. Transportation Science, 52(4), 965–981. 33 Othman, M. S., Shurbevski, A., Karuno, Y., Nagamochi, H. (2017). Routing of Carrier-vehicle Systems with Dedicated Last-stretch Delivery Vehicle and Fixed Carrier Route. Journal of Information Processing, 25, 655–666. Bouman, P., Agatz, N., Schmidt, M. (2017). Dynamic programming approaches for the traveling salesman problem with drone. ERIM Report Series Research in Management (ERS-2017-011– LIS). Bouman, P., Agatz, N., Schmidt M. (2018). Instances for the TSP with drone (and some solutions). https://doi.org/10.5281/zenodo.1204676. Boysen, N., Schwerdfeger, S., Weidinger, F. (2018). Scheduling last-mile deliveries with truckbased autonomous robots. European Journal of Operational Research, 271(3), 1085–1099. Campbell, J. F., Sweeney, D., Zang, J. (2017). Strategic design for delivery with trucks and drones. Technical Report. SCMA-2017-0201. College of Business Administration. University of Missouri. St. Louis. MO USA. Chang, Y. S., Lee, H. J. (2018). Optimal delivery routing with wider drone-delivery areas along a shorter truck-route. Expert Systems with Applications. 104. 307–317. Cheng, C., Adulyasak, Y., Rousseau. L. M. (2018). Formulations and Exact Algorithms for Drone Routing Problem. Technical report. 2018-2018-31. CIRRELT. Canada. Chowdhury, S., Emelogu, A., Marufuzzaman, M., Nurre, S. G. (2017). Drones for disaster response and relief operations: A continuous approximation model. International Journal of Production Economics, 188, 167–184. Dayarian, I., Savelsbergh, M., Clarke, J. P. (2017). Same-day delivery with drone resupply, Optimization Online, http://www.optimization-online.org/DB_FILE/2017/09/6206.pdf. Dijkstra, E. W., (1959). A note on two problems in connection with graphs, Numer. Math., 1, 83– 89. Dorling, K., Heinrichs, J., Messier, G. G., Magierowski, S. (2017). Vehicle Routing Problems for Drone Delivery, IEEE Transactions on Systems, Man. and Cybernetics, 47(1), 70–85. Ferrandez, S. M., Harbison, T., Wieber, T., Stnuges, R., Rich, R. (2016). Optimization of a truckdrone in tandem delivery network using k-means and genetic algorithm, Journal of Industrial Engineering and Management, 9(2), 374–388. Garone, E., Naldi, R., Casavola, A. (2011). Traveling salesman problem for a class of carriervehicle systems, Journal of Guidance, Control, and Dynamics, 34(4), 1272–1276. Greenwood, F. (2015). Above and Beyond: Humanitarian Uses of Drones, World Politics Review, Available from: http://www.worldpoliticsreview.com/articles/16750/above-and-beyond- humanitarian-uses-of-drones. [Accessed on December 15. 2018]. Ha, Q. M., Deville, Y., Pham, Q. D., Hà, M. H. (2018). On the min-cost Traveling Salesman Problem with Drone, Transportation Research Part C: Emerging Technologies, 86, 597–621. Ham, A. M., 2018. Integrated scheduling of m-truck, m-drone, and m-depot constrained by timewindow, drop-pickup, and m-visit using constraint programming, Transportation Research Part C: Emerging Technologies, 91, 1−14. Jeong, H. Y., Lee, S., Song, B. D. (2019). Truck-Drone Hybrid Delivery Routing: Payload-Energy dependency and No-Fly Zones, International Journal of Production Economics (In Press). Johnson, D. S., McGeoch, L. A. (1997). The Traveling Salesman Problem: A Case Study in Local Optimization. In: Aarts, E. H. L. and Lenstra, J. K.. Eds., Local Search in Combinatorial Optimization, John Wiley & Sons, New York, 215−310. Karak, A., Abdelghany, K. (2019). The hybrid vehicle-drone routing problem for pick-up and delivery services, Transportation Research Part C: Emerging Technologies, 102, 427−449. 34 Kitjacharoenchai, P., Ventresca, M., Moshref-Javadi, M., Lee, S., Tanchoco, J. M. A., Brunese, P. A. (2019). Multiple traveling salesman problem with drones: Mathematical model and heuristic approach, Computers & Industrial Engineering, 129, 14−30. Luo, Z., Liu, Z., Shi, J. (2017). A two-echelon cooperated routing problem for a ground vehicle and its carried unmanned aerial vehicle, Sensors, 17, 1144. Manyam, S. G., Sundar, K., Casbeer, D. W. (2018). Cooperative Routing for an Air-Ground Vehicle Team--Exact Algorithm, Transformation Method, and Heuristics, arXiv preprint arXiv:1804.09546. Marinelli, M., Caggiani, L., Ottomauelli, M., Dell`Orco, M. (2017). En-route truck-drone parcel delivery for optimal vehicle routing strategies, IEE Intelligent Transport Systems, l2(4), 253–261. Mathew, N., Smith, S. L., Waslander, S. L. (2015). Planning Paths for Package Delivery in Heterogeneous Multirobot Teams, IEEE Transactions on Automation Science and Engineering, 12, 1298–1308. Morgan, P. (2005). Carbon Fibers and Their Composites, London: Taylor and Francis. Moshref-Javadi, M., Lee, S. (2017). Using drones to latency in distribution systems, IIE Annual Conference, Proceedings, 235−240, Institute of Industrial and Systems Engineers (IISE). Murray, C. C., Chu, A. G. (2015). The flying sidekick traveling salesman problem: Optimization of drone-assisted parcel delivery, Transportation Research Part C: Emerging Technologies, 54, 86–109. Murray, C., Raj, R. (2020). The multiple flying sidekicks traveling salesman problem: Parcel Delivery with Multiple Drones, Transportation Research Part C: Emerging Technologies, 110, 368–398. Poikonen, S., Wang, X., Golden, B. (2017). The vehicle routing problem with drones: Extended models and connections, Networks, 70(1), 34–13. Ponza, A. (2016). Optimization of drone-assisted parcel delivery. Master Thesis. University of Padova. Italy. Pugliese, L. D. P., Guerriero, F. (2017). Last-mile delivery by using drones and classical vehicles, International Conference on Optimization and Decision Science, 557–565 (Springer). Rabta, B., Wankmüller, C., Reiner, G. (2018). A drone fleet model for last-mile distribution in disaster relief operations, International Journal of Disaster Risk Reduction, 28, 107–112. Reddy, T. B. (2010). Chapter 14: Lithium primary batteries. In: Linden’s Handbook of Batteries, 4th edition, 14.1–14.90 (McGraw-Hill Professional. New York). Ruiz, R., Stützle, T. (2007). A simple and effective iterated greedy algorithm for the permutation flowshop scheduling problem, European Journal of Operational Research, 177, 2033–2049. Sacramento, D., Pisinger, D., Ropke, S. (2019). An adaptive large neighborhood search metaheuristic for the vehicle routing problem with drones, Transportation Research Part C: Emerging Technologies, 102, 289–315. Schermer, D., Moeini, M., Wendt, O. (2019). A matheuristic for the vehicle routing problem with drones and its variants, Transportation Research Part C: Emerging Technologies, 106, 166–204. Ulmer, M. W., Thomas, B. W. (2018). Same‐day delivery with heterogeneous fleets of drones and vehicles, Networks, 72, 475–505. Wang, X., Poikonen, S., Golden, B. (2017). The vehicle routing problem with drones: Several worst-case results, Optimization Letters, 11(4), 679–697. Wang, Z., Sheu, J. B. (2019). Vehicle routing problem with drones, Transportation Research Part B: Methodological, 122, 350–364. 35 Wen, T., Zhang, Z., Wong, K. K. L. (2016). Multi-Objective Algorithm for Blood Supply via Unmanned Aerial Vehicles to the Wounded in an Emergency Situation, PloS one, 11(5), 1–22. Yurek, E. E., Ozmutlu, H. C. (2018). A decomposition-based iterative optimization algorithm for traveling salesman problem with drone, Transportation Research Part C: Emerging Technologies, 91, 249–262.