基于Matlab的ER50A型拆捆機器人運動特性分析

仝浩源，馬立東，孟智娟

(太原科技大學 機械工程學院，太原 030024)

當前，在特鋼棒材精整生產線上需要對管棒材進行打捆和拆捆操作。打捆已經有現成的專機，而拆捆需要依賴人工操作，勞動強度大，成本高。因此利用一種拆捆機器人，完成智能化自動拆捆具有十分重大的研究意義。將拆捆裝置安裝在工業機器人末端法蘭盤上，就可以將其組裝為性價比高的智能拆捆機器人。由于機器人造價高、工業生產環境復雜，對機器人進行拆捆作業前可先進行仿真實驗，以便讓研究人員提前了解機器人的運行軌跡和工作情況，掌握機器人運動特性的相關數據[1]。對機器人的運動特性進行相關研究，可以為之后的機器人技術分析和運動控制打好一定的基礎，是開展機器人應用的必經之路，對后續工程的研發起著至關重要的作用[2]。

本文通過MDH參數法建立ER50A拆捆機器人模型，根據機器人的DH參數和連接桿之間的位姿變換矩陣推導出ER50A運動學方程，通過五次多項式插值算法對拆捆作業中任意兩點進行軌跡規劃，在MATLAB軟件中安裝機器人工具箱Robotics-Toolbox[3]，利用工具箱中的機器人仿真功能建立機器人模型。利用建立的機器人模型對運動學方程進行了仿真驗證，根據蒙特卡洛數學算法對拆捆機器人的工作空間進行仿真模擬，通過軌跡規劃仿真得出機器人關節運動曲線，仿真結果表明了理論推導的準確性，為下一步拆捆作業的進行提供了可靠的參考和指導。

1 ER50A機器人建模

1.1 ER50A拆捆機器人實物圖

如圖1所示，ER50A拆捆機器人由拆捆末端和機器人本體構成，該拆捆機器人本體具有6個自由度，機器人負載為50 kg，重復定位精度可達到0.08 mm.

圖1 ER50A拆捆機器人結構模型

1.2 確定ER50A機器人的D-H參數

1955年，Denavit和Hartenberg二人一起提出了一種經典的機器人建模方法，被人們稱為D-H法[4]，由于該方法坐標系序號和對應關節軸的序號不一致，故有關學者對其進行改造，改造后的方法被稱為MDH參數法。在MDH參數法中，各關節坐標系建立方法為：zi軸指向關節i的軸線方向，逆時針方向為正；xi軸沿zi和zi+1軸的公垂線方向朝向zi+1軸的方向為正；xi軸和zi軸確定后即可確定yi軸[5]，由此建立如圖2所示的D-H坐標系。

圖2 ER50A機器人D-H坐標系

在圖2中，{0}坐標系建立在基座上，{1}坐標系建立在關節1處，依次類推。其中各連桿參數分別代表如下：

ai—連桿長度，zi和zi+1軸之間的最小距離；

θi—繞zi軸由xi-1轉向xi的關節角；

di—連桿距離，即ai和ai-1之間的距離；

ai—連桿扭角，即zi和zi+1軸之間的夾角；

通過圖2建立的ER50A機器人模型確定表1中的數據，如表1所示。

表1 ER50A拆捆機器人d-H參數表

在空間中兩個物體之間的位置關系可以通過移動和轉動來描述，MDH參數法的主要作用就是通過矩陣方程來確定兩個連桿之間的位置關系，根據矩陣依次相乘運算可得出機械臂末端和基坐標系之間的齊次變化矩陣關系。

從連桿坐標系{i-1}到{i}，被定義為基本平移和旋轉，經歷的變換為：繞zi軸旋轉θi、沿zi軸平移di、沿xi軸平移ai-1和繞xi軸旋轉ai-1，用公式表述如下：

i-1iT=Rz(θi)×Dz(di)×Dx(ai-1)×Rx(ai-1)

(1)

它可展開為：

i-1iT=

1.3 ER50A型機器人仿真模型的建立

根據表1中的D-H數據，在Matlab中利用Robotic-Toolbox提供的Link函數建立ER50A拆捆機器人結構模型，link函數里面的五個參數為D-H表中的數據，程序如下：

startup_rvc;

% theta d a alpha offset

L(1)= Link([0，0.304，0，0，0]，'modified');

L(2)= Link([0，0，0.220，-pi/2，0]，'modified');

L(3)= Link([0，0，0.900，0，0]，'modified');

L(4)= Link([0，0.2535，0.160，-pi/2，0]，'modified');

L(5)= Link([0，0，0，pi/2，0]，'modified');

L(6)= Link([0，0.200，0，-pi/2，0]，'modified');

ER50A = SerialLink(L，'name'，'ER50A機械臂')；

ER50A.qlim=[-pi，pi;-13*pi/18，7*pi/18;-4*pi/9，35*pi/36;-2*pi，2*pi;-23*pi/36，23*pi/36;-2.5*pi，2.5*pi]

ER50A.display();%顯示機器人DH參數表

ER50A.teach;%各關節驅動指令

theta=[0，-pi/2，0，0，0，0];

ER50A.plot(theta);%SerialLink類函數

thy = isrevolute(L2);%驗證關節2是否為旋轉關節

如圖3所示為仿真所建機器人模型：

圖3 ER50A機器人仿真模型

在圖3中，界面左上方表示機器人末端在基坐標系下的坐標值(笛卡爾空間)，左下為各關節位置調節滑塊，右側為ER50A三維模型，通過左側滑塊的來回移動或者直接輸入數值可以驅動機器人在三維空間中運動。

2 機器人運動學

2.1 機器人正運動學研究

(2)

公式(2)中各個參數具體計算如下：

r11=c1[c23(c4c5c6-s4s5)-s22s5c5]+s1(s4c5c6+c4s6)

r21=s1[c22(c4c5c6-s4s6)-s23s5c6]+c1(s4c5c6+c4s6)

r31=-s23(c4c5c6-s4s6)-c23s5s6

r12=c1[c23(-c4c5c6-s4s6)+s23s5c6]+s1(c5c6-s4c5s6)

r22=s1[s23(c4c5c6-s4s6)+s23s5c6]+c1(c4c6-s4c5s6)

r32=-s23(-c4c5c6-s4s6)+s23s5c6

r13=-c1(c23c4c5+s23c5)-s1s4s5

r23=-s1(c23c4c5+s23c5)+c1c4s5s5

r33=s23c4s5-c23c5

PX=c1[a2c2+a3c23-d4s22-d3s1

PY=s1[a2c2+a3c23-d4s22-d3c1

PZ=-a3s22-a2s2-c23d4

現給定機器人6個關節的旋轉角度分別為：θ1=0,θ1=-pi/3,θ3=pi/4,θ4=0,θ5=-pi/3,θ6=-pi/4.代入公式(2)中根據Matlab編寫程序解得：

T=0.183 0 0.183 0 0.965 9 1.083

0.707 1 -0.707 1 0 0

0.683 0 0.683 0-0.258 8 0.828 2

0 0 0 1

利用Matlab中的fkine函數對ER50A機器人進行正運動學仿真，設初始六個關節變量為theta，則有向量：

theta=[0，-pi/3，pi/4，0，-pi/3，-pi/4]

調用fkine函數：

T=ER50A.fkine(theta)%正運動位姿矩陣；

在相同旋轉關節角下，通過Matlab編程得到的機器人運動學正解和使用fkine函數解的結果完全一致，由此可以證明所建正運動學方程準確、合理。

2.2 機器人逆運動學分析和仿真驗證

在實際拆捆過程中，往往是已知末端裝置的空間位姿來求關節的旋轉變量。對機械臂進行控制以及在任務空間進行相應的軌跡規劃操作，都需要提前知道實際的關節轉動變量，而機器人逆解運算就是已知末端裝置位姿反求關節轉動量[6]。

封閉解法是機器人求逆解的方法之一，一般使用其中的代數解法對機器人進行逆解運算，其主要原理就是通過一步步的矩陣變換，在矩陣中找到一一對應的三角函數關系來求解各個變量[7]。

(7)

對于逆運動學的解，可分為以下三種情況：

(1)無解。作業空間中的位置不在機器人工作空間內，機器人無法到達該位置；再者，實際關節角都有一定的轉動范圍，超出這個范圍機器人也無法到達指定位姿，因此不存在解。

(2)唯一解。機械臂只能從一個方向到達指定位姿，此時存在唯一解。

(3)多個解。在各自關節的運動范圍內，機械臂可以從多個方向到達同一個位姿，此時存在多解。這時就要考慮選擇選擇其中一組最合適的解：一要考慮“路程最短”；二要考慮避障。

本次仿真使用的已知位姿矩陣為2.1節中的T矩陣，結合表1中的D-H參數和上述逆解公式，通過MATLAB編程計算，在限制ER50A機器人各關節運動范圍的情況下，可以得到八組解。然后在Matlab中通過機械臂仿真ikine逆解函數對ER50A機器人6關節進行求解計算，程序如下：

q1=ER50A.ikine(T);

所得結果為：

q1=[0.000 0 -1.042 7 0.785 4 0.000 0

-1.047 2 -0.785 4];

用MATLAB編程得到的一組解與使用fkine函數仿真所得q1的值相同。由此可以驗證所建機器人模型的正確性以及逆解算法的有效性。

3 ER50A機器人工作空間分析和仿真

機器人的工作空間決定了機器人末端能否到達任務規劃中所需的位置，他是末端在笛卡爾空間中點的集合所構成的云圖，對其進行比必要的分析有利于為之后的軌跡規劃提供參考[8]。

串聯型機械臂的工作空間大小主要通過前三個關節的位置轉動進行變化，本文所用的ER50A機械臂就是一款普通的串聯機器人，所以本文利用蒙特卡洛算法模擬出機械臂前三軸可行的工作空間模型圖。在實際生活中，我們經常使用蒙特卡洛算法去解決一些數學問題，它是一種隨機抽樣的數值方法，在機械臂研究中，通常被用來分析機器人的工作空間范圍。

如圖4所示是在MATLAB中用蒙特卡羅算法模擬ER50A三維工作空間的點云圖像，實驗代碼如下：

圖4 ER50A三維工作空間

A=unifrnd(-pi，pi，[1，3 000]);%第一關節變量限位

B=unifrnd(-13*pi/18，7*pi/18，[1，3 000]);%第二關節限位

C=unifrnd(-4*pi/9，35*pi/36，[1，3 000]);%第三關節限位

G= cell(30 000，3);%建立元胞數組

for n = 1:30 000

G{n} =[A(n)B(n)C(n)];

end %產生30 000組隨機點

H1=cell2mat(G)；%將元胞數組轉化為矩陣

T=double(robot.fkine(H1)); %機械臂正解

在圖4中灰色部分表示末端執行器能夠到達的空間位置，該仿真結果與機器人廠商提供的機器人活動范圍一致，由此可以證明所建機器人模型的準確性，同時可以為下一步機器人軌跡規劃的活動范圍提供參考。

4 機器人軌跡規劃和仿真

4.1 軌跡規劃

根據實際生產過程中的拆捆工況需求對機器人進行路徑規劃。軌跡規劃是在已有路徑的基礎上對機器人運行路徑點插值擬合，然后根據插值函數計算出機器人工作過程中的速度和加速度變化情況，力求讓機器人行走過程更加平滑[9]。目前，軌跡規劃的方法很多，其中多項式插值法應用最為廣泛[10]。在多項式插值中，三次多項式插值相對簡單，但生成的軌跡加速度是不連續的[11]。故本文采用高階五次多項式插值，用θ(t)表示執行器隨時間的運動軌跡，通過一階導數和二階導數可以得到終端速度和加速度的變化趨勢，如公式(8)、公式(9)、公式(10)所示：

θ(t)=a0+a1t+a2t2+a3t3+a4t4+a5r5

(8)

(9)

(10)

公式(8)可知五次多項式中共有6個未知數(a0，a1，a2，a3,a4,a5)，求解方程需要有6個約束條件，故需要給出機器人始末位置、始末速度和始末加速度約束，將其分別帶入公式(8)、公式(9)和公式(10)中，約束方程如下：

θ0、θF為已知始末位置，滿足如下約束條件：

θ0=a0

(11)

θf=a0+a1t+a2t2+a3t3+a4t4+a5t5

(12)

(13)

(14)

(15)

(16)

由以上6個方程可以解出6個未知變量，由此可將機器人軌跡運動的五次多項式求解出來。

4.2 機器人軌跡規劃仿真

根據第4.1節所述，為了直觀看出仿真結果是否與預想效果一致，不讓機器人第四個關節轉動，將機器人啟動位置q0設置為[0，0，0，0，0，0]，目的地q1位置坐標[pi/2，-pi/3，pi/4，0，pi/3，-pi/4]，運動時間t為[0:0.5:16].在Matlab中輸入[q，qd，qdd]=jtraj(qA，qAB，t)代碼可以得到如圖5所示的關節轉動圖，對機器人運動過程中的速度、加速度曲線以及在笛卡爾空間中的軌跡進行仿真，得到如圖6-圖8所示的圖形。

圖5 角度變化圖

圖6 角速度變化曲線

圖7 角加速度變化曲線

圖8 笛卡爾空間軌跡曲線

從圖8可以清晰的看到機器人從q0到q1的軌跡變化曲線；圖5為使用五次多項式插值后機械臂關節角度從出發點到目的地隨時間變化的運動過程；圖6、圖7分別為各個關節在轉動過程中速度、加速度平滑變化圖，從圖5-圖7中可以看到第四關節一直為0，各曲線變化平滑、無拐點和突變，證明機器人在拆捆作業過程中結構穩定，不會產生大的抖動，同時綜合圖5、圖6、圖7可以看到第四關節位移、速度、加速度均為零，最終證明了仿真結果的正確性。

5 結論

(1)以ER50A機器人為研究對象，通過MDH參數法建立了機器人結構模型，根據模型建立兩桿之間的矩陣變化方程，為下一步運動學分析提供了理論指導。

(2)通過建立的ER50A數學模型，求出了機器人的正解、逆解，建立了末端與基座標系之間的矩陣變化關系，并通過MATLAB仿真分析證明了其準確性。

(3)根據蒙特卡洛算法模擬出機械臂可行的工作空間點云圖，證明了ER50A機器人的工作空間能夠滿足一定的任務需求。

(4)采用高階五次多項式插值算法對機器人末端進行軌跡規劃，利用仿真軟件得到各關節的實時變化規律，節約了研究成本，同時仿真結果說明了機器人結構參數的合理性，為下一步拆捆作業的進行提供了可靠的參考和指導。