Lambda sweep method

The λ-sweep method

Since September 2026. The method behind the pharmacy results on the ServiceAccess branch: the model, how one λ picks one point and a sweep traces the frontier, how LP fractions become a real set of pharmacies, what the scenarios are, and what the model deliberately excludes. Code references are to lp_run.jl and lambda_sweep_simplex.jl — see Pharmacy service access.

1. The model — uncapacitated facility location with soft coverage

For clients i (populated 1 km² cells, weight pop_i = residents) and candidate cells j, over the OD pairs the road network delivers (t_ij ≤ 120 min, within the choice set):

min   Σ_i pop_i · [ Σ_j c(t_ij)·x_ij  +  BIG·(1 − Σ_j x_ij) ]   +   λ · Σ_j y_j

s.t.  Σ_j x_ij ≤ 1            a client may be left unserved
      x_ij ≤ y_j              only to open facilities
      0 ≤ y_j ≤ 1,  x_ij ≥ 0  (LP relaxation of y_j ∈ {0,1})

• x is assignment and y is facility openness — the p-median convention (ReVelle & Swain 1970).
• Soft coverage. A client the model cannot or will not route is priced at BIG — 120 min under LINEAR, 1.0 (the saturation value) under LOGISTIC — and at the same price in the baseline, so the star and the frontier are directly comparable. Writing the per-client cost as Σ_k x_k·(c(t_k) − BIG)·pop_k + BIG·Σ_i pop_i, the constant drops out and the x-coefficients are fixed across λ; only the y-coefficients change per sweep point. The LP may now choose to strand a remote client when a facility serving it would cost more than the BIG it saves, so the frontier extends below the full-coverage floor and every area brackets S1 and S2.
• λ = w · €100 000 is the price of a location. A linear facility cost a + b·q reduces to λ·#open while (nearly) all demand is served. The €100 000 is a placeholder, not an estimate: the frontier, the scenario deltas and the ranking of areas are invariant to it; only the € labels move with the true a.
• Why the relaxation. With y_j ∈ {0,1} this is the NP-hard uncapacitated facility-location problem in its strong disaggregated form (one x_ij ≤ y_j per OD pair), whose LP relaxation is known to be nearly integral; the sweeps confirm it (frac_y stays small). The LP optimum is a certified lower bound; the rounding below gives a feasible upper bound.
• Relation to the p-median. Imposing Σ_j y_j = p instead of pricing it gives the p-median problem; S1 at the baseline count is a p-median instance in soft-coverage form. The λ-sweep is its Lagrangian relaxation (Cornuéjols, Fisher & Nemhauser 1977): it reaches only the p's on the lower convex envelope of the value function, so p-values in non-convex gaps need a fixed-count solve.

Travel-cost functions

	c(t), t in minutes	BIG	Notes
LINEAR	t	120	person-minutes; unbounded, stranding expensive
LOGISTIC	1 / (1 + e^{−(t−25)/10})	1.0	adapted logit, midpoint 25 / scale 10: c(0) = 0.08, c(15) = 0.27, c(45) = 0.88, c(60) = 0.97. Lower below ~10 min (indifferent to minor relocations), saturating by ~45 min (caps remote weight)

Both are swept over the same 1-2-5/decade λ-grid. The logistic cost spans 0.08–1 where the linear spans 0–120, so its useful λ-range is compressed (many λ values map to one solution). A rescaling so that c(0) and c(60) match the linear curve is an open group decision.

2. The sweep — one LP per λ, warm-started

build_lp_warmstart builds the model once (HiGHS, presolve on). solve_at_w! rewrites only the y objective coefficients and re-solves with dual simplex from the previous optimal basis, which makes each step along the grid cheap. The grid is 1-2-5 per decade, w ∈ [1e-4, 5], common to every area, extended upward (×2.5, plateau-detected) until S1 brackets; SWEEP_MULTS gives a denser grid. LP_TIME_LIMIT (1 h) turns a pathological point into a skipped point. SOLVER=ipm forces interior point + crossover — basis-distance-immune, useful on the very largest LINEAR regions at high λ, but it fails on soft-coverage LOGISTIC (the coefficient range is too wide).

Editor's note (15 Sep 2026, against commit e25416b). The upward extension (×2.5, plateau-detected) was removed in commit f93aea1 (9 Jul 2026) because the top of the common grid was intractable; the grid now runs from 1e-4 up to SWEEP_WMAX (default 5.0 LINEAR, 0.5 LOGISTIC) and the driver prints a hint to raise SWEEP_WMAX if S1 is still not bracketed. See The lambda sweep and the Pareto frontier.

How one λ picks one point. We minimise travel + λ·#open; every solution with the same objective lies on a line of slope −λ. Slide the line down until it last touches the feasible set — that tangency is the λ-optimum. A steeper slope buys fewer, busier pharmacies. Many λ give many tangencies, and their lower-left envelope is the Pareto frontier. The frontier itself needs no λ; λ only says which point on it a decision-maker would pick.

The sweep read along λ. The frontier charts put λ on a secondary axis; the chart below puts it on the x-axis, with the two travel-cost bounds on the left and the facility count on the right, for the 41-area aggregate. The count collapses over two decades of λ while the bounds stay together; they separate only near S2. See Pharmacy results September 2026 for the reading and the Netherlands version; doc/lambda_axis_chart.py <REGION> draws any area.

Fractional sum_y is informative, not noise. Three clients, three candidates, unit weights, λ = 2: any one facility costs 1+1+3+2 = 7, any two cost 1+1+1+4 = 7, but y = (½,½,½) with each client split over its two near pharmacies costs 3 + 3 = 6 < 7. The LP strictly beats the integer optimum; geographies that resist a clean integer covering show the gap.

3. From fractions to pharmacies — the multistart rounding

An LP-guided matheuristic (multistart_round), seconds per point — the LP solve dominates:

1. Fix the count p = round(Σ_j y*_j), the integer count nearest the LP's fractional total, so every rounded point stays on the same axis as the relaxation.
2. Seed 10 candidate sets: top-p by y*; the lazy-greedy set (greedy_round, CELF — marginal travel gains are submodular, so after each pick only candidates sharing a now-better-served client are re-evaluated); and 8 y*-weighted random draws (Efraimidis–Spirakis A-Res). Must-open facilities (y* ≈ 1) always stay. Top-p alone starves regions the LP covered with many low-y substitutes; the greedy seed fixes that.
3. Polish by swaps (swap_round!): best-improving open↔closed swaps over the fractional pool — the classic p-median vertex-substitution search in its fast-interchange form (Resende & Werneck 2007). The swap profit is the exact travel change (gain − loss + interaction term), so polishing is monotone; ≤ 12 rounds per seed.
4. Score coverage-honestly, keep the best (travel_of): every client at its nearest open pharmacy at c(t); a client with no open pharmacy in its OD is priced at BIG, the same price as in the LP and the baseline, so abandoning remote clusters is never free. The unserved count is reported next to every figure.

Because the rounded set is a feasible point of the same soft-coverage problem, LP relaxation ≤ integer optimum ≤ multistart holds rigorously — a certified gap on total cost, not on travel and facility cost separately. Deterministic (fixed RNG seed: MS_RESTARTS = 10, MS_ROUNDS = 12, MS_SEED). Measured gaps: +0–18 % under LINEAR, +0–4.7 % under LOGISTIC.

4. Scenarios

	Definition	How it is read off
Baseline ★	Today's pharmacies; each client to its nearest by road time; unreachable clients priced at BIG (coverage-consistent).	baseline_metrics
S1 — same count, less travel	Travel minimised at today's number of locations.	Chord-interpolated on the multistart frontier at the baseline count.
S2 — same travel, fewer locations	Count minimised at today's travel cost.	Interpolated where the frontier crosses the baseline travel level.
S3 — the whole frontier	Every tangency between and beyond S1 and S2.	The sweep itself.

Gains are stated in native units first — locations, and person-minutes (LINEAR) or dimensionless logistic cost with the mean minutes beside it — which need no λ; the λ at which the sweep reaches each point prices them second. The exported S1/S2 location files snap to the nearest swept grid point; in five areas (FRC, FRH, FRJ, FRK, SE2) both scenarios snap to the same point. The interpolated tables are unaffected; an exact fixed-count solve is the planned refinement.

Editor's note (15 Sep 2026, against commit e25416b). Since commits 29fcf7f and e25416b (10 Sep 2026, issue #52) the driver pins S1 and S2 by bisection inside the sweep the moment their bracket closes (REFINE_TOL 0.2 % or one facility, REFINE_ITERS 8), so the exported rows sit at the target rather than on the coarse grid; the five-area coincidence above describes the pre-#52 snap. The three S1/S2 definitions in circulation (snapped, chord-interpolated, pinned) and which artefact uses which are tabulated on Scenarios S1 S2 S3.

5. Scope rules — what the model deliberately excludes

Under soft coverage every resident the model cannot route is priced at BIG. That is right for someone genuinely remote and wrong when the network or the supply data is simply missing, because it turns a data gap into apparent policy headroom. Two rules separate the cases.

1. Road network per landbody. The largest strongly-connected network is kept for every separate landbody, not just the single largest in the study area. Before, one component survived per area, so island networks — and the ferry links inside them — were discarded; Sicilia and Sardegna are separate landbodies, so at most one could ever be kept. No source data changed. ITG: mean travel 22.3 → 4.9 min, unreachable residents 983 169 → 591.
2. Regions the model cannot serve at all (issue #49). If more than 50 % of a NUTS region's inhabitant locations are absent from the OD matrix, the whole region — population and candidate locations — is excluded. Level NUTS3 where populated, else NUTS2, else NUTS1; judged on the existing network so baseline and sweep agree on who is in scope. Excluded today: PT200 Azores (919 cells, 212 855 residents) and PT300 Madeira (419 cells, 244 867 residents) — population but no pharmacy anywhere in the source data. Portugal's figures now cover the mainland only (baseline travel −62 %). No other area comes near the threshold.

After both rules, residents still unreachable within 120 min are 0.012 % of the total; 34 of 41 areas have none, the largest remainder is SE2 (165 cells, 14 210 residents, 0.32 %).

6. The choice set, and where it binds

Each client's OD row holds every candidate out to the road time of its 5th-nearest existing pharmacy (max_nr_facilities_per_client = 5, or 120 min). That radius holds ~70–90 candidates on average — it is not "5 choices", it is a radius fixed by today's five, and when the optimiser removes facilities the radius does not grow. Measured at S2: SE2 has 1 730 cells / 449 261 residents (10.2 %) with exactly one open facility left inside their radius, so the optimiser must keep it open even where a 6th-nearest pharmacy would serve them; ITF has 9.0 % already stranded and a further 23.9 % one closure away. The frontier is under-estimated wherever the radius, not the geography, decides who can be served. Widening the candidate-side radius (keeping the existing side at 5 so the baseline stays comparable) is the first remaining roadmap item, and a precondition for RSSV candidate reduction (which samples candidate subsets and needs every client to keep a reachable one).

7. Aggregation and metrics

• Aggregated frontier over the 41 disjoint areas (13 countries + FR/IT/SE/PL NUTS-1; country-level Poland dropped in favour of its 7 NUTS-1): at a common λ the sum of the regional optima is the combined optimum (separability). Summed at the union of swept w-values inside the range every area covers (68 points); an area without that exact λ is log-interpolated between its adjacent sweep points.
• Improvement rectangle per area (doc/frontier_metrics.py): baseline ↔ S1 (vertical) ↔ S2 (horizontal); its area, raw and relative to baseline count × travel; the crossing of the baseline→corner diagonal with the frontier; and λ at that crossing — the price per location at which the balanced improvement is the optimum. Areas are ranked on these, coloured by pharmacy-regulation typology (doc/policy_typology.csv).

8. Not pursued — the catchment-cap ladder

cap_scenario.jl (drivers s1_cap.jl, s2_cap.jl) implements the cost-function-free rungs of the team's A→D ladder: S1-A/B/C — a fixed count under a max-catchment cap, multiple or one pharmacy per cell, urban cells held fixed; S2-A/B — a minimum count under min + max caps, globally or outside urban centres. Run on the Netherlands (S1-A: +13.3 % travel vs today at a 35 000 cap, multiple per cell; S1-B: −6.5 % at one per cell) and dropped on 4 September 2026: observed catchments vary so widely that realistic min/max bounds cannot be set. The λ-sweep is the general method.

References

ReVelle & Swain 1970 (p-median) · Cornuéjols, Fisher & Nemhauser 1977 (Lagrangian relaxation) · Leskovec et al. 2007 (CELF) · Efraimidis & Spirakis 2006 (A-Res weighted sampling) · Resende & Werneck 2007 (fast interchange) · Albuquerque, Figueiredo & Genre-Grandpierre 2026, SSRN 7133060 (RSSV candidate reduction, territorial coverage constraints — the planned next steps).