Battery Electric Vehicle Traveling Salesman Problem with Drone

Abstract

The idea of deploying electric vehicles and unmanned aerial vehicles (UAVs), also known as drones, to deliver packages in logistics operations has attracted increasing attention in the past few years. In this paper, we propose an innovative problem where a battery electric vehicle (BEV) paired with drone is utilized to deliver first-aid items in a rural area. This problem is termed battery electric vehicle traveling salesman problem with drone (BEVTSPD). In BEVTSPD, the BEV and the drone perform delivery tasks coordinately while the BEV can serve as a drone hub. The BEV can also refresh its battery energy to full capacity in battery-swap stations available in the network. An arc-based mixed-integer programming model defined in a multigraph is presented for BEVTSPD. An exact branch-and-price (BP) algorithm and a Variable Neighborhood Search (VNS) 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 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 BEV’s driving range.

Appendix

1.1 Handling Additional Side Constraints with MILP Formulation

In this section, we aim to relax some of the operational assumptions of BEVTSPD 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).

1.1.1 Fixed Charging Time or Service Time

Assume that each node $i \in C \cup S$, is associated with a fixed time cost $f_{i}$. If i is a customer node, $f_{i}$ represents the fixed service time. If i is a charging station node, $f_{i}$ 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 i to the arc cost of the node’s incoming or outgoing arcs. For example, define the modified travel times $c^{mT}_{ij} = c^{T}_{ij} + f_{i}, \forall \ i \in C \cup S$. We can handle these side constraints by replacing the travel times $c^{T}_{ij}$ with the modified travel times $c^{mT}_{ij}$.

1.1.2 Loops

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 $(i-j-i)$ is allowed which is also called a loop.

Define a new binary decision variable $z_{ij}$, which equals to 1 if the drone performs loop $( i-j-i)$ and 0 otherwise. To account for the impact of the loop, the objective function is changed to:

$$\begin{aligned} t^{*} = min \ t_{0'} + \sum _{j}(c^{D}_{0j}+c^{D}_{j0^{'}})z_{0j} + \sum _{(i,j):i \ne j}(c^{D}_{ij}+c^{D}_{ji})z_{ij} \end{aligned}$$

(29)

Constraints (7) are replaced by:

$$\begin{aligned} y^{T}_{i}+ y^{D}_{i}+ y^{C}_{i}+\sum _{j}z_{ji} =1, \ \forall \ i \in C_{0} \end{aligned}$$

(30)

$$\begin{aligned} z_{ij} \le y^{C}_{i}, \ \forall \ i,j \in C_{0}: i \ne j \end{aligned}$$

(31)

Constraints (16) are replaced by:

$$\begin{aligned} b^{d2}_{i} \le b^{a2}_{i} - \sum _{j}\sum _{k} y_{ijk}d_{ijk} - \sum _{j:j \ne i}(e^{D}_{ij}+e^{D}_{ji})z_{ij}, \ \forall \ i \in C_{0} \end{aligned}$$

(32)

The new objective (29) considers the extra waiting time incurred by the loops. Constraints (30) indicate that a node i can be a truck node, combined node, or a drone node that is visited by a drone sortie or a loop. Constraints (31) state that a node can perform a loop only when it is a combined node. Constraints (32) are the updated shared energy constraints that consider the existence of loops.

1.1.3 Incompatible Customers

An incompatible customer is a customer that can only be served by the truck and not by the drone. Denote $C^{D} \in C$ 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:

$$\begin{aligned} y^{D}_{i} =0, \ \forall \ i \in C \backslash C^{D} \end{aligned}$$

(33)

$$\begin{aligned} z_{ji} = 0, \ \forall \ i \in C \backslash C^{D},j \in C_{0}: i \ne j \end{aligned}$$

(34)

1.1.4 Launch and Retrieve Time

Assume that it takes the truck $S_{L}$ 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 $l_{i}, i \in C_{0}$, which equals one if node i is a drone node and not directly visited by the drone from the depot. Then the objective (1) is replaced by:

$$\begin{aligned} t^{*} = min \ t_{0'} + \sum _{j}(c^{D}_{0j}+c^{D}_{j0^{'}})z_{0j} + \sum _{(i,j):i \ne j}(c^{D}_{ij}+c^{D}_{ji}+S_{L})z_{ij} + \sum _{i \in C_{0}}S_{L}l_{i} \end{aligned}$$

(35)

Besides, we need to add constraints:

$$\begin{aligned} l_{i} \ge y^{D}_{i} - x^{D}_{0i} \end{aligned}$$

(36)

The new objective (35) considers the extra launch time of the loops and at each launch node. The constraints (36) state that $l_{i}$ equals one only when node i is a drone node and not directly visited by the drone from the depot.

Similarly, assume that it takes the truck $S_{R}$ 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:

$$\begin{aligned} t^{*} = min \ t_{0'} + \sum _{j}(c^{D}_{0j}+c^{D}_{j0^{'}}+S_{R})z_{0j} + \sum _{(i,j):i \ne j}(c^{D}_{ij}+c^{D}_{ji}+S_{R})z_{ij} + \sum _{i \in C_{0}}S_{R}y^{D}_{i} \end{aligned}$$

(37)

1.2 Handling Additional Side Constraints in Dynamic Programming

1.2.1 Launch and Retrieve Times

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 $f(ng, k, i^{C}, i^{T}, \tau , b^{C}, b^{T}, w) - u_{j} + S_{L}$. When adding a drone leg, the new state value is $f(ng, k, i^{C}, i^{T}, \tau , b^{C}, b^{T}, False) - u_{j} + c(i^{T},j,p) + S_{R}$.

1.2.2 Maximum Number of Customers Per Truck Leg

In dynamic programming, we can add a truck arc when $i^{T} != i^{C}$ 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 $n^{T}$ indicating the current number of customers visited in the current truck leg. In this way, the new state is

$$\begin{aligned} f(ng, k, i^{C}, i^{T}, \tau , b^{C}, b^{T}, w, n^{T}). \end{aligned}$$

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 $\bar{n}$ is small (e.g. $\bar{n} \le 3$), the influence of the new state variable is minimal. Besides, as $\bar{n}$ 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.

1.2.3 Weight-Dependent Drone Flying Range

The weight-dependent drone flying range assumption can be enforced in DP by changing the original constraint $e^{D}_{(i^{C}, j, 0)} + e^{D}_{(j, i^{T}, 0)} \le Q^{d}$ when adding drone leg to the new constraint $f^{D}_{(i^{C}, j, 0)} + f^{D}_{(j, i^{T}, 0)} \le Q^{d}$ where $f^{D}_{(i, j, p)}$ is the energy consumption function for drones when traversing arc (i, j, p).

1.2.4 Loops

When the loop is allowed in the operation process, in DP, an extra action of adding a loop is needed. For current state $f(ng, k, i^{C}, i^{T}, \tau , b^{C}, b^{T}, w)$, the add loop action is only available when $i^{C} = i^{T}, \tau =0, w=False, b^{C}=b^{T}$, which indicates that both the truck and the drone are at node $i^{C}$ at the current state, and $b^{C} \ge e^{D}_{i^{T}j}+e^{D}_{ji^{T}}$ which indicates that the remaining energy is sufficient for the drone to perform sortie $\langle i^{T}, j, i^{T} \rangle$. After the add loop action, the resulting new state is $f((ng \cup \{i^{T}, i^{C}\}) \cap N_{j}, k+1, i^{T}, i^{T}, 0, b^{T}, b^{T}, False)$ and the new state value is $f(ng, k, i^{C}, i^{T}, 0, b^{C}, b^{T}, w) + t^{D}_{(i^{C}, j, 0)} + t^{D}_{(j, i^{T}, 0)} - u_{j}$.

1.3 Proof of Statement 3

Proof

Define a reverse ng-route that starts with the destination depot while back-propagating to the origin depot. Thus we could also define a function $f(ng_{b}, n-k, i^{C}_{b}, i^{T}, \tau _{b}, b^{C}_{b}, b^{T}, w) = f^{b}$ as the representation of this reverse partial ng-route. In this way, a full ng-route could be achieved by combining $f(ng, k, i^{C}, i^{T}, \tau , b^{C}, b^{T}, w)$ with $f(ng_{b}, n-k, i^{C}_{b}, i^{T}, \tau _{b}, b^{C}_{b}, b^{T}, w)$.

If all the above conditions are satisfied, all the potential propagation that could be achieved by propagating $f(ng_{2}, k, i^{C}, i^{T}, \tau , b^{C}_{2}, b^{T}_{2}, w)$ can also be achieved by propagating $f(ng_{1}, k, i^{C}, i^{T}, \tau , b^{C}_{1}, b^{T}_{1}, w)$.

By combining $f^{1}$ with a backward partial route, we can derive a full ng-route. Assuming a route performs drone leg $i^{C}-j-i^{C}_{b}$ where j is the drone node, denote $t_{cb} = t^{D}_{i^{C}j} + t^{D}_{ji^{C}_{b}} + f_{j}$ as the needed time for the drone to finish this drone leg. Then the reduced cost of the full route is:

$$\begin{aligned} \hat{d_{1}} = f_{1}+f^{b}+u_{i^{T}}-u_{j}+max\{\tau _{1}+\tau _{b}, t_{cb}\} \end{aligned}$$

Similarly, the reduced cost of the full route by combining $f^{2}$ with $f^{b}$ is

$$\begin{aligned} \hat{d_{2}} = f_{2}+f^{b}+u_{i^{T}}-u_{j}+max\{\tau _{2}+\tau _{b}, t_{cb}\} \end{aligned}$$

If all the conditions from (28a) to (28f) are satisfied, we need to prove that $\hat{d_{1}} \le \hat{d_{2}}$. This can be illustrated in three cases:

• If $t_{cb} \ge max\{\tau _{1}+\tau _{b}, \tau _{2}+\tau _{b}\}$, then $\hat{d_{1}} = f_{1}+f^{b}+u_{i^{T}}-u_{j}+ t_{cb}$ and $\hat{d_{2}} = f_{2}+f^{b}+u_{i^{T}}-u_{j}+ t_{cb}$. Thus $\hat{d_{1}} \le \hat{d_{2}}$ because of condition (28a).

• If $t_{cb} \le min\{\tau _{1}+\tau _{b}, \tau _{2}+\tau _{b}\}$, then $\hat{d_{1}} = f_{1}+f^{b}+u_{i^{T}}-u_{j}+ \tau _{1}+\tau _{b}$ and $\hat{d_{2}} = f_{2}+f^{b}+u_{i^{T}}-u_{j}+\tau _{2}+\tau _{b}$. Thus $\hat{d_{1}} \le \hat{d_{2}}$ because of condition (28f).

• If $\tau _{1}+\tau _{b} \le t_{cb} \le \tau _{2}+\tau _{b}$, then $\hat{d_{1}} = f_{1}+f^{b}+u_{i^{T}}-u_{j}+ t_{cb}$ and $\hat{d_{2}} = f_{2}+f^{b}+u_{i^{T}}-u_{j}+ \tau _{2}+\tau _{b}$. Thus $\hat{d_{1}} \le \hat{d_{2}}$ because of condition (28a).

Thus, $\hat{d_{1}} \le \hat{d_{2}}$ in all the scenarios and the dominance rule stands.