读文献primal dual PBD
Posted beidou111
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了读文献primal dual PBD相关的知识，希望对你有一定的参考价值。
文章目录
基于最优化的时间积分（变分法）
首先写出equation of motion，实际上就是动量方程或者f=ma
这里M是质量矩阵。他是对角阵，每一项分别是对应顶点的质量。
- u+代表下一时刻开始时的速度（或者是这一时刻结束时的速度），
- u-代表这一时刻开始时的速度。
- u~代表惯性作用下的下时刻速度，或者说，无约束情况下的下一时刻速度，或者说，仅仅在外力作用下的下一时刻速度。
- q代表广义坐标。如果是直角坐标系，就是位置x。之所以用广义坐标代替是因为这里还考虑了例如摆线的极坐标。和速度一样，＋代表下一时刻。
我们能够通过变分原理将动量方程转换为能量最小化问题。目标函数定义为
- 其中Ui代表势能。如果是最简单情况我们就可以认为是弹性力对应的弹性势能。
补充：为什么势能的导数是保守力？见附录
最小化问题就是求使得目标函数g最小的下一时刻速度。
求解最小化问题的时候，需要用到目标函数的梯度（例如梯度下降法）。求导得g的梯度为
其中广义力为
（把f带进去试试，会发现求导公式就是公式4）
其中G是kinematic map。它是广义速度与广义坐标的导数之间的转换矩阵。
显然，假如在笛卡尔坐标系下，G就是单位阵。
如果用最简单的显示时间积分，这一时刻的广义坐标与下一时刻的广义坐标的关系就是
梯度下降法
一个求解最小化问题（公式3）的最简单的方法就是使用梯度下降法。
首先我们有了梯度d，如公式4所示。
然后我们向着梯度d的方向下降一小步，步长是alpha
如前所述，使用最简单的显示时间积分，那么下一时刻的广义位置就是
梯度下降法的缺陷是很慢。我们可以通过预处理来进行加速。为此我们需要找到预处理矩阵P。然后用P乘以d即可。于是公式5变为。
一个对于矩阵P的简单的选择是：选取Hessian矩阵的逆的近似值。
补充：为什么预处理矩阵P应该选择Hessian的逆呢？ 见附录
加了预处理的梯度下降，此时就相当于是牛顿法了。
二次势能(primal space)
对于目标函数g（公式2）来说，可以分为惯性项和势能项。对于势能来说，使用基于约束的公式，最简单的是二次势能。
k是刚度系数。C是约束。约束可以是向量也可以是标量（后面我们主要用标量）。
根据能量最小化原理，势能的导数就是保守力。
补充：为什么会有一个G? 见附录
这里J是Constraint的Jacobian
如前所述，使用Newton style的preconditioner，那个preconditioner取Hessian矩阵。
上式中，M是质量矩阵。他是对角阵，每个元素是顶点的质量。剩下的唯一不确定的就是f对u的导数。也就是force Jacobian。 f = − ∑ i G T ∂ U i T ∂ q + \\mathbff=-\\sum_i \\mathbfG^T \\frac\\partial U_i^T\\partial \\mathbfq^+ f=−∑iGT∂q+∂UiT。对每个单元来说，我们对u求导得到
推导我们放到附录
其中第二项叫做几何刚度，geometric stiffness。 忽略几何刚度，只考虑第一项，则H的逆，也就是预处理矩阵，近似为
(GN for Gauss-Newton)
这就是Gauss Newton法迭代（一种quasi-Newton 法）的预处理矩阵。
由于求逆会耗费大量时间，所以一般用对角的倒数矩阵代替
其中下标d代表degree of freedom，也就是节点（三维每个坐标摊平了）。
附录
附录1：为什么势能的导数就是保守力。
首先我们说保守力的概念。
保守力就是做功与路径无关的力。只和起点和重点有关。
因此保守力积分就是个定积分。积分上下限分别是终点和起点。
于是积分出来的原函数，就是功。我们随便取一个零点，功就是能量。这种保守力积分出来的能量叫做势能。也就是说，势能就是保守力积分。
所以反过来，势能求导，就是保守力。
附录2：为什么预处理矩阵P应该选择Hessian的逆呢？
关于这点，请看
https://blog.csdn.net/weixin_43940314/article/details/121125847
简单来说，就是求 a r g m i n g ( u ) argmin g(u) argming(u)相当于求解其导数得0的点（驻点）。令其导数就恰好取为某个Au=b的线性方程组。也就是
g ′ ( u ) = A u − b = 0 g'(u)=Au-b = 0 g′(u)=Au−b=0
那么求解Au=b就得到了期望的u。
而求解该方程最直接的办法就是直接对A求逆。也就是
u
=
A
−
1
b
u = A^-1b
u=A−1b
就直接得到了解。
于是梯度下降法对于g来说，相当于二阶导数（第一阶导数是求驻点，第二阶导数来自于梯度）
附录3：为什么会有一个G？
这是因为f为U对x求导。所以根据链式法则，等于
f
=
∂
U
/
∂
x
=
(
∂
U
/
∂
q
)
(
∂
q
/
∂
x
)
f = \\partial U/ \\partial x = (\\partial U / \\partial q) (\\partial q / \\partial x)
f=∂U/∂x=(∂U/∂q)(∂q/∂x)
速度u与广义坐标导数的关系为
q
˙
=
G
u
\\dot q = G u
q˙=Gu
上市可以认为是广义坐标系下的速度（实际上是广义坐标的时间导数）与直角坐标系下的速度u的关系。又因为我们采用了最简单的显式时间积分，速度和位置之间只不过是乘以个dt的关系，而dt是个常数，因此是线性的。
Δ q / Δ t = G u \\Delta q/\\Delta t = Gu Δq/Δt=Gu
所以广义坐标和直角坐标系位置的转换矩阵也是矩阵G。
Δ q = G Δ x \\Delta q = G \\Delta \\mathbf x Δq=GΔx
所以 ( ∂ q / ∂ x ) (\\partial q / \\partial x) (∂q/∂x)这一项就是G。而转置不转置，只是为了内积的写法。因为内积是标量，所以最后写出来都是一样的。
由于 q ˙ = G u \\dot q = G u q˙=Gu所以显式时间积分后 Δ q = G u \\Delta q = G u Δq=Gu。因为零点对导数没影响，我们这里就认为零点是0了，所以 q = G u d t q = G u dt q=Gudt。
附录4：推导f对u求导
目标是推导出
∂ f ∂ u = − Δ t k [ J T J + ∂ J ∂ u C ] \\frac\\partial \\mathbff\\partial \\mathbfu=-\\Delta t k\\left[\\mathbfJ^T \\mathbfJ+\\frac\\partial \\mathbfJ\\partial \\mathbfu C\\right] ∂u∂f=−Δtk[JTJ+∂u∂JC]
其中
- f是保守力
- u是速度
- C是constraint function: C(q)
- J是constraint Jacobian， J = ( ∂ C / ∂ q ) G J=(\\partial C/\\partial q) G J=(∂C/∂q)G
- k是stiffness
- Δ t \\Delta t Δt是时间步长.
已知
f
=
−
G
T
∂
U
T
∂
q
\\mathbff=-\\mathbfG^T \\frac\\partial U^T\\partial \\mathbfq
f=−GT∂q∂UT
J
=
∂
C
∂
q
G
J=\\frac\\partial C\\partial q G
J=∂q∂CG
Δ
q
=
G
u
d
t
\\Delta q = G u dt
Δq=Gudt
U
=
1
2
k
C
2
U = \\frac12kC^2
U=21kC2
G是kinmatic map, 用于转换广义坐标系与直角坐标系
我们这里不加上标的q都默认是q+. 因为q-是已知量（上一时刻的广义位置），而只有下一时刻的位置q+是未知量。
由于
q
˙
=
G
u
\\dot q = G u
q˙=Gu所以显式时间积分后
Δ
q
=
G
u
\\Delta q = G u
Δq=Gu。因为零点对导数没影响，我们这里就认为零点是0了，所以
q
=
G
u
d
t
q = G u dt
q=Gudt。
∂
f
/
∂
u
=
−
∂
(
G
T
∂
U
T
∂
q
)
∂
u
=
−
G
∂
2
U
∂
q
2
∂
q
∂
u
=
−
G
2
d
t
∂
2
U
∂
q
2
=
−
G
2
d
t
k
C
∂
C
∂
q
∂
q
=
−
G
2
d
t
k
(
(
∂
C
∂
q
)
2
+
C
∂
2
C
∂
q
2
)
=
−
G
2
d
t
k
(
(
G
−
1
J
)
2
+
(
G
−
1
)
2
∂
J
∂
q
C
)
=
−
d
t
k
(
J
2
+
∂
J
∂
q
C
)
\\beginalign* \\partial f/\\partial u&=- \\frac\\partial(\\mathbfG^T \\frac\\partial U^T\\partial \\mathbfq)\\partial u \\\\ &=-G\\frac\\partial^2 U\\partial q^2 \\frac\\partial q\\partial u\\\\ &=-G^2 dt \\frac\\partial^2 U\\partial q^2\\\\ &=-G^2 dt \\frackC\\frac\\partial C\\partial q\\partial q\\\\ &=-G^2 dt \\space k \\left( (\\frac\\partial C\\partial q)^2 + C\\frac\\partial^2 C\\partial q^2\\right)\\\\ &=-G^2 dt \\space k \\left((G^-1J)^2 + (G^-1)^2 \\frac\\partial J\\partial qC\\right)\\\\ &=- dt \\space k \\left(J^2 + \\frac\\partial J\\partial qC\\right) \\endalign*
∂f/∂u=−∂u∂(GT∂q∂UT)=−G∂q2∂2U∂u∂q 以上是关于读文献primal dual PBD的主要内容，如果未能解决你的问题，请参考以下文章 primal and dual linear problem primal and dual linear problem 使用Python实现Primal-Dual Interior Point Method 使用Python实现Primal-Dual Interior Point Method