论文笔记_SIGGRAPH2019会前课程：An Introduction to Physics-Based Animation_3

VI. Constraints

discuss the use of constraints约束 in physics-based animation

• the difference between bilateral 双边的(equality) and unilateral 单边的(inequality) constraints and
  the resulting implications由此产生的影响
• discuss the difference between soft and hard constraints
• penalty forces, Lagrange multipliers拉格朗日乘数, and generalized coordinates广义坐标
• how constraints are enforced in practical constrained rigid body systems
• non-penetration constraints非穿透性约束（result from collision碰撞 and contact接触）

除了牛顿第二定律描述的物体在力的作用下的运动，在更多的场景中定义物体的位置或速度需要满足的条件比定义作用在物体上的力更自然和方便，这些条件就是约束。
e.g.静置在桌子上的物体可能受重力或其他一些力，约束它不要穿透桌子，物理上桌子给一个支持力，与物体给桌子的力大小相等方向相反，最终对物体的作用效果保证了不穿透，与非穿透约束相关的这个力就正是维持这个约束必需的力。
在更多的情景中，约束力是系统中为了维持约束即时产生的对其他力的反应，定义约束和必须的力比定义约束力更自然和方便。
约束有很多种，约束动力学问题也有很多种解法。

5.1 Bilateral and Unilateral Constraints

自由度degrees of freedom (DOF)，是系统中完全定义系统配置的独立参数的数量。e.g.三维空间，单个粒子有三个自由度，单个刚体有六个自由度。

双边约束/等式约束

一个完整的约束是系统中位置自由度的关系，形式为：$g(x_1,x_2,...,x_n,t)=0$
它可能取决于时间。由于这个关系必须被满足，它就会隐式地去掉了一些自由度。因此若n个中有n-1个已知，那另一个就可以由上边的等式定义。
上式定义了一个n维空间的manifold流形（流形是高维空间中曲线、曲面概念的拓广），这个约束是双边或等式约束，使得位置恰在流形上，而非任一边。
e.g.二维空间的两个刚体，两端相连，可关于接触点旋转，如图，刚体自由度为${\mathbf{x}}_1,{\theta }_1,{\mathbf{x}}_2,{\theta }_2$，其约束:
$g({\mathbf{x}}_1,{\theta }_1,{\mathbf{x}}_2,{\theta }_2)=({\mathbf{x}}_1+{\mathbf{r}}_1)-({\mathbf{x}}_2+{\mathbf{r}}_2)=0$
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-tEXJNctS-1616637060925)(5.%E7%BA%A6%E6%9D%9F_md_files/image_20200704010003.png?v=1&type=image&token=V1:kx0El9b2SZo3CWwnotbBRyiCjfs5TDHkp1VqmvGnY2E)]
𝐠 ∶ ℝ6 →ℝ2,隐式消除了两个自由度，故此约束有6-2=4个自由度。也可以认为是整个系统在二维空间的位置有两个自由度，加上两刚体各自的角度共四个自由度。

单边约束/不等式约束

形式：$g(x_1,x_2,...,x_n,t)\ge 0$
物体与另一物体相接触会产生这种形式的约束。这种情况下物体可以分开但不可以穿透。
如图，球与平面接触，可在平面上运动，可从平面上拿起，蓝色点为接触点，$y_1$是球心的高，约束为：
$g(x_1,y_1,{\theta }_1)=y_1-h\ge 0$
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-whAKGRtM-1616637060926)(5.%E7%BA%A6%E6%9D%9F_md_files/image_20200704011525.png?v=1&type=image&token=V1:SMsF1Rpn8dLu-pCZGcsU0aoZh50Rz0eca77xuA3DXhs)]
球可以自由滚动，侧面位置$x_1$和朝向${\theta }_1$是独立的，若进一步限制旋转，则在接触中这两个量不会再独立变化。
其他的单侧约束的情况，比如限定铰接的刚体连接角的法定范围。

5.2 Soft vs. Hard Constraints

hard constraint

必须恰好满足。维持它的力必须克服任何违反约束的抵抗的力。

soft constraint

不需要恰好满足。可以容忍一定的错误。维持它的力不像前者那么强横，约束的抵抗力可以使其被违反，直到这种力消退。

5.3 Solution Methods: Penalty Forces, Lagrange Multipliers,Generalized Coordinates

约束动力学问题的仿真的几种方式。被广泛地分为两类：使用maximal coordinates,最大的坐标（这种坐标无视约束，利用辅助条件定义约束）；使用generalized coordinates广义坐标（有时指减少的的坐标，可能因为与最大坐标相比，它们的数量减少了），其参数代表计算约束后的自由度，由此可推出最大坐标。
Penalty methods罚函数法和拉格朗日乘数法是使用最大坐标的例子。

Penalty Methods

通过使用类似弹簧的恢复的力量来惩罚约束违反以维持约束。
e.g.考虑下图场景，用弹簧连接两刚体，弹簧静息长度为0，刚体1在接触点受到的力：
$f_{21}({\mathbf{x}}_1,{\theta }_1,{\mathbf{x}}_2,{\theta }_2)=-k(({\mathbf{x}}_2+{\mathbf{r}}_2({\theta }_2))-({\mathbf{x}}_1+{\mathbf{r}}_1({\theta }_1)))$
刚体2：$f_{12}({\mathbf{x}}_1,{\theta }_1,{\mathbf{x}}_2,{\theta }_2)=-f_{21}({\mathbf{x}}_1,{\theta }_1,{\mathbf{x}}_2,{\theta }_2)$
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-DrWXcUho-1616637060927)(5.%E7%BA%A6%E6%9D%9F_md_files/image_20200704014434.png?v=1&type=image&token=V1:gXEl57a8LLsunzXcsJNpuJDizCig7Dl9Kpvr-t1zufw)]
惩罚力可以根据与约束的违反相关的能量来表示，详见[Faure et al. 2008]
包含取决于约束的导数的阻尼项的惩罚力，详见[Bender et al. 2014]
约束违反的惩罚项实现相对简单，但存在一些缺陷，详见[Drumwright 2008; Bender et al. 2014].
和前面例子中的参数k相同，惩罚力项包含一个或更多强度参数，多个约束的竞争力的复杂系统中，Tuning stiffness parameters优化刚度参数可能是个挑战。
k较大，可以更好地维持约束。
但这样也向系统中引入了一个僵硬的力，可能会导致严格的时间步长限制或必需的隐式积分。也可能导致不期待的振动，使得静止接触的物体的场景难以实现，详见[Bender et al. 2014]。
尽管如此，罚函数法仍被广泛应用，一些致力于减轻其劣势的研究详见[Harmon et al. 2009; Tang et al. 2012;
Xu et al. 2014]

Lagrange multipliers

在运动方程中显式包含约束力的方式。与惩罚力不同，其约束力就是抵消任何其他破坏约束的竞争力的力。
从数学上讲，这足以维持一个最初满足的约束条件，在实际应用中，需要附加的稳定技术来抵消约束满足中的漂移。
考虑约束的等式：$g(\mathbf{x})=0$,其中$\mathbf{x}=(x_1,x_2,...,x_n)$,这是一个定义了所有合法的位置为x的隐式曲面方程。约束的梯度为向量
$▽g(\mathbf{x})=\left[\begin{array}{c}\frac{\partial g}{\partial x_1}(\mathbf{x})\\ \frac{\partial g}{\partial x_2}(\mathbf{x})\\ ...\\ \frac{\partial g}{\partial x_n}(\mathbf{x})\end{array}\right]$
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-r7Dfdaka-1616637060928)(5.%E7%BA%A6%E6%9D%9F_md_files/image_20200704230508.png?v=1&type=image&token=V1:i2W8xtvHriAbA37EyIwd3nIqTypDLpK7FzuCTLlYWqA )]

如图，梯度向量垂直于面，向量z与约束曲面相切与点x，满足${\mathbf{z}}^T▽g(\mathbf{x})=0$
考虑任意力$\mathbf{F}$作用于有约束的粒子上,当它欲使粒子离开约束曲面时会产生一个与之相抵消的约束力${\mathbf{F}}_c$，除了摩擦力，约束力仅作用在$▽g(\mathbf{x})$方向，因此workless,详见[Lanczos 2012]。考虑一个虚拟位移$\delta \mathbf{r}$,与约束一致，${\mathbf{F}}_c^T\delta \mathbf{r}=0$,故$\delta \mathbf{r}=\alpha \mathbf{z}$,那么$\delta {\mathbf{r}}^T▽g(\mathbf{x})=0$,所以${\mathbf{F}}_c=\lambda ▽g$,$\lambda$被称为约束的拉格朗日乘数，是标量，决定了约束力的强度。
更多的情况下，考虑m个约束,等式为：
$g(\mathbf{x})=\left[\begin{array}{c}{g}_1(\mathbf{x})\\ g_2(\mathbf{x})\\ ...\\ g_m(\mathbf{x})\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ ...\\ 0\end{array}\right]=0$
约束的雅可比矩阵：
$J(\mathbf{x})=\left[\begin{array}{cccc}\frac{\partial g_1}{\partial x_1}(\mathbf{x})& \frac{\partial g_1}{\partial x_2}(\mathbf{x})& ...& \frac{\partial g_1}{\partial x_n}(\mathbf{x})\\ \frac{\partial g_2}{\partial x_1}(\mathbf{x})& \frac{\partial g_2}{\partial x_2}(\mathbf{x})& ...& \frac{\partial g_2}{\partial x_n}(\mathbf{x})\\ ...& & & \\ \frac{\partial g_m}{\partial x_1}(\mathbf{x})& \frac{\partial g_m}{\partial x_2}(\mathbf{x})& ...& \frac{\partial g_m}{\partial x_n}(\mathbf{x})\end{array}\right]=\left[\begin{array}{c}▽g_1(\mathbf{x}{)}^T\\ ▽g_2(\mathbf{x}{)}^T\\ ...\\ ▽g_m(\mathbf{x}{)}^T\end{array}\right]$
那么总的约束力由m项组成，
${\mathbf{F}}_c={\lambda }_1▽g_1+{\lambda }_2▽g_2+...+{\lambda }_m▽g_m=J^T\lambda$
这里，$\lambda =({\lambda }_1,...{\lambda }_m)$
包含作用力和约束力的运动公式写为：
$M\mathbf{a}=\mathbf{F}+{\mathbf{F}}_c=\mathbf{F}+J^T\lambda$
进行类似[Baraff 1996; Witkin 2001]的推导，未知数$\lambda$可以通过要求系统加速度与约束兼容来定义。约束等式两边关于时间t作微分；$0=\dot{g}(\mathbf{x})=\frac{\partial g^T}{\partial \mathbf{x}}\dot{\mathbf{x}}(t)=J(t)\mathbf{v}(t)$
这表明。合法的速度正交于约束的梯度，平行于上图中的向量z。
再关于t进行一次微分：$0=\ddot{g}(\mathbf{x})=\dot{J}\mathbf{v}(t)+J\dot{\mathbf{v}}(t)=\dot{J}\mathbf{v}(t)+J\mathbf{a}(t)$
这是加速度必须满足的条件，以保持位置和速度水平约束。
从运动公式得出加速度并代入：
$J{M}^{-1}(\mathbf{F}+J^T\lambda )=-\dot{J}\mathbf{v}(t)$,
把带$\lambda$的项单独放左边：$(J{M}^{-1}{J}^T)\lambda =-J{M}^{-1}\mathbf{F}-\dot{J}\mathbf{v}(t)$
实际应用中，约束往往只耦合系统中的几个自由度，因此J和M一般是稀疏矩阵，但$(J{M}^{-1}{J}^T)$一般不稀疏。
重写上式应用于稀疏KKT系统：
$\left[\begin{array}{cc}M& -J^T\\ -J& 0\end{array}\right]\left[\begin{array}{c}\mathbf{y}\\ \lambda \end{array}\right]=\left[\begin{array}{c}0\\ -\mathbf{b}\end{array}\right]$
此处，$\mathbf{y}=\mathbf{a}-M^{-1}\mathbf{F}$,$\mathbf{b}=-J{M}^{-1}\mathbf{F}-\dot{J}\mathbf{v}(t)$
为为它的解决方案描述了一个线性时间算法。应用于许多实际情景。
KKT系统常出现于各种约束动力学问题的解决方法中。

Generalized coordinates

约束动力学的另一种方法是基于寻找独立坐标的最小集——广义坐标$q_1,...,q_N$,用这些坐标和它们的导数来表示运动方程。用广义坐标表示的运动方程是变分力学的一部分，变分力学由欧拉、拉格朗日、汉密尔顿等人在牛顿之后发展而来，详见[Lanczos 2012]。$q_i$的选择不唯一，但有的情况可能有自然而方便的某个选择。
广义坐标的一个最出名的例子是：一个刚体可以看作是粒子的集合，每一对粒子被约束保持固定的距离。三维空间中6个自由度（3个位置和3个朝向）的最小集用以表示一个刚体，物体上任意点的最大坐标可以用2.2.1节中已知的公式由广义坐标计算出来。
n（n>N）个最大坐标作为广义坐标的函数被给出:
$x_1=x_1(q_1,q_2,...,q_N)$
$x_2=x_2(q_1,q_2,...,q_N)$
…
$x_n=x_n(q_1,q_2,...,q_N)$
结合所有的极大坐标和广义坐标，定义这个$\mathbf{x}=\mathbf{x}(\mathbf{q})$,这个方程组也可能取决于时间，详见[Murray et al. 1986; Lanczos 2012]，此处不做讨论。
使用$\mathbf{q}$或者$\dot{\mathbf{q}}$,这些表达式可用于定义系统的动能和势能、广义的速度。
系统的速度为$\dot{\mathbf{x}}(t)=\frac{\partial \mathbf{x}}{\partial \mathbf{q}}\dot{\mathbf{q}}$
令$J(\mathbf{q})=\frac{\partial \mathbf{x}}{\partial \mathbf{q}}$为坐标函数的雅可比矩阵，
动能为：$T=\frac{1}{2}{\dot{\mathbf{x}}}^{T}M\dot{\mathbf{x}}=\frac{1}{2}(J\dot{\mathbf{q}}{)}^{T}MJ\dot{\mathbf{q}}=\frac{1}{2}{\dot{\mathbf{q}}}^T{J}^{T}MJ\dot{\mathbf{q}}$
类似于极大坐标系中的力，可以定义广义坐标系中的广义力，例如刚体中的广义力是力的六个组成部分和作用在物体上的力矩，如果两组坐标下的两组力在符合约束条件的虚位移上做相同的虚功，那么将产生等效的动力学，详见[Lanczos 2012]。故任何系统中的广义力可通过使两种力在与约束条件一致的虚位移（virtual displacements）上所作的虚功（virtual work）相等得出，即：$\mathbf{F}\cdot \delta \mathbf{x}=\mathbf{G}\cdot \delta \mathbf{q}$.
这里$\mathbf{G}\in$ℝ${}^N$,表示广义力，位移$\delta \mathbf{x}$服从约束，$\delta \mathbf{x}=J\delta \mathbf{q}$。由此可得：
${\mathbf{F}}^T\delta \mathbf{x}={\mathbf{G}}^T\delta \mathbf{q}$;
${\mathbf{F}}^{T}J\delta \mathbf{x}={\mathbf{G}}^T\delta \mathbf{q}$
由于坐标$\delta \mathbf{q}$不受约束，是任意的，那就必然有：
${\mathbf{F}}^{T}J={\mathbf{G}}^T$或$\mathbf{G}=J^T\mathbf{F}$

验证这个结果——
对于刚体，粒子速度：$\dot{\mathbf{x}}=\left[\begin{array}{cc}I& {\mathbf{r}}^{\ast T}\end{array}\right]\left[\begin{array}{c}\mathbf{v}\\ \omega \end{array}\right]$.
所以$J=\left[\begin{array}{cc}I& {\mathbf{r}}^{\ast T}\end{array}\right]$.
力F作用于粒子上，得到广义力$\mathbf{G}=J^T\mathbf{F}=\left[\begin{array}{c}I\\ {\mathbf{r}}^{\ast }\end{array}\right]\mathbf{F}=\left[\begin{array}{c}\mathbf{F}\\ \mathbf{r}\times \mathbf{F}\end{array}\right]$。(类似于力和力矩的表达式)

拉格朗日运动方程，类似于牛顿第二定律，可写为：$\frac{d}{dt}(\frac{\partial T}{\partial \dot{\mathbf{q}}})=\mathbf{G}$;
$\mathbf{G}$是广义力。定义拉格朗日算符$L=T-V$，V是和力有关的势能，将上式表示为欧拉-拉格朗日方程
$\frac{d}{dt}(\frac{\partial L}{\partial \dot{\mathbf{q}}})-\frac{\partial L}{\partial \mathbf{q}}=0$。
附加的力，考虑与广义坐标无关的非保守力或约束力可以被包括在方程的右侧。
e.g.考虑如下图所示的一个二维的摆。单个粒子质量m，最大坐标系下坐标 $\mathbf{x}=(x，y)$,位置的坐标有两个，但由于摆的约束，系统自由度为1，使用摆和垂直方向的角度$\theta$作为广义坐标，将连接点作为原点，最大坐标为：
$x(\theta )=lsin\theta$;
$y(\theta )=-lcos\theta$。
速度为：$\left[\begin{array}{c}\dot{x}(\theta )\\ \dot{y}(\theta )\end{array}\right]=\left[\begin{array}{c}lcos\theta \dot{\theta }\\ lsin\theta \dot{\theta }\end{array}\right]=J\dot{\theta }$;
动能：$T=\frac{1}{2}m({\dot{x}}^2+{\dot{y}}^2)=\frac{1}{2}m{l}^2{\dot{\theta }}^2$.
那么$\frac{\partial T}{\partial \dot{\theta }}=m{l}^2\dot{\theta }$,则$\frac{d}{dt}\frac{\partial T}{\partial \dot{\theta }}=m{l}^2\ddot{\theta }$。
唯一需要明确解释的力是重力，在极大坐标系下表示为$(0,-mg{)}^T$,由此，广义力为：$J^T\left[\begin{array}{c}0\\ -mg\end{array}\right]=\left[\begin{array}{cc}lcos\theta & lsin\theta \end{array}\right]\left[\begin{array}{c}0\\ -mg\end{array}\right]=-lmgsin\theta$
将此式上一个等式合并，得到$m{l}^2\ddot{\theta }=-lmgsin\theta$或者著名的钟摆运动方程$\ddot{\theta }+\frac{g}{l}sin\theta =0$

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-oWBmd58H-1616637060929)(5.%E7%BA%A6%E6%9D%9F_md_files/image_20200705183325.png?v=1&type=image&token=V1:JCbdKCRa8Jy5Nlw6ULk4yfbErNifXRW5NNXxlG_lOqE)]

5.4 Practical Constrained Rigid Body Systems 实用的约束刚体系统

实际的刚体仿真系统必须处理各种约束条件，包括关节约束，或连接(如球和球座、转动关节、棱柱)，碰撞，静止和滑动接触。布娃娃物理系统(ragdoll physics),对于关节/衔接约束规定的连通性是树的连通性(没有闭环)，广义坐标通常被有效使用，详见[Bender et al. 2014]，比如Featherstone’s algorithm，详见[Featherstone 1983]。对于其他的约束，约束力被广泛应用，对拉格朗日乘数的公式进行了一些修正，在5.3节阐述。

5.3节中使用约束方程的二阶导，给出了约束满足的加速级条件。然后通过将由牛顿第二定律得出的加速度带代入到这个条件，对拉格朗日乘数的线性系统求导，并且解出$\lambda$.有效决定了维持合法加速所需的力.假设位置和速度都已合法，合法的加速度会随时间保持位置和速度水平的约束。在实践中，加速水平约束有一些话题，使其不适合作为刚体系统的通用方法。牛顿第二定律描述了速度的连续变化，然而刚体的碰撞会导致速度的不连续变化，比如一个坚硬的石块撞击地面，然后立即弹起。在模拟摩擦接触场景时，约束条件的加速水平解法也存在问题，比如把一个物体推过一个粗糙的表面(就像Painleve悖论一样)。

基于这些原因，很多系统在速度水平模拟刚体约束，使用描述有限速度变化的牛顿方程积分，详见[Mirtich 1996; Guendelman et al. 2003; Weinstein et al. 2006]。速度水平的公式中，冲量不连续地改变物体的速度。冲量定义为力对时间的积分：
$\mathbf{j}={\int }_{t_1}^{t_2}\mathbf{f}dt$.

因此它有动量的单位,$N\cdot s=kgm/s$,一个脉冲作用于质量m的粒子会导致粒子动量的变化。质量为m的粒子关于时间对牛顿第二定律积分，得到动量-脉冲方程：
${\int }_{t_1}^{t_2}m\mathbf{a}dt={\int }_{t_1}^{t_2}\mathbf{f}dt$
$=>m(\mathbf{v}(t_2)-v(t_1))=\mathbf{j}$

为了在速度水平重写应用于稀疏KKT系统的：
$\left[\begin{array}{cc}M& -J^T\\ -J& 0\end{array}\right]\left[\begin{array}{c}\mathbf{y}\\ \lambda \end{array}\right]=\left[\begin{array}{c}0\\ -\mathbf{b}\end{array}\right]$。（$\mathbf{y}=\mathbf{a}-M^{-1}\mathbf{F}$,$\mathbf{b}=-J{M}^{-1}\mathbf{F}-\dot{J}\mathbf{v}(t)$）
考虑离散化牛顿定律：$M{\mathbf{V}}^{n+1}=M{\mathbf{V}}^n+\Delta t\mathbf{F}+\Delta t{J}^T{\lambda }^{n+1}$。
速度的约束为：$J{\mathbf{V}}^{n+1}=0$,定义${\mu }^{n+1}=\Delta t{\lambda }^{n+1}$,将未知项${\mathbf{V}}^{n+1}$和${\mu }^{n+1}$挪到等式左边，结合两个等式得到线性系统，详见[Tournier et al. 2015]：
$\left[\begin{array}{cc}M& -J^T\\ -J& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}^{n+1}\\ {\mu }^{\mathbf{n}+1}\end{array}\right]=\left[\begin{array}{c}M{\mathbf{V}}^n+\Delta t\mathbf{F}\\ 0\end{array}\right]$.
这个线性系统类似于前边的，但是用单位动量代替了力，用冲量形式的$\mu =\Delta t\lambda$代替$\lambda$作为未知量。从等式中消掉${\mathbf{V}}^{n+1}$,得到：$J{M}^{-1}{J}^T{\mu }^{n+1}=-J{\mathbf{V}}^n-\Delta tJ{M}^{-1}\mathbf{F}$。这是耦合了约束冲量的解的方程，发生在物体速度被调整到适应约束的情况下，可能会导致不同约束间的冲突。必须同时考虑系统中的所有约束以得到全部都满足的解决方案。

实践中很多系统用基于冲量的方法（impulse-based approach，比如[Mirtich 1996; Guendelman et al. 2003; Weinstein et al. 2006]），使用有限次高斯-赛德尔迭代（Gauss-Seidel iterations）计算方程的近似解,在迭代中每个约束使用其他约束的最新计算状态独立地更新。相互作用的应用可能会严格地限制允许迭代的次数，并且需要类似warmstarting the solution暖化解决方案(使用前一个时间步中的解)的一些技术来获得可接受的结果，见[Catto 2014].当一个轻物体与两个重物体相互作用时，高斯-塞德尔方法的收敛速度会非常慢，因为在一次成对迭代中，轻物不会成功地移动重物很多。这种messenger problem可通过同时解一些约束部分缓解，是一种解全局系统和本地成对约束的折中，见[Shinar et al. 2008; Catto 2014].

另一个在实践中产生的问题是由于离散化和数值误差产生的漂移。在5.3节中加速度水平的方程式中，初始假设位置和速度满足约束，之后仅要求加速度合法，经过一段时间会造成速度和位置从合法状态漂移，产生显著的不自然现象，比如错误的交叠或裂缝。本节目前为止介绍的速度水平的方程式有类似的问题，原因是没有明确地纠正约束中的任何位置误差。稳定化技术被用于抑制漂移，见[Bender et al. 2014].一个广泛应用的方法是Baumgarte stabilization，见[Baumgarte 1972],

其加速度约束方程$\ddot{\mathbf{g}}(\mathbf{x})=0$，被换为[Bender et al. 2014]中的：
$\ddot{\mathbf{g}}(\mathbf{x})+2\alpha \dot{\mathbf{g}}(\mathbf{x})+{\beta }^2\mathbf{g}(\mathbf{x})=0$;
参数$\alpha 和\beta$可调整。

对速度水平的方程，其速度约束方程$\dot{\mathbf{g}}(\mathbf{x})=0$被换为：
$\dot{\mathbf{g}}(\mathbf{x})+\gamma \mathbf{g}(\mathbf{x})=0$;
参数$\gamma$可调整。

当位置水平的约束$\mathbf{g}(\mathbf{x})=0$已经满足，会得到和之前相同的方程，如果位置约束存在误差，计算的速度将使误差为零。

许多实际的刚体系统向硬约束添加了依从性，通过将$J\mathbf{v}=0$换为
$J\mathbf{v}+\frac{\beta }{\Delta t}\mathbf{g}(\mathbf{x})+\gamma \lambda =0$.
(e.g., Bullet [Coumans 2010], ODE [Smith et al. 2005], and Box2D [Catto 2011a]);
$\beta 和\gamma$为常量，见[Catto 2011b].线性方程组变为：
$\left[\begin{array}{cc}M& -J^T\\ -J& {\gamma }^I\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}^{n+1}\\ {\mu }^{n+1}\end{array}\right]=\left[\begin{array}{c}M{\mathbf{V}}^n+\Delta t\mathbf{F}\\ -\frac{\beta }{\Delta t}\mathbf{g}({\mathbf{X}}^n)\end{array}\right]$

这样使得约束稍缓和，同时抵制漂移。另外，新的对角块正则化了方程使其更易使用数值方法求解，在冗余或接近冗余约束的情况下，可能制约能力变差，详见[Catto 2011b;Tournier et al. 2015; Lacoursiere 2007]。

5.5 Non-penetration Constraints, Collisions,and Contact

许多刚性或可变形固体的动态场景通常会导致物体间的许多碰撞和接触。Soft bodies还有自我碰撞，即物体的一部分变形并且撞上了另一部分。特别地，布料仿真会经历许多同时发生的自我碰撞，要解决这个问题是很有挑战性的。所有这些情况，求解器的计算成本的显著成分是在检测和解决碰撞和接触（解决碰撞之前要先检测到）。

Collision detection.碰撞检测

通过检查物体几何的相互渗透来检测。碰撞检测算法使用了大量的几何表示。
多边形几何常用于表示刚体和可变形体。许多算法已经发展为寻找多边形的交，特别是凸多边形（内部任意两点连线仍在多边形内部）。

分离轴定理（Separating Axis Theorem）

认为如果在两个凸形状之间可以找到一个分离的超平面，则它们是不相交的，因此不会贯穿。潜在的碰撞几何面自然成为分离超平面的候选。凸分解（Convex decompositions）常用于非凸物体碰撞检测的计算。使用多边形几何要考虑大量相交的情况，例如点-面交叉和边-边交叉。

有向距离场the signed distance field

另一种在碰撞检测中常用的几何表示。物体的有向距离场$\varphi (\mathbf{x})$将从$\mathbf{x}$到物体表面最近的点的有向距离分配到空间$\mathbf{x}$的每个点。物体内部的点$\varphi <0$,物体外部$\varphi >0$,或者反过来，取决于约定。
表面就隐式地定义为$\varphi$的零水平集，点$\mathbf{x}$使$\varphi (\mathbf{x})=0$.通常$\varphi$在包含物体的结构化网格上取样。另一物体表面上的点$\mathbf{p}$可搜索$\varphi (\mathbf{p})$来确定一定时间内它在内部还是外部。尽管内部/外部测试很快，但是由于对象很多，在存储距离字段时可能会有很大的内存开销。尽管在预处理步骤中对于刚体可以计算一次$\varphi$,对于可变形几何，可能需要间歇地重新计算。用带符号的距离表示的另一优势是，解碰撞和接触的表面法线容易使用有向距离场计算出：$\mathbf{n}(\mathbf{x})=▽\varphi (\mathbf{x})$.

加速结构

n个物体有$O(n^2)$对，为使碰撞检测应用于许多相互作用的物体的仿真实践，需要使用加速结构。常用方法是为几何对象创建包围盒bounding volumes，特别是那些可进行快速相交测试的。包围盒可用于快速减去可能相交的集合，原因是两物体不相交，其包围盒也不相交，常用球、axis-aliged bounding boxes轴对齐包围盒、oriented bounding boxes。另外还可以计算包围盒的等级结构，潜在地使得更有效地从大量的候选中挑选交集。

除了包围盒，空间划分还可以通过指示在空间的特定区域中同时存在哪些对象来加速碰撞检测。

连续碰撞检测

在离散碰撞检测中，考虑场景中几何对象的状态处于固定时刻，不考虑轨迹。对于大体积对象往往足够了。但对于薄的物体，如布，或非常小的物体可能在一个时间步长中经历相交，这在时间步长开始或结束的固定时刻时观察起来并不明显。比如，一个小物体被扔到一块布上，可能在时间步开始时完全在布的一边，结束时完全在另一边，但在两个时刻之间物体经过了布料，所以原则上是应该检测并解决碰撞。连续碰撞检测就在检测相交时考虑时间步内的运动轨迹。
求解交叉问题的分析方法在某些情况下存在。也用Ray casting光线投射和数值求解器，比离散碰撞检测更复杂计算量更大，有些系统会通过限制时间步长来避免连续碰撞检测。

Collision response (1D).

碰撞反应中，求解器解决了一个之前检测到的冲突，通常是通过计算和应用一对冲量到碰撞的物体上使其更新的速度是合法的。

e.g.一维空间的两个粒子，质量为$m_1,m_2$，以速度$v_1,v_2$相撞，如下图。根据牛顿第三定律，它们会相互施加大小相等方向相反的力，所以有一对冲量$j和-j$影响粒子动量变化。令碰撞后速度为$v_1^{\prime },v_2^{\prime }$,在系统中加入冲量，得到：
$m_1{v}_1^{\prime }=m_1{v}_1+j$;
$m_2{v}_2^{\prime }=m_2{v}_2-j$.

碰撞是非弹性的，或者说粘性的，碰撞后两粒子共同运动，i.e.$v_1^{\prime }=v_2^{\prime }$.
于是我们得到了三个未知数、三个等式，解出碰撞后的速度并带入等式：
$\frac{1}{m_1}(m_1{v}_1+j)=\frac{1}{m_2}(m_2{v}_2-j)$
解出j:$(\frac{1}{m_1}+\frac{1}{m_2})j=v_2-v_1$;
=>$j=(\frac{1}{m_1}+\frac{1}{m_2}{)}^{-1}(v_2-v_1)$.

如果$v_1=v,v_2=-v$;$m_1=m_2=m$,由上得，$j=-mv$,碰撞后$v_1^{\prime }=v_2^{\prime }=0$;
如果$v_1=v,v_2=-v$;$m_1=3,m_2=1$,由上得，碰撞后$v_1^{\prime }=v_2^{\prime }=\frac{1}{2}v$,则这些粒子继续沿着质量更大的粒子运动的方向，以减小的速度一起运动;

用于弹性碰撞的情况，特别是刚体仿真，粒子碰撞后的相对速度不再为0，而是：
$(v_2^{\prime }-v_1^{\prime })=-ϵ(v_2-v_1)$,
这里$ϵ\in [0,1]$是the coefficient of restitution回弹系数。当$ϵ=0$是就是刚刚非弹性碰撞的例子。带入表达式求解j:
$j=(\frac{1}{m_1}+\frac{1}{m_2}{)}^{-1}(1+ϵ)(v_2-v_1)$

如果$v_1=v,v_2=-v$;$m_1=m_2=m$,对于$ϵ=1$的完全弹性碰撞，$j=-2mv$,$v_1^{\prime }=-v,v_2^{\prime }=v$

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-SjGKygk1-1616637060930)(5.约束_md_files/一维质点碰撞.png)]

Deformable object collisions可变形物体碰撞

从物理上讲，碰撞产生的弹跳是由于碰撞物体在碰撞区域的变形造成的,物体碰撞，压缩，存储弹性能，能量的释放导致碰撞后的反弹。因此，对可变形物体可以使用上面描述的非弹性碰撞方程，将一维的情况直接推广到2、3维：
$m_1{\mathbf{v}}_1^{\prime }=m_1{\mathbf{v}}_1+\mathbf{j}$;
$m_2{\mathbf{v}}_2^{\prime }=m_2{\mathbf{v}}_2-\mathbf{j}$;
${\mathbf{v}}_2^{\prime }={\mathbf{v}}_1^{\prime }$
二维的情况是6个方程和未知数，三维的情况是9个。速度是在碰撞点处的，冲量作用于碰撞点，当物体变形时，可能会有许多碰撞点。

Rigid body collisions

与可变性物体相比，二者物理上碰撞的过程是一样的，但理想刚体不会有任何形变，刚体碰撞的模拟需要不同的方式。
碰撞定律常用于联系碰撞前后的量，见[Chatterjee and Ruina 1998].其推导与上面类似，但还要考虑速度和动量的角分量。在运动的切向方向上的摩擦处理大大复杂化了刚体碰撞。首先考虑粘滞和无摩擦碰撞的简单情况。

根据[Ruina and Pratap 2009],考虑两质量为$m_1,m_2$的刚体的碰撞，惯性张量$I_1,I_2$,线速度${\mathbf{v}}_1,{\mathbf{v}}_2$，角速度${\omega }_1,{\omega }_2$，碰撞后的速度和前面一样用素数（prime）表示.场景如下图，碰撞时与物体相切的平面单位法向量为$\mathbf{n}$。
碰撞前后的线动量为：
$m_1{\mathbf{v}}_{C_1}^{\prime }=m_1{\mathbf{v}}_{C_1}+\mathbf{j}$ ;
$m_2{\mathbf{v}}_{C_2}^{\prime }=m_2{\mathbf{v}}_{C_2}-\mathbf{j}$.
角动量：
$I_1{\omega }_1^{\prime }=I_1{\omega }_1+{\mathbf{r}}_1\times \mathbf{j}$ ;
$I_2{\omega }_2^{\prime }=I_2{\omega }_2-{\mathbf{r}}_2\times \mathbf{j}$

对于粘性碰撞，其他方程表明物体在碰撞后在接触点立即以相同的速度运动:
${\mathbf{v}}_{P1}^{\prime }={\mathbf{v}}_{P2}^{\prime }$.
这里，
${\mathbf{v}}_{P1}^{\prime }={\mathbf{v}}_{C1}^{\prime }+{\omega }_1^{\prime }\times {\mathbf{r}}_1$ ;
${\mathbf{v}}_{P2}^{\prime }={\mathbf{v}}_{C2}^{\prime }+{\omega }_2^{\prime }\times {\mathbf{r}}_2$ .
是第一个和第二个物体在碰撞点各自的速度。

方程/未知数个数，二维为8，三维为15.

对于无摩擦弹性碰撞，设$\mathbf{j}=j\mathbf{n}$，
$({\mathbf{v}}_{P2}^{\prime }-{\mathbf{v}}_{P1}^{\prime })\cdot \mathbf{n}=-ϵ({\mathbf{v}}_{P2}-{\mathbf{v}}_{P1})\cdot \mathbf{n}$.
这里，$ϵ\in [0,1]$。已知物体在碰撞前的线速度和角速度、碰撞平面，可解方程确定碰撞后物体的法向冲量和速度。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Iwb65sGY-1616637060931)(5.约束_md_files/刚体碰撞.png)]

Frictional collisions摩擦碰撞

在摩擦碰撞中,冲量$\mathbf{j}$由法向分量和切向分量组成,分别阻止或抑制物体碰撞点的相对切向速度。
Coulomb friction model库伦摩擦模型
描述在粘着或滑动接触中产生的摩擦力的经验模型。切向力的大小不能超过$\mu \Vert {\mathbf{f}}_n\Vert$,接触点黏着的情况下可能会更小。切向力的方向与运动的切向方向正好相反，系数$\mu$为摩擦系数。

滑动接触，切向力处于其可能值的最大值：
${\mathbf{f}}_t=-\mu \Vert {\mathbf{f}}_n\Vert \mathbf{t}$。
这里，$\mathbf{t}$是运动方向。

当接触点粘着时，切向力刚好足以抵消其他切向力。

总结：接触点相对速度为正，摩擦力为$-\mu f_n\mathbf{t}$;相对速度为负，摩擦力为$\mu f_n\mathbf{t}$，为0，介于二者之间。（这种情景改编自[Ruina and Pratap 2009]）

库伦摩擦模型力和冲量的约束可以用the friction cone摩擦锥来表示，定义不等式:
$\Vert {\mathbf{f}}_t\Vert \le \mu \Vert {\mathbf{f}}_n\Vert$.
如下图：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Jo1g7SAg-1616637060932)(5.约束_md_files/摩擦锥.png)]

在[Guendelman et al. 2003]中,通过首先假定在切向方向粘滞向刚体碰撞中引入切向摩擦冲量，即：
$(I-\mathbf{n}{\mathbf{n}}^T)({\mathbf{v}}_{P2}^{\prime }-{\mathbf{v}}_{P1}^{\prime })=0$.
将它与前面的等式一同解出，若$\mathbf{j}$的解在摩擦锥中，那就接受粘滞的假设，使用冲量$\mathbf{j}$；若不在，则施加滑动摩擦。在实践中，类似这个的代数碰撞定律可以产生一些非物理的结果，对几个碰撞定律更详细的分析见[Chatterjee and Ruina 1998]