文章目錄

• • 基于最優(yōu)化的時間積分（變分法）
  • 梯度下降法
  • 二次勢能(primal space)
  • 附錄
  • • 附錄1：為什么勢能的導(dǎo)數(shù)就是保守力。
    • 附錄2：為什么預(yù)處理矩陣P應(yīng)該選擇Hessian的逆呢？
    • 附錄3：為什么會有一個G？
    • 附錄4：推導(dǎo)f對u求導(dǎo)

基于最優(yōu)化的時間積分（變分法）

首先寫出equation of motion，實(shí)際上就是動量方程或者f=ma

這里M是質(zhì)量矩陣。他是對角陣，每一項分別是對應(yīng)頂點(diǎn)的質(zhì)量。

• u+代表下一時刻開始時的速度（或者是這一時刻結(jié)束時的速度），
• u-代表這一時刻開始時的速度。
• u~代表慣性作用下的下時刻速度，或者說，無約束情況下的下一時刻速度，或者說，僅僅在外力作用下的下一時刻速度。
  
• q代表廣義坐標(biāo)。如果是直角坐標(biāo)系，就是位置x。之所以用廣義坐標(biāo)代替是因為這里還考慮了例如擺線的極坐標(biāo)。和速度一樣，＋代表下一時刻。

我們能夠通過變分原理將動量方程轉(zhuǎn)換為能量最小化問題。目標(biāo)函數(shù)定義為

• 其中Ui代表勢能。如果是最簡單情況我們就可以認(rèn)為是彈性力對應(yīng)的彈性勢能。

補(bǔ)充：為什么勢能的導(dǎo)數(shù)是保守力？見附錄

最小化問題就是求使得目標(biāo)函數(shù)g最小的下一時刻速度。

求解最小化問題的時候，需要用到目標(biāo)函數(shù)的梯度（例如梯度下降法）。求導(dǎo)得g的梯度為

其中廣義力為

（把f帶進(jìn)去試試，會發(fā)現(xiàn)求導(dǎo)公式就是公式4）

其中G是kinematic map。它是廣義速度與廣義坐標(biāo)的導(dǎo)數(shù)之間的轉(zhuǎn)換矩陣。

顯然，假如在笛卡爾坐標(biāo)系下，G就是單位陣。

如果用最簡單的顯示時間積分，這一時刻的廣義坐標(biāo)與下一時刻的廣義坐標(biāo)的關(guān)系就是

梯度下降法

一個求解最小化問題（公式3）的最簡單的方法就是使用梯度下降法。

首先我們有了梯度d，如公式4所示。

然后我們向著梯度d的方向下降一小步，步長是alpha

如前所述，使用最簡單的顯示時間積分，那么下一時刻的廣義位置就是

梯度下降法的缺陷是很慢。我們可以通過預(yù)處理來進(jìn)行加速。為此我們需要找到預(yù)處理矩陣P。然后用P乘以d即可。于是公式5變?yōu)椤?/p>

一個對于矩陣P的簡單的選擇是：選取Hessian矩陣的逆的近似值。

補(bǔ)充：為什么預(yù)處理矩陣P應(yīng)該選擇Hessian的逆呢？ 見附錄

加了預(yù)處理的梯度下降，此時就相當(dāng)于是牛頓法了。

二次勢能(primal space)

對于目標(biāo)函數(shù)g（公式2）來說，可以分為慣性項和勢能項。對于勢能來說，使用基于約束的公式，最簡單的是二次勢能。

k是剛度系數(shù)。C是約束。約束可以是向量也可以是標(biāo)量（后面我們主要用標(biāo)量）。

根據(jù)能量最小化原理，勢能的導(dǎo)數(shù)就是保守力。

補(bǔ)充：為什么會有一個G? 見附錄

這里J是Constraint的Jacobian

如前所述，使用Newton style的preconditioner，那個preconditioner取Hessian矩陣。

上式中，M是質(zhì)量矩陣。他是對角陣，每個元素是頂點(diǎn)的質(zhì)量。剩下的唯一不確定的就是f對u的導(dǎo)數(shù)。也就是force Jacobian。$\mathbf{f}=?{\sum }_i{\mathbf{G}}^T\frac{?U_i^T}{?{\mathbf{q}}^{+}}$。對每個單元來說，我們對u求導(dǎo)得到

推導(dǎo)我們放到附錄

其中第二項叫做幾何剛度，geometric stiffness。 忽略幾何剛度，只考慮第一項，則H的逆，也就是預(yù)處理矩陣，近似為

(GN for Gauss-Newton)

這就是Gauss Newton法迭代（一種quasi-Newton 法）的預(yù)處理矩陣。

由于求逆會耗費(fèi)大量時間，所以一般用對角的倒數(shù)矩陣代替

其中下標(biāo)d代表degree of freedom，也就是節(jié)點(diǎn)（三維每個坐標(biāo)攤平了）。

附錄

附錄1：為什么勢能的導(dǎo)數(shù)就是保守力。

首先我們說保守力的概念。

保守力就是做功與路徑無關(guān)的力。只和起點(diǎn)和重點(diǎn)有關(guān)。

因此保守力積分就是個定積分。積分上下限分別是終點(diǎn)和起點(diǎn)。

于是積分出來的原函數(shù)，就是功。我們隨便取一個零點(diǎn)，功就是能量。這種保守力積分出來的能量叫做勢能。也就是說，勢能就是保守力積分。

所以反過來，勢能求導(dǎo)，就是保守力。

附錄2：為什么預(yù)處理矩陣P應(yīng)該選擇Hessian的逆呢？

關(guān)于這點(diǎn)，請看

https://blog.csdn.net/weixin_43940314/article/details/121125847

簡單來說，就是求$argming(u)$相當(dāng)于求解其導(dǎo)數(shù)得0的點(diǎn)（駐點(diǎn)）。令其導(dǎo)數(shù)就恰好取為某個Au=b的線性方程組。也就是

$$g^{\prime }(u)=Au?b=0$$

那么求解Au=b就得到了期望的u。

而求解該方程最直接的辦法就是直接對A求逆。也就是
$$u=A^{?1}b$$
就直接得到了解。

于是梯度下降法對于g來說，相當(dāng)于二階導(dǎo)數(shù)（第一階導(dǎo)數(shù)是求駐點(diǎn)，第二階導(dǎo)數(shù)來自于梯度）

附錄3：為什么會有一個G？

這是因為f為U對x求導(dǎo)。所以根據(jù)鏈?zhǔn)椒▌t，等于
$$f=?U/?x=(?U/?q)(?q/?x)$$

速度u與廣義坐標(biāo)導(dǎo)數(shù)的關(guān)系為
$$\dot{q}=Gu$$
上市可以認(rèn)為是廣義坐標(biāo)系下的速度（實(shí)際上是廣義坐標(biāo)的時間導(dǎo)數(shù)）與直角坐標(biāo)系下的速度u的關(guān)系。又因為我們采用了最簡單的顯式時間積分，速度和位置之間只不過是乘以個dt的關(guān)系，而dt是個常數(shù)，因此是線性的。

$$\Delta q/\Delta t=Gu$$

所以廣義坐標(biāo)和直角坐標(biāo)系位置的轉(zhuǎn)換矩陣也是矩陣G。

$$\Delta q=G\Delta \mathbf{x}$$

所以$(?q/?x)$這一項就是G。而轉(zhuǎn)置不轉(zhuǎn)置，只是為了內(nèi)積的寫法。因為內(nèi)積是標(biāo)量，所以最后寫出來都是一樣的。

由于$\dot{q}=Gu$所以顯式時間積分后$\Delta q=Gu$。因為零點(diǎn)對導(dǎo)數(shù)沒影響，我們這里就認(rèn)為零點(diǎn)是0了，所以$q=Gudt$。

附錄4：推導(dǎo)f對u求導(dǎo)

目標(biāo)是推導(dǎo)出

$$\frac{?\mathbf{f}}{?\mathbf{u}}=?\Delta tk\left[{\mathbf{J}}^T\mathbf{J}+\frac{?\mathbf{J}}{?\mathbf{u}}C\right]$$

其中

• f是保守力
• u是速度
• C是constraint function: C(q)
• J是constraint Jacobian， $J=(?C/?q)G$
• k是stiffness
• $\Delta t$是時間步長.

已知

$$\mathbf{f}=?{\mathbf{G}}^T\frac{?U^T}{?\mathbf{q}}$$
$$J=\frac{?C}{?q}G$$
$$\Delta q=Gudt$$
$$U=\frac{1}{2}k{C}^2$$

G是kinmatic map, 用于轉(zhuǎn)換廣義坐標(biāo)系與直角坐標(biāo)系

我們這里不加上標(biāo)的q都默認(rèn)是q+. 因為q-是已知量（上一時刻的廣義位置），而只有下一時刻的位置q+是未知量。
由于$\dot{q}=Gu$所以顯式時間積分后$\Delta q=Gu$。因為零點(diǎn)對導(dǎo)數(shù)沒影響，我們這里就認(rèn)為零點(diǎn)是0了，所以$q=Gudt$。

$$\begin{array}{rl}?f/?u& =?\frac{?({\mathbf{G}}^T\frac{?U^T}{?\mathbf{q}})}{?u}\\ & =?G\frac{{?}^{2}U}{?q^2}\frac{?q}{?u}\\ & =?G^{2}dt\frac{{?}^{2}U}{?q^2}\\ & =?G^{2}dt\frac{kC\frac{?C}{?q}}{?q}\\ & =?G^{2}dtk\left((\frac{?C}{?q}{)}^2+C\frac{{?}^{2}C}{?q^2}\right)\\ & =?G^{2}dtk\left((G^{?1}J{)}^2+(G^{?1}{)}^2\frac{?J}{?q}C\right)\\ & =?dtk\left(J^2+\frac{?J}{?q}C\right)\end{array}$$

總結(jié)

以上是生活随笔為你收集整理的【读文献】primal dual PBD的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。