iLQR and Differential Dynamic Programming

A boat approaching a quay must arrive in the right place, face along the berth, and lose its remaining momentum. A sequence of thrust commands determines all three outcomes. The shooting formulation turns this task into an optimization over those commands: simulate a candidate sequence, measure its cost, and improve it. How can an optimizer use the time structure of the simulation when computing that improvement?

Linearizing the dynamics along the current trajectory and quadratically approximating its costs gives a problem whose controls can be eliminated backward, one time step at a time. The resulting correction is then tested through the nonlinear simulator. Repeating these operations gives the iterative linear-quadratic regulator, or iLQR. Keeping additional curvature from the dynamics gives differential dynamic programming, or DDP. Both methods compute corrections around a nominal trajectory. Their backward recursions follow from substitution and quadratic minimization.

Learning Goals¶

After working through the derivations and examples, you should be able to:

• Construct a local quadratic optimal-control problem around a feasible rollout.

• Eliminate its controls backward and recover an affine correction at each step.

• Implement nonlinear trial rollouts with line search and regularization.

• Identify the dynamics-curvature terms that distinguish DDP from iLQR.

• Assess a docking plan using its cost, arrival velocity, and hull clearance.

• Explain how thermal load and nonlinear power loss shape a bounded cooling schedule.

Prerequisites¶

The chapter uses the finite-horizon optimal-control formulation, single shooting, and first and second derivatives. The sequential-methods discussion introduces general SQP subproblems, and the initial-value-problem appendix reviews numerical integration. Neither dynamic programming nor a previous derivation of LQR is required.

A Docking Problem¶

Docking requires the boat to reach the berth, align with the quay, and shed its momentum. A trajectory that reaches the right position with appreciable translational or angular velocity will carry the boat away again.

Consider a two-metre boat moving in the horizontal plane. Its state and control are

The position and velocity are measured in fixed world coordinates, $\psi$ is the heading, and $\omega$ is the angular velocity. Each control is a fraction of maximum thrust. Positive commands push forward and negative commands push backward. Unequal thrusts turn the boat. There is no directly commanded sideways force, so an approach with sideways momentum requires a maneuver.

Let $\mathbf e(\psi)=(\cos\psi,\sin\psi)^\top$ point along the boat and $\mathbf n(\psi)=(-\sin\psi,\cos\psi)^\top$ point across it. The simplified inertial model is

The drag opposes velocity along each boat axis, with stronger resistance to sideways motion. Thrust changes acceleration, so releasing the controls does not stop the boat immediately. These equations are a synthetic teaching model for calm water, with neither contact forces nor waves.

The quay occupies $p_y\leq0$, and the desired arrival state is $\mathbf x_\star=(0,1,0,0,0,0)$. The boat should lie parallel to the quay with its center one metre from the edge. The angled approach starts from $(-8,5,-0.5,0.2,-0.1,0)$. The sideways-drift case changes only $v_y$ to $-0.5\ \mathrm{m/s}$. Angles are in radians. Both experiments use $T-1=200$ piecewise-constant controls, each held for $h=0.1\ \mathrm s$, for a total of 20 seconds. Thus $\mathbf x_1$ is the initial state and $\mathbf x_T$ is the arrival state, as in the preceding chapters. One fourth-order Runge--Kutta step gives $\mathbf x_{t+1}=\mathbf f_t(\mathbf x_t,\mathbf u_t)$ for $t=1,\ldots,T-1$.

To discourage the hull from approaching the quay, let $y_j(\mathbf x)$ be the world-coordinate height of rectangular corner $j$. With a desired buffer $a=0.2\ \mathrm m$ and smoothing length $s=0.08\ \mathrm m$, define

This is a smooth penalty for violating the buffer. Its finite weight permits violations, so the final trajectories will also undergo a separate geometric clearance check. The full rectangle is used for that check, including the area outside the pointed boat icon in the figures.

For $\mathbf e_p=\mathbf p-(0,1)^\top$, the running and terminal costs are

These numerical weights apply to the stated SI coordinates and dimensionless thrust commands. The periodic heading penalty assigns the same cost to angles that differ by a full turn. The simulator retains an unwrapped angle; it does not insert a discontinuity at $\pm\pi$ into the differentiated transition.

The initial control sequence is zero, so the boat coasts and slows under drag. It stops far from the berth. Successive control updates can first bring it closer and then adjust the arrival heading and velocity.

Taylor Models of Dynamics and Cost¶

The local model describes the change in cost caused by a small change in the thrust sequence. Start from the finite-horizon problem

Here $\mathbf U=(\mathbf u_1,\ldots,\mathbf u_{T-1})$ is the complete control sequence. Its rollout determines every later state.

Set aside the control bounds while deriving the local correction. Roll out a nominal control sequence $\bar{\mathbf U}$ to obtain states satisfying $\bar{\mathbf x}_{t+1}=\mathbf f_t(\bar{\mathbf x}_t,\bar{\mathbf u}_t)$. For a nearby trajectory, write $\delta\mathbf x_t=\mathbf x_t-\bar{\mathbf x}_t$ and $\delta\mathbf u_t=\mathbf u_t-\bar{\mathbf u}_t$, and stack them as $\delta\mathbf z_t=(\delta\mathbf x_t,\delta\mathbf u_t)$. Taylor expansion of the discrete transition at the nominal state and control gives

Subtracting $\bar{\mathbf x}_{t+1}$ and dropping the quadratic remainder gives the linearized dynamics $\delta\mathbf x_{t+1}=A_t\delta\mathbf x_t+B_t\delta\mathbf u_t$. Both rollouts start from the same state, so $\delta\mathbf x_1=\mathbf0$. There is no constant defect because the nominal trajectory satisfies the discrete dynamics. Independently chosen state guesses would generally leave such a defect. We write $F_t=[A_t\ B_t]$ for the combined Jacobian.

Taylor expansion of the running cost to second order gives

Here $\bar c_t=c_t(\bar{\mathbf x}_t,\bar{\mathbf u}_t)$ and $c_{x,t}=\nabla_{\mathbf x}c_t(\bar{\mathbf x}_t,\bar{\mathbf u}_t)$; the other gradient and Hessian blocks are evaluated at the same nominal pair. We denote the stacked gradient by $c_{z,t}$ and the displayed block Hessian by $c_{zz,t}$. Gradients are column vectors, matching the convention of the preceding chapters. We retain the linear terms: individual stage gradients need not vanish even at an optimal trajectory because the dynamics couple stages, and at a non-optimal nominal trajectory they also drive its correction. The $c_{xu,t}$ block measures how the marginal cost of a control changes with the state. Even if the running cost has no such cross term, eliminating later controls can create one in the quadratic tail.

The terminal cost has the same expansion, with no control coordinate:

Dropping the remainders in these Taylor expansions yields a quadratic cost subject to linearized transitions. The state at one step couples only to its neighbors. This is the local problem solved by backward elimination. A general SQP method also forms local quadratic subproblems, but its exact Lagrangian Hessian includes terms from dynamics curvature. Those terms enter below when we develop DDP.

Backward Elimination¶

The last control affects only its own stage and the terminal state when its starting state is fixed. Eliminating it first leaves a shorter quadratic problem of the same form. Repeating that step works backward through the horizon.

A two-step calculation¶

Consider a scalar example with $T=3$, initial state $x_1=0$, dynamics $x_{t+1}=x_t+u_t$, and cost

Treat $x_2$ as given while eliminating $u_2$. Substitution of $x_3=x_2+u_2$ gives

Differentiating with respect to $u_2$ gives $2u_2+x_2-1=0$, so $u_2=(1-x_2)/2$. Substituting this expression back into the last two terms of the cost leaves $(x_2-1)^2/4$. Since $x_2=u_1$, the remaining problem is

Its derivative is $u_1+(u_1-1)/2$, giving $u_1=1/3$. Forward substitution then gives $x_2=1/3$, $u_2=1/3$, $x_3=2/3$, and $J=1/6$. The final position is short of one because the objective trades terminal error against effort; arrival was penalized, not imposed as an equality.

The expression $u_2=(1-x_2)/2$ also carries information that a single number would discard: if the state reaching the last step changes, the minimizing last control changes with it. This dependence becomes the matrix feedback correction in the general calculation.

Eliminating a vector control¶

Suppose all controls after time $t$ have already been eliminated from the local problem. Write the remaining quadratic tail as

At the terminal step, $S_T$ is the quadratic Taylor polynomial of $c_T$, so $s_T=c_T(\bar{\mathbf x}_T)$ and its other coefficients are the $p_T,P_T$ defined above. At earlier steps, $S_{t+1}$ is obtained by elimination from a finite optimization problem. This construction does not require a function over all states of the nonlinear system.

Substitute $\delta\mathbf x_{t+1}=A_t\delta\mathbf x_t+B_t\delta\mathbf u_t$ into this tail and add the stage cost. With $F_t=[A_t\ B_t]$, the new linear and quadratic coefficients are

We use $M_t$ for the local Hessian because $H_t$ denoted the Hamiltonian in the Pontryagin chapter. Partition $q_t$ and $M_t$ according to the state and control coordinates. Apart from a constant, the expression to minimize is

The time index is suppressed within this one-step calculation. When $M_{uu}$ is positive definite, setting the control derivative to zero gives

where

Restoring the time index gives the affine correction $\delta\mathbf u_t=k_t+K_t\delta\mathbf x_t$. The vector $k$ changes the nominal command even at zero state deviation. The matrix $K$ adjusts that change for a different state arriving from the earlier steps. In code, these formulas are linear solves with $M_{uu}$; an explicit matrix inverse is unnecessary.

Substituting the minimizing control back into the quadratic expression gives the coefficients needed by the preceding step:

The subtracted term in $P_t$ is the Schur complement associated with eliminating the control block. Repeat this operation from $T-1$ to one, starting with the terminal coefficients. A subsequent forward substitution recovers the state and control changes. For fixed state and control dimensions, the number of these elimination steps grows linearly with the horizon.

For linear dynamics and quadratic costs, this procedure solves the original unconstrained problem exactly when its control minimizations are well posed. That problem is called the linear-quadratic regulator problem, and the matrix recursion for $P_t$ is a Riccati recursion. Along a nonlinear trajectory, $A_t$, $B_t$, and the cost derivatives describe only a local approximation. They must be recomputed after the trajectory changes.

Nonlinear Rollouts and iLQR¶

The backward pass minimizes a quadratic approximation. A large correction can make its linearized state prediction inaccurate, so a proposed control sequence must be evaluated through the original discrete simulator before it is accepted.

Starting from the same initial state, form a nonlinear trial trajectory using

The step size $\alpha$ scales the feedforward change. The feedback term responds to the state deviation that actually occurs during this trial. Its gain is not multiplied by $\alpha$. At $\alpha=0$, induction through these equations recovers the nominal trajectory, because each state deviation is then zero. For an exact linear-quadratic problem, $\alpha=1$ gives the complete solution derived by elimination.

A line search first tries $\alpha=1$, then successively smaller values. Each trial has its own nonlinear states and cost. The implementation uses $1,1/2,\ldots,2^{-11}$ and accepts a strictly lower finite cost with a small sufficient-decrease check. If

the backward pass estimates a reduction $\widehat\Delta J(\alpha)=-\alpha d_1-\alpha^2d_2$. Acceptance requires $J(\bar{\mathbf U})-J(\mathbf U^+)\geq10^{-4}\max(\widehat\Delta J(\alpha),0)$ as well as a strict decrease. This estimate guides acceptance; it does not replace evaluation of the original objective.

The control curvature may be indefinite or nearly singular. Replace the matrix used in the control solve by $M_{uu}+\mu I$, where $\mu>0$. Increasing $\mu$ makes the correction more conservative and can make the solve positive definite. The code checks a Cholesky factorization, increases $\mu$ after a failed backward pass or line search, and decreases it after an accepted step. The numerical defaults start at 10^-4, multiply by ten after failure, and divide by three after acceptance, with a floor of 10^-9.

When regularization changes $k$ and $K$, the simplified Schur-complement formulas above no longer describe substitution into the original quadratic model. The implementation therefore uses the full expressions

Here $M_{uu}$ is the original curvature block; regularization was used to choose the correction. These expressions also remain applicable when some controls are fixed at their bounds.

Alternating local linear-quadratic approximation, backward elimination, and nonlinear rollout gives iLQR Li & Todorov, 2004. The iteration uses second derivatives of the costs but only first derivatives of the dynamics. It is often described as a Gauss--Newton approximation to DDP. That terminology should not obscure what is approximated: the omitted terms come from dynamics curvature. If a nonlinear cost is itself a sum of squared residuals, replacing its Hessian by a residual-Jacobian product would be a further approximation; the docking code uses the exact cost Hessian.

The returned gains describe a neighborhood of the final nominal trajectory. The docking animation below executes one stored plan. The gains used inside the optimizer do not imply that the animation replans as the boat moves.

Dynamics Curvature and DDP¶

iLQR drops the quadratic remainder in the transition’s first-order Taylor expansion. For a turning boat, a change in heading changes the direction of thrust, so products of state and control changes can affect the next state. DDP retains these second-order terms in the backward calculation.

For component $i$ of the transition, Taylor expansion one order further gives

Here $f^i_{zz,t}$ is the Hessian of the $i$th component of $\mathbf f_t$ with respect to the stacked state and control, evaluated at the nominal pair. Substitute this expression into the same quadratic tail $S_{t+1}$. Its linear term $p_{t+1}^\top\delta\mathbf x_{t+1}$ now contributes an additional quadratic term. Its quadratic term contributes $\tfrac12\delta\mathbf z_t^\top F_t^\top P_{t+1}F_t\delta\mathbf z_t$; products involving the second-order part of the transition there are of order three or higher and are discarded. Consequently,

The extra contraction weights curvature of each next-state component by its linear coefficient in the eliminated tail. At a given backward stage, this is the second-order chain rule for composing the transition with that tail. In block form, the differences from iLQR are

Setting those three sums to zero recovers the iLQR construction. Retaining them and using the same quadratic elimination and nonlinear line search gives DDP Tassa et al., 2014. The two methods share an algorithmic structure; the curvature choice can change both the direction of a step and the amount of regularization needed to accept it. This derivation does not require solving a subproblem with quadratic dynamics constraints.

For the boat, even the continuous acceleration contains a state--control cross derivative: $\partial^2\dot v_x/(\partial\psi\,\partial u_L) =-(F_{\max}/m)\sin\psi$. The discrete derivatives used by the solver also include how heading and velocity change within the RK4 step. Automatic differentiation is applied through that complete step. A differential equation that is affine in its controls need not yield a numerical transition whose second control derivatives vanish.

The extra curvature costs computation and can be indefinite. DDP therefore still needs regularization and nonlinear acceptance tests. Neither method guarantees the globally best docking maneuver. For affine dynamics the transition Hessians are zero, so the two backward constructions coincide when all other choices agree.

Thruster Limits and Stopping Criteria¶

The boat’s controls must remain in $[-1,1]^2$, even when the unconstrained correction asks for more thrust. At zero state deviation, the feedforward step is the solution of the box-constrained quadratic problem

For two controls, each coordinate is either free, fixed at its lower bound, or fixed at its upper bound. The teaching implementation checks the nine combinations, minimizes over each face, and keeps the best feasible candidate. Positive-definite curvature makes this a strictly convex QP. Larger systems can use an active-set solver rather than enumerate all combinations Tassa et al., 2014.

The feedback rows of saturated controls are zero while their active set stays fixed. For the free coordinates $\mathcal F$, solve $(M_{uu}+\mu I)_{\mathcal F\mathcal F}K_{\mathcal F} =-(M_{ux})_{\mathcal F}$. Use the full substitution formulas to update the quadratic tail. With bounds, this tail describes a neighborhood with the selected active sets; across active-set changes the minimized expression is generally piecewise quadratic. During the nonlinear trial, project the resulting control onto its box before simulation. This enforces the physical bounds when a trial leaves the neighborhood represented by the backward pass. The nonlinear cost check then decides whether that trial should be retained.

The stopping test is separate from docking success. Let $g=\nabla_{\mathbf U}J$ be the gradient of the original shooting objective, computed by the adjoint recursion along the current nonlinear rollout. With normalized controls, the solver stops when

At an interior optimum this reduces to a small control gradient. At a bound it also allows a gradient that points toward a forbidden improvement. The experiment permits at most 120 accepted updates. If regularization exceeds 10¹² without an acceptable trial, the solver retains its last accepted trajectory and reports failure to make progress. A small residual certifies approximate first-order stationarity, not a globally optimal path or a successful arrival.

Docking Trajectories and Recorded Futures¶

Do the two curvature choices produce the same maneuver from the same starting plan? Both methods receive the coasting initialization, identical costs and control limits, and the same stopping settings. The following tables report their final nonlinear costs and arrival checks. Each final control sequence is also replayed with four RK4 substeps per control interval, giving a $0.025\ \mathrm s$ validation grid.

The arrival checks below use the finer replay.

All four runs satisfy the stopping test and the docking checks: position error below $0.2\ \mathrm m$, heading error below $3^\circ$, speed below $0.05\ \mathrm{m/s}$, angular speed below $1^\circ/\mathrm s$, bounded thrust, and positive hull clearance. The maximum position difference between the planning trajectory and the finer replay is below $10^{-6}\ \mathrm m$ in these runs. Clearance is checked on the finer grid; this remains a sampled validation, not a proof of continuous-time collision avoidance.

In the angled approach, iLQR reaches cost 8.0515 in 18 accepted updates, while DDP reaches a different stationary trajectory with cost 11.9879 in 34 updates. The DDP trajectory makes two full turns near the berth. Its final unwrapped heading is approximately $720^\circ$, which the periodic heading cost treats as the same orientation as zero. Reaching this stationary trajectory does not remove the extra motion accumulated on the way there. In the sideways-drift case, both reach cost 7.6674, with 13 accepted updates for iLQR and 45 for DDP. These are outcomes of two particular local solves, not a general ranking of the algorithms. The difference in the first case also shows why comparing iteration counts alone can hide a difference in the solutions obtained.

Within one selected plan, simulation time has a different meaning from optimizer iteration. At two seconds, for example, the boat occupies the two-second state of that rollout. Its remaining states give a prediction of where it will go if the rest of the stored controls are applied. Outlines spaced equally in time spread apart while the boat is moving quickly and cluster as it comes to rest. Reverse thrust helps remove momentum before arrival.

The replay reads recorded simulations, so exploring it does not rerun an optimizer. Recomputing a plan from newly observed states is the additional step introduced in receding-horizon control.

Pulling Down a Thermoacoustic Refrigerator¶

A standing-wave thermoacoustic refrigerator has a loudspeaker at one end of a sealed tube and a stack of thin plates inside it. Oscillating gas parcels exchange heat with the plates as they compress and expand, carrying heat from the cold end of the stack toward the hot end. Heat exchangers connect those ends to a payload and to ambient air. How should the loudspeaker amplitude vary when acoustic energy and the payload temperature at a fixed deadline both matter?

We represent the cold payload and hot exchanger by two lumped temperatures, $\mathbf x_t=(T_{c,t},T_{h,t})$, and control the normalized driver amplitude $u_t\in[0,1]$. Let $\Delta T=T_h-T_c$ and $\eta(\Delta T)=1-\Delta T/\Delta T_0$. The continuous model is

Here $\dot Q_c$ is heat removed from the cold side, $\dot W$ is work supplied by the driver, and $UA$ is the hot exchanger’s conductance to ambient. At positive work, the cooling coefficient of performance (COP) is $\dot Q_c/\dot W$. The short-stack approximation motivates the $u^2$ terms and their dependence on the temperature span Swift, 1988Swift, 2017. The cubic term represents an additional loss at high amplitude. These coefficients are synthetic teaching values, not measurements of a particular device.

Both temperatures begin at $20\ ^\circ\mathrm C$. A fourth-order Runge--Kutta step of length $h=1\ \mathrm s$ defines the discrete transition $\mathbf x_{t+1}=\mathbf f_t(\mathbf x_t,u_t)$ for $T-1=300$ control intervals. The cost uses the same running-plus-terminal form as the docking problem:

There is no running temperature-tracking term: the target applies at the end of the five-minute pull-down, while the running cost measures energy use. The terminal penalty is soft, so a plan can trade some arrival error for lower energy. Both methods start from the constant sequence $u_t=0.5$ and use the same control box and stopping settings.

Costs are evaluated on the final accepted nonlinear trajectory. The iteration-limit row is a retained local plan, not a converged solve.

For comparison, constant full amplitude in the baseline case ends at 2.20 °C, uses 1426 J, and has cost 149.89.

The baseline plan drives relatively hard at first, eases off, then ramps to full amplitude near the deadline. Early driving uses the initially small temperature span, but it also warms the hot side and weakens later pumping. The fixed load adds the same total heat over every 300-second plan. Control timing still changes the temperature span and therefore the final cold temperature. Running at full amplitude throughout uses about $1426\ \mathrm J$ and reaches $2.20\ ^\circ\mathrm C$; its nonlinear cost is 149.886, versus 51.224 for the baseline local plan. The soft terminal cost need not favor the coldest possible final state.

Without the parasitic load, a nearly constant amplitude is economical. To see the tendency, temporarily hold the temperature span fixed and write $v=u^2$. Heat pumped is then proportional to $v$, while the cubic loss is proportional to $v^{3/2}$, a convex function for $v\geq0$. For a fixed total amount of pumping, Jensen’s inequality favors spreading $v$ through time. The actual schedule is not exactly constant because the span and hot-side temperature still respond to the controls. With $C_h=200\ \mathrm{J/K}$ and $UA=0.2\ \mathrm{W/K}$, the amplitude rises through most of the horizon, then turns off for the final two intervals. Both thermal storage and heat rejection changed in that experiment, so it does not isolate the effect of either parameter.

The DDP curvature term also has a direct physical interpretation here. For the continuous cold-side equation,

Increasing the temperature span weakens pumping, and this mixed derivative measures how that effect changes with amplitude. The solver differentiates the complete RK4 transition, whose Hessian also contains effects from the intermediate integration stages.

When $k_3=0$, both pumping and work are quadratic in amplitude at a fixed temperature span. Their instantaneous ratio then no longer depends on amplitude, leaving weaker preferences among some schedules. In the recorded run, iLQR reaches its iteration limit while DDP satisfies the stopping test; both reach nearly the same final cost. For the sluggish hot side, DDP reaches a different, slightly lower-cost local plan than iLQR. Lowering the target to $0\ ^\circ\mathrm C$ produces a plan saturated at $u_t=1$ throughout the fixed horizon, exercising the box-constrained backward step. These outcomes depend on the initial sequence. Starting the baseline from constant amplitudes 0.1 and 0.9 leads to other local plans, so neither result certifies a global minimum.

Exercises¶

Computational Sources¶

The implementation separates the generic backward recursion from the physical models and artifact builders:

• Shared iLQR and DDP solver

• Boat dynamics, costs, and validation

• Experiment and figure builder

• Numerical diagnostics and source hashes

• Thermoacoustic model and costs

• Thermoacoustic experiment builder

• Thermoacoustic diagnostics

Run uv run python scripts/build_boat_docking_artifacts.py from the repository root to reproduce the solves, figures, results table, and browser data. Normal book builds read those artifacts. Tests independently check the derivatives, the composed DDP curvature, the scalar calculation, and agreement between backward elimination and a dense quadratic solve.

Run uv run python scripts/build_thermoacoustic_pulldown_artifacts.py to regenerate the refrigerator results table, figures, and diagnostics.

Summary and Outlook¶

A feasible nonlinear rollout supplies the point around which dynamics and costs are approximated. Backward elimination of a local quadratic problem then supplies both a feedforward change and a state-dependent correction. iLQR repeats that calculation with linearized dynamics; DDP includes the second-order chain-rule terms from the transition. Nonlinear rollouts, regularization, and acceptance tests connect these local calculations to a decreasing sequence of actual trajectory costs. The refrigerator applies the same calculation to a bounded driver amplitude, with energy use and terminal temperature competing in the objective.

The docking and refrigerator experiments both use discrete transitions obtained by integrating ordinary differential equations. Continuous-time transcription and collocation develops other ways to represent continuous trajectories inside a finite optimization problem. Later, finite-horizon dynamic programming gives a broader interpretation of the quadratic tail: in the linear-quadratic case, the function produced by elimination is the optimal cost-to-go.