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

How to Maximise Coverage with a Limited Number of Sites

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.8% of residents, against a ceiling of 99.5% with every site open.
  • spopt's own coverage figure said 85.3%. Its perc_cov attribute counts demand points, not people.
  • The radius chose the sites. A plan optimised for 15 minutes shared only 2 of its 8 sites with the 10-minute plan and put just 68.7% of residents within 10 minutes.
  • Committing two sites in advance cost 2.6 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()

C holds drive times in seconds from each block to each candidate, and w holds each block's population. Then solve again at other radii and other numbers of sites before recommending anything.

Flow from a travel-time matrix to the coverage ceiling, maximal coverage solved for a range of site counts, coverage measured in people, and the radius tested.
The solve is one step of five; the other four decide whether its answer can be trusted.

Step-by-step solution

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.5% of residents were within 10 minutes of one, and 92.0% 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.

3. Weight by population

Without weights, MCLP maximises the number of demand points covered. On the county's blocks the unweighted plan covered 86.5% of blocks but 92.8% of residents; the weighted plan covered 85.3% of blocks and 93.8% 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.32 and 329 uncovered blocks. Those 329 blocks held 10,473 people, 6.2% of residents, so the plan covered 93.8% 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.6%  66.3%  73.4%  79.6%  84.6%  88.4%  91.5%  93.8%  95.0%  96.0%  97.1%  97.9%
gain     50.6   15.8    7.0    6.2    5.0    3.8    3.2    2.3    1.2    1.0    1.0    0.9

Each solve took under a second. 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.

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.9%          85.8%              0 of 8
 8 min               86.9%          91.9%              2 of 8
10 min               93.8%          93.8%              8 of 8
12 min               97.4%          86.3%              2 of 8
15 min               99.6%          68.7%              2 of 8

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.

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, both in the dense centre — and letting the model choose 6 more covered 91.2% within 10 minutes, 2.6 points below the free optimum. Committing candidate 47, the best single site, covered 92.3%: 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.47 points with 8 sites and by 5.52 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.

Bar chart of the share of residents within 10 minutes for 8-site plans optimised for radii of 5, 8, 10, 12 and 15 minutes.
A plan optimised for 15 minutes leaves nearly a third of residents beyond 10.

Code examples

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)

It printed all 106 candidates open: 99.5%, then the curve in step 5, from p= 1 covered 50.6% gain 50.6 pts to p=12 covered 97.9% gain 0.9 pts.

Example 2 — how much the radius changes the plan 

def radius_sensitivity(C, w, p, radii, standard):
    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")


radius_sensitivity(C, w, p=8, radii=[300, 480, 600, 720, 900], standard=600)
optimised for  5 min: 65.9% within 5 min, 85.8% within the 10-min standard
optimised for  8 min: 86.9% within 8 min, 91.9% within the 10-min standard
optimised for 10 min: 93.8% within 10 min, 93.8% within the 10-min standard
optimised for 12 min: 97.4% within 12 min, 86.3% within the 10-min standard
optimised for 15 min: 99.6% within 15 min, 68.7% within the 10-min standard

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

free, free_sites = solve_mclp(C, w, radius=600, p=8)
fixed, fixed_sites = solve_mclp(C, w, radius=600, p=8, fixed=[49, 86])
print(f"spopt perc_cov: {free.perc_cov:.1f}% of demand points ({free.n_cli_uncov} uncovered)")
print(f"population covered: {covered_share(C, w, free_sites, 600):.1%} free, "
      f"{covered_share(C, w, fixed_sites, 600):.1%} with 49 and 86 committed; "
      f"{len(set(free_sites) & set(fixed_sites))} sites in common")
spopt perc_cov: 85.3% of demand points (329 uncovered)
population covered: 93.8% free, 91.2% with 49 and 86 committed; 4 sites in common

Explanation

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; the p-median model on the same matrix took 99.7 seconds.

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

Two panels comparing an unweighted maximal coverage plan, which covers more blocks but fewer people, with a population-weighted plan.
Weights decide what is being maximised; perc_cov reports the unweighted figure either way.

Edge cases or notes

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

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.8% 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.3%, the share of blocks; weighted by population it covered 93.8%.

Do I need to weight the demand points?

Yes, unless every point represents the same demand. An unweighted plan covered more blocks, 86.5%, but fewer residents, 92.8%, 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 68.7% of residents within 10.

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.8% to 91.2%.

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.5% of residents; the rest were beyond every candidate.