质量矩阵与科氏力Mass Matrix & Coriolis Terms
机械臂的动力学方程 \(M(\mathbf q)\ddot{\mathbf q}+C(\mathbf q,\dot{\mathbf q})\dot{\mathbf q}+\mathbf g(\mathbf q)=\boldsymbol\tau\) 里,三项各管一件事:手臂此刻"有多沉"(质量矩阵),一个关节的运动怎样拽动别的关节(科氏力与离心力),重力在拉什么。这一节把它们拆开来看:质量矩阵是衡量动能的"尺子",科氏项是尺子随姿态变化的代价,再看工具端到底"感觉有多重",最后用它们做计算力矩控制和参数辨识。 In the manipulator equation \(M(\mathbf q)\ddot{\mathbf q}+C(\mathbf q,\dot{\mathbf q})\dot{\mathbf q}+\mathbf g(\mathbf q)=\boldsymbol\tau\) each term answers one question: how heavy the arm feels right now (the mass matrix), how motion in one joint tugs on the others (Coriolis and centrifugal terms), and what gravity is pulling on. This lesson takes them apart: the mass matrix is a ruler for kinetic energy, the Coriolis term is the price of that ruler changing with posture, the tool-tip then "feels" a very different mass in different directions, and finally computed-torque control and parameter identification put all three to work.
1为什么机器人需要它Why robots need this
8.1 节用拉格朗日方法、8.2 节用牛顿–欧拉递推,都得到了同一个方程:
Lesson 8.1 (Lagrange) and Lesson 8.2 (recursive Newton–Euler) arrive at the same equation:
递推算法回答"给定运动,要多大力矩",但它是一个黑盒。控制器要的往往是方程的结构:
- 重力补偿 / 手拖示教:只需要 \(\mathbf g(\mathbf q)\)。让关节输出 \(\boldsymbol\tau=\mathbf g(\mathbf q)\),手臂就"悬浮"在空中,你可以用手推着它走。
- 计算力矩控制(逆动力学控制):用 \(M\)、\(C\)、\(\mathbf g\) 把非线性的手臂"抵消"成一堆互不干扰的积分器,再用简单的 PD 去控制。
- 阻抗 / 操作空间控制:想让手爪像一个"质量–弹簧–阻尼",就得知道手爪端的等效质量 \(\Lambda(\mathbf q)\),它由 \(M\) 和雅可比算出。
- 仿真与全身控制:物理引擎每一步都在解 \(\ddot{\mathbf q}=M^{-1}(\boldsymbol\tau-C\dot{\mathbf q}-\mathbf g)\);足式机器人的全身控制把 \(M\) 放进二次规划的等式约束里(第四部分)。
- 参数辨识与 sim-to-real:方程对质量、质心、惯量是线性的,可以用最小二乘从数据里估计它们;仿真里惯量写错,策略迁移到真机就会"手感不对"。
这一节的目标是:看清每一项长什么样、为什么长这样、在控制里干什么。我们用同一个平面两连杆臂贯穿全篇:\(l_1=1\) m,\(l_2=0.8\) m,质量 \(m_1=2\) kg,\(m_2=1.5\) kg,两根连杆都是均匀杆(质心在中点,绕质心的转动惯量 \(I=ml^2/12\)),\(\theta_i\) 从水平方向量起(\(\mathbf q=(\theta_1,\theta_2)\)),重力沿 \(-y\)。
The recursive algorithm answers "given a motion, what torque?" but it is a black box. Controllers usually want the structure of the equation:
- Gravity compensation / hand-guiding: needs only \(\mathbf g(\mathbf q)\). Command \(\boldsymbol\tau=\mathbf g(\mathbf q)\) and the arm "floats" where you leave it, so you can lead it by hand.
- Computed-torque (inverse-dynamics) control: use \(M\), \(C\), \(\mathbf g\) to cancel the nonlinear arm into a set of independent double integrators, then control those with a plain PD law.
- Impedance / operational-space control: to make the gripper behave like a mass–spring–damper you must know the apparent mass \(\Lambda(\mathbf q)\) at the tool, computed from \(M\) and the Jacobian.
- Simulation and whole-body control: a physics engine solves \(\ddot{\mathbf q}=M^{-1}(\boldsymbol\tau-C\dot{\mathbf q}-\mathbf g)\) at every step; whole-body control for legged robots puts \(M\) in an equality constraint of a quadratic program (Part IV).
- Identification and sim-to-real: the equation is linear in the masses, centres of mass and inertias, so least squares can estimate them from data; wrong inertias in simulation make a policy "feel off" on the real robot.
The goal of this lesson is to see what each term looks like, why it has that form, and what it does in a controller. One planar two-link arm runs through the whole page: \(l_1=1\) m, \(l_2=0.8\) m, masses \(m_1=2\) kg, \(m_2=1.5\) kg, both links uniform rods (centre of mass at the midpoint, inertia about it \(I=ml^2/12\)), angles \(\theta_i\) measured from the horizontal (\(\mathbf q=(\theta_1,\theta_2)\)), gravity along \(-y\).
2直觉:\(M\) 是衡量动能的"尺子"Intuition: \(M\) is a ruler for kinetic energy
推一扇很重的门,推门轴附近的地方和推门的外沿,手感完全不同。手臂也是这样:关节 1 转起来"有多沉",取决于质量离关节 1 有多远。手臂伸直时,质量都甩在外面,转动惯量大;折起来时,质量缩在身边,转动惯量小。所以手臂的惯性不是一个数,而是随姿态变化的一个矩阵 \(M(\mathbf q)\)。
怎么把"有多沉"量化?看动能。给定关节速度 \(\dot{\mathbf q}\),手臂的总动能是一个关于 \(\dot{\mathbf q}\) 的二次型:
Push a heavy door near its hinge or at its outer edge and it feels completely different. An arm is the same: how "heavy" joint 1 feels depends on how far the mass sits from joint 1. Stretched out, the mass is flung outward and the rotational inertia is large; folded, the mass hugs the shoulder and it is small. So the arm's inertia is not a number but a matrix that changes with posture, \(M(\mathbf q)\).
How do we quantify "how heavy"? Look at the kinetic energy. For joint velocities \(\dot{\mathbf q}\), the arm's total kinetic energy is a quadratic form in \(\dot{\mathbf q}\):
1.7 节的二次型,1.11–1.12 节的特征值与椭圆在这里直接派上用场:\(M\) 对称正定,等能量面 \(\{\dot{\mathbf q}:\tfrac12\dot{\mathbf q}^\top M\dot{\mathbf q}=E\}\) 是关节速度平面里的一个椭圆。它的半轴是 \(\sqrt{2E/\lambda_i}\),沿 \(M\) 的特征向量:特征值大的方向很"沉",同样的能量只能转得很慢,椭圆在那里很窄;特征值小的方向很"轻",椭圆很长。下图左边拖动手臂的姿态,右边拖动琥珀色的速度手柄 \(\dot{\mathbf q}\),看动能怎么变。
The quadratic forms of Lesson 1.7 and the eigenvalues and ellipses of Lessons 1.11–1.12 pay off right here: \(M\) is symmetric positive definite, and the level set \(\{\dot{\mathbf q}:\tfrac12\dot{\mathbf q}^\top M\dot{\mathbf q}=E\}\) is an ellipse in the joint-velocity plane. Its semi-axes are \(\sqrt{2E/\lambda_i}\) along the eigenvectors of \(M\): a large eigenvalue is a "heavy" direction (the same energy buys little speed, so the ellipse is narrow), a small eigenvalue is a "light" direction (the ellipse is long). In the figure drag the posture on the left and the amber velocity handle \(\dot{\mathbf q}\) on the right, and watch the kinetic energy.
- 不加负载,把姿态调到让椭圆最圆(条件数 \(\kappa(M)=\lambda_{\max}/\lambda_{\min}<5.6\))。这时肘部大概弯了多少?上图曲线里哪个位置最圆?With no payload, pose the arm so the ellipse is as round as it gets (condition number \(\kappa(M)=\lambda_{\max}/\lambda_{\min}<5.6\)). Roughly how far is the elbow bent? Where on the lower plot is that?
- 把速度手柄拖到椭圆长轴的顶端(沿 \(M\) 最小特征值的方向,动能与 \(E\) 相差不到 10%)。同样的动能,关节速度的大小和在最"沉"的方向上相比差了多少倍?Drag the velocity handle to the tip of the ellipse's long axis (along the eigenvector of \(M\)'s smallest eigenvalue, kinetic energy within 10% of \(E\)). For the same energy, how many times larger is the joint speed than along the "heaviest" direction?
\(M(\mathbf q)\) 不是"质量",而是关节空间里衡量动能的度量(尺子):同一个关节速度,在不同姿态下的动能不同;同一个姿态下,沿不同方向的"沉"也不同。它只依赖手臂的形状:对这个平面臂,\(M\) 只与 \(\theta_2\) 有关,与 \(\theta_1\)(整条手臂怎么摆)无关。上图读数里"\(\theta_1\) 改变后 \(M\) 的变化"始终是 0。
\(M(\mathbf q)\) is not "the mass" but a metric (a ruler) for kinetic energy on joint space: the same joint velocity carries different energy in different postures, and in one posture some directions are heavier than others. It depends only on the arm's shape: for this planar arm \(M\) depends on \(\theta_2\) only, not on \(\theta_1\) (how the whole arm is swung). The readout's "change of \(M\) when \(\theta_1\) changes" is always 0.
3动手玩:离心力与科氏力是从哪来的Play: where centrifugal and Coriolis torques come from
"科氏力"听上去神秘,其实只是牛顿第二定律。让关节 1 以恒定角速度转动,连杆 2 的质心在画圆,圆周运动必须有指向圆心的加速度 \(\mathbf a=-\omega^2\mathbf r\),也就必须有一个指向圆心的合力 \(m\mathbf a\)。这个力要由关节 2 的销钉和电机提供,它对关节 2 的力矩就是电机必须出的力矩。如果电机不出,连杆就会"甩出去":这就是离心的感觉。
下图保持两个关节速度 \(\dot\theta_1,\dot\theta_2\) 恒定(\(\ddot{\mathbf q}=0\),先忽略重力),琥珀色箭头是维持这种运动,每根连杆质心所需的合力 \(m_i\mathbf a_i\)。右侧的条形图是关节力矩里与速度有关的部分 \(C(\mathbf q,\dot{\mathbf q})\dot{\mathbf q}\),按来源拆开:
"Coriolis force" sounds mysterious, but it is just Newton's second law. Spin joint 1 at a constant rate: link 2's centre of mass travels on a circle, and circular motion needs an acceleration toward the centre, \(\mathbf a=-\omega^2\mathbf r\), hence a net force \(m\mathbf a\) toward the centre. Joint 2's pin and motor must supply that force, and its moment about joint 2 is the torque the motor has to deliver. If the motor did not, the link would "fly outward": that is the feeling of centrifugal force.
The figure keeps both joint rates \(\dot\theta_1,\dot\theta_2\) constant (\(\ddot{\mathbf q}=0\), gravity ignored for now). The amber arrows are the net force \(m_i\mathbf a_i\) each link's centre of mass needs to keep this motion going. The bars on the right are the velocity-dependent part of the joint torques, \(C(\mathbf q,\dot{\mathbf q})\dot{\mathbf q}\), split by origin:
- 离心项(含 \(\dot q_i^2\)):某个关节自己的速度平方,例如 \(\tau_2\ni h\dot\theta_1^2\)。
- Centrifugal terms (containing \(\dot q_i^2\)): a joint's own rate squared, e.g. \(\tau_2\ni h\dot\theta_1^2\).
- 科氏项(含 \(\dot q_i\dot q_j,\ i\ne j\)):两个关节速度的乘积,例如 \(\tau_1\ni-2h\dot\theta_1\dot\theta_2\)。
- Coriolis terms (containing \(\dot q_i\dot q_j,\ i\ne j\)): a product of two joint rates, e.g. \(\tau_1\ni-2h\dot\theta_1\dot\theta_2\).
其中 \(h=m_2l_1l_{c2}\sin\theta_2\)。按下"播放"让手臂真的转起来(\(\theta_2\) 会跟着变),或者拖动滑块与姿态手柄。
Here \(h=m_2l_1l_{c2}\sin\theta_2\). Press "Play" to actually spin the arm (\(\theta_2\) will change with it), or drag the sliders and posture handles.
- 只转关节 1(\(\dot\theta_2=0\),\(|\dot\theta_1|>3\) rad/s)。关节 2 明明没有加速,为什么要出 \(|\tau_2|>3\) N·m?把 \(\theta_2\) 拉到 0 或 π(手臂伸直 / 折回),这个力矩怎么变?Spin only joint 1 (\(\dot\theta_2=0\), \(|\dot\theta_1|>3\) rad/s). Joint 2 is not accelerating, so why does it need \(|\tau_2|>3\) N·m? Pull \(\theta_2\) to 0 or π (arm straight or folded): what happens to the torque?
- 让两个关节都在转(\(|\dot\theta_1|>1\),\(|\dot\theta_2|>1\)),而关节 1 的速度相关力矩几乎为零(\(|\tau_1|<0.3\) N·m,\(|\sin\theta_2|>0.35\))。提示:\(2\dot\theta_1\dot\theta_2+\dot\theta_2^2=0\)。Make both joints spin (\(|\dot\theta_1|>1\), \(|\dot\theta_2|>1\)) while joint 1's velocity-dependent torque is nearly zero (\(|\tau_1|<0.3\) N·m, \(|\sin\theta_2|>0.35\)). Hint: \(2\dot\theta_1\dot\theta_2+\dot\theta_2^2=0\).
- 在右上角切换到"另一种 \(C\)":\(C\dot{\mathbf q}\) 的数值不变,但"斜对称误差"不再是 0。这说明了什么?Switch the top-right selector to "alternative \(C\)": the numbers \(C\dot{\mathbf q}\) do not change, but the "skew-symmetry error" is no longer 0. What does that tell you?
离心项和科氏项不是新的物理力,而是用关节角作坐标描述牛顿定律时多出来的"速度乘速度"项:关节角是曲线坐标,\(M\) 随坐标变化,\(\dot{\mathbf q}\) 变化的同时尺子自己也在变,\(\tfrac{d}{dt}(M\dot{\mathbf q})\) 里便冒出了 \(\dot M\dot{\mathbf q}\)。如果手臂的某个姿态变化不改变 \(M\)(\(h=0\),伸直或折回时),这些项就消失。
Centrifugal and Coriolis terms are not new physical forces; they are the "velocity-times-velocity" terms that appear when Newton's law is written in joint coordinates. Joint angles are curvilinear coordinates, \(M\) varies with them, and as \(\dot{\mathbf q}\) changes the ruler itself changes: a \(\dot M\dot{\mathbf q}\) pops out of \(\tfrac{d}{dt}(M\dot{\mathbf q})\). Where a posture change leaves \(M\) unchanged (\(h=0\), arm straight or folded), these terms vanish.
4第二、第三视角:工具端"感觉"有多重Second and third views: how heavy does the tool feel?
在关节空间里 \(M\) 说得清楚,可是人和任务关心的是手爪端:用手推一下末端,它"像多重的东西"?在方向 \(\mathbf u\) 上推一个力 \(f\),手爪的加速度由 \(\ddot{\mathbf x}=J\ddot{\mathbf q}\)(忽略速度项与重力)和 \(\ddot{\mathbf q}=M^{-1}J^\top\mathbf F\)(3.4 节的 \(\boldsymbol\tau=J^\top\mathbf F\))给出:
\(M\) is crisp in joint space, but people and tasks care about the tool tip: push on it, and what mass does it feel like? Push with force \(\mathbf F\) in direction \(\mathbf u\): the tip acceleration follows from \(\ddot{\mathbf x}=J\ddot{\mathbf q}\) (velocity terms and gravity ignored) and \(\ddot{\mathbf q}=M^{-1}J^\top\mathbf F\) (the \(\boldsymbol\tau=J^\top\mathbf F\) of Lesson 3.4):
\(\Lambda\) 叫操作空间惯量(Khatib,1987),\(m(\mathbf u)\) 是沿 \(\mathbf u\) 推时的等效质量:推 1 N 沿 \(\mathbf u\) 方向得到的加速度分量是 \(1/m\)。它和 3.4 节的可操作度椭圆是亲戚:若 \(M=I\),\(\Lambda^{-1}=JJ^\top\),等效质量就是 \(1/\sigma_i^2\),手爪"容易动"的方向(椭圆长轴)等效质量小,"难动"的方向(短轴)等效质量大。一般的 \(M\) 只是把关节速度的"尺子"换成了动能:速度椭圆变成"单位动能 \(\tfrac12\dot{\mathbf x}^\top\Lambda\dot{\mathbf x}=\tfrac12\) 的手爪速度集合"。
在下图拖动手臂,再拖动手爪上的琥珀色小圆改变推的方向(同时画出加速度响应)。右图是极坐标图,半径 = 该方向上的等效质量。
\(\Lambda\) is called the operational-space inertia (Khatib, 1987) and \(m(\mathbf u)\) the apparent mass along \(\mathbf u\): pushing with 1 N along \(\mathbf u\) accelerates the tip by \(1/m\) in that direction. It is a relative of the manipulability ellipse of Lesson 3.4: if \(M=I\), \(\Lambda^{-1}=JJ^\top\) and the apparent masses are \(1/\sigma_i^2\), so directions where the tip moves easily (the long axis) are light, and stiff directions (the short axis) are heavy. A general \(M\) only swaps the "ruler" on joint velocities for kinetic energy: the velocity ellipse becomes the set of tip velocities with kinetic energy \(\tfrac12\dot{\mathbf x}^\top\Lambda\dot{\mathbf x}=\tfrac12\).
Drag the arm in the figure, then drag the small amber circle at the tip to change the push direction (the acceleration response is drawn too). The right panel is a polar plot whose radius is the apparent mass in that direction.
- 找出让等效质量最各向同性的姿态(\(\kappa(\Lambda)<4.6\))。它是"伸直"还是"弯成约 90°"?能不能做到 \(\kappa=1\)?Find the posture where the apparent mass is the most isotropic (\(\kappa(\Lambda)<4.6\)). Is it "stretched" or "bent around 90°"? Can you reach \(\kappa=1\)?
- 把手臂拉到接近伸直,让最大的等效质量超过 20 kg。手臂总共才 3.5 kg,这 20 kg 是哪里来的?Stretch the arm nearly straight until the largest apparent mass exceeds 20 kg. The whole arm weighs only 3.5 kg, so where do the 20 kg come from?
矩阵视角:\(M=\partial^2T/\partial\dot{\mathbf q}^2\) 就是动能对速度的 Hessian(3.5 节),因为 \(T\) 是二次的,Hessian 就是 \(M\) 本身。几何视角:\(M(\mathbf q)\) 给关节空间定义了一个随位置变化的"长度尺子" \(\|\dot{\mathbf q}\|_M^2=\dot{\mathbf q}^\top M\dot{\mathbf q}\)(黎曼度量)。不受力、无重力时手臂走的是这个度量下的"直线"(测地线),克里斯托费尔符号就是描述"直线"如何弯曲的联络,所以 \(C\) 恰好由 \(\partial M/\partial\mathbf q\) 构造(第 5 节)。操作空间视角:把这把尺子经 \(J\) 投影到手爪端,就是 \(\Lambda\)。
Matrix view: \(M=\partial^2T/\partial\dot{\mathbf q}^2\) is the Hessian of the kinetic energy with respect to the velocity (Lesson 3.5); since \(T\) is quadratic, that Hessian is \(M\) itself. Geometric view: \(M(\mathbf q)\) gives joint space a position-dependent "length ruler" \(\|\dot{\mathbf q}\|_M^2=\dot{\mathbf q}^\top M\dot{\mathbf q}\) (a Riemannian metric). With no torque and no gravity the arm follows "straight lines" (geodesics) of this metric, and the Christoffel symbols are the connection describing how those lines bend, which is why \(C\) is built from \(\partial M/\partial\mathbf q\) (Section 5). Operational-space view: project that ruler through \(J\) to the tool tip and you get \(\Lambda\).
5再看数学The math
提示:把鼠标悬停在公式中带底色的符号上,对应的图元素会被高亮(FIG 1 的椭圆、FIG 2 的条形图、FIG 3 的等效质量椭圆)。Tip: hover the highlighted symbols in the formulas; the matching element lights up (FIG 1's ellipse, FIG 2's bars, FIG 3's apparent-mass ellipse).
从动能得到 \(M\),并证明它对称正定The mass matrix from kinetic energy, and why it is symmetric positive definite
第 \(i\) 根连杆的质心速度与角速度都是关节速度的线性函数:\(\mathbf v_{ci}=J_{vi}(\mathbf q)\dot{\mathbf q}\),\(\boldsymbol\omega_i=J_{\omega i}(\mathbf q)\dot{\mathbf q}\)(7.4 节的雅可比,对象换成各连杆的质心)。刚体动能是平动加转动:
The centre-of-mass velocity and angular velocity of link \(i\) are linear in the joint velocities: \(\mathbf v_{ci}=J_{vi}(\mathbf q)\dot{\mathbf q}\), \(\boldsymbol\omega_i=J_{\omega i}(\mathbf q)\dot{\mathbf q}\) (the Jacobians of Lesson 7.4, taken at each link's centre of mass). A rigid body's kinetic energy is translation plus rotation:
其中 \({}^{0}I_i=R_iI_iR_i^\top\) 是连杆 \(i\) 的惯量张量在基座系中的表示。方括号就是 \(M(\mathbf q)\)。用 5.10 节的 twist 写,它更紧凑:\(M=\sum_i J_{bi}^\top\mathcal G_iJ_{bi}\),\(J_{bi}\) 是连杆 \(i\) 质心系的 \(6\times n\) 体雅可比,\(\mathcal G_i=\mathrm{diag}(I_i,\,m_i\mathbf 1_3)\) 是空间惯量(角部分在前,和本课程的 twist 顺序一致)。
where \({}^{0}I_i=R_iI_iR_i^\top\) is link \(i\)'s inertia tensor expressed in the base frame. The bracket is \(M(\mathbf q)\). With the twists of Lesson 5.10 it is more compact: \(M=\sum_i J_{bi}^\top\mathcal G_iJ_{bi}\), with \(J_{bi}\) the \(6\times n\) body Jacobian of link \(i\)'s centre-of-mass frame and \(\mathcal G_i=\mathrm{diag}(I_i,\,m_i\mathbf 1_3)\) the spatial inertia (angular part first, matching the course's twist order).
对称正定的证明。每一项 \(m_iJ_{vi}^\top J_{vi}\) 与 \(J_{\omega i}^\top{}^{0}I_iJ_{\omega i}\) 都是对称半正定的(\(m_i>0\),\(I_i\) 正定)。对任何 \(\dot{\mathbf q}\ne0\),\(T=0\) 要求每根连杆的质心速度和角速度都为零,也就是整台机器人静止,而关节速度非零必然使至少一根连杆在动(链式结构中第一个有非零关节速度的连杆本身就在转或移动)。所以 \(\dot{\mathbf q}^\top M\dot{\mathbf q}=2T>0\),\(M\) 正定。又因为对转动关节的臂,\(M\) 的元素是 \(\mathbf q\) 的有界连续(三角)函数,存在与姿态无关的常数 \(0<\underline m\le\lambda_{\min}(M)\le\lambda_{\max}(M)\le\overline m\):这是控制器稳定性证明里反复用到的"有界性"。
Why \(M\) is symmetric positive definite. Every term \(m_iJ_{vi}^\top J_{vi}\) and \(J_{\omega i}^\top{}^{0}I_iJ_{\omega i}\) is symmetric positive semi-definite (\(m_i>0\), \(I_i\) positive definite). For any \(\dot{\mathbf q}\ne0\), \(T=0\) would require every link's centre-of-mass and angular velocity to vanish, i.e. the whole robot to be at rest, but a non-zero joint velocity necessarily moves at least one link (the first link with a non-zero joint rate is itself rotating or sliding). Hence \(\dot{\mathbf q}^\top M\dot{\mathbf q}=2T>0\) and \(M\) is positive definite. For a revolute-joint arm the entries of \(M\) are bounded continuous (trigonometric) functions of \(\mathbf q\), so there are posture-independent constants \(0<\underline m\le\lambda_{\min}(M)\le\lambda_{\max}(M)\le\overline m\): the "boundedness" used again and again in controller stability proofs.
对本课的平面两连杆臂,展开 \(T\)(质心 2 的速度平方为 \(l_1^2\dot\theta_1^2+l_{c2}^2(\dot\theta_1+\dot\theta_2)^2+2l_1l_{c2}\cos\theta_2\,\dot\theta_1(\dot\theta_1+\dot\theta_2)\))并读出系数:
For this lesson's planar two-link arm, expand \(T\) (the squared speed of centre 2 is \(l_1^2\dot\theta_1^2+l_{c2}^2(\dot\theta_1+\dot\theta_2)^2+2l_1l_{c2}\cos\theta_2\,\dot\theta_1(\dot\theta_1+\dot\theta_2)\)) and read off the coefficients:
只有 \(\cos\theta_2\) 出现:\(\theta_2=0\)(伸直)时 \(M_{11}\) 最大;\(\theta_2=\pi\)(折回)时最小,FIG 1 下方的曲线就是这三个式子。\(M_{22}\) 是常数:关节 2 单独转动时,连杆 2 绕关节 2 的转动惯量不随姿态变,但是它和关节 1 的耦合 \(M_{12}\) 变了。
Only \(\cos\theta_2\) appears: \(M_{11}\) is largest at \(\theta_2=0\) (straight) and smallest at \(\theta_2=\pi\) (folded); the curves under FIG 1 are exactly these three formulas. \(M_{22}\) is constant: link 2's inertia about joint 2 does not change with posture, but its coupling \(M_{12}\) to joint 1 does.
克里斯托费尔符号:从 \(\partial M/\partial\mathbf q\) 造出 \(C\)Christoffel symbols: building \(C\) from \(\partial M/\partial\mathbf q\)
8.1 节的欧拉–拉格朗日方程 \(\tfrac{d}{dt}\tfrac{\partial T}{\partial\dot q_k}-\tfrac{\partial T}{\partial q_k}+\tfrac{\partial V}{\partial q_k}=\tau_k\) 对二次型 \(T=\tfrac12\sum M_{ij}\dot q_i\dot q_j\) 求导,得到
Apply Lesson 8.1's Euler–Lagrange equation \(\tfrac{d}{dt}\tfrac{\partial T}{\partial\dot q_k}-\tfrac{\partial T}{\partial q_k}+\tfrac{\partial V}{\partial q_k}=\tau_k\) to the quadratic form \(T=\tfrac12\sum M_{ij}\dot q_i\dot q_j\):
括号里的系数对 \(i,j\) 不对称,但 \(\dot q_i\dot q_j\) 对称,所以可以只保留其对称部分,这样得到的系数叫第一类克里斯托费尔符号:
The coefficient in brackets is not symmetric in \(i,j\) while \(\dot q_i\dot q_j\) is, so only its symmetric part matters; the result is the Christoffel symbol of the first kind:
两连杆算例。\(M\) 只依赖 \(q_2=\theta_2\),所以只有 \(\partial M/\partial q_2\) 非零:\(\partial M_{11}/\partial q_2=-2h\),\(\partial M_{12}/\partial q_2=-h\),\(\partial M_{22}/\partial q_2=0\),其中 \(h=m_2l_1l_{c2}\sin\theta_2\)。代入得 \(c_{211}=-h\),\(c_{121}=c_{221}=-h\),\(c_{112}=h\),其余为零,即
Two-link example. \(M\) depends only on \(q_2=\theta_2\), so only \(\partial M/\partial q_2\) is non-zero: \(\partial M_{11}/\partial q_2=-2h\), \(\partial M_{12}/\partial q_2=-h\), \(\partial M_{22}/\partial q_2=0\), with \(h=m_2l_1l_{c2}\sin\theta_2\). Substituting gives \(c_{211}=-h\), \(c_{121}=c_{221}=-h\), \(c_{112}=h\), all others zero, i.e.
这正是 FIG 2 四根条形图的来源:含 \(\dot\theta_i^2\) 的是离心项,含 \(\dot\theta_1\dot\theta_2\) 的是科氏项。代码里的"克里斯托费尔 − 闭式"读数,就是对 \(M\) 数值求偏导后用上面的公式算出 \(C\),与这个闭式的差(约 \(10^{-9}\))。
This is where FIG 2's four bars come from: terms with \(\dot\theta_i^2\) are centrifugal, the term with \(\dot\theta_1\dot\theta_2\) is Coriolis. The "Christoffel − closed form" readout numerically differentiates \(M\), builds \(C\) with the formula above, and compares it with this closed form (the difference is about \(10^{-9}\)).
\(C\) 不唯一,\(C\dot{\mathbf q}\) 唯一;斜对称性与能量守恒\(C\) is not unique, \(C\dot{\mathbf q}\) is; skew-symmetry and energy balance
运动方程里真正出现的只有向量 \(C\dot{\mathbf q}\),它由 \(c_{ijk}\) 的对称部分唯一确定。但"哪个矩阵乘 \(\dot{\mathbf q}\) 得到它"有很多种分法:加上任何满足 \(S\dot{\mathbf q}=0\) 的 \(S\) 都行。例如 \(C_b=\begin{bmatrix}0&-h(2\dot\theta_1+\dot\theta_2)\\ h\dot\theta_1&0\end{bmatrix}\) 给出同样的 \(C\dot{\mathbf q}\)。克里斯托费尔的分法的特别之处是下面的性质:
Only the vector \(C\dot{\mathbf q}\) appears in the equation of motion, and it is fixed uniquely by the symmetric part of \(c_{ijk}\). But there are many ways to write it as "matrix times \(\dot{\mathbf q}\)": add any \(S\) with \(S\dot{\mathbf q}=0\). For example \(C_b=\begin{bmatrix}0&-h(2\dot\theta_1+\dot\theta_2)\\ h\dot\theta_1&0\end{bmatrix}\) gives the same \(C\dot{\mathbf q}\). What is special about the Christoffel choice is this property:
(把 \(\dot M_{kj}=\sum_i\partial_iM_{kj}\dot q_i\) 与 \(2C_{kj}=\sum_i(\partial_iM_{kj}+\partial_jM_{ki}-\partial_kM_{ij})\dot q_i\) 相减即得;交换 \(k\leftrightarrow j\) 并利用 \(M\) 的对称性,符号反号。)斜对称矩阵的二次型为零,所以对任何满足 \(C\dot{\mathbf q}\) 相同的分法都有 \(\dot{\mathbf q}^\top(\dot M-2C)\dot{\mathbf q}=0\),进而得到功率平衡:
(Subtract \(2C_{kj}=\sum_i(\partial_iM_{kj}+\partial_jM_{ki}-\partial_kM_{ij})\dot q_i\) from \(\dot M_{kj}=\sum_i\partial_iM_{kj}\dot q_i\); swapping \(k\leftrightarrow j\) and using the symmetry of \(M\) flips the sign.) A skew-symmetric matrix has a vanishing quadratic form, so \(\dot{\mathbf q}^\top(\dot M-2C)\dot{\mathbf q}=0\) for every factorization with the same \(C\dot{\mathbf q}\), and we get the power balance:
科氏力与离心力不做功:它们只在关节之间重新分配动能,总动能的变化只来自电机输出功率 \(\dot{\mathbf q}^\top\boldsymbol\tau\) 和重力做功。(对 \(n=2\),克里斯托费尔分法是同时满足"给出 \(C\dot{\mathbf q}\)"和"\(\dot M-2C\) 斜对称"的唯一矩阵;\(n\ge3\) 时还有别的。自适应控制的证明用到的是斜对称性,所以工程上默认用克里斯托费尔分法,Pinocchio 的 computeCoriolisMatrix 给出的 \(C\) 也是如此。)
Coriolis and centrifugal terms do no work: they only shuffle kinetic energy between joints, and the total changes only through the motor power \(\dot{\mathbf q}^\top\boldsymbol\tau\) and gravity. (For \(n=2\) the Christoffel factorization is the unique matrix with both properties, "gives \(C\dot{\mathbf q}\)" and "\(\dot M-2C\) skew"; for \(n\ge3\) there are others. Adaptive-control proofs use the skew-symmetry, so the Christoffel choice is the engineering default, and the \(C\) returned by Pinocchio's computeCoriolisMatrix follows it.)
重力项与"对参数线性"The gravity term and "linear in the parameters"
重力来自势能:\(V=\sum_i m_ig\,h_i(\mathbf q)\)(\(h_i\) 是质心高度),\(\mathbf g(\mathbf q)=\partial V/\partial\mathbf q\)。平面臂里 \(g_1=(m_1l_{c1}+m_2l_1)g\cos\theta_1+m_2l_{c2}g\cos(\theta_1+\theta_2)\),\(g_2=m_2l_{c2}g\cos(\theta_1+\theta_2)\):每一项都是 \(\partial(\text{高度})/\partial\theta\) 乘以重量,也就是重力对该关节的力臂。把前面所有的式子合起来,每一项都是"一个只和 \((\mathbf q,\dot{\mathbf q},\ddot{\mathbf q})\) 有关的已知函数 × 一个只含物理参数的数",所以整个方程可以写成
Gravity comes from the potential energy: \(V=\sum_im_ig\,h_i(\mathbf q)\) (\(h_i\) the height of centre \(i\)) and \(\mathbf g(\mathbf q)=\partial V/\partial\mathbf q\). For the planar arm \(g_1=(m_1l_{c1}+m_2l_1)g\cos\theta_1+m_2l_{c2}g\cos(\theta_1+\theta_2)\) and \(g_2=m_2l_{c2}g\cos(\theta_1+\theta_2)\): each term is a weight times the lever arm of gravity about that joint. Putting all the formulas together, every term is "a known function of \((\mathbf q,\dot{\mathbf q},\ddot{\mathbf q})\) times a number made of physical parameters", so the whole equation can be written
\(Y\) 叫回归矩阵。两根连杆有 6 个惯性参数(\(m_i,l_{ci},I_i\)),但只有 5 个组合 \(a,\dots,e\) 能被力矩数据区分,称为基参数(Gautier–Khalil,1990)。\(\boldsymbol\pi\) 的辨识就是一个线性最小二乘问题(1.13 节):把 \(N\) 个时刻的 \(Y_k\) 叠成大矩阵 \(A\),解 \(\hat{\boldsymbol\pi}=\arg\min\|A\boldsymbol\pi-\boldsymbol\tau\|^2\);\(A\) 的条件数决定了测量噪声被放大多少,因此要选"激励充分"的轨迹(FIG 5)。
\(Y\) is the regressor matrix. The two links have six inertial parameters (\(m_i,l_{ci},I_i\)) but only five combinations \(a,\dots,e\) can be told apart from torque data: the base parameters (Gautier–Khalil, 1990). Identifying \(\boldsymbol\pi\) is a linear least-squares problem (Lesson 1.13): stack the \(Y_k\) of \(N\) instants into a big \(A\) and solve \(\hat{\boldsymbol\pi}=\arg\min\|A\boldsymbol\pi-\boldsymbol\tau\|^2\); the condition number of \(A\) decides how much measurement noise is amplified, hence the need for a "sufficiently exciting" trajectory (FIG 5).
正动力学、逆动力学,以及怎样算 \(M\)Inverse and forward dynamics, and how to compute \(M\)
8.2 节的牛顿–欧拉递推(RNEA)一次 \(O(n)\) 地算出逆动力学,不需要显式的 \(M\) 和 \(C\)。要得到 \(M\) 本身有一个小技巧:令 \(\dot{\mathbf q}=0\)、重力为零,那么 \(\tau=M\ddot{\mathbf q}\),取 \(\ddot{\mathbf q}=\mathbf e_j\),RNEA 的输出就是 \(M\) 的第 \(j\) 列。更高效的是 Featherstone 的组合刚体算法(CRBA),而正动力学则有 \(O(n)\) 的关节体算法(ABA)。偏置项 \(\mathbf b\) 是令 \(\ddot{\mathbf q}=0\) 的 RNEA。MuJoCo 里的 qfrc_bias 就是这个 \(\mathbf b\)。
Lesson 8.2's recursive Newton–Euler algorithm (RNEA) computes inverse dynamics in \(O(n)\) without forming \(M\) or \(C\). To get \(M\) itself there is a small trick: set \(\dot{\mathbf q}=0\) and gravity to zero, so \(\tau=M\ddot{\mathbf q}\); with \(\ddot{\mathbf q}=\mathbf e_j\) the RNEA output is column \(j\) of \(M\). More efficient is Featherstone's composite-rigid-body algorithm (CRBA), and forward dynamics has the \(O(n)\) articulated-body algorithm (ABA). The bias term \(\mathbf b\) is RNEA with \(\ddot{\mathbf q}=0\); MuJoCo's qfrc_bias is exactly this \(\mathbf b\).
带浮动基座的机器人(足式、人形,9.2 节)把广义坐标扩成 \(\mathbf q=(\text{基座位姿},\mathbf q_{\text{joints}})\),\(M\) 变成 \(\begin{bmatrix}M_{bb}&M_{bj}\\M_{jb}&M_{jj}\end{bmatrix}\),基座那 6 行没有电机,\(\boldsymbol\tau\) 的对应行为零,由接触力 \(J_c^\top\mathbf f_c\) 提供:形式不变,只是多了接触。
For floating-base robots (legged, humanoid, Lesson 9.2) the coordinates grow to \(\mathbf q=(\text{base pose},\mathbf q_{\text{joints}})\) and \(M\) becomes \(\begin{bmatrix}M_{bb}&M_{bj}\\M_{jb}&M_{jj}\end{bmatrix}\); the six base rows have no motors, so the corresponding rows of \(\boldsymbol\tau\) are zero and contact forces \(J_c^\top\mathbf f_c\) supply them: the form is unchanged, there is just contact.
计算力矩控制与能量(无源性)论证Computed torque and the energy (passivity) argument
令 \(\mathbf e=\mathbf q_d-\mathbf q\)。计算力矩控制律用模型 \(\hat M,\hat C,\hat{\mathbf g}\) 去抵消手臂:
Let \(\mathbf e=\mathbf q_d-\mathbf q\). The computed-torque law uses a model \(\hat M,\hat C,\hat{\mathbf g}\) to cancel the arm:
模型精确时代入动力学方程,\(M(\ddot{\mathbf q}-\mathbf a)=0\),而 \(M\) 可逆,所以 \(\ddot{\mathbf e}+K_d\dot{\mathbf e}+K_p\mathbf e=0\):误差满足线性、各关节互相解耦的二阶系统,\(K_p=\omega_n^2\mathbf 1\)、\(K_d=2\omega_n\mathbf 1\) 时每个关节都是临界阻尼,指数收敛。模型不精确时,同样的推导给出
With an exact model, substituting into the dynamics gives \(M(\ddot{\mathbf q}-\mathbf a)=0\), and since \(M\) is invertible, \(\ddot{\mathbf e}+K_d\dot{\mathbf e}+K_p\mathbf e=0\): the error obeys a linear, joint-by-joint decoupled second-order system, critically damped in every joint for \(K_p=\omega_n^2\mathbf 1\), \(K_d=2\omega_n\mathbf 1\), converging exponentially. With an imperfect model the same derivation gives
稳定的线性系统被有界的输入 \(\boldsymbol\eta\) 驱动,误差有界,且 \(\|\boldsymbol\eta\|\) 正比于模型误差:这就是"优雅退化"。增大 \(\omega_n\) 能迅速压低误差(本课仿真里大致按 \(1/\omega_n^2\) 下降),但会放大噪声、触及电机饱和与采样延迟,所以实际要靠更好的模型(FIG 5)或自适应、鲁棒控制(Slotine–Li,第四部分)。
A stable linear system driven by the bounded input \(\boldsymbol\eta\) has a bounded error, and \(\|\boldsymbol\eta\|\) scales with the model error: that is "graceful degradation". Raising \(\omega_n\) shrinks the error quickly (roughly as \(1/\omega_n^2\) in this lesson's simulation) but amplifies noise and runs into motor saturation and sampling delay, so in practice one needs a better model (FIG 5) or adaptive and robust control (Slotine–Li, Part IV).
只用重力补偿的 PD 为什么也稳定?用 \(\tilde{\mathbf q}=\mathbf q-\mathbf q^*\),控制律 \(\boldsymbol\tau=\mathbf g(\mathbf q)-K_p\tilde{\mathbf q}-K_d\dot{\mathbf q}\)(\(K_p,K_d\) 正定),取能量函数 \(W=T+\tfrac12\tilde{\mathbf q}^\top K_p\tilde{\mathbf q}\)。由功率平衡,
Why is PD with gravity compensation stable? Take \(\tilde{\mathbf q}=\mathbf q-\mathbf q^*\), the law \(\boldsymbol\tau=\mathbf g(\mathbf q)-K_p\tilde{\mathbf q}-K_d\dot{\mathbf q}\) (\(K_p,K_d\) positive definite) and the energy function \(W=T+\tfrac12\tilde{\mathbf q}^\top K_p\tilde{\mathbf q}\). By the power balance,
能量只减不增;\(\dot W=0\) 时 \(\dot{\mathbf q}=0\),进而 \(\ddot{\mathbf q}=0\),方程给出 \(K_p\tilde{\mathbf q}=0\),所以 LaSalle 不变原理保证 \(\mathbf q\to\mathbf q^*\)。关键一步用到了斜对称性:科氏项在 \(\dot W\) 里被消掉了。这就是"无源性"论证:手臂是一个以 \(\dot{\mathbf q}^\top\boldsymbol\tau\) 为功率输入的耗能系统。
The energy never increases; \(\dot W=0\) forces \(\dot{\mathbf q}=0\), hence \(\ddot{\mathbf q}=0\), and the equation then gives \(K_p\tilde{\mathbf q}=0\), so LaSalle's invariance principle gives \(\mathbf q\to\mathbf q^*\). The key step used the skew-symmetry: the Coriolis term cancelled in \(\dot W\). This is the "passivity" argument: the arm is a system whose power input is \(\dot{\mathbf q}^\top\boldsymbol\tau\) and which does not create energy by itself.
操作空间惯量的推导Deriving the operational-space inertia
令手爪受力 \(\mathbf F\),关节力矩 \(\boldsymbol\tau=J^\top\mathbf F\)。由 \(\ddot{\mathbf x}=J\ddot{\mathbf q}+\dot J\dot{\mathbf q}\) 和 \(\ddot{\mathbf q}=M^{-1}(J^\top\mathbf F-C\dot{\mathbf q}-\mathbf g)\) 消去 \(\ddot{\mathbf q}\):
Let the tool feel a force \(\mathbf F\) and the joints carry \(\boldsymbol\tau=J^\top\mathbf F\). Eliminate \(\ddot{\mathbf q}\) using \(\ddot{\mathbf x}=J\ddot{\mathbf q}+\dot J\dot{\mathbf q}\) and \(\ddot{\mathbf q}=M^{-1}(J^\top\mathbf F-C\dot{\mathbf q}-\mathbf g)\):
手爪端也满足"质量×加速度 + 速度项 + 重力 = 力",只是质量变成了 \(\Lambda\)。它还有一个最小动能的解释:满足 \(J\dot{\mathbf q}=\dot{\mathbf x}\) 的所有关节速度里,动能最小的一个是 \(\dot{\mathbf q}=M^{-1}J^\top\Lambda\dot{\mathbf x}\),最小动能就是 \(\tfrac12\dot{\mathbf x}^\top\Lambda\dot{\mathbf x}\)(FIG 3 左边的椭圆正是它的等值线)。当 \(J\) 不满行秩时(奇异位形,或任务维数大于关节数),\(JM^{-1}J^\top\) 奇异,\(\Lambda\) 在丢失的方向上趋于无穷:该方向上手爪根本无法被关节加速。
The tool end obeys the same form "mass × acceleration + velocity terms + gravity = force", with \(\Lambda\) as the mass. There is also a minimum-energy reading: among all joint velocities with \(J\dot{\mathbf q}=\dot{\mathbf x}\), the one with the least kinetic energy is \(\dot{\mathbf q}=M^{-1}J^\top\Lambda\dot{\mathbf x}\), and that minimum is \(\tfrac12\dot{\mathbf x}^\top\Lambda\dot{\mathbf x}\) (the ellipse on the left of FIG 3 is its level set). When \(J\) loses row rank (a singular posture, or a task with more dimensions than joints), \(JM^{-1}J^\top\) is singular and \(\Lambda\) tends to infinity along the lost direction: the joints simply cannot accelerate the tool that way.
6机器人中的它:计算力矩控制与参数辨识In robotics: computed-torque control and identification
下面是一个真实的仿真:手爪要在 1 秒量级的周期里画一个圆(或 8 字),关节用 RK4 以 1 ms 步长积分,控制器以 1 kHz 更新(零阶保持)。三种控制器共用同一个初始误差:(a) 纯 PD;(b) PD + 重力补偿;(c) 计算力矩。PD 的增益按沿参考轨迹的平均惯量缩放,使三者名义带宽 \(\omega_n\) 相同。"模型质量误差 \(\delta\)"让控制器以为所有质量都是真实值的 \(1+\delta\) 倍,"关节摩擦 \(b\)"是控制器不知道的粘滞摩擦 \(-b\dot{\mathbf q}\)。
把 \(\delta\) 和 \(b\) 都拖到 0,看计算力矩的误差塌缩到仿真精度;再拖开,看它如何"退化"。把"误差"切换成对数坐标,指数收敛是一条直线。
Here is a real simulation: the tool must trace a circle (or a figure-eight) in a period of about a second; the joints are integrated with RK4 at 1 ms and the controller updates at 1 kHz (zero-order hold). Three controllers start from the same initial error: (a) plain PD; (b) PD + gravity compensation; (c) computed torque. The PD gains are scaled by the average inertia along the reference so that all three share the nominal bandwidth \(\omega_n\). The "model mass error \(\delta\)" makes the controller believe every mass is \(1+\delta\) times the true one, and "joint friction \(b\)" is viscous friction \(-b\dot{\mathbf q}\) the controller does not know about.
Drag both \(\delta\) and \(b\) to 0 and watch the computed-torque error collapse to simulation precision; then move away and watch it "degrade". Switch the error plot to log scale: exponential convergence is a straight line.
- 把模型质量误差 \(\delta\) 和摩擦 \(b\) 都调到 0:计算力矩的 RMS 降到 1 mm 以下,而纯 PD 仍超过 5 cm。Set the model mass error \(\delta\) and friction \(b\) both to 0: computed torque's RMS falls below 1 mm while plain PD stays above 5 cm.
- 把模型质量误差调到 \(-30\%\)(模型偏轻)。计算力矩的误差上升了,但仍然比"PD + 重力补偿"小吗?再试 \(+30\%\),结果一样吗?Set the model mass error to \(-30\%\) (model too light). Computed torque's error rises, but is it still below "PD + gravity compensation"? Try \(+30\%\): is the result the same?
- 让模型几乎准确(\(|\delta|\le5\%\),\(b\le0.2\)),把 \(\omega_n\) 降到 7 rad/s 以下:纯 PD 的 RMS 超过 30 cm,计算力矩仍在 3 cm 以内。是谁在做主要的工作?Keep the model nearly right (\(|\delta|\le5\%\), \(b\le0.2\)) and lower \(\omega_n\) below 7 rad/s: plain PD's RMS exceeds 30 cm while computed torque stays within 3 cm. Who is doing most of the work?
模型从哪里来?厂家给的 CAD 质量参数往往不准,负载(夹爪、工件)又会变。方法是辨识:让手臂跑一条"激励充分"的轨迹,记录 \((\mathbf q,\dot{\mathbf q},\ddot{\mathbf q},\boldsymbol\tau)\),解线性最小二乘 \(A\boldsymbol\pi\approx\boldsymbol\tau\)。下图中,数据是上面的真实动力学加上高斯力矩噪声。选三种激励:丰富(两个关节都大幅运动)、只动关节 1(关节 2 冻结)、缓慢小幅运动,看 \(A\) 的条件数和各参数的相对误差。
Where does the model come from? Manufacturer CAD masses are often inaccurate, and payloads (grippers, workpieces) change. The answer is identification: run a "sufficiently exciting" trajectory, record \((\mathbf q,\dot{\mathbf q},\ddot{\mathbf q},\boldsymbol\tau)\) and solve the linear least squares \(A\boldsymbol\pi\approx\boldsymbol\tau\). In the figure the data are the true dynamics above plus Gaussian torque noise. Pick one of three excitations, rich (both joints move a lot), joint 1 only (joint 2 frozen), or slow and small, and watch the condition number of \(A\) and the relative error of each parameter.
- 选"丰富"激励,噪声 \(\sigma\ge0.5\) N·m,让五个参数的相对误差都低于 10%。Choose the "rich" excitation with noise \(\sigma\ge0.5\) N·m and get all five relative errors below 10%.
- 选"缓慢小幅"激励,噪声 \(\sigma\ge0.5\) N·m:最差的参数误差超过 100%,\(\kappa(A)\) 是多少?为什么手臂缓慢运动时惯量几乎看不出来?Choose the "slow and small" excitation with noise \(\sigma\ge0.5\) N·m: the worst parameter is off by more than 100%. What is \(\kappa(A)\)? Why are the inertias nearly invisible when the arm moves slowly?
重力补偿与手拖示教:协作臂(UR 的 freedrive、KUKA LBR iiwa、Franka 的手拖模式)读出关节力矩传感器或电流,只补 \(\mathbf g(\mathbf q)\)(加一点摩擦补偿),手臂就能悬浮、被人推着走。力矩控制臂:Franka(libfranka 的 franka::Model 提供 mass、coriolis、gravity)和 iiwa 以 1 kHz 在内环里用 \(M,C,\mathbf g\) 做计算力矩与阻抗控制。动力学库:Pinocchio(crba、rnea、aba、computeCoriolisMatrix)、Drake(CalcMassMatrix、CalcBiasTerm)、MuJoCo(mj_fullM、qfrc_bias)都在做本节的计算。操作空间 / 全身控制:Khatib 的操作空间框架用 \(\Lambda\) 做任务空间控制;四足与人形的全身控制器(第四部分)把 \(M\ddot{\mathbf q}+\mathbf b=S^\top\boldsymbol\tau+J_c^\top\mathbf f_c\) 当作二次规划的等式约束。辨识:Atkeson 等(1986)开创了用最小二乘估计连杆与负载的惯性参数,如今商用臂出厂时和换夹爪后都会做负载辨识。sim-to-real:仿真里的惯量、电机转子惯量、摩擦与真机不符,是策略迁移失败的常见原因,所以强化学习常做动力学随机化(Peng 等,2018);也有人用神经网络直接学 \(M\)、\(C\),同时保持"\(M\) 正定"之类结构(Lutter 等,Deep Lagrangian Networks,2019)。
Gravity compensation and hand-guiding: collaborative arms (UR freedrive, KUKA LBR iiwa, Franka's hand-guiding) read joint torque sensors or motor currents, compensate only \(\mathbf g(\mathbf q)\) (plus a little friction), and the arm floats and follows a push. Torque-controlled arms: the Franka (libfranka's franka::Model exposes mass, coriolis, gravity) and the iiwa use \(M,C,\mathbf g\) in a 1 kHz inner loop for computed-torque and impedance control. Dynamics libraries: Pinocchio (crba, rnea, aba, computeCoriolisMatrix), Drake (CalcMassMatrix, CalcBiasTerm) and MuJoCo (mj_fullM, qfrc_bias) all perform this lesson's computations. Operational-space and whole-body control: Khatib's framework uses \(\Lambda\) for task-space control; whole-body controllers for quadrupeds and humanoids (Part IV) treat \(M\ddot{\mathbf q}+\mathbf b=S^\top\boldsymbol\tau+J_c^\top\mathbf f_c\) as an equality constraint of a quadratic program. Identification: Atkeson et al. (1986) pioneered least-squares estimation of link and payload inertial parameters, and commercial arms run payload identification at the factory and after a gripper change. Sim-to-real: mismatched inertias, motor rotor inertia and friction between simulation and hardware are a common reason policies fail to transfer, which is why reinforcement learning often uses dynamics randomization (Peng et al., 2018); others learn \(M\) and \(C\) with neural networks while keeping structure such as "\(M\) is positive definite" (Lutter et al., Deep Lagrangian Networks, 2019).
7常见误区Common pitfalls
FIG 1 里 \(M\) 的条件数从伸直时的 45 左右降到弯曲时的 5 左右,\(M_{11}\) 差了近三倍。用固定增益的关节 PD,有效带宽 \(\sqrt{K_p/M}\) 随姿态变化:同一组增益在一个姿态下临界阻尼,在另一个姿态下就欠阻尼振荡。直觉在哪里失灵?对工业臂的每个关节单独调参常常"够用",是因为大减速比让电机转子惯量压过了 \(M\) 的变化(见误区 4);换成直驱或低减速比的臂(力矩控制的协作臂、四足的腿)就不行。
In FIG 1 the condition number of \(M\) falls from about 45 when the arm is stretched to about 5 when it is bent, and \(M_{11}\) changes by almost a factor of three. With a fixed-gain joint PD the effective bandwidth \(\sqrt{K_p/M}\) varies with posture: the same gains are critically damped in one posture and ring in another. Where does intuition fail? Tuning each joint of an industrial arm separately often "works" because a large gear ratio lets motor rotor inertia swamp the variation of \(M\) (see Pitfall 4); with direct-drive or low-ratio arms (torque-controlled cobots, legs of quadrupeds) it does not.
这里的科氏项与离心项是坐标的产物:一个不受力的质点,在笛卡尔坐标里 \(\ddot x=0\),一点速度项都没有;换成极坐标,\(\ddot r-r\dot\phi^2=0\) 里立刻冒出"离心项"。机械臂的 \(C\dot{\mathbf q}\) 同理,是用关节角描述运动的代价。判断标准是它们不做功(斜对称性):只改变动能在关节间的分配,不改变总动能。另外它不同于地球自转的科氏力,不要把两者混为一谈。
The Coriolis and centrifugal terms here are artifacts of the coordinates: a force-free particle in Cartesian coordinates has \(\ddot x=0\) with no velocity terms at all; switch to polar coordinates and \(\ddot r-r\dot\phi^2=0\) suddenly shows a "centrifugal term". The arm's \(C\dot{\mathbf q}\) is the same: the price of describing motion in joint angles. The test is that they do no work (skew-symmetry): they redistribute kinetic energy between joints without changing the total. Do not confuse them with the Coriolis force of the Earth's rotation either.
唯一的是向量 \(C\dot{\mathbf q}\),不是矩阵 \(C\)。FIG 2 里两种 \(C\) 给出同样的力矩,却一个满足 \(\dot M-2C\) 斜对称、一个不满足。需要做无源性或自适应控制的证明,或者要用 \(\dot{\mathbf q}^\top(\dot M-2C)\dot{\mathbf q}=0\) 之外的性质(如 \(\dot M=C+C^\top\))时,必须说明用的是克里斯托费尔写法;只算力矩时,任何一种都行,不必纠结。
What is unique is the vector \(C\dot{\mathbf q}\), not the matrix \(C\). In FIG 2 both versions give the same torques, but one makes \(\dot M-2C\) skew-symmetric and the other does not. For passivity or adaptive-control proofs, or anything beyond \(\dot{\mathbf q}^\top(\dot M-2C)\dot{\mathbf q}=0\) (such as \(\dot M=C+C^\top\)), you must say you use the Christoffel version; when you only compute torques, any version will do.
本节的 \(M\) 只含连杆。真实关节还有转子惯量 \(I_m\),经减速比 \(N\) 折算到关节侧是 \(N^2I_m\)("反射惯量"):\(N=100\) 时,它可以比连杆的惯量大一到两个数量级,使 \(M\) 近似为常数对角阵,粘滞、库仑和 Stribeck 摩擦也常常与 \(C\dot{\mathbf q}\) 一样大。FIG 4 里没建模的摩擦会让计算力矩的误差升到毫米以上。真实的辨识要把转子惯量和摩擦也放进回归矩阵(9.4 节)。
This lesson's \(M\) contains only the links. A real joint also has rotor inertia \(I_m\), which a gear ratio \(N\) reflects to the joint side as \(N^2I_m\) ("reflected inertia"): at \(N=100\) it can exceed the link inertia by one or two orders of magnitude, making \(M\) nearly constant and diagonal, and viscous, Coulomb and Stribeck friction are often as large as \(C\dot{\mathbf q}\). In FIG 4 the unmodelled friction pushes computed torque's error above a millimetre. Real identification puts rotor inertia and friction into the regressor too (Lesson 9.4).
FIG 3 中把手臂拉直,\(J\) 奇异,\(\Lambda^{-1}=JM^{-1}J^\top\) 出现零特征值,沿丢失方向的 \(\Lambda\) 趋于无穷。这不是数值问题,而是物理事实:那个方向上手爪根本动不了。写代码时别直接对 \(JM^{-1}J^\top\) 求逆,而要先检查 \(\sigma_{\min}(J)\),必要时加阻尼或降维。冗余臂(\(J\) 行满秩而列多)没有这个问题,但关节空间的逆应该用动力学一致的广义逆 \(\bar J=M^{-1}J^\top\Lambda\),而不是 \(J^+\)。还有,位置和姿态混在一个 6 维 \(\Lambda\) 里时单位是 kg 与 kg·m²,对它求特征值没有意义(同 3.4 节的误区 4)。
In FIG 3 stretch the arm and \(J\) becomes singular: \(\Lambda^{-1}=JM^{-1}J^\top\) gets a zero eigenvalue and \(\Lambda\) tends to infinity along the lost direction. That is not a numerical glitch but a physical fact: the tool cannot move that way at all. In code, do not blindly invert \(JM^{-1}J^\top\); check \(\sigma_{\min}(J)\) first and damp or reduce the task dimension if needed. A redundant arm (\(J\) of full row rank with extra columns) has no such problem, but the joint-space inverse should be the dynamically consistent generalized inverse \(\bar J=M^{-1}J^\top\Lambda\), not \(J^+\). And when position and orientation share one 6D \(\Lambda\), the units are kg and kg·m², so its eigenvalues are meaningless (the same as Pitfall 4 of Lesson 3.4).
8小练习Check yourself
9延伸阅读Further reading
- M. W. Spong, S. Hutchinson & M. Vidyasagar, Robot Modeling and Control, 2nd ed., Wiley, 2020, Ch. 7–8 — 动力学(质量矩阵、克里斯托费尔符号、斜对称性、线性参数化)与多变量控制(逆动力学控制)Dynamics (mass matrix, Christoffel symbols, skew-symmetry, linear parameterization) and multivariable control (inverse-dynamics control)
- J.-J. E. Slotine & W. Li, Applied Nonlinear Control, Prentice Hall, 1991, Ch. 9 — 机器人动力学的性质、无源性证明、自适应控制Properties of robot dynamics, passivity proofs, adaptive control
- R. M. Murray, Z. Li & S. S. Sastry, A Mathematical Introduction to Robotic Manipulation, CRC Press, 1994, Ch. 4 — 机械臂动力学与控制,用旋量写法推导 \(M\)、\(C\)Manipulator dynamics and control, with \(M\) and \(C\) derived in screw notation
- K. Lynch & F. Park, Modern Robotics, Ch. 8 and Ch. 11 — 开链动力学(拉格朗日、牛顿–欧拉、空间惯量)与机器人控制(计算力矩、任务空间运动控制)Dynamics of open chains (Lagrangian, Newton–Euler, spatial inertia) and robot control (computed torque, task-space motion control)
- B. Siciliano, L. Sciavicco, L. Villani & G. Oriolo, Robotics: Modelling, Planning and Control, Springer, 2009, Ch. 7–8 — 动力学(含参数线性化)与运动控制Dynamics (including linear parameterization) and motion control
- O. Khatib, A Unified Approach for Motion and Force Control of Robot Manipulators: The Operational Space Formulation, IEEE J. Robotics and Automation 3(1), 1987 — 操作空间惯量 \(\Lambda\) 的出处The origin of the operational-space inertia \(\Lambda\)
- R. Featherstone, Rigid Body Dynamics Algorithms, Springer, 2008 — RNEA、CRBA、ABA 的权威参考The reference for RNEA, CRBA and ABA
- C. G. Atkeson, C. H. An & J. M. Hollerbach, Estimation of Inertial Parameters of Manipulator Loads and Links, Int. J. Robotics Research 5(3), 1986; M. Gautier & W. Khalil, Direct Calculation of Minimum Set of Inertial Parameters of Serial Robots, IEEE Trans. Robotics and Automation 6(3), 1990 — 惯性参数辨识与基参数Inertial-parameter identification and base parameters
- J.-J. E. Slotine & W. Li, On the Adaptive Control of Robot Manipulators, Int. J. Robotics Research 6(3), 1987 — 用斜对称性证明的自适应控制Adaptive control proved with the skew-symmetry property
- M. Lutter, C. Ritter & J. Peters, Deep Lagrangian Networks: Using Physics as Model Prior for Deep Learning, ICLR 2019 — 用神经网络学习保持物理结构的 \(M\)、\(C\)Learning \(M\) and \(C\) with neural networks while preserving physical structure