Combinatorial Optimization for Fleets

The Capacitated Vehicle Routing Problem (CVRP) is a combinatorial constraint, not a geometric one. 10 nodes produce ~3.6 million possible routes. 50 nodes produce ~10⁶⁴. The number of permutations grows factorially with N, which puts the problem in NP-Hard territory, no known algorithm guarantees optimality in polynomial time.

Vehicle routing plan for a taxi fleet serving distributed demand

The standard formulation minimizes total travel cost across a fleet of capacity-constrained vehicles. Given a depot node $0$ and customer nodes $1\dots N$, each with demand $q_i$, and $K$ identical vehicles each with capacity $Q$:

$$\text{Minimize }\sum _{k=1}^K\sum _{i=0}^N\sum _{j=0}^N{d}_{ij}\cdot x_{ijk}$$

Subject to:

• Every customer is visited exactly once: ${\sum }_k{\sum }_j{x}_{ijk}=1$ for all $i>0$
• Flow conservation: ${\sum }_i{x}_{ipk}-{\sum }_j{x}_{pjk}=0$ for all $p,k$
• Capacity: ${\sum }_i{q}_i{\sum }_j{x}_{ijk}\le Q$ for all $k$

Where $x_{ijk}$ is 1 if vehicle $k$ traverses edge $(i,j)$, and 0 otherwise. The decision variables are binary, the constraints are linear, the cost depends on $d_{ij}$ – and that distance matrix is where the engineering starts.

You can compute $d_{ij}$ two ways. Haversine distance uses lat/lng directly and runs in O(N²) time – roughly 2ms for 100 nodes. It ignores roads, one-way streets, and traffic, but it works for any coordinate pair with zero API calls. The alternative is a routing engine like OSRM, which accounts for road network topology. It also runs $O(N^2)$, but each query is an HTTP request: at ~50ms per pair, 100 nodes takes over 8 minutes to build the matrix. Haversine is a prototyping tool; OSRM is a production cost.

With the cost matrix built, the solver layer handles the combinatorial search. OR-Tools provides a production-grade CVRP solver with exact and heuristic strategies built in:

from ortools.constraint_solver import routing_enums_pb2, pywrapcp
import math
import numpy as np
 
def solve_cvrp(locations, demands, num_vehicles, vehicle_capacity):
    n = len(locations)
    dist = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            if i != j:
                lat1, lon1 = locations[i]['lat'], locations[i]['lon']
                lat2, lon2 = locations[j]['lat'], locations[j]['lon']
                dphi = math.radians(lat2 - lat1)
                dlambda = math.radians(lon2 - lon1)
                a = math.sin(dphi/2)**2 + math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.sin(dlambda/2)**2
                dist[i][j] = 6371 * 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
 
    manager = pywrapcp.RoutingIndexManager(n, num_vehicles, 0)
    routing = pywrapcp.RoutingModel(manager)
 
    def distance_callback(from_index, to_index):
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return dist[from_node][to_node]
 
    transit_callback_index = routing.RegisterTransitCallback(distance_callback)
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
 
    def demand_callback(from_index):
        from_node = manager.IndexToNode(from_index)
        return demands[from_node]
 
    demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)
    routing.AddDimensionWithVehicleCapacity(
        demand_callback_index, 0, [vehicle_capacity] * num_vehicles, True, 'Capacity')
 
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = \
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARCH
    search_parameters.local_search_metaheuristic = \
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH
    search_parameters.time_limit.seconds = 10
 
    solution = routing.SolveWithParameters(search_parameters)
    return solution, routing, manager

The time limit of 10 seconds is deliberate, it forces the solver into the heuristic regime for any non-trivial instance. The boundary between exact and approximate depends on N and the number of vehicles. Benchmarks on standard CVRP instances (Augerat, Christofides) give a rough rule of thumb:

N	Vehicles	OR-Tools solve time	Regime
10	2–3	<1s	Exact (MIP/CP)
50	5–10	5–30s	Exact-to-heuristic boundary
100	10–20	30s–2m	Heuristic dominated
500	20–50	10s (limited)	Guided Local Search only
1000+	50+	10s (limited)	GLS, solution quality varies

Below ~50 nodes, you can reliably find provably optimal routes. Above that, the cost matrix quality (Haversine vs. road distance) dominates solution quality more than the solver strategy. Spending engineering effort on the distance model pays off before tuning solver hyperparameters does.