本发明属于数控机床可靠性评估技术领域，涉及一种考虑工况条件差异的数控机床可靠性评估的方法，包括以下步骤：1、采集工况信息、故障数据及运行数据，进行数据整理；2、考虑工况条件差异的数控机床变参数幂律分布可靠性建模；3、筛选工况因素；4、模型参数估计；5、考虑工况条件差异的数控机床混合变参数幂律分布可靠性建模；6、对混合变参数幂律分布进行假设检验。本发明利用混合变参数幂律分布对同一工况水平中多种工况条件下的数控机床可靠性函数进行拟合，同时，能够综合的考虑多种工况因素对可靠性水平评估的影响；相比于传统的可靠性建模与评估方法，本发明能够准确的对不同工况水平中多种工况条件下的数控机床的可靠性进行评估。
1.一种考虑工况差异的数控机床可靠性评估的方法，其特征在于，包括以下步骤：步骤1：采集工况信息、故障数据及运行数据，进行数据整理；
步骤2：考虑工况差异的数控机床变参数幂律分布可靠性建模：
步骤5：考虑工况差异的数控机床混合变参数幂律分布可靠性建模：
2.根据权利要求1所述的一种考虑工况差异的数控机床可靠性评估的方法，其特征在于：步骤1中所述的工况信息是指在数控机床在发生故障或截尾时所处的工况条件中的各个工况因素，包括：数控机床信息、工件信息、刀具信息、切削参数及环境因素；
所述故障数据记录时应明确每一个故障的故障现象、故障部位、故障原因及故障处理内容；
所述运行数据主要记录数控机床在实际运行过程中的运行时间、停机时间和维修时间数据。
3.根据权利要求1所述的一种考虑工况差异的数控机床可靠性评估的方法，其特征在于：所述数据整理是指：
假设共有K台数控机床，在数控机床用户现场可靠性试验过程中总共采集了n条故障数据，其中包括r条完整的故障数据及n-r条截尾数据，每台数控机床最多存在一个截尾数据，r＝n1+n2+…+nq+…+nK；n-r＝n1+n2+…+na+…+nK；
第j个故障发生时间为tj，则第q台数控机床的第j个故障间隔工作时间tq,j＝tj-tj-1；
ta为第a台数控机床的可靠性试验结束时间， 为第q台数控机床在可靠性试验周期内的最后一个故障发生的时间，当 不等于ta时，则第a台数控机床存在截尾数据，即na＝1，且第a台数控机床的截尾时间为采集的第q台数控机床的第j条故障数据对应的工况信息GKq,j＝[Z1bj Z2bj ... Zibj ... Zkbj]；
采集的第a台数控机床的截尾数据对应的工况信息GKa＝[Z1ba Z2ba ... Ziba ... Zkba]；
其中Zibj表示第j个故障间隔工作时间所对应第b种工况条件下的第i个工况因素的工况协变量；
Ziba表示第a台数控机床的截尾时间所对应第b种工况条件下的第i个工况因素的工况协变量；
将现场跟踪试验采集的工况信息进行整理，最终确定了数控机床的m种工况条件，及每种工况条件对应的k种工况因素；
采集的第b种工况条件对应的工况信息GKb＝[Z1b Z2b ... Zib ... Zkb]，其中Zib表示数控机床第b种工况条件下的第i个工况因素的工况协变量。
4.根据权利要求1、2或3所述的一种考虑工况差异的数控机床可靠性评估的方法，其特征在于：步骤2中所述考虑工况差异的数控机床变参数幂律分布可靠性建模，具体包括以下步骤：
假定数控机床的故障间隔工作时间t服从两参数的幂律分布，其故障率函数λ(t)为：其中：β为幂律分布的形状参数；
即而得到数控机床的可靠度函数R(t)为：概率密度函数f(t)为：
f(t)＝-dR(t)/dt (3)累积分布函数F(t)为：
考虑m种工况条件下k种工况因素对数控机床可靠性水平的影响，对幂律分布中的尺度参数η建立变参数幂律分布模型；
选择第1种工况条件作为基准工况条件，则数控机床在第b种工况条件下幂律分布的尺度参数ηb为：其中：第1种工况条件下第i个工况因素的工况协变量Zi1为基准工况条件下第i个工况因素的工况协变量；
Zib为第b种工况条件下第i个工况因素的工况协变量；
η1为数控机床在基准工况条件下幂律分布的尺度参数；
将式(6)代入式(1)得到第b种工况条件下故障率函数λb(t)：将式(7)代入式(2)中，得到数控机床在第b种工况条件下考虑k种工况因素的可靠度函数Rb(t)为：将式(8)代入式(3)中，得到数控机床在第b种工况条件下考虑k种工况因素的概率密度函数fb(t)为：将式(7)代入式(4)中，得到数控机床在第b种工况条件下考虑k种工况因素的累积分布函数Fb(t)为：将式(9)代入式(5)中，得到数控机床在第b种工况条件下考虑k种工况因素的数控机床的平均故障间隔工作时间MTBFb为：其中：Γ(1+1/(β+1))——Gamma函数。
5.根据权利要求4所述的一种考虑工况差异的数控机床可靠性评估的方法，其特征在于：步骤3中所述筛选工况因素，具体包括以下步骤：
所述比例模型为第b种工况条件与基准工况条件故障率函数的比例关系；
将基准工况条件下幂律分布的尺度参数η1代入式(1)得到基准工况条件下故障率函数λ1(t)，即而得到式(7)中第b种工况条件下故障率函数与基准工况条件下故障率函数的比例关系为：其中：Zi1为第1种工况条件下第i个工况因素的工况协变量；
Zib为第b种工况条件下第i个工况因素的工况协变量；
θi为第i种工况因素的比例系数，且θi＝ξi(1+β)；
3.2、采用偏似然函数估计比例系数：为了估计式(12)中的第i种工况因素的比例系数θi，使用偏似然函数L(θi)：其中：Zibj为第j个故障间隔工作时间所对应第b种工况条件的第i个工况因素的工况协变量；
将式(14)对比例系数θi求偏导数并令其等于0得：将采集的数据整理后代入式(15)，采用Newton-raphson数值算法进行求解，即可得出第i种工况因素的比例系数θi；
3.3、采用卡方检验法筛选工况因素：第i种工况因素的卡方检验统计量χ2(k)i为：其中：θi为第i种工况因素的比例系数；
为在所有工况条件下第i种工况因素的最小工况协变量；
Zibj为第j个故障间隔工作时间所对应第b种工况条件下的第i个工况因素的工况协变量。
6.根据权利要求5所述的一种考虑工况差异的数控机床可靠性评估的方法，其特征在于：步骤4中所述模型参数估计，具体包括以下步骤：在选择了合适的工况因素的基础上，将采集的数据整理后代入式(8)和式(9)采用极大似然估计方法估计在基准工况条件下的故障率函数参数尺度参数η1和形状参数β；极大似然估计函数L(η1,β)为：其中：tq,j为第q台数控机床的第j个故障间隔工作时间；
对式(18)的形状参数β和尺度参数η1分别求偏导数并令其等于零得：根据式(15)求得的第i种工况因素的比例系数θi，将ξi＝θi/(1+β)代入式(19)和式(20)，采用Newton-raphson数值算法得到基准工况条件下的故障率函数的尺度参数η1与形状参数β的估计值，同时得到第i种工况因素的系数ξi的值，再根据式(6)得到第b种工况条件下幂律分布的尺度参数ηb的估计值。
7.根据权利要求6所述的一种考虑工况差异的数控机床可靠性评估的方法，其特征在于：步骤5中所述考虑工况差异的数控机床混合变参数幂律分布可靠性建模，具体包括以下步骤：
为了建立考虑m种工况条件差异的数控机床混合变参数幂律分布可靠性模型，需依据斯特奇斯经验公式将m种工况条件划分成d组：d＝[1+3.32lg(m)] (21)其中：[]代表取整；
令δb为第b种工况条件下的工况分组系数，则工况分组系数δb为：其中：ξi为第i种工况因素的系数；
设第e组内有me种工况条件，且 根据工况分组组距Δδ的大小，得到第e组工况分组系数的区间为(δmax-(e-1)Δδ,δmax-eΔδ)，再根据第e组工况分组系数的区间得到第e组内的me种工况条件，即得到第e种工况水平下的me种工况条件，利用混合变参数幂律分布对同一工况水平下工况条件的数控机床可靠性函数进行拟合，同时，能够综合的考虑多种工况因素对可靠性水平评估的影响；
5.2建立混合变参数幂律分布可靠性模型：本发明中混合变参数幂律分布模型适用于数控机床存在一种或多种工况条件差异的情况；根据分组结果，对第e种工况水平下的me种工况条件建立混合变参数幂律分布模型；
令第b种工况条件下的混合权系数为pb，同时 且pb＞0,b∈me，则第e种工况水平对应的混合累积分布函数Fe(t)为：第e种工况水平对应的混合概率密度分布函数fe(t)为：第e种工况水平对应的混合可靠度函数Re(t)为：第e种工况水平对应的混合平均故障间隔工作时间MTBFe为：将式(10)、(9)、(8)和(11)代入式(24)、(25)、(26)和(27)得到考虑工况差异的数控机床混合变参数幂律分布累积分布函数Fe(t)为：概率密度分布函数fe(t)为：
设第e组内工况协变量Zib的标准差σib(Z)为：其中： 为第e种工况水平下的所有工况条件中第i种工况因素的工况协变量的平均值；
则混合变参数幂律分布的相对权系数Pb为：将相对权系数进行归一化处理，得到混合权系数pb为：
8.根据权利要求7所述的一种考虑工况差异的数控机床可靠性评估的方法，其特征在于：步骤6中所述假设检验是指：
采用Kolmogorov-Sminov假设检验方法，即K-S假设检验，对混合变参数幂律分布进行假设检验；
设第e个混合变参数幂律分布的K-S检验统计量De为:De＝max|Fn(t)-Fe(t)|,e＝1,2,…,d (35)其中：Fn(t)——样本经验分布累计分布函数；
其中，样本经验累计分布函数Fn(t)为：在给定K-S检验显著性水平α下的临界值Dn,α，当由试验数据计算得到De大于临界值Dn,α时，可判断原假设不成立；当由试验数据计算得到De小于临界值Dn,α时，则可判断原假设成立。
一种考虑工况差异的数控机床可靠性评估的方法
技术领域
[0001] 本发明属于数控机床可靠性评估技术领域，涉及一种考虑工况差异的数控机床可靠性评估的方法，更确切地说，本发明提出了一种考虑工况差异的数控机床混合变参数幂律分布的可靠性建模与评估的方法。
背景技术
[0002] 数控机床作为装备制造业的工作母机，其用户遍及各个行业和地区，不同行业由于加工产品类型的不同使得数控机床加工过程中所处工况条件中的工况因素(如切削力、扭矩、切削速度等)存在很大的差异。而传统的数控机床可靠性建模与评估方法主要是基于故障数据的建模，没有考虑数控机床所处工况条件的差异对其可靠性的影响，这样就会降低数控机床可靠性建模与评估的准确性。研究表明：同型号同批次的数控机床由于所处的工况条件的不同，呈现出的可靠性水平差异很大。目前，在研究数控机床故障分布规律时，常用的故障分布模型是两参数或三参数的威布尔分布模型，其前提条件是故障间隔工作时间独立同分布，然而，数控机床是集机、电、液、气于一体的复杂可修系统，在数控机床发生故障之后，往往很难做到修复如新。因此，不应该将它们拟合成为一个威布尔分布模型。
发明内容
[0003] 为了更好的对数控机床进行可靠性建模与评估，本发明利用大量的用户现场可靠性试验数据，提出了一种考虑工况差异的数控机床混合变参数幂律分布的可靠性建模与评估的方法。本发明是采用如下技术方案实现的，结合附图说明如下：
[0004] 一种考虑工况差异的数控机床可靠性评估的方法，包括以下步骤：
[0005] 步骤1：采集工况信息、故障数据及运行数据，进行数据整理；
[0006] 步骤2：考虑工况差异的数控机床变参数幂律分布可靠性建模；
[0007] 步骤3：筛选工况因素；
[0008] 步骤4：模型参数估计；
[0009] 步骤5：考虑工况差异的数控机床混合变参数幂律分布可靠性建模；
[0010] 步骤6：对混合变参数幂律分布进行假设检验。
[0011] 技术方案步骤1中所述的工况信息是指在数控机床在发生故障或截尾时所处的工况条件中的各个工况因素，包括：数控机床信息、工件信息、刀具信息、切削参数及环境因素；
[0012] 所述故障数据是指每一个故障的故障现象、故障部位、故障原因及故障处理内容；
[0013] 所述运行数据是指数控机床在实际运行过程中的运行时间、停机时间和维修时间数据。
[0014] 技术方案中所述数据整理是指：
[0015] 假设共有K台数控机床，在数控机床用户现场可靠性试验过程中总共采集了n条故障数据，其中包括r条完整的故障数据及n-r条截尾数据，每台数控机床最多存在一个截尾数据， r＝n1+n2+...+nq+...+nK；n-r＝n1+n2+...+na+...+nK；
[0016] nq表示第q台数控机床共发生的故障数量；
[0017] na＝1表示第a台数控机床存在截尾数据；
[0018] na＝0表示第a台数控机床不存在截尾数据；
[0019] 第j个故障发生时间为tj，则第q台数控机床的第j个故障间隔工作时间tq，j＝tj-tj-1；
[0020] ta为第a台数控机床的可靠性试验结束时间， 为第q台数控机床在可靠性试验周期内的最后一个故障发生的时间，当 不等于ta时，则第a台数控机床存在截尾数据，即na＝1，且第a台数控机床的截尾时间为
[0021] 采集的第q台数控机床第j条故障数据对应的工况信息GKq，j＝ [Z1bj Z2bj ... Zibj ... Zkbj]；
[0022] 采集的第a台数控机床的截尾数据对应的工况信息GKa＝[Z1ba Z2ba ... Ziba ... Zkba]；
[0023] 其中Zibj表示第j个故障间隔工作时间所对应第b种工况条件下的第i个工况因素的工况协变量；
[0024] Ziba表示第a台数控机床的截尾时间所对应第b种工况条件下的第i个工况因素的工况协变量；
[0025] 将现场跟踪试验采集的工况信息进行整理，最终确定了数控机床的m种工况条件，及每种工况条件对应的k种工况因素；
[0026] 采集的第b种工况条件对应的工况信息GKb＝[Z1b Z2b ... Zib ... Zkb]，其中Zib表示数控机床第b种工况条件下的第i个工况因素的工况协变量。
[0027] 技术方案步骤2中所述考虑工况差异的数控机床变参数幂律分布可靠性建模，具体包括以下步骤：
[0028] 2.1、建立两参数幂律分布模型：
[0029] 假定数控机床的故障间隔工作时间t服从两参数的幂律分布，其故障率函数λ(t)为：
[0030]
[0031] 其中：β为幂律分布的形状参数；
[0032] η为幂律分布的尺度参数；
[0033] 即而得到数控机床的可靠度函数R(t)为：
[0034]
[0035] 概率密度函数f(t)为：
[0036]
[0037] 累积分布函数F(t)为：
[0038]
[0039] 平均故障间隔工作时间MTBF为：
[0040]
[0041] 2.2、建立变参数幂律分布模型：
[0042] 考虑m种工况条件下k种工况因素对数控机床可靠性水平的影响，对幂律分布中的尺度参数η建立变参数幂律分布模型；
[0043] 选择第1种工况条件作为基准工况条件，则数控机床在第b种工况条件下幂律分布的尺度参数ηb为：
[0044]
[0045] 其中：Zi1为第1种工况条件下第i个工况因素的工况协变量，即基准工况条件下第i个工况因素的工况协变量；
[0046] Zib为第b种工况条件下第i个工况因素的工况协变量；
[0047] η1为在基准工况条件下幂律分布的尺度参数；
[0048] ξi为第i种工况因素的系数；
[0049] 将式(6)代入式(1)得到第b种工况条件下故障率函数λb(t)为：
[0050]
[0051] 将式(7)代入式(2)中，得到数控机床在第b种工况条件下考虑k种工况因素的可靠度函数Rb(t)为：
[0052]
[0053] 将式(8)代入式(3)中，得到数控机床在第b种工况条件下考虑k种工况因素的概率密度函数fb(t)为：
[0054]
[0055] 将式(7)代入式(4)中，得到数控机床在第b种工况条件下考虑k种工况因素的累积分布函数Fb(t)为：
[0056]
[0057] 将式(9)代入式(5)中，得到数控机床在第b种工况条件下考虑k种工况因素的数控机床的平均故障间隔工作时间MTBFb为：
[0058]
[0059] 其中：Γ(1+1/(1+β))——Gamma函数。
[0060] 技术方案步骤3中所述筛选工况因素，具体包括以下步骤：
[0061] 3.1、建立比例模型：
[0062] 所述比例模型为第b种工况条件与基准工况条件故障率函数的比例关系；
[0063] 将基准工况条件下幂律分布的尺度参数η1代入式(1)得到基准工况条件下故障率函数λ1(t)，即而得到式(7)中第b种工况条件下故障率函数与基准工况条件下故障率函数的比例关系为：
[0064]
[0065] 其中：Zi1为第1种工况条件下第i个工况因素的工况协变量；
[0066] Zib为第b种工况条件下第i个工况因素的工况协变量；
[0067] θi为第i种工况因素的比例系数，且θi＝ξi(1+β)；
[0068] 3.2、采用偏似然函数估计比例系数：
[0069] 为了估计式(12)中的第i种工况因素的比例系数θi，使用偏似然函数L(θi)：
[0070]
[0071] 其中：Zibj为第j个故障间隔工作时间所对应第b种工况条件的第i个工况因素的工况协变量；
[0072] 将式(13)化简取对数后得：
[0073]
[0074] 将式(14)对比例系数θi求偏导数并令其等于0得：
[0075]
[0076] 将采集的数据整理后代入式(15)，采用Newton-raphson(牛顿迭代法)数值算法进行求解，即可得出第i种工况因素的比例系数θi；
[0077] 3.3、采用卡方检验法筛选工况因素：
[0078] 第i种工况因素的卡方检验统计量χ2(k)i为：
[0079]
[0080] 其中：θi为第i种工况因素的比例系数；
[0081] 为在所有工况条件下第i种工况因素的最小工况协变量；
[0082] Zibj为第j个故障间隔工作时间所对应第b种工况条件下的第i个工况因素的工况协变量。
[0083] 技术方案步骤4中所述模型参数估计，具体包括以下步骤：
[0084] 在选择了合适的工况因素的基础上，将采集的数据整理后代入式(8)和式(9)采用极大似然估计方法估计在基准工况条件下的故障率函数参数尺度参数η1和形状参数β；极大似然估计函数L(η1，β)为：
[0085]
[0086] 其中：tq，j为第q台数控机床的第j个故障间隔工作时间；
[0087] 为第a台数控机床的截尾时间；
[0088] 将式(17)化简取对数后得：
[0089]
[0090] 对式(18)的形状参数β和尺度参数η1分别求偏导数并令其等于零得：
[0091]
[0092]
[0093] 根据式(15)求得的第i种工况因素的比例系数θi，将ξi＝θi/(1+β)代入式(19)和式 (20)，采用Newton-raphson(牛顿迭代法)数值算法得到基准工况条件下的故障率函数的尺度参数η1与形状参数β的估计值，同时得到第i种工况因素的系数ξi的值，再根据式(6) 得到第b种工况条件下幂律分布的尺度参数ηb的估计值。
[0094] 技术方案步骤5中所述考虑工况差异的数控机床混合变参数幂律分布可靠性建模，具体包括以下步骤：
[0095] 5.1、对工况条件进行分组：
[0096] 为了建立考虑m种工况条件差异的数控机床混合变参数幂律分布可靠性模型，需依据斯特奇斯经验公式将m种工况条件划分成d组：
[0097] d＝[1+3.32lg(m)] (21)
[0098] 其中：[]代表取整。
[0099] 令δb为第b种工况条件下的工况分组系数，则工况分组系数δb为：
[0100]
[0101] 其中：ξi为第i种工况因素的系数；
[0102] 设δmax为工况分组系数中的最大值，δmin为工况分组系数中的最小值，则工况分组组距Δδ为：
[0103]
[0104] 设第e组内有me种工况条件，且 根据工况分组组距Δδ的大小，得到第e组工况分组系数的区间为(δmax-(e-1)Δδ，δmax-eΔδ)，再根据第e组工况分组系数的区间得到第e组内的me种工况条件，即得到第e种工况水平下的me种工况条件，利用混合变参数幂律分布对同一工况水平下工况条件的数控机床可靠性函数进行拟合，同时，能够综合的考虑多种工况因素对可靠性水平评估的影响。
[0105] 5.2建立混合变参数幂律分布可靠性模型
[0106] 本发明中混合变参数幂律分布模型适用于数控机床存在一种或多种不同工况条件的情况。根据工况条件分组结果，对第e种工况水平下的me种工况条件建立混合变参数幂律分布模型。令第b种工况条件下的混合权系数为pb，同时 且pb＞0，b∈me。则第e种工况水平对应的混合累积分布函数Fe(t)为：
[0107]
[0108] 第e种工况水平对应的混合概率密度分布函数fe(t)为：
[0109]
[0110] 第e种工况水平对应的混合可靠度函数Re(t)为：
[0111]
[0112] 第e种工况水平对应的混合平均故障间隔工作时间MTBFe为：
[0113]
[0114] 将式(10)、(9)、(8)和(11)代入式(24)、(25)、(26)和(27)得到考虑工况差异的数控机床混合变参数幂律分布累积分布函数Fe(t)为：
[0115]
[0116] 概率密度分布函数fe(t)为：
[0117]
[0118] 可靠度函数Re(t)为：
[0119]
[0120] 平均故障间隔工作时间MTBFe为：
[0121]
[0122] 5.3估计混合权系数pb
[0123] 设第e组内工况协变量Zib的标准差σib(Z)为：
[0124]
[0125] 其中：Zib表示数控机床第b种工况条件下的第i个工况因素的工况协变量；
[0126] 为第e种工况水平下的所有工况条件中第i种工况因素的工况协变量的平均值。
[0127] 则混合变参数幂律分布的相对权系数Pb为：
[0128]
[0129] 将相对权系数进行归一化处理，得到混合权系数pb为：
[0130]
[0131] 技术方案步骤6中所述假设检验是指：
[0132] 采用Kolmogorov-Sminov(柯尔莫诺夫-斯米尔诺夫)假设检验方法，即K-S假设检验，对混合变参数幂律分布进行假设检验；
[0133] 设第e个混合变参数幂律分布的K-S检验统计量De为：
[0134] De＝max|Fn(t)-Fe(t)|，e＝1，2，…，d (35)
[0135] 其中：Fn(t)——样本经验分布累计分布函数；
[0136] 其中，样本经验累计分布函数Fn(t)为：
[0137]
[0138] 在给定K-S检验显著性水平α下的临界值Dn，α，当由试验数据计算得到De大于临界值 Dn，α时，可判断原假设不成立。当由试验数据计算得到De小于临界值Dn，α时，则可判断原假设成立。
[0139] 与现有技术相比本发明的有益效果是：
[0140] 1、数控机床在实际切削加工过程中，由于处在不同的行业和加工不同类型的零件，导致数控机床所受的切削载荷相差很大，然而传统的数控机床可靠性建模与评估方法，由于未考虑工况条件的差异，导致同批次的数控机床可靠性评估的结果可能相差很大，不能够准确的评估处于不同工况条件下的数控机床的真实可靠性水平。本发明提出了一种考虑工况差异的数控机床混合变参数幂律分布可靠性建模与评估方法，即幂律分布中的尺度参数随工况条件的变化而变化，相比于传统的可靠性建模与评估方法，此方法能够准确的对不同工况条件下的数控机床的可靠性水平进行准确评估。
[0141] 2、现有的数控机床可靠性建模与评估方法是依据数控机床的故障数据拟合成运用威布尔模型建立数控机床可靠性模型，然而使用威布尔分布模型建立可靠性模型的前提条件是故障间隔工作时间独立同分布，也就是说数控机床每次发生故障之后都需要“修复如初”。然而，数控机床是集机、电、液、气于一体的复杂可修系统，在数控机床发生故障之后，往往很难做到修复如初，所以，本发明专利采用混合变参数幂律分布建立数控机床可靠性模型，与现有的方法相比更加符合实际情况。
[0142] 3、在幂律分布中，不同工况条件对应的尺度参数是不同的，同时，为了能够准确评价工况条件差异对数控机床可靠性水平的影响，本发明专利采用混合分布，将同一工况水平下的多种工况条件对应的幂律分布整合成一个分布函数，混合分布中混合权系数的大小表示在工况水平下的工况条件对数控机床可靠性水平的影响程度。从而能够准确的评价不同工况水平下的数控机床的可靠性水平。
[0143] 4、数控机床在实际切削加工过程中，影响数控机床可靠性的工况因素种类很多，而对数控机床可靠性水平影响程度却相差较大，为了能够找出对数控机床可靠性水平影响较大的工况因素，本发明专利采用卡方检验法进行工况因素筛选。
[0144] 5、为了能够准确的评价多种工况条件下的数控机床的可靠性水平，本发明专利制定了《数控机床故障记录表》、《数控机床工况信息记录表》和《数控机床运行记录表》，通过对数控机床进行数据采集，获取数控机床的故障数据、工况信息及运行数据。以便能够准确真实的反应不同工况条件下的数控机床可靠性水平。
[0145] 6、本发明提出的一种考虑工况差异的可靠性建模与评估方法，不仅局限于数控机床的可靠性建模与评估的研究，也可用于其他机电液复杂产品的可靠性建模与评估。
附图说明
[0146] 下面结合附图对本发明专利作进一步的说明：
[0147] 图1是本发明所述考虑工况差异的数控机床可靠性评估方法流程图；
[0148] 图2是卡方检验法的流程图；
[0149] 图3是数控机床的可靠度函数曲线示意图；
[0150] 图4是数控机床的概率密度函数曲线示意图。
具体实施方式
[0151] 参阅图1，本发明专利所述的数控机床可靠性建模与评估方法，通过数控机床用户现场可靠性试验，获取大量的数控机床故障数据、工况信息和运行数据的基础上，建立了考虑工况差异的数控机床混合变参数幂律分布可靠性模型。首先，建立了幂律分布尺度参数随工况条件的变化而变化的变参数幂律分布模型，其次，建立比例模型，并利用偏似然函数估计比例模型第i种工况因素比例系数θi，采用卡方检验法筛选出对数控机床可靠性水平影响较大的工况因素，利用极大似然法估计基准工况条件下幂律分布尺度参数η1和形状参数β，进而得到第i种工况因素的系数ξi和第b种工况条件下幂律分布尺度参数ηb；利用工况分组系数δ对工况条件进行分组，再根据分组的结果，对第b种工况条件下混合权系数pb进行估计，得到考虑工况差异的混合变参数幂律分布可靠性模型，最后，对混合变参数幂律分布进行 Kolmogorov-Sminov(柯尔莫诺夫-斯米尔诺夫)假设检验。这种方法较好的解决了考虑工况差异的数控机床可靠性建模与评估的问题。
[0152] 1.数据采集与整理
[0153] 1.1采集故障信息
[0154] 针对数控机床开展用户现场可靠性试验，在可靠性试验过程中详细记录数控机床的故障信息，包括故障数据、工况信息及运行数据等，为后期的可靠性评估等提供基础数据。
[0155] 故障数据记录时应明确每一个故障的故障现象、故障部位、故障原因及故障处理等内容，并详细填写《数控机床故障记录表》，详见表1。
[0156] 工况信息是指在数控机床在发生故障或截尾时所处的工况条件中的各个工况因素，包括：数控机床信息(数控机床型号、数控机床编号等)、工件信息(工件材料、工件硬度等)、刀具信息(刀具名称、刀具材料等)、切削参数(背吃刀量、进给速度、切削速度等)及环境因素(环境温度、湿度等)等。并详细填写《数控机床工况信息记录表》，详见表2。
[0157] 运行数据主要记录数控机床在实际运行过程中的运行时间、停机时间和维修时间等数据。并详细填写《数控机床运行记录表》，详见表3。
[0158] 采集的工况信息、故障数据及运行数据，为后期可靠性建模与评估方法等提供基础数据。
[0159] 表1数控机床故障记录表
[0160]
[0161] 表2数控机床工况信息记录表
[0162]
[0163]
[0164] 表3数控机床运行记录表
[0165]
[0166]
[0167] 1.2数据整理
[0168] 假设共有K台数控机床，在数控机床用户现场可靠性试验过程中总共采集了n条故障数据，其中包括r(r＝n1+n2+...+nq+...+nK)条完整的故障数据及n-r(n-r＝n1+n2+...+na+...+nK) 条截尾数据(每台数控机床最多存在一个截尾数据)，nq表示第q台数控机床共发生的故障数量，na＝1表示第a台数控机床存在截尾数据，na＝0表示第a台数控机床不存在截尾数据，ta为第a台数控机床的可靠性试验结束时间，且第j个故障发生时间为tj，则第q台数控机床的第j个故障间隔工作时间tq，j＝tj-tj-1， 为第q台数控机床在可靠性试验周期内的最后一个故障发生的时间，当 不等于ta时，则第a台数控机床存在截尾数据，即na＝1，且第a台数控机床的截尾时间为 采集的第j条故障数据对应的工况信息 GKq，j＝[Z1bj Z1bj … Zibj … Zkbj]，采集的第a台数控机床的截尾数据对应的工况信息 GKa＝[Z1ba Z1ba … Ziba … Zkba]，其中Zibj表示第j个故障间隔工作时间所对应第b种工况条件下的第i个工况因素的工况协变量，Ziba表示第a台数控机床的截尾时间所对应第b 种工况条件下的第i个工况因素的工况协变量。
[0169] 将现场跟踪试验采集的工况信息进行整理，最终确定了数控机床的m种工况条件，及每种工况条件对应的k种工况因素。采集的第b种工况条件对应的工况信息 GKb＝[Z1b Z2b … Zib … Zkb]，其中Zib表示数控机床第b种工况条件下的第i个工况因素的工况协变量。
[0170] 2.考虑工况差异的数控机床变参数幂律分布可靠性建模
[0171] 2.1建立两参数幂律分布模型
[0172] 本发明假定数控机床的故障间隔工作时间t服从两参数的幂律分布，其故障率函数λ(t) 为：
[0173]
[0174] 其中：β为幂律分布的形状参数；
[0175] η为幂律分布的尺度参数；
[0176] 即而得到数控机床的可靠度函数R(t)为：
[0177]
[0178] 概率密度函数f(t)为：
[0179]
[0180] 累积分布函数F(t)为：
[0181]
[0182] 平均故障间隔工作时间MTBF为：
[0183]
[0184] 2.2建立变参数幂律分布模型
[0185] 考虑m种工况条件下k种工况因素对数控机床可靠性水平的影响，对幂律分布中的尺度参数η建立了变参数幂律分布模型。选择第1种工况条件作为基准工况条件，则数控机床在第b种工况条件下幂律分布的尺度参数ηb为：
[0186]
[0187] 其中：Zi1为第1种工况条件下第i个工况因素的工况协变量，即在基准工况条件下第i 个工况因素的工况协变量；
[0188] Zib为第b种工况条件下第i个工况因素的工况协变量；
[0189] η1为数控机床在基准工况条件下幂律分布的尺度参数；
[0190] ξi为第i种工况因素的系数；
[0191] 将式(6)代入式(1)可以得到第b种工况条件下故障率函数λb(t)为：
[0192]
[0193] 将式(7)代入式(2)中，得到数控机床在第b种工况条件下考虑k种工况因素的可靠度函数Rb(t)为：
[0194]
[0195] 将式(8)代入式(3)中，得到数控机床在第b种工况条件下考虑k种工况因素的概率密度函数fb(t)为：
[0196]
[0197] 将式(7)代入式(4)中，得到数控机床在第b种工况条件下考虑k种工况因素的累积分布函数Fb(t)为：
[0198]
[0199] 将式(9)代入式(5)中，得到数控机床在第b种工况条件下考虑k种工况因素的数控机床的平均故障间隔工作时间MTBFb为：
[0200]
[0201] 其中：Γ(1+1/(β+1))——Gamma函数；
[0202] 3.筛选工况因素
[0203] 3.1建立比例模型
[0204] 比例模型为第b种工况条件与基准工况条件故障率函数的比例关系。
[0205] 将基准工况条件下幂律分布的尺度参数η1代入式(1)得到基准工况条件下故障率函数λ1(t)，即而得到第b种工况条件下故障率函数(式(7))与基准工况条件下故障率函数的比例关系为：
[0206]
[0207] 其中：Zi1为第1种工况条件下第i个工况因素的工况协变量；
[0208] Zib为第b种工况条件下第i个工况因素的工况协变量；
[0209] θi为第i种工况因素的比例系数，且θi＝ξi(1+β)。
[0210] 3.2采用偏似然函数估计比例系数
[0211] 为了估计式(12)中的第i种工况因素的比例系数θi，使用偏似然函数L(θi)：
[0212]
[0213] 其中：Zibj为第j个故障间隔工作时间所对应第b种工况条件的第i个工况因素的工况协变量；
[0214] 将式(13)化简取对数后得：
[0215]
[0216] 将式(14)对比例系数θi求偏导数并令其等于0得：
[0217]
[0218] 将采集的数据整理后代入式(15)，采用Newton-raphson数值算法进行求解，即可得出第i种工况因素的比例系数θi。
[0219] 3.3采用卡方检验法筛选工况因素
[0220] 由于数控机床在实际切削加工过程中，处在复杂的工况条件下，为了去除对数控机床可靠性影响较小或无影响的工况因素，采用卡方检验法对工况因素进行筛选，工况因素筛选流程见附图2。第i种工况因素的卡方检验统计量χ2(k)i为：
[0221]
[0222] 其中：θi为第i种工况因素的比例系数；
[0223] 为在所有工况条件下第i种工况因素的最小工况协变量；
[0224] Zibj为第j个故障间隔工作时间所对应第b种工况条件下第i个工况因素的工况协变量；
[0225] 对于给定显著性水平α下，参照表4，查询 的值，根据附图2进行筛选，判断当 i＝1，2，…，k时χ2(k)i是否全部大于 当且仅当k个χ2(k)i全部大于 时，k种工况因素对数控机床可靠性水平影响较大。
[0226] 表4卡方分布临界值表
[0227]
[0228]
[0229]
[0230] 4.模型参数估计
[0231] 在选择了合适的工况因素的基础上，将采集的数据整理后代入式(8)和式(9)采用极大似然估计方法估计在基准工况条件下的故障率函数参数尺度参数η1和形状参数β。则极大似然估计函数L(η1，β)为：
[0232]
[0233] 其中：tq，j为第q台数控机床的第j个故障间隔工作时间；
[0234] 为第a台数控机床的截尾时间；
[0235] 将式(17)化简取对数后得：
[0236]
[0237] 对式(18)的参数β和η1分别求偏导数并令其等于零得：
[0238]
[0239]
[0240] 根据式(15)求得的第i种工况因素的比例系数θi，将ξi＝θi/(1+β)代入式(19)和式 (20)，采用Newton-raphson数值算法得到基准工况条件下的故障率函数的尺度参数η1与形状参数β的估计值，同时得到第i种工况因素的系数ξi的值，再根据式(6)得到第b种工况条件下幂律分布的尺度参数ηb的估计值。
[0241] 5.考虑工况差异的数控机床混合变参数幂律分布可靠性建模
[0242] 5.1对工况条件进行分组
[0243] 为了建立考虑m种工况条件差异的数控机床混合变参数幂律分布可靠性模型，需依据斯特奇斯经验公式将m种工况条件划分成d组：
[0244] d＝[1+3.32lg(m)] (21)
[0245] 其中：[]代表取整。
[0246] 令δb为第b种工况条件下的工况分组系数，则工况分组系数δb为：
[0247]
[0248] 其中：ξi为第i种工况因素的系数；
[0249] 设δmax为工况分组系数的最大值，δmin为工况分组系数的最小值，则工况分组组距Δδ为：
[0250]
[0251] 设第e组内有me种工况条件，且 根据工况分组组距Δδ的大小，得到第e组工况分组系数的区间为(δmax-(e-1)Δδ，δmax-eΔδ)，再根据第e组工况分组系数的区间得到第e组内的me种工况条件，即得到第e种工况水平下的me种工况条件，利用混合变参数幂律分布对同一工况水平下工况条件的数控机床可靠性函数进行拟合，同时，能够综合的考虑多种工况因素对可靠性水平评估的影响。
[0252] 5.2建立混合变参数幂律分布可靠性模型
[0253] 本发明中混合变参数幂律分布模型适用于数控机床存在一种或多种工况条件差异的情况。根据分组结果，对第e种工况水平下的me种工况条件建立混合变参数幂律分布模型。令第b种工况条件下的混合权系数为pb，同时 且pb＞0，b∈me。则第e种工况水平对应的混合累积分布函数Fe(t)为：
[0254]
[0255] 第e种工况水平对应的混合概率密度分布函数fe(t)为：
[0256]
[0257] 第e种工况水平对应的混合可靠度函数Re(t)为：
[0258]
[0259] 第e种工况水平对应的混合平均故障间隔工作时间MTBFe为：
[0260]
[0261] 将式(10)、(9)、(8)和(11)代入式(24)、(25)、(26)和(27)得到考虑工况差异的数控机床混合变参数幂律分布累积分布函数Fe(t)为：
[0262]
[0263] 概率密度分布函数fe(t)为：
[0264]
[0265] 可靠度函数Re(t)为：
[0266]
[0267] 平均故障间隔工作时间MTBFe为：
[0268]
[0269] 5.3估计混合权系数pb
[0270] 设第e组内工况协变量Zib的标准差σib(Z)为：
[0271]
[0272] 其中： 为第e种工况水平下的所有工况条件中第i种工况因素的工况协变量的平均值。
[0273] 则混合变参数幂律分布的相对权系数Pb为：
[0274]
[0275] 将相对权系数进行归一化处理，得到混合权系数pb为：
[0276]
[0277] 6.假设检验
[0278] 本发明采用Kolmogorov-Sminov假设检验方法，即K-S假设检验，对混合变参数幂律分布进行假设检验。
[0279] 设第e个混合变参数幂律分布的K-S检验统计量De为：
[0280] De＝max|Fn(t)-Fe(t)|，e＝1，2，…，d (35)
[0281] 其中：Fn(t)——样本经验分布累计分布函数；
[0282] 其中，样本经验累计分布函数Fn(t)为：
[0283]
[0284] 给定K-S检验显著性水平α下，查询表5得到K-S检验的临界值Dn，α。当由试验数据计算得到De＞Dn，α时，可判断原假设不成立。当De＜Dn，α时，则可判断原假设成立。
[0285] 表5 K-S检验临界值表
[0286]
[0287]
[0288] 下面结合实例对数控机床进行可靠性建模与评估，对本发明中的具体施加方式进一步说明。进行了6个多月的数控机床用户现场跟踪试验，收集到数控机床在20种工况条件下的 151条故障数据。具体建模与评估过程如下：
[0289] (1)故障数据、工况信息及运行数据的获取
[0290] 对受试数控机床进行现场跟踪试验，记录数控机床的故障数据、工况信息及运行数据并填写表1、表2及表3。
[0291] (2)数据整理
[0292] 根据该批数控机床用户现场可靠性试验的具体情况，初步考虑了切削速度、扭矩、切削力、温度及换刀次数等五种工况因素对数控机床进行可靠性建模，第b种工况条件下第1种工况因素为切削速度，第b种工况条件下第2种工况因素为扭矩，第b种工况条件下第3种工况因素为切削力，第b种工况条件下第4种工况因素为温度，第b种工况条件下第5种工况因素为换刀次数，通过计算整理得到数控机床在不同工况条件下的故障信息，如表6所示。
[0293] 表6数控机床在不同工况条件下的故障信息
[0294]
[0295]
[0296]
[0297] (3)筛选工况因素
[0298] 根据表6中的数据利用式(15)得到5种工况因素的比例系数θi，分别是：θ1＝0.0050，θ2＝0.0350，θ3＝0.0528，θ4＝0.0319，θ5＝0.0251。
[0299] 为了选择对数控机床可靠性水平影响较大的工况因素，以工况条件1作为基准工况条件，初步选择5种工况因素(切削速度、扭矩、切削力、温度及换刀次数)，即k＝5，通过式(16) 进行计算出5种工况因素的卡方检验统计量，计算结果如表7所示，发现χ2(5)4与χ2(5)5均小于 且换刀次数的卡方检验统计量χ2(5)5最小，故忽略换刀次数对数控机床可靠性水平的影响，再对剩余其他4种工况因素重新进行卡方检验，此时k＝4，再计算出4种2 2
工况因素的卡方检验统计量，发现χ(4)4小于 且温度的卡方检验统计量χ(4)4最小，故忽略温度对数控机床可靠性水平的影响，再对剩余其他3种工况因素重新进行卡方检验，此时k＝3，再计算出3种工况因素的卡方检验统计量，计算结果如表2所示，发现χ2(3)1、χ2(3)2与χ2(3)3均大于 可以看出工况因素中扭矩、切削力与切削速度对数控机床可靠性水平影响较大，温度与换刀次数对数控机床可靠性水平影响较小，故对数控机床进行可靠性评估时，需要考虑扭矩、切削力与切削速度三种工况因素对数控机床可靠性水平的影响。
[0300] 表7卡方检验法
[0301]
[0302] (4)模型参数估计
[0303] 在选择对数控机床可靠性水平影响较大的工况因素后，根据式(19)和(20)，利用极大似然法估计出在基准工况条件在幂律分布的尺度参数为η1＝642.4，形状参数为β＝0.2865。即而求得工况因素切削速度的系数ξ1＝0.0038，工况因素扭矩的系数ξ2＝0.0262，工况因素切削力的系数ξ3＝0.0041；根据式(16)得到第b种工况条件下尺度参数ηb，详见表
8。
[0304] (5)考虑工况差异的数控机床混合变参数幂律分布的可靠性建模
[0305] 根据式(21-23)对工况条件进行分组，利用式(32-34)估计第b种工况条件的混合权系数pb，详见表8；利用式(30-31)得到第e种工况水平下的混合变参数幂律分布可靠度模型Re(t) 及平均故障间隔工作时间MTBFe，详见表8。
[0306] (6)假设检验
[0307] 假设该数控机床的故障间隔工作时间服从混合变参数幂律分布。
[0308] 利用式(35-36)，计算得到第e个K-S检验统计量De为：
[0309] D1＝0.1143，D2＝0.0884，D3＝0.0458，D4＝0.0656，D5＝0.1192。
[0310] 取K-S检验显著性水平α＝0.01，当n＝151时，查表5得：Dn，α＝0.1326。
[0311] 由于De均小于Dn，α，因此，该数控机床的故障间隔工作时间服从混合变参数幂律分布。
[0312] 表8混合变参数幂律分布模型求解
[0313]
[0314]
