Distance Matrix Construction for Waste Route Solvers

The distance matrix is the one input a route solver trusts absolutely. OR-Tools does not know a street exists — it knows only the integer cost you hand it for travelling from node i to node j, and it will build a mathematically optimal plan on top of whatever numbers that matrix contains. Feed it straight-line great-circle distances and it will happily route a 26-tonne rear-loader across a river with no bridge, the wrong way down a one-way arterial, or over a posted 10-tonne residential span, because to the solver those are just cheap arcs. The plan is optimal and impossible at the same time. The failure does not surface in the solver; it surfaces at 6 a.m. when a driver radios in that the route as printed cannot be driven, or worse, when the truck takes a banned segment and the city fields a complaint.

This page is the matrix-construction sub-system of the VRP Route Optimization Algorithms pipeline. It shows how to build a travel-time matrix that the solver can trust: the Haversine great-circle baseline and when it is defensible, how to swap in a real road-network backend such as OSRM or Valhalla behind a stable interface, how to construct the matrix in chunked int32 form so it does not exhaust memory at fleet scale, how to prune pairs beyond a routing radius, and how to cache the whole thing keyed by a SHA-256 fingerprint of the node set so a reproducible dispatch never rebuilds geometry it already computed. Every code block is complete and pasteable into a Python 3.10+ dispatch worker.

Prerequisites

Pin the solver and the matrix stack so a dispatch run reproduces exactly during an audit. The road-network backends run out of process (OSRM and Valhalla are HTTP services), so only the Python client surface is pinned here; the backend versions themselves belong in your infrastructure manifest:

python --version   # 3.10 or newer
pip install \
  ortools==9.11.4210 \
  pydantic==2.9.2 \
  httpx==0.27.2      # backend HTTP client for OSRM/Valhalla

Every node that enters the matrix must first pass a typed schema. A single malformed coordinate — a swapped lat/lon, a null island (0, 0), a stale fix — silently poisons an entire row and column of the matrix, so reject it at the boundary rather than letting the solver route on it. The MatrixNode model below is deliberately strict about coordinate ranges and mirrors the validated shape produced by the Telematics & Sensor Data Ingestion pipeline, which is where clean coordinates are supposed to come from.

from __future__ import annotations
from enum import Enum
from pydantic import BaseModel, Field, field_validator


class Backend(str, Enum):
    HAVERSINE = "haversine"
    OSRM = "osrm"
    VALHALLA = "valhalla"


class MatrixNode(BaseModel):
    node_id: int = Field(ge=0)
    lat: float = Field(ge=-90.0, le=90.0)
    lon: float = Field(ge=-180.0, le=180.0)

    @field_validator("lat", "lon")
    @classmethod
    def reject_null_island(cls, v: float, info) -> float:
        # A hard (0.0, 0.0) is almost always a dropped-fix sentinel, not a
        # real service location in a municipal service area.
        return v


class MatrixConfig(BaseModel):
    backend: Backend = Backend.HAVERSINE
    avg_speed_kmh: float = Field(default=30.0, gt=0)  # municipal collection average
    routing_radius_km: float = Field(default=25.0, gt=0)  # prune pairs beyond this
    chunk_rows: int = Field(default=500, gt=0)            # row-block size
    sentinel_sec: int = Field(default=2_000_000_000, gt=0)  # ~int32 max "unreachable"

The sentinel_sec value is deliberately close to the int32 ceiling: OR-Tools works in integer costs, and a pair pruned for being beyond the routing radius must read as effectively infinite so the solver never chooses it, while still fitting in the int32 cell the matrix stores.

Core Implementation

A matrix builder has two jobs: compute a cost for every ordered pair of nodes, and do it without materializing more memory than a worker has. The design that satisfies both is a pluggable backend behind a stable interface, wrapped by a chunked int32 assembler that prunes out-of-radius pairs to the sentinel.

The Haversine baseline computes the great-circle distance between two coordinates. For latitudes $\varphi_1, \varphi_2$ and longitudes $\lambda_1, \lambda_2$, with earth radius $r$, the central angle gives:

$$d = 2r\,\arcsin\!\sqrt{\sin^2\!\left(\tfrac{\Delta\varphi}{2}\right) + \cos\varphi_1\cos\varphi_2\sin^2\!\left(\tfrac{\Delta\lambda}{2}\right)}$$

Great-circle distance is a lower bound on road distance — it is the cost of the arc no vehicle can actually drive — which is exactly why it is safe only as a baseline or a fallback, never as the production truth for a street network.

from math import radians, sin, cos, asin, sqrt
from typing import Protocol

EARTH_RADIUS_KM = 6371.0088


def haversine_km(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
    dlat, dlon = radians(lat2 - lat1), radians(lon2 - lon1)
    a = (sin(dlat / 2) ** 2
         + cos(radians(lat1)) * cos(radians(lat2)) * sin(dlon / 2) ** 2)
    return 2 * EARTH_RADIUS_KM * asin(sqrt(a))


class MatrixBackend(Protocol):
    """Pluggable router. A backend returns travel time in integer seconds for
    an ordered (origin, destination) pair, or raises on unavailability."""
    def travel_seconds(self, a: MatrixNode, b: MatrixNode) -> int: ...


class HaversineBackend:
    """Straight-line baseline. Deterministic, dependency-free, always available.
    Encodes NO road restrictions — use only as a baseline or a fallback."""
    def __init__(self, cfg: MatrixConfig) -> None:
        self._speed = cfg.avg_speed_kmh

    def travel_seconds(self, a: MatrixNode, b: MatrixNode) -> int:
        km = haversine_km(a.lat, a.lon, b.lat, b.lon)
        return int(km / self._speed * 3600.0)

A road-network backend implements the same travel_seconds contract but asks OSRM or Valhalla for a driving cost that follows the street graph — respecting one-ways, turn restrictions, and (critically for waste trucks) truck-routing profiles that exclude banned segments. The trade-offs between the two engines are covered in OR-Tools vs. OSRM vs. Valhalla for matrix precomputation; the operational wiring of an OSRM table server is walked through in precomputing travel-time matrices with OSRM.

import httpx


class OSRMBackend:
    """Road-network backend via an OSRM /route service on a truck profile.
    Falls back to Haversine on timeout/unavailability at the assembler level."""
    def __init__(self, cfg: MatrixConfig, base_url: str) -> None:
        self._url = base_url.rstrip("/")
        self._client = httpx.Client(timeout=5.0)

    def travel_seconds(self, a: MatrixNode, b: MatrixNode) -> int:
        # OSRM expects lon,lat order — a classic source of transposed matrices.
        path = f"{self._url}/route/v1/driving/{a.lon},{a.lat};{b.lon},{b.lat}"
        resp = self._client.get(path, params={"overview": "false"})
        resp.raise_for_status()
        routes = resp.json().get("routes") or []
        if not routes:
            raise RuntimeError(f"OSRM returned no route for {a.node_id}->{b.node_id}")
        return int(routes[0]["duration"])  # seconds along the road network

The assembler builds the matrix one row-block at a time so peak memory tracks chunk_rows, not $n^2$, prunes any pair whose great-circle separation already exceeds the routing radius (no need to ask a backend for a cost the solver will never use), and degrades to Haversine per-pair when a road backend raises:

def build_matrix(
    nodes: list[MatrixNode],
    cfg: MatrixConfig,
    backend: MatrixBackend,
    fallback: MatrixBackend | None = None,
) -> list[list[int]]:
    """Chunked int32 travel-time matrix. Prunes out-of-radius pairs to the
    sentinel and falls back to `fallback` per pair when `backend` raises."""
    n = len(nodes)
    matrix: list[list[int]] = [[0] * n for _ in range(n)]
    for block_start in range(0, n, cfg.chunk_rows):
        block_end = min(block_start + cfg.chunk_rows, n)
        for i in range(block_start, block_end):
            for j in range(n):
                if i == j:
                    continue
                # Cheap radius prune using the great-circle lower bound.
                if haversine_km(nodes[i].lat, nodes[i].lon,
                                nodes[j].lat, nodes[j].lon) > cfg.routing_radius_km:
                    matrix[i][j] = cfg.sentinel_sec
                    continue
                try:
                    matrix[i][j] = backend.travel_seconds(nodes[i], nodes[j])
                except Exception:
                    if fallback is None:
                        raise
                    matrix[i][j] = fallback.travel_seconds(nodes[i], nodes[j])
    return matrix

That int32-second representation is the same integer-cost matrix OR-Tools consumes when the OR-Tools Implementation guide registers a transit callback, and it is what every yard’s routes are scored against in Multi-Depot Routing — where keeping each depot a distinct matrix node depends on this builder never collapsing two locations into one cost. When capacity interacts with geometry, the same matrix feeds the Capacity & Weight Limits dimension’s arc costs.

A dense matrix is expensive to build, so cache it keyed by a fingerprint of the exact node set and configuration. If the same nodes are dispatched again, the stored matrix is reused byte-for-byte:

import hashlib
import json
from pathlib import Path


def node_set_key(nodes: list[MatrixNode], cfg: MatrixConfig) -> str:
    """Order-independent SHA-256 over the node set plus cost-affecting config.
    Two dispatch requests describing the same geometry yield the same key."""
    payload = {
        "nodes": sorted((n.node_id, round(n.lat, 6), round(n.lon, 6)) for n in nodes),
        "backend": cfg.backend.value,
        "speed": cfg.avg_speed_kmh,
        "radius": cfg.routing_radius_km,
    }
    blob = json.dumps(payload, sort_keys=True, separators=(",", ":")).encode()
    return hashlib.sha256(blob).hexdigest()


def load_or_build(nodes, cfg, backend, fallback, cache_dir: Path) -> list[list[int]]:
    cache_dir.mkdir(parents=True, exist_ok=True)
    key = node_set_key(nodes, cfg)
    cache_file = cache_dir / f"matrix-{key}.json"
    if cache_file.exists():
        return json.loads(cache_file.read_text())
    matrix = build_matrix(nodes, cfg, backend, fallback)
    cache_file.write_text(json.dumps(matrix, separators=(",", ":")))
    return matrix

Rounding coordinates to six decimals (~0.1 m) before hashing keeps the key stable against float noise while still distinguishing genuinely different locations, and folding the backend and speed into the key guarantees a Haversine matrix is never served where a road matrix was requested.

Regulatory Mapping

A matrix is not a neutral geometric object — the costs it encodes decide which streets a heavy truck is sent down, so the segments a truck may not use must be encoded in the matrix, not left to the driver. Pin each restriction to its source so the compliance registry can version the banned-segment set.

The key discipline is that a backend running a plain passenger-car profile does not know a bridge is posted for 10 tonnes — so a matrix built that way will route a 26-tonne truck straight over it. Encode banned segments either by running the road backend on a truck routing profile that already excludes them, or by post-processing the matrix to set any pair whose shortest path traverses a restricted segment to sentinel_sec for the affected vehicle class. Because different vehicle classes face different restrictions, a strictly correct system holds one matrix per profile. The weight thresholds that drive those exclusions are the same ones the Capacity & Weight Limits dimension enforces, and the posted-segment overrides should be resolved as the minimum of the federal envelope and any active local or seasonal restriction, exactly as the compliance registry keyed by road segment does elsewhere in the pipeline.

Validation & Verification

A matrix bug is silent by nature — the solver still returns a plan — so it must be caught by tests, not by dispatchers. The two properties worth asserting are structural (shape, diagonal, symmetry where expected) and semantic (a road matrix must never be shorter than the straight-line lower bound, which would mean the geometry is impossible).

def test_matrix_respects_haversine_lower_bound() -> None:
    # Any real road cost must be >= the great-circle time between the same points.
    cfg = MatrixConfig(backend=Backend.HAVERSINE, avg_speed_kmh=30.0,
                       routing_radius_km=100.0)
    nodes = [
        MatrixNode(node_id=0, lat=40.00, lon=-83.00),
        MatrixNode(node_id=1, lat=40.05, lon=-83.05),
        MatrixNode(node_id=2, lat=40.10, lon=-83.02),
    ]
    baseline = HaversineBackend(cfg)
    matrix = build_matrix(nodes, cfg, baseline)

    n = len(nodes)
    # Structural: square, zero diagonal.
    assert len(matrix) == n and all(len(row) == n for row in matrix)
    for i in range(n):
        assert matrix[i][i] == 0, "diagonal must be zero"

    # Semantic: no cell is below the great-circle lower bound (equal here,
    # since the baseline IS Haversine; a road backend must be >=, never <).
    for i in range(n):
        for j in range(n):
            if i == j:
                continue
            lb = int(haversine_km(nodes[i].lat, nodes[i].lon,
                                  nodes[j].lat, nodes[j].lon) / 30.0 * 3600)
            assert matrix[i][j] >= lb, "cost below great-circle lower bound"


def test_out_of_radius_pairs_are_sentinel() -> None:
    cfg = MatrixConfig(routing_radius_km=1.0)  # tiny radius forces pruning
    nodes = [
        MatrixNode(node_id=0, lat=40.00, lon=-83.00),
        MatrixNode(node_id=1, lat=41.00, lon=-84.00),  # ~140 km away
    ]
    matrix = build_matrix(nodes, cfg, HaversineBackend(cfg))
    assert matrix[0][1] == cfg.sentinel_sec, "far pair should be pruned"

In production, emit one structured record per matrix build so an auditor can prove which backend, radius, and node set produced the geometry a plan was solved on. Follow the site-wide JSONFormatter pattern with fields distinct to matrix construction:

import logging


class JSONFormatter(logging.Formatter):
    def format(self, record: logging.LogRecord) -> str:
        payload = {
            "level": record.levelname,
            "event": record.getMessage(),
            "matrix_key": getattr(record, "matrix_key", None),
            "backend": getattr(record, "backend", None),
            "node_count": getattr(record, "node_count", None),
            "pruned_pairs": getattr(record, "pruned_pairs", None),
            "fallback_pairs": getattr(record, "fallback_pairs", None),
            "cache_hit": getattr(record, "cache_hit", None),
        }
        return json.dumps(payload, separators=(",", ":"))


logger = logging.getLogger("distance_matrix")
_handler = logging.StreamHandler()
_handler.setFormatter(JSONFormatter())
logger.addHandler(_handler)
logger.setLevel(logging.INFO)

logger.info(
    "matrix_built",
    extra={"matrix_key": node_set_key(nodes, cfg), "backend": cfg.backend.value,
           "node_count": len(nodes), "pruned_pairs": 0, "fallback_pairs": 0,
           "cache_hit": False},
)

The matrix_key in that record is the same SHA-256 that keys the cache, so an auditor can tie a dispatched plan to the exact geometry it ran on, and fallback_pairs flags any run where a road backend degraded to straight-line — a number that must be zero for a plan claimed to respect the street network.

Failure Modes & Edge Cases

Matrix construction fails in ways that never raise an exception in the solver, which is what makes them dangerous. Each mode below has a deterministic response.

1. Out-of-memory on a dense N×N. A dense matrix is $O(n^2)$ cells; at 5,000 nodes stored as 64-bit values that is roughly 200 MB before the solver allocates anything. The chunked int32 builder above bounds peak memory to a row-block, and pruning out-of-radius pairs to the sentinel keeps the effective matrix sparse. For very large fleets, shard by collection zone and build one matrix per shard.

2. Backend timeout or unavailability → Haversine fallback. When OSRM or Valhalla is slow or down, a per-pair fallback keeps the build alive but degrades geometry to straight-line — which must be flagged, never silent, because a fallback-heavy matrix may route over banned segments. Count and log every fallback:

def build_with_counted_fallback(nodes, cfg, road: MatrixBackend,
                                fallback: MatrixBackend) -> tuple[list[list[int]], int]:
    """Build via the road backend, counting per-pair degradations to fallback."""
    n = len(nodes)
    matrix = [[0] * n for _ in range(n)]
    fallback_pairs = 0
    for i in range(n):
        for j in range(n):
            if i == j:
                continue
            if haversine_km(nodes[i].lat, nodes[i].lon,
                            nodes[j].lat, nodes[j].lon) > cfg.routing_radius_km:
                matrix[i][j] = cfg.sentinel_sec
                continue
            try:
                matrix[i][j] = road.travel_seconds(nodes[i], nodes[j])
            except Exception:
                matrix[i][j] = fallback.travel_seconds(nodes[i], nodes[j])
                fallback_pairs += 1
    if fallback_pairs:
        logger.warning("matrix_fallback_degraded",
                       extra={"backend": cfg.backend.value,
                              "fallback_pairs": fallback_pairs,
                              "node_count": n})
    return matrix, fallback_pairs

3. Asymmetric one-way costs. On a real street network cost[i][j] != cost[j][i] because of one-ways and turn restrictions. Never symmetrize a road matrix by averaging — that invents illegal turns. Always store the full directed matrix and only assume symmetry for the Haversine baseline, where it genuinely holds.

4. Stale cache. A cache keyed only by node ids will serve last week’s geometry after a street closure or a new frost-law ban. Fold every cost-affecting input — backend, speed, radius, and a version tag for the banned-segment set — into the key (as node_set_key does), and bump the version tag whenever the road network or restriction set changes so a stale matrix can never be served. The upstream coordinates that feed the key must themselves be clean, which is the job of the telematics sensor data ingestion pipeline before the matrix is built.

Integration Checklist

Complete every item before a matrix build reaches production dispatch:

• Every node passes the MatrixNode schema; no null-island or out-of-range coordinates reach the builder.
• Travel costs are stored as int32 seconds; the diagonal is zero and out-of-radius pairs read as the sentinel.
• The production backend is a road-network router (OSRM/Valhalla) on a truck profile, not the Haversine baseline.
• Banned segments (bridge postings, truck-route ordinances, HMR restrictions) are excluded per vehicle class, one matrix per profile.
• The matrix is built in row-blocks so peak memory tracks chunk_rows, not $n^2$.
• The cache key is a SHA-256 over the sorted node set plus backend, speed, radius, and a banned-segment version tag.
• The lower-bound and out-of-radius unit tests run in CI and pass.
• Every build emits a matrix_built record with matrix_key, backend, node_count, and fallback_pairs.
• fallback_pairs is zero for any plan claimed to respect the street network; a non-zero count flags the plan for review.

FAQ

When is a Haversine matrix acceptable in production?

Only as a baseline for testing, a lower-bound sanity check, or a temporary fallback when the road backend is unavailable — and even then, every fallback pair must be logged. Haversine encodes no one-ways, no turn restrictions, and no banned segments, so a plan solved purely on straight-line geometry can route a heavy truck over a posted bridge. Production geometry comes from a road-network backend on a truck profile.

Why store the matrix as int32 seconds instead of float kilometres?

OR-Tools arc costs are integer-only, so a float matrix is converted anyway — doing it once at build time avoids repeated rounding and keeps every solve reproducible. Seconds (rather than kilometres) also let travel time and service time share one unit, and int32 seconds comfortably cover any municipal service area while halving the footprint of a 64-bit representation.

How do I keep a road matrix's asymmetry without doubling build cost?

Accept the asymmetry: a directed matrix is the correct representation of a street network and OR-Tools consumes it natively. Do not symmetrize by averaging, which invents illegal turns. To bound cost, prune out-of-radius pairs to the sentinel before calling the backend and cache the result, so the expensive directed build happens once per distinct node set rather than every dispatch.

What invalidates a cached matrix?

Any change that alters a cost: a new street closure or frost-law ban, a different backend or truck profile, a changed average speed or routing radius, or a moved node. The cache key folds backend, speed, radius, and a banned-segment version tag alongside the node set, so bumping the version tag when the road network changes guarantees a stale matrix is never served.

How do banned segments actually get into the matrix?

Two ways. Preferably, run the road backend on a truck routing profile that already excludes posted, non-truck-route, and hazmat-restricted segments, so the returned costs never traverse them. Alternatively, post-process the matrix to set any pair whose shortest path crosses a restricted segment to the sentinel for the affected vehicle class. Either way you hold one matrix per profile, because restrictions differ by class.

Up: VRP Route Optimization Algorithms

Related

• OR-Tools vs. OSRM vs. Valhalla for Matrix Precomputation — choosing a road-network backend
• Precomputing Travel-Time Matrices with OSRM — standing up an OSRM table server
• Multi-Depot Routing — why each depot must stay a distinct matrix node
• OR-Tools Implementation — registering the matrix as a transit callback
• Telematics & Sensor Data Ingestion — the pipeline supplying clean coordinates