Maximizing 20-Turn Cumulative Population with Python

November 4, 2026

How Much Should a Region Spend on Disaster Preparedness?

Floods and earthquakes destroy houses and factories. Every turn, a planner must decide how much of the budget goes to prevention, how much to stockpiles, and how much to keep as repair materials. In this article I build a small simulation of this problem, optimize the budget split, and look at the result.

1. Problem Setup

A region starts with the following state.

Item	Symbol	Initial value
Population	$P_0$	100
Housing capacity	$H_0$	130
Production capacity	$K_0$	70

Each turn, output is limited by labor and production capacity:

$$Y_t = \min\left(K_t,\ \lambda P_t\right), \qquad \lambda = 0.6$$

A fixed share of output becomes the budget:

$$B_t = s,Y_t, \qquad s = 0.25$$

The budget is split into four uses.

Share	Use	Effect
$a$	Disaster defenses (levees, seismic retrofits)	Reduce damage
$b$	Relief stockpile (food, water, tents)	Reduces deaths among the homeless
$c$	Repair materials	Used to rebuild after a disaster
$d = 1-a-b-c$	Growth investment	Adds housing and production capacity

A parameter $w$ decides which assets get the repair materials first: a share $w$ goes to housing and $1-w$ to production facilities.

Stocks and depreciation. Defenses $D_t$, relief supplies $S_t$, and repair materials $M_t$ decay a little each turn:

$$D_{t+1} = 0.97,D_t + a B_t,\quad S_{t+1} = 0.95,S_t + b B_t,\quad M_{t+1} = 0.98,M_t + c B_t$$

Disasters. A flood occurs with probability 0.15 per turn and an earthquake with probability 0.04. Each has a random intensity $u \in [0.3, 1]$. Defenses reduce the damage rate with diminishing returns:

$$\delta^{F} = 0.45,u_F\left(1 - 0.85\left(1 - e^{-D/30}\right)\right),\qquad \delta^{E} = 0.70,u_E\left(1 - 0.60\left(1 - e^{-D/50}\right)\right)$$

Combined, the fraction of assets lost is

$$\phi = 1 - (1-\delta^{F})(1-\delta^{E})$$

Housing falls by $\phi H$ and production capacity by $0.8,\phi K$. The lost amount becomes repair debt, which repair materials can rebuild. Direct deaths are $P,(0.04,\delta^{F} + 0.12,\delta^{E})$.

Homeless people. People without a home number $\max(0, P - H)$. They need supplies, and unmet needs raise mortality:

$$\text{mortality} = 0.01 + 0.08\times(\text{unmet share of needs})$$

Recovery. Population grows by up to 2.5% per turn, limited by free housing.

Objective. The goal is to choose $(a, b, c, w)$ to maximize the expected cumulative population:

$$\max_{a,b,c,w}\ \ \mathbb{E}\left[\sum_{t=1}^{20} P_t\right]$$

The expectation is estimated by Monte Carlo over random disaster scenarios, and optimized with differential evolution, a global optimizer that needs no gradients.

2. Source Code

The simulator advances all scenarios at once with NumPy arrays instead of looping over scenarios. This makes the optimizer run in a few seconds.

import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
from scipy.optimize import differential_evolution

plt.rcParams["figure.dpi"] = 110

# ============================================================
# 1. Parameters
# ============================================================
T = 20                       # number of turns
P0, H0, K0 = 100.0, 130.0, 70.0   # initial population, housing capacity, production capacity
LABOR = 0.6                  # output per person (capped by production capacity)
SAVE_RATE = 0.25             # share of output available as budget each turn

P_FLOOD, P_QUAKE = 0.15, 0.04          # per-turn disaster probabilities
FLOOD_DMG, QUAKE_DMG = 0.45, 0.70      # asset damage at intensity 1.0
FLOOD_MORT, QUAKE_MORT = 0.04, 0.12    # direct mortality at intensity 1.0

D_DECAY, S_DECAY, M_DECAY = 0.97, 0.95, 0.98   # depreciation of defenses, supplies, materials
D_FLOOD, D_QUAKE = 30.0, 50.0          # saturation scales of defense effect
FLOOD_MAX, QUAKE_MAX = 0.85, 0.60      # maximum damage reduction by defenses

NEED = 0.15                  # supply need per homeless person per turn
HL_MORT, HL_MORT_UNMET = 0.01, 0.08    # mortality of homeless people
GROWTH = 0.025               # natural growth (limited by free housing)

# ============================================================
# 2. Disaster scenarios (generated once, vectorized)
# ============================================================
def make_scenarios(n, seed):
    rng = np.random.default_rng(seed)
    flood = (rng.random((n, T)) < P_FLOOD) * rng.uniform(0.3, 1.0, (n, T))
    quake = (rng.random((n, T)) < P_QUAKE) * rng.uniform(0.3, 1.0, (n, T))
    return flood, quake

# ============================================================
# 3. Simulator (all scenarios advance simultaneously with NumPy)
# ============================================================
def simulate(a, b, c, w, sc, path=False):
    """a: defense, b: relief stock, c: repair materials, rest d: growth investment,
    w: share of repair materials sent to housing."""
    flood, quake = sc
    n = flood.shape[0]
    d = max(0.0, 1.0 - a - b - c)
    P = np.full(n, P0); H = np.full(n, H0); K = np.full(n, K0)
    D = np.zeros(n); S = np.zeros(n); M = np.zeros(n)
    Hd = np.zeros(n); Kd = np.zeros(n)
    cum = np.zeros(n)
    hist = [P.copy()]
    for t in range(T):
        Y = np.minimum(K, LABOR * P)
        B = SAVE_RATE * Y
        D = D * D_DECAY + a * B
        S = S * S_DECAY + b * B
        M = M * M_DECAY + c * B
        H = H + 0.5 * d * B
        K = K + 0.5 * d * B
        # repair with stocked materials
        rh = np.minimum(Hd, w * M)
        rk = np.minimum(Kd, (1 - w) * M)
        left = M - rh - rk
        xh = np.minimum(Hd - rh, left); left = left - xh
        xk = np.minimum(Kd - rk, left); left = left - xk
        H = H + rh + xh; Hd = Hd - rh - xh
        K = K + rk + xk; Kd = Kd - rk - xk
        M = left
        # disasters
        fp = FLOOD_MAX * (1 - np.exp(-D / D_FLOOD))
        qp = QUAKE_MAX * (1 - np.exp(-D / D_QUAKE))
        fh = FLOOD_DMG * flood[:, t] * (1 - fp)
        qh = QUAKE_DMG * quake[:, t] * (1 - qp)
        frac = 1 - (1 - fh) * (1 - qh)
        lossH = frac * H; lossK = 0.8 * frac * K
        H = H - lossH; Hd = Hd + lossH
        K = K - lossK; Kd = Kd + lossK
        P = P - P * (FLOOD_MORT * fh + QUAKE_MORT * qh)
        # homeless people and relief supplies
        hl = np.maximum(0.0, P - H)
        need = NEED * hl
        used = np.minimum(S, need)
        S = S - used
        unmet = np.where(need > 0, 1 - used / np.maximum(need, 1e-9), 0.0)
        P = P - hl * (HL_MORT + HL_MORT_UNMET * unmet)
        # recovery growth
        P = P + np.minimum(GROWTH * P, np.maximum(0.0, H - P))
        cum += P
        hist.append(P.copy())
    return (cum, np.array(hist)) if path else cum

def decode(x):
    a = x[0]
    b = (1 - a) * x[1]
    c = (1 - a - b) * x[2]
    return a, b, c, x[3]

# ============================================================
# 4. Optimization (differential evolution on Monte Carlo mean)
# ============================================================
train = make_scenarios(1500, seed=0)
test = make_scenarios(5000, seed=1)

t0 = time.time()
res = differential_evolution(
    lambda x: -simulate(*decode(x), train).mean(),
    bounds=[(0, 1)] * 4, seed=0, maxiter=40, popsize=12, tol=1e-6, polish=False)
a, b, c, w = decode(res.x)
d = 1 - a - b - c
print(f"Optimization time: {time.time() - t0:.1f} s")
print("=== Optimal allocation ===")
print(f"Defense a={a:.3f} | Relief stock b={b:.3f} | Repair materials c={c:.3f} | Growth d={d:.3f}")
print(f"Share of repair materials to housing w={w:.3f}")

# ============================================================
# 5. Strategy comparison on unseen scenarios
# ============================================================
strategies = {
    "No disaster spend": (0.0, 0.0, 0.0, 0.5),
    "Reactive only": (0.0, 0.0, 0.30, 0.5),
    "Balanced": (0.15, 0.10, 0.20, 0.5),
    "Optimized": (a, b, c, w),
}
results, paths = {}, {}
print("\n=== Cumulative population over 20 turns (5,000 test scenarios) ===")
print(f"{'Strategy':<20}{'Mean':>10}{'P10':>10}{'P90':>10}")
for name, params in strategies.items():
    cum, hist = simulate(*params, test, path=True)
    results[name], paths[name] = cum, hist
    print(f"{name:<20}{cum.mean():>10.1f}{np.percentile(cum, 10):>10.1f}{np.percentile(cum, 90):>10.1f}")

calm = (np.zeros((1, T)), np.zeros((1, T)))
print(f"\nDisaster-free benchmark (growth-only spend): {simulate(0, 0, 0, 0.5, calm)[0]:.1f}")
gain = results["Optimized"].mean() / results["No disaster spend"].mean() - 1
print(f"Improvement over 'No disaster spend': {gain * 100:.1f}%")

# ============================================================
# 6. Grid scan for the 3D surface (defense a vs repair materials c)
# ============================================================
grid = np.linspace(0, 0.9, 25)
Z = np.full((len(grid), len(grid)), np.nan)
t0 = time.time()
for i, ai in enumerate(grid):
    for j, cj in enumerate(grid):
        if ai + b + cj <= 1.0:
            Z[j, i] = simulate(ai, b, cj, w, train).mean()
print(f"\nGrid scan time: {time.time() - t0:.1f} s")
A, C = np.meshgrid(grid, grid)
Zm = np.ma.masked_invalid(Z)

# ============================================================
# 7. Visualization (single figure)
# ============================================================
names = list(strategies.keys())
colors = {"No disaster spend": "#9e9e9e", "Reactive only": "#e39b2d",
          "Balanced": "#3b82c4", "Optimized": "#c0392b"}
turns = np.arange(T + 1)

fig = plt.figure(figsize=(18, 16))
gs = GridSpec(3, 2, figure=fig, hspace=0.32, wspace=0.22)

ax1 = fig.add_subplot(gs[0, 0])
for nm in names:
    h = paths[nm]
    ax1.plot(turns, h.mean(axis=1), color=colors[nm], lw=2.4, label=nm)
    if nm in ("No disaster spend", "Optimized"):
        ax1.fill_between(turns, np.percentile(h, 10, axis=1),
                         np.percentile(h, 90, axis=1), color=colors[nm], alpha=0.15)
ax1.set_title("Population trajectory (mean and 10-90% band)")
ax1.set_xlabel("Turn"); ax1.set_ylabel("Population"); ax1.legend(); ax1.grid(alpha=0.3)

ax2 = fig.add_subplot(gs[0, 1])
means = [results[n].mean() for n in names]
p10 = [np.percentile(results[n], 10) for n in names]
bars = ax2.bar(names, means, color=[colors[n] for n in names])
ax2.scatter(names, p10, color="black", zorder=3, label="10th percentile")
for bar, m in zip(bars, means):
    ax2.text(bar.get_x() + bar.get_width() / 2, m + 15, f"{m:.0f}", ha="center", fontweight="bold")
ax2.set_ylim(0, max(means) * 1.15)
ax2.set_title("Cumulative population over 20 turns")
ax2.set_ylabel("Person-turns"); ax2.legend(); ax2.grid(axis="y", alpha=0.3)

ax3 = fig.add_subplot(gs[1, 0], projection="3d")
surf = ax3.plot_surface(A, C, Zm, cmap="viridis", edgecolor="none", alpha=0.95)
ax3.scatter([a], [c], [simulate(a, b, c, w, train).mean()], color="red", s=80, depthshade=False)
ax3.set_xlabel("Defense share a"); ax3.set_ylabel("Repair-material share c")
ax3.set_zlabel("Mean cumulative population")
ax3.set_title("3D surface: defense vs repair materials (red = optimum)")
ax3.view_init(elev=28, azim=-125)
fig.colorbar(surf, ax=ax3, shrink=0.55, pad=0.08)

ax4 = fig.add_subplot(gs[1, 1])
cs = ax4.contourf(A, C, Zm, levels=20, cmap="viridis")
ax4.scatter([a], [c], color="red", s=90, label="Optimum", zorder=3)
ax4.set_xlabel("Defense share a"); ax4.set_ylabel("Repair-material share c")
ax4.set_title("Contour view of the same surface"); ax4.legend()
fig.colorbar(cs, ax=ax4)

ax5 = fig.add_subplot(gs[2, 0])
for nm in ("No disaster spend", "Reactive only", "Optimized"):
    ax5.hist(results[nm], bins=60, alpha=0.55, color=colors[nm], label=nm)
ax5.set_title("Distribution of outcomes across disaster scenarios")
ax5.set_xlabel("Cumulative population"); ax5.set_ylabel("Scenarios")
ax5.legend(); ax5.grid(alpha=0.3)

ax6 = fig.add_subplot(gs[2, 1], projection="3d")
sel = np.linspace(0, 1, 40)
sorted_idx = np.argsort(results["Optimized"])
pick = sorted_idx[(sel * (len(sorted_idx) - 1)).astype(int)]
for k, idx in enumerate(pick):
    ax6.plot(turns, np.full(T + 1, k), paths["Optimized"][:, idx], color=plt.cm.plasma(k / 39), lw=1.2)
ax6.set_xlabel("Turn"); ax6.set_ylabel("Scenario rank (low to high)")
ax6.set_zlabel("Population")
ax6.set_title("3D trajectories of 40 scenarios (Optimized)")
ax6.view_init(elev=24, azim=-60)

fig.suptitle("Disaster preparedness and population maintenance", fontsize=17, y=0.94)
plt.show()

3. Code Walkthrough

Parameters (Section 1)

All constants live at the top, so you can change the disaster frequency or damage level and rerun. FLOOD_MAX and QUAKE_MAX cap how much defenses can help: even a very large investment cannot erase damage entirely.

Scenario generation (Section 2)

make_scenarios draws a Bernoulli event (flood or earthquake) for every scenario and turn, then multiplies it by a uniform intensity in $[0.3, 1]$. The result is two arrays of shape (n_scenarios, T). Two independent sets are created:

train (1,500 scenarios) is used inside the optimizer.
test (5,000 scenarios) is used only for the final comparison.

Evaluating on scenarios the optimizer never saw avoids overstating the benefit of the optimal plan.

The simulator (Section 3)

The state variables (P, H, K, D, S, M, Hd, Kd) are NumPy arrays with one entry per scenario. Each turn performs these steps:

Production and budget. Output is min(K, LABOR * P), and a fixed share becomes the budget B.
Stock updates. Defenses, relief supplies, and repair materials decay and then receive their share of the budget. Growth investment is split equally between housing and production capacity.
Repairs. Materials are sent first to the assets that the parameter $w$ prioritizes. Any leftover material is used on the other asset. Repaired assets reduce the repair debt Hd and Kd.
Disasters. The defense stock determines the damage reduction through $1 - e^{-D/\text{scale}}$. Housing and production capacity are reduced, and the loss is added to the repair debt.
Homelessness and relief. People beyond housing capacity are homeless. They consume relief supplies, and the unmet share raises their mortality.
Recovery. Population grows at up to 2.5%, but only into free housing.

Because every operation uses np.minimum, np.maximum, and array arithmetic, there is only one Python loop over the 20 turns. This is the speed-up: a naive version would loop over scenarios as well, which would be roughly 1,500 times slower per evaluation.

Parameter decoding (`decode`)

The optimizer works in the unit cube $[0,1]^4$. The function converts these to valid shares:

$$a = x_0,\quad b = (1-a),x_1,\quad c = (1-a-b),x_2$$

This guarantees $a + b + c \le 1$ without needing explicit constraints.

Optimization (Section 4)

differential_evolution evolves a population of candidate allocations. It handles the non-smooth min/max structure and the noisy Monte Carlo objective well. polish=False skips the final gradient-based refinement, which is unhelpful for a noisy objective.

Strategy comparison (Section 5)

Four strategies are evaluated on the same 5,000 test scenarios:

No disaster spend: everything goes to growth.
Reactive only: 30% goes to repair materials and nothing to prevention.
Balanced: a hand-made mix of all four.
Optimized: the optimizer’s answer.

The disaster-free benchmark (no disasters at all) gives an upper bound on what any strategy can reach.

3D surface scan (Section 6)

With $b$ and $w$ fixed at their optimal values, the code scans a $25 \times 25$ grid of defense share $a$ and repair-material share $c$. Cells where $a + b + c > 1$ are masked out, which is why the contour plot is triangular.

4. Execution Results

Optimization time: 9.7 s
=== Optimal allocation ===
Defense a=0.453 | Relief stock b=0.004 | Repair materials c=0.346 | Growth d=0.197
Share of repair materials to housing w=1.000

=== Cumulative population over 20 turns (5,000 test scenarios) ===
Strategy                  Mean       P10       P90
No disaster spend       2186.8    1684.8    2591.6
Reactive only           2192.4    1693.6    2589.0
Balanced                2310.8    1857.6    2604.7
Optimized               2432.8    2079.5    2610.2

Disaster-free benchmark (growth-only spend): 2618.3
Improvement over 'No disaster spend': 11.3%

Grid scan time: 1.5 s

5. Reading the Results

Optimal allocation

The optimizer found a clear pattern:

Use	Share
Defenses $a$	about 45%
Relief stockpile $b$	almost 0%
Repair materials $c$	about 35%
Growth investment $d$	about 20%
Housing share of repair materials $w$	about 100%

Three points stand out.

Prevention and repair work together. About 80% of the budget goes to defenses and repair materials, while growth gets only about 20%. Putting everything into growth is poor because a single large disaster erases years of expansion.

Housing comes first. The optimum sends nearly all repair materials to housing. The reason is the homelessness penalty: a person without a home dies at a higher rate and cannot be supported by population growth, while lost production capacity only reduces the budget. Rebuilding houses protects the population directly.

Relief stockpiles are almost unused. In this model, defenses and fast rebuilding keep the number of homeless people small, so large stockpiles are rarely drawn on. This depends on the parameters: if rebuilding were slow or homelessness widespread, $b$ would rise.

Strategy comparison

Strategy	Mean	10th percentile	90th percentile
No disaster spend	2186.8	1684.8	2591.6
Reactive only	2192.4	1693.6	2589.0
Balanced	2310.8	1857.6	2604.7
Optimized	2432.8	2079.5	2610.2

The optimized plan improves the mean by about 11.3% over spending nothing on disasters, and it lifts the 10th percentile by more than 20%. The gain is therefore concentrated in the bad scenarios, which is what preparedness is for. The disaster-free benchmark is about 2618, so the optimized plan recovers most of the gap between the unprepared case and a world with no disasters.

Spending only on repair materials (“Reactive only”) barely helps. Without defenses, damage is large and materials are used up quickly.

Figure guide

Population trajectory (top left). The unprepared region’s lower band drifts below the starting population after repeated disasters, while the optimized path grows steadily and its band stays above 100.
Cumulative population bars (top right). Mean values rise from left to right, and the black dots (10th percentile) rise even faster, showing improved worst-case protection.
3D surface (middle left). The surface is a smooth hill. Too little defense leaves the region exposed; too much starves growth and repair. The red dot marks the peak.
Contour view (middle right). The same hill seen from above. The high plateau is fairly wide, which means the plan is robust: allocations near the optimum perform almost as well.
Outcome distribution (bottom left). Without preparation, outcomes spread to the left with a long tail of poor scenarios. The optimized histogram is concentrated near the best outcomes.
3D trajectories (bottom right). Forty scenarios sorted from worst to best. Even the worst runs under the optimized plan avoid collapse, and most climb steadily.

6. Conclusion

In this model, a balanced allocation of roughly 45% to defenses, 35% to repair materials, and 20% to growth maximizes the expected 20-turn cumulative population, with nearly all repair materials directed to housing. The main lesson is that disaster spending is not a cost that competes with growth. It protects the very assets that growth builds, and its biggest payoff appears in the worst scenarios.

To extend the model, you could add more disaster types, make rebuilding speed depend on budget, or let the allocation change over time as a policy rather than a constant.

Maximize Its Population in 30 Years

November 3, 2026

How Much Should a Nation Invest in Irrigation and Grain Reserves to Maximize Its Population in 30 Years?

Rainfall changes from year to year, and when a drought hits, agricultural output drops sharply. If we plan only around the average harvest, one bad stretch of weather can set a population back for years. In this post we build a small model that puts two kinds of preparation on the table at once:

Irrigation investment, which permanently reduces how much a drought hurts the harvest.
Grain reserves, which are consumed in lean years but lose some value to spoilage every year.

We then use Monte Carlo simulation to find the policy that maximizes the expected population after 30 years.

1. The Model

Rainfall. The rainfall index in year $t$ follows a Gamma distribution with mean 1 and coefficient of variation $\mathrm{CV}=0.45$:

$$
R_t \sim \mathrm{Gamma}\left(k,\ \tfrac{1}{k}\right),\qquad k=\frac{1}{\mathrm{CV}^2}
$$

Harvest. Let $I_t\in[0,1]$ be the irrigation coverage. Irrigation covers a fraction $\varepsilon$ of the rainfall shortfall, and keeping it running costs $u$ per unit of coverage:

$$
y_t = \min(1,R_t) + \bigl(1-\min(1,R_t)\bigr),\varepsilon, I_t,\qquad
F_t = F_{\max}, y_t - u, I_t
$$

Surplus and deficit. With population $P_t$ (food demand is one unit per person):

$$
\text{Surplus}_t=\max(F_t-P_t,,0),\qquad \text{Deficit}_t=\max(P_t-F_t,,0)
$$

Policy. Two parameters define the policy:

$a\in[0,1]$ is the share of the surplus invested in irrigation.
$s\ge 0$ is the target reserve, measured in years of consumption.

Irrigation dynamics. Here $C$ is the cost of covering the whole country:

$$
I_{t+1}=\min!\left(1,\ I_t+\frac{\mathrm{Inv}_t}{C}\right),\qquad \mathrm{Inv}_t=\min\bigl(a\cdot \text{Surplus}_t,\ (1-I_t),C\bigr)
$$

Reserve dynamics. In a deficit year the reserve is drawn down first. Surplus left after the irrigation investment is stored up to the target $sP_t$. A spoilage rate $\delta$ is applied each year:

$$
S_{t+1}=(1-\delta)\Bigl(S_t-\text{Draw}_t+\text{Store}_t\Bigr)
$$

Population growth. If a shortfall remains after drawing on the reserve, the population shrinks in proportion to the unmet demand. Otherwise it grows at up to $g$, depending on how much food is left over:

$$
P_{t+1}=P_t,(1+r_t),\qquad
r_t=
\begin{cases}
-m,\dfrac{\text{Shortage}_t}{P_t} & \text{if Shortage}_t>0\[3mm]
g\cdot\min!\left(1,\ \dfrac{\text{Leftover}_t}{0.1,P_t}\right) & \text{otherwise}
\end{cases}
$$

Objective. We search for the policy that maximizes the expected population at year 30:

$$
\max_{a,,s}\ \ \mathbb{E}\bigl[P_{30}\bigr]
$$

The trade-off is that surplus spent on irrigation or stored in the reserve is surplus that cannot go toward population growth today. The right balance between the two is what we want to find.

2. Full Source Code

The search uses 3,000 rainfall scenarios, and the chosen policy is then evaluated on 5,000 fresh scenarios it has never seen.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import time

# ---------------------------------------------------------------
# 1. Parameters
# ---------------------------------------------------------------
T = 30                # planning horizon (years)
N_TRAIN = 3000        # scenarios used to search for the best policy
N_TEST = 5000         # fresh scenarios used to evaluate the chosen policy
P0 = 100.0            # initial population (millions)
F_MAX = 130.0         # food output with perfect water supply (million "person-years" of food)
G = 0.025             # maximum annual population growth rate
M = 0.5               # mortality coefficient when food is short
SPOIL = 0.05          # annual storage loss rate
COST_FULL = 350.0     # food units needed to irrigate 100% of the land
EFF = 0.5             # share of the rainfall shortfall that full irrigation can cover
UPKEEP = 2.5          # annual upkeep per 100% irrigation coverage (food units)
CV = 0.45             # coefficient of variation of annual rainfall
DROUGHT = 0.70        # rainfall index below this value is called a drought

def make_rain(n, seed):
    rng = np.random.default_rng(seed)
    k = 1.0 / CV**2
    return rng.gamma(shape=k, scale=1.0 / k, size=(n, T))   # mean = 1

# ---------------------------------------------------------------
# 2. Vectorised simulator (all scenarios are advanced at once)
# ---------------------------------------------------------------
def simulate(a, s, R, track=False):
    n = R.shape[0]
    P = np.full(n, P0)
    I = np.zeros(n)
    S = np.zeros(n)
    if track:
        hP = np.zeros((T + 1, n)); hI = np.zeros((T + 1, n)); hS = np.zeros((T + 1, n))
        hP[0] = P
    for t in range(T):
        rr = np.minimum(1.0, R[:, t])
        y = rr + (1.0 - rr) * EFF * I
        prod = F_MAX * y - UPKEEP * I
        surplus = np.maximum(prod - P, 0.0)
        deficit = np.maximum(P - prod, 0.0)
        draw = np.minimum(S, deficit)
        S = S - draw
        inv = np.minimum(a * surplus, (1.0 - I) * COST_FULL)
        I = np.minimum(1.0, I + inv / COST_FULL)
        room = np.maximum(s * P - S, 0.0)
        store = np.minimum(surplus - inv, room)
        S = S + store
        leftover = surplus - inv - store
        shortage = deficit - draw
        r = np.where(shortage > 0.0,
                     -M * shortage / P,
                     G * np.minimum(1.0, leftover / (0.1 * P)))
        S = S * (1.0 - SPOIL)
        P = P * (1.0 + r)
        if track:
            hP[t + 1] = P; hI[t + 1] = I; hS[t + 1] = S
    if track:
        return P, hP, hI, hS
    return P

# ---------------------------------------------------------------
# 3. Grid search over (a, s)
# ---------------------------------------------------------------
R_train = make_rain(N_TRAIN, seed=1)
R_test = make_rain(N_TEST, seed=2)

a_grid = np.linspace(0.0, 1.0, 21)
s_grid = np.linspace(0.0, 2.0, 17)
Z = np.zeros((len(a_grid), len(s_grid)))

t0 = time.time()
for i, a in enumerate(a_grid):
    for j, s in enumerate(s_grid):
        Z[i, j] = simulate(a, s, R_train).mean()
elapsed = time.time() - t0

i_best, j_best = np.unravel_index(np.argmax(Z), Z.shape)
a_best, s_best = a_grid[i_best], s_grid[j_best]

# ---------------------------------------------------------------
# 4. Out-of-sample evaluation
# ---------------------------------------------------------------
P_base, hP_base, hI_base, hS_base = simulate(0.0, 0.0, R_test, track=True)
P_opt, hP_opt, hI_opt, hS_opt = simulate(a_best, s_best, R_test, track=True)
P_irr = simulate(a_best, 0.0, R_test)
P_stk = simulate(0.0, s_best, R_test)

drought_prob = (R_test < DROUGHT).mean()

print("=== Settings ===")
print(f"Horizon: {T} years, scenarios (search/test): {N_TRAIN}/{N_TEST}")
print(f"Drought probability per year (rain < {DROUGHT}): {drought_prob:.3f}")
print(f"Grid search: {len(a_grid) * len(s_grid)} policies in {elapsed:.2f} s")
print()
print("=== Optimal policy ===")
print(f"Share of surplus invested in irrigation a* = {a_best:.2f}")
print(f"Target reserve s* = {s_best:.3f} years of consumption")
print()
print("=== Population after 30 years (millions, test scenarios) ===")
print(f"{'Policy':<28}{'Mean':>9}{'5th pct':>10}{'Median':>9}{'Collapse<50':>13}")
rows = [("No preparation (a=0, s=0)", P_base),
        ("Irrigation only", P_irr),
        ("Stockpile only", P_stk),
        ("Optimal combination", P_opt)]
for name, p in rows:
    print(f"{name:<28}{p.mean():>9.1f}{np.percentile(p, 5):>10.1f}"
          f"{np.median(p):>9.1f}{(p < 50).mean():>12.1%}")
gain = P_opt.mean() / P_base.mean() - 1
print()
print(f"Gain of optimal policy over no preparation: {gain:+.1%}")
print(f"Mean irrigation coverage at year 30: {hI_opt[-1].mean():.1%}")
print(f"Mean reserve at year 30: {hS_opt[-1].mean():.1f} million food units")

# ---------------------------------------------------------------
# 5. One figure with every graph
# ---------------------------------------------------------------
fig = plt.figure(figsize=(20, 12))
gs = GridSpec(2, 3, figure=fig, hspace=0.3, wspace=0.25)
years = np.arange(T + 1)

# (1) rainfall distribution
ax1 = fig.add_subplot(gs[0, 0])
ax1.hist(R_test.ravel(), bins=60, color="steelblue", alpha=0.8)
ax1.axvline(DROUGHT, color="red", ls="--", lw=2, label=f"Drought threshold = {DROUGHT}")
ax1.set_title("(1) Annual rainfall index")
ax1.set_xlabel("Rainfall index (mean = 1)")
ax1.set_ylabel("Frequency")
ax1.legend()

# (2) 3D surface
ax2 = fig.add_subplot(gs[0, 1], projection="3d")
Sg, Ag = np.meshgrid(s_grid, a_grid)
ax2.plot_surface(Sg, Ag, Z, cmap="viridis", edgecolor="none", alpha=0.9)
ax2.scatter([s_best], [a_best], [Z[i_best, j_best]], color="red", s=80)
ax2.set_title("(2) Expected population at year 30")
ax2.set_xlabel("Target reserve s")
ax2.set_ylabel("Irrigation share a")
ax2.set_zlabel("Population (millions)")
ax2.view_init(elev=28, azim=-125)

# (3) heat map
ax3 = fig.add_subplot(gs[0, 2])
im = ax3.imshow(Z, origin="lower", aspect="auto", cmap="viridis",
                extent=[s_grid[0], s_grid[-1], a_grid[0], a_grid[-1]])
ax3.scatter([s_best], [a_best], color="red", marker="*", s=250,
            label=f"Optimum (s={s_best:.2f}, a={a_best:.2f})")
ax3.set_title("(3) Heat map of the same surface")
ax3.set_xlabel("Target reserve s")
ax3.set_ylabel("Irrigation share a")
ax3.legend(loc="lower right")
fig.colorbar(im, ax=ax3, label="Expected population (millions)")

# (4) population paths
ax4 = fig.add_subplot(gs[1, 0])
for hP, color, label in [(hP_base, "gray", "No preparation"),
                         (hP_opt, "crimson", "Optimal policy")]:
    lo, med, hi = np.percentile(hP, [5, 50, 95], axis=1)
    ax4.fill_between(years, lo, hi, color=color, alpha=0.2)
    ax4.plot(years, med, color=color, lw=2.5, label=f"{label} (median, 5-95% band)")
    ax4.plot(years, hP.mean(axis=1), color=color, lw=1.5, ls="--")
ax4.set_title("(4) Population paths (dashed = mean)")
ax4.set_xlabel("Year")
ax4.set_ylabel("Population (millions)")
ax4.legend()
ax4.grid(alpha=0.3)

# (5) final distribution
ax5 = fig.add_subplot(gs[1, 1])
bins = np.linspace(0, max(P_opt.max(), P_base.max()), 60)
ax5.hist(P_base, bins=bins, color="gray", alpha=0.6, label=f"No preparation (mean {P_base.mean():.0f})")
ax5.hist(P_opt, bins=bins, color="crimson", alpha=0.6, label=f"Optimal policy (mean {P_opt.mean():.0f})")
ax5.set_title("(5) Distribution of final population")
ax5.set_xlabel("Population at year 30 (millions)")
ax5.set_ylabel("Number of scenarios")
ax5.legend()

# (6) irrigation and reserve under the optimal policy
ax6 = fig.add_subplot(gs[1, 2])
ax6.plot(years, hI_opt.mean(axis=1) * 100, color="royalblue", lw=2.5, label="Irrigation coverage (%)")
ax6.set_xlabel("Year")
ax6.set_ylabel("Irrigation coverage (%)", color="royalblue")
ax6.grid(alpha=0.3)
ax6b = ax6.twinx()
ax6b.plot(years, hS_opt.mean(axis=1), color="darkorange", lw=2.5, label="Reserve (million food units)")
ax6b.set_ylabel("Reserve (million food units)", color="darkorange")
ax6.set_title("(6) Mean irrigation and reserve (optimal policy)")
h1, l1 = ax6.get_legend_handles_labels()
h2, l2 = ax6b.get_legend_handles_labels()
ax6.legend(h1 + h2, l1 + l2, loc="center right")

plt.show()

3. Code Walkthrough

Section 1: Parameters

Symbol	Variable	Value	Meaning
$T$	`T`	30	planning horizon in years
$P_0$	`P0`	100	initial population (millions)
$F_{\max}$	`F_MAX`	130	food output under perfect rainfall, which acts as a carrying capacity
$g$	`G`	0.025	maximum annual growth rate
$m$	`M`	0.5	how strongly shortages shrink the population
$\delta$	`SPOIL`	0.05	annual spoilage of stored grain
$C$	`COST_FULL`	350	cost of irrigating the entire country
$\varepsilon$	`EFF`	0.5	share of the rainfall shortfall that full irrigation can offset
$u$	`UPKEEP`	2.5	annual upkeep per unit of irrigation coverage
CV	`CV`	0.45	year-to-year variability of rainfall

make_rain draws a matrix of shape (scenarios, years) from a Gamma distribution. Rainfall cannot be negative, and the Gamma distribution has a long tail on the wet side and a fat mass on the dry side. With the shape parameter $k=1/\mathrm{CV}^2$ and scale $1/k$, the mean is exactly 1.

Section 2: The Simulator

The key design choice is vectorization. P, I and S are arrays with one element per scenario, so a single pass over the 30 years advances thousands of scenarios at once through NumPy operations. A plain Python loop over scenarios would call the inner year-by-year logic millions of times across the policy grid. With vectorization, the whole grid search finishes in a couple of seconds, so no further speed-up is needed.

Inside the yearly loop, the order of events matters:

Harvest. rr is the capped rainfall index $\min(1,R_t)$. The harvest is rr + (1 - rr) * EFF * I, which is exactly $y_t$ above, and prod subtracts the irrigation upkeep.
Surplus or deficit. Exactly one of surplus and deficit is positive in any year.
Draw the reserve. In a deficit year, draw moves grain from the reserve to cover the gap.
Invest in irrigation. A share a of the surplus goes to inv, capped so coverage never exceeds 100%.
Store grain. room is the gap between the target reserve $sP_t$ and the current reserve, and store fills it from the remaining surplus.
Population update. Anything still unmet after the reserve is the shortage, which drives a decline. Otherwise, growth is proportional to the leftover food, saturating at $g$ once the leftover reaches 10% of the population.
Spoilage. The reserve shrinks by the factor $1-\delta$.

When track=True, the function also records the full history of population, irrigation and reserves, which we need for the time-series plots.

Section 3: Grid Search

We evaluate every combination of 21 irrigation shares $a$ and 17 reserve targets $s$, which is 357 policies. Each is scored by the mean final population over the 3,000 training scenarios and stored in the matrix Z. np.argmax followed by np.unravel_index recovers the best pair of grid indices.

Section 4: Out-of-Sample Evaluation

Picking the winner on the same scenarios used to score it would flatter the policy. So the chosen policy is re-run on 5,000 new scenarios generated with a different seed. We compare four cases: no preparation, irrigation only, stockpile only, and the optimal combination. For each case we report the mean, the 5th percentile (a bad-luck measure), the median, and the share of scenarios in which the population falls below 50 million, which we call a collapse.

Section 5: Graphs

All six panels are drawn into a single figure through GridSpec, so one plt.show() produces everything.

4. Results

Console Output

=== Settings ===
Horizon: 30 years, scenarios (search/test): 3000/5000
Drought probability per year (rain < 0.7): 0.275
Grid search: 357 policies in 1.90 s

=== Optimal policy ===
Share of surplus invested in irrigation a* = 0.80
Target reserve s* = 0.250 years of consumption

=== Population after 30 years (millions, test scenarios) ===
Policy                           Mean   5th pct   Median  Collapse<50
No preparation (a=0, s=0)        72.5      50.7     72.0        4.4%
Irrigation only                 103.4      91.2    103.3        0.0%
Stockpile only                   86.4      63.2     86.3        0.3%
Optimal combination             108.5      97.2    108.8        0.0%

Gain of optimal policy over no preparation: +49.6%
Mean irrigation coverage at year 30: 100.0%
Mean reserve at year 30: 20.3 million food units

Graphs

5. Reading the Results

The comparison table

In a run of this model, the drought probability per year (rainfall index below 0.7) came out at roughly 27%, which is more than one year in four. Under those conditions the picture was clear:

No preparation left the population at about 72 million on average, down from the initial 100 million, and roughly 4% of scenarios collapsed below 50 million.
Stockpile only raised the average to about 86 million. Reserves smooth out single bad years but cannot fix a structural problem: when droughts come back to back, the reserve runs dry and spoilage eats into it the whole time.
Irrigation only did better, at about 103 million, because it permanently raises the harvest in every dry year rather than only compensating after the fact.
The optimal combination reached about 108 million on average, roughly a 50% improvement over doing nothing, and it also had the best 5th percentile. Preparation does not just raise the average, it lifts the floor.

The optimum found by the search was a high irrigation share (about 80% of surplus) together with a modest reserve of a quarter of a year of consumption. Irrigation does the heavy lifting, and a small reserve covers the gaps while irrigation is still being built up and the shocks that irrigation cannot absorb.

Graph (1): Rainfall distribution

The histogram shows the right-skewed rainfall index with the red drought threshold at 0.7. The large mass to the left of the line is the reason this problem is hard: bad years are common, not rare events.

Graph (2): 3D surface of expected population

The horizontal axes are the reserve target $s$ and the irrigation share $a$, and the vertical axis is the expected population at year 30. Moving along the $a$ axis produces a steep rise, which shows how valuable irrigation is. Moving along the $s$ axis gives a much gentler hill that falls off at large $s$, because holding too much grain costs spoilage and growth. The red point marks the optimum.

Graph (3): Heat map

The same surface seen from above. The bright region is a broad plateau along high irrigation shares with small reserves. The plateau is wide, which means the policy is forgiving: you do not need to hit the exact optimum to capture most of the benefit.

Graph (4): Population paths

The median and the 5th-to-95th percentile band for both policies. Without preparation, the band widens over time and its lower edge drops, showing how one unlucky run of droughts can drive a population down. Under the optimal policy the whole band sits higher and is narrower, so the outcome is both better and more predictable.

Graph (5): Distribution of final population

The gray histogram (no preparation) is wide and has a long left tail, which is the collapse risk. The red histogram (optimal policy) is shifted right and concentrated. Shrinking that left tail is the real prize of drought preparation.

Graph (6): Irrigation coverage and reserve over time

Irrigation coverage climbs steadily and approaches 100% as surplus is converted into infrastructure, while the average reserve rises early and then levels off at its target. This confirms the strategy: build irrigation first, and keep a modest buffer for the shocks that remain.

6. Summary

Planning for the average year is not enough when rainfall varies strongly, because the population is hurt more by bad stretches than it gains from good ones.
Irrigation investment changes the structure of the system, while grain reserves only buffer it. In this example irrigation was the more valuable tool.
The best policy combined both, and its advantage showed up most strongly in the worst-case outcomes.
Vectorizing the simulation over scenarios makes it practical to search hundreds of policies in seconds, and testing the winner on fresh scenarios keeps the comparison honest.

Natural extensions include making the irrigation share depend on the current reserve, adding a cost for the reserve itself, or replacing the grid search with a dynamic-programming or reinforcement-learning approach to find a policy that reacts to the state of the system each year.

Maximizing Population Capacity Under a Budget

November 2, 2026

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

import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

np.random.seed(0)

# ------------------------------------------------------------
# 1. Problem data
# ------------------------------------------------------------
cities = ["Aoba", "Bell", "Cedar", "Dune", "Elm", "Fjord", "Glen", "Haven", "Iris", "Jade"]
coords = np.array([[1, 8], [4, 9], [8, 8], [2, 5], [5, 5],
                   [9, 5], [1, 2], [4, 1.5], [7, 2], [9.5, 1]], dtype=float)
population = np.array([120, 85, 60, 95, 150, 70, 55, 110, 80, 45], dtype=float)  # thousands

# (name, type, cost in $M, served city indices)
segments = [
    ("R1", "Road",  30, [0, 1, 3]),
    ("R2", "Road",  25, [1, 2, 4]),
    ("R3", "Road",  35, [3, 4, 6, 7]),
    ("R4", "Road",  28, [4, 5, 8]),
    ("R5", "Road",  22, [7, 8, 9]),
    ("G1", "Grid",  40, [0, 1, 3, 4]),
    ("G2", "Grid",  32, [2, 5, 4]),
    ("G3", "Grid",  38, [3, 6, 7, 4]),
    ("G4", "Grid",  26, [7, 8, 9]),
    ("G5", "Grid",  18, [1, 2]),
    ("W1", "Water", 20, [0, 3, 6]),
    ("W2", "Water", 24, [4, 7, 8]),
    ("W3", "Water", 16, [2, 5]),
    ("W4", "Water", 21, [1, 4, 8, 9]),
]
types = ["Road", "Grid", "Water"]
weights = {"Road": 0.30, "Grid": 0.40, "Water": 0.30}
BUDGET = 120

n_seg, n_city = len(segments), len(cities)
seg_names = [s[0] for s in segments]
seg_type = [s[1] for s in segments]
seg_cost = np.array([s[2] for s in segments], dtype=float)

# coverage tensor: A[t][j, i] = 1 if segment j (of type t) serves city i
A = {t: np.zeros((n_seg, n_city)) for t in types}
for j, (nm, tp, c, srv) in enumerate(segments):
    A[tp][j, srv] = 1.0

# ------------------------------------------------------------
# 2. Exhaustive search (vectorized over all 2^14 subsets)
# ------------------------------------------------------------
masks = np.arange(2 ** n_seg)
bits = ((masks[:, None] >> np.arange(n_seg)) & 1).astype(float)   # (16384, 14)
cost_all = bits @ seg_cost

pop_by_type = np.zeros((len(masks), len(types)))   # covered population per type
cover_all = {}
for k, t in enumerate(types):
    cov = (bits @ A[t]) > 0                         # (16384, 10) boolean
    cover_all[t] = cov
    pop_by_type[:, k] = cov @ population

w_vec = np.array([weights[t] for t in types])
capacity_all = pop_by_type @ w_vec

feasible = cost_all <= BUDGET
best_idx = np.argmax(np.where(feasible, capacity_all, -1.0))
best_cap = capacity_all[best_idx]
best_sel = [j for j in range(n_seg) if (best_idx >> j) & 1]

# ------------------------------------------------------------
# 3. Greedy baseline (marginal gain per cost)
# ------------------------------------------------------------
def capacity_of(sel):
    cov = {t: np.zeros(n_city, dtype=bool) for t in types}
    for j in sel:
        cov[seg_type[j]] |= A[seg_type[j]][j] > 0
    return sum(weights[t] * population[cov[t]].sum() for t in types)

def greedy(budget):
    sel, spent = [], 0.0
    while True:
        base = capacity_of(sel)
        best_j, best_ratio = None, 0.0
        for j in range(n_seg):
            if j in sel or spent + seg_cost[j] > budget:
                continue
            ratio = (capacity_of(sel + [j]) - base) / seg_cost[j]
            if ratio > best_ratio:
                best_j, best_ratio = j, ratio
        if best_j is None:
            break
        sel.append(best_j)
        spent += seg_cost[best_j]
    return sel

greedy_sel = greedy(BUDGET)
greedy_cap = capacity_of(greedy_sel)

# ------------------------------------------------------------
# 4. Console report
# ------------------------------------------------------------
print("=" * 60)
print(f"Budget: {BUDGET} $M | Candidates: {n_seg} | Subsets searched: {2**n_seg:,}")
print("=" * 60)
print("\n[Optimal selection]")
print("Segments :", ", ".join(seg_names[j] for j in best_sel))
print(f"Cost     : {cost_all[best_idx]:.0f} $M")
print(f"Capacity : {best_cap:.2f} (thousand-weighted persons)")
print("\n[Greedy selection]")
print("Segments :", ", ".join(seg_names[j] for j in greedy_sel))
print(f"Cost     : {seg_cost[greedy_sel].sum():.0f} $M")
print(f"Capacity : {greedy_cap:.2f}")
print(f"\nGreedy / Optimal = {greedy_cap / best_cap * 100:.2f} %")

cov_best = {t: cover_all[t][best_idx] for t in types}
table = pd.DataFrame(
    {"Population(k)": population.astype(int),
     **{t: cov_best[t].astype(int) for t in types}},
    index=cities)
table["Weighted capacity"] = [
    population[i] * sum(weights[t] * cov_best[t][i] for t in types) for i in range(n_city)]
print("\n[City-level coverage under the optimal plan]")
print(table.round(1).to_string())

# ------------------------------------------------------------
# 5. Sensitivity data for graphs
# ------------------------------------------------------------
budgets = np.linspace(0, cost_all.max(), 41)
opt_curve = np.array([capacity_all[cost_all <= b].max() for b in budgets])
gre_curve = np.array([capacity_of(greedy(b)) for b in budgets])

B_grid = np.linspace(40, 250, 22)
G_grid = np.linspace(0.1, 0.8, 22)
Z = np.zeros((len(G_grid), len(B_grid)))
for a, wg in enumerate(G_grid):
    wr = (1 - wg) / 2
    cap = pop_by_type @ np.array([wr, wg, wr])
    for b, B in enumerate(B_grid):
        Z[a, b] = cap[cost_all <= B].max()
Bm, Gm = np.meshgrid(B_grid, G_grid)

# ------------------------------------------------------------
# 6. One combined figure
# ------------------------------------------------------------
type_color = {"Road": "#e07b39", "Grid": "#d6c21f", "Water": "#2a7fbf"}
fig = plt.figure(figsize=(20, 15))
gs = GridSpec(2, 2, figure=fig, hspace=0.25, wspace=0.18)

# (a) map
ax1 = fig.add_subplot(gs[0, 0])
offsets = {"Road": -0.12, "Grid": 0.0, "Water": 0.12}
for j in range(n_seg):
    sel = j in best_sel
    pts = coords[segments[j][3]]
    center = pts.mean(axis=0)
    for p in pts:
        ax1.plot([center[0], p[0]], [center[1] + offsets[seg_type[j]], p[1]],
                 color=type_color[seg_type[j]],
                 lw=3.0 if sel else 0.8,
                 alpha=0.95 if sel else 0.25,
                 ls="-" if sel else "--", zorder=1)
    if sel:
        ax1.text(center[0], center[1] + offsets[seg_type[j]], seg_names[j],
                 fontsize=9, ha="center", va="center", fontweight="bold",
                 bbox=dict(boxstyle="round,pad=0.2", fc="white",
                           ec=type_color[seg_type[j]]), zorder=3)
ax1.scatter(coords[:, 0], coords[:, 1], s=population * 6, c="#444444",
            alpha=0.85, zorder=2, edgecolor="white", linewidth=1.5)
for i, nm in enumerate(cities):
    ax1.annotate(f"{nm}\n{int(population[i])}k", coords[i], textcoords="offset points",
                 xytext=(0, 18), ha="center", fontsize=9)
for t in types:
    ax1.plot([], [], color=type_color[t], lw=3, label=t)
ax1.set_title("(a) Optimal Infrastructure Network (solid = built, dashed = skipped)", fontsize=13)
ax1.set_xlim(-0.5, 11); ax1.set_ylim(0, 10.5)
ax1.legend(loc="lower left"); ax1.set_xlabel("x"); ax1.set_ylabel("y")
ax1.grid(alpha=0.2)

# (b) 3D bar: city x type
ax2 = fig.add_subplot(gs[0, 1], projection="3d")
for k, t in enumerate(types):
    xs = np.arange(n_city)
    zs = population * weights[t] * cov_best[t]
    ax2.bar3d(xs, np.full(n_city, k), np.zeros(n_city), 0.6, 0.6, zs,
              color=type_color[t], alpha=0.9, shade=True)
ax2.set_xticks(np.arange(n_city) + 0.3)
ax2.set_xticklabels(cities, rotation=40, fontsize=8)
ax2.set_yticks(np.arange(len(types)) + 0.3)
ax2.set_yticklabels(types, fontsize=9)
ax2.set_zlabel("Weighted capacity")
ax2.set_title("(b) 3D Capacity by City and Infrastructure Type", fontsize=13)
ax2.view_init(elev=28, azim=-60)

# (c) cost vs capacity
ax3 = fig.add_subplot(gs[1, 0])
ax3.plot(budgets, opt_curve, "-o", ms=4, color="#c0392b", label="Optimal (exhaustive)")
ax3.plot(budgets, gre_curve, "-s", ms=4, color="#2c7fb8", label="Greedy")
ax3.axvline(BUDGET, color="gray", ls="--")
ax3.scatter([BUDGET], [best_cap], s=140, color="#c0392b", zorder=5)
ax3.annotate(f"Budget {BUDGET}\nOptimal {best_cap:.1f}", (BUDGET, best_cap),
             textcoords="offset points", xytext=(15, -45), fontsize=10,
             arrowprops=dict(arrowstyle="->"))
ax3.set_xlabel("Budget ($M)"); ax3.set_ylabel("Capacity")
ax3.set_title("(c) Budget vs. Achievable Capacity", fontsize=13)
ax3.legend(); ax3.grid(alpha=0.3)

# (d) 3D surface
ax4 = fig.add_subplot(gs[1, 1], projection="3d")
surf = ax4.plot_surface(Bm, Gm, Z, cmap="viridis", edgecolor="none", alpha=0.95)
ax4.set_xlabel("Budget ($M)"); ax4.set_ylabel("Grid weight")
ax4.set_zlabel("Optimal capacity")
ax4.set_title("(d) Optimal Capacity Surface: Budget × Grid Weight", fontsize=13)
ax4.view_init(elev=26, azim=-125)
fig.colorbar(surf, ax=ax4, shrink=0.55, pad=0.08)

fig.suptitle("Infrastructure Prioritization under a Budget Constraint", fontsize=18, y=0.95)
plt.show()

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.
bits expands 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_cost gives the cost of every subset at once.
For each type, (bits @ A[t]) > 0 is 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 @ population converts coverage to covered population per type, and pop_by_type @ w_vec produces 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_curve is the best achievable capacity for 41 budget levels. Because the capacity of every subset is already computed, each point is just a masked max.
Z is 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$. Since pop_by_type is 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.

Optimal Resource Allocation Among Three Competing Regions

November 1, 2026

Maximizing Total Population with Python

Imagine a government that must divide a fixed development budget among three regions. Every unit of budget spent in one region improves living conditions there and lets its population grow, but that same unit is no longer available to the other two regions. How should the budget be split so that the total population ten years from now is as large as possible?

This is a classic constrained optimization problem. In this article we build a concrete model, solve it numerically with Python, verify the answer by brute-force search, and visualize the result in 3D.

1. Problem Setting

Three regions $i = 1, 2, 3$ (Region A, B, C) share a total budget $R$. Let $x_i \ge 0$ be the amount invested in region $i$.

$$
\max_{x_1,x_2,x_3}\ \sum_{i=1}^{3} P_i(T;,x_i)
\quad \text{subject to} \quad \sum_{i=1}^{3} x_i \le R,\quad x_i \ge 0
$$

Population model

Population in each region follows logistic growth:

$$
\frac{dP_i}{dt} = r_i P_i\left(1-\frac{P_i}{K_i}\right),\qquad P_i(0)=P_{i,0}
$$

The closed-form solution is

$$
P_i(t) = \frac{K_i P_{i,0}}{P_{i,0} + \left(K_i - P_{i,0}\right)e^{-r_i t}}
$$

Investment raises the carrying capacity of the region, but with diminishing returns:

$$
K_i(x_i) = K_{i,0} + \Delta K_i\left(1 - e^{-x_i/s_i}\right)
$$

Here $K_{i,0}$ is the capacity without investment, $\Delta K_i$ is the maximum additional capacity that investment can create, and $s_i$ is the scale of investment at which the benefit starts to saturate.

Optimality conditions

Because $P_i(T)$ is increasing and concave in $K_i$, and $K_i$ is increasing and concave in $x_i$, each $P_i(T;x_i)$ is an increasing concave function of $x_i$. The problem is therefore a concave maximization over a simplex, so any local optimum is the global optimum. Introducing a Lagrange multiplier $\lambda$ for the budget constraint gives the KKT conditions:

$$
\frac{\partial P_i(T;x_i)}{\partial x_i} = \lambda \quad (x_i > 0),\qquad
\frac{\partial P_i(T;x_i)}{\partial x_i} \le \lambda \quad (x_i = 0)
$$

In words: at the optimum, the marginal population gain per unit of budget is equal across all funded regions. The multiplier $\lambda$ is the shadow price of the budget, that is, the extra population obtained from one more unit of budget.

2. Example Problem

Parameter	Region A	Region B	Region C
Initial population $P_{i,0}$ (thousand)	50	30	20
Capacity without investment $K_{i,0}$ (thousand)	80	60	50
Maximum capacity gain $\Delta K_i$ (thousand)	120	100	90
Investment scale $s_i$	40	25	20
Growth rate $r_i$ (per year)	0.30	0.45	0.60

The total budget is $R = 100$ units and the planning horizon is $T = 10$ years.

3. Source Code

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# ------------------------------------------------------------
# Model parameters
# ------------------------------------------------------------
names = ["Region A", "Region B", "Region C"]
P0 = np.array([50.0, 30.0, 20.0])     # initial population (thousand)
K0 = np.array([80.0, 60.0, 50.0])     # capacity without investment
dK = np.array([120.0, 100.0, 90.0])   # maximum capacity gain from investment
s = np.array([40.0, 25.0, 20.0])      # investment scale (saturation)
r = np.array([0.30, 0.45, 0.60])      # growth rate per year
T = 10.0                              # planning horizon (years)
R = 100.0                             # total budget

region_colors = ["#4cc9f0", "#f72585", "#ffd166"]


# ------------------------------------------------------------
# Population model (vectorized)
# ------------------------------------------------------------
def final_population(x):
    """Population of each region at time T. x has shape (..., 3)."""
    x = np.asarray(x, dtype=float)
    K = K0 + dK * (1.0 - np.exp(-x / s))
    q = np.exp(-r * T)
    return K * P0 / (P0 * (1.0 - q) + K * q)


def total_population(x):
    return np.sum(final_population(x), axis=-1)


def population_region(xi, i):
    """Final population of a single region i for investment xi (array allowed)."""
    xi = np.asarray(xi, dtype=float)
    K = K0[i] + dK[i] * (1.0 - np.exp(-xi / s[i]))
    q = np.exp(-r[i] * T)
    return K * P0[i] / (P0[i] * (1.0 - q) + K * q)


def trajectory(x, t):
    """Population over time for allocation x. Returns array of shape (len(t), 3)."""
    x = np.asarray(x, dtype=float)
    K = K0 + dK * (1.0 - np.exp(-x / s))
    tt = np.asarray(t, dtype=float)[:, None]
    return K * P0 / (P0 + (K - P0) * np.exp(-r * tt))


# ------------------------------------------------------------
# Optimization with SLSQP (budget constraint: sum(x) = R)
# ------------------------------------------------------------
def objective(x):
    return -total_population(x)


constraints = [{"type": "eq", "fun": lambda x: R - np.sum(x)}]
bounds = [(0.0, R)] * 3
x_init = np.full(3, R / 3.0)

res = minimize(
    objective,
    x_init,
    method="SLSQP",
    bounds=bounds,
    constraints=constraints,
    options={"ftol": 1e-12, "maxiter": 500},
)
x_opt = res.x
total_opt = total_population(x_opt)

# ------------------------------------------------------------
# Brute-force verification on a simplex grid (step = 1 unit)
# ------------------------------------------------------------
g = np.arange(0, int(R) + 1)
G1, G2 = np.meshgrid(g, g, indexing="ij")
mask = (G1 + G2) <= R
x1 = G1[mask].astype(float)
x2 = G2[mask].astype(float)
x3 = R - x1 - x2
alloc_grid = np.stack([x1, x2, x3], axis=1)
total_grid = total_population(alloc_grid)
k_best = int(np.argmax(total_grid))
x_grid_best = alloc_grid[k_best]
total_grid_best = total_grid[k_best]

# ------------------------------------------------------------
# Comparison strategies
# ------------------------------------------------------------
strategies = {
    "Optimal": x_opt,
    "Equal split": np.full(3, R / 3.0),
    "Proportional to initial pop.": R * P0 / P0.sum(),
    "No investment": np.zeros(3),
}
strategy_totals = {k: float(total_population(v)) for k, v in strategies.items()}

# ------------------------------------------------------------
# Marginal returns (shadow price check)
# ------------------------------------------------------------
eps = 1e-6
marginal = np.zeros(3)
for i in range(3):
    xp = x_opt.copy()
    xm = x_opt.copy()
    xp[i] += eps
    xm[i] = max(xm[i] - eps, 0.0)
    marginal[i] = (total_population(xp) - total_population(xm)) / (xp[i] - xm[i])

# ------------------------------------------------------------
# Console output
# ------------------------------------------------------------
print("=" * 64)
print("Optimal resource allocation among three competing regions")
print("=" * 64)
print(f"Solver status          : {res.message}")
print(f"Total budget           : {R:.1f}")
print(f"Planning horizon       : {T:.0f} years")
print("-" * 64)
print(f"{'Region':<12}{'Investment':>14}{'Final pop.':>14}{'Marginal gain':>16}")
final_opt = final_population(x_opt)
for i in range(3):
    print(f"{names[i]:<12}{x_opt[i]:>14.3f}{final_opt[i]:>14.3f}{marginal[i]:>16.5f}")
print("-" * 64)
print(f"Total population (SLSQP)      : {total_opt:.4f} thousand")
print(f"Total population (grid search): {total_grid_best:.4f} thousand")
print(f"Grid-search allocation        : "
      f"A={x_grid_best[0]:.0f}, B={x_grid_best[1]:.0f}, C={x_grid_best[2]:.0f}")
print("-" * 64)
print(f"{'Strategy':<32}{'Total pop.':>14}{'vs Optimal':>14}")
for k, v in strategy_totals.items():
    ratio = 100.0 * v / strategy_totals["Optimal"]
    print(f"{k:<32}{v:>14.3f}{ratio:>13.2f}%")
print("=" * 64)

# ------------------------------------------------------------
# Visualization (all panels in a single figure)
# ------------------------------------------------------------
plt.style.use("dark_background")
fig = plt.figure(figsize=(21, 12))
gs = fig.add_gridspec(2, 3)
fig.suptitle("Resource Allocation Among Three Competing Regions",
             fontsize=20, fontweight="bold")

# (1) 3D surface of total population on the budget simplex
ax1 = fig.add_subplot(gs[0, 0], projection="3d")
surf = ax1.plot_trisurf(x1, x2, total_grid, cmap="viridis",
                        linewidth=0.0, antialiased=True, alpha=0.95)
ax1.scatter([x_opt[0]], [x_opt[1]], [total_opt], color="red", s=140,
            depthshade=False, label="Optimum")
ax1.set_xlabel("Investment in A")
ax1.set_ylabel("Investment in B")
ax1.set_zlabel("Total population (thousand)")
ax1.set_title("Total population surface\n(Investment in C = R - A - B)")
ax1.view_init(elev=28, azim=-125)
ax1.legend(loc="upper left")
fig.colorbar(surf, ax=ax1, shrink=0.6, pad=0.1)

# (2) Contour map on the simplex
ax2 = fig.add_subplot(gs[0, 1])
cf = ax2.tricontourf(x1, x2, total_grid, levels=30, cmap="viridis")
ax2.tricontour(x1, x2, total_grid, levels=10, colors="white",
               linewidths=0.4, alpha=0.6)
eq = strategies["Equal split"]
ax2.scatter([x_opt[0]], [x_opt[1]], color="red", s=140, marker="*",
            label="Optimal", zorder=5)
ax2.scatter([eq[0]], [eq[1]], color="white", s=80, marker="o",
            label="Equal split", zorder=5)
ax2.set_xlabel("Investment in A")
ax2.set_ylabel("Investment in B")
ax2.set_title("Contours of total population")
ax2.legend(loc="upper right")
fig.colorbar(cf, ax=ax2)

# (3) Allocation comparison
ax3 = fig.add_subplot(gs[0, 2])
labels3 = ["Optimal", "Equal split", "Proportional to initial pop."]
width = 0.25
pos = np.arange(3)
bar_colors = ["#06d6a0", "#8d99ae", "#ef476f"]
for j, lab in enumerate(labels3):
    vals = strategies[lab]
    bars = ax3.bar(pos + (j - 1) * width, vals, width,
                   label=lab, color=bar_colors[j])
    for b, v in zip(bars, vals):
        ax3.text(b.get_x() + b.get_width() / 2.0, v + 0.8, f"{v:.1f}",
                 ha="center", va="bottom", fontsize=9)
ax3.set_xticks(pos)
ax3.set_xticklabels(names)
ax3.set_ylabel("Investment")
ax3.set_title("Budget allocation by strategy")
ax3.legend()

# (4) Population trajectories
ax4 = fig.add_subplot(gs[1, 0])
t = np.linspace(0.0, T, 200)
traj_opt = trajectory(x_opt, t)
traj_none = trajectory(np.zeros(3), t)
for i in range(3):
    ax4.plot(t, traj_opt[:, i], color=region_colors[i], lw=2.5,
             label=f"{names[i]} (optimal)")
    ax4.plot(t, traj_none[:, i], color=region_colors[i], lw=1.5, ls="--",
             alpha=0.8, label=f"{names[i]} (no investment)")
ax4.set_xlabel("Years")
ax4.set_ylabel("Population (thousand)")
ax4.set_title("Population trajectories")
ax4.legend(fontsize=8, ncol=2)
ax4.grid(alpha=0.25)

# (5) Marginal returns and the shadow price
ax5 = fig.add_subplot(gs[1, 1])
xs = np.linspace(0.0, R, 400)
for i in range(3):
    m = np.gradient(population_region(xs, i), xs)
    ax5.plot(xs, m, color=region_colors[i], lw=2.5, label=names[i])
    m_at = float(np.interp(x_opt[i], xs, m))
    ax5.scatter([x_opt[i]], [m_at], color=region_colors[i], s=120,
                edgecolor="white", zorder=5)
ax5.axhline(float(np.max(marginal)), color="white", ls=":", lw=1.5,
            label="Shadow price level")
ax5.set_xlabel("Investment in a region")
ax5.set_ylabel("Marginal population gain per unit")
ax5.set_title("Marginal returns (dots = optimal allocation)")
ax5.legend()
ax5.grid(alpha=0.25)

# (6) Total population by strategy
ax6 = fig.add_subplot(gs[1, 2])
keys = list(strategy_totals.keys())
vals6 = [strategy_totals[k] for k in keys]
colors6 = ["#06d6a0", "#8d99ae", "#ef476f", "#6c757d"]
bars6 = ax6.bar(np.arange(len(keys)), vals6, color=colors6)
for b, v in zip(bars6, vals6):
    ax6.text(b.get_x() + b.get_width() / 2.0, v + 2.0, f"{v:.1f}",
             ha="center", va="bottom", fontsize=10)
ax6.set_xticks(np.arange(len(keys)))
ax6.set_xticklabels(["Optimal", "Equal\nsplit", "Proportional\nto initial pop.",
                     "No\ninvestment"])
ax6.set_ylabel("Total population after T years (thousand)")
ax6.set_title("Total population by strategy")
ax6.set_ylim(0, max(vals6) * 1.15)

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

4. Code Walkthrough

Parameters and the population model

The arrays P0, K0, dK, s, and r hold the table of parameters from Section 2, one entry per region. The function final_population implements

$$
P_i(T) = \frac{K_i P_{i,0}}{P_{i,0}\left(1-e^{-r_iT}\right) + K_i, e^{-r_iT}}
$$

which is the logistic solution evaluated at $t=T$. It is written with NumPy broadcasting, so it accepts either a single allocation of shape (3,) or a whole batch of shape (N, 3). This vectorization is what makes the brute-force search below essentially instantaneous: about five thousand candidate allocations are evaluated in a single call without any Python loop. trajectory evaluates the same formula at many time points to draw growth curves, and population_region evaluates one region at a time for the marginal-return plot.

Numerical optimization

The optimization uses SLSQP, a sequential quadratic programming method suited to smooth problems with equality constraints and bounds. Because SciPy minimizes, the objective is the negative of the total population. The constraint $R - \sum_i x_i = 0$ forces the whole budget to be spent, which loses nothing because population is increasing in every $x_i$. The starting point is an equal split. Since the problem is concave, SLSQP converges to the global optimum regardless of the starting point.

Brute-force verification

To confirm that the solver has not been trapped somewhere unexpected, the code enumerates every integer allocation $(x_1, x_2, x_3)$ with $x_1 + x_2 \le R$ and $x_3 = R - x_1 - x_2$. The mask G1 + G2 <= R keeps only points on the budget simplex, and np.argmax picks the best one. The grid optimum should agree with the SLSQP optimum to within the one-unit grid spacing, and the SLSQP total should be equal to or slightly higher than the grid total.

Marginal returns

After the optimum is found, the code estimates $\partial P_{\text{total}}/\partial x_i$ by central-style finite differences. If the theory is correct, every region that receives funding shows the same marginal gain, and that common value is the shadow price $\lambda$ of the budget.

Comparison strategies

Three reference policies are evaluated alongside the optimal plan: splitting equally, splitting in proportion to the initial population (a natural “fairness” rule), and investing nothing. They show how much the optimization actually gains over intuitive rules of thumb.

5. Results

Result images

Console output

================================================================
Optimal resource allocation among three competing regions
================================================================
Solver status          : Optimization terminated successfully
Total budget           : 100.0
Planning horizon       : 10 years
----------------------------------------------------------------
Region          Investment    Final pop.   Marginal gain
Region A            35.986       137.355         0.95683
Region B            33.593       128.951         0.95683
Region C            30.420       118.858         0.95683
----------------------------------------------------------------
Total population (SLSQP)      : 385.1636 thousand
Total population (grid search): 385.1561 thousand
Grid-search allocation        : A=36, B=34, C=30
----------------------------------------------------------------
Strategy                            Total pop.    vs Optimal
Optimal                                385.164       100.00%
Equal split                            384.874        99.92%
Proportional to initial pop.           379.521        98.53%
No investment                          186.835        48.51%
================================================================

6. Interpreting the Graphs

Top left: 3D surface of total population. The horizontal axes are the investments in Regions A and B, and the investment in Region C is whatever remains of the budget. The surface is a single smooth dome with no competing peaks, which is the visual signature of the concavity proved earlier. The red marker sits at the top of the dome, the best allocation the budget allows.

Top center: contour map. The same surface seen from above. The triangular domain is the set of all feasible budget splits. The star marks the optimal plan and the white dot marks the equal split. The distance between them shows how far an intuitive “split it evenly” policy sits from the true optimum, and the contour lines show how quickly the total falls off as the allocation moves away from the best point.

Top right: allocation by strategy. Side-by-side bars compare the optimal split with the equal split and the split proportional to initial population. Differences in bar height show how the optimal plan deliberately favors regions whose investment converts into population most efficiently, rather than regions that are simply large today.

Bottom left: population trajectories. Solid curves show each region’s population over the ten years under the optimal plan, and dashed curves show what would happen with no investment. The gap between solid and dashed curves of the same color is the benefit of funding that region. Regions with higher growth rates $r_i$ approach their capacity sooner, while slower-growing regions are still climbing at $t = T$.

Bottom center: marginal returns. Each curve shows the extra population gained from one more unit of budget as a function of the amount already invested in that region. All curves decline because of diminishing returns. The dots mark the optimal allocation, and for every funded region they lie at the same height, the dotted horizontal line representing the shadow price $\lambda$. This is the KKT condition made visible: moving a unit of budget from one region to another cannot help, because the gain in one is exactly cancelled by the loss in the other.

Bottom right: total population by strategy. The optimal plan produces the largest population, with the equal split and the proportional rule below it and the no-investment baseline far behind. The gap between the optimal plan and the intuitive rules is the practical value of solving the allocation problem properly.

7. Conclusion

A competing allocation problem looks complicated because every investment decision affects every other region, yet its structure is simple. With diminishing returns in each region, the total population is a concave function of the allocation, so a single numerical solve finds the global optimum, and the optimal plan is characterized by one clean principle: equalize the marginal gain per unit of budget across all funded regions. The shadow price $\lambda$ then tells decision makers exactly how much additional population one extra unit of budget would buy, which is valuable information when negotiating the size of the budget itself.

The same framework extends naturally to more regions, additional constraints such as minimum guaranteed funding per region, or investment that is spread over time as a dynamic optimization problem.

Optimizing Population Migration Between Cities

October 31, 2026

Maximizing the 10-Year Total with Python

The Problem

A country has several cities with very different characteristics. Industrial cities offer plenty of jobs but suffer from a housing shortage. Agricultural cities produce abundant food but offer few jobs. We can relocate residents between cities, and we must decide how many people to move where so that the total population after 10 years is as large as possible.

A Concrete Example

We model three cities (unit: thousand people).

City	Initial population	Housing (year 0)	Housing added per year	Jobs (year 0)	Jobs added per year	Food production (people fed)
Industrial	500	520	30	900	30	400
Agricultural	300	500	4	200	3	900
Commercial	200	330	8	380	10	300

Food is a national resource shared by everyone and grows by 1% per year.

The Mathematical Model

For city $i$ in year $t$, with population $P_{i,t}$, housing capacity $H_{i,t}$, and job capacity $J_{i,t}$, we define three satisfaction factors:

$$
h_{i,t}=\min!\left(1,\frac{H_{i,t}}{P_{i,t}}\right),\qquad
j_{i,t}=\min!\left(1,\frac{J_{i,t}}{P_{i,t}}\right),\qquad
f_{t}=\min!\left(1,\frac{\sum_i F_{i,t}}{\sum_i P_{i,t}}\right)
$$

The overall quality of life is

$$
q_{i,t}=h_{i,t},\bigl(0.6+0.4,j_{i,t}\bigr),f_t
$$

and the yearly population update is

$$
P_{i,t+1}=P_{i,t}\Bigl(1+b,q_{i,t}-d\Bigr),\qquad b=0.04,; d=0.01
$$

Let $m_{ij}$ be the number of people moving from city $i$ to city $j$ at the start. A share $\lambda=0.04$ of the arriving migrants is lost while settling, so the starting population becomes

$$
P_{j,0}’=P_{j,0}-\sum_{k\neq j} m_{jk}+(1-\lambda)\sum_{i\neq j} m_{ij}
$$

The optimization problem is

$$
\max_{m_{ij}\ge 0}\ \sum_{i=1}^{3}P_{i,10}\qquad \text{subject to}\qquad \sum_{j\neq i}m_{ij}\le 0.6,P_{i,0}
$$

To make the search unconstrained, each city gets a leaving ratio and a split ratio between its two destinations, giving six decision variables. Because the objective is non-convex and full of $\min$ kinks, we use differential evolution, a global optimizer.

Source Code

The simulation is fully vectorized with NumPy, so the 3,721 scenarios for the 3D surface are evaluated in a single call and the whole script finishes in a few seconds.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.optimize import differential_evolution

np.set_printoptions(precision=1, suppress=True)

# ------------------------------------------------------------
# 1. Problem settings (population unit: thousand people)
# ------------------------------------------------------------
CITIES = ["Industrial", "Agricultural", "Commercial"]
N_CITY = 3
YEARS = 10

P0 = np.array([500.0, 300.0, 200.0])      # initial population
H0 = np.array([520.0, 500.0, 330.0])      # housing capacity at year 0
dH = np.array([30.0, 4.0, 8.0])           # housing added per year
J0 = np.array([900.0, 200.0, 380.0])      # job capacity at year 0
dJ = np.array([30.0, 3.0, 10.0])          # jobs added per year
F0 = np.array([400.0, 900.0, 300.0])      # food production (people fed)
F_GROWTH = 0.01                           # annual food growth

BIRTH = 0.04
DEATH = 0.01
MOVE_LOSS = 0.04                          # share of migrants lost while settling

# ------------------------------------------------------------
# 2. Simulation (vectorized: pop_start may have shape (N, 3))
# ------------------------------------------------------------
def simulate(pop_start, years=YEARS, record=False):
    P = np.array(pop_start, dtype=float)
    history = [P.copy()]
    for t in range(years):
        H = H0 + dH * t
        J = J0 + dJ * t
        F_total = F0.sum() * (1.0 + F_GROWTH) ** t
        food = np.minimum(1.0, F_total / P.sum(axis=-1, keepdims=True))
        housing = np.minimum(1.0, H / P)
        jobs = np.minimum(1.0, J / P)
        quality = housing * (0.6 + 0.4 * jobs) * food
        growth = BIRTH * quality - DEATH
        P = P * (1.0 + growth)
        history.append(P.copy())
    if record:
        return np.array(history)
    return P

# ------------------------------------------------------------
# 3. Migration model
#    x = [out_I, out_A, out_C, split_I, split_A, split_C]
#    out_i   : share of city i that leaves
#    split_i : share of leavers sent to the first destination
# ------------------------------------------------------------
DEST = {0: (1, 2), 1: (0, 2), 2: (0, 1)}

def flow_matrix(x):
    x = np.asarray(x, dtype=float)
    M = np.zeros((N_CITY, N_CITY))
    for i in range(N_CITY):
        leavers = P0[i] * x[i]
        d1, d2 = DEST[i]
        M[i, d1] = leavers * x[3 + i]
        M[i, d2] = leavers * (1.0 - x[3 + i])
    return M

def after_migration(M):
    out = M.sum(axis=1)
    inflow = M.sum(axis=0)
    return P0 - out + inflow * (1.0 - MOVE_LOSS)

def total_after_10y(x):
    return simulate(after_migration(flow_matrix(x))).sum()

def objective(x):
    return -total_after_10y(x)

# ------------------------------------------------------------
# 4. Optimization
# ------------------------------------------------------------
bounds = [(0.0, 0.6)] * 3 + [(0.0, 1.0)] * 3
result = differential_evolution(
    objective, bounds, seed=42, popsize=30, maxiter=300, tol=1e-10, polish=True
)
x_opt = result.x
M_opt = flow_matrix(x_opt)
P_start_opt = after_migration(M_opt)

hist_base = simulate(P0, record=True)
hist_opt = simulate(P_start_opt, record=True)

# Simple rule-of-thumb strategy for comparison: move 20% of the agricultural city to the industrial city
x_rule = np.array([0.0, 0.20, 0.0, 1.0, 1.0, 1.0])
M_rule = flow_matrix(x_rule)
hist_rule = simulate(after_migration(M_rule), record=True)

# ------------------------------------------------------------
# 5. Console report
# ------------------------------------------------------------
print("=" * 64)
print(" Population Migration and Urban Growth (unit: thousand people)")
print("=" * 64)
print(f"{'City':<14}{'Initial':>10}{'No move':>12}{'Rule':>12}{'Optimal':>12}")
for i, name in enumerate(CITIES):
    print(f"{name:<14}{P0[i]:>10.1f}{hist_base[-1, i]:>12.1f}"
          f"{hist_rule[-1, i]:>12.1f}{hist_opt[-1, i]:>12.1f}")
print("-" * 64)
print(f"{'Total':<14}{P0.sum():>10.1f}{hist_base[-1].sum():>12.1f}"
      f"{hist_rule[-1].sum():>12.1f}{hist_opt[-1].sum():>12.1f}")
print()
print("Optimal migration flows (rows: from, columns: to, thousand people)")
print(f"{'':<14}" + "".join(f"{c:>14}" for c in CITIES))
for i, name in enumerate(CITIES):
    print(f"{name:<14}" + "".join(f"{M_opt[i, j]:>14.1f}" for j in range(N_CITY)))
print()
gain_base = hist_opt[-1].sum() - hist_base[-1].sum()
print(f"Gain over no migration : {gain_base:.1f} thousand "
      f"({gain_base / hist_base[-1].sum() * 100:.2f} %)")
gain_rule = hist_opt[-1].sum() - hist_rule[-1].sum()
print(f"Gain over rule of thumb: {gain_rule:.1f} thousand "
      f"({gain_rule / hist_rule[-1].sum() * 100:.2f} %)")

# ------------------------------------------------------------
# 6. Surface: Agricultural -> Industrial (a) vs Agricultural -> Commercial (b)
#    All other flows are fixed at the optimal values (fully vectorized)
# ------------------------------------------------------------
n_grid = 61
a_vals = np.linspace(0.0, P0[1], n_grid)
b_vals = np.linspace(0.0, P0[1], n_grid)
A_grid, B_grid = np.meshgrid(a_vals, b_vals)
valid = (A_grid + B_grid) <= P0[1]

out_I = M_opt[0].sum()
out_C = M_opt[2].sum()
inflow_I = M_opt[2, 0] + A_grid
inflow_A = M_opt[0, 1] + M_opt[2, 1]
inflow_C = M_opt[0, 2] + B_grid
start = np.stack([
    P0[0] - out_I + inflow_I * (1 - MOVE_LOSS),
    P0[1] - A_grid - B_grid + inflow_A * (1 - MOVE_LOSS),
    P0[2] - out_C + inflow_C * (1 - MOVE_LOSS),
], axis=-1).reshape(-1, N_CITY)
Z = simulate(start.copy()).sum(axis=1).reshape(A_grid.shape)
Z = np.where(valid, Z, np.nan)

best_idx = np.nanargmax(Z)
bi, bj = np.unravel_index(best_idx, Z.shape)
print()
print(f"Grid best (A->I, A->C): ({A_grid[bi, bj]:.1f}, {B_grid[bi, bj]:.1f}) "
      f"-> total {Z[bi, bj]:.1f} thousand")

# ------------------------------------------------------------
# 7. One combined figure
# ------------------------------------------------------------
colors = ["#d62728", "#2ca02c", "#1f77b4"]
years = np.arange(YEARS + 1)

fig = plt.figure(figsize=(20, 12))
fig.suptitle("Population Migration and Urban Growth: Maximizing the 10-Year Total",
             fontsize=18, fontweight="bold")

# (1) 3D surface
ax1 = fig.add_subplot(2, 3, 1, projection="3d")
surf = ax1.plot_surface(A_grid, B_grid, Z, cmap=cm.viridis, edgecolor="none", alpha=0.92)
ax1.scatter([A_grid[bi, bj]], [B_grid[bi, bj]], [Z[bi, bj]], color="red", s=80,
            depthshade=False, label="Best")
ax1.set_xlabel("Agri → Industrial (k)")
ax1.set_ylabel("Agri → Commercial (k)")
ax1.set_zlabel("Total pop. in 10 yrs (k)")
ax1.set_title("3D: Total Population vs Migration")
ax1.view_init(elev=28, azim=-125)
ax1.legend(loc="upper left")
fig.colorbar(surf, ax=ax1, shrink=0.55, pad=0.1)

# (2) Contour
ax2 = fig.add_subplot(2, 3, 2)
cf = ax2.contourf(A_grid, B_grid, Z, levels=25, cmap="viridis")
ax2.scatter([A_grid[bi, bj]], [B_grid[bi, bj]], color="red", s=90, marker="*",
            label="Best", zorder=5)
ax2.set_xlabel("Agri → Industrial (k)")
ax2.set_ylabel("Agri → Commercial (k)")
ax2.set_title("Contour Map of the Same Surface")
ax2.legend()
fig.colorbar(cf, ax=ax2)

# (3) Flow heatmap
ax3 = fig.add_subplot(2, 3, 3)
im = ax3.imshow(M_opt, cmap="Oranges")
ax3.set_xticks(range(N_CITY))
ax3.set_yticks(range(N_CITY))
ax3.set_xticklabels(CITIES)
ax3.set_yticklabels(CITIES)
ax3.set_xlabel("To")
ax3.set_ylabel("From")
ax3.set_title("Optimal Migration Flows (k)")
for i in range(N_CITY):
    for j in range(N_CITY):
        ax3.text(j, i, f"{M_opt[i, j]:.1f}", ha="center", va="center", fontsize=12)
fig.colorbar(im, ax=ax3)

# (4) City populations: no migration vs optimal
ax4 = fig.add_subplot(2, 3, 4)
for i, name in enumerate(CITIES):
    ax4.plot(years, hist_base[:, i], "--", color=colors[i], label=f"{name} (no move)")
    ax4.plot(years, hist_opt[:, i], "-", color=colors[i], linewidth=2.5,
             label=f"{name} (optimal)")
ax4.set_xlabel("Year")
ax4.set_ylabel("Population (k)")
ax4.set_title("City Populations Over 10 Years")
ax4.grid(alpha=0.3)
ax4.legend(fontsize=8, ncol=2)

# (5) Total population by strategy
ax5 = fig.add_subplot(2, 3, 5)
ax5.plot(years, hist_base.sum(axis=1), "--", color="gray", linewidth=2, label="No migration")
ax5.plot(years, hist_rule.sum(axis=1), "-.", color="#ff7f0e", linewidth=2,
         label="Rule of thumb")
ax5.plot(years, hist_opt.sum(axis=1), "-", color="#9467bd", linewidth=3, label="Optimal")
ax5.set_xlabel("Year")
ax5.set_ylabel("Total population (k)")
ax5.set_title("Total Population by Strategy")
ax5.grid(alpha=0.3)
ax5.legend()

# (6) Housing pressure (population / housing capacity)
ax6 = fig.add_subplot(2, 3, 6)
H_t = np.array([H0 + dH * t for t in range(YEARS + 1)])
for i, name in enumerate(CITIES):
    ax6.plot(years, hist_base[:, i] / H_t[:, i], "--", color=colors[i])
    ax6.plot(years, hist_opt[:, i] / H_t[:, i], "-", color=colors[i], linewidth=2.5,
             label=name)
ax6.axhline(1.0, color="black", linewidth=1)
ax6.set_xlabel("Year")
ax6.set_ylabel("Population / housing capacity")
ax6.set_title("Housing Pressure (dashed: no move, solid: optimal)")
ax6.grid(alpha=0.3)
ax6.legend()

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

Code Walkthrough

1. Problem settings

All city parameters are stored as NumPy arrays of length three, so every formula in the model is applied to all cities at once. P0 is the initial population, H0/dH describe housing and its yearly growth, J0/dJ do the same for jobs, and F0 is food production. MOVE_LOSS is the settling loss $\lambda$.

2. The simulation function

simulate implements the update rule $P_{i,t+1}=P_{i,t}(1+bq_{i,t}-d)$. Three details matter:

Vectorization. pop_start may be a single vector of shape (3,) or a matrix of shape (N, 3) holding $N$ scenarios. The food factor uses sum(axis=-1, keepdims=True), which works for both shapes. This is the key to the speed-up: instead of looping over thousands of scenarios in Python, NumPy evaluates them all in one pass.
Min functions. np.minimum(1.0, H / P) caps each satisfaction factor at one. Extra housing above the population does not make life better, but a shortage reduces growth.
Time-varying capacity. Housing and jobs grow linearly, and food grows geometrically, so a good migration plan must anticipate the future, not just today’s situation.

3. The migration model

The six decision variables are encoded as [out_I, out_A, out_C, split_I, split_A, split_C]. flow_matrix converts them into a $3\times3$ matrix $M$, where $M_{ij}$ is the number of people moving from city $i$ to $j$. Because a city can never send away more than 60% of its residents, every candidate is automatically feasible, and no constraint handling is needed. after_migration subtracts the leavers, adds the arrivals after the settling loss, and returns the new starting population.

4. Optimization

differential_evolution is a population-based global optimizer. It is well suited here because the objective has kinks from the $\min$ operations and may have several local optima. A fixed seed makes the result reproducible, and polish=True refines the final answer with a local search.

5. Comparison strategies

Two baselines put the optimum in context: no migration, and a rule of thumb that simply moves 20% of the agricultural population to the industrial city, following the intuition that jobs are where the industry is.

6. The 3D surface

To visualize the landscape, we vary two flows, Agricultural → Industrial and Agricultural → Commercial, over a $61\times61$ grid while keeping the other flows at their optimal values. All 3,721 scenarios are stacked into one (3721, 3) array and simulated in a single call. Combinations that would move more people than the agricultural city has are masked with NaN.

7. One combined figure

Six panels are drawn on a single canvas, so one plt.show() produces the entire visual summary.

Console Output

================================================================
 Population Migration and Urban Growth (unit: thousand people)
================================================================
City             Initial     No move        Rule     Optimal
Industrial         500.0       672.0       745.8       684.5
Agricultural       300.0       381.1       312.1       241.3
Commercial         200.0       268.8       268.8       411.4
----------------------------------------------------------------
Total             1000.0      1321.8      1326.7      1337.3

Optimal migration flows (rows: from, columns: to, thousand people)
                  Industrial  Agricultural    Commercial
Industrial               0.0           0.0           0.0
Agricultural             9.7           0.0         110.6
Commercial               0.0           0.0           0.0

Gain over no migration : 15.4 thousand (1.17 %)
Gain over rule of thumb: 10.6 thousand (0.80 %)

Grid best (A->I, A->C): (20.0, 100.0) -> total 1337.3 thousand

Graph Output

Interpreting the Results

Overall comparison

Doing nothing yields about 1,322 thousand people after 10 years. The simple rule of thumb of sending 20% of the farmers to the industrial city reaches about 1,327 thousand, while the optimized plan reaches about 1,337 thousand. The optimum beats the rule of thumb by roughly 0.8%. In a model where the baseline already grows by about 32%, that is a meaningful extra gain obtained purely by relocating people at the start.

The optimal plan

The optimizer moves almost nobody into the industrial city. Instead, it sends roughly 111 thousand people from the agricultural city to the commercial city and only about 10 thousand to the industrial city. This is the opposite of the naive intuition. The industrial city has many jobs, but its housing is already tight, so extra residents there are penalized by the housing factor $h_{i,t}$. The commercial city strikes the best balance between jobs and housing, so it absorbs most of the newcomers.

The 3D surface and contour map

The surface has a clear peak, and the red marker sits on top of it. Moving too few people leaves the agricultural city’s spare housing unused, while moving too many pushes the destination into overcrowding, where growth collapses. The ridge runs roughly along the line where the total number of migrants stays constant. This shows that the total number of migrants matters more than how they are split, with the split fine-tuning the result. The contour map shows the same shape from above and makes the optimal region easy to read.

Flow heatmap

The heatmap confirms that the optimal plan is essentially a single-corridor flow from the agricultural city to the commercial city, with all other flows near zero.

City populations

Compared with the dashed no-migration curves, the commercial city grows far faster under the optimal plan, while the agricultural city shrinks. The industrial city ends up close to its no-migration level, which is exactly the point: it was already near its comfortable capacity.

Total population by strategy

The optimal curve starts slightly lower because of the settling loss, then overtakes the other strategies as the better housing match compounds year after year. This is a classic trade-off: a small up-front cost for a larger long-term gain.

Housing pressure

The ratio of population to housing capacity shows why the plan works. Under no migration, the industrial city is already above 1, meaning a housing shortage, and keeps growing worse. Under the optimal plan, the pressure is spread more evenly across cities, which keeps the quality factor $q_{i,t}$ high and growth strong everywhere.

Conclusion

Maximizing long-term population is not about sending people to where the jobs are, nor to where the food is. It is about matching people to places that can both employ and house them, while accounting for how capacity changes over time. By combining a simple growth model, a vectorized simulation, and differential evolution, we found a migration plan that beats both doing nothing and a plausible rule of thumb. The same framework extends easily to more cities, yearly migration decisions, or additional constraints such as transport capacity and budgets.

Optimal Water Allocation Across Six Villages with Linear Programming and Python

October 30, 2026

Upstream vs. Downstream

Introduction

Picture a river that flows through six villages. The three villages at the top (upstream) get first access to the water. The three at the bottom (downstream) get whatever is left, plus a tributary that joins the river midway. Every liter pulled out upstream is a liter that downstream cannot use, yet farming and drinking water are both essential to every village.

This article formulates the problem as an optimization model, solves it exactly with linear programming, and compares the result with an uncoordinated “first come, first served” policy. We then use 3D and 2D visualizations to show how the optimal allocation moves as the upstream intake limit and the river inflow change.

Problem Setting

Villages are indexed $i = 1,\dots,6$ from upstream to downstream. Villages $1$–$3$ are upstream, and villages $4$–$6$ are downstream.
The river inflow at the top is $Q_0$. A tributary with flow $Q_T$ joins just after village 3.
Village $i$ withdraws $a_i$ for agriculture and $w_i$ for domestic use, in thousand m³/day.
A fraction $c$ of every withdrawal is consumed. The remaining $1-c$ returns to the river.
At the river mouth, an environmental minimum flow $Q_{\mathrm{env}}$ must remain.
The goal is to maximize the total population supported across all villages.

Mathematical Formulation

Let $P_i$ be the population (in thousands) that village $i$ can sustain. Food production saturates because arable land is limited:

$$
P_i^{\mathrm{food}}(a_i) = P_i^{\max}\left(1 - e^{-a_i/s_i}\right)
$$

Here $P_i^{\max}$ is the maximum population the land of village $i$ can feed, and $s_i$ is the water scale of its irrigation. Domestic water supports

$$
P_i^{\mathrm{dom}}(w_i) = \frac{w_i}{d}
$$

where $d$ is the per-capita domestic demand. The population is limited by the scarcer of the two:

The flow arriving at village $i$ is

$$
Q_i = Q_0 + \mathbb{1}[i \ge 4],Q_T - c\sum_{j<i}\left(a_j + w_j\right)
$$

The complete optimization problem is

$$
\begin{aligned}
\max_{a,,w,,P}\quad & \sum_{i=1}^{6} P_i \
\text{s.t.}\quad & P_i \le P_i^{\max}\left(1 - e^{-a_i/s_i}\right), \
& d,P_i \le w_i, \
& a_i + w_i \le Q_i, \
& Q_0 + Q_T - c\sum_{j=1}^{6}\left(a_j + w_j\right) \ge Q_{\mathrm{env}}, \
& a_i,\ w_i,\ P_i \ge 0
\end{aligned}
$$

The function $P_i^{\mathrm{food}}$ is concave, so the first constraint defines a convex set. Every other constraint is linear. We can therefore replace the concave curve by its tangent lines at $K$ points $\alpha_k$:

$$
P_i \le P_i^{\mathrm{food}}(\alpha_k) + \frac{dP_i^{\mathrm{food}}}{da}\bigg|_{\alpha_k}\left(a_i - \alpha_k\right), \qquad k = 1,\dots,K
$$

with

$$
\frac{dP_i^{\mathrm{food}}}{da} = \frac{P_i^{\max}}{s_i},e^{-a/s_i}
$$

The minimum of these tangents approximates the concave curve from above, and the error shrinks as $K$ grows. The problem becomes a pure linear program that HiGHS solves in milliseconds, so the 3D parameter sweep later in the article stays fast.

Example Parameters

Village	Zone	$P_i^{\max}$ (thousand people)	$s_i$ (thousand m³/day)
U1	Upstream	6	25
U2	Upstream	8	30
U3	Upstream	5	20
D1	Downstream	9	28
D2	Downstream	7	35
D3	Downstream	10	40

The remaining parameters are:

Domestic demand: $d = 0.15$ thousand m³/day per thousand people (150 L per person per day)
Consumptive fraction: $c = 0.6$
River inflow: $Q_0 = 100$ thousand m³/day
Tributary flow: $Q_T = 30$ thousand m³/day
Environmental minimum flow: $Q_{\mathrm{env}} = 40$ thousand m³/day

As a baseline we use an uncoordinated policy. Each village, starting from the top, takes enough water to support 92% of its maximum population (domestic water first). It is limited only by the water physically arriving and by the total consumable budget $Q_0 + Q_T - Q_{\mathrm{env}}$.

Source Code

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import gridspec
from scipy.optimize import linprog

# ------------------------------------------------------------
# Parameters
# ------------------------------------------------------------
NAMES = ["U1", "U2", "U3", "D1", "D2", "D3"]
N = len(NAMES)
N_UP = 3
P_MAX = np.array([6.0, 8.0, 5.0, 9.0, 7.0, 10.0])
S = np.array([25.0, 30.0, 20.0, 28.0, 35.0, 40.0])
D = 0.15
C = 0.6
Q0 = 100.0
Q_TRIB = 30.0
Q_ENV = 40.0
K = 80
BIG = 1.0e6


def food_pop(a):
    return P_MAX * (1.0 - np.exp(-a / S))


def population(a, w):
    return np.minimum(food_pop(a), w / D)


# ------------------------------------------------------------
# Linear program (variables: a[0:N], w[N:2N], P[2N:3N])
# ------------------------------------------------------------
rows_tan = []
rhs_tan = []
for i in range(N):
    alpha = np.linspace(0.0, 6.0 * S[i], K)
    f = P_MAX[i] * (1.0 - np.exp(-alpha / S[i]))
    m = P_MAX[i] / S[i] * np.exp(-alpha / S[i])
    b = f - m * alpha
    for mk, bk in zip(m, b):
        row = np.zeros(3 * N)
        row[2 * N + i] = 1.0
        row[i] = -mk
        rows_tan.append(row)
        rhs_tan.append(bk)
A_tan = np.array(rows_tan)
b_tan = np.array(rhs_tan)

A_dom = np.zeros((N, 3 * N))
for i in range(N):
    A_dom[i, 2 * N + i] = D
    A_dom[i, N + i] = -1.0

A_flow = np.zeros((N, 3 * N))
for i in range(N):
    A_flow[i, i] += 1.0
    A_flow[i, N + i] += 1.0
    for j in range(i):
        A_flow[i, j] += C
        A_flow[i, N + j] += C

A_env = np.zeros((1, 3 * N))
for j in range(N):
    A_env[0, j] = C
    A_env[0, N + j] = C

A_cap = np.zeros((1, 3 * N))
for j in range(N_UP):
    A_cap[0, j] = 1.0
    A_cap[0, N + j] = 1.0

A_UB = np.vstack([A_tan, A_dom, A_flow, A_env, A_cap])
COST = np.concatenate([np.zeros(2 * N), -np.ones(N)])
BOUNDS = [(0.0, None)] * (3 * N)


def solve_lp(q0, cap=BIG):
    q_in = np.where(np.arange(N) < N_UP, q0, q0 + Q_TRIB)
    b_ub = np.concatenate([b_tan, np.zeros(N), q_in,
                           [q0 + Q_TRIB - Q_ENV], [cap]])
    res = linprog(COST, A_ub=A_UB, b_ub=b_ub, bounds=BOUNDS, method="highs")
    if not res.success:
        return None
    x = np.clip(res.x, 0.0, None)
    return x[:N], x[N:2 * N], x[2 * N:]


# ------------------------------------------------------------
# Uncoordinated (first come, first served) baseline
# ------------------------------------------------------------
def greedy_allocation(q0, target=0.92):
    budget = q0 + Q_TRIB - Q_ENV
    cum = 0.0
    a = np.zeros(N)
    w = np.zeros(N)
    for i in range(N):
        q_in = q0 if i < N_UP else q0 + Q_TRIB
        arrival = q_in - cum
        w_want = D * target * P_MAX[i]
        a_want = -S[i] * np.log(1.0 - target)
        limit = max(min(arrival, (budget - cum) / C), 0.0)
        take = min(w_want + a_want, limit)
        w[i] = min(w_want, take)
        a[i] = take - w[i]
        cum += C * take
    return a, w


def river_flow(a, w, q0):
    cum = np.cumsum(a + w) * C
    trib = np.where(np.arange(N) >= N_UP, Q_TRIB, 0.0)
    return q0 + trib - cum


# ------------------------------------------------------------
# Base case
# ------------------------------------------------------------
sol = solve_lp(Q0)
if sol is None:
    raise RuntimeError("The linear program could not be solved.")
a_opt, w_opt, _ = sol
pop_opt = population(a_opt, w_opt)

a_gr, w_gr = greedy_allocation(Q0)
pop_gr = population(a_gr, w_gr)

flow_opt = river_flow(a_opt, w_opt, Q0)
flow_gr = river_flow(a_gr, w_gr, Q0)

df = pd.DataFrame({
    "Village": NAMES,
    "Zone": ["Upstream"] * N_UP + ["Downstream"] * (N - N_UP),
    "Agri_opt": a_opt,
    "Dom_opt": w_opt,
    "Pop_opt": pop_opt,
    "Agri_greedy": a_gr,
    "Dom_greedy": w_gr,
    "Pop_greedy": pop_gr,
}).round(2)

print("=== Allocation (thousand m3/day) and population (thousand people) ===")
print(df.to_string(index=False))
print()
print(f"Total population (optimal)       : {pop_opt.sum():.2f}")
print(f"Total population (uncoordinated) : {pop_gr.sum():.2f}")
print(f"Improvement                      : {100.0 * (pop_opt.sum() / pop_gr.sum() - 1.0):.1f} %")
print(f"Upstream population (optimal)    : {pop_opt[:N_UP].sum():.2f}")
print(f"Downstream population (optimal)  : {pop_opt[N_UP:].sum():.2f}")
print(f"Upstream population (uncoord.)   : {pop_gr[:N_UP].sum():.2f}")
print(f"Downstream population (uncoord.) : {pop_gr[N_UP:].sum():.2f}")
print(f"River mouth flow (optimal)       : {flow_opt[-1]:.2f}")
print(f"River mouth flow (uncoordinated) : {flow_gr[-1]:.2f}")

# ------------------------------------------------------------
# Sweep: upstream intake cap x river inflow
# ------------------------------------------------------------
caps = np.linspace(0.0, 160.0, 33)
q0_grid = np.linspace(60.0, 160.0, 21)
TOT = np.full((len(q0_grid), len(caps)), np.nan)
UPP = np.full((len(q0_grid), len(caps)), np.nan)
DNP = np.full((len(q0_grid), len(caps)), np.nan)

for qi, q in enumerate(q0_grid):
    for ci, cap in enumerate(caps):
        res_sweep = solve_lp(q, cap)
        if res_sweep is None:
            continue
        p = population(res_sweep[0], res_sweep[1])
        TOT[qi, ci] = p.sum()
        UPP[qi, ci] = p[:N_UP].sum()
        DNP[qi, ci] = p[N_UP:].sum()

X, Y = np.meshgrid(caps, q0_grid)
qi_base = int(np.argmin(np.abs(q0_grid - Q0)))
up_use_opt = float((a_opt + w_opt)[:N_UP].sum())

# ------------------------------------------------------------
# Visualization (single figure)
# ------------------------------------------------------------
plt.style.use("dark_background")
fig = plt.figure(figsize=(21, 12))
gs = gridspec.GridSpec(2, 3, figure=fig, hspace=0.32, wspace=0.28,
                       top=0.92, bottom=0.06, left=0.04, right=0.97)

# (1) 3D surface
ax1 = fig.add_subplot(gs[0, 0], projection="3d")
surf = ax1.plot_surface(X, Y, TOT, cmap="viridis", edgecolor="none", alpha=0.95)
ax1.set_xlabel("Upstream intake cap")
ax1.set_ylabel("River inflow Q0")
ax1.set_zlabel("Total population")
ax1.set_title("Optimal total population surface")
ax1.view_init(elev=28, azim=-125)
fig.colorbar(surf, ax=ax1, shrink=0.55, pad=0.08)

# (2) Contour map
ax2 = fig.add_subplot(gs[0, 1])
cs = ax2.contourf(X, Y, TOT, levels=20, cmap="viridis")
fig.colorbar(cs, ax=ax2)
ax2.scatter([up_use_opt], [Q0], marker="*", s=300, color="#ff4d6d",
            edgecolor="white", zorder=5, label="Base-case optimum")
ax2.set_xlabel("Upstream intake cap")
ax2.set_ylabel("River inflow Q0")
ax2.set_title("Top view: total population")
ax2.legend(loc="lower right")

# (3) Upstream vs downstream trade-off
ax3 = fig.add_subplot(gs[0, 2])
ax3.stackplot(caps, UPP[qi_base, :], DNP[qi_base, :],
              labels=["Upstream population", "Downstream population"],
              colors=["#4cc9f0", "#f72585"], alpha=0.85)
ax3.axvline(up_use_opt, color="white", linestyle="--", linewidth=1.2,
            label="Unconstrained optimum")
ax3.set_xlabel("Upstream intake cap")
ax3.set_ylabel("Population (thousand)")
ax3.set_title("Upstream vs downstream trade-off (Q0 = 100)")
ax3.legend(loc="lower right")

# (4) Population by village
ax4 = fig.add_subplot(gs[1, 0])
xs = np.arange(N)
ax4.bar(xs - 0.2, pop_opt, width=0.4, color="#4cc9f0", label="Optimal")
ax4.bar(xs + 0.2, pop_gr, width=0.4, color="#f72585", label="Uncoordinated")
ax4.set_xticks(xs)
ax4.set_xticklabels(NAMES)
ax4.set_ylabel("Population (thousand)")
ax4.set_title("Population by village")
ax4.legend()

# (5) River flow profile
ax5 = fig.add_subplot(gs[1, 1])
ax5.plot(xs, flow_opt, "o-", color="#4cc9f0", linewidth=2, label="Optimal")
ax5.plot(xs, flow_gr, "s--", color="#f72585", linewidth=2, label="Uncoordinated")
ax5.axhline(Q_ENV, color="#ffd60a", linestyle=":", linewidth=1.5,
            label="Environmental minimum")
ax5.axvline(N_UP - 0.5, color="gray", linestyle="--", linewidth=1)
ax5.text(N_UP - 0.45, ax5.get_ylim()[1] * 0.97, "Tributary joins",
         color="lightgray", va="top")
ax5.set_xticks(xs)
ax5.set_xticklabels(NAMES)
ax5.set_ylim(bottom=0.0)
ax5.set_ylabel("River flow after village")
ax5.set_title("River flow profile")
ax5.legend(loc="lower left")

# (6) 3D bars of agricultural withdrawal
ax6 = fig.add_subplot(gs[1, 2], projection="3d")
ax6.bar3d(np.arange(N), np.zeros(N), np.zeros(N), 0.6, 0.5,
          np.asarray(a_opt, dtype=float), color="#4cc9f0", alpha=0.9)
ax6.bar3d(np.arange(N), np.ones(N), np.zeros(N), 0.6, 0.5,
          np.asarray(a_gr, dtype=float), color="#f72585", alpha=0.9)
ax6.set_xticks(np.arange(N) + 0.3)
ax6.set_xticklabels(NAMES)
ax6.set_yticks([0.25, 1.25])
ax6.set_yticklabels(["Optimal", "Uncoord."])
ax6.set_zlabel("Agricultural withdrawal")
ax6.set_title("Agricultural water by village")
ax6.view_init(elev=26, azim=-60)

fig.suptitle("Upstream and Downstream Water Allocation: Population Maximization",
             fontsize=18)
plt.show()

Code Walkthrough

Parameters and helper functions

The first block defines the village names, the number of upstream villages (N_UP = 3), and the physical parameters. food_pop evaluates the saturating food curve $P_i^{\max}(1 - e^{-a_i/s_i})$ for all villages at once with NumPy broadcasting. population takes the elementwise minimum of the food-limited and water-limited populations, which is exactly $P_i = \min{P_i^{\mathrm{food}}, w_i/d}$.

Building the linear program

The decision vector is laid out as $[a_1..a_6,\ w_1..w_6,\ P_1..P_6]$, which gives 18 variables.

Tangent constraints. For each village, K = 80 tangent points are placed from $0$ to $6 s_i$. Each tangent becomes one row $P_i - m_k a_i \le b_k$, where $m_k$ is the slope and $b_k = f(\alpha_k) - m_k \alpha_k$ is the intercept. Together these 480 rows approximate the concave food curve from above.
Domestic constraints. The rows $d P_i - w_i \le 0$ encode $P_i \le w_i/d$.
Flow constraints. Row $i$ encodes $a_i + w_i + c\sum_{j<i}(a_j + w_j) \le Q_i^{\mathrm{in}}$. This is the arrival-flow condition: a village cannot withdraw more than what reaches it after all upstream consumption.
Environmental constraint. A single row enforces $c\sum_j (a_j + w_j) \le Q_0 + Q_T - Q_{\mathrm{env}}$.
Upstream intake cap. A final row limits total upstream withdrawal to cap. It is set to a huge value (BIG) in the base case and varied during the sweep.

All of these matrices are built once outside the solver function. Only the right-hand side changes between solves, so the 693 linear programs in the parameter sweep run quickly. linprog(..., method="highs") returns the global optimum, and the minus signs in COST turn maximization into minimization. The function returns None if the solver fails, and the main code raises a clear error in that case.

Uncoordinated baseline

greedy_allocation simulates the “first come, first served” policy. Walking downstream, each village wants enough water for 92% of its maximum population (domestic water first, then irrigation). It receives the smaller of its want, the water that physically arrives, and the remaining consumable budget. Because upstream villages are served first, they can exhaust the whole budget before the water reaches the lower valley.

River profile

river_flow computes the flow after each village as $Q_0$, plus the tributary for downstream villages, minus the cumulative consumption $c\sum(a_j + w_j)$. This is the quantity compared against the environmental minimum.

Parameter sweep

The double loop solves the LP for 21 river inflows (60 to 160) and 33 upstream intake caps (0 to 160). The results are stored in three matrices holding the total, upstream, and downstream populations. These matrices feed the 3D surface, the contour map, and the stacked trade-off plot.

Visualization

All six panels are drawn in one figure with a dark theme.

A 3D surface of the optimal total population versus intake cap and river inflow.
A top-down contour view of the same surface, with the base-case optimum marked by a star.
A stacked area chart of upstream and downstream population versus the intake cap at $Q_0 = 100$.
Grouped bars of the population of each village under both policies.
The river flow profile along the valley, with the environmental minimum and the tributary junction marked.
A 3D bar chart comparing agricultural withdrawals in each village under both policies.

Execution Results

Result Image

Console Output

=== Allocation (thousand m3/day) and population (thousand people) ===
Village       Zone  Agri_opt  Dom_opt  Pop_opt  Agri_greedy  Dom_greedy  Pop_greedy
     U1   Upstream     19.92     0.49     3.30        63.14        0.83        5.52
     U2   Upstream     26.19     0.70     4.66        60.51        1.10        6.94
     U3   Upstream     15.94     0.41     2.75        23.72        0.69        3.47
     D1 Downstream     30.82     0.90     6.01         0.00        0.00        0.00
     D2 Downstream     21.43     0.48     3.21         0.00        0.00        0.00
     D3 Downstream     31.88     0.82     5.49         0.00        0.00        0.00

Total population (optimal)       : 25.41
Total population (uncoordinated) : 15.93
Improvement                      : 59.5 %
Upstream population (optimal)    : 10.70
Downstream population (optimal)  : 14.71
Upstream population (uncoord.)   : 15.93
Downstream population (uncoord.) : 0.00
River mouth flow (optimal)       : 40.00
River mouth flow (uncoordinated) : 40.00

Interpreting the Results

The 3D surface and contour map

The surface shows the maximum achievable total population for every combination of upstream intake cap and river inflow. Along the river-inflow axis, the surface rises as more water enters the system. Along the intake-cap axis, it rises and then flattens into a plateau. Beyond the plateau, relaxing the cap further adds nothing, because the optimizer would not use that extra upstream water anyway. In the contour view, the plateau appears as vertical contour lines. The star marks the base-case optimum at the point where the plateau begins.

The upstream versus downstream trade-off

The stacked area chart makes the conflict visible. Raising the upstream intake cap increases upstream population, but the downstream layer shrinks because less water arrives. The total is largest when the two effects are balanced. Since only 60% of each withdrawal is consumed, the return flow softens the conflict, and the optimizer takes advantage of it by directing water to the villages with the best population gain per unit of water.

Village populations and agricultural water

Under the uncoordinated policy, the upstream villages fill their demand first. The consumable budget is exhausted by the time the water reaches the lower valley, so the downstream villages receive close to nothing. The optimal plan does the opposite of “fill the top first”. Because the food curve saturates, the last units of water are worth little to a village that is already well supplied. The optimizer therefore holds back upstream intake once the marginal gain falls and sends that water downstream, where it supports far more people. The 3D bar chart shows this clearly: the pattern of agricultural withdrawals under the optimal plan is spread across the whole valley, while the uncoordinated plan is heavily concentrated at the top.

The river flow profile

Both policies respect the environmental minimum at the river mouth, but they use the water very differently. The uncoordinated policy drains the river to the environmental limit within the upstream section. The optimal policy keeps more water flowing into the lower valley, so the flow declines gradually along the whole river before reaching the same floor at the mouth.

Conclusion

We modeled a river shared by upstream and downstream villages and recast the nonlinear food-production curves as a linear program using tangent-line approximation. The result is an exact and extremely fast solver that finds the globally optimal split between agricultural and domestic water. Compared with a first-come-first-served policy, the optimal allocation supports substantially more people, because it recognizes that water is more valuable where it is scarce.

The same framework extends naturally to more villages, seasonal flow variation, canal capacity limits, or economic terms such as crop revenue. Any of these can be added as extra linear constraints or objective terms, and the solver remains fast.

Maximizing Total Population with Linear Programming

October 29, 2026

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

import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
from scipy.optimize import linprog
from scipy.sparse import coo_matrix

# ---------------------------------------------------------------
# 1. Problem data
# ---------------------------------------------------------------
T = 20
G_RATE = 0.35
FOOD_PER_CAP = 1.0
STOCK0 = 800.0
REFILL = 80.0
BUDGET = 6000.0
EPS = 1e-4
FLOOR = 0.6

names = ["Warehouse", "Aria", "Borea", "Cinder", "Dune", "Ember", "Fjord", "Glen"]
N = len(names)
NC = N - 1
pop0 = np.array([50, 80, 30, 60, 40, 20, 70], dtype=float)
prod = np.array([30, 40, 10, 20, 50, 5, 25], dtype=float)
K = np.array([120, 150, 100, 130, 90, 80, 140], dtype=float)
pos = {0: (0.0, 0.0), 1: (-2.0, 1.5), 2: (0.0, 2.0), 3: (2.0, 1.2),
       4: (-1.5, 3.5), 5: (1.5, 3.2), 6: (0.0, 4.6), 7: (2.8, 4.2)}

roads = [(0, 1, 90, 2), (0, 2, 120, 3), (0, 3, 80, 2),
         (1, 2, 30, 1), (1, 4, 40, 2), (2, 4, 50, 1), (2, 5, 45, 2),
         (3, 5, 40, 2), (3, 7, 25, 3), (4, 6, 35, 3), (5, 6, 30, 1),
         (5, 7, 25, 3), (6, 7, 30, 2)]

arcs = []
for u, v, cap, cost in roads:
    arcs.append((u, v, cap, cost))
    if u != 0:
        arcs.append((v, u, cap, cost))
nA = len(arcs)
arc_from = np.array([a[0] for a in arcs])
arc_to = np.array([a[1] for a in arcs])
arc_cap = np.array([a[2] for a in arcs], dtype=float)
arc_cost = np.array([a[3] for a in arcs], dtype=float)

# ---------------------------------------------------------------
# 2. Build the LP once (sparse matrices)
# ---------------------------------------------------------------
oF = nA * T
oP = oF + NC * T
oS = oP + NC * (T + 1)
nv = oS + T + 1


def ix(e, t): return e * T + t
def iF(j, t): return oF + j * T + t
def iP(j, t): return oP + j * (T + 1) + t
def iS(t): return oS + t


def build_model():
    er, ec, ev, eb = [], [], [], []
    r = 0
    for j in range(NC):
        for t in range(T):
            er += [r, r, r]
            ec += [iP(j, t + 1), iP(j, t), iF(j, t)]
            ev += [1.0, -(1.0 - G_RATE), -G_RATE / FOOD_PER_CAP]
            eb.append(0.0)
            r += 1
    for t in range(T):
        er += [r, r]
        ec += [iS(t + 1), iS(t)]
        ev += [1.0, -1.0]
        for e in range(nA):
            if arc_from[e] == 0:
                er.append(r)
                ec.append(ix(e, t))
                ev.append(1.0)
        eb.append(REFILL)
        r += 1
    A_eq = coo_matrix((ev, (er, ec)), shape=(r, nv)).tocsr()
    b_eq = np.array(eb)

    ur, uc, uv, ub = [], [], [], []
    r = 0
    for j in range(NC):
        node = j + 1
        for t in range(T):
            ur.append(r); uc.append(iF(j, t)); uv.append(1.0)
            for e in range(nA):
                if arc_to[e] == node:
                    ur.append(r); uc.append(ix(e, t)); uv.append(-1.0)
                if arc_from[e] == node:
                    ur.append(r); uc.append(ix(e, t)); uv.append(1.0)
            ub.append(prod[j])
            r += 1
    for e in range(nA):
        for t in range(T):
            ur.append(r); uc.append(ix(e, t)); uv.append(arc_cost[e])
    ub.append(BUDGET)
    r += 1
    A_ub = coo_matrix((uv, (ur, uc)), shape=(r, nv)).tocsr()
    b_ub = np.array(ub)

    c = np.zeros(nv)
    for j in range(NC):
        c[iP(j, T)] = -1.0
    for e in range(nA):
        for t in range(T):
            c[ix(e, t)] = EPS * arc_cost[e]
    return A_eq, b_eq, A_ub, b_ub, c


A_EQ, B_EQ, A_UB, B_UB, C_VEC = build_model()


def solve(budget=BUDGET, stock0=STOCK0, cap_scale=1.0):
    lb = np.zeros(nv)
    ubd = np.full(nv, np.inf)
    for e in range(nA):
        for t in range(T):
            ubd[ix(e, t)] = arc_cap[e] * cap_scale
    for j in range(NC):
        for t in range(T + 1):
            ubd[iP(j, t)] = K[j]
            lb[iP(j, t)] = FLOOR * pop0[j]
        lb[iP(j, 0)] = pop0[j]
        ubd[iP(j, 0)] = pop0[j]
    lb[iS(0)] = stock0
    ubd[iS(0)] = stock0
    b_ub = B_UB.copy()
    b_ub[-1] = budget
    res = linprog(C_VEC, A_ub=A_UB, b_ub=b_ub, A_eq=A_EQ, b_eq=B_EQ,
                  bounds=np.column_stack([lb, ubd]), method="highs")
    if res.status != 0:
        raise RuntimeError("LP failed: " + res.message)
    x = res.x
    flow = x[:oF].reshape(nA, T)
    F = x[oF:oP].reshape(NC, T)
    P = x[oP:oS].reshape(NC, T + 1)
    S = x[oS:]
    return dict(flow=flow, F=F, P=P, S=S, total=P[:, T].sum(),
                cost=float((arc_cost[:, None] * flow).sum()))


# ---------------------------------------------------------------
# 3. Solve and report
# ---------------------------------------------------------------
t0 = time.time()
opt = solve()
t_single = time.time() - t0

print("=" * 60)
print("Optimal plan (relay roads allowed)")
print("=" * 60)
print(f"Solve time                : {t_single:.3f} s")
print(f"Initial total population  : {pop0.sum():.1f} thousand")
print(f"Final total population    : {opt['total']:.1f} thousand")
print(f"Transport cost used       : {opt['cost']:.1f} / {BUDGET:.0f}")
print(f"Warehouse stock at turn 20: {opt['S'][-1]:.1f} t")
print("-" * 60)
print(f"{'City':<8}{'Initial':>10}{'Final':>10}{'Capacity':>10}{'Growth %':>10}")
for j in range(NC):
    p_end = opt["P"][j, T]
    print(f"{names[j+1]:<8}{pop0[j]:>10.1f}{p_end:>10.1f}{K[j]:>10.1f}"
          f"{(p_end / pop0[j] - 1) * 100:>10.1f}")

# Sensitivity sweep (3D graph)
budgets = np.linspace(3500, 9000, 8)
cap_scales = np.linspace(0.8, 2.0, 8)
Z = np.zeros((len(cap_scales), len(budgets)))
t0 = time.time()
for a, cs in enumerate(cap_scales):
    for b, bd in enumerate(budgets):
        Z[a, b] = solve(budget=bd, cap_scale=cs)["total"]
print("-" * 60)
print(f"Sweep of {Z.size} LPs finished in {time.time() - t0:.2f} s")

# ---------------------------------------------------------------
# 4. Visualization (one figure)
# ---------------------------------------------------------------
fig = plt.figure(figsize=(22, 13))
gs = fig.add_gridspec(2, 3, hspace=0.30, wspace=0.38)
colors = plt.cm.tab10(np.arange(NC))

# (1) network
ax1 = fig.add_subplot(gs[0, 0])
tot_flow = opt["flow"].sum(axis=1)
util = tot_flow / (arc_cap * T)
max_f = max(tot_flow.max(), 1.0)
for e in range(nA):
    if tot_flow[e] < 1e-6:
        continue
    x0, y0 = pos[arc_from[e]]
    x1, y1 = pos[arc_to[e]]
    ax1.annotate("", xy=(x1, y1), xytext=(x0, y0),
                 arrowprops=dict(arrowstyle="-|>", lw=1 + 7 * tot_flow[e] / max_f,
                                 color=cm.viridis(util[e]), alpha=0.85,
                                 shrinkA=18, shrinkB=18,
                                 connectionstyle="arc3,rad=0.08"))
    ax1.text((x0 + x1) / 2, (y0 + y1) / 2, f"{tot_flow[e]:.0f}", fontsize=8,
             ha="center", va="center",
             bbox=dict(boxstyle="round,pad=0.15", fc="white", ec="none", alpha=0.8))
for n in range(N):
    x, y = pos[n]
    if n == 0:
        ax1.scatter(x, y, s=1200, c="gold", edgecolors="k", zorder=3, marker="s")
    else:
        ax1.scatter(x, y, s=300 + 6 * opt["P"][n - 1, T], c=[colors[n - 1]],
                    edgecolors="k", zorder=3)
    ax1.text(x, y - 0.38, names[n], ha="center", fontsize=9, fontweight="bold")
sm = cm.ScalarMappable(norm=Normalize(0, 1), cmap="viridis")
sm.set_array([])
fig.colorbar(sm, ax=ax1, fraction=0.04, pad=0.02, label="Road utilization")
ax1.set_title("Total food shipped per road (20 turns)")
ax1.axis("off")

# (2) trajectories
ax2 = fig.add_subplot(gs[0, 1])
turns = np.arange(T + 1)
for j in range(NC):
    ax2.plot(turns, opt["P"][j], color=colors[j], lw=2, label=names[j + 1])
    ax2.axhline(K[j], color=colors[j], lw=0.6, ls=":")
ax2.set_xlabel("Turn")
ax2.set_ylabel("Population (thousand)")
ax2.set_title("Population trajectories (dotted = capacity)")
ax2.legend(ncol=2, fontsize=8)
ax2.grid(alpha=0.3)

# (3) 3D population surface
ax3 = fig.add_subplot(gs[0, 2], projection="3d")
Xg, Yg = np.meshgrid(turns, np.arange(NC))
ax3.plot_surface(Xg, Yg, opt["P"], cmap="viridis", edgecolor="k", lw=0.3, alpha=0.9)
ax3.set_xlabel("Turn")
ax3.set_yticks(np.arange(NC))
ax3.set_yticklabels([names[j + 1] for j in range(NC)], fontsize=7)
ax3.set_zlabel("Population")
ax3.set_title("3D: population by city and turn")
ax3.view_init(elev=28, azim=-60)

# (4) heatmap
ax4 = fig.add_subplot(gs[1, 0])
used = [e for e in range(nA) if opt["flow"][e].sum() > 1e-6]
heat = np.array([opt["flow"][e] / arc_cap[e] for e in used])
im = ax4.imshow(heat, aspect="auto", cmap="magma", vmin=0, vmax=1)
ax4.set_yticks(range(len(used)))
ax4.set_yticklabels([f"{names[arc_from[e]]}→{names[arc_to[e]]}" for e in used], fontsize=8)
ax4.set_xticks(range(0, T, 2))
ax4.set_xticklabels(range(1, T + 1, 2))
ax4.set_xlabel("Turn")
ax4.set_title("Road utilization heatmap (flow / capacity)")
fig.colorbar(im, ax=ax4, fraction=0.04, pad=0.02)

# (5) stock and cost
ax5 = fig.add_subplot(gs[1, 1])
ax5.plot(turns, opt["S"], color="tab:orange", lw=2.5, label="Warehouse stock (t)")
ax5.set_xlabel("Turn")
ax5.set_ylabel("Stock (t)", color="tab:orange")
ax5.grid(alpha=0.3)
ax5b = ax5.twinx()
cum_cost = np.concatenate([[0], np.cumsum((arc_cost[:, None] * opt["flow"]).sum(axis=0))])
ax5b.plot(turns, cum_cost, color="tab:red", lw=2.5, ls="--", label="Cumulative cost")
ax5b.axhline(BUDGET, color="gray", ls=":")
ax5b.set_ylabel("Cumulative cost", color="tab:red")
h1, l1 = ax5.get_legend_handles_labels()
h2, l2 = ax5b.get_legend_handles_labels()
ax5.legend(h1 + h2, l1 + l2, loc="center right")
ax5.set_title("Warehouse stock and transport cost")

# (6) 3D sensitivity
ax6 = fig.add_subplot(gs[1, 2], projection="3d")
Bg, Cg = np.meshgrid(budgets, cap_scales)
ax6.plot_surface(Bg, Cg, Z, cmap="plasma", edgecolor="k", lw=0.3, alpha=0.9)
ax6.set_xlabel("Budget")
ax6.set_ylabel("Road capacity multiplier")
ax6.set_zlabel("Final population")
ax6.set_title("3D: final population vs. budget and road capacity")
ax6.view_init(elev=26, azim=-125)

fig.suptitle("Resource allocation to multiple cities: 20-turn LP optimization",
             fontsize=17, fontweight="bold", y=0.96)
plt.show()

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 c contains $-1$ for each final-turn population, because linprog minimizes, 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.

Optimal Immigration and Birth Support Over 30 Turns

October 28, 2026

Steering Population Growth:

Introduction

Every growing society faces the same dilemma. Grow too fast, and the food supply collapses under the weight of new mouths. Grow too slowly, and the workforce is too small to produce enough food to sustain a prosperous economy.

In this article we turn that dilemma into a concrete optimization problem. At every turn, a planner chooses two quantities:

the number of immigrants to accept, and
the birth-support budget to spend.

The goal is to maximize the cumulative population over 30 turns. We build a small simulation model, solve it numerically in Python, and visualize the result with a 3D landscape and a 3D trajectory plot.

Problem Setting

Let $t = 0, 1, \dots, T-1$ with $T = 30$. The state consists of the population $P_t$ and the food stock $S_t$. The controls are normalized intensities $x_t, y_t \in [0, 1]$:

$$
I_t = I_{\max}, x_t, \qquad B_t = B_{\max}, y_t
$$

where $I_t$ is the requested number of immigrants and $B_t$ is the requested birth-support spending.

Production and food balance

Only a fraction $\lambda$ of the population works. Production has diminishing returns in labor:

$$
L_t = \lambda P_t, \qquad Y_t = A, L_t^{\beta}, \qquad 0 < \beta < 1
$$

Each person needs $c$ units of food per turn. The satisfaction ratio and the next food stock (with spoilage rate $\delta$ and storage capacity $S_{\max}$) are

$$
s_t = \min\left(1,\ \frac{S_t + Y_t}{c,P_t}\right), \qquad
S_{t+1} = \min\left(S_{\max},\ (1-\delta)\max\left(0,\ S_t + Y_t - c,P_t\right)\right)
$$

Budget constraint

The treasury collects a tax $\tau$ on production. Each immigrant costs $c_I$ (housing and integration). If the requested spending exceeds the budget, both requests are scaled down proportionally:

$$
\kappa_t = \min\left(1,\ \frac{\tau Y_t}{c_I I_t + B_t}\right), \qquad
\tilde I_t = \kappa_t I_t, \qquad \tilde B_t = \kappa_t B_t
$$

Demographics

Birth support has diminishing returns per capita, while food shortage raises mortality:

$$
b_t = b_0 + b_1\left(1 - \exp\left(-\frac{\tilde B_t}{\sigma P_t}\right)\right), \qquad
d_t = \min\left(1,\ d_0 + d_1 (1 - s_t)\right)
$$

$$
P_{t+1} = \max\left(1,\ P_t,(1 + b_t - d_t) + \tilde I_t\right)
$$

Objective

$$
\max_{x_t,, y_t \in [0,1]} \ J = \sum_{t=1}^{T} P_t
$$

Parameters

Symbol	Meaning	Value
$T$	Number of turns	30
$P_0,\ S_0$	Initial population, initial food stock	1000, 1500
$\lambda$	Working share of the population	0.6
$A,\ \beta$	Productivity, labor elasticity	18, 0.7
$c$	Food need per person	1.0
$\delta,\ S_{\max}$	Spoilage rate, storage capacity	0.05, 3000
$\tau$	Tax rate on production	0.10
$c_I,\ I_{\max},\ B_{\max}$	Cost per immigrant, maximum requests	2.0, 150, 100
$b_0,\ b_1,\ \sigma$	Base birth rate, birth-support gain, per-capita scale	0.015, 0.030, 0.05
$d_0,\ d_1$	Base death rate, shortage mortality coefficient	0.010, 0.60

With these values, the steady-state carrying capacity (where $Y = cP$) is roughly 4,600 people. The initial population is far below it, so there is room to grow, but a fast surge overshoots the capacity and triggers shortage deaths.

Solution Strategy

Landscape analysis. Evaluate every constant policy $(x, y)$ on a $41 \times 41$ grid to see the shape of the objective in 3D.
Full optimization. Optimize all 60 decision variables (30 turns × 2 controls) with multi-start L-BFGS-B, seeded with the best constant policy.
Local polish. Refine the best solution with a sparse random local search that only accepts improvements.
Comparison. Compare against four baseline policies.

Because the objective is evaluated thousands of times, the simulator uses only scalar Python arithmetic with the math module, which avoids NumPy’s per-call overhead on scalars and keeps one simulation in the tens-of-microseconds range.

Full Source Code

import math
import time

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# ----------------------------------------------------------------------
# Model parameters
# ----------------------------------------------------------------------
T = 30
P0 = 1000.0
S0 = 1500.0
LABOR_SHARE = 0.6
A_PROD = 18.0
BETA = 0.7
NEED_PER_CAPITA = 1.0
SPOIL = 0.05
S_MAX = 3000.0
TAX_RATE = 0.10
COST_IMM = 2.0
I_MAX = 150.0
B_MAX = 100.0
BASE_BIRTH = 0.015
BIRTH_GAIN = 0.030
BIRTH_SCALE = 0.05
BASE_DEATH = 0.010
SHORT_DEATH = 0.60


def clip01(v):
    if v < 0.0:
        return 0.0
    if v > 1.0:
        return 1.0
    return v


# ----------------------------------------------------------------------
# Fast scalar simulator
# ----------------------------------------------------------------------
def simulate(x, y, trace=False):
    P = P0
    S = S0
    total = 0.0
    rec = None
    if trace:
        rec = {"P": [P0], "S": [S0], "sat": [], "imm": [], "bsp": [],
               "budget": [], "births": [], "deaths": [], "bind": []}
    for t in range(T):
        xt = clip01(float(x[t]))
        yt = clip01(float(y[t]))

        labor = LABOR_SHARE * P
        prod = A_PROD * labor ** BETA
        need = NEED_PER_CAPITA * P
        avail = S + prod
        if avail >= need:
            left = avail - need
            sat = 1.0
        else:
            left = 0.0
            sat = avail / need

        budget = TAX_RATE * prod
        imm = I_MAX * xt
        bsp = B_MAX * yt
        cost = COST_IMM * imm + bsp
        binding = 0.0
        if cost > budget and cost > 0.0:
            kappa = budget / cost
            imm *= kappa
            bsp *= kappa
            binding = 1.0

        birth_rate = BASE_BIRTH + BIRTH_GAIN * (1.0 - math.exp(-bsp / (BIRTH_SCALE * P)))
        death_rate = min(1.0, BASE_DEATH + SHORT_DEATH * (1.0 - sat))
        births = P * birth_rate
        deaths = P * death_rate

        P = max(P + births - deaths + imm, 1.0)
        S = min(left * (1.0 - SPOIL), S_MAX)
        total += P

        if trace:
            rec["P"].append(P)
            rec["S"].append(S)
            rec["sat"].append(sat)
            rec["imm"].append(imm)
            rec["bsp"].append(bsp)
            rec["budget"].append(budget)
            rec["births"].append(births)
            rec["deaths"].append(deaths)
            rec["bind"].append(binding)

    if trace:
        return total, {key: np.array(val) for key, val in rec.items()}
    return total


def objective(z):
    return -simulate(z[:T], z[T:])


# ----------------------------------------------------------------------
# Step 1: landscape of constant policies
# ----------------------------------------------------------------------
t_start = time.time()
grid = np.linspace(0.0, 1.0, 41)
Z = np.zeros((grid.size, grid.size))
for i, yv in enumerate(grid):
    for j, xv in enumerate(grid):
        Z[i, j] = simulate(np.full(T, xv), np.full(T, yv))
iy, ix = np.unravel_index(int(np.argmax(Z)), Z.shape)
best_cx = float(grid[ix])
best_cy = float(grid[iy])
best_const_val = float(Z[iy, ix])
t_landscape = time.time() - t_start

# ----------------------------------------------------------------------
# Step 2: multi-start L-BFGS-B over all 60 variables
# ----------------------------------------------------------------------
t_start = time.time()
rng = np.random.default_rng(42)
n_var = 2 * T
bounds = [(0.0, 1.0)] * n_var

best_z = np.concatenate([np.full(T, best_cx), np.full(T, best_cy)])
best_val = best_const_val

starts = [best_z.copy(),
          np.full(n_var, 0.5),
          np.full(n_var, 0.2),
          np.full(n_var, 0.8)]
starts += [rng.uniform(0.0, 1.0, n_var) for _ in range(5)]

for z0 in starts:
    res = minimize(objective, z0, method="L-BFGS-B", bounds=bounds,
                   options={"maxiter": 150, "eps": 1e-4})
    val = -float(res.fun)
    if val > best_val:
        best_val = val
        best_z = np.clip(np.asarray(res.x, dtype=float), 0.0, 1.0)

# ----------------------------------------------------------------------
# Step 3: sparse random local search (accept improvements only)
# ----------------------------------------------------------------------
cur_z = best_z.copy()
cur_val = best_val
for it in range(6000):
    sigma = 0.05 if it < 3000 else 0.01
    cand = cur_z.copy()
    idx = rng.integers(0, n_var, size=5)
    cand[idx] += rng.normal(0.0, sigma, size=5)
    cand = np.clip(cand, 0.0, 1.0)
    v = simulate(cand[:T], cand[T:])
    if v > cur_val:
        cur_z = cand
        cur_val = v

opt_x = cur_z[:T]
opt_y = cur_z[T:]
t_opt = time.time() - t_start

# ----------------------------------------------------------------------
# Step 4: policy comparison
# ----------------------------------------------------------------------
zeros = np.zeros(T)
ones = np.ones(T)
policies = {
    "No control": (zeros, zeros),
    "Max immigration": (ones, zeros),
    "Max birth support": (zeros, ones),
    "Best constant": (np.full(T, best_cx), np.full(T, best_cy)),
    "Optimized": (opt_x, opt_y),
}
colors = {
    "No control": "#9aa0a6",
    "Max immigration": "#ff7043",
    "Max birth support": "#66bb6a",
    "Best constant": "#42a5f5",
    "Optimized": "#ffd54f",
}

results = {}
for name, (px, py) in policies.items():
    tot, rec = simulate(px, py, trace=True)
    results[name] = (tot, rec)

# ----------------------------------------------------------------------
# Console report
# ----------------------------------------------------------------------
print("=" * 72)
print("Population growth control: results")
print("=" * 72)
print(f"Landscape evaluation time   : {t_landscape:.2f} s (41 x 41 constant policies)")
print(f"Optimization time           : {t_opt:.2f} s")
print(f"Best constant policy        : x = {best_cx:.3f}, y = {best_cy:.3f}")
print("-" * 72)
print(f"{'Policy':<20s}{'Cumulative pop.':>18s}{'Final pop.':>14s}{'Min satisf.':>14s}")
for name, (tot, rec) in results.items():
    print(f"{name:<20s}{tot:>18,.0f}{rec['P'][-1]:>14,.0f}{rec['sat'].min():>14.3f}")
gain = 100.0 * (results["Optimized"][0] / results["Best constant"][0] - 1.0)
print("-" * 72)
print(f"Gain of optimized policy over best constant policy: {gain:+.2f} %")
print(f"Turns where the budget constraint was binding     : "
      f"{int(results['Optimized'][1]['bind'].sum())} / {T}")
print("-" * 72)
print("Optimized policy, turn by turn")
print(f"{'Turn':>4s}{'Population':>13s}{'Food stock':>13s}{'Immigrants':>13s}"
      f"{'Birth supp.':>13s}{'Satisf.':>10s}")
rec_opt = results["Optimized"][1]
for t in range(T):
    print(f"{t + 1:>4d}{rec_opt['P'][t + 1]:>13,.0f}{rec_opt['S'][t + 1]:>13,.0f}"
          f"{rec_opt['imm'][t]:>13,.1f}{rec_opt['bsp'][t]:>13,.1f}"
          f"{rec_opt['sat'][t]:>10.3f}")

# ----------------------------------------------------------------------
# Single combined figure
# ----------------------------------------------------------------------
plt.style.use("dark_background")
fig = plt.figure(figsize=(21, 12))
turns0 = np.arange(T + 1)
turns1 = np.arange(1, T + 1)

# (1) 3D landscape of constant policies
ax1 = fig.add_subplot(2, 3, 1, projection="3d")
Xg, Yg = np.meshgrid(grid, grid)
surf = ax1.plot_surface(Xg, Yg, Z / 1e3, cmap="viridis", alpha=0.92,
                        linewidth=0, antialiased=True)
ax1.scatter([best_cx], [best_cy], [best_const_val / 1e3], color="red", s=70,
            depthshade=False, label="Best constant policy")
ax1.set_xlabel("Immigration level x")
ax1.set_ylabel("Birth-support level y")
ax1.set_zlabel("Cumulative pop. (thousands)")
ax1.set_title("3D landscape of constant policies")
ax1.view_init(elev=28, azim=-125)
ax1.legend(loc="upper left")
fig.colorbar(surf, ax=ax1, shrink=0.55, pad=0.08)

# (2) Population trajectories
ax2 = fig.add_subplot(2, 3, 2)
for name, (tot, rec) in results.items():
    lw = 3.0 if name == "Optimized" else 1.8
    ax2.plot(turns0, rec["P"], label=name, color=colors[name], linewidth=lw)
ax2.set_xlabel("Turn")
ax2.set_ylabel("Population")
ax2.set_title("Population trajectories")
ax2.grid(alpha=0.25)
ax2.legend()

# (3) Optimal controls
ax3 = fig.add_subplot(2, 3, 3)
ax3.bar(turns1, rec_opt["imm"], color="#ff7043", alpha=0.85, label="Immigrants accepted")
ax3.set_xlabel("Turn")
ax3.set_ylabel("Immigrants per turn", color="#ff7043")
ax3b = ax3.twinx()
ax3b.plot(turns1, rec_opt["bsp"], color="#66bb6a", linewidth=2.5, marker="o",
          markersize=4, label="Birth-support spending")
ax3b.set_ylabel("Birth-support spending", color="#66bb6a")
ax3.set_title("Optimized controls (after budget scaling)")
h1, l1 = ax3.get_legend_handles_labels()
h2, l2 = ax3b.get_legend_handles_labels()
ax3.legend(h1 + h2, l1 + l2, loc="upper right")
ax3.grid(alpha=0.25)

# (4) Food stock and satisfaction under the optimized policy
ax4 = fig.add_subplot(2, 3, 4)
ax4.plot(turns0, rec_opt["S"], color="#ffd54f", linewidth=2.5, label="Food stock")
ax4.set_xlabel("Turn")
ax4.set_ylabel("Food stock", color="#ffd54f")
ax4b = ax4.twinx()
ax4b.plot(turns1, rec_opt["sat"], color="#4dd0e1", linewidth=2.5, marker="s",
          markersize=4, label="Satisfaction ratio")
ax4b.set_ylim(0.0, 1.05)
ax4b.set_ylabel("Satisfaction ratio", color="#4dd0e1")
ax4.set_title("Optimized policy: food stock and satisfaction")
h1, l1 = ax4.get_legend_handles_labels()
h2, l2 = ax4b.get_legend_handles_labels()
ax4.legend(h1 + h2, l1 + l2, loc="lower left")
ax4.grid(alpha=0.25)

# (5) 3D trajectories in (turn, population, food stock)
ax5 = fig.add_subplot(2, 3, 5, projection="3d")
for name, (tot, rec) in results.items():
    lw = 3.0 if name == "Optimized" else 1.8
    ax5.plot(turns0, rec["P"] / 1e3, rec["S"] / 1e3, color=colors[name],
             linewidth=lw, label=name)
ax5.set_xlabel("Turn")
ax5.set_ylabel("Population (thousands)")
ax5.set_zlabel("Food stock (thousands)")
ax5.set_title("3D state-space trajectories")
ax5.view_init(elev=24, azim=-60)
ax5.legend(loc="upper left")

# (6) Cumulative population by policy
ax6 = fig.add_subplot(2, 3, 6)
names = list(results.keys())
vals = [results[n][0] / 1e3 for n in names]
bars = ax6.bar(range(len(names)), vals, color=[colors[n] for n in names])
ax6.set_xticks(range(len(names)))
ax6.set_xticklabels(names, rotation=20, ha="right")
ax6.set_ylabel("Cumulative population (thousands)")
ax6.set_title("Objective value by policy")
for b, v in zip(bars, vals):
    ax6.text(b.get_x() + b.get_width() / 2.0, v, f"{v:,.1f}",
             ha="center", va="bottom", fontsize=10)
ax6.grid(axis="y", alpha=0.25)

fig.suptitle("Population Growth Control over 30 Turns", fontsize=18)
fig.tight_layout(rect=[0, 0, 1, 0.95])
plt.show()

Detailed Code Walkthrough

Parameters and helper

All constants sit at the top of the file and map one-to-one onto the table above. clip01 guarantees that any control passed to the simulator is a valid intensity in $[0, 1]$, so the optimizer can never push the model into an undefined region.

The simulator

simulate is the heart of the program. In each turn it performs these steps in order:

Production. It computes $Y_t = A(\lambda P_t)^{\beta}$ and the available food $S_t + Y_t$.
Satisfaction. If the available food covers the need $cP_t$, then $s_t = 1$ and the surplus is carried over. Otherwise the stock is emptied and $s_t$ equals the covered fraction.
Budget scaling. The requested spending $c_I I_t + B_t$ is compared with the tax revenue $\tau Y_t$. If it is too large, both requests are shrunk by the same factor $\kappa_t$. This is why requesting more than the treasury can afford is simply wasteful, never infeasible.
Demographics. Births use a saturating function of per-capita support, and deaths rise linearly with the shortage $1 - s_t$.
State update. The new population is floored at 1 to keep the divisions well defined, and the new stock is capped by spoilage and storage capacity.

When trace=True, the function also records every intermediate quantity for plotting. The optimization calls use trace=False, so no list operations slow the inner loop.

Landscape analysis

The double loop evaluates $41 \times 41 = 1681$ constant policies. Besides powering the 3D surface, it supplies the best constant policy, which is both a strong baseline and a safe starting point for the optimizer.

Full optimization

The decision vector z concatenates the 30 immigration intensities and the 30 birth-support intensities. L-BFGS-B handles the box constraints natively. Gradients come from finite differences with a step of $10^{-4}$, which is large enough to cope with the kinks created by min, max, and the budget scaling. Nine starting points (the best constant policy, three flat policies, and five random ones) reduce the risk of getting stuck in a poor local optimum. Since best_val is initialized with the best constant policy’s value, the final answer can never be worse than that baseline.

Local polish

The random local search perturbs only five randomly chosen variables at a time and accepts a candidate only if it improves $J$. The step size shrinks from 0.05 to 0.01 halfway through. This cheaply repairs the small imperfections that finite-difference gradients leave behind on a non-smooth objective.

Reporting and plotting

The console report prints the timing, the best constant policy, a comparison table, and a turn-by-turn log of the optimized policy. The figure is assembled from six panels in a single plt.show() call, so all graphs appear as one image.

Execution Results

Result Image

Console Output

========================================================================
Population growth control: results
========================================================================
Landscape evaluation time   : 0.08 s (41 x 41 constant policies)
Optimization time           : 11.94 s
Best constant policy        : x = 1.000, y = 0.625
------------------------------------------------------------------------
Policy                 Cumulative pop.    Final pop.   Min satisf.
No control                      32,441         1,161         1.000
Max immigration                 91,495         5,423         1.000
Max birth support               47,467         2,223         1.000
Best constant                   94,805         5,891         1.000
Optimized                       96,760         6,119         1.000
------------------------------------------------------------------------
Gain of optimized policy over best constant policy: +2.06 %
Turns where the budget constraint was binding     : 17 / 30
------------------------------------------------------------------------
Optimized policy, turn by turn
Turn   Population   Food stock   Immigrants  Birth supp.   Satisf.
   1        1,085        1,981         74.6          9.2     1.000
   2        1,174        2,445         78.9          9.9     1.000
   3        1,270        2,892         83.3         10.8     1.000
   4        1,370        3,000         87.9         11.5     1.000
   5        1,476        3,000         92.6         12.5     1.000
   6        1,588        3,000         97.4         13.5     1.000
   7        1,707        3,000        102.3         14.4     1.000
   8        1,831        3,000        107.4         15.6     1.000
   9        1,962        3,000        112.7         16.7     1.000
  10        2,100        3,000        118.1         17.8     1.000
  11        2,244        3,000        123.6         19.2     1.000
  12        2,396        3,000        129.2         20.6     1.000
  13        2,555        3,000        135.2         21.7     1.000
  14        2,722        3,000        141.3         23.0     1.000
  15        2,897        3,000        147.3         24.9     1.000
  16        3,079        3,000        150.0         33.7     1.000
  17        3,269        3,000        150.0         48.2     1.000
  18        3,467        3,000        150.0         63.2     1.000
  19        3,672        3,000        150.0         78.4     1.000
  20        3,885        3,000        150.0         94.0     1.000
  21        4,101        3,000        150.0        100.0     1.000
  22        4,319        2,997        150.0        100.0     1.000
  23        4,539        2,937        150.0        100.0     1.000
  24        4,760        2,819        150.0        100.0     1.000
  25        4,983        2,644        150.0        100.0     1.000
  26        5,207        2,412        150.0        100.0     1.000
  27        5,433        2,123        150.0        100.0     1.000
  28        5,660        1,779        150.0        100.0     1.000
  29        5,889        1,379        150.0        100.0     1.000
  30        6,119          924        150.0        100.0     1.000

Reading the Graphs

Panel 1: 3D landscape of constant policies

The surface height is the cumulative population achieved when the same immigration level $x$ and birth-support level $y$ are used in all 30 turns. Two forces shape it. Raising either control adds people, but pushing them too high drives the population past the carrying capacity, where shortage mortality eats into the gains. The budget scaling also flattens the surface where requests exceed what the treasury can pay for. The red marker shows the best constant policy. A broad ridge rather than a sharp peak means that the exact mix matters less than staying out of the overshoot region.

Panel 2: Population trajectories

This panel shows how each policy moves the population over time. “No control” shows pure natural dynamics. “Max immigration” and “Max birth support” show what happens when a single lever is pushed to its limit. A curve that rises steeply and then bends downward signals overshoot: the stock was drained and shortage deaths took over. The optimized curve is the one that keeps the area under the curve largest, since $J$ is exactly that area in discrete form.

Panel 3: Optimized controls

The bars show the immigrants actually accepted after budget scaling, and the line shows the actual birth-support spending. Because tax revenue grows with production, the affordable budget changes over time, and the shape of the optimal controls follows the economy’s capacity to pay. Where the two controls move in opposite directions, the planner is trading one lever for the other. The two levers differ in cost-effectiveness: an immigrant costs a fixed amount and arrives immediately, while birth support has diminishing returns per capita.

Panel 4: Food stock and satisfaction

This panel exposes the resource side of the trade-off. The food stock acts as a buffer that lets the population temporarily exceed what current production can feed. A satisfaction ratio of 1 means no one is short of food. Any dip below 1 means shortage mortality is active. Seeing exactly when the planner lets the buffer run down, and how close the satisfaction ratio gets to the danger zone, explains why the optimal policy is a controlled use of the stock rather than a strict avoidance of shortage.

Panel 5: 3D state-space trajectories

Here each policy is drawn as a curve in $(\text{turn},\ P_t,\ S_t)$ space. A healthy trajectory climbs in population while the food stock stays positive. A trajectory that heads toward the floor of the food-stock axis while the population is still high is living off its reserves and is exposed to a sudden shortage. Comparing the optimized curve with the extreme policies shows at a glance how the planner balances population growth against reserve depletion.

Panel 6: Objective value by policy

The bars compare the final objective $J$ for all five policies. The gap between the best constant policy and the optimized policy measures the value of adapting the controls over time rather than holding them fixed. The gaps to the single-lever policies show how costly it is to ignore one of the two instruments.

Discussion

Three structural features of the model explain the behavior of the optimal policy.

Diminishing returns on labor. Because $\beta < 1$, doubling the population does not double the food supply. Per-capita food falls as the population grows, which creates a natural carrying capacity and makes unlimited growth self-defeating.

A buffer with a cost. The food stock lets the society run a temporary deficit, but spoilage makes holding a large stock expensive. The optimizer therefore treats the stock as a resource to be spent wisely rather than a quantity to be maximized.

A finite horizon. The objective counts only the first 30 turns. Consequences that would arrive after the horizon do not affect $J$, so the optimal policy can legitimately behave differently near the final turns than in the middle of the horizon. This end-of-horizon effect is a general property of finite-horizon optimal control and is worth keeping in mind when interpreting the optimized controls.

The model is intentionally simple, yet it already contains the essential tension of the problem: growth is valuable in itself, but unchecked growth erodes the resource base that supports it. Extensions are straightforward. Introducing age structure with a maturation delay would make birth support slower to pay off, a stochastic harvest would turn the stock into a true insurance, and a closed-loop controller that reacts to the observed state could outperform the open-loop schedule found here.

Conclusion

We formulated population growth control as a 60-variable constrained optimization problem and solved it with a landscape scan, multi-start L-BFGS-B, and a sparse local search, all built on a fast scalar simulator. The 3D landscape revealed the trade-off between the two levers, and the 3D state-space trajectories showed how the optimized policy uses the food buffer. The key lesson is that the best way to maximize cumulative population is neither to grow as fast as possible nor to grow cautiously, but to steer the growth rate so that population, workforce, and food supply rise together.

How to Maximize a Nation's Population 50 Years from Now

October 27, 2026

Investing in Health vs. Machines

Every government faces the same awkward question: do we spend on hospitals today, or on factories?

Medical spending saves lives, but the money cannot also be used to build the capital that produces tomorrow’s output. At the same time, the population is aging. Children, workers and the elderly all change in size, and workers pay for everything.

In this article I build a small simulation of this dilemma and solve it numerically in Python. The goal is to maximize the total population 50 years from now, and to find out what the best investment path looks like.

1. The Model

Time is split into 10 periods of 5 years each, so $t = 0, 1, \dots, 9$ covers 50 years. The population has three age groups: children $C_t$, workers $W_t$ and the elderly $E_t$.

Production

Output follows a Cobb–Douglas function of capital $K_t$ and workers $W_t$:

$$
Y_t = A , K_t^{\alpha} , W_t^{1-\alpha}
$$

A fixed share $s$ of output is invested each period, so the investment budget is $I_t = s,Y_t,\Delta t$. The decision variable $u_t \in [0,1]$ is the share of that budget that goes to medical investment. The remainder goes to productive capital:

$$
K_{t+1} = (1-\delta_K),K_t + (1-u_t),I_t
$$

$$
H_{t+1} = (1-\delta_H),H_t + \frac{u_t,I_t}{h}
$$

Here $H_t$ is the accumulated health stock, which depreciates at rate $\delta_H$, and $h$ is a scale constant.

Mortality

The baseline mortality rates $m_C, m_W, m_E$ are scaled by a common multiplier:

$$
\mu_t = \Big[,m_{\min} + (1-m_{\min}),e^{-H_t}\Big]\cdot\Big(\frac{y_t}{y_0}\Big)^{-\theta},
\qquad y_t = \frac{Y_t}{N_t}
$$

The first factor says that medical investment lowers mortality, with diminishing returns and a floor $m_{\min}$. The second factor says that higher income per capita, which is an effect of the production investment, also lowers mortality through nutrition and sanitation. This second term is what creates a real trade-off. Without it, putting everything into health would always win.

Population dynamics

With birth rate $f$ and aging rates $g_{CW}$ (child to worker) and $g_{WE}$ (worker to elderly):

$$
C_{t+1} = C_t,(1 - g_{CW} - \mu_t m_C) + f,W_t
$$

$$
W_{t+1} = W_t,(1 - g_{WE} - \mu_t m_W) + g_{CW},C_t
$$

$$
E_{t+1} = E_t,(1 - \mu_t m_E) + g_{WE},W_t
$$

The optimization problem

$$
\max_{u_0,\dots,u_9 \in [0,1]} ; N_{10} = C_{10} + W_{10} + E_{10}
$$

Example setup

Item	Value
Initial population (millions)	children 25, workers 60, elderly 15
Initial capital $K_0$	150
Baseline mortality per 5 years $(m_C, m_W, m_E)$	0.05, 0.10, 0.45
Capital share $\alpha$, saving rate $s$	0.35, 0.25
Birth rate $f$	0.11 per worker per period

With no health investment at all, this society shrinks. That is the baseline we want to beat.

2. The Code

The simulator is fully vectorized. It evaluates hundreds of investment plans in a single call, so the 3D surface (1,681 scenarios) is computed in a fraction of a second and no extra speed-up is needed. The optimization uses SLSQP from 20 random starting points to avoid poor local solutions.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.optimize import minimize

plt.rcParams["figure.dpi"] = 110
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.3

# ------------------------------------------------------------
# 1. Parameters (one period = 5 years, 10 periods = 50 years)
# ------------------------------------------------------------
T, DT = 10, 5
ALPHA, A = 0.35, 1.0
S_RATE = 0.25
D_K, D_H = 0.20, 0.15
F_BIRTH = 0.11
G_CW, G_WE = 1 / 4, 1 / 9
M0 = np.array([0.05, 0.10, 0.45])
M_FLOOR = 0.15
H_SCALE = 60.0
THETA = 0.8
C0, W0, E0, K0, H0 = 25.0, 60.0, 15.0, 150.0, 0.0
Y_INIT = A * K0 ** ALPHA * W0 ** (1 - ALPHA)
N_INIT = C0 + W0 + E0

# ------------------------------------------------------------
# 2. Vectorized simulator: U has shape (batch, T)
# ------------------------------------------------------------
def simulate(U):
    U = np.atleast_2d(np.asarray(U, dtype=float))
    B = U.shape[0]
    C = np.full(B, C0); W = np.full(B, W0); E = np.full(B, E0)
    K = np.full(B, K0); H = np.full(B, H0)
    out = np.zeros((T + 1, 6, B))
    out[0] = [C, W, E, K, H, A * K ** ALPHA * W ** (1 - ALPHA)]
    for t in range(T):
        Y = A * K ** ALPHA * W ** (1 - ALPHA)
        I = S_RATE * Y * DT
        N = C + W + E
        ratio = np.clip((Y / N) / (Y_INIT / N_INIT), 0.5, 2.0)
        mult = (M_FLOOR + (1 - M_FLOOR) * np.exp(-H)) * ratio ** (-THETA)
        mc, mw, me = M0[0] * mult, M0[1] * mult, M0[2] * mult
        u = U[:, t]
        Cn = C * (1 - G_CW - mc) + F_BIRTH * W
        Wn = W * (1 - G_WE - mw) + G_CW * C
        En = E * (1 - me) + G_WE * W
        Kn = (1 - D_K) * K + (1 - u) * I
        Hn = (1 - D_H) * H + u * I / H_SCALE
        C, W, E, K, H = Cn, Wn, En, Kn, Hn
        out[t + 1] = [C, W, E, K, H, A * K ** ALPHA * W ** (1 - ALPHA)]
    return out

def final_population(U):
    return simulate(U)[-1, :3].sum(axis=0)

# ------------------------------------------------------------
# 3. Optimization (the control of the last period has no effect
#    on the 50-year population, so only T-1 controls are optimized)
# ------------------------------------------------------------
def objective(x):
    return -final_population(np.append(x, x[-1]))[0]

rng = np.random.default_rng(42)
best = None
for _ in range(20):
    x0 = rng.uniform(0, 1, T - 1)
    res = minimize(objective, x0, bounds=[(0, 1)] * (T - 1), method="SLSQP")
    if best is None or res.fun < best.fun:
        best = res
u_opt = np.append(best.x, best.x[-1])

# Best constant share
grid_u = np.linspace(0, 1, 101)
N_const = final_population(np.tile(grid_u[:, None], (1, T)))
u_const = grid_u[np.argmax(N_const)]

strategies = {
    "No medical investment (u=0)": np.zeros(T),
    "All-in medical (u=1)": np.ones(T),
    "Balanced (u=0.5)": np.full(T, 0.5),
    f"Best constant (u={u_const:.2f})": np.full(T, u_const),
    "Optimal time-varying": u_opt,
}
results = {name: simulate(u)[:, :, 0] for name, u in strategies.items()}

# ------------------------------------------------------------
# 4. Console report
# ------------------------------------------------------------
years = np.arange(T + 1) * DT
print("=== Optimal medical investment share u_t (per 5-year period) ===")
for t in range(T):
    print(f"  Year {years[t]:2d}-{years[t+1]:2d}: u = {u_opt[t]:.3f}")
print()
print("=== Population after 50 years (millions) ===")
header = f"{'Strategy':34s}{'Children':>10s}{'Workers':>10s}{'Elderly':>10s}{'Total':>10s}{'Output Y':>10s}"
print(header)
print("-" * len(header))
for name, r in results.items():
    c, w, e, k, h, y = r[-1]
    print(f"{name:34s}{c:10.1f}{w:10.1f}{e:10.1f}{c+w+e:10.1f}{y:10.1f}")
print()
r = results["Optimal time-varying"]
print("=== Age composition under the optimal policy ===")
for t in [0, 5, 10]:
    c, w, e = r[t, :3]
    n = c + w + e
    print(f"  Year {years[t]:2d}: children {c/n:5.1%} | workers {w/n:5.1%} | elderly {e/n:5.1%} | total {n:6.1f}")
gain = results["Optimal time-varying"][-1, :3].sum() / results["No medical investment (u=0)"][-1, :3].sum() - 1
print()
print(f"Gain of the optimal policy over no medical investment: {gain:+.1%}")

# ------------------------------------------------------------
# 5. Data for the 3D surfaces
# ------------------------------------------------------------
ug = np.linspace(0, 1, 41)
UE, UL = np.meshgrid(ug, ug)
Ugrid = np.concatenate(
    [np.repeat(UE.reshape(-1, 1), 5, axis=1), np.repeat(UL.reshape(-1, 1), 5, axis=1)], axis=1
)
Z = final_population(Ugrid).reshape(UE.shape)
iz = np.unravel_index(np.argmax(Z), Z.shape)

us = np.linspace(0, 1, 41)
traj = simulate(np.tile(us[:, None], (1, T)))[:, :3, :].sum(axis=1)  # (T+1, 41)
TT, UU = np.meshgrid(years, us)

# ------------------------------------------------------------
# 6. One combined figure
# ------------------------------------------------------------
fig = plt.figure(figsize=(21, 12))
col = {"C": "#4C9BE8", "W": "#F2A541", "E": "#7A5C9E"}

ax1 = fig.add_subplot(2, 3, 1)
ax1.step(years[:-1], u_opt, where="post", color="#D1495B", lw=2.5)
ax1.set_ylim(0, 1.05)
ax1.set_xlabel("Year"); ax1.set_ylabel("Medical share u (optimal)", color="#D1495B")
ax2 = ax1.twinx()
ax2.plot(years, r[:, 3], color="#2E86AB", lw=2, label="Capital K")
ax2.plot(years, r[:, 4] * 50, color="#3BB273", lw=2, ls="--", label="Health stock H (x50)")
ax2.set_ylabel("Stock"); ax2.legend(loc="upper right"); ax2.grid(False)
ax1.set_title("Optimal investment path and stocks")

ax3 = fig.add_subplot(2, 3, 2)
ax3.stackplot(years, r[:, 0], r[:, 1], r[:, 2],
              labels=["Children", "Workers", "Elderly"],
              colors=[col["C"], col["W"], col["E"]], alpha=0.9)
ax3.set_xlabel("Year"); ax3.set_ylabel("Population (millions)")
ax3.set_title("Age structure under the optimal policy"); ax3.legend(loc="upper left")

ax4 = fig.add_subplot(2, 3, 3)
styles = ["#999999", "#E4572E", "#17BEBB", "#76B041", "#111111"]
for (name, res), c in zip(results.items(), styles):
    ax4.plot(years, res[:, :3].sum(axis=1), lw=3 if "Optimal" in name else 2, color=c, label=name)
ax4.set_xlabel("Year"); ax4.set_ylabel("Total population (millions)")
ax4.set_title("Total population by strategy"); ax4.legend(fontsize=8)

ax5 = fig.add_subplot(2, 3, 4, projection="3d")
surf = ax5.plot_surface(UE, UL, Z, cmap=cm.viridis, alpha=0.9, linewidth=0)
ax5.scatter([UE[iz]], [UL[iz]], [Z[iz]], color="red", s=80, depthshade=False)
ax5.set_xlabel("u (years 0-25)"); ax5.set_ylabel("u (years 25-50)"); ax5.set_zlabel("Population at year 50")
ax5.set_title(f"Early vs late investment (max = {Z[iz]:.1f})")
ax5.view_init(elev=28, azim=-125)
fig.colorbar(surf, ax=ax5, shrink=0.55, pad=0.1)

ax6 = fig.add_subplot(2, 3, 5, projection="3d")
surf2 = ax6.plot_surface(TT, UU, traj.T, cmap=cm.plasma, alpha=0.9, linewidth=0)
ax6.set_xlabel("Year"); ax6.set_ylabel("Constant medical share u"); ax6.set_zlabel("Total population")
ax6.set_title("Population over time for each constant share")
ax6.view_init(elev=26, azim=-130)
fig.colorbar(surf2, ax=ax6, shrink=0.55, pad=0.1)

ax7 = fig.add_subplot(2, 3, 6)
labels = ["Year 0", "Year 25", "Year 50\n(no medical)", "Year 50\n(optimal)"]
snap = [r[0, :3], r[5, :3], results["No medical investment (u=0)"][-1, :3], r[-1, :3]]
shares = np.array([s / s.sum() for s in snap]) * 100
bottom = np.zeros(len(snap))
for i, (nm, c) in enumerate(zip(["Children", "Workers", "Elderly"], [col["C"], col["W"], col["E"]])):
    ax7.bar(labels, shares[:, i], bottom=bottom, color=c, label=nm)
    for j, v in enumerate(shares[:, i]):
        ax7.text(j, bottom[j] + v / 2, f"{v:.0f}%", ha="center", va="center", color="white", fontsize=10)
    bottom += shares[:, i]
ax7.set_ylim(0, 118)
ax7.set_ylabel("Share of population (%)"); ax7.set_title("Demographic shares"); ax7.legend(loc="upper center", ncol=3, fontsize=9)

plt.tight_layout()
plt.show()

3. Code Walkthrough

Parameters. All constants are gathered at the top. Time steps are 5 years, so aging rates are per-period fractions: children spend about 20 years as children ($1/4$ per period) and workers about 45 years as workers ($1/9$ per period). Y_INIT and N_INIT store initial output and population, which are needed as the reference point for income per capita.

simulate(U). This is the core of the program. It takes a matrix U of shape (batch, 10), where every row is one complete investment plan. All state variables are NumPy arrays of length batch, so a single loop over 10 periods advances every plan simultaneously. Inside the loop, the code does the following:

Computes output $Y_t$ and the investment budget $I_t$.
Computes the mortality multiplier $\mu_t$ from the health stock and from income per capita. The ratio is clipped to $[0.5, 2]$ so that extreme values cannot make mortality explode or vanish.
Updates the three age groups with the transition equations.
Splits $I_t$ between capital $K$ and health $H$ according to $u_t$.

The result has shape $(11, 6, \text{batch})$: 11 time points, six variables ($C, W, E, K, H, Y$) and the batch dimension.

final_population. It returns the sum of the three age groups at year 50. This is the objective function.

Optimization. SLSQP is a gradient-based method that handles the bounds $0 \le u_t \le 1$. Because the landscape could in principle have several local optima, I run it from 20 random starts and keep the best. The control of the last period cannot influence the population at year 50 (its effect would only appear in the next period), so I optimize only nine variables and repeat the last value for display purposes.

Benchmarks. Four reference strategies are compared with the optimum: no medical spending, all-in medical spending, a balanced 50/50 split, and the best constant share found by scanning 101 values.

3D data. For the first surface, I split the 50 years into two halves and let $u$ differ between them: UE for the first 25 years and UL for the last 25 years. The 41 × 41 grid is flattened into 1,681 plans and evaluated in one vectorized call. For the second surface, I evaluate 41 constant shares over the whole time axis.

Figure. Seven panels are drawn on one canvas, including two 3D surfaces.

4. Results

Console output

=== Optimal medical investment share u_t (per 5-year period) ===
  Year  0- 5: u = 1.000
  Year  5-10: u = 0.763
  Year 10-15: u = 0.500
  Year 15-20: u = 0.485
  Year 20-25: u = 0.475
  Year 25-30: u = 0.470
  Year 30-35: u = 0.468
  Year 35-40: u = 0.468
  Year 40-45: u = 0.470
  Year 45-50: u = 0.470

=== Population after 50 years (millions) ===
Strategy                            Children   Workers   Elderly     Total  Output Y
------------------------------------------------------------------------------------
No medical investment (u=0)             14.1      29.1      12.5      55.8      78.7
All-in medical (u=1)                    19.9      43.6      35.9      99.4      30.8
Balanced (u=0.5)                        20.1      44.4      40.2     104.7      80.4
Best constant (u=0.64)                  20.3      44.8      40.7     105.8      70.6
Optimal time-varying                    20.5      45.3      41.6     107.4      80.3

=== Age composition under the optimal policy ===
  Year  0: children 25.0% | workers 60.0% | elderly 15.0% | total  100.0
  Year 25: children 22.2% | workers 48.8% | elderly 29.1% | total  100.5
  Year 50: children 19.1% | workers 42.2% | elderly 38.7% | total  107.4

Gain of the optimal policy over no medical investment: +92.6%

Graph output

5. Interpreting the Results

Optimal investment path (top left). The best policy is front-loaded. In the first 5 years it puts everything into medicine ($u \approx 1$), then tapers to about 0.76, and from year 10 onward it settles at a steady share of roughly 0.47. The logic is that early health investment has diminishing returns. The first dollars cut mortality sharply, but once the health stock is built, further spending buys less. At that point the money is better used for capital, which raises income and keeps mortality low through the second channel. Notice that capital dips at first because of the sacrifice, then recovers and grows steadily.

Age structure (top middle). Even under the optimal policy, the population ages. The elderly grow steadily while the worker group shrinks. This is the cost of saving lives: the people we keep alive become elderly, and births alone cannot offset the shift. In the optimal run, the elderly make up 15% of the population at the start and nearly 39% at year 50.

Strategy comparison (top right). Without medical investment the population falls from 100 to about 56 million. Every health-oriented strategy reverses this decline. The notable point is the all-in strategy. It starts well but flattens and even declines late, because starving capital reduces income per capita, which pushes mortality back up. The final output $Y$ makes this vivid. All-in medicine ends with an output of about 31, compared with about 80 for the optimal plan, so the economy has been hollowed out. The optimal time-varying path beats even the best constant share, though only slightly. The gap between the best constant plan and the optimum is small, which shows that the level of investment matters more than its exact timing, while front-loading adds a modest extra gain.

3D surface: early vs. late (bottom left). The horizontal axes are the medical shares in the first and second halves of the period, and the height is the population at year 50. The surface is a smooth hill with an interior peak (marked in red). Corners are poor: spending nothing early or late is clearly bad, and so is going to 100% in both halves. This confirms that the optimum is a genuine balance, not a corner solution. The ridge also shows that a high early share paired with a moderate late share does best, consistent with front-loading.

3D surface: constant share over time (bottom middle). This plot shows how total population evolves for each constant $u$. At $u=0$ the surface slopes sharply downward over time. As $u$ increases, the decline turns into growth. For very high $u$, growth flattens at later dates, which reflects the capital starvation described above. The ridge of the surface lies at intermediate values of $u$.

Demographic shares (bottom right). Without medical investment, the age mix at year 50 looks healthy on paper: 52% workers and only 22% elderly. But that “healthy” mix comes from a population that has shrunk by almost half because people died. The optimal policy produces a much larger population, but with only 42% workers and 39% elderly. This is the central tension in the problem: saving lives and keeping a young age structure pull in opposite directions.

6. Takeaways

Medical investment sharply raises the long-run population (about +93% versus doing nothing in this example), even though it reduces short-term production.
The optimum is an interior solution, not “all medicine”: the productive economy has to be kept alive, because it funds the very health system we are building.
Front-loading is worthwhile. Build the health stock early, then shift part of the budget back to capital.
A larger population is an older population. A policy that maximizes headcount will push the elderly share up, so pension and care systems need to be planned alongside medical investment.

The model is intentionally simple. A natural next step is to add a fertility response to income, age-specific health returns, or a utility function that values workers more than the elderly. Each of these would change the optimal path, and the same code structure can handle them.

Optimal Food Storage Under Spoilage and Storage Costs

October 26, 2026

Eat Now or Store for Later?

Every village, colony, and game-economy designer eventually faces the same dilemma. A harvest has just come in. Do you feed everyone generously today, which makes the population grow, or do you lock part of the grain away for a future that might be hungrier? Storage is not free, and some of the grain will rot. In this article we turn that dilemma into a precise optimization problem, solve it in Python with an exact gradient (so it runs fast), and visualize the answer in a single dashboard that includes 3D plots.

1. The Problem

A settlement lives for $T = 20$ turns. In every turn:

A harvest $H_t$ arrives. It varies with a seasonal cycle on top of a slowly growing trend.
The settlement decides how much of its available food to eat now and how much to store.
Eaten food determines population growth. Stored food loses a fraction to spoilage and costs a fraction in upkeep.

The objective is to maximize the total population summed over all 20 turns.

The tension is easy to state:

Eating more now grows the population sooner, and a larger population counts in every later turn too.
But a larger population also needs more food later.
Eating too much in a bumper-harvest turn wastes food, because growth saturates.
Starving in a lean turn is punished heavily, because growth collapses.
Stored food is a buffer against lean turns, but it decays.

2. Mathematical Model

State and decision. In turn $t = 0, 1, \dots, T-1$ the state is the population $P_t$ and the food stock $S_t$. The decision is the eating fraction $u_t \in [0, 1]$.

Food balance.

$$
A_t = S_t + H_t, \qquad c_t = u_t A_t, \qquad r_t = \frac{c_t}{P_t}
$$

Here $A_t$ is the available food, $c_t$ is the food eaten, and $r_t$ is the per-capita ration, where $r_t = 1$ means exactly one unit of food per person.

Population growth. The growth factor depends on the ration through a smooth saturating curve:

$$
P_{t+1} = P_t , g(r_t), \qquad g(r) = g_{\mathrm{lo}} + \left(g_{\mathrm{hi}} - g_{\mathrm{lo}}\right)\left(1 - e^{-(r/r_c)^2}\right)
$$

With no food the population shrinks by the factor $g_{\mathrm{lo}}$. With plentiful food it grows by at most $g_{\mathrm{hi}}$. The squared exponent makes the curve S-shaped, so mild shortages are tolerable but severe shortages are disastrous, and extra food beyond the knee of the curve brings little benefit.

Storage with spoilage and cost. Each stored unit independently spoils with probability $\delta$ per turn. For a large stock, the law of large numbers makes the expected stock an excellent description, so spoilage removes the fraction $\delta$ of what is stored. Storage upkeep costs $\kappa$ food units per stored unit per turn. Writing $\rho = 1 - \delta - \kappa$ for the survival factor:

$$
S_{t+1} = \rho , (1 - u_t) , A_t
$$

Objective.

$$
\max_{u_0, \dots, u_{T-1} \in [0,1]} ; J = \sum_{t=0}^{T-1} P_{t+1}
$$

Food left in storage after the final turn is worthless, so the model automatically pushes the settlement to run its stock down toward the end.

3. Concrete Example

Quantity	Symbol	Value
Number of turns	$T$	20
Initial population	$P_0$	100
Initial stock	$S_0$	40
Growth factor with no food	$g_{\mathrm{lo}}$	0.70
Maximum growth factor	$g_{\mathrm{hi}}$	1.12
Ration scale	$r_c$	1.0
Harvest	$H_t$	$(110 + 5t)\left(1 + a \sin\frac{2\pi t}{5}\right)$
Seasonal amplitude	$a$	0.5
Spoilage probability	$\delta$	0.15
Storage cost	$\kappa$	0.03

The harvest grows linearly and oscillates with a period of five turns, so there are rich turns and lean turns. We compare three policies:

Optimal: the best eating fractions found by the optimizer.
Eat-all: consume everything every turn ($u_t = 1$).
Constant 80%: always eat 80% of the available food and store the rest.

We also repeat the optimization over a grid of spoilage rates $\delta$ and seasonal amplitudes $a$ to see when storage is worth it.

4. Solution Method and Speed-Up

The problem has 20 continuous decision variables with simple bounds, and it is not convex. The growth curve is S-shaped, so a single starting point may end in a poor local optimum. We therefore use bounded L-BFGS-B from several different starting policies and keep the best result.

The naive approach lets the optimizer estimate gradients by finite differences, which costs 21 full simulations per gradient. We instead compute the exact gradient with a backward (adjoint) recursion. Let $J_t = \sum_{k=t}^{T-1} P_{k+1}$ be the population total from turn $t$ onward, and define the sensitivities

$$
\mu_t = \frac{\partial J_t}{\partial P_t}, \qquad \nu_t = \frac{\partial J_t}{\partial S_t}, \qquad \mu_T = \nu_T = 0 .
$$

With $a_t = 1 + \mu_{t+1}$ and $g’(r) = (g_{\mathrm{hi}} - g_{\mathrm{lo}}), e^{-(r/r_c)^2}, \dfrac{2r}{r_c^2}$, the recursion runs from $t = T-1$ down to $t = 0$:

$$
\frac{\partial J}{\partial u_t} = a_t, g’(r_t), A_t - \nu_{t+1}, \rho, A_t
$$

$$
\mu_t = a_t \left( g(r_t) - r_t, g’(r_t) \right), \qquad
\nu_t = a_t, g’(r_t), u_t + \nu_{t+1}, \rho, (1 - u_t)
$$

One forward pass plus one backward pass gives the whole gradient. That is roughly twenty times cheaper than finite differences, which makes the 121-point parameter sweep in Section 7 finish quickly.

5. Full Source Code

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
from scipy.optimize import minimize

# ----------------------------------------------------------------------
# Problem parameters
# ----------------------------------------------------------------------
T = 20
P0 = 100.0
S0 = 40.0
G_LO, G_HI, R_C = 0.70, 1.12, 1.0
BASE, TREND, PERIOD = 110.0, 5.0, 5.0
AMP0, DELTA0, KAPPA = 0.5, 0.15, 0.03


def harvest(amp):
    t = np.arange(T)
    return (BASE + TREND * t) * (1.0 + amp * np.sin(2.0 * np.pi * t / PERIOD))


def growth(r):
    return G_LO + (G_HI - G_LO) * (1.0 - np.exp(-(r / R_C) ** 2))


def growth_with_slope(r):
    e = math.exp(-(r / R_C) ** 2)
    g = G_LO + (G_HI - G_LO) * (1.0 - e)
    dg = (G_HI - G_LO) * e * 2.0 * r / R_C ** 2
    return g, dg


# ----------------------------------------------------------------------
# Forward simulation of a given policy u (eating fraction per turn)
# ----------------------------------------------------------------------
def simulate(u, H, rho):
    P = np.empty(T + 1)
    S = np.empty(T + 1)
    R = np.empty(T)
    C = np.empty(T)
    P[0] = P0
    S[0] = S0
    for t in range(T):
        A = S[t] + H[t]
        c = u[t] * A
        r = c / P[t]
        g, _ = growth_with_slope(r)
        P[t + 1] = P[t] * g
        S[t + 1] = rho * (A - c)
        R[t] = r
        C[t] = c
    return P, S, R, C


# ----------------------------------------------------------------------
# Objective and exact gradient (forward pass + adjoint backward pass)
# ----------------------------------------------------------------------
def neg_value_and_grad(u, H, rho):
    P = np.empty(T + 1)
    S = np.empty(T + 1)
    G = np.empty(T)
    Gd = np.empty(T)
    A = np.empty(T)
    R = np.empty(T)
    P[0] = P0
    S[0] = S0
    for t in range(T):
        A[t] = S[t] + H[t]
        R[t] = u[t] * A[t] / P[t]
        G[t], Gd[t] = growth_with_slope(R[t])
        P[t + 1] = P[t] * G[t]
        S[t + 1] = rho * (1.0 - u[t]) * A[t]
    J = P[1:].sum()

    grad = np.empty(T)
    wP = 0.0
    wS = 0.0
    for t in range(T - 1, -1, -1):
        aP = 1.0 + wP
        grad[t] = aP * Gd[t] * A[t] - wS * rho * A[t]
        wP_new = aP * (G[t] - Gd[t] * R[t])
        wS_new = aP * Gd[t] * u[t] + wS * rho * (1.0 - u[t])
        wP, wS = wP_new, wS_new
    return -J, -grad


# ----------------------------------------------------------------------
# Multi-start bounded optimization
# ----------------------------------------------------------------------
STARTS = [np.full(T, 1.0), np.full(T, 0.8), np.full(T, 0.5)]


def solve(H, rho, starts):
    best = None
    for u_init in starts:
        res = minimize(
            neg_value_and_grad,
            u_init,
            args=(H, rho),
            jac=True,
            method="L-BFGS-B",
            bounds=[(0.0, 1.0)] * T,
            options={"maxiter": 500, "ftol": 1e-12, "gtol": 1e-9},
        )
        if best is None or res.fun < best.fun:
            best = res
    return best


# ----------------------------------------------------------------------
# Base case: optimal policy versus two simple baselines
# ----------------------------------------------------------------------
H_base = harvest(AMP0)
rho0 = 1.0 - DELTA0 - KAPPA

res0 = solve(H_base, rho0, STARTS)
u_opt = res0.x
P_opt, S_opt, R_opt, C_opt = simulate(u_opt, H_base, rho0)
J_opt = P_opt[1:].sum()

P_eat, S_eat, R_eat, C_eat = simulate(np.ones(T), H_base, rho0)
J_eat = P_eat[1:].sum()

P_con, S_con, R_con, C_con = simulate(np.full(T, 0.8), H_base, rho0)
J_con = P_con[1:].sum()

print("=" * 78)
print(
    "Optimal policy  (spoilage {:.0%}, storage cost {:.0%}, "
    "seasonal amplitude {:.1f})".format(DELTA0, KAPPA, AMP0)
)
print("=" * 78)
print(
    f"{'Turn':>4} {'Harvest':>9} {'Stock':>9} {'Eaten':>9} "
    f"{'Eat frac':>9} {'Ration':>7} {'Pop(end)':>9}"
)
for t in range(T):
    print(
        f"{t + 1:>4d} {H_base[t]:>9.1f} {S_opt[t]:>9.1f} {C_opt[t]:>9.1f} "
        f"{u_opt[t]:>9.3f} {R_opt[t]:>7.3f} {P_opt[t + 1]:>9.1f}"
    )
print("-" * 78)
print(f"Total population over {T} turns")
print(f"  Optimal policy : {J_opt:>10.1f}")
print(f"  Eat-all policy : {J_eat:>10.1f}")
print(f"  Constant 80%   : {J_con:>10.1f}")
print(f"  Gain of optimal over eat-all: {100.0 * (J_opt / J_eat - 1.0):.2f} %")

# ----------------------------------------------------------------------
# Parameter sweep: spoilage rate x seasonal amplitude
# ----------------------------------------------------------------------
deltas = np.linspace(0.0, 0.5, 11)
amps = np.linspace(0.0, 0.8, 11)
Z_opt = np.zeros((len(amps), len(deltas)))
Z_eat = np.zeros((len(amps), len(deltas)))

for i, amp in enumerate(amps):
    H_i = harvest(amp)
    for j, delta in enumerate(deltas):
        rho_ij = 1.0 - delta - KAPPA
        res = solve(H_i, rho_ij, STARTS)
        Z_opt[i, j] = -res.fun
        P_e, _, _, _ = simulate(np.ones(T), H_i, rho_ij)
        Z_eat[i, j] = P_e[1:].sum()

Z_gain = 100.0 * (Z_opt - Z_eat) / Z_eat
imax, jmax = np.unravel_index(np.argmax(Z_gain), Z_gain.shape)
print("-" * 78)
print(
    f"Sweep: largest gain over eat-all = {Z_gain[imax, jmax]:.2f} % "
    f"(spoilage {deltas[jmax]:.0%}, amplitude {amps[imax]:.1f})"
)
print(f"Sweep: smallest gain over eat-all = {Z_gain.min():.2f} %")

# ----------------------------------------------------------------------
# One combined dashboard
# ----------------------------------------------------------------------
plt.style.use("dark_background")
C_OPT, C_EAT, C_CON = "#00e5ff", "#ff6e6e", "#ffd166"

fig = plt.figure(figsize=(22, 12))
fig.suptitle(
    "Eat Now or Store for Later?  Optimal Food Policy over 20 Turns",
    fontsize=18,
    fontweight="bold",
)

turns0 = np.arange(T + 1)
turns1 = np.arange(1, T + 1)

# (1) Population trajectories
ax1 = fig.add_subplot(2, 3, 1)
ax1.plot(turns0, P_opt, color=C_OPT, lw=2.5, label=f"Optimal (sum {J_opt:.0f})")
ax1.plot(turns0, P_eat, color=C_EAT, lw=2, ls="--", label=f"Eat-all (sum {J_eat:.0f})")
ax1.plot(turns0, P_con, color=C_CON, lw=2, ls=":", label=f"Constant 80% (sum {J_con:.0f})")
ax1.set_title("Population")
ax1.set_xlabel("Turn")
ax1.set_ylabel("Population")
ax1.grid(alpha=0.25)
ax1.legend()

# (2) Harvest, consumption, stock under the optimal policy
ax2 = fig.add_subplot(2, 3, 2)
ax2.bar(turns1, H_base, color="#4a90d9", alpha=0.45, label="Harvest")
ax2.plot(turns1, C_opt, "o-", color=C_OPT, lw=2, label="Food eaten")
ax2.plot(turns1, S_opt[:T], "s--", color=C_CON, lw=2, label="Stock at start of turn")
ax2.set_title("Optimal policy: harvest, consumption, stock")
ax2.set_xlabel("Turn")
ax2.set_ylabel("Food units")
ax2.grid(alpha=0.25)
ax2.legend()

# (3) 3D trajectories: turn - population - stock
ax3 = fig.add_subplot(2, 3, 3, projection="3d")
ax3.plot(turns0, P_opt, S_opt, color=C_OPT, lw=2.5, label="Optimal")
ax3.plot(turns0, P_eat, S_eat, color=C_EAT, lw=2, label="Eat-all")
ax3.plot(turns0, P_con, S_con, color=C_CON, lw=2, label="Constant 80%")
ax3.set_title("3D trajectories")
ax3.set_xlabel("Turn")
ax3.set_ylabel("Population")
ax3.set_zlabel("Stock")
ax3.view_init(elev=24, azim=-60)
ax3.legend()

# (4) 3D surface: optimal total population
X, Y = np.meshgrid(deltas, amps)
ax4 = fig.add_subplot(2, 3, 4, projection="3d")
surf4 = ax4.plot_surface(X, Y, Z_opt, cmap="viridis", edgecolor="none", alpha=0.95)
ax4.scatter([DELTA0], [AMP0], [J_opt], color="red", s=70, depthshade=False)
ax4.set_title("Optimal total population")
ax4.set_xlabel("Spoilage rate")
ax4.set_ylabel("Seasonal amplitude")
ax4.set_zlabel("Total population")
ax4.view_init(elev=28, azim=-125)
fig.colorbar(surf4, ax=ax4, shrink=0.55, pad=0.08)

# (5) 3D surface: percentage gain over eat-all
ax5 = fig.add_subplot(2, 3, 5, projection="3d")
surf5 = ax5.plot_surface(X, Y, Z_gain, cmap="plasma", edgecolor="none", alpha=0.95)
ax5.set_title("Value of storing: gain over eat-all (%)")
ax5.set_xlabel("Spoilage rate")
ax5.set_ylabel("Seasonal amplitude")
ax5.set_zlabel("Gain (%)")
ax5.view_init(elev=28, azim=-125)
fig.colorbar(surf5, ax=ax5, shrink=0.55, pad=0.08)

# (6) Growth curve with the operating points of the optimal policy
ax6 = fig.add_subplot(2, 3, 6)
r_max = max(2.5, float(R_opt.max()) * 1.1)
rr = np.linspace(0.0, r_max, 300)
ax6.plot(rr, growth(rr), color="white", lw=2)
sc = ax6.scatter(
    R_opt, growth(R_opt), c=np.arange(1, T + 1), cmap="plasma", s=60, zorder=3
)
ax6.axhline(1.0, color="gray", ls="--")
ax6.axvline(1.0, color="gray", ls="--")
ax6.set_title("Growth factor g(r) and optimal operating points")
ax6.set_xlabel("Per-capita ration r")
ax6.set_ylabel("Growth factor")
ax6.grid(alpha=0.25)
fig.colorbar(sc, ax=ax6, label="Turn")

fig.subplots_adjust(
    left=0.05, right=0.97, bottom=0.07, top=0.90, wspace=0.28, hspace=0.30
)
plt.show()

6. Code Walkthrough

Parameters and helper functions. The constants at the top mirror the table in Section 3. harvest(amp) builds the harvest sequence, a linear trend multiplied by a sine wave whose size is controlled by amp. growth(r) is the vectorized growth curve used for plotting. growth_with_slope(r) is the scalar version that returns both $g(r)$ and $g’(r)$, which the gradient needs. It uses math.exp, which is faster than NumPy for single numbers and quietly returns zero for very negative arguments, so extreme rations cannot cause an overflow.

simulate. This is a literal translation of the model. For each turn it computes the available food, the amount eaten, and the ration, then updates the population with the growth factor and the stock with the survival factor $\rho$. It also returns the rations and consumption so we can plot them later.

neg_value_and_grad. This is the heart of the speed-up. The first loop is the same simulation, but it stores the intermediate growth factors, their slopes, the available food, and the rations. The second loop walks backward from the last turn and carries the two sensitivities wP and wS, the code versions of $\mu$ and $\nu$. At each turn it writes the exact derivative of the total population with respect to that turn’s eating fraction. The function returns the negated value and gradient because SciPy minimizes.

solve. Because the growth curve is S-shaped, the problem can have several local optima. The function therefore runs L-BFGS-B from three starting policies (eat everything, eat 80%, eat half) and keeps the best result. All eating fractions are bounded to $[0, 1]$ through the bounds argument. Passing jac=True tells SciPy that our function returns the value and gradient together, so no finite differences are needed.

Base case. The script solves the 20-turn problem for the base parameters, then simulates the two baseline policies for comparison. It prints a turn-by-turn table (harvest, stock, food eaten, eating fraction, ration, and ending population) and the three totals.

Parameter sweep. Two nested loops solve the full optimization on an 11 × 11 grid of spoilage rates (0% to 50%) and seasonal amplitudes (0 to 0.8). For every grid point we record the optimal total and the eat-all total, then convert the difference into a percentage gain. The exact gradient is what keeps this loop fast: 121 optimizations, each with three starts, finish in a short time.

Dashboard. All six panels are drawn in a single figure so that one image captures the whole story. The 3D panels use projection="3d", and view_init sets the camera angle so that the surfaces are easy to read.

7. Results

Console Output

==============================================================================
Optimal policy  (spoilage 15%, storage cost 3%, seasonal amplitude 0.5)
==============================================================================
Turn   Harvest     Stock     Eaten  Eat frac  Ration  Pop(end)
   1     110.0      40.0     136.6     0.910   1.366     105.5
   2     169.7      11.0     130.4     0.722   1.237     108.5
   3     155.3      41.2     121.3     0.617   1.117     108.5
   4      88.3      61.7     110.4     0.736   1.018     105.3
   5      68.2      32.4     100.6     1.000   0.955     100.2
   6     135.0       0.0     135.0     1.000   1.347     105.4
   7     206.6       0.0     144.6     0.700   1.372     111.3
   8     187.6      50.8     139.9     0.587   1.257     115.0
   9     105.9      80.8     133.0     0.712   1.156     116.1
  10      81.3      44.1     125.4     1.000   1.080     114.9
  11     160.0       0.0     160.0     1.000   1.393     121.7
  12     243.5       0.0     170.5     0.700   1.401     129.1
  13     220.0      59.8     164.9     0.589   1.277     134.0
  14     123.6      94.2     155.9     0.716   1.164     135.6
  15      94.4      50.7     145.1     1.000   1.071     133.7
  16     185.0       0.0     185.0     1.000   1.383     141.5
  17     280.4       0.0     232.7     0.830   1.645     154.5
  18     252.3      39.1     224.4     0.770   1.452     165.2
  19     141.2      54.9     196.2     1.000   1.188     168.1
  20     107.5       0.0     107.5     1.000   0.640     141.3
------------------------------------------------------------------------------
Total population over 20 turns
  Optimal policy :     2515.4
  Eat-all policy :     2329.4
  Constant 80%   :     2303.2
  Gain of optimal over eat-all: 7.98 %
------------------------------------------------------------------------------
Sweep: largest gain over eat-all = 39.83 % (spoilage 0%, amplitude 0.8)
Sweep: smallest gain over eat-all = 0.00 %

Result Image

8. How to Read the Results

Population panel (top left). The three curves show how the same village evolves under three different philosophies. The legend reports each policy’s 20-turn sum, which is exactly the quantity we maximized. The gap between the cyan optimal curve and the red eat-all curve is the value created purely by timing the consumption. Whatever the exact numbers, the optimal curve cannot lie below the eat-all curve in total, because eat-all is one of the starting points of the search.

Harvest, consumption, and stock panel (top middle). The blue bars are the harvest, the cyan line is the food actually eaten, and the yellow dashed line is the stock at the beginning of each turn. Where the cyan line is smoother than the bars, the settlement is turning a spiky harvest into steady meals by moving food through the warehouse. The stock line should rise after rich harvests and fall during lean turns. Near the final turns it should shrink toward zero, because leftover grain at the end of turn 20 earns nothing.

3D trajectory panel (top right). Each policy is a path through the space of turn, population, and stock. The eat-all path hugs the floor of the stock axis because it never stores anything. The optimal path lifts off the floor whenever it builds a buffer and returns to it as the horizon approaches. Rotating the view in a notebook makes it easy to see how the buffer and the population interact.

Optimal total population surface (bottom left). This surface shows the best achievable total as a function of the spoilage rate and the seasonal amplitude. Higher spoilage means a leakier warehouse, which lowers what an optimal planner can achieve. The red dot marks the base case from Section 3 so the single run can be located on the surface.

Value of storing surface (bottom middle). This is the most decision-relevant picture. It plots the percentage improvement of the optimal policy over eat-all. The surface answers the practical question of when storage is worth the trouble. Storage should matter most where seasons swing widely and spoilage is low, since that is where moving food in time pays the most. It should matter least where the harvest is steady or spoilage is severe. The sweep summary lines in the console print the largest and smallest gains on the grid.

Growth curve panel (bottom right). The white curve is the S-shaped growth factor $g(r)$, and the colored dots are the rations the optimal policy actually chooses, colored by turn. Dots clustered near the knee of the curve reveal the central idea of the solution: because the curve flattens at high rations and falls steeply at low rations, an optimal planner avoids both overeating and starvation and keeps the ration in the productive middle zone. Any dots that stray far from that zone typically appear at the very end of the horizon, where the planner spends the last of its stock.

9. Takeaways

The model is small, but it captures the core logic of storage decisions:

Smoothing is valuable because the growth response is nonlinear. Because starvation hurts more than abundance helps, shifting food from rich turns to lean turns raises the total even when part of the stored food is lost.
Spoilage and upkeep set the price of smoothing. Each unit moved through the warehouse is discounted by $\rho = 1 - \delta - \kappa$, so storage is worthwhile only when the benefit of smoothing exceeds that loss.
The horizon matters. With a finite horizon, the stock has no value at the end, so the optimal policy ends by consuming everything.
Exact gradients make the exploration cheap. The adjoint recursion turned a 121-point parameter sweep into a task that finishes quickly, which lets us study the whole landscape instead of a single case.

Change $\delta$, $\kappa$, the harvest cycle, or the shape of the growth curve and rerun the script. Each change tells a different story about when it pays to eat now and when it pays to wait.

How Much Should a Region Spend on Disaster Preparedness?

1. Problem Setup

2. Source Code

3. Code Walkthrough

Parameters (Section 1)

Scenario generation (Section 2)

The simulator (Section 3)

Parameter decoding (decode)

Optimization (Section 4)

Strategy comparison (Section 5)

3D surface scan (Section 6)

4. Execution Results

5. Reading the Results

Optimal allocation

Strategy comparison

Figure guide

6. Conclusion

How Much Should a Nation Invest in Irrigation and Grain Reserves to Maximize Its Population in 30 Years?

1. The Model

2. Full Source Code

3. Code Walkthrough

Section 1: Parameters

Section 2: The Simulator

Section 3: Grid Search

Section 4: Out-of-Sample Evaluation

Section 5: Graphs

4. Results

Console Output

Graphs

5. Reading the Results

The comparison table

Graph (1): Rainfall distribution

Graph (2): 3D surface of expected population

Graph (3): Heat map

Graph (4): Population paths

Graph (5): Distribution of final population

Graph (6): Irrigation coverage and reserve over time

6. Summary

Choosing Which Roads, Power Lines, and Water Mains to Build

The Problem

Concrete Example

Mathematical Formulation

Source Code

Code Walkthrough

1. Problem Data

2. Vectorized Exhaustive Search

3. Greedy Baseline

4. Console Report

5. Sensitivity Data

6. Combined Figure

Execution Results

Interpreting the Results

Optimal Plan vs. Greedy

Panel (a): Network Map

Panel (b): 3D Capacity by City and Type

Panel (c): Budget vs. Capacity

Panel (d): 3D Surface of Budget × Grid Weight

Summary

Maximizing Total Population with Python

1. Problem Setting

Population model

Optimality conditions

2. Example Problem

3. Source Code

4. Code Walkthrough

Parameters and the population model

Numerical optimization

Brute-force verification

Marginal returns

Comparison strategies

5. Results

Result images

Console output

6. Interpreting the Graphs

7. Conclusion

Maximizing the 10-Year Total with Python

The Problem

A Concrete Example

The Mathematical Model

Source Code

Parameter decoding (`decode`)