Finite-Horizon Dynamic Programming

MPC creates feedback by solving a new trajectory problem from each measured state. Can the decisions for an entire family of states be computed in advance instead? Dynamic programming decomposes a finite-horizon control problem into cost-to-go functions indexed by time and state.

The deterministic discrete-time optimal control problem provides the starting point. Stochastic successors will enter in the next chapter through conditional expectations.

Consider a typical DOCP of Bolza type:

Rather than considering only the total cost from the initial time to the final time, dynamic programming introduces the concept of cost from an arbitrary point in time to the end. This leads to the definition of the “cost-to-go” or “value function” $J_k(\mathbf{x}_k)$:

This function represents the total cost incurred from stage $k$ onwards to the end of the time horizon, given that the system is initialized in state $\mathbf{x}_k$ at stage $k$. Suppose the problem has been solved from stage $k+1$ to the end, yielding the optimal cost-to-go $J_{k+1}^\star(\mathbf{x}_{k+1})$ for any state $\mathbf{x}_{k+1}$ at stage $k+1$. The question then becomes: how does this information inform the decision at stage $k$?

Given knowledge of the optimal behavior from $k+1$ onwards, the task reduces to determining the optimal action $\mathbf{u}_k$ at stage $k$. This control should minimize the sum of the immediate cost $c_k(\mathbf{x}_k, \mathbf{u}_k)$ and the optimal future cost $J_{k+1}^\star(\mathbf{x}_{k+1})$, where $\mathbf{x}_{k+1}$ is the resulting state after applying action $\mathbf{u}_k$. Mathematically, this is expressed as:

This equation is known as Bellman’s equation, named after Richard Bellman, who formulated the principle of optimality:

An optimal policy has the property that whatever the previous state and decision, the remaining decisions must constitute an optimal policy with regard to the state resulting from the previous decision.

In other words, any sub-path of an optimal path, from any intermediate point to the end, must itself be optimal. This principle is the basis for the backward induction procedure which computes the optimal value function and provides closed-loop control capabilities without having to use an explicit NLP solver.

Dynamic programming can handle nonlinear systems and non-quadratic cost functions naturally. It provides a global optimal solution, when one exists, and can incorporate state and control constraints with relative ease. These advantages come from a change in what is being computed. Trajectory optimization takes the initial state $\mathbf{x}_1$ as given and returns one sequence of states and controls, typically a local optimum. With state dimension $n$ and control dimension $m$, the resulting nonlinear program has on the order of $T(n+m)$ decision variables, so its size grows only linearly with the horizon and the dimension. Parametric optimization widened the question to a family of problems indexed by the initial state, but still answered it locally: sensitivity updates carry the solution of one instance to a nearby one, and learned controllers are fitted only on states sampled from the operating domain. Backward recursion asks for the optimal cost-to-go and the optimal decision from every state at every stage. The answer is a function on the whole state space rather than a vector of decision variables.

Representing such a function is what becomes expensive as the dimension of the state space increases, a difficulty Bellman termed the “curse of dimensionality.” If each of the $n$ state coordinates is discretized into $N$ levels, the table storing $J_t^\star$ has $N^n$ entries, and each entry requires its own minimization over the controls. Computation and memory therefore grow exponentially with the state dimension, rendering direct application of dynamic programming intractable for high-dimensional problems. The curse is the price of asking for a solution everywhere. Trajectory optimization gives up the global, closed-loop answer to keep its cost polynomial in the dimension; dynamic programming obtains that answer and pays for it with exponential growth.

Fortunately, learning-based methods offer efficient tools to combat the curse of dimensionality on two fronts: by using function approximation (e.g., neural networks) to avoid explicit discretization, and by leveraging randomization through Monte Carlo methods inherent in the learning paradigm. Most of this course is dedicated to those ideas.

Backward Recursion¶

How does the terminal cost propagate backward into an optimal action and value for every earlier state?

The principle of optimality provides a methodology for solving optimal control problems. Beginning at the final time horizon and working backwards, at each stage the local optimization problem given by Bellman’s equation is solved. This process, termed backward recursion or backward induction, constructs the optimal value function stage by stage.

Upon completion of this backward pass, we now have access to the optimal control to take at any stage and in any state. Furthermore, we can simulate optimal trajectories from any initial state by applying the optimal policy at each stage to generate the optimal trajectory.

The proof uses a “cut-and-paste” argument. Any control sequence can be cut after its first step and its tail replaced by an optimal continuation without increasing its cost. Comparing first steps, each followed by its optimal continuation, then gives a sequence that no other sequence can beat. Figure 1 follows these two moves with concrete costs.

Example: Seasonal Operation of a Hydroelectric Reservoir¶

On a northern river, most of the water arrives in a few weeks of spring, when the snowpack melts. This is the freshet. Electricity is worth the most in winter, when heating demand peaks. A reservoir behind a dam moves water from one season to the other. Each week the operator decides how much water to send through the turbines. Water released now is sold at this week’s price. Water kept in the reservoir can be sold later, and it also raises the water level, so every cubic metre released afterwards falls from a greater height and produces more energy. Keeping water has a risk as well. When the freshet reaches a full reservoir, the excess goes over the spillway and earns nothing.

The hydro cascade in the trajectory optimization chapter computed one plan from one initial state. An operator needs more than that. Inflows never follow the forecast exactly, so the reservoir is rarely at the level the plan assumed, and the operator needs a release for whatever level it is at. Dynamic programming provides this: one backward pass gives a release for every week and every storage level. Reservoir operators call the resulting curves rule curves, and dynamic programming has been used to compute them for decades Yeh, 1985.

Model. The stage is one week, with $T=52$ weeks starting in January. Volumes are measured in cubic hectometres ($1\,\mathrm{hm}^3 = 10^6\,\mathrm{m}^3$).

The inflows and prices are stylized curves, not measurements: a winter base flow of 100 m³/s, a freshet that peaks in mid-May at about four times the mean flow, an autumn rise from rain, and a price that is highest in early February and lowest in August. Water above capacity spills, so storage evolves as

Releasing $u$ with head $h(x)$ produces $\kappa\,h(x)\,u$ megawatt-hours, where $\kappa = \rho g \eta = 1000 \cdot 9.81 \cdot 0.9$ joules per cubic metre per metre, converted to MWh per hm³. The head grows with storage, which is what makes keeping water valuable beyond its later sale price. Water left at the end of the year is credited at $p_T = 50$ $/MWh, so the terminal value is $J_T^\star(x) = p_T\,\kappa\,h(x)\,(x - x_{\min})$. Since the objective is a revenue to maximize rather than a cost, Bellman’s equation takes a maximum:

Call the expression inside the maximum

the total revenue from storage $x$ in week $t$ if we release $u$ now and act optimally afterwards. It weighs the revenue earned by releasing $u$ this week against the revenue still obtainable from the storage that remains, and Bellman’s equation reads $J_t^\star(x) = \max_u Q_t(x, u)$.

Code. The storage grid has a step of 10 hm³, and the inflows are rounded to multiples of 10 hm³. Every successor $x + w_t - u$ is therefore a grid point, and Algorithm 1 applies as stated: $J_{t+1}^\star$ is looked up, never interpolated. The next section removes this convenience. The loop over states and actions of the algorithm is written with NumPy broadcasting: x is a column of all storage levels and u a row of all releases, so each line evaluates every (storage, release) pair at once.

import numpy as np

# Calendar: 52 weekly decisions; week 0 is the first week of January.
T = 52
weeks = np.arange(T)

# Inflow forecast w_t (hm³ per week): winter base flow, spring freshet, autumn rain.
w = 60 + 340 * np.exp(-((weeks - 19) / 3.5) ** 2) + 60 * np.exp(-((weeks - 43) / 4) ** 2)
w = 10 * np.round(w / 10)            # multiples of the grid step

# Price p_t ($/MWh): winter peak, summer trough.
p = 60 + 30 * np.cos(2 * np.pi * (weeks - 5) / 52)

# Reservoir and plant.
x_min, x_max = 200, 2000             # dead storage and capacity (hm³)
u_min, u_max = 30, 250               # ecological minimum and turbine capacity (hm³/week)
kappa = 2.45                         # MWh per hm³ per metre of head
p_T = 50                             # $/MWh credited to water left at the end of the year

def head(x):                         # metres of head; grows with storage
    return 40 + 40 * np.sqrt(x / x_max)

def terminal_value(x):
    return p_T * kappa * head(x) * (x - x_min)

dx = 10
X = np.arange(x_min, x_max + dx, dx)   # storage grid: 181 levels
U = np.arange(u_min, u_max + dx, dx)   # release grid: 23 levels

def index(x):                          # position of an on-grid storage level
    return np.round((x - x_min) / dx).astype(int)

J = np.zeros((T + 1, len(X)))
policy = np.zeros((T, len(X)))
J[T] = terminal_value(X)

for t in reversed(range(T)):
    x, u = X[:, None], U[None, :]                  # all (storage, release) pairs
    x_next = np.minimum(x + w[t] - u, x_max)       # water above capacity spills
    feasible = x + w[t] - u >= x_min               # cannot draw below dead storage
    value = p[t] * kappa * head(x) * u + J[t + 1][index(np.maximum(x_next, x_min))]
    value = np.where(feasible, value, -np.inf)
    J[t] = value.max(axis=1)
    policy[t] = U[value.argmax(axis=1)]

def simulate(x0):
    storage = [x0]
    for t in range(T):
        u = policy[t, index(storage[-1])]
        storage.append(min(storage[-1] + w[t] - u, x_max))
    return np.array(storage)

print(f"optimal revenue from x0 = 1100 hm³: {J[0, index(1100)] / 1e6:.1f} M$")

optimal revenue from x0 = 1100 hm³: 74.1 M$

The backward pass evaluates $52 \times 181 \times 23 \approx 2.2 \times 10^5$ (week, storage, release) triples and takes a fraction of a second. Its output is not one trajectory but two tables, $J_t^\star(x)$ and $u_t^\star(x)$, with one entry for every week and storage level.

from pathlib import Path
import sys

code_directory = Path.cwd() / "code"
if str(code_directory) not in sys.path:
    sys.path.insert(0, str(code_directory))

from hydro_reservoir import (
    format_comparison,
    make_interpolation_figure,
    make_policy_figure,
    make_threshold_figure,
    simulate_rule,
)

water_value = np.gradient(J, dx, axis=1)       # dJ_t/dx on the grid, $ per hm³
rollouts = [simulate(x0) for x0 in (300, 700, 1100, 1500, 1900)]
reservoir = dict(w=w, p=p, head=head, terminal_value=terminal_value, x_min=x_min,
                 x_max=x_max, u_min=u_min, u_max=u_max, kappa=kappa)

make_policy_figure(w, p, X, policy, water_value[:T], rollouts)

The left map is the policy. It has two regions: release at turbine capacity or release the ecological minimum, with a thin boundary between them. The white storage paths show the policy the way an operator would read it. Whatever the level in January, the reservoir is drawn down to dead storage by late March, refilled by the freshet and then more slowly through the summer until it is full in mid-September, held full until late November, and drawn down at turbine capacity in December, when prices are high again. The initial storage affects only the first few weeks, because every path reaches the same rule curve before the freshet. A fuller reservoir in January simply means more water to sell during the winter drawdown.

The right map is the derivative of the value function with respect to storage, which hydro operators call the water value. It is lowest, close to zero, at full storage during the freshet: the reservoir cannot hold more, so an extra hm³ would spill and earn nothing. It is highest in winter at low storage, when every hm³ can be sold at the winter price and none is at risk of spilling.

Why the policy is all or nothing. Fix the week $t$ and the storage $x$, and look at $Q_t(x, u)$ as a function of the release $u$ alone. It has two terms, and $u$ enters each one once.

The first term, $p_t\,\kappa\,h(x)\,u$, is linear in $u$. The head $h(x)$ depends on the storage at the start of the week, which is fixed, so each extra hm³ released earns the same amount, $p_t\,\kappa\,h(x)$.

The second term depends on $u$ only through next week’s storage $x_{t+1} = x + w_t - u$. Releasing one more hm³ lowers $x_{t+1}$ by one hm³, so by the chain rule this term changes at the rate $-\partial J_{t+1}^\star/\partial x$ evaluated at $x_{t+1}$, minus the water value. (If the reservoir overflows, $x_{t+1} = \bar x$ does not depend on $u$, and this rate is zero: the extra hm³ would have spilled anyway.)

Adding the two rates gives the slope of $Q_t$ in $u$:

We maximize $Q_t$ over the interval $[u_{\min}, u_{\max}]$, and the sign of this slope tells us where the maximum is:

• If the slope is positive for every $u$ in the interval, $Q_t$ increases from one end to the other, and the maximum is at $u = u_{\max}$.

• If it is negative for every $u$, $Q_t$ decreases, and the maximum is at $u = u_{\min}$.

• Only if the slope changes sign inside the interval can the maximum be strictly between the bounds, at a release where the two values are equal.

The first term of the slope does not depend on $u$ at all. The second depends on $u$ only through the water value, and over the storage levels reachable in one week the water value changes little compared with the gap between it and the price. For most $(t, x)$ the slope therefore keeps one sign over the whole interval, and the release is at a bound. This is the two-region map of Figure 2. The interior case is the thin boundary between the regions: there the policy releases exactly enough water to bring $x_{t+1}$ to the level at which the water value equals the price.

Two physical limits also produce releases between $u_{\min}$ and $u_{\max}$, and both show up in the figures. At low storage the constraint $x + w_t - u \ge x_{\min}$ caps the release below turbine capacity, so “release as much as possible” means emptying the reservoir to dead storage. At full storage the slope jumps at the release $u = x + w_t - \bar x$: below it the extra water would spill, the second term is zero and the slope is positive; above it the slope is the price minus the water value at $\bar x$. When that is negative, the best release is the kink itself, which keeps the reservoir exactly full and releases the inflow. This is what happens from September to November in Figure 3.

The decision rule is therefore a price comparison: run the turbines when the price exceeds the value of the water in the reservoir, both expressed per hm³. In hydro-dominated power systems such as Norway’s, producers compute water values with dynamic programming and use them to set the prices at which they offer their energy.

The same structure appears in continuous time, where the slope is called a switching function and a control that sits at its bounds is called bang-bang. Whenever the reward is linear in a bounded control, the optimum is at a bound except where the switching function is zero, the singular case.

x0 = 1100
comparison = {
    "dynamic programming": simulate_rule(lambda t, x: policy[t, index(x)], x0, **reservoir),
    "run of river (release = inflow)": simulate_rule(lambda t, x: w[t], x0, **reservoir),
    "full release when price > mean": simulate_rule(
        lambda t, x: u_max if p[t] > p.mean() else u_min, x0, **reservoir),
}
print(format_comparison(comparison))

optimal = comparison["dynamic programming"]
marginal_release = p * kappa * head(optimal["storage"][:-1])
marginal_storage = np.array([np.interp(optimal["storage"][t + 1], X, water_value[t + 1])
                             for t in range(T)])
make_threshold_figure(p, optimal["storage"], optimal["release"], marginal_release,
                      marginal_storage, u_min, u_max)

The table above the figure compares the policy with two simpler operating rules over the same year. Running the river as it comes, releasing each week’s inflow, earns about 10% less because it sells most of the water during the freshet, when prices are low. Releasing fully whenever the price is above its annual mean earns even less and spills 900 hm³, because the rule ignores the storage level. It empties the reservoir by early February and has only the small winter inflow to sell for the rest of the high-price season. In May, when the price falls below its mean, it switches to the minimum release in the middle of the freshet. The reservoir is full by late June, and from July to November every hm³ above the minimum flow spills. Dynamic programming accounts for the price, the head and the remaining capacity at the same time, through one function of storage.

Explore the solution. The page below re-solves the same backward recursion in your browser whenever a parameter changes. Backward pass animates the recursion: the maps fill in from week 52 back to week 1, one column per Bellman update, and the panel beside the map shows the value function of the column being computed. Play the year then operates the reservoir forward from the chosen initial storage. The page compares the two terms of the derivative per MWh rather than per hm³: dividing both by $\kappa\,h(x_t)$ turns the first into the price $p_t$ and the second into a water value in $/MWh.

Question. Set the price of water left at the end of the year to zero. What does the policy do in November and December? Is this what an operator who will still run the dam next year should do?

Question. Starting from 1100 hm³, the policy keeps the minimum release through January, although the price is already above 80 $/MWh, and releases at full capacity from mid-February to mid-March. Why does it wait?

Question. Suppose the river had three such reservoirs in series, each discretized with the same 181 storage levels and 23 releases. How many (storage, release) combinations would one Bellman update have to evaluate at each week?

Handling Continuous Spaces with Interpolation¶

When a successor falls between the states stored in the table, which value should the recursion use for it, and what problem does the recursion then solve?

Theorem 1 places no restriction on the state space. A program, however, has to store $J_{t+1}^\star$ as finitely many numbers. The reservoir code stores one value per 10 hm³ of storage, and it avoided the question by construction: its inflows were rounded to multiples of the grid step, so every successor $x + w_t - u$ landed on a grid point. Measured inflows are not multiples of anything. In a week with 63.7 hm³ of inflow and a 30 hm³ release, a reservoir at 1000 hm³ moves to 1033.7 hm³, between two grid points, and Bellman’s equation then asks for $J_{t+1}^\star(1033.7)$, which the table does not contain. Choosing a value for it is a representation problem that the collocation chapter already treated for functions of time.

Nodal Representation of the Cost-to-Go¶

In collocation, a trajectory piece was stored by its values $y_j$ at support nodes $\sigma_j$ and reconstructed between them by cardinal functions, $p(\tau)=\sum_j y_j\ell_j(\tau)$ with $\ell_j(\sigma_i)=\delta_{ij}$. A cost-to-go table is the same kind of object, with the nodes placed in the state space instead of on the time axis. Choose grid states $\mathbf{x}^{(1)},\dots,\mathbf{x}^{(N)}$ and let $v_{t,i}$ denote the value stored for $\mathbf{x}^{(i)}$ at stage $t$. A family of cardinal functions on the grid turns these numbers into a function of the state:

The array J in the reservoir code holds the nodal coordinates $v_{t,i}$. An interpolation rule is a choice of cardinal functions, which the array alone does not specify.

When the collocation chapter chose how to represent a control between its support values, it listed three options: a constant command on each interval, a line between endpoint commands, or a higher-degree polynomial. The same three options are available for the cost-to-go. On a one-dimensional grid $x^{(1)}<\dots<x^{(N)}$, such as the storage levels of the reservoir, they are:

1. Rounding (nearest-neighbor interpolation). The cardinal function $\ell_i$ equals one on the states closer to $x^{(i)}$ than to any other node and zero elsewhere, so $\hat J_t(x)$ is the value stored at the nearest node. It is constant on each cell, like a zero-order-hold control, except that the cell is centered on its node.

2. Linear interpolation. Each $\ell_i$ is a “hat” that rises from zero at $x^{(i-1)}$ to one at $x^{(i)}$ and returns to zero at $x^{(i+1)}$. Between two adjacent nodes only two hats are nonzero, and $\hat J_t(x)=(1-\lambda)\,v_{t,i}+\lambda\,v_{t,i+1}$ with $\lambda=(x-x^{(i)})/(x^{(i+1)}-x^{(i)})\in[0,1]$. These are the continuous piecewise-linear controls with shared endpoint values of the trapezoidal transcription.

3. Cubic spline. Cubic pieces are joined with continuous first and second derivatives. Each cardinal function is then nonzero over the whole grid, with lobes of alternating sign that decay away from its node.

In all three cases the weights $\ell_i(x)$ sum to one, because each rule reproduces a constant function exactly. They differ in sign. Rounding and linear interpolation have nonnegative weights, so $\hat J_t(x)$ is a convex combination of stored values and stays within their range. The collocation chapter used the same property to bound a piecewise-linear control by bounding its support values. The spline has negative weights, like the quadratic Lagrange basis of the classroom exercise whose weights at $\tau=1/4$ were $(3/8,3/4,-1/8)$. It can therefore leave the range of the stored values, as the bottom row of Figure 5 shows on data with a jump. In exchange, it is far more accurate than linear interpolation on smooth functions.

from finite_horizon_dp_figures import make_cardinal_functions_figure

make_cardinal_functions_figure()

Substituting (15) into Bellman’s equation and evaluating it only at the grid states gives a recursion on $N$ numbers per stage:

Compared with Algorithm 1, the loop over states becomes a loop over grid states, and the lookup of $J_{t+1}^\star$ at the successor becomes an evaluation of $\hat J_{t+1}$. Each stage still consists of independent one-step problems, one per grid state. For the reservoir, which maximizes revenue, the minimum becomes a maximum. The numbers $v_{t,i}$ are approximations of $J_t^\star(\mathbf{x}^{(i)})$, and the cardinal functions determine how good they are.

Grid Models Behind the Interpolated Recursion¶

The weights $\ell_j(\mathbf{f}_t(\mathbf{x}^{(i)},\mathbf{u}))$ in (16) depend only on where the successor lands relative to the grid. Read as a description of where the system goes next, they show which problem the interpolated recursion solves.

With rounding, exactly one weight equals one, the weight of the node nearest to the successor. The recursion is then the backward recursion of Theorem 1 for a different deterministic system, which lives on the grid and moves to the node nearest to $\mathbf{f}_t(\mathbf{x},\mathbf{u})$. The theorem applies to this rounded system without change, so the recursion returns its optimal values and decisions. The approximation error lies entirely in the replacement of the model, and that replacement can be qualitatively wrong. Any change smaller than half a grid step disappears in the rounded system. In the example below, the weekly gain of storage in summer is such a change, and the rounded reservoir cannot fill.

With linear interpolation, the successor’s weight is split between the two nodes that bracket it, as $1-\lambda$ and $\lambda$. These numbers are nonnegative and sum to one, so they can be read as probabilities. The bracket in (16) becomes the stage cost plus the expected stored value at a random next node, and the interpolated recursion is the exact backward recursion of a finite-state system with random successors: from $x^{(i)}$ under $u$, the system moves to $x^{(j)}$ with probability $\ell_j(f_t(x^{(i)},u))$. In the reservoir, a successor of 1033.7 hm³ on the 10 hm³ grid becomes 1030 hm³ with probability 0.63 and 1040 hm³ with probability 0.37. Hat functions reproduce linear functions, $\sum_j\ell_j(x)\,x^{(j)}=x$, so the random next node equals the true successor on average. Linear interpolation thus replaces the deterministic model by a Markov chain on the grid. Matching the mean motion of such a chain to the continuous dynamics is the local-consistency condition at the base of the Markov chain approximation method of Kushner & Dupuis (2001) for continuous-time stochastic control. In $n$ dimensions, multilinear interpolation on a rectangular grid multiplies one-dimensional hat weights. Its $2^n$ weights on the corners of the enclosing cell are again nonnegative and sum to one, so the same reading holds. The next chapter makes the inflow itself random and interpolates each of its successors linearly. The inflow probabilities and the interpolation weights then multiply, and the result is again a distribution over grid states.

A cubic spline places negative weights on some nodes, and negative numbers cannot be read as probabilities. No finite-state problem has (16) as its exact recursion. The interpolated value $\hat J_{t+1}$ can also leave the range of the stored values, and the minimization over $\mathbf{u}$ is free to steer toward successors where it does.

Error Propagation Across Stages¶

Each stage of (16) uses stored values that already contain the error of later stages. The sign of the weights decides whether that error is passed on unchanged or amplified. To make this precise, define the error of the nodal values and the interpolation error of the exact cost-to-go,

where the supremum ranges over the successors reached from the grid. The constant $\Lambda=\sup_{\mathbf{y}}\sum_j|\ell_j(\mathbf{y})|$ bounds how much the interpolant can magnify errors in the stored values. Because the weights sum to one, $\Lambda\ge 1$, with equality exactly when they are nonnegative.

With nonnegative weights, every stage adds its own interpolation error and no earlier error is amplified, so the error grows at most linearly with the horizon. With negative weights, the bound allows each earlier error to be multiplied by $\Lambda>1$ at every subsequent stage. The size of $\delta_t$ depends on the rule and the grid spacing $\Delta$. For a cost-to-go with bounded derivatives on a one-dimensional grid, rounding gives $\delta\le\tfrac{\Delta}{2}\max|J'|$ and linear interpolation gives $\delta\le\tfrac{\Delta^2}{8}\max|J''|$, so halving the spacing halves the first bound and quarters the second. A cubic spline converges faster still on smooth functions, which is why it can be the most accurate of the three in practice even though its $\Lambda$ exceeds one.

The proof compares the two minima through the largest difference between their brackets, not an average one, and the example below shows why that is the right comparison. An interpolation error that looks unbiased across controls becomes biased once a control is selected, because the minimization prefers controls whose successors happen to be evaluated favorably. On the coarse reservoir grid, rounding moves the feasible successors by -0.3 hm³ on average, over all weeks, grid states and releases. The successors that the maximization selects are moved by +45 hm³ on average, out of a possible 50, and $94\%$ of them are rounded up. The same selection effect reappears as overestimation bias in fitted Q-iteration.

These bounds do not remove the curse of dimensionality from the opening of this chapter. A grid with $N$ levels per coordinate still has $N^n$ nodes, each with its own one-step problem, and multilinear interpolation touches $2^n$ stored values per evaluation. A more accurate rule reduces the number of levels needed for a given error, not the exponential dependence on $n$. The approximate value methods replace the grid by parametric approximators and face the same question about weights. Interpolation with nonnegative weights that sum to one is an averager in Gordon’s sense and keeps value iteration stable, while least-squares fits and neural networks carry no such guarantee.

Decisions Between Grid Points¶

The recursion returns decisions at the grid states only, but the reservoir will be measured between them. One option is to interpolate the table of decisions $\boldsymbol{\pi}_{t,i}$. The reservoir’s release table is all or nothing, as Figure 2 shows, so in each week it jumps from the ecological minimum to turbine capacity at a switching storage. Linear interpolation of that table stays within the release limits but returns intermediate releases that the recursion never evaluated. A cubic spline through it overshoots at the switch, as in the bottom row of Figure 5, and commands releases below the ecological minimum and above turbine capacity.

The other option keeps the decision inside an optimization. At the measured state, solve the one-step problem with the interpolant standing in for the future:

At a grid state this reproduces the stored decision. Between grid states it returns an admissible control that is consistent with the values the recursion computed. The rule (20) is a receding-horizon controller with a prediction horizon of one step and terminal cost $\hat J_{t+1}$. In the receding-horizon chapter, the terminal cost stood in for the cost-to-go beyond the prediction horizon and had to be supplied by the designer. Here the interpolated backward recursion computes it. This sharpens the answer to the question that opened the chapter. What can be computed in advance for a whole family of states is the cost-to-go, stored at grid states. The decision at a particular state still comes from an online minimization, now over a single step. A longer online horizon makes the decision depend less on $\hat J$, at the price of a larger online problem.

Example: Reservoir Operation with Unrounded Inflows¶

We return to the reservoir with the inflow forecast left unrounded, so that successors fall between grid points, and solve it on a coarse grid of 19 storage levels spaced 100 hm³ apart. Two versions of Algorithm 2 differ only in the interpolation rule passed to them: nearest neighbor or linear (np.interp). A reference solution uses linear interpolation on a grid twenty times finer. Each value function is then used the way an operator would use it: the reservoir is simulated with the exact, unrounded dynamics, and each week’s release is chosen by the one-step lookahead (20) with the interpolated $\hat J_{t+1}$.

# Unrounded inflows: successors now fall between grid points.
w_raw = 60 + 340 * np.exp(-((weeks - 19) / 3.5) ** 2) + 60 * np.exp(-((weeks - 43) / 4) ** 2)
U_fine = np.linspace(u_min, u_max, 45)

def nearest(x_query, grid, values):
    step = grid[1] - grid[0]
    i = np.clip(np.round((x_query - grid[0]) / step).astype(int), 0, len(grid) - 1)
    return values[i]

def linear(x_query, grid, values):
    return np.interp(x_query, grid, values)

def backward_recursion(grid, interpolate):
    J = np.zeros((T + 1, len(grid)))
    J[T] = terminal_value(grid)
    for t in reversed(range(T)):
        x, u = grid[:, None], U_fine[None, :]
        x_next = np.clip(x + w_raw[t] - u, x_min, x_max)
        value = p[t] * kappa * head(x) * u + interpolate(x_next, grid, J[t + 1])
        J[t] = np.where(x + w_raw[t] - u >= x_min, value, -np.inf).max(axis=1)
    return J

def operate(grid, J, interpolate, x0):
    """Simulate the true dynamics, choosing each release by one Bellman step."""
    storage, revenue = [x0], 0.0
    for t in range(T):
        x = storage[-1]
        x_next = np.clip(x + w_raw[t] - U_fine, x_min, x_max)
        value = p[t] * kappa * head(x) * U_fine + interpolate(x_next, grid, J[t + 1])
        u = U_fine[np.argmax(np.where(x + w_raw[t] - U_fine >= x_min, value, -np.inf))]
        revenue += p[t] * kappa * head(x) * u
        storage.append(min(x + w_raw[t] - u, x_max))
    return np.array(storage), revenue + terminal_value(storage[-1])

coarse = np.arange(x_min, x_max + 100, 100)   # 19 levels
fine = np.arange(x_min, x_max + 5, 5)         # 361 levels
schemes = {
    "fine grid, linear (reference)": (fine, linear),
    "coarse grid, nearest neighbor": (coarse, nearest),
    "coarse grid, linear": (coarse, linear),
}
outcomes = {}
for name, (grid, interpolate) in schemes.items():
    J_scheme = backward_recursion(grid, interpolate)
    storage, revenue = operate(grid, J_scheme, interpolate, 1100)
    outcomes[name] = (grid, J_scheme[0], storage)
    print(f"{name:<32} predicted {linear(1100, grid, J_scheme[0]) / 1e6:6.2f} M$"
          f"   realized {revenue / 1e6:6.2f} M$")

fine grid, linear (reference)    predicted  74.16 M$   realized  74.17 M$
coarse grid, nearest neighbor    predicted  98.46 M$   realized  71.06 M$
coarse grid, linear              predicted  74.02 M$   realized  74.15 M$

make_interpolation_figure(T, outcomes, "fine grid, linear (reference)")

Linear interpolation on 19 points predicts a revenue within 0.2% of the 361-point reference and realizes one within 0.03%. Rounding on the same grid fails in the two ways anticipated above, and both come from the size of a typical weekly change. Its grid model is wrong in summer. The reservoir gains about 30 hm³ per week at the minimum release, less than half the grid step, so the rounded successor equals the current grid point. In the rounded model, storing water gains nothing, so the policy releases the summer inflow at summer prices instead of keeping it for winter, and the reservoir never fills. Its value function is also too optimistic, by about a third, because the maximization selects the releases whose successors round up and so adds water that does not exist. Linear interpolation keeps the fractional part of every step, so small changes accumulate correctly over the year.

Question. Replace linear by a cubic spline, lambda x_query, grid, values: CubicSpline(grid, values)(x_query) with CubicSpline from scipy.interpolate, and rerun the coarse-grid recursion. Is the value function more or less accurate than with linear interpolation? Does the answer contradict Proposition 1?

Question. Instead of solving the one-step problem at the measured storage, an operator interpolates the coarse-grid release table between grid states. Record the maximizing release in backward_recursion for the linear rule, then operate the reservoir from 1100 hm³ by interpolating that table, first linearly and then with a cubic spline. What goes wrong in each case?

Linear Quadratic Regulator via Dynamic Programming¶

What form do the value function and feedback law take when the dynamics are linear and every cost is quadratic?

Linear dynamics and quadratic costs give a backward recursion that can be solved in closed form. The value function remains quadratic at every stage, and the optimal policy is a linear feedback law. No state grid, interpolation scheme, or general function approximator is needed. The recursion tracks a finite sequence of matrices.

Consider a discrete-time linear system:

where $\mathbf{x}_t \in \mathbb{R}^n$ is the state and $\mathbf{u}_t \in \mathbb{R}^m$ is the control input. The matrices $A_t \in \mathbb{R}^{n \times n}$ and $B_t \in \mathbb{R}^{n \times m}$ describe the system dynamics at time $t$.

The cost function to be minimized is quadratic:

where $Q_T \succeq 0$ (positive semidefinite), $Q_t \succeq 0$, and $R_t \succ 0$ (positive definite) are symmetric matrices of appropriate dimensions. The positive definiteness of $R_t$ ensures the minimization problem is well-posed.

A quadratic terminal cost implies a quadratic value function at every earlier stage. Suppose the optimal cost-to-go at stage $t$ has the form

for some positive semidefinite matrix $P_t$. At the terminal time, this is true by definition: $P_T = Q_T$.

Backward induction verifies the hypothesis. Assume $J_{t+1}^\star(\mathbf{x}_{t+1}) = \frac{1}{2}\mathbf{x}_{t+1}^\top P_{t+1} \mathbf{x}_{t+1}$. Bellman’s equation at stage $t$ is

Substituting the dynamics $\mathbf{x}_{t+1} = A_t\mathbf{x}_t + B_t\mathbf{u}_t$ and the quadratic form for $J_{t+1}^\star$:

Expanding the last term:

The expression inside the minimization becomes:

Collecting terms involving $\mathbf{u}_t$ and calling the result $G_t$:

This is a quadratic function of $\mathbf{u}_t$, and $J_t^\star(\mathbf{x}_t) = \min_{\mathbf{u}_t} G_t(\mathbf{x}_t,\mathbf{u}_t)$. To find the minimizer, we set the gradient of $G_t$ with respect to $\mathbf{u}_t$ to zero:

Since $R_t + B_t^\top P_{t+1} B_t$ is positive definite (both $R_t$ and $P_{t+1}$ are positive semidefinite with $R_t$ strictly positive), we can solve for the optimal control:

Define the gain matrix:

so that $\mathbf{u}_t^\star = -K_t\mathbf{x}_t$. This is a linear feedback policy: the optimal control is simply a linear function of the current state.

Substituting $\mathbf{u}_t^\star$ back into the cost-to-go expression and simplifying (by completing the square), we obtain:

where $P_t$ satisfies the discrete-time Riccati equation:

The resulting backward recursion is:

The iLQR and DDP chapter applies this recursion to nonlinear problems. Linearizing the dynamics and quadratically approximating the costs along a nominal trajectory gives a local linear-quadratic problem, whose Riccati recursion yields a correction to the trajectory. Repeating that step gives iLQR and differential dynamic programming.

Local Stabilization of the Cart-Pole¶

The cart-pole in the trajectory-optimization chapter started from the downward configuration and required a large nonlinear maneuver to reach the top. Once it is near the upright equilibrium, a smaller problem remains: reject local deviations by moving the cart in response to the measured state. The nonlinear state and input are the same as before,

with $\theta=0$ at upright. One step of fourth-order Runge--Kutta integration with step size $h$ defines the discrete map $\mathbf{x}_{k+1}=F_h(\mathbf{x}_k,u_k)$, as in the trajectory-optimization chapter. The upright configuration with zero input is a fixed point of this map: $F_h(\mathbf{x}^\star,u^\star)=\mathbf{x}^\star$ with $(\mathbf{x}^\star,u^\star)=(0,0)$.

Write the state and input as deviations from this point, $\mathbf{x}_k=\mathbf{x}^\star+\delta\mathbf{x}_k$ and $u_k=u^\star+\delta u_k$. A first-order Taylor expansion of $F_h$ around $(\mathbf{x}^\star,u^\star)$ gives

where $A$ and $B$ are the Jacobians of the discrete map at the equilibrium:

Because $F_h(\mathbf{x}^\star,u^\star)=\mathbf{x}^\star$, subtracting $\mathbf{x}^\star$ from both sides and dropping the second-order remainder leaves the linear deviation dynamics

The approximation is good only while $\delta\mathbf{x}_k$ and $\delta u_k$ stay small, which is what the third rollout below tests. We differentiate the RK4 step itself, not the continuous-time vector field, so $A$ and $B$ already account for the sampling at step size $h$ and need no separate discretization. The code computes both Jacobians with central finite differences of rk4_step.

The experiment uses $h=0.02$ s, $Q=\operatorname{diag}(2,0.2,80,3)$, and $R=0.15$. For this time-invariant infinite-horizon case, the Riccati recursion converges to a fixed matrix $P$ that satisfies the discrete algebraic Riccati equation. The corresponding policy is $u_k=-K\delta\mathbf{x}_k$.

The unconstrained linear closed loop is asymptotically stable when every eigenvalue of $A-BK$ lies inside the unit disk. The physical implementation adds two constraints that are absent from that eigenvalue calculation: acceleration is clipped to $|u|\leq 8\;\mathrm{m\,s^{-2}}$, and the cart must remain inside a 2.4 m rail. Three deterministic nonlinear rollouts distinguish the claims supported by the linearization:

1. An uncontrolled pole begins $5^\circ$ from upright.

2. The LQR controller begins from the same $5^\circ$ displacement.

3. The same controller begins $45^\circ$ from upright with the same actuator and rail limits.

from pathlib import Path
import sys

code_directory = Path.cwd() / "code"
if str(code_directory) not in sys.path:
    sys.path.insert(0, str(code_directory))

from cartpole_control import (
    CartPoleParameters,
    design_lqr,
    format_lqr_metrics,
    make_lqr_animation,
    make_lqr_figure,
    run_lqr_cases,
)

cartpole_parameters = CartPoleParameters()
lqr_design = design_lqr(cartpole_parameters)
lqr_cases = run_lqr_cases(cartpole_parameters)

print(format_lqr_metrics(lqr_cases, lqr_design))
make_lqr_figure(lqr_cases, cartpole_parameters)

The closed-loop eigenvalues certify asymptotic stability of the unconstrained linearized model, not every trajectory of the nonlinear constrained plant. The $5^\circ$ rollout remains in a region where the linear approximation supplies useful actions. In the $45^\circ$ rollout, the controller immediately commands more acceleration than the actuator can provide, and the cart then reaches the end of the rail. Increasing the entries of $Q$ cannot remove those physical limits.

from IPython.display import HTML, display
import matplotlib.pyplot as plt

lqr_animation = make_lqr_animation(lqr_cases, cartpole_parameters)
display(HTML(lqr_animation.to_jshtml()))
plt.close(lqr_animation._fig)

Balancing a pen on a finger motivates the choice of cart acceleration as the input, because the finger stabilizes the object by moving its base. The cart-pole model replaces the contact by a planar frictionless hinge. A real pen can slip, detach, flex, and rotate out of the plane, while sensing and hand motion introduce delays. The calculation establishes local stabilization for the stated rigid-body model; the classroom demonstration shares its instability and feedback mechanism, not all of its equations.

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
def linearize_upright(
    parameters: CartPoleParameters = CartPoleParameters(),
    *,
    step_size: float = 0.02,
    epsilon: float = 1e-6,
) -> tuple[np.ndarray, np.ndarray]:
    """Linearize the exact discrete RK4 update at the upright equilibrium."""

    equilibrium = np.zeros(4)
    state_matrix = np.empty((4, 4), dtype=float)
    for column in range(4):
        perturbation = np.zeros(4)
        perturbation[column] = epsilon
        state_matrix[:, column] = (
            rk4_step(equilibrium + perturbation, 0.0, parameters, step_size)
            - rk4_step(equilibrium - perturbation, 0.0, parameters, step_size)
        ) / (2.0 * epsilon)
    input_matrix = (
        rk4_step(equilibrium, epsilon, parameters, step_size)
        - rk4_step(equilibrium, -epsilon, parameters, step_size)
    )[:, None] / (2.0 * epsilon)
    return state_matrix, input_matrix


def design_lqr(
    parameters: CartPoleParameters = CartPoleParameters(),
    *,
    step_size: float = 0.02,
) -> LQRDesign:
    """Solve the discrete algebraic Riccati equation at the upright state."""

    state_matrix, input_matrix = linearize_upright(parameters, step_size=step_size)
    cost_matrix = np.diag([2.0, 0.2, 80.0, 3.0])
    control_cost = np.array([[0.15]])
    riccati_matrix = solve_discrete_are(
        state_matrix,
        input_matrix,
        cost_matrix,
        control_cost,
    )
    gain = np.linalg.solve(
        control_cost + input_matrix.T @ riccati_matrix @ input_matrix,
        input_matrix.T @ riccati_matrix @ state_matrix,
    )
    eigenvalues = np.linalg.eigvals(state_matrix - input_matrix @ gain)
    return LQRDesign(
        state_matrix,
        input_matrix,
        cost_matrix,
        control_cost,
        riccati_matrix,
        gain,
        eigenvalues,
    )

Download the shared nonlinear cart-pole and LQR source.

Summary and Outlook¶

Backward recursion replaces one trajectory optimization with a sequence of state-indexed subproblems. Storing the cost-to-go at grid states and choosing cardinal functions to interpolate between them makes the recursion a finite computation. With nonnegative weights, the interpolated recursion is the exact recursion of a grid model, deterministic for rounding and stochastic for linear interpolation, and its errors add across stages instead of multiplying. A decision between grid states comes from a one-step lookahead against the interpolated cost-to-go. The linear-quadratic case needs no grid: the value function stays quadratic, the recursion becomes the Riccati recursion, and the policy is linear state feedback.

The iLQR and DDP chapter uses this Riccati recursion as a local step for nonlinear trajectory optimization. It linearizes the dynamics around a nominal trajectory, solves the resulting linear-quadratic problem backward, and tests the correction through the nonlinear simulator.

Linear interpolation already replaced each deterministic successor by a distribution over two grid states. How do the value and policy change when the model itself prescribes a distribution of possible successors? Stochastic dynamic programming replaces the deterministic continuation value by a conditional expectation.