文章目錄

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

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

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

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

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

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

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

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

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

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

其中廣義力為

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

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

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

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

梯度下降法

一個(gè)求解最小化問(wèn)題（公式3）的最簡(jiǎn)單的方法就是使用梯度下降法。

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

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

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

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

一個(gè)對(duì)于矩陣P的簡(jiǎn)單的選擇是：選取Hessian矩陣的逆的近似值。

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

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

二次勢(shì)能(primal space)

對(duì)于目標(biāo)函數(shù)g（公式2）來(lái)說(shuō)，可以分為慣性項(xiàng)和勢(shì)能項(xiàng)。對(duì)于勢(shì)能來(lái)說(shuō)，使用基于約束的公式，最簡(jiǎn)單的是二次勢(shì)能。

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

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

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

這里J是Constraint的Jacobian

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

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

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

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

(GN for Gauss-Newton)

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

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

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

附錄

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

首先我們說(shuō)保守力的概念。

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

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

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

所以反過(guò)來(lái)，勢(shì)能求導(dǎo)，就是保守力。

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

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

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

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

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

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

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

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

附錄3：為什么會(huì)有一個(gè)G？

這是因?yàn)閒為U對(duì)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)的時(shí)間導(dǎo)數(shù)）與直角坐標(biāo)系下的速度u的關(guān)系。又因?yàn)槲覀儾捎昧俗詈?jiǎn)單的顯式時(shí)間積分，速度和位置之間只不過(guò)是乘以個(gè)dt的關(guān)系，而dt是個(gè)常數(shù)，因此是線性的。

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

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

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

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

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

附錄4：推導(dǎo)f對(duì)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$是時(shí)間步長(zhǎng).

已知

$$\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+. 因?yàn)閝-是已知量（上一時(shí)刻的廣義位置），而只有下一時(shí)刻的位置q+是未知量。
由于$\dot{q}=Gu$所以顯式時(shí)間積分后$\Delta q=Gu$。因?yàn)榱泓c(diǎn)對(duì)導(dǎo)數(shù)沒(méi)影響，我們這里就認(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的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

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