Sage Journals: Discover world-class research

Abstract

Drones and electric vehicles (EVs) represent promising technologies for enhancing the efficiency and sustainability of last-mile delivery services. This paper focuses on the optimization of customer deliveries through the integration of a plug-in hybrid electric vehicle (PHEV) and a drone. Our model, named the plug-in hybrid electric vehicle traveling-salesman problem with drone (PHEVTSPD), assumes the PHEV can be recharged, either fully or partially, at charging stations, while the drone can be launched or retrieved from the EV. Both the EV and drone are capable of independently serving the customer. In comparison with traditional truck-only or drone-only delivery models, the hybrid EV–drone model overcomes the limitations of drone payload capacity and EV service area, thereby significantly improving delivery efficiency and reducing greenhouse-gas emissions. This research presents a three-index mixed-integer linear program (MILP) formulation of PHEVTSPD. Additionally, a linear or piecewise linear approximation of the concave time–state-of-charge (SoC) function is adopted in the model. To solve the proposed problem, we introduce an adaptive large-neighborhood search (ALNS) metaheuristic. Numerical analysis results reveal that the proposed ALNS method outperforms variable neighborhood search (VNS) with an average optimality gap of approximately 3% when solving instances with 10 nodes. Furthermore, a piecewise linear function with a six-line-segment approximation demonstrates an average of 10.8% lower cost compared with a linear approximation.

Keywords

aviation urban air mobility freight systems urban freight transportation logistics

In the United States, the transportation sector generates 28% of the national greenhouse-gas emissions ( 1 ). Many local governments and corporate policies aim to promote transportation modes with lower pollution and greenhouse-gas emissions. Electric vehicles (EVs) are an emerging alternative to internal combustion engines, and several companies have started using them in their operations. For example, in 2018, FedEx announced a fleet expansion and added 1000 electric delivery vehicles to operate commercial and residential pickups and delivery services in the United States ( 2 ). Switching to electric fleets not only has long-term effects on mitigating the impact of climate change but may also have immediate financial benefits, as fuel cost accounts for 39% to 60% of operating costs in the trucking sector ( 3 ). Compared with conventional internal combustion engines and petroleum-fuel-powered vehicles, EVs are much more energy-efficient and require less maintenance, which indicates potential savings to freight and logistics companies ( 4 , 5 ). Another new trend in recent years is the integration of unmanned aerial vehicles (UAVs), or drones, into e-commerce and on-demand item delivery. Amazon, Google, and DHL announced plans to use UAVs to deliver small packages. The past few years have witnessed a dramatic increase in UAV applications ( 6 ). According to Teal Group’s prediction, commercial use of UAVs will grow eightfold over the decade to reach U.S. $7.3 billion in 2027 ( 7 ). For simplicity, for the rest of the paper, the terms “UAV” and “drone” are used interchangeably.

While EVs and drones exhibit considerable potential for last-mile delivery services, they are not without their limitations. Despite recent strides in battery technology allowing private electric vehicles and trucks to achieve a driving range of 300 to 400 mi per charge, the practical driving range for EVs dedicated to commercial use, such as delivery vans, hovers around 150 mi in real-world scenarios ( 8 ). Consequently, en-route charging becomes a necessity, particularly when navigating suburban or rural areas or executing cross-town routes. Simultaneously, the application of drones faces a significant constraint in their payload capacity, notably for flat-wing drones commonly deployed in delivery operations. For instance, drones utilized in Amazon Prime Air typically carry payloads of around 5 lb ( 9 ).

In response to these challenges, this paper introduces a hybrid EV–drone delivery model designed to overcome the aforementioned drawbacks. Similar to the traveling-salesman problem with drone (TSPD) proposed in Murray and Chu ( 10 ) and Agatz et al. ( 11 ), this model pairs an EV with a drone for efficient delivery tasks. In contrast to traditional truck-only or drone-only delivery approaches, the hybrid EV–drone model adeptly overcomes the restricted payload capacity of drones and the limited service area of EVs, substantially enhancing delivery efficiency and mitigating greenhouse-gas emissions. Notably, both the EV and the drone operate independently to serve customers while the EV serves as a drone hub, facilitating the recharge and launch/retrieval of the drone. This setup allows the drone to be deployed multiple times en route, reaching hard-to-access customers by flying directly without being hindered by ground congestion. The battery energy can be replenished en route at charging stations within the network. This novel problem is termed the plug-in hybrid electric vehicle traveling-salesman problem with drone, or PHEVTSPD.

This article is based on Chapter 4 of the author’s PhD thesis ( 12 ) and technical report ( 13 ) and is an extension of Zhu et al. ( 14 ). Specifically, the main contributions of this paper are:

Integration of EV and drone routing: The proposed problem, PHEVTSPD, is one of the few pieces of research that address the EV routing problem and drone routing problem simultaneously, providing a comprehensive framework for optimizing hybrid transportation systems.

Partial recharge with nonlinear state-of-charge (SoC) function: PHEVTSPD considers partial recharge of EVs, incorporating a nonlinear concave time–SoC function. This function models the relationship between charging time and energy charge amount using a piecewise linear approximation, similar to Montoya et al. ( 15 ).

MILP formulation: A mixed-integer linear programming (MILP) formulation for the PHEVTSPD, defined in an augmented network, is proposed. This formulation captures the complexity of both vehicle and drone routing, as well as the charging dynamics of EVs.

Adaptive large-neighborhood search (ALNS): A specially designed ALNS method, which incorporates constraint-programming (CP) modeling as a subproblem, is presented. Test results indicate that this approach is more efficient than the variable neighborhood search (VNS) method, which is widely used to solve the TSPD.

Real-world case study: A real-world case study based on the downtown Austin network demonstrates that the final completion time increases with the accuracy of the piecewise linear approximation. This highlights the practical applicability and effectiveness of the proposed methods in realistic scenarios.

The rest of the paper is organized as follows. The next section presents the literature review. The MILP formulation for PHEVTSPD is described after that, followed by the proposed ALNS method. The subsequent section compares the computational performance of the MILP formulation and ALNS algorithm; it also includes sensitivity analysis of several key parameters. The final section presents the conclusion and potential future research directions.

Literature Review

The vehicle routing problem (VRP) and the traveling-salesman problem (TSP), first proposed by Dantzig and Ramser ( 16 ), are among the most well-studied optimization problems in operations research. Researchers have considered many VRP and TSP variants—incorporating service time windows, capacities, maximum route lengths, distinguishing pickups and deliveries, fleet inhomogeneities, and so forth. Various exact and heuristic methods have been proposed to solve the problem ( 17 – 21 ). References ( 22 – 27 ) and ( 28 ) provides a thorough literature review of VRP variants and solution-algorithm families.

Erdoğan and Miller-Hooks ( 29 ) introduced the green VRP (G-VRP), in which the goal is to route a fleet of alternative-fueled vehicles (AFV) to serve a set of customers within a time limit while respecting the driving range of the vehicles. Conrad and Figliozzi ( 30 ) developed the recharging VRP in which vehicles can recharge at particular customer locations while considering customer time windows and fleet capacity constraints. Schneider et al. ( 31 ) also modeled customer time windows and fleet capacity constraints in the electric VRP (E-VRP) problem while making the recharging times dependent on the remaining charge levels. Montoya et al. ( 15 ) extended the previous E-VRP models to consider nonlinear charging functions using piecewise linear approximations. The solution methods for E-VRP variants are diverse and range from exact methods such as branch and bound, branch and cut ( 32 ), and branch and price ( 31 , 33 ) to heuristic methods such as the modified savings method of Clarke and Wright with density-based clustering ( 29 ), local improvement based on neighborhood swap ( 31 , 34 ), and metaheuristics such as simulated annealing and tabu search ( 35 – 37 ). Pelletier et al. ( 38 ) and Erdelić and Carić ( 39 ) provide a comprehensive survey of the different E-VRP variants and associated solution algorithms. Unlike the research mentioned above, this paper focuses on the traveling-salesman variant. Doppstadt et al. ( 40 ) formulated the TSP for hybrid EVs considering four modes of operation—combustion, electric, charging, and boost. An iterated tabu search with local search operators, which switch route structure and operating modes, was used to solve real-world instances. Doppstadt et al. ( 41 ) extend Doppstadt et al.’s model ( 40 ) by considering customer time windows and proposing a new VNS-based solution method. Liao et al. ( 42 ) provided an efficient dynamic-programming-based polynomial-time algorithm for the EV shortest-travel-time path problem and approximation algorithms for the PHEV touring problems. The algorithms incorporated battery capacity constraints and battery swaps. Roberti and Wen ( 43 ) developed a MILP formulation for the EV TSP problem with time windows for both full- and partial-recharge policies. While there has been a significant amount of work on E-VRP and TSP, except for Zhu et al. ( 14 ), this has yet to consider an integrated delivery system with drones.

Meanwhile, many studies have investigated the efficiency of delivery systems that deploy UAVs. Otto et al. ( 44 ) provides a detailed review of the various civil applications of drones in domains such as agriculture, monitoring, transport, security, and so forth. Murray and Chu ( 10 ) introduced the flying-sidekick TSP (FSTSP), which assumes that a truck can launch its UAV at the depot or customer node and remain on its route. At the same time, the UAV delivers one small parcel to another customer before meeting again at a rendezvous location (another customer node on the truck’s route). Murray and Chu ( 10 ) proposed a two-stage route and reassigned heuristic wherein the first stage, a truck TSP tour that visits all customers, is determined. In the second stage, selected customers are reassigned to the UAV based on cost savings. Murray and Chu ( 10 ) also introduced the parallel-drone-scheduling TSP (PDSTSP), in which multiple drones and a truck originating from a depot serve a set of customers. Mbiadou Saleu et al. ( 45 ) developed an iterative two-stage heuristic involving customer partitioning and routing optimization for the PDSTSP. Agatz et al. ( 11 ) presented an integer programming formulation for a variant of FSTSP called TSPD and a “route-first, cluster-second” heuristic. Ha et al. ( 46 ) studied a variant of Murray and Chu ( 10 ) to minimize operating- and waiting-time costs rather than completion time. The authors propose two heuristics—a modification of Murray and Chu’s ( 10 ) heuristic to minimize costs and greedy randomized adaptive-search procedure (GRASP). Yurek and Ozmutlu ( 47 ) developed an iterative two-stage algorithm to solve the FSTSP of Murray and Chu ( 10 ), which was referred to as TSPD. The truck and drone routes are determined sequentially in two stages. Bouman et al. ( 48 ) modified the Bellman–Held–Karp dynamic-programming algorithm for the TSP to develop an exact solution approach for TSPD, whereas Poikonen et al. ( 49 ) used a branch-and-bound method. De Freitas and Penna ( 50 , 51 ) developed a randomized variable neighborhood descent heuristic, which modified an initial TSP solution obtained from the Concorde solver to solve the FSTSP. Boysen et al. ( 52 ) focused on scheduling single and multiple drone deliveries launched from a truck with a fixed route. Jeong et al. ( 53 ) modified Murray and Chu’s model ( 10 ) to include the impact of payload on energy consumption and no-fly zones and proposed a two-stage construction and search heuristic. Dayarian et al. ( 54 ) focused on a new variant in which drones are used to resupply a truck making deliveries. Kim and Moon ( 55 ) studied a variant termed traveling-salesman problem with drone station (TSPDS), in which a truck is used to resupply a drone station, which is different from a depot. Multiple drones then make deliveries to customers from drone stations. The truck will also make deliveries to customers after supplying the drone station. A two-phase solution algorithm is developed to solve the problem. Roberti and Ruthmair ( 56 ) proposed a compact mixed-integer programming formulation for TSPD and discussed how this formulation could be extended to address some common additional side constraints. In addition, it presented an exact branch-and-price algorithm in which the subproblem is solved using a dynamic-programming approach. More recently, Freitas et al. ( 57 ) proposed another mixed-integer programming formulation and general VNS metaheuristic for the same problem.

While there has been a significant body of work on integrating drones into existing routing frameworks since 2015, none consider EVs and their associated range constraints. Recently, the authors studied a similar problem in the context of battery EVs and drones ( 14 ). In Zhu et al. ( 14 ), we assume batteries are swapped and do not consider partial charging. Moreover, the ALNS heuristic developed in this current research outperforms the VNS metaheuristic developed in Zhu et al. ( 14 ). In addition, this paper also assumes a nonlinear relationship between the amount of energy charged and charging time, which significantly increases the problem’s difficulty but also indicates that the proposed model has great potential to be used in a real-world application.

Problem Description and Formulations

In the PHEVTSPD, a logistics company aims to deliver a specified number of packages to customers using a single PHEV paired with a commercial drone. In this context, The PHEV is assumed to be a small electric van specifically designed for delivery tasks, with a battery size similar to the van built by Alke, an Italian company), specifically designed for delivery tasks. Each customer is to be served exactly once by either the PHEV or the drone. The objective is to determine a coordinated route for the PHEV and the drone that minimizes the total route completion time.

Additional operational and energy assumptions for the proposed problem are summarized below:

Initially, both the PHEV and the drone are located at the depot, with the PHEV fully charged.

Both the PHEV and the drone travel at constant speeds within the network, and altitude differences are ignored.

Similar to the TSPD, the PHEV acts as a drone hub, launching the drone at customer node $i$ and retrieving it at customer node $k$ while the drone independently serves customer $j$ . Both nodes $i$ and $k$ must also be visited by the PHEV. This $i - j - k$ sequence is referred to as a drone sortie. This operation must satisfy the drone’s loading capacity, flight range, and energy requirements. Each drone launch or retrieval incurs a fixed time cost.

The proposed problem does not explicitly model customer service time, for simplicity, although service time can be incorporated as arc-traversing time.

As the PHEV traverses an arc in the network, its energy decreases at a known constant rate. The PHEV can recharge its battery at charging stations during operations. The amount of energy recharged and the charging time are decision variables.

The PHEV’s battery can also recharge the drone when it is onboard, implying that drone charging time is not separately considered. This is mainly because the energy consumption of the drone sortie is known before it is launched, so the drone does not need to be fully charged as long as it has enough energy to finish the sortie. With the drone sortie, the remaining energy of the PHEV’s battery decreases proportionally. The ratio of decreased PHEV battery energy and drone consumed energy is denoted $γ : γ > 0$ . This assumption is called the shared energy assumption.

Figure 1 illustrates a coordinated PHEV–drone route. Compared with the TSPD or the FSTSP with Drones (FSTSPD), the proposed problem follows a similar operational pattern but substitutes a conventional truck with a PHEV. A key difference is that the PHEV can recharge its battery at charging stations. Furthermore, the shared energy assumption implies that the PHEV and drone must coordinate not only operationally but also for energy management. The authors believe this assumption addresses the practical issue of the drone’s limited power supply and enhances the realism of the proposed problem scenario.

Figure 1.

A plug-in hybrid electric vehicle (PHEV)–drone coordinated route.

For the rest of the section, we introduce a MILP model and then explain how to consider partial charge and incorporate nonlinear time–SoC functions.

MILP Model for PHEVTSPD

The proposed model extends the formulation proposed in Murray and Chu ( 10 ). The original network is a directed, connected graph $G = (V, A)$ . The node set $V$ consists of the customer set $C = {1, 2, . . ., c}$ , the depot set ${0, 0'}$ , and charging-station (CS) node set $S = {c + 1, c + 2, . . ., c + s}$ . Thus $V = {0, 0}^{'} \cup C \cup S$ and $| V | = c + s + 1$ . In the arc set $A = {(i, j) : i, j \in V}$ , each arc $(i, j)$ is associated with two non-negative travel time costs $c_{ij}^{T}$ and $c_{ij}^{D}$ , which corresponds to the travel time needed for the PHEV and drone to travel from node $i$ to node $j$ , respectively. In addition, each arc is also associated with two non-negative energy costs $e_{ij}^{T}$ and $e_{ij}^{D}$ for the PHEV and the drone, respectively. Both $e_{ij}^{T}$ and $e_{ij}^{D}$ are also measured in time units.

To permit multiple visits to the charging-station nodes while requiring exactly one visit to the customer nodes, the original network $G$ is augmented into $m$ copies to create $G^{a} = (N, E^{a})$ . $m$ should be set as small as possible to reduce the network size but large enough not to restrict multiple beneficial visits. In Erdoğan and Miller-Hooks ( 29 ), this value is set to 2 or 3.

The sets used in the model are:

$C$ : Set of all customers in the problem, $C = {1, 2, . . ., c}$ and $| C | = c$ .

$C^{'}$ : Subset of customers that are available to drone delivery service, $C' \subset C$ .

$S$ : Set of all charging stations in the original network, $S = {c + 1, c + 2, . . ., c + s}$ and $| S | = s$ .

$S^{a}$ : Augmented set of all charging stations which includes a total of $m$ copies of set $S$ . $S^{a} = {c + 1, . ., c + s, c + s + 1, . . ., c + ms}$ and $| S^{a} | = ms$ . Note that in the augmented network, each node can be visited at most once.

$N$ : Set of all nodes in the augmented network, $N = S' \cup C \cup {0, 0'}$ , where $0$ and $0'$ both represent the depot and $| N | = c + ms + 2$ .

$N_{d}$ : Set of nodes from which a vehicle may depart in the augmented network. $N_{d} = {0, 1, . . .,} c \cup S'$

$N_{a}$ : Set of nodes to which a vehicle may arrive in the augmented network. $N_{a} = {1, 2, . . . c} \cup S' \cup {v_{c + ms + 1}}$

$N'$ : Set of all customer nodes and charging-station nodes in the augmented network. $N' = C \cup S^{a}$

$A$ : Set of all feasible arcs in the augmented network. $A = A_{t} \cup A_{d}$

$A_{d}$ : Set of all arcs in the augmented network that can be traversed by drone. $A_{d} = {(i, j) : i, j \in V, e_{ij}^{D} \leq Q^{D}}$

$A_{t}$ : Set of all arcs in the augmented network that can be traversed by PHEV. $A_{t} = {(i, j) : i, j \in V, e_{ij}^{T} \leq Q^{T}}$

$D$ : Set of drone’s feasible sortie $〈 i, j, k 〉$ , where the drone is launched from node $i$ , travels to node $j$ and returns to node $k$ . $D = {〈 i, j, k 〉 : i \in N_{d}, j \in C^{'}, k \in N_{a}, i \neq j, j \neq k, i \neq k, e_{ij}^{D} + e_{jk}^{D} \leq Q^{D}}$ .

A summary of parameters is given below:

$c_{ij}^{T}$ : Time cost of PHEV when traversing arc $(i, j) \in A_{t}$ .

$c_{ij}^{D}$ : Time cost for the drone when traversing arc $(i, j) \in A_{d}$ .

$d_{ijk}$ : Travel time cost for the drone to accomplish drone sortie $〈 i, j, k 〉 \in D$ .

$s_{i}$ : Service time for customer $i \in C$ .

$e_{ij}^{T}$ : Energy cost for the PHEV when traversing arc $(i, j) \in A_{t}$ .

$e_{ij}^{D}$ : Energy cost for the drone when traversing arc $(i, j) \in A_{d}$ .

$e_{ijk}^{D}$ : Energy cost for the drone to accomplish sortie $〈 i, j, k 〉 \in D$ .

$Q^{D}$ : Flight-energy capacity of the drone.

$Q^{T}$ : Driving-energy capacity of the PHEV.

$γ$ : Drone–PHEV energy consumption ratio, $γ > 0$ .

$M$ : A positive large number.

The decision variables are listed below:

$x_{ij}$ : Equals 1 if the PHEV travels arc $(i, j) \in A_{t}$ and 0 otherwise, where $i \neq j$ and $i \in N_{d}, j \in N_{a}$ .

$y_{ijk}$ : Equals 1 if the drone performs sortie $〈 i, j, k 〉 \in D$ and 0 otherwise.

$p_{ij}$ : Equals 1 if customer node $i$ is visited before customer $j$ in the PHEV’s path and 0 otherwise.

$u_{i}$ : Position of node $i$ in the PHEV’s path.

$b_{i}^{a}$ : Remaining battery charge of the PHEV on arrival at node $i$ .

$b_{i}^{d 1}$ : Remaining battery charge of the PHEV on departure from node $i$ before launching drone.

$b_{i}^{d 2}$ : Remaining battery charge of the PHEV on departure from node $i$ after launching drone.

$t_{j}^{a}$ : Time when the PHEV arrives at node $j$ .

$t_{j}^{d}$ : Time when the PHEV departs from node $j$ .

$t'_{j}$ : Time when the drone arrives at node $j$ .

Objective:

\min t_{c + m s + 1}

(1)

Routing constraints:

\begin{array}{l} \sum_{\begin{matrix} i \in N_{d} \\ i \neq j \end{matrix}} x_{i j} + \sum_{\begin{matrix} i \in N_{d} \\ i \neq j \end{matrix}} \sum_{\begin{matrix} k \in N_{a} \\ 〈 i, j, k 〉 \in D \end{matrix}} y_{i j k} = 1 \\ \forall j \in C \end{array}

(2)

\sum_{j \in N_{a}} x_{0 j} = 1

(3)

\sum_{i \in N_{d}} x_{i, 0'} = 1

(4)

\begin{matrix} u_{i} - u_{j} + 1 \leq (c + ms + 2) (1 - x_{ij}) \\ \forall i \in N', j \in N_{a}, j \neq i \end{matrix}

(5)

\begin{matrix} \sum_{\binom{i \in N_{d}}{i \neq j}} x_{ij} = \sum_{\binom{k \in N_{a}}{k \neq j}} x_{jk} \end{matrix}

(6)

\begin{matrix} \sum_{j \in C'} \sum_{\binom{k \in N_{a}}{〈 i, j, k 〉 \in D}} y_{ijk} \leq 1 \\ \forall j \in N' \\ \forall i \in N_{d} \end{matrix}

(7)

\begin{matrix} \sum_{i \in N_{d}} \sum_{\binom{j \in C'}{〈 i, j, k 〉 \in D}} y_{ijk} \leq 1 \\ \forall k \in N_{a} \end{matrix}

(8)

\begin{matrix} 2 y_{ijk} \leq \sum_{\binom{h \in N_{d}}{h \neq i}} x_{hi} + \sum_{\binom{l \in N_{d}}{l \neq k}} x_{lk} \\ \forall i \in N_{d}, j \in C', k \in N_{a}, 〈 i, j, k 〉 \in D \end{matrix}

(9)

\begin{matrix} y_{0 jk} \leq \sum_{\binom{h \in N_{d}}{h \neq k}} x_{hk} \\ \forall j \in C^{'}, k \in N_{a}, 〈 0, j, k 〉 \in D \end{matrix}

(10)

\begin{matrix} 1 - (c + ms + 2) (1 - \sum_{\binom{j \in C^{'}}{〈 i, j, k 〉 \in D}} y_{ijk}) \leq u_{k} - u_{i} \\ \forall i \in N_{d}, k \in N_{a}, k \neq i \end{matrix}

(11)

The objective (1) is to minimize the operational time. Constraints 2 to 11 are associated with the routing of the two vehicles. Constraint 2 guarantees that each customer node is visited once by either the PHEV or the drone. Constraint 3 and 4 state that the PHEV must start from and return to the depot. Constraint 5 is a sub-tour elimination constraint for the PHEV. Constraint 6 indicates that if the PHEV visits node $j$ , it must depart from node $j$ . Constraints 7 and 8 state that each node can launch or retrieve the drone at most once. Constraint 9 ensures that if the drone sortie $〈 i, j, k 〉$ is executed, then the PHEV should visit node $i$ and $k$ . Constraint 10 states that if the drone is launched from the depot and returned to node $k$ , then node $k$ should be visited by the PHEV. Constraint 11 shows sub-tour elimination constraints for the drone.

Energy constraints:

\begin{matrix} b_{j}^{a} \leq (b_{i}^{d 2} - e_{ij}^{T}) x_{ij} + M (1 - x_{ij}) \\ \forall i \in N_{d}, j \in N_{a}, i \neq j \end{matrix}

(12)

b_{0}^{a} = Q^{T}

(13)

b_{i}^{a} = b_{i}^{d 1} \forall i \in C_{0}

(14)

b_{i}^{a} \leq b_{i}^{d 1} \forall i \in S

(15)

b_{i}^{d 2} = b_{i}^{d 1} - γ \times \sum_{\binom{j \in C^{'}}{j \neq i}} \sum_{\binom{k \in N_{a}}{〈 i, j, k 〉 \in D}} y_{ijk} e_{ijk}^{D} \forall i \in N_{d}

(16)

0 \leq b_{i}^{a}, b_{i}^{d 1}, b_{i}^{d 2} \leq Q^{T} \forall i \in N

(17)

Constraints 12 to 17 are associated with the battery energy level. In particular, Constraint 12 states that if the PHEV traverses arc $(i, j)$ , then the energy level before arriving at node $j$ is $e_{ij}^{T}$ less than the energy level on leaving node $i$ . Constraint 13 ensures that the PHEV is fully charged when it departs from the depot. Constraint 14 indicate that if the PHEV visits a customer node $i$ , the remaining energy on departure at node $i$ is the same as the energy on arrival at node $i$ . Similarly, Constraint 15 enforces the principle that if the PHEV visits a charging station $i$ , then the remaining energy on departure at node $i$ is no less than the energy on arrival at node $i$ . Constraint 16 models the energy change at each node $i \in N_{d}$ because of the launch of the drone. Constraint 17 gives the domain constraints for $b_{i}^{a}, b_{i}^{d 1}, b_{i}^{d 2}$ .

Timing constraints:

t'_{i} \geq t_{i}^{d} - M (1 - \sum_{\binom{j \in C^{'}}{j \neq i}} \sum_{\binom{k \in N_{a}}{〈 i, j, k 〉 \in D}} y_{ijk}) \forall i \in N_{d}

(18)

t'_{i} \leq t_{i}^{d} + M (1 - \sum_{\binom{j \in C^{'}}{j \neq i}} \sum_{\binom{k \in N_{a}}{〈 i, j, k 〉 \in D}} y_{ijk}) \forall i \in N_{d}

(19)

t'_{k} \geq t_{k}^{a} - M (1 - \sum_{\binom{i \in N_{d}}{i \neq k}} \sum_{\binom{j \in C'}{〈 i, j, k 〉 \in D}} y_{ijk}) \forall k \in N_{a}

(20)

t'_{k} \leq t_{k}^{a} + M (1 - \sum_{\binom{i \in N_{d}}{i \neq k}} \sum_{\binom{j \in C^{'}}{〈 i, j, k 〉 \in D}} y_{ijk}) \forall k \in N_{a}

(21)

\begin{matrix} t_{k}^{a} \geq t_{h}^{d} + c_{hk}^{T} - M (1 - x_{hk}) \\ \forall h \in N_{d}, k \in N_{a}, k \neq h \end{matrix}

(22)

t_{j}^{d} \geq t_{j}^{a} + s_{i} \forall j \in C_{0}

(23)

t_{j}^{d} \geq t_{j}^{a} + f (b_{j}^{d 1}, b_{j}^{a}) \forall j \in S

(24)

\begin{matrix} {t'}_{j} \geq {t'}_{j} + c_{ij}^{D} + s_{j} - M (1 - \sum_{\binom{k \in N_{a}}{〈 i, j, k 〉 \in D}} y_{ijk}) \\ \forall j \in C', i \in N_{d}, i \neq j \end{matrix}

(25)

\begin{matrix} {t'}_{k} \geq {t'}_{j} + c_{jk}^{D} - M (1 - \sum_{\binom{i \in N_{d}}{〈 i, j, k 〉 \in D}} y_{ijk}) \\ \forall j \in C^{'}, k \in N_{a}, k \neq j \end{matrix}

(26)

Constraints 18 to 26 are associated with the travel time of the two vehicles. In particular, Constraints 18 to 21 ensure that the travel time and drone range limit are correctly handled. To be more specific, if drone sortie $〈 i, j, k 〉$ is executed, then $t'_{i} = t_{i}^{d}$ and $t'_{k} = t_{k}^{a}$ . Constraint 22 indicates that if the PHEV travels from node $h$ to node $k$ where $h \in N_{d}, k \in N_{a}$ , its arrival time at node $k$ must incorporate its arrival time at node $h$ , travel time from node $h$ to node $k$ , the drone’s launch time at node $h$ , and retrieve time at node $k$ . This constraint is not binding if the EV does not travel from node $h$ to node $k$ . Constraints 25 and 26 are associated with the drone’s arrival time. Suppose there is a drone route of $〈 i, j, k 〉$ , then the drone’s arrival time of node $j$ and node $k$ should be related to the drone’s travel time between $i$ and $j$ , and $j$ and $k$ , and the drone’s retrieve time at node $k$ .

Ordering constraints:

\begin{matrix} e_{ij}^{D} + e_{jk}^{D} \leq Q^{D} + M (1 - y_{ijk}) \\ \forall k \in N_{a}, j \in C^{'}, i \in N_{d}, 〈 i, j, k 〉 \in D \end{matrix}

(27)

\begin{matrix} u_{i} - u_{j} \leq - 1 + (c + ms + 2) (1 - p_{ij}) \\ \forall i, j \in N', i \neq j \end{matrix}

(28)

\begin{matrix} p_{ij} + p_{ji} = 1 \\ \forall i, j \in N^{'}, i \neq j \end{matrix}

(29)

\begin{matrix} {t'}_{l} - {t'}_{k} \geq - M (3 - \sum_{\binom{j \in C'}{\binom{〈 i, j, k 〉 \in D}{j \neq l}}} y_{ijk} - \sum_{\begin{matrix} m \in C' \\ m \neq i \\ m \neq k \\ m \neq l \end{matrix}} \sum_{\begin{matrix} n \in N_{a} \\ 〈 l, m, n 〉 \in D \\ n \neq i \\ n \neq k \end{matrix}} y_{lmn} - p_{il}) \\ \forall i, l \in N_{d}, k \in N_{a}, i \neq k \neq l, 〈 i, j, k 〉, 〈 l, m, n 〉 \in D \end{matrix}

(30)

Constraints 27 to 30 are associated with ordering the two vehicles. Constraint 27 ensures that the drone route should be within the drone’s flight range. Constraint 28 is a sub-tour elimination constraint, and Constraint 29 ensures the correct ordering of two different nodes. Constraint 30 indicates that if there exist two drone route deliveries $〈 i, j, k 〉$ and $〈 l, m, n 〉$ and node $i$ is visited before node $l$ by the PHEV, then node $l$ must be visited after node $k$ .

Domain constraints:

t_{0} = 0

(31)

p_{0 j} = 1 \forall j \in N_{a}

(32)

x_{ij} \in {0, 1} \forall i \in N_{d}

(33)

y_{ijk} \in {0, 1} \forall i \in N_{d}

(34)

1 \leq u_{i} \leq c + ms + 2 \forall i \in N_{a}

(35)

t_{i} \geq 0 \forall i \in N

(36)

t'_{i} \geq 0 \forall i \in N

(37)

p_{ij} \in {0, 1} \forall i \in N_{d}, j \in N_{a}, j \neq i

(38)

The proposed formulation is also capable of accommodating common side constraints, which are discussed in detail in the Appendix. We acknowledge the existence of various models addressing the TSPD, which are similar to PHEVTSPD. Examples include the two-index formulation originally proposed by Roberti and Ruthmair ( 56 ), the formulation by Luo et al. ( 58 ), and the one proposed in Freitas et al. ( 57 ). Although these MILP models are efficient, they do not apply to our problem because these models do not model the energy balance at the charging-station nodes explicitly since there are no charging-station nodes involved.

Modeling Partial-Recharge Policies

In the literature related to EV routing, a full recharge assumption is commonly adopted because of its simplicity—for example, Schneider et al. ( 31 ). This assumption is also valid when the EV is assumed to be a battery-swappable EV, where its battery can be exchanged for a fully charged one at battery swapping stations ( 14 , 59 ).

However, full recharge can be time-consuming, depending on the SoC level, available charging technology, and battery capacity. In time-sensitive routing scenarios, such as logistics companies aiming to deliver parcels within specific time constraints, a full-charge policy might yield suboptimal solutions. For instance, if the PHEV visits a charging station near the end of its route, it does not need a full charge to return to the depot. This situation also arises when the PHEV visits consecutive charging stations, where full charging is unnecessary for the journey to the next station. Additionally, in respect of vehicle maintenance, deep charge/discharge cycles are detrimental to the lifespan of lithium batteries, which are commonly used in EVs, and significantly accelerate battery degradation ( 60 ). Research indicates that a complete charge/discharge cycle hastens battery capacity fading and reduces the number of dynamic-stress-test cycles, as shown in Figure 2.

Figure 2.

Battery capacity for dynamic-stress-test cycles.

Given these drawbacks of full recharge, partial recharge, where the battery is charged only enough to complete the route, is often preferable in real-world applications. Notable research adopting this assumption includes Montoya et al. ( 15 ), Desaulniers et al. ( 61 ), Bruglieri et al. ( 62 ), and Froger et al. ( 63 ). Surveys on the E-VRP, its variants, and different charging strategies are provided by Erdelić and Carić ( 39 ) and more recently by Abid et al. ( 64 ).

In the proposed problem, we adopt a partial-recharge assumption with a nonlinear time–SoC function. We assume that the charging stations are connected to a plug-in electric grid, and the PHEV can be charged with any amount of energy at CS nodes. Both the amount of energy charged and the charging time, which depends on the initial SoC level and the amount of energy charged, are decision variables.

In the plug-in charging stations, the actual plug-in time–SoC curve is a nonlinear concave function, as shown in Figure 3. In the figure, $u$ and $i$ represent terminal voltage and current, respectively. The time–SoC curve is nonlinear because the $i$ and $u$ change during the charging process. In the first stage, SoC increases linearly with charging time as the current $i$ is constant. After SoC reaches a specific value, the current $i$ decreases exponentially while the voltage $u$ is kept constant to avoid damage to the battery. As a result, the SoC grows concavely with time in the second phase. This charging curve indicates that charging efficiency decreases drastically when the SoC level is high.

Figure 3.

A typical charging curve.

In the field of EV routing that adopts a partial-recharge policy, a common strategy assumes a linear relationship between charging time and the charged amount. Research such as that in ( 35 , 37 , 62 , 66 , 67 ) often uses a conservative estimate of the actual time–SoC function, ensuring that for any charging time $t$ , the actual SoC level after charging is higher than the linear approximation to maintain solution feasibility.

One major drawback of the linear approximation of the time–SoC function is its accuracy. The relationship between the amount of charged energy and the charging time is not linear; it depends on the initial time–SoC level before charging, as depicted in Figure 4. Montoya et al. ( 15 ) were the first to model the concave time–SoC function using a piecewise linear function. This approach is more precise than the linear counterpart, as shown in Figure 5, where $| R |$ represents the number of segments in the piecewise linear function. Montoya et al. ( 15 ) considered four additional categories of decision variables in the MILP formulation. The total cost when using the nonlinear charging function was found to be 2.70% cheaper than that of the linear function on average. Zuo et al. ( 68 ) further simplified this piecewise linear approximation, requiring only one additional decision variable. This research adopts the same technique as Zuo et al. ( 68 ) to model the piecewise linear time–SoC function.

Figure 4.

Relationship between charged time and charged energy amount.

Figure 5.

Actual time–SoC function and piecewise linear function approximation with $| R | = 4$ .

To deploy this technique, some necessary changes in the notations are:

$e_{ij}^{T} \in [0, 1]$ : Energy cost relative to PHEV’s energy capacity when traversing arc $(i, j)$ .

$e_{ijk}^{D}$ : Energy cost relative to PHEV’s energy capacity when performing sortie $(i, j, k)$ .

$b_{i}^{a} \in [0, 1]$ : Remaining energy relative to PHEV’s energy capacity on arriving at node $i$ .

$b_{i}^{d 1} \in [0, 1]$ : Remaining energy relative to PHEV’s energy capacity on departure at node $i$ while before launch operation.

$b_{i}^{d 2} \in [0, 1]$ : Remaining energy relative to PHEV’s energy capacity on departure at node $i$ after launch operation.

For the new notation, all the energy values are now represented as their ratio to the PHEV battery energy capacity. Let $f (t)$ denote the concave energy-charging function at charging time $t$ . Let $K_{r}$ and $B_{r}$ denote the slope and intercept of each secant line $r \in R$ . Let $S_{i}$ represent the energy-charging level and $t_{i}$ the required charging time for the battery to be charged from $SoC = 0 %$ to $SoC = S_{i}$ at charging station $i$ . Thus, it is obvious that

S_{i} \leq K_{r} t_{i} + B_{r}, \forall r \in R

(39)

Let $(t_{i}, s_{i}), i = 1, 2, . . ., k + 1$ be the breakpoints of the secant lines of the SoC charging curve. Since there are a total of $| R | = k$ secant lines, there are $k + 1$ breakpoints. Before charging, the current SoC level at charging-station node $i \in S$ is $b_{i}^{a}$ . Denote $t_{i}^{ba}$ as the time needed for the battery to be charged from $SoC = 0 %$ to $SoC = b_{i}^{a}$ . Let $Δ S$ be the amount of energy that needs to be charged and $Δ t$ be the time required to charge the battery from $SoC = b_{i}^{a}$ to $SoC = b_{i}^{a} + Δ S$ . Then, with the same principle of Equation 39, we can derive the following equation:

b_{i}^{a} + Δ S \leq K_{r} (t_{i}^{ba} + Δ t) + B_{r}, \forall r \in R

(40)

Note that in our MILP formulation $Δ t = t_{i}^{d} - t_{i}^{a}$ and $Δ S = b_{i}^{d 1} - b_{i}^{a}$ . So the above equation is equivalent to

b_{i}^{d 1} \leq K_{r} (t_{i}^{ba} + t_{i}^{d} - t_{i}^{a}) + B_{r}, \forall r \in R, i \in S

(41)

Equation 41 models the relationship between $b_{i}^{d 1}$ , $t^{d}$ , and $t^{a}$ . To linearize it, the following equations are needed. Let $α_{ri}$ represent a binary variable that equals 1 if a particular secant line $r$ is used to represent the relation between $b_{i}^{a}$ and $t_{i}^{ba}$ . Thus,

b_{i}^{a} \geq s_{r} + M (α_{ri} - 1), \forall r \in R, i \in S

(42)

b_{i}^{a} \leq s_{r + 1} + M (1 - α_{ri}), \forall r \in R, i \in S

(43)

\sum_{r \in R} α_{ri} = 1, \forall i \in S

(44)

where Constraints 42 and 43 ensure that if $α_{r} = 1$ then $b_{i}^{a}$ is bounded by $s_{r}$ and $s_{r + 1}$ and Constraint 44 ensures that only one secant line is selected. Now, we establish the relationship between existing variables with $t^{ba}$ , which can be expressed as

b_{i}^{a} \geq K_{r} t_{i}^{ba} + B_{r} - M (1 - α_{ri}), \forall r \in R, i \in S

(45)

b_{i}^{a} \leq K_{r} t_{i}^{ba} + B_{r} + M (1 - α_{ri}), \forall r \in R, i \in S

(46)

With a new reformulation and the introduction of new auxiliary variables $α_{ri}, t_{i}^{ba}, \forall i \in S, r \in R$ , Constraint 24 is dropped and the complete formulation for PHEVTSPD is Objective ( 1 ) with constraints ( 2 ) to ( 23 ), constraints ( 25 ) to ( 39 ), and constraints ( 42 ) to ( 47 ).

ALNS Solution Method

This section proposes a metaheuristic solution method for PHEVTSPD. In the past, multiple researchers have shown that VNS produces some of the best results when solving TSPD ( 51 , 69 ). However, when applying the neighborhood structures to a feasible solution $s$ , the resulting solution $s^{t}$ might be infeasible as a result of the PHEV’s driving-range-limit constraint. Thus, we propose an alternative metaheuristic, ALNS, first proposed in Shaw ( 70 ) and Ropke and Pisinger ( 71 ). Compared with VNS, the number of iterations and newly found solutions in Large Neighborhood Search (LNS) is much smaller as it adopts more costly methods for obtaining new solutions.

Overview of ALNS

Typically, in LNS, two kinds of operators exist: destroy methods and repair methods. The destroy method aims to “destroy” the current solution by eliminating a target number of nodes and the repair method tries to restore a feasible solution. In ALNS, both methods are chosen based on their performance score, with the roulette-wheel selection rule. Denote $A$ as the set that contains all the available destroy methods and $B$ as the repair-method set. At the start of iteration $j$ , the performance score of method $i$ is denoted as $w_{ij}$ , then the probability of method $i$ being chosen is

p_{ij} = \frac{w_{ij}}{\sum_{i} w_{ij}}

(47)

After each iteration, $w_{ij}$ will be updated based on their performance in the iteration as follows:

w_{i, j + 1} = ρ w_{ij} + \frac{λ}{τ} (1 - ρ)

(48)

In the formula, $ρ$ is the score-shrinking rate, and $λ / τ$ is the score for the method in iteration $j$ , where $λ$ is an indication of the fitness of the newly found solution and $τ$ is an indication of the compute time of the current iteration. Intuitively, the value of $λ / τ$ is higher if we find a better solution with a lower cost or enter a search space that has not been visited before in a short amount of time.

We also control the acceptance probability of a new solution with a global time-varying parameter $T$ called temperature. This concept is frequently used in simulated annealing ( 72 ). In this research, if a better solution is found, it is always accepted. If the new solution has a higher cost, the probability of it being accepted is

p (x^{t}) = e^{\frac{f (s) - f (s^{t})}{T}}

(49)

where $f (s)$ denote the objective value of solution $s$ . As the solution space is explored, the value of $T$ gradually decreases from the initial value to zero. Given a prespecified time limit for the search phase, the value of $T$ is updated as

T = T (1 - \frac{t}{t_{m}})

(50)

where $t$ is the elapsed search time and $t_{m}$ is the preset time limit. To avoid the value of the denominator decreasing to zero in the formula (50), the search terminates when the value of $T$ is sufficiently small. In this research, the initial temperature is set to be 1000. The pseudo-code for the ALNS algorithm is given in

Algorithm 1 ALNS
Input: PHEVTSPD instance and preset parameters Output: EVTSP-D solution 1: $s \leftarrow InitialSolution ()$ 2: $s^{} \leftarrow s$ 3: While $t \leq t_{m}$ do 4: Choose a pair of destroy and repair method $(d (s), r (s))$ based on their performance score 5: Get a temporary solution $s^{temp} \leftarrow r (d (s))$ 6: Get new feasible solution $s^{t}$ . If $s^{temp}$ is feasible, then $s^{t} \leftarrow s^{temp}$ . Otherwise, $s^{t} \leftarrow g (s^{temp})$ where $g (s)$ is the feasibility restoring function 7: $T = T_{0} (1 - t / t_{m})$ 8: if $random (0, 1) < e^{\frac{f (s) - f (s^{t})}{T}}$ then 9: $s \leftarrow s^{t}$ 10: end If 11: Update performance score for destroy/repair methods 12: end while 13: return $s^{}$

Algorithm 1 ALNS

Input: PHEVTSPD instance and preset parameters
Output: EVTSP-D solution
1:

s \leftarrow InitialSolution ()

s^{*} \leftarrow s

3: While

t \leq t_{m}

do
4: Choose a pair of destroy and repair method

(d (s), r (s))

based on their performance score
5: Get a temporary solution

s^{temp} \leftarrow r (d (s))

6: Get new feasible solution

s^{t}

. If

s^{temp}

is feasible, then

s^{t} \leftarrow s^{temp}

. Otherwise,

s^{t} \leftarrow g (s^{temp})

where

g (s)

is the feasibility restoring function
7:

T = T_{0} (1 - t / t_{m})

8: if

random (0, 1) < e^{\frac{f (s) - f (s^{t})}{T}}

then
9:

s \leftarrow s^{t}

10: end If
11: Update performance score for destroy/repair methods
12: end while
13: return

s^{*}

Initial Solution

The modified Clarke–Wright savings algorithm (MCWS), a greedy maximum-saving algorithm proposed in Erdoğan and Miller-Hooks ( 29 ), is used to construct a PHEV route that satisfies the battery constraint. The resulting initial solution has all the customers served by the PHEV and has no drone sorties.

Removal Methods

In each iteration, the removal operator removes some customer nodes from the current solution $s$ and returns a partial solution $s^{p}$ and a node set $P$ . In this research, two different destruction methods are considered. The first one is random removal, which randomly removes a target number of customer nodes from the current solution until a target number is reached. The second is cluster removal, which removes a random customer node and progressively removes its nearby nodes until a target number is reached.

Repair Methods

In each iteration, with a given partial solution $s^{p}$ and node set $P$ , the repair method aims to insert nodes in $P$ into $s^{p}$ and return a new feasible solution $s^{t}$ . In general, the repair method is a two-phase modification process. It first tries to construct a feasible solution and then improves it. This paper considers three different repair methods, two of which are inspired by Sacramento et al. ( 73 ).

Greedy Repair

The greedy repair method is the “cheapest” of all three repair methods. In the first phase, the feasible routes are constructed by inserting all the nodes in $P$ into the truck’s routes greedily. Later, the method improves the solution by deploying the drone as much as possible. The procedure is shown in Algorithm 2.

Algorithm 2 Greedy repair
Input: partial solution $s^{p}$ and picked node set $P$ Output: new feasible solution of $s^{t}$ 1: node $i$ in $P$ 2: Find the insert location with least-extra-cost; 3: If the resulting route is feasible, continue. Otherwise, insert CS to restore its feasibility. 4: end for 5: for node $i$ in $r^{T}$ _do 6: if $i$ could be served by drone with the previous node as launch node and next node as retrieve node then 7: Add sortie $s$ to the current route; 8: Delete $i$ from the truck’s route. 9: end if 10: end for 11: return $s^{p}$

Algorithm 2 Greedy repair

Input: partial solution

s^{p}

and picked node set

P

Output: new feasible solution of

s^{t}

1: node

i

P

2: Find the insert location with least-extra-cost;
3: If the resulting route is feasible, continue. Otherwise,
insert CS to restore its feasibility.
4: end for
5: for node

i

r^{T}

_do
6: if

i

could be served by drone with the previous node as launch node and next node as retrieve node then
7: Add sortie

s

to the current route;
8: Delete

i

from the truck’s route.
9: end if
10: end for
11: return

s^{p}

Nearby Repair

The first phase of the nearby repair is similar to that of the greedy repair. In the second phase, the nearby repair aims to explore all possible drone sorties (by exploring different launch and retrieve nodes) and pick the sortie with the largest time-savings. This process is repeated until there is no potential drone sortie with positive time-saving. The detailed procedure of the nearby repair is shown in Algorithm 3.

Algorithm 3 Nearby repair
Input: partial solution $s^{p}$ and picked node set $P$ Output: new feasible solution of $s^{t}$ 1: node $i$ in $P$ 2: Insert node into $s^{p}$ before a customer node that is closest to the chosen node 3: end for 4: for each feasible segment $f$ in $r^{T}$ do 5: for node $i$ in $f$ do 6: Examine every feasible drone sortie that serves $i$ ; 7: Select sortie $s$ with the greatest saving; 8: If resulting route is feasible, add sortie $s$ into $s^{p}$ ; 9: Break $f$ into two remaining smaller segments and add them into segment set. 10: end for 11: end for 12: return $s^{p}$

Algorithm 3 Nearby repair

Input: partial solution

s^{p}

and picked node set

P

Output: new feasible solution of

s^{t}

1: node

i

P

2: Insert node into

s^{p}

before a customer node that is closest to the chosen node
3: end for
4: for each feasible segment

f

r^{T}

do
5: for node

i

f

do
6: Examine every feasible drone sortie that serves

i

;
7: Select sortie

s

with the greatest saving;
8: If resulting route is feasible, add sortie

s

into

s^{p}

;
9: Break

f

into two remaining smaller segments and add them into segment set.
10: end for
11: end for
12: return

s^{p}

CP Repair

Constraint programming is a paradigm for solving combinatorial problems by logically expressing the constraints. This method is specifically to derive a feasible solution and has been widely used in tackling schedule-related problems. In the community of vehicle routing, especially drone-related vehicle routing, CP has been proven to be an effective method for obtaining solutions of good quality—for example, Montemanni and Dell’Amico ( 74 , 75 ) and Ham ( 76 ). In this section, it is an essential part of the algorithm to get a feasible complete solution from the partial solution. The CP repair method is the most computationally demanding of all three techniques. It will eliminate all these existing drone sorties and re-insert these customer nodes into the partial solution, by solving a CP scheduling problem with the CPLEX CP optimizer. The CP repair method is described in Algorithm 4. A detailed description of the CP method and its equivalent MILP formulation, which is first proposed in Yurek and Ozmutlu ( 47 ) are shown in the Appendix ( 74 – 76 ).

Algorithm 4 CP repair
Input: partial solution $s^{p}$ and picked node set $P$ Output: an improved feasible solution of $s^{t}$ 1 : for node $i$ in $P$ do 2: insert node into $s^{p}$ with greedy method. 3: end for 4: $p' \leftarrow$ null set; 5: Delete all drone sorties in $r_{i}^{D}$ ,add drone nodes into $p'$ ; 6: Randomly select truck nodes in $r_{i}^{T}$ and add them into $p'$ ; 7: $r_{i}^{T'} \leftarrow$ remaining truck route; 8: Solve the CP model or MILP model with $r_{i}^{T'}$ and $p'$ . 9: return $s^{p}$

Algorithm 4 CP repair

Input: partial solution

s^{p}

and picked node set

P

Output: an improved feasible solution of

s^{t}

1 : for node

i

P

do
2: insert node into

s^{p}

with greedy method.
3: end for
4:

p' \leftarrow

null set;
5: Delete all drone sorties in

r_{i}^{D}

,add drone nodes into

p'

;
6: Randomly select truck nodes in

r_{i}^{T}

and add them into

p'

;
7:

r_{i}^{T'} \leftarrow

remaining truck route;
8: Solve the CP model or MILP model with

r_{i}^{T'}

and

p'

.
9: return

s^{p}

Feasibility-Resorting Function and Cost-Evaluation Function

At each iteration, after the removal methods and repair methods are applied to the current solution $s$ , a temporary (possibly infeasible) solution $s^{t}$ is returned. If $s^{t}$ is feasible, then its quality is ready to be evaluated. If $s^{t}$ is infeasible, the feasibility-resorting function $g (s)$ needs to be called. If $g (s)$ fails to restore the feasibility of $s^{t}$ , the solution is then discarded.

The feasibility-resorting function $g (s)$ aims to insert charging-station nodes into the EV’s route of $s^{t}$ to restore its feasibility. Note that this process will not affect the feasibility of drone sorties. If the PHEV’s driving range constraint is violated at arc $(n_{i}, n_{i + 1})$ , the feasibility-resorting function $g (s)$ need to be called:

Step 1: Identify the violating arc $(n_{i}, n_{i + 1})$ . If none is found, the solution $s$ is feasible. Otherwise, move to Step 2.

Step 2: Set current arc $a^{cur}$ as arc $(n_{i}, n_{i + 1})$ . Identify the CS node $s_{i}$ that is closest to node $n_{i}$ and insert it into $a^{cur}$ . If after inserting $s_{i}$ the PHEV could traverse the current arc successfully, return to Step 1. Otherwise, move to Step 3.

Step 3: Set current arc $a^{cur}$ as arc $(n_{i - 1}, n_{i})$ , which is before the current arc. Identify the CS node $s_{i - 1}$ that is closest to node $n_{i - 1}$ . Try to insert $s_{i - 1}$ OR $s_{i}$ into the current arc. If after inserting $s_{i}$ the PHEV could traverse the current arc successfully, return to Step 1. Otherwise, repeat Step 3, until a preset iteration limit is reached.

Computational Experiments

This section compares the performance of ALNS and VNS, and the computational efficiency and accuracy of ALNS are tested against a commercial solver. A sensitivity analysis of the number of line segments used in the time–SoC function approximation process is also performed.

Experimental Setting

Although publicly available instances of the TSPD exist, converting these instances into those suitable for PHEVTSPD proves challenging. As a result, numerical analyses are conducted on randomly generated instances. For each instance, the numbers of customers and charging stations are specified beforehand. Then, the coordinates of customers and charging stations are randomly generated within a 40 km × 40 km square, while the depot is always located at the center of the square. The distance matrix for the PHEV employs the Manhattan metric, while the Euclidean metric is applied to the drone, assuming that the drone can fly directly between two nodes in the network. Each instance is treated as an electric vehicle traveling-salesman problem (EVTSP) and the MCWS algorithm applied before it is chosen to ensure its feasibility. In this study, the PHEV is assumed to have a similar battery configuration as the Alke ATX340E, and the drone is the DJI MATRICE 600 PRO. A summary of configurations of the PHEV and drone is shown in Table 1; these are primarily derived from the companies’ official websites ( 77 , 78 ).

Table 1.

List of Default Parameters Used in Tests

Parameter	Value
PHEV driving speed	40 km/h
PHEV maximum driving range	100 km
Drone flying speed	60 km/h
Drone maximum flight range	20 km
Drone launch time	100 s
Drone retrieve time	29 s
Drone–PHEV energy consumption rate $γ$	0.4
Drone–PHEV speed ratio $β$	1.5

Note: PHEV = plug-in hybrid electric vehicle.

The proposed MILP formulation is implemented in Pyomo and solved using ILOG’s CPLEX Concert Technology solver (version 12.6.3). Simultaneously, the ALNS algorithm is coded in Python. The CP repair method is solved with the CPLEX CP optimizer (version 12.6.3). All experiments are conducted on a 3.6 GHz Intel Core i7 desktop with 32 GB RAM.

Estimation of the Concave Time–SoC Function

Multiple SoC charging functions of different batteries under different charging techniques exist. For example, a detailed time–SoC function provided by TeoTeslaFan in the Tesla forum is publicly available ( 79 ), which indicates that for the Tesla Model S, the battery could be fully charged in 110 min with a Tesla Supercharger. For our problem, the PHEV has a charging time of less than 90 min, according to Alke’s official website. Thus, necessary changes are made to the function posted in TeoTeslaFan (80). This concave SoC function is estimated with piecewise linear approximation with different segments, shown in Figure 6.

Figure 6.

Time–SoC function and linear/piecewise linear approximation with various $R$ values: (a) linear approximation, (b) piecewise linear approximately R = 2, (c) piecewise linear approximately R = 4, and (d) piecewise linear approximately R = 6.

Performance Comparison Between ALNS and VNS

The performance of the ALNS is compared with the VNS presented in Campuzano et al. ( 69 ). The test results are shown in Figure 7, where the displayed numbers are the average value of 20 randomly generated instances.

Figure 7.

Comparison between ALNS and VNS.

As can be seen, when the size of the tested instances is small $(| C | < = 10)$ , both heuristics have similar performance. However, as the instance size gets bigger $(| C | > 10)$ , ALNS continuously outperforms VNS, and the gap in the final cost is about 15%, which indicates that ALNS performs better than VNS for our problem.

Comparison Between MILP Formulation and ALNS Method

The performance of the ALNS is also compared with the CPLEX solver, whose results are shown in Table 2. In the table, “R1” represents the method approximating the time–SoC function with $| R | = 1$ while “R2” stands for results with $| R | = 2$ . In the last Gap columns, the optimality gap of ALNS is calculated as

\frac{t_{A} - t_{C}}{t_{C}} \times 100 %

(51)

where $t_{A}$ and $t_{C}$ represent the average completion time returned by ALNS and CPLEX.

Table 2.

Performance Comparison Between ALNS and CPLEX on Base Instances

					Opt				Running time				Gap (%)
					CPLEX		ALNS		CPLEX		ALNS		ALNS
N	n	$n_{cs}$	$β$	No. ins	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2
8	4	2	1.5	20	796.58	772.20	796.58	792.20	17.20	10.80	5	5	0.00	2.59
	4	2	2	20	743.25	741.24	743.25	741.24	17.20	10.80	5	5	0.00	0.00
	4	2	2.5	20	705.69	699.21	705.69	699.21	17.10	10.80	5	5	0.00	0.00
9	5	2	1.5	20	959.54	930.61	969.26	951.07	139.83	113.25	5	5	1.01	2.20
	5	2	2	20	902.21	867.30	907.65	870.78	140.41	114.21	5	5	0.60	0.40
	5	2	2.5	20	848.14	824.64	849.61	825.89	139.00	114.30	5	5	0.17	0.15
10	6	2	1.5	10	1156.80	1122.78	1197.47	1157.64	2259.33	2804.00	20	20	3.52	3.10
	6	2	2	10	1052.39	1019.85	1084.68	1040.02	2305.68	2814.27	20	20	3.07	1.98
	6	2	2.5	10	980.32	960.29	1000.26	990.95	2345.17	2845.67	20	20	2.03	3.19
14	8	3	1.5	10	na	na	1349.78	1295.78	3600	3600	20	20	na	na
	8	3	2	10	na	na	1186.20	1136.56	3600	3600	20	20	na	na
	8	3	2.5	10	na	na	921.53	902.60	3600	3600	20	20	na	na
20	10	5	1.5	10	na	na	1657.97	1484.31	3600	3600	20	20	na	na
	10	5	2	10	na	na	1431.35	1318.85	3600	3600	20	20	na	na
	10	5	2.5	10	na	na	1237.06	1112.13	3600	3600	20	20	na	na
36	20	8	1.5	10	na	na	1898.38	1634.4	3600	3600	20	20	na	na
	20	8	2	10	na	na	1675.70	1424.23	3600	3600	20	20	na	na
	20	8	2.5	10	na	na	1528.29	1396.48	3600	3600	20	20	na	na

Note: Opt = optimal solution; CPLEX = CPLEX solver (software); ALNS = adaptive large-neighborhood search; no. ins = number of instances; R1, R2 = results of time–SoC (state-of-charge) function when R = 1 or 2; na = not applicable (CPLEX failed to converge within 1 h).

We also conduct additional tests on some PHEVTSPD variants described below, and the results are shown in Table 3.

MaxLeg: the problem variant that limits the number of customers per truck leg

LRT: the variant considers the launch/retrieve time

Range: the variant that assumes the drone’s flight range depends on the parcel’s weight

Loop: the variant which allows the self-loop

Table 3.

Performance Comparison Between ALNS and CPLEX on Instances with Side Constraints

$\| N \|$	n	$n_{cs}$	$β$	No. ins	MaxLeg		LRT		Range		Loop
$\| N \|$	n	$n_{cs}$	$β$	No. ins	$R = 1$	$R = 2$	$R = 1$	$R = 2$	$R = 1$	$R = 2$	$R = 1$	$R = 2$
8	4	2	1.5	20	0.01	0.00	0.00	0.00	0.00	0.00	na	na
	4	2	2	20	0.00	0.00	0.00	0.00	0.00	0.00	na	na
	4	2	2.5	20	0.00	0.00	0.00	0.00	0.00	0.00	na	na
9	5	2	1.5	20	1.02	1.45	1.75	1.62	1.32	1.04	na	na
	5	2	2	20	1.69	1.24	1.56	1.46	0.79	0.78	na	na
	5	2	2.5	20	1.30	1.41	1.32	1.48	1.74	1.35	na	na
10	6	2	1.5	10	3.05	3.07	2.98	3.45	3.14	3.21	na	na
	6	2	2	10	3.04	3.14	3.62	3.11	3.12	3.23	na	na
	6	2	2.5	10	2.97	3.25	3.75	3.18	3.29	3.37	na	na

Note: ALNS = adaptive large-neighborhood search; CPLEX = CPLEX Optimization Studio (software); MaxLeg = the problem variant that limits the number of customers per truck leg; LRT = launch retrieval time; no. ins = number of instances; na = not applicable (CPLEX failed to converge within 1 h).

Table 2 highlights the performance disparities between the commercial solver and the ALNS algorithm for base cases. Notably, the commercial solver struggles to solve instances with six customers and two charging stations within a 1 h timeframe. In these instances, with an augmented network size of 10 nodes, ALNS demonstrates an average optimality gap of less than 3.5% when allotted a computational time of 5 s. To be more specific, the ALNS method obtained the optimal solution in 48 out of all the 60 instances tested with different $β$ and $R$ values. For the remaining 12 instances where the ALNS solution is suboptimal, a comparison with the optimal CPLEX solution indicates that the suboptimality of ALNS solutions is mostly a result of the charging-station ordering related to the feasibility-resorting function. As instances expand to more than 10 nodes in the augmented network, CPLEX fails to converge within the 1 h limit.

An insightful comparison of the final costs for PHEVTSPD instances reveals an average 7% lower cost compared with the base cases, with this margin seemingly widening as instance sizes increase. This underscores the efficacy of a piecewise linear approximation of the time–SoC function, particularly with a high $R$ value, in yielding higher-quality results. Additionally, the drone’s speed emerges as a critical factor influencing the final cost. Specifically, the final delivery time for PHEVTSPD with $β = 2.5$ is approximately 12% lower than that with $β = 1.5$ .

The results from Table 3 highlight similar optimality gap trends for PHEVTSPD variants as observed in the base problem. However, for the Loop variant, CPLEX fails to terminate within 1 h, even for a modest network with four customers and two charging stations.

Comparison of PHEVTSPD with Different $R$ Values

In this subsection, we compare the final cost of PHEVTSPD with different time–SoC function approximations. The actual time–SoC function and its approximation with different $| R |$ values are shown in Figure 6. Obviously, as the $| R |$ value increases, the accuracy of the approximation also increases.

The final delivery time with different $| R |$ values is shown in Table 4. As expected, the average final cost of PHEVTSPD instances decreases as the $| R |$ value increases. The average final cost with $| R | = 6$ is about 15,728.3 s, while that with linear approximation is about 17,643.2 s, 10.88% higher than when $| R | = 6$ . These results indicate that using piecewise linear function approximation with a low $R$ value might generate a suboptimal solution.

Table 4.

Performance Comparison of PHEVTSPD with Different $| R |$ Values

N	C	#CS	#Ins	Final cost (s)
N	C	#CS	#Ins	R = 1	R = 2	R = 4	R = 6
10	6	2	10	11,315.2	10,890.1	10,848.2	10,840.9
20	10	5	10	16,300.2	15,373.8	15,078.0	14,943.2
30	16	7	10	18,155.7	17,126.0	16,777.4	16,640.2
40	20	10	10	19,973.4	17,573.5	16,979.3	16,906.0
50	26	12	10	22,471.6	19,350.2	19,342.4	19,311.1
Average				17,643.2	16,062.7	15,805.1	15,728.3

In addition, we also test the final delivery time of the downtown Austin network with various $| R |$ values. A total of 35 centroids of ZIP Code Tabulated Areas (ZCTAs) in the Austin metro area are chosen as nodes in the network. The resulting instance is tested with a simple heuristic method to guarantee its feasibility. Out of these 35 nodes, 10 nodes are selected as the charging stations, and the remaining 25 nodes are considered customer nodes and the depot. This network is shown in Figure 8. Two different layouts are tested. In the first central layout, a node that lies in the central area of the network is chosen as the depot. In the side layout, the depot is the node at the network’s bottom.

Figure 8.

Downtown Austin network.

The sensitivity analysis is shown in Figure 9. The central-layout instance generally has less completion time than the side-layout instance. For both network layouts, the final delivery time decreases gradually as $| R |$ increases, and the gap between $| R | = 1$ and $| R | = 6$ is about 15.0% and 12.7% for central and side layouts, respectively.

Figure 9.

Final cost of Austin network with various $| R |$ values.

Conclusions and Future Work

This paper focuses on optimizing the cooperative routing between a PHEV and a drone for efficient deliveries to a designated set of customers. The electric vehicle can recharge at charging-station nodes with flexible energy levels, and the charging time is contingent on the SoC level. Both the charge-energy amount and charging time are treated as decision variables. To address this complex problem, a MILP formulation is proposed within an augmented network to solve PHEVTSPD.

This research utilizes a piecewise linear approximation to incorporate the nonlinear time–SoC function into the model. To tackle instances of practical sizes, a specially designed adaptive neighborhood-search metaheuristic method is introduced. Throughout the search process, the performance of each destroy and repair method is continually evaluated and updated. Additionally, a novel repair method based on CP, explicitly designed for PHEVTSPD, is proposed. The numerical analysis section conducts a comprehensive comparison between the MILP formulation, ALNS, and VNS. Results indicate that ALNS outperforms VNS when the number of nodes in the augmented network exceeds 10, with ALNS maintaining an average optimality gap of less than 3.5%. Moreover, test results with varying $| R |$ values show that a piecewise function with six segments exhibits an average 10% cost reduction compared with a linear function approximation.

This research suggests potential extensions in several directions. First, obtaining exact solutions for PHEVTSPD instances with practical sizes $(| C | > 20)$ remains a challenge. Future research could also explore uncertainties and randomness in the problem, such as variations in customer demands or the drone’s flight time/speed.

Supplemental Material

sj-pdf-1-trr-10.1177_03611981241278351 – Supplemental material for Plug-In Hybrid Electric Vehicle Traveling-Salesman Problem with Drone

Supplemental material, sj-pdf-1-trr-10.1177_03611981241278351 for Plug-In Hybrid Electric Vehicle Traveling-Salesman Problem with Drone by Tengkuo Zhu, Stephen D. Boyles and Avinash Unnikrishnan in Transportation Research Record

Footnotes

Author Contributions

The authors confirm contribution to the paper as follows: study conception and design: Tengkuo Zhu, Stephen Boyles, Avinash Unnikrishnan; data collection: Tengkuo Zhu; analysis and interpretation of results: Tengkuo Zhu, Stephen Boyles, Avinash Unnikrishnan; draft manuscript preparation: Stephen Boyles, Avinash Unnikrishnan. All authors reviewed the results and approved the final version of the manuscript.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: 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).

ORCID iD

Tengkuo Zhu

Supplemental Material

Supplemental material for this article is available online.

References

United States Environmental Protection Agency (EPA). Sources of Greenhouse Gas Emissions. 2021. https://www.epa.gov/ghgemissions/sources-greenhouse-gas-emissions.

FedEx. FedEx Continues Advancing Fleet Electrification Goals. 2022. https://newsroom.fedex.com/.

Sahin

Yilmaz

Ust

Guneri

A. F.

Gulsun

An Approach for Analysing Transportation Costs and a Case Study. European Journal of Operational Research, Vol. 193, No. 1, 2009, pp. 1–11.

Howey

D. A.

Martinez-Botas

Cussons

Lytton

Comparative Measurements of the Energy Consumption of Electric, Hybrid and Internal Combustion Engine Vehicles. Transportation Research Part D: Transport and Environment, Vol. 16, No. 6, 2011, pp. 459–464.

Balthasar

Tait

Riera-Palou

Harrison

A New Comparison Between the Life Cycle Greenhouse Gas Emissions of Battery Electric Vehicles and Internal Combustion Vehicles. Energy Policy, Vol. 44, 2012, pp. 160–173.

DroneZon. Drones for Deliveries from Medicine to Post, Packages and Pizza. 2019. https://www.dronezon.com/drones-for-good/drone-parcel-pizza-delivery-service/.

Teal Group. Teal Group Predicts Worldwide Civil Drone Production Will Soar Over the Next Decade. 2018. https://www.tealgroup.com/.

Loads of Vans. Electric Van Range Guide 2023: Optimise Your EV. 2024. https://www.loadsofvans.com/blog/electric-van-range-guide.

AmazonAir. AmazonPrimeAir. 2024. https://www.aboutamazon.com/news/transportation/amazon-prime-air-drone-delivery-mk30-photos.

10.

Murray

C. C.

Chu

A. G.

The Flying Sidekick Traveling-Salesman Problem: Optimization of Drone-Assisted Parcel Delivery. Transportation Research Part C: Emerging Technologies, Vol. 54, 2015, pp. 86–109.

11.

Agatz

Bouman

Schmidt

Optimization Approaches for the Traveling Salesman Problem with Drone. Transportation Science, Vol. 52, No. 4, 2018, pp. 965–981.

12.

Zhu

Incorporating UAVs into Urban Logistics Combinatorial Optimization Problems. PhD thesis. The University of Texas at Austin, 2022.

13.

Boyles

S. D.

Zhu

Assessment of Parcel Delivery Systems Using Unmanned Aerial Vehicles. University of North Carolina at Charlotte. Center for Advanced Multimodal Mobility Solutions and Education, 2020.

14.

Zhu

Boyles

S. D.

Unnikrishnan

Battery Electric Vehicle Traveling Salesman Problem with Drone. Networks and Spatial Economics, Vol. 24, No. 1, 2024, pp. 49–97.

15.

Montoya

Guéret

Mendoza

J. E.

Villegas

J. G.

The Electric Vehicle Routing Problem with Nonlinear Charging Function. Transportation Research Part B: Methodological, Vol. 103, 2017, pp. 87–110.

16.

Dantzig

G. B.

Ramser

J. H.

The Truck Dispatching Problem. Management Science, Vol. 6, No. 1, 1959, pp. 80–91.

17.

Baldacci

Mingozzi

Roberti

Recent Exact Algorithms for Solving the Vehicle Routing Problem Under Capacity and Time Window Constraints. European Journal of Operational Research, Vol. 218, No. 1, 2012, pp. 1–6.

18.

Laporte

Fifty Years of Vehicle Routing. Transportation Science, Vol. 43, No. 4, 2009, pp. 408–416.

19.

Desaulniers

Desrosiers

Spoorendonk

The Vehicle Routing Problem with Time Windows: State-of-the-Art Exact Solution Methods. In Wiley Encyclopedia of Operations Research and Management Science ( Cochran

J. J.

, ed.), John Wiley & Sons, Inc., 2010.

20.

Osman

I. H.

Metastrategy Simulated Annealing and Tabu Search Algorithms for the Vehicle Routing Problem. Annals of Operations Research, Vol. 41, 1993, pp. 421–451.

21.

Gendreau

Hertz

Laporte

New Insertion and Postoptimization Procedures for the Traveling Salesman Problem. Operations Research, Vol. 40, No. 6, 1992, pp. 1086–1094.

22.

Braekers

Ramaekers

Van Nieuwenhuyse

The Vehicle Routing Problem: State of the Art Classification and Review. Computers & Industrial Engineering, Vol. 99, 2016, pp. 300–313.

23.

Montoya-Torres

J. R.

Franco

J. L.

Isaza

S. N.

Jiménez

H. F.

Herazo-Padilla

A Literature Review on the Vehicle Routing Problem with Multiple Depots. Computers & Industrial Engineering, Vol. 79, 2015, pp. 115–129.

24.

Pillac

Gendreau

Guéret

Medaglia

A. L.

A Review of Dynamic Vehicle Routing Problems. European Journal of Operational Research, Vol. 225, No. 1, 2013, pp. 1–11.

25.

Toth

Vigo

Vehicle Routing: Problems, Methods, and Applications. Society for Industrial and Applied Mathematics, 2014.

26.

Eksioglu

Vural

A. V.

Reisman

The Vehicle Routing Problem: A Taxonomic Review. Computers & Industrial Engineering, Vol. 57, No. 4, 2009, pp. 1472–1483.

27.

Golden

B. L.

Raghavan

Wasil

E. A.

The Vehicle Routing Problem: Latest Advances and New Challenges, Vol. 43. Springer Science & Business Media, New York, NY, 2008.

28.

Toth

Vigo

The Vehicle Routing Problem. SIAM, 2002.

29.

Erdoğan

Miller-Hooks

A Green Vehicle Routing Problem. Transportation Research Part E: Logistics and Transportation Review, Vol. 48, No. 1, 2012, pp. 100–114.

30.

Conrad

R. G.

Figliozzi

M. A.

The Recharging Vehicle Routing Problem. Proc., 2011 Industrial Engineering Research Conference, IISE Norcross, GA, 2011, p. 8.

31.

Schneider

Stenger

Goeke

The Electric Vehicle-Routing Problem with Time Windows and Recharging Stations. Transportation Science, Vol. 48, No. 4, 2014, pp. 500–520.

32.

Koç

Ç.

Karaoglan

The Green Vehicle Routing Problem: A Heuristic Based Exact Solution Approach. Applied Soft Computing, Vol. 39, 2016, pp. 154–164.

33.

Hiermann

Puchinger

Ropke

Hartl

R. F.

The Electric Fleet Size and Mix Vehicle Routing Problem with Time Windows and Recharging Stations. European Journal of Operational Research, Vol. 252, No. 3, 2016, pp. 995–1018.

34.

Masmoudi

M. A.

Hosny

Demir

Genikomsakis

K. N.

Cheikhrouhou

The Dial-a-Ride Problem with Electric Vehicles and Battery Swapping Stations. Transportation Research Part E: Logistics and Transportation Review, Vol. 118, 2018, pp. 392–420.

35.

Keskin

Çatay

Partial Recharge Strategies for the Electric Vehicle Routing Problem with Time Windows. Transportation Research Part C: Emerging Technologies, Vol. 65, 2016, pp. 111–127.

36.

Goeke

Schneider

Routing a Mixed Fleet of Electric and Conventional Vehicles. European Journal of Operational Research, Vol. 245, No. 1, 2015, pp. 81–99.

37.

Felipe

Á.

Ortuño

M. T.

Righini

Tirado

A Heuristic Approach for the Green Vehicle Routing Problem with Multiple Technologies and Partial Recharges. Transportation Research Part E: Logistics and Transportation Review, Vol. 71, 2014, pp. 111–128.

38.

Pelletier

Jabali

Laporte

Veneroni

Battery Degradation and Behaviour for Electric Vehicles: Review and Numerical Analyses of Several Models. Transportation Research Part B: Methodological, Vol. 103, 2017, pp. 158–187.

39.

Erdelić

Carić

A Survey on the Electric Vehicle Routing Problem: Variants and Solution Approaches. Journal of Advanced Transportation, Vol. 2019, 2019, p. 5075671.

40.

Doppstadt

Koberstein

Vigo

The Hybrid Electric Vehicle Traveling Salesman Problem. European Journal of Operational Research, Vol. 253, No. 3, 2016, pp. 825–842.

41.

Doppstadt

Koberstein

Vigo

The Hybrid Electric Vehicle Traveling Salesman Problem with Time Windows. European Journal of Operational Research, Vol. 284, No. 2, 2020, pp. 675–692.

42.

Liao

C.-S.

S.-H.

Shen

Z.-J. M.

The Electric Vehicle Touring Problem. Transportation Research Part B: Methodological, Vol. 86, 2016, pp. 163–180.

43.

Roberti

Wen

The Electric Traveling Salesman Problem with Time Windows. Transportation Research Part E: Logistics and Transportation Review, Vol. 89, 2016, pp. 32–52.

44.

Otto

Agatz

Campbell

Golden

Pesch

Optimization Approaches for Civil Applications of Unmanned Aerial Vehicles (UAVs) or Aerial Drones: A Survey. Networks, Vol. 72, No. 4, 2018, pp. 411–458.

45.

Mbiadou Saleu

R. G.

Deroussi

Feillet

Grangeon

Quilliot

An Iterative Two-Step Heuristic for the Parallel Drone Scheduling Traveling Salesman Problem. Networks, Vol. 72, No. 4, 2018, pp. 459–474.

46.

Q. M.

Deville

Pham

Q. D.

Há

M. H.

On the Min-Cost Traveling Salesman Problem with Drone. Transportation Research Part C: Emerging Technologies, Vol. 86, 2018, pp. 597–621.

47.

Yurek

E. E.

Ozmutlu

H. C.

A Decomposition-Based Iterative Optimization Algorithm for Traveling Salesman Problem with Drone. Transportation Research Part C: Emerging Technologies, Vol. 91, 2018, pp. 249–262.

48.

Bouman

Agatz

Schmidt

Dynamic Programming Approaches for the Traveling Salesman Problem with Drone. Networks, Vol. 72, No. 4, 2018, pp. 528–542.

49.

Poikonen

Golden

Wasil

E. A.

A Branch-and-Bound Approach to the Traveling Salesman Problem with a Drone. Informs Journal on Computing, Vol. 31, No. 2, 2019, pp. 335–346.

50.

De Freitas

J. C.

Penna

P. H. V.

A Variable Neighborhood Search for Flying Sidekick Traveling Salesman Problem. International Transactions in Operational Research, Vol. 27, No. 1, 2020, pp. 267–290.

51.

De Freitas

J. C.

Penna

P. H. V.

A Randomized Variable Neighborhood Descent Heuristic to Solve the Flying Sidekick Traveling Salesman Problem. Electronic Notes in Discrete Mathematics, Vol. 66, 2018, pp. 95–102.

52.

Boysen

Briskorn

Fedtke

Schwerdfeger

Drone Delivery from Trucks: Drone Scheduling for Given Truck Routes. Networks, Vol. 72, No. 4, 2018, pp. 506–527.

53.

Jeong

H. Y.

Song

B. D.

Lee

Truck-Drone Hybrid Delivery Routing: Payload-Energy Dependency and No-Fly Zones. International Journal of Production Economics, Vol. 214, 2019, pp. 220–233.

54.

Dayarian

Savelsbergh

Clarke

J.-P.

Same-Day Delivery with Drone Resupply. Transportation Science, Vol. 54, No. 1, 2020, pp. 229–249.

55.

Kim

Moon

Traveling Salesman Problem with a Drone Station. IEEE Transactions on Systems, Man, and Cybernetics: Systems, Vol. 49, No. 1, 2018, pp. 42–52.

56.

Roberti

Ruthmair

Exact Methods for the Traveling Salesman Problem with Drone. Transportation Science, Vol. 55, No. 2, 2021, pp. 315–335.

57.

Freitas

J. C.

Penna

P. H. V.

Toffolo

T. A.

Exact and Heuristic Approaches to Truck Drone Delivery Problems. EURO Journal on Transportation and Logistics, Vol. 12, 2023, p. 100094.

58.

Luo

Poon

Zhang

Liu

Lim

The Multi-Visit Traveling Salesman Problem with Multi-Drones. Transportation Research Part C: Emerging Technologies, Vol. 128, 2021, p. 103172.

59.

Chen

Miao

The Electric Vehicle Routing Problem with Time Windows and Battery Swapping Stations. Proc., 2016 IEEE International Conference on Industrial Engineering and Engineering Management (IEEM), Bali, Indonesia, IEEE, New York, 2016, pp. 712–716.

60.

Battery University. How-to-Prolong-Lithium-Based-Batteries. 2020. https://batteryuniversity.com/article/bu-808-how-to-prolong-lithium-based-batteries.

61.

Desaulniers

Errico

Irnich

Schneider

Exact Algorithms for Electric Vehicle-Routing Problems with Time Windows. Operations Research, Vol. 64, No. 6, 2016, pp. 1388–1405.

62.

Bruglieri

Pezzella

Pisacane

Suraci

A Variable Neighborhood Search Branching for the Electric Vehicle Routing Problem with Time Windows. Electronic Notes in Discrete Mathematics, Vol. 47, 2015, pp. 221–228.

63.

Froger

Mendoza

J. E.

Jabali

Laporte

Improved Formulations and Algorithmic Components for the Electric Vehicle Routing Problem with Nonlinear Charging Functions. Computers & Operations Research, Vol. 104, 2019, pp. 256–294.

64.

Abid

Tabaa

Chakir

Hachimi

Routing and Charging of Electric Vehicles: Literature Review. Energy Reports, Vol. 8, 2022, pp. 556–578.

65.

Hõimoja

Rufer

Dziechciaruk

Vezzini

An Ultrafast EV Charging Station Demonstrator. Proc., International Symposium on Power Electronics Power Electronics, Electrical Drives, Automation and Motion, Sorrento, Italy, IEEE, New York, 2012, pp. 1390–1395.

66.

Sassi

Cherif

W. R.

Oulamara

Vehicle Routing Problem with Mixed Fleet of Conventional and Heterogenous Electric Vehicles and Time Dependent Charging Costs. ORCHIDS - Operations research for Complex HybrId Decision Sytems, 2014.

67.

Schiffer

Walther

The Electric Location Routing Problem with Time Windows and Partial Recharging. European Journal of Operational Research, Vol. 260, No. 3, 2017, pp. 995–1013.

68.

Zuo

Xiao

You

Kaku

A New Formulation of the Electric Vehicle Routing Problem with Time Windows Considering Concave Nonlinear Charging Function. Journal of Cleaner Production, Vol. 236, 2019, p. 117687.

69.

Campuzano

Lalla-Ruiz

Mes

A Multi-start VNS Algorithm for the TSP-D with Energy Constraints. Proc., International Conference on Computational Logistics, Springer, Cham, 2021, pp. 393–409.

70.

Shaw

Using Constraint Programming and Local Search Methods to Solve Vehicle Routing Problems. Proc., International Conference on Principles and Practice of Constraint Programming, Springer, Berlin, Heidelberg, 1998, pp. 417–431.

71.

Ropke

Pisinger

An Adaptive Large Neighborhood Search Heuristic for the Pickup and Delivery Problem with Time Windows. Transportation Science, Vol. 40, No. 4, 2006, pp. 455–472.

72.

Kirkpatrick

Gelatt

C. D.

Jr. Vecchi

M. P.

Optimization by Simulated Annealing. Science, Vol. 220, No. 4598, 1983, pp. 671–680.

73.

Sacramento

Pisinger

Ropke

An Adaptive Large Neighborhood Search Meta-Heuristic for the Vehicle Routing Problem with Drones. Transportation Research Part C: Emerging Technologies, Vol. 102, 2019, pp. 289–315.

74.

Montemanni

Dell’Amico

Constraint Programming Models for the Parallel Drone Scheduling Vehicle Routing Problem. EURO Journal on Computational Optimization, Vol. 11, 2023, p. 100078.

75.

Montemanni

Dell’Amico

Solving the Parallel Drone Scheduling Traveling Salesman Problem via Constraint Programming. Algorithms, Vol. 16, No. 1, 2023, p. 40.

76.

Ham

A. M.

Integrated Scheduling of m-Truck, m-Drone, and m-Depot Constrained by Time Window, Drop-Pickup, and m-Visit Using Constraint Programming. Transportation Research Part C: Emerging Technologies, Vol. 91, 2018, pp. 1–14.

77.

ALKE. Electric Delivery Van. 2021. https://www.alke.com/electric-delivery-van.

78.

DJI. DJI. 2021. https://www.dji.com/.

79.

TeoTeslaFan. Model S Supercharging Times Compared S60, S70D, S85, P85D, S85D. 2015. https://www.tesla.com/support/.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.18 MB

Plug-In Hybrid Electric Vehicle Traveling-Salesman Problem with Drone

Abstract

Keywords

Literature Review

Problem Description and Formulations

MILP Model for PHEVTSPD

Modeling Partial-Recharge Policies

ALNS Solution Method

Overview of ALNS

Initial Solution

Removal Methods

Repair Methods

Greedy Repair

Nearby Repair

CP Repair

Feasibility-Resorting Function and Cost-Evaluation Function

Computational Experiments

Experimental Setting

Estimation of the Concave Time–SoC Function

Performance Comparison Between ALNS and VNS

Comparison Between MILP Formulation and ALNS Method

Comparison of PHEVTSPD with Different R Values

Conclusions and Future Work

Supplemental Material

sj-pdf-1-trr-10.1177_03611981241278351 – Supplemental material for Plug-In Hybrid Electric Vehicle Traveling-Salesman Problem with Drone

Footnotes

Author Contributions

Declaration of Conflicting Interests

Funding

ORCID iD

Supplemental Material

References

Supplementary Material

Comparison of PHEVTSPD with Different $R$ Values