文章目录
基于最优化的时间积分（变分法）
首先写出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 + \mathbf{f}=-\sum_i \mathbf{G}^T \frac{\partial U_i^T}{\partial \mathbf{q}^{+}} 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^{-1}b
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 \mathbf{f}}{\partial \mathbf{u}}=-\Delta t k\left[\mathbf{J}^T \mathbf{J}+\frac{\partial \mathbf{J}}{\partial \mathbf{u}} 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
\mathbf{f}=-\mathbf{G}^T \frac{\partial U^T}{\partial \mathbf{q}}
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 = \frac{1}{2}kC^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 ) \begin{align*} \partial f/\partial u&=- \frac{\partial(\mathbf{G}^T \frac{\partial U^T}{\partial \mathbf{q}})}{\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 \frac{kC\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^{-1}J)^2 + (G^{-1})^2 \frac{\partial J}{\partial q}C\right)\\ &=- dt \space k \left(J^2 + \frac{\partial J}{\partial q}C\right) \end{align*} ∂f/∂u=−∂u∂(GT∂q∂UT)=−G∂q2∂2U∂u∂q=−G2dt∂q2∂2U=−G2dt∂qkC∂q∂C=−G2dt k((∂q∂C)2+C∂q2∂2C)=−G2dt k((G−1J)2+(G−1)2∂q∂JC)=−dt k(J2+∂q∂JC)
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
