Fixing a Location-Allocation Model That Never Solves
Problem statement
A p-median or coverage model that worked on a test area runs for an hour on the real one, returns "Not Solved", reports "Infeasible", or fails before the solver starts. The formulation is usually correct; the model is too large, the travel-time matrix contains something the solver cannot use, or the question has no answer. Each has a different fix, and a longer time limit fixes none of them.
Measured with spopt and PuLP's CBC solver, choosing 8 sites for the 2,241 populated census blocks of Chittenden County, Vermont:
- With 200 candidate sites the p-median model had 448,400 variables and 450,442 constraints. Given a 300-second limit, CBC returned after 508 seconds with no solution, and spopt raised
RuntimeError: Model is not solved: Not Solved. - Accepting a 5% optimality gap did not help: the same error after 505 seconds.
- Keeping only block–site pairs within 20 minutes — 55.6% of the variables — solved the restricted model to optimality in 227.8 s. A 15-minute cutoff made the same model infeasible, and on 106 candidates the 20-minute model returned sites 6.7% worse than the unrestricted optimum.
- One unreachable block in the matrix failed in 0.3 seconds with
PulpError: Cannot multiply variables with NaN/inf values.
Quick answer
Work out which failure you have before changing anything:
symptom cause fix
RuntimeError ... Not Solved (after the limit) too many variables for the solver aggregate demand, sparsify, prune
RuntimeError ... Infeasible no set of sites satisfies constraints check the radius, the cutoff, the candidates
PulpError: Cannot multiply ... NaN/inf unreachable or missing travel times fix the matrix before building
"Optimal" but slow or poor time limit hit with an incumbent read sol_status; shrink the model
p-centre far slower than p-median minimax objective solve p-median first, or limit candidates
Then check the model size — demand points × candidates — against what your solver handles: here 237,652 variables took 100 seconds and 448,400 did not finish.
Step-by-step solution
1. Read the status, not just the traceback
spopt checks PuLP's status after every solve and raises for anything but optimal. Not Solved means the solver stopped before finding any feasible solution — usually the time limit. Infeasible means it proved there is none. A PulpError before solving means the model could not be built. When the time limit stops CBC after it has a solution, PuLP reports Optimal anyway; model.problem.sol_status of 2 instead of 1 shows the solution was not proved optimal.
2. Check the matrix before building
Every travel time must be a finite number. A single block that cannot reach any candidate — a inf from the shortest-path search — made PuLP refuse to build the objective in 0.29 s, and the same with NaN in 0.31 s. Find such rows, fix the network or snapping, or drop the demand point knowingly (Example 1).
3. Check that a feasible answer exists
Coverage and set-covering models can be impossible. With a 10-minute standard, set covering on the county's 106 school sites raised Model is not solved: Infeasible, because the worst-placed block was 15.5 minutes from every school. The smallest feasible radius is the largest value, over demand points, of each point's nearest candidate time — compute it before choosing a standard.
4. Estimate the size before you build
A p-median model has one assignment variable per demand–candidate pair plus one per candidate, and as many "assign only to open sites" constraints. On the same 2,241 blocks the solve time grew faster than the model:
candidates variables solve status
25 56,050 7.6 s Optimal
50 112,100 19.8 s Optimal
100 224,200 62.3 s Optimal
106 237,652 101.6 s Optimal
200 448,400 >300 s Not Solved
Building the 200-candidate model took 5.5 s and 924 MB; the run's peak memory reached 2.5 GB.
5. Do not expect the time limit to rescue a model
The 300-second limit returned after 508 seconds of wall time — writing and reading the model file and the solver's own start-up count too — and with nothing to return. A 5% relative gap (gapRel=0.05) changes when CBC stops after it has a solution; it does nothing if CBC never finds one. Shrink the model instead.
6. Aggregate the demand
Demand points are the larger dimension. Solving the 106-school model on 118 block groups instead of 2,241 blocks took 0.34 s instead of 101.6 s, and the chosen sites were 0.41% worse when evaluated on the blocks — see choosing facility locations with a p-median model.
7. Drop assignments nobody would make
A block will never be served by a site an hour away when eight sites are open. Creating assignment variables only for pairs within a cutoff, plus each block's nearest candidate, cut the 200-candidate model to 249,041 variables at 20 minutes, and it solved to optimality in 227.8 s with a mean trip of 5.37 minutes. At 15 minutes it had 165,294 variables and was infeasible: some blocks' only nearby candidates could not all be served by 8 open sites. That optimum is optimal for the restricted model only: a cutoff forbids any answer that leaves a block beyond it. On the 106 schools, where the unrestricted optimum left 29 blocks beyond 20 minutes, the 20-minute model chose worse sites (Example 2). Choose a cutoff comfortably above the travel times you expect in the answer, and re-solve with a larger one if the result changes.
8. Prune candidates only where it removes many
Keeping only candidates that were the nearest site for at least one block removed 6 of 200 here, and the model still did not solve in 395 s. Pruning helps when candidate lists are dense — every road intersection, every parcel — and many are dominated by neighbours.
9. Choose the objective with solve time in mind
On the first 50 candidates, p-median solved in 20.5 s and p-centre in 117.9 s; on the 106 schools, 99.7 s and 916.4 s. Minimising the worst trip gives the solver a far weaker bound to prune with. Solve p-median or maximal coverage first, then p-centre on a reduced candidate list if the worst trip matters — see location-allocation explained.
Code examples
Example 1 — preflight a cost matrix
import numpy as np
def preflight(C, radius=None, cutoff=None):
"""Report size, bad values and feasibility before building a location model (C: demand x candidates, seconds)."""
n, k = C.shape
bad = ~np.isfinite(C)
print(f"{n:,} demand x {k} candidates -> p-median variables {n * k + k:,}")
if bad.any():
rows = np.flatnonzero(bad.all(axis=1))
print(f" {int(bad.sum()):,} non-finite costs; {len(rows)} demand points reach no candidate: {rows[:10].tolist()}")
nearest = np.where(bad, np.inf, C).min(axis=1)
print(f" smallest radius that covers everyone: {nearest.max():.0f} s")
if radius is not None:
print(f" radius {radius} s: {int((nearest > radius).sum())} demand points cannot be covered")
if cutoff is not None:
kept = int((C <= cutoff).sum() + (np.where(bad, np.inf, C) > cutoff).all(axis=1).sum())
print(f" cutoff {cutoff} s keeps {kept:,} assignment variables ({kept / (n * k):.1%})")
preflight(C, radius=600, cutoff=1200)
broken = C.copy()
broken[0, :] = np.inf # one block that reaches no school
preflight(broken)2,241 demand x 106 candidates -> p-median variables 237,652
smallest radius that covers everyone: 929 s
radius 600 s: 24 demand points cannot be covered
cutoff 1200 s keeps 121,855 assignment variables (51.3%)
2,241 demand x 106 candidates -> p-median variables 237,652
106 non-finite costs; 1 demand points reach no candidate: [0]
smallest radius that covers everyone: inf s
C is the travel-time matrix in seconds from the 2,241 blocks to the county's 106 schools, as built in measuring network distance for many pairs, and w holds each block's population.
Example 2 — a p-median model with assignments only within a cutoff
import time
import pulp
def sparse_p_median(C, w, p, cutoff, time_limit=600):
"""p-median with assignment variables only where C <= cutoff, plus each demand point's nearest candidate."""
n, k = C.shape
prob = pulp.LpProblem("sparse_p_median", pulp.LpMinimize)
open_ = [pulp.LpVariable(f"y{j}", cat="Binary") for j in range(k)]
assign = {}
for i in range(n):
allowed = np.flatnonzero(C[i] <= cutoff)
if allowed.size == 0:
allowed = np.array([int(np.argmin(C[i]))])
for j in allowed:
assign[i, int(j)] = pulp.LpVariable(f"x{i}_{j}", lowBound=0, upBound=1)
prob += pulp.lpSum(w[i] * C[i, j] * x for (i, j), x in assign.items())
by_demand = {}
for (i, j), x in assign.items():
prob += x <= open_[j]
by_demand.setdefault(i, []).append(x)
for xs in by_demand.values():
prob += pulp.lpSum(xs) == 1
prob += pulp.lpSum(open_) == p
start = time.perf_counter()
prob.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=time_limit))
sites = [j for j, y in enumerate(open_) if y.value() and y.value() > 0.5]
print(f"cutoff {cutoff / 60:.0f} min: {len(assign):,} assignment variables, {time.perf_counter() - start:.1f} s, "
f"status {pulp.LpStatus[prob.status]} (sol_status {prob.sol_status}), sites {sites}")
return sites
sparse_sites = sparse_p_median(C, w, 8, 1200)cutoff 20 min: 121,855 assignment variables, 132.2 s, status Optimal (sol_status 1), sites [11, 24, 52, 62, 71, 72, 77, 99]
Assignment variables are continuous: with fixed open sites, sending each demand point to its nearest open site is always optimal, so they need not be declared binary.
On the 106 schools this model was not a shortcut. The full model's optimum, from Example 3, left 29 blocks more than 20 minutes from their nearest open school — the worst at 26.9 minutes — so the 20-minute model could not choose it. Its own optimum was a different set of sites with a 6.7% higher population-weighted travel time (mean trip 5.61 minutes against 5.26), and it took 132.2 s against 99.4 s for the full model. A cutoff removes options as well as variables; it pays off only when the model is too large to solve otherwise, and its answer should be checked against a larger cutoff.
Example 3 — solve with a status you can trust
from spopt.locate import PMedian
def solve_and_report(C, w, p, time_limit):
start = time.perf_counter()
model = PMedian.from_cost_matrix(C, w, p_facilities=p)
try:
model.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=time_limit))
except RuntimeError as error:
print(f"{time.perf_counter() - start:.1f} s: {error}")
return None
proved = model.problem.sol_status == pulp.LpSolutionOptimal
sites = [j for j, v in enumerate(model.fac_vars) if v.value() > 0.5]
print(f"{time.perf_counter() - start:.1f} s: {pulp.LpStatus[model.problem.status]}, "
f"{'proved optimal' if proved else 'best found within the limit'}; sites {sites}")
return sites
solve_and_report(C, w, 8, time_limit=20)
dense_sites = solve_and_report(C, w, 8, time_limit=900)82.2 s: Model is not solved: Not Solved. See ``pulp.constants.LpStatus`` for more information.
99.4 s: Optimal, proved optimal; sites [24, 49, 51, 63, 68, 70, 72, 77]
The 20-second limit counts CBC's processing time, not the wall-clock time around it, which is why the first call returned after 82.2 s — with nothing.
Explanation
Why p-median models blow up
The assignment formulation needs a variable and a constraint for every demand–candidate pair. Doubling candidates doubles both, and the linear relaxation the solver works from grows harder to tighten. Branch and bound has to rule out combinations of open sites; with 200 candidates choosing 8 there are 5.5 × 10¹³ of them, and a weak relaxation prunes few.
Why a time limit returned nothing
CBC can only return what it has found. On a large p-median model it may spend the whole limit on the root relaxation and preprocessing before its first feasible integer solution. The time limit then expires with no incumbent, and spopt, correctly, refuses to report sites. A relative gap stops a search that has a solution; it cannot create one.
Why a cutoff can make a model infeasible
Removing assignment options removes feasible solutions. If a demand point's only allowed candidates are all far from the few that can be opened, no choice of p sites serves every point. The fallback of always keeping each point's nearest candidate does not guarantee feasibility for small p, which is what happened at 15 minutes.
Why aggregation is so effective
It shrinks the larger dimension of the matrix by an order of magnitude or more, and it removes near-duplicate demand points that make the relaxation hard. The price is a small, measurable loss of quality, which you can check by evaluating the chosen sites on the original demand.
Edge cases or notes
- Temporary model files are written by PuLP for CBC; a 448,400-variable model produced files of hundreds of megabytes, and a full temporary directory makes the solve fail.
- Commercial and HiGHS solvers handle larger models; switching solver changes the limits, not the causes.
- Warm starts from a heuristic solution give CBC an incumbent, so a time limit can return something.
- Capacity constraints make models much harder and can make them infeasible when total capacity is short.
- Duplicate candidates at the same network node add variables without adding options; deduplicate by node.
- Fixed existing facilities reduce the search; add them as constraints rather than candidates to choose.
- Coverage models stay small because they need only which pairs are within the radius; maximal coverage solved in 0.19 s where p-median took 99.7 s.
Internal links
- How to choose facility locations with a p-median model in Python — the model, sizes and aggregation
- Location-allocation explained: p-median, coverage and what each optimises — objectives and solve times
- How to maximise coverage with a limited number of sites — the model that solves quickly
- How to measure network distance for many origin–destination pairs — building a matrix without gaps
- Fixing an OSMnx graph that is disconnected — where infinite travel times come from
- How to find the areas nobody can reach within a travel time — the demand that makes covering infeasible
- Census geographies explained: blocks, tracts, output areas and why they nest — aggregation levels
- How to measure distance to the nearest facility for every home — nearest-candidate times for preflight
FAQ
Why does my PuLP location model say "Not Solved"?
The solver stopped, usually at the time limit, before finding any feasible solution. With 448,400 variables CBC found none in 300 seconds; shrink the model rather than extending the limit.
What does "Model is not solved: Infeasible" mean in spopt?
No choice of sites satisfies the constraints. A 10-minute set-covering model was infeasible because some blocks were 15.5 minutes from every candidate.
Why does PuLP say it cannot multiply variables with NaN or inf values?
The cost matrix contains infinite or missing travel times, often from demand points that cannot reach any candidate. Fix or remove them before building the model.
Does a larger optimality gap make a model solve?
Only if the solver has already found a solution. A 5% gap on the 200-candidate p-median still ended with no solution after 505 seconds.
How can I make a p-median model smaller?
Aggregate demand, remove assignment variables beyond a sensible travel-time cutoff, and deduplicate candidates. A 20-minute cutoff made the 200-candidate model solvable in 227.8 s, but a cutoff can exclude the true optimum, so check the answer with a larger one.
Why is p-centre so much slower than p-median?
Minimising the single worst trip gives the solver a weak bound to prune with. On 50 candidates p-centre took 117.9 s against 20.5 s for p-median.