工業(yè)機器人廣泛應用于工業(yè)生產(chǎn)過(guò)程中的搬運、裝配、涂裝等方面,并且正向著(zhù)高速、高精度、高智能化的方向發(fā)展,加之產(chǎn)業(yè)需求的不斷擴大與提升,諸如新型的激光焊接、激光切割、多機器人協(xié)作等任務(wù)對機器人的控制精度提出了更高要求[1,2]1-2]。設計基于機器人動(dòng)力學(xué)模型的控制器是實(shí)現高速高精度運動(dòng)控制的有效途徑[3]3],但該類(lèi)控制器的設計需要已知機器人的動(dòng)力學(xué)參數。一般情況下,機器人的動(dòng)力學(xué)參數難以直接測量,需采用特定的辨識實(shí)驗獲得[4]4]。Atkeson等[5]5]選擇關(guān)節空間的五次多項式作為激勵軌跡,采用小二乘法作為參數估計的方法。Swevers等[6]6]首次提出了采用基于傅里葉級數的激勵軌跡進(jìn)行辨識實(shí)驗,該方法產(chǎn)生較大影響,在隨后的辨識研究中,基本均采用傅里葉級數作為激勵軌跡[7]7]。傅里葉級數是指由不同幅值和周期的三角函數構成的非線(xiàn)性函數,該函數具有連續的多階導數,其二階導數具有較大幅值,能充分激發(fā)系統的動(dòng)態(tài)特性;并且通過(guò)合理設置參數,能保證函數起始和結束時(shí)刻具有相同值,便于機器人重復執行激勵軌跡。但傳統傅里葉級數形式的激勵軌跡無(wú)法保證各關(guān)節的角速度和角加速度在軌跡的起始和停止時(shí)刻為零,導致機器人在執行激勵軌跡時(shí)存在抖動(dòng)。文獻[8]在此基礎上提出了一種基于改進(jìn)的傅里葉級數的辨識方法,該方法考慮速度、加速度邊界條件,保證激勵軌跡起始階段的穩定性,提高了辨識精度。目前,確定激勵軌跡參數有兩種應用較多的方法,一種是基于回歸矩陣的條件數小化方法[9]9];另一種是基于對數行列式的優(yōu)化方法[10]10]。參數估計的方法也有很多,其中應用較多的兩種方法為小二乘法[11]11]和大似然估計法[8]8]。其他的方法包括擴展卡爾曼濾波法[12]12],加權小二乘法[12]12],非線(xiàn)性小二乘法[13]13]等。在機器人控制方面,Slotine等做了大量深入的研究,提出了直接自適應、間接自適應、復合自適應等控制算法[14,15,16,17]14-17];針對機器人工作環(huán)境限制、負載變化等問(wèn)題,提出將自適應理論結合阻抗控制,并證明了控制系統的穩定性[18]18];焦曉紅等[19]19]將自適應結合魯棒控制,針對誤差擾動(dòng)上界確定和不確定兩種情況提出控制律,并通過(guò)仿真證明了控制算法的可行性。文獻[20]提出了自適應神經(jīng)網(wǎng)絡(luò )控制算法。利用徑向基神經(jīng)網(wǎng)絡(luò )來(lái)調整滑?？刂圃鲆?消除滑?？刂戚斎攵秳?dòng)并優(yōu)化系統性能。

筆者以ZCR07S六自由度工業(yè)機器人為實(shí)驗研究對象,基于Codesys平臺完成機器人動(dòng)力學(xué)參數辨識實(shí)驗與自適應控制實(shí)驗。辨識出一組動(dòng)力學(xué)參數并實(shí)現了自適應控制算法。

n連桿剛性機器人動(dòng)力學(xué)方程可由二階非線(xiàn)性微分方程[21]21]描述為

H(q)q¨+C(q,q˙)q˙+G(q)+τf=τ$\mathit{{\rm H}}(\mathit{q})\ddot{\mathit{q}}+\mathit{C}(\mathit{q},\dot{\mathit{q}})\dot{\mathit{q}}+\mathit{G}(\mathit{q})+\mathit{\tau }{}_f=\mathit{\tau }$(1)

式中:q,q˙$\dot{\mathit{q}}$,q¨$\ddot{\mathit{q}}$∈Rn×1分別代表關(guān)節角位置、角速度和角加速度向量;τ∈Rn×1為關(guān)節驅動(dòng)力矩;H(q)∈Rn×n為慣性矩陣;C(q,q˙$\dot{\mathit{q}}$)∈Rn×1為哥氏力和離心力矩陣;G(q)∈Rn×1為重力項;τf∈Rn×1為摩擦力項,本研究中摩擦力模型為τf=fvq˙+fcsgn(q˙)$\mathit{\tau }{}_f=\mathit{f}{}_v\dot{\mathit{q}}+\mathit{f}{}_c\mathrm{sgn}(\dot{\mathit{q}})$,其中fv,fc分別為黏性摩擦系數與庫倫摩擦系數,sgn代表符號函數。

機器人系統的動(dòng)力學(xué)結構特性如下:

特性1 可適當定義C(q,q˙),H˙?2C$\mathit{C}(\mathit{q},\dot{\mathit{q}}),\dot{\mathit{{\rm H}}}-2\mathit{C}$為斜對稱(chēng)矩陣,即對任意x∈Rn均有

xT[H˙?2C]x=0$\mathit{x}{}^{\text{Τ}}[\dot{\mathit{{\rm H}}}-2\mathit{C}]\mathit{x}=0$ (2)

特性2 存在一個(gè)依賴(lài)于機械臂參數的參數向量a,使得H(q),C(q,q˙$\dot{\mathit{q}}$),G(q),τf滿(mǎn)足線(xiàn)性關(guān)系

H(q)q¨+C(q,q˙)q˙+G(q)+τf=Y(q,q˙,q¨)a?????????(3)$\begin{array}{l}\mathit{{\rm H}}(\mathit{q})\ddot{\mathit{q}}+\mathit{C}(\mathit{q},\dot{\mathit{q}})\dot{\mathit{q}}+\mathit{G}(\mathit{q})+\\ \mathit{\tau }{}_f=\mathit{Y}(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})\mathit{a}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}(3)\end{array}$

式中:Y(q,q˙$\dot{\mathit{q}}$,q¨$\ddot{\mathit{q}}$)為n×m的已知關(guān)節變量的回歸矩陣,它是機器人廣義坐標及其各階導數的已知函數矩陣,a為m×k描述機器人質(zhì)量特性的未知定常參數向量。即機器人的動(dòng)力學(xué)特性是機器人慣性參數的線(xiàn)性函數。

式(3)中Y的有些列恒為零,有些列間存在常值線(xiàn)性關(guān)系,所以不是所有慣性參數都對力矩有影響,通過(guò)線(xiàn)性關(guān)系重組標準慣性參數[11,22]或數值解析方法[23]獲得可辨識的小慣性參數集amin。式(3)可重寫(xiě)為

H(q)q¨+C(q,q˙)q˙+G(q)+τf=Y(q,q˙,q¨)a=Y′(q,q˙,q¨)amin?????????(4)$\begin{array}{l}\mathit{{\rm H}}(\mathit{q})\ddot{\mathit{q}}+\mathit{C}(\mathit{q},\dot{\mathit{q}})\dot{\mathit{q}}+\mathit{G}(\mathit{q})+\mathit{\tau }{}_f=\\ \mathit{Y}(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})\mathit{a}={\mathit{Y}}^{\prime }(\mathit{q},\dot{\mathit{q}},\ddot{\mathit{q}})\mathit{a}{}_{\min }\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}(4)\end{array}$

式中Y′(q,q˙$\dot{\mathit{q}}$,q¨$\ddot{\mathit{q}}$)各列線(xiàn)性無(wú)關(guān)。六自由度串聯(lián)機器人小慣性參數如表1所示,慣性參數重組關(guān)系為

Izz1r=(m3+m4+m5+m6)(d22${}_2^2$+d23${}_3^2$+r23${}_3^2$)+

Ia1+d22${}_2^2$m2+2mz3r3+Iyy2+Iyy3+Izz1

Ixx2r=-(m3+m4+m5+m6)·

(d23${}_3^2$+r22${}_2^2$-r23${}_3^2$)+Ixx2-Iyy2

Ixz2r=Ixz2-d3(m3+m4+m5+m6+mz3)r3

Izz2r=Ia2+d23${}_3^2$(m3+m4+m5+m6)+Izz2

mx2r=d3(m3+m4+m5+m6)+mx2

Ixx3r=2mz4r4+(m4+m5+m6)r24${}_4^2$+

Ixx3-Iyy3+Iyy4

Izz3r=2mz4r4+(m4+m5+m6)r24${}_4^2$+

Iyy4-Izz3

my3r=my3-mz4-(m4+m5+m6)r4

Ixx4r=Ixx4-Iyy4+Iyy5

Izz4r=Iyy5+Izz4

my4r=Imy4+Imz5

Ixx5r=Ixx5-Iyy5+Iyy6

Izz5r=Iyy6+Izz5

my5r=Imy5-Imz6

Ixx6r=Ixx6-Iyy6

為實(shí)現基于動(dòng)力學(xué)模型的控制,需要對結構已知的模型參數進(jìn)行辨識,獲取機器人的一組小慣性參數值,即式(4)中的amin。機器人的參數辨識是一個(gè)復雜的系統工程,包括運動(dòng)學(xué)與動(dòng)力學(xué)建模、激勵軌跡的選取、實(shí)驗設計、數據采集與處理、參數估計、模型驗證等環(huán)節,其流程圖如圖1所示。

ZCR07S機器人6 個(gè)關(guān)節均為旋轉關(guān)節,本體如圖2所示,運動(dòng)學(xué)模型的建立采用文獻[11]中的方法,連桿坐標系如圖3所示,MDH參數如表2所示。

動(dòng)力學(xué)模型由SYMORO+[24]軟件推導,以自定義符號表達式形式在robot-name.dyn文件中給出,由于篇幅限制,文中不再給出。

各關(guān)節激勵軌跡為N次諧波的正弦和余弦函數有限項和,激勵軌跡關(guān)節角位置表達式為

qi(t)=∑l=1Nalωflsin(ωflt)?blωflcos(ωflt)+qi0?????????(5)$\begin{array}{l}q{}_i(t)={\displaystyle \sum _{l=1}^{{\rm N}}\frac{a{}_l}{\omega {}_{f}l}}\sin (\omega {}_{f}lt)-\\ \frac{b{}_l}{\omega {}_{f}l}\cos (\omega {}_{f}lt)+q{}_{i0}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}(5)\end{array}$

關(guān)節角速度、角加速度分別為一階、二階導數。實(shí)驗中:基頻為0.1 Hz;周期為10 s;N為5;qi0為位置偏移量;al,bl確定了三角函數的幅值。

確定激勵軌跡參數al,bl,qi0的問(wèn)題可描述為非線(xiàn)性約束優(yōu)化問(wèn)題:

優(yōu)化目標為

q*(t)=argmin(J) (6)

約束條件為

式中:J為待確定的目標函數;qmax,q˙$\dot{q}$max,q¨$\ddot{q}$max分別為關(guān)節限制大位置、速度、加速度。實(shí)驗所用機器人的關(guān)節位置、速度限位如表3所示。

式(7)中后面的3 個(gè)約束條件保證了激勵軌跡起始和終止位置、速度、加速度為零。根據Hadamard不等式,正定矩陣的行列式小于等于其對角項的乘積,可計算det(YTY)上界,可小化矩陣條件數。目標函數J[25]為

J=1∏j=1n∑i=1nY2ij$\mathit{J}=\frac{1}{{\displaystyle \prod _{j=1}^n{\displaystyle \sum _{i=1}^n\mathit{Y}}}{}_{ij}^2}$

式中Y為m×n的回歸矩陣,Yij為矩陣Y的第i行第j列元素。通過(guò)MatlabFunction函數計算目標函數參數,將目標函數與式(7)中的約束條件代入到Matlab fmincon函數計算激勵軌跡參數。激勵軌跡如圖4所示。

辨識實(shí)驗設計與數據采集通過(guò)Codesys平臺完成,同時(shí)采集關(guān)節角位置、速度和力矩,采樣頻率為250 Hz。由Codesys采集到的力矩經(jīng)過(guò)如下轉換關(guān)系換算為機械臂端力矩:機械臂端力矩=采集的力矩相對值×電機額定力矩×減速器減速比/1 000。實(shí)驗所用電機額定力矩和減速器減速比如表4所示。

實(shí)驗中采集的關(guān)節位置軌跡光滑,無(wú)需濾波處理,關(guān)節速度、加速度通過(guò)中心差分算法計算得到,并通過(guò)Matlab smooth函數對關(guān)節加速度、力矩進(jìn)行平滑處理。

小慣性參數值的估計采用標準小二乘法,去除掉一些相對標準偏差較大(RSD>30%)的慣性參數后,得到29 個(gè)描述6 個(gè)關(guān)節機器人動(dòng)力學(xué)模型的小慣性參數,如表5所示。

通過(guò)比較實(shí)測力矩與由辨識模型所估計的力矩的均方根誤差來(lái)判斷所辨識的模型是否符合精確度要求。圖5為重新生成的一組驗證軌跡。

圖6,7分別代表激勵軌跡和驗證軌跡下的實(shí)測力矩與估計力矩。

表6為激勵軌跡和驗證軌跡下的實(shí)測力矩與估計力矩的均方根誤差,誤差值較小,辨識的模型可用。

選擇控制律和自適應律[14]為

τ=H?(q)q¨r+C?(q,q˙)q˙r+G?(q)+F?(q˙)?KDs?????????(8)$\mathit{\tau }=\widehat{\mathit{{\rm H}}}(\mathit{q})\ddot{\mathit{q}}{}_r+\widehat{\mathit{C}}(\mathit{q},\dot{\mathit{q}})\dot{\mathit{q}}{}_r+\widehat{\mathit{G}}(\mathit{q})+\widehat{\mathit{F}}(\dot{\mathit{q}})-\mathit{{\rm K}}{}_D\mathit{s}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}\text{?}(8)$

a?˙=?Γ?1YT(q,q˙,q˙r,q¨r)s$\dot{\widehat{\mathit{a}}}=-\mathit{\Gamma }{}^{-1}\mathit{Y}{}^{\text{Τ}}(\mathit{q},\dot{\mathit{q}},\dot{\mathit{q}}{}_r,\ddot{\mathit{q}}{}_r)\mathit{s}$ (9)

式中:H?$\widehat{\mathit{{\rm H}}}$(q),C?$\widehat{\mathit{C}}$(q,q˙$\dot{\mathit{q}}$),G?$\widehat{\mathit{G}}$(q),F?$\widehat{\mathit{F}}$(q˙$\dot{\mathit{q}}$)分別為慣性矩陣、離心力和哥氏力、重力、摩擦力的估計值,其中f?$\widehat{\mathit{f}}$v,f?$\widehat{\mathit{f}}$c分別為辨識的黏性摩擦系數與庫倫摩擦系數;KD,Γ分別為PD控制器增益、自適應增益,均選為正定對角常數矩陣;q˙$\dot{\mathit{q}}$r,q¨$\ddot{\mathit{q}}$r分別為參考速度、參考加速度,q˙r=q˙??Λq?,q¨r=q¨??Λq?˙,$\dot{\mathit{q}}{}_r=\dot{\mathit{q}}{}^{*}-\mathit{\Lambda }\stackrel{?}{\mathit{q}},\ddot{\mathit{q}}{}_r=\ddot{\mathit{q}}{}^{*}-\mathit{\Lambda }\dot{\stackrel{?}{\mathit{q}}},$其中q˙?,q¨?$\dot{\mathit{q}}{}^{*},\ddot{\mathit{q}}{}^{*}$分別為期望速度、期望加速度;q,q˙$\dot{\mathit{q}}$分別為實(shí)際位置、實(shí)際速度,其中q?=q?q?,q?˙=q˙?q˙?$\stackrel{?}{\mathit{q}}=\mathit{q}-\mathit{q}{}^{*},\dot{\stackrel{?}{\mathit{q}}}=\dot{\mathit{q}}-\dot{\mathit{q}}{}^{*}$;s為位置誤差和速度誤差的加權和,其表達式為s=q˙?q˙r=q?˙+Λq?$\mathit{s}=\dot{\mathit{q}}-\dot{\mathit{q}}{}_r=\dot{\stackrel{?}{\mathit{q}}}+\mathit{\Lambda }\stackrel{?}{\mathit{q}}$。由式(4)知:控制律式(8)可寫(xiě)為

τ=Y(q,q˙,q˙r,q¨r)a?min?KDs$\mathit{\tau }=\mathit{Y}(\mathit{q},\dot{\mathit{q}},\dot{\mathit{q}}{}_r,\ddot{\mathit{q}}{}_r)\widehat{\mathit{a}}{}_{\min }-\mathit{{\rm K}}{}_D\mathit{s}$ (10)

式中:a?$\widehat{\mathit{a}}$min為由辨識得到的小慣性參數;回歸矩陣Y可直接由符號計算軟件SYMORO+得出,其中矩陣元素Yij以自定義形式在文件robot-name.dim中給出,對應于robot-name.dim中的DGiKj。則由得到的小慣性參數a?$\widehat{\mathit{a}}$min與回歸矩陣Y即可實(shí)現自適應控制器。自適應控制器結構如圖8所示。

研究中使用的六軸工業(yè)機器人控制系統采用一主多從的控制模式,主站采用ARM+Linux控制方案,從站采用帶EtherCAT接口的高創(chuàng )CDHD系列伺服驅動(dòng),主站與從站、從站與從站之間的通訊采用工業(yè)以太網(wǎng)EtherCAT技術(shù),軟件平臺以Codesys作為開(kāi)發(fā)環(huán)境,以ST和CFC作為主要編程語(yǔ)言實(shí)現機器人動(dòng)力學(xué)控制,機器人控制系統結構如圖9所示。

實(shí)驗的期望軌跡為

?????????????????????????????????q?1=5t2?109t3?????0≤t<3q?2=3t2?23t3???0≤t<3q?3=5t2?109t3???0≤t<3q?4=?10t2+209t3???0≤t<3q?5=10t2?209t3???0≤t<3q?6=10t2?209t3???0≤t<3q?i=q?i(3)????3≤t,i=1,2,?,6$\{\begin{array}{l}q{}_1^{*}=5t{}^2-\frac{10}{9}t{}^3\text{?}\text{?}\text{?}\text{?}\text{?}0\le t<3\\ q{}_2^{*}=3t{}^2-\frac{2}{3}t{}^3\text{?}\text{?}\text{?}0\le t<3\\ q{}_3^{*}=5t{}^2-\frac{10}{9}t{}^3\text{?}\text{?}\text{?}0\le t<3\\ q{}_4^{*}=-10t{}^2+\frac{20}{9}t{}^3\text{?}\text{?}\text{?}0\le t<3\\ q{}_5^{*}=10t{}^2-\frac{20}{9}t{}^3\text{?}\text{?}\text{?}0\le t<3\\ q{}_6^{*}=10t{}^2-\frac{20}{9}t{}^3\text{?}\text{?}\text{?}0\le t<3\\ q{}_i^{*}=q{}_i^{*}(3)\text{?}\text{?}\text{?}\text{?}3\le t,i=1,2,?,6\end{array}$

程序運行采樣周期為4 ms,自適應增益均為1,在實(shí)際的調節過(guò)程中,根據系統的穩定性與跟蹤性能之間的權衡,確定KD=diag[400,350,130,22,9,2],Λ=diag[10,10,10,30,30,25](KP=ΛKD)。實(shí)驗結果如圖10所示,圖10中曲線(xiàn)由上至下分別代表期望位置y1,實(shí)際位置y2,期望速度y3,實(shí)際速度y4,位置誤差y5,速度誤差y6。跟蹤過(guò)程中,關(guān)節1至關(guān)節6大位置誤差分別為-0.06°,0.22°,0.34°,-0.29°,0.45°,0.85°。關(guān)節1至關(guān)節6位置定位誤差分別為-0.06°,0.02°,0.04°,-0.11°,-0.11°,0.09°。

采用一種計算效率高的優(yōu)化準則設計辨識實(shí)驗的激勵軌跡,通過(guò)Codesys平臺完成辨識實(shí)驗,去掉相對標準偏差較大的慣性參數后,得到29 個(gè)描述機器人動(dòng)力學(xué)特性的小慣性參數,并通過(guò)驗證軌跡驗證了辨識參數的精確性;接著(zhù)基于所辨識的小慣性參數和機器人動(dòng)力學(xué)模型,將一種自適應控制算法應用于六自由度工業(yè)機器人軌跡跟蹤控制,通過(guò)實(shí)驗驗證了所研究方法的有效性。主要體現為應用創(chuàng )新:第一個(gè)創(chuàng )新點(diǎn)是辨識出了ZCR07S工業(yè)機器人動(dòng)力學(xué)參數;第二個(gè)創(chuàng )新點(diǎn)是將自適應控制算法應用于ZCR07S工業(yè)機器人,并取得了良好的控制效果。