RV(rotational vector)减速器是一种效率高、体积小、重量轻、扭转刚度大、传动精度高的新型传动机构，在工业机器人、数控机床、医药化工设备等领域应用广泛。其设计参数众多，约束条件复杂，传动性能相互耦合，传统设计方法难以获得最优解。随着优化理论逐渐成熟，多目标优化算法逐渐应用于摆线类减速器设计。Wang等^[1]对K-H-V摆线减速器进行了优化，提升了传动效率并缩小了体积。Wang等^[2]改进了NSGA-Ⅱ算法，增强了种群的分布性，并用于摆线针轮减速器的优化。Jat等^[3]使用NSGA-Ⅱ对深沟球轴承的基本额定动载荷和弹流动态最小膜厚度进行了优化。

扭转刚度是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-Ⅱ算法，并将其应用于RV减速器的结构参数优化。

1 结构优化的数学模型

RV减速器的结构如图 1所示，为保证与现有RV减速器的兼容性，基本参数由设计人员给出，包括输入功率P(W)、输入转速n(r/min)、行星轮个数n_p、输出扭矩T₂(mN·m)，表示为

$ \mathit{\boldsymbol{C}} = {\{ P,n,{n_{\rm{p}}},{T_2}\} ^{\rm{T}}}。$

(1)

RV减速器的中心轮齿数z₁、行星轮齿数z₂、摆线轮齿数z_g、针齿中心圆直径D_z、针齿直径d_z和短幅系数K₁直接影响机械效率，z₁、z₂、D_z、d_z及行星轮齿宽b、渐开线齿轮模数m、摆线轮齿宽B及转臂轴承节圆直径D_m、转臂轴承圆柱滚子直径D_r、转臂轴承圆柱滚子数量Z、转臂轴承圆柱滚子有效长度L、支承轴承节圆直径D′_m、支承轴承圆柱滚子直径D′_r、支承轴承圆柱滚子数量Z′、支承轴承圆柱滚子有效长度L′直接影响减速器的体积和扭转刚度。设计变量表示为

$ \mathit{\boldsymbol{X}} = {\left\{ {{z_1},{z_2},b,m,{z_{\rm{g}}},{D_{\rm{z}}},{d_{\rm{z}}},B,{K_1},{D_{\rm{m}}},{D_{\rm{r}}},Z,L,D_{\rm{m}}^\prime ,D_{\rm{r}}^\prime ,{Z^\prime },{L^\prime }} \right\}^{\rm{T}}}。$

(2)

1.1 RV减速器的性能指标

由于RV减速器的体积与扭转刚度的耦合关系较强，需同时考虑体积和扭转刚度目标，传动效率作为重要性能指标也应纳入优化目标。

1.1.1 体积

在满足RV减速器所需的输入功率和输出转矩的基础上，将减小外形尺寸作为设计目标。如图 2所示，体积目标函数为

$ \begin{array}{l} V = \frac{{\rm{ \mathsf{ π} }}}{4}\left\{ {{{\left( {{d_3} + 2{\tau _1}} \right)}^2}\left( {{h_1} + {h_2} + {h_3}} \right) + \left[ {{{\left( {{d_1} + 2\left( {{\tau _1} + {\tau _2}} \right)} \right)}^2} - {{\left( {{d_3} + 2{\tau _1}} \right)}^2}} \right]{h_2} + } \right.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {\left[ {{{\left( {{d_3} + 2\left( {{\tau _1} + {\tau _2} + {\tau _3}} \right)} \right)}^2} - {{\left( {{d_3} + 2{\tau _1}} \right)}^2}} \right]{h_3}} \right\}, \end{array} $

(3)

式中：τ₁, τ₂, τ₃表示针齿壳各部分的厚度；h₁, h₂, h₃表示针齿壳各部分的长度，h₃=δ₁+b+δ₂+2L′+2B-h₁-h₂，d₃=d₁+2d₂。δ₁和δ₂分别表示输出端盘端面到行星轮端面的距离，以及行星轮端面到支承轴承端面的距离。

1.1.2 扭转刚度

影响RV减速器刚度的元件主要有输入轴、渐开线行星传动机构、转臂轴承、支承轴承、曲柄轴和摆线针轮传动机构。

1) 在输入扭矩T₁作用下，输入轴的扭转角为

$ {\theta _{{\rm{si}}}} = {T_1}\sum\nolimits_{i = 1}^3 {\frac{{{l_i}}}{{G{I_{\rm{P}}}}}} , $

式中：l_i是输入轴各部分的长度；G为切变模量；I_P=πd_v⁴/32，其中d_v为等效直径：

$ {d_{\rm{v}}} = \sqrt[4]{{\left. {\sum\nolimits_{i = 1}^3 {{l_i}} /\sum\nolimits_{i = 1}^3 {{l_i}} /d_i^4} \right)}}。$

式中d_i为输入轴各部分的直径。

2) 通过数值弹性力学求解，渐开线行星传动机构总啮合刚度为

$ {c_{\rm{r}}} = \left( {0.75{\varepsilon _\alpha } + 0.25} \right){c^\prime }, $

(4)

式中：ε_α为端面重合度；c′为单对齿刚度。

3) 转臂轴承为一体化滚子轴承，以圆柱滚子为滚动体，曲柄轴为内滚道，摆线轮坐标孔为外滚道，材料为GCr15。采用单列短圆柱滚子轴承刚度经验公式^[12]计算转臂轴承的径向刚度为

$ {K_{\rm{r}}} = 0.340 \times {10^4}(|R{|^{0.1}}{Z^{0.9}}{L^{0.8}}{\rm{co}}{{\rm{s}}^{1.9}}\beta ), $

(5)

式中：转臂轴承径向合力^[13]$ \left| R \right| = \left| {{R_{\rm{t}}}} \right| + \sqrt {{{\left| {{R_x}} \right|}^2} + {{\left| {{R_y}} \right|}^2}} $，Meng^[6]给出了计算转臂轴承切向受力R_t、x向受力R_x和y向受力R_y的方法；一体化轴承滚动体接触角β=0°。

4) 支承轴承同为一体化滚子轴承，考虑变形协调条件，输出端盘刚性较大，视为刚体。n_p对支承轴承在输出端盘扭矩的作用下将产生相等的径向形变，单个支撑轴承的径向力为

$ |{R^\prime }| = \frac{{{T_2}}}{{2{n_{\rm{p}}}{a_{{\rm{sp}}}}}}, $

(6)

式中a_sp为渐开线齿轮中心距。将|R′|代入式(5)，可得支承轴承径向刚度K′_r。

5) 根据文献[14]，曲柄轴的总变形为

$ \delta = \frac{{|R| \cdot 2L_2^2L_3^2}}{{6EI({L_2} + {L_3})}} - \frac{{{F_{\rm{t}}} \cdot {L^2}1{\kern 1pt} {\kern 1pt} {L_2}}}{{2EI}}, $

(7)

式中：L_i为曲柄轴各段的长度；I为惯性矩；E为弹性模量; F_t为切向分力。

6) 由文献[4]可知，摆线针轮传动处于低速级，对减速器扭转刚度的影响较大，故采用CAE仿真计算摆线轮与针齿间啮合刚度K″，以提高模型的精度。为提高CAE模型建模效率，开发了参数化建模软件，图 3为典型零件图形化建模界面。

用户界面采用三维建模软件NX12的Block UI Styler实现，通过C++完成用户界面编辑及尺寸参数定义，采用Visual Studio对程序进行编译与链接。尺寸驱动的RV减速器模型将创建三维模型的过程分解，分别对特征、对象、实体的操作进行函数化处理。在建模界面输入Kriging样本点对应的尺寸参数，驱动软件生成三维零件图。

将图 3所示摆线针轮传动模型导入仿真软件中，对针齿壳施加固定约束，在摆线轮端面施加额定扭矩T_g，其中单个摆线轮传递的扭矩T_g=0.55T₂。对摆线轮和针齿壳进行自动网格划分，对涉及接触作用的区域局部细分，采用Abaqus二次开发实现上述流程的自动化，如图 4所示。通过脚本接口建立图形用户界面GUI与内核的通信，提取摆线轮在扭矩T_g作用下当前角度的角位移β_H。

将输入轴与针齿壳固定，在输出轴上施加额定扭矩T₂，各弹性元件将引起相应的弹性转角θ_i如表 1所示，RV减速器的总刚度K′^[7]为

$ {K^\prime } = \frac{{{T_2}}}{{\sum\limits_{i = 1}^6 {{\theta _i}} }}。$

(8)

表中渐开线传动机构的位移$ {\delta _{{\rm{sp}}}} = \frac{{{F_{\rm{t}}}}}{{{c_{\rm{r}}}\cos \alpha }} $，为渐开线齿轮副的啮合角；转角$ {\theta _{{\rm{sp}}}} = \frac{{{\delta _{{\rm{sp}}}}\cos \alpha }}{{{r_{\rm{s}}}}} = \frac{{{F_{\rm{t}}}}}{{{c_{\rm{r}}}{r_{\rm{s}}}}} $r_s是节圆半径；z_b为针齿数；转臂轴承的位移δ_cb=R/K_r；支承轴承的位移$ {\delta _{{\rm{sb}}}} = \left| {R'} \right|/{K'_{\rm{r}}};i_{65}^1 = - i_{15}^6/i_{16}^5, i_{15}^6 = {z_{\rm{b}}}{z_2}/{z_1};i_{16}^5 $为当针齿固定时，输入轴相对于输出轴的转速比。

1.1.3 传动效率

RV减速器的主要效率损失来自于齿轮啮合摩擦损失和轴承摩擦损失^[15]，表示为

$ \eta = {\eta _{16}}{\kern 1pt} {\kern 1pt} {\eta _{\rm{B}}}{\rm{ }}。$

(9)

式中：η_B为轴承总效率，η₁₆为封闭差动齿轮传动效率，

$ {\eta _{16}} = \frac{{\left( {\frac{{{z_{\rm{b}}}}}{{{z_{\rm{g}}}}} - 1} \right)\left( {\frac{{{z_{\rm{b}}}}}{{{z_{\rm{g}}}}} - \eta _6^{6,2} + \frac{{{z_2}}}{{{z_1}}}\frac{{{z_{\rm{b}}}}}{{{z_{\rm{g}}}}}\eta _1^6} \right)}}{{\left( {\frac{{{z_{\rm{b}}}}}{{{z_{\rm{g}}}}} - \eta _6^{6,2}} \right)\left( {\frac{{{z_{\rm{b}}}}}{{{z_{\rm{g}}}}} - 1 + \left( {\frac{{{z_2}{z_{\rm{b}}}}}{{{z_1}{z_{\rm{g}}}}}} \right)} \right)}}。$

式中：η₁⁶为渐开线齿轮啮合效率，

$ \eta _1^6 = 1 - {\rm{ \mathsf{ π} }}{\kern 1pt} {\mu ^\prime }\left( {\frac{1}{{{z_1}}} + \frac{1}{{{z_2}}}} \right)\left( {\varepsilon _1^2 + \varepsilon _2^2 + 1 - {\varepsilon _1} - {\varepsilon _2}} \right); $

η₆^{6, 2}为摆线针轮啮合效率，

$ \eta _6^{6,2} = 1 - \frac{{fC}}{{{K_1}}}\left( {1 - \frac{{{d_{\rm{z}}}}}{{{D_{\rm{z}}}}}} \right)。$

式中：f为滑动摩擦因数；滑动特性系数$ C = \frac{4}{{\rm{ \mathsf{ π} }}}{z_{\rm{g}}} $。

1.2 结构参数的约束条件

1) 考虑相邻行星轮齿顶不干涉(g₄)、行星齿轮装配条件(h₁)、弯曲接触强度(g₅, g₆)等，如表 2所示。

表中齿数比u=z₂/z₁；齿顶高系数h_a^*=1；K为载荷系数；σ_F为齿根弯曲应力；[σ_F]为许用齿根弯曲应力；σ_H为齿面接触应力；[σ_H]为许用齿面接触应力。齿廓系数Y_Fa、应力修正系数Y_Sa、节点区域系数Z_H、弹性系数Z_E通过调用齿轮校核子程序计算。

2) 考虑圆柱滚子安全接触应力(g₈)^[3]、转臂轴承几何约束(g₁₀)、轴承寿命(g₁₁)等，如表 3所示。

表中的Q_max=4.08×10^-3|R|/Z；s为安全系数；[σ_c]为许用压应力；D_o为行星轮传动轴直径；e为偏心距；n_b为转臂轴承内外圈相对转速；寿命指数ε=10/3；p为平均当量动载荷；基本额定动载荷C_d为

$ {{C_{\rm{d}}} = {f_{\rm{c}}}{L^{\frac{7}{9}}}{Z^{\frac{3}{4}}}D_{\rm{r}}^{\frac{{29}}{{27}}}}。$

平均当量动载荷p为

$ {p = |R|\left[ {{{\left( {\sqrt {{{\left| {{R_x}} \right|}^2} + {{\left| {{R_y}} \right|}^2}} /|R| - 0.5} \right)}^2} + 0.75} \right]}。$

支承轴承的约束条件相似，故不再赘述。

3) 考虑针齿分布圆直径(g₁₂)、摆线轮不根切条件(g₁₆)、摆线针轮接触强度(g₁₇)等，如表 4所示。

表中$ {K_2} = {D_{\rm{z}}}\sin \left\{ {\frac{{\rm{ \mathsf{ π} }}}{{{z_{\rm{g}}}}}} \right\}/{d_{\rm{z}}} $；Y_1max是关于K₁、K₂和z_g的系数；K₂为针径系数。

4) 考虑两级传动的尺寸均衡(g₁₈)，并保证最大输出转矩(g₁₉)，如表 5所示。

2 多目标优化的理论基础 2.1 MP-NSGA-Ⅱ算法

ridPSO和ridDE^[10-11]是适用于单目标MINLP问题的主流算法，不具备多目标寻优能力。笔者结合多目标进化算法NSGA-Ⅱ提出MP-NSGA-Ⅱ。该算法继承了NSGA-Ⅱ解集分布性良好、时间复杂度低、收敛速度快等优点, 并且扩展了处理实数、整数及离散变量的能力。

求解优化问题的初始步骤是根据自变量范围矩阵随机生成编码种群，设计变量中$ {\left\{ {b, {D_{\rm{z}}}, {d_{\rm{z}}}, B, {K_1}, {D_{\rm{m}}}, {{D'}_{\rm{m}}}, {D_{\rm{r}}}, {{D'}_{\rm{r}}}, L, L'} \right\}^{\rm{T}}} $为连续变量，$ {\left\{ {{z_1}, {z_2}, Z, Z'} \right\}^{\rm{T}}} $为整数变量，$ {\left\{ {{z_{\rm{g}}}, m} \right\}^{\rm{T}}} $为离散变量。在进化算法库Geatpy的基础上改进，如图 5所示。

处理混合种群的具体步骤为：

1) 产生初始种群：利用设计变量的取值范围生成自变量范围矩阵，根据前m₁个连续变量和后m₂个离散变量在矩阵中的位置将其分割为2个子矩阵M_r和M_i。针对子矩阵M_r，使用rand函数生成一个实数值的初始种群P_r。针对子矩阵M_i，使用rand函数和四舍五入法生成一个十进制整数的初始种群P_i。

2) 混合种群交叉：种群P_r和P_i均是行数为种群规模，列数为设计变量数的矩阵，将P_i转化为浮点数矩阵，两矩阵水平合并，生成混合种群P_m, 采用两点交叉实现个体间染色体的重组。

3)子种群变异：将混合种群分割为子种群P_r和P_i，P_i转化为整数矩阵。采用实数值高斯变异算子对P_r进行变异，整数值变异算子实现矩阵P_i中个体突变，并将两矩阵合成混合种群P_m。

由于截断法^[9]只能处理连续的整数变量，为处理MIP(mixed integer programming)问题中的离散变量，提出一种基于数组索引的离散变量通用编码方案如图 6。将离散变量Z_d的n_d个可取值设为数组Z_arr, 数组索引编码为取值(0, 1, …, n_d-1)的整数子种群Z_num。在评估函数值和限制条件的阶段，根据数组Z_arr和数组索引种群Z_num解码出离散子种群Z_dis。

为增强种群的分布性并降低计算代价，引入考虑拥挤距离的非支配排序。由于计算欧氏距离的效率较低，采用差分计算目标值的偏移量占比代表拥挤距离。根据混合种群个体目标值由小到大排序，两相邻个体x_i和x_j间的拥挤距离定义为

$ {d_{\rm{c}}} = \frac{{f({x_j}) - f({x_i})}}{{f{{(x)}_{\max }} - f{{(x)}_{{\rm{min}}}}}} $

(10)

式中：f(x)是个体x的目标值，f(x)_max和f(x)_min分别是混合种群中最大和最小的目标值。根据拥挤距离更新个体的适应度V_fit，用于选择下一代个体。

2.2 双加点准则Kriging代理

利用拉丁超立方试验设计在设计变量空间内采样，得到Kriging代理模型的初始样本点，样本点满足使下式取得最小值：

$ \sum\nolimits_{i = 1}^{N - 1} {\sum\nolimits_{j = i + 1}^{N - 1} {\frac{1}{{{{\left\| {{x_i} - {x_j}} \right\|}^2}}}} } , $

(11)

式中：$ \left\| {{x_i} - {x_j}} \right\| $为两样本点间距离，N为样本点的数量。

Kriging代理模型定义了设计变量x与预测值y的关系，表达式^[8]为

$ y(x) = F(\beta ,x) + z(x), $

(12)

式中：F(β, x)为设计变量空间的全局模型，z(x)为按N(0, σ₂)随机分布的局部偏差, z(x)的统计特征为

$ \left\{ {\begin{array}{*{20}{l}} {E[z(x)] = 0,}\\ {{\mathop{\rm Var}\nolimits} [z(x)] = \sigma _{\rm{z}}^2,}\\ {{\mathop{\rm cov}} \left( {z\left( {{x_i}} \right),z\left( {{x_j}} \right)} \right) = \sigma _{\rm{z}}^2R\left( {{x_i},{x_j}} \right)}。\end{array}} \right. $

(13)

R(x_i, x_j)是用于刻画样本点x_i和x_j间关联程度的相关模型，通常采用高斯相关模型

$ R\left( {{x_i},{x_j}} \right) = \exp \left( { - \sum\limits_{k = 1}^{{n_{\rm{v}}}} {{\mathit{\boldsymbol{\theta }}_k}} {{\left| {x_i^k - x_j^k} \right|}^2}} \right), $

(14)

式中：n_v是设计变量的维数；θ是待定的相关参数向量；x_i^k、x_j^k和θ_k分别是x_i、x_j和θ的第k个分量。

利用线性加权插值方法得到Kriging模型在预测点x处的响应值和预测方差：

$ {\hat y(x) = F(\beta ,x) + {\mathit{\boldsymbol{r}}^{\rm{T}}}(x){\mathit{\boldsymbol{R}}^{ - 1}}(\mathit{\boldsymbol{g}} - \hat \beta F),} $

(15)

$ {{{\hat e}^2}(x) = {\sigma ^2}\left[ {1 - {\mathit{\boldsymbol{r}}^{\rm{T}}}\mathit{\boldsymbol{Rr}} + \left[ {\frac{{{{\left( {1 - {\mathit{\boldsymbol{q}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{r}}} \right)}^2}}}{{{\mathit{\boldsymbol{q}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{q}}}}} \right]} \right],} $

(16)

式中：R是相关模型矩阵；r(x)是x点与样本点间的相关模型向量；g是样本点响应的向量；q是元素均为1且个数为n_v的单位列向量。

在优化过程中添加样本点可有效提高代理模型的精度，通常选择在期望提高(EI)或均方误差(MSE)较大处加点。为有效利用进化算法优化过程中的信息，将Pareto最优集引入加点准则。在迭代过程中根据Pareto集的拥挤距离选择分散样本点更新Kriging模型。单次迭代中均方误差较大处加点数S_E和Pareto最优解处加点数S_P分别为

$ {S_{\rm{E}}} = \left[ {\frac{{g - {g_{\rm{c}}}}}{g} \cdot {C_{{\rm{E}}1}} + {C_{{\rm{T}}1}}} \right],{S_{\rm{P}}} = \left[ {\frac{{{g_{\rm{c}}}}}{g} \cdot {C_{{\rm{E}}2}} + {C_{{\rm{T}}2}}} \right], $

(17)

式中：g是迭代总次数；g_c是当前迭代数；C_E和C_T分别是调整系数和最小加点数。

2.3 熵权法Pareto选优

熵权法是一种客观赋权方法，利用决策指标的熵计算熵权值，在i个性能指标X₁, X₂, …, X_i，j个Pareto最优解的评价问题中，第p个性能指标X_p={x₁, x₂, …, x_j}，X_p的熵H_p定义为

$ {H_p} = - {k_j}\sum\nolimits_{q = 1}^j {{f_{pq}}} \ln {f_{pq}}(p = 1,2, \cdots ,i), $

(18)

式中：$ {f_{pq}} = \frac{{{r_{pq}}}}{{\sum\nolimits_{q = 1}^j {{r_{pq}}} }}, {k_j} = \frac{1}{{\ln j}} $。

归一化后的指标$ {r_{pq}} = \frac{{{X_{pq}} - \min \left( {{X_p}} \right)}}{{\max \left( {{X_p}} \right) - \min \left( {{X_p}} \right)}} $。指标的熵越大，熵权值越小，熵权值w_p定义为

$ {w_p} = \frac{{1 - {H_p}}}{{i - \sum\limits_{p = 1}^i {{H_p}} }}。$

(19)

3 优化方法的验证及实现 3.1 MP-NSGA-Ⅱ评估

为评估MP-NSGA-Ⅱ算法的有效性，对文献[11]中的14个MINLP问题进行仿真，已知最优源自MDE、MDELS、MDEIHS、ridPSO和ridDE。

表 6均为最小化目标问题，MP-NSGA-Ⅱ的种群规模为1 000。为体现算法的收敛特性，进化代数设为10 000。MP-NSGA-Ⅱ单次运行最优解与已知最优差距小于0.1%，并改进了3个已知最优解：P₅(x=1.374 823 1, y=1)，P₉(x=[27, 27, 27], y=[88, 44])，P₁₂(x=[0.902 19, 0.887 75, 0.949 18, 0.848 72], y=[5, 5, 4, 6])。由于P₇中已知最优的自变量x₂=0时目标函数分母为0，原解不成立。

3.2 Kriging近似模型

选择影响摆线针轮间刚度K″的{z_g, D_z, d_z, K₁, B}^T为设计变量，利用拉丁超立方试验设计(LHS)在5个设计变量的取值范围内随机选取32个样本点。利用Abaqus计算各样本点对应的仿真值，采用pyKriging模块中的train()和predict()函数构建变量与预测值间的Kriging代理模型。图 7是以K₁和z_g为自变量的Kriging可视化模型，可见Kriging响应面提供了梯度信息，降低了目标函数的非线性，避免了陷入局部最优。

采用复相关系数R²检验拟合模型的精度为

$ {R^2} = \frac{{\sum\limits_{u = 1}^v {{{\left( {{{\hat y}_u} - \bar y} \right)}^2}} }}{{\sum\limits_{u = 1}^v {{{\left( {{y_u} - \bar y} \right)}^2}} }}, $

(20)

式中：v是检验模型的样本点数量；$ {\hat y_u} $和y_u分别是第u个样本点的预测值与仿真值；y是仿真值的平均值。

复相关系数越接近1，近似模型的精度越高。选择16个样本点检验Kriging模型的精度，如图 8所示。检验得复相关系数R²为0.920 8，说明该代理模型的精度能够满足结构优化的需要。

3.3 结构优化子程序

RV减速器的优化目标为{η, V, K′}^T，约束函数为几何及应力约束，优化模型为

$ \left\{ {\begin{array}{*{20}{c}} {\min \left[ { - \eta (\mathit{\boldsymbol{X}},\mathit{\boldsymbol{C}}),V(\mathit{\boldsymbol{X}},\mathit{\boldsymbol{C}}), - {K^\prime }(\mathit{\boldsymbol{X}},\mathit{\boldsymbol{C}})} \right],}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\quad g(\mathit{\boldsymbol{X}},\mathit{\boldsymbol{C}}) \le 0,}\\ {h(\mathit{\boldsymbol{X}},\mathit{\boldsymbol{C}}) = 0}。\end{array}} \right. $

(21)

式中：g(X, C)为不等式约束；h(X, C)为等式约束。

由上述的优化模型及算法原理，以PySide2为开发框架，编写一体化结构RV减速器设计软件。其优化子程序如图 9所示，由基本参数及设计变量设置、优化算法参数设置及优化结果输出组成。

根据RV减速器的实际工况给出RV减速器的动力学及结构基本参数，如图 9左上。根据RV减速器的设计要求及BAJ-25E减速器的设计参数，初算设计变量范围，如图 9左下。

设置MP-NSGA-Ⅱ算法的种群规模为3 000，遗传代数为1 000，代沟为0.5，交叉概率为90%，变异概率为10%。选择方式为轮盘赌选择，重组方式为两点交叉。Pareto前沿如图 10所示，可知优化目标间相互制约，需从解集中优选理想解。

生成RV减速器的结构及性能参数表如图 9右下所示。设定传动效率及扭转刚度下限和体积上限，以缩减Pareto最优解集的规模。

4 结果与分析

各优化目标通过结构参数相互耦合，利用两目标的Pareto前沿分析优化目标间的耦合关系，为后续的Pareto选优提供依据，如图 11所示。

从图 11(a)和(b)可看出扭转刚度和体积变化时，传动效率维持在85.2%~85.6%，与其余优化目标的耦合关系不显著。11(c)表明体积与扭转刚度成显著正相关。因此，主要考虑体积和扭转刚度的设计要求，初步筛选5个设计方案如表 7所示。

采用熵权法，将性能指标决策矩阵归一化，求解各指标的信息熵H_p，并计算各目标熵权值：

$ {w_p} = [0.293{\kern 1pt} {\kern 1pt} {\kern 1pt} 2,0.318{\kern 1pt} {\kern 1pt} {\kern 1pt} 6,0.388{\kern 1pt} {\kern 1pt} {\kern 1pt} 2],p = 1,2,3。$

由熵权值w_p及归一化指标r_pq对各方案评分，为

$ {Z_q} = {r_{1q}}{w_1} - {r_{2q}}{w_2} + {r_{3q}}{w_3}\quad (q = 1,2,3,4,5)。$

(22)

根据Z_q排序，选择理想方案2。BAJ-25E的结构参数与优化后的RV减速器参数如表 8所示。

由表 8知，经MP-NSGA-Ⅱ算法优化可直接获得符合约束条件的优化值，无需圆整处理。与初始值相比，优化解的传动效率提升了1.24%，体积减小了1.69%，扭转刚度增大了53.83%。

在制造过程中，尺寸参数可能出现1%左右的偏差^[3]，对性能指标的敏感性分析有助于指导实际生产中的误差控制。研究连续变量变化±1%对优化后的减速器扭转刚度的影响，如表 8。可见D_z、K₁、L和L′对扭转刚度的影响均超过0.5%，其它参数波动的影响较小，均低于0.2%。

5 结论

1) 分析了影响RV减速器性能的结构参数，综合17个设计变量，3个目标函数和21个约束条件建立了结构优化的数学模型。

2) 提出MP-NSGA-Ⅱ算法，利用MINLP基准问题测试了改进算法的性能，并改进了3个已知最优解，证明了该算法的有效性。

3) 结合Kriging与MP-NSGA-Ⅱ得到Pareto前沿，利用熵权法完成Pareto选优，有效提高了RV减速器的综合性能。