昆明理工大学-计算机仿真技术课程报告

1、计算机仿真技术课程设计报告计算机仿真技术课程设计报告学 院： 机电工程学院 专 业： 机械工程及自动化 姓 名： 学 号： 专业班级： 指导教师： 刘 孝 保 2015年 6 月15日目录1 设计背景12 设计任务13 数学模型建立13.1正运动求解13.2动力学方程建立33.3拉格朗日方程求解43.3.1惯性力43.3.2离心力和科氏力53.3.3重力64 MATLAB建模仿真及验证64.1 Matlab求解流程64.2 Matlab求解关节雅克比矩阵74.3 Matlab求解其它相关矩阵84.4 Matlab求解运动参数84.5 结果对比85 结论96 课程体会9参考文献11附录121 设

2、计背景过去几十年，中国制造业迅速发展，成为以劳动密集型产业为主导的“世界工厂”。然而，随着近几年国内劳动力成本的大幅度上升，以及面临着印度等国外低廉劳动力价格竞争，迫使中国制造业产业必须采取成本更低、更具竞争力的机器人生产制作一体化模式。码垛机器人是一种利用工业机器人技术与机电一体化相结合的自动化产品，它可将粉粒状原料生产过程中的称量、装袋、码垛等一系列生产线任务实现无人工干预的自动化作业。随着物流和食品等行业的不断发展，以及工业自动化理念的推广和深入，码垛机器人发挥着越来越重要的作用，并且已逐渐受到诸多加工生产企业的青睐。码垛机器人它不仅可以准确、高效地完成码垛作业，而且可以降低工人的劳动强

3、度，提高生产效率，保障工人安全的同时，减少了物料的破损和浪费。码垛机器人的整机动态特性对于其负载、运行速度和稳定性具有重要意义，而动态性能提高的基础是其动力学模型的准确性。目前，国内机器人方面的研究多侧重于机构设计以及静态运动学分析1，对码垛机器人动力学相关研究较少，而码垛机器人动力学模型的建立对其整机动态特性研究、尺度参数优化以及控制系统研究都具有非常重要的意义。 国内对码垛机器人的研究刚刚起步，成熟的产业化码垛机器人产品还未出现2。随着自动化技术的发展和进步，工业码垛机器人在物流自动化领域的应用日趋广泛，工业码垛机器人相关技术的研究也得以深入展开。为了能够更好地实现控制，提高码垛机器人的性

4、能，以及让搬运码垛机器人安全、稳定的工作，必须要对其进行动力学分析。2 设计任务本文以ABB公司的IRB6200型号码垛机器人进行动力学仿真分析。利用上述机器人相关参数在Matlab软件上编程的方式运用拉格朗日方程得到了和SolidWorks软件运动仿真相似的结果，验证了此方法的正确性。3 数学模型建立3.1正运动求解由于拉格朗日方程求解中需要运用到关节的坐标数据，所以此处利用D-H参数3计算每个关节间和基坐标系间的变换矩阵，确定坐标间的参数关系。首先运用SoildWorks建模功能对目标码垛机器人IRB6200建模，对确定好的模型测量长度参数，从而确定机器人的D-H参数，如图1和表1所示。图

5、1: 码垛机器人表 1:D-H参数ii-1ai-1dii10001(0) 290320130.52(90)3097595.53(0)4902008874(0)590005(0)6900-1636(0)表中：i-1绕xi-1轴转到从zi-1到zi轴的角度ai-1沿xi-1轴方向从zi-1到zi轴的位移di沿zi轴方向从xi-1到xi轴的位移i绕zi轴转到从xi-1到xi轴的角度得到的广义坐标结构简图如图2所示。图2:广义坐标结构简图已知每个关节的广义转角，运用正运动求解得到相邻坐标间的变换公式。 （1）表示从坐标到坐标的变换公式4，其中第四列的前三行分别表示转换的三维方向位移变化。 （2）表示第

6、个关节相对基坐标系的三维坐标。由于Matlab编程不能模拟质量均匀分布的物体属性，所以建立的变换公式必须采用重心的坐标从而等效为重心运动，保证结果的正确性。3.2动力学方程建立 根据拉格朗日方程5,6有： （3）式中：L拉格朗日算子，总动能K和总势能U之差 （4）第i个关节的广义加速度第i个关节的广义速度第i个关节的广义位移第i个关节的广义力自由度数（关节数）因为势能U是关于高度（位置）的函数，与时间无关，故 （5）拉格朗日方程经简化后得到：（6）对于多体系统而言，为广义惯性力向量，为重力向量，用表示。 （7）式中：惯性矩阵利用式(7)简化拉格朗日方程 （8） （9）（10）式（10）中为惯性

7、力；为离心力和科氏力，用表示。将式(10)代入(6)得到动力学方程的一般形式： （11）3.3拉格朗日方程求解3.3.1惯性力由动能定义知总动能K等于每个杆件关节的直线运动动能和旋转运动动能之和：（12）式中：第i个关节质心绝对速度第i个关节质心绝对角速度第i个关节质量第i个关节质心旋转惯量在每个关节自身的坐标系下的运动速度称为广义速度，而相对机架的基坐标系则为绝对速度，两者之间的转换需要利用雅克比矩阵7。 （13） （14） （15）（16）其中：基坐标系中末端位移基坐标系中末端z轴方向位移低副类型 （17）将式(13)(15)代入式(12)得：（18）将式(7)代入上式得： （19）3.3

8、.2离心力和科氏力将分解为离心力与科氏力之和，分别由和表示。（20）其中 （21） （22）为矩阵第i行j列元素，代表惯量大小。 （23）由于已由式(19)求出，故已知，能够进一步求解离心力和科氏力。3.3.3重力 （24）重力由雅克比矩阵和重量矩阵相乘即可得到。4 MATLAB建模仿真及验证4.1 Matlab求解流程在实际工业运用中，由于码垛机器人涉及多个曲面，用SolidWorks对机器人进行准确建模耗时费力，而结构设计中的一些次要结构如连接孔等在对模型进行仿真分析时，为了提高运算效率，又不得不进行必要的简化，导致分析结果与实际存在一定误差，且简化的越多误差越大。而利用数学模型建模求解的

9、方法只需要知道机器人的重心位置、质量和转动惯量便可以得到最后的驱动力矩结果，避免了模型简化导致的分析误差、缩短了研究时间以及简化了工作流程，具有实际运用意义。Matlab软件具有强大的计算能力，由于涉及较多的微分，可利用机器人D-H参数和建立好的数学模型，在Matlab平台上进行编程求解。由SolidWorks模型可以得到码垛机器人每个关节的质量以及转动惯量。根据上述模型推导可知，只要求出雅克比矩阵，模型其余力矩皆可被求出。具体求解步骤流程图如图3所示。图3:求解驱动力矩流程图首先确定码垛机器人的正运动关节角度，基于D-H参数可以建立每个关节的变换矩阵。4.2 Matlab求解关节雅克比矩阵其

10、次根据每个关节的变换公式和低副类型可以求得雅克比矩阵。以下是计算第六个关节雅克比矩阵的程序：joint=6;syms Jw6 Jv6%第六关节%Jw6(1,1)=0;Jw6(3,6)=0;Jv6(1,1)=0;Jv6(3,6)=0;for i=1:joint Jw6(:,i)=T0ii(1:3,3);endfor i=1:jointJv6(:,i)=diff(T0ijoint(1:3,4),qq(i);endJ6=Jv6;Jw6;4.3 Matlab求解其它相关矩阵得到所有关节的雅克比矩阵之后，基于式(19)、(20)、(23)、(24)求解可以得到质量、惯性力、科氏力以及重力矩阵。例如求解质

11、量矩阵的程序如下：freedom=6;%计算质量矩阵%for i=1:freedomMi=m(i)*Jvi*Jvi+Jwi*Ic(i)*Jwi;endMM=M1+M2+M3+M4+M5+M6;4.4 Matlab求解运动参数运用Matlab自带的jtraj()函数可以求得规划路径的广义位移、速度以及加速度等运动参数。码垛机器人的起始广义坐标向量：qA=0,pi/2,0,0,0,0 ；终止点广义坐标向量：qB=-pi/4,pi/4,0,-pi/4,pi/4,-pi/2；运用语句：q,qd,qdd=jtraj(qA,qB,t);可以得到广义位移q、速度qd和加速度qdd向量，其结果均为jtraj函

12、数内部利用五次多项式插值求解得到，所以得到的运动规划路径具有先加速再减速的特性。4.5 结果对比最后基于式(11)得出每个关节的驱动力矩。SolidWorks软件中的“Motion 分析”模块可对码垛机器人添加马达进行运动仿真，由于SolidWorks马达力矩仿真结果恒大于等于零，需对Matlab得到的力矩数据取绝对值。由于篇幅有限，仅对较为重要的第一关节和第二关节力矩大小仿真对比，数据如表2、表3所示，表中横坐标数据单位为秒，力矩单位为Nm。分析如图4和图5所示。由图可见，两种方法仿真的误差较小，验证了Matlab拉格朗日方程仿真的可靠性。图4:关节1仿真曲线对比图5:关节2仿真曲线对比5

13、结论本文提出的利用雅克比矩阵求解机器人运动方程的方法，实际上是对拉格朗日方程简化求解从而得到码垛机器人的动力学参数结果，实现了在MATLAB软件上快速准确的建立动力学模型。与SolidWorks运动仿真结果对比可以看出，雅克比矩阵求解拉格朗日方程的算法具有较大的可行性，在精度要求不是特别高的场合，例如机器人重要的前臂和大臂的设计控制，此方法可以缩短实物的设计时间，减小样机的开发成本，进一步优化现有码垛机器人的结构参数，具有较强的实际意义。6 课程体会在科技高度发展的今天，计算机在人们之中的作用越来越突出。而matlab作为一种计算机工程软件，我们学习它，有助于我们更好的了解计算机，与计算机进行

14、交流，因此，matlab的学习对我们尤其重要。自从大二开始参加课外比赛接触matlab以来，我们学习matlab的时间已经一年多了，对matlab也算得上半个入门者，期间也写过一些课程或课外竞赛编过一些代码，但是也仅仅是停留在入门的范围，对里面的好多东西还是很陌生，更多的在运用起来的时候还是感到很棘手。所以巩固加深学习matlab，不管是在专业运用层面还是对于以后学习更高级的语言来说重要性都不言而喻。在这一个学期的课程学习中，老师会在课堂上要讲授许多关于matlab的语法规则，包括GUI、Simulink模块功能。在课堂上感觉已经听懂的指令或使用方法，但到完成作业过程中由发现还是没掌握，这便是

15、由于我们不是十分深层次的了解接触程序设计，缺乏程序设计的实践所致。没办法，就只能将老师课堂上讲解的知识点找出来对应着编写代码。后来随着上课的深入，我们会在课堂的授课中，将自己许多不理解的内容点以及函数名记录下来，在课后通过软件自带的help教程以及网络搜索努力为自己的知识盲点漏洞打补丁。通过多次课后练习，对于语法知识有了较感性的认识，在理解的基础上就自然而然地掌握了许多matlab的语法规定。并且对于一些内容自己认为在课堂上听懂了，但在自己的实际练习操作中会发现许多的问题和偏差，学习matlab不能停留在学习它的语法规则，而是要利用学到的知识编写matlab程序，解决实际问题。即把matlab

16、作为工具，描述解决实际问题的步骤，由计算机帮助我们解题。只有通过实践编程才能检验自己是否掌握matlab。上机实验是学习程序设计语言必不可少的实践环节，特别是matlab灵活、简洁，更需要通过编程的实践来真正掌握它。对于程序设计语言的学习目的，可以概括为学习语法规定、掌握程序设计方法、提高程序开发能力，这些都必须通过充分的实际上机操作才能完成。因为学时所限，课程不能安排过多的统一上机实验，所以我会有效地利用课后时间积极尝试以及运用，尽快掌握用matlab开发程序的能力，为今后的继续学习打下一个良好的基础。完成程序的编写，决不意味着万事大吉。你认为万无一失的程序，实际运行时可能不断出现麻烦。如编

17、译程序检测出一大堆错误。有时程序本身不存在语法错误，也能够顺利运行，但是运行结果显然是错误的。就如老师上课常说的一句话，程序实际上是调出来的。我们在资料书上了解到开发环境所提供的编译系统无法发现这种程序逻辑错误，只能靠自己的经验分析判断错误所在。程序的调试是一个技巧性很强的工作，对于初学者来说，尽快掌握程序调试方法是非常重要的。有时候一个消耗你几个小时时间的小小错误。所以发现问题解决问题也是学习编程能力的一个非常好的方法。自己编好程序调试运行时，可能有很多你想不到的情况发生，通过解决这些问题，可以逐步提高自己对matlab的理解和程序开发能力。通过一个学期的matlab软件计算机实践学习，我们

18、了解了一些关于matlab的知识，理解巩固了我们matlab的理论知识，对我们将来到社会工作将会有莫大的帮助。在我看来无论学习什么语言，最重要的是掌握习编程思想，不仅仅是编程，其他任何方面思想是最为重要的，思维模式的养成和成功有着深深的联系，所以这一点对于处于学生阶段的我们来说非常重要。在这期间，我们也查找了不少资料，发现很多知识是我们从没有接触过的，我们并不了解，所以借此增长了自己的见识，了解了更多关于它的应用。在编写程序代码时，需要什么函数，需要什么模块得一步一步的来，慢慢地研究，每天积累点，不能急于求成，有问题就应该去查资料，讨论。这段经历让我知道，只要你努力，任何东西都不会太难。所以感

19、谢刘老师通过这一门课程让我们班级对matlab软件从最开始的接触到了解最后直至对知识的深入，更是通过任务的合理安排与布置，让我们掌握知识的同时，更是提高了我们能力，培养了我们思维模式的养成，相信这一次课程经历会在今后生活的各方面起着重要作用！参考文献1 张良安, 单家正, 马寅东, 刘凤臣, 解安东, 刘黎明. PT1300型局部闭链式码垛机器人动力学分析J. 山东科技大学学报(自然科学报), 2013, 32(2)2 李晓刚, 刘晋浩. 码垛机器人的研究与应用现状、问题及对策J. 包装工程, 2011, 32(2):96-102.3 干敏耀, 马骏骑, 陈永星, 范启印. 基于MATLAB的

20、PUMA机器人运动学仿真J. 昆明理工大学学报(理工版), 2003, 28(6)4 聂晓波, 陈亚峰. 焊接机器人运动学分析J,自动化仪表. 2013: 34(9).5 焦磊涛, 张伟, 仲军, 樊继壮, 王猛, 仿青蛙跳跃机器人的动力学分析J. 哈尔滨商业大学学报(自然科学版), 2013, 29(6)6 朱斌, 张楠, 胡南, 李宏峰, 煤矿救援机器人的多体动力学J. 辽宁工程技术大学学报(自然科学版), 2013,32(10)7 米士彬, 金振林. 基于雅克比矩阵求解并联机器人位置正解方法J. 燕山大学学报, 2011, 35(5).8 王田苗, 陶永, 我国工业机器人技术现状与产业化

21、发展战略J. 机械工程学报, 2014, 50(9)7 解本铭, 王永浩, 打磨机操作臂运动仿真研究J. 中国名航大学学报, 2009, 27(1)8 罗家佳, 胡国清, 基于MATLAB的机器人运动仿真研究J. 厦门大学学报(自然科学版), 2005, 44(5)11 翟敬梅, 康博, 张铁, 六自由度喷涂机器人动力学分析及仿真J. 机械设计与制造, 2012(1)附录%L=LINK(alpha A theta D Sigma，modified)，%连杆的前四个元素依次为扭角，长度，夹角（关节角），偏置（距离）%最后一个Sigma为0（转动关节）或1（移动关节） %IRB 6620 modc

22、lear %清除变量 L1=link( 0 0 0 0 0,mod); %1真实关节1L2=link(pi/2 320 pi/2 130.5 0,mod); %3真实关节2L3=link( 0 975 0 95.5 0,mod); %4真实关节3L4=link(pi/2 200 0 887 0,mod); %5辅助关节L5=link(pi/2 0 0 0 0,mod); %8真实关节5L6=link(pi/2 0 0 -163 0,mod);%9真实关节6L7=link(0 0 0 -19 0,mod);%9真实关节6r=robot(L1 L2 L3 L4 L5 L6 L7);%构建机器人r.

23、name=IRB6620;%命名close alldrivebot(r); freedom=6; %自由度m=220.5 110.7 134.3 97 16.7 2.4;Ic=12.2 11 3.7 1.3 0.12 0.011; g=0 0 -9.8; qA=0 pi/2 0 0 0 0 ; %起始点关节空间矢量qB=-pi/4 , pi/4 , 0 , -pi/4 , pi/4 , -pi/2 ; %终止点关节空间矢量%plot(r,qA)%plot(r,qB)t=0:0.076:2; % 产生时间向量 tt=0:0.04:2;q,qd,qdd=jtraj(qA,qB,t); plot(r

24、,q) % %Ti-1 i= cos(qi) -sin(qi) 0 Ai-1% sin(qi)cos(ai-1) cos(qi)cos(ai-1) -sin(ai-1) -sin(ai-1)Di% sin(qi)sin(ai-1) cos(qi)sin(ai-1) cos(ai-1) cos(ai-1)Di% 0 0 0 1; q1= sym(q1,real);q2= sym(q2,real);q3= sym(q3,real);q4= sym(q4,real);q5= sym(q5,real);q6= sym(q6,real);qq=q1 q2 q3 q4 q5 q6;T01=cos(q1)

25、-sin(q1) 0 0 sin(q1) cos(q1) 0 0 0 0 1 0 0 0 0 1 ;T12=cos(q2) -sin(q2) 0 0.32 0 0 -1 -0.1305 sin(q2) cos(q2) 0 0 0 0 0 1 ;T23= cos(q3) -sin(q3) 0 0.975 sin(q3) cos(q3) 0 0 0 0 1 0.0955 0 0 0 1 ;T34=cos(q4) -sin(q4) 0 0.2 0 0 -1 -0.887 sin(q4) cos(q4) 0 0 0 0 0 1; T555=cos(q4) -sin(q4) 0 0.2 0 0 -1 -

26、0.573 sin(q4) cos(q4) 0 0 0 0 0 1;T45=cos(q5) -sin(q5) 0 0 0 0 -1 0 sin(q5) cos(q5) 0 0 0 0 0 1 ;T56=cos(q6) -sin(q6) 0 0 0 0 -1 0.163 sin(q6) cos(q6) 0 0 0 0 0 1 ;T67= cos(0) -sin(0) 0 0 sin(0) cos(0) 0 0 0 0 1 -0.02 0 0 0 1; T12T=cos(q2) -sin(q2) 0 0.149 0 0 -1 -0.033 sin(q2) cos(q2) 0 0 0 0 0 1 ;

27、 T23T= cos(q3) -sin(q3) 0 0.416 sin(q3) cos(q3) 0 0 0 0 1 0.088 0 0 0 1 ; T34T=cos(q4) -sin(q4) 0 0.118 0 0 -1 0.048 sin(q4) cos(q4) 0 0 0 0 0 1; T45T=cos(q5) -sin(q5) 0 0 0 0 -1 0 sin(q5) cos(q5) 0 0 0 0 0 1 ; T56T=cos(q6) -sin(q6) 0 0 0 0 -1 0.36 sin(q6) cos(q6) 0 0 0 0 0 1 ; T67T= cos(0) -sin(0)

28、0 0 sin(0) cos(0) 0 0 0 0 1 -0.019 0 0 0 1;T02=T01*T12;T03=T02*T23;T04=T03*T34;T05=T04*T45;T06=T05*T56;T07=T06*T67; T02T=T01*T12T; T03T=T02*T23T; T04T=T03*T34T; %T05T=T04*T45T; T05T=T03*T555*T45T; T06T=T05*T56T; T07T=T06*T67T;T0i=T02T T03T T04T T05T T06T T07T; %计算雅克比% %对于第一关节%joint=1;syms Jw1 Jv1Jw1

29、(1,1)=0;Jw1(3,freedom)=0;Jv1(1,1)=0;Jv1(3,freedom)=0;for i=1:joint Jw1(:,i)=T0ii(1:3,3);endfor i=1:joint Jv1(:,i)=diff(T0ijoint(1:3,4),qq(i);endJ1=Jv1;Jw1;%J1=simple(J1); %对于第二关节%joint=2;syms Jw2 Jv2Jw2(1,1)=0;Jw2(3,freedom)=0;Jv2(1,1)=0;Jv2(3,freedom)=0;for i=1:joint Jw2(:,i)=T0ii(1:3,3);endfor i=1

30、:joint Jv2(:,i)=diff(T0ijoint(1:3,4),qq(i);endJ2=Jv2;Jw2;%J2=simple(J2); %对于第三关节%joint=3;syms Jw3 Jv3Jw3(1,1)=0;Jw3(3,freedom)=0;Jv3(1,1)=0;Jv3(3,freedom)=0;for i=1:joint Jw3(:,i)=T0ii(1:3,3);endfor i=1:joint Jv3(:,i)=diff(T0ijoint(1:3,4),qq(i);endJ3=Jv3;Jw3;%J3=simple(J3); %对于第四关节%joint=4;syms Jw4

31、Jv4Jw4(1,1)=0;Jw4(3,freedom)=0;Jv4(1,1)=0;Jv4(3,freedom)=0;for i=1:joint Jw4(:,i)=T0ii(1:3,3);endfor i=1:joint Jv4(:,i)=diff(T0ijoint(1:3,4),qq(i);endJ4=Jv4;Jw4;%J4=simple(J4); %对于第五关节%joint=5;syms Jw5 Jv5Jw5(1,1)=0;Jw5(3,freedom)=0;Jv5(1,1)=0;Jv5(3,freedom)=0;for i=1:joint Jw5(:,i)=T0ii(1:3,3);endf

32、or i=1:joint Jv5(:,i)=diff(T0ijoint(1:3,4),qq(i);endJ5=Jv5;Jw5;%J5=simple(J5); %对于第六关节%joint=6;syms Jw6 Jv6Jw6(1,1)=0;Jw6(3,freedom)=0;Jv6(1,1)=0;Jv6(3,freedom)=0;for i=1:joint Jw6(:,i)=T0ii(1:3,3);endfor i=1:joint Jv6(:,i)=diff(T0ijoint(1:3,4),qq(i);endJ6=Jv6;Jw6;%J6=simple(J6); Jv=Jv1 Jv2 Jv3 Jv4

33、Jv5 Jv6;Jw=Jw1 Jw2 Jw3 Jw4 Jw5 Jw6;J=J1 J2 J3 J4 J5 J6; freedom=6;%计算质量矩阵%for i=1:freedom Mi=m(i)*Jvi*Jvi+Jwi*Ic(i)*Jwi;endMM=M1+M2+M3+M4+M5+M6;%m ijk%for ii=1:freedom for jj=1:freedom for kk=1:freedom MMMii,jj,kk=diff(MM(ii,jj),qq(kk); end endend%b ijk%for ii=1:freedom for jj=1:freedom for kk=1:fre

34、edom bii,jj,kk=(MMMii,jj,kk+MMMii,kk,jj-MMMjj,kk,ii)/2; end endend%C(q)%for i=1:freedom for j=1:freedom Cq(i,j)=bi,j,j; endend %B(q)%for i=1:freedom jj=1; for j=1:freedom for k=j+1:freedom Bq(i,jj)=2*bi,j,k; jj=jj+1; end endend%G(q)%syms G G=-Jv1*g*m(1)+Jv2*g*m(2)+Jv3*g*m(3)+Jv4*g*m(4)+Jv5*g*m(5)+Jv6*g*m(6); %*循环开始*for xunhuan=1:length(q) qture=q(xunhuan,:); qdture=qd(xunhuan,:); qddture=qdd(xunhuan,:);