3. 多体动力学方程建立

多体系统动力学的核心任务，是把由多个刚体（或柔性体）通过铰链、滑移副、球铰等连接而成的机械系统，转化为一组可求解的数学方程。不同建模路线得到的方程形式不同，但描述的是同一个物理系统。理解它们的推导逻辑、适用场景与相互转换关系，是后续参数化优化的前提——因为优化器最终要反复求解这些方程，方程形式直接决定计算效率与灵敏度信息的获取方式。

一、广义坐标与约束：方程建立的共同起点

无论采用哪种方法，第一步都是选择描述系统位形的坐标。对含 n_b 个刚体的平面系统，若每个刚体用质心坐标 (x, y) 和姿态角 θ 描述，则笛卡尔坐标数为 3n_b。这些坐标并非独立，铰链会施加限制，例如转动副要求两刚体连接点重合，即两个标量约束方程。设系统共有 m 个独立约束，则系统自由度为 DOF = 3n_b − m。

约束一般写成位置约束向量形式 Φ(q, t) = 0。对时间求导得到速度约束 Φ_q q̇ = −Φ_t，再求导得到加速度约束 Φ_q q̈ = γ，其中 γ 汇集了速度二次项与时间显式项。这个 Φ_q 就是约束雅可比矩阵，它同时是后续微分代数方程和灵敏度分析的关键矩阵。需要强调：本文统一约定约束以 Φ = 0 表示，雅可比按 Φ_q = ∂Φ/∂q 定义，全文与图示保持一致。

二、拉格朗日方程：能量视角的建模路线

拉格朗日方法从标量能量出发，避开铰链内部约束力的求解。对无约束的广义坐标 q，方程写作：

d/dt(∂T/∂q̇)ᵀ − (∂T/∂q)ᵀ + Φqᵀ λ = Q
Φ(q, t) = 0

其中 T 为动能，Q 为广义力（含弹簧、阻尼、驱动），λ 为拉格朗日乘子，物理意义是约束反力。乘子项 Φ_q^Tλ 正是把约束力引入方程的方式。若改用独立广义坐标（如各铰链转角），可消去乘子，得到纯微分方程，维数更低，但推导独立坐标变换对复杂拓扑系统较为繁琐。

拉格朗日路线的优势在于：动能与势能都是标量，便于用符号推导，也便于把柔性体、驱动元件统一纳入能量表达。在参数化优化中，若设计变量是杆长、质心位置或刚度，能量形式对参数的导数往往更易解析求取。

三、牛顿-欧拉方法：逐体受力与递推

牛顿-欧拉方法对每个刚体分别写出牛顿方程与欧拉方程：

m_i a_i = F_i^ext + F_i^joint
J_i α_i + ω_i × (J_i ω_i) = M_i^ext + M_i^joint

它显式包含铰链约束力与力矩，物理图像直观，适合处理接触、摩擦、碰撞等需要真实反力的场景。工程实现中常配合递推算法：先由基座向外做运动学递推，算出各体速度、加速度；再由外向内做动力学递推，传递力与力矩。这种 O(n) 复杂度算法在链式或树状机器人上效率很高。

代价是方程数量多、约束力作为未知量一并求解，且对闭环系统需要额外处理。因此实际软件常把牛顿-欧拉用于开链高效计算，把拉格朗日或 DAE 形式用于含闭环的通用建模。

四、约束方程与微分代数方程形式

当保留笛卡尔坐标并显式保留约束时，运动方程与约束方程联立，形成微分代数方程（DAE）。最常用的指标-3 形式为：

[ M   Φqᵀ ] [ q̈ ]   [ Q − Φqᵀ... 的显式项 ]
[ Φq   0  ] [ λ  ] = [ γ            ]

这是关于 q̈ 和 λ 的线性方程组，每步积分需求解一个鞍点结构的线性系统。它的优点是约束漂移可通过 Baumgarte 稳定化或投影修正控制；缺点是代数约束使系统指标升高，直接积分易出现数值不稳定。若把约束二次求导后只积分 q̈，位置约束会随时间漂移，因此工程上常配合约束稳定项：

Φq q̈ = γ − 2α(Φq q̇ + βΦ) − β²Φ

其中 α、β 为稳定参数，典型取值使约束误差按指数衰减。另一种思路是把 DAE 降为常微分方程（ODE），即用独立坐标消去乘子，但需要处理坐标分区与拓扑变化。

五、三种路线的对比与选用

六、面向参数化优化的实现要点

在优化循环中，设计变量 p（如杆长、质量、刚度）会改变 M、Q、Φ_q 和 γ。若采用 DAE 形式，每步需求解线性系统，且需对 p 求导以获得灵敏度。常用做法是：先由方程对 p 求偏导，得到灵敏度 DAE，再与状态方程同步积分。这样可避免有限差分的重复求解，显著降低优化成本。

综上，拉格朗日方程提供能量视角的紧凑表达，牛顿-欧拉方法提供逐体受力的直观递推，约束方程与 DAE 形式则把两者统一到可数值求解的框架中。理解它们的等价性与差异，才能在参数化优化中根据系统拓扑、计算规模与灵敏度需求，选择最合适的方程形式。