CN107273609A

CN107273609A - 一种基于Kriging模型齿轮传动可靠性评估方法

Info

Publication number: CN107273609A
Application number: CN201710446255.4A
Authority: CN
Inventors: 黄洪钟; 钱正坤; 黄土地; 李彦锋; 付国忠; 米金华; 郭骏宇
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2017-06-14
Filing date: 2017-06-14
Publication date: 2017-10-20

Abstract

本发明公开一种基于Kriging模型齿轮传动可靠性评估方法，应用于可靠性评估领域。本申请首先根据初始的样本建立Kriging模型，然后根据建立的Kriging模型，判断变异系数是否满足要求，通过增加样本点的方式，重新计算新的变异系数；大大降低了传统方法中用有限元方法求应力的计算量；本申请采用蒙特卡罗方法计算结构可靠度，同时结合学习函数以较少的样本点得到高精度的代理模型，从而减少了试验次数。

Description

一种基于Kriging模型齿轮传动可靠性评估方法

技术领域

本发明属于齿轮传动可靠性评估领域，特别涉及一种齿轮传动可靠性评估技术。

背景技术

齿轮传动广泛应用于工业生产和日常生活中，如汽车、风力发电、机床、航空航天等领域。在实际运用过程中，由于各种不同的外部因素影响，使得齿轮失效形式多种多样，如齿面点蚀、齿面磨损、轮齿折断、齿面胶合和塑性变形等。当前机械产品朝着高可靠性方向发展，如果不针对这些失效机理开展可靠性研究，产品势必会失去竞争力。

齿轮传动可靠度计算主要是通过构建结构功能函数，求解结构强度和应力(如接触应力和弯曲应力等)，然后根据一次可靠度、二次可靠度或蒙特卡罗方法求解结构可靠度。为了得到高精度的可靠度，一般采用蒙特卡罗方法进行可靠度计算。但结构强度和应力均是随机变量，且服从不同分布类型。影响齿轮应力的因素较多，主要有三个方面：内部激励、外部激励和系统参数不确定。这就导致使用蒙特卡罗方法计算可靠度时计算量较大，有时无法进行计算(如有些应力需要通过试验方式获得)。

当前，为解决上述问题通常采用Kriging代理模型近似结构功能函数，然后通过蒙特卡罗方法求解可靠度。但传统的采样方法很难得到较好的样本点，使得构建的Kriging模型精度较差。通过增加样本点可以提高代理模型的精度，但这样会增加试验成本。可靠度的计算精度很大程度上取决于代理模型的精度，因此在保证较低的试验成本情况下获得高精度的代理模型，对齿轮传动可靠度的计算十分重要。目前，国内对基于Kriging模型齿轮传动可靠性评估研究较少。

发明内容

本发明为解决上述技术问题，提出了一种基于Kriging模型齿轮传动可靠性评估方法。

本发明采用的技术方案是：一种基于Kriging模型齿轮传动可靠性评估方法，包括：

S1、根据齿轮传动结构功能函数，确定随机变量以及分布类型；

S2、对每个随机变量进行蒙特卡罗随机抽样，生成N个符合该变量分布类型的随机数，所述N个随机数作为对应随机变量的候选样本点；

每个随机变量对应的候选样本点构成候选样本点集合；

S3、采用拉丁超立方抽样法从每个随机变量的候选样本点中抽取相同数量的初始训练样本点，根据各个随机变量对应的初始训练样本点得到初始训练样本矩阵；

S4、将初始训练样本矩阵带入结构功能函数，得到真实响应值矩阵；

S5、根据当前的总的初始样本以及真实响应值矩阵，构造初始Kriging模型；

S6、根据步骤S5的初始Kriging模型，计算每一个随机变量对应的候选样本点的学习函数值，并找出该随机变量对应的最小学习函数值；若存在某个最小学习函数值小于2；则将各随机变量对应的候选样本点中最小学习函数值对应的候选样本点，作为该随机变量的最佳样本点；并计算各最佳样本点的真实响应值；然后将各随机变量的最佳样本点加入到初始训练样本中，得到更新后的初始训练样本；并且将各最佳样本点对应的真实响应值加入到真实响应值矩阵中，得到更新后的真实响应值矩阵；然后返回步骤S5；否则执行步骤S7；

S7、根据当前Kriging模型以及候选样本点集合，计算构件的失效率和变异系数，若当前的变异系数小于或等于预设阈值，则输出当前失效率和变异系数；否则执行步骤S8；

S8、重新对每个随机变量进行蒙特卡罗随机抽样，生成M个符合该变量分布类型的随机数，根据各个变量对应的M个随机数，替换当前的候选样本点集合；然后返回步骤S7；

所述M＝N+q*E；E为常数，q为迭代次数。

进一步地，所述步骤S6具体为：

S61、根据步骤S5构造的初始Kriging模型，计算候选样本点集合中所有候选样本的预测值μ_y(x_i)和标准差σ_y(x_i)；i为候选样本点集合中候选样本点的序号；

所述预测值μ_y(x_i)计算式为：

μ_y(x_i)＝f^Tβ^*+r^TR^-1(Y-Fβ^*)

其中，β^*为中间变量；F表示多项式函数矩阵；f表示多项式函数矩阵F中的元素；Y表示函数真实响应值；上标-1表示求逆矩阵；R表示候选样本点集合中任意两个候选样本点的相关性；

所述标准差σ_y(x_i)的计算式为：

其中，r表示检测点与候选样本点之间的相关性；u、σ均为中间变量；

S62、根据步骤S61得到预测值μ_y(x_i)和标准差σ_y(x_i)，计算每一个随机变量对应的候选样本点的学习函数值U(x_i)；

S63、根据步骤S62计算得到的U(x_i)，找出该随机变量对应的最小学习函数值；并判断是否存在最小的学习函数值小于2，若是则将各随机变量对应的候选样本点中最小学习函数值对应的候选样本点，作为该随机变量的最佳样本点；并计算各最佳样本点的真实响应值；然后将各随机变量的最佳样本点加入到初始训练样本中，得到更新后的初始训练样本；并且将各最佳样本点对应的真实响应值加入到真实响应值矩阵中，得到更新后的真实响应值矩阵，并返回步骤S5；否则执行步骤S7。

进一步地，步骤S7所述失效率和变异系数的计算公式如下：

其中，表示失效率，表示变异系数，I_f(y(x_i))为示性函数，当y(x_i)≤0时I_f(y(x_i))＝1；当y(x_i)＞0时I_f(y(x_i))＝0；N_y≤0为落入失效区域内样本点总数；N_MC为当前候选样本点集合中候选样本点的总数量。

本发明的有益效果：本发明的方法，首先根据初始的样本建立Kriging模型，然后根据建立的Kriging模型，判断变异系数是否满足要求，通过增加样本点的方式，重新计算新的变异系数；本申请采用蒙特卡罗方法计算结构可靠度，同时结合学习函数以较少的样本点得到高精度的代理模型，减少了试验次数，大大降低了计算量。

附图说明

图1为本发明的方案流程图。

具体实施方式

为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。

本发明的技术方案如图1所示；一种基于Kriging模型齿轮传动可靠性评估方法，包括以下步骤：

S1：建立齿轮传动结构功能函数，确定函数中的随机变量S＝[x₁,x₂,…,x_n]^T和分布类型。

结构功能函数为：

G(x)＝D₀-D(x) (1)

其中，D₀结构强度，D(x)结构应力；这里的结构功能函数也可采用现有的Kriging模型来近似得到，以简化运算。

各随机变量分布类型通过数据收集和假设检验得到。

S2：对每个随机变量进行蒙特卡罗随机抽样，生成N个符合该变量分布类型的随机数。由于这些样本点用来更新试验样本点，故称这些样本点为候选样本总体。

即在本步骤中，每个随机变量对应了一个包含N个样本点的候选样本。本实施例首次对每个随机变量进行蒙特卡罗随机抽样时，生成1000个随机数。

S3：采用拉丁超立方抽样法对每个随机变量抽取相同数量的初始训练样本点(Design of Experiment，DOE)构成该随机变量的初始训练样本，由各个随机变量对应的初始训练样本构成初始训练样本矩阵，用S＝[s₁,s₂,…,s_m]^T表示；

比如现有3个随机变量，第一个随机变量蒙特卡罗抽样随机数为：1 5 8 9共四个；第二个随机变量蒙特卡罗抽样随机数为：5 6 8 0共四个；第三个随机变量蒙特卡罗抽样随机数为：4 5 8 0共四个；那么S矩阵为

其中，s₁表示第一个随机变量对应的初始训练样本；同理，s₂表示第二个随机变量对应的初始训练样本；同理，s_m表示第m个随机变量对应的初始训练样本；上标T表示转置。

所述的随机变量的初始训练样本，即从该随机变量的候选样本中，选择一部分样本点作为初始样本点组成该随机变量的初始训练样本，一般当N＝1000时，会选择12个左右的样本点作为初始训练样本。

S4、将初始训练样本矩阵代入结构功能函数求得函数真实响应值Y_DOE，这里的Y_DOE通过有限元方法得到。

S5：根据S和Y_DOE构造初始Kriging模型。

根据S和Y_DOE构造的初始Kriging模型为：

其中，β表示线性回归系数；f(x)表示变量X的多项式函数；x表示随机变量X的具体数据；n表示f_i(x)的个数；z(x)是服从正态分布的随机变量；

z(x)的协方差矩阵为：

函数为样本空间中任意两个样本点ω和x的相关函数，θ为相关函数的参数，函数具体表达式为：

若初始样本为S＝[s₁,s₂,…,s_m]^T，则多项式函数矩阵F可表示为：

F＝[f(s₁),f(s₂),…,f(s_m)]^T (6)

初始样本的相关性R_ij为：

若x为检测点，S为样本点，则它们之间的相关性可用r(x)表示：

初始样本响应值Y为：

Y＝[y₁,y₂,…,y_m] (9)

S6、根据步骤S5的初始Kriging模型，计算候选样本中所有点的预测值μ_y(x_i)和标准差σ_y(x_i)，并计算学习函数U(x_i)找出最小值minU(x)。

预测值μ_y(x_i)和标准差σ_y(x_i)的计算过程为：

若用线性组合来近似一个检测点的响应，则

近似值和真实响应值之间的偏差为：

其中，Z＝[z₁,z₂,…,z_m]^T为初始样本点误差。

为了使预测无偏，则F^Tc-f(x)＝0。偏差均方差可表示为：

采用拉格朗日方法，求最小均方差可得c，约束函数为F^Tc-f(x)＝0。

L(c,λ)＝σ²(1+c^TRc-2c^Tr)-λ^T(F^Tc-f) (13)

参数c的梯度为：

L'_c(c,λ)＝2σ²(Rc-r)-Fλ (14)

由优化一阶必要条件，可得：

把代入式(15)，得：

将式(16)代入式(10)，可得：

由得：

β^*＝(F^TR^-1F)^-1F^TR^-1Y (18)

把式(18)代入式(17)，得到Kriging模型预测值为：

Kriging模型均方误差为：

其中，r表示检测点与候选样本点之间的相关性，具体计算式参考式(8)；u、σ均为中间变量；且

u＝F^TR^-1r-f，

学习函数表达式如下：

其中，μ_y(x)表示Kriging代理模型预测均值；σ_y(x)表示Kriging代理模型预测标准差。

当存在minU(x)＜2时，得到各随机变量学习函数最小值minU(x)对应的该随机变量候选样本中的最佳样本点x^*，并求得该最佳样本点的真实响应值y(x^*)，将该最佳样本点加入到初始训练样本向量中，得到更新后的初始训练样本向量，并将该最佳样本点对应的真实响应值加入到真实响应值矩阵中，更新真实响应值矩阵，然后返回步骤S5；否则进行步骤S7。

S7、根据当前Kriging模型以及候选样本点集合，计算构件的失效率和变异系数若当前的变异系数小于或等于预设阈值δ，则输出当前失效率和变异系数；否则执行步骤S8；

和的计算公式如下：

其中，I_f(y(x_i))为示性函数，当y(x_i)≤0时I_f(y(x_i))＝1；当y(x_i)＞0时I_f(y(x_i))＝0；N_y≤0为落入失效区域内样本点总数；N_MC为当前候选样本点集合中候选样本点的总数量；MC表示蒙特卡洛抽样。

一般情况下，取δ＝0.03。

所述M＝N+q*E；E为常数，q为迭代次数。

当首次生成的N＝1000个样本点计算得到的变异系数不满足要求，可以采用每次增加500个样本的方式，即第一次迭代不满足要求时，对每个随机变量生成1500个蒙特卡罗随机样本，用于计算结构可靠度；来使得变异系数满足要求。第二次迭代不满足要求则生成2000个蒙特卡罗随机样本，用于计算结构可靠度；直至当前的变异系数小于或等于预设阈值δ。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims

1.一种基于Kriging模型齿轮传动可靠性评估方法，其特征在于，包括：

每个随机变量对应的候选样本点构成候选样本点集合；

S6、根据步骤S5的初始Kriging模型，计算每一个随机变量对应的候选样本点的学习函数值，并找出该随机变量对应的最小学习函数值；若存在某个最小学习函数值小于2，则将各随机变量对应的候选样本点中最小学习函数值对应的候选样本点，作为该随机变量的最佳样本点；并计算各最佳样本点的真实响应值；然后将各随机变量的最佳样本点加入到初始训练样本中，得到更新后的初始训练样本；并且将各最佳样本点对应的真实响应值加入到真实响应值矩阵中，得到更新后的真实响应值矩阵；然后返回步骤S5；否则执行步骤S7；

所述M＝N+q*E；E为常数，q为迭代次数。

2.根据权利要求1所述的一种基于Kriging模型齿轮传动可靠性评估方法，其特征在于，所述步骤S6具体为：

所述预测值μ_y(x_i)计算式为：

μ_y(x_i)＝f^Tβ^*+r^TR^-1(Y-Fβ^*)

所述标准差σ_y(x_i)的计算式为：

<mrow> <msubsup> <mi>&sigma;</mi> <mi>y</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> <mo>&lsqb;</mo> <mn>1</mn> <mo>+</mo> <msup> <mi>u</mi> <mi>T</mi> </msup> <msup> <mrow> <mo>(</mo> <msup> <mi>F</mi> <mi>T</mi> </msup> <msup> <mi>R</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mi>F</mi> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mi>u</mi> <mo>-</mo> <msup> <mi>r</mi> <mi>T</mi> </msup> <msup> <mi>R</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mi>r</mi> <mo>&rsqb;</mo> </mrow>

<mrow> <mi>U</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mo>|</mo> <mfrac> <mrow> <msub> <mi>&mu;</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>&sigma;</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>|</mo> </mrow>

S63、根据步骤S62计算得到的U(x_i)，找出该随机变量对应的最小学习函数值；并判断是否存在最小的学习函数值大于2，若是则将各随机变量对应的候选样本点中最小学习函数值对应的候选样本点，作为该随机变量的最佳样本点；并计算各最佳样本点的真实响应值；然后将各随机变量的最佳样本点加入到初始训练样本中，得到更新后的初始训练样本；并且将各最佳样本点对应的真实响应值加入到真实响应值矩阵中，得到更新后的真实响应值矩阵，并返回步骤S5；否则执行步骤S7。

3.根据权利要求1所述的一种基于Kriging模型齿轮传动可靠性评估方法，其特征在于，步骤S7所述失效率和变异系数的计算公式如下：

<mrow> <msub> <mi>P</mi> <mi>f</mi> </msub> <mo>&ap;</mo> <msub> <mover> <mi>P</mi> <mo>^</mo> </mover> <mi>f</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <msub> <mi>N</mi> <mrow> <mi>M</mi> <mi>C</mi> </mrow> </msub> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>M</mi> <mi>C</mi> </mrow> </msub> </munderover> <msub> <mi>I</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mi>y</mi> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <msub> <mi>N</mi> <mrow> <mi>y</mi> <mo>&le;</mo> <mn>0</mn> </mrow> </msub> <msub> <mi>N</mi> <mrow> <mi>M</mi> <mi>C</mi> </mrow> </msub> </mfrac> </mrow>

<mrow> <msub> <mi>&delta;</mi> <msub> <mover> <mi>P</mi> <mo>^</mo> </mover> <mi>f</mi> </msub> </msub> <mo>=</mo> <msqrt> <mfrac> <mrow> <mn>1</mn> <mo>-</mo> <msub> <mover> <mi>P</mi> <mo>^</mo> </mover> <mi>f</mi> </msub> </mrow> <mrow> <mo>(</mo> <msub> <mi>M</mi> <mrow> <mi>M</mi> <mi>C</mi> </mrow> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> <msub> <mover> <mi>P</mi> <mo>^</mo> </mover> <mi>f</mi> </msub> </mrow> </mfrac> </msqrt> </mrow>