辛軍煒, 鄧援超, 宋志成, 蔡宇飛
(湖北工業大學 機械工程學院,武漢 430068,E-mail:yiyuansdiy@163.com)
平面五桿機構是桿數最少的二自由度機構[1],由于運動副都為低副,其承載能力大,因此在工程上應用廣泛。相對于四桿機構,五桿機構具有較高的柔性,其中RRPRR型五桿機構的可動性不受桿長條件的約束[2],末端工作點具有較高的柔性,可實現較復雜的運動軌跡。
在連桿軌跡生成的研究中,由于傅里葉描述符對軌跡有很好的描述,且軌跡曲線的平移、縮放、旋轉對曲線的特征參數沒有影響,通過曲線特征參數擬合的方法得到了廣泛應用。主流的方法有通過數值圖譜法將曲線的特征參數存儲在計算機中后使用模糊識別算法進行匹配,以及采用擬合算法通過曲線的特征參數進行連桿參數擬合。褚金奎等[3]通過諧波特征參數法,對連桿機構的尺寸,機構安裝位置,連桿上工作點的位置參數進行綜合,建立了數值圖譜庫;馮立艷等[4]通過建立機構實際尺寸和安裝參數與傅里葉特征諧波參數之間的關系,結合遺傳算法對五桿機構的特征尺寸進行優化擬合;李學剛等[5]通過將五桿機構進行拆分,提出采用傅里葉級數法對平面五桿機構軌跡代數求解的方法;丁健生等[6]基于最小二乘法,利用傅里葉級數對位置點近似逼近,實現了少位置離散點特征提取的方法;Uesaka等[7]提出了一種將曲線總曲率的復指數函數通過擴展的傅里葉系數展開,從而得到的新的傅里葉描述符,使之適用于開式曲線和閉式曲線的擬合。Yue等[8]構建了由常用平面連桿機構生成的開、閉平面曲線庫,從而將連桿曲線生成問題轉化為庫搜索問題和局部優化問題。
目前采用傅里葉級數法的連桿軌跡曲線綜合中,由于傅里葉級數的性質,多是針對整周采樣點進行的連桿曲線擬合。工程實際中,往往只給出了部分軌跡的采樣點,如在包裝機械中,往往只需要機構滿足一部分軌跡曲線,對整周的軌跡并沒有要求,針對這一問題,研究借鑒了傅里葉描述符的特點,通過與最小二乘法結合,借助遺傳算法,對開式曲線進行擬合。
在傅里葉級數理論中,任何周期函數都可以由正弦函數以及余弦函數組成的無窮級數表示。通過在復平面中建立坐標系,用傅里葉級數的復數形式表示由連桿運動軌跡曲線所形成的復函數,且此過程可逆。
本篇論文的開式曲線的傅里葉描述符采用最小二乘法,通過定義傅里葉描述符對任務軌跡曲線進行擬合,當只給出了部分的連桿運動軌跡時,首先建立該連桿機構任務點的傅里葉級數的復數方程,之后將經擬合后的傅里葉系數與復數方程的系數進行匹配,根據復數方程中桿件的相互關系得出連桿機構的各個參數。
對于平面連桿機構的運動軌跡可以看做一個復函數在復平面上形成的曲線,當連桿曲柄轉動一周時,工作點所形成的曲線可以看做周期函數的采樣,這些采樣點可表示為復數形式z(t)=x(t)+jy(t),一個等時間間隔采樣的點數為N的離散傅里葉逆變換(IDFT)可表示為
(1)
此時z為原函數,An為傅里葉系數,其中n=0,1,…,N-1。則傅里葉系數An可表示為
(2)
以上公式為封閉曲線的傅里葉展開式。在開式曲線中,由于其采樣點未包含整周采樣點信息,我們無法直接使用傅里葉變換。在此,我們使用最小二乘法去近似擬合曲線的傅里葉系數,稱之為開式傅里葉描述符。它們的形式為
(3)
其中:p為諧波成分的階數,且k∈[-p,p];系數αk為傅里葉描述符的系數,代表每個圓周運動的半徑(為|αk|)和初始相位(為∠αk)的復矢量,通過首尾相接,分別以角速度kω0運動從而得到z(t)[9],這些旋轉矢量又叫做諧波分量。
根據離散傅里葉級數理論,低階諧波成分用來描述曲線的輪廓,高階諧波成分用于描述曲線細節。故階數越高,傅里葉表達式的項數越多,對曲線的描述也就越精確。當諧波成分為3階時,對曲線軌跡點的擬合誤差已經小于1%(下文有例證),可以滿足實際需求,若需要進一步提高擬合的精度,則可以通過提高擬合的階數[10]。
最小二乘法又稱最小平方法,通過最小化平方誤差得到與目標數據點匹配的最佳函數匹配,常用于曲線擬合。采用構建誤差函數進行最小二乘法擬合來求解系數αk,通過最小化Δ值,從而得到最佳的系數值αk
(4)
式(4)中的z(t)為原始曲線,可將其轉為復數矩陣形式求解
GX=Y
(5)
對式(5)使用LU分解進行運算[11]
(6)
其中:θi=ω0ti;k,m∈[-p,p];且p為諧波成分的階數。至此,經過解矩陣方程,我們可以得到矩陣X,矩陣X中各個αm即為軌跡曲線各階次的傅里葉系數αk,進而得到擬合曲線的函數表達式z(t)。
▲圖1 一般位置下滑塊五桿機構運動簡圖
如圖1所示,該RRPRR型五桿機構有兩個自由度,以曲柄AB與曲柄DC為主動件,轉動角速度分別為ω1和ω3,與機架AD的初始轉角分別為δ1和δ3,并設δ1和δ3的初始的相位差為δh,其中角β為機架AD與x軸之間的夾角,角θ為連桿BC與機架AD間的夾角,并隨著連桿機構運動的時間而變化。A點坐標為(xa,ya),BC間的距離l2隨滑塊C的運動而變化,滑塊C上固連展桿CP的長度為l0,展桿CP與桿BC之間的夾角為α,當曲柄AB和曲柄CD轉動時,工作點P運動所行成的軌跡即為五桿機構的軌跡曲線。
如圖1,設曲柄AB與曲柄DC轉動的角速度比ω1∶ω3=h∶v,且h和v均為任意整數。δ1和δ3的關系為
δ1=δh+δ3
(7)
(8)
通過ABCDA形成的閉環,得到對應的矢量方程為
(9)
式(9)的復數形式為
(10)
將式(10)取共軛后得到
(11)
將式(9)和式(10)相除消l2后,即可得到以連桿擺動角θ(t)為變量的轉角函數ejθ(t):
(12)
式中:“±”號表示同一機構的兩種構型,取“+”號表示該機構為正裝機構,“-”號表示該機構為反裝機構。本篇論文采用正裝機構,式(12)的轉角函數取正號。
在RRPRR型雙自由度五桿機構中,影響機構工作點P軌跡曲線的因素有:兩原動件角速度比ω1∶ω3以及10個位置參數l1、l3、l4、l0、α、xa、ya、β、δ1、δ3[12]。其中l1、l3、l4為基礎尺寸,l2為關聯尺寸與l1、l3、l4的取值有關,在保證桿長l2滿足連桿持續運動條件下,l2的取值不影響連桿機構運動軌跡。當基礎尺寸給定時,展桿參數l0和α決定工作點P的運動軌跡的形狀和大小,xa、ya和β決定連桿軌跡的位置和方向,δ1的取值與δ3和δh有關,且其中δ1和δ3決定了工作點P的初始位置。
當曲柄AB和曲柄DC轉動方向相同時,滿足該五桿機構能夠連續轉動的桿長條件為[13]
(13)
當連桿機構整周轉動時,轉角函數ejθ(t)的傅里葉變換為
(14)
式(14)中:ck為第k次諧波分量的幅值;φk為第k次諧波分量的初相位;θ(t)是連桿BC相對于機架AD的擺角θ的函數。
當給出的擬合點少于一個周期,即為一個開式曲線時,轉角函數ejθ(t)則不是一個周期函數,不能采用傅里葉變化得到系數ckejφk的值,但可通過第一節中最小二乘法擬合的方法獲得。
如圖2(a)所示,滑塊五桿機構在標準位置下,A點位于坐標原點O點,機架AD與x軸重合,且曲柄DC與機架AD共線,此時曲柄AB與機架AD的初始夾角為δh。曲柄AB與曲柄DC轉動的角速度比ω1∶ω3=h∶v,且h和v均為任意整數。
▲圖2 兩種位置下滑塊五桿機構運動簡圖
代入轉角函數ejθ(t)的傅里葉級數表達式后,此時工作點P點的軌跡曲線則通過傅里葉描述符的形式表達,具體為
rP(t)=rAD+rDC(t)+rCP(t)=l4+l3ejvωt+
l0ejα(c0ejφ0+c1ejφ1ejωt+c-1ejφ-1e-jωt+c2ejφ2ej2ωt+
c-2ejφ-2e-j2ωt+…+ckejφkejkωt)
(15)
當曲柄DC以恒定角速度vω轉動δ3的角位移時,P點的傅里葉級數展開式為
rP(t)′=l4+l3ejδ3ejvωt+l0ej[α+θ(t+δ3/v)]=
l4+l3ejδ3ejvωt+l0ejα[c0ejφ0+c1ej(φ1+δ3/v)ejωt+
c-1ej(φ-1-δ3/v)e-jωt+c2ej(φ2+2δ3/v)ej2ωt+
c-2ej(φ-2-2δ3/v)e-j2ωt+…+ckej(φk+kδ3/v)ejkωt]
(16)
如圖2(b)所示,一般位置下五桿機構的配置中:A點相對于原點O有R的位置偏移以及μ的角偏移,其坐標的向量表示為Rejμ,且機架AD相對于x軸有β的角度;當曲柄DC相對于機架AD以恒定角速度vω轉動δ3的角位移時,曲柄AB以hω的角速度轉動,并相對于機架AD有δ1=δh+hδ3/v的角位移。五桿機構P點軌跡曲線的傅里葉表達式為
rP(t)′=[Rejμ+l4ejβ+l0c0ej(α+β+φ0)]+
l0c1ej(α+β+φ1+δ3/v)ejωt+l0c-1ej(α+β+φ-1-δ3/v)e-jωt+
[l3ej(δ3+β)+l0cvej(α+β+φv+δ3)]ejvωt…+
l0ckej(α+β+φk+kδ3/v)ejkωt+…
(17)
將式(17)中0次和v次諧波分量的系數分別進行合并,經整理得到rp(t)′的傅里葉系數
(18)
其中:Ck=ckejφk,k為諧波成分的次數,且k∈[-p,p],p為正整數。
首先,針對給定的期望軌跡點,傅里葉描述符理論上有無窮多項,但對于實際問題,我們采用有限階數的傅里葉描述符[14],即可較好近似定義任務軌跡曲線R
(19)
其中:p為諧波成分的階數,x(t)+jy(t)即為期望軌跡點R在復平面上的坐標,期望軌跡點傅里葉系數Rk的值可由第一節中的最小二乘法解出。
連桿機構工作點軌跡P通過有限階數的傅里葉級數可以近似表示為
(20)
其中:p為諧波成分的階數,Pk為連桿機構軌跡的傅里葉系數。當任務曲線R與連桿曲線P擬合時,滿足
Rk=Pk,k∈[-p,p]
(21)
于是對于任務曲線R的傅里葉系數,滿足
(22)
其中:Ck=ckejφk,從R的傅里葉系數我們觀察到,R0、Rv、Rk均包含展桿參數l0ej(α+β),定義
(23)
即
K1+jK2=l0ej(α+β)
(24)
針對RRPRR型五桿機構的10個參數,我們設解集
M={l1、l3、l4、xa、ya、δ3、δh、K1、K2、β}
(25)
連桿機構各桿尺寸按比例縮放不影響連桿轉角,為便于程序化求解,將解集M分別定義為解集
(26)
針對多參數的非線性優化問題,目前已經有粒子群算法、蟻群算法、遺傳算法等算法,其中遺傳算法是基于生物進化理論而產生的隨機全局優化算法,該算法適合較復雜的優化問題,相比于傳統優化算法,遺傳算法可以獲得很好的效果[15]。采用遺傳算法計算系數Rk的最小二乘解,通過最小化誤差函數I1,使任務曲線盡可能穿過所提供的目標點
(27)
(28)
即誤差函數的取值與K1、K2無關,由此
(29)
通過優化算法對l3/l1、l4/l1、δ3、δh的值進行搜索后,當I1在許可誤差范圍內時,可得到解集M1{l3/l1、l4/l1、δ3、δh、K1、K2},同時也得到了Ckejkδ3/v的值,由式(22)中的v次諧波成分
Rv=l3ej(δ3+β)+Cvl0ej(α+β)ejδ3/v
(30)
將已知量帶入式(30),分離實部和虛部后相除,消去l3即可獲得β,后代入β即可獲得l3。得到l3后代入解集M1即可得到l1和l4。
由式(23),根據已知K1、K2、β的值,可解出l0和α。
由式(22)中的0次諧波成分
R0=Rejμ+l4ejβ+C0l0ej(α+β)
(31)
通過代入已知參數,將式(31)分離實部虛部可得到Rejμ的值,因
Rejμ=xa+jya
(32)
由式(32)即可得xa和ya。
到此已解出l1、l4、l0、α,解集M2{l3、xa、ya、β}的全部參數,RRPRR型五桿機構的連桿參數已經全部解出。
為進一步驗證開式傅里葉描述符在連桿曲線擬合過程中的實際效果,給出了兩個應用實例,分別為整周采樣點下的擬合與開式曲線下的擬合。
整周采樣點下的擬合相當于擬合角度為360°,對于開式傅里葉描述符同樣適用。給定32個期望軌跡點,如圖3所示,對應位置的坐標見表1。
表1 期望軌跡的取樣點坐標
▲圖3 期望軌跡點
設定實例中曲柄AB和曲柄DC傳動比ω1∶ω3為1∶1。根據給定的坐標點,使用MATLAB編寫程序,利用遺傳算法對坐標點進行擬合,擬合誤差為1.0850e-05,擬合效果如圖4所示。
▲圖4 整周采樣點下的擬合
從32個取樣點中,計算不同階數下的擬合誤差,隨機選取其中7個擬合點為例,列出不同階數下軌跡點的擬合誤差,如表2。
表2 各個位置點在不同階數下的擬合誤差
從表2中可以看出,當階數k為3階時,坐標點的擬合誤差已經小于1%;當階數為5時,擬合誤差已經減少到0.5%以下。
在某紙盒生產過程中,需要將紙張從紙張定位裝置移送到紙堆,從而完成紙張的送料,如圖5所示。
▲圖5 紙張送料示意圖
根據預定的紙張移動軌跡,從軌跡中選取16個點作為期望的軌跡點,并設曲柄AB和曲柄DC的傳動比ω1∶ω3為1∶1,坐標點如表3。
表3 給定的擬合坐標點
因擬合的路徑是開式曲線,我們需要設定其擬合的角度,分別為110°、120°、130°、140°,擬合效果如圖6所示。
根據擬合的結果,當擬合角度為120°時,擬合誤差最小,擬合效果滿足預期,得出的數據如下:
通過分離實部和虛部,可以解出解集M2={l3、xa、ya、β},其中
▲圖6 不同擬合角度的擬合誤差
由K1、K2以及β的值,求解出α以及l0的值。
α=1.109 0,l0=3.810 6
由以上的公式,可以得知該擬合的RRPRR型五桿機構的10個參數,匯總如下:
機架原點A的坐標(xa,ya)為:(41.982 4,11.947 6);
曲柄AB的桿長為17.287 9 mm;
曲柄CD的桿長為13.200 1mm;
機架AD的長度為17.050 2 mm;
機架AD與x軸夾角β為0.460 6 rad;
曲柄AB初始角度δ1為1.137 4 rad;
曲柄CD初始角度δ3為0.209 3 rad;
連桿CP的長度為3.810 6 mm;
連桿CP與BC間夾角α為1.109 0 rad。
因此,本實例選用擬合角度擬合角度為120°時的機構作為紙張輸送機構。
研究采用基于傅里葉描述符的最小二乘法擬合連桿軌跡曲線,并給出了該方法在RRPRR型五桿機構開式曲線的擬合過程,取得了較好的擬合效果,為后續連桿機構開式曲線的擬合提供了可借鑒的方法。在綜合案例中,采用的遺傳算法通過給定不同的擬合角度對其進行擬合求解,從而獲得滿足預期結果的連桿機構。經過改造的傅里葉描述符克服了原有必須取樣整周期采樣點的連桿曲線擬合,更加符合實際工程應用。