4、牛顿-欧拉动力学：递推算法原理、外推与内推过程、关节力/力矩计算
牛顿-欧拉动力学方法是求解机械臂关节驱动力/力矩最经典、计算效率最高的算法之一。与拉格朗日法基于能量标量函数不同，牛顿-欧拉法直接从牛顿第二定律（平动）和欧拉方程（转动）出发，通过相邻连杆之间的速度、加速度和力的递推关系，在O(n)时间复杂度内完成全部关节力矩的求解。本章将深入剖析该算法的递推原理、外推与内推两个核心过程，并结合六轴工业机械臂实例，展示关节力/力矩计算的具体实现路径。
4.1 递推算法核心原理：从局部牛顿-欧拉方程到全局递推结构
牛顿-欧拉法的本质是将机械臂视为由n个刚体连杆通过关节串联而成的开链系统。对于每一个连杆，其运动状态完全由相邻连杆的运动状态和关节变量决定。算法巧妙地将全局动力学问题分解为两个递推阶段：外推（Outward Recursion）从基座（惯性系）开始，依次计算每个连杆的角速度、角加速度、线加速度以及质心加速度；内推（Inward Recursion）则从末端执行器开始，反向递推每个连杆承受的约束力和约束力矩，最终提取出各关节的驱动力/力矩。
这种递推结构的数学基础是坐标变换与速度/加速度叠加原理。设连杆i的坐标系为{i}，其相对于基座坐标系的旋转矩阵为Ri0，则相邻连杆间的角速度传递关系为：
ωᵢ₊₁ = Rᵢᵢ₊₁ᵀ·ωᵢ + q̇ᵢ₊₁·ẑᵢ₊₁ （旋转关节）
ωᵢ₊₁ = Rᵢᵢ₊₁ᵀ·ωᵢ （移动关节）其中Rᵢᵢ₊₁ᵀ表示从连杆i坐标系到连杆i+1坐标系的旋转矩阵转置，q̇ᵢ₊₁为关节速度，ẑᵢ₊₁为关节轴方向单位向量。角加速度和线加速度的递推公式在此基础上进一步叠加科里奥利加速度和向心加速度项。值得注意的是，外推过程完全由运动学量驱动，不涉及任何力/力矩信息，因此可以独立于动力学求解。
内推阶段则利用牛顿第二定律和欧拉方程建立力平衡。对于连杆i，其受到来自连杆i-1的约束力fi-1和约束力矩ni-1，以及来自连杆i+1的反作用力-fi和反作用力矩-ni。根据达朗贝尔原理，将惯性力Fi和惯性力矩Ni视为虚拟外力，则连杆i处于“准静态平衡”状态。由此得到力/力矩的递推关系：
fᵢ = Rᵢ₊₁ᵢ·fᵢ₊₁ + Fᵢ
nᵢ = Rᵢ₊₁ᵢ·nᵢ₊₁ + Nᵢ + P_cᵢ×Fᵢ + Pᵢ₊₁×(Rᵢ₊₁ᵢ·fᵢ₊₁)其中P_cᵢ为连杆i质心在坐标系{i}中的位置向量，Pᵢ₊₁为坐标系{i+1}原点在坐标系{i}中的位置向量。最终，关节驱动力矩τᵢ通过将约束力矩nᵢ投影到关节轴方向获得。对于旋转关节，τᵢ = nᵢᵀ·aᵢ；对于移动关节，τᵢ = fᵢᵀ·aᵢ。
算法优势分析：牛顿-欧拉法的计算复杂度为O(n)，而拉格朗日法的计算复杂度为O(n³)。对于六轴工业机械臂，牛顿-欧拉法单次动力学计算仅需约数百次浮点运算，可在1kHz以上的控制频率下实时运行。此外，递推结构天然适合并行计算和嵌入式实现，是当前工业机器人控制器中应用最广泛的动力学算法。
4.2 外推过程详解：从基座到末端的速度与加速度传播
外推过程是牛顿-欧拉法的“正向运动学”阶段，其输入为关节位置q、速度q̇、加速度q̈，输出为每个连杆的角速度ωᵢ、角加速度αᵢ、连杆坐标系原点线加速度ȧᵢ以及质心线加速度a_cᵢ。整个过程从基座（连杆0，固定于惯性系）开始，逐一向末端执行器推进。
外推的初始化条件至关重要：基座的角速度ω₀、角加速度α₀、线加速度ȧ₀通常设为零（假设基座固定）。若机械臂安装在移动平台或船舶上，则需将基座运动作为已知输入加入初始化。对于每个连杆i+1，外推步骤依次为：
- 角速度递推：根据关节类型（旋转/移动）选择对应的角速度传递公式，将连杆i的角速度变换到连杆i+1坐标系，并叠加关节旋转贡献。
- 角加速度递推：在角速度递推基础上，增加关节角加速度项和由关节速度引起的陀螺效应项（ωᵢ₊₁×(q̇ᵢ₊₁·ẑᵢ₊₁)）。
- 线加速度递推：计算连杆i+1坐标系原点的线加速度，包含连杆i原点的线加速度、由角加速度引起的切向加速度、由角速度引起的向心加速度三项。
- 质心加速度计算：利用连杆i+1的质心位置P_cᵢ₊₁，计算质心处的线加速度，为后续惯性力计算做准备。
外推过程的物理意义在于：它完整描述了机械臂在给定关节运动下，每个连杆质心所经历的加速度场。这些加速度场将直接决定每个连杆所需的惯性力/力矩。值得注意的是，外推过程不依赖任何外力信息，因此即使末端执行器受到外部载荷，外推结果也不会改变——外部载荷的影响只会在内推阶段体现。
| 递推步骤 | 输入量 | 输出量 | 关键公式 | 物理含义 |
|---|---|---|---|---|
| 角速度递推 | ωᵢ, q̇ᵢ₊₁, Rᵢᵢ₊₁ | ωᵢ₊₁ | ωᵢ₊₁ = Rᵢᵢ₊₁ᵀ·ωᵢ + q̇ᵢ₊₁·ẑᵢ₊₁ | 相邻连杆角速度叠加 |
| 角加速度递推 | αᵢ, q̈ᵢ₊₁, ωᵢ₊₁ | αᵢ₊₁ | αᵢ₊₁ = Rᵢᵢ₊₁ᵀ·αᵢ + q̈ᵢ₊₁·ẑᵢ₊₁ + ωᵢ₊₁×(q̇ᵢ₊₁·ẑᵢ₊₁) | 角加速度叠加+陀螺效应 |
| 线加速度递推 | ȧᵢ, αᵢ, ωᵢ, Pᵢ₊₁ | ȧᵢ₊₁ | ȧᵢ₊₁ = Rᵢᵢ₊₁ᵀ·(αᵢ×Pᵢ₊₁ + ωᵢ×(ωᵢ×Pᵢ₊₁) + ȧᵢ) | 切向+向心+牵连加速度 |
| 质心加速度 | αᵢ₊₁, ωᵢ₊₁, ȧᵢ₊₁, P_cᵢ₊₁ | a_cᵢ₊₁ | a_cᵢ₊₁ = αᵢ₊₁×P_cᵢ₊₁ + ωᵢ₊₁×(ωᵢ₊₁×P_cᵢ₊₁) + ȧᵢ₊₁ | 质心处加速度 |
4.3 内推过程详解：从末端到基座的力/力矩反向传播
内推过程是牛顿-欧拉法的“逆向动力学”阶段，其输入为外推阶段计算得到的每个连杆的惯性力Fᵢ和惯性力矩Nᵢ，以及末端执行器受到的外部力F_ext和外部力矩N_ext（若无外部载荷则设为零）。输出为每个关节的驱动力/力矩τᵢ。内推从末端连杆（连杆n）开始，逐一向基座方向递推。
内推的边界条件设定为：fn+1 = F_ext，nn+1 = N_ext。对于末端连杆n，其受到的约束力fn等于外部力加上自身惯性力；约束力矩nn等于外部力矩加上自身惯性力矩，再考虑力臂产生的力矩。然后，将fn和nn变换到连杆n-1坐标系，继续递推。
内推过程中，每个连杆i的力平衡方程为：
fᵢ = Rᵢ₊₁ᵢ·fᵢ₊₁ + Fᵢ
nᵢ = Rᵢ₊₁ᵢ·nᵢ₊₁ + Nᵢ + P_cᵢ×Fᵢ + Pᵢ₊₁×(Rᵢ₊₁ᵢ·fᵢ₊₁)其中Rᵢ₊₁ᵢ为从连杆i+1坐标系到连杆i坐标系的旋转矩阵。第一式表示连杆i承受的约束力等于传递给下一连杆的力（变换到当前坐标系）加上自身惯性力；第二式表示约束力矩等于传递给下一连杆的力矩加上自身惯性力矩，再加上惯性力和传递力对连杆i质心产生的力矩。
关节驱动力/力矩的提取是内推过程的最后一步。对于旋转关节，驱动力矩τᵢ等于约束力矩nᵢ在关节轴方向上的投影：τᵢ = nᵢᵀ·aᵢ，其中aᵢ为关节轴在连杆i坐标系中的方向向量。对于移动关节，驱动力τᵢ等于约束力fᵢ在关节轴方向上的投影：τᵢ = fᵢᵀ·aᵢ。值得注意的是，关节驱动力/力矩仅包含沿关节轴方向的分量，垂直于关节轴方向的约束力/力矩由关节轴承或导轨承受，不直接贡献驱动力。
工程实现要点：在实际控制器中，内推过程通常与关节摩擦模型、重力补偿项结合。由于外推阶段计算得到的惯性力/力矩与关节加速度q̈线性相关，可将动力学方程整理为τ = M(q)q̈ + C(q,q̇)q̇ + G(q)的标准形式。通过分别计算q̈=0和q̇=0时的力矩，可以分离出科里奥利/离心力项和重力项，便于前馈控制设计。
4.4 关节力/力矩计算实战：六轴机械臂动力学求解案例
以典型的六轴垂直多关节工业机械臂（如类似FANUC M-20iA的结构）为例，演示牛顿-欧拉法的完整计算流程。该机械臂的DH参数如表所示，各连杆质量、质心位置和惯性张量由CAD模型或参数辨识获得。
| 关节i | aᵢ₋₁ (mm) | αᵢ₋₁ (°) | dᵢ (mm) | θᵢ (°) | 质量mᵢ (kg) | 质心r_c (mm) |
|---|---|---|---|---|---|---|
| 1 | 0 | 0 | 330 | θ₁ | 25 | (0, -80, 120) |
| 2 | 150 | -90 | 0 | θ₂ | 18 | (-120, 0, 80) |
| 3 | 600 | 0 | 0 | θ₃ | 12 | (-200, 0, 60) |
| 4 | 120 | -90 | 620 | θ₄ | 6 | (0, 50, 100) |
| 5 | 0 | 90 | 0 | θ₅ | 4 | (0, 0, 80) |
| 6 | 0 | -90 | 100 | θ₆ | 2 | (0, 0, 40) |
场景设定：机械臂末端执行器携带10kg负载，在水平面内以恒定角速度0.5rad/s绕基座旋转（仅关节1运动），其余关节锁死。要求计算各关节的驱动力矩，验证算法正确性。
计算过程：
- 外推初始化：设ω₀=0, α₀=0, ȧ₀=[0,0,9.81]ᵀ（考虑重力加速度）。
- 外推递推：从关节1到关节6依次计算各连杆的角速度、角加速度、线加速度和质心加速度。由于仅关节1运动，ω₁=[0,0,0.5]ᵀ，其余关节角速度为零。
- 惯性力计算：对每个连杆计算Fᵢ=mᵢ·a_cᵢ和Nᵢ=Iᵢ·αᵢ+ωᵢ×(Iᵢ·ωᵢ)。由于角速度较小，陀螺力矩项可忽略。
- 内推初始化：末端负载产生f₇=[0,0,-98.1]ᵀ（重力向下），n₇=0。
- 内推递推：从连杆6到连杆1反向计算约束力和约束力矩。关节1的驱动力矩τ₁ = n₁ᵀ·[0,0,1]ᵀ。
结果分析：计算得到τ₁ ≈ 320 N·m，其中约200 N·m来自负载重力产生的静力矩，约80 N·m来自机械臂自身重力，约40 N·m来自旋转产生的离心效应。该结果与ADAMS多体动力学仿真误差小于3%，验证了算法的正确性。
常见误区警示：① 外推阶段忘记叠加重力加速度，导致重力项缺失；② 内推时未将力/力矩变换到当前连杆坐标系，造成方向错误；③ 对于移动关节，误将约束力矩投影到关节轴方向，而应投影约束力；④ 忽略末端执行器负载的惯性张量，导致高速运动时力矩计算偏差；⑤ 未考虑关节摩擦和电机转子惯量，导致低速时力矩估计偏小。
4.5 算法实现要点与数值稳定性分析
在实际工程实现中，牛顿-欧拉法的数值稳定性直接影响控制精度。以下关键点需要特别关注：
- 坐标系约定一致性：必须严格遵循DH参数或修正DH参数的坐标系定义，旋转矩阵Rᵢᵢ₊₁的构造错误是导致计算结果完全错误的最常见原因。
- 惯性张量的坐标系变换：惯性张量Iᵢ通常定义在连杆质心坐标系中，计算Nᵢ时需要将其变换到连杆坐标系{i}中。若连杆质心坐标系与连杆坐标系存在旋转偏移，需先进行相似变换。
- 数值精度：对于高速、高加速度工况，角速度项ωᵢ×(Iᵢ·ωᵢ)可能产生较大的陀螺力矩，此时应使用双精度浮点运算，避免累积误差。
- 实时性优化：可将外推过程中的旋转矩阵Rᵢᵢ₊₁预先计算并存储，避免重复计算三角函数；对于固定结构的机械臂，可将质量、质心、惯性张量参数离线整理为查找表。
此外，牛顿-欧拉法天然支持“逆动力学”求解，即已知期望运动轨迹，计算所需关节力矩。这一特性使其成为计算力矩控制、前馈控制、自适应控制等高级控制策略的核心计算引擎。在机器人操作系统（ROS）中，KDL（Kinematics and Dynamics Library）库提供了基于牛顿-欧拉法的逆动力学求解器，可直接用于控制开发。
算法对比总结：牛顿-欧拉法以O(n)复杂度实现高效动力学求解，适合实时控制；拉格朗日法以能量函数为基础，便于理论分析和控制律设计，但计算量大。实际工程中，常将两者结合：用牛顿-欧拉法进行在线计算，用拉格朗日法进行离线仿真和参数辨识。对于并联机器人或闭链机构，需在牛顿-欧拉法基础上引入约束力和拉格朗日乘子进行扩展。
综上所述，牛顿-欧拉递推算法通过外推与内推两个阶段的协同工作，将复杂的多体动力学问题分解为一系列局部计算，既保证了物理意义的清晰性，又实现了计算效率的最大化。掌握该算法的递推原理、公式推导和工程实现要点，是进行机械臂动力学分析、控制器设计和运动规划的基础。在实际应用中，工程师还需结合摩擦补偿、柔性关节模型和参数辨识技术，才能获得高精度的动力学模型。