朱漢子, 盧劍偉, 陳新法, 尤 昕

(1.合肥工業(yè)大學(xué) 機(jī)械工程學(xué)院,安徽 合肥 230009; 2.合肥工業(yè)大學(xué) 汽車與交通工程學(xué)院,安徽 合肥 230009; 3.合肥泰禾光電科技股份有限公司,安徽 合肥 231200)

隨著機(jī)械臂應(yīng)用領(lǐng)域的增加,高速、高精度是機(jī)械臂追求的重要技術(shù)性能[1],全面掌握其動(dòng)力學(xué)性能并進(jìn)行精準(zhǔn)控制是高性能機(jī)械臂研發(fā)的關(guān)鍵。其中,機(jī)械臂動(dòng)力學(xué)參數(shù)的辨識(shí)是掌握其動(dòng)力學(xué)行性能并進(jìn)行控制的重要前提。但多關(guān)節(jié)機(jī)械臂動(dòng)力學(xué)模型自身較為復(fù)雜,不是所有的參數(shù)都能夠被辨識(shí)出來[2]。因此,從標(biāo)準(zhǔn)動(dòng)力學(xué)參數(shù)中篩選確定最小動(dòng)力學(xué)參數(shù)集,不僅可以降低動(dòng)力學(xué)參數(shù)辨識(shí)的復(fù)雜度,還可以提高基于動(dòng)力學(xué)模型的機(jī)械臂控制的快速性和魯棒性[3]。

國外研究人員較早開始關(guān)注多關(guān)節(jié)機(jī)械臂最小動(dòng)力學(xué)參數(shù)的集合問題[4]。文獻(xiàn)[5]將動(dòng)力學(xué)參數(shù)融入機(jī)械臂動(dòng)力學(xué)模型中構(gòu)建新的公式,利用遞歸的閉環(huán)關(guān)系法求解所需的動(dòng)力學(xué)參數(shù)組合;文獻(xiàn)[6]通過對(duì)機(jī)械臂動(dòng)力學(xué)參數(shù)回歸矩陣的分析,總結(jié)相應(yīng)的動(dòng)力學(xué)參數(shù)組合規(guī)則,并根據(jù)機(jī)械臂關(guān)節(jié)間的運(yùn)動(dòng)關(guān)系,按照組合規(guī)則確定最小動(dòng)力學(xué)參數(shù)集[7]。但是上述獲得最小動(dòng)力學(xué)參數(shù)集的方法還有不足:一方面需要較長時(shí)間分析動(dòng)力學(xué)參數(shù)的回歸矩陣來確定動(dòng)力學(xué)參數(shù)的組合關(guān)系;另一方面在實(shí)際應(yīng)用過程中根據(jù)組合規(guī)則確定最小參數(shù)集的過程繁瑣。

為此,本文在考慮關(guān)節(jié)摩擦和關(guān)節(jié)電機(jī)轉(zhuǎn)動(dòng)慣量基礎(chǔ)上,首先利用MDH(modified Denavit-Hartenberg)方法建立多關(guān)節(jié)機(jī)械臂運(yùn)動(dòng)學(xué)模型,采用遞推式牛頓歐拉方法建立其動(dòng)力學(xué)模型;然后利用隨機(jī)數(shù)模擬機(jī)械臂關(guān)節(jié)的運(yùn)動(dòng)學(xué)參數(shù),設(shè)計(jì)了一個(gè)列變換旋轉(zhuǎn)矩陣,通過該矩陣重組動(dòng)力學(xué)參數(shù)的廣義回歸矩陣,確定回歸矩陣各列的獨(dú)立性和相關(guān)性;最后對(duì)重組的動(dòng)力學(xué)參數(shù)的回歸矩陣進(jìn)行正交三角(QR)分解,從標(biāo)準(zhǔn)動(dòng)力學(xué)參數(shù)中確定最小動(dòng)力學(xué)參數(shù)集及其對(duì)應(yīng)的最簡回歸矩陣。

1 機(jī)械臂動(dòng)力學(xué)模型

六關(guān)節(jié)機(jī)械臂簡圖如圖1所示。

圖1 六關(guān)節(jié)機(jī)械臂MDH簡圖

本文利用MDH方法建立其模型,并采用牛頓歐拉方法建立包含關(guān)節(jié)電機(jī)慣性量和關(guān)節(jié)摩擦在內(nèi)的機(jī)械臂動(dòng)力學(xué)模型[8]。

(1)

各關(guān)節(jié)電機(jī)轉(zhuǎn)動(dòng)慣量Ia的力矩τdriver如下:

(2)

庫侖-黏滯摩擦所引起的力矩τfriction如下:

(3)

其中:fv為黏滯摩擦系數(shù);fc為庫侖摩擦系數(shù)。

以上三者是線性化的,通過線性疊加獲得機(jī)器人關(guān)節(jié)力矩τall的公式[9]，具體如下:

τall=τ+τdriver+τfriction

(4)

整個(gè)機(jī)器人動(dòng)力學(xué)參數(shù)φ為:

φ=[φ1φ2…φn]T

(5)

第i連桿標(biāo)準(zhǔn)動(dòng)力學(xué)參數(shù)φi如下:

φi=[LxxiLxyiLxziLyyiLyzi

lxilyilzimiIaifcifvi]

(6)

其中

Lxyi=Ixyi-mirxiryi;

Lxzi=Ixzi-mirxirzi;

Lyzi=Iyzi-miryirzi;

lxi=mirxi;lyi=miryi;lzi=mirzi。

連桿i與MDH坐標(biāo)系的i關(guān)系如圖2所示。

圖2 連桿i與MDH坐標(biāo)系的i關(guān)系

機(jī)械臂的第i連桿質(zhì)心Ci的轉(zhuǎn)動(dòng)慣量為[IxxiIxyiIxziIxxiIyyiIyziIzzi],機(jī)械臂的第i連桿質(zhì)心Ci相對(duì)第i連桿MDH坐標(biāo)原點(diǎn)oi的質(zhì)心位置為[rxi,ryi,rzi]。采用平行軸定理消去機(jī)械臂動(dòng)力學(xué)的非線性項(xiàng),由(6)式的變量定義可得機(jī)械臂的第i連桿的質(zhì)心Ci的轉(zhuǎn)動(dòng)慣量相對(duì)第i連桿MDH坐標(biāo)原點(diǎn)oi的轉(zhuǎn)動(dòng)慣量的轉(zhuǎn)化量為[LxxiLxyiLxziLxxiLyyiLyziLzzi]和質(zhì)心位置的轉(zhuǎn)化量為[lxilyilzi],mi為第i連桿的質(zhì)量。

2 動(dòng)力學(xué)參數(shù)分析

機(jī)械臂最小動(dòng)力學(xué)參數(shù)集是從標(biāo)準(zhǔn)動(dòng)力學(xué)參數(shù)中篩選組合出來的[10],由于并不是每一個(gè)標(biāo)準(zhǔn)動(dòng)力學(xué)參數(shù)都對(duì)機(jī)器人關(guān)節(jié)力矩有影響,并且機(jī)械臂動(dòng)力學(xué)參數(shù)對(duì)關(guān)節(jié)力矩的貢獻(xiàn)值也有大小之別。機(jī)械臂標(biāo)準(zhǔn)動(dòng)力學(xué)參數(shù)的可篩選組合性表明了其動(dòng)力學(xué)參數(shù)對(duì)應(yīng)回歸矩陣空間是可以進(jìn)行壓縮和精簡的[11]。

(7)

ψ∈Rm×n(m≥n)是非列滿秩的矩陣,ψ的最小線性無關(guān)列組對(duì)應(yīng)標(biāo)準(zhǔn)參數(shù)φ的基參數(shù)空間。最小線性無關(guān)列數(shù)也是ψ的空間維度k,而n-k列是需要?jiǎng)h除和線性整合的。ψ的空間維度k對(duì)應(yīng)標(biāo)準(zhǔn)參數(shù)φ的基參數(shù)個(gè)數(shù),此k列貢獻(xiàn)度最大。篩選出的n-k列對(duì)應(yīng)標(biāo)準(zhǔn)參數(shù)φ中貢獻(xiàn)度幾乎為0的參數(shù)和貢獻(xiàn)度是相關(guān)性的參數(shù),其中,貢獻(xiàn)度為0的參數(shù)需要?jiǎng)h除,貢獻(xiàn)度是相關(guān)性的參數(shù)則需要基于基參數(shù)進(jìn)行線性組合。

3 列變換矩陣的求取

標(biāo)準(zhǔn)的矩陣QR分解如下:

(8)

其中:Q為m×m的正交矩陣;R為m×n的上三角矩陣。

通過對(duì)ψ進(jìn)行分解求解動(dòng)力學(xué)參數(shù)貢獻(xiàn)度。由于多關(guān)節(jié)機(jī)械臂動(dòng)力學(xué)參數(shù)的回歸矩陣ψ是非列滿秩的矩陣,對(duì)其標(biāo)準(zhǔn)QR分解不具有唯一性。但是通過Householder方式的 QR分解能夠解決ψ的標(biāo)準(zhǔn)QR分解不唯一的問題[12]。于是通過列變換的QR分解使得ψ具有QR分解的唯一性,即

(9)

列變換矩陣Π∈Rn×n顯示ψ各列的獨(dú)立性與相關(guān)性,其中Ak、Bn-k為R的分塊矩陣。

r=rank(ψ),Q為正則矩陣,Ak為k×k的上三角非奇異矩陣,rank(ψΠ)=rank(Ak)。由Householder矩陣H1,…,Hk和變換矩陣Π1,…,Πk可得:

(H1,…,Hk)ψ(Π1,…,Πk)=R=

(10)

令C(m-k)(n-k)=[zk+1,…,zn],其中參數(shù)zj為C的各列,令‖zp‖2=max(‖zk+1‖2,…,‖zn‖2),其中j=k+1,…,n。如果rank(Ak)=k,那么‖zp‖2就等于0。于是,該變換矩陣Π的設(shè)計(jì)方法如下:

廣義矩陣ψ初次QR分解,對(duì)初次分解獲得R矩陣的對(duì)角線數(shù)值的非零值對(duì)應(yīng)的列位置和零值對(duì)應(yīng)的位置進(jìn)行標(biāo)記,根據(jù)標(biāo)記位置構(gòu)造列變換原矩陣Θ,然后對(duì)n階單位矩陣進(jìn)行Θ列變換來構(gòu)造列變換矩陣Π。

(11)

令E為m×n的單位矩陣,該隨機(jī)數(shù)下的廣義矩陣的旋轉(zhuǎn)矩陣p為:

p=ΘE

(12)

4 基參數(shù)及參數(shù)線性組合

利用(12) 式列變換矩陣p可得:

(13)

其中:φ1為基參數(shù);φ2為基于基參數(shù)的可線性組合的參數(shù)。

將ψ進(jìn)行列變換,可得k個(gè)獨(dú)立的列矩陣ψ1和n-k個(gè)需刪除、組合的列矩陣ψ2:

ψrandomp=[ψ1ψ2]

(14)

對(duì) (14) 式隨機(jī)廣義矩陣ψrandom進(jìn)行列變換后,對(duì)ψrandomp進(jìn)行QR分解,可得:

[Q1RkQ2Rk×(n-k)]

(15)

φmin=φ1+βφ2

(16)

最小參數(shù)集φmin及其對(duì)應(yīng)回歸矩陣ψ1與標(biāo)準(zhǔn)動(dòng)力學(xué)參數(shù)φ及其對(duì)應(yīng)回歸矩陣ψ關(guān)系如下:

ψφ=ψ1φmin

(17)

5 求解與驗(yàn)證

通過數(shù)值算例驗(yàn)證本文提出方法求解的最小動(dòng)力學(xué)參數(shù)集的正確性和有效性,機(jī)械臂的MDH參數(shù)見表1所列。

表1 機(jī)械臂MDH參數(shù)

應(yīng)用本方法得到的機(jī)械臂最小動(dòng)力學(xué)參數(shù)集的符號(hào)組合解如下:

Lzz1+Lyy2+Lyy3+0.487lz2+0.300 2lz3+

0.592 9m2+0.209(m3+m4+m5+m6),

Lxx2-Lyy2-0.186(m3-m4-m5-m6),

Lxy2,

Lxz2-0.431 8lz3+0.040 3(m3+m4+m5+m6),

Lyz2,

Lzz2+0.186(m3+m4+m5+m6),

lx2+0.431 8(m3+m4+m5+m6),

ly2,

Lxx3-Lyy3+Lyy4+0.866 2lz4+

0.187 1(m4+m5+m6),

Lxy3-0.020 3l4z-0.008 8(m4-m5-m6),

Lxz3,Lyz3,

Lzz3+Lyy4+0.866 3lz4+0.188(m4+m5+m6),

lx3-0.020 3(m4-m5-m6),

ly3-lz4-0.433 1(m4-m5-m6),

Lxx4-Lyy4+Lyy5,

Lxy4,Lxz4,Lyz4,

Lzz4+Lyy5,lx4,

ly4+lz5,

Lxx5-Lyy5+Lyy6,Lxy5,

Lyz5,Lzz5+Lyy6,lx5,

ly5-lz6,

Lxx6-Lyy6,Lxy6,Lxz6,

Lyz6,Lzz6,lx6,ly6。

因?yàn)槟Σ另?xiàng)和電機(jī)轉(zhuǎn)動(dòng)慣量項(xiàng)參數(shù)不能精簡,為了保證機(jī)械臂動(dòng)力學(xué)模型的完整性,本文將其納入建模之中,所以原有各連桿參數(shù)有[LxxiLxyiLxziLyyiLyziLzzilxilyilzimi],其中i=1,…,6,參數(shù)總量為60個(gè)。本文方法去除連桿1中的Lxx1,Lxy1,Lxz1,Lyy1,Lyz1,lx1,ly1,lz1和連桿6中的lz6,最后通過線性組合成上述的最小動(dòng)力學(xué)參數(shù)集的36個(gè)符號(hào)組合解結(jié)果,同時(shí)在CORE I5 CPU 上求解該六關(guān)節(jié)機(jī)械臂的最小參數(shù)集用時(shí)為3.64 s。

為驗(yàn)證該方法的最小參數(shù)集在動(dòng)力學(xué)上的準(zhǔn)確性與有效性,利用Matlab軟件里robot tools模塊建立該機(jī)械臂模型如圖3所示,設(shè)計(jì)的圓弧軌跡如圖4所示。

圖3 機(jī)械臂結(jié)構(gòu)

圖4 機(jī)械臂軌跡

該軟件模塊獲得圓弧軌跡下各關(guān)節(jié)運(yùn)動(dòng)變量,其中機(jī)械臂第3關(guān)節(jié)的角度、角速度和角加速度如圖5所示。對(duì)比簡化前與簡化后機(jī)械臂的各個(gè)關(guān)節(jié)的力矩值,如圖6所示。

圖5 機(jī)械臂的角位移、角速度和角加速度

圖6 模型簡化前與簡化后力矩對(duì)比

通過Adams建立的機(jī)械臂仿真模型,采用有限項(xiàng)傅里葉級(jí)數(shù)獲得辨識(shí)激勵(lì)軌跡,利用最小二乘法辨識(shí)六關(guān)節(jié)機(jī)械臂的前3個(gè)關(guān)節(jié)的最小動(dòng)力學(xué)參數(shù)集，并計(jì)算前3個(gè)關(guān)節(jié)的關(guān)節(jié)力矩,與Adams仿真出的該激勵(lì)軌跡下的前3個(gè)關(guān)節(jié)力矩對(duì)比如圖7所示;3個(gè)關(guān)節(jié)的仿真力矩與辨識(shí)力矩誤差值如圖8所示。

圖7 關(guān)節(jié)力矩Adams仿真與辨識(shí)計(jì)算力矩

圖8 關(guān)節(jié)力矩Adams仿真值與辨識(shí)計(jì)算值誤差值

通過以上2種方式對(duì)比,本文的最小動(dòng)力學(xué)參數(shù)集計(jì)算方法完整反映出原有機(jī)械臂動(dòng)力學(xué)特性;同時(shí)通過辨識(shí)結(jié)果分析可知,該方法計(jì)算的最小動(dòng)力學(xué)參數(shù)集有較好的辨識(shí)性。

6 結(jié) 論

為了求取多關(guān)節(jié)機(jī)械臂最小動(dòng)力學(xué)參數(shù)集,本文利用生成的隨機(jī)數(shù)模擬多關(guān)節(jié)機(jī)械臂運(yùn)動(dòng)參數(shù),設(shè)計(jì)一種回歸矩陣的變換矩陣,從而在較短時(shí)間內(nèi)確定最小動(dòng)力學(xué)參數(shù)集。以六關(guān)節(jié)機(jī)械臂為例,通過數(shù)值算例對(duì)所提出的方法進(jìn)行了驗(yàn)證,結(jié)果表明,本文提出的方法獲取的最小動(dòng)力學(xué)參數(shù)集準(zhǔn)確有效,所提出的最小動(dòng)力學(xué)參數(shù)計(jì)算方法穩(wěn)定、可靠。