Allocating Food Across Cities Over 20 Turns
Imagine a small realm with a single central warehouse and seven cities. Each city has a different population, a different local food production capacity, and a different carrying capacity. Food can only travel along roads, and every road has a transport capacity per turn and a cost per unit shipped. The warehouse is refilled a little every turn, and the ruler has a fixed transport budget for the whole 20-turn horizon.
Goal: decide how much food to ship along every road in every turn so that the total population after 20 turns is as large as possible.
This is a classic multi-period network flow problem with population dynamics attached. Everything in it is linear, so it can be solved exactly as a linear program.
1. The Concrete Example
There is one warehouse and seven cities: Aria, Borea, Cinder, Dune, Ember, Fjord, and Glen.
| City | Initial population (thousand) | Local production (t / turn) | Carrying capacity (thousand) |
|---|---|---|---|
| Aria | 50 | 30 | 120 |
| Borea | 80 | 40 | 150 |
| Cinder | 30 | 10 | 100 |
| Dune | 60 | 20 | 130 |
| Ember | 40 | 50 | 90 |
| Fjord | 20 | 5 | 80 |
| Glen | 70 | 25 | 140 |
- The warehouse starts with 800 t and receives 80 t of new food every turn.
- Only Aria, Borea, and Cinder have a direct road to the warehouse. Dune, Ember, Fjord, and Glen can only be supplied by relaying food through other cities.
- The total transport budget is 6000 cost units.
- No city may fall below 60% of its initial population.
2. Mathematical Model
Notation. Let $T=20$ be the number of turns. Let $i$ index the cities and $e$ index the directed roads. Let $P_i(t)$ be the population of city $i$ at turn $t$, and $F_i(t)$ the food consumed there. Let $x_e(t)$ be the food shipped on road $e$, and $S(t)$ the warehouse stock.
Population dynamics. Population moves toward the level that the available food can support, at adjustment rate $g$ (with $c$ tons of food needed per thousand people):
$$
P_i(t+1) = (1-g),P_i(t) + \frac{g}{c},F_i(t)
$$
Food balance in each city. Consumption cannot exceed local production plus inflow minus outflow (relaying through a city is allowed):
$$
F_i(t) \le \pi_i + \sum_{e \in \text{in}(i)} x_e(t) - \sum_{e \in \text{out}(i)} x_e(t)
$$
Warehouse stock. With refill $R$ per turn:
$$
S(t+1) = S(t) - \sum_{e \in \text{out}(0)} x_e(t) + R
$$
Road capacity, population limits, and budget. With road capacity $u_e$, unit cost $w_e$, carrying capacity $K_i$, floor ratio $\rho$, and budget $B$:
$$
0 \le x_e(t) \le u_e, \qquad \rho,P_i(0) \le P_i(t) \le K_i, \qquad \sum_{t=0}^{T-1}\sum_e w_e,x_e(t) \le B
$$
Objective. A tiny cost penalty $\varepsilon$ removes pointless circular shipments without changing the optimal population:
$$
\max ;; \sum_i P_i(T) ;-; \varepsilon \sum_{t}\sum_e w_e,x_e(t)
$$
The model has about 700 variables, all linear, so the optimum is global. We solve it with SciPy’s HiGHS solver.
3. Full Source Code
1 | import time |
4. Code Walkthrough
4.1 Problem data
T, G_RATE, FOOD_PER_CAP, STOCK0, REFILL, BUDGET, and FLOOR correspond to $T$, $g$, $c$, the initial stock, $R$, $B$, and $\rho$ in the model. The NumPy arrays pop0, prod, and K hold the initial populations, the local production $\pi_i$, and the carrying capacities $K_i$.
Each tuple in roads is (node_a, node_b, capacity, cost), where node 0 is the warehouse. Roads between cities are bidirectional, so the loop that builds arcs creates two directed arcs for each of them. Roads leaving the warehouse are one-way, because food never flows back into the warehouse. Splitting the data into arc_from, arc_to, arc_cap, and arc_cost lets later steps work with vectorized lookups.
4.2 Variable layout
The LP has four families of variables, packed into one long vector:
- $x_e(t)$ for every arc and turn, starting at offset 0.
- $F_i(t)$, the consumed food, starting at
oF. - $P_i(t)$ for $t=0,\dots,T$, starting at
oP. - $S(t)$, the warehouse stock, starting at
oS.
The tiny helpers ix, iF, iP, and iS convert a (row, turn) pair into a column index. All the constraint-building code goes through these four functions, which keeps indexing mistakes out.
4.3 Building the model (speed-up)
build_model() assembles the constraint matrices as sparse COO triplets (row, column, value) and converts them to CSR format.
- Equality constraints encode the population dynamics for each city and turn, and the warehouse stock balance for each turn.
- Inequality constraints encode the per-city food balance, plus one final row for the total budget.
- The objective vector
ccontains $-1$ for each final-turn population, becauselinprogminimizes, plus the small cost penalty $\varepsilon w_e$ on every shipment variable.
The key speed-up is that this model is built only once. The function solve() then only changes the bounds and the single budget entry in b_ub. Bounds also encode the initial populations, the initial stock, and the floor $\rho P_i(0)$. This is why the 64-LP sensitivity sweep later finishes in well under a second. Rebuilding the matrices for each scenario with Python loops would be many times slower. Because the matrices are sparse, HiGHS also solves each LP in milliseconds.
4.4 The solver function
solve() accepts three knobs: the budget, the initial stock, and a road-capacity multiplier. It sets the lower and upper bounds, calls linprog(..., method="highs"), raises a clear error if the model is infeasible, and unpacks the solution vector into flow, F, P, and S. It also returns the total final population and the total transport cost.
4.5 Reporting
The script solves the base scenario, measures the solve time, and prints a summary table. It then runs the sensitivity sweep over 8 budgets and 8 road-capacity multipliers (64 LPs) and stores the final population of each scenario in the matrix Z.
4.6 One figure with six panels
All plots go into a single figure drawn on a 2×3 grid:
- Network map: arrow thickness is the total food shipped over 20 turns, arrow color is road utilization, and node size is the final population.
- Population trajectories: one line per city, with dotted lines for the carrying capacities.
- 3D population surface: city on one axis, turn on the other, population as height.
- Road-utilization heatmap: flow divided by capacity, for every road that is used, on every turn.
- Warehouse stock and cumulative cost: twin axes showing how the stock drains and the budget is consumed.
- 3D sensitivity surface: final population as a function of budget and road-capacity multiplier.
5. Execution Results
============================================================ Optimal plan (relay roads allowed) ============================================================ Solve time : 0.011 s Initial total population : 350.0 thousand Final total population : 443.2 thousand Transport cost used : 6000.0 / 6000 Warehouse stock at turn 20: 529.6 t ------------------------------------------------------------ City Initial Final Capacity Growth % Aria 50.0 109.9 120.0 119.8 Borea 80.0 132.9 150.0 66.1 Cinder 30.0 73.8 100.0 145.9 Dune 60.0 36.0 130.0 -40.0 Ember 40.0 36.7 90.0 -8.3 Fjord 20.0 12.0 80.0 -40.0 Glen 70.0 42.0 140.0 -40.0 ------------------------------------------------------------ Sweep of 64 LPs finished in 0.71 s
6. Interpreting the Results
Overall outcome
Starting from 350 thousand people, the optimal plan reaches about 443 thousand after 20 turns. The whole transport budget of 6000 is consumed, so the budget is the binding constraint. About 530 t of food is still sitting in the warehouse at the end, which tells us that money, not food, is the scarce resource.
Who grows and who does not
The cities with a direct road to the warehouse grow the most: Aria (+120%), Borea (+66%), and Cinder (+146%). These cities offer the cheapest way to turn one ton of food into population. Dune, Fjord, and Glen are pushed all the way down to the 60% floor, and Ember ends slightly below its start. Reaching them requires relaying through other cities at a higher total cost per ton, so the optimizer spends the budget where each unit of cost buys the most population.
The population trajectories
The trajectory plot shows a distinct three-phase shape. Populations first fall to their floors within the first few turns. They then stay flat for roughly ten turns, with total population hovering around 210 thousand. In the last five or six turns they surge upward. The objective only rewards the population at turn 20, so the optimal strategy is to save transport budget during the middle of the horizon and spend it in an end-game burst. The population dynamics have a lag, which is why the surge begins a few turns before the end instead of on the last turn.
The 3D population surface
The 3D surface makes this pattern easy to see. The “plateau” is the long flat stretch where every city sits at its floor. The sharp ridges on the right edge are the cities that receive the end-game shipments.
The road-utilization heatmap
The heatmap confirms the end-game story. The roads from the warehouse to Aria, Borea, and Cinder are almost idle for most of the horizon, and from around turn 16 they run at 100% of their capacity. The relay roads toward Glen and Fjord run at a steady moderate level, because those cities must be kept alive above the floor all the time.
Stock and cost
The warehouse stock climbs steadily while shipments are small, peaks around turn 15, and then collapses as the big shipments begin. The cumulative-cost curve is nearly flat at first and then rises steeply to hit the budget line exactly at turn 20.
The 3D sensitivity surface
The final 3D graph shows the final population for 64 combinations of budget and road capacity. The worst case is only about 311 thousand, and the best is about 635 thousand. Two regimes are visible:
- At low budgets the surface is almost flat in the road-capacity direction. Money is the bottleneck, so widening the roads does not help.
- At high budgets the surface rises steeply with road capacity. Money is no longer the bottleneck, and the roads become the limiting factor.
This tells a planner when investing in new roads pays off and when it is better to simply fund more transport.
7. Summary
- A multi-city, multi-turn food distribution problem can be written as a linear program with flow variables on roads, population variables in cities, and stock variables at the warehouse.
- The sparse model is built once and re-solved with different bounds, so even a 64-scenario sensitivity analysis takes less than a second.
- The optimal plan concentrates resources on cities with cheap access to the warehouse and saves most of the budget for an end-game surge, because only the final population counts.
- The sensitivity surface shows whether money or road capacity is the real bottleneck, which is exactly the kind of insight that makes this type of model useful for planning.
