Location-Allocation Explained: p-Median, Coverage and What Each Optimises

On this page

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 167,770 residents of the 2,239 populated census blocks with a road route to a school, 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.72 min        26.99 min         89.0%          87.6 s
p-centre                7.74 min      13.92 min        29.69 min         75.1%         919.2 s (at its 900 s limit)
maximal coverage        6.38 min      10.95 min        29.36 min         93.8%           0.7 s

p-median and maximal coverage each won the measure they optimise. The p-centre run won nothing: it stopped at its 900-second time limit, PuLP still reported Optimal, and its longest trip — the one number it exists to minimise — was 2.7 minutes worse than p-median's. The three sets of 8 sites shared at most 2 sites between any pair, and the solve times differed by a factor of more than 1,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.

Two models fix the number of sites and optimise service; one fixes the service standard and minimises sites.

Demand, candidates, costs, then four objectives

1. Assemble demand, candidates and costs

Demand is population in small units — here 2,239 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,239 × 106 entries, after dropping the two blocks whose nearest road node reaches no school at all. 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 88.8% of residents within 10 minutes. Its longest trip was 26.99 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, and it is by far the hardest to solve, because a minimax objective gives the solver little to prune with. Here CBC was still searching when its 900-second limit arrived. PuLP reported Optimal anyway — the status it gives whenever the solver stops holding a feasible solution — and the solution it held had a longest trip of 29.69 minutes, worse than p-median's 26.99, with a mean of 7.74 minutes and only 75.1% of residents within 10. The clock ends that search, so the sites vary between runs. A p-centre answer that loses to p-median on the longest trip is not an optimum; check model.problem.sol_status (1 is proved, 2 is not) and give the solver longer, or a smaller candidate list, before reporting one.

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. Its longest trip was 29.36 minutes, longer than p-median's, because anyone already beyond reach contributes nothing to its objective. It solved in under a second.

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.4%  73.5%  79.7%  84.5%  88.3%  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 no 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.

Half the county from one site, and under one percentage point from the twelfth.

Four location models solved and compared

The data used on this page

Chittenden County, Vermont, from the 2020 Census: TIGER/Line 2020 blocks for Vermont, with each populated block reduced to a representative point carrying its POP20. The candidates are every OpenStreetMap amenity=school feature in the county — 106 of them. 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 the nearest nodes. Two blocks (553 people) whose nearest node has no road route to any school are dropped here, before any model sees them — the fix the troubleshooting guide describes. Solved with spopt 0.7 and PuLP 3.3's bundled CBC; spopt fails on PuLP 4.0, so install spopt "pulp<4". Overpass via the maps.mail.ru mirror while overpass-api.de was refusing connections. The OpenStreetMap answers are read as the Overpass mirror gave them on 10 October 2026.

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

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

# OpenStreetMap as Overpass returned it on 10 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-location-allocation-explained-2026-10-10.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
homes = vt_blocks[(vt_blocks.COUNTYFP20 == "007") & (vt_blocks.POP20 > 0)].copy()
homes["geometry"] = homes.geometry.representative_point()
homes = homes.reset_index(drop=True)

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)
index = {n: i for i, n in enumerate(G.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)
T = sp.csr_array((list(best.values()), (rows, cols)), shape=(len(index), len(index)))
s_idx = [index[n] for n in ox.distance.nearest_nodes(G, schools.geometry.x, schools.geometry.y)]
h_all = np.array([index[n] for n in ox.distance.nearest_nodes(G, homes.geometry.x, homes.geometry.y)])
C_all = dijkstra(T.T, directed=True, indices=s_idx)[:, h_all].T

reachable = np.isfinite(C_all).all(axis=1)
print(f"dropping {(~reachable).sum()} blocks ({homes.POP20[~reachable].sum():,} people) with no road route")
homes, h_idx = homes[reachable].reset_index(drop=True), list(h_all[reachable])
C, w = C_all[reachable], homes["POP20"].to_numpy(float)
print(f"C: {C.shape[0]:,} blocks x {C.shape[1]} schools, {w.sum():,.0f} residents")

dropping 2 blocks (553 people) with no road route
C: 2,239 blocks x 106 schools, 167,770 residents

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)
print(f"sites shared: p-median/p-centre {len(set(pm) & set(pc))}, p-median/MCLP {len(set(pm) & set(mc))}, "
      f"p-centre/MCLP {len(set(pc) & set(mc))}")

p-median: Optimal, 87.64s, 8 sites | mean 5.26 | p95 13.72 | max 26.99 | within 10 min 89.0%
p-centre: Optimal, 919.20s, 8 sites | mean 7.74 | p95 13.92 | max 29.69 | within 10 min 75.1%
MCLP: Optimal, 0.65s, 8 sites | mean 6.38 | p95 10.95 | max 29.36 | within 10 min 93.8%
LSCP: Model is not solved: Infeasible. See ``pulp.constants.LpStatus`` for more information.
sites shared: p-median/p-centre 0, p-median/MCLP 1, p-centre/MCLP 2

The printed times include building each model: 87.6 s for p-median, 919.2 s for p-centre — stopped by its limit, whatever the status says — and 0.65 s for MCLP. 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.

The p-centre run here passed the first check and failed the third.

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")

p= 1  optimal 50.6%  greedy 50.6%  gap 0.00 pts
p= 2  optimal 66.4%  greedy 60.9%  gap 5.54 pts
p= 3  optimal 73.5%  greedy 68.6%  gap 4.84 pts
p= 4  optimal 79.7%  greedy 75.6%  gap 4.17 pts
p= 5  optimal 84.5%  greedy 81.8%  gap 2.70 pts
p= 6  optimal 88.3%  greedy 87.1%  gap 1.19 pts
p= 7  optimal 91.5%  greedy 89.6%  gap 1.88 pts
p= 8  optimal 93.8%  greedy 91.3%  gap 2.48 pts
p= 9  optimal 95.0%  greedy 92.7%  gap 2.24 pts
p=10  optimal 96.0%  greedy 94.0%  gap 2.05 pts
p=11  optimal 97.1%  greedy 95.1%  gap 1.97 pts
p=12  optimal 97.9%  greedy 96.1%  gap 1.81 pts

At a 10-minute radius the greedy heuristic trailed the optimum by 5.58 points with 2 sites and by 2.48 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.

Example 4 — the other figures quoted above

from math import comb

every = evaluate(C, w, range(C.shape[1]))
print(f"combinations of 8 from {C.shape[1]}: {comb(C.shape[1], 8):.1e}")
print(f"all {C.shape[1]} open: mean {every['mean_min']:.2f} min, longest {every['max_min']:.1f} min, "
      f"within 10 min {every['within_r']:.1%}")
for radius in (600, 300):
    ok = C.min(axis=1) <= radius
    m = LSCP.from_cost_matrix(C[ok], service_radius=radius).solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=300))
    print(f"set covering at {radius // 60} min over the {w[ok].sum() / w.sum():.1%} who can be covered: "
          f"{len(chosen(m))} sites")

combinations of 8 from 106: 3.0e+11
all 106 open: mean 2.23 min, longest 15.5 min, within 10 min 99.5%
set covering at 10 min over the 99.5% who can be covered: 26 sites
set covering at 5 min over the 92.0% who can be covered: 54 sites

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,440 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 times of 0.65 s, 101.7 s and — without proving anything — 895.6 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.

No model wins every column, and a time-limited p-centre can lose even its own.

Time limits, candidate lists and capacity

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; only problem.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.99. The p-centre run stopped at its 900-second limit with a longest trip of 29.69 — a reminder that a time-limited minimax solve can report Optimal without being it.

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.58 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,239 × 106 matrix, maximal coverage took 0.65 s, p-median 87.6 s and p-centre hit a 900-second limit without proving an optimum.