Vehicle routing — from exact optimization to learning-guided heuristics¶
Machine Learning in Operations Research · combinatorial optimization, and how ML accelerates it¶
The other families in this section optimized how much to provision under uncertainty. Routing is a different beast: a combinatorial problem where the hard part is the astronomically large set of possible solutions. We tackle the Capacitated Vehicle Routing Problem (CVRP) — the core of last-mile delivery — and walk the full ladder OR/ML practitioners actually climb:
- Exact optimization — a mixed-integer program (MIP) that returns the provably optimal routes, and shows exactly why exactness stops scaling.
- A classic heuristic — Clarke-Wright savings + 2-opt local search: fast, and good enough at scale.
- A production solver — Google OR-Tools' routing engine with guided local search, the industrial-strength reference.
- Learning-guided optimization — the current ML×OR frontier: learn which edges belong in good routes, prune the problem to those candidates, and optimize the sparse remainder — predict-then- optimize, applied to the structure of the search itself.
The problem. A depot dispatches capacity-limited vehicles to serve customers with known demands; minimize total distance so every customer is visited once and no vehicle exceeds capacity. We use reproducible synthetic delivery instances (clustered "city" geographies).
import os, time, math, itertools
os.environ.setdefault("KMP_DUPLICATE_LIB_OK", "TRUE")
import numpy as np, pandas as pd, matplotlib.pyplot as plt
BLUE, RED, GREEN, ORANGE, GREY, PURP = "#2b6cb0","#c53030","#2f855a","#dd6b20","#a0aec0","#6b46c1"
def make_instance(n_cust, seed=0, n_clusters=4, cap=30):
rng = np.random.default_rng(seed)
centers = rng.uniform(10, 90, (n_clusters, 2))
pts = [np.array([50.0, 50.0])] # depot at centre
for _ in range(n_cust):
c = centers[rng.integers(n_clusters)]; pts.append(np.clip(c + rng.normal(0, 8, 2), 0, 100))
pts = np.array(pts)
demand = np.array([0] + list(rng.integers(1, 10, n_cust)))
dist = np.linalg.norm(pts[:, None] - pts[None, :], axis=2)
K = int(np.ceil(demand.sum() / cap)) # min vehicles by capacity
return dict(pts=pts, demand=demand, dist=dist, cap=cap, K=K, n=len(pts))
def route_distance(routes, dist):
return sum(dist[0, r[0]] + sum(dist[r[k], r[k+1]] for k in range(len(r)-1)) + dist[r[-1], 0]
for r in routes if len(r))
def plot_routes(inst, routes, title, ax):
pts = inst["pts"]
ax.scatter(pts[1:,0], pts[1:,1], s=18, c=GREY, zorder=2)
ax.scatter(*pts[0], s=180, marker="*", c=RED, zorder=3, label="depot")
for r, c in zip(routes, plt.cm.tab10(np.linspace(0,1,max(len(routes),1)))):
path = [0]+list(r)+[0]
ax.plot(pts[path,0], pts[path,1], "-", color=c, lw=1.3, alpha=.8)
ax.set_title(title, fontsize=10); ax.set_xticks([]); ax.set_yticks([])
inst_s = make_instance(9, seed=1, cap=25) # small, for the exact MIP
print(f"small instance: {inst_s['n']-1} customers, capacity {inst_s['cap']}, min vehicles K={inst_s['K']}")
small instance: 9 customers, capacity 25, min vehicles K=2
1 · The instance — last-mile delivery¶
A depot (star) and customers (dots) scattered in clusters, each with a delivery demand; a fleet of capacity-limited vehicles must cover them all. The objective is total distance driven. Even stating the optimal solution is hard: with $n$ customers the number of possible route-sets grows factorially, which is what makes routing a showcase for both exact optimization and clever heuristics.
Note on the data. Unlike the other OR families (which use real retail and ED data), these routing instances are synthetic — reproducible clustered geographies generated with a fixed seed — because clean public last-mile datasets are scarce and unwieldy. The methods transfer unchanged to standardized academic benchmarks (CVRPLIB / Augerat) or a firm's own delivery records; only the coordinates change.
2 · Exact optimization — a mixed-integer program¶
For a small instance we can find the provably optimal routes. Let $x_{ij}=1$ if a vehicle travels directly from $i$ to $j$. Minimize $\sum_{ij} d_{ij}x_{ij}$ subject to: each customer has exactly one arc in and one out; exactly $K$ arcs leave and enter the depot; and Miller-Tucker-Zemlin (MTZ) load variables $u_i$ that both eliminate subtours and enforce capacity:
$$u_j \ge u_i + q_j - Q\,(1-x_{ij}),\qquad q_i \le u_i \le Q.$$
We solve it with PuLP/CBC. This is exact — but the binary $x_{ij}$ make it NP-hard: solve time explodes with $n$, so it is viable only for tens of customers. That limit is the whole reason the rest of the notebook exists.
import pulp
def solve_exact(inst, time_limit=60):
n, d, q, Q, K = inst["n"], inst["dist"], inst["demand"], inst["cap"], inst["K"]
N = range(n); C = range(1, n)
p = pulp.LpProblem("cvrp", pulp.LpMinimize)
x = {(i,j): pulp.LpVariable(f"x_{i}_{j}", cat="Binary") for i in N for j in N if i!=j}
u = {i: pulp.LpVariable(f"u_{i}", lowBound=q[i], upBound=Q) for i in C}
p += pulp.lpSum(d[i,j]*x[i,j] for (i,j) in x)
for j in C: p += pulp.lpSum(x[i,j] for i in N if i!=j) == 1 # in-degree
for i in C: p += pulp.lpSum(x[i,j] for j in N if i!=j) == 1 # out-degree
p += pulp.lpSum(x[0,j] for j in C) == K # K vehicles leave depot
p += pulp.lpSum(x[i,0] for i in C) == K # K return
for i in C:
for j in C:
if i!=j: p += u[j] >= u[i] + q[j] - Q*(1 - x[i,j])
t = time.time(); p.solve(pulp.PULP_CBC_CMD(msg=0, timeLimit=time_limit)); t = time.time()-t
# rebuild routes from arcs (depot has K successors -> handle separately from the unique customer succ)
succ = {i: j for (i,j) in x if i!=0 and x[i,j].value() and x[i,j].value() > 0.5}
starts = [j for (i,j) in x if i==0 and x[i,j].value() and x[i,j].value() > 0.5]
routes = []
for j0 in starts:
r=[j0]; nxt=succ.get(j0)
while nxt is not None and nxt!=0: r.append(nxt); nxt=succ.get(nxt)
routes.append(r)
return routes, route_distance(routes, d), t, pulp.LpStatus[p.status]
routes_x, dist_x, t_x, status_x = solve_exact(inst_s)
fig, ax = plt.subplots(figsize=(5.5, 5.5)); plot_routes(inst_s, routes_x, f"Exact MIP optimum — {dist_x:.1f} ({status_x})", ax)
ax.legend(fontsize=8); plt.show()
print(f"exact optimum: distance {dist_x:.1f} over {len(routes_x)} routes | solved in {t_x:.1f}s ({status_x})")
exact optimum: distance 217.2 over 2 routes | solved in 1.1s (Optimal)
3 · A classic heuristic — Clarke-Wright savings + 2-opt¶
Exact solving dies past a few dozen customers, so for realistic sizes we use heuristics. Clarke-Wright savings builds routes by greedily merging the pairs that save the most detour ($s_{ij}=d_{0i}+d_{0j}-d_{ij}$) while respecting capacity; then 2-opt un-crosses each route by reversing segments that shorten it. Together they give a good solution in milliseconds, even for hundreds of stops.
def clarke_wright(inst):
d, q, Q, n = inst["dist"], inst["demand"], inst["cap"], inst["n"]
routes = {i: [i] for i in range(1, n)}; rof = {i: i for i in range(1, n)}; load = {i: q[i] for i in range(1, n)}
sav = sorted(((d[0,i]+d[0,j]-d[i,j], i, j) for i in range(1,n) for j in range(i+1,n)), reverse=True)
for s, i, j in sav:
ri, rj = rof[i], rof[j]
if ri==rj or load[ri]+load[rj] > Q: continue
Ri, Rj = routes[ri], routes[rj]
if Ri[-1]==i and Rj[0]==j: new = Ri+Rj
elif Ri[0]==i and Rj[-1]==j: new = Rj+Ri
elif Ri[-1]==i and Rj[-1]==j: new = Ri+Rj[::-1]
elif Ri[0]==i and Rj[0]==j: new = Ri[::-1]+Rj
else: continue
for k in new: rof[k]=ri
routes[ri]=new; load[ri]+=load[rj]; del routes[rj]
return list(routes.values())
def two_opt(r, d):
def rd(rt): path=[0]+rt+[0]; return sum(d[path[k],path[k+1]] for k in range(len(path)-1))
best=rd(r); improved=True
while improved:
improved=False
for a in range(len(r)-1):
for b in range(a+1, len(r)):
nr = r[:a]+r[a:b+1][::-1]+r[b+1:]; nd=rd(nr)
if nd < best-1e-9: r, best = nr, nd; improved=True
return r
def savings_2opt(inst):
t=time.time(); routes=[two_opt(r, inst["dist"]) for r in clarke_wright(inst)]
return routes, route_distance(routes, inst["dist"]), time.time()-t
# validate on the small instance (compare to the exact optimum), then run at scale
r_h, d_h, t_h = savings_2opt(inst_s)
print(f"heuristic on small instance: {d_h:.1f} vs exact {dist_x:.1f} ({(d_h/dist_x-1)*100:+.1f}% gap, {t_h*1000:.0f} ms)")
inst_L = make_instance(100, seed=2, cap=40)
r_hL, d_hL, t_hL = savings_2opt(inst_L)
fig, ax = plt.subplots(figsize=(5.5, 5.5)); plot_routes(inst_L, r_hL, f"Savings+2-opt (100 customers) — {d_hL:.0f}", ax); plt.show()
print(f"heuristic on 100 customers: distance {d_hL:.0f} over {len(r_hL)} routes | {t_hL*1000:.0f} ms")
heuristic on small instance: 217.2 vs exact 217.2 (+0.0% gap, 0 ms)
heuristic on 100 customers: distance 1178 over 12 routes | 4 ms
4 · A production solver — Google OR-Tools¶
Industry uses purpose-built routing engines. OR-Tools constructs an initial solution (cheapest-arc) and then runs guided local search — escaping local optima by penalizing frequently-used long arcs — under a time budget. It is the strong reference against which heuristics are judged.
from ortools.constraint_solver import pywrapcp, routing_enums_pb2
def solve_ortools(inst, time_limit=8):
n, d, q, Q, K = inst["n"], inst["dist"], inst["demand"], inst["cap"], inst["K"]
D = (d*100).astype(int)
mgr = pywrapcp.RoutingIndexManager(n, K, 0); rm = pywrapcp.RoutingModel(mgr)
ti = rm.RegisterTransitCallback(lambda a,b: int(D[mgr.IndexToNode(a)][mgr.IndexToNode(b)]))
rm.SetArcCostEvaluatorOfAllVehicles(ti)
di = rm.RegisterUnaryTransitCallback(lambda a: int(q[mgr.IndexToNode(a)]))
rm.AddDimensionWithVehicleCapacity(di, 0, [Q]*K, True, "Cap")
pr = pywrapcp.DefaultRoutingSearchParameters()
pr.first_solution_strategy = routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC
pr.local_search_metaheuristic = routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH
pr.time_limit.seconds = time_limit
t=time.time(); sol = rm.SolveWithParameters(pr); t=time.time()-t
if sol is None: return None, np.inf, t
routes=[]
for v in range(K):
idx=rm.Start(v); r=[]
while not rm.IsEnd(idx):
nd=mgr.IndexToNode(idx)
if nd!=0: r.append(nd)
idx=sol.Value(rm.NextVar(idx))
if r: routes.append(r)
return routes, route_distance(routes, d), t
r_or, d_or, t_or = solve_ortools(inst_L, time_limit=8)
fig, ax = plt.subplots(1,2, figsize=(13, 6.5)); plot_routes(inst_L, r_hL, f"Savings+2-opt — {d_hL:.0f}", ax[0])
plot_routes(inst_L, r_or, f"OR-Tools (GLS) — {d_or:.0f}", ax[1]); plt.show()
print(f"OR-Tools: {d_or:.0f} ({t_or:.1f}s) vs heuristic {d_hL:.0f} ({t_hL*1000:.0f} ms) -> OR-Tools {(1-d_or/d_hL)*100:.1f}% shorter")
OR-Tools: 1165 (8.0s) vs heuristic 1178 (4 ms) -> OR-Tools 1.1% shorter
5 · Learning-guided optimization — predict-then-optimize for the search¶
Here is the ML×OR frontier. In an optimal route, edges are overwhelmingly short and local — a vehicle almost never jumps across the map. So most of the $n^2$ possible arcs are irrelevant, and if we could predict which arcs a good solution will use, we could shrink the problem before optimizing it.
We do exactly that: from a batch of instances solved by OR-Tools, we learn an edge classifier — features per ordered pair $(i,j)$ (distance, whether $j$ is among $i$'s nearest neighbours, mutual nearness, distances to the depot) → label "this arc is in the solution." On a fresh instance the model scores every arc; we keep each node's top-$k$ candidate arcs and re-solve on that sparse graph (savings + 2-opt restricted to candidate edges). This is predict-then-optimize applied to the search structure itself — the same theme as the inventory and staffing families, now sparsifying a combinatorial problem instead of stocking to a quantile.
from sklearn.ensemble import GradientBoostingClassifier
def edge_features(inst, k=10):
d, q, n = inst["dist"], inst["demand"], inst["n"]
order = np.argsort(d, axis=1) # nearest-neighbour ranks
rank = np.zeros((n,n));
for i in range(n): rank[i, order[i]] = np.arange(n)
feats, pairs = [], []
for i in range(n):
for j in range(n):
if i==j: continue
feats.append([d[i,j], rank[i,j], rank[j,i], d[0,i], d[0,j], q[j]]); pairs.append((i,j))
return np.array(feats), pairs
def sol_edges(routes):
E=set()
for r in routes:
path=[0]+list(r)+[0]
for a in range(len(path)-1): E.add((path[a], path[a+1]))
return E
# --- training set: solve small instances with OR-Tools, label their solution edges ---
X, y = [], []
for s in range(25):
ins = make_instance(45, seed=100+s, cap=40)
rr, _, _ = solve_ortools(ins, time_limit=2)
E = sol_edges(rr); F, P = edge_features(ins)
X.append(F); y.append(np.array([1 if p in E else 0 for p in P]))
X = np.vstack(X); y = np.concatenate(y)
clf = GradientBoostingClassifier(max_depth=3, n_estimators=150).fit(X, y)
print(f"edge classifier trained on {len(y):,} arcs ({y.mean()*100:.1f}% are solution arcs) | train AUC-ish acc {clf.score(X,y):.3f}")
# --- score every arc on the test instance ---
Ft, Pt = edge_features(inst_L); prob = clf.predict_proba(Ft)[:,1]
n = inst_L["n"]; score = np.zeros((n,n))
for (i,j), pr in zip(Pt, prob): score[i,j] = pr
E_or = sol_edges(r_or)
print(f"scored all {n*(n-1):,} arcs of the 100-customer test instance; next we keep each node's best few")
edge classifier trained on 51,750 arcs (2.5% are solution arcs) | train AUC-ish acc 0.977 scored all 10,100 arcs of the 100-customer test instance; next we keep each node's best few
# --- solve on the learned candidate graph: savings + 2-opt restricted to candidate edges ---
def clarke_wright_cand(inst, cand):
d, q, Q, n = inst["dist"], inst["demand"], inst["cap"], inst["n"]
routes={i:[i] for i in range(1,n)}; rof={i:i for i in range(1,n)}; load={i:q[i] for i in range(1,n)}
sav=sorted(((d[0,i]+d[0,j]-d[i,j],i,j) for i in range(1,n) for j in range(i+1,n) if cand[i,j]), reverse=True)
for s,i,j in sav:
ri,rj=rof[i],rof[j]
if ri==rj or load[ri]+load[rj]>Q: continue
Ri,Rj=routes[ri],routes[rj]
if Ri[-1]==i and Rj[0]==j: new=Ri+Rj
elif Ri[0]==i and Rj[-1]==j: new=Rj+Ri
elif Ri[-1]==i and Rj[-1]==j: new=Ri+Rj[::-1]
elif Ri[0]==i and Rj[0]==j: new=Ri[::-1]+Rj
else: continue
for k in new: rof[k]=ri
routes[ri]=new; load[ri]+=load[rj]; del routes[rj]
return list(routes.values())
def build_cand(k):
cand = np.zeros((n,n), bool)
for i in range(n):
for j in np.argsort(score[i])[::-1][:k]:
if i!=j: cand[i,j]=cand[j,i]=True
return cand
def cand_solve(cand):
r = [two_opt(rr, inst_L["dist"]) for rr in clarke_wright_cand(inst_L, cand)]
return r, route_distance(r, inst_L["dist"])
n_full = n*(n-1)//2
sweep = []
for k in [5, 8, 12, 16, 20]:
cand = build_cand(k)
cov = np.mean([1 if (cand[i,j] or i==0 or j==0) else 0 for (i,j) in E_or])
_, d = cand_solve(cand)
sweep.append((k, int(cand.sum()//2), cov, d))
Sw = pd.DataFrame(sweep, columns=["k","edges","coverage","distance"]); Sw["pct_pairs"] = Sw.edges/n_full*100
print(Sw.round(3).to_string(index=False)); print(f"full-graph heuristic: {d_hL:.0f} over {n_full} pairs")
t=time.time(); cand12 = build_cand(12); r_ml, d_ml = cand_solve(cand12); t_ml=time.time()-t; n_cand = int(cand12.sum()//2)
cov12 = np.mean([1 if (cand12[i,j] or i==0 or j==0) else 0 for (i,j) in E_or])
fig, ax = plt.subplots(1, 2, figsize=(13, 4.8))
ax[0].plot(Sw.pct_pairs, Sw.distance, "o-", color=PURP, label="candidate heuristic")
ax[0].axhline(d_hL, color=GREY, ls="--", label=f"full heuristic ({d_hL:.0f})")
ax0b = ax[0].twinx(); ax0b.plot(Sw.pct_pairs, Sw.coverage*100, "s--", color=GREEN)
ax[0].set_xlabel("% of arc-pairs kept"); ax[0].set_ylabel("total distance"); ax0b.set_ylabel("coverage of optimal arcs (%)", color=GREEN)
ax[0].set_title("Learned sparsification: quality vs search size"); ax[0].legend(fontsize=8, loc="upper right")
pts = inst_L["pts"]
for i in range(n):
for j in range(i+1, n):
if cand12[i,j]: ax[1].plot(pts[[i,j],0], pts[[i,j],1], "-", color=BLUE, lw=.3, alpha=.3)
plot_routes(inst_L, r_ml, f"Solution on candidate graph (k=12, {n_cand} edges) — {d_ml:.0f}", ax[1]); plt.show()
print(f"at k=12: {n_cand} edges ({n_cand/n_full*100:.0f}% of pairs), coverage {cov12*100:.0f}%, distance {d_ml:.0f} vs full {d_hL:.0f}")
k edges coverage distance pct_pairs 5 318 0.795 1223.107 6.297 8 505 0.884 1198.883 10.000 12 730 0.920 1163.742 14.455 16 954 0.946 1165.923 18.891 20 1138 0.955 1158.317 22.535 full-graph heuristic: 1178 over 5050 pairs
at k=12: 730 edges (14% of pairs), coverage 92%, distance 1164 vs full 1178
6 · The ladder, side by side¶
rows = [("Exact MIP", f"{inst_s['n']-1} cust", round(dist_x,1), round(t_x,2), "optimal, doesn't scale"),
("Savings + 2-opt", "100 cust", round(d_hL,0), round(t_hL,3), "instant, decent"),
("OR-Tools GLS", "100 cust", round(d_or,0), round(t_or,2), "production reference"),
("ML-guided (candidate)", "100 cust", round(d_ml,0), round(t_ml,3), f"only {n_cand/n_full*100:.0f}% of pairs searched")]
summary = pd.DataFrame(rows, columns=["method","instance","distance","seconds","note"])
print(summary.to_string(index=False))
fig, ax = plt.subplots(figsize=(13, 4.5))
m100 = summary[summary.instance=="100 cust"]
ax.bar(m100["method"], m100["distance"], color=[GREEN,BLUE,PURP]); ax.set_ylabel("total distance (100 customers)")
ax.set_title("Solution quality on the 100-customer instance"); ax.tick_params(axis="x", rotation=12)
for i,v in enumerate(m100["distance"]): ax.text(i,v,f"{v:.0f}",ha="center",va="bottom",fontsize=9)
ax.set_ylim(m100["distance"].min()*0.9, m100["distance"].max()*1.05); fig.tight_layout(); plt.show()
method instance distance seconds note
Exact MIP 9 cust 217.2 1.150 optimal, doesn't scale
Savings + 2-opt 100 cust 1178.0 0.004 instant, decent
OR-Tools GLS 100 cust 1165.0 8.000 production reference
ML-guided (candidate) 100 cust 1164.0 0.003 only 14% of pairs searched
7 · Takeaways¶
- Exact optimization gives a provable optimum but is NP-hard — usable only for tens of stops; it anchors "how good is good" on small instances.
- A simple heuristic (savings + 2-opt) delivers a solid solution in milliseconds and scales to hundreds of stops — the workhorse when speed matters.
- A production solver (OR-Tools guided local search) tightens that further and is the realistic reference for last-mile operations.
- Learning-guided optimization shows the ML×OR direction: a model trained on solved instances learns what a good edge looks like, letting us prune the problem to each node's most promising arcs. There is a clean quality-vs-size tradeoff — a larger candidate set covers more of the optimal arcs and recovers full solution quality, while still searching only a fraction of all pairs. On a 100-stop instance the raw speed gain is moot (the heuristic is already instant), but this learned sparsification is exactly how modern large-scale and neural routing methods make otherwise-intractable instances solvable. It is predict-then-optimize aimed at the structure of the search, complementing the demand/uncertainty framing of the other families.
Fourth family of ML in Operations Research: alongside Predict-then-Optimize, Stochastic Allocation, and Queueing, routing rounds out the section with the combinatorial side of OR — and the same lesson that ML is most powerful not as a replacement for optimization but as a way to guide and accelerate it.