How to Choose Facility Locations with a p-Median Model in Python

On this page

Problem statement

You need to open a fixed number of facilities — clinics, depots, libraries — and want the choice that makes the average trip as short as possible for everyone. The p-median model finds that choice exactly. It is also the location model that grows fastest: it has a variable for every pairing of a demand point with a candidate site, and at realistic sizes the solver runs for minutes or never returns.

Measured by choosing 8 of 106 school sites in Chittenden County, Vermont, for the 168,323 residents of 2,241 census blocks, with drive times on the road network:

The optimal 8 sites gave a mean trip of 5.26 minutes, against 2.23 minutes if every candidate were open.
The model had 237,440 variables and took 88 seconds to solve.
Doubling the candidates from 100 to 200 multiplied the solve time by more than six. 200 road junctions took 257 seconds; a different 200 candidates in the troubleshooting guide produced no solution at all in 5 minutes.
Solving on 118 block groups took under a second, and its sites were only 0.81% worse when measured on the blocks.

Quick answer

Build a demand × candidate travel-time matrix, solve with a time limit, then evaluate the chosen sites on the full demand:

import pulp
from spopt.locate import PMedian

model = PMedian.from_cost_matrix(C, w, p_facilities=8)
model.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=900))
sites = [j for j, v in enumerate(model.fac_vars) if v.value() > 0.5]
nearest = C[:, sites].min(axis=1)
mean_minutes = (nearest * w).sum() / w.sum() / 60

C holds travel times in seconds from each demand point to each candidate, and w the population of each demand point. If C has more than a few hundred thousand cells, aggregate the demand first.

The size check comes before the solve: the variable count predicts whether the solver will return.

What you get — from The data used on this page

dropping 2 blocks (553 people) with no road route

Weighted demand, candidates, a matrix, a size check

1. Represent demand as weighted points

Use the smallest areas you have population for, each as one point with its population as the weight. Here that is 2,241 populated census blocks, each represented by a point inside the block. Population-weighted points are what the objective averages over, so an empty block adds nothing and a dense one counts in proportion.

2. Choose the candidates deliberately

The model can only pick from the list it is given. The 106 school sites are a plausible list of publicly owned land. As a comparison, 100 road junctions drawn at random from the county gave a slightly higher optimal mean, 5.28 minutes against 5.26, and 200 of them a lower one, 5.10. Which sites are on the list matters as much as how many.

3. Build the travel-time matrix

Run one shortest-path search per candidate on the reversed road graph: every demand point's time to that candidate falls out of a single search. For 2,239 blocks and 106 schools the matrix took about 60 milliseconds. Check it for infinite values before solving — a demand point that cannot reach any candidate makes the model meaningless (Example 1).

4. Check the size before solving

The p-median model has one assignment variable for every demand–candidate pair plus one variable per candidate, so its size is n × k + k. On the same 2,239 blocks, with an 8-site budget and random road junctions as candidates (Example 4):

candidates   variables   solve time   status          mean trip
     25         56,000        5.7 s   Optimal          6.07 min
     50        112,000       11.8 s   Optimal          5.90 min
    100        224,000       39.0 s   Optimal          5.28 min
    106 schools 237,440      88.2 s   Optimal          5.26 min
    200        448,000      256.9 s   Optimal          5.10 min

Solve time doubled from 25 to 50 junctions, more than tripled from 50 to 100, and grew six-and-a-half-fold from 100 to 200. That last solve finished inside the 300-second limit, but only just: one run while other jobs shared the machine was stopped at the limit with no solution; the troubleshooting guide's 200 candidates — the 106 schools plus 94 blocks — did not, and spopt raised RuntimeError: Model is not solved: Not Solved. The variable count predicts the trouble; it does not decide it.

5. Solve with a time limit and read the status

Always pass timeLimit, and read pulp.LpStatus[model.problem.status] afterwards. When the limit stops CBC after it has found a solution, PuLP still reports Optimal; model.problem.sol_status of 2 rather than 1 is the sign that optimality was not proved.

6. Evaluate the sites on more than the objective

The optimal sites were [24, 49, 51, 63, 68, 70, 72, 77] among the 106 candidates. Measured on the blocks, the mean trip was 5.26 minutes, the 95th percentile 13.72 minutes (computed in location-allocation explained), the longest trip 26.99 minutes, and 89.0% of residents were within 10 minutes. The model minimises the mean and nothing else; the longest trip is whatever that choice leaves.

7. Aggregate demand when the model is too big

Replace blocks with block groups or tracts, represented by population-weighted centroids, solve on those, and evaluate the chosen sites back on the blocks:

demand units   variables   solve    own-level mean   mean on blocks   sites shared with block solution
2,239 blocks     237,440    88.2 s      5.261 min       5.261 min            8 of 8
118 groups        12,614     0.58 s     5.107 min       5.302 min (+0.81%)   6 of 8
40 tracts          4,346     0.13 s     4.391 min       5.423 min (+3.10%)   2 of 8

The block-group solution was over a hundred times faster and only 0.81% worse. The tract solution looked best on its own terms — 4.39 minutes — and was worst in reality, because a tract's centroid hides how far its residents actually travel.

8. Compare with other objectives before deciding

p-median sites sit near population centres and tolerate long trips at the edges. If the longest trip or a service standard matters, compare with p-centre and maximal coverage on the same matrix, as in location-allocation explained.

Block groups kept almost all of the quality at a fraction of the cost; tracts did not.

Solving and describing a p-median plan

The data used on this page

Chittenden County, Vermont. Demand is the 2020 Census: the TIGER/Line 2020 blocks for Vermont, each populated block as a representative point with its POP20. The candidates are every OpenStreetMap amenity=school feature in the county — 106 of them — and, for the size test in Example 4, road intersections drawn with a fixed seed. Drive times come from OpenStreetMap's drive network for the county plus 10.5 km, with OSMnx's imputed free-flow speeds, road-only between nearest nodes. Two blocks (553 people) whose node reaches no school are dropped before any model sees them. Solved with spopt 0.7 and PuLP 3.3's bundled CBC (spopt fails on PuLP 4.0: install spopt "pulp<4"). Overpass via the maps.mail.ru mirror.

import io
import os
import requests
import urllib.request
import zipfile

import geopandas as gpd
import numpy as np
import osmnx as ox
import scipy.sparse as sp
from scipy.sparse.csgraph import dijkstra

# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-p-median-facility-location-python-2026-10-09.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.overpass_rate_limit = False
ox.settings.http_user_agent = "spatialworkflow.io guide examples (https://www.spatialworkflow.io/)"
CRS = 32145                                               # NAD83 / Vermont, metres

if not os.path.exists("tl_2020_50_tabblock20.zip"):
    urllib.request.urlretrieve("https://www2.census.gov/geo/tiger/TIGER2020/TABBLOCK20/"
                               "tl_2020_50_tabblock20.zip", "tl_2020_50_tabblock20.zip")
vt_blocks = gpd.read_file("tl_2020_50_tabblock20.zip").to_crs(CRS)
county = vt_blocks[vt_blocks.COUNTYFP20 == "007"].union_all()           # Chittenden
blocks = vt_blocks[(vt_blocks.COUNTYFP20 == "007") & (vt_blocks.POP20 > 0)].copy()
blocks["geometry"] = blocks.geometry.representative_point()

county_ll = gpd.GeoSeries([county], crs=CRS).to_crs(4326).iloc[0]
schools = ox.features_from_polygon(county_ll, tags={"amenity": "school"}).to_crs(CRS)
schools["geometry"] = schools.geometry.representative_point()
schools = schools.sort_index().reset_index()[["element", "id", "name", "geometry"]]

reach = gpd.GeoSeries([county.buffer(10_500)], crs=CRS).to_crs(4326).iloc[0]
G = ox.add_edge_travel_times(ox.add_edge_speeds(ox.graph_from_polygon(reach, network_type="drive")))
G = ox.project_graph(G, to_crs=CRS)
nodes = list(G.nodes)
index = {n: i for i, n in enumerate(nodes)}
node_xy = np.array([(G.nodes[n]["x"], G.nodes[n]["y"]) for n in nodes])
best = {}
for u, v, d in G.edges(data=True):                      # shortest of any parallel edges
    best[(index[u], index[v])] = min(best.get((index[u], index[v]), np.inf), d["travel_time"])
rows, cols = zip(*best)
graph = sp.csr_array((list(best.values()), (rows, cols)), shape=(len(nodes), len(nodes)))
school_nodes = np.array([index[n] for n in ox.distance.nearest_nodes(G, schools.geometry.x, schools.geometry.y)])
all_block_nodes = np.array([index[n] for n in ox.distance.nearest_nodes(G, blocks.geometry.x, blocks.geometry.y)])

reachable = np.isfinite(dijkstra(graph.T.tocsr(), directed=True, indices=np.unique(school_nodes),
                                 min_only=True)[all_block_nodes])
print(f"dropping {(~reachable).sum()} blocks ({blocks.POP20[~reachable].sum():,} people) with no road route")
blocks, block_nodes = blocks[reachable].reset_index(drop=True), all_block_nodes[reachable]
w = blocks["POP20"].to_numpy(float)

dropping 2 blocks (553 people) with no road route

Example 1 — a travel-time matrix from a road graph

import numpy as np
from scipy.sparse.csgraph import dijkstra


def travel_time_matrix(graph, demand_nodes, candidate_nodes):
    """Seconds from each demand node to each candidate node on a directed sparse graph."""
    unique, inverse = np.unique(candidate_nodes, return_inverse=True)
    to_candidates = dijkstra(graph.T.tocsr(), directed=True, indices=unique)   # candidate x node
    C = np.ascontiguousarray(to_candidates[inverse][:, demand_nodes].T)
    stranded = int(np.isinf(C).all(axis=1).sum())
    if stranded:
        raise ValueError(f"{stranded} demand points cannot reach any candidate")
    return C


import time

start = time.perf_counter()
C = travel_time_matrix(graph, block_nodes, school_nodes)
print(f"{C.shape[0]:,} x {C.shape[1]} matrix in {(time.perf_counter() - start) * 1000:.0f} ms")

2,239 x 106 matrix in 61 ms

graph is a SciPy sparse matrix of travel times between road nodes, and the node arguments are row indices into it. Reversing the graph with .T turns "from each candidate" into "to each candidate", which matters on one-way streets. Called with the 2,239 routable block nodes and 106 school nodes, it returned a 2,239 × 106 matrix in about 60 ms.

Example 2 — solve, then describe the answer

import time

import pulp
from spopt.locate import PMedian


def solve_p_median(C, w, p, time_limit=900):
    n, k = C.shape
    print(f"{n:,} demand points x {k} candidates = {n * k + k:,} variables")
    started = time.perf_counter()
    model = PMedian.from_cost_matrix(C, w, p_facilities=p)
    model.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=time_limit))
    sites = [j for j, v in enumerate(model.fac_vars) if v.value() > 0.5]
    print(f"{pulp.LpStatus[model.problem.status]} in {time.perf_counter() - started:.1f} s: sites {sites}")
    return sites


def describe(C, w, sites):
    minutes = C[:, sites].min(axis=1) / 60
    load = np.bincount(C[:, sites].argmin(axis=1), weights=w, minlength=len(sites))
    print(f"mean {(minutes * w).sum() / w.sum():.2f} min, longest {minutes.max():.2f} min, "
          f"within 10 min {w[minutes <= 10].sum() / w.sum():.1%}")
    print("people per site:", ", ".join(f"{x:,.0f}" for x in load))


sites = solve_p_median(C, w, p=8)
describe(C, w, sites)

2,239 demand points x 106 candidates = 237,440 variables
Optimal in 88.2 s: sites [24, 49, 51, 63, 68, 70, 72, 77]
mean 5.26 min, longest 26.99 min, within 10 min 89.0%
people per site: 14,086, 47,785, 28,626, 13,405, 14,081, 9,695, 23,604, 16,488

The solve time includes building the model and varies from run to run. The loads matter because the model has no capacities: the busiest site served 47,785 people, nearly five times the 9,695 of the quietest.

Example 3 — solve on block groups, evaluate on blocks

import geopandas as gpd
from scipy.spatial import cKDTree


def aggregate_demand(blocks, key_length=12):
    """Population-weighted centroids of block groups (the first 12 characters of a block GEOID)."""
    b = blocks.assign(unit=blocks.GEOID20.str[:key_length],
                      wx=blocks.geometry.x * blocks.POP20, wy=blocks.geometry.y * blocks.POP20)
    g = b.groupby("unit").agg(POP20=("POP20", "sum"), wx=("wx", "sum"), wy=("wy", "sum")).reset_index()
    return gpd.GeoDataFrame(g[["unit", "POP20"]], crs=blocks.crs,
                            geometry=gpd.points_from_xy(g.wx / g.POP20, g.wy / g.POP20))


groups = aggregate_demand(blocks)
_, group_nodes = cKDTree(node_xy).query(np.column_stack([groups.geometry.x, groups.geometry.y]))
C_groups = travel_time_matrix(graph, group_nodes, school_nodes)
group_sites = solve_p_median(C_groups, groups.POP20.to_numpy(float), p=8)
describe(C, w, group_sites)

118 demand points x 106 candidates = 12,614 variables
Optimal in 0.4 s: sites [24, 32, 51, 63, 68, 70, 77, 87]
mean 5.30 min, longest 26.99 min, within 10 min 89.6%
people per site: 14,068, 30,755, 26,785, 13,952, 15,475, 9,695, 15,819, 41,221

blocks holds one point per populated block, node_xy the road node coordinates in the same projected CRS, and describe is called with the block matrix C, so the comparison is on the same footing as Example 2.

Example 4 — the size test and the aggregation table

def mean_minutes(C, w, sites):
    return (C[:, sites].min(axis=1) * w).sum() / w.sum() / 60


# candidates: road intersections (three or more streets) drawn with a fixed seed
junctions = np.array([index[n] for n, d in G.nodes(data=True) if d.get("street_count", 0) >= 3])
rng = np.random.default_rng(0)
picked = rng.choice(junctions, 200, replace=False)
C_int = dijkstra(graph.T.tocsr(), directed=True, indices=picked)[:, block_nodes].T
C_int[np.isinf(C_int)] = C_int[np.isfinite(C_int)].max() * 10       # a junction some blocks cannot reach

for k in (25, 50, 100, 200):
    started = time.perf_counter()
    model = PMedian.from_cost_matrix(C_int[:, :k], w, p_facilities=8)
    try:
        model.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=300))
        status, chosen = pulp.LpStatus[model.problem.status], [j for j, v in enumerate(model.fac_vars) if v.value() > 0.5]
        mean = f"{mean_minutes(C_int[:, :k], w, chosen):.2f} min"
    except RuntimeError:
        status, mean = "Not Solved", "-"
    print(f"{k:4d} junctions {C.shape[0] * k + k:>9,} variables {time.perf_counter() - started:7.1f} s  {status:11} {mean}")

# the same model on blocks, block groups and tracts, evaluated on the blocks
for label, length in (("block groups", 12), ("tracts", 11)):
    units = aggregate_demand(blocks, key_length=length)
    _, unit_nodes = cKDTree(node_xy).query(np.column_stack([units.geometry.x, units.geometry.y]))
    C_units = travel_time_matrix(graph, unit_nodes, school_nodes)
    started = time.perf_counter()
    model = PMedian.from_cost_matrix(C_units, units.POP20.to_numpy(float), p_facilities=8)
    model.solve(pulp.PULP_CBC_CMD(msg=False))
    s_units = [j for j, v in enumerate(model.fac_vars) if v.value() > 0.5]
    own = mean_minutes(C_units, units.POP20.to_numpy(float), s_units)
    on_blocks = mean_minutes(C, w, s_units)
    print(f"{len(units):4d} {label:12} {C_units.size + C_units.shape[1]:>7,} variables {time.perf_counter() - started:5.2f} s  "
          f"own-level mean {own:.3f}  on blocks {on_blocks:.3f} "
          f"({on_blocks / mean_minutes(C, w, sites) - 1:+.2%}), shared sites {len(set(s_units) & set(sites))} of 8")

  25 junctions    56,000 variables     5.7 s  Optimal     6.07 min
  50 junctions   112,000 variables    11.8 s  Optimal     5.90 min
 100 junctions   224,000 variables    39.0 s  Optimal     5.28 min
 200 junctions   448,000 variables   256.9 s  Optimal     5.10 min
 118 block groups  12,614 variables  0.58 s  own-level mean 5.107  on blocks 5.302 (+0.81%), shared sites 6 of 8
  40 tracts         4,346 variables  0.13 s  own-level mean 4.391  on blocks 5.423 (+3.10%), shared sites 2 of 8

What the model is and why it slows down

What the model is

Each candidate j has a 0/1 variable saying whether it opens, and each demand point i has an assignment variable for each candidate. The constraints say that every demand point is assigned exactly once, only to an open candidate, and that exactly p candidates open. The objective is the sum over pairs of population × travel time × assignment. With no capacities, each demand point ends up assigned to its nearest open site.

Why it gets slow

Both the assignment variables and the "only to an open site" constraints number n × k: 237,440 of each for the county. PuLP builds every one of them as a Python object and writes the model to a file for CBC to read and CBC then has to prove that no other combination of sites is better. Doubling the candidates doubles the model and more than doubles the search.

Why aggregation mostly works

Moving each person to their block group's centroid changes their travel time a little in a dense area and a lot in a large rural one. The p-median objective averages over everyone, so the small errors in dense places dominate and the chosen sites barely change. Tracts are big enough that the errors stop cancelling: their centroids understated the mean trip by about a minute and pulled the sites towards the wrong places.

What the model leaves out

Capacity, the longest trip, and the fact that people do not always use the nearest site. The median model is the right tool when average travel is the honest objective; it is the wrong one when the promise is that nobody travels too far.

Each doubling of the model roughly tripled the solve time, and the last one more than quintupled it.

Units, finite times and forced sites

Units do not change the answer. Seconds or minutes scale the objective, not the choice.
Travel times must be finite. Validate the matrix; see fixing a location-allocation model that never solves for what infinities do.
Ties exist. Different site sets can give the same objective; the solver returns one.
Existing facilities can be forced open by fixing their variables, which also shrinks the search.
Evaluate on the finest demand you have, never on the aggregated units the model was solved on.
Snap distances are not in the matrix. Add them if homes or candidates sit far from the road.
A random candidate list is not a plan. It is useful for testing scale, not for choosing land.

Internal links

Location-allocation explained: p-median, coverage and what each optimises — p-median against other objectives
How to maximise coverage with a limited number of sites — the coverage alternative
Fixing a location-allocation model that never solves — slow and failed solves
How to measure network distance for many origin–destination pairs — building the matrix at scale
How to measure distance to the nearest facility for every home — evaluating existing sites
How to count the population within reach of each facility — per-site demand
How to find the areas nobody can reach within a travel time — where the median model leaves people behind
Census geographies explained: blocks, tracts, output areas and why they nest — the units used for aggregation

FAQ

What does a p-median model optimise?

The population-weighted total, or average, travel time from every demand point to its nearest chosen site. For 8 of 106 school sites in Chittenden County the optimum was a mean drive of 5.26 minutes.

How big a p-median model can I solve?

It depends on the solver and machine, but the size is demand points × candidates. With CBC here, 237,440 variables took 88 seconds, and 448,000 took 257 seconds for one set of candidates and produced no solution within 300 seconds for another.

How do I make a p-median model solve faster?

Aggregate the demand, prune the candidates, or both. Solving on 118 block groups instead of 2,239 blocks took under a second and gave sites only 0.81% worse on the blocks.

Why did spopt raise "Model is not solved: Not Solved"?

The time limit ran out before CBC found any feasible integer solution. spopt checks the PuLP status after solving and raises for anything other than Optimal.

Should I aggregate to census tracts?

Usually not. Tract centroids made the solution look better than it was, 4.39 minutes on its own terms, while its real mean on blocks was 5.42 minutes, the worst of the three levels tested.

Does the p-median model account for facility capacity?

No. Every demand point goes to its nearest open site regardless of load; in the county solution the busiest site served nearly five times as many people as the quietest.