精华内容
下载资源
问答
  • 自平衡立方体机器人动力学建模.pdf
  • 多运动模式机器人动力学建模及仿真
  • 便携式自主水下机器人动力学建模方法研究.pdf
  • 研究一种3-RRRT新型高速搬运机器人动力学建模及正向动力学求解方法。以多体系统理论和Kane方法建立了该并联机器人的运动学与动力学模型,运用违约修正约束稳定法对动力学模型进行求解。并利用Matlab软件平台进行了...
  • 利用有限元理论研究柔性并联机器人动力学建模的理论和方法,分析各支链的弹性位移及其耦合关系,建立柔性并联机器人系统的运动约束条件和动力约束条件。综合考虑构件的分布质量、集中质量以及杆件的剪切变形、弯曲...
  • 面向全滚动的球形机器人动力学建模方法,北京邮电大学孙汉旭老师
  • 机器人动力学建模之牛顿欧拉法推导

    千次阅读 多人点赞 2019-07-24 21:23:31
    机器人动力学建模之牛顿欧拉法推导 最近在研究多连杆的机器人建模,发现许多使用牛顿欧拉法的建模方式,都直接采用了如下的一条公式: Fextj=MjVj˙+βj F_{extj}=M_j \dot{V_j}+\beta_j Fextj​=Mj​Vj​˙​+βj...

    机器人动力学建模之牛顿欧拉法推导

    最近在研究多连杆的机器人建模,发现许多使用牛顿欧拉法的建模方式,都直接采用了如下的一条公式:
    F e x t j = M j V j ˙ + β j (0) F_{extj}=M_j \dot{V_j}+\beta_j \tag{0} Fextj=MjVj˙+βj(0)
    其中:
    V j = [ U j Ω j ] , 这 是 连 杆 j 相 对 于 惯 性 系 的 速 度 和 角 速 度 组 成 的 六 维 向 量 ( 表 示 在 j 系 下 ) 。 V_j=\left[ \begin{matrix} U_j\\ \Omega_j \end{matrix} \right],这是连杆j相对于惯性系的速度和角速度组成的六维向量(表示在j系下)。 Vj=[UjΩj],jj
    j H j + 1 , 代 表 的 就 是 一 个 6 维 的 变 换 矩 阵 , 将 j + 1 系 下 的 力 变 换 到 j 系 下 。 {\vphantom{}}^{j}H_{j+1}, 代表的就是一个6维的变换矩阵,将{j+1}系下的力变换到{j}系下。 jHj+1,6j+1j
    F e x t j , 连 杆 受 到 的 外 力 、 外 力 矩 在 j 系 下 的 表 示 。 F_{extj},连杆受到的外力、外力矩在j系下的表示。 Fextj,j
    M j = [ m I 3 × 3 − m r j ^ m r j ^ I j ] M_j=\left[ \begin{matrix} m I_{3\times3} & -m\hat{r_j} \\ m\hat{r_j} & I_j\\ \end{matrix} \right] Mj=[mI3×3mrj^mrj^Ij]
    β j = [ m Ω j × U j + m Ω j × ( Ω j × r j ) m r j × ( Ω j × U j ) + Ω j × I j Ω j ] \beta_j=\left[ \begin{matrix} m\Omega_j \times U_j + m\Omega_j \times (\Omega_j \times r_j)\\ m r_j \times (\Omega_j \times U_j) + \Omega_j \times I_j \Omega_j \end{matrix} \right] βj=[mΩj×Uj+mΩj×(Ωj×rj)mrj×(Ωj×Uj)+Ωj×IjΩj]
    x ^ = [ 0 − x 3 x 2 x 3 0 − x 1 − x 2 x 1 0 ] , x ^ p = x × p \hat{x}=\left[ \begin{matrix} 0 &-x_3 & x_2 \\ x_3 & 0 & -x_1 \\ -x_2 & x_1 & 0 \end{matrix} \right] ,\hat{x}p=x\times p x^=0x3x2x30x1x2x10,x^p=x×p

    这个公式的形式很漂亮,用起来也比较方便。但是,如果从牛顿-欧拉方程来看,很难直接推出这个公式。看了很多文献,也没有找到推导过程,于是索性自己推导了一番。

    推导过程

    在这里插入图片描述

    1 问题描述

    图中蓝色的物体表示一个刚体,我们选择刚体上的某一点,定义了一个与刚体固连的坐标系B,即本体系。
    我们假设刚体可被视作为很多微元的组合,并且假设刚体是均匀的。
    最终,我们希望得到外力与刚体相对惯性系的速度和加速度的关系式。

    2 符号说明

    T T T:代表刚体的动能
    P b P_b Pb:代表刚体的动量,在{B}系下的表示
    Π b \Pi_b Πb:代表刚体的角动量,在{B}系下的表示
    U b U_b Ub:代表刚体相对于惯性系的速度,在{B}系下的表示
    U s U_s Us:代表刚体相对于惯性系的速度,在{S}系下的表示
    V i b V_{ib} Vib:代表刚体微元 i 相对于惯性系的速度,在{B}系下的表示
    Ω b \Omega_b Ωb:代表刚体相对于惯性系的角速度,在{B}系下的表示
    Ω s \Omega_s Ωs:代表刚体相对于惯性系的角速度,在{S}系下的表示
    p b p_b pb:代表p向量在{B}系下的表示
    p s p_s ps:代表p向量在{S}系下的表示
    r c b r_{cb} rcb:代表 r c r_c rc向量在{B}系下的表示
    r i b r_{ib} rib:代表 r i r_i ri向量在{B}系下的表示
    m m m:刚体质量
    I b I_b Ib:刚体的惯量矩阵,相对于{B}系的
    F e x t b F_{extb} Fextb:刚体所受合外力在{B}系的表示
    M e x t b M_{extb} Mextb:刚体所受合外力矩在{B}系的表示

    3 基本原理

    第一步,计算刚体的动能。

    动能是个标量,在任何坐标系下表示都相等,但这个动能必须是相对惯性系而言的。
    也就是说,速度和角速度必须是相对与惯性系的。
    

    第二步,动能对速度、角速度求导,得到动量、角动量。

    速度可以是表示在{B}系下的,也可以是{S}系下的。如果是相对于表示在{B}系下的速度求导,那么求出来的动量也是表示在{B}系下的。反之亦然。
    

    第三步,动量和角动量对时间求导,得到外力。

    这一步的求导需要格外注意。我们需要在惯性系下对时间求导,这样出来的才会是外力。牛顿定律只在惯性系下有效。
    如何在惯性系下求导,参考博文[关于机器人运动学与动力学建模的几点领悟](https://blog.csdn.net/handsome_for_kill/article/details/96473701)
    

    链接:关于机器人运动学与动力学建模的几点领悟

    4 好戏开始

    4.1 计算刚体动能

    刚体的动能会等于所有刚体微元的动能之和。

    T = ∫ 1 2 V i b 2   d m (1) T=\int \frac{1}{2}V_{ib}^2\ dm \tag{1} T=21Vib2 dm(1)

    我们知道, V i b V_{ib} Vib可以被表示为:
    V i b = U b + Ω b × r i b (2) V_{ib}=U_b+\Omega_b \times r_{ib} \tag{2} Vib=Ub+Ωb×rib(2)

    所以,我们把公式(2)代入公式(1),可以得到:
    T = ∫ 1 2 V i b 2 d m = ∫ 1 2 ( U b + Ω b × r i b ) ⋅ ( U b + Ω b × r i b ) d m = ∫ 1 2 ( U b ⋅ U b + 2 U b ⋅ Ω b × r i b + Ω b × r i b ⋅ Ω b × r i b ) d m = ∫ 1 2 U b 2 d v + ∫ U b ⋅ Ω b × r i b d m + ∫ 1 2 Ω b × r i b ⋅ Ω b × r i b d m = 1 2 m U b 2 + U b ⋅ Ω b × ∫ r i b d m + 1 2 Ω b ∫ r i b × Ω b × r i b d m = 1 2 m U b 2 + m U b ⋅ Ω b × r c b + 1 2 Ω b I b Ω b = 1 2 m U b 2 + m U b ⋅ Ω b × r c b + 1 2 I b Ω b 2 (3) \begin{aligned} T &= \int \frac{1}{2} V_{ib}^2 dm \\ &= \int \frac{1}{2} (U_b+\Omega_b \times r_{ib})\cdot(U_b+\Omega_b \times r_{ib}) dm \\ &= \int \frac{1}{2} (U_b \cdot U_b+2U_b\cdot\Omega_b \times r_{ib}+ \Omega_b \times r_{ib}\cdot\Omega_b \times r_{ib}) dm \\ &= \int \frac{1}{2} U_b^2 dv + \int U_b\cdot\Omega_b \times r_{ib} dm + \int \frac{1}{2} \Omega_b \times r_{ib}\cdot\Omega_b \times r_{ib} dm \\ &=\frac{1}{2} m U_b^2 +U_b\cdot\Omega_b \times \int r_{ib} dm+\frac{1}{2} \Omega_b \int r_{ib} \times \Omega_b \times r_{ib} dm \\ &=\frac{1}{2} m U_b^2 +mU_b\cdot\Omega_b \times r_{cb}+\frac{1}{2} \Omega_b I_b \Omega_b\\ &=\frac{1}{2} m U_b^2 +mU_b\cdot\Omega_b \times r_{cb}+\frac{1}{2} I_b \Omega_b^2\\ \end{aligned} \tag{3} T=21Vib2dm=21(Ub+Ωb×rib)(Ub+Ωb×rib)dm=21(UbUb+2UbΩb×rib+Ωb×ribΩb×rib)dm=21Ub2dv+UbΩb×ribdm+21Ωb×ribΩb×ribdm=21mUb2+UbΩb×ribdm+21Ωbrib×Ωb×ribdm=21mUb2+mUbΩb×rcb+21ΩbIbΩb=21mUb2+mUbΩb×rcb+21IbΩb2(3)

    在这部分的推导中,用到了很多公式:
    Ω b × r i b = − r i b × Ω b U b ⋅ Ω b × r i b = U b × Ω b ⋅ r i b m = ∫ 1 d m r c b = ∫ r i b d m ∫ 1 d m ∫ r i b × Ω b × r i b d m = ∫ [ r i b y 2 + r i b z 2 − r i b y r i b x − r i b z r i b x − r i b x r i b y r i b z 2 + r i b x 2 − r i b z r i b y − r i b x r i b z − r i b y r i b z r i b x 2 + r i b y 2 ] d m ⋅ Ω b = I b Ω b \begin{aligned} &\Omega_b \times r_{ib} =-r_{ib} \times \Omega_b \\ &U_b\cdot\Omega_b \times r_{ib} =U_b\times \Omega_b \cdot r_{ib} \\ &m={\int 1 dm} \\ &r_{cb}=\frac{\int r_{ib} dm}{\int 1 dm} \\ & \int r_{ib} \times \Omega_b \times r_{ib} dm=\int \left[ \begin{matrix} r_{iby}^2+r_{ibz}^2 & -r_{iby}r_{ibx}& -r_{ibz}r_{ibx}\\ -r_{ibx}r_{iby}&r_{ibz}^2+r_{ibx}^2 & -r_{ibz}r_{iby}\\ -r_{ibx}r_{ibz}& -r_{iby}r_{ibz} & r_{ibx}^2+r_{iby}^2 \\ \end{matrix} \right] dm \cdot \Omega_b = I_b\Omega_b \end{aligned} Ωb×rib=rib×ΩbUbΩb×rib=Ub×Ωbribm=1dmrcb=1dmribdmrib×Ωb×ribdm=riby2+ribz2ribxribyribxribzribyribxribz2+ribx2ribyribzribzribxribzribyribx2+riby2dmΩb=IbΩb

    最终,我们得到了动能的表达式:

    T = 1 2 m U b 2 + m U b ⋅ Ω b × r c b + 1 2 I b Ω b 2 (4) T =\frac{1}{2} m U_b^2 +mU_b\cdot\Omega_b \times r_{cb}+\frac{1}{2} I_b \Omega_b^2 \tag{4} T=21mUb2+mUbΩb×rcb+21IbΩb2(4)

    4.2 计算刚体动量

    动能对速度求导,得到动量:
    P b = ∂ T ∂ U b = m U b + m Ω b × r c b (5) P_b=\frac{\partial T}{\partial U_b}=mU_b+m\Omega_b \times r_{cb} \tag{5} Pb=UbT=mUb+mΩb×rcb(5)
    动能对角速度求导,得到角动量:
    Π b = ∂ T ∂ Ω b = m r c b × U b + I b Ω b (6) \Pi_b=\frac{\partial T}{\partial \Omega_b}=m r_{cb} \times U_b+I_b\Omega_b \tag{6} Πb=ΩbT=mrcb×Ub+IbΩb(6)

    4.3 计算与合外力的关系

    首先,我们要分析,上一步算出来的刚体动量和动量矩是什么?

    动量的定义没有歧义,但是动量矩有。

    我们知道,根据动量矩定理,只有刚体对固定点(惯性系原点S)的动量矩对时间求导等于刚体上的合外力对该固定点的合外力矩。即:

    d d t Π s = M s \frac{{\rm d}}{{\rm dt}}\Pi_s=M_s dtdΠs=Ms

    而我们上面计算的是对动点b的绝对动量矩,它和 M e x t b M_{extb} Mextb是什么关系呢?

    这篇文章中就不详细推导了,这部分推导可以参见文章:机器人动力学建模之刚体动力学基础学习

    这里直接给出结论:

    d d t Π b = M e x t b − U b × P b (7) \frac{{\rm d}}{{\rm dt}}\Pi_b=M_{extb}-U_b\times P_b \tag{7} dtdΠb=MextbUb×Pb(7)

    接下来,进入正题:
    注意,之前提到,这里需要在惯性系下,将动量和角动量对时间求导:

    1、动量对时间求导,得到合外力:
    F e x t b = d s P b d t = m d s U b d t + m d s ( Ω b × r c b ) d t = m ( U b ˙ + Ω b × U b ) + m Ω ˙ b × r c b + m Ω b × ( Ω b × r c b ) = m U b ˙ + m Ω b × U b + m Ω ˙ b × r c b + m Ω b × ( Ω b × r c b ) (8) \begin{aligned} F_{extb} =\frac{d^s P_b}{dt}&=m\frac{d^s U_b}{dt}+m\frac{d^s(\Omega_b\times r_{cb})}{dt} \\ &=m(\dot{U_b}+\Omega_b \times U_b)+m\dot\Omega_b\times r_{cb}+m\Omega_b \times (\Omega_b\times r_{cb}) \\ &=m\dot{U_b}+m\Omega_b \times U_b+m\dot\Omega_b\times r_{cb}+m\Omega_b \times (\Omega_b\times r_{cb}) \\ \end{aligned} \tag{8} Fextb=dtdsPb=mdtdsUb+mdtds(Ωb×rcb)=m(Ub˙+Ωb×Ub)+mΩ˙b×rcb+mΩb×(Ωb×rcb)=mUb˙+mΩb×Ub+mΩ˙b×rcb+mΩb×(Ωb×rcb)(8)

    2、根据公式(7),角动量对时间求导,得到合外力矩:
    M e x t b = d s Π b d t + U b × P b = m d s ( r c b × U b ) d t + d s ( I b Ω b ) d t + m U b × ( Ω b × r c b ) = m ( r c b × U ˙ b + Ω b × ( r c b × U b ) ) + I b Ω ˙ b + Ω b × ( I b Ω b ) + m U b × ( Ω b × r c b ) = m r c b × U b ˙ + I b Ω ˙ b + Ω b × ( I b Ω b ) + m [ Ω b × ( r c b × U b ) + U b × ( Ω b × r c b ) ] = m r c b × U b ˙ + I b Ω ˙ b + Ω b × ( I b Ω b ) + m r c b × ( Ω b × U b ) (9) \begin{aligned} M_{extb} =\frac{d^s \Pi_b}{dt}+U_b \times P_b&=m \frac{d^s (r_{cb}\times U_b)}{dt}+ \frac{d^s (I_b\Omega_b)}{dt}+mU_b \times (\Omega_b\times r_{cb}) \\ &=m (r_{cb}\times\dot U_b+\Omega_b\times(r_{cb}\times U_b))+ I_b\dot\Omega_b+\Omega_b\times(I_b\Omega_b)+mU_b \times (\Omega_b\times r_{cb})\\ &=m r_{cb}\times \dot{U_b}+ I_b \dot{\Omega}_b+\Omega_b\times(I_b\Omega_b)+m [\Omega_b\times(r_{cb}\times U_b)+U_b \times (\Omega_b\times r_{cb})]\\ &=m r_{cb}\times \dot{U_b}+ I_b \dot{\Omega}_b+\Omega_b\times(I_b\Omega_b)+m r_{cb}\times (\Omega_b\times U_b)\\ \end{aligned} \tag{9} Mextb=dtdsΠb+Ub×Pb=mdtds(rcb×Ub)+dtds(IbΩb)+mUb×(Ωb×rcb)=m(rcb×U˙b+Ωb×(rcb×Ub))+IbΩ˙b+Ωb×(IbΩb)+mUb×(Ωb×rcb)=mrcb×Ub˙+IbΩ˙b+Ωb×(IbΩb)+m[Ωb×(rcb×Ub)+Ub×(Ωb×rcb)]=mrcb×Ub˙+IbΩ˙b+Ωb×(IbΩb)+mrcb×(Ωb×Ub)(9)

    整理一下公式(7)和公式(8),得到:
    F e x t b = m U b ˙ + m Ω b × U b + m Ω ˙ b × r c b + m Ω b × ( Ω b × r c b ) M e x t b = m r c b × U b ˙ + I b Ω ˙ b + Ω b × ( I b Ω b ) + m r c b × ( Ω b × U b ) (10) \begin{aligned} &F_{extb}=m\dot{U_b}+m\Omega_b \times U_b+m\dot\Omega_b\times r_{cb}+m\Omega_b \times (\Omega_b\times r_{cb})\\ &M_{extb}=m r_{cb}\times \dot{U_b}+ I_b \dot{\Omega}_b+\Omega_b\times(I_b\Omega_b)+m r_{cb}\times (\Omega_b\times U_b)\\ \end{aligned} \tag{10} Fextb=mUb˙+mΩb×Ub+mΩ˙b×rcb+mΩb×(Ωb×rcb)Mextb=mrcb×Ub˙+IbΩ˙b+Ωb×(IbΩb)+mrcb×(Ωb×Ub)(10)

    整理成矩阵形式:
    [ F e x t b M e x t b ] = [ m I 3 × 3 − m r ^ c b m r ^ c b I b ] [ U b ˙ Ω ˙ b ] + [ m Ω b × U b + m Ω b × ( Ω b × r c b ) m r c b × ( Ω b × U b ) + Ω b × ( I b Ω b ) ] (11) \left[ \begin{matrix} F_{extb}\\ M_{extb} \end{matrix} \right]= \left[ \begin{matrix} mI_{3\times 3} &-m\hat{r}_{cb}\\ m \hat{r}_{cb} & I_b \end{matrix} \right] \left[ \begin{matrix} \dot{U_b}\\ \dot\Omega_b \end{matrix} \right]+ \left[ \begin{matrix} m{\Omega}_b \times U_b+m\Omega_b \times (\Omega_b\times r_{cb})\\ m r_{cb}\times (\Omega_b\times U_b)+ \Omega_b\times(I_b\Omega_b) \end{matrix} \right] \tag{11} [FextbMextb]=[mI3×3mr^cbmr^cbIb][Ub˙Ω˙b]+[mΩb×Ub+mΩb×(Ωb×rcb)mrcb×(Ωb×Ub)+Ωb×(IbΩb)](11)
    可以看到,这个式子和公式(0)是完全一摸一样的,推到这里我们就得到了答案。

    值得注意,这里的推导用到了下面的求导方法,其中 r c b r_{cb} rcb是常数向量。
    Ω b × r c b = [ Ω b x , Ω b y , Ω b z ] × [ r x , r y , r z ] = [ Ω b y r z − Ω b z r y , Ω b z r x − Ω b x r z , Ω b x r y − Ω b y r x ] \Omega_b\times r_{cb}=[\Omega_{bx},\Omega_{by},\Omega_{bz}]\times[r_x,r_y,r_z]=[\Omega_{by}r_z-\Omega_{bz}r_y,\Omega_{bz}r_x-\Omega_{bx}r_z,\Omega_{bx}r_y-\Omega_{by}r_x] Ωb×rcb=[Ωbx,Ωby,Ωbz]×[rx,ry,rz]=[ΩbyrzΩbzry,ΩbzrxΩbxrz,ΩbxryΩbyrx]
    d s ( Ω b × r c b ) d t = d s d t [ ( Ω b y r z − Ω b z r y ) ⋅ i b + ( Ω b z r x − Ω b x r z ) ⋅ j b + ( Ω b x r y − Ω b y r x ) ⋅ k b ] = d s ( Ω b y r z − Ω b z r y ) d t ⋅ i b + d s ( Ω b z r x − Ω b x r z ) d t ⋅ j b + d s ( Ω b x r y − Ω b y r x ) d t ⋅ k b         + ( Ω b y r z − Ω b z r y ) ⋅ d s i b d t + ( Ω b z r x − Ω b x r z ) ⋅ d s j b d t + ( Ω b x r y − Ω b y r x ) ⋅ d s k b d t = ( Ω ˙ b y r z − Ω ˙ b z r y ) ⋅ i b + ( Ω ˙ b z r x − Ω ˙ b x r z ) ⋅ j b + ( Ω ˙ b x r y − Ω ˙ b y r x ) ⋅ k b         + Ω b × ( Ω b y r z − Ω b z r y ) ⋅ i b + Ω b × ( Ω b z r x − Ω b x r z ) ⋅ j b + Ω b × ( Ω b x r y − Ω b y r x ) ⋅ k b = Ω ˙ b × r c b + Ω b × ( Ω b × r c b ) \begin{aligned} \frac{d^s(\Omega_b\times r_{cb})}{dt}&=\frac{d^s}{dt}\Big[ (\Omega_{by}r_z-\Omega_{bz}r_y)\cdot i_b + (\Omega_{bz}r_x-\Omega_{bx}r_z) \cdot j_b + (\Omega_{bx}r_y-\Omega_{by}r_x) \cdot k_b \Big] \\ &=\frac{d^s (\Omega_{by}r_z-\Omega_{bz}r_y)}{dt}\cdot i_b+\frac{d^s(\Omega_{bz}r_x-\Omega_{bx}r_z)}{dt} \cdot j_b+\frac{d^s(\Omega_{bx}r_y-\Omega_{by}r_x)}{dt} \cdot k_b \\ &\ \ \ \ \ \ \ +(\Omega_{by}r_z-\Omega_{bz}r_y)\cdot \frac{d^s i_b}{dt} + (\Omega_{bz}r_x-\Omega_{bx}r_z) \cdot\frac{d^s j_b}{dt} + (\Omega_{bx}r_y-\Omega_{by}r_x) \cdot \frac{d^s k_b}{dt} \\ &=( \dot{\Omega}_{by}r_z- \dot\Omega_{bz}r_y)\cdot i_b + ( \dot\Omega_{bz}r_x- \dot\Omega_{bx}r_z) \cdot j_b + ( \dot\Omega_{bx}r_y- \dot\Omega_{by}r_x) \cdot k_b \\ &\ \ \ \ \ \ \ +\Omega_{b}\times(\Omega_{by}r_z-\Omega_{bz}r_y)\cdot i_b +\Omega_{b}\times (\Omega_{bz}r_x-\Omega_{bx}r_z) \cdot j_b + \Omega_{b}\times(\Omega_{bx}r_y-\Omega_{by}r_x) \cdot k_b\\ &=\dot\Omega_b\times r_{cb}+\Omega_b \times (\Omega_b\times r_{cb}) \end{aligned} dtds(Ωb×rcb)=dtds[(ΩbyrzΩbzry)ib+(ΩbzrxΩbxrz)jb+(ΩbxryΩbyrx)kb]=dtds(ΩbyrzΩbzry)ib+dtds(ΩbzrxΩbxrz)jb+dtds(ΩbxryΩbyrx)kb       +(ΩbyrzΩbzry)dtdsib+(ΩbzrxΩbxrz)dtdsjb+(ΩbxryΩbyrx)dtdskb=(Ω˙byrzΩ˙bzry)ib+(Ω˙bzrxΩ˙bxrz)jb+(Ω˙bxryΩ˙byrx)kb       +Ωb×(ΩbyrzΩbzry)ib+Ωb×(ΩbzrxΩbxrz)jb+Ωb×(ΩbxryΩbyrx)kb=Ω˙b×rcb+Ωb×(Ωb×rcb)
    同理,我们有:
    d s ( r c b × U b ) d t = r c b × U ˙ b + Ω b × ( r c b × U b ) \frac{d^s (r_{cb}\times U_b)}{dt}=r_{cb}\times\dot U_b+\Omega_b\times(r_{cb}\times U_b) dtds(rcb×Ub)=rcb×U˙b+Ωb×(rcb×Ub)
    此外,还有:
    Ω b × ( r c b × U b ) + U b × ( Ω b × r c b ) = r c b × ( Ω b × U b ) \Omega_b\times(r_{cb}\times U_b)+U_b \times (\Omega_b\times r_{cb})=r_{cb}\times (\Omega_b\times U_b) Ωb×(rcb×Ub)+Ub×(Ωb×rcb)=rcb×(Ωb×Ub)

    另外,由于对于{b}系而言, I b I_b Ib是常数矩阵,所以
    d s ( I d Ω b ) d t = d s d t ( [ I x x I x y I x z I y x I y y I y z I z x I z y I z z ] [ Ω b x Ω b y Ω b z ] ) = d s d t ( I x x Ω b x + I x y Ω b y + I x z Ω b z I y x Ω b x + I y y Ω b y + I y z Ω b z I z x Ω b x + I z y Ω b y + I z z Ω b z ) = ( d s d t ( I x x Ω b x ⋅ i b + I x y Ω b y ⋅ j b + I x z Ω b z ⋅ k b ) d s d t ( I y x Ω b x ⋅ i b + I y y Ω b y ⋅ j b + I y z Ω b z ⋅ k b ) ) d s d t ( I z x Ω b x ⋅ i b + I z y Ω b y ⋅ j b + I z z Ω b z ⋅ k b ) ) ) = ( I x x Ω ˙ b x ⋅ i b + I x y Ω ˙ b y ⋅ j b + I x z Ω ˙ b z ⋅ k b + Ω b × ( I x x Ω b x ⋅ i b + I x y Ω b y ⋅ j b + I x z Ω b z ⋅ k b ) I y x Ω ˙ b x ⋅ i b + I y y Ω ˙ b y ⋅ j b + I y z Ω ˙ b z ⋅ k b + Ω b × ( I y x Ω b x ⋅ i b + I y y Ω b y ⋅ j b + I y z Ω b z ⋅ k b ) I z x Ω ˙ b x ⋅ i b + I z y Ω ˙ b y ⋅ j b + I z z Ω ˙ b z ⋅ k b + Ω b × ( I z x Ω b x ⋅ i b + I z y Ω b y ⋅ j b + I z z Ω b z ⋅ k b ) ) = I b Ω ˙ b + Ω b × ( I b Ω b ) \begin{aligned} \frac{d^s (I_d\Omega_b)}{dt}&=\frac{d^s}{dt} \left(\left[ \begin{matrix} I_{xx} & I_{xy} & I_{xz}\\ I_{yx} & I_{yy} & I_{yz}\\ I_{zx} & I_{zy} & I_{zz}\\ \end{matrix} \right] \left[ \begin{matrix} \Omega_{bx}\\ \Omega_{by}\\ \Omega_{bz} \end{matrix} \right]\right) \\ &=\frac{d^s}{dt} \left( \begin{matrix} I_{xx}\Omega_{bx}+I_{xy}\Omega_{by}+I_{xz}\Omega_{bz}\\ I_{yx}\Omega_{bx}+I_{yy}\Omega_{by}+I_{yz}\Omega_{bz}\\ I_{zx}\Omega_{bx}+I_{zy}\Omega_{by}+I_{zz}\Omega_{bz}\\ \end{matrix} \right)\\ &=\left( \begin{matrix} \frac{d^s}{dt}(I_{xx}\Omega_{bx}\cdot i_b+I_{xy}\Omega_{by}\cdot j_b+I_{xz}\Omega_{bz}\cdot k_b)\\ \frac{d^s}{dt}(I_{yx}\Omega_{bx}\cdot i_b+I_{yy}\Omega_{by}\cdot j_b+I_{yz}\Omega_{bz}\cdot k_b))\\ \frac{d^s}{dt}(I_{zx}\Omega_{bx}\cdot i_b+I_{zy}\Omega_{by}\cdot j_b+I_{zz}\Omega_{bz}\cdot k_b)) \end{matrix} \right)\\ &=\left( \begin{matrix} I_{xx}\dot \Omega_{bx}\cdot i_b+I_{xy}\dot \Omega_{by}\cdot j_b+I_{xz}\dot \Omega_{bz}\cdot k_b+\Omega_b\times(I_{xx}\Omega_{bx}\cdot i_b+I_{xy}\Omega_{by}\cdot j_b+I_{xz}\Omega_{bz}\cdot k_b)\\ I_{yx}\dot \Omega_{bx}\cdot i_b+I_{yy}\dot \Omega_{by}\cdot j_b+I_{yz}\dot \Omega_{bz}\cdot k_b+\Omega_b\times(I_{yx}\Omega_{bx}\cdot i_b+I_{yy}\Omega_{by}\cdot j_b+I_{yz}\Omega_{bz}\cdot k_b)\\ I_{zx}\dot \Omega_{bx}\cdot i_b+I_{zy}\dot \Omega_{by}\cdot j_b+I_{zz}\dot \Omega_{bz}\cdot k_b+\Omega_b\times(I_{zx}\Omega_{bx}\cdot i_b+I_{zy}\Omega_{by}\cdot j_b+I_{zz}\Omega_{bz}\cdot k_b) \end{matrix} \right)\\ &=I_b\dot\Omega_b+\Omega_b\times(I_b\Omega_b) \end{aligned} dtds(IdΩb)=dtdsIxxIyxIzxIxyIyyIzyIxzIyzIzzΩbxΩbyΩbz=dtdsIxxΩbx+IxyΩby+IxzΩbzIyxΩbx+IyyΩby+IyzΩbzIzxΩbx+IzyΩby+IzzΩbz=dtds(IxxΩbxib+IxyΩbyjb+IxzΩbzkb)dtds(IyxΩbxib+IyyΩbyjb+IyzΩbzkb))dtds(IzxΩbxib+IzyΩbyjb+IzzΩbzkb))=IxxΩ˙bxib+IxyΩ˙byjb+IxzΩ˙bzkb+Ωb×(IxxΩbxib+IxyΩbyjb+IxzΩbzkb)IyxΩ˙bxib+IyyΩ˙byjb+IyzΩ˙bzkb+Ωb×(IyxΩbxib+IyyΩbyjb+IyzΩbzkb)IzxΩ˙bxib+IzyΩ˙byjb+IzzΩ˙bzkb+Ωb×(IzxΩbxib+IzyΩbyjb+IzzΩbzkb)=IbΩ˙b+Ωb×(IbΩb)

    小结

    最近学习建模耽误了很多时间,其实这么老的方法,直接把论文里的公式拿来用就OK了,但是我就是控制不住自己,想要把来龙去脉搞清楚。查了很多资料都没有查到,最后只好自己动手。

    当然这里也借鉴了很多地方的资料,参考如下:

    惯性矩、惯性积、转动惯量、惯性张量
    关于机器人运动学与动力学建模的几点领悟
    牛顿—欧拉方程
    UNDERWATER GLIDERS: DYNAMICS, CONTROL AND DESIGN

    最后,希望本文能够对大家有所帮助!

    展开全文
  • 应用位形流形最小嵌入模型,对带有二阶非完整约束的欠驱动机器人动力学建模和控制进行研究。采用位形流形最小嵌入模型,简化了动力学方程,为欠驱动机器人动力学建模和控制的进一步研究奠定了基础。首先用一个齐次...
  • 柔性两轮移动机器人动力学建模及模型分析,阮晓钢,任红格,设计了柔性两轮移动机器人,针对腰部柔性的结构特点,以柔性多体系统动力学理论为基础,并根据假设模态法,将其抽象为平面柔性三
  • 机器人动力学建模,可以借鉴文中提及的工具箱,进行三维的仿真建模。
  • 两轮自平衡机器人动力学建模及其平衡控制.
  • 原理和等价树状结构的混联支路并联机器人动力学建模方法, 推导出封闭形式的混联支路并联机器人逆动力学和正动力学模 型, 动力学模型的维数等于机器人的自由度数。利用最小惯性参数和递推牛顿一欧拉法可以减少支路逆...
  • 行业分类-外包设计-基于多体系统离散时间传递矩阵法的机器人动力学建模方法.zip
  • 二、工业机器人动力学 机器人动力学描述的是关节力矩、动力学参数及关节运动的关系,用于机器人动力学建模的方法很多,如牛顿-欧拉方法、拉格朗日方法、凯恩方法、算子代数方法等。对于同一个机器人,无论采用何种...

     

    二、工业机器人动力学

            机器人动力学描述的是关节力矩、动力学参数及关节运动的关系,用于机器人动力学建模的方法很多,如牛顿-欧拉方法、拉格朗日方法、凯恩方法、算子代数方法等。对于同一个机器人,无论采用何种建模方法,最终得到的动力学模型都是等价的,可以表示为:

                                         \tiny \tau _{dyn}=D\left ( q \right )\ddot{q}+C\left ( q,\dot{q} \right )+G\left ( q \right )(2-1)                                         (2-1)

            其中\tiny D\left ( q \right )为惯性项,\tiny C\left ( q,\dot{q} \right )为科氏力及离心力项,\tiny G\left ( q \right )为重力项,每一项都是机器人惯性参数与关节运动参数的函数。机器人的10个惯性参数可表示为向量的形式:\tiny P_{iner}=\left ( I_{xx}, I_{xy}, I_{xz}, I_{yy}, I_{yz}, I_{zz}, H_{x}, H_{y}, H_{z}, m\right ),其中,参数\tiny I_{xx}\sim I_{zz}为机器人惯量矩阵\tiny I中的6个参数,\tiny H_{x}\sim H_{z}\tiny H=m\times \vec{r_{c}}=m\left ( r_{cx}, r_{cy}, r_{cz}\right )的3个分量(\tiny \vec{r_{c}}为质心向量)。上述9个量均包含(2-1)在\tiny D\left ( q \right )\tiny C\left ( q,\dot{q} \right )项内,\tiny m表示质量,包含在\tiny G\left ( q \right )项。

            基于模型的控制方案主要包括计算力矩控制、动力学前馈控制等,要想通过这些控制方法实现对轨迹的完全精确跟踪,控制方法中的动力学模型必须与机器人实际的动态特性相符。而各种典型建模方法所得到的动力学模型(2-1),只是在理想情况下的结果。实际情况中,影响机器人动力学的因素很多,如加工、装配、材料分布不均等引起的偏差;关节弹性引起的变形所带来的运动学参数偏差;关节摩擦引起的摩擦力矩;由传动方案所引起的不同关节间的运动耦合等。这些因素中,很多无法进行精确建模。为了不增加动力学模型的复杂性,理想的动力学建模方法并未完全考虑这些因素的作用,因此所得到的动力学模型(1-2)与实际的机器人动力学特性是有偏差的。动力学模型的偏差映射到控制方案中,就会引起轨迹的跟踪误差。

    2.1  动力学参数辨识

            完整的动力学参数辨识主要包括动力学建模、动力学模型的线性化(辨识模型),辨识轨迹优化、辨识算法构造、参数采集与处理、试验验证等几个方面。不同辨识方案在建模、线性化、轨迹优化及试验验证方面没有太大区别,区别主要体现在辨识算法和采集方面。

            就辨识算法而言,目前已有神经网络辨识,遗传算法辨识,最大似然估计辨识,卡尔曼滤波算法辨识,最小二乘法辨识等。

            数据采集的区别主要体现在力矩的采集上,关节运动参数一般都是通过安装在电机上的编码器测得关节转角,再对关节转角进行微分得到角速度、角加速度。而力矩的采集大致分为两类,即力传感器直接测量及通过电机电流间接测量。对于直接测量,需要在几机器人关节安装力传感器,一般选在末端或基座处,一是因为其他关节在装配后没有安装力传感器的空间,另一方面,每个关节都安装测力传感器势必会大大增加辨识成本;间接测量是通过电机电流测测量值简介计算出电机的驱动力矩值,电机电流与驱动力矩满足:

                                                                               \tiny \tau _{in}=k\cdot i_{c} 

            \tiny \tau为关节力矩,\tiny k为电机转矩常数,\tiny i_{c}电机驱动电流。不同测量方案对可辨识的参数类型及参数辨识过程有影响,而不同辨识算法仅对辨识精度有影响。

    2.2  牛顿-欧拉动力学建模

            牛顿-欧拉动力学方法基于两个基本方程,即力平衡方程以及力矩平衡方程,分别为:
                                                                                 \tiny f_{c}=ma_{c}

                                                                    \tiny n_{c}=I_{c}\cdot \alpha +\omega \times \left ( I_{c}\cdot \omega \right )

            \tiny f_{c}表示作用于机械臂质心处的合力,\tiny a_{c}表示机械臂质心的线加速度,\tiny n_{c}表示作用于机械臂质心处的合力矩,\tiny I_{c}表示相对机械臂质心表示的机械臂惯量矩阵,\tiny \alpha表示机械臂角的加速度,\tiny \omega表示机械臂角的角速度。

            牛顿-欧拉动力学建模方法包括两部分,即正向运动学递推及反向动力学递推:

    (1)正向运动学递推

            角速度递推:

            角加速度递推:

            线加速度递推:


            质心处线加速度:

            其中,\tiny _{i}^{i+1}\textrm{R}表示第\tiny i与第\tiny i+1坐标系间的姿态转换矩阵,\tiny ^{i}\textrm{p}_{i+1}表示第\tiny i与第\tiny i+1坐标原点间的距离向量,\tiny ^{i+1}\textrm{z}_{i+1}表示\tiny i+1关节的轴线方向,\tiny \dot{q}_{i+1},\ddot{q}_{i+1}均为关节变量,分别表示关节的角速度和角加速度,其他各符号含义同前,其中左上标代表参数在哪个坐标系表示,右下标表示参数所隶属的机械臂。
    (2)反向运动学递推

            \tiny i机械臂质心处的合力:

            \tiny i机械臂关节处的作用力:

            \tiny i机械臂质心处的合力矩:

            \tiny i机械臂关节处的合力矩:

            最后,将作用于关节\tiny i处的合力矩向\tiny i关节轴线方向投影,得到关节\tiny i的驱动力矩:

            这样,得到工业机器人的动力学方程,一般可写为标准形式:

            \tiny \tau为机器人的驱动力矩向量,满足:

            其中, \tiny \tau _{_{i}}表示第\tiny i关节的驱动力矩。

            \tiny D\left ( q\right )称为机械臂的质量矩阵,是一个对称阵:

                                                         \tiny D\left ( q \right )=\begin{bmatrix} D_{11} &D_{12} & .... & D_{1n}\\ D_{12} &D_{22} & .... & D_{2n}\\ \vdots&\vdots&\vdots&\vdots \\ D_{1n} &D_{2n} & .... & D_{nn} \end{bmatrix}
            对角项\tiny D_{ii}通过第\tiny i关节的角加速度\tiny \ddot{q}_{i}产生对第 \tiny i关节力矩\tiny \tau _{i}的力矩分量,非对角项\tiny D_{ik}通过第\tiny k关节的角加速度\tiny \ddot{q}_{k}产生对第\tiny i关节力矩\tiny \tau _{i}的力矩分量。

            \tiny C\left ( q,\dot{q} \right )为科氏力及离心力项,满足:

                                                                               \tiny C\left ( q,\dot{q} \right )=\left ( c_{1},c_{2},...,c_{n} \right )

            其中,

                                                                                    \tiny c_{i}=\dot{q}^{T}\cdot C_{i}\cdot \dot{q}

                                                           \tiny C_{i}=\begin{bmatrix} C_{i11} &C_{i12} & .... & C_{i1n}\\C_{i12} &C_{i22} & .... & C_{i2n}\\ \vdots&\vdots&\vdots&\vdots \\ C_{i1n} &C_{i2n} & .... & C_{inn} \end{bmatrix}

            \tiny c_{ijk}的后两个下标\tiny j,k表示此力矩分量与\tiny j,k关节的速度有关,他们的动态力相互作用在关节\tiny i处产生反作用力(力矩),标号\tiny i总表示“感受”到速度引起的反作用力(力矩)的关节编号。当\tiny j=k时,\tiny c_{ijk}与关节\tiny i“感受”到的关节\tiny k的角速度产生的离心力有关;而当\tiny j\neq k时,\tiny c_{ijk}与关节\tiny i“感受”到的关节\tiny j\tiny k的速度产生的科氏力有关。

    2.3 牛顿-欧拉动力学编程

        利用Newton-Euler方程建立动力学模型,在Matlab中利用m语言编写程序,完成动力学模型正反解验证。

    (1)建立牛顿-欧拉方程逆动力学模型。输入是六个关节的期望角度、速度和加速度,输出是六个关节力矩。

    (2)建立凯恩方程正动力学模型。输入是上一级的六个关节的力矩,输出是实际期望角度、速度和加速度。

                                                              图 1 基于牛顿-欧拉方程和凯恩方程建立的机器人动力学模型simulink仿真图

                                                                                                图2 六个关节输出力矩

                                                                            图3 六个关节实际输出角度、速度和加速度 

    展开全文
  • 刚体动力学基础学习 1 符号 rrr:刚体上某个质量微元对固定点O的位置矢径 rPr_PrP​:刚体上一点P对固定点O的位置矢径 ρ\rhoρ:刚体上某个质量微元对基点P的位置矢径 ρc\rho_cρc​:刚体质心C对P的位置矢径 mmm...

    刚体动力学基础学习

    在这里插入图片描述

    1 符号

    r r r:刚体上某个质量微元对固定点O的位置矢径

    r P r_P rP:刚体上一点P对固定点O的位置矢径

    ρ \rho ρ:刚体上某个质量微元对基点P的位置矢径

    ρ c \rho_c ρc:刚体质心C对P的位置矢径

    m m m:质量

    v k v_k vk:点k的速度

    ω \omega ω:角速度

    a k a_k ak:点k的加速度

    Q Q Q:动量

    G k G_k Gk:关于点k的绝对动量矩

    G k ′ G_k^{'} Gk:关于点k的相对动量矩

    F F F:合外力

    L k L_k Lk:对点k的合外力矩

    2 动量

    (2-1) Q = ∫ m r ˙   d m = m r ˙ c = m v c Q=\int_m \dot r \ {\rm dm}=m\dot r_c=mv_c\tag{2-1} Q=mr˙ dm=mr˙c=mvc(2-1)

    3 动量矩

    动量矩的这部分内容相对复杂,因为要分成好多种情况进行讨论。分别是:

    • 对固定点的动量矩
    • 对动点的绝对动量矩
    • 对动点的相对动量矩

    首先,要对动点和固定点做区分。动点就是在动坐标系上的点,固定点就是在固定坐标系(或者说是惯性系)下的点。

    其次,要对绝对动量矩和相对动量矩做区分。动量矩的计算公式是矢径和速度的乘积的积分,所谓绝对动量矩就是用绝对速度求解的动量矩,而相对动量矩就是用相对速度(相对动坐标系的速度)来求解的动量矩。由此可知,对固定点求动量矩,都没有绝对动量矩和相对动量矩的说法,因为对固定点都是绝对动量矩。而对于动点求动量矩,就有二者的区分。

    另外,特别值得注意的是,动点中有一个点,非常特殊,那就是质心,它是占据有非常特殊地位的一个动点。对于质心而言,绝对动量矩和相对动量矩是相等的。而对于一般的动点,这一条完全不成立。

    3.1 对固定点O的动量矩

    (3-1) G O = ∫ m r × r ˙   d m = ∫ m ( r P + ρ ) × ( r ˙ P + ρ ˙ )   d m = ∫ m ( r P + ρ ) × ( v P + ω × ρ )   d m = r P × v c m + ρ c × v P m + J P ⋅ ω \begin{aligned} G_O&=\int_m r\times\dot r \ {\rm dm}\\ &=\int_m (r_P+\rho)\times (\dot r_P + \dot \rho) \ {\rm dm}\\ &=\int_m (r_P+\rho)\times (v_P + \omega\times\rho) \ {\rm dm}\\ &=r_P\times v_c m+\rho_c\times v_P m+J_P\cdot \omega\\ \end{aligned}\tag{3-1} GO=mr×r˙ dm=m(rP+ρ)×(r˙P+ρ˙) dm=m(rP+ρ)×(vP+ω×ρ) dm=rP×vcm+ρc×vPm+JPω(3-1)

    3.2 对动点P的相对动量矩

    公式(3-1)右端的第三项,实际上就是刚体在相对P点平动坐标系运动中对P点的动量矩。
    (3-2) G P ′ = J P ⋅ ω G_P^{'}=J_P\cdot\omega\tag{3-2} GP=JPω(3-2)
    它的意思就是说,假设有一个坐标系固连在P点上,这个坐标系相对固定点O在平动,然后,计算刚体相对于P点的动量矩。因此,上式也可以写成:
    (3-3) G P ′ = ∫ m ρ × ρ ˙   d m = ∫ m ( ρ c + ρ ′ ) × ( ρ ˙ c + ρ ˙ ′ )   d m = ρ c × ρ ˙ c m + J c ⋅ ω = J P ⋅ ω \begin{aligned} G_P^{'}&=\int_m \rho\times\dot\rho \ {\rm dm}\\ &=\int_m (\rho_c+\rho^{'})\times(\dot\rho_c+\dot\rho^{'})\ {\rm dm}\\ &=\rho_c\times\dot\rho_c m+J_c\cdot \omega=J_P\cdot\omega \end{aligned}\tag{3-3} GP=mρ×ρ˙ dm=m(ρc+ρ)×(ρ˙c+ρ˙) dm=ρc×ρ˙cm+Jcω=JPω(3-3)

    3.3 对动点P的绝对动量矩

    公式(3-1)右端的第二和第三项,实际上是刚体对P点的绝对动量矩。
    (3-4) G P = ρ c × v P m + J P ⋅ ω G_P=\rho_c\times v_P m+J_P\cdot \omega\tag{3-4} GP=ρc×vPm+JPω(3-4)
    而上式又可以写作:
    (3-5) G P = ∫ m ρ × r ˙   d m = ∫ m ( ρ c + ρ ′ ) × r ˙   d m = ρ c × v c m + J c ⋅ ω \begin{aligned} G_P&=\int_m \rho\times\dot r \ {\rm dm}\\ &=\int_m (\rho_c+\rho^{'})\times\dot r\ {\rm dm}\\ &=\rho_c\times v_c m+J_c\cdot \omega \end{aligned}\tag{3-5} GP=mρ×r˙ dm=m(ρc+ρ)×r˙ dm=ρc×vcm+Jcω(3-5)

    3.4 对质心C的相对动量矩和绝对动量矩

    相对动量矩
    (3-6) G c ′ = J c ⋅ ω G_c^{'}=J_c\cdot\omega\tag{3-6} Gc=Jcω(3-6)
    刚体在绝对运动(相对于固定点O的运动)中,对刚体质心的动量矩,等于,刚体在相对运动(相对基点P的运动)中对质心的动量矩。即:
    (3-7) G c = G c ′ G_c=G_c^{'}\tag{3-7} Gc=Gc(3-7)

    3.5 固定点的动量矩与对动点的动量矩之间的关系

    对固定点的动量矩与对动点P的绝对动量矩的关系为:
    (3-8) G O = r P × v c m + G P = r P × v c m + ρ c × v c m + J c ⋅ ω G_O=r_P\times v_c m+G_P=r_P\times v_c m+\rho_c\times v_c m+J_c\cdot \omega\tag{3-8} GO=rP×vcm+GP=rP×vcm+ρc×vcm+Jcω(3-8)
    当P点就是C点时,
    (3-9) G O = r c × v c m + J c ⋅ ω = r c × v c m + G c G_O=r_c\times v_c m+J_c\cdot \omega=r_c\times v_c m+G_c\tag{3-9} GO=rc×vcm+Jcω=rc×vcm+Gc(3-9)

    4 动量定理

    (4-1) d Q d t = F = m r ¨ c = m a c \frac{{\rm d}Q}{{\rm d}t}=F=m\ddot r_c=ma_c\tag{4-1} dtdQ=F=mr¨c=mac(4-1)

    上式便是牛顿动力学方程。

    5 动量矩定理

    5.1 动量矩基本定理(无需证明的)

    刚体对固定点O的动量矩定理:
    (5-1) d d t G O = L O \frac{{\rm d}}{{\rm d}t}G_O=L_O \tag{5-1} dtdGO=LO(5-1)
    刚体对质心的动量矩定理:
    (5-2) d d t G c = L c \frac{{\rm d}}{{\rm d}t}G_c=L_c \tag{5-2} dtdGc=Lc(5-2)

    5.2 刚体对动点的绝对动量矩定理

    首先申明,刚体对动点的绝对动量矩就是前文推导的 G P G_P GP。对动点的绝对动量矩的意思就是,对平移坐标系的原点求动量矩,但是求解时用的是绝对速度。

    所以,
    (5-3) G P = ∫ m ρ × r ˙   d m = ρ c × v c m + G c G_P=\int_m \rho\times\dot r \ {\rm dm}=\rho_c\times v_c m+G_c \tag{5-3} GP=mρ×r˙ dm=ρc×vcm+Gc(5-3)
    所以求导,
    (5-4) d d t G P = ρ ˙ c × m v c + ρ c × m a c + d d t G c = ( r ˙ c − r ˙ P ) × m v c + ρ c × F + L c = − v P × m v c + L P \begin{aligned} \frac{{\rm d}}{{\rm d}t}G_P &=\dot \rho_c\times mv_c+\rho_c\times ma_c+\frac{{\rm d}}{{\rm d}t}G_c\\ &=(\dot r_c-\dot r_P)\times mv_c+\rho_c\times F+L_c\\ &=-v_P\times mv_c+L_P \end{aligned}\tag{5-4} dtdGP=ρ˙c×mvc+ρc×mac+dtdGc=(r˙cr˙P)×mvc+ρc×F+Lc=vP×mvc+LP(5-4)
    其中, L P L_P LP代表了刚体所受的外力系对动点P的主矩,它的定义为:
    (5-5) L P = ρ c × F + L c L_P=\rho_c\times F+L_c\tag{5-5} LP=ρc×F+Lc(5-5)
    最终,刚体对动点的绝对动量矩定理的表达形式为:
    (5-6) d d t G P = − v P × m v c + L P \frac{{\rm d}}{{\rm d}t}G_P=-v_P\times mv_c+L_P\tag{5-6} dtdGP=vP×mvc+LP(5-6)

    5.3 刚体对动点的相对动量矩定理

    首先申明,刚体对动点的相对动量矩就是前文推导的 G P ′ ′ G_P^{''} GP。对动点的相对动量矩的意思就是,对平移坐标系的原点求动量矩,但是求解时用的是相对速度。
    (5-7) G P ′ = ∫ m ρ × ρ ˙   d m = ρ c × ρ ˙ c m + G c G_P^{'}=\int_m \rho\times\dot\rho \ {\rm dm}=\rho_c\times\dot\rho_c m+G_c\tag{5-7} GP=mρ×ρ˙ dm=ρc×ρ˙cm+Gc(5-7)
    求导:
    (5-8) d d t G P ′ = ρ c × m ρ ¨ c + d d t G c = ρ c × m ( r ¨ c − r ¨ P ) + L c = ρ c × m a c − ρ c × m a P + L c = ρ c × F − ρ c × m a P + L c = − ρ c × m a P + L P \begin{aligned} \frac{{\rm d}}{{\rm d}t}G_P^{'} &=\rho_c\times m \ddot\rho_c +\frac{{\rm d}}{{\rm d}t} G_c\\ &=\rho_c\times m (\ddot r_c-\ddot r_P) +L_c\\ &=\rho_c\times m a_c-\rho_c\times ma_P+L_c\\ &=\rho_c\times F-\rho_c\times ma_P+L_c\\ &=-\rho_c\times ma_P+L_P\\ \end{aligned}\tag{5-8} dtdGP=ρc×mρ¨c+dtdGc=ρc×m(r¨cr¨P)+Lc=ρc×macρc×maP+Lc=ρc×Fρc×maP+Lc=ρc×maP+LP(5-8)
    又因为, G P ′ = J P ⋅ ω G_P^{'}=J_P\cdot\omega GP=JPω,其中 J P J_P JP是一个张量在某组坐标系(基)下的坐标矩阵,它的取值和基的选取非常相关。如果 J P J_P JP表示在固定坐标系下,那么由于刚体的旋转, J P J_P JP的取值时时刻刻都会发生变化。而如果 J P J_P JP表示在与刚体固连的动坐标系下,那么其取值可以为一个常量。另外,必须注意,当 J P J_P JP表示在动坐标系下时,此时公式中相应的 ω \omega ω也是表示在动坐标系下的。所以:
    (5-9) d d t G P ′ = d d t ( J P ⋅ ω ) = J P ⋅ ω ˙ + ω × J P ⋅ ω \begin{aligned} \frac{{\rm d}}{{\rm d}t}G_P^{'}=\frac{{\rm d}}{{\rm d}t}(J_P\cdot\omega)={J_P}\cdot\dot\omega+\omega\times J_P\cdot\omega \end{aligned}\tag{5-9} dtdGP=dtd(JPω)=JPω˙+ω×JPω(5-9)

    所以,若把公式(5-9)带入公式(5-8),可得:
    (5-10) J P ⋅ ω ˙ + ω × J P ⋅ ω + ρ c × m a P = L P {J_P}\cdot\dot\omega+\omega\times J_P\cdot\omega+\rho_c\times ma_P=L_P \tag{5-10} JPω˙+ω×JPω+ρc×maP=LP(5-10)
    若P点就是质心,则 ρ c = 0 \rho_c=0 ρc=0,由公式(5-10),可得:
    (5-11) J c ⋅ ω ˙ + ω × J c ⋅ ω = L c {J_c}\cdot\dot\omega+\omega\times J_c\cdot\omega=L_c \tag{5-11} Jcω˙+ω×Jcω=Lc(5-11)
    上式便是鼎鼎大名的欧拉动力学方程了。

    最终,刚体对动点的相对动量矩定理的表达形式为:
    (5-12) d d t G P ′ = − ρ c × m a P + L P \frac{{\rm d}}{{\rm d}t}G_P^{'}=-\rho_c\times ma_P+L_P\tag{5-12} dtdGP=ρc×maP+LP(5-12)

    6 牛顿—欧拉公式

    (6-1) F = m a c L c = J c ⋅ ω ˙ + ω × J c ⋅ ω \begin{aligned} &F=ma_c\\ &L_c={J_c}\cdot\dot\omega+\omega\times J_c\cdot\omega \end{aligned}\tag{6-1} F=macLc=Jcω˙+ω×Jcω(6-1)

    以上,是大名鼎鼎的牛顿欧拉公式,看着形式很简单,但在应用时,有一些细节务必注意。

    • a c a_c ac是质心处的绝对加速度
    • J c J_c Jc是质心处定义的在动坐标系下表示的惯性张量
    • ω \omega ω是刚体在动坐标系下表示的角速度,通常也就是体坐标系
    • L c L_c Lc是对质心的合外力矩

    7 利用牛顿—欧拉公式推导常见的多连杆机器人动力学方程

    假设所有量都是表示在动坐标系下的,也就是刚体的随体坐标系{b}系,另外{b}系的原点设置在b点,b点与刚体质心不重合。此时的牛顿欧拉方程为:
    (7-1) F b = m a c b L c b = J c ⋅ ω ˙ b + ω b × J c ⋅ ω b \begin{aligned} &F^b=ma_c^b\\ &L_c^b={J_c}\cdot\dot\omega^b+\omega^b\times J_c\cdot\omega^b \end{aligned}\tag{7-1} Fb=macbLcb=Jc