⚙️ How-To Location & accessibility Networks & routing Performance & scale Location & Accessibility · step 13 of 20 Updated

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

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,652 variables, took 101.6 seconds to solve and about 1 GB of memory.
  • With 200 candidates instead of 106, the solver found no solution at all in 5 minutes.
  • Solving on 118 block groups took 0.34 seconds, and its sites were only 0.41% 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.

Flow from demand points and candidate sites to a travel-time matrix, a size check, a time-limited p-median solve and an evaluation of the chosen sites on the full demand.
The size check comes before the solve: the variable count predicts whether the solver will return.

Step-by-step solution

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 intersections drawn at random from the county gave a slightly lower optimal mean, 5.15 minutes against 5.26. 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,241 blocks and 106 schools the matrix took 66 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,241 blocks, with an 8-site budget and random intersections as candidates:

candidates   variables   solve time   status          mean trip
     25         56,050        7.6 s   Optimal          5.62 min
     50        112,100       19.8 s   Optimal          5.24 min
    100        224,200       62.3 s   Optimal          5.15 min
    106 schools 237,652     101.6 s   Optimal          5.26 min
    200        448,400      > 300 s   Not Solved       —

Solve time roughly tripled each time the model doubled, and at 448,400 variables CBC had no integer solution after its 300-second limit. spopt then raised RuntimeError: Model is not solved: Not Solved.

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.73 minutes, the longest trip 26.94 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,241 blocks     237,652   101.6 s      5.261 min       5.261 min            8 of 8
118 groups        12,614     0.34 s     5.100 min       5.282 min (+0.41%)   5 of 8
40 tracts          4,346     0.09 s     4.380 min       5.410 min (+2.83%)   2 of 8

The block-group solution was 300 times faster and only 0.41% worse. The tract solution looked best on its own terms — 4.38 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.

Bar chart of the mean trip on census blocks for p-median sites chosen using blocks, block groups and tracts, with the solve time for each.
Block groups kept almost all of the quality at a fraction of the cost; tracts did not.

Code examples

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

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,241 block nodes and 106 school nodes, it returned a 2,241 × 106 matrix in 66 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,241 demand points x 106 candidates = 237,652 variables
Optimal in 113.3 s: sites [24, 49, 51, 63, 68, 70, 72, 77]
mean 5.26 min, longest 26.94 min, within 10 min 89.0%
people per site: 14,086, 48,014, 28,950, 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 48,014 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, 63, 68, 70, 76, 77, 98]
mean 5.28 min, longest 26.94 min, within 10 min 90.8%
people per site: 14,068, 30,755, 13,952, 15,594, 8,605, 12,846, 15,819, 56,684

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.

Explanation

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,546 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 — peak memory for the run was about 1 GB — 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.

Bar chart of p-median solve time as candidates grow from 25 to 200, ending with no solution found within the 300-second limit.
Each doubling of the model roughly tripled the solve time, until the solver stopped returning.

Edge cases or notes

  • 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.

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, 237,652 variables took 101.6 seconds and 448,400 variables produced no solution within 300 seconds.

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,241 blocks took 0.34 seconds and gave sites only 0.41% 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.38 minutes on its own terms, while its real mean on blocks was 5.41 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.