Choosing Which Roads, Power Lines, and Water Mains to Build
The Problem
A regional government has several candidate infrastructure segments: roads, power-grid lines, and water mains. A single segment can serve several cities at once, but each one costs money. With a fixed budget, which segments should be built so that the total population capacity is as large as possible?
Concrete Example
- There are 10 cities with populations between 45k and 150k.
- There are 14 candidate segments: 5 roads (R1–R5), 5 grid lines (G1–G5), and 4 water mains (W1–W4).
- Each segment has a cost (in $M) and a list of the cities it serves.
- The budget is $B = 120$ ($M).
- A city’s capacity depends on which infrastructure types reach it. The weights are $w_{\text{road}}=0.30$, $w_{\text{grid}}=0.40$, $w_{\text{water}}=0.30$.
- Two segments of the same type reaching the same city do not add up. Coverage counts once per type.
Mathematical Formulation
Let $x_j \in {0,1}$ indicate whether segment $j$ is built, $J_t$ the set of segments of type $t$, and $S_j$ the set of cities served by segment $j$. Let $y_{i,t}$ indicate whether city $i$ is covered by type $t$.
$$
\max_{x,,y} \quad \sum_{i=1}^{10} p_i \sum_{t \in {\text{road},,\text{grid},,\text{water}}} w_t , y_{i,t}
$$
$$
\text{s.t.} \quad \sum_{j} c_j x_j \le B
$$
$$
y_{i,t} \le \sum_{j \in J_t :, i \in S_j} x_j, \qquad y_{i,t} \le 1, \qquad x_j \in {0,1}, \quad y_{i,t}\in{0,1}
$$
This is a budgeted maximum coverage problem, which is NP-hard in general. With only 14 candidates there are $2^{14} = 16{,}384$ subsets, so exhaustive search is exact and fast, provided it is vectorized with NumPy. A greedy method based on marginal gain per cost is included as a baseline.
Source Code
1 | import numpy as np |
Code Walkthrough
1. Problem Data
coords and population define the 10 cities. segments is a list of tuples (name, type, cost, served cities), which keeps the whole problem in one readable table. The coverage tensor A[t][j, i] is 1 when segment $j$ of type $t$ serves city $i$. Keeping one matrix per type lets us test “is city $i$ covered by type $t$?” with a single matrix product.
2. Vectorized Exhaustive Search
The key to speed is avoiding a Python loop over $2^{14}$ subsets.
masks = np.arange(2**n_seg)lists every subset as an integer.bitsexpands each integer into a 0/1 row, so row $m$ tells which segments are built in subset $m$. Its shape is $(16384, 14)$.cost_all = bits @ seg_costgives the cost of every subset at once.- For each type,
(bits @ A[t]) > 0is a Boolean matrix of shape $(16384, 10)$. Entry $(m, i)$ is true if at least one built segment of type $t$ serves city $i$. This is exactly the “coverage counts once” rule, because the count is thresholded at zero. cov @ populationconverts coverage to covered population per type, andpop_by_type @ w_vecproduces the weighted capacity for all subsets:
$$
\text{Cap}(m) = \sum_{t} w_t \sum_{i} p_i \cdot \mathbb{1}\Big[\textstyle\sum_{j \in J_t,, i\in S_j} x_j^{(m)} \ge 1\Big]
$$
Infeasible subsets (cost above the budget) are masked with $-1$, and argmax picks the best one. The bit pattern of best_idx is then decoded into the list of chosen segments. The whole search takes a fraction of a second.
3. Greedy Baseline
capacity_of evaluates any selection directly. greedy repeatedly adds the segment with the highest ratio
$$
\frac{\Delta \text{Cap}(j)}{c_j}
$$
among those that still fit in the remaining budget, and stops when nothing improves. Greedy is intuitive and fast, but it can lock in cheap high-ratio segments that later block a better combination.
4. Console Report
The script prints the optimal and greedy selections, their costs, the greedy-to-optimal ratio, and a per-city table showing which of the three infrastructure types reach each city and how much weighted capacity it contributes.
5. Sensitivity Data
opt_curveis the best achievable capacity for 41 budget levels. Because the capacity of every subset is already computed, each point is just a maskedmax.Zis a surface over two parameters: the budget and the weight of the power grid. The road and water weights share the remainder equally, $w_{\text{road}} = w_{\text{water}} = (1-w_{\text{grid}})/2$. Sincepop_by_typeis reused, changing the weights costs only a matrix-vector product, so the 22 × 22 surface is computed without re-searching subsets.
6. Combined Figure
A single GridSpec figure contains four panels:
- (a) A map of cities (marker size proportional to population) with built segments drawn as thick solid lines and skipped segments as thin dashed lines. Each segment is drawn as a hub-and-spoke from its centroid to the cities it serves.
- (b) A 3D bar chart with city on one axis, infrastructure type on another, and weighted capacity as height.
- (c) The budget–capacity curve for the optimal and greedy methods, with the chosen budget marked.
- (d) A 3D surface of optimal capacity over budget and grid weight.
Execution Results
============================================================
Budget: 120 $M | Candidates: 14 | Subsets searched: 16,384
============================================================
[Optimal selection]
Segments : R3, G1, W1, W4
Cost : 116 $M
Capacity : 492.00 (thousand-weighted persons)
[Greedy selection]
Segments : W4, G1, W1, G4
Cost : 107 $M
Capacity : 463.00
Greedy / Optimal = 94.11 %
[City-level coverage under the optimal plan]
Population(k) Road Grid Water Weighted capacity
Aoba 120 0 1 1 84.0
Bell 85 0 1 1 59.5
Cedar 60 0 0 0 0.0
Dune 95 1 1 1 95.0
Elm 150 1 1 1 150.0
Fjord 70 0 0 0 0.0
Glen 55 1 0 1 33.0
Haven 110 1 0 0 33.0
Iris 80 0 0 1 24.0
Jade 45 0 0 1 13.5
Interpreting the Results
Optimal Plan vs. Greedy
With a budget of $120M, the exhaustive search selects R3, G1, W1, and W4 at a cost of $116M, reaching a weighted capacity of 492.0. The greedy method selects W4, G1, W1, and G4 at $107M and reaches only 463.0, about 94.1 % of the optimum.
The reason is visible in the selection. Greedy fills the budget with segments that look efficient one at a time, including the cheap grid line G4. The exhaustive search instead spends the remaining money on the road R3. R3 reaches four cities, including Elm (150k), Dune (95k), and Haven (110k), so it unlocks a large block of new coverage that no single cheap segment can match. This is the classic weakness of ratio-based greedy rules under a hard budget: the best final combination is not the sequence of locally best choices.
Panel (a): Network Map
The solid lines show where the money goes. The grid line G1 and the water main W1 both pull in the northwestern cluster around Aoba, Bell, and Dune. W4 stretches toward Elm, Iris, and Jade, while R3 links the central and southern cities. Cities such as Cedar and Fjord, which are far from the cheap clusters and have modest populations, receive nothing. This is the practical consequence of a budget constraint: capacity-per-dollar decides who gets served.
Panel (b): 3D Capacity by City and Type
The tallest stacks belong to Elm and Dune, which are covered by all three types and therefore receive their full weighted capacity. Cedar and Fjord show no bars at all. Cities with a single infrastructure type, such as Iris and Jade (water only) or Haven (road only), contribute small bars. The chart makes it clear that the plan concentrates investment where multiple infrastructure types overlap, because every additional type multiplies the benefit for a populous city.
Panel (c): Budget vs. Capacity
Both curves rise steeply at low budgets, where the cheapest high-impact segments are bought, and flatten at high budgets, where remaining segments mostly duplicate existing coverage. This diminishing-returns shape is typical of coverage problems. The optimal curve is never below the greedy curve, and the gap appears at specific budgets where greedy commits to the wrong segment. The marked point at $120M is one of those budgets, which is why it was chosen for this example.
Panel (d): 3D Surface of Budget × Grid Weight
The surface rises monotonically along the budget axis, reflecting the diminishing-returns curve from panel (c). Along the grid-weight axis the optimal capacity also increases, because a larger share of the weight goes to the grid, whose segments (G1, G3) cover the most populous cities. The slope is steepest in the low-budget region, which means the value of an extra dollar is highest when the budget is tight. This surface can be used for policy questions such as “how much would an additional $20M change the outcome if power reliability is valued more?” Such questions can be answered by reading the surface instead of rerunning a full optimization.
Summary
- The infrastructure prioritization problem is a budgeted maximum coverage problem with a coverage-once rule per infrastructure type.
- Vectorizing all $2^{14}$ subsets with NumPy bit-expansion turns an exponential-looking search into a near-instant matrix computation.
- A greedy ratio rule is quick but can fall several percent short of the optimum, as the 94.1 % result shows.
- Budget and weight sensitivity can be explored cheaply by reusing the precomputed per-subset coverage, which produces the 3D surface in a few lines.
