Structured Simulators and Neural Operators › Differentiable Simulators

A simulator is a world model you wrote yourself. It takes a state $x_t$ and an action $a_t$ and produces a next state through hand-coded physics:

$$x_{t+1} = f(x_t, a_t; \theta),$$

where:

• $x_t$, $a_t$: the state and the action at time $t$
• $x_{t+1}$: the next state produced by the hand-coded physics
• $f$: the hand-coded transition function that advances the state one step
• $\theta$: the physical parameters, collecting masses, lengths, friction coefficients, stiffnesses

The distinction between the structure of $f$ and the values of $\theta$ is the whole design. The structure encodes your physical understanding: which forces act, how they compose, what integration scheme advances the state. The parameters encode what you do not know precisely: the true mass of the payload, the actual coefficient of friction on the table, the real stiffness of a contact. Many legacy or externally wrapped simulators expose only forward outputs, even though classical sensitivity machinery may exist. A differentiable simulator additionally exposes derivatives of the implemented computation through analytic, automatic, implicit, adjoint, or custom derivative rules. For the rollout above, the quantities of interest are the following partial derivatives:

$$\frac{\partial x_T}{\partial \theta}, \qquad \frac{\partial x_T}{\partial x_0}, \qquad \frac{\partial x_T}{\partial a_{0:T-1}}.$$

where:

• $T$: the rollout horizon, the number of simulation steps taken from the initial state
• $x_T$: the simulator state after $T$ steps, the output whose sensitivity we compute
• $\partial x_T / \partial \theta$: how the final state changes with each physical parameter
• $\partial x_T / \partial x_0$: how the final state changes with the initial condition
• $\partial x_T / \partial a_{0:T-1}$: how the final state changes with the entire sequence of actions applied over the rollout

Each of these three gradients powers a different capability. The parameter gradient drives system identification, the initial-condition gradient drives sensitivity analysis and state estimation, and the action gradient drives trajectory optimization. One derivative mechanism enables all three, which makes differentiability valuable for a world model.

Why does this matter? Three reasons dominate.

System identification becomes gradient descent. Part II, Chapter 3 treated system identification as an optimization problem and surveyed gradient-free and linear methods. Those methods are workhorses, but they can scale poorly with dimension and nonlinear dynamics. When the simulator is differentiable, identification of $\theta$ reduces to minimizing a loss $L(x_{1:T}^{\text{data}}; \theta)$ by backpropagation. The method avoids choosing parameter updates from forward evaluations alone; each iteration still simulates the trajectory, compares it with data, backpropagates through the rollout, and updates the parameters. The gradient combines information from each step of the rollout.

Planning becomes gradient descent. If you can differentiate the accumulated reward or a target cost through the rollout, you can optimize an action sequence directly. This is the modern descendant of shooting methods from classical optimal control (Part II, Chapter 4), and we will do exactly this in the worked example. Automatic differentiation replaces much of the system-specific derivative transcription with loss.backward(), provided every operation on the executed path has a valid derivative. The resulting gradient is exact for the implemented computation, not necessarily for omitted or idealized physics.

Learning can wrap around the simulator. Gradients flowing through simulation let you fit residual neural networks on top of physics, fit simulator parameters jointly with policy parameters, or backpropagate a learned perception module into a physical model. This pattern appears in some hybrid differentiable-simulation pipelines: the simulator supplies a structured backbone, while learned components model residual effects and end-to-end derivatives coordinate the trainable pieces. It is one design option for sim-to-real work, not a description of every sim-to-real pipeline.

The rollout is a deep computation graph: applying $f$ once is one layer, so a horizon-$T$ rollout is a $T$-layer network. Everything you know about training deep networks applies, with the twist that the layers are physics steps and the depth is dictated by the horizon you care about, not by an architecture choice. There are two standard ways to get gradients through it.

Unrolled backpropagation treats the whole trajectory as a feedforward graph and applies reverse-mode autodiff. This is what torch.autograd does automatically. It is exact and simple, and because reverse-mode autodiff computes the gradient of a scalar loss with respect to all inputs in a single backward pass whose cost is comparable to the forward pass, it is remarkably efficient per unit of information. Its weakness is memory: without checkpointing or recomputation, it stores intermediate states and memory grows linearly with $T$. For a long rollout of a high-dimensional state, memory can become the binding constraint, depending on state size, per-step compute, hardware, and the checkpointing or recomputation strategy.

The adjoint method avoids storing the trajectory by solving an augmented system backward in time. For continuous dynamics $\dot{x} = f(x, t; \theta)$ and a terminal cost $J = \Phi(x_T)$, which measures the terminal state, define the adjoint state $\lambda(t)$ by

$$\dot{\lambda} = -\left(\frac{\partial f}{\partial x}\right)^{\!\top} \lambda, \qquad \lambda(T) = \frac{\partial \Phi}{\partial x_T},$$

where:

• $\lambda(t)$: the adjoint state, a costate variable that accumulates sensitivity information as it is integrated backward in time from $t = T$ to $t = 0$
• $\partial f / \partial x$: the Jacobian of the vector field with respect to the state, evaluated along the forward trajectory
• $\lambda(T)$: the terminal condition, equal to the gradient of the terminal cost $\Phi$ with respect to the final state $x_T$

and then the gradient is an integral along the backward pass:

$$\frac{dJ}{d\theta} = -\int_T^0 \lambda^\top \frac{\partial f}{\partial \theta}\, dt.$$

where:

• $dJ/d\theta$: the gradient of the terminal cost with respect to the parameters
• $\lambda^\top$: the adjoint state (transposed), propagated backward from $T$ to $0$ along the reverse integration
• $\partial f / \partial \theta$: the Jacobian of the vector field with respect to the parameters
• The integral runs from $T$ down to $0$ because the adjoint system is solved backward in time; the minus sign comes from that reversed direction of integration

The intuition behind the adjoint state is worth pausing on: $\lambda(t)$ answers the question "if you perturb the state slightly at time $t$, how does the final cost change?" It starts at $T$ holding exactly the cost's sensitivity to the terminal state, and as it integrates backward, each step multiplies it by the transpose Jacobian of the dynamics, which converts sensitivity one step later into sensitivity one step earlier. The parameter gradient then follows by accumulating, along the backward sweep, how each parameter perturbation at each moment feeds into that running sensitivity. In the ideal continuous-adjoint formulation, storage can be independent of the number of accepted steps when the trajectory can be reconstructed. Practical implementations retain solver workspace or checkpoints and trade storage for recomputation. There is a pleasing unity here: discrete unrolled backpropagation is in fact an adjoint method for the discretized system; the distinction is practical, not fundamental. Which side wins depends on state dimension, horizon, numerical stability, and whether the integrator can reconstruct or restart the trajectory reliably.

Physics is not everywhere smooth, and this is the central engineering difficulty. The clean calculus of the previous section assumes $f$ is differentiable, and the most interesting parts of physical interaction are precisely where it is not.

• Contacts and impacts create hybrid maps whose event activation, timing, contact set, and normal can change non-smoothly with state or parameters. Within a fixed impact mode, a restitution map such as $v_n^+ = -e v_n^-$ is linear and continuous in the incoming normal velocity; the derivative difficulty comes from switching between modes and moving event times.
• Friction and stiction involve complementarity conditions (either the object sticks with force below a threshold, or it slips with force at the threshold), which are non-smooth by nature. The switch between sticking and slipping is a set-valued force law, and set-valued maps have no single Jacobian at the switching surface.
• Chaotic dynamics amplify perturbations along unstable directions approximately as $e^{\lambda T}$ over the regime where linearization applies. After several Lyapunov times, long-horizon prediction and gradient-based optimization can become badly conditioned. The derivative can still be mathematically correct; the problem is that it describes extreme sensitivity to initial conditions, parameters, and numerical error.

Practical responses depend on the engine and contact pipeline. Some differentiable simulators use compliant or softened contact, while others differentiate constraint, convex, or position-based formulations; frameworks such as DiffTaichi and Warp support more than one simulator formulation. Any surrogate changes the relationship between the implemented gradient and ideal rigid-contact physics. Against chaotic sensitivity, practitioners often differentiate over shorter windows and re-fit or re-plan frequently, exactly the receding-horizon pattern we will see in the worked example. An adjoint changes memory and computation, not the underlying sensitivity.

We can isolate one mechanism by computing parameter sensitivity in two linear systems. One transition contracts perturbations and the other amplifies them. This is not a chaos experiment; it shows how repeated Jacobian products can make an additive parameter sensitivity saturate or grow with rollout length.