Solving the Multi-Destination Shortest Path Problem
Drone delivery is rapidly becoming a practical reality — companies like Amazon, Zipline, and Wing are already operating commercial drone fleets. At the heart of every efficient drone delivery system lies a classic combinatorial optimization problem: given a set of delivery destinations, what is the shortest route that visits all of them and returns to the depot?
This is a variant of the Travelling Salesman Problem (TSP), and in this article we’ll model a concrete drone delivery scenario, solve it with both exact (ILP) and heuristic methods, and visualize the results in 2D and 3D.
Problem Setup
Imagine a drone depot located in a city center. The drone must deliver packages to 10 customer locations scattered across the city, then return to the depot. We want to minimize total flight distance.
Let the depot be node $0$ and delivery locations be nodes $1, 2, \ldots, n$. The drone visits each node exactly once. The decision variable is:
$$x_{ij} \in {0, 1}, \quad x_{ij} = 1 \text{ if the drone flies directly from node } i \text{ to node } j$$
The objective is:
$$\min \sum_{i=0}^{n} \sum_{j=0}^{n} d_{ij} \cdot x_{ij}$$
subject to:
$$\sum_{j \neq i} x_{ij} = 1 \quad \forall i \quad \text{(leave each node exactly once)}$$
$$\sum_{i \neq j} x_{ij} = 1 \quad \forall j \quad \text{(enter each node exactly once)}$$
$$u_i - u_j + n \cdot x_{ij} \leq n - 1 \quad \forall i \neq j, ; i,j \geq 1 \quad \text{(MTZ subtour elimination)}$$
where $d_{ij}$ is the Euclidean distance between nodes $i$ and $j$, and $u_i$ are auxiliary ordering variables enforcing a single connected tour (Miller–Tucker–Zemlin constraints).
The Euclidean distance between nodes is:
$$d_{ij} = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}$$
For the 3D visualization we also assign a flight altitude to each leg, reflecting real drone operations where altitude varies based on urban obstacle avoidance.
Full Source Code
1 | # ============================================================ |
Code Walkthrough
Section 0 — Environment Setup
PuLP is installed quietly. NumPy, Matplotlib (including mpl_toolkits.mplot3d), and itertools are imported. All randomness is seeded at 42 for reproducibility.
Section 1 — Problem Instance
We place the depot at coordinates $(5, 5)$ and scatter 10 customer nodes uniformly at random over a $10 \times 10$ km grid. Each node also receives a random flight altitude between 30 m and 120 m (the depot sits at 0 m), simulating real urban airspace constraints. The full $11 \times 11$ Euclidean distance matrix $D$ is precomputed once using np.linalg.norm.
Section 2 — Exact ILP Solver (MTZ)
The TSP is cast as an Integer Linear Program:
- Decision variables $x_{ij} \in {0,1}$ encode whether the drone travels from node $i$ to $j$.
- Degree constraints enforce exactly one departure and one arrival per node.
- MTZ constraints $u_i - u_j + (n-1)x_{ij} \leq n-2$ eliminate subtours by assigning an ordering integer $u_i$ to each non-depot node.
PuLP’s CBC solver finds the provably optimal solution. For $n = 11$ nodes this takes under a few seconds on Colab.
Section 3 — Nearest-Neighbour Heuristic
A greedy $O(n^2)$ baseline: always fly to the closest unvisited customer. Fast but suboptimal — it commonly produces routes 20–40% longer than optimal on random instances.
Section 4 — 2-opt Improvement
Starting from the NN tour, 2-opt repeatedly reverses a segment of the route if doing so reduces total distance:
$$\Delta = d_{i-1,i} + d_{j,j+1} - d_{i-1,j} - d_{i,j+1} < 0$$
Iterations continue until no improving swap exists. This reliably closes most of the gap between NN and ILP.
Section 5 — Running All Solvers
All three methods are timed independently. The ILP gap versus 2-opt is printed to quantify how close the heuristic came to optimality.
Visualization Panels
Panel A — ILP Optimal Route (2D). The complete optimal tour is drawn with directional arrows over the city map. The depot (gold star) and customers (cyan dots) are labelled. This is the ground-truth solution.
Panel B — Route Comparison (2D). All three routes are overlaid: grey (NN), green (2-opt), and red (ILP). The visual difference between NN and the optimized routes immediately reveals where greedy choices go wrong — typically long backtracking legs that 2-opt and ILP eliminate.
Panel C — 3D Flight Path. The ILP route is lifted into 3D using the per-node altitude assignments. Each flight leg is coloured by average altitude using the plasma colormap — warmer colours indicate higher segments. Vertical dashed poles connect each delivery point to ground level. Quiver arrows show flight direction. This panel captures the true nature of drone operations, where altitude is a continuous variable.
Panel D — Distance Bar Chart. A direct comparison of total tour length across the three methods. The difference between NN and ILP quantifies the cost of greedy decisions; the 2-opt bar shows how much of that gap is recoverable without integer programming.
Panel E — Per-leg Distances (ILP). Each individual flight leg in the ILP-optimal tour is shown as a bar, coloured by the average altitude of that leg. This reveals which segments dominate total flight time and energy consumption.
Panel F — Altitude Profile. The altitude of the drone at each stop in the ILP route is plotted as a filled line chart, giving a longitudinal view of the vertical profile of the entire delivery mission.
Results
=== Drone Delivery Route Optimization === Depot : node 0 at [5. 5.] Customers : nodes 1–10 Distance matrix (km) shape: (11, 11) [ILP] Route : [0, 10, 8, 3, 9, 4, 6, 1, 5, 2, 7, 0] [ILP] Total distance : 31.5408 km (solved in 0.18s) [NN] Route : [0, 9, 10, 8, 3, 7, 2, 5, 1, 4, 6, 0] [NN] Total distance : 34.6317 km [2opt] Route : [0, 9, 4, 6, 1, 5, 2, 7, 10, 3, 8, 0] [2opt] Total distance : 31.8671 km (solved in 0.0003s) ILP gap vs 2-opt: 0.3263 km
Figure saved.
Key Takeaways
The MTZ-based ILP formulation guarantees global optimality for small-to-medium fleets (up to ~15 nodes in seconds on Colab). The 2-opt heuristic almost always matches or comes within 1–2% of the ILP optimum at a fraction of the compute cost, making it the practical choice for real-time dispatch. For fleets with 50+ destinations, metaheuristics (genetic algorithms, ant colony optimization) or branch-and-cut solvers like Google OR-Tools become necessary.
The 3D altitude model is a reminder that real drone routing is not flat: energy consumption scales with $\sqrt{\Delta x^2 + \Delta y^2 + \Delta z^2}$, and future extensions should incorporate battery capacity constraints, no-fly zones, and wind fields into the objective.