CN112001087A - 一种旋转关节型工业机器人非线性动力学建模分析方法 - Google Patents

一种旋转关节型工业机器人非线性动力学建模分析方法 Download PDF

Info

Publication number
CN112001087A
CN112001087A CN202010875365.4A CN202010875365A CN112001087A CN 112001087 A CN112001087 A CN 112001087A CN 202010875365 A CN202010875365 A CN 202010875365A CN 112001087 A CN112001087 A CN 112001087A
Authority
CN
China
Prior art keywords
point
industrial robot
formula
matrix
hinge
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010875365.4A
Other languages
English (en)
Other versions
CN112001087B (zh
Inventor
李波
田威
崔光裕
廖文和
李宇飞
张霖
胡俊山
石瀚斌
梁爽
王俊
张苇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Changjiang Industrial Technology Research Institute Co ltd
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN202010875365.4A priority Critical patent/CN112001087B/zh
Publication of CN112001087A publication Critical patent/CN112001087A/zh
Application granted granted Critical
Publication of CN112001087B publication Critical patent/CN112001087B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Manipulator (AREA)

Abstract

本发明公开了一种旋转关节型工业机器人非线性动力学建模分析方法，包括如下步骤：将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的系统动力学模型；选取工业机器人各关节转角值作为系统广义坐标，确定积分变量x(ti)及其初始条件x(t0)；计算工业机器人体元件和铰元件的传递矩阵和传递方程；根据元件传递矩阵和机器人系统拓扑结构，计算总传递矩阵，结合边界条件，求解总传递方程，获得边界点状态矢量中的未知状态变量；根据动力学模型中各连接点处的加速度、角加速度，对动力学模型进行运动学分析，获得广义坐标的二阶导数；结合数值积分方法，获得ti+1时刻的积分变量x(ti+1)，当到达期望的计算时间则结束，否则返回重新确定初始条件并再次分析计算。

Description

一种旋转关节型工业机器人非线性动力学建模分析方法
技术领域
本发明属于工业机器人动力学建模技术领域，具体涉及一种旋转关节型工业机器人非线 性动力学建模分析方法。
背景技术
工业机器人在机械加工、汽车制造、航空航天制造等自动化生产领域的广泛应用，对其 运动控制精度和动态特性等核心性能提出了更高要求。传统的基于运动学的控制与设计研究 已经无法满足应用需求，必须从动力学层面入手，通过对机器人加工系统进行动力学研究， 获取机器人运动过程中的动力响应，作为机器人高精度控制和优化设计的基础。机器人的动 力学研究对于提高机器人运动精度与动态特性至关重要，一方面必须建立能够准确描述机器 人系统动态特性的动力学方程，另一方面在保证动力学建模精度的前提下实现动力学快速计 算，满足实时控制要求。
严格来说工业机器人动力学系统是一个非线性多体系统，其动力学方程相当复杂。动力 学研究方法主要有牛顿-欧拉法、拉格朗日法、约丹速度变分法、Kane法、Wittenburg法、 Schiehlen法、Shabana法等。已公开发明专利“一种工业机器人整体动力学建模及动力学参数 辨识方法，CN 110539302 A”采用拉格朗日方法建立工业机器人整体动力学模型，并进行了动 力学参数辨识。现有技术存在的问题是：现有动力学建模方法在计算工业机器人动力学响应 时，需要建立系统的总体动力学方程，虽然能够较为准确地描述工业机器人的实际动态特性， 但其系统矩阵阶次随着自由度的增加不断提高，导致计算速度下降、难以满足机器人运动过 程中的实时计算与控制需求。
发明内容
本发明所要解决的技术问题是针对上述现有技术的不足，提供一种旋转关节型工业机器 人非线性动力学建模分析方法。
为实现上述技术目的，本发明采取的技术方案为：
一种旋转关节型工业机器人非线性动力学建模分析方法，其中：包括如下步骤：
步骤S1：将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体 系统动力学模型，并建立拓扑结构图；
步骤S2：选取工业机器人各关节转角值作为系统广义坐标，确定积分变量x(ti)及其初始 参数；
步骤S3：对于ti时刻，根据积分变量进行对铰元件和体元件进行运动学分析，得到各体 元件的位置、姿态、速度、角速度变量；
步骤S4：计算工业机器人体元件和铰元件的传递矩阵和传递方程；
步骤S5：根据元件传递矩阵和机器人系统拓扑结构，计算总传递矩阵，结合边界条件， 求解总传递方程，获得边界点状态矢量中的未知状态变量；
步骤S6：根据动力学模型中各连接点处的加速度、角加速度，对动力学模型进行运动学 分析，获得广义坐标的二阶导数，结合数值积分方法，获得ti+1时刻的积分变量x(ti+1)，当到 达期望的计算时间，则结束，否则返回步骤S3继续计算。
为优化上述技术方案，采取的具体措施还包括：
进一步地，步骤S1具体为：
S11：根据旋转型工业机器人的几何结构、装配关系，以及各个部件间的连接方式，将旋 转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体系统动力学模型；
S12：将工业机器人的连杆和连接件等价于体元件，将工业机器人用于连接的关节等价于 铰元件，其中体元件包含质量，铰元件不包含质量；
S13：将工业机器人的基座边界编号为0，铰元件编号为奇数，体元件编号为偶数。
进一步地，步骤S2具体为：选取工业机器人各关节转角值作为系统广义坐标，即 q＝[θ1 θ3 θ5 θ7 θ9 θ11]T，确定积分变量
Figure BDA0002652520310000021
和每个铰元件的初始转动角度和角速度。
进一步地，步骤S3具体为：
S31：以工业机器人的基座为原点建立oxyz惯性坐标系，以体元件的力的输入点I为中心 建立IxIyIzI连体坐标系，建立各部件连接点的状态矢量z为
Figure BDA0002652520310000022
其中，
Figure BDA0002652520310000023
为连接点的绝对加速度在惯性坐标系oxyz的坐标列阵，
Figure BDA0002652520310000024
为连接点处对应坐标系的绝对角加速度在惯性系oxyz的坐标列阵， m＝[mx my mz]T和q＝[qx qy qz]T分别为连接点处内力矩、内力在惯性系oxyz的坐标列阵；
S32：对体元件进行运动学分析
由体元件的输入点I和输出点O的绝对角速度、角加速度相等得：
Figure BDA0002652520310000031
式中，
Figure BDA0002652520310000032
表示输入点I的绝对角速度和绝对角加速度，
Figure BDA0002652520310000033
表示输出点O的 绝对角速度和绝对角加速度；
输出点O相对于惯性坐标系原点的位置矢量在惯性坐标系中的投影rO表示为
rO＝rI+rIO＝rI+AIlIO (3)
ΩI＝AIωI
式中，rI为输入点I相对于惯性坐标系原点的位置矢量，rIO和lIO分别表示输出点相对于 输入点的位置矢量在惯性坐标系和连体坐标系IxIyIzI的投影，AI为输入点连体系IxIyIzI相对于 惯性系oxyz的方向余弦矩阵，ΩI、ωI分别为输入点连体坐标系在惯性坐标系中的角速度矢量 在惯性坐标系oxyz和输入端连体坐标系IxIyIzI中的投影；
计算体元件的速度
Figure BDA0002652520310000034
S33：对铰元件进行动力学分析，则铰元件输入点与输出点状态矢量为
Figure BDA0002652520310000035
式中，Aj,O为铰元件输出点连体系相对于惯性坐标系的方向余弦矩阵，
Figure BDA0002652520310000036
表 示铰元件外接刚体相对于内接刚体的角速度。
进一步地，步骤S4具体为：对式(3)关于时间求导可得
Figure BDA0002652520310000037
Figure BDA0002652520310000038
Figure BDA0002652520310000039
Figure BDA00026525203100000310
其中
Figure BDA00026525203100000311
表示任一向量a对应的反对称矩阵；
将式(9)代入式(6)可得
Figure BDA0002652520310000041
同理，对于质心C，有
Figure BDA0002652520310000042
根据刚体质心运动定律，可得体元件随输入点平动的动力学方程为
Figure BDA0002652520310000043
式中，m为体元件质量，qI,qO分别为输入点I、输出点O的内力，fC为作用在质心C处的外力；
根据体元件相对于任意动点的相对动量矩定理，选择输入点I为动点，并投影至惯性坐 标系，可得
Figure BDA0002652520310000044
式中，JI为体元件相对于输入点的转动惯量在输入点连体坐标系IxIyIzI中的矩阵表示，mC表示作用在体元件质心处的外力矩；
将式(12)代入式(13)，并整理式(10)、(4)、(13)、(12)可得
Figure BDA0002652520310000045
将式(14)写成矩阵的形式
zO＝UzI (15)
可得体元件的传递矩阵为
Figure BDA0002652520310000046
式中，I3为3阶单位矩阵，O3为3×3全零矩阵，O3×1为3×1全零矩阵，O1×3为1×3全零矩 阵；
得到各体元件的传递方程为
Figure BDA0002652520310000051
式中，Ui(i＝2,4,6,8,10,12)为体元件的传递矩阵，z2,1,z2,3,L,z12,11为工业机器人系统内各 连接点状态矢量，z12,0为工业机器人系统自由端边界点状态矢量；
S42：计算铰元件的传递矩阵和传递方程
由铰元件的输入点与输出点状态矢量式(5)第四项等式对时间求一阶导数可得
Figure BDA0002652520310000052
式(18)两边同时乘以
Figure BDA0002652520310000053
得
Figure BDA0002652520310000054
式中，
Figure BDA0002652520310000055
Figure BDA0002652520310000056
为方向余弦矩阵Aj,O的转置。
铰元件绕z轴方向的内力矩
Figure BDA0002652520310000057
可得
Figure BDA0002652520310000058
式中，H2＝[0 0 1]。
当外接体元件j+1后再接一个铰元件，可得
Figure BDA0002652520310000059
由于外接体元件j+1的传递方程和传递矩阵分别为
zj+1,O＝Uj+1zj+1,I (22)
Figure BDA00026525203100000510
式中，u3,1、u3,2、u3,3、u3,4、u3,5为铰元件j外接体元件j+1传递矩阵的分块矩阵中的元素。
则有
Figure BDA0002652520310000061
等式两边同时左乘
Figure BDA0002652520310000062
并结合式(21)，整理可得
Figure BDA0002652520310000063
由于第j个铰元件输出点状态矢量与第j+1个体元件输入点状态矢量相等，即
Figure BDA0002652520310000064
将式(26)代入式(25)，得
Figure BDA0002652520310000065
联立式(19)与式(27)，可得铰元件j的传递矩阵为
Figure BDA0002652520310000066
Figure BDA0002652520310000067
Figure RE-GDA0002705925310000068
则各铰元件的传递方程为
Figure BDA0002652520310000069
式中，Uj(j＝1,3,5,7,9,11)为铰元件的传递矩阵，z1,0为工业机器人系统固定端边界点状态 矢量。
进一步地，步骤S5具体为：
根据各体元件传递方程(17)和各铰元件传递方程(29)，得到机器人系统总传递方程
z12,0＝U1-12z1,0 (30)
式中，U1-12＝U12…U2U1为系统总传递矩阵；
工业机器人固定端边界点和自由端边界点的边界条件为
Figure BDA0002652520310000071
将边界条件(31)代入式(30)，得到分块矩阵形式为
Figure BDA0002652520310000072
Figure BDA0002652520310000073
优化式(33)，计算边界点状态矢量中的未知状态变量
Figure BDA0002652520310000074
根据计算得到的未知状态变量获得系统固定端状态矢量，然后依次利用各元件的传递方 程(17)和(29)，获得系统中各连接点的状态矢量。
进一步地，所述步骤S6中铰元件θj(j＝1,3,5,7,9,11)的广义坐标的二阶导数的计算过程为：
在式(18)两边同时左乘
Figure BDA0002652520310000075
得到
Figure BDA0002652520310000076
本发明的有益效果：
本发明提供了一种旋转关节型工业机器人非线性动力学建模分析方法，无需系统总体动 力学方程、系统矩阵阶次不随系统自由度的增加而增加，可大幅度提高机器人空间大运动非 线性动力学建模与计算效率。
附图说明
图1是本发明的旋转关节型工业机器人多体系统动力学模型示意图；
图2是本发明的旋转关节型工业机器人动力学模型拓扑结构示意图；
图3是本发明的空间运动刚体及其坐标系示意图；
图4是本发明的空间运动柱铰及其坐标系示意图；
图5是本发明的旋转关节型工业机器人多体系统动力学数值计算流程图。
具体实施方式
以下结合附图对本发明的实施例作进一步详细描述，本发明采用KUKA KR500-2830MT 型工业机器人多体系统动力学模型。
如图1所示，本发明为一种旋转关节型工业机器人非线性动力学建模分析方法，其中： 包括如下步骤：
步骤S1：将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体 系统动力学模型，并建立拓扑结构图；
S11：根据旋转型工业机器人的几何结构、装配关系，以及各个部件间的连接方式，将旋 转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体系统动力学模型；
S12：将工业机器人的连杆和连接件等价于体元件，将工业机器人用于连接的关节等价于 铰元件，其中体元件包含质量，铰元件不包含质量；
S13：将工业机器人的基座边界编号为0，铰元件编号为奇数，体元件编号为偶数。
机器人地基编号为0；6个关节依次编号为1、3、5、7、9、11；6个连杆依次编号为2、4、6、8、10、12。需要指出的是，KUKA机器人的后三个电机在第3连杆上，因此体元件6 的质量包括第3连杆及后三个电机伺服机构的质量。对于其他型号的机器人，只需根据其电 机的具体布局方式，选择对应的铰接形式即可。
将机器人各连杆视为空间运动刚体，各关节视为空间光滑柱铰，电机产生的柱铰轴方向 的控制力矩分别作用在相邻的两个连杆上，且锁止机构工作时的弹性阻尼作用等效为与柱铰 轴方向并联的扭簧和旋转阻尼器。将机器人加工时外界对末端执行器的力处理为作用在机器 人系统上的外力。
如图2所示，KUKA KR500-2830MT型工业机器人多体动力学模型即为：在电机驱动力 和外界负载加工力作用下的由6个空间运动刚体和6个空间柱铰组成的多刚体系统。根据 KUKA KR500-2830MT型工业机器人多体系统动力学模型，建立其动力学模型拓扑结构，其 中封闭曲线代表体元件；箭头代表铰元件，箭头方向代表机器人系统连接点状态矢量的传递 方向；数字为元件编号；边界统一编号为0。
步骤S2：选取工业机器人各关节转角值作为系统广义坐标，确定积分变量x(ti)及其初始 参数；
步骤S2具体为：选取工业机器人各关节转角值作为系统广义坐标，即 q＝[θ1 θ3θ5 θ7 θ9 θ11]T，确定积分变量
Figure BDA0002652520310000091
和每个铰元件的初始转动角度和角速 度。
步骤S3：于ti时刻，根据积分变量进行对铰元件和体元件进行运动学分析，得到各体元 件的位置、姿态、速度、角速度变量。
步骤S3具体为：
S31：以工业机器人的基座为原点建立oxyz惯性坐标系，以体元件的力的输入点I为中 心建立IxIyIzI连体坐标系，建立各部件连接点的状态矢量z为
Figure BDA0002652520310000092
其中，
Figure BDA0002652520310000093
为连接点的绝对加速度在惯性坐标系oxyz的坐标列阵，
Figure BDA0002652520310000094
为连接点处对应坐标系的绝对角加速度在惯性系oxyz的坐标列阵， m＝[mx my mz]T和q＝[qx qy qz]T分别为连接点处内力矩、内力在惯性系oxyz的坐标列阵；
S32：对体元件进行运动学分析
由体元件的输入点I和输出点O的绝对角速度、角加速度相等得：
Figure BDA0002652520310000095
式中，
Figure BDA0002652520310000096
表示输入点I的绝对角速度和绝对角加速度，
Figure BDA0002652520310000097
表示输出点O的 绝对角速度和绝对角加速度；
输出点O相对于惯性坐标系原点的位置矢量在惯性坐标系中的投影rO表示为
rO＝rI+rIO＝rI+AIlIO (3)
ΩI＝AIωI
式中，rI为输入点I相对于惯性坐标系原点的位置矢量，rIO和lIO分别表示输出点相对于 输入点的位置矢量在惯性坐标系和连体坐标系IxIyIzI的投影，AI为输入点连体系IxIyIzI相对于 惯性系oxyz的方向余弦矩阵，ΩI、ωI分别为输入端连体坐标系在惯性坐标系中的角速度矢量 在惯性坐标系oxyz和输入端连体坐标系IxIyIzI中的投影；
计算体元件的速度
Figure BDA0002652520310000101
S33：对铰元件进行动力学分析，则铰元件输入点与输出点状态矢量为
Figure BDA0002652520310000102
式中，Aj,O为铰元件输出点连体系相对于惯性坐标系的方向余弦矩阵，
Figure BDA0002652520310000103
表 示铰元件外接刚体相对于内接刚体的角速度。
步骤S4：计算工业机器人体元件和铰元件的传递矩阵和传递方程；
步骤S4具体为：
S41：计算体元件的传递矩阵和传递方程
对式(3)关于时间求导可得
Figure BDA0002652520310000104
Figure BDA0002652520310000105
Figure BDA0002652520310000106
Figure BDA0002652520310000107
其中
Figure BDA0002652520310000108
表示任一向量a对应的反对称矩阵；
将式(9)代入式(6)可得
Figure BDA0002652520310000109
同理，对于质心C，有
Figure BDA00026525203100001010
根据刚体质心运动定律，可得体元件随输入点平动的动力学方程为
Figure BDA00026525203100001011
式中，m为体元件质量，qI,qO分别为输入点I、输出点O的内力，fC为作用在质心C处的外力；
根据体元件相对于任意动点的相对动量矩定理，选择输入点I为动点，并投影至惯性坐 标系，可得
Figure BDA0002652520310000111
式中，JI为体元件相对于输入点的转动惯量在输入点连体坐标系IxIyIzI中的矩阵表示，mC表示作用在体元件质心处的外力矩；
将式(12)代入式(13)，并整理式(10)、(4)、(13)、(12)可得
Figure BDA0002652520310000112
将式(14)写成矩阵的形式
zO＝UzI (15)
可得体元件的传递矩阵为
Figure BDA0002652520310000113
式中，I3为3阶单位矩阵，O3为3×3全零矩阵，O3×1为3×1全零矩阵，O1×3为1×3全零矩 阵；
得到各体元件的传递方程为
Figure BDA0002652520310000114
式中，Ui(i＝2,4,6,8,10,12)为体元件的传递矩阵，z2,1,z2,3,L,z12,11为工业机器人系统内各 连接点状态矢量，z12,0为工业机器人系统自由端边界点状态矢量；
S42：计算铰元件的传递矩阵和传递方程
由铰元件的输入点与输出点状态矢量式(5)第四项等式对时间求一阶导数可得
Figure BDA0002652520310000115
式(18)两边同时乘以
Figure BDA0002652520310000116
得
Figure BDA0002652520310000121
式中，
Figure BDA0002652520310000122
Figure BDA0002652520310000123
为方向余弦矩阵Aj,O的转置；
铰元件绕z轴方向的内力矩
Figure BDA0002652520310000124
可得
Figure BDA0002652520310000125
式中，H2＝[0 0 1]。
当外接体元件j+1后再接一个铰元件，可得
Figure BDA0002652520310000126
由于外接体元件j+1的传递方程和传递矩阵分别为
zj+1,O＝Uj+1zj+1,I (22)
Figure BDA0002652520310000127
式中，u3,1、u3,2、u3,3、u3,4、u3,5为铰元件j外接体元件j+1传递矩阵的分块矩阵中的元素；
则有
Figure BDA0002652520310000128
等式两边同时左乘
Figure BDA0002652520310000129
并结合式(21)，整理可得
Figure BDA00026525203100001210
由于第j个铰元件输出点状态矢量与第j+1个体元件输入点状态矢量相等，即
Figure BDA00026525203100001211
将式(26)代入式(25)，得
Figure BDA00026525203100001212
联立式(19)与式(27)，可得铰元件j的传递矩阵为
Figure BDA0002652520310000131
Figure BDA0002652520310000132
Figure RE-GDA0002705925310000133
则各铰元件的传递方程为
Figure BDA0002652520310000134
式中，Uj(j＝1,3,5,7,9,11)为铰元件的传递矩阵，z1,0为工业机器人系统固定端边界点状态 矢量。
步骤S5：根据元件传递矩阵和机器人系统拓扑结构，计算总传递矩阵，结合边界条件， 求解总传递方程，获得边界点状态矢量中的未知状态变量；
步骤S5具体为:
根据各体元件传递方程(17)和各铰元件传递方程(29)，得到机器人系统总传递方程
z12,0＝U1-12z1,0 (30)
式中，U1-12＝U12…U2U1为系统总传递矩阵；
工业机器人固定端边界点和自由端边界点的边界条件为
Figure BDA0002652520310000135
将边界条件(31)代入式(30)，得到分块矩阵形式为
Figure BDA0002652520310000141
Figure BDA0002652520310000142
优化式(33)，计算边界点状态矢量中的未知状态变量
Figure BDA0002652520310000143
根据计算得到的未知状态变量获得系统固定端状态矢量，然后依次利用各元件的传递方 程(17)和(29)，获得系统中各连接点的状态矢量。
步骤S5：根据动力学模型中各连接点处的加速度、角加速度，对动力学模型进行运动学 分析，获得广义坐标的二阶导数；
所述步骤S6中铰元件θj(j＝1,3,5,7,9,11)的广义坐标的二阶导数的计算过程为：在式(18) 两边同时左乘
Figure BDA0002652520310000144
得到
Figure BDA0002652520310000145
以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于 本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术 人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (7)

1.一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于，包括如下步骤：
步骤S1：将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体系统动力学模型，并建立拓扑结构图；
步骤S2：选取工业机器人各铰元件作为系统广义坐标，确定积分变量x(ti)及其初始条件；
步骤S3：对于ti时刻，根据积分变量进行对铰元件和体元件进行运动学分析，得到各体元件的位置、姿态、速度、角速度变量；
步骤S4：计算工业机器人体元件和铰元件的传递矩阵和传递方程；
步骤S5：根据元件传递矩阵和机器人系统拓扑结构，计算总传递矩阵，结合边界条件，求解总传递方程，获得边界点状态矢量中的未知状态变量；
步骤S6：根据动力学模型中各连接点处的加速度、角加速度，对动力学模型进行运动学分析，获得广义坐标的二阶导数,结合数值积分方法，获得ti+1时刻的积分变量x(ti+1)，当到达期望的计算时间，则结束，否则返回步骤S3继续计算。
2.根据权利要求1所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S1具体为：
S11：根据旋转型工业机器人的几何结构、装配关系，以及各个部件间的连接方式，将旋转型工业机器人三维实体模型转化为由体元件和铰元件组成的机器人多体系统动力学模型；
S12：将工业机器人的连杆和连接件等价于体元件，将工业机器人用于连接的关节等价于铰元件，其中体元件包含质量，铰元件不包含质量；
S13：将工业机器人的基座边界编号为0，铰元件编号为奇数，体元件编号为偶数。
3.根据权利要求2所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S2具体为：选取工业机器人各关节转角值作为系统广义坐标，即q＝[θ1θ3 θ5 θ7 θ9 θ11]T，确定积分变量
Figure FDA0002652520300000011
和每个铰元件的初始转动角度和角速度。
4.根据权利要求3所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S3具体为：
S31：以工业机器人的基座为原点建立oxyz惯性坐标系，以体元件的力的输入点I为中心建立IxIyIzI连体坐标系，建立各部件连接点的状态矢量z为
Figure FDA0002652520300000012
其中，
Figure FDA0002652520300000021
为连接点的绝对加速度在惯性坐标系oxyz的坐标列阵，
Figure FDA0002652520300000022
为连接点处对应坐标系的绝对角加速度在惯性系oxyz的坐标列阵，m＝[mx my mz]T和q＝[qx qy qz]T分别为连接点处内力矩、内力在惯性系oxyz的坐标列阵；
S32：对体元件进行运动学分析
由体元件的输入点I和输出点O的绝对角速度、角加速度相等得：
Figure FDA0002652520300000023
式中，ΩI、
Figure FDA0002652520300000024
表示输入点I的绝对角速度和绝对角加速度，ΩO、
Figure FDA0002652520300000025
表示输出点O的绝对角速度和绝对角加速度；
输出点O相对于惯性坐标系原点的位置矢量在惯性坐标系中的投影rO表示为
rO＝rI+rIO＝rI+AIlIO (3)
ΩI＝AIωI
式中，rI为输入点I相对于惯性坐标系原点的位置矢量，rIO和lIO分别表示输出点相对于输入点的位置矢量在惯性坐标系和连体坐标系IxIyIzI的投影，AI为输入点连体系IxIyIzI相对于惯性系oxyz的方向余弦矩阵，ΩI、ωI分别为输入端连体坐标系在惯性坐标系中的角速度矢量在惯性坐标系oxyz和输入端连体坐标系IxIyIzI中的投影；
计算体元件的速度
Figure FDA0002652520300000026
式中，
Figure FDA0002652520300000027
表示输出点O的速度，
Figure FDA0002652520300000028
表示输入点I的速度，
Figure FDA0002652520300000029
表示AI对时间的一阶导数。
S33：对铰元件进行动力学分析，则铰元件输入点与输出点状态矢量为
Figure FDA00026525203000000210
式中，
Figure FDA00026525203000000211
mj,O、qj,O、Ωj,O分别为第j个铰元件输出点的绝对加速度、内力矩、内力和绝对角速度的坐标列阵；
Figure FDA00026525203000000212
mj,I、qj,I、Ωj,I分别为第j个铰元件输入点的绝对加速度、内力矩、内力和绝对角速度的坐标列阵；Aj,O为铰元件输出点连体系相对于惯性坐标系的方向余弦矩阵，
Figure FDA0002652520300000031
表示铰元件外接刚体相对于内接刚体的角速度，θj为旋转角度。
5.根据权利要求4所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S4具体为：
S41：计算体元件的传递矩阵和传递方程
对式(3)关于时间求导可得
Figure FDA0002652520300000032
Figure FDA0002652520300000033
Figure FDA0002652520300000034
Figure FDA0002652520300000035
其中
Figure FDA0002652520300000036
表示任一向量a对应的反对称矩阵；
将式(9)代入式(6)可得
Figure FDA0002652520300000037
同理，对于质心C，有
Figure FDA0002652520300000038
根据刚体质心运动定律，可得体元件随输入点平动的动力学方程为
Figure FDA0002652520300000039
式中，m为体元件质量，qI,qO分别为输入点I、输出点O的内力，fC为作用在质心C处的外力；
根据体元件相对于任意动点的相对动量矩定理，选择输入点I为动点，并投影至惯性坐标系，可得
Figure FDA00026525203000000310
式中，JI为体元件相对于输入点的转动惯量在输入点连体坐标系IxIyIzI中的矩阵表示，mC表示作用在体元件质心处的外力矩；
将式(12)代入式(13)，并整理式(10)、(4)、(13)、(12)可得
Figure FDA0002652520300000041
将式(14)写成矩阵的形式
zO＝UzI (15)
可得体元件的传递矩阵为
Figure FDA0002652520300000042
式中，I3为3阶单位矩阵，O3为3×3全零矩阵，O3×1为3×1全零矩阵，O1×3为1×3全零矩阵；
得到各体元件的传递方程为
Figure FDA0002652520300000043
式中，Ui(i＝2,4,6,8,10,12)为体元件的传递矩阵，z2,1,z2,3,L,z12,11为工业机器人系统内各连接点状态矢量，z12,0为工业机器人系统自由端边界点状态矢量；
S42：计算铰元件的传递矩阵和传递方程
由铰元件的输入点与输出点状态矢量式(5)中第四项等式对时间求一阶导数可得
Figure FDA0002652520300000044
式(18)两边同时乘以
Figure FDA0002652520300000045
得
Figure FDA0002652520300000046
式中，
Figure FDA0002652520300000047
Figure FDA0002652520300000048
为方向余弦矩阵Aj,O的转置。
铰元件绕z轴方向的内力矩
Figure FDA0002652520300000049
可得
Figure FDA0002652520300000051
式中，H2＝[0 0 1]。
当外接体元件j+1后再接一个铰元件，可得
Figure FDA0002652520300000052
由于外接体元件j+1的传递方程和传递矩阵分别为
zj+1,O＝Uj+1zj+1,I (22)
Figure FDA0002652520300000053
式中，u3,1、u3,2、u3,3、u3,4、u3,5为铰元件j外接体元件j+1传递矩阵的分块矩阵中的元素。
则有
Figure FDA0002652520300000054
等式两边同时左乘
Figure FDA0002652520300000055
并结合式(21)，整理可得
Figure FDA0002652520300000056
由于第j个铰元件输出点状态矢量与第j+1个体元件输入点状态矢量相等，即
Figure FDA0002652520300000057
将式(26)代入式(25)，得
Figure FDA0002652520300000058
联立式(19)与式(27)，可得铰元件j的传递矩阵为
Figure FDA0002652520300000059
Figure FDA00026525203000000510
Figure DEST_PATH_FDA0002705925300000061
则各铰元件的传递方程为
Figure FDA0002652520300000062
式中，Uj(j＝1,3,5,7,9,11)为铰元件的传递矩阵，z1,0为工业机器人系统固定端边界点状态矢量。
6.根据权利要求5所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S5具体为：
根据各体元件传递方程(17)和各铰元件传递方程(29)，得到机器人系统总传递方程
z12,0＝U1-12z1,0 (30)
式中，U1-12＝U12…U2U1为系统总传递矩阵；
工业机器人固定端边界点和自由端边界点的边界条件为
Figure FDA0002652520300000063
将边界条件(31)代入式(30)，得到分块矩阵形式为
Figure FDA0002652520300000064
Figure FDA0002652520300000065
优化式(33)，计算边界点状态矢量中的未知状态变量
Figure FDA0002652520300000066
根据计算得到的未知状态变量获得系统固定端状态矢量，然后依次利用各元件的传递方程(17)和(29)，获得系统中各连接点的状态矢量。
7.根据权利要求6所述的一种旋转关节型工业机器人非线性动力学建模分析方法，其特征在于：所述步骤S6中铰元件θj(j＝1,3,5,7,9,11)的广义坐标的二阶导数的计算过程为：在式(18)两边同时左乘
Figure FDA0002652520300000071
得到
Figure FDA0002652520300000072
CN202010875365.4A 2020-08-27 2020-08-27 一种旋转关节型工业机器人非线性动力学建模分析方法 Active CN112001087B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010875365.4A CN112001087B (zh) 2020-08-27 2020-08-27 一种旋转关节型工业机器人非线性动力学建模分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010875365.4A CN112001087B (zh) 2020-08-27 2020-08-27 一种旋转关节型工业机器人非线性动力学建模分析方法

Publications (2)

Publication Number Publication Date
CN112001087A true CN112001087A (zh) 2020-11-27
CN112001087B CN112001087B (zh) 2021-10-19

Family

ID=73470459

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010875365.4A Active CN112001087B (zh) 2020-08-27 2020-08-27 一种旋转关节型工业机器人非线性动力学建模分析方法

Country Status (1)

Country Link
CN (1) CN112001087B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113954079A (zh) * 2021-11-23 2022-01-21 北京邮电大学 一种同质模块化机器人通用数学表征方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104793497A (zh) * 2015-04-20 2015-07-22 天津理工大学 基于多体系统离散时间传递矩阵法的机器人动力学建模方法
CN106708078A (zh) * 2017-02-21 2017-05-24 西北工业大学 一种适用于空间机器人执行器故障下的快速姿态稳定方法
CN109543264A (zh) * 2018-11-12 2019-03-29 天津理工大学 一种基于多维度重构校正的柔性多体机器人建模与求解方法
CN110340888A (zh) * 2018-10-30 2019-10-18 大连理工大学 一种空间机器人抓捕控制系统、强化学习方法及动力学建模方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104793497A (zh) * 2015-04-20 2015-07-22 天津理工大学 基于多体系统离散时间传递矩阵法的机器人动力学建模方法
CN106708078A (zh) * 2017-02-21 2017-05-24 西北工业大学 一种适用于空间机器人执行器故障下的快速姿态稳定方法
CN110340888A (zh) * 2018-10-30 2019-10-18 大连理工大学 一种空间机器人抓捕控制系统、强化学习方法及动力学建模方法
CN109543264A (zh) * 2018-11-12 2019-03-29 天津理工大学 一种基于多维度重构校正的柔性多体机器人建模与求解方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
XIAOTING RUI: "DEVELOPMENTS IN TRANSFER MATRIX METHOD FOR MULTIBODY SYSTEMS (RUI METHOD) AND ITS APPLICATIONS", 《IDETC/CIE 2018》 *
YUE GUO, WEI CHEN AND XINHUA ZHAO: "Transfer Matrix Method for Multi-body Systems Dynamics Modeling of Rescue Robot", 《IEEE》 *
郭月: "基于多刚体系统传递矩阵法的救援机器人动力学建模与规划", 《中国优秀博硕士学位论文全文数据库(硕士).信息科技辑》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113954079A (zh) * 2021-11-23 2022-01-21 北京邮电大学 一种同质模块化机器人通用数学表征方法

Also Published As

Publication number Publication date
CN112001087B (zh) 2021-10-19

Similar Documents

Publication Publication Date Title
CN103495977B (zh) 一种6r型工业机器人负载识别方法
CN111993417B (zh) 一种基于rbf神经网络的机械臂自适应阻抗控制方法
CN104723340A (zh) 基于连接和阻尼配置的柔性关节机械臂的阻抗控制方法
Zhao et al. Dynamics analysis and characteristics of the 8-PSS flexible redundant parallel manipulator
CN109657282B (zh) 一种基于拉格朗日动力学的h型运动平台建模方法
CN109634111B (zh) 一种高速重载机器人动态变形计算方法
CN104723341A (zh) 基于连接和阻尼配置的柔性关节机械臂的位置控制方法
CN102962838A (zh) 具有封闭式运动学正解的六自由度并联机构及解析方法
CN112001087B (zh) 一种旋转关节型工业机器人非线性动力学建模分析方法
Xie et al. Dynamic modeling and performance analysis of a new redundant parallel rehabilitation robot
CN110909438A (zh) 一种基于动力学模型的轻载关节型并联机器人控制方法
CN111597657B (zh) 一种旋转关节型工业机器人模态参数和振动响应计算方法
Noppeney et al. Task-space impedance control of a parallel Delta robot using dual quaternions and a neural network
CN110774286B (zh) 一种基于刚柔耦合动力学的五自由度机械手的控制方法
CN114888812A (zh) 一种基于刚度性能优化的移动机械臂站位规划方法
Zaher et al. Joints flexibility effect on the dynamic performance of robots
Liu et al. N-PD cross-coupling synchronization control based on adjacent coupling error analysis
Yang et al. Research on gravity compensation in motion control of multi-joint robot
Wang et al. Newton-Euler method for dynamic modeling and control of parallel polishing manipulator
Xu et al. Kinematic analysis of a novel type of 6-DOF hybrid serial-parallel manipulator
CN117697767A (zh) 空间站弹性关节空间机械臂控制方法及装置
CN109571483B (zh) 一种空间机械臂任务轨迹规划域构建方法
Staicu et al. Spatial parallel robots
Hovland et al. Benchmark of the 3-dof gantry-tau parallel kinematic machine
Zheng et al. Dynamic analysis and simulation of six-axis cooperative robot based on screw theory

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20220130

Address after: No. 29, Qinhuai District, Qinhuai District, Nanjing, Jiangsu

Patentee after: Nanjing University of Aeronautics and Astronautics Asset Management Co.,Ltd.

Address before: No. 29, Qinhuai District, Qinhuai District, Nanjing, Jiangsu

Patentee before: Nanjing University of Aeronautics and Astronautics

TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20220223

Address after: 211599 No. 59, Wangqiao Road, Xiongzhou street, Liuhe District, Nanjing, Jiangsu

Patentee after: Nanjing Changjiang Industrial Technology Research Institute Co.,Ltd.

Address before: No. 29, Qinhuai District, Qinhuai District, Nanjing, Jiangsu

Patentee before: Nanjing University of Aeronautics and Astronautics Asset Management Co.,Ltd.

TR01 Transfer of patent right