How to Maximise Coverage with a Limited Number of Sites

On this page

Problem statement

A service with a standard — an ambulance within ten minutes, a clinic within a short drive — and a budget for a fixed number of sites needs the maximal covering location problem (MCLP): choose the sites that put the most people within the standard. The model itself is quick to solve. The mistakes happen around it.

Measured by choosing 8 of Chittenden County's 106 school sites for a 10-minute drive standard, serving the 168,323 residents of 2,241 census blocks:

The best 8 sites covered 93.4% of residents, against a ceiling of 99.2% with every site open.
spopt's own coverage figure said 85.2%. Its perc_cov attribute counts demand points, not people.
The radius chose the sites. A plan optimised for 15 minutes shared only 1 or 2 of its 8 sites with the 10-minute plan and put just 65.2% to 68.3% of residents within 10 minutes.
Committing two sites in advance cost 2.5 percentage points. Forcing even the best single site into the plan cost 1.5.

Quick answer

Build a demand × candidate travel-time matrix, pass population as the weights, and compute coverage from the chosen sites yourself:

import pulp
from spopt.locate import MCLP

model = MCLP.from_cost_matrix(C, w, service_radius=600, p_facilities=8)
model.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=300))
sites = [j for j, v in enumerate(model.fac_vars) if v.value() > 0.5]
covered = w[(C[:, sites] <= 600).any(axis=1)].sum() / w.sum()
print(f"{len(sites)} sites cover {covered:.1%} of residents within 10 minutes")

8 sites cover 93.4% of residents within 10 minutes

C holds drive times in seconds from each block to each candidate, and w holds each block's population. spopt 0.7 needs PuLP 3 — pip install spopt "pulp<4" — because PuLP 4.0 changed the LpVariable signature and spopt's model builder fails with TypeError: ... unexpected keyword argument 'lowBound'. Then solve again at other radii and other numbers of sites before recommending anything.

The solve is one step of five; the other four decide whether its answer can be trusted.

A standard, a ceiling, weights, then solve

1. Fix the standard and build the matrix

Coverage is binary: a block is covered when its travel time to some chosen site is within the radius. Measure the radius in the units of the matrix — 600 seconds here — and use network travel time rather than straight-line distance if the standard is about travel. The 2,241 × 106 matrix of drive times took well under a second to build from one reversed shortest-path search per candidate.

2. Find the ceiling before optimising

With all 106 candidates open, 99.2% of residents were within 10 minutes of one, and 91.7% within 5 minutes. No choice of sites can beat those figures, and the blocks beyond them — 15.5 minutes from the nearest candidate at worst — need a different candidate list, not a better model. Two blocks, 553 people, had no road route to any candidate at all: their nearest node sits on a fragment of the graph. They count as uncovered in every figure here.

3. Weight by population

Without weights, MCLP maximises the number of demand points covered. On the county's blocks the unweighted plan covered 86.4% of blocks but 92.5% of residents; the weighted plan covered 85.2% of blocks and 93.4% of residents. They shared 5 of 8 sites. Pass the quantity you care about — population, households, calls — as weights.

4. Compute coverage in people yourself

After solving, spopt sets perc_cov and n_cli_uncov. Both count demand points: the weighted 8-site plan reported perc_cov 85.2 and 332 uncovered blocks. Those 332 blocks held 11,053 people, 6.6% of residents, so the plan covered 93.4% of the population. Report the population share, and state what the demand points are.

5. Solve for a range of site counts

Coverage rises quickly and then flattens:

sites      1      2      3      4      5      6      7      8      9     10     11     12
covered  50.4%  66.2%  73.2%  79.5%  84.2%  88.0%  91.2%  93.4%  94.7%  95.7%  96.8%  97.6%
gain     50.4   15.8    7.0    6.2    4.8    3.8    3.2    2.3    1.2    1.0    1.0    0.9

Each solve took about a second (Example 1). The curve is what a budget discussion needs: the ninth site bought 1.2 points, the twelfth 0.9. Covering everyone who can be covered is a different model — set covering — and needed 26 sites (Example 4).

6. Test the radius

The same 8-site budget optimised for different radii, each solution then measured against the 10-minute standard:

optimised for     own radius   within 10 min   sites shared with the 10-min plan
 5 min               65.5%          85.4%              0 of 8
 8 min               86.6%          91.5%              2 of 8
10 min               93.4%          93.4%              8 of 8
12 min               97.0%     72.6% or 86.0%          2 of 8
15 min               99.2%     68.3% or 65.2%       2 or 1 of 8

The 12-minute and 15-minute rows have two answers because the model has tied optima. Two runs returned different sets of 8 sites, both covering exactly 97.0% within 12 minutes, one putting 86.0% of residents within 10 minutes and the other 72.6%; at 15 minutes both covered 99.2% and put 68.3% or 65.2% within 10. The solver is not wrong — the objective cannot tell them apart — but a plan chosen from a single solve can be one of several equally good plans that differ on everything the objective ignores.

A longer radius lets sites spread out to reach remote blocks, and the dense centre then falls outside the tighter standard. If the standard is uncertain, choose sites that do well across the plausible radii, not the optimum for one.

Tied optima: equal on the objective, 13 points apart on the measure a reader would check next.

7. Fix sites that are already committed

Existing facilities, or sites already promised, go in through predefined_facilities_arr, a 0/1 array over candidates. Committing candidates 49 and 86 — the best pair for a 5-minute standard, Mater Christi School in Burlington and South Burlington High School — and letting the model choose 6 more covered 90.9% within 10 minutes, 2.5 points below the free optimum. Committing candidate 47, the best single site, covered 91.9%: the best site to open first is not part of the best set of 8.

8. Compare with a greedy plan

Adding at each step the site that covers the most still-uncovered people is easy to explain but not optimal. At 10 minutes it trailed the optimum by 2.45 points with 8 sites and by 5.53 points with 2. The optimal plans are not nested — the best 2 sites did not include the best single site — which is exactly what one-at-a-time selection cannot discover.

A plan optimised for 15 minutes leaves about a third of residents beyond 10.

Coverage curves and how the radius changes the plan

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, points and outlines, by representative point, in a fixed order. Drive times come from OpenStreetMap's drive network for the county plus 10.5 km, with OSMnx's imputed free-flow speeds, and are road-only: from each block's nearest node to each school's. 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 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-maximal-coverage-location-python-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)
w = homes["POP20"].to_numpy(float)

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)))
h_idx = [index[n] for n in ox.distance.nearest_nodes(G, homes.geometry.x, homes.geometry.y)]
s_idx = [index[n] for n in ox.distance.nearest_nodes(G, schools.geometry.x, schools.geometry.y)]

t = time.perf_counter()
C = dijkstra(T.T, directed=True, indices=s_idx)[:, h_idx].T      # blocks x candidates, seconds
print(f"{w.sum():,.0f} residents in {len(homes):,} blocks; {len(schools)} candidate schools; "
      f"matrix {C.shape[0]:,} x {C.shape[1]} in {time.perf_counter() - t:.2f} s")
nearest = C.min(axis=1)
stranded = ~np.isfinite(nearest)
print(f"within 10 min of any candidate: {w[nearest <= 600].sum() / w.sum():.1%}; "
      f"within 5 min: {w[nearest <= 300].sum() / w.sum():.1%}; "
      f"worst reachable block {nearest[~stranded].max() / 60:.1f} min; "
      f"blocks with no road route to any candidate: {stranded.sum()} ({w[stranded].sum():,.0f} people)")

168,323 residents in 2,241 blocks; 106 candidate schools; matrix 2,241 x 106 in 0.07 s
within 10 min of any candidate: 99.2%; within 5 min: 91.7%; worst reachable block 15.5 min; blocks with no road route to any candidate: 2 (553 people)

The examples share two helpers:

import numpy as np
import pulp
from spopt.locate import MCLP


def covered_share(C, w, sites, radius):
    """Share of the weight within radius of at least one chosen site."""
    return w[(C[:, sites] <= radius).any(axis=1)].sum() / w.sum()


def solve_mclp(C, w, radius, p, fixed=None):
    committed = None
    if fixed:
        committed = np.zeros(C.shape[1], dtype=int)
        committed[fixed] = 1
    model = MCLP.from_cost_matrix(C, w, service_radius=radius, p_facilities=p,
                                  predefined_facilities_arr=committed)
    model.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=300))
    return model, [j for j, v in enumerate(model.fac_vars) if v.value() > 0.5]

Example 1 — a coverage curve with its ceiling

def coverage_curve(C, w, radius, max_sites):
    ceiling = covered_share(C, w, list(range(C.shape[1])), radius)
    print(f"all {C.shape[1]} candidates open: {ceiling:.1%}")
    previous = 0.0
    for p in range(1, max_sites + 1):
        _, sites = solve_mclp(C, w, radius, p)
        share = covered_share(C, w, sites, radius)
        print(f"p={p:2d}  covered {share:.1%}  gain {100 * (share - previous):4.1f} pts")
        previous = share


coverage_curve(C, w, radius=600, max_sites=12)

all 106 candidates open: 99.2%
p= 1  covered 50.4%  gain 50.4 pts
p= 2  covered 66.2%  gain 15.8 pts
p= 3  covered 73.2%  gain  7.0 pts
p= 4  covered 79.5%  gain  6.2 pts
p= 5  covered 84.2%  gain  4.8 pts
p= 6  covered 88.0%  gain  3.8 pts
p= 7  covered 91.2%  gain  3.2 pts
p= 8  covered 93.4%  gain  2.3 pts
p= 9  covered 94.7%  gain  1.2 pts
p=10  covered 95.7%  gain  1.0 pts
p=11  covered 96.8%  gain  1.0 pts
p=12  covered 97.6%  gain  0.9 pts

Example 2 — how much the radius changes the plan

def radius_sensitivity(C, w, p, radii, standard):
    _, reference = solve_mclp(C, w, standard, p)
    for radius in radii:
        _, sites = solve_mclp(C, w, radius, p)
        print(f"optimised for {radius // 60:2d} min: {covered_share(C, w, sites, radius):.1%} within {radius // 60} min, "
              f"{covered_share(C, w, sites, standard):.1%} within the {standard // 60}-min standard, "
              f"{len(set(sites) & set(reference))} of {p} sites shared")


radius_sensitivity(C, w, p=8, radii=[300, 480, 600, 720, 900], standard=600)

optimised for  5 min: 65.5% within 5 min, 85.4% within the 10-min standard, 0 of 8 sites shared
optimised for  8 min: 86.6% within 8 min, 91.5% within the 10-min standard, 2 of 8 sites shared
optimised for 10 min: 93.4% within 10 min, 93.4% within the 10-min standard, 8 of 8 sites shared
optimised for 12 min: 97.0% within 12 min, 86.0% within the 10-min standard, 2 of 8 sites shared
optimised for 15 min: 99.2% within 15 min, 65.2% within the 10-min standard, 1 of 8 sites shared

Example 3 — committed sites, and the figure not to report

_, pair = solve_mclp(C, w, radius=300, p=2)              # the best pair for a 5-minute standard
_, (best_one,) = solve_mclp(C, w, radius=600, p=1)       # the best single site at 10 minutes
free, free_sites = solve_mclp(C, w, radius=600, p=8)
fixed, fixed_sites = solve_mclp(C, w, radius=600, p=8, fixed=pair)
one, one_sites = solve_mclp(C, w, radius=600, p=8, fixed=[best_one])

uncovered = ~(C[:, free_sites] <= 600).any(axis=1)
print(f"spopt perc_cov: {free.perc_cov:.1f}% of demand points ({free.n_cli_uncov} uncovered)")
print(f"those {uncovered.sum()} blocks hold {w[uncovered].sum():,.0f} people ({w[uncovered].sum() / w.sum():.1%})")
print(f"committed pair {pair}: {', '.join(schools.loc[pair, 'name'])}")
print(f"population covered: {covered_share(C, w, free_sites, 600):.1%} free, "
      f"{covered_share(C, w, fixed_sites, 600):.1%} with the pair committed "
      f"({len(set(free_sites) & set(fixed_sites))} sites in common), "
      f"{covered_share(C, w, one_sites, 600):.1%} with site {best_one} committed; "
      f"site {best_one} in the free plan: {best_one in free_sites}")
trip = C[:, free_sites].min(axis=1)
print(f"longest trip in the free plan: {trip[np.isfinite(trip)].max() / 60:.2f} min")

spopt perc_cov: 85.2% of demand points (331 uncovered)
those 331 blocks hold 11,026 people (6.6%)
committed pair [49, 86]: Mater Christi School, South Burlington High School
population covered: 93.4% free, 90.9% with the pair committed (4 sites in common), 91.9% with site 47 committed; site 47 in the free plan: False
longest trip in the free plan: 29.36 min

Example 4 — weights, set covering and a greedy plan

from spopt.locate import LSCP

# unweighted: every block counts once
_, unweighted = solve_mclp(C, np.ones(len(w)), 600, 8)
_, weighted = solve_mclp(C, w, 600, 8)
blocks = np.ones(len(w))
for label, sites in (("unweighted", unweighted), ("weighted", weighted)):
    print(f"{label:10}: {covered_share(C, blocks, sites, 600):.1%} of blocks, "
          f"{covered_share(C, w, sites, 600):.1%} of residents")
print(f"sites shared: {len(set(unweighted) & set(weighted))} of 8")

# covering everyone who can be covered
coverable = C.min(axis=1) <= 600
lscp = LSCP.from_cost_matrix(C[coverable], service_radius=600)
lscp.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=300))
print(f"set covering: {sum(v.value() > 0.5 for v in lscp.fac_vars)} sites to cover all "
      f"{coverable.sum():,} coverable blocks")


def greedy(C, w, radius, p):
    """Add, one at a time, the site that covers the most still-uncovered weight."""
    cover, done, chosen = C <= radius, np.zeros(len(w), bool), []
    for _ in range(p):
        j = int(np.argmax((cover & ~done[:, None]).T @ w))
        chosen.append(j)
        done |= cover[:, j]
    return chosen


print(f"{(C <= 600).mean():.1%} of block-school pairs are within 10 minutes")
_, wrong = solve_mclp(C, w, 10, 8)                       # 10 meant as minutes, read as seconds
print(f"service_radius=10 on a matrix in seconds: {covered_share(C, w, wrong, 10):.1%} covered")

for p in (1, 2, 8):
    _, best = solve_mclp(C, w, 600, p)
    g = greedy(C, w, 600, p)
    print(f"p={p}: optimal {covered_share(C, w, best, 600):.2%}, greedy {covered_share(C, w, g, 600):.2%}, "
          f"gap {100 * (covered_share(C, w, best, 600) - covered_share(C, w, g, 600)):.2f} pts; optimal sites {sorted(best)}")

unweighted: 86.4% of blocks, 92.5% of residents
weighted  : 85.2% of blocks, 93.4% of residents
sites shared: 5 of 8
set covering: 26 sites to cover all 2,215 coverable blocks
16.7% of block-school pairs are within 10 minutes
service_radius=10 on a matrix in seconds: 1.2% covered
p=1: optimal 50.42%, greedy 50.42%, gap 0.00 pts; optimal sites [47]
p=2: optimal 66.20%, greedy 60.68%, gap 5.53 pts; optimal sites [34, 42]
p=8: optimal 93.45%, greedy 90.98%, gap 2.47 pts; optimal sites [15, 24, 36, 42, 50, 62, 69, 104]

Why coverage models solve quickly

What the model maximises

MCLP has a 0/1 variable for each candidate and one for each demand point. A demand point may count as covered only if some open candidate is within the radius, exactly p candidates may open, and the objective is the weighted sum of covered demand points. It never looks at how far inside the radius anyone is, or how far outside.

Why coverage models solve quickly

The only thing the model needs from the matrix is which pairs are within the radius. Here 16.7% of block–school pairs were within 10 minutes, so each coverage constraint involves few variables, and the linear relaxation is close to the integer answer. Each of the solves above took under a second; a p-median model on the same matrix, with a variable for every block–school pair, is far slower.

Why the radius moves every site

A block counts fully when it is 1 second inside the radius and not at all when 1 second outside. Changing the radius changes which blocks are just inside or just outside each site's reach, so the best combination can change completely — as it did between 5 and 10 minutes, with no sites in common.

What coverage leaves out

Everyone beyond the radius is worth nothing to the objective, whether they are 11 minutes away or 29. The 10-minute plan's longest trip was 29.36 minutes, the worst of the three model types compared in location-allocation explained. If the far tail matters, report it next to the coverage share, or constrain it.

Weights decide what is being maximised; perc_cov reports the unweighted figure either way.

Binary coverage, ties and radius units

Coverage is binary. Partial or distance-decayed coverage needs a different model.
Backup coverage is separate. A block covered by one site loses service when that site is busy; spopt's LSCPB models backup coverage for set covering.
Ties are common. Several site sets can cover the same population; the solver returns one.
Solutions depend on the candidate list. Blocks no candidate can reach stay uncovered whatever the budget.
Capacity is ignored. A site covering 20,000 people counts the same as one covering 200.
A radius in the wrong unit is silent. A matrix in seconds with service_radius=10 still solves; it covered 1.2% of residents, the blocks within 10 seconds of a site.
Solutions are not nested. Do not phase a rollout by taking the first k sites of a larger plan without checking coverage at each phase.

Internal links

Location-allocation explained: p-median, coverage and what each optimises — how coverage compares with other objectives
How to choose facility locations with a p-median model in Python — minimising average travel instead
Fixing a location-allocation model that never solves — infeasible covering models and slow solves
How to find the areas nobody can reach within a travel time — measuring coverage for existing sites
How to measure network distance for many origin–destination pairs — building the matrix
How to count the population within reach of each facility — the coverage of each site on its own
Accessibility measures explained: nearest, cumulative and gravity — measures that are not binary
How to measure distance to the nearest facility for every home — nearest-site travel times

FAQ

What is the maximal covering location problem?

A model that chooses a fixed number of sites to maximise the demand within a service standard, such as a travel time. With 8 of 106 school sites and a 10-minute drive, the best plan covered 93.4% of the county's residents.

Why is spopt's perc_cov lower than my population coverage?

Because perc_cov is the percentage of demand points covered, not of the weights. The plan above reported 85.2%, the share of blocks; weighted by population it covered 93.4%.

Do I need to weight the demand points?

Yes, unless every point represents the same demand. An unweighted plan covered more blocks, 86.4%, but fewer residents, 92.5%, and shared only 5 of 8 sites with the weighted plan.

How sensitive are the sites to the service radius?

Very. The 5-minute and 10-minute plans had no sites in common, and a plan optimised for 15 minutes put only 65–68% of residents within 10, depending on which of its tied optima the solver returned.

How do I keep existing facilities in the solution?

Pass a 0/1 array over the candidates as predefined_facilities_arr. Committing two central sites reduced the best 8-site coverage from 93.4% to 90.9%.

How many sites do I need to cover everyone?

Use a set covering model instead, restricted to demand that some candidate can reach. Within 10 minutes that took 26 sites for 99.2% of residents; the rest were beyond every candidate.