这个系列是GAMES103：基于物理的计算机动画入门(GAMES 103: Intro to Physics Based Animation)的同步课程笔记。本课程是计算机动画的入门课程，着重介绍各种基于物理的动画仿真模拟技术。本节课主要介绍布料模拟中的约束方法。
PBD and Strain Limiting
布料的力学特征之一是当它拉伸到一定程度后抵抗拉伸的刚度会变得非常大，因此在进行仿真时需要增加弹簧的刚度。但仅仅增加弹簧刚度会带来新的问题：对于显式积分而言过大的刚度会导致数值不稳定性，而对于隐式积分则容易造成系统的病态(ill-conditioned)。为了解决这些问题我们往往需要缩短时间步长或是增加迭代次数，但这种做法又增加了计算的复杂度，降低了仿真的效率。
A Single Spring
我们可以通过对弹簧的位移进行约束来解决这样的问题。以单根弹簧为例，假设它的初始长度为\(L\)且刚度无穷大，则弹簧两端的位移需要满足约束：
\[\phi (\mathbf{x}) = \Vert \mathbf{x}_i - \mathbf{x}_j \Vert - L = 0\]当弹簧不满足约束时我们需要调整两端的位置来保证弹簧的长度仍然等于原长\(L\)：
这样我们就可以构造一个优化问题：
\[\begin{aligned} \{ \mathbf{x}_i^\text{new}, \mathbf{x}_j^\text{new} \} &= \arg \min \frac{1}{2} \bigg( m_i \Vert \mathbf{x}_i^\text{new} - \mathbf{x}_i \Vert^2 + m_j \Vert \mathbf{x}_j^\text{new} - \mathbf{x}_j \Vert^2 \bigg) \\ \text{s. t.} \ \ \phi (\mathbf{x}) &= 0 \end{aligned}\]从几何的角度上来看，求解这个优化问题相当于把点\(\mathbf{x} = \{ \mathbf{x}_i, \mathbf{x}_j \}\)投影到满足约束的区域\(\Omega = \{ \mathbf{x}: \phi(\mathbf{x}) = 0 \}\)上，并且使得投影过程的能量消耗最小。
求解这个优化问题可以得到更新后的端点位置：
\[\mathbf{x}_i^\text{new} = \mathbf{x}_i - \frac{m_j}{m_i + m_j} (\Vert \mathbf{x}_i - \mathbf{x}_j \Vert - L) \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert}\] \[\mathbf{x}_j^\text{new} = \mathbf{x}_j + \frac{m_i}{m_i + m_j} (\Vert \mathbf{x}_i - \mathbf{x}_j \Vert - L) \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert}\]Multiple Springs
当需要处理多根弹簧组成的系统时，我们需要通过迭代的方式来求解。最简单的方法是对于每根弹簧单独求解，然后依次更新弹簧两端的位置并不断迭代，这种求解方法称为Gauss-Seidel方法(Gauss-Seidel Approach)：
通常情况下我们无法保证迭代后一定能够满足所有的约束，但随着迭代次数的增加系统的约束情况会不断改善。需要注意的是Gauss-Seidel方法的效果与弹簧的顺序有很大的关系，弹簧的顺序不仅会影响收敛速度，还会影响最终的优化结果。
我们可以使用Jacobi方法(Jacobi Approach)来改善Gauss-Seidel方法依赖于弹簧顺序的问题。在Jacobi方法中会先保存每个顶点的更新量，然后对所有更新量取平均再更新弹簧顶点的位置：
Jacobi方法的缺点在于它的收敛速度要小于Gauss-Seidel方法，但它不依赖于弹簧的计算顺序而且容易并行，因此Jacobi方法的应用更为广泛。
Position Based Dynamics
在此基础上我们就可以设计算法来解决布料的仿真问题。我们首先来介绍position based dynamics(PBD)方法，它的做法非常直观：我们首先按照一般的质点系统来更新质点的位置和速度，然后利用Gauss-Seidel方法或是Jacobi方法考虑系统的约束，最后利用更新后的位置来更新质点的运动状态：
PBD的缺陷在于它不是一个完全基于物理的方法，整个系统的弹性实际上是由迭代的次数以及网格的密度来决定：迭代次数越多、网格越密系统的弹性就越大。
PBD算法的优势在于它容易实现且易于并行，在网格分辨率比较低的情况下有非常快的收敛速度，而且除了布料模拟外也可以应用在其他类型的约束上；但它的缺点在于它不是基于物理的方法，因此得到的不是一个完全准确的解，同时在网格分辨率较高时计算效率会比较低。
Strain Limiting
为了提高PBD的计算效率，我们可以在投影前的仿真计算中先施加一定的约束来防止出现过大的变形，这样的做法称为strain limiting。对于单根弹簧，我们在仿真时施加限制：
\[\sigma^\text{min} \leq \frac{1}{L} \Vert \mathbf{x}_i - \mathbf{x}_j \Vert \leq \sigma^\text{max}\]不难发现此时我们的约束条件与PBD计算投影时的约束非常相似，只不过是把等式约束放宽为了不等式约束。求解优化问题可以得到：
\[\sigma = \frac{1}{L} \Vert \mathbf{x}_i - \mathbf{x}_j \Vert\] \[\sigma_0 = \min (\max (\sigma, \sigma^\text{min}), \sigma^\text{max})\] \[\mathbf{x}_i^\text{new} = \mathbf{x}_i - \frac{m_j}{m_i + m_j} (\Vert \mathbf{x}_i - \mathbf{x}_j \Vert - \sigma_0 L) \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert}\] \[\mathbf{x}_j^\text{new} = \mathbf{x}_j + \frac{m_i}{m_i + m_j} (\Vert \mathbf{x}_i - \mathbf{x}_j \Vert - \sigma_0 L) \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert}\]在弹簧的基础上我们可以在三角形上定义类似的优化问题：
\[\{ \mathbf{x}_i^\text{new}, \mathbf{x}_j^\text{new}, \mathbf{x}_k^\text{new} \} = \arg \min \frac{1}{2} \bigg( m_i \Vert \mathbf{x}_i^\text{new} - \mathbf{x}_i \Vert^2 + m_j \Vert \mathbf{x}_j^\text{new} - \mathbf{x}_j \Vert^2 + m_k \Vert \mathbf{x}_k^\text{new} - \mathbf{x}_k \Vert^2 \bigg) \\\]求解后可以发现更新后的三角形只是在原始三角形的位置上进行了缩放：
strain limiting在物理仿真中有非常丰富的应用。它可以提高仿真过程中的数值稳定性，同时避免出现大变形的问题。除此之外strain limiting很适合描述布料在形变时的非线性行为，甚至对布料的锁死问题也有一定的缓解。
Projective Dynamics
projective dynamics(PD)与PBD非常相似，但在PD中我们不会直接更新质点的位置，而是利用位置构造出系统能量：
\[E(\mathbf{x}) = \sum_{e=\{i, j\}} \frac{k}{2} \Vert (\mathbf{x}_i - \mathbf{x}_j) - (\mathbf{x}_{e, i}^\text{new} - \mathbf{x}_{e, j}^\text{new}) \Vert^2\] \[\{ \mathbf{x}_{e, i}^\text{new}, \mathbf{x}_{e, j}^\text{new} \} = \text{Projection}_e (\mathbf{x}_i , \mathbf{x}_j)\]上式的意义是直接令投影后的质点位置\(\mathbf{x}^\text{new}\)满足约束从而构造出能量\(E\)。由于投影后弹簧两端的质点仍然在弹簧定义的直线上，因此可以对能量项进行变形得到：
\[\begin{aligned} E(\mathbf{x}) &= \sum_{e=\{i, j\}} \frac{k}{2} \Vert (\mathbf{x}_i - \mathbf{x}_j) - (\mathbf{x}_{e, i}^\text{new} - \mathbf{x}_{e, j}^\text{new}) \Vert^2 \\ &= \sum_{e=\{i, j\}} \frac{k}{2} \bigg\Vert (\mathbf{x}_i - \mathbf{x}_j) - L_e \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert} \bigg\Vert^2 \\ &= \sum_{e=\{i, j\}} \frac{k}{2} \bigg\Vert \Vert \mathbf{x}_i - \mathbf{x}_j \Vert \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert} - L_e \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert} \bigg\Vert^2 \\ &= \sum_{e=\{i, j\}} \frac{k}{2} ( \Vert \mathbf{x}_i - \mathbf{x}_j \Vert - L_e )^2 \end{aligned}\]其中\(L_e\)为弹簧变形后的长度，不难发现此时系统的能量与弹簧系统的能量是相同的。假设更新后的位置与初始状态无关，对能量项求偏导可以得到质点受到的外力：
\[\begin{aligned} \mathbf{f}_i &= -k \frac{\partial E}{\partial \mathbf{x}_i} = - \sum_{e=\{i, j\}} (\mathbf{x}_i - \mathbf{x}_j) - (\mathbf{x}_{e, i}^\text{new} - \mathbf{x}_{e, j}^\text{new}) \\ &= -k \sum_{e=\{i, j\}} (\mathbf{x}_i - \mathbf{x}_j) - L_e \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert} \\ &= -k \sum_{e=\{i, j\}} ( \Vert \mathbf{x}_i - \mathbf{x}_j \Vert - L_e ) \frac{\mathbf{x}_i - \mathbf{x}_j}{\Vert \mathbf{x}_i - \mathbf{x}_j \Vert} \end{aligned}\]不难发现此时每个质点收到的外力与推导弹簧系统时弹簧端点受到的外力完全相同。但需要注意的是如此定义的系统能量和直接求解弹簧系统相比具有不同的Hessian矩阵，以下图定义的弹簧网格为例，此时系统的Hessian矩阵为常数：
因此使用PD进行仿真时我们的计算流程与隐式积分完全相同，只需要使用固定的Hessian矩阵来代替原始系统的Hessian矩阵即可：
使用PD来代替隐式积分的优势在于此时计算状态更新时的系数矩阵\(\mathbf{A} = \frac{1}{\Delta t^2} \mathbf{M} + \mathbf{H}\)始终为常数，因此我们可以在计算更新量时通过事先的矩阵分解来直接计算更新量。从解线性系统的角度上讲PD的本质是使用一个常数矩阵来近似系统真实的Hessian矩阵，当这个近似的矩阵有比较好的性质(正定)时我们仍然可以得到线性系统的近似解。
Constrained Dynamics
PBD和PD方法的局限性在于它们难以处理系统刚度非常大的情况，而且如果想要严格地保证约束成立就需要进行非常多次的迭代。constrained dynamics则是一种可以严格保证约束成立的仿真技术，对于弹簧质点系统当弹簧刚度无穷大时系统的约束可以表示为：
\[\phi_e (\mathbf{x}) = \Vert \mathbf{x}_{e,i} - \mathbf{x}_{e,j} \Vert - L_e\]我们可以把系统的能量改写成关于约束的函数：
\[\begin{aligned} E(\mathbf{x}) &= \sum_{e=\{i, j\}} \frac{k}{2} ( \Vert \mathbf{x}_i - \mathbf{x}_j \Vert - L_e )^2 \\ &= \frac{1}{2} \boldsymbol{\phi}^T \mathbf{C}^{-1} \boldsymbol{\phi} \end{aligned}\]其中\(\mathbf{C} = \frac{1}{k} \mathbf{I}\)称为系统的柔度矩阵(compliant matrix)。对能量求偏导得到质点受到的力：
\[\mathbf{f} (\mathbf{x}) = -\nabla E = - \bigg( \frac{\partial E}{\partial \boldsymbol{\phi}} \frac{\partial \boldsymbol{\phi}}{\partial \mathbf{x}} \bigg)^T = - \mathbf{J}^T \mathbf{C}^{-1} \boldsymbol{\phi} = \mathbf{J}^T \boldsymbol{\lambda}\]其中\(\mathbf{J}\)为约束关于质点位置的Jacobi矩阵，\(\boldsymbol{\lambda} = -\mathbf{C}^{-1} \boldsymbol{\phi}\)称为Lagrangian乘子。接下来我们对隐式积分公式进行变形：
\[\mathbf{v}^{\text{new}} = \mathbf{v} + \Delta t \mathbf{M}^{-1} \mathbf{f}^{\text{new}} \iff \mathbf{M} \mathbf{v}^{\text{new}} - \Delta t \mathbf{f}^{\text{new}} = \mathbf{M} \mathbf{v}\]不难发现隐式积分实际上就是动量守恒。接着带入力的公式得到：
\[\mathbf{M} \mathbf{v}^{\text{new}} - \Delta t \mathbf{J}^T \boldsymbol{\lambda}^{\text{new}} = \mathbf{M} \mathbf{v}\]同时，利用\(\boldsymbol{\lambda}\)的定义有：
\[\mathbf{C} \boldsymbol{\lambda}^{\text{new}} = - \boldsymbol{\phi}^{\text{new}} \approx - \boldsymbol{\phi} - \mathbf{J} (\mathbf{x}^{\text{new}} - \mathbf{x}) \approx - \boldsymbol{\phi} - \Delta t \mathbf{J} \mathbf{v}^{\text{new}}\]将两个式子联立起来可以得到矩阵形式的方程组：
\[\begin{bmatrix} \mathbf{M} & -\Delta t \mathbf{J}^T \\ \Delta t \mathbf{J}^T & \mathbf{C} \end{bmatrix} \begin{bmatrix} \mathbf{v}^{\text{new}} \\ \boldsymbol{\lambda}^{\text{new}} \end{bmatrix} = \begin{bmatrix} \mathbf{M} \mathbf{v} \\ - \boldsymbol{\phi} \end{bmatrix}\]因此我们可以通过求解这个方程组来实现系统状态的更新，常用的求解方法包括直接解或是利用Schur补来进行消元求解。
constrained dynamics最大的优势在于它可以处理系统刚度无穷大的情况，此时只需要将\(\mathbf{C}\)矩阵设置为零矩阵即可。