游戏模拟——Position based dynamics

目录

• Verlet积分
  • 基本积分方法
  • Verlet 算位置
  • Verlet 算速度
• PBD
  • 基于力的方法解碰撞
    • 过冲问题
  • 基于位置的方法解碰撞
  • 算法流程
  • 求解器借用的思想
  • 关于动量守恒
  • 约束投影
  • 简单约束举例
  • 刚度体现

计算机图形中动态系统模拟最流行的方法是基于力的。累积内部和外部力量，根据牛顿的第二个运动定律计算加速度。然后使用时间积分方法来更新速度，最后是对象的位置。
一些模拟方法（大多数刚性体模拟器）使用基于冲量的方法并直接操纵速度。
PBD是一种省略了速度层的直接作用于位置的控制方法，方法的主要优点是它的可控性、性能、数值稳定，在游戏开发领域应用很广，近年的UE chaos底层也采用了PBD方法。

Verlet积分#

verlet积分在游戏开发领域应用也很广泛，并且PBD论文中提到了PBD的思想和verlet有相近之处，这里也顺便介绍一下。

基本积分方法#

积分方法的阶数一般是看泰勒展开中导数的最高项。n阶方法的误差是O(hn+1)" role="presentation">O(hn+1)$O(h^{n+1})$。
值得一提的是，explicit euler、semi-implicit euler都是一阶方法。
积分方法讨论的通常是给定一个函数S(t)，然后我们没有办法得到S(t)的解析式（假定这里是和时间有关的函数），此时我们通过一些方法来对S(t) 做近似，达到一个我们给定一个时刻的数据，可以推算其他时刻数据的效果，这里一切的理论基础都是泰勒展开，然后进行各种化简。
这里用图来表示，小红色的矩形代表积分过程，左边是explicit euler，右边是semi-implicit euler。

左图用的矩形的左边来表示一个小时间步里整体的情况，意思就是时间步刚开始的值。左图用的矩形的右边来表示一个小时间步里整体的情况，意思就是时间步结束时的值。

Verlet 算位置#

verlet积分是一种牛顿运动方程的积分方法，他的积分方法有着良好的数值稳定性，对比简单的欧拉方法没有显著的额外计算成本。

xt+Δt=xt+v(t)Δt+12a(t)Δt2+16b(t)Δt3+O(Δt4)" role="presentation">xt+Δt=xt+v(t)Δt+12a(t)Δt2+16b(t)Δt3+O(Δt4)$$x_{t+\mathrm{\Delta }t}=x_t+v(t)\mathrm{\Delta }t+\frac{1}{2}a(t)\mathrm{\Delta }{t}^2+\frac{1}{6}b(t)\mathrm{\Delta }{t}^3+O(\mathrm{\Delta }{t}^4)$$

xt−Δt=xt−v(t)Δt+12a(t)Δt2−16b(t)Δt3+O(Δt4)" role="presentation">xt−Δt=xt−v(t)Δt+12a(t)Δt2−16b(t)Δt3+O(Δt4)$$x_{t-\mathrm{\Delta }t}=x_t-v(t)\mathrm{\Delta }t+\frac{1}{2}a(t)\mathrm{\Delta }{t}^2-\frac{1}{6}b(t)\mathrm{\Delta }{t}^3+O(\mathrm{\Delta }{t}^4)$$

两者求和，消去一阶和三阶项得到

xt+Δt=2xt−xt−Δt+a(t)Δt2+O(Δt4)" role="presentation">xt+Δt=2xt−xt−Δt+a(t)Δt2+O(Δt4)$$x_{t+\mathrm{\Delta }t}=2x_t-x_{t-\mathrm{\Delta }t}+a(t)\mathrm{\Delta }{t}^2+O(\mathrm{\Delta }{t}^4)$$

我们可以将t+Δt" role="presentation">t+Δt$t+\mathrm{\Delta }t$ 当成当前帧所在的时刻，那么t就是上一帧，t−Δt" role="presentation">t−Δt$t-\mathrm{\Delta }t$就是上上一帧，我们只需要两个位置信息和上一帧的力信息，就可以得到当前帧的位置信息，并且这个位置信息是三阶精度的，误差四阶精度，wiki百科中提到，这个方法这里的误差是局部误差，并且公式中的加速度是根据精确解计算的，全局误差（即累积误差）与时间步长的平方（Δt^2）成正比。虽然局部离散化误差为 4 阶，但由于微分方程的2阶，推导一番后全局误差为 2 阶，我们将它定义为2阶方法。

有的书上也说这个公式也可以从二阶中心差分推导出，所以是二阶方法。

Verlet 算速度#

接下来是求速度信息，
两者求差，消去二阶项得到

x(t+Δt)−x(t−Δt)=2v(t)Δt+13b(t)Δt3" role="presentation">x(t+Δt)−x(t−Δt)=2v(t)Δt+13b(t)Δt3$$x(t+\mathrm{\Delta }t)-x(t-\mathrm{\Delta }t)=2v(t)\mathrm{\Delta }t+\frac{1}{3}b(t)\mathrm{\Delta }{t}^3$$

这样可以推算上一帧的速度信息

v(t)=x(t+Δt)−x(t−Δt)2Δt+b(t)Δt33Δt" role="presentation">v(t)=x(t+Δt)−x(t−Δt)2Δt+b(t)Δt33Δt$$v(t)=\frac{x(t+\mathrm{\Delta }t)-x(t-\mathrm{\Delta }t)}{2\mathrm{\Delta }t}+\frac{b(t)\mathrm{\Delta }{t}^3}{3\mathrm{\Delta }t}$$

v(t)=x(t+Δt)−x(t−Δt)2Δt+O(Δt2)" role="presentation">v(t)=x(t+Δt)−x(t−Δt)2Δt+O(Δt2)$$v(t)=\frac{x(t+\mathrm{\Delta }t)-x(t-\mathrm{\Delta }t)}{2\mathrm{\Delta }t}+O(\mathrm{\Delta }{t}^2)$$

这个算上一帧速度的方法是一阶精度二阶误差的。
在他的原始方法里推算当帧的速度用的中值定理+近似，导致速度一定是不准的，得出速度公式为

v(t+Δt)=x(t+Δt)−x(t)Δt+O(Δt)" role="presentation">v(t+Δt)=x(t+Δt)−x(t)Δt+O(Δt)$$v(t+\mathrm{\Delta }t)=\frac{x(t+\mathrm{\Delta }t)-x(t)}{\mathrm{\Delta }t}+O(\mathrm{\Delta }t)$$

这个速度的准度可以说几乎没有，零阶精度，一阶误差，推导可以看verlet一开始 x(t+Δt)" role="presentation">x(t+Δt)$x(t+\mathrm{\Delta }t)$的展开将xt移到左边再除以步长得到：

x(t+Δt)−x(t)Δt=v(t)+12a(t)Δt+16b(t)Δt2" role="presentation">x(t+Δt)−x(t)Δt=v(t)+12a(t)Δt+16b(t)Δt2$$\frac{x(t+\mathrm{\Delta }t)-x(t)}{\mathrm{\Delta }t}=v(t)+\frac{1}{2}a(t)\mathrm{\Delta }t+\frac{1}{6}b(t)\mathrm{\Delta }{t}^2$$

所以如果12a(t)Δt+16b(t)Δt2" role="presentation">12a(t)Δt+16b(t)Δt2$\frac{1}{2}a(t)\mathrm{\Delta }t+\frac{1}{6}b(t)\mathrm{\Delta }{t}^2$作为误差，误差就是O(Δt)" role="presentation">O(Δt)$O(\mathrm{\Delta }t)$的。

wiki百科里则提供了一种更加高精度的速度推算方法叫 velocity verlet ，利用梯形法（头尾值加起来求平均，之前的小矩形可以想象成梯形来看）算的加速度来推算速度，让速度的精度到达二阶，至于为什么是二阶，因为梯形法的结果是二阶精度，这个精度貌似可以传递。

v(t+Δt)=v(t)+a(t)+a(t+Δt)2Δt" role="presentation">v(t+Δt)=v(t)+a(t)+a(t+Δt)2Δt$$v(t+\mathrm{\Delta }t)=v(t)+\frac{a(t)+a(t+\mathrm{\Delta }t)}{2}\mathrm{\Delta }t$$

a的t+Δt" role="presentation">t+Δt$t+\mathrm{\Delta }t$是一般来说是一个和位置有关的函数，先求了位置，接下来这里面都是已知量了。

PBD#

基于力的方法解碰撞#

首先一个直观印象关于如何处理动力学中的穿插问题，基于力的方法往往是计算穿插深度，然后给根据深度和穿插的方向给予力，这过程中有一个刚度系数，表示这个物体的坚硬程度，调的越大，两个物体越难穿插，你就会看到物体看起来很硬，实际航就是穿插产生的力的缩放系数，刚度越大实际上要求时间步长越小，否则很容易出现穿插过深弹出一个很大的力这种情况，其实就是数值解误差太大。

过冲问题#

显式积分法，也称为开环法或单步法，广泛用于求解科学和工程各个领域的常微分方程 (ODE)。然而，这些方法可能会遇到过冲问题，这会导致不准确或不稳定的解决方案。
产生原因：
1、两刚体相交过深导致计算出来的力太大，求解的初始条件残差太大
2、时间步长太大
3、高频振荡：在求解涉及高频振荡的 ODE 时，显式方法可能难以准确捕获这些振荡，从而导致超调。在这种情况下，使用对刚性系统更稳定的隐式积分方法会有所帮助。
4、误差累积：使用显式方法时，误差会随着时间的推移而累积，从而导致显着的超调。

要解决显式积分中的过冲问题，请考虑使用更小的时间步长、自适应时间步长、高阶方法，或者在适当的时候切换到隐式或半隐式方法。此外，优化初始条件并使用灵敏度分析可以提高解决方案的准确性和稳定性。

基于位置的方法解碰撞#

在PBD方法中，当检测到两个物体发生穿透时，直接根据约束修正物体位置到一个不会碰撞的位置，然后更新速度信息。他其实是一个反向的过程，虽然这中间力的计算不明确了，但是表现是正常的，这是我们期待的，这就是pbd论文提到的visually plausible 视觉正确，并且避免了之前提到的过冲问题。

算法流程#

我们表达一个物体，使用N个顶点，M个约束，一个顶点i对应的质量是mi，位置则是xi，速度 vi
约束的索引一般使用j
我们对于其中一个约束有：

• 定义一个约束影响的顶点数量是nj，可以理解为第 j个约束所影响的顶点数目为nj
• 定义约束函数C 顶点数目为3nj，位置信息为xyz，所有的数据总量就是3nj，约束函数就是读取这个3nj的数据算出一个实数
• 约束函数影响的顶点的索引，nj个索引
• 定义这个约束的刚度值
• 定义这个约束是等式约束和不等式约束
  等式约束比较强硬，要满足约束函数等于0，不等式约束则是约束函数大于等于0就行，kj定义了这个约束的刚度，刚度越大这个约束要越快被满足。
  

算法最开始是先初始化x和v，还有质量的倒数
每个循环开始，先通过力算速度，在通过算出的速度算位置，计算出一个预测的速度、位置。
然后在根据这个预测的位置情况进行碰撞约束的生成。
然后走一个固定的迭代次数，做约束投影，让p，其实就是位置，落在正确的 满足约束函数的地方，因为有多个约束，其实这里就是一个多个约束方程的方程组。
约束投影完了只有得到正确的位置，然后根据投影后的位置和上一次迭代的位置算一下速度（这里和verlet简单的算速度方法很像），然后把投影后的位置更新到位置数据里。
最后在(16)行根据摩擦(friction)和恢复(restitution)系数修改速度。

这个方案是无条件稳定的，因为他不会像一般的explicit方法那样对未来进行外插，而是将顶点移动到正确的位置。
他的稳定性不取决于时间步长，而是取决于约束函数的样子。
这个求解器并不像以往的显式或者隐式积分器，而是处于一种中间态，如果每个时间步只跑一个迭代他看起来像explict，而添加迭代次数，可以更像隐式积分方案。

求解器借用的思想#

需要注意的是约束函数是通常都是非线性的，比如距离约束

C=|p1−p2|−d" role="presentation">C=|p1−p2|−d$$C=|p_1-p_2|-d$$

这里的非线性，其实指的是他这个方程组的未知数不是一次的，有绝对值的话，可以理解成平方再开方，一般的迭代法解线性方程组使用jacobi、gauss-seidel迭代方法，他们的未知数有多个，但是都是一次的。
但是在这里作者是借用了GS方法的思想，来解一个非线性方程组。
原gauss-seidel迭代法是，我们一整个方程组，然后从第一条方程开始，代入一个初值，然后方程左边算出第一个未知数X，然后这个X算出来后会直接替换掉下面方程组求解时用到X的地方，第二条方程算出Y，也是一样的，先满足一条方程，然后将它的影响带到下一条方程，经过固定次数的迭代后，系统的未知数们会收敛，就求解成功了。
作者称这个方法为Non-Linear Gauss-Seidel方法，本质上就是把所有的约束方程列出来，然后带入一个初值，然后先求一条约束方程，得到一个delta p，然后把这个delta p应用，然后带入第二条约束方程继续求，直到最后一条方程，如此迭代多次，即使是非线性方程，最后方程组也会收敛。
值得注意的是，GS方案比较依赖方程的组的满足顺序，过度约束 且 没有保证求解顺序会导致结果振荡。

关于动量守恒#

PBD方法里提到了动量守恒的话题，PBD本身方法的流程他是想说明按照他的流程走，只要约束函数定义时候是动量守恒的，那么结果是动量守恒的。
如果你这个系统最后输出的结果是没有违背动量守恒的基础，那么将看到一些不期待的运动，这里坐着定义为 Ghost force，看着就像有力量拖着物体 旋转物体一样。

后面作者提到，对于一个系统来说只有内部约束需要考虑动量守恒，因为外力作用会作用在全局而产生动量。内部约束不会改变刚体运动状态，即刚体的旋转、位移，作者称之为 与刚体模态 正交。
因此我们沿着约束函数梯度的方向去影响整个系统他最后满足约束了约束函数（满足这个内部力的性质），他的最后就是动量守恒的。
这里我需要说的详细的一点，因为我自己也对动量守恒这块比较困惑：
我的理解是PBD本身其实和动量守恒不沾边，主要是约束函数的定义是动量守恒的，你沿着约束函数的梯度去影响系统最后达到了约束函数，那么约束函数的最终状态是动量守恒的，我们拿这个状态去做表现，表现就是动量守恒的，但是如果我们因为很多约束，最后求解没成功，有很多约束没满足，那么最后的表现就是不动量守恒的。
如果我们定义了一个动量不守恒的C，最后求解成功了，那么表现出来也就是动量不守恒的，如果胡乱定义了一个约束还要求动量守恒 其实没有讨论的必要。

约束投影#

这块公式相当多，
首先对于一个等式约束我们要求的就是这个状态代入约束函数后输出一个0
即

C(p+Δp)=0" role="presentation">C(p+Δp)=0$$C(p+\mathrm{\Delta }p)=0$$

而这个函数可以一阶泰勒展开得到

C(p+Δp)≈C(p)+∇pC(p)⋅Δp" role="presentation">C(p+Δp)≈C(p)+∇pC(p)⋅Δp$$C(p+\mathrm{\Delta }p)\approx C(p)+{\mathrm{\nabla }}_{p}C(p)\cdot \mathrm{\Delta }p$$

我们要把delta P限制在约束函数的梯度方向即

Δp=λ∇pC(p)" role="presentation">Δp=λ∇pC(p)$$\mathrm{\Delta }p=\lambda {\mathrm{\nabla }}_{p}C(p)$$

这个式子代入上一个式子得到λ
得到

Δp=−C(p)|∇pC(p)|2∇pC(p)" role="presentation">Δp=−C(p)|∇pC(p)|2∇pC(p)$$\mathrm{\Delta }p=-\frac{C(p)}{|{\mathrm{\nabla }}_{p}C(p){|}^2}{\mathrm{\nabla }}_{p}C(p)$$

这就是对于单个约束函数，整个系统接下来要做的改变。
而对于单个点p_i来说，前面的系数λ是一样的，后面的梯度方向则变成了这个delta P在单个点的分量，其实是梯度的某几个维度，也就是关于pi的偏导数
前面的系数用s表示，

Δpi=−s∇piC(p)" role="presentation">Δpi=−s∇piC(p)$$\mathrm{\Delta }{p}_i=-s{\mathrm{\nabla }}_{p_i}C(p)$$

s=C(p)∑j|∇pjC(p)|2" role="presentation">s=C(p)∑j|∇pjC(p)|2$$s=\frac{C(p)}{\sum _j|{\mathrm{\nabla }}_{p_j}C(p){|}^2}$$

其实这里下面的分母就是梯度的模长，梯度的模长等于各个偏导数分量的平方和。
特别注意下维度：
这里的p是指所有的点的位置，维度是点的数量 n * 3 ，其中delta P会保证和约束函数的梯度一样的方向，约束函数的梯度也是一个 n * 3的向量 。

delta P1则表示delta P向量的第一点分量，是一个3分量的向量，其实就是第一个点的位移，他应该保证和约束函数在这个第一个点的偏导数保证方向相同。

简单约束举例#

这里就不举例了，可以查看参考文章里面的距离约束
https://zhuanlan.zhihu.com/p/48737753

刚度体现#

公式里还有一个k参数，这个参数的作用意会一下就是，我们有一条约束方程，要以怎么样的速度去满足这条方程，如果k等于1则说明直接满足，不留余地，等于0-1的值，则说明可能是跑几轮才能比较接近精确解。
刚度越大，越快满足约束方程，也就和刚度体现的越硬对应了。
一个问题是，他自己指出材料的刚度依赖时间步长，定步长就不太会有这些问题，verlet积分的推断也是依赖定步长。
多次迭代下误差为

Δp(1−k)ns" role="presentation">Δp(1−k)ns$$\mathrm{\Delta }p(1-k{)}^{n_s}$$

另外，为了消除多次迭代中残差随着次数指数级增长的问题，他这里用的一个非常怪的式子反推了，然后误差变为线性依赖于k而不依赖于ns,ns是迭代次数。

k′=1−(1−k)1ns" role="presentation">k′=1−(1−k)1ns$$k^{\prime }=1-(1-k{)}^{\frac{1}{n_s}}$$

最终误差来到

Δp(1−k′)ns=Δp(1−k)" role="presentation">Δp(1−k′)ns=Δp(1−k)$$\mathrm{\Delta }p(1-k^{\prime }{)}^{n_s}=\mathrm{\Delta }p(1-k)$$

flyingziming
2023.4.7