<figure><img src="image_1.png" /></figure>

HAL
open science

# A Branch-Cut-And-Price approach for the Two-Echelon Vehicle Routing Problem with Drones

Sylvain Lichau, Ruslan Sadykov, Julien François, Rémy Dupas

► To cite this version:

Sylvain Lichau, Ruslan Sadykov, Julien François, Rémy Dupas. A Branch-Cut-And-Price approach for the Two-Echelon Vehicle Routing Problem with Drones. 2024. hal-04584429

HAL Id: hal-04584429

https://hal.science/hal-04584429

Preprint submitted on 23 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# A Branch-Cut-And-Price approach for the Two-Echelon Vehicle Routing Problem with Drones

Sylvain Lichau, Ruslan Sadykov, Julien François, Rémy Dupas

► To cite this version:

Sylvain Lichau, Ruslan Sadykov, Julien François, Rémy Dupas. A Branch-Cut-And-Price approach for the Two-Echelon Vehicle Routing Problem with Drones. Computers and Operations Research, 2025, 173, pp.106869. 10.1016/j.cor.2024.106869. hal-04584429v2

HAL Id: hal-04584429

https://hal.science/hal-04584429v2

Submitted on 7 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# A Branch-Cut-And-Price approach for the Two-Echelon Vehicle Routing Problem with Drones

**Sylvain Lichau**

Univ. Bordeaux, CNRS, Bordeaux INP, IMS, UMR 5218, F-33400 Talence, France,
Email: sylvain.lichau@u-bordeaux.fr

**Ruslan Sadykov**

Univ. Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251, F-33400 Talence, France

**Julien François**

Univ. Bordeaux, CNRS, Bordeaux INP, IMS, UMR 5218, F-33400 Talence, France,

**Rémy Dupas**

Univ. Bordeaux, CNRS, Bordeaux INP, IMS, UMR 5218, F-33400 Talence, France

June 2023

## Abstract

In this paper, we propose a new set-partitioning model for the two-echelon vehicle routing problem with drones (2E-VRP-D), where partial routes corresponding to drone movements are enumerated using an efficient dynamic program. To solve the model we use an exact branch-cut-and-price algorithm, and a labelling algorithm for the pricing problem both based on state-of-the art literature. We also propose an adaptation of the well-known rounded capacity cuts for this problem, as well as pre-processing methods to reduce the size of the problem. In addition, this paper presents an effective heuristic branch-cut-and-price, based on the exact branch-cut-and-price algorithm. Extensive computational experiments are presented, showing that the exact algorithm can solve all the instances from the literature on the 2E-VRP-D, and almost multiply by four the size of the clustered instances solved optimally. Sensitivity analysis is also conducted for the proposed improvements.

# 1 Introduction

In recent years, the explosion in e-commerce and the ever-growing concentration of populations in cities have led to a drastic increase of demand for home deliveries, overloading last mile logistics. To overcome this challenge, incorporating drones in deliveries is a promising solution. In addition to paralleling deliveries, enabling customers to receive their packages faster, using drones would also decrease costs for the provider (Wang, 2015). It would also decrease the pollution emitted by truck routes, drones being less polluting than trucks. By reducing
# A Branch-Cut-And-Price Approach for the Two-Echelon
Vehicle Routing Problem with Drones

Sylvain Lichau<sup>a,b,*</sup>, Ruslan Sadykov<sup>b,d</sup>, Julien François<sup>c</sup>, Rémy Dupas<sup>c</sup>

<sup>a</sup>Univ. Bordeaux, CNRS, Bordeaux INP, IMS, UMR 5218, Talence, F-33400, France

<sup>b</sup>Univ. Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR
5251, Talence, F-33400, France

<sup>c</sup>Univ. Bordeaux, CNRS, Bordeaux INP, IMS, UMR 5218, Talence, F-33400, France

<sup>d</sup> Atopima, 16 place Sainte-Eulalie, Bordeaux, F-33000, France

## Abstract

We present a new set-partitioning model for the two-echelon vehicle routing problem with drones (2E-VRP-D). In this model, partial routes corresponding to drone movements are enumerated by an efficient dynamic program. A state-of-the-art exact branch-cut-and-price algorithm solves the model using a labeling algorithm for the pricing problems. We further propose an adaptation of the well-known rounded capacity cuts, specifically tailored to address this problem, as well as pre-processing methods designed to enhance the algorithm's efficiency. In addition, we propose an effective heuristic branch-cut-and-price based on the exact algorithm. Computational experiments demonstrate that the exact algorithm solves all instances from the 2E-VRP-D literature and significantly increases the size of optimally solved clustered instances. Finally, a sensitivity analysis is conducted to assess the impact of the proposed enhancements.

## 1. Introduction

In recent years, the quick growth of e-commerce and the continuous urbanization of populations have resulted in a significant increase in demand for home delivery services, overloading last-mile logistics. To address this challenge, the integration of drones in delivery operations presents a promising solution. In addition to paralleling deliveries, allowing customers to receive their packages faster, the use of drones would also decrease costs for the provider (Wang, 2015). The pollution emitted by the deliveries would diminish because drones are less polluting than trucks. Moreover, traffic congestion could be alleviated substantially by decreasing the volume of truck-based deliveries in urban areas (Figliozzi, 2017).

Numerous pilot projects involving drone deliveries have been successful, and

<footnote>*Corresponding author
Email address: sylvain.lichau@u-bordeaux.fr (Sylvain Lichau)</footnote>
the number of customers delivered by trucks in cities, and in particular on small roads, the traffic congestion would decrease significantly.

Numerous drone delivery pilot projects have been successfully completed and many delivery companies are trying to integrate drones to deliveries. One of the first successful implementation of a couple drone/truck delivery, to our knowledge, has been established in France in 2016, to a difficult-to-reach area (DPDgroup, 2016). In urban area, in 2017 a pilot project named Vans & drones was conducted in Zürich (Roca-Riu and Menendez, 2019), which combined drones and vans to test on-demand delivery. In 10 days, 50 packages were successfully delivered, and the delay was always lower than two hours. Swiss Post has permanently established the use of drones to transfer blood samples between hospitals in Lugano, Zürich and Bern (Post, 2023). More recently, Amazon is claiming that they will start using drone deliveries in Italy, the UK and the US by the end of 2024 (Amazon, 2023).

Despite the positive impact of drone-assisted delivery and the existence of some real-word applications, certain difficulties remain in terms of implementation as it might bring this new technology very close to the public. Drone use raises security issues as the drone must be completely foolproof and ready to face any problem, such as extreme weather, technical default, breakdown, signal loss or even attacks. Drones are also subject to prohibited air zone and other legal restrictions (Lin et al., 2018). Finally, public risk beliefs, such as privacy invasions or injuring people due to malfunction, are also an important obstacle to the implementation of drone deliveries on large scales, in particular in urban areas (Zhu, 2019). Garg et al. (2023) provides a literature review over the potential benefits, challenges, and limitations of the usage of drones in last mile delivery.

Due to the potential evolution of public opinion, regulation and technology, it is of first importance to understand how deliveries with drones will impact cities and how to design them in the best possible way. In line with these future developments, the field of operations research is tackling several variants of drone delivery problems, as reviewed in (Macrina et al., 2020) and (Chung et al., 2020). Researchers are considering different truck-drone cooperation. In the classic vehicle routing problem with drones (VRP-D), each truck is paired with a single drone that delivers a single customer leaving and merging only when the truck is parked at a customer. Many different cooperation systems exist, leading to a large variety of problems.

In this paper we study the Two-Echelon Vehicle Routing Problem with Drones (2E-VRP-D). In this problem a set of vehicles and a set of drones have to serve a set of customers from a depot. The aim is to minimize the total duration of the routes. We have chosen to focus on this problem because it is sufficiently generic to correspond to several real-world situations, whereas encompassing the following set of requirements:
many logistics companies are actively working on integrating drones into their existing delivery systems. One of the first successful implementations of a drone/truck delivery system, to our knowledge, was established in France in 2016, specifically targeting a remote area (DPDgroup, 2016). In 2017, a pilot project named Vans & drones combining these two types of vehicles to test on-demand delivery in urban areas was conducted in Zürich (Roca-Riu and Menendez, 2019). Over a period of ten days, they successfully delivered 50 packages, with delays not exceeding two hours. Since 2017, Swiss Post has permanently deployed drones to transfer blood samples between hospitals in Lugano, Zürich, and Bern (Post, 2023). More recently, Amazon has announced the launch of drone deliveries in Italy, the UK, and the US by the end of 2024 (Amazon, 2023).

Despite the promising impact of drone-assisted delivery and the existence of real-world applications, several challenges arise from its large-scale adoption, particularly regarding the close interaction between this new technology and the citizens. Drone use raises security issues as drones must be completely fool-proof and ready to manage various problems, such as extreme weather, technical malfunctions, breakdowns, signal loss, or potential attacks. Drones are also subject to prohibited air zones and other legal restrictions (Lin et al., 2018). Finally, the public's risk beliefs, such as privacy invasions or injuring people due to malfunction, are also significant obstacles to large-scale implementations of drone deliveries. (Zhu, 2019). Garg et al. (2023) present a literature review of the potential benefits, challenges, and limitations of using drones in last-mile delivery.

Due to the potential evolution of factors such as public opinion, regulation, and technology, it is particularly important to understand how deliveries with drones will impact cities and how to design them in the best possible way. In line with these future developments, the field of operations research is tackling several variants of drone delivery problems, as reviewed in (Macrina et al., 2020) and (Chung et al., 2020). Researchers are considering different truck-drone co-operation systems. In the classic vehicle routing problem with drones (VRP-D), each truck is paired with a single drone. The truck must be parked for the drone to launch and land. The truck continues its route while the drone delivers packages to nearby customers. Many different cooperation systems exist, leading to a large variety of problems.

We study the Two-Echelon Vehicle Routing Problem with Drones (2E-VRP-D). In this problem, a set of vehicles and autonomous drones must deliver packages to customers from a depot. The objective is to minimize the total duration of the routes. This problem is general enough to apply to various real-world situations. It can be summarized as follows:

* Each truck carries an adjustable number of drones. Each drone and its necessary equipment reduce the associated truck payload capacity. Indeed, the technological gap between implementing a truck that carries one drone and a truck that carries several is relatively small, yet it significantly
* Each truck contains an adjustable number of drones. The useful capacity of the trucks depends on the number of carried drones, each drone and its necessary equipment having a weight. Carrying several drones on a truck rather than just one can be easily implemented, while having an important impact on the delivery.

* Full synchronization exists between the truck and the drone. Indeed, a model where drones can change truck en route will be hard to implement on a real-world application, for security and responsibility reasons. A drone can fly only when its assigned truck is stationed at a customer, so as to achieve a compromise between simplicity and efficiency in the proposed model. Furthermore, the drones can operate multiple round trips from the same customer. Choosing the number of drones carried for each vehicle enables us to balance the total number of drones, the weight taken by the drones in a truck and the efficiency of the route.

* Each drone has a limited capacity and can only carry one package per flight. One package per drone is realistic as it is more secure than allowing a drone to carry several packages at the same time. It also reduces the payload, enabling longer flight. With removable batteries, drones can be instantaneously recharged when reaching the truck. Carrying a sufficient number of batteries for deliveries does not entail any particular transport constraints. Static truck when drones are flying allows close range monitoring in case of problems. A battery limits the maximal flight distance, dependent of the payload. It is assumed that drones are fully recharged upon reaching the truck.

* Each customer is characterized by a demand, a deadline, its availability for drone delivery and a service time dependent on the means of transport. The service at a customer delivered by a vehicle happens in parallel with the deliveries made by drones from this customer.

This is not the usual cooperation system between trucks and drones, but it can be argued that it enables an operator to remain near the drones to monitor them, or in case of an accident. For these reasons we think that the 2E-VRP-D is a highly topical and realistic problem as it is easy to implement in the delivery process. It is related to the two-echelon vehicle routing problem (2E-VRP), see (Sluijk et al., 2023), as we can consider the deliveries of the trucks as a first echelon, and the deliveries of the drones as a second echelon. In this context, Zhou et al. (2023) recently proposed a branch-and-price approach that optimally solves instances of up to 30 customers in less than 3 hours for the 2E-VRP-D. Given that real-world last mile delivery problems commonly contain a higher number of customers to serve it is necessary to be able to solve larger instances in a reasonable time. In this paper we propose a branch-cut-and-price algorithm for the 2E-VRP-D using recent advances from the well-known capacitated vehicle routing problem (CVRP). The main contributions of this paper are as follows:
improves delivery efficiency.

* Each truck and its drone(s) are fully synchronized, i.e., drones are assigned to a truck. A model where drones could change trucks en route would be complex to implement on a real-world application for security and responsibility reasons. To achieve a compromise between implementation simplicity and delivery efficiency, a drone can only fly when its associated truck is parked at a customer's location. Furthermore, the drones can operate multiple round trips from the same parking spot. The selection of each vehicle's number of carried drones is made individually.

* Each drone has a limited capacity and, for security reasons, can carry only one package per flight. Drones are considered autonomous. Before each launch, a setup time is required to load the package and change the battery. Drones, utilizing removable batteries, are assumed to be fully recharged upon leaving the truck. Carrying sufficient batteries for deliveries does not introduce any significant transport constraints. A stationary truck during drone flights allows a quick reaction in case of an incident. The maximum flight distance is constrained by the drone's battery capacity and the weight of the payload.

* Each customer is characterized by a demand, a deadline, and a service time dependent on the means of transport (i.e., truck or drone). Some customers cannot be served by drone. As drones are considered autonomous, the driver can serve the customer where the truck is parked, while the drones independently make other deliveries from the same customer.

Therefore, the 2E-VRP-D is a highly relevant and pragmatic problem as it is relatively straightforward to implement in the delivery process. It is related to the two-echelon vehicle routing problem (2E-VRP), see (Sluijk et al., 2023), because truck deliveries can be considered as a first echelon, while drone deliveries represent a second echelon. In this context, Zhou et al. (2023) proposed a branch-and-price approach that optimally solves instances of up to 30 customers in less than 3 hours for the 2E-VRP-D. Given that real-world last-mile delivery problems commonly contain a higher number of customers to serve, it is necessary to be able to solve larger instances in a reasonable time. We propose a branch-cut-and-price algorithm for the 2E-VRP-D using recent advances from the Capacitated Vehicle Routing Problem (CVRP). The main contributions of this paper are the following:

* A new branch-cut-and-price approach with state-of-the-art components based on enumerating partial routes corresponding to drone operations in the subproblems. The enumeration is modeled as a scheduling problem and solved using an efficient dynamic program. An adaptation of the rounded capacity cuts strengthens the algorithm cuts to better suit this problem, by considering drones' weight in the truck. Additionally, we propose preprocessing methods to eliminate dominated and infeasible enumerated partial routes.
* A new branch-cut-and-price approach based on the enumeration of the partial routes corresponding to drone movements in the subproblems and modelled as scheduling problems, and solved using an efficient dynamic program and state-of-the-art components. The algorithm is strengthened by an adaptation of the rounded capacity cuts to better suit this problem, by taking into account the drone's weight in the truck. We also propose pre-processing methods to eliminate dominated and infeasible the enumerated partial routes to efficiently reduce the subproblems size without excluding optimal solutions.

* A heuristic branch-cut-and-price based on the restriction of the subproblem solution. This heuristic is particularly efficient for instances where customers are close to each other. In this setting, drones flying from a parked customer *c* can reach many customers. We say *c* has a large drone-neighborhood. Moreover, in this setting, the drone-neighborhoods are highly interconnected, meaning that a customer can be reached by drone by many different customers.

* Extensive experiments are conducted on the instances and show that our methods can solve clustered instances with up to 120 customers, increasing greatly the size of the solved instances from the literature. We also show the benefit of the improvements proposed. Six new sets of medium-size instances with 40, 100, 120, 130, 140 and 150 customers extending the dataset of small-size instances proposed in (Zhou et al., 2023) and based on (Chen et al., 2021), (Solomon, 1987) and (Gehring and Homberger, 1999).

In section 2, we present a literature review mainly focusing on exact methods to solve vehicle routing problems with drones. In section 3, we define the problem we are addressing. In section 4, we present the pricing problem and how we solve it. In section 5, we present our approach, providing details on the branch-and-cut algorithm and on improvements we use. In section 6, we present the heuristic branch-cut-and-price. Finally, in section 7, we provide computational experiments highlighting the strength of our method.

# 2 Litterature review

The vehicle routing problem with drones (VRPD) class of problem is a variant of the classic vehicle routing problem (VRP) that integrates drones in the deliveries. This section is separated in two parts, the first one will cover the variants of the travelling salesman problem (TSP) (i.e. single tour), and the second the variants of the VRP (i.e. multiple tours).

## 2.1 Drones integrated in TSP

The interest to incorporate drones into routing problem emerged with the work of Murray and Chu (2015), in which a drone and a truck cooperate to solve a
* A heuristic branch-cut-and-price based on the restriction of the subproblem solution space. This heuristic demonstrates particular efficiency in instances where customers are located randomly. In such cases, drones launched from a vehicle parked at a customer can serve a large number of nearby customers, creating a large "drone-reachable neighborhood".

* We introduce six new sets of medium-size instances with 40, 100, 120, 130, 140 and 150 customers extending the dataset of small-size instances proposed in (Zhou et al., 2023) and based on (Chen et al., 2021), (Solomon, 1987) and (Gehring and Homberger, 1999). Extensive experiments on the instances show that our methods can solve clustered instances with up to 120 customers, significantly expanding the size of the solved instances compared to the existing literature. We also conduct sensitivity analyses on the adapted rounded capacity cuts and the preprocessing methods.

In section 2, we present a literature review with a primary focus on exact methods for solving vehicle routing problems with drones. In section 3, we formally define the problem being addressed. In section 5, we present the pricing problem and propose our solution approach. In section 4, we present our methodology, detailing the branch-cut-and-price algorithm along with the proposed improvements. In section 6, we present the heuristic branch-cut-and-price. Finally, in section 7, we provide computational experiments highlighting the strength of our method.

## 2. Literature review

The vehicle routing problem with drones (VRPD) class of problem is a variant of the classic vehicle routing problem (VRP) that integrates drones into the deliveries. The first part of this section covers the variants with drones of the traveling salesman problem (TSP) (i.e., single tour), and the second part covers the variants with drones of the VRP (i.e., multiple tours).

### 2.1. Drones integrated in TSP

The interest to incorporate drones into the routing problem emerged with the work of Murray and Chu (2015), in which a drone and a truck cooperate to solve a TSP-like problem called the flying sidekick traveling salesman problem (FSTSP). They presented a MIP model and two heuristics to solve instances with up to 10 customers. Various extensions were proposed, as in (Agatz et al., 2018) that allowed trucks to wait for the drones at the nodes and proposed a MILP and two heuristics based on the route-first, cluster-second method to solve the problem. The authors also experimentally analyze the time saving achieved by using a truck paired with a drone versus a truck alone. Bouman et al. (2018) used a dynamic program to solve the TSP-D instances with up to 20 customers. They demonstrated that reducing the number of customers visited by truck during drone flights significantly reduced the computation time without diminishing the solution quality. Yang et al. (2023) consider the robust drone-truck
TSP-like problem called the flying sidekick travelling salesman problem (FSTSP). They presented a MIP model and two heuristics to solve the problem with up to 10 customers. Various extensions were proposed, as in (Agatz et al., 2018) that allowed trucks to wait for the drones at the nodes and proposed a MILP and two heuristics based on the route-first, cluster-second method to solve the problem. The authors also analyse experimentally the time saving achieved by the use of a truck paired with a drone versus a truck alone. Bouman et al. (2018) used a dynamic program to solve the TSP-D up to 20 customers. They also showed that decreasing the number of customers visited by a truck while the drone was flying significantly decreased the computation time, without decreasing the quality of the solution significantly. Yang et al. (2023) consider the robust drone-truck delivery problem, where a truck and a drone have to deliver customers to maximize a profit function under truck routes uncertainty. The uncertainty is considered on the duration of the truck arcs, due to network congestion. They propose a branch and price algorithm to solve instances with up to 40 customers.

To improve the efficiency of the cooperation even further, researchers have studied the use of several drones on the truck. In (Karak and Abdelghany, 2019) the authors propose the “Hybrid Vehicle-Drone Routing Problem” (HVDRP), where a single vehicle carries a fleet of drones performing pickup and delivery operations. They propose a MIP formulation along with 3 heuristics to solve the problem. Murray and Raj (2020) introduced the “Multiple Flying Sidekicks travelling Salesman Problem” (mFSTSP), which is an extension with several vehicles of the FSTSP. It considers an arbitrary number of heterogeneous drones that may be deployed from the depot or from the delivery truck. They provide a MIP formulation solving instances with at most 8 customers, as well as a heuristic approach. Kang and Lee (2021) introduce the Heterogeneous Drone-Truck Routing Problem (HDTRP), where drones with different speeds and batteries are setup on a truck, and can fly only when the truck is parked at a customer. They propose a MILP formulation as well as a branch-and-cut algorithm to solve up to 50 customers instances optimally.

## 2.2 Drones integrated in VRP

The use of drones with a set of trucks was first introduced in (Wang et al., 2017), which presents the Vehicle Routing Problem with Drones (VRPD), where multiple trucks are equipped with drones to serve a set of customers. In this problem, the drones carried by a truck can deliver one parcel per flight, and the truck can move between the launch and the retrieval of its drones. The authors analyse different worst-case scenarios and provide upper bounds on the savings induced by the use of drone. The work is extended in (Wang and Sheu, 2019), but the drones can operate multiple deliveries per trip and the drone-truck assignment is open. The authors propose a branch-and-price algorithm to solve up to 13 customers. They also provide a sensitivity analysis over the maximal authorized flying duration of the drones. Poikonen et al. (2017) consider the battery of the drones, a different metric for the drones and the trucks,
delivery problem, where a truck and a drone have to serve customers to maxi-
mize a profit function under truck route uncertainty. They consider truck travel
duration uncertainty is considered due to network congestion. They propose a
branch and price algorithm to solve instances with up to 40 customers.

Researchers have studied the cooperation of several drones with the truck to im-
prove efficiency further. In (Karak and Abdelghany, 2019), the authors propose
the "Hybrid Vehicle-Drone Routing Problem" (HVDRP), where a single vehicle
carries a fleet of drones performing pickup and delivery operations. They pro-
pose a MIP formulation along with three heuristics to solve the problem. Murray
and Raj (2020) introduced the "Multiple Flying Sidekicks Traveling Salesman
Problem" (mFSTSP), which is an extension of the FSTSP with several vehicles.
It considers an arbitrary number of heterogeneous drones that can be deployed
from the depot or the delivery truck. They provide a MIP formulation solving
instances with at most 8 customers and a heuristic approach. Kang and Lee
(2021) introduce the Heterogeneous Drone-Truck Routing Problem (HDTRP),
where drones with different speeds and batteries are set up on a truck and can
fly only when the truck is parked at a customer. They propose a MILP formula-
tion and a branch-and-cut algorithm to solve instances with up to 50 customers
optimally.

## 2.2. Drones integrated in VRP

The use of drones with a set of trucks was first introduced in (Wang et al., 2017), which presents the Vehicle Routing Problem with Drones (VRPD), where multiple trucks are equipped with drones to serve a set of customers. In this problem, the drones carried by a truck can deliver one parcel per flight, and the truck can move between the launch and the retrieval of its drones. The authors analyze different worst-case scenarios and provide upper bounds on the savings induced by using drones. The work is extended in (Wang and Sheu, 2019), but the drones can operate multiple deliveries per trip, and the drone-truck assignment is open. The authors propose a branch-and-price algorithm to solve instances with up to 13 customers. They also provide a sensitivity analysis of the drones' maximal authorized flying duration. Poikonen et al. (2017) consider the battery of the drones, a different metric for the drones and the trucks, as well as a cost objective function. They provide theoretical results regarding the VRPD.

Di Puglia Pugliese and Guerriero (2017) present the first MIP model for an extension of the VRP-D called the vehicle-drone routing problem with time windows (VDRPTW), which also considers multiple drones per truck, time windows, and truck capacity. They also provide a numerical study comparing VRP to VRP-D. They concluded that drones are only advantageous when the drone's transportation cost is lower than the truck's. This work is extended in (Di Puglia Pugliese et al., 2020) where the authors present a MIP model solving up to 15 customer instances. They also compare VRP, VRP-D, and routing problem with drones (RPD), where the drones serve customers directly from the depot, on their CO2 emission and traffic congestion. Sacramento et al. (2019) address an extension of the VRP-D where the drones can return to the depot
as well as a cost objective function. They provide theoretical results regarding the VRPD. Di Puglia Pugliese and Guerriero (2017) present the first MIP model for an extension of the VRP-D called the vehicle-drone routing problem with time windows (VDRPTW), which also considers multiple drones per truck, time windows and truck capacity. They also provide a numerical study comparing VRP to VRP-D. They analyse that drone is only useful when the drone's transportation cost is lower than the truck's transportation cost. This work is extended in (Di Puglia Pugliese et al., 2020) where the authors present a MIP model solving up to 15 customers instances. They also compare VRP, VRP-D and routing problem with drones (RPD), where the drones deliver customers directly from the depot, on their CO2 emission and traffic congestion. Sacramento et al. (2019) address an extension of the VRP-D where the drones can return to the depot directly. They present an adaptive large neighborhood search and with a MIP model. In (Schermer et al., 2019), the authors consider the Vehicle Routing Problem and En Route Operation (VRPDERO), where a drone can be launched and retrieved on a moving truck. They propose a MIP model and a heuristic. In (Kitjacharoenchai et al., 2020) the objective function is to minimize the total completion time. A truck can only launch or retrieve a single drone upon reaching a node. A drone can deliver multiple customers before returning to its assigned truck. The authors propose a mixed integer non-linear program (MINLP) and solve instances up to 9 customers optimally. Tamke and Buscher (2021) proposed a branch-and-cut algorithm for the VRP-D without truck capacity solving up to 30 customers instances. The work is extended in (Tamke and Buscher, 2023), where the speed at which a drone performs a flight must be selected from a discrete set, and influence the energy consumption. The authors propose a MIP model along with valid inequalities and preprocessing methods to accelerate the solving. They solve instances with 20 customers optimally, and find feasible solutions on 50 customers instances with a heuristic version of their MIP.

Although the VRPD and its extensions are the most studied, other works related to the integration of drones in vehicle routing problems are also considered. Ulmer and Thomas (2018) consider a dynamic variant, the Same-Day Delivery Routing Problem with Heterogeneous Fleets of drones and vehicles (SDDPHF), where customers requests arrive while the vehicles are already on the road. Drones make deliveries from the depot, serving a single customer before coming back to the depot. The authors present a policy function approximation, to decide whether a request is fulfilled by a drone or a vehicle. Kitjacharoenchai et al. (2019) introduce the multiple travelling salesman problem with drones (mTSPD), where multiple vehicles and drones deliver a set of customers from a depot. They do not limit the number of drones on a vehicle and the truck-drone assignment is open (i.e. the launch and retrieve operations can be operated by different trucks). The authors propose a MIP model solving instances with up to 9 customers instances and a heuristic. Dukkanci et al. (2021) proposes the Energy Minimizing and Range Constrained Drone Delivery Problem (ERDDP) where trucks move from the depot to parking spots. Drones can fly from a parking spot or from the depot to deliver a single customer. The objective function
directly. They present an Adaptive Large Neighborhood Search (ALNS) and a MIP model. In (Schermer et al., 2019), the authors consider the Vehicle Routing Problem and En Route Operation (VRPDERO), where a drone can be launched and retrieved on a moving truck. They propose a MIP model and a heuristic. In (Kitjacharoenchai et al., 2020), the objective function is to minimize the total completion time. A truck can only launch or retrieve a single drone upon reaching a node. A drone can deliver parcels to multiple customers before returning to its assigned truck. The authors propose a mixed integer non-linear program (MINLP) that optimally solves instances with up to 9 customers. Tamke and Buscher (2021) proposed a branch-and-cut algorithm for the VRP-D without truck capacity solving up to 30 customer instances. The work is extended in (Tamke and Buscher, 2023), where the speed at which a drone performs a flight must be selected from a discrete set and influence the energy consumption. The authors propose a MIP model, valid inequalities, and preprocessing methods to accelerate the solving. They solve instances with 20 customers optimally and find feasible solutions on 50 customer instances with a heuristic version of their MIP.

Although the VRPD and its extensions are the most studied, other works related to the integration of drones in vehicle routing problems are also considered. Ulmer and Thomas (2018) consider a dynamic variant, the Same-Day Delivery Routing Problem with Heterogeneous Fleets of drones and vehicles (SDDPHF), where customers requests arrive while the vehicles are already on the road. Drones make deliveries from the depot, serving a single customer before coming back to the depot. The authors present a policy function approximation to decide whether a request is fulfilled by a drone or a vehicle. Kitjacharoenchai et al. (2019) introduce the multiple traveling salesman problem with drones (mTSPD), where multiple vehicles and drones serve a set of customers from a depot. They do not limit the number of drones on a vehicle, and the truck-drone assignment is open (i.e., the launch and retrieve operations can be operated by different trucks). The authors propose a MIP model solving instances with up to 9 customers and a heuristic. Dukkanci et al. (2021) propose the Energy Minimizing and Range Constrained Drone Delivery Problem (ERDDP), where trucks move from the depot to parking spots. Drones can fly from a parking spot or from the depot to deliver a single customer. The objective function minimizes a cost function over the energy spent and the operational costs. The drone's speed directly impacts energy consumption and is considered a decision variable. The authors conduct a parametric analysis of the number of vehicles used. Meng et al. (2023) introduce the multi-visit drone routing problem for pickup and delivery services (MDRP-PD). They propose a MIP and a two-stage heuristic based on simulated annealing. Zhou et al. (2023) present the two-echelon vehicle routing problem with drones (2E-VRP-D). We recall that in this problem, the number of drones on a truck is a decision variable, which allows a better use of the drone fleet. The drones can only fly when their assigned truck is parked at a nearby client. To solve the problem, they presented a MIP, a branch and price, and a heuristic. They managed to solve up to 30 clients to optimality. A similar problem is studied in (Yin et al., 2024), which considers routing
minimizes a cost function over the energy spent, and the operational costs. The drone speed directly impacts its energy consumption, and is considered a decision variable. The authors conduct a parametric analysis over the number of vehicles used. Meng et al. (2023) introduce the multi-visit drone routing problem for pickup and delivery services (MDRP-PD), propose a MIP and a two-stage heuristic based on simulated annealing. Zhou et al. (2023) presented the two-echelon vehicle routing problem with drones (2E-VRP-D). We recall that in this problem, the number of drones on a truck is a decision variable, which is allowing for a better use of the drone fleet. The drones can only fly when their assigned truck is parked at a nearby client. To solve the problem, they presented a MIP, a branch and price and a heuristic. They managed to solve up to 30 clients to optimality. A very similar problem is studied in (Yin et al., 2024), which consider the routing of vehicle-and-drone cooperative to deliver blood products to hospitals. The main difference between the two problems is the objective function, as Yin et al. (2024) penalise the squared weighted arrival time at each hospital. They use a Benders-decomposition column-generation algorithm to solve instances up to 45 customers.

To give an overview of the literature regarding exact approaches for drones integrated in VRP, we present a summary table. Let us first recall notation for the vehicle routing problems with drones introduced in (Tamke and Buscher, 2021). The criteria are the following: number of trucks (#t), number of drones per truck (#d), drone-truck assignment (dta) i.e. whether the drone has an assigned truck or not, launch and retrieve operations (lro), drone range (ran), drone capacity (cap), and objective function (obj). Most criteria are self-explanatory. In the table 1, the keyword 'return' refers to the ability of a truck to return to an already visited customer. The keyword 'loops' refers to the ability of a drone to operate loop operation, in which the start and end location of the flight is the same.
vehicle-and-drone cooperatives to deliver blood products to hospitals. The primary difference between the two problems lies in the objective function, as Yin et al. (2024) penalize the weighted, squared arrival time at each hospital. They use a Benders-decomposition column-generation algorithm to solve instances of up to 45 customers.

To give an overview of the literature regarding exact approaches for drones integrated in VRP, we present a summary table. Let us first recall the notation for the vehicle routing problems with drones introduced in (Tamke and Buscher, 2021). The criteria are the following: number of trucks (#t), number of drones per truck (#d), drone-truck assignment (dta) i.e., whether the drone has an assigned truck or not, launch and retrieve operations (lro), drone range (ran), drone capacity (cap), and objective function (obj). Most criteria are self-explanatory. In Table 1, the keyword 'return' refers to the ability of a truck to return to an already visited customer. The keyword 'loops' refers to the ability of a drone to operate loop operations, in which the start and end locations of the flight are the same.

<table><caption>Table 1: Classification criteria for delivery systems with truck-drone cooperation. (Tamke and Buscher, 2021)</caption><thead><tr><th>Criterion</th><th>Value</th><th>Description</th><th>Criterion</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td rowspan="2">#t</td><td>l</td><td>single truck</td><td rowspan="4">ran</td><td>t</td><td>time</td></tr><tr><td>m</td><td>multiple trucks</td><td>d</td><td>distance</td></tr><tr><td rowspan="2">#d</td><td>l</td><td>single drone</td><td>n</td><td>number of nodes</td></tr><tr><td>m</td><td>multiple drones</td><td>e</td><td>energy consumption</td></tr><tr><td>u</td><td>unlimited number</td><td>u</td><td>unrestricted</td></tr><tr><td rowspan="2">dta</td><td>f</td><td>fixed</td><td rowspan="2">cap</td><td>m</td><td>multiple modes</td></tr><tr><td>o</td><td>open</td><td>l</td><td>single customer</td></tr><tr><td rowspan="3">lro</td><td>a</td><td>return and loops allowed</td><td rowspan="3"></td><td>m</td><td>multiple customers</td></tr><tr><td>p</td><td>return and loops prohibited</td><td></td><td></td></tr><tr><td>l</td><td>only loops allowed</td><td></td><td></td></tr><tr><td rowspan="6">obj</td><td>min ct</td><td>minimize completion time</td><td></td><td></td><td></td></tr><tr><td>min tc</td><td>minimize total costs</td><td></td><td></td><td></td></tr><tr><td>min wt</td><td>minimize customer waiting times</td><td></td><td></td><td></td></tr><tr><td>bi-obj</td><td>minimize total costs and completion time</td><td></td><td></td><td></td></tr><tr><td>min max ct</td><td>minimize maximum completion time</td><td></td><td></td><td></td></tr><tr><td>min total ct</td><td>minimize total completion time</td><td></td><td></td><td></td></tr></tbody></table>

Using the notation defined in Table 1, we present Table 2, providing information about the problem studied, the exact resolutions method and the results of different paper considering the integration of drones to vehicle routing problems. The algorithms are abbreviated as follows: BP stands for branch-and-price, BC for branch-and-cut and BCP for branch-and-cut-and-price.
<table><caption>Table 1: Classification criteria for delivery systems with truck-drone cooperation. (Tamke and Buscher, 2021)</caption><thead><tr><th>Criterion</th><th>Value</th><th>Description</th><th>Criterion</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td rowspan="2">#t</td><td>l</td><td>single truck</td><td rowspan="4">ran</td><td>t</td><td>time</td></tr><tr><td>m</td><td>multiple trucks</td><td>d</td><td>distance</td></tr><tr><td rowspan="2">#d</td><td>l</td><td>single drone</td><td>n</td><td>number of nodes</td></tr><tr><td>m</td><td>multiple drones</td><td>e</td><td>energy consumption</td></tr><tr><td></td><td>u</td><td>unlimited number</td><td>u</td><td>unrestricted</td></tr><tr><td rowspan="2">dta</td><td>f</td><td>fixed</td><td rowspan="4">cap</td><td>m</td><td>multiple modes</td></tr><tr><td>o</td><td>open</td><td>l</td><td>single customer</td></tr><tr><td rowspan="2">lro</td><td>a</td><td>return and loops allowed</td><td>m</td><td>multiple customers</td></tr><tr><td>p</td><td>return and loops prohibited</td><td></td><td></td></tr><tr><td></td><td>l</td><td>only loops allowed</td><td></td><td></td></tr><tr><td rowspan="5">obj</td><td>min ct</td><td>minimize completion time</td><td></td><td></td><td></td></tr><tr><td>min tc</td><td>minimize total costs</td><td></td><td></td><td></td></tr><tr><td>min wt</td><td>minimize customer waiting times</td><td></td><td></td><td></td></tr><tr><td>bi-obj</td><td>minimize total costs and completion time</td><td></td><td></td><td></td></tr><tr><td>min max ct</td><td>minimize maximum completion time</td><td></td><td></td><td></td></tr><tr><td></td><td>min total ct</td><td>minimize total completion time</td><td></td><td></td><td></td></tr></tbody></table>

Using the notation defined in Table 1, we present Table 2, providing information about the problem studied, the exact resolutions method and the results of different paper considering the integration of drones to vehicle routing problems. The Approach are abbreviated as follows: BP stands for branch-and-price, BC for branch-and-cut and BCP for branch-and-cut-and-price.
<table><caption>Table 2: Overview of exact approaches for drones integrated in VRP</caption><thead><tr><th>Reference</th><th>Problem</th><th>#t</th><th>#d</th><th>dta</th><th>lro</th><th>ran</th><th>cap</th><th>obj</th><th>Approach</th><th>Size</th><th>Note</th></tr></thead><tbody><tr><td>Sacramento et al. (2019)</td><td>VRP-D</td><td>⟨m</td><td>1</td><td>f</td><td>p</td><td>t</td><td>1</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>12</td><td></td></tr><tr><td>Sacramento et al. (2019)</td><td>VRP-D</td><td>⟨m</td><td>1</td><td>f</td><td>p</td><td>t</td><td>1</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>12</td><td></td></tr><tr><td>Di Puglia Pugliese and Guerriero (2017)</td><td>VDRPTW</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>d</td><td>1</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>10</td><td>Time Window</td></tr><tr><td>Meng et al. (2023)</td><td>MDRP-PD</td><td>⟨m</td><td>1</td><td>f</td><td>p</td><td>e</td><td>m</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>11</td><td>Pickup and Delivery</td></tr><tr><td>Wang and Sheu (2019)</td><td>VRPD</td><td>⟨m</td><td>m</td><td>o</td><td>p</td><td>t</td><td>m</td><td>min tc</td><td>⟩</td><td>BP</td><td>13</td><td>Docking Hubs</td></tr><tr><td>Schermer et al. (2019)</td><td>VRPDERO</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>d</td><td>1</td><td>min max ct</td><td>⟩</td><td>MIIP</td><td>10</td><td>En-route Operation</td></tr><tr><td>Kitjacharoenchai et al. (2019)</td><td>mTSPD</td><td>⟨m</td><td>u</td><td>o</td><td>p</td><td>u</td><td>1</td><td>min max ct</td><td>⟩</td><td>MIIP</td><td>9</td><td></td></tr><tr><td>Kitjacharoenchai et al. (2020)</td><td>VRPD</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>d</td><td>m</td><td>min total ct</td><td>⟩</td><td>MINLP</td><td>9</td><td></td></tr><tr><td>Di Puglia Pugliese et al. (2020)</td><td>VRP-D</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>d</td><td>1</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>15</td><td>Time Window</td></tr><tr><td>Dukkanci et al. (2021)</td><td>ERDDP</td><td>⟨m</td><td>m</td><td>f</td><td>1</td><td>e</td><td>1</td><td>min tc</td><td>⟩</td><td>SOCP</td><td>37</td><td>Drone Speed Variable</td></tr><tr><td>Tamke and Buscher (2021)</td><td>VRPD</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>d</td><td>1</td><td>min max ct</td><td>⟩</td><td>BC</td><td>20</td><td></td></tr><tr><td>Zhou et al. (2023)</td><td>2E-VRP-D</td><td>⟨m</td><td>m</td><td>f</td><td>1</td><td>e</td><td>1</td><td>min total ct</td><td>⟩</td><td>BP</td><td>30</td><td>Variable nb Drone per Truck</td></tr><tr><td>Tamke and Buscher (2023)</td><td>VRPD-DSS</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>e</td><td>1</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>20</td><td>Drone Speed Variable</td></tr><tr><td>Yin et al. (2024)</td><td>BPDRP</td><td>⟨m</td><td>m</td><td>f</td><td>1</td><td>t</td><td>1</td><td>min tc</td><td>⟩</td><td>BP</td><td>45</td><td>Quadratic objective</td></tr><tr><td>Our</td><td>2E-VRP-D</td><td>⟨m</td><td>m</td><td>f</td><td>1</td><td>e</td><td>1</td><td>min total ct</td><td>⟩</td><td>BCP</td><td>120</td><td>Variable Drone per Truck</td></tr></tbody></table>
<table><caption>Table 2: Overview of exact approaches for drones integrated in VRP</caption><thead><tr><th>Reference</th><th>Problem</th><th>#t</th><th>#d</th><th>dta</th><th>ho</th><th>ran</th><th>cap</th><th>obj</th><th>Approach</th><th>Size</th><th>Note</th></tr></thead><tbody><tr><td>Sacramento et al. (2019)</td><td>VRP-D</td><td>⟨m</td><td>1</td><td>f</td><td>p</td><td>t</td><td>1</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>12</td><td></td></tr><tr><td>Sacramento et al. (2019)</td><td>VRP-D</td><td>⟨m</td><td>1</td><td>f</td><td>p</td><td>t</td><td>1</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>12</td><td></td></tr><tr><td>Di Puglia Pugliese and Guerriero (2017)</td><td>VDRPTW</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>d</td><td>1</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>10</td><td>Time Window</td></tr><tr><td>Meng et al. (2023)</td><td>MDRP-PD</td><td>⟨m</td><td>1</td><td>f</td><td>p</td><td>e</td><td>m</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>11</td><td>Pickup and Delivery</td></tr><tr><td>Wang and Sheu (2019)</td><td>VRPD</td><td>⟨m</td><td>m</td><td>o</td><td>p</td><td>t</td><td>m</td><td>min tc</td><td>⟩</td><td>BP</td><td>13</td><td>Docking Hubs</td></tr><tr><td>Schermer et al. (2019)</td><td>VRPDERO</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>d</td><td>1</td><td>min max ct</td><td>⟩</td><td>MIIP</td><td>10</td><td>En-route Operation</td></tr><tr><td>Kitjacharoenchai et al. (2019)</td><td>mTSPD</td><td>⟨m</td><td>u</td><td>o</td><td>p</td><td>u</td><td>1</td><td>min max ct</td><td>⟩</td><td>MIIP</td><td>9</td><td></td></tr><tr><td>Kitjacharoenchai et al. (2020)</td><td>VRPD</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>d</td><td>m</td><td>min total ct</td><td>⟩</td><td>MINLP</td><td>9</td><td></td></tr><tr><td>Di Puglia Pugliese et al. (2020)</td><td>VRP-D</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>d</td><td>1</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>15</td><td>Time Window</td></tr><tr><td>Dukkanci et al. (2021)</td><td>ERDDP</td><td>⟨m</td><td>m</td><td>f</td><td>1</td><td>e</td><td>1</td><td>min tc</td><td>⟩</td><td>SOCP</td><td>37</td><td>Drone Speed Variable</td></tr><tr><td>Tamke and Buscher (2021)</td><td>VRPD</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>d</td><td>1</td><td>min max ct</td><td>⟩</td><td>BC</td><td>20</td><td></td></tr><tr><td>Zhou et al. (2023)</td><td>2E-VRP-D</td><td>⟨m</td><td>m</td><td>f</td><td>1</td><td>e</td><td>1</td><td>min total ct</td><td>⟩</td><td>BP</td><td>30</td><td>Variable nb Drone per Truck</td></tr><tr><td>Tamke and Buscher (2023)</td><td>VRPD-DSS</td><td>⟨m</td><td>m</td><td>f</td><td>p</td><td>e</td><td>1</td><td>min tc</td><td>⟩</td><td>MIIP</td><td>20</td><td>Drone Speed Variable</td></tr><tr><td>Yin et al. (2024)</td><td>BPDRP</td><td>⟨m</td><td>m</td><td>f</td><td>1</td><td>t</td><td>1</td><td>min tc</td><td>⟩</td><td>BP</td><td>45</td><td>Quadratic objective</td></tr><tr><td>Our</td><td>2E-VRP-D</td><td>⟨m</td><td>m</td><td>f</td><td>1</td><td>e</td><td>1</td><td>min total ct</td><td>⟩</td><td>BCP</td><td>120</td><td>Variable Drone per Truck</td></tr></tbody></table>
### 3. Formulation

We now introduce the addressed problem. The 2E-VRP-D was first defined in (Zhou et al., 2023). In this problem, $k_v$ vehicles and $k_d$ drones have to deliver a set of customers $V = \{1, ..., n\}$ from a depot denoted 0. Each truck can carry at most $D$ drones.

The following assumptions are made:

(a) Both drones and vehicles deliver packages, but some customers have to be served by truck.

(b) A drone is limited by its capacity and battery but is fully recharged upon reaching the vehicle.

(c) A drone can only carry one package per flight and is synchronized with a single assigned vehicle.

(d) A drone can only fly when its vehicle is parked at a customer but can take off multiple times from the same customer's location.

(e) Each drone and its necessary equipment have a weight $q_1^d$.

(f) The total capacity of a vehicle, denoted by $Q$, limits the weight of packages and drones it can carry.

(g) The service time at a customer delivered by a vehicle and the drone deliveries from this customer occur simultaneously.

We define $N_i$, $i \in V$, the drone-reachable neighborhood of $i$, as the set of customers that can be served by drone from $i$. For completeness, a formal definition of $N_i$, $i \in V$ can be found in Appendix A. We define a drone schedule as an ordered set of round trips for each drone originating from a truck parked at a customer location. The problem is defined by the complete symmetric directed graph $G = (V_+, A)$, where $V_+ = \{0\} \cup V$ is the set of vertices and $A = \{(i, j) | i, j \in V_+, i \neq j\}$ is the set of arcs. A customer $i \in V$ is defined by a demand $q_i$, a vehicle service time $s_i^v$, a drone service time $s_i^d$, an autonomous drone setup time $o_i$, and a deadline $\bar{d}_i$, which is the latest date a vehicle can arrive to serve the customer $i$. Travelling time between two nodes $i, j \in V_+$ is noted $c_{i,j}^v$ for a vehicle and $c_{i,j}^d$ for a drone.

We now present our mathematical formulation of the problem. We define a route as the combination of a truck route for the first echelon and a set of drone schedules for the second echelon. Let $\mathcal{R}^d$, $d \in \{0, ..., D\}$, be the set of all feasible routes of a vehicle carrying $d$ drones. For every $d \in \{0, ..., D\}$, every $i \in V$, and every subset $K \subseteq N_i$, let $S_{iK}^d$ be the set of all optimal drone schedules serving the set of customers $K$ with $d$ drones from the customer $i$. Set $S_{iK}^d$ may contain several schedules because the arrival date of the truck at customer $i$ may vary, and the customers have deadlines. Hence, the optimal solution may not be feasible for a given truck arrival date. Note that set $S_{iK}^d$ contains the Pareto front, minimizing the maximal drone schedule duration and maximizing the maximal truck arrival date. We note $m_s$, $s \in S_{iK}^d$ the makespan of the drone schedule $s$. Let $b_{i,j}^r$ be equal to one if and only if the vehicle moves from $i$ to $j$ in route $r \in \mathcal{R}^d$. Let $a_s^r$ be equal to one if and only if the drones follow drone schedule $s \in S_{iK}^d$ when the vehicle is parked at customer $i \in V$ in route $r$. Let variable
## 3 Formulation

We now introduce the addressed problem. The 2E-VRP-D was first defined in (Zhou et al., 2023). In this problem a set of vehicles denoted $K_v = \{1, ..., k_v\}$ and a set of drones $K_d = \{1, ..., k_d\}$ have to serve a set of customers $V = \{1, ..., n\}$ from a depot 0. Let $D = \{0, ..., \Gamma\}$ be the set of all possible numbers of drones a vehicle can carry. The following assumptions are made:

(a) Both drones and vehicles deliver packages, but some customers can not be delivered by drones.

(b) A drone is limited by its capacity and its battery but is fully recharged upon reaching the vehicle.

(c) A drone can only carry one package per flight, and is synchronised with a single assigned vehicle.

(d) A drone can fly only when its vehicle is parked at a customer, but can take off multiple times from the same customer.

(e) The vehicle total capacity is noted $Q$. The vehicle payload is impacted by the amount of carried drones, each drone and its necessary equipment having a weight $q_1^d$.

(f) The vehicle service time at a customer and the drone deliveries from this customer happen in parallel.

We define $K_i, i \in V$, the drone-reachable neighborhood of $i$, as the set of customers that can be delivered by drone from $i$. We call a drone schedule an ordered set of round trips of the drones from a customer. The problem is defined by the complete symmetric directed graph $G = (V_+, A)$, where $V_+ = \{0\} \cup V$ is the set of vertices and $A = \{(i,j) | i,j \in V_+, i \neq j\}$ is the set of arcs. A customer $i \in V$ has a demand $q_i$, a vehicle service time $s_i^v$, a drone service time $s_i^d$ and a deadline $\bar{d}_i$, which is the latest date a vehicle can arrive to serve the customer $i$. Travelling time between two nodes $i, j \in V_+$ is noted $c_{i,j}^v$ for a vehicle and $c_{i,j}^d$ for a drone.

We now present our mathematical formulation of the problem. Let $\mathcal{R}^d, d \in D$, be the set of all feasible routes of a vehicle carrying $d$ drones. For every $d \in D$, every $i \in V$, and every subset $K \subseteq K_i$, let $S_{iK}^d$ be the set of all optimal drone schedule serving the set of customers $K$ with $d$ drones from the customer $i$. Set $S_{iK}^d$ may contain several schedules, because the arrival date of the truck at customer $i$ may vary and the customers have deadlines, hence the optimal solution may not be feasible for a given truck arrival date. Note that set $S_{iK}^d$ contains the Pareto front minimizing the maximal drone schedule duration, and maximizing the maximal truck arrival date. Let $b_{i,j}^r$ be equal to one if and only if the vehicle moves from $i$ to $j$ in route $r \in \mathcal{R}^d$. Let $a_s^r$ be equal to one if and only if the drones follow drone schedule $s \in S_{iK}^d$ when the vehicle is parked at customer $i \in V$ in route $r$. Let variable $\lambda_r, r \in \mathcal{R}^d, d \in D$ be equal to 1 if the route $r$ is used in the solution and 0 otherwise. We also define the following additional variables :

$x_{ij}^d$ — an integer variable which is equal to the number of times a vehicle
$\lambda_r$, $r \in \mathbb{R}^d$, $d \in \{0, ..., D\}$ be equal to 1 if the route $r$ is used in the solution and 0 otherwise. We also define the following additional variables :

$x_{ij}^d$ — an integer variable equal to the number of times a vehicle with $d \in \{0, ..., D\}$ drones travels between node $i \in V_+$ and node $j \in V_+$, $i < j$, regardless of the order of visit.

$y^d$ — a continuous variable equal to the total truck immobility time in routes of vehicles with $d$ drones.

$z_j^d$ — a binary variable equal to one if and only if a vehicle with $d$ drones sends a drone to visit customer $j \in V$.

$u^d$ — the integer variable equal to the number of used vehicles with $d$ drones.

Variables with superscript $d$ belong to subproblem $d \in \{0, ..., D\}$.

The set partitioning formulation is presented in two parts, using the previous notations. The first part (i.e., equations (2) to (4)) uses the additional variables to model the problem. The second part (i.e., equations (5) to (8)) maps the variable $\lambda_r$ to the additional variables.

$$ \min \sum_{d \in \{0, ..., D\}} \left( \sum_{i \in V_+} \sum_{j \in V_+} c_{i,j}^v \times x_{ij}^d + y^d \right) \qquad (1) $$

$$ \text{s.t.} \quad \sum_{d \in \{0, ..., D\}} \left( \sum_{j=1}^{i-1} \frac{1}{2} x_{ji}^d + \sum_{j=i+1}^{n} \frac{1}{2} x_{ij}^d + z_i^d \right) = 1, \quad \forall i \in V, \qquad (2) $$

$$ \sum_{d \in \{0, ..., D\}} u^d \leq k_v, \qquad (3) $$

$$ \sum_{d \in \{0, ..., D\}} d \times u^d \leq k_d, \qquad (4) $$

$$ x_{ij}^d = \sum_{r \in \mathcal{R}^d} (b_{ij}^r + b_{ji}^r) \times \lambda_r, \quad \forall d \in \{0, ..., D\}, i, j \in V_+, i < j, \qquad (5) $$

$$ y^d = \sum_{r \in \mathcal{R}^d} \sum_{j \in V} \sum_{K \in \mathcal{K}_j} \sum_{s \in S_{jK}^d} (m_s \times a_s^r) \times \lambda_r, \quad \forall d \in \{0, ..., D\}, \forall i \in V, \qquad (6) $$

$$ z_i^d = \sum_{r \in \mathcal{R}^d} \sum_{j \in V} \sum_{K \in \mathcal{K}_j} \sum_{i \in K} a_s^r \times \lambda_r, \quad \forall d \in \{0, ..., D\}, \forall i \in V \qquad (7) $$

$$ u^d = \sum_{r \in \mathcal{R}^d} \sum_{j \in V} b_{0j}^r \times \lambda_r, \quad \forall d \in \{0, ..., D\}, \qquad (8) $$

$$ \lambda_r \in \{0, 1\}, \quad \forall d \in \{0, ..., D\}, \forall r \in \mathcal{R}^d. \qquad (9) $$
with $d \in D$ drones travels between node $i \in V_+$ and node $j \in V_+$, $i < j$, regardless of the order of visit.

$y^d$ — a continuous variable which is equal to the total truck immobility time in routes of vehicles with $d$ drones.

$z_j^d$ — a binary variable which is equal to one if and only if a vehicle with $d$ drones sends a drone to visit customer $j \in V$.

$u^d$ — the integer variable which is equal to the number of used vehicles with $d$ drones.

Variables with superscript $d$ belong to subproblem $d \in D$.

Using the previous notations, a set partitioning formulation is presented in two parts. The first part (i.e. equations (2) to (4)) uses the additional variables to model the problem. The second part (i.e. equations (5) to (8)) maps the variable $\lambda_r$ to the additional variables.

$$ \min \sum_{d \in D} \left( \sum_{i \in V_+} \sum_{j \in V_+} c_{i,j}^v \times x_{ij}^d + y^d \right) \qquad (1) $$

$$ \text{s.t.} \quad \sum_{d \in D} \left( \sum_{j=1}^{i-1} \frac{1}{2} x_{ji}^d + \sum_{j=i+1}^{n} \frac{1}{2} x_{ij}^d + z_i^d \right) = 1, \quad \forall i \in V, \qquad (2) $$

$$ \sum_{d \in D} u^d \leq k_v, \qquad (3) $$

$$ \sum_{d \in D} du^d \leq k_d, \qquad (4) $$

$$ x_{ij}^d = \sum_{r \in \mathcal{R}^d} (b_{ij}^r + b_{ji}^r) \times \lambda_r, \quad \forall d \in D, i, j \in V_+, i < j, \qquad (5) $$

$$ y^d = \sum_{r \in \mathcal{R}^d} \sum_{j \in V} \sum_{K \in \mathcal{K}_j} \sum_{s \in S_{jK}^d} (m_s \times a_s^r) \times \lambda_r, \quad \forall d \in D, \forall i \in V, \qquad (6) $$

$$ z_i^d = \sum_{r \in \mathcal{R}^d} \sum_{j \in V} \sum_{K \in \mathcal{K}_j: s \in S_{jK}^d} a_s^r \times \lambda_r, \quad \forall d \in D, \forall i \in V \qquad (7) $$

$$ u^d = \sum_{r \in \mathcal{R}^d} \sum_{j \in V} b_{0j}^r \times \lambda_r, \quad \forall d \in D, \qquad (8) $$

$$ \lambda_r \in \{0, 1\}, \quad \forall d \in D, \forall r \in \mathcal{R}^d. \qquad (9) $$

The objective function (1) and constraints (2)–(4) belong to the master. The objective function (1) minimizes the total duration of the routes. The constraints (2) ensure that each customer is visited exactly once. The constraints (3) and (4) verify that the solution uses a feasible number of trucks and drones. The constraints (5)–(8) define the mapping between the additional variables and the route variables. Note that additional variables are used mainly for the clarity
The objective function (1) and constraints (2)-(4) belong to the master. The objective function (1) minimizes the total duration of the routes. The constraints (2) ensure that each customer is visited exactly once. The constraints (3) and (4) verify that the solution uses a feasible number of trucks and drones. The constraints (5)-(8) define the mapping between the additional variables and the route variables. Note that additional variables are used mainly for the clarity of the formulation and can be eliminated from the formulation using relations (5)-(8).

# 4. Branch-cut-and-price algorithm

Because it would be impractical to enumerate the set $\mathcal{R}^d, d \in \{0, ..., D\}$ due to its exponential size over the number of customers, we solve this formulation using a column and cut generation approach. We use a similar branch-and-cut-and-price (BCP) as in (Pessoa et al., 2020), with the following features. The BCP algorithm iteratively generates new variables (columns) to improve the linear programming (LP) relaxation of an integer program. Based on the LP solution, branching decisions are made, and cutting planes are added to strengthen the relaxation. Two primal heuristics are used in branch and bounds nodes. The first heuristic solves the restricted master problem using a MIP solver and a time limit. The second heuristic is called the diving heuristic, where columns with the most significant fractional values are iteratively set to one, and the remaining problem is solved using column generation.

## 4.1. Column generation

We use three-stage column generation. The pricing is solved using heuristic labeling algorithms in the two first stages as in (Sadykov et al., 2021), and using the exact labeling dynamic algorithm as in (Pessoa et al., 2020) and (Sadykov et al., 2021). The first two stages generate at most 30 columns, and the last stage at most 150. Automatic dual price smoothing stabilization is used, as in (Pessoa et al., 2018), to help the convergence of the process. The bucket-arc elimination procedure is also used, removing arcs that have proven not to belong to any optimal solution. After each call, we use a route enumeration technique similar to (Baldacci et al., 2008), which tries to enumerate all improving routes for a subproblem. If it succeeds, the future pricing problems will be solved by inspection of the enumerated routes. Furthermore, if the number of enumerated routes in all subproblems is less than a certain threshold (we use 5000 as the threshold), all the routes are added to the restricted master problem and solved by the MIP solver.

## 4.2. Cutting planes

### 4.2.1. Adapted rounded capacity cuts

The rounded capacity cuts were first introduced in (Laporte and Nobert, 1983) for the constrained vehicle routing problem (CVRP). Let $C \subset V$ and let $h_r$ be
of the formulation and can be eliminated from the formulation using relations
(5)–(8).

4 Pricing subproblem

We design the pricing subproblems as Resource Constrained Shortest Path Prob-
lem (RCSPP), where a feasible solution is a feasible route for the 2E-VRP-D.
The objective of this problem is to find the shortest path between a source and
a sink in a graph where each arc consumes resources and requires that the re-
maining resources lie between lower and upper limits.

We solve the subproblems using the bucket graph labelling algorithm presented
in (Sadykov et al., 2021). The algorithm tries to extend partial paths called
labels to form complete paths. Given their set of accumulated resources, labels
are stored and extended according to so-called buckets, helping greatly with
dominance checks. This algorithm supports ng-paths relaxation as well as the
presence of limited rank-1 cuts.

We model a RCSP problem for each $d \in D$, in a graph defined as follows. Let $\mathcal{G}^d = (\mathcal{V}^d, \mathcal{A}^d)$ be a directed graph, where $\mathcal{V}^d = V_+ \cup V_+'$, with $V_+'$ a copy of $V_+$. Set $\mathcal{A}^d = A^{dT} \cup A^{dD}$ with $A^{dT}$ the set of arcs representing a truck delivery possible with $d$ drones and $A^{dD}$ the set of arcs representing possible drone deliveries with $d$ drones. The source of the graph is node $0'$ and the sink is node 0. Arcs in $A^{dD}$ go from a vertex $i \in V$ to the vertex $i' \in V'$, which is the copy of $i$. Let $S_{iK}^d$ be the set of drone routes with $d$ drone, delivering the set of customers $K \subseteq \mathcal{K}_i$ from $i$. For each $i \in V$, each possible subset $K \subseteq \mathcal{K}_i$, and each drone routes $s \in S_{iK}^d$, we define an arc in set $A^{dD}$, labelled $(i, i')^s$. Let $a_s^r$ be equal to one if and only if route $r \in \mathcal{R}^d$ follows this arc, i.e., the drones follow drone route $s \in S_{iK}^d$ in route $r$ when the vehicle is positioned at customer $i \in V$. Arcs in $A^{dT}$ go from a vertex $i' \in V_+'$ to a vertex $j \in V^+$. Let $b_{i,j}^r$ be equal to one if and only if route $r \in \mathcal{R}^d$ follows arc $(i', j)$. Note that the prime is omitted in the notation without ambiguity, as there is no arc between $i$ and $j$ with $i, j \in V_+$.
the number of times $r \in \mathbb{R}^d$ visits C i.e.

$$h_r = \frac{1}{2} \times \sum_{i \in V_+ \setminus C} \sum_{j \in C} b_{\min(i,j), \max(i,j)}^r + \sum_{i \in V_+ \setminus C} \sum_{K \in \mathcal{N}_i: K \cap C \neq \emptyset} \sum_{s \in S_{iK}^d} a_s^r \quad \forall d \in D, r \in \mathbb{R}^d \quad (10)$$

These valid inequalities can be formulated as follows:

$$\sum_{d \in D} \sum_{r \in \mathcal{R}^d} \lambda_r \times h_r \geq \left\lfloor \frac{\sum_{i \in C} q_i}{Q} \right\rfloor, \quad C \subset V \quad (11)$$

However, in this problem, drones occupy a significant portion of the truck's capacity (about 20% per drone in the instances), which must be taken into consideration. Let $k = \left\lfloor \frac{Q}{q_1^d} \right\rfloor$, and we adapt the rounded capacity cuts as follows:

$$\sum_{d \in D} \sum_{r \in \mathcal{R}^d} \lambda_r \times h_r \times (k-d) \geq \left\lfloor \frac{\sum_{i \in C} q_i}{q_1^d} \right\rfloor, \quad C \subset V \quad (12)$$

We will refer to these as adapted rounded capacity cuts (ARCC). We use the greedy construction heuristic from (Lysgaard et al., 2004) to separate the ARCC.

<caption>**Algorithm 1** Separation for the adapted rounded capacity cuts</caption>

**Require:** $t \in \mathbb{R}$
1: **for** $i \in V$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;C $\leftarrow \{i\}$
3: &nbsp;&nbsp;&nbsp;&nbsp;**while** C $\neq$ V **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;C $\leftarrow$ C $\cup$ $\{\arg\max_{j \in V \setminus C} (\text{violation}(C \cup \{j\}))\}$
5: &nbsp;&nbsp;&nbsp;&nbsp;**if** violation(C) > t **then**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;addCut(C)

We describe the separation algorithm in Algorithm (1), where t is a user-defined threshold, the function violation(C) calculates the difference between the right-hand side and the left-hand side in constraints (12) for set C and addCut(C) add the valid inequality to the problem. The algorithm can be explained as follows: starting with a single customer in the set, while not all customers are in the set, we add a customer that is not already in the set, and that maximizes the violation. If the violation of the set exceeds the threshold, we add the corresponding cut to the linear program. We generate at most 100 cuts per round.

### 4.2.2. Limited memory rank-one cuts

We use limited memory rank-one cuts (lm-R1Cs) similar to (Pecin et al., 2017a), (Pessoa et al., 2020). The R1Cs uses a Chvátal-Gomory rounding on a weighted sum of a subset of constraints (2). When the coefficients are all equal, these cuts are called Subset-Row cuts (Jepsen et al., 2008).
<figure><img src="image_3.png" /><figcaption>Figure 1: Example of a subproblem graph $\mathcal{G}^d$</figcaption></figure>

Figure (1) shows an example of a subproblem graph. We represented only three arcs from $i \in V$ to $i' \in V'$ for clarity, but there should be an arc for all $s \in S_{iK}^d, \forall K \subseteq K_i$.

## 4.1 Resources for the pricing subproblem

In our model, we use 2 resources, the duration $\beta_d$ of a route, and the capacity $\beta_c$ of a truck, and we define a pricing subproblem for each $d \in D$. The objective function is the minimization of the reduced cost. The reduced cost of a route $r$ with a cost $c_r$ is given by: $\bar{c}_r = c_r - \sum_{i \in V} \alpha_{ir} \times \pi_i - \lambda_0^v - d_r \times \lambda_0^d$ where $\alpha_{ir}$ is equal to one if and only if customer $i$ is visited in the route $r$, $d_r$ the number of drones used in route $r$, $\pi_i$ are the dual variables associated with constraints (2), $\lambda_0^v$ is the dual variable associated with constraint (3) and $\lambda_0^d$ is dual variable associated with constraint (4).

Let $[l_{a,r}, u_{a,r}]$ be the resource feasibility interval, for each $a \in \mathcal{A}^d$ and $r \in \{\beta_c, \beta_d\}$, i.e. the accumulated resource $r$ has to belong lay in the interval $[l_{a,r}, u_{a,r}]$ to travel through arc $a$. Let $q_{a,r}$ be the consumption of resource $r$ on arc $a$. For the resource capacity we set:

$$l_{a,\beta_c} = 0 \qquad \forall a \in \mathcal{A}^d \qquad (10)$$

$$u_{a,\beta_c} = Q - d \times q_1^d \qquad \forall a \in \mathcal{A}^d \qquad (11)$$

$$q_{a,\beta_c} = \begin{cases} 0 & \forall a = (i, 0), i \in V' \\ q_j & \forall a = (i, j), i \in V'_+, j \in V \\ \sum_{k \in K} q_k & \forall a = (i, i')^s, i \in V, i' \in V', K \in K_i, s \in S_{iK}^d \end{cases} \qquad (12)$$

The truck starts with a remaining capacity equal to $Q - d \times q_1^d$. The remaining capacity of a truck must lie between 0 and $Q - d \times q_1^d$ on every arcs of the graph.
The R1Cs are defined as follows. Let $C \subset V_+$, with multipliers $0 < \rho_i < 1 \quad \forall i \in C$ and let $h_i^r$ be the number of time the vertex $i$ is visited in route $r$:

$$h_i^r = \frac{1}{2} \sum_{j \in C} b_{\min(i,j), \max(i,j)}^r + \sum_{j \in V_+ \setminus C} \sum_{K \in \mathcal{K}_j: s \in S_j^d} a_s^r \quad \forall d \in D, r \in \mathcal{R}^d \quad (13)$$

The following equation is a valid inequality:

$$\sum_{d \in D} \sum_{r \in \mathcal{R}^d} \left[ \rho_i \sum_{i \in C} h_i^r \right] \lambda_r \le \left[ \sum_{i \in C} \rho_i \right] \quad (14)$$

Because the R1Cs are non-robust constraints, we use lm-R1Cs instead, where the set $C$ is paired with a memory set $M$ such that $C \subseteq M \subseteq V_+$. $\left[ \rho_i \sum_{i \in C} h_i^r \right]$ in equation (14) is replaced by $\alpha(C, M, \rho, P)$ and is calculated during the pricing computation as in algorithm 2.

**Algorithm 2** $\alpha(C, M, \rho, P = (v_0 = 0, v_1, ..., v_k, v_k + 1 = 0'))$

1: $\alpha \leftarrow 0, s \leftarrow 0$

2: **for** $i = 1$ **to** $k$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;**if** $(v_{i-1}, v_i) \notin M$ **then**

4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$s \leftarrow 0$

5: &nbsp;&nbsp;&nbsp;&nbsp;**if** $v_i \in C$ **then**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$s \leftarrow s + \rho_{v_i}$

7: &nbsp;&nbsp;&nbsp;&nbsp;**if** $s \ge 1$ **then**

8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$s \leftarrow s - 1, \alpha \leftarrow \alpha + 1$

9: **return** $\alpha$

If the partial path $P$ leaves the set $M$, all previous visits are forgotten, leading to an easier-to-compute but smaller coefficient. Note that if $M = V$, lm-R1Cs are equivalent to R1Cs. To separate those cuts, we find a violated constraint and calculate its minimal memory set such that the limited memory R1Cs have the same violation.

### 4.3. Branching

We select a collection of branching candidates from the total number of drones ($\sum_{d \in \{0, ..., D\}} d \times u^d$), the total number of trucks ($\sum_{d \in \{0, ..., D\}} u^d$), the number of trucks with $d$ drones (variables $u$) and the truck travel variables (variables $x$).

We use a multi-phase strong branching procedure to choose the most promising branching candidate, similar to (Pecin et al., 2017a). The idea is to spend
The consumption of remaining capacity is set to the demand of the customer served, or the sum of the demand of the customers delivered if the arc represents a drone route (i.e. $a = (i, i')i \in V, i' \in V'$), or 0 if the destination of the arc is node 0.

Let $\bar{a}_s$ be the maximal starting time of a solution $s$ that will be formally defined in section (4.2.1). For the resource duration we set:

$$l_{a,\beta_d} = 0 \qquad \forall a \in \mathcal{A}^d \qquad (13)$$

$$u_{a,\beta_d} = \begin{cases} \bar{d}_j & \forall a = (i,j), i \in V'_+, j \in V_+ \\ \bar{a}_s & \forall a = (i, i')^s, i \in V, i' \in V', K \in K_i, s \in S_{iK}^d \end{cases} \qquad (14)$$

$$q_{a,\beta_d} = \begin{cases} c_{i,j}^v & \forall a = (i,j), i \in V'_+, j \in V_+ \\ \max(s_i^v, m_s) & \forall a = (i, i')^s, i \in V, i' \in V', K \in K_i, s \in S_{iK}^d \end{cases} \qquad (15)$$

The truck starts at time 0. A truck arc delivering a customer can only be used if the accumulated duration lie between 0 and the customer's deadline. A drone arc delivering a set of customers can only be used if the accumulated duration lie between 0 and the drone route's maximum starting time. The duration consumption of a drone arc is the maximum between the service time in truck at customer $i$ and the duration of the drone route, which are operated in parallel.

## 4.2 Drone routes enumeration

We enumerate all the feasible drone routes that might belong to an optimal solution. Therefore we must calculate the optimal drone routes from all customers $i \in V$ such that $\bar{f}_i^d = 1$, for any number of drones on the truck $d \in [1, ..., \Gamma]$, and for all sets of available customers, $K \subseteq K_i$. Because of customers deadline, the optimal drone route delivering the set $K$ depends on the truck arrival time at $i$. Let $T_i = [c_{0,i}^v, ..., \bar{d}_i]$ be the set of all possible arrival times of a truck at customer $i$.

### 4.2.1 Scheduling problem

Since the drones can only deliver one customer before going back to the truck, and they are automatically fully recharged when reaching the truck, we can model the drone routing problem for a truck with $d$ drones arrived at a customer $i$ at time $a \in T_i$ and delivering set $K \subseteq K_i$ customers as a scheduling problem with $d$ parallel identical machines, $|K|$ jobs with deadlines, a unique global release date $a$ and minimizing the makespan. The notation of the scheduling problem is $P|\bar{d}_j|C_{max}$ (Graham et al., 1979). Note that the global release date $a$ does not appear in the notation, since the release date is unique for all the jobs, the machines will start working at time $a$. Hence, with a simple trick on the deadlines, we can get rid of the global release date.

In this scheduling problem, if there exists an optimal solution of the problem, there always exists an optimal solution without waiting time. Hence increasing
more time evaluating branching variables in the lowest depth of the branch-and-bound tree where each selection has a greater impact on the overall time and spend less time as the depth increases, taking advantage of the history of previous evaluations.

In the first phase, half the candidates are chosen from history from previous calls of the strong branching procedure. The remaining half is selected in a balanced manner from all branching strategies. The selection method is based on the distance between the variable's value and the nearest integer. The candidates with the largest distances are selected.

In the second phase, selected candidates from the first phase are evaluated by solving the current restricted master linear program modified for each candidate without column generation. We select the candidates using the product rule (Achterberg, 2007).

In the third phase, selected candidates from the second phase are evaluated by solving the restricted master linear program modified for each candidate, with heuristic column generation but without cut generation. We select the best candidate using the product rule.

## 5. Pricing subproblem

We design the pricing subproblems as a Resource Constrained Shortest Path Problem (RCSP), where a feasible solution is a feasible route for the 2E-VRP-D. The objective of this problem is to determine a shortest path between a source and a sink in a graph where each arc consumes resources and requires that the remaining resources lie between lower and upper limits.

We solve the subproblems using the bucket graph labeling algorithm presented in (Sadykov et al., 2021). The algorithm tries to extend partial paths called labels to form complete paths. Given their accumulated resources, labels are stored and extended according to so-called buckets, helping greatly with dominance checks. This algorithm supports ng-paths relaxation (Baldacci et al., 2011) as well as the presence of limited rank-1 cuts (Pecin et al., 2017b).

We model a RCSP problem for each $d \in \{0, ..., D\}$ in a graph defined as follows. Let $G^d = (\mathcal{V}^d, \mathcal{A}^d)$ be a directed graph, where $\mathcal{V}^d = V_+ \cup V'_+$, with $V'_+$ a copy of $V_+$. Set $\mathcal{A}^d = A^{dT} \cup A^{dD}$ with $A^{dT}$ the set of arcs representing a truck delivery possible with $d$ drones and $A^{dD}$ the set of arcs representing possible drone deliveries with $d$ drones. The source of the graph is node 0', and the sink is node 0. Arcs in $A^{dD}$ go from a vertex $i \in V$ to the vertex $i' \in V'$, which is the copy of $i$. Let $S_{iK}^d$ be the set of drone schedules with $d$ drone, serving the set of customers $K \subseteq \mathcal{N}_i$ from $i$. For each $i \in V$, each possible subset $K \subseteq \mathcal{N}_i$, and each drone schedules $s \in S_{iK}^d$, we define an arc in set $A^{dD}$, labelled $(i, i')^s$. Let $a_s^r$ be equal to one if and only if route $r \in \mathcal{R}^d$ follows this arc, i.e., the drones follow drone schedule $s \in S_{iK}^d$ in route $r$ when the vehicle is positioned at customer $i \in V$. Arcs in $A^{dT}$ go from a vertex $i' \in V'_+$ to a vertex $j \in V^+$. Let $b_{i,j}^r$ be equal to one
a can only increase the makespan of a solution. Therefore, if an optimal solution $s^*$ at time $a$ remains feasible at time $a+\Delta$, with $\Delta > 0$ then $s^*$ is also an optimal solution at time $a+\Delta$. In other word, an optimal schedule can be optimal for a range of elements of $T_i$. We model the problem as follows:

Let $J_K^i$ be the set of jobs where each job represents the delivery of a customer $K \subseteq K_i$ by drone from $i \in V$. We note $p_t$ the processing time of the job $t \in J_K^i$, $v_t$ the customer it delivers and $\bar{d}_t$ its deadline. We set $p_t = c_{i,v_t}^d + s_{v_t}^d + c_{v_t,i}^d$. To keep the consistency with the definition of deadline in the 2E-VRP-D, we set $\bar{d}_t = \bar{d}_{v_t} + s_{v_t}^d + c_{v_t,i}^d - a$.

As mentioned earlier, the set of solutions for a given $i \in V, d \in D, K \subseteq K_i$ is noted $S_{iK}^d$. For $s \in S_{iK}^d$, let $m_s$ be its makespan, i.e. its completion time minus its starting time, $C_t^s$ be the completion time of a job $t \in J_K^i$ and $\bar{a}_s$ be its maximal starting time i.e. $\bar{a}_s = \min_{t \in J_K^i} (\bar{d}_t - C_t^s)$. Hence, $s$ remains optimal and feasible for all $a' \in [a, \bar{a}_s]$ and we don't have to calculate other optimal solutions in the interval.

### 4.2.2 Scheduling algorithm

We use a dynamic program to find all the solutions in a $S_{iK}^d$, $i \in V, d \in D, K \subseteq K_i$. Because we solve a whole set $S_{iK}^d$, we have to find the solution for all $a \in T_i$. Hence we must find the Pareto front minimizing the makespan and maximizing the maximal starting time. The idea of the dynamic program is to check all possible combinations of solution of problems with fewer machines and subsets of jobs.

First we solve $S_{iK}^1$ $\forall i \in V_+$, $K \in K_i$ using the Jackson's rule, also known as the Earliest Deadline First (Jackson, 1955), (Stankovic et al., 1998), and checking that the schedule is feasible; otherwise there is no solution in $S_{iK}^1$. Note that $|S_{iK}^1| \le 1$. Let $K \subseteq K_i$, $s \in S_{iK}^d$ and $Z$ be the set of jobs done on one of the machines of $s$. $\exists u \in S_{iZ}^1$ and $\exists v \in S_{iK \setminus Z}^{d-1}$ such that by combining $u$ and $v$, we obtain $s$. We write $s = (u, v)$. We have, $m_{(u,v)} = \max(m_u, m_v)$ and $\bar{a}_{(u,v)} = \min(\bar{a}_u, \bar{a}_v)$.

$$S_{iK}^d = \{(u, v) | u \in S_{iZ}^1, v \in S_{iK \setminus Z}^{d-1}, \forall s \in S_{iK}^d, m_s < m_{(u,v)} \Rightarrow \bar{a}_s < \bar{a}_{(u,v)}(1), \forall Z \subseteq K\} \quad (16)$$

The condition (1) ensure that $S_{iK}^d$ contains the Pareto front minimizing $m_s$ and maximizing $\bar{a}_s$. To avoid symmetries, we set the machine running set $Z$ to be the bottleneck machine inducing the makespan. This enable us to only check sets $Z$ such that $\sum_{t \in Z} p_t \ge \frac{1}{d} \times \sum_{t \in K} p_t$, because the makespan of a feasible solution can not be lower than $\frac{1}{d} \times \sum_{t \in K} p_t$.

$|\{Z \subseteq K\}| = 2^{|K|}$, hence the complexity of this dynamic program is exponential in $|K|$. The algorithm is output sensitive, as the output size can be exponential over the input size. The complexity of this algorithm is polynomial, $O(n^3)$ where $n$ is the size of its output. ($n$ schedules to calculate, for each we check
if and only if route $r \in \mathcal{R}^d$ follows arc $(i', j)$. Note that the prime is omitted in the notation without ambiguity, as there is no arc between $i$ and $j$ with $i, j \in V_+.

<figure><img src="image_2.png" /><figcaption>Figure 1: Example of a subproblem graph $\mathcal{G}^d$</figcaption></figure>

Figure (1) shows an example of a subproblem graph. We represented only three arcs from $i \in V$ to $i' \in V'$ for clarity, but there should be an arc for all $s \in S_{iK}^d, \forall K \subseteq N_i$.

### 5.1. Resources for the pricing subproblem

In our model, we use 2 resources, the duration $\beta_d$ of a route and the capacity $\beta_c$ of a truck, and we define a pricing subproblem for each $d \in \{0, ..., D\}$. The objective function is the minimization of the reduced cost. The reduced cost of a route $r$ with a cost $c_r$ is given by: $\bar{c}_r = c_r - \sum_{i \in V} \alpha_{ir} \times \pi_i - \lambda_0^v - d_r \times \lambda_0^d$

where $\alpha_{ir}$ is equal to one if and only if customer $i$ is visited in the route $r$, $d_r$ the number of drones used in route $r$, $\pi_i$ are the dual variables associated with constraints (2), $\lambda_0^v$ is the dual variable associated with constraint (3) and $\lambda_0^d$ is dual variable associated with constraint (4).

Let $[l_{a,r}, u_{a,r}]$ be the resource feasibility interval, for each $a \in \mathcal{A}^d$ and $r \in \{\beta_c, \beta_d\}$, i.e. the accumulated resource $r$ has to lie within the interval $[l_{a,r}, u_{a,r}]$ to travel through arc $a$. Let $q_{a,r}$ be the consumption of resource $r$ on arc $a$. For the resource capacity, we set:
each combination of schedule with 1 machine and with n-1 machine: $O(n^2)$)

### 4.2.3 Criteria for valid arcs elimination

In this section, we introduce two criteria that allow us to perform preprocessing to eliminate infeasible and dominated arcs in $A^{dD}$, from a subproblem graphs $G$, $d \in D$.

* Any arc in $S_{iK}^d, d \in D, i \in V, K \subseteq K_i$ is infeasible if the total payload is higher than the truck capacity, i.e. $q_i + d \times q_1^d + \sum_{k \in K} q_k > Q$ because the scheduling solution is not feasible in the routing problem, thus the arc representing this schedule is infeasible.

* Let a truck with $d$ drones be parked at customer $i$. If there is a faster route (combining truck and drones) delivering set of customers $K$ than delivering all of them by drone from $i$ (i.e. using an appropriate drone schedule in $S_{iK}^d$), then we know that the corresponding solution in $S_{iK}^d$ will not be used in any optimal solution and we can remove it safely.

<figure><img src="image_4.png" /><figcaption>Figure 2: Example of an application of SPR criterion</figcaption></figure>

In the figure 2 we see an example of an application of the SPR criterion. In the left side a single drone delivers customers $j,k,l$ from customer $l$. In the right side a drone delivers $j$ from $i$ then the truck moves to customer $k$ and then the drone delivers $l$ from $k$. The duration of the left partial solution is longer than the duration of the right side. It is not sufficient to remove $S_{i\{jkl\}}^1$, indeed we also need to verify that the right solution is feasible regardless of where the vehicle moves. This will be explained in detail in the following section. The SPR criterion generalize and tackle such solutions.
$$l_{a,\beta_c} = 0 \qquad \forall a \in \mathcal{A}^d \qquad (15)$$

$$u_{a,\beta_c} = Q - d \times q_1^d \qquad \forall a \in \mathcal{A}^d \qquad (16)$$

$$q_{a,\beta_c} = \begin{cases} 0 & \forall a = (i, 0), i \in V' \\ q_j & \forall a = (i, j), i \in V'_+, j \in V \\ \sum_{k \in K} q_k & \forall a = (i, i')^s, i \in V, i' \in V', K \in \mathcal{N}_i, s \in S_{iK}^d \end{cases} \qquad (17)$$

The truck starts with a remaining capacity of $Q-d \times q_1^d$. The remaining capacity of a truck have to lie between 0 and $Q-d \times q_1^d$ on every arc of the graph. The remaining capacity consumption of an arc is set to the sum of the demands of the customers served, or 0 if the destination is the depot.

Let $\Delta_s$ be the maximal starting time of a solution $s$ that will be formally defined in section (5.2.1). For the resource duration, we set:

$$l_{a,\beta_d} = 0 \qquad \forall a \in \mathcal{A}^d \qquad (18)$$

$$u_{a,\beta_d} = \begin{cases} \bar{d}_j & \forall a = (i, j), i \in V'_+, j \in V_+ \\ \Delta_s & \forall a = (i, i')^s, i \in V, i' \in V', K \in \mathcal{N}_i, s \in S_{iK}^d \end{cases} \qquad (19)$$

$$q_{a,\beta_d} = \begin{cases} c_{i,j}^v & \forall a = (i, j), i \in V'_+, j \in V_+ \\ \max(s_i^v, m_s) & \forall a = (i, i')^s, i \in V, i' \in V', K \in \mathcal{N}_i, s \in S_{iK}^d \end{cases} \qquad (20)$$

The truck starts at time 0. A truck arc serving a customer can only be used if the accumulated duration lies between 0 and the customer's deadline. A drone arc serving a set of customers can only be used if the accumulated duration lies between 0 and the drone schedule's maximum starting time. The duration consumption of a drone arc is the maximum between the truck's service time at customer $i$ and the duration of the drone schedule, as they are operated in parallel.

## 5.2. Drone schedules enumeration

The objective is to enumerate all optimal drone schedules from all customers $i \in V$, for all number of drones on the truck $d \in \{1, ..., D\}$, and for all sets of available customers, $K \subseteq \mathcal{N}_i$. Note that the number of subsets of $\mathcal{N}_i$ is exponential. Because of customers' deadlines, the optimal drone schedule serving the set of customers $K$ depends on the truck arrival time at $i$. We note $T_i = [c_{0,i}^v, \bar{d}_i]$ the interval containing all possible truck arrival times at customer $i$. Note that our approach does not require time discretization.
More formally, for a given $i \in V_+$, $d \in D$, $K \in \mathcal{K}_i$ and $j \in V \setminus K \cup \{i\}$, let SPR$(d, i, K, j)$ be the duration of the shortest partial route from $i$ to $j$ using a truck equipped with $d$ drones and delivering either by drone or by truck the customers of set $K$. If for any feasible $j$, SPR$(d, i, K, j) < m_s + c_{i,j}^v$, $s$ will not be used in an optimal solution. We consider the following set of feasible $j$:

$$ \left\{ j \in V_+ \setminus K \cup \{i\} \mid (c_{0,i}^v + \max(s_i^v, m_s) + c_{i,j}^v < \bar{d}_j) \wedge \left( q_i + d \times q_1^d + q_j + \sum_{k \in K} q_k \le Q \right) \right\} \quad (17) $$

In our computational experiment and because checking this criterion would be too time consuming, we only consider special cases where the truck is only allowed to visit a single customer of $K$, and the others must be delivered by drone. This criterion let us remove up to 90% of the arcs in $\mathcal{G}^1$, 61% in $\mathcal{G}^2$ and 6% in $\mathcal{G}^3$.

# 5 Branch-cut-and-price algorithm

Because it would be impractical to enumerate the set $\mathcal{R}^d, d \in D$ due to its exponential size over the number of customers we solve this formulation using a column and cut generation approach. We use a similar branch-and-cut-and-price (BCP) as in (Pessoa et al., 2020), with the following features. The BCP algorithm iteratively generates new variables (columns) to improve a linear programming (LP) relaxation of an integer program. Branching decisions are made based on the LP solution and cutting planes are added to strengthen the relaxation. Two primal heuristics are used in branch and bounds nodes. The first heuristic consists in solving the restricted master problem using a MIP solver and a time limit. The second heuristic is called the diving heuristic, where columns with largest fractional values are iteratively set to one and the remaining problem is solved using column generation.

## 5.1 Column generation

We use three-stage column generation. The pricing is solved using heuristic labelling algorithms in the two first stages as in (Sadykov et al., 2021), and using the exact labelling dynamic algorithm as in (Pessoa et al., 2020) and (Sadykov et al., 2021). The first two stages generate at most 30 columns, and the last stage at most 150. Automatic dual price smoothing stabilization is used, as in (Pessoa et al., 2018) to help the convergence of the process. The bucket-arc elimination procedure is also used, removing arcs proven not to belong to any optimal solution. After each call, we use a route enumeration technique similar to (Baldacci et al., 2008), which tries to enumerate all improving routes for a subproblem, and if it succeeds, the future pricing problems are solved by inspection of the enumerated routes. Furthermore, if the number of enumerated routes in all subproblem is less than a certain threshold (we use 5000 as threshold), all
### 5.2.1. Scheduling problem

Each drone can only deliver a single customer per round-trip and is automatically fully recharged when reaching the truck. Consequently, we can model the drone routing problem for a truck with $d$ drones arriving at customer $i$ at time $t_a \in T_i$ and serving set $K \subseteq \mathcal{N}_i$ customers as a scheduling problem with $d$ parallel identical machines, $|K|$ jobs with deadlines and minimizing the makespan. The arrival time $t_a$ can be treated as a common release date for all jobs and eliminated by subtracting $t_a$ from each job's deadline. This operation is included in the job deadline formula 22. Therefore, the notation of this scheduling problem is $P|\bar{d}_j|C_{\max}$ (Graham et al., 1979). In this problem, waiting times between jobs can be eliminated without compromising the quality of a solution. Therefore, solutions involving waiting time can be omitted.

We model the problem as follows:

Let $J_K^i$ be the set of jobs where each job corresponds to the delivery of a customer in $K \subseteq \mathcal{N}_i$ by drone from $i \in V$. Each job corresponds to a customer's delivery. The processing time of the job is defined by the period of drone unavailability caused by the delivery, which includes the setup time, travel times, and service time. The job's deadline is derived from the customer's delivery deadline, though it is adjusted to align with the definition of processing time. For $t \in J_K^i$, let:

* $p_t$ be its processing time

* $\bar{d}_t$ be its deadline

* $v_t$ be the customer whose delivery it represents.

Formally, we have the following equations:

$$p_t = o_{v_t} + c_{i,v_t}^d + s_{v_t}^d + c_{v_t,i}^d \quad (21)$$

$$\bar{d}_t = \bar{d}_{v_t} + s_{v_t}^d + c_{v_t,i}^d - t_a \quad (22)$$

To solve this problem, a distinction must be made depending on the number of machines.

* When considering only one machine, this scheduling problem can be solved in polynomial time using the Jackson's rule, also known as the Earliest Deadline First (Jackson, 1955), (Stankovic et al., 1998). Therefore, processing the jobs ordered by increasing deadlines yields an optimal solution. Modifying $t_a$ does not affect the sequence of jobs sorted by increasing deadlines. Consequently, for a given $t_a$, either this solution is optimal or no feasible solution exists.

* When considering more than one machine, $P||C_{\max}$ is known to be NP-hard Błażewicz et al. (2007), therefore, $P|\bar{d}_j|C_{\max}$ is NP-hard.

### 5.2.2. Dynamic program computing Pareto front

As mentioned earlier, the set of schedules for a given $i \in V$, $d \in \{0, ..., D\}$, $K \subseteq \mathcal{N}_i$ is noted $S_{iK}^d$. To be more precise, $S_{iK}^d$ is the Pareto front minimizing the
the routes are added to the restricted master problem, and solved by the MIP solver.

## 5.2 Cutting planes

### 5.2.1 Adapted rounded capacity cuts

The rounded capacity cuts were first introduced in (Laporte and Nobert, 1983) for the constrained vehicle routing problem (CVRP). Let $C \subset V$ and let $h_r$ be the number of times $r \in \mathbb{R}^d$ visits C i.e.

$$h_r = \frac{1}{2} \times \sum_{i \in V_+ \setminus C} \sum_{j \in C} b_{\min(i,j), \max(i,j)}^r + \sum_{i \in V_+ \setminus C} \sum_{K \in \mathcal{K}_i: K \cap C \neq \emptyset} \sum_{s \in S_{iK}^d} a_s^r \quad \forall d \in D, r \in \mathbb{R}^d \quad (18)$$

These valid inequalities can be formulated as follows:

$$\sum_{d \in D} \sum_{r \in \mathbb{R}^d} \lambda_r \times h_r \geq \left\lfloor \frac{\sum_{i \in C} q_i}{Q} \right\rfloor, \quad C \subset V \quad (19)$$

However, because in this problem the drones take a large amount of the truck's capacity (about 20% per drone in the instances), we must take them into account. Let $k = \left\lfloor \frac{Q}{q_1^d} \right\rfloor$, and adapt the rounded capacity cuts as follows:

$$\sum_{d \in D} \sum_{r \in \mathbb{R}^d} \lambda_r \times h_r \times (k-d) \geq \left\lfloor \frac{\sum_{i \in C} q_i}{q_1^d} \right\rfloor, \quad C \subset V \quad (20)$$

We will refer to these cuts as adapted rounded capacity cuts (ARCC). To separate the ARCC we use the greedy construction heuristic from (Lysgaard et al., 2004).

<caption>**Algorithm 1** Separation for the adapted rounded capacity cuts</caption>

Require: $t \in \mathbb{R}$

for $i \in V$ do
  $C \leftarrow \{i\}$
  while $C \neq V$ do
    $C \leftarrow C \cup \{\arg \max_{j \in V \setminus C} (\text{violation}(C \cup \{j\}))\}$
    if $\text{violation}(C) > t$ then
      addCut($C$)

We describe the separation algorithm in Algorithm (1), where $t$ is a user-defined threshold, the function $\text{violation}(C)$ calculates the difference between the right-hand side and the left-hand side in constraints (20) for set $C$ and $\text{addCut}(C)$ add the valid inequality to the problem. The algorithm can be explained as follows: starting with a single customer in the set, while not all customers are
makespan and maximizing the maximal starting time.

This section presents a dynamic program that solves for the Pareto front corresponding to a set $S_{iK}^d$.

For $s \in S_{iK}^d$, let:

* $m_s$ be its makespan.

* $C_t^s$ be the completion time of job $t \in J_K^i$.

* $\Delta_s$ be its maximal starting time, i.e., $\Delta_s = \min_{t \in J_K^i} (\bar{d}_t - C_t^s)$.

The algorithm computes $S_{iK}^1$ for all $i \in V$ and for all $K \in \mathcal{N}_i$. As explained in the previous section, when using a single machine, we only have to consider the solution given by Jackson's rule. Consequently, solving the problem solely for truck arrival time $t_a = c_{0,i}^v$, the earliest time a truck can arrive at customer $i$, is sufficient to obtain $S_{iK}^1$. Note that if $\Delta_s$ is negative, $s$ is unfeasible and does not belong to $S_{iK}^1$.

The dynamic programming method exploits the fact that we need to compute solutions for all $K \subseteq \mathcal{N}_i$ and all $d \in \{1, ..., D\}$. Thus, calculated optimal solutions can be combined to generate feasible solutions for problems involving additional machines.

For instance, let $X, Z \subseteq \mathcal{N}_i$ such that $X \cap Z = \emptyset$. Let $d \in \{2, ..., D\}$. Let $x^* \in S_{iX}^1$ and $z^* \in S_{iZ}^{d-1}$ be already calculated optimal solutions. Then, assigning $x^*$ on one machine and $z^*$ on $d-1$ machines, we obtain a feasible solution for the scheduling problem with $d$ machines and the set of jobs $X \cup Z$. We note a solution obtained from combining two solutions $x^*$ and $z^*$ as $(x^*, z^*)$. Although solution $(x^*, z^*)$ might not be an optimal solution for the problem, combining solutions in such a way ensures that we enumerate only promising feasible solutions.

More generally, we write $s = (u, v)$, the solution obtained from combining solution $u$ and solution $v$. We have, $m_{(u,v)} = \max(m_u, m_v)$ and $\Delta_{(u,v)} = \min(\Delta_u, \Delta_v)$.

Let $H_{iK}^d$ be the set of all such combinations. To avoid symmetries, we assume the single-machine solution to be the bottleneck machine inducing the makespan. Hence, we have:

$$H_{iK}^d = \{(u, v) | u \in S_{iZ}^1, v \in S_{iK \setminus Z}^{d-1}, m_u \ge m_v\}$$

Minor improvements can reduce computation time. For example, the single-machine solution $u$ will not appear in $H_{iK}^d$ if it does not verify the equations $m_u \ge \frac{1}{d} \times \sum_{t \in J_K^i} p_t$ and $m_u \ge \max_{t \in J_K^i} p_t$, as the condition $m_u \ge m_v$ would never be verified.

A solution $s$ dominates a solution $s'$ if its makespan is lower and its maximum
in the set, we add a customer that is not already in the set and that maximizes the violation. If the violation of the set exceeds the threshold, we add the corresponding cut to the linear program. We generate at most 100 cuts per round.

### 5.2.2 Limited memory rank-one cuts

We use limited memory rank-one cuts (lm-R1C) similar to (Pecin et al., 2017), (Pessoa et al., 2020). The R1C uses a Chvátal-Gomory rounding on a weighted sum of a subset of constraints (2). When the coefficients are all equal, these cuts are called Subset-Row cuts (Jepsen et al., 2008).

The R1C are defined as follows. Let $C \subset V_+$, with multipliers $0 < \rho_i < 1 \quad \forall i \in C$ and let $h_i^r$ be the number of time the vertex $i$ is visited in route $r$:

$$h_i^r = \frac{1}{2} \sum_{j \in C} b_{\min(i,j), \max(i,j)}^r + \sum_{j \in V_+ \setminus C} \sum_{K \in \mathcal{K}_j: s \in S_{jK}^d} a_s^r \quad \forall d \in D, r \in \mathcal{R}^d \quad (21)$$

The following equation is a valid inequalities:

$$\sum_{d \in D} \sum_{r \in \mathcal{R}^d} \left[ \rho_i \sum_{i \in C} h_i^r \right] \lambda_r \le \left[ \sum_{i \in C} \rho_i \right] \quad (22)$$

Because the R1C are non-robust constraints, we use the lm-R1C instead, where the set $C$ is paired with a memory set $M$ such that $C \subseteq M \subseteq V_+$. $\left[ \rho_i \sum_{i \in C} h_i^r \right]$ in equation (22) is replaced by $\alpha(C, M, \rho, P)$ and is calculated during the pricing computation as in algorithm 2.

**Algorithm 2** $\alpha(C, M, \rho, P = (v_0 = 0, v_1, ..., v_k, v_k + 1 = 0'))$
$\alpha \leftarrow 0, s \leftarrow 0$
**for** $i = 1$ **to** $k$ **do**
&nbsp;&nbsp;&nbsp;&nbsp;**if** $(v_{i-1}, v_i) \notin M$ **then**
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**s** $\leftarrow 0$
&nbsp;&nbsp;&nbsp;&nbsp;**if** $v_i \in C$ **then**
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**s** $\leftarrow s + \rho_{v_i}$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $s \ge 1$ **then**
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**s** $\leftarrow s - 1, \alpha \leftarrow \alpha + 1$
&nbsp;&nbsp;&nbsp;&nbsp;**return** $\alpha$

If the partial path $P$ leaves the set $M$, all previous visits are forgotten, leading to an easier to compute but smaller coefficient. Note that if $M = V$, the lm-R1C are equivalent to the R1C. To separate those cuts, we find a violated constraint, and we calculate its minimal memory set such that the limited memory R1C has the same violation.
starting time is higher. Formally, we define the solutions dominance as follows:

$$s \succeq s' \Leftrightarrow (m_s \le m_{s'} \land \Delta_s \ge \Delta_{s'})$$

$$s = s' \Leftrightarrow (m_s = m_{s'} \land \Delta_s = \Delta_{s'})$$

$$s \succ s' \Leftrightarrow (s \succeq s' \land s \neq s')$$

We now describe the dynamic program's formula to calculate a set $S_{iK}^d$.

$$S_{iK}^d = \{s \in H_{iK}^d, \forall s' \in S_{iK}^d \setminus \{s\}, s' \not\succeq s\} \quad (23)$$

We recall that $S_{iK}^d$ contains the Pareto front minimizing $m_s$ and maximizing $\Delta_s$. When two solutions on the Pareto front have identical objective values, either can be chosen without loss of generality. We remove one arbitrarily. Algorithm 3 describes the procedure to calculate $S_{iK}^d$.

**Algorithm 3 computeS(d ∈ {2, ..., D}, i ∈ V, K ⊆ N_i)**

1: $S \leftarrow \{\}$

2: **for** $X \subset K$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;**for** $x \in S_{iX}^1$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Z $\leftarrow K \setminus X$

5: &nbsp;&nbsp;&nbsp;&nbsp;**for** $z \in S_{iZ}^{d-1}$ **do**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $\nexists s \in S | s \succ (x, z)$ **then**

7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;S $\leftarrow S \cup \{(x, z)\}$

8: **return** S

Note that with specific jobs' deadlines and processing times, the number of solutions in a set $S_{iK}^d$ can also be exponential. In practice, the size of the sets remained manageable throughout all processed instances.

The algorithm is exponential in the size of K. However, as the output size can be exponentially larger than the input size, the algorithm is output-sensitive. The complexity of the algorithm computing algorithm 3 for each $i \in V, d \in \{2, ..., D\}$ and $K \subseteq N_i$ is polynomial in the size of its output.

### 5.2.3. Pareto front enumeration example

In this section, we present a simple example of enumeration using the dynamic program. Figure 2 represents a customer and its drone-reachable neighborhood. We consider a truck parked at customer i. Customers 1, 2, and 3 belong to the drone-reachable neighborhood of i ($N_i = \{1, 2, 3\}$). Note that for readability reasons, different notations are used for the customers in the drone-reachable neighborhood and for the customer where the truck is parked. Figure 2 represents the one-way flight duration on each arc. Deadlines are provided above each customer node in $N_i$. We consider $c_{0,i}^v$ to be zero, the service time for drones to be one, and the setup time to be zero for all customers.
## 5.3 Branching

We obtain branching candidates from: the total number of drones ($\sum_{d \in D} d \times u^d$), the total number of trucks ($\sum_{d \in D} u^d$), the number of trucks with $d$ drones (variables $u$) and the truck travel variables (variables $x$).

We use a multi-phase strong branching procedure, similar to (Pecin et al., 2017), to choose the most promising branching candidate. The idea is to spend more time evaluating branching variables in the lowest depth of the branch-and-bound tree where each selection has a greater impact on the overall time, and spend less time as the depth increases, taking advantage of the history of previous evaluations.

In the first phase, half the candidates are chosen from history from previous calls of the strong branching procedure, while the others are chosen in a balanced way from all the branching strategies based on the distance between its value and the closest integer value. The candidates with the largest distances are selected.

In the second phase, selected candidates from the first phase are evaluated by solving the current restricted master linear program modified for each candidate, without column generation. We select the candidates using the product rule (Achterberg, 2007).

In the third phase, selected candidates from the second phase are evaluated by solving the restricted master linear program modified for each candidate, with heuristic column generation but without cut generation. We select the best candidate using the product rule.

# 6 Heuristic branch-cut-and-price

In this section, we propose a heuristic branch-cut-and-price based on the reduction of the size of the subproblems graphs. It has four parameters :

$\Gamma_h$ — The maximum number of drones a truck can carry.

$w_h$ — The maximum number of customers that can be delivered by drone from a customer in a scheduling.

$k_h$ — The maximum number of drone-reachable neighbors kept. Only the $k_h$ closest drone-reachable neighbors are kept.

$p_h$ — The number of schedule with earliest deadline kept for each $S_{iK}^d$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Example customer *i* and its drone-reachable neighborhood</figcaption></figure>

In this example, we calculate $S_{i\{1,2,3\}}^2$ using formula 23. In other words, we calculate the optimal drone schedules delivering a set of customers $\{1, 2, 3\}$ using two drones from a truck parked at customer *i* for any truck arrival date.

We note a schedule with one machine as the order of the processed jobs between parentheses. For instance: $(1, 3)$ is the schedule with one machine, processing job 1, then job 3. We note a schedule with more than one machine as the order of each machine's processed jobs separated by a "|". For instance, $(3|1, 2)$ is the schedule with two machines, where the first machine processes job 3 and the second machine processes job 1 and then job 2.

The first step of the algorithm is to model the problem as a scheduling problem. To convert the drone travel from customer *i* to customer $v_t$ to a job *t*, we use equation 21 and equation 22. In table 3, the instance's data is summarized, and jobs are presented.

<table><caption>Table 3: Drone round-trips conversion to jobs</caption><thead><tr><th>t</th><th>v<sub>t</sub></th><th>o<sub>v<sub>t</sub></sub></th><th>c<sup>d</sup><sub>i,v<sub>t</sub></sub> = c<sup>d</sup><sub>v<sub>t</sub>,i</sub></th><th>s<sub>v<sub>t</sub></sub></th><th>d&#773;<sub>v<sub>t</sub></sub></th><th>p<sub>t</sub></th><th>d&#773;<sub>t</sub></th></tr></thead><tbody><tr><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>3</td><td>3</td><td>5</td></tr><tr><td>2</td><td>2</td><td>0</td><td>1.5</td><td>1</td><td>4.5</td><td>4</td><td>7</td></tr><tr><td>3</td><td>3</td><td>0</td><td>2</td><td>1</td><td>7</td><td>5</td><td>10</td></tr></tbody></table>

The second step is to sort the jobs by increasing deadlines in the ordered set $K \subseteq N_i$. For clarity, this order is the lexicographic order here. The third step is to solve $S_{iK}^1$ with the Earliest Deadline First rule.
$\Gamma_h$ decreases the number of subproblems. $w_h$ and $k_h$ reduce the number of drones arcs in each subproblems graphs. Note that the difference between $w_h$ and $k_h$ is that $k_h$ control which customer can be delivered by drone from a parked vehicle, while $w_h$ controls how many customers can be delivered by drone from a parked vehicle. Finally, $p_h$ reduces number of schedule and therefore the number of arcs in subproblems with multiple drones.

Let $\mathcal{K}_i^k$ be the set of the $k$ closest to $i$ customers in $\mathcal{K}_i$. More formally, $\mathcal{K}_i^k \subseteq \mathcal{K}_i$, $|\mathcal{K}_i^k| = k$ and $\forall k \in \mathcal{K}_i^k, k' \in \mathcal{K}_i \setminus \mathcal{K}_i^k, c_{ik}^d \le c_{ik'}^d$. Let $S_{iK}^{dp}$ be the set of the $p_h$ solutions from $S_{iK}^d$ with earliest deadlines. To apply these rules, we transform each subproblems graphs $\mathcal{G}^d = (\mathcal{V}^d, \mathcal{A}^d), d \in [1, ..., \Gamma]$ with $\mathcal{A}^d = A^T \cup A^D$ from the exact branch-cut-and-price subproblems graphs as follows: we create the graphs $\mathcal{G}_h^d = (\mathcal{V}_h^d, \mathcal{A}_h^d), d \in [1, ..., \Gamma_h]$, where $\mathcal{V}_h^d = \mathcal{V}^d$ and $\mathcal{A}_h^d = A^T \cup \{(i, i')^s \in A^D, \forall K \subseteq \mathcal{K}_i^{kh}, s \in S_{iK}^{dp_h} : |K| \le w_h\}$

# 7 Computational experiments

## 7.1 Configuration

Our algorithms such as the dynamic program for the enumeration, the separation algorithm for the ARCC and the SPR criterion are coded in C++. We also used:

* IBM ILOG CPLEX 20.1 as the LP solver in the relaxed restricted master problem and as the solver for the enumerated MIPs.

* The generic branch-cut-and-price code from (Sadykov and Vanderbeck, 2021), with automatic dual price smoothing stabilization from (Pessoa et al., 2018), primal heuristics from (Sadykov et al., 2019). We also use the labelling algorithm code from (Sadykov et al., 2021).

The experiments presented in this paper were carried out:

* Using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr).

* Running on 2.6 GHz 2x 18-core Cascade Lake Intel Xeon Skylake Gold 6240 CPU using a single thread with 5.3 Go memory.

## 7.2 Instances

To measure the performance of our algorithm, we used and extended the sets of instances from (Zhou et al., 2023), which is composed of modified instances from (Chen et al., 2021), and from (Solomon, 1987). In all these instances, the trucks and the drones follow the same distance metric. We also used clustered instances from (Gehring and Homberger, 1999) to experiment on larger instances, as they
Formulas for the makespan and the maximal starting time are as follows:

$$m_{s_K} = \sum_{t \in J_K^i} p_t$$
$$\Delta_{s_K} = \min_{J_K^i} \left( \bar{d}_t - \sum_{\substack{t' \in J_K^i \\ C_{t'}^{s_K} \le C_t^{s_K}}} p_{t'} \right)$$

$$S_{i\{1\}}^{1} = \{(1) : m_{(1)} = 3, \Delta_{(1)} = 2\}$$

$$S_{i\{2\}}^{1} = \{(2) : m_{(2)} = 4, \Delta_{(2)} = 3\}$$

$$S_{i\{3\}}^{1} = \{(3) : m_{(3)} = 5, \Delta_{(3)} = 5\}$$

$$S_{i\{1,2\}}^{1} = \{(1,2) : m_{(1,2)} = 7, \Delta_{(1,2)} = 0\}$$

$$S_{i\{1,3\}}^{1} = \{(1,3) : m_{(1,3)} = 8, \Delta_{(1,3)} = 2\}$$

$$S_{i\{2,3\}}^{1} = \{(2,3) : m_{(2,3)} = 9, \Delta_{(2,3)} = 1\}$$

$$S_{i\{1,2,3\}}^{1} = \{(1,2,3) : m_{(1,2,3)} = 12, \Delta_{(1,2,3)} = -2\}$$

$S_{i\{1,2,3\}}^{1}$ is unfeasible as $\Delta_{1,2,3} < 0$.

The fourth step is to compute $S_{i\{1,2,3\}}^{2}$ using the formula 23.

$$\{(u,v)|u \in S_{iZ}^{1}, v \in S_{i\{1,2,3\} \setminus Z}^{1} \quad \forall Z \subseteq \{1,2,3\}\} = \{(1), (2,3), \\ ((2), (1,3)), \\ ((3), (1,2))\}$$

With this set, each combination of solutions is calculated using the values in step four. For instance ((1), (2, 3)), 1 is processed on the first machine, while 2 and 3 are processed on the remaining machine, as $d - 1 = 1$. Their makespan and maximal starting time are calculated as follows:

$$m_{(1|2,3)} = \max(m_{1},m_{2,3}) = 9$$

$$\Delta_{(1|2,3)} = \min(\Delta_{1}, \Delta_{2,3}) = 1$$

$$m_{(2|1,3)} = \max(m_{2},m_{1,3}) = 8$$

$$\Delta_{(2|1,3)} = \min(\Delta_{2}, \Delta_{1,3}) = 2$$

$$m_{(3|1,2)} = \max(m_{3},m_{1,2}) = 7$$

$$\Delta_{(3|1,2)} = \min(\Delta_{3}, \Delta_{1,2}) = 0$$
are an extension of the Solomon (1987) instances. The maximal drone capacity is set to 10 kg as in (Chen et al., 2021), the service times with drone are set to $s_i^d = [0.5s_i^v]$, $i \in V$.

In the instances proposed in (Chen et al., 2021), an additional $[0.2 \times |V|]$ customers are randomly selected to be unavailable for drone delivery. Table 3 describes the different data sets. For each data set, we give the name of the instances (Instances) where $X$ is the id of a particular instance of the set, the number of instances in the set (Size), the customer positioning (Pos), the distance metric (Dist), the original authors of the instances (Authors), the capacity of the vehicle ($Q$), the maximum number of drones allowed on a vehicle ($\Gamma$), the total number of vehicles ($k_v$) and the total number of drones ($k_d$). In column (Pos), "rand" stands for "random", and "clust" for "clustered". In column (Dist), "eucl" stands for "euclidean".

<table><caption>Table 3: Instances sets and their characteristics</caption><thead><tr><th>Instances</th><th>|V|</th><th>Size</th><th>Pos</th><th>Dist</th><th>Authors</th><th>Q</th><th>Γ</th><th>k<sub>v</sub></th><th>k<sub>d</sub></th></tr></thead><tbody><tr><td colspan="10" style="text-align:center; font-weight:bold;">Small-Size Instances</td></tr><tr><td>Cardiff_10_X</td><td>10</td><td>20</td><td>rand</td><td>real</td><td>Chen et al. (2021)*</td><td>170</td><td>2</td><td>4</td><td>4</td></tr><tr><td>Cardiff_15_X</td><td>15</td><td>20</td><td>rand</td><td>real</td><td>Chen et al. (2021)</td><td>170</td><td>2</td><td>5</td><td>5</td></tr><tr><td>Cardiff_25_X</td><td>25</td><td>20</td><td>rand</td><td>real</td><td>Chen et al. (2021)</td><td>200</td><td>3</td><td>5</td><td>10</td></tr><tr><td>c10X_35</td><td>35</td><td>9</td><td>clust</td><td>eucl</td><td>Solomon (1987)</td><td>250</td><td>3</td><td>5</td><td>10</td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">Medium-Size Instances</td></tr><tr><td>Cardiff_40_X</td><td>40</td><td>20</td><td>rand</td><td>real</td><td>Chen et al. (2021)*</td><td>200</td><td>3</td><td>5</td><td>10</td></tr><tr><td>c10X_100</td><td>100</td><td>9</td><td>clust</td><td>eucl</td><td>Solomon (1987)</td><td>250</td><td>3</td><td>14</td><td>28</td></tr><tr><td>120_C1_2_X</td><td>120</td><td>10</td><td>clust</td><td>eucl</td><td>Gehring and Homberger (1999)</td><td>250</td><td>3</td><td>40</td><td>60</td></tr><tr><td>130_C1_2_X</td><td>130</td><td>10</td><td>clust</td><td>eucl</td><td>Gehring and Homberger (1999)</td><td>250</td><td>3</td><td>40</td><td>60</td></tr><tr><td>140_C1_2_X</td><td>140</td><td>10</td><td>clust</td><td>eucl</td><td>Gehring and Homberger (1999)</td><td>250</td><td>3</td><td>40</td><td>60</td></tr><tr><td>150_C1_2_X</td><td>150</td><td>10</td><td>clust</td><td>eucl</td><td>Gehring and Homberger (1999)</td><td>250</td><td>3</td><td>40</td><td>60</td></tr></tbody></table>

Note that the set with 10 customers and the set with 40 customers are marked with '*' as they are not in (Chen et al., 2021), but are derived from sets with larger number of customers from (Chen et al., 2021), namely 15 and 50, by keeping only the 10 and 40 first customers.

The instances $c10X_N$ are the 100-customers clustered instances called $c10X$ in (Solomon, 1987), restricted to the $N$ first customers. We categorize instances with 35 or fewer customers as small, while instances with more than 35 customers will be referred to as medium. The instances $N_C1_2_X$ are the 200-customer clustered instances called $C1_2_X$ in (Gehring and Homberger, 1999), restricted to the first $N$ customers.

## 7.3 Results

In this section we present all the results obtained by setting the computation time limit at three hours per instance. Unless specified otherwise, the ARCC
<figure><img src="image_4.png" /><figcaption>Figure 3: Pareto front for set $S_{i\{1,2,3\}}^2$</figcaption></figure>

Figure 3 represents the Pareto front for $S_{i\{1,2,3\}}^2$ in the objective space. The red arrows represent the optimization direction for each objective. We observe that $(1|2, 3)$ is dominated by $(2|1, 3)$ as $m_{(1|2,3)} > m_{(2|1,3)}$ and $\Delta_{(1|2,3)} < \Delta_{(2|1,3)}$. Therefore, $(1|2, 3)$ does not belong in $S_{i\{1,2,3\}}^2$. The two remaining solutions are not dominated and belong to $S_{i\{1,2,3\}}^2$. Finally, we have

$$S_{i\{1,2,3\}}^2 = \{(2|1, 3), (3|1, 2)\}$$

In conclusion, to deliver customers 1, 2, and 3 from $i$ using two drones:

* If the truck arrives at most at time 0, the optimal solution is to use a drone to deliver customer 3, and the other drone to deliver customer 1 and then customer 2. In that case, the makespan of the deliveries is 7.

* If the truck arrives at most at time 2, the optimal solution is to use a drone to deliver customer 2, and the other drone to deliver customer 1 and then customer 3. In that case, the makespan of the deliveries is 8.

* If the truck arrives after 2, there is no feasible solution.

A final remark is that if the truck travel time between the depot and customer $i$ is strictly greater than 0 ($c_{(0,i)}^t > 0$), then the drone schedule $(1, 2)$, in which a single drone serves customer 1 followed by customer 2, would not be included in $S_{iK}^1$, as the truck could not arrive before the schedule's maximal starting time ($\Delta_{(1,2)} = 0$). In that case, $((3), (1, 2))$ would not be computed in the set of possible combinations for $S_{i\{1,2,3\}}^2$, resulting in $(3|1, 2)$ not belonging in $S_{i\{1,2,3\}}^2$.

### 5.2.4. Criteria for valid arcs elimination

In this section, we introduce two criteria used in preprocessing to eliminate in-feasible and dominated arcs in $A^{dD}$, from a subproblem graph $G$, $d \in \{0, ..., D\}$.
and the SPR criterion are used. To measure the difficulty of an instance more precisely, we introduce $\hat{K} = \max_{i \in V} |K_i|$ the largest drone-reachable neighborhood of an instance, and $|\mathcal{A}^{\max}| = \max_{d \in D} |\mathcal{A}^d|$ the maximum number of arcs in a subproblem graph. $|\mathcal{A}^{\max}|$ is extremely dependent on $\hat{K}$, as the number of arcs in a subproblem graph is exponential in the size of the drone-reachable neighborhoods. Gap is calculated as follows: $\frac{\text{PB}-\text{DB}}{\text{DB}}$, with PB the primal bound and DP the dual bound.

We present in the first subsection a comparison with the state-of-the-art results on the small-size instances. In the second subsection, we present the performance of our algorithm on medium size instances. In a third section, we present the results of the heuristic branch-cut-and-price on instances with large drone-reachable neighborhoods.

### 7.3.1 Performance of the exact BCP on the small-size instances

Table 4 shows the summary of the performance of our algorithm on the small-size instances and a comparison with the results obtained in (Zhou et al., 2023). The comparison is fair, as we use a computer with similar speed. For each data set, we give the name of the set (Instances), the number of instances solved optimally over the total number of instances in the dataset (Solved), the average time in seconds required to solve the instances of the dataset (Time), the average number of nodes processed in the branch-cut-and-price algorithm (Nodes), and the maximum number of arcs in the subproblems graphs ($|\mathcal{A}^{\max}|$).

<table><caption>Table 4: Average exact performance comparison given small-size instances</caption><thead><tr><th rowspan="2">Instances</th><th colspan="2">Zhou et al.</th><th colspan="4">Our</th></tr><tr><th>Solved</th><th>Time (s)</th><th>Solved</th><th>Time (s)</th><th>Nodes</th><th>|A<sup>max</sup>|</th></tr></thead><tbody><tr><td>Cardiff_10_X</td><td>20/20</td><td>0.4</td><td>20/20</td><td>0.57</td><td>1</td><td>217</td></tr><tr><td>Cardiff_15_X</td><td>20/20</td><td>0.7</td><td>20/20</td><td>1.33</td><td>1</td><td>993</td></tr><tr><td>Cardiff_25_X</td><td>15/20</td><td>3561.5</td><td>20/20</td><td>86.25</td><td>2.5</td><td>5207</td></tr><tr><td>c10X_35</td><td>4/9</td><td>7408</td><td>9/9</td><td>556.30</td><td>8.12</td><td>2490</td></tr></tbody></table>

Our branch-price-and-cut algorithm can solve the small-size instances very efficiently, indeed it solves all instances featured in (Zhou et al., 2023), including 10 solved for the first time. The computation time is improved by a factor 13 for the largest instances of the set. The number of nodes processed is very low. Our algorithm solves any of the small-size instances in less than half an hour. Therefore we increase the difficulty of the instances, by adding customers.

### 7.3.2 Performance of the exact BCP on the medium-size instances

Table 5 shows the summary of the performance of our algorithm on the medium-size instances. For each data set, we give the name of the set (Instances), the
The first criterion is based on a maximum weight argument. Any schedule in $S_{iK}^d, d \in \{0, ..., D\}, i \in V, K \subseteq N_i$ does not correspond to a feasible drone schedule if the total payload is higher than the truck capacity, i.e. $q_i+d \times q_1^d + \sum_{k \in K} q_k > Q$ because the scheduling solution is not feasible in the routing problem, thus the arc representing this schedule is infeasible.

The second criterion is based on a replacement argument. Let a truck with $d$ drones be parked at customer $i$. Suppose there is a faster route (combining truck and drones) serving a set of customers $K$ than serving all of them by drone from $i$ (i.e., using an appropriate drone schedule in $S_{iK}^d$). In that case, we know that the corresponding solution in $S_{iK}^d$ will not be used in any optimal solution, and we can remove it safely. More formally, for a given $i \in V, d \in \{0, ..., D\}, K \in N_i$ and $j \in V \setminus K \cup \{i\}$, let SPR$(d, i, K, j)$ be the duration of the shortest partial route from $i$ to $j$ using a truck equipped with $d$ drones and serving either by drone or by truck the customers of set $K$. If for any feasible $j$, SPR$(d, i, K, j) < m_s + c_{i,j}^v$, $s$ will not be used in an optimal solution.
We consider the following set of feasible $j$:

$$ \left\{ j \in V_+ \setminus K \cup \{i\} \mid (c_{0,i}^v + \max(s_i^v, m_s) + c_{i,j}^v < \bar{d}_j) \wedge \left( q_i + d \times q_1^d + q_j + \sum_{k \in K} q_k \le Q \right) \right\} \quad (24) $$

<figure><img src="image_5.png" /><figcaption>Figure 4: Context for an application of SPR criterion with one drone.</figcaption></figure>

Figure 4 presents the intuition behind the SPR criterion. On the left side, a single drone serves customers $j,k,l$ from customer $i$. On the right side, a drone serves customer $j$ from customer $i$, then the truck moves to customer $k$, and then the drone serves customer $l$ from customer $k$. It may be possible that the total duration of the left-side solution is longer than that of the right-side
maximum number of arcs in the subproblems graphs (|A<sup>max</sup>|) of all instances, the average maximal size of drone-reachable neighborhood (K), the number of instances solved optimally over the total number instances (Solved), the number of instances for which no feasible solution has been found (No Sol), the average Gap of the instances for which a feasible solution has been found, the average time in seconds required to solve the instances of the dataset (Time) and the average number of nodes processed in the branch-cut-and-price algorithm (Nodes). The detailed results can be found in table 11, table 12 and table 13 in Appendix A.

<table><caption>Table 5: Average exact performance comparison given medium-size instances</caption><thead><tr><th>Instances</th><th>|A<sup>max</sup>|</th><th>K&#770;</th><th>Solved</th><th>No Sol</th><th>Gap</th><th>Time (s)</th><th>Nodes</th></tr></thead><tbody><tr><td>Cardiff_40_X</td><td>78422.55</td><td>13</td><td>9/20</td><td>6</td><td>3.12</td><td>6241.96</td><td>12.2</td></tr><tr><td>c10X_100</td><td>19009.23</td><td>11</td><td>8/9</td><td>0</td><td>0.14</td><td>3134.03</td><td>30.4</td></tr><tr><td>120_C1_2_X</td><td>14931.50</td><td>3</td><td>7/10</td><td>3</td><td>0</td><td>4430.94</td><td>58.8</td></tr><tr><td>130_C1_2_X</td><td>17615.80</td><td>4</td><td>5/10</td><td>4</td><td>0.07</td><td>7081.01</td><td>77.4</td></tr><tr><td>140_C1_2_X</td><td>20367.80</td><td>4</td><td>3/10</td><td>4</td><td>0.27</td><td>8229.01</td><td>84.3</td></tr><tr><td>150_C1_2_X</td><td>23332.80</td><td>4</td><td>2/10</td><td>4</td><td>0.45</td><td>9384.44</td><td>107.8</td></tr></tbody></table>

The size of the drone-reachable neighborhood has a heavy impact on the number of arcs of the subproblems graphs, which strongly affect calculation time. Therefore, the size of the drone-reachable neighborhood should be considered when measuring the size of an instance. Our algorithm is particularly efficient to solve clustered instances: we solve most clustered instances with 120 customers and some with up to 150 customers. It is worth mentioning that in Gehring and Homberger (1999) instances, due to the weight of the packages, less than 50% of the customers are available for drone delivery, leading to very small drone-reachable neighborhoods. Increasing the number of customers in Gehring and Homberger (1999) instances slightly increases the gap, the number of processed nodes, and decreases the number of instances solved to optimality. However with randomized instances, the performance of our method deteriorates much faster. Already with 40 customers the algorithm can only find optimal solutions for half of the Chen et al. (2021) instances and cannot find any feasible solution for 6 instances. This can be explained by the large and highly interconnected drone-reachable neighborhoods in the randomized instances. Note that instances *Cardiff_40_X* in addition of having a high number of arcs in average, also have great disparity among the instances of the set, as the standard deviation of the maximum number of arc in a subproblem graph is 61705.8. In *Cardiff_40_X* instances, we only solve less than half instances. Due to the extremely high number of arcs, solving a node is slow. To obtain good quality solutions, we propose to use the heuristic defined in section 6.
solution. In that case, if any solution using the left-side deliveries can replace it with the right-side deliveries without diminishing the objective value while remaining feasible, the left-side deliveries do not have to be considered, and we can remove $S_{i\{jkl\}}^1$.

# 6. Heuristic branch-cut-and-price

In this section, we propose a heuristic branch-cut-and-price based on reducing the size of the subproblems graphs. It has four parameters :

$D_h$ — The maximum number of drones a truck can carry.

$w_h$ — The maximum number of customers that can be served by drone from a customer in a schedule.

$k_h$ — The maximum number of drone-reachable neighbors. Only the $k_h$ closest drone-reachable neighbors are kept.

$p_h$ — The number of schedule with earliest deadline kept for each $S_{iK}^d$.

$D_h$ decreases the number of subproblems. $w_h$ and $k_h$ reduce the number of drone arcs in each subproblem graph. Note that the difference between $w_h$ and $k_h$ is that $k_h$ controls which customer can be served by drone from a parked vehicle, while $w_h$ controls how many customers can be served by drone from a parked vehicle. Finally, $p_h$ reduces the number of schedules and, therefore, the number of arcs in subproblems with multiple drones.

Let $\mathcal{N}_i^k$ be the set of the $k$ closest to $i$ customers in $\mathcal{N}_i$. More formally, $\mathcal{N}_i^k \subseteq \mathcal{N}_i$, $|\mathcal{N}_i^k| = k$ and $\forall k \in \mathcal{N}_i^k, k' \in \mathcal{N}_i \setminus \mathcal{N}_i^k, c_{ik}^d \le c_{ik'}^d$. Let $S_{iK}^{dp}$ be the set of the $p_h$ solutions from $S_{iK}^d$ with earliest deadlines. To apply these rules, we transform each subproblems graphs $\mathcal{G}^d = (\mathcal{V}^d, \mathcal{A}^d), d \in \{1, ..., D\}$ with $\mathcal{A}^d = A^T \cup A^D$ from the exact branch-cut-and-price subproblems graphs as follows: we create the graphs $\mathcal{G}_h^d = (\mathcal{V}_h^d, \mathcal{A}_h^d), d \in \{1, ..., D_h\}$, where $\mathcal{V}_h^d = \mathcal{V}^d$ and $\mathcal{A}_h^d = A^T \cup \{(i, i')^s \in A^D, \forall K \subset \mathcal{N}_i^{k_h}, s \in S_{iK}^{dp_h} : |K| \le w_h\}$

# 7. Computational experiments

## 7.1. Configuration

Our algorithms, such as the dynamic program for the enumeration, the separation algorithm for the ARCC, and the SPR criterion, are coded in C++. We also used:

* IBM ILOG CPLEX 20.1 as the LP solver in the relaxed restricted master problem and as the solver for the enumerated MIPs.
### 7.3.3 Performance of the heuristic BCP on the random medium-size instances

To evaluate the performance of the proposed heuristic, we tested several different settings, which are reported in Table 6. The last row of the table corresponds to the performance of the exact BCP algorithm. In the absence of results in the literature, we compare all the results obtained by both the exact BCP and the heuristic BCP, and retain the BKS from one of these methods. The average gap between the heuristic solution and the BKS is then calculated as follows: $\frac{\text{PB-BKS}}{\text{BKS}}$, with PB the primal bound and BKS the best known solution value (Gap BKS).

Table 6 shows the summary of the performance of our heuristic branch-cut-and-price on the Chen et al. (2021) instances, with different sets of parameters. For each set of parameters $\Gamma_h$, $w_h$, $k_h$ and $p_h$ we give the average number of arcs in the largest subproblem of each instance ($|\mathcal{A}^{\max}|$), the average number of nodes processed in the branch-cut-and-price algorithm (Nodes), the average gap between the best known solution and the solution found in percentage for the instances where a solution was found (Gap BKS), the number of instances for which the solver reached the time limit (Time lim), the number of instances for which the solver was not able to find a feasible solution (No sol) and the average time solving time in seconds (Time). When the value of $p_h$ is inf, it means the parameter is left unrestricted. The detailed results can be found in table 14 in Appendix A.

Note that the maximum authorized number of drones per truck in these instances is 3, and the largest drone-neighborhood is 15, hence setting the related parameters to these values is similar to leaving them unrestricted.

<table><caption>Table 6: Heuristic average performances given parameters on the Chen et al. (2021) instances with 40 customers</caption><thead><tr><th>Γ<sub>h</sub></th><th>w<sub>h</sub></th><th>k<sub>h</sub></th><th>p<sub>h</sub></th><th>|A<sup>max</sup>|</th><th>Time lim</th><th>No sol</th><th>Gap BKS</th><th>Time (s)</th><th>Nodes</th></tr></thead><tbody><tr><td>2</td><td>2 2</td><td>1</td><td>1796.30</td><td>0</td><td>0</td><td>4</td><td>392.72</td><td>5</td><td></td></tr><tr><td>2</td><td>3 3</td><td>1</td><td>1937.20</td><td>0</td><td>0</td><td>2.09</td><td>510.60</td><td>7.1</td><td></td></tr><tr><td>2</td><td>5 5</td><td>1</td><td>2641.55</td><td>0</td><td>0</td><td>1.05</td><td>360.44</td><td>5.6</td><td></td></tr><tr><td>2</td><td>8 5</td><td>1</td><td>6250.15</td><td>6</td><td>0</td><td>0.74</td><td>3738.86</td><td>19.4</td><td></td></tr><tr><td>3</td><td>5 3</td><td>1</td><td>2490.50</td><td>0</td><td>0</td><td>0.85</td><td>681.09</td><td>9.7</td><td></td></tr><tr><td>3</td><td>5 5</td><td>1</td><td>2664.95</td><td>0</td><td>0</td><td>0.67</td><td>381.38</td><td>5.4</td><td></td></tr><tr><td>3</td><td>5 5</td><td>inf</td><td>2773.70</td><td>0</td><td>0</td><td>0.67</td><td>348.23</td><td>5.1</td><td></td></tr><tr><td>3</td><td>8 5</td><td>1</td><td>6465.15</td><td>4</td><td>1</td><td>0.31</td><td>3496.81</td><td>25.2</td><td></td></tr><tr><td>3</td><td>8 5</td><td>inf</td><td>7767.95</td><td>6</td><td>2</td><td>0.26</td><td>3878.18</td><td>20.6</td><td></td></tr><tr><td>3</td><td>8 8</td><td>1</td><td>7088.75</td><td>6</td><td>2</td><td>0.18</td><td>4052.90</td><td>23.2</td><td></td></tr><tr><td>3</td><td>15 10</td><td>1</td><td>15222.95</td><td>5</td><td>0</td><td>0.04</td><td>4125.97</td><td>29.3</td><td></td></tr><tr><td>3</td><td>15 10</td><td>inf</td><td>78192.80</td><td>11</td><td>8</td><td>2.49</td><td>6203.80</td><td>12.5</td><td></td></tr><tr><td>3</td><td>15 15</td><td>1</td><td>43094.35</td><td>9</td><td>7</td><td>0.44</td><td>5625.45</td><td>13.5</td><td></td></tr><tr><td>3</td><td>15 15</td><td>inf</td><td>78422.55</td><td>11</td><td>6</td><td>1.89</td><td>6241.96</td><td>13</td><td></td></tr></tbody></table>
* The generic branch-cut-and-price code from (Sadykov and Vanderbeck, 2021), with automatic dual price smoothing stabilization from (Pessoa et al., 2018), primal heuristics from (Sadykov et al., 2019). We also use the labeling algorithm code from (Sadykov et al., 2021).

The experiments presented in this paper were carried out:

* Using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux. Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr).

* Running on 2.6 GHz 2x 18-core Cascade Lake Intel Xeon Skylake Gold 6240 CPU using a single thread with 5.3 Go memory.

## 7.2. Instances

To measure the performance of our algorithm, we used and extended the sets of instances from (Zhou et al., 2023), which is composed of modified instances from (Chen et al., 2021), and from (Solomon, 1987). In all these instances, the trucks and the drones follow the same distance metric. We also used clustered instances from (Gehring and Homberger, 1999) to experiment on larger instances, as they are an extension of the Solomon (1987) instances. The maximal drone capacity is set to 10 kg as in (Chen et al., 2021), the drone service times are set to $s_i^d = [0.5s_i^v]$, $i \in V$.

In the instances proposed in (Chen et al., 2021), an additional $[0.2 \times |V|]$ customers are randomly selected to be unavailable for drone delivery. Table 4 describes the different data sets. For each data set, we give the name of the instances (Instances) where $X$ is the id of a particular instance of the set, the number of instances in the set (Size), the customer positioning (Pos), the distance metric (Dist), the original authors of the instances (Authors), the capacity of the vehicle ($Q$), the maximum number of drones allowed on a vehicle ($D$), the total number of vehicles ($k_v$) and the total number of drones ($k_d$). In column (Pos), "Rand" stands for "Random", and "Clust" for "Clustered". In column (Dist), "Eucl" stands for "Euclidean."
We can observe in table 6 that reducing the number of arcs in $\mathcal{A}^{dD}$ with $d \in [1, ..., \Gamma_h]$ by decreasing $w_h$ or $k_h$ can significantly decrease computation time, without decreasing the solutions quality much. For these instances, reducing $\Gamma_h$ has a heavy negative impact on the gap BKS, without having a positive impact on the computation time. It is worth mentioning that the heuristic with parameters $\Gamma_h = 3$, $w_h = 15$, $k_h = 10$, $p_h = 1$ gives a major improvement compared to the exact method, as both the average gap and the average computation time largely decrease. The heuristic with parameters $\Gamma_h = 3$, $w_h = 5$, $k_h = 5$, $p_h = 1$ despite having a slightly larger gap has a very fast average solving time of less than 10 minutes. Unrestricted parameter $p_h$ has a major negative impact on the number of arcs and on the computation time, while having no positive impact on the gap. That is why we set it to one in most cases.

## 7.4 Sensitivity analysis

In this section we first present the impact of the ARCC and then we present the impact of the SPR criterion. We use the concept of cut off values. Cut off values are primal bound that the user provides to the algorithm, in order to prune the branch and bound tree. We provide the best known solution as cut off value. It can be use to prove the optimality of a solution, but here the goal is to avoid random performance of the primal heuristics, so we can fairly compare the performance with and without the tested improvement. Unless specified otherwise, the ARCC and the SPR criterion are used. We recall that the gap is calculated as follows: $\frac{PB-DB}{DB}$, with PB the primal bound and DP the dual bound.

### 7.4.1 Impact of the adapted rounded capacity cuts

In this section, we present the impact of the adapted rounded capacity cuts. The capacity constraint in the Chen et al. (2021) instances is not tight, as a large proportion of the trucks capacity is not used in optimal solution, therefore the ARCC are not useful on these instances, and the instances with 100 customers are all solved within 10 minutes, most of them solved at root node. For these reasons, we focus on the instances from Gehring and Homberger (1999) in this section. As shown in table 5, within the three-hour limit, the exact branch-cut-and-price does not find feasible solution for several instances of the Gehring and Homberger (1999) datasets. To find accurate cut off values, we run the exact BCP on these instances with a ten-hour time limit. However, we did not find feasible solution for all instances with 150 customers. Therefore, we present the results on the Gehring and Homberger (1999) instances with 120, 130 and 140 customers in table 7. We compare the performances with and without the ARCC. Table 7 shows the summary of the experiments to measure the impact of the ARCC. For each data set, we give the name of the set (Instances), whether the ARCC were used or not (ARCC), the average gap obtained on instances that are not solved optimally within the time limit (Gap), the average time in
<table><caption>Table 4: Instances sets and their characteristics</caption><thead><tr><th>Instances</th><th>|V|</th><th>Size</th><th>Pos</th><th>Dist</th><th>Authors</th><th>Q</th><th>D</th><th>k<sub>v</sub></th><th>k<sub>d</sub></th></tr></thead><tbody><tr><td colspan="10" style="text-align: center; font-weight: bold;">Small-Size Instances</td></tr><tr><td>Cardiff_10_X</td><td>10</td><td>20</td><td>Rand</td><td>Real</td><td>Chen et al. (2021)*</td><td>170</td><td>2</td><td>4</td><td>4</td></tr><tr><td>Cardiff_15_X</td><td>15</td><td>20</td><td>Rand</td><td>Real</td><td>Chen et al. (2021)</td><td>170</td><td>2</td><td>5</td><td>5</td></tr><tr><td>Cardiff_25_X</td><td>25</td><td>20</td><td>Rand</td><td>Real</td><td>Chen et al. (2021)</td><td>200</td><td>3</td><td>5</td><td>10</td></tr><tr><td>c10X_35</td><td>35</td><td>9</td><td>Clust</td><td>Eucl</td><td>Solomon (1987)</td><td>250</td><td>3</td><td>5</td><td>10</td></tr><tr><td colspan="10" style="text-align: center; font-weight: bold;">Medium-Size Instances</td></tr><tr><td>Cardiff_40_X</td><td>40</td><td>20</td><td>Rand</td><td>Real</td><td>Chen et al. (2021)*</td><td>200</td><td>3</td><td>5</td><td>10</td></tr><tr><td>c10X_100</td><td>100</td><td>9</td><td>Clust</td><td>Eucl</td><td>Solomon (1987)</td><td>250</td><td>3</td><td>14</td><td>28</td></tr><tr><td>120_C1_2_X</td><td>120</td><td>10</td><td>Clust</td><td>Eucl</td><td>Gehring and Homberger (1999)</td><td>250</td><td>3</td><td>40</td><td>60</td></tr><tr><td>130_C1_2_X</td><td>130</td><td>10</td><td>Clust</td><td>Eucl</td><td>Gehring and Homberger (1999)</td><td>250</td><td>3</td><td>40</td><td>60</td></tr><tr><td>140_C1_2_X</td><td>140</td><td>10</td><td>Clust</td><td>Eucl</td><td>Gehring and Homberger (1999)</td><td>250</td><td>3</td><td>40</td><td>60</td></tr><tr><td>150_C1_2_X</td><td>150</td><td>10</td><td>Clust</td><td>Eucl</td><td>Gehring and Homberger (1999)</td><td>250</td><td>3</td><td>40</td><td>60</td></tr></tbody></table>

Note that the set with 10 customers and the set with 40 customers are marked with '*' as they are not in (Chen et al., 2021), but are derived from sets with a larger number of customers from (Chen et al., 2021), namely 15 and 50, by keeping only the 10 and 40 first customers.

The instances `c10X_N` are the 100-customers clustered instances called `c10X` in (Solomon, 1987), restricted to the *N* first customers. We categorize instances with 35 or fewer customers as small, while instances with more than 35 customers will be referred to as medium. The instances `N_C1_2_X` are the 200-customer clustered instances called `C1_2_X` in (Gehring and Homberger, 1999), restricted to the first *N* customers.

### 7.3. Results

This section presents all the results obtained by setting the computation time limit to three hours per instance. Unless specified otherwise, the ARCC and the SPR criterion are used. To measure the difficulty of an instance more precisely, we introduce $\hat{\mathcal{K}} = \max_{i \in V} |\mathcal{N}_i|$ the largest drone-reachable neighborhood of an instance, and $|\mathcal{A}^{\max}| = \max_{d \in D} |\mathcal{A}^d|$ the maximum number of arcs in a subproblem graph. $|\mathcal{A}^{\max}|$ is strongly dependent on $\hat{\mathcal{K}}$, as the number of arcs in a subproblem graph is exponential in the size of the drone-reachable neighborhoods. The gap is calculated as follows: $\frac{\text{PB}-\text{DB}}{\text{DB}}$, with PB the primal bound and DP the dual bound.

In the first subsection, we present a comparison with the state-of-the-art results on the small-size instances. In the second subsection, we present the performance of our algorithm on medium-sized instances. In the third subsection, we present the results of the heuristic branch-cut-and-price on instances with large drone-reachable neighborhoods.
seconds required to solve the instances of the dataset (Time) and the average number of nodes processed in the branch-cut-and-price algorithm (Nodes). Note that the number of optimal solutions found is omitted as it is the same for both. The detailed results can be found in the table 15 in Appendix A.

<table><caption>Table 7: Average impact of the adapted rounded capacity cuts on the (Gehring and Homberger, 1999) instances</caption><thead><tr><th>Instances</th><th>ARCC</th><th>Gap</th><th>Time (s)</th><th>Nodes</th></tr></thead><tbody><tr><td>120_C1_2_X</td><td>yes</td><td>0.12</td><td>2437.30</td><td>31.6</td></tr><tr><td>120_C1_2_X</td><td>no</td><td>0.12</td><td>2451.13</td><td>40.6</td></tr><tr><td>130_C1_2_X</td><td>yes</td><td>0.28</td><td>5784.12</td><td>88.2</td></tr><tr><td>130_C1_2_X</td><td>no</td><td>0.30</td><td>5744.47</td><td>67.6</td></tr><tr><td>140_C1_2_X</td><td>yes</td><td>0.56</td><td>5955.80</td><td>75.4</td></tr><tr><td>140_C1_2_X</td><td>no</td><td>0.61</td><td>5925.10</td><td>72.6</td></tr></tbody></table>

We can observe on table 7 that the ARCC have a positive impact of the gap, growing with the size of the instances. The ARCC have a negligible impact on the computation time, this is due to the fact that the capacity constraints are not very tight for these instances.

### 7.4.2 Impact of the shorter partial route criterion

In this section we study the impact of the shorter partial route criterion both in terms of efficiency of the exact algorithm, and in terms of graph size. As the instances from (Gehring and Homberger, 1999) have a small drone-reachable neighborhood, this criterion does not have a significant impact on the number of arcs. For this reason, we perform the sensitivity analysis only for the medium-size instances from (Chen et al., 2021) and (Solomon, 1987). We compare the performances of the algorithm with and without the SPR criterion. We recall that $G^d$ is the directed graph representing the delivery possibilities for a truck with $d$ drones. Table 8 shows the summary of the experiments to measure the impact of the SPR criterion on the medium-size instances. For each data set, we give the name of the set (Instances), the number of arcs in subproblem graph and the proportion of arcs compared to without the SPR criterion (|$A^1$|), (|$A^2$|), (|$A^3$|), the average Gap of the instances for which a feasible solution has been found, the average time in seconds required to solve the instances of the dataset (Time) and the average number of nodes processed in the branch-cut-and-price algorithm (Nodes).
### 7.3.1. Performance of the exact BCP on the small-size instances

Table 5 shows the summary of the performance of our algorithm on the small-size instances and a comparison with the results obtained in (Zhou et al., 2023). The comparison is fair, as we use a computer with similar speed. For each data set, we give the name of the set (Instances), the number of instances solved optimally over the total number instances in the dataset (Solved), the average time in seconds required to solve the instances of the dataset (Time), the average number of nodes processed in the branch-cut-and-price algorithm (Nodes), and the maximum number of arcs in the subproblems graphs ($|\mathcal{A}^{\max}|$).

<table><caption>Table 5: Average exact performance comparison given small-size instances</caption><thead><tr><th rowspan="2">Instances</th><th colspan="2">Zhou et al.</th><th colspan="4">Our</th></tr><tr><th>Solved</th><th>Time (s)</th><th>Solved</th><th>Time (s)</th><th>Nodes</th><th>|A<sup>max</sup>|</th></tr></thead><tbody><tr><td>Cardiff_10_X</td><td>20/20</td><td>0.4</td><td>20/20</td><td>0.57</td><td>1</td><td>217</td></tr><tr><td>Cardiff_15_X</td><td>20/20</td><td>0.7</td><td>20/20</td><td>1.33</td><td>1</td><td>993</td></tr><tr><td>Cardiff_25_X</td><td>15/20</td><td>3561.5</td><td>20/20</td><td>86.25</td><td>2.5</td><td>5207</td></tr><tr><td>c10X_35</td><td>4/9</td><td>7408</td><td>9/9</td><td>556.30</td><td>8.12</td><td>2490</td></tr></tbody></table>

Our branch-price-and-cut algorithm can solve the small-size instances very efficiently. Indeed, it solves all instances featured in (Zhou et al., 2023), including 10 solved for the first time. The computation time is improved by a factor of 13 for the largest instances of the set. The number of nodes processed is very low. Our algorithm solves any small instance in less than half an hour. Therefore, we increase the difficulty of the instances by adding customers.

### 7.3.2. Performance of the exact BCP on the medium-size instances

Table 6 shows the summary of the performance of our algorithm on the medium-size instances. For each data set, we give the name of the set (Instances), the maximum number of arcs in the subproblems graphs ($|\mathcal{A}^{\max}|$) of all instances, the average maximal size of the drone-reachable neighborhoods ($\hat{\mathcal{K}}$), the number of instances solved optimally over the total number instances (Solved), the number of instances for which no feasible solution has been found (No Sol), the average Gap of the instances for which a feasible solution has been found, the average time in seconds required to solve the instances of the dataset (Time) and the average number of nodes processed in the branch-cut-and-price algorithm (Nodes). The detailed results can be found in Table 11, Table12, Table 13, and Table 14 in Appendix B.
<table><caption>Table 8: Sensitivity analysis over SPR criterion</caption><thead><tr><th>Instances</th><th>SPR</th><th>|A<sup>1</sup>|</th><th>|A<sup>2</sup>|</th><th>|A<sup>3</sup>|</th><th>Gap</th><th>Time (s)</th><th>Nodes</th></tr></thead><tbody><tr><td>Cardiff-40-X</td><td>yes</td><td>4449.15 (20%)</td><td>62169.85 (54%)</td><td>72709.35 (97%)</td><td>6.84</td><td>6376.47</td><td>10.7</td></tr><tr><td>Cardiff-40-X</td><td>no</td><td>22597.45 (100%)</td><td>116217.80 (100%)</td><td>75056.20 (100%)</td><td>7.26</td><td>6318.14</td><td>8.8</td></tr><tr><td>c10X_100</td><td>yes</td><td>12902 (68%)</td><td>18298.67 (95%)</td><td>19009.23 (100%)</td><td>0</td><td>206.81</td><td>5.23</td></tr><tr><td>c10X_100</td><td>no</td><td>19079 (100%)</td><td>19326.56 (100%)</td><td>19009.23 (100%)</td><td>0</td><td>274.63</td><td>5.89</td></tr></tbody></table>

Note that using cut off values, the algorithm solves all *c10X_100* instances to optimality. The criterion removes a very large number of arcs in *G*¹ in both cases, less in *G*² and almost none in *G*³. It removes more arcs when instances have larger drone-reachable neighborhood. For both data sets, the criterion accelerates most solved instances. The unsolved instances gaps are decreased of 0.5% in average. Note that although the criterion seems to have a negative impact on computation time in the *Cardiff-40-X* instances, in fact, the time decreases for 7 out of 9 instances solved optimally. Also the criterion decreases the number of nodes in the *c10X_100* instances. It is a positive impact, as it improves convergence speed. The criterion increases the number of nodes in instances *Cardiff-40-X*, which is also a positive result, as it is mainly due to the number of nodes in instances where the optimal solution was not found, demonstrating a faster node processing.

### 7.4.3 Impact of the restricted master heuristic

In this section, we present the impact of the restricted master heuristic. In this experiment, cut off values are not provided to the solver.

<table><caption>Table 9: Sensitivity analysis over restricted master heuristic.</caption><thead><tr><th>Instances</th><th>RMH</th><th>Solved</th><th>No Sol</th><th>Gap</th><th>Time</th><th>Nodes</th></tr></thead><tbody><tr><td>Cardiff-40-X</td><td>yes</td><td>9/20</td><td>6</td><td>3.12</td><td>6241.96</td><td>12.2</td></tr><tr><td>Cardiff-40-X</td><td>no</td><td>11/20</td><td>2</td><td>3.64</td><td>5621.62</td><td>39.0</td></tr><tr><td>c10X_100</td><td>yes</td><td>8/9</td><td>0</td><td>0.14</td><td>3134.03</td><td>30.4</td></tr><tr><td>c10X_100</td><td>no</td><td>8/9</td><td>1</td><td>0.00</td><td>4184.41</td><td>57.0</td></tr><tr><td>120_C1_2_X</td><td>yes</td><td>7/10</td><td>3</td><td>0</td><td>4430.94</td><td>58.8</td></tr><tr><td>120_C1_2_X</td><td>no</td><td>6/10</td><td>2</td><td>0.19</td><td>6116.54</td><td>132.40</td></tr></tbody></table>

We can observe that the number of nodes increases drastically when the restricted master heuristics are deactivated.

## 8 Conclusion

In this paper, we introduced a new branch-cut-and-price algorithm for the two-echelon vehicle routing problem with drones, based on the enumeration of the
<table><caption>Table 6: Average exact performance comparison given medium-size instances</caption><thead><tr><th>Instances</th><th>|A<sup>max</sup>|</th><th>Ĉ</th><th>Solved</th><th>No</th><th>Sol</th><th>Gap</th><th>Time (s)</th><th>Nodes</th></tr></thead><tbody><tr><td>Cardiff_40_X</td><td>78422.55</td><td>13</td><td>9/20</td><td>6</td><td>3.12</td><td>6241.96</td><td>12.2</td><td></td></tr><tr><td>c10X_100</td><td>19009.23</td><td>11</td><td>8/9</td><td>0</td><td>0.14</td><td>3134.03</td><td>30.4</td><td></td></tr><tr><td>120_C1_2_X</td><td>14931.50</td><td>3</td><td>7/10</td><td>3</td><td>0</td><td>4430.94</td><td>58.8</td><td></td></tr><tr><td>130_C1_2_X</td><td>17615.80</td><td>4</td><td>5/10</td><td>4</td><td>0.07</td><td>7081.01</td><td>77.4</td><td></td></tr><tr><td>140_C1_2_X</td><td>20367.80</td><td>4</td><td>3/10</td><td>4</td><td>0.27</td><td>8229.01</td><td>84.3</td><td></td></tr><tr><td>150_C1_2_X</td><td>23332.80</td><td>4</td><td>2/10</td><td>4</td><td>0.45</td><td>9384.44</td><td>107.8</td><td></td></tr></tbody></table>

The size of the drone-reachable neighborhood significantly impacts the number of arcs of the subproblems graphs, which strongly affects calculation time. Therefore, the size of the drone-reachable neighborhood should be considered when measuring the size of an instance. Our algorithm is particularly efficient for solving clustered instances: it solves most clustered instances with 120 customers and some with up to 150 customers. It is worth mentioning that in Gehring and Homberger (1999) instances, due to the weight of the packages, less than 50% of customers are available for drone delivery, resulting in limited size drone-reachable neighborhoods. Increasing the number of customers in Gehring and Homberger (1999) instances slightly increases the gap and the number of processed nodes and decreases the number of instances solved to optimality. However, with randomized instances, the performance of our method deteriorates much faster. With 40 customers, the algorithm can only find optimal solutions for half of the Chen et al. (2021) instances and cannot find any feasible solution for 6 instances. The modest performance can be attributed to the large and highly interconnected drone-reachable neighborhoods in the randomized instances, leading to an important number of set $S_{iK}^d$. However, even in these instances, the number of solutions in each $S_{iK}^d$ remained tractable. Note that instances Cardiff_40_X, in addition to having a high number of arcs on average, also have significant disparity among the instances of the set, as the standard deviation of the maximum number of arcs in a subproblem graph is 61705.8. We only solve less than half instances of Cardiff_40_X dataset. Due to the extremely high number of arcs, solving a node is slow. We propose to use the heuristic defined in section 6 to obtain good quality solutions.

### 7.3.3. Performance of the heuristic BCP on the random medium-size instances

To evaluate the performance of the proposed heuristic, we tested several different settings, which are reported in Table 7. The last row of the table corresponds to the performance of the exact BCP algorithm. In the absence of results in the literature, we compare all the results obtained by both the exact BCP and the heuristic BCP and retain the BKS from one of these methods. The average gap between the heuristic solution and the BKS is then calculated as follows: $\frac{\text{PB-BKS}}{\text{BKS}}$, with PB the primal bound and BKS the best-known solution value (Gap BKS).

Table 7 shows the summary of the performance of our heuristic branch-cut-
drones round trips between customers and the vehicle parked at another customer. We presented ways to improve the method by reducing the number of drone trips enumerated, and by adapting the well-known rounded capacity cuts, supported by sensitivity analysis. We also proposed an efficient heuristic to tackle dense instances.

Extensive numerical experiments demonstrated the high efficiency of our method. It allows for solving all instances of literature containing up to 35 customers improving significantly computation times by a factor higher than 10 and optimally solving 10 of them for the first time. We also proposed 6 new sets of instances to promote research on VRP-D type problems. The experiments on our method shows that we solve all but one clustered instance with 100 customers, and more than half of the clustered instances with 120 customers. In the random dataset with 40 customers and large drone-reachable neighborhood, 11 instances remain open. To tackle these instances we proposed an heuristic branch-cut-and-price based on the elimination of certain arcs in the subproblems, and using the exact approach for solving. This approach shows great results, as the gaps with the best known solutions are very small. We also show the improvement caused by our adaptation of the rounded capacity cuts, and by the SPR criterion.

Testing these new instances highlighted some limitations of our branch-cut-and-price approach. On the largest instances, the size of the sub-problem graph becomes untractable due to the number of parallel arcs representing the drones round trips possibilities which is exponential given the number of customers that can be delivered by drone from a customer. A possible remedy is to model the scheduling problem as a part of the RCSP, to avoid the explicit enumeration of the scheduling solutions. This would require a special care for the resources, as only one drone is responsible for the duration of the makespan. In the same vein, it would be useful to find other rules of dominance for the drones arcs. In order to increase the efficiency of our algorithm, another possible enhancement concerns the primal heuristic to quickly find good quality primal bounds. It would also be very valuable to create instances more suitable for truck/drone deliveries, with different metrics for the drones and for the trucks and realistic features regarding last miles delivery service. One could also consider uncertainty for the travel with the truck. Comparing the gain using the 2E-VRP-D and the VRPD, on different sets of instances and objective functions could also be very interesting. In addition, most public authorities at national or metropolitan level are currently moving towards tighter regulation of deliveries in city centers. In particular, more restrictions are being imposed on access times for vehicles, and therefore also for customer deliveries. We could integrate time window constraints with non-zero ready times into our model. In particular, this extension will require a special care for the drones operations, as drones might need to wait between round trips. Another perspective is to apply this type of enumeration strategy to problems with a similar structure, such as the park-and-loop vehicle routing problem for instance, in which the driver of a vehicle parked at a client can make a walking tour to deliver nearby customers, to avoid network congestion and use of gas.
and-price on the Chen et al. (2021) instances with different sets of parameters. For each set of parameters $D_h$, $w_h$, $k_h$, and $p_h$, we give the average number of arcs in the largest subproblem of each instance ($|\mathcal{A}^{\max}|$), the average number of nodes processed in the branch-cut-and-price algorithm (Nodes), the average gap between the best-known solution and the solution found in percentage for the instances where a solution was found (Gap BKS), the number of instances for which the solver reached the time limit (Time lim), the number of instances for which the solver was not able to find a feasible solution (No sol) and the average solving time in seconds (Time). When the value of $p_h$ is inf, it means the parameter is left unrestricted. The detailed results can be found in table 15 in Appendix B.

Note that the maximum authorized number of drones per truck in these instances is 3, and the largest drone-reachable neighborhood is 15; hence, setting the related parameters to these values is similar to leaving them unrestricted.

<table><caption>Table 7: Heuristic average performances given parameters on the Chen et al. (2021) instances with 40 customers</caption><thead><tr><th>D<sub>h</sub></th><th>w<sub>h</sub></th><th>k<sub>h</sub></th><th>p<sub>h</sub></th><th>|A<sup>max</sup>|</th><th>Time lim</th><th>No sol</th><th>Gap BKS</th><th>Time (s)</th><th>Nodes</th></tr></thead><tbody><tr><td>2</td><td>2 2</td><td>1</td><td>1796.30</td><td>0</td><td>0</td><td>4</td><td>392.72</td><td>5</td><td></td></tr><tr><td>2</td><td>3 3</td><td>1</td><td>1937.20</td><td>0</td><td>0</td><td>2.09</td><td>510.60</td><td>7.1</td><td></td></tr><tr><td>2</td><td>5 5</td><td>1</td><td>2641.55</td><td>0</td><td>0</td><td>1.05</td><td>360.44</td><td>5.6</td><td></td></tr><tr><td>2</td><td>8 5</td><td>1</td><td>6250.15</td><td>6</td><td>0</td><td>0.74</td><td>3738.86</td><td>19.4</td><td></td></tr><tr><td>3</td><td>5 3</td><td>1</td><td>2490.50</td><td>0</td><td>0</td><td>0.85</td><td>681.09</td><td>9.7</td><td></td></tr><tr><td>3</td><td>5 5</td><td>1</td><td>2664.95</td><td>0</td><td>0</td><td>0.67</td><td>381.38</td><td>5.4</td><td></td></tr><tr><td>3</td><td>5 5</td><td>inf</td><td>2773.70</td><td>0</td><td>0</td><td>0.67</td><td>348.23</td><td>5.1</td><td></td></tr><tr><td>3</td><td>8 5</td><td>1</td><td>6465.15</td><td>4</td><td>1</td><td>0.31</td><td>3496.81</td><td>25.2</td><td></td></tr><tr><td>3</td><td>8 5</td><td>inf</td><td>7767.95</td><td>6</td><td>2</td><td>0.26</td><td>3878.18</td><td>20.6</td><td></td></tr><tr><td>3</td><td>8 8</td><td>1</td><td>7088.75</td><td>6</td><td>2</td><td>0.18</td><td>4052.90</td><td>23.2</td><td></td></tr><tr><td>3</td><td>15 10</td><td>1</td><td>15222.95</td><td>5</td><td>0</td><td>0.04</td><td>4125.97</td><td>29.3</td><td></td></tr><tr><td>3</td><td>15 10</td><td>inf</td><td>78192.80</td><td>11</td><td>8</td><td>2.49</td><td>6203.80</td><td>12.5</td><td></td></tr><tr><td>3</td><td>15 15</td><td>1</td><td>43094.35</td><td>9</td><td>7</td><td>0.44</td><td>5625.45</td><td>13.5</td><td></td></tr><tr><td>3</td><td>15 15</td><td>inf</td><td>78422.55</td><td>11</td><td>6</td><td>1.89</td><td>6241.96</td><td>13</td><td></td></tr></tbody></table>

We can observe in table 7 that reducing the number of arcs in $\mathcal{A}^{DD}$ with $d \in \{1, ..., D_h\}$ by decreasing $w_h$ or $k_h$ can significantly decrease computation time, without substantially decreasing the quality of the solutions. Reducing $D_h$ has a significant negative impact on the gap BKS without decreasing the computation time for these instances. It is worth mentioning that the heuristic with parameters $D_h = 3$, $w_h = 15$, $k_h = 10$, and $p_h = 1$ gives a major improvement compared to the exact method, as both the average gap and the average computation time largely decrease. The heuristic with parameters $D_h = 3$, $w_h = 5$, $k_h = 5$, $p_h = 1$, despite having a slightly larger gap, exhibits a swift average solving time of less than 10 minutes. Unrestricted parameter $p_h$ significantly increases the number of arcs and the computation time, while having no positive
# Appendix A Detailed computational experiments

In the next tables we present the detailed results of our experiments. We present different information among: the name of the set (Instances), the maximum number of arcs in the subproblem $G^\Gamma$ ($|A^\Gamma|$) of all instances, the average maximal size of drone-reachable neighborhood ($\hat{K}$), the number of instances solved optimally over the total number instances (Solved), the number of instances for which no feasible solution have been found (No Sol), the average Gap of the instances for which a feasible solution have been found, the average time in seconds required to solve the instances of the dataset (Time) and the average number of nodes processed in the branch-cut-and-price algorithm (Nodes).

<table><caption>Table 10: Detailed performance of the exact algorithm on the small-size instances</caption><thead><tr><th>Instance</th><th>|A<sup>max</sup>|</th><th>Root Gap</th><th>Root Time</th><th>Nodes</th><th>DB</th><th>PB</th><th>Gap</th><th>Time (s)</th></tr></thead><tbody><tr><td>Cardiff_10_01</td><td>217</td><td>0</td><td>0.72</td><td>1</td><td>1223.62</td><td>1223.62</td><td>0</td><td>0.79</td></tr><tr><td>Cardiff_10_02</td><td>134</td><td>0</td><td>0.447</td><td>1</td><td>1905.08</td><td>1905.08</td><td>0</td><td>0.49</td></tr><tr><td>Cardiff_10_03</td><td>147</td><td>0</td><td>0.543</td><td>1</td><td>1412.67</td><td>1412.67</td><td>0</td><td>0.6</td></tr><tr><td>Cardiff_10_04</td><td>136</td><td>0</td><td>0.572</td><td>1</td><td>1633.32</td><td>1633.32</td><td>0</td><td>0.62</td></tr><tr><td>Cardiff_10_05</td><td>136</td><td>0</td><td>0.436</td><td>1</td><td>1682.93</td><td>1682.93</td><td>0</td><td>0.48</td></tr><tr><td>Cardiff_10_06</td><td>155</td><td>0</td><td>0.591</td><td>1</td><td>1221.75</td><td>1221.75</td><td>0</td><td>0.65</td></tr><tr><td>Cardiff_10_07</td><td>129</td><td>0</td><td>0.473</td><td>1</td><td>1931.75</td><td>1931.75</td><td>0</td><td>0.52</td></tr><tr><td>Cardiff_10_08</td><td>143</td><td>0</td><td>0.615</td><td>1</td><td>1578.15</td><td>1578.15</td><td>0</td><td>0.67</td></tr><tr><td>Cardiff_10_09</td><td>126</td><td>0</td><td>0.456</td><td>1</td><td>2552.03</td><td>2552.03</td><td>0</td><td>0.49</td></tr><tr><td>Cardiff_10_10</td><td>125</td><td>0</td><td>0.46</td><td>1</td><td>2045.18</td><td>2045.18</td><td>0</td><td>0.5</td></tr><tr><td>Cardiff_10_11</td><td>137</td><td>0</td><td>0.424</td><td>1</td><td>2095.25</td><td>2095.25</td><td>0</td><td>0.46</td></tr><tr><td>Cardiff_10_12</td><td>161</td><td>0</td><td>0.578</td><td>1</td><td>1737.25</td><td>1737.25</td><td>0</td><td>0.63</td></tr><tr><td>Cardiff_10_13</td><td>136</td><td>0</td><td>0.5</td><td>1</td><td>1485.88</td><td>1485.88</td><td>0</td><td>0.57</td></tr><tr><td>Cardiff_10_14</td><td>125</td><td>0</td><td>0.441</td><td>1</td><td>1926.62</td><td>1926.62</td><td>0</td><td>0.5</td></tr><tr><td>Cardiff_10_15</td><td>202</td><td>0</td><td>0.479</td><td>1</td><td>1814.33</td><td>1814.33</td><td>0</td><td>0.53</td></tr><tr><td>Cardiff_10_16</td><td>148</td><td>0</td><td>0.602</td><td>1</td><td>1140.05</td><td>1140.05</td><td>0</td><td>0.69</td></tr><tr><td>Cardiff_10_17</td><td>129</td><td>0</td><td>0.449</td><td>1</td><td>1733.83</td><td>1733.83</td><td>0</td><td>0.5</td></tr><tr><td>Cardiff_10_18</td><td>137</td><td>0</td><td>0.456</td><td>1</td><td>1960.93</td><td>1960.93</td><td>0</td><td>0.51</td></tr><tr><td>Cardiff_10_19</td><td>136</td><td>0</td><td>0.478</td><td>1</td><td>1533.67</td><td>1533.67</td><td>0</td><td>0.53</td></tr><tr><td>Cardiff_10_20</td><td>135</td><td>0</td><td>0.502</td><td>1</td><td>1366.42</td><td>1366.42</td><td>0</td><td>0.59</td></tr><tr><td>Cardiff_15_01</td><td>993</td><td>0</td><td>3.59</td><td>1</td><td>1332</td><td>1332</td><td>0</td><td>3.66</td></tr><tr><td>Cardiff_15_02</td><td>383</td><td>0</td><td>0.773</td><td>1</td><td>1970.62</td><td>1970.62</td><td>0</td><td>0.83</td></tr><tr><td>Cardiff_15_03</td><td>292</td><td>0</td><td>1.099</td><td>1</td><td>1971.38</td><td>1971.38</td><td>0</td><td>1.19</td></tr><tr><td>Cardiff_15_04</td><td>320</td><td>0</td><td>0.799</td><td>1</td><td>2188.32</td><td>2188.32</td><td>0</td><td>0.84</td></tr><tr><td>Cardiff_15_05</td><td>311</td><td>0</td><td>0.727</td><td>1</td><td>1831.27</td><td>1831.27</td><td>0</td><td>0.79</td></tr><tr><td>Cardiff_15_06</td><td>408</td><td>0</td><td>1.395</td><td>1</td><td>1695.17</td><td>1695.17</td><td>0</td><td>1.5</td></tr><tr><td>Cardiff_15_07</td><td>318</td><td>0</td><td>0.755</td><td>1</td><td>2050.67</td><td>2050.67</td><td>0</td><td>0.82</td></tr><tr><td>Cardiff_15_08</td><td>318</td><td>0</td><td>1.091</td><td>1</td><td>1689.55</td><td>1689.55</td><td>0</td><td>1.17</td></tr><tr><td>Cardiff_15_09</td><td>322</td><td>0</td><td>0.687</td><td>1</td><td>2765.68</td><td>2765.68</td><td>0</td><td>0.74</td></tr></tbody><tfoot><tr><td colspan="9" style="text-align: center;">Continued on next page</td></tr></tfoot></table>
impact on the gap. That is why we set it to one in most cases.

### 7.4. Sensitivity analysis

In this section, we first present the impact of the ARCC and then present the impact of the SPR criterion. We use the concept of cut-off values. Cut-off values are primal bounds that the user provides to the algorithm to prune the branch and bound tree. We use the best-known solution as a cut-off value. Although these values can be used to prove solution optimality, the goal here is to prevent random performance of the primal heuristics, enabling a fair comparison of performance with and without the tested improvement. Unless specified otherwise, the ARCC and the SPR criterion are used. We recall that the gap is calculated as follows: $\frac{PB-DB}{DB}$, with PB the primal bound and DP the dual bound.

#### 7.4.1. Impact of the adapted rounded capacity cuts

This section presents the impact of the adapted rounded capacity cuts. The capacity constraint in the Chen et al. (2021) instances is not tight, as a substantial proportion of the truck's capacity remains unused in an optimal solution. Therefore, the ARCC are not helpful for Chen et al. (2021) instances. In addition, the instances with 100 customers are all solved within 10 minutes, most of them solved at the root node. For these reasons, we focus on the instances from Gehring and Homberger (1999) in this section. As shown in table 6, within the three-hour limit, the exact branch-cut-and-price does not find a feasible solution for several instances of the Gehring and Homberger (1999) datasets. We have run the exact BCP on these instances with a ten-hour time limit to find accurate cut-off values. However, we did not find a feasible solution for all instances with 150 customers. Therefore, we present the results on the Gehring and Homberger (1999) instances with 120, 130, and 140 customers in table 8. We compare the performances with and without the ARCC. Table 8 shows the experiments' summary to measure the ARCC's impact. For each data set, we give the name of the set (Instances), whether the ARCC were used or not (ARCC), the average gap obtained on instances that are not solved optimally within the time limit (Gap), the average time in seconds required to solve the instances of the dataset (Time) and the average number of nodes processed in the branch-cut-and-price algorithm (Nodes). Note that the number of optimal solutions found is omitted as the ARCC does not impact this number. The detailed results can be found in the table 16 in Appendix B.
<table><caption>Table 10: Detailed performance of the exact algorithm on the small-size instances</caption><thead><tr><th>Instance</th><th>A<sup>max</sup></th><th>Root Gap</th><th>Root Time</th><th>Nodes</th><th>DB</th><th>PB</th><th>Gap</th><th>Time (s)</th></tr></thead><tbody><tr><td>Cardiff_15_10</td><td>271</td><td>0</td><td>0.9</td><td>1</td><td>2984.98</td><td>2984.98</td><td>0</td><td>0.94</td></tr><tr><td>Cardiff_15_11</td><td>319</td><td>0</td><td>0.824</td><td>1</td><td>2242.98</td><td>2242.98</td><td>0</td><td>0.89</td></tr><tr><td>Cardiff_15_12</td><td>379</td><td>0</td><td>0.991</td><td>1</td><td>2155.83</td><td>2155.83</td><td>0</td><td>1.06</td></tr><tr><td>Cardiff_15_13</td><td>369</td><td>0</td><td>1.512</td><td>1</td><td>1606.55</td><td>1606.55</td><td>0</td><td>1.59</td></tr><tr><td>Cardiff_15_14</td><td>300</td><td>0</td><td>4.598</td><td>1</td><td>2369.88</td><td>2369.88</td><td>0</td><td>4.76</td></tr><tr><td>Cardiff_15_15</td><td>728</td><td>0</td><td>1.228</td><td>1</td><td>1883.45</td><td>1883.45</td><td>0</td><td>1.29</td></tr><tr><td>Cardiff_15_16</td><td>522</td><td>0</td><td>1.129</td><td>1</td><td>1756.85</td><td>1756.85</td><td>0</td><td>1.22</td></tr><tr><td>Cardiff_15_17</td><td>293</td><td>0</td><td>0.583</td><td>1</td><td>1858.17</td><td>1858.17</td><td>0</td><td>0.64</td></tr><tr><td>Cardiff_15_18</td><td>330</td><td>0</td><td>0.893</td><td>1</td><td>2288.73</td><td>2288.73</td><td>0</td><td>0.99</td></tr><tr><td>Cardiff_15_19</td><td>303</td><td>0</td><td>0.69</td><td>1</td><td>1762.17</td><td>1762.17</td><td>0</td><td>0.76</td></tr><tr><td>Cardiff_15_20</td><td>306</td><td>0</td><td>0.871</td><td>1</td><td>1697.1</td><td>1697.1</td><td>0</td><td>0.95</td></tr><tr><td>Cardiff_25_01</td><td>1622</td><td>0</td><td>78.333</td><td>1</td><td>2233.87</td><td>2233.87</td><td>0</td><td>78.81</td></tr><tr><td>Cardiff_25_02</td><td>1688</td><td>0.3</td><td>47.381</td><td>3</td><td>1979.95</td><td>1979.95</td><td>0</td><td>64.68</td></tr><tr><td>Cardiff_25_03</td><td>985</td><td>0.7</td><td>28.089</td><td>3</td><td>2389.7</td><td>2389.7</td><td>0</td><td>40.9</td></tr><tr><td>Cardiff_25_04</td><td>1280</td><td>5.8</td><td>31.931</td><td>3</td><td>2797.63</td><td>2797.63</td><td>0</td><td>66.87</td></tr><tr><td>Cardiff_25_05</td><td>1481</td><td>0</td><td>27.47</td><td>1</td><td>2062.97</td><td>2062.97</td><td>0</td><td>27.72</td></tr><tr><td>Cardiff_25_06</td><td>1140</td><td>0</td><td>24.194</td><td>1</td><td>2075.48</td><td>2075.48</td><td>0</td><td>24.52</td></tr><tr><td>Cardiff_25_07</td><td>1212</td><td>1.3</td><td>31.748</td><td>3</td><td>2112.6</td><td>2112.6</td><td>0</td><td>44.93</td></tr><tr><td>Cardiff_25_08</td><td>1119</td><td>0</td><td>20.74</td><td>1</td><td>2605.48</td><td>2605.48</td><td>0</td><td>20.91</td></tr><tr><td>Cardiff_25_09</td><td>1634</td><td>1.4</td><td>33.46</td><td>3</td><td>2055.22</td><td>2055.22</td><td>0</td><td>48.93</td></tr><tr><td>Cardiff_25_10</td><td>1605</td><td>2.7</td><td>75.233</td><td>5</td><td>2217.12</td><td>2217.12</td><td>0</td><td>296.04</td></tr><tr><td>Cardiff_25_11</td><td>5207</td><td>0.7</td><td>85.909</td><td>3</td><td>2051.12</td><td>2051.12</td><td>0</td><td>99.0</td></tr><tr><td>Cardiff_25_12</td><td>1245</td><td>7</td><td>18.103</td><td>5</td><td>4775.93</td><td>4775.93</td><td>0</td><td>157.68</td></tr><tr><td>Cardiff_25_13</td><td>2730</td><td>0</td><td>57.655</td><td>1</td><td>2248.63</td><td>2248.63</td><td>0</td><td>58.04</td></tr><tr><td>Cardiff_25_14</td><td>900</td><td>4</td><td>57.477</td><td>3</td><td>2530.68</td><td>2530.68</td><td>0</td><td>257.38</td></tr><tr><td>Cardiff_25_15</td><td>1240</td><td>0.7</td><td>35.323</td><td>3</td><td>2871.6</td><td>2871.6</td><td>0</td><td>51.14</td></tr><tr><td>Cardiff_25_16</td><td>1553</td><td>1.5</td><td>29.461</td><td>3</td><td>2503.27</td><td>2503.27</td><td>0</td><td>40.64</td></tr><tr><td>Cardiff_25_17</td><td>3780</td><td>0</td><td>65.629</td><td>1</td><td>2275.83</td><td>2275.83</td><td>0</td><td>65.9</td></tr><tr><td>Cardiff_25_18</td><td>922</td><td>0.6</td><td>56.438</td><td>3</td><td>2167.9</td><td>2167.9</td><td>0</td><td>125.93</td></tr><tr><td>Cardiff_25_19</td><td>1796</td><td>1.7</td><td>55.092</td><td>3</td><td>2257.8</td><td>2257.8</td><td>0</td><td>71.74</td></tr><tr><td>Cardiff_25_20</td><td>2449</td><td>0</td><td>82.744</td><td>1</td><td>2816.63</td><td>2816.63</td><td>0</td><td>83.14</td></tr><tr><td>c101_35</td><td>2460</td><td>0.7</td><td>20.556</td><td>3</td><td>1896.5</td><td>1896.5</td><td>0</td><td>113.09</td></tr><tr><td>c102_35</td><td>2460</td><td>0.5</td><td>27.488</td><td>5</td><td>1892.1</td><td>1892.1</td><td>0</td><td>729.18</td></tr><tr><td>c103_35</td><td>2460</td><td>0.6</td><td>17.262</td><td>3</td><td>1889.2</td><td>1889.2</td><td>0</td><td>32.16</td></tr><tr><td>c104_35</td><td>2472</td><td>0.8</td><td>50.63</td><td>25</td><td>1888.6</td><td>1888.6</td><td>0</td><td>1406.2</td></tr><tr><td>c105_35</td><td>2464</td><td>0.5</td><td>38.436</td><td>3</td><td>1896.5</td><td>1896.5</td><td>0</td><td>257.74</td></tr><tr><td>c106_35</td><td>2460</td><td>0.6</td><td>34.21</td><td>5</td><td>1896.5</td><td>1896.5</td><td>0</td><td>241.3</td></tr><tr><td>c107_35</td><td>2490</td><td>0.6</td><td>37.053</td><td>3</td><td>1896.5</td><td>1896.5</td><td>0</td><td>162.77</td></tr><tr><td>c108_35</td><td>2460</td><td>0.7</td><td>28.446</td><td>5</td><td>1893</td><td>1893</td><td>0</td><td>408.04</td></tr><tr><td>c109_35</td><td>2490</td><td>0.5</td><td>121.423</td><td>21</td><td>1893</td><td>1893</td><td>0</td><td>1656.22</td></tr></tbody></table>
<table><caption>Table 8: Average impact of the adapted rounded capacity cuts on the (Gehring and Homberger, 1999) instances</caption><thead><tr><th>Instances</th><th>ARCC</th><th>Gap</th><th>Time (s)</th><th>Nodes</th></tr></thead><tbody><tr><td>120_C1_2_X</td><td>yes</td><td>0.12</td><td>2437.30</td><td>31.6</td></tr><tr><td>120_C1_2_X</td><td>no</td><td>0.12</td><td>2451.13</td><td>40.6</td></tr><tr><td>130_C1_2_X</td><td>yes</td><td>0.28</td><td>5784.12</td><td>88.2</td></tr><tr><td>130_C1_2_X</td><td>no</td><td>0.30</td><td>5744.47</td><td>67.6</td></tr><tr><td>140_C1_2_X</td><td>yes</td><td>0.56</td><td>5955.80</td><td>75.4</td></tr><tr><td>140_C1_2_X</td><td>no</td><td>0.61</td><td>5925.10</td><td>72.6</td></tr></tbody></table>

We can observe in table 8 that the ARCC have a positive impact on the gap, growing with the size of the instances. The ARCC have a negligible impact on the computation time because the capacity constraints are not very tight for these instances.

### 7.4.2. Impact of the shorter partial route criterion

In this section, we study the impact of the shorter partial route criterion in terms of the exact algorithm's efficiency and graph size. In our computational experiment, and because checking this criterion would be too time-consuming, we only consider special cases where the truck can only visit a single customer of K, and the others must be delivered by drone. As the instances from (Gehring and Homberger, 1999) have a small drone-reachable neighborhood, this criterion does not significantly impact the number of arcs. For this reason, we perform the sensitivity analysis only for the medium-size instances from (Chen et al., 2021) and (Solomon, 1987). We compare the performances of the algorithm with and without the SPR criterion. We recall that $G^d$ is the directed graph representing the delivery possibilities for a truck with d drones. Table 9 shows the experiments' summary to measure the SPR criterion's impact on the medium-size instances. For each data set, we give the name of the set (Instances), the number of arcs in subproblem graph and the proportion of arcs compared to without the SPR criterion (|$A^1$|), (|$A^2$|), (|$A^3$|), the average Gap of the instances for which a feasible solution has been found, the average time in seconds required to solve the instances of the dataset (Time) and the average number of nodes processed in the branch-cut-and-price algorithm (Nodes).

<table><caption>Table 9: Sensitivity analysis over SPR criterion</caption><thead><tr><th>Instances</th><th>SPR</th><th>$|A^1|$</th><th>$|A^2|$</th><th>$|A^3|$</th><th>Gap</th><th>Time (s)</th><th>Nodes</th></tr></thead><tbody><tr><td>Cardiff_40_X</td><td>yes</td><td>4449.15 (20%)</td><td>62169.85 (54%)</td><td>72709.35 (97%)</td><td>6.84</td><td>6376.47</td><td>10.7</td></tr><tr><td>Cardiff_40_X</td><td>no</td><td>22597.45 (100%)</td><td>116217.80 (100%)</td><td>75056.20 (100%)</td><td>7.26</td><td>6318.14</td><td>8.8</td></tr><tr><td>c10X_100</td><td>yes</td><td>12902 (68%)</td><td>18298.67 (95%)</td><td>19009.23 (100%)</td><td>0</td><td>206.81</td><td>5.23</td></tr><tr><td>c10X_100</td><td>no</td><td>19079 (100%)</td><td>19326.56 (100%)</td><td>19009.23 (100%)</td><td>0</td><td>274.63</td><td>5.89</td></tr></tbody></table>

Note that with cut-off values, the algorithm solves all c10X_100 instances to optimality. The criterion removes a substantial number of arcs in $G^1$ in both
Note that the instance *Cardiff_40_11* is not feasible under our assumptions.

<table><caption>Table 11: Detailed performance of the exact algorithm on instances adapted from (Chen et al., 2021)</caption><thead><tr><th>Instance</th><th>|A<sup>max</sup>|</th><th>K̂</th><th>Root Gap</th><th>Root Time</th><th>Nodes</th><th>DB</th><th>PB</th><th>Gap</th><th>Time (s)</th></tr></thead><tbody><tr><td>Cardiff_40_01</td><td>13034</td><td>10</td><td>0.3</td><td>214.127</td><td>3</td><td>2931.88</td><td>2931.88</td><td>0</td><td>306.84</td></tr><tr><td>Cardiff_40_02</td><td>95883</td><td>14</td><td>20.1</td><td>92.803</td><td>37</td><td>3044.72</td><td>3217.42</td><td>5.7</td><td>10800</td></tr><tr><td>Cardiff_40_03</td><td>30318</td><td>12</td><td>2.7</td><td>168.228</td><td>3</td><td>2682.5</td><td>2682.5</td><td>0</td><td>351.02</td></tr><tr><td>Cardiff_40_04</td><td>169680</td><td>15</td><td>-</td><td>401.227</td><td>9</td><td>4292.31</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>Cardiff_40_05</td><td>25801</td><td>11</td><td>1.7</td><td>329.68</td><td>11</td><td>3770.9</td><td>3770.9</td><td>0</td><td>2141.71</td></tr><tr><td>Cardiff_40_06</td><td>26451</td><td>12</td><td>12.5</td><td>203.151</td><td>19</td><td>2909</td><td>2909</td><td>0</td><td>1633.04</td></tr><tr><td>Cardiff_40_07</td><td>89523</td><td>15</td><td>-</td><td>674.039</td><td>3</td><td>2093.42</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>Cardiff_40_08</td><td>9510</td><td>9</td><td>0.6</td><td>191.295</td><td>3</td><td>3834.02</td><td>3834.02</td><td>0</td><td>358.6</td></tr><tr><td>Cardiff_40_09</td><td>81734</td><td>14</td><td>35.1</td><td>381.086</td><td>11</td><td>3580.1</td><td>4484.18</td><td>25.3</td><td>10800</td></tr><tr><td>Cardiff_40_10</td><td>55913</td><td>14</td><td>1.3</td><td>171.288</td><td>3</td><td>2692.08</td><td>2692.08</td><td>0</td><td>573.32</td></tr><tr><td>Cardiff_40_11</td><td>48327</td><td>12</td><td>0</td><td>29.336</td><td>1</td><td>unfeasible</td><td>unfeasible</td><td>0</td><td>30.12</td></tr><tr><td>Cardiff_40_12</td><td>38128</td><td>14</td><td>0.7</td><td>200.695</td><td>3</td><td>2642.87</td><td>2642.87</td><td>0</td><td>276.22</td></tr><tr><td>Cardiff_40_13</td><td>121649</td><td>14</td><td>21.9</td><td>196.029</td><td>23</td><td>2723.1</td><td>2913.13</td><td>7</td><td>10800</td></tr><tr><td>Cardiff_40_14</td><td>247124</td><td>15</td><td>-</td><td>1716.28</td><td>1</td><td>2236.37</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>Cardiff_40_15</td><td>33041</td><td>12</td><td>20.3</td><td>181.949</td><td>39</td><td>2606.74</td><td>2707.12</td><td>3.9</td><td>10800</td></tr><tr><td>Cardiff_40_16</td><td>153756</td><td>15</td><td>-</td><td>785.406</td><td>5</td><td>2927.97</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>Cardiff_40_17</td><td>24463</td><td>13</td><td>3.2</td><td>278.477</td><td>23</td><td>2733.69</td><td>2781.87</td><td>1.8</td><td>10800</td></tr><tr><td>Cardiff_40_18</td><td>93744</td><td>14</td><td>-</td><td>224.516</td><td>21</td><td>2707.07</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>Cardiff_40_19</td><td>10767</td><td>11</td><td>1.9</td><td>198.707</td><td>3</td><td>2784.05</td><td>2784.05</td><td>0</td><td>368.26</td></tr><tr><td>Cardiff_40_20</td><td>85341</td><td>14</td><td>-</td><td>214.132</td><td>23</td><td>2649.34</td><td>-</td><td>-</td><td>10800</td></tr></tbody></table>

<table><caption>Table 12: Detailed performance of the exact algorithm on 100-customer instances</caption><thead><tr><th>Instance</th><th>|A<sup>max</sup>|</th><th>Root Gap</th><th>Root Time</th><th>Nodes</th><th>DB</th><th>PB</th><th>Gap</th><th>Time (s)</th></tr></thead><tbody><tr><td>c101_100</td><td>18845</td><td>0.3</td><td>88.031</td><td>5</td><td>5298.76</td><td>5298.76</td><td>0</td><td>643.78</td></tr><tr><td>c102_100</td><td>18845</td><td>0.4</td><td>86.379</td><td>7</td><td>5298.36</td><td>5298.36</td><td>0</td><td>751.62</td></tr><tr><td>c103_100</td><td>18898</td><td>0.4</td><td>104.734</td><td>43</td><td>5296.46</td><td>5296.46</td><td>0</td><td>4008.0</td></tr><tr><td>c104_100</td><td>20064</td><td>1.4</td><td>141.282</td><td>121</td><td>5274.95</td><td>5339.06</td><td>1.2</td><td>10800</td></tr><tr><td>c105_100</td><td>18849</td><td>0.3</td><td>90.378</td><td>19</td><td>5298.76</td><td>5298.76</td><td>0</td><td>3249.86</td></tr><tr><td>c106_100</td><td>18844</td><td>0.4</td><td>90.421</td><td>9</td><td>5298.76</td><td>5298.76</td><td>0</td><td>1144.6</td></tr><tr><td>c107_100</td><td>18875</td><td>0.3</td><td>80.05</td><td>5</td><td>5298.76</td><td>5298.76</td><td>0</td><td>769.38</td></tr><tr><td>c108_100</td><td>18845</td><td>0.4</td><td>93.02</td><td>15</td><td>5295.36</td><td>5295.36</td><td>0</td><td>1712.0</td></tr><tr><td>c109_100</td><td>19018</td><td>0.4</td><td>125.412</td><td>49</td><td>5295.16</td><td>5295.16</td><td>0</td><td>5127.0</td></tr></tbody></table>
cases, less in $G^2$, and almost none in $G^3$. It removes more arcs when instances have larger drone-reachable neighborhoods.

However, while the computation efficiency improvement is mild, the SPR criterion slightly accelerates computation time in most solved instances for both data sets. The unsolved instances gaps are decreased by 0.5% on average. Although the criterion seems to have a negative impact on computation time in the *Cardiff_40_X* instances, the time decreases for 7 out of 9 instances solved optimally. Also, the criterion decreases the number of nodes in the *c10X_100* instances. The impact is positive, as it improves convergence speed. The criterion increases the number of nodes in instances *Cardiff_40_X*, which is also a positive result, as it is mainly due to the number of nodes in instances where the optimal solution was not found, demonstrating a faster node processing.

### 7.4.3. Impact of the restricted master heuristic

In this section, we present the impact of the restricted master heuristic. In this experiment, cut-off values are not provided to the solver.

<table><caption>Table 10: Sensitivity analysis over restricted master heuristic.</caption><thead><tr><th>Instances</th><th>RMH</th><th>Solved</th><th>No Sol</th><th>Gap</th><th>Time</th><th>Nodes</th></tr></thead><tbody><tr><td>Cardiff_40_X yes</td><td>yes</td><td>9/20</td><td>6</td><td>3.12</td><td>6241.96</td><td>12.2</td></tr><tr><td>Cardiff_40_X no</td><td>no</td><td>11/20</td><td>2</td><td>3.64</td><td>5621.62</td><td>39.0</td></tr><tr><td>c10X_100 yes</td><td>yes</td><td>8/9</td><td>0</td><td>0.14</td><td>3134.03</td><td>30.4</td></tr><tr><td>c10X_100 no</td><td>no</td><td>8/9</td><td>1</td><td>0.00</td><td>4184.41</td><td>57.0</td></tr><tr><td>120_C1_2_X yes</td><td>yes</td><td>7/10</td><td>3</td><td>0</td><td>4430.94</td><td>58.8</td></tr><tr><td>120_C1_2_X no</td><td>no</td><td>6/10</td><td>2</td><td>0.19</td><td>6116.54</td><td>132.40</td></tr></tbody></table>

It can be observed that the absence of the restricted master heuristic leads to faster node processing. Removing the restricted master heuristic positively impacts the results on the instances *Cardiff_40_X*. This can be explained by the fact that elementary columns are challenging to obtain in these instances due to the large drone-reachable neighborhoods and the lack of well-defined clusters. Consequently, the restricted master heuristic's solution space may be very limited or even empty, potentially leading to negative effects.

### 7.5. Large instance solution illustration

In this section, we display an optimal solution for the *c100_02.txt* instance.
<table><caption>Table 13: Detailed performance of the exact algorithm on instances from (Gehring and Homberger, 1999)</caption><thead><tr><th>Instance</th><th>A<sup>max</sup></th><th>K</th><th>Root Gap</th><th>Root Time</th><th>Nodes</th><th>DB</th><th>PB</th><th>Gap</th><th>Time (s)</th></tr></thead><tbody><tr><td>120.C1.2.1</td><td>14931</td><td>3</td><td>0.4</td><td>45.306</td><td>7</td><td>9090</td><td>9090</td><td>0</td><td>251.78</td></tr><tr><td>120.C1.2.2</td><td>14931</td><td>3</td><td>0.3</td><td>87.072</td><td>29</td><td>8920.5</td><td>8920.5</td><td>0</td><td>1827.92</td></tr><tr><td>120.C1.2.3</td><td>14931</td><td>3</td><td>1035</td><td>185.251</td><td>143</td><td>8815.06</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>120.C1.2.4</td><td>14931</td><td>3</td><td>1046.7</td><td>163.474</td><td>153</td><td>8728.17</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>120.C1.2.5</td><td>14931</td><td>3</td><td>0.5</td><td>58.236</td><td>19</td><td>8986.2</td><td>8986.2</td><td>0</td><td>1454.68</td></tr><tr><td>120.C1.2.6</td><td>14931</td><td>3</td><td>0.9</td><td>91.461</td><td>17</td><td>8988.4</td><td>8988.4</td><td>0</td><td>1899.7</td></tr><tr><td>120.C1.2.7</td><td>14931</td><td>3</td><td>0.5</td><td>67.444</td><td>29</td><td>8958.6</td><td>8958.6</td><td>0</td><td>2438.2</td></tr><tr><td>120.C1.2.8</td><td>14931</td><td>3</td><td>0.3</td><td>166.804</td><td>9</td><td>8817.9</td><td>8817.9</td><td>0</td><td>642.36</td></tr><tr><td>120.C1.2.9</td><td>14931</td><td>3</td><td>0.3</td><td>171.819</td><td>51</td><td>8790.62</td><td>8790.62</td><td>0</td><td>3394.73</td></tr><tr><td>120.C1.2.10</td><td>14931</td><td>3</td><td>1044.6</td><td>193.899</td><td>131</td><td>8741.16</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>130.C1.2.1</td><td>17614</td><td>4</td><td>0.6</td><td>61.628</td><td>11</td><td>9325.92</td><td>9325.92</td><td>0</td><td>744.67</td></tr><tr><td>130.C1.2.2</td><td>17614</td><td>4</td><td>0.3</td><td>105.671</td><td>11</td><td>9141.2</td><td>9141.2</td><td>0</td><td>787.67</td></tr><tr><td>130.C1.2.3</td><td>17614</td><td>4</td><td>1005.5</td><td>147.702</td><td>105</td><td>9054.61</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>130.C1.2.4</td><td>17617</td><td>4</td><td>1015.2</td><td>162.767</td><td>125</td><td>8973.53</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>130.C1.2.5</td><td>17614</td><td>4</td><td>0.5</td><td>87.648</td><td>29</td><td>9209.62</td><td>9209.62</td><td>0</td><td>3294.54</td></tr><tr><td>130.C1.2.6</td><td>17614</td><td>4</td><td>0.8</td><td>119.465</td><td>109</td><td>9222.4</td><td>9222.4</td><td>0</td><td>9784.0</td></tr><tr><td>130.C1.2.7</td><td>17614</td><td>4</td><td>0.6</td><td>122.121</td><td>21</td><td>9197.92</td><td>9197.92</td><td>0</td><td>2199.14</td></tr><tr><td>130.C1.2.8</td><td>17614</td><td>4</td><td>0.7</td><td>148.914</td><td>153</td><td>9053.06</td><td>9091.7</td><td>0.4</td><td>10800</td></tr><tr><td>130.C1.2.9</td><td>17614</td><td>4</td><td>1008.9</td><td>177.902</td><td>107</td><td>9029.23</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>130.C1.2.10</td><td>17614</td><td>4</td><td>1013.6</td><td>215.358</td><td>103</td><td>8985.44</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>140.C1.2.1</td><td>20366</td><td>4</td><td>0.6</td><td>82.444</td><td>51</td><td>9970.42</td><td>9970.42</td><td>0</td><td>2840.77</td></tr><tr><td>140.C1.2.2</td><td>20366</td><td>4</td><td>0.5</td><td>135.393</td><td>129</td><td>9786.12</td><td>9817.38</td><td>0.3</td><td>10800</td></tr><tr><td>140.C1.2.3</td><td>20366</td><td>4</td><td>931.4</td><td>192.591</td><td>125</td><td>9710.4</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>140.C1.2.4</td><td>20369</td><td>4</td><td>939</td><td>181.705</td><td>87</td><td>9632.6</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>140.C1.2.5</td><td>20366</td><td>4</td><td>0.4</td><td>94.988</td><td>25</td><td>9848.82</td><td>9848.82</td><td>0</td><td>1981.57</td></tr><tr><td>140.C1.2.6</td><td>20366</td><td>4</td><td>1.1</td><td>121.619</td><td>89</td><td>9812.82</td><td>9891.3</td><td>0.8</td><td>10800</td></tr><tr><td>140.C1.2.7</td><td>20366</td><td>4</td><td>0.3</td><td>152.451</td><td>21</td><td>9836.32</td><td>9836.32</td><td>0</td><td>1867.67</td></tr><tr><td>140.C1.2.8</td><td>20366</td><td>4</td><td>0.6</td><td>193.816</td><td>103</td><td>9707.63</td><td>9753.78</td><td>0.5</td><td>10800</td></tr><tr><td>140.C1.2.9</td><td>20366</td><td>4</td><td>933</td><td>215.832</td><td>107</td><td>9683.85</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>140.C1.2.10</td><td>20366</td><td>4</td><td>939.2</td><td>211.671</td><td>105</td><td>9636.41</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>150.C1.2.1</td><td>23331</td><td>4</td><td>0.4</td><td>95.102</td><td>35</td><td>10699.3</td><td>10699.3</td><td>0</td><td>1539.33</td></tr><tr><td>150.C1.2.2</td><td>23331</td><td>4</td><td>0.6</td><td>202.898</td><td>103</td><td>10551.2</td><td>10599.8</td><td>0.5</td><td>10800</td></tr><tr><td>150.C1.2.3</td><td>23331</td><td>4</td><td>1.6</td><td>160.827</td><td>137</td><td>10453.6</td><td>10611.5</td><td>1.5</td><td>10800</td></tr><tr><td>150.C1.2.4</td><td>23334</td><td>4</td><td>860.4</td><td>221.244</td><td>91</td><td>10425.2</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>150.C1.2.5</td><td>23331</td><td>4</td><td>0.6</td><td>107.043</td><td>123</td><td>10582.5</td><td>10628.2</td><td>0.4</td><td>10800</td></tr><tr><td>150.C1.2.6</td><td>23331</td><td>4</td><td>0.5</td><td>170.548</td><td>159</td><td>10561.4</td><td>10595.1</td><td>0.3</td><td>10800</td></tr><tr><td>150.C1.2.7</td><td>23331</td><td>4</td><td>0.5</td><td>131.694</td><td>113</td><td>10602.5</td><td>10602.5</td><td>0</td><td>5905.0</td></tr><tr><td>150.C1.2.8</td><td>23331</td><td>4</td><td>856.6</td><td>215.009</td><td>105</td><td>10468.6</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>150.C1.2.9</td><td>23331</td><td>4</td><td>860.5</td><td>236.564</td><td>105</td><td>10416.2</td><td>-</td><td>-</td><td>10800</td></tr><tr><td>150.C1.2.10</td><td>23331</td><td>4</td><td>865</td><td>210.899</td><td>107</td><td>10381</td><td>-</td><td>-</td><td>10800</td></tr></tbody></table>

Continued on next page
<figure><img src="image_6.png" /><figcaption>Figure 5: Drone routes (left) and Drone schedule (right) representations</figcaption></figure>

Black arrows represent truck movement. Green dashed arrows represent a drone delivery. Each path of green arrows represents the delivery order of a single drone from the first customer of the path. Multiple green arrows from a customer illustrate multiple parallel drone deliveries.

Figure 5 visually explains the drone schedule representation. The diagram on the left depicts the sequence of movements along the drone route: a truck arrives at customer 1, a drone makes a first round-trip to deliver customer 2 (labeled 1 in green), then the same drone makes a second round-trip to deliver customer 3 (labeled 2 in green), and finally, the truck leaves customer 1. The diagram on the right depicts the drone schedule representation of the same sequence. The path of green arrows represents the sequence of round trips made from its source by a drone. Therefore, a single drone performs a round trip to customer 2 and then to customer 3.

We use this scheduling representation because it allows the representation of schedules without further annotations and significantly simplifies the visualization.