基于MATLAB 六轴喷涂机器人运动仿真研究

近年来，随着工业机器人技术的成熟与发展，工业机器人应用技术迈向了一个新的领域，而喷涂机器人占据其中的重要一环。

与传统的喷涂设备相比，使用机器人不仅可以极大地减少涂料对人体的伤害，还可以提高喷涂效率，且喷涂质量改善明显。

因此，许多公司都在研发适应各种环境的喷涂机器人。

然而，随着制造业的快速发展，一方面用户对工件表面的喷涂质量提出了更高的要求，另一方面需要喷涂的工件种类越来越多。

为了满足各种工件的喷涂要求，需要对机器人进行更加精准的控制。

此外，机器人运动学的研究也是为了更好地对机器人进行实时控制，以便于机器人的轨迹规划。

目前，对机器人运动学研究的学者很多，研究运动学的方法也多种多样，但是本质都是分析机器人运动轨迹、各个关节转角和其末端位姿间的逻辑关系[1-3]。

本文研究机器人运动学，以GR6150型六轴机器人为研究目标，采用坐标法建立运动学方程，探究机器人的位姿和关节转角之间的关系，运用现代计算机技术及MATLAB 辅助工具进行建模仿真与验证。

1　喷涂机器人结构分析与建立数学模型
1.1　GR6150喷涂机器人结构分析
GR6150型机器人是一款六轴喷涂机器人，机器人本体结构模型如图1所示。

该机器人是由若干活动连杆和转动关节组成，每个关节都有一组电机和减速器负责驱动。

大臂和小臂（2轴和3轴）结合处安装有3组电机，减速器小臂内部通过3个链条传动，手腕部分（4轴、5轴、6轴）由3组啮合的齿轮组成。

图1　GR6150机器人模型
该种类型的机器人整体结构比较紧凑，特别是手腕部分结构。

这种设计概念是由于考虑喷涂机器人的喷涂效率和效果与机器人的运动轨迹有关[4]。

因此，为了增大手腕的关节转动范围和减少手腕末端关节部位出现涂料漏喷，GR6150在设计时调整手腕结构的紧凑度，然后让小臂的旋转关节靠近手腕，缩短旋转关节的衔接长度，大大的改善喷涂质量。

1.2　数学模型建立
1.2.1　喷涂机器人
D-H坐标系的建立
由于GR6150机器人每个关节都是转动关节，根据关节转动的方向，采用D-H坐标法如图2所示。

图2　D-H坐标系
根据齐次坐标变换法得到相邻杆件之间的位姿矩阵，推导过程如下。

建立两个相邻连杆i-1和i的连杆坐标系{i-1}与{i}。

假设T(i-1)i表示第i连杆相对于第i-1[5]由齐次坐标变换原理可知，从i-1
齐次变换可以描述为：
()
()()(
1
000
i i i
i i
T Rot z,θTrans,,d Trans a,
−
=
即：
摘　要：喷涂机器人的运动学研究是动力学控制和轨迹规划的理论基础。

为对喷涂机器人的运动学进行具象化研究，以埃夫特公司提供的GR6150型喷涂机器人为研究对象。

GR6150型机器人是一款六轴喷涂机器人，也是应用非常广泛的多自由度喷涂机器人。

根据机器人外形结构，采用D-H法建立坐标系，对GR6150型机器人进行正、逆运动学求解和分析。

利用Matlab对GR6150型机器人建立机器人六轴连杆模型，同时对机器人进行正运动学仿真分析，从而验证了运动方程的合理性和连杆模型的正确性。

该研究为验证六轴机器人的运动学研究提供了有效方法，也为动力学研究奠定了理论基础。

关键词：喷涂机器人 运动学 仿真
（2）根据式（2），带入连杆参数，解出相邻连杆位姿矩阵T01,T12,T23,T(i-1)i,…,T(n-1)n，将它们按顺序相乘，可以求出T0n，得：
T0n=T01T12T23…T(i-1)i T(n-1)n（3）其中，n为机器人关节数，T0n为机器人手腕末端相对于基座的位姿矩阵。

1.2.2　喷涂机器人运动数学模型的建立
如图3所示，为了更好地建立模型且方便求解运动学方程，让机器人处于一个特殊位姿状态即初始位姿。

喷涂机器人在初始位姿状态下每个关节上都建立有一个空间坐标系，且坐标零点处于如图3所示的工作截面上。

在这种位姿下，在各个关节上建立D-H坐标系，使得运动方程推导过程得到较大简化，运动学求解更方便。

图3　喷涂机器人初始姿态
根据上述机器人初始位姿，建立喷涂机器人的坐标系如图4所示。

图4　喷涂机器人坐标系
喷涂机器人的D-H的参数如表1所示。

其中θi为变量，表示每个关节的转动范围。

根据GR6150实际设计尺寸，表中其余各参数值为：a3=1400mm，d4=1500mm，d5=108.7mm，d6=87mmd，β=70°。

2　GR6150机器人运动学分析
根据上述公式（2）求得运动方程以及机器人工具箱解出机器人末端位姿方程。

2.1　运动学正解
根据上述连杆坐标系和连杆参数，将表1中参数带入式（2），得：
11
11
01
00
00
0010
0001
c s
s c
T
−
=
（4）
22
22
12
00
00
0100
0001
c s
s c
T
−
=
−
（5）
3333
3333
23
0010
0001
c s a c
s c a s
T
−
=
（6）
44
44
34
4
00
00
010
0001
c s
s c
T
d
−
=
−
（7）
555
555
45
5
0001
ββ
ββ
ββ
c s c s s
s c c c s
T
s c d
−
−
=
（8）
666
666
56
6
0001
ββ
ββ
ββ
c s c s s
s c c c s
T
s c d
−−
=
−
（9）
式中，c i=cosθi，s i=sinθi，cβ=cosβ，sβ=sinβ。

根据式（3）末端位姿运动方程以及上面计算得到的6个相邻连杆的位姿矩阵，得到GR6150
末端相对于基座标系的变换位姿矩阵T06
06011223344556
000
x x x
y y y
z z z
n o a
n o a
T T T T T T T
n o a
==
经计算，矩阵中的各个参数如下：
()()()12341454112345123564112351412345123561234141236 y ββββββββ
c s c s c s c c s c s c c s s c s s c s c s c c s s s c s s =++−− +−+−− −++ （12）
()()()2345235234562345234452356234236 z ββββββββ
n s s s c c s s s c c c s c c s s c c c c c s s c s s s c c s s =−−++−++ （13）
()()()()()1234145
14123451235614123451412345123561412341236 x ββββββββ
o c c c s s c s c c c s s c c s s s s s s c c c s s c c c s c c c s c s c c s c c c s s c s c c s =−++−−+ −+−−+ −+
（14）
()()()()()1234145
12341451235612341451234145123561231234146 y ββββββββ
o s c c c s c s c s c c s c s s s s s s c c c s s s c s c c c c s s c c c c s s c s c s c c s c s =−−−+−−
+++++ −+
（15）
(
)()()4545232356
4545
23
2356234236 z βββ
β
ββ
β
β
o s s c c c s c s s s s c c c c s c c s c c s s s c c c s
=−−−+
+−+
+ （16）
()()()141234514123451235123414123 x ββββββa s s c c c s s c c c s c c c s s c s c c s s c s c s c c =−+−−
−++
（17）()()())123414512341451235123414123 y ββββββ
a s c c c s s s c s c c c c s s c s s s c s c c s s s c c =−+−+− ++−
（18）()()2345234523523423 z ββββββ
a s c s s s c c c c s s s s s c c c =+−−+ （19）
()()()()12341451412345123561234141236141234612351234123 x ββββββββp c c c s s s s s c c s c c c s c s d s c c s s c s c s c d c s c c c s d s c s d c c s d c c a =−++−−
+−−− −−−+ （20）()()()()14123451412345123551412341236123414612351234123
y ββββββββp c s s c c s c c s c s c c s s c s d s c c s c s s s s c d c s c s c c d s s s d c s s d s c a =−+−− ++−+
+−−+ （21）
()()()454523235623234234523523423623
z ββββββββp c s c s c s c s s d s s s c s d s s d s c d c s s s c c d c s a =+−+
−−−−+− （22）
为了计算方便，其中c 23=cos(θ2+θ3)，s 23=sin(θ2+θ3)。

060112233445560
1y
y y y z z z z n o a p T T T T T T T n o a p ==
（23）
[]2
26
3
6
2q pi /pi /pi /pi /pi /pi /=− （24）
()1
q ikine r,T ,q =
（25）
其中n x 、n y 、n z 、o x 、o y 、o z 、a x 、a y 、a z 、p x 、p y 、p z 在正解过程中已经求出，是已知的；最右边6个相邻连杆位姿矩阵取决于各个关节转角θ1…θ6的大小，求解时，从T 06开始求解关节变量，然后根据坐标变换使T 06表达式的各元素与式（23）中间表达式形式相等，矩阵左右元素对应相等，即可确定θ1。

在求得θ1后，式（23）两边左乘T 0-11，得：
T 0-11T 06=T 12T 23T 34T 45T 56
（26）
式中，左边为θ1和T 06各元素的函数，可以得到变量关节。

同理，重复以上求解过程，可以依次求得θ3、θ4、θ5、θ6等关节变量。

3　机器人运动学仿真3.1　建立运动模型
由第2节得出机器人运动学方程，根据机器人的连杆参数和关节转角，使用Matlab 机器人工具箱建立连杆模
型，然后输入关节变量，验证运动方程的准确性。

图5为在Matlab 中的连杆模型。

3.2　运动方程的验证
在之前通过理论分析得出运动学方程，为了验证运动方程的正确性，利用Matlab 软件中Robotics Toolbox 机器人函数进行实例验算[3]。

3.2.1　正运动学方程验证
根据关节角范围，假设θ1=π/2，θ2=π/2，θ3=-π/6，θ4=π/3， θ5=π/6，θ6=π/2，代入式（10），利用Matlab 编写程序得到正向运动学解为：
060661204893056871160054850202008114155805119084840134854150
001......T ....−
−− =
−−
（27）
调用Matlab 中Robotics Toolbox 工具箱的fkine 函数，对末端空间位姿进行求解：
T =fkine (r ,q )
其中，r 指本文研究对象GR6150-pi /6 pi /3 pi /6 pi /2]3.2.2　逆运动学方程验证
前面已知逆解求解公式，用Matlab 工具箱中q 1=ikine (r ,T ,q )函数，如表2所示[6]。

由表2可知，计算值与仿真值最大差值0.0008，误差为0.05%，在误差允许范围内，所以利用Matlab工具箱进行仿真结果与理论计算得到的逆解基本一致，证明了逆运动学方程基本正确。

3.3　末端轨迹规划仿真及分析
轨迹规划是根据作业任务的要求，事先规定机器人的操作顺序和动作过程。

轨迹规划仿真可以更详尽直观地描述工业机器人的运动过程[7-11]。

根据上述建立的仿真模型，在Matlab中调用jtraj函数构建末端从空间一点到另一点的轨迹曲线，假设在工作空间中某一点A点为初始位置，此点的位姿关节向量用q A=[0 0 0 0 0 0]表示，B点为空间中异于A的一点，此点的位姿q B=[pi/2 pi/2 -pi/6 pi/3 pi/6 pi/2]表示。

根据关节A和B的关节变量参数，现在用Matlab工具箱对机器人末端进行从A点到B点的轨迹输出，仿真结果如图6所示，得到由起点到终点的运动轨迹曲线以及位置、速度、加速度变化曲线。

图6　末端运动轨迹
在仿真采样点为50步时，机器人从A到B的空间轨迹输出图形如图6所示。

随着步数的变化，得到机器人手腕末端的位置、速度、加速度3个随着步数变化的连续光滑的曲线。

其中，由位置曲线图可以看出，机器人末端位移由零逐渐变化约1.57rad。

由速度曲线可以得到机器人末端在开始和结束位置速度均为零，在接近中间位置时速度达到最大值，此时加速度为零。

由加速曲线看出，在开始和结束位置加速度均为0，随着机器人连续运动，加速度由正值逐渐到负值，出现了正负两个极值。

该机器人末端运动位置由起点到终点的过程速度和加速度曲线连续平滑且没有断层，说明速度和驱动力没有突变。

因此，在这个过程中，机器人可以相对平稳运行，没有出现较大的振动，从而验证了GR6150型六轴机器人结构设计的合理性和运动模型的有效性。

4　结语
本文以希美埃公司生产的GR6150型六轴机器人为研究对象，用D-H方法建立坐标系进行正解和逆解的分析，求解出一种六自由度机器人的运动算法，最后用Matlab软件中的Robotics toolbox对GR6150进行建模仿真，验证了运动学算法的有效性。

通过仿真结果，证明了GR6150型六轴机器人运动学方程的有效性，可为机器人末端的动力学研究与控制提供理论基础。

参考文献
[1]王晓强，王帅军，刘建亭.基于MATLAB的IRB2400工业机
器人运动学分析[J].机床与液压，2014，（3）：54-57，40. [2]熊有伦.机器人技术基础[M].武汉：华中理工大学出版社，
1996.
[3]蔡自兴.机器人[M].北京：清华大学出版社，2000.
[4]陈伟，赵德安，梁震.喷涂机器人的喷枪轨迹优化设计与实验[J].
中国机械工程，2011，（17）：
[5]孙亮，马江，阮晓钢.
控制工程，2010，（3）：388-391.
[6]苏学满.基于matlab
设计与制造，2013，（1）：78-80.
[7]舒志兵，张海荣.
器人技术与应用，2006，（6）：39-42.
图5　滑块图和Matlab下的三维模型图
表2　仿真值与计算值
关节变量θ1θ2θ3θ4θ5θ6计算值 1.57 1.57-0.5233 1.04670.5233 1.57仿真值-1.5708-1.5708-0.5236 1.04720.5236 1.5708
[10]CORKE P I.Robotics Toolbox for Matlab（Release7.1）[EB/
OL].（2002-04-01）[2020-02-13].http://www.cat.csiro.au/cmst/ staff/pic/robot.
[11]刑迪雄，张琦.基于CATIA V5的工业机器人运动学仿真研
究[J].机械，2011，（1）：57-61.
Based on MATLAB Six Axis Spraying Robot Motion Simulation Research
GAO Fangjun1, ZHENG Lei2, XU Qiang2, TAN Xiaohong1 (1.School of Mechanical Engineering, Hefei University of Technology, Hefei 230009;2.EFORT Intelligent Equipment Co., spray robot that is widely used,According to the shape structure of the robot, the kinematical equation were established by using the modified form of D-H notation，and the forward kinematics and The inverse kinematics problem was analyzed by using the coordinate frame. Matlab was used to build a six-axis connecting rod model for GR6150 robot, and forward kinematics simulation was conducted for the robot. Thus, the rationality of the geometric model design and the correctness of the mathematical model are verified. This study provides an effective method to verify the kinematics of six-axis robot ,and lay a theoretical foundation for dynamics research.
Key words: painting robot, kinematics, simulation
（上接第52页）横截面上的风量；
（3）控制调节风扇的转速，使腔内压力达到设定值；
（4）控制主风扇，使腔内风速达到要求；
（5）待系统风速稳定后，根据排气的风速变化，再次调节排气风扇，从（2）开始新一轮调整。

该系统执行单元和控制对象均有多个。

根据上述的控制循环，周期性调整各个风机的速度，最终使4个被控参数落在要求的范围内。

由于空气流动响应的速度较慢，控制周期选取不宜过小，目前为10s调整一个风扇。

实际应用时，根据一组常用的机器参数要求手动调整优化，得到一组风扇的转速参数。

设备正常启动时，从这一组默认的参数开始运行，根据操作人员输入参数进行动态调整。

3　测试效果及总结
该设备目前已经整机运行，经苏州大学卫生与环境技术研究所进行真机检测，所测结果均符合设计要求，即前门风速0.35～0.5m/s，腔内风速0.35～0.5m/s，腔内气压±20～50bar，排气风量占比大于30%[1]，配药腔内洁净度符合GB 50457-2008中100级的要求，前门操作区域清洁度满足10000级要求[2]，且使用烟雾检查前门操作台气流方向，气流内外的交换未影响到操作区输送装置。

配药腔及前门菌落检测均为0。

该循环方式满足了化疗药物配制时对洁净度的要求，同时对人员、药物、环境做到了充分的保护。

70%内循环的方式最大限度地降低了能耗，在保证安全的前提下实现了经济效益的最大化。

参考文献
[1]国家食品药品监督管理局.YY 0569-2011Ⅱ级生物安全柜[S].
北京：中国标准出版社，2011.
[2]中华人民共和国住房和城乡建设部.GB50457-2008医药工业洁
净厂房设计规范[S].北京：中国计划出版社，2009.
Design of Air Flow System for Chemotherapy Drugs Automatic Compounding Machine
CUI Mingke
(Shanghai Swisslog Healthcare Co., Ltd., Shanghai 200127)
Abstract: Most of the chemotherapy drugs are toxic to human and their preparation process has always been harmful to nurses. The traditional solution for the nurses is to get them prepared in a biosafety cabinet in PIVAS (Pharmacy Intravenous Admixture Service), but this does not get contact with the drugs avoided completely. The company where the author works has
Key words
air purification。