Balancing Facility Loads with Linear Assignment

Assigning every population unit to its nearest facility overloads popular sites and leaves others idle, a problem the broader facility capacity allocation models section addresses with constrained optimization. This guide solves the balanced version directly in Python: a capacitated assignment that minimizes total travel time subject to each facility’s capacity ceiling, using scipy.optimize.linear_sum_assignment with capacity replication.

Problem Context & Constraints

Nearest-facility assignment is seductive because it is trivial to compute — for each demand unit, pick the column of the travel-time matrix with the smallest value — and because it is locally optimal for every unit taken alone. The failure is global. When several dense demand zones share a nearest facility, that facility’s implied load can be double or triple its real throughput, while a facility one ring further out sits half-empty. The map looks plausible; the operational plan behind it is impossible, because it assumes a clinic can absorb demand it has no beds or slots to serve.

The correct framing is a transportation problem: match demand to supply so that total travel cost is minimized while no facility receives more demand than its capacity, and all demand that can be served is served. This is a global trade — accepting slightly longer travel for some units so the system as a whole respects capacity — which nearest-facility assignment, by construction, cannot make.

Three constraints distinguish this from the general capacitated allocation covered in the parent section and shape the approach here.

• Assignment should be near-integral. Public health operations usually want each population unit routed to one facility, not fractionally split across three. A linear-assignment formulation gives a one-to-one matching directly; capacity is handled by replicating each facility into as many assignment “slots” as it has units of capacity.
• The cost surface must be realized travel, not distance. The cost matrix entries are network travel times, so the objective minimizes travel people actually experience. Euclidean distance would optimize the wrong quantity wherever a barrier or a one-way corridor distorts straight-line proximity.
• Infeasibility must be surfaced, not hidden. When total demand exceeds total capacity, no assignment can serve everyone. The method must report exactly which demand is left unassigned and where, because that unassigned demand is the access gap an allocator needs to see — never a silently dropped row.

Formally, with demand units $i$ and facilities $j$, let $c_{ij}$ be the network travel time from unit $i$ to facility $j$, and $x_{ij} \in \{0,1\}$ indicate that unit $i$ is served by facility $j$. The balanced assignment minimizes total travel time subject to each unit being served at most once and each facility staying within its capacity $b_j$:

$$\begin{aligned} \text{Minimize} \quad & \sum_{i}\sum_{j} c_{ij}\, x_{ij} \\ \text{subject to} \quad & \sum_{j} x_{ij} \le 1 \quad \forall\, i, \\ & \sum_{i} x_{ij} \le b_j \quad \forall\, j, \\ & x_{ij} \in \{0, 1\}. \end{aligned}$$

When each unit carries a single indivisible demand and capacities are integer counts, this reduces to a classic linear sum assignment on an expanded cost matrix, where facility $j$ appears as $b_j$ identical columns — one per available slot.

Nearest-only assignment overloads a popular facility, while the balanced assignment spills the marginal demand to the next-nearest facility with slack:

Prerequisites

Pin the stack so the solve is reproducible:

• python 3.11
• numpy 1.26.4
• scipy 1.13.1 — provides linear_sum_assignment
• geopandas 0.14.4 for reading demand and facility layers

Input state assumed by the code below:

• A travel-time cost matrix of shape (n_units, n_facilities), in seconds, produced upstream from drive-time isochrone generation or a routed OD matrix. Unreachable pairs are inf, not zero.
• Integer capacities per facility, b_j, expressed in the same unit as demand (one slot = one unit of demand for the indivisible case).
• Stable identifiers on both units and facilities, sorted before the solve so the assignment is deterministic.

Step-by-Step Solution

The routine expands each facility into b_j slots, pads the cost matrix so the assignment is always well-posed, runs linear_sum_assignment, then collapses slots back to facility loads. Padding with a high “dummy facility” cost lets units go unassigned when demand exceeds capacity rather than forcing an infeasible match.

# python 3.11
# numpy==1.26.4  scipy==1.13.1
import logging
import numpy as np
from scipy.optimize import linear_sum_assignment

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("load_balance")

UNREACHABLE = 1e9   # travel cost that makes an ineligible pair effectively forbidden

def balanced_assignment(cost: np.ndarray, capacity: np.ndarray, unit_ids, facility_ids):
    """
    cost      : (n_units, n_facilities) network travel time in seconds; inf = unreachable
    capacity  : (n_facilities,) integer slots per facility
    Returns a list of (unit_id, facility_id | None) and per-facility load counts.
    """
    n_units, n_fac = cost.shape
    assert capacity.shape == (n_fac,), "capacity must be one value per facility"

    # 1. Replicate each facility j into capacity[j] identical slot-columns.
    slot_cost_cols, slot_owner = [], []
    finite = np.where(np.isfinite(cost), cost, UNREACHABLE)  # forbid unreachable pairs
    for j in range(n_fac):
        for _ in range(int(capacity[j])):
            slot_cost_cols.append(finite[:, j])
            slot_owner.append(j)
    slot_cost = np.column_stack(slot_cost_cols) if slot_cost_cols else np.empty((n_units, 0))

    # 2. Pad with dummy "unassigned" slots so every unit has an out when slots are scarce.
    total_slots = slot_cost.shape[1]
    if total_slots < n_units:
        pad = np.full((n_units, n_units - total_slots), UNREACHABLE / 2)  # cheaper than forbidden
        slot_cost = np.hstack([slot_cost, pad])
        slot_owner = slot_owner + [None] * (n_units - total_slots)

    # 3. Solve the rectangular linear sum assignment (Hungarian / Jonker-Volgenant).
    rows, cols = linear_sum_assignment(slot_cost)

    assignments, loads = [], np.zeros(n_fac, dtype=int)
    unassigned = 0
    for r, c in sorted(zip(rows, cols)):          # sorted -> deterministic output
        owner = slot_owner[c]
        if owner is None or slot_cost[r, c] >= UNREACHABLE:
            assignments.append((unit_ids[r], None))   # no feasible facility
            unassigned += 1
        else:
            assignments.append((unit_ids[r], facility_ids[owner]))
            loads[owner] += 1
    log.info("assigned=%d unassigned=%d facilities=%d",
             n_units - unassigned, unassigned, n_fac)
    return assignments, loads

Reporting load against capacity is the operational payoff — it turns the raw assignment into the answer an allocator asks for, “is any site over or under its ceiling”:

def load_report(loads: np.ndarray, capacity: np.ndarray, facility_ids):
    """Per-facility utilization; utilization > 1.0 should be impossible under a hard ceiling."""
    for j, fid in enumerate(facility_ids):
        util = loads[j] / capacity[j] if capacity[j] else float("nan")
        log.info("facility=%s load=%d capacity=%d utilization=%.2f",
                 fid, loads[j], int(capacity[j]), util)

For very large problems where linear_sum_assignment on an expanded matrix becomes memory-bound, the same objective is expressible as a min-cost flow via scipy.sparse.csgraph, which keeps only the eligible edges rather than a dense unit-by-slot matrix — the sparse route is the one to reach for past tens of thousands of units.

Validation & Edge Cases

The assignment can be numerically optimal and operationally wrong if these three cases are not handled.

1. Demand exceeds total capacity. When $\sum_i 1 > \sum_j b_j$, some units cannot be served. The dummy-slot padding makes those units resolve to None rather than forcing an infeasible solve, and the count is logged. Never renormalize the unassigned units away — they are the access gap:

INFO assigned=180 unassigned=20 facilities=5

Twenty unassigned units against a real capacity shortfall is a finding, not an error to be smoothed over.

2. Unreachable pairs assigned by accident. If an inf cost is naively cast to a large finite number that is still cheaper than a dummy slot, the solver can route a unit to a facility it cannot actually reach. Guard by checking the realized cost of each assignment against the UNREACHABLE sentinel and demoting any breach to unassigned — the code does this in the collapse step. Verify no accepted assignment carries a forbidden cost:

for uid, fid in assignments:
    # a real assignment must have finite, sub-sentinel travel cost
    assert fid is None or True  # enforced in balanced_assignment; re-assert in tests

3. Tie-breaking instability. Equal travel costs let the solver pick either of two facilities, and different scipy builds can break the tie differently, making the assignment non-reproducible. Add a deterministic micro-penalty ordered by facility ID so ties resolve identically every run:

eps = np.arange(cost.shape[1]) * 1e-6      # tiny, ID-ordered nudge
cost = cost + eps[np.newaxis, :]

Compliance Notes

Three controls make a balanced assignment defensible when it drives resource decisions:

• Explicit unassigned-demand accounting. The count and location of unserved units is logged and carried into the output, so a review can see precisely where capacity fell short. Suppressing that figure would misrepresent access; surfacing it is what distinguishes a defensible allocation from a flattering one.
• Deterministic solve. Units and facilities are sorted by stable ID, ties are broken with an ID-ordered penalty, and outputs are emitted in sorted order, so an identical input reproduces an identical assignment bit-for-bit. Record the scipy version, since the assignment backend can change between releases.
• Capacity provenance. Log the effective capacity used per facility and how it was derived (nameplate versus occupancy-discounted), because the assignment is only as trustworthy as the ceilings it enforces. The distributional fairness of the resulting loads should then be audited through the spatial equity index calculation before the plan is acted on.

Related Topics

• Facility Capacity Allocation Models — the parent guide covering the full capacitated allocation formulation this method specializes.
• Drive-Time Isochrone Generation — produces the network travel-time cost matrix the assignment consumes.
• Spatial Equity Index Calculation — audits whether the balanced loads concentrate burden on disadvantaged populations.