Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Classical Planning and Optimal Control › Linear-Quadratic Control

Michael BrenndoerferJune 10, 202613 min read

World Models HandbookAll chapters

Part 2 of 12Inference, Decisions, and Control

Chapter 10Classical Planning and Optimal Control › Linear-Quadratic Control

Reading section 5 of 10 · Linear-Quadratic Control · 13 minSection 5 of 10 · 13 min

Section 5 of 10 · 13 min read

Classical Planning and Optimal Control › Linear-Quadratic ControlOne section · Sections in this chapter5 / 10 ↓

Sections in this chapter10

Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Finite tabular DP needs an enumerable state set. Continuous-state dynamic programming also exists, but representing an arbitrary value function is difficult. The linear-quadratic regulator (LQR) is a foundational, widely used special case in which the Bellman equation reduces to Riccati equations. Iterative LQR and differential dynamic programming, in turn, build local trajectory updates from LQR-like subproblems.

The setting is:

Linear dynamics (from Chapter 3's system identification, or first-principles modeling):
$x_{t+1} = A x_t + B u_t,$
where $x_t$ is the state, $u_t$ the control input, $A$ the state-transition matrix describing how the current state carries into the next, and $B$ the input matrix describing how the control changes the next state.
Quadratic cost, penalizing deviation from zero and control effort:
$J = \sum_{t=0}^{T-1} \left( x_t^\top Q x_t + u_t^\top R u_t \right) + x_T^\top Q_f x_T,$
where:
- $x_t^\top Q x_t$ : the per-step cost in state directions weighted by $Q$ ; it is positive for every nonzero deviation only when $Q \succ 0$
- $u_t^\top R u_t$ : the per-step control cost, penalizing control effort
- $x_T^\top Q_f x_T$ : the terminal cost on the final state $x_T$
- $Q$ , $R$ , $Q_f$ : symmetric cost matrices with $Q, Q_f \succeq 0$ and $R \succ 0$ , so state costs are nonnegative while control is strictly penalized

Why quadratic? Squared error is a natural penalty for regulation, and, together with linear dynamics, it gives a value-function family that is closed under Bellman backups: if $V(x) = x^\top P x$ , the next backup is quadratic too. A probabilistic interpretation is available when the quadratic terms are proportional to negative log densities of Gaussian state preferences and control priors, but that interpretation needs those modeling assumptions; it is not intrinsic to every LQR problem. The closure property is the computational trick: we solve for the finitely many entries of $P$ rather than for an arbitrary function.

Deriving the Riccati equationLink Copied

Guess that the value function is quadratic, $V(x) = x^\top P x$ , and run one step of the Bellman backup:

x^\top P' x = \min_u \left[ x^\top Q x + u^\top R u + (Ax + Bu)^\top P (Ax + Bu) \right].

where:

$P$ : the matrix of the current quadratic value estimate $V(x) = x^\top P x$
$P'$ : the updated matrix after one backup, defining the next value iterate
$Ax + Bu$ : the model's prediction of the next state under the linear dynamics
$\min_u$ : optimization over all candidate controls $u$ , selecting the one with the smallest one-step cost plus successor value

The expression inside the min is a quadratic in $u$ , so it can be minimized analytically by setting its gradient to zero:

\begin{aligned} 0 &= \nabla_u \left[ u^\top R u + (Ax + Bu)^\top P (Ax + Bu) \right] && \text{(first-order condition; } x^\top Q x \text{ does not involve } u \text{)} \\ &= 2 R u + 2 B^\top P (Ax + Bu) && \text{(differentiate each quadratic term)} \\ &= 2 (R + B^\top P B) u + 2 B^\top P A x, && \text{(collect the terms multiplying } u \text{)} \end{aligned}

Solving this linear equation for $u$ gives the optimal feedback law

u^*(x) = -\underbrace{(R + B^\top P B)^{-1} B^\top P A}_{K}\, x.

Two things just happened that are worth savoring. First, the optimal policy came out linear in the state, with the gain matrix $K$ as the only unknown. We did not restrict the policy class; we allowed arbitrary (measurable) functions from state to action, and the optimizer chose a linear one because the problem's structure made it optimal. Second, the minimizing $u$ can be computed in closed form, no search over actions needed. Compare this with the grid world, where every action had to be enumerated and evaluated at every state. Here the entire maximization over a continuum of controls collapses to solving a small linear system. Structure bought us the same thing the heuristic bought A*: a way to skip most of the work. Now substitute $u^*(x) = -Kx$ back into the backup, so that $Ax + Bu^* = (A - BK)x$ :

\begin{aligned} x^\top P' x &= x^\top Q x + x^\top K^\top R\, K x + x^\top (A - B K)^\top P (A - B K) x && \text{(insert } u^* = -Kx \text{)} \\ &= x^\top \left[ Q + A^\top P A - A^\top P B (R + B^\top P B)^{-1} B^\top P A \right] x. && \text{(expand and collect terms using } K = (R + B^\top P B)^{-1} B^\top P A \text{)} \end{aligned}

Reading off the matrix and imposing the fixed-point condition $P' = P$ yields

P = Q + A^\top P A - A^\top P B (R + B^\top P B)^{-1} B^\top P A.

This is the discrete algebraic Riccati equation (DARE). The finite-horizon version is a backward recursion in time: start with $P_T = Q_f$ and step the recursion backward to $t=0$ . The infinite-horizon equation has a unique positive semidefinite stabilizing solution under the usual stabilizability and detectability conditions, and the resulting modeled closed-loop matrix $A-BK$ is Schur stable. Along its trajectories,

V(x_{t+1})-V(x_t)=-x_t^\top(Q+K^\top R K)x_t.

This is nonpositive in general and strictly negative for every nonzero $x_t$ when $Q+K^\top RK$ is positive definite. Under those conditions, $V(x)=x^\top P x$ is a Lyapunov certificate for the modeled closed loop. The theorem does not automatically certify a different true plant or an approximate nonlinear system.

Certainty equivalence

LQR assumes known $A$ and $B$ . A common practice when they are estimated is certainty equivalence: fit the model, then design the controller as if the estimate were exact. In the standard partially observed LQG setting, the separation principle makes an optimal controller from a Kalman state estimate followed by the corresponding LQR feedback. That theorem does not cover generic parameter-estimation error or nonlinear learned models. In those cases certainty equivalence is a heuristic: an optimizer may seek regions where the fitted model is optimistic, so robustness margins and evaluation on the true plant or a credible uncertainty set matter.

Building an LQR on an identified modelLink Copied

Let's do this end to end. Following Chapter 3, we fit $A$ and $B$ of a double integrator from noisy transition data, solve the Riccati equation on the estimated model, and deploy the resulting gain on the true system. This identify-plan-deploy chain is one common model-based-control pattern. Later systems may also change the representation, uncertainty treatment, planner, data-collection loop, adaptation rule, and safety architecture.

First, generate data and identify the model by least squares.

In[15]:

Code

rng = np.random.default_rng(3)
dt = 0.1
A_true = np.array([[1.0, dt], [0.0, 1.0]])  # position and velocity update
B_true = np.array([[0.0], [dt]])  # control is acceleration

n = 400
X = rng.normal(scale=1.0, size=(n, 2))
U = rng.normal(scale=1.5, size=(n, 1))
noise = rng.normal(scale=0.01, size=(n, 2))
X_next = X @ A_true.T + U @ B_true.T + noise

Z = np.hstack([X, U])  # regressors: state and action
Theta, *_ = np.linalg.lstsq(Z, X_next, rcond=None)
A_hat, B_hat = Theta[:2, :].T, Theta[2:, :].T

Out[16]:

Console

True A:
 [[1.  0.1]
 [0.  1. ]] 
Estimated A:
 [[ 1.0005e+00  9.9900e-02]
 [-5.0000e-04  1.0000e+00]]

True B:
 [[0. ]
 [0.1]] 
Estimated B:
 [[0.0004]
 [0.1003]]

The estimates land within a few thousandths of the true matrices in this seeded experiment. The regression succeeds because the synthetic state-action design matrix $Z=[X,U]$ is full rank and well conditioned; white-noise controls alone would not establish identifiability of the state columns in $A$ . In a rollout dataset, input design and state visitation jointly determine the regressor's rank and conditioning. We will judge whether certainty-equivalent control worked here by checking the true closed loop below, rather than by treating a particular parameter-error magnitude as a universal threshold.

Now solve the Riccati equation by fixed-point iteration, which is nothing more than repeatedly applying the Bellman backup for quadratic value functions until $P$ stops moving.

In[17]:

Code

Q = np.eye(2)  # penalize state deviation
R = 0.5 * np.eye(1)  # penalize control effort
P = Q.copy()
riccati_converged = False
for _ in range(5000):
    K = np.linalg.solve(R + B_hat.T @ P @ B_hat, B_hat.T @ P @ A_hat)
    P_new = Q + A_hat.T @ P @ (A_hat - B_hat @ K)
    if np.max(np.abs(P_new - P)) < 1e-10:
        P = P_new
        riccati_converged = True
        break
    P = P_new
if not riccati_converged:
    raise RuntimeError("Riccati iteration did not converge")
K = np.linalg.solve(R + B_hat.T @ P @ B_hat, B_hat.T @ P @ A_hat)

Compare this loop with value iteration. Both initialize a value representation and iterate a Bellman-derived update toward a fixed point. Here the representation is a two-by-two matrix rather than a table, and analytic minimization replaces action enumeration. They are structurally analogous iterations with different operators, domains, and convergence assumptions.

In[18]:

Code

closed_loop = A_true - B_true @ K
eig_mags = np.abs(np.linalg.eigvals(closed_loop))

Out[19]:

Console

Gain K: [1.27244878 2.09436387]
Closed-loop eigenvalue magnitudes: [0.8963 0.8963]
Stabilizing: True

The gain came out positive on both state components, meaning the controller pushes acceleration opposite to position and velocity. We evaluated the gain on the true matrices and found both eigenvalue magnitudes below one. That is an empirical check for this identified example; the Riccati theorem itself guarantees stability for the estimated model on which the gain was designed. Finally, deploy the gain on the true dynamics, starting from a state with nonzero velocity, and compare against doing nothing.

In[20]:

Code

T = 50
x0 = np.array([0.0, 1.5])
x_un = np.zeros((T + 1, 2))
x_cl = np.zeros((T + 1, 2))
x_un[0], x_cl[0] = x0, x0

for t in range(T):
    u_t = -K @ x_cl[t]
    x_cl[t + 1] = A_true @ x_cl[t] + B_true @ u_t
    x_un[t + 1] = A_true @ x_un[t]

state_cost_un = float(sum(x @ Q @ x for x in x_un))
state_cost_cl = float(sum(x @ Q @ x for x in x_cl))

Out[21]:

Console

Uncontrolled final state:  position +7.500, velocity +1.500
LQR final state:           position +0.011, velocity -0.016
Accumulated state cost, uncontrolled: 1080.6
Accumulated state cost, LQR:          10.87

The uncontrolled cart coasts with its initial velocity over this noiseless simulation, and its accumulated squared-state cost dwarfs the regulated one. This illustrates why feedback matters: the open-loop system has repeated eigenvalues at one, so velocity persists and position drifts, while the two-component gain $K$ makes this tested closed loop forget its initial condition exponentially.

Out[22]:

Visualization

Line chart with regulated position rising briefly and then decaying, while uncontrolled position rises linearly. — Position of the double integrator under LQR versus no control. Initial velocity first raises the regulated position, after which it decays toward zero; the uncontrolled cart drifts linearly.

Line chart with velocity decaying to zero under control and staying constant without. — Velocity of the same two runs. LQR bleeds the initial velocity off smoothly; the uncontrolled velocity never changes.

The controller was designed on the estimated matrices but deployed on the true system. The explicit true-plant eigenvalue check shows a comfortable stability margin for this experiment. Certainty equivalence prescribes the design step; that numerical check supplies the robustness evidence here. On other plants, especially away from a valid local model, parameter error can change stability or performance substantially. The gap between the model used for planning and the world used for execution therefore needs direct evaluation.

The role of the R weight in that trade-off can be made concrete by re-solving the Riccati equation across a range of control penalties.

In[23]:

Code

import numpy as np

## Re-solve the infinite-horizon Riccati equation on the identified double
## integrator model (A_hat, B_hat from above) for a sweep of R values.
## Deterministic fixed-point iteration, tolerance 1e-9, 20000-iteration cap.
r_values = np.logspace(-2, 1.5, 12)
gain_norms = []
eig_mags_best = []
Q_sweep = np.eye(2)
for r_val in r_values:
    R_sweep = r_val * np.eye(1)
    P_s = Q_sweep.copy()
    for _ in range(20000):
        K_s = np.linalg.solve(
            R_sweep + B_hat.T @ P_s @ B_hat, B_hat.T @ P_s @ A_hat
        )
        P_new = Q_sweep + A_hat.T @ P_s @ (A_hat - B_hat @ K_s)
        if np.max(np.abs(P_new - P_s)) < 1e-9:
            P_s = P_new
            break
        P_s = P_new
    gain_norms.append(np.linalg.norm(K_s))
    eig_mags_best.append(
        np.max(np.abs(np.linalg.eigvals(A_true - B_true @ K_s)))
    )
gain_norms = np.array(gain_norms)
eig_mags_best = np.array(eig_mags_best)

Out[24]:

Visualization

Gain magnitude rises as control cost R decreases, while the worst closed-loop eigenvalue magnitude follows a nonmonotone curve below one in this experiment. — LQR gain magnitude (solid) and worst true-plant closed-loop eigenvalue magnitude (dashed) versus the control penalty R. The horizontal axis is logarithmic. Cheaper control increases the gain; the spectral-radius response is not monotone over the full sweep.

As $R$ shrinks, control becomes cheap and the gain grows. The worst eigenvalue magnitude of the true plant is nonmonotone over this sweep, so gain size alone should not be read as a universal speed measure. Every tested true-plant eigenvalue remains below the dotted unit line, but that is a numerical result for this small identification error. The Riccati guarantee applies directly to each estimated-model design; transfer to the true plant requires a robustness margin or an explicit check such as this one.

When the structure is approximately rightLink Copied

LQR's limits mirror its strengths: its exact solution needs linear dynamics and quadratic costs. Iterative LQR (iLQR) and differential dynamic programming approximate nonlinear dynamics and costs around a nominal trajectory, compute local feedback updates, roll out a new nominal, and repeat. Model predictive control (MPC) solves a finite-horizon problem at each step, applies the first action, observes the next state, and replans. That feedback can correct some disturbances and model errors before they accumulate across many horizons, but it does not erase within-horizon error or guarantee robustness. Local approximation and receding-horizon replanning form a bridge to the planning agents of Part VII.