主要参考学习资料及插图来源:
《机器人学导论(第4版)》John J.Craig著
台大机器人学之动力学——林沛群
前置知识:线性代数-理论力学
码字不易,求点赞收藏加关注(´•ω•̥`)
有问题欢迎评论区讨论~
目录
- 操作臂动力学
- 刚体的加速度
- 线加速度
- 角加速度
- 质量分布
- 牛顿-欧拉方程递推动力学方程
- 牛顿方程和欧拉方程
- 计算速度和加速度的外推法
- 作用在连杆上的力和力矩
- 计算力和力矩的内推法
- 牛顿-欧拉递推动力学算法
- 考虑重力的动力学算法
- 操作臂动力学方程的结构
- 迭代形式与封闭形式
- 状态空间方程
- 位形空间方程
- 操作臂动力学的拉格朗日方程
- 笛卡尔空间中的操作臂动力学
- 笛卡尔状态空间方程
- 笛卡尔位形空间中的力矩方程
- 考虑非刚体影响
操作臂动力学
刚体的加速度
对刚体的线速度和角速度进行求导分别得到线加速度和角加速度:
B V ˙ Q = d d t B V Q = lim Δ t → 0 B V Q ( t + Δ t ) − B V Q ( t ) Δ t ^B\dot{V}_Q=\frac{\mathrm{d}}{\mathrm{d}t}^BV_Q=\lim_{\Delta t\rightarrow0}\frac{^BV_Q(t+\Delta t)-^BV_Q(t)}{\Delta t} BV˙Q=dtdBVQ=limΔt→0ΔtBVQ(t+Δt)−BVQ(t)
A Ω ˙ B = d d t A Ω B = lim Δ t → 0 A Ω B ( t + Δ t ) − A Ω B ( t ) Δ t ^A\dot{\Omega}_B=\frac{\mathrm{d}}{\mathrm{d}t}^A\Omega_B=\lim_{\Delta t\rightarrow0}\frac{^A\Omega_B(t+\Delta t)-^A\Omega_B(t)}{\Delta t} AΩ˙B=dtdAΩB=limΔt→0ΔtAΩB(t+Δt)−AΩB(t)
当刚体所处的瞬时参考坐标系为世界坐标系 { U } \{U\} {U}时,可用下列符号表示:
v ˙ a = U V ˙ A O R G \dot{v}_a=^U\dot{V}_AORG v˙a=UV˙AORG
ω ˙ A = U Ω ˙ A \dot{\omega}_A=^U\dot{\Omega}_A ω˙A=UΩ˙A
线加速度
应用第五章所得公式:
A V Q = A V B O R G + B A R B V Q + A Ω B × B A R B Q ^AV_Q=^AV_{BORG}+^A_BR^BV_Q+^A\Omega_B\times^A_BR^BQ AVQ=AVBORG+BARBVQ+AΩB×BARBQ
当 { A } \{A\} {A}与 { B } \{B\} {B}原点重合:
A V Q = d d t ( B A R B Q ) = B A R B V Q + A Ω B × B A R B Q ^AV_Q=\frac{\mathrm{d}}{\mathrm{d}t}(^A_BR^BQ)=^A_BR^BV_Q+^A\Omega_B\times^A_BR^BQ AVQ=dtd(BARBQ)=BARBVQ+AΩB×BARBQ①
对上式求导:
A V ˙ Q = d d t ( B A R B Q ) + A Ω ˙ B × B A R B Q + A Ω B × d d t ( B A R B Q ) ^A\dot{V}_Q=\frac{\mathrm{d}}{\mathrm{d}t}(^A_BR^BQ)+^A\dot{\Omega}_B\times^A_BR^BQ+^A\Omega_B\times\frac{\mathrm{d}}{\mathrm{d}t}(^A_BR^BQ) AV˙Q=dtd(BARBQ)+AΩ˙B×BARBQ+AΩB×dtd(BARBQ)
= B A R B V ˙ Q + A Ω B × B A R B V Q + A Ω ˙ B × B A R B Q + A Ω B × ( B A R B V Q + A Ω B × B A R B Q ) =^A_BR^B\dot{V}_Q+^A\Omega_B\times^A_BR^BV_Q+^A\dot{\Omega}_B\times^A_BR^BQ+^A\Omega_B\times(^A_BR^BV_Q+^A\Omega_B\times^A_BR^BQ) =BARBV˙Q+AΩB×BARBVQ+AΩ˙B×BARBQ+AΩB×(BARBVQ+AΩB×BARBQ)
= B A R B V ˙ Q + 2 A Ω B × B A R B V Q + A Ω ˙ B × B A R B Q + A Ω B × ( A Ω B × B A R B Q ) =^A_BR^B\dot{V}_Q+2^A\Omega_B\times^A_BR^BV_Q+^A\dot{\Omega}_B\times^A_BR^BQ+^A\Omega_B\times(^A\Omega_B\times^A_BR^BQ) =BARBV˙Q+2AΩB×BARBVQ+AΩ˙B×BARBQ+AΩB×(AΩB×BARBQ)
再考虑原点不重合的情况,加上 { B } \{B\} {B}原点的加速度:
A V ˙ Q = A V ˙ B O R G + B A R B V ˙ Q + 2 A Ω B × B A R B V Q + A Ω ˙ B × B A R B Q + A Ω B × ( A Ω B × B A R B Q ) ^A\dot{V}_Q=^A\dot{V}_{BORG}+^A_BR^B\dot{V}_Q+2^A\Omega_B\times^A_BR^BV_Q+^A\dot{\Omega}_B\times^A_BR^BQ+^A\Omega_B\times(^A\Omega_B\times^A_BR^BQ) AV˙Q=AV˙BORG+BARBV˙Q+2AΩB×BARBVQ+AΩ˙B×BARBQ+AΩB×(AΩB×BARBQ)②
上式为求解加速度的一般公式。
当计算旋转关节操作臂连杆线加速度时, B Q ^BQ BQ为常量,即:
B V Q = B V ˙ Q = 0 ^BV_Q=^B\dot{V}_Q=0 BVQ=BV˙Q=0
此时②简化为:
A V ˙ Q = A V ˙ B O R G + A Ω ˙ B × B A R B Q + A Ω B × ( A Ω B × B A R B Q ) ^A\dot{V}_Q=^A\dot{V}_{BORG}+^A\dot{\Omega}_B\times^A_BR^BQ+^A\Omega_B\times(^A\Omega_B\times^A_BR^BQ) AV˙Q=AV˙BORG+AΩ˙B×BARBQ+AΩB×(AΩB×BARBQ)③
角加速度
假设 { B } \{B\} {B}以 A Ω B ^A\Omega_B AΩB相对于 { A } \{A\} {A}转动,同时 { C } \{C\} {C}以 B Ω C ^B\Omega_C BΩC相对于 { B } \{B\} {B}转动,则:
A Ω C = A Ω B + B A R B Ω C ^A\Omega_C=^A\Omega_B+^A_BR^B\Omega_C AΩC=AΩB+BARBΩC
对上式求导:
A Ω ˙ C = A Ω ˙ B + d d t ( B A R B Ω C ) ^A\dot{\Omega}_C=^A\dot{\Omega}_B+\frac{\mathrm{d}}{\mathrm{d}t}(^A_BR^B\Omega_C) AΩ˙C=AΩ˙B+dtd(BARBΩC)
最后一项用 B Ω C ^B\Omega_C BΩC代入①中 B Q ^BQ BQ:
A Ω ˙ C = A Ω ˙ B + B A R B Ω ˙ C + A Ω B × B A R B Ω C ^A\dot{\Omega}_C=^A\dot{\Omega}_B+^A_BR^B\dot{\Omega}_C+^A\Omega_B\times^A_BR^B\Omega_C AΩ˙C=AΩ˙B+BARBΩ˙C+AΩB×BARBΩC④
上式用于计算操作臂连杆的角加速度。
质量分布
对一个可以在三维空间自由移动的刚体来说可能存在无穷个旋转轴。在一个刚体绕任意轴做旋转运动时,我们需要一种能够表征刚体质量分布的方式。我们引入惯性张量,可以看做是对一个物体惯性矩的广义度量。
惯性张量可以在任何坐标系中定义,但一般在刚体自身坐标系中定义,用左上标表明已知惯性张量所在的坐标系。 { A } \{A\} {A}中的惯性张量可用 3 × 3 3\times3 3×3矩阵表示如下:
A I = [ I x x − I x y − I x z − I x y I y y − I y z − I x z − I y z I z z ] ^AI=\begin{bmatrix}I_{xx}&-I_{xy}&-I_{xz}\\-I_{xy}&I_{yy}&-I_{yz}\\-I_{xz}&-I_{yz}&I_{zz}\end{bmatrix} AI= Ixx−Ixy−Ixz−IxyIyy−Iyz−Ixz−IyzIzz
矩阵中各元素为:
I x x = ∭ V ( y 2 + z 2 ) ρ d v I_{xx}=\iiint_V(y^2+z^2)\rho\mathrm{d}v Ixx=∭V(y2+z2)ρdv
I y y = ∭ V ( x 2 + z 2 ) ρ d v I_{yy}=\iiint_V(x^2+z^2)\rho\mathrm{d}v Iyy=∭V(x2+z2)ρdv
I z z = ∭ V ( x 2 + y 2 ) ρ d v I_{zz}=\iiint_V(x^2+y^2)\rho\mathrm{d}v Izz=∭V(x2+y2)ρdv
I x y = ∭ V x y ρ d v I_{xy}=\iiint_Vxy\rho\mathrm{d}v Ixy=∭Vxyρdv
I x z = ∭ V x z ρ d v I_{xz}=\iiint_Vxz\rho\mathrm{d}v Ixz=∭Vxzρdv
I y z = ∭ V y z ρ d v I_{yz}=\iiint_Vyz\rho\mathrm{d}v Iyz=∭Vyzρdv
式中刚体由单元体 d v \mathrm{d}v dv组成,单元体的密度为 ρ \rho ρ,每个单元体的位置由矢量 A P = ( x , y , z ) T ^AP=(x,y,z)^T AP=(x,y,z)T确定。
I x x I_{xx} Ixx、 I y y I_{yy} Iyy和 I z z I_{zz} Izz称为惯性矩,其余三个交叉项称为惯量积。对一个刚体来说,该六个相互独立的参量取决于所在坐标系的位置和姿态。当选择的坐标系姿态使刚体的惯量积为零时,该坐标系的轴被称为主轴,而相应的惯量矩被称为主惯性矩。
平行移轴定理描述了一个以刚体质心为原点的坐标系平移到另一个坐标系时惯性张量的变换关系。假设 { C } \{C\} {C}是以刚体质心为原点的坐标系, { A } \{A\} {A}为任意平移后的坐标系,则平行移轴定理可以表示为:
A I z z = C I z z + m ( x c 2 + y c 2 ) ^AI_{zz}=^CI_{zz}+m(x^2_c+y^2_c) AIzz=CIzz+m(xc2+yc2)
A I x y = C I x y − m x c y c ^AI_{xy}=^CI_{xy}-mx_cy_c AIxy=CIxy−mxcyc
式中矢量 P c = ( x c , y c , z c ) T P_c=(x_c,y_c,z_c)^T Pc=(xc,yc,zc)T表示刚体质心在 { A } \{A\} {A}中的位置,其余轴下标的变换公式同理。
平行移轴定理的矢量-矩阵形式为:
A I = C I + m ( P c T P c I 3 − P c P c T ) ^AI=^CI+m(P^T_cP_cI_3-P_cP_c^T) AI=CI+m(PcTPcI3−PcPcT)
式中 I 3 I_3 I3为 3 × 3 3\times3 3×3单位矩阵。
惯性张量的其他性质:
- 如果坐标系的两个坐标轴构成的平面为刚体质量分布的对称平面,则第三个坐标轴与这两个坐标轴的惯性积为零。
- 三个惯量矩的和与参考坐标系的姿态无关。
- 惯性张量的特征值为刚体的主惯性矩,特征矢量为主轴。
大多数操作臂连杆的几何形状及结构组成都比较复杂,因而很难直接应用公式求解,一般使用测量装置来测量每个连杆的惯性矩。
牛顿-欧拉方程递推动力学方程
牛顿方程和欧拉方程
牛顿方程以及描述旋转运动的欧拉方程描述了力、惯量和加速度之间的关系。
牛顿方程
F = m v ˙ c F=m\dot{v}_c F=mv˙c
F F F:作用在质心上的力;
m m m:刚体的总质量;
v ˙ c \dot{v}_c v˙c:刚体质心的加速度。
欧拉方程
N = C I ω ˙ + ω × C I ω N=^CI\dot{\omega}+\omega\times^CI\omega N=CIω˙+ω×CIω
N N N:作用在刚体上的力矩;
C I ^CI CI:刚体在 { C } \{C\} {C}中的惯性张量, { C } \{C\} {C}的原点在质心;
ω \omega ω和 ω ˙ \dot{\omega} ω˙:刚体旋转的角速度和角加速度。
计算速度和加速度的外推法
为了计算作用在连杆上的惯性力,需要计算操作臂每个连杆质心在某一时刻的角速度、线加速度和角加速度。可应用递推方法完成这些计算,首先对连杆 1 1 1进行计算,接着计算下一个连杆,一直外推到连杆 n n n。
第五章机器人连杆的运动中我们给出了角速度在连杆之间的传播方程:
i + 1 ω i + 1 = i i + 1 R i ω i + θ ˙ i + 1 i + 1 Z ^ i + 1 ^{i+1}\omega_{i+1}=^{i+1}_{i}R^i\omega_i+\dot{\theta}_{i+1}{}^{i+1}\hat{Z}_{i+1} i+1ωi+1=ii+1Riωi+θ˙i+1i+1Z^i+1
本章④给出了连杆之间角加速度变换的方程:
i + 1 ω ˙ i + 1 = i i + 1 R i ω ˙ i + i i + 1 R i ω i × θ ˙ i + 1 Z ^ i + 1 + θ ¨ i + 1 i + 1 Z ^ i + 1 ^{i+1}\dot{\omega}_{i+1}=^{i+1}_iR^i\dot{\omega}_i+^{i+1}_iR^i\omega_i\times\dot{\theta}_{i+1}\hat{Z}_{i+1}+\ddot{\theta}_{i+1}{}^{i+1}\hat{Z}_{i+1} i+1ω˙i+1=ii+1Riω˙i+ii+1Riωi×θ˙i+1Z^i+1+θ¨i+1i+1Z^i+1
当第 i + 1 i+1 i+1个关节为移动关节时,上式可简化为:
i + 1 ω ˙ i + 1 = i i + 1 R i ω ˙ i ^{i+1}\dot{\omega}_{i+1}=^{i+1}_{i}R^i\dot{\omega}_i i+1ω˙i+1=ii+1Riω˙i
本章③可得到每个连杆坐标系原点的线加速度:
i + 1 v ˙ i + 1 = i i + 1 R [ i ω ˙ i × i P i + 1 + i ω i × ( i ω i × i P i + 1 ) + i v ˙ i ] ^{i+1}\dot{v}_{i+1}=^{i+1}_iR[^i\dot{\omega}_i\times^iP_{i+1}+^i\omega_i\times(^i\omega_i\times^iP_{i+1})+^i\dot{v}_i] i+1v˙i+1=ii+1R[iω˙i×iPi+1+iωi×(iωi×iPi+1)+iv˙i]
当第 i + 1 i+1 i+1个关节为移动关节时,上式变为②:
i + 1 v ˙ i + 1 = i i + 1 R [ i ω ˙ i × i P i + 1 + i ω i × ( i ω i × i P i + 1 ) + i v ˙ i ] + 2 i + 1 ω i + 1 × d ˙ i + 1 i + 1 Z ^ i + 1 + d ¨ i + 1 i + 1 Z ^ i + 1 ^{i+1}\dot{v}_{i+1}=^{i+1}_iR[^i\dot{\omega}_i\times^iP_{i+1}+^i\omega_i\times(^i\omega_i\times^iP_{i+1})+^i\dot{v}_i]+2^{i+1}\omega_{i+1}\times\dot{d}_{i+1}{}^{i+1}\hat{Z}_{i+1}+\ddot{d}_{i+1}{}^{i+1}\hat{Z}_{i+1} i+1v˙i+1=ii+1R[iω˙i×iPi+1+iωi×(iωi×iPi+1)+iv˙i]+2i+1ωi+1×d˙i+1i+1Z^i+1+d¨i+1i+1Z^i+1
同理用③可得到每个连杆质心的线加速度:
i v ˙ C i = i ω ˙ i × i P C i + i ω i × ( i ω i × i P C i ) + i v ˙ i ^i\dot{v}_{C_i}=^i\dot{\omega}_i\times^iP_{C_i}+^i\omega_i\times(^i\omega_i\times^iP_{C_i})+^i\dot{v}_i iv˙Ci=iω˙i×iPCi+iωi×(iωi×iPCi)+iv˙i
其中 { C i } \{C_i\} {Ci}原点位于连杆质心,姿态与 { i } \{i\} {i}相同。
作用在连杆上的力和力矩
计算出每个连杆质心的线加速度和角加速度之后,运用牛顿-欧拉公式可以计算出作用在连杆质心上的惯性力和力矩:
F i = m v ˙ C i F_i=m\dot{v}_{C_i} Fi=mv˙Ci
N i = C i I ω ˙ i + ω i × C i I ω i N_i=^{C_i}I\dot{\omega}_i+\omega_i\times^{C_i}I\omega_i Ni=CiIω˙i+ωi×CiIωi
计算力和力矩的内推法
计算出每个连杆上的作用力和力矩之后,需要计算这些产生施加在连杆上的力和力矩所对应的关节力矩。
第五章操作臂的静力一节中定义了:
f i = f_i= fi=连杆 i − 1 i-1 i−1作用在连杆 i i i上的力;
n i = n_i= ni=连杆 i − 1 i-1 i−1作用在连杆 i i i上的力矩。
则作用在连杆 i i i上的合力:
i F i = i f i − i + 1 i R i + 1 f i + 1 ^iF_i=^if_i-^i_{i+1}R^{i+1}f_{i+1} iFi=ifi−i+1iRi+1fi+1⑤
作用在连杆质心上的合力矩:
i N i = i n i − i n i + 1 + ( − i P C i ) × i f i − ( i P i + 1 − i P C i ) × i f i + 1 ^iN_i=^in_i-^in_{i+1}+(-^iP_{C_i})\times^if_i-(^iP_{i+1}-^iP_{C_i})\times^if_{i+1} iNi=ini−ini+1+(−iPCi)×ifi−(iPi+1−iPCi)×ifi+1
利用⑤和旋转矩阵变换上式可写成:
i N i = i n i − i + 1 i R i + 1 n i + 1 − i P C i × i F i − i P i + 1 × i + 1 i R i + 1 f i + 1 ^iN_i=^in_i-^i_{i+1}R^{i+1}n_{i+1}-^iP_{C_i}\times^iF_i-^iP_{i+1}\times^i_{i+1}R^{i+1}f_{i+1} iNi=ini−i+1iRi+1ni+1−iPCi×iFi−iPi+1×i+1iRi+1fi+1
最后重新排列力和力矩方程形成相邻连杆从高序号向低序号的迭代关系:
i f i = i + 1 i R i + 1 f i + 1 + i F i ^if_i=^i_{i+1}R^{i+1}f_{i+1}+^iF_i ifi=i+1iRi+1fi+1+iFi
i n i = i N i + i + 1 i R i + 1 n i + 1 + i P C i × i F i + i P i + 1 × i + 1 i R i + 1 f i + 1 ^in_i=^iN_i+^i_{i+1}R^{i+1}n_{i+1}+^iP_{C_i}\times^iF_i+^iP_{i+1}\times^i_{i+1}R^{i+1}f_{i+1} ini=iNi+i+1iRi+1ni+1+iPCi×iFi+iPi+1×i+1iRi+1fi+1
应用这些方程从连杆 n n n开始向内递推到机器人基座依次求解。
最终由第五章操作臂的静力一节给出关节驱动力为:
τ i = i n i T i Z ^ i \tau_i=^in_i^T{}^i\hat{Z}_i τi=iniTiZ^i
对于移动关节:
τ i = i f i T i Z ^ i \tau_i=^if_i^T{}^i\hat{Z}_i τi=ifiTiZ^i
牛顿-欧拉递推动力学算法
由关节运动计算关节力矩的完整算法由两部分组成:第一部分是对每个连杆应用牛顿-欧拉方程从连杆 1 1 1到连杆 n n n向外递推计算连杆的速度和加速度;第二部分是从连杆 n n n到连杆 1 1 1跌倒计算连杆间的相互作用力和力矩以及关节驱动力矩。对于转动关节,该算法归纳如下:
外推:
i + 1 ω i + 1 = i i + 1 R i ω i + θ ˙ i + 1 i + 1 Z ^ i + 1 ^{i+1}\omega_{i+1}=^{i+1}_{i}R^i\omega_i+\dot{\theta}_{i+1}{}^{i+1}\hat{Z}_{i+1} i+1ωi+1=ii+1Riωi+θ˙i+1i+1Z^i+1
i + 1 ω ˙ i + 1 = i i + 1 R i ω ˙ i + i i + 1 R i ω i × θ ˙ i + 1 Z ^ i + 1 + θ ¨ i + 1 i + 1 Z ^ i + 1 ^{i+1}\dot{\omega}_{i+1}=^{i+1}_iR^i\dot{\omega}_i+^{i+1}_iR^i\omega_i\times\dot{\theta}_{i+1}\hat{Z}_{i+1}+\ddot{\theta}_{i+1}{}^{i+1}\hat{Z}_{i+1} i+1ω˙i+1=ii+1Riω˙i+ii+1Riωi×θ˙i+1Z^i+1+θ¨i+1i+1Z^i+1
i + 1 v ˙ i + 1 = i i + 1 R [ i ω ˙ i × i P i + 1 + i ω i × ( i ω i × i P i + 1 ) + i v ˙ i ] ^{i+1}\dot{v}_{i+1}=^{i+1}_iR[^i\dot{\omega}_i\times^iP_{i+1}+^i\omega_i\times(^i\omega_i\times^iP_{i+1})+^i\dot{v}_i] i+1v˙i+1=ii+1R[iω˙i×iPi+1+iωi×(iωi×iPi+1)+iv˙i]
i + 1 v ˙ C i + 1 = i + 1 ω ˙ i + 1 × i + 1 P C i + 1 + i + 1 ω i + 1 × ( i + 1 ω i + 1 × i + 1 P C i + 1 ) + i + 1 v ˙ i + 1 ^{i+1}\dot{v}_{C_{i+1}}=^{i+1}\dot{\omega}_{i+1}\times^{i+1}P_{C_{i+1}}+^{i+1}\omega_{i+1}\times(^{i+1}\omega_{i+1}\times^{i+1}P_{C_{i+1}})+^{i+1}\dot{v}_{i+1} i+1v˙Ci+1=i+1ω˙i+1×i+1PCi+1+i+1ωi+1×(i+1ωi+1×i+1PCi+1)+i+1v˙i+1
F i + 1 = m v ˙ C i + 1 F_{i+1}=m\dot{v}_{C_{i+1}} Fi+1=mv˙Ci+1
N i + 1 = C i + 1 I ω ˙ i + 1 + ω i + 1 × C i + 1 I ω i + 1 N_{i+1}=^{C_{i+1}}I\dot{\omega}_{i+1}+\omega_{i+1}\times^{C_{i+1}}I\omega_{i+1} Ni+1=Ci+1Iω˙i+1+ωi+1×Ci+1Iωi+1
内推:
i f i = i + 1 i R i + 1 f i + 1 + i F i ^if_i=^i_{i+1}R^{i+1}f_{i+1}+^iF_i ifi=i+1iRi+1fi+1+iFi
i n i = i N i + i + 1 i R i + 1 n i + 1 + i P C i × i F i + i P i + 1 × i + 1 i R i + 1 f i + 1 ^in_i=^iN_i+^i_{i+1}R^{i+1}n_{i+1}+^iP_{C_i}\times^iF_i+^iP_{i+1}\times^i_{i+1}R^{i+1}f_{i+1} ini=iNi+i+1iRi+1ni+1+iPCi×iFi+iPi+1×i+1iRi+1fi+1
τ i = i n i T i Z ^ i \tau_i=^in_i^T{}^i\hat{Z}_i τi=iniTiZ^i
考虑重力的动力学算法
令 0 v ˙ 0 = G ^0\dot{v}_0=G 0v˙0=G即可将作用在连杆上的重力因素包括到动力学方程中去,其中 G G G与重力矢量等大反向。
操作臂动力学方程的结构
迭代形式与封闭形式
上述动力学方程主要应用于数值计算(迭代形式)或作为分析方法用于符号方程的推导(封闭形式)。
对于数值计算,只要将待求操作臂的惯性张量、连杆质量、连杆质心位置矢量和相邻连杆的旋转矩阵代入这些方程中即可计算出任何运动情况下的关节力矩。
当需要对方程的结构进行研究,例如重力项的形式、重力和惯性力的影响效果等,就要给出封闭形式的动力学方程。
一个如图所示的RR操作臂,为简单起见假设每个连杆的质量都集中在连杆的末端,应用上述方程得到的封闭形式的动力学方程如下:
τ 1 = m 2 l 2 2 ( θ ¨ 1 + θ ¨ 2 ) + m 2 l 1 l 2 c 2 ( 2 θ ¨ 1 + θ ¨ 2 ) + ( m 1 + m 2 ) l 1 2 θ ¨ 1 − m 2 l 1 l 2 s 2 θ ˙ 2 2 − 2 m 2 l 1 l 2 s 2 θ ˙ 1 θ ˙ 2 + m 2 l 2 g c 12 + ( m 1 + m 2 ) l 1 g c 1 \tau_1=m_2l_2^2(\ddot{\theta}_1+\ddot{\theta}_2)+m_2l_1l_2c_2(2\ddot{\theta}_1+\ddot{\theta}_2)+(m_1+m_2)l^2_1\ddot{\theta}_1-m_2l_1l_2s_2\dot{\theta}_2^2-2m_2l_1l_2s_2\dot{\theta}_1\dot{\theta}_2+m_2l_2gc_{12}+(m_1+m_2)l_1gc_1 τ1=m2l22(θ¨1+θ¨2)+m2l1l2c2(2θ¨1+θ¨2)+(m1+m2)l12θ¨1−m2l1l2s2θ˙22−2m2l1l2s2θ˙1θ˙2+m2l2gc12+(m1+m2)l1gc1
τ 2 = m 2 l 1 l 2 c 2 θ ¨ 1 + m 2 l 1 l 2 s 2 θ ˙ 1 2 + m 2 l 2 g c 12 + m 2 l 2 2 ( θ ¨ 1 + θ ¨ 2 ) \tau_2=m_2l_1l_2c_2\ddot{\theta}_1+m_2l_1l_2s_2\dot{\theta}^2_1+m_2l_2gc_{12}+m_2l^2_2(\ddot{\theta}_1+\ddot{\theta}_2) τ2=m2l1l2c2θ¨1+m2l1l2s2θ˙12+m2l2gc12+m2l22(θ¨1+θ¨2)
如此复杂的函数表达式描述的只是一个最简单的操作臂,可见一个封闭形式的六自由度操作臂的动力学方程是相当复杂的。
状态空间方程
忽略一个方程中的某些细节而仅显示方程的某些结构可以很方便地表示操作臂的动力学方程。
用牛顿-欧拉方程对操作臂进行分析时,动力学方程可写成如下形式:
τ = M ( Θ ) Θ ¨ + V ( Θ , Θ ˙ ) + G ( Θ ) \tau=M(\Theta)\ddot{\Theta}+V(\Theta,\dot{\Theta})+G(\Theta) τ=M(Θ)Θ¨+V(Θ,Θ˙)+G(Θ)
式中 M ( Θ ) M(\Theta) M(Θ)为操作臂的 n × n n\times n n×n质量矩阵, V ( Θ , Θ ˙ ) V(\Theta,\dot{\Theta}) V(Θ,Θ˙)为 n × 1 n\times1 n×1的离心力和科氏力矢量, G ( Θ ) G(\Theta) G(Θ)是 n × 1 n\times1 n×1重力矢量。上式称为状态空间方程,因为矢量 V ( Θ , Θ ˙ ) V(\Theta,\dot{\Theta}) V(Θ,Θ˙)取决于位置和速度。
以上面给出的RR操作臂为例:
τ = [ τ 1 τ 2 ] , Θ ¨ = [ θ ¨ 1 θ ¨ 2 ] \tau=\begin{bmatrix}\tau_1\\\tau_2\end{bmatrix},\ddot{\Theta}=\begin{bmatrix}\ddot{\theta}_1\\\ddot{\theta}_2\end{bmatrix} τ=[τ1τ2],Θ¨=[θ¨1θ¨2]
质量矩阵 M ( Θ ) M(\Theta) M(Θ)的所有各项均为 Θ \Theta Θ的函数并与 Θ ¨ \ddot{\Theta} Θ¨相乘:
M ( Θ ) = [ l 2 2 + 2 l 1 l 2 m 2 c 2 + l 1 2 ( m 1 + m 2 ) l 2 2 m 2 + l 1 l 2 m 2 c 2 l 2 2 m 2 + l 1 l 2 m 2 c 2 l 2 2 m 2 ] M(\Theta)=\begin{bmatrix}l^2_2+2l_1l_2m_2c_2+l_1^2(m_1+m_2)&l^2_2m_2+l_1l_2m_2c_2\\l^2_2m_2+l_1l_2m_2c_2&l^2_2m_2\end{bmatrix} M(Θ)=[l22+2l1l2m2c2+l12(m1+m2)l22m2+l1l2m2c2l22m2+l1l2m2c2l22m2]
速度项 V ( Θ , Θ ˙ ) V(\Theta,\dot{\Theta}) V(Θ,Θ˙)包含了所有与关节速度有关的项:
V ( Θ , Θ ˙ ) = [ − m 2 l 1 l 2 s 2 θ ˙ 2 2 − 2 m 2 l 1 l 2 s 2 θ ˙ 1 θ ˙ 2 m 2 l 1 l 2 s 2 θ ˙ 1 2 ] V(\Theta,\dot{\Theta})=\begin{bmatrix}-m_2l_1l_2s_2\dot{\theta}^2_2-2m_2l_1l_2s_2\dot{\theta}_1\dot{\theta}_2\\m_2l_1l_2s_2\dot{\theta}^2_1\end{bmatrix} V(Θ,Θ˙)=[−m2l1l2s2θ˙22−2m2l1l2s2θ˙1θ˙2m2l1l2s2θ˙12]
− m 2 l 1 l 2 s 2 θ ˙ 1 2 -m_2l_1l_2s_2\dot{\theta}^2_1 −m2l1l2s2θ˙12是与离心力有关的项,因为它是关节速度的平方;
− 2 m 2 l 1 l 2 s 2 θ ˙ 1 θ ˙ 2 -2m_2l_1l_2s_2\dot{\theta}_1\dot{\theta}_2 −2m2l1l2s2θ˙1θ˙2是与科氏力有关的项,因为它总是包含两个不同关节速度的乘积。
重力项 G ( Θ ) G(\Theta) G(Θ)包含了所有与重力加速度 g g g有关的项:
G ( Θ ) = [ m 2 l 2 g c 12 + ( m 1 + m 2 ) l 1 g c 1 m 2 l 2 g c 12 ] G(\Theta)=\begin{bmatrix}m_2l_2gc_{12}+(m_1+m_2)l_1gc_1\\m_2l_2gc_{12}\end{bmatrix} G(Θ)=[m2l2gc12+(m1+m2)l1gc1m2l2gc12]
位形空间方程
将动力学方程中速度项 V ( Θ , Θ ˙ ) V(\Theta,\dot{\Theta}) V(Θ,Θ˙)离心力和科氏力拆开来,并提取出关节速度项,写成另外一种形式如下:
τ = M ( Θ ) Θ ¨ + B ( Θ ) ( Θ ˙ Θ ˙ ) + C ( Θ ) ( Θ ˙ 2 ) + G ( Θ ) \tau=M(\Theta)\ddot{\Theta}+B(\Theta)(\dot{\Theta}\dot{\Theta})+C(\Theta)(\dot{\Theta}^2)+G(\Theta) τ=M(Θ)Θ¨+B(Θ)(Θ˙Θ˙)+C(Θ)(Θ˙2)+G(Θ)⑥
式中 B ( Θ ) B(\Theta) B(Θ)为 n × n ( n − 1 ) / 2 n\times n(n-1)/2 n×n(n−1)/2阶科氏力系数矩阵, ( Θ ˙ Θ ˙ ) (\dot{\Theta}\dot{\Theta}) (Θ˙Θ˙)是 n ( n − 1 ) / 2 × 1 n(n-1)/2\times1 n(n−1)/2×1阶关节速度矢量,即:
( Θ ˙ Θ ˙ ) = ( θ ˙ 1 θ ˙ 2 , θ ˙ 1 θ ˙ 3 , ⋯ , θ ˙ n − 1 θ ˙ n ) T (\dot{\Theta}\dot{\Theta})=(\dot{\theta}_1\dot{\theta}_2,\dot{\theta}_1\dot{\theta}_3,\cdots,\dot{\theta}_{n-1}\dot{\theta}_n)^T (Θ˙Θ˙)=(θ˙1θ˙2,θ˙1θ˙3,⋯,θ˙n−1θ˙n)T
C ( Θ ) C(\Theta) C(Θ)是 n × n n\times n n×n阶离心力系数矩阵, ( Θ ˙ 2 ) (\dot{\Theta}^2) (Θ˙2)是 n × 1 n\times1 n×1阶矢量,即:
( Θ ˙ 2 ) = ( θ ˙ 1 2 , θ ˙ 2 2 , ⋯ , θ ˙ n 2 ) T (\dot{\Theta}^2)=(\dot{\theta}^2_1,\dot{\theta}^2_2,\cdots,\dot{\theta}^2_n)T (Θ˙2)=(θ˙12,θ˙22,⋯,θ˙n2)T
⑥称为位形空间方程,因为它的系数矩阵仅是操作臂位置的函数。
以上面给出的RR操作臂为例:
( Θ ˙ Θ ˙ ) = ( θ ˙ 1 θ ˙ 2 ) , ( Θ ˙ 2 ) = [ θ ˙ 1 2 θ ˙ 2 2 ] (\dot{\Theta}\dot{\Theta})=(\dot{\theta}_1\dot{\theta}_2),(\dot{\Theta}^2)=\begin{bmatrix}\dot{\theta}^2_1\\\dot{\theta}^2_2\end{bmatrix} (Θ˙Θ˙)=(θ˙1θ˙2),(Θ˙2)=[θ˙12θ˙22]
B ( Θ ) = [ − 2 m 2 l 1 l 2 s 2 0 ] B(\Theta)=\begin{bmatrix}-2m_2l_1l_2s_2\\0\end{bmatrix} B(Θ)=[−2m2l1l2s20]
C ( Θ ) = [ 0 − m 2 l 1 l 2 s 2 m 2 l 1 l 2 s 2 0 ] C(\Theta)=\begin{bmatrix}0&-m_2l_1l_2s_2\\m_2l_1l_2s_2&0\end{bmatrix} C(Θ)=[0m2l1l2s2−m2l1l2s20]
操作臂动力学的拉格朗日方程
牛顿-欧拉公式是一种解决动力学问题的力平衡方法,而拉格朗日公式则是一种基于能量的动力学方法。对于同一个操作臂来说,两种方法得到的动力学方程式相同的。
第 i i i个连杆的动能 k i k_i ki可以表示为:
k i = 1 2 m i v C i T v C i + 1 2 i ω i T C i I i i ω i k_i=\frac{1}{2}m_iv^T_{C_i}v_{C_i}+\frac{1}{2}{}^i\omega^T_i{}^{C_i}I_i{}^i\omega_i ki=21mivCiTvCi+21iωiTCiIiiωi
式中第一项是连杆质心线速度动能,第二项是连杆角速度动能。整个操作臂的动能是各个连杆动能之和:
k = Σ i = 1 n k i k=\Sigma^n_{i=1}k_i k=Σi=1nki
式中 v C i v_{C_i} vCi和 i ω i ^i\omega_i iωi是 Θ \Theta Θ和 Θ ˙ \dot{\Theta} Θ˙的函数,因此操作臂的动能 k k k可描述为关节位置和速度的标量函数:
k ( Θ , Θ ˙ ) = 1 2 Θ ˙ T M ( Θ ) Θ ˙ k(\Theta,\dot{\Theta})=\frac{1}{2}\dot{\Theta}^TM(\Theta)\dot{\Theta} k(Θ,Θ˙)=21Θ˙TM(Θ)Θ˙
这里 M ( Θ ) M(\Theta) M(Θ)为状态空间方程中介绍的 n × n n\times n n×n质量矩阵。
第 i i i个连杆的势能 u i u_i ui可以表示为:
u i = − m i 0 g T 0 P C i + u r e f i u_i=-m_i{}^0g^T{}^0P_{C_i}+u_{ref_i} ui=−mi0gT0PCi+urefi
这里 0 g ^0g 0g是 3 × 1 3\times1 3×1重力矢量, 0 P C i ^0P_{C_i} 0PCi是位于第 i i i个连杆质心的矢量, u r e f i u_{ref_i} urefi是使 u i u_i ui的最小值为零的常数(即最小取到势能零点,实际上可以相对于任意一个参考零点)。操作臂的总势能为各个连杆势能之和:
u = Σ i = 1 n u i u=\Sigma^n_{i=1}u_i u=Σi=1nui
式中 0 P C i ^0P_{C_i} 0PCi是 Θ \Theta Θ的函数,因此操作臂的势能 u u u可以描述为关节位置的标量函数。
拉格朗日力学公式给出了一种从标量函数推导动力学方程的方法,该标量函数称为拉格朗日函数,即一个机械系统的动能和势能的差值:
L ( Θ , Θ ˙ ) = k ( Θ , Θ ˙ ) − u ( Θ ) \mathcal{L}(\Theta,\dot{\Theta})=k(\Theta,\dot{\Theta})-u(\Theta) L(Θ,Θ˙)=k(Θ,Θ˙)−u(Θ)
则操作臂的运动方程为:
d d t ∂ L ∂ Θ ˙ − ∂ L ∂ Θ = d d t ∂ k ∂ Θ ˙ − ∂ k ∂ Θ + ∂ u ∂ Θ = τ \frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial\mathcal{L}}{\partial\dot{\Theta}}-\frac{\partial\mathcal{L}}{\partial\Theta}=\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial k}{\partial\dot{\Theta}}-\frac{\partial k}{\partial\Theta}+\frac{\partial u}{\partial\Theta}=\tau dtd∂Θ˙∂L−∂Θ∂L=dtd∂Θ˙∂k−∂Θ∂k+∂Θ∂u=τ
笛卡尔空间中的操作臂动力学
笛卡尔状态空间方程
已知关节变量建立的状态空间方程为:
τ = M ( Θ ) Θ ¨ + V ( Θ , Θ ˙ ) + G ( Θ ) \tau=M(\Theta)\ddot{\Theta}+V(\Theta,\dot{\Theta})+G(\Theta) τ=M(Θ)Θ¨+V(Θ,Θ˙)+G(Θ)⑦
在操作臂控制当中,有时我们希望用笛卡尔变量的一般形式建立操作臂动力学方程:
F = M x ( Θ ) χ ¨ + V x ( Θ , Θ ˙ ) + G x ( Θ ) \mathcal{F}=M_x(\Theta)\ddot{\chi}+V_x(\Theta,\dot{\Theta})+G_x(\Theta) F=Mx(Θ)χ¨+Vx(Θ,Θ˙)+Gx(Θ)
式中 F \mathcal{F} F是作用于机器人末端执行器上的力-力矩矢量, χ \chi χ是表达末端执行器位置和姿态的笛卡尔矢量。
与关节空间相对应, M x ( Θ ) M_x(\Theta) Mx(Θ)是笛卡尔质量矩阵, V x ( Θ , Θ ˙ ) V_x(\Theta,\dot{\Theta}) Vx(Θ,Θ˙)是笛卡尔空间中的速度项矢量, G x ( Θ ) G_x(\Theta) Gx(Θ)是笛卡尔空间中的重力项矢量。
在第五章力域中的雅可比一节中有:
τ = J T ( Θ ) F \tau=J^T(\Theta)\mathcal{F} τ=JT(Θ)F⑧
其中 J ( Θ ) J(\Theta) J(Θ)与 F \mathcal{F} F和 χ \chi χ在同一个笛卡尔坐标系下。
接下来推导关节变量与笛卡尔变量的状态空间方程之间的关系:
对⑦两边同乘 J − T J^{-T} J−T并结合⑧得到:
J − T τ = F = J − T M ( Θ ) Θ ¨ + J − T V ( Θ , Θ ˙ ) + J − T G ( Θ ) J^{-T}\tau=\mathcal{F}=J^{-T}M(\Theta)\ddot{\Theta}+J^{-T}V(\Theta,\dot{\Theta})+J^{-T}G(\Theta) J−Tτ=F=J−TM(Θ)Θ¨+J−TV(Θ,Θ˙)+J−TG(Θ)⑨
又由雅可比矩阵的定义得:
χ ˙ = J Θ ˙ \dot{\chi}=J\dot{\Theta} χ˙=JΘ˙
两边求导:
χ ¨ = J ˙ Θ ˙ + J Θ ¨ \ddot{\chi}=\dot{J}\dot{\Theta}+J\ddot{\Theta} χ¨=J˙Θ˙+JΘ¨
Θ ¨ = J − 1 χ ¨ − J − 1 J ˙ Θ ˙ \ddot{\Theta}=J^{-1}\ddot{\chi}-J^{-1}\dot{J}\dot{\Theta} Θ¨=J−1χ¨−J−1J˙Θ˙
将上式代入⑨得:
F = J − T M ( Θ ) J − 1 χ ¨ − J − T M ( Θ ) J − 1 J ˙ Θ ˙ + J − T V ( Θ , Θ ˙ ) + J − T G ( Θ ) \mathcal{F}=J^{-T}M(\Theta)J^{-1}\ddot{\chi}-J^{-T}M(\Theta)J^{-1}\dot{J}\dot{\Theta}+J^{-T}V(\Theta,\dot{\Theta})+J^{-T}G(\Theta) F=J−TM(Θ)J−1χ¨−J−TM(Θ)J−1J˙Θ˙+J−TV(Θ,Θ˙)+J−TG(Θ)
由此得出笛卡尔空间动力学方程中各项的表达式:
M x ( Θ ) = J − T ( Θ ) M ( Θ ) J − 1 ( Θ ) M_x(\Theta)=J^{-T}(\Theta)M(\Theta)J^{-1}(\Theta) Mx(Θ)=J−T(Θ)M(Θ)J−1(Θ)
V x ( Θ , Θ ˙ ) = J − T ( Θ ) ( V ( Θ , Θ ˙ ) − M ( Θ ) J − 1 ( Θ ) J ˙ ( Θ ) Θ ˙ ) V_x(\Theta,\dot{\Theta})=J^{-T}(\Theta)(V(\Theta,\dot{\Theta})-M(\Theta)J^{-1}(\Theta)\dot{J}(\Theta)\dot{\Theta}) Vx(Θ,Θ˙)=J−T(Θ)(V(Θ,Θ˙)−M(Θ)J−1(Θ)J˙(Θ)Θ˙)
G x ( Θ ) = J − T ( Θ ) G ( Θ ) G_x(\Theta)=J^{-T}(\Theta)G(\Theta) Gx(Θ)=J−T(Θ)G(Θ)
笛卡尔位形空间中的力矩方程
用笛卡尔空间动力学方程写出等价的关节力矩:
τ = J T ( Θ ) ( M x ( Θ ) χ ¨ + V x ( Θ , Θ ˙ ) + G x ( Θ ) ) \tau=J^T(\Theta)(M_x(\Theta)\ddot{\chi}+V_x(\Theta,\dot{\Theta})+G_x(\Theta)) τ=JT(Θ)(Mx(Θ)χ¨+Vx(Θ,Θ˙)+Gx(Θ))
按与⑥类似的思路将上式改写为位形空间方程得:
τ = J T ( Θ ) M x ( Θ ) χ ¨ + B x ( Θ ) ( Θ ˙ Θ ˙ ) + C x ( Θ ) ( Θ ˙ 2 ) + G ( Θ ) \tau=J^T(\Theta)M_x(\Theta)\ddot{\chi}+B_x(\Theta)(\dot{\Theta}\dot{\Theta})+C_x(\Theta)(\dot{\Theta}^2)+G(\Theta) τ=JT(Θ)Mx(Θ)χ¨+Bx(Θ)(Θ˙Θ˙)+Cx(Θ)(Θ˙2)+G(Θ)
式中 B x ( Θ ) B_x(\Theta) Bx(Θ)是 n × n ( n − 1 ) / 2 n\times n(n-1)/2 n×n(n−1)/2阶科氏力系数矩阵, ( Θ ˙ Θ ˙ ) (\dot{\Theta}\dot{\Theta}) (Θ˙Θ˙)是 n ( n − 1 ) / 2 × 1 n(n-1)/2\times1 n(n−1)/2×1阶关节速度矢量,即:
( Θ ˙ Θ ˙ ) = ( θ ˙ 1 θ ˙ 2 , θ ˙ 1 θ ˙ 3 , ⋯ , θ ˙ n − 1 θ ˙ n ) T (\dot{\Theta}\dot{\Theta})=(\dot{\theta}_1\dot{\theta}_2,\dot{\theta}_1\dot{\theta}_3,\cdots,\dot{\theta}_{n-1}\dot{\theta}_n)^T (Θ˙Θ˙)=(θ˙1θ˙2,θ˙1θ˙3,⋯,θ˙n−1θ˙n)T
C x ( Θ ) C_x(\Theta) Cx(Θ)是 n × n n\times n n×n阶离心力系数矩阵, ( Θ ˙ 2 ) (\dot{\Theta}^2) (Θ˙2)是 n × 1 n\times1 n×1阶矢量,即:
( Θ ˙ 2 ) = ( θ ˙ 1 2 , θ ˙ 2 2 , ⋯ , θ ˙ n 2 ) T (\dot{\Theta}^2)=(\dot{\theta}^2_1,\dot{\theta}^2_2,\cdots,\dot{\theta}^2_n)T (Θ˙2)=(θ˙12,θ˙22,⋯,θ˙n2)T
G ( Θ ) G(\Theta) G(Θ)与关节空间方程中的相同,但一般情况下 B x ( Θ ) ≠ B ( Θ ) , C x ( Θ ) ≠ C ( Θ ) B_x(\Theta)\neq B(\Theta),C_x(\Theta)\neq C(\Theta) Bx(Θ)=B(Θ),Cx(Θ)=C(Θ)。
考虑非刚体影响
我们推导出的动力学方程没有包含摩擦力。
最简单的摩擦力模型为粘性摩擦力,摩擦力矩与关节运动速度成正比:
τ f r i c t i o n = v θ ˙ \tau_{friction}=v\dot{\theta} τfriction=vθ˙
式中 v v v是粘性摩擦系数。
有时应用另一个简单的摩擦力模型,就是库伦摩擦。库伦摩擦是一个常数,它的符号取决于关节速度:
τ f r i c t i o n = c s g n ( θ ˙ ) \tau_{friction}=c\mathrm{sgn}(\dot{\theta}) τfriction=csgn(θ˙)
式中 c c c是库伦摩擦系数。当 θ ˙ = 0 \dot{\theta}=0 θ˙=0, c c c一般取 1 1 1,称为静摩擦系数;当 θ ˙ ≠ 0 \dot{\theta}\neq0 θ˙=0, c < 1 c<1 c<1,称为动摩擦系数。
(书中给出的公式及其解释博主未能理解,但根据定义,库伦摩擦即静摩擦定律和滑动摩擦定律的结合)
比较合理的模型是二者兼顾:
τ f r i c t i o n = c s g n ( θ ˙ ) + v θ ˙ \tau_{friction}=c\mathrm{sgn}(\dot{\theta})+v\dot{\theta} τfriction=csgn(θ˙)+vθ˙
在许多操作臂关节中,摩擦力也与节点位置有关。主要原因是齿轮不是理想圆,齿轮的偏心会导致摩擦力随关节位置而变化,因此一个比较复杂的摩擦力模型为:
τ f r i c t i o n = f ( θ , θ ˙ ) \tau_{friction}=f(\theta,\dot{\theta}) τfriction=f(θ,θ˙)
将其附加到刚体力学模型的动力学方程中得到一个更完整的模型:
τ = M ( Θ ) Θ ¨ + V ( Θ , Θ ˙ ) + G ( Θ ) + F ( Θ , Θ ˙ ) \tau=M(\Theta)\ddot{\Theta}+V(\Theta,\dot{\Theta})+G(\Theta)+F(\Theta,\dot{\Theta}) τ=M(Θ)Θ¨+V(Θ,Θ˙)+G(Θ)+F(Θ,Θ˙)
还有其他一些影响因素,例如连杆弯曲效应引起谐振等,但是这些因素建模十分复杂,暂不讨论。
本章完