結合Kriging與改進NSGA-II的RV減速器優化

2021-03-13 08:37:52繆嘉成李朝陽陳兵奎

重慶大學學報 2021年2期

繆嘉成，李朝陽，陳兵奎

(重慶大學機械機械傳動國家重點實驗室，重慶 400044)

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

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

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

筆者以BAJ-25E為研究對象，建立RV減速器效率和體積的目標函數以及部分剛度的Kriging代理模型，提出一種能夠同時處理離散種群、整數種群與實數種群的改進NSGA-II算法，并將其應用于RV減速器的結構參數優化。

1 結構優化的數學模型

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

圖1 RV減速器結構原理圖Fig.1 Structure of the RV reducer

(1)

(2)

1.1 RV減速器的性能指標

由于RV減速器的體積與扭轉剛度的耦合關系較強，需同時考慮體積和扭轉剛度目標，傳動效率作為重要性能指標也應納入優化目標。

1.1.1 體積

在滿足RV減速器所需的輸入功率和輸出轉矩的基礎上，將減小外形尺寸作為設計目標。如圖2所示，體積目標函數為

圖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表示針齒殼各部分的長度，h3=δ1+b+δ2+2L′+2B-h1-h2，d3=d1+2d2。δ1和δ2分別表示輸出端盤端面到行星輪端面的距離，以及行星輪端面到支承軸承端面的距離。

1.1.2 扭轉剛度

影響RV減速器剛度的元件主要有輸入軸、漸開線行星傳動機構、轉臂軸承、支承軸承、曲柄軸和擺線針輪傳動機構。

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

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

2)通過數值彈性力學求解，漸開線行星傳動機構總嚙合剛度為

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

(4)

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

3)轉臂軸承為一體化滾子軸承，以圓柱滾子為滾動體，曲柄軸為內滾道，擺線輪坐標孔為外滾道，材料為GCr15。采用單列短圓柱滾子軸承剛度經驗公式[12]計算轉臂軸承的徑向剛度為

(5)

4)支承軸承同為一體化滾子軸承，考慮變形協調條件，輸出端盤剛性較大，視為剛體。np對支承軸承在輸出端盤扭矩的作用下將產生相等的徑向形變，單個支撐軸承的徑向力為

(6)

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

5)根據文獻[14]，曲柄軸的總變形為

(7)

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

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

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

用戶界面采用三維建模軟件NX12的Block UI Styler實現，通過C++完成用戶界面編輯及尺寸參數定義，采用Visual Studio對程序進行編譯與鏈接。尺寸驅動的RV減速器模型將創建三維模型的過程分解，分別對特征、對象、實體的操作進行函數化處理。在建模界面輸入Kriging樣本點對應的尺寸參數，驅動軟件生成三維零件圖。

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

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

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

表1 各元件引起的彈性轉角Table 1 Elastic angle caused by each component

(8)

1.1.3 傳動效率

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

η=η16ηB。

(9)

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

1.2 結構參數的約束條件

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

表2 漸開線行星傳動機構的約束條件Table 2 Constraints of the involute planetary transmission

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

表3 轉臂軸承的約束條件Table 3 Constraints of turning arm bearing

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

平均當量動載荷p為

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

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

表4 擺線針輪傳動機構的約束條件Table 4 Constraints of the cycloidal pin transmission

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

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

2 多目標優化的理論基礎

2.1 MP-NSGA-II算法

ridPSO和ridDE[10-11]是適用于單目標MINLP問題的主流算法，不具備多目標尋優能力。筆者結合多目標進化算法NSGA-II提出MP-NSGA-II。該算法繼承了NSGA-II解集分布性良好、時間復雜度低、收斂速度快等優點,并且擴展了處理實數、整數及離散變量的能力。

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

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

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

2)混合種群交叉：種群Pr和Pi均是行數為種群規模，列數為設計變量數的矩陣，將Pi轉化為浮點數矩陣，兩矩陣水平合并，生成混合種群Pm,采用兩點交叉實現個體間染色體的重組。

3)子種群變異：將混合種群分割為子種群Pr和Pi，Pi轉化為整數矩陣。采用實數值高斯變異算子對Pr進行變異，整數值變異算子實現矩陣Pi中個體突變，并將兩矩陣合成混合種群Pm。

由于截斷法[9]只能處理連續的整數變量，為處理MIP(mixed integer programming)問題中的離散變量，提出一種基于數組索引的離散變量通用編碼方案如圖6。將離散變量Zd的nd個可取值設為數組Zarr,數組索引編碼為取值(0,1,…,nd-1)的整數子種群Znum。在評估函數值和限制條件的階段，根據數組Zarr和數組索引種群Znum解碼出離散子種群Zdis。

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

為增強種群的分布性并降低計算代價，引入考慮擁擠距離的非支配排序。由于計算歐氏距離的效率較低，采用差分計算目標值的偏移量占比代表擁擠距離。根據混合種群個體目標值由小到大排序，兩相鄰個體xi和xj間的擁擠距離定義為

(10)

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

2.2 雙加點準則Kriging代理

利用拉丁超立方試驗設計在設計變量空間內采樣，得到Kriging代理模型的初始樣本點，樣本點滿足使下式取得最小值：

(11)

Kriging代理模型定義了設計變量x與預測值y的關系，表達式[8]為

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

(12)

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

(13)

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

(14)

利用線性加權插值方法得到Kriging模型在預測點x處的響應值和預測方差：

(15)

(16)

式中：R是相關模型矩陣；r(x)是x點與樣本點間的相關模型向量；g是樣本點響應的向量；q是元素均為1且個數為nv的單位列向量。

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

(17)

式中：g是迭代總次數；gc是當前迭代數；CE和CT分別是調整系數和最小加點數。

2.3 熵權法Pareto選優

熵權法是一種客觀賦權方法，利用決策指標的熵計算熵權值，在i個性能指標X1,X2,…,Xi，j個Pareto最優解的評價問題中，第p個性能指標Xp={x1,x2,…,xj}，Xp的熵Hp定義為

(18)

(19)

3 優化方法的驗證及實現

3.1 MP-NSGA-II評估

為評估MP-NSGA-II算法的有效性，對文獻[11]中的14個MINLP問題進行仿真，已知最優源自MDE、MDELS、MDEIHS、ridPSO和ridDE。

表6均為最小化目標問題，MP-NSGA-II的種群規模為1 000。為體現算法的收斂特性，進化代數設為10 000。MP-NSGA-II單次運行最優解與已知最優差距小于0.1%，并改進了3個已知最優解：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中已知最優的自變量x2=0時目標函數分母為0，原解不成立。