摘要

我们研究两类求解器，用于计算由大量刚体构成、并通过带摩擦的接触相互作用的系统的时间演化。接触通过互补条件建模；摩擦则表述为一个变分问题。系统动力学由一组微分代数方程 (DAE) 与微分变分不等式 (DVI) 耦合构成。在时间离散之后，将施加在速度层面的互补条件适当放松，从而得到一个锥互补问题 (CCP)。CCP 的求解成为仿真的瓶颈，并通过最小化一个带锥约束的等价二次优化问题来实现。本文针对这一受约束优化问题考察两大类求解器。

投影 Gauss–Jacobi (PGJ)、投影 Gauss–Seidel (PGS) 以及加速投影梯度下降 (APGD) 方法代表了第一类求解器。它们属于一阶方法，只利用目标函数值及其梯度信息。第二类求解器由对称锥内点 (SCIP) 方法和原–对偶内点 (PDIP) 方法所代表。这类二阶方法依靠牛顿步来确定下降方向，并通过线搜索来计算步长。所有五种方法都利用图形处理单元 (GPU) 上的并行计算；其中牛顿步使用稀疏并行 GPU 线性求解器实现。

我们通过两类数值算例——填充 (filling) 与犁土/切削 (drafting)——来评估这五种求解策略在收敛速度、精度和计算成本方面的表现。为保证可比性，所有数值实验均在同一套开源代码的基础上进行，该代码经过修改以集成上述五种方法。由 Elsevier B.V. 出版。

1. 引言

本研究试图回答下面这个问题：

在用微分变分方法来表述刚体摩擦接触动力学问题时，采用二阶优化方法（即使用 Hessian 信息的方法）是否比一阶方法（只使用梯度信息的方法）更高效？

这个问题在多体动力学问题中非常关键，特别是刚体之间通过摩擦和接触力相互作用的场景。这类多体系统的代表之一就是颗粒物料体系。正如文献 [1] 中指出的，工业中超过 50% 的被加工材料都是颗粒形式，理解它们的动力学行为，对非常广泛的实际应用都至关重要，比如：增材制造、土力学（terramechanics）、纳米颗粒自组装、雪崩动力学、复合材料、火山碎屑流、行星和小行星的形成、陨石撞击坑的形成，以及制药、食品加工、农业、采矿等诸多行业。

在本文中，多体动力学问题是置于一个微分变分框架下进行建模的。一方面，这种做法的好处在于：它有坚实的解析理论基础，并且能够处理诸如 Painlevé 悖论 [2] 这样的经典难题。另一方面，这种建模方式的数值求解过程要求在每一个时间步上解一个规模巨大的优化问题，其解给出了每个接触点处起作用的摩擦力和接触力。 当系统中存在大量活跃接触（例如在颗粒动力学中，这是常态）时，这个优化问题就会成为整个数值求解过程中的主要瓶颈。

我们想要弄清的是：采用二阶优化算法是否能加速整个数值求解过程？

直观上看，可以预期：如果在求解过程中引入更多的关于问题本身的信息——也就是说，把 Hessian 融入求解过程——就能够得到需要更少迭代次数的算法来求解上述优化问题。

然而，二阶方法同时也需要解线性方程组，而在多体动力学问题中，这些线性系统往往非常庞大。问题的关键在于：迭代次数的减少，能否抵消每次迭代计算成本增加所带来的开销？

研究表明，对于相对较小规模的问题，当二阶优化方法配合快速并行稀疏线性求解器（用于预条件）时，这种基于二阶方法的方案往往是更优选择。

但是，随着多体动力学问题规模的增长，在线性代数上的计算时间预期会逐渐占据主导地位，并最终抵消甚至超过由迭代次数减少带来的好处。

上述这些定性的论断，在后文会通过定量分析加以验证。这需要一些预备知识，而这些预备内容会在第 2.1 节中给出。第 2.1 节中，我们将建立这样一个微分变分问题：其解描述了由刚体组成、并通过摩擦与接触相互作用的多体系统的时间演化。

用于求解摩擦力和接触力的一阶和二阶优化方法会在第 3 节中介绍。第 4 节中将对这些方法进行规模（缩放）分析，并结合一组在数值求解上已知较具挑战性的基准算例进行评估。

第 5 节给出了验证研究的结果，我们在其中将数值结果与实验数据进行了对比。第 6 节对全文进行总结并归纳主要经验和结论。最后，附录中汇总了一些推导，这些推导对于论证的完整性来说很重要，但在全文“叙事主线”中相对次要。

2. 预备知识：带摩擦与接触的刚体动力学

2.1 运动方程的建立

一组数量为 ($n_b$) 的刚体通过摩擦和接触相互作用，其随时间的演化在本文中用与每个刚体 (j)（其中 (1 \le j \le n_b)）相关联的笛卡尔坐标来描述。 广义坐标数组

[ \mathbf{q} = [\mathbf{r}1^T, \boldsymbol{\varepsilon}1^T, \ldots, \mathbf{r}{n_b}^T, \boldsymbol{\varepsilon}{n_b}^T]^T \in \mathbb{R}^{n_r n_b}, ]

以及它的时间导数

[ \dot{\mathbf{q}} = [\dot{\mathbf{r}}1^T, \dot{\boldsymbol{\varepsilon}}1^T, \ldots, \dot{\mathbf{r}}{n_b}^T, \dot{\boldsymbol{\varepsilon}}{n_b}^T]^T \in \mathbb{R}^{n_r n_b}, ]

用来表示系统的状态。其中，对刚体 (j) 而言，(\mathbf{r}_j) 和 (\boldsymbol{\varepsilon}_j) 分别是质心的绝对位置和该刚体的欧拉参数（姿态）。

欧拉参数的时间导数 (\dot{\boldsymbol{\varepsilon}}) 可以被另一组未知量所替代；也就是说，可以使用在刚体局部坐标中表示的角速度 (\bar{\boldsymbol{\omega}})。这样可以构造广义速度向量

[ \mathbf{v} = [\dot{\mathbf{r}}1^T, \bar{\boldsymbol{\omega}}1^T, \ldots, \dot{\mathbf{r}}{n_b}^T, \bar{\boldsymbol{\omega}}{n_b}^T]^T, ]

它可以通过文献 [3] 中给出的映射关系，和 (\dot{\mathbf{q}}) 互相对应：

[ \dot{\mathbf{q}} = \mathbf{T}(\mathbf{q}),\mathbf{v}. \tag{1} ]

对应 LaTeX： \dot{\mathbf{q}} = \mathbf{T}(\mathbf{q})\,\mathbf{v}.

运动方程可以写成矩阵形式：

[ \mathbf{M},\dot{\mathbf{v}} = \mathbf{f}(t,\mathbf{q},\mathbf{v}), \tag{2} ]

其中 (\mathbf{M} \in \mathbb{R}^{6 n_b \times 6 n_b}) 是广义质量矩阵，(\mathbf{f}(t,\mathbf{q},\mathbf{v})) 是施加到系统上的广义力向量 ([3,4])。

LaTeX： \mathbf{M}\,\dot{\mathbf{v}} = \mathbf{f}(t,\mathbf{q},\mathbf{v})

系统中存在双边约束（bilateral constraints），例如球铰链、转动关节、移动关节等等，这些约束限制了两个刚体之间的相对运动。从数学上看，每一个这样的物理关节都会引入一组完整约束（holonomic）代数方程：

[ i \in \mathcal{B}:\quad g_i(\mathbf{q}, t)=0, \tag{3} ]

其中 (\mathcal{B}) 是所有双边约束的集合。

这些物理关节在速度层面同样会产生运动学约束，这可以通过对式 (3) 中的位置约束对时间求导得到：

[ \dot{g}_i(\mathbf{q}, t) = \bar{\mathbf{G}}_i(\mathbf{q}),\mathbf{v} + \frac{\partial g_i}{\partial t} = 0. ]

要满足这些约束，就需要引入一个约束反力（constraint reaction force）

[ \bar{\mathbf{G}}i^T(\mathbf{q}),\bar{\boldsymbol{\gamma}}{i,b}, ]

它被加到式 (2) 的右端项中。该反力依赖于一个未知的拉格朗日乘子 (\bar{\boldsymbol{\gamma}}_{i,b}) ([3,4])。

相关 LaTeX：

  1. \dot{g}_i(\mathbf{q}, t) = \bar{\mathbf{G}}_i(\mathbf{q})\,\mathbf{v} + \frac{\partial g_i}{\partial t} = 0
  2. \bar{\mathbf{G}}_i^T(\mathbf{q})\,\bar{\boldsymbol{\gamma}}_{i,b}

顺带帮你理一理这些符号在干嘛（方便你写 APGD / SOCCP）

按你的背景我用“数值仿真/求解器视角”再解释一遍：

  1. 状态变量层次

    • (\mathbf{q})：广义坐标，堆的是每个刚体的平移 (\mathbf{r}_j) 和姿态参数 (\boldsymbol{\varepsilon}_j)（这里用的是四元数/欧拉参数）。
    • (\dot{\mathbf{q}})：它的时间导数。

    • 但做刚体动力学时，姿态部分其实更喜欢用 局部角速度 (\bar{\boldsymbol{\omega}}_j)，于是定义了 [ \mathbf{v} = [\dot{\mathbf{r}}_1^T, \bar{\boldsymbol{\omega}}_1^T, \ldots]^T ] 然后通过一个几何矩阵 (\mathbf{T}(\mathbf{q})) 把 (\mathbf{v}) 转成 (\dot{\mathbf{q}})：

      [ \dot{\mathbf{q}} = \mathbf{T}(\mathbf{q}) \mathbf{v}. ]

      你以后在 Chrono/Verosim 里看到的各种“Gamma”、“Gyration matrix”其实都和这个映射有亲戚关系。

  2. 动力学层次

    • (\mathbf{M})：把所有刚体的质量、惯量拼成一个大的块对角矩阵，所以尺寸是 (6n_b\times 6n_b)（每个刚体 3 个平移 + 3 个转动自由度）。
    • (\mathbf{f}(t,\mathbf{q},\mathbf{v}))：所有外力/广义力（重力、驱动力、接触力、阻尼……）做完映射之后的统称。
    • 不带约束时就是经典的 [ \mathbf{M}\dot{\mathbf{v}}=\mathbf{f}. ]

  3. 几何约束 / 关节

    • 物理关节（球铰、转动、平动……）在位置层面写成 [ g_i(\mathbf{q},t)=0, ] 这是一个标量或小向量方程，典型 holonomic 约束。
    • 对时间求导，得到速度层面的约束： [ \dot{g}_i = \bar{\mathbf{G}}_i(\mathbf{q})\mathbf{v} + \frac{\partial g_i}{\partial t}=0. ] 这里 (\bar{\mathbf{G}}_i) 就是你在 LCP/SOCCP 里会看到的 Jacobian 行。在 Chrono 的记号里大概对应于 (C_q) 那一行。

  4. 拉格朗日乘子和反力

    • 为了强制 (\dot{g}i=0)，对每个约束引入拉格朗日乘子 (\bar{\boldsymbol{\gamma}}{i,b})（下标 b=bilateral）。
    • 广义力里就多了一项 [ \mathbf{f}{\text{constraint}} = \sum{i\in\mathcal{B}} \bar{\mathbf{G}}i^T(\mathbf{q}),\bar{\boldsymbol{\gamma}}{i,b}. ]
    • 从“数值求解器视角”，(\bar{\boldsymbol{\gamma}}_{i,b}) 就是你在 LCP/CCP 中求的 λ / γ。
    • 后面再加上接触摩擦，就会形成一个大的 KKT / VI / SOCCP 系统，然后才轮到 APGD、PGS、PCG 这些朋友登场。

你可以把这一节当成：

“如何从几何约束 + 刚体质量/惯量，把问题写成一个带拉格朗日乘子的矩阵方程，为后面的补充问题（LCP/CCP/SOCCP）铺路。”

后面你读到构造 Schur complement (N) 或 (Z) 的时候，这些 (\mathbf{M})、(\bar{\mathbf{G}}i)、(\bar{\boldsymbol{\gamma}}{i,b}) 都会直接进到你现在在写的 APGD 求解器里。