Managing Renewable and Finite Resources Over 100 Turns with Python

November 7, 2026

Planning Long-Term Civilization Survival

A growing population produces more, but it also consumes more. When some resources regenerate and others run out for good, and when droughts, earthquakes, and swings in demand arrive at random, the question becomes: how should a civilization plan so that its population after 100 turns is as large as possible, without ever dying out along the way?

In this article we build a concrete model, turn it into an optimization problem, and solve it with a vectorized evolutionary algorithm in Python.

1. The Example Problem

The civilization starts with a population $P_0 = 100$, a renewable stock $R_0 = 800$ (forests, fisheries, fertile soil), a finite stock $F_0 = 2000$ (ore, fossil fuel), and a resource efficiency $E_0 = 1$.

Each turn $t = 0, 1, \dots, 99$ the planner chooses three controls:

$u_t \in [0, 0.30]$: the fraction of the renewable stock harvested
$v_t \in [0, 0.06]$: the fraction of the finite stock mined
$w_t \in [0, 0.50]$: the share of output invested in resource efficiency

To keep the search space small, the controls are piecewise constant over 10 blocks of 10 turns, which gives $3 \times 10 = 30$ decision variables.

Three random shocks hit each scenario:

Drought: with probability 0.05 per turn, the renewable stock is multiplied by $0.55$.
Earthquake: with probability 0.03 per turn, the population is multiplied by $0.80$ and efficiency is cut by 10%.
Demand fluctuation: per-capita demand is multiplied by $\eta_t = 1 + 0.15\sin(2\pi t/25 + \phi) + 0.05,\xi_t$, where $\xi_t \sim \mathcal{N}(0,1)$.

2. Mathematical Model

Let $\delta_t$ and $q_t$ be the drought and earthquake multipliers (equal to 1 when no disaster occurs). After shocks, $R_t’ = \delta_t R_t$ and $P_t’ = q_t P_t$.

Renewable resource (logistic regrowth minus harvest):

Finite resource:

$$
M_t = v_t F_t, \qquad F_{t+1} = F_t - M_t
$$

Production (labor and resources, with diminishing returns):

$$
Y_t = A,(P_t’)^{\alpha},\left(H_t + \gamma M_t\right)^{1-\alpha}
$$

Demand and supply ratio:

$$
D_t = \frac{d_0,\eta_t,P_t’}{E_t}, \qquad s_t = \frac{(1 - w_t),Y_t}{D_t}
$$

Population growth (decline is faster than growth):

$$
g_t = \begin{cases} 0.05\tanh\bigl(3(s_t - 1)\bigr) & s_t \ge 1 \ 0.12\tanh\bigl(3(s_t - 1)\bigr) & s_t < 1 \end{cases}, \qquad P_{t+1} = P_t’,(1 + g_t)
$$

If $P_{t+1} < 1$, the civilization is extinct and the population is set to 0 permanently.

Efficiency growth:

The parameters are $r = 0.12$, $K = 1000$, $A = 1.35$, $\alpha = 0.6$, $\gamma = 1.5$, $d_0 = 1$, $\kappa = 0.05$, and $E_{\max} = 4$.

Optimization problem. Over $S$ random scenarios, with the controls $\boldsymbol{\theta} = (u_t, v_t, w_t)$ as decision variables, we solve

with $P_{\mathrm{safe}} = 30$ and $\lambda = 20$. The first term maximizes the expected final population. The second term is a penalty that forces the population to stay away from extinction in every scenario, not just on average.

3. Solution Strategy

The model is nonlinear, has hard constraints (extinction is absorbing), and contains random shocks, so gradient-based methods are a poor fit. We use differential evolution, a population-based global optimizer.

The simulation of 100 turns is cheap, but an evolutionary algorithm needs tens of thousands of evaluations. The code below is therefore written to be fast:

The simulator is vectorized over candidates and scenarios simultaneously, so one call evaluates a whole generation of 360 candidate plans on 40 scenarios with only 100 loop iterations.
SciPy’s vectorized=True and updating="deferred" options pass the entire population to the objective function at once.
Plans are optimized on 40 training scenarios and then validated on 400 unseen scenarios, so the reported performance is out-of-sample.

4. Source Code

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

# ----------------------------------------------------------------------
# 1. Problem setting
# ----------------------------------------------------------------------
T, NB = 100, 10
BL = T // NB
P0, R0, F0, E0 = 100.0, 800.0, 2000.0, 1.0
K, REGEN = 1000.0, 0.12
A, ALPHA, GAMMA, D0 = 1.35, 0.6, 1.5, 1.0
U_MAX, V_MAX, W_MAX = 0.30, 0.06, 0.50
E_CAP, E_GROW = 4.0, 0.05
P_SAFE, PENALTY = 30.0, 20.0

def make_scenarios(seed, n):
    rng = np.random.default_rng(seed)
    t = np.arange(T)[None, :]
    phase = rng.uniform(0.0, 2.0 * np.pi, size=(n, 1))
    return {
        "drought": np.where(rng.random((n, T)) < 0.05, 0.55, 1.0),
        "quake": np.where(rng.random((n, T)) < 0.03, 0.80, 1.0),
        "demand": 1.0 + 0.15 * np.sin(2.0 * np.pi * t / 25.0 + phase)
                  + 0.05 * rng.standard_normal((n, T)),
    }

# ----------------------------------------------------------------------
# 2. Vectorized simulator: (candidates M) x (scenarios S) at once
# ----------------------------------------------------------------------
def simulate(theta, scen, record=False):
    theta = np.atleast_2d(theta)
    M = theta.shape[0]
    S = scen["drought"].shape[0]
    u = np.repeat(theta[:, 0:NB] * U_MAX, BL, axis=1)
    v = np.repeat(theta[:, NB:2 * NB] * V_MAX, BL, axis=1)
    w = np.repeat(theta[:, 2 * NB:3 * NB] * W_MAX, BL, axis=1)
    P = np.full((M, S), P0)
    R = np.full((M, S), R0)
    F = np.full((M, S), F0)
    E = np.full((M, S), E0)
    min_p = P.copy()
    if record:
        hP = np.zeros((T + 1, M, S))
        hR = np.zeros((T + 1, M, S))
        hF = np.zeros((T + 1, M, S))
        hP[0], hR[0], hF[0] = P, R, F
    for t in range(T):
        dr = scen["drought"][:, t]
        qk = scen["quake"][:, t]
        dm = scen["demand"][:, t]
        R = R * dr
        P = P * qk
        E = E * (1.0 - 0.5 * (1.0 - qk))
        ut = u[:, t:t + 1]
        vt = v[:, t:t + 1]
        wt = w[:, t:t + 1]
        H = ut * R
        Mn = vt * F
        X = H + GAMMA * Mn
        Y = A * np.power(P, ALPHA) * np.power(X, 1.0 - ALPHA)
        D = np.maximum(D0 * dm * P / E, 1e-9)
        s = (1.0 - wt) * Y / D
        g = np.where(s >= 1.0, 0.05, 0.12) * np.tanh(3.0 * (s - 1.0))
        R = np.clip(R + REGEN * R * (1.0 - R / K) - H, 0.0, None)
        F = F - Mn
        E = np.minimum(E * (1.0 + E_GROW * wt), E_CAP)
        P = P * (1.0 + g)
        P = np.where(P < 1.0, 0.0, P)
        min_p = np.minimum(min_p, P)
        if record:
            hP[t + 1], hR[t + 1], hF[t + 1] = P, R, F
    if record:
        return P, min_p, hP, hR, hF
    return P, min_p

def objective(x, scen):
    p_fin, min_p = simulate(x.T, scen)
    shortfall = np.maximum(0.0, P_SAFE - min_p).sum(axis=1)
    return -p_fin.mean(axis=1) + PENALTY * shortfall

# ----------------------------------------------------------------------
# 3. Optimization with differential evolution
# ----------------------------------------------------------------------
train = make_scenarios(2024, 40)
test = make_scenarios(777, 400)

t0 = time.time()
result = differential_evolution(
    objective, bounds=[(0.0, 1.0)] * (3 * NB), args=(train,),
    vectorized=True, updating="deferred", maxiter=300, popsize=12,
    tol=1e-9, seed=1, polish=False, mutation=(0.5, 1.0), recombination=0.8,
)
elapsed = time.time() - t0
theta_opt = result.x

theta_greedy = np.concatenate([np.ones(NB), np.ones(NB), np.zeros(NB)])
theta_cons = np.concatenate([np.full(NB, 0.2), np.full(NB, 0.15), np.full(NB, 0.2)])
names = ["Greedy", "Conservative", "Optimized"]
colors = ["#ff5c5c", "#ffd166", "#06d6a0"]
thetas = np.vstack([theta_greedy, theta_cons, theta_opt])

p_fin, min_p, hP, hR, hF = simulate(thetas, test, record=True)

# ----------------------------------------------------------------------
# 4. Console report
# ----------------------------------------------------------------------
print(f"Optimization time: {elapsed:.1f} s  ({result.nit} generations)")
print(f"Best training objective: {-result.fun:.2f}")
print()
print(f"{'Policy':<14}{'Mean P(100)':>13}{'Median':>10}{'5th pct':>10}{'Extinct %':>12}{'Min<30 %':>11}")
for k, name in enumerate(names):
    print(f"{name:<14}{p_fin[k].mean():>13.1f}{np.median(p_fin[k]):>10.1f}"
          f"{np.percentile(p_fin[k], 5):>10.1f}{100 * np.mean(min_p[k] < 1.0):>12.1f}"
          f"{100 * np.mean(min_p[k] < P_SAFE):>11.1f}")
print()
print("Optimized schedule per 10-turn block")
print(f"{'Block':<8}{'Harvest u':>11}{'Mining v':>11}{'Invest w':>11}")
for b in range(NB):
    print(f"{b * BL + 1:>3}-{(b + 1) * BL:<4}{theta_opt[b] * U_MAX:>11.3f}"
          f"{theta_opt[NB + b] * V_MAX:>11.4f}{theta_opt[2 * NB + b] * W_MAX:>11.3f}")

# ----------------------------------------------------------------------
# 5. Policy landscape for the 3D surface
# ----------------------------------------------------------------------
G = 26
gu = np.linspace(0.0, 1.0, G)
gv = np.linspace(0.0, 1.0, G)
GU, GV = np.meshgrid(gu, gv)
grid_theta = np.zeros((G * G, 3 * NB))
grid_theta[:, 0:NB] = GU.ravel()[:, None]
grid_theta[:, NB:2 * NB] = GV.ravel()[:, None]
grid_theta[:, 2 * NB:] = theta_opt[None, 2 * NB:]
grid_p, _ = simulate(grid_theta, train)
Z = grid_p.mean(axis=1).reshape(G, G)
best = np.unravel_index(np.argmax(Z), Z.shape)

# ----------------------------------------------------------------------
# 6. One combined figure
# ----------------------------------------------------------------------
plt.style.use("dark_background")
fig = plt.figure(figsize=(22, 13))
gs = GridSpec(2, 3, figure=fig)
turns = np.arange(T + 1)

ax1 = fig.add_subplot(gs[0, 0])
for k in range(3):
    lo, med, hi = np.percentile(hP[:, k, :], [10, 50, 90], axis=1)
    ax1.fill_between(turns, lo, hi, color=colors[k], alpha=0.18)
    ax1.plot(turns, med, color=colors[k], lw=2.2, label=names[k])
ax1.set_title("Population (median and 10-90% band)")
ax1.set_xlabel("Turn")
ax1.set_ylabel("Population")
ax1.legend(loc="upper left")
ax1.grid(alpha=0.2)

ax2 = fig.add_subplot(gs[0, 1])
rm = np.percentile(hR[:, 2, :], [10, 50, 90], axis=1)
fm = np.percentile(hF[:, 2, :], [10, 50, 90], axis=1)
ax2.fill_between(turns, rm[0], rm[2], color="#4cc9f0", alpha=0.2)
ax2.plot(turns, rm[1], color="#4cc9f0", lw=2.2, label="Renewable stock R")
ax2.set_xlabel("Turn")
ax2.set_ylabel("Renewable stock", color="#4cc9f0")
ax2.grid(alpha=0.2)
ax2b = ax2.twinx()
ax2b.plot(turns, fm[1], color="#f72585", lw=2.2, label="Finite stock F")
ax2b.set_ylabel("Finite stock", color="#f72585")
ax2.set_title("Resource stocks under the optimized plan")
h1, l1 = ax2.get_legend_handles_labels()
h2, l2 = ax2b.get_legend_handles_labels()
ax2.legend(h1 + h2, l1 + l2, loc="center right")

ax3 = fig.add_subplot(gs[0, 2])
tt = np.arange(T)
ax3.step(tt, np.repeat(theta_opt[0:NB], BL), where="post", color="#06d6a0", lw=2.2,
         label=f"Harvest rate u (max {U_MAX:.2f})")
ax3.step(tt, np.repeat(theta_opt[NB:2 * NB], BL), where="post", color="#f72585", lw=2.2,
         label=f"Mining rate v (max {V_MAX:.2f})")
ax3.step(tt, np.repeat(theta_opt[2 * NB:], BL), where="post", color="#ffd166", lw=2.2,
         label=f"Efficiency investment w (max {W_MAX:.2f})")
ax3.set_title("Optimized control schedule (fraction of maximum)")
ax3.set_xlabel("Turn")
ax3.set_ylabel("Normalized control")
ax3.set_ylim(-0.05, 1.05)
ax3.legend(loc="center left", fontsize=9)
ax3.grid(alpha=0.2)

ax4 = fig.add_subplot(gs[1, 0], projection="3d")
ax4.plot_surface(GU * U_MAX, GV * V_MAX, Z, cmap="viridis", edgecolor="none", alpha=0.92)
ax4.scatter([gu[best[1]] * U_MAX], [gv[best[0]] * V_MAX], [Z[best]], color="#ff5c5c", s=80,
            depthshade=False)
ax4.set_title("Mean final population vs constant (u, v)")
ax4.set_xlabel("Harvest rate u")
ax4.set_ylabel("Mining rate v")
ax4.set_zlabel("Mean P(100)")
ax4.view_init(elev=28, azim=-125)

ax5 = fig.add_subplot(gs[1, 1], projection="3d")
for k in range(3):
    ax5.plot(hP[:, k, :].mean(axis=1), hR[:, k, :].mean(axis=1), hF[:, k, :].mean(axis=1),
             color=colors[k], lw=2.4, label=names[k])
    ax5.scatter([hP[0, k, :].mean()], [hR[0, k, :].mean()], [hF[0, k, :].mean()],
                color="white", s=30, depthshade=False)
ax5.set_title("Mean trajectory in (P, R, F) space")
ax5.set_xlabel("Population")
ax5.set_ylabel("Renewable R")
ax5.set_zlabel("Finite F")
ax5.legend(loc="upper left")
ax5.view_init(elev=22, azim=-60)

ax6 = fig.add_subplot(gs[1, 2])
bp = ax6.boxplot([p_fin[k] for k in range(3)], patch_artist=True, showfliers=False)
for patch, c in zip(bp["boxes"], colors):
    patch.set_facecolor(c)
    patch.set_alpha(0.6)
for med in bp["medians"]:
    med.set_color("white")
ax6.set_xticks([1, 2, 3])
ax6.set_xticklabels(names)
for k in range(3):
    ax6.text(k + 1, ax6.get_ylim()[1] * 0.95,
             f"extinct {100 * np.mean(min_p[k] < 1.0):.1f}%", ha="center", color=colors[k])
ax6.set_title("Final population over 400 unseen scenarios")
ax6.set_ylabel("Population at turn 100")
ax6.grid(alpha=0.2, axis="y")

fig.suptitle("Long-term civilization survival plan", fontsize=18)
plt.tight_layout(rect=[0, 0, 1, 0.97])
plt.show()

5. Code Walkthrough

Section 1: problem setting. All model constants are defined once at the top. make_scenarios draws, for every scenario, a drought multiplier array, an earthquake multiplier array, and a demand multiplier array, each of shape (n, T). A fixed random seed makes the experiment reproducible. Two independent scenario sets are created later: 40 for training and 400 for testing.

Section 2: the vectorized simulator. simulate is the heart of the program. The decision vector theta has 30 entries in $[0,1]$: the first 10 are harvest intensities, the next 10 mining intensities, and the last 10 investment shares. They are scaled by U_MAX, V_MAX, and W_MAX and expanded to one value per turn with np.repeat.

State variables P, R, F, and E are arrays of shape (M, S): $M$ candidate plans times $S$ scenarios. The time loop applies the shocks, harvests and mines, computes production $Y$, demand $D$, and the supply ratio $s$, and then updates every stock. Because a single loop advances all candidates and all scenarios at once, the cost per generation is nearly independent of population size. The line P = np.where(P < 1.0, 0.0, P) makes extinction absorbing, and min_p tracks the lowest population ever reached in each run. With record=True the full histories are kept for plotting.

Section 3: the objective and the optimizer. objective receives an array of shape (D, M), which is the format SciPy uses when vectorized=True, so it is transposed before simulating. It returns the negative mean final population plus the survival penalty, because SciPy minimizes. Differential evolution uses popsize=12, so each generation contains $12 \times 30 = 360$ candidate plans, and it runs for 300 generations. updating="deferred" is required for vectorized evaluation. The whole search takes only seconds, whereas a plain Python loop over every candidate and scenario would be orders of magnitude slower.

Baselines. Two hand-made policies serve as a reference. The greedy policy always harvests and mines at the maximum rate and never invests in efficiency. The conservative policy harvests and mines at low rates and invests a small constant share.

Section 4: the console report. The optimized plan and both baselines are evaluated on the 400 unseen scenarios. The table reports the mean, median, and 5th percentile of the final population, the extinction rate, and the share of scenarios in which the population dipped below the safety level of 30. The optimized control schedule is printed block by block.

Section 5: the policy landscape. To visualize how sensitive the outcome is to the two resource-extraction controls, we evaluate a $26 \times 26$ grid of constant harvest rates $u$ and mining rates $v$, holding the optimized investment schedule fixed. Each grid point is one more row in the same vectorized call.

Section 6: the combined figure. All six panels are drawn into one figure. The panels show the population bands for the three policies, the resource stocks under the optimized plan, the control schedule, the 3D landscape of final population against $(u, v)$, the 3D trajectory in $(P, R, F)$ space, and the box plots of final population on unseen scenarios.

6. Execution Results

Console Output

Optimization time: 35.6 s  (300 generations)
Best training objective: 793.55

Policy          Mean P(100)    Median   5th pct   Extinct %   Min<30 %
Greedy                  0.2       0.0       0.0        82.8      100.0
Conservative           64.1      60.1      24.2         0.0       20.0
Optimized             846.8     833.7     381.8         0.0        5.8

Optimized schedule per 10-turn block
Block     Harvest u   Mining v   Invest w
  1-10        0.098     0.0029      0.499
 11-20        0.043     0.0246      0.492
 21-30        0.051     0.0218      0.493
 31-40        0.055     0.0251      0.493
 41-50        0.042     0.0278      0.436
 51-60        0.045     0.0267      0.353
 61-70        0.031     0.0315      0.272
 71-80        0.042     0.0333      0.152
 81-90        0.038     0.0566      0.000
 91-100       0.153     0.0597      0.001

Graph Output

7. Interpreting the Results

Population trajectories. The greedy plan looks attractive at first: harvesting and mining at full speed produces a short boom in population. But the renewable stock collapses below the level where it can regenerate, the finite stock is burned through, and the population falls toward zero. In the large majority of unseen scenarios the greedy civilization goes extinct. The conservative plan never collapses, but it never grows either: the population stagnates at a low level because nothing is invested in efficiency, and a noticeable share of scenarios still dip below the safety level after a disaster. The optimized plan first dips slightly, then grows slowly, and finally accelerates sharply. Its final population is many times larger than the conservative one, with no extinctions in the test set.

Resource stocks. Under the optimized plan the renewable stock is quickly drawn down from its initial value to a moderate level and then held steady in a band that allows strong regeneration, close to the maximum-sustainable-yield point of the logistic curve. This is the signature of sustainable management: harvest the flow, not the stock. The finite stock declines smoothly and is deliberately spent over the whole horizon, ending near exhaustion at turn 100.

Control schedule. The schedule has a clear three-phase structure.

Investment phase (early turns). Mining is almost off, harvesting is modest, and the investment share $w_t$ sits at its upper bound. Output is sacrificed to build efficiency $E_t$, which lowers the per-capita demand $D_t$ for the rest of the game.
Steady phase (middle turns). Harvesting stays near the sustainable level, mining rises gradually, and investment is phased out as efficiency approaches its cap.
Drawdown phase (final turns). With the horizon about to end, there is no reason to preserve stock for the future. Investment drops to zero, mining reaches its maximum, and harvesting jumps. This is the classic end-of-horizon effect in finite-horizon optimal control.

3D landscape. The surface of mean final population against constant $(u, v)$ is far from flat. Pushing the harvest rate too high creates a deep valley: the renewable stock is destroyed and the civilization collapses. Pushing it too low wastes production. The good region is a narrow ridge, and the optimizer has to find it in a 30-dimensional space, where an exhaustive grid search would be hopeless.

3D trajectory. In $(P, R, F)$ space, the greedy trajectory dives toward the zero-population axis after a short excursion. The conservative trajectory barely moves in population. The optimized trajectory first loses renewable stock to a sustainable level, then climbs steadily in population while the finite stock is consumed, and ends at a high population.

Robustness on unseen scenarios. The box plots confirm that the result is not an artifact of the 40 training scenarios. On 400 fresh scenarios the optimized plan keeps a clearly higher median final population than the alternatives, with a zero extinction rate. The penalty term $\lambda \sum_s \max(0, P_{\mathrm{safe}} - \min_t P_t^{(s)})$ is what buys this safety: it makes plans that gamble on avoiding disasters unattractive, even when their average outcome would look good.

8. Conclusion

We turned a vague question about civilization survival into a concrete optimization problem with renewable and finite resources, population-dependent production, and random disasters. A vectorized simulator combined with differential evolution found a plan that is both ambitious (a final population far above the cautious baseline) and safe (no extinction on unseen scenarios). The structure of the solution is intuitive: invest early in efficiency, harvest the renewable resource at its sustainable rate, spend the finite resource steadily, and draw everything down only when the horizon ends.

The same framework extends naturally to richer models: more resource types, age-structured populations, finer control blocks, or risk measures such as conditional value-at-risk in place of the simple survival penalty.

Optimizing Equipment Investment to Maximize Cumulative Population Over 20 Turns

November 6, 2026

Automation vs. Scarcity

Automation is a classic lever for growth. Every machine you add lets each worker produce more food, which supports a larger population. But machines are not free. They need materials to build, energy to run, and they wear out. Meanwhile, the food spent on construction is food your people cannot eat today.

In this article we turn that tension into a concrete optimization problem, solve it with a vectorized evolutionary search plus local polishing, and visualize the result with a nine-panel dark-themed figure that includes two 3D plots.

1. The Problem

A settlement starts with 100 people. Each turn (think of it as one year) the government decides how many automation units $u_t$ to build. The goal is to maximize the total population summed over 20 turns:

$$
\max_{u_0,\dots,u_{19}} ;; \sum_{t=1}^{20} P_t \qquad \text{subject to } 0 \le u_t \le 12 .
$$

Model

Workers and effective equipment. A fixed fraction $\eta$ of the population works. The energy supply $E$ can power only $K_{\max} = E/e$ units of equipment, so any equipment above that ceiling sits idle:

$$
L_t = \eta P_t, \qquad K_t^{\mathrm{eff}} = \min!\left(K_t,; K_{\max}\right), \qquad K_{\max} = \frac{E}{e}.
$$

Food production. Land is fixed, so labor has diminishing returns ($\beta<1$). Automation multiplies output, with a saturating effect that depends on equipment per worker:

$$
Y_t = A,L_t^{\beta}\left(1+\gamma,\frac{K_t^{\mathrm{eff}}}{K_t^{\mathrm{eff}}+\kappa L_t}\right).
$$

Construction is limited by materials and food. Building one unit costs $m$ units of material and $c$ units of food. The number actually built is capped by what is available:

$$
I_t=\min!\left(u_t,;\frac{M_t}{m},;\frac{Y_t}{2c}\right), \qquad F_t = Y_t - c,I_t .
$$

Population dynamics. Let $\nu$ be the food needed per person. The ratio $q_t = F_t/(\nu P_t)$ drives growth. Surplus food ($q_t>1$) raises the population, a shortage ($q_t<1$) shrinks it, and the rate is bounded by $r$:

$$
P_{t+1}=P_t\left[1+r\tanh!\big(2,(q_t-1)\big)\right].
$$

Stocks. Equipment depreciates at rate $\delta$. Materials are consumed by construction and partially replenished by recycling at a constant rate $\rho$:

$$
K_{t+1}=(1-\delta)K_t+I_t, \qquad M_{t+1}=M_t-m,I_t+\rho .
$$

Parameters

Symbol	Meaning	Value
$P_0$	Initial population	100
$\eta$	Working fraction	0.6
$\nu$	Food need per person	1.0
$A,\ \beta$	Land productivity, labor exponent	9.5, 0.6
$\gamma,\ \kappa$	Automation gain, saturation scale	1.2, 0.5
$E,\ e$	Energy supply per turn, energy per unit	30, 0.5 (so $K_{\max}=60$)
$K_0$	Initial equipment	5
$M_0,\ \rho,\ m$	Initial materials, recycling per turn, material per unit	400, 4, 6
$\delta$	Depreciation	0.06
$c$	Food cost per unit built	2
$r$	Maximum growth rate per turn	0.12

The material budget is the crux. The stock of 400 plus 4 per turn of recycling supports roughly 80 units over the whole horizon. Filling the energy ceiling of 60 units would eat almost all of it, leaving nothing to replace depreciation. The optimizer must find the right balance.

2. Source Code

The entire program is a single script. The simulator is fully vectorized, so it can evaluate hundreds of candidate investment schedules in one call. This makes the global search finish in seconds instead of minutes.

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

plt.style.use("dark_background")

# ------------------------------------------------------------------
# Model parameters
# ------------------------------------------------------------------
T = 20
P0 = 100.0
WORKER_FRAC = 0.6
NEED = 1.0
A_LAND = 9.5
BETA = 0.6
GAIN = 1.2
KAPPA = 0.5
ENERGY_SUPPLY = 30.0
ENERGY_PER_UNIT = 0.5
K_CAP = ENERGY_SUPPLY / ENERGY_PER_UNIT
K0 = 5.0
M0 = 400.0
M_RECYCLE = 4.0
M_PER_UNIT = 6.0
DEPREC = 0.06
FOOD_PER_UNIT = 2.0
R_MAX = 0.12
I_MAX = 12.0

HIST_KEYS = ("P", "K", "Keff", "M", "Y", "I", "F", "ypw")


# ------------------------------------------------------------------
# Vectorized simulator: each row of I_req is one investment schedule
# ------------------------------------------------------------------
def simulate(I_req):
    I_req = np.atleast_2d(np.asarray(I_req, dtype=float))
    n = I_req.shape[0]
    P = np.full(n, P0)
    K = np.full(n, K0)
    M = np.full(n, M0)
    total = np.zeros(n)
    hist = {key: np.zeros((n, T)) for key in HIST_KEYS}

    for t in range(T):
        L = WORKER_FRAC * P
        K_eff = np.minimum(K, K_CAP)
        Y = A_LAND * L ** BETA * (1.0 + GAIN * K_eff / (K_eff + KAPPA * L))

        I_cap = np.minimum(M / M_PER_UNIT, 0.5 * Y / FOOD_PER_UNIT)
        I = np.minimum(np.maximum(I_req[:, t], 0.0), I_cap)

        F = Y - FOOD_PER_UNIT * I
        q = F / (NEED * P)
        P_new = P * (1.0 + R_MAX * np.tanh(2.0 * (q - 1.0)))

        K = (1.0 - DEPREC) * K + I
        M = M - M_PER_UNIT * I + M_RECYCLE
        P = P_new
        total += P

        hist["P"][:, t] = P
        hist["K"][:, t] = K
        hist["Keff"][:, t] = K_eff
        hist["M"][:, t] = M
        hist["Y"][:, t] = Y
        hist["I"][:, t] = I
        hist["F"][:, t] = F
        hist["ypw"][:, t] = Y / L

    return total, hist


# ------------------------------------------------------------------
# Optimization: batch differential evolution + local polishing
# ------------------------------------------------------------------
def neg_total_batch(X):
    X = np.asarray(X, dtype=float)
    if X.ndim == 1:
        X = X.reshape(-1, 1)
    total, _ = simulate(X.T)
    return -total


def neg_total_scalar(x):
    total, _ = simulate(x)
    return -float(total[0])


bounds = [(0.0, I_MAX)] * T
best_x, best_val = None, np.inf

t_start = time.time()
for seed in (0, 1, 2, 3):
    res = differential_evolution(
        neg_total_batch,
        bounds,
        vectorized=True,
        updating="deferred",
        popsize=30,
        maxiter=800,
        tol=1e-10,
        mutation=(0.4, 1.0),
        recombination=0.85,
        seed=seed,
        polish=False,
    )
    pol = minimize(neg_total_scalar, res.x, method="L-BFGS-B", bounds=bounds)
    if pol.fun <= res.fun:
        cand_x, cand_val = pol.x, pol.fun
    else:
        cand_x, cand_val = res.x, res.fun
    print(f"seed {seed}: total population = {-cand_val:,.2f}")
    if cand_val < best_val:
        best_x, best_val = cand_x.copy(), cand_val
elapsed = time.time() - t_start
x_opt = np.clip(best_x, 0.0, I_MAX)

# ------------------------------------------------------------------
# Baseline strategies
# ------------------------------------------------------------------
const_grid = np.linspace(0.0, I_MAX, 241)
tot_const, _ = simulate(np.tile(const_grid[:, None], (1, T)))
c_best = float(const_grid[np.argmax(tot_const)])

strategies = {
    "Optimal": x_opt,
    "No automation": np.zeros(T),
    "Best constant": np.full(T, c_best),
    "Max effort": np.full(T, I_MAX),
}

results = {}
for name, x in strategies.items():
    total, h = simulate(x)
    results[name] = (float(total[0]), {k: v[0] for k, v in h.items()})

h_opt = results["Optimal"][1]
base_total = results["No automation"][0]

print()
print(f"Search finished in {elapsed:.1f} s")
print(f"Best constant investment per turn: {c_best:.2f} units")
print()
header = f"{'Strategy':<16}{'Total pop.':>14}{'Final pop.':>14}{'Final K':>10}{'Materials left':>16}{'vs. none':>11}"
print(header)
print("-" * len(header))
for name, (tot, h) in results.items():
    gain = 100.0 * (tot / base_total - 1.0)
    print(f"{name:<16}{tot:>14,.1f}{h['P'][-1]:>14,.1f}{h['K'][-1]:>10.1f}{h['M'][-1]:>16.1f}{gain:>10.1f}%")
print()
print("Optimal investment schedule (units built per turn):")
print(np.array2string(h_opt["I"], precision=2, floatmode="fixed", max_line_width=100))

# ------------------------------------------------------------------
# Two-parameter slice for the 3D surface:
# constant investment in early turns vs. constant investment in late turns
# ------------------------------------------------------------------
EARLY_END = 6
grid = np.linspace(0.0, I_MAX, 49)
GE, GL = np.meshgrid(grid, grid)
turn_idx = np.arange(T)
sched = np.where(turn_idx[None, None, :] < EARLY_END, GE[:, :, None], GL[:, :, None])
Z = simulate(sched.reshape(-1, T))[0].reshape(GE.shape)
best_idx = np.unravel_index(np.argmax(Z), Z.shape)

# ------------------------------------------------------------------
# Visualization: one figure, nine panels
# ------------------------------------------------------------------
COL = {
    "Optimal": "#00e5ff",
    "No automation": "#ff6b6b",
    "Best constant": "#ffd166",
    "Max effort": "#b388ff",
}
t_axis = np.arange(1, T + 1)
t_full = np.arange(0, T + 1)

fig = plt.figure(figsize=(21, 17))
gs = fig.add_gridspec(3, 3, hspace=0.38, wspace=0.27)
fig.subplots_adjust(left=0.05, right=0.97, top=0.92, bottom=0.05)
fig.suptitle("Automation vs. Scarcity: Maximizing Cumulative Population over 20 Turns",
             fontsize=21, fontweight="bold")

# (1) Population trajectories
ax1 = fig.add_subplot(gs[0, 0])
for name, (tot, h) in results.items():
    ax1.plot(t_full, np.r_[P0, h["P"]], color=COL[name],
             lw=3.4 if name == "Optimal" else 2.0,
             label=f"{name} (total {tot:,.0f})")
ax1.set_title("Population trajectory", fontsize=14)
ax1.set_xlabel("Turn")
ax1.set_ylabel("Population")
ax1.grid(alpha=0.25)
ax1.legend(fontsize=9, loc="upper left")

# (2) Optimal investment and materials
ax2 = fig.add_subplot(gs[0, 1])
bars = ax2.bar(t_axis, h_opt["I"], color=COL["Optimal"], alpha=0.85)
ax2.set_title("Optimal construction schedule and materials", fontsize=14)
ax2.set_xlabel("Turn")
ax2.set_ylabel("Equipment built per turn")
ax2.grid(alpha=0.25)
axm = ax2.twinx()
(line_m,) = axm.plot(t_axis, h_opt["M"], color="#ffa94d", marker="o", ms=4, lw=2.2)
axm.set_ylabel("Materials stock")
axm.grid(False)
ax2.legend([bars, line_m], ["Equipment built", "Materials stock"], loc="upper right", fontsize=9)

# (3) Equipment stock vs. energy ceiling
ax3 = fig.add_subplot(gs[0, 2])
for name, (tot, h) in results.items():
    ax3.plot(t_full, np.r_[K0, h["K"]], color=COL[name],
             lw=3.4 if name == "Optimal" else 1.8, label=name)
ax3.axhline(K_CAP, color="white", ls="--", lw=1.5, label="Energy-supported ceiling")
ax3.set_ylim(0, K_CAP * 1.12)
ax3.set_title("Equipment stock vs. energy ceiling", fontsize=14)
ax3.set_xlabel("Turn")
ax3.set_ylabel("Equipment units")
ax3.grid(alpha=0.25)
ax3.legend(fontsize=9, loc="upper left")

# (4) Food balance under the optimal plan
ax4 = fig.add_subplot(gs[1, 0])
P_start = np.r_[P0, h_opt["P"][:-1]]
demand = NEED * P_start
ax4.plot(t_axis, h_opt["Y"], color="#69f0ae", lw=2.6, label="Production")
ax4.plot(t_axis, h_opt["F"], color=COL["Optimal"], lw=2.6, label="Net food after construction")
ax4.plot(t_axis, demand, color="#ff8a80", lw=2.2, ls="--", label="Demand")
ax4.fill_between(t_axis, h_opt["F"], demand, where=h_opt["F"] >= demand,
                 interpolate=True, color="#69f0ae", alpha=0.18)
ax4.set_title("Food balance (optimal plan)", fontsize=14)
ax4.set_xlabel("Turn")
ax4.set_ylabel("Food units")
ax4.grid(alpha=0.25)
ax4.legend(fontsize=9, loc="upper left")

# (5) Output per worker
ax5 = fig.add_subplot(gs[1, 1])
for name, (tot, h) in results.items():
    ax5.plot(t_axis, h["ypw"], color=COL[name],
             lw=3.4 if name == "Optimal" else 1.8, label=name)
ax5.set_title("Food output per worker", fontsize=14)
ax5.set_xlabel("Turn")
ax5.set_ylabel("Output per worker")
ax5.grid(alpha=0.25)
ax5.legend(fontsize=9, loc="best")

# (6) Totals by strategy
ax6 = fig.add_subplot(gs[1, 2])
names = list(results.keys())
vals = [results[n][0] for n in names]
ax6.bar(range(len(names)), vals, color=[COL[n] for n in names], alpha=0.9)
for i, v in enumerate(vals):
    ax6.text(i, v, f"{v:,.0f}", ha="center", va="bottom", fontsize=11, fontweight="bold")
ax6.set_xticks(range(len(names)))
ax6.set_xticklabels(names, rotation=12)
ax6.set_ylim(0, max(vals) * 1.12)
ax6.set_title("Cumulative population by strategy", fontsize=14)
ax6.set_ylabel("Sum of population over 20 turns")
ax6.grid(alpha=0.25, axis="y")

# (7) 3D surface: early vs. late investment level
ax7 = fig.add_subplot(gs[2, 0], projection="3d")
ax7.plot_surface(GE, GL, Z, cmap="viridis", edgecolor="none", alpha=0.95, antialiased=True)
ax7.scatter([GE[best_idx]], [GL[best_idx]], [Z[best_idx]], color="white", s=70, depthshade=False)
ax7.set_title(f"Total population surface\n(early = turns 1-{EARLY_END}, late = turns {EARLY_END + 1}-{T})",
              fontsize=13)
ax7.set_xlabel("Early investment")
ax7.set_ylabel("Late investment")
ax7.set_zlabel("Total population")
ax7.view_init(elev=28, azim=-125)

# (8) 3D trajectories: time, equipment, population
ax8 = fig.add_subplot(gs[2, 1], projection="3d")
for name, (tot, h) in results.items():
    ax8.plot(t_full, np.r_[K0, h["K"]], np.r_[P0, h["P"]], color=COL[name],
             lw=3.4 if name == "Optimal" else 2.0, label=name)
ax8.set_title("State-space trajectories", fontsize=13)
ax8.set_xlabel("Turn")
ax8.set_ylabel("Equipment K")
ax8.set_zlabel("Population P")
ax8.view_init(elev=24, azim=-60)
ax8.legend(loc="upper left", fontsize=9)

# (9) Contour view of the same slice
ax9 = fig.add_subplot(gs[2, 2])
cf = ax9.contourf(GE, GL, Z, levels=24, cmap="viridis")
fig.colorbar(cf, ax=ax9, pad=0.02, label="Total population")
ax9.scatter([GE[best_idx]], [GL[best_idx]], marker="*", s=380, color="white",
            edgecolor="black", zorder=5, label="Best on slice")
ax9.plot([0, I_MAX], [0, I_MAX], ls="--", color="white", alpha=0.6, lw=1.2,
         label="Constant investment")
ax9.set_title("Slice contour: early vs. late investment", fontsize=14)
ax9.set_xlabel("Early investment (units per turn)")
ax9.set_ylabel("Late investment (units per turn)")
ax9.legend(fontsize=9, loc="upper right")

plt.show()

3. Console Output

seed 0: total population = 4,037.14
seed 1: total population = 4,037.14
seed 2: total population = 4,036.60
seed 3: total population = 4,037.14

Search finished in 32.7 s
Best constant investment per turn: 10.95 units

Strategy            Total pop.    Final pop.   Final K  Materials left   vs. none
---------------------------------------------------------------------------------
Optimal                4,037.1         269.0      31.5            16.0      50.1%
No automation          2,688.8         144.0       1.5           480.0       0.0%
Best constant          4,032.5         269.1      33.6             4.0      50.0%
Max effort             4,032.3         268.2      33.1             4.0      50.0%

Optimal investment schedule (units built per turn):
[1.20e+01 1.20e+01 1.19e+01 1.11e+01 9.79e+00 7.86e+00 5.13e+00 1.56e+00 6.67e-01 6.67e-01 6.67e-01
 6.67e-01 6.67e-01 6.67e-01 6.67e-01 6.67e-01 6.67e-01 1.24e-07 2.39e-08 5.05e-09]

4. Result Figure

5. Code Walkthrough

5.1 The vectorized simulator

simulate is the heart of the program. It accepts a matrix whose rows are candidate investment schedules and advances all of them through the 20 turns at the same time using NumPy array operations. Each turn follows the model exactly:

Compute workers $L_t=\eta P_t$ and the powered equipment $K_t^{\mathrm{eff}}=\min(K_t,K_{\max})$.
Compute production $Y_t$ with the saturating automation multiplier.
Convert the requested investment $u_t$ into the investment actually built. It is limited by the materials in stock and by a rule that construction can consume at most half of the food produced. Together these keep the model feasible for every schedule, so no candidate can produce negative stocks or a negative population.
Subtract the food cost of construction, compute the food ratio $q_t$, and update the population with the $\tanh$ growth law. The $\tanh$ keeps the growth rate within $\pm r$, so the population stays strictly positive.
Update the equipment and materials stocks and accumulate the objective.

Because every operation acts on whole arrays, evaluating 600 schedules costs only slightly more than evaluating one. That is the speed-up trick of this article.

5.2 Why a global search

The objective is non-convex and has kinks. The min operations switch regimes, for example when materials run out or when equipment exceeds the energy ceiling. Gradient-based methods started from a single point can easily get stuck. Differential evolution maintains a population of 600 schedules and recombines them, which explores the 20-dimensional space robustly.

Two details matter for speed and reliability:

vectorized=True with updating="deferred" makes SciPy pass the whole population to the objective at once. neg_total_batch receives an array of shape (parameters, solutions), transposes it, and calls the batch simulator.
polish=False disables SciPy’s internal polishing, because that step would call the batch objective with a single point. Instead we polish ourselves with L-BFGS-B using the scalar objective neg_total_scalar. The script runs four random seeds and keeps the best result, which guards against an unlucky run.

5.3 Baseline strategies

To judge the optimum, we compare it with three simple policies:

No automation. Never build anything. The initial equipment simply depreciates.
Best constant. Build the same amount every turn. The amount is found by scanning 241 values at once with the batch simulator.
Max effort. Build as much as possible (12 units) every turn until materials run dry.

5.4 The 3D slice

The full decision space has 20 dimensions, so we cannot plot it directly. Instead we study a two-dimensional slice: a constant investment level for turns 1 to 6 (the $x$-axis) and another constant level for turns 7 to 20 (the $y$-axis). One call to the batch simulator evaluates all 49 × 49 = 2,401 grid points, which become the surface in panel 7 and the contour map in panel 9.

6. Reading the Graphs

Panel 1 – Population trajectory. The optimal curve sits above all three baselines and, importantly, the legend shows the cumulative totals that serve as our objective. Even with no automation, the population grows at first because the land can feed more than 100 people. It then stalls as diminishing returns to labor set in. Automation lifts that ceiling, which is why the optimal curve keeps climbing after the other curves flatten.

Panel 2 – Construction schedule and materials. The bars show what was actually built, after the material and food caps. The orange line shows how the material stock is drawn down. Aggressive building early on converts the material stock into productive capacity while the payoff still has many turns to accumulate. As the stock approaches its limit, construction shrinks to roughly the level that recycling can afford.

Panel 3 – Equipment versus the energy ceiling. The dashed line is the maximum number of units the energy supply can power. Equipment above it would be wasted. The comparison shows the cost of the max-effort policy: it burns through materials quickly and then cannot replace depreciating machines, whereas the optimal plan paces its construction to what materials allow.

Panel 4 – Food balance. The gap between production and net food is the food diverted to construction. During heavy-building turns, net food falls and may drop below demand, which temporarily slows growth. This is the short-term price of investing, and the optimizer accepts it only when the later gains exceed the loss.

Panel 5 – Output per worker. This directly answers the opening question. Automation raises output per worker, but diminishing returns to labor pull it down as the population grows. The optimal policy keeps productivity per worker higher than the other strategies for most of the horizon.

Panel 6 – Totals. The bar chart summarizes the ranking, and the console table quantifies the gain over doing nothing.

Panel 7 – 3D surface. The height is the cumulative population as a function of early and late investment. The surface is not a simple ramp: more early investment helps only up to a point, because the food diverted to construction slows growth and the material stock is finite. The white dot marks the best point of this restricted slice.

Panel 8 – State-space trajectories. Each strategy traces a path through (turn, equipment, population). The optimal path rises steeply in both equipment and population, while the no-automation path stays low in equipment, and the max-effort path climbs in equipment but bends back as machines wear out without enough materials to replace them.

Panel 9 – Contour map. The same slice seen from above. The dashed diagonal marks constant-investment policies. If the best point lies off this diagonal, it shows that front-loading or back-loading beats spreading construction evenly.

7. Takeaways

Automation is an investment with a payback period. It reduces food available today in exchange for higher output later, so timing matters as much as quantity.
Materials act as a budget. With a fixed stock and slow recycling, every unit built early competes with maintenance later.
Energy sets a hard ceiling on how much equipment is useful, so the right question is not how much automation to build, but how much can be sustained.
A vectorized simulator combined with differential evolution and a local polish solves a non-smooth, 20-dimensional dynamic planning problem in seconds, and the same structure extends naturally to richer models with several resources, stochastic shocks, or multiple regions.

Optimizing Land Use for Population and Environmental Value Over 50 Years

November 5, 2026

Forest or Farmland?

Clearing forest for crops is tempting. Farmland feeds more people right away, and a growing population needs food. But forests also act as a sponge for rainfall: they recharge groundwater and keep rivers flowing. Remove too much of them and the water that irrigates those fields starts to disappear decades later.

In this article we turn that dilemma into a concrete optimization problem and solve it with Python. We search for the land-use policy that maximizes a combination of population and environmental resource value 50 years from now.

1. The Problem Setting

Consider a region with a total area of $L = 1000\ \mathrm{km^2}$. At the start:

Forest: $F_0 = 700\ \mathrm{km^2}$
Farmland: $A_0 = 300\ \mathrm{km^2}$
Population: $P_0 = 30$ thousand people
Water stock (aquifer and reservoirs): $W_0 = 800$ million $\mathrm{m^3}$, with capacity $W_{\max} = 1000$ million $\mathrm{m^3}$

Every year we may convert some forest into farmland. We restrict the decision to a simple two-parameter policy:

$u$: the fraction of the remaining forest converted each year
$T_c$: the number of years during which conversion continues, after which it stops

The conversion in year $t$ is

$$
c_t =
\begin{cases}
u,F_t & (t < T_c)\
0 & (t \ge T_c)
\end{cases}
$$

and the land-use dynamics are

$$
F_{t+1} = F_t - c_t, \qquad A_{t+1} = A_t + c_t .
$$

2. The Mathematical Model

Food production. Farmland yields $y_0$ thousand people’s worth of food per $\mathrm{km^2}$ at full water supply, scaled by a water factor:

$$
Y_t = y_0, A_t \sqrt{\frac{W_t}{W_{\max}}}
$$

The food ratio $r_t$ is production divided by population:

$$
r_t = \frac{Y_t}{P_t}
$$

Population dynamics. Population grows when $r_t > 1$ and shrinks when $r_t < 1$:

$$
P_{t+1} = P_t\left[1 + g_{\max}\tanh\bigl(2,(r_t - 1)\bigr)\right]
$$

Water dynamics. Forest increases recharge, while farmland irrigation and domestic use withdraw water. A fraction $\lambda$ is lost naturally each year:

$$
W_{t+1} = \mathrm{clip}\Bigl(W_t + a_0 + a_1 F_t - w_a A_t - w_p P_t - \lambda W_t,\ 0,\ W_{\max}\Bigr)
$$

This is the key mechanism. Converting forest cuts recharge $a_1 F_t$ and raises irrigation demand $w_a A_t$ at the same time, so the water stock declines slowly and the damage shows up late.

Objective function. We want to maximize the value at the horizon $T = 50$:

$$
J = w,\frac{P_T}{P_0} + (1-w)\left[\kappa,\frac{F_T}{F_0} + (1-\kappa),\frac{W_T}{W_{\max}}\right]
$$

The first term is the normalized population. The bracket is the normalized environmental value, a mix of forest area and water stock. We use $w = 0.5$ and $\kappa = 0.5$. The optimization problem is

$$
\max_{u,,T_c}\ J(u, T_c)
$$

Parameter values used in the code:

Symbol	Value	Meaning
$y_0$	0.11	People fed per $\mathrm{km^2}$ of farmland (thousand)
$a_0$	25	Base recharge (million $\mathrm{m^3}$/yr)
$a_1$	0.30	Extra recharge per $\mathrm{km^2}$ of forest
$\lambda$	0.04	Natural water loss rate
$w_a$	0.25	Irrigation withdrawal per $\mathrm{km^2}$ of farmland
$w_p$	0.60	Domestic withdrawal per thousand people
$g_{\max}$	0.03	Maximum population growth rate

3. Solution Strategy

The decision space has only two dimensions, so we can evaluate it exhaustively with a grid search:

$u$: 81 values from 0% to 10% per year
$T_c$: 51 values from 0 to 50 years

That gives $81 \times 51 = 4131$ policies. Running a 50-year simulation for each one in a Python loop would be slow. Instead, we write the simulation so that every state variable is a NumPy array of shape $(81, 51)$. A single 50-step loop then advances all 4131 policies at the same time.

We also run a sensitivity analysis on the population weight $w$. Because the final states of all policies are stored, we can re-evaluate $J$ for any $w$ without simulating again.

4. Source Code

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

plt.style.use("dark_background")

# ------------------------------------------------------------
# 1. Parameters
# ------------------------------------------------------------
L = 1000.0          # total land area [km^2]
F0 = 700.0          # initial forest [km^2]
A0 = L - F0         # initial farmland [km^2]
P0 = 30.0           # initial population [thousand people]
W0 = 800.0          # initial water stock [million m^3]
WMAX = 1000.0       # water storage capacity [million m^3]
T = 50              # planning horizon [years]

Y0 = 0.11           # people fed per km^2 of farmland at full water [thousand/km^2]
RECH0 = 25.0        # base recharge [million m^3/yr]
RECH1 = 0.30        # extra recharge per km^2 of forest [million m^3/yr/km^2]
LOSS = 0.04         # natural water loss rate [1/yr]
W_IRR = 0.25        # irrigation withdrawal per km^2 of farmland
W_DOM = 0.60        # domestic withdrawal per thousand people
GMAX = 0.03         # maximum population growth rate [1/yr]

W_POP = 0.5         # weight of population in the objective
KAPPA = 0.5         # weight of forest versus water inside environmental value

# ------------------------------------------------------------
# 2. Simulation (vectorized over arrays of policies)
# ------------------------------------------------------------
def simulate(u, tc, keep_path=False):
    u = np.asarray(u, dtype=float)
    tc = np.asarray(tc, dtype=float)
    u, tc = np.broadcast_arrays(u, tc)
    F = np.full(u.shape, F0)
    A = np.full(u.shape, A0)
    P = np.full(u.shape, P0)
    W = np.full(u.shape, W0)
    if keep_path:
        Fh, Ah, Ph, Wh = [F.copy()], [A.copy()], [P.copy()], [W.copy()]
    for t in range(T):
        conv = np.where(t < tc, u * F, 0.0)
        food = Y0 * A * np.sqrt(np.clip(W, 0.0, WMAX) / WMAX)
        ratio = food / np.maximum(P, 1e-9)
        P_next = P * (1.0 + GMAX * np.tanh(2.0 * (ratio - 1.0)))
        recharge = RECH0 + RECH1 * F
        withdraw = W_IRR * A + W_DOM * P
        W_next = np.clip(W + recharge - withdraw - LOSS * W, 0.0, WMAX)
        F = F - conv
        A = A + conv
        P = P_next
        W = W_next
        if keep_path:
            Fh.append(F.copy()); Ah.append(A.copy())
            Ph.append(P.copy()); Wh.append(W.copy())
    if keep_path:
        return np.array(Fh), np.array(Ah), np.array(Ph), np.array(Wh)
    return F, A, P, W

def objective(F, W, P, w_pop=W_POP):
    env = KAPPA * F / F0 + (1.0 - KAPPA) * W / WMAX
    return w_pop * P / P0 + (1.0 - w_pop) * env

# ------------------------------------------------------------
# 3. Grid search over (conversion rate u, conversion period tc)
# ------------------------------------------------------------
u_grid = np.linspace(0.0, 0.10, 81)
tc_grid = np.arange(0, T + 1, 1)
U, TC = np.meshgrid(u_grid, tc_grid, indexing="ij")
Ff, Af, Pf, Wf = simulate(U, TC)
J = objective(Ff, Wf, Pf)

idx = np.unravel_index(np.argmax(J), J.shape)
u_best, tc_best = u_grid[idx[0]], tc_grid[idx[1]]
J_best = J[idx]
print("=== Optimal land-use policy ===")
print(f"Annual conversion rate u*      : {u_best*100:.2f} % of remaining forest per year")
print(f"Conversion period tc*          : {tc_best} years")
print(f"Forest after 50 years          : {Ff[idx]:.1f} km^2 ({Ff[idx]/L*100:.1f} % of land)")
print(f"Farmland after 50 years        : {Af[idx]:.1f} km^2")
print(f"Population after 50 years      : {Pf[idx]:.1f} thousand")
print(f"Water stock after 50 years     : {Wf[idx]:.1f} million m^3")
print(f"Objective J*                   : {J_best:.4f}")

policies = {
    "No conversion": (0.0, 0),
    "Optimal": (u_best, tc_best),
    "Moderate (3%/yr, 20 yr)": (0.03, 20),
    "Aggressive (8%/yr, 50 yr)": (0.08, 50),
}
paths = {}
print("\n=== Policy comparison (year 50) ===")
print(f"{'Policy':<28}{'Forest':>9}{'Water':>9}{'Pop':>9}{'J':>9}")
for name, (uu, tt) in policies.items():
    Fp, Ap, Pp, Wp = simulate(uu, tt, keep_path=True)
    paths[name] = (Fp, Ap, Pp, Wp)
    Jp = objective(Fp[-1], Wp[-1], Pp[-1])
    print(f"{name:<28}{float(Fp[-1]):>9.1f}{float(Wp[-1]):>9.1f}{float(Pp[-1]):>9.1f}{float(Jp):>9.4f}")

# ------------------------------------------------------------
# 4. Sensitivity to the weight on population
# ------------------------------------------------------------
w_list = np.linspace(0.0, 1.0, 41)
opt_forest = []
opt_pop = []
for wv in w_list:
    Jw = objective(Ff, Wf, Pf, wv)
    k = np.unravel_index(np.argmax(Jw), Jw.shape)
    opt_forest.append(Ff[k] / L * 100.0)
    opt_pop.append(Pf[k])
opt_forest = np.array(opt_forest)
opt_pop = np.array(opt_pop)

# ------------------------------------------------------------
# 5. Visualization (single figure)
# ------------------------------------------------------------
years = np.arange(T + 1)
colors = {
    "No conversion": "#4fc3f7",
    "Optimal": "#ffd54f",
    "Moderate (3%/yr, 20 yr)": "#81c784",
    "Aggressive (8%/yr, 50 yr)": "#ef5350",
}

fig = plt.figure(figsize=(20, 15))
gs = GridSpec(3, 3, figure=fig, hspace=0.35, wspace=0.28)

# (a) 3D objective surface
ax1 = fig.add_subplot(gs[0, 0], projection="3d")
ax1.plot_surface(U * 100, TC, J, cmap="viridis", edgecolor="none", alpha=0.92)
ax1.scatter([u_best * 100], [tc_best], [J_best], color="#ff5252", s=90, depthshade=False)
ax1.set_xlabel("Conversion rate u [%/yr]")
ax1.set_ylabel("Conversion period [yr]")
ax1.set_zlabel("Objective J")
ax1.set_title("(a) Objective surface")
ax1.view_init(elev=28, azim=-125)

# (b) Contour map
ax2 = fig.add_subplot(gs[0, 1])
cs = ax2.contourf(U * 100, TC, J, levels=30, cmap="viridis")
fig.colorbar(cs, ax=ax2, label="Objective J")
ax2.plot(u_best * 100, tc_best, "o", color="#ff5252", markersize=10, label="Optimum")
ax2.set_xlabel("Conversion rate u [%/yr]")
ax2.set_ylabel("Conversion period [yr]")
ax2.set_title("(b) Contour map of J")
ax2.legend(loc="upper right")

# (c) 3D state-space trajectories (forest, water, population)
ax3 = fig.add_subplot(gs[0, 2], projection="3d")
for name, (Fp, Ap, Pp, Wp) in paths.items():
    ax3.plot(Fp, Wp, Pp, color=colors[name], linewidth=2.5, label=name)
    ax3.scatter([Fp[-1]], [Wp[-1]], [Pp[-1]], color=colors[name], s=50, depthshade=False)
ax3.set_xlabel("Forest [km$^2$]")
ax3.set_ylabel("Water [Mm$^3$]")
ax3.set_zlabel("Population [thousand]")
ax3.set_title("(c) State-space trajectories")
ax3.view_init(elev=24, azim=-60)
ax3.legend(loc="upper left", fontsize=8)

# (d) Population
ax4 = fig.add_subplot(gs[1, 0])
for name, (Fp, Ap, Pp, Wp) in paths.items():
    ax4.plot(years, Pp, color=colors[name], linewidth=2.5, label=name)
ax4.set_xlabel("Year")
ax4.set_ylabel("Population [thousand]")
ax4.set_title("(d) Population")
ax4.grid(alpha=0.25)
ax4.legend(fontsize=8)

# (e) Forest
ax5 = fig.add_subplot(gs[1, 1])
for name, (Fp, Ap, Pp, Wp) in paths.items():
    ax5.plot(years, Fp, color=colors[name], linewidth=2.5, label=name)
ax5.set_xlabel("Year")
ax5.set_ylabel("Forest area [km$^2$]")
ax5.set_title("(e) Forest area")
ax5.grid(alpha=0.25)

# (f) Water
ax6 = fig.add_subplot(gs[1, 2])
for name, (Fp, Ap, Pp, Wp) in paths.items():
    ax6.plot(years, Wp, color=colors[name], linewidth=2.5, label=name)
ax6.set_xlabel("Year")
ax6.set_ylabel("Water stock [Mm$^3$]")
ax6.set_title("(f) Water stock")
ax6.grid(alpha=0.25)

# (g) 3D bars: final population for coarse policy grid
ax7 = fig.add_subplot(gs[2, 0], projection="3d")
u_coarse = np.array([0.0, 0.02, 0.04, 0.06, 0.08, 0.10])
tc_coarse = np.array([0, 10, 20, 30, 40, 50])
Uc, Tcc = np.meshgrid(u_coarse, tc_coarse, indexing="ij")
_, _, Pc, _ = simulate(Uc, Tcc)
xpos = Uc.ravel() * 100
ypos = Tcc.ravel().astype(float)
zpos = np.zeros_like(xpos)
dz = Pc.ravel()
norm = (dz - dz.min()) / (dz.max() - dz.min() + 1e-12)
ax7.bar3d(xpos - 0.35, ypos - 1.5, zpos, 0.7, 3.0, dz, color=plt.cm.plasma(norm), shade=True)
ax7.set_xlabel("u [%/yr]")
ax7.set_ylabel("Period [yr]")
ax7.set_zlabel("Population [thousand]")
ax7.set_title("(g) Final population")
ax7.view_init(elev=30, azim=-130)

# (h) Trade-off between population and forest
ax8 = fig.add_subplot(gs[2, 1])
sc = ax8.scatter(opt_pop, opt_forest, c=w_list, cmap="plasma", s=60)
fig.colorbar(sc, ax=ax8, label="Population weight w")
ax8.set_xlabel("Optimal final population [thousand]")
ax8.set_ylabel("Remaining forest [% of land]")
ax8.set_title("(h) Trade-off frontier")
ax8.grid(alpha=0.25)

# (i) Optimal forest share against weight
ax9 = fig.add_subplot(gs[2, 2])
ax9.plot(w_list, opt_forest, color="#ffd54f", linewidth=2.5)
ax9.fill_between(w_list, opt_forest, color="#ffd54f", alpha=0.2)
ax9.set_xlabel("Population weight w")
ax9.set_ylabel("Forest kept at optimum [% of land]")
ax9.set_title("(i) Optimal forest share vs. weight")
ax9.grid(alpha=0.25)

fig.suptitle("Forest Conversion Policy: Population versus Environmental Value over 50 Years",
             fontsize=18, y=0.95)
plt.show()

5. Code Walkthrough

Section 1: Parameters. All physical constants and weights are defined at the top, so every experiment in this article can be reproduced by changing a single line. The area is in $\mathrm{km^2}$, population in thousands, and water in millions of $\mathrm{m^3}$.

Section 2: The simulate function. This is the core of the program.

np.broadcast_arrays(u, tc) makes the policy parameters share one shape. We can pass either scalars (a single policy) or a 2D grid (thousands of policies) and use the same code.
Each state variable (F, A, P, W) is an array of that shape, so one for t in range(T) loop simulates every policy at once. This is the speed-up: it replaces roughly four thousand 50-step Python loops with one.
np.where(t < tc, u * F, 0.0) implements the rule “convert a fraction $u$ of the forest while $t < T_c$, then stop”, elementwise for every policy.
The food ratio uses np.sqrt(np.clip(W, 0.0, WMAX) / WMAX). The clip keeps the square root safe, and np.maximum(P, 1e-9) avoids division by zero.
The water update is the equation from Section 2. Forest adds recharge, farmland and people subtract withdrawals, and the whole stock is clipped to $[0, W_{\max}]$.
Recharge and withdrawal are computed from the land state before this year’s conversion is applied, so the effect of deforestation reaches the water stock with a delay.
With keep_path=True, the function also returns the full time series. We use this later to draw trajectories.

Section 3: Grid search. np.meshgrid(..., indexing="ij") builds the policy grid, and a single call to simulate(U, TC) returns the final state of every policy. np.argmax followed by np.unravel_index finds the best cell. The program then prints the optimal policy and compares it against three reference policies: no conversion, a moderate one (3% per year for 20 years), and an aggressive one (8% per year for 50 years).

Section 4: Sensitivity analysis. The objective is cheap to re-evaluate because the final states Ff, Wf, Pf are already stored. For 41 values of the population weight $w$ from 0 to 1, we find the best policy and record the remaining forest share and the population it achieves.

Section 5: Visualization. Nine panels are combined into one figure with GridSpec:

(a) A 3D surface of the objective $J$ over the policy space, with the optimum marked.
(b) The same function as a contour map.
(c) A 3D trajectory of each policy in (forest, water, population) space.
(d)–(f) Time series of population, forest area, and water stock.
(g) A 3D bar chart of final population over a coarse policy grid.
(h)–(i) The trade-off between population and forest as the weight $w$ changes.

6. Execution Results

=== Optimal land-use policy ===
Annual conversion rate u*      : 6.50 % of remaining forest per year
Conversion period tc*          : 4 years
Forest after 50 years          : 535.0 km^2 (53.5 % of land)
Farmland after 50 years        : 465.0 km^2
Population after 50 years      : 49.8 thousand
Water stock after 50 years     : 997.5 million m^3
Objective J*                   : 1.2700

=== Policy comparison (year 50) ===
Policy                         Forest    Water      Pop        J
No conversion                   700.0   1000.0     32.8   1.0475
Optimal                         535.0    997.5     49.8   1.2700
Moderate (3%/yr, 20 yr)         380.7      0.0     23.3   0.5238
Aggressive (8%/yr, 50 yr)        10.8      0.0     14.5   0.2460

7. Interpreting the Results

The optimal policy is an early, short burst of conversion. The grid search selects a conversion rate of about 6.5% per year for only 4 years. That converts roughly 165 $\mathrm{km^2}$ of forest, so about 53% of the land stays forested at year 50. Farmland grows from 300 to about 465 $\mathrm{km^2}$, and the population reaches about 49.8 thousand, against about 32.8 thousand if nothing is converted. Meanwhile the water stock stays essentially full, at about 997 million $\mathrm{m^3}$.

Panel (a) and (b): a narrow ridge, not a single peak. The objective surface shows a bright ridge running along the curve where the total amount of forest removed, roughly $u \times T_c$, is moderate. Many combinations of $u$ and $T_c$ that clear about the same area score about the same. The decision that matters is how much forest to remove in total, and the exact schedule matters much less. The surface then falls steeply toward the far corner, where large amounts of forest are cleared.

Panel (d)–(f): the delayed collapse. The aggressive and moderate policies show the trap clearly. Their populations climb faster than the optimal one for the first 15 to 25 years, because more farmland means more food. Their forest keeps shrinking, so recharge falls and irrigation demand keeps rising. Eventually the water stock crashes to zero, in panel (f) after roughly 15 years for the aggressive policy and 30 years for the moderate one. The water factor in food production then collapses, and the population falls, ending below even the no-conversion case. In the end, the short-term gain is far outweighed by the long-term loss.

The comparison at year 50 makes this concrete:

Policy	Forest ($\mathrm{km^2}$)	Water (Mm$^3$)	Population (thousand)	$J$
No conversion	700.0	1000.0	32.8	1.0475
Optimal	535.0	997.5	49.8	1.2700
Moderate (3%/yr, 20 yr)	380.7	0.0	23.3	0.5238
Aggressive (8%/yr, 50 yr)	10.8	0.0	14.5	0.2460

Panel (c): state-space trajectories. In the 3D view, the optimal and no-conversion paths stay near the high-water plane. The moderate and aggressive paths first move toward lower forest while population rises, then turn and plunge to zero water. The picture shows that the destructive policies pass through a high-population state that cannot be sustained.

Panel (g): a cliff, not a slope. The final-population bars show that the population is not a smooth function of the policy. Cautious policies give moderate populations, a narrow band of mild conversion gives the tallest bars, and heavy conversion causes a sudden drop. This is typical of systems with a tipping point, where the water stock hits its lower bound.

Panel (h) and (i): how values change the answer. As the population weight $w$ approaches 0, the optimum keeps about 70% of the land as forest, which is simply the do-nothing policy. At around $w \approx 0.15$ to $0.2$ the optimum switches to a limited conversion that leaves a little over half the land forested, and it stays there for all larger $w$. Even a society that values only population should not clear everything, because a policy that destroys the water supply also destroys the population. This is why the environmental and population objectives are less opposed than they first appear, at least up to the tipping point.

8. Summary

We modeled forest-to-farmland conversion with land, water, food, and population dynamics, and searched for the policy that maximizes $J$ over 50 years.
Vectorizing the simulation lets a single 50-step loop evaluate 4131 policies at once.
The best policy is a brief, early conversion of about 165 $\mathrm{km^2}$ of forest, which raises the final population by about 50% while keeping the water stock intact.
Policies that clear more forest look better for 15 to 25 years and then collapse as the water supply fails.
The timing of damage matters as much as its size. Because deforestation affects water with a delay, the true cost of a decision only becomes visible long after it is made.

The same framework extends naturally to richer models, for example adding soil erosion, a climate-driven rainfall trend, or a continuous-time control law solved with dynamic programming.

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.

Planning Long-Term Civilization Survival

1. The Example Problem

2. Mathematical Model

3. Solution Strategy

4. Source Code

5. Code Walkthrough

6. Execution Results

Console Output

Graph Output

7. Interpreting the Results

8. Conclusion

Automation vs. Scarcity

1. The Problem

Model

Parameters

2. Source Code

3. Console Output

4. Result Figure

5. Code Walkthrough

5.1 The vectorized simulator

5.2 Why a global search

5.3 Baseline strategies

5.4 The 3D slice

6. Reading the Graphs

7. Takeaways

Forest or Farmland?

1. The Problem Setting

2. The Mathematical Model

3. Solution Strategy

4. Source Code

5. Code Walkthrough

6. Execution Results

7. Interpreting the Results

8. Summary

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

Parameter decoding (`decode`)