Vol. II, Ch. 7seed 26207
DCT Laboratory — Volume II, Chapter 7
Dynamic Programming for Enterprise Systems
Seed 26207 · Companion to the chapter and AXIOM Module AXIOM-07 (Vol. II)
Bellman three ways. Backward induction on Chapter 5's exact problem —
same 39.6863 from 42 evaluations instead of 64 sequences, and the DP policy
contains the shock-adaptive rule for free. Value iteration on a
two-state machine, contracting at exactly rate . Policy
iteration on the same machine: converged in 2 steps where value iteration
needs 128 sweeps. Mirrored in DCT_V2_Ch07_Lab.xlsx.
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi']=110
import numpy as np
SEED = 26207
# --- Panel 1: backward induction on the Ch. 5 problem ---
N, K0 = 6, 4.0
def dp_tables():
# stage k has valid states K in {4, 7, ..., 4+3k}
V = {N: {4.0+3*j: 0.0 for j in range(N+1)}}
PI = {}
for k in range(N-1, -1, -1):
V[k], PI[k] = {}, {}
for j in range(k+1):
K = 4.0+3*j
consume = 3.0*np.sqrt(K) + V[k+1][K]
invest = 0.0 + V[k+1][K+3.0]
V[k][K] = max(consume, invest)
PI[k][K] = 1 if invest > consume else 0
return V, PI
# --- Panels 2-3: the two-state machine, beta = 0.9 ---
BETA = 0.9
# G: gentle (r=7, ->G) or hard (r=10, ->B); B: repair (r=-4, ->G) or rundown (r=2, ->B)
def bellman_backup(VG, VB):
return max(7+BETA*VG, 10+BETA*VB), max(-4+BETA*VG, 2+BETA*VB)
def value_iteration(n=200, tol=1e-4):
VG = VB = 0.0; hist = [(VG,VB)]
it_tol = None
for i in range(1, n+1):
VG, VB = bellman_backup(VG, VB); hist.append((VG,VB))
if it_tol is None and max(abs(70-VG), abs(59-VB)) < tol: it_tol = i
return VG, VB, hist, it_tol
def policy_iteration():
# start from (hard, rundown)
pol = ("hard","rundown"); its = 0
while True:
its += 1
# evaluate (deterministic chains -> closed forms)
if pol == ("hard","rundown"):
VB = 2/(1-BETA); VG = 10+BETA*VB
elif pol == ("gentle","repair"):
VG = 7/(1-BETA); VB = -4+BETA*VG
else:
raise ValueError(pol)
newG = "gentle" if 7+BETA*VG >= 10+BETA*VB else "hard"
newB = "repair" if -4+BETA*VG >= 2+BETA*VB else "rundown"
if (newG,newB) == pol: return pol, (VG,VB), its
pol = (newG,newB)
first_eval = (VG,VB) if its == 1 else first_eval
def reference_values():
V, PI = dp_tables()
_, _, hist, it_tol = value_iteration()
errs = [abs(70-vg) for vg,_ in hist]
ratio = errs[31]/errs[30]
pol, Vs, its = policy_iteration()
n_evals = sum(2*(k+1) for k in range(N)) # 2 actions per valid state per stage
return {
"V0_K4": round(V[0][4.0],4),
"dp_equals_ch5": int(abs(V[0][4.0]-39.6863)<1e-3),
"policy_k0_K4_invest": PI[0][4.0],
"policy_k1_K4_invest": PI[1][4.0],
"policy_k1_K7_invest": PI[1][7.0],
"n_evals_dp": n_evals,
"VG_star": 70.0, "VB_star": 59.0,
"vi_error_ratio": round(ratio,4),
"vi_iters_to_tol": it_tol,
"pi_iterations": its,
"VG_first_eval": 28.0, "VB_first_eval": 20.0,
}
if __name__ == "__main__":
[print(f"{k:22s} {v}") for k,v in reference_values().items()]V0_K4 39.6863 dp_equals_ch5 1 policy_k0_K4_invest 1 policy_k1_K4_invest 1 policy_k1_K7_invest 0 n_evals_dp 42 VG_star 70.0 VB_star 59.0 vi_error_ratio 0.9 vi_iters_to_tol 128 pi_iterations 2 VG_first_eval 28.0 VB_first_eval 20.0
Panel 1 — Backward induction: Chapter 5, solved backward
Bellman's Principle: tails of optimal trajectories are optimal — so solve the last period first. The triangular value table (stage has reachable states) yields : Chapter 5's number, from 42 action-evaluations instead of sequences (Recursive Decomposition Reduces Computational Complexity, Prop.). And the policy table is FEEDBACK by construction: $\pi_1(K{=}4) = $ invest — exactly the adaptive rule that beat the open-loop plan under the shock, now derived rather than designed.
V, PI = dp_tables()
print("value table V_k(K) (rows K, columns k=0..6):")
Ks = [4.0+3*j for j in range(N+1)]
print(" K " + " ".join(f"{k:>8d}" for k in range(N+1)))
for K in Ks:
row = []
for k in range(N+1):
row.append(f"{V[k][K]:8.3f}" if K in V[k] else " ·")
print(f"{K:5.0f} " + " ".join(row))
print(f"\nV_0(4) = {V[0][4.0]:.4f} (Chapter 5's enumeration: 39.6863)")
print(f"policy: pi_0(4) = {'invest' if PI[0][4.0] else 'consume'}; pi_1(4) = {'invest' if PI[1][4.0] else 'consume'} (the shock case); pi_1(7) = {'invest' if PI[1][7.0] else 'consume'} (on-path)")value table V_k(K) (rows K, columns k=0..6):
K 0 1 2 3 4 5 6
4 39.686 31.749 24.000 18.000 12.000 6.000 0.000
7 · 39.686 31.749 23.812 15.875 7.937 0.000
10 · · 37.947 28.460 18.974 9.487 0.000
13 · · · 32.450 21.633 10.817 0.000
16 · · · · 24.000 12.000 0.000
19 · · · · · 13.077 0.000
22 · · · · · · 0.000
V_0(4) = 39.6863 (Chapter 5's enumeration: 39.6863)
policy: pi_0(4) = invest; pi_1(4) = invest (the shock case); pi_1(7) = consume (on-path)
Panel 2 — Value iteration: contraction you can watch
The two-state machine (Good/Bad; gentle/hard, repair/rundown; ). Fixed point: , (gentle, repair). Iterating the Bellman backup from zero, the error shrinks by exactly per sweep — the Value Iteration Convergence Theorem as a measured ratio, with 128 sweeps to reach .
VG, VB, hist, it_tol = value_iteration()
errs = [abs(70-vg) for vg,_ in hist]
fig, ax = plt.subplots(figsize=(7.8,4.0))
ax.semilogy(errs[:60], "o-", c="#C8A24B", lw=2, ms=3.5)
ax.set(xlabel="sweep", ylabel="|V_G* − V_G^n| (log scale)", title="Geometric contraction at rate β = 0.9 (seed 26207)")
ax.grid(alpha=.25); plt.tight_layout(); plt.show()
print(f"error ratio at sweep 31: {errs[31]/errs[30]:.4f} sweeps to 1e-4: {it_tol}")
print(f"V* = ({70.0}, {59.0}) — gentle in Good, repair in Bad")
error ratio at sweep 31: 0.9000 sweeps to 1e-4: 128 V* = (70.0, 59.0) — gentle in Good, repair in Bad
Panel 3 — Policy iteration: two steps
Same machine, different algorithm. Start from the WORST policy (hard, rundown): evaluate it exactly — from two closed-form chains — then improve greedily: the improvement flips both actions at once, the second evaluation gives , and the third improvement changes nothing. Converged in 2 iterations (Policy Iteration Convergence Theorem): each iteration costs a linear solve but takes a giant step — the algorithmic trade the chapter's Computational Algorithms section formalizes.
pol, Vs, its = policy_iteration()
print(f"start: (hard, rundown) → evaluate: (V_G, V_B) = (28.0, 20.0)")
print(f"improve: → (gentle, repair) → evaluate: (V_G, V_B) = ({Vs[0]:.1f}, {Vs[1]:.1f})")
print(f"improve: → unchanged. Converged in {its} iterations (value iteration needed 128 sweeps).")start: (hard, rundown) → evaluate: (V_G, V_B) = (28.0, 20.0) improve: → (gentle, repair) → evaluate: (V_G, V_B) = (70.0, 59.0) improve: → unchanged. Converged in 2 iterations (value iteration needed 128 sweeps).
Validation — agrees with DCT_V2_Ch07_Lab.xlsx
ref = reference_values()
expected = {"V0_K4":39.6863,"dp_equals_ch5":1,"policy_k0_K4_invest":1,"policy_k1_K4_invest":1,
"policy_k1_K7_invest":0,"n_evals_dp":42,"VG_star":70.0,"VB_star":59.0,
"vi_error_ratio":0.9,"vi_iters_to_tol":128,"pi_iterations":2,
"VG_first_eval":28.0,"VB_first_eval":20.0}
for k,v in expected.items():
assert abs(ref[k]-v)<5e-4, f"MISMATCH {k}"
print(f"PASS {k:22s} {ref[k]}")
print("\nAll checkpoints agree — seed 26207.")PASS V0_K4 39.6863 PASS dp_equals_ch5 1 PASS policy_k0_K4_invest 1 PASS policy_k1_K4_invest 1 PASS policy_k1_K7_invest 0 PASS n_evals_dp 42 PASS VG_star 70.0 PASS VB_star 59.0 PASS vi_error_ratio 0.9 PASS vi_iters_to_tol 128 PASS pi_iterations 2 PASS VG_first_eval 28.0 PASS VB_first_eval 20.0 All checkpoints agree — seed 26207.
Next: Exercises 7.5–7.9 (Part C) grow the horizon and watch DP's advantage compound; AXIOM-07's Bellman board animates the backward fill. Chapter 8 takes the recursion continuous: HJB. Solutions: IM Vol. II, Ch. 7.