多体系统LCP模型转换为QP的严格性分析

1.原始约束与LCP→QP转化的条件

原始的多体接触/摩擦约束满足线性互补问题(LCP)的等价QP条件需要相当严格的假设。 特别是要求系统矩阵对称正定和约束集凸。

对于无摩擦的接触（纯单侧约束），接触力可通过对称的Delassus矩阵 $N = D^T M^{-1} D$ 构造正定二次型；此时LCP可视为某个QP的KKT条件 （线性互补问题（LCP）其实就是某类二次规划（QP）的最优性条件，特别是在约束是盒约束（或不等式）时）， 满足等价转换条件。

然而，引入库仑摩擦后，法向与切向约束耦合导致非线性互补问题（NCP），不再满足直接转化为凸QP的充分条件。 论文指出基于库仑摩擦的DVI离散化会形成一个非凸的互补问题，需要特殊处理。 因此，原始约束并未天然满足将LCP严格等价为QP的条件。 必须通过假设或修改（如对称化、凸化）使问题满足凸QP条件。

正如论文所述，他们通过引入决策(D3)对NCP进行凸化，使之成为圆锥互补问题(CCP)，

这依赖于接触矩阵$N$的对称半正定性并将摩擦锥约束视为凸集Υ_i。 只有经过这样的处理，才能将原始问题转换为在凸约束下求解的QP。

总之，在未作特殊处理前，原始库仑摩擦约束不具备LCP到QP严格等价的充分条件，

需通过对模型加以假设和修正（如将NCP凸化为CCP）才能建立等价的QP模型。

？？？

我的问题是：你的意思是LCP直接进行提取schur补矩阵和右端项转换为QP是不准确的吗？我这个代码有进行凸化转化为CCP吗？

考证后：

所以这个我在求解器设计中，没有直接从 LCP 生硬转 QP，而是通过 projectFriction 做了摩擦锥投影， 这个操作已经是凸化，所以求解的是 CCP（Convex Cone Complementarity Problem）。这种 CCP 是可以等价成二阶锥规划（SOCP）或凸 QP 来解的，所以你的方法是符合 QP 求解假设的（但前提是摩擦锥投影正确实现）。

如果要并行，最好从建模阶段就按 CCP 的接触块结构组织数据，这样能像 Chrono 那样实现并行求解.

2.APGD求解的QP与原始LCP（KKT）的关系

APGD方法实际求解的是由接触互补问题的KKT条件形成的QP。

论文明确指出，经过上述凸化处理后，摩擦接触问题可表述为带圆锥约束的二次优化问题 (Equation 11)，

\[\min_{\lambda \in \mathcal{K}} \quad f(\lambda) = \frac{1}{2} \lambda^T Z \lambda + \lambda^T d\]

其一阶最优性条件正是圆锥互补形式。

也就是说，该QP的KKT条件等价于求解松弛后接触问题的互补条件(CCP)。

APGD（加速投影梯度下降）利用Nesterov方法，在每个仿真步求解这个带摩擦锥约束的QP， 以求出法向和摩擦方向的约束力。

实际上，APGD方法等价于求解上述QP的对偶形式： 正如论文所述，通过放松原始问题得到的凸优化问题的对偶形式“实质上就是本文采用的方法”。

因此，APGD并非直接求原始LCP，而是求解来自（松弛后）LCP的KKT条件对应的QP。 这一QP的目标函数一般是$ \frac{1}{2}\gamma^T N \gamma + r^T \gamma$， 其中$\gamma$为未知冲量/约束力，$N$为接触的Schur补矩阵，$r$为外力和上一时刻状态形成的项。

QP的约束$\gamma_i\in Υ_i$即各接触力需位于库仑摩擦圆锥内。 APGD通过迭代$\gamma$使其满足KKT：梯度步相当于计算$N\gamma + r$，投影步将$\gamma$投影回摩擦锥可行域。

综上，APGD求解的正是松弛后接触问题KKT条件对应的QP，从数值角度保证了等价于该松弛互补问题的解。 需要注意，这里的KKT对应的是经过凸化/松弛处理后的问题，而非未经松弛的原始NCP模型。

3.摩擦锥投影、互补松弛与伪变量的近似处理 （后续深入研究）

在该方法中，库仑摩擦约束通过圆锥投影精确处理，而法向-切向互补条件则经过松弛处理，引入了一定近似。 具体而言：

相关知识点

1.系统矩阵对称正定

“系统矩阵对称正定”和“约束集凸”是将**线性互补问题（LCP）转换为二次规划问题（QP）**的两个重要数学条件。下面我详细讲解它们各自的含义，以及为什么这两者对从 LCP 到 QP 的等价转化至关重要。

系统矩阵对称正定（Symmetric Positive Definite, SPD）定义：

一个矩阵 $M \in \mathbb{R}^{n \times n}$ 被称为对称正定矩阵，当且仅当满足以下两个条件： * 对称性：$M = M^T$ * 正定性：对任意非零向量 $z \in \mathbb{R}^n$，有：

\[z^T M z > 0\]

🧩 在多体动力学中，它是谁？

在APGD 和 Dantzig 所求解的 LCP 形式如下：

\[w = M z + q,\quad w \geq 0,\; z \geq 0,\; w^T z = 0\]

其中：

所以如果：

则 $J M^{-1} J^T$ 就是对称正定的。

✅ 重要性：

如果系统矩阵是对称正定的，那么原始的 LCP 就可以等价地转化为一个凸二次规划（QP）问题：

\[\min_{z \geq 0} \;\; \tfrac{1}{2} z^T M z + q^T z\]

这一步的等价性是精确的（数学上严格成立）。

2.约束集凸

约束集凸（Feasible Set is Convex）

✅ 定义：

所谓“约束集凸”是指 QP 的可行解集合（也就是变量 $z$ 的取值范围）必须构成一个凸集。即：

🧩 在本问题中，它是谁？

对于无摩擦 LCP：

对于含摩擦的 LCP：

✅ 为什么它重要？

二次规划只有在目标函数是凸函数（由 SPD 保证）且约束集是凸集时，才能保证：

🔁 总结一下

条件 数学意义 对 QP 转化影响
系统矩阵对称正定 $M \succ 0$ 目标函数是凸函数 保证凸优化问题
约束集是凸集 $z \in \mathcal{C}$ 可行域是凸 保证最优解存在与唯一
二者共同作用 保证 LCP ⇄ QP 完全等价 可安全使用 APGD 求解

3.Delassus矩阵

\[J \, M^{-1} \, J^T \, \lambda = b\]

这里：

含义

4.NCP

这是因为 摩擦约束本质上是非线性的，而且它会把法向约束和切向约束耦合在一起，从而让原本的 LCP（线性互补问题）变成 NCP（非线性互补问题）。

我给你拆开讲一下原因：

1️⃣ 先看无摩擦的情况（纯 LCP）

在无摩擦接触中：

力学离散后，KKT 约束是线性的：

\[w_n = N_{nn} \lambda_n + r_n \ge 0\]

这就是一个 LCP，因为 $w_n$ 是 $\lambda_n$ 的线性函数。

2️⃣ 加入摩擦（法向和切向耦合）

当你有摩擦时，切向约束来自 库仑摩擦锥：

\[\|\lambda_t\| \le \mu \lambda_n\]

其中：

这个约束 不是线性：

3️⃣ 为什么导致 NCP

在动力学离散方程中，法向和切向力的计算公式是：

\[w_n = N_{nn} \lambda_n + N_{nt} \lambda_t + r_n\] \[w_t = N_{tn} \lambda_n + N_{tt} \lambda_t + r_t\]

摩擦锥条件会把 $\lambda_t$ 和 $\lambda_n$ 绑在一起：

因为这个关系是 乘法形式 + 范数非线性，所以方程系统变成了 非线性互补问题（NCP）。

4️⃣ 如何处理 NCP

在数值计算中，常见的处理方式有：

  1. 线性化摩擦锥（变成 LCP）

    • 用多边形近似圆锥（切向摩擦力方向离散成若干条边）
    • 这样 $|\lambda_t| \le \mu \lambda_n$ 变成若干线性不等式

  2. 直接解 NCP

    • 用非线性互补求解器（如 PATH、semismooth Newton）
    • 代价高，但更准确

  3. 投影法（PGS/APGD 等）

    • 在每一步梯度下降后，把解投影到摩擦锥上
    • 这是 Chrono、你代码里的 projectFriction() 这种方法

💡 一句话总结： 摩擦约束中的 $|\lambda_t| \le \mu \lambda_n$ 是一个非线性几何约束，它把法向力和切向力耦合起来，使得系统从 LCP 变成 NCP。

Reference