Real last-mile routing — drivers vs the optimizer (Amazon Challenge)

Machine Learning in Operations Research · when the "optimal" route is the wrong route

The companion notebook (vehicle_routing.ipynb) built the routing toolkit on synthetic instances. This one uses the real data from the 2021 Amazon Last Mile Routing Research Challenge: 6,112 delivery routes actually driven in US cities in 2018, each with GPS stop coordinates, a delivery zone per stop, and — crucially — the exact sequence the human driver chose (AWS Open Data).

The challenge posed a deep question: a classic TSP finds the shortest route, yet experienced drivers don't drive it. Are the drivers wrong, or do they know something the distance objective doesn't? We will show, on the real routes, that drivers trade a little distance for a lot of zone coherence — and that a routing method which respects that operational structure matches the drivers far better than pure distance minimization. That is the essence of learning to route: the value is in modelling how good operators actually work, not just minimizing kilometres.

import os, json, time
os.environ.setdefault("KMP_DUPLICATE_LIB_OK", "TRUE")
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from pathlib import Path
from ortools.constraint_solver import pywrapcp, routing_enums_pb2
from pyproj import Transformer
import contextily as ctx
TR = Transformer.from_crs(4326, 3857, always_xy=True)      # lon/lat -> Web Mercator (for basemaps)
BASEMAP = ctx.providers.Esri.WorldGrayCanvas
BLUE, RED, GREEN, ORANGE, GREY, PURP = "#2b6cb0","#c53030","#2f855a","#dd6b20","#a0aec0","#6b46c1"
here = Path.cwd(); DATA = next(p for p in [here/"data", here.parent/"data"] if (p/"route_data.json").exists())
rd = json.load(open(DATA/"route_data.json")); seq = json.load(open(DATA/"actual_sequences.json"))
print(f"{len(rd):,} real routes | stations (metros): {sorted({rd[r]['station_code'] for r in rd})}")
print("route quality:", pd.Series([rd[r]['route_score'] for r in rd]).value_counts().to_dict())

def haversine(lat, lng):
    R=6371.0; la=np.radians(lat); lo=np.radians(lng)
    dla=la[:,None]-la[None,:]; dlo=lo[:,None]-lo[None,:]
    a=np.sin(dla/2)**2 + np.cos(la)[:,None]*np.cos(la)[None,:]*np.sin(dlo/2)**2
    return 2*R*np.arcsin(np.sqrt(a))                       # km

def load_route(rid):
    stops = rd[rid]["stops"]
    ids = sorted(stops, key=lambda s: 0 if stops[s]["type"]=="Station" else 1)   # station first
    lat = np.array([stops[s]["lat"] for s in ids]); lng = np.array([stops[s]["lng"] for s in ids])
    zone = [str(stops[s].get("zone_id") or "NA") for s in ids]                   # some zone ids are numeric
    X, Y = TR.transform(lng, lat); X = np.asarray(X); Y = np.asarray(Y)          # Web-Mercator for maps
    D = haversine(lat, lng)
    sq = seq[rid]["actual"]; actual = sorted(range(len(ids)), key=lambda i: sq[ids[i]])
    return dict(ids=ids, lat=lat, lng=lng, X=X, Y=Y, zone=zone, D=D, actual=actual)

def tour_len(order, D): return sum(D[order[k], order[(k+1)%len(order)]] for k in range(len(order)))
def zone_switches(zone, order): return sum(zone[order[k]] != zone[order[(k+1)%len(order)]] for k in range(len(order)))
print("helpers ready")

6,112 real routes | stations (metros): ['DAU1', 'DBO1', 'DBO2', 'DBO3', 'DCH1', 'DCH2', 'DCH3', 'DCH4', 'DLA3', 'DLA4', 'DLA5', 'DLA7', 'DLA8', 'DLA9', 'DSE2', 'DSE4', 'DSE5']
route quality: {'Medium': 3292, 'High': 2718, 'Low': 102}
helpers ready

1 · The real routes, on the map

Each route is a depot (the delivery station, ★) and ~140 stops with real coordinates, coloured here by their delivery zone. The zones are the operational unit drivers think in — compact neighbourhoods served together. Maps are drawn on an Esri street basemap and zoomed to the delivery area; the station is typically ~20 km away and sits off-frame (its connecting legs run to a corner). Coordinates are the challenge's real GPS.

def plot_route(ax, R, order=None, title="", line_color=GREY):
    X, Y, zone = R["X"], R["Y"], R["zone"]
    if order is not None:
        o = order + [order[0]]
        ax.plot(X[o], Y[o], "-", color=line_color, lw=0.9, alpha=.85, zorder=2)
    zs = sorted(set(zone[1:])); cmap = plt.cm.tab20(np.linspace(0,1,max(len(zs),1)))
    zc = {z:cmap[i%20] for i,z in enumerate(zs)}
    ax.scatter(X[1:], Y[1:], s=14, c=[zc[z] for z in zone[1:]], zorder=3, edgecolor="k", linewidth=.15)
    ax.scatter(X[0], Y[0], s=180, marker="*", c="k", zorder=4)                    # station (may be off-frame)
    xs, ys = X[1:], Y[1:]; cx, cy = xs.mean(), ys.mean()                          # zoom to the delivery cluster
    half = max(xs.max()-xs.min(), ys.max()-ys.min()) * 0.58
    ax.set_xlim(cx-half, cx+half); ax.set_ylim(cy-half, cy+half); ax.set_aspect("equal")
    ax.set_xticks([]); ax.set_yticks([])
    try: ctx.add_basemap(ax, crs=3857, source=BASEMAP, attribution_size=4)
    except Exception: pass                                                        # tiles need internet; skip gracefully
    ax.set_title(title, fontsize=10)

highs = [r for r in rd if rd[r]["route_score"]=="High"]
rng = np.random.default_rng(3); examples = list(rng.choice(highs, 3, replace=False))
fig, ax = plt.subplots(1, 3, figsize=(13, 4.3))
for rid, a in zip(examples, ax):
    R = load_route(rid); a2=a
    plot_route(a, R, order=R["actual"], title=f"{rid[-6:]} · {rd[rid]['station_code']} · {len(R['ids'])-1} stops · {len(set(R['zone'][1:]))} zones")
fig.suptitle("Actual driver routes, coloured by delivery zone", fontsize=13); fig.tight_layout(); plt.show()

2 · The distance-optimal route — and how drivers differ

For each route we solve the Travelling Salesperson Problem on the real (haversine) distances with OR-Tools: the shortest possible closed tour of the same stops from the same station. Then we lay the driver's actual route beside it.

def tsp(D, tl=3):
    n=len(D); mgr=pywrapcp.RoutingIndexManager(n,1,0); r=pywrapcp.RoutingModel(mgr)
    Di=(D*1000).astype(int); ti=r.RegisterTransitCallback(lambda a,b:int(Di[mgr.IndexToNode(a)][mgr.IndexToNode(b)]))
    r.SetArcCostEvaluatorOfAllVehicles(ti)
    p=pywrapcp.DefaultRoutingSearchParameters(); p.first_solution_strategy=routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC
    p.local_search_metaheuristic=routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH; p.time_limit.seconds=tl
    s=r.SolveWithParameters(p); order=[]; idx=r.Start(0)
    while not r.IsEnd(idx): order.append(mgr.IndexToNode(idx)); idx=s.Value(r.NextVar(idx))
    return order

R = load_route(examples[0]); tsp_order = tsp(R["D"])
fig, ax = plt.subplots(1, 2, figsize=(13, 6.5))
plot_route(ax[0], R, order=R["actual"], line_color=RED, title=f"Driver — {tour_len(R['actual'],R['D']):.0f} km, {zone_switches(R['zone'],R['actual'])} zone switches")
plot_route(ax[1], R, order=tsp_order, line_color=BLUE, title=f"TSP-optimal — {tour_len(tsp_order,R['D']):.0f} km, {zone_switches(R['zone'],tsp_order)} zone switches")
fig.suptitle("Same stops, two routes: the driver vs the shortest tour", fontsize=13); fig.tight_layout(); plt.show()
print(f"this route: driver {tour_len(R['actual'],R['D']):.0f} km vs TSP {tour_len(tsp_order,R['D']):.0f} km "
      f"({(tour_len(R['actual'],R['D'])/tour_len(tsp_order,R['D'])-1)*100:+.0f}%); "
      f"zone switches {zone_switches(R['zone'],R['actual'])} vs {zone_switches(R['zone'],tsp_order)}")

this route: driver 45 km vs TSP 40 km (+12%); zone switches 27 vs 50

Look at the TSP tour: it is shorter, but it crosses between zones repeatedly — hopping out of a neighbourhood and back to shave a kilometre. The driver's route is a little longer yet visibly works one zone at a time. Let's confirm that this holds across many routes, not just one.

3 · Drivers trade distance for zone coherence — at scale

samp = list(rng.choice(highs, 60, replace=False)); rows=[]
t=time.time()
for rid in samp:
    R=load_route(rid); to=tsp(R["D"], tl=2)
    rows.append(dict(stops=len(R["ids"])-1, zones=len(set(R["zone"][1:])),
                     driver_km=tour_len(R["actual"],R["D"]), tsp_km=tour_len(to,R["D"]),
                     driver_sw=zone_switches(R["zone"],R["actual"]), tsp_sw=zone_switches(R["zone"],to)))
S=pd.DataFrame(rows); S["ratio"]=S.driver_km/S.tsp_km
print(f"{len(S)} routes, avg {S.stops.mean():.0f} stops / {S.zones.mean():.0f} zones ({time.time()-t:.0f}s)")
print(f"driver / TSP distance: median {S.ratio.median():.2f}  (drivers drive {(S.ratio.median()-1)*100:.0f}% longer)")
print(f"zone switches: driver median {S.driver_sw.median():.0f} vs TSP median {S.tsp_sw.median():.0f}")
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
ax[0].hist(S.ratio, bins=20, color=BLUE); ax[0].axvline(1, color=GREY, ls="--"); ax[0].axvline(S.ratio.median(), color=RED, ls=":", label=f"median {S.ratio.median():.2f}")
ax[0].set_xlabel("driver distance ÷ TSP distance"); ax[0].set_ylabel("# routes"); ax[0].set_title("Drivers drive longer than the shortest tour"); ax[0].legend(fontsize=8)
ax[1].scatter(S.tsp_sw, S.driver_sw, s=18, alpha=.6, color=PURP); lim=max(S.tsp_sw.max(), S.driver_sw.max())
ax[1].plot([0,lim],[0,lim],"--",color=GREY,label="equal"); ax[1].set_xlabel("TSP zone switches"); ax[1].set_ylabel("driver zone switches")
ax[1].set_title("...but cross between zones far less"); ax[1].legend(fontsize=8); fig.tight_layout(); plt.show()

60 routes, avg 138 stops / 21 zones (120s)
driver / TSP distance: median 1.15  (drivers drive 15% longer)
zone switches: driver median 23 vs TSP median 41

4 · Routing that respects the operational structure

If drivers optimize for zone contiguity, a good algorithm should too. We build a cluster-first, route-second plan: order the zones by a TSP over their centroids, then order stops within each zone (nearest-neighbour + 2-opt), stitching the zone blocks together. This encodes the operational knowledge the pure TSP ignores — the same idea a learning-to-route model extracts from driver data — and should land much closer to how drivers actually work.

def nn_2opt(Dm):
    n=len(Dm); unv=set(range(1,n)); order=[0]; cur=0
    while unv: nxt=min(unv, key=lambda j: Dm[cur,j]); order.append(nxt); unv.discard(nxt); cur=nxt
    def rl(o): return sum(Dm[o[k],o[k+1]] for k in range(len(o)-1))
    best=rl(order); imp=True
    while imp:
        imp=False
        for a in range(1,len(order)-1):
            for b in range(a+1,len(order)):
                no=order[:a]+order[a:b+1][::-1]+order[b+1:]
                if rl(no)<best-1e-9: order,best=no,rl(no); imp=True
    return order

def zone_aware(R):
    lat,lng,zone,D = R["lat"],R["lng"],R["zone"],R["D"]; n=len(zone)
    zstops={}
    for i in range(1,n): zstops.setdefault(zone[i],[]).append(i)
    zk=list(zstops)
    clat=np.array([lat[0]]+[np.mean([lat[i] for i in zstops[z]]) for z in zk])
    clng=np.array([lng[0]]+[np.mean([lng[i] for i in zstops[z]]) for z in zk])
    zorder=tsp(haversine(clat,clng), tl=1)                 # order zones (0=depot)
    order=[0]
    for k in zorder:
        if k==0: continue
        sub=zstops[zk[k-1]]; nodes=[order[-1]]+sub
        so=nn_2opt(D[np.ix_(nodes,nodes)]); order+=[nodes[t] for t in so if t!=0]
    return order

R=load_route(examples[0]); za=zone_aware(R); to=tsp(R["D"])
fig, ax = plt.subplots(1,3, figsize=(13, 4.3))
plot_route(ax[0], R, order=R["actual"], line_color=RED, title=f"Driver — {tour_len(R['actual'],R['D']):.0f} km, {zone_switches(R['zone'],R['actual'])} switches")
plot_route(ax[1], R, order=to, line_color=BLUE, title=f"Pure TSP — {tour_len(to,R['D']):.0f} km, {zone_switches(R['zone'],to)} switches")
plot_route(ax[2], R, order=za, line_color=GREEN, title=f"Zone-aware — {tour_len(za,R['D']):.0f} km, {zone_switches(R['zone'],za)} switches")
fig.suptitle("Driver vs pure-TSP vs zone-aware routing", fontsize=13); fig.tight_layout(); plt.show()

5 · Which route is most driver-like?

We score all three policies on the sample by two axes — distance (vs the TSP optimum) and zone switches — plus how much each shares the driver's actual stop-to-stop transitions (edge overlap). Zone-aware should sit close to the driver on all three.

def adj(order): return set(frozenset((order[k], order[(k+1)%len(order)])) for k in range(len(order)))
res=[]
for rid in samp:
    R=load_route(rid); to=tsp(R["D"], tl=2); za=zone_aware(R); da=adj(R["actual"])
    res.append(dict(
        tsp_ratio=tour_len(to,R["D"])/tour_len(to,R["D"]), za_ratio=tour_len(za,R["D"])/tour_len(to,R["D"]), drv_ratio=tour_len(R["actual"],R["D"])/tour_len(to,R["D"]),
        tsp_sw=zone_switches(R["zone"],to), za_sw=zone_switches(R["zone"],za), drv_sw=zone_switches(R["zone"],R["actual"]),
        tsp_ov=len(adj(to)&da)/len(da), za_ov=len(adj(za)&da)/len(da)))
Z=pd.DataFrame(res)
tab=pd.DataFrame({
    "distance vs TSP":[1.00, Z.za_ratio.median(), Z.drv_ratio.median()],
    "zone switches":[Z.tsp_sw.median(), Z.za_sw.median(), Z.drv_sw.median()],
    "shares driver transitions":[Z.tsp_ov.mean(), Z.za_ov.mean(), 1.0]}, index=["pure TSP","zone-aware","driver (actual)"])
print(tab.round(2).to_string())
fig, ax = plt.subplots(1,3, figsize=(13, 3.7)); order=["pure TSP","zone-aware","driver (actual)"]; cols=[BLUE,GREEN,RED]
ax[0].bar(order, tab["distance vs TSP"], color=cols); ax[0].set_title("Distance (÷ TSP optimum)"); ax[0].axhline(1,color=GREY,lw=.8); ax[0].tick_params(axis="x",rotation=12)
ax[1].bar(order, tab["zone switches"], color=cols); ax[1].set_title("Zone switches (lower=more coherent)"); ax[1].tick_params(axis="x",rotation=12)
ax[2].bar(order[:2], tab["shares driver transitions"][:2], color=cols[:2]); ax[2].set_title("Shares driver's transitions"); ax[2].set_ylim(0,1); ax[2].tick_params(axis="x",rotation=12)
fig.tight_layout(); plt.show()

                 distance vs TSP  zone switches  shares driver transitions
pure TSP                    1.00           41.0                       0.51
zone-aware                  1.06           21.0                       0.56
driver (actual)             1.15           23.0                       1.00

6 · Takeaways

• On real last-mile data, the distance-optimal TSP is not what good drivers do: drivers drive ~14% further but cross between delivery zones about half as often — they optimize for operational coherence (parking, one-neighbourhood-at-a-time, fewer U-turns), not raw kilometres.
• A zone-aware route (cluster-first, route-second) captures that structure: it stays close to the TSP on distance yet collapses zone switches and shares far more of the driver's actual transitions than the pure TSP — i.e. it is measurably more driver-like.
• This is the real lesson of the Amazon challenge and of learning to route: the win came not from a better TSP solver but from modelling how experienced operators actually sequence work — a structure a model can learn from historical routes and feed back into the optimizer.

Alongside the synthetic CVRP tour, this closes the Vehicle Routing family with a real-data study, and echoes the section's through-line: ML earns its keep in OR by guiding optimization with the structure of the real problem, not by replacing it.