Spatial Epi Public Health GIS

Siting New Clinics With P-Median Optimization

This guide, part of Facility Capacity Allocation Models, covers location–allocation for a common public health decision: where to put a small number of new clinics.

Problem Context & Constraints

A county has funding for three new community health clinics and a shortlist of 40 candidate sites — county-owned buildings, leased storefronts, school-based options. Planners want the sites that bring the most people closest to care. Picking by intuition tends to favour visible neighbourhoods and politically active ones; picking by optimisation makes the objective explicit and the trade-offs visible.

The p-median problem chooses p sites from the candidates to minimise the total demand-weighted travel time from every demand point to its nearest open site. With demand weighted by the population in need — uninsured residents, Medicaid enrollees, or a composite need index — the solution puts clinics where they shorten the most need-weighted trips. It is a classic integer program: a binary variable for each candidate (open or not) and a binary variable for each demand–site assignment.

Three constraints shape a usable model. Existing facilities must be included as fixed open sites, or the model will re-site them. Coverage floors — no one in need more than a set time from a clinic — express equity goals that minimising the average alone will not guarantee. And capacity, where a single site cannot absorb unlimited demand, turns the problem into a capacitated p-median, which assigns demand to sites under their capacity rather than simply to the nearest.

Choosing Three Sites From Forty A county map with demand areas shaded by need, two existing clinics shown as filled squares and forty candidate sites shown as small open circles. The p-median solution opens three candidates, shown as large filled circles: one in a high-need area far from existing clinics, one on the edge of a dense neighbourhood, and one serving a rural corridor. Each demand area is assigned to its nearest open clinic, reducing total need-weighted travel time by thirty-eight percent. Existing clinics fixed, three new sites chosen high need, far dense edge rural corridor squares: existing clinics; large dots: new sites; shading: need-weighted demand

Prerequisites

  • python 3.11
  • pulp 2.8.0 with the bundled CBC solver, numpy 1.26.4, pandas 2.2.2
  • Demand points — block groups or tracts — with a need weight such as uninsured population
  • Candidate sites vetted for feasibility: zoning, lease availability, transit access
  • Existing clinic locations
  • A travel time matrix from every demand point to every candidate and existing site, built as in Caching OSRM Route Matrices for Large Batches

Step-by-Step Solution

Formulate the p-median with existing sites fixed open, add an optional maximum travel time constraint, solve with CBC, and compare the result against the status quo and against simple heuristics.

# P-median clinic siting with fixed existing sites and a coverage floor.
# Pinned: pulp==2.8.0, numpy==1.26.4, pandas==2.2.2
import logging

import numpy as np
import pandas as pd
import pulp

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


def solve_pmedian(T: pd.DataFrame, weight: pd.Series, existing: list[str], p_new: int,
                  max_minutes: float | None = None, time_limit: int = 600) -> dict:
    """T: minutes, rows = demand ids, columns = site ids (candidates + existing)."""
    I, J = list(T.index), list(T.columns)
    candidates = [j for j in J if j not in existing]

    # Prune assignments that can never be optimal or are forbidden by the floor.
    allowed = {(i, j) for i in I for j in J
               if max_minutes is None or T.at[i, j] <= max_minutes}
    uncoverable = [i for i in I if not any((i, j) in allowed for j in J)]
    if uncoverable:
        raise ValueError(f"{len(uncoverable)} demand points have no site within {max_minutes} min")

    m = pulp.LpProblem("p_median", pulp.LpMinimize)
    y = pulp.LpVariable.dicts("open", J, cat="Binary")
    x = pulp.LpVariable.dicts("assign", list(allowed), lowBound=0, upBound=1)

    m += pulp.lpSum(weight[i] * T.at[i, j] * x[(i, j)] for (i, j) in allowed)
    for i in I:
        m += pulp.lpSum(x[(i, j)] for j in J if (i, j) in allowed) == 1
    for (i, j) in allowed:
        m += x[(i, j)] <= y[j]
    for j in existing:
        m += y[j] == 1
    m += pulp.lpSum(y[j] for j in candidates) == p_new

    status = m.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=time_limit))
    opened = [j for j in candidates if y[j].value() > 0.5]
    assign = {i: max((j for j in J if (i, j) in allowed), key=lambda j: x[(i, j)].value() or 0)
              for i in I}
    obj = pulp.value(m.objective)
    log.info("status=%s new_sites=%s objective=%.0f", pulp.LpStatus[status], opened, obj)
    return {"opened": opened, "assign": assign, "objective": obj,
            "status": pulp.LpStatus[status]}


def evaluate(T, weight, open_sites) -> dict:
    nearest = T[open_sites].min(axis=1).to_numpy()
    w = (weight.reindex(T.index) / weight.sum()).to_numpy()
    order = np.argsort(nearest)
    p90 = nearest[order][np.searchsorted(np.cumsum(w[order]), 0.9)]
    return {"mean_min": float((nearest * w).sum()),
            "p90_min": float(p90),
            "share_over_30": float(w[nearest > 30].sum())}

Assignment variables are continuous here, not binary. For the uncapacitated p-median, the optimal solution assigns each demand point wholly to its nearest open site anyway, so relaxing the assignment variables is exact and makes the problem much faster to solve. For a capacitated version, where demand can be split across sites, the assignments must be made binary if split assignment is not acceptable.

The evaluate function reports more than the objective. The mean travel time is what the p-median minimises; the 90th percentile and the share of need beyond 30 minutes show what happens to the worst-served, which is where public health equity goals usually sit.

Efficiency and Equity Pull Apart Three solutions compared. The status quo has a need-weighted mean travel time of twenty-two minutes and a ninetieth percentile of forty-eight minutes, with nineteen percent of need beyond thirty minutes. The unconstrained p-median cuts the mean to fourteen minutes and the ninetieth percentile to thirty-six, with nine percent beyond thirty. Adding a thirty-minute floor raises the mean slightly to fifteen minutes but brings the ninetieth percentile to twenty-eight and leaves no one beyond thirty minutes. A small efficiency cost buys the equity floor Mean 90th pct Beyond 30 min Status quo 22 min 48 min 19% P-median 14 min 36 min 9% P-median + 30-min floor 15 min 28 min 0% need-weighted travel times; present all three to decision-makers

Validation & Edge Cases

Solver status. Always check the status. “Optimal” means the solution is proven best; “Not Solved” with a time limit means CBC stopped with the best solution found so far and an optimality gap. Report the gap. For a few hundred demand points and 40 candidates, CBC normally proves optimality in seconds; thousands of demand points with hundreds of candidates may need a longer limit or aggregation.

Infeasible coverage floors. If some demand points have no candidate or existing site within the floor, the model is infeasible. The code detects this before solving and lists them; the planning answer is either to add candidates in that area or to relax the floor there explicitly.

Sensitivity to p. Solve for p from one to five and plot the objective. The marginal gain of each additional clinic usually falls steeply; showing the curve helps decision-makers see whether a fourth site is worth seeking funding for.

Near-optimal alternatives. Several site combinations are often within one or two percent of the optimum. Those alternatives matter because a planner may reject the optimal site for reasons outside the model — a lease falls through, a neighbourhood objects. Enumerate near-optimal solutions by adding a constraint that excludes the optimum and re-solving, and present the top few.

Demand weights. Weighting by total population and by uninsured population produce different answers. The choice encodes the policy goal, so it should be agreed with program leadership before the model runs, not chosen afterwards to match a preferred site.

Travel mode. Clinics for low-income populations are often reached by transit or on foot. A drive-time matrix will favour sites near highways. Build the matrix in the mode the target population uses, or in a mix weighted by vehicle availability, using Building GTFS Travel-Time Matrices for Clinic Access.

Candidate lists that pre-decide the answer. The optimiser can only choose among the sites it is given. If candidates cluster in one part of the county because that is where the facilities department owns buildings, the solution will too. Review the candidate map for gaps before solving, and where a high-need area has no candidate, add a notional site there to learn how much the missing option costs — even if it cannot be built, the difference is useful evidence for acquiring one.

Capacity. Without a capacity constraint, a site in a dense area may be assigned more patients than it can see. Add capacities per site and move to a capacitated formulation, as in Balancing Facility Loads With Linear Assignment.

INFO pmedian status=Optimal new_sites=['C07', 'C22', 'C31'] objective=412880
INFO pmedian floor=30 min: status=Optimal new_sites=['C07', 'C19', 'C31'] objective=438215
What Each Additional Clinic Buys A curve of total need-weighted travel time against the number of new clinics, from zero to five. The first clinic reduces total travel time by about twenty percent, the second by about twelve percent more, the third by seven percent, and the fourth and fifth by three and two percent. The curve flattens after three, which is the budgeted number, suggesting a fourth clinic would bring much smaller gains. Diminishing returns from each new site budgeted: 3 0 1 2 3 4 5 total weighted travel time number of new clinics

Compliance Notes

Record the candidate list and how it was vetted, the demand geography and need weight, the travel time matrix source, mode and date, the constraints, solver version, status and gap, the chosen sites, near-optimal alternatives, and the evaluation metrics for each.

Siting decisions are public decisions. Publish the objective, the weights and the constraints in plain language with the result, so residents can see that the choice followed stated goals. Where the final decision departs from the model’s recommendation, record why; the model supports the decision, it does not make it.