本文档属于 Robotics Tutorial 项目，作者：Pengfei Guo，达妙科技。采用 CC BY 4.0 协议，转载请注明出处。

S04 可微分仿真理论：让梯度穿过物理，也知道梯度何时会失真

本章定位：可微分仿真把仿真器从“只能前向滚动的黑盒”变成“可以对状态、控制、物理参数和策略参数求导的计算图”。这让系统辨识、轨迹优化、策略学习和控制器调参可以使用一阶梯度。但机器人不是普通神经网络，接触、摩擦、碰撞、积分器稳定性和长时域链式法则会让梯度爆炸、消失、偏置甚至方向错误。本章的目标不是把可微分仿真神化，而是建立一套可判断、可调试、可落地的理论框架。

前置依赖：S01 MuJoCo 核心引擎与接触模型，S03 MJX/GPU 生态，足式/70 腿足简化模型理论中的离散化与 MPC，复合/30 多模态 MPC 中的接触模式与约束激活，任一自动微分基础章节中的前向模式与反向模式。

关键问题：什么时候应该相信仿真的梯度？什么时候只应该把它当成近似方向？什么时候应该回到 PPO、MPPI、CEM 或有限差分？

---

S04.0 前置自测

📋 答不出 2 题以上时，建议先回到前置章节复习。

#	问题	来源	预期关键词
1	写出离散时间动力学 \(x_{k+1}=F(x_k,u_k,p)\) 中 \(A_k=\partial F/\partial x\) 和 \(B_k=\partial F/\partial u\) 的物理含义。	足式/70 简化模型、MPC 基础	状态灵敏度、输入灵敏度、线性化
2	显式欧拉和隐式欧拉的区别是什么？为什么隐式方法对刚性系统更稳定？	数值积分基础	当前斜率、下一时刻斜率、稳定域
3	接触互补条件 \(\phi(q)\ge 0,\lambda\ge 0,\phi(q)\lambda=0\) 分别表达什么？	S01 MuJoCo 核心、接触力学	不穿透、法向力非负、非接触或有力
4	自动微分、有限差分、解析导数各自最容易出错的地方是什么？	CppAD/Ceres/JAX 基础	计算图、步长误差、推导和实现错误
5	多模态 MPC 中 mode schedule 为什么会改变约束集合？	复合/30 多模态 MPC	支撑/摆动、接触/自由、激活项

自测标准：

正确数量	阅读建议
5/5	可以直接阅读本章，并重点做练习与代码验证。
3-4/5	可以阅读本章，但每遇到线性化、接触或积分器时要回看前置章节。
0-2/5	建议先补离散动力学、接触约束和自动微分，否则后面的梯度偏差会很难理解。

---

S04.0.1 本章目标

学完本章后，你应该能够：

1. **解释**为什么机器人仿真需要求导，以及可微分仿真相对 PPO、有限差分、轨迹优化的优势和边界。
2. **推导**离散时间动力学的单步灵敏度、多步链式法则、伴随递推和长 horizon 梯度病态。
3. **区分**显式积分、半隐式积分、隐式积分在可微分仿真中的稳定性、导数计算方式和工程取舍。
4. **分析**接触非光滑性、互补约束、软接触和平滑硬接触对梯度质量的影响。
5. **选择**自动微分、解析导数、有限差分、隐函数定理和伴随法的适用场景。
6. **设计**梯度检查实验，定位梯度为 NaN、爆炸、消失、符号错误、模式切换错误等常见问题。
7. **判断**什么时候可微分仿真适合系统辨识、控制优化和策略学习，什么时候应退回零阶方法或混合方法。

---

S04.0.2 知识地图

可微分仿真不是一个单独技巧，而是一条从动力学到优化的链路：

连续动力学
  dx/dt = f(x,u,p)
      │
      ▼
数值积分
  x_{k+1}=F(x_k,u_k,p)
      │
      ▼
局部导数
  A_k=∂F/∂x_k, B_k=∂F/∂u_k, P_k=∂F/∂p
      │
      ▼
多步反传
  dL/dθ 经过 x_0 → x_1 → ... → x_N
      │
      ▼
优化用途
  system ID / trajectory optimization / SHAC / differentiable MPC
      │
      ▼
失效来源
  接触切换 / 长 horizon / stiff dynamics / reward 非光滑 / 仿真近似

模块	本章对应小节	核心判断
为什么求导	S04.1	梯度减少采样，但不保证方向正确
离散动力学	S04.2	真正被求导的是积分后的映射 \(F\)
显式/隐式积分	S04.3	稳定性和可微性要同时考虑
接触非光滑	S04.4	模式边界处导数可能不存在
软接触与互补约束	S04.5	光滑性、物理准确度、步长依赖之间权衡
AD vs 解析导数	S04.6	省力不等于可靠，解析不等于便宜
伴随法/反传	S04.7	反向递推避免显式构造大 Jacobian
梯度检查	S04.8	不检查梯度就等于盲目优化
适用边界	S04.9	可微分仿真是工具，不是默认答案

---

S04.1 为什么要对仿真求导 ⭐

这一节解决什么问题：传统仿真只能回答“给定动作会发生什么”，可微分仿真还要回答“为了让结果更好，动作或参数应该往哪个方向改”。

S04.1.1 从黑盒仿真到带导数的物理层

传统仿真器的接口像这样：

\[ x_{k+1}=F(x_k,u_k,p) \]

其中：

符号	含义	例子
\(x_k\)	当前状态	关节角、速度、基座位姿、接触状态
\(u_k\)	当前控制	关节力矩、期望关节位置、肌肉激活
\(p\)	物理参数	质量、惯量、摩擦系数、阻尼、接触刚度
\(F\)	一步仿真映射	碰撞检测、约束求解、积分器的组合

普通仿真只能给出 \(x_{k+1}\)。

可微分仿真还给出：

\[ A_k=\frac{\partial F}{\partial x_k},\quad B_k=\frac{\partial F}{\partial u_k},\quad P_k=\frac{\partial F}{\partial p} \]

这三个矩阵告诉优化器：

导数	工程含义	常见用途
\(A_k\)	当前状态扰动会如何传播到下一步	LQR/iLQR、稳定性分析、状态估计
\(B_k\)	控制输入扰动会如何影响下一步	MPC、轨迹优化、策略梯度
\(P_k\)	物理参数扰动会如何影响下一步	system ID、real-to-sim、域参数学习

类比到腿足简化模型：

在足式/70 中，LIPM 离散化后得到 \(x_{k+1}=A_dx_k+B_du_k\)。

那里 \(A_d\) 和 \(B_d\) 是手推出来的解析矩阵。

可微分仿真做的是同一件事，只是对象从 LIPM 扩展到完整多体动力学、接触、摩擦和执行器模型。

区别在于：

LIPM/MPC 中的导数	可微分仿真中的导数
模型低维，公式短	模型高维，包含积分和接触求解
多数情况下光滑	接触切换处非光滑
可以手推	常需要 AD、隐函数定理或混合方法
主要用于 MPC	可用于 MPC、RL、system ID、设计优化

本质洞察：可微分仿真的核心不是“仿真器能被神经网络训练”，而是“把物理步进变成优化器可以感知的局部线性模型”。只要能得到可信的局部线性模型，就能把控制、辨识和学习放进同一套链式法则里。

S04.1.2 如果不用导数会怎样

没有导数时，优化器只能靠试探。

例如要调一个落脚点 \(u\)，让质心最终位置 \(x_N\) 接近期望值 \(x^\star\)：

\[ L(u)=\frac{1}{2}\|x_N(u)-x^\star\|^2 \]

零阶方法不知道 \(\frac{dL}{du}\)。

它只能采样很多个 \(u_i\)：

u_1 → rollout → L_1
u_2 → rollout → L_2
u_3 → rollout → L_3
...
选择 loss 较小的方向

这在维度低时可行。

当控制变量变成一条轨迹 \(u_0,\ldots,u_{N-1}\) 时，维度会迅速变大。

如果每步 12 个关节力矩，horizon 为 50：

\[ \dim(U)=12\times 50=600 \]

在 600 维空间里靠随机采样寻找下降方向，样本效率很低。

一阶方法则直接利用：

\[ \nabla_U L = \left[ \frac{\partial L}{\partial u_0}, \frac{\partial L}{\partial u_1}, \ldots, \frac{\partial L}{\partial u_{N-1}} \right] \]

每个控制量都能得到“往哪里改”的信息。

这就是可微分仿真的吸引力。

但如果接触导致 \(\nabla_U L\) 有偏差，下降方向可能是错的。

所以本章始终强调两句话：

直觉	更精确的说法
梯度让学习更快	当梯度可信时，一阶优化比零阶采样更省样本
可微分仿真很强	当系统光滑、horizon 适中、目标函数连续时最强
接触也能求导	接触可以被平滑、近似或隐式求导，但硬接触边界处导数不一定存在
有导数就比 PPO 好	一阶梯度可能低方差但有偏，PPO/REINFORCE 可能高方差但估计目标不同

S04.1.3 四类最常见用途

用途一：系统辨识

真实机器人给出一段轨迹：

\[ \tau_{\text{real}}=\{x_0^{r},u_0,x_1^{r},u_1,\ldots,x_N^{r}\} \]

仿真器使用参数 \(p\) 生成预测轨迹：

\[ x_{k+1}^{s}=F(x_k^{s},u_k,p) \]

定义误差：

\[ L(p)=\sum_{k=1}^{N}\|x_k^{s}(p)-x_k^{r}\|_Q^2 \]

可微分仿真给出：

\[ \nabla_pL \]

于是质量、摩擦、阻尼、关节延迟等参数可以通过梯度下降拟合。

这比网格搜索更高效。

但它也更容易过拟合局部轨迹。

如果训练轨迹只包含慢速直线行走，学到的摩擦参数未必能解释跳跃或侧向滑动。

用途二：轨迹优化

给定初始状态 \(x_0\) 和目标状态 \(x^\star\)，优化控制序列：

\[ \min_{u_0,\ldots,u_{N-1}} \sum_{k=0}^{N-1}\ell(x_k,u_k)+\ell_f(x_N) \]

约束：

\[ x_{k+1}=F(x_k,u_k) \]

如果 \(F\) 可微，iLQR、DDP、SQP 和梯度下降都可以使用仿真器导数。

传统轨迹优化常把动力学写成解析方程。

可微分仿真则让“仿真器本身”成为动力学方程。

这对复杂接触和高维机器人很有吸引力。

用途三：策略学习

策略写成：

\[ u_k=\pi_\theta(o_k) \]

rollout 损失：

\[ L(\theta)=\sum_{k=0}^{N-1}\ell(x_k,\pi_\theta(o_k)) \]

可微分仿真允许：

\[ \frac{dL}{d\theta} = \sum_k \frac{\partial L}{\partial x_k} \frac{\partial x_k}{\partial \theta} + \frac{\partial L}{\partial u_k} \frac{\partial u_k}{\partial \theta} \]

这就是“梯度穿过物理步进”的含义。

SHAC、AHAC 等方法把长 horizon 截短，前半段用可微分仿真的解析梯度，后半段用 critic 估计终值。

用途四：控制器参数学习

不是所有可学习对象都必须是神经网络。

MPC 权重、阻抗控制增益、步态相位参数、落脚点修正系数都可以看作参数 \(\theta\)。

例如落脚点控制器：

\[ p_{\text{foot}}=p_0 + k_v(v-v^\star) \]

如果仿真可微，就能直接学习 \(k_v\)：

\[ \frac{dL}{dk_v} = \frac{\partial L}{\partial p_{\text{foot}}} \frac{\partial p_{\text{foot}}}{\partial k_v} \]

这种用法经常比端到端策略更稳健。

因为结构仍然由控制理论给出，学习只负责调少量关键参数。

S04.1.4 可微分仿真与几种常见方法的关系

方法	用不用仿真模型	用不用仿真导数	梯度来源	典型优点	典型风险
PPO/REINFORCE	用	不用	采样回报估计	鲁棒、适合非光滑奖励	方差高、样本多
MPPI/CEM	用	不用	控制序列采样	容易处理非凸和约束	高维控制采样昂贵
有限差分轨迹优化	用	近似用	扰动仿真差值	实现简单	维度高时代价大
解析模型 MPC	用	用	手推动力学导数	快、可解释	模型简化明显
可微分仿真	用	用	AD/解析/隐式导数	样本效率高	接触和长时域梯度病态
混合 actor-critic	用	部分用	短时域仿真梯度 + critic	兼顾效率和稳定性	critic 偏差和仿真梯度偏差叠加

本质洞察：可微分仿真不是 PPO 的简单替代品。PPO 在估计随机目标的零阶梯度，可微分仿真在估计给定仿真模型上的一阶灵敏度。二者回答的问题相近但不相同；当模型非光滑或目标不连续时，一阶灵敏度可能非常快地给出一个错误方向。

S04.1.5 小代码：一维系统中梯度为什么省样本

下面的代码不依赖机器人库，只用一个一维离散系统说明“有梯度”和“靠采样”的区别。

import numpy as np


def rollout(u, x0=0.0, dt=0.1):
    """一维双积分器：x 是位置，v 是速度，u 是加速度序列。"""
    x = x0
    v = 0.0
    for a in u:
        # 半隐式欧拉：先更新速度，再更新位置
        v = v + dt * a
        x = x + dt * v
    return x


def loss(u, target=1.0):
    """终端位置误差。"""
    xN = rollout(u)
    return 0.5 * (xN - target) ** 2


def finite_difference_grad(u, eps=1e-5):
    """中心差分梯度。维度越高，需要的 rollout 次数越多。"""
    g = np.zeros_like(u)
    for i in range(len(u)):
        up = u.copy()
        um = u.copy()
        up[i] += eps
        um[i] -= eps
        g[i] = (loss(up) - loss(um)) / (2.0 * eps)
    return g


def analytic_grad(u, target=1.0, dt=0.1):
    """手推梯度：每个加速度会影响后续所有位置。"""
    xN = rollout(u, dt=dt)
    dL_dxN = xN - target
    g = np.zeros_like(u)
    n = len(u)
    for i in range(n):
        # 第 i 个加速度先影响速度，再影响从 i 到 N-1 的位置
        dxN_dui = dt * dt * (n - i)
        g[i] = dL_dxN * dxN_dui
    return g


u = np.zeros(20)
g_fd = finite_difference_grad(u)
g_an = analytic_grad(u)
print("有限差分前 5 项:", g_fd[:5])
print("解析梯度前 5 项:", g_an[:5])
print("最大误差:", np.max(np.abs(g_fd - g_an)))

这段程序的教学意义：

观察	解释
有限差分需要 \(2n\) 次 rollout	每个维度都要正负扰动一次
解析梯度一次就能得到所有方向	链式法则复用中间结构
维度越高，差距越大	机器人控制轨迹很容易有几百到几千维
一维系统梯度很可靠	因为系统光滑、没有接触切换

S04.1.6 常见陷阱

类型	错误想法	现象	根本原因	正确做法
概念误区	“能求导就一定能优化好”	loss 下降一段后突然变差	导数只描述局部，接触切换后局部模型失效	监控梯度质量，并与有限差分方向对比
编程陷阱	把仿真器输出的全部状态都放进 loss	梯度被无关状态主导	loss 量纲混乱，角度、速度、力混在一起	对每类误差做归一化和权重设计
思维陷阱	用可微分仿真替代所有 RL	接触丰富任务训练不稳定	一阶梯度可能低方差但有偏	对照 PPO/MPPI 基线，按任务选择
工程陷阱	只看最终 reward	训练看似变好但动作不可部署	梯度可能利用仿真漏洞	同时记录约束裕度、能耗、接触冲击

S04.1.7 练习

#	练习	难度
1	用上面的双积分器代码把 horizon 从 20 改到 200，比较有限差分和解析梯度的计算次数。解释为什么高维控制序列中一阶信息更重要。	⭐
2	设计一个 system ID 损失，用来同时拟合质量 \(m\)、粘性阻尼 \(b\) 和库仑摩擦 \(f_c\)。写出哪些轨迹能激发这些参数，哪些轨迹会导致不可辨识。	⭐⭐
3	比较“直接学习神经网络策略”和“学习 MPC 权重”的可解释性差异。要求从安全约束、调试日志和部署风险三个角度回答。	⭐⭐

---

S04.2 离散时间动力学与灵敏度 ⭐⭐⭐

这一节解决什么问题：可微分仿真真正求导的对象不是连续方程 \(\dot{x}=f(x,u)\)，而是经过积分器和接触求解器之后的一步映射 \(x_{k+1}=F(x_k,u_k,p)\)。

S04.2.1 从连续系统到离散映射

机器人动力学常从连续形式开始：

\[ \dot{x}=f(x,u,p) \]

但计算机仿真必须离散化。

采样时间为 \(h\) 时，一步仿真写成：

\[ x_{k+1}=F_h(x_k,u_k,p) \]

下标 \(h\) 提醒我们：

同一个连续动力学 \(f\)，使用不同积分器和不同步长，会得到不同的离散映射 \(F_h\)。

例如显式欧拉：

\[ F_h(x_k,u_k,p)=x_k+h f(x_k,u_k,p) \]

四阶 Runge-Kutta：

\[ \begin{aligned} k_1 &= f(x_k,u_k,p)\\ k_2 &= f(x_k+\frac{h}{2}k_1,u_k,p)\\ k_3 &= f(x_k+\frac{h}{2}k_2,u_k,p)\\ k_4 &= f(x_k+h k_3,u_k,p)\\ x_{k+1} &= x_k+\frac{h}{6}(k_1+2k_2+2k_3+k_4) \end{aligned} \]

这两个 \(F_h\) 的导数不同。

所以可微分仿真中的“动力学导数”必须明确指的是：

离散时间步进函数的导数，而不是连续微分方程右端的导数。

这点很重要。

如果你用 \(\partial f/\partial x\) 近似 \(\partial F/\partial x\)，就默认了显式欧拉和小步长。

在刚性接触或大步长场景下，这个近似可能严重错误。

S04.2.2 单步线性化

对一步映射 \(F\) 在当前轨迹 \((\bar{x}_k,\bar{u}_k)\) 附近做一阶 Taylor 展开：

\[ x_{k+1} \approx F(\bar{x}_k,\bar{u}_k) + A_k(x_k-\bar{x}_k) + B_k(u_k-\bar{u}_k) \]

其中：

\[ A_k=\left.\frac{\partial F}{\partial x}\right|_{\bar{x}_k,\bar{u}_k}, \quad B_k=\left.\frac{\partial F}{\partial u}\right|_{\bar{x}_k,\bar{u}_k} \]

如果还要对物理参数 \(p\) 求导：

\[ P_k=\left.\frac{\partial F}{\partial p}\right|_{\bar{x}_k,\bar{u}_k,\bar{p}} \]

扰动传播写成：

\[ \delta x_{k+1}=A_k\delta x_k+B_k\delta u_k+P_k\delta p \]

这就是所有可微分仿真优化的局部骨架。

矩阵	维度	说明
\(A_k\)	\(n_x\times n_x\)	状态误差如何传播
\(B_k\)	\(n_x\times n_u\)	控制扰动如何进入状态
\(P_k\)	\(n_x\times n_p\)	参数扰动如何改变下一步
\(\delta x_k\)	\(n_x\)	当前状态小扰动
\(\delta u_k\)	\(n_u\)	当前控制小扰动
\(\delta p\)	\(n_p\)	参数小扰动

从 MPC 视角看，\(A_k,B_k\) 是优化器的预测模型。

从策略学习视角看，它们是反向传播中的 Jacobian。

从 system ID 视角看，\(P_k\) 是参数可辨识性的局部证据。

S04.2.3 多步链式法则

考虑两步：

\[ x_1=F_0(x_0,u_0) \]

\[ x_2=F_1(x_1,u_1) \]

若要计算 \(x_2\) 对 \(u_0\) 的导数：

\[ \frac{\partial x_2}{\partial u_0} = \frac{\partial F_1}{\partial x_1} \frac{\partial x_1}{\partial u_0} = A_1B_0 \]

三步：

\[ \frac{\partial x_3}{\partial u_0} = A_2A_1B_0 \]

推广到 \(N\) 步：

\[ \frac{\partial x_N}{\partial u_i} = A_{N-1}A_{N-2}\cdots A_{i+1}B_i \]

这条公式解释了长 horizon 梯度的两个典型病态：

情况	数学表现	工程现象
转移矩阵乘积的主奇异值长期衰减	\(A_{N-1}\cdots A_{i+1}\) 趋近 0	梯度消失，早期动作学不到
转移矩阵乘积的主奇异值长期增长	\(A_{N-1}\cdots A_{i+1}\) 放大扰动	梯度爆炸，更新一步就发散
\(A_k\) 在接触切换处跳变	乘积不稳定	梯度符号忽正忽负
\(B_i\) 接近零	控制不影响状态	参数或动作不可辨识

这和深度神经网络中的梯度消失/爆炸非常像。

区别在于这里每一层不是线性层或激活函数，而是一次物理仿真步。

工程日志里常记录每一步的 \(\|A_k\|\) 或谱半径，但这只是启发式指标。真正决定长时域梯度的是矩阵乘积的方向性和奇异值增长；非正规矩阵甚至可能在单步范数不大的情况下产生短时放大。

因此可微分仿真也常使用神经网络中的工程技巧：

技巧	在可微分仿真中的作用
短 horizon	减少 \(A\) 连乘长度
gradient clipping	限制爆炸梯度
loss normalization	避免某些物理量主导梯度
checkpointing	在内存和重算之间折中
critic terminal value	用值函数替代远期反传
curriculum	先训练平滑简单任务，再进入接触复杂任务

S04.2.4 线性化不是全局模型

线性化只在局部有效。

如果控制更新太大：

\[ u_k^{new}=u_k+\alpha\Delta u_k \]

即使 \(\Delta u_k\) 是下降方向，也可能因为 \(\alpha\) 太大跨过接触模式边界。

一旦跨过模式边界，原来的 \(A_k,B_k\) 就不再描述新的动力学。

这就是为什么 iLQR、DDP、SQP 和梯度策略学习都需要 line search、trust region 或步长控制。

优化器	控制线性化失效的方法
iLQR/DDP	forward rollout line search
SQP	trust region / merit function
PPO	clipped objective 限制策略更新
SHAC/AHAC	短 horizon + 梯度裁剪 + critic
system ID	小学习率 + 参数物理边界

本质洞察：\(A_k,B_k\) 不是“真实世界的公式”，而是“当前轨迹附近的一阶显微镜”。显微镜能看清局部纹理，却不能告诉你跨过接触边界后世界会怎样。

S04.2.5 参数灵敏度与可辨识性

对参数 \(p\) 的多步导数满足递推：

\[ \frac{\partial x_{k+1}}{\partial p} = A_k\frac{\partial x_k}{\partial p} + P_k \]

初始条件通常是：

\[ \frac{\partial x_0}{\partial p}=0 \]

如果初始状态也依赖参数，则要把对应项加入。

这个递推告诉我们：

参数不是“只影响当前一步”。

一个质量参数的改变会影响加速度，进而影响速度，再影响位置，最后影响接触时机。

可辨识性取决于两件事：

条件	含义
参数真的影响轨迹	\(P_k\) 不能长期为零
轨迹把影响放大到可观测量	\(A\) 递推后的误差能在传感器中出现

例如：

参数	需要什么激励	如果缺少激励
质量	加减速、跳跃、外力扰动	仅匀速运动中难以辨识
粘性阻尼	不同速度段	速度范围太窄时与摩擦混淆
库仑摩擦	正反向低速运动	只单向运动时符号不充分
接触刚度	明显接触压缩	无冲击或无接触时不可辨识
关节延迟	快速变化的控制输入	平滑慢动作中难以估计

S04.2.6 C++ 示例：离散双积分器的 \(A,B\) 矩阵

#include <Eigen/Dense>
#include <iostream>

struct DoubleIntegrator {
    double dt;

    // 状态 x = [position, velocity]，输入 u = acceleration
    Eigen::Vector2d step(const Eigen::Vector2d& x, double u) const {
        Eigen::Vector2d next;
        // 半隐式欧拉：速度先变，位置再用新速度更新
        next(1) = x(1) + dt * u;
        next(0) = x(0) + dt * next(1);
        return next;
    }

    void linearize(Eigen::Matrix2d& A, Eigen::Vector2d& B) const {
        // 对 step 函数手推导数：
        // v_next = v + dt*u
        // q_next = q + dt*v_next = q + dt*v + dt^2*u
        A.setIdentity();
        A(0, 1) = dt;
        B << dt * dt, dt;
    }
};

int main() {
    DoubleIntegrator sys{0.01};
    Eigen::Matrix2d A;
    Eigen::Vector2d B;
    sys.linearize(A, B);
    std::cout << "A =\n" << A << "\n";
    std::cout << "B =\n" << B << "\n";
    return 0;
}

这段代码对应的数学式是：

\[ \begin{bmatrix} q_{k+1}\\ v_{k+1} \end{bmatrix} = \begin{bmatrix} 1 & h\\ 0 & 1 \end{bmatrix} \begin{bmatrix} q_k\\ v_k \end{bmatrix} + \begin{bmatrix} h^2\\ h \end{bmatrix} u_k \]

注意这里的 \(B\) 与显式欧拉不同。

显式欧拉若先用旧速度更新位置：

\[ q_{k+1}=q_k+h v_k \]

则：

\[ B= \begin{bmatrix} 0\\h \end{bmatrix} \]

同一个连续系统，不同积分器给出的离散导数不同。

这就是本节最重要的工程结论。

S04.2.7 常见陷阱

类型	错误做法	现象	根本原因	正确做法
概念误区	把 \(\partial f/\partial x\) 当成 \(\partial F/\partial x\)	MPC 预测和仿真不一致	连续导数没有包含积分器	对实际 step 函数求导
编程陷阱	改了积分器但没改解析 Jacobian	梯度检查失败	\(F\) 变了，\(A,B\) 没变	Jacobian 与 step 同步测试
思维陷阱	认为线性化适用于任意大更新	line search 失败、接触模式乱跳	一阶近似只局部成立	使用小步长、trust region
工程陷阱	只检查单步导数	多步训练发散	单步正确不代表链式乘积稳定	同时检查 1、10、50 步梯度

S04.2.8 练习

#	练习	难度
1	手推显式欧拉、半隐式欧拉下双积分器的 \(A,B\)。解释为什么位置对当前输入的导数不同。	⭐
2	对一个 LIPM 离散系统，计算 \(\partial x_N/\partial u_0=A^{N-1}B\)。改变 CoM 高度 \(h\)，观察自然频率 \(\omega=\sqrt{g/h}\) 如何影响梯度大小。	⭐⭐
3	设计一个实验，比较单步梯度正确但 100 步梯度爆炸的情况。要求记录每一步的 \(\|A_k\|_2\) 和累计梯度范数。	⭐⭐⭐

---

S04.3 显式积分、隐式积分与可微性 ⭐⭐⭐

这一节解决什么问题：积分器不只是数值细节。它决定仿真稳定性，也决定导数的形式、代价和可靠性。

S04.3.1 为什么积分器会影响梯度

连续系统：

\[ \dot{x}=f(x,u) \]

积分器给出离散映射：

\[ x_{k+1}=F_h(x_k,u_k) \]

可微分仿真要反传的是 \(F_h\)。

所以积分器改变三件事：

影响	说明
前向轨迹	同一输入下得到的 \(x_{k+1}\) 不同
局部导数	\(A_k=\partial F_h/\partial x_k\) 不同
优化地形	loss 对输入和参数的曲率不同

举一个最小例子：

\[ \dot{x}=a x \]

显式欧拉：

\[ x_{k+1}=x_k+h a x_k=(1+ha)x_k \]

导数：

\[ \frac{\partial x_{k+1}}{\partial x_k}=1+ha \]

隐式欧拉：

\[ x_{k+1}=x_k+h a x_{k+1} \]

整理：

\[ x_{k+1}=\frac{1}{1-ha}x_k \]

导数：

\[ \frac{\partial x_{k+1}}{\partial x_k}=\frac{1}{1-ha} \]

当 \(a<0\) 且系统很 stiff 时，显式欧拉可能因为 \(|1+ha|>1\) 而数值发散。

隐式欧拉则可能保持稳定。

这说明：

前向仿真不稳定时，梯度通常也不可靠。

S04.3.2 显式欧拉：最简单但稳定域小

显式欧拉：

\[ x_{k+1}=x_k+h f(x_k,u_k) \]

导数：

\[ A_k=I+h\frac{\partial f}{\partial x}(x_k,u_k) \]

\[ B_k=h\frac{\partial f}{\partial u}(x_k,u_k) \]

优点：

优点	解释
实现简单	只需要当前状态的斜率
导数简单	直接使用 \(f\) 的 Jacobian
适合教学	公式短，容易做梯度检查

缺点：

缺点	后果
稳定域小	stiff 系统需要很小步长
能量行为差	机械系统可能虚假增能
接触中容易穿透	大步长下错过碰撞事件

显式欧拉适合学习推导。

它不适合作为高刚度接触系统的默认仿真器。

如果你把接触刚度调大，但步长不变，显式欧拉经常产生爆炸梯度。

S04.3.3 半隐式欧拉：机器人仿真常见折中

许多机械系统把状态分成位置 \(q\) 和速度 \(v\)：

\[ \dot{q}=v \]

\[ \dot{v}=a(q,v,u) \]

显式欧拉：

\[ v_{k+1}=v_k+h a(q_k,v_k,u_k) \]

\[ q_{k+1}=q_k+h v_k \]

半隐式欧拉：

\[ v_{k+1}=v_k+h a(q_k,v_k,u_k) \]

\[ q_{k+1}=q_k+h v_{k+1} \]

区别只是一行。

但对机械系统很重要。

半隐式欧拉先更新动量，再用新速度更新位置。

它通常比显式欧拉更好地保持机械系统能量结构。

对双积分器，半隐式欧拉让 \(q_{k+1}\) 直接依赖当前输入 \(u_k\)：

\[ \frac{\partial q_{k+1}}{\partial u_k}=h^2 \]

显式欧拉中这一项是 0。

这不是数学矛盾。

这是离散化选择不同。

可微分仿真中必须接受这个事实：

梯度是离散仿真器的梯度，不是理想连续世界的梯度。

S04.3.4 RK4：高精度但导数图更长

Runge-Kutta 4 阶方法：

\[ \begin{aligned} k_1 &= f(x_k,u_k)\\ k_2 &= f(x_k+\frac{h}{2}k_1,u_k)\\ k_3 &= f(x_k+\frac{h}{2}k_2,u_k)\\ k_4 &= f(x_k+h k_3,u_k)\\ x_{k+1} &= x_k+\frac{h}{6}(k_1+2k_2+2k_3+k_4) \end{aligned} \]

RK4 的导数通过四个中间斜率传播。

如果用自动微分，计算图会自然包含 \(k_1,k_2,k_3,k_4\)。

如果手推 Jacobian，需要逐层套链式法则。

优点：

优点	说明
平滑系统精度高	同步长下截断误差小
适合无接触轨迹优化	飞行器、机械臂自由空间常用
梯度通常更平滑	因为数值误差更小

缺点：

缺点	说明
每步多次动力学计算	前向更贵，反向更贵
不自动解决 stiff 问题	高刚度接触仍可能不稳定
中间状态可能穿透	接触事件仍需要特殊处理

S04.3.5 隐式欧拉：稳定但需要解方程

隐式欧拉：

\[ x_{k+1}=x_k+h f(x_{k+1},u_k) \]

右端使用未知的 \(x_{k+1}\)。

因此每一步都要解非线性方程：

\[ G(x_{k+1},x_k,u_k)=x_{k+1}-x_k-hf(x_{k+1},u_k)=0 \]

这看起来更麻烦。

但它对 stiff 系统更稳定。

刚性弹簧、阻尼接触、约束稳定化都常常需要隐式思想。

隐式步的导数不能只对公式直接读出来。

对 \(G=0\) 两边对 \(x_k\) 求导：

\[ \frac{\partial G}{\partial x_{k+1}} \frac{\partial x_{k+1}}{\partial x_k} + \frac{\partial G}{\partial x_k}=0 \]

于是：

\[ \frac{\partial x_{k+1}}{\partial x_k} = - \left(\frac{\partial G}{\partial x_{k+1}}\right)^{-1} \frac{\partial G}{\partial x_k} \]

对隐式欧拉：

\[ \frac{\partial G}{\partial x_{k+1}} = I-h\frac{\partial f}{\partial x}(x_{k+1},u_k) \]

\[ \frac{\partial G}{\partial x_k}=-I \]

所以：

\[ A_k= \left( I-h f_x(x_{k+1},u_k) \right)^{-1} \]

对输入：

\[ \frac{\partial G}{\partial u_k} = -h f_u(x_{k+1},u_k) \]

\[ B_k= \left( I-h f_x(x_{k+1},u_k) \right)^{-1} h f_u(x_{k+1},u_k) \]

这就是隐式积分可微性的核心。

它不是“反传穿过 Newton 求解的每一次迭代”这一种方式。

更常见、更高效的方式是用隐函数定理直接求解线性系统。

S04.3.6 隐式导数的工程含义

隐式导数中出现矩阵求解：

\[ \left( I-h f_x \right)^{-1}v \]

工程上不应该显式求逆。

应该解线性方程：

\[ \left( I-h f_x \right)y=v \]

原因：

做法	问题
显式求逆	数值不稳定，计算贵
分解后求解	稳定，可复用前向求解的分解
迭代线性求解	大规模稀疏系统更合适

这与可微分 MPC 中 KKT 隐式求导非常相似。

两者都把“求解器的输出”看作隐函数。

前向计算求出解。

反向计算不必展开所有迭代，而是解一个伴随线性系统。

S04.3.7 积分器选型表

积分器	前向代价	稳定性	导数难度	适用场景
显式欧拉	低	差	低	教学、小步长、平滑系统
半隐式欧拉	低	中	低-中	刚体仿真常见折中
RK4	中-高	中	中	无接触或弱接触高精度轨迹
隐式欧拉	高	好	高	stiff、软接触、约束稳定
变分积分器	中-高	好	高	保结构机械系统
事件驱动积分	高	依赖事件处理	高	碰撞时刻需要精确处理

决策流程：

系统是否有强接触或 stiff 弹簧？
  ├─ 否 → 需要高精度？
  │      ├─ 是 → RK4
  │      └─ 否 → 半隐式欧拉
  └─ 是 → 是否能接受每步解方程？
         ├─ 是 → 隐式欧拉或隐式接触求解
         └─ 否 → 降低刚度、减小步长、改用软化模型或短 horizon

S04.3.8 Python 示例：显式与隐式欧拉的梯度差异

import numpy as np


def explicit_step(x, a, h):
    """显式欧拉：x_{k+1}=(1+h*a)*x_k。"""
    return (1.0 + h * a) * x


def implicit_step(x, a, h):
    """隐式欧拉：x_{k+1}=x_k/(1-h*a)。"""
    return x / (1.0 - h * a)


def rollout(step_fn, x0, a, h, n):
    x = x0
    for _ in range(n):
        x = step_fn(x, a, h)
    return x


a = -50.0      # stiff 衰减系统
h = 0.05       # 对显式欧拉来说步长偏大
n = 20

xe = rollout(explicit_step, 1.0, a, h, n)
xi = rollout(implicit_step, 1.0, a, h, n)

# 对线性系统，末状态对初值的梯度就是单步放大因子的 n 次方。
grad_explicit = (1.0 + h * a) ** n
grad_implicit = (1.0 / (1.0 - h * a)) ** n

print("显式欧拉末状态:", xe)
print("隐式欧拉末状态:", xi)
print("显式欧拉 dx_N/dx_0:", grad_explicit)
print("隐式欧拉 dx_N/dx_0:", grad_implicit)

这个例子中连续系统 \(\dot{x}=-50x\) 本应快速衰减。

显式欧拉的放大因子是：

\[ 1+h a=1-2.5=-1.5 \]

绝对值大于 1。

所以数值解会震荡并发散。

隐式欧拉的放大因子是：

\[ \frac{1}{1-ha}=\frac{1}{3.5} \]

绝对值小于 1。

所以数值解稳定衰减。

这说明前向积分器本身就能决定优化是否有意义。

如果前向仿真已经数值不稳定，再精确的自动微分也只是在精确地反传错误的数值轨迹。

S04.3.9 常见陷阱

类型	错误做法	现象	根本原因	正确做法
概念误区	认为高阶积分器一定更适合可微分仿真	RK4 在接触场景仍然发散	高阶精度不等于处理非光滑事件	接触问题优先考虑事件/约束处理
编程陷阱	隐式步反传穿过所有 Newton 迭代	内存大、梯度慢、迭代数变化导致不稳定	计算图包含求解过程而非解映射	使用隐函数定理或自定义 backward
工程陷阱	接触刚度变大但步长不变	穿透变小但梯度爆炸	stiff 时间尺度没有被解析	减小步长、隐式化、降低刚度
思维陷阱	把积分误差当成模型误差	system ID 学出奇怪参数	参数在补偿数值积分偏差	先做步长收敛实验

S04.3.10 练习

#	练习	难度
1	对 \(\dot{x}=ax\)，分别推导显式欧拉、隐式欧拉和中点法的一步导数。画出不同 \(ha\) 下的稳定区域。	⭐⭐
2	用双质量弹簧系统比较显式欧拉和隐式欧拉。逐步增大弹簧刚度 \(k\)，记录前向能量和梯度范数。	⭐⭐⭐
3	思考为什么半隐式欧拉在游戏物理和机器人仿真中常见，但高精度轨迹优化仍常用 RK4 或 collocation。	⭐⭐

---

S04.4 接触非光滑性：梯度最容易失真的地方 ⭐⭐⭐⭐

这一节解决什么问题：机器人可微分仿真最难的不是刚体动力学，而是接触模式切换。接触让状态转移从一个光滑函数变成分段函数，分段边界处的导数可能不存在。

S04.4.1 从一个球撞地面开始

考虑一维小球高度 \(y\) 和速度 \(v\)。

没有接触时：

\[ \dot{y}=v,\quad \dot{v}=-g \]

如果一步之后仍在空中，离散映射是平滑的。

但如果小球撞到地面，需要施加冲量。

最简单的弹性碰撞模型：

\[ v^+ = -e v^- \]

其中：

符号	含义
\(v^-\)	碰撞前速度
\(v^+\)	碰撞后速度
\(e\)	恢复系数，\(0\le e\le 1\)

当初始高度或初始速度稍微变化时，可能出现两种情况：

情况	一步内是否碰撞	映射
A	不碰撞	继续自由落体
B	碰撞	速度反向并缩放

这就是模式切换。

在模式内部，映射可以很光滑。

在模式边界，映射可能不光滑。

边界由“刚好碰到地面”定义：

\[ \phi(q)=0 \]

其中 \(\phi(q)\) 是间隙函数。

S04.4.2 接触映射为什么不是普通函数

硬接触常用互补条件描述：

\[ \phi(q)\ge 0 \]

\[ \lambda \ge 0 \]

\[ \phi(q)\lambda=0 \]

含义是：

条件	物理含义
\(\phi(q)\ge 0\)	物体不能穿透
\(\lambda\ge 0\)	法向接触力只能推开，不能吸住
\(\phi\lambda=0\)	要么有间隙无力，要么接触有力

第三条是非光滑性的核心。

因为它允许两种互斥模式：

非接触模式：phi > 0, lambda = 0
接触模式：  phi = 0, lambda >= 0

模式切换不是一个连续可微的变化。

当 \(\phi\) 从正数接近 0，\(\lambda\) 可以从 0 突然变成正值。

如果还包含摩擦，切向力还会出现 sticking/sliding 模式切换。

接触现象	数学结构	梯度风险
触地/离地	活动约束切换	导数跳变
弹性碰撞	冲量映射	速度不连续
静摩擦/滑动	摩擦锥边界切换	切向导数不连续
接触点替换	最近点不连续	接触 Jacobian 跳变
多接触冗余	力分配不唯一	解映射集合值

S04.4.3 分段光滑与边界项

一个最小数学例子是 Heaviside 函数：

\[ H(z)= \begin{cases} 0,&z<0\\ 1,&z\ge 0 \end{cases} \]

在 \(z\ne 0\) 的地方，导数为 0。

在 \(z=0\) 的地方，普通意义下导数不存在。

如果某个 loss 是：

\[ L(\theta)=H(\theta) \]

那么在每个分段内部看，导数几乎处处为 0。

但从期望角度看，如果 \(\theta\) 控制随机变量跨过边界的概率，真实目标的变化来自边界移动。

这就是 Suh 等人在刚性/不连续系统中讨论的一阶批量梯度偏差的直觉来源。

可微分仿真沿着单条样本轨迹反传时，通常只看到当前分段内部的导数。

它可能忽略“模式边界移动”带来的贡献。

用一句话说：

分段内部的梯度不等于分段函数期望的梯度。

S04.4.4 FoBG 偏差的直觉

考虑目标：

\[ J(\theta)=\mathbb{E}_{\xi}[L(F(\theta,\xi))] \]

一阶批量梯度常做：

\[ \widehat{g}_{\text{FoBG}} = \frac{1}{N}\sum_{i=1}^N \nabla_{\theta}L(F(\theta,\xi_i)) \]

如果 \(F\) 光滑，交换梯度和期望通常比较自然：

\[ \nabla_{\theta}\mathbb{E}[L] \approx \mathbb{E}[\nabla_{\theta}L] \]

但如果 \(F\) 有跳变，边界贡献不能简单忽略。

可微分仿真给出的梯度可能：

表现	含义
方差很低	每条轨迹内部导数稳定
偏差很大	忽略了跨模式边界的概率变化
更新很快	梯度计算效率高
方向错误	快速走向错误解

这也是本章反复强调的判断：

可微分仿真的梯度不是“天然更真”。

它更像是在当前模式内部给出的局部灵敏度。

当优化目标主要由模式切换决定时，局部灵敏度会失去代表性。

S04.4.5 与多模态 MPC 的联系

复合/30 多模态 MPC 中强调：

模式决定动力学、约束和代价。

可微分仿真中的接触切换也是模式问题。

区别是：

多模态 MPC	可微分仿真
mode schedule 常显式给定	接触模式常由仿真器自动决定
优化器知道每段激活哪些约束	AD 可能只看到当前执行路径
切换可加入过渡模式	切换可能表现为非光滑事件
调试看模式表	调试看接触事件和 Jacobian

如果把接触模式完全隐藏在仿真器内部，梯度调试会变难。

一个实用做法是记录每步的接触模式：

时间	接触点数量	接触法向	sticking/sliding	接触力	梯度范数
\(k\)	\(n_c\)	\(\hat{n}_i\)	每个接触点状态	\(\lambda_i\)	\(\|\partial L/\partial x_k\|\)

当梯度突然爆炸或符号翻转时，常能在同一时刻看到接触模式切换。

S04.4.6 Python 示例：软化 Heaviside 观察梯度

import numpy as np


def hard_step(z):
    """硬切换：除了 z=0，数值导数几乎处处为 0。"""
    return 1.0 if z >= 0.0 else 0.0


def soft_step(z, alpha):
    """用 sigmoid 软化切换，alpha 越大越接近硬切换。"""
    return 1.0 / (1.0 + np.exp(-alpha * z))


def soft_step_grad(z, alpha):
    """sigmoid 的解析导数。"""
    s = soft_step(z, alpha)
    return alpha * s * (1.0 - s)


for alpha in [1.0, 10.0, 100.0]:
    zs = np.linspace(-0.1, 0.1, 5)
    grads = [soft_step_grad(z, alpha) for z in zs]
    print("alpha =", alpha, "grads =", grads)

这个例子展示了平滑化的权衡：

\(\alpha\)	函数形状	梯度形状	工程含义
小	很软	梯度宽而小	优化容易，但物理边界模糊
中	适度平滑	梯度集中	常用折中
大	接近硬切换	梯度尖峰	物理更硬，但训练不稳定

软接触模型也有类似权衡。

接触越硬，越接近真实刚性约束。

但梯度越尖锐，数值越难稳定。

S04.4.7 常见陷阱

类型	错误想法	现象	根本原因	正确做法
概念误区	“接触力能算出来，所以接触可微”	梯度在触地瞬间跳变	解存在不代表解映射可微	分析活动集和模式边界
编程陷阱	不记录接触模式	无法解释梯度尖峰	只看状态看不出约束切换	记录接触数量、力、法向和模式
思维陷阱	把接触软化当成完全解决	sim2real 中接触行为变差	软化改变了物理模型	同时评估前向物理和梯度质量
工程陷阱	在高冲击任务上长 horizon 反传	loss 震荡、梯度 NaN	多次碰撞导致 Jacobian 连乘病态	短 horizon、事件分段、混合零阶方法

S04.4.8 练习

#	练习	难度
1	用一维小球撞地模型写出“未碰撞”和“已碰撞”两段映射，分别手推对初始速度的导数。说明边界处为什么不可微。	⭐⭐
2	复现 sigmoid 软化例子，逐步增大 \(\alpha\)，画出函数和梯度。解释这对应接触刚度调大时的什么现象。	⭐⭐
3	设计一个接触日志格式，用于定位四足机器人训练中梯度爆炸的时间点。字段至少包含接触脚、法向力、切向力比例和梯度范数。	⭐⭐⭐

---

S04.5 软接触、互补约束与平滑硬接触 ⭐⭐⭐⭐

这一节解决什么问题：处理接触梯度有多条路线。软接触让函数变光滑，但改变物理；互补约束保留硬接触结构，但求导更难；平滑硬接触试图在二者之间折中。

S04.5.1 penalty-based 软接触

最常见的软接触是弹簧阻尼模型。

设间隙函数为 \(\phi(q)\)。

当 \(\phi<0\) 表示穿透。

穿透深度：

\[ d=\max(0,-\phi(q)) \]

法向力可以写成：

\[ \lambda_n=k d + c \dot{d} \]

并加上非负截断：

\[ \lambda_n=\max(0,kd+c\dot{d}) \]

优点：

优点	说明
实现简单	不需要解互补问题
容易放进 AD	力是状态的显式函数
对学习友好	接触过渡可以变得平滑

缺点：

缺点	后果
允许穿透	物理上不是真正刚性约束
刚度依赖步长	\(k\) 大时需要小步长或隐式积分
参数难调	\(k,c\) 同时影响前向行为和梯度
非负截断仍非光滑	max 处导数不连续

如果把外层非负截断替换为 softplus：

\[ \text{softplus}_\beta(z)=\frac{1}{\beta}\log(1+\exp(\beta z)) \]

则外层力输出变成光滑近似：

\[ \lambda_n=\text{softplus}_\beta(kd+c\dot{d}) \]

但这只平滑了“压缩力不能为负”的外层裁剪。如果穿透深度仍然定义为 \(d=\max(0,-\phi(q))\)，接触开启边界 \(\phi(q)=0\) 仍然非光滑。想让几何接触边界也变得平滑，需要进一步使用 \(d=\text{softplus}_\beta(-\phi(q))\) 这类近似，并重新检查法向速度项、穿透惩罚和物理偏差。\(\beta\) 也会引入新的平滑参数：太小会让接触过软，太大又接近原来的非光滑边界。

S04.5.2 MuJoCo 软接触的直觉

MuJoCo 的接触求解不是简单弹簧模型。

它把约束、接触、摩擦和阻抗式软化放进统一优化框架中。

从可微分仿真的学习角度，重要的是直觉：

参数方向	物理效果	梯度效果
接触更软	允许更明显的变形/穿透	梯度更平滑
接触更硬	更接近刚性接触	梯度更尖锐
阻尼更大	冲击更少	可能抑制高频梯度
步长更小	事件解析更细	代价更高但导数更可靠

这也是为什么 MuJoCo/MJX 常被用于可微分物理实验。

软接触让许多接触场景在工程上可以反传。

但软接触不是对 Suh 2022 结论的否定。

当接触参数趋向硬接触极限时，非光滑问题仍然回来。

S04.5.3 LCP/NCP 互补约束

硬接触可写成互补问题。

最简法向互补：

\[ 0\le \lambda \perp \phi(q)\ge 0 \]

这里的 \(\perp\) 表示：

\[ \lambda^T\phi(q)=0 \]

带摩擦时，问题更复杂。

常见形式包括：

形式	说明
LCP	线性互补问题，常来自线性化接触
NCP	非线性互补问题，保留非线性几何
SOCP	二阶锥摩擦约束
QP	最小化加速度或冲量误差并满足锥近似

互补约束的困难在于活动集。

哪些接触激活，哪些不激活，哪些摩擦方向 stick，哪些 slide，都可能变化。

在固定活动集内部，可以求导。

活动集切换处，导数可能跳变。

S04.5.4 隐函数定理求接触梯度

假设接触求解器输出 \(z^\star\)。

它满足一组残差方程：

\[ R(z^\star,x,u,p)=0 \]

其中 \(z^\star\) 可能包含接触力、冲量、下一步速度或 KKT 乘子。

如果 \(R\) 对 \(z\) 的 Jacobian 非奇异，则隐函数定理给出：

\[ \frac{\partial z^\star}{\partial x} = - \left(\frac{\partial R}{\partial z}\right)^{-1} \frac{\partial R}{\partial x} \]

接着把 \(z^\star\) 代入 \(x_{k+1}=H(x,u,z^\star)\)：

\[ \frac{\partial x_{k+1}}{\partial x} = H_x+H_z\frac{\partial z^\star}{\partial x} \]

这就是 Dojo、可微优化层、可微 MPC 和许多隐式仿真方法背后的共同思想。

它的优势是：

不需要把求解器迭代过程完整展开。

它的限制是：

条件	失败表现
活动集稳定	活动集切换时不满足
Jacobian 非奇异	冗余接触或退化几何可能奇异
残差足够光滑	原始互补条件不光滑
求解器收敛	前向解不准时反向也不准

S04.5.5 解析平滑硬接触

平滑硬接触的目标是：

保留硬接触结构，同时避免活动集切换的尖锐不连续。

常见手段：

手段	思路	代价
barrier	用对数障碍让变量留在可行域内部	障碍参数影响精度
NCP smoothing	把互补函数替换为平滑近似	需要调平滑参数
entropy regularization	让离散/锥选择更平滑	物理解释更弱
interior-point sensitivity	在内点路径上求灵敏度	实现复杂

平滑参数通常记为 \(\mu\) 或 \(\tau\)。

当 \(\mu\) 较大：

效果	说明
梯度更平滑	优化更容易
物理更软	与硬接触差异更大

当 \(\mu\to 0\)：

效果	说明
接近硬接触	前向物理更准确
梯度更尖锐	活动集切换重新出现

所以平滑参数不是越小越好。

它应该随训练阶段、任务冲击强度和部署目标共同选择。

S04.5.6 四类接触梯度策略对比

策略	前向物理	梯度质量	实现难度	推荐场景
penalty 软接触	中	平滑但有模型偏差	低	教学、软体、低冲击学习
硬接触 + 活动集导数	高	模式内准确，边界跳变	高	轨迹优化、低切换接触
平滑互补/内点	高-中	较平滑，可控偏差	高	接触丰富研究、sim2real 前沿
代理模型	前向高，反向近似	依赖代理质量	中	四足 locomotion、工程折中

代理模型的思想是：

前向使用高保真仿真。

反向使用更平滑或更低维的近似模型。

这不是偷懒。

它承认“用于预测的模型”和“用于给出优化方向的模型”可以不同。

在腿足控制中，SRBD 也有类似精神：

真实机器人是高维多体系统。

MPC 用的是简化模型。

只要简化模型给出的方向足够有用，闭环反馈就能补偿剩余误差。

S04.5.7 常见陷阱

类型	错误做法	现象	根本原因	正确做法
概念误区	“软接触越软越好学”	学出的策略依赖穿透	软接触改变可行动作集合	同时限制穿透和接触力峰值
编程陷阱	对 max/clip 不做检查	梯度突然为 0 或 NaN	截断产生非光滑点	使用平滑近似或记录活动区间
工程陷阱	用平滑参数补偿模型错误	参数越调越不物理	平滑不是万能拟合项	分开调前向物理和梯度稳定
思维陷阱	只比较训练速度	sim2real 失败	快速优化可能利用软模型漏洞	加入高保真验证和扰动测试

S04.5.8 练习

#	练习	难度
1	写出 penalty 接触力 \(\lambda=k\max(0,-\phi)\) 对 \(\phi\) 的导数。说明 \(\phi=0\) 处的问题。	⭐
2	将 max 换成 softplus，推导其导数。比较 \(\beta=5,50,500\) 下的接触力曲线。	⭐⭐
3	对一个箱子推墙任务，列出哪些接触适合软化，哪些约束不应过度软化。要求解释物理后果。	⭐⭐⭐

---

S04.6 自动微分 vs 解析导数 vs 有限差分 ⭐⭐⭐

这一节解决什么问题：可微分仿真并不等于“全部交给自动微分”。导数有多种来源，选错会带来性能、内存和正确性问题。

S04.6.1 三种导数来源

方法	核心思想	优点	缺点
有限差分	扰动输入，比较输出差值	最容易实现，适合检查	慢，步长敏感
解析导数	手推公式或符号推导	快，可解释，可利用结构	推导和实现容易错
自动微分	按计算图机械套链式法则	通用，和代码同步	内存大，控制流/非光滑需小心

有限差分公式：

\[ \frac{\partial f}{\partial x_i} \approx \frac{f(x+\epsilon e_i)-f(x-\epsilon e_i)}{2\epsilon} \]

解析导数直接写出：

\[ f(x)=x^2 \Rightarrow f'(x)=2x \]

自动微分则把程序分解成基础操作：

y = x * x
dy/dx = dy/dy * d(x*x)/dx = 1 * 2x

自动微分不是数值差分。

它没有有限差分的截断误差。

但它仍然会忠实地反传程序中的非光滑、分支、clip 和近似。

如果程序中的物理模型错了，AD 会精确地给出错误模型的导数。

S04.6.2 前向模式与反向模式

设函数：

\[ y=f(x),\quad x\in\mathbb{R}^n,\quad y\in\mathbb{R}^m \]

Jacobian：

\[ J=\frac{\partial y}{\partial x}\in\mathbb{R}^{m\times n} \]

前向模式擅长计算：

\[ Jv \]

反向模式擅长计算：

\[ J^T w \]

选择依据：

场景	推荐
输入维度小，输出维度大	前向模式
输入维度大，输出标量 loss	反向模式
MPC 需要完整 \(A,B\)	前向批量或解析导数
策略网络训练	反向模式
system ID 参数少	前向模式很合适
只要方向导数	JVP/VJP 不必显式构造 Jacobian

机器人仿真中常见维度：

对象	维度特点	更自然的模式
物理参数 \(p\)	十几个到几十个	前向模式
控制轨迹 \(U\)	几百到几千	反向模式
状态 Jacobian \(A\)	方阵，可能稀疏	解析/前向批量
标量训练 loss	输出 1 维	反向模式

S04.6.3 解析导数仍然重要

自动微分很方便，但解析导数没有过时。

原因有四个：

原因	说明
稀疏结构	解析导数能直接填 sparse block
物理解释	能看出每一项来自哪个力或约束
性能	可避免构建巨大计算图
稳定性	可对奇异点做专门处理

例如刚体动力学中，质量矩阵、Coriolis 项、接触 Jacobian 都有结构。

Pinocchio、Drake、MuJoCo 这类库不会简单依赖黑盒 AD。

它们通常混合使用：

部分	常见导数方式
运动学	解析 Jacobian
质量矩阵和动力学	解析递推或模板 AD
接触求解	解析/隐式/数值近似
高层 loss	反向模式 AD
策略网络	反向模式 AD

S04.6.4 有限差分的真正价值：检查

有限差分不适合大规模优化。

但它适合做梯度检查。

例如检查方向导数：

\[ g_{\text{AD}}^T v \quad \text{vs} \quad \frac{L(\theta+\epsilon v)-L(\theta-\epsilon v)}{2\epsilon} \]

方向导数检查比逐维检查更便宜。

如果 \(\theta\) 有 10000 维，逐维中心差分要 20000 次 rollout。

方向检查只要 2 次 rollout。

常用流程：

步骤	做法
1	采样一个随机方向 \(v\) 并归一化
2	用 AD 算 \(g=\nabla_\theta L\)
3	计算 \(g^Tv\)
4	用中心差分算方向导数
5	比较相对误差

相对误差：

\[ \text{relerr} = \frac{|d_{\text{AD}}-d_{\text{FD}}|} {\max(1,|d_{\text{AD}}|,|d_{\text{FD}}|)} \]

S04.6.5 常见陷阱

类型	错误做法	现象	根本原因	正确做法
概念误区	“AD 梯度一定正确”	与有限差分不一致	程序包含 detach、分支或非光滑操作	对最小例子做方向导数检查
编程陷阱	在反传路径中原地修改状态	梯度丢失或报错	计算图被破坏	使用函数式 step，避免隐式全局状态
工程陷阱	显式构造巨大 Jacobian	内存爆炸	不需要完整矩阵	使用 JVP/VJP 或伴随递推
思维陷阱	用有限差分当训练梯度	计算极慢	维度太高	有限差分只做校验和小规模基线

S04.6.6 练习

#	练习	难度
1	对 \(f(x)=[x_1x_2,\sin x_1]\) 手推 Jacobian，并用有限差分验证。	⭐
2	对一个 100 维控制序列，只做 5 个随机方向导数检查。解释为什么这比逐维检查更实用。	⭐⭐
3	选择一个接触 loss，说明哪些操作会让 AD 梯度和有限差分差异变大。	⭐⭐⭐

---

S04.7 伴随法与反向传播：长时域梯度如何计算 ⭐⭐⭐⭐

这一节解决什么问题：多步仿真不应该显式构造所有 \(\partial x_i/\partial u_j\)。伴随法用反向递推高效计算 loss 对全部控制和参数的梯度。

S04.7.1 损失函数设定

考虑离散系统：

\[ x_{k+1}=F_k(x_k,u_k,\theta) \]

损失：

\[ L= \ell_N(x_N,\theta) + \sum_{k=0}^{N-1}\ell_k(x_k,u_k,\theta) \]

目标是计算：

\[ \frac{dL}{d\theta} \]

或：

\[ \frac{dL}{du_k} \]

直接展开会产生大量链式项。

伴随法定义：

\[ \lambda_k=\frac{dL}{dx_k} \]

它表示：

从第 \(k\) 步状态开始，后续全部损失对当前状态的敏感度。

S04.7.2 反向递推

终端条件：

\[ \lambda_N=\frac{\partial \ell_N}{\partial x_N} \]

对中间步：

\[ \lambda_k = \frac{\partial \ell_k}{\partial x_k} + A_k^T\lambda_{k+1} \]

其中：

\[ A_k=\frac{\partial F_k}{\partial x_k} \]

控制梯度：

\[ \frac{dL}{du_k} = \frac{\partial \ell_k}{\partial u_k} + B_k^T\lambda_{k+1} \]

参数梯度：

\[ \frac{dL}{d\theta} = \sum_{k=0}^{N-1} \left( \frac{\partial \ell_k}{\partial \theta} + P_k^T\lambda_{k+1} \right) + \frac{\partial \ell_N}{\partial \theta} \]

这三条公式就是反向传播穿过物理仿真的核心。

它们和神经网络反传完全同构：

神经网络	可微分仿真
layer activation	状态 \(x_k\)
layer weights	控制/参数 \(u_k,\theta\)
layer Jacobian	\(A_k,B_k,P_k\)
backprop signal	伴随变量 \(\lambda_k\)

S04.7.3 为什么不显式构造大 Jacobian

如果要把所有状态对所有控制的导数写出来，会得到块下三角矩阵：

\[ \begin{bmatrix} \frac{\partial x_1}{\partial u_0} & 0 & \cdots\\ \frac{\partial x_2}{\partial u_0} & \frac{\partial x_2}{\partial u_1} & \cdots\\ \vdots & \vdots & \ddots \end{bmatrix} \]

块数随 \(N^2\) 增长。

伴随法只需反向扫一遍。

复杂度随 \(N\) 线性增长。

这就是长 horizon 反传的基本可行性来源。

但内存仍然是问题。

反向时需要前向轨迹上的 \(x_k,u_k\) 和可能的中间量。

方法	内存	计算	说明
保存全部中间量	高	低	标准反传
checkpointing	中	中	保存部分，反向时重算
只保存状态	中	高	重新计算局部 Jacobian
伴随 ODE 连续法	低	可能不稳定	神经 ODE 常见，但接触系统需谨慎

S04.7.4 SHAC 的短 horizon 思想

长 horizon 问题：

\[ L=\sum_{k=0}^{H-1}\gamma^k r_k \]

如果 \(H=1000\)，反传需要穿过 1000 个物理步。

接触、摩擦、积分误差和 Jacobian 连乘都会放大病态。

SHAC 的做法是截断：

\[ L_h= \sum_{k=0}^{h-1}\gamma^k r_k + \gamma^h V_\psi(x_h) \]

其中 \(h\) 常远小于 \(H\)。

前 \(h\) 步用可微分仿真反传。

\(h\) 步之后用 critic 估计。

这不是简单近似。

它是在承认：

短时域物理梯度更可信。

长时域价值由统计学习补足。

组件	作用
可微分仿真	提供短时域低方差梯度
critic	估计截断后的长期收益
短 horizon	避免 Jacobian 连乘病态
并行仿真	保持样本覆盖

S04.7.5 反向传播中的 detach 与 stop-gradient

有时故意阻断梯度是正确的。

例如：

场景	为什么阻断
critic target	避免 bootstrap 目标同时追自己
高风险接触事件	避免错误梯度污染前面策略
reset 状态	episode 边界不应反传
随机采样	需要 reparameterization 或 score function
safety clamp	clamp 外梯度可能误导优化

但无意间阻断梯度是常见 bug。

典型表现：

现象	可能原因
loss 变化但梯度全 0	中途 detach 或转 numpy
只有策略最后一层有梯度	观测构造断图
物理参数梯度为 0	参数没有进入 step 函数
接触参数梯度偶尔为 NaN	接触求解失败或除零

S04.7.6 伪代码：短 horizon 可微 rollout

def differentiable_rollout(policy, dynamics_step, x0, horizon):
    """短时域可微 rollout，所有变量保持在计算图中。"""
    xs = [x0]
    us = []
    rewards = []
    x = x0
    for k in range(horizon):
        # 策略根据当前状态输出动作
        u = policy(x)
        # 物理步进必须是可微函数，不能转成普通数组
        x_next = dynamics_step(x, u)
        # 奖励也应尽量保持连续可导
        r = reward_fn(x, u)
        xs.append(x_next)
        us.append(u)
        rewards.append(r)
        x = x_next
    return xs, us, rewards


def shac_like_loss(policy, critic, dynamics_step, x0, horizon, gamma):
    """短 horizon + critic 终值的损失。"""
    xs, us, rewards = differentiable_rollout(policy, dynamics_step, x0, horizon)
    total = 0.0
    discount = 1.0
    for r in rewards:
        total = total + discount * r
        discount = discount * gamma
    # critic 估计截断之后的价值
    total = total + discount * critic(xs[-1])
    # 最大化回报等价于最小化负回报
    return -total

关键检查：

代码位置	检查内容
policy(x)	是否保留计算图
dynamics_step(x,u)	是否包含不可微分支或外部黑盒
reward_fn	是否使用硬阈值奖励
critic(xs[-1])	target 是否需要 stop-gradient
loss backward	梯度范数是否有限

S04.7.7 常见陷阱

类型	错误做法	现象	根本原因	正确做法
概念误区	horizon 越长越好	梯度爆炸或消失	Jacobian 连乘病态	从短 horizon 开始
编程陷阱	rollout 中转成 numpy	前面层无梯度	计算图断开	保持张量类型一致
工程陷阱	不做梯度裁剪	偶发接触导致参数飞掉	接触尖峰放大反传信号	记录并裁剪梯度范数
思维陷阱	critic 只是一项小修补	截断后长期行为差	terminal value 决定远期方向	对 critic 做单独验证

S04.7.8 练习

#	练习	难度
1	从离散系统 \(x_{k+1}=F(x_k,u_k)\) 出发，手推伴随递推 \(\lambda_k=\ell_{x,k}+A_k^T\lambda_{k+1}\)。	⭐⭐
2	用一个线性系统实现短 horizon rollout，比较 horizon=10、50、200 时梯度范数变化。	⭐⭐
3	设计一个 SHAC 风格实验：只改变截断长度 \(h\)，记录收敛速度、梯度范数和最终策略质量。	⭐⭐⭐

---

S04.8 梯度检查与故障定位 ⭐⭐⭐

这一节解决什么问题：可微分仿真中的梯度必须被检查。没有梯度检查，训练失败时很难区分是优化问题、物理问题、代码问题还是接触问题。

S04.8.1 最小梯度检查清单

一个可微分仿真实验至少检查四层：

层级	检查对象	方法
单步 step	\(x_{k+1}=F(x_k,u_k)\)	有限差分对比 AD
多步 rollout	\(x_N(U)\)	方向导数检查
loss	\(L(\theta)\)	AD vs 中心差分
优化方向	更新后 loss 是否下降	小步长 line search

不要一开始就在完整四足策略上查梯度。

应该从最小系统开始。

推荐顺序：

双积分器
  ↓
无接触摆杆
  ↓
带软接触小球
  ↓
单腿接触
  ↓
完整四足或人形

每一层通过后再进入下一层。

S04.8.2 中心差分步长怎么选

中心差分误差有两部分：

误差	来源
截断误差	\(\epsilon\) 太大时 Taylor 近似不准
舍入误差	\(\epsilon\) 太小时浮点消减

经验上 double 精度常从：

\[ \epsilon=10^{-6} \]

或：

\[ \epsilon=10^{-4}\max(1,|\theta_i|) \]

开始。

但接触系统中不能只试一个 \(\epsilon\)。

应该扫描：

\[ \epsilon\in\{10^{-2},10^{-3},10^{-4},10^{-5},10^{-6}\} \]

如果系统光滑，会出现一段相对误差较低的平台。

如果没有平台，常见原因是：

现象	可能原因
小 \(\epsilon\) 误差大	浮点噪声或求解器容差太粗
大 \(\epsilon\) 误差大	非线性或跨接触模式
所有 \(\epsilon\) 都差	AD 断图、解析导数错误、非确定性
正负扰动接触模式不同	点在模式边界附近

S04.8.3 方向导数检查代码

import numpy as np


def directional_check(loss_fn, theta, grad, eps_list):
    """比较 AD 梯度方向导数和中心差分方向导数。"""
    rng = np.random.default_rng(0)
    v = rng.normal(size=theta.shape)
    v = v / (np.linalg.norm(v) + 1e-12)

    ad_dir = float(np.dot(grad.reshape(-1), v.reshape(-1)))
    rows = []

    for eps in eps_list:
        lp = loss_fn(theta + eps * v)
        lm = loss_fn(theta - eps * v)
        fd_dir = float((lp - lm) / (2.0 * eps))
        rel = abs(ad_dir - fd_dir) / max(1.0, abs(ad_dir), abs(fd_dir))
        rows.append((eps, ad_dir, fd_dir, rel))

    return rows


# 使用方式：
# theta 是待检查参数，grad 是自动微分得到的梯度。
# loss_fn 必须是确定性的：同一个 theta 调用多次得到同一个 loss。

检查结果建议打印成表：

\(\epsilon\)	AD 方向导数	FD 方向导数	相对误差	判断
\(10^{-2}\)	...	...	...	可能有非线性
\(10^{-3}\)	...	...	...	观察是否下降
\(10^{-4}\)	...	...	...	常用判断区间
\(10^{-5}\)	...	...	...	观察舍入误差

S04.8.4 梯度健康指标

训练时每轮至少记录：

指标	含义	异常信号
loss	优化目标	突然变 NaN 或震荡
reward 分项	行为来源	某项主导全部梯度
\(\|\nabla_\theta L\|\)	总梯度范数	爆炸或长期为 0
max gradient	局部尖峰	单个参数异常
contact count	接触数量	与梯度尖峰同步变化
penetration	穿透深度	软接触被滥用
solver residual	接触/约束求解精度	前向解不可靠
step size	优化步长	更新过大跨模式

一个实用经验：

把梯度范数、接触数量、最大接触力画在同一张时间轴上。

很多梯度故障会和接触事件对齐。

如果梯度尖峰不对应接触，往往是 reward、归一化或代码断图问题。

S04.8.5 故障排查表

症状	可能原因	排查步骤	修复方向
梯度全 0	detach、clip 饱和、loss 与参数无关	检查每层参数梯度，移除 clip 做最小测试	保持计算图，重写连续 loss
梯度 NaN	除零、sqrt 负数、接触求解失败	打印首个 NaN 时间步和状态	加 epsilon、限制状态、提高求解容差
有限差分不匹配	步长不合适或非光滑	扫描 \(\epsilon\)，记录接触模式	避开边界点或使用平滑近似
单步正确，多步错误	Jacobian 连乘病态	比较 1/10/50 步误差	缩短 horizon、梯度裁剪
loss 下降但行为变差	loss 权重错误或利用仿真漏洞	分项 reward 和物理裕度同图查看	归一化、加入约束和验证场景
system ID 参数不物理	参数不可辨识或补偿数值误差	做激励分析和步长收敛	增加轨迹多样性、加物理边界

S04.8.6 常见陷阱

类型	错误做法	现象	根本原因	正确做法
概念误区	“训练能跑就说明梯度对”	训练到奇怪动作	错误梯度也能降低错误 loss	做方向导数检查
编程陷阱	梯度检查时保留随机噪声	FD 每次不同	loss 非确定性	固定随机种子或关闭随机项
工程陷阱	只在无接触状态检查	接触后失败	检查点覆盖不足	分别检查空中、接触前、接触后
思维陷阱	相对误差必须极小	接触附近误差不可避免	非光滑系统不能按光滑标准判	结合模式一致性解释误差

S04.8.7 练习

#	练习	难度
1	给双积分器写方向导数检查，要求相对误差在 \(10^{-6}\) 量级附近出现平台。	⭐
2	对 softplus 接触模型做同样检查，比较 \(\beta\) 不同时的误差平台。	⭐⭐
3	设计一份训练日志，使它能同时解释 loss、梯度、接触和约束残差之间的因果关系。	⭐⭐⭐

---

S04.9 可微分仿真的适用边界 ⭐⭐⭐

这一节解决什么问题：可微分仿真是强工具，但不是默认答案。判断边界比掌握 API 更重要。

S04.9.1 适合可微分仿真的任务

任务	适合原因	注意事项
system ID	参数维度低，轨迹误差可导	需要丰富激励和物理边界
低接触轨迹优化	动力学较光滑	注意积分器和约束
软体/流体/MPM	本身连续介质更平滑	计算量大
阻抗参数学习	参数少，结构强	不要让软接触漏洞主导
短 horizon 策略学习	梯度链短	需要 critic 或终端代价
可微分 MPC 参数学习	优化层结构明确	KKT/活动集需要平滑处理

S04.9.2 不适合作为首选的任务

任务	风险	更稳妥方案
高频硬碰撞 manipulation	梯度偏差大	PPO/MPPI + 局部可微模块
稀疏奖励探索	梯度没有有效信号	PPO/SAC/CEM、课程学习
长 horizon locomotion 直接反传	梯度爆炸/消失	SHAC/AHAC、PPO、分层训练
接触模式极多的任务	活动集频繁切换	mode schedule、混合优化
前向模型不可信	优化会利用模型误差	先做 system ID 和 sim2sim 验证
强离散逻辑控制	分支不可微	显式模式建模或零阶搜索

S04.9.3 工程决策流程

目标是否需要对物理参数或控制轨迹求精确局部灵敏度？
  ├─ 否 → PPO / MPPI / CEM / 常规 MPC 可能更简单
  └─ 是
      │
      ├─ 系统是否大部分时间光滑或软接触？
      │    ├─ 是 → 可微分仿真值得尝试
      │    └─ 否 → 继续判断
      │
      ├─ 是否能把 horizon 截短并用 critic/终端代价补足？
      │    ├─ 是 → SHAC/AHAC 风格
      │    └─ 否 → 零阶或混合方法更稳妥
      │
      ├─ 是否能记录并控制接触模式？
      │    ├─ 是 → 可做平滑硬接触或代理模型
      │    └─ 否 → 先完善日志和可解释性
      │
      └─ 是否有可靠梯度检查？
           ├─ 是 → 进入训练
           └─ 否 → 先做最小系统验证

S04.9.4 与框架生态的连接

框架/方向	可微机制	适合用途	边界
MJX	JAX 计算图	MuJoCo 生态内的可微物理、system ID、短 rollout	接触梯度仍需检查
Brax	JAX 反向模式	大规模并行、教学和策略学习	物理细节与 MuJoCo 不完全一致
Drake AutoDiffXd	C++ 前向模式	轨迹优化、MPC、解析模型验证	大规模反向训练不如 JAX 方便
DiffTaichi	source-to-source AD	MPM、软体、流体教学	刚体接触生态较弱
Dojo 类隐式引擎	隐函数定理	接触丰富研究	工程生态和部署成本
NVIDIA Warp	tape-based AD	自定义 CUDA 物理和渲染	需要手动管理内核和内存
Genesis	内置可微分仿真接口	多物理统一场景的可微原型和教学	具体机器人任务的梯度质量需逐项验证

框架选择不应只看“是否可微”。

还要看：

问题	判断
前向物理是否可信	仿真器首先要像真实系统
梯度是否可检查	没有检查就不能用于控制结论
是否支持批量并行	策略学习通常需要大量环境
是否能记录接触细节	接触任务必须可调试
是否与部署链路一致	MJCF/URDF/执行器模型要能迁移

S04.9.5 章末综合练习

#	练习	难度
1	用 MJX 或任一可微环境实现 cart-pole 的 \(\partial x_{t+1}/\partial u_t\) 检查，并与中心差分对比。	⭐⭐
2	用软接触小球复现“接触越硬，梯度越尖”的现象。记录穿透深度、法向力、梯度范数。	⭐⭐⭐
3	设计一个简化 SHAC：短 horizon 为 16、32、64 三档，比较梯度范数和训练稳定性。	⭐⭐⭐
4	做一个可微 system ID：用已知质量和摩擦生成轨迹，再从错误参数开始优化。说明哪些参数能恢复，哪些参数互相混淆。	⭐⭐⭐
5	跨章综合题：结合足式/70 的 SRBD、复合/30 的模式表和本章的接触梯度，设计“四足臂推门”的可微分仿真实验。要求说明哪些部分用可微梯度，哪些部分用 mode schedule 或零阶搜索。	⭐⭐⭐⭐

---

本章小结

核心公式速查

主题	公式	含义
离散动力学	\(x_{k+1}=F(x_k,u_k,p)\)	真正被求导的一步仿真映射
单步线性化	\(\delta x_{k+1}=A_k\delta x_k+B_k\delta u_k+P_k\delta p\)	局部扰动传播
多步控制导数	\(\partial x_N/\partial u_i=A_{N-1}\cdots A_{i+1}B_i\)	长 horizon 梯度病态来源
隐式步残差	\(G(x_{k+1},x_k,u_k)=0\)	隐式积分或接触求解
隐函数导数	\(dz^\star/dx=-(R_z)^{-1}R_x\)	不展开求解器迭代
伴随递推	\(\lambda_k=\ell_{x,k}+A_k^T\lambda_{k+1}\)	反向传播穿过物理
控制梯度	\(dL/du_k=\ell_{u,k}+B_k^T\lambda_{k+1}\)	每步控制的优化方向
互补约束	\(0\le\lambda\perp\phi(q)\ge0\)	硬接触数学表达
方向导数检查	\(g^Tv\approx[L(\theta+\epsilon v)-L(\theta-\epsilon v)]/(2\epsilon)\)	高维梯度检查

一句话记忆

概念	一句话
可微分仿真	把物理步进变成可反传的局部线性模型
离散导数	求的是积分器和求解器组合后的 \(F\)，不是只求连续 \(f\)
显式积分	简单便宜，但 stiff 接触下容易炸
隐式积分	稳定性更好，但导数要解线性系统
接触非光滑	模式内部可微，模式边界可能不可微
软接触	用物理偏差换梯度平滑
互补约束	用数学约束表达“不穿透和只推不拉”
自动微分	精确反传程序，但不保证程序物理正确
伴随法	用反向递推避免构造巨大 Jacobian
梯度检查	可微分仿真从研究玩具变成工程工具的门槛

方法选型总结

场景	推荐方法
光滑系统短 horizon	可微分仿真直接反传
长 horizon locomotion	SHAC/AHAC 或 PPO + 可微模块
高冲击硬接触	零阶方法、mode schedule、平滑硬接触研究路线
system ID	可微分仿真 + 多激励轨迹 + 参数边界
MPC 参数学习	可微分 MPC 或可微仿真外层调参
部署前验证	高保真非可微仿真、sim2sim、真实日志回放

---

累积项目：Mini-DiffSim 梯度实验台

本章新增一个可复用实验台，用于后续 S05 可微分 MPC 和策略学习实验。

阶段 1：光滑动力学基线

任务	验收标准
实现双积分器和摆杆 step	单步 AD 与有限差分相对误差有稳定平台
实现显式、半隐式、RK4 三种积分器	输出不同离散导数并能解释差异
记录多步梯度范数	能观察 horizon 增长带来的变化

阶段 2：软接触模块

任务	验收标准
实现 penalty 接触和 softplus 接触	接触力曲线可视化
扫描刚度和平滑参数	穿透、力峰值、梯度范数同时记录
做方向导数检查	接触边界附近能解释误差

阶段 3：短 horizon 策略优化

任务	验收标准
实现短 rollout loss	horizon=8/16/32 可切换
加入梯度裁剪	接触尖峰不导致参数 NaN
加入 terminal value 接口	为 S05 的可微 MPC/critic 连接预留

阶段 4：调试面板

记录项	用途
loss 分项	判断优化目标是否被单项主导
梯度范数	识别爆炸和消失
接触数量	对齐模式切换
穿透深度	发现软接触漏洞
接触力峰值	评估冲击和稳定性
有限差分误差	验证梯度可信度

🔧 故障排查手册

症状	可能原因	排查步骤	相关小节
单步梯度与有限差分不一致	step 中有分支、clip、积分器导数写错	固定状态，扫描 \(\epsilon\)，打印 \(A,B\)	S04.2, S04.8
多步梯度爆炸	\(A_k\) 连乘范数过大或接触尖峰	记录每步 \(\|A_k\|\)、接触力、梯度范数	S04.2, S04.4
梯度消失	horizon 太长、reward 太远、控制影响弱	缩短 horizon，加入中间 loss，检查 \(B_k\)	S04.2, S04.7
训练利用穿透	软接触过软或 loss 没惩罚穿透	记录最大穿透深度，增加物理裕度项	S04.5
system ID 得到不物理参数	激励不足或参数互相补偿	做参数敏感度分析，增加物理边界	S04.2, S04.9
loss 下降但策略不稳	梯度优化了错误代理目标	做高保真验证，检查接触和能耗指标	S04.1, S04.9
接触瞬间 NaN	除零、奇异法向、求解器残差大	打印首个 NaN 状态，降低步长或加正则	S04.4, S04.8
隐式导数很慢	展开了求解器迭代或显式求逆	改用线性求解和自定义 backward	S04.3, S04.5

---

章末回看：应该形成的判断力

可微分仿真的第一层能力，是会写：

\[ x_{k+1}=F(x_k,u_k,p) \]

并对它求导。

第二层能力，是知道这个导数来自积分器、接触求解器、平滑参数和计算图，而不是来自理想物理世界本身。

第三层能力，是能在训练失败时把问题拆开：

问题类别	典型证据
模型问题	前向轨迹与真实或高保真仿真不一致
积分问题	步长变化导致结果大变
接触问题	梯度尖峰与接触模式切换对齐
优化问题	小步长下降，大步长发散
代码问题	AD 与有限差分在光滑点不一致
目标问题	loss 下降但物理裕度恶化

最终要形成的工程判断是：

可微分仿真适合用来给出短时域、局部、结构化的优化方向。

它不适合被当成“所有物理学习问题的默认答案”。

当系统光滑、目标连续、horizon 适中、梯度经过检查时，它能显著提高样本效率。

当任务由硬接触切换、稀疏奖励和长时域探索主导时，零阶方法、显式模式建模、代理模型和高保真验证仍然不可替代。

⚠️ 章末陷阱：把可导当成梯度可信

一段程序能被自动微分系统反传，并不意味着梯度适合控制。接触分支、裁剪、投影和求解器容差都会让梯度反映局部数值路径，而不是稳定的物理趋势。

⚠️ 章末陷阱：只在无接触轨迹上做梯度检查

无接触轨迹通常很光滑，有限差分容易通过。但机器人控制的难点往往出现在接触建立、脱离和滑动切换附近。梯度检查必须覆盖这些边界点，才有教学价值。

⚠️ 章末陷阱：用长时域梯度解释短时域模型

可微分仿真在短 horizon 内能提供清晰灵敏度；horizon 过长时，接触模式切换和数值误差会不断放大。把长时域失败全部归因于模型参数，往往会掩盖积分器、目标函数和探索不足的问题。

练习：验证一个梯度是否值得相信

A 型：构造一维弹簧阻尼系统，写出解析导数、自动微分导数和有限差分导数，比较三者在不同步长下的误差。

B 型：构造一个带 max(0, x) 接触惩罚的离散系统，分别在接触边界两侧做有限差分。解释为什么边界附近的梯度会对扰动方向敏感。

C 型：选择一个短 horizon 四足落脚任务，记录每步接触数量、最大穿透、梯度范数和 loss 分项。说明梯度尖峰是否与接触切换对齐。

---

延伸阅读

#	资料	难度	推荐理由
1	Suh, Simchowitz, Zhang, Tedrake, "Do Differentiable Simulators Give Better Policy Gradients?", ICML 2022	⭐⭐⭐⭐	理解刚性/不连续系统下一阶梯度偏差的核心论文
2	Xu, Makoviychuk 等, "Accelerated Policy Learning with Parallel Differentiable Simulation", ICLR 2022	⭐⭐⭐	SHAC 的主要来源，短 horizon + critic 的代表
3	Freeman 等, "Brax: A Differentiable Physics Engine for Large Scale Rigid Body Simulation", 2021	⭐⭐	JAX 可微物理和大规模并行的代表
4	Hu 等, "DiffTaichi: Differentiable Programming for Physical Simulation", ICLR 2020	⭐⭐⭐	source-to-source AD 在物理仿真中的经典工作
5	de Avila Belbute-Peres 等, "End-to-End Differentiable Physics for Learning and Control", NeurIPS 2018	⭐⭐⭐	把 LCP 作为可微层嵌入学习和控制的早期工作
6	Howell 等, "Dojo: A Differentiable Physics Engine for Robotics", 2022	⭐⭐⭐⭐	隐式接触与平滑互补的研究路线
7	Song, Kim, Scaramuzza, "Learning Quadruped Locomotion Using Differentiable Simulation", CoRL 2024	⭐⭐⭐	四足可微仿真 sim-to-real 的重要案例
8	MuJoCo / MJX 官方文档	⭐⭐	了解 MuJoCo 软接触、MJX 和 JAX 工作流