Vehicle Routing Problem in Python — Deep Dive

The Vehicle Routing Problem in production is less about elegant math and more about encoding messy business rules into a solver that returns usable routes under time pressure. This guide covers full CVRPTW implementation, multi-depot extensions, and tuning for fleet sizes up to 1,000 stops.

Full CVRPTW implementation

from ortools.constraint_solver import routing_enums_pb2, pywrapcp

def solve_cvrptw(
    distance_matrix: list[list[int]],
    demands: list[int],
    time_matrix: list[list[int]],
    time_windows: list[tuple[int, int]],
    vehicle_capacities: list[int],
    num_vehicles: int,
    depot: int = 0,
    time_limit_seconds: int = 30,
):
    manager = pywrapcp.RoutingIndexManager(
        len(distance_matrix), num_vehicles, depot
    )
    routing = pywrapcp.RoutingModel(manager)

    # Distance callback
    def distance_cb(from_idx, to_idx):
        return distance_matrix[manager.IndexToNode(from_idx)][manager.IndexToNode(to_idx)]

    dist_cb_idx = routing.RegisterTransitCallback(distance_cb)
    routing.SetArcCostEvaluatorOfAllVehicles(dist_cb_idx)

    # Capacity dimension
    def demand_cb(idx):
        return demands[manager.IndexToNode(idx)]

    demand_cb_idx = routing.RegisterUnaryTransitCallback(demand_cb)
    routing.AddDimensionWithVehicleCapacity(
        demand_cb_idx, 0, vehicle_capacities, True, "Capacity"
    )

    # Time dimension
    def time_cb(from_idx, to_idx):
        return time_matrix[manager.IndexToNode(from_idx)][manager.IndexToNode(to_idx)]

    time_cb_idx = routing.RegisterTransitCallback(time_cb)
    routing.AddDimension(time_cb_idx, 60, 1440, False, "Time")
    time_dim = routing.GetDimensionOrDie("Time")

    for loc in range(len(time_windows)):
        index = manager.NodeToIndex(loc)
        time_dim.CumulVar(index).SetRange(*time_windows[loc])

    # Minimize waiting + transit (soft span cost)
    for v in range(num_vehicles):
        time_dim.SetSpanCostCoefficientForVehicle(1, v)

    search_params = pywrapcp.DefaultRoutingSearchParameters()
    search_params.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PARALLEL_CHEAPEST_INSERTION
    )
    search_params.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH
    )
    search_params.time_limit.FromSeconds(time_limit_seconds)

    solution = routing.SolveWithParameters(search_params)
    if not solution:
        return None

    routes = []
    for v in range(num_vehicles):
        route = []
        index = routing.Start(v)
        while not routing.IsEnd(index):
            node = manager.IndexToNode(index)
            t = solution.Value(time_dim.CumulVar(index))
            route.append({"stop": node, "arrival": t})
            index = solution.Value(routing.NextVar(index))
        route.append({"stop": manager.IndexToNode(index),
                       "arrival": solution.Value(time_dim.CumulVar(index))})
        routes.append(route)
    return routes, solution.ObjectiveValue()

Drop penalties for infeasible stops

When not every stop can fit within constraints, allow the solver to skip stops at a cost:

penalty = 10_000  # high enough to prefer serving, low enough to allow drops
for node in range(1, len(distance_matrix)):
    routing.AddDisjunction([manager.NodeToIndex(node)], penalty)

Dropped stops appear as unvisited in the solution. The business layer can then assign them to a next-day run or a partner carrier.

Multi-depot setup

OR-Tools supports per-vehicle start and end depots:

starts = [0, 0, 5, 5]  # vehicles 0-1 start at depot 0, vehicles 2-3 at depot 5
ends = [0, 0, 5, 5]
manager = pywrapcp.RoutingIndexManager(num_locations, num_vehicles, starts, ends)

This is critical for companies with regional warehouses. Each vehicle returns to its origin depot, and the solver avoids assigning stops far from a vehicle’s home base because the return leg would inflate cost.

Heterogeneous fleets

Real fleets mix vehicle types — vans, trucks, bikes. Model this with per-vehicle capacity and cost:

vehicle_capacities = [100, 100, 500, 500]  # vans vs trucks
vehicle_costs = [10, 10, 25, 25]           # cost per km

for v in range(num_vehicles):
    routing.SetFixedCostOfVehicle(vehicle_costs[v] * 100, v)  # fixed dispatch cost

The solver now prefers cheaper vehicles unless their capacity forces using a truck.

Solver strategy comparison

Benchmarked on a 300-stop, 15-vehicle CVRPTW instance (synthetic, uniform distribution):

Strategy	First solution	Metaheuristic	10 s result	30 s result
A	PATH_CHEAPEST_ARC	GUIDED_LOCAL_SEARCH	142,380	138,200
B	PARALLEL_CHEAPEST_INSERTION	GUIDED_LOCAL_SEARCH	139,100	136,800
C	SAVINGS	SIMULATED_ANNEALING	145,900	141,500
D	CHRISTOFIDES	TABU_SEARCH	143,200	139,900

Strategy B consistently wins for CVRPTW. For pure CVRP without time windows, SAVINGS + GUIDED_LOCAL_SEARCH often edges ahead.

Scaling beyond 1,000 stops

At scale, a single OR-Tools invocation struggles with both memory and solution quality. Decomposition strategies:

1. Geographic clustering — K-Means or DBSCAN splits stops into zones matching roughly one vehicle each. Solve each zone independently.
2. Time-based partitioning — morning deliveries and afternoon deliveries solved separately, sharing vehicles across shifts.
3. Hierarchical solve — first assign stops to depots (assignment problem), then solve VRP per depot.

Parallelizing with concurrent.futures.ProcessPoolExecutor across clusters yields near-linear speedup.

from concurrent.futures import ProcessPoolExecutor

def solve_cluster(cluster_data):
    return solve_cvrptw(**cluster_data)

with ProcessPoolExecutor(max_workers=4) as pool:
    results = list(pool.map(solve_cluster, cluster_configs))

Integration with dispatch systems

A production pipeline typically:

1. Pulls new orders from a database at a scheduling cutoff time (e.g., 5 AM).
2. Builds the distance matrix via OSRM or cached historical data.
3. Solves VRP with a 30-60 second budget.
4. Writes routes to a dispatch table that drivers access via a mobile app.
5. Monitors driver GPS and triggers re-optimization when delays exceed a threshold.

The re-optimization loop is the hardest part. Removing completed stops, inserting new orders, and respecting partially-driven routes requires careful state management.

Pitfalls from real deployments

• Service time omission — forgetting the 5 minutes per stop to park, walk, and deliver. Over 30 stops, that is 2.5 hours missing from the model.
• Symmetric matrix for asymmetric roads — one-way streets and highway ramps mean A→B ≠ B→A.
• Fixed vehicle count — letting the solver decide how many vehicles to use (via fixed costs) often outperforms forcing an exact number.
• Overfit to average day — seasonal demand spikes break static route templates. Re-solve at least weekly.

The one thing to remember: Production VRP combines capacity, time windows, and fleet heterogeneity into a single model that OR-Tools solves heuristically — and decomposition plus re-routing keep it viable at scale.