key: cord-0102283-vmfj2tfx
authors: Zhu, Tengkuo; Boyles, Stephen D.; Unnikrishnan, Avinash
title: Electric Vehicle Traveling Salesman Problem with Drone with Fixed-time-full-charge Policy
date: 2022-01-20
journal: nan
DOI: nan
sha: 41b69f642ee5bec872dd7578f7bd781cedc5ec33
doc_id: 102283
cord_uid: vmfj2tfx

The idea of deploying electric vehicles and unmanned aerial vehicles (UAVs), also known as drones, to perform"last-mile"delivery in logistics operations has attracted increasing attention in the past few years. In this paper, we propose the electric vehicle traveling salesman problem with drone (EVTSPD), in which the electric vehicle (EV) and the drone perform delivery tasks coordinately while the electric vehicle may need to visit charging stations occasionally to recharge. We further assume that the EV can refresh its energy to full battery capacity with fixed time at charging stations. Thus, the proposed problem is termed EVTSPD-FF. In this paper, an arc-based mixed-integer programming model defined in a multigraph is presented for EVTSPD-FF. An exact branch-and-price (BP) algorithm and a variable neighborhood search heuristic are developed to solve instances with up to 25 customers in one minute. Numerical experiments show that the heuristic is much more efficient than solving the arc-based model using the ILOG CPLEX solver and BP algorithm. A real-world case study on the Austin network and the sensitivity analysis of different parameters are also conducted and presented. The results indicate that drone speed has a more significant effect on delivery time than the EV's driving range.

Although e-commerce was growing fast before COVID-19 hit, the ongoing pandemic has accelerated the shift towards a more digital world and pushed more consumers worldwide online. Based on the United Nations conference on trade and development, the e-commerce sector saw a significant rise from 16 percent to 19 percent in 2020 in its share of all retail sales. According to the US commerce department, in the United States, online sales hit $791.70 billion in 2020, up 32.4% from $598.02 billion in the prior year. With consumers increasingly turning to e-commerce for their shopping needs, efficient order fulfillment and parcel distribution is the expectation of every online shopping experience. As a result, businesses in this area have begun racing to develop new technologies and experimental supply chain models to increase parcel volume and expedite deliveries.

One of the most significant expenses and challenges in the item shipping process is same-day, last-mile delivery, where the goods are transported from a distribution hub to the customer's front door. The last mile of delivery cost accounts for about half of the total shipping costs, as the final leg of shipment typically involves multiple stops with low drop sizes, especially in urban areas. In recent years, a large variety of alternative last-mile transportation modes appeared in both academia and industry, such as electric vehicles (Lin et al., 2016) , autonomous robots (Jennings and Figliozzi, 2019) and unmanned aerial vehicles, or drones . Some new delivery concepts are also proposed, such as crowdsourcing delivery (Rougès and Montreuil, 2014) . The problem researched in this paper involves the operation of two new techniques mentioned above: electric vehicles and drones.

Meanwhile, political and practical considerations in recent years, such as concern about greenhouse gas emissions, are driving the move away from internal combustion engines towards trucks and vans using alternative power sources. The deployment of a new generation of electric vehicles is critical for reducing greenhouse gas emissions. This concept has also been extended to logistics. The new invention of electric delivery trucks or vans is increasingly replacing the traditional delivery truck as an approach for logistics companies to "boost their bottom lines" while fighting climate change and urban pollution. UPS, Amazon, and DHL all announced their initiatives to partially or fully replace their current pickup and delivery fleet with battery-powered vehicles.

At the same time, auto manufacturers compete with start-ups to produce the most efficient and "smart" electric delivery vehicles. One of the latest inventions is small electric delivery vans, a picture of which is shown in Figure 1 . It is an electric, lightweight van designed especially for small or medium-sized package delivery. Compared to the traditional bulky delivery truck, this small van is significantly more eco-friendly and cost-effective. This van is compact in dimensions (140-200 cm of width), can carry a load up to 1600-2000 kilograms, and has a maximum driving range of 60-200 km with a lithium battery. This small-sized and flexible electric vehicle might represent the future landscape of logistics.

Delivery via unmanned aerial vehicles, or drones, has also shown to have great potential to satisfy future customer expectations for delivery service quality as efficient, environmentally friendly, low maintenance cost delivery modes. The use of drones as a delivery method has been applied in several areas such as logistics, military operations, public security, traffic surveillance, monitoring, and humanitarian relief. In particular, the use of drones in logistics has sparked continuous interest for the general public ever since online retail giant Amazon announced its "Prime Air" initiative back in 2013. Multiple companies, such as Google, DHL, UPS, and Walmart, have revealed and tested drone delivery services in various countries. Compared to ground transportation delivery, drone delivery has the clear advantage of greater traveling speed, as it can travel directly from the launch point to the customer's location without being affected by the congested ground traffic. This property is crucial for several time-sensitive applications such as blood and medicine delivery. Besides, based on the results reported in (Goodchild and Toy, 2018) , drones tend to have a 2 emissions advantage over trucks in service zones that are either closer to the depot or have smaller numbers of customers. The past few years have witnessed a dramatic increase in UAV applications (DroneZon, 2019) . According to Teal Group's prediction, commercial use of UAVs will grow eightfold over the decade to reach US $7.3 billion in 2027 (Teal Group, 2018) . For the rest of the paper, the term "UAV" and "drone" are used interchangeably.

However, these new techniques do not come without drawbacks. For electric vehicles, although their driving range has been dramatically improved recently due to the advancement in the power grid system and lithium battery, this value ranges typically from 150 miles to 300 miles. The driving range is significantly smaller than electric passenger cars for normal-sized electric trucks or vans for commercial use. For example, the expected driving range of electric delivery vans built by Rivian, an American electric vehicle automaker, is about 150 miles (Insideevs, 2021) . The Mercedes-Benz eSprinter for Amazon has a maximum range of 104 miles with a 55 kWh battery and 70 miles with a 41 kWh battery (?). As for small-sized electric vans, their driving range is about 60 miles to 190 miles with a 10-20 kWh lithium battery. Thus, for electric vans that can only cover less than 100 miles per charge, their success and fluent operation require them to recharge their battery en route when necessary. For drone delivery, apart from the limitation of related regulations of the Federal Aviation Administration and privacy issues in the urban area, its drawback mainly lies in a limited capacity and flight range. A rotary-winged drone, such as the one used in project "Amazon prime air," typically has a capacity of up to five pounds, with a maximum flight range of 15 miles. These values also indicate that a launched drone can only serve one customer before it is retrieved and recharged in a practical scenario.

Several new delivery concepts have been proposed in academia to mitigate the disadvantages of drone delivery, one of which is the in-tandem delivery method where the trucks are teamed up with drones to accomplish delivery tasks. Initially presented in (Murray and Chu, 2015) , this idea requires that in operation, the truck carries a large number of parcels and charging facilities for the drones and serves as a drone hub while delivering parcels independently. The truck could launch the drone and retrieve the drone later before the drone's energy is depleted. This setting leads to a combinatorial optimization problem that requires the coordination of both truck and drone routes.

This paper describes a similar problem while replacing the traditional delivery truck with a small electric truck/van. This problem is called the electric vehicle traveling salesmen problem with drone (EVTSPD). A single truck/van equipped with a drone performs the delivery task in this problem. On its route, it needs to visit charging stations in the network to recharge its battery energy when necessary. Besides, we further assume that the truck/van shares its energy with its equipped drone. When the drone is launched from the truck/van to perform a delivery task, the energy required for drone delivery should be deducted from the remaining energy of the truck. It turns out that the addition of this assumption significantly increases the difficulty of the problem and needs specially constructed techniques to handle it. From a practical point of view, by combining two newly emerged techniques -EV and drone -we can gain some insight into the future landscape of logistics by extending current solution techniques to more complicated problem variants we may face later.

In general, the main contributions of this paper are:

• The paper investigates the combination of EV and drone routing and propose the electric vehicle traveling salesman problem with drone. In EVTSPD, the proposed shared energy assumption ensures that the coordination of two vehicles should be considered throughout the operation. • In terms of the charging policy, this paper assumes that the EV could refresh its battery to full capacity with fixed amount of time every time it visits the charging station nodes. This proposed problem, EVTSPD with fixed-time-full-charge policy, is termed as EVTSPD-FF. • A mixed integer/linear programming formulation of EVTSPD-FF is presented. This formulation is defined in a constructed multigraph, which contains no charging station nodes at the expense of an increased number of arcs. It is estimated that when three charging stations exist in the instance, the constructed multigraph has approximately 3-5 times more arcs compared to the original network. • We present a set partitioning formulation for EVTSPD-FF and propose an exact branch-and-price (BP) solution method based on the multigraph. A dynamic programming (DP) approach is used to derive the least-cost route in the pricing problem of the branch-and-price algorithm, based on the ng-route relaxation of the problem. This approach proved to be efficient in identifying the columns which would enter the master problem of the BP algorithm. • The numerical analysis results indicate that the BP method can solve instances containing up to 10 customers to proven optimality in one hour, which is significantly more efficient than solving the problem via a commercial solver. • A variable neighborhood search (VNS) heuristic method is also proposed in this paper to solve EVTSPD-FF of practical size. • Extensive computational runs are conducted to test the efficiency of the proposed MILP model, BP algorithm, and VNS heuristic. A case study based on a real-world Austin network is also performed, along with the sensitivity analysis of several critical parameters of EV and UAV. The rest of the paper is organized as follows. Section 2 provides a thorough literature review of the related and recent research on EV routing and drone routing problems. Section 3 describes the EVTSPD-FF and its MILP formulation. Section 4 proposes the set partitioning formulation and the exact BP algorithm used to solve the problem. Specifically, it introduces the DP method used to solve the pricing problem of BP. Besides, a VNS method is also proposed in this section. Section 5 presents the computational results of the BP algorithm compared to MILP models, the performance evaluation of VNS, and a real-world case study. Finally, section 6 concludes the work and gives insight into potential future research directions.

The number of publications related to EVs or drones has grown significantly in recent years. A thorough and comprehensive on both areas are not the main subject of this section. Thus we only include recent papers that focus specifically on EV and drone routing.

The electric vehicle routing problem is a recent variant of traditional VRP. More details about VRP and its variants could be seen in (Braekers et al., 2016 , Montoya-Torres et al., 2015 . EVRP is a branch of the alternative-fuel powered vehicle routing problem where the conventional vehicles are replaced by alternative-fuel powered vehicles (AFVs), which have limited fuel tank capacity or battery. Thus, EVRP is also a branch of the green vehicle routing problem (GVRP). In fact, EVRP is called GVRP in the early papers of this topic, such as (Erdoğan and Miller-Hooks, 2012) . Readers can refer to (Asghari et al., 2020) as the latest literature review on GVRP. For a review paper that specially focuses on the electric vehicle routing problem and its variants, readers can refer to (Erdelić and Carić, 2019) .

There are two basic configurations of EVs: the battery electric vehicle (BEV) and the hybrid electric vehicle (HEV). The former one is powered exclusively by batteries inside the vehicle, while the latter one can be powered by batteries and an internal combustion engine. Consequently, BEV has much less driving range compared to HEV. Similar to a majority of EV routing-related papers regarding BEV as EV, in this problem, we also assume that the electric van is BEV. For the rest of the paper, we assume that EV is referred to as BEV if not clearly stated otherwise. An EV could renew its battery energy at charging stations in two ways: by swapping empty batteries with fully charged ones or by plug-in charging via an electric power source. The battery swapping process could be finished within minutes, while the plug-in charging time depends on the desired state of charge level and charging function.

The early research on the routing of electric fleet includes (Gonçalves et al., 2011, Conrad and Figliozzi, 2011) . The former focuses on VRPPD with mixed feet of EVs and conventional internal combustion engine vehicles where the charging time is based on the total distance traveled, while the latter one assumes that vehicles are allowed to recharge at customer locations to recover up to 80% of its charge capacity. The numerical analysis results of the latter paper imply that customer time windows have a greater impact on the tour distance recharging time is long. (Erdoğan and Miller-Hooks, 2012) formally formulated GVRP in which a fleet of AFVs are involved. The AFVs can refresh their driving electricity at charging station nodes with fixed charging time. No customer time window and vehicle capacity are considered in the problem. The MILP formulation is proposed based on the augmented network. Two problem-specific heuristics are presented: one is based on the Clarke and Wright Savings algorithm, while the other one is based on the clustering algorithm. (Schneider et al., 2014) proposed a problem that considers both the service time window and vehicle capacity. This problem assumes a constant and linear charging rate at the charging station node, while the recharging time depends on the remaining energy level of the vehicle when arriving at charging station nodes. A hybrid variable neighborhood search heuristic/tabu search method is used to solve the problem. In the previously discussed paper, the charging station nodes are copied multiple times to create the augmented network to enforce that each node could visit at most one time. A different approach is adopted in (Bruglieri et al., 2016) where for each pair of customers, the charging stations are pre-computed to ensure that the vehicle can travel through the two customers via the computed charging stations. A new MILP model is proposed based on this network. However, the limitation of this approach is that only one charging station node can be traversed between two consecutive customer visits, which could be sub-optimal on some occasions.

All the above-mentioned papers use the battery swapping scheme and assume fixed charging time, which is not accurate in practice. In fact, the recharging of the battery employs a multi-phase charging scheme (Gao et al., 2002, Tremblay and Dessaint, 2009) to protect the battery from over-charging, which makes the function of charging time nonlinear. (Montoya et al., 2017) extended the previous E-VRP models to consider nonlinear charging functions and used a piece-wise linear approximation to simplify the process. The numerical analysis in the paper indicates that neglecting nonlinear charging might lead to an infeasible route. A local search meta-heuristic is used to solve the problem. (Froger et al., 2019) proposed two new formulations for nonlinear charging EVRP: one is the arc-based formulation, and another one is path-based formulation. The numerical analysis results show that the path-based model outperforms the one in (Montoya et al., 2017) . (Koç et al., 2019) extended the problem by considering shared charging stations. This new problem aims to minimize the total charging station construction cost and driving cost. A similar problem is proposed in (Kullman et al., 2021) , in which the author considers EVRP with public charging stations. In the public charging station, if the demand is full, the waiting time is modeled with queuing theory. A decomposition approach is used to derive the bounds of the problem.

With the complex constraints of EVRP, most of the papers use heuristics/meta-heuristics to obtain a solution. Exact methods for the EVRP and its variants are rather scarce. (Koç and Karaoglan, 2016) presented a new formulation while assuming that no more than one refueling stop is made by any vehicle when traveling between two customers or the depot and one customer, which is similar to the approach used in (Bruglieri et al., 2016) . A branch-and-cut algorithm is used to obtain the lower bounds, while a heuristic based on simulated annealing is adopted to obtain upper bounds. This solution method can solve instances containing up to 20 customers. (Desaulniers et al., 2016) adopts an exact branch-price-and-cut algorithm to solve four variants of EVRP. This algorithm uses customized monodirectional and bidirectional labeling algorithms for generating feasible vehicle routes and can solve instances with 100 customers and 21 recharging stations. (Hiermann et al., 2016) also reports an exact branch-and-price algorithm and adaptive large neighborhood search method for EVRP with time window and various fleet sizes. More recently, (Andelmin and Bartolini, 2017) proposed an exact branch-price-and-cut algorithm for solving EVRP with maximum duration constraints based on the creation of a multigraph. This approach avoids the limitation of the approach used in (Koç and Karaoglan, 2016) and is shown to have excellent performance, which can solve instances with 110 customers to optimality. This multigraph approach is also used in this paper. (Lee, 2021) presents an exact branch-and-price approach for the EVRP with nonlinear charging functions by introducing an extended charging stations network. This extended charging stations network is very similar to the refueling path network in this paper. However, in its branch-and-price scheme, the subproblem only obtains route segment and not complete routes, which is different from our approach.

The number of recent publications related to the application of drones is growing extremely fast. (Otto et al., 2018) review several important and potential civil drone applications in the area of agriculture, monitoring, transport, security, etc. This section does not include a detailed review of the drone-related application, so we only consider the papers that involve drone routing-related optimization problems. A recent review on this subject can be seen in (Macrina et al., 2020) .

Several works emerged in recent years that focus on serving customers with a fleet that is composed of drones only. Most of these drone-only routing problems pay special attention to the energy consumption, limited capacity, and flight range of drones. Unlike the truck-drone in-tandem routing problem where a single drone is assumed to serve only one customer per trip, in several drone-only problems, it is normally assumed to be capable of serving multiple customers in one trip. (Dorling et al., 2016) proposes two models of drone-only VRP with different objectives. The first model minimizes costs under the delivery time limit, while the second one minimizes the delivery time under budget constraints. A simulated annealing (SA) heuristic is used to solve the problems. (Coelho et al., 2017) introduced a new problem named the green UAV routing problem, where multiple charging station nodes are introduced. This model also has multiple objective functions, which include total travel distance/time, the number of drones used, and so on. A heuristic called the multi-objective smart pool search (MOSPOOLS) is used to obtain the solution. (Liu, 2019) present a model where drones are used to perform on-demand meal delivery services. It is a pick-up and delivery problem where the meal orders have to be fulfilled in a dynamic order. A time-discretized MIP formulation and a heuristic are proposed in the paper. Additional researches related to this topic includes (Troudi and Osman, 2019) and (Choudhury et al., 2021) .

Another class of drone routing problem is the one that focuses on truck-drone in-tandem delivery systems. One obvious advantage of this delivery method is the drone can greatly extend their service radius because now the truck could be served as a drone hub that launch/retrieve the drone. Besides, it also can be seen as a natural transitional method from the traditional truck-only delivery to fully automated drone fleet delivery. This newly emerged delivery method is firstly proposed in (Murray and Chu, 2015) , where it introduces the flying sidekick TSP (FSTSP), where a single truck and a single drone aim to serve a set of customers. The truck and drone can perform delivery tasks independently, while the drone could be launched from the truck and retrieved later at a different location. The drone can only carry one parcel per trip, and constant launching/retrieving time is considered. The problem aims to minimize the coordinated route duration, including the waiting time caused by the synchronization of both vehicles. A MILP formulation, along with two heuristics, is also proposed in the paper. Since then, various research has been conducted, either trying to improve the solving efficiency of FSTSP or to explore new ways of combining multiple trucks and drones. (Ferrandez et al., 2016 , Es Yurek and Ozmutlu, 2018 , De Freitas and Penna, 2018 , Ha et al., 2018 are some of the examples that aim to study the variant of FSTSP or to propose new solution methods. In particular, (Ferrandez et al., 2016) study the effectiveness of this new delivery system and adopts K-means algorithms to determine optimal launch locations and a genetic algorithm to obtain truck route between launch locations. It finds that the new system is more efficient in terms of energy and time. Similar analyses are also performed in (Carlsson and Song, 2018) where it assumes that drones can be launched/retrieved at some points which are not necessarily customer locations. ) study a variant of FSTSP where the drone can be launched/retrieved at the same customer location. This variant is named the traveling salesman problem with drones (TSPD). Two route-first, cluster-second heuristics are proposed to solve TSPD. present an exact dynamic programming algorithm to solve the same problem proposed in . (Es Yurek and Ozmutlu, 2018) presents an iterative method to solve FSTSP by pre-deciding the customers served by the truck or the drone, where for each case, the truck route and drone route are obtained in sequence. This method performs well when the number of customers is small but unable to solve problems containing more than ten customers. (De Freitas and Penna, 2018) propose a variable neighborhood search method to solve both FSTSP and TSPD. (Ha et al., 2018 ) study a variant of FSTSP where the objective function aims at minimizing total operational costs, which includes transportation cost and waiting cost. A heuristic method based on a greedy randomized adaptive search procedure (GRASP) is used to solve the problem. More recently, (Salama and Srinivas, 2020 ) studied a variant of FSTSP where multiple drones operate in conjunction with a single truck. In this problem, the truck needs to launch several drones at a focal point in each customer cluster. A three-step heuristic is presented, where the first step is to find the focal point and then assign customers to focal points, and lastly, the routing of a truck among these focal points. Besides, (Raj and Murray, 2020) also studies FSTSP with multiple drones and considers different travel speeds, payload capacities, service times, and flight range. A MILP formulation is proposed, and a heuristic approach consisting of three subproblems is proposed.

As a generalization of FSTSP and TSPD, the vehicle routing problem with drones (VRPD) is another attempt to combine trucks and drones in a more complex way. (Wang et al., 2017) introduced VRPD, where a fleet of trucks, each equipped with multiple drones, aims to perform the delivery task. In this problem, drones can be launched/retrieved at any of the customer locations. An analysis of several worst-case scenarios is conducted, from which the best possible time savings can be derived. (Schermer et al., 2018) propose two route-first cluster-second heuristics to solve VRPD, where the improving phase is based on local search. (Wang and Sheu, 2019) extend the problem by proposing an arc-based MILP formulation. In this problem, there exist several drone hub nodes that serve as a hub that can refresh/reload drones. A branch-and-price algorithm is used to solve this problem. (Schermer et al., 2019) considers an extension of VRPD, called vehicle routing problem with drones and en route operations (VRPDERO). In VRPDERO, drone operations might also start and end at some discrete points on the arcs traversed by the vehicles, and the objective is to minimize the makespan. A MILP formulation is proposed, and a variable neighborhood search is used to obtain a solution to the problem. (Sacramento et al., 2019) presents an adaptive large neighborhood search method to solve a VRPD variant where each truck is equipped with exactly one drone. (Kitjacharoenchai et al., 2020 ) considers a two-echelon VRPD where truck routing and drone routing are separate. A similar problem is also seen in (Liu et al., 2020) .

In this paper, we propose and analyze the the routing problem that involves both the EV and the drone, which is named as EVTSPD. In this problem, the shared energy assumption indicates that the synchronization between the EV and the drone has to be considered. Specifically, we focus on EVTSPD-FF which further assumes that the EV could refresh its battery energy with fixed amount of time every time it visits the charging station nodes. This assumption allows us to redefine the problem in a constructed multigraph. A MILP formulation is proposed based on the idea of multigraph, along with an exact branch-and-price scheme and an efficient variable neighbourhood search algorithm to solve the problem of different size.

In EVTSPD, a central planner aims to serve all the customers in the network using a single EV equipped with a single drone, which is initially parked at the depot. In the operation process, the EV and the drone perform delivery tasks in tandem, where the EV can serve as the drone hub that can launch, retrieve and reload the drone, and both the EV and the drone can serve customers independently. This operation pattern is similar to the flying sidekick traveling salesman problem introduced in (Murray and Chu, 2015) and the traveling salesman problem with drone studied in . However, the most significant difference between EVTSPD and the previous TSPD research is that the EV is associated with driving energy which decreases gradually as it traverses arcs in the network. Thus, in EVTSPD, the EV needs to visit specialized charging station nodes in the network to recharge the energy. The objective of this problem is to obtain an EV-drone coordinated route that has minimal operational time.

In the modeling process, the following assumptions or simplifications are adopted, most of which are in analogy to existing VRPs and EVRPs:

1) The drone is assumed to be a rotary-wing UAV so that it can take-off and land vertically on the EV.

2) Both the EV's and the drone's speed are assumed to be constant, and the altitude differences in the network are ignored.

3) The consumed energy of traveling is assumed to be a linear function of traveled distance. 4) At the beginning of the planning horizon, both the EV and the drone are located at the depot and the EV is assumed to be fully-charged. 5) The EV can recharge its energy at the charging stations in the network. The charging process is assumed to start immediately as the EV arrives at the station. It is also assumed that all charging stations have unlimited capacities. 6) All charging stations are considered to be level-3 fast charging stations so that it would not lead to unacceptable high charging times. The EVTSPD is defined on an undirected, complete graph = ( , ), with a node set consisting of the customer set = {1, 2, ..., }, the depot set {0, 0 }, and a set of charging stations (CS) nodes = { + 1, + 2, ..., + }. The node set is thus = {0, 0 } ∪ ∪ and | | = + + 1. . The edge set = {( , ) : , ∈ , < } contains the edges connecting vertices of . Each edge ( , ) is associated with two non-negative travel time cost and , which corresponds to the travel time needed for the truck and UAV to travel from node to node , respectively. In general, as UAV travels faster than the truck, is normally no less than . Besides, each edge is also associated with two non-negative energy cost and for the truck and the UAV, respectively. Note that both and are also measured in time unit. To simplify the problem, we further assume that for all the edges ≤ and ≤ , where and are energy capacity for the truck and drone, respectively. In this paper, as we assume that the truck is an electric vehicle. The term "truck", "electric vehicle", or "electric truck" are used interchangeably unless with further instruction. Similarly, the term "UAV" and "drone" are also used interchangeably. The planner aims to serve all the customers in the network using the truck or the drone. In the operational stage, after they depart from the depot, both the truck and the drone can serve customers independently. Besides, the truck serves as the drone hub to launch/retrieve the drone at a customer node. After launching the drone, the truck continues its planned route and serves customers before retrieving the drone at another customer node. Meanwhile, the drone will also perform delivery tasks, after which it will be retrieved by the truck. We assume that the launching/retrieving time is zero (this assumption will be relaxed later in this section). In practice, this can be done if the truck carries multiple drones and only launches the fully loaded/charged ones after retrieving the drone in operation. We define the process of the drone's operation from the launch to the next retrieve as "drone sortie" or "UAV sortie." In each drone sortie, we assume only one customer can be served, considering the low loading capacity of the drone. This indicates that each drone sortie contains exactly three nodes: launch node, customer node (also called "drone node"), and retrieve node. Besides, waiting time might be generated for each drone sortie when the truck and the drone arrive at the retrieve node at different times.

The goal of EVTSPD is to find a coordinated tour that starts and ends at the depot and visits a subset of vertices (including charging stations when necessary) such that the total delivery time is minimized. All customers must be served once, either by the EV or by the drone.

The final feasible solution of EVTSPD consists of the EV route and several non-overlapping UAV routes, which start and end on the EV route. A simple example is shown in Figure 2 

We define "drone node" as the node that is solely visited by the drone, "truck node" as the node that is solely visited by the truck, and "combined node" as the node that is visited by both the truck and the drone.

Below is a summary of parameters that are available in the original network :

: Travel time for the EV when traversing arc ( , ) : Travel time for the drone when traversing arc ( , ) : Travel time cost for the drone to launch from node , serves node and return to node : Driving energy cost for the EV when traversing arc ( , ) : Flight energy cost for the drone when traversing arc ( , ) : Flight energy cost for the drone to launch from node , serves node and return to node : Flight energy capacity of the drone : Driving energy capacity of the EV : A positive large number

One special assumption adopted in this research, distinct from previous VRP or TSPD research, is the shared energy assumption between the EV and the drone. To be more specific, we assume that there is a Lithium-Ion battery of limited capacity on the EV that can be used by both EV and drone. This battery not only serves as the energy source for the EV, it can also recharge the battery on the drone if it is on board. During the operation, if the drone is launched to serve a customer independently, the energy consumed by this drone sortie (which can be calculated beforehand) should be deducted from the remaining energy level of the battery.

To better illustrate this assumption, denote as the remaining energy of EV upon arrival at node and as the remaining energy of EV upon departure from node . Using the network in Figure 3 , assume the electric vehicle launches the UAV at node 1, travels to node 3, and retrieves the UAV at node 4 while the UAV serves the customer 2. If 1 = 100, which indicates that EV's remaining battery level upon arrival at node 1 is 100, then the EV's battery level upon departure from node 1 could be calculated as 1 = 100 − 124 = 100 − 20 = 80, where 124 represents the required energy of the UAV to be launched at node 1, serves customer 2 and returns to node 4. The battery level of the other nodes is also shown in Figure 3 .

There are two major reasons to make this assumption. First of all, unlike bigger ones used in military applications, the smaller drones, mainly used in logistics, usually have an inadequate power supply. An EV-drone coordinated delivery system can solve this problem as the battery on EV can also be used by the drone.

Second, simply ignoring the consumed energy during the drone's operation would lead to inaccurate EV driving range and overall result, as the drone's energy consumption rate is comparable to that of an EV. The average energy consumption rate for EVs on the market is around 200 wh/km (Electric Vehicle Database, 2022). For the specialized electric delivery van, whose driving range is about 86 km with 10 kWh LiFePO4 battery and 200 km with 20 kWh battery, the energy consumption rate is about 100-116 wh/km (Alke, 2022) . In terms of commercial drones, for example, DJI's MATRICE 600 PRO, when equipped with six TB47S batteries (4500 mAh, 22.2 V), can fly at most 16 minutes with a 6 kg payload, according to its official website (DJI, 2022) . Assuming this situation, it can fly at a maximum speed of about 60 km/h. The drone can travel up to 16 km before the batteries are depleted, which indicates that the energy consumption rate is about 37.5 wh/km, about 20% of a typical EV and 40% of a smaller electric van. The shared energy assumption aims to account for this extra 20%-40% percent of energy consumption from the drone. Besides, to consider the difference between EVs and drones in their energy consumption rate, we introduce another parameter to represent the energy consumption ratio between the drone and the EV. In the numerical analysis section, this value is chosen as = 0.4. In terms of the drone's charging time, in reality, it takes at least 30 minutes for the zero-energy-level battery to be fully charged. In practice, this time could be bypassed by swapping the used battery with the alternative backup batteries, which could be charged regardless of whether the drone is on-board or not.

In EVTSPD, the truck is assumed to be an electric vehicle associated with driving energy. As the truck traverses the arcs in the network, its remaining energy decreases gradually. Thus, the truck needs to visit charging station nodes when necessary to avoid its remaining energy being depleted. Note that no limit is set on the number of stops made for recharging at any charging stations. Thus, the final solution might contain more than one visit to some specific charging station nodes, while some other charging station nodes might never be visited at all.

In this paper, we assume that the electric vehicle is BEV and the charging stations are battery swapping stations (BSSs). For this setting, the BEV could renew its battery energy at battery swapping stations (BSSs) by swapping empty batteries with fully charged ones, which could be finished within minutes (Sarker et al., 2014) . Thus, the EV is expected to refresh its battery to full capacity with a fixed time every time it visits a charging station. This fixed-time-full-charge policy for EV is proposed in the early EV-related research and is commonly adopted in the literature (Erdoğan and Miller-Hooks, 2012, Conrad and Figliozzi, 2011) . EVTSPD which adopts this assumption is named as EVTSPD-FF.

Although this fixed-time-full-charge assumption might not be widely adopted in practice, this paper later shows that this setting enables us to create a multigraph that only contains customer nodes. The EVTSPD-FF can be defined on this multigraph instead of the original network. In other words, this assumption enables us to achieve significant computational efficiency that is otherwise impossible if this assumption is discarded. Later in this paper, we will show how to exploit this property to solve the problem exactly with a branch-and-price algorithm.

Before we present our MILP model for EVTSPD-FF, one thing to note is that this proposed MILP formulation is not defined on the original network but a constructed multigraph . We will explain the motivations and reasons for constructing this new network before the detailed steps of the multigraph construction.

Compared to the normal TSPD, the major complication in EVTSPD-FF is the extra EV driving range constraints: the planner needs to decide when and which charging station nodes should be visited before the driving energy is depleted. However, generating such a feasible route based on the original network is difficult. To enable multiple visits to the same charging station nodes, one needs to create several "dummy" copies of the charging station node set and insert these "dummy" nodes back into to construct an augmented network 0 . In (Andelmin and Bartolini, 2017) , the authors propose an efficient way to overcome this difficulty: reformulate the original network into a multigraph. In this research, we also adopt this technique to construct a multigraph for two purposes: it serves as an alternative network that EVTSPD-FF could be defined on. The multigraph is also used to generate feasible routes in the BP algorithm introduced later. For the rest of this subsection, we will explain the steps of constructing this multigraph. Denote = ( , ) as the subgraph containing the charging station node set , and arc set . This arc set contains all the arcs between the CS nodes, i.e., = {( , ) : , ∈ , ≤ }. Thus, network is a graph that only contains CS nodes and CS arcs. Let 0 = ∪ {0} be the set containing all customers and the single depot. Meanwhile, we define refuel path as a path that starts at ∈ 0 , ends at ∈ 0 , ≠ and traverses a path in graph . Besides, if a vehicle travels from node to node directly by traversing arc ( , ), the path is called a "direct path". Therefore, a path between two customer nodes is either a direct path or a refuel path. For every direct path we assign a index 0 to it so that a direct path between and is represented as ( , , 0). We also assign a unique index > 0 to every refuel path so that a refuel path between and is represented as ( , , ).

We can now define an extended arc set ℰ that contains all the refuel paths and direct arcs ( , ), ∀ , ∈ 0 . The refuel path that starts at and ends at is represented as arc ( , , ) in ℰ and the direct arc that starts at and ends at is represented as arc ( , , 0) in ℰ. In this way, the arc set ℰ contains all the paths that starts at a node ∈ 0 , end at another node ∈ 0 and visits no other customers along the route, regardless how many CS nodes are visited between and .

At this point, we can define a multigraph = ( 0 , ℰ) where node-set 0 = ∪ {0} is the set containing all customers and the single depot and ℰ is the extended arc set. Each path from to in original network will be represented as an arc from to in the multigraph . Compared to the original network , in the multigraph , the charging station nodes can be eliminated as they are inherently represented by each arc in . In this way, any original route in the original network can be represented as an alternative route in the multigraph , which has the same property as the original route.

Keep in mind that the main reason for constructing this multigraph is to create an efficient alternative network to avoid the construction of the augmented network, where the charging station nodes are duplicated. Also, note that as the multigraph contains all the direct arcs between nodes in 0 in the original network and we prohibit the drones from being launched or retrieved at the charging station nodes, any feasible drone sorties in the original network corresponds to the same drone sortie in . For this reason, both the EV's route and drone's sorties can be constructed on the multigraph . This indicates that EVTSPD-FF can be modeled on multigraph .

Next, we will introduce the properties associated with each arc in ℰ. Each arc ( , , ), > 0 that represents a refuel path is associated with three properties: minimum energy requirement ( , , ) , remaining energy ( , , ) , and travel time cost ( , , ) . Minimum energy requirement of arc ( , , ) is the minimum energy required at node to traverse this arc. Remaining energy of arc ( , , ) indicates the remaining energy of the truck after it traverses this arc. Travel time cost of arc ( , , ) is the time cost of EV when traversing this arc. To calculate these properties, denote 1 ( , , ) and 2 ( , , ) as the energy consumption of the first and the last arc in original network traversed by refuel path . Apparently, ( , , ) = 1 ( , , ) and ( , , ) = − 2 ( , , ) . As for ( , , ) , it is simply the sum of the travel time of all the original direct arcs in .

As for each arc ( , , 0) that represents a direct arc ( , ) in original network , it only has two properties: minimum energy requirement and travel time cost, which are calculated in the same way as described above. Keep in mind that as arc ( , , 0) ∈ corresponds to a direct arc ( , ) in , its first and last arc are the same. Also note that this arc has no associated remaining energy ( , , ) . This is because the truck visits no CS nodes when traversing direct arc. Thus, the remaining energy of EV after traversing direct arc ( , , 0) is not a fixed value, and thus no remaining energy is associated with this arc.

To better explain these properties in ℰ and , a simple case is shown in Figure 4 . In this figure, the first column represents the refuel path or direct arc in the original network , and the second column shows their corresponding arcs in multigraph .

Theoretically, the number of arcs in is prohibitive as the number of refuel paths between any two nodes is potentially infinite. However, a lot of these refuel paths can be eliminated from as they are dominated by other paths and will not appear in the optimal solution of the problem. Below we introduce the dominance rule among the refuel paths.

Theorem 1 Let = ( , , ..., , ) be a refuel path from ∈ to ∈ . This refuel path is represented as arc ( , , ) in . Let = ( , , ..., , ) be another refuel path from to . This refuel path is represented as arc ( , , ) in . Refuel path is said to be dominated by if the following conditions are satisfied: (i) ( , , ) ( , , ) and at least one of the three inequalities is strict.

The above dominance rule is easy to understand and prove, which indicates that a path is better than (or "dominates") another one if it has less energy requirement, greater remaining energy, and less travel time. As there exists at least one optimal solution to EVTSPD-FF whose traversed refuel paths are non-dominated, only non-dominated refuel paths are necessary when constructing the multigraph . For the rest of this section, the steps of the multigraph construction are explained in detail.

In general, the construction process is divided into two phases: the network elimination phase and the construction phase. In the elimination phase, some arcs that would not appear in the feasible solution of the problem are eliminated from the original network . The first step is straightforward. The second step was proposed in (Schneider et al., 2014) which involved eliminating arcs that would violate the fuel energy constraint when visiting the two closest CS nodes. 1) Remove from each arc ( , ) such that > 2) Remove from each arc ( , ) such that ∈ , ∈ , and + + > , ∀ , ∈ ∪ {0}; After the elimination phase, Algorithm 1 describes the the construction phase. In line 5, we derive all the feasible paths from to by checking three situations: the vehicle travels from to directly, the vehicle travels from to while traversing only one CS node, and the vehicle travels from to while traversing two or more CS nodes using . In line 6, we use dominance rule 1 to check if these feasible paths are non-dominated. Besides, in line 8, we need to assign a unique index to each non-dominated path to derive the route in the original network. Also, note that when adding arcs to , we need also store the arc's energy requirement, time cost, and remaining energy (if any).

It is worth noting that an upper bound on the number of arcs of is ( + 1) 2 2 where is the maximum number of nodes of the shortest path in between any two stations, and is the number of CS nodes. However, the number of arcs reduces significantly in practice with dominance rule 1 and Algorithm 1. Based on the numerical analysis, we Add all the remaining paths in to the multigraph as arc ( , , ), where is a unique index so that we can keep track of the exact route in the original network 9: end for 10: return Multigraph found that the number of arcs in is about 3-8 times that in the original network (only considering the arcs that traverse two customers nodes). (Schneider et al., 2014) provided a comprehensive analysis of the multigraph where the multigraph is compared with the extended graph, which assumes that the vehicle can make at most one refueling stop when traveling between two customers.

With the construction of , now we are ready to present the arc-based MILP formulation of EVTSPD-FF. Below is a summary of parameters or sets used in the model: : Equals one if customer node is a truck customer ∈ {0, 1} : Equals one if customer node is a drone customer ∈ {0, 1}

: Equals one if customer node is a combined customer : Arrival time of the truck or the drone or both at node ∈ 0 1 ≥ 0 : Theoretical remaining driving energy charge of the EV upon arrival at node , which is measured in time units. 2 ≥ 0 : Actual remaining driving energy charge of the EV upon arrival at node , which is measured in time units ≥ 0

: Remaining battery charge of the EV upon departure from node , which is measured in time units

Here we introduce two driving energy charges upon arrival at node : 1 and 2 . The first one is a theoretical value which aims to check if the remaining battery charge at previous node is sufficient to traverse the arc ( , , ). If 1 < 0, this indicates that the truck cannot traverse the arc ( , , ) from node . 2 is the actual remaining charge upon arrival at node . If arc ( , , ) : > 0 is traversed, then 2 = . If arc ( , , 0) is traversed, then 2 = 1 .

Objective 

In the above formulation, objective function (1) aims to minimize the final completion time. Constraints (2) are truck route flow balance constraints at node . Constraints (3) state that if node is visited by the truck then it is either a truck node or a combined node. Constraints (4) indicate that the truck must depart from and return to the depot. Constraints (5) are the drone route flow balance constraints at node . Constraints (6) state that if a drone traverses arc ( , , 0), then node is either a drone node or combined node. Constraints (7) indicate that the drone must departs from and returns to the depot. Constraints (8) state that node needs to be either truck node, drone node or combined node. Constraints (9) and (10) are the time coordination constraints for the truck and drone, respectively. Constraints (11) ensure that the drone travels from to if and only if the truck visits at least one of the two customers and . Constraints (12) prohibit sortie of type (0, , 0 ). Constraints (13) is the driving energy constraints of the truck. Constraints (14) and (15) calculate the actual arrival battery level if different arcs are traversed. Constraints (16) indicate that the truck departs from the depot with fully-charged battery. Constraints (17) are the shared energy constraints. Constraints (18) -(21) are the coordination constraints between and other decision variables.

In this section, we aim to relax some of the operational assumptions of EVTSPD-FF and explore how the arc-based model can handle these additional side constraints that are commonly seen in practice and literature. In this section, we only discuss some of the most common ones. For further discussion on this topic, the reader can refer to (Cavani et al., 2021) .

Assume that each node ∈ ∪ , is associated with a fixed time cost . If is a customer node, represents the fixed service time. If is a charging station node, represents the fixed charging time.

To account for this fixed node cost in the first model, we can simply add this fixed cost of node to the arc cost of the node's incoming or outgoing arcs. For example, define the modified travel times = + , ∀ ∈ ∪ . We can handle these side constraints by replacing the travel times with the modified travel times .

In this case, the launch node and the retrieve node of a drone sortie could be the same node. It indicates that the truck is allowed to wait at the launch node until the drone returns. With this assumption the drone sortie of the form ( − − ) is allowed which is also called a loop.

Define a new binary decision variable , which equals to 1 if the drone performs loop ( − − ) and 0 otherwise. To account for the impact of the loop, the objective function is changed to:

Constraints (7) are replaced by:

Constraints (16) are replaced by:

The new objective (22) considers the extra waiting time incurred by the loops. Constraints (23) indicate that a node can be a truck node, combined node, or a drone node that is visited by a drone sortie or a loop. Constraints (24) states that a node can perform loop only when it is a combined node. Constraints (25) is the updated shared energy constraints that consider the existence of loops.

A incompatible customer is a customer that can only be served by the truck and not by the drone. Denote ∈ be the set of customers that can be served by the drone. We can add the following constraints into the first model to account for incompatible customers:

= 0, ∀ ∈ \ , ∈ 0 : ≠

Assume that it takes the truck time to launch the drone, and during this process, the truck and the drone do not move. This assumption adds extra time cost to each launch node. We further assume that this launch time is not incurred if the launch node is the depot. To model launch time, we define a new binary variable , ∈ 0 , which equals one if node is a drone node and not directly visited by the drone from the depot. Then the objective (1) is replaced by:

Besides, we need to add constraints:

The new objective (28) considers the extra launch time of the loops and at each launch node. The constraints (29) states that equals one only when node is a drone node and not directly visited by the drone from the depot. Similarly, assume that it takes the truck time to retrieve the drone, and during this process, the truck and the drone do not move. We can model this retrieve time by replacing objective (1) with the new objective:

Solving the EVTSPD-FF of practical size with a commercial solver using the MILP formulations mentioned in the last section is prohibitive, considering the enormous number of constraints and decision variables involved. This section first introduces an exact branch-and-price solution algorithm. This solution method is similar to the work presented in (Roberti and Ruthmair, 2021) but with extra complications to consider the fuel energy constraint and the multigraph. The branch-and-price algorithm is based on the set partitioning formulation of the problem. However, in the case of EVTSPD-FF, solving the pricing problem exactly, namely, obtaining a coordinated route with the least reduced cost that visits each customer exactly once while satisfying all the driving/flight range constraints, would be as difficult as solving the original problem. Thus, a relaxation of the EVTSPD-FF route based on the well-known ng-route relaxation (Baldacci et al., 2011) is implemented. Using this relaxation, we might derive a "worse" lower bound than solving the pricing problem exactly. However, great computational benefits could be achieved, significantly increasing the algorithm's efficiency. This section presents the set partitioning formulation of EVTSPD-FF, followed by a detailed explanation of the ng-route relaxation adopted when solving the pricing problem using the dynamic programming technique. Then, the branch-and-bound strategy is introduced.

Although the BP algorithm is much more efficient than the commercial solver, it can only solve instances containing up to 10 customers within one hour, significantly less than its TSPD counterpart. Thus, at the end of this section, we also introduce a variable neighborhood search (VNS) heuristic, inspired by the work proposed in (Campuzano et al., 2021) , that aims to solve EVTSPD-FF instance of practical size.

An alternative formulation that has been widely used to model vehicle routing problems and its variants is that based on Set Partitioning (SP) or Set Covering (SC), which was initially proposed by (Balinski and Quandt, 1964) . This formulation uses a possibly exponential number of binary variables for each elementary circuit starting from the depot, traversing each customer, and returning to the depot. The use of SP is widely used to tackle any variants of the vehicle routing problem because all the additional side constraints could be defined in a feasible route.

The set partitioning formulation of EVTSPD-FF is also defined in the multigraph . Denote the coordinated route as , which includes a truck route and drone route . The truck route starts from the origin depot 0, traverses arcs = ( , +1 , ), ∈ , = 0, 1, ..., and returns to the destination depot. The drone route include several drone sorties that has launch node and retrieve node, which are traversed by the truck route.

The travel time cost of this route is ( ) = ( ) + , namely the sum of the travel times of the arcs it traverses and the extra waiting times at each (retrieve) node. The driving energy ( ) of this route at each vertex is:

In equation (31), when the truck/drone departs origin depot 0, its remaining energy is . If the truck traverses an arc ( −1 , , 0), then the remaining energy at node is the remaining energy at node −1 minus the consumed energy along the arc. If the truck traverses an arc ( −1 , , −1 ) where −1 ≠ 0, it means that the truck traverses refuel path −1 from −1 to . Thus the remaining energy at node is a fixed value ( −1 , , −1 ) associated with the traversed arc.

Meanwhile, a route is feasible if the following conditions are satisfied: 1) 0 ≤ ( ) ≤ , ∀ = 0, 1, ..., + 1 2) ≥ ( −1 , , −1 ) , ∀ = 0, 1, ..., These conditions indicate that whenever it aims to traverse an arc, the remaining energy of the truck should be no less than the required energy of the arc.

Let ℛ be the index set of all the feasible routes of EVTSPD-FF in . Thus, each feasible route consists of an ordered sequence of operations and combined legs (see their definitions in Section 3). Each such route ∈ ℛ is characterized by coefficients that indicate the number of times customer ∈ is visited in route and by that represents the total time cost of route . Let be the binary decision variable that takes the value 1 if and only if route is in the solution. The EVTSPD-FF could be modeled as the following set partitioning problem:

The objective function (32) aims at minimizing the duration of the selected route. Constraint (33) indicates the only one of such route can be selected. Constraints (34) guarantees each customer be visited exactly once and constraint (35) is the domain constraint for decision variables. Solving the formulation SP via a commercial solver would be prohibitive as there is an exponential number of feasible routes in the set ℛ. The most common approach to solving SP is to resort to the column generation (CG) or branch-and-price algorithm and generate feasible routes dynamically. By replacing the full set ℛ with the dynamically computed restricted route set ℛ , where is the index of the iterations, the linear relaxation of SP formulation is called the master problem, or in short, MP. The process of finding a possible route and adding it to ℛ is called a pricing problem. The basic idea of CG or BP is to find any basic feasible solutions in the pricing problem using the dual variable of the master problem and add such solutions to the restricted route set ℛ of the master problem until this process can be no longer continued. By replacing the full set ℛ with the restricted route set ℛ and solving MP, we get a lower bound of the EVTSPD-FF as it is a relaxation of the original problem. Then the branch-and-bound process is carried on to obtain the binary solution of SP and hence the optimal solution of EVTSPD-FF.

The critical step of this solution algorithm is to generate basic feasible solutions in the pricing problem, i.e., a route that visits each customer exactly once and satisfies all the drone and truck's related energy capacity constraints. Sadly, finding such a route is as complex as solving the original problem considering all the constraints involved. Thus, a relaxation technique is required to easily generate feasible routes at the cost of increasing the number of iterations.

One of the most efficient route relaxations in solving routing problems and its variants is the ng-route relaxation proposed in (Baldacci et al., 2011 ). An ng-route is not necessary an elementary route that visits each customer exactly once. In an ng-route, a customer could be visited multiple times, while some may not be visited at all. When a partial ng-route propagates to the next customer, one cannot visit if this customer lies in a dynamically computed set called ng-set. In other words, the ng-route can visit a customer multiple times if, in between these visits, the route visits another customer that is far away from .

Using the ng-route relaxation in the pricing problem of BP indicates that some ng-route routes that are not necessarily elementary are obtained in each iteration, which would be added back into the restricted route set ℛ of the master problem. Solving the master problem of SP enables us to get a lower bound of the original problem. Compared to solving the pricing problem exactly without using any relaxation technique, ng-route relaxation would result in a "bad" route which might lead to deteriorated lower bounding in each iteration of BP. However, the drawback of using this relaxation is not significant because most ng-route is typically cost-inefficient and unattractive and thus unlikely to appear in the MP solution, and only valuable ng-route would appear in the MP solution. Meanwhile, it is much easier to find such least-reduced-cost ng-routes than finding the least-reduced-cost elementary EVTSPD-FF route in the pricing problem, as few state variables are needed, and some strong dominance rules could be applied to reduce the state that needs to be checked. In general, adopting ng-route relaxation can lead to a more efficient algorithm. The rest of this section will illustrate how this technique solves EVTSPD-FF.

For each customer ∈ , define the ng-set of customer as the set that contains all the neighbors of customer , i.e., all the customers near . This ng-set indicates all the customers who cannot visit between two consecutive visit customer . The size of the ng-set determines the propagating direction of ng-routes. A larger ng-set indicates greater sub-tour allowed and better lower bounding from solving the master problem. Meanwhile, finding a least-cost ng-route with a greater ng-set size is more complicated. On the other hand, a small ng-set would result in an easy-to-solve pricing problem but may lead to a relatively worse lower bound. Based on the research conducted in (Roberti and Ruthmair, 2021) , setting the size of ng-set as 5 is a good trade-off between the quality of bounding and the difficulty of solving the pricing problem. This research tested three different ng-set sizes and reported their computational results in the numerical analysis section.

One critical step in the branch-and-price algorithm is to obtain the least-reduced-cost route given the dual information of the master problem solution. As mentioned above, the ng-route relaxation technique is used so that the non-elementary route is allowed to enter the restricted route set ℛ of the master problem. As a result, in the pricing problem of the BP algorithm, we aim to find a least-reduced-cost ng-route. This goal is accomplished by the dynamic programming (DP) algorithm introduced in this section. This algorithm combines elements of the solution method proposed in (Roberti and Ruthmair, 2021) and (Andelmin and Bartolini, 2017) with necessary modifications made to fit our purpose. The difference between these three DP algorithms will be discussed later in this section.

The problem description section demonstrates how an EVTSPD-FF solution route can be constructed to concatenate operations and combined legs. Remember that each such operation and combined leg serve their own customers. We define a feasible ng-route as a route whose total number of customer visits across all operations and combined legs is where is the total number of customers in the instance while satisfying the drone and truck's flight/driving energy limit constraints. The DP recursion used in this research is based on the idea of propagating the current partial route to the next available customer until the total customer visits equal to , and then the route returns to the destination depot. Such route is constructed in the multigraph and can propagate to the next customer by three possible actions: adding a truck arc, a drone leg, or a combined arc. By deploying these three actions, all the feasible ng-routes could be generated, and the least-cost route could be identified and added to ℛ . In this process, some strong dominance rules could be used to eliminate unpromising and dominated routes and thus facilitate the procedure.

Let 0 ∈ R and ∈ R, ∈ be the dual variables associated with constraints (33) and (34), respectively of the linear relaxation of SP formulation. Thus, in the pricing problem, the objective function is − − ∈ where is the route duration of route and is a parameter indicating the number of times customer is visited in route . By solving the pricing problem, we aim to find an ng-route with the least objective value. In DP, we define a function ( , , , , , , , ) as the representation of a partial ng-route. The elements of this function are described as follows:

• is the ng-set of the route, which indicates the set of customers that cannot be visited in the next propagation. It is a subset of customer set and also a subset of , ng-set of customer . • represents the number of customer visits performed so far. When equals , we either terminate the route or force the route to return to the depot. • ∈ is the last combined node visited by the truck in the partial route. • ∈ is the last customer node visited by the truck (regardless of whether the drone is on-board or not) in the partial route.

• ≥ 0 represents the truck's travel time since the last time the truck visits .

• represents the remaining driving energy of the truck at node . • represents the remaining driving energy of the truck at node . • is a binary indicator that represents whether is a fixed value when adding a truck arc. If the truck visit some CS nodes between and , then when adding another truck arc, the value will not change and =

. If no CS nodes exist between and , then = and = .

As we can see, we need to keep track of the remaining driving energy of the truck at both and . If = , it means that the partial route finishes an operation or traverses a combined arc. In this situation, equals and represents the maximum driving energy without charging. When ≠ , it means that the route has an ongoing operation and ≠ . If the route visits next customer by truck while traversing arc ( , , ) ∈ , we need to check if remaining energy at is sufficient to traverse this arc, that is, check if ≥ ( , , ) . Otherwise, if the route aims to visit customer by drone, it indicates that a drone leg − − would be performed. Because of the shared energy assumption, we need to check if ≥ + . Meanwhile, the binary parameter helps us to keep track of the value of .

The value of ( , , , , , , , ) represents the minimal duration of any partial route that starts from the depot and can be characterized by the tuple ( , , , , , , , ). In the beginning of the DP, the route is initialized as (∅, 0, 0, 0, 0, , , ) = − 0 . In the propagation process, given the current state ( , , , , , , , ) , there are three different actions can be taken based on the current partial route state: 1) Add truck arc: the current partial route will visit the next customer ∈ /( ∪ { , }) by truck only. However, this action can be taken if and only if the following conditions are satisfied:

• There exist an arc ( , , ) in to connect and . 2) Add combined arc: if = and = 0, it means that the current partial route just finishes an operation or combined arc. In this situation, the current partial route can visit the next customer ∈ \( ∪ { , } by adding a combined arc, if and only if the following conditions are satisfied:

• There exist an arc in to connect and . 3) Add drone leg: if ≠ and consequently ≠ 0, the current partial route has an ongoing unfinished operation.

In this situation, the current partial route can visit the next customer ∈ \( ∪ { , } by adding a drone A simple example is shown in Figure 5 to illustrate the three actions and their effect on the partial route introduced above. Note that in this figure, a double line arrow represents a combined arc and a single line arrow represents a truck arc. The three actions are deployed at every partial node to progress to the next customer node until the partial route visits exact customers. If the partial route has to visit customers, we also need to check if the route has returned to the depot, and an extra arc needs to be added if necessary. In this way, by exploring all the possible partial routes, we obtain an ng-route with the least reduced cost given the dual information of the master problem. However, the number of possible states still grows exponentially with the instance size, especially considering that multiple arcs exist between two nodes in the multigraph, and three possible actions exist. Thus, the following dominance rules are used to limit the number of states in the algorithm. During the algorithm process, these two dominance rules are constantly called to check if the new state dominates or is dominated by others. Let ( 1 , , , , , 1 , 1 , ) = 1 and ( 2 , , , , , 2 , 2 , ) = 2 . If the following conditions are satisfied:

Then ( 1 , , , , , 1 , 1 , ) dominates ( 2 , , , , , 2 , 2 , )

Proof: If all the five above conditions are satisfied, all the potential states that could be achieved by propagating ( 2 , , , , , 2 , 2 , ) can also be achieved by propagating ( 1 , , , , , 1 , 1 , ) with a lower cost. Thus the dominance rule stands.

Theorem 3 Let ( 1 , , , , , 1 , 1 , ) = 1 and ( 2 , , , , , 2 , 2 , ) = 2 . If the following conditions are satisfied:

Then ( 1 , , , , , 1 , 1 , ) dominates ( 2 , , , , , 2 , 2 , )

Proof: Define a reverse ng-route that starts with the destination depot while back propagating to origin depot. Thus we could also define a function ( , − , , , , , , ) = as the representation of this reverse partial ng-route. In this way, a full ng-route could be achieved by combining ( , , , , , , , ) with ( , − , , , , , , ). If all the above conditions are satisfied, all the potential propagation that could be achieved by propagating ( 2 , , , , , 2 , 2 , ) can also be achieved by propagating ( 1 , , , , , 1 , 1 , ). By combining 1 with backward partial route, we can derive a full ng-route. Assuming a route performs drone leg − − where is the drone node, denote = + + as the needed time for drone to finish the this drone leg. Then the reduced cost of the full route is:

Similarly, the reduced cost of the full route by combining 2 with iŝ

If all the conditions from (37a) to (37f) are satisfied, we need to prove thatˆ1 ≤ˆ2. This can be illustrated in three cases:

• If ≥ { 1 + , 2 + }, thenˆ1 = 1 + + − + andˆ2 = 2 + + − + . Thusˆ1 ≤ˆ2 because of condition (37a).

• If ≤ { 1 + , 2 + }, thenˆ1 = 1 + + − + 1 + andˆ2 = 2 + + − + 2 + . Thuŝ 1 ≤ˆ2 because of condition (37f). • If 1 + ≤ ≤ 2 + , thenˆ1 = 1 + + − + andˆ2 = 2 + + − + 2 + . Thusˆ1 ≤ˆ2 because of condition (37a). Thus,ˆ1 ≤ˆ2 in all the scenarios and the dominance rule stands.

The launch and retrieve time assumption can be incorporated in the DP model by changing the state value calculation to account for the launch and retrieve times. When adding a truck arc, the new state value is ( , , , , , , , )− + . When adding a drone leg, the new state value is ( , , , , , , , )− + ( , , ) + .

In dynamic programming, we can add a truck arc when ! = and there is enough energy to traverse the next truck arc. When the maximum number of customers per truck leg assumption is enforced, we need to check an additional condition, that is, the current number of customers in the current truck leg, to ensure a feasible route. We need an extra element to record this value when constructing states in DP to do so. Define a new state variable indicating the current number of customers visited in the current truck leg. In this way, the new state is ( , , , , , , , , ) .

Note that it might seem non-trivial to add a state variable in DP as this would lead to expanded state space. However, in practice, when the upper bound¯is small (e.g.¯≤ 3), the influence of the new state variable is minimal. Besides, as places an upper bound of the number of consecutive "add truck arc" actions, the propagation process and hence the overall computation process is shortened, as suggested in the numerical analysis section.

The weight-dependent drone flying range assumption can be enforced in DP by changing the original constraint ( , ,0) + ( , ,0) ≤ when adding drone leg to the new constraint ( , ,0) + ( , ,0) ≤ where ( , , ) is the energy consumption function for drones when traversing arc ( , , ).

When loop is allowed in the operation process, in DP, an extra action of adding loop is needed. For current state ( , , , , , , , ), the add loop action is only available when = , = 0, = , = , which indicates that both the truck and the drone are at node at current state, and ≥ + which indicates that the remaining energy is sufficient for drone to perform sortie < , , >. After the add loop action, the resulting new state is (( ∪ { , }) ∩ , + 1, , , 0, , , ) and the new state value is ( , , , , 0, , , ) + ( , ,0) + ( , ,0) − .

With the ng-route relaxation technique, the columns (routes) entered into the master problem of BP might not be elementary routes. As a result, the final solution of the master problem might be fractional. If this situation happens, a branch-and-bound procedure is carried out to obtain a feasible integer solution, and thus a branching strategy guiding this process is needed. For a node in the B&B tree, given an optimal solution * of the master problem of SP formulation and its corresponding column set (ng-route set), we follow the routine introduced in (Roberti and Ruthmair, 2021) and perform a binary search based on three different types of decisions: (1) * , ∈ , which represents the times each customer is served by the drone alone in the optimal solution, (2) * , ( , , ) ∈ , which represents the times each arc ( , , ) is traversed by the truck only in the optimal solution, and (3) * , ( , ) ∈ , which represents the times each arc ( , ) is traversed by the drone only in the optimal solution. Note that the value of * , * , * should all lie in the interval [0, 1] because of constraints (35).

With the three kinds of optimal solution values, we can perform branching based on these values in a hierarchical manner. First of all, we check if * , ∈ is binary for each customer . If there exists some customer whose * value is not binary, we then find the one whose value is closest to 0.5, that is, * = ∈ | * − 0.5| and branch on this customer. If all the values of * are binary, we can move to the values of * , ( , , ) ∈ and similarly branch on the truck arc whose value is closest to 0.5, that is, ( , , ) * = ( , , ) ∈ | * − 0.5|. If all the values of * are either zero or one, we then move to the last layer of branching decisions * , ( , ) ∈ and branch on the drone arc whose value is the closest to 0.5, that is, ( , ) * = ( , ) ∈ | * − 0.5|. Note that we end up with two child nodes with corresponding branching decision values at zero or one for each of these branching decisions. After performing all three layers of branching procedures, we can obtain all the feasible integral solutions to the original EVTSPD-FF and pick the least-cost one as the optimal solution. To increase the efficiency of this branching process, one can apply a depth-first-search method to obtain an integral solution as soon as possible so that it can be used as an upper bound to prune the remaining nodes in the later process.

As the solution EVTSPD-FF solely consists of truck arcs and drone legs, it seems unnecessary to branch on the value of * , ∈ . However, as explained in (Roberti and Ruthmair, 2021) , only branching on * and * might lead to an unbalanced search tree as it is more restricting to set the value of a specific * or * to one than setting the value to zero. Besides, based on the numerical results conducted in this paper, we find that for most small-sized cases, merely branching on * , ∈ is sufficient to obtain an integral solution. For these reasons, the first layer of the branching decision is on * , ∈ .

VNS scheme is a classical meta-heuristic that aims to escape local optima by changing the way of selecting neighbourhood structure during the process. This method is first proposed in (Mladenović and Hansen, 1997) and later improved in (Hansen and Mladenović, 2001) . Over the past few decades, this approach has been successfully applied in solving multiple problems, especially in the field of optimization. In the VRP community, it can be used to solve different variant of the problem such as TSP, VRP with time window and VRP with multi-depot. (Campuzano et al., 2021) compare the performance of VNS with the MILP formulation proposed in (Cavani et al., 2021) to solve TSPD. Thus, in this section, we also present a similar VNS scheme to solve EVTSPD-FF instance of larger size. The general framework of VNS is shown in Algorithm 2.

Unlike TSPD where any permutation of the customers nodes is a feasible route, in EVTSPD-FF, finding an feasible initial solution is non-trivial, especially when the EV's energy constraints are tight. A simple initial solution of EVTSPD-FF is using the EV to serve all the customers and ignore the drone completely. However, it is still non-trivial to obtain a feasible EV route. In the past literature, there exist various ways of finding a feasible solution. Some examples are the Modified Clarke-Wright (MCWS) algorithm proposed in (Erdoğan and Miller-Hooks, 2012) , dynamic programming approach introduced in BP algorithm and PseudoGreedy algorithm, a greedy method proposed in (Felipe et al., 2014) . In this research, the MCWS method is used to generate a feasible EV route in VNS due to its simplicity.

As a critical component of the VNS algorithm, the neighborhood structure is used to escape the local optima so that the algorithm can explore the new space in the feasible region. The performance and effectiveness of the neighborhood operator directly affect the algorithm's overall performance. In this research, we use a variety of neighborhood operators whose performance has been tested in the past literature. All the used operators used in this research are shown in Table  1 : Swap location of three consecutive nodes with other two node (Campuzano et al., 2021) ReinsertCS Delete one random CS node from EV's route and insert another one This paper

As mentioned previously, the major complication of EVTSPD-FF compared to its TSPD counterpart is that it is difficult for EVTSPD-FF to maintain the solution's feasibility during the search process. In VNS, a saving-based greedy method is used to insert charging stations into found solutions so that its feasibility can be maintained.

Note that maintaining solution feasibility is a costly process in VNS. Checking the solution's feasibility is of computational complexity O ( ) where is the number of nodes in the current EV's route, while inserting one CS node into the current solution is of computational complexity O ( | |) where | | is the number of CS nodes in the network. Meanwhile, most of the neighborhood structure operator requires at most quadratic time O ( 2 ). This indicates that, compared to TSPD, for each new solution found in the searching process, the number of operations needed to maintain the solution's feasibility is at least doubled. In practice, the increased complexity is even more significant because the inserting CS function might be called several times to obtain a feasible solution.

In this section, the effectiveness of the proposed MILP formulation, the efficiency of the exact BP method and the performance of the VNS heuristic are all compared and tested. Then, a real world case study is conducted on the downtown Austin network to demonstrate the practicality of the proposed model.

While benchmark instances are available for drone related routing problems , without the location of charging station nodes, these instances are not guaranteed to be feasible in the case of EVTSPD-FF with randomly generated charging stations. Thus, in this paper, we conduct a numerical analysis on randomly generated instances with realistic parameter setting.

In this test, for all the generated instances, the single depot is located at (0, 0), and the coordinates of the customer nodes and the charging station nodes are uniformly distributed between -20 km and 20 km. The number of customers | | is an input parameter while the number of charging station nodes is roughly a third of | |. The distance matrix for the electric vehicle is calculated using the Manhattan metric, while the Euclidean metric is used for the UAV to reflect its greater mobility options. For each randomly generated instance, the MCWS algorithm is used to solve the instance as electrical vehicle TSP to guarantee the instance is feasible for EVTSPD-FF.

In this test, we assume the electric van is Alke model ATX340E and the drone is DJI MATRICE 600 PRO. The parameter setting associated with the electric van and drone are mostly adopted from the company's official website (Alke, 2022 , DJI, 2022 . As per the Federal Aviation Administration's regulation, the maximum allowed altitude is 400 feet (122 meters) above ground level. In this research, the drone's cruise altitude is estimated to be 50 meters. The drone's launch time is estimated to be the cruise altitude divided by the taking off vertical speed. In the table, the launch time also includes the time need for swapping battery and reloading parcel, which indicate that the launch time is greater than the retrieve time. The proposed branch-and-price algorithm and ALNS are coded in Python, while the MILP formulation is implemented in Pyomo and solved with ILOG's CPLEX Concert Technology solver (version 12.6.3). All experiments are run on a 3.6 GHz Intel Core i7 desktop with 32 GB RAM.

We first compare the performance between the two exact solution methods: solving the MILP formulation with commercial solver (CPLEX) and the proposed BP method. The computational results are shown in Table 3 , which is a summary for the average test results for instances with different size (| | = 5, 6, 7, 8, 9, 10) and different drone-truck travel speed ratio ( = 1.5, 2, 3). In the table, Column No.ins reports the number of instances tested. Arc represents the result from the arc-based MILP model. As mentioned before, BP with three different ng-set sizes is implemented. In the column Opt, the value indicates the number of instances solved by different approaches. Note that the exact optimal solution is not available for large instances ( ≥ 8). In this case, the solution with the lowest total cost is considered the optimal solution among all the returned solutions of all the four solution methods. Column RunningTime reveals the average operational time of different approaches over the solved instances, which is measured in seconds. Column BB nodes presents the number of branch-and-bound nodes searched to obtain the optimal solution. The last column, GapRoot reveals the average gap at the root node between the optimal solution cost and the lower bound returned by the root node. Note that only the instance that is solved by the corresponding approach is considered for each approach. As can be seen from Table 3 , it is obvious that the branch-and-price algorithm is more efficient than the proposed MILP models defined on the multigraph, regardless of the value of the parameter . Specifically, the arc-based model can solve instances containing 8 customers while the BP method can solve instances containing up to 10 customers in one hour.

As we mentioned earlier, BP-5 is the most efficient in solving the pricing problem among all three BP methods, while BP-7 can obtain the best lower bounding at the root node. This property is also demonstrated by the results in Table 3 . According to the table, BP-7 solves the most cases among the three. This is mainly because BP-7 needs to examine fewer nodes before termination in the branch-and-bound procedure by obtaining a better lower bound at the root node. This result is also illustrated in the "BB nodes" column, as the average BB nodes for BP-7 is much lower than the other two methods.

Besides, comparing the performance of the three BP methods, it seems that when ≤ 7, BP-is the most efficient method and BP-5 is more efficient than the other two methods when ≥ 8. In the first case, when ≤ 7, the branch-and-bound tree for method BP-only contains the root node, which indicates that the problem could be solved in one iteration by solving the root node. However, when ≥ 8, all three BP methods need to examine other branch-and-bound nodes other than the root node to terminate the algorithm. As the results suggest, as the number of customers increases, the average BB nodes in the BP method increase accordingly, from less than 3 when ≤ 7 to around 8 when = 10. Meanwhile, when ≥ 8, the average computational time for BP-5 is much lower than the other two methods, as it is simpler to solve the pricing problem for BP-5 than BP-6 and BP-7.

Furthermore, the test results with various drone-truck travel speed ratio values indicate that the value of parameter has little impact on the overall computational efficiency of all three solution methods. This is likely because the value of only affects the number of feasible drone sorties in the instances whose impact on the total number of decision variables is minimal.

Besides, we also conduct additional tests on the performance of MILP model and BP5 method on some commonly seen variants of EVTSPD-FF. The following EVTSPD-FF variants are considered:

• Base: the base case of EVTSPD-FF without any side constraints • LRT: the EVTSPD-FF variant considers the launch/retrieve time, whose value are shown in Table 2 • Range: the EVTSPD-FF variant which assumes the drone's flight range is dependent on the weight of parcel it carries. In this research, we adopts a simplified assumption that its flight range is doubled when no parcel is on boarded and that the range decreases linearly with the weight of the parcel until it reaches the base range value with full-capacity payload. • Loop: the EVTSPD-FF variant which allows the self-loop The computational results are shown in Table 4 , where the unit for each value is second. It can been seen that the BP method is very reliable in handling additional side constraints as it is fairly simple to consider these extra side constraints in dynamic programming propagation. For the MILP model, the computational time of maximum leg variant and self-loop variant is about 35% higher than normal EVTSPD-FF. 

In this section, we analyze the performance of the proposed variable neighourhood search heuristic and compare it with the exact branch-and-cut method. The parameter setting of the tested instances are the same as previously described. The test results are shown in Table 5 . In the table, the unit for the final cost columns is 10 seconds and the optimality gap of VNS is calculated as the difference between final results of two methods divided by the optimal cost.

As can be seen from the table, the performance of VNS is reasonably acceptable as the optimality gap is less than 3.5% for instances containing less than 8 customers and 3 charging stations while the running time of VNS is significantly less than that of BP. For instances with more than 10 customers, the BP fails to terminate in 1 hour while the VNS could solve instances with 25 customers in about 1 minute. Although we expect the optimality gap of VNS to grow as the size of the instances increases, the current test results indicate the VNS is competitive to exact solution methods. 

In this section, a real-world case study is conducted to illustrate the effectiveness of the proposed algorithm in solving EVTSPD-FF with practical size. This study also analyzes the effects of several key parameters on the final delivery time. In this case study, the downtown Austin network is analyzed. A total of 35 nodes are chosen randomly from the network, most of which are the centroids of ZIP Code Tabulated Areas (ZCTAs) in Austin metro area, while ensuring that the resulting EVTSPD-FF problem is feasible. Out of these 35 nodes, 10 nodes are selected as the charging stations and the remaining 25 nodes are considered as customer nodes or depot. Two different layouts are test. In the first layout we pick a node that lies in the central area of the network as the depot while in the second layout the depot is chosen as a node that lies in the suburban district of Austin. These two layouts are named as central layout and side layout, respectively. Thus, the final network contains one depot, 24 customers and 10 charging stations. All customers can be served by either the EV or the UAV. The travel time of the EV from one node to another node is estimated using the Google map Python API for the peak hour of a typical Monday. Thus, the resulting travel time matrix is asymmetric. Besides, the travel time of the UAV from one node to another is calculated as the direct distance between the two nodes divided by the drone's speed, which is 60 km/h. All the other parameters are the same as the previously described. A representation of the network is shown in Figure 6 .

With the default setting, the completion delivery time obtained via the proposed VNS heuristic is 16040 seconds (approximately 4.5 hours), and the final route consists of 14 customers being served by the UAV and ten customers being served by the EV.

In the rest of this section, we study the impact of characteristics of different parameters on the performance of this EV-UAV delivery system, from which we will gain insights into this new mode of transportation. All the results are the best-found solution based on ten independent runs of VNS.

The UAV is typically faster than the EV as it can take shortcuts without being affected by ground traffic congestion. The X-Wing project initiated by Google declares their delivery UAVs can reach the speed of 120 km/h when tested in a suburban area in Australia. This section explores the effect of the drone's travelling speed on the final delivery cost.

Note that the default travel speed is 60km/h and three additional speed are tested, namely, 50 km/h, 70 km/h and 80 km/h. The final solution cost of different drone travelling speed is shown in Figure 7 .

As shown, the delivery completion time decreases as the travelling speed of the UAV increases, for both layouts. As the UAV speed increases, the total delivery time, the EV travel time, and the total waiting time reduce gradually. For central layout, the final route completion cost decreases about 16.67% when drone's speed increases from 50 km/h to 80 km/h. For side layout this value about 11.1%. This result indicates that increasing the UAV speed is an effective approach to reduce the risk and operational cost incurred by the waiting period during the delivery.

The driving range of an EV depends on the size of the battery it carries. However, as battery size increases, the operational and maintenance cost also grows, so logistics companies deploying electric vehicles must make a trade-off between these two factors. In EVTSPD-FF, the driving range of the EV affects the final delivery completion time in that the low driving range indicates EV has to visit charging stations frequently. Tighter energy constraints also make it difficult for the heuristic algorithm to find a "good" feasible solution. Denote the EV's driving range as , which is measured in seconds. This section presents the result of how EV's driving range affects the solution route cost, as illustrated in Figure 8 . As presented, for side layout the final delivery time decreases with the EV driving range increases. However, this trend is not obvious for central layout. 

This paper introduces the electric vehicle traveling salesman problem with the drone, a new delivery concept combining electric vehicle and drone. In EVTSPD, the electric vehicle and drone operate in a cooperative way, where the EV serves as a drone hub that launches and retrieves the drone similar to the FSTSP introduced in (Murray and Chu, 2015) and TSPD introduced in . Meanwhile, the unique features of the EVTSPD, compared to the past research, are two-fold. First, as the electric vehicle is modeled in the problem, the truck needs to visit charging station nodes to avoid depleting its energy. Second, we propose a novel shared energy assumption, which indicates that the EV and drone shared their energy in the operation process. EVTSPD aims to find a coordinated EV-drone route that minimizes total travel time and serves a set of customers while incorporating stops at charging stations in route plans to ensure sufficient charge.

In this paper, the EV is assumed to be fully-charged with fixed time at charging station nodes. We present a novel arc-based mixed integer/linear programming formulation for EVTSPD-FF defined in a constructed multigraph. In this multigraph, all the feasible route between two non-charging station nodes are pre-calculated and stored in the network as an arc between the two nodes. In this way, the dimension of the multigraph is significantly smaller than the original network which needs to be augmented to enable multiple visit to one single charging station. Later in this paper we show that computational efficiency could be achieved by utilizing this multigraph and an exact branch-and-price algorithm is proposed to solve the EVTSPD-FF. In this BP algorithm, a dynamic programming method is used to propagate in the multigraph to obtain a feasible ng-route of the problem. This research also presents a variable neighbourhood search heuristic to solve the problem.

In the numerical analysis section, the computational efficiency of the MILP formulation, exact BP algorithm, and the proposed VNS heuristic are tested and compared. The test results indicate that the BP method is much more efficient than solving the MILP model via a commercial solver. The VNS is the fastest of all three method with reasonable optimality gap. Considering the enormous amount of constraints in the problem, the exact method can only be used to solve small-sized instances and for instances containing more than ten customers the VNS is the only viable method. Besides, a real-world case study based on downtown Austin network is performed to show that the VNS could solve the problem with practical size within a reasonable computational time. It also conducts several sensitivity analyses to illustrate how key parameters affect the final integrated route. The test results demonstrate that the UAV speed has a more significant influence on the final delivery time compared to the EV driving range.

The EVTSPD-FF formulations, along with these solution techniques, will aid organizations with EV fleets in overcoming difficulties resulting from limited recharging infrastructure. The new delivery concept of using UAV and EV to perform last-mile delivery would result in financial and environmental benefits when considering the reduced operation cost of fueling and switching to drone, which does not require a costly human pilot. Besides, the research in the newly-merged approach will provide intuition for the future development of a more sophisticated delivery service.

There remain several practical challenges for drone delivery regarding payload capacity, safety, and public acceptance (Watkins et al., 2020) . It remains to be seen whether battery-swapping stations (like the type we assume) or fast-charging stations ultimately make more economic sense for logistics fleets. As the number of customers increases, it will become essential to consider a multi-vehicle version of the current problem, perhaps with heterogeneous EV and drone range and capacity or alternative (non full-charge) policies. All of these should be addressed in future research.

Electric vehicle routing problem

Study of sidewalk autonomous delivery robots and their potential impacts on freight efficiency and travel

Persistent UAV delivery logistics: MILP formulation and efficient heuristic

Crowdsourcing delivery: New interconnected business models to reinvent delivery

Delivery by drone: An evaluation of unmanned aerial vehicle technology in reducing CO2 emissions in the delivery service industry

Drones For Deliveries From Medicine To Post

Teal Group Predicts Worldwide Civil Drone Production Will Soar Over the Next Decade

Rivian Electric Van Could Be Called RCV

The flying sidekick traveling salesman problem: Optimization of drone-assisted parcel delivery

The vehicle routing problem: State of the art classification and review

A literature review on the vehicle routing problem with multiple depots

A green vehicle routing problem

Green vehicle routing problem: A state-of-the-art review

A survey on the electric vehicle routing problem: variants and solution approaches

Optimization of a distribution network using electric vehicles: A VRP problem

The recharging vehicle routing problem

The electric vehicle-routing problem with time windows and recharging stations

A new mathematical programming model for the green vehicle routing problem

Dynamic lithium-ion battery model for system simulation

Experimental validation of a battery dynamic model for EV applications

The electric vehicle routing problem with nonlinear charging function

Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions

The electric vehicle routing problem with shared charging stations

Electric vehicle routing with public charging stations

The green vehicle routing problem: A heuristic based exact solution approach

Exact algorithms for electric vehicle-routing problems with time windows

The electric fleet size and mix vehicle routing problem with time windows and recharging stations

An exact algorithm for the green vehicle routing problem

An exact algorithm for the electric-vehicle routing problem with nonlinear charging time

Optimization approaches for civil applications of unmanned aerial vehicles (UAVs) or aerial drones: A survey

Drone-aided routing: A literature review

Vehicle routing problems for drone delivery

A multi-objective green UAV routing problem

An optimization-driven dynamic vehicle routing algorithm for on-demand meal delivery using drones

Analysis and Design of Band-Pass Filter Based on Metamaterial

Efficient large-scale multi-drone delivery using transit networks

Optimization of a truck-drone in tandem delivery network using k-means and genetic algorithm

Optimization approaches for the traveling salesman problem with drone

Dynamic programming approaches for the traveling salesman problem with drone

A decomposition-based iterative optimization algorithm for traveling salesman problem with drone

A randomized variable neighborhood descent heuristic to solve the flying sidekick traveling salesman problem

On the min-cost traveling salesman problem with drone

Coordinated logistics with a truck and a drone

Joint optimization of customer location clustering and drone-based routing for last-mile deliveries

The multiple flying sidekicks traveling salesman problem with variable drone speeds

The vehicle routing problem with drones: several worst-case results

Algorithms for solving the vehicle routing problem with drones

Vehicle routing problem with drones

A hybrid VNS/Tabu search algorithm for solving the vehicle routing problem with drones and en route operations

An adaptive large neighborhood search metaheuristic for the vehicle routing problem with drones

Two echelon vehicle routing problem with drones in last mile delivery

Two-echelon routing problem for parcel delivery by cooperated truck and drone

Energy consumption of full electric vehicles

Electric Van Alke

DJI MATRICE 600 PRO

Optimal operation and services scheduling for an electric vehicle battery swapping station

Exact methods for the traveling salesman problem with multiple drones

Exact methods for the traveling salesman problem with drone

New route relaxation and pricing strategies for the vehicle routing problem

A Multi-start VNS Algorithm for the TSP-D with Energy Constraints

On an integer program for a delivery problem

Variable neighborhood search

Variable neighborhood search: Principles and applications

A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges

A variable neighborhood search for flying sidekick traveling salesman problem

Ten questions concerning the use of drones in urban environments

This research is based on work supported by the National Science Foundation under Grant No. 1826230, 1562109/1562291, 1562109, 1826337, 1636154, and 1254921. This work is also supported by the Center for Advanced Multimodal Mobility Solutions and Education (CAMMSE).