๐Ÿšš Solving the Delivery Route Selection Problem with Python

What is the Delivery Route Selection Problem?

The Delivery Route Selection Problem (a variant of the Vehicle Routing Problem / VRP) asks:

Given a set of delivery locations and a depot, find the shortest total route that visits all locations exactly once and returns to the depot.

This is essentially the classic Traveling Salesman Problem (TSP), and itโ€™s NP-hard โ€” meaning brute force becomes infeasible as the number of stops grows.


Problem Setup

Imagine a delivery company with 1 depot and 10 delivery points. We want to find the optimal route that:

  • Starts and ends at the depot
  • Visits all 10 locations exactly once
  • Minimizes total travel distance

Mathematical Formulation

Let $n$ be the number of nodes (depot + delivery points), and $d_{ij}$ be the distance between node $i$ and node $j$.

We want to minimize the total route cost:

$$\min \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} 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)}$$

$$x_{ij} \in {0, 1}$$

Where $x_{ij} = 1$ if the route goes directly from node $i$ to node $j$.


Solution Strategy: Nearest Neighbor Heuristic + 2-opt Improvement

We use two algorithms:

  1. Nearest Neighbor Heuristic โ€” Fast greedy construction of an initial route
  2. 2-opt Local Search โ€” Iteratively reverse sub-routes to improve the solution

The 2-opt improvement works by checking if swapping two edges reduces total distance:

$$\text{If } d_{i,i+1} + d_{j,j+1} > d_{i,j} + d_{i+1,j+1} \Rightarrow \text{reverse the segment}$$


Full Source Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
# ============================================================
# Delivery Route Selection Problem
# Solver: Nearest Neighbor Heuristic + 2-opt Local Search
# Visualization: 2D route map + 3D distance surface
# ============================================================

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
import warnings
warnings.filterwarnings('ignore')

# โ”€โ”€ 1. Problem Definition โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
np.random.seed(42)

NUM_LOCATIONS = 11 # 1 depot + 10 delivery points
DEPOT_INDEX = 0

# Generate random (x, y) coordinates in [0, 100]
coords = np.random.randint(0, 100, size=(NUM_LOCATIONS, 2)).astype(float)
coords[DEPOT_INDEX] = [50.0, 50.0] # Fix depot at center

labels = ['Depot'] + [f'Stop {i}' for i in range(1, NUM_LOCATIONS)]

print("=" * 50)
print(" Delivery Route Optimization")
print("=" * 50)
print(f"Number of delivery stops : {NUM_LOCATIONS - 1}")
print(f"Depot location : {coords[DEPOT_INDEX]}")
print()

# โ”€โ”€ 2. Distance Matrix โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
def build_distance_matrix(coords):
"""Euclidean distance matrix between all node pairs."""
n = len(coords)
D = np.zeros((n, n))
for i in range(n):
for j in range(n):
diff = coords[i] - coords[j]
D[i, j] = np.sqrt(diff @ diff)
return D

D = build_distance_matrix(coords)

def total_distance(route, D):
"""Total route distance including return to depot."""
return sum(D[route[i], route[i+1]] for i in range(len(route) - 1))

# โ”€โ”€ 3. Nearest Neighbor Heuristic โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
def nearest_neighbor(D, start=0):
"""Greedy construction: always move to the closest unvisited node."""
n = len(D)
visited = [False] * n
route = [start]
visited[start] = True

for _ in range(n - 1):
current = route[-1]
nearest = None
nearest_dist = float('inf')
for j in range(n):
if not visited[j] and D[current, j] < nearest_dist:
nearest = j
nearest_dist = D[current, j]
route.append(nearest)
visited[nearest] = True

route.append(start) # return to depot
return route

# โ”€โ”€ 4. 2-opt Local Search โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
def two_opt(route, D, max_iter=1000):
"""
Improve route by reversing sub-segments.
Stops when no improvement is found or max_iter reached.
"""
best = route[:]
best_dist = total_distance(best, D)
improved = True
iteration = 0

while improved and iteration < max_iter:
improved = False
iteration += 1
for i in range(1, len(best) - 2):
for j in range(i + 1, len(best) - 1):
new_route = best[:i] + best[i:j+1][::-1] + best[j+1:]
new_dist = total_distance(new_route, D)
if new_dist < best_dist - 1e-10:
best = new_route
best_dist = new_dist
improved = True
return best, best_dist, iteration

# โ”€โ”€ 5. Solve โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
t0 = time.time()

nn_route = nearest_neighbor(D, start=DEPOT_INDEX)
nn_dist = total_distance(nn_route, D)
print(f"[Nearest Neighbor] Distance = {nn_dist:.2f} Route: {nn_route}")

opt_route, opt_dist, iters = two_opt(nn_route, D)
elapsed = time.time() - t0

print(f"[2-opt Optimized ] Distance = {opt_dist:.2f} Route: {opt_route}")
print(f"\nImprovement : {nn_dist - opt_dist:.2f} ({(nn_dist - opt_dist)/nn_dist*100:.1f}%)")
print(f"2-opt iters : {iters}")
print(f"Elapsed time: {elapsed*1000:.1f} ms")

# โ”€โ”€ 6. Visualization โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
fig = plt.figure(figsize=(18, 14))
fig.patch.set_facecolor('#0d1117')

# ---- 6-A. 2D: Nearest Neighbor Route โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
ax1 = fig.add_subplot(2, 2, 1)
ax1.set_facecolor('#161b22')

def draw_route(ax, route, coords, labels, title, color, dist):
for i in range(len(route) - 1):
a, b = coords[route[i]], coords[route[i+1]]
ax.annotate("", xy=b, xytext=a,
arrowprops=dict(arrowstyle='->', color=color,
lw=1.8, mutation_scale=15))
for idx, (x, y) in enumerate(coords):
c = '#ff6b6b' if idx == DEPOT_INDEX else '#4ecdc4'
s = 180 if idx == DEPOT_INDEX else 100
ax.scatter(x, y, s=s, c=c, zorder=5, edgecolors='white', linewidths=0.8)
offset = (3, 3) if idx != DEPOT_INDEX else (3, -8)
ax.annotate(labels[idx], (x, y), xytext=offset,
textcoords='offset points', color='white',
fontsize=7.5, fontweight='bold')
ax.set_title(f'{title}\nTotal Distance: {dist:.2f}',
color='white', fontsize=11, pad=8)
ax.tick_params(colors='#8b949e')
for spine in ax.spines.values():
spine.set_edgecolor('#30363d')
ax.set_xlim(-5, 108)
ax.set_ylim(-5, 108)

draw_route(ax1, nn_route, coords, labels, 'Nearest Neighbor Route', '#f7c59f', nn_dist)

# ---- 6-B. 2D: Optimized Route โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
ax2 = fig.add_subplot(2, 2, 2)
ax2.set_facecolor('#161b22')
draw_route(ax2, opt_route, coords, labels, '2-opt Optimized Route', '#a8ff78', opt_dist)
ax2.set_title(f'2-opt Optimized Route\nTotal Distance: {opt_dist:.2f} '
f'(โ†“{(nn_dist-opt_dist)/nn_dist*100:.1f}%)',
color='#a8ff78', fontsize=11, pad=8)

# ---- 6-C. Bar chart: edge-by-edge distance comparison โ”€โ”€โ”€โ”€โ”€โ”€โ”€
ax3 = fig.add_subplot(2, 2, 3)
ax3.set_facecolor('#161b22')

nn_edges = [D[nn_route[i], nn_route[i+1]] for i in range(len(nn_route)-1)]
opt_edges = [D[opt_route[i], opt_route[i+1]] for i in range(len(opt_route)-1)]
x_pos = np.arange(len(nn_edges))
width = 0.38

ax3.bar(x_pos - width/2, nn_edges, width, label='Nearest Neighbor',
color='#f7c59f', edgecolor='#0d1117', linewidth=0.5)
ax3.bar(x_pos + width/2, opt_edges, width, label='2-opt Optimized',
color='#a8ff78', edgecolor='#0d1117', linewidth=0.5)

ax3.set_xlabel('Edge Index', color='#8b949e')
ax3.set_ylabel('Distance', color='#8b949e')
ax3.set_title('Edge-by-Edge Distance Comparison', color='white', fontsize=11)
ax3.legend(facecolor='#21262d', edgecolor='#30363d', labelcolor='white', fontsize=9)
ax3.tick_params(colors='#8b949e')
for spine in ax3.spines.values():
spine.set_edgecolor('#30363d')

# ---- 6-D. 3D Distance Surface โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€
ax4 = fig.add_subplot(2, 2, 4, projection='3d')
ax4.set_facecolor('#0d1117')

n = NUM_LOCATIONS
I_idx = np.arange(n)
J_idx = np.arange(n)
II, JJ = np.meshgrid(I_idx, J_idx)
ZZ = D[II, JJ]

surf = ax4.plot_surface(II, JJ, ZZ, cmap='plasma',
alpha=0.85, edgecolor='none',
linewidth=0, antialiased=True)

for k in range(len(opt_route) - 1):
i, j = opt_route[k], opt_route[k+1]
ax4.plot([i, j], [j, i], [D[i,j], D[i,j]],
'w-o', lw=2.0, ms=4, alpha=0.9, zorder=10)

# โ† labelpad ใ‚’ๅ‰Š้™คใ—ใ€ใƒฉใƒ™ใƒซ่‰ฒใ‚’ๅˆฅ้€”่จญๅฎš
cbar = fig.colorbar(surf, ax=ax4, shrink=0.5, aspect=8, label='Distance')
cbar.ax.yaxis.label.set_color('#8b949e')

ax4.set_xlabel('Node i', color='#8b949e', labelpad=6)
ax4.set_ylabel('Node j', color='#8b949e', labelpad=6)
ax4.set_zlabel('Distance', color='#8b949e', labelpad=6)
ax4.set_title('3D Distance Matrix Surface\n(white line = optimized route)',
color='white', fontsize=10, pad=10)
ax4.tick_params(colors='#8b949e', labelsize=7)
ax4.xaxis.pane.fill = False
ax4.yaxis.pane.fill = False
ax4.zaxis.pane.fill = False
ax4.xaxis.pane.set_edgecolor('#30363d')
ax4.yaxis.pane.set_edgecolor('#30363d')
ax4.zaxis.pane.set_edgecolor('#30363d')
ax4.view_init(elev=30, azim=-55)

plt.suptitle('Delivery Route Optimization โ€” TSP Solver',
color='white', fontsize=15, fontweight='bold', y=0.98)
plt.tight_layout(rect=[0, 0, 1, 0.97])
plt.savefig('route_optimization.png', dpi=150,
bbox_inches='tight', facecolor='#0d1117')
plt.show()
print("\n[Plot saved as route_optimization.png]")

Code Walkthrough

1. Problem Definition

1
2
coords = np.random.randint(0, 100, size=(NUM_LOCATIONS, 2)).astype(float)
coords[DEPOT_INDEX] = [50.0, 50.0]

We create 11 nodes (1 depot + 10 stops) placed randomly on a 100ร—100 grid. The depot is fixed at the center (50, 50). np.random.seed(42) ensures the results are reproducible every run.


2. Distance Matrix

1
2
def build_distance_matrix(coords):
D[i, j] = np.sqrt(diff @ diff)

All pairwise Euclidean distances are computed once and stored in an $n \times n$ matrix $D$. Looking up $D[i,j]$ throughout the algorithm is $O(1)$, making repeated distance queries essentially free.


3. Nearest Neighbor Heuristic

1
2
def nearest_neighbor(D, start=0):
nearest = argmin D[current, unvisited]

Starting at the depot, the algorithm greedily jumps to the closest unvisited node at each step, then appends the depot again at the end to close the loop. Time complexity is $O(n^2)$. The expected tour length is approximately:

$$L_{NN} \approx 0.7124 \sqrt{n \cdot A}$$

where $A$ is the area of the bounding region. This gives a decent starting solution, but often produces โ€œcrossingโ€ edges that can be improved.


1
new_route = best[:i] + best[i:j+1][::-1] + best[j+1:]

For every pair of edges $(i \to i+1)$ and $(j \to j+1)$, the algorithm checks whether reversing the segment between them reduces total distance. The improvement condition is:

$$\Delta = d_{i,j} + d_{i+1,j+1} - d_{i,i+1} - d_{j,j+1} < 0$$

When $\Delta < 0$, the reversal is accepted and the search restarts. The loop continues until no improving swap exists. Each pass is $O(n^2)$, and in practice only a handful of passes are needed.


5. Four-Panel Visualization

Panel What it shows
Top-left Nearest Neighbor route with directed arrows
Top-right 2-opt optimized route with improvement %
Bottom-left Bar chart comparing each edgeโ€™s distance before/after optimization
Bottom-right 3D surface of the full distance matrix with the optimized route overlaid as a white line

The 3D surface is particularly insightful. The $z$-axis shows $D[i,j]$ for every node pair โ€” โ€œpeaksโ€ represent long, expensive edges. The white path overlaid on the surface traces exactly which edges the optimizer chose, visually confirming that it navigates around the high-cost peaks.


Execution Result

==================================================
  Delivery Route Optimization
==================================================
Number of delivery stops : 10
Depot location           : [50. 50.]

[Nearest Neighbor]  Distance = 383.42  Route: [0, 9, 7, 1, 10, 8, 4, 3, 5, 2, 6, 0]
[2-opt Optimized ]  Distance = 333.11  Route: [0, 4, 3, 5, 8, 1, 10, 7, 9, 6, 2, 0]

Improvement : 50.31 (13.1%)
2-opt iters : 3
Elapsed time: 2.6 ms


[Plot saved as route_optimization.png]

Key Insights

Why 2-opt works: The Nearest Neighbor heuristic frequently produces routes with crossing edges. By the triangle inequality, any two crossing edges can always be uncrossed to yield a shorter total tour. 2-opt systematically eliminates all such crossings.

Scalability: For $n \leq 20$ this runs in milliseconds. For $n \sim 100$ it remains fast. For $n > 1{,}000$ you would want Lin-Kernighan, Google OR-Tools, or simulated annealing as the $O(n^2)$ per-pass cost becomes significant.

Real-world extensions would add time windows per stop, vehicle capacity limits, and multiple vehicles โ€” transforming this into a full CVRPTW (Capacitated VRP with Time Windows). Even in those advanced solvers, 2-opt remains a standard improvement subroutine.