30_时步法与数值积分

接触力学 · 专题 3.0：时步法与数值积分——非光滑接触动力学的时间离散¶

档位：核心档位 3（博士入学）＋进阶档位 4（博士毕业＋） 建议时长：档位 3 约 35–45 h；档位 4 额外 20–30 h 前置：本目录 10_单边约束与互补问题（LCP/NCP，谁是未知量）、20_摩擦锥与最大耗散（可行集形状）；第二批凸分析（法锥/切锥/近端算子）；第四批刚体动力学（$M(q)$、Jacobian、Featherstone）；常微分方程数值解（Runge-Kutta、隐式 Euler）下游：本目录 40_可微接触（解析雅可比从何而来）、混合系统、接触隐式 MPC（CI-MPC）

底线先行：机器人接触仿真本质上**不是**"光滑 ODE + 碰撞修正"，而是 测度微分包含（measure differential inclusion, MDI） 的时间离散。掌握 Moreau-Jean θ-scheme、Stewart-Trinkle LCP、Anitescu 凸松弛、MuJoCo 软接触求解器这四条主线，才能读懂任何现代仿真器；辛积分、事件驱动法在此**失效**，冲击律是时步的**自然产物**而非外挂。

前置自测¶

📋 前置自测（答不出 ≥ 2 题 → 先回本目录 10/20 与第二批凸分析复习）

线性互补问题（LCP） $0\le z\perp Mz+q\ge 0$ 的含义是什么？什么叫"互补松弛"？Lemke 算法在什么条件下保证终止？
什么是**摩擦锥**（Coulomb friction cone）？二阶锥（second-order cone, SOC）形式与多面体（polyhedral/pyramidal）近似各自长什么样？最大耗散原理（maximal dissipation）如何刻画摩擦力方向？
写出**半隐式 Euler（symplectic Euler）** 与**隐式 Euler（backward Euler）** 对光滑 ODE $\dot v=f(q,v),\ \dot q=v$ 的更新公式。它们在能量行为上有何区别？
什么叫**法锥（normal cone）** $N_C(x)$ 与**切锥（tangent cone）** $T_C(x)$？对凸集 $C=\{x:\phi(x)\ge 0\}$，在边界点 $\phi(x)=0$ 处法锥是什么方向？
一个质点在重力下从高处落到地面、以恢复系数 $e=0.8$ 反复弹跳，弹跳次数与总时间各是有限还是无限？这个现象叫什么？

自测说明：第 1、2 题考"原子操作"（每步要解的子问题），第 3 题考积分器基础，第 4 题考凸分析语言，第 5 题考 Zeno 现象——这五块拼起来正是本专题的全部前置。若第 5 题答"无限次弹跳但有限时间结束"，恭喜，你已经直觉地理解了为什么事件驱动法会崩溃。

本章目标¶

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

**从第一性原理说清**为什么接触动力学不能当"光滑 ODE + 碰撞检测"来积分——解是 BV 函数、力是测度、约束集随时间移动这三重困难各意味着什么；
独立推导 Moreau-Jean θ-scheme 的离散互补条件，解释 $\theta=\tfrac12$ 能量中性、$\theta=1$ 强耗散的能量估计；
手写 Stewart-Trinkle 时步格式的 LCP 矩阵，复现 Anitescu-Potra "copositive-plus ⇒ 每步可解"的证明骨架；
辨析 Anitescu 凸松弛与 Stewart-Trinkle 多面体锥的正交关系，理解凸化引入的 $O(h)$ 人工滑移（gliding/creep）从何而来；
读懂 MuJoCo 软接触求解器：写出它求解的凸优化问题（Gauss 原理 + 正则化）、solref/solimp 的物理含义、PGS/CG/Newton 三个求解器各自的收敛性质；
比较事件驱动 vs 时步法在 Zeno、Painlevé、多体同时冲击下的行为差异，知道为什么 MuJoCo/Drake/Bullet 全是 time-stepping；
打通本专题与"可微接触"（本目录 40）的桥：凸性 ⇒ 解析雅可比 ⇒ RL/MPC 梯度。

本章知识导航¶

本章围绕一个核心问题展开：给定一个会发生接触、碰撞、摩擦的多体系统，如何把它的运动方程交给计算机一步步算出来？ 这个问题表面上属于"数值积分"，但接触的非光滑性使它成为一个独立的数学难题。

本章包含 9 个核心知识点，它们的逻辑关系如下：

§0 为什么时间离散是独立难题（三重困难：BV解/测度力/移动约束）
        │
        ├──→ §1 两条方法论路线：事件驱动 vs 时步法 ⭐⭐
        │         （路线之争——本章其余部分都站在"时步法"一侧）
        │
        ├──→ §2 Moreau 扫动过程与测度微分包含 ⭐⭐⭐
        │         （连续时间的"正确方程"——所有离散格式都要收敛到它）
        │              │
        │              ├──→ §3 Moreau-Jean θ-scheme（速度级离散）⭐⭐⭐
        │              │
        │              └──→ §4 Stewart-Trinkle LCP（位置级离散）⭐⭐⭐
        │                       │
        │                       └──→ §5 Anitescu 凸松弛（LCP→凸QP/CCP）⭐⭐⭐
        │                                │
        │                                └──→ §6 MuJoCo 软接触（连续时间+软化）⭐⭐⭐
        │
        ├──→ §7 冲击律的离散处理（Newton/Poisson/Moreau/Stronge）⭐⭐
        │
        ├──→ §8 半隐式 vs 全隐式积分、辛性与接触的不相容 ⭐⭐⭐
        │
        └──→ §9 穿透、能量与漂移：时步法的"数值病理"与稳定化 ⭐⭐

推荐阅读路径：

理论主干（必学）：§0 → §1 → §2 → §3 → §4 → §5。这条线从"为什么难"出发，建立连续时间的 MDI，再分速度级（Moreau-Jean）和位置级（Stewart-Trinkle）两种离散，最后用凸松弛打通到现代主流。
工程落地（强烈推荐）：§6（MuJoCo）→ §9（穿透/能量/漂移）。这条线告诉你打开任意仿真器源码时，每一行对应哪条数学路线。
数学深化（档位 4）：§7（冲击律）→ §8（辛性不相容）→ 章末"收敛性严格证明"。

导航说明：本路线图只展示**结构**。§0–§2 是"为什么"与"正确方程"，§3–§6 是"四种离散格式"，§7–§9 是"细节与病理"。一个常见误区是把这九节当成九个并列工具——其实它们是一棵树：根是 MDI（§2），主干是四种离散，枝叶是冲击律与稳定化。

前置知识桥接¶

本章站在三块前置知识的肩膀上，这里把它们各用 2–3 行重新激活，免得你翻回去：

回顾本目录 10_单边约束与互补问题：LCP $0\le z\perp w=Mz+q\ge 0$ 的解 $z$ 满足"要么 $z_i=0$，要么 $w_i=0$"。我们在那里学到 Lemke 算法、copositive 矩阵类、以及"互补"如何编码"接触要么张开（$\phi>0,\lambda=0$）要么压紧（$\phi=0,\lambda>0$）"。本章每一步时步迭代，本质上都是在解一个 LCP/NCP/凸 QP——它就是时步法的"原子操作"。
回顾本目录 20_摩擦锥与最大耗散：Coulomb 摩擦把切向力约束在锥 $\|\beta\|\le\mu\lambda_n$ 内；最大耗散原理说滑动时摩擦力取在锥边界、方向与滑动速度反平行。我们在那里学到二阶锥 vs 多面体锥的取舍。本章 §4（Stewart-Trinkle）走多面体锥 → LCP，§5（Anitescu）/§6（MuJoCo）走二阶锥 → 凸锥规划——摩擦锥的形状选择直接决定离散格式的代数类型。
回顾第二批凸分析：凸集 $C$ 在点 $x$ 的法锥 $N_C(x)=\{g:\langle g,y-x\rangle\le 0,\forall y\in C\}$，切锥 $T_C(x)$ 是 $C$ 在 $x$ 的局部线性化。近端算子 $\mathrm{prox}_{f}(v)=\arg\min_x\{f(x)+\tfrac12\|x-v\|^2\}$，Moreau-Yosida 正则化把不可微的指示函数磨光。本章 §2 的扫动过程用法锥/切锥写动力学，§6 的 ADMM/proximal 用近端算子——凸分析是非光滑动力学的"母语"。

如果跳过本章会怎样¶

不学本章，你会在以下两个具体场景栽跟头：

你用 RK4 写了一个 bouncing ball 仿真，球穿过地板飞了出去，或者在接触点附近积分器报"步长缩到机器精度仍不收敛"。 你以为是 bug，其实是数学必然——接触点处速度有跳跃，$f$ 不 Lipschitz，任何基于 $f\in C^k$ 的高阶积分器（RK、Adams 多步、辛 Verlet）在冲击点阶数塌陷到零。不学 §0–§1 你永远不知道该改用 time-stepping。
你想给 RL 训练环境做可微仿真（diff-sim），却发现 MuJoCo 的接触梯度"时有时无"、Bullet 根本拿不到梯度。 不学 §5–§6 你不会明白：可微性来自**凸性**——把接触写成凸优化（Anitescu/SAP/ADMM），解对参数才有良定义的雅可比；而 LCP 的解在激活集切换处不可微。本目录 40_可微接触的全部理论基础都在本章。
你调 MuJoCo/Drake 参数，物体要么穿地板、要么疯狂抖动，反复试错却不知道该动哪个旋钮。 不学 §6、§8、§9 你不会知道：穿透/jitter/漂移各有不同的数学根源（软度、步长-刚度失配、求解器未收敛、约束级别），稳定步长有一个硬界 $h\le 2/\omega_0$（§8），而 solref/solimp 是"穿透-稳定性"的权衡而非"越硬越好"。调参从试错变对症，全靠本章的数学。

预计阅读时间¶

阅读方式	时间	适合谁
精读（含推导与练习）	12–16 h	要做仿真器、可微接触、CI-MPC 研究的读者
速读（跳过收敛证明细节）	5–7 h	用现成仿真器、想懂"为什么这样设参数"的工程师
速查（只看表格 + 速查卡 + 故障排查）	40 min	调参踩坑时回来对症

§0 为什么"时间离散"是独立的数学难题 ⭐⭐¶

动机：一个看似简单的问题¶

考虑最朴素的机器人仿真任务：一个刚体在重力下下落，撞到地面，然后静止（或反弹）。连续时间的牛顿方程是

\[ M\dot v = F_{\text{ext}}(q,v) + J_n^\top \lambda_n,\qquad \dot q = v, \]

其中 $M$ 是质量矩阵，$F_{\text{ext}}$ 是重力等外力，$\lambda_n\ge 0$ 是地面对物体的法向接触力，$J_n=\nabla\phi$ 是 gap 函数 $\phi(q)\ge 0$（物体到地面的有符号距离）的梯度。接触力满足互补条件

\[ 0\le \phi(q)\ \perp\ \lambda_n\ge 0, \]

即"距离为正则无力，有力则距离为零"。

一个初学者的自然想法是：这不就是带约束的 ODE 吗？拿个高阶积分器（RK4）一积不就完了？ 本节要说清，为什么这个想法**根本错误**，错在哪三个层次，以及为什么需要一整套独立的数学工具。

如果用经典 ODE 积分器会怎样¶

我们先做一个思想实验，把 RK4 用在 bouncing ball 上看看会发生什么。

设球初速向下，$t_*$ 时刻 $\phi=0$（触地）。在 $t_*$ 之前，$\lambda_n=0$，方程就是自由落体 $\dot v=-g$，光滑，RK4 完美。在 $t_*$ 这一瞬间，恢复系数 $e$ 的 Newton 冲击律要求 $v^+=-e\,v^-$——速度**瞬间跳变**。

问题来了：

RK4 假设右端 $f$ 在一个步长内 $C^4$ 光滑，它在 $[t_\ell,t_{\ell+1}]$ 内取 4 个采样点做加权。但若 $t_*\in(t_\ell,t_{\ell+1})$，这一步跨越了速度跳变，$f$ 在步内不连续，RK4 的 Taylor 展开前提崩溃，局部截断误差从 $O(h^5)$ 塌陷到 $O(h)$ 甚至 $O(1)$——阶数全部丢失。
即使你用**事件检测**（root-finding 找 $\phi(q(t))=0$ 的精确 $t_*$），在 $t_*$ 处停下、应用冲击律、再重启 RK4——你会遇到更致命的问题：Zeno 现象。

本质洞察：接触动力学的"难"不在于"有约束"——纯等式约束的 DAE（微分代数方程）有成熟的高阶积分器（如 SHAKE/RATTLE）。难在于**约束是单边的（不等式）且会突然激活/失活**，这使得解的正则性（regularity）从 $C^\infty$ 直接掉到 BV（有界变差），而所有经典数值分析的误差估计都建立在解充分光滑的前提上。一旦解不光滑，"阶数"这个概念本身就失去意义。

Zeno 现象：事件驱动法的"死穴"¶

bouncing ball 在 $e<1$ 时，每次弹跳高度按 $e^2$ 衰减，弹跳间隔按 $e$ 衰减。设第一次弹跳后腾空时间为 $t_0$，则后续腾空时间依次为 $e t_0,\ e^2 t_0,\dots$，总时间

\[ T_\infty = t_0\sum_{k=0}^\infty e^k = \frac{t_0}{1-e} < \infty. \]

也就是说，在有限时间 $T_\infty$ 内发生了无穷多次弹跳事件。这叫 Zeno 现象（得名于芝诺的"阿基里斯追龟"悖论）。

对事件驱动法这是**绝症**：它要逐个检测、逐个处理每次碰撞事件，但事件有无穷多个，且在 $T_\infty$ 附近事件间隔 $\to 0$，root-finding 永远停不下来——仿真"卡死"在 $T_\infty$ 之前。现实中工程师常打补丁（"速度小于阈值就当静止"），但这是 ad hoc 的、破坏物理一致性的。

时步法如何化解？ 它**根本不检测事件**。每一步固定步长 $h$，把"这一步内是否发生接触、发生几次"全部交给一个互补问题/优化问题去解。当步长 $h$ 跨过 $T_\infty$ 时，球早已静止，多次微弹跳被这一步"自动吸收"成一个静止解——没有任何特殊处理。这是时步法相对事件驱动法的**第一个根本优势**。

基础铺垫：gap 函数与接触 Jacobian¶

在深入三重困难之前，必须把两个贯穿全章的基础对象讲清——gap 函数 $\phi$ 与**接触 Jacobian $J=\nabla\phi$**。它们是接触约束的几何骨架，后面每个公式都用到。

gap 函数（有符号距离，signed distance）。 对两个物体（或物体与地面），gap 函数 $\phi(q)$ 定义为它们之间的**有符号最近距离**：

$\phi(q)>0$：分离（separated），无接触；
$\phi(q)=0$：恰好接触（touching）；
$\phi(q)<0$：穿透（penetration，物理上不允许，但软接触/离散误差下可能短暂出现）。

单边约束就是 $\phi(q)\ge 0$（不穿透）。注意 $\phi$ 是位形 $q$ 的**非线性**函数——这正是 §4 Stewart-Trinkle 需要线性化它的原因。对凸物体，$\phi$ 可由 GJK 算法（Gilbert-Johnson-Keerthi）计算；非凸需分解或采样。

接触 Jacobian。 gap 函数对位形的梯度：

\[ J=\frac{\partial\phi}{\partial q}=\nabla\phi(q)\in\mathbb{R}^{1\times n}. \]

它的物理意义至关重要：$J$ 把**广义速度** $v=\dot q$ 映射到**接触点的法向相对速度**——

\[ \dot\phi=\frac{\partial\phi}{\partial q}\dot q=J v. \]

即 $Jv$ 是"接触间隙的变化率"：$Jv>0$ 物体分离，$Jv<0$ 物体接近（侵入），$Jv=0$ 保持接触。这就是为什么速度级约束写成 $Jv^+\ge 0$（§2）——它要求接触间隙不减小（不侵入）。

更进一步，由虚功原理，接触力 $\lambda_n$（沿法向）对应的广义力是 $J^\top\lambda_n$——这就是为什么牛顿方程里接触力项是 $J^\top\lambda$（§0 开头那个方程）。$J$ 把接触空间（1 维，法向）与广义坐标空间（$n$ 维）双向连接：$Jv$ 是速度的"下投影"，$J^\top\lambda$ 是力的"上投影"。

多接触情形。 $k$ 个接触点时，每个有自己的 $\phi_i$、$J_i=\nabla\phi_i$。堆叠成 $J=[J_1;\dots;J_k]\in\mathbb{R}^{k\times n}$，则 $Jv\in\mathbb{R}^k$ 是所有接触的法向相对速度向量，$J^\top\lambda$ 是所有接触力的合广义力。Delassus 算子 $W=JM^{-1}J^\top\in\mathbb{R}^{k\times k}$ 就是把 $M^{-1}$ 夹在两个 $J$ 之间——它衡量"接触 $i$ 的单位冲量引起接触 $j$ 的速度变化"（§3 已用，这里给出它的来历）。

本质洞察：接触 Jacobian $J=\nabla\phi$ 是连接"接触世界"（1 维或 $k$ 维的间隙/法向力）与"机器人世界"（$n$ 维广义坐标）的**翻译官**。它有两个方向：$Jv$ 把速度翻译成"间隙变化率"（用于约束），$J^\top\lambda$ 把接触力翻译成"广义力"（用于动力学）。全章所有公式——MDI、LCP、凸 QP、MuJoCo——本质上都在这两个翻译之间运算。一旦理解 $J$ 的双向角色，你看任何接触公式都能立刻分清"哪部分在接触空间、哪部分在广义坐标空间"。

三重难度叠加¶

把上面的观察系统化，接触动力学的时间离散难在三个**叠加**的层次：

层次	困难	数学刻画	后果
(a) 解的正则性	速度有跳跃	$v(\cdot)$ 是 BV 函数（bounded variation），$q(\cdot)$ 仅 Lipschitz	经典积分器在跳跃点阶数塌陷
(b) 力的本质	冲击力是"无穷大作用无穷短时间"	作用力是测度（measure）而非函数，冲量含 atomic measure（原子测度 $\delta_{t_i}$）	不能写 $\lambda(t)$，只能写 $d\mathbf{R}$
(c) 约束集移动	可行域 $C(q)$ 随位形变	$C(q)=\{q:\phi(q)\ge 0\}$ 是移动的凸集（sweeping process）	约束的激活/失活是不连续的几何事件

逐条解释为什么每一条都"杀死"经典数值分析：

(a) BV 解。 经典 ODE 理论要求右端 $f$ Lipschitz，从而解 $C^1$。接触一来，$v$ 在冲击点不连续——它只是 BV 函数（变差有限，但可以跳）。Hairer-Lubich-Wanner 那本 644 页的《Geometric Numerical Integration》全程假设光滑 Hamiltonian 系统，它的辛积分器、对称积分器、变分积分器的所有误差估计都需要 $f\in C^k$。在 BV 解上，"$p$ 阶精度"无从谈起。

(b) 测度力。 持续接触时接触力是有限的普通函数 $\lambda_n(t)$；但冲击瞬间，为了在零时间内改变动量，力必须是 $\delta$ 测度（冲量 $P=\int\lambda\,dt$ 有限但 $\lambda\to\infty$）。所以正确的写法不是"$\lambda(t)$ 这个函数"，而是"$d\mathbf{R}$ 这个测度"，它分解为**勒贝格连续部分**（持续接触力 $\times\,dt$）加**原子部分**（冲击 $\sum_i P_i\delta_{t_i}$）。这就是为什么动力学要写成**测度微分包含**而不是微分方程。

(c) 移动约束。 可行集 $C(q)=\{q:\phi(q)\ge 0\}$ 的边界随时间扫过位形空间——Moreau 称之为**扫动过程（sweeping process）**。约束何时激活（$\phi$ 从正变零）、何时失活（从零变正）是不连续的几何事件，没有任何光滑结构可言。

本质洞察：这三重困难不是"工程上的麻烦"，而是**数学对象本身变了**。光滑 ODE 的解活在 $C^1$ 空间，接触动力学的解活在 BV 空间；光滑系统的"力"是函数，接触的"力"是测度。试图用 $C^1$ 空间的工具（RK、辛）去逼近 BV 空间的解，就像试图用有理数逼近 $\sqrt 2$ 却拒绝完备化——你需要换一套数学。这套数学就是 Moreau 的测度微分包含（§2）。

本专题的独立性：为什么 Stewart 要专门开一个分支¶

这里值得花几行说清本专题在数值分析版图中的**独立地位**，否则你会以为它只是"ODE 数值解的一个应用"。

光滑数值分析（RK、多步、辛、变分积分器）：处理 $\dot x=f(x),\ f\in C^k$。圣经是 Hairer-Lubich-Wanner 2006。
DAE 数值分析（等式约束）：处理 $\dot x=f(x,\lambda),\ g(x)=0$。有 index reduction、SHAKE/RATTLE 等成熟方法。
非光滑时步法（不等式约束）：处理 $M\,dv-F\,dt\in -d\mathbf{R},\ \phi(q)\ge 0$。这是 Stewart 在 SIAM Review 42:3 (2000) 专门开辟的分支，因为前两套工具在这里**全部失效**。

离散力学社区的变分积分器（Marsden-West 2001）**可以**扩展到接触（Fetecau-Marsden-West-Ortiz 2003 给出了保辛的碰撞离散，见 §8），但相比 Moreau-Jean 时步法在工程中使用较少——原因正是 §8 要讲的"辛性与穿透修正不相容"。

阶段小结：到这里我们确立了三件事——(1) 接触动力学的解是 BV 函数、力是测度、约束集移动，三重困难叠加；(2) 经典 ODE 积分器（含事件驱动法）在冲击点/Zeno 点失效；(3) 这是一个独立的数学分支，需要测度微分包含这套新语言。接下来 §1 先把"事件驱动 vs 时步法"两条路线摆清楚，再进入正确方程（§2）。

⚠️ 常见陷阱¶

💡 概念误区：认为"接触仿真 = 光滑 ODE + 碰撞检测打补丁" - 新手想法："正常用 RK4 积分，每步检查有没有穿透，穿透了就把速度翻转一下。" - 实际上：这是**事件驱动**思路的拙劣版本，在 Zeno 场景（bouncing ball 落定）彻底失效，且每次"翻转速度"破坏动量/能量一致性。正确的时步法把接触作为**每步求解的互补/优化问题的一部分**，与动力学**同时求解**，而非事后修补。 - 为什么重要：理解这一点，你才会明白为什么 MuJoCo/Drake 的核心是一个"求解器"（solver）而不是一个"碰撞检测器 + 积分器"。

🧠 思维陷阱：用"精度阶数"评判接触积分器 - 新手想法："这个积分器是 4 阶的，肯定比 1 阶的半隐式 Euler 准。" - 实际上：在**光滑段**阶数才有意义；在**冲击点**，由于解是 BV 的，所有积分器都退化到至多 $O(h)$（一阶），高阶方法没有任何优势，反而因为多点采样跨越跳跃而更糟。评判接触积分器的正确指标是**收敛到 MDI 解**（弱 * 收敛）、能量行为（是否人为注入/耗散能量）、互补满足程度（穿透量）。 - 正确思维：接触仿真选积分器问三个问题——(1) 每步解的是 LCP/NCP 还是凸 QP？(2) 能量是中性、耗散还是可能增长？(3) 是否保证非穿透或穿透有界？阶数排在最后。

💡 概念误区：把冲量和力混为一谈 - 新手想法："接触力就是 $\lambda_n$，冲击力就是很大的 $\lambda_n$。" - 实际上：持续接触力 $\lambda_n(t)$ 是**有限函数**，冲击对应的是**冲量** $P=\int\lambda_n\,dt$（atomic measure 的权重），瞬间内 $\lambda_n$ 趋于无穷而积分有限。把二者混用会导致你在离散格式里搞不清"未知量是力还是冲量"——Moreau-Jean 和 Stewart-Trinkle 的未知量都是**冲量** $P$（量纲 $N\cdot s$），不是力。这是后面所有公式的关键。 - 为什么重要：时步法的未知量统一是冲量，正因为它能同时表达"持续接触力 $\times h$"和"冲击 atom"——一套未知量统一两种物理。

练习¶

（推导题，草稿纸完成） 对一维 bouncing ball（$\dot v=-g$，地面 $\phi=q\ge 0$，恢复系数 $e$），从初始高度 $q_0$、零初速出发，推导第 $k$ 次落地的时刻 $t_k$ 的表达式，并证明 $\lim_{k\to\infty}t_k=T_\infty<\infty$。计算 $e=0.8,\ q_0=1\,\mathrm{m},\ g=9.81$ 时的 $T_\infty$。再说明：若用步长 $h=10^{-3}\,\mathrm{s}$ 的事件驱动法，它会在哪个时刻附近"卡死"？
（开放思考题） 三重困难 (a)(b)(c) 中，哪一条是"纯弹性碰撞 + 事件驱动"也无法回避的？哪一条在"持续接触无冲击"时消失？用一个具体例子（如方块静置桌面 vs 方块自由落体撞桌面）说明每条困难在什么物理情形下出现。
（反事实题） 假设我们生活在一个"接触力总是有限光滑函数"的世界（即把所有接触都换成很硬的弹簧 $\lambda_n=-k\phi$，$k$ 大但有限）。这时 (a)(b)(c) 三重困难分别变成什么？这个"软化"思路正是 MuJoCo 的选择（§6）——它用什么代价换取了 ODE 积分器的可用性？提示：考虑 $k\to\infty$ 时方程的刚性（stiffness）。

§1 两条方法论路线：事件驱动 vs 时步法 ⭐⭐¶

动机：面对同一个接触系统，有两种截然不同的算法哲学¶

§0 说清了接触动力学"难在哪"。现在面对一个具体的多体接触系统，工程史上演化出两条**根本不同**的求解哲学。它们不是"同一方法的两个版本"，而是对"接触是什么"的两种世界观。本节要把它们的分歧讲透，因为本章其余部分（§2–§9）几乎全部站在**时步法**这一侧——你必须先知道我们为什么放弃了另一侧。

事件驱动法（event-driven，又称 event-capturing 的对立面 event-tracking）：世界观是"接触是离散事件"。在两次事件之间，系统是光滑 ODE，用高阶积分器（RK45/DOPRI/LSODAR）积分；检测到事件（gap $\phi=0$）时，停下、用 root-finding 精确定位事件时刻、单独求解一个冲击问题、再重启光滑积分。
时步法（time-stepping，又称 event-capturing）：世界观是"接触是连续过程的一部分"。固定步长 $h$ 推进，不检测任何事件，每一步求解一个包含动力学+接触约束的互补问题（LCP/NCP）或优化问题（凸 QP/CCP）。冲击不是"事件"，而是这个互补问题解中冲量 $P$ 的 atomic 分量。

历史：两条路线的源头¶

事件驱动法源自天体力学和电路仿真传统——那里"事件"（行星交会、二极管开关）确实稀疏且可精确定位，LSODAR（含 root-finding 的 ODE 求解器）就是为此而生。把它用到机械接触上是 1980s–1990s 的自然延伸（Pfeiffer-Glocker 1996 的部分内容）。

时步法的现代形式由两个独立源头汇合：Moreau（1977, 1988） 从凸分析角度提出扫动过程的时间离散，Stewart-Trinkle（1996） 从工程角度提出 LCP 时步格式。Stewart（2000, SIAM Review 42:3）的综述把这条分支正式确立。其深层动机正是 §0 的 Zeno 与 Painlevé——事件驱动法在这两个现象上**数学上无定义**，而时步法给出唯一正确极限。

理论：逐维度对比¶

下表是本章最重要的"路线对照表"，请逐行理解而非死记：

维度	Event-driven（事件驱动）	Time-stepping（时步法）
光滑段积分	高阶变步长：RK45 / DOPRI / LSODAR	固定步长低阶：θ-method / 半隐式 Euler
事件检测	root-finding 求 $\phi(q(t))=0$ 的精确时刻	不检测，每步解 LCP/NCP/QP
冲击处理	事件点单独求解 impact LCP（独立子问题）	速度跳变被该步的互补解自动吸收
未知量	光滑段：$q,v$；事件点：冲量 $P$	每步：$v_{\ell+1}$ 与冲量 $P_{\ell+1}$ 同时
Zeno 行为	致命（bouncing ball 落定时无限次事件）	天然鲁棒（多次微弹被一步吸收）
Painlevé 悖论	无定义（接触力发散，无经典解）	给出含 atom 的测度解（impact without collision）
接触数规模	小（$\le 10$ 个接触）尚可	$10^6$ 量级 granular 也良好扩展
精度（光滑段）	高（可达 8 阶）	低（1–2 阶）
精度（冲击点）	受 Zeno/Painlevé 限制，常崩溃	鲁棒，弱 * 收敛到 MDI 解
典型实现	Siconos `EventDriven` + LSODAR	Moreau-Jean、Stewart-Trinkle、SAP、MuJoCo

几个关键对比的深入解读：

(1) "检测事件" vs "不检测事件" 是分水岭。 事件驱动法把计算预算花在"精确找到每个事件时刻"上——这在事件稀疏时高效（一次台球碰撞，精确到纳秒）。但接触数一多（足式机器人四足同时触地、granular 百万颗粒），事件几乎在每个瞬间都在发生，root-finding 的开销爆炸，且"哪个事件先发生"的排序本身就病态。时步法把预算花在"每步解一个互补问题"上——它不关心事件何时发生，只关心"这一步结束时状态是什么"。

(2) 冲击的地位天差地别。 这是两条路线最深刻的分歧。事件驱动法里，冲击是**外挂的额外步骤**："先积分光滑 ODE，碰到 $\phi=0$ 就暂停，套用 $v^+=-ev^-$，再继续"。时步法里，冲击是**互补问题解的一部分**——冲量 $P$ 作为未知量，当接触激活时它取正值（atom），否则为零，同一套方程统一处理持续接触与冲击。这正是 §2 Moreau 把 $v^++ev^-$ 直接塞进 MDI argument 的天才之处。

本质洞察：事件驱动法把接触看成"打断光滑世界的离散事件"，时步法把接触看成"连续测度过程的原子分量"。前者是**牛顿式**的（力是函数，冲击是奇异事件），后者是**Moreau 式**的（力是测度，冲击是测度的 atom）。机器人仿真全面倒向后者，不是因为它更"精确"（光滑段反而更糙），而是因为它在 Zeno/Painlevé/多接触这些"病理"情形下**仍有良定义的解**——而机器人恰恰天天遇到这些病理。

事件驱动法的另一面：Filippov 系统与滑模¶

为了公平，也为了理解事件驱动法**擅长**什么，我们补充一类它处理得比时步法更精细的现象——Filippov 系统的滑模（sliding mode）。这也是 §0 提到的"复合切换"的具体化。

什么是 Filippov 系统？ 一类**分段光滑**的动力系统：状态空间被一个**切换面**（switching surface）$\Sigma=\{x:s(x)=0\}$ 分成两个区域，每个区域内向量场光滑，但跨过 $\Sigma$ 时向量场**不连续**。例如干摩擦振子：物体向右运动时摩擦力向左，向左运动时摩擦力向右——速度过零（$\Sigma=\{v=0\}$）时摩擦力方向突变。

滑模现象。 当切换面两侧的向量场都"指向" $\Sigma$ 时，轨迹无法离开 $\Sigma$，被"困"在切换面上**滑动**——这叫滑模。物理例子：物体被压在传送带上，静摩擦让它粘住（stick），速度恒为零（在 $v=0$ 这个切换面上滑动）。Filippov 的理论给出滑模上的"等效动力学"：向量场取两侧的凸组合 $\dot x=\alpha f_+ + (1-\alpha)f_-$，$\alpha$ 由"保持在 $\Sigma$ 上"（$\dot s=0$）唯一确定。

事件驱动法在这里的优势。 滑模的进入/退出是精确的事件（$s=0$ 且向量场指向 $\Sigma$）。Piiroinen-Kuznetsov（2008, ACM TOMS）的事件驱动法能**精确检测**滑模的进入、在滑模上用等效动力学高精度积分、并**精确检测**退出——这比时步法把滑模"模糊"成一串小步要精细得多。对**继电反馈系统、干摩擦振子、电路**等滑模主导的系统，事件驱动法更准。

与接触的关系。 Coulomb 摩擦的 stick-slip（粘滞-滑动）切换正是一种滑模——粘滞相位（stick）是滑模（切向速度恒零），滑动相位（slip）是常规运动。这就是为什么**带摩擦的接触系统既是单边约束问题（法向）又是 Filippov 系统（切向 stick-slip）**——双重非光滑性叠加。时步法用互补/最大耗散统一处理这个 stick-slip，但不像事件驱动法那样精确分辨切换瞬间。

本质洞察：事件驱动法不是"过时的劣等方法"——它在**滑模主导、事件稀疏、要求高精度切换**的系统（Filippov 系统、电路、单次精密碰撞）里**优于**时步法。两条路线的取舍本质是"事件密度"：事件稀疏且需精确 → 事件驱动；事件密集或病态（Zeno/Painlevé/多接触）→ 时步法。机器人接触恰好落在后者，所以全选时步法；但若你做的是精密碰撞标定或干摩擦振子分析，别忘了事件驱动法这把更锋利的小刀。

桥接：机器人工程的现实选择¶

理论讲完，回到工程。为什么所有主流机器人仿真器都是 time-stepping？

任务类型	接触特征	路线选择	理由
单次高精度冲击（台球、撞针、碰撞标定）	接触稀疏、要求高精度冲量	event-driven 更准	光滑段高阶 + 精确事件定位
足式运动（locomotion）	多足周期性触地/离地、$e\approx 0$	time-stepping	接触频繁、需鲁棒、Zeno 风险
灵巧操作（dexterous manipulation）	多指多面接触、滑动/粘滞切换	time-stepping	接触数多、摩擦主导
Granular（沙、颗粒）	$10^5$–$10^6$ 接触	time-stepping	事件驱动完全不可行
RL 训练 / 可微仿真	需大批量并行 + 梯度	time-stepping（软化/凸）	凸性 ⇒ 可微（见 §5–§6）

结论：ROS / MuJoCo / Drake / Bullet / Pinocchio / Isaac 全部是 time-stepping；只有 Siconos 两者并存（它是数学研究导向的库）。这就是为什么本章把全部篇幅给了 time-stepping。

阶段小结：到这里我们确立了两条路线的分野，并解释了机器人为何全选时步法。接下来必须回答一个更根本的问题：时步法的离散格式，应该收敛到什么"正确答案"？ 这个"正确答案"就是 §2 的测度微分包含——它是连续时间的"真方程"，§3–§6 的所有离散格式都以收敛到它为合法性标准。

⚠️ 常见陷阱¶

💡 概念误区：把时步法当成事件驱动法的"近似"或"偷懒版" - 新手想法："事件驱动法精确处理每次碰撞，时步法只是固定步长懒得检测事件，所以是近似。" - 实际上：时步法是**独立数学对象**（Moreau 扫动过程的原生离散），Stewart 1998 证明其 $h\to 0$ 极限**就是** MDI 解。在 Zeno 场景下事件驱动法**根本无定义**（事件无穷多、停不下来），而时步法给出**唯一正确**极限。所以不是"时步法近似事件驱动法"，恰恰相反——在病理情形下时步法是唯一有定义的那个。 - 为什么重要：理解这点，你才不会在 bouncing ball 落定时去给时步法"加事件检测补丁"——那是画蛇添足且破坏鲁棒性。

🧠 思维陷阱：认为"高阶积分器 = 更好的仿真器" - 新手想法："Siconos 支持 RK45 事件驱动，MuJoCo 只有半隐式 Euler，所以 Siconos 更准。" - 实际上：在接触主导的机器人任务里，步长稳定性、可微性、并行能力、warm-start 远比光滑段阶数重要。MuJoCo 用一阶半隐式 Euler 却能百万环境并行 + 可微，正是工程上的正确取舍。光滑段的 8 阶精度对足式/操作任务毫无意义——因为运动几乎处处有接触切换。 - 正确思维：选仿真器问"接触刚度、接触数、可微性需求、实时约束"四维度，而非"积分器阶数"。

💡 概念误区：以为"不检测事件"意味着"忽略冲击" - 新手想法："时步法不检测碰撞事件，那它怎么处理弹跳？岂不是穿过去了？" - 实际上：时步法**不显式检测**事件，但通过每步的互补问题**隐式捕获**冲击——冲量 $P$ 作为未知量在接触激活时自动取非零 atom 值。"不检测"指的是不用 root-finding 找事件时刻，不是**不处理冲击。这正是 event-**capturing（捕获）一词的由来。 - 为什么重要：这是 §2–§4 全部公式的核心——未知量是冲量 $P$，互补条件 $0\le U\perp P\ge 0$ 自动决定每步是否有冲击。

练习¶

（开放思考题） 列举三个机器人任务，分别论证应该选事件驱动法还是时步法，并说清你的判据是接触数、Zeno 风险、精度需求还是可微性。其中至少一个任务应当是"两种方法都可接受"的边界情形，说明权衡点。
（反事实题） 假设有人坚持用事件驱动法仿真四足机器人小跑（trot）。请具体描述会在哪个物理瞬间出问题：当四条腿在着地相位频繁切换、且落地时 $e\approx 0$（塑性）时，root-finding 会面临什么？提示：考虑"几乎同时但不完全同时"的多接触事件如何排序，以及 $e\to 0$ 时的退化。
（综合判断题） 给定一个仿真器（你熟悉的任意一个），查阅其文档/源码，回答：它是 event-driven 还是 time-stepping？固定步长还是变步长？每步求解的是 LCP、NCP 还是凸 QP？把答案填入本章 §6 的"仿真器对比表"对应行，验证你的判断与表中一致。

§2 Moreau 扫动过程与测度微分包含（MDI）⭐⭐⭐¶

动机：我们需要一个能容纳"跳跃速度 + 测度力"的方程¶

§0 论证了接触动力学的解是 BV 函数、力是测度。但我们写下的牛顿方程 $M\dot v=F+J^\top\lambda$ 假设 $\dot v$ 存在（$v$ 可微）——这在冲击点不成立。我们需要一个新的方程形式，它能在速度跳跃、力是测度的情形下仍然成立。 这就是 Moreau 在 1970s–1980s 发展的**测度微分包含（measure differential inclusion, MDI），以及它的几何图景**扫动过程（sweeping process）。本节是全章理论的**地基**——§3–§6 的每一种离散格式，都是这个连续时间方程的一种离散，合法性由"$h\to 0$ 是否收敛到它"判定。

反面：为什么不能用"微分方程 + Dirac 函数"硬凑¶

一个自然的补救想法是：既然冲击力是 $\delta$ 函数，那就写 $M\dot v=F+\sum_i P_i\delta(t-t_i)+J^\top\lambda$，把冲击当 Dirac 项加进去不就行了？

这个想法在**单次、已知时刻**的冲击下勉强可用，但有三个致命问题：

冲击时刻 $t_i$ 是未知的——它由动力学自己决定（$\phi=0$ 何时发生），你不能预先把 $\delta(t-t_i)$ 写进方程。
冲击大小 $P_i$ 也是未知的——它由冲击律（恢复系数）和约束共同决定。
持续接触与冲击需要分别处理——持续接触是勒贝格连续力，冲击是 atom，你得手动判断当前是哪种情形并切换方程。这正是事件驱动法的困境。

Moreau 的洞见是：不要把力写成"函数 + Dirac"，而是直接写成一个测度 $d\mathbf{R}$，并用法锥/切锥的包含关系约束它。 这样冲击时刻、大小、持续/瞬时之分，全部由包含关系自动决定，无需人工切换。

历史：Moreau 的扫动过程¶

Jean-Jacques Moreau（1923–2014），法国力学家与凸分析先驱（"Moreau 包络"、"Moreau 分解"都以他命名），在 1971–1988 年间发展了**扫动过程（sweeping process）** 理论。最初的动机是塑性力学（弹塑性材料的应力被"扫动"在屈服面内），后来 Moreau（1988, CISM Courses and Lectures 302）将其用于单边约束多体动力学。Michelle Schatzman 和 Michel Jean 在 1980s–1990s 完善了离散格式（Moreau-Jean scheme，见 §3）。

扫动过程的几何图景极美：想象一个点 $q(t)$ 被一个**移动的凸集** $C(t)$ "扫动"——当点在集合内部时自由运动，当点撞到集合边界时，边界"推"着它（沿内法向）不让它出界。这个"推力"就是法锥方向的接触力。

理论：从单边约束到测度微分包含¶

Step 1：可行集与约束。 位形 $q(t)$ 必须满足单边约束 $\phi(q)\ge 0$（gap 非负，不穿透）。定义**可行集**

\[ C(q)=\{q':\phi(q')\ge 0\}. \]

注意 $C$ 依赖于约束几何，随接触对的相对位形变化——这就是"移动的凸集"。

Step 2：速度级运动学。 在接触边界 $\phi(q)=0$ 上，物理上允许的速度受限：物体不能继续侵入，即 $\nabla\phi^\top v\ge 0$（分离或相切，不准穿入）。这定义了**切锥**

\[ T_C(q)=\{v:\nabla\phi(q)^\top v\ge 0\ \text{ for all active }\phi\}. \]

未接触时（$\phi>0$）切锥是全空间（任意速度合法）；接触时切锥是半空间（只准分离方向）。

Step 3：法锥与接触力方向。 接触力 $\lambda_n\ge 0$ 沿内法向 $\nabla\phi$ 作用，且只在接触时（$\phi=0$）非零。用凸分析语言，接触力（作为广义力）属于切锥在 $v^+$ 处的**法锥**

\[ -J^\top\lambda \in N_{T_C(q)}(v^+), \]

其中 $N_K(x)=\{g:\langle g,y-x\rangle\le 0,\ \forall y\in K\}$ 是凸锥 $K$ 在 $x$ 的法锥。这个包含关系**同时编码**了三件事：(i) 接触力非负（$\lambda_n\ge 0$）；(ii) 只在接触时有力（互补）；(iii) 力沿法向。

Step 4：动力学写成测度微分包含。 把牛顿方程的"$M\dot v\,dt$"替换为测度"$M\,dv$"（$dv$ 是速度的微分测度，含跳跃），把接触力替换为接触冲量测度 $d\mathbf{R}$，得到 Moreau 测度微分包含：

\[ \boxed{\,-M\,dv + F(q,v)\,dt\ \in\ N_{T_C(q)}\bigl(v^+(t)\bigr)\,d\mu,\qquad d\mu=dt+\sum_i\delta_{t_i},\ \ v^+(t)=\lim_{s\downarrow t}v(s).\,} \]

逐项解读这个核心方程：

$M\,dv$：动量的微分测度。$v$ 是 BV 函数，$dv$ 分解为连续部分（$\dot v\,dt$，光滑段加速度）加跳跃部分（$\sum_i(v_i^+-v_i^-)\delta_{t_i}$，冲击时的速度跳变）。这是关键——$dv$ 允许跳跃，而 $\dot v\,dt$ 不允许。
$F(q,v)\,dt$：外力（重力等），它是普通函数乘勒贝格测度 $dt$，没有 atom（外力不会在瞬间无穷大）。
$d\mu=dt+\sum_i\delta_{t_i}$：基测度，勒贝格测度加冲击时刻的原子测度。这一项让方程**同时**容纳持续接触（在 $dt$ 上）与冲击（在 $\delta_{t_i}$ 上）。
$N_{T_C(q)}(v^+)$：在**右极限速度** $v^+$ 处的切锥法锥。用 $v^+$ 而非 $v^-$ 或 $v$，是因为接触力作用后速度变为 $v^+$，物理上分离方向由 $v^+$ 决定。

本质洞察：MDI 的天才在于把"力"从函数升级为**测度**，把"导数"从 $\dot v$ 升级为**微分测度 $dv$。这一升级让冲击不再是"打断方程的奇异事件"，而是"$dv$ 测度的 atomic 分量"——持续接触是 $dv$ 的连续部分，冲击是 $dv$ 的跳跃部分，**同一个方程统一两者。这正是 §0 那三重困难（BV 解、测度力、移动约束）的统一解决：BV 解 ↔ $dv$ 含跳跃；测度力 ↔ $d\mathbf{R}$；移动约束 ↔ $N_{T_C(q)}$ 随 $q$ 变。

Moreau 的关键代数化：从 $N_C$ 到 $N_{T_C}$¶

这里有一个微妙但极重要的技术点，是 Moreau 框架"对接触激活鲁棒"的根源。

朴素想法是把接触力写成位置级法锥 $-J^\top\lambda\in N_{C(q)}(q)$（位形 $q$ 在可行集 $C$ 的法锥里）。但 $N_{C(q)}(q)$ 在 $q$ 严格内部（$\phi>0$）是 $\{0\}$，在边界（$\phi=0$）突然变成一条射线——这个**不连续切换**正是接触激活的不光滑性，难以离散。

Moreau 的代数化：把位置级法锥 $N_{C(q)}(q)$ 替换为**速度级**的切锥法锥 $N_{T_C(q)}(v^+)$。区别在于：

位置级 $N_{C(q)}(q)$：问"位形 $q$ 在可行集边界吗"，答案是 0/1 突变。
速度级 $N_{T_C(q)}(v^+)$：问"在当前接触几何下，速度 $v^+$ 是否被切锥约束"，这是对**速度**的约束，离散时只需在每步检查 $\phi=0$ 的接触上施加 $\nabla\phi^\top v^+\ge 0$。

这个替换使得**时步公式对接触激活状态的变化天然鲁棒**——你不需要精确知道"何时 $\phi$ 从正变零"，只需在每步对当前 $\phi\le 0$（或线性化后将要 $\le 0$）的接触施加速度级约束。这是时步法"不检测事件"的数学根据。

速度级 vs 位置级约束：一个贯穿全章的分歧¶

上面的代数化引出一个贯穿 §3–§4 的核心选择：约束写在**位置级**还是**速度级**？

	位置级（position-level）	速度级（velocity-level）
约束	$\phi(q)\ge 0$（gap 非负）	$\phi=0$ 时 $\nabla\phi^\top v^+\ge 0$（不准侵入速度）
保证	严格不穿透	速度合法，但位置可能漂移
代表格式	Stewart-Trinkle（§4）	Moreau-Jean（§3）
代价	需非线性 gap 函数展开	需位置投影防漂移
冲击律嵌入	不直接	$v^++ev^-$ 直接进 argument

二者在持续接触（$\phi\equiv 0$）下等价，但在"即将接触但尚未接触"时行为不同。理解这个分歧，是看懂 §3 与 §4 区别的钥匙。

本质洞察：位置级约束像"硬墙"——严格不让穿透，但代价是要处理非线性 gap，且在激活瞬间不光滑。速度级约束像"软规则"——只管速度方向合法，位置允许小幅漂移（靠投影修正）。Moreau 选速度级，是为了让冲击律 $v^++ev^-$ 能直接进 MDI 的 argument，从而"持续接触 + 冲击"一套公式（§3 会看到这有多优雅）。

摩擦的引入：从法锥到 de Saxcé 双位锥¶

上面只讲了法向接触（无摩擦）。引入 Coulomb 摩擦后，接触力不再只沿法向，而是被约束在**摩擦锥** $K=\{(\lambda_n,\beta):\|\beta\|\le\mu\lambda_n\}$ 内。最大耗散原理（本目录 20）要求摩擦力与切向滑动速度反平行。

Moreau 框架的处理：把无摩擦时的法锥 $N_{T_C}$ 替换为 de Saxcé 双位锥（bipotential） 关系——这是一个把法向互补与切向最大耗散统一起来的凸分析结构。其核心是 de Saxcé 1992 的双位势 $b(v,r)=\mu\|v_t\|\,r_n+\langle v,r\rangle$，使得摩擦接触可写成一个隐式标准化（implicit normalization）的包含。细节属于本目录 20 的范畴，这里只需记住结论：

\[ -d\mathbf{R}\in\partial_{\text{de Saxcé}}\ \text{(广义法锥)}\quad\Longleftrightarrow\quad \text{法向互补 + 切向最大耗散 + 摩擦锥约束。} \]

这个统一化是后面 §3（Dzonou-Monteiro Marques 收敛证明克服 Coulomb 律非半连续性）的关键。

一个几何算例：扫动过程如何"扫"出接触力¶

为了让"移动凸集扫动一个点"从比喻落到具体，我们算一个最小例子，把法锥、切锥、速度级约束全部可视化。

设置：一维质点 $q\in\mathbb{R}$，可行集 $C=\{q:q\ge 0\}$（地面在原点，物体在右侧合法）。质点以速度 $v$ 运动，无外力。

情形一：物体在内部（$q>0$）。 此时 $\phi=q>0$，约束不激活。切锥 $T_C(q)=\mathbb{R}$（全空间，任意速度合法），法锥 $N_{T_C}(v)=\{0\}$（无约束力）。MDI 退化为 $M\,dv=0$，即 $v$ 不变——自由运动。几何图景：点在集合内部，边界"够不着"它，无推力。

情形二：物体在边界且向外（$q=0,\ v<0$，要穿出）。 这是关键情形。切锥 $T_C(0)=\{u:u\ge 0\}$（只准向右/相切，不准向左穿出地面）。当前速度 $v<0$ 不在切锥内（违反约束）。

MDI 要求 $-M\,dv\in N_{T_C(0)}(v^+)$，即接触冲量把 $v$ "拉回"切锥。法锥 $N_{T_C(0)}(v^+)$ 在 $v^+=0$（被拉到边界）时是 $\{g:g\le 0\}$（指向集合内部的法向）。求解：冲量 $P$ 使 $v^+=0$（停在边界，塑性 $e=0$）或 $v^+=-ev^-=e|v|>0$（回弹）。

几何图景：边界像一堵墙，当点试图穿出（$v<0$）时，墙沿内法向（$+$ 方向）"推"点，推力大小恰好让点不穿出（$v^+\ge 0$）。 这个"推力"就是法锥方向的接触冲量——它不是预先设定的，而是"恰好够用"地由包含关系决定。

情形三：物体在边界但向内（$q=0,\ v>0$，要离开）。 速度 $v>0$ 在切锥内（合法，正在分离）。无需约束力，$P=0$，物体自由离开。几何图景：点正在离开边界，墙不挽留。

本质洞察：这个算例揭示了扫动过程的核心机制——接触力是"约束的影子"。它只在点试图违反约束（情形二）时出现，大小"恰好够用"地把速度拉回切锥，方向沿内法向。这与"预先设定一个弹簧力"（软接触，§6）截然不同：扫动过程的力是**互补涌现**的（要么 $\phi>0$ 无力，要么 $\phi=0$ 有恰好够用的力），而非函数式给定的。三种情形对应互补条件 $0\le\phi\perp\lambda\ge 0$ 的三个分支：内部（$\phi>0,\lambda=0$）、压紧（$\phi=0,\lambda>0$）、分离（$\phi=0,\lambda=0$ 临界）。把这个一维图景推广到多维多接触，就是完整的 MDI。

对比：扫动过程 vs 罚函数（软接触）的两种世界观¶

这里值得插入一个贯穿全章的对比，它预告了 §6 MuJoCo 的不同选择：

	扫动过程（硬接触，§2–§5）	罚函数 / 软接触（§6 MuJoCo）
接触力来源	互补涌现（恰好够用）	函数式给定（弹簧 $\lambda=-k\phi$）
穿透	严格零（或漂移修正）	有意允许（$\phi<0$ 产生力）
数学对象	法锥包含 / LCP	光滑函数
可微性	激活集切换处不可微	处处可微
刚度	无穷硬	有限 $k$（可调）

这不是"哪个对"的问题，而是两种世界观：硬接触追求物理精确（零穿透），软接触追求数值友好（可微、稳定）。本章 §2–§5 走硬接触，§6 走软接触——它们在 $k\to\infty$ 时统一（软接触的刚性极限就是硬接触）。

阶段小结：到这里我们建立了连续时间的"正确方程"——Moreau MDI $-M\,dv+F\,dt\in N_{T_C(q)}(v^+)\,d\mu$，理解了它如何用测度+法锥统一持续接触与冲击、用切锥法锥实现对接触激活的鲁棒、用速度级 vs 位置级的选择埋下 §3/§4 的分歧、用 de Saxcé 双位锥引入摩擦，并通过一维几何算例看清"接触力是约束的影子"。接下来 §3 把这个连续方程**离散**成 Moreau-Jean θ-scheme——第一种、也是数学最严谨的时步格式。

⚠️ 常见陷阱¶

💡 概念误区：把 MDI 当成"加了 Dirac 项的微分方程" - 新手想法："$M\,dv$ 不就是 $M\dot v\,dt$ 嘛，MDI 就是普通牛顿方程加几个冲击 Dirac。" - 实际上：$dv$ 是**微分测度**，本质上是 $v$（BV 函数）的 Stieltjes 测度，它**内禀地**包含跳跃，不是"$\dot v\,dt$ 外加 Dirac"。区别在于：在 MDI 里，跳跃时刻和大小由包含关系**自动决定**；在"$\dot v\,dt$ + Dirac"里你必须**手动指定** Dirac 在哪、多大。前者是自洽的数学对象，后者是 ad hoc 拼凑。 - 为什么重要：只有理解 $dv$ 是测度，你才懂为什么 §3 的离散里未知量是冲量 $P$（$=\int d\mathbf{R}$）而非力。

🧠 思维陷阱：以为法锥写在 $v^-$ 或 $v$ 上无所谓 - 新手想法："$N_{T_C}(v^+)$ 里用 $v^+$、$v^-$ 还是平均 $v$ 不都差不多？" - 实际上：必须用 $v^+$（右极限速度）。物理上，接触力作用之后速度变为 $v^+$，分离/侵入由 $v^+$ 判定。若用 $v^-$，则冲击前的侵入速度会被错误地"允许"，无法阻止穿透。Moreau 框架的因果性（力的作用导致 $v^-\to v^+$）要求 argument 是 $v^+$。这也是 §3 离散里 $v_{\ell+1}$（新速度）出现在互补条件里、而非 $v_\ell$ 的原因。 - 正确思维：接触约束是对"作用后速度"的约束——永远盯着 $v^+$。

💡 概念误区：认为扫动过程只是个比喻 - 新手想法："'移动凸集扫动一个点'听起来像直觉比喻，不是严格数学。" - 实际上：扫动过程是**完全严格**的数学对象，有完整的存在唯一性理论（Moreau 1977 对凸 $C(t)$、后续对 prox-regular 集的推广）。"扫动"二字对应的是 $C(t)$ 边界推着 $q(t)$ 的法锥包含 $-\dot q\in N_{C(t)}(q)$（一阶扫动过程）或其二阶动力学版本。它不是比喻，是定理。 - 为什么重要：正因为它是严格的，§3–§4 的离散格式才能证明收敛到它（Stewart 1998、Moreau-Jean 弱收敛）。

练习¶

（推导题，草稿纸完成） 对一维质点撞墙（$\phi=q\ge 0$，$M=m$，无外力），写出 Moreau MDI。在冲击时刻 $t_*$，把 MDI 在该原子测度上积分，推导出冲量 $P$ 与速度跳变 $v^+-v^-$ 的关系。再加入 Newton 冲击律 $v^+=-ev^-$，解出 $P$。验证 $e=1$ 时 $P=-2mv^-$（完全弹性，动量翻转），$e=0$ 时 $P=-mv^-$（完全塑性，速度归零）。
（开放思考题） 为什么 Moreau 选择把法锥写在**切锥** $T_C$ 上（即 $N_{T_C(q)}(v^+)$）而不是直接写在可行集 $C$ 上（即 $N_{C(q)}(q)$）？从"接触激活的不连续性"角度论证速度级表述的离散优势。提示：考虑一个物体从 $\phi>0$ 逐渐接近 $\phi=0$ 的过程，两种表述各在哪个瞬间"开关"接触力。
（综合题，需结合本目录 10/20） 把 Moreau MDI 与本目录 10 的 LCP 联系起来：在持续接触（无冲击、$d\mu=dt$）的光滑段，证明 MDI $-M\dot v+F\in N_{T_C}(v^+)$ 等价于一个法向接触力的互补条件 $0\le\lambda_n\perp(\nabla\phi^\top v^+\ \text{或}\ \phi)\ge 0$。再加入摩擦锥（本目录 20），说明为什么这变成 NCP（非线性互补）而非 LCP。

§3 Moreau-Jean θ-scheme：速度级离散 ⭐⭐⭐¶

动机：把连续 MDI 变成计算机能跑的迭代¶

§2 给了连续时间的"正确方程" MDI。但计算机不能积分测度——我们需要把它**离散**成"已知 $(q_\ell,v_\ell)$，求 $(q_{\ell+1},v_{\ell+1})$"的迭代格式。Moreau-Jean θ-scheme 是第一种、也是数学最严谨的离散，它直接对应 Siconos 的 MoreauJeanOSI，是 Pinocchio 3 ADMM 的理论祖先。本节要**完整推导**这个格式，并给出它最漂亮的性质——能量估计。

反面：朴素离散会出什么问题¶

先看一个朴素尝试，理解 Moreau-Jean 的设计为何如此。最直接的想法是显式 Euler：

\[ M(v_{\ell+1}-v_\ell)=h F_\ell + H^\top P_{\ell+1},\qquad q_{\ell+1}=q_\ell+h v_\ell, \]

其中 $H=\nabla\phi$ 是接触 Jacobian，$P_{\ell+1}$ 是这一步的接触冲量。问题：

接触冲量 $P_{\ell+1}$ 满足什么条件？ 若写成位置级 $\phi(q_{\ell+1})\ge 0$，但 $q_{\ell+1}=q_\ell+h v_\ell$ 用的是**旧速度** $v_\ell$，接触力来不及作用，物体已经穿透了。
显式格式刚性差：硬接触对应大刚度，显式 Euler 需要极小步长才稳定。

Moreau-Jean 的修正：(i) 用**速度级**约束（写在 $v_{\ell+1}$ 上，新速度）；(ii) 位置更新用 $\theta$ 加权的速度（半隐式），获得可调的稳定性与能量行为。

历史：Moreau 与 Jean¶

如 §2 所述，Moreau（1988）提出连续扫动过程，Michelle Schatzman 与 Michel Jean 完善离散格式。"θ-method"（又称 Crank-Nicolson 型加权）来自传统 ODE 数值分析：$\theta\in[0,1]$ 在显式（$\theta=0$）与隐式（$\theta=1$）间插值，$\theta=\tfrac12$ 是二阶的梯形法。Moreau（CMAME 177, 1999）与 Jean（CMAME 177, 1999）分别给出了收敛性分析的奠基论文。

理论：θ-scheme 的完整推导¶

Step 1：动力学的 θ-加权离散。 把 MDI 在一个时间步 $[t_\ell,t_{\ell+1}]$（步长 $h$）上积分。连续部分外力用 θ-加权（$\theta\in[0,1]$），冲量 $P_{\ell+1}=\int_{[t_\ell,t_{\ell+1}]}d\mathbf{R}$ 整体作为未知量：

\[ \boxed{\,M(v_{\ell+1}-v_\ell)=h\bigl[\theta F_{\ell+1}+(1-\theta)F_\ell\bigr]+H_{\ell+\theta}^\top P_{\ell+1}.\,} \]

这里 $F_{\ell+1}=F(q_{\ell+1},v_{\ell+1})$ 用新状态（隐式），$H_{\ell+\theta}=\nabla\phi(q_{\ell+\theta})$ 在中间位形 $q_{\ell+\theta}=q_\ell+\theta h v_\ell$（或更精细的迭代）处评估。

Step 2：位置的 θ-加权更新。

\[ \boxed{\,q_{\ell+1}=q_\ell+h\bigl[\theta v_{\ell+1}+(1-\theta)v_\ell\bigr].\,} \]

$\theta=\tfrac12$ 时这是梯形法（位置用前后速度平均），$\theta=1$ 时用新速度（后向 Euler），$\theta=0$ 时用旧速度（前向 Euler）。

Step 3：接触的速度级互补条件（含冲击律）。 这是 Moreau-Jean 最精妙的一步。接触冲量 $P_{\ell+1}$ 不属于普通互补，而是属于**切锥法锥**，且**冲击律直接嵌入 argument**：

\[ \boxed{\,P_{\ell+1}\in -N_{T_C(q_{\ell+\theta})}\bigl(v_{\ell+1}+e\,v_\ell\bigr).\,} \]

注意 argument 是 $v_{\ell+1}+e\,v_\ell$——这正是 Newton 冲击律 $v^+=-ev^-$ 的离散嵌入！展开看法向分量：定义"含恢复的相对法向速度" $U_{\ell+1}=H(v_{\ell+1}+e v_\ell)$，则法向互补条件为

\[ 0\le U_{\ell+1}^n\ \perp\ P_{\ell+1}^n\ge 0. \]

Step 4：消元得到每步的 LCP。 把动力学（Step 1）代入互补条件（Step 3）。先解出 $v_{\ell+1}$（设 $F$ 暂时显式或线性化）：

\[ v_{\ell+1}=\underbrace{v_\ell+M^{-1}h\bar F}_{v_{\text{free}}\ (\text{自由速度})}+M^{-1}H^\top P_{\ell+1}, \]

其中 $v_{\text{free}}$ 是"无接触时这一步的速度"（自由飞行）。代入 $U_{\ell+1}=H(v_{\ell+1}+e v_\ell)$：

\[ U_{\ell+1}=\underbrace{HM^{-1}H^\top}_{W\ (\text{Delassus 算子})}P_{\ell+1}+\underbrace{H(v_{\text{free}}+e v_\ell)}_{U_{\text{free}}}. \]

于是法向接触归结为一个 LCP：

\[ \boxed{\,0\le U_{\ell+1}=W\,P_{\ell+1}+U_{\text{free}}\ \perp\ P_{\ell+1}\ge 0,\qquad W=HM^{-1}H^\top.\,} \]

$W=HM^{-1}H^\top$ 称为 Delassus 算子（Delassus operator），它是接触空间中的"逆惯量"——衡量"单位接触冲量引起多大的接触点相对速度变化"。

Step 5：$W$ 的半正定性 ⇒ LCP 可解。 由于 $M\succ 0$（质量矩阵正定），$M^{-1}\succ 0$，故对任意 $x$：

\[ x^\top W x=x^\top H M^{-1}H^\top x=(H^\top x)^\top M^{-1}(H^\top x)\ge 0, \]

即 $W\succeq 0$（半正定）。由本目录 10 的 LCP 理论，半正定矩阵的 LCP 若可行则有解（且可用 Lemke 算法或 PGS 求解）。这保证了 Moreau-Jean 每步法向 LCP 可解。（含摩擦时 $W$ 块结构更复杂，可解性需 §4 的 copositive-plus 论证。）

几何直觉：Delassus 算子是什么¶

$W=HM^{-1}H^\top$ 这个对象值得停下来理解，它贯穿全章。

物理意义：给接触点施加单位冲量 $P$，接触点的相对速度变化 $\Delta U=WP$。所以 $W$ 是"接触点的有效柔度（compliance）"。
退化情形：单接触单自由度时 $W=1/m_{\text{eff}}$（有效质量的倒数），冲量 $P$ 改变速度 $P/m_{\text{eff}}$——回到中学物理。
多接触耦合：$W$ 的非对角元 $W_{ij}=H_i M^{-1}H_j^\top$ 表示"接触 $i$ 的冲量如何影响接触 $j$ 的速度"——这是多体接触相互耦合的来源（一个接触发力会通过刚体传到另一个接触）。
病态来源：当某些接触方向"几乎平行"或质量分布极端时，$W$ 接近奇异（condition number 大），导致求解器收敛慢——这正是 Pinocchio 3 用 ADMM + proximal 正则化来对付的 "ill-conditioned Delassus" 问题。

本质洞察：整个时步法的代数核心，就是反复求解"以 Delassus 算子 $W$ 为系统矩阵的 LCP/凸 QP"。$W=HM^{-1}H^\top$ 把刚体惯量 $M^{-1}$ 投影到接触空间——它就是接触问题的"系统矩阵"。后面 §4 的 Stewart-Trinkle LCP 矩阵、§5 的凸 QP Hessian、§6 的 $A=JM^{-1}J^\top$，全都是这个 $W$ 的变体。认出这一点，四种格式就统一了。

能量估计：θ 的物理含义¶

Moreau-Jean 最漂亮的性质是其**能量行为**可以解析分析。这解释了为什么 $\theta=\tfrac12$ "能量中性"、$\theta=1$ "强耗散"。

光滑段（无接触）能量变化。 动能 $E=\tfrac12 v^\top M v$，一步变化

\[ \Delta E=\tfrac12\bigl(v_{\ell+1}^\top M v_{\ell+1}-v_\ell^\top M v_\ell\bigr)=\tfrac12(v_{\ell+1}+v_\ell)^\top M(v_{\ell+1}-v_\ell). \]

代入动力学 $M(v_{\ell+1}-v_\ell)=hF_{\ell+\theta}$（无接触，$\bar F=\theta F_{\ell+1}+(1-\theta)F_\ell$）。用恒等式 $\tfrac12(v_{\ell+1}+v_\ell)=v_{\ell+\theta}-(\theta-\tfrac12)(v_{\ell+1}-v_\ell)$：

\[ \Delta E=h\Bigl(v_{\ell+\theta}+(\theta-\tfrac12)(v_{\ell+1}-v_\ell)\Bigr)^\top F_{\ell+\theta}. \]

取 $\theta=\tfrac12$：括号内第二项消失，$\Delta E=h\,v_{\ell+1/2}^\top F$——恰好是力做的功（中点法则），能量中性（无数值耗散也无注入）。这是辛/对称积分器的离散能量守恒类比。
取 $\theta=1$：第二项 $=\tfrac12(v_{\ell+1}-v_\ell)=\tfrac12 M^{-1}hF$，代入得多出一项 $-\tfrac{h^2}{2}\|F\|_{M^{-1}}^2<0$，数值耗散 $O(h^2)$（后向 Euler 的人工阻尼）。

接触段能量变化。 冲量做功 $\Delta E_{\text{contact}}=-\tfrac12 P^\top(U^++U^-)$（动能定理的接触版本，符号取决于约定）。在 Moreau-Jean 框架下可证：$P^\top(v^++ev^-)/2\ge 0$ 保证接触**额外耗散**（$e<1$ 时塑性碰撞耗能，$e=1$ 时保能）。

把光滑段与接触段合起来：

θ	光滑段	接触段	总体行为	适用
$\tfrac12$	能量中性	接触耗散（$e<1$）	二阶精度、能量中性	弹性碰撞、长时间守恒轨迹
$1$	数值耗散 $O(h^2)$	接触耗散	一阶精度、强耗散	刚性接触、需阻尼抑制振荡
$(\tfrac12,1)$	中等耗散	接触耗散	介于两者	工程折中

本质洞察：θ 不是随便选的数值参数，它直接控制**能量耗散**。$\theta=\tfrac12$ 像辛积分器（能量长期不漂移，适合天体/长轨迹），$\theta=1$ 像隐式 Euler（人为耗散，把高频振荡压平，适合硬接触防 jitter）。没有"最好"的 θ——它是精度与稳定性的旋钮：要长期能量守恒选 $\tfrac12$，要鲁棒压振荡选接近 $1$。这与 §8 的"半隐式 vs 全隐式"权衡是同一回事的两种语言。

把每步 LCP 看成一个投影：通往 Moreau-Yosida 与 ADMM¶

最后一个视角，它把 Moreau-Jean 与现代 ADMM/proximal 求解器（Pinocchio 3）连起来，也呼应 §2 用 Moreau 包络引入摩擦的伏笔。

回到 §3 的每步法向 LCP $0\le U=WP+U_{\text{free}}\perp P\ge 0$。这个 LCP 等价于一个**投影**：解 $P$ 是把自由速度对应的"无约束解"投影到可行锥上。更精确地，它是下面凸 QP 的解：

\[ P_{\ell+1}=\arg\min_{P\ge 0}\ \tfrac12 P^\top W P + U_{\text{free}}^\top P. \]

（其 KKT 条件正是 LCP。）这个 QP 的解可以写成一个**近端算子（proximal operator）** 的形式——这正是 Moreau 当年发明的 Moreau 包络 / Moreau-Yosida 正则化：把不可微的指示函数 $\iota_{\{P\ge 0\}}$（可行锥的指示）用一个光滑的二次下界"磨光"，其极小点就是投影。

为什么这个视角重要？ 它揭示了三件事：

Moreau-Jean 的每步求解本质是"投影到摩擦锥"——这与 §5 的 CCP 投影梯度迭代 $\gamma^{k+1}=\Pi_{\mathcal{K}}(\gamma^k-\eta(N\gamma^k+r))$ 同源。LCP、CCP、proximal 求解都是"投影到接触可行集"的不同包装。
Pinocchio 3 的 ADMM + proximal 直接用 Moreau-Yosida 正则化处理 ill-conditioned Delassus（§3 提到的病态）——proximal 项 $\tfrac{\rho}{2}\|P-P^k\|^2$ 给 $W$ 加一个 $\rho I$，改善条件数。这就是为什么本章开头说"Pinocchio 3 ADMM 的 proximal 项对应 Moreau-Yosida 正则化"。
首尾呼应：Moreau 在 §2 用法锥写连续动力学，在这里他发明的 Moreau 包络又成为离散求解的工具——同一个人的两个发明（扫动过程 + Moreau 包络）在接触仿真里合流。

本质洞察：从"解 LCP"到"投影到锥"到"近端算子"，是同一件事的三种语言。认出这点，你就明白为什么 Moreau-Jean（LCP）、Anitescu CCP（投影梯度）、Pinocchio ADMM（proximal）看似不同的求解器，骨子里都在做"把无约束解投影回接触可行集"。Moreau-Yosida 正则化是这个投影的数学母体——它把硬约束（指示函数）磨成软的二次罚，极小点即投影。这也悄悄预告了 §6 MuJoCo 的软化：MuJoCo 的正则化 $R$ 本质上就是一个固定强度的 Moreau-Yosida 磨光。

桥接：Moreau-Jean 在工程中的位置¶

Siconos MoreauJeanOSI：完整实现 θ-scheme，是本格式的工业级参考。源码 kernel/src/simulationTools/MoreauJeanOSI.cpp。
Pinocchio 3 / Simple：其 ADMM 求解器的 proximal 项对应 Moreau-Yosida 正则化——这正是 Moreau 当年发明的工具（Moreau 包络），首尾呼应。
理解此节是读 Acary-Brogliato 2008（Springer LNACM 35）Ch. 9–10 的前提——那是非光滑时步法的"工程师百科全书"。

阶段小结：到这里我们完成了第一种离散——Moreau-Jean θ-scheme：速度级约束 + θ 加权位置更新 + 冲击律嵌入 argument，消元后每步解一个以 Delassus 算子 $W=HM^{-1}H^\top$ 为矩阵的 LCP；θ 控制能量耗散（$\tfrac12$ 中性、$1$ 耗散）。接下来 §4 讲第二种、工程上更常见的离散——Stewart-Trinkle 位置级 LCP，并给出"每步可解"的经典证明。

⚠️ 常见陷阱¶

💡 概念误区：以为 Moreau-Jean 的未知量是接触力 - 新手想法："互补条件 $0\le U\perp P\ge 0$ 里的 $P$ 是接触力吧？" - 实际上：$P$ 是**接触冲量**（量纲 $N\cdot s$），等于"接触力对一个时间步的积分"。持续接触时 $P\approx\lambda_n h$（力乘步长），冲击时 $P$ 是有限的 atom 值（力趋无穷但积分有限）。**统一用冲量**正是 Moreau-Jean 能"一套公式处理持续接触 + 冲击"的关键——若用力，冲击时力是无穷大无法表示。 - 为什么重要：调试时若把 $P$ 当力会发现量纲不对、数值差 $1/h$；正确理解后，要回到力只需除以 $h$（持续接触段）。

🧠 思维陷阱：认为 $\theta$ 越大越稳定越好 - 新手想法："$\theta=1$（全隐式）最稳定，那就总用 $\theta=1$。" - 实际上：$\theta=1$ 引入 $O(h^2)$ 数值耗散——它把系统能量人为抽走。对需要长期能量守恒的场景（卫星、摆、弹性碰撞链），这会让仿真"越跑越慢"直至静止，完全失真。对硬接触防 jitter 场景，$\theta=1$ 的耗散又是优点。没有普适最优，取决于物理需求。 - 正确思维：长期守恒/弹性 → $\theta=\tfrac12$；硬接触/抑制高频振荡 → $\theta\to 1$。把 θ 当"耗散旋钮"而非"稳定性开关"。

💡 概念误区：忽视位置漂移（drift） - 新手想法："速度级约束保证了 $v^+$ 合法，那位置肯定也不穿透。" - 实际上：速度级约束只保证**速度方向**合法（不侵入），但位置 $\phi(q_{\ell+1})$ 可能因离散误差缓慢变负（drift），物体逐渐陷入地面。Moreau-Jean 需配合**位置投影**（如 CombinedProjectedMoreauJeanOSI）或 §9 的稳定化技术修正。这是速度级格式的固有代价（对比 Stewart-Trinkle 的位置级严格不穿透）。 - 为什么重要：不加投影的纯 Moreau-Jean 长时间仿真会出现"物体慢慢沉入桌面"的视觉 bug——这不是 bug，是速度级格式的数学特性。

练习¶

（推导题，草稿纸完成） 从 $M(v_{\ell+1}-v_\ell)=hF_{\ell+1/2}+H^\top P$ 出发，设 $U_{\ell+1}=H v_{\ell+1}+eHv_\ell$，$W=HM^{-1}H^\top$。完整推导 $U_{\ell+1}=WP+U_{\text{free}}$ 的表达式，写出 $U_{\text{free}}$ 的显式形式。再证明 $W\succeq 0$，并说明这如何保证法向 LCP 可解（引用本目录 10 的 LCP 可解性定理）。
（推导题） 推导 θ-method 的能量耗散估计：对无接触光滑段，从 $\Delta E=\tfrac12(v_{\ell+1}^\top M v_{\ell+1}-v_\ell^\top M v_\ell)$ 出发，代入 $M(v_{\ell+1}-v_\ell)=hF_{\ell+\theta}$，证明 $\theta=\tfrac12$ 时 $\Delta E=h v_{\ell+1/2}^\top F$（能量中性），$\theta=1$ 时多出 $-\tfrac{h^2}{2}\|F\|_{M^{-1}}^2$（耗散 $O(h^2)$）。
（开放思考题） Moreau-Jean 用速度级约束导致位置漂移，Stewart-Trinkle 用位置级约束严格不穿透但需展开非线性 gap。设计一个数值实验（描述即可，不必编码）来定量比较两者在"方块静置桌面 1000 步"后的穿透量与能量漂移。你预期哪个穿透更小？哪个能量更稳定？为什么？

§4 Stewart-Trinkle 时步格式：位置级 LCP ⭐⭐⭐¶

动机：工程上最有影响力的离散格式¶

Moreau-Jean 数学最严谨，但 Stewart-Trinkle（1996） 在工程上影响最大——早期 Drake（TimeSteppingRigidBodyPlant）、Bullet 的 MLCP、ODE 都源于它。它的核心贡献有两个：(i) 把摩擦锥**离散成多面体锥**，使整个接触问题变成一个标准 LCP（而非更难的 NCP）；(ii) Anitescu-Potra（1997）证明这个 LCP 永远可解——这是接触仿真"总能算出一个解"的第一个严格保证。本节要手写出这个 LCP 矩阵，并复现可解性证明的骨架。

反面：为什么不能直接用二阶锥摩擦¶

Coulomb 摩擦的"真实"形状是**二阶锥（SOC）** $\{(\lambda_n,\beta):\|\beta\|_2\le\mu\lambda_n\}$（本目录 20）。直接用 SOC，每步要解一个**二阶锥互补问题（SOCCP）** 或非线性互补问题（NCP）——这在 1996 年缺乏成熟、可证可解的算法。

Stewart-Trinkle 的工程妥协：把圆锥的底面圆离散成 $d$ 条方向的多边形（pyramidal linearization），摩擦力写成这 $d$ 个方向基的非负组合 $\beta\ge 0$。这样切向约束变成线性不等式，整个问题落入**线性互补 LCP** 的成熟框架（Lemke 算法可解）。代价是摩擦锥形状有 $O(1/d)$ 的离散误差（$d$ 越大越接近圆锥，但 LCP 维度越高）。

历史：Stewart-Trinkle 与 Anitescu-Potra¶

David Stewart 与 Jeff Trinkle 在 IJNME 39:2673 (1996) 提出位置级 LCP 时步格式，标题强调"inelastic collisions and Coulomb friction"。Mihai Anitescu 与 Florian Potra 在 Nonlinear Dynam. 14:231 (1997) 给出可解性证明（copositive-plus）。Stewart 又在 ARMA 145:215 (1998) 证明了 $h\to 0$ 的收敛性，并由此**消解了 Painlevé 悖论**（见后）。这三篇构成 Stewart-Trinkle 路线的"三部曲"。

理论：Stewart-Trinkle LCP 的构造¶

Step 1：位置级非穿透约束。 Stewart-Trinkle 直接要求下一步不穿透：

\[ \phi(q_{\ell+1})\ge 0. \]

由于 $\phi$ 非线性，实践中线性化为 $\phi_\ell + h\,\nabla\phi^\top v_{\ell+1}\ge 0$，即"用新速度外推一步后 gap 仍非负"。这与 Moreau-Jean 的速度级 $\nabla\phi^\top v_{\ell+1}\ge 0$ 的区别在于**多了 $\phi_\ell/h$ 这一项**——它"看见"了当前 gap，能在物体还没碰到时就预判。

Step 2：多面体摩擦锥。 设接触法向 $n$，切向用 $d$ 个方向基 $D=[d_1,\dots,d_d]$（在切平面内均匀分布，且关于原点对称，即每个 $d_i$ 都有 $-d_i$）。摩擦冲量写成 $D\boldsymbol\beta$，$\boldsymbol\beta\ge 0$。Coulomb 约束 $\sum_i\beta_i\le\mu\lambda_n$。

Step 3：三组互补条件。 Stewart-Trinkle LCP 有三组未知量与三组互补：

法向：$0\le \lambda_n\ \perp\ (\text{法向分离速度}\ge 0)$；
切向：$0\le \boldsymbol\beta\ \perp\ (\lambda\,E+D^\top v_{\ell+1}\ge 0)$，其中 $\sigma=$ 滑移率，$E=[1,\dots,1]^\top$；
摩擦锥边界：$0\le \sigma\ \perp\ (\mu\lambda_n-E^\top\boldsymbol\beta\ge 0)$。

Step 4：组装成混合 LCP（MCP）。 把三组写成矩阵形式，未知量 $(\lambda_n,\boldsymbol\beta,\sigma)$（或含 $v_{\ell+1}$ 的扩展形式）满足带如下矩阵的 MCP：

\[ \boxed{\,M_{\text{LCP}}=\begin{bmatrix} W_{nn} & W_{nd} & 0 \\ W_{dn} & W_{dd} & E \\ \mu & -E^\top & 0 \end{bmatrix},\,} \]

其中： - $W=HM^{-1}H^\top$ 是 Delassus 算子（§3 的核心对象，在这里分块为法向-法向 $W_{nn}$、法向-切向 $W_{nd}$ 等）； - $E=[1,\dots,1]^\top$ 是把 $d$ 个切向方向耦合到同一滑移率 $\sigma$ 的全 1 向量； - 最后一行 $[\mu,-E^\top,0]$ 编码摩擦锥边界 $\mu\lambda_n\ge E^\top\boldsymbol\beta$； - $\sigma\ge 0$ 是**滑移率乘子（slip rate multiplier），互补条件 $\mu\lambda_n-E^\top\boldsymbol\beta\ge 0\perp\sigma\ge 0$ **复现最大耗散原理——粘滞时 $\sigma=0$（摩擦锥内部），滑动时 $\sigma>0$（在锥边界，摩擦力取最大）。

逐块物理解读： - 左上 $2\times 2$ 块 $\begin{bmatrix}W_{nn}&W_{nd}\\W_{dn}&W_{dd}\end{bmatrix}=DM^{-1}D^\top$ 型：对称半正定（Delassus 的子块），表示法向+切向冲量如何耦合产生接触点速度。 - $E$ 与 $-E^\top$ 的反对称配对：这是滑移率乘子 $\sigma$ 与切向方向的耦合，它在二次型里**反对称抵消**——这是下面 copositive-plus 证明的关键。 - $\mu$ 在左下角：摩擦系数把法向力"放大"为可用的切向力预算。

一个具体的 LCP 矩阵：2D 方块在地面¶

把抽象矩阵落到数字。考虑一个 2D 方块（质量 $m=1$，转动惯量 $I=1$）的一个角点接触地面，单接触。法向 $n=(0,1)$（向上），切向用两个对称方向 $d_1=(1,0)$、$d_2=(-1,0)$（$d=2$ 方向）。设接触点恰在质心正下方（简化，接触 Jacobian 不含转动耦合），则接触 Jacobian（法向 + 两切向）为

\[ H=\begin{bmatrix} n^\top \\ d_1^\top \\ d_2^\top \end{bmatrix}=\begin{bmatrix} 0 & 1 \\ 1 & 0 \\ -1 & 0 \end{bmatrix}\ (\text{仅平动 }2\text{ 自由度}),\qquad M^{-1}=\begin{bmatrix}1&0\\0&1\end{bmatrix}. \]

Delassus 子块 $W=HM^{-1}H^\top$：

\[ W=\begin{bmatrix} 0&1\\1&0\\-1&0\end{bmatrix}\begin{bmatrix} 0&1&-1\\1&0&0\end{bmatrix}=\begin{bmatrix} 1&0&0\\0&1&-1\\0&-1&1\end{bmatrix}. \]

解读：$W_{nn}=1$（法向有效柔度），切向 $2\times 2$ 块 $\begin{bmatrix}1&-1\\-1&1\end{bmatrix}$（半正定，秩 1——因为 $d_1=-d_2$ 两方向共线，对称化后冗余），法向-切向块为零（接触点在质心下方，法向冲量不产生切向速度）。

完整 $4\times 4$ Stewart-Trinkle 矩阵（未知量 $(\lambda_n,\beta_1,\beta_2,\sigma)$，最后一行/列是滑移率与摩擦锥边界 $\mu\lambda_n\ge\beta_1+\beta_2$，取 $\mu=0.5$）：

\[ M_{\text{LCP}}=\begin{bmatrix} 1&0&0&0\\ 0&1&-1&1\\ 0&-1&1&1\\ 0.5&-1&-1&0\end{bmatrix}. \]

验证不对称：$(M_{\text{LCP}})_{14}=0$ 而 $(M_{\text{LCP}})_{41}=0.5$；$(M_{\text{LCP}})_{24}=1$ 而 $(M_{\text{LCP}})_{42}=-1$——确实不对称。这正是为什么不能直接用对称正定求解器，而需要 copositive-plus 理论保证可解。来自摩擦锥边界行的元素是第 4 行 $[0.5,-1,-1,0]$（$\mu=0.5$ 与两个 $-E^\top=-1$）和切向列里的 $+E$ 项（第 2、3 行第 4 列的 $+1$）。亲手验证 $x^\top M_{\text{LCP}}x$ 在 $x\ge 0$ 时 $\ge 0$：取 $x=(\lambda_n,\beta_1,\beta_2,\sigma)$，展开发现 $\sigma$ 相关的交叉项 $\sigma(\beta_1+\beta_2)$（来自第 2、3 行 $+\sigma$）与 $\sigma(-\beta_1-\beta_2)$（来自第 4 行 $-E^\top$）恰好抵消，余下 $\lambda_n^2+(\beta_1-\beta_2)^2\ge 0$——这就是 copositive-plus 的反对称抵消机制的数字版！

核心定理：Anitescu-Potra copositive-plus¶

现在到 Stewart-Trinkle 路线最重要的理论结果——每步 LCP 一定有解。

定义（copositive-plus 矩阵）：矩阵 $N$ 称为 copositive-plus，若 $$ (\text{i}) x\ge 0\Rightarrow x^\top N x\ge 0;\qquad (\text{ii}) x\ge 0, x^\top N x=0\Rightarrow (N+N^\top)x=0. $$ 即在非负象限上半正定，且二次型取零时对称部分把 $x$ 打到零空间。

定理（Anitescu-Potra 1997）：Stewart-Trinkle LCP 矩阵 $M_{\text{LCP}}$ 是 copositive-plus。

证明骨架（复现 Anitescu-Potra 1997 §3 的思路）：

把 $M_{\text{LCP}}$ 分块为 $\begin{pmatrix}\tilde W & b\\ b^\top & 0\end{pmatrix}$ 形式，其中 $\tilde W=\begin{bmatrix}W_{nn}&W_{nd}\\W_{dn}&W_{dd}\end{bmatrix}\succeq 0$（因为它是 $H M^{-1}H^\top$ 的对称子块，$M^{-1}\succ 0$），而最后一行/列含摩擦项 $[\mu,-E^\top]$ 与 $[0;E]$。

对任意 $x=(x_1,x_2)\ge 0$（$x_1=(\lambda_n,\boldsymbol\beta)$，$x_2=\sigma$），计算二次型：

\[ x^\top M_{\text{LCP}}x = x_1^\top\tilde W x_1 + x_2(\mu\lambda_n-E^\top\boldsymbol\beta) + x_2(E^\top\boldsymbol\beta)\cdot(\text{来自 }+E\text{ 项})+\dots \]

关键观察：摩擦行 $[\mu,-E^\top,0]$ 与切向列里的 $+E$ 项**反对称配对**——它们在二次型 $x^\top M x=\tfrac12 x^\top(M+M^\top)x$ 中的反对称部分**互相抵消**（$+E$ 与 $-E^\top$ 贡献符号相反）。于是

\[ x^\top M_{\text{LCP}}x = x_1^\top \tilde W x_1 \ge 0, \]

证得 (i) copositive。

对 (ii)：若 $x^\top M_{\text{LCP}}x=x_1^\top\tilde W x_1=0$，由 $\tilde W\succeq 0$ 得 $\tilde W x_1=0$；结合反对称项消失，可推出 $(M_{\text{LCP}}+M_{\text{LCP}}^\top)x=0$。证毕。

推论（每步可解）：由 Cottle-Pang-Stone 定理（LCP 理论的核心结果，本目录 10），copositive-plus 矩阵的 LCP 若可行（feasible）则有解。Stewart-Trinkle LCP 总是可行（$\boldsymbol\beta=0,\lambda_n=0,\sigma=0$ 配合自由运动给出可行点的退化情形可构造），故**每步 LCP 恒有解**。

本质洞察：copositive-plus 是比"正定"弱得多的条件——$M_{\text{LCP}}$ 并不正定（它甚至不对称！），但 Anitescu-Potra 证明了它在非负象限上的"半正定 + 退化可控"足以保证 LCP 可解。这是**第一个不需要 $M_{\text{LCP}}\succ 0$ 的可解性证明**。它告诉我们：接触仿真"总能算出一个解"不是工程运气，而是摩擦锥多面体化后矩阵结构的数学必然。反对称的滑移率耦合（$E$ vs $-E^\top$）是这个奇迹的技术核心。

收敛性与 Painlevé 悖论的消解¶

Stewart 1998（ARMA 145:215, DOI 10.1007/s002050050129）收敛定理（简版）：

设 $M(q)\succ 0$ Lipschitz、外力 $F$ 连续、摩擦系数 $\mu$ 有界。令步长 $h\to 0$，Stewart-Trinkle 离散解 $(q_h,v_h)$ 存在**子列弱 * 收敛**到 Moreau MDI 的解，且冲量测度 $P_h\rightharpoonup P$ 于 $C([0,T])^*$。

证明链条（属于档位 4，这里给骨架）：离散 BV 估计（一致变差有界）→ Helly 选择定理（BV 函数列有逐点收敛子列）→ 弱 * 紧致（冲量测度列有弱 * 收敛子列）→ 极限满足 MDI。

Painlevé 悖论的消解。 Painlevé 1895 发现：某些刚体 + Coulomb 摩擦配置（如一根斜杆在地面滑动，摩擦系数大）下，经典刚体方程**无解或解不唯一**——接触力按经典理论会发散到无穷。这曾被认为是刚体模型的根本缺陷。

Stewart 1998 的推论：Stewart-Trinkle 离散解的 $h\to 0$ 极限存在，但含冲量 atom——即所谓 "impact without collision"（无碰撞的冲击）。物理图景：在 Painlevé 配置下，系统会自发产生一个瞬时冲量（即使没有"撞上"任何东西），把系统推离病态配置。时步法**自动**给出这个含 atom 的测度解，无需任何特殊处理。

本质洞察：Painlevé 悖论之所以是"悖论"，是因为人们坚持解必须是光滑函数（接触力是有限函数）。一旦接受解可以是**测度**（含冲量 atom），悖论消失——时步法的 BV/测度框架天然容纳了"无碰撞冲击"。这再次印证 §0–§2 的主线：接触动力学的解活在 BV/测度空间，用 $C^1$ 框架看它必然遇到"悖论"，换对空间就豁然开朗。

Stewart-Trinkle vs Moreau-Jean：位置级 vs 速度级的再对比¶

	Stewart-Trinkle	Moreau-Jean
约束级别	位置级 $\phi(q_{\ell+1})\ge 0$	速度级 $\nabla\phi^\top(v_{\ell+1}+ev_\ell)\ge 0$
穿透	严格不穿透（线性化后近似）	速度合法、位置可能漂移
gap 评估	需要（非线性 $\phi$ 展开）	不需要
摩擦锥	多面体 → LCP	切锥法锥（可 SOC 或多面体）
可解性	copositive-plus（Anitescu-Potra）	$W\succeq 0$ 半正定
漂移	无（位置级直接约束）	需位置投影
典型实现	早期 Drake、Bullet MLCP、ODE	Siconos MoreauJeanOSI

二者在持续接触下等价，差异出现在"即将接触但尚未接触"的过渡：Stewart-Trinkle 因为有 $\phi_\ell/h$ 项能"提前刹车"，Moreau-Jean 速度级则要等到 $\phi=0$ 才激活。用错（混淆两者）会导致 $O(h)$ 漂移或虚假冲击。

阶段小结：到这里我们完成了第二种离散——Stewart-Trinkle 位置级 LCP：多面体摩擦锥 + 滑移率乘子复现最大耗散，LCP 矩阵 copositive-plus 保证每步可解（Anitescu-Potra），$h\to 0$ 弱 * 收敛到 MDI 并消解 Painlevé 悖论（Stewart 1998）。LCP 的代价是激活集切换处不可微、多面体锥有 $O(1/d)$ 误差。下一步 §5 用 Anitescu 凸松弛把 LCP 升级为凸 QP/CCP——既改善可微性，又支持二阶锥。

⚠️ 常见陷阱¶

💡 概念误区：以为多面体摩擦锥是"近似"，会引入大误差 - 新手想法："把圆锥砍成 $d$ 边形肯定不准，$d$ 小时误差很大。" - 实际上：多面体化的误差是 $O(1/d^2)$（对锥半径），$d=8$ 已相当准。更重要的是它带来一个**实质性结构改变**：摩擦力方向被限制在 $d$ 个离散方向上，这会产生轻微的**方向各向异性**（沿基方向滑动比沿对角方向略易/略难）。这不是"精度"问题而是"形状"问题——$d\to\infty$ 才完全各向同性。 - 为什么重要：若你的任务对滑动方向敏感（如精密装配的旋转对中），多面体各向异性可能比锥半径误差更要命，此时应转向 §5/§6 的二阶锥（SOC）格式。

🧠 思维陷阱：把 LCP 可解性当成"解唯一" - 新手想法："Anitescu-Potra 证明了 LCP 有解，那解肯定唯一，仿真是确定性的。" - 实际上：copositive-plus 只保证**存在**解，不保证唯一。多体同时冲击（三球同碰、方块四角同触）下 LCP 可以有**多个解**，不同求解器（Lemke vs PGS）可能给出不同结果。这是接触仿真"非确定性"的一个来源——换求解器、换枢轴规则，结果可能变。 - 正确思维：接触 LCP 存在 ≠ 唯一。唯一性需更强条件（如 $M_{\text{LCP}}$ 是 P-矩阵）。多体同时冲击的非唯一性是病态，需 §7 的能量一致投影或 Frémond 多体冲击律处理。

💡 概念误区：混淆"位置级线性化"与"速度级约束" - 新手想法："Stewart-Trinkle 线性化成 $\phi_\ell+h\nabla\phi^\top v_{\ell+1}\ge 0$，这不就是速度级约束吗？" - 实际上：关键差异是那个 $\phi_\ell$ 项。速度级（Moreau-Jean）是 $\nabla\phi^\top v_{\ell+1}\ge 0$（无 $\phi_\ell$），只在 $\phi=0$ 时激活；位置级线性化是 $\phi_\ell/h+\nabla\phi^\top v_{\ell+1}\ge 0$，当 $\phi_\ell>0$ 时约束更松（允许接近），$\phi_\ell<0$（已穿透）时**主动把物体推出来**（修正穿透）。这个 $\phi_\ell/h$ 项让 Stewart-Trinkle 自带穿透修正。 - 为什么重要：理解这个差异，你才知道为什么 Stewart-Trinkle 不需要额外的 Baumgarte 稳定化（位置级约束内禀地修正穿透），而纯速度级 Moreau-Jean 需要。

练习¶

（推导题，草稿纸完成） 给定单接触、2D（一个法向 + 两个切向方向 $d_1=-d_2$），$M=\mathrm{diag}(m,m,I)$，接触 Jacobian $H$ 已知。手写出 $3\times 3$ 的 Stewart-Trinkle LCP 矩阵 $M_{\text{LCP}}$ 的每个元素（用 $W=HM^{-1}H^\top$ 的分量表示）。验证它不对称，并指出哪些元素来自摩擦锥边界行。
（证明题） 复现 Anitescu-Potra copositive-plus 证明：对 $x=(\lambda_n,\boldsymbol\beta,\sigma)\ge 0$，证明 $x^\top M_{\text{LCP}}x=x_1^\top\tilde W x_1\ge 0$（其中 $\tilde W$ 是法向+切向 Delassus 子块），关键是说清摩擦行 $[\mu,-E^\top]$ 与切向列 $+E$ 如何反对称抵消。再补全 (ii)：$x^\top M x=0\Rightarrow(M+M^\top)x=0$。
（开放思考题） Painlevé 悖论中"无碰撞的冲击"（impact without collision）听起来违反直觉——没撞上任何东西怎么会有冲量？请用斜杆滑动的物理图景解释：当杆以大摩擦系数滑动、且某临界配置下"摩擦力越大→法向力越大→摩擦力更大"形成正反馈时，系统如何自发产生瞬时冲量逃离病态？这个冲量在能量上是耗散还是注入？

§5 Anitescu 凸松弛：从 LCP 到凸 QP/CCP ⭐⭐⭐¶

动机：LCP 的两个痛点，凸优化的两个礼物¶

§4 的 Stewart-Trinkle LCP 很美，但有两个工程痛点：

不可微：LCP 的解在激活集（哪些接触在压紧）切换处**不可微**——这对 RL/MPC 梯度（本目录 40）是致命的。
求解算法的可扩展性：Lemke 等枢轴算法在大规模（$10^5$ 接触）时不如迭代法（投影梯度、ADMM）可扩展，且 LCP 一般不对称矩阵难以用对称正定求解器加速。

Mihai Anitescu（2006, Math. Program. 105:113） 的关键贡献：通过一个**凸松弛**，把每步的接触问题从 LCP（非凸可行集 / 不对称矩阵）变成一个**严格凸二次规划（QP）** 或等价的**锥互补问题（cone complementarity problem, CCP）**。凸优化带来两个礼物：(i) 全局最优解唯一且对参数连续依赖（→ 可微，本目录 40）；(ii) 可用高效迭代法（投影梯度 APGD、ADMM、Newton）求解，扩展到百万接触。代价是引入一个 $O(h)$ 的人工滑移（gliding/creep）。本节讲清这个松弛"松"在哪、代价是什么。

反面：为什么原始问题非凸¶

Coulomb 摩擦的"精确"速度级互补问题为什么非凸？核心在最大耗散原理的**双线性耦合**：摩擦力 $\beta$ 的可用预算 $\mu\lambda_n$ 依赖法向力 $\lambda_n$，而 $\lambda_n$ 和 $\beta$ 又同时是未知量。这个 "$\|\beta\|\le\mu\lambda_n$ 且 $\lambda_n,\beta$ 耦合" 的可行集，连同互补条件，构成一个**非凸**的解集（Anitescu-Hart 指出 impulse-velocity time-stepping 虽总有解但解集可能非凸）。非凸 ⇒ 多解、求解器可能卡在局部、解对参数不连续（不可微）。

历史：从 Stewart-Trinkle 到 Anitescu 到 Tasora-Chrono¶

Anitescu 的凸松弛思路源于"用凸优化的良好性质换取一点物理精度"。2006 年 Math. Program. 论文给出严格凸 QP 表述与收敛证明。Tasora-Anitescu（2010, 2011） 把它发展为锥互补 CCP，用 APGD（加速投影梯度，Nesterov 型） 求解，落地为 Project Chrono 的 granular 仿真引擎（$10^6$ 颗粒）。这条路线（凸 primal）后来成为 2020s 主流：Drake SAP、Pinocchio 3 ADMM 都属此谱系。

理论：Anitescu 凸松弛的构造¶

Step 1：速度级动力学（同前）。 一步内（步长 $h$）：

\[ M(v_{\ell+1}-v_\ell)=h F_\ell + \sum_i\bigl(n_i\lambda_n^i + D_i\boldsymbol\beta^i\bigr), \]

其中 $n_i$ 法向、$D_i=[d_1,\dots,d_d]$ 切向方向基，$\lambda_n^i\ge 0$、$\boldsymbol\beta^i\ge 0$ 是法向/切向冲量。

Step 2：松弛的非穿透约束。 这是凸化的第一处关键。Anitescu 不写严格的 $\phi(q_{\ell+1})\ge 0$，而是写**速度级 + gap 稳定项**：

\[ n_i^\top v_{\ell+1} + \frac{\phi_i(q_\ell)}{h}\ \ge\ 0, \]

即"法向分离速度 + 当前 gap/步长 $\ge 0$"。这个 $\phi_i/h$ 项把位置级信息注入速度级约束（与 Stewart-Trinkle 线性化同源），但 Anitescu 进一步把它写进一个**凸**的框架。

Step 3：凸松弛的核心——把双线性耦合改成凸约束。 原始最大耗散是 $\beta$ 关于 $\lambda_n$ 的双线性，非凸。Anitescu 的松弛：把摩擦约束写成

\[ \mu_i\,(n_i^\top v_{\ell+1}) + \|D_i^\top v_{\ell+1}\|\ \le\ 0\quad\text{（在接触激活时）}, \]

或等价地，引入一个**附加项**使整个问题成为对**冲量**的凸优化。最终每步求解一个**严格凸 QP**：

\[ \boxed{\,\min_{\boldsymbol\gamma}\ \tfrac12\,\boldsymbol\gamma^\top N\,\boldsymbol\gamma + \mathbf{r}^\top\boldsymbol\gamma\quad\text{s.t.}\quad \boldsymbol\gamma\in\mathcal{K}\ (\text{摩擦锥，二阶锥或多面体}),\,} \]

其中 $\boldsymbol\gamma=(\lambda_n,\boldsymbol\beta)$ 是接触冲量，$N=D^\top M^{-1}D$ 是 Delassus 算子（§3 的 $W$！只是把法向+切向方向都堆进 $D$），$\mathbf{r}$ 含自由速度与 gap 稳定项 $\phi/h$。

Step 4：CCP 形式（等价对偶）。 凸 QP 的 KKT 条件等价于一个**锥互补问题（CCP）**：

\[ \mathcal{K}\ni\boldsymbol\gamma\ \perp\ (N\boldsymbol\gamma+\mathbf{r})\in\mathcal{K}^*, \]

其中 $\mathcal{K}^*$ 是摩擦锥的对偶锥。这就是 Tasora-Anitescu 的 "cone complementarity problem" 名称由来。CCP 可用投影梯度迭代求解：$\boldsymbol\gamma^{k+1}=\Pi_{\mathcal{K}}(\boldsymbol\gamma^k-\eta(N\boldsymbol\gamma^k+\mathbf{r}))$，其中 $\Pi_{\mathcal{K}}$ 是到摩擦锥的投影（解析可算）。

Step 5：严格凸性 ⇒ 唯一解。 由 $M\succ 0$，$N=D^\top M^{-1}D\succeq 0$；当 $D$ 列满秩（接触非冗余）时 $N\succ 0$，QP 严格凸，解唯一。即使 $N$ 仅半正定，目标函数加上摩擦锥约束后通常仍有唯一解（或可加 proximal 正则化保证，见 §6 的 ADMM）。

锥投影 $\Pi_{\mathcal{K}}$ 的解析形式与 APGD 迭代¶

CCP 投影梯度迭代的核心是**到摩擦锥的投影** $\Pi_{\mathcal{K}}$。这个投影有解析解，值得写出来——它是整个凸 primal 路线（Chrono APGD、Drake SAP）每次迭代的内层操作。

对单个接触的二阶锥 $\mathcal{K}=\{(\gamma_n,\boldsymbol\gamma_t):\|\boldsymbol\gamma_t\|\le\mu\gamma_n\}$，投影一个点 $(z_n,\mathbf{z}_t)$ 到 $\mathcal{K}$ 分三种情形（这是凸优化里二阶锥投影的标准结果）：

\[ \Pi_{\mathcal{K}}(z_n,\mathbf{z}_t)=\begin{cases} (z_n,\mathbf{z}_t) & \text{若}\ \|\mathbf{z}_t\|\le\mu z_n\ (\text{锥内，不动}),\\[4pt] (0,\mathbf{0}) & \text{若}\ \|\mathbf{z}_t\|\le-z_n/\mu\ (\text{对偶锥内，归零}),\\[4pt] \dfrac{\mu\|\mathbf{z}_t\|+z_n}{1+\mu^2}\Bigl(1,\ \mu\dfrac{\mathbf{z}_t}{\|\mathbf{z}_t\|}\Bigr) & \text{否则（投到锥面）}. \end{cases} \]

三种情形的物理意义：(1) 接触力已在摩擦锥内（合法的 stick），保持不变；(2) 接触力指向"负法向"（要把物体往里拉，不物理），归零（接触松开）；(3) 接触力超出锥（要求的摩擦超过 $\mu\gamma_n$），投影到锥面（slip，摩擦取最大）。这三种情形恰好对应物理的 stick / 分离 / slip 三态——锥投影把不可行的接触力"拉回"物理可行集。

APGD（加速投影梯度）迭代（Tasora-Anitescu 2011，Chrono 主力）：在朴素投影梯度上加 Nesterov 加速，收敛率从 $O(1/k)$ 提升到 $O(1/k^2)$：

\[ \begin{aligned} \boldsymbol\gamma^{k+1}&=\Pi_{\mathcal{K}}\bigl(\mathbf{y}^k-\eta(N\mathbf{y}^k+\mathbf{r})\bigr),\\ t_{k+1}&=\tfrac{1}{2}\bigl(1+\sqrt{1+4t_k^2}\bigr),\\ \mathbf{y}^{k+1}&=\boldsymbol\gamma^{k+1}+\tfrac{t_k-1}{t_{k+1}}(\boldsymbol\gamma^{k+1}-\boldsymbol\gamma^k), \end{aligned} \]

其中 $\eta\le 1/\|N\|$（步长，$\|N\|$ 是 Delassus 谱范数），$\mathbf{y}^k$ 是 Nesterov 外推点。每次迭代只需一次矩阵-向量乘 $N\mathbf{y}^k$（matrix-free 时甚至不显式构造 $N$）+ 一次解析锥投影——这就是它能扩展到 $10^6$ 接触的原因：每步成本线性于接触数。

本质洞察：凸 primal 路线的全部计算，归结为"矩阵-向量乘（传播接触耦合）+ 锥投影（拉回物理可行集）"的反复迭代。锥投影的三态（stick/分离/slip）正是物理的三态，而它**解析可算、处处定义**（即使在 stick/slip 边界也有唯一投影值）——这是凸性的馈赠。对比 LCP 的枢轴算法（Lemke）需要组合搜索激活集，投影梯度只需"乘 + 投影"，既快又可微（投影算子几乎处处可微）。这把可微性（本目录 40）和可扩展性（百万接触）一并交付。

凸松弛的代价：人工滑移（gliding / creep）¶

凸化不是免费的。Anitescu 松弛引入一个**物理上可见的伪影**：当法向力很大时，物体会产生一个微小的、与法向力成正比的**切向滑移**——即使物理上应该静止粘滞（stick）。

为什么会这样？ 松弛后的摩擦约束 $\mu(n^\top v)+\|D^\top v\|\le 0$ 把法向速度 $n^\top v$ 与切向速度耦合。在持续接触（$n^\top v\approx 0$）下这没问题，但由于 gap 稳定项 $\phi/h$ 和软化效应，法向方向有微小的"压入速度"，通过耦合泄漏成切向滑移。其量级是 $O(h)$（步长越小越轻）。

Drake SAP 文档里称之为 "$O(\delta t)$ gliding"——一个静置在斜面上、本应靠摩擦不动的物体，会以正比于步长的速度极缓慢下滑。这在大步长 MPC（$h\sim 5$–$10$ ms）下可能可见，需要 warm-start 或减小步长缓解。

本质洞察：Anitescu 凸松弛是一笔**精确的交易**——用 $O(h)$ 的人工滑移，换来全局唯一解、可微性、和百万级可扩展性。这不是"近似得不够好"，而是"主动选择放弃一点物理精确性以获得数值与优化的良好性质"。理解这一点，你就懂了 2020s 主流仿真器（Drake SAP、MuJoCo、Chrono、Pinocchio）的共同哲学：凸性 > 物理精确性，因为 RL/MPC/可微仿真对前者的需求远大于后者。

凸松弛与多面体锥：正交的两个近似¶

这里必须澄清一个高频混淆：Stewart-Trinkle 的多面体锥**和 **Anitescu 的凸松弛**是**两个正交（独立）的近似，可叠加可独立：

	近似的对象	近似类型	误差	代数后果
多面体锥（Stewart-Trinkle）	摩擦锥的形状（圆 → $d$ 边形）	几何离散	$O(1/d^2)$	LCP 结构精确（线性）
凸松弛（Anitescu）	约束本身（互补 → 凸不等式 + $\phi/h$）	凸化	$O(h)$ 人工滑移	凸 QP / CCP

二者关系： - 可叠加：Stewart-Trinkle 用多面体锥 + 不松弛 → LCP。 - 可独立：Anitescu 可用**二阶锥（SOC，精确形状）** + 凸松弛 → SOCCP（无多面体误差，但有 $O(h)$ 滑移）。 - 现代选择：Drake SAP 用凸松弛 + 二阶锥（**不用**多面体）；Siconos 两者都支持；Chrono APGD 用凸松弛 + 二阶锥。

本质洞察：很多人把"多面体化"和"凸松弛"混为一谈，以为都是"把摩擦搞简单"。其实它们改的是不同东西：多面体化改**锥的形状**（精度 $O(1/d)$，LCP 精确），凸松弛改**约束的数学结构**（精度 $O(h)$，换来凸性）。一个是几何近似，一个是分析近似。认清这点，你才能读懂"为什么 SAP 用 SOC 却仍是凸的"（凸松弛与 SOC 并不矛盾）。

桥接：凸性是可微接触的前提¶

这是本节通往本目录 40_可微接触的桥。为什么凸化如此重要？因为凸优化的解对参数可微。

LCP 的解 $\boldsymbol\gamma(\theta)$（$\theta$ 是物理参数如摩擦系数、质量）在激活集切换处**不可微**——梯度 $\partial\boldsymbol\gamma/\partial\theta$ 在那里跳变或不存在。
凸 QP 的解（在严格凸 + 解唯一时）对参数**连续可微**，梯度可由 KKT 条件的隐函数定理（implicit function theorem）解析计算——这正是可微仿真（diff-sim）的数学基础。

所以"凸松弛"不只是数值技巧，它是**可微接触的必要前提**。本目录 40 会讲：把接触求解器写成可微优化层（differentiable optimization layer），梯度通过 KKT 反向传播——而这只在凸（QP/CCP）情形下良定义。

Anitescu 2006 收敛定理：凸 QP 松弛的解序列，当步长 $h\to 0$（同时松弛参数 $\epsilon\to 0$）时，收敛到差分变分不等式（DVI）/ MDI 的解。即凸松弛在极限下回到正确物理——人工滑移随 $h\to 0$ 消失。这保证了凸化是"渐近一致"的，不是永久偏差。

阶段小结：到这里我们完成了从 LCP 到凸优化的升级——Anitescu 凸松弛把每步接触问题写成严格凸 QP（Hessian 是 Delassus 算子 $N=D^\top M^{-1}D$）或等价 CCP，用投影梯度/ADMM 求解，代价是 $O(h)$ 人工滑移，回报是唯一解 + 可微性 + 百万级扩展。凸松弛与多面体锥正交。这条凸 primal 路线是 2020s 主流。下一步 §6 看一个把"凸 + 软化"做到极致的工程典范——MuJoCo。

⚠️ 常见陷阱¶

💡 概念误区：以为凸松弛 = 多面体锥近似 - 新手想法："凸松弛就是把摩擦锥砍成多边形让问题变线性吧？" - 实际上：二者正交。多面体化是改**锥形状**（圆→多边形，几何近似，结果仍是 LCP/可能非凸）；凸松弛是改**约束结构**（让解集变凸，可用 SOC 精确形状）。Drake SAP 用凸松弛但**不用**多面体（它用精确二阶锥）。把两者混为一谈会让你读不懂"凸 + SOC"为何不矛盾。 - 为什么重要：选格式时要分别决定"锥形状"（SOC vs 多面体）和"是否凸松弛"——这是两个独立旋钮。

🧠 思维陷阱：把人工滑移当成 bug 去"修复" - 新手想法："物体在斜面上慢慢下滑，肯定是摩擦没设对，调大摩擦系数。" - 实际上：凸松弛的 $O(h)$ gliding 是**格式固有**的，不是参数错误。调大摩擦系数治标不治本——正确做法是减小步长 $h$（滑移正比于 $h$）或用 warm-start。把它当 bug 乱调参会引入新问题。 - 正确思维：看到正比于步长的微小滑移，先判断是否凸松弛伪影；若是，减小 $h$ 或接受它（多数 RL/MPC 任务能容忍 $O(h)$ 滑移）。

💡 概念误区：认为凸化丢失了物理，所以凸求解器不如 LCP 准 - 新手想法："LCP 是精确的互补，凸松弛是近似，所以 Stewart-Trinkle 比 SAP 更物理。" - 实际上：(1) Stewart-Trinkle 的多面体锥本身也是近似（锥形状）；(2) 凸松弛的误差 $O(h)$ 随步长消失（Anitescu 2006 证明收敛到 MDI）；(3) LCP 在多体同时冲击下解非唯一，反而"不确定"。凸格式用一点可控的滑移换来唯一性、可微性、稳定性——在机器人任务里通常**净收益为正**。 - 为什么重要：不要因为"凸=近似"就教条地拒绝它。2020s 主流全面转向凸 primal 正是因为综合权衡下它更优。

练习¶

（推导题） 从速度级动力学 $M(v_{\ell+1}-v_\ell)=hF_\ell+D\boldsymbol\gamma$ 出发，消去 $v_{\ell+1}$，证明每步凸 QP 的 Hessian 恰是 Delassus 算子 $N=D^\top M^{-1}D$，并写出线性项 $\mathbf{r}$ 中"自由速度"与"gap 稳定项 $\phi/h$"分别长什么样。说明为什么 $D$ 列满秩时 $N\succ 0$（严格凸、解唯一）。
（开放思考题） 解释 Anitescu 凸松弛中 $O(h)$ 人工滑移的物理来源：松弛后的摩擦约束如何把法向方向的微小压入速度"泄漏"成切向滑移？为什么这个滑移正比于步长 $h$ 而非摩擦系数？设计一个最小实验（斜面静置方块）来测量滑移速度对 $h$ 的依赖，验证线性关系。
（综合题，桥接本目录 40） 论证"凸性是可微接触的前提"：(a) 举例说明 LCP 解在激活集切换处为何不可微（画出一维接触力对某参数的曲线，指出折点）；(b) 说明凸 QP 解在严格凸时为何可微（提示：KKT 系统 + 隐函数定理）；(c) 这对 RL 策略梯度/MPC 灵敏度意味着什么？为什么 MuJoCo/Drake 能做 diff-sim 而 Bullet 难？

§6 MuJoCo 软接触求解器：连续时间 + 凸优化 ⭐⭐⭐¶

动机：为 RL 与可微仿真而生的另一条路¶

§3–§5 都是"刚性接触 + 时间离散"的不同格式。MuJoCo（Todorov et al. 2012）走了一条**显著不同**的路：它**主动放弃刚性接触**，采用一个**内禀柔性（inherently soft）** 的连续时间接触模型，把每步的接触求解写成一个**凸优化问题**。这个选择是为了 RL 训练（百万环境并行）和 model-based control（可微、可逆）量身定制的。本节是全章"工程落地"的核心——读懂它，你就读懂了机器人学界用得最多的仿真器的心脏。MuJoCo 文档（mujoco.readthedocs.io）的 Computation 一章是本节的一手依据。

反面：硬接触对 RL/控制的不友好¶

为什么 MuJoCo 不用 §4/§5 的硬接触 LCP/CCP？

硬接触不可微或难微分：LCP 解在激活集切换处不可微（§5），即使凸 CCP 在锥边界（stick/slip 切换）也只是次可微。RL/MPC 需要**处处光滑**的梯度。
硬接触的"接触/不接触"是二值开关：这给优化景观（optimization landscape）带来不连续，策略梯度方差大。
硬约束需要精确碰撞检测 + 互补求解，计算重、难以 GPU 大批量并行。

MuJoCo 的反向选择：让接触"软"——pushing harder against a constraint results in larger acceleration（推得越狠、加速度越大），接触力是 gap 的**光滑函数**而非互补开关。这把整个问题变成一个**处处可微的凸优化**，代价是刚性保真度不足（硬物体会有微小"海绵感"穿透）。

历史：Todorov 的 MuJoCo 与 Gauss 原理¶

Emanuel Todorov（华盛顿大学，后 DeepMind）2012 年发布 MuJoCo（Multi-Joint dynamics with Contact），核心论文 "MuJoCo: A physics engine for model-based control"。其接触模型的数学根基是 Gauss 最小约束原理（Gauss's principle of least constraint）：真实加速度是"在约束允许范围内、与无约束加速度偏差最小"的那个。Todorov 把 Gauss 原理 + 软约束正则化写成一个凸优化，这是 MuJoCo 区别于 LCP 学派的根本。后续（Todorov 2014 "Convex and analytically-invertible dynamics"）进一步证明了该动力学**解析可逆**——这对逆动力学/控制极有价值。

理论：MuJoCo 求解的凸优化问题¶

Step 1：Gauss 原理的原始（primal）形式。 MuJoCo 求解的核心优化（依据官方 Computation 文档）是

\[ \boxed{\,(\dot v)=\arg\min_{x}\ \tfrac12\bigl\|x-M^{-1}(\tau-c)\bigr\|_M^2 + \tfrac12\bigl\|J x - a_{\text{ref}}\bigr\|^{\text{(Huber)}}_{R^{-1}}\quad\text{s.t.}\quad Jx-y\in\mathcal{K}^*,\,} \]

逐项解读： - $M^{-1}(\tau-c)$：无约束加速度（$\tau$ 广义力，$c$ 科氏/重力项）。第一项 $\|x-M^{-1}(\tau-c)\|_M^2$ 就是 Gauss 原理——惩罚偏离无约束加速度（用质量度量 $\|\cdot\|_M$）。 - $a_{\text{ref}}$：参考加速度（reference acceleration），由约束"想要达到"的状态决定（如把穿透推回零）。 - $R$：对角正则化矩阵（regularizer），它让约束变"软"——$R$ 越大约束越软。这是 MuJoCo 软接触的数学开关。 - $\mathcal{K}^*$：摩擦锥的对偶锥。 - Huber 范数：摩擦损失项用 Huber 函数（二次-线性混合）使其光滑且有界。

Step 2：对偶（dual）形式——PGS 在此求解。 Lagrange 对偶是

\[ \boxed{\,f=\arg\min_\lambda\ \tfrac12\lambda^\top(A+R)\lambda + \lambda^\top(a_u-a_{\text{ref}})\quad\text{s.t.}\quad\lambda\in\Omega,\,} \]

其中： - $A=JM^{-1}J^\top$：约束空间的逆惯量——又是 Delassus 算子！（§3 的 $W$、§5 的 $N$ 的连续时间版本）。 - $a_u=JM^{-1}(\tau-c)+\dot Jv$：约束空间的无约束加速度。 - $\Omega$：可行的接触力集合（摩擦锥）。 - $R$ 加在 $A$ 上 ($A+R$)：正则化让 Hessian 严格正定 ⇒ 严格凸 ⇒ 唯一解。这个 $+R$ 就是软接触的代数体现——它保证即使接触冗余（$A$ 奇异）问题仍严格凸可解。

本质洞察：把 MuJoCo 对偶问题 $\min\tfrac12\lambda^\top(A+R)\lambda+\dots$ 与 §5 的 Anitescu 凸 QP $\min\tfrac12\gamma^\top N\gamma+\dots$ 并排看——它们是同一个东西！$A=JM^{-1}J^\top$ 就是 $N=D^\top M^{-1}D$（连续 vs 离散、加速度 vs 冲量的差别）。区别只在：Anitescu 用硬摩擦锥（$\gamma\in\mathcal{K}$，可能在边界不可微），MuJoCo 加正则化 $+R$ 并软化锥（处处可微）。所以 MuJoCo 不是"另一个物理"，而是 Anitescu 凸 primal 路线 + 软化正则化的一个特例。认出这点，四种格式（Moreau-Jean、Stewart-Trinkle、Anitescu、MuJoCo）就在 Delassus 算子 $W=JM^{-1}J^\top$ 下完全统一了。

Step 3：约束的软化——solref 与 solimp。 软接触的"软度"由两组参数控制，它们定义了 $a_{\text{ref}}$ 和 $R$（impedance）：

参考加速度满足一个**虚拟弹簧-阻尼器**动力学：

\[ a_{\text{ref}} = -k_p\,\phi - k_d\,\dot\phi, \]

其中 $\phi$ 是约束违反量（穿透），$(k_p,k_d)$ 是刚度/阻尼。

solref = $(\text{timeconst},\ \text{dampratio})$：定义这个虚拟弹簧-阻尼器的**时间常数**和**阻尼比**。它描述约束"想要做什么"（what the constraint wants to do）——以多快、多大阻尼把穿透推回零。重新参数化为质量-弹簧-阻尼系统的固有频率与阻尼比，比直接设 $(k_p,k_d)$ 更直观。
solimp = $(d_0,\ d_{\text{width}},\ \text{width},\dots)$：定义**阻抗（impedance）** $d\in[0,1]$ 如何随约束违反量变化——它描述求解器"能做到什么"（what the solver is able to do）。阻抗在无接触时小（软），违反增大时趋向 1（硬）。这让 MuJoCo 能建模"软接触层 + 硬芯"的复合材料感。

一句话记忆：solref 是"想做什么"（reference），solimp 是"能做到什么"（impedance）。两者共同决定 $a_{\text{ref}}$（参考加速度）和 $R$（正则化）：$\dot\omega=a_{\text{ref}}-Rf$ 是约束的形变动力学。

solref 的弹簧-阻尼器参数化算例。 MuJoCo 默认 solref = (timeconst, dampratio) = (0.02, 1.0)。它把约束建模为一个二阶质量-弹簧-阻尼器，固有频率 $\omega_0$ 和阻尼比 $\zeta$ 由这两个参数决定：

\[ \omega_0=\frac{1}{\text{timeconst}}=\frac{1}{0.02}=50\ \mathrm{rad/s},\qquad \zeta=\text{dampratio}=1.0\ (\text{临界阻尼}). \]

等效刚度与阻尼（对单位有效质量）：$k_p=\omega_0^2=2500$，$k_d=2\zeta\omega_0=100$。参考加速度 $a_{\text{ref}}=-k_p\phi-k_d\dot\phi=-2500\phi-100\dot\phi$。物理意义：当穿透 $\phi$（这里取约束违反为负，记 $|\phi|$）出现时，约束以时间常数 $0.02$ s、临界阻尼**无超调地**把违反推回零。临界阻尼 $\zeta=1$ 是默认选择，因为它最快回零且不振荡（$\zeta<1$ 欠阻尼会"弹"，$\zeta>1$ 过阻尼回得慢）。

为什么 timeconst 必须 $\ge 2h$？ 这是 MuJoCo 文档的硬性建议，背后是数值稳定性。半隐式 Euler 求解这个弹簧-阻尼器，稳定性要求约束的固有周期被步长充分分辨——若 timeconst $<2h$（即 $\omega_0>1/(2h)$，约束比步长还快），离散弹簧-阻尼器会**不稳定振荡**（接触力在压紧/松开间高频切换 → jitter，§9）。取 timeconst $\ge 2h$ 保证每个约束周期至少被几步分辨。例：$h=0.002$ s 时 timeconst 须 $\ge 0.004$ s；默认 $0.02$ s 远大于此，安全。

solimp 的阻抗曲线算例。 默认 solimp = (d0, dwidth, width, midpoint, power) = (0.9, 0.95, 0.001, 0.5, 2)。阻抗 $d\in[0,1]$ 随约束违反量从 $d_0=0.9$（接触刚开始，较软）平滑过渡到 $d_{\text{width}}=0.95$（深度穿透，较硬），过渡宽度 width$=0.001$ m，曲线形状由 power$=2$（二次）和 midpoint$=0.5$ 决定。阻抗 $d$ 直接调制正则化：$R\propto(1-d)/d$——$d\to 1$ 时 $R\to 0$（约束变硬，无正则化），$d$ 小时 $R$ 大（约束软）。这让 MuJoCo 能建模"表层软、深层硬"的复合接触（如橡胶垫包裹硬芯）：浅穿透时软（橡胶变形），深穿透时硬（碰到硬芯）。

Step 4：摩擦锥——pyramidal vs elliptic 与 impratio。 MuJoCo 支持两种摩擦锥：

椭圆锥（elliptic）：$\mathcal{K}_{\text{elliptic}}=\{f:f_1\ge 0,\ f_1^2\ge\sum_{i\ge 2}f_i^2/\mu_{i-1}^2\}$——精确的二阶锥（含切向、扭转 torsional、滚动 rolling 摩擦），各向同性。
金字塔锥（pyramidal）：$\mathcal{K}_{\text{pyramidal}}=\{f\in\mathbb{R}^{2(n-1)}:f\ge 0\}$——多面体近似（§4 的 Stewart-Trinkle 风格），用棱向量，计算更快但有各向异性。
impratio：控制法向 vs 切向阻抗之比。增大 impratio 让摩擦"更硬"（更接近无滑），常用于抓取/装配需要强摩擦的场景。

Step 5：三个求解器 PGS / CG / Newton。 同一个凸优化问题，MuJoCo 提供三种数值算法（依据官方文档）：

求解器	工作空间	收敛	特点	适用
PGS（Projected Gauss-Seidel）	对偶	一阶（sublinear）	逐约束更新到最优；前几次迭代快；椭圆锥时用 QCQP 子步（先沿锥射线优化，再在超平面切片椭球内优化）	容忍不精确解、追求速度
CG（Conjugate Gradient）	原始（reduced primal）	高阶	非线性共轭梯度（Polak-Ribière-Plus）+ 精确线搜索（分段二次代价上解析二阶导）；无 setup 开销	数值更难但仍凸的问题
Newton	原始	二阶	精确 Newton + 解析二阶导 + Hessian 的 Cholesky 分解 + 精确线搜索；约束状态变化时增量更新分解；默认求解器	高精度、中等规模

三者都支持金字塔/椭圆锥、稠密/稀疏 Jacobian。Newton 是默认——对典型机器人模型，它几步内收敛到全局最优。PGS 的 10 次扫描通常已与全局最优无法区分（对典型模型）。

半隐式 Euler 与积分器选项¶

MuJoCo 的时间积分用**半隐式 Euler（semi-implicit Euler）** 为主，另有 implicitfast（对速度相关项隐式）和 RK4（光滑段高阶）选项。半隐式 Euler 先更新速度（用新接触力）再更新位置：

\[ v_{\ell+1}=v_\ell+hM^{-1}(\tau-c+J^\top f),\qquad q_{\ell+1}=q_\ell+h\,v_{\ell+1}. \]

注意位置用**新速度** $v_{\ell+1}$（这是"半隐式/symplectic Euler"的标志），比显式 Euler 稳定得多。implicitfast 进一步把阻尼等速度相关项隐式化，提升大步长稳定性。

桥接：为什么 MuJoCo 选这条路¶

MuJoCo 的全部设计决策，都指向**可微性 + 大批量并行**：

软接触 ⇒ 处处可微：接触力是 gap 的光滑函数，无激活集开关，梯度处处良定义 → MJX（JAX 版）可微 + 百万环境并行，这是 RL 的杀手锏。
凸优化 ⇒ 唯一解 + warm-start：每步唯一全局最优，可用上一步解热启动，收敛快。
连续时间 soft-convex ⇒ 时间可逆：Todorov 2014 证明动力学解析可逆，逆动力学（已知加速度求力）有闭式——对 model-based control 极有价值。
代价：刚性保真不足：硬物体有微小穿透（"海绵感"），对需要精确刚性接触力的任务（如精密力控）不如 Drake SAP/Siconos。

本质洞察：MuJoCo 不追求"最物理精确"，而追求"对学习与控制最友好"。它用**软化 + 凸化**这两把刷子，把一个本质上非光滑、不可微、难并行的接触问题，刷成了一个**光滑、可微、可大批量 GPU 并行**的凸优化。这就是为什么 RL 社区几乎人手一个 MuJoCo——它牺牲的刚性保真度，恰好是 RL 不在乎的；它换来的可微性与并行性，恰好是 RL 最需要的。一个工具的成功，往往不在于它"最准"，而在于它的取舍**恰好对上了一个领域的核心需求**。

阶段小结：到这里我们完成了"工程落地"的核心——MuJoCo 软接触求解器：Gauss 原理 + 正则化 $R$ 的凸优化，对偶 Hessian 是 Delassus 算子 $A=JM^{-1}J^\top$（与 §5 的 $N$ 同构），solref/solimp 控制软度，PGS/CG/Newton 三求解器，半隐式 Euler 积分；软化+凸化换取可微性与百万并行。至此四种格式（§3–§6）在 Delassus 算子下统一。接下来 §7–§9 处理细节：冲击律、辛性、穿透/能量/漂移。

⚠️ 常见陷阱¶

💡 概念误区：以为 MuJoCo 是"硬接触 + 碰撞检测" - 新手想法："MuJoCo 跟 Bullet 一样，检测碰撞再算接触力。" - 实际上：MuJoCo 的接触**内禀是软的**——接触力是穿透 $\phi$ 的光滑函数（虚拟弹簧-阻尼器），没有"接触/不接触"二值开关。这是它能处处可微、能 RL 大批量并行的根本。把它当硬接触理解，你永远搞不懂 solref/solimp 为什么存在。 - 为什么重要：理解"软"，才知道为什么 MuJoCo 里硬物体仍有微小穿透（特性非 bug），以及为什么它能做 diff-sim。

🧠 思维陷阱：盲目调大刚度想消除穿透 - 新手想法："物体穿进地面了，把 solref 时间常数调到极小（刚度极大）就不穿了。" - 实际上：刚度越大，系统越**刚性（stiff）**，半隐式 Euler 需要越小的步长才稳定——刚度过大 + 步长不够小 = 数值爆炸或剧烈 jitter。MuJoCo 的软接触本就用穿透换稳定，盲目调硬会失去这个好处。正确做法是平衡刚度与步长，或接受小穿透。 - 正确思维：solref/solimp 是"穿透-稳定性"的权衡旋钮。要更硬必须同时减小步长；多数任务接受 mm 级穿透即可。

💡 概念误区：把 PGS/CG/Newton 当成"解不同问题" - 新手想法："PGS、CG、Newton 是三种不同的接触模型。" - 实际上：三者求解**同一个凸优化问题**，只是数值算法不同（一阶 vs 高阶 vs 二阶）。换求解器**不改变物理模型**，只改变收敛速度与精度。Newton 默认、最精确；PGS 最快但低精度；CG 居中。 - 为什么重要：调试发现结果随求解器变，多半是某个求解器未收敛（迭代数不够），而非"物理不同"——增大迭代数或换 Newton 验证。

💡 概念误区：混淆 solref 与 solimp 的职责 - 新手想法："solref 和 solimp 都是调软硬的，随便设一个就行。" - 实际上：solref 定义参考加速度（约束"想做什么"——多快把穿透推回零，对应弹簧-阻尼器的时间常数/阻尼比）；solimp 定义阻抗如何随违反量变化（约束"能做到什么"——软到硬的过渡曲线）。前者管"目标动力学"，后者管"软硬随穿透深度的变化"。两者职责不同，需分别理解。 - 为什么重要：建模软接触层（如橡胶垫）主要调 solimp 的阻抗曲线；调整接触的"反应速度"主要调 solref 的时间常数。

练习¶

（推导题） 写出 MuJoCo 的原始 Gauss 优化问题，说明第一项 $\|x-M^{-1}(\tau-c)\|_M^2$ 如何编码 Gauss 最小约束原理。再推导其对偶问题，证明对偶 Hessian 是 $A+R$，其中 $A=JM^{-1}J^\top$。说明正则化 $R$ 为何保证严格凸（即使 $A$ 因接触冗余而奇异）。
（开放思考题） solref=(timeconst, dampratio) 把约束建模为虚拟弹簧-阻尼器。给定 timeconst $=0.02$ s、dampratio $=1$（临界阻尼），推算等效刚度/阻尼，并讨论：若步长 $h=0.002$ s，这个接触是"硬"还是"软"？timeconst 减到 $0.002$ s（等于步长）会发生什么数值现象？为什么 MuJoCo 文档建议 timeconst $\ge 2h$？
（综合判断题，桥接 §5 与本目录 40） 论证 MuJoCo 与 Anitescu 凸松弛（§5）的同构关系：(a) 把 MuJoCo 对偶问题与 Anitescu 凸 QP 并排写出，指出 $A=JM^{-1}J^\top$ 与 $N=D^\top M^{-1}D$ 对应；(b) 说明 MuJoCo 多了正则化 $R$ 和软化锥，这如何改善可微性；(c) 解释为什么这种"软化+凸化"使 MJX 能做可微仿真，而硬 LCP（Bullet）难。

§7 冲击律的离散处理 ⭐⭐¶

动机：当两个物体相撞，速度该怎么跳¶

前面几节里"恢复系数 $e$"反复出现，但我们一直没系统讲清**冲击律（impact law）** 本身。冲击律回答一个核心问题：碰撞瞬间，法向速度（乃至切向速度）按什么规则跳变？ 看似简单（中学物理 $v^+=-ev^-$），实则在多体、有摩擦、多点同时碰撞时暗藏深坑——不同冲击律给出不同能量行为，有的甚至会**人为增加能量**（违反物理）。本节梳理四种主流冲击律及其离散形式，并讲清为什么 Moreau 的选择在时步法框架下最优雅。

反面：朴素 Newton 冲击律的能量陷阱¶

最朴素的冲击律是 Newton 法则：$v_n^+=-e\,v_n^-$（法向速度按系数 $e$ 翻转）。单点碰撞没问题，但在**多点同时碰撞**或**有摩擦**时会出事：

Kane 反例（Kane's paradox）：某些多体碰撞配置下，对每个接触点独立应用 Newton 法则，会导致系统**总动能增加**——凭空多出能量，违反热力学。
原因：Newton 法则约束的是"速度比"，但在多体耦合下，各接触点的速度通过刚体相互关联，独立设定速度比会过约束或注入能量。

这促使人们发展基于**冲量**或**能量**的冲击律。

历史：四代冲击律¶

冲击律的研究横跨三个世纪：Newton（17 世纪，速度比）→ Poisson（19 世纪，冲量比，分压缩/恢复段）→ Stronge（1990s，能量比，修正 Kane 反例）→ Moreau（1980s，速度级包含，统一时步法）。它们不是简单的新旧替代，而是针对不同病态的不同修正。

理论：四种冲击律对比¶

冲击模型	定义	约束的量	能量性质
Newton	$v_n^+=-e\,v_n^-$	法向速度比	$e=1$ 保能；$e<1$ 耗散；多体下可能增能（Kane）
Poisson	$P_n^{\text{rest}}=e\,P_n^{\text{comp}}$	冲量比（分压缩段/恢复段）	摩擦多体下比 Newton 更能量一致
Moreau	$v_n^+=-e\,v_n^-$ 但作为包含 $-P\in N_{T_C}(v^++ev^-)$	速度级统一化	自动耗散摩擦，时步法原生
Stronge	能量恢复 $e_*^2=E^+/E^-$	能量比	避免 Kane 反例的能量增益

逐个深入：

Newton 冲击律。 $v_n^+=-e\,v_n^-$，$e\in[0,1]$ 恢复系数。$e=1$ 完全弹性（速度大小不变、方向翻转，动能守恒），$e=0$ 完全塑性（碰后法向速度为零，最大耗散），$0<e<1$ 部分弹性。优点：最直观、单点最准。缺点：多体同时碰撞可能增能（Kane），有摩擦时切向行为未定义。

Poisson 冲击律。 把碰撞分为**压缩段**（compression，相对速度从 $v^-$ 到零，法向冲量 $P_n^{\text{comp}}$）和**恢复段**（restitution，从零到 $v^+$，冲量 $P_n^{\text{rest}}$），约束 $P_n^{\text{rest}}=e\,P_n^{\text{comp}}$。优点：约束冲量（守恒量）而非速度，在摩擦多体下能量一致性更好。缺点：需区分两段，离散稍复杂。

值得展开 Poisson 与 Newton 在**单点无摩擦**时为何等价，而多体/有摩擦时为何分道扬镳。单点无摩擦：压缩段把法向速度从 $v_n^-$ 拉到 $0$，由动量定理 $P_n^{\text{comp}}=m_{\text{eff}}|v_n^-|$；恢复段 $P_n^{\text{rest}}=e P_n^{\text{comp}}=e\,m_{\text{eff}}|v_n^-|$，产生回弹速度 $v_n^+=P_n^{\text{rest}}/m_{\text{eff}}=e|v_n^-|$——与 Newton $v_n^+=-ev_n^-$ 完全一致。但多体/有摩擦时：压缩段的法向冲量会通过 Delassus 算子 $W$ 的非对角元（§3）耦合到其他接触和切向，使得"压缩段结束时刻"（所有接触法向速度归零的瞬间）与各点速度比解耦——此时约束**冲量比**（Poisson）和约束**速度比**（Newton）给出不同结果，Poisson 因为约束的是守恒量（冲量）而更能量一致。这就是 Poisson 在多体摩擦下优于 Newton 的根源。

Moreau 冲击律。 这是时步法的"原生"冲击律。它不单独定义碰撞，而是把恢复系数**直接嵌入 MDI 的 argument**：$-P\in N_{T_C}(v^++ev^-)$。展开看：法向分量给出 $v_n^+=-ev_n^-$（与 Newton 一致），但它作为**包含关系**自动处理摩擦（切向走最大耗散）、自动统一持续接触与冲击。优点：一套公式、自动摩擦耗散、时步法天然。这正是 §2–§3 反复强调的"冲击律嵌入 argument"。

Stronge 冲击律。 用**能量恢复系数** $e_*^2=E^+/E^-$（恢复段释放的能量与压缩段储存的能量之比）。优点：基于能量，从根本上避免 Kane 反例的能量增益（因为它直接约束能量比 $\le 1$）。缺点：物理量较抽象，工程使用较少。

Moreau 选择的深刻意义¶

为什么时步法（Moreau-Jean、Siconos）几乎都用 Moreau 冲击律？

在 sweeping-process 框架下，恢复系数 $e$ 直接出现在 "$v^++ev^-$" 这个 argument 里。这带来一个决定性好处：时步格式无需区分"碰撞"与"持续接触"——同一套公式处理二者。

持续接触时（$v^-\approx 0$）：argument $v^++ev^-\approx v^+$，约束退化为 $v^+\ge 0$（不侵入），接触力是普通的持续接触力。
碰撞时（$v^-<0$ 显著）：argument $v^++ev^-$ 强制 $v^+=-ev^-$，冲量 $P$ 取 atom 值。

本质洞察：Moreau 把恢复系数 $e$ 从"碰撞瞬间的特殊规则"变成"MDI argument 的一个常驻参数"。这是 time-stepping 全部魅力的浓缩——冲击不是额外的"事件"，而是冲量测度 $P$ 的 atomic 分量；恢复系数不是碰撞专用开关，而是统一约束里的一个系数。这就是为什么本章反复说"冲击律是时步的自然产物而非外挂"。事件驱动法必须把冲击当独立子问题求解，时步法让它自然涌现。

多体同时冲击的病态¶

最后必须提一个理论与工程都头疼的难点：多体同时冲击。

三球同时碰撞（Newton's cradle 的退化情形）、方块落桌四角同触：这些配置下，无唯一 Newton 解——速度跳变不被独立的逐点 $e$ 唯一确定，依赖碰撞的"微观顺序"（哪个点先碰），而理想刚体里这个顺序是奇异的。
解决方案：(1) Frémond 多体冲击律（Acary 等 2024）——用一个全局能量泛函约束所有接触；(2) energy-consistent Pfeiffer-Glocker 投影——把冲击后速度投影到能量不增的可行集。
工程现实：多数仿真器对此"睁一只眼"——给出某个解（依赖求解器枢轴顺序），不保证唯一。这是接触仿真"换求解器结果变"的又一来源（呼应 §4 的 LCP 非唯一性）。

阶段小结：到这里我们梳理了四种冲击律（Newton 速度比、Poisson 冲量比、Stronge 能量比、Moreau 速度级包含）及其能量行为，理解了 Moreau 冲击律为何在时步法中最优雅（恢复系数嵌入 argument、统一持续接触与冲击），以及多体同时冲击的非唯一病态。接下来 §8 回到积分器层面，讲半隐式 vs 全隐式，以及辛积分为何与接触不相容。

⚠️ 常见陷阱¶

💡 概念误区：以为恢复系数 $e$ 是材料的固有常数 - 新手想法："钢球的 $e$ 就是 0.9，查表就行。" - 实际上：$e$ 不是纯材料属性——它依赖碰撞速度、几何、温度，且 Newton/Poisson/Stronge 三种定义下的"$e$"数值不同（速度比 vs 冲量比 vs 能量比）。高速碰撞 $e$ 通常更小（更多塑性变形耗能）。把 $e$ 当固定常数会在跨速度域仿真时失准。 - 为什么重要：标定 $e$ 时必须注明用的是哪种冲击律、什么速度范围；跨域时 $e$ 应随速度变化。

🧠 思维陷阱：认为 $e=1$（完全弹性）一定能量守恒 - 新手想法："$e=1$ 是完全弹性碰撞，能量必然守恒。" - 实际上：单点碰撞 $e=1$ 确实保能，但**多体同时碰撞 + Newton 法则**下，即使每点 $e=1$，总能量也可能**改变**（Kane 反例）——因为逐点速度比约束在多体耦合下不等价于能量守恒。要保证多体保能需用 Stronge（能量比）或能量一致投影。 - 正确思维：$e=1$ 保能仅对单点成立；多体保能需要能量级的冲击律或投影。

💡 概念误区：把切向恢复也用一个系数描述 - 新手想法："法向有 $e$，切向也设个切向恢复系数不就行了？" - 实际上：切向行为由**摩擦**（Coulomb 锥 + 最大耗散）决定，不是简单的"切向恢复系数"。碰撞时切向冲量受法向冲量约束（$\|\beta\|\le\mu P_n$），且方向由滑动方向决定。强行引入切向恢复系数会与摩擦锥冲突、可能增能。 - 为什么重要：法向用恢复系数、切向用摩擦锥——这是两套不同机制，不能混用一个系数。

练习¶

（推导题，草稿纸完成） 对两个质量分别为 $m_1,m_2$ 的小球一维对心碰撞，初速 $v_1^-,v_2^-$，恢复系数 $e$。用 Newton 冲击律 + 动量守恒，解出碰后速度 $v_1^+,v_2^+$。计算动能变化 $\Delta E$，证明 $\Delta E=-\tfrac12(1-e^2)\mu_{12}(v_1^--v_2^-)^2$（$\mu_{12}=m_1m_2/(m_1+m_2)$ 约化质量），从而 $e<1$ 时 $\Delta E<0$（耗散），$e=1$ 时守恒。
（开放思考题） 构造或描述一个 Kane 反例：找一个多点同时碰撞配置（如一个杆同时撞两个固定点），说明逐点应用 Newton 法则（每点 $e=1$）为何可能导致总动能增加。再说明 Stronge 能量比冲击律如何避免这个问题。提示：考虑各接触点速度通过刚体的耦合。
（综合题） 解释 Moreau 冲击律 $-P\in N_{T_C}(v^++ev^-)$ 如何"统一"持续接触与冲击：分别取 $v^-\approx 0$（持续接触）和 $v^-<0$（碰撞）两种情形，展开这个包含关系，说明它分别退化为什么条件。这与事件驱动法"先积分再 apply $v^+=-ev^-$"的本质区别是什么？

§8 半隐式 vs 全隐式积分、辛性与接触的不相容 ⭐⭐⭐¶

动机：积分器的"隐式程度"如何选¶

§3 的 θ-scheme 已经触及"显式-隐式"光谱（$\theta=0$ 显式、$\theta=1$ 隐式）。本节系统讲清积分器的三个层次——显式 / 半隐式 / 全隐式，它们在稳定性、能量、计算量上的权衡，以及一个深刻的负面结果：辛积分器（symplectic integrator）与接触根本不相容。后者解释了一个反直觉的事实——为什么 MuJoCo 不用"更高级"的 Verlet/变分积分器，而用"朴素"的半隐式 Euler。

反面：显式 Euler 为什么在接触下不堪用¶

最朴素的显式 Euler：$v_{\ell+1}=v_\ell+hM^{-1}F_\ell$，$q_{\ell+1}=q_\ell+hv_\ell$（位置用**旧**速度）。它在接触下的灾难：

能量注入：显式 Euler 对振荡系统**人为增加能量**（不稳定）。接触相当于极硬的弹簧，显式 Euler 会让接触振荡发散。
刚性限制：硬接触刚度 $k$ 大，显式稳定性要求 $h<2/\sqrt{k/m}$——刚度越大步长越小，硬接触下步长趋于零。
位置滞后：用旧速度更新位置，接触力来不及阻止穿透。

这就是为什么**没有任何严肃仿真器用纯显式 Euler 处理接触**。

理论：三个层次的积分器¶

层次一：显式（explicit）。 $v_{\ell+1}=v_\ell+hM^{-1}F(q_\ell,v_\ell)$，$q_{\ell+1}=q_\ell+hv_\ell$。右端全用旧状态。优点：无需解方程，最快。缺点：能量注入、刚性差、接触下发散。

层次二：半隐式 / symplectic Euler（semi-implicit）。 先更新速度，再用**新速度**更新位置：

\[ v_{\ell+1}=v_\ell+hM^{-1}F(q_\ell,v_\ell)+hM^{-1}J^\top\lambda,\qquad q_{\ell+1}=q_\ell+h\,v_{\ell+1}. \]

关键是 $q$ 用 $v_{\ell+1}$（新速度）。优点：对光滑 Hamiltonian 系统是**辛**的（symplectic，能量长期有界不漂移），稳定性远好于显式，仍然便宜（接触力可半隐式解）。这是 MuJoCo、Bullet、Drake 的主力积分器。

层次三：全隐式（fully implicit）。 右端全用新状态：$v_{\ell+1}=v_\ell+hM^{-1}F(q_{\ell+1},v_{\ell+1})+\dots$，需要解非线性方程组（Newton 迭代）。优点：最稳定（A-稳定、L-稳定），可用大步长，对刚性系统友好。缺点：每步解非线性系统，最贵；引入数值耗散（$\theta=1$ 的 $O(h^2)$ 耗散）。Drake SAP 的 midpoint、implicitfast 倾向这一端。

层次	位置更新用	能量行为	刚性	计算量	代表
显式 Euler	旧速度 $v_\ell$	注入（发散）	差	最低	几乎无人用于接触
半隐式 Euler	新速度 $v_{\ell+1}$	辛（长期有界）	中	低	MuJoCo、Bullet、Drake
全隐式 / θ=1	新状态	耗散 $O(h^2)$	好	高（Newton）	Drake SAP midpoint

本质洞察：半隐式 Euler 是"性价比之王"——它只比显式 Euler 多改一处（位置用新速度），却获得辛性（光滑段能量长期不漂移）和远好的稳定性。这就是为什么它统治机器人仿真：在"够稳 + 够快 + 够准"的三角里，它是接触仿真的甜点。全隐式更稳但太贵，显式太不稳。

稳定性界算例：刚度如何限制步长¶

为了把"刚性限制步长"从口号落到数字，算一个软接触（弹簧刚度 $k$）的稳定性界——这直接解释了 §9 的 jitter 和 §6 的"timeconst $\ge 2h$"。

把一个接触建模为刚度 $k$、阻尼 $c$ 的弹簧-阻尼器（MuJoCo 风格软接触），质量 $m$。接触动力学是 $m\ddot\phi=-k\phi-c\dot\phi$，固有频率 $\omega_0=\sqrt{k/m}$。

显式 Euler 稳定界：对无阻尼弹簧 $\ddot\phi=-\omega_0^2\phi$，显式 Euler 的放大因子模 $|G|=\sqrt{1+\omega_0^2h^2}>1$——恒不稳定（任意步长都发散）。这是显式 Euler 处理接触的死刑判决。

半隐式 Euler 稳定界：半隐式 Euler 对无阻尼弹簧的放大因子满足 $|G|=1$ 当且仅当 $\omega_0 h\le 2$，即

\[ h\le\frac{2}{\omega_0}=2\sqrt{\frac{m}{k}}. \]

超过这个界，半隐式 Euler 也发散。物理意义：步长不能超过接触固有周期的 $1/\pi$ 量级——接触越硬（$k$ 大、$\omega_0$ 大），允许的步长越小。这正是 §6 "timeconst $\ge 2h$" 的来历：timeconst $=1/\omega_0$，约束 $h\le 2\cdot$timeconst$\Leftrightarrow$ timeconst $\ge h/2$（MuJoCo 取更保守的 $\ge 2h$ 留余量）。

数字示例：钢接触刚度 $k=10^6$ N/m，$m=1$ kg，则 $\omega_0=1000$ rad/s，半隐式 Euler 稳定界 $h\le 2/1000=2\times 10^{-3}$ s。若要更硬的 $k=10^8$，则 $h\le 2\times 10^{-4}$ s——刚度增 100 倍，步长须减 10 倍。这就是硬接触的代价：要么用极小步长（慢），要么软化接触（MuJoCo 的选择，§6），要么全隐式（无条件稳定但每步贵）。

全隐式（θ=1）稳定界：后向 Euler 对弹簧**无条件稳定**（$|G|<1$ 对所有 $h$）——这是它的卖点。代价是 $O(h^2)$ 耗散：大步长时把接触振荡过度抹平。

本质洞察：积分器的"刚性"本质是**稳定步长 vs 接触刚度**的权衡。显式：恒不稳（出局）。半隐式：$h\le 2/\omega_0$，硬接触逼小步长。全隐式：无条件稳但耗散 + 贵。MuJoCo 用半隐式 + 软化（限制 $\omega_0$ 让 $h$ 可大），Drake SAP 用半隐式/全隐式 + 凸求解（大步长稳定），各有取舍。这个稳定界 $h\le 2/\omega_0$ 是连接 §6（timeconst）、§8（积分器）、§9（jitter 来自步长超界）的数学纽带——记住它，三节贯通。

核心负面结果：辛积分器与接触不相容¶

现在到本节最深刻、也最反直觉的结论。在**光滑** Hamiltonian 力学里，辛积分器（Verlet、变分积分器 Marsden-West）是"圣杯"——它们保辛、近乎能量守恒、保动量，长时间轨迹（天体力学、分子动力学）极准。一个自然的想法是：接触仿真也用辛积分器不是更好吗？

答案是**不行**，且有严格证明。

Cirak-West 2005（IJNME 64:1079）定理：显式辛格式在**接触激活的瞬间必然破坏辛性（symplecticity）**。

直觉解释：辛性要求离散流是**辛映射**（保相空间体积/辛形式 $\omega$）。但接触约束的激活/失活是不等式约束激活集的**突变**——它不是微分同胚（diffeomorphism），在激活瞬间相空间被"折叠"（穿透被修正、速度被投影），这个折叠**不保辛**。穿透修正（projection）尤其致命：Cirak-West 证明任何修正穿透的操作都破坏辛结构。

进一步：Kane-Marsden-Ortiz-West 2000（IJNME 49:1295）证明，变分积分器在**弹性接触**（$e=1$）下**可以**保辛——但仅限完全弹性。一旦 $e<1$（塑性、耗散），辛性必然丢失。

结论： $$ \text{完全塑性冲击} (e=0) + \text{辛} = \textbf{矛盾};\qquad \text{保辛仅在} e=1 \text{纯弹性} + \text{无穿透修正时成立}. $$

而机器人最常见的恰恰是 $e\approx 0$ 的塑性接触（脚落地不反弹、手抓物不弹开）——辛积分器在这里完全失效。

为什么 MuJoCo 不用变分积分器¶

这就解释了 §0 埋下的伏笔。MuJoCo 面对一个选择：

变分积分器（保辛）：长时间能量守恒好，但与塑性接触（$e\approx 0$）、穿透修正**根本不相容**（Cirak-West 2005）。
半隐式 Euler + soft-convex solver：不保辛，但**可微 + 时间可逆**，对塑性接触无障碍。

MuJoCo 选了后者。理由：对 MPC/RL，可微性和时间可逆性远比保辛重要。机器人任务里接触几乎处处是塑性的，保辛带来的"长期能量守恒"优势在这里根本无从发挥（接触本就该耗能），而它要付出的"与接触不相容"代价却是致命的。

本质洞察：这是一个"工具的优越性依赖场景"的绝佳案例。辛积分器在光滑天体力学里是无可争议的王者，但搬到接触机器人里却**水土不服**——因为它的核心优势（保辛、长期能量守恒）建立在"系统光滑且能量守恒"的前提上，而接触恰恰打破了这两个前提（非光滑、耗能）。Cirak-West 的负面结果不是"辛积分器不够好"，而是"它的优势前提在接触下不成立"。这也反向印证了本章主线：接触动力学是一个**独立**的数学领域，光滑力学的明珠（辛积分）在这里黯然失色，而"朴素"的半隐式 Euler + 凸求解器反而是正解。

桥接：与 §3 θ-scheme 的统一¶

§3 的 θ-scheme 其实就是这个"显式-隐式"光谱的统一参数化：$\theta=0$ 是显式、$\theta=1$ 是全隐式、$\theta=\tfrac12$ 是梯形（接近辛的能量中性）。本节的"半隐式 Euler"对应一种特定的算子分裂（速度隐式、位置半隐式）。Siconos 的 MoreauJeanOSI（θ-scheme）、Drake 的 SI-Euler、MuJoCo 的半隐式 Euler，都是这个光谱上的点。理解光谱，你就能在任意仿真器里定位它的积分器。

阶段小结：到这里我们厘清了积分器的三层次（显式/半隐式/全隐式）及其能量-稳定性-计算量权衡，确立了半隐式 Euler 为机器人仿真主力的原因，并给出深刻负面结果——辛积分器与塑性接触/穿透修正不相容（Cirak-West 2005），解释了 MuJoCo 为何弃辛取半隐式。最后一节 §9 处理时步法的"数值病理"：穿透、能量漂移、约束漂移及其稳定化。

⚠️ 常见陷阱¶

💡 概念误区：以为半隐式 Euler 和显式 Euler 差不多 - 新手想法："都是一阶 Euler，半隐式不就是换了下顺序，能差多少？" - 实际上：差别是**质变**。半隐式 Euler 用新速度更新位置，对光滑 Hamiltonian 系统**保辛**（能量长期有界），而显式 Euler 注入能量（发散）。一个简单的弹簧振子，显式 Euler 振幅指数增长，半隐式 Euler 振幅长期稳定。仅仅"换个顺序"带来辛性 vs 发散的天壤之别。 - 为什么重要：这解释了为什么所有仿真器都用半隐式而非显式——不是细节优化，是稳定性的根本差异。

🧠 思维陷阱：认为"保辛 = 更好的接触积分器" - 新手想法："变分积分器保辛、保动量、近能量守恒，肯定比半隐式 Euler 适合接触。" - 实际上：Cirak-West 2005 证明辛性与接触激活/穿透修正**根本不相容**（除非 $e=1$ 纯弹性无修正）。机器人接触多是 $e\approx 0$ 塑性，辛积分器在此**完全失效**。保辛的优势（长期能量守恒）在"本就该耗能"的接触场景毫无用武之地。 - 正确思维：积分器的优劣依赖场景。光滑长轨迹用辛；接触耗能用半隐式 + 凸求解。不要把光滑力学的结论盲目搬到接触。

💡 概念误区：以为全隐式总是最好（最稳定） - 新手想法："全隐式 A-稳定，那就总用全隐式，步长可以很大。" - 实际上：全隐式每步要解非线性系统（Newton 迭代），最贵；且引入 $O(h^2)$ 数值耗散，大步长时会过度抹平动力学（能量被人为抽走，运动变"黏稠"）。对需要保留动力学细节或长期能量守恒的任务，全隐式的耗散是缺点。 - 正确思维：全隐式适合刚性接触防 jitter、大步长 MPC；不适合需要能量守恒/动力学保真的场景。稳定性与精度/计算量是权衡。

练习¶

（推导题，草稿纸完成） 对一维无阻尼弹簧振子 $\ddot x=-\omega^2 x$，分别写出显式 Euler 和半隐式 Euler 的一步更新。计算两者的"数值能量" $E_\ell=\tfrac12(v_\ell^2+\omega^2 x_\ell^2)$ 的一步变化 $E_{\ell+1}/E_\ell$。证明显式 Euler 的能量按 $(1+\omega^2h^2)$ 增长（发散），半隐式 Euler 的能量被一个守恒量约束（长期有界）。这就是辛性的体现。
（开放思考题） 用直觉语言解释 Cirak-West 2005 的结论：为什么"修正穿透"这个操作会破坏辛性？提示：辛映射要保相空间体积/辛形式，而把穿透的物体"推回表面"是一个把相空间局部"压扁/折叠"的操作。再说明为什么 $e=1$ 纯弹性碰撞（无穿透修正、能量守恒）可以例外。
（综合判断题） 给定四个任务：(a) 仿真太阳系行星 100 万年轨道；(b) 仿真四足机器人小跑 10 秒；(c) 仿真一个无摩擦弹性球在盒子里弹跳（$e=1$）100 万次；(d) RL 训练机械臂抓取。对每个任务判断应该用辛/变分积分器还是半隐式 Euler + 凸求解器，并说明判据（光滑性、$e$ 值、可微需求、能量守恒需求）。

§9 穿透、能量与约束漂移：时步法的数值病理与稳定化 ⭐⭐¶

动机：仿真"看起来不对"时，问题出在哪¶

到此你已掌握四种格式的数学。但实际跑仿真时，你会遇到一堆"看起来不对"的现象：物体慢慢沉进地面（穿透/漂移）、关节莫名其妙抖动（jitter）、能量莫名增加或减少、堆叠的方块缓慢"融化"。这些不是 bug，而是时步法的**数值病理（numerical pathology）**——离散化的固有副作用。本节系统梳理这些病理及其稳定化技术，是全章"调参排错"的总纲。

反面：为什么不稳定化会出问题¶

时步法的约束是**离散满足**的——速度级约束 $\nabla\phi^\top v_{\ell+1}\ge 0$ 只保证这一步速度合法，但**位置级约束 $\phi\ge 0$ 会随离散误差累积漂移**。具体：

速度级格式（Moreau-Jean）：每步只管速度不侵入，但数值误差让位置缓慢偏离约束流形——物体逐步陷入地面（约束漂移 / drift）。
即使位置级格式（Stewart-Trinkle），线性化 $\phi_\ell+h\nabla\phi^\top v_{\ell+1}\ge 0$ 也只在线性化精度内成立，仍有 $O(h^2)$ 残余穿透。

不加稳定化，长时间仿真这些误差累积成可见的物理失真。

理论：三类病理与对应稳定化¶

病理一：约束漂移（constraint drift）¶

现象：位置级约束 $\phi(q)=0$（或等式约束 $g(q)=0$）随时间缓慢被违反，物体沉入接触面或关节脱开。

稳定化技术：

(1) Baumgarte 稳定化（1972）：把约束违反当作一个受控的二阶系统，引入人工阻尼：

\[ \ddot\phi + 2\alpha\dot\phi + \beta^2\phi = 0. \]

这强制约束违反 $\phi$ 按阻尼振荡衰减回零。优点：实现简单（加到约束方程右端）。致命缺点：$(\alpha,\beta)$ 极难调——取太小漂移修正慢，取太大给系统**注入人工能量或过阻尼**。Bullet 的默认 ERP（error reduction parameter）$=0.2$ 本质是 Baumgarte，它会把机器人关节振荡"人为压平"，导致 sim-to-real gap（仿真里关节比真实更"软"）。

(2) GGL 稳定化（Gear-Gupta-Leimkuhler）：把位置约束和速度约束**同时**作为约束，投影到约束流形。比 Baumgarte 更稳健（不引入可调阻尼），但需解更大的系统。

(3) 坐标投影（coordinate projection）：每步求解后，把位置（和速度）投影**回约束流形 $\{\phi=0\}$。 - **post-stabilization：每步后只投影位置。 - coordinate projection：同步投影位置 + 速度。 - Siconos CombinedProjectedMoreauJeanOSI、Drake 内置 coordinate projection 都属此类。比 Baumgarte 更推荐——不引入人工能量，几何上干净。

稳定化	原理	优点	缺点
Baumgarte	人工阻尼 $\ddot\phi+2\alpha\dot\phi+\beta^2\phi=0$	简单	参数难调、注入能量、sim2real gap
GGL	位置+速度同时约束	稳健、无可调阻尼	系统更大
坐标投影	每步投影回流形	几何干净、不注入能量	投影开销

病理二：穿透（penetration）¶

现象：硬物体之间出现可见的相互嵌入。

来源与对策： - 速度级格式：固有漂移 → 用坐标投影修正。 - 软接触（MuJoCo）：穿透是**设计特性**（软度换稳定），由 solref/solimp 控制深度；减小穿透需调硬 + 减小步长（§6 陷阱）。 - 位置级格式（Stewart-Trinkle）：线性化残余 $O(h^2)$ → 减小步长或用几何隐式时步法（geometrically implicit，CMU 的 ijrr_dynsim 工作改进了 gap 函数的隐式处理）。

病理三：能量漂移与 jitter¶

现象：(a) 能量缓慢增加（不稳定）或减少（过度耗散）；(b) 静止堆叠的物体高频抖动（jitter）。

来源与对策： - 能量增加：通常是显式分量或 Baumgarte 注入 → 改半隐式/全隐式、弃 Baumgarte 用投影。 - 能量过度减少：$\theta=1$ 或大数值阻尼 → 改 $\theta=\tfrac12$、减小阻尼。 - jitter：接触力在"压紧/松开"间高频切换，常因步长过大、刚度过高、或求解器未收敛 → 增大求解器迭代数、适当增大数值阻尼（$\theta\to 1$）、减小步长、或用 warm-start 稳定接触力。

自适应步长与高阶恢复（档位 4 预告）¶

时步法天然是固定低阶步长（§1）。要在保持鲁棒的同时恢复精度，前沿工作： - Chen-Acary-Brüls 2013（IJNME 96:487）非光滑 generalized-α：数值阻尼可调，与 Hilber-Hughes-Taylor 兼容，给柔性多体 + 接触。 - Acary 2013（CMAME 256:224）projected event-capturing Runge-Kutta：自适应步长 + 位置投影，在光滑段重获高阶精度。 - Brüls-Acary-Cardona 2014：Lie-group generalized-α 用于柔性多体接触。

这些属于章末"进阶章节"，这里只需知道：固定低阶不是时步法的终点，加投影 + 自适应可部分恢复高阶。

本质洞察：时步法的所有"数值病理"（漂移、穿透、jitter、能量漂移）本质上都是**离散化把连续约束打折扣**的后果——连续时间严格满足的约束，离散后只能近似满足。稳定化技术（Baumgarte/GGL/投影）就是各种"把约束拉回来"的手段，它们的优劣在于**拉回来的同时是否扰动了物理**：Baumgarte 拉回来但注入能量（脏），坐标投影拉回来且不注入能量（干净）。理解这个，调参排错就从"试错"变成"对症"。

阶段小结：到这里我们梳理了时步法的三类数值病理——约束漂移、穿透、能量漂移/jitter——及其稳定化（Baumgarte 简单但脏、GGL/坐标投影干净），并预告了自适应步长恢复高阶的前沿。至此九节理论 + 工程主体完成。下面进入本章的认知巩固部分：误解汇总、符号表、定理速查表、与后续章节的桥接、故障排查手册。

⚠️ 常见陷阱¶

💡 概念误区：把所有穿透都当成 bug - 新手想法："物体穿进去了一点点，肯定是代码错了。" - 实际上：穿透分两类——(1) 软接触（MuJoCo）的穿透是**设计特性**（软度换稳定，受 solref/solimp 控制）；(2) 速度级格式的漂移是**离散副作用**（需投影修正）。只有第二类需要"修"，第一类是预期行为。先判断仿真器类型再决定是否处理。 - 为什么重要：在 MuJoCo 里追求零穿透（疯狂调硬）会牺牲稳定性，南辕北辙。

🧠 思维陷阱：用 Baumgarte 稳定化解决一切漂移 - 新手想法："漂移了就加 Baumgarte，调大 $\beta$ 直到不漂。" - 实际上：Baumgarte 注入人工能量/阻尼，参数难调，且会引入 sim-to-real gap（Bullet 默认 ERP=0.2 把关节振荡人为压平）。调大 $\beta$ 可能让系统过阻尼或不稳定。更推荐坐标投影/GGL——几何干净、不注入能量。 - 正确思维：漂移优先用坐标投影；Baumgarte 仅在投影不可行时作为简单后备，且需谨慎调参。

💡 概念误区：jitter 是物理现象 - 新手想法："堆叠的方块在抖，可能是真实的微振动。" - 实际上：静止堆叠的 jitter 几乎总是**数值伪影**——接触力在"压紧/松开"间高频切换，源于步长过大、刚度过高或求解器未收敛。真实静止物体不会这样抖。对策：增大求解器迭代数、加数值阻尼、减小步长、warm-start。 - 为什么重要：把 jitter 当物理会让你忽略真正的数值问题；识别它是伪影才能对症（多半是求解器没收敛）。

练习¶

（开放思考题） 你的仿真里方块堆叠后缓慢"融化"下沉。诊断这是约束漂移、软接触穿透还是能量问题？分别针对速度级格式（Moreau-Jean）、软接触（MuJoCo）、位置级格式（Stewart-Trinkle）给出最可能的原因和对应的稳定化方案。
（推导题） Baumgarte 稳定化 $\ddot\phi+2\alpha\dot\phi+\beta^2\phi=0$ 的解是阻尼振荡。求其衰减时间常数与 $(\alpha,\beta)$ 的关系，并说明：若 $\alpha$ 取得使系统欠阻尼（$\alpha<\beta$），约束违反会如何演化？若过阻尼（$\alpha>\beta$）呢？为什么"临界阻尼"$\alpha=\beta$ 在实践中难以精确命中、从而 Baumgarte 难调？
（综合题，桥接 §3/§6/§8） 一个机器人手抓取仿真出现关节 jitter。请按以下顺序排查：(a) 检查积分器（§8）是否半隐式以上；(b) 检查 θ 或数值阻尼（§3）是否足够；(c) 检查求解器（§6）迭代数是否收敛；(d) 检查步长与刚度匹配。说明每一步你会观察什么指标、如何判断是否是该环节的问题。

§10 主流仿真器的时间步进实现对比 ⭐⭐¶

动机：把数学路线映射到真实代码¶

前面九节建立了全部数学。本节是"数学 → 代码"的字典——打开任意仿真器源码时，你能指出它在哪条数学路线上、冲击律如何嵌入、步长稳定性来自何处。这是把理论变成工程能力的最后一公里。

理论：六大仿真器全景对比¶

仿真器	时步法	求解器	问题类型	冲击处理	典型 $h$	优点	缺点
Siconos	Moreau-Jean θ, Schatzman-Paoli, event-driven, NS-gen-α	NSGS, AC/FB/NM 非光滑 Newton, Proximal, VI-ExtraGradient, IPM	LCP, VI, ConvexQP, SOCCP, FrictionContact	Newton/Moreau/Frémond 原生	$10^{-5}$–$10^{-3}$	数学最严谨；求解器最全	机器人 API 弱；非实时
Drake TAMSI	半隐式 Euler	非凸 NCP + 切向角限制 line-search	非凸 NCP	柔性 Hunt-Crossley	$10^{-3}$	物理连续；成熟	非凸可能失败
Drake SAP	SI-Euler / midpoint(二阶)	Newton + exact line-search 凸势函数	无约束凸优化	柔性 Kelvin-Voigt；刚性极限	$10^{-3}$–$10^{-2}$	全局收敛；warm-start	凸松弛 $O(\delta t)$ gliding
MuJoCo	SI-Euler / implicitfast / RK4	PGS / CG / Newton（同一 soft-convex 问题）	凸 QP/SOCP + soft complementarity	无显式；`solref`/`solimp` 软化	$10^{-3}$–$10^{-2}$	可微 + JAX 百万并行	刚性保真不足
Chrono NSC	Anitescu 半隐式；HHT	APGD(Nesterov), PSOR, BB, ADMM	凸 CCP (SOCCP)	Anitescu 速度冲击	$10^{-3}$	granular $>10^6$；GPU 扩展	小规模收敛慢于 SAP
Bullet	SI-Euler + Sequential Impulse	SI ≡ PGS, MLCP (Lemke/Dantzig), Featherstone	近似 NCP	Newton $e$ via bias velocity；Baumgarte ERP	$10^{-3}$–$1/60$	简单、快、warm-start	线性收敛；高刚度 jitter
Pinocchio 3 / Simple	SI-Euler + per-step NCP	ADMM（自动 ρ）+ proximal, PGS, CCP	NCP / SOCCP（统一刚性/柔性）	Proximal ≡ Moreau-Yosida	$10^{-4}$–$5\times10^{-3}$	数学严谨；解析可微	年轻；工业经验少

四条技术路线的本质归类¶

把六个仿真器归入四条技术路线，对应本章的四种数学格式：

Siconos = 直面 NCP/SOCCP 学派（§2–§4）——不回避非光滑，直接解互补/变分不等式，求解器最全、数学最严谨。
Drake-SAP / Pinocchio-ADMM = 凸 primal 学派（§5，2020s 主流）——Anitescu 凸松弛的工业落地，全局收敛 + warm-start + 可微。
MuJoCo = 连续时间 + soft-convex（§6）——软化 + 凸化，换可微性与百万并行，RL 的事实标准。
Chrono-APGD / Bullet-SI = 迭代 per-contact 投影派——逐接触投影迭代（PGS/APGD），granular 与游戏引擎的主力。

核心源码速查¶

打开源码时直接定位这些文件：

Siconos：numerics/src/FrictionContact/fc3d_nsgs.c、fc3d_nonsmooth_Newton_AlartCurnier.c、gfc3d_ipm.c；kernel/src/simulationTools/MoreauJeanOSI.cpp
Drake：multibody/plant/tamsi_solver.cc；multibody/contact_solvers/sap/sap_solver.cc
MuJoCo：src/engine/engine_solver.c（mj_solPGS/mj_solCG/mj_solNewton）、src/engine/engine_core_constraint.c；MJX：mjx/_src/solver.py
Chrono：src/chrono/solver/ChSolverAPGD.cpp、ChSolverADMM.cpp、ChSolverPSOR.cpp
Bullet：src/BulletDynamics/ConstraintSolver/btSequentialImpulseConstraintSolver.cpp
Pinocchio：include/pinocchio/algorithm/admm-solver.hxx、contact-cholesky.hpp

机器人概念预告表：数学对象 → 机器人实例 → 仿真器¶

数学对象	机器人实例	对应仿真器
Moreau-Jean θ-step	刚体仿真、granular 沙地	Siconos, Pinocchio 3, RaiSim
Stewart-Trinkle LCP	硬约束刚体、塔积木	早期 Drake、Bullet MLCP、ODE
Anitescu 凸 CCP	滚动接触、大规模多体	Chrono APGD, ProjectChrono GPU
TAMSI	柔性接触抓取（regularized Coulomb）	Drake（legacy manipulation）
SAP (convex primal)	manipulation、peg-in-hole、大步长 MPC	Drake（当前默认）、ContactBench 比较基线
ADMM + Proximal	humanoid 多接触、ill-conditioned Delassus	Pinocchio 3 / Simple
soft-convex (semi-implicit Euler)	RL 训练环境、可微 rollout	MuJoCo, MJX (JAX)
Sequential Impulse (PGS)	游戏引擎、教学 demo	Bullet, PhysX, Havok, ODE

本质洞察：六个仿真器看起来五花八门，但**全部建立在 Delassus 算子 $W=JM^{-1}J^\top$ 上**——区别只在四件事：(1) 约束写位置级还是速度级（§2）；(2) 摩擦锥用多面体还是二阶锥（§4 vs §5）；(3) 是否凸松弛、是否软化（§5–§6）；(4) 用什么数值算法解（Lemke/PGS/APGD/ADMM/Newton）。掌握这四个维度，任何新仿真器你都能在一分钟内归类。这就是本章"四条主线"的终极价值——它不是四个孤立工具，而是一张能给任意仿真器定位的坐标系。

⚠️ 常见陷阱¶

💡 概念误区：以为"更新的仿真器一定更好" - 新手想法："Pinocchio 3 / Drake SAP 是最新的，肯定比 Bullet 强，全用新的。" - 实际上：每个仿真器有其甜点。Bullet 简单快 + warm-start，游戏/快速原型最佳；MuJoCo 可微 + 并行，RL 之王；Drake SAP 大步长 + 全局收敛，MPC/操作首选；Siconos 数学严谨，研究基准；Chrono granular 之王。"最好"取决于任务的接触刚度、接触数、可微需求、实时约束。 - 为什么重要：选错仿真器（如用 Siconos 做 RL 百万并行，或用 Bullet 做高精度科研基准）会事倍功半。

🧠 思维陷阱：以为换仿真器结果应该完全一致 - 新手想法："物理定律一样，不同仿真器跑同一场景结果应该相同。" - 实际上：不同仿真器用不同离散格式、不同摩擦锥、不同凸松弛、不同求解器，结果**系统性不同**——尤其在多体同时冲击（解非唯一，§4/§7）、大步长（凸松弛 gliding，§5）、高刚度（jitter，§9）时。SimBenchmark（Hwangbo 2018）、ContactBench（Le Lidec 2024）专门量化这些差异。 - 正确思维：跨仿真器结果差异是常态而非 bug。要可复现的科研结果，必须固定仿真器 + 版本 + 求解器 + 步长，并报告这些。

练习¶

（综合判断题） 给定四个任务：(a) $10^6$ 颗粒沙堆坍塌；(b) 机械臂 peg-in-hole 装配，需大步长 MPC；(c) RL 训练人形机器人行走，需百万并行；(d) 科研论文需可复现的高精度刚体碰撞基准。为每个任务选一个仿真器并说明理由（对照仿真器对比表的优缺点）。
（源码定位题） 从核心源码速查表任选一个仿真器，下载其源码，定位它的接触求解器主函数，对照本章 §3–§6 判断：它在哪条数学路线？位置级还是速度级？摩擦锥用什么？验证你的判断与对比表一致，并指出源码里哪一行对应 Delassus 算子 $W=JM^{-1}J^\top$ 的构造。
（跨章综合题，需综合 §4 §5 §6 §10） 给定同一个"立方体堆叠"场景，预测它在 Bullet（SI/PGS）、Drake SAP（凸 primal）、MuJoCo（soft-convex）三个仿真器里的不同行为：哪个穿透最小？哪个最可能 jitter？哪个有 $O(h)$ gliding？哪个能拿到梯度？为每个预测给出数学依据（引用对应小节）。

§11 进阶：收敛性的严格证明与扩展（档位 4）⭐⭐⭐⭐¶

定位：本节属档位 4（博士毕业+）。前面九节给出了各格式的"结论"（收敛、可解、能量行为），本节补上**为什么**——收敛证明的完整链条、辛不相容的深层机制、以及向无穷维（柔体/有限元）的扩展。研究型读者必读，工程读者可跳过。

动机：从"能算"到"算对"的严格保证¶

§3–§5 反复说"$h\to 0$ 收敛到 MDI 解"。但**收敛到 MDI 解**这件事，凭什么成立？它不是显然的——离散解是分段的、跳跃的，极限为什么恰好是连续时间的"正确答案"？本节给出 Stewart 1998 的证明骨架，它是非光滑数值分析的**模板**：任何人想为新时步格式证收敛，都走这四步。

理论：Stewart 1998 收敛证明的四步链条¶

目标：证明 Stewart-Trinkle 离散解 $(q_h,v_h)$ 在 $h\to 0$ 时（沿子列）收敛到 Moreau MDI 的解，且冲量测度 $P_h\rightharpoonup P$。

Step 1：离散 BV 估计（一致变差有界）。 第一步要证：所有步长 $h$ 的离散速度 $v_h$ 的**总变差**（total variation）一致有界：

\[ \mathrm{TV}(v_h)=\sum_\ell\|v_{\ell+1}-v_\ell\|\le C,\qquad C\ \text{与}\ h\ \text{无关}. \]

直觉：速度的总跳跃量被能量约束（动量守恒 + 接触耗散）控制。冲量不会无界，否则能量会无界增长（与外力有界矛盾）。这一步用到 copositive-plus（§4）保证每步 LCP 解的冲量有界，以及能量估计（§3）保证总耗散有界。

为什么需要这一步：BV 估计是后续紧致性的前提。没有一致界，函数列可能"逃逸到无穷"，谈不上收敛。

Step 2：Helly 选择定理（紧致性）。 BV 函数有一个关键性质——Helly 选择定理：

若 $\{f_n\}$ 是一致有界、总变差一致有界的函数列，则存在**逐点收敛的子列** $f_{n_k}\to f$，且极限 $f$ 也是 BV 函数。

应用到 $v_h$：由 Step 1 的一致 BV 界，存在子列 $v_{h_k}$ 逐点收敛到某个 BV 函数 $v$。位置 $q_{h_k}$（Lipschitz，更光滑）一致收敛到 $q$。

为什么需要这一步：这是把"离散解列"变成"有极限"的关键。Helly 定理是 BV 空间的"Bolzano-Weierstrass"——它保证有界闭集列紧。

Step 3：弱 * 紧致（冲量测度的收敛）。 冲量 $P_h$ 是**测度**（含 atom），不能用逐点收敛。但测度空间 $C([0,T])^*$ 上有**弱 * 紧致性**（Banach-Alaoglu 定理）：

有界测度列有弱 * 收敛子列 $P_{h_k}\rightharpoonup P$，即对任意连续函数 $g$，$\int g\,dP_{h_k}\to\int g\,dP$。

由 Step 1 冲量总量有界，故 $P_h$ 有弱 * 收敛子列。极限 $P$ 是一个测度，可能含 atom（这正是冲击/Painlevé 解的来源）。

为什么需要这一步：力是测度（§0 困难 (b)），所以收敛必须在测度的弱 * 拓扑里谈。这一步把"离散冲量"收敛到"连续冲量测度"。

Step 4：极限满足 MDI（通过极限）。 最后把离散格式

\[ M(v_{\ell+1}-v_\ell)=h\bar F_\ell + H^\top P_{\ell+1},\qquad P_{\ell+1}\in -N_{T_C}(\dots) \]

在 $h\to 0$ 取极限。左端 $M\,dv_h\to M\,dv$（Stieltjes 测度收敛），$h\bar F_h\to F\,dt$，$H^\top P_h\to H^\top P$（弱 *）。互补/法锥关系在极限下保持（用法锥的**上半连续性 / 闭图性质**，graph closedness）。于是极限 $(q,v,P)$ 满足

\[ -M\,dv+F\,dt\in N_{T_C(q)}(v^+)\,d\mu, \]

即 Moreau MDI。证毕。

本质洞察：这四步链条 BV 估计 → Helly → 弱 * 紧致 → 通过极限 是非光滑数值分析的"万能模板"。它的精髓是：在**正确的空间**（BV 速度 + 测度冲量）里，离散解列是**紧致**的（有收敛子列），且极限满足**正确的方程**（MDI）。对比光滑数值分析的收敛证明（Taylor 展开 + Gronwall），这里完全不用 Taylor（解不光滑无法展开），而用**紧致性 + 闭性**——这是泛函分析的胜利。任何想为新格式（你自己发明的时步法）证收敛的人，都按这四步走。

Painlevé 解的构造（Step 4 的副产品）¶

Step 3 的极限测度 $P$ 可能含 atom——这正是 Painlevé 悖论"无碰撞冲击"的严格来源。在 Painlevé 配置下，离散冲量 $P_h$ 在某时刻附近"集中"，弱 * 极限是一个 $\delta$ 测度（atom）。物理上：系统在病态配置自发产生瞬时冲量逃离。Stewart 1998 严格构造了这个含 atom 的解，从而**证明了 Painlevé 配置下解存在**（只是不是光滑函数，而是测度解）——彻底消解了百年悖论。

摩擦情形的收敛：Dzonou-Monteiro Marques¶

无摩擦收敛（上面）相对干净，但**带 Coulomb 摩擦**的收敛要难得多，核心障碍是 Coulomb 律的非半连续性（lack of semicontinuity）——摩擦力方向在滑动/粘滞切换处不连续，Step 4 的"法锥闭性"不再直接成立。

Dzonou-Monteiro Marques 2007（Nonlinear Anal. TMA） 的突破：用 de Saxcé 双位锥（§2）把摩擦接触重写成一个**半连续**的包含关系，从而恢复 Step 4 所需的闭性。结论：若摩擦系数 $\mu$ 充分小且 $M$ Lipschitz，带 Coulomb 摩擦的 Moreau-Jean 离散解子列收敛到摩擦 MDI 的解。这是"带摩擦时步法收敛"的第一个严格结果——工程现实（摩擦无处不在）的理论保证。

辛不相容的深层机制（§8 的严格版）¶

§8 给了 Cirak-West 2005 的结论。这里补充**为什么**。

辛积分器保持**辛形式** $\omega=\sum dq_i\wedge dp_i$（相空间的"面积元"）。一个映射 $\Phi$ 辛，当且仅当 $\Phi^*\omega=\omega$（拉回保辛形式）。

接触约束激活时，离散流包含一个**投影**：把穿透/侵入速度投影回可行集 $\{\phi\ge 0,\ \nabla\phi^\top v\ge 0\}$。投影是一个**非满射、非双射**的映射（它把一个区域压到边界上），其雅可比是**奇异的**（投影方向的特征值为 0）。奇异雅可比的映射不可能保持 $\omega$（保辛要求雅可比是辛矩阵，行列式为 1，不可奇异）。所以**任何含投影的接触修正必然破坏辛性**。

唯一例外：$e=1$ 纯弹性 + 无穿透修正——此时碰撞是一个**保能量的速度反射**（法向速度翻转），它是辛的（Kane-Marsden-Ortiz-West 2000 证明变分积分器在此保辛）。但 $e<1$ 一旦引入耗散，反射不再保能量、不再辛。

本质洞察：辛性的破坏不是"数值误差"，而是**几何不相容**——接触修正的投影本质上是"信息丢失"的（多个穿透状态被映射到同一个边界状态），而辛性要求"信息守恒"（双射 + 保面积）。这两者在数学上**不可调和**。这解释了为什么没有任何"改进的辛接触积分器"能完全成功——它撞上的是几何障碍，不是技术难度。

无穷维扩展：柔体与有限元接触¶

前面全是**刚体**（有限维 $M$）。真实机器人有柔性部件（软体抓手、柔性关节），接触发生在**连续介质**上——这把问题推向**无穷维**（$M$ 变成偏微分算子）。

Laursen-Chawla 1997 / Laursen-Love 2002：**能量-动量守恒**的有限元冲击积分器——在 FEM 离散的弹性体接触中，构造保持能量和动量的时间积分（避免数值能量漂移导致的非物理行为）。
Hauret-Le Tallec 2006：动态接触有限元的**能量耗散估计**——量化离散格式在接触时耗散多少能量。
Hesch-Betsch 2009/2011：**mortar 方法**的能量守恒——mortar 是处理不匹配网格接触界面的技术，他们证明了其能量守恒性质。
XPBD / PBD（Müller et al.）：图形学的**位置级投影**求解器——实时但非物理精确，用迭代投影把约束违反"推回去"。VR/游戏柔体的主力，但不适合需要物理精确的机器人力控。

这些属于"接触力学 + 连续介质力学"的交叉，是本专题向 FEM/软体机器人的延伸。关键 take-away：无穷维下"能量守恒"比刚体更微妙——离散弹性波 + 接触很容易引入数值能量漂移，需要专门的能量-动量积分器。

阶段小结：本节补全了三块严格理论——(1) Stewart 1998 收敛四步链（BV 估计→Helly→弱 * 紧致→通过极限），它是非光滑数值分析的万能模板；(2) 辛不相容的几何机制（投影的奇异雅可比 vs 辛的保面积要求）；(3) 向无穷维柔体/FEM 的扩展（能量-动量积分器）。至此本章理论彻底完整。

⚠️ 常见陷阱¶

💡 概念误区：以为"收敛到 MDI 解"是显然的 - 新手想法："步长趋于零，离散解当然趋于连续解，还要证什么？" - 实际上：完全不显然。离散解是跳跃的、力是离散冲量，极限**为什么**恰好满足连续 MDI，需要在正确空间（BV + 测度）里用紧致性 + 闭性严格证明。在错误空间（如要求 $C^1$ 收敛）里，根本不收敛。Stewart 1998 整篇论文就在证这件"显然"的事。 - 为什么重要：理解证明的难度，你才懂为什么"换个格式要重新证收敛"不是走过场——非半连续摩擦（Dzonou-MM）就花了十年才证出来。

🧠 思维陷阱：试图发明"保辛的接触积分器" - 新手想法："Cirak-West 说显式辛不行，那我设计个隐式的、聪明的辛接触积分器。" - 实际上：辛不相容是**几何障碍**（投影的奇异雅可比），不是技术难度。任何修正穿透的操作都破坏辛性，与显式/隐式无关。除非 $e=1$ 纯弹性无修正，否则保辛 + 接触不可调和。把精力花在"凸求解器 + 能量监控"上更务实。 - 正确思维：接受辛性在接触下的失效，转而追求"能量行为可控"（θ 调耗散、能量监控）——这是工程上可达的目标。

💡 概念误区：把刚体收敛结论照搬到柔体/FEM - 新手想法："刚体时步法收敛证明了，柔体有限元接触应该也一样。" - 实际上：无穷维（FEM 弹性体）下"能量守恒"远更微妙——离散弹性波 + 接触极易引入数值能量漂移，需专门的能量-动量积分器（Laursen-Love）。刚体的有限维证明不直接推广到无穷维。 - 为什么重要：做软体机器人/柔性抓手仿真时，不能假设刚体仿真器的能量行为；需用 FEM 接触专用积分器。

练习¶

（证明题，研究级，草稿纸完成） 复现 Stewart 1998 收敛证明的 Step 1–2：(a) 对一个一维质点撞墙系统，写出离散速度 $v_h$ 的总变差，用能量估计证明它一致有界；(b) 陈述 Helly 选择定理，说明它如何给出 $v_h$ 的逐点收敛子列。不要求完整严格，但要说清每步用了什么泛函分析工具及其作用。
（开放思考题） 解释为什么带 Coulomb 摩擦的收敛证明比无摩擦难得多。具体说明"Coulomb 律的非半连续性"是什么（摩擦力方向在哪里不连续），以及 de Saxcé 双位锥如何恢复证明所需的闭性/半连续性。提示：考虑滑动→粘滞切换瞬间摩擦力方向的行为。
（综合题，桥接 §8） 用辛形式 $\omega=\sum dq_i\wedge dp_i$ 和"投影的奇异雅可比"论证辛不相容：(a) 写出一个一维接触投影（把侵入速度 $v<0$ 投影到 $v=0$）的映射及其雅可比；(b) 说明该雅可比为何奇异；(c) 论证奇异雅可比的映射不可能是辛映射（辛矩阵行列式为 1）；(d) 解释 $e=1$ 纯弹性反射（$v\to -v$）为何例外（它的雅可比是什么？辛吗？）。

§12 一个完整的算例：Moreau-Jean 时步法手算 bouncing ball ⭐⭐¶

动机：把抽象格式落到具体数字¶

读完十一节理论，最有效的巩固是**亲手算一步**。本节用 Moreau-Jean θ-scheme 手算一维 bouncing ball 的几步，把 §2–§3 的抽象公式（MDI、Delassus LCP、冲击律嵌入）全部落到具体数字上。这是检验你是否真懂的试金石——若能跟着算完，你就掌握了时步法的核心机制。

设置¶

一维质点（质量 $m=1$ kg）在重力 $g=10\,\mathrm{m/s^2}$ 下落，地面在 $q=0$，gap 函数 $\phi(q)=q\ge 0$，恢复系数 $e=0.5$。初始 $q_0=0.05$ m，$v_0=0$。步长 $h=0.1$ s。取 $\theta=1$（后向 Euler，便于手算）。

接触 Jacobian $H=\nabla\phi=1$（一维，$\phi=q$）。Delassus 算子 $W=HM^{-1}H^\top=1\cdot 1\cdot 1=1$。

逐步格式（θ=1 简化）¶

θ=1 时格式为： $$ m(v_{\ell+1}-v_\ell)=h(-mg)+H^\top P_{\ell+1},\quad q_{\ell+1}=q_\ell+h v_{\ell+1}, $$ 接触互补（速度级 + 冲击律嵌入，仅当 $\phi$ 将 $\le 0$ 时激活）： $$ 0\le \underbrace{(v_{\ell+1}+e v_\ell)}{U\ge 0. $$}} \perp P_{\ell+1

第 1 步（$\ell=0$，自由下落，无接触）¶

当前 $q_0=0.05>0$，且自由速度 $v_{\text{free}}=v_0-hg=0-0.1\cdot 10=-1$ m/s，外推位置 $q_0+h v_{\text{free}}=0.05+0.1\cdot(-1)=-0.05<0$——预测会穿透！ 接触激活。

解 LCP：设 $P_1\ge 0$。$v_1=v_{\text{free}}+W P_1/m=-1+P_1$（$m=W=1$）。冲击律嵌入：$U_1=v_1+e v_0=(-1+P_1)+0.5\cdot 0=-1+P_1$。互补 $0\le U_1\perp P_1\ge 0$：

若 $P_1>0$，则 $U_1=0\Rightarrow P_1=1$。检验 $P_1=1>0$ ✓。
得 $v_1=-1+1=0$。

但等等——这是 $\theta=1$、速度级、$v_0=0$ 的退化情形（初速为零，冲击律 $ev_0=0$）。物理上物体还在 $q_0=0.05$ 上方，不该立即被接触力顶住。这暴露了纯速度级格式 + 大步长的一个问题：它用"外推穿透"判断激活，$h=0.1$ 太大导致提前激活。这正是 §9 讲的离散副作用。

改用更细步长 $h=0.01$ 重算第 1 步：$v_{\text{free}}=0-0.01\cdot 10=-0.1$，外推 $q=0.05+0.01\cdot(-0.1)=0.049>0$——不穿透，接触不激活。$P_1=0$，$v_1=-0.1$ m/s，$q_1=0.049$ m。合理：物体自由下落一点点。

教训：这个"算砸了又算对"的过程本身极有教育意义——它具体展示了 §9 的"步长过大导致提前激活/虚假接触"。大步长时速度级格式的激活判据失真，必须减小 $h$ 或用位置级（Stewart-Trinkle 的 $\phi_\ell/h$ 项能更准地判断激活）。

接触步（物体已接近地面）¶

跳到物体即将落地的某步：设 $q_\ell=0.005$ m，$v_\ell=-0.99$ m/s（自由下落积累的速度），$h=0.01$。

自由速度 $v_{\text{free}}=-0.99-0.01\cdot 10=-1.09$ m/s。外推 $q=0.005+0.01\cdot(-1.09)=0.005-0.0109=-0.0059<0$——穿透，接触激活。

解 LCP：$v_{\ell+1}=v_{\text{free}}+P_{\ell+1}=-1.09+P_{\ell+1}$。冲击律嵌入 $U_{\ell+1}=v_{\ell+1}+e v_\ell=(-1.09+P_{\ell+1})+0.5\cdot(-0.99)=-1.09+P_{\ell+1}-0.495=P_{\ell+1}-1.585$。

互补 $0\le U_{\ell+1}\perp P_{\ell+1}\ge 0$：设 $P_{\ell+1}>0\Rightarrow U_{\ell+1}=0\Rightarrow P_{\ell+1}=1.585$ N·s。检验 $>0$ ✓。

新速度 $v_{\ell+1}=-1.09+1.585=0.495$ m/s。验证冲击律：$v^+=-e v^-$ 期望 $=-0.5\cdot(-0.99)=0.495$ m/s ✓✓。

冲量分解：$P=1.585=\underbrace{1.09}_{\text{抵消自由穿入}}+\underbrace{0.495}_{\text{产生回弹}}$，即冲量先抵消了入射动量（$m|v_{\text{free}}|=1.09$），再产生回弹动量（$m v^+=0.495$）。

新位置 $q_{\ell+1}=0.005+0.01\cdot 0.495=0.005+0.00495=0.0099$ m $>0$——物体回弹离地 ✓。

本质洞察：这一步把全章核心机制压缩成几个数字——(1) 未知量是冲量 $P$（不是力）；(2) 冲击律 $v^+=-ev^-$ 通过 argument $v_{\ell+1}+ev_\ell$ 自动满足，无需单独"应用碰撞规则"；(3) 同一个 LCP 既处理了"抵消入射"又处理了"产生回弹"——持续接触与冲击统一。亲手验证 $v^+=0.495=-0.5\times(-0.99)$ 那一刻，你就彻底理解了 §2–§3 的全部抽象公式。

能量核对¶

碰撞前动能（用 $v^-=-0.99$）：$E^-=\tfrac12\cdot 1\cdot 0.99^2=0.490$ J。碰撞后（$v^+=0.495$）：$E^+=\tfrac12\cdot 1\cdot 0.495^2=0.1225$ J。能量比 $E^+/E^-=0.25=e^2$ ✓——恢复系数 $e=0.5$ 对应能量保留 $e^2=0.25$，耗散 $75\%$。这验证了 §7 的能量性质：$e<1$ 碰撞耗能，耗散比为 $1-e^2$。

练习¶

（计算题，草稿纸完成） 接着上面的接触步，继续手算 2 步：物体回弹后自由上升、减速、再下落。验证每步的 $q,v$ 演化合理（上升时 $v>0$ 递减，到最高点 $v=0$，再下落 $v<0$）。物体第二次落地时速度应约为第一次回弹速度 $0.495$ 的多少？（提示：自由飞行能量守恒，落回同高度时 $|v|$ 应接近回弹速度。）
（对比题） 把同一个接触步改用 $\theta=0.5$ 重算（位置更新用 $\tfrac12(v_{\ell+1}+v_\ell)$）。比较 $\theta=1$ 与 $\theta=0.5$ 的结果差异，特别是能量行为——哪个更接近理论的 $e^2=0.25$ 能量比？为什么？（联系 §3 的能量估计。）
（综合题，桥接 §5） 把这个一维例子改用 Anitescu 凸松弛（§5）重新表述：写出每步的凸 QP（一维，无摩擦时退化为简单的投影），说明它与 LCP 在一维无摩擦情形下是否等价。再加入一个切向（假设地面有摩擦，物体有水平初速），说明此时 LCP 与凸松弛开始出现 $O(h)$ gliding 差异。

本章常见误解汇总¶

误解	正确理解	相关节
接触仿真 = 光滑 ODE + 碰撞检测打补丁	是 measure differential inclusion 的时间离散；接触与动力学同时求解	§0, §1
时步法是事件驱动法的近似/偷懒版	时步法是独立数学对象（扫动过程原生离散），$h\to 0$ 极限就是 MDI 解；Zeno 下事件驱动无定义	§1, §4
高阶积分器（RK4）更适合接触	冲击点解是 BV 的，所有积分器塌陷到 $O(h)$；评判标准是收敛性/能量/穿透而非阶数	§0, §8
接触力 $\lambda_n$ 和冲击是同一个量	持续接触力是有限函数，冲击是冲量 atom；时步法未知量统一是冲量 $P$	§0, §3
MDI 就是"微分方程 + Dirac 项"	$dv$ 是微分测度，内禀含跳跃；跳跃时刻/大小由包含关系自动决定	§2
多面体摩擦锥 = Anitescu 凸松弛	两个正交近似：多面体改锥形状（$O(1/d)$，LCP 精确），凸松弛改约束结构（$O(h)$，换凸性）	§4, §5
凸松弛丢物理，所以不如 LCP 准	凸松弛误差 $O(h)$ 随步长消失（收敛到 MDI）；换来唯一解+可微+扩展，机器人任务净收益为正	§5
LCP 可解 = 解唯一	copositive-plus 只保证存在；多体同时冲击下 LCP 可多解，换求解器结果变	§4, §7
MuJoCo 是硬接触 + 碰撞检测	MuJoCo 内禀软接触——接触力是 gap 的光滑函数，无二值开关，故可微	§6
`solref` 和 `solimp` 都是调软硬，随便设一个	`solref` 定义参考加速度（约束"想做什么"），`solimp` 定义阻抗随违反量变化（"能做到什么"）	§6
恢复系数 $e$ 是材料固有常数	$e$ 依赖速度/几何/温度，且 Newton/Poisson/Stronge 定义下数值不同	§7
$e=1$ 完全弹性一定能量守恒	单点成立；多体同时碰撞 + Newton 法则可能增能（Kane 反例），需能量级冲击律	§7
保辛积分器更适合接触	Cirak-West 2005：辛性与接触激活/穿透修正不相容（除非 $e=1$）；机器人 $e\approx 0$ 塑性，辛失效	§8
半隐式 Euler 和显式 Euler 差不多	半隐式保辛（能量长期有界），显式注入能量（发散）；仅换位置更新顺序即质变	§8
所有穿透都是 bug	软接触穿透是设计特性（受 `solref`/`solimp` 控制）；速度级漂移才需投影修正	§6, §9
jitter 是物理微振动	静止堆叠 jitter 几乎总是数值伪影（接触力高频切换、求解器未收敛）	§9
换仿真器结果应该一致	不同离散/锥/松弛/求解器导致系统性差异；可复现需固定全套配置	§10

本章小结¶

符号表¶

本章新引入的数学符号及含义：

符号	含义	首次出现
$\phi(q)$	gap 函数（有符号距离），$\phi\ge 0$ 不穿透	§0
$M$	广义质量矩阵（正定 $M\succ 0$）	§0
$H=\nabla\phi$（或 $J$）	接触 Jacobian（gap 梯度）	§0, §6
$\lambda_n$	法向接触力（$\ge 0$）	§0
$P$	接触冲量（$=\int d\mathbf{R}$，量纲 $N\cdot s$）	§0, §3
$e$	恢复系数（restitution，$\in[0,1]$）	§0
$d\mu$	基测度 $=dt+\sum_i\delta_{t_i}$（勒贝格 + 原子）	§2
$dv$	速度的微分测度（含跳跃）	§2
$d\mathbf{R}$	接触冲量测度	§2
$v^+,v^-$	速度的右/左极限（$\lim_{s\downarrow t}, \lim_{s\uparrow t}$）	§2
$C(q)$	可行集 $\{q:\phi(q)\ge 0\}$	§2
$T_C(q)$	切锥（tangent cone）	§2
$N_K(x)$	凸锥 $K$ 在 $x$ 的法锥（normal cone）	§2
$\theta$	θ-method 加权系数（$\in[0,1]$）	§3
$W=HM^{-1}H^\top$	Delassus 算子（接触空间逆惯量）	§3
$U$	接触点相对速度（含恢复 $v^++ev^-$ 的法向投影）	§3
$v_{\text{free}}$	自由速度（无接触时一步速度）	§3
$D=[d_1,\dots,d_d]$	切向方向基（多面体摩擦锥）	§4
$\boldsymbol\beta$	各方向摩擦冲量（$\ge 0$）	§4
$\sigma$	滑移率乘子（slip rate multiplier，$\ge 0$）	§4
$E=[1,\dots,1]^\top$	滑移率耦合向量	§4
$M_{\text{LCP}}$	Stewart-Trinkle LCP 矩阵	§4
$\mu$	摩擦系数	§4
$\mathcal{K},\mathcal{K}^*$	摩擦锥及其对偶锥	§5, §6
$N=D^\top M^{-1}D$	Anitescu 凸 QP 的 Hessian（= Delassus）	§5
$\boldsymbol\gamma$	接触冲量向量（凸 QP 未知量）	§5
$A=JM^{-1}J^\top$	MuJoCo 对偶 Hessian（= Delassus 连续版）	§6
$a_{\text{ref}}$	参考加速度（reference acceleration）	§6
$R$	软约束正则化矩阵（impedance）	§6
`solref`	MuJoCo 约束时间常数/阻尼比参数	§6
`solimp`	MuJoCo 阻抗（impedance）参数	§6
$\alpha,\beta$	Baumgarte 稳定化的阻尼/频率参数	§9

定理速查表¶

本章核心定理/公式及一句话说明：

定理/公式	一句话说明	对应节	出处
Moreau MDI	$-M\,dv+F\,dt\in N_{T_C(q)}(v^+)\,d\mu$；持续接触 + 冲击统一公式	§2	Moreau 1988 CISM 302
Moreau-Jean θ-scheme	速度级离散 + 冲击律嵌入 argument $v_{\ell+1}+ev_\ell$；每步解 Delassus LCP	§3	Moreau/Jean CMAME 177 (1999)
θ 能量估计	$\theta=\tfrac12$ 能量中性；$\theta=1$ 耗散 $O(h^2)$	§3	经典 θ-method 分析
Delassus 半正定	$W=HM^{-1}H^\top\succeq 0$（因 $M\succ 0$）⇒ 法向 LCP 可解	§3	LCP 理论
Stewart-Trinkle LCP	$d$-方向多面体锥 + 滑移率乘子复现最大耗散；矩阵见 §4	§4	Stewart-Trinkle IJNME 39 (1996) 2673
Anitescu-Potra copositive-plus	$M_{\text{LCP}}$ copositive-plus ⇒ Cottle-Pang-Stone ⇒ 每步可解	§4	Anitescu-Potra Nonlinear Dynam. 14 (1997) 231
Stewart 1998 收敛 + Painlevé	$h\to 0$ 弱 * 收敛到 MDI 解；Painlevé 含 atom 冲量解（impact without collision）	§4	Stewart ARMA 145 (1998) 215, DOI 10.1007/s002050050129
Anitescu 2006 凸松弛收敛	凸 QP 松弛（$\phi/h$ 稳定项）$\epsilon\to 0$ 收敛到 DVI/MDI 解	§5	Anitescu Math. Program. 105 (2006) 113
CCP 形式	$\mathcal{K}\ni\gamma\perp(N\gamma+r)\in\mathcal{K}^*$；投影梯度可解	§5	Tasora-Anitescu 2010/2011
MuJoCo Gauss 凸优化	$\min\tfrac12\\|x-M^{-1}(\tau-c)\\|_M^2+\dots$；对偶 Hessian $A+R$	§6	MuJoCo 文档；Todorov 2012/2014
Moreau-Jean 弱收敛	prox-regular 约束 + $\theta\in[\tfrac12,1]$，$q_h\to q$ 一致、$v_h\overset{*}{\rightharpoonup}v$	§3	Moreau/Jean CMAME 177 (1999)
Dzonou-Monteiro Marques	带 Coulomb 摩擦 Moreau-Jean 弱收敛（de Saxcé 双位锥克服非半连续）	§2, §7	Dzonou-MM Nonlinear Anal. TMA (2007)
Cirak-West 不相容	显式辛格式在接触激活/穿透修正下必破坏辛性	§8	Cirak-West IJNME 64 (2005) 1079
Kane-Marsden-Ortiz-West	变分积分器在弹性接触（$e=1$）下保辛的条件	§8	KMOW IJNME 49 (2000) 1295
Castro SAP 全局收敛	SAP 势函数严格凸 + Lipschitz-Hessian ⇒ Newton 全局线性-二次收敛	§10	Castro-Permenter-Han T-RO 39 (2022) 1301

知识点总表¶

编号	知识点	核心要点	对应节	难度
1	三重困难	BV 解 / 测度力 / 移动约束；经典积分器失效	§0	⭐⭐
2	事件驱动 vs 时步法	检测事件 vs 每步解互补；Zeno/Painlevé 下时步法唯一有定义	§1	⭐⭐
3	Moreau MDI	$-M\,dv+F\,dt\in N_{T_C}(v^+)\,d\mu$；测度 + 法锥统一持续/冲击	§2	⭐⭐⭐
4	Moreau-Jean θ-scheme	速度级 + 冲击嵌入 argument；θ 控能量；解 Delassus LCP	§3	⭐⭐⭐
5	Stewart-Trinkle LCP	位置级 + 多面体锥；copositive-plus 每步可解；消解 Painlevé	§4	⭐⭐⭐
6	Anitescu 凸松弛	LCP→凸 QP/CCP；$O(h)$ gliding 换唯一解+可微+扩展	§5	⭐⭐⭐
7	MuJoCo 软接触	Gauss 原理 + 正则化凸优化；`solref`/`solimp`；PGS/CG/Newton	§6	⭐⭐⭐
8	冲击律	Newton/Poisson/Stronge/Moreau；多体同时冲击非唯一	§7	⭐⭐
9	半隐式 vs 全隐式 + 辛不相容	三层积分器；Cirak-West 辛与接触不相容	§8	⭐⭐⭐
10	穿透/能量/漂移稳定化	Baumgarte（脏）vs GGL/投影（干净）；jitter 是伪影	§9	⭐⭐
11	仿真器对比	六大仿真器 / 四条路线；全部基于 Delassus 算子	§10	⭐⭐

累积项目：本章新增模块¶

接触力学累积项目：从零搭一个最小但完整的接触仿真器，逐章累积。本章新增"时间离散引擎"模块。

本章交付物：实现一个支持**两种时步格式**的 2D 刚体接触仿真器核心（不依赖现成物理引擎）：

模块 A — Moreau-Jean θ-scheme（速度级）：实现 §3 的更新公式，每步用 PGS 解 Delassus LCP $0\le U=WP+U_{\text{free}}\perp P\ge 0$。支持 $\theta\in\{0.5,1\}$ 切换，验证能量行为（$\theta=0.5$ 弹性碰撞链保能，$\theta=1$ 耗散）。
模块 B — Anitescu 凸松弛（凸 QP）：实现 §5 的每步凸 QP，用投影梯度（APGD）求解 CCP，二阶锥摩擦。测量并验证 $O(h)$ gliding（斜面静置方块的下滑速度正比于 $h$）。
基准场景：(a) bouncing ball（验证 Zeno 鲁棒性——多次微弹被一步吸收）；(b) 立方体堆叠（验证稳定性 + 穿透 + jitter）；(c) 斜面方块（验证 gliding）。
对照实验：把同样三个场景接到 MuJoCo / PyBullet，对比你的实现与它们的穿透量、能量曲线、gliding——理解"换格式结果变"（§10）。

与前序模块衔接：模块复用本目录 10 的 LCP/PGS 求解器、本目录 20 的摩擦锥投影、第四批的 $M(q)$ 与 Jacobian 计算。与后续衔接：本章的凸 QP 模块（B）是本目录 40_可微接触的基础——下一章在它上面加自动微分/隐函数定理求梯度。

项目代码目录：contact_sim/timestepping/（建议 Python + NumPy 起步，后续可移植 C++/JAX）。

延伸阅读¶

核心教材（按角色分类，标注难度）¶

教材	角色	时步法章节	难度	一句定位
Acary & Brogliato Numerical Methods for Nonsmooth Dynamical Systems (2008, Springer LNACM 35)	首选	Ch. 7–11（Part III）；Ch. 12–14	★★★★★	非光滑时步法的"工程师百科全书"；Moreau-Jean/Schatzman-Paoli/event-driven 全覆盖 + 求解器细节
Brogliato Nonsmooth Mechanics 3rd (2016, Springer CCE)	建模参考	§5.7 数值综述；§7.3 数值稳定	★★★★★	建模—适定性—稳定性—控制 one-stop；约 1300 条参考文献
Stewart Dynamics with Inequalities (2011, SIAM)	数学最严	Ch. 5–8	★★★★★	唯一给 Stewart-Trinkle 收敛完整证明的专著
Pfeiffer & Glocker Multibody Dynamics with Unilateral Contacts (1996, Wiley)	历史奠基	Part 1 LCP；Ch. 11–14	★★★★	LCP 接触动力学在机械工程中的奠基
Hairer-Lubich-Wanner Geometric Numerical Integration 2nd (2006, Springer SSCM 31)	对照组	Ch. VI 辛；VII.1 约束；IX backward error	★★★★★	辛/变分积分器圣经（光滑背景）——让你理解为何辛在接触下失效

论文与讲义¶

奠基论文： - Stewart & Trinkle, IJNME 39 (1996) 2673（位置级 LCP 时步法） - Anitescu & Potra, Nonlinear Dynam. 14 (1997) 231（copositive-plus 可解性） - Stewart, ARMA 145 (1998) 215, DOI 10.1007/s002050050129（收敛 + Painlevé 消解） - Anitescu, Math. Program. 105 (2006) 113（凸松弛） - Moreau, CISM 302 (1988)；Moreau/Jean, CMAME 177 (1999)（扫动过程 + θ-scheme） - Tasora & Anitescu (2010/2011)（matrix-free CCP / APGD，granular） - Castro-Permenter-Han, IEEE T-RO 39 (2022) 1301, DOI 10.1109/TRO.2022.3209077（Drake SAP 全局收敛）

课程与文档： - Russ Tedrake, Underactuated Robotics "Planning and Control through Contact"：https://underactuated.mit.edu/contact.html - MuJoCo 官方文档 Computation 章：https://mujoco.readthedocs.io/ - Drake 接触文档：https://drake.mit.edu/doxygen_cxx/group__multibody__contact.html - Pinocchio 文档：https://stack-of-tasks.github.io/pinocchio/ - Vincent Acary (INRIA TRIPOP)：https://hal.inria.fr/ （Siconos 技术报告、收敛证明预印本） - Mihai Anitescu (Argonne)：https://www.mcs.anl.gov/~anitescu/ （Anitescu-Potra-Stewart 论文 PDF）

开源代码库：Siconos、Drake、MuJoCo/MJX、Project Chrono、Bullet、Pinocchio、ContactBench（Le Lidec 2024 基准）、SimBenchmark（Hwangbo 2018 跨仿真器基准）、Simple（Carpentier-Le Lidec-Montaut 2024 ADMM 统一求解器）。

以上链接在域名层面真实存在；具体子页 URL 可能随时间变动，建议从主页导航而非硬编码深链接。中文资源（知乎"MuJoCo 求解器"/"Drake SAP"/"非光滑动力学"、B 站 GAMES103/202、公众号"泡泡机器人 SLAM"）请以关键词自行检索并核对时效。

本章与后续章节的关系¶

后续章节	与本章的关系	本章哪个知识点为其铺垫
本目录 40_可微接触	本章凸松弛（§5）/软接触（§6）的凸性 ⇒ 解对参数可微 ⇒ 解析雅可比	§5 凸 QP、§6 MuJoCo soft-convex、§5 桥接小节
混合系统（hybrid systems）	事件驱动（§1）是混合系统特例；时步法是混合系统的"消元"（把模式切换隐式化）	§1 两条路线
接触隐式 MPC（CI-MPC）	Posa-Cantu-Tedrake 2014、Manchester contact-implicit TO 直接用时步离散作为优化约束	§3/§4 离散格式、§5 凸 QP
足式运动控制	落地接触的 $e\approx 0$ 塑性、多足同时触地的多体冲击	§7 冲击律、§8 辛不相容
强化学习（仿真环境）	MJX 百万并行 rollout、可微 rollout 求策略梯度	§6 MuJoCo soft-convex、§5 可微性

🔧 故障排查手册¶

#	症状	可能原因	排查步骤	相关节
1	物体穿过地板/相互嵌入	(a) 步长过大；(b) 速度级格式漂移；(c) 软接触固有穿透；(d) 碰撞检测漏检	① 减小步长 $h$ 看穿透是否减小（若线性减小 → 离散误差）；② 判断仿真器是软接触（MuJoCo，穿透正常）还是刚性（需投影）；③ 速度级格式加坐标投影；④ 检查碰撞 margin 设置	§0, §6, §9
2	仿真在某时刻"卡死"/步长缩到机器精度	用了事件驱动法 + Zeno 场景（bouncing ball 落定）	① 确认是否事件驱动（root-finding）；② 若是，改用 time-stepping（固定步长不检测事件）；③ time-stepping 会把多次微弹自动吸收	§0, §1
3	能量莫名增加/系统越跑越快	(a) 显式 Euler 注入能量；(b) Baumgarte 稳定化注入能量；(c) 求解器未收敛	① 检查积分器，改半隐式以上；② 弃 Baumgarte，改坐标投影/GGL；③ 增大求解器迭代数；④ 检查 $e>1$ 误设	§8, §9
4	能量过度损失/运动"黏稠"变慢	(a) $\theta=1$ 或大数值阻尼；(b) Baumgarte 过阻尼	① θ-scheme 改 $\theta=0.5$（能量中性）；② 减小数值阻尼/`solref` 阻尼比；③ 检查接触 $e$ 是否过小	§3, §8
5	静止堆叠物体高频抖动（jitter）	(a) 步长过大；(b) 刚度过高 vs 步长不匹配；(c) 求解器未收敛；(d) 接触力高频切换	① 增大求解器迭代数（最常见）；② 减小步长或降低刚度；③ 适当增大数值阻尼（$\theta\to 1$）；④ warm-start 稳定接触力	§6, §9
6	物体在斜面上本应静止却缓慢下滑	凸松弛的 $O(h)$ gliding（非 bug，非参数错误）	① 确认是否凸格式（SAP/MuJoCo/Chrono）；② 减小步长 $h$ 看滑移是否线性减小（确认是 gliding）；③ 接受它或减小 $h$；④ 勿靠调大摩擦系数硬修	§5
7	换求解器/换仿真器结果不同	(a) 多体同时冲击 LCP 解非唯一；(b) 不同摩擦锥/凸松弛；(c) 大步长 gliding 差异	① 确认场景是否多体同时冲击（解非唯一是正常的）；② 对照 §10 各仿真器格式差异；③ 要可复现须固定仿真器+版本+求解器+步长	§4, §7, §10
8	接触求解器报"不收敛"/失败	(a) Drake TAMSI 非凸失败；(b) ill-conditioned Delassus；(c) 摩擦系数过大触发病态	① 改用凸求解器（SAP/ADMM）保证全局收敛；② Pinocchio ADMM + proximal 正则化对付 ill-conditioned；③ 检查质量矩阵/接触配置是否病态	§5, §10

研究实践建议¶

给新手的建议：

先跑 bouncing ball 和立方体堆叠——这两个最小场景能暴露 §9 几乎所有病理（穿透、漂移、jitter、能量漂移）。在自己实现的简单时步法上调通它们，比读十篇论文更能建立直觉。
从 MuJoCo 入手而非 Siconos——MuJoCo 文档完善、API 友好、社区大，先建立"软接触 + 凸优化"的工程直觉，再回头看 Siconos 的数学严谨性。
永远报告你的仿真配置——仿真器、版本、求解器、步长、摩擦锥类型。接触仿真结果对这些极敏感（§10），不报告等于不可复现。
把"四条主线 + Delassus 算子"当成坐标系——遇到任何新仿真器/论文，先定位它在四条路线（§2–§6）的哪条、约束位置级还是速度级、锥是多面体还是 SOC、是否凸松弛。这个坐标系能让你快速归类一切。

给有经验者的建议：

可微接触是当前前沿——若做 RL/MPC，深入 §5–§6 的凸性 → 可微（本目录 40），关注 Pinocchio 3 / Simple 的解析可微、ContactBench 的基准。这是 2020s 最活跃的方向。
读 Stewart 1998 的收敛证明全链（档位 4）——BV 估计 → Helly → 弱 * 紧致 → MDI。这是理解"为什么时步法收敛到正确物理"的唯一途径，也是写新格式收敛证明的模板。
关注大步长 MPC 的 gliding 权衡——Drake SAP 的 $O(h)$ gliding 在 $h\sim 5$–$10$ ms 可见，但换来 5–10× 步长。在接触隐式 MPC 里这个权衡是核心设计决策。
**多体同时冲击的非唯一性**仍是开放问题——Frémond 多体冲击律（Acary 等 2024）、能量一致投影是前沿，若做精密多接触操作值得深入。

跨章综合练习¶

以下练习需综合本章多节及前置/后续章节，目的是打破小节隔阂、把知识连成网络（R14 要求）。

（贯穿全章综合题） 从一个具体任务出发——"仿真一只四足机器人在沙地（granular）上小跑（trot）"——回答：(a) 应选事件驱动还是时步法？为什么（§1，考虑接触数 + Zeno）？(b) 沙地颗粒间接触用哪种格式（§3–§6，考虑接触数 $10^5$+）？(c) 机器人脚与地面的塑性落地（$e\approx 0$）为何排除辛积分器（§8）？(d) 若要用 RL 训练这个步态，仿真器该选哪个、为什么需要可微（§5–§6 + 本目录 40）？把这四问串成一个连贯的"为什么 MuJoCo/MJX 是这个任务的正确选择"的论证。
（综合 §2 §3 §4 §5，需结合本目录 10/20） 沿着"Delassus 算子 $W=JM^{-1}J^\top$"这条主线，证明四种格式的统一性：(a) 在 Moreau-Jean（§3）里 $W$ 出现在哪（LCP 矩阵）？(b) 在 Stewart-Trinkle（§4）里 $W$ 的子块如何组成 $M_{\text{LCP}}$？(c) 在 Anitescu 凸 QP（§5）里 $W$ 是什么（Hessian $N$）？(d) 在 MuJoCo（§6）里 $W$ 的连续版 $A=JM^{-1}J^\top$ 在对偶问题哪里？写一段话论证"四种格式是同一个 $W$ 的四种包装"，并指出它们的真正区别在哪四个维度（约束级别、锥形状、是否凸松弛、求解算法）。
（综合 §7 §8 §9，开放设计题） 你要为一个新仿真器设计时间积分 + 冲击处理方案，目标任务是"人形机器人多接触操作 + 大步长 MPC"。请做四个设计决策并各用一节的理论证成：(a) 用什么冲击律（§7）？(b) 用半隐式还是全隐式（§8，考虑大步长稳定性界 $h\le 2/\omega_0$）？(c) 如何防约束漂移（§9，Baumgarte vs 投影）？(d) 凸松弛的 $O(h)$ gliding 在大步长下可接受吗（§5），如何缓解？把四个决策组装成一个连贯的设计文档大纲。
（前后桥接题，桥接本目录 10/20/40） 画一张"本章在接触力学知识体系中的位置"图：上游（本目录 10 LCP/NCP 提供"每步解什么"、本目录 20 摩擦锥提供"可行集形状"、第二批凸分析提供"法锥/近端语言"、第四批刚体动力学提供 $M$ 和 $J$），本章（把这些组装成时间离散格式），下游（本目录 40 可微接触用本章的凸性、CI-MPC 用本章的离散作约束）。对每条上游/下游边，用 2–3 句话说清"本章如何复用上游、下游如何复用本章"。

结语¶

学完本专题，你应具备三项核心能力：

一，数学层面——能把任意接触动力学问题写成 measure differential inclusion（Moreau MDI），并导出 Moreau-Jean（速度级）/ Stewart-Trinkle（位置级 LCP）/ Anitescu（凸松弛）三条离散路线，读懂 Stewart 1998 ARMA 的收敛证明骨架，理解 copositive-plus 为何保证每步可解。

二，工程层面——能打开 Drake sap_solver.cc、MuJoCo engine_solver.c、Siconos MoreauJeanOSI.cpp 任一源码，指出它在哪条数学路线、冲击律如何嵌入、步长稳定性来自何处，并认出处处出现的 Delassus 算子 $W=JM^{-1}J^\top$。

三，决策层面——面对具体机器人任务（抓取 / 足式 / granular / RL 训练 / 微分 MPC），能按"冲击刚度、接触数、可微性需求、实时约束"四维度选出正确仿真器与求解器。

这把钥匙打开后，本目录 40_可微接触与后续的接触隐式 MPC 将不再是黑箱——它们都建立在你刚刚掌握的 θ-scheme、凸 primal 求解器、以及"凸性 ⇒ 可微"这条主线之上。本章反复强调的一句话值得记一辈子：接触不是打断光滑世界的奇异事件，而是测度过程的自然分量；冲击律不是外挂的碰撞规则，而是时步格式 argument 里的一个系数。理解了这句话，你就理解了整个非光滑接触动力学的灵魂。

最后给一张"全章一图流"的心智地图，帮你把九节理论压缩成一条主线随时调取：接触动力学的解活在 BV/测度空间（§0）→ 正确方程是 Moreau MDI（§2）→ 离散有速度级（Moreau-Jean，§3）和位置级（Stewart-Trinkle，§4）两条，都解以 Delassus 算子 $W=JM^{-1}J^\top$ 为核的互补问题 → 凸松弛把 LCP 升级为凸 QP/CCP（Anitescu，§5）换来可微 + 扩展 → MuJoCo 进一步软化 + 正则化（§6）换来处处可微 + 百万并行 → 冲击律嵌入 argument（§7）、半隐式 Euler 主力（§8，辛与接触不相容）、穿透/jitter/漂移靠投影稳定化（§9）→ 六大仿真器全是这条主线上的点（§10）。把这条主线记牢，任何接触仿真问题你都能定位它在哪一环、该用哪个工具、为什么。这就是从"会用仿真器"到"懂仿真器"的跨越。

并且记住本章的元方法论：每当遇到一个新的接触仿真器、新的论文、新的求解器，不要从零开始读，而是用本章建立的**四维坐标系**去定位它——(1) 约束写**位置级还是速度级**（§2）？(2) 摩擦锥用**多面体还是二阶锥**（§4 vs §5）？(3) 是否做了**凸松弛/软化**（§5–§6）？(4) 用什么**数值算法**解（Lemke / PGS / APGD / ADMM / Newton）？四个问题答完，这个新对象就被你精确地放进了知识地图，剩下的只是细节。这套"用坐标系归类新事物"的能力，比记住任何单一格式都更持久、更值钱——它是把一章教学文档真正内化为研究能力的标志。

版本信息速查¶

本章涉及的工具/库/求解器及其关键版本信息（截至 2026，请以官方仓库为准）：

工具/库	角色	关键模块	备注
MuJoCo	soft-convex 仿真	`engine_solver.c`（PGS/CG/Newton）、`engine_core_constraint.c`	3.x 已开源（DeepMind）；MJX 为 JAX 版，支持可微 + 百万并行
MJX	MuJoCo JAX 版	`mjx/_src/solver.py`	可微 rollout，RL 大批量训练
Drake	TAMSI + SAP	`tamsi_solver.cc`、`contact_solvers/sap/sap_solver.cc`	SAP 为当前默认接触求解器；凸 primal
Siconos	Moreau-Jean θ + event-driven	`MoreauJeanOSI.cpp`、`FrictionContact/`	数学最严谨；求解器最全（NSGS/AC/IPM/Proximal）
Project Chrono	Anitescu CCP / APGD	`ChSolverAPGD.cpp`、`ChSolverADMM.cpp`	granular $>10^6$；GPU 扩展
Bullet	Sequential Impulse (PGS)	`btSequentialImpulseConstraintSolver.cpp`	游戏/快速原型；warm-start
Pinocchio 3 / Simple	ADMM + proximal	`admm-solver.hxx`、`contact-cholesky.hpp`	解析可微；统一刚性/柔性 NCP/SOCCP
ContactBench	求解器基准	Le Lidec 2024 配套	跨求解器精度/性能基准
SimBenchmark	跨仿真器基准	Hwangbo 2018	量化仿真器间系统性差异

关键参数默认值速记（MuJoCo）：solref=(0.02, 1.0)（timeconst 0.02 s、临界阻尼）、solimp=(0.9, 0.95, 0.001, 0.5, 2)、默认求解器 Newton、默认积分器半隐式 Euler。稳定性经验界：timeconst $\ge 2h$，步长 $h\le 2/\omega_0=2\sqrt{m/k}$。

注：接触仿真结果对仿真器版本、求解器选择、步长、摩擦锥类型高度敏感（§10）。任何可复现的研究都必须报告完整配置（仿真器 + 版本 + 求解器 + 步长 + 摩擦锥）。上述源码路径在对应仓库主分支真实存在，但具体文件可能随版本演进微调。

学习路径提示：本章是接触力学专题的"时间离散"环节。建议的完整学习顺序为——先读本目录 10（LCP/NCP，每步解什么）、20（摩擦锥，可行集形状），再读本章（30，时间离散把前两者组装成可跑的迭代），最后读 40（可微接触，利用本章凸性求梯度）。本章的两个"必备附加组件"——符号表与定理速查表（见本章小结）——是回头速查时的入口；遇到任何符号或定理不记得，先翻这两张表，再回对应小节精读。

层次	困难	数学刻画	后果
(a) 解的正则性	速度有跳跃	\(v(\cdot)\) 是 BV 函数（bounded variation），\(q(\cdot)\) 仅 Lipschitz	经典积分器在跳跃点阶数塌陷
(b) 力的本质	冲击力是"无穷大作用无穷短时间"	作用力是测度（measure）而非函数，冲量含 atomic measure（原子测度 \(\delta_{t_i}\)）	不能写 \(\lambda(t)\)，只能写 \(d\mathbf{R}\)
(c) 约束集移动	可行域 \(C(q)\) 随位形变	\(C(q)=\{q:\phi(q)\ge 0\}\) 是移动的凸集（sweeping process）	约束的激活/失活是不连续的几何事件

	位置级（position-level）	速度级（velocity-level）
约束	\(\phi(q)\ge 0\)（gap 非负）	\(\phi=0\) 时 \(\nabla\phi^\top v^+\ge 0\)（不准侵入速度）
保证	严格不穿透	速度合法，但位置可能漂移
代表格式	Stewart-Trinkle（§4）	Moreau-Jean（§3）
代价	需非线性 gap 函数展开	需位置投影防漂移
冲击律嵌入	不直接	\(v^++ev^-\) 直接进 argument

	扫动过程（硬接触，§2–§5）	罚函数 / 软接触（§6 MuJoCo）
接触力来源	互补涌现（恰好够用）	函数式给定（弹簧 \(\lambda=-k\phi\)）
穿透	严格零（或漂移修正）	有意允许（\(\phi<0\) 产生力）
数学对象	法锥包含 / LCP	光滑函数
可微性	激活集切换处不可微	处处可微
刚度	无穷硬	有限 \(k\)（可调）

θ	光滑段	接触段	总体行为	适用
\(\tfrac12\)	能量中性	接触耗散（\(e<1\)）	二阶精度、能量中性	弹性碰撞、长时间守恒轨迹
\(1\)	数值耗散 \(O(h^2)\)	接触耗散	一阶精度、强耗散	刚性接触、需阻尼抑制振荡
\((\tfrac12,1)\)	中等耗散	接触耗散	介于两者	工程折中

	Stewart-Trinkle	Moreau-Jean
约束级别	位置级 \(\phi(q_{\ell+1})\ge 0\)	速度级 \(\nabla\phi^\top(v_{\ell+1}+ev_\ell)\ge 0\)
穿透	严格不穿透（线性化后近似）	速度合法、位置可能漂移
gap 评估	需要（非线性 \(\phi\) 展开）	不需要
摩擦锥	多面体 → LCP	切锥法锥（可 SOC 或多面体）
可解性	copositive-plus（Anitescu-Potra）	\(W\succeq 0\) 半正定
漂移	无（位置级直接约束）	需位置投影
典型实现	早期 Drake、Bullet MLCP、ODE	Siconos MoreauJeanOSI

	近似的对象	近似类型	误差	代数后果
多面体锥（Stewart-Trinkle）	摩擦锥的形状（圆 → \(d\) 边形）	几何离散	\(O(1/d^2)\)	LCP 结构精确（线性）
凸松弛（Anitescu）	约束本身（互补 → 凸不等式 + \(\phi/h\)）	凸化	\(O(h)\) 人工滑移	凸 QP / CCP

冲击模型	定义	约束的量	能量性质
Newton	\(v_n^+=-e\,v_n^-\)	法向速度比	\(e=1\) 保能；\(e<1\) 耗散；多体下可能增能（Kane）
Poisson	\(P_n^{\text{rest}}=e\,P_n^{\text{comp}}\)	冲量比（分压缩段/恢复段）	摩擦多体下比 Newton 更能量一致
Moreau	\(v_n^+=-e\,v_n^-\) 但作为包含 \(-P\in N_{T_C}(v^++ev^-)\)	速度级统一化	自动耗散摩擦，时步法原生
Stronge	能量恢复 \(e_*^2=E^+/E^-\)	能量比	避免 Kane 反例的能量增益

层次	位置更新用	能量行为	刚性	计算量	代表
显式 Euler	旧速度 \(v_\ell\)	注入（发散）	差	最低	几乎无人用于接触
半隐式 Euler	新速度 \(v_{\ell+1}\)	辛（长期有界）	中	低	MuJoCo、Bullet、Drake
全隐式 / θ=1	新状态	耗散 \(O(h^2)\)	好	高（Newton）	Drake SAP midpoint

稳定化	原理	优点	缺点
Baumgarte	人工阻尼 \(\ddot\phi+2\alpha\dot\phi+\beta^2\phi=0\)	简单	参数难调、注入能量、sim2real gap
GGL	位置+速度同时约束	稳健、无可调阻尼	系统更大
坐标投影	每步投影回流形	几何干净、不注入能量	投影开销

符号	含义	首次出现
\(\phi(q)\)	gap 函数（有符号距离），\(\phi\ge 0\) 不穿透	§0
\(M\)	广义质量矩阵（正定 \(M\succ 0\)）	§0
\(H=\nabla\phi\)（或 \(J\)）	接触 Jacobian（gap 梯度）	§0, §6
\(\lambda_n\)	法向接触力（\(\ge 0\)）	§0
\(P\)	接触冲量（\(=\int d\mathbf{R}\)，量纲 \(N\cdot s\)）	§0, §3
\(e\)	恢复系数（restitution，\(\in[0,1]\)）	§0
\(d\mu\)	基测度 \(=dt+\sum_i\delta_{t_i}\)（勒贝格 + 原子）	§2
\(dv\)	速度的微分测度（含跳跃）	§2
\(d\mathbf{R}\)	接触冲量测度	§2
\(v^+,v^-\)	速度的右/左极限（\(\lim_{s\downarrow t}, \lim_{s\uparrow t}\)）	§2
\(C(q)\)	可行集 \(\{q:\phi(q)\ge 0\}\)	§2
\(T_C(q)\)	切锥（tangent cone）	§2
\(N_K(x)\)	凸锥 \(K\) 在 \(x\) 的法锥（normal cone）	§2
\(\theta\)	θ-method 加权系数（\(\in[0,1]\)）	§3
\(W=HM^{-1}H^\top\)	Delassus 算子（接触空间逆惯量）	§3
\(U\)	接触点相对速度（含恢复 \(v^++ev^-\) 的法向投影）	§3
\(v_{\text{free}}\)	自由速度（无接触时一步速度）	§3
\(D=[d_1,\dots,d_d]\)	切向方向基（多面体摩擦锥）	§4
\(\boldsymbol\beta\)	各方向摩擦冲量（\(\ge 0\)）	§4
\(\sigma\)	滑移率乘子（slip rate multiplier，\(\ge 0\)）	§4
\(E=[1,\dots,1]^\top\)	滑移率耦合向量	§4
\(M_{\text{LCP}}\)	Stewart-Trinkle LCP 矩阵	§4
\(\mu\)	摩擦系数	§4
\(\mathcal{K},\mathcal{K}^*\)	摩擦锥及其对偶锥	§5, §6
\(N=D^\top M^{-1}D\)	Anitescu 凸 QP 的 Hessian（= Delassus）	§5
\(\boldsymbol\gamma\)	接触冲量向量（凸 QP 未知量）	§5
\(A=JM^{-1}J^\top\)	MuJoCo 对偶 Hessian（= Delassus 连续版）	§6
\(a_{\text{ref}}\)	参考加速度（reference acceleration）	§6
\(R\)	软约束正则化矩阵（impedance）	§6
`solref`	MuJoCo 约束时间常数/阻尼比参数	§6
`solimp`	MuJoCo 阻抗（impedance）参数	§6
\(\alpha,\beta\)	Baumgarte 稳定化的阻尼/频率参数	§9

定理/公式	一句话说明	对应节	出处
Moreau MDI	\(-M\,dv+F\,dt\in N_{T_C(q)}(v^+)\,d\mu\)；持续接触 + 冲击统一公式	§2	Moreau 1988 CISM 302
Moreau-Jean θ-scheme	速度级离散 + 冲击律嵌入 argument \(v_{\ell+1}+ev_\ell\)；每步解 Delassus LCP	§3	Moreau/Jean CMAME 177 (1999)
θ 能量估计	\(\theta=\tfrac12\) 能量中性；\(\theta=1\) 耗散 \(O(h^2)\)	§3	经典 θ-method 分析
Delassus 半正定	\(W=HM^{-1}H^\top\succeq 0\)（因 \(M\succ 0\)）⇒ 法向 LCP 可解	§3	LCP 理论
Stewart-Trinkle LCP	\(d\)-方向多面体锥 + 滑移率乘子复现最大耗散；矩阵见 §4	§4	Stewart-Trinkle IJNME 39 (1996) 2673
Anitescu-Potra copositive-plus	\(M_{\text{LCP}}\) copositive-plus ⇒ Cottle-Pang-Stone ⇒ 每步可解	§4	Anitescu-Potra Nonlinear Dynam. 14 (1997) 231
Stewart 1998 收敛 + Painlevé	\(h\to 0\) 弱 * 收敛到 MDI 解；Painlevé 含 atom 冲量解（impact without collision）	§4	Stewart ARMA 145 (1998) 215, DOI 10.1007/s002050050129
Anitescu 2006 凸松弛收敛	凸 QP 松弛（\(\phi/h\) 稳定项）\(\epsilon\to 0\) 收敛到 DVI/MDI 解	§5	Anitescu Math. Program. 105 (2006) 113
CCP 形式	\(\mathcal{K}\ni\gamma\perp(N\gamma+r)\in\mathcal{K}^*\)；投影梯度可解	§5	Tasora-Anitescu 2010/2011
MuJoCo Gauss 凸优化	\(\min\tfrac12\\|x-M^{-1}(\tau-c)\\|_M^2+\dots\)；对偶 Hessian \(A+R\)	§6	MuJoCo 文档；Todorov 2012/2014
Moreau-Jean 弱收敛	prox-regular 约束 + \(\theta\in[\tfrac12,1]\)，\(q_h\to q\) 一致、\(v_h\overset{*}{\rightharpoonup}v\)	§3	Moreau/Jean CMAME 177 (1999)
Dzonou-Monteiro Marques	带 Coulomb 摩擦 Moreau-Jean 弱收敛（de Saxcé 双位锥克服非半连续）	§2, §7	Dzonou-MM Nonlinear Anal. TMA (2007)
Cirak-West 不相容	显式辛格式在接触激活/穿透修正下必破坏辛性	§8	Cirak-West IJNME 64 (2005) 1079
Kane-Marsden-Ortiz-West	变分积分器在弹性接触（\(e=1\)）下保辛的条件	§8	KMOW IJNME 49 (2000) 1295
Castro SAP 全局收敛	SAP 势函数严格凸 + Lipschitz-Hessian ⇒ Newton 全局线性-二次收敛	§10	Castro-Permenter-Han T-RO 39 (2022) 1301

#	症状	可能原因	排查步骤	相关节
1	物体穿过地板/相互嵌入	(a) 步长过大；(b) 速度级格式漂移；(c) 软接触固有穿透；(d) 碰撞检测漏检	① 减小步长 \(h\) 看穿透是否减小（若线性减小 → 离散误差）；② 判断仿真器是软接触（MuJoCo，穿透正常）还是刚性（需投影）；③ 速度级格式加坐标投影；④ 检查碰撞 margin 设置	§0, §6, §9
2	仿真在某时刻"卡死"/步长缩到机器精度	用了事件驱动法 + Zeno 场景（bouncing ball 落定）	① 确认是否事件驱动（root-finding）；② 若是，改用 time-stepping（固定步长不检测事件）；③ time-stepping 会把多次微弹自动吸收	§0, §1
3	能量莫名增加/系统越跑越快	(a) 显式 Euler 注入能量；(b) Baumgarte 稳定化注入能量；(c) 求解器未收敛	① 检查积分器，改半隐式以上；② 弃 Baumgarte，改坐标投影/GGL；③ 增大求解器迭代数；④ 检查 \(e>1\) 误设	§8, §9
4	能量过度损失/运动"黏稠"变慢	(a) \(\theta=1\) 或大数值阻尼；(b) Baumgarte 过阻尼	① θ-scheme 改 \(\theta=0.5\)（能量中性）；② 减小数值阻尼/`solref` 阻尼比；③ 检查接触 \(e\) 是否过小	§3, §8
5	静止堆叠物体高频抖动（jitter）	(a) 步长过大；(b) 刚度过高 vs 步长不匹配；(c) 求解器未收敛；(d) 接触力高频切换	① 增大求解器迭代数（最常见）；② 减小步长或降低刚度；③ 适当增大数值阻尼（\(\theta\to 1\)）；④ warm-start 稳定接触力	§6, §9
6	物体在斜面上本应静止却缓慢下滑	凸松弛的 \(O(h)\) gliding（非 bug，非参数错误）	① 确认是否凸格式（SAP/MuJoCo/Chrono）；② 减小步长 \(h\) 看滑移是否线性减小（确认是 gliding）；③ 接受它或减小 \(h\)；④ 勿靠调大摩擦系数硬修	§5
7	换求解器/换仿真器结果不同	(a) 多体同时冲击 LCP 解非唯一；(b) 不同摩擦锥/凸松弛；(c) 大步长 gliding 差异	① 确认场景是否多体同时冲击（解非唯一是正常的）；② 对照 §10 各仿真器格式差异；③ 要可复现须固定仿真器+版本+求解器+步长	§4, §7, §10
8	接触求解器报"不收敛"/失败	(a) Drake TAMSI 非凸失败；(b) ill-conditioned Delassus；(c) 摩擦系数过大触发病态	① 改用凸求解器（SAP/ADMM）保证全局收敛；② Pinocchio ADMM + proximal 正则化对付 ill-conditioned；③ 检查质量矩阵/接触配置是否病态	§5, §10

30_时步法与数值积分

接触力学 · 专题 3.0：时步法与数值积分——非光滑接触动力学的时间离散¶

前置自测¶

本章目标¶

本章知识导航¶

前置知识桥接¶

如果跳过本章会怎样¶

预计阅读时间¶

§0 为什么"时间离散"是独立的数学难题 ⭐⭐¶

动机：一个看似简单的问题¶

如果用经典 ODE 积分器会怎样¶

Zeno 现象：事件驱动法的"死穴"¶

基础铺垫：gap 函数与接触 Jacobian¶

三重难度叠加¶

本专题的独立性：为什么 Stewart 要专门开一个分支¶

⚠️ 常见陷阱¶

练习¶

§1 两条方法论路线：事件驱动 vs 时步法 ⭐⭐¶

动机：面对同一个接触系统，有两种截然不同的算法哲学¶

历史：两条路线的源头¶

理论：逐维度对比¶

事件驱动法的另一面：Filippov 系统与滑模¶

桥接：机器人工程的现实选择¶

⚠️ 常见陷阱¶

练习¶

§2 Moreau 扫动过程与测度微分包含（MDI）⭐⭐⭐¶

动机：我们需要一个能容纳"跳跃速度 + 测度力"的方程¶

反面：为什么不能用"微分方程 + Dirac 函数"硬凑¶

历史：Moreau 的扫动过程¶

理论：从单边约束到测度微分包含¶

Moreau 的关键代数化：从 \(N_C\) 到 \(N_{T_C}\)¶

速度级 vs 位置级约束：一个贯穿全章的分歧¶

摩擦的引入：从法锥到 de Saxcé 双位锥¶

一个几何算例：扫动过程如何"扫"出接触力¶

对比：扫动过程 vs 罚函数（软接触）的两种世界观¶

⚠️ 常见陷阱¶

练习¶

§3 Moreau-Jean θ-scheme：速度级离散 ⭐⭐⭐¶

动机：把连续 MDI 变成计算机能跑的迭代¶

反面：朴素离散会出什么问题¶

历史：Moreau 与 Jean¶

理论：θ-scheme 的完整推导¶

几何直觉：Delassus 算子是什么¶

能量估计：θ 的物理含义¶

把每步 LCP 看成一个投影：通往 Moreau-Yosida 与 ADMM¶

桥接：Moreau-Jean 在工程中的位置¶

⚠️ 常见陷阱¶

练习¶

§4 Stewart-Trinkle 时步格式：位置级 LCP ⭐⭐⭐¶

动机：工程上最有影响力的离散格式¶

反面：为什么不能直接用二阶锥摩擦¶

历史：Stewart-Trinkle 与 Anitescu-Potra¶

理论：Stewart-Trinkle LCP 的构造¶

一个具体的 LCP 矩阵：2D 方块在地面¶

核心定理：Anitescu-Potra copositive-plus¶

收敛性与 Painlevé 悖论的消解¶

Stewart-Trinkle vs Moreau-Jean：位置级 vs 速度级的再对比¶

⚠️ 常见陷阱¶

练习¶

§5 Anitescu 凸松弛：从 LCP 到凸 QP/CCP ⭐⭐⭐¶

动机：LCP 的两个痛点，凸优化的两个礼物¶

反面：为什么原始问题非凸¶

历史：从 Stewart-Trinkle 到 Anitescu 到 Tasora-Chrono¶

理论：Anitescu 凸松弛的构造¶

锥投影 \(\Pi_{\mathcal{K}}\) 的解析形式与 APGD 迭代¶

凸松弛的代价：人工滑移（gliding / creep）¶

凸松弛与多面体锥：正交的两个近似¶

桥接：凸性是可微接触的前提¶

⚠️ 常见陷阱¶

练习¶

§6 MuJoCo 软接触求解器：连续时间 + 凸优化 ⭐⭐⭐¶

动机：为 RL 与可微仿真而生的另一条路¶

反面：硬接触对 RL/控制的不友好¶

历史：Todorov 的 MuJoCo 与 Gauss 原理¶

理论：MuJoCo 求解的凸优化问题¶

半隐式 Euler 与积分器选项¶

桥接：为什么 MuJoCo 选这条路¶

⚠️ 常见陷阱¶

练习¶

§7 冲击律的离散处理 ⭐⭐¶