Location-Allocation Explained: p-Median, Coverage and What Each Optimises
Problem statement
Choosing one site is a ranking problem: score every candidate and take the best. Choosing eight is not, because the best second site depends on where the first one went, and the number of combinations explodes — there are about 3.0 × 10¹¹ ways to choose 8 sites from 106. Location-allocation models solve that combinatorial problem exactly, for a stated objective.
The objective is the whole decision. Measured by choosing 8 of Chittenden County's 106 school sites to serve the 168,323 residents of its 2,241 populated census blocks, with drive times on the road network:
model mean trip 95th percentile longest trip within 10 min solve time
p-median 5.26 min 13.73 min 26.94 min 89.0% 99.7 s
p-centre 6.37 min 13.41 min 24.12 min 84.3% 916.4 s
maximal coverage 6.38 min 10.90 min 29.30 min 93.8% 0.2 s
Each model won the measure it optimises and lost at least one other. The three sets of 8 sites shared at most 2 sites between any pair. And the models differed in solve time by a factor of nearly 5,000.
Quick answer
Pick the model from the promise the service makes:
the service promises... model
the shortest trips on average p-median
nobody travels longer than necessary p-centre
as many people as possible within a standard, with k sites maximal coverage (MCLP)
everyone within a standard, with as few sites as possible set covering (LSCP)
With spopt and PuLP, each is a few lines on a demand × candidate cost matrix:
import pulp
from spopt.locate import LSCP, MCLP, PCenter, PMedian
solver = pulp.PULP_CBC_CMD(msg=False, timeLimit=900)
pmedian = PMedian.from_cost_matrix(C, population, p_facilities=8).solve(solver)
coverage = MCLP.from_cost_matrix(C, population, service_radius=600, p_facilities=8).solve(solver)
chosen = [j for j, v in enumerate(pmedian.fac_vars) if v.value() > 0.5]C holds travel times in seconds from each block to each candidate site; service_radius=600 is ten minutes.
Step-by-step solution
1. Assemble demand, candidates and costs
Demand is population in small units — here 2,241 populated blocks. Candidates are places a facility could go — here 106 school sites, a plausible stand-in for publicly owned land. The cost matrix is the drive time from every demand point to every candidate: 2,241 × 106 entries, with no unreachable pairs. Every model below reads the same three inputs, so any error in them is shared by every answer.
2. Know the best possible service first
With all 106 candidates open, the population-weighted mean trip would be 2.23 minutes and the longest 15.5 minutes, with 99.5% of residents within 10 minutes. Eight sites cannot beat that, and the gap between it and each 8-site solution is the price of having only eight.
3. p-median: minimise total travel
The p-median model chooses p sites and assigns every demand point to its nearest open site, minimising population-weighted travel time. It found the lowest mean trip, 5.26 minutes, and served 89.0% of residents within 10 minutes. Its longest trip was 26.94 minutes: efficiency tolerates a few remote blocks paying heavily.
4. p-centre: minimise the longest trip
The p-centre model minimises the single largest travel time. It cut the longest trip to 24.12 minutes, at the cost of a longer average (6.37 minutes) and the lowest share within 10 minutes (84.3%). It was also by far the hardest to solve: 916 seconds, against 100 for p-median on the same matrix, because a minimax objective gives the solver little to prune with.
5. Maximal coverage: maximise people within a standard
The maximal covering location problem (MCLP) chooses p sites to maximise the population within a service radius, ignoring how far inside or outside it anyone is. With a 10-minute radius it covered 93.8% of residents — the most of any model — and had the best 95th percentile. It also had the worst longest trip, 29.30 minutes, because anyone already beyond reach contributes nothing to its objective. It solved in 0.19 seconds.
6. Set covering: the fewest sites to cover everyone
The location set covering problem (LSCP) turns the question round: what is the smallest number of sites that puts everyone within the standard? At 10 minutes it is infeasible here: the solver raised
RuntimeError: Model is not solved: Infeasible. See ``pulp.constants.LpStatus`` for more information.
because some blocks are 15.5 minutes from every candidate school. Restricted to the blocks that some school can reach within 10 minutes — 99.5% of residents — it needed 26 sites; at a 5-minute standard, 54 sites to cover the 92.0% who can be covered at all.
7. Expect diminishing returns
Adding sites one at a time under maximal coverage with a 10-minute standard:
sites 1 2 3 4 5 6 8 10 12
covered 50.6% 66.3% 73.4% 79.6% 84.6% 88.4% 93.8% 96.0% 97.9%
The first site covered half the county; the twelfth added 0.9 percentage points. The curve is the argument for — or against — a budget.
8. Compare the solutions, not only the objectives
The three 8-site answers shared almost nothing: p-median and p-centre had 1 site in common, p-median and maximal coverage 1, p-centre and maximal coverage 2. Evaluate every solution on every measure, as in the table above, and show the trade-off to whoever decides.
Code examples
Example 1 — evaluate any set of sites on every objective
import numpy as np
def weighted_quantile(values, weights, q):
order = np.argsort(values)
cum = np.cumsum(weights[order]) / weights.sum()
return float(values[order][np.searchsorted(cum, q)])
def evaluate(C, w, sites, radius=600):
"""Population-weighted outcomes when demand uses its nearest chosen site (C in seconds)."""
t = C[:, list(sites)].min(axis=1)
return {
"mean_min": (t * w).sum() / w.sum() / 60,
"p95_min": weighted_quantile(t, w, 0.95) / 60,
"max_min": t.max() / 60,
"within_r": w[t <= radius].sum() / w.sum(),
}
def chosen(model):
return sorted(j for j, v in enumerate(model.fac_vars) if v.value() is not None and v.value() > 0.5)Scoring every model's sites with the same function is what makes the comparison fair. Each model's own objective value is not comparable across models.
Example 2 — solve and compare four models
import time
import pulp
from spopt.locate import LSCP, MCLP, PCenter, PMedian
def run(name, build, C, w, limit=900, radius=600):
started = time.perf_counter()
try:
model = build().solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=limit))
except RuntimeError as err:
print(f"{name}: {err}")
return None
sites = chosen(model)
e = evaluate(C, w, sites, radius)
print(f"{name}: {pulp.LpStatus[model.problem.status]}, {time.perf_counter() - started:.2f}s, {len(sites)} sites | "
f"mean {e['mean_min']:.2f} | p95 {e['p95_min']:.2f} | max {e['max_min']:.2f} | "
f"within {radius // 60} min {e['within_r']:.1%}")
return sites
P, R = 8, 600
pm = run("p-median", lambda: PMedian.from_cost_matrix(C, w, p_facilities=P), C, w)
pc = run("p-centre", lambda: PCenter.from_cost_matrix(C, p_facilities=P), C, w)
mc = run("MCLP", lambda: MCLP.from_cost_matrix(C, w, service_radius=R, p_facilities=P), C, w)
ls = run("LSCP", lambda: LSCP.from_cost_matrix(C, service_radius=R), C, w, limit=300)On the county's matrix the solver itself took 99.7 s for p-median, 916.4 s for p-centre and 0.19 s for MCLP; the printed times also include building each model, 104.86 s and 0.70 s for p-median and MCLP in one run. The outcomes are those in the table at the top. The set-covering call printed LSCP: Model is not solved: Infeasible. followed by spopt's pointer to pulp.constants.LpStatus, rather than returning sites, because spopt raises when the solver reports a status other than optimal.
Example 3 — a coverage curve, optimal against greedy
def greedy_coverage(cover, w, p):
"""Add the site that covers the most still-uncovered people, p times."""
covered = np.zeros(cover.shape[0], dtype=bool)
shares = []
for _ in range(p):
gain = (cover & ~covered[:, None]).T @ w
covered |= cover[:, int(np.argmax(gain))]
shares.append(w[covered].sum() / w.sum())
return shares
cover = C <= 600
greedy = greedy_coverage(cover, w, 12)
for p in range(1, 13):
m = MCLP.from_cost_matrix(C, w, service_radius=600, p_facilities=p).solve(pulp.PULP_CBC_CMD(msg=False))
optimal = w[cover[:, chosen(m)].any(axis=1)].sum() / w.sum()
print(f"p={p:2d} optimal {optimal:.1%} greedy {greedy[p - 1]:.1%} gap {100 * (optimal - greedy[p - 1]):.2f} pts")At a 10-minute radius the greedy heuristic trailed the optimum by 5.52 points with 2 sites and by 2.47 points with 8. The optimal 2-site solution did not contain the optimal single site at all — which is exactly what a greedy method cannot discover.
Explanation
Why objectives pull sites apart
Minimising total travel is dominated by where most people live, so p-median sites sit near population centres. Minimising the longest trip is dominated by the most remote block, so p-centre spreads sites towards the edges. Maximising coverage is dominated by the population just outside the radius of other sites, and it abandons anyone who cannot be brought inside. Each objective is a value judgement expressed as arithmetic; there is no neutral choice.
Why optimal sets are not nested
The best single site sits in the densest part of the county. With two sites, it is better to split that area between two sites slightly off-centre than to keep the first and add a second elsewhere. So the optimal k-site solution is not the optimal (k − 1)-site solution plus one. Greedy methods, which never undo a choice, cannot find those rearrangements, which is why they trailed the optimum by several points.
Why solve times differ so much
Coverage and set-covering models only need to know whether each block is within the radius of each site, so their constraint matrix is sparse and their relaxations are strong. The p-median model has an assignment variable for every block–site pair — 237,546 here — and p-centre adds a minimax objective whose relaxation gives little guidance, so the solver searches far longer to prove optimality. The same matrix produced solve times of 0.19 s, 99.7 s and 916.4 s.
Why infeasibility is information
A set-covering model that cannot be solved is saying something true: no choice of candidates puts everyone within the standard. The useful response is not to force a solution but to report the smallest feasible standard — 15.5 minutes here — or the share of demand that can be covered, and to decide whether to add candidate sites or relax the promise.
Edge cases or notes
- A time-limited run can still say Optimal. PuLP's CBC interface labels a run that stopped on time with a solution in hand as
Optimal; onlyproblem.sol_status(2, integer feasible, rather than 1) shows that optimality was never proved. - Candidate lists drive everything. A model can only choose from the sites it is offered.
- Demand aggregation changes answers. Blocks, block groups and tracts can choose different sites; test on the finest data you have.
- Capacity is not modelled here. Real facilities fill up; capacitated variants exist and are harder to solve.
- Existing facilities can be fixed. Force current sites open and let the model add to them.
- Infeasible set covering raises in spopt. Catch the error and report the smallest feasible radius.
- Equity and efficiency trade off. Show p-median and p-centre side by side rather than choosing silently.
Internal links
- How to choose facility locations with a p-median model in Python — the p-median workflow in full
- How to maximise coverage with a limited number of sites — coverage models in practice
- Fixing a location-allocation model that never solves — when the solver runs forever or reports infeasibility
- Location analysis explained: catchments, accessibility and site selection — where allocation fits
- How to find the areas nobody can reach within a travel time — measuring coverage before optimising it
- How to measure network distance for many origin–destination pairs — building the cost matrix
- Site suitability explained: constraints, factors and weights — producing candidate sites
- Accessibility measures explained: nearest, cumulative and gravity — evaluating the outcome
FAQ
What is location-allocation?
Choosing a set of facility sites and assigning demand to them together, so that a stated objective — total travel, the longest trip, or population within a standard — is optimised for the whole set rather than one site at a time.
What is the difference between p-median and p-centre?
p-median minimises total or average travel; p-centre minimises the longest trip. For 8 of 106 school sites, p-median gave a mean trip of 5.26 minutes and a longest of 26.94, p-centre 6.37 and 24.12.
When should I use maximal coverage instead?
When the service has a standard, such as ten minutes, and a fixed number of sites. It covered 93.8% of residents within ten minutes with 8 sites, more than either median or centre models, but left its worst-served residents further away.
Why did my set covering model fail?
Because no set of candidates can put every demand point within the radius. Here the smallest feasible standard was 15.5 minutes, so a 10-minute model was infeasible.
Is a greedy approach good enough?
It is a reasonable start but not optimal. At a 10-minute radius greedy selection trailed the optimum by up to 5.52 percentage points, because the best sets of sites are not built by adding one site at a time.
How long do these models take to solve?
It depends on the model more than the data. On the same 2,241 × 106 matrix, maximal coverage took 0.19 s, p-median 99.7 s and p-centre 916.4 s.