Fixing a Location-Allocation Model That Never Solves

On this page

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,000 variables. Given a 300-second limit, CBC returned after 313 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 312 seconds.
Keeping only block–site pairs within 20 minutes — 50.3% of the variables — let CBC solve the model inside a 600-second limit: proved optimal in 360 seconds on an idle machine, while an earlier run on a busy one stopped at the limit holding a solution it had not proved. A 15-minute cutoff made the same model infeasible, and on 106 candidates the 20-minute model returned sites 6.8% worse than the unrestricted optimum.
One unreachable block in the matrix failed in half a second 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,440 variables took 160 seconds and 448,000 did not finish.

Three different errors, three different causes; only one of them is about speed.

Read the status, check the matrix, estimate the size

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.48 s, and the same with NaN in 0.49 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,000     5.1 s     Optimal
    50         112,000    18.5 s     Optimal
   100         224,000    64.5 s     Optimal
   106         237,440    88.6 s     Optimal
   200         448,000   >300 s      Not Solved

Building the 200-candidate model took 5.4 s before the solver started. The candidates beyond the 106 schools are 94 populated blocks drawn with a fixed seed, and the times, like all those in Example 4, were measured at a load average of 5 to 9 on a 24-thread machine.

5. Do not expect the time limit to rescue a model

The wall time exceeds the limit because writing the model and starting the solver count too.

The 300-second limit returned after 483 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.43 s instead of 89 s, and the chosen sites were 0.81% 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 225,299 variables at 20 minutes, and CBC found a solution within its 600-second limit, reported with sol_status 2 — feasible, not proved optimal. At 15 minutes it had 146,363 variables and was infeasible: some blocks' only nearby candidates could not all be served by 8 open sites. Even a proved optimum would be 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 15 of 200 here, which still leaves a model of about 414,000 variables — too many. 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 18.5 s and p-centre in 52 s on an idle machine; on a busy one an earlier run took 24 s and 988 s, so the penalty for the minimax objective ranges from three-fold to forty-fold with the conditions. On the 106 schools p-centre hit a 900-second limit without proving anything, while p-median needed 89 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.

Only reducing the assignment variables turned "Not Solved" into an answer.

Preflight checks and a cutoff-limited model

The data used on this page

Chittenden County, Vermont, from the 2020 Census: TIGER/Line 2020 blocks for Vermont, with each populated block reduced to a representative point carrying its POP20. The candidates are every OpenStreetMap amenity=school feature in the county — 106 of them. Drive times come from OpenStreetMap's drive network for the county plus 10.5 km, with OSMnx's imputed free-flow speeds, road-only between the nearest nodes. Two blocks (553 people) whose nearest node has no road route to any school are dropped here, before any model sees them — the fix step 2 describes. Solved with spopt 0.7 and PuLP 3.3's bundled CBC; spopt fails on PuLP 4.0, so install spopt "pulp<4". Overpass via the maps.mail.ru mirror while overpass-api.de was refusing connections. The OpenStreetMap answers are read as the Overpass mirror gave them on 10 October 2026.

import io
import os
import requests
import time
import urllib.request
import zipfile

import geopandas as gpd
import numpy as np
import osmnx as ox
import pandas as pd
import scipy.sparse as sp
from scipy.sparse.csgraph import dijkstra

# OpenStreetMap as Overpass returned it on 10 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-location-model-never-solves-2026-10-10.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.overpass_rate_limit = False
ox.settings.http_user_agent = "spatialworkflow.io guide examples (https://www.spatialworkflow.io/)"
CRS = 32145                                               # NAD83 / Vermont, metres

if not os.path.exists("tl_2020_50_tabblock20.zip"):
    urllib.request.urlretrieve("https://www2.census.gov/geo/tiger/TIGER2020/TABBLOCK20/"
                               "tl_2020_50_tabblock20.zip", "tl_2020_50_tabblock20.zip")
vt_blocks = gpd.read_file("tl_2020_50_tabblock20.zip").to_crs(CRS)
county = vt_blocks[vt_blocks.COUNTYFP20 == "007"].union_all()           # Chittenden
homes = vt_blocks[(vt_blocks.COUNTYFP20 == "007") & (vt_blocks.POP20 > 0)].copy()
homes["geometry"] = homes.geometry.representative_point()
homes = homes.reset_index(drop=True)

county_ll = gpd.GeoSeries([county], crs=CRS).to_crs(4326).iloc[0]
schools = ox.features_from_polygon(county_ll, tags={"amenity": "school"}).to_crs(CRS)
schools["geometry"] = schools.geometry.representative_point()
schools = schools.sort_index().reset_index()[["element", "id", "name", "geometry"]]

reach = gpd.GeoSeries([county.buffer(10_500)], crs=CRS).to_crs(4326).iloc[0]
G = ox.add_edge_travel_times(ox.add_edge_speeds(ox.graph_from_polygon(reach, network_type="drive")))
G = ox.project_graph(G, to_crs=CRS)
index = {n: i for i, n in enumerate(G.nodes)}
best = {}
for u, v, d in G.edges(data=True):                      # shortest of any parallel edges
    best[(index[u], index[v])] = min(best.get((index[u], index[v]), np.inf), d["travel_time"])
rows, cols = zip(*best)
T = sp.csr_array((list(best.values()), (rows, cols)), shape=(len(index), len(index)))
s_idx = [index[n] for n in ox.distance.nearest_nodes(G, schools.geometry.x, schools.geometry.y)]
h_all = np.array([index[n] for n in ox.distance.nearest_nodes(G, homes.geometry.x, homes.geometry.y)])
C_all = dijkstra(T.T, directed=True, indices=s_idx)[:, h_all].T

reachable = np.isfinite(C_all).all(axis=1)
print(f"dropping {(~reachable).sum()} blocks ({homes.POP20[~reachable].sum():,} people) with no road route")
homes, h_idx = homes[reachable].reset_index(drop=True), list(h_all[reachable])
C, w = C_all[reachable], homes["POP20"].to_numpy(float)
print(f"C: {C.shape[0]:,} blocks x {C.shape[1]} schools, {w.sum():,.0f} residents")

dropping 2 blocks (553 people) with no road route
C: 2,239 blocks x 106 schools, 167,770 residents

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,239 demand x 106 candidates -> p-median variables 237,440
  smallest radius that covers everyone: 928 s
  radius 600 s: 24 demand points cannot be covered
  cutoff 1200 s keeps 121,710 assignment variables (51.3%)
2,239 demand x 106 candidates -> p-median variables 237,440
  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 county's blocks to its 106 schools, built in the setup above, 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,710 assignment variables, 43.5 s, status Optimal (sol_status 1), sites [11, 24, 52, 62, 71, 72, 77, 98]

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 27.0 minutes — so the 20-minute model could not choose it. Its own optimum was a different set of sites with a 6.8% higher population-weighted travel time (mean trip 5.62 minutes against 5.26 — Example 3 prints both). It was faster, 43.5 s against 87.7 s for the full model, and that is the trap: 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.

The cutoff made the model smaller by forbidding the answer the full model found.

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)

72.4 s: Model is not solved: Not Solved. See ``pulp.constants.LpStatus`` for more information.
87.7 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 long after 20 s — with nothing.

def mean_trip(C, w, sites):
    return np.average(C[:, sites].min(axis=1), weights=w) / 60


beyond = C[:, dense_sites].min(axis=1) > 1200
print(f"full optimum: mean trip {mean_trip(C, w, dense_sites):.2f} min; {beyond.sum()} blocks beyond 20 minutes, "
      f"worst {C[:, dense_sites].min(axis=1).max() / 60:.1f} min")
print(f"20-minute model: mean trip {mean_trip(C, w, sparse_sites):.2f} min "
      f"({mean_trip(C, w, sparse_sites) / mean_trip(C, w, dense_sites) - 1:+.1%})")

full optimum: mean trip 5.26 min; 29 blocks beyond 20 minutes, worst 27.0 min
20-minute model: mean trip 5.62 min (+6.8%)

Example 4 — the measurements behind the steps

These runs take about an hour between them; the timings depend on the machine and on what else it is doing, and so do the outcomes of the solves that run into a time limit.

import pulp
from spopt.locate import LSCP, PCenter

# step 2: a matrix with an unreachable row, or a gap
for label, value in (("inf", np.inf), ("NaN", np.nan)):
    bad = C.copy()
    bad[0, :] = value
    start = time.perf_counter()
    try:
        PMedian.from_cost_matrix(bad, w, p_facilities=8).solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=60))
    except Exception as error:
        print(f"{label}: {type(error).__name__}: {error} ({time.perf_counter() - start:.2f} s)")

# step 3: set covering at 10 minutes
lscp = LSCP.from_cost_matrix(C, service_radius=600)
try:
    lscp.solve(pulp.PULP_CBC_CMD(msg=False))
except RuntimeError as error:
    print(f"set covering, 10 min: {error}")

# 200 candidates: the 106 schools, then 94 populated blocks drawn with a fixed seed
extra = np.random.default_rng(0).choice(len(homes), 94, replace=False)
C200 = np.hstack([C, dijkstra(T.T, directed=True, indices=[h_idx[b] for b in extra])[:, h_idx].T])


def pmedian_run(C, w, p=8, time_limit=300, gap=None, label=""):
    start = time.perf_counter()
    model = PMedian.from_cost_matrix(C, w, p_facilities=p)
    built = time.perf_counter() - start
    try:
        model.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=time_limit, gapRel=gap))
        status = pulp.LpStatus[model.problem.status]
    except RuntimeError as error:
        status = str(error).split(". See")[0]
    n, k = C.shape
    print(f"{label:>14} {n * k + k:>9,} variables  built {built:5.1f} s  total {time.perf_counter() - start:6.1f} s  {status}")


# step 4: the size sweep
for k in (25, 50, 100, 106, 200):
    pmedian_run(C200[:, :k], w, label=f"{k} candidates")
pmedian_run(C200, w, gap=0.05, label="200, gap 5%")                     # step 5

# step 7: assignments only within a cutoff
for cutoff in (1200, 900):
    sparse_p_median(C200, w, 8, cutoff, time_limit=600)

# step 8: prune candidates that are nobody's nearest
nearest_somewhere = np.unique(C200.argmin(axis=1))
print(f"pruning keeps {len(nearest_somewhere)} of 200 candidates")

# step 9: p-centre against p-median on the first 50 candidates
start = time.perf_counter()
centre = PCenter.from_cost_matrix(C200[:, :50], p_facilities=8)
centre.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=900))
print(f"p-centre, 50 candidates: {time.perf_counter() - start:.1f} s, {pulp.LpStatus[centre.problem.status]}")

# step 6: aggregate the blocks to block groups, solve, and evaluate on the blocks
bg = homes["GEOID20"].str[:12].to_numpy()
groups = pd.DataFrame({"bg": bg, "x": homes.geometry.x * w, "y": homes.geometry.y * w, "w": w}).groupby("bg").sum()
gx, gy = groups["x"] / groups["w"], groups["y"] / groups["w"]          # population-weighted centres
g_idx = [index[n] for n in ox.distance.nearest_nodes(G, gx, gy)]
C_bg = dijkstra(T.T, directed=True, indices=s_idx)[:, g_idx].T
start = time.perf_counter()
agg = PMedian.from_cost_matrix(C_bg, groups["w"].to_numpy(), p_facilities=8)
agg.solve(pulp.PULP_CBC_CMD(msg=False))
agg_sites = [j for j, v in enumerate(agg.fac_vars) if v.value() > 0.5]
print(f"{len(groups)} block groups: {time.perf_counter() - start:.2f} s; on the blocks, mean trip "
      f"{mean_trip(C, w, agg_sites):.2f} min against {mean_trip(C, w, dense_sites):.2f} "
      f"({mean_trip(C, w, agg_sites) / mean_trip(C, w, dense_sites) - 1:+.2%})")

inf: PulpError: Cannot multiply variables with NaN/inf values (0.54 s)
NaN: PulpError: Cannot multiply variables with NaN/inf values (0.56 s)
set covering, 10 min: Model is not solved: Infeasible. See ``pulp.constants.LpStatus`` for more information.
 25 candidates    56,000 variables  built   0.4 s  total    5.1 s  Optimal
 50 candidates   112,000 variables  built   0.9 s  total   18.5 s  Optimal
100 candidates   224,000 variables  built   1.7 s  total   64.5 s  Optimal
106 candidates   237,440 variables  built   1.9 s  total   88.6 s  Optimal
200 candidates   448,000 variables  built   3.5 s  total  313.2 s  Model is not solved: Not Solved
   200, gap 5%   448,000 variables  built   3.4 s  total  312.0 s  Model is not solved: Not Solved
cutoff 20 min: 225,313 assignment variables, 360.2 s, status Optimal (sol_status 1), sites [24, 52, 85, 86, 115, 169, 172, 197]
cutoff 15 min: 146,389 assignment variables, 36.9 s, status Infeasible (sol_status -1), sites [16, 21, 42, 61, 62, 70, 71, 76, 123, 134, 183]
pruning keeps 185 of 200 candidates
p-centre, 50 candidates: 52.1 s, Optimal
118 block groups: 0.43 s; on the blocks, mean trip 5.30 min against 5.26 (+0.81%)

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.

The same sites, the same matrix: minimising the worst trip took three times longer on 50 candidates, and on 106 did not finish.

Model files, other solvers and warm starts

Temporary model files are written by PuLP for CBC; a 448,000-variable model produces 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 solves in well under a second where p-median took 160 s here.

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
NetworkX Says No Path Exists Between Two Nodes
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,000 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 312 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 let CBC find a solution for the 200-candidate model, 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 52 s against 18.5 s for p-median on an idle machine, and 988 s against 24 s on a busy one.