论文阅读笔记：Complementary Dynamics

论文题目：Complementary Dynamics 。It is a SIGGRAPH Asia 2020 paper by Jiayi Eris Zhang, Seungbae Bang, David I.W. Levin and Alec Jacobson. More details is in the Project Page

---

0. 摘要

论文提出了一种

1. 背景

为了实现****

2. 方法初探

首先，定义仿真所涉及的变量。定义骨架参数 rig parameter 为 \(\boldsymbol{p} \in \mathbb{R}^{m}\) ；基于骨架蒙皮结构，由骨架偏移实现的网格节点位移为 \(\boldsymbol{u}^{r}\) ；基于本文所提出的方法，计算得到的 complementary displacement 网格节点位移为 \(\boldsymbol{u}^{c}\) ；在 \(t\) 时刻，网格节点的位移是两者相加的和，即

\[\boldsymbol{u}_{t} = \boldsymbol{u}^{r}_{t} + \boldsymbol{u}^{c}_{t} \]

在软体对象的形变仿真中，可基于优化隐式积分方法进行求解计算。 \(t\) 时刻网格节点的位移为如下非线性最优化问题 non-linear optimization problem 的解，即

\[\boldsymbol{u}_{t} = argmin E_{t}(\boldsymbol{u}_{t}) \]

其中，\(E_{t}\) 是能量函数，包括了弹性势能 potential energy 、动能 momentum term ，外力做功 external work 等，即

\[E_{t}(\boldsymbol{u}_{t}) = \Phi(\boldsymbol{u}_{t}) + \frac{1}{2h^2} \ddot{\boldsymbol{u}_{t}}^T\boldsymbol{M}\ddot{\boldsymbol{u}_{t}} - \boldsymbol{u}_{t}^T\boldsymbol{f}(\boldsymbol{u}_{t}) \]

其中，\(h > 0\) 是仿真时间步长；\(\boldsymbol{M} \in \mathbb{R}^{dn\times dn}\) 是质量矩阵 mass matrix 。

那么，在本论文的仿真计算中，complementary displacement 网格节点位移 \(\boldsymbol{u}^{c}_{t}\) 可通过求解上述优化问题得到，（注：\(\boldsymbol{u}^{r}_{t}\) 是已知量），即

\[\boldsymbol{u}^{c}_{t} = argmin E_{t}(\boldsymbol{u}^{r}_{t} + \boldsymbol{u}^{c}_{t}) \]

思考：对上述问题直接求解，得到的结果将会是 \(\boldsymbol{u}^{c}_{t} = -\boldsymbol{u}^{r}_{t}\) ，也就是会 undo the rig displacements ，显然是不对的。因此，在仿真计算中，应加入适当约束，对 complementary displacement \(\boldsymbol{u}^{c}_{t}\) 施加一定的限制。

3. 骨架正交约束 Rig Orthogonality Constraints

接下来思考如何定义约束。换一个视角，如果 \(\boldsymbol{u}^{c}_{t}\) 和 \(\boldsymbol{u}^{r}_{t}\) 互不相关，那么，给定当前的骨架参数 \(\boldsymbol{p}_{t}\) ，对于任意有效的 \(\boldsymbol{u}^{c}_{t}\) 都应满足如下标准，即

\[argmin \frac{1}{2}||\boldsymbol{u}^{r}(\boldsymbol{p}_t) + \boldsymbol{u}^{c}_{t} - \boldsymbol{u}^{r}(\boldsymbol{p})||^2_{\boldsymbol{M}} = \boldsymbol{p}_t \]

其中，\(||\boldsymbol{x}||^2_{\boldsymbol{M}} = \boldsymbol{x}^T\boldsymbol{M}\boldsymbol{x}\) 。该约束将保证 \(\boldsymbol{u}^{c}_{t}\) 和 \(\boldsymbol{u}^{r}_{t}\) 互不相关。
（上面这一段的表述，还需要再整理，感觉没讲清楚）

有了上述条件，我们取优化目标函数的一阶导数，即

\[\frac{\partial\frac{1}{2}||\boldsymbol{u}^{r}(\boldsymbol{p}_t) + \boldsymbol{u}^{c}_{t} - \boldsymbol{u}^{r}(\boldsymbol{p})||^2_{\boldsymbol{M}}}{\partial\boldsymbol{p}}|_{\boldsymbol{p}_t} = 0 \]

将上述公式拆分，得到

\[\frac{\partial(\boldsymbol{u}^{r}(\boldsymbol{p}_t) + \boldsymbol{u}^{c}_{t} - \boldsymbol{u}^{r}(\boldsymbol{p}))}{\partial\boldsymbol{p}}|^T_{\boldsymbol{p}_t} \boldsymbol{M}(\boldsymbol{u}^{r}(\boldsymbol{p}_t) + \boldsymbol{u}^{c}_{t} - \boldsymbol{u}^{r}(\boldsymbol{p})) = 0 \]

由于 \(\boldsymbol{u}^{r}(\boldsymbol{p}_t)\) 以及 \(\boldsymbol{u}^{c}_{t}\) 均与 \(\boldsymbol{p}\) 无关，对其偏导数为零，上述公式进一步简化为

\[\frac{\partial\boldsymbol{u}^{r}}{\partial\boldsymbol{p}}|^{T}_{\boldsymbol{p}_t} \boldsymbol{M} \boldsymbol{u}^{c}_{t} = 0 \]

其中，\(\boldsymbol{J}^T_t = \frac{\partial\boldsymbol{u}^{r}}{\partial\boldsymbol{p}}|^{T}_{\boldsymbol{p}_t}\) 为骨架的雅可比矩阵 rig Jacobian matrix 。

上述约束意味着，\(\boldsymbol{u}^{c}\) 与骨架蒙皮的位移正交互补: This constraint can read as requiring that \(\boldsymbol{u}^{c}\) lies in the orthogonal complement of the rig, first order.

4. 约束仿真 Constrainted Simulation

基于上述骨架正交约束，可以重新给出形变仿真的优化问题，即

\[\boldsymbol{u}^{c}_{t} = argmin E_{t}(\boldsymbol{u}^{r}_{t} + \boldsymbol{u}^{c}) \quad \text{subject to} \quad \boldsymbol{J}^T_t \boldsymbol{M} \boldsymbol{u}^{c} = 0 \]

对于上述优化问题的求解，通常可采用拉格朗日乘子法，即上述优化问题可转化为如下线性系统：

\[\begin{bmatrix} \boldsymbol{Q} & \boldsymbol{C}^{T}\\ \boldsymbol{C} & \boldsymbol{0} \end{bmatrix} \begin{bmatrix} \boldsymbol{x}\\ \lambda \end{bmatrix} = \begin{bmatrix} \boldsymbol{l}\\ \boldsymbol{0} \end{bmatrix} \]

其中，\(\boldsymbol{Q} = \partial^2\Phi/\partial{\boldsymbol{u}^{C}}^2 + h^{-2}\boldsymbol{M}\)，是与 Hessian 矩阵相关的项；\(\boldsymbol{l} = -\partial\Phi/\partial\boldsymbol{u}^{C} + \frac{\boldsymbol{M}}{h^2}(\boldsymbol{u}_{t}^{r} - \boldsymbol{u}_{t-h} - h\dot{\boldsymbol{u}}_{t-h}) + \boldsymbol{f}_{t}\)，是 Gradient 相关的项；\(\boldsymbol{C} = \boldsymbol{J}^T_t\boldsymbol{M}\)，是骨架 Jacobian 相关的项。

5. 优化问题求解 Newton-Raphson方法

采用拉格朗日乘子法对上述优化问题求解，最终转化为非线性方程组 \(\boldsymbol{Ax} = \boldsymbol{b}\) 。在本文中，采用了 Newton-Raphson 方法对该非线性方程组进行求解。结合 Newton-Raphson方法求解多元非线性方程组 中给出的介绍，求解流程如下：

1、初始化。计算 rig 分量 \(\boldsymbol{u}_t^{r}\)，计算初始位移 \(\boldsymbol{u}^{(0)}\)，complementary 分量为零，即 \(\boldsymbol{u}^{c(0)} = \boldsymbol{0}\) ；最大迭代次数为 \(N = 20\) ；置 \(k = 0\)；精度要求为 \(\epsilon = 1e^{-6}\)

2、求解方程组 \(\boldsymbol{A}(\boldsymbol{x}^{(k)})\boldsymbol{d} = -\boldsymbol{F}(\boldsymbol{x}^{(k)})\)

注：Newton-Raphson方法求解多元非线性方程组参见 Newton-Raphson方法求解多元非线性方程组