The cost of carrier consistency: Last-mile delivery by vehicle and drone for subscription-based orders

Abstract

We investigate a deterministic last-mile delivery problem in which a firm delivering scheduled orders determines a subset of customers to be prioritized for drone delivery, relegating the rest to vehicle delivery. The vehicle and drone deliver in tandem over a planning horizon. The problem is amenable to a multi-day Traveling Salesman Problem with Drone, with consistency in the assignment of customers to each carrier across all days. We develop a mixed-integer program (MIP) that yields provably optimal solutions for instances having up to 40 customers over a 6-day horizon. We also devise a heuristic underpinned by the MIP to effectively and accurately prioritize customers for drone delivery. A computational study involving instances with up to 150 customers ordering over 6 days shows the matheuristic often matching the results of the proposed MIP, and also reveals a logistical burden of less than 0.5% on average over a solution lacking carrier consistency..

Keywords:

• Last-mile delivery
• drone
• routing
• traveling salesman
• optimization

1. Introduction

By the time Wing expanded their drone delivery service to Dallas, Texas in April 2022, the company had overseen more than 200,000 successful flights (French, 2022). Describing his company’s strategy for the expansion, Wing CTO Adam Woodworth explained: “We’re going to invite customers in groups to make sure everyone has a good first experience with drone delivery” (French, 2022). DroneUp CTO John Verson, whose company has been delivering parcels via drone in partnership with Walmart for two years prior to expansion plans announced in May 2022, described a similar strategy: “We are starting small with the idea that we can deliver to a certain set of customers” (Barefoot, 2022; Repko, 2022). Rather than deploying drone delivery en masse, these firms have purposefully limited drone delivery to a select group of customers to start, knowing that the consistency and reliability of drone delivery for these initial customers is the key determinant of whether the service might expand. In the words of Vernon, CTO of DroneUp: “As we continue to improve a certain set of operations, we will be expanding” (Barefoot, 2022; Michael, 2022). Firms initially limit drone deliveries to ensure their logistical feasibility and success. That consistency also serves customers, many of whom who may be experiencing drone delivery for the first time (Figure 1).

Figure 1. The drones operated by both Wing and DroneUp are designed to lower parcels via cables to customers who are physically present for delivery. (a) A DroneUp drone carrying a parcel (Michael, 2022). (b) A Wing drone customer retrieving his parcel (French, 2022).

This work addresses two questions central to drone delivery: given deterministic customer delivery schedules, which subset of customers should be prioritized for drone delivery initially? And what is the logistical cost associated with consistent drone service for those customers? This article introduces the Traveling Salesman Problem with Subscription-Based Orders (TSP-SBO), which investigates this new notion of carrier consistency. Motivated by the potential to capture a larger market, a supplier selects a set of customers who may be offered drone delivery, seeking to serve those customers consistently via drone and thereby to demonstrate the viability of the related operations to the remaining customers. Customers in the TSP-SBO place subscription orders, in which they express recurring demand met on a pre-determined delivery schedule. For example, a customer may order coffee every morning, bananas every other day, or prescriptions once a month. Depending on the recurrence of these subscription orders, the delivery schedule in each period of the TSP-SBO features a subset of customers requiring delivery. Unlike other multi-period problems, the TSP-SBO does not permit shifts of deliveries across periods. The key requirement is that customers must be served by the same carrier (i.e., either the vehicle or the drone) across periods. Figure 2 illustrates a toy-sized example for how delivery schedules impact carrier assignment consistency for six customers across three periods.

Figure 2. The delivery schedule determines the periods in which each customer requires service, which in turn influences the carriers’ assignments and routes across periods. (a) A toy-sized delivery schedule. (b) A legend for the illustrated routes. (c) Period 1 (d) Period 2 (e) Period 3.

In each period illustrated in Figure 2(c)–(e), a routing solution is determined to serve customers with scheduled deliveries. A customer be must be served by the same carrier, across periods, whenever a delivery is scheduled. For example, Customer 3 requires delivery in Periods 1 and 2, and is served by the vehicle in both periods. Consistent service requires the presence of nearby customers being served by the vehicle in the same period. For example, Customer 4 is served by the drone in Period 2, with a flight sequence of (3,4,0), where 0 corresponds to the depot. In Period 3, however, Customer 3 does not have a scheduled delivery, while Customer 4 does. As such, the drone flight sequence in Period 3 for Customer 4 is (5,4,0).

As contributions, firstly this work presents a novel mixed-integer programming (MIP) formulation for the TSP-SBO that underpins a matheuristic approach for generating high-quality solutions. The matheuristic solves a multi-dimensional knapsack sub-problem to initialize a solution for the TSP-SBO, then uses the proposed MIP formulation to dynamically improve that solution through search procedures. Secondly, over a variety of benchmark instances, the efficacy of the matheuristic compares favorably against that of the MIP formulation, yielding solutions that consistently deviate from their optimal counterparts by less than 0.5% on average with substantial computational savings. The computational study also explores the logistical cost of serving a subset of customers via drone consistently before rolling out the service to a larger customer base. The results suggest the logistical cost of consistent service via drone is limited to less than 0.5% as compared with solutions that do not require consistency.

The remainder of the paper is organized as follows. Section 2 presents a review of related literature. Section 3 introduces the formal problem statement alongside the proposed MIP formulation. Section 4 discusses the matheuristic. Section 5 presents a computational study comparing exact and matheuristic results across a variety of instances. Section 6 concludes this work with a summary of our findings.

2. Literature review

Before introducing two streams of literature that relate to the notions of carrier synchronization and multi-period demand in the TSP-SBO, this section begins by outlining common operational assumptions for vehicle-and-drone routing problems.

2.1. Operational assumptions in vehicle-and-drone routing problems

There is a growing literature focused on the single-period problem of routing a vehicle and its companion drone through a network of customer deliveries at minimal cost or duration, which is known as the Traveling Salesman Problem with Drone (TSP-D). Like the TSP-SBO, the TSP-D allows for collaborative deliveries between the carriers, and therefore requires that the vehicle and drone travel in sync to accommodate en-route rendezvous between them. Due to the novelty of vehicle-and-drone deliveries, a diverse set of operational assumptions has been considered in the TSP-D literature.

Table 1 shows that vehicle-and-drone routing problems are governed by a variety of assumptions about the number of periods considered in the problem, the drone’s flight time and speed, the locations at which it may rendezvous with the vehicle, the number of drones considered (some references allow for 1 or more drones), and the number of customers the vehicle may visit while the drone is in flight (the intervening vehicle visits column in Table 1). Although this diversity in operational assumptions has led to a divergence in problem settings for the TSP-D, that is at least in part due to the ever-evolving nature of the real-world problem. According to a recent review by Boysen et al. (2021), “there may be some elaborate problem versions not yet investigated, but seeing that drones are not yet operating in mass markets and many operational details are still unclear …scientific research should rather concentrate on the identification of the most promising drone-delivery concept, which may vary for last-mile delivery tasks.” In that vein, the subsections that follow review vehicle-and-drone routing literature with common assumptions aligned with those in this study.

Table 1. A summary of recent works.

Article	Planning horizon (Periods)	Max drone flight time	Drone count	Intervening vehicle visits	Drone-eligible customers	Drone speed
Murray and Chu (2015)	Single	Fixed	1	Yes	100%	Fixed
Agatz et al. (2018)	Single	Unlimited	1	Yes	100%	Fixed
Bouman et al. (2018)	Single	Unlimited	1	No	N/A	Fixed
Boysen et al. (2018)	Single	Unlimited	$\ge$1	No	100%	Fixed
Carlsson and Song (2018)	Single	Unlimited	1	Yes	100%	Fixed
El-Adle et al. (2021)	Single	Fixed	1	Yes	100%	Fixed
Raj and Murray (2020)	Single	Weight-dependent	$\ge$1	Yes	85%	Variable
Roberti and Ruthmair (2021)	Single	Unlimited	$\ge$1	Yes	100%	Fixed
This work	Multiple	Fixed	1	No	100%	Fixed

2.2. Single-period variants

Each period of a TSP-SBO instance, when separated from the remainder of the problem, is in fact an instance of the TSP-D. There is a growing literature for the TSP-D, spanning exact optimization (Murray & Chu, 2015; Carlsson & Song, 2018; Agatz et al., 2018; El-Adle et al., 2021; Poikonen et al., 2019; Roberti & Ruthmair, 2021) and heuristic approaches (Bouman et al., 2018; Schermer et al., 2019; Dell’Amico et al., 2021) that have solved networks having up to 39 nodes optimally and up to 200 nodes heuristically. For recent reviews, the interested reader is referred to Macrina et al. (2020) and Boysen et al. (2021).

Distinctly from the TSP-D, the TSP-SBO requires carrier consistency over a multi-day horizon. Firstly, the delivery schedule in each period of the TSP-SBO may feature distinct customers. Thus, solving one instance of the TSP-SBO is akin to solving a series of instances of the TSP-D, each of which may have a distinct set of customers who require a delivery in that period. Since the TSP-SBO requires that the same carrier consistently serve the same customer across all periods, optimal customer assignments thus minimize the routing cost of serving a customer across all periods in the horizon. Even if carrier assignments in the TSP-SBO were pre-determined as part of a solution approach, there would remain to be solved in each in each time period a distinct TSP-D instance.

2.3. Multi-period demand & carrier consistency in routing problems

Although few references exist in the routing literature for subscription-based ordering (Belavina et al., 2017), the notions of multi-period demand and carrier consistency have been investigated, with problem definitions that set them apart from the TSP-SBO. Recent examples of the former include Archetti et al. (2015), in which a fleet of homogeneous capacitated vehicles is routed for delivery in each period. Demand from one period may be postponed to a later date, subject to inventory and penalty costs, and the objective is to minimize those costs alongside the overall transportation cost. Lagan et al. (2021) investigated multi-period routing in the context of parcel delivery with dynamic demand, such that customer orders are revealed over time, and the decision-maker selects strategies and priorities for delivery schedules. Distinctly from these works, the TSP-SBO requires orders to be delivered according to a deterministic delivery schedule, without the ability to delay delivery. Just as important, the TSP-SBO features two heterogeneous carriers with synchronized movements.

Consistency in routing problems can take several forms, including vehicle (i.e., using the same vehicle/driver) or temporal (i.e., delivering in a consistent time window) consistency, with changes allowed to the vehicle’s route in each period depending upon demand (Mourgaya & Vanderbeck, 2007). The TSP-SBO features a type of carrier consistency, in that customers chosen for vehicle or drone delivery are always served by the same carrier. Nevertheless, since the drone is launched and retrieved at customer nodes (distinct from the customer who receives the drone delivery), consistency in the TSP-SBO depends upon the synchronous movement of the carriers as well as customers’ co-expression of demand.

The combination of carrier consistency over a multi-period horizon, the requirement to deliver orders in the period in which they are due, as well as the synchronized movements of the vehicle and drone, set the TSP-SBO apart from the extant literature as a uniquely unified, computationally-challenging problem.

3. Problem statement & MIP formulation

This section describes the TSP-SBO, including details of drone delivery operations that justify the problem settings, then proposes a MIP formulation.

3.1. Problem statement

Consider a network of customer locations represented by n nodes, a network of roadways for the vehicle as well as a network of flight paths for the drone whose edges connect those nodes, and a $\mathcal{T}$-period horizon over which orders are delivered via a pre-determined schedule. Periods refer to generic windows of time in which customers express demand. For the application discussed in this article, one period represents one day. In this deterministic setting, the subset of customers requiring delivery in each period is known; likewise, for any customer, the subset of periods in which that customer requires delivery is also known, as illustrated in Figure 2(a). In each period, the drone and vehicle traverse an optimized tour that begins and ends at a central depot while collaboratively serving all customers with orders in that period. Every customer is served consistently by the same carrier across all periods, and deliveries may not be postponed for later periods. A decision-maker determines in advance a maximal share of the customer base that may be served via drone in the initial roll-out of the service. In summary, the TSP-SBO thus seeks an assignment of customers to the vehicle and drone that achieves consistent service by the same carrier to each customer while respecting the decision-maker’s limit on the number of customers to be served via drone. The objective is to determine assignments that minimize the overall return time of both carriers to the depot while feasibly delivering all orders across the planning horizon.

3.2. Drone delivery operations

As outlined in Section 2, the TSP-SBO adopts common operational assumptions from the literature on vehicle-and-drone delivery problems. To clarify the practical relevance of these assumptions, this section describes various elements of drone delivery.

3.2.1. Drone flight range

DroneUp drones are designed to carry parcels weighing up to 10 pounds, and to deliver them in as quickly as 30 min (Guggina, 2022). The drones are battery powered, exhausting their charge as they fly. After a flight, the drone’s battery pack may be swapped, or there may be multiple drones such that the drone used for the most recent flight charges while another delivers new orders (Wing, 2022). In either case, the drone replenishes its battery capacity after each flight. Accordingly, in the TSP-SBO, the drone is permitted to serve one customer per flight and the duration of any single flight must not exceed F, the drone’s flight range. When retrieved by the vehicle, we assume the drone is loaded with a new parcel, and that the drone may again embark on a new flight with a duration of less than F. When the drone is not flying, it is transported aboard the vehicle.

For customer k to receive a drone delivery, the vehicle may launch the drone from a customer i served by the vehicle, and then the vehicle may retrieve the drone from customer j who is also served by the vehicle, where all nodes are distinct. Several studies have examined the impact of allowing a drone to launch and land at the same node, known in the literature as a drone cycle (Agatz et al., 2018; Schermer et al., 2019; Roberti & Ruthmair, 2021). Most of the literature cited in Section 2, however, assume distinct launching and landing nodes for the drone. Hereafter, (i, j) are referred to as rendezvous pairs for k.

3.2.2. Simultaneous carrier movement and waiting

While a commercial drone is in flight, regulations in the US generally require the pilot to track “the position, altitude, attitude, and movement” of the drone unaided (i.e., using only the pilot’s sight) (FAA, 2022). Accordingly, in this problem the vehicle is permitted to move during a drone flight, but is not permitted to make intervening deliveries to other customers. That is, given the rendezvous pair (i, j) used to serve customer k via drone, the vehicle must travel directly from i to j.

Finally, although launching the drone requires both carriers to be present simultaneously at a node, the drone’s landing and the vehicle’s arrival to a node may be asynchronous. A key distinction related to waiting time is the number of parcels with each carrier. After dropping a parcel during a delivery, the drone may alter its flight speed so that it arrives at a rendezvous node simultaneously with the vehicle (Raj & Murray, 2020). If the drone exhausts its flight capacity, it may await the vehicle on the ground. In either case, since the drone is no longer carrying a parcel, its waiting time does not affect the overall completion time of the deliveries. Should the vehicle await the drone at a landing site, that waiting time counts towards the duration of the tour.

3.2.3. Consistent availability of rendezvous pairs

A central challenge in the TSP-SBO is that delivery via drone necessitates the co-expression of customer orders. As illustrated in Figure 3a, customers D1 and V1 both have at least two rendezvous pairs available in the first period of demand. In a second period shown in Figure 3b, however, customer V1 lacks two nearby rendezvous pairs, which eliminates the possibility of drone service in that period. And by extension, V1 must be served by the vehicle in all periods with an active subscription. Note also that D1 does not use the same rendezvous pairs in Period 2 as those in Period 1. Nevertheless, as long as at least one pair of rendezvous pairs may be used for drone service whenever a given customer requires delivery, then the customer may be served by the drone in all periods. This phenomenon is hereafter referenced as drone eligibility: the availability of rendezvous pairs with co-expressed subscriptions in the same periods as the node of interest. Whenever a customer lacks at least one rendezvous node pair in any period with an active subscription, that customer is deemed to be drone-ineligible.

Figure 3. An illustration of the availability of rendezvous pairs for customers in distinct periods. (a) In Period 1, rendezvous pairs are available for both D1 and V1. (b) In Period 2, only D1 has rendezvous pairs, which is why V1 must be served by the vehicle.

3.3. MIP formulation

The parameters and variables for the proposed formulation are as follows:

3.3.1. Input parameters

• $T=\{1,\dots ,\mathcal{T}\}:$ Set of service periods.

• $V=\{0,1,\dots ,n,n+1\}:$ Set of service nodes in the network.

• $V^{+}=V\backslash \{0,n+1\}:$ Set of customer nodes, which excludes the central depot (represented by dummy nodes 0 and n + 1).

• $A=\{(i,j):i,j\in V,i\ne j\}:$ Set of arcs connecting distinct nodes in the network.

• α: Ratio of drone to vehicle speed.

• τ_ij: Vehicle travel time for $(i,j)\in A.$

• f_ij: Drone flight time for $(i,j)\in A.$

• $d_j^t:$ Binary indicator such that $d_j^t=1\iff$ customer j requires parcel delivery in period t, $\forall j\in V^{+},t\in T.$

• F: Maximum flight duration for the drone to make a single delivery then return to the vehicle.

• $C^t\subseteq V^{+}$ $=\{j\in V^{+}$ $|d_j^t=1\}:$ Set of customers requiring parcel delivery in period $t\in T.$

• $R^t\subseteq V=C^t\cup \{0,n+1\}:$ Set of nodes requiring carrier visits in period $t\in T.$

• ${\widehat{f}}_{\mathit{\text{ikj}}}=\mathrm{\text{max}}(0,f_{\mathit{\text{ik}}}+f_{\mathit{\text{kj}}}-{\tau }_{\mathit{\text{ij}}}):$ The vehicle’s waiting time for a drone flight originating from i, delivering to k, and landing at j, $\forall$ $(i,j)\in A,k\in D,$ with distinct i, j, k.

• $E_k^t=\{(i,j)\in R^t$ $|f_{\mathit{\text{ik}}}+f_{\mathit{\text{kj}}}\le F,$ $i,j\in R^t$$\}:$ All pairs of rendezvous pairs permitting a drone delivery to customer k in period t, with distinct i, j, k.

• $D\subseteq V^{+}=\{j\in V^{+}$ $|$ $d_j^t>0$ $\wedge$ $\left|E_j^t\right|>0$ $\forall t\in T\}:$ Set of nodes for which drone delivery is feasible by virtue of having rendezvous pairs available in each period in which the customer requires delivery.

• $\Gamma =$ The maximal share of customers for whom drone delivery may be offered.

• $D_{\mathrm{\text{max}}}^t:$ The objective value of an initial feasible solution for period $t\in T.$

• $M_{\mathit{\text{ij}}}^t:$ A sufficiently large scalar associated with period $t\in T,(i,j)\in A.$

The proposed formulation makes use of Miller-Tucker-Zemlin sub-tour elimination constraints and a set of “big M” scalars. The procedure for obtaining both the scalar M values and a starting feasible solution, $D_{\mathrm{\text{max}}}^t,$ may be found in Appendix A.

3.3.2. Decision variables

• $x_{\mathit{\text{ij}}}^t\in \{0,1\}:$ $x_{\mathit{\text{ij}}}^t=1$ $\iff$ the vehicle travels from i to j in period $t\in T,$ $i,j\in R^t,i\ne j.$

• $y_{\mathit{\text{ikj}}}^t\in \{0,1\}:$ $y_{\mathit{\text{ikj}}}^t=1$ $\iff$ the drone launches from i, delivers to k, and lands at j in period t, $\forall$ $t\in T,$ $k\in C^t,\hspace{0.17em}(i,j)\in E_k^t,$ with distinct i, k, j.

• ${\delta }_j\in [0,1]:$ ${\delta }_j=1$ $\iff$ the drone delivers to customer j across all periods in which j expresses demand, $\forall j\in V^{+}.$ Although the δ-variables are declared continuous and bounded between 0 and 1, they assume binary values due to their relationship with the binary x-variables in the model.

• $b_j^t\in [0,D_{\mathrm{\text{max}}}^t]:$ The departure time of the vehicle at node j in period $t\in T,j\in R^t.$

The formulation for the Traveling Salesman Problem with Subscription-Based Orders, Model TSPSBO, is presented as follows: (1a) $$\mathbf{\text{TSPSBO}}:\text{Minimize}\mathrm{ }\hspace{1em}\sum _{t\in T}\sum _{i\in R^t}\sum _{j\in R^t|j\ne i}{\tau }_{\mathit{\text{ij}}}{x}_{\mathit{\text{ij}}}^t+\sum _{t\in T}\sum _{k\in (C^t\cap D)}\sum _{(i,j)\in E_k^t}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ijk}}}^t$$(1a) (1b) $$\mathrm{s}.\mathrm{t}.\hspace{1em}\sum _{j\in C^t}{x}_{0j}^t=1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,$$(1b) (1c) $$\sum _{j\in C^t}{x}_{\mathit{\text{jn}}+1}^t=1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,$$(1c) (1d) $$\sum _{j\in R^t}{x}_{\mathit{\text{ij}}}^t-\sum _{j\in R^t}{x}_{\mathit{\text{ji}}}^t=0,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in C^t,$$(1d) (1e) $$\sum _{j\in R^t}{x}_{\mathit{\text{kj}}}^t=1-{\delta }_k,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,k\in C^t,$$(1e) (1f) $$\sum _{(i,j)\in E_k^t}{y}_{\mathit{\text{ikj}}}^t+\sum _{i\in R^t}{x}_{\mathit{\text{ik}}}^t=1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,k\in C^t,$$(1f) (1g) $$\begin{array}{c}{y}_{\mathit{\text{ikj}}}^t\le x_{\mathit{\text{ij}}}^t,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,k\in (C^t\cap D),\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}(i,j)\in E_k^t,\end{array}$$(1g) (1h) $$\sum _{k\in C^t|(i,j)\in E_k^t}{y}_{\mathit{\text{ikj}}}^t\le 1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i,j\in R^t,$$(1h) (1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) (1j) $$b_{n+1}^t\le \sum _{i\in R^t}\sum _{l\in R^t}{\tau }_{\mathit{\text{il}}}{x}_{\mathit{\text{il}}}^t+\sum _{k\in (C^t\cap D)}\sum _{(i,l)\in E_k^t}{\widehat{f}}_{\mathit{\text{ikl}}}{y}_{\mathit{\text{ikl}}}^t,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,$$(1j) (1k) $$b_j^t\le D_{\mathrm{\text{max}}}^t,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}\forall t\in T,j\in R^t,$$(1k) (1l) $$\sum _{j\in V^{+}}{\delta }_j\le \lfloor \Gamma \left|V^{+}\right|\rfloor ,$$(1l) (1m) $$b,\delta ,\ge 0\hspace{1em}x,y\mathrm{ }\text{binary}.\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}$$(1m)

For each time period, the objective (1a)(1a) $$\mathbf{\text{TSPSBO}}:\text{Minimize}\mathrm{ }\hspace{1em}\sum _{t\in T}\sum _{i\in R^t}\sum _{j\in R^t|j\ne i}{\tau }_{\mathit{\text{ij}}}{x}_{\mathit{\text{ij}}}^t+\sum _{t\in T}\sum _{k\in (C^t\cap D)}\sum _{(i,j)\in E_k^t}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ijk}}}^t$$(1a) minimizes the return time to the depot for all carriers: the summation of x-variables minimizes the travel time of the vehicle, while the summation of y-variables accounts for the duration of the total vehicle idle time, which accrues when the vehicle awaits the return of the drone from a delivery it is completing. Thus (1a)(1a) $$\mathbf{\text{TSPSBO}}:\text{Minimize}\mathrm{ }\hspace{1em}\sum _{t\in T}\sum _{i\in R^t}\sum _{j\in R^t|j\ne i}{\tau }_{\mathit{\text{ij}}}{x}_{\mathit{\text{ij}}}^t+\sum _{t\in T}\sum _{k\in (C^t\cap D)}\sum _{(i,j)\in E_k^t}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ijk}}}^t$$(1a) is the makespan of a tour in which both carriers return to the depot. Owing to Constraints (1b)(1b) $$\mathrm{s}.\mathrm{t}.\hspace{1em}\sum _{j\in C^t}{x}_{0j}^t=1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,$$(1b) and (1c)(1c) $$\sum _{j\in C^t}{x}_{\mathit{\text{jn}}+1}^t=1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,$$(1c) , the vehicle’s path begins and ends at the depot (represented as nodes 0 and n + 1, respectively), while Constraint (1d)(1d) $$\sum _{j\in R^t}{x}_{\mathit{\text{ij}}}^t-\sum _{j\in R^t}{x}_{\mathit{\text{ji}}}^t=0,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in C^t,$$(1d) in combination with the subtour elimination constraints discussed below ensures a Hamiltonian path for the vehicle through the intervening customer nodes for each time period. Across time periods, Constraint (1e)(1e) $$\sum _{j\in R^t}{x}_{\mathit{\text{kj}}}^t=1-{\delta }_k,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,k\in C^t,$$(1e) ensures the consistent assignment of a single carrier to each customer, while Constraint (1f)(1f) $$\sum _{(i,j)\in E_k^t}{y}_{\mathit{\text{ikj}}}^t+\sum _{i\in R^t}{x}_{\mathit{\text{ik}}}^t=1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,k\in C^t,$$(1f) ensures that each customer is delivered to either by the vehicle or by the drone. Constraint (1 g) disallows intervening vehicle deliveries during a drone flight, and Constraint (1h)(1h) $$\sum _{k\in C^t|(i,j)\in E_k^t}{y}_{\mathit{\text{ikj}}}^t\le 1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i,j\in R^t,$$(1h) allows at most one drone flight to launch from each rendezvous node.

Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) determines the departure time of the vehicle from nodes, with the secondary purpose of a subtour elimination constraint of the Miller-Tucker-Zemlin- (MTZ-) type. Since MTZ-based formulations are known to yield relatively weak linear relaxations, Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) has been enhanced by scaling M, as suggested in Remark 1 and detailed in the Appendix. Constraint (1j)(1j) $$b_{n+1}^t\le \sum _{i\in R^t}\sum _{l\in R^t}{\tau }_{\mathit{\text{il}}}{x}_{\mathit{\text{il}}}^t+\sum _{k\in (C^t\cap D)}\sum _{(i,l)\in E_k^t}{\widehat{f}}_{\mathit{\text{ikl}}}{y}_{\mathit{\text{ikl}}}^t,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,$$(1j) provides an upper bound for the value of the arrival time at the depot in each time period, since the latest possible departure time from any node must be less than the makespan of the entire tour in the same period. Constraint (1k)(1k) $$b_j^t\le D_{\mathrm{\text{max}}}^t,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}\forall t\in T,j\in R^t,$$(1k) generalizes this upper bound to all variables associated with departure time from any node. Constraint (1 l) limits the number of customers (with the total customer base represented by the cardinality of V⁺) served via drone according to the firm’s early roll-out preference (a share of customers represented by Γ). Finally, Constraint (1 m) enforces non-negativity and binary restrictions on the decision variables.

Remark 1.

As shown in Appendix A, the M scalar in Constraints (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) may be set to:$$M_{\mathit{\text{ij}}}^t=D_{\mathrm{\text{max}}}^t-{\tau }_{\mathit{\text{in}}+1}+{\tau }_{\mathit{\text{ij}}}-{\tau }_{0j}\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in R^t\backslash \{n+1\},j\in R^t\backslash \left\{0\right\}|i\ne j.$$

4. Matheuristic for TSP-SBO

This section introduces a matheuristic, underpinned by Model TSPSBO, to tractably solve large instances. The matheuristic is comprised of three phases: initially, a customer-centric multi-dimensional knapsack problem yields a starting feasible assignment of customers to carriers; in two subsequent phases of neighbourhood searches dictated by customer subscription patterns and geo-spatial analyses, the matheuristic improves that initial solution. The diversification schemes and termination criteria for the heuristic are also highlighted thereafter. Figure 4 provides an overview of the proposed matheuristic procedure.

Figure 4. An overview of the matheuristic for the TSP-SBO.

In what follows, customer assignments refer to values of the binary δ-variables, indicating which of the vehicle (${\delta }_j=0$) or drone (${\delta }_j=1$) will serve customer $j\in V^{+}$ across all periods. In contrast, routing assignments refer to values of the binary x- and y-variables that determine the path of the vehicle and drone in each period, respectively. For customers that are drone-eligible, a cluster refers to the nearby nodes within flight range of the node. The incumbent refers to the best-known feasible solution for Model TSPSBO.

4.1. Initialization

The matheuristic begins by generating feasible customer assignments for Model TSPSBO by solving a multi-dimensional knapsack problem (MKP). As discussed in Section 3.2.3, the feasibility of drone service depends on the consistent availability of rendezvous pairs near the customer to be served via drone. As a result, there exists an upper limit on the total number of customers that may be served via drone: that upper bound is $\lfloor \frac{(m-1)}{2}\rfloor ,$ where m is the number of nodes in the cluster of interest. In fact, the TSP-SBO may be characterized by two types of upper bounds on drone assignments: the aforementioned global upper bound for the entire network; and localized cluster-specific bounds based on the availability of rendezvous pairs therein. Combined, these bounds present a special structure in the TSP-SBO. The following notation is introduced to exploit that structure:

4.1.1. Input parameters

• s_k: The minimal savings associated with serving a drone-eligible customer $k\in D$ via drone (as compared with vehicle service) from the nearest rendezvous pairs available in each period, defined as $s_k=\sum _{t\in T}{\mathrm{\text{min}}}_{(i,j)\in E_k^t}({\tau }_{\mathit{\text{ik}}}+{\tau }_{\mathit{\text{kj}}})\alpha .$ Note α is the ratio of the drone over the vehicle’s speed.

• $c_{\mathrm{\text{max}}}:$ The maximum number of nodes that may be served by drone in the network, where $c_{\mathrm{\text{max}}}=\lfloor \frac{((n+2)-1)}{2}\rfloor =\lfloor \frac{(n+1)}{2}\rfloor .$ Note that n + 2 is the total number of nodes in the network since n is the number of customer nodes, and the depot is represented by 2 nodes.

• ${\omega }_{\mathit{\text{tk}}}=\underset{(i,j)\in E_k^t}{\cup }\{i,j\}:$ The union of the cluster-specific nodes that may serve as rendezvous node pairs for drone-eligible customer $k\in D$ in period $t\in T.$

• $p_{\mathit{\text{tk}}}^1:$ The maximum number of cluster-specific nodes that are ineligible for drone service and may be used as rendezvous points for node $k\in D$ in period $t\in T,$ where $p_{\mathit{\text{tk}}}^1=\mathrm{\text{max}}(0,\lfloor \frac{(|{\omega }_{\mathit{\text{tk}}}\backslash D|-1)}{2}\rfloor ).$

• $p_{\mathit{\text{tk}}}^2:$ The maximum number of cluster-specific nodes that are drone-eligible and may be used as rendezvous points for node $k\in D$ in period $t\in T,$ where $p_{\mathit{\text{tk}}}^2=\mathrm{\text{min}}(1,\mathrm{\text{max}}(0,\lfloor \frac{(|{\omega }_{\mathit{\text{tk}}}\cap D|-1)}{2}\rfloor ).$

The parameters $p_{\mathit{\text{tk}}}^1$ and $p_{\mathit{\text{tk}}}^2$ partition the cluster for customer k based on the feasibility for drone service whenever there is more than one node eligible for such service. Figure 5 illustrates why both parameters are needed; the figure uses polar coordinates, with radii corresponding to half of the drone’s flight range. All nodes within one radial arc (and 60 degrees) from node k form the cluster for k. Both Nodes 6 and 16 are drone-eligible in Figure 5, and Node 6 may be served by the following rendezvous pairs: $(0,16);(16,3);(16,4);(16,11);(16,13);$ as well as their reciprocals (e.g., (3, 16)). Note that all the rendezvous pairs associated with Node 6 are dependent on the availability of Node 16 as part of a rendezvous pair. Thus if Node 16 were assigned to the drone, then Node 6 would no longer be eligible for drone service.

Figure 5. A radial plot of carrier assignments to customers in a densely-populated neighbourhood.

To enforce this mutually exclusive relationship, $p_{\mathit{\text{tk}}}^2$ allows no more than one node per cluster to be served via drone using the $\mathrm{\text{max}}$ operator. Considering Node 6, for example, ${\omega }_{2,6}=\{0,3,4,11,13,16,21\},$ where Nodes 0 and 21 are dummy nodes representing the depot. Of those nodes, only $({\omega }_{2,6}\cap D)=\{3,4,11,13,16\}$ are drone-eligible, meaning that in general up to $\lfloor \frac{(5-1)}{2}\rfloor =2$ nodes may be assigned for drone service. By construction, however, $p_{\mathit{\text{tk}}}^2=\mathrm{\text{min}}(1,\mathrm{\text{max}}(0,2))=1.$ The inner $\mathrm{\text{max}}$ operator is required in the case $|({\omega }_{2,6}\cap D)|=0,$ such that $p_{\mathit{\text{tk}}}^2$ is bounded below by 0.

Similarly, the $\mathrm{\text{max}}$ operator is required to enforce a lower bound of 0 for $p_{\mathit{\text{tk}}}^1,$ which considers drone-ineligible nodes in the cluster. Using Node 6 in Figure 5 as an example, only one drone-ineligible node (the depot) exists in that cluster. Since drone-ineligible nodes must always be served by the vehicle, the only upper bound on the number of such nodes that may be used as rendezvous points is $\lfloor \frac{(1-1)}{2}\rfloor =0.$ Thus for Node 6 in period 2, $p_{2,6}^1=0,\hspace{0.17em}{p}_{2,6}^2=1,$ meaning no more than one node may be assigned to drone service in the cluster of node 6. To enforce such a constraint, Model MKP is presented as follows, with binary z-variables such that $z_k=1$ $\iff$ $k\in D$ is served by the drone: (2a) $$\mathbf{\text{MKP}}:\text{Maximize}\mathrm{ }\hspace{1em}\sum _{k\in D}{s}_k{z}_k$$(2a) (2b) $$\mathrm{s}.\mathrm{t}.\hspace{1em}\sum _{k\in D}{z}_k\le c_{\mathrm{\text{max}}},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}$$(2b) (2c) $$z_k+\sum _{{}_{j\in ({\omega }_{\mathit{\text{tk}}}\cap D)}}{z}_j\le p_{\mathit{\text{tk}}}^1+p_{\mathit{\text{tk}}}^2,\forall t\in T,k\in C^t,$$(2c) $$z\mathrm{ }\text{binary}.\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}$$

Model MKP focuses on assigning customers for service via the drone, without regard for corresponding routing assignments. The result is a limited number of decision variables, on the order of the cardinality of D, the set of drone-eligible nodes. Across all drone-eligible nodes, objective (2a)(2a) $$\mathbf{\text{MKP}}:\text{Maximize}\mathrm{ }\hspace{1em}\sum _{k\in D}{s}_k{z}_k$$(2a) in Model MKP maximizes the savings associated with avoiding vehicular travel time in favor of using the drone. Constraint (2b)(2b) $$\mathrm{s}.\mathrm{t}.\hspace{1em}\sum _{k\in D}{z}_k\le c_{\mathrm{\text{max}}},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}$$(2b) sets $c_{\mathrm{\text{max}}}$ as the global maximum for the number of nodes that may be served via drone. At clusters associated with each $k\in D,$ Constraint (2c)(2c) $$z_k+\sum _{{}_{j\in ({\omega }_{\mathit{\text{tk}}}\cap D)}}{z}_j\le p_{\mathit{\text{tk}}}^1+p_{\mathit{\text{tk}}}^2,\forall t\in T,k\in C^t,$$(2c) enforces localized knapsacks. The left-hand side of Constraint (2c)(2c) $$z_k+\sum _{{}_{j\in ({\omega }_{\mathit{\text{tk}}}\cap D)}}{z}_j\le p_{\mathit{\text{tk}}}^1+p_{\mathit{\text{tk}}}^2,\forall t\in T,k\in C^t,$$(2c) captures all drone-eligible nodes in the cluster of k, including k. The right-hand side of Constraint (2c)(2c) $$z_k+\sum _{{}_{j\in ({\omega }_{\mathit{\text{tk}}}\cap D)}}{z}_j\le p_{\mathit{\text{tk}}}^1+p_{\mathit{\text{tk}}}^2,\forall t\in T,k\in C^t,$$(2c) specifies the maximum number of drone assignments according to the cluster partitioning scheme described previously. Given the drone assignments from Model MKP, the remaining customer assignments in Model TSPSBO may be fixed for the vehicle.

The feasible customer assignments in Model MKP come at the expense of simplifying (namely, by over-constraining) a key feature of the TSP-SBO problem statement. Specifically, in Model MKP the parameter s_k assumes that the nearest rendezvous pair is always used to serve customer node k via drone. Under the original problem statement, any rendezvous pair may be used in any period. This simplification is restored in subsequent phases of the matheuristic.

Given carrier assignments from Model MKP, determining the carrier paths in each period is a separable problem: Model TSPSBO decomposes into $\mathcal{T}$ independent instances of the TSP-D for each period (with carrier assignments fixed across the multi-period horizon). Hereafter, for a given set of customer assignments, let the operational completion be a conforming set of feasible routing decisions. Appendix B shows a single-period variant of Model TSPSBO, which treats customer assignments (i.e., the δ-variables) as input parameters and yields a corresponding operational completion. To accelerate the search for improved solutions, the following stages of the matheuristic are underpinned by this single-period variant of the formulation as well.

4.2. Subscription-based improvement

Phase 1 of the matheuristic focuses on improving period-independent customer assignments (as opposed to period-dependent routing decisions). As such, Phase 1 solves Model TSPSBO iteratively, with the variables associated with certain customer clusters fixed, while certain other customer clusters (and their associated variables) are targeted for re-optimization.

Customer clusters are targeted based on their subscription schedules, since those with more frequent deliveries have a correspondingly relatively large number of routing variables associated with them. Those variables factor directly into Objective (1a)(1a) $$\mathbf{\text{TSPSBO}}:\text{Minimize}\mathrm{ }\hspace{1em}\sum _{t\in T}\sum _{i\in R^t}\sum _{j\in R^t|j\ne i}{\tau }_{\mathit{\text{ij}}}{x}_{\mathit{\text{ij}}}^t+\sum _{t\in T}\sum _{k\in (C^t\cap D)}\sum _{(i,j)\in E_k^t}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ijk}}}^t$$(1a) . Indirectly, since these clusters require frequent delivery, their assignments are likely to influence the assignments of other customers with less frequent subscriptions. For instance, erroneously assigning the drone to deliver to a node with frequent subscriptions may preclude the neighbouring nodes from being served via drone. The following notation is relevant to Phase 1:

4.2.1. Phase 1 notation:

• $P^t=\{j\in V^{+}|\sum _{l\in T}{d}_j^l=t\}:$ Sets of customers organized by the frequency of their subscription (e.g., P³ corresponds to all customers requiring delivery in 3 distinct periods).

• $Z_{\text{incumbent}}^{*}:$ The best-known solution identified by the matheuristic.

• NII: The number of non-improving iterations deployed in Phase 1 before termination.

• $\lambda \in \{1,\dots \mathcal{T}\}:$ An index for the number of time periods in the planning horizon.

At the initiation of Phase 1, $Z_{\text{incumbent}}^{*}$ is the operational completion of the solution from Model MKP, NII $=0,$ and $\lambda =\mathcal{T}.$ Phase 1 proceeds as follows:

4.2.2. Phase 1 procedure:

• Step 1. Unfix the values of all variables in Model TSPSBO.

• Step 2. $\forall i\notin P^{\lambda },$ fix δ_i to the values in $Z_{\text{incumbent}}^{*}.$

• This step ensures only a subset of customer nodes corresponding to $P^{\lambda }$ is targeted for carrier assignment re-optimization.

• Step 3. $\forall i\notin P^{\lambda },$ fix x, y variables associated with node i to the values in $Z_{\text{incumbent}}^{*}.$

• This step ensures that the routing assignments to be re-optimized will focus on nodes contained within $P^{\lambda }.$

• Step 4. Solve Model TSPSBO.

• Step 5. Update $Z_{\text{incumbent}}^{*}$ if an improved solution has been found. Else, $\mathit{\text{NII}}=\mathit{\text{NII}}+1.$

• Step 6. Let $\lambda =\lambda -1.$

• Step 7. If λ = 0 or the threshold for NII has been exceeded, terminate.

Although re-optimization of both period independent (δ) and period dependent (x, y) variables is possible through this procedure, the focus is on the former. The next phase of the matheuristic provides the primary framework for the improvement of the latter.

4.3. Geo-spatial Improvement

Re-optimizing a solution on the basis of subscriptions in Phase 1 disregards an essential component of the problem: the geo-spatial proximity necessary to find rendezvous pairs for drone service. Indeed, F, the flight range of the drone, circumscribes a customer cluster wherein period independent and period dependent variables are intertwined. Distinctly from Phase 1, Phase 2 targets customer clusters and nodes adjacent to that cluster as part of its iterative improvement procedure. Moreover, Phase 2 prioritizes for re-optimization those clusters with The following notation introduces Phase 2:

4.3.1. Phase 2 notation:

• ${\omega }_{\mathit{\text{tk}}}=\underset{(i,j)\in E_k^t}{\cup }\{i,j\}:$ The union of the rendezvous pairs for drone-eligible customer $k\in D$ in period $t\in T.$

• ${\Omega }_k=\{t\in T|k\in C^t\}:$ The set of all periods in which customer k requires delivery.

• ${\omega }_k=\underset{t\in {\Omega }_k}{\cup }{\omega }_{\mathit{\text{tk}}}:$ The set of rendezvous pairs for customer k across all time periods.

• $Z_{\text{incumbent}}^{*}:$ The best-known solution identified by the matheuristic.

• NII: The number of non-improving iterations deployed in Phase 1 before termination.

Phase 2 of the matheuristic also relies on several data structures to track elements of the solution that have already been explored.

• exploredNodes $\subset V^{+}:$ A set comprised of customer nodes $k\in V^{+}.$

• exploredClusters $\subset V^{+}:$ A set comprised of customer nodes $k\in V^{+}.$

• expandedCluster $\subset V^{+}:$ A set comprised of customer nodes $k\in V^{+}.$

At the initiation of Phase 2, $Z_{\text{incumbent}}^{*}$ is the best-known solution from the previous phase, NII $=0,$ and exploredNodes $=$ exploredClusters $=$ expandedClusters $=\varnothing .$ Phase 2 proceeds as follows:

4.3.2. Phase 2 procedure:

• Step 1. $\forall i,j\in V^{+},$ if ${\omega }_i\subset {\omega }_j$ then,

• $\Rightarrow \mathit{\text{exploredNodes}}=\mathit{\text{exploredNodes}}\cup \{l\in V^{+}|l\in {\omega }_i\}$

• $\Rightarrow \mathit{\text{exploredClusters}}=\mathit{\text{exploredClusters}}\cup \left\{i\right\}.$

  For distinct drone-eligible nodes j and k, if ${\omega }_j\subset {\omega }_k,$ re-optimizing the cluster of k ensures the cluster for node j is re-optimized as well. This step eliminates dominated clusters from consideration.

• Step 2. $\widehat{k}=\text{arg} {\mathrm{\text{max}}}_{\mathrm{k}\in ({\mathrm{V}}^{+}\backslash \mathit{\text{exploredClusters}})}(\text{card}({\omega }_{\mathrm{k}}\backslash \mathit{\text{exploredNodes}}))$

  Select a drone-eligible node $\widehat{k}$ such that the cluster associated with $\widehat{k}$ includes the greatest number of nodes that have yet to be re-optimized.

• Step 3. $\mathit{\text{expandedCluster}}=\{l\in V^{+}|l\in {\omega }_k\hspace{0.17em}\vee \hspace{0.17em}x>0\hspace{0.17em}\vee \hspace{0.17em}y>0\}$

  Expand the cluster by augmenting ω_k with nodes from the carriers’ path in the incumbent solution.

• Step 4. Unfix all variable values in Model TSPSBO.

• Step 5. $\forall i,j,l\in \mathit{\text{expandedCluster}},$ fix ${\delta }_i,x_{\mathit{\text{ij}}},y_{\mathit{\text{ijl}}}$ to the values from $Z_{\text{incumbent}}^{*}.$

  To focus the re-optimization effort on only those nodes relevant to the targeted cluster, Phase 2 fixes the value of all variables associated with nodes excluded from expandedCluster.

• Step 6. Solve Model TSPSBO.

• Step 7. Update $Z_{\text{incumbent}}^{*}$ if an improved solution has been found. Else, $\mathit{\text{NII}}=\mathit{\text{NII}}+1.$

• Step 8. If $V^{+}\backslash \mathit{\text{exploredClusters}}=\varnothing$ or NII threshold exceeded, then terminate.

The special structure of the TSP-SBO permits optimal assignments in the local cluster of a drone-eligible customer to be in fact globally optimal for the instance. By virtue of the limited flight range of the drone, the customer assignments within each cluster are independent from the remainder of the network. In fact, exhaustively improving customer and routing assignments in all clusters would inexorably produce a globally optimal solution. This insight motivates the diversification scheme in Phase 2, which distinctly from Phase 1, avoids repeated re-optimization of the same dominated clusters.

4.4. Diversification schemes & termination criteria

The interplay between the two phases serves a valuable role in diversifying the search for a solution. Improved solutions from each phase may cascade into the next, such that alternating between the two phases, rather than using either one in isolation, enhances the search to yield higher-quality solutions. Moreover, depending on the demand or geospatial profile of instances being solved, a decision-maker may elect to execute a phase of the heuristic multiple times. The choice of how frequently to deploy each phase, and when to terminate the matheuristic’s search, is left to the decision-maker; Section 5.3 discusses our computational experience.

5. Computational study

This section introduces instances on which the proposed exact optimization and matheuristic approaches are tested, followed by computational and analytical insights.

5.1. Description of problem instances

Instances are characterized by two main features: (i) a network configuration detailing the geographical location of customers; and (ii) a delivery profile for the schedule of subscription-based deliveries for customers over the entire multi-period planning horizon.

5.1.1. Network configurations

Five network configurations were generated, each featuring the depot at the origin. Each configuration includes 150 customer nodes distributed uniformly in a 360${}^{°}$ pattern (spanning $2.5F$) around the origin. For each instance discussed below where n < 150, the first n customer nodes from each configuration were chosen. The network is connected, with τ_ij (the vehicle’s travel time between nodes) equal to the Euclidean distance between i and j, and f_ij (the drone’s travel time) equal to $\alpha \times {\tau }_{\mathit{\text{ij}}},$ where $\alpha =0.5$ such that the drone travels twice as fast as the vehicle, and the drone’s flight range F = 30 units. These settings match those used in other recent TSP-D works (Agatz et al., 2018; Schermer et al., 2019; El-Adle et al., 2021). In particular, the instances have been designed such that the flight range of the drone matches that of existing drone deliveries conducted in the United States for 30 min per delivery (Guggina, 2022).

5.1.2. Delivery profiles

The increasing volume and frequency of global parcel delivery discussed in Section 1 has motivated multiple customer categories in this study: $C_A,$ C_B, and C_C whose subscriptions recur with likelihood ρ_A, ρ_B, ρ_C, respectively. Each customer requires delivery in at least one time period in the planning horizon and is assigned to exactly one of the disjoint customer categories. Let 0 $\le$ ρ_C $\le$ ρ_B $\le$ ρ_A $\le$ 1, such that customers in category C_A have the highest probability of requiring delivery in a particular period.

There is an infinite spectrum of possibilities for both the proportion of customer categories and the subscription frequency associated with each, which in turn shape TSP-SBO solutions. It is impractical to investigate all of these profiles. Certain profiles, however, may be more amenable to analytical interpretation. For example, setting ${\rho }_C={\rho }_B={\rho }_A=1$ ensures that all customers have active subscriptions in each period, and thus the TSP-SBO reduces to the simpler (but still challenging) TSP-D. Alternatively, setting one or more of ρ_C, ρ_B, and ρ_A to 0 would be eliminate that category from the delivery network entirely. The customer categories and delivery profiles used in this study balance the recurrence of customer subscriptions in each period to simulate realistic conditions. Specifically, Appendix C summarizes parcel delivery data from several countries that informed the selection of the delivery profiles outlined as follows.

• Customer Category A: 10% share, with ${\rho }_A=\frac{2}{3},$ and $C_A=\{1,\dots ,\lceil 0.1\left|V^{+}\right|\rceil \}.$

• Customer Category B: 30% share, with ${\rho }_B=0.35,$ and $C_B=\{\lceil 0.1\left|V^{+}\right|\rceil +1,\dots ,\lceil 0.4\left|V^{+}\right|\rceil \}.$

• Customer Category C: 60% share, with ${\rho }_C=0.1,$ and $C_C=\{V^{+}\backslash (C_A\cup C_B)\}.$

5.1.3. Instances

We examine datasets¹ S, M, and L, which involve small-, medium-, and large-sized instances, respectively. For each dataset, five distinct network configurations are examined and, for each of these, a number of distinct delivery profiles (based on the demand parameters specified above) are generated. With planning horizons of $\mathcal{T}=$ 3 or $\mathcal{T}=$ 6 time-periods (typically days), the computational study includes:

• Dataset S: 150 instances; ($n=$ 20, 30, or 40) × ($\mathcal{T}=$ 3 or 6) × 5 network configurations × 5 delivery profiles;

• Dataset M: 40 instances; ($n=$ 50 or 75) × ($\mathcal{T}=$ 3 or 6) × 5 network configurations × 2 delivery profiles;

• Dataset L: 20 instances; ($n=$ 100 or 150) × ($\mathcal{T}=$ 3 or 6) × 5 network configurations × 1 delivery profile.

Based on the instance generation scheme described above, Table 2 summarizes several key statistics which may be helpful to characterize the instances explored.Table 2 shows that, on average in a given period, 36-50% of customers require delivery. Moreover, the flight range of the drone permits, on average, 47-98% of customer nodes to serve as part of a rendezvous pair. Note that the former statistic is period-dependent, while the latter is period-independent. Both the average subscriptions and number of drone-eligible customers increase with the size (and customer density) of the network. Finally, note that under $\mathcal{T}=3,$ there were generally more drone-eligible customers as compared with those eligible under $\mathcal{T}=6.$ Although adding time periods (and nodes) might increase the overall number of subscriptions, adding customers with infrequent subscriptions may sparsify the subscription schedule and diminish the viability of drone service.

Table 2. Instance profile statistics.

Instance size (nodes)	$\mathcal{T}$	Average active subscriptions per period (nodes)	Average drone-eligible customers per instance (nodes)
20	3	10.0	12.4
	6	7.3	9.4
30	3	14.9	23.6
	6	11.3	18.7
40	3	19.9	32.9
	6	14.4	25.2
50	3	24.5	42.6
	6	17.9	40.4
75	3	37.2	71.1
	6	27.9	66.4
100	3	49.5	96.9
	6	37.1	91.4
150	3	71.0	147.6
	6	54.0	144.6

5.1.4. Computational settings and parameter tuning

Models TSPSBO, TSPSBO-SP (see Appendix B), and MKP as well as the matheuristic were implemented using AMPL and solved using GUROBI version 9.0.2, which was in our computational experience the most effective solver for the formulations and problem instances presented. The runs were performed on a computer with an Intel i7-7700K processor and 32 GB of RAM. Using that platform, the following solver settings and algorithmic parameters yielded the best computational efficiency in our experience. For an initial solution using Model MKP, the solver was limited to 300 CPU seconds, with the vast majority of instances yielding an optimal solution in a small fraction of that time. Each phase of the matheuristic relies upon an incumbent solution from the previous phase, and the solver begins with those incumbent solutions as a warm start. The initial feasible solution used in Phase 1 is generated using the assignments yielded by Model MKP and routes yielded by the greedy algorithm described in Appendix A. In Phase 1, we set NII = 3, with a solver limit of 120 CPU seconds in each iteration. In Phase 2, we set NII $=3$ across all instances, with a solver time limit that varied by instance size: 30 CPU seconds for instances with $n\le 40;$ 40 CPU seconds for $n\in [50,100];$ and 120 CPU seconds for n = 150. The aforementioned time limits were selected based on the performance of the solver with Model TSPSBO during the exact solution approach. More generally, the solver settings were adjusted for the matheuristic, with each re-optimization terminating with a relative gap of 5% or larger, and a focus on finding feasible solutions rather than improving the lower bound of the solution at hand. Finally, the termination criteria (setting NII $=3$) in each phase were determined empirically on the basis of the matheuristic’s performance on the S and M datasets. We aimed to balance computational efficiency with the matheuristic’s efficacy.

5.2. Base results

The following statistics are introduced to assess the quality of the proposed matheuristic. In what follows, let $Z_{\text{MIP}}^{*}$ be the best solution yielded by the exact optimization approach of Model TSPSBO (within a time limit of 10,000 s), let $Z_{\text{MH}}^{*}$ be the best solution yielded by the matheuristic, and let $Z_{\text{MKP}}^{*}$ be the best operational completion of the customer assignments generated from Model MKP.

5.2.1. Summary Statistics

• Δ: The deviation between the matheuristic solution and that yielded by the solver using the proposed (time-limited) exact optimization approach, with $\Delta =\frac{(Z_{\text{MH}}^{*}-Z_{\text{MIP}}^{*})}{Z_{\text{MIP}}^{*}}.$

• ${\Delta }_{\mu }:$ The average of Δ over a set of instances.

• ${\Delta }_{\mathrm{\text{max}}}:$ The maximum of Δ over a set of instances.

• ${\Delta }^{\text{MKP}}:$ The deviation between the matheuristic solution and that yielded by a (time-limited) completion for Model MKP, with ${\Delta }^{\text{MKP}}=\frac{(Z_{\text{MH}}^{*}-Z_{\text{MKP}}^{*})}{Z_{\text{MKP}}^{*}}.$

• ${\text{CPU}}_{\text{MH}}:$ The matheuristic’s average computational effort across a set of instances.

• ${\text{CPU}}_{\text{MIP}}:$ The average computational effort (capped at 10,000 seconds) required by the exact optimization of Model TSPSBO across a set of instances.

Table 3 summarizes benchmark testing across Dataset S, in which all instances were solved optimally using Model TSPSBO. The Opt. and ${\Delta }_{\mathrm{\text{max}}}$ columns show the matheuristic found optimal solutions for more than 73% of instances, with a deviation of at most 2.8% from the optimal solutions. Across Dataset S, the average deviation between matheuristic solutions and their optimal counterparts was 0.5% or less.

Table 3. Computational results for dataset S, networks having 40 nodes or fewer.

n	$\mathcal{T}$	${\Delta }_{\mu }(\%)$	${\Delta }_{\mathrm{\text{max}}}(\%)$	Opt.	${\text{CPU}}_{\text{MH}}$ (s)	${\text{CPU}}_{\text{MIP}}$ (s)	${\Delta }_{\text{MKP}}(\%)$
20	3	0.1	1.2	21/25	3.9	1.7	−19.4
	6	0.2	1.4	19/25	4.6	1.5	−13.0
30	3	0.0	0.2	24/25	34.4	54.8	−29.0
	6	0.1	1.1	20/25	36.4	30.3	−19.6
40	3	0.5	2.8	10/25	139.7	920.4	−33.0
	6	0.2	0.9	16/25	158.2	1,759.2	−24.2

As the size of the network and number of periods increase, Table 3 also shows that the matheuristic’s computational performance scales better than that of the exact optimization approach. Although the exact optimization approach has a slight computational advantage for $n=$ 20 and 30 under $\mathcal{T}=$ 3, this is due to the fixed computational expense associated with the matheuristic’s progression through its phases and the number of non-improving iterations the heuristic considered for convergence. For instances of practical relevance (i.e., $n\ge 40$), however, the matheuristic has a clear computational advantage. Finally, the ${\Delta }^{\text{MKP}}$ column highlights the relative savings achieved by the matheuristic by comparing its final solution (after all phases) with its starting solution (the time-limited routing completion for customer assignments from Model MKP). The magnitude of this metric is dependent on both the quality of the starting solution, and the ability of the matheuristic to improve upon it. For comparative purposes, it is included as part of the assessments in all datasets.

In the absence of optimal solutions for Datasets M and L (exceptions are highlighted in table footnotes), Tables 4 and 5 introduce two new columns using both upper and lower bounding approaches to analyze the quality of the matheuristic solution. ${\text{LB}}_{\mu }$ compares the solutions yielded by the matheuristic against the lower bound available at the termination of the time-limited exact optimization approach; values close to 0 indicate a solution that is provably near-optimal. Δ_L compares the deviation between the upper bounds produced by the exact optimization approach over the same computational effort consumed by the matheuristic; negative values indicate the matheuristic outperformed the exact optimization approach when allowed a comparable computational time.

Table 4. Computational results for dataset M, networks having 50–75 nodes.

n	$\mathcal{T}$	${\Delta }_{\mu }(\%)$	${\Delta }_{\mathrm{\text{max}}}(\%)$	${\Delta }_t(\%)$	${\text{LB}}_{\mu }(\%)$	${\text{CPU}}_{\text{MH}}$ (s)	${\text{CPU}}_{\text{MIP}}$ (s)	${\Delta }_{\text{MKP}}(\%)$
50	3	0.4^3	1.7	−0.3	−0.6	333.6	3,922.8	−38.7
	6	0.6^4	2.0	−3.9	−1.9	504.6	6,860.1	−29.0
75	3	0.8	6.4	−6.0	−5.7	692.8	10,000.0	−42.4
	6	−1.2	0.7	−0.5	−6.5	1,074.9	10,000.0	−35.3

³ Opt. = 2/8.

⁴ Opt. = 0/5.

Table 5. Computational results for dataset L, networks having 100–150 nodes.

n	$\mathcal{T}$	${\Delta }_{\mu }(\%)$	${\Delta }_{\mathrm{\text{max}}}(\%)$	${\Delta }_t(\%)$	${\text{LB}}_{\mu }(\%)$	${\text{CPU}}_{\text{MH}}$ (s)	${\text{CPU}}_{\text{MIP}}$ (s)	${\Delta }_{\text{MKP}}(\%)$
100	3	−0.1	1.1	−5.6	−8.8	818.1	10,000.0	−46.7
	6	−1.8	0.5	−9.7	−11.5	1,250.8	10,000.0	−40.6
150	3	−4.3	−1.2	−18.3	−16.9	1,620.5	10,000.0	−44.1
	6	−9.3	−4.3	−21.2	−16.2	2,174.3	10,000.0	−41.6

In line with the results for Dataset S, ${\Delta }_{\mu }$ remains consistently below 0.8% for Dataset M, as well as exhibiting consistently negative Δ_L values between −0.3% to nearly −10%. The difference between Δ_L and ${\Delta }_{\mu }$ shows that even with the additional computational effort granted to the exact optimization approach, the solver achieved only marginal improvements over the matheuristic solutions. For example, for instances having 75 nodes under $\mathcal{T}=$ 3, the matheuristic on average required about 700 s. Limited to the same amount of time, the exact approach yielded solutions that were 6.0% inferior, on average. In the more than 9,000 s permitted thereafter to the exact optimization approach, the new solution improved upon that yielded by the matheuristic by 0.8%, on average. Indeed, the difference between Δ_L and ${\Delta }_{\mu }$ suggests that most of the computational effort expended by the exact optimization approach improves the incumbent lower bound, rather than discovering an improved feasible solution. To that end, negative ${\text{LB}}_{\mu }$ values between −0.6 to −10.6% demonstrate the tight lower bounds on the matheuristic’s solutions.

For instances having 100 or more nodes, the quality of both upper and lower bounds generated by the exact optimization approach greatly deteriorated. With increasingly negative values of ${\Delta }_{\mu }$ and ${\text{LB}}_{\mu },$ Table 5 indicates a growing gap between the matheuristic solution and the upper and lower bounds, respectively, yielded by the exact optimization approach in 10,000 CPU seconds. Moreover, the negative values in the Δ_t column affirm that as instance size scales, the exact optimization approach becomes vastly less effective compared with the matheuristic. Finally, the reported ${\Delta }_{\text{MKP}}$ values from the largest instances conform to the same trend from Datasets S and M, which further evidences the efficacy of the matheuristic.

The ${\Delta }_{\mathrm{\text{max}}}$ column shows at least one instance (in Table 4 with $\mathcal{T}=3$ under n = 75) where the gap between the matheuristic and exact solutions was as large as 6.4%. But that instance is an outlier: of the 210 instances in the test-bed, the matheuristic matched or improved upon the exact solution in 66.1% of instances ($\Delta \le$ 0%); 90.9% of instances exhibited $\Delta \le$ 1%; and all but 3/210 instances exhibited $\Delta \le$ 2.5%.

5.3. Matheuristic progression & solution insights

Although the matheuristic did not always yield optimal solutions, for a large number of instances the procedure yielded solutions that were very nearly optimal. For example, Figure 6a and b summarize results from instances having 40 nodes under $\mathcal{T}=$ 3, for which the matheuristic yielded 10/25 optimal solutions. Overall, however, the average deviation between the matheuristic and optimal solutions was 0.2% over the same instances. Moreover, Figure 6a, which tracks the overlap in customer assignments in each phase as compared with the optimal solution, shows that the overlap remained below 90% in the final phase. For the delivery profiles explored in this work, this pattern suggests the presence of consistent near-optimal customer assignments, even when sub-optimal, may yet yield high quality routing decisions and overall near-optimal solutions.

Figure 6. For each phase of the matheuristic, Figure 6a tracks the overlap of customers assigned to the drone in the matheuristic solution (as compared to the optimal solution). Meanwhile Figure 6b compares the quality of the incumbent routing solution to the optimal routes possible under that set of assignments. (a) Quality of Customer Assignments by Phase (b) Quality of Routing Assignments by Phase.

In the earliest phase of the matheuristic, sub-optimal customer assignments do adversely affect the quality of the solution. For example, Figure 6a plots customer assignment overlap (with the optimal solution) near 78% in the initialization phase. Correspondingly, the heuristic yielded solutions that were up to 14% sub-optimal on average during that phase, as shown in Figure 6b. In the first iteration of Phase 1, however, the customer assignment overlap worsened to 74%, even as the corresponding solution quality improved to be 2% sub-optimal on average. This trend suggests that a sub-optimal set of customer assignments may be amenable to a routing completion that yields a near-optimal solution for the instance.

These figures highlight the importance of the interplay between phases, and by extension the search criteria they use, in yielding high-quality solutions efficiently. As alluded to in Section 4.4, the needs of the decision-maker may dictate the sequence and number of iterations allotted to each phase of the heuristic. In our computational experience, two main loops through the matheuristic were the most effective setting. For larger instances especially, the second main loop served to significantly improve solution quality. At the same time, the matheuristic’s design leaves for users the freedom to terminate the search according to their own needs for solution quality and computational efficiency.

Several trends persisted in the best-known solutions yielded by the matheuristic: a relatively large share of customers served by drone; and relatively short waiting time by the vehicle (as captured by the $\widehat{f}$-variables, defined as the vehicle’s idle time during a drone flight). Across the S, M, and L datasets, on average there were 37.0%, 42.8%, and 45.0% of customers served by the drone, respectively. Those averages approach the upper bound of 50% of customers being served by the drone, which is imposed by the problem statement. Moreover, the average vehicle waiting time, expressed as a share of the objective value (essentially the share of the route’s duration that was spent idling) was 1.5%, 2.7%, and 3.4% across the datasets in the aforementioned order, respectively. Collectively, these trends indicate for the datasets in this study large networks may benefit from drone service, and that as the share of customers served by the drone increases, the vehicle idling time also increases, but remains a small share of less than 4% of the overall completion time of deliveries.

5.4. Carrier consistency insights

As previewed by the title and introduction of this article, serving customers consistently by the same carrier carries a logistical cost. Indeed, there are at least two ways in which consistent service may increase the logistical cost of drone delivery:

1. Carrier Consistency: Either the drone or the vehicle must exclusively deliver to a customer across the multi-day horizon;

2. Retention Consistency: As alluded to in Section 1, several drone delivery firms offer their service progressively, starting with a small group of customers that eventually expands. Retention consistency ensures that any customers served by the drone would retain service by that carrier, even as that service expands to include new customers.

Section 1 outlines several reasons, including the requirement for drone deliveries to be attended and received via a smartphone application, as to why the first type of consistency discussed above is a key feature of the TSPSBO. This discussion focuses on the logistical implications of the second type of consistency, in which a firm progressively increases the share of customers, Γ, for whom drone delivery is offered. Suppose for example that a decision-maker plans to offer drone delivery to 30% of customers in a particular network ultimately. To begin, the decision-maker may elect to offer the service to only 10% of customers. In other words, a decision-maker may plan to set $\Gamma =0.3$ eventually, but set $\Gamma =0.1$ in an early adoption period to ensure the reliability of the service first. In this case, which customers should be offered drone service first?

This study presents logistical insights associated with two strategies to address this critical question:

• Bounded Strategy: As Γ (the share of customers for whom the firm elects to provide drone delivery coverage) increases, a decision-maker is free to offer or to renege drone delivery service for customers. Given Γ, the decision-maker selects an optimal set of customers for assignment to the drone, without regard for whether they had been offered drone delivery previously.

• Decrementing Strategy: As Γ increases, a decision-maker may only incorporate new customers into drone delivery service. The implication is that once a customer has been offered drone delivery, their assignment to that carrier may not be revoked even if Γ later increases.

The Decrementing Strategy is so named because it begins with the projected maximal value of Γ chosen by the decision-maker, then moves backwards to designate the customers assigned to the drone for smaller values of Γ. Decrementing is comparable, mathematically, to appending an additional constraint to fix certain period-independent carrier assignments across repeated runs of the formulation. In other words, the Bounded Strategy is a relaxation of the Retention Consistency requirement in the Decrementing Strategy. Table 6 below compares the performance of the two strategies across values of Γ, starting with $\Gamma =0.1$ and proceeding incrementally until $\Gamma =0.5.$ The objective values and computational effort correspond to solutions from Model TSPSBO, which was used to solve instances from Dataset S. For each instance size, both the average objective value and computational effort are reported for the Bounded Strategy. For the Decrementing Strategy, the relative percent change from the Bounded Strategy values are shown. For example, for instances having 30 nodes under $\Gamma =0.10,$ the Decrementing strategy yielded objective values that were 0.4% higher, on average, than those of the Bounded Strategy.

Table 6. A comparison of managerial strategies for consistent drone service.

Bounded strategy
	$\Gamma =0.10$		$\Gamma =0.20$		$\Gamma =0.30$		$\Gamma =0.40$		$\Gamma =0.50$
n	$Z^{*}(\%)$	CPU (s)	$Z^{*}(\%)$	CPU (s)	$Z^{*}(\%)$	CPU (s)	$Z^{*}(\%)$	CPU (s)	$Z^{*}(\%)$	CPU (s)
20	100.0	1.4	100.0	1.8	100.0	1.5	100.0	1.9	100.0	1.4
30	100.0	14.6	100.0	22.6	100.0	36.7	100.0	26.8	100.0	23.8
40	100.0	338.9	100.0	753.3	100.0	813.9	100.0	804.8	100.0	475.8
Decrementing strategy
	$\Gamma =0.10$		$\Gamma =0.20$		$\Gamma =0.30$		$\Gamma =0.40$		$\Gamma =0.50$
n	$Z^{*}(\%)$	CPU (s)	$Z^{*}(\%)$	CPU (s)	$Z^{*}(\%)$	CPU (s)	$Z^{*}(\%)$	CPU (s)	$Z^{*}(\%)$	CPU (s)
20	100.1	0.5	100.0	0.6	100.0	0.4	100.1	0.3	100.0	1.6
30	100.4	5.1	100.2	2.8	100.2	2.4	100.0	1.2	100.0	23.5
40	100.2	34.0	100.3	38.8	100.2	17.6	100.0	9.8	100.0	473.4

Table 6 shows that for larger values of Γ, the gap in both objective values and computational effort between the two strategies disappears. That is because at $\Gamma =0.5,$ the Retention Consistency constraints become non-binding, since no more than 50% of customers may be served by the drone due to the lack of available rendezvous pairs. For smaller Γ values, however, there is a small “cost” of carrier consistency, no larger than 0.5%, for retaining consistent service for drone customers in the Decrementing Strategy. At the same time, the Decrementing Strategy is computationally advantageous, as might be expected. For instances having 40 nodes, for example, the Decrementing Strategy required an average of 17.6 CPU seconds for a solution under $\Gamma =0.3,$ while the Bounded Strategy required more than 800 CPU seconds, on average. For a marginal logistical cost as compared with the Bounded Strategy, decision-makers can thus offer customers both carrier and retention consistency through the Decrementing Strategy, with the added benefit of computational efficiency for the decision-maker.

6. Conclusion

This work investigates a last-mile delivery problem for subscription-based orders in which a firm must consistently assign customers to deliveries by either a vehicle or its companion drone, and to determine accompanying vehicle-drone routing decisions across a planning horizon. The proposed novel MIP formulation optimally solves instances having up to 40 customers over 6 days. For larger instances having up to 150 nodes over the same horizon, a proposed matheuristic tractably and consistently yields near-optimal solutions. The matheuristic exploits two aspects of the problem: the co-expression of customer deliveries and the geo-spatial dispersion of customer locations. From generating an initial solution to the subsequent cluster search procedures, each phase of the matheuristic is guided by analytical features of the problem that must be carefully balanced in a high-quality solution. The availability of rendezvous pairs for drone launch/recovery as well as the trade-offs involved in assigning carriers to dense customer clusters are all intertwined features of the TSP-SBO analyzed by the matheuristic. Based on these analytics, the matheuristic targets certain customer clusters for re-optimization using a variety of iterations of the proposed formulation.

In our computational study, the proposed matheuristic yielded optimal solutions for 110/150 instances with $n\le 40$ customers and up to 6 time-periods; for instances with $50\le n\le 75$ customers, it remarkably produced solutions proven by both upper- and lower-bounding strategies to be near-optimal; and for $n\ge 100$ customers, it produced comparable statistical patterns, suggesting high-quality solutions in manageable computational efforts (within 1,800 s at most). For the largest instances, the solver’s solutions to the proposed MIP within a time limit of 10,000 s were non-competitive.

This study advocates for the use of data and decision analytics jointly, in a manner that exploits the underlying data and problem features to guide algorithmic frameworks. In this context, our proposed MIP is first used as a standalone methodology for small-sized instances and then as a backbone for our matheuristic where cluster searches are driven by subscription patterns and location analytics. In light of the growing interest and need for unmanned vehicles such as aerial drones, this confluence of analytics and computational optimization is, in our opinion and experience, a promising research paradigm that may deliver numerical performance, and insights into managerial policies, for challenging logistical problems. A dynamic setting for this problem may be of future interest, to address the inclusion of new customers into the schedule and delivery plans.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are freely available in Gituhub at https://bit.ly/3l2phXV.

Notes

1 All instances as well as best-known solutions available at: https://bit.ly/3l2phXV

Appendix A.

Derivation of tightened big “M” values

In what follows, valid values are derived for the $M_{\mathit{\text{ij}}}^t$ parameter deployed in Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) in Section 3.2: (1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}\widehat{f_{\mathit{\text{ikj}}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in R^t\backslash \{n+1\},\\ j\in R^t\backslash \left\{0\right\}|i\ne j\end{array}$$(1i)

In each time period $t\in T,$ Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) sets the departure time, $b_j^t,$ from a node $j\in R^t$ visited by the vehicle ($\sum _{i\in R^t}{x}_{\mathit{\text{ij}}}^t=1$). If j is not visited by the vehicle, $M_{\mathit{\text{ij}}}^t$ must be large enough to retain the validity of Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) . The following inequalities determine the earliest and latest departure times for customer nodes (i.e., excluding the depot) visited by the vehicle, respectively: (1) $$b_i^t\ge {\tau }_{0i},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in C^t,$$(1) (2) $$b_i^t\le D_{\mathrm{\text{max}}}^t-{\tau }_{\mathit{\text{in}}+1},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in C^t$$(2)

Inequality (1)(1) $$b_i^t\ge {\tau }_{0i},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in C^t,$$(1) holds by construction, since the earliest arrival time for a customer i visited by the vehicle is the duration of a trip from the depot to that node (${\tau }_{0i}$). For Inequality (2)(2) $$b_i^t\le D_{\mathrm{\text{max}}}^t-{\tau }_{\mathit{\text{in}}+1},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in C^t$$(2) , $D_{\mathrm{\text{max}}}^t$ is the objective value associated with a feasible tour for period $t\in T.$ The vehicle’s latest departure from a customer node i may not exceed the difference between $D_{\mathrm{\text{max}}}^t$ and the vehicle’s travel time from that node to the depot, ${\tau }_{\mathit{\text{in}}+1}.$ As specified in Section 3.1, nodes 0 and n + 1 represent source and sink nodes, with shared coordinates at the central depot location. We also assume symmetrical inter-node distances.

To scale $M_{\mathit{\text{ij}}}^t,$ suppose $x_{\mathit{\text{ij}}}^t=x_{\mathit{\text{ji}}}^t=0,$ which implies $\sum _{k\in (C^t\cap D)}\widehat{f_{\mathit{\text{ikj}}}}{y}_{\mathit{\text{ikj}}}^t=0$ by Constraint (1f)(1f) $$\sum _{(i,j)\in E_k^t}{y}_{\mathit{\text{ikj}}}^t+\sum _{i\in R^t}{x}_{\mathit{\text{ik}}}^t=1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,k\in C^t,$$(1f) . Thus Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) simplifies to: $$b_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}-M_{\mathit{\text{ij}}}^t\Rightarrow M_{\mathit{\text{ij}}}^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}-b_j^t,$$

Consider the following cases:

1. Node i is the depot

1. Whenever the vehicle departs the source node representing the depot, its departure time is equal to 0 ($b_0^t=0$ $\forall t\in T$). With i = 0, simplifying Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) yields $$M_{0j}^t\ge {\tau }_{0j}-b_j^t$$ For a customer node j, by Inequality (1)(1) $$b_i^t\ge {\tau }_{0i},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in C^t,$$(1) , the earliest departure time from j must be at least equal to the duration of the vehicle’s trip to j, $b_j^t\ge {\tau }_{0j}.$ Further simplifying Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) yields $$M_{0j}^t\ge 0$$

2. Since the vehicle should never depart the sink node representing the depot, Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) ensures $i\ne n+1.$ Therefore this case is irrelevant for the problem at hand.

2. Node j is the depot

1. Since the vehicle’s tour should begin at the source node representing the depot, Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) ensures $j\ne 0.$ Therefore this case is irrelevant for the problem at hand.

2. At the sink node representing the depot, the vehicle’s “departure” time is in fact the duration of a feasible tour for that period ($b_{n+1}^t=D_{\mathrm{\text{max}}}^t$ $\forall t\in T$). With $j=n+1,$ simplifying Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) yields

$$M_{\mathit{\text{in}}+1}^t\ge b_i^t+{\tau }_{\mathit{\text{in}}+1}-D_{\mathrm{\text{max}}}^t$$

Given the duration of a feasible tour, for any customer node $i\in V^{+},$ the latest time by which the vehicle must depart i is $b_i^t\le D_{\mathrm{\text{max}}}^t-{\tau }_{\mathit{\text{in}}+1}$ according to Inequality (2)(2) $$b_i^t\le D_{\mathrm{\text{max}}}^t-{\tau }_{\mathit{\text{in}}+1},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in C^t$$(2) . With $i\in V^{+}$ and $j=n+1,$ simplifying Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) yields $$M_{\mathit{\text{in}}+1}^t\ge D_{\mathrm{\text{max}}}^t-{\tau }_{\mathit{\text{in}}+1}+{\tau }_{\mathit{\text{in}}+1}-D_{\mathrm{\text{max}}}^t\Rightarrow M_{\mathit{\text{in}}+1}^t\ge 0$$

• Nodes $i,j\in V^{*}$: To assume the worst-case bound, let the departure time from node i be as large as possible, and the departure time from node j be as small as possible. Then by Inequalities (1)(1) $$b_i^t\ge {\tau }_{0i},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in C^t,$$(1) and (2)(2) $$b_i^t\le D_{\mathrm{\text{max}}}^t-{\tau }_{\mathit{\text{in}}+1},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in C^t$$(2) , $b_i^t=D_{\mathrm{\text{max}}}^t-{\tau }_{\mathit{\text{in}}+1}$ and $b_j^t={\tau }_{0j}.$ With both i and j representing customer nodes, simplifying Constraint (1i)(1i) $$\begin{array}{c}{b}_j^t\ge b_i^t+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^t\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^t-M_{\mathit{\text{ij}}}^t(1-x_{\mathit{\text{ij}}}^t),\mathrm{\hspace{1em}}\forall t\in T,\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^t\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^t\backslash \left\{0\right\}|i\ne j,\end{array}$$(1i) yields

$$M_{\mathit{\text{ij}}}^t\ge D_{\mathrm{\text{max}}}^t-{\tau }_{\mathit{\text{in}}+1}+{\tau }_{\mathit{\text{ij}}}-{\tau }_{0j},$$

The valid values of M are as follows: $$M_{\mathit{\text{ij}}}^t=D_{\mathrm{\text{max}}}^t-{\tau }_{\mathit{\text{in}}+1}+{\tau }_{\mathit{\text{ij}}}-{\tau }_{0j},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall t\in T,i\in R^t\backslash \{n+1\},j\in R^t\backslash \left\{0\right\}|i\ne j\hspace{1em}\square$$

Appendix B.

Single-period variant of model TSPSBO

This section introduces a single-period variant of the formulation from Section 3, which underpins key stages of the matheuristic. For convenience, the parameter and variable definitions from the original formulation in Section 3.3 are reproduced below. A key distinction between the single-period variant and the original model is that variables from the latter affect decisions across multiple periods (i.e., they are period independent) and are therefore parameterized in this single-period variant of the model.

In what follows, we refer to the single period variant of Model TSPSBO as Model TSPSBO-SP. Without loss of generality, we designate $\overline{t}\in T$ as the single period of reference for Model TSPSBO-SP. Finally, parameterized period-independent variables are denoted by a bar (e.g., $\overline{\delta }$) for clarity.

Input parameters

• $T=\{1,\dots ,\mathcal{T}\}:$ Set of service periods.

• $V=\{0,1,\dots ,n,n+1\}:$ Set of service nodes in the network.

• $V^{+}=V\backslash \{0,n+1\}:$ Set of customer nodes, which excludes the central depot (represented by dummy nodes 0 and n + 1).

• $A=\{(i,j):i,j\in V,i\ne j\}:$ Set of arcs connecting distinct nodes in the network.

• α: Ratio of drone to vehicle speed.

• τ_ij: Vehicle travel time for $(i,j)\in A.$

• f_ij: Drone flight time for $(i,j)\in A.$

• $d_j^{\overline{t}}:$ Binary indicator such that $d_j^{\overline{t}}=1\iff$ customer j requires parcel delivery in period $\overline{t},\hspace{0.17em}\forall j\in V^{+}.$

• F: Maximum flight duration for the drone to make a single delivery then return to the vehicle.

• $C^{\overline{t}}$ $=\{j\in V^{+}$ $|d_j^{\overline{t}}=1\}:$ Set of customers requiring parcel delivery in period $\overline{t}.$

• $R^{\overline{t}}=C^{\overline{t}}\cup \{0,n+1\}:$ Set of nodes requiring carrier visits in period $\overline{t}.$

• ${\widehat{f}}_{\mathit{\text{ikj}}}=\mathrm{\text{max}}(0,f_{\mathit{\text{ik}}}+f_{\mathit{\text{kj}}}-{\tau }_{\mathit{\text{ij}}}):$ The vehicle’s waiting time for a drone flight originating from i, delivering to k, and landing at j, $\forall$ $(i,j)\in A,k\in D,$ with distinct i, j, k.

• $E_k^{\overline{t}}=\{(i,j)\in R^{\overline{t}}$ $|f_{\mathit{\text{ik}}}+f_{\mathit{\text{kj}}}\le F,$ $i,j\in R^{\overline{t}}$$\}:$ All pairs of rendezvous pairs permitting a drone delivery to customer k in period $\overline{t},$ with distinct i, j, k.

• $D=\{j\in V^{+}$ $|$ $d_j^{\overline{t}}>0$ $\wedge$ $\left|E_j^{\overline{t}}\right|>0$$\}:$ Set of nodes for which drone delivery is feasible by virtue of having rendezvous pairs available in each period in which the customer requires delivery.

• $D_{\mathrm{\text{max}}}^{\overline{t}}:$ The objective value of a feasible solution in period $\overline{t}$ in which all nodes $i\in R^{\overline{t}}$ are served.

• $M_{\mathit{\text{ij}}}^{\overline{t}}:$ A sufficiently large scalar associated with period $\overline{t},(i,j)\in A.$

• $\overline{{\delta }_j}\in [0,1]:$ $\overline{{\delta }_j}=1$ $\iff$ the drone delivers to customer j across all periods in which j expresses demand, $\forall j\in V^{+}.$ Since δ is a period independent variable in the multi-period formulation, it is treated as a parameter in this single-period variant of the formulation.

Decision variables

• $x_{\mathit{\text{ij}}}^{\overline{t}}\in \{0,1\}:$ $x_{\mathit{\text{ij}}}^{\overline{t}}=1$ $\iff$ the vehicle travels from i to j, $\forall$ $i,j\in R^{\overline{t}},i\ne j.$

• $y_{\mathit{\text{ikj}}}^{\overline{t}}\in \{0,1\}:$ $y_{\mathit{\text{ikj}}}^{\overline{t}}=1$ $\iff$ the drone launches from i, delivers to k, and lands at j, $\forall$ $k\in C^{\overline{t}},\hspace{0.17em}(i,j)\in E_k^{\overline{t}},$ with distinct i, k, j.

• $b_j^{\overline{t}}\in [0,D_{\mathrm{\text{max}}}^{\overline{t}}]:$ The departure time of the vehicle at node j, for $j\in R^{\overline{t}}.$

The Single-Period variant of the formulation for the Traveling Salesman Problem with Subscription-Based Orders, Model TSPSBO-SP, is presented as follows: (B.1a) $$\mathbf{\text{TSPSBO}}\mathbf{-}\mathbf{\text{SP}}:\text{Minimize}\mathrm{ }\hspace{1em}\sum _{i\in R^{\overline{t}}}\sum _{j\in R^{\overline{t}}|j\ne i}{\tau }_{\mathit{\text{ij}}}{x}_{\mathit{\text{ij}}}^{\overline{t}}+\sum _{k\in (C^{\overline{t}}\cap D)}\sum _{(i,j)\in E_k^{\overline{t}}}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{ijk}}^{\overline{t}}$$(B.1a) (B.1b) $$\mathrm{s}.\mathrm{t}.\hspace{1em}\sum _{j\in C^{\overline{t}}}{x}_{0j}^{\overline{t}}=1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}$$(B.1b) (B.1c) $$\sum _{j\in C^{\overline{t}}}{x}_{\mathit{\text{jn}}+1}^{\overline{t}}=1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}$$(B.1c) (B.1d) $$\sum _{j\in R^{\overline{t}}}{x}_{\mathit{\text{ij}}}^{\overline{t}}-\sum _{j\in R^{\overline{t}}}{x}_{\mathit{\text{ji}}}^{\overline{t}}=0,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}}\forall i\in C^{\overline{t}},$$(B.1d) (B.1e) $$\sum _{j\in R^{\overline{t}}}{x}_{\mathit{\text{kj}}}^{\overline{t}}=1-\overline{{\delta }_k},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}}\forall k\in C^{\overline{t}},$$(B.1e) (B.1f) $$\sum _{(i,j)\in E_k^{\overline{t}}}{y}_{\mathit{\text{ikj}}}^{\overline{t}}+\sum _{i\in R^{\overline{t}}}{x}_{\mathit{\text{ik}}}^{\overline{t}}=1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall k\in C^{\overline{t}},$$(B.1f) (B.1g) $$\begin{array}{c}{y}_{\mathit{\text{ikj}}}^{\overline{t}}\le x_{\mathit{\text{ij}}}^{\overline{t}},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}\forall k\in (C^{\overline{t}}\cap D),\\ \mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{\hspace{1em}}(i,j)\in E_k^{\overline{t}},\end{array}$$(B.1g) (B.1h) $$\sum _{k\in C^{\overline{t}}|(i,j)\in E_k^{\overline{t}}}{y}_{\mathit{\text{ikj}}}^{\overline{t}}\le 1,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\forall i,j\in R^{\overline{t}},$$(B.1h) (B.1i) $$\begin{array}{c}{b}_j^{\overline{t}}\ge b_i^{\overline{t}}+{\tau }_{\mathit{\text{ij}}}+\sum _{k\in (C^{\overline{t}}\cap D)}{\widehat{f}}_{\mathit{\text{ikj}}}{y}_{\mathit{\text{ikj}}}^{\overline{t}}-M_{\mathit{\text{ij}}}^{\overline{t}}(1-x_{\mathit{\text{ij}}}^{\overline{t}}),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall i\in R^{\overline{t}}\backslash \{n+1\},\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\forall j\in R^{\overline{t}}\backslash \left\{0\right\}|i\ne j,\end{array}$$(B.1i) (B.1j) $$b_{n+1}^{\overline{t}}\le \sum _{i\in R^{\overline{t}}}\sum _{l\in R^{\overline{t}}}{\tau }_{\mathit{\text{il}}}{x}_{\mathit{\text{il}}}^{\overline{t}}+\sum _{k\in (C^{\overline{t}}\cap D)}\sum _{(i,l)\in E_k^{\overline{t}}}{\widehat{f}}_{\mathit{\text{ikl}}}{y}_{\mathit{\text{ikl}}}^{\overline{t}},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}$$(B.1j) (B.1k) $$b_{n+1}^{\overline{t}}\le D_{\mathrm{\text{max}}}^{\overline{t}},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}$$(B.1k) (B.1l) $$b\ge 0\hspace{1em}x,y\mathrm{ }\text{binary}.\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}$$(B.1l)

Appendix C.

Demand profile generation for instances

Table C1 summarizes commercial parcel delivery trends for households around the world, some of which may be too heavy or bulky for drone delivery. As outlined in Section 1, subscriptions to different product types may lead households to receive more or fewer parcels than the average in their country; creating three separate customer categories helps account for this reality. Specifically, customer categories A and B, which collectively comprise 40% of all customers, receive subscriptions with an average 50% likelihood, which matches the upper range of average daily parcel deliveries in countries like the United States, the United Kingdom, and Japan from Table C1. The category C customers, whose subscriptions recur with 10% likelihood, account for subscriptions which may require far less frequent deliveries.

Table C1. International parcel delivery statistics, 2020.

Country	Parcels^a (Billions)	Households (Millions)	average parcels delivered per household (Annually)	Average parcels delivered per household (Daily)
Japan	9.1	53.33^b	168.8	0.47
United Kingdom	5.0	27.79^c	136.6	0.49
United States	20.2	120.75^d	122.8	0.46
Australia	1.0	7.76^e	120.4	0.35
Germany	4.1	41.51^f	89.1	0.27

^aSpadafora (2021).

^bJapanese Statistics Bureau (2019).

^cBritish Office for National Statistics (2021).

^dUnited States Census Bureau (2019).

^eAustralian Institute of Family Studies (2019).

^fGerman Destatis (2019).