Optimal Network Segmentation for Zero Trust

Balancing Security, Latency, and Operational Load

Zero Trust and micro-segmentation promise a network in which a single compromised host cannot reach everything else. Yet every additional boundary has a price. Traffic that crosses a boundary passes through a policy enforcement point and picks up latency. Every boundary also adds firewall rules that someone has to write, review, and audit.

Push segmentation to the extreme, one segment per asset, and the attack surface for lateral movement nearly disappears. The rule base then explodes, and east-west traffic pays a latency toll on almost every request. Leave the network flat and operations are trivial, but one phishing click can reach the crown jewels.

In this article we turn that tension into an optimization problem and solve it in Python for a concrete 30-asset network.

1. Problem Setting

We have $N = 30$ assets in six tiers: Web (5), App (8), DB (4), Admin (3), IoT (6), and Corp (4). Each asset $i$ has:

  • a compromise probability $p_i$ (how exposed it is),
  • an asset value $v_i$ (how much damage its loss causes).

Assets communicate according to a symmetric traffic matrix $T_{ij}$ in kreq/s. A link exists when $T_{ij} > 0$, and we write $L_{ij} = \mathbb{1}[T_{ij} > 0]$.

A design is an assignment $a_i \in {0, \dots, K_{\max}-1}$ of every asset to a segment. Define the cross-segment indicator

$$
x_{ij} = \mathbb{1}[a_i \neq a_j].
$$

Risk

Inside a segment the network is flat, so a compromised asset can reach every other member. Across segments, only the allowed flows leak, with a small residual factor $\varepsilon$:

$$
\mathrm{Risk}(a) = \sum_{s} \Big(\sum_{i \in s} p_i\Big)\Big(\sum_{j \in s} v_j\Big) - \sum_i p_i v_i ;+; \varepsilon \sum_{i,j} p_i v_j , L_{ij}, x_{ij}.
$$

Latency

Each cross-segment hop passes an enforcement point that adds $\delta$ ms. The traffic-weighted average added latency is

$$
\mathrm{Lat}(a) = \delta \cdot \frac{\sum_{i<j} T_{ij}, x_{ij}}{\sum_{i<j} T_{ij}}.
$$

Operational load

Operations cost grows with the number of segments to manage, $K_{\mathrm{used}}$, and with the number of firewall rules, one per communicating pair that crosses a boundary:

$$
\mathrm{Ops}(a) = c_{\mathrm{seg}}, K_{\mathrm{used}} + c_{\mathrm{rule}} \sum_{i<j} L_{ij}, x_{ij}.
$$

Normalized objective

To make the weights meaningful, each term is divided by its natural scale: the flat-network risk $R_0$, the maximum latency $\delta$, and the cost $\mathrm{Ops}_{\mathrm{micro}}$ of full micro-segmentation. We minimize

$$
J(a) = w_r \frac{\mathrm{Risk}(a)}{R_0} + w_l \frac{\mathrm{Lat}(a)}{\delta} + w_o \frac{\mathrm{Ops}(a)}{\mathrm{Ops}_{\mathrm{micro}}}, \qquad w_r + w_l + w_o = 1.
$$

The parameters are $\delta = 1.8$ ms, $\varepsilon = 0.08$, $c_{\mathrm{seg}} = 6$, $c_{\mathrm{rule}} = 1$, and the main scenario uses $(w_r, w_l, w_o) = (0.5, 0.3, 0.2)$.

This is a combinatorial problem: the number of ways to partition 30 assets is astronomically large, so brute force is out. We use simulated annealing, with one important speed-up described below.

2. Full Source Code

The whole experiment is a single script.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
import math
import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib.gridspec import GridSpec

plt.style.use("dark_background")

# ------------------------------------------------------------
# 1. Problem definition: 30 assets in 6 tiers
# ------------------------------------------------------------
rng = np.random.default_rng(2026)
TIERS = ["Web", "App", "DB", "Admin", "IoT", "Corp"]
COUNTS = np.array([5, 8, 4, 3, 6, 4])
N = int(COUNTS.sum())
tier_of = np.repeat(np.arange(len(TIERS)), COUNTS)

names = []
seen = [0] * len(TIERS)
for t in tier_of:
seen[t] += 1
names.append(f"{TIERS[t]}{seen[t]}")

p_base = np.array([0.30, 0.15, 0.05, 0.08, 0.35, 0.25])
v_base = np.array([2.0, 4.0, 10.0, 8.0, 1.0, 3.0])
p = p_base[tier_of] * rng.uniform(0.8, 1.2, N)
v = v_base[tier_of] * rng.uniform(0.8, 1.2, N)

link_prob = np.array([
[0.30, 0.90, 0.00, 0.10, 0.00, 0.60],
[0.90, 0.45, 0.85, 0.10, 0.30, 0.40],
[0.00, 0.85, 0.80, 0.10, 0.00, 0.00],
[0.10, 0.10, 0.10, 0.20, 0.10, 0.10],
[0.00, 0.30, 0.00, 0.10, 0.25, 0.00],
[0.60, 0.40, 0.00, 0.10, 0.00, 0.40],
])
rate_mean = np.array([
[1.0, 8.0, 0.0, 0.3, 0.0, 2.0],
[8.0, 3.0, 10.0, 0.3, 1.5, 1.5],
[0.0, 10.0, 6.0, 0.3, 0.0, 0.0],
[0.3, 0.3, 0.3, 0.5, 0.3, 0.3],
[0.0, 1.5, 0.0, 0.3, 0.5, 0.0],
[2.0, 1.5, 0.0, 0.3, 0.0, 1.0],
])

T = np.zeros((N, N))
for i in range(N):
for j in range(i + 1, N):
ti, tj = tier_of[i], tier_of[j]
if rng.random() < link_prob[ti, tj]:
T[i, j] = rate_mean[ti, tj] * rng.uniform(0.5, 1.5)
T[j, i] = T[i, j]
L = (T > 0).astype(float)

DELTA = 1.8
EPS = 0.08
C_SEG = 6.0
C_RULE = 1.0

R0 = float(p.sum() * v.sum() - (p * v).sum())
T_TOT = float(T.sum() / 2.0)
N_LINKS = int(L.sum() / 2)
OPS_MICRO = C_SEG * N + C_RULE * N_LINKS

# ------------------------------------------------------------
# 2. Cost model
# ------------------------------------------------------------
def evaluate(a):
a = np.asarray(a, dtype=int)
K = int(a.max()) + 1
Pseg = np.bincount(a, weights=p, minlength=K)
Vseg = np.bincount(a, weights=v, minlength=K)
cross = (a[:, None] != a[None, :]).astype(float)
risk_in = float((Pseg * Vseg).sum() - (p * v).sum())
risk_x = EPS * float((np.outer(p, v) * L * cross).sum())
risk = risk_in + risk_x
lat = DELTA * 0.5 * float((T * cross).sum()) / T_TOT
rules = 0.5 * float((L * cross).sum())
k_used = int(np.count_nonzero(np.bincount(a, minlength=K)))
ops = C_SEG * k_used + C_RULE * rules
return {"risk": risk, "lat": lat, "rules": rules, "ops": ops, "k": k_used,
"risk_n": risk / R0, "lat_n": lat / DELTA, "ops_n": ops / OPS_MICRO}

def total_cost(m, w):
return w[0] * m["risk_n"] + w[1] * m["lat_n"] + w[2] * m["ops_n"]

# ------------------------------------------------------------
# 3. Fast simulated annealing (O(1) cross-term update per move)
# ------------------------------------------------------------
def anneal(w, kmax, iters=15000, seed=0, t_start=0.05, t_end=1e-4):
wr, wl, wo = w
if kmax == 1:
return np.zeros(N, dtype=int)
rs = np.random.default_rng(seed)
C = (wr * EPS * (np.outer(p, v) + np.outer(v, p)) * L / R0
+ wl * T / T_TOT + wo * C_RULE * L / OPS_MICRO)
np.fill_diagonal(C, 0.0)
seg_cost = wo * C_SEG / OPS_MICRO
k_in = wr / R0

a = rs.integers(0, kmax, N)
onehot = np.zeros((N, kmax))
onehot[np.arange(N), a] = 1.0
S = C @ onehot
Pseg = np.bincount(a, weights=p, minlength=kmax).astype(float)
Vseg = np.bincount(a, weights=v, minlength=kmax).astype(float)
cnt = np.bincount(a, minlength=kmax).astype(int)

J = total_cost(evaluate(a), w)
best_J, best_a = J, a.copy()

nodes = rs.integers(0, N, iters)
shifts = rs.integers(1, kmax, iters)
us = rs.random(iters)
temps = t_start * (t_end / t_start) ** (np.arange(iters) / max(iters - 1, 1))

for it in range(iters):
i = int(nodes[it])
ai = int(a[i])
b = (ai + int(shifts[it])) % kmax
pi = p[i]
vi = v[i]
d = k_in * (pi * (Vseg[b] - Vseg[ai]) + vi * (Pseg[b] - Pseg[ai]) + 2.0 * pi * vi)
d += S[i, ai] - S[i, b]
if cnt[ai] == 1:
d -= seg_cost
if cnt[b] == 0:
d += seg_cost
if d <= 0.0 or us[it] < math.exp(-d / temps[it]):
a[i] = b
S[:, ai] -= C[:, i]
S[:, b] += C[:, i]
Pseg[ai] -= pi
Pseg[b] += pi
Vseg[ai] -= vi
Vseg[b] += vi
cnt[ai] -= 1
cnt[b] += 1
J += d
if J < best_J - 1e-12:
best_J = J
best_a = a.copy()
return best_a

def best_of(w, kmax, restarts=3, iters=15000, base_seed=0):
best_a, best_J = None, float("inf")
for r in range(restarts):
cand = anneal(w, kmax, iters=iters, seed=base_seed + 101 * r + kmax)
j = total_cost(evaluate(cand), w)
if j < best_J:
best_J, best_a = j, cand
return best_a, best_J

# ------------------------------------------------------------
# 4. Experiments
# ------------------------------------------------------------
W_MAIN = (0.5, 0.3, 0.2)
t0 = time.time()

k_list = list(range(1, 13))
sweep_a, sweep_J, sweep_m = {}, {}, {}
for k in k_list:
a_k, j_k = best_of(W_MAIN, k, restarts=4, iters=15000)
sweep_a[k], sweep_J[k], sweep_m[k] = a_k, j_k, evaluate(a_k)

k_best = min(k_list, key=lambda k: sweep_J[k])
a_opt = sweep_a[k_best]

scenarios = {
"Flat (1 segment)": np.zeros(N, dtype=int),
"Tier-based (6)": tier_of.copy(),
"Micro (1 per asset)": np.arange(N),
"Optimized (SA)": a_opt,
}
metrics = {name: evaluate(a) for name, a in scenarios.items()}

wl_grid = np.linspace(0.05, 0.70, 8)
Z = np.zeros((len(wl_grid), len(k_list)))
pareto = []
for iy, wl in enumerate(wl_grid):
w = (0.8 - wl, wl, 0.2)
for ix, k in enumerate(k_list):
a_k, j_k = best_of(w, k, restarts=2, iters=12000, base_seed=7)
Z[iy, ix] = j_k
m = evaluate(a_k)
pareto.append((m["risk_n"], m["lat_n"], m["ops_n"], wl))
pareto = np.array(pareto)

elapsed = time.time() - t0

# ------------------------------------------------------------
# 5. Console report
# ------------------------------------------------------------
print("=" * 78)
print("Network segmentation design: problem summary")
print("=" * 78)
print(f"Assets: {N} | communicating pairs: {N_LINKS} | total traffic: {T_TOT:.1f} kreq/s")
print(f"Weights (risk, latency, ops) = {W_MAIN}")
print()
print(f"{'Scenario':<22}{'K':>4}{'Risk':>10}{'Lat[ms]':>10}{'Rules':>8}{'Ops':>9}{'J':>9}")
print("-" * 72)
for name, m in metrics.items():
print(f"{name:<22}{m['k']:>4}{m['risk']:>10.2f}{m['lat']:>10.3f}"
f"{m['rules']:>8.0f}{m['ops']:>9.1f}{total_cost(m, W_MAIN):>9.4f}")
print()
print("Optimal number of segments versus the upper limit K_max")
print(f"{'K_max':>6}{'K_used':>8}{'J':>10}")
for k in k_list:
print(f"{k:>6}{sweep_m[k]['k']:>8}{sweep_J[k]:>10.4f}")
print()
print("Optimized segments")
used = sorted(np.unique(a_opt), key=lambda s: (-np.sum(a_opt == s), s))
rank_of = {int(s): r for r, s in enumerate(used)}
node_rank = np.array([rank_of[int(x)] for x in a_opt])
for idx, s in enumerate(used, start=1):
members = [names[i] for i in range(N) if a_opt[i] == s]
print(f" Segment {idx} ({len(members)} assets): {', '.join(members)}")
print()
red_vs_flat = 100 * (1 - metrics["Optimized (SA)"]["risk"] / metrics["Flat (1 segment)"]["risk"])
ops_vs_micro = 100 * (1 - metrics["Optimized (SA)"]["ops"] / metrics["Micro (1 per asset)"]["ops"])
print(f"Risk reduction vs flat network : {red_vs_flat:.1f} %")
print(f"Operational load reduction vs micro : {ops_vs_micro:.1f} %")
print(f"Total computation time : {elapsed:.1f} s")

# ------------------------------------------------------------
# 6. Visualization (single figure)
# ------------------------------------------------------------
fig = plt.figure(figsize=(21, 13))
gs = GridSpec(2, 3, figure=fig)
tier_colors = plt.cm.tab10(np.arange(len(TIERS)))

# (a) Sorted traffic matrix
ax1 = fig.add_subplot(gs[0, 0])
order = np.lexsort((tier_of, node_rank))
Ts = T[np.ix_(order, order)]
im = ax1.imshow(Ts, cmap="magma", interpolation="nearest")
sizes = [int(np.sum(node_rank == r)) for r in range(len(used))]
pos = 0
for sz in sizes[:-1]:
pos += sz
ax1.axhline(pos - 0.5, color="cyan", lw=1.2)
ax1.axvline(pos - 0.5, color="cyan", lw=1.2)
ax1.set_xticks(range(N))
ax1.set_yticks(range(N))
ax1.set_xticklabels([names[i] for i in order], rotation=90, fontsize=6)
ax1.set_yticklabels([names[i] for i in order], fontsize=6)
ax1.set_title("(a) Traffic matrix sorted by optimized segment")
fig.colorbar(im, ax=ax1, fraction=0.046, pad=0.04, label="kreq/s")

# (b) Segment map
ax2 = fig.add_subplot(gs[0, 1])
seg_ids = list(used)
Ks = len(seg_ids)
coords = np.zeros((N, 2))
for idx, s in enumerate(seg_ids):
members = np.where(a_opt == s)[0]
ang = 2 * np.pi * idx / Ks
center = np.array([5.0 * np.cos(ang + 0.6), 5.0 * np.sin(ang + 0.6)])
rad = 0.0 if len(members) == 1 else 0.30 + 0.07 * len(members)
for q, m_i in enumerate(members):
th = 2 * np.pi * q / len(members)
coords[m_i] = center + rad * np.array([np.cos(th), np.sin(th)])
ax2.add_patch(plt.Circle(center, rad + 0.45, fill=False, ec="white", ls="--", lw=0.9))
ax2.text(center[0], center[1] - rad - 0.75, f"Seg {idx + 1}", ha="center", fontsize=8, color="white")
for i in range(N):
for j in range(i + 1, N):
if L[i, j] > 0:
same = a_opt[i] == a_opt[j]
ax2.plot([coords[i, 0], coords[j, 0]], [coords[i, 1], coords[j, 1]],
color="#888888" if same else "#ff4d4d",
alpha=0.35 if same else 0.45,
lw=min(0.3 + 0.25 * T[i, j], 2.5), zorder=1)
ax2.scatter(coords[:, 0], coords[:, 1], c=tier_colors[tier_of], s=70,
edgecolors="white", linewidths=0.5, zorder=3)
handles = [Line2D([0], [0], marker="o", color="none", markerfacecolor=tier_colors[t],
markersize=8, label=TIERS[t]) for t in range(len(TIERS))]
handles.append(Line2D([0], [0], color="#ff4d4d", lw=2, label="Cross-segment flow"))
ax2.legend(handles=handles, loc="upper right", fontsize=7, ncol=2)
ax2.set_aspect("equal")
ax2.set_xlim(-7.4, 7.4)
ax2.set_ylim(-7.4, 7.4)
ax2.axis("off")
ax2.set_title("(b) Optimized segmentation map")

# (c) Cost breakdown
ax3 = fig.add_subplot(gs[0, 2])
labels = list(metrics.keys())
comp = np.array([[W_MAIN[0] * metrics[n]["risk_n"],
W_MAIN[1] * metrics[n]["lat_n"],
W_MAIN[2] * metrics[n]["ops_n"]] for n in labels])
bottom = np.zeros(len(labels))
for c_idx, (lab, col) in enumerate(zip(["Risk", "Latency", "Operations"],
["#ff6b6b", "#4dd0e1", "#ffd54f"])):
ax3.bar(range(len(labels)), comp[:, c_idx], bottom=bottom, color=col, label=lab, edgecolor="black")
bottom += comp[:, c_idx]
for x, tot in enumerate(bottom):
ax3.text(x, tot + 0.01, f"{tot:.3f}", ha="center", fontsize=9)
ax3.set_xticks(range(len(labels)))
ax3.set_xticklabels([l.replace(" (", "\n(") for l in labels], fontsize=8)
ax3.set_ylabel("Weighted cost J")
ax3.set_ylim(0, bottom.max() * 1.15)
ax3.legend(fontsize=8)
ax3.set_title("(c) Cost breakdown by design")

# (d) K_max sweep
ax4 = fig.add_subplot(gs[1, 0])
ax4.plot(k_list, [sweep_J[k] for k in k_list], "o-", color="white", lw=2, label="Total J")
ax4.plot(k_list, [W_MAIN[0] * sweep_m[k]["risk_n"] for k in k_list], "s--", color="#ff6b6b", label="Risk term")
ax4.plot(k_list, [W_MAIN[1] * sweep_m[k]["lat_n"] for k in k_list], "^--", color="#4dd0e1", label="Latency term")
ax4.plot(k_list, [W_MAIN[2] * sweep_m[k]["ops_n"] for k in k_list], "d--", color="#ffd54f", label="Operations term")
ax4.axvline(k_best, color="lime", ls=":", lw=1.5)
ax4.text(k_best + 0.15, ax4.get_ylim()[1] * 0.9, f"best K_max = {k_best}", color="lime", fontsize=9)
ax4.set_xlabel("Upper limit of segments K_max")
ax4.set_ylabel("Weighted cost")
ax4.set_xticks(k_list)
ax4.grid(alpha=0.25)
ax4.legend(fontsize=8)
ax4.set_title("(d) Trade-off versus the number of segments")

# (e) 3D surface
ax5 = fig.add_subplot(gs[1, 1], projection="3d")
KX, WY = np.meshgrid(np.array(k_list), wl_grid)
ax5.plot_surface(KX, WY, Z, cmap="viridis", edgecolor="none", alpha=0.9)
kb_idx = np.argmin(Z, axis=1)
ax5.scatter(np.array(k_list)[kb_idx], wl_grid, Z[np.arange(len(wl_grid)), kb_idx],
color="red", s=35, depthshade=False, label="Optimal K_max")
for _ax in (ax5,):
for _axis in (_ax.xaxis, _ax.yaxis, _ax.zaxis):
_axis.set_pane_color((0.08, 0.08, 0.10, 1.0))
ax5.set_xlabel("K_max")
ax5.set_ylabel("Latency weight")
ax5.set_zlabel("Optimal J")
ax5.view_init(elev=28, azim=-125)
ax5.legend(fontsize=8)
ax5.set_title("(e) Optimal cost surface")

# (f) 3D Pareto scatter
ax6 = fig.add_subplot(gs[1, 2], projection="3d")
sc = ax6.scatter(pareto[:, 0], pareto[:, 1], pareto[:, 2], c=pareto[:, 3], cmap="plasma",
s=22, alpha=0.8, depthshade=False)
markers = {"Flat (1 segment)": "*", "Tier-based (6)": "P", "Micro (1 per asset)": "X", "Optimized (SA)": "D"}
for name, m in metrics.items():
ax6.scatter([m["risk_n"]], [m["lat_n"]], [m["ops_n"]], marker=markers[name], s=170,
color="white" if name != "Optimized (SA)" else "lime",
edgecolors="black", depthshade=False, label=name)
for _axis in (ax6.xaxis, ax6.yaxis, ax6.zaxis):
_axis.set_pane_color((0.08, 0.08, 0.10, 1.0))
ax6.set_xlabel("Risk (normalized)")
ax6.set_ylabel("Latency (normalized)")
ax6.set_zlabel("Ops (normalized)")
ax6.view_init(elev=22, azim=-60)
ax6.legend(fontsize=7, loc="upper left")
fig.colorbar(sc, ax=ax6, fraction=0.03, pad=0.08, label="Latency weight")
ax6.set_title("(f) Pareto cloud in risk-latency-ops space")

fig.suptitle("Optimal Network Segmentation: Security vs. Latency vs. Operational Load", fontsize=16)
fig.subplots_adjust(left=0.04, right=0.95, top=0.93, bottom=0.06, wspace=0.22, hspace=0.25)
plt.show()

3. Code Walkthrough

Section 1: building the network

The 30 assets are generated tier by tier. Every tier has a base compromise probability and asset value, with ±20% jitter per host. Databases are valuable but rarely compromised ($p=0.05$, $v=10$). IoT devices are the opposite: very exposed and nearly worthless on their own ($p=0.35$, $v=1$), which makes them ideal pivot points for an attacker.

The traffic matrix is sampled from two tier-level tables. link_prob gives the chance that two assets talk, and rate_mean gives their mean rate. Web–App and App–DB are dense and heavy, as in a classic three-tier architecture. Admin hosts touch everything at a trickle, and IoT devices only talk to a few App servers.

Section 2: the cost model

evaluate() computes the three terms directly from the formulas above. The in-segment risk uses a neat identity: $\sum_{i,j \in s,, i \ne j} p_i v_j = P_s V_s - \sum_{i\in s} p_i v_i$, where $P_s$ and $V_s$ are segment totals. A single np.bincount therefore replaces a double loop. The cross-segment terms use the mask cross, which is 1 where two assets sit in different segments.

evaluate() is the reference implementation. It is $O(N^2)$ per call and is used for reporting and for the initial state of the annealer.

Section 3: the fast simulated annealing

A naive annealer would call evaluate() at every move, which costs $O(N^2)$. Instead we compute the change in cost $\Delta J$ incrementally. All pairwise cross-segment costs (risk leak, latency, rules) are folded into one symmetric matrix $C_{ij}$ so that

$$
J_{\mathrm{cross}} = \sum_{i<j} C_{ij}, x_{ij}.
$$

We maintain a table $S_{i,s} = \sum_{j:,a_j = s} C_{ij}$, the total cross cost between asset $i$ and everything currently in segment $s$. Moving asset $i$ from segment $a$ to segment $b$ then changes the cross cost by exactly

$$
\Delta J_{\mathrm{cross}} = S_{i,a} - S_{i,b}.
$$

Neighbors left behind in $a$ become cross-segment, and neighbors in $b$ stop being cross-segment. That is an $O(1)$ lookup. After an accepted move, only two columns of $S$ are updated, which is $O(N)$.

The in-segment risk term is quadratic in the segment totals. Moving asset $i$ changes it by

$$
\Delta J_{\mathrm{in}} = \frac{w_r}{R_0}\Big[p_i (V_b - V_a) + v_i (P_b - P_a) + 2 p_i v_i\Big],
$$

where $P$ and $V$ are segment totals before the move. The segment-count penalty is added or removed when a segment becomes empty or newly occupied. A move is accepted with the Metropolis rule $\min(1, e^{-\Delta J / \tau})$, with the temperature $\tau$ cooled geometrically from 0.05 to $10^{-4}$.

best_of() runs several independent restarts with different seeds and keeps the best, which makes the result robust against unlucky starts.

Section 4: experiments

Three experiments are run:

  1. $K_{\max}$ sweep. For $K_{\max} = 1,\dots,12$ we search for the best partition using at most $K_{\max}$ segments, with the main weights $(0.5, 0.3, 0.2)$. The best result over all $K_{\max}$ is our optimized design.
  2. Baselines. A flat network, a classic tier-based design with six zones, and full micro-segmentation with one segment per asset.
  3. Weight sweep. For eight latency weights $w_l \in [0.05, 0.70]$, with $w_o = 0.2$ fixed and $w_r = 0.8 - w_l$, we repeat the $K_{\max}$ sweep. This produces the cost surface and the cloud of Pareto-style trade-off points.

Sections 5 and 6: report and figure

The console report prints a comparison table, the $K_{\max}$ sweep, the member list of each optimized segment, and the headline improvements. All six panels are drawn in one figure, so a single plt.show() produces everything.

4. Console Output

==============================================================================
Network segmentation design: problem summary
==============================================================================
Assets: 30 | communicating pairs: 130 | total traffic: 689.6 kreq/s
Weights (risk, latency, ops) = (0.5, 0.3, 0.2)

Scenario                 K      Risk   Lat[ms]   Rules      Ops        J
------------------------------------------------------------------------
Flat (1 segment)         1    746.75     0.000       0      6.0   0.5039
Tier-based (6)           6     90.72     1.651     109    145.0   0.4295
Micro (1 per asset)     30     15.80     1.800     130    310.0   0.5106
Optimized (SA)           5    253.95     0.166      44     74.0   0.2454

Optimal number of segments versus the upper limit K_max
 K_max  K_used         J
     1       1    0.5039
     2       2    0.3019
     3       3    0.2564
     4       4    0.2460
     5       5    0.2454
     6       5    0.2454
     7       5    0.2454
     8       5    0.2454
     9       5    0.2454
    10       5    0.2454
    11       5    0.2454
    12       5    0.2454

Optimized segments
  Segment 1 (17 assets): Web1, Web2, Web3, Web4, Web5, App1, App2, App3, App4, App5, App6, App7, App8, DB1, DB2, DB3, DB4
  Segment 2 (5 assets): IoT1, IoT2, IoT3, IoT4, IoT6
  Segment 3 (3 assets): IoT5, Corp3, Corp4
  Segment 4 (3 assets): Admin1, Admin2, Admin3
  Segment 5 (2 assets): Corp1, Corp2

Risk reduction vs flat network      : 66.0 %
Operational load reduction vs micro : 76.1 %
Total computation time              : 24.2 s

5. Result Figure

6. Interpreting the Results

Panel (a): the sorted traffic matrix

Sorting the traffic matrix by optimized segment makes the structure visible. Heavy traffic (bright cells) collapses into one large diagonal block: Web, App, and DB. Those three tiers exchange so much traffic that cutting between them would be costly in latency and rules. The remaining segments (IoT, Admin, Corp) have only faint traffic to the core, so boundaries around them are cheap. The optimizer found exactly the cut lines that sever the least traffic.

Panel (b): the segmentation map

The red edges are the only flows that cross a boundary. They are thin and sparse compared with the dense gray mass inside the big segment. Notably, the optimizer does not follow the tier structure blindly. It groups the IoT devices together, keeps the Admin hosts in their own segment, and splits the corporate machines depending on how they connect to the rest.

Panel (c): cost breakdown

Each design pays in a different currency:

  • Flat: almost the entire cost is risk, since one compromised host exposes everything.
  • Tier-based: risk is low, but the latency term is large because many Web–App and App–DB flows now cross a boundary.
  • Micro-segmentation: risk nearly vanishes, but latency and operations dominate. Every one of the 130 communicating pairs needs its own rule.
  • Optimized: a much smaller total than any baseline, because it accepts a moderate amount of in-segment risk in exchange for very little latency and operational cost.

Panel (d): the number of segments

Going from one segment to four, total cost drops steeply. Risk falls fast while latency and operations rise only slightly. Beyond five segments the curve is flat. Raising $K_{\max}$ further adds no benefit because the optimizer simply leaves the extra segments empty: any additional boundary would cost more in latency and rules than it saves in risk. The marginal-benefit logic of segmentation is visible here. The first few cuts are very valuable, and later cuts are not worth their price.

Panel (e): the 3D cost surface

The surface shows the optimal cost as a function of both the segment limit and the weight we place on latency. When latency is cheap (small $w_l$), the surface has a deep valley at moderate $K_{\max}$. As latency weighs more, the best designs shift toward fewer segments. The red markers trace the optimal $K_{\max}$ for each weight, showing how the right architecture depends on the business priorities rather than on a universal rule.

Panel (f): the Pareto cloud

Every dot is a solution found during the weight sweep, placed in the space of normalized risk, latency, and operations. The baselines sit at the extreme corners: flat at high risk, micro-segmentation at high operations. The optimized design lies near the point closest to the origin. Solutions on the surface between the extremes represent honest trade-offs. Moving along it, you can only buy lower risk by paying in latency or operations.

7. Conclusion

Network segmentation is not a matter of “more is safer”. With a simple quantitative model, simulated annealing finds a design with a handful of segments that is far better than both the flat network and the textbook tier-based layout under the stated weights. Two ingredients made it practical:

  • the incremental cost update, which reduces each annealing step from a full matrix evaluation to a constant-time lookup plus a column update,
  • multiple restarts, which guard against poor local optima.

The same framework extends naturally: add per-flow criticality, use measured latency per enforcement point, or constrain certain assets to share a segment for compliance. The objective keeps its shape, and only the cost matrix $C_{ij}$ changes.

Log Retention Optimization

Balancing Storage Cost, Incident Investigation, and Compliance

How long should you keep your logs? Keep them too briefly and you can’t investigate an incident or pass an audit. Keep them too long and storage costs balloon. This post turns the question into an optimization problem and solves it with Python.

The Problem

A company collects six log sources. Each has a daily volume $v_s$ (GB/day), a compliance minimum retention $T_{\min,s}$, and a weight $w_s$, the probability that the source is needed in an incident investigation.

Source GB/day Compliance minimum (days) Weight $w_s$
Firewall 400 90 0.5
DNS 250 90 0.4
Web/App 300 90 0.6
Auth/IdP 40 365 0.9
EDR 150 180 0.8
Cloud Audit 60 365 0.8

Each source has two decision variables:

  • $h$: days kept in hot storage (fast, searchable, expensive)
  • $T$: total retention in days, with the remaining $T-h$ days in cold archive (cheap, slow to restore)

The Model

Storage cost. With hot price $p_h$ and cold price $p_c$ in dollars per GB-month, the steady-state annual cost is

$$C_{\text{stor}}(h,T)=12,v,\bigl[p_h,h+p_c,(T-h)\bigr]$$

Dwell time. The time between a compromise and its detection, $D$, is modeled as log-normal:

$$F(d)=P(D\le d)=\Phi!\left(\frac{\ln d-\mu}{\sigma}\right),\qquad \mu=\ln(\text{median})$$

Investigation risk. With $\lambda$ relevant incidents per year, a restore penalty $r$ when the needed log is only in cold storage, and a loss $L$ when the log has already been deleted:

$$R(h,T)=\lambda,w,\Bigl[r,\bigl(F(T)-F(h)\bigr)+L,\bigl(1-F(T)\bigr)\Bigr]$$

Optimization problem.

$$\min_{h,T}; C(h,T)=C_{\text{stor}}(h,T)+R(h,T)\quad\text{s.t.}\quad T_{\min}\le T\le T_{\max},;; h_{\min}\le h\le T$$

Here $T_{\max}$ is a data-minimization cap and $h_{\min}$ is the shortest hot window needed for day-to-day monitoring.

Optimality conditions. Setting the partial derivatives to zero, with $f=F’$ the dwell-time density, gives

$$f(T^*)=\frac{12,v,p_c}{\lambda,w,(L-r)},\qquad f(h^*)=\frac{12,v,(p_h-p_c)}{\lambda,w,r}$$

The first condition says to keep logs in cold storage until the marginal storage cost equals the marginal reduction in expected loss. The second says the same for the hot window, using the restore penalty.

Speed-up. The cost separates into $C(h,T)=a(h)+b(T)$ with

$$a(h)=12,v,(p_h-p_c),h-\lambda,w,r,F(h),\qquad b(T)=12,v,p_c,T+\lambda,w,\bigl[(r-L)F(T)+L\bigr]$$

For each $T$, the best $h$ is the prefix minimum $\min_{h\le T}a(h)$, which one cumulative-minimum pass computes. This removes the two-dimensional search.

The Source Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.special import ndtr

plt.rcParams["font.size"] = 10
pd.set_option("display.width", 250)
pd.set_option("display.max_columns", None)
pd.set_option("display.float_format", lambda x: f"{x:,.1f}")

# ---------------------------------------------------------------
# 1. Parameters
# ---------------------------------------------------------------
SOURCES = pd.DataFrame({
"source": ["Firewall", "DNS", "Web/App", "Auth/IdP", "EDR", "Cloud Audit"],
"gb_per_day": [400, 250, 300, 40, 150, 60],
"t_min": [90, 90, 90, 365, 180, 365],
"weight": [0.5, 0.4, 0.6, 0.9, 0.8, 0.8],
}).set_index("source")

P_HOT, P_COLD = 0.10, 0.004 # $ per GB-month
LAMBDA, LOSS, RESTORE = 10, 60_000, 8_000
MEDIAN_DWELL, SIGMA = 45, 1.2 # log-normal dwell time
MU = np.log(MEDIAN_DWELL)
H_MIN, H_MAX, T_MAX = 7, 365, 730 # days

# ---------------------------------------------------------------
# 2. Model
# ---------------------------------------------------------------
def dwell_cdf(d, mu=MU, sigma=SIGMA):
d = np.asarray(d, dtype=float)
return np.where(d > 0, ndtr((np.log(np.maximum(d, 1e-9)) - mu) / sigma), 0.0)

def cost_surface(v, w, H, T, loss=LOSS, mu=MU, sigma=SIGMA):
Fh, Ft = dwell_cdf(H, mu, sigma), dwell_cdf(T, mu, sigma)
storage = 12 * v * (P_HOT * H + P_COLD * (T - H))
risk = LAMBDA * w * (RESTORE * (Ft - Fh) + loss * (1 - Ft))
return storage, risk

def optimize_bruteforce(v, w, t_min):
H = np.arange(0, H_MAX + 1, dtype=float)[:, None]
T = np.arange(t_min, T_MAX + 1, dtype=float)[None, :]
s, r = cost_surface(v, w, H, T)
tot = np.where((H < H_MIN) | (H > T), np.inf, s + r)
i, j = np.unravel_index(np.argmin(tot), tot.shape)
return int(H[i, 0]), int(T[0, j]), float(tot[i, j])

def optimize_fast(v, w, t_min, loss=LOSS, mu=MU, sigma=SIGMA):
h = np.arange(0, H_MAX + 1, dtype=float)
t = np.arange(t_min, T_MAX + 1, dtype=float)
Fh, Ft = dwell_cdf(h, mu, sigma), dwell_cdf(t, mu, sigma)
a = 12 * v * (P_HOT - P_COLD) * h - LAMBDA * w * RESTORE * Fh
a[h < H_MIN] = np.inf
b = 12 * v * P_COLD * t + LAMBDA * w * ((RESTORE - loss) * Ft + loss)
pm = np.minimum.accumulate(a)
arg = np.maximum.accumulate(np.where(a == pm, np.arange(h.size), 0))
k = np.minimum(t, H_MAX).astype(int)
tot = pm[k] + b
j = int(np.argmin(tot))
return int(h[arg[k[j]]]), int(t[j]), float(tot[j])

def evaluate(v, w, h, T):
s, r = cost_surface(v, w, np.array(float(h)), np.array(float(T)))
return float(s), float(r)

# ---------------------------------------------------------------
# 3. Solve and compare
# ---------------------------------------------------------------
t0 = time.perf_counter()
bf = {n: optimize_bruteforce(r.gb_per_day, r.weight, r.t_min) for n, r in SOURCES.iterrows()}
t_bf = time.perf_counter() - t0
t0 = time.perf_counter()
fs = {n: optimize_fast(r.gb_per_day, r.weight, r.t_min) for n, r in SOURCES.iterrows()}
t_fs = time.perf_counter() - t0
for n in SOURCES.index:
assert np.isclose(bf[n][2], fs[n][2]), n
print(f"Brute force: {t_bf*1000:.1f} ms | Fast: {t_fs*1000:.1f} ms | Speed-up: {t_bf/t_fs:.1f}x")

rows = []
for n, r in SOURCES.iterrows():
h, T, tot = fs[n]
s, rk = evaluate(r.gb_per_day, r.weight, h, T)
rows.append([n, r.gb_per_day, r.t_min, h, T, s, rk, s + rk])
opt = pd.DataFrame(rows, columns=["source", "GB/day", "T_min", "hot_days", "total_days",
"storage_$/yr", "risk_$/yr", "total_$/yr"]).set_index("source")
print(opt)
print(opt[["storage_$/yr", "risk_$/yr", "total_$/yr"]].sum())

policies = {
"A: Compliance min\n(all hot)": lambda r: (r.t_min, r.t_min),
"B: 365 d\n(all hot)": lambda r: (max(365, r.t_min), max(365, r.t_min)),
"C: Compliance min\n(30 d hot)": lambda r: (min(30, r.t_min), r.t_min),
"D: 365 d\n(30 d hot)": lambda r: (30, max(365, r.t_min)),
"E: Optimized\ntiered": None,
}
pol_rows = []
for name, f in policies.items():
st = rk = 0.0
for n, r in SOURCES.iterrows():
h, T = (fs[n][0], fs[n][1]) if f is None else f(r)
s, k = evaluate(r.gb_per_day, r.weight, h, T)
st += s
rk += k
pol_rows.append([name.replace("\n", " "), st, rk, st + rk])
pol = pd.DataFrame(pol_rows, columns=["policy", "storage_$/yr", "risk_$/yr", "total_$/yr"]).set_index("policy")
print(pol)

# Sensitivity: median dwell time x loss per blind incident
medians = np.linspace(15, 120, 8)
losses = np.linspace(20_000, 200_000, 8)
sens = np.zeros((medians.size, losses.size))
for i, m in enumerate(medians):
for j, L in enumerate(losses):
sens[i, j] = sum(optimize_fast(r.gb_per_day, r.weight, r.t_min, loss=L, mu=np.log(m))[2]
for _, r in SOURCES.iterrows())

# ---------------------------------------------------------------
# 4. Visualization (single figure)
# ---------------------------------------------------------------
fw = SOURCES.loc["Firewall"]
hg = np.arange(0, 121, 2, dtype=float)
tg = np.arange(fw.t_min, T_MAX + 1, 5, dtype=float)
HH, TT = np.meshgrid(hg, tg, indexing="ij")
s, r = cost_surface(fw.gb_per_day, fw.weight, HH, TT)
Z = np.where((HH < H_MIN) | (HH > TT), np.nan, s + r)
h_o, t_o, c_o = fs["Firewall"]

fig = plt.figure(figsize=(21, 14))

ax1 = fig.add_subplot(2, 3, 1, projection="3d")
ax1.plot_surface(HH, TT, Z / 1e3, cmap=cm.viridis, edgecolor="none", alpha=0.9)
ax1.scatter([h_o], [t_o], [c_o / 1e3], color="red", s=80, label="Optimum")
ax1.set_xlabel("Hot days h"); ax1.set_ylabel("Total days T"); ax1.set_zlabel("Cost (k$/yr)")
ax1.set_title("Firewall: total annual cost surface"); ax1.view_init(elev=28, azim=-125); ax1.legend()

ax2 = fig.add_subplot(2, 3, 2, projection="3d")
MM, LL = np.meshgrid(medians, losses / 1e3, indexing="ij")
ax2.plot_surface(MM, LL, sens / 1e3, cmap=cm.plasma, edgecolor="none")
ax2.set_xlabel("Median dwell time (days)"); ax2.set_ylabel("Loss per blind incident (k$)")
ax2.set_zlabel("Optimal total (k$/yr)"); ax2.set_title("Sensitivity of the optimal total cost")
ax2.view_init(elev=25, azim=-130)

ax3 = fig.add_subplot(2, 3, 3)
cs = ax3.contourf(HH, TT, Z / 1e3, levels=25, cmap="viridis")
fig.colorbar(cs, ax=ax3, label="k$/yr")
ax3.plot(h_o, t_o, "r*", markersize=16, label=f"Optimum (h={h_o}, T={t_o})")
ax3.set_xlabel("Hot days h"); ax3.set_ylabel("Total days T")
ax3.set_title("Firewall: contour view"); ax3.legend(loc="upper right")

ax4 = fig.add_subplot(2, 3, 4)
tt = np.arange(fw.t_min, T_MAX + 1, dtype=float)
s4, r4 = cost_surface(fw.gb_per_day, fw.weight, np.array(float(h_o)), tt)
ax4.plot(tt, s4 / 1e3, label="Storage")
ax4.plot(tt, r4 / 1e3, label="Expected risk")
ax4.plot(tt, (s4 + r4) / 1e3, "k", lw=2, label="Total")
ax4.axvline(t_o, color="red", ls="--", label=f"T* = {t_o}")
ax4.set_xlabel("Total retention T (days)"); ax4.set_ylabel("k$/yr")
ax4.set_title(f"Firewall: cost trade-off at h = {h_o}"); ax4.legend(); ax4.grid(alpha=0.3)

ax5 = fig.add_subplot(2, 3, 5)
x = np.arange(len(opt)); wd = 0.27
ax5.bar(x - wd, opt["hot_days"], wd, label="Hot days h*")
ax5.bar(x, opt["total_days"], wd, label="Total days T*")
ax5.bar(x + wd, opt["T_min"], wd, label="Compliance minimum")
ax5.set_xticks(x); ax5.set_xticklabels(opt.index, rotation=20)
ax5.set_ylabel("Days"); ax5.set_title("Optimal retention per source"); ax5.legend(); ax5.grid(alpha=0.3, axis="y")

ax6 = fig.add_subplot(2, 3, 6)
xp = np.arange(len(pol))
ax6.bar(xp, pol["storage_$/yr"] / 1e3, label="Storage")
ax6.bar(xp, pol["risk_$/yr"] / 1e3, bottom=pol["storage_$/yr"] / 1e3, label="Expected risk")
ax6.set_xticks(xp); ax6.set_xticklabels(list(policies.keys()), fontsize=8)
ax6.set_ylabel("k$/yr"); ax6.set_title("Policy comparison (all sources)"); ax6.legend(); ax6.grid(alpha=0.3, axis="y")
for i, v in enumerate(pol["total_$/yr"] / 1e3):
ax6.text(i, v, f"{v:,.0f}", ha="center", va="bottom", fontsize=9)

plt.tight_layout()
plt.show()

Code Walkthrough

1. Parameters

SOURCES holds the six log sources with daily volume, compliance minimum, and investigation weight. The constants define the economics and the threat model:

  • P_HOT and P_COLD are storage prices in dollars per GB-month.
  • LAMBDA is the number of relevant incidents per year.
  • RESTORE is the extra cost when the needed log must be rehydrated from cold storage (delay and engineering time).
  • LOSS is the cost of a “blind” incident, where the logs were already deleted.
  • MEDIAN_DWELL and SIGMA define the log-normal dwell-time distribution. The heavy right tail is what makes long retention valuable.
  • H_MIN, H_MAX, and T_MAX are the search bounds. T_MAX acts as a data-minimization cap, so the optimizer can never recommend keeping data forever.

2. Model

dwell_cdf is the log-normal CDF, computed with scipy.special.ndtr, the standard normal CDF. The np.where guard returns 0 for non-positive days, so the function accepts scalars, vectors, and 2D grids without special cases.

cost_surface implements the two cost components from the formulas above. It relies on NumPy broadcasting: pass H as a column and T as a row, and it returns the whole $(h,T)$ cost grid at once with no loops.

optimize_bruteforce evaluates every feasible $(h,T)$ pair and masks infeasible ones (h < H_MIN or h > T) with infinity. It is simple and obviously correct, so it serves as the reference.

optimize_fast uses the separable structure $C=a(h)+b(T)$:

  • np.minimum.accumulate(a) gives the best hot window for every possible $T$ in one pass.
  • The arg line recovers which $h$ achieved that minimum. Because a == pm is true wherever a new running minimum appears, taking a cumulative maximum over those indices tracks the latest minimizer.
  • The final argmin over pm[k] + b is a one-dimensional search.

The assert in the solving section checks that both methods find the same minimum cost for every source, so the fast version is verified rather than assumed.

evaluate returns the storage and risk components for any chosen policy, which makes policy comparisons easy.

3. Solving and Comparing

The script times both solvers and prints the speed-up. It then builds the optimal-policy table and compares five strategies:

  • A: the compliance minimum, all hot.
  • B: one year for everything, all hot.
  • C: the compliance minimum with a 30-day hot window.
  • D: one year with a 30-day hot window.
  • E: the optimized tiered policy.

The last block re-solves the whole problem for an 8×8 grid of median dwell times and blind-incident losses. This is only practical because the fast solver is so cheap.

4. Visualization

All six panels are drawn into one figure and shown with a single plt.show(). They are explained in the results section below.

Execution Results

Brute force: 29.9 ms | Fast: 3.6 ms | Speed-up: 8.3x
             GB/day  T_min  hot_days  total_days  storage_$/yr  risk_$/yr  total_$/yr
source                                                                               
Firewall      400.0   90.0        26         535      22,252.8   32,133.2    54,386.0
DNS           250.0   90.0        36         586      17,400.0   21,734.6    39,134.6
Web/App       300.0   90.0        46         634      25,027.2   27,937.7    52,964.9
Auth/IdP       40.0  365.0       218         730      11,447.0   11,522.7    22,969.8
EDR           150.0  180.0        99         730      22,363.2   20,565.4    42,928.6
Cloud Audit    60.0  365.0       168         730      13,714.6   12,922.6    26,637.2
storage_$/yr   112,204.8
risk_$/yr      126,816.3
total_$/yr     239,021.1
dtype: float64
                              storage_$/yr  risk_$/yr  total_$/yr
policy                                                           
A: Compliance min (all hot)      178,800.0  354,460.0   533,260.0
B: 365 d (all hot)               525,600.0   97,315.1   622,915.1
C: Compliance min (30 d hot)      48,624.0  509,527.0   558,151.0
D: 365 d (30 d hot)               62,496.0  286,668.1   349,164.1
E: Optimized tiered              112,204.8  126,816.3   239,021.1

Interpreting the Results

Optimal policy per source. The optimizer recommends very different policies for different sources:

  • Firewall, DNS, and Web/App are high-volume with modest weights. Their hot windows are short (about 26, 36, and 46 days), but their total retention is far beyond the 90-day minimum, at roughly 535, 586, and 634 days. The cold tier costs so little that extending retention is cheap insurance against the long tail of slow detections.
  • Auth/IdP, EDR, and Cloud Audit are low-volume, high-weight sources. They hit the 730-day cap, and their hot windows are much longer (about 218, 99, and 168 days) because a restore delay on these logs is expensive and keeping them hot is cheap in absolute terms.

Policy comparison. The optimized tiered policy costs about $239k per year in total. The alternatives are much worse:

Policy Total ($/yr)
A: compliance minimum, all hot about 533k
B: 365 days, all hot about 623k
C: compliance minimum, 30 days hot about 558k
D: 365 days, 30 days hot about 349k
E: optimized tiered about 239k

Policies A and C look cheap on storage but carry very large expected risk. Policy B cuts risk but pays heavily for hot storage. Even policy D, which is a reasonable-sounding rule of thumb, costs roughly 46% more than the optimum. The lesson is that neither “keep less” nor “keep everything hot” wins. The savings come from matching each source’s tiering and retention to its volume, value, and the shape of your detection-time distribution.

3D cost surface (top left) and contour view (top right). For the firewall, the surface is a valley with a clear minimum, marked in red.

  • Along the $T$ axis, cost falls steeply at first as the risk of blind incidents drops, then rises slowly as cold storage accumulates.
  • Along the $h$ axis, cost rises gently with hot storage but also falls as restore penalties are avoided, which produces the shallow bowl around the optimum.
  • Cutting off the left side of the surface is the minimum hot window, and cutting off the diagonal is the rule $h\le T$.
  • The contour view shows the same minimum from above. The valley is much flatter along $T$ than along $h$, so a retention period near the optimum is forgiving while the hot window deserves more careful tuning.

Cost trade-off curve (bottom left). Storage cost grows linearly with $T$, while expected risk decays as the log-normal tail is covered. Their sum is the black curve, whose minimum sits at $T^*$. This is a direct picture of the condition $f(T^*)=12,v,p_c/\bigl(\lambda,w,(L-r)\bigr)$: the optimum is where the density of “late discoveries” no longer justifies another day of archive.

Sensitivity surface (top middle). The optimal total cost rises with both the median dwell time and the loss per blind incident.

  • Slower detection pushes more incidents into the tail, so more retention and more cost are needed.
  • Larger blind-incident losses make long retention more valuable, so the optimal spend grows.
  • Reducing your median detection time lowers the cost of the optimal policy, which is a business case for investing in detection engineering.

Per-source bar chart (bottom middle). The compliance minimum is a floor, not a target. For every source, the optimal total retention is at or above it, and for most it is far above. The chart also shows how much the hot window varies, from under a month to over seven months.

Stacked policy chart (bottom right). Each bar splits storage from expected risk. Policies at the extremes are dominated by one component, and the optimized policy balances the two.

Conclusion

Log retention is not a single number to apply across the board. By modeling detection delay statistically and pricing storage tiers explicitly, the problem becomes a small, well-behaved optimization with a clean closed-form intuition. The separable structure also lets you solve it fast enough to run sensitivity analyses over the assumptions that matter most.

Optimizing the Scope of MFA Enforcement

Balancing Usability and Security with Python

Multi-factor authentication (MFA) is one of the most effective defenses against credential theft. But “MFA for everyone, everywhere” is not free. Every prompt costs seconds of employee time, increases help-desk tickets, and pushes frustrated users toward workarounds. Meanwhile, some user groups hold almost nothing worth stealing.

In this article we turn the question “Who should be forced to use MFA?“ into a small optimization problem and solve it exactly in Python.

1. The Problem Setting

A company has 14 user segments, such as executives, IT administrators, finance staff, contractors, and cafeteria staff. For each segment $i$ we know:

  • $n_i$: number of users
  • $p_i$: annual probability that a user’s credentials get compromised
  • $L_i$: average loss per compromise (dollars)
  • $e_i$: the fraction of compromises that MFA prevents (effectiveness)
  • $f_i$: annual usability cost per user (lost time, support tickets, workarounds)
  • $D_i$: one-time-per-year deployment and operation cost for the segment

The decision variable is $x_i \in {0,1}$, where $x_i = 1$ means MFA is enforced for segment $i$.

Benefit (expected annual loss avoided):

$$
B_i = n_i , p_i , L_i , e_i
$$

Cost (usability plus deployment):

$$
C_i = n_i , f_i + D_i
$$

Problem 1: unconstrained optimum. Maximize the net value:

$$
\max_{x \in {0,1}^N} \sum_{i=1}^{N} x_i ,(B_i - C_i)
$$

Problem 2: budget-constrained optimum. Maximize the avoided loss while keeping the total friction and deployment cost within a budget $K$:

$$
\max_{x \in {0,1}^N} \sum_{i=1}^{N} x_i B_i \quad \text{s.t.} \quad \sum_{i=1}^{N} x_i C_i \le K
$$

Problem 2 is a classic 0-1 knapsack problem. With $N = 14$ there are only $2^{14} = 16{,}384$ possible policies, so we can enumerate all of them and obtain the provably exact optimum.

The residual risk after applying a policy $x$ is:

$$
R(x) = \sum_{i=1}^{N} n_i p_i L_i ;-; \sum_{i=1}^{N} x_i B_i
$$

For the 3D visualization, we also use the per-user view. Let $r_i = p_i L_i e_i$ be the expected loss avoided per user. Enforcing MFA on one user is worthwhile when:

$$
v = r - f > 0
$$

2. Source Code

The entire program is a single script. To avoid a slow Python loop over $2^{14}$ policies, the enumeration is fully vectorized: every policy is encoded as a bitmask, expanded into a binary matrix $M \in {0,1}^{2^N \times N}$, and evaluated with two matrix-vector products.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # noqa: F401

np.set_printoptions(suppress=True)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)

# ------------------------------------------------------------------
# 1. Problem data: 14 user segments
# ------------------------------------------------------------------
rows = [
# name, users, p, loss, friction/user/yr, deploy cost, MFA effectiveness
("Executives", 40, 0.060, 250000, 60, 1500, 0.97),
("IT administrators", 60, 0.050, 400000, 70, 2000, 0.98),
("Finance staff", 120, 0.035, 90000, 55, 3000, 0.95),
("HR staff", 80, 0.030, 60000, 50, 2500, 0.95),
("Contractors", 150, 0.050, 30000, 90, 4000, 0.92),
("Engineers", 400, 0.020, 40000, 45, 5000, 0.95),
("Customer support", 300, 0.030, 15000, 40, 4000, 0.90),
("Sales", 250, 0.025, 12000, 50, 3500, 0.90),
("Marketing", 100, 0.015, 6000, 45, 2500, 0.90),
("Legal", 30, 0.030, 80000, 55, 1200, 0.95),
("Back-office clerks", 350, 0.008, 2500, 60, 4500, 0.88),
("Warehouse kiosk staff", 500, 0.004, 800, 35, 6000, 0.80),
("Interns", 80, 0.010, 1500, 65, 2000, 0.88),
("Cafeteria staff", 120, 0.003, 500, 30, 2500, 0.75),
]
df = pd.DataFrame(rows, columns=["segment", "users", "p", "loss", "friction", "deploy", "eff"])
N = len(df)

# Benefit B_i = n_i * p_i * L_i * e_i (expected annual loss avoided)
# Cost C_i = n_i * f_i + D_i (usability cost + deployment cost)
df["benefit"] = df.users * df.p * df.loss * df.eff
df["cost"] = df.users * df.friction + df.deploy
df["net"] = df.benefit - df.cost
df["ratio"] = df.benefit / df.cost
df["risk_per_user"] = df.p * df.loss * df.eff

B = df.benefit.values
C = df.cost.values

# ------------------------------------------------------------------
# 2. Exact optimization by vectorized enumeration of all 2^N policies
# ------------------------------------------------------------------
codes = np.arange(2 ** N, dtype=np.int64)
M = ((codes[:, None] >> np.arange(N)) & 1).astype(np.float64) # (2^N, N)
tot_benefit = M @ B
tot_cost = M @ C

# (a) Unconstrained optimum: maximize sum(B - C)
best_free = int(np.argmax(tot_benefit - tot_cost))

# (b) Budget-constrained optimum: maximize benefit s.t. cost <= budget
BUDGET = 150000.0
feasible = tot_cost <= BUDGET
best_budget = int(np.argmax(np.where(feasible, tot_benefit, -np.inf)))

# (c) Greedy heuristic by benefit/cost ratio under the same budget
order = np.argsort(-df.ratio.values)
greedy = np.zeros(N)
spent = 0.0
for i in order:
if spent + C[i] <= BUDGET and B[i] > C[i]:
greedy[i] = 1
spent += C[i]
greedy_benefit = float(greedy @ B)

df["enforce_free"] = M[best_free].astype(bool)
df["enforce_budget"] = M[best_budget].astype(bool)
df["enforce_greedy"] = greedy.astype(bool)

# (d) Baselines
all_mask = np.ones(N)
total_loss_no_mfa = float((df.users * df.p * df.loss).sum())

def summarize(name, mask):
b = float(mask @ B)
c = float(mask @ C)
return {"policy": name, "segments": int(mask.sum()),
"users_covered": int(mask @ df.users.values),
"benefit": b, "cost": c, "net": b - c,
"residual_loss": total_loss_no_mfa - b}

summary = pd.DataFrame([
summarize("No MFA", np.zeros(N)),
summarize("MFA for everyone", all_mask),
summarize("Unconstrained optimum", M[best_free]),
summarize(f"Budget optimum (<= ${BUDGET:,.0f})", M[best_budget]),
summarize("Greedy by ratio (same budget)", greedy),
])

print("=== Segment data and per-segment economics ===")
print(df[["segment", "users", "benefit", "cost", "net", "ratio",
"enforce_free", "enforce_budget"]].round(2).to_string(index=False))
print("\n=== Policy comparison ===")
print(summary.round(0).to_string(index=False))
print("\nGreedy / exact benefit ratio under budget: "
f"{greedy_benefit / tot_benefit[best_budget]:.4f}")
print("Enforced (budget optimum):", ", ".join(df.segment[df.enforce_budget]))
print("Not enforced (budget optimum):", ", ".join(df.segment[~df.enforce_budget]))

# ------------------------------------------------------------------
# 3. Pareto frontier: best benefit for each friction budget
# ------------------------------------------------------------------
budgets = np.linspace(0, tot_cost.max(), 120)
frontier = np.array([tot_benefit[tot_cost <= b].max() for b in budgets])
residual_frontier = total_loss_no_mfa - frontier

# ------------------------------------------------------------------
# 4. Sensitivity: number of segments worth enforcing
# ------------------------------------------------------------------
fric_mult = np.linspace(0.5, 3.0, 26)
eff_scale = np.linspace(0.5, 1.0, 26)
count_grid = np.zeros((len(eff_scale), len(fric_mult)))
for a, es in enumerate(eff_scale):
for b_, fm in enumerate(fric_mult):
eff = np.minimum(df.eff.values * es, 0.99)
ben = df.users.values * df.p.values * df.loss.values * eff
cst = df.users.values * df.friction.values * fm + df.deploy.values
count_grid[a, b_] = np.sum(ben > cst)

# ------------------------------------------------------------------
# 5. Surface: net value per user v(r, f) = r - f, r = expected loss avoided
# ------------------------------------------------------------------
r_axis = np.linspace(0, 400, 60)
f_axis = np.linspace(0, 400, 60)
R, F = np.meshgrid(r_axis, f_axis)
V = R - F

# ------------------------------------------------------------------
# 6. All figures in one image
# ------------------------------------------------------------------
fig = plt.figure(figsize=(21, 13))
fig.suptitle("Optimizing the Scope of MFA Enforcement", fontsize=20, fontweight="bold", y=0.98)

# (1) Net value per segment
ax1 = fig.add_subplot(2, 3, 1)
d1 = df.sort_values("net")
colors = ["#2a9d8f" if v > 0 else "#e76f51" for v in d1.net]
ax1.barh(d1.segment, d1.net / 1000, color=colors)
ax1.axvline(0, color="black", lw=1)
ax1.set_xlabel("Net value of enforcing MFA (k$ / year)")
ax1.set_title("(1) Net value per segment")
ax1.grid(axis="x", alpha=0.3)

# (2) Benefit vs cost
ax2 = fig.add_subplot(2, 3, 2)
c2 = np.where(df.enforce_budget, "#2a9d8f", "#e76f51")
ax2.scatter(df.cost / 1000, df.benefit / 1000, s=df.users / 2 + 30, c=c2,
alpha=0.8, edgecolors="black")
lim = max(df.cost.max(), df.benefit.max()) / 1000
ax2.plot([0.5, lim], [0.5, lim], "k--", lw=1, label="Benefit = Cost")
for _, r in df.iterrows():
ax2.annotate(r.segment, (r.cost / 1000, r.benefit / 1000), fontsize=7,
xytext=(4, 4), textcoords="offset points")
ax2.set_xscale("log"); ax2.set_yscale("log")
ax2.set_xlabel("Cost: friction + deployment (k$ / year)")
ax2.set_ylabel("Benefit: loss avoided (k$ / year)")
ax2.set_title("(2) Benefit vs cost (green = enforced under budget)")
ax2.legend(loc="upper left"); ax2.grid(alpha=0.3, which="both")

# (3) Pareto frontier
ax3 = fig.add_subplot(2, 3, 3)
ax3.plot(budgets / 1000, residual_frontier / 1000, color="#264653", lw=2,
label="Optimal frontier")
for name, mask, mk, col in [("Unconstrained optimum", M[best_free], "*", "#e9c46a"),
("Budget optimum", M[best_budget], "o", "#2a9d8f"),
("Greedy", greedy, "s", "#f4a261"),
("MFA for everyone", all_mask, "X", "#e76f51")]:
ax3.scatter(mask @ C / 1000, (total_loss_no_mfa - mask @ B) / 1000,
s=220, marker=mk, color=col, edgecolors="black", zorder=5, label=name)
ax3.axvline(BUDGET / 1000, color="gray", ls=":", label="Friction budget")
ax3.set_xlabel("Total MFA cost (k$ / year)")
ax3.set_ylabel("Residual expected loss (k$ / year)")
ax3.set_title("(3) Security vs usability trade-off")
ax3.legend(fontsize=8); ax3.grid(alpha=0.3)

# (4) 3D surface
ax4 = fig.add_subplot(2, 3, 4, projection="3d")
ax4.plot_surface(R, F, V, cmap="viridis", alpha=0.85, linewidth=0)
ax4.plot_surface(R, F, np.zeros_like(V), color="gray", alpha=0.25)
ax4.set_xlabel("Loss avoided r ($/user/yr)")
ax4.set_ylabel("Friction f ($/user/yr)")
ax4.set_zlabel("Net value v = r - f")
ax4.set_title("(4) Net value per user (gray plane: v = 0)")
ax4.view_init(elev=25, azim=-60)

# (5) 3D scatter of segments
ax5 = fig.add_subplot(2, 3, 5, projection="3d")
ec = np.where(df.enforce_free, "#2a9d8f", "#e76f51")
ax5.scatter(df.risk_per_user, df.friction, df.users, c=ec, s=90,
edgecolors="black", depthshade=True)
for _, r in df.iterrows():
ax5.text(r.risk_per_user, r.friction, r.users, r.segment, fontsize=6)
ax5.set_xlabel("Loss avoided r ($/user/yr)")
ax5.set_ylabel("Friction f ($/user/yr)")
ax5.set_zlabel("Users")
ax5.set_title("(5) Segments in 3D (green = enforce)")
ax5.view_init(elev=22, azim=-50)

# (6) Sensitivity heatmap
ax6 = fig.add_subplot(2, 3, 6)
im = ax6.imshow(count_grid, origin="lower", aspect="auto", cmap="YlGnBu",
extent=[fric_mult[0], fric_mult[-1], eff_scale[0], eff_scale[-1]])
cs = ax6.contour(fric_mult, eff_scale, count_grid, colors="white", linewidths=0.8)
ax6.clabel(cs, fmt="%d", fontsize=8)
ax6.scatter([1.0], [1.0], color="red", s=80, zorder=5, label="Baseline")
ax6.set_xlabel("Friction multiplier")
ax6.set_ylabel("MFA effectiveness scale")
ax6.set_title("(6) Segments worth enforcing")
ax6.legend(loc="lower left")
fig.colorbar(im, ax=ax6, label="Number of segments")

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

3. Code Walkthrough

Step 1: Defining the data

Each tuple in rows describes one segment. The values are deliberately realistic in shape: executives and IT administrators have a high compromise probability and an enormous loss per incident, while cafeteria and warehouse kiosk staff have tiny loss values. Contractors are interesting because they have a high compromise probability but also the highest friction ($90 per user per year), as they use unmanaged devices and many different tools.

The two core formulas are computed directly as pandas columns: benefit is $n_i p_i L_i e_i$ and cost is $n_i f_i + D_i$. We also derive net, ratio (benefit divided by cost), and risk_per_user, which is $r_i = p_i L_i e_i$ and is used later in the 3D plots.

Step 2: Exact optimization by vectorized enumeration

The key line is:

1
M = ((codes[:, None] >> np.arange(N)) & 1).astype(np.float64)

codes holds the integers $0, 1, \dots, 2^N - 1$. Shifting each integer right by $0, 1, \dots, N-1$ bits and masking with & 1 extracts every bit, so row $k$ of M is the binary representation of $k$. That row is exactly one enforcement policy $x$. The resulting matrix has $16{,}384$ rows and $14$ columns.

Then M @ B and M @ C compute $\sum_i x_i B_i$ and $\sum_i x_i C_i$ for all policies at once. This replaces a Python loop of $16{,}384 \times 14$ iterations with two BLAS calls that finish almost instantly, which is the speed-up for the time-consuming part.

  • best_free maximizes tot_benefit - tot_cost, which solves Problem 1.
  • best_budget masks out infeasible policies with $-\infty$ and maximizes the benefit, which solves Problem 2.

Because this is a brute-force search over the entire space, the answer is guaranteed to be the global optimum.

Step 3: Greedy heuristic

For larger organizations (say, 60 or more segments), $2^N$ enumeration is impossible. The standard fallback is a greedy rule: sort segments by benefit-to-cost ratio and add them while the budget allows. We include it as a comparison to show how close the cheap heuristic comes to the exact optimum. I also require $B_i > C_i$, so the heuristic never buys a segment that destroys value.

Step 4: Baselines and summary

summarize evaluates any policy mask and returns the number of segments, covered users, total benefit, total cost, net value, and residual expected loss $R(x)$. We compare five policies: no MFA, MFA for everyone, the unconstrained optimum, the budget optimum, and the greedy policy.

Step 5: Pareto frontier

For 120 budget levels from $0$ up to the cost of covering everyone, we find the maximum achievable benefit among all feasible policies. Converting this to residual loss gives the best possible security level for each amount of friction we are willing to accept.

Step 6: Sensitivity analysis

Real estimates of friction and effectiveness are uncertain. We sweep a friction multiplier (0.5x to 3x) and an effectiveness scale (0.5x to 1.0x) and count how many segments still satisfy $B_i > C_i$. This shows how robust the recommended scope is.

Step 7: Visualization

All six panels are drawn into one figure, so a single plt.show() produces the whole picture. Panels 4 and 5 are 3D plots created with projection="3d".

4. Console Output

=== Segment data and per-segment economics ===
              segment  users   benefit  cost       net  ratio  enforce_free  enforce_budget
           Executives     40  582000.0  3900  578100.0 149.23          True            True
    IT administrators     60 1176000.0  6200 1169800.0 189.68          True            True
        Finance staff    120  359100.0  9600  349500.0  37.41          True            True
             HR staff     80  136800.0  6500  130300.0  21.05          True            True
          Contractors    150  207000.0 17500  189500.0  11.83          True            True
            Engineers    400  304000.0 23000  281000.0  13.22          True            True
     Customer support    300  121500.0 16000  105500.0   7.59          True            True
                Sales    250   67500.0 16000   51500.0   4.22          True            True
            Marketing    100    8100.0  7000    1100.0   1.16          True            True
                Legal     30   68400.0  2850   65550.0  24.00          True            True
   Back-office clerks    350    6160.0 25500  -19340.0   0.24         False            True
Warehouse kiosk staff    500    1280.0 23500  -22220.0   0.05         False           False
              Interns     80    1056.0  7200   -6144.0   0.15         False            True
      Cafeteria staff    120     135.0  6100   -5965.0   0.02         False            True

=== Policy comparison ===
                       policy  segments  users_covered   benefit     cost       net  residual_loss
                       No MFA         0              0       0.0      0.0       0.0      3167980.0
             MFA for everyone        14           2580 3039031.0 170850.0 2868181.0       128949.0
        Unconstrained optimum        10           1530 3030400.0 108550.0 2921850.0       137580.0
 Budget optimum (<= $150,000)        13           2080 3037751.0 147350.0 2890401.0       130229.0
Greedy by ratio (same budget)        10           1530 3030400.0 108550.0 2921850.0       137580.0

Greedy / exact benefit ratio under budget: 0.9976
Enforced (budget optimum): Executives, IT administrators, Finance staff, HR staff, Contractors, Engineers, Customer support, Sales, Marketing, Legal, Back-office clerks, Interns, Cafeteria staff
Not enforced (budget optimum): Warehouse kiosk staff

5. Graph Output

6. Reading the Results

What the numbers say

Without any MFA, the expected annual loss across the whole company is about $3.17 million. Enforcing MFA for everyone avoids about $3.04 million of that, at a total cost of about $171 thousand, leaving a residual expected loss of about $129 thousand.

The unconstrained optimum is different. It enforces MFA in 10 of the 14 segments, covering 1,530 of 2,580 users. Its cost is only about $109 thousand, and its net value is about $2.92 million. That beats “MFA for everyone” (net about $2.87 million) by roughly $54 thousand per year, even though it leaves slightly more residual loss (about $138 thousand). The four skipped segments are back-office clerks, warehouse kiosk staff, interns, and cafeteria staff. For these groups the expected loss avoided is smaller than the friction and deployment cost, so forcing MFA destroys value.

Why the budget optimum enforces more segments

With a budget of $150 thousand, the exact optimum enforces 13 segments and skips only the warehouse kiosk staff. This looks surprising at first, because Problem 2 maximizes avoided loss, not net value. When the budget is not tight, adding a low-value segment still increases total benefit a little, so the solver keeps adding until the budget is nearly exhausted (about $147 thousand spent). If your goal is cost-effectiveness rather than maximum risk reduction, Problem 1 is the right formulation. If your goal is “reduce risk as much as the friction budget allows,” Problem 2 is the right one. Making this distinction explicit is an important part of a policy discussion.

The greedy heuristic reaches the same set as the unconstrained optimum, and its benefit is 99.76% of the exact budget optimum. For a large number of segments, this makes greedy a practical and nearly lossless substitute.

Panel (1): Net value per segment

The bar chart makes the asymmetry obvious. IT administrators and executives generate net values in the hundreds of thousands to over a million dollars. Engineers, finance staff, and contractors form a solid middle tier. Marketing sits barely above zero (a net of about $1.1 thousand, benefit-to-cost ratio of 1.16), which marks the boundary. The four red bars at the bottom are the segments where MFA enforcement costs more than it saves.

Panel (2): Benefit versus cost on a log-log scale

Points above the dashed diagonal are worth enforcing. The bubble size is proportional to the number of users. IT administrators and executives sit far above the line, with benefit-to-cost ratios of roughly 190 and 149. The low-risk, high-headcount groups (warehouse kiosk, back-office, cafeteria) sit below it. Marketing hugs the diagonal, which shows how sensitive the decision is for borderline groups.

Panel (3): The security-versus-usability frontier

The curve shows the lowest achievable residual loss for each cost level. It drops steeply at first, because the cheapest segments (administrators, executives, finance) remove most of the risk, and then flattens into a long tail. Beyond roughly $100 thousand, each additional dollar buys almost no extra risk reduction. The markers confirm this: the unconstrained optimum, the greedy policy, and the budget optimum all lie almost on the flat part of the curve, while “MFA for everyone” is the farthest to the right with essentially no extra security to show for it. This chart is the best one to show to management, since it turns the “security versus convenience” argument into a concrete curve.

Panel (4): 3D surface of per-user net value

The surface is the plane $v = r - f$ over the plane of expected loss avoided $r$ and friction $f$, with a translucent gray plane at $v = 0$. The decision boundary is the diagonal line $r = f$ where the surface crosses the gray plane. Everything where the surface is above the plane is worth protecting. The geometry shows that MFA scope should be set by the ratio of risk to friction, not by an organization-wide blanket rule.

Panel (5): Segments in 3D

Each segment is a point in the space of loss avoided per user, friction per user, and headcount. Green points are recommended for enforcement and red ones are not. The green points are spread along the high-risk axis, with administrators and executives at the far end, while red points collect near the origin of the risk axis, including the largest groups by headcount such as warehouse kiosk staff. The third axis makes a practical point: the red segments are not small, so skipping them removes a large amount of user-facing friction.

Panel (6): Sensitivity heatmap

The horizontal axis multiplies the friction estimate and the vertical axis scales MFA effectiveness. The color and white contour labels show how many segments still pass the $B_i > C_i$ test. The baseline (red dot) corresponds to the 10 segments found above. Moving right (friction becomes more painful) or down (MFA becomes less effective) shrinks the enforced scope, and the core group of high-risk segments stays enforced in almost every scenario. This means the headline recommendation, “always enforce MFA for administrators, executives, finance, HR, legal, and contractors,” is robust, while the decision for borderline groups such as marketing and sales depends heavily on how well you can estimate friction.

7. Takeaways

Casting MFA scope as a 0-1 optimization gives three practical lessons.

  1. Risk-based scoping beats blanket enforcement. Skipping only four low-risk segments raised the net value by about $54 thousand per year in this example while removing friction for 1,050 users.
  2. The marginal value curve is steep and then flat. Most of the risk reduction comes from a small number of high-value segments, so these should be enforced first and with the strongest factors (for example, phishing-resistant authenticators).
  3. Borderline groups need better data. The sensitivity map shows that the uncertain decisions are the ones near $B_i \approx C_i$. Measuring real friction costs through help-desk tickets and login time is more valuable than refining the risk model for groups that are already far from the boundary.

For larger organizations, replace the brute-force enumeration with a dynamic-programming knapsack solver or an integer programming library, and add constraints such as regulatory requirements that force MFA for certain segments regardless of economics.

Maximizing Population Between Upstream and Downstream Villages with Linear Programming

Sharing a River Fairly

When a river runs through several settlements, every liter taken upstream is a liter missing downstream. Each village also has to split its share between drinking water and irrigation for the food its people eat. In this post, I build a small but realistic model of this dilemma, solve it exactly with linear programming in Python, and visualize the results, including a 3D surface that shows how river flow and upstream limits shape the maximum sustainable population.

1. The Example Scenario

A single river passes through six settlements in this order:

Order Settlement Max population (thousand) Irrigation water per capita
1 Up-1 Highland 20 1.2
2 Up-2 Terrace 15 1.4
3 Up-3 Forest 10 1.6
4 Down-1 Plain 25 1.5
5 Down-2 Delta 20 1.8
6 Down-3 Coast 15 2.0

Other conditions:

  • Water is measured in abstract flow units, and population in thousands of people.
  • The river flow at the source is $Q = 100$.
  • Domestic water need per capita is $d = 0.2$.
  • A fraction $r = 0.3$ of the withdrawn water returns to the river, so only $(1-r)$ of each withdrawal is truly consumed.
  • An environmental minimum flow $E = 10$ must remain in the river at every point.
  • Every village must keep at least 20% of its maximum population.
  • Optionally, the total upstream withdrawal can be capped by a limit $L$.

2. Mathematical Formulation

Let $P_k$ be the population, $W^{\mathrm{dom}}_k$ the domestic withdrawal, and $W^{\mathrm{agr}}_k$ the irrigation withdrawal of settlement $k$, with total withdrawal $W_k = W^{\mathrm{dom}}_k + W^{\mathrm{agr}}_k$.

$$
\max \sum_{k=1}^{6} P_k
$$

subject to

$$
d,P_k \le W^{\mathrm{dom}}_k, \qquad a_k P_k \le W^{\mathrm{agr}}_k \qquad (k=1,\dots,6)
$$

$$
\sum_{j<k} (1-r),W_j + W_k \le Q - E \qquad (k=1,\dots,6)
$$

$$
\sum_{k \in \text{upstream}} W_k \le L
$$

$$
0.2,P_k^{\max} \le P_k \le P_k^{\max}, \qquad W^{\mathrm{dom}}_k,; W^{\mathrm{agr}}_k \ge 0
$$

The river-balance constraint says that at settlement $k$, the water left in the river after all upstream consumption must still cover this village’s withdrawal while preserving the environmental flow. Because every constraint is linear, the whole problem is a linear program that can be solved exactly.

3. Complete Python Code

The code is a single script. For speed, the constraint matrix is built only once, and only the right-hand side vector changes between the 960 linear programs solved for the sweep, using the fast HiGHS solver.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from scipy.optimize import linprog

# ------------------------------------------------------------
# 1. Problem data
# ------------------------------------------------------------
names = ["Up-1 Highland", "Up-2 Terrace", "Up-3 Forest",
"Down-1 Plain", "Down-2 Delta", "Down-3 Coast"]
n = len(names) # number of settlements
n_up = 3 # first three are upstream

P_max = np.array([20.0, 15.0, 10.0, 25.0, 20.0, 15.0]) # capacity (thousand people)
P_min = 0.2 * P_max # minimum guarantee (20 %)
d = 0.2 # domestic water per capita
a = np.array([1.2, 1.4, 1.6, 1.5, 1.8, 2.0]) # irrigation water per capita
r = 0.3 # return-flow ratio
E = 10.0 # environmental minimum flow

# ------------------------------------------------------------
# 2. Linear program (matrix built once, only the right-hand side changes)
# x = [P_1..P_6, W_dom_1..W_dom_6, W_agr_1..W_agr_6]
# ------------------------------------------------------------
def build_constraint_matrix():
A = np.zeros((3 * n + 1, 3 * n))
row = 0
for k in range(n): # domestic: d * P_k <= W_dom_k
A[row, k] = d
A[row, n + k] = -1.0
row += 1
for k in range(n): # food: a_k * P_k <= W_agr_k
A[row, k] = a[k]
A[row, 2 * n + k] = -1.0
row += 1
for k in range(n): # river balance at settlement k
for j in range(k):
A[row, n + j] = 1.0 - r
A[row, 2 * n + j] = 1.0 - r
A[row, n + k] = 1.0
A[row, 2 * n + k] = 1.0
row += 1
for k in range(n_up): # total upstream withdrawal <= L
A[row, n + k] = 1.0
A[row, 2 * n + k] = 1.0
return A

A_UB = build_constraint_matrix()
C = np.concatenate([-np.ones(n), np.zeros(2 * n)])
BOUNDS = [(P_min[k], P_max[k]) for k in range(n)] + [(0, None)] * (2 * n)

def solve(Q, L):
b = np.zeros(A_UB.shape[0])
b[2 * n:3 * n] = Q - E
b[3 * n] = L
return linprog(C, A_ub=A_UB, b_ub=b, bounds=BOUNDS, method="highs")

# ------------------------------------------------------------
# 3. Base case: Q = 100, no upstream limit
# ------------------------------------------------------------
Q0, L0 = 100.0, 1e6
res = solve(Q0, L0)
x = res.x
P, W_dom, W_agr = x[:n], x[n:2 * n], x[2 * n:]
W_tot = W_dom + W_agr

df = pd.DataFrame({
"Population (thousand)": P,
"Capacity (thousand)": P_max,
"Domestic water": W_dom,
"Irrigation water": W_agr,
"Total withdrawal": W_tot,
}, index=names).round(2)

print("=== Base case (Q = 100, no upstream limit) ===")
print(df.to_string())
print(f"\nTotal population : {P.sum():.2f} thousand")
print(f"Upstream : {P[:n_up].sum():.2f} thousand")
print(f"Downstream : {P[n_up:].sum():.2f} thousand")

# river flow profile
flow = [Q0]
for k in range(n):
flow.append(flow[-1] - (1 - r) * W_tot[k])
flow = np.array(flow)
print("\nRiver flow after each settlement:", np.round(flow, 2))

# ------------------------------------------------------------
# 4. Sweep: river flow Q x upstream withdrawal limit L
# ------------------------------------------------------------
L_min = sum((d + a[k]) * P_min[k] for k in range(n_up))
Q_grid = np.linspace(60, 160, 30)
L_grid = np.linspace(L_min, 100, 30)
Z = np.full((len(L_grid), len(Q_grid)), np.nan)

t0 = time.time()
for i, L in enumerate(L_grid):
for j, Q in enumerate(Q_grid):
rr = solve(Q, L)
if rr.status == 0:
Z[i, j] = -rr.fun
elapsed = time.time() - t0
print(f"\nSolved {Z.size} linear programs in {elapsed:.2f} s")

# trade-off curve at Q = 100
L_curve = np.linspace(L_min, 100, 60)
up_pop, down_pop, tot_pop = [], [], []
for L in L_curve:
rr = solve(Q0, L)
up_pop.append(rr.x[:n_up].sum())
down_pop.append(rr.x[n_up:n].sum())
tot_pop.append(-rr.fun)
up_pop, down_pop, tot_pop = map(np.array, (up_pop, down_pop, tot_pop))

best = int(np.argmax(tot_pop))
print(f"Smallest feasible upstream limit : {L_min:.1f}")
print(f"Best total population : {tot_pop[best]:.2f} thousand "
f"(upstream limit >= {L_curve[best]:.1f})")
print(f"Population at strictest limit : {tot_pop[0]:.2f} thousand")

# ------------------------------------------------------------
# 5. One combined figure
# ------------------------------------------------------------
fig = plt.figure(figsize=(21, 12))
gs = GridSpec(2, 3, figure=fig, hspace=0.35, wspace=0.28)
colors = ["#2a9d8f"] * n_up + ["#e76f51"] * (n - n_up)

# (a) river flow profile
ax1 = fig.add_subplot(gs[0, 0])
labels = ["Source"] + [s.split()[0] for s in names]
ax1.step(range(n + 1), flow, where="post", color="#264653", lw=2.5)
ax1.scatter(range(n + 1), flow, color="#264653", zorder=3)
ax1.axhline(E, color="red", ls="--", label="Environmental minimum flow")
ax1.axvspan(-0.3, n_up + 0.5, color="#2a9d8f", alpha=0.12, label="Upstream")
ax1.axvspan(n_up + 0.5, n + 0.3, color="#e76f51", alpha=0.12, label="Downstream")
ax1.set_xticks(range(n + 1))
ax1.set_xticklabels(labels, rotation=30)
ax1.set_ylabel("River flow")
ax1.set_title("(a) River flow profile")
ax1.legend(loc="lower left")
ax1.grid(alpha=0.3)

# (b) population vs capacity
ax2 = fig.add_subplot(gs[0, 1])
idx = np.arange(n)
ax2.bar(idx - 0.2, P_max, width=0.4, color="lightgray", label="Capacity")
ax2.bar(idx + 0.2, P, width=0.4, color=colors, label="Optimal population")
ax2.set_xticks(idx)
ax2.set_xticklabels([s.split()[0] for s in names], rotation=30)
ax2.set_ylabel("Population (thousand)")
ax2.set_title("(b) Optimal population per settlement")
ax2.legend()
ax2.grid(alpha=0.3, axis="y")

# (c) water allocation
ax3 = fig.add_subplot(gs[0, 2])
ax3.bar(idx, W_dom, color="#457b9d", label="Domestic")
ax3.bar(idx, W_agr, bottom=W_dom, color="#a8dadc", label="Irrigation")
ax3.set_xticks(idx)
ax3.set_xticklabels([s.split()[0] for s in names], rotation=30)
ax3.set_ylabel("Water withdrawal")
ax3.set_title("(c) Water allocation")
ax3.legend()
ax3.grid(alpha=0.3, axis="y")

# (d) 3D surface
ax4 = fig.add_subplot(gs[1, 0:2], projection="3d")
QQ, LL = np.meshgrid(Q_grid, L_grid)
surf = ax4.plot_surface(QQ, LL, Z, cmap="viridis", edgecolor="none", alpha=0.95)
ax4.set_xlabel("River flow Q")
ax4.set_ylabel("Upstream limit L")
ax4.set_zlabel("Max total population (thousand)")
ax4.set_title("(d) Maximum population vs river flow and upstream limit")
ax4.view_init(elev=28, azim=-125)
fig.colorbar(surf, ax=ax4, shrink=0.6, pad=0.08)

# (e) trade-off curve
ax5 = fig.add_subplot(gs[1, 2])
ax5.plot(L_curve, up_pop, color="#2a9d8f", lw=2.5, label="Upstream")
ax5.plot(L_curve, down_pop, color="#e76f51", lw=2.5, label="Downstream")
ax5.plot(L_curve, tot_pop, color="black", lw=2.5, ls="--", label="Total")
ax5.set_xlabel("Upstream withdrawal limit L")
ax5.set_ylabel("Population (thousand)")
ax5.set_title("(e) Upstream-downstream trade-off (Q = 100)")
ax5.legend()
ax5.grid(alpha=0.3)

plt.show()

4. Code Walkthrough

Section 1: Problem data

The arrays P_max, P_min, and a hold each village’s capacity, guaranteed minimum population, and irrigation water requirement per person. The scalars d, r, and E are the domestic need, return-flow ratio, and environmental flow. The first three villages are upstream, which is encoded by n_up = 3.

Section 2: Building the linear program

The decision vector stacks three blocks of six variables: populations, domestic withdrawals, and irrigation withdrawals, 18 variables in total. The constraint matrix A_UB has $3n+1 = 19$ rows:

  • Rows 1 to 6 encode $d P_k - W^{\mathrm{dom}}_k \le 0$.
  • Rows 7 to 12 encode $a_k P_k - W^{\mathrm{agr}}_k \le 0$.
  • Rows 13 to 18 are the river-balance constraints. For settlement $k$, all upstream withdrawals are multiplied by $(1-r)$ because part of that water flows back, while the settlement’s own withdrawal counts fully.
  • Row 19 is the upstream cap $L$.

Since linprog minimizes, the objective vector C is $-1$ for every population variable and $0$ elsewhere, which maximizes total population. The helper solve(Q, L) only rewrites the right-hand side vector b, so the expensive matrix construction happens once. This is the speed-up that makes the 900-point sweep finish in a few seconds.

Section 3: Base case

We solve the problem with $Q=100$ and an effectively unlimited upstream cap (L0 = 1e6). The result is unpacked into populations and withdrawals, shown in a table, and used to compute the river flow after each village with

$$
Q_{k} = Q_{k-1} - (1-r),W_k .
$$

Section 4: Parameter sweep

L_min is the smallest upstream cap that still satisfies the 20% minimum guarantee for upstream villages. We then solve the LP on a $30 \times 30$ grid of river flow $Q \in [60,160]$ and cap $L \in [L_{\min},100]$, storing the optimal population in Z. Infeasible cases would remain NaN, which keeps the code safe. A second loop traces the trade-off curve at $Q=100$.

Section 5: Visualization

All five charts are drawn into one figure using GridSpec. Panels (a) to (c) describe the base case, panel (d) is the 3D surface, and panel (e) is the trade-off curve between upstream and downstream population.

5. Where to Paste the Results

==================================================================
Budget: $100k | Total expected annual loss without controls: $1750.0k
==================================================================

[Optimal allocation]
  Firewall               $  14.40k  ( 14.4%)
  IDS/IPS                $  20.26k  ( 20.3%)
  Employee Training      $  27.17k  ( 27.2%)
  Penetration Testing    $   7.55k  (  7.5%)
  EDR                    $  15.62k  ( 15.6%)
  Patch Management       $  15.00k  ( 15.0%)
  Residual expected loss : $770.35k
  Risk reduction         : 55.98%

[Strategy comparison]
  Equal split              residual $  791.65k | reduction  54.76%
  Risk-weighted split      residual $  782.48k | reduction  55.29%
  Random search (50,000)   residual $  770.78k | reduction  55.96%
  Greedy ($1k steps)       residual $  770.38k | reduction  55.98%
  SLSQP optimum            residual $  770.35k | reduction  55.98%

[Marginal loss reduction per extra $1k at the optimum]
  Firewall                 4.6685
  IDS/IPS                  4.6685
  Employee Training        4.6685
  Penetration Testing      4.6685
  EDR                      4.6685
  Patch Management         4.6685

[Budget sweep]
  Budget $ 20k -> residual $ 1404.18k ( 19.8% reduction)
  Budget $ 50k -> residual $ 1085.86k ( 38.0% reduction)
  Budget $ 80k -> residual $  874.84k ( 50.0% reduction)
  Budget $110k -> residual $  726.10k ( 58.5% reduction)
  Budget $140k -> residual $  618.23k ( 64.7% reduction)
  Budget $170k -> residual $  538.56k ( 69.2% reduction)
  Budget $200k -> residual $  478.78k ( 72.6% reduction)

/tmp/ipykernel_962/2916647809.py:224: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all Axes decorations.
  fig.tight_layout(rect=[0, 0, 1, 0.96])

6. Interpreting the Results

Base case: who gets the water?

The optimal plan supports a total of about 76.2 thousand people: roughly 43.8 thousand upstream and 32.4 thousand downstream. Three observations stand out.

  1. Upstream villages are filled first. Highland and Terrace reach their full capacity, and Forest reaches about 8.8 thousand of its 10 thousand. This is not favoritism. Because 30% of upstream water returns to the river, one unit withdrawn upstream costs the system only $(1-r)=0.7$ units of consumption, while the village’s need per person is among the lowest. Using water where it is cheapest to turn into population maximizes the total.
  2. Down-1 Plain also fills completely, thanks to its relatively low irrigation need of 1.5 per person, while Delta and Coast, with needs of 1.8 and 2.0, end up near the lower end of their allowed range (about 4.4 and 3.0 thousand). The 20% guarantee is what keeps Coast at 3.0 thousand rather than zero.
  3. The river is nearly exhausted at the mouth. Panel (a) shows the flow dropping step by step from 100 to about 12, just above the environmental minimum of 10. Water is the binding resource, and the plan uses it almost completely.

Panels (b) and (c)

Panel (b) compares the optimal population with each village’s capacity. Panel (c) shows that irrigation takes the overwhelming share of withdrawals in every village, roughly six to eight times the domestic use. This is typical of real river basins, where agriculture dominates water demand, and it explains why differences in irrigation efficiency drive the optimal allocation far more than differences in drinking-water needs.

3D surface: panel (d)

The surface shows the maximum total population as a function of river flow $Q$ and the upstream withdrawal limit $L$.

  • Along the $Q$ axis, population rises almost linearly at first, because every additional unit of water directly supports more people, and then flattens out as villages hit their capacity. The maximum on the surface is about 105 thousand, which is the sum of all capacities.
  • Along the $L$ axis, a very strict upstream limit lowers the maximum population, because water that could be used efficiently upstream is forced to flow down to less efficient villages.
  • Beyond a certain limit, increasing $L$ no longer changes anything, which shows up as a flat plateau. At that point the upstream limit is no longer binding.

Trade-off curve: panel (e)

At $Q=100$, the total population is about 64.2 thousand under the strictest limit and climbs to about 76.2 thousand once the limit reaches roughly 68. Upstream population grows with $L$ while downstream population shrinks, and the total keeps rising until the upstream limit stops being binding. This illustrates a useful policy insight: a moderate upstream cap costs the basin very little in total population, but it shifts a meaningful number of people from the upstream to the downstream side. Policy makers who care about fairness can therefore choose a cap in the range where the total curve is still nearly flat, trading a small loss in total population for a much more balanced distribution.

7. Conclusion

By expressing the upstream-downstream water conflict as a linear program, we obtained a provably optimal allocation of domestic and agricultural water, an exact measure of how much each unit of upstream limitation costs in total population, and a clear picture of how river flow and return-flow efficiency shape the maximum sustainable population. The same framework scales easily: adding more villages, seasonal flows, or different return-flow ratios only requires changing the data arrays and the constraint matrix.

Optimizing Fraud Detection Thresholds

Balancing False Positives and False Negatives

Every fraud detection system eventually hits the same question: at what score should we raise an alert? Set the threshold too low and investigators drown in false alarms. Set it too high and real fraud slips through. In this post we build a synthetic payment-fraud problem, find the cost-minimizing threshold, check it against Bayes decision theory, and visualize how the optimum shifts as the business context changes.

1. The Problem Setup

We score 200,000 transactions with a model that outputs a fraud score $s \in [0,1]$. The setup is:

  • Fraud prevalence: $\pi = 0.01$ (1% of transactions)
  • Cost of a false positive (manual review of a legitimate transaction): $C_{FP} = 5$ USD
  • Cost of a false negative (undetected fraud, average loss): $C_{FN} = 300$ USD
  • Legitimate scores follow $\text{Beta}(2, 8)$, and fraudulent scores follow $\text{Beta}(6, 3)$

A transaction raises an alert if $s \ge t$. For threshold $t$, the error rates are

$$
\text{FPR}(t) = \frac{FP(t)}{N_{0}}, \qquad \text{FNR}(t) = \frac{FN(t)}{N_{1}}
$$

where $N_0$ and $N_1$ are the numbers of legitimate and fraudulent transactions. The total operating cost is

$$
\text{Cost}(t) = C_{FP}\cdot FP(t) + C_{FN}\cdot FN(t)
$$

and the goal is

$$
t^{*} = \arg\min_{t \in [0,1]} \text{Cost}(t)
$$

A Theoretical Check

Bayes decision theory says to raise an alert whenever the expected cost of alerting is lower than the expected cost of staying silent:

$$
(1-p(s)),C_{FP} \le p(s),C_{FN} \quad\Longleftrightarrow\quad p(s) \ge \frac{C_{FP}}{C_{FP}+C_{FN}}
$$

where the posterior fraud probability is

$$
p(s) = \frac{\pi f_1(s)}{\pi f_1(s) + (1-\pi) f_0(s)}
$$

Here $f_1$ and $f_0$ are the score densities of fraudulent and legitimate transactions. We will compare the empirical optimum with this theoretical one.

Sensitivity to the Cost Ratio

Define $r = C_{FN}/C_{FP}$. The cost in units of $C_{FP}$ becomes

$$
\frac{\text{Cost}(t)}{C_{FP}} = FP(t) + r \cdot FN(t)
$$

We will explore this as a surface over $(t, r)$ in 3D.

2. The Full Source Code

The naive approach re-evaluates every transaction for every threshold, which costs $O(N \cdot T)$. The fast approach sorts the scores once and uses binary search, which costs $O(N \log N + T \log N)$. Both are included and compared.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
import time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # noqa: F401
from scipy.stats import beta as beta_dist

# ------------------------------------------------------------
# 1. Synthetic data
# ------------------------------------------------------------
rng = np.random.default_rng(42)
N = 200_000
PI = 0.01
C_FP = 5.0
C_FN = 300.0

is_fraud = rng.random(N) < PI
scores = np.where(is_fraud, rng.beta(6, 3, N), rng.beta(2, 8, N))

n1 = int(is_fraud.sum())
n0 = N - n1

thresholds = np.linspace(0.0, 1.0, 1001)

# ------------------------------------------------------------
# 2. Naive evaluation (slow)
# ------------------------------------------------------------
def naive_counts(scores, is_fraud, thresholds):
fp = np.zeros(len(thresholds), dtype=np.int64)
fn = np.zeros(len(thresholds), dtype=np.int64)
for i, t in enumerate(thresholds):
alert = scores >= t
fp[i] = np.sum(alert & ~is_fraud)
fn[i] = np.sum(~alert & is_fraud)
return fp, fn

# ------------------------------------------------------------
# 3. Fast evaluation (sort once + binary search)
# ------------------------------------------------------------
s0 = np.sort(scores[~is_fraud])
s1 = np.sort(scores[is_fraud])

def fast_counts(th):
fp = len(s0) - np.searchsorted(s0, th, side="left")
fn = np.searchsorted(s1, th, side="left")
return fp.astype(np.int64), fn.astype(np.int64)

t0 = time.perf_counter()
fp_naive, fn_naive = naive_counts(scores, is_fraud, thresholds)
time_naive = time.perf_counter() - t0

t0 = time.perf_counter()
fp, fn = fast_counts(thresholds)
time_fast = max(time.perf_counter() - t0, 1e-9)

assert np.array_equal(fp, fp_naive) and np.array_equal(fn, fn_naive)

# ------------------------------------------------------------
# 4. Cost minimization
# ------------------------------------------------------------
fpr = fp / n0
fnr = fn / n1
cost = C_FP * fp + C_FN * fn

best = int(np.argmin(cost))
t_best = thresholds[best]
t_default = 0.5
default = int(np.argmin(np.abs(thresholds - t_default)))

# Theoretical (Bayes) threshold
s_grid = np.linspace(0.001, 0.999, 9999)
f1 = beta_dist.pdf(s_grid, 6, 3)
f0 = beta_dist.pdf(s_grid, 2, 8)
posterior = PI * f1 / (PI * f1 + (1 - PI) * f0)
c_star = C_FP / (C_FP + C_FN)
t_bayes = s_grid[np.argmax(posterior >= c_star)]

# ------------------------------------------------------------
# 5. Console report
# ------------------------------------------------------------
print("=" * 62)
print("Fraud Detection Threshold Optimization")
print("=" * 62)
print(f"Transactions: {N:,} (fraud: {n1:,} / legitimate: {n0:,})")
print(f"C_FP = {C_FP:.0f} USD, C_FN = {C_FN:.0f} USD")
print("-" * 62)
print(f"{'':24s}{'t = 0.50':>16s}{'optimal':>16s}")
print(f"{'Threshold':24s}{thresholds[default]:>16.3f}{t_best:>16.3f}")
print(f"{'False positives':24s}{fp[default]:>16,d}{fp[best]:>16,d}")
print(f"{'False negatives':24s}{fn[default]:>16,d}{fn[best]:>16,d}")
print(f"{'FPR':24s}{fpr[default]:>16.4%}{fpr[best]:>16.4%}")
print(f"{'FNR':24s}{fnr[default]:>16.4%}{fnr[best]:>16.4%}")
print(f"{'Total cost (USD)':24s}{cost[default]:>16,.0f}{cost[best]:>16,.0f}")
saving = cost[default] - cost[best]
print("-" * 62)
print(f"Saving vs. t = 0.50: {saving:,.0f} USD ({saving / cost[default]:.1%})")
print(f"Empirical optimal threshold : {t_best:.3f}")
print(f"Bayes theoretical threshold : {t_bayes:.3f}")
print("-" * 62)
print(f"Naive loop : {time_naive:.4f} sec")
print(f"Fast method: {time_fast:.6f} sec")
print(f"Speedup : {time_naive / time_fast:,.0f}x")

# ------------------------------------------------------------
# 6. Sensitivity analysis over the cost ratio r = C_FN / C_FP
# ------------------------------------------------------------
ratios_line = np.logspace(0, np.log10(300), 80)
cost_mat = fp[None, :] + ratios_line[:, None] * fn[None, :]
t_opt_emp = thresholds[np.argmin(cost_mat, axis=1)]

t_opt_bayes = np.zeros(len(ratios_line))
for i, r in enumerate(ratios_line):
c = 1.0 / (1.0 + r)
t_opt_bayes[i] = s_grid[np.argmax(posterior >= c)]

# 3D surface data
thr_3d = np.linspace(0.2, 0.95, 76)
ratios_3d = np.linspace(5, 200, 60)
fp3, fn3 = fast_counts(thr_3d)
Tm, Rm = np.meshgrid(thr_3d, ratios_3d)
Z = (fp3[None, :] + ratios_3d[:, None] * fn3[None, :]) * C_FP / N * 1000.0
Z = np.minimum(Z, 1500.0)
best_idx = np.argmin(Z, axis=1)
line_t = thr_3d[best_idx]
line_z = Z[np.arange(len(ratios_3d)), best_idx]

# ------------------------------------------------------------
# 7. All plots in one figure
# ------------------------------------------------------------
fig = plt.figure(figsize=(20, 12))

ax1 = fig.add_subplot(2, 3, 1)
ax1.hist(scores[~is_fraud], bins=60, alpha=0.6, density=True, label="Legitimate")
ax1.hist(scores[is_fraud], bins=60, alpha=0.6, density=True, label="Fraud")
ax1.axvline(t_best, color="red", ls="--", label=f"Optimal t = {t_best:.2f}")
ax1.set_title("Score Distributions")
ax1.set_xlabel("Fraud score")
ax1.set_ylabel("Density")
ax1.legend()

ax2 = fig.add_subplot(2, 3, 2)
ax2.plot(thresholds, fpr, label="FPR")
ax2.plot(thresholds, fnr, label="FNR")
ax2.axvline(t_best, color="red", ls="--", label="Optimal t")
ax2.set_title("FPR and FNR vs. Threshold")
ax2.set_xlabel("Threshold t")
ax2.set_ylabel("Rate")
ax2.legend()

ax3 = fig.add_subplot(2, 3, 3)
ax3.plot(thresholds, C_FP * fp, label="False positive cost")
ax3.plot(thresholds, C_FN * fn, label="False negative cost")
ax3.plot(thresholds, cost, color="black", lw=2, label="Total cost")
ax3.scatter([t_best], [cost[best]], color="red", zorder=5, s=80, label="Minimum")
ax3.set_title("Cost Curve")
ax3.set_xlabel("Threshold t")
ax3.set_ylabel("Cost (USD)")
ax3.legend()

ax4 = fig.add_subplot(2, 3, 4)
ax4.plot(fpr, 1 - fnr, lw=2)
ax4.plot([0, 1], [0, 1], "k:", alpha=0.5)
ax4.scatter([fpr[best]], [1 - fnr[best]], color="red", s=90, zorder=5, label="Optimal")
ax4.scatter([fpr[default]], [1 - fnr[default]], color="orange", s=90, zorder=5, label="t = 0.50")
ax4.set_title("ROC Curve with Operating Points")
ax4.set_xlabel("False positive rate")
ax4.set_ylabel("True positive rate")
ax4.legend()

ax5 = fig.add_subplot(2, 3, 5, projection="3d")
ax5.plot_surface(Tm, Rm, Z, cmap="viridis", alpha=0.85, linewidth=0)
ax5.plot(line_t, ratios_3d, line_z, color="red", lw=3, label="Optimal threshold path")
ax5.set_title("Cost Surface: Threshold x Cost Ratio")
ax5.set_xlabel("Threshold t")
ax5.set_ylabel("Cost ratio r")
ax5.set_zlabel("Cost per 1,000 txns (USD)")
ax5.view_init(elev=28, azim=-125)
ax5.legend()

ax6 = fig.add_subplot(2, 3, 6)
ax6.plot(ratios_line, t_opt_emp, lw=2, label="Empirical optimum")
ax6.plot(ratios_line, t_opt_bayes, "--", lw=2, label="Bayes theoretical")
ax6.axvline(C_FN / C_FP, color="gray", ls=":", label="Our case (r = 60)")
ax6.set_xscale("log")
ax6.set_title("Optimal Threshold vs. Cost Ratio")
ax6.set_xlabel("Cost ratio r = C_FN / C_FP")
ax6.set_ylabel("Optimal threshold")
ax6.legend()

plt.tight_layout()
plt.show()

3. Code Walkthrough

Data generation

We draw a Bernoulli label with probability $\pi = 0.01$ for each transaction, then draw a score from $\text{Beta}(6,3)$ for fraud and $\text{Beta}(2,8)$ for legitimate transactions. The two distributions overlap, which is exactly why a trade-off exists. If they were perfectly separable, a single threshold would produce zero errors.

Naive evaluation

naive_counts loops over 1,001 thresholds. In each iteration it builds a Boolean alert vector across all 200,000 transactions and counts false positives (alert & ~is_fraud) and false negatives (~alert & is_fraud). It is easy to read, but it touches every transaction for every threshold.

Fast evaluation

fast_counts sorts the legitimate and fraudulent scores once. For any threshold $t$:

  • The number of legitimate transactions with $s \ge t$ is $N_0$ minus the count below $t$, which is exactly len(s0) - searchsorted(s0, t).
  • The number of fraudulent transactions with $s < t$ is searchsorted(s1, t).

Using side="left" makes the boundary behave identically to the naive scores >= t rule, so the assert confirms the two methods agree exactly. The sorting cost is paid once, and each threshold then costs only $O(\log N)$.

Cost minimization

With fp and fn arrays in hand, cost = C_FP * fp + C_FN * fn is a single vectorized expression, and np.argmin returns the best threshold. We also report the result at the conventional default of $t = 0.5$ to quantify the improvement.

Bayes threshold

We compute the posterior $p(s)$ on a fine grid using the exact Beta densities, then find the first score where $p(s) \ge C_{FP}/(C_{FP}+C_{FN})$. Because the likelihood ratio $f_1/f_0$ is monotonically increasing in $s$, the posterior crosses the cutoff exactly once, so this is the theoretical optimum. The empirical minimum should land very close to it.

Sensitivity analysis

The 2D analysis uses 80 cost ratios on a log scale from 1 to 300 and reuses the same fp and fn arrays, so the entire ratio-by-threshold matrix is built with broadcasting and no extra passes over the data. The 3D surface uses a restricted threshold range of $[0.2, 0.95]$ and caps the cost at 1,500 USD per 1,000 transactions so that extreme values at very low thresholds do not flatten the interesting region.

4. Execution Results

Console Output

==============================================================
Fraud Detection Threshold Optimization
==============================================================
Transactions: 200,000  (fraud: 2,065 / legitimate: 197,935)
C_FP = 5 USD, C_FN = 300 USD
--------------------------------------------------------------
                                t = 0.50         optimal
Threshold                          0.500           0.450
False positives                    3,971           7,784
False negatives                      317             188
FPR                              2.0062%         3.9326%
FNR                             15.3511%         9.1041%
Total cost (USD)                 114,955          95,320
--------------------------------------------------------------
Saving vs. t = 0.50: 19,635 USD (17.1%)
Empirical optimal threshold : 0.450
Bayes theoretical threshold : 0.442
--------------------------------------------------------------
Naive loop : 1.4651 sec
Fast method: 0.000544 sec
Speedup    : 2,692x

Graph Output

5. Interpreting the Graphs

Score Distributions (top left). The legitimate scores pile up near the low end, while the fraud scores concentrate at higher values. The overlap in the middle is the region where every threshold choice necessarily trades one error type for the other. The red dashed line marks the cost-optimal threshold, which sits well to the left of where the two histograms would cross if the classes were equally likely, because a missed fraud is far more expensive than a wasted review.

FPR and FNR vs. Threshold (top center). As $t$ increases, FPR falls and FNR rises. This is the core trade-off. Raising the threshold reduces investigator workload but lets more fraud through.

Cost Curve (top right). The false-positive cost decreases with $t$, the false-negative cost increases, and their sum forms a U-shaped curve. The minimum is the optimal operating point. Notice how flat the curve is near the bottom. Small deviations from $t^*$ cost little, but moving far away, especially toward high thresholds, is expensive because every missed fraud costs 60 times as much as a review.

ROC Curve (bottom left). The optimal operating point is not the point closest to the top-left corner, nor the point that balances FPR and FNR. It sits where the slope of the ROC curve equals the cost-and-prevalence-weighted ratio $\dfrac{(1-\pi),C_{FP}}{\pi,C_{FN}}$. In practice this puts the operating point at a noticeably higher true positive rate than the default $t = 0.5$, accepting more false alarms in exchange for catching more fraud.

3D Cost Surface (bottom center). The two horizontal axes are the threshold $t$ and the cost ratio $r$, and the height is the cost per 1,000 transactions. The surface forms a valley whose floor shifts as $r$ grows. The red line traces the valley floor, the optimal threshold for each ratio. As fraud becomes relatively more expensive, the valley moves toward lower thresholds, and the walls on the high-threshold side rise steeply. This view makes it clear that the optimal threshold is not a universal constant. It is a function of the business cost structure.

Optimal Threshold vs. Cost Ratio (bottom right). The empirical optimum, found by brute-force search over the data, follows the Bayes-theoretical curve closely across the whole range. The small deviations come from sampling noise in the 2,000 or so fraud cases. The gray dotted line marks our scenario of $r = 60$. As $r$ grows, the optimal threshold decreases steadily, which means a business facing larger fraud losses should be willing to review more alerts.

6. Conclusion

Threshold selection in fraud detection is a cost-minimization problem, not an accuracy-maximization problem. With heavily imbalanced classes and asymmetric costs, the default threshold of 0.5 is almost never optimal. Three lessons stand out:

  1. Define the costs explicitly. The optimal threshold depends on $\pi$, $C_{FP}$, and $C_{FN}$, and changes whenever any of them change.
  2. Use the theory as a sanity check. The Bayes rule $p(s) \ge C_{FP}/(C_{FP}+C_{FN})$ gives a principled target that empirical search should match.
  3. Vectorize the search. Sorting once and using binary search makes evaluating thousands of candidate thresholds essentially free, which makes sensitivity analysis like the 3D surface above practical even on large datasets.

In production, you can extend this framework with per-transaction loss amounts, a limit on the number of alerts analysts can review per day, or threshold schedules by customer segment. The same cost-minimization logic carries over directly.

Managing Renewable and Finite Resources Over 100 Turns with Python

Planning Long-Term Civilization Survival

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

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

1. The Example Problem

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

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

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

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

Three random shocks hit each scenario:

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

2. Mathematical Model

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

Renewable resource (logistic regrowth minus harvest):

Finite resource:

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

Production (labor and resources, with diminishing returns):

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

Demand and supply ratio:

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

Population growth (decline is faster than growth):

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

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

Efficiency growth:

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

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

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

3. Solution Strategy

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

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

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

4. Source Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from scipy.optimize import differential_evolution

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5. Code Walkthrough

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

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

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

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

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

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

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

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

6. Execution Results

Console Output

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

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

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

Graph Output

7. Interpreting the Results

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

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

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

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

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

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

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

8. Conclusion

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

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

Optimizing Equipment Investment to Maximize Cumulative Population Over 20 Turns

Automation vs. Scarcity

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

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


1. The Problem

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

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

Model

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

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

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

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

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

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

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

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

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

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

Parameters

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

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


2. Source Code

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
import time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # noqa: F401
from scipy.optimize import differential_evolution, minimize

plt.style.use("dark_background")

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

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


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

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

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

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

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

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

return total, hist


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


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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

plt.show()

3. Console Output

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

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

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

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

4. Result Figure


5. Code Walkthrough

5.1 The vectorized simulator

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

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

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

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

Two details matter for speed and reliability:

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

5.3 Baseline strategies

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

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

5.4 The 3D slice

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


6. Reading the Graphs

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

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

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

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

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

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

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

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

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


7. Takeaways

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

Optimizing Land Use for Population and Environmental Value Over 50 Years

Forest or Farmland?

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

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

1. The Problem Setting

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

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

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

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

The conversion in year $t$ is

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

and the land-use dynamics are

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

2. The Mathematical Model

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

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

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

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

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

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

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

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

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

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

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

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

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

Parameter values used in the code:

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

3. Solution Strategy

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

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

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

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

4. Source Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

plt.style.use("dark_background")

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5. Code Walkthrough

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

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

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

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

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

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

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

6. Execution Results

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

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

7. Interpreting the Results

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

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

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

The comparison at year 50 makes this concrete:

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

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

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

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

8. Summary

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

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

Maximizing 20-Turn Cumulative Population with Python

How Much Should a Region Spend on Disaster Preparedness?

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

1. Problem Setup

A region starts with the following state.

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

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

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

A fixed share of output becomes the budget:

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

The budget is split into four uses.

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

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

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

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

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

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

Combined, the fraction of assets lost is

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

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

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

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

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

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

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

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

2. Source Code

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from mpl_toolkits.mplot3d import Axes3D # noqa: F401
from scipy.optimize import differential_evolution

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3. Code Walkthrough

Parameters (Section 1)

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

Scenario generation (Section 2)

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

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

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

The simulator (Section 3)

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

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

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

Parameter decoding (decode)

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

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

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

Optimization (Section 4)

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

Strategy comparison (Section 5)

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

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

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

3D surface scan (Section 6)

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

4. Execution Results

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

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

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

Grid scan time: 1.5 s

5. Reading the Results

Optimal allocation

The optimizer found a clear pattern:

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

Three points stand out.

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

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

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

Strategy comparison

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

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

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

Figure guide

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

6. Conclusion

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

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

Maximize Its Population in 30 Years

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

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

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

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

1. The Model

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

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

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

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

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

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

Policy. Two parameters define the policy:

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

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

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

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

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

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

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

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

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

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

2. Full Source Code

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import time

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

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

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

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

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

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

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

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

drought_prob = (R_test < DROUGHT).mean()

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

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

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

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

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

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

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

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

plt.show()

3. Code Walkthrough

Section 1: Parameters

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

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

Section 2: The Simulator

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

Inside the yearly loop, the order of events matters:

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

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

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

Section 4: Out-of-Sample Evaluation

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

Section 5: Graphs

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

4. Results

Console Output

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

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

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

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

Graphs

5. Reading the Results

The comparison table

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

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

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

Graph (1): Rainfall distribution

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

Graph (2): 3D surface of expected population

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

Graph (3): Heat map

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

Graph (4): Population paths

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

Graph (5): Distribution of final population

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

Graph (6): Irrigation coverage and reserve over time

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

6. Summary

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

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