← Past problems · 2013 set

2013 · Problem A — Emergency Medical Response

Maximum Coverage Location Integer programming Facility location Sensitivity analysis

Read the official problem page →

The prompt, restated

A small regional EMS agency operates across six zones with known populations and a known 6×6 matrix of zone-to-zone average ambulance travel times (in minutes). The agency's service standard is to reach an emergency within 8 minutes of dispatch. Teams are asked: (1) Choose the placement of three ambulances — one per zone, with at most one ambulance per zone — that maximizes the population reachable within the 8-minute standard. (2) Re-run the analysis with only two ambulances, then with one, and report which zones lose coverage as the fleet shrinks. (3) Discuss catastrophic-event preparedness: what happens when one ambulance is already committed to a multi-casualty incident, or when a hospital diversion forces longer turn-around times? (4) Write a one-page memo to the Emergency Service Coordinator summarizing the recommended placement, the coverage gaps, and any operational caveats.

Note the prompt is deliberately small-N — six zones is well inside brute force. The points come from the defense: how you weight equity vs. raw coverage, what you do when ties appear, and how cleanly you communicate the result.

Key modeling idea

This is the classical Maximum Coverage Location Problem (MCLP) (Church & ReVelle, 1974). Let $a_{ij} = 1$ if zone $i$ can reach zone $j$ within 8 minutes (i.e., $t_{ij} \le 8$). Let $x_j \in \{0,1\}$ indicate whether an ambulance is stationed in zone $j$, and $y_i \in \{0,1\}$ indicate whether demand zone $i$ is covered. The integer program is $$\max \sum_i p_i\, y_i \quad \text{s.t.}\quad y_i \le \sum_j a_{ij} x_j,\;\; \sum_j x_j = k,\;\; x_j,y_i\in\{0,1\}$$ with $p_i$ the population of zone $i$ and $k\in\{1,2,3\}$ the fleet size. With only six zones, $\binom{6}{3}=20$ placements — enumerate all of them.

Suggested approach

  • Step 1 — Build the coverage matrix. From the official travel-time table, compute $a_{ij}=\mathbb{1}[t_{ij}\le 8]$. Sanity-check that every zone covers itself ($t_{ii}=0$) and that the matrix is roughly symmetric.
  • Step 2 — Enumerate placements. For each $k\in\{1,2,3\}$, iterate over all $\binom{6}{k}$ subsets, compute total covered population, and rank (see technique 11 — optimization).
  • Step 3 — Break ties with equity. When two placements cover the same total population, prefer the one whose minimum zone coverage probability is higher (max-min fairness). Document the tiebreaker.
  • Step 4 — Degrade the fleet. Report the marginal value of each successive ambulance: coverage jumps from one to two units, then two to three. Diminishing returns is the headline finding.
  • Step 5 — Catastrophe scenario. Model one ambulance as "busy with probability $\rho$" and recompute expected covered population. Tie $\rho$ to a Poisson call-arrival rate using the Erlang-B / loss formula (technique 9 — Markov chains & queueing).

Data sources to consider

SourceWhat you get
NFPA 1710 (Standard for Career Fire Departments)The 8-minute / 90th-percentile EMS response benchmark used in U.S. cities
NEMSIS (National EMS Information System) public reportsReal call-volume and response-time distributions by region
Church & ReVelle (1974), "The Maximal Covering Location Problem"The canonical MCLP formulation and proof of NP-hardness in general
U.S. Census tract population dataPlug-in for $p_i$ if you want to extend beyond the toy six-zone case
OSRM / OpenStreetMap travel-time matricesRealistic $t_{ij}$ if you want to validate the given table against road geometry

Common pitfalls

  • Greedy without checking optimality. Greedy MCLP is a $(1-1/e)\approx 63\%$ approximation in general — on six zones, it often gets the wrong answer. Just enumerate.
  • Confusing coverage with response. A zone is "covered" if some ambulance can reach it in 8 minutes — not that the response will actually happen in 8 minutes once queueing and dispatch overhead are added. State this gap explicitly.
  • Treating the travel-time matrix as deterministic. Real travel times have ~30% coefficient of variation [illustrative]; convert the 8-minute threshold into a probability $P(T_{ij}\le 8)$ and report expected covered population.
  • Forgetting equity. A placement that covers 85% of population but leaves one entire low-density zone with zero coverage is politically untenable. Always report the minimum per-zone coverage alongside the total.
  • No catastrophe analysis. The prompt asks for it directly. Don't just hand-wave — model busy probability with a one-server M/M/1/0 loss formula or a Monte Carlo (technique 10).

Python sketch

Brute-force MCLP for $k=1,2,3$ with a busy-ambulance Monte Carlo extension.

import numpy as np
from itertools import combinations

# --- Inputs (replace with the official table) ---
# travel times t_ij in minutes between the six zones [illustrative]
T = np.array([
    [ 0,  4,  6, 11,  9, 14],
    [ 4,  0,  5, 10,  8, 12],
    [ 6,  5,  0,  7,  6,  9],
    [11, 10,  7,  0,  5,  6],
    [ 9,  8,  6,  5,  0,  7],
    [14, 12,  9,  6,  7,  0],
])
pop = np.array([21000, 35000, 15000, 30000, 20000, 18000])  # [illustrative]
A   = (T <= 8).astype(int)              # coverage matrix

def covered(placement):
    reachable = A[list(placement), :].max(axis=0)  # zone i covered if ANY chosen j covers it
    return int((reachable * pop).sum()), reachable

# --- Step 2: enumerate placements for k = 1, 2, 3 ---
for k in (1, 2, 3):
    best = max(combinations(range(6), k), key=lambda s: covered(s)[0])
    total, mask = covered(best)
    print(f"k={k}: stations at {best}, covered = {total:,} of {pop.sum():,}")
    print(f"       per-zone covered? {mask.tolist()}")

# --- Step 5: catastrophe — one ambulance busy with prob rho ---
rng  = np.random.default_rng(0)
rho  = 0.35                              # call-volume-driven busy fraction
trials = 20_000
placement = (1, 2, 4)                    # best k=3 plan from above [illustrative]
totals = []
for _ in range(trials):
    available = [j for j in placement if rng.random() > rho]
    totals.append(covered(tuple(available))[0] if available else 0)
print(f"E[covered | rho={rho}] = {np.mean(totals):,.0f}")

Sensitivity & validation checklist

  • Sweep the coverage threshold from 6 to 10 minutes — the optimal placement should be stable across at least ±1 minute.
  • Perturb each $t_{ij}$ by Gaussian noise ($\sigma = 1$ min) and recompute the optimal set 1,000 times; report the placement-stability frequency.
  • Vary populations $\pm 20\%$ per zone — if the optimal flips on small population changes, the model is fragile and you need a regret-minimization formulation.
  • Cross-check the enumeration against a tiny ILP solver (e.g., PuLP / OR-Tools) to confirm both methods agree.
  • Validate the catastrophe model: $\rho$ should be consistent with $\lambda \bar{s}$, where $\lambda$ is call arrival rate and $\bar{s}$ is mean service time per call.

Related pages