Spatial Epi Public Health GIS

Facility Capacity Allocation Models

This guide sits within Healthcare Access & Network Analysis Automation and covers how to convert static facility inventories into demand-responsive, audit-ready service footprints. Facility capacity allocation models are the computational core that matches population demand to bounded facility throughput under real road-network accessibility — replacing one-off mapping with version-controlled pipelines that reconcile real-time patient demand, bed and staffing ceilings, and equity constraints while maintaining a defensible audit trail.

Facility Capacity Allocation Pipeline Five stages connected left to right: Harmonize geospatial data and enforce schema, Build network-constrained service footprints, Apply constraint-based allocation algorithms, Validate statistics and audit equity, and Deploy the production pipeline with monitoring. Harmonize schema enforce, CRS align Service Footprints network- constrained Allocate constraint-based optimization Validate statistics & equity audit Deploy production pipeline Capacity allocation — harmonized supply and demand are matched under network and capacity limits

Concept & Epidemiological Alignment

A facility capacity allocation model answers a single surveillance-relevant question: given the demand surface of a population and the bounded throughput of each facility, which demand belongs to which facility once travel cost and saturation are both binding? The output is not a map of nearest facilities but an assignment — a set of demand-to-facility links whose total load never exceeds any facility’s ceiling and whose aggregate weighted travel time is minimized. This distinguishes it from raw accessibility scoring, which assumes infinite capacity and silently overcounts access wherever a popular facility is already full.

Use a capacity-constrained allocation model when supply is genuinely scarce relative to demand and when the cost of misallocation is clinical — surge periods, mass-vaccination logistics, dialysis or oncology referral planning, or EMS receiving-hospital assignment. Where capacity is abundant and the question is purely “can people reach a facility within t minutes,” a simpler two-step floating catchment or an accessibility index is the right tool and a constrained optimizer is over-engineering. The allocation model earns its complexity precisely when the capacity constraint binds.

Three assumptions must hold for the results to be defensible in a public health context. First, the demand surface must be a credible estimate of need, not just population — age-standardized incidence or service-specific utilization rates, not raw census counts, or the allocation optimizes for the wrong objective. Second, the impedance surface must reflect realized travel, so the model consumes drive-time isochrone generation output rather than Euclidean distance. Third, capacity must be expressed in the same temporal unit as demand (beds-per-day, slots-per-clinic-session); mixing an annual demand figure with an instantaneous bed count produces a meaningless constraint.

Method Selection

The right allocation formulation depends on whether capacity binds, whether assignment must be all-or-nothing, and how large the problem is. Pick the row that matches your operational constraints before writing any solver code.

Method Capacity ceiling Assignment Best when Avoid when
Two-step floating catchment (2SFCA) Soft (ratio) Fractional Screening access deserts; capacity rarely binds A facility’s saturation must hard-block further demand
Gravity / Huff allocation None Probabilistic Modeling realistic choice spillover across many sites You need a guaranteed, auditable single assignment per zone
Capacitated p-median (LP/MILP) Hard Binary or fractional Scarce supply, surge planning, referral routing N is very large and a fractional answer is acceptable
Min-cost flow / transportation Hard Fractional Tens of thousands of zones; divisible demand Each zone must route to exactly one facility

For the production walkthrough below the working assumption is a binding hard ceiling with a near-single assignment per demand zone, which lands on the capacitated p-median formulation solved as a mixed-integer program; the hospital bed capacity model is the surge-specialized instance of exactly this formulation.

Spatial Data Prerequisites

Allocation is only as trustworthy as the geometry and attributes feeding it. Three input layers must be prepared and gated before optimization:

  • Facility points — a point geometry per facility carrying a stable facility_id, a capacity value in the same temporal unit as demand, and a service_type. Snap each point to the routable road graph; an unsnapped facility silently drops out of every catchment.
  • Demand zones — polygon geometries (census tracts, hex bins, or service areas) with a demand attribute derived from need, not raw population, plus the socioeconomic covariates needed for the later equity audit.
  • Impedance surface — the network travel-time matrix or isochrone polygons that quantify reachable cost between each demand zone centroid and each facility.

Before any distance calculation, enforce a canonical projected coordinate reference system for public health — an equal-area projection such as EPSG:5070 (CONUS Albers) so areas and network distances are metric and undistorted. Run topology cleaning (buffer(0) to repair self-intersections, drop zero-area slivers and invalid multipart features) and verify a minimum viable problem size: at least one reachable facility per demand zone within the impedance ceiling, or the LP is infeasible. Validate input attributes against an explicit schema so a missing capacity or a negative demand is caught at ingestion rather than surfacing as a silent solver failure. These data-type and format expectations follow the site’s spatial data types and formats conventions.

The ingestion gate harmonizes CRS, repairs geometry, and validates the schema deterministically, logging provenance for HIPAA and GDPR audit:

# pinned: geopandas==0.14.4, shapely==2.0.4, pyproj==3.6.1, pydantic==2.7.1
import geopandas as gpd
from pydantic import BaseModel, ValidationError, field_validator
import logging, json

log = logging.getLogger("alloc.ingest")

class FacilitySchema(BaseModel):
    facility_id: str
    capacity: int            # same temporal unit as demand (e.g. beds/day)
    service_type: str
    operational_status: str

    @field_validator("capacity")
    @classmethod
    def _positive(cls, v: int) -> int:
        if v <= 0:
            raise ValueError("capacity must be positive")
        return v

def harmonize_and_validate(gdf: gpd.GeoDataFrame, target_crs: str = "EPSG:5070") -> gpd.GeoDataFrame:
    gdf = gdf.to_crs(target_crs)                      # equal-area, metric distances
    gdf["geometry"] = gdf.geometry.buffer(0)         # repair self-intersections
    gdf = gdf[gdf.geometry.is_valid & ~gdf.geometry.is_empty].copy()

    valid = []
    for _, row in gdf.sort_values("facility_id").iterrows():   # deterministic order
        try:
            FacilitySchema(**row.drop("geometry").to_dict())
            valid.append(row)
        except ValidationError as e:
            log.warning(json.dumps({"facility_id": row.get("facility_id"),
                                     "event": "schema_reject", "errors": e.errors()}))
    out = gpd.GeoDataFrame(valid, crs=target_crs)
    log.info(json.dumps({"event": "ingest_summary",
                         "input_rows": len(gdf), "valid_rows": len(out)}))
    return out

Network-Constrained Service Footprints

Capacity allocation cannot rely on Euclidean distance: real accessibility depends on road-network topology, time-of-day speed profiles, and modality-specific routing. The allocation pipeline consumes the isochrone output described under drive-time isochrone generation to establish service boundaries from actual travel time rather than straight-line proximity. In production this means interfacing with OSRM, Valhalla, or ArcGIS Network Analyst via REST or local Docker containers, caching the network graph, applying time-dependent speed profiles, and handling seasonal closures and detours as explicit edge cases.

Containerized Valhalla instances give deterministic, auditable results without external rate limits. Generate isochrone polygons at standardized intervals (10, 20, 30 minutes), intersect them with the population demand surface to compute reachable catchment volumes, and persist the resulting zone-to-facility travel-time matrix as the cost input for the optimizer. That matrix — not the geometry — is what the allocation solver actually reads.

Production Implementation: Constraint-Based Allocation

Once network accessibility is quantified, allocation becomes a mathematical optimization. The standard formulation is a capacitated p-median assignment that distributes demand across facilities while respecting throughput ceilings: every demand zone is matched to a facility under that facility’s capacity limit, so the optimizer trades nearest-facility preference against saturation — the moment a facility fills, surrounding demand spills to the next-reachable site.

Capacitated Allocation: Demand-to-Facility Assignment Under Throughput Limits Four demand zones on the left, each with a demand volume, are assigned by minimum network travel time to three facilities on the right, each with a bed capacity. Zone A and Zone B route to Facility 1; Zone C routes to Facility 2; Zone D would prefer Facility 2 but it is saturated, so its demand spills to Facility 3 — illustrating how the capacity constraint overrides nearest-facility preference. Demand zones (dᵢ) Facilities (capacity cⱼ) Zone A demand 80 Zone B demand 60 Zone C demand 90 Zone D demand 70 Facility 1 cap 160 · load 140 Facility 2 cap 90 · load 90 (full) Facility 3 cap 120 · load 70 nearest, but cⱼ exhausted spill assignment Capacity constraint Σ dᵢ·xᵢⱼ ≤ cⱼ overrides minimum-travel-time preference, redirecting saturated demand

The objective minimizes total weighted travel time while enforcing each facility’s capacity ceiling and assigning every demand zone exactly once:

Minimize:∑i∑jdemandi⋅distanceij⋅xijSubject to:∑jxij=1,∀ i∑idemandi⋅xij≤capacityj,∀ jxij∈{0,1}\begin{aligned} \text{Minimize:} \quad & \sum_{i}\sum_{j} \text{demand}_i \cdot \text{distance}_{ij} \cdot x_{ij} \\ \text{Subject to:} \quad & \sum_{j} x_{ij} = 1, \quad \forall\, i \\ & \sum_{i} \text{demand}_i \cdot x_{ij} \le \text{capacity}_j, \quad \forall\, j \\ & x_{ij} \in \{0, 1\} \end{aligned}

The implementation below builds the assignment in PuLP over the network travel-time matrix, sorts every index by a stable ID so the model is deterministic, and logs the solver’s convergence status and objective value for the audit trail:

# pinned: pulp==2.8.0, pandas==2.2.2, numpy==1.26.4
import pulp, pandas as pd, hashlib, json, logging

log = logging.getLogger("alloc.solve")

def allocate(demand: pd.Series, capacity: pd.Series, tt: pd.DataFrame,
             penalty: float | None = None) -> pd.DataFrame:
    """
    demand  : zone_id -> demand units
    capacity: facility_id -> ceiling (same temporal unit as demand)
    tt      : DataFrame [zone_id x facility_id] of network travel time (minutes)
    penalty : per-unit cost for soft over-capacity (None = hard constraint)
    """
    zones = sorted(demand.index)          # deterministic ordering
    facs  = sorted(capacity.index)

    prob = pulp.LpProblem("capacitated_allocation", pulp.LpMinimize)
    x = {(i, j): pulp.LpVariable(f"x_{i}_{j}", cat="Binary") for i in zones for j in facs}

    over = {}
    if penalty is not None:               # soft-constraint surge slack
        over = {j: pulp.LpVariable(f"over_{j}", lowBound=0) for j in facs}

    prob += pulp.lpSum(demand[i] * tt.at[i, j] * x[(i, j)] for i in zones for j in facs) \
            + pulp.lpSum(penalty * over[j] for j in over)

    for i in zones:                       # each zone assigned exactly once
        prob += pulp.lpSum(x[(i, j)] for j in facs) == 1
    for j in facs:                        # capacity ceiling (hard or soft)
        load = pulp.lpSum(demand[i] * x[(i, j)] for i in zones)
        prob += load <= capacity[j] + (over[j] if j in over else 0)

    status = prob.solve(pulp.PULP_CBC_CMD(msg=False))
    rows = [{"zone_id": i, "facility_id": j,
             "demand": demand[i], "travel_min": tt.at[i, j]}
            for i in zones for j in facs if x[(i, j)].value() == 1]
    result = pd.DataFrame(rows).sort_values(["zone_id", "facility_id"])

    cfg = json.dumps({"zones": zones, "facilities": facs, "penalty": penalty}, sort_keys=True)
    log.info(json.dumps({"event": "solve_complete",
                         "status": pulp.LpStatus[status],
                         "objective": pulp.value(prob.objective),
                         "config_sha256": hashlib.sha256(cfg.encode()).hexdigest()}))
    return result

Parameter Selection & Tuning

Three parameter families govern allocation behavior and each must be set deliberately, not left at defaults. The capacity ceiling per facility should be expressed at the same temporal granularity as demand and discounted for realistic occupancy — a hospital’s licensed bed count is not its usable surge capacity once diversion thresholds and staffing ratios are applied. Encode the effective ceiling, not the nameplate figure, and log which discount was used.

The impedance ceiling (the maximum travel time at which a zone-facility pair is even eligible) is the lever that controls feasibility versus equity. Set it too tight and rural zones become infeasible; too loose and the optimizer routes patients past clinically reasonable limits. Service-specific thresholds matter: ICU and trauma assignments warrant stricter ceilings than general medical-surgical, and pediatric or stroke-ready routing carries its own weights.

The soft-constraint penalty governs surge behavior. With a hard ceiling the LP is infeasible the instant aggregate demand exceeds aggregate capacity; introduce the over_j slack variable with a per-unit penalty to let the model overflow gracefully and report where and by how much capacity is breached. Tune the penalty and the impedance ceiling together with an automated grid or Bayesian sweep, scoring each combination on total weighted travel time, the count of over-capacity facilities, and the equity spread, so the chosen configuration trades efficiency against resilience explicitly rather than by accident.

Statistical Validation & Equity Auditing

Allocation outputs require statistical validation before operational deployment. Compute spatial autocorrelation (Global and Local Moran’s I) on assignment residuals to detect clustering of underserved zones, and run residual analysis between predicted and observed utilization to confirm calibration. For distributional fairness across socioeconomic strata, integrate the spatial equity index calculation — Gini coefficients and Lorenz curves quantify whether the allocation concentrates burden on already-disadvantaged populations.

Validation pipelines should execute automated hypothesis testing:

  • Kolmogorov–Smirnov tests for travel-time distribution alignment against the prior period
  • Chi-square goodness-of-fit for facility utilization patterns versus expected load
  • Monte Carlo simulation to stress-test the assignment under stochastic demand shocks

Version-control every validation metric alongside the model parameters that produced it, so a regulatory reviewer can reproduce both the allocation and the evidence that it passed.

Edge Cases & Failure Modes

Production allocation pipelines fail in characteristic ways, and each needs an explicit guard:

  • Infeasible LP (demand exceeds capacity). A hard-constraint model returns Infeasible the moment Σ demand > Σ reachable capacity. Detect the status, fall back to the soft-penalty formulation, and surface the over-capacity facilities rather than emitting an empty assignment.
  • Island zones with no reachable facility. A demand zone whose nearest facility lies beyond the impedance ceiling makes the Σ_j x_ij = 1 constraint impossible. Pre-screen for orphan zones, log them as access gaps, and relax the ceiling for those zones only — never silently drop them.
  • Transboundary CRS drift. Facilities and demand zones ingested in different projections produce nonsense travel distances. Enforce one canonical projected CRS at ingestion and assert it on every layer before the solver runs.
  • Zero or null demand zones. Zones with zero demand inflate the variable count without affecting the objective; filter them before building the model. Null demand must error at the schema gate, never default to zero.
  • Memory and runtime for large N. A binary p-median has |zones| × |facilities| variables; past roughly 50,000 zones the MILP becomes intractable. Switch to a min-cost-flow / transportation formulation with fractional assignment, or pre-cluster demand into representative bins, and cap the solver wall-clock with a logged time limit so a runaway solve cannot stall the pipeline.

Compliance & Audit Controls

Operational deployment requires a fault-tolerant, audit-ready architecture. Orchestrate batch runs via Apache Airflow or Prefect with explicit retry logic and exponential backoff for routing-API failures, and run every spatial operation in isolated Docker containers with pinned dependency versions to prevent environment drift.

Determinism is a compliance requirement, not a convenience: sort zones and facilities by stable IDs before building the model (as the code above does) so re-running identical inputs reproduces an identical assignment bit-for-bit. Hash the full configuration — zone set, facility set, capacity ceilings, impedance ceiling, penalty — with SHA-256 and record the digest with every run, so any allocation decision can be tied back to the exact parameters that produced it.

Emit structured JSON logs at every stage capturing input dataset checksums and row counts, CRS transformation and validation flags, routing-engine version and graph build timestamp, solver convergence status and objective value, and equity-metric outputs with any threshold violations. Persist the assignment output with an explicit, documented schema — zone_id, facility_id, demand, travel_min, run_sha256, model_version — and attach ISO 19115 lineage metadata (source, transformation steps, parameter configuration) so the data product is self-describing. Ship logs to a centralized SIEM or object store with immutable retention to satisfy both internal governance and external HIPAA / GDPR audit.