这是“GNN × XPBD 混合布料”系列的第一篇。后续文章会让图网络预测布料状态或修正量，但在此之前，必须先把物理解算器说清楚：网络究竟在预测什么，哪些量仍由约束保证，训练时的时间步又为什么不能随便改变。

本文用两种标签区分信息来源：论文结论表示原论文明确推导或实验支持的结论；工程建议表示为了实现、调试和扩展混合系统而给出的实践取舍，不能冒充论文结论。

1. 从半隐式积分开始

把布料离散成粒子位置 $\mathbf{x}_i$、速度 $\mathbf{v}_i$、质量 $m_i$，并记逆质量 $w_i=1/m_i$。固定点直接令 $w_i=0$。对长度为 $h$ 的一个时间步，常见的半隐式 Euler 预测是：

$$\mathbf{v}_i^*=\mathbf{v}_i^n+h\,w_i\mathbf{f}_i^{\mathrm{ext}}, \qquad \widetilde{\mathbf{x}}_i=\mathbf{x}_i^n+h\mathbf{v}_i^*.$$

如果到此为止，弹性力也需要显式积分；刚度一高，允许的 $h$ 就会迅速变小。PBD 改变了问题：先用外力预测位置，再直接把预测位置投影回满足约束的区域，最后由位置差重建速度：

$$\mathbf{v}_i^{n+1}=\frac{\mathbf{x}_i^{n+1}-\mathbf{x}_i^n}{h}.$$

因此，“位置修正”并非只影响这一帧的形状，它还会通过重建速度影响下一步运动。

论文结论： Müller 等人的 PBD 原论文把一般约束统一到位置层，并强调位置投影的可控性以及处理碰撞穿透的便利性。它不是把传统弹簧力换个写法，而是把约束求解放在速度重建之前。

2. PBD 投影就是局部线性化

设一个标量等式约束连接若干粒子：

$$C(\mathbf{x}_1,\ldots,\mathbf{x}_k)=0.$$

在当前位置附近做一阶展开，并把各点修正限制在约束梯度方向，可得一次质量加权投影：

$$\Delta\lambda= \frac{-C(\mathbf{x})} {\sum_j w_j\left\|\nabla_{\mathbf{x}_j}C(\mathbf{x})\right\|^2}, \qquad \Delta\mathbf{x}_j=w_j\nabla_{\mathbf{x}_j}C\,\Delta\lambda.$$

这相当于对单个非线性约束做一次 Newton 式局部修正。按顺序逐约束更新位置，就是非线性的 Projected Gauss-Seidel：后面的约束会立即看到前面已经修改过的位置。碰撞等不等式约束 $C(\mathbf{x})\ge 0$ 则只在被违反时进入活动集。

论文结论： 对平移和旋转不变的内部约束，沿质量加权梯度分配修正可保持线性、角动量，不会凭空产生拖拽物体的“幽灵力”。固定、碰撞等约束代表外界作用，本来就允许改变物体总动量。这个区分也提醒我们：以后若让网络直接修改粒子位置，就不能默认它仍具有同样的守恒性质。

传统 PBD 常再乘一个 $k\in[0,1]$ 作为“刚度”。问题是，同一约束迭代 $N$ 次后的有效修正近似为 $1-(1-k)^N$，所以增加迭代次数也会让材料变硬。原论文给出 $k'=1-(1-k)^{1/N}$ 的补偿，但也明确指出时间步依赖仍然存在；多个耦合约束下，单约束推导也不能保证全局完全不变。

3. XPBD：把刚度改写成柔度

XPBD 论文从弹性能出发。对约束向量 $\mathbf{C}(\mathbf{x})$，用柔度（compliance）$\boldsymbol{\alpha}$ 表示逆刚度：

$$U(\mathbf{x})=\frac{1}{2}\mathbf{C}(\mathbf{x})^T \boldsymbol{\alpha}^{-1}\mathbf{C}(\mathbf{x}), \qquad \widetilde{\boldsymbol{\alpha}}=\frac{\boldsymbol{\alpha}}{h^2}.$$

对单个标量约束，第 $r$ 次迭代的核心更新为：

$$\Delta\lambda= \frac{-C(\mathbf{x}^{(r)})-\widetilde{\alpha}\lambda^{(r)}} {\sum_j w_j\left\|\nabla_{\mathbf{x}_j}C\right\|^2+\widetilde{\alpha}},$$

$$\lambda^{(r+1)}=\lambda^{(r)}+\Delta\lambda, \qquad \Delta\mathbf{x}_j=w_j\nabla_{\mathbf{x}_j}C\,\Delta\lambda.$$

$\lambda$ 不再只是推导中用完即弃的比例系数，而是在一个时间步的多次约束迭代中累计的总拉格朗日乘子。每个新时间步的基线算法把它初始化为零；跨时间步 warm start 是额外策略，不应和原始算法混在一起。$\alpha=0$ 时，公式退化为无限刚约束的 PBD 投影。

论文结论： XPBD 将约束与明确的弹性势能对应，并用 $\widetilde{\alpha}=\alpha/h^2$ 处理时间尺度，使材料参数不再像 PBD 的投影比例那样直接绑定迭代次数和时间步。论文还给出由总乘子估计约束力的途径。

工程建议： “解耦”不等于“任意少迭代都得到同一结果”。柔度定义了收敛目标，但局部 Gauss-Seidel 仍需要传播修正；长链、规则网格和强耦合约束在预算不足时仍有残差。应同时记录 $C$ 的分布、迭代次数和帧时间，而不是只看一张静帧。

4. 一份可落地的 XPBD 伪代码

下面省略碰撞、摩擦和阻尼，只保留布料约束主干：

function stepXPBD(particles, constraints, h, iterations):
    # 1. 保存时间步起点，并用外力预测
    for p in particles:
        x_old[p] = x[p]
        if invMass[p] > 0:
            v[p] += h * invMass[p] * externalForce[p]
            x[p] += h * v[p]

    # 2. 总乘子在本时间步内累计
    for c in constraints:
        lambda[c] = 0

    # 3. 非线性 Gauss-Seidel 投影
    repeat iterations times:
        for c in constraints:
            C, gradients = evaluate(c, x)
            alphaTilde = compliance[c] / (h * h)
            denom = alphaTilde
            for each particle p touched by c:
                denom += invMass[p] * dot(gradients[p], gradients[p])

            if denom is too small:
                continue

            deltaLambda = (-C - alphaTilde * lambda[c]) / denom
            lambda[c] += deltaLambda

            for each particle p touched by c:
                x[p] += invMass[p] * gradients[p] * deltaLambda

    # 4. 由投影后的位置重建速度
    for p in particles:
        v[p] = (x[p] - x_old[p]) / h

工程建议： 先用串行、固定约束顺序的版本建立正确性基线，再做着色并行或 Jacobi。并行方式改变了“本轮能看到多新的位置”，收敛速度和材料观感都会改变；这不是简单的数据布局优化。

5. 拉伸约束：最小而完整的例子

对边 $(i,j)$，静止长度为 $L_{ij}$：

$$C_s=\|\mathbf{x}_i-\mathbf{x}_j\|-L_{ij}, \qquad \nabla_i C_s=\mathbf{n},\quad \nabla_j C_s=-\mathbf{n},$$

其中 $\mathbf{n}=(\mathbf{x}_i-\mathbf{x}_j)/\|\mathbf{x}_i-\mathbf{x}_j\|$。于是：

$$\Delta\lambda_s= \frac{-C_s-\widetilde{\alpha}_s\lambda_s} {w_i+w_j+\widetilde{\alpha}_s}.$$

这条式子已经覆盖固定端点：固定点的 $w=0$，所有修正自然落到可动端。实现时要为接近零长度的边设置退化保护。

工程建议： 不要把“每条边相同的 $\alpha_s$”直接解释为网格无关的连续材料。分辨率、三角形形状和约束密度改变后，宏观刚度可能变化。若要跨网格比较，应从面积、边长及目标本构关系出发做离散标定。

6. 弯曲约束：同一求解器，不同几何量

两相邻三角形共享边 $(1,2)$，对侧顶点为 $3,4$。令单位法线为 $\mathbf{n}_1,\mathbf{n}_2$，静止二面角为 $\phi_0$，原始 PBD 布料模型使用：

$$C_b=\arccos\!\left(\operatorname{clamp}(\mathbf{n}_1\cdot\mathbf{n}_2,-1,1)\right)-\phi_0.$$

记 $d=\mathbf{n}_1\cdot\mathbf{n}_2$。在三角形非退化且 $d\in(-1,1)$ 时，四个顶点的梯度都遵循同一条链式法则：

$$\nabla_{\mathbf{x}_i}C_b= -\frac{\nabla_{\mathbf{x}_i}(\mathbf{n}_1\cdot\mathbf{n}_2)}{\sqrt{1-d^2}}, \qquad i\in\{1,2,3,4\}.$$

其中 $\nabla(\mathbf{n}_1\cdot\mathbf{n}_2)$ 必须包含单位法线归一化的导数；不能把法线当常量，也不能照搬 stretch 的 $\pm\mathbf n$。平移不变性还要求四个位置梯度之和为零，这是一条很实用的实现检查。解析梯度可按 PBD 原论文给出的四点 stencil 展开，随后仍原样代入 XPBD 的 $\Delta\lambda$ 公式，只需为弯曲约束指定独立柔度 $\alpha_b$。原论文指出，这种二面角约束不依赖边长，因此比“给两个对侧顶点加距离弹簧”更容易把拉伸与弯曲分开控制。

工程建议： acos 在法线近乎平行、三角形接近退化时数值敏感。至少应 clamp 点积、拒绝零面积三角形，并对梯度做有限差分测试。若改用带符号 atan2 角度，也必须重新验证梯度、绕序和翻面行为。

7. 迭代还是 Small Steps

设渲染帧步长为 $\Delta t_f$，预算允许执行 $N$ 次约束遍历，有两种安排：

1. 一个大步长，预测一次，然后做 $N$ 次 XPBD 迭代；
2. 分成 $N$ 个 $h=\Delta t_f/N$ 的子步，每个子步只做一次 XPBD 遍历。

Small Steps in Physics Simulation报告，在其链条、布料和流体实验中，第二种安排以相近约束遍历数得到更低的约束误差与数值阻尼。直觉原因是外力预测造成的位置误差含有 $h^2$ 尺度；先把时间离散变细，后续每个局部问题也更容易。论文的算法在每个子步令乘子从零开始，因此单次遍历时，$-\widetilde\alpha\lambda$ 项在求解中可以省略。

function stepSmallSteps(frameDt, substeps):
    h = frameDt / substeps
    contacts = detectSweptContactsForWholeFrame()
    repeat substeps times:
        predictPositions(h)
        resetLambdasToZero()
        solveEachConstraintOnce(h, contacts)
        rebuildVelocities(h)

论文结论： Small Steps 比较的是“相同数量的约束遍历”，并在论文场景中观察到子步方案显著降低残差；这不是对所有碰撞系统、硬件或并行实现的普遍速度保证。

工程建议： 子步越多，预测与速度更新次数也越多；若每个子步都重建碰撞、更新动画蒙皮或调用神经网络，成本可能远高于多做一次约束迭代。论文通过整帧扫掠检测并复用接触集摊薄碰撞成本，但高速运动和轨迹变化可能让接触集过时。实际系统应分别测量 broad phase、narrow phase、约束求解和网络推理，并按风险选择接触刷新频率。小步长还会减少隐式数值耗散，必要时应显式建模约束阻尼。

8. 在接入网络前先通过四组验证

一个“看起来会飘动”的方布不足以证明求解器正确。建议先固定随机种子、约束顺序和碰撞场景，用四组小实验建立数值基线：

实验	只改变什么	应记录什么	主要排查对象
单边悬挂	质量、柔度、时间步	静态伸长量、$\lambda$	符号、单位、$h^2$ 缩放
粒子链	迭代数或子步数	最大残差、末端轨迹	修正传播与求解顺序
双三角折页	静止角、弯曲柔度	角度误差、梯度误差	法线绕序、翻面、退化面
小块悬布	网格分辨率	stretch/bend 分位数、能量趋势	离散尺度与并行差异

有限差分梯度检查尤其重要。对任意坐标分量 $x_k$，比较解析梯度和中心差分：

$$\frac{\partial C}{\partial x_k} \approx \frac{C(\mathbf{x}+\varepsilon\mathbf{e}_k)-C(\mathbf{x}-\varepsilon\mathbf{e}_k)}{2\varepsilon}.$$

误差应在远离退化构型时随 $\varepsilon$ 先下降，随后才受浮点舍入影响。若不同坐标分量呈现系统性反号，通常是梯度方向或三角形绕序错误；若只在接近平面或零面积时爆炸，则应检查角度函数和退化保护。

工程建议： 把时间步减半、迭代数翻倍、Gauss-Seidel 改为 Jacobi 都应作为独立测试轴，不能一次全改。除最终位置外，还应输出每类约束的 $|C|$ 均值、P95、最大值以及乘子范围。这样在加入 GNN 后，才能判断网络是降低了整体误差，还是只把少数严重违约藏进平均值。

9. 给 GNN 混合系统留下清晰边界

到这里可以得到一个适合作为学习系统底座的接口：

• 状态：$\mathbf{x}^n,\mathbf{v}^n$、质量、拓扑与静止几何；
• 可学习量：预测位置、速度增量、约束修正、柔度或外力先验；
• 物理校正边界：固定点，以及接触、stretch、bend 的最终校正仍交给 XPBD；有限预算下仍须检查约束残差，不能把“经过投影”等同于“严格满足”；
• 监督信号：位置误差之外，再记录约束残差、乘子/力估计和长序列漂移；
• 数值协议：训练和部署必须声明 $h$、子步数、迭代数、求解顺序及乘子重置边界。

工程建议： 第一版混合系统不要让网络同时预测“所有东西”。先固定 XPBD 为可信校正层，只让 GNN 提供一个定义清楚、可消融的输入；否则视觉改善来自网络、子步还是约束预算，很难回答。

下一篇将从粒子网格转向图表示，讨论消息传递网络如何表达布料的局部相互作用：从布料网格到图网络。

参考文献

• Matthias Müller, Bruno Heidelberger, Marcus Hennix, John Ratcliff, Position Based Dynamics, VRIPhys 2006.
• Miles Macklin, Matthias Müller, Nuttapong Chentanez, XPBD: Position-Based Simulation of Compliant Constrained Dynamics, MIG 2016.
• Miles Macklin et al., Small Steps in Physics Simulation, SCA 2019.

贡献者: Astroite