在现实世界中，几乎没有物体是完全自由的。机械系统通过铰链连接，车轮沿着轨道滚动，机器人的关节限制了运动的自由度。这些限制——我们称之为约束——是物理仿真中最具挑战性也是最重要的课题之一。本章将深入探讨如何在数值仿真中准确、稳定、高效地处理各种约束，这对于构建可靠的机器人仿真环境至关重要。
约束动力学的核心挑战在于：如何在保持系统物理真实性的同时，处理数值积分带来的约束违背？这个看似简单的问题，实际上涉及了深刻的数学理论和工程权衡。从拉格朗日在18世纪提出的优雅数学框架，到现代物理引擎中的实用算法，约束处理技术经历了漫长的演化。
对于强化学习应用而言，约束的准确建模直接影响到：
本章的学习目标：
在经典力学中，约束限制了系统的运动自由度。考虑一个包含 $n$ 个刚体的系统，每个刚体有6个自由度（3个平移+3个旋转），系统总共有 $N = 6n$ 个自由度。约束可以表示为关于系统状态的方程：
完整约束（Holonomic Constraints）： \(\mathbf{C}(\mathbf{q}, t) = \mathbf{0}\)
其中 $\mathbf{q} \in \mathbb{R}^N$ 是广义坐标，$\mathbf{C}: \mathbb{R}^N \times \mathbb{R} \rightarrow \mathbb{R}^m$ 是 $m$ 个约束方程。
非完整约束（Nonholonomic Constraints）： \(\mathbf{A}(\mathbf{q}, t)\dot{\mathbf{q}} + \mathbf{b}(\mathbf{q}, t) = \mathbf{0}\)
这类约束涉及速度，无法积分为位置约束。
对于完整约束，通过对时间求导，我们得到速度级约束： \(\dot{\mathbf{C}} = \mathbf{J}\dot{\mathbf{q}} + \frac{\partial \mathbf{C}}{\partial t} = \mathbf{0}\)
其中 $\mathbf{J} = \frac{\partial \mathbf{C}}{\partial \mathbf{q}}$ 是约束雅可比矩阵。
进一步求导得到加速度级约束： \(\ddot{\mathbf{C}} = \mathbf{J}\ddot{\mathbf{q}} + \dot{\mathbf{J}}\dot{\mathbf{q}} + \frac{\partial^2 \mathbf{C}}{\partial t^2} = \mathbf{0}\)
达朗贝尔原理是约束动力学的基石。它指出，在任何时刻，作用在系统上的主动力与惯性力的合力在任何与约束相容的虚位移上做的虚功为零：
\[(\mathbf{F} - \mathbf{M}\ddot{\mathbf{q}})^T \delta\mathbf{q} = 0\]其中 $\mathbf{F}$ 是主动力，$\mathbf{M}$ 是质量矩阵，$\delta\mathbf{q}$ 是虚位移。
关键洞察：约束力 $\mathbf{F}_c$ 垂直于约束流形，因此可以表示为： \(\mathbf{F}_c = \mathbf{J}^T \boldsymbol{\lambda}\)
其中 $\boldsymbol{\lambda} \in \mathbb{R}^m$ 就是著名的拉格朗日乘子。
结合牛顿第二定律和约束条件，我们得到约束系统的运动方程：
\[\begin{cases} \mathbf{M}\ddot{\mathbf{q}} = \mathbf{F} + \mathbf{J}^T \boldsymbol{\lambda} \\ \mathbf{J}\ddot{\mathbf{q}} = -\dot{\mathbf{J}}\dot{\mathbf{q}} \end{cases}\]这就是拉格朗日第一类方程（Lagrange’s equations of the first kind）。将其写成矩阵形式：
\[\begin{bmatrix} \mathbf{M} & -\mathbf{J}^T \\ \mathbf{J} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{q}} \\ \boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} \mathbf{F} \\ -\dot{\mathbf{J}}\dot{\mathbf{q}} \end{bmatrix}\]这是一个混合的微分-代数方程组（DAE），其中：
几何解释： 约束将系统限制在配置空间的一个子流形上。拉格朗日乘子的作用是提供垂直于这个流形的力，使系统保持在流形上运动。
约束流形
╱╲
╱ ╲
╱ ╲ ← 无约束轨迹
╱ ╲
╱ •────→ F (主动力)
╱ ↓
╱ λJ^T (约束力)
╱ ↓
╱ •────→ 约束轨迹
拉格朗日乘子 $\boldsymbol{\lambda}$ 具有明确的物理意义：它表示维持约束所需的力的大小。考虑一个简单的例子：
例：摆的约束力 一个质量为 $m$ 的质点通过长度为 $l$ 的刚性杆连接到原点。约束方程为： \(C = x^2 + y^2 - l^2 = 0\)
约束雅可比： \(\mathbf{J} = [2x, 2y]\)
约束力： \(\mathbf{F}_c = \lambda \mathbf{J}^T = 2\lambda [x, y]^T\)
这正是指向原点的向心力，$\lambda$ 的值决定了杆的张力。
在实际实现中，我们需要求解线性系统来获得加速度和拉格朗日乘子。常用的方法包括：
方法1：直接求解（适用于小规模系统） 直接求解混合系统： \(\boldsymbol{\lambda} = (\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^T)^{-1}(\mathbf{J}\mathbf{M}^{-1}\mathbf{F} + \dot{\mathbf{J}}\dot{\mathbf{q}})\) \(\ddot{\mathbf{q}} = \mathbf{M}^{-1}(\mathbf{F} + \mathbf{J}^T\boldsymbol{\lambda})\)
这里的 $\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^T$ 称为Delassus算子，在接触力学中起重要作用。
方法2：Schur补法（更稳定） 利用Schur补，我们可以先求解约束力：
稀疏性利用： 对于关节体系统，约束雅可比矩阵通常是稀疏的。例如，铰链约束只影响相连的两个刚体：
\[\mathbf{J} = \begin{bmatrix} \mathbf{J}_1 & -\mathbf{J}_2 & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{J}_3 & -\mathbf{J}_4 & \cdots & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots \end{bmatrix}\]利用这种稀疏结构可以大大提高计算效率。
即使我们在加速度级精确满足约束，数值积分仍会导致约束违背。这是因为：
考虑简单的显式欧拉积分： \(\mathbf{v}_{n+1} = \mathbf{v}_n + h\ddot{\mathbf{q}}_n\) \(\mathbf{q}_{n+1} = \mathbf{q}_n + h\mathbf{v}_{n+1}\)
即使 $\ddot{\mathbf{q}}n$ 满足加速度约束，$\mathbf{q}{n+1}$ 也可能违背位置约束。这种违背会随时间累积，最终导致系统完全偏离约束流形。
漂移的数学分析： 设约束误差为 $\boldsymbol{\epsilon} = \mathbf{C}(\mathbf{q})$，其动力学为： \(\ddot{\boldsymbol{\epsilon}} = \mathbf{J}\ddot{\mathbf{q}} + \dot{\mathbf{J}}\dot{\mathbf{q}}\)
如果我们只强制 $\ddot{\boldsymbol{\epsilon}} = 0$，则误差满足： \(\boldsymbol{\epsilon}(t) = \boldsymbol{\epsilon}_0 + t\dot{\boldsymbol{\epsilon}}_0\)
这是线性增长！误差会无限制地增大。
Baumgarte在1972年提出了一个优雅的解决方案：不是强制 $\ddot{\mathbf{C}} = 0$，而是强制：
\[\ddot{\mathbf{C}} + 2\alpha\dot{\mathbf{C}} + \beta^2\mathbf{C} = \mathbf{0}\]这将约束误差的动力学改变为一个稳定的二阶系统。选择合适的参数 $\alpha$ 和 $\beta$，误差将指数衰减而不是线性增长。
控制论视角： Baumgarte方法本质上是一个PD控制器：
修改后的约束方程变为： \(\mathbf{J}\ddot{\mathbf{q}} = -\dot{\mathbf{J}}\dot{\mathbf{q}} - 2\alpha\mathbf{J}\dot{\mathbf{q}} - \beta^2\mathbf{C}\)
参数选择对系统性能至关重要。考虑误差动力学的特征方程： \(s^2 + 2\alpha s + \beta^2 = 0\)
根的位置决定了系统行为：
临界阻尼（$\alpha = \beta$）： \(s_{1,2} = -\alpha\) 误差最快收敛，无振荡。推荐用于大多数应用。
过阻尼（$\alpha > \beta$）： \(s_{1,2} = -\alpha \pm \sqrt{\alpha^2 - \beta^2}\) 收敛较慢，但更稳定。
欠阻尼（$\alpha < \beta$）： \(s_{1,2} = -\alpha \pm i\sqrt{\beta^2 - \alpha^2}\) 有振荡，可能导致视觉artifacts。
实用参数选择指南：
频率域分析： Baumgarte参数影响系统的频率响应。高频约束（如碰撞）需要更大的 $\beta$ 值，而低频约束（如关节限制）可以使用较小的值。
从现代控制理论角度，Baumgarte方法可以推广为更一般的误差反馈形式：
\[\ddot{\mathbf{C}}_{\text{desired}} = -\mathbf{K}_p\mathbf{C} - \mathbf{K}_d\dot{\mathbf{C}} - \mathbf{K}_i\int\mathbf{C}dt\]其中：
自适应Baumgarte： 根据约束违背的大小动态调整参数： \(\alpha = \alpha_0 \cdot \min(1, \frac{\|\mathbf{C}\|}{\epsilon_{\text{tol}}})\)
这在约束违背较大时提供更强的纠正，而在接近满足时减少干扰。
多约束系统的耦合问题： 当系统有多个约束时，不同约束的Baumgarte参数可能相互影响。解决方案：
投影方法是处理约束漂移的另一种重要技术。基本思想是：先进行无约束积分，然后将结果投影回约束流形。
速度投影算法：
投影步骤等价于求解优化问题： \(\min_{\Delta\mathbf{v}} \frac{1}{2}\Delta\mathbf{v}^T\mathbf{M}\Delta\mathbf{v} \quad \text{s.t.} \quad \mathbf{J}(\mathbf{v}^* + \Delta\mathbf{v}) = \mathbf{0}\)
使用拉格朗日乘子法，解为： \(\Delta\mathbf{v} = -\mathbf{M}^{-1}\mathbf{J}^T(\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^T)^{-1}\mathbf{J}\mathbf{v}^*\)
几何解释： 速度投影是将速度向量正交投影到约束切空间（零空间）上：
v*（无约束速度）
↓
• ← 投影
╱│
╱ │Δv
╱ │
╱ ↓
╱ v_proj（投影后速度）
╱________________
约束切空间
位置投影直接纠正配置空间中的约束违背：
非线性投影（精确但昂贵）： 求解非线性方程组： \(\mathbf{C}(\mathbf{q} + \Delta\mathbf{q}) = \mathbf{0}\)
使用牛顿-拉夫逊迭代： \(\Delta\mathbf{q}_k = -\mathbf{J}^+\mathbf{C}(\mathbf{q}_k)\)
其中 $\mathbf{J}^+ = \mathbf{J}^T(\mathbf{J}\mathbf{J}^T)^{-1}$ 是雅可比的伪逆。
线性化投影（快速近似）： 假设约束在当前点附近是线性的： \(\mathbf{C}(\mathbf{q} + \Delta\mathbf{q}) \approx \mathbf{C}(\mathbf{q}) + \mathbf{J}\Delta\mathbf{q} = \mathbf{0}\)
解为： \(\Delta\mathbf{q} = -\mathbf{J}^+\mathbf{C}(\mathbf{q})\)
质量加权投影： 考虑系统的惯性特性： \(\Delta\mathbf{q} = -\mathbf{M}^{-1}\mathbf{J}^T(\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^T)^{-1}\mathbf{C}\)
这保持了系统的动量特性。
后稳定化（Post-stabilization）在每个时间步结束后应用投影：
完整的积分流程：
1. 预测步：
v* = v_n + h*M^(-1)*F
q* = q_n + h*v*
2. 速度投影：
λ_v = solve(J*M^(-1)*J^T*λ = J*v*)
v** = v* - M^(-1)*J^T*λ_v
3. 位置投影：
λ_q = solve(J*M^(-1)*J^T*λ = C(q*)/h)
q_{n+1} = q* - h*M^(-1)*J^T*λ_q
v_{n+1} = v** - M^(-1)*J^T*λ_q
分离式vs耦合式投影：
耦合投影需要求解更大的系统： \(\begin{bmatrix} \mathbf{M} & \mathbf{J}^T \\ \mathbf{J} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \Delta\mathbf{v} \\ \boldsymbol{\mu} \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ -\mathbf{J}\mathbf{v}^* - \frac{1}{h}\mathbf{C} \end{bmatrix}\)
从微分几何角度，约束定义了配置空间中的一个子流形 $\mathcal{M}$。投影方法的本质是找到流形上最近的点。
黎曼度量下的投影： 配置空间的自然度量由质量矩阵 $\mathbf{M}$ 诱导： \(d(\mathbf{q}_1, \mathbf{q}_2)^2 = (\mathbf{q}_1 - \mathbf{q}_2)^T\mathbf{M}(\mathbf{q}_1 - \mathbf{q}_2)\)
投影是在这个度量下的最近点问题： \(\mathbf{q}_{\text{proj}} = \arg\min_{\mathbf{q} \in \mathcal{M}} d(\mathbf{q}, \mathbf{q}^*)^2\)
切空间和法空间： 在约束流形上的点 $\mathbf{q}$，切空间 $T_{\mathbf{q}}\mathcal{M}$ 由： \(T_{\mathbf{q}}\mathcal{M} = \{\mathbf{v} : \mathbf{J}(\mathbf{q})\mathbf{v} = \mathbf{0}\}\)
法空间 $N_{\mathbf{q}}\mathcal{M}$ 由约束梯度张成： \(N_{\mathbf{q}}\mathcal{M} = \text{span}(\mathbf{J}^T)\)
投影就是沿法方向移动到流形上。
投影的迭代改进： 对于大的约束违背，单次投影可能不够准确。可以使用迭代投影：
for i = 1 to max_iterations:
C = compute_constraint(q)
if ||C|| < tolerance:
break
J = compute_jacobian(q)
Δq = -J^+ * C
q = q + α * Δq // α是步长因子
步长因子 $\alpha \in (0,1]$ 可以提高收敛的鲁棒性。
约束违背在实际仿真中是不可避免的，主要来源包括：
数值来源：
系统来源：
定量分析： 设真实解为 $\mathbf{q}^$，数值解为 $\mathbf{q}$，误差 $\mathbf{e} = \mathbf{q} - \mathbf{q}^$。 约束违背的泰勒展开： \(\mathbf{C}(\mathbf{q}) = \mathbf{C}(\mathbf{q}^*) + \mathbf{J}\mathbf{e} + \frac{1}{2}\mathbf{e}^T\mathbf{H}\mathbf{e} + O(\|\mathbf{e}\|^3)\)
其中 $\mathbf{H}$ 是约束的Hessian张量。
硬约束（Hard Constraints）： 必须精确满足的约束，如关节连接。
实现方式：
while ||C(q)|| > ε_hard:
λ = solve_constraint_force()
q = update_with_constraint(q, λ)
软约束（Soft Constraints）： 允许小幅违背的约束，通过罚函数实现。
罚函数方法： \(\mathbf{F}_{\text{penalty}} = -k_p\mathbf{C} - k_d\dot{\mathbf{C}}\)
其中 $k_p$ 是刚度系数，$k_d$ 是阻尼系数。
混合方法： 根据约束类型和违背程度动态切换：
if ||C|| < threshold_soft:
use_soft_constraint()
elif ||C|| < threshold_hard:
use_baumgarte_stabilization()
else:
use_projection_method()
当系统过约束或约束冲突时，需要优先级机制：
优先级分层：
分层求解算法：
for priority in [0, 1, 2, 3]:
constraints = get_constraints_at_priority(priority)
null_space = compute_null_space(higher_priority_constraints)
λ = solve_in_null_space(constraints, null_space)
apply_constraint_forces(λ)
零空间投影： 确保低优先级约束不影响高优先级： \(\mathbf{J}_{\text{low}} = \mathbf{J}_{\text{low}}(\mathbf{I} - \mathbf{J}_{\text{high}}^+\mathbf{J}_{\text{high}})\)
自适应容差： 根据系统状态动态调整容差： \(\epsilon_{\text{tol}} = \epsilon_0 \cdot (1 + \alpha\|\mathbf{v}\| + \beta\|\mathbf{F}\|)\)
约束正则化： 添加正则化项避免数值问题： \((\mathbf{J}\mathbf{M}^{-1}\mathbf{J}^T + \epsilon\mathbf{I})\boldsymbol{\lambda} = \mathbf{b}\)
其中 $\epsilon$ 是正则化参数（Tikhonov正则化）。
约束缩放： 归一化不同量纲的约束： \(\tilde{\mathbf{C}}_i = \frac{\mathbf{C}_i}{w_i}, \quad w_i = \max(\|\mathbf{J}_i\|, \epsilon_{\text{scale}})\)
紧急恢复机制： 当约束违背超过阈值时的处理：
铰链关节（Revolute Joint）： 约束两个刚体只能绕固定轴相对旋转。 约束方程（5个约束）： \(\mathbf{C} = \begin{bmatrix} \mathbf{r}_1 + \mathbf{R}_1\mathbf{s}_1 - \mathbf{r}_2 - \mathbf{R}_2\mathbf{s}_2 \\ \mathbf{n}_1^T\mathbf{n}_2 \\ \mathbf{n}_1^T\mathbf{b}_2 \end{bmatrix} = \mathbf{0}\)
其中 $\mathbf{s}_i$ 是局部锚点，$\mathbf{n}_i$ 是旋转轴，$\mathbf{b}_2$ 是垂直于 $\mathbf{n}_2$ 的向量。
球关节（Spherical Joint）： 允许绕点的任意旋转（3个约束）： \(\mathbf{C} = \mathbf{r}_1 + \mathbf{R}_1\mathbf{s}_1 - \mathbf{r}_2 - \mathbf{R}_2\mathbf{s}_2 = \mathbf{0}\)
滑动关节（Prismatic Joint）： 只允许沿固定轴平移（5个约束）： \(\mathbf{C} = \begin{bmatrix} (\mathbf{r}_1 - \mathbf{r}_2) \times \mathbf{n} \\ \mathbf{R}_1^T\mathbf{n}_1 - \mathbf{R}_2^T\mathbf{n}_2 \end{bmatrix} = \mathbf{0}\)
关节限制（Joint Limits）： 使用不等式约束： \(\theta_{\min} \leq \theta \leq \theta_{\max}\)
转换为互补性条件： \(\lambda_{\min} \geq 0, \quad \theta - \theta_{\min} \geq 0, \quad \lambda_{\min}(\theta - \theta_{\min}) = 0\)
在强化学习中，约束参数的选择直接影响策略学习的效果：
约束刚度的影响：
优化策略： \(k_{\text{optimal}} = \arg\min_{k} \mathbb{E}_{\pi}[L_{\text{task}} + \alpha L_{\text{constraint}}(k)]\)
约束噪声的作用： 添加适度的约束噪声可以提高策略的鲁棒性： \(\mathbf{C}_{\text{noisy}} = \mathbf{C} + \mathcal{N}(0, \sigma^2)\)
训练时逐步减小噪声（课程学习）： \(\sigma(t) = \sigma_0 \cdot \exp(-t/\tau)\)
约束松弛调度： 早期训练阶段放松约束，促进探索：
if episode < warmup_episodes:
constraint_strength = 0.1
elif episode < curriculum_episodes:
constraint_strength = 0.1 + 0.9 * (episode - warmup_episodes) / curriculum_episodes
else:
constraint_strength = 1.0
动态约束松弛： 根据策略性能动态调整约束强度： \(k(t) = k_{\min} + (k_{\max} - k_{\min}) \cdot \sigma(\text{performance})\)
其中 $\sigma$ 是基于性能的sigmoid函数。
约束感知的探索： 在action space中考虑约束可行域： \(a_{\text{feasible}} = \text{proj}_{\mathcal{F}}(a_{\text{raw}} + \epsilon)\)
其中 $\mathcal{F}$ 是约束定义的可行域，$\epsilon$ 是探索噪声。
分层控制架构：
约瑟夫-路易·拉格朗日（Joseph-Louis Lagrange, 1736-1813）在1788年出版的《分析力学》（Mécanique Analytique）彻底改变了力学的研究方法。他自豪地宣称：”在这本书中没有一张图”，将力学完全建立在分析（微积分）的基础上。
革命性贡献：
拉格朗日方法的优雅在于其普适性。无论是天体运动还是机械系统，都可以用同一套数学框架描述。这种统一性对现代物理仿真至关重要。
拉格朗日的工作奠定了变分力学的基础，其影响延续至今：
最小作用量原理： \(\delta S = \delta \int_{t_1}^{t_2} L(\mathbf{q}, \dot{\mathbf{q}}, t) dt = 0\)
这个原理不仅适用于经典力学，还贯穿了整个物理学：
现代计算方法的基础：
非完整约束是速度相关的约束，无法积分为位置约束。典型例子包括：
数学表示： \(\mathbf{A}(\mathbf{q})\dot{\mathbf{q}} = \mathbf{0}\)
其中 $\mathbf{A}$ 不是某个函数的梯度。
可积性条件（Frobenius定理）： 约束可积（完整）当且仅当： \(\mathbf{A}_i \cdot (\nabla \times \mathbf{A}_j) = 0, \quad \forall i,j\)
处理非完整系统有两种主要方法：
Nonholonomic方法（物理正确）： 直接应用达朗贝尔原理： \((\mathbf{F} - \mathbf{M}\ddot{\mathbf{q}})^T\delta\mathbf{q} = 0\) 其中 $\delta\mathbf{q}$ 满足 $\mathbf{A}\delta\mathbf{q} = 0$。
Vakonomic方法（变分一致）： 先变分后约束： \(\delta \int L dt = 0, \quad \text{subject to} \quad \mathbf{A}\dot{\mathbf{q}} = 0\)
两种方法通常给出不同的运动方程！Nonholonomic是物理正确的。
轮式机器人： 无滑动约束： \(\dot{x}\sin\theta - \dot{y}\cos\theta = 0\)
这导致机器人不能侧向移动，必须通过操纵实现平行泊车。
空间机器人： 角动量守恒（无外力矩）： \(\mathbf{L} = \sum_i \mathbf{I}_i\boldsymbol{\omega}_i = \text{const}\)
通过内部运动改变姿态（猫的翻正反射）。
软体机器人： 不可伸展约束： \(\|\frac{\partial \mathbf{r}}{\partial s}\| = 1\)
其中 $s$ 是弧长参数。
约束动力学是物理仿真的核心，本章系统介绍了处理约束的主要方法：
关键概念：
核心公式汇总：
约束系统运动方程： \(\begin{bmatrix}\mathbf{M} & -\mathbf{J}^T \\ \mathbf{J} & \mathbf{0}\end{bmatrix}\begin{bmatrix}\ddot{\mathbf{q}} \\ \boldsymbol{\lambda}\end{bmatrix} = \begin{bmatrix}\mathbf{F} \\ -\dot{\mathbf{J}}\dot{\mathbf{q}}\end{bmatrix}\)
Delassus算子：$\mathbf{D} = \mathbf{J}\mathbf{M}^{-1}\mathbf{J}^T$
投影算子：$\mathbf{P} = \mathbf{I} - \mathbf{J}^+\mathbf{J}$
方法选择指南：
题1：推导双摆系统的约束方程和雅可比矩阵。
题2：证明Baumgarte方法在临界阻尼条件下（$\alpha = \beta$）误差按 $e^{-\alpha t}$ 衰减。
题3：说明为什么投影方法保持系统的对称性。
题4：分析罚函数方法中刚度系数 $k$ 对系统频率的影响。
题5：设计一个自适应Baumgarte参数调整算法，根据约束违背的历史自动调整 $\alpha$ 和 $\beta$。
题6：对于一个包含 $n$ 个刚体和 $m$ 个约束的系统，分析不同约束处理方法的计算复杂度。
题7：证明非完整约束系统不满足能量守恒，即使没有外力。
题8：设计一个混合约束求解器，对不同类型的约束使用不同的处理方法。
陷阱1：Baumgarte参数设置过大
陷阱2：约束雅可比奇异
陷阱3：投影顺序依赖
陷阱4：过约束系统
陷阱5：约束初始化不一致
陷阱6：忽略约束的非线性
陷阱7：密集矩阵运算
陷阱8：过度迭代