第10章：可微分物理仿真

开篇导言

可微分物理仿真是近年来物理仿真领域最重要的突破之一，它将传统的前向动力学仿真转变为可以计算梯度的计算图，从而支持基于梯度的优化方法。这一技术革新不仅为轨迹优化、系统辨识和控制器设计提供了新的工具，更重要的是，它将物理仿真无缝集成到深度学习框架中，实现了端到端的学习。本章将深入探讨可微分仿真的数学基础、实现技术以及在强化学习中的应用，特别关注如何处理接触等非光滑现象带来的挑战。

学习目标

完成本章学习后，您将能够：

理解自动微分与符号微分的区别，掌握前向和反向自动微分的实现原理
运用隐式函数定理计算约束系统的梯度，理解KKT条件的可微分性
实现接触力的光滑近似方法，处理摩擦锥的可微分松弛
诊断和解决梯度爆炸与消失问题，设计稳定的可微分仿真器
将可微分仿真集成到强化学习训练流程中，实现基于梯度的轨迹优化

10.1 自动微分与解析梯度

10.1.1 自动微分的基本原理

自动微分（Automatic Differentiation, AD）是计算函数导数的一种算法技术，它既不同于符号微分的表达式操作，也不同于数值微分的有限差分近似。AD的核心思想是将复杂函数分解为基本运算的组合，然后通过链式法则精确计算导数。

考虑一个简单的例子，计算函数 $f(x_1, x_2) = \sin(x_1 x_2) + x_1^2$ 的梯度：

计算图表示：
    x1 ──┬──> (*) ──> sin ──> (+) ──> f
         │     ↑                ↑
         └──> x2              x1^2

前向模式AD从输入开始，逐层传播导数：

设 $\dot{x}_1 = 1, \dot{x}_2 = 0$（计算 $\partial f/\partial x_1$）
$v_1 = x_1 x_2$，$\dot{v}_1 = \dot{x}_1 x_2 + x_1 \dot{x}_2 = x_2$
$v_2 = \sin(v_1)$，$\dot{v}_2 = \cos(v_1) \dot{v}_1 = x_2 \cos(x_1 x_2)$
$v_3 = x_1^2$，$\dot{v}_3 = 2x_1 \dot{x}_1 = 2x_1$
$f = v_2 + v_3$，$\dot{f} = \dot{v}_2 + \dot{v}_3 = x_2 \cos(x_1 x_2) + 2x_1$

反向模式AD（即反向传播）从输出开始，逐层传播梯度：

设 $\bar{f} = 1$
$\bar{v}_2 = \bar{f} \cdot 1 = 1$，$\bar{v}_3 = \bar{f} \cdot 1 = 1$
$\bar{v}_1 = \bar{v}_2 \cdot \cos(v_1) = \cos(x_1 x_2)$
$\bar{x}_1 = \bar{v}_1 \cdot x_2 + \bar{v}_3 \cdot 2x_1 = x_2 \cos(x_1 x_2) + 2x_1$
$\bar{x}_2 = \bar{v}_1 \cdot x_1 = x_1 \cos(x_1 x_2)$

10.1.2 物理仿真中的自动微分实现

在物理仿真中，我们需要对整个仿真轨迹计算梯度。设仿真的离散时间演化为：

$$\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}$$ 这是一个线性递推关系，可以通过前向或反向模式高效计算。

10.1.3 解析梯度 vs 自动微分

虽然自动微分提供了通用的梯度计算方法，但在某些情况下，手工推导的解析梯度具有优势：

计算效率：对于特定结构（如稀疏雅可比），解析梯度可以利用结构信息
数值稳定性：某些操作（如归一化）的解析梯度可以避免数值问题
内存效率：不需要存储中间计算图

例如，对于刚体旋转的四元数表示 $\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计算更高效且稳定。

10.1.4 混合微分策略

实践中，我们通常采用混合策略：

对复杂的数值积分器使用AD
对简单的几何变换使用解析梯度
对不可微操作（如碰撞检测）使用光滑近似

10.2 隐式函数定理的应用

10.2.1 约束优化问题的可微分性

物理仿真中的许多问题可以表示为约束优化： $$\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$ 是可微的。

10.2.2 隐式积分器的梯度计算

考虑隐式欧拉方法： $$\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}$。

10.2.3 LCP求解器的可微分性

线性互补问题（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$ 提供了正则化，使函数处处可微。

10.2.4 接触力的隐式微分

对于接触约束 $\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}$ 是活跃约束集。

10.3 接触的可微分处理

10.3.1 接触检测的光滑化

传统的接触检测是二值函数： $$\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$ 控制过渡区域的宽度。

10.3.2 法向接触力的软约束模型

硬接触约束 $\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})$$

10.3.3 摩擦力的可微分模型

库仑摩擦模型 $|\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}$。

10.3.4 接触雅可比的计算

对于点接触，接触力对位置的雅可比包含几何和力学两部分： $$\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}}$$ 其中：

$\partial d/\partial \mathbf{q}$ 是距离函数的梯度（通常是法向量）
$\partial \mathbf{n}/\partial \mathbf{q}$ 涉及曲率信息

对于凸形状，可以使用GJK算法的中间结果计算这些导数。

10.4 梯度爆炸与消失问题

10.4.1 长序列仿真中的梯度传播

考虑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$，梯度指数衰减（消失）。

10.4.2 梯度裁剪与归一化

梯度裁剪： $$\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$ 是可学习参数。

10.4.3 时间反向传播的截断

对于长序列，可以使用截断的时间反向传播（TBPTT）：

每K步截断梯度流
保持前向状态连续性
在截断点重新初始化梯度

实现伪代码：

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()]  # 截断梯度

10.4.4 自适应时间步长

使用可学习或自适应的时间步长可以缓解梯度问题： $$h_t = h_0 \cdot \sigma(\mathbf{w}^T \phi(\mathbf{x}_t))$$ 其中 $\phi$ 是特征函数，$\mathbf{w}$ 是可学习参数。这允许系统在稳定区域使用大步长，在复杂区域使用小步长。

RL应用：基于梯度的轨迹优化

应用场景1：模型预测控制

可微分仿真使得我们可以直接优化控制序列： $$\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$$ 其中梯度通过可微分仿真器反向传播计算。

应用场景2：策略梯度增强

将可微分仿真集成到策略梯度中： $$\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}}$ 是通过可微分仿真计算的辅助损失，如预测误差或约束违背。

应用场景3：系统辨识与域适应

使用真实数据优化仿真参数： $$\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$ 通过可微分仿真计算，实现快速参数调优。

历史人物：Wengert与自动微分的诞生

Robert Edwin Wengert在1964年的博士论文中首次系统地提出了自动微分的概念。他认识到，任何计算机程序都可以视为基本运算的序列，而每个基本运算的导数都是已知的。通过系统地应用链式法则，可以自动计算整个程序的导数。

Wengert的关键洞察是引入"带号表"（tape）的概念，记录计算过程中的所有中间变量及其依赖关系。这个想法直接导致了现代深度学习框架中计算图的概念。虽然Wengert的工作在当时没有立即得到广泛应用（部分因为计算资源限制），但它为后来的自动微分工具（如ADIFOR、TAF）和现代深度学习框架（如TensorFlow、PyTorch）奠定了理论基础。

有趣的是，Wengert最初的动机是解决工程优化问题，特别是化工过程的参数优化。他没有预见到60年后，他的方法会成为人工智能革命的基石。这再次证明了基础研究的长远价值——看似抽象的数学工具可能在未来产生革命性的应用。

高级话题：神经ODE与连续深度模型

连续时间动力学的参数化

神经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）

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$ 是参数化的哈密顿量。这保证了辛结构和长期稳定性。

本章小结

可微分物理仿真通过将传统的前向动力学仿真转化为可计算梯度的系统，为轨迹优化、系统辨识和端到端学习开辟了新的可能。本章的关键要点包括：

核心概念

自动微分：通过计算图和链式法则精确计算梯度，避免了符号微分的表达式膨胀和数值微分的精度问题
隐式函数定理：为约束系统和隐式积分器的梯度计算提供了理论基础
接触光滑化：通过软约束和光滑近似处理不可微的接触现象
梯度流控制：通过裁剪、归一化和截断等技术管理长序列的梯度传播

关键公式

反向模式AD递推：$\bar{\mathbf{x}}_t = (\partial f/\partial \mathbf{x}_t)^T \bar{\mathbf{x}}_{t+1} + (\partial l_t/\partial \mathbf{x}_t)^T$
隐式函数梯度：$d\mathbf{x}^*/d\theta = -(\partial g/\partial \mathbf{x})^{-1}(\partial g/\partial \theta)$
光滑接触力：$\lambda_n = k_n \text{smax}(0, -d)^{1.5} + d_n \text{smax}(0, -\dot{d})$
Fischer-Burmeister函数：$\phi(a,b) = a + b - \sqrt{a^2 + b^2 + \epsilon^2}$

实践要点

混合使用自动微分和解析梯度以平衡效率和实现复杂度
仔细选择光滑化参数 $\epsilon$ 来平衡可微分性和物理真实性
监控梯度范数并采用适当的正则化技术
在长序列仿真中使用梯度截断或检查点技术管理内存

练习题

基础题

练习10.1 手工计算函数 $f(x,y) = e^{xy} + \sin(x+y)$ 在点 $(1,0)$ 处的梯度，分别使用前向和反向模式自动微分。验证两种方法得到相同结果。

提示

构建计算图，识别中间变量，然后应用链式法则。前向模式需要两次遍历（分别对x和y），反向模式只需一次。

答案

计算图：$v_1 = xy$, $v_2 = e^{v_1}$, $v_3 = x+y$, $v_4 = \sin(v_3)$, $f = v_2 + v_4$

在点 $(1,0)$：$v_1 = 0$, $v_2 = 1$, $v_3 = 1$, $v_4 = \sin(1)$

前向模式（对x）：$\dot{v}_1 = y = 0$, $\dot{v}_2 = e^{v_1}\dot{v}_1 = 0$, $\dot{v}_3 = 1$, $\dot{v}_4 = \cos(1)$, $\partial f/\partial x = \cos(1)$

反向模式：$\bar{f} = 1$, $\bar{v}_2 = 1$, $\bar{v}_4 = 1$, $\bar{v}_3 = \cos(1)$, $\bar{v}_1 = e^{v_1} = 1$ $\bar{x} = \bar{v}_1 \cdot y + \bar{v}_3 = 0 + \cos(1) = \cos(1)$ $\bar{y} = \bar{v}_1 \cdot x + \bar{v}_3 = 1 + \cos(1)$

梯度：$\nabla f = [\cos(1), 1+\cos(1)]$

练习10.2 考虑弹簧-质量系统 $m\ddot{x} + c\dot{x} + kx = 0$，使用隐式欧拉方法离散化。推导位置 $x_{n+1}$ 对参数 $k$ 的导数。

提示

先写出隐式欧拉的离散方程，然后应用隐式函数定理。

答案

隐式欧拉：$x_{n+1} = x_n + hv_{n+1}$, $v_{n+1} = v_n - h(c v_{n+1} + k x_{n+1})/m$

整理得：$(1 + hc/m)v_{n+1} + (hk/m)x_{n+1} = v_n$, $x_{n+1} - hv_{n+1} = x_n$

矩阵形式：$\begin{bmatrix} 1+hc/m & hk/m \\ -h & 1 \end{bmatrix} \begin{bmatrix} v_{n+1} \\ x_{n+1} \end{bmatrix} = \begin{bmatrix} v_n \\ x_n \end{bmatrix}$

应用隐式函数定理：$\frac{\partial x_{n+1}}{\partial k} = -\mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial k} \begin{bmatrix} v_{n+1} \\ x_{n+1} \end{bmatrix}$

其中 $\partial \mathbf{A}/\partial k = \begin{bmatrix} 0 & h/m \\ 0 & 0 \end{bmatrix}$

练习10.3 实现一维接触的光滑化模型。给定距离函数 $d(x)$，设计一个处处可微的接触力函数，满足：当 $d > \epsilon$ 时力约为0，当 $d < -\epsilon$ 时力约为 $k|d|$。

提示

可以使用sigmoid函数或双曲正切函数进行光滑过渡。

答案

一种可能的实现： $$f(d) = \frac{k}{2}\left(1 - \tanh\left(\frac{d}{\epsilon}\right)\right) \cdot \max(0, -d - \epsilon/2)$$ 或使用更光滑的版本： $$f(d) = k \cdot \epsilon \cdot \log\left(1 + e^{-d/\epsilon}\right) \cdot \left(1 + \tanh\left(\frac{-d-\epsilon}{\epsilon}\right)\right)/2$$ 这确保了：

$d > \epsilon$ 时，$f \approx 0$
$d < -\epsilon$ 时，$f \approx k|d|$
在 $[-\epsilon, \epsilon]$ 区间内光滑过渡

挑战题

练习10.4 考虑二维刚体碰撞，设计一个可微分的碰撞响应模型。刚体A和B的状态为 $(x_A, y_A, \theta_A, v_{xA}, v_{yA}, \omega_A)$ 和类似的B状态。当检测到碰撞时，如何计算可微分的冲量？

提示

关键是将离散的碰撞事件转化为连续的接触力，考虑使用罚函数方法或基于LCP的光滑化。

答案

可微分碰撞响应的关键步骤：

渗透深度和法向量：$d = \text{signed_distance}(A, B)$, $\mathbf{n} = \nabla d$
相对速度：在接触点 $\mathbf{p}$： $$\mathbf{v}_{\text{rel}} = (\mathbf{v}_A + \omega_A \times \mathbf{r}_A) - (\mathbf{v}_B + \omega_B \times \mathbf{r}_B)$$
光滑接触力： $$f_n = k \cdot \text{smax}(0, -d)^{1.5} + b \cdot \text{smax}(0, -\mathbf{v}_{\text{rel}} \cdot \mathbf{n})$$
摩擦力（光滑库仑模型）： $$\mathbf{f}_t = -\mu f_n \frac{\mathbf{v}_t}{|\mathbf{v}_t| + \epsilon}$$ 其中 $\mathbf{v}_t = \mathbf{v}_{\text{rel}} - (\mathbf{v}_{\text{rel}} \cdot \mathbf{n})\mathbf{n}$
力和力矩： $$\mathbf{F}_A = (f_n \mathbf{n} + \mathbf{f}_t), \quad \boldsymbol{\tau}_A = \mathbf{r}_A \times \mathbf{F}_A$$ 所有操作都是可微的，梯度可以通过自动微分计算。

练习10.5 分析为什么在可微分仿真中，显式积分器比隐式积分器更容易出现梯度爆炸。设计一个自适应的混合积分策略。

提示

考虑雅可比矩阵的谱半径，以及显式和隐式方法的稳定性区域。

答案

分析：

显式欧拉：$\mathbf{x}_{n+1} = \mathbf{x}_n + h\mathbf{f}(\mathbf{x}_n)$，雅可比 $\mathbf{J} = \mathbf{I} + h\mathbf{J}_f$
如果 $\mathbf{J}_f$ 有特征值 $\lambda$，则 $\mathbf{J}$ 的特征值为 $1 + h\lambda$
对于刚性系统，$|\lambda|$ 很大，导致 $|1 + h\lambda| > 1$，梯度指数增长

隐式欧拉：$\mathbf{J} = (\mathbf{I} - h\mathbf{J}_f)^{-1}$，特征值为 $1/(1 - h\lambda)$，当 $h\lambda$ 很大时趋于0，导致梯度消失但不爆炸。

自适应混合策略：

def adaptive_integrator(x, f, h, stiffness_threshold=10):
    J_f = compute_jacobian(f, x)
    eigenvalues = compute_eigenvalues(J_f)
    stiffness = max(abs(eigenvalues)) * h

    if stiffness < stiffness_threshold:
        # 非刚性：使用显式方法（更快）
        return x + h * f(x)
    else:
        # 刚性：使用隐式方法（更稳定）
        return solve_implicit(lambda x_new: x_new - x - h*f(x_new))

或使用IMEX（隐式-显式）分裂： $$\mathbf{x}_{n+1} = \mathbf{x}_n + h[\mathbf{f}_{\text{explicit}}(\mathbf{x}_n) + \mathbf{f}_{\text{implicit}}(\mathbf{x}_{n+1})]$$

练习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}$ 的导数？

提示

使用Fischer-Burmeister函数或其他NCP函数将LCP转化为非线性方程组，然后应用隐式函数定理。

答案

使用Fischer-Burmeister函数 $\phi(a,b) = a + b - \sqrt{a^2 + b^2 + \epsilon^2}$：

转化为方程组： $$\mathbf{F}(\mathbf{z}, \mathbf{q}) = \begin{bmatrix} \phi(z_1, (\mathbf{M}\mathbf{z} + \mathbf{q})_1) \\ \vdots \\ \phi(z_n, (\mathbf{M}\mathbf{z} + \mathbf{q})_n) \end{bmatrix} = 0$$
计算雅可比： $$\frac{\partial \phi}{\partial a} = 1 - \frac{a}{\sqrt{a^2 + b^2 + \epsilon^2}}$$ $$\frac{\partial \phi}{\partial b} = 1 - \frac{b}{\sqrt{a^2 + b^2 + \epsilon^2}}$$
应用隐式函数定理： $$\frac{\partial \mathbf{z}^*}{\partial \mathbf{q}} = -\left(\frac{\partial \mathbf{F}}{\partial \mathbf{z}}\right)^{-1} \frac{\partial \mathbf{F}}{\partial \mathbf{q}}$$ 其中： $$\frac{\partial \mathbf{F}}{\partial \mathbf{z}} = \text{diag}(\partial \phi/\partial a) + \text{diag}(\partial \phi/\partial b) \mathbf{M}$$ $$\frac{\partial \mathbf{F}}{\partial \mathbf{q}} = \text{diag}(\partial \phi/\partial b)$$ 实现时需要注意数值稳定性，特别是在 $a \approx 0$ 或 $b \approx 0$ 附近。

练习10.7（开放性问题）设计一个可微分的绳索仿真系统。绳索由N个质点通过约束连接，需要处理拉伸、弯曲和自碰撞。讨论主要的技术挑战和可能的解决方案。

提示

考虑使用Position Based Dynamics (PBD)或XPBD框架，它们天然适合可微分化。

答案

主要挑战：

约束耦合：绳索的约束高度耦合，导致雅可比矩阵稠密
自碰撞：组合爆炸的碰撞对，不连续的接触图
数值刚性：不可拉伸约束导致的刚性问题

解决方案：

XPBD框架：

# 位置更新
x_new = x + dt * v + dt^2 * f/m
# 约束投影（可微分）
for constraint in constraints:
    lambda_i = solve_constraint(x_new, constraint)
    x_new += M^(-1) * J^T * lambda_i

距离约束的可微分投影： $$\Delta \mathbf{x} = -\frac{C(\mathbf{x})}{\nabla C^T \mathbf{M}^{-1} \nabla C + \alpha/h^2} \mathbf{M}^{-1} \nabla C$$ 其中 $C = |\mathbf{x}_i - \mathbf{x}_j| - L_0$ 是约束函数。
弯曲约束： $$C_{\text{bend}} = \arccos(\mathbf{e}_1 \cdot \mathbf{e}_2) - \theta_0$$ 使用光滑的arccos近似避免奇异性。
自碰撞的空间哈希加速： - 使用可微分的空间哈希函数 - 软化碰撞检测：$w_{ij} = \exp(-d_{ij}^2/\sigma^2)$ - 加权平均碰撞响应
多重网格求解： - 粗网格捕获全局行为 - 细网格处理局部细节 - 可微分的插值算子

这个系统的梯度计算开销很大，建议使用检查点技术和混合精度训练。

练习10.8 证明使用辛积分器的可微分仿真保持梯度的辛结构。这对长期轨迹优化有什么影响？

提示

从哈密顿系统的变分原理出发，证明梯度动力学也是哈密顿的。

答案

证明概要：

设哈密顿系统：$\dot{\mathbf{q}} = \partial H/\partial \mathbf{p}$, $\dot{\mathbf{p}} = -\partial H/\partial \mathbf{q}$

辛积分器保持：$\mathbf{J}^T \mathbf{\Omega} \mathbf{J} = \mathbf{\Omega}$，其中 $\mathbf{J}$ 是雅可比，$\mathbf{\Omega} = \begin{bmatrix} 0 & \mathbf{I} \\ -\mathbf{I} & 0 \end{bmatrix}$

对于参数化轨迹 $(\mathbf{q}(\theta), \mathbf{p}(\theta))$，变分方程： $$\frac{d}{dt}\begin{bmatrix} \delta \mathbf{q} \\ \delta \mathbf{p} \end{bmatrix} = \begin{bmatrix} \partial^2 H/\partial \mathbf{p} \partial \mathbf{q} & \partial^2 H/\partial \mathbf{p}^2 \\ -\partial^2 H/\partial \mathbf{q}^2 & -\partial^2 H/\partial \mathbf{q} \partial \mathbf{p} \end{bmatrix} \begin{bmatrix} \delta \mathbf{q} \\ \delta \mathbf{p} \end{bmatrix}$$

这仍是哈密顿形式，对应的哈密顿量是原哈密顿量的二阶展开。

影响：

长期稳定性：梯度不会人为增长或衰减
守恒量：梯度流保持某些积分不变量
可逆性：前向和反向传播数值上一致
优化景观：避免了梯度的病态累积，优化问题更良态

实践意义：对于长期轨迹优化（如空间任务规划），使用辛积分器可以避免梯度的数值漂移，提高优化的可靠性。

常见陷阱与错误

陷阱1：过度光滑化

问题：为了可微分性使用过大的 $\epsilon$，导致物理行为失真。解决：使用退火策略，训练初期用大 $\epsilon$，逐步减小。

陷阱2：梯度穿越不连续点

问题：即使使用光滑近似，在某些点（如 $v = 0$ 的摩擦力）梯度仍可能不连续。解决：在这些点附近使用次梯度或随机平滑。

陷阱3：内存爆炸

问题：存储整个仿真轨迹用于反向传播消耗大量内存。解决：使用梯度检查点或递归反向传播。

陷阱4：数值精度损失

问题：长链的梯度传播累积浮点误差。解决：使用混合精度训练，关键部分用双精度。

陷阱5：局部最优陷阱

问题：可微分仿真的优化景观复杂，容易陷入局部最优。解决：结合全局优化方法（如CMA-ES）或使用多起点策略。

最佳实践检查清单

设计阶段

[ ] 识别系统中的不可微操作并设计光滑近似
[ ] 评估是否需要解析梯度还是自动微分足够
[ ] 设计梯度流的监控和诊断机制
[ ] 考虑内存与计算的权衡

实现阶段

[ ] 实现梯度检查（有限差分验证）
[ ] 添加梯度裁剪和归一化
[ ] 实现检查点机制管理内存
[ ] 为关键操作编写自定义反向传播

调试阶段

[ ] 监控梯度范数的演化
[ ] 检查数值稳定性（NaN/Inf）
[ ] 验证物理守恒量（能量、动量）
[ ] 比较不同积分器的梯度质量

优化阶段

[ ] 调整光滑化参数 $\epsilon$
[ ] 实验不同的优化器（Adam, L-BFGS等）
[ ] 尝试梯度截断长度的不同选择
[ ] 评估计算图简化的可能性

部署阶段

[ ] 评估实时性能需求
[ ] 考虑固定某些梯度路径加速计算
[ ] 实现自适应精度控制
[ ] 准备梯度调试和可视化工具