Elastic Locomotion with Mixed Second-order Differentiation

一句话总结

提出一种把反向自动微分与复步长有限差分融合的”混合二阶微分”工具，能高效计算带接触屏障的逆仿真损失函数的 Hessian，从而用牛顿法直接求解软体角色的肌肉激活，让用户通过高层运动描述驱动弹性体完成爬行、跳跃、翻转等生动的运动控制。

研究背景

很多角色依赖刚性骨架来产生协调运动，但自然界中大量生物没有内骨骼，它们靠肌肉收缩或舒张来与环境交互。给这类软体角色做运动控制天然困难：软体表面柔韧且接触面积大，会产生高维接触问题，涉及一组不等式约束。接触模型虽已被充分研究（可简化为线性互补规划），但其组合本质使问题属于 NP 难，实际计算仍很棘手。

同时，真实材料大多是非线性的。仅仅给定外力做正向仿真就已经很昂贵，因为成千上万个自由度是双向耦合的。如果用户想直接、直观地控制仿真结果（例如让身体沿指定路径运动），目标函数与控制参数之间只能通过正向过程隐式关联，没有闭式表达式，还夹杂接触约束。这种高非线性叠加高维度的组合，超出了现有仿真技术的能力。

本文受增量势接触（IPC）启发，用对数屏障惩罚把不等式约束重写为无约束优化，从而回避互补规划的组合爆炸，能处理更高维的接触。但 IPC 会进一步加剧正向仿真的非线性，导致优化地形崎岖不平、充满鞍点与局部极小，只依赖梯度信息的现有可微/逆仿真技术几乎无能为力。作者据此提出可高效评估二阶导数的混合微分，让牛顿法这类高阶优化工具得以直接使用。

方法

整体框架

弹性运动被建模为一个标准的逆仿真问题：给定嵌入了可收缩/可伸展肌肉纤维的软体，在用户指定的高层运动描述（如身体轨迹）下，反解出驱动身体自然协调运动所需的肌肉激活力。方法把高维时空优化分解为逐帧的非线性规划：每一帧用一个混合微分模态高效返回损失函数的一阶和二阶梯度信息，且该微分与具体采用的仿真算法、接触交互数量无关。二阶下降方向对这种强非线性逆问题至关重要。

flowchart LR
    A[用户高层运动输入] --> B[运动学损失 + 正则项]
    C[正向仿真<br/>非线性 FEM + 隐式欧拉 + IPC 屏障] --> B
    B --> D[反向 AD<br/>求梯度]
    D --> E[CSFD 提升<br/>沿虚轴扰动]
    E --> F[损失 Hessian]
    F --> G[线搜索牛顿法]
    G --> H[最优肌肉激活]
    H --> I[进入下一帧]
    I --> C

正向动力学用非线性 FEM 加隐式欧拉写成力平衡形式：

\[\boldsymbol{x}_{t+1} = \boldsymbol{x}_t + \Delta t\, \boldsymbol{v}_{t+1}, \quad \boldsymbol{v}_{t+1} = \boldsymbol{v}_t + \Delta t\, \boldsymbol{M}^{-1}(\boldsymbol{f}_{int} + \boldsymbol{f}_{ext} + \boldsymbol{f}_d + \boldsymbol{f}_c + \boldsymbol{f}_m).\]

其中弹性力采用 stable Neo-Hookean 能量，外力为重力，阻尼用 Rayleigh 阻尼；接触力遵循 IPC，当表面三角形足够接近碰撞体时激活对数屏障势能 \(B(d, \hat{d})\)，\(\hat{d}\) 是约束容差（实现取 \(\hat{d}=1\mathrm{E}{-3}\)），\(\kappa\) 控制屏障初始刚度，摩擦力则用上一迭代碰撞力信息做滞后（lagged）计算。

肌肉纤维建模为带 \(M\) 段的多段折线，每段是一个只能沿自身方向收缩或伸展、不能弯曲的弹簧，施加激活量 \(a\)。一段激活会影响附近多个单元，第 \(i\) 个单元与第 \(j\) 段肌肉的影响权重基于二者在网格上的测地距离 \(g_{ij}\) 用高斯核给出：

\[w_{ij} = \exp\!\left(-\frac{g_{ij}^2}{c^2}\right).\]

权重只依赖静止形状，可预计算。所有肌肉力用一个位姿相关的激活矩阵编码为 \(\boldsymbol{f}_m = \boldsymbol{A}(\boldsymbol{x})\,\boldsymbol{a}\)，其中 \(\boldsymbol{a}\) 是所有肌肉段的激活向量。作者假定肌肉布置是给定的（纵向、环形、径向、斜向/螺旋等不同排布支持缩短、弯曲、变细拉长、扭转等动作）。

优化目标写为运动学目标项与正则项的加权和：

\[\arg\min_{\boldsymbol{a}} L(\boldsymbol{x}, \boldsymbol{a}) = \sum_i w_i G_i(\boldsymbol{x}_{t+1}, \boldsymbol{a}) + \sum_j \lambda_j R_j(\boldsymbol{x}_{t+1}, \boldsymbol{a}).\]

由于 \(\boldsymbol{x}_{t+1}\) 通过运动方程隐式依赖 \(\boldsymbol{a}\) 并涉及动态接触与摩擦，问题强非线性。求解时先做几次梯度下降提供稳健初值（对初始选择不敏感），再切换到牛顿法利用其在局部极小附近的二阶收敛；上一帧的激活 \(\boldsymbol{a}_t\) 作为初始猜测。

关键设计：混合二阶微分（CSFD-AD）

求解流程本身是标准的线搜索非线性优化，真正的难点在于损失 Hessian 的评估。作者结合两种微分模态：

• 反向 AD：损失是”多输入到单标量”的映射，反向 AD 通过前向存值、反向传播伴随量（adjoint）高效得到梯度。但若单靠反向 AD 求 Hessian，需要两遍 AD 并保存整张计算图与所有中间量，在这种屏障在环 FEM 过程里内存开销不可承受，且反复套用 AD 求高阶导数在数值上不稳定。

• 复步长有限差分（CSFD）：把实值函数提升到复域，沿虚轴施加扰动来估计导数，\(f'(x_0) \approx \dfrac{\mathrm{Im}\, f^{*}(x_0 + h i)}{h}\)。分子没有相减项，因此不受传统有限差分的相减相消影响，扰动 \(h\) 可取到极小（如 \(1\times 10^{-16}\)），精度可达机器精度，媲美解析导数。

核心思路是把反向 AD 当作一个”把输入映射到一阶导数”的通用黑盒函数，再用 CSFD 去扰动它，一次求值就自然得到二阶导数，作者称之为 CSFD-AD。它沿用 AD 算法，但把每个中间变量与伴随量都提升为复数。Hessian 的第 \(k\) 列由

\[[\nabla^2 f]_k \approx \frac{\mathrm{Im}\big(\nabla f(\boldsymbol{x} + h i\cdot \boldsymbol{e}_k)\big)}{h}\]

给出，因此只需 \(M\) 次扰动（而非广义高阶 CSFD 的 \(M^2\) 次）即可组装整个 Hessian，且这 \(M\) 次扰动可并行。

为让 CSFD-AD 在真实计算流程中可用，作者仔细处理了运算符与函数的复数重载与可微性：关系运算（如 \(>\)）只比较实部，保证与原算法分支一致；条件分支因扰动沿虚轴（与实轴正交）总是返回 \(\mathrm{Re}(x_0^{*})\) 所在分支的导数，无需特殊安全处理；大多数初等函数在复域解析（满足 Cauchy–Riemann 条件）可直接提升，而绝对值 abs 因标准复数模长不解析，需重新推导其解析形式与伴随量。此外还为矩阵/向量运算做了专门实现来避免计算图节点膨胀（如 \(\|\boldsymbol{x}\|^2\) 用专用伴随 \(\bar{\boldsymbol{x}} = 2y\boldsymbol{x}\)），并支持 trace、determinant、inverse、奇异值等线性代数函数。SVD 因复矩阵奇异值恒为实数而不解析，作者转而对迭代式数值 SVD（隐式移位对称 QR SVD）逐操作做复数重载；稀疏线性求解 \(\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}\) 则借鉴伴随法，前向存分解、反向用同一分解求 \(\bar{\boldsymbol{b}} = \boldsymbol{A}^{-\top}\bar{\boldsymbol{x}}\)，避开对矩阵求逆求伴随导致的计算图指数膨胀，也不依赖预设的迭代次数（这对接触迭代数剧烈波动的场景很关键）。

运动控制器

框架支持一系列直观的高层控制器，用户无需手动推导闭式导数：线性运动控制器约束关键点（质心 COM、局部质心、相对位置）的位置、速度、加速度；角运动控制器通过操纵角动量变化 \(G_{angular} = \|\dot{\boldsymbol{L}}_{angular}(\boldsymbol{x}) - \dot{\boldsymbol{L}}^{\star}\|^2\) 控制旋转与扭转，或通过控制转动惯量（MOI）\(G_{MOI}\) 让角色收缩身体从而在空中转得更快（跳水运动员常用技巧）；还能控制弹性能变化 \(G_{elastic}\)（如通过蹦床势能控制起跳后的动能高度）以及投影底面积 \(G_{proj}\)（落地时张开底面保持平衡）。正则项包括约束总肌肉能量 \(R_{energy}=\frac{1}{2}k\|\boldsymbol{a}\|^2\) 以及惩罚激活变化率超过阈值的 \(R_{change}\)，用来保证目标函数 Hessian 正定、增强优化稳定性（类似 LM 方法）。

实验结果

实现基于 Stan Math 库（反向 AD）与 Eigen 库（复数线性代数），用 Intel TBB 并行计算 Hessian，运行于 Intel i9-13900KF、128GB 内存的台式机，所有实验时间步 \(\Delta t = 1/40\) 秒。下表汇报各示例的规模与耗时（秒/时间步）：

软体	单元数	激活维度	最大接触数	仿真	梯度	Hessian	总求解
Caterpillar	2.6K	62	140	0.8	0.5	5.2	15.4
Starfish	2.8K	66	192	0.9	0.5	6.2	20.3
Chess	8.4K	72	224	1.4	1.0	12.8	36.6
Lamp	5.2K	68	166	1.1	0.9	8.8	19.5
Stool	6.8K	40	246	1.0	0.6	4.8	21.6
“T”	1.2K	51	25	0.2	0.2	3.4	5.5
Puffer fish	27.6K	32	171	4.7	2.8	27.1	68.6

与优化器对比：在”投篮”实验中，一个 I 形软体把球推入篮筐，仅用梯度信息（迭代 1000 次）时球无法沿预定轨迹运动、甚至到不了篮筐；本文的二阶优化让球准确穿过篮圈；即便只在释放球前最后一步切到梯度下降，轨迹仍不准（球从筐沿弹回）。在多种损失组合的代表性帧上，与梯度下降及 L-BFGS（同一 warm start 后比较）相比，本文方法在所有例子中都表现出强二阶收敛；当帧内摩擦接触很多时收敛差距被拉大，而在无接触场景下 L-BFGS 也能给出不错收敛，梯度下降则始终不理想。

与基于 LCP 的接触求解对比：Tan 等人 2012 的 QPCC 求解器每个接触点有 10 种状态、最小值需在 \(10^N\) 的解空间中搜索，最坏情况要穷举，对中等分辨率软体不可行（其实现只支持四个接触面片，牺牲了精度）。在有 140 个接触点、接触状态逐帧变化的软糖毛毛虫爬行示例中，QPCC 无法在合理时间内得到可用解、毛毛虫几乎不动；本文借助 IPC 与线搜索牛顿对接触点数量不敏感，能稳健前进。

更多结果覆盖：爬行（海星、模拟蠕动波的软糖毛毛虫，含纵向与环形肌肉的收缩舒张）、跳跃（物理化的软体 Luxo Jr. 台灯、连续跳跃的软体象棋，起跳/空中/落地分阶段用不同控制器）、蹦床背翻（台灯与可变形蹦床之间的高维复杂接触，用弹性势能目标控制每次跳跃高度）、行走（驮书的四足软体凳子，处理凳-地与凳-书接触）、滚动与旋转（经典 T 形角色跳落旋转，可复现 Tan 等人 2012 的结果；带四个环形肌肉的刺豚鱼左右翻滚、跳上平台取星，QPCC 无法生成此运动）。

亮点与局限

亮点：

• 提出通用的混合二阶微分 CSFD-AD，把黑盒反向 AD 用 CSFD 沿虚轴提升，一次求值得到二阶导，仅需 \(M\) 次可并行扰动即可组装 Hessian，避免了双遍 AD 的内存爆炸与数值不稳定。
• 使强非线性、接触在环的逆仿真第一次能直接用牛顿法求解，二阶收敛明显优于梯度下降与 L-BFGS，尤其在大量摩擦接触时优势显著。
• 与具体仿真算法、接触数量解耦，且对控制器高度友好——用户可用复杂目标函数而无需手推导数，展示了从爬行、跳跃、翻转到行走、滚动的丰富运动。
• 工程上完整处理了关系运算、条件分支、abs、SVD、稀疏线性求解等在复域的可微性与实现细节。

局限：

• 组装大规模 Hessian 昂贵，复杂度至少 \(O(N^2)\)。
• 作者指出可考虑把 CSFD-AD 与不需要显式 Hessian 的优化算法（如 L-BFGS）结合。
• 因损失是多对一映射而选用反向 AD；若是低维输入映射到高维输出，则更适合用前向 AD 版本的 CSFD-AD。

延伸思考

该方法把”解析微分负责梯度、数值微分负责升阶”这一分工做成了通用工具，其价值可能超出软体运动控制本身——凡是需要在强非线性、约束在环的仿真流程里拿到可靠 Hessian 的逆问题（如机器人控制、可制造性设计、参数辨识）都可能受益，作者也明确提到方法与形变仿真细节无本质绑定。值得关注的开放问题是如何降低 Hessian 组装的 \(O(N^2)\) 代价：是否能用无矩阵（matrix-free）的 Hessian-向量积配合共轭梯度型牛顿-CG，在保留二阶信息的同时避免显式组装，从而把方法推向更高分辨率的角色与更复杂的多体接触场景。