Differentiable Planning

In Sampling-Based Planning and Model Predictive Control, we planned by proposing action sequences, rolling each one through the world model, scoring the result, and keeping the best. Derivative-free population methods remain usable when model gradients are unavailable or unreliable, but their performance depends on the proposal, objective, and evaluation budget. Random shooting with a twenty-step horizon and six action dimensions searches a 120-dimensional space. If each added coordinate has conditional probability at most a fixed $q<1$ of landing in its acceptable subset under the proposal, then one sample's success probability is at most $q^d$ after $d$ such coordinates. Independent uniform coordinates with a fixed acceptable range fraction $q$ are one example. This is the relevant curse of dimensionality, not a guarantee that every high-dimensional planning problem is hard.

The derivative-free random-shooting planner used in the previous chapter does not compute or use an action-space derivative telling it which infinitesimal change lowers cost. Random shooting relies on sampled candidates, while CEM uses elite-sample statistics to concentrate its proposal on regions found promising. If the proposal assigns probability $p$ to near-good plans, $N$ independent candidates contain at least one with probability $1-(1-p)^N$. When adding coordinates makes $p$ shrink exponentially, maintaining the same success probability requires exponentially more candidates. On smooth problems with reliable gradients, ignoring those derivatives can therefore require many more model evaluations to locate a useful update. The advantage is conditional: unstable, biased, or locally misleading gradients can erase it.

Differentiable planning asks a different question. Instead of "which of these candidates is best?", it asks "in which local direction should I change this candidate?" If the rollout and cost are differentiable along the executed trajectory, reverse-mode automatic differentiation gives the entire action-gradient vector at a small constant multiple of the cost of one objective evaluation, independent of the number of action coordinates. For cost minimization, whenever $\nabla_a J\neq0$, the negative gradient is a local descent direction. With a suitable step size and constraint handling, any nonzero feasible negative-gradient or projected-gradient direction reduces cost locally and updates the action blocks with nonzero components. This direct sensitivity calculation can replace or complement population sampling.

The recipe is short enough to fit in a paragraph. Start with a guess for the action sequence. Roll the model forward from the current state, recording every intermediate state. Evaluate the cost of that trajectory. Backpropagate the cost through the rollout to obtain $\partial J / \partial a_t$ for every time step $t$. Take a gradient step on the action sequence. Repeat until a declared stopping rule or computation budget is met, such as a maximum iteration count paired with an objective-improvement tolerance and patience window. Then apply the first action and replan from the new state. Planning has become optimization.

This reframing has strict prerequisites and important failure modes:

• The action-to-objective computation must provide a usable action derivative or derivative estimator. Ordinary reverse mode needs usable local pullbacks along the executed path. Many neural latent-dynamics models for continuous control use differentiable operations on this path, while stochastic objectives can instead use a valid derivative estimator. Hard contact and friction switches can introduce kinks or jumps. Discrete mode changes can also produce branch-dependent derivatives.
• Decision-relevant objective terms on the differentiated action-to-objective path must contribute a usable action derivative or valid derivative estimator; separately handled terms need not. Quadratic tracking costs are differentiable, although their action signal can vanish or be poorly scaled. A hard indicator cost ("one if you miss the target") has a zero, and therefore uninformative, pathwise gradient away from its boundary; at the boundary it is generally nondifferentiable. A smooth surrogate can supply a pathwise gradient. For an expected cost under a suitable stochastic formulation, a score-function estimator can be valid under the regularity conditions developed later, though it often has high variance and may need variance reduction; if it is not usable, choose another valid sensitivity route. If differentiation through the hard term is bypassed, preserve its decision effect through a constraint, full-objective candidate evaluation, or an outer discrete or hybrid search.
• The gradients must remain numerically usable across the planning horizon. A long rollout contains long Jacobian products. Those products can shrink, grow, remain bounded, oscillate, or show transient amplification; strong contraction, expansion, and ill-conditioning are the failure modes that concern us here.

Each condition requires its own tools. This chapter builds the derivative machinery, examines how costs and terminal values shape the objective, and then focuses on numerical gradient usability, a practical failure mode over long horizons.

We derive backpropagation through dynamics as an adjoint recursion, translate it into action optimization, examine what makes a cost or value function usable as a planning objective, and confront the numerical realities of long rollout derivatives. Along the way we build a complete differentiable planner for a double integrator and then isolate its gradient chain under contracting and expanding scalar dynamics.

The previous chapter used $s_t$, $a_t$, $f$, $H$, $\ell$, and $\ell_f$. Here we keep $a_t$, $f$, and $H$, but write the state as $z_t$ to emphasize that it may be latent, and write the stage and terminal costs as $c_t$ and $c_H$. Thus $z_t$ corresponds to the previous chapter's $s_t$, while $c_t$ and $c_H$ correspond to $\ell$ and $\ell_f$.

We work with a deterministic, differentiable dynamics model

$$z_{t+1} = f(z_t, a_t), \qquad z_t \in \mathbb{R}^n,\; a_t \in \mathbb{R}^m,$$

where $z_t$ is the model state: it may be a fully observed physical configuration, a latent or belief state, or an observation when that observation is sufficient for the modeled transition. The action $a_t$ is continuous. The initial state $z_0$ is given and fixed. Over a horizon of $H$ transitions we accumulate the cost

$$J(a_0, \dots, a_{H-1}) = \sum_{t=0}^{H-1} c_t(z_t, a_t) + c_H(z_H),$$

where $c_t$ is the stage cost at time $t$ and $c_H$ is an optional terminal cost. The actions are the inputs to a recurrence, and $J$ is a scalar function of those inputs. This is exactly the structure of a recurrent neural network unrolled over time, except that we are differentiating with respect to the inputs rather than the weights.

To make the chain rule concrete, trace the computational paths from one action $a_t$. It has a direct path to the stage cost $c_t$ and a path through the next state $z_{t+1}=f(z_t,a_t)$ to later states and costs. The corresponding derivatives can still be zero: the cost may be locally flat, the action Jacobian may vanish, or later dynamics may erase the perturbation. The structural point is that a current action can contribute directly to the present cost and indirectly to future costs through the intervening transitions.

Define the state Jacobian and action Jacobian at time $t$ as

$$F_t = \frac{\partial f}{\partial z}\bigg|_{(z_t, a_t)} \in \mathbb{R}^{n \times n}, \qquad G_t = \frac{\partial f}{\partial a}\bigg|_{(z_t, a_t)} \in \mathbb{R}^{n \times m}.$$

The sensitivity of a future state to a past action is then a product of state Jacobians:

$$\frac{\partial z_k}{\partial a_t} = F_{k-1} F_{k-2} \cdots F_{t+1} \, G_t, \qquad k > t.$$

Read this product from right to left: the action perturbs $z_{t+1}$ through $G_t$, and each subsequent Jacobian transforms that perturbation. A subsequent Jacobian may mix or rotate directions and may shear a neighborhood of perturbations. Whether it amplifies or attenuates a particular perturbation depends on that direction and the chosen norm. If all relevant $F_j$ have operator norm at most a common $\rho<1$, the sensitivity decays at least geometrically with separation. Growth requires corresponding expansion along the propagated directions; an operator norm above one alone does not guarantee it.

Computing every action gradient by expanding its Jacobian products separately would repeat shared factors. For example, the state sensitivities $\partial z_H/\partial z_{H-1}=F_{H-1}$ and $\partial z_H/\partial z_{H-2}=F_{H-1}F_{H-2}$ share the same suffix. A naive scalar implementation that forms each horizon-length product independently performs $O(H^2)$ transition multiplications; one backward sweep performs $O(H)$. The adjoint (co-state) recursion is reverse-mode autodiff organized to reuse those suffixes.

Define the adjoint vector

$$\lambda_t = \frac{\partial J}{\partial z_t} \in \mathbb{R}^n,$$

interpreted as a column vector. The adjoint is a sensitivity of the objective to the state: it says how much the total cost would change if we could nudge $z_t$ directly. Because $z_t$ feeds into every later state and cost through the dynamics, this sensitivity already accumulates everything that happens downstream of time $t$. At the end of the horizon the adjoint is just the terminal cost gradient:

$$\lambda_H = \frac{\partial c_H}{\partial z_H}.$$

Now walk backwards. The state $z_t$ influences the objective in exactly two ways: through the stage cost $c_t(z_t, a_t)$, and through the next state $z_{t+1} = f(z_t, a_t)$, whose adjoint we already have. The chain rule gives

$$\lambda_t = \frac{\partial c_t}{\partial z_t} + F_t^\top \lambda_{t+1}, \qquad t = H-1, H-2, \dots, 0.$$

The transpose appears because $z_{t+1}=f(z_t,a_t)$ depends on $z_t$ through $F_t$, and $\lambda$ is a column vector. Under the local linearization, $\delta z_{t+1}=F_t\delta z_t+o(\lVert\delta z_t\rVert)$. Pulling a gradient back through the transition therefore multiplies by $F_t^\top$, converting next-state sensitivity into current-state sensitivity. For this discrete recurrence, the recursion is exactly the BPTT backward pass. It is the discrete-time analogue of the continuous ODE adjoint, which instead satisfies a backward differential equation; differentiating a chosen ODE discretization produces a discrete adjoint whose details depend on that discretization.

The action gradient then follows immediately. The action $a_t$ enters the objective only through $c_t(z_t, a_t)$ and through $z_{t+1} = f(z_t, a_t)$:

$$\boxed{\;\frac{\partial J}{\partial a_t} = \frac{\partial c_t}{\partial a_t} + G_t^\top \lambda_{t+1}.\;}$$

where:

• $\partial c_t / \partial a_t$: the direct sensitivity of the stage cost at time $t$ to the action taken at time $t$
• $G_t = \partial f / \partial a \big|_{(z_t, a_t)}$: the action Jacobian, mapping an action perturbation into a next-state perturbation
• $\lambda_{t+1} = \partial J / \partial z_{t+1}$: the adjoint of the next state, which aggregates the sensitivity of every future cost

One backward sweep computes $\lambda_H,\lambda_{H-1},\dots,\lambda_1$ and every action gradient. Its arithmetic cost is a small constant multiple of the forward evaluation. A standard uncheckpointed implementation uses $O(H)$ activation memory when per-step size is fixed; it retains the local information needed for the pullbacks, not necessarily explicit Jacobian matrices. Autodiff computes vector-Jacobian products such as $F_t^\top\lambda_{t+1}$ and $G_t^\top\lambda_{t+1}$ directly. Checkpointing or recomputation can reduce retained activations at additional compute cost.

The next two cells compute and plot analytic scalar multiplication counts for two organizations of the same action-sensitivity calculation. They do not benchmark Python exponentiation. Independently forming $\rho^{H-1-t}$ by repeated multiplication uses $\sum_{r=0}^{H-1}r=H(H-1)/2$ multiplications by $\rho$. With the terminal seed $\lambda_H$ available, the adjoint sweep needs $H-1$ multiplications to produce $\lambda_{H-1},\dots,\lambda_1$; together, $\lambda_1,\dots,\lambda_H$ are the adjoints required by all $H$ action gradients.

This is the cheap gradient principle: reverse mode computes a full $H\times m$ action gradient at a small constant multiple of one scalar objective evaluation, whereas coordinate-wise finite differences require one or more extra evaluations per action coordinate. This evaluation-cost asymmetry is one important reason differentiable planning is attractive in high dimensions. Directional and simultaneous-perturbation estimators have different tradeoffs, so the comparison here is specifically with coordinate-wise finite differences.

Substituting the adjoint recursion into itself gives the explicit form

$$\lambda_t = \sum_{k=t}^{H-1} \Big(\prod_{j=t}^{k-1} F_j^\top\Big) \frac{\partial c_k}{\partial z_k} \;+\; \Big(\prod_{j=t}^{H-1} F_j^\top\Big)\frac{\partial c_H}{\partial z_H},$$

and therefore

$$\frac{\partial J}{\partial a_t} = \frac{\partial c_t}{\partial a_t} + G_t^\top \Bigg[\sum_{k=t+1}^{H-1} \Big(\prod_{j=t+1}^{k-1} F_j^\top\Big)\frac{\partial c_k}{\partial z_k} + \Big(\prod_{j=t+1}^{H-1} F_j^\top\Big)\frac{\partial c_H}{\partial z_H}\Bigg].$$

Every indirect future-cost contribution contains an action Jacobian and a possibly empty product of state-Jacobian transposes applied to a cost gradient. This formula motivates the stability analysis below, and it explains a subtlety that is easy to miss: the gradient with respect to an early action is a sum of many terms, one for each future cost, each weighted by a generally different, time-indexed Jacobian product. Terms can reinforce or cancel. A vanishing product for one future step does not necessarily mean the total gradient vanishes. Conversely, a single large product term can dominate the sum even when the surrounding terms are small, which is one reason the gradient profile across actions can look erratic in practice.

Ordinary reverse mode needs a valid local pullback for every operation on the executed path. At a switch, an autodiff implementation may select a branch or subgradient convention; a planner may instead install a surrogate or custom pullback; and some operations provide no useful derivative. Backpropagation-based planners therefore need a usable sensitivity route for the optimized path: ordinary pullbacks, event-aware or implicit sensitivities, a relaxation or surrogate, or an estimator. Otherwise they must bypass or search the hard component separately.

• Contact and collision. Rigid-contact dynamics are often smooth within a fixed mode but hybrid across impacts. Classical derivatives can jump or fail to exist when contact timing or the active mode changes, while event-aware sensitivities can be valid away from grazing and other degenerate events. Compliant differentiable contact laws or relaxed complementarity conditions can smooth selected transitions, but they change the physical model and do not automatically smooth friction or contact activation. Another local approximation differentiates the active branch while holding the contact set fixed. This is locally exact only while perturbations preserve the full contact and friction mode; across a mode change it omits event sensitivity, although an optimizer may still use it heuristically.
• Discrete actions and argmax operations. With a one-hot encoding, hard $a=\arg\max_i q_i$ has zero Jacobian away from ties and is nondifferentiable at ties, so it supplies no informative ordinary pathwise gradient for changing the selected category. A temperature-softmax relaxes a deterministic argmax; for categorical sampling, Gumbel-softmax relaxes the one-hot Gumbel-max sample. A straight-through estimator keeps the hard forward value but substitutes a surrogate backward derivative. These introduce approximation choices. For background on discrete and latent action representations, see Actions, Inverse Dynamics, and Latent Actions.
• Stochastic dynamics. A sample $z_{t+1} \sim p_\theta(z_{t+1} \mid z_t, a_t)$ does not automatically provide a pathwise derivative with respect to the action. For a reparameterizable continuous distribution, write the sample as a deterministic function of the action and independent noise, for example $z_{t+1} = \mu_\theta(z_t, a_t) + \sigma_\theta(z_t, a_t) \odot \epsilon$ with $\epsilon\sim\mathcal N(0,I)$. Holding the sampled noise fixed then gives a pathwise gradient for a Monte Carlo objective. When reparameterization is unavailable, differentiating $\mathbb E_{z\sim p_a}[L(a,z)]$ yields the ordinary term $\mathbb E[\partial_a L]$ plus the likelihood-ratio term $\mathbb E[L\,\partial_a\log p_a(z)]$, provided differentiation can be interchanged with integration and the required terms are integrable. The likelihood-ratio term also requires action-independent support or vanishing boundary contributions; nonzero boundary terms from action-dependent support need a separate valid estimator. If the score itself is estimated, its conditional expectation given $z$ must equal the true score, or the weighted term must otherwise be jointly unbiased. Score-function estimators often have high variance. See Stochasticity and Uncertainty for the probabilistic background.
• Rendering and perception. At fixed inputs and parameters, a differentiable deterministic image encoder defines a repeatable mathematical gradient, although it may be expensive or badly scaled. Implementation-level repeatability additionally requires deterministic forward and backward kernels in a fixed software and hardware environment. Stochastic rendering, noisy observations, and sampled layers can add gradient variance. Surrogate derivatives can add bias or mismatch; randomized surrogates can also add variance. Hard thresholding, such as a binarized mask or hard attention selection, needs a relaxation, estimator, or custom pullback.

Simulation systems such as DiffTaichi and the published Brax physics engine provide automatic differentiation for substantial parts of a physics computation. So do MJX-JAX and supported NVIDIA Warp kernels. Coverage is backend- and feature-dependent: current MJX-Warp is not differentiable, and contact smoothness depends on the chosen formulation and settings. Planning by backpropagation may instead use analytical dynamics, a learned model, event-aware sensitivities, or an implicitly differentiated solver.

Once we can compute $\nabla_a J \in \mathbb{R}^{H \times m}$, action-sequence planning is a constrained finite-dimensional optimization problem:

$$\min_{a_0, \dots, a_{H-1}} \; J(a) \quad \text{subject to} \quad a_t \in \mathcal{A} \;\; \forall t,$$

where $\mathcal{A}$ is the admissible action set. Independent limits form a box: with shared scalar limits it is $[a_{\min},a_{\max}]^m$, and in general it is $\prod_i[a_{\min,i},a_{\max,i}]$. Coupled actuator limits need non-box sets; rate constraints need sequence-coupled constraints or state augmentation; state-dependent limits need state-action or path constraints such as $a_t\in\mathcal{A}(z_t)$. With continuous decision variables, the problem is a nonlinear program when its objective or explicit constraints are nonlinear. Discrete decisions instead lead to a discrete or mixed-integer formulation. The linear-quadratic example below is a convex quadratic program.

The simplest solver is gradient descent with a fixed scalar step: $a\leftarrow a-\eta\nabla_aJ$. Adaptive methods such as Adam are convenient in implementations like the worked example below. Away from its numerical stabilizer and interventions such as clipping or weight decay, Adam's moment normalization makes updates less sensitive to a uniform positive rescaling of persistent gradients. It is not fully scale-free, and it does not remove poor conditioning. For a smooth positive-curvature local model, one global step size must be stable for the direction of largest curvature, which can force slow progress in flatter positive-curvature directions. Adaptive, accelerated, and preconditioned first-order methods may reduce that penalty but do not eliminate it in general.

Two simple approaches are especially easy to implement.

Projected gradient descent takes a step and then clips each coordinate back into the box:

$$a \leftarrow \mathrm{clip}\big(a - \eta \nabla_a J,\; a_{\min},\; a_{\max}\big).$$

This enforces the box constraint exactly. At a boundary, projection removes only the component of a proposed step that points outside the box; a later descent step that points inward can still move the coordinate. Repeated outward steps can pin a coordinate at the boundary, while changing or noisy directions can make iterates oscillate along active faces. The projection operator also introduces a kink that must be handled in convergence analysis near active constraints.

Reparameterization removes the constraint entirely. Write

$$a_t = \frac{a_{\min} + a_{\max}}{2} + \frac{a_{\max} - a_{\min}}{2} \tanh(u_t)$$

and optimize the unconstrained variable $u$. In exact real arithmetic, every finite $u$ maps into the strict interior of the box. In floating-point arithmetic, $\tanh(u)$ can round to $\pm1$, and midpoint or half-width arithmetic can cancel, overflow, or underflow. The affine formula alone therefore does not universally guarantee even closed-box feasibility. With well-scaled bounds and overflow-free arithmetic, it typically lands in the closed box; apply a final clip when exact numerical box feasibility is required. The other cost is saturation: when $|u_t|$ is large, $\mathrm{d}a_t/\mathrm{d}u_t$ becomes exponentially small. The action lies very near a boundary and may become effectively stuck in finite precision, which is harmless only while the desired action remains there and harmful when an inward recovery is needed. A logistic parameterization has the same limitation.

Both methods are used in simple first-order planners. Reparameterization guarantees strict-interior feasibility in exact arithmetic but can compound vanishing gradients when it saturates. In finite precision, overflow-safe bound arithmetic followed by a final clamp, or another interval-preserving implementation, is needed to guarantee a result in the closed box. Projection can represent an exact boundary optimum directly, whereas real-arithmetic $\tanh$ approaches a boundary only asymptotically. Overall convergence still depends on scaling, step size, curvature, and active constraints.

The trajectory transcription changes both sparsity and conditioning. A key distinction is which trajectory quantities become free variables and how the dynamics are enforced.

• Single shooting optimizes only the actions and eliminates states by forward simulation. Each action-gradient block contains direct terms plus sums weighted by products of later transition Jacobians. This is the action-only formulation established above.
• Multiple shooting also promotes segment-boundary states to decision variables and adds continuity constraints. This exposes sparse local continuity Jacobians to a constrained solver and often avoids the worst full-horizon sensitivity products, although the resulting KKT system can still be ill-conditioned. It is a standard option for difficult trajectory-optimization and nonlinear-MPC problems.
• Direct collocation represents states and controls at mesh or collocation nodes and enforces integration defects. A-stable or L-stable implicit collocation schemes, such as appropriate Radau formulas, can be more stable for stiff modes than explicit integration at the same step size, although their nonlinear constraint Jacobians can still be poorly conditioned.

For sufficiently short, well-conditioned receding-horizon problems with modest path constraints, single shooting is a practical baseline with less transcription bookkeeping. As horizon-induced sensitivity products become problematic, multiple shooting can justify its additional variables and constraints; the crossover depends on dynamics, time step, dimensions, constraints, scaling, and solver. Its local continuity blocks shorten sensitivity chains but are not guaranteed to be well-conditioned.

Gradient descent uses only first-order information, and planning objectives can be badly conditioned when long-horizon sensitivities create strongly unequal curvature across action directions. Curvature-aware methods attack that local geometry directly.

Differential Dynamic Programming (DDP) and its Gauss-Newton variant iLQR compute a local quadratic approximation of the cost-to-go. Around a nominal pair $(\bar z_t,\bar a_t)$, define $\delta z=z-\bar z_t$ and $\delta a=a-\bar a_t$; let $q_z,q_a$ denote the corresponding gradient blocks of the local action cost-to-go $\mathcal Q_t$, and let $Q_{zz},Q_{za},Q_{az},Q_{aa}$ denote its Hessian blocks. When the unconstrained local control subproblem is well posed, typically after regularization makes $Q_{aa}$ positive definite, the backward algebra solves that local model exactly; it only approximates the original nonlinear problem. The backward pass augments the first-order adjoint, or value gradient, with a value Hessian:

$$\mathcal Q_t(z, a) \approx \text{const} + \tfrac{1}{2}\begin{bmatrix} \delta z \\ \delta a \end{bmatrix}^\top \begin{bmatrix} Q_{zz} & Q_{za} \\ Q_{az} & Q_{aa} \end{bmatrix} \begin{bmatrix} \delta z \\ \delta a \end{bmatrix} + \begin{bmatrix} q_z \\ q_a \end{bmatrix}^\top \begin{bmatrix} \delta z \\ \delta a \end{bmatrix},$$

and, when $Q_{aa}$ is positive definite, the local feedforward perturbation is $\delta a^\star=-Q_{aa}^{-1}q_a$. At each time step the sweep forms $Q_{aa}$, $Q_{az}$, and the local linear terms from the stage-cost derivatives, $F_t$, $G_t$, and the next value gradient and Hessian, then propagates the resulting value gradient and Hessian. Full DDP also retains contractions with second derivatives of the dynamics, while iLQR drops those dynamics terms but keeps the first and second cost derivatives. The resulting finite-horizon, generally time-varying Riccati recursion is related to, but not identical to, the discrete algebraic Riccati equation used for the infinite-horizon linear time-invariant example below.

On smooth problems with a useful local quadratic model, DDP or iLQR can converge in fewer iterations than first-order descent. The advantage is not guaranteed, especially when $Q_{aa}$ is indefinite or poorly stabilized rollouts make the local model inaccurate. DDP retains dynamics second derivatives, while iLQR drops them; both typically need regularization and a line search or trust region. With a descent-producing regularization, backtracking can reject an aggressive step until an acceptance condition is met. It can still stall or fail when no useful local direction exists.

In a smoothly evolving MPC loop, consecutive problems share structure and usually differ modestly unless disturbances, references, constraints, or modes change abruptly. Under those conditions, shifting the unused action suffix left, appending a terminal guess, and checking or restoring feasibility for the new state and constraints often gives a useful full-horizon initial sequence. The appended value can be the previous last action or zero only when it satisfies the new box, rate, coupled, and state-dependent constraints. This warm start alone does not guarantee recursive feasibility; that requires additional MPC assumptions. Warm starts can reduce the iterations needed per update, but whether this reaches a real-time deadline depends on the model, hardware, control rate, and accuracy target. The same shift-and-pad pattern is a common warm-starting discipline for CEM-based MPC, as introduced in Sampling-Based Planning and Model Predictive Control.

A useful hybrid is to evaluate a population of CEM candidates for a few distribution-update rounds to find a promising basin, then refine the best sequence with gradients. The sampling phase explores without requiring differentiability; the gradient phase performs local refinement. The outcome still depends on population size, initialization, and the objective's basin structure.

There is a second, distinct meaning of "differentiable planning" to separate from the first.

So far we have used gradients to perform the planning. We can also make the planning procedure a differentiable layer so gradients flow into upstream cost or dynamics parameters. One route solves the optimization problem in the forward pass and differentiates its KKT conditions implicitly. A standard sufficient regularity regime for classical implicit differentiation uses a differentiable local KKT residual, a strongly regular or locally unique solution branch, and a nonsingular linearized optimality system. A locally stable active set commonly permits a fixed-active-set derivative only when the problem data are differentiable and the corresponding reduced KKT system is nonsingular, conditions typically supported by constraint qualification and suitable second-order assumptions. Active-set changes may make the solution map nondifferentiable, but they do not do so in every problem. These conditions enable the classical inverse-system sensitivity formula; they are not universal necessities for every differentiable selected optimizer map. Unrolling a fixed number of solver iterations is a different route with different memory and approximation tradeoffs.

This distinction matters in learning-to-plan systems and differentiable optimization layers. It is separate from the next part of this chapter: using differentiable values inside a planning objective. The planned Policies Learned in Imagination chapter will instead examine policy learning from model-generated rollouts.

The planning objective is a consequential design choice because it defines what the optimizer is asked to prefer. With accurate gradients and a suitable step-size or line-search rule, a gradient-based planner follows local descent directions in that objective, even if it has a convenient minimum in the wrong place.

A standard stage cost for smooth continuous control and trajectory tracking is

$$c_t(z_t, a_t) = (z_t - z^\star)^\top Q (z_t - z^\star) + a_t^\top R a_t,$$

with $z^\star$ the target state, $Q\succeq0$ weighting state error, and $R\succ0$ penalizing effort. For physical coordinates with documented units and suitably structured, calibrated $Q$ and $R$, the terms can be interpreted as weighted state-error and action-effort contributions. For a learned latent state, the first term is a $Q$-weighted squared discrepancy—a squared seminorm if $Q$ is singular—and need not align with task semantics. Without a bound on decoder sensitivity, numerically close latents can still decode to substantially different observations. For a decoder $d(z_t)$ and target $x^\star$ in $\mathbb R^p$, a decoded cost $(d(z_t)-x^\star)^\top Q_x(d(z_t)-x^\star)$ with $Q_x\succeq0$ in $\mathbb R^{p\times p}$ measures error in observation coordinates and costs a decoder evaluation per step. It is reliable as a task-semantic cost only if the observation metric is appropriate and decoder error over planned states does not materially distort task-relevant costs, rankings, or gradients. See Sufficient State and State Abstraction for the representation question.

Other common forms:

• Smooth surrogates for sparse costs. Away from its boundary, a binary goal term has zero pathwise gradient and gives an action optimizer no local signal from that term. For $\beta>0$, a smooth alternative is $c(z)=\beta^{-1}\log\!\left(1+\exp\left(\beta(\lVert z-z^\star\rVert^2-\epsilon^2)\right)\right)$, a softplus approximation to the hinge violation $\max(0,\lVert z-z^\star\rVert^2-\epsilon^2)$ rather than to a binary indicator. A Gaussian goal bonus, $c(z)=-\exp(-\lVert z-z^\star\rVert^2/2\sigma^2)$, is another option. At fixed distance $r=\lVert z-z^\star\rVert>0$, its state-gradient magnitude is $r\sigma^{-2}\exp(-r^2/2\sigma^2)$ and is largest at $\sigma=r/\sqrt{2}$. In the tail $r>\sqrt{2}\sigma$, shrinking $\sigma$ further weakens the gradient; for $\sigma>r/\sqrt{2}$, shrinking it can strengthen the gradient. Tune $\sigma$ to the candidate-distance distribution, and narrow it as candidates approach only while preserving a useful signal.
• Logarithmic barriers. For a barrier weight $\kappa>0$, an inequality $g(z)>0$ on a strictly feasible trajectory can be represented by $c(z)=-\kappa\log g(z)$. The term inherits $g$'s regularity in the strict interior: it is smooth there when $g$ is smooth, and merely differentiable there when $g$ is merely differentiable. It diverges to $+\infty$ as $g(z)\to0^+$. At any finite positive $\kappa$ it is an interior approximation, not an exact replacement for a general hard constraint, and it cannot be evaluated at an infeasible point.
• Penalty terms. For a penalty weight $\rho>0$, softer alternatives add $\rho \max(0, -g(z))^2$. Penalties allow small violations in exchange for lower cost elsewhere; barriers do not.
• Learned costs. A neural cost can be learned from demonstrations through inverse optimal control or from explicit preferences through reward modeling. Smooth operations make it differentiable everywhere; standard ReLU networks provide gradients almost everywhere. Approximation and generalization errors, especially off distribution, can become exploitable when they make reachable plans look artificially attractive.
• Energy-based and contrastive costs. If a model provides a differentiable energy whose low-energy configurations align with task success, that energy can serve as a planning cost. Predictive compatibility alone does not guarantee control utility; see the chapter on predictive embeddings and energy-based models.

Cost shaping requires caution. A generic dense auxiliary term can change the optimal plan. For this finite-horizon, undiscounted cost convention, one invariant construction is

$$c'_t=c_t+\Phi_t(z_t)-\Phi_{t+1}(z_{t+1}),\qquad c'_H=c_H+\Phi_H(z_H).$$

The added terms telescope to the fixed constant $\Phi_0(z_0)$, so they do not change the ordering of feasible action sequences. Omitting the terminal correction leaves the residual $-\Phi_H(z_H)$, which is generally action-dependent and can change the optimum when $\Phi_H(z_H(a))$ varies across feasible plans. Discounted reward shaping uses a different sign and includes the discount factor; the reinforcement-learning result does not transfer verbatim without matching those conventions.

When a finite horizon approximates a longer or infinite-horizon objective, it omits a tail whose value can vary across terminal states. A terminal cost can approximate that continuation:

$$J_H^V(a) = \sum_{t=0}^{H-1} c_t(z_t, a_t) + V(z_H), \qquad V(z) \approx \inf_{\pi} \; \mathbb{E}_{\pi}\Big[\sum_{k=0}^{\infty} c(z_k, a_k) \;\Big|\; z_0 = z\Big],$$

assuming the continuation cost is finite and uses the same dynamics, constraints, costs, and information pattern as the planner. If $V$ is the exact optimal continuation cost, the subproblems attain their minima, and $H\geq1$, dynamic programming makes every minimizing prefix begin with an action in the same infinite-horizon optimal first-action set. Without a unique optimum or consistent tie-breaking, different horizons may select different actions from that set. In MPC, suitable terminal costs—often together with terminal sets and other conditions—can support stability arguments. Structurally, the terminal value is the planning analogue of an $n$-step temporal-difference bootstrap: both replace an omitted tail with a value estimate, but a planning terminal value guides action optimization only when it distinguishes reachable terminal states, while the bootstrap supplies a learning target.

For an equilibrium-centered, unconstrained infinite-horizon LTI linear-quadratic regulator, under the usual stabilizability and detectability assumptions, the value has the quadratic form

$$V(z) = (z - z^\star)^\top P (z - z^\star),$$

where $P$ solves the discrete algebraic Riccati equation

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

The target and input must be shifted to a feasible equilibrium; general affine tracking can add linear and constant terms. We will solve this DARE numerically for the exactly specified LQ model and test the resulting terminal value in a short-horizon planner.

Outside special cases such as LQR, a nonlinear problem's value function need not have this quadratic form and may lack a tractable closed form. It can be approximated, for example by a learned $V_\theta$. When it preserves decision-relevant value differences or gradients over candidate terminal states, it can support a shorter explicit rollout. Shortening that rollout reduces the maximum transition-Jacobian chain length and may mitigate horizon-induced vanishing or exploding sensitivities, although the terminal value can introduce its own conditioning problems. The model and value become coupled when joint training links them through shared parameters, representations, losses, gradients, or targets. Coupling also arises when value targets depend on model rollouts. An optimizer can push candidate terminal states outside the value function's training support. Bellman-style training may impose consistency constraints, but it does not remove approximation error. A learned terminal value therefore trades part of the omitted-tail error for value-approximation error while the shorter action problem still has to be optimized.

Differentiable value functions support at least three planning patterns.

As a terminal cost. Backpropagation starts from $\nabla V_\theta(z_H)$. That terminal-state adjoint propagates through the rollout and contributes a gradient to every earlier action.

For direct action optimization. If $Q_\theta(z,a)$ denotes a reward or return to maximize and is differentiable in a continuous action, one ascent iteration is

$$a^{(k+1)} = a^{(k)} + \eta \nabla_a Q_\theta\!\left(z,a^{(k)}\right),$$

with appropriate constraint handling for constrained actions, such as projection or reparameterization for simple action bounds. This $Q$ uses a reward convention, unlike the continuation cost $V$ above; for an action cost-to-go, the corresponding update uses the negative gradient. Each iteration differentiates one $Q$ evaluation rather than an explicit model rollout. Finding the global greedy action still requires solving a generally non-convex problem and may need restarts; differentiability alone guarantees neither a global maximizer nor even progress under an arbitrary step. For the displayed one-step update, only errors in the action gradient at the current iterate alter the step. Errors in action rankings can still affect candidate selection or the eventual greedy action.

As a tail estimate after finite lookahead. A terminal value can make short rollouts effective by estimating what lies beyond them. TD-MPC scores short latent rollouts with predicted rewards and a learned terminal value during sampling-based trajectory optimization, whereas MuZero backs up a value from a search-tree leaf. The planned chapters on MuZero and Value-Equivalent Models and TD-MPC and Control-Centric Representations will develop those distinct uses.

An imperfect learned component can create an exploitation surface when its error varies across planner-reachable candidates and changes objective values or derivatives. A vanilla optimizer sees the supplied values and derivatives, not an independent test of whether the surrogate remains valid. When a reachable model error creates a smooth direction toward lower predicted cost, gradient optimization can exploit it rapidly. It is not guaranteed to find every error, nor is it uniformly faster than sampling: local basins, flat regions, constraints, and initialization still matter.

Concretely, this shows up as:

• Plans that drive the latent state to a region where a learned cost is artificially low, decoding to something nonsensical.
• Terminal cost-to-go estimates that are too low in unexplored regions, making unsupported terminal states look artificially attractive.
• Reachable dynamics-model errors that make impossible trajectories look favorable. Directed gradients can refine such an exploit quickly once it lies in the current basin, although the comparison with sampling depends on the error landscape and budget. This is the model-exploitation failure introduced in Sampling-Based Planning and Model Predictive Control.

Several defenses against model exploitation have analogues for learned costs and values, but their uncertainty or support measures must be defined for the component they protect:

• Uncertainty-aware objectives. An ensemble penalty can discourage plans in regions of high disagreement. When disagreement is a validated proxy for decision-relevant epistemic uncertainty, high-disagreement regions can also carry elevated exploitation risk. Agreement can still hide shared systematic bias. PETS uses probabilistic ensembles and particles to propagate uncertainty during trajectory evaluation rather than adding this exact penalty; the planned chapter on probabilistic and ensemble model-based RL will compare these mechanisms.
• Smoothness constraints. Penalizing cost-gradient magnitude over a training distribution can discourage steep local variation there. It neither guarantees the removal of spurious minima nor, by itself, guarantees control of gradients at planner-reachable states outside that distribution.
• Conservative objectives. Some conservative offline-RL objectives, such as COMBO-style value regularizers, penalize Q-values under distributions chosen to emphasize out-of-support model-rollout tuples, aiming for conservative estimates there. Other conservative methods act through reward penalties or policy constraints.
• Trust regions in action space. A conventional trust region bounds each action-sequence update around the current optimizer iterate. To bound cumulative drift from an initially trusted plan, impose a constraint relative to that fixed reference or an equivalent episode-wide deviation budget.
• Evaluation on the real system. No amount of surrogate engineering substitutes for checking whether the plan works on the real system.

A useful mental model is that ordinary optimization against a fixed learned surrogate can behave like an adversarial search for errors inside the optimizer's basin, even though it is not a formal min-max problem. With other factors controlled, gradients that align better with a locally valid descent direction and have usable scale and variance can make local surrogate optimization more effective. That can improve task performance or expose smooth model errors more aggressively. Optimization quality and real-system performance must therefore be evaluated separately.

We now turn to the central technical difficulty. Everything so far assumed that $\nabla_aJ$ is usable. Over long horizons it can fail when accumulated sensitivities are strongly contractive, expansive, ill-conditioned, discontinuous, or noisy. A long horizon alone does not force failure.

Return to the sensitivity formula

$$\frac{\partial z_k}{\partial a_t} = F_{k-1} F_{k-2} \cdots F_{t+1} \, G_t.$$

If the state Jacobians were constant, $F_j = F$ for all $j$, this product would be $F^{k-t-1} G_t$. For a scalar $F=\rho$, or for a normal matrix whose relevant modes are excited, its asymptotic scale follows the eigenvalue magnitudes:

• If $|\rho|<1$, the product decays geometrically. The influence of early actions on late states vanishes.
• If $|\rho|>1$ and the action excites the scalar mode with a nonzero gain, the sensitivity magnitude grows geometrically. The influence of early actions on late states then becomes unbounded as the separation grows.
• If $|\rho|\approx1$, the modal magnitude changes little per transition. Over $p$ transitions it remains near one only when $p|\log|\rho||$ is small, and negative $\rho$ alternates sign.

For a general matrix, spectral radius alone is insufficient over a finite horizon: non-normal matrices can exhibit large transient growth even when every eigenvalue lies inside the unit disk. Time-varying nonlinear rollouts require the singular values and directions of the actual product $F_{k-1}\cdots F_{t+1}$. The scalar spectral-radius example is therefore intuition, not a complete stability test. A real system's rates come from its physical dynamics. The implemented model's state and action Jacobians supply the planner's rollout-sensitivity factors and match the real rates only to the extent that the model is locally faithful. The full objective gradient also depends on cost or value derivatives and on the coordinates used to parameterize actions. Constraints separately determine feasible or projected directions unless they are encoded through a barrier, penalty, or reparameterization. For nonlinear dynamics, Jacobians depend on the evaluation point and can vary along a rollout.

Numerical overflow and underflow. Float32's largest finite magnitude is about $3.4\times10^{38}$, its smallest positive normal value is about $1.2\times10^{-38}$, and subnormal values extend to about $1.4\times10^{-45}$. With $\rho = 1.2$ and $H = 100$, $\rho^{H} \approx 8.3 \times 10^7$; with $\rho = 1.5$, the product is about $4.1\times10^{17}$ at $H=100$ and $1.7\times10^{35}$ at $H=200$, still finite, but it overflows near $H=219$. In the other direction, $0.9^{400}\approx5.0\times10^{-19}$ is still normal; it reaches the normal/subnormal boundary only near $H=829$. These scalar thresholds do not predict the conditioning of a matrix product, but they show when finite-precision range becomes an additional failure mode. Float64 postpones overflow and underflow without fixing vanishing or exploding sensitivity.

Ill-conditioning. Ill-conditioning can remain severe even when overflow and underflow are absent because it concerns relative directional scales rather than representable range. For a symmetric positive-definite quadratic Hessian whose smallest and largest eigenvalues are $\mu$ and $L$, respectively, $\kappa_2=L/\mu$; if $[\mu,L]$ is only an enclosing interval, then $\kappa_2\leq L/\mu$. Gradient descent with step $1/L$ contracts Euclidean iterate error by at most $1-1/\kappa_2$ per step when $L$ is the largest eigenvalue; a constant-factor reduction can therefore take on the order of $\kappa_2$ steps. In linear-quadratic problems the action Hessian is built from rollout Jacobians and cost curvature. Exact nonlinear Hessians also contain second derivatives of the dynamics and costs. This spread motivates regularized curvature-aware methods such as DDP, but they do not eliminate conditioning problems in general.

Unequal update scales. Consider the terminal-only objective for $z_{t+1}=0.85z_t+a_t$ with $H=40$. The final-state sensitivities are $\partial z_{40}/\partial a_0=0.85^{39}\approx0.00177$ and $\partial z_{40}/\partial a_{39}=1$. Away from a stationary terminal cost, the common terminal derivative cancels in their ratio, so the late-action gradient is about $0.85^{-39}\approx566$ times larger. With plain gradient descent using one global learning rate, early-action updates are proportionally smaller. Global-norm clipping rescales the entire vector and preserves this ratio; it bounds magnitude but does not repair the imbalance. Per-coordinate moment normalization can make persistent nonzero components more similar in update magnitude, except where optimizer history or $\epsilon$ dominates, and may over-amplify directions with little physical influence.

In a strongly contracting system, weak early-action influence on the terminal state can be physically correct. Early actions may still matter for stage costs. At a given candidate rollout, relative action-gradient weights depend on the initial state and visited state-action points through the dynamics, costs, horizon, discounting, and parameterization. Those weights may or may not suit the practical objective. If contraction is a model artifact, the sensitivity is misleading. If it is physical, weak terminal sensitivity is faithful, although its optimization consequences still depend on the rest of the formulation.

Gradient variance under stochastic dynamics. Under regularity conditions that permit differentiation under the expectation, reparameterization gives $\nabla_a\mathbb{E}_\epsilon[L(a,\epsilon)]=\mathbb{E}_\epsilon[\nabla_aL(a,\epsilon)]$. In practice, we estimate it by averaging pathwise gradients from sampled noise realizations. Jacobian chains can amplify or attenuate each sample's influence, so variance may grow, shrink, or plateau with horizon. When measured Monte Carlo variance is large relative to the mean gradient, use more trajectories or validated control variates. For paired comparisons, common random numbers can help when the induced covariance is favorable.

The mitigations fall into five families.

(A) Shorten the effective horizon.

• Receding-horizon MPC with a finite, task-appropriate planning horizon, the central deployment pattern in Sampling-Based Planning and Model Predictive Control.
• A sufficiently decision-faithful terminal value, which estimates the omitted tail without another long transition chain.
• Truncated backpropagation at selected recurrent-state boundaries. One scheme detaches the state $K$ transitions before the end and differentiates only the final window; this saves work but generally biases the full planning gradient unless the omitted chain-rule contributions are zero.
• Discounting with $0\leq\gamma<1$, which geometrically downweights late terms and shortens the effective horizon. It changes the objective and does not guarantee a better condition number.

(B) Control the gradient scale.

• Global-norm clipping: rescale the whole vector when its norm exceeds a calibrated threshold. This preserves direction, unlike coordinate-wise clipping, and is preferable when direction preservation is the goal. It does not correct component imbalance.
• Cost normalization: choose a practical scale from representative loss and gradient measurements, then tune the learning rate, optimizer $\epsilon$, and clipping threshold together. An order-one norm is a convention, not a requirement.
• Higher precision when diagnostics show overflow, underflow, or harmful roundoff. Float64 extends numerical range without fixing true sensitivity decay or conditioning. FP16 mixed-precision optimization may need dynamic loss scaling when backward gradients would otherwise underflow, together with finite-value checks on overflow-prone paths. BF16 and inference-only mixed precision do not inherently require loss scaling; safeguards should match the format and measured numerical range. Sensitive reductions and optimizer state still need adequate precision.

(C) Improve conditioning.

• Second-order and quasi-Newton methods (DDP, iLQR, Gauss-Newton, L-BFGS) can reduce sensitivity to eigenvalue spread when the local curvature model is accurate enough and suitably regularized or globalized. Approximate or indefinite curvature can produce an unhelpful step, while regularization choices and finite precision can leave the numerical problem ill-conditioned.
• A compatible, well-chosen preconditioner, for example one based on a positive diagonal approximation to the Gauss-Newton Hessian, can reduce the transformed spectral spread. A poor, stale, or mismatched preconditioner can leave conditioning unchanged or worsen it.
• Architecture choices inside the transition can help, but must be checked on planner-relevant rollout Jacobians. Layer normalization offers no conditioning guarantee. An upper spectral bound can limit expansion without preventing vanishing directions. A residual form $f(z,a)=z+g(z,a)$ stays near identity only when $\partial g/\partial z$ is controlled, not merely when $g$ has small values.
• For a linear component known to be asymptotically stable, Schur or diagonal parameterizations can constrain eigenvalues inside the unit disk. Orthogonal recurrences instead preserve norms with eigenvalues on the unit circle. Neither eigenvalue constraint controls finite-horizon transient growth in a non-normal system; monitor singular values or finite-horizon products separately.

(D) Break the chain.

• Multiple shooting or direct collocation, which lifts state variables and exposes sparse local dynamics constraints instead of fully condensing every dependency into the action objective. The constrained system can still be ill-conditioned.
• A square-root checkpoint schedule, which reduces activation memory from $O(H)$ to $O(\sqrt{H})$ for a uniform chain at the cost of recomputation. Other schedules make different tradeoffs. Checkpointing addresses activation-memory pressure rather than Jacobian conditioning; it is worthwhile when the value of lower memory use outweighs its recomputation cost.
• Explicit Jacobian or finite-horizon sensitivity regularizers over planner-relevant trajectories. Multi-step prediction losses expose rollout error but do not by themselves guarantee well-scaled derivative products.

(E) Avoid differentiating through hard components.

• Straight-through estimators for quantization or discretization: use a hard value in the forward pass and substitute a surrogate backward derivative. In straight-through Gumbel-softmax, that surrogate is the soft relaxation's gradient.
• Score-function estimators when pathwise reparameterization is unavailable. For $\mathbb E_{z\sim p_a}[L(a,z)]$, use the ordinary $\partial_a L$ term together with the likelihood-ratio term when differentiation and integration can be interchanged, the required terms are integrable, and the support is action-independent or its boundary contributions vanish. Nonzero boundary terms need a separate valid estimator. An estimated score must be conditionally unbiased given $z$, or must otherwise make the weighted term jointly unbiased. Score-function estimators often have higher variance than a well-behaved pathwise estimator, so control variates are commonly useful.
• Carefully chosen smooth surrogates: differentiable compliant contact for selected contact transitions, an interior barrier for a differentiable inequality whose initialization and every evaluated iterate remain strictly feasible, or a finite-temperature relaxation of a hard categorical choice. Each changes the original problem and must be validated.

Conventional supervised dynamics models are often fit with one-step prediction error, sometimes augmented with multi-step rollout losses. A standard value-only prediction objective does not explicitly control or guarantee the state-Jacobian spectrum. A model can therefore predict well on its training distribution and still be a poor substrate for gradient planning.

Two remedies target different aspects of the problem. Training with a multi-step prediction objective exposes the model to rollout errors caused by its repeated transitions, but it does not by itself guarantee well-conditioned Jacobians. An explicit regularizer can instead penalize selected Jacobian norms, singular values, or finite-horizon products; a guaranteed bound requires an enforced constraint or suitable parameterization. Suppressing physically real expansion may bias the dynamics. These choices connect to Predictive and Generative Objectives.

The deeper point is that one-step prediction quality and closed-loop decision usefulness are different properties. A model can score well on a one-step benchmark and still produce unusable plans. Part XI will revisit this distinction in evaluation. The architecture and learned parameters determine the realized derivatives; data, objective, regularization, initialization, and optimization influence which parameters are learned. Those derivatives are one major factor in how well a gradient planner works, alongside model accuracy, cost design, constraints, horizon, and optimizer.

The numerical values in the worked example are not task-independent defaults. A measurement-driven starting procedure is:

1. Nondimensionalize or rescale each state and action coordinate using documented characteristic scales so their numerical magnitudes are comparable, and transform costs and bounds consistently. Then measure representative objective and gradient scales; coordinate normalization alone does not make gradients order one.
2. Choose a physical lookahead from the system's controllability and task timescales, sampling interval, terminal-value quality, and compute budget. Expressing the horizon only as a number of steps hides these dependencies.
3. Tune the optimizer and learning rate from measured objective decrease. If you use clipping, set its threshold from the observed unclipped gradient distribution and log how often it activates.
4. When consecutive receding-horizon problems evolve smoothly, shift the unused suffix, append a terminal action, and check or restore feasibility for the newly measured state and constraints before using the full-length sequence as a warm start. Reset or use additional starts after abrupt disturbances, reference changes, constraint changes, or mode switches.
5. Log gradients by time step. A large span is a diagnostic trigger, not a verdict: inspect contraction or expansion, action gains, cost placement, cancellations, constraints, units, and horizon. For curvature, use matrix-free Hessian-vector products to estimate the most negative and largest algebraic eigenvalues of a symmetric action Hessian, or Jacobian-vector and vector-Jacobian products to estimate extreme singular values of a nonsymmetric rollout Jacobian or sensitivity product. Form the matrix directly only when the action problem is small enough.
6. If progress stalls, first check gradient correctness, scale, saturation, active constraints, smoothness, and initialization. Consider regularized iLQR or DDP when curvature is the limiting factor and a useful local quadratic model exists.

We use a problem whose action gradient and optimum have closed matrix forms that we can verify numerically. The continuous dynamics are $\dot{x}=v$ and $\dot{v}=a$, with acceleration held constant over each sampling interval. Their exact zero-order-hold discretization at time step $\Delta t$ is

$$z_t = \begin{bmatrix} x_t \\ v_t \end{bmatrix}, \qquad z_{t+1} = \underbrace{\begin{bmatrix} 1 & \Delta t \\ 0 & 1 \end{bmatrix}}_{A} z_t + \underbrace{\begin{bmatrix} \tfrac{1}{2}\Delta t^2 \\ \Delta t \end{bmatrix}}_{B} a_t.$$

The cost is quadratic with target $g$:

$$J(a) = \sum_{t=0}^{H} (z_t - g)^\top Q (z_t - g) + \sum_{t=0}^{H-1} R a_t^2.$$

This linear-quadratic structure makes the action gradient affine and the unconstrained optimum available from a linear solve. It is therefore a useful test case: we can compare autodiff with an independently assembled quadratic gradient at a declared action vector.

Because the dynamics are linear, the state at time $k$ is an affine function of the whole action sequence:

$$z_k = A^k z_0 + \sum_{j=0}^{k-1} A^{k-1-j} B \, a_j = A^k z_0 + G_k a,$$

where $a \in \mathbb{R}^H$ collects the actions and $G_k \in \mathbb{R}^{2 \times H}$ has columns

$$[G_k]_{:,j} = \begin{cases} A^{k-1-j} B, & j \le k-1, \\ 0, & j \ge k. \end{cases}$$

Substituting into the cost and collecting terms gives a quadratic function of $a$:

$$J(a) = a^\top M a + 2 a^\top b + \text{const}, \qquad M = \sum_{k=0}^{H} G_k^\top Q_k G_k + R I_H, \qquad b = \sum_{k=0}^{H} G_k^\top Q_k \big(A^k z_0 - g\big),$$

with $Q_k=Q$ for all $k$ in this example. The gradient is affine in the actions:

$$\nabla_a J = 2 M a + 2 b.$$

Three consequences follow. The conditioning and nonlinear-gradient lessons extend beyond this toy problem, while strict convexity and the unique global optimum rely on this linear-quadratic structure with $R>0$.

The objective is strictly convex. $M$ is a sum of positive semidefinite terms plus $RI_H$, so $M\succ0$ when $R>0$. The unique local minimum is the global minimum,

$$a^\star=-M^{-1}b.$$

With exact gradients and a suitable step size, gradient descent converges to it; an excessive fixed step can still diverge.

For a suitably scaled fixed step, the convergence bound is governed by the condition number of $M$. Near the optimum, equal displacements along eigenvectors with small eigenvalues raise the cost less than displacements along steep eigenvectors. A single stable scalar step then contracts slowly along the flat directions. For this double integrator, the eigenvalue spread grows over the horizons we test; that trend is not universal.

For nonlinear multi-step dynamics, the gradient is generally not affine in the actions. Pendulum, cart-pole, and neural-dynamics planning objectives are often non-convex, although special nonlinear cases can remain convex. If the objective is bounded below and has an $L$-Lipschitz gradient, fixed-step gradient descent with $0<\alpha<2/L$ drives the gradient norm toward zero. Any convergent iterate sequence therefore converges to a stationary point, which need not be a local minimum. The result can depend strongly on initialization, and warm starts can help retain a useful basin when consecutive problems change smoothly.

We now build the planner. Each experiment isolates one point: autograd matches the analytic gradient at a tested action vector, Adam improves an open-loop plan, a Riccati terminal value improves a two-step MPC controller under a fixed optimization budget, and scalar dynamics expose contracting and expanding gradient chains.

The chapter-wide book setup runs before the first executable example so its plotting helpers are available throughout. Seeding controls the common random streams and improves repeatability within a fixed software and hardware environment; it does not by itself guarantee identical results across platforms or library versions.

Next we define the double integrator, the quadratic cost, and the horizon. Keeping the constants in one place makes it easy to vary the horizon later.

The two gradients agree to numerical precision at this tested point. We check the gradient norm explicitly because agreement between two negligible vectors would be a weak test. Here the all-zero diagnostic vector has a substantial nonzero gradient.

Now we hand the action sequence to Adam. The "parameters" are actions rather than model weights; the loop still zeroes gradients, computes the loss, backpropagates, and takes an update.

The objective decreases throughout the plotted run. Under this fixed budget, the trajectory approaches and overshoots the target position, ending near $(x,v)=(1.05,-0.28)$ rather than settling exactly at $(1,0)$. In this run, the optimizer produces rising and then falling velocity while trading finite-horizon position and velocity tracking against action effort. No separately designed feedback law is used: given the model, cost, initialization, and optimizer, backpropagation through the known dynamics computes the action gradients, and Adam uses those gradients to produce the open-loop action sequence.

We now assemble $M$ for increasing horizons in this fixed double-integrator problem. The true action Hessian is $2M$, which has the same condition number as $M$. Horizon growth need not worsen every system, but it does here.

We solve the discrete algebraic Riccati equation numerically, then compare three receding-horizon controllers under the same 100-iteration Adam budget per planning call: a declared two-step horizon without a terminal value, the same two-step horizon with the Riccati value, and a twelve-step reference without a terminal value. Finite optimizer budgets mean this is an empirical comparison of these implementations, not an exact Bellman-equivalence test.