Optimal Immigration and Birth Support Over 30 Turns

October 28, 2026

Steering Population Growth:

Introduction

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

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

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

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

Problem Setting

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

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

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

Production and food balance

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

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

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

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

Budget constraint

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

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

Demographics

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

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

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

Objective

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

Parameters

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

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

Solution Strategy

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

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

Full Source Code

import math
import time

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

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


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


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

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

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

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

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

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

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


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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Detailed Code Walkthrough

Parameters and helper

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

The simulator

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

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

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

Landscape analysis

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

Full optimization

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

Local polish

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

Reporting and plotting

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

Execution Results

Result Image

Console Output

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

Reading the Graphs

Panel 1: 3D landscape of constant policies

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

Panel 2: Population trajectories

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

Panel 3: Optimized controls

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

Panel 4: Food stock and satisfaction

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

Panel 5: 3D state-space trajectories

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

Panel 6: Objective value by policy

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

Discussion

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

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

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

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

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

Conclusion

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

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

October 27, 2026

Investing in Health vs. Machines

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

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

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

1. The Model

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

Production

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

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

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

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

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

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

Mortality

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

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

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

Population dynamics

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

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

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

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

The optimization problem

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

Example setup

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

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

2. The Code

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

plt.tight_layout()
plt.show()

3. Code Walkthrough

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

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

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

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

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

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

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

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

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

4. Results

Console output

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

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

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

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

Graph output

5. Interpreting the Results

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

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

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

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

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

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

6. Takeaways

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

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

Optimal Food Storage Under Spoilage and Storage Costs

October 26, 2026

Eat Now or Store for Later?

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

1. The Problem

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

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

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

The tension is easy to state:

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

2. Mathematical Model

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

Food balance.

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

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

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

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

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

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

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

Objective.

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

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

3. Concrete Example

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

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

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

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

4. Solution Method and Speed-Up

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

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

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

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

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

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

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

5. Full Source Code

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

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


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


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


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


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


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

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


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


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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6. Code Walkthrough

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

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

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

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

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

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

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

7. Results

Console Output

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

Result Image

8. How to Read the Results

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

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

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

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

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

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

9. Takeaways

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

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

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

Optimizing Construction Timing over 20 Turns

October 25, 2026

When Should You Build the Farm?

In strategy games and in real capacity planning, the same dilemma keeps showing up. A production facility pays off only after a delay, and the resources spent building it are taken from the very thing you want to grow. A farm that starts producing food three turns after construction begins is a clean example. Build too early and the construction cost starves the young population. Build too late and the population outgrows its food supply.

In this article we turn that dilemma into a concrete optimization problem, solve it exactly in Python, and visualize why the optimum looks the way it does.

1. Problem Setting

We plan over $T = 20$ turns. A village starts with $P_0 = 50$ people and an initial food stock $S_0 = 80$. At most $K = 5$ farms can be built, and each farm has three properties.

It costs $C = 90$ food units in total, paid in equal installments over the three turns of construction.
It starts producing $F = 16$ food per turn exactly $3$ turns after construction begins.
Construction may start at any turn $b_k \in {1, \dots, 20}$, or never (we encode “never” as $b_k = 21$).

A build plan is therefore the vector $\mathbf{b} = (b_1, \dots, b_K)$.

State transition

The number of farms producing at turn $t$ and the number under construction are

$$
A_t(\mathbf{b}) = \sum_{k=1}^{K} \mathbb{1}\left[b_k + 3 \le t\right], \qquad
B_t(\mathbf{b}) = \sum_{k=1}^{K} \mathbb{1}\left[b_k \le t < b_k + 3\right].
$$

Food available in turn $t$, after paying construction installments, is

$$
V_t = S_{t-1} + \left(M + F,A_t\right) - \frac{C}{3},B_t ,
$$

where $M = 40$ is the output of the existing fields. The demand is $D_t = c,P_{t-1}$ with $c = 0.8$ food per person. The satisfaction ratio is $\rho_t = \max\left(V_t / D_t,\ 0\right)$.

Population grows when food is plentiful and shrinks under shortage:

$$
P_t = P_{t-1}\left[1 + r,\mathrm{clip}!\left(\frac{\rho_t - 1}{0.3},,0,,1\right)\left(1 - \frac{P_{t-1}}{N_{\max}}\right) - \delta,\mathrm{clip}\left(1-\rho_t,,0,,1\right)\right]
$$

with $r = 0.14$, $\delta = 0.5$ and $N_{\max} = 600$. Leftover food is stored with $10%$ spoilage:

$$
S_t = 0.9,\max\left(V_t - D_t,\ 0\right).
$$

Objective

$$
\max_{\mathbf{b}\in{1,\dots,21}^K} ; P_{20}(\mathbf{b})
$$

Because farms are interchangeable, only the multiset of build turns matters. The number of distinct plans is

$$
\binom{21 + 5 - 1}{5} = 53{,}130 ,
$$

which is small enough to enumerate exhaustively and obtain a provably optimal plan.

2. Source Code

The simulator is vectorized over all plans at once. Instead of calling a Python loop 53,130 times, it runs a single 20-step loop on arrays of length 53,130. This is what makes the full search finish in a fraction of a second.

import itertools
import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

plt.style.use("dark_background")

# ------------------------------------------------------------
# Parameters
# ------------------------------------------------------------
T = 20            # planning horizon (turns)
DELAY = 3         # a farm starts producing 3 turns after construction begins
KMAX = 5          # maximum number of farms
NEVER = T + 1     # code meaning "this farm is never built"

P0 = 50.0         # initial population
S0 = 80.0         # initial food stock
R = 0.14          # maximum population growth rate per turn
DEATH = 0.5       # starvation death rate coefficient
BASE = 40.0       # food produced per turn by the existing fields
F = 16.0          # food produced per turn by one farm
CONS = 0.8        # food consumed per person per turn
COST = 90.0       # total construction cost of one farm (paid over DELAY turns)
SPOIL = 0.9       # fraction of leftover food that survives to the next turn
CAP = 600.0       # carrying capacity of the village

# ------------------------------------------------------------
# Vectorized simulator: evaluates many build plans at once
# ------------------------------------------------------------
def simulate(plans, record=False):
    plans = np.atleast_2d(plans)
    n = plans.shape[0]
    P = np.full(n, P0)
    S = np.full(n, S0)
    if record:
        hist = {k: [] for k in ("P", "S", "prod", "need", "active")}
        hist["P"].append(P.copy())
        hist["S"].append(S.copy())
    for t in range(1, T + 1):
        active = ((plans + DELAY) <= t).sum(axis=1)
        building = ((plans <= t) & (t < plans + DELAY)).sum(axis=1)
        prod = BASE + F * active
        avail = S + prod - (COST / DELAY) * building
        need = CONS * P
        ratio = np.clip(avail / need, 0.0, None)
        growth = R * np.clip((ratio - 1.0) / 0.3, 0.0, 1.0) * (1.0 - P / CAP)
        deaths = DEATH * np.clip(1.0 - ratio, 0.0, 1.0)
        P = P * (1.0 + growth - deaths)
        S = np.maximum(avail - need, 0.0) * SPOIL
        if record:
            hist["P"].append(P.copy())
            hist["S"].append(S.copy())
            hist["prod"].append(prod.copy())
            hist["need"].append(need.copy())
            hist["active"].append(active.copy())
    if record:
        return P, {k: np.array(v) for k, v in hist.items()}
    return P

# ------------------------------------------------------------
# Exhaustive search over all build plans
# ------------------------------------------------------------
start = time.time()
plans = np.array(
    list(itertools.combinations_with_replacement(range(1, NEVER + 1), KMAX)),
    dtype=int,
)
final_pop = simulate(plans)
elapsed = time.time() - start

n_farms = (plans != NEVER).sum(axis=1)
top = final_pop.max()
candidates = np.where(final_pop >= top - 1e-6)[0]
best_idx = int(candidates[np.argmin(n_farms[candidates])])
best_plan = plans[best_idx]
best_value = final_pop[best_idx]

print(f"Number of candidate plans : {len(plans):,}")
print(f"Search time               : {elapsed:.3f} sec")
print(f"Optimal build turns       : {[int(x) for x in best_plan if x != NEVER]}")
print(f"Final population (turn {T}): {best_value:.1f}")

# Best plan for each number of farms
best_by_k = []
for k in range(KMAX + 1):
    mask = (plans == NEVER).sum(axis=1) == (KMAX - k)
    idx = np.where(mask)[0][np.argmax(final_pop[mask])]
    best_by_k.append((k, plans[idx], final_pop[idx]))
    built = [int(x) for x in plans[idx] if x != NEVER]
    print(f"Best with {k} farm(s): turns {built}  ->  population {final_pop[idx]:.1f}")

# Reference strategies
ref_plans = {
    "Optimal": best_plan,
    "All at turn 1": np.array([1] * KMAX),
    "Evenly spaced (1,4,7,10,13)": np.array([1, 4, 7, 10, 13]),
    "All late (12-16)": np.array([12, 13, 14, 15, 16]),
    "No farms": np.array([NEVER] * KMAX),
}
ref_results = {}
for name, plan in ref_plans.items():
    P_end, h = simulate(plan, record=True)
    ref_results[name] = (plan, float(P_end[0]), h)
    print(f"{name:30s}: final population {P_end[0]:7.1f}")

# ------------------------------------------------------------
# Surface: first two farms (b1, b2) with the other three at optimum
# ------------------------------------------------------------
fixed = [int(x) for x in best_plan[2:]]
grid = np.arange(1, T + 1)
B1, B2 = np.meshgrid(grid, grid)
surf_plans = np.stack(
    [B1.ravel(), B2.ravel()] + [np.full(B1.size, v) for v in fixed], axis=1
)
Z = simulate(surf_plans).reshape(B1.shape)

# ------------------------------------------------------------
# Visualization (single figure)
# ------------------------------------------------------------
fig = plt.figure(figsize=(20, 12))
fig.suptitle("Optimal Farm Construction Timing over 20 Turns", fontsize=20, y=0.98)

colors = ["#00e5ff", "#ff5252", "#ffd740", "#b388ff", "#9e9e9e"]

# (1) 3D surface
ax1 = fig.add_subplot(2, 3, 1, projection="3d")
surf = ax1.plot_surface(B1, B2, Z, cmap=cm.viridis, edgecolor="none", alpha=0.95)
ax1.scatter(
    [best_plan[0]], [best_plan[1]], [best_value],
    color="red", s=80, depthshade=False, label="Optimum",
)
ax1.set_xlabel("Build turn of farm 1")
ax1.set_ylabel("Build turn of farm 2")
ax1.set_zlabel("Final population")
ax1.set_title("Final population vs. first two build turns")
ax1.view_init(elev=28, azim=-60)
ax1.legend(loc="upper left")
fig.colorbar(surf, ax=ax1, shrink=0.55, pad=0.1)

# (2) 3D trajectories
ax2 = fig.add_subplot(2, 3, 2, projection="3d")
turns = np.arange(0, T + 1)
for (name, (plan, val, h)), c in zip(ref_results.items(), colors):
    ax2.plot(turns, h["P"][:, 0], h["S"][:, 0], color=c, lw=2.2, label=name)
ax2.set_xlabel("Turn")
ax2.set_ylabel("Population")
ax2.set_zlabel("Food stock")
ax2.set_title("3D trajectory: turn, population, food stock")
ax2.view_init(elev=24, azim=-55)
ax2.legend(fontsize=8, loc="upper left")

# (3) Population trajectories
ax3 = fig.add_subplot(2, 3, 3)
for (name, (plan, val, h)), c in zip(ref_results.items(), colors):
    ax3.plot(turns, h["P"][:, 0], color=c, lw=2.4, label=f"{name} ({val:.0f})")
ax3.set_xlabel("Turn")
ax3.set_ylabel("Population")
ax3.set_title("Population over time")
ax3.grid(alpha=0.25)
ax3.legend(fontsize=8)

# (4) Best value by number of farms
ax4 = fig.add_subplot(2, 3, 4)
ks = [k for k, _, _ in best_by_k]
vals = [v for _, _, v in best_by_k]
bars = ax4.bar(ks, vals, color=cm.plasma(np.linspace(0.2, 0.9, len(ks))))
for b, v in zip(bars, vals):
    ax4.text(b.get_x() + b.get_width() / 2, v + 3, f"{v:.0f}", ha="center")
ax4.set_xlabel("Number of farms built")
ax4.set_ylabel("Best final population")
ax4.set_title("Marginal value of each additional farm")
ax4.grid(alpha=0.25, axis="y")

# (5) Food balance of optimal plan
ax5 = fig.add_subplot(2, 3, 5)
h = ref_results["Optimal"][2]
tt = np.arange(1, T + 1)
ax5.plot(tt, h["prod"][:, 0], color="#76ff03", lw=2.4, label="Food production")
ax5.plot(tt, h["need"][:, 0], color="#ff5252", lw=2.4, label="Food demand")
ax5.fill_between(
    tt, h["prod"][:, 0], h["need"][:, 0],
    where=h["prod"][:, 0] >= h["need"][:, 0], color="#76ff03", alpha=0.2,
)
ax5.fill_between(
    tt, h["prod"][:, 0], h["need"][:, 0],
    where=h["prod"][:, 0] < h["need"][:, 0], color="#ff5252", alpha=0.2,
)
for b in best_plan:
    if b != NEVER:
        ax5.axvline(b, color="#ffd740", ls=":", alpha=0.7)
        if b + DELAY <= T:
            ax5.axvline(b + DELAY, color="#00e5ff", ls="--", alpha=0.5)
ax5.set_xlabel("Turn  (yellow: construction start, cyan: production start)")
ax5.set_ylabel("Food per turn")
ax5.set_title("Food balance under the optimal plan")
ax5.grid(alpha=0.25)
ax5.legend()

# (6) Distribution of all plans
ax6 = fig.add_subplot(2, 3, 6)
ax6.hist(final_pop, bins=60, color="#40c4ff", alpha=0.85)
ax6.axvline(best_value, color="red", lw=2.5, label=f"Optimum = {best_value:.0f}")
ax6.set_xlabel("Final population")
ax6.set_ylabel("Number of plans")
ax6.set_title("Outcomes of all 53,130 build plans")
ax6.grid(alpha=0.25)
ax6.legend()

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

3. Code Walkthrough

Parameters. All model constants live at the top, so you can change the cost, the delay or the farm output and immediately see how the optimal schedule shifts. NEVER = T + 1 is a sentinel meaning “this farm is never built”. With it, plans with fewer than five farms are represented in the same fixed-length array.

The simulator. simulate accepts an array of plans of shape $(n, K)$ and advances all of them through the 20 turns simultaneously.

active counts farms whose production has started, which is the condition $b_k + 3 \le t$.
building counts farms in their three construction turns, which is the condition $b_k \le t < b_k + 3$. Each pays one third of the cost per turn.
avail is the stock plus production minus the construction installment. This is where the cost of building early enters the model.
ratio is the satisfaction ratio $\rho_t$. Growth switches on smoothly once food exceeds $100%$ of demand and saturates at $130%$. Deaths grow linearly as the ratio falls below $100%$.
S carries leftover food forward with $10%$ spoilage. Spoilage matters because it stops the model from solving the problem by simply hoarding a huge stockpile early.

Exhaustive search. itertools.combinations_with_replacement generates every multiset of five build turns drawn from ${1,\dots,21}$, giving exactly $53{,}130$ rows. A single call to simulate evaluates all of them, so the optimum is guaranteed, not heuristic.

Tie-breaking. If several plans reach the same final population, the code picks the one with the fewest farms. A farm that adds nothing but still costs resources should not appear in the answer.

Best plan per farm count. The loop over k finds the best schedule for exactly $k$ farms. This isolates the marginal value of each additional farm.

Reference strategies. Four intuitive strategies, namely building everything immediately, building at an even pace, building everything late and building nothing, are simulated with record=True to keep their full histories for plotting.

Surface data. The last three build turns are frozen at their optimal values, while the first two sweep a $20 \times 20$ grid. This produces a smooth landscape in which the effect of timing is easy to see.

4. Visualization of the Results

How to read the six panels

Panel 1, the 3D surface. The two horizontal axes are the build turns of the first and second farms, and the height is the final population. The surface rises steeply when both farms are built early and then falls off as either one is pushed late. The red dot marks the optimum. There is a pronounced ridge, not a flat plateau, which shows that timing matters a great deal. Because farms are interchangeable, the surface is also symmetric about the diagonal.

Panel 2, the 3D trajectories. Each curve traces turn, population and food stock together. The optimal plan keeps the stock low in the early turns, because every unit of food is being invested in construction, and then lets the stock climb as the farms come online. The all-late strategy collapses toward zero population, and the no-farm strategy stays pinned at the starting point.

Panel 3, population over time. This is the most direct comparison. The optimal plan overtakes the evenly spaced plan and ends clearly above it. The strategy that builds everything at turn 1 is the surprise. Paying for five farms at once leaves nothing to feed the population, and it ends almost where it started.

Panel 4, marginal value of each farm. The final population rises sharply with each of the first few farms and then flattens. The fourth farm still helps, while a fifth adds nothing, because by then the village already grows at its maximum rate and extra food is wasted to spoilage. This is why the optimal plan contains four farms and not five.

Panel 5, food balance. Production is a staircase that rises each time a farm completes, and demand is a smooth exponential curve. Yellow dotted lines mark construction starts and cyan dashed lines mark production starts. Each new farm is scheduled to arrive just before demand would overtake supply, so shortages are avoided without paying for capacity long before it is needed.

Panel 6, distribution of all plans. Most of the $53{,}130$ plans cluster at poor outcomes, and the optimum sits at the far right tail. A randomly chosen schedule is very unlikely to be good, which shows why a systematic search is worth the effort.

5. Console Output

Number of candidate plans : 53,130
Search time               : 0.177 sec
Optimal build turns       : [1, 3, 6, 9]
Final population (turn 20): 179.3
Best with 0 farm(s): turns []  ->  population 50.0
Best with 1 farm(s): turns [5]  ->  population 114.2
Best with 2 farm(s): turns [1, 2]  ->  population 149.8
Best with 3 farm(s): turns [1, 3, 6]  ->  population 168.6
Best with 4 farm(s): turns [1, 3, 6, 9]  ->  population 179.3
Best with 5 farm(s): turns [1, 3, 6, 9, 11]  ->  population 179.3
Optimal                       : final population   179.3
All at turn 1                 : final population    53.9
Evenly spaced (1,4,7,10,13)   : final population   167.1
All late (12-16)              : final population     3.3
No farms                      : final population    50.0

6. Discussion

The exhaustive search shows that the best schedule is a staggered one, with farms started at turns $1, 3, 6$ and $9$. Three lessons come out of it.

Build early, but not all at once. The first farm goes up at turn 1 because its payback is long and the population is still small enough to absorb the cost. Dumping every farm into the first turn is catastrophic, since the combined cost starves the population before any farm produces.
Spacing follows the delay. The gaps between builds (2, 3 and 3 turns) are close to the 3-turn production delay. Each farm’s construction cost is paid while the previous farm is still not producing, and the schedule is tuned so that the food dip never becomes deadly.
Stop at the right number. The marginal value of capacity drops to zero once growth is limited by the maximum growth rate and not by food. Building a fifth farm only burns resources.

The same method scales to richer settings. You can raise KMAX, add several facility types, or make the cost depend on the current population. When the plan space becomes too large to enumerate, dynamic programming or a genetic algorithm can replace the brute-force search, but the vectorized simulator can be reused unchanged as the evaluation function.

Optimizing Mining and Reinvestment to Maximize the Final Population

October 24, 2026

How Much Ore Should You Dig?

Ore deposits are finite. Dig too slowly and the colony never grows. Dig too fast and the mine runs dry while people still depend on it. Selling ore buys food and machines. Food keeps people alive and growing. Machines raise the extraction capacity. This post turns that trade-off into a concrete optimization problem and solves it in Python.

1. The Problem

A colony lives on a single ore deposit. Each year it decides:

how many tons to extract ($E$),
how much of the revenue to reinvest in machines (the rest buys food).

The goal is to maximize the population at the end of year $T$ without wasting the deposit or running it dry too early.

2. Mathematical Model

Let $t = 0, 1, \dots, T-1$ be the year, $R_t$ the remaining ore, $K_t$ the number of machines and $P_t$ the population.

Extraction is limited by the target, the remaining ore and the machine capacity:

$$
e_t = \min\left(E,; R_t,; c,K_t\right)
$$

Revenue and its split between reinvestment and food:

$$
\text{rev}_t = p,e_t, \qquad I_t = s_t,\text{rev}_t, \qquad F_t = (1-s_t),\text{rev}_t
$$

The reinvestment share is $s_t = s$ during the first $\tau$ years and $s_t = 0$ afterwards (all revenue goes to food).

State transitions:

$$
R_{t+1} = R_t - e_t, \qquad K_{t+1} = (1-\delta)K_t + \frac{I_t}{m}
$$

Food sufficiency and population growth:

$$
z_t = \frac{F_t}{n,P_t}, \qquad P_{t+1} = P_t\left(1 + 2g_{\max},\operatorname{clip}(z_t - 1,,-1,,0.5)\right)
$$

With $z_t = 1$ the population is stable, $z_t = 0$ means starvation ($-10%$ per year) and $z_t \ge 1.5$ gives the maximum growth $+g_{\max}$.

Optimization problem:

$$
\max_{E,,s,,\tau}; P_T \quad \text{s.t.}\quad \sum_{t} e_t \le R_0,;; 0 \le s \le 1,;; 0 \le \tau \le T
$$

The deposit constraint $\sum_t e_t \le R_0$ is enforced automatically by $e_t \le R_t$.

Parameters used in this example:

Symbol	Meaning	Value
$T$	Horizon	30 years
$R_0$	Initial ore	1000 t
$P_0$	Initial population	500
$K_0$	Initial machines	10
$p$	Ore price	40
$c$	Capacity per machine	2 t/year
$m$	Machine price	400
$\delta$	Depreciation	5%
$n$	Food need per person	1
$g_{\max}$	Max growth rate	5%

3. Solution Strategy

The decision space $(E, s, \tau)$ has only three dimensions, so the problem is solved in two stages:

Exhaustive search over a grid of $71 \times 51 \times 31 = 112{,}251$ scenarios. All scenarios are simulated at the same time with NumPy broadcasting, so the year loop runs only 30 times instead of $30 \times 112{,}251$ times.
Continuous refinement of the best grid point with Nelder–Mead.

The optimum is then compared with three baseline policies:

(a) never reinvest,
(b) always reinvest 80% of revenue,
(c) over-mine with $E = 80$ t/year.

4. Source Code

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

# ------------------------------------------------------------
# 1. Model parameters
# ------------------------------------------------------------
T = 30                    # planning horizon (years)
R0, P0, K0 = 1000.0, 500.0, 10.0   # initial ore [t], population, machines
PRICE = 40.0              # ore price per ton
CAP = 2.0                 # tons per machine per year
MACH = 400.0              # price of one machine
DEP = 0.05                # machine depreciation rate
NEED = 1.0                # food units required per person per year
GMAX = 0.05               # max population growth rate per year


def simulate(E, s, tau, record=False):
    """Vectorised simulation. E, s, tau may be scalars or broadcastable arrays."""
    E, s, tau = np.broadcast_arrays(np.asarray(E, float),
                                    np.asarray(s, float),
                                    np.asarray(tau, float))
    R = np.full(E.shape, R0)
    P = np.full(E.shape, P0)
    K = np.full(E.shape, K0)
    keys = ("R", "K", "P", "ore", "inv", "food")
    hist = {k: [] for k in keys}
    for t in range(T):
        ore = np.minimum(np.minimum(E, R), CAP * K)       # extraction
        rev = PRICE * ore                                  # sales revenue
        inv = np.where(t < tau, s, 0.0) * rev              # reinvestment
        food = rev - inv                                   # food budget
        z = food / (np.maximum(P, 1e-9) * NEED)            # food sufficiency
        rate = 2.0 * GMAX * np.clip(z - 1.0, -1.0, 0.5)    # growth rate
        if record:
            for k, v in zip(keys, (R, K, P, ore, inv, food)):
                hist[k].append(float(v))
        R = R - ore
        K = (1.0 - DEP) * K + inv / MACH
        P = np.maximum(P * (1.0 + rate), 0.0)
    if record:
        for k, v in zip(("R", "K", "P"), (R, K, P)):
            hist[k].append(float(v))
        return float(P), {k: np.array(v) for k, v in hist.items()}
    return P


# ------------------------------------------------------------
# 2. Exhaustive search over (E, s, tau) in one vectorised call
# ------------------------------------------------------------
E_grid = np.linspace(10, 80, 71)
s_grid = np.linspace(0, 1, 51)
tau_grid = np.arange(0, T + 1, dtype=float)

t0 = time.time()
PF = simulate(E_grid[:, None, None], s_grid[None, :, None], tau_grid[None, None, :])
elapsed = time.time() - t0
iE, iS, iT = np.unravel_index(np.argmax(PF), PF.shape)


# ------------------------------------------------------------
# 3. Continuous refinement with Nelder-Mead
# ------------------------------------------------------------
def objective(x):
    E = np.clip(x[0], 1.0, 200.0)
    s = np.clip(x[1], 0.0, 1.0)
    tau = np.clip(x[2], 0.0, T)
    return -float(simulate(E, s, tau))


res = minimize(objective, [E_grid[iE], s_grid[iS], tau_grid[iT]],
               method="Nelder-Mead",
               options={"xatol": 1e-3, "fatol": 1e-4, "maxiter": 2000})
Eo = float(np.clip(res.x[0], 1.0, 200.0))
so = float(np.clip(res.x[1], 0.0, 1.0))
to = float(np.clip(res.x[2], 0.0, T))
Po, ho = simulate(Eo, so, to, record=True)

# ------------------------------------------------------------
# 4. Baseline policies for comparison
# ------------------------------------------------------------
# (a) never invest, feed everyone
a_E, _ = np.unravel_index(np.argmax(PF[:, 0, :]), PF[:, 0, :].shape)
Pa, ha = simulate(E_grid[a_E], 0.0, 0.0, record=True)
# (b) reinvest 80% of revenue for the whole horizon
j_s = int(np.argmin(np.abs(s_grid - 0.8)))
b_E = int(np.argmax(PF[:, j_s, T]))
Pb, hb = simulate(E_grid[b_E], s_grid[j_s], float(T), record=True)
# (c) over-mining: extract as fast as capacity allows (E = 80), best s and tau
c_E = len(E_grid) - 1
c_s, c_t = np.unravel_index(np.argmax(PF[c_E]), PF[c_E].shape)
Pc, hc = simulate(E_grid[c_E], s_grid[c_s], tau_grid[c_t], record=True)

# ------------------------------------------------------------
# 5. Console report
# ------------------------------------------------------------
print("=" * 70)
print(f"Scenarios evaluated : {PF.size:,} (vectorised, {elapsed:.3f} s)")
print(f"Grid optimum        : E={E_grid[iE]:.1f}, s={s_grid[iS]:.2f}, "
      f"tau={tau_grid[iT]:.0f}, P_T={PF[iE, iS, iT]:.1f}")
print(f"Refined optimum     : E={Eo:.2f} t/yr, s={so:.3f}, tau={to:.2f} yr")
print(f"                      P_T={Po:.1f}, ore left={ho['R'][-1]:.1f} t, "
      f"machines={ho['K'][-1]:.1f}")
print("-" * 70)
print(f"{'Policy':<34}{'Final pop.':>12}{'Ore left [t]':>14}")
for name, p, h in (("Optimal", Po, ho),
                   ("(a) No reinvestment", Pa, ha),
                   ("(b) Reinvest 80% always", Pb, hb),
                   ("(c) Over-mining (E=80)", Pc, hc)):
    print(f"{name:<34}{p:>12.1f}{h['R'][-1]:>14.1f}")
print("=" * 70)

# ------------------------------------------------------------
# 6. All figures in one image
# ------------------------------------------------------------
years = np.arange(T + 1)
yr = np.arange(T)
policies = ((ho, "Optimal", "tab:red"), (ha, "(a) No reinvestment", "tab:blue"),
            (hb, "(b) Reinvest 80%", "tab:green"), (hc, "(c) Over-mining", "tab:orange"))

fig = plt.figure(figsize=(21, 13))
fig.suptitle("Mining, Reinvestment and Final Population", fontsize=20, y=0.995)

ax1 = fig.add_subplot(2, 3, 1, projection="3d")
EE, SS = np.meshgrid(E_grid, s_grid, indexing="ij")
ax1.plot_surface(EE, SS, PF[:, :, iT], cmap=cm.viridis, linewidth=0, antialiased=True, alpha=0.95)
ax1.scatter([Eo], [so], [Po], color="red", s=80)
ax1.set_xlabel("Extraction E [t/yr]"); ax1.set_ylabel("Reinvestment share s")
ax1.set_zlabel("Final population")
ax1.set_title(f"3D surface: E vs s (tau = {tau_grid[iT]:.0f})")
ax1.view_init(elev=28, azim=-125)

ax2 = fig.add_subplot(2, 3, 2, projection="3d")
E2, T2 = np.meshgrid(E_grid, tau_grid, indexing="ij")
ax2.plot_surface(E2, T2, PF[:, iS, :], cmap=cm.plasma, linewidth=0, antialiased=True, alpha=0.95)
ax2.set_xlabel("Extraction E [t/yr]"); ax2.set_ylabel("Reinvestment years tau")
ax2.set_zlabel("Final population")
ax2.set_title(f"3D surface: E vs tau (s = {s_grid[iS]:.2f})")
ax2.view_init(elev=28, azim=-125)

ax3 = fig.add_subplot(2, 3, 3)
im = ax3.imshow(PF[:, :, iT].T, origin="lower", aspect="auto", cmap="viridis",
                extent=[E_grid[0], E_grid[-1], s_grid[0], s_grid[-1]])
ax3.plot(Eo, so, "r*", markersize=18, label="Optimum")
ax3.set_xlabel("Extraction E [t/yr]"); ax3.set_ylabel("Reinvestment share s")
ax3.set_title("Heatmap of final population"); ax3.legend()
fig.colorbar(im, ax=ax3, label="Final population")

ax4 = fig.add_subplot(2, 3, 4)
for h, label, col in policies:
    ax4.plot(years, h["R"], lw=2.2, label=label, color=col)
ax4.set_xlabel("Year"); ax4.set_ylabel("Remaining ore [t]")
ax4.set_title("Ore reserve"); ax4.grid(alpha=0.3); ax4.legend()

ax5 = fig.add_subplot(2, 3, 5)
for h, label, col in policies:
    ax5.plot(years, h["P"], lw=2.2, label=label, color=col)
ax5.set_xlabel("Year"); ax5.set_ylabel("Population")
ax5.set_title("Population trajectory"); ax5.grid(alpha=0.3); ax5.legend()

ax6 = fig.add_subplot(2, 3, 6)
ax6.bar(yr, ho["food"], color="tab:green", label="Food spending")
ax6.bar(yr, ho["inv"], bottom=ho["food"], color="tab:purple", label="Equipment spending")
ax6.set_xlabel("Year"); ax6.set_ylabel("Revenue allocation")
ax6.set_title("Optimal policy: how revenue is spent")
ax6.grid(alpha=0.3); ax6.legend()

plt.tight_layout()
plt.show()

5. Code Walkthrough

5.1 Parameters and the `simulate` function

The constants at the top mirror the table in Section 2. The core is simulate, which runs the 30-year economy forward in time.

np.broadcast_arrays makes E, s and tau share one shape. Passing arrays of shapes (71,1,1), (1,51,1) and (1,1,31) therefore evaluates the whole $71\times51\times31$ cube at once. This is the speed-up: the Python loop runs over the 30 years only, and each iteration is a NumPy operation over 112,251 scenarios.
ore = min(E, R, CAP*K) implements $e_t=\min(E,R_t,cK_t)$. The term R guarantees that no more ore than the deposit holds can ever be extracted, so the deposit constraint is always satisfied.
np.where(t < tau, s, 0.0) switches reinvestment off after year $\tau$, which is the “invest first, harvest later” structure.
The growth rate 2*GMAX*clip(z-1, -1, 0.5) gives $-10%$ at $z=0$, $0%$ at $z=1$, and saturates at $+5%$ for $z\ge1.5$. Food beyond $1.5$ times the need is wasted, which is why over-feeding is never optimal.
With record=True, the function stores the yearly trajectories of every state and returns them for plotting. Without it, only the final population array is returned, which keeps the grid search cheap.

5.2 Exhaustive search

PF is a three-dimensional array whose entry PF[i, j, k] is the final population for E_grid[i], s_grid[j] and tau_grid[k]. np.argmax followed by np.unravel_index recovers the best triple. Because the decision space is small, this brute-force approach is guaranteed to find the global optimum of the grid and, unlike a gradient-based method, is not troubled by the non-smooth min and clip operations.

The grid has finite resolution, so scipy.optimize.minimize with Nelder–Mead polishes the best grid point. Nelder–Mead needs no gradients, which suits this piecewise-defined model. The np.clip calls inside objective keep the search inside the feasible box $0\le s\le1$, $0\le\tau\le T$.

5.4 Baseline policies

The three baselines show what the optimization is worth:

(a) No reinvestment: the best $E$ with $s=0$. The colony never expands its machine fleet.
(b) Reinvest 80% always: the colony pours almost everything into equipment for the whole horizon and starves its own people.
(c) Over-mining: the extraction target is set to the maximum of the grid ($80$ t/year), and only $s$ and $\tau$ are optimized.

5.5 Console report

The report prints the number of scenarios and the time taken, the grid optimum, the refined optimum, and a comparison table of final population and ore left for all four policies.

6. Console Output

======================================================================
Scenarios evaluated : 112,251 (vectorised, 0.087 s)
Grid optimum        : E=44.0, s=0.46, tau=21, P_T=1080.0
Refined optimum     : E=45.35 t/yr, s=0.459, tau=21.00 yr
                      P_T=1082.0, ore left=0.5 t, machines=14.9
----------------------------------------------------------------------
Policy                              Final pop.  Ore left [t]
Optimal                                 1082.0           0.5
(a) No reinvestment                      296.6         686.9
(b) Reinvest 80% always                  287.0           1.7
(c) Over-mining (E=80)                  1077.5           0.8
======================================================================

How to Read the Console Output

In my test run, the whole $112{,}251$-scenario search finished in a fraction of a second, and the refined optimum was close to

$$
E^* \approx 45.4\ \text{t/year}, \qquad s^* \approx 0.46, \qquad \tau^* \approx 21\ \text{years}, \qquad P_T \approx 1{,}082 .
$$

Almost no ore is left at the end ($\approx 0.5$ t), so the deposit is used up almost exactly at the horizon. The machine fleet grows from 10 to about 15 units.

The comparison table shows the following:

No reinvestment ends at roughly $297$ people. The fleet decays through depreciation and the capacity $cK_t$ falls below the ore flow needed, so about $687$ t of ore stays in the ground unsold. This is the opposite failure to depletion: the resource is wasted.
Reinvesting 80% always ends at roughly $287$ people. The ore is used up, but the colony spent too much on machines, starved, and never recovered.
Over-mining with $E=80$ reaches about $1{,}077$, almost as good as the optimum. The reason is that the real throttle is the machine capacity, not the target $E$: once $E$ exceeds $cK_t$, raising it changes nothing. What separates good from bad policies is therefore the timing and size of the investment, not the extraction target alone.

7. Graphs

How to Read the Graphs

Top left: 3D surface of $E$ and $s$ (at the optimal $\tau$). The surface is a ridge rather than a sharp peak. For small $E$, the colony sells too little ore and the final population is low whatever $s$ is. For intermediate $E$ there is a clear maximum in $s$: with too little reinvestment, capacity is the bottleneck, and with too much, food runs short. The red dot marks the refined optimum sitting on top of the ridge.

Top center: 3D surface of $E$ and $\tau$ (at the optimal $s$). This shows the role of the investment window. Too short a window leaves the fleet too small to mine the whole deposit in time. Too long a window starves the colony because revenue goes to machines instead of food. The best $\tau$ lies in the second half of the horizon, so the colony invests for most of its life and switches to pure food spending near the end.

Top right: heatmap of final population. This is the top view of the first surface. The bright region is narrow along the $s$ axis and wide along the $E$ axis. In other words, the result is sensitive to the reinvestment share but forgiving about the extraction target, because capacity limits extraction anyway. The red star marks the optimum.

Bottom left: ore reserve. The optimal curve (red) descends steadily and touches zero right at year 30. The no-reinvestment curve (blue) stays high, which shows unused wealth. The over-mining curve (orange) drains the deposit fast, but only as fast as machine capacity allows, so it behaves much like the optimum.

Bottom center: population trajectory. The optimal population dips or stagnates early, because part of the revenue goes to machines, and then grows faster later as the larger fleet raises revenue. This is the classic invest now, benefit later pattern. The always-reinvest policy (green) collapses because the colony starves, and the no-reinvestment policy (blue) stalls because capacity decays.

Bottom right: allocation of revenue under the optimal policy. Each bar is one year’s revenue, split into food (green) and equipment (purple). The purple part is present for roughly the first 21 years and vanishes afterwards, which is exactly the switching structure of the optimum. The bars also grow over time as the machine fleet expands.

8. Conclusion

A finite deposit turns the colony’s problem into a timing problem: what matters is when to invest and when to switch to consumption.
The optimum extracts the deposit almost exactly as the horizon ends. It leaves no ore unused, and it does not exhaust the mine early and starve the people.
Both extremes are bad. Too little reinvestment wastes the resource, and too much starves the population.
Because the decision space is small, vectorized brute force plus a Nelder–Mead polish is fast and finds the global optimum without worrying about non-smooth model terms.

For richer variants, one can let the population supply the mining labor, make the ore price depend on extraction, or replace the single switching time $\tau$ with a fully time-dependent share $s_t$ optimized by dynamic programming.

Maximizing the Population After 15 Years with Linear Programming

October 23, 2026

How Should a Colony Split Its Budget Between Housing and Food?

Every growing settlement faces the same dilemma. Building houses raises the number of people the town can accept, but every unit of budget spent on construction is a unit not spent on farmland, irrigation, or storage. Houses also take years to complete, while food investment pays off right away. So how should the money be split each year to maximize the settled population 15 years from now?

In this article we model the problem as a linear program, solve it with SciPy, and compare the optimal plan against simple “fixed ratio” strategies. We also explore how the construction period and the budget scale change the answer, using 3D surfaces.

1. A Concrete Example

Consider a small colony with the following setting.

Item	Value
Planning horizon	15 years
Initial population	800
Initial housing capacity	1,000 people
Initial food capacity	1,000 people
Housing gained per budget unit	25 people
Food gained per budget unit	15 people
Base annual budget	6 units
Extra budget per person fed (tax base)	0.012 units
Maximum annual population growth	30 %
Construction period	3 years

Housing is a delayed investment: money spent in year $s$ only produces homes in year $s+L$. Food investment takes effect in the next year. Because the budget grows as more people are fed, early food investment also pays dividends later.

2. Mathematical Formulation

Let $h_t$ and $f_t$ be the budget spent on housing and food in year $t$ ($t=0,\dots,T-1$), and let $P_t$ be the population in year $t$.

Housing capacity (with construction period $L$):

$$
H_t = H_0 + c_h \sum_{s=0}^{t-L} h_s
$$

Food capacity:

$$
F_t = F_0 + c_f \sum_{s=0}^{t-1} f_s
$$

Annual budget, which grows with the number of people that can be fed:

$$
B_t = B_0 + g,F_t
$$

Budget constraint:

$$
h_t + f_t \le B_t
;;\Longleftrightarrow;;
h_t + f_t - g,c_f \sum_{s=0}^{t-1} f_s \le B_0 + g,F_0
$$

The right-hand form is linear in the decision variables, which is what makes this a linear program.

Population constraints:

$$
P_t \le H_t,\qquad P_t \le F_t,\qquad P_t \le (1+\gamma),P_{t-1}
$$

Objective:

$$
\max ; P_{T}
\quad \text{subject to} \quad h_t,\ f_t \ge 0
$$

The population can never exceed the smaller of housing and food capacity, so the optimal plan must balance the two at the end of the horizon while handling the construction delay.

3. Source Code

Everything, including the solver, the strategy comparison, the sensitivity analysis, and all graphs, is in a single cell.

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

# ---------------- Parameters ----------------
T = 15          # planning horizon (years)
P0 = 800.0      # initial population
H0 = 1000.0     # initial housing capacity (people)
F0 = 1000.0     # initial food capacity (people fed)
CH = 25.0       # people housed per budget unit
CF = 15.0       # people fed per budget unit
B0 = 6.0        # base annual budget
G = 0.012       # budget per person fed (tax base)
GAMMA = 0.30    # max annual population growth
L_BASE = 3      # construction period (years)


def solve_lp(L=L_BASE, scale=1.0):
    nh, nf, npop = T, T, T
    n = nh + nf + npop
    ih = lambda t: t
    jf = lambda t: nh + t
    kp = lambda t: nh + nf + (t - 1)
    A, b = [], []
    for t in range(T):  # budget constraint
        row = np.zeros(n)
        row[ih(t)] = 1.0
        row[jf(t)] = 1.0
        for s in range(t):
            row[jf(s)] -= scale * G * CF
        A.append(row)
        b.append(scale * (B0 + G * F0))
    for t in range(1, T + 1):  # housing capacity
        row = np.zeros(n)
        row[kp(t)] = 1.0
        for s in range(0, t - L + 1):
            row[ih(s)] -= CH
        A.append(row)
        b.append(H0)
    for t in range(1, T + 1):  # food capacity
        row = np.zeros(n)
        row[kp(t)] = 1.0
        for s in range(t):
            row[jf(s)] -= CF
        A.append(row)
        b.append(F0)
    for t in range(1, T + 1):  # growth limit
        row = np.zeros(n)
        row[kp(t)] = 1.0
        if t == 1:
            A.append(row)
            b.append((1 + GAMMA) * P0)
        else:
            row[kp(t - 1)] = -(1 + GAMMA)
            A.append(row)
            b.append(0.0)
    c = np.zeros(n)
    c[kp(T)] = -1.0
    res = linprog(c, A_ub=np.array(A), b_ub=np.array(b),
                  bounds=[(0, None)] * n, method="highs")
    if not res.success:
        raise RuntimeError(res.message)
    h = res.x[:nh]
    f = res.x[nh:nh + nf]
    return h, f, -res.fun


def simulate(L=L_BASE, scale=1.0, h_arr=None, f_arr=None, alpha=None):
    H = np.zeros(T + 1)
    Fd = np.zeros(T + 1)
    P = np.zeros(T + 1)
    B = np.zeros(T)
    h = np.zeros(T)
    f = np.zeros(T)
    P[0] = P0
    for t in range(T + 1):
        H[t] = H0 + CH * h[:max(0, t - L + 1)].sum()
        Fd[t] = F0 + CF * f[:t].sum()
        if t >= 1:
            P[t] = min(H[t], Fd[t], (1 + GAMMA) * P[t - 1])
        if t < T:
            B[t] = scale * (B0 + G * Fd[t])
            if h_arr is not None:
                h[t], f[t] = h_arr[t], f_arr[t]
            else:
                h[t], f[t] = alpha * B[t], (1 - alpha) * B[t]
    return H, Fd, P, B, h, f


# ---------------- Solve the base case ----------------
h_opt, f_opt, P_final_lp = solve_lp()
H, Fd, P, B, h_s, f_s = simulate(h_arr=h_opt, f_arr=f_opt)

print("=== Base case: optimal investment plan ===")
print(f"Maximum population after {T} years: {P[T]:,.1f} people")
print(f"Total budget spent on housing: {h_opt.sum():8.1f}")
print(f"Total budget spent on food   : {f_opt.sum():8.1f}")
print()
print("Year | Budget | Housing | Food  | Housing cap | Food cap | Population")
for t in range(T):
    print(f"{t:4d} | {B[t]:6.1f} | {h_opt[t]:7.1f} | {f_opt[t]:5.1f} | "
          f"{H[t]:11.0f} | {Fd[t]:8.0f} | {P[t]:10.0f}")
print(f"{T:4d} |    -   |    -    |   -   | {H[T]:11.0f} | {Fd[T]:8.0f} | {P[T]:10.0f}")

# ---------------- Fixed-share strategies ----------------
alphas = np.linspace(0, 1, 41)
fixed_pop = np.array([simulate(alpha=a)[2][T] for a in alphas])
best_i = int(np.argmax(fixed_pop))
print()
print("=== Fixed housing share strategies ===")
print(f"Best fixed share: {alphas[best_i]:.2f} -> population {fixed_pop[best_i]:,.1f}")
print(f"Optimal dynamic plan: population {P[T]:,.1f} "
      f"(+{P[T] - fixed_pop[best_i]:,.1f} vs best fixed share)")

# ---------------- Sensitivity: construction period x budget ----------------
Ls = np.arange(1, 7)
scales = np.linspace(0.6, 1.6, 11)
Z = np.zeros((len(Ls), len(scales)))
for i, L in enumerate(Ls):
    for j, sc in enumerate(scales):
        Z[i, j] = solve_lp(L=int(L), scale=float(sc))[2]
print()
print("=== Max population by construction period (rows) and budget scale (columns) ===")
print("L \\ scale " + " ".join(f"{s:6.2f}" for s in scales))
for i, L in enumerate(Ls):
    print(f"L={L}      " + " ".join(f"{Z[i, j]:6.0f}" for j in range(len(scales))))

# ---------------- Fixed share x construction period ----------------
alpha_grid = np.linspace(0, 1, 21)
Zf = np.zeros((len(Ls), len(alpha_grid)))
for i, L in enumerate(Ls):
    for j, a in enumerate(alpha_grid):
        Zf[i, j] = simulate(L=int(L), alpha=float(a))[2][T]

# ---------------- Figure ----------------
fig = plt.figure(figsize=(20, 12))
years = np.arange(T + 1)

ax1 = fig.add_subplot(2, 3, 1)
ax1.plot(years, H, "s--", label="Housing capacity")
ax1.plot(years, Fd, "^--", label="Food capacity")
ax1.plot(years, P, "o-", lw=3, color="k", label="Population")
ax1.set_xlabel("Year")
ax1.set_ylabel("People")
ax1.set_title("(a) Capacity and population under the optimal plan")
ax1.grid(alpha=0.3)
ax1.legend()

ax2 = fig.add_subplot(2, 3, 2)
ax2.bar(np.arange(T), f_opt, label="Food investment", color="tab:green")
ax2.bar(np.arange(T), h_opt, bottom=f_opt, label="Housing construction", color="tab:orange")
ax2.plot(np.arange(T), B, "k.-", label="Available budget")
ax2.set_xlabel("Year")
ax2.set_ylabel("Budget units")
ax2.set_title("(b) Optimal annual budget allocation")
ax2.grid(alpha=0.3)
ax2.legend()

ax3 = fig.add_subplot(2, 3, 3)
ax3.plot(alphas, fixed_pop, "o-", label="Fixed housing share")
ax3.axhline(P[T], color="r", ls="--", label="Optimal dynamic plan")
ax3.set_xlabel("Housing share of the budget")
ax3.set_ylabel(f"Population at year {T}")
ax3.set_title("(c) Fixed share vs optimal plan")
ax3.grid(alpha=0.3)
ax3.legend()

ax4 = fig.add_subplot(2, 3, 4, projection="3d")
SX, LY = np.meshgrid(scales, Ls)
surf = ax4.plot_surface(SX, LY, Z, cmap="viridis", edgecolor="k", linewidth=0.3)
ax4.set_xlabel("Budget scale")
ax4.set_ylabel("Construction period (years)")
ax4.set_zlabel("Population")
ax4.set_title("(d) 3D: optimal population")
ax4.view_init(elev=28, azim=-120)
fig.colorbar(surf, ax=ax4, shrink=0.6, pad=0.1)

ax5 = fig.add_subplot(2, 3, 5, projection="3d")
AX, LY2 = np.meshgrid(alpha_grid, Ls)
surf2 = ax5.plot_surface(AX, LY2, Zf, cmap="plasma", edgecolor="k", linewidth=0.3)
ax5.set_xlabel("Housing share")
ax5.set_ylabel("Construction period (years)")
ax5.set_zlabel("Population")
ax5.set_title("(e) 3D: fixed-share strategies")
ax5.view_init(elev=28, azim=-60)
fig.colorbar(surf2, ax=ax5, shrink=0.6, pad=0.1)

ax6 = fig.add_subplot(2, 3, 6)
im = ax6.imshow(Z, origin="lower", aspect="auto", cmap="viridis",
                extent=[scales[0], scales[-1], Ls[0] - 0.5, Ls[-1] + 0.5])
ax6.set_xlabel("Budget scale")
ax6.set_ylabel("Construction period (years)")
ax6.set_title("(f) Heatmap of (d)")
fig.colorbar(im, ax=ax6, label="Population")

plt.tight_layout()
plt.show()

4. Code Walkthrough

4.1 Parameters

The constants at the top mirror the table in Section 1. CH and CF are the “efficiency” of one budget unit in housing and food. G is the tax-base coefficient that links food capacity to the next year’s budget. GAMMA caps annual population growth, which prevents the model from jumping instantly to the final population.

4.2 Building the linear program: `solve_lp`

The decision vector has $3T$ entries, laid out as

Three small helper lambdas (ih, jf, kp) convert a year index into a position in this vector. Four blocks of inequality rows $A x \le b$ are then built:

Budget: $h_t + f_t - g,c_f\sum_{s<t} f_s \le B_0 + gF_0$. The inner loop subtracts the tax revenue created by earlier food investments. The scale argument multiplies the whole budget so we can run sensitivity studies.
Housing: $P_t - c_h\sum_{s\le t-L} h_s \le H_0$. The range range(0, t - L + 1) encodes the construction delay. If $t<L$ the range is empty, so no new homes exist yet.
Food: $P_t - c_f\sum_{s<t} f_s \le F_0$.
Growth: $P_1 \le (1+\gamma)P_0$ and $P_t - (1+\gamma)P_{t-1}\le 0$ for later years.

Since linprog minimizes, the objective vector c has $-1$ on $P_T$, and the sign is flipped back when returning -res.fun. The HiGHS solver is used, and a failure raises an explicit error instead of silently returning garbage.

4.3 Forward simulation: `simulate`

This function replays a plan year by year with plain arithmetic. It is used for two purposes:

To recompute $H_t$, $F_t$, $P_t$ and $B_t$ from the optimal $h_t, f_t$ so they can be tabulated and plotted.
To evaluate fixed-share strategies, where a constant fraction alpha of each year’s budget goes to housing and the rest to food.

The population update min(H[t], Fd[t], (1 + GAMMA) * P[t - 1]) is the exact counterpart of the three population constraints.

4.4 Sensitivity analysis and speed

The double loop calls solve_lp for 6 construction periods × 11 budget scales, which is 66 independent linear programs. Each LP has only 45 variables, and HiGHS solves it in milliseconds, so the whole analysis finishes within a few seconds. This is a major advantage of the linear formulation: a brute-force search over year-by-year allocations would face a continuous space of dimension 30 and would be hopelessly slow. Likewise, the fixed-share sweep is only $6 \times 21$ cheap simulations, so no further acceleration is needed.

4.5 Plotting

A single figure with six panels is drawn:

(a) housing capacity, food capacity and population over time
(b) stacked bars of the optimal allocation against the available budget
(c) fixed-share strategies compared with the optimal plan
(d) a 3D surface of the optimal population versus budget scale and construction period
(e) a 3D surface of fixed-share strategies versus housing share and construction period
(f) a heatmap view of (d)

5. Results

Console output

=== Base case: optimal investment plan ===
Maximum population after 15 years: 8,202.7 people
Total budget spent on housing:    288.1
Total budget spent on food   :    480.2

Year | Budget | Housing | Food  | Housing cap | Food cap | Population
   0 |   18.0 |     0.0 |  18.0 |        1000 |     1000 |        800
   1 |   21.2 |     0.0 |  21.2 |        1000 |     1270 |       1000
   2 |   25.1 |     0.0 |  25.1 |        1000 |     1589 |       1000
   3 |   29.6 |     0.0 |  29.6 |        1000 |     1965 |       1000
   4 |   34.9 |     0.2 |  34.7 |        1000 |     2408 |       1000
   5 |   41.1 |    12.1 |  29.1 |        1000 |     2928 |       1000
   6 |   46.4 |    15.7 |  30.7 |        1000 |     3364 |       1000
   7 |   51.9 |    20.4 |  31.5 |        1006 |     3825 |       1006
   8 |   57.6 |    26.5 |  31.1 |        1307 |     4297 |       1307
   9 |   63.2 |    34.5 |  28.7 |        1699 |     4763 |       1699
  10 |   68.3 |    44.8 |  23.5 |        2209 |     5193 |       2209
  11 |   72.6 |    59.0 |  13.6 |        2872 |     5546 |       2872
  12 |   75.0 |    75.0 |   0.0 |        3734 |     5750 |       3734
  13 |   75.0 |     0.0 |  75.0 |        4854 |     5750 |       4854
  14 |   88.5 |     0.0 |  88.5 |        6328 |     6875 |       6310
  15 |    -   |    -    |   -   |        8203 |     8203 |       8203

=== Fixed housing share strategies ===
Best fixed share: 0.43 -> population 5,800.3
Optimal dynamic plan: population 8,202.7 (+2,402.4 vs best fixed share)

=== Max population by construction period (rows) and budget scale (columns) ===
L \ scale   0.60   0.70   0.80   0.90   1.00   1.10   1.20   1.30   1.40   1.50   1.60
L=1        4188   5196   6393   7805   9455  11371  13573  16054  18848  21931  25289
L=2        4046   4985   6082   7354   8823  10505  12389  14504  16827  19352  22068
L=3        3894   4762   5759   6902   8203   9647  11257  13020  14925  16967  19113
L=4        3733   4523   5419   6437   7562   8809  10163  11612  13163  14781  16437
L=5        3560   4267   5058   5939   6898   7939   9044  10223  11443  12704  13969
L=6        3389   3994   4668   5396   6182   7002   7870   8752   9656  10618  11551

Graph image

6. Interpreting the Results

The optimal plan: food first, houses later, food again at the end

With a 3-year construction period, the optimal plan reaches a final population of about 8,203 people. It uses roughly 480 budget units on food and 288 on housing in total, and the timing is what matters.

Years 0–3: all-in on food. Housing capacity (1,000) is the binding limit, but the plan still spends everything on food. Why? Food capacity increases the tax base, so each early food unit generates extra budget for every subsequent year. This compounding is the reason to invest in food first, even though the population temporarily stalls at 1,000.
Years 4–12: a gradually rising housing share. Construction starts around year 5, and the housing share of the budget climbs steadily until year 12, when the entire budget goes into houses. Because homes take 3 years to appear, this spending shows up as capacity in years 8 through 15, and the population starts rising again from year 7.
Years 13–14: back to food. Housing money spent now would arrive after year 15 and be worthless for the objective, so the last two years switch entirely to food to push food capacity up to meet the housing capacity.

In panel (a), the housing and food capacity curves converge to the same value, about 8,203, at year 15. This is the signature of the optimum: any surplus on one side would be wasted budget, because population is limited by the smaller of the two.

Why a fixed ratio is not good enough

Panel (c) shows that the best fixed housing share is about 43 %, giving about 5,800 people. The optimal time-varying plan beats it by roughly 2,400 people (about 41 %). A fixed ratio cannot express “invest in food now, houses later, and food again at the end”. The delay and the tax-base feedback make the best ratio depend on time, and this is exactly what the linear program captures.

The 3D sensitivity surface

Panel (d) shows how the optimal population depends on the budget scale and the construction period.

Along the budget axis, the population rises monotonically. At $L=3$ it goes from about 3,900 at scale 0.6 to about 19,100 at scale 1.6. The growth is faster than linear because extra budget also compounds through the tax base.
Along the construction-period axis, a shorter period is always better. At the baseline budget (scale 1.0), the optimum falls from about 9,455 at $L=1$ to about 8,203 at $L=3$ and about 6,182 at $L=6$.
The two effects interact. At low budgets the construction period matters little in absolute terms, but at high budgets it matters a lot: at scale 1.6 the gap between $L=1$ and $L=6$ is more than 13,000 people. In other words, shortening construction is most valuable precisely when the budget is large.

The heatmap in panel (f) shows the same information from above. The contour-like bands slope diagonally, so one extra year of construction time costs about as much as a noticeable cut in the budget.

Panel (e) shows that for each construction period there is a single “ridge” in the housing share. As the construction period grows, the ridge shifts and the peak height drops, confirming that no single fixed ratio is robust across conditions, while the optimization automatically adapts.

7. Takeaways

Housing and food are not simple substitutes. Housing is a delayed investment, while food is an immediate investment that also expands the budget.
The optimal plan is a time-varying path: build the economic base first, ramp up housing, and spend the last years on whatever can still pay off before the deadline.
Writing the problem as a linear program makes it solvable exactly and instantly, and it makes sensitivity analysis over many scenarios essentially free.
The same framework can be extended with extra constraints such as maintenance costs, multiple building types, or random harvest shocks.

Optimal Water Allocation Over 20 Turns

October 22, 2026

Balancing Agriculture, Households, and Industry

Water is the one resource that every sector needs and no sector can replace. Farms need it to feed people, households need it to keep people alive, and factories turn it into the tax revenue that pays for tomorrow’s infrastructure. When a drought arrives, every cubic meter given to one sector is taken from another.

In this article we formulate a concrete water allocation problem as a 20-turn dynamic optimization, solve it numerically with Python, and visualize the result with a single combined figure that includes two 3D plots.

1. Problem Statement

A region receives a limited amount of water $W_t$ in each turn $t = 0, 1, \dots, 19$. The planner splits it among three sectors:

$$
\alpha_t^{\mathrm{agr}} + \alpha_t^{\mathrm{dom}} + \alpha_t^{\mathrm{ind}} = 1, \qquad \alpha_t^{\bullet} \ge 0 .
$$

The goal is to maximize the sum of the final population and the final asset value:

$$
\max_ ; J = P_{20} + w , A_{20}, \qquad w = 1 .
$$

Here $P$ is the population in thousands of people and $A$ is the asset value in millions of USD.

2. Mathematical Model

2.1 Water supply

The base supply drops sharply in a drought centered on turn 9. Infrastructure built from past tax revenue reduces leakage and raises the effective supply:

$$
W_t^{\mathrm{base}} = W_0\left(1 - 0.55, e^{-\left(\frac{t-9}{3}\right)^2}\right), \qquad
W_t = W_t^{\mathrm{base}}\left(1 + \eta_{\max}\frac{A_t}{A_t + A_{1/2}}\right).
$$

2.2 Satisfaction ratios

The water delivered to each sector is $W_t^{\mathrm{agr}} = \alpha_t^{\mathrm{agr}} W_t$, and similarly for the other two. The demand is proportional to the population, with per-capita coefficients $q_a$ and $q_d$:

$$
r_t^{\mathrm{agr}} = \frac{W_t^{\mathrm{agr}}}{q_a P_t}, \qquad
r_t^{\mathrm{dom}} = \frac{W_t^{\mathrm{dom}}}{q_d P_t}.
$$

The satisfaction level is a smooth version of $\min(1, r)$, which keeps the objective differentiable for gradient-based optimization:

$$
s(r) = 1 - \frac{1}{\beta}\log!\left(1 + e^{\beta (1-r)}\right).
$$

2.3 Population dynamics

Water shortage raises mortality, and a food shortage also suppresses births. The logistic factor represents the carrying capacity $K$:

$$
m_t = m_0 + m_d,(1 - s_t^{\mathrm{dom}})^2 + m_a,(1 - s_t^{\mathrm{agr}})^2,
$$

$$
P_{t+1} = P_t\left[1 + b, s_t^{\mathrm{agr}} s_t^{\mathrm{dom}}\left(1 - \frac{P_t}{K}\right) - m_t\right].
$$

2.4 Industry, tax revenue, and assets

Industrial output saturates with industrial water and scales with the square root of the workforce. A fraction $\kappa$ of the output is collected as tax and reinvested, while assets depreciate at rate $\delta$:

$$
Y_t = Y_{\max}\left(1 - e^{-W_t^{\mathrm{ind}}/k_i}\right)\sqrt{\frac{P_t}{P_0}}, \qquad
A_{t+1} = (1-\delta)A_t + \kappa Y_t .
$$

This creates the central trade-off. Industrial water increases future assets and future water efficiency, while household and agricultural water keeps people alive and growing.

2.5 Concrete example

Symbol	Value	Meaning
$P_0$, $A_0$	100, 200	Initial population (thousand) and asset value (million USD)
$W_0$	85	Base water supply per turn
$q_a$, $q_d$	0.50, 0.30	Per-capita agricultural and domestic demand
$K$, $b$	260, 0.09	Carrying capacity and birth coefficient
$m_0, m_d, m_a$	0.010, 0.25, 0.20	Mortality coefficients
$Y_{\max}$, $k_i$	60, 25	Industrial output scale and saturation constant
$\kappa$, $\delta$	0.5, 0.03	Reinvestment rate and depreciation
$\eta_{\max}$, $A_{1/2}$	0.4, 400	Infrastructure efficiency gain

3. Solution Approach

We compare five policies:

Equal split, with one third of the water to each sector.
Agriculture first, with shares (0.60, 0.30, 0.10).
Industry first, with shares (0.25, 0.20, 0.55).
Best constant, found by a grid search over all fixed allocations.
Optimized (dynamic), where the allocation can change every turn.

The dynamic policy has $20 \times 3 = 60$ logits $z_t$, mapped to valid shares by a softmax:

$$
\alpha_t^{k} = \frac{e^{z_t^{k}}}{\sum_{j} e^{z_t^{j}}} .
$$

The softmax removes the simplex constraint, so we can use L-BFGS-B directly. Because the objective is non-convex, we run six starting points and keep the best result. The simulator uses plain scalar arithmetic inside the 20-step loop, so the full optimization finishes in about a second and no further acceleration is needed.

4. Source Code

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

plt.style.use("dark_background")

# ----------------------------------------------------------------------
# 1. Model parameters
# ----------------------------------------------------------------------
T = 20
P0, A0 = 100.0, 200.0
W0 = 85.0
Q_A, Q_D = 0.50, 0.30
K_POP = 260.0
BIRTH = 0.09
M0, M_D, M_A = 0.010, 0.25, 0.20
Y_MAX, K_I = 60.0, 25.0
KAPPA, DELTA = 0.5, 0.03
EFF_MAX, EFF_HALF = 0.4, 400.0
WEIGHT = 1.0
BETA = 40.0

def base_supply(t):
    return W0 * (1.0 - 0.55 * math.exp(-((t - 9.0) / 3.0) ** 2))

def smooth_sat(r):
    x = 1.0 - r
    sp = np.logaddexp(0.0, BETA * x) / BETA
    return float(min(1.0, max(0.0, 1.0 - sp)))

def simulate(alloc):
    P, A = P0, A0
    Ps, As, Sa, Sd, Ws = [P], [A], [], [], []
    for t in range(T):
        eff = 1.0 + EFF_MAX * A / (A + EFF_HALF)
        W = base_supply(t) * eff
        wa, wd, wi = alloc[t, 0] * W, alloc[t, 1] * W, alloc[t, 2] * W
        s_a = smooth_sat(wa / (Q_A * P))
        s_d = smooth_sat(wd / (Q_D * P))
        mort = M0 + M_D * (1.0 - s_d) ** 2 + M_A * (1.0 - s_a) ** 2
        growth = BIRTH * s_a * s_d * max(0.0, 1.0 - P / K_POP)
        P_new = max(P * (1.0 + growth - mort), 1.0)
        out = Y_MAX * (1.0 - math.exp(-wi / K_I)) * math.sqrt(P / P0)
        A_new = (1.0 - DELTA) * A + KAPPA * out
        P, A = P_new, A_new
        Ps.append(P); As.append(A); Sa.append(s_a); Sd.append(s_d); Ws.append(W)
    J = P + WEIGHT * A
    return J, np.array(Ps), np.array(As), np.array(Sa), np.array(Sd), np.array(Ws)

def constant_policy(a, d, i):
    return np.tile(np.array([a, d, i], dtype=float), (T, 1))

def softmax_rows(z):
    z = z - z.max(axis=1, keepdims=True)
    e = np.exp(z)
    return e / e.sum(axis=1, keepdims=True)

def neg_objective(zflat):
    alloc = softmax_rows(zflat.reshape(T, 3))
    return -simulate(alloc)[0]

# ----------------------------------------------------------------------
# 2. Baseline policies and best constant policy
# ----------------------------------------------------------------------
policies = {
    "Equal split": constant_policy(1 / 3, 1 / 3, 1 / 3),
    "Agriculture first": constant_policy(0.60, 0.30, 0.10),
    "Industry first": constant_policy(0.25, 0.20, 0.55),
}

step = 0.025
grid = np.arange(0.0, 1.0 + 1e-9, step)
gx, gy, gz = [], [], []
best_const, best_const_J = None, -1e18
for a in grid:
    for i in grid:
        d = 1.0 - a - i
        if d < -1e-9:
            continue
        d = max(d, 0.0)
        J = simulate(constant_policy(a, d, i))[0]
        gx.append(a); gy.append(i); gz.append(J)
        if J > best_const_J:
            best_const_J, best_const = J, (a, d, i)
policies["Best constant"] = constant_policy(*best_const)

# ----------------------------------------------------------------------
# 3. Dynamic optimization (multi-start L-BFGS-B)
# ----------------------------------------------------------------------
t0 = time.time()
rng = np.random.default_rng(42)
starts = [np.zeros(T * 3)]
bc = np.log(np.maximum(np.array(best_const), 1e-3))
starts.append(np.tile(bc, T))
for _ in range(4):
    starts.append(rng.normal(0.0, 0.7, T * 3))

best_res = None
for z0 in starts:
    res = minimize(neg_objective, z0, method="L-BFGS-B",
                   bounds=[(-8.0, 8.0)] * (T * 3),
                   options={"maxiter": 300})
    if best_res is None or res.fun < best_res.fun:
        best_res = res
elapsed = time.time() - t0
policies["Optimized (dynamic)"] = softmax_rows(best_res.x.reshape(T, 3))

# ----------------------------------------------------------------------
# 4. Evaluation
# ----------------------------------------------------------------------
results = {name: simulate(al) for name, al in policies.items()}
print("Optimization time: %.2f s" % elapsed)
print("Best constant allocation (agri, domestic, industry): "
      "(%.3f, %.3f, %.3f)" % best_const)
print()
print("%-22s %10s %12s %10s" % ("Policy", "Population", "Asset value", "J"))
for name, r in results.items():
    print("%-22s %10.2f %12.2f %10.2f" % (name, r[1][-1], r[2][-1], r[0]))
print()
opt_alloc = policies["Optimized (dynamic)"]
print("Optimized allocation by turn (agri, domestic, industry):")
for t in range(T):
    print("turn %2d  W=%6.1f  (%.3f, %.3f, %.3f)" %
          (t + 1, results["Optimized (dynamic)"][5][t],
           opt_alloc[t, 0], opt_alloc[t, 1], opt_alloc[t, 2]))

# ----------------------------------------------------------------------
# 5. Visualization
# ----------------------------------------------------------------------
colors = {
    "Equal split": "#4cc9f0",
    "Agriculture first": "#90be6d",
    "Industry first": "#f9c74f",
    "Best constant": "#b388eb",
    "Optimized (dynamic)": "#ff4d6d",
}
turns = np.arange(1, T + 1)
fig = plt.figure(figsize=(20, 12))
gs = GridSpec(2, 3, figure=fig)

ax1 = fig.add_subplot(gs[0, 0])
ax1.stackplot(turns, opt_alloc[:, 0], opt_alloc[:, 1], opt_alloc[:, 2],
              labels=["Agriculture", "Domestic", "Industry"],
              colors=["#90be6d", "#4cc9f0", "#f9c74f"], alpha=0.9)
ax1.set_title("Optimized water allocation share")
ax1.set_xlabel("Turn"); ax1.set_ylabel("Share of available water")
ax1.set_xlim(1, T); ax1.set_ylim(0, 1)
ax1b = ax1.twinx()
ax1b.plot(turns, results["Optimized (dynamic)"][5], color="white", lw=2, ls="--",
          label="Effective supply")
ax1b.set_ylabel("Effective supply")
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax1b.get_legend_handles_labels()
ax1.legend(h1 + h2, l1 + l2, loc="lower left", fontsize=8)

ax2 = fig.add_subplot(gs[0, 1])
for name, r in results.items():
    ax2.plot(np.arange(0, T + 1), r[1], lw=2.5 if "Optimized" in name else 1.6,
             color=colors[name], label=name)
ax2.set_title("Population trajectory")
ax2.set_xlabel("Turn"); ax2.set_ylabel("Population (thousand)")
ax2.legend(fontsize=8); ax2.grid(alpha=0.2)

ax3 = fig.add_subplot(gs[0, 2])
for name, r in results.items():
    ax3.plot(np.arange(0, T + 1), r[2], lw=2.5 if "Optimized" in name else 1.6,
             color=colors[name], label=name)
ax3.set_title("Asset value trajectory")
ax3.set_xlabel("Turn"); ax3.set_ylabel("Asset value (million USD)")
ax3.legend(fontsize=8); ax3.grid(alpha=0.2)

ax4 = fig.add_subplot(gs[1, 0], projection="3d")
surf = ax4.plot_trisurf(np.array(gx), np.array(gy), np.array(gz),
                        cmap="plasma", linewidth=0.1, antialiased=True, alpha=0.95)
ax4.scatter([best_const[0]], [best_const[2]], [best_const_J], color="white",
            s=60, label="Best constant")
ax4.scatter([best_const[0]], [best_const[2]], [results["Optimized (dynamic)"][0]],
            color="#ff4d6d", s=80, marker="*", label="Optimized (dynamic)")
ax4.set_title("Objective J vs constant allocation")
ax4.set_xlabel("Agriculture share"); ax4.set_ylabel("Industry share")
ax4.set_zlabel("J")
ax4.view_init(elev=28, azim=-125)
ax4.legend(fontsize=8)
fig.colorbar(surf, ax=ax4, shrink=0.55, pad=0.08, label="J")

ax5 = fig.add_subplot(gs[1, 1], projection="3d")
for name, r in results.items():
    ax5.plot(r[1], r[2], np.arange(0, T + 1), lw=2.5 if "Optimized" in name else 1.5,
             color=colors[name], label=name)
    ax5.scatter([r[1][-1]], [r[2][-1]], [T], color=colors[name], s=30)
ax5.set_title("State-space trajectory")
ax5.set_xlabel("Population (thousand)"); ax5.set_ylabel("Asset value")
ax5.set_zlabel("Turn")
ax5.view_init(elev=22, azim=-60)
ax5.legend(fontsize=8)

ax6 = fig.add_subplot(gs[1, 2])
names = list(results.keys())
pops = [results[n][1][-1] for n in names]
assets = [WEIGHT * results[n][2][-1] for n in names]
xpos = np.arange(len(names))
ax6.bar(xpos, pops, color="#4cc9f0", label="Population")
ax6.bar(xpos, assets, bottom=pops, color="#f9c74f", label="Asset value")
for k, n in enumerate(names):
    ax6.text(k, pops[k] + assets[k] + 3, "%.1f" % results[n][0],
             ha="center", fontsize=9)
ax6.set_xticks(xpos)
ax6.set_xticklabels([n.replace(" ", "\n") for n in names], fontsize=8)
ax6.set_ylim(0, max(p + a for p, a in zip(pops, assets)) * 1.12)
ax6.set_title("Final objective J = population + asset value")
ax6.set_ylabel("J"); ax6.legend(fontsize=8); ax6.grid(alpha=0.2, axis="y")

fig.suptitle("Water allocation optimization over 20 turns", fontsize=18)
fig.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

5. Code Walkthrough

5.1 Parameters and helper functions

The constants in the first block correspond one-to-one with the table in Section 2.5. base_supply implements the Gaussian-shaped drought, which cuts the base supply by up to 55% around turn 9.

smooth_sat implements the smoothed satisfaction function $s(r)$. It uses np.logaddexp(0, x), which computes $\log(1+e^{x})$ without numerical overflow. The result is clipped to $[0, 1]$ to remove tiny floating-point excursions.

5.2 The simulator

simulate runs the 20-turn dynamics for any allocation matrix of shape (T, 3). In each turn it:

Computes the effective supply, including the infrastructure bonus from current assets.
Splits the supply into agricultural, domestic, and industrial water.
Converts each delivery into a satisfaction level.
Updates the population from the mortality and growth formulas.
Updates the assets from the industrial output.

It returns the objective $J$ together with all time series, which are reused later for plotting.

5.3 Softmax parametrization

softmax_rows converts a $20 \times 3$ matrix of logits into valid shares. Every row is non-negative and sums to one, so the optimizer works in an unconstrained space. Subtracting the row maximum before exponentiation is a standard stabilization step.

5.4 Baselines and the constant-policy grid

The three hand-made policies are the reference points. The double loop then scans every fixed allocation $(\alpha^{\mathrm{agr}}, \alpha^{\mathrm{ind}})$ on a grid of step 0.025, with $\alpha^{\mathrm{dom}} = 1 - \alpha^{\mathrm{agr}} - \alpha^{\mathrm{ind}}$. These roughly 800 evaluations serve two purposes. They give the best constant policy, and they provide the data for the 3D surface.

5.5 Dynamic optimization

The decision vector has 60 entries. We run L-BFGS-B from six starting points:

the uniform allocation,
the best constant allocation,
four random perturbations with a fixed seed.

The smooth satisfaction function matters here. A hard min would produce kinks that confuse gradient-based search, while the softplus version has a well-defined gradient everywhere. Bounds of $\pm 8$ on the logits keep the softmax numerically tame.

5.6 Evaluation and printing

Every policy is re-simulated, and the script prints a comparison table followed by the optimized allocation for each turn. The console output is the numerical counterpart of the figure.

5.7 The combined figure

All plots are drawn in a single $2 \times 3$ figure:

Top left: the optimized allocation as a stacked area, with the effective supply as a dashed line.
Top middle and right: population and asset trajectories for all five policies.
Bottom left (3D): the objective $J$ as a surface over constant allocations, with the best constant policy and the dynamic optimum marked.
Bottom middle (3D): the trajectory of each policy in (population, assets, turn) space.
Bottom right: a stacked bar chart that splits the final $J$ into its population and asset components.

6. Results

Combined figure

Console output

Optimization time: 4.39 s
Best constant allocation (agri, domestic, industry): (0.350, 0.250, 0.400)

Policy                 Population  Asset value          J
Equal split                 90.11       390.24     480.36
Agriculture first          115.21       238.17     353.38
Industry first              68.17       415.54     483.71
Best constant               89.24       410.27     499.50
Optimized (dynamic)         90.79       414.54     505.33

Optimized allocation by turn (agri, domestic, industry):
turn  1  W=  96.3  (0.308, 0.254, 0.438)
turn  2  W=  97.0  (0.309, 0.250, 0.441)
turn  3  W=  97.4  (0.309, 0.246, 0.445)
turn  4  W=  97.2  (0.307, 0.244, 0.449)
turn  5  W=  95.3  (0.303, 0.244, 0.454)
turn  6  W=  89.9  (0.295, 0.247, 0.458)
turn  7  W=  79.4  (0.282, 0.258, 0.460)
turn  8  W=  64.7  (0.267, 0.280, 0.453)
turn  9  W=  50.9  (0.272, 0.313, 0.415)
turn 10  W=  45.1  (0.324, 0.333, 0.342)
turn 11  W=  51.0  (0.375, 0.310, 0.314)
turn 12  W=  65.0  (0.379, 0.268, 0.353)
turn 13  W=  80.2  (0.362, 0.235, 0.403)
turn 14  W=  91.3  (0.348, 0.219, 0.433)
turn 15  W=  97.5  (0.344, 0.216, 0.440)
turn 16  W= 100.1  (0.347, 0.219, 0.433)
turn 17  W= 101.2  (0.353, 0.225, 0.422)
turn 18  W= 101.6  (0.358, 0.232, 0.410)
turn 19  W= 101.9  (0.362, 0.239, 0.399)
turn 20  W= 102.1  (0.367, 0.245, 0.388)

7. Discussion

7.1 The one-sided policies fail in opposite ways

The agriculture-first policy keeps the population large, since farmers and households are almost never short of water. But it starves industry, so the asset value barely moves from its starting level and the final $J$ is by far the lowest.

The industry-first policy shows the opposite failure. Assets grow quickly, but during the drought the population falls to its lowest level. Because industrial output scales with $\sqrt{P_t}$, the shrinking workforce eventually drags on the very thing the policy was trying to maximize.

7.2 The 3D surface has a single broad peak

The surface of $J$ over constant allocations is smooth and unimodal. Its maximum lies at roughly 35% agriculture, 25% domestic use, and 40% industry. The surface falls off steeply when the agricultural share drops very low, because that region produces famine-level mortality, and also when the industrial share approaches zero, because the asset engine stalls. The practical lesson is that a moderately balanced allocation already captures most of the achievable value.

7.3 Dynamic reallocation adds a further gain

The optimized policy beats the best constant allocation. The stacked area chart shows why:

Before the drought, the optimizer gives industry its largest share. Early assets compound through reinvestment and also raise the effective water supply, so water spent on industry pays back twice.
During the drought, industry is cut back while the domestic and agricultural shares rise. The effective supply is at its lowest, and protecting the population from the quadratic mortality penalty is worth more than the lost output.
After the drought, the industrial share recovers and the allocation drifts back toward its pre-drought pattern.

In other words, the optimizer treats industry as an investment to be front-loaded and treats household and farm water as insurance to be strengthened exactly when the supply is scarce.

7.4 State-space view

In the 3D trajectory plot, every policy starts at the same point. The agriculture-first path moves toward high population and low assets, and the industry-first path moves toward high assets and low population. The optimized path stays close to the efficient frontier of the other policies. It accepts a temporary population dip during the drought but ends with the highest combined value.

7.5 Final composition of the objective

The stacked bars make the trade-off explicit. Populations and asset values sit on different scales, so the weight $w$ matters: increasing $w$ shifts the optimum toward industry, and decreasing it shifts the optimum toward households and agriculture. Because the code exposes WEIGHT as a single constant, the policy preference can be explored by changing one line.

8. Conclusion

We modeled the allocation of scarce water as a 20-turn dynamic optimization problem with three coupled effects: shortage-driven mortality, industrial tax revenue that compounds into assets, and infrastructure that improves the effective supply. A smooth satisfaction function made the problem suitable for gradient-based optimization, and a softmax parametrization removed the allocation constraints.

The main findings were:

Extreme policies are dominated. Neglecting either industry or agriculture sharply lowers the final objective.
A well-chosen constant split already performs well, and the objective surface around it is smooth and forgiving.
A time-varying policy performs better still. Invest in industry while water is plentiful, then shift water to people and food during the drought.

The same framework extends naturally to stochastic droughts, multiple regions sharing a river, or water-saving technology investments as an additional decision variable.

Optimizing Farmland vs. Housing with Python

October 21, 2026

How Many People Can a Village Feed and House?

A village of 100 people owns a fixed amount of land. Every hectare turned into a field produces food, and every hectare turned into housing raises the population cap. Ten years from now, how many people can the village support, and what split of land gets it there?

This is a small problem, but it shows how two limiting factors, food and housing, pull against each other. In this post I model the village, solve it by brute-force search, and plot the result in 2D and 3D.

1. The Problem

The village has:

Item	Value
Total land $L$	60 ha
Initial population $P_0$	100 people
Planning horizon	10 years
Farming productivity $A$	8.0
Diminishing-returns exponent $\alpha$	0.85
Housing density $D$	10 people/ha
Maximum growth rate $r$	0.35 per year
Starvation severity $m$	0.60
Food needed per person $c$	1.0 per year

The decision variable is $x$, the hectares of farmland. The remaining $L - x$ hectares become residential land. The goal is to maximize the population after 10 years.

2. The Model

Food production shows diminishing returns, since extra fields need more labor, are less fertile and so on:

$$
F(x) = A,x^{\alpha}
$$

Housing capacity is linear in residential land:

$$
H(x) = D,(L - x)
$$

Each year the population grows logistically toward the housing capacity, and it shrinks if food runs short. The food sufficiency ratio is

$$
q_t = \frac{F(x)}{c,P_t}, \qquad s_t = \max\bigl(0,\ 1 - q_t\bigr)
$$

and the population update is

$$
P_{t+1} = P_t + r,P_t\left(1 - \frac{P_t}{H(x)}\right) - m,P_t,s_t
$$

The optimization problem is

$$
\max_{0 \le x \le L}\ P_{10}(x)
$$

Why is there an interior optimum? Too little farmland starves the village. Too little housing caps the population. At a long-run equilibrium with food shortage, the growth term and the starvation term balance:

so the best split balances the two constraints rather than maximizing either one.

3. The Complete Source Code

Everything is in one script. The search evaluates thousands of candidate land splits simultaneously with NumPy array operations instead of a Python loop over candidates.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401  (registers the 3D projection)

# ------------------------------------------------------------
# 1. Parameters of the village
# ------------------------------------------------------------
L = 60.0        # total usable land [ha]
P0 = 100.0      # initial population
YEARS = 10      # planning horizon
A = 8.0         # farming productivity coefficient
ALPHA = 0.85    # diminishing returns on farmland (0 < ALPHA <= 1)
D = 10.0        # people housed per hectare of residential land
R = 0.35        # maximum annual growth rate
M = 0.60        # starvation severity
C = 1.0         # food needed per person per year


# ------------------------------------------------------------
# 2. Vectorised simulator (x, A, D may be scalars or arrays)
# ------------------------------------------------------------
def simulate(x, A=A, D=D):
    """Return the population history, shape (YEARS + 1, *x.shape)."""
    x, A_, D_ = np.broadcast_arrays(
        np.asarray(x, dtype=float),
        np.asarray(A, dtype=float),
        np.asarray(D, dtype=float),
    )
    food = A_ * x ** ALPHA                       # annual food production
    cap = np.maximum(D_ * (L - x), 1e-9)         # housing capacity
    P = np.full(x.shape, P0)
    hist = [P.copy()]
    for _ in range(YEARS):
        q = food / (C * np.maximum(P, 1e-9))     # food sufficiency ratio
        shortage = np.clip(1.0 - q, 0.0, 1.0)
        P = P + R * P * (1.0 - P / cap) - M * P * shortage
        P = np.maximum(P, 0.0)
        hist.append(P.copy())
    return np.array(hist)


# ------------------------------------------------------------
# 3. Optimisation by exhaustive search on a fine grid
# ------------------------------------------------------------
x_grid = np.linspace(0.0, L, 6001)               # step = 0.01 ha
final_pop = simulate(x_grid)[-1]
i_best = int(np.argmax(final_pop))
x_best, p_best = x_grid[i_best], final_pop[i_best]

print("=" * 60)
print("Village land-allocation problem")
print("=" * 60)
print(f"Total land L = {L:.0f} ha, initial population = {P0:.0f}, horizon = {YEARS} years")
print(f"Optimal farmland   : {x_best:.2f} ha ({x_best / L:.1%} of the land)")
print(f"Optimal housing    : {L - x_best:.2f} ha")
print(f"Population in year {YEARS}: {p_best:.2f}")
print("-" * 60)
print(f"{'Farmland [ha]':>14} {'Housing [ha]':>13} {'Pop. (year 10)':>15}")
for xv in [10, 15, 20, 25, 30, 35, 40, 45, 50]:
    pv = simulate(float(xv))[-1]
    print(f"{xv:>14.0f} {L - xv:>13.0f} {float(pv):>15.2f}")
print("-" * 60)

best_hist = simulate(x_best)
food_best = A * x_best ** ALPHA
cap_best = D * (L - x_best)
print(f"Food production at optimum : {food_best:.2f} (feeds {food_best / C:.1f} people)")
print(f"Housing capacity at optimum: {cap_best:.2f} people")
print("Year-by-year population at the optimum:")
for t, p in enumerate(best_hist):
    print(f"  year {t:>2}: {float(p):7.2f}")

# ------------------------------------------------------------
# 4. Data for the figures
# ------------------------------------------------------------
years = np.arange(YEARS + 1)
sample_x = [10, 20, 30, float(np.round(x_best, 2)), 45, 50]
x_3d = np.linspace(0.0, L - 1.0, 70)
T3, X3 = np.meshgrid(years, x_3d)
Z3 = simulate(x_3d)[:, :].T                       # shape (len(x_3d), YEARS + 1)

a_vals = np.linspace(4.0, 12.0, 60)
XX, AA = np.meshgrid(np.linspace(0.0, L - 1.0, 60), a_vals)
ZZ = simulate(XX, A=AA)[-1]
best_x_for_a = XX[0][np.argmax(ZZ, axis=1)]

# ------------------------------------------------------------
# 5. One combined figure (2 x 3 panels)
# ------------------------------------------------------------
fig = plt.figure(figsize=(21, 12))

# (a) final population vs farmland
ax1 = fig.add_subplot(2, 3, 1)
ax1.plot(x_grid, final_pop, color="tab:green", lw=2)
ax1.axvline(x_best, color="red", ls="--", label=f"optimum: {x_best:.2f} ha")
ax1.scatter([x_best], [p_best], color="red", zorder=5)
ax1.set_xlabel("Farmland [ha]")
ax1.set_ylabel("Population in year 10")
ax1.set_title("(a) Final population vs. farmland")
ax1.grid(alpha=0.3)
ax1.legend()

# (b) trajectories
ax2 = fig.add_subplot(2, 3, 2)
for xv in sample_x:
    ax2.plot(years, simulate(float(xv)), marker="o", ms=4, label=f"farmland = {xv:g} ha")
ax2.set_xlabel("Year")
ax2.set_ylabel("Population")
ax2.set_title("(b) Population trajectories")
ax2.grid(alpha=0.3)
ax2.legend(fontsize=8)

# (c) the two limiting factors at year 0 conditions
ax3 = fig.add_subplot(2, 3, 3)
food_curve = A * x_grid ** ALPHA / C
house_curve = D * (L - x_grid)
ax3.plot(x_grid, food_curve, color="tab:orange", lw=2, label="People fed by farmland")
ax3.plot(x_grid, house_curve, color="tab:blue", lw=2, label="People housed by residential land")
ax3.plot(x_grid, np.minimum(food_curve, house_curve), color="black", ls=":", lw=2,
         label="Bottleneck = min(food, housing)")
ax3.axvline(x_best, color="red", ls="--")
ax3.set_xlabel("Farmland [ha]")
ax3.set_ylabel("People")
ax3.set_title("(c) Food limit vs. housing limit")
ax3.grid(alpha=0.3)
ax3.legend(fontsize=8)

# (d) 3D surface: time x farmland -> population
ax4 = fig.add_subplot(2, 3, 4, projection="3d")
ax4.plot_surface(T3, X3, Z3, cmap=cm.viridis, edgecolor="none", alpha=0.95)
ax4.set_xlabel("Year")
ax4.set_ylabel("Farmland [ha]")
ax4.set_zlabel("Population")
ax4.set_title("(d) 3D: population over time and farmland")
ax4.view_init(elev=28, azim=-130)

# (e) 3D surface: farmland x productivity -> final population
ax5 = fig.add_subplot(2, 3, 5, projection="3d")
ax5.plot_surface(XX, AA, ZZ, cmap=cm.plasma, edgecolor="none", alpha=0.95)
ax5.set_xlabel("Farmland [ha]")
ax5.set_ylabel("Productivity A")
ax5.set_zlabel("Population in year 10")
ax5.set_title("(e) 3D: final population vs. farmland and productivity")
ax5.view_init(elev=28, azim=-125)

# (f) heat map with the optimal ridge
ax6 = fig.add_subplot(2, 3, 6)
cs = ax6.contourf(XX, AA, ZZ, levels=25, cmap="plasma")
fig.colorbar(cs, ax=ax6, label="Population in year 10")
ax6.plot(best_x_for_a, a_vals, color="white", lw=2.5, label="Optimal farmland")
ax6.set_xlabel("Farmland [ha]")
ax6.set_ylabel("Productivity A")
ax6.set_title("(f) Optimal farmland as productivity changes")
ax6.legend(loc="upper right")

plt.suptitle("Farmland vs. Housing: Maximizing the Village Population After 10 Years",
             fontsize=17, y=1.0)
plt.tight_layout()
plt.show()

4. Code Walkthrough

Section 1: Parameters

All constants of the model are defined at the top, so you can experiment by changing a single number. ALPHA < 1 encodes diminishing returns on farmland. M controls how brutally the population shrinks when food runs short.

Section 2: The vectorised simulator

simulate(x, A, D) implements the update rule

$$
P_{t+1} = P_t + r,P_t\left(1 - \frac{P_t}{H}\right) - m,P_t,s_t
$$

Key details:

np.broadcast_arrays gives x, A and D a common shape. You can pass a single farmland value, an array of 6,001 candidates, or a 2D grid of farmland × productivity, and the same function handles all of them.
np.maximum(..., 1e-9) guards against division by zero. When farmland is 0 hectares, food is 0, and when it is 60 hectares, housing is 0. Without the guard you would get NaN values or warnings.
np.clip(1.0 - q, 0.0, 1.0) implements $s_t = \max(0, 1-q_t)$ and also caps the shortage at 100%, so the population cannot drop by more than $m$ in a year.
P = np.maximum(P, 0.0) prevents negative populations.
The loop runs only 10 times (one per year). The “loop” over candidate land splits is handled by NumPy internally.

Section 3: Optimisation

The search space is one-dimensional, so exhaustive search is both simple and reliable. np.linspace(0, L, 6001) creates 6,001 candidates with 0.01 ha resolution. A single call to simulate evaluates them all at once, and np.argmax picks the best. A scalar optimizer could be fooled by a flat plateau or an unusual shape, while a dense grid cannot.

The code then prints a summary, a table of selected land splits, and the year-by-year trajectory at the optimum. It also prints the food production and housing capacity at the optimum, which shows which constraint is binding.

Section 4: Data for the figures

sample_x contains the farmland values whose trajectories will be compared in panel (b).
T3, X3, Z3 build the grid for the first 3D surface. simulate(x_3d) returns an array of shape (11, 70) (years × candidates), so it is transposed to (70, 11) to match the meshgrid.
XX, AA, ZZ build the second surface by evaluating the final population over a grid of farmland and productivity A. This works only because simulate accepts array-valued A.
best_x_for_a finds, for each productivity level, the farmland that maximizes the final population, which gives the “optimal ridge” for the heat map.

Section 5: One combined figure

All six panels are drawn into a single fig with a 2 × 3 layout, so one plt.show() displays everything:

(a) the objective function $P_{10}(x)$ with the optimum marked,
(b) population trajectories for different land splits,
(c) the food limit versus the housing limit,
(d) a 3D surface of population over time and farmland,
(e) a 3D surface of final population over farmland and productivity,
(f) a heat map of (e) with the optimal ridge overlaid.

5. Console Output

============================================================
Village land-allocation problem
============================================================
Total land L = 60 ha, initial population = 100, horizon = 10 years
Optimal farmland   : 36.81 ha (61.4% of the land)
Optimal housing    : 23.19 ha
Population in year 10: 191.08
------------------------------------------------------------
 Farmland [ha]  Housing [ha]  Pop. (year 10)
            10            50          105.00
            15            45          135.00
            20            40          157.78
            25            35          174.43
            30            30          185.40
            35            25          190.66
            40            20          189.70
            45            15          148.77
            50            10          100.00
------------------------------------------------------------
Food production at optimum : 171.46 (feeds 171.5 people)
Housing capacity at optimum: 231.90 people
Year-by-year population at the optimum:
  year  0:  100.00
  year  1:  119.91
  year  2:  140.17
  year  3:  159.58
  year  4:  177.00
  year  5:  188.34
  year  6:  190.59
  year  7:  191.00
  year  8:  191.06
  year  9:  191.08
  year 10:  191.08

Reading the numbers

The search finds that the best allocation is about 36.8 hectares of farmland (61.4% of the land) and 23.2 hectares of housing, giving a population of about 191 people after 10 years, nearly double the starting 100.

The comparison table shows the shape of the trade-off:

With 10 ha of farmland, food is the bottleneck and the population barely moves (about 105).
With 50 ha of farmland, only 10 ha remain for housing. The cap of $10 \times 10 = 100$ people is exactly the starting population, so the village never grows (100.00).
The best results sit in a narrow band around 35–40 ha, where the final population is about 190.

At the optimum the village produces food for about 171 people but has housing for about 232. The binding constraint is food, but housing must be large enough that it does not choke growth. Substituting the equilibrium population $P^{*} \approx 191$ into the balance condition from Section 2 gives about $0.0616$ on both sides, so the simulation has settled at the equilibrium predicted by the formula.

The year-by-year trajectory rises quickly at first (100 → 120 → 140 → 160 → 177 → 188) and then flattens out around 191 from year 6 on. In other words, the village reaches its carrying capacity well before the ten years are up, and the optimization is really about choosing which carrying capacity it reaches.

6. Graph Output

Reading the graphs

(a) Final population vs. farmland. The curve rises steeply while food is scarce, peaks near 37 ha, and then falls sharply as housing shrinks. The drop on the right is much steeper than the climb on the left, which shows that over-farming is more dangerous than under-farming: housing vanishes linearly, and a population cannot live in fields it has no room to inhabit.

(b) Population trajectories. Each line is a different land split. Splits near the optimum (30–45 ha) all climb toward about 190. Small farmland shares flatten at low levels because of famine. The 50 ha line stays flat at 100 because the housing cap equals the starting population.

(c) Food limit vs. housing limit. The orange curve (people the farmland can feed) rises with diminishing returns. The blue line (people the housing can hold) falls linearly. The dotted black curve is the minimum of the two, which is the theoretical ceiling. The two curves cross at roughly 41 ha, but the red dashed optimum sits to the left of that crossing, at about 37 ha. The reason is that a population does not suddenly stop at its ceiling. Near the housing cap, logistic growth slows down, so giving up a little food to gain more room is worth it.

(d) 3D surface: time × farmland. The ridge of high population runs along the farmland axis near 35–40 ha and grows with time, then flattens by year 6 or so. The cliff on the high-farmland side marks housing exhaustion, and the low plateau on the low-farmland side marks famine. This is the whole dynamic of the problem in one picture.

(e) 3D surface: farmland × productivity. As productivity $A$ rises, the whole surface lifts and the peak moves toward less farmland. The land needed to feed a given number of people shrinks when each hectare is more productive, so the freed land can become housing.

(f) Heat map with optimal ridge. The white line is the optimal farmland as a function of productivity. It bends toward smaller farmland as $A$ grows, confirming the 3D picture. Better farming technology is not a reason to farm more land. It is a reason to build more houses.

7. Conclusion

With the 60 ha village in this example, the best strategy is to devote roughly 61% of the land to farming and 39% to housing, which brings the population from 100 to about 191 in ten years. The central lesson is that the answer is neither “maximize food” nor “maximize housing”. Growth is capped by the scarcer of the two resources, and the optimum balances them, tilted slightly toward food because starvation is more punishing than crowding. As farming becomes more productive, the best allocation shifts toward housing.

Designing Optimal Probes in Quantum Metrology

October 20, 2026

Maximizing Quantum Fisher Information with Python

How precisely can you measure a phase with $N$ qubits? The answer depends on two design choices: which input state you prepare and which measurement you perform. In this post we work through a concrete example, an $N$-qubit phase estimator, and solve the probe design problem numerically, both without noise and under local dephasing.

1. The Problem

A parameter $\varphi$ is imprinted on a probe state through the unitary

$$
U_\varphi = e^{-i\varphi G}, \qquad G=\frac12\sum_{i=1}^{N}\sigma_z^{(i)} .
$$

The quantum Cramér–Rao bound limits the variance of any unbiased estimator after $\nu$ repetitions:

$$
\mathrm{Var}(\hat\varphi);\ge;\frac{1}{\nu,F_Q[\rho_\varphi]} .
$$

The quantum Fisher information (QFI) $F_Q$ depends only on the probe state. For a pure state it is

$$
F_Q = 4\left(\langle G^2\rangle-\langle G\rangle^2\right).
$$

For a mixed state $\rho=\sum_i\lambda_i|i\rangle\langle i|$ it becomes

$$
F_Q = 2\sum_{i,j}\frac{(\lambda_i-\lambda_j)^2}{\lambda_i+\lambda_j},\big|\langle i|G|j\rangle\big|^2 .
$$

Step 1: Choosing the input state

The eigenvalues of $G$ range from $-N/2$ to $+N/2$, so the variance is bounded by $\mathrm{Var}(G)\le N^2/4$ and therefore

$$
F_Q\le N^2 \quad(\text{Heisenberg limit}).
$$

A product state such as $|+\rangle^{\otimes N}$ only reaches $F_Q=N$ (the standard quantum limit). The state that saturates the bound is the GHZ state

$$
|\mathrm{GHZ}\rangle=\frac{|0\rangle^{\otimes N}+|1\rangle^{\otimes N}}{\sqrt2}.
$$

Step 2: Choosing the measurement

The QFI is an upper bound. A concrete measurement with outcome probabilities $p_k(\varphi)$ yields the classical Fisher information

$$
F_C(\varphi)=\sum_k\frac{\big(\partial_\varphi p_k\big)^2}{p_k};\le;F_Q .
$$

We measure the same tilted Pauli observable on every qubit,

$$
\sigma_\vartheta=\sin\vartheta,\sigma_x+\cos\vartheta,\sigma_z ,
$$

where $\vartheta=0$ is the computational basis (blind to the phase of a GHZ state) and $\vartheta=\pi/2$ is the $x$-basis.

Step 3: Adding noise

Independent dephasing with probability $p$ on each qubit multiplies every density-matrix element by a factor that depends on the Hamming distance $d_H(i,j)$ between the two bit strings:

$$
\rho_{ij};\longrightarrow;(1-2p)^{,d_H(i,j)},\rho_{ij}.
$$

Because dephasing commutes with $U_\varphi$, we can compute the QFI of the noisy state directly. For the two reference probes the closed forms are

$$
F_Q^{\mathrm{GHZ}}=N^2(1-2p)^{2N},\qquad F_Q^{\mathrm{prod}}=N(1-2p)^2 .
$$

GHZ is better only while $(1-2p)^{2N-2}>1/N$. For $N=3$ the crossover is at $p=\tfrac12\left(1-3^{-1/4}\right)\approx 0.120$. Beyond that point the optimal probe must be found numerically.

2. Full Python Code

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

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(7)

# ------------------------------------------------------------------
# 1. Basic building blocks
# ------------------------------------------------------------------
def popcount_table(N):
    return np.array([bin(k).count("1") for k in range(2**N)])

def generator_diag(N):
    # G = (1/2) * sum_i Z_i is diagonal; |0> -> +1/2, |1> -> -1/2
    return 0.5 * (N - 2 * popcount_table(N))

def dephasing_kernel(N, p):
    # Local dephasing: rho_ij -> rho_ij * (1-2p)^{Hamming distance(i, j)}
    pc = popcount_table(N)
    idx = np.arange(2**N)
    ham = pc[idx[:, None] ^ idx[None, :]]
    return (1.0 - 2.0 * p) ** ham

def pure_qfi(psi, g):
    prob = np.abs(psi) ** 2
    return 4.0 * (np.sum(prob * g**2) - np.sum(prob * g) ** 2)

def mixed_qfi(rho, g):
    lam, U = np.linalg.eigh(rho)
    Gm = U.conj().T @ (g[:, None] * U)
    S = lam[:, None] + lam[None, :]
    D = lam[:, None] - lam[None, :]
    W = np.zeros_like(S)
    mask = S > 1e-12
    W[mask] = D[mask] ** 2 / S[mask]
    return 2.0 * np.sum(W * np.abs(Gm) ** 2)

def ghz_state(N):
    psi = np.zeros(2**N, dtype=complex)
    psi[0] = psi[-1] = 1 / np.sqrt(2)
    return psi

def product_state(N):
    return np.ones(2**N, dtype=complex) / np.sqrt(2**N)

def to_state(x):
    d = len(x) // 2
    psi = x[:d] + 1j * x[d:]
    return psi / np.linalg.norm(psi)

# ------------------------------------------------------------------
# 2. Numerical search for the optimal probe state
# ------------------------------------------------------------------
def optimize_pure(N, n_restarts=4):
    g = generator_diag(N)
    best_val, best_psi = -1.0, None
    for _ in range(n_restarts):
        x0 = rng.normal(size=2 * 2**N)
        res = minimize(lambda x: -pure_qfi(to_state(x), g), x0, method="L-BFGS-B")
        if -res.fun > best_val:
            best_val, best_psi = -res.fun, to_state(res.x)
    return best_val, best_psi

def optimize_noisy(N, p, n_restarts=3, x_warm=None):
    g = generator_diag(N)
    K = dephasing_kernel(N, p)
    def objective(x):
        psi = to_state(x)
        rho = np.outer(psi, psi.conj()) * K
        return -mixed_qfi(rho, g)
    starts = [rng.normal(size=2 * 2**N) for _ in range(n_restarts)]
    if x_warm is not None:
        starts.append(x_warm)
    best_val, best_x = -1.0, None
    for x0 in starts:
        res = minimize(objective, x0, method="L-BFGS-B")
        if -res.fun > best_val:
            best_val, best_x = -res.fun, res.x
    return best_val, best_x

# ------------------------------------------------------------------
# 3. Classical Fisher information for a local projective measurement
# ------------------------------------------------------------------
def classical_fisher(psi, g, Mfull, phi):
    psi_phi = np.exp(-1j * phi * g) * psi
    dpsi = -1j * g * psi_phi
    a = Mfull @ psi_phi
    da = Mfull @ dpsi
    prob = np.abs(a) ** 2
    dprob = 2.0 * np.real(np.conj(a) * da)
    mask = prob > 1e-12
    return np.sum(dprob[mask] ** 2 / prob[mask])

def measurement_matrix(N, theta):
    c, s = np.cos(theta / 2), np.sin(theta / 2)
    Mt = np.array([[c, s], [-s, c]])
    M = Mt
    for _ in range(N - 1):
        M = np.kron(M, Mt)
    return M

# ------------------------------------------------------------------
# 4. Experiments
# ------------------------------------------------------------------
t0 = time.time()

# (A) Noiseless: optimal pure probe for N = 1..5
Ns = np.arange(1, 6)
qfi_opt_A = []
for N in Ns:
    val, _ = optimize_pure(int(N))
    qfi_opt_A.append(val)
qfi_opt_A = np.array(qfi_opt_A)

print("=== (A) Noiseless case ===")
print(f"{'N':>3} {'product':>9} {'GHZ':>9} {'optimized':>11} {'N^2':>6}")
for N, v in zip(Ns, qfi_opt_A):
    g = generator_diag(int(N))
    print(f"{N:>3} {pure_qfi(product_state(int(N)), g):>9.4f} "
          f"{pure_qfi(ghz_state(int(N)), g):>9.4f} {v:>11.4f} {N**2:>6d}")

# (B) Dephasing: N = 3, optimize the probe for each p (warm start)
N3 = 3
g3 = generator_diag(N3)
p_grid = np.linspace(0.0, 0.45, 10)
qfi_opt_B, qfi_ghz_B, qfi_prod_B, x_store = [], [], [], []
x_warm = None
for p in p_grid:
    K = dephasing_kernel(N3, p)
    val, x_warm = optimize_noisy(N3, p, n_restarts=3, x_warm=x_warm)
    qfi_opt_B.append(val)
    x_store.append(x_warm.copy())
    qfi_ghz_B.append(mixed_qfi(np.outer(ghz_state(N3), ghz_state(N3).conj()) * K, g3))
    qfi_prod_B.append(mixed_qfi(np.outer(product_state(N3), product_state(N3).conj()) * K, g3))
qfi_opt_B, qfi_ghz_B, qfi_prod_B = map(np.array, (qfi_opt_B, qfi_ghz_B, qfi_prod_B))

print("\n=== (B) Local dephasing, N = 3 ===")
print(f"{'p':>6} {'product':>9} {'GHZ':>9} {'optimized':>11}")
for p, a, b, c in zip(p_grid, qfi_prod_B, qfi_ghz_B, qfi_opt_B):
    print(f"{p:>6.3f} {a:>9.4f} {b:>9.4f} {c:>11.4f}")

# (C) Measurement: tilted local measurement for GHZ and product probes
theta_grid = np.linspace(0.05, np.pi - 0.05, 60)
phi_grid = np.linspace(0.011, np.pi / N3 - 0.011, 80)
psi_ghz, psi_prod = ghz_state(N3), product_state(N3)
FC_ghz = np.zeros((len(theta_grid), len(phi_grid)))
FC_prod = np.zeros_like(FC_ghz)
for i, th in enumerate(theta_grid):
    M = measurement_matrix(N3, th)
    for j, ph in enumerate(phi_grid):
        FC_ghz[i, j] = classical_fisher(psi_ghz, g3, M, ph)
        FC_prod[i, j] = classical_fisher(psi_prod, g3, M, ph)
ratio_ghz = FC_ghz.mean(axis=1) / N3**2
ratio_prod = FC_prod.mean(axis=1) / N3

i_mid = len(theta_grid) // 2
print("\n=== (C) Measurement ===")
print(f"GHZ, computational-basis-like tilt (theta ~ 0): mean F_C = {FC_ghz[0].mean():.4f}")
print(f"GHZ, near-equatorial tilt (theta ~ pi/2):       mean F_C = {FC_ghz[i_mid].mean():.4f}  (F_Q = {N3**2})")
print(f"Product, near-equatorial tilt (theta ~ pi/2):   mean F_C = {FC_prod[i_mid].mean():.4f}  (F_Q = {N3})")

# (D) Entanglement advantage over (N, p)
N_list = np.arange(1, 7)
p_list = np.linspace(0.0, 0.45, 31)
ratio_surface = np.zeros((len(N_list), len(p_list)))
max_dev = 0.0
for a, N in enumerate(N_list):
    N = int(N)
    g = generator_diag(N)
    pg, pp = ghz_state(N), product_state(N)
    for b, p in enumerate(p_list):
        K = dephasing_kernel(N, p)
        fg = mixed_qfi(np.outer(pg, pg.conj()) * K, g)
        fp = mixed_qfi(np.outer(pp, pp.conj()) * K, g)
        ratio_surface[a, b] = fg / fp
        max_dev = max(max_dev, abs(fg - N**2 * (1 - 2*p)**(2*N)),
                      abs(fp - N * (1 - 2*p)**2))
p_star = 0.5 * (1 - 3 ** (-1 / 4))
print("\n=== (D) Closed-form check ===")
print(f"max |numerical - analytic| over the (N, p) grid = {max_dev:.2e}")
print(f"N = 3: GHZ stops beating the product state near p = {p_star:.4f}")
print(f"Total computation time: {time.time() - t0:.1f} s")

# ------------------------------------------------------------------
# 5. One figure with all results
# ------------------------------------------------------------------
fig = plt.figure(figsize=(20, 12))
gs = GridSpec(2, 3, figure=fig, hspace=0.32, wspace=0.28)

ax1 = fig.add_subplot(gs[0, 0])
Nf = np.linspace(1, 5, 100)
ax1.plot(Nf, Nf, "b--", label=r"Standard quantum limit $N$")
ax1.plot(Nf, Nf**2, "r-", label=r"Heisenberg limit $N^2$")
ax1.plot(Ns, qfi_opt_A, "ko", ms=9, mfc="none", mew=2, label="Numerically optimized probe")
ax1.set_xlabel("Number of qubits $N$")
ax1.set_ylabel(r"Quantum Fisher information $F_Q$")
ax1.set_title("(a) Noiseless optimal probe")
ax1.legend()
ax1.grid(alpha=0.3)

ax2 = fig.add_subplot(gs[0, 1])
ax2.plot(p_grid, qfi_prod_B, "b--o", label=r"Product state $|+\rangle^{\otimes 3}$")
ax2.plot(p_grid, qfi_ghz_B, "r-s", label="GHZ state")
ax2.plot(p_grid, qfi_opt_B, "k-", lw=2.5, label="Optimized probe")
ax2.set_xlabel("Dephasing probability $p$")
ax2.set_ylabel(r"$F_Q$")
ax2.set_title("(b) Probe optimization under dephasing ($N=3$)")
ax2.legend()
ax2.grid(alpha=0.3)

ax3 = fig.add_subplot(gs[0, 2])
labels = [format(k, "03b") for k in range(2**N3)]
width = 0.27
for off, ip in zip((-width, 0.0, width), (0, 3, 6)):
    probs = np.abs(to_state(x_store[ip])) ** 2
    ax3.bar(np.arange(2**N3) + off, probs, width=width, label=f"$p={p_grid[ip]:.2f}$")
ax3.set_xticks(np.arange(2**N3))
ax3.set_xticklabels([f"|{s}>" for s in labels])
ax3.set_ylabel("Population of the optimal probe")
ax3.set_title("(c) Structure of the optimal input state")
ax3.legend()
ax3.grid(alpha=0.3, axis="y")

ax4 = fig.add_subplot(gs[1, 0], projection="3d")
PH, TH = np.meshgrid(phi_grid, theta_grid)
surf = ax4.plot_surface(PH, TH, FC_ghz, cmap="viridis", edgecolor="none", alpha=0.95)
ax4.set_xlabel(r"Phase $\phi$")
ax4.set_ylabel(r"Measurement tilt $\vartheta$")
ax4.set_zlabel(r"$F_C$", labelpad=6)
ax4.set_title("(d) GHZ probe: classical Fisher information (3D)")
ax4.view_init(elev=28, azim=-60)
fig.colorbar(surf, ax=ax4, shrink=0.55, pad=0.08)

ax5 = fig.add_subplot(gs[1, 1])
ax5.plot(theta_grid, ratio_ghz, "r-", lw=2.5, label=r"GHZ: $\bar F_C / N^2$")
ax5.plot(theta_grid, ratio_prod, "b--", lw=2.5, label=r"Product: $\bar F_C / N$")
ax5.axvline(np.pi / 2, color="gray", ls=":")
ax5.set_xlabel(r"Measurement tilt $\vartheta$ (rad)")
ax5.set_ylabel(r"$\bar F_C / F_Q$ (averaged over $\phi$)")
ax5.set_title("(e) How much of $F_Q$ does the measurement extract?")
ax5.legend()
ax5.grid(alpha=0.3)

ax6 = fig.add_subplot(gs[1, 2], projection="3d")
PP, NN = np.meshgrid(p_list, N_list)
surf2 = ax6.plot_surface(PP, NN, ratio_surface, cmap="plasma", edgecolor="none", alpha=0.95)
ax6.plot_surface(PP, NN, np.ones_like(ratio_surface), color="gray", alpha=0.25)
ax6.set_xlabel("Dephasing $p$")
ax6.set_ylabel("Qubits $N$")
ax6.set_zlabel("GHZ / product", labelpad=6)
ax6.set_title("(f) Entanglement advantage (gray plane = 1)")
ax6.view_init(elev=26, azim=-125)
fig.colorbar(surf2, ax=ax6, shrink=0.55, pad=0.08)

fig.suptitle("Optimal probe design in quantum metrology", fontsize=18, y=0.98)
plt.show()

3. Code Walkthrough

Building blocks

Since $G=\tfrac12\sum_i\sigma_z^{(i)}$ is diagonal in the computational basis, it is stored as a vector g. The basis state with $n_1$ ones has eigenvalue $(N-2n_1)/2$, and popcount_table counts those ones.

dephasing_kernel is the main speed-up. Instead of applying $3N$ Kraus operators one at a time, it builds the matrix $K_{ij}=(1-2p)^{d_H(i,j)}$ once. The Hamming distance comes from the bitwise XOR of the indices, idx[:,None] ^ idx[None,:]. The noisy state is then a single element-wise product, rho * K.

pure_qfi implements $4,\mathrm{Var}(G)$ using only the populations $|\psi_k|^2$. mixed_qfi implements the eigendecomposition formula. The mask S > 1e-12 removes the $0/0$ terms that appear when two eigenvalues vanish, which would otherwise produce NaN.

Parametrizing and optimizing the probe

to_state maps a real vector of length $2\cdot 2^N$ to a normalized complex state. This turns the constrained search over unit vectors into an unconstrained problem that L-BFGS-B can handle.

optimize_pure minimizes $-F_Q$ from several random starting points and keeps the best result. It uses the cheap pure-state formula, which is why $N=5$ (a 32-dimensional state) finishes almost instantly.

optimize_noisy does the same for the noisy QFI, but with two additions that speed things up. It evaluates $F_Q$ on rho * K, and it accepts a warm start x_warm. As $p$ increases in small steps, the previous optimum is an excellent initial guess, so only three random restarts per $p$ are needed.

Measurement model

measurement_matrix builds $M^{\dagger\otimes N}$, where $M=R_y(\vartheta)$ has the eigenvectors of $\sigma_\vartheta$ as columns. Measuring $\sigma_\vartheta$ on every qubit means the outcome probabilities are $p_k=|(M^{\dagger\otimes N}\psi_\varphi)_k|^2$.

classical_fisher computes the derivative analytically,

$$
\partial_\varphi p_k = 2,\mathrm{Re}!\left[a_k^*,\partial_\varphi a_k\right],\qquad \partial_\varphi\psi_\varphi=-iG\psi_\varphi ,
$$

so no finite-difference error enters. Outcomes with $p_k\le10^{-12}$ are masked out to avoid division by zero.

Experiments (A)–(D)

(A) For $N=1,\dots,5$ the code optimizes over all pure states and compares the result with $N$ and $N^2$.
(B) For $N=3$ and ten values of $p$ it compares the product state, GHZ and the numerically optimized probe.
(C) It sweeps the measurement tilt $\vartheta$ and the phase $\varphi$ on a $60\times80$ grid. The phase range $0<\varphi<\pi/3$ covers one period of the GHZ signal $\cos 3\varphi$. The grids are offset slightly so they avoid the exact points where outcome probabilities vanish.
(D) It evaluates the ratio $F_Q^{\rm GHZ}/F_Q^{\rm prod}$ for $N=1,\dots,6$ and 31 values of $p$. It also reports the largest deviation from the closed-form expressions, which serves as a built-in correctness check of mixed_qfi.

The figure

All six panels are drawn in a single GridSpec figure. Panels (d) and (f) are 3D surfaces.

4. Results

=== (A) Noiseless case ===
  N   product       GHZ   optimized    N^2
  1    1.0000    1.0000      1.0000      1
  2    2.0000    4.0000      4.0000      4
  3    3.0000    9.0000      9.0000      9
  4    4.0000   16.0000     16.0000     16
  5    5.0000   25.0000     25.0000     25

=== (B) Local dephasing, N = 3 ===
     p   product       GHZ   optimized
 0.000    3.0000    9.0000      9.0000
 0.050    2.4300    4.7830      4.7830
 0.100    1.9200    2.3593      2.6425
 0.150    1.4700    1.0588      1.7240
 0.200    1.0800    0.4199      1.1678
 0.250    0.7500    0.1406      0.7770
 0.300    0.4800    0.0369      0.4867
 0.350    0.2700    0.0066      0.2711
 0.400    0.1200    0.0006      0.1201
 0.450    0.0300    0.0000      0.0300

=== (C) Measurement ===
GHZ, computational-basis-like tilt (theta ~ 0): mean F_C = 0.0000
GHZ, near-equatorial tilt (theta ~ pi/2):       mean F_C = 8.7309  (F_Q = 9)
Product, near-equatorial tilt (theta ~ pi/2):   mean F_C = 2.8973  (F_Q = 3)

=== (D) Closed-form check ===
max |numerical - analytic| over the (N, p) grid = 2.13e-14
N = 3: GHZ stops beating the product state near p = 0.1201
Total computation time: 2.1 s

5. Discussion of the Results

(a) Noiseless optimum: the Heisenberg limit emerges from the optimizer

The black circles are the best QFI values found by the optimizer for $N=1,\dots,5$. They sit on the red curve $N^2$ and well above the blue line $N$. The optimizer started from random states and was never told about GHZ, yet it reproduced the bound $F_Q\le(g_{\max}-g_{\min})^2=N^2$. For $N=3$ the values are $3$ for the product state and $9$ for both GHZ and the optimized probe.

(b) Under dephasing, GHZ is not always best

At $p=0$ the optimized probe equals GHZ ($F_Q=9$). GHZ coherence decays as $(1-2p)^{N}$, so its QFI collapses like $(1-2p)^{2N}$ while the product state’s QFI decays only like $(1-2p)^2$. The curves cross near $p\approx0.12$. In the intermediate regime ($p\approx0.1$–$0.3$) the optimized probe is strictly better than both reference states. At $p=0.10$, for example, the optimized value is about $2.64$, against roughly $2.36$ for GHZ and $1.92$ for the product state. For very large $p$ the optimum merges with the product state. No fixed textbook state is optimal across the whole range, which is the point of probe optimization.

(c) Structure of the optimal state

At $p=0$ the optimizer produces a GHZ state, with weight $0.5$ on $|000\rangle$ and $|111\rangle$ and nothing elsewhere. As noise grows, population leaks into the other six basis states. The weight on the two extreme states drops while the rest rises toward a uniform distribution, which is the signature of a product state $|+\rangle^{\otimes3}$ (every population equal to $1/8$). The optimal probe therefore interpolates between GHZ and a product state, shifting toward the more robust end as noise increases.

(d) 3D view: the measurement matters as much as the state

This surface shows the classical Fisher information $F_C(\varphi,\vartheta)$ of the GHZ probe. Near $\vartheta=\pi/2$, the $x$-basis, there is a high plateau close to $F_Q=9$. As $\vartheta\to0$ or $\pi$ the surface falls to zero. In the computational basis the outcome statistics of a GHZ state do not depend on $\varphi$ at all, so a perfect probe combined with the wrong measurement carries no information. The dips at the edges of the plateau, near $\varphi\to0$ and $\varphi\to\pi/3$, appear because some outcome probabilities approach zero there.

(e) How much of the QFI the measurement extracts

Averaged over $\varphi$, both the GHZ and product curves peak near $\vartheta=\pi/2$ at roughly $0.9$–$1.0$ of their respective $F_Q$. The GHZ curve is narrower, so entangled probes demand more precise measurement alignment. Entanglement raises the achievable sensitivity but makes it more fragile with respect to measurement errors.

(f) 3D view: where entanglement pays off

The surface is the ratio $F_Q^{\rm GHZ}/F_Q^{\rm prod}=N(1-2p)^{2N-2}$, and the gray plane marks ratio $=1$. In the noiseless limit the ratio equals $N$, so it grows linearly with the number of qubits. For any $p>0$ the exponent $2N-2$ makes larger GHZ states lose their advantage faster. The surface drops below the plane in a region that expands as $N$ increases, so the best probe size depends on the noise level. The console check confirms that the numerical QFI agrees with the closed forms to machine precision.

6. Takeaways

The QFI bounds the best possible precision, so maximizing it over input states is the first step of probe design. In the noiseless case this gives the Heisenberg-limited GHZ state.
The measurement must be matched to the probe. A GHZ state read out in the computational basis yields $F_C=0$, while an $x$-basis readout approaches $F_Q$.
With noise, the optimal state is neither GHZ nor a product state. It is a numerically found compromise, and a short gradient-based optimization over $2^N$ amplitudes is enough to locate it.
The same pipeline carries over to other generators, noise channels and qubit numbers. Only generator_diag and dephasing_kernel need to change.

Variational Quantum Algorithms in Practice

October 19, 2026

Solving a Hydrogen-Molecule VQE and a MaxCut QAOA from Scratch in Python

Variational quantum algorithms (VQE and QAOA) split the work between two machines. A parameterized quantum circuit prepares a state $|\psi(\boldsymbol{\theta})\rangle$ and measures a cost function, and a classical optimizer updates $\boldsymbol{\theta}$ to minimize it. This hybrid loop is the standard strategy for today’s noisy, small-scale quantum devices.

In this article we solve two concrete problems with nothing but NumPy and SciPy:

VQE finds the ground-state energy of a two-qubit hydrogen-molecule Hamiltonian.
QAOA solves a MaxCut problem on a 6-node graph, with circuit depths $p = 1, 2, 3, 4$.

Everything is simulated with an exact state-vector, so no quantum SDK is required.

1. Theory

1.1 VQE (Variational Quantum Eigensolver)

Given a Hamiltonian $H$, VQE minimizes the expectation value

$$
E(\boldsymbol{\theta}) = \langle \psi(\boldsymbol{\theta}) | H | \psi(\boldsymbol{\theta}) \rangle ;\ge; E_0 ,
$$

where $E_0$ is the true ground-state energy. The inequality is the variational principle: any trial state gives an upper bound on $E_0$, so minimizing over $\boldsymbol{\theta}$ pushes the energy down toward the ground state.

The Hamiltonian is written as a sum of Pauli strings,

$$
H = \sum_k c_k P_k, \qquad P_k \in {I, X, Y, Z}^{\otimes n},
$$

and on hardware each $\langle P_k\rangle$ is estimated separately from measurements. Our simulator computes the full expectation exactly.

1.2 QAOA (Quantum Approximate Optimization Algorithm)

For MaxCut on a graph $G=(V,E)$, the cost operator is

$$
C = \sum_{(u,v)\in E} \frac{1 - Z_u Z_v}{2},
$$

whose eigenvalue on a computational basis state $|z\rangle$ equals the number of cut edges. A depth-$p$ QAOA state is

$$
|\boldsymbol{\gamma},\boldsymbol{\beta}\rangle = \prod_{l=1}^{p} e^{-i\beta_l B}, e^{-i\gamma_l C}, |+\rangle^{\otimes n}, \qquad B = \sum_{j} X_j .
$$

The classical optimizer maximizes

$$
F_p(\boldsymbol{\gamma},\boldsymbol{\beta}) = \langle \boldsymbol{\gamma},\boldsymbol{\beta} | C | \boldsymbol{\gamma},\boldsymbol{\beta} \rangle ,
$$

and the quality of the result is measured by the approximation ratio

$$
r_p = \frac{F_p^{*}}{C_{\max}} \in [0,1].
$$

As $p\to\infty$, QAOA is guaranteed to approach $r=1$, so a deeper circuit should never do worse.

2. Problem Setup

Problem A: the H$_2$ molecule (VQE)

At a bond length of $0.735,\text{Å}$, after parity mapping and two-qubit reduction, the electronic Hamiltonian of H$_2$ becomes

$$
H = g_0, II + g_1, IZ + g_2, ZI + g_3, ZZ + g_4, XX,
$$

with

$$
g_0=-1.0524,\quad g_1=0.3979,\quad g_2=-0.3979,\quad g_3=-0.0113,\quad g_4=0.1809 \ \text{(Hartree)}.
$$

Our ansatz uses four parameters:

$$
|\psi(\boldsymbol{\theta})\rangle = \big(R_y(\theta_2)\otimes R_y(\theta_3)\big);\mathrm{CNOT}_{0\to1};\big(R_y(\theta_0)\otimes R_y(\theta_1)\big),|00\rangle ,
$$

where

$$
R_y(\theta)=\begin{pmatrix}\cos\frac{\theta}{2} & -\sin\frac{\theta}{2}\[2pt] \sin\frac{\theta}{2} & \cos\frac{\theta}{2}\end{pmatrix}.
$$

Because $H$ is a real matrix, its ground state is real. Every real two-qubit state can be written as a local real rotation applied to $\cos a,|00\rangle+\sin a,|11\rangle$, so this ansatz is expressive enough to reach the exact ground state.

Problem B: MaxCut on a 6-node graph (QAOA)

The graph has 6 nodes and 8 edges: a ring $0!-!1!-!2!-!3!-!4!-!5!-!0$ plus the two chords $(0,3)$ and $(1,4)$. We find the best parameters for $p=1,2,3,4$ and compare the result with the brute-force optimum.

3. Full Source Code

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

plt.style.use('dark_background')
rng = np.random.default_rng(7)
t_start = time.time()

# =====================================================================
# Part A : VQE for the two-qubit H2 Hamiltonian
# =====================================================================
I2 = np.eye(2, dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)


def kron(*mats):
    out = np.array([[1.0 + 0j]])
    for m in mats:
        out = np.kron(out, m)
    return out


def ry(t):
    c, s = np.cos(t / 2.0), np.sin(t / 2.0)
    return np.array([[c, -s], [s, c]], dtype=complex)


CNOT = np.array([[1, 0, 0, 0],
                 [0, 1, 0, 0],
                 [0, 0, 0, 1],
                 [0, 0, 1, 0]], dtype=complex)

g0, g1, g2, g3, g4 = (-1.052373245772859, 0.39793742484318045,
                      -0.39793742484318045, -0.01128010425623538,
                      0.18093119978423156)

H = (g0 * kron(I2, I2) + g1 * kron(I2, Z) + g2 * kron(Z, I2)
     + g3 * kron(Z, Z) + g4 * kron(X, X))

E_exact = float(np.linalg.eigvalsh(H)[0])
psi_init = np.array([1, 0, 0, 0], dtype=complex)


def vqe_state(th):
    u1 = kron(ry(th[0]), ry(th[1]))
    u2 = kron(ry(th[2]), ry(th[3]))
    return u2 @ (CNOT @ (u1 @ psi_init))


def vqe_energy(th):
    psi = vqe_state(th)
    return float(np.real(np.vdot(psi, H @ psi)))


histories, results = [], []
for run in range(5):
    hist = []

    def objective(th, hist=hist):
        e = vqe_energy(th)
        hist.append(e)
        return e

    x0 = rng.uniform(-np.pi, np.pi, 4)
    res = minimize(objective, x0, method='COBYLA', tol=1e-10,
                   options={'maxiter': 400, 'rhobeg': 0.5})
    histories.append(np.array(hist))
    results.append(res)

best_run = int(np.argmin([r.fun for r in results]))
th_opt = results[best_run].x
E_vqe = results[best_run].fun

print("=" * 60)
print("Part A : VQE (H2, bond length 0.735 A)")
print("=" * 60)
print(f"Exact ground energy (diagonalization) : {E_exact:.10f} Ha")
print(f"VQE energy (best of 5 restarts)       : {E_vqe:.10f} Ha")
print(f"Absolute error                        : {abs(E_vqe - E_exact):.3e} Ha")
print("Optimal parameters (rad)              :", np.round(th_opt, 5))

# 3D landscape: vary theta_0, theta_1 with theta_2, theta_3 fixed at optimum
grid_n = 61
t0_grid = np.linspace(-np.pi, np.pi, grid_n)
t1_grid = np.linspace(-np.pi, np.pi, grid_n)
T0, T1 = np.meshgrid(t0_grid, t1_grid)
E_land = np.zeros_like(T0)
for i in range(grid_n):
    for j in range(grid_n):
        E_land[i, j] = vqe_energy(np.array([T0[i, j], T1[i, j],
                                            th_opt[2], th_opt[3]]))

w0 = float(np.arctan2(np.sin(th_opt[0]), np.cos(th_opt[0])))
w1 = float(np.arctan2(np.sin(th_opt[1]), np.cos(th_opt[1])))
E_mark = vqe_energy(np.array([w0, w1, th_opt[2], th_opt[3]]))

# =====================================================================
# Part B : QAOA for MaxCut on a 6-node graph
# =====================================================================
n = 6
dim = 2 ** n
edges = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 0), (0, 3), (1, 4)]

idx_all = np.arange(dim)
cut = np.zeros(dim)
for (u, v) in edges:
    cut += (((idx_all >> u) & 1) != ((idx_all >> v) & 1)).astype(float)
max_cut = float(cut.max())
n_opt = int(np.sum(cut == max_cut))


def apply_rx_all(psi, beta):
    c, s = np.cos(beta), np.sin(beta)
    U = np.array([[c, -1j * s], [-1j * s, c]], dtype=complex)
    t = psi.reshape([2] * n)
    for ax in range(n):
        t = np.moveaxis(np.tensordot(U, t, axes=([1], [ax])), 0, ax)
    return t.reshape(-1)


def qaoa_state(params, p):
    gammas, betas = params[:p], params[p:]
    psi = np.full(dim, 1.0 / np.sqrt(dim), dtype=complex)
    for g, b in zip(gammas, betas):
        psi = np.exp(-1j * g * cut) * psi
        psi = apply_rx_all(psi, b)
    return psi


def qaoa_expect(params, p):
    psi = qaoa_state(params, p)
    return float(np.real(np.sum(np.abs(psi) ** 2 * cut)))


def neg_expect(params, p):
    return -qaoa_expect(params, p)


P_MAX = 4
N_RESTART = 10
best_params, best_val = {}, {}
for p in range(1, P_MAX + 1):
    best = None
    for r in range(N_RESTART):
        x0 = np.concatenate([rng.uniform(0, np.pi, p),
                             rng.uniform(0, np.pi / 2, p)])
        res = minimize(neg_expect, x0, args=(p,), method='L-BFGS-B')
        if best is None or res.fun < best.fun:
            best = res
    best_params[p] = best.x
    best_val[p] = -best.fun

ps = np.arange(1, P_MAX + 1)
ratios = np.array([best_val[p] / max_cut for p in ps])

print()
print("=" * 60)
print("Part B : QAOA (MaxCut, 6 nodes, 8 edges)")
print("=" * 60)
print(f"Brute-force maximum cut : {max_cut:.0f}  ({n_opt} optimal bitstrings)")
print(f"{'p':>3} | {'<C>':>10} | {'ratio':>8}")
print("-" * 30)
for p in ps:
    print(f"{p:>3} | {best_val[p]:>10.5f} | {best_val[p] / max_cut:>8.5f}")

psi_final = qaoa_state(best_params[P_MAX], P_MAX)
probs = np.abs(psi_final) ** 2
order = np.argsort(probs)[::-1][:4]
print()
print(f"Top bitstrings at p = {P_MAX} (qubit 0 is the rightmost bit):")
for k in order:
    print(f"  |{k:0{n}b}>  prob = {probs[k]:.4f}   cut = {cut[k]:.0f}")
print(f"Probability of sampling an optimal cut : {probs[cut == max_cut].sum():.4f}")

# 3D landscape for p = 1
gam_grid = np.linspace(0, 2 * np.pi, 61)
bet_grid = np.linspace(0, np.pi, 61)
GG, BB = np.meshgrid(gam_grid, bet_grid)
F_land = np.zeros_like(GG)
for i in range(GG.shape[0]):
    for j in range(GG.shape[1]):
        F_land[i, j] = qaoa_expect(np.array([GG[i, j], BB[i, j]]), 1)

g1_opt = float(best_params[1][0] % (2 * np.pi))
b1_opt = float(best_params[1][1] % np.pi)
F_mark = qaoa_expect(np.array([g1_opt, b1_opt]), 1)

# =====================================================================
# Combined figure
# =====================================================================
fig = plt.figure(figsize=(21, 12))
fig.suptitle("Variational Quantum Algorithms: VQE (H$_2$) and QAOA (MaxCut)",
             fontsize=20, y=0.995)

# (1) VQE convergence
ax1 = fig.add_subplot(2, 3, 1)
for k, h in enumerate(histories):
    err = np.maximum(np.minimum.accumulate(h) - E_exact, 1e-12)
    if k == best_run:
        ax1.semilogy(err, color='cyan', lw=2.5, label='best run')
    else:
        ax1.semilogy(err, color='gray', lw=1.0, alpha=0.7)
ax1.axhline(1.6e-3, color='orange', ls='--', lw=1.5,
            label='chemical accuracy (1.6 mHa)')
ax1.set_xlabel("Cost-function evaluations")
ax1.set_ylabel(r"$E(\theta) - E_{\mathrm{exact}}$  [Ha]")
ax1.set_title("(1) VQE convergence (5 restarts)")
ax1.legend(loc='upper right')
ax1.grid(alpha=0.25)

# (2) VQE 3D landscape
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
ax2.plot_surface(T0, T1, E_land, cmap='viridis', alpha=0.92,
                 linewidth=0, antialiased=True)
ax2.scatter([w0], [w1], [E_mark], color='red', s=90, depthshade=False,
            label='VQE optimum')
ax2.set_xlabel(r"$\theta_0$")
ax2.set_ylabel(r"$\theta_1$")
ax2.set_zlabel("Energy [Ha]")
ax2.set_title(r"(2) VQE energy landscape ($\theta_2,\theta_3$ fixed)")
ax2.legend(loc='upper left')
ax2.view_init(elev=28, azim=-60)

# (3) Graph with best cut
ax3 = fig.add_subplot(2, 3, 3)
best_idx = int(np.argmax(probs))
bits = [(best_idx >> k) & 1 for k in range(n)]
ang = 2 * np.pi * np.arange(n) / n
pos = np.column_stack([np.cos(ang), np.sin(ang)])
for (u, v) in edges:
    is_cut = bits[u] != bits[v]
    ax3.plot([pos[u, 0], pos[v, 0]], [pos[u, 1], pos[v, 1]],
             color='lime' if is_cut else 'gray',
             ls='-' if is_cut else '--', lw=3 if is_cut else 1.5, zorder=1)
for k in range(n):
    ax3.scatter(pos[k, 0], pos[k, 1], s=900,
                color='deepskyblue' if bits[k] == 0 else 'orange',
                edgecolor='white', zorder=2)
    ax3.text(pos[k, 0], pos[k, 1], str(k), ha='center', va='center',
             color='black', fontsize=14, fontweight='bold', zorder=3)
ax3.set_aspect('equal')
ax3.axis('off')
ax3.set_title(f"(3) Most probable cut at p={P_MAX}: {cut[best_idx]:.0f} of "
              f"{max_cut:.0f} edges")

# (4) QAOA p=1 3D landscape
ax4 = fig.add_subplot(2, 3, 4, projection='3d')
ax4.plot_surface(GG, BB, F_land, cmap='magma', alpha=0.92,
                 linewidth=0, antialiased=True)
ax4.scatter([g1_opt], [b1_opt], [F_mark], color='cyan', s=90,
            depthshade=False, label='p=1 optimum')
ax4.set_xlabel(r"$\gamma$")
ax4.set_ylabel(r"$\beta$")
ax4.set_zlabel(r"$\langle C \rangle$")
ax4.set_title(r"(4) QAOA landscape at $p=1$")
ax4.legend(loc='upper left')
ax4.view_init(elev=30, azim=-55)

# (5) Approximation ratio vs depth
ax5 = fig.add_subplot(2, 3, 5)
ax5.plot(ps, ratios, 'o-', color='cyan', lw=2.5, ms=10)
ax5.axhline(1.0, color='orange', ls='--', lw=1.5, label='optimum')
for p, r in zip(ps, ratios):
    ax5.text(p, r - 0.012, f"{r:.3f}", ha='center', va='top', fontsize=11)
ax5.set_xticks(ps)
ax5.set_ylim(min(ratios) - 0.06, 1.03)
ax5.set_xlabel("QAOA depth $p$")
ax5.set_ylabel("Approximation ratio")
ax5.set_title("(5) Approximation ratio vs depth")
ax5.legend(loc='lower right')
ax5.grid(alpha=0.25)

# (6) Output distribution at p = P_MAX
ax6 = fig.add_subplot(2, 3, 6)
colors = plt.cm.plasma(cut / max_cut)
ax6.bar(np.arange(dim), probs, color=colors, width=1.0)
sm = plt.cm.ScalarMappable(cmap='plasma', norm=plt.Normalize(0, max_cut))
sm.set_array([])
fig.colorbar(sm, ax=ax6, label="Cut value")
ax6.set_xlabel("Bitstring index $z$")
ax6.set_ylabel("Probability")
ax6.set_title(f"(6) Output distribution at $p={P_MAX}$")
ax6.grid(alpha=0.25)

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

print()
print(f"Total elapsed time: {time.time() - t_start:.1f} s")

4. Code Walkthrough

4.1 Part A: VQE

Hamiltonian construction. The helper kron chains Kronecker products, so kron(I2, Z) is $I\otimes Z$. The first factor is qubit 0 (the most significant bit). The five Pauli terms are assembled into a $4\times4$ matrix H. We also diagonalize H exactly with eigvalsh, giving E_exact as a reference to judge the VQE result.

Ansatz. vqe_state applies the four-parameter circuit directly to the vector $|00\rangle$: a layer of $R_y$ gates, a CNOT, and a second layer of $R_y$ gates. vqe_energy evaluates $\langle\psi|H|\psi\rangle$ with np.vdot, which conjugates its first argument. The result is real because $H$ is Hermitian.

Optimization loop. We run COBYLA from five random starting points, because variational landscapes can have local minima. COBYLA is gradient-free, which mirrors how shot-noise-limited hardware experiments are often run. Each run stores its full history of energy evaluations, so we can draw convergence curves later. The run with the lowest final energy is kept.

Landscape scan. To visualize the cost function, we vary $\theta_0$ and $\theta_1$ over $[-\pi,\pi]$ while holding $\theta_2,\theta_3$ at their optimal values. The energy has period $2\pi$ in every parameter, so the optimum is wrapped back into the plotted range with arctan2(sin, cos) before it is marked on the surface.

4.2 Part B: QAOA

Cost vector. The cost operator $C$ is diagonal in the computational basis. We therefore precompute cut, a vector of length $2^6=64$ holding the number of cut edges for every bitstring, using bit masks. The phase-separation step $e^{-i\gamma C}$ then becomes one element-wise multiplication, np.exp(-1j*g*cut) * psi.

Mixer. The mixer $e^{-i\beta B}=\bigotimes_j e^{-i\beta X_j}$ acts on each qubit independently. We reshape the 64-element vector into a tensor of shape $(2,2,2,2,2,2)$ and contract a $2\times2$ matrix into one axis at a time with tensordot. Building the full $64\times64$ mixer matrix is never necessary.

Optimization. For each depth $p$ we maximize $F_p$ by minimizing its negative with L-BFGS-B from 10 random starting points. Random restarts matter here because the QAOA landscape has many equivalent and inequivalent local optima.

Post-processing. The code evaluates the optimized state at $p=4$ and extracts the probability distribution $|\psi_z|^2$. It also computes the total probability of sampling a maximum cut, which is the quantity that matters when you actually measure the circuit.

4.3 Why this code is fast

The slow way to write this simulation is to build every gate as a full $2^n\times2^n$ matrix and multiply, which costs $O(4^n)$ per gate. Three choices avoid that:

the diagonal cost is stored as a vector, so the phase step costs $O(2^n)$;
the mixer is applied qubit by qubit through tensordot, costing $O(n,2^n)$ instead of $O(4^n)$;
the Pauli matrices for VQE are tiny ($4\times4$), so exact evaluation is essentially free.

With these choices both landscape scans (3,721 evaluations each) and all optimizations finish quickly on a standard Colab CPU.

5. Execution Results

============================================================
Part A : VQE (H2, bond length 0.735 A)
============================================================
Exact ground energy (diagonalization) : -1.8572750302 Ha
VQE energy (best of 5 restarts)       : -1.8572750302 Ha
Absolute error                        : 1.332e-15 Ha
Optimal parameters (rad)              : [-2.30741  1.87472 -3.95238  1.77846]

============================================================
Part B : QAOA (MaxCut, 6 nodes, 8 edges)
============================================================
Brute-force maximum cut : 8  (2 optimal bitstrings)
  p |         |    ratio
------------------------------
  1 |    5.63245 |  0.70406
  2 |    7.02977 |  0.87872
  3 |    7.04608 |  0.88076
  4 |    7.97840 |  0.99730

Top bitstrings at p = 4 (qubit 0 is the rightmost bit):
  |101010>  prob = 0.4968   cut = 8
  |010101>  prob = 0.4968   cut = 8
  |010001>  prob = 0.0007   cut = 6
  |101110>  prob = 0.0007   cut = 6
Probability of sampling an optimal cut : 0.9935


Total elapsed time: 17.3 s

6. Interpreting the Figure

(1) VQE convergence. Each curve shows the running minimum of the energy error $E(\boldsymbol{\theta})-E_{\text{exact}}$ on a logarithmic axis. By the variational principle this error can never be negative, so the curves decrease monotonically toward zero. The orange dashed line marks chemical accuracy, $1.6\times10^{-3}$ Ha, the threshold at which a quantum-chemistry result becomes practically useful. A run that crosses this line and keeps falling has found the correct basin. Runs that stall higher have been trapped in a local minimum, which is exactly why several restarts are used.

(2) VQE energy landscape (3D). This surface is the energy as a function of $\theta_0$ and $\theta_1$ with the other two angles frozen. Its smooth, periodic shape reflects the trigonometric dependence of the energy on rotation angles. The red marker sits at a minimum of the surface. The picture shows why gradient-free and gradient-based optimizers work well on small, shallow ansätze: the valleys are broad and well-defined. It also hints at the difficulty of larger circuits, where many such valleys coexist.

(3) Best cut found by QAOA. The graph is drawn with its most probable bitstring from the $p=4$ state. Node colors (blue or orange) indicate which side of the partition each node belongs to. Green solid edges are cut, and gray dashed edges are not. A good partition alternates colors around the ring and uses the two chords effectively. Because the cost is symmetric under flipping all bits, the complementary coloring is exactly as probable.

(4) QAOA landscape at $p=1$ (3D). The surface of $\langle C\rangle$ over $(\gamma,\beta)$ shows a periodic ridge-and-valley structure. The cyan marker is the best $p=1$ point found by the optimizer, and it sits on one of the highest ridges. The several equivalent peaks come from the symmetries of the problem, namely the periodicity of the mixer in $\beta$ and the integer-valued cost in $\gamma$. This explains why random restarts consistently find the same maximum value, even from different starting angles.

(5) Approximation ratio versus depth. The ratio $r_p$ should rise as $p$ increases and approach the dashed line at 1. Even $p=1$ already gives a non-trivial guarantee, and each additional layer adds two variational parameters and therefore strictly more expressive power. The gain per layer typically shrinks, which is the usual trade-off between circuit depth, which costs noise on real hardware, and solution quality.

(6) Output distribution. Each bar is the probability of measuring one of the 64 bitstrings, colored by its cut value. After optimization at $p=4$, the probability mass concentrates on bright, high-cut bitstrings, while low-cut strings are suppressed. Because measurements sample from this distribution, the key practical quantity is the total probability on optimal cuts: repeating the circuit a handful of times is usually enough to see a maximum cut at least once.

7. Conclusion

The two examples show the same pattern. A short parameterized circuit prepares a state, a cost function is evaluated, and a classical optimizer improves the parameters.

VQE turned the variational principle into a practical minimization and reproduced the exact ground-state energy of H$_2$ with a four-parameter circuit.
QAOA encoded a combinatorial problem into a phase operator and a mixer, and its solution quality improved with circuit depth $p$.

The same structure scales to larger molecules and graphs. The hard parts there are noise, barren plateaus, and the cost of measurement, but the hybrid quantum-classical loop is unchanged.

Steering Population Growth:

Introduction

Problem Setting

Production and food balance

Budget constraint

Demographics

Objective

Parameters

Solution Strategy

Full Source Code

Detailed Code Walkthrough

Parameters and helper

The simulator

Landscape analysis

Full optimization

Local polish

Reporting and plotting

Execution Results

Result Image

Console Output

Reading the Graphs

Panel 1: 3D landscape of constant policies

Panel 2: Population trajectories

Panel 3: Optimized controls

Panel 4: Food stock and satisfaction

Panel 5: 3D state-space trajectories

Panel 6: Objective value by policy

Discussion

Conclusion

Investing in Health vs. Machines

1. The Model

Production

Mortality

Population dynamics

The optimization problem

Example setup

2. The Code

3. Code Walkthrough

4. Results

Console output

Graph output

5. Interpreting the Results

6. Takeaways

Eat Now or Store for Later?

1. The Problem

2. Mathematical Model

3. Concrete Example

4. Solution Method and Speed-Up

5. Full Source Code

6. Code Walkthrough

7. Results

Console Output

Result Image

8. How to Read the Results

9. Takeaways

When Should You Build the Farm?

1. Problem Setting

State transition

Objective

2. Source Code

3. Code Walkthrough

4. Visualization of the Results

How to read the six panels

5. Console Output

6. Discussion

How Much Ore Should You Dig?

1. The Problem

2. Mathematical Model

3. Solution Strategy

4. Source Code

5. Code Walkthrough

5.1 Parameters and the simulate function

5.2 Exhaustive search

5.3 Continuous refinement

5.4 Baseline policies

5.5 Console report

6. Console Output

How to Read the Console Output

7. Graphs

How to Read the Graphs

5.1 Parameters and the `simulate` function

4.2 Building the linear program: `solve_lp`

4.3 Forward simulation: `simulate`