科學(xué)計(jì)算選講優(yōu)化算法在并聯(lián)機(jī)器人尺度綜合中的應(yīng)用董成林天津大學(xué)

1、優(yōu)化算法在并聯(lián)機(jī)器人尺度綜合中的應(yīng)用董成林天津大學(xué)機(jī)械工程學(xué)院10132021601工程問(wèn)題如圖1所示，是一種具有1T2R 3自由度的3-UPU位置型并聯(lián)機(jī)構(gòu)，該并聯(lián)機(jī)構(gòu)動(dòng)平臺(tái)與3自由度轉(zhuǎn)頭串接可組成6自由度混聯(lián)機(jī)器人，具有工作空間/裝備占地比大、精度和剛度高等優(yōu)點(diǎn)圖1 3-TPT機(jī)器人圖2 3-TPT機(jī)器人結(jié)構(gòu)簡(jiǎn)圖在并聯(lián)機(jī)構(gòu)的設(shè)計(jì)中，尺度綜合是核心內(nèi)容之一，其涉及的主要問(wèn)題是：在 滿足一定機(jī)構(gòu)約束的條件下，利用各種性能評(píng)價(jià)指標(biāo)構(gòu)造目標(biāo)函數(shù)， 通過(guò)優(yōu)化尺 度參數(shù)使機(jī)構(gòu)在其工作空間內(nèi)獲得優(yōu)良的操作性能。通常可采用的性能指標(biāo)有： 雅可比矩陣的條件數(shù)(condition number)、最小奇異值(

2、minimum singular value)、 可操作性(manipulability)、各向同性(isotropy)等。在對(duì)該并聯(lián)機(jī)構(gòu)進(jìn)行尺度綜合 之前，要先求得機(jī)構(gòu)的無(wú)量綱雅可比矩陣，因此要對(duì)機(jī)構(gòu)進(jìn)行位置分析和速度分析。具體計(jì)算過(guò)程這里不作詳述，只給出最后求得的雅可比矩陣:Ja J c Jcv(1)式中,J cvJ av =然后確定該機(jī)器人的工作空間及設(shè)計(jì)變圖3 3-TPT并聯(lián)機(jī)構(gòu)的任務(wù)空間qimin-(Cx S3,1 T-(F S3,2 )T (C3 匯 S3,3 )IT(q +q%i 嚴(yán) S5)T一(C2 * q2S3,2 產(chǎn) S5,2 )TI (C3 +q3S3,3 )漢 S5,3

3、 )量與運(yùn)動(dòng)學(xué)性能評(píng)價(jià)指標(biāo)。參見(jiàn)圖3,根據(jù)3-TPT機(jī)構(gòu)結(jié)構(gòu)特點(diǎn)，定義其任務(wù)空間 W為 一半徑為R ,高為h的圓柱體，且令點(diǎn)B至W 下平面距離為H。顯見(jiàn)，在初始位形(即系 A-xiyiZi與B-xyz平行)下，點(diǎn)P處于W的 中軸線上。利用位置逆解模型在 Wt邊界上搜索，可確定支鏈桿長(zhǎng)極限qimin和qimax(i =123)做如下歸化處理，即a2bHR印de'a2 =，飛，H，R，a，d，e -b2b?b?b?a?a 2a 2約束條件：'a,= 1 _ ' d扎.=扎 ba,，e = 27 =qimax_qimin 乞%，i =1,2,3為了便于機(jī)構(gòu)尺度的參數(shù)化設(shè)計(jì)，用

4、b2對(duì)、b和R，及用a2對(duì)a、d和e性能評(píng)價(jià)指標(biāo):式中，V表示W(wǎng)的體積，m為無(wú)量綱雅克比矩陣J Pv的最小奇異值。為了控制操作速度上界，應(yīng)使cm ' max(7)至此，3-TPT并聯(lián)機(jī)構(gòu)的尺度綜合可歸結(jié)為：給定 帯二hR、上二及相 應(yīng)的約束條件，確定無(wú)量綱參數(shù)機(jī)構(gòu)的獨(dú)立設(shè)計(jì)變量為、'a和冷2，以使其操作性能在給定的工作空間中最優(yōu)。即以 時(shí)為目標(biāo)函數(shù)，以支鏈桿長(zhǎng)比及動(dòng)靜平臺(tái)形狀比等為約束的非線性多約束單目標(biāo)優(yōu)化問(wèn)題：:r min(8)CTm ( X )+ >GflXi-Yrtr,上坷=I _-d ,上q =仏s.t.片蘭 , i =1,2,30.5 W <1 , 0.

5、25蘭扎a?蘭0.5, 1.5蘭盯 2.5T *_X ='色'a2 b1' d ' H ，二二 max；二 m X2非線性多約束單(多)目標(biāo)優(yōu)化在工程中，存在這么一類(lèi)優(yōu)化問(wèn)題：f(x) > min(9)gj(x0j Is.t.Wj(x) =0 j = L其中，x =(為，x2xn)T Rn， f : Rn R， gj : Rn ； R(j I U L)， I =1,2,.m， L =m 1,m 2.m p。該問(wèn)題即為非線性約束優(yōu)化問(wèn)題。到目前為止，還沒(méi)有特別有效的方法直接 得到最優(yōu)解，人們普遍采用迭代的方法求解：首先選擇一個(gè)初始點(diǎn)，利用當(dāng)前迭 代點(diǎn)的或已產(chǎn)

6、生的迭代點(diǎn)的信息，產(chǎn)生下一個(gè)迭代點(diǎn)，一步一步逼近最優(yōu)解，進(jìn)而得到一個(gè)迭代點(diǎn)列，這樣便構(gòu)成求解(9)的迭代算法。針對(duì)本文第一節(jié)中所提出的問(wèn)題，可以借助于matlab優(yōu)化工具包中的fmincon函數(shù)，該函數(shù)主要用于解決非線性多約束單目標(biāo)優(yōu)化問(wèn)題。對(duì)于大規(guī)模 運(yùn)算，該函數(shù)采用信賴域法進(jìn)行求解；中小規(guī)模預(yù)算，那么采用序列二次規(guī)劃算法。 下面對(duì)這兩種算法進(jìn)行介紹：2.1信賴域算法信賴域方法是一類(lèi)很新的方法，它和線搜索法并列為目前求解非線性規(guī)劃問(wèn)題的兩類(lèi)主要數(shù)值方法。信賴域方法思想新穎，算法可靠，具有很強(qiáng)的收斂性.它不僅能很快地解決良態(tài)問(wèn)題，而且也能有效地求解病態(tài)(ill-conditioned)的優(yōu)化問(wèn)

7、題。下面給出信賴域算法的一般步驟：Step1.初始化選取初始點(diǎn)x0，對(duì)稱陣B。 Rn n，最大信賴域半徑厶 0 ;取初始信賴域半 徑禮* (0)，容許誤差；- 0，參數(shù)0,1/ 4);置k = 0。Step2.收斂性檢驗(yàn)假設(shè)|八廣|2乞；，那么Xk為近似最優(yōu)解，算法停止；否那么，求解k T 1 Tmin mk(d) = fk C f ) d qd Bkds.t |d|25得到可能的移動(dòng)量dk，在求解(10)時(shí)，采用一種基于共軛梯度法的算法，即Steilhaug:1) .初始化：取初始最優(yōu)解d。=0, g0，初始搜索方向?yàn)閞0=-g°，容 許誤差M 0，置k = 0 o2) .收斂性檢

8、驗(yàn)：假設(shè)|gk |，那么取近似解d = dk，算法停止。假設(shè)|rk|B = (rk)TBr0那么沿著直線d =dk .rk求解(10)，得到一個(gè)滿足|d |.:的近似解，算法停止。 其他情形下，令% 咄gk|2/|rk|B，dk+ = dkrk假設(shè)|dk1|2-：那么尋找步長(zhǎng) 一 0 ,使得|dk"rk|2“并且令近似解d =dk rk算法停止；否那么，轉(zhuǎn)3) o3) .迭代改良：令 gk41 = gk+c(kBrk，假設(shè)|g心 |2« |g0|那么近似解d二dk1，算法停止；否那么，令k =|gk1|2/|gk|2，rk1gk1置 k 1 > k，轉(zhuǎn) 2)oStep

9、3.迭代改良根據(jù)實(shí)際下降量：AReck 二 f(xk) - f (xk dk)預(yù)期下降量PReck 二 mJ0)-mk(dk)得到下降量比匚二 ARe4 /PReck令新的迭代點(diǎn)和信賴域半徑分別為口dk,pk>nX沖=2-xk,Pk 0和r"kl|dk|/4,Pk 小4Ak+=<min2Ak,&, Pk >3/4。Ak,elseStep4.校正矩陣令fk" Jf(xk 1),通過(guò)修正矩陣Bk ,得到新的對(duì)稱矩陣Bki ；置k k , 轉(zhuǎn) Step2.2.2序列二次規(guī)劃法序列二次規(guī)劃法，簡(jiǎn)稱SQP方法，亦稱約束變尺度法。下面就問(wèn)題(9)的SQP 算法

10、作簡(jiǎn)要介紹。SQP算法的一個(gè)重要優(yōu)點(diǎn)就是：如果沿著二次規(guī)劃子問(wèn)題的解方向搜索最終 能得到步長(zhǎng)為1,那么算法就是超線性收斂的。下面給出一般的SQP算法的根本步驟。Stepl.給定初始點(diǎn)X。 Rn，選取正定矩陣Bo，令k = 0.Step2.求解子規(guī)劃Q(Xk,Bk)，即1min i f (xk)T ddTBkd2gj(xQ '、gj(Xk)T0, j Is.t.Tgj(xQ '、gj(xQ d =0, j l得到dk.Step3.令Xk 廠Xk ： kdk( ： k -0).Step4.修正Bk，使Bk 1正定，令 k 1，返回Step3.3工程問(wèn)題求解設(shè)定R =1.2、0.3、

11、 =1，調(diào)用MATLAB優(yōu)化工具包中的基于以上算 法的有約束非線性規(guī)劃fmineon函數(shù)，可得到如下一組使得> min的最優(yōu)解答0.6331 '，二=0.35810.25000.63310.3669-1.7381圖4示出了在取最優(yōu)尺度參數(shù)時(shí) 二m在Wt中上、中、下三個(gè)截面的分布。由圖可見(jiàn)，二m關(guān)于X軸對(duì)稱，在Wt中軸線上取得最大值，且沿徑向和Z坐標(biāo)增大單調(diào)減小；機(jī)構(gòu)在整個(gè)任務(wù)空間內(nèi)具有良好的全域運(yùn)動(dòng)學(xué)性能1.510.5 y-0 5-1-1.50.40330320.310.5y-0.5z三 H 0.5h010.380.360.340.320.30 281.510.340.5y0.3

12、-0.50.28-10.32z = H h10 12圖4二m在Wt的上、中、下三個(gè)截面的分布附錄：matlab程序代碼% 主函數(shù)“ f_TPT_slove ：x0=0.75,0.3,0.75,0.25,2;1=0.5,0.25,0.5,0,1.5;u=1,0.5,1,0.5,2.5;optio ns=optimset('Maxlter',15);xo2,yo2=fmi neon ('f_TPT_mi nsi ngn um',xO,l,u,'f_TPT_ nlco n',optio ns) %子函數(shù)1 (目標(biāo)函數(shù))“ f_TPT_minsingnu

13、m RR=1.2*b2;hh=0.3*RR;s051=1;0;0; va01=a1*0;-1;0; vb1= b1*0;-1;0;s052=0;1;0; va02=a2*1;0;0; vb2=b2*1;0;0;s053=-1;0;0;va03=a1*0;1;0;vb3=b1*0;1;0;n r=15; n h=4; n p=20;dr=RR/(nr-1); dh=hh/(nh-1);dp=2*pi/(np-1);for i=1: nhh=(i-1)*dh+H;%工作空間軸向掃描for j=1: npalfa=(j-1)*dp;%任務(wù)空間周向掃描；for k=1: nrr=(k-1)*dr;%任

14、務(wù)空間徑向掃描；% Kin ematic an alysis xp=r*cos(alfa)+d; yp=r*si n( alfa); zp=h; rp=xp;yp;zp;thita 1=0.1;phi1=0.1;psi1=0;if xp+d=0thita 1=0;elsethita1=2*ata n( (-1*zp+sqrt(zp*zp+xp*xp-d*d)/(xp+d); endphi1=-1*ata n(yp/(xp-b2)*si n(thita1)+zp*cos(thita1);R(:,:,3)=Rotatio n(1,psi1)*Rotatio n(2,thita1)*Rotatio

15、n(1,phi1);vx1(:,3)=R(:,1,3);vy1(:,3)=R(:,2,3);vz1(:,3)=R(:,3,3);%支鏈桿長(zhǎng)求解va仁 R(:,:,3)*va01;va2=R(:,:,3)*va02;va3=R(:,:,3)*va03;vc1=-va1+d*vx1(:,3)+e*vz1(:,3);vc2=-va2+d*vx1(:,3)+e*vz1(:,3);vc3=-va3+d*vx1(:,3)+e*vz1(:,3);q1= norm(rp-vbl-vcl);q2=no rm(rp_vb2_vc2);q3=no rm(rp-vb3-vc3);si 1=1;0;0; s31=(rp

16、-vb1-vc1)/q1; s41=s31;s12=0;1;0; s32=(rp-vb2-vc2)/q2;s42=s32;s13=-1;0;0;s33=(rp-vb3-vc3)/q3;s43=s33;s21=cross(s31,s11)/norm(cross(s31,s11);s22=cross(s32,s12)/norm(cross(s32,s12);s23=cross(s33,s13)/norm(cross(s33,s13);s53=R(:,:,3)*sO53;s5 仁R(:,:,3)*s051;s52=R(:,:,3)*sO52;Jav=s31's32's33'J

17、aw=-(cross(vc1,s31)'-(cross(vc2,s32)'-(cross(vc3,s33)' Jcv=s51's52's53'Jcw=-(cross(vc1+q1*s31,s51)'-(cross(vc2+q2*s32,s52)'-(cross(vc3+q3*s33,s53)'Jpv=Jav-Jaw* in v(Jcw)*Jcv; ss=svd(Jpv);minsin gularit ynu mber(i,j,k)=min( abs(ss);end end end mea n=sum(sum(sum(mi

18、nsi ngularitynu mber)/( nr)* (n p)*( nh); stdev=sqrt(sum(sum(sum(mi nsi ngularitynu mber-mea n).A2)/( nr)*( np)*( nh);ff(1)=-mea n;ff(2)=stdev;%子函數(shù)約束函數(shù)“ f_TPT_nlcon function c,ceq=f_TPT_ nlcon(x) la1=x(1); la2=x(2); lb 1=x(3); ld=x(4); lH=x(5);b2=1; a2=b2*la2; b1=b2*lb1; a1=a2*la1; d=a2*ld;H=b2*IH;

19、e=2*a2;RR=1.2*b2; hh=0.3*RR;s051=1;0;0; va01=a1*0;-1;0; s052=0;1;0; va02=a2*1;0;0; s053=-1;0;0; va03=a1*0;1;0;vb仁 b1*0;-1;0;vb2=b2*1;0;0;vb3=b1*0;1;0;nr=15; nh=2; np=30;dr=RR/(nr-1); dh=hh/(nh-1);dp=2*pi/(np-1);for i=1: nhh=(i-1)*dh+H;%工作空間軸向掃描；for j=1: npalfa=(j-1)*dp;%任務(wù)空間周向掃描；for k=1: nrr=(k-1)*dr;%任務(wù)空間徑向掃描；% Kin ematic an alysis xp=r*cos(alfa)+d; yp=r*si n( alfa); zp=h;rp=xp;yp;zp;thita 1=0.1;phi1=0.1;psi1=0;if xp+d=0thita 1=0;elsethita1=2*ata n( (-1*zp+sqrt(zp*zp+xp*xp-d*d)/(xp+d); endphi1=-1*ata n(yp/(xp-b2)*si n(thita1)+zp*cos(th