一種改進的SCARA機器人動力學參數辨識方法

嚴 浩 白瑞林 吉 峰

1.江南大學輕工過程先進控制教育部重點實驗室，無錫，2141222.無錫信捷電氣股份有限公司，無錫，214072

一種改進的SCARA機器人動力學參數辨識方法

嚴 浩1白瑞林1吉 峰2

1.江南大學輕工過程先進控制教育部重點實驗室，無錫，2141222.無錫信捷電氣股份有限公司，無錫，214072

針對SCARA機器人動力學參數辨識問題，提出了一種基于優化改進傅里葉級數的辨識方法。根據SCARA機器人完整動力學方程，推導得到動力學模型的線性形式。采用改進傅里葉級數作為機器人關節的激勵軌跡，使得關節角度滿足連續周期性，并且關節角速度和角加速度在軌跡起始和停止時刻為零。為進一步提高辨識精度，以SCARA機器人觀測矩陣的條件數為目標函數，采用基于排擠機制的小生境遺傳算法對激勵軌跡的系數進行優化?？紤]到測量噪聲的影響，采用加權最小二乘法(WLS)作為參數估計方法。實驗結果表明，采用所提方法能準確辨識出SCARA機器人的動力學參數，兩關節力矩測量值和預測值的殘差均方根分別減小了11.50%和26.35%。

機器人；動力學；小生境遺傳算法；激勵軌跡；參數辨識

0 引言

隨著工業機器人應用于生命醫療、激光焊接、汽車電子等高精密領域[1-2]，機器人技術向高速高精度方向發展。目前，國產工業機器人大部分仍采用傳統PID控制策略，難以實現高精度控制。設計基于機器人動力學模型的控制策略是實現高速高精度運動控制的有效方法[3-4]，但該類控制策略要以機器人的精準動力學參數為基礎。

機器人動力學參數的獲取方法可分為解體測量法、CAD法以及整體辨識法。機器人結構復雜，很多參數難以通過直接測量得到。采用計算機建模的CAD法忽略了機器人的裝配誤差和材料的分布特性；整體辨識法能夠考慮機器人實際工作中的各種因素影響，因而受到了廣泛關注[5-6]。

機器人動力學參數辨識包括6個步驟：動力學建模、模型線性化、激勵軌跡優化、數據采樣與處理、參數估計方法、試驗驗證[7]。ATKESON等[8]構造5次多項式作為關節空間的激勵軌跡，并采用最小二乘法作為動力學參數估計方法。SWEVER等[9]提出基于周期性傅里葉級數的辨識方法，使得機器人可以循環執行激勵軌跡。但傳統傅里葉級數無法保證各關節的角速度和角加速度在軌跡的起始和停止時刻為零，導致機器人在執行激勵軌跡時存在抖動。同時由于采樣誤差的存在，由未經優化的傅里葉級數得到的觀測矩陣會放大采樣誤差的影響[10]。丁亞東等[11]針對多關節機器人的動力學參數辨識，提出一種分步辨識方法，降低了辨識方程復雜性。陳恩偉等[12]探討了人工神經網絡在動力學參數辨識中的一般規律，使神經網絡的結構與權值具有明確物理意義，但僅分析了機器人末端連桿慣性參數的辨識問題。

為了提高SCARA機器人的動力學參數辨識精度，本文對傳統傅里葉級數進行改進，使得關節角速度和角加速度在軌跡起始和停止時刻為零。同時，以觀測矩陣的最小條件數作為目標函數，采用基于排擠機制的小生境遺傳算法對傅里葉級數的系數進行優化，降低采樣誤差對辨識精度的影響。為進一步降低采樣誤差的影響，采用加權最小二乘法對動力學參數進行辨識。

1 動力學模型及其線性化

對于n關節的機器人系統，考慮關節摩擦的影響，其動力學方程可以用二階非線性微分方程來表示[13]：

(1)

SCARA機器人第3、第4關節共用一個連桿，第2、3、4關節的伺服電機、同步帶、滾珠絲杠、升降臺等都安裝在關節2的連桿上，因此，關節1、關節2決定SCARA機器人的平面定位精度和軌跡跟蹤精度，關節3決定豎直方向上的運動。由于關節3與前兩關節垂直，所以前兩關節的運動不會對關節3產生影響。由此，僅對關節1和關節2進行動力學建模及動力學參數辨識。式(1)的標量表達式如下：

(2)

根據SCARA機器人結構，利用拉格朗日法對SCARA機器人進行動力學建模[14]，得到式(2)中的各項參數：

式中，m1、m2分別為連桿1和連桿2的質量；l1、l2分別為連桿1和連桿2的長度；(x1,y1,z1)、(x2,y2,z2)分別為前兩連桿的質心在坐標系1和坐標系2中的位置坐標；IZz1、IZz2分別為連桿1和連桿2沿Z軸的轉動慣量；fci和fvi分別表示第i關節的庫侖摩擦因數和黏性摩擦因數。

文獻[15]已經證明存在一依賴于機器人參數的參數向量，使得機器人動力學方程可以線性化，即

(3)

根據式(2)的SCARA機器人動力學方程具體形式，得到線性化方程的顯式表達如下：

(4)

式中各項參數如下：

機器人動力學參數辨識實驗過程為：機器人各關節跟蹤特定的激勵軌跡，通過采樣或間接計算得到關節的輸入力矩、關節角度、角速度、角加速度，估計動力學參數φ。機器人各關節跟蹤激勵軌跡時，在采樣時刻t1,t2,…,tN對關節輸入力矩等進行N次直接或間接采樣，則可以得到如下方程：

YNφ=τN

(5)

即

(6)

2 激勵軌跡優化設計

采用周期性的傅里葉級數作為激勵軌跡，能充分激勵機器人的動力學特性，并且多次采樣取均值可以提高信噪比。但在軌跡的起始和停止時刻，傳統傅里葉級數不能保證關節角速度和角加速度為零，使機器人無法平滑運行。

對于n階非齊次方程組YNφ=τN，當由N組觀測數據構成的觀測矩陣YN和N組測量得到的輸出力矩τN各有擾動δYN和δτN時，則添加擾動后的非齊次方程組為

(YN+δYN)(φ+δφ)=(τN+δτN)

(7)

由此解得

(8)

因YNφ=τN，則有‖τN‖≤‖YN‖‖φ‖，則式(8)可轉化為

(9)

則動力學參數φ解的相對誤差為

(10)

(11)

2.1改進的傅里葉級數

傳統的M階傅里葉級數定義如下：

(12)

式中，i為第i個關節；ωf為傅里葉級數的基頻，每個關節的基頻相同；ai,k、bi,k為傅里葉級數的系數；θ0i為偏移量；M為傅里葉級數的階數，決定軌跡的帶寬；ai,k、bi,k、θ0i為自由系數。

相應的關節運行速度和加速度為

(13)

機器人動力學參數辨識過程中，為保證運行的連續周期性，以及起始和結束時刻的平穩性，需要激勵軌跡滿足如下約束條件：

(14)

t0=0ttf=2π/(2πωf)

式中，t0、ttf分別為周期運動的起始和結束時刻；θi_init為第i關節的起始和結束時刻角度。

傳統的傅里葉級數難以保證上述約束條件，因此，對傳統傅里葉級數的激勵軌跡進行改進，將常數項θ0i用5次多項式來代替，得到改進傅里葉級數的激勵軌跡：

(15)

解得5次多項式的系數為

2.2激勵軌跡優化

確定激勵軌跡的形式之后，以一定頻率對改進傅里葉級數進行采樣，得到對應的關節角度、角速度和角加速度，共采集N組數據，從而構造觀測矩陣YN。以觀測矩陣YN的條件數作為目標函數，即適應度函數，為增強全局尋優能力，采用基于排擠機制的小生境遺傳算法對傅里葉級數的系數進行優化，獲得最小條件數時對應的系數ai,k、bi,k。

基于排擠機制的小生境遺傳算法尋優步驟如下：

(1)選擇編碼策略。取M=5，改進的傅里葉級數共有20個待優化參數，因此，將20個參數進行實數編碼，構成可行解的集合。

(2)確定適應度函數。根據式(6)計算得到觀測矩陣YN，并以其條件數作為目標函數，即適應度函數。

(3)建立排擠機制。記第i個個體為pi，第j個個體為pj，則個體pi和pj之間的歐氏距離為‖pi-pj‖，取一較小整數ξ，f(·)表示個體的適應度值。在對每一代進行選擇、交叉、變異操作之前，比較種群中任意兩個個體的相似性，若‖pi-pj‖lt;ξ，則表明pi和pj的相似程度較大。通過對適應度較小者施加一較強的懲罰函數，使其適應度變得極小，從而排擠掉其中適應度較小者，即若f(pi)gt;f(pj)，則令f(pi)=δf(pj)(δ為一很小的正數)，則在以后的進化中pj會以極大的概率被淘汰掉。

(4)確定遺傳算法參數。共進化400代，設置種群規模為50，交叉概率為0.4，變異概率為0.1。

(5)選擇操作。種群中個體的選擇操作采用輪盤賭法，并且將種群中個體適應度值取倒數，使得在輪盤中適應度值越小的個體在輪盤中被選中的概率越大。選擇操作控制進化方向，使得下一代種群優于上一代種群。

(6)交叉操作。以一定概率對種群中的個體進行染色體交叉操作，被選中染色體的交叉位置相同。

(7)變異操作。對種群中染色體以一定概率進行變異操作，得到新的種群個體。

(8)判斷是否已完成預定的進化代數，達到則停止進化。

采用改進的傅里葉級數作為激勵軌跡，并用對改進傅里葉級數的20個參數進行尋優，優化得到的傅里葉級數參數如表1所示。觀測矩陣YN的條件數隨進化代數的變化過程如圖1所示。

表1 優化后的傅里葉級數參數

圖1 觀測矩陣條件數Fig.1 Condition number of the observation matrix

得到的SCARA機器人兩關節位置輸入如圖2所示，對應的兩關節角速度和角加速度分別如圖3、圖4所示。從圖中可以看出，各關節的起始和結束時刻位置相同，保證了軌跡運行的周期性，并且起始和結束時刻的關節角速度和角加速度為零，保證了機器人平滑運行。

圖2 激勵軌跡角度Fig.2 Angle of the excitation trajectory

圖3 激勵軌跡角速度Fig.3 Angular velocity of the excitation trajectory

圖4 激勵軌跡角加速度Fig.4 Angular acceleration of the excitation trajectory

3 采樣數據處理與參數估計

3.1采樣數據處理

采用優化后的改進傅里葉級數作為機器人各關節激勵軌跡，運行時對關節的輸入力矩、關節角度和角速度進行采樣。由于采樣得到的數據易受測量噪聲的影響，所以多次重復執行激勵軌跡，對采樣得到的關節角度、角速度和關節力矩取均值以提高信噪比：

為避免速度測量噪聲給加速度計算帶來的影響，文獻[9]根據采樣得到的各關節角度，反向擬合成傅里葉級數形式，再微分得到各關節角速度和角加速度。目前伺服電機一般采用高精度編碼器，根據關節角度求角速度也經過濾波處理，得到的角速度信號比較準確，所以本文僅將采樣得到的關節角速度擬合成傅里葉級數形式，關節速度的傅里葉級數

采用線性最小二乘法得到參數ai,k、bi,k，再對關節速度的傅里葉級數進行微分，得到關節角加速度信號，避免了兩次微分帶來的傳遞誤差。

3.2加權最小二乘估計

機器人在重復跟蹤激勵軌跡時，在N個時間節點t1,t2,…,tN，通過高精度光電編碼器對機器人關節角度和角速度進行采集，由于電機輸出力矩不易直接測量，所以根據電機反饋輸入電流和電機轉矩系數計算得到關節輸入力矩。最終得到式(5)中的觀測矩陣YN和關節力矩向量τN。

式(5)的求解可采用加權最小二乘法[16]。在數據的采樣中，關節角度和角速度是通過高分辨率的光電編碼器經過濾波處理后測得的，相比關節力矩的測量噪聲可以忽略不計。由力矩測量噪聲的方差構成加權矩陣，從而區分準確數據和不準確數據，以此進一步提高辨識精度。此時，動力學參數φ的加權最小二乘估計

(16)

力矩測量噪聲方差通過下式進行估算：

式中，Ν為一個周期內采樣的點數；τij(k)為第j周期內，第i關節在第k次采樣得到的力矩輸出。

4 實驗設計與結果分析

辨識實驗在4自由度SCARA機器人的前兩關節進行，如圖5所示。機器人各個關節采用永磁同步伺服電機驅動，關節位置通過伺服電機的17位光電編碼器測量，速度通過對位置信號濾波后計算得到，力矩信號通過伺服電機的電樞反饋電流Iq和電機轉矩系數計算得到。機器人各個關節由伺服驅動器進行獨立PID控制，在NI實時系統上通過NI板卡進行激勵軌跡的位置給定。

圖5 SCARA機器人前兩關節Fig.5 The first two joints of SCARA robot

改進的激勵軌跡為5階傅里葉級數，即M=5，則兩個關節共有20個參數需要進行優化，優化后的關節激勵軌跡角度、角速度和角加速度分別如圖2、圖3和圖4所示。傅里葉級數基頻ωf=0.1 Hz，軌跡周期為10 s，重復執行激勵軌跡周期數L=5。

采用文獻[9]中提出的方法和本文提出的方法分別進行辨識，得到辨識結果如表2所示。

表2 辨識結果

機器人動力學參數的實際值是無法提前測得的，不能將兩種方法辨識得到的值與實際值進行對比。所以，選取一條不同于激勵軌跡的可執行軌跡作為驗證軌跡，使機器人跟蹤驗證軌跡，同時對關節輸入力矩信號進行采樣，再根據辨識出的動力學參數φ和期望的運動軌跡，用式(3)計算出關節力矩的預測值。為評價參數辨識的精度，采用力矩測量值和預測值的殘差均方根εRMS作為評價指標：

式中，τm(k)、τpred(k)分別為跟蹤驗證軌跡時，第k個采樣點處的采樣力矩和力矩預測值；K為每次運行的采樣點總數。

為滿足速度和加速度的邊界條件，仍選取式(15)所示形式的軌跡作為驗證軌跡，如圖6所示。

圖6 SCARA機器人關節驗證軌跡Fig.6 The joint validation trajectory of the SCARA robot

將上述驗證軌跡作為機器人兩關節的位置給定，根據兩種方法辨識得到的動力學參數，計算得到相應的各關節預測力矩輸出。同時采集各個關節的實際輸出力矩，得到兩關節測量力矩和預測力矩如圖7所示。

(a)關節1

(b)關節2圖7 關節力矩測量值和預測值Fig.7 The measured values and predicted values ofthe joint torque

根據式(16)計算得到兩種方法對應的各關節力矩殘差均方根，如表3所示?？梢钥闯觯簩Ω道锶~級數進行改進，使得角速度和角加速度在軌跡起始和停止時刻為零，并且使用基于排擠機制的小生境遺傳算法優化觀測矩陣條件數的方法，可以明顯提高辨識精度。采用所提方法，兩關節力矩測量值和預測值的殘差均方根分別降低了11.50%和26.35%。

表3 兩關節力矩殘差均方根

5 結論

根據SCARA機器人的完整動力學方程，推導得到動力學模型的線性形式。采用改進的傅里葉級數作為機器人關節的激勵軌跡，使得關節角度滿足連續周期性，并且關節角速度和角加速度在軌跡起始和停止時刻為零。為提高辨識精度，采用基于排擠機制的小生境遺傳算法對激勵軌跡的系數進行優化，以最小化SCARA機器人觀測矩陣的條件數。考慮到測量噪聲的影響，采用加權最小二乘法(WLS)作為參數估計方法。實驗結果表明，采用所提方法能準確辨識出SCARA機器人的動力學參數，并且相比于文獻[9]，兩關節力矩測量值和預測值的殘差均方根分別降低了11.50%和26.35%。

[1] 王田苗, 陶永. 我國工業機器人技術現狀與產業化發展戰略 [J]. 機械工程學報, 2014, 50(9): 1-13.

WANG Tianmiao, TAO Yong. Research Status and Industrialization Development Strategy of Chinese Industrial Robot [J]. Journal of Mechanical Engineering, 2014, 50(9): 1-13.

[2] BROGARDH T. Robot Control Overview: an Industrial Perspective [J]. Modeling Identification and Control, 2009, 30(3): 167-180.

[3] WU J, WANG J, YOU Z. An Overview of Dynamic Parameter Identification of Robots [J]. Robotics and Computer-Integrated Manufacturing, 2010, 26(5): 414-419.

[4] BONA B, INDRI M, SMALDONE N. Rapid Prototyping of a Model-based Control with Friction Compensation for a Direct-drive Robot [J]. IEEE/ASME Transactions on Mechatronics, 2006, 11(5): 576-584.

[5] 耿令波. 工業機器人動力學參數辨識方法研究 [D]. 南京：南京航空航天大學, 2013.

GENG Lingbo. Study on Method of Dynamic Parameter Identification of Industrial Robots [D].Nanjing: Nanjing University of Aeronautics and Astronautics, 2013.

[6] BENIMELI F, MATA V, VALERO F. A Comparison between Direct and Indirect Dynamic Parameter Identification Methods in Industrial Robots [J]. Robotica, 2006, 24(5): 579-590.

[7] QIN Z, BARON L, BIRGLEN L. A New Approach to the Dynamic Parameter Identification of Robotic Manipulators [J]. Robotica, 2010, 28(4): 539-547.

[8] ATKESON C G, AN C H, HOLLERBACH J M. Estimation of Inertial Parameters of Manipulator Loads and Links [J]. International Journal of Robotics Research, 1986, 5(3): 101-119.

[10] 黎柏春, 王振宇, DEMIN A,等. 一種改進的機器人動力學參數辨識方法 [J]. 中國工程機械學報, 2015, 13(5): 381-387.

LI Baichun, WANG Zhenyu,DEMIN A, et al. Improved Parametric Identification Method Based on Robot Dynamics [J]. Chinese Journal of Construction Machinery, 2015, 13(5): 381-387.

[11] 丁亞東, 陳柏, 吳洪濤,等. 一種工業機器人動力學參數的辨識方法 [J]. 華南理工大學學報: 自然科學版, 2015, 43(3): 49-56.

DING Yadong, CHEN Bai, WU Hongtao, et al. An Identification Method of Industrial Robot’s Dynamic Parameters [J]. Journal of South China University of Technology (Natural Science Edition), 2015,43(3): 49-56.

[12] 陳恩偉, 劉正士, 干方建. 機器人末端臂慣性參數辨識的人工神經網絡方法 [J]. 中國機械工程, 2006, 17(3): 268-271.

CHEN Enwei, LIU Zhengshi, GAN Fangjian. Application of ANN Method in Inertial Parameter Identification of the End Effector of Robot [J]. China Mechanical Engineering, 2006, 17(3): 268-271.

[13] 蔡自興.機器人學[M].北京：清華大學出版社，2005: 87-94.

CAI Zixing. Robotics [M]. Beijing:Tsinghua University Press, 2005: 87-94.

[14] 陸金偉, 李躍民. SCARA機器人的完整動力學顯方程的推導[J]. 振動、測試與診斷, 1995(3): 23-29.

LU Jinwei, LI Yuemin. The Derivation of the Complete Dynamic Equation of the SCARA Robot [J]. Journal of Vibration and Measurement and Diagnosis, 1995(3): 23-29.

[15] LIN S K. Minimal Linear Combinations of the Inertia Parameters of a Manipulator [J]. IEEE Transactions on Robotics and Automation, 1995, 11(3): 360-373.

[16] CALANCA A, CAPISANI L M, FERRARA A, et al. MIMO Closed Loop Identification of an Industrial Robot [J]. IEEETransactions on Control Systems Technology, 2011, 19(5): 1214-1224.

(編輯王旻玥)

AnImprovedDynamicsParameterIdentificationMethodforSCARARobots

YAN Hao1BAI Ruilin1JI Feng2

1.Key Laboratory of Advanced Process Control for Light Industry,Jiangnan University,Wuxi,Jiangsu,214122 2.Xinjie Electronic Co.,Ltd.，Wuxi，Jiangsu，214072

In order to identify the dynamics parameters of SCARA robots, a new identification method was proposed based on improved Fourier series. According to the complete dynamics equations of the SCARA robots, a linear form of the dynamics model was derived. The improved Fourier series was used as the excitation trajectory of the robot joints, and the joint angles satisfied the continuous periodicity, and the angular velocities and angular accelerations were as zero at the beginning and end of the trajectory. In order to further improve the identification accuracy, the condition numbers of the observation matrix of the SCARA robots were used as the objective function, and based on crowding mechanism the niche genetic algorithm was used to optimize the coefficients of the trajectory. Considering the influences of measurement noises, the weighted least square method (WLS) was used as the parameter estimation method. The experimental results show that the proposed method may accurately identify the dynamics parameters of the SCARA robots, torque measurements and predictive values of two joints, the residual root mean square values are reduced by 11.50% and 26.35% respectively.

robot; dynamics; niche genetic algorithm; excitation trajectory; parameter identification

TP242

10.3969/j.issn.1004-132X.2017.22.010

2016-12-08

江蘇省產學研前瞻性聯合研究項目(BY2015019-38)；江蘇高校優勢學科建設工程項目(PAPD)

嚴浩，男，1992年生。江南大學物聯網工程學院碩士研究生。主要研究方向為檢測技術與自動化裝置、機器人控制技術。E-mail:jiangnan_yh@163.com。白瑞林，男，1955年生。江南大學物聯網工程學院教授、博士研究生導師。吉峰，男，1981年生。無錫信捷電氣股份有限公司工程師。