The facility location model explained

The facility location model explained

In one sentence: the optimisation problem the Julia side solves at each price λ — which candidate cells to open and where each resident goes — in words before symbols, so you can read a sweep log or a frontier chart and know what it can and cannot claim.

Where this sits: GeoDMS has produced the OD matrix (GeoDMS OD matrix and choice set); this page is the model lp_run.jl builds from it, which the sweep (The lambda sweep and the Pareto frontier), the rounding (From LP fractions to real pharmacies) and the scenarios (Scenarios S1 S2 S3) assume. Formal statement: Lambda sweep method §1; OR depth: Background theory.

The question the model answers

Suppose a pharmacy location cost a fixed amount λ (lambda — in the code w times a placeholder €100 000; Glossary, The lambda sweep and the Pareto frontier), and every minute a resident travels to reach one also had a price. Which cells would you open, and where would everybody go? The ingredients:

• Clients i: every populated 1 km² cell, weighted by its residents pop_i (Glossary → client).
• Candidates j: every cell with at least 50 inhabitants (MinPopSizeForPharmacyLocation) plus every cell holding a pharmacy today; why 50 and what it costs: GeoDMS OD matrix and choice set.
• Pairs (i, j) with road time t_ij in minutes: only the pairs the OD matrix contains — the client's choice set, the candidates no farther by road than its 5th-nearest existing pharmacy, capped at 120 min (Glossary → choice set; where the radius binds: Scope rules and data gaps). A pair not in the OD does not exist for the model.

Two kinds of decision: y (is it open?) and x (who goes where?)

• y_j — is candidate j open? In reality 1 or 0.
• x_ij — what share of client i is served by candidate j? In reality also 1 or 0 (one pharmacy per client).

These are the letters of the p-median literature (below), adopted on 24 July 2026 (commit 4cee238). Before that date the code, logs and deck used the opposite letters — see Pitfalls.

The cost, term by term

The model minimises one total with three parts:

1. Travel. A resident served by j costs c(t_ij) — minutes under LINEAR, a 0-to-1 S-curve under LOGISTIC (Travel cost functions) — times pop_i · x_ij.
2. Stranding. A resident served by nobody costs a flat BIG; the unserved share of client i is 1 − Σ_j x_ij.
3. Locations. Every open location costs λ.

In symbols (Lambda sweep method §1):

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

s.t.  Σ_j x_ij ≤ 1            each client is served at most once; the rest is unserved
      x_ij ≤ y_j              you can only use an open share of a candidate
      0 ≤ x_ij ≤ 1,  0 ≤ y_j ≤ 1

In words: travel of the served, plus the stranding price of the unserved, plus the bill for the open locations; nobody is served twice; nobody is sent to a closed door. A large λ buys few, busy pharmacies and a small λ many; sweeping λ traces the frontier, the least travel achievable for each number of open locations (Glossary). λ is a pure per-location price: the school-era per-client cost term is constant once everybody is served, so facility_cost() in settings.jl is never called and the € label is a placeholder (Alternatives not pursued).

Why population weights? Without pop_i a hamlet of 12 and a suburb of 8 000 would pull equally — cell-minutes, not person-minutes. The column is total_pop (env CLIENT_WEIGHT); clients in an excluded NUTS region keep weight 0 (Scope rules and data gaps).

Why a client may be left unserved (soft coverage), and what BIG does

The coverage constraint is Σ_j x_ij ≤ 1, not = 1. With = 1 (hard coverage, still in the legacy run_scenario) every client must be assigned, so every area has a full-coverage floor — the fewest locations that reach every client — below which no λ can push the optimiser. In several areas that floor sat above today's pharmacy count (Netherlands: 1 899 against 1 615 cells, commit acd59a6), so the frontier could not reach S1, and in four areas not S2 either.

Soft coverage (agreed 10 July 2026, commit 4ec47a9) removes the floor: a client may stay unserved at BIG per resident, so the optimiser can abandon a remote hamlet when the location that would serve it costs more than the stranding it saves. The frontier then runs past today's count and every area brackets S1 and S2 (Glossary; Scenarios S1 S2 S3).

BIG is the price of one unserved resident, big_cost() in settings.jl (around line 97; env BIG_TRAVELTIME_MIN): c(120) = 120 person-minutes under LINEAR, 1.0 under LOGISTIC. Every OD pair has c(t_ij) ≤ BIG, so serving is never worse than stranding: the model strands only when the location needed is not worth it. Why those two values, and how BIG went from a per-region maximum to a fixed c(120): Travel cost functions.

The baseline ★ (baseline_metrics in lambda_sweep_simplex.jl, around line 167) charges the same BIG to residents the Existing-side OD cannot reach — coverage-consistent since 206923c; before that they were silently dropped, so the ★ could sit below the LP lower bound (Scenarios S1 S2 S3).

The "constant drops out" trick

Multiply out the per-client bracket:

pop_i · [ Σ_j c(t_ij)·x_ij + BIG·(1 − Σ_j x_ij) ]
   =  Σ_j (c(t_ij) − BIG)·pop_i·x_ij   +   BIG·pop_i

The last term is what you would pay if everybody were stranded — the same for every answer, so it cannot change which answer is smallest. build_lp_warmstart leaves it out and gives every x_ij the coefficient (c(t_ij) − BIG)·pop_i (around line 420): each unit of assignment earns back the difference between stranding and serving. No "unserved" variable is needed — the slack under the ≤ 1 is the unserved share — and no x-coefficient depends on λ: between two sweep points only the coefficients of the y_j change (set_objective_coefficient in solve_at_w!, around lines 466–469), which is what lets HiGHS warm-start each grid point from the previous answer (The lambda sweep and the Pareto frontier).

The trap: the objective HiGHS prints in the log is Σ x·(c − BIG)·pop + λ·Σy, a large negative number under LINEAR, not the total cost. solve_at_w! re-adds the constant when it computes travel_relax (around lines 502–513). Read travel_relax and sum_y from the sweep row, never the HiGHS objective (Log lines and output files).

Why solve the relaxation, and what "nearly integral" means

In the real problem y_j ∈ {0, 1}. That integer problem is NP-hard, and at up to fifteen million OD rows per area and 15–40 solves per sweep it is out of reach. build_lp_warmstart solves the LP relaxation instead (Glossary): 0 ≤ y_j ≤ 1 (around line 406), so a candidate may be half open and a client half served by two pharmacies — an ordinary linear programme.

The formulation is "strong"

x_ij ≤ y_j is written once per OD pair (around line 426), not once per candidate as the weak form Σ_i x_ij ≤ |I|·y_j would. Under the weak form a candidate serving one client in a hundred needs only y_j = 0.01; under the strong form a candidate that fully serves anybody must be fully open. The cost is size: one row per OD pair.

The reward is that the relaxation is nearly integral: most y_j come out at exactly 0 or 1 anyway (Background theory §3 says why). frac_y in every log row counts the candidates strictly between 0 and 1: few against the tens of thousands of candidates, but ~10 % of the open count near S1/S2 (median 9 % over the 148 S1/S2 points in doc/deck_data.json, September 2026) and more than the open count at high λ — which is why the rounding step matters.

The relaxed optimum is a certified lower bound

Every integer answer is also allowed by the relaxation, so the relaxation can only do as well or better; the multistart set the rounding step builds from the fractions (From LP fractions to real pharmacies) is a real network, so it can only do as well or worse. In words: the LP's travel plus λ times its fractional count (sum_y) never exceeds the best real answer, which never beats the rounded set's travel (travel_multi) plus λ times its count (n_open):

travel_relax + λ·sum_y   ≤   integer optimum at λ   ≤   travel_multi + λ·n_open

The sandwich holds on total cost — travel plus λ times count — not on travel alone; the deck's travel-only "multi vs relax" column can therefore be off by up to λ/2, even negative. Worked case and measured gaps: From LP fractions to real pharmacies.

A toy where the fraction beats every whole-number answer

Fractional y can be strictly better than any real set of pharmacies — the example Lambda sweep method §2 quotes, with the matrix it omits: three clients of one resident each, three candidates on a ring, each client between two of them.

c(t_ij)	cand 1	cand 2	cand 3
A	1	1	3
B	3	1	1
C	1	3	1

Take λ = 2 and BIG = 120, so stranding never pays. Total = travel + λ · open. One candidate: two clients at 1, one at 3: travel 5 + locations 2 = 7. Two candidates: everybody at 1: 3 + 4 = 7. All three: 3 + 6 = 9.

The best real answer costs 7. The LP sets y = (½, ½, ½) and splits every client half-and-half over its two near candidates — feasible because x_ij = ½ ≤ y_j and Σ_j x_ij = 1. Travel 3; locations 2 × 1.5 = 3; total 6 < 7. HiGHS confirms it (sum_y = 1.5, frac_y = 3). Rounding fixes p = round(1.5) = 2 (Julia rounds half to even) and returns a pair at 7: a certified gap of 1 on a bound of 6.

With "today's count" = 2 the same toy shows why S1 can need more than a sweep: the LP reports 1.5 for λ < 4 and one site above, never 2 — the integrality gap, not a concave dent — so only a fixed-count solve (Σ_j y_j = 2) would land on 2. Rounding still returns a pair, so "today's count" is read from the rounded set.

Relation to UFL and p-median, briefly

• UFL (uncapacitated facility location, Glossary): fixed cost per open site plus assignment cost; the count is an output. This model is a UFL with the same λ everywhere and an outside option (stranding at BIG).
• p-median: minimise weighted travel with exactly p sites; the count is an input. S1 is a p-median instance in soft-coverage form.

Pricing locations at λ instead of fixing p is a Lagrangian relaxation of the count constraint (Glossary): it can only land on counts on the lower convex envelope of the integer frontier and skips any concave dent — smallest example and the mechanism of issue #54 on The lambda sweep and the Pareto frontier, proof sketch in Background theory §5. Hence S1/S2 are bracketed between two sweep points and bisected on λ, not solved directly; at scale dents are shallow (Scenarios S1 S2 S3).

In the code

Open lp_run.jl at build_lp_warmstart (around line 397, permalink) and read the comment above the objective; then solve_at_w! (around line 460, permalink). The data come from load_from in lambda_sweep_simplex.jl (around line 71), not settings.jl → load_country, which lacks the NUTS exclusion and the candidate subsample. Function by function: Code walkthrough lp_run; every env var: Code walkthrough settings. Always set TRAVEL_FUNC: its default QUADRATIC is a likely school-era leftover (inference) nobody reports.

Wiki symbol	In the code
client i, candidate j	locations[i] = client's OD row indices; j ∈ facilities, sorted facility_rel ids
pair (i, j), t_ij	OD row k: clients_col[k], facilities_col[k], t_ij_col[k] (od[:t_ij] ./ 60)
pop_i	wpop[k] per OD row; client_pop[i] per client
x_ij, y_j	x[k], one variable per OD row; y[j], keyed by facility id
c(t), BIG, λ	c(t), big_cost() in settings.jl; λ = w * FACILITY_MIN_COSTS
Σ_j y*_j, fractional count, LP travel	sum_y; n_frac (log: frac_y); travel_relax
Σ_j x_ij ≤ 1, x_ij ≤ y_j	sum(x[k] for k in rows) <= 1 per client; x[k] <= y[facilities_col[k]] per OD row

Pitfalls and open questions

• The HiGHS objective in the log is not the cost. Read travel_relax and sum_y.
• The certified gap is on total cost, not on travel alone.
• x and y were swapped before 24 July 2026 (commit 4cee238): sum_x meant the number of open locations in older logs and issue #47; a human reading an old log must swap them.
• mean_t is asymmetric — sweep rows count stranded residents at 0 min, the baseline at 120 min; compare cost_c, not mean_t (Log lines and output files).
• The choice-set radius caps the model. A client whose OD holds one open candidate must keep it or pay BIG; where the radius rather than the geography decides, the frontier is under-estimated (Scope rules and data gaps).
• Open: an exact soft-coverage p-median MIP to pin S1/S2 (doc/todo.md B5); the LOGISTIC rescaling (Travel cost functions); the mean_t asymmetry; which HiGHS termination statuses are accepted (Code walkthrough lp_run).

As of September 2026 (commit e25416b).

See also

• Lambda sweep method · Background theory · Travel cost functions · The lambda sweep and the Pareto frontier · From LP fractions to real pharmacies · Scenarios S1 S2 S3 · Scope rules and data gaps · Glossary