結(jié)合Kriging與改進(jìn)NSGA-II的RV減速器優(yōu)化

繆嘉成，李朝陽(yáng)，陳兵奎

(重慶大學(xué) 機(jī)械機(jī)械傳動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，重慶 400044)

RV(rotational vector)減速器是一種效率高、體積小、重量輕、扭轉(zhuǎn)剛度大、傳動(dòng)精度高的新型傳動(dòng)機(jī)構(gòu)，在工業(yè)機(jī)器人、數(shù)控機(jī)床、醫(yī)藥化工設(shè)備等領(lǐng)域應(yīng)用廣泛。其設(shè)計(jì)參數(shù)眾多，約束條件復(fù)雜，傳動(dòng)性能相互耦合，傳統(tǒng)設(shè)計(jì)方法難以獲得最優(yōu)解。隨著優(yōu)化理論逐漸成熟，多目標(biāo)優(yōu)化算法逐漸應(yīng)用于擺線類減速器設(shè)計(jì)。Wang等[1]對(duì)K-H-V擺線減速器進(jìn)行了優(yōu)化，提升了傳動(dòng)效率并縮小了體積。Wang等[2]改進(jìn)了NSGA-II算法，增強(qiáng)了種群的分布性，并用于擺線針輪減速器的優(yōu)化。Jat等[3]使用NSGA-II對(duì)深溝球軸承的基本額定動(dòng)載荷和彈流動(dòng)態(tài)最小膜厚度進(jìn)行了優(yōu)化。

扭轉(zhuǎn)剛度是RV減速器的關(guān)鍵性能指標(biāo)之一[4]。中外學(xué)者針對(duì)擺線及RV減速器剛度特性進(jìn)行了深入研究[4-7]，但未見針對(duì)其剛度優(yōu)化方法的報(bào)道。RV減速器剛度分析通常采用數(shù)值方法或ANSYS有限元仿真[4,6]，前者難以精確反映減速器結(jié)構(gòu)參數(shù)與扭轉(zhuǎn)剛度間的非線性關(guān)系，后者的計(jì)算量難以滿足優(yōu)化算法的要求。為減少耗時(shí)的CAE(computer aided engineering)模擬，需結(jié)合理論模型與CAD (computer aided design)二次開發(fā)，建立部分扭轉(zhuǎn)剛度的Kriging代理模型[8]。

RV減速器結(jié)構(gòu)優(yōu)化的關(guān)鍵是解決多目標(biāo)混合整數(shù)非線性規(guī)劃(MOMINLP, multi-objectives mixed integer non-liner programming)問(wèn)題，求解小樣本問(wèn)題常用分支界定法、割平面法等精確算法。由于精確方法求解高維問(wèn)題的時(shí)間復(fù)雜度極高，中外學(xué)者對(duì)進(jìn)化算法加以改進(jìn)[9]，部分研究基于實(shí)數(shù)編碼的粒子群算法(PSO, particle swarm optimization)或差分進(jìn)化算法(DE, differential evolution algorithm )，利用三角函數(shù)、Sigmod函數(shù)等建立實(shí)數(shù)與整數(shù)的映射關(guān)系[10-11]。

筆者以BAJ-25E為研究對(duì)象，建立RV減速器效率和體積的目標(biāo)函數(shù)以及部分剛度的Kriging代理模型，提出一種能夠同時(shí)處理離散種群、整數(shù)種群與實(shí)數(shù)種群的改進(jìn)NSGA-II算法，并將其應(yīng)用于RV減速器的結(jié)構(gòu)參數(shù)優(yōu)化。

1 結(jié)構(gòu)優(yōu)化的數(shù)學(xué)模型

RV減速器的結(jié)構(gòu)如圖1所示，為保證與現(xiàn)有RV減速器的兼容性，基本參數(shù)由設(shè)計(jì)人員給出，包括輸入功率P(W)、輸入轉(zhuǎn)速n(r/min)、行星輪個(gè)數(shù)np、輸出扭矩T2(mN·m)，表示為

圖1 RV減速器結(jié)構(gòu)原理圖Fig.1 Structure of the RV reducer

(1)

(2)

1.1 RV減速器的性能指標(biāo)

由于RV減速器的體積與扭轉(zhuǎn)剛度的耦合關(guān)系較強(qiáng)，需同時(shí)考慮體積和扭轉(zhuǎn)剛度目標(biāo)，傳動(dòng)效率作為重要性能指標(biāo)也應(yīng)納入優(yōu)化目標(biāo)。

1.1.1 體積

在滿足RV減速器所需的輸入功率和輸出轉(zhuǎn)矩的基礎(chǔ)上，將減小外形尺寸作為設(shè)計(jì)目標(biāo)。如圖2所示，體積目標(biāo)函數(shù)為

圖2 一體化RV減速器外形尺寸Fig.2 Dimensions of the integrated RV reducer

[(d3+2(τ1+τ2+τ3))2-(d3+2τ1)2]h3},

(3)

式中：τ1,τ2,τ3表示針齒殼各部分的厚度；h1,h2,h3表示針齒殼各部分的長(zhǎng)度，h3=δ1+b+δ2+2L′+2B-h1-h2，d3=d1+2d2。δ1和δ2分別表示輸出端盤端面到行星輪端面的距離，以及行星輪端面到支承軸承端面的距離。

1.1.2 扭轉(zhuǎn)剛度

影響RV減速器剛度的元件主要有輸入軸、漸開線行星傳動(dòng)機(jī)構(gòu)、轉(zhuǎn)臂軸承、支承軸承、曲柄軸和擺線針輪傳動(dòng)機(jī)構(gòu)。

1)在輸入扭矩T1作用下，輸入軸的扭轉(zhuǎn)角為

式中di為輸入軸各部分的直徑。

2)通過(guò)數(shù)值彈性力學(xué)求解，漸開線行星傳動(dòng)機(jī)構(gòu)總嚙合剛度為

cr=(0.75εα+0.25)c′，

(4)

式中：εα為端面重合度；c′為單對(duì)齒剛度。

3)轉(zhuǎn)臂軸承為一體化滾子軸承，以圓柱滾子為滾動(dòng)體，曲柄軸為內(nèi)滾道，擺線輪坐標(biāo)孔為外滾道，材料為GCr15。采用單列短圓柱滾子軸承剛度經(jīng)驗(yàn)公式[12]計(jì)算轉(zhuǎn)臂軸承的徑向剛度為

(5)

4)支承軸承同為一體化滾子軸承，考慮變形協(xié)調(diào)條件，輸出端盤剛性較大，視為剛體。np對(duì)支承軸承在輸出端盤扭矩的作用下將產(chǎn)生相等的徑向形變，單個(gè)支撐軸承的徑向力為

(6)

式中asp為漸開線齒輪中心距。將|R′|代入式(5)，可得支承軸承徑向剛度K′r。

5)根據(jù)文獻(xiàn)[14]，曲柄軸的總變形為

(7)

式中：Li為曲柄軸各段的長(zhǎng)度；I為慣性矩；E為彈性模量;Ft為切向分力。

6)由文獻(xiàn)[4]可知，擺線針輪傳動(dòng)處于低速級(jí)，對(duì)減速器扭轉(zhuǎn)剛度的影響較大，故采用CAE仿真計(jì)算擺線輪與針齒間嚙合剛度K″，以提高模型的精度。為提高CAE模型建模效率，開發(fā)了參數(shù)化建模軟件，圖3為典型零件圖形化建模界面。

圖3 RV減速器參數(shù)化建模軟件Fig.3 Parametric modeling software of the RV reducer

用戶界面采用三維建模軟件NX12的Block UI Styler實(shí)現(xiàn)，通過(guò)C++完成用戶界面編輯及尺寸參數(shù)定義，采用Visual Studio對(duì)程序進(jìn)行編譯與鏈接。尺寸驅(qū)動(dòng)的RV減速器模型將創(chuàng)建三維模型的過(guò)程分解，分別對(duì)特征、對(duì)象、實(shí)體的操作進(jìn)行函數(shù)化處理。在建模界面輸入Kriging樣本點(diǎn)對(duì)應(yīng)的尺寸參數(shù)，驅(qū)動(dòng)軟件生成三維零件圖。

將圖3所示擺線針輪傳動(dòng)模型導(dǎo)入仿真軟件中，對(duì)針齒殼施加固定約束，在擺線輪端面施加額定扭矩Tg，其中單個(gè)擺線輪傳遞的扭矩Tg=0.55T2。對(duì)擺線輪和針齒殼進(jìn)行自動(dòng)網(wǎng)格劃分，對(duì)涉及接觸作用的區(qū)域局部細(xì)分，采用Abaqus二次開發(fā)實(shí)現(xiàn)上述流程的自動(dòng)化，如圖4所示。通過(guò)腳本接口建立圖形用戶界面GUI與內(nèi)核的通信，提取擺線輪在扭矩Tg作用下當(dāng)前角度的角位移βH。

圖4 擺線針輪有限元分析程序Fig.4 FEM analysis program of the cycloid-pin

將輸入軸與針齒殼固定，在輸出軸上施加額定扭矩T2，各彈性元件將引起相應(yīng)的彈性轉(zhuǎn)角θi如表1所示，RV減速器的總剛度K′[7]為

表1 各元件引起的彈性轉(zhuǎn)角Table 1 Elastic angle caused by each component

(8)

1.1.3 傳動(dòng)效率

RV減速器的主要效率損失來(lái)自于齒輪嚙合摩擦損失和軸承摩擦損失[15]，表示為

η=η16ηB。

(9)

式中：ηB為軸承總效率，η16為封閉差動(dòng)齒輪傳動(dòng)效率，

1.2 結(jié)構(gòu)參數(shù)的約束條件

1)考慮相鄰行星輪齒頂不干涉(g4)、行星齒輪裝配條件(h1)、彎曲接觸強(qiáng)度(g5,g6)等，如表2所示。

表2 漸開線行星傳動(dòng)機(jī)構(gòu)的約束條件Table 2 Constraints of the involute planetary transmission

2)考慮圓柱滾子安全接觸應(yīng)力(g8)[3]、轉(zhuǎn)臂軸承幾何約束(g10)、軸承壽命(g11)等，如表3所示。

表3 轉(zhuǎn)臂軸承的約束條件Table 3 Constraints of turning arm bearing

表中的Qmax=4.08×10-3|R|/Z；s為安全系數(shù)；[σc]為許用壓應(yīng)力；Do為行星輪傳動(dòng)軸直徑；e為偏心距；nb為轉(zhuǎn)臂軸承內(nèi)外圈相對(duì)轉(zhuǎn)速；壽命指數(shù)ε=10/3；p為平均當(dāng)量動(dòng)載荷；基本額定動(dòng)載荷Cd為

平均當(dāng)量動(dòng)載荷p為

支承軸承的約束條件相似，故不再贅述。

3)考慮針齒分布圓直徑(g12)、擺線輪不根切條件(g16)、擺線針輪接觸強(qiáng)度(g17)等，如表4所示。

表4 擺線針輪傳動(dòng)機(jī)構(gòu)的約束條件Table 4 Constraints of the cycloidal pin transmission

4)考慮兩級(jí)傳動(dòng)的尺寸均衡(g18)，并保證最大輸出轉(zhuǎn)矩(g19)，如表5所示。

表5 整體約束條件Table 5 Overall constraints

2 多目標(biāo)優(yōu)化的理論基礎(chǔ)

2.1 MP-NSGA-II算法

ridPSO和ridDE[10-11]是適用于單目標(biāo)MINLP問(wèn)題的主流算法，不具備多目標(biāo)尋優(yōu)能力。筆者結(jié)合多目標(biāo)進(jìn)化算法NSGA-II提出MP-NSGA-II。該算法繼承了NSGA-II解集分布性良好、時(shí)間復(fù)雜度低、收斂速度快等優(yōu)點(diǎn),并且擴(kuò)展了處理實(shí)數(shù)、整數(shù)及離散變量的能力。

圖5 MP-NSGA-II算法流程圖Fig.5 MP-NSGA-II algorithm flow chart

處理混合種群的具體步驟為：

1)產(chǎn)生初始種群：利用設(shè)計(jì)變量的取值范圍生成自變量范圍矩陣，根據(jù)前m1個(gè)連續(xù)變量和后m2個(gè)離散變量在矩陣中的位置將其分割為2個(gè)子矩陣Mr和Mi。針對(duì)子矩陣Mr，使用rand函數(shù)生成一個(gè)實(shí)數(shù)值的初始種群Pr。針對(duì)子矩陣Mi，使用rand函數(shù)和四舍五入法生成一個(gè)十進(jìn)制整數(shù)的初始種群Pi。

2)混合種群交叉：種群Pr和Pi均是行數(shù)為種群規(guī)模，列數(shù)為設(shè)計(jì)變量數(shù)的矩陣，將Pi轉(zhuǎn)化為浮點(diǎn)數(shù)矩陣，兩矩陣水平合并，生成混合種群Pm,采用兩點(diǎn)交叉實(shí)現(xiàn)個(gè)體間染色體的重組。

3)子種群變異：將混合種群分割為子種群Pr和Pi，Pi轉(zhuǎn)化為整數(shù)矩陣。采用實(shí)數(shù)值高斯變異算子對(duì)Pr進(jìn)行變異，整數(shù)值變異算子實(shí)現(xiàn)矩陣Pi中個(gè)體突變，并將兩矩陣合成混合種群Pm。

由于截?cái)喾╗9]只能處理連續(xù)的整數(shù)變量，為處理MIP(mixed integer programming)問(wèn)題中的離散變量，提出一種基于數(shù)組索引的離散變量通用編碼方案如圖6。將離散變量Zd的nd個(gè)可取值設(shè)為數(shù)組Zarr,數(shù)組索引編碼為取值(0,1,…,nd-1)的整數(shù)子種群Znum。在評(píng)估函數(shù)值和限制條件的階段，根據(jù)數(shù)組Zarr和數(shù)組索引種群Znum解碼出離散子種群Zdis。

圖6 離散變量通用編碼方案Fig.6 Discrete variable universal coding scheme

為增強(qiáng)種群的分布性并降低計(jì)算代價(jià)，引入考慮擁擠距離的非支配排序。由于計(jì)算歐氏距離的效率較低，采用差分計(jì)算目標(biāo)值的偏移量占比代表?yè)頂D距離。根據(jù)混合種群個(gè)體目標(biāo)值由小到大排序，兩相鄰個(gè)體xi和xj間的擁擠距離定義為

(10)

式中：f(x)是個(gè)體x的目標(biāo)值，f(x)max和f(x)min分別是混合種群中最大和最小的目標(biāo)值。根據(jù)擁擠距離更新個(gè)體的適應(yīng)度Vfit，用于選擇下一代個(gè)體。

2.2 雙加點(diǎn)準(zhǔn)則Kriging代理

利用拉丁超立方試驗(yàn)設(shè)計(jì)在設(shè)計(jì)變量空間內(nèi)采樣，得到Kriging代理模型的初始樣本點(diǎn)，樣本點(diǎn)滿足使下式取得最小值：

(11)

Kriging代理模型定義了設(shè)計(jì)變量x與預(yù)測(cè)值y的關(guān)系，表達(dá)式[8]為

y(x)=F(β,x)+z(x)，

(12)

式中：F(β,x)為設(shè)計(jì)變量空間的全局模型，z(x)為按N(0,σ2)隨機(jī)分布的局部偏差,z(x)的統(tǒng)計(jì)特征為

(13)

R(xi,xj)是用于刻畫樣本點(diǎn)xi和xj間關(guān)聯(lián)程度的相關(guān)模型，通常采用高斯相關(guān)模型

(14)

利用線性加權(quán)插值方法得到Kriging模型在預(yù)測(cè)點(diǎn)x處的響應(yīng)值和預(yù)測(cè)方差：

(15)

(16)

式中：R是相關(guān)模型矩陣；r(x)是x點(diǎn)與樣本點(diǎn)間的相關(guān)模型向量；g是樣本點(diǎn)響應(yīng)的向量；q是元素均為1且個(gè)數(shù)為nv的單位列向量。

在優(yōu)化過(guò)程中添加樣本點(diǎn)可有效提高代理模型的精度，通常選擇在期望提高(EI)或均方誤差(MSE)較大處加點(diǎn)。為有效利用進(jìn)化算法優(yōu)化過(guò)程中的信息，將Pareto最優(yōu)集引入加點(diǎn)準(zhǔn)則。在迭代過(guò)程中根據(jù)Pareto集的擁擠距離選擇分散樣本點(diǎn)更新Kriging模型。單次迭代中均方誤差較大處加點(diǎn)數(shù)SE和Pareto最優(yōu)解處加點(diǎn)數(shù)SP分別為

(17)

式中：g是迭代總次數(shù)；gc是當(dāng)前迭代數(shù)；CE和CT分別是調(diào)整系數(shù)和最小加點(diǎn)數(shù)。

2.3 熵權(quán)法Pareto選優(yōu)

熵權(quán)法是一種客觀賦權(quán)方法，利用決策指標(biāo)的熵計(jì)算熵權(quán)值，在i個(gè)性能指標(biāo)X1,X2,…,Xi，j個(gè)Pareto最優(yōu)解的評(píng)價(jià)問(wèn)題中，第p個(gè)性能指標(biāo)Xp={x1,x2,…,xj}，Xp的熵Hp定義為

(18)

(19)

3 優(yōu)化方法的驗(yàn)證及實(shí)現(xiàn)

3.1 MP-NSGA-II評(píng)估

為評(píng)估MP-NSGA-II算法的有效性，對(duì)文獻(xiàn)[11]中的14個(gè)MINLP問(wèn)題進(jìn)行仿真，已知最優(yōu)源自MDE、MDELS、MDEIHS、ridPSO和ridDE。

表6均為最小化目標(biāo)問(wèn)題，MP-NSGA-II的種群規(guī)模為1 000。為體現(xiàn)算法的收斂特性，進(jìn)化代數(shù)設(shè)為10 000。MP-NSGA-II單次運(yùn)行最優(yōu)解與已知最優(yōu)差距小于0.1%，并改進(jìn)了3個(gè)已知最優(yōu)解：P5(x=1.374 823 1,y=1)，P9(x=[27,27,27],y=[88,44])，P12(x=[0.902 19,0.887 75,0.949 18, 0.848 72],y=[5,5,4,6])。由于P7中已知最優(yōu)的自變量x2=0時(shí)目標(biāo)函數(shù)分母為0，原解不成立。

表6 MINLP基準(zhǔn)問(wèn)題仿真Table 6 MINLP benchmark simulation

3.2 Kriging近似模型

圖7 Kriging可視化模型Fig.7 Kriging visual model

采用復(fù)相關(guān)系數(shù)R2檢驗(yàn)擬合模型的精度為

(20)

復(fù)相關(guān)系數(shù)越接近1，近似模型的精度越高。選擇16個(gè)樣本點(diǎn)檢驗(yàn)Kriging模型的精度，如圖8所示。檢驗(yàn)得復(fù)相關(guān)系數(shù)R2為0.920 8，說(shuō)明該代理模型的精度能夠滿足結(jié)構(gòu)優(yōu)化的需要。

圖8 Kriging模型精度檢驗(yàn)Fig.8 Precision validation of the Kriging model

3.3 結(jié)構(gòu)優(yōu)化子程序

RV減速器的優(yōu)化目標(biāo)為{η,V,K′}T，約束函數(shù)為幾何及應(yīng)力約束，優(yōu)化模型為

(21)

式中：g(X,C)為不等式約束；h(X,C)為等式約束。

由上述的優(yōu)化模型及算法原理，以PySide2為開發(fā)框架，編寫一體化結(jié)構(gòu)RV減速器設(shè)計(jì)軟件。其優(yōu)化子程序如圖9所示，由基本參數(shù)及設(shè)計(jì)變量設(shè)置、優(yōu)化算法參數(shù)設(shè)置及優(yōu)化結(jié)果輸出組成。

根據(jù)RV減速器的實(shí)際工況給出RV減速器的動(dòng)力學(xué)及結(jié)構(gòu)基本參數(shù)，如圖9左上。根據(jù)RV減速器的設(shè)計(jì)要求及BAJ-25E減速器的設(shè)計(jì)參數(shù)，初算設(shè)計(jì)變量范圍，如圖9左下。

圖9 一體化結(jié)構(gòu)RV減速器設(shè)計(jì)軟件界面Fig.9 Integrated RV reducer design software interface

設(shè)置MP-NSGA-II算法的種群規(guī)模為3 000，遺傳代數(shù)為1 000，代溝為0.5，交叉概率為90%，變異概率為10%。選擇方式為輪盤賭選擇，重組方式為兩點(diǎn)交叉。Pareto前沿如圖10所示，可知優(yōu)化目標(biāo)間相互制約，需從解集中優(yōu)選理想解。

圖10 Pareto前沿Fig.10 Pareto frontier

生成RV減速器的結(jié)構(gòu)及性能參數(shù)表如圖9右下所示。設(shè)定傳動(dòng)效率及扭轉(zhuǎn)剛度下限和體積上限，以縮減Pareto最優(yōu)解集的規(guī)模。

4 結(jié)果與分析

各優(yōu)化目標(biāo)通過(guò)結(jié)構(gòu)參數(shù)相互耦合，利用兩目標(biāo)的Pareto前沿分析優(yōu)化目標(biāo)間的耦合關(guān)系，為后續(xù)的Pareto選優(yōu)提供依據(jù)，如圖11所示。

圖11 各優(yōu)化目標(biāo)間的耦合關(guān)系Fig.11 Coupling relationships between optimization objectives

從圖11(a)和(b)可看出扭轉(zhuǎn)剛度和體積變化時(shí)，傳動(dòng)效率維持在85.2%～85.6%，與其余優(yōu)化目標(biāo)的耦合關(guān)系不顯著。11(c)表明體積與扭轉(zhuǎn)剛度成顯著正相關(guān)。因此，主要考慮體積和扭轉(zhuǎn)剛度的設(shè)計(jì)要求，初步篩選5個(gè)設(shè)計(jì)方案如表7所示。

表7 初選設(shè)計(jì)方案Table 7 Preliminary design selection

采用熵權(quán)法，將性能指標(biāo)決策矩陣歸一化，求解各指標(biāo)的信息熵Hp，并計(jì)算各目標(biāo)熵權(quán)值：

wp=[0.293 2,0.318 6,0.388 2],p=1,2,3。

由熵權(quán)值wp及歸一化指標(biāo)rpq對(duì)各方案評(píng)分，為

Zq=r1qw1-r2qw2+r3qw3(q=1,2,3,4,5)。

(22)

根據(jù)Zq排序，選擇理想方案2。BAJ-25E的結(jié)構(gòu)參數(shù)與優(yōu)化后的RV減速器參數(shù)如表8所示。

表8 優(yōu)化參數(shù)對(duì)比及敏感性分析Table 8 Comparison of optimized parameters and sensitivity analysis

由表8知，經(jīng)MP-NSGA-II算法優(yōu)化可直接獲得符合約束條件的優(yōu)化值，無(wú)需圓整處理。與初始值相比，優(yōu)化解的傳動(dòng)效率提升了1.24%，體積減小了1.69%，扭轉(zhuǎn)剛度增大了53.83%。

在制造過(guò)程中，尺寸參數(shù)可能出現(xiàn)1%左右的偏差[3]，對(duì)性能指標(biāo)的敏感性分析有助于指導(dǎo)實(shí)際生產(chǎn)中的誤差控制。研究連續(xù)變量變化±1%對(duì)優(yōu)化后的減速器扭轉(zhuǎn)剛度的影響，如表8?？梢奃z、K1、L和L′對(duì)扭轉(zhuǎn)剛度的影響均超過(guò)0.5%，其它參數(shù)波動(dòng)的影響較小，均低于0.2%。

5 結(jié) 論

1)分析了影響RV減速器性能的結(jié)構(gòu)參數(shù)，綜合17個(gè)設(shè)計(jì)變量，3個(gè)目標(biāo)函數(shù)和21個(gè)約束條件建立了結(jié)構(gòu)優(yōu)化的數(shù)學(xué)模型。

2)提出MP-NSGA-II算法，利用MINLP基準(zhǔn)問(wèn)題測(cè)試了改進(jìn)算法的性能，并改進(jìn)了3個(gè)已知最優(yōu)解，證明了該算法的有效性。

3)結(jié)合Kriging與MP-NSGA-II得到Pareto前沿，利用熵權(quán)法完成Pareto選優(yōu)，有效提高了RV減速器的綜合性能。