基于Delaunay三角劃分網格計算多連桿機器人工作空間

吳 虎， 李鑫寧， 楊先海， 薛 鵬， 譚 帥

(1. 山東理工大學機械工程學院， 山東淄博 255000； 2.國家采油裝備工程技術研究中心， 山東東營 257091)

多連桿機器人是由多個活動關節(jié)將剛性桿件組成的連桿機構，具有工作空間范圍大、靈活性好、機械結構簡單、避障性能好、可避免機構奇異性等優(yōu)點，因此廣泛應用于工業(yè)生產、工程實踐、醫(yī)療手術等領域[1-2]。多連桿機器人為機械裝置的主要執(zhí)行部件，其靈活程度對于整個機構系統(tǒng)的流暢運行起著至關重要的作用。

工作空間的基本概念[3]及靈活性[4]、 邊界的確定[5]、靈活角的分析[6]、 工作空間的解析關系[7]等研究領域的探索， 為多連桿機器人工作空間的研究提供了理論基礎。 目前， 多連桿機器人工作空間的求解方法比較多， 主要可以歸納為圖解法(幾何法)、 解析法和數(shù)值積分法三大類。 圖解法[8]利用幾何繪圖的原理求得工作空間的邊界， 該方法直觀性強。 García-Murillo等[9]提出一種用于分析2(3-RRPS, R為轉動副, P為移動副, S為球副)并聯(lián)機器人工作空間的幾何方法， 通過計算機輔助設計(CAD)軟件進行實例分析計算。 劉懷周等[10]利用圖解法分析深水作業(yè)串聯(lián)機械臂的工作空間邊界曲線， 運用遺傳算法求得工作空間最優(yōu)解。 解析法[11]通過多次求解邊界曲面(線)的包絡問題， 確定工作空間的邊界方程。 Chen等[12]運用解析法確定混合刻板繪圖儀的工作空間邊界， 并采用優(yōu)化算法計算了工作空間的體積。 于紅英等[13]利用解析法建立了并聯(lián)機構各支鏈工作空間的邊界方程。 何世凱等[14]提出矩陣零空間的方法求解力旋量工作空間邊界。 利用數(shù)值積分法[15]計算多連桿機器人的工作空間時， 通過正運動學方程計算多連桿機器人末端執(zhí)行機構端點的坐標， 從而形成工作空間邊界曲線或曲面。 隨著計算機軟、 硬件技術的不斷發(fā)展， 數(shù)值積分法可以分析任意形式多連桿機器人的工作空間， 因此得到了越來越廣泛的應用， 其中比較典型的是基于蒙特卡羅法的工作空間求解。 與圖解法與解析法相比， 蒙特卡羅法省略了復雜的數(shù)學推導過程， 可編程性強， 容易實現(xiàn)可視化， 被廣泛地應用在各種領域機器人的工作空間求解中， 但是傳統(tǒng)的蒙特卡羅法計算工作空間不夠精確。 He等[16]提出了一種基于增加邊界點密度的蒙特卡羅法計算工作空間。 Zhao等[17]利用傳統(tǒng)的蒙特卡羅法生成種子工作空間， 根據(jù)正態(tài)分布擴展種子工作空間， 設置精度閾值動態(tài)調整標準差， 提出了一種體素算法來求出工作空間的體積， 通過不斷細化邊界部分來減小結果誤差。 Bader等[18]混合計算自運動流形范圍估計工作空間包絡和蒙特卡羅積分估計方向體積的算法， 針對高冗余機器人提出了一種容錯工作空間的計算方法。 Hamida等[4]將遺傳算法和蒙特卡羅法相結合， 針對機器人結構緊湊和靈巧性最大化進行優(yōu)化。 樊桂菊等[19]建立了以可達工作空間性能和結構緊湊為指標的優(yōu)化模型， 利用遺傳算法求最優(yōu)參數(shù)， 采用網格劃分法求解工作空間體積， 達到較高的精度。 雖然上述方法在計算精度、 邊界確定等方面取得了很大的進步， 但同時也存在直觀性不強、 表達式繁瑣、 對于多自由度的機器人運算量過大等問題。

為了更準確地描述多連桿機器人的工作空間，在工作空間密度函數(shù)[20]和自適應非結構化網格生成理論[21-22]的基礎上，本文中提出基于Delaunay三角劃分網格的多連桿機器人工作空間計算方法(簡稱本文方法)，通過與數(shù)值積分法、 網格劃分法進行對比， 驗證本文方法的準確性。

1 靈活工作空間與工作空間

可達工作空間是在參考坐標系中末端執(zhí)行機構所能達到的所有位置坐標的集合，記作W；靈活工作空間是在總的可達工作空間內，末端執(zhí)行機構可以以任意姿態(tài)達到的坐標點所構成的工作空間，記作W1；次工作空間是工作空間中去除靈活工作空間所剩余的空間，記作W2。三者的關系為

W=W1+W2

。

(1)

圖1所示為平面3R機器人的運動模型。 該機構由長度分別為l1、l2、l3的連桿組成， 令l1>l2+l3，l2>l3， 機械臂驅動轉角θi∈[-π,π](i=1, 2, 3)，以此例來說明靈活工作空間與工作空間的關系。

li(i=1,2,3)—連桿長度； θi—機械臂驅動轉角； xoy—平面直角坐標系。圖1 平面3R(R為轉動副)機器人的運動模型

多連桿機器人的工作空間與靈活工作空間如圖2所示。圓C1半徑R1=l1-l2-l3，圓C2半徑R2=l1-(l2-l3)=l1-l2+l3，圓C3半徑R3=l1+(l2-l3)=l1+l2-l3，圓C4半徑R4=l1+l2+l3。顯然，工作空間W是圓C1和C4圍成的圓環(huán)區(qū)域，靈活工作空間W1是圓C2和C3圍成的圓環(huán)區(qū)域，其他剩余的工作空間為W2。

li(i=1,2,3)—連桿長度； Ch(h=1,2,3,4)—圓心為坐標原點o的不同半徑的圓。圖2 多連桿機器人的工作空間與靈活工作空間

2 工作空間密度函數(shù)

將機械機構的工作空間均勻地劃分成等尺寸的單元，每個單元內末端執(zhí)行機構可以到達的點數(shù)與整個工作空間內到達點數(shù)比值除以單元的面積或體積即為工作空間密度[23]。工作密度函數(shù)值越大，證明機械機構的靈活性越大， 同時可以更加直觀地表達多連桿機器人的工作空間。 蒙特卡羅法屬于實驗數(shù)學的一個分支， 利用隨機數(shù)進行統(tǒng)計實驗， 將均值、 概率等作為待解決問題的數(shù)值解， 也稱為蒙特卡羅模擬。 通過研究多連桿機器人的工作空間密度， 采用蒙特卡羅法，從概率統(tǒng)計的角度借助計算機技術進行仿真， 可以直觀地描述工作空間的形狀、 特征及性質。

2.1 蒙特卡羅法計算步驟

工作空間解算的蒙特卡羅算法原理[6]是

W={M(k)∶k∈G}

，

(2)

式中：k為廣義關節(jié)變量；M(k)為正運動學參數(shù)；G為關節(jié)空間集合。

采用蒙特卡羅法求解多連桿機器人工作空間的主要思想如下：通過計算機仿真隨機產生大量的關節(jié)空間點，映射至工作空間，得到一定數(shù)量的末端位置隨機點，進而統(tǒng)計解算工作空間的邊界。蒙特卡羅算法的步驟如下：

1)根據(jù)靈活工作空間的定義，對各個關節(jié)轉角的樣本進行劃分取點，構建樣本空間G；

2)依據(jù)正運動學參數(shù)M(k)解算機械臂末端執(zhí)行機構所能到達的各個位置；

3)輸出機械臂末端執(zhí)行機構的位置，構建靈活工作空間W，樣本容量越大，所得到的幾何輪廓越清晰，靈活工作空間W越精確。

平面3R機器人正運動計算模型如圖3所示。根據(jù)圖中多連桿機器人末端執(zhí)行機構位置與連桿長度及關節(jié)轉角的幾何關系，用Denavit-Hartenberg(D-H)法列出其正運動學方程。

li(i=1,2,3)—連桿長度； θi—機械臂驅動角； P—末端執(zhí)行機構位置點； o、 A、 B—關節(jié)中心點； x0oy0、 x1oy1、 x2Ay2、 x2By2—平面直角坐標系。圖3 平面3R(R為轉動副)機器人正運動計算模型

假設坐標軸z0、z1、z2垂直于紙面向外，從空間中的點(x0,y0,z0)到(x1,y1,z1)的齊次旋轉變換矩陣為

(3)

從點(x1,y1,z1)到(x2,y2,z2)的齊次旋轉變換矩陣為

(4)

從點(x2,y2,z2)到(x3,y3,z3)的齊次旋轉變換矩陣為

(5)

從點(x3,y3,z3)到(x0,y0,z0)的齊次旋轉變換矩陣為

(6)

可得末端執(zhí)行機構P點的位置矢量為

(7)

式中(xP,yP,zP)為末端執(zhí)行機構位置點P點的坐標。

基于以上公式推導，可以得到兩桿機器人的坐標如下：

(8)

(9)

式中l(wèi)i為第i個連桿的長度。

以上運動學方程是在直角坐標系下構建的，要得到極坐標系下對應的表達式，可利用直角坐標與極坐標的關系，即

(10)

(11)

式中：r為極徑；φ為極角。

2.2 密度函數(shù)的構建

蒙特卡羅法計算多連桿機器人工作空間采用的是統(tǒng)計學的方法。 Sigmoid函數(shù)是統(tǒng)計學中常用的數(shù)學模型。 根據(jù)工作空間密度函數(shù)符合隨機密度函數(shù)[20]的特性， 通過建立數(shù)學模型來描述工作空間密度。

設N連桿機械臂桿件的密度函數(shù)為fN(p)，其中p∈G，N連桿機械臂的工作空間密度函數(shù)為

f1,2,…,N(p)=(f1*f2*…*fN)(p)

,

(12)

式中*為在G上的卷積運算符。

ρ1,2,…,N(x)=(ρ1⊙ρ2⊙…⊙ρN)(x)

,

(13)

工作空間是連續(xù)的空間，其特點是隨機變量可以連續(xù)地取值，所有的連續(xù)隨機變量都有概率密度函數(shù)。u為多維向量，μ為均值，σ2為方差，多元正態(tài)分布密度函數(shù)的形式為

(14)

對式(14)進行嚴格泛化，得

(15)

式中：Σ為N×N型對稱正定矩陣； |·|為矩陣的行列式。μ、Σ的求解如下：

(16)

(17)

(18)

連續(xù)的概率密度函數(shù)不局限于1，因此在工作空間的研究中，所有連續(xù)隨機變量都是可測的，所有連續(xù)分布都有密度。對于N連桿機械臂，當每個連桿轉角有M個狀態(tài)時，連桿機械臂工作空間內的樣本個數(shù)為K=MN，此樣本空間的方差S1,2,…,N(K)為

(19)

式中XK為空間樣本方差的矩陣。

當M足夠大時，實際工作空間密度函數(shù)的方差為

(20)

式中In為n×n型單位矩陣。

根據(jù)卷積迭代及極值定理，Sigmoid函數(shù)構建的工作空間密度函數(shù)為

exp[-(u-μ)T(u-μ)Σ-1] 。

(21)

設桿長l1=1.8，l2=2.4，l3=3.1，在轉角θ1、θ2、θ3不受限制的情況下，取樣本點個數(shù)M=8 000，根據(jù)式(9)運算并繪制出工作空間密度圖，如圖4所示。

2.3 計算仿真

根據(jù)構建的密度函數(shù)，建立五自由度多連桿機器人模型，如圖5所示，其中l(wèi)1=0.42，l2=0.375，l3=0.42，l1、l2、l3與底座轉軸不在一個平面內(l1、l2、l3與底座轉軸沿轉軸同方向偏移0.1、 0.2、 0.3)，轉角θ1、θ2、θ3以及底座轉角θ4、 執(zhí)行機構轉角θ5不受限制， 所取的樣本點個數(shù)M=20 000。 圖6所示為五自由度多連桿機器人工作空間仿真視圖，4個視圖全面描述了工作空間密度。

xy—坐標平面。圖4 3R(R為轉動副)機器人工作空間密度云圖

xyz—空間直角坐標系。圖5 五自由度多連桿機器人模型

(a)工作空間三維視圖

(b)xy坐標平面

(c)xz坐標平面

(d)yz坐標平面圖6 五自由度多連桿機器人工作空間仿真視圖

六自由度多連桿機器人模型如圖7所示，其中l(wèi)1=0.42，l2=0.375，l3=0.42，l4=0.42， 桿l1、l2、l3、l4與底座轉軸均在一個平面內，轉角θ1、θ2、θ3、θ4以及底座轉角θ5， 執(zhí)行機構轉角θ6不受限制，所取的樣本點個數(shù)M=20 000。六自由度多連桿機器人工作空間仿真視圖如圖8所示。

圖7 六自由度多連桿機器人模型

3 工作空間計算

根據(jù)連桿機器人運動學及工作空間密度研究，提出基于Delaunay三角劃分網格計算多連桿機器人工作空間的方法， 原理如圖9所示。 本文方法是一種非結構化網格生成方法，首先建立一個包含所有點的初始三角網，然后將剩余的點逐點插入，并且確保初始三角網成為Delaunay三角網，同時，每個Delaunay三角形的外接圓不包含其他結點；若出現(xiàn)外接圓內有其他結點，則需修改局部網格，直至滿足每個Delaunay三角形的外接圓不包含其他結點這一條件。

(a)工作空間三維視圖

(b)xy坐標平面

(c)xz坐標平面

(d)yz坐標平面圖8 六自由度多連桿機器人工作空間仿真視圖

圖9 基于Delaunay三角劃分網格的 多連桿機器人工作空間計算方法原理

1)建立初始網格。首先找到包含該密度空間點集的外接矩形，通常保證外接矩形的邊長是點集最大邊界長度的3倍。然后任意連接矩形的對角線，形成2個三角形，作為初始Delaunay三角網格。

2)點的逐次插入。在保證計算區(qū)域內工作空間密度函數(shù)均勻的情況下，插入點必在三角形內或三角形邊上，3種插入點的情況如圖10所示。插入點P′與三角形的3個頂點U、V、W相連，三角形及其外接圓形成Delaunay腔，然后連接插入點與空腔的每個頂點，每次插入一個新點之后，重復連接插入點與空腔的每個頂點操作。

U、 V、 W—三角形的3個頂點； V1—三角形的鄰近頂點； P′—插入點。圖10 點的逐次插入中插入點的3種情況

3)構造外接圓。首先對U、V、W這3個點構造外接圓T1，把所有結點包含進去，之后每次只引入一個結點，生成相應的圓T2，T3，……，直至所有結點都進入Delaunay三角劃分網格。

4)刪除最近邊。將三角形中離新結點最近的一條邊刪除，對網格進行重生成，于是新結點與老結點形成新的三角形，如圖11所示。

圖11 通過刪除最近邊對網格進行重生成所得新三角形

6)工作空間面積的求和。工作空間面積S的求解公式為

(22)

式中：n1為內部空間子網格個數(shù)；n2為外部邊界空間子網格個數(shù)。

以圖3為模型，機器人連桿參數(shù)如表1所示。

表1 平面3R(R為轉動副)機器人連桿長度與機械臂驅動轉角

用數(shù)值積分法、 網格劃分法和本文方法分別計算工作空間面積， 并與理論計算值進行比較， 如表2所示。 從表中可以看出， 本文方法的計算結果較準確， 與其他2種方法相比， 誤差率分別降低了0.99%和0.50%， 提高了工作空間的計算精度。

表2 不同方法計算的工作空間面積比較

4 結論

為了更準確地描述多連桿機器人的工作空間，在對多連桿機器人進行運動學分析的基礎上，本文中提出基于Delaunay三角劃分網格計算連桿機器人工作空間，得出如下結論：

1)應用構建的Sigmoid工作空間密度函數(shù)，不僅可以直觀地表達工作空間的工作范圍，而且可以定量地得出工作空間體積的數(shù)值，避免了繁瑣的表達式及過量的運算，為工程中的實際應用提供了有效的計算方法，同時為工作空間后續(xù)的深入研究提供了思路。

2)通過本文方法與數(shù)值積分法、 網格劃分發(fā)法等其他計算方法的對比，驗證了本文方法的準確性。

3)在工作空間的研究中，正運動學的D-H參數(shù)算法隨著連桿機器人自由度的增加，需要在每個關節(jié)增加坐標系，計算量隨之增大。如何減少計算量、 提高計算效率是需要進一步研究的問題。