Exact Methods for the Traveling Salesman Problem with Multiple Drones

Drone delivery is drawing increasing attention in last-mile delivery. Effective solution methods to solve decision-making problems arising in drone delivery allow to run and assess drone delivery systems. In this paper, we focus on delivery systems with a single traditional vehicle and multiple drones working in tandem to fulﬁll customer requests. We address the Traveling Salesman Problem with Multiple Drones (TSP-MD) and investigate the modeling challenges posed by the presence of multiple drones, which have proven to be hard to handle in the literature. We propose a compact Mixed-Integer Linear Programming (MILP) model to formulate the TSP-MD and several families of valid inequalities. Moreover, we illustrate an exact decomposition approach based on the compact MILP and a branch-and-cut algorithm. We show that this exact approach can solve instances with up to 24 customers to proven optimality, improving upon existing exact methods that can solve similar problems with up to ten customers only.


Introduction
The idea of delivering parcels with drones is becoming more and more popular. Since 2013, major transport and delivery companies, such as DHL and UPS, have been carrying out experiments to test the feasibility of parcel drone delivery, see, e.g., DHL Parcelcopter project 1 and UPS Flight Forward 1 http://www.dpdhl.com/en/media-relations/specials/dhl-parcelcopter.html project 2 . Through its Prime Air system, 3 Amazon is also working on a system that can deliver parcels in less than 30 minutes. Moreover, Google with its Alphabet section has designed X-Wing, 4 an unmanned aerial vehicle specifically designed for parcel delivery. Replacing traditional vehicles with drones to deliver parcels can allow to reduce transportation costs, decrease delivery times and carbon footprint, as well as reach areas inaccessible with traditional vehicles (see, e.g., Joerss et al. (2016), Lee et al. (2016), Kellermann et al. (2020)).
Parcel delivery with drones is also drawing increasing attention in the research community, particularly among researchers who develop solution methods for distribution problems. Distribution problems arising in parcel drone delivery can be classified into three main categories depending on how drones interact (if any) with traditional vehicles (i.e., trucks): (1) drones operate on their own, (2) drones cooperate with trucks, and (3) drones work in tandem with trucks. When drones operate on their own, they perform all deliveries without the support of other vehicles; in this type of delivery systems, the main decisions concern locating drone's facilities, such as warehouses and recharging points (see, e.g., Coelho et al. (2017), Paradiso et al. (2020), Park et al. (2020)). When drones cooperate with trucks, they work independently, in the sense that drones travel their own routes and serve some of the customers while trucks travel other routes to serve other customers; the cooperation between the two types of vehicles consists of making sure that all customers are served efficiently. An example of cooperation between drones and trucks without the two being synchronized is the Parallel Drone Scheduling Traveling Salesman Problem investigated by Murray & Chu (2015). In the last type of problems, where drones work in tandem with trucks, drones and trucks are synchronized; trucks both visit customers to deliver parcels but also act as moving stations and recharging facilities, where drones can take off and land to pick up and deliver parcels to serve their own customers and can recharge their batteries (see, e.g., Murray & Chu (2015), Agatz et al. (2018), Poikonen et al. (2019), Ha et al. (2020)).
In this paper, we focus on the last type of delivery systems, where trucks and drones work in tandem, and we focus on solution methods to design optimal combined routes for the two types of vehicles. In particular, we investigate systems where there are a single truck and multiple drones. This topic is more and more popular in the scientific literature also due to the growing interest shown by different organizations, as shown by Otto et al. (2018). To the best of our knowledge, in the scientific literature currently available, no exact methods that can solve to proven optimality instances with more than ten customers have been proposed. Several authors (see, e.g., Schermer et al. (2019), 2 http://www.ups.com/us/en/services/shipping-services/flight-forward-drones.page 3 https://www.amazon.com/Amazon-Prime-Air/b?ie=UTF8&node=8037720011 4 https://x.company/projects/wing/ Hà et al. (2020), Murray & Raj (2020), Poikonen & Golden (2020)) propose Mixed Integer Linear Programming (MILP) models to solve small instances of different distribution problems involving a truck and multiple drones, but the main contributions of these papers lie in the proposal of heuristics and metaheuristics that are suitable to solve medium and large instances of the corresponding problems.
The problem we investigate can be described as follows. A set of customer locations is given. Each customer requires the delivery of a single package. Packages are initially located at a depot. A truck and one or more drones work in tandem to deliver all packages to the customers. The truck and the drones start and end their route at the depot. The truck can serve customers along the route. At the depot and at all the customer locations visited by the truck, the drones can take off to fly towards a customer to serve and can land after delivering a package to a customer. A single customer can be served by a drone in between each take-off/landing operation. While the drones are airborne, the truck moves and can fulfill other deliveries. When the drones are not airborne, they are housed on the truck. We call this problem the Traveling Salesman Problem with Multiple-Drones (TSP-MD). Figure   1 shows a feasible solution of a TSP-MD with seven customers, a single drone, and a truck. The depot is represented with the rectangle. Customers are represented with circles. White customers are served by the truck, and grey customers are served by the drone. Straight lines represent the movements of the truck while dashed lines represent the movements of the drone. The truck leaves the depot while the drone takes off to serve Customer 2. While the drone is airborne, the truck serves Customers 1 and 3. At Customer 3, the two vehicles rejoin. The drone is launched to serve Customer 5 while the truck serves Customers 4 and 6. At Customer 6, the drone lands on the truck and is launched to serve Customer 7. The two vehicles rejoin at the depot, where their route ends. In this work, we aim at contributing to the literature on exact methods to solve distribution problems where trucks and drones work in tandem and multiple drones are used. The existing literature shows that the presence of multiple drones makes the exact solutions of these problems particularly challenging, so we study the TSP-MD as previously described, and we do not consider other side constraints that have been studied in the literature, such as limited drone flying range or the impossibility of serving some customers with the drones. As commonly done in the literature, the objective of the TSP-MD is to minimize the time required to complete all the deliveries and return to the depot. The main contributions of this paper are: (a) we formulate the TSP-MD with a compact MILP model; (b) we propose numerous sets of valid inequalities for the compact MILP; (c) we propose an exact solution method that decomposes the TSP-MD into easier problems and solves each of these problems with a branch-and-cut method; (d ) we present extensive computational experiments that show that the proposed decomposition method can solve instances with up to 24 customers to proven optimality.
The rest of the paper is organized as follows. Section 2 provides an overview of the relevant literature on the TSP-MD. Section 3 formally introduces the TSP-MD. Section 4 presents the compact MILP formulation for the TSP-MD. In Section 5, we present several families of valid inequalities for the MILP presented in Section 4. In Section 6, we describe an exact decomposition approach, based on the MILP and the valid inequalities presented in Sections 4 and 5, to solve the TSP-MD to proven optimality. Section 7 discusses the computational results achieved by the exact decomposition approach. Finally, Section 8 summarizes the main insights of our paper and outlines potential future research avenues.

Literature Review
As shown by Macrina et al. (2020), the academic literature on decision-making problems for last-mile delivery with combined trucks and drones has steadily increased since the seminal paper of Murray & Chu (2015). For an in-depth analysis of the literature on the topic, we refer the reader to the recent exhaustive reviews of Otto et al. (2018), Chung et al. (2020), and Macrina et al. (2020). In the next two sections, we review only the literature on exact solution techniques for last-mile delivery problems where drones and trucks work in tandem and are synchronized (Section 2.1) and papers investigating synchronized routing problems with one or more trucks and multiple drones (Section 2.2).

Exact Solution Techniques
Whereas the literature on heuristic methods to solve drone-related last-mile delivery problems where drones and trucks work in tandem and are synchronized is rich, few contributions on exact solution techniques are available. The main contributions, listed in chronological order, are Bouman et al. and Vásquez et al. (2020). All of these contributions address problems featuring a single truck and a single drone but different operational constraints and solution approaches, ranging from dynamic programming recursions to branch-and-bound, branch-and-cut, and branch-and-price methods. Bouman et al. (2018) investigate the Traveling Salesman Problem with Drone (TSP-D), where a   single truck and a single drone are available to serve a set of customers. The TSP-D is the singledrone version of the TSP-MD. They propose dynamic programming recursions to solve the TSP-D to optimality and test them on instances with up to 20 nodes. The results indicate that instances with 16 nodes can be solved to proven optimality and that the addition of the restriction that the truck can visit a limited number of customers while the drone is airborne significantly reduces the solution time.
Dell'Amico et al. (2019) study the Flying Sidekick TSP introduced by Murray & Chu (2015). They propose a three-index and a two-index formulation and different sets of valid inequalities that are embedded in branch-and-cut algorithms. These algorithms are tested on the 72 10-customer instances introduced by Murray & Chu (2015). The results show that the two-index formulation is the better formulation as it can solve 59 out of the 72 benchmark instances with an average computing time of less than 20 minutes.
El-Adle et al. (2019) study the TSP-D of Bouman et al. (2018) with some additional features, such as the battery life of the drone. They illustrate a MILP model with binary variables and present a combination of valid inequalities, pre-processing techniques, and other bound tightening strategies.
The computational results demonstrate that optimal solutions of instances with up to 24 nodes can be found and encouraging results can be achieved on larger instances with 28 and 32 nodes.
The TSP-D is investigated by Poikonen et al. (2019), who also consider additional constraints such as drone battery capacity and incompatibilities between customers and drones. They propose a heuristic branch-and-bound approach that provides a provably optimal solution of the problem under some circumstances. At each node of the search tree, a bound corresponding to a potential order to deliver a subset of packages is computed by solving a dynamic program. Even though the main contribution of the paper is represented by four heuristics based on branch-and-bound, the computational results show that instances with ten nodes can be solved to proven optimality.
Another contribution on the TSP-D is owed to Roberti & Ruthmair (2020), who consider a basic TSP-D and show how to model and address different side constraints, such as drone battery capacity, incompatibilities between drones and customers, launch and rendezvous times, etc. The authors present an exact branch-and-price algorithm based on a set-partitioning formulation of the problem, where the pricing problem is solved by a dynamic programming recursion that builds upon the ngroute relaxation introduced by Baldacci et al. (2011). Computational results show that benchmark instances with up to 39 customers can be solved to proven optimality. Vásquez et al. (2020) propose a MILP formulation for the TSP-D and an exact two-stage decomposition method. In the first stage, the sequence of customers served by the truck is defined. In the second stage, the operations of the drone are defined based on the solution of the first stage. The authors design a Benders-like decomposition algorithm strengthened by different families of valid inequalities. The computational results on benchmark instances demonstrate that the optimal solution of instances with up to 25 nodes can be found and that the drone speed strongly affects the computational performance of the proposed method -the faster the drone, the lower the computing time of the exact method.

Solution Approaches for Multi-Drone Problems
The literature on problems featuring multiple drones is diverse both in terms of problem settings and solution methods. The papers that address problems that are closer to the TSP-MD are owed to Schermer et al. (2019), Kitjacharoenchai et al. (2020), Murray & Raj (2020), and Poikonen & Golden (2020). We revise these four contributions in the following. The problem addresses two echelons of delivery: in the first level, trucks are routed from the depot to serve some customers, and, in the second level, drones perform routes starting and ending at the trucks. The authors propose a MILP formulation to solve small instances and two heuristics to solve large-sized instances -a route construction heuristic and a Large Neighborhood Search.
Murray & Raj (2020) develop a MILP formulation for the TSP-MD with additional side constraints.
In particular, they consider that if multiple drones are launched and retrieved at the same customer location, they must be synchronized in order to avoid collisions. The MILP they propose can solve small instances only. Therefore, the authors propose a three-phase heuristic approach to solve instances with more than eight customers. First, the set of customers is partitioned into two subsets (the customers served by the truck, and those served by the drones), and a feasible TSP-tour for the truck is built. Then, the paths of the drones are added to the TSP-tour. In the third phase, they determine the correct coordination of drones at each truck node considering launch, recovery, and service times. Poikonen & Golden (2020) consider the k-Multi-visit Drone Routing Problem, where a truck and multiple drones work in tandem. Each drone is allowed to carry multiple packages whose weight is taken into account to design the corresponding routes. They propose a constructive heuristic approach based on the solution of the TSP for the truck and the concept of operation, which is defined as a set of actions that start with the truck and all the drones being at the same (launching) location and end with all the drones rejoining the truck.
Even though the literature on solution approaches to solve problems featuring truck(s) and multiple drones working in tandem is richer and richer, to the best of our knowledge, no method is available to optimally solve instances with more than ten customers. The main contribution of the available papers lie in the proposition of heuristic and metaheuristic techniques. We aim at partly filling this gap by proposing a decomposition method based on the branch-and-cut paradigm to solve instances of the TSP-MD with up to 24 customers.

Problem Description
The TSP-MD can be formally described as follows. A complete directed graph G = (V, A) is given. The vertex set V is defined as V = N ∪ {0, 0 }, where N represents a set of n customers to serve and 0 and 0 are two copies of the depot, indicating the initial and final vertex of the vehicles' route, respectively -in the following, we also use the notation N 0 = N ∪ {0} and N 0 = N ∪ {0 }. The arc set A is defined as A truck is located at the depot. The truck is equipped with m identical drones. The truck and the drones are used to serve all customers. The time the truck takes to traverse arc (i , j ) ∈ A is indicated by t T i j , and the time the drones take to traverse arc (i , j ) ∈ A is indicated by t D i j . We assume the drones to be faster than the truck, so t D i j ≤ t T i j holds for each arc (i , j ) ∈ A. As commonly done in the literature, we also assume that the triangle inequality holds for both truck travel times, t T i j , and drone travel times, t D i j . We also assume that both truck travel times and drone travel times are symmetric, i.e., for each arc Each customer can be visited by the truck while some or all drones are airborne and the others are on-board. As commonly done in the literature, we neglect service times at the customers, but they can be easily included, if needed, by updating the drones and vehicle traveling times. The drones have limited carrying capacity, so they can serve a single customer before rejoining the truck.
The drones can take off and rejoin the truck at customer locations and the depot only. While drones are airborne, the truck is traveling, so a drone cannot take off and land at the same location while the truck stands still. The goal of the TSP-MD is to find the route of minimum duration to serve all customers by the truck and the drones while considering the synchronization between the truck and the drones. Customer 5 -when traveling from 0 to 8, the second drone is on the truck. The third drone traverses the dash-dotted path 0 − 2 − 6 − 1 − 0 and serves Customers 2 and 1. Numbers above or below the dashed, dotted, and dash-dotted lines represent the arc travel times, t D i j . To compute the duration of the solution, the departure time (Dep.T ) at each node visited by the truck has to be computed. We assume that all vehicles leave the depot at time 0. The departure time from Customer 8 is 10 because the truck arrives at Customer 8 at time 9 but has to wait for the first drone to rejoin after serving Customer 3. The departure time from Customer 6 is 15 because the truck arrives there at time 15 and does not have to wait for the third drone, which arrives at Customer 6 at time 13. The departure time from Customer 9 is 19 because the truck arrives there at time 16 but has to wait on the first drone to arrive after serving Customer 7. Finally the duration of the solution is 35 units because the truck returns to 0 at time 32, the first drone at time 33, the second drone at time 35, and the last drone at time 33.

Compact Formulation
In this section, we introduce a compact MILP to formulate the TSP-MD. This MILP is the basis of the exact approach that is presented in Section 6. Let x i j ∈ {0, 1} be a binary variable equal to 1 if the truck traverses arc (i , j ) ∈ A (no matter if there is any drone on-board or all of them are airborne), and let y i ∈ {0, 1} be a binary variable equal to 1 if the truck visits node i ∈ N . Moreover, let us call drone leg a sequence of three nodes 〈i , j , k〉 (i ∈ N 0 , j ∈ N , k ∈ N 0 ), where i is the node where a drone takes off from the truck, j is a customer served by the drone on its own, and k is the node where the drone rejoins the truck; for example, in Figure 2, the first drone performs three drone legs: 〈0, 3, 8〉, 〈8, 7, 9〉, and 〈9, 4, 0 〉.
the set of all feasible drone legs, and let z i j k ∈ {0, 1} be a binary variable equal to 1 if one of the drones performs drone leg 〈i , j , k〉 ∈ L to serve customer j . Finally, let a i ∈ R + be the departure time of the truck from node i ∈ V (a 0 is actually the maximum arrival time at 0 over all vehicles), and let w i ∈ Z + be a commodity variable representing the number of drones airborne when the truck leaves node i ∈ V . Then, the TSP-MD can be formulated as follows: The objective function (1a) aims at minimizing the total duration of the route that serves all customers. Constraints (1b) are flow conservation constraints of the truck. Constraints (1c) link the arc variables with the y-variables. Constraints (1d) ensure that the truck leaves from and returns to the depot. Constraints (1e) guarantee that each customer is served exactly once.
Constraints (1f) ensure that drone legs start and end only at locations visited by the truck. Notice that the last two terms of the right-hand side are not necessary for the correctness of the constraints.
guarantee that a drone leg is selected if and only if its corresponding launching and landing locations are visited by the truck. If the sum of the last two terms of (1f) is strictly positive, y i must be one and the first two terms of the right-hand side must be zero. Thus, constraints (1f) is a valid lifting of Constraint (1g)  Single-Drone Case. Being developed for the TSP-MD, model (1) is also valid for the special case where there is a single drone (i.e., m = 1). However, in this special case, the n × (n − 1) constraints (1f) can be replaced by the following 2n constraints Constraints (2a) (constraints (2b), resp.) guarantee that the drone can perform a drone leg starting from (ending at, resp.) a given customer i ∈ N only if customer i is visited by the truck (i.e., y i = 1).
Additional Notation. The following additional notation is used in the rest of the paper. Let A(S) ⊂ A, S ⊆ N , denote the set of arcs whose endpoints belong to the set S, i.e., other words, A(S) is the set of arcs induced by the set of customers S. Moreover, let A(S 1 , S 2 ) be the set of arcs starting from any node of the set S 1 ⊂ V and ending at any node of the set S 2 ⊂ V , i.e., We also denote by x(Â) the sum of the x-variables of the arcs belonging to the setÂ ⊆ A, i.e., x(Â) = (i , j )∈Â x i j . Finally, let y(S) be the sum of the y-variables of the customers of the set S ⊆ N , i.e., y(S) = i ∈S y i .

Valid Inequalities
Our computational experience about using formulation (1) to solve TSP-MD shows that its linear relaxation is as weak as it can be. Indeed, the lower bound provided by the linear relaxation of (1) is 0 in all the test instances we have used. Indeed, the bigM-coefficients in constraints (1m) make it feasible to set all variables a equal to 0 for highly fractional solutions of the x and z variables. Therefore, in this section, we present different sets of valid inequalities that strengthen the linear relaxation of formulation (1).

Minimum Number of Customers to Serve with the Truck
We can observe that, while the truck traverses an arc, m customers can be served with the drones if all of them are airborne. Moreover, all but the last arc traversed by the truck end at a customer served by the truck. Therefore, we can compute the minimum number of customers that the truck has to serve, n T , as n T = n−m m+1 . The following inequality forces the truck to serve at least n T customers and is valid for (1):

Symmetry Breaking
Because of the assumption that both truck and drone travel times are symmetric (see Section 3), we can observe that every feasible solution of (1) has a corresponding symmetric solution obtained by traversing its arcs in the opposite direction. The two solutions have the same duration. The following Symmetry-Breaking Constraints are valid for (1): Constraints (4) halve the set of feasible solutions of (1) by forcing the index of the first customer served by the truck not to be greater than the index of the last customer visited by the truck; Constraints (5) force the index of the first customer served by the truck to be strictly lower than the index of the last customer served by the truck.

Lower Bound to a 0
Schermer et al. (2019) observe that the total duration of any feasible route, a 0 , is bounded from below by the total travel time of the truck and the average travel time of the drones. Therefore, the following inequalities are valid for (1): Valid inequality (6a) can be strengthened by considering (part of) the time the truck has to wait for the drones. In particular, we lift valid inequality (6a) by considering the average delay, if any, caused by the drones performing drone legs 〈i * , k, j * 〉, for a given arc (i * , j * ) ∈ A, when the truck traverses arc (i * , j * ) ∈ A. Figure 2 shows two examples of this case. The first case is when the truck traverses arc (0, 8) (i.e., x 08 = 1) and the first drone performs drone leg 〈0, 3, 8〉 (i.e., z 038 = 1); as t D 03 = 4, t D 38 = 6, and t T 08 = 9, the unit of waiting time given by t D 03 + t D 38 − t T 08 = 4 + 6 − 9 = 1 can safely be added to the right-hand side of (6a). Similarly, another unit of waiting time can be added to the right-hand side of (6a) by considering the last arc traversed by the truck, (9, 0 ) and the last drone leg performed by the first drone, 〈9, 4, 0 〉.
To strengthen inequality (6a) by considering this waiting time, let us introduce a non-negative continuous variable d i ∈ R + representing a lower bound to the minimum time the truck has to wait for the drone(s) after traversing an arc leaving from node i ∈ N 0 . Moreover, let L i j ⊂ L be the set of drone legs that start from node i ∈ N 0 , end at node j ∈ N 0 , and make the truck wait for the drones if the truck traverses arc (i , j ) ∈ A, i.e., The following O(|A| + 1) inequalities can be added to (1) to lift (6a): where M i j is a large enough value that can be computed, for example, as (6a) by adding the set of d -variables to its right-hand side. Constraints (7b) set d i , i ∈ N 0 , equal to the average delay caused by the drones performing drone legs of type 〈i , k, j 〉 if x i j = 1 (and, consequently, y i = 1). Notice that constraints (7b) are meaningful, for each arc (i , j ) ∈ A, only if there exists at least one drone leg 〈i , k, j 〉 that makes the truck waiting on the drones after traversing arc (i , j ) ∈ A (i.e., only if |L i j | ≥ 1).
By observing that, for each arc (i , j ) ∈ A, the following inequality holds

Maximum Drone Legs Leaving from a Node
We can establish the following relationship between the number of drone legs 〈i , j , k〉 ∈ L departing from node i ∈ N 0 and the variables w i and y i 〈i , j ,k〉∈L Constraints (9a) stipulate that the number of drones airborne w i when the truck leaves node i ∈ N 0 is greater than or equal to the number of selected drone legs departing from node i . Constraints (9b) are on-off constraints stating that one or more drone legs departing from node i ∈ N can be selected if and only if i is served by the truck (i.e., y i = 1).

Relationship between x(V ) and y(N )
We can easily observe that the number of arcs traversed by the truck must be equal to the number of customers served by the truck plus 1. Therefore, the following inequality holds for (1)

Drone Legs Incompatible with First and Last Truck Arcs
The next two sets of valid inequalities are based on the relationship between the first (or the last) two arcs traversed by the truck and the drone legs that are incompatible with them. The two following sets of inequalities are valid for (1) x 0i + x i j + 1 m 〈s,r, j 〉∈L : Constraints (11a) relate the first two arcs traversed by the truck and the incompatible drone legs whereas constraints (11b) are about the last two arcs traversed by the truck.
Proof. We prove the validity of (11a) by considering all possible cases; the validity of (11b) can be proven with similar argumentation.
To prove the validity of (11a), we consider four cases determined by all possible values of its righthand side: 1. y i = 0 and y j = 0: as the truck does not visit customers i and j , then (a) the truck cannot traverse arc (0, i ) (i.e., x 0i = 0), (b) it cannot traverse arc (i , j ) (i.e., x i j = 0), (c) no drone legs starting from a node s ∈ N 0 \ {0, i } and ending at j can be selected (i.e., 〈s,r, j 〉∈L : s =i ,0 z sr j = 0), and (d ) drone leg 〈0, i , j 〉 cannot be selected. Therefore, whenever y i + y j = 0, the left-hand side cannot be strictly positive.
2. y i = 1 and y j = 0: as the truck serves customer i (i.e., y i = 1), it can traverse arc (0, i ) (i.e., x 0i = 1); however, because y j = 0, then (a) the truck cannot traverse arc (i , j ) (i.e., x i j = 0), and (b) no drone legs ending at j can be selected, so both 〈s,r, j 〉∈L : s =i ,0 z sr j and z 0i j must be 0. Therefore, if y i = 1 and y j = 0, the highest value the left-hand side can take is 1, which corresponds to x 0i = 1.
3. y i = 0 and y j = 1: as the truck does not serve customer i (i.e., y i = 0), it cannot traverse arcs (0, i ) and (i , j ) (i.e., x 0i = 0 and x i j = 0); moreover, because y j = 1, the number of drone legs ending at j cannot be greater then the number of drones, so 〈s,r, j 〉∈L : s =i ,0 z sr j + z 0i j ≤ m · y j . Therefore, if y i = 0 and y j = 1, the highest value the left-hand side can take is 1.
4. y i = 1 and y j = 1: as the truck serves both i and j , we show that the highest value the left-hand side can take in any feasible solution is 2, which happens in the following three cases: • x 0i = 1 and x i j = 1: as the first two arcs traversed by the truck are (0, i ) and (i , j ), no drone legs starting at a node different from 0 and i can land at j , so 〈s,r, j 〉∈L : s =i ,0 z sr j = 0. Furthermore, as the truck serves i , z 0i j must be zero.
• x 0i = 1, x i j = 0, and 1 m 〈s,r, j 〉∈L : s =i ,0 z sr j = 1: this happens if the first arc traversed by the truck is (0, i ), j is served by the truck later, and m drone legs ending at j are selected.
However, because x 0i = 1, z 0i j must be zero, so the left-hand side is equal to 2.
• x 0i = 0, x i j = 1, and 1 m 〈s,r, j 〉∈L : s =i ,0 z sr j = 1: this happens when the truck does not traverse arc (0, i ), but it traverses arc (i , j ) and m drone legs starting from a node different from 0 and i and ending at j are selected. Because x i j = 1, so z 0i j must be zero, and the left-hand side is equal to 2.

Too-Short Paths Cuts
The idea behind Too-Short Paths Cuts (TSPC) is that, for any subset of customers S ⊂ N such that |S| < n T , at most |S| x i j variables such that i , j ∈ S ∪ {0, 0 }, i = j , can be selected; indeed, as the truck must serve at least n T customers, if more than |S| of these x i j variables are selected, the truck makes a close path serving less than n T customers.
To formally describe TSPC, let us refer to the fractional solution depicted in Figure 3. The TSP-MD instance features six customers (n = 6) and a single drone (m = 1). Solid lines show strictly positive x i j variables, namely, x 02 = x 23 = x 30 = 1 and x 45 = x 56 = x 64 = 2 3 . The values of the x i j variables also determine the values of the y-variables, which are y 2 = y 3 = 1 and y 4 = y 5 = y 6 = 2 3 . Dashed lines represent strictly positive drone legs, namely, z 043 = z 053 = z 063 = 1 3 , z 310 = 1. Because n = 6 and m = 1, the truck must visit at least three customer, indeed, n T = n−m m+1 = 3. By looking at variables x 02 , x 23 , and x 30 (all of which are equal to 1), we can see that they form a close path 0 − 2 − 3 − 0 that serves two customers only. Such a path can be cut off by adding inequality x 02 + x 23 + x 30 ≤ y 2 + y 3 , which can be lifted as x 02 + x 03 + x 23 + x 32 + x 20 + x 30 ≤ y 2 + y 3 .
In general terms, TSPC can be described as follows:

Lifted Generalized Subtour-Elimination Constraints
Generalized Subtour-Elimination Constraints (GSEC) have been extensively used as valid inequalities for variants of the TSP where not all customers are to be visited, e.g., in the Undirected Selective TSP (Gendreau et al. (1998)) and the Orienteering Problem (Fischetti et al. (1998)). GSEC can be stated, among others, in the following form As in the TSP-MD, the truck is not supposed to serve all customers, constraints (13) are valid also for formulation (1). However, constraints (13) can be lifted as follows.
Let n D = n − n T be the maximum number of customers that can be served by the drones alone.

x(A(S \{r })) and y(S) = y(S \{r }). If we define S = S \{r }, constraint x(A(S))+1 ≤ y(S) reduces to
x(A(S )) ≤ y(S )−1, which corresponds to the well-known SEC for the set S and is valid because y 3 + y 5 + y 6 + y 7 + y 8 + y 9 1.6 − y 6 0.2 which is violated by 0.2. Notice that the corresponding constraint (13) obtained by removing z 063 from the left-hand side is not violated.

Maximum Outgoing Drone Legs Cuts
Maximum Outgoing Drone Legs Cuts (MODLC) limit the number of drone legs that take off within a given subset of customers S ⊂ N and land at nodes not contained in S. We explain the idea behind MODLC with the numerical example depicted in Figure  In any feasible solution where x 12 = 1 (or x 21 = 1), y 1 and y 2 must be equal to 1. However, we can also observe that, if x 12 = 1 (or x 21 = 1), there can be at most two drone legs that take off at customers 1 or 2 and end at nodes served by the truck after serving 1 or 2. Therefore, the following inequality holds for (1) and is violated by the current fractional solution:
Our computational experience shows that it is computationally convenient to add all these valid inequalities to (1)  (1) with a branch-and-cut algorithm. However, more clever separation procedures (even potentially running in polynomial time) can be devised based on the rich literature on algorithms to separate GSEC (13). Again, the most violated LGSEC (14) cut is added (if any).
As to the separation of MODLC (15)-(17b), we have noticed that they are rarely violated for sets of customers of cardinality greater than three, so we separate them, by inspection, for subsets of customers S that satisfy two conditions: (a) 1 ≤ |S| ≤ 3, and (b) S ⊆ N 0.05 (ŷ). MODLC (15)-(17b) are separated only if no violated LGSEC (14) was found, and the most violated MODLC is added (if any).

Solving the TSP-MD by Decomposition
Finding an optimal solution of the TSP-MD by solving formulation (1) with the different valid inequalities presented in Section 5 is challenging even for instances with 15-20 customers. Therefore, in this section, we propose a decomposition method that finds an optimal TSP-MD solution by solving m + 1 MILPs similar to (1), where each MILP is easier to solve than (1), and the overall required computational effort is usually lower than solving (1) as it is.
The key observation of this decomposition approach is that the set of feasible solutions of the TSP-MD can be partitioned into m + 1 sets, where each index set γ (γ = 0, 1, . . . , m) is characterized by the fact that exactly γ out of m drones serve at least two customers each. For example, if m = 2, the feasible TSP-MD solutions can be partitioned into three sets: in the first set of solutions, both drones serve at least two customers each (i.e., γ = 2); in the second set of solutions, one drone serves at least two customers, and the other drone serves no more than one customer (i.e., γ = 1); in the third and last set of solutions, both drones serve no more than one customer each (i.e., γ = 0).
Let t * (γ) be an optimal solution cost of the TSP-MD where exactly γ (γ = 0, 1, . . . , m) drones serve at least two customers each and the other m − γ drones serve at most one customer each. Then, the optimal TSP-MD cost can be expressed as In our computational experience, most of the times t * coincides with t * (m). Therefore, in the following (see 6.3), we describe an exact method that first computes t * (m) by solving a MILP similar to (1) (see Section 6.1) and then computes a lower bound to each value t * (γ), γ = 0, 1, . . . , m − 1, by solving another MILP derived from (1) (see 6.2).
6.1. Computing t * (m) Vásquez et al. (2020) present some properties of the optimal solution(s) of a problem similar to the TSP-MD but with a single drone. Inspired by these properties, we present some properties of the optimal solution(s) of the TSP-MD under the assumption that each drone is supposed to serve at least two customers. We use these properties to derive a MILP to compute t * (m) by adjusting (1).

Proposition 3. There exists an optimal solution of the TSP-MD where all the drones are airborne when
the truck is moving if the following three conditions hold: (i ) truck travel times, t T i j , and drone travel times, t D i j , satisfy the triangle inequality; (i i ) the truck is not faster than the drones, i.e., t T i j ≤ t D i j for each (i , j ) ∈ A; and (i i i ) each drone serves at least two customers.
Proof. We prove Proposition 3 by showing the example displayed in Figure 6 with seven customers and a drone. Panel (a) shows a feasible solution involving two drone legs 〈1, 6, 2〉 and 〈3, 7, 5〉. The duration of solution (a), t a , is t a = t T 01 +max{t T 12 , t D 16 +t D 62 }+t T 23 +max{t T 34 +t T 45 , t D 37 +t D 75 }+t T 50 . Let us assume the two drone legs 〈1, 6, 2〉 and 〈3, 7, 5〉 are replaced by drone legs 〈0, 6, 3〉 and 〈3, 7, 0 〉, respectively (see panel (b)  In general, whenever a drone performs a drone leg that is preceded or followed by an arc traversed by the drone on-board of the truck, such a drone leg can be replaced by another drone leg such that the drone takes off and lands at the same node (in two subsequent drone legs) without increasing the duration of the whole solution. This is the case no matter how many drones are available. Notice that this is true if and only if each drone serves at least two customers, as imposed by Condition (i i i ).
Because of Proposition 3, to compute t * (m), the following constraints can be added to (1): Constraints (18a) state that if customer i ∈ N is visited by the truck, then m drones must be airborne when the truck leaves i . Constraint (18b) ensures that the truck does not visit more than n − 2m customers as each drone is supposed to serve at least two customers. Constraints (18c) ensure that all drones take off and land at the depot at the beginning and the end of the route. Constraints (18d) are flow conservation constraints for the drones guaranteeing that all drones are airborne when the truck leaves each served customer.
If constraints (18) are added to (1), because of constraints (18a), the w-variables can be removed from the model and constraints (1g)-(1k) become redundant. Therefore, t * (m) can be computed as follows: Notice that all valid inequalities introduced in Section 5 are still valid and can be added to (19).

Computing a Lower Bound to t
To compute a lower bound, lb(γ) (γ = 0, 1, . . . , m − 1), to t * (γ), we solve a MILP derived from (19) by allowing m − γ customers not to be visited at all and γ drones to serve at least two customers each.
The MILP to compute lb(γ) is the following: Constraints (20b) replace constraints (1e) by stating that every customer is served at most once. Constraint (20c) guarantees that n − (m − γ) customers are served. Constraint (20d) replaces constraint (18b) and ensures that at most n − (m − γ) − 2 · γ = n − m − γ customers are served by the truck. Constraints (20e) replace constraints (18c) to ensure that exactly γ drones take off and land at the depot at the beginning and the end of the route. Moreover, notice that, in constraints (1f) and (1m), m can be replaced with max 1, γ to tighten the formulation.
As to the valid inequalities presented in Section 5, constraints (7a) and (10), are valid for (20) and can be separated without any changes. Constraints (6) and (9) are valid for (20), too, and can be separated, but m can also be replaced by γ to tighten the formulation. Constraints (8), (11), and MODLC (15)-(17b) are also valid for (20), but m can also be replaced by max 1, γ . Inequality (3) must be changed as the minimum number of customers that must be served by the truck if γ drones are used out of m (let us call such a value n T (γ)) depends on γ and can be computed as n T (γ) = n−m γ+1 ; thus, constraint (3) is replaced by y(N ) ≥ n T (γ). Similarly, in Symmetry Breaking constraints (4), (5) and in TSPC (12), n T must be replaced by n T (γ). Finally, in LGSEC (14), n D must be replaced by n − n T (γ).
Optimality Check. In general, lb(γ) is just a lower bound to t * (γ), and there is no guarantee that the two values coincide. However, it can be the case that lb(γ) = t * (γ) and is possible to derive an optimal solution corresponding to lb(γ) from the optimal solution of (20). We show this case in Figure 7. The top part of Figure 7 shows an optimal solution of formulation (20) for an instance with eight customers (n = 8), two drones (m = 2), and γ = 1. Customer 8 is not served in this optimal solution.
The truck serves Customers 1, 2, 3, 4, 5, in this sequence, and the only available drone performs drone legs 〈0, 6, 2〉 and 〈2, 7, 0 〉. Let a * 1 and a * 5 be the departure time of the truck from Customers 1 and 5 in this solution. If at least one of the three inequalities t D 08 . Indeed, one of the two drone legs 〈0, 8, 5〉 and 〈1, 8, 0 〉 can be added to the solution without affecting the duration of the tour.
In general terms, let (x * , y * , z * , a * ) be an optimal solution of formulation (20) of cost lb(γ). Let i f and i be the first and the last customer served by the truck in such a solution. Let a * i f (a * i , resp.) be the minimum (maximum, resp.) departure time of the truck from node i f (i , resp.) such that the duration of the tour is lb(γ) computed as Moreover, letN * be the set of m − γ customers that are not served in (x * , y * , z * , a * ). If, for all customers i ∈N * , the following inequality holds then lb(γ) = t * (γ).

Exact Solution Method
The exact method we propose to solve the TSP-MD consists of the following steps: Step 1. Compute an upper bound to t * . Solve formulation (19) with the addition of constraint y(N ) = n − 2 · m and the valid inequalities described in Section 5 to compute an upper bound to t * .
Let ub be such an upper bound.
Step 2. Compute t * (m). Solve formulation (19) with the addition of the valid inequalities described in Section 5 and by replacing M with ub. Use the optimal solution of Step 1 to hot-start the solver.
Step 3. Compute lb(γ) with γ = 0, 1, . . . , m − 1. For each value of γ between 0 and m − 1, solve formulation (20) with the addition of the valid inequalities presented in Sections 5 and 6.2. If the corresponding optimal solution does not pass the optimality check described in Section 6.2, then go to Step 5.
Return t * , and stop.
Step 5. Compute t * by using formulation (1). Solve formulation (1) with the addition of the valid inequalities described in Section 5, and by replacing M with ub. Return t * .
In our computational experience, Step 5 has proved to be necessary only in a few small instances with nine customers and at least two drones. However, in these cases, the optimal solution is quickly found by solving (1) because of the small size of the instances.

Computational Results
In this section, we discuss the computational results achieved by the exact solution method described in Section 6, provide managerial insights, and report a sensitivity analysis on key parameters and components of the proposed solution method. This section is organized as follows. In Section 7.1, we describe the test instances used in the computational experiments. Section 7.2 summarizes the results obtained on the baseline TSP-MD test instances we have selected. In Section 7.3, we show how the speed of the drone affects the performance of the solution method, and we offer some managerial insight to investigate the importance of the drone speed. Section 7.4 compares the performance of our exact solution method with the performance of the compact formulation (1) presented in Section 4. Section 7.5 shows how TSPC (12), LGSEC (14), and MODLC (15)-(17b) affect the performance of the solution method.
Both the exact solution method and the compact formulation (1) are coded in C, compiled with g++ 7.5.0, and solved with CPLEX 12.8. All experiments are conducted on a single thread of an Intel Xeon E3-1245v5 machine running at 3.50 GHz. We use the default parameters of CPLEX, except for CPX_PARAM_THREADS (set to one to use a single thread) and CPX_PARAM_EPINT (set to 1e-10 to decrease integrality tolerance). A time limit of two hours is imposed to solve each instance. All computing times reported in this section are in seconds.

Test Instances
We use a set of 40 baseline instances with 10, 15, 20, and 25 nodes (ten instances per size) derived from the instances introduced by Poikonen et al. (2019). The instances are created by randomly locating the depot and the n customers on a 50-by-50 grid, where node coordinates are uniformly distributed in the two dimensions. The geographical x-y coordinates of each node i are given as (x i , y i ). As commonly done in the literature, for each arc (i , j ) ∈ A, we compute truck travel times, t T i j , according to the Manhattan metric (i.e., t T i j = |x i − x j | + |y i − y j | ) and drone travel times, t D i j , according to the Euclidean metric (i.e., where α is a parameter that represents the ratio between the speed of the truck and the speed of the drones and is set equal to 0.5. Travel times are rounded up to the nearest integer to ensure that the triangle inequality holds. Each instance is tested with one, two, and three drones, leading to a baseline set of 120 instances. These 120 instances are used in Section 7.2. In Section 7.3, where we investigate the importance of parameter α, we use 240 more instances, still derived from the baseline instances, but by setting α = 0.333 and α = 1.  Step 3 provides a solution that is better than the optimal solution of Step 2, and T is the average computing time. About

Computational Results on the Baseline Test Instances
Step 5, Call indicates how many times Step 5 is called, and T is the average computing time.

Overall
Step 1 Step 2 Step 3 Step 5 Table 1 shows that the exact method finds an optimal solution for 97 out of the 120 instances, with an average computing time of less than ten minutes. We can see that increasing the number of customers and/or drones makes the TSP-MD more difficult to solve; this is due to higher average gaps (see column %Root) in Step 2. By looking at the difference between the gaps of columns %Root and %LP, we can also appreciate the significant contribution of separating LGSEC, TSPC, and MODLC -further details about the effect of adding these cuts are provided in Section 7.5. The right-most columns of the table show that Step 3 never finds a better solution than Step 2, and Step 5 is necessary only once (on an instance with nine customers and three drones).

Computational Results with Different Drone Speeds
In this section, we show how the speed of the drone affects the performance of the exact solution method and provide some managerial insights on the importance of this parameter.
Tables 2 and 3 summarize the results achieved by the exact solution method on the test instances when setting α = 0.333 (i.e., drones are three times faster than the truck) and α = 1 (i.e., drones are as fast as the truck), respectively. The format of these two tables is the same of Table 1.

Overall
Step 1 Step 2 Step 3 Table 2 shows that 109 instances can be solved to optimality when α = 0.333, i.e., 12 more instances than when α = 0.5, and the average computing time is less than five minutes. Moreover, by comparing Tables 1 and 2, we can see that all instances with up to 19 customers can be solved to optimality when α = 0.333, which is not the case when α = 0.5. As in Table 1, increasing the number of customers and/or drones makes the instances harder to solve due to higher gaps in Step 2.
Step 3 provides the optimal solution only once, and Step 5 is necessary only on an instance with nine customers and three drones. All in all, we can conclude that the higher the speed of the drones, the higher the chance that an instance can be solved to optimality within lower computing times.

Overall
Step 1 Step 2 Step 3   In Figure 8, we show how the average optimal solution cost changes when the number of available drones and their speed change. In the left panel of Figure 8, we compare the average solution cost over the 9-customer instances solved to optimality for any combination of values of α = 0.333, 0.5, 1 and m = 1, 2, 3; values are ratios over the average cost with m = 1 and α = 1. A similar comparison is provided in the right panel for instances with 14 customers. Figure 8 shows that the higher the number of drones and their speed, the lower the optimal solution costs. When α = 0.5, the optimal solution cost decreases by 17 to 26% compared to when α = 1, and, when α = 0.333, there is a further decrease of 6 to 10%. Having two drones instead of one allows to decrease the solution cost by 13 to 21%, and adding a third drone allows to further decrease it by 9-21%. In terms of solution cost, having a single fast drone (with α = 0.333) is similar to having three slow drones (with α = 1).

Comparison between the Exact Solution Method and the Compact Formulation
In this section, we compare the performance of the compact formulation (1) (hereafter called CF) presented in Section 4, when solved with CPLEX with the addition of all the valid inequalities presented in in Section 5, with the performance of the exact decomposition method (hereafter called Dec) presented in Section 6. As a benchmark set, we use the 40 baseline instances with two drones. Table 4 reports, for both CF and Dec, the number of instances solved to optimality (Opt), the average computing time (T), and the average computing time over the instances solved to optimality by both CF and Dec (T both ).  Table 4 shows that CF and Dec have similar performances on 9-customer instances. However, the benefits of Dec over CF stand out on larger instances with 14 and 19 customers. Indeed, both methods can solve all 14-customer instances, but Dec is nine times faster than CF. Moreover, Dec can solve three 19-customer instances more than CF and is almost three times faster. On 24-customer instances, the two approaches have similar performances. These results indicate that Dec, on average, outperforms CF.

Computational Results of the Exact Solution Method without LGSEC, TSPC, and MODLC
This section shows how separating LGSEC, TSPC, and MODLC affects the performance of the exact solution method. As we have done in Section 7.4, we use the 40 baseline instances with two drones as a benchmark set.   Table 5 shows that, on instances with up to 14 customers, the exact solution method features similar computational performance with or without cuts. On larger instances with 19 or 24 customers, removing one of the three families of valid inequalities results in fewer instances solved to optimality, larger gaps for the open instances, and usually larger computing times to find an optimal solution.
These results illustrate the benefits of separating the three families of valid inequalities.

Conclusions and Discussion
In this paper, we have investigated the Traveling Salesman Problem with Multiple Drones (TSP-MD).
We have proposed a compact formulation to model the problem and several sets of valid inequalities to improve the continuous relaxation of this compact formulation. Finding an optimal solution of test instances with just 15 nodes might take more than an hour of computing time if the compact formulation is solved as it is with an off-the-shelf solver. Therefore, we have proposed an exact decomposition approach that solves the TSP-MD to optimality by decomposing the problem into m + 1 simpler problems, where m stands for the number of drones available. This exact approach provides encouraging results and allows to solve instances with up to 24 customers and three drones in less than two hours of computing time, more than doubling the size of solvable instances of related problems with similar approaches from the literature. We have also conducted a sensitivity analysis to shed light on the advantages and computational implications of varying the number of drones available and their speed. The main insights we have offered are that (a) increasing the number of drones and their speed can significantly decrease the time to serve all customers and (b) the faster the drones, the higher the chances that the exact method can find a provably optimal solution.
Our study has been aimed at investigating the implications and challenges of having multiple drones in a truck-and-drones distributing setting, so we have focused our attention on a neat problem, without considering side constraints that arise in real-life applications and have been recently explored in the scientific literature. We expect that the proposed solution method can be adjusted to handle many side constraints, such as the flying range of the drones (possibly weight-dependent), the unavailability of customers to be served by drones, launch and rendezvous times for the drones to take off and land, etc. As many of such side constraints could be modeled by properly defining the set of feasible drone legs and the number of variables of the formulation we have proposed would be limited by these side constraints, we are confident that our exact solution method could provide comparable performance to those illustrated in our paper on instances of similar size.
We envision several future research directions to extend our study beyond the scope of this paper. First, richer problems with different side constraints could be considered, and the effectiveness of our solution approach could be investigated. Then, other objective functions, for example taking environmental considerations into account, could be examined, potentially in a multi-objective optimization setting. Finally, uncertainty could be considered to explore settings that are as close as possible to real-life settings.