Solving the Parallel Drone Scheduling Traveling Salesman Problem via Constraint Programming

: Drones are currently seen as a viable way of improving the distribution of parcels in urban and rural environments, while working in coordination with traditional vehicles, such as trucks. In this paper, we consider the parallel drone scheduling traveling salesman problem, where a set of customers requiring a delivery is split between a truck and a ﬂeet of drones, with the aim of minimizing the total time required to serve all the customers. We propose a constraint programming model for the problem, discuss its implementation and present the results of an experimental program on the instances previously cited in the literature to validate exact and heuristic algorithms. We were able to decrease the cost (the time required to serve customers) for some of the instances and, for the ﬁrst time, to provide a demonstrated optimal solution for all the instances considered. These results show that constraint programming can be a very effective tool for attacking optimization problems with traveling salesman components, such as the one discussed.


Introduction
E-commerce has experienced a boom in recent decades and, according to statistics provided by Statista (https://www.statista.com/statistics/379046/worldwide-retail-ecommerce-sales/,accessed on 8 January 2023), an enormous increase in e-commerce sales is occurring worldwide.One consequence is a huge increase in home-delivered parcels, which can be readily observed in our own cities.It has been reported [1] that many e-commerce and parcel delivery companies are offering ever faster delivery, such as sameday and instant delivery.The report suggests that 20 to 25% of consumers are willing to pay more to receive their parcel on the same day, and 2% would want instant delivery.
Companies may be able to offer this kind of delivery through application of cutting-edge technology, such as the use of different types of vehicles, including aerial drones.In [2], the authors forecast that autonomous vehicles, including drones, will deliver about 80% of all parcels in the coming decades.Among the different aspects of the use of aerial drones for last mile deliveries, we focus on operational planning, aiming at optimized delivery strategies.There are several advantages associated with the use of aerial drones: they do not need to follow the road network and can fly in (almost) straight lines, and they are not affected by traffic congestion.These properties also have an important positive impact on sustainability, making these innovative solutions of interest to companies (through reduction in costs), to customers (through more efficient delivery services) and for the whole of society (through more sustainable logistics).
The seminal work of [3] pioneered a new routing optimization problem in which a truck and a drone collaborate to make deliveries.From an operational research perspective, the authors present two new prototypical variants developed from the traditional traveling salesman problem (TSP) called flying sidekick TSP (FSTSP) and parallel drone scheduling TSP (PDSTSP).In both cases a truck and drones collaborate to deliver parcels, the difference being that, in the former, model drones can be launched and collected from the truck during its tour, while, in the latter, the drones are operated directly from the central depot, while the truck executes a traditional delivery tour.In the remainder of the paper, we will focus on the second problem, addressing the interested reader to [4] for full details and some solution strategies for the FSTSP.
More specifically, in the PDSTSP, there is a truck that can leave the depot, serve a set of customers, and return to the depot, and a set of drones, each of which can leave the depot, serve a customer, and return to the depot before serving other customers.Not all the customers can be served by the drones, either due to their location or the characteristics of their parcel.The objective of the approach is to minimize the time needed to serve all the customers and to get the last vehicle back to the depot.
A first mixed-integer linear programming (MILP) model for the PDSTSP and some heuristic methods [5,6] are proposed in [3].A more refined model, and the first metaheuristic method, based on a two-step strategy embedding a dynamic programming-based component, are discussed in [7] (later extended in [8] for a variation of the problem using several trucks).A similar two-step approach, but based on matheuristics concepts [9] is presented in [10].A hybrid ant colony optimization metaheuristic is discussed in [11].An improved variable neighbour search metaheuristic is presented in [12].Most of the work published to date has focused on heuristic methods, while exact methods have been able to solve systematically only instances with up to 20 customers [10].Two other recent studies have dealt with variations of the problem where several trucks are employed, and also present results for the basic PDSTSP.In [13], three MILP models are proposed, one of which is arc-based and the other two are set-covering-based, together with a branch-and-price approach using one of the set-covering-based models.A heuristic version of one of the solvers is also discussed, leading to a matheuristic approach developed to target large instances.A more realistic variation of the PDSTSP is introduced in [14], where concepts such as the capacity of the vehicles, load balancing and decoupling of costs and times are taken into account.The authors propose an MILP model and a ruin-and-recreate metaheuristic for the problem.Several other PDSTSP variants and extensions have also been introduced and investigated in the literature.We refer the interested reader to [15,16] for a complete survey.The reader interested in technological details about hardware design for UAVs can refer to [17,18].Solutions for autonomous flying for drones are described in [19].Finally, policies and regulations for flying drones are analyzed in [20].
In this paper, we explore the potential of constraint programming (CP) for the PDSTSP.The use of CP is motivated by recent developments in black-box solvers for such a paradigm.Although still not widely used in the optimization community, these solvers can now benefit greatly from parallel computation, and offer a new perspective for optimization on modern personal computers.We discuss how to use CP efficiently and show how the paradigm can be used to optimally solve (mostly for the first time) instances traditionally used to benchmark exact and heuristic PDSTSP approaches in the literature.In earlier research studies attention was focused mainly on heuristic methods, since initial studies indicated that exact algorithms were not likely to be suitable for large instances.In this paper, we show that CP can be a "game-changer" for PDSTSP and potentially for other optimization problems with similar characteristics.
The rest of the paper is organized as follows: PDSTSP is formally defined in Section 2; in Section 3, a CP model for the PDSTSP, that represents the main contribution of the paper, together with the results, is proposed; Section 4 presents extensive experimental results relating to validation of the proposed approach; our conclusions are drawn in Section 5.

The Parallel Drone Scheduling Traveling Salesman Problem
The PDSTSP can be represented on a complete directed graph G = (V, A), where the node set V = {0, 1, ..., n} includes the depot (node 0) and a set of customers C = {1, ..., n} to be served.The set A contains ordered pairs of nodes, representing connections.A truck and a set D of identical drones are available to deliver parcels to the customers.The truck starts from the depot 0, visits a subset of the customers, and returns back to the depot.The drones operate back-and-forth trips from the depot to customers, delivering one parcel per trip.Not all the customers can be served by a drone due to the weight of the parcel, excessive distance of the customer location from the depot, or adverse morphological configurations of the terrain.Let C D ⊆ C denote the set of customers that can be served by drones.These customers are referred to as drone-eligible in the remainder of the paper.The travel time incurred by the truck to go from node i to node j is denoted as t T ij , while the time required by a drone to serve a customer i (back-and-forth trip) is denoted as t U i .The truck and the drones start from the depot at time 0; the objective of the PDSTSP is to minimize the time required to complete all the deliveries and to have the truck and all the drones back at the depot.Note that, since the truck and drones work in parallel, the objective function translates into minimizing the maximum time required by the vehicles.In the remainder of the paper, we refer to this quantity as cost.
An example of a PDSTSP instance is provided in Figure 1.An example of PDSTSP is provided on the left; we assume that two drones are available (travel times are omitted for the sake of simplicity).A solution is provided on the right.The black arcs represent the tour of the truck that visits nodes 1 and 3 before going back to the depot.The red arcs depict the missions of the first drone (that serves nodes 2 and 6), while the blue arcs depict the missions of the first drone (that serves nodes 4 and 5).

A Constraint Programming Model
The PDSTSP can be formally described by the following CP model.The variables are defined as follows: A binary variable y ij , with i, j ∈ V, takes value 1 if node i is visited right before node j by the truck during its tour, and value 0 otherwise.Note that, in case a customer i ∈ C D is visited by a drone, then y ii is set to 1, 0 otherwise.A binary variable x di , with d ∈ D and i ∈ C D , takes value 1 if the drone-eligible customer i is visited by the drone d, and value 0 otherwise.Finally, a variable α is used to capture the objective function value, which is the time taken by the last vehicle to go back to the depot after completion of its tasks.
(PDSTSP) min α (1) The objective function (1) minimizes α that will be assigned the time of the longest tasks among the vehicles.Constraint (2) implements a circuit that imposes y variables to take values in such a way as to form a feasible truck tour, and eventually self-loops for the variables associated with customers not visited by the truck.Constraint (3) requires that, if a customer i ∈ C D is not visited by the truck (which means y ii = 1), then a drone d must visit it.Constraint (4) sets α to be at least as large as the length of the truck tour.Constraint (5) requires α to take a value larger than or equal to the time spent by each drone d to execute the tasks assigned to it.Finally, constraints ( 6) and ( 7) define the domain of the binary variables.

Computational Experiments
The constraint programming model described in Sections 3 was implemented in Python 3.9 and solved via the CP-SAT solver of Google OR-Tools (https://developers. google.com/optimization,accessed on 8 January 2023) 9.4.Note that this CP solver, like most modern versions, relies heavily on bounds from linear programming during computations.The experiments were run on a laptop computer equipped with 8GB of RAM, an Apple M1 chip and using all eight CPU cores available-four performance cores running up to 3.204 GHz and four efficiency cores running up to 2.064 GHz-unless otherwise indicated.
The outcomes of the extensive experimental program that was undertaken are discussed in the remainder of this section.
A quite vast literature related to solving methods for the PDSTSP exists.We compare our results with the exact methods (namely, MILPs attacked via commercial solvers) presented in [3,10,13], and, for the larger instances, with the heuristic methods discussed in [7,[10][11][12]14].The results are organized by benchmark sets, after some preliminary considerations are presented and experiments on solver-related aspects are described.The optimal solutions obtained in the study are available in the Supplementary Materials.

Preliminary Experiments
In this section, we report the results of some preliminary experiments aimed at tuning a parameter of the solving procedure implemented and understanding the behaviour of the CP-SAT solver while running in a multi-threading environment.

Scaling and Precision
For technical reasons, the CP solver requires the use of integer coefficients in the constraints; then, scaling is necessary to represent floating point values for the coefficients that appear in constraints ( 4) and (5) of our model.In practice, the coefficients need to be multiplied by a constant F and, then, only the integer part is considered.Once the optimal solution is found, the value of α has to be scaled back by dividing by F. The choice of the scaling factor F has a double impact on the performance of the solver: high values for the factor guarantee higher numerical precision in the results (once scaled back), but, at the same time, lower numerical precision leads to faster execution.Therefore, a trade-off has to be found.
We carried out a preliminary study to identify a suitable value for the factor F that would then be adopted for the rest of the experiments.A representative instance to illustrate our results was identified in berlin52_0_80_2_2_1 (originally proposed in [7] and with an optimal solution of cost 5290.652).The results obtained by considering F ∈ {1; 10; 100; 1000; 10, 000; 100, 000; 1, 000, 000} are reported in Figure 2. In the first chart, we report the computation time (in seconds) required by the solver for different values of F; in the second chart, the value of the optimal cost affected by precision errors (estimated), and the value of the real cost of the same solution, are plotted.From Figure 2, it can be observed that smaller scaling factors correspond to shorter computing times, both to retrieve the optimal solution and to confirm its optimality.On the other hand, the precision associated with a smaller scaling factor is not sufficient, as shown by the poor quality of the estimation of the cost of the tentative optimal solution retrieved.The results presented are illustrative and representative.Based on the results, we chose to use F = 10, 000 for the remaining experiments.This value guarantees acceptable computation times and sufficiently safe precision for floating point numbers.Note that, from the results, the choice can appear as over-conservative, since values such as 100 and 1000 for F have already been shown to be good approximations; however, since the computation time associated with 10, 000 is only marginally increased, we preferred it, to be sure to avoid any numerical instability problems in our results.

Multi-Threading Computation
CP solvers notoriously benefit greatly from the use of multi-threading computation, typically substantially more than MILP solvers, due to the different nature of the solving routines.On the other hand, nowadays, home-computers and laptops offer multi-core architectures, making parallel computation affordable.In this section, we assess the effect of multi-threading on the performance of solvers while running on the PDTSP.
For the test reported, we selected the instance eil101_0_0_1_2_1 first proposed in [7] since it reduces to a classic traveling salesman problem, with no drone-eligible customers.This also enables previously proposed MILP-based solvers for the PDSTSP to close the instance.A maximum running time of 180 s was set.We solve the CP model described in Section 3 with the CP-SAT solver and the MILP model described in [10] attacked by the Gurobi (https://www.gurobi.com,accessed on 5 January 2023) 9.1.1solver.For both cases, we allow the use of a number of parallel threads between 1 and 8; in Figure 3, we report on the time required by the solvers (180 is reported when the time limit is reached).The results reported in Figure 3 indicate that the CP approach is able to benefit substantially from the use of multi-threading computation.This result was expected, given the general characteristics of routines solving constraint programming and their scalability properties.More surprising is the drastic change in performance between four and five threads.This behaviour was observed in several other instances (although the jump was not always between four and five) and is likely a consequence of activation of crucial sub-solvers on newly available parallel cores.In general, the chart also testifies to the importance of the heuristic routines used to distribute the computation economically and efficiently and to select the alternative solvers that are eventually run in parallel.On the other hand, the MILP solver does not seem to benefit from multi-threading in these short tests, possibly due to the overhead of distributing tasks and collecting results and to the reduced variety of the methods run in parallel.It should be noted that the situation might be less extreme when more challenging instances are considered, but this is difficult to observe, since, even with longer computation times (hours) allowed and eight threads, generic instances were not closed.
In conclusion, multi-threading appears to be the basis of the success of modern CP solvers; therefore, all the experiments reported in Section 4.2 were run on eight cores.

Main Experimental Results
In this section, we consider the instances commonly adopted in the previous PDSTSP literature.We also present optimal solution costs and the computation times required by the CP model discussed in Section 3 to retrieve these solutions.
It should be noted that the optimal solution costs were only already published for the 720 small instances discussed in Section 4.2.1 and for two of the instances addressed in Section 4.2.2.For all remaining 96 instances, optimality was demonstrated here for the first time.

Instances Proposed in [3]
The first instances analysed were those originally proposed in [3].A total of 120 configurations with 10 customers and 120 configurations with 20 customers were created; these instances were solved with a single truck and either one, two, or three drones, resulting in 720 test instances.
In Table 1, we compare the results published in [3,10,13] with the new CP-based approach.The rows give the average values over the 360 instances with 10 customers and the 360 instances with 20 customers.For each instance size, and for each algorithm, we give the average computing time in seconds required to solve the instances and the standard deviations.It should be noted that the experiments were run on different hardware configurations, so the comparison might not be fair.In addition, the experiments for the methods based on MILP models discussed in [10,13] were run on a single core, while those from [3] and our CP model were run using eight cores.This might make the comparison even less fair, but it has also to be considered that, according to the experiments reported in Section 4.1.2,it is very unlikely that the MILP solvers could achieve a real substantial advantage by running in parallel for the problem considered.This is not true for the CP solver, which is designed to exploit as much as possible from any multi-core architecture (in our case, a laptop computer).The results presented in Table 1 highlight the substantial improvements in exact methods for the PDSTSP that have occurred in the last few years.The CP model discussed in Section 3 (bold entries) is the fastest, with a clear advantage (especially in terms of standard deviation) for the instances with 20 customers.

Instances Proposed in Mbiadou Saleu et al. in [7]
A second program of experiments to test the new CP model was carried out on the instances originated from TSPLIB [21] in [7].For each original TSPLIB instance considered, several PDSTSP instances were created; we refer the reader to [7,14] for the details.The results cover instances with 48 to 229 customers (first number in the name of each instance); the results are presented for several methods (we show the best results disclosed in a single publication where alternative methods are discussed) together with those of the new CP model.It should be noted that the latter is the first effective exact solving method used for these instances.The results are grouped by original TSPLIB instance and are summarized in Tables 2-7.For each heuristic method, the best upper bound is reported, while, for the CP method, the optimal solution cost, the time required to retrieve the optimal solution (Sec found), and the time to prove its optimality (Sec opt) are annotated.In the tables, we report in italics any suboptimal result, and, in bold, newly improved upper bounds.In this case, the results are taken from the original papers and, therefore, are likely to have been obtained on different hardware configurations.
Examining the results reported in Tables 2-7, several observations can be made.The first is that the heuristic methods are extremely effective.The last methods particularly retrieved almost all the optimal solutions (without proving their optimality).The CP solver was able to marginally improve one solution over the previously best known results.It should be noted that this result indirectly validates the quality of the heuristic methods previously proposed.The times required to solve instances to optimality varied from a few seconds for the smaller instances to more than 17,000 s for the largest and most difficult instance.

Instances Proposed in [13]
A third set of experiments to test the model presented in Section 3 was carried out on the instances originated from TSPLIB [21] in [13].Since this dataset partially overlaps with that discussed in Section 4.2.2,only the results on the uncovered instances are reported here.We rename the instances according to the logic followed in previous publications for the instances presented in [7].We refer the reader to [13] for full details about the instances (the mapping of the naming notation should also be straightforward).The results again cover instances with 48 to 229 customers served with six drones.We report figures for the methods previously applied on these instances, namely a MILP, a branch-and-price approach (B&P) and a matheuristic algorithm (matheuristic), all discussed in [13].The results are summarized in Table 8.
The results show that CP was able to find three new optimal solutions, with improvements of up to 15% over the previous solution values.For the hardest instances, however, significant computing time was required-more than 32 h of computation.
In Figure 4 the optimality gaps are reported for all the methods available for these instances.For each method M, the gap is calculated in terms of the percentage deviation from the optimal solution as cost(M)−cost(CP) cost(CP) • 100, where cost(M) is the cost of the best solution retrieved by method M. From the chart, it is possible to see how the methods previously introduced in [13] have difficulty finding high quality solutions for the larger instances, with substantial gaps with respect to CP in some cases.From the results on these instances, the best method among those discussed in [13] appears to be the MILP-based method, although it performed poorly on one of the instances.

Conclusions
The aim of the investigation was to show the potential of constraint programming solvers on optimization problems that-like the parallel drone scheduling traveling salesman problem-extend the concepts of a traditional traveling salesman problem.The new constraint programming model we propose leads-by exploiting multi-threading computation-to a very effective exact algorithm for the problem considered.The method was able to find the optimal solution for all the instances normally considered in the literature, and this was the first optimality proof for most of the instances.During the process, four of the previously best-known solution costs were also improved.
Finally, we believe that straightforward adaptations of the model we propose could lead to state-of-the-art results for many of the several extensions proposed to the basic problem considered here, although this is beyond the scope of the present paper.From a more general perspective, we also believe that constraint programming has potential for use in effectively attacking many other optimization problems arising in different fields.

Figure 1 .
Figure 1.An example of PDSTSP is provided on the left; we assume that two drones are available (travel times are omitted for the sake of simplicity).A solution is provided on the right.The black arcs represent the tour of the truck that visits nodes 1 and 3 before going back to the depot.The red arcs depict the missions of the first drone (that serves nodes 2 and 6), while the blue arcs depict the missions of the first drone (that serves nodes 4 and 5).

Figure 2 .
Figure 2. Times and optimal objective function cost for different values of the scaling factor Finstance berlin52_0_80_2_2_1 from [7].

Figure 3 .
Figure 3.Time required by the MILP solver and by the CP solver with different numbers of threads (180 s maximum time)-instance eil101_0_0_1_2_1 from [7].

Figure 4 .
Figure 4. Optimality gap for the different methods tested on the instances proposed in [13].