可微分物理仿真是近年来物理仿真领域最重要的突破之一，它将传统的前向动力学仿真转变为可以计算梯度的计算图，从而支持基于梯度的优化方法。这一技术革新不仅为轨迹优化、系统辨识和控制器设计提供了新的工具，更重要的是，它将物理仿真无缝集成到深度学习框架中，实现了端到端的学习。本章将深入探讨可微分仿真的数学基础、实现技术以及在强化学习中的应用，特别关注如何处理接触等非光滑现象带来的挑战。
完成本章学习后，您将能够：
自动微分（Automatic Differentiation, AD）是计算函数导数的一种算法技术，它既不同于符号微分的表达式操作，也不同于数值微分的有限差分近似。AD的核心思想是将复杂函数分解为基本运算的组合，然后通过链式法则精确计算导数。
考虑一个简单的例子，计算函数 $f(x_1, x_2) = \sin(x_1 x_2) + x_1^2$ 的梯度：
计算图表示：
x1 ──┬──> (*) ──> sin ──> (+) ──> f
│ ↑ ↑
└──> x2 x1^2
前向模式AD从输入开始，逐层传播导数：
反向模式AD（即反向传播）从输出开始，逐层传播梯度：
在物理仿真中，我们需要对整个仿真轨迹计算梯度。设仿真的离散时间演化为：
\[\mathbf{x}_{t+1} = f(\mathbf{x}_t, \mathbf{u}_t, \theta)\]其中 $\mathbf{x}t$ 是状态，$\mathbf{u}_t$ 是控制输入，$\theta$ 是系统参数。对于损失函数 $L = \sum{t=0}^T l(\mathbf{x}_t, \mathbf{u}_t)$，我们需要计算：
\[\frac{\partial L}{\partial \theta} = \sum_{t=0}^T \frac{\partial l}{\partial \mathbf{x}_t} \frac{\partial \mathbf{x}_t}{\partial \theta}\]关键在于递归计算 $\partial \mathbf{x}_t/\partial \theta$：
\[\frac{\partial \mathbf{x}_{t+1}}{\partial \theta} = \frac{\partial f}{\partial \mathbf{x}_t} \frac{\partial \mathbf{x}_t}{\partial \theta} + \frac{\partial f}{\partial \theta}\]这是一个线性递推关系，可以通过前向或反向模式高效计算。
虽然自动微分提供了通用的梯度计算方法，但在某些情况下，手工推导的解析梯度具有优势：
例如，对于刚体旋转的四元数表示 $\mathbf{q} = [\cos(\theta/2), \sin(\theta/2)\mathbf{n}]$，其对旋转角 $\theta$ 的导数为：
\[\frac{\partial \mathbf{q}}{\partial \theta} = \frac{1}{2}[-\sin(\theta/2), \cos(\theta/2)\mathbf{n}]\]直接使用这个解析表达式比通过AD计算更高效且稳定。
实践中，我们通常采用混合策略：
物理仿真中的许多问题可以表示为约束优化：
\[\min_{\mathbf{x}} f(\mathbf{x}, \theta) \quad \text{s.t.} \quad g(\mathbf{x}, \theta) = 0\]其KKT条件为： \(\begin{cases} \nabla_x f + \lambda^T \nabla_x g = 0 \\ g(\mathbf{x}, \theta) = 0 \end{cases}\)
根据隐式函数定理，如果雅可比矩阵： \(\mathbf{J} = \begin{bmatrix} \nabla_{xx}^2 f + \lambda^T \nabla_{xx}^2 g & \nabla_x g^T \\ \nabla_x g & 0 \end{bmatrix}\)
是非奇异的，则解 $(\mathbf{x}^, \lambda^)$ 关于参数 $\theta$ 是可微的。
考虑隐式欧拉方法： \(\mathbf{x}_{t+1} = \mathbf{x}_t + h \mathbf{f}(\mathbf{x}_{t+1}, t+h)\)
这定义了一个隐式关系 $\mathbf{g}(\mathbf{x}{t+1}, \mathbf{x}_t) = \mathbf{x}{t+1} - \mathbf{x}t - h\mathbf{f}(\mathbf{x}{t+1}) = 0$。
应用隐式函数定理： \(\frac{\partial \mathbf{x}_{t+1}}{\partial \mathbf{x}_t} = -\left(\frac{\partial \mathbf{g}}{\partial \mathbf{x}_{t+1}}\right)^{-1} \frac{\partial \mathbf{g}}{\partial \mathbf{x}_t} = (\mathbf{I} - h\mathbf{J}_f)^{-1}\)
其中 $\mathbf{J}f = \partial \mathbf{f}/\partial \mathbf{x}{t+1}$。
线性互补问题（LCP）：找到 $\mathbf{z} \geq 0, \mathbf{w} \geq 0$ 使得： \(\mathbf{w} = \mathbf{M}\mathbf{z} + \mathbf{q}, \quad \mathbf{w}^T\mathbf{z} = 0\)
可以重新表示为光滑的非线性方程组： \(\mathbf{\Phi}(\mathbf{z}, \mathbf{w}) = \begin{bmatrix} \mathbf{w} - \mathbf{M}\mathbf{z} - \mathbf{q} \\ \phi(\mathbf{w}, \mathbf{z}) \end{bmatrix} = 0\)
其中 $\phi$ 是光滑互补函数，如Fischer-Burmeister函数： \(\phi(a, b) = a + b - \sqrt{a^2 + b^2 + \epsilon^2}\)
参数 $\epsilon > 0$ 提供了正则化，使函数处处可微。
对于接触约束 $\mathbf{c}(\mathbf{q}) \geq 0$，接触力通过以下优化问题确定：
\[\min_{\lambda \geq 0} \frac{1}{2}\lambda^T \mathbf{A} \lambda + \lambda^T \mathbf{b} \quad \text{s.t.} \quad \mathbf{c}(\mathbf{q}) + \mathbf{J}_c \Delta \mathbf{v} \geq 0\]其中 $\mathbf{J}_c = \partial \mathbf{c}/\partial \mathbf{q}$ 是接触雅可比。解的导数可以通过隐式函数定理计算：
\[\frac{\partial \lambda^*}{\partial \mathbf{q}} = -(\mathbf{A}_{\mathcal{A}\mathcal{A}})^{-1} \mathbf{J}_{c,\mathcal{A}}^T \frac{\partial \mathbf{c}}{\partial \mathbf{q}}\]其中 $\mathcal{A}$ 是活跃约束集。
传统的接触检测是二值函数： \(\text{contact}(\mathbf{x}) = \begin{cases} 1 & \text{if } d(\mathbf{x}) \leq 0 \\ 0 & \text{otherwise} \end{cases}\)
为了可微分性，我们使用光滑近似： \(\text{contact}_\epsilon(\mathbf{x}) = \sigma(-d(\mathbf{x})/\epsilon)\)
其中 $\sigma$ 是sigmoid函数，$\epsilon$ 控制过渡区域的宽度。
硬接触约束 $\mathbf{c} \geq 0, \lambda \geq 0, \lambda^T\mathbf{c} = 0$ 可以用软约束近似：
\[\lambda_n = k_n \max(0, -d)^p + d_n \max(0, -\dot{d})\]其中 $k_n$ 是刚度，$d_n$ 是阻尼，$p$ 是指数（通常取1.5模拟赫兹接触）。
为了保证处处可微，使用光滑max函数： \(\text{smax}(x, \epsilon) = \epsilon \log(1 + e^{x/\epsilon})\)
| 库仑摩擦模型 $ | \mathbf{f}_t | \leq \mu f_n$ 是不可微的。常用的光滑近似包括： |
粘性摩擦模型： \(\mathbf{f}_t = -\mu f_n \tanh(\mathbf{v}_t / v_0)\)
光滑摩擦锥： \(\mathbf{f}_t = -\mu f_n \frac{\mathbf{v}_t}{|\mathbf{v}_t| + \epsilon}\)
LuGre模型（考虑粘滑转换）： \(\begin{aligned} \dot{z} &= v - \sigma_0 \frac{|v|}{g(v)} z \\ f &= \sigma_0 z + \sigma_1 \dot{z} + \sigma_2 v \end{aligned}\)
| 其中 $g(v) = \mu_c + (\mu_s - \mu_c)e^{- | v/v_s | ^2}$。 |
对于点接触，接触力对位置的雅可比包含几何和力学两部分：
\[\frac{\partial \mathbf{f}_c}{\partial \mathbf{q}} = \frac{\partial \mathbf{f}_c}{\partial d} \frac{\partial d}{\partial \mathbf{q}} + \frac{\partial \mathbf{f}_c}{\partial \mathbf{n}} \frac{\partial \mathbf{n}}{\partial \mathbf{q}}\]其中：
对于凸形状，可以使用GJK算法的中间结果计算这些导数。
考虑T步仿真的梯度： \(\frac{\partial L}{\partial \mathbf{x}_0} = \sum_{t=0}^T \frac{\partial l_t}{\partial \mathbf{x}_t} \prod_{k=0}^{t-1} \frac{\partial \mathbf{x}_{k+1}}{\partial \mathbf{x}_k}\)
| 如果 $ | \partial \mathbf{x}_{k+1}/\partial \mathbf{x}_k | > 1$，梯度指数增长（爆炸）； |
| 如果 $ | \partial \mathbf{x}_{k+1}/\partial \mathbf{x}_k | < 1$，梯度指数衰减（消失）。 |
梯度裁剪： \(\mathbf{g}_{\text{clip}} = \begin{cases} \mathbf{g} & \text{if } |\mathbf{g}| \leq \theta \\ \theta \frac{\mathbf{g}}{|\mathbf{g}|} & \text{otherwise} \end{cases}\)
层归一化（应用于状态转换）： \(\mathbf{x}_{t+1} = \gamma \frac{f(\mathbf{x}_t) - \mu}{\sigma + \epsilon} + \beta\)
其中 $\mu, \sigma$ 是均值和标准差，$\gamma, \beta$ 是可学习参数。
对于长序列，可以使用截断的时间反向传播（TBPTT）：
实现伪代码：
for epoch in epochs:
states = [x0]
for t in range(T):
states.append(forward(states[-1]))
if (t+1) % K == 0:
loss = compute_loss(states[-K:])
backward(loss)
states = [states[-1].detach()] # 截断梯度
使用可学习或自适应的时间步长可以缓解梯度问题：
\[h_t = h_0 \cdot \sigma(\mathbf{w}^T \phi(\mathbf{x}_t))\]其中 $\phi$ 是特征函数，$\mathbf{w}$ 是可学习参数。这允许系统在稳定区域使用大步长，在复杂区域使用小步长。
可微分仿真使得我们可以直接优化控制序列：
\[\min_{\mathbf{u}_{0:T}} \sum_{t=0}^T l(\mathbf{x}_t, \mathbf{u}_t) \quad \text{s.t.} \quad \mathbf{x}_{t+1} = f(\mathbf{x}_t, \mathbf{u}_t)\]使用梯度下降： \(\mathbf{u}_{0:T}^{(k+1)} = \mathbf{u}_{0:T}^{(k)} - \alpha \nabla_{\mathbf{u}} L\)
其中梯度通过可微分仿真器反向传播计算。
将可微分仿真集成到策略梯度中：
\[\nabla_\theta J = \mathbb{E}_{\pi_\theta}\left[\sum_{t=0}^T \nabla_\theta \log \pi_\theta(\mathbf{u}_t|\mathbf{x}_t) R_t + \lambda \nabla_\theta L_{\text{sim}}\right]\]其中 $L_{\text{sim}}$ 是通过可微分仿真计算的辅助损失，如预测误差或约束违背。
使用真实数据优化仿真参数：
\[\theta^* = \arg\min_\theta \sum_{i=1}^N |\mathbf{x}_i^{\text{real}} - \mathbf{x}_i^{\text{sim}}(\theta)|^2\]梯度 $\partial \mathbf{x}^{\text{sim}}/\partial \theta$ 通过可微分仿真计算，实现快速参数调优。
Robert Edwin Wengert在1964年的博士论文中首次系统地提出了自动微分的概念。他认识到，任何计算机程序都可以视为基本运算的序列，而每个基本运算的导数都是已知的。通过系统地应用链式法则，可以自动计算整个程序的导数。
Wengert的关键洞察是引入”带号表”（tape）的概念，记录计算过程中的所有中间变量及其依赖关系。这个想法直接导致了现代深度学习框架中计算图的概念。虽然Wengert的工作在当时没有立即得到广泛应用（部分因为计算资源限制），但它为后来的自动微分工具（如ADIFOR、TAF）和现代深度学习框架（如TensorFlow、PyTorch）奠定了理论基础。
有趣的是，Wengert最初的动机是解决工程优化问题，特别是化工过程的参数优化。他没有预见到60年后，他的方法会成为人工智能革命的基石。这再次证明了基础研究的长远价值——看似抽象的数学工具可能在未来产生革命性的应用。
神经ODE将离散的层替换为连续的动力学：
\[\frac{d\mathbf{x}}{dt} = f_\theta(\mathbf{x}(t), t)\]其中 $f_\theta$ 是神经网络。这与物理仿真有天然的联系——两者都是求解ODE。
计算梯度时，神经ODE使用伴随方法避免存储整个轨迹：
\[\frac{d\mathbf{a}}{dt} = -\mathbf{a}^T \frac{\partial f}{\partial \mathbf{x}}, \quad \mathbf{a}(T) = \frac{\partial L}{\partial \mathbf{x}(T)}\]然后 $\partial L/\partial \mathbf{x}(0) = \mathbf{a}(0)$。这与可微分物理仿真的反向模式AD等价。
PINN将物理定律作为损失函数的一部分：
\[L = L_{\text{data}} + \lambda L_{\text{physics}}\]其中： \(L_{\text{physics}} = \int_\Omega \left|\frac{\partial u}{\partial t} - \mathcal{L}[u]\right|^2 dx\)
$\mathcal{L}$ 是物理算子（如Navier-Stokes方程）。
保持能量守恒的神经网络架构：
\[\begin{aligned} \frac{d\mathbf{q}}{dt} &= \frac{\partial H_\theta}{\partial \mathbf{p}} \\ \frac{d\mathbf{p}}{dt} &= -\frac{\partial H_\theta}{\partial \mathbf{q}} \end{aligned}\]其中 $H_\theta$ 是参数化的哈密顿量。这保证了辛结构和长期稳定性。
可微分物理仿真通过将传统的前向动力学仿真转化为可计算梯度的系统，为轨迹优化、系统辨识和端到端学习开辟了新的可能。本章的关键要点包括：
练习10.1 手工计算函数 $f(x,y) = e^{xy} + \sin(x+y)$ 在点 $(1,0)$ 处的梯度，分别使用前向和反向模式自动微分。验证两种方法得到相同结果。
练习10.2 考虑弹簧-质量系统 $m\ddot{x} + c\dot{x} + kx = 0$，使用隐式欧拉方法离散化。推导位置 $x_{n+1}$ 对参数 $k$ 的导数。
| 练习10.3 实现一维接触的光滑化模型。给定距离函数 $d(x)$，设计一个处处可微的接触力函数，满足：当 $d > \epsilon$ 时力约为0，当 $d < -\epsilon$ 时力约为 $k | d | $。 |
练习10.4 考虑二维刚体碰撞，设计一个可微分的碰撞响应模型。刚体A和B的状态为 $(x_A, y_A, \theta_A, v_{xA}, v_{yA}, \omega_A)$ 和类似的B状态。当检测到碰撞时，如何计算可微分的冲量？
练习10.5 分析为什么在可微分仿真中，显式积分器比隐式积分器更容易出现梯度爆炸。设计一个自适应的混合积分策略。
练习10.6 推导并实现可微分的LCP求解器。给定 $\mathbf{M}\mathbf{z} + \mathbf{q} \geq 0$, $\mathbf{z} \geq 0$, $\mathbf{z}^T(\mathbf{M}\mathbf{z} + \mathbf{q}) = 0$，如何计算解 $\mathbf{z}^*$ 对参数 $\mathbf{q}$ 的导数？
练习10.7（开放性问题）设计一个可微分的绳索仿真系统。绳索由N个质点通过约束连接，需要处理拉伸、弯曲和自碰撞。讨论主要的技术挑战和可能的解决方案。
练习10.8 证明使用辛积分器的可微分仿真保持梯度的辛结构。这对长期轨迹优化有什么影响？
问题：为了可微分性使用过大的 $\epsilon$，导致物理行为失真。 解决：使用退火策略，训练初期用大 $\epsilon$，逐步减小。
问题：即使使用光滑近似，在某些点（如 $v = 0$ 的摩擦力）梯度仍可能不连续。 解决：在这些点附近使用次梯度或随机平滑。
问题：存储整个仿真轨迹用于反向传播消耗大量内存。 解决：使用梯度检查点或递归反向传播。
问题：长链的梯度传播累积浮点误差。 解决：使用混合精度训练，关键部分用双精度。
问题：可微分仿真的优化景观复杂，容易陷入局部最优。 解决：结合全局优化方法（如CMA-ES）或使用多起点策略。