繩驅準連續型機械臂的動力學研究

黃玉平，賈龍飛，郭亞星，陶云飛，鄭繼貴

（1.北京精密機電設備研究所，北京，100076；2.航天伺服驅動與傳動技術實驗室，北京，100076）

0 引 言

在諸多領域中，定期的監測和維護是確保設備安全工作的重要方式，但由于此類設備具有復雜的結構，使得可供監測及維護作業的空間非常狹小，維修人員不易進行操作。為了順利完成狹小空間中設備的監測和維護，機械臂被用來代替人類進行作業。此類機械臂需要具備穿越狹窄環境、回避障礙物的能力，并具有足夠大的靈巧操作空間。

準連續型機械臂具有多個自由度，可呈現多種運動模式，可在狹小、復雜及非結構化的環境中執行任務[1]，且具有多目標優化能力，廣泛應用于航天[2]、航空[3]、特種加工、核電站監測[4]、醫療[5]和災害救援[6]等領域。由于準連續型機械臂采用繩索驅動的方式，具有多個自由度，而且驅動繩和臂桿之間存在摩擦力，為實現準連續型機械臂的精確控制，故需建立準確的動力學模型。

常用的動力學建模方法包括牛頓歐拉法[7]、Schiehlen 法[8]、Kine 法、拉格朗日法[9]等，其中拉格朗日法借助廣義坐標，可極大簡化建模過程，因此在實際工程中應用廣泛。但如果模型中的自由度超過6 個或更多，則會大大增加拉格朗日方法的運算量，在實際應用中無法滿足實時控制的要求。遞歸牛頓歐拉法在多個自由度的問題上，具有較小的運算量，故廣泛用于建立準連續型機械臂的動力學建模。

針對準連續型機械臂的動力學研究，Yuan[10]建立了一種連續型機械臂的靜力學模型，驗證了在機械臂偏轉角度較大的情況下，摩擦力對機械臂的運動精度影響較大。Rone[11]利用虛功率原理建立了繩驅連續型機械臂的靜態動力學摩擦模型，在不考慮慣性效應的情況下推導了靜力平衡模型。馬曙光[12,13]利用空間算子代數方法推導了機械臂的動力學方程，對動力學模型中的繩索模型進行改進，通過增大仿真時間步長臨界值來提高動力學模型的數值仿真效率，可實現實時仿真，但未給出具體的關節力矩與驅動繩拉力的轉換關系，并未對最優化問題展開分析。劉天亮[14]針對繩驅動超冗余機械臂建立了動力學模型，并給出力矩與驅動繩拉力的轉換關系，但建立的動力學模型較為復雜，且利用該方法得到的驅動繩拉力變化曲線不光滑。

準連續型機械臂具有多個自由度，因此本文采用遞歸牛頓歐拉法進行動力學建模分析。通過確定機械臂運動過程中，各變量之間的轉換關系，在考慮摩擦力的前提下，建立了完整的動力學模型。并針對關節力矩與驅動繩拉力之間的轉換，優化了最優化模型中的目標函數，采用置信域算法解決了24 根驅動繩拉力值曲線的不光滑問題，并通過降元進一步縮減了計算時間。

1 動力學建模

1.1 機械結構

本文所研究的準連續型機械臂如圖1a 所示，主要包括機械臂本體結構、驅動機構、導軌。機械臂本體結構由8 個剛性臂桿串聯而成，臂桿之間通過十字軸萬向節連接，驅動機構由伺服電機、聯軸器、滾珠絲杠、滑塊組成。伺服電機帶動絲杠轉動來控制滑塊移動，從而拉動驅動繩，進一步可通過控制驅動繩產生的位移及拉力來調整機械臂的姿態。

每一組驅動機構控制一根驅動繩，與驅動繩數量相同的24組驅動機構以圓周陣列的形式均布于驅動機構中的基座上。穿過孔ija、孔ijb、孔ijc 的3 根驅動繩驅動第j（j=1,2,…,8 且j≥i）個臂桿，在第i個臂桿的圓形截面上取驅動繩穿過孔ija、孔ijb、孔ijc 的圓心與該圓形截面圓心連成的直線與Zi軸之間的夾角為該點的位置角度φija、φijb、φijc。如圖1c 中的孔41a、孔47a、孔43c 的位置角度分別為0、90、270 。

圖1 準連續型機械臂結構示意Fig.1 Schematic Diagram of Quais-continuous Manipulator

續圖1

1.2 運動學分析

準連續型機械臂模型基于以下假設：

a）每節臂桿為剛體；

b）驅動繩不可伸長，且質量和慣量很小。

圖2 為基座的本體坐標系{O0}以及8 個臂桿的本體坐標系{O1}～{O8}。本文建立相對于地面固定不變的世界坐標系{Oxyz}，當8 個臂桿的軸線在同一水平方向時，以8 個臂桿的軸線方向為X軸方向，以重力的反方向為Z軸方向，Y軸方向可根據右手定則確定。

圖2 機械臂的D-H 坐標系Fig.2 The D-H Coordinate System of Manipulator

第i個臂桿具有2 個自由度，分別是相對于第i-1個臂桿的偏航角αi和俯仰角θi，坐標系{Oi-1}先沿著Xi-1軸平移L，再繞Zi-1軸旋轉αi，最后繞Yi軸旋轉θi轉換為坐標系{Oi}。可得到轉換矩陣為

在轉換矩陣的基礎上可推導出各驅動繩長度、關節角度、臂桿末端點等數據之間的轉換關系，即可建立準連續型機械臂完整的運動學模型。

1.3 動力學分析

1.3.1 單臂桿受力分析

下面在第i個臂桿的本體坐標系{Oi}中，對其進行受力分析，得到如圖3 所示的單臂桿受力分析圖。

圖3 單臂桿受力分析Fig.3 Force Analysis of the Single Arm

圖3 中的力和力矩都為第i個臂桿本體坐標系{Oi}中的力和力矩，根據力平衡及力矩平衡，可得：

圖4 單臂桿簡化受力分析示意及局部放大圖Fig.4 Simplified Force Analysis of the Single Arm and Local Enlarged View

合并式（2）和式（3），可得式（4）。

基座運動的力由直線導軌提供，用符號F0表示。由于有末端外力Fe及自身重量Gi的影響，需要24 根驅動繩提供驅動力，進而轉換為驅動力矩，來維持整個系統的平衡。fiia、fiib、fiic為第i個臂桿的3 根驅動繩的拉力，合力為fii。fija、fijb、fijc為第i個關節處驅動第j個臂桿的3 根驅動繩的拉力，合力為fij。為驅動第i+1 至第8 個臂桿的繩子在第i個關節處的合力。通過第i個關節的3×(9-i)根驅動繩在第i個關節處的合力為fi。以上各力可通過式（8）進行轉換。niiy、niiz為驅動繩拉力fii等效至Yi軸、Zi軸的力矩，為等效至Yi軸、Zi軸的力矩。niix、niiy、niiz的合力矩為nii，的合力矩為，通過第i個關節的所有驅動繩在第i個關節處的合力fi平移至坐標系{Oi}原點產生的合力矩。以上各力矩可通過式（9）進行轉換。

1.3.2 數學模型建立

牛頓歐拉法是一種基于牛頓方程及歐拉方程分析矢量力學的方法，可采用遞歸方法計算準連續型機械臂的動力學參數。牛頓歐拉法的遞推由兩部分組成：第1 部分為正向遞推，從基座到第1 個臂桿，然后到第2 個臂桿，依次往后，直至第8 個臂桿，由式（10）可計算每個臂桿的角速度、角加速度及質心加速度；第2 部分為反向遞推，從第8 個臂桿到基座，由式（11）反向計算作用力和力矩及關節驅動力矩。本文需要將關節驅動力矩等效轉換為驅動繩的驅動力矩。

式中mi和Ici分別為第i個臂桿的的質量和質心的慣性張量。

綜上可得準連續型機械臂各參數轉化示意，如圖5所示。圖5 中包括運動學與動力學2 大模塊，運動學模塊主要分為3 個小模塊，分別為求解旋轉矩陣及平移矩陣、求解驅動繩在每個位置的方向以及求解驅動繩與臂桿軸線的夾角及驅動繩相對于臂桿的速度。動力學模塊主要分為4 個小模塊：利用牛頓方程求解質心加速度等參數、利用歐拉方程求解合力矩等參數、利用最優化算法將二力矩轉換為三繩拉力、利用摩擦力公式計算繩子在不同位置處的拉力大小。

圖5 準連續型機械臂各參數轉化示意Fig.5 Schematic of Conversion of Parameters

2 動力學仿真

在建立動力學模型的基礎上，借助Matlab 進行準連續型機械臂的動力學仿真分析，得出準連續型機械臂中24 根驅動繩的拉力大小與其運動狀態的關系。圖6 為準連續型機械臂的運動軌跡。為保證對應的軌跡規劃下各臂桿的位置、速度及加速度變化曲線光滑且連續，令機械臂由水平狀態開始，圍繞x0軸線做圓周運動，0～3 s 由水平位置過渡到“畫圓”的位置，3～15 s 每個臂桿的末端點圍繞一個圓轉動，15～18 s 由“畫圓”的位置過渡到水平位置。

圖6 準連續機械臂的運動軌跡Fig.6 The Trajectory of Manipulator

規定第i個臂桿的起始點為(xi，yi，zi)，坐標的變化規律如式（13）所示。其中T=8 s 為臂桿末端繞圓轉一周的時間，采用七次多項式求解g(t)，得出式（14）。

根據圖6 對應的運動軌跡及前面建立的動力學模型，仿真可得8 個關節處的所有驅動繩的等效力矩在Yi軸、Zi軸方向上的分量大小niy、niz，如圖7 所示。初始狀態下等效力矩在Yi軸上的分量為niy(0)，從第1～8個關節分別為-64mgl、-49mgl、-36mgl、-25mgl、-16mgl、-9mgl、-4mgl、-1mgl（其中，m為臂桿的質量；g為重力加速度；l為臂桿的長度）。從圖7 中可以看出當臂桿末端繞圓形軌跡轉動一周時，等效力矩變化了 3個周期，等效力矩曲線光滑且連續。

圖7 各個關節處的等效力矩變化Fig.7 Equivalent Torque Variation of Each Joint

3 非光滑時變最優化問題研究

通過動力學模型可推導出各個關節處所需的等效力矩，初始狀態下在第i個關節處的所有驅動繩的等效力矩在Yi軸上的分量niy(0)=-(9-i)2mgl，進而可推出在第i個關節處的驅動第i個臂桿的驅動繩的等效力矩在Yi軸上的分量niiy(0)=(2i-17)mgl。如何將關節兩個方向的力矩niiy、niiz轉換為3 根驅動繩上的拉力尤為重要。

由2 個方程求解3 個未知量可有無窮解，首先建立對應的最優化問題的數學模型，如式（15）所示。

上述問題的目標函數J1(f)為變量fiia、fiib、fiic的線性函數，約束條件也為變量fiia、fiib、fiic的線性函數，可利用Matlab優化工具箱中的linprog函數求解該線性規劃問題，通過仿真可得到24 根驅動繩拉力隨時間的變化曲線，以驅動第3 個臂桿的繩子為例，3 根驅動繩拉力隨時間的變化曲線如圖8 所示。

圖8 線性規劃下的驅動繩拉力隨時間的變化曲線Fig.8 The Variation Curve of the Tension of Driving Cables with Time under Linear Programming

從圖8 中可以看出驅動繩拉力變化曲線連續但有些曲線不光滑、導數不連續，始終有1 根驅動繩的拉力約為零。若按照圖8 中的變化規律控制驅動繩拉力，則容易拉斷驅動繩，影響整個系統的運動平穩性。針對上述問題，將線性目標函數J1(f)改為如式（16）所示的非線性目標函數J2(f)，并設定了如式（17）所示的的fiia，fiib，fiic的下限。

非線性規劃問題的求解方法有許多種，本文分別采用遺傳算法、模擬退火算法、粒子群算法、置信域算法進行非線性規劃問題的求解，圖9 為4 種算法下的驅動第3 個臂桿的3 根驅動繩拉力在9～10 s 的局部放大圖。

圖9 4 種算法下的驅動第3 個臂桿繩拉力的局部放大示意Fig.9 Local Enlarged View of the Tension of Cables Driving the Third Link under Four Algorithms

從圖9 中可以看出，若利用前3 種算法進行最優化求解得到的曲線相對不光滑且計算時間較長，在某些地方容易陷入局部最小值，而且前3 種算法計算速度較慢且易出錯，故不適用于計算24 根驅動繩拉力的實時期望數值。

故本文采用置信域算法進行求解非線性規劃問題，可將計算時間大大縮短，通過仿真可得出24 根驅動繩拉力隨時間的變化曲線，以驅動第3 個臂桿的3根驅動繩為例，其拉力變化曲線如圖10 所示。

圖10 非線性規劃下的驅動繩拉力隨時間的變化Fig.10 The Variation Curve of the Tension of Driving Cables with Time under Nonlinear Programming

分析仿真結果發現，任何時刻在第i個關節處驅動第i個臂桿的3 根驅動繩的拉力的等效力矩nii都等于本關節所有驅動繩的拉力的等效力矩ni減去驅動第i+1至第8 個臂桿的繩子在第i個關節處的等效合力矩，從而驗證了動力學模型的正確性。

對比圖8、圖10，可以看出雖然線性規劃下的驅動繩拉力的最大值（例如圖8 中max(f33)=21.14mgl/r）小于非線性規劃下的驅動繩拉力的最大值（例如圖10 中max(f33)=38.93mgl/r），但非線性規劃下的24 根驅動繩的拉力變化曲線連續且光滑，在機械臂的實時運動過程中，可提高整個系統的平穩性，有效地降低抖動。

該優化問題中含有3 個變量，為進一步降低計算量、減少計算時間，可利用式（7）中的兩力矩niiy、niiz與3 根驅動繩拉力fiia、fiib、fiic的等式關系，將3 個變量轉化為1 個變量，再采用置信域算法針對非線性目標函數求取最優解。最終得到的仿真結果與圖10相同，但仿真時間從108 s 縮短至12 s，大大提高了兩力矩與3 繩拉力之間的最優化轉換速度。

4 結束語

本文利用牛頓歐拉法建立了準連續型機械臂的動力學模型，精準分析了運動過程中24 根驅動繩在不同位置的方向、速度，以及摩擦力對繩子拉力大小的影響，將每個關節處的驅動繩拉力及支持力兩個未知量合并為一個未知量進行求解，從而實現了機械臂各參數之間的轉化，并通過仿真驗證了動力學模型的正確性。

同時針對關節力矩與驅動繩拉力之間的轉換關系，提出了一種非線性目標函數，采用置信域算法消除了驅動繩拉力曲線的不光滑問題，提高了整個系統的平穩性，有效降低了抖動，并通過降元進一步縮減了計算時間，可滿足準連續型機械臂的實時控制需求。

后續工作中將在本文建立的動力學模型的基礎上，在準連續性機械臂的控制模型中增加驅動繩拉力的反饋信號，采用PID 或其他方法實時調節24 根驅動繩拉力，從而提高機械臂運動的平穩性及精確性。