CN102981452A - 数控机床三类功能部件的可靠性建模与可靠性评估方法 - Google Patents

Abstract

本发明涉及一种数控机床三类功能部件的可靠性建模与可靠性评估方法，该方法包括下述步骤：根据数控机床可靠性现场试验数据绘制故障间隔时间的WPP图；如果WPP图趋向于一条直线，则对功能部件寿命进行二参数威布尔分布建模并对其可靠性进行评估；如果WPP图具有明显的拐点，则采用二重威布尔混合分布建模并对其参数进行估计；对二重威布尔混合模型进行可靠性评估。本发明基于运行状态的数控机床三类功能部件的可靠性现场试验数据进行功能部件的可靠性建模和评估，采用多重威布尔分布模型揭示这些功能部件产品寿命周期过程中故障分布规律，并进行综合参数估计，较好地解决了功能部件这类复杂系统的可靠性建模和评估问题。

Description

数控机床三类功能部件的可靠性建模与可靠性评估方法

技术领域

本发明属于数控机床关键功能部件的可靠性工程领域，涉及一种高速主轴单元、刀库及机械手、动力伺服刀架等三种关键功能部件的可靠性建模与可靠性评估的方法。

背景技术

数控机床功能部件的可靠性对整机的可靠性有重要影响，特别是高速主轴单元、刀库及机械手、动力伺服刀架对主机可靠性的影响尤为明显，是国产数控机床可靠性的薄弱环节，是影响整机可靠性水平的关键部件。科技查新表明，以往虽然已有关于数控车床功能部件的可靠性方面的报道，但这些报道主要集中在功能部件的故障现象分析，以及功能部件的故障对数控机床主机可靠性的影响，大多数都是研究数控机床整机可靠性时，指出某一功能部件的故障问题，在研究功能部件故障分布规律时，建模常用的方法是将功能部件的失效数据拟合成简单的二参数或三参数威布尔分布，在确定出分布参数之后,对失效数据进行可靠性评估和预测。由于数控机床功能部件是由多个零部件构成的复杂系统，零部件之间具有故障相关性，而且其故障可能是在多种失效机理共同作用下发生的,并且在不同的寿命阶段,不同的失效机理对系统的失效起到不同的主导作用。例如，由于涉及原材料、制造工艺等方面的差异，在不同时期生产的功能部件产品可能有着不同的寿命分布。随着数控装备功能的日趋先进以及一些尚未成熟的新技术的应用，日趋复杂的数控产品对可靠性工程提出了更高要求，如果单纯使用标准威布尔分布来描述,会出现较大的误差。

发明内容

本发明要解决的技术问题是提供一种基于运行状态的数控机床动力伺服刀架、刀库及机械手、高速主轴单元等功能部件的可靠性现场试验数据进行功能部件的可靠性建模和评估的数控机床三类功能部件的可靠性建模与可靠性评估方法，该方法采用多重威布尔分布模型揭示这些功能部件产品寿命周期过程中故障分布规律，并进行综合参数估计，较好地解决了功能部件这类复杂系统的可靠性建模和评估问题。

为了解决上述技术问题，本发明的数控机床三类功能部件的可靠性建模与可靠性评估方法包括下述步骤：

步骤一：

根据数控机床可靠性现场试验数据绘制故障间隔时间的WPP图(Weibull ProbabilityPaper)，WPP图横坐标x_i及纵坐标y_i如下：

$\left\{\begin{array}{c}{x}_i=\ln \left(t_i\right)\\ y_i=\ln [-\ln (1-\frac{i-0.3}{n+0.4})]\end{array}\right.---\left(1\right)$

式（1）中，t_i(i=1,2,…,n)为关键功能部件第i次故障相对于考核起始时刻t₀的时刻；n是关键功能部件发生故障的总次数；

步骤二：

判断故障间隔时间WPP图是否趋向于一条直线；如果WPP图趋向于一条直线，则转步骤三；如果WPP图具有明显的拐点，则转步骤四；

步骤三：

对功能部件寿命进行二参数威布尔分布建模并对其可靠性进行评估；

步骤四：

采用二重威布尔混合分布建模并对其参数进行估计；所述二重威布尔混合模型的可靠度函数R(t)如下：

$R\left(t\right)=p{R}_1\left(t\right)+q{R}_2\left(t\right)$

$=\mathrm{pexp}[-{\left(\frac{t}{{\mathrm{\α}}_1}\right)}^{{\mathrm{\β}}_1}]+\mathrm{qexp}[-{\left(\frac{t}{{\mathrm{\α}}_2}\right)}^{{\mathrm{\β}}_2}]---\left(10\right)$

式(10)中：

t——时间变量；

α₁和β₁——第一个子分布的形状参数和尺度参数；

α₂和β₂——第二个子分布的形状参数和尺度参数；

p和q——分别为第一个和第二个子分布在二重威布尔混合模型中的比重，并且0≤p≤1，0≤q≤1，p+q＝1；

步骤五：

对二重威布尔混合模型进行可靠性评估，关键功能部件的平均无故障时间MTBF(MeanTime Between Failures)观测值按式（15）进行计算：

$\mathrm{MTBF}=\frac{1}{N_0}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{t}_i=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{t}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{r}_i}---\left(15\right)$

式（15）中：

N₀—在评定周期内关键功能部件累计故障频数；

n—关键功能部件抽样台数；

t_i—在评定周期内第i台关键功能部件的实际工作时间，单位是小时；

r_i—在评定周期内第i台关键功能部件出现的故障频数；

MTBF的点估计值按式(24)计算：

$\mathrm{MTBF}=p\·{\mathrm{\α}}_1\·\mathrm{\Γ}(1+\frac{1}{{\mathrm{\β}}_1})+q\·{\mathrm{\α}}_2\·\mathrm{\Γ}(1+\frac{1}{{\mathrm{\β}}_2})---\left(24\right)$

式(24)中：

α₁和β₁——第一个子分布的形状参数和尺度参数；

α₂和β₂——第二个子分布的形状参数和尺度参数；

p和q——分别为第一个和第二个子分布在二重威布尔混合模型中的比重，并且0≤p≤1，0≤q≤1，p+q=1；

采用定时截尾的试验指标区间估计，获得MTBF的置信水平为1-α的双侧置信区间为：

其中α为显著性水平，r为关键功能部件发生故障的次数，T^*为定时截尾试验总试验时间。

所述步骤三中二参数威布尔分布建模方法如下：

威布尔分布的概率密度函数为：

$f\left(t\right)=\frac{\mathrm{\β}}{\mathrm{\α}}{\left(\frac{t}{\mathrm{\α}}\right)}^{\mathrm{\β}-1}\exp [-{\left(\frac{t}{\mathrm{\α}}\right)}^{\mathrm{\β}}],t\≥0---\left(4\right)$

累积故障分布函数为：

$F\left(t\right)=1-\exp [-{\left(\frac{t}{\mathrm{\α}}\right)}^{\mathrm{\β}}],t\≥0---\left(5\right)$

式（4）、式（5）中，t为时间变量，β为形状参数，β>0；α为尺度参数，α>0；

利用式（7）通过迭代法求出α和β的估计值

和

将和

值代入式（4）、式（5）得到二参数威布尔分布的概率密度函数f(t)和累积故障分布函数F（t）；

$\left\{\begin{array}{c}\frac{1}{\mathrm{\β}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{t}_i^{\mathrm{\β}}\ln t_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{t}_i^{\mathrm{\β}}}-\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\ln t_i\\ {\mathrm{\α}}^{\mathrm{\β}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{t}_i^{\mathrm{\β}}\end{array}\right.---\left(7\right)$

（7）中，t_i(i=1,2,…,n)为关键功能部件第i次故障相对于考核起始时刻t₀的时刻；n是关键功能部件发生故障的总次数；

所述步骤三中在二参数威布尔分布建模后还需对二参数威布尔分布模型进行假设检验。

二参数威布尔分布模型假设检验采用d检验法，具体如下：

将D_n与临界值D_n,α进行比较，满足式（8），则接受原假设，否则拒绝原假设；

$D_n=\underset{-\&infin;<t<+\&infin;}{\sup }|F_n\left(t_i\right)-F_0\left(t_i\right)|=\max \left\{d_i\right\}<D_{n,\mathrm{\&alpha;}}---\left(8\right)$

式(8)中：

F₀(t_i)——原假设分布函数，即式（5）中的累积故障分布函数F(t)；

F_n(t)——样本大小为n的经验分布函数，采用中位秩法，有

$F_n\left(t\right)=\left\{\begin{array}{cc}0,& t<t_i\\ \frac{i-0.3}{n+0.4},& t_i\&le;t<t_{i+1}\\ 1,& t\&GreaterEqual;t_n\end{array}\right.---\left(9\right)$

D_n,α——临界值；

式(8)、式(9)中t_i是将n个实验数据由小到大的排列序列，i=1,2,3...(n-1),n为样本大小，即关键功能部件发生故障的总次数。

所述步骤三中，根据似然比（Likelihood Ratio）理论对二参数威布尔分布模型的可靠性进行评估，具体如下：

二参数威布尔分布的对数似然函数为：

$\ln L(\mathrm{\&alpha;},\mathrm{\&beta;})=n\ln \mathrm{\&beta;}-\mathrm{n\&beta;}\ln \mathrm{\&alpha;}+(\mathrm{\&beta;}-1)\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\ln t_i-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\frac{t_i}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}}---\left(17\right)$

式（17）中，t_i(i=1,2,…,n)为关键功能部件第i次故障相对于考核起始时刻t₀的时刻；n是关键功能部件发生故障的总次数；

令原假设H₀:β=β₀，备择假设H:β≠β₀，则似然比为：

$l=\frac{L(\widetilde{\mathrm{\&alpha;}},{\mathrm{\&beta;}}_0)}{L(\widehat{\mathrm{\&alpha;}},\widehat{\mathrm{\&beta;}})}---\left(18\right)$

式（18）中：

和

是α和β的极大似然估计值；

是在假设H₀:β=β₀成立，参数α的条件极大似然估计，它由对lnL(α,β₀)关于α施行极大化来求得：

$\widetilde{\mathrm{\&alpha;}}={\left(\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i^{\mathrm{\&beta;}}\right)}^{\frac{1}{\mathrm{\&beta;}}}---\left(19\right)$

同时

服从自由度为1的χ²分布；设置信水平为γ，令X=-2lnl；解不等式

得到β₀值的集合，即在置信水平为γ时，形状参数β的置信区间；

通过公式(20)得到从形状参数β角度考虑的MTBF的置信区间；

$\mathrm{MTBF}=\widetilde{\mathrm{\&alpha;}}\mathrm{\&Gamma;}(1+\frac{1}{{\mathrm{\&beta;}}_0})---\left(20\right)$

令原假设H₀:α=α₀，备择假设H:α≠α₀，则似然比为：

$l=\frac{L({\mathrm{\&alpha;}}_0,\widetilde{\mathrm{\&beta;}})}{L(\widehat{\mathrm{\&alpha;}},\widehat{\mathrm{\&beta;}})}---\left(21\right)$

式(21)中：

和

是α和β的极大似然估计值；

是在假设H₀:α=α₀成立，参数β的条件极大似然估计，它由对lnL(α,β₀)关于β施行极大化来求得，令

得到方程：

$\frac{n}{\mathrm{\&beta;}}-n\ln {\mathrm{\&alpha;}}_0+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\ln t_i-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\frac{t_i}{{\mathrm{\&alpha;}}_0}\right)}^{\mathrm{\&beta;}}\ln \left(\frac{t_i}{{\mathrm{\&alpha;}}_0}\right)=0---\left(22\right)$

式（22）中，t_i(i=1,2,…,n)为关键功能部件第i次故障相对于考核起始时刻t₀的时刻；n是关键功能部件发生故障的总次数；

解方程(22)得到α₀对应的

同时

服从自由度为1的χ²分布;令X=-2lnl，解不等式

得到α₀值的集合，即在置信水平为γ时，形状参数α的置信区间；

通过公式(23)得到从尺度参数α角度考虑的MTBF的置信区间；

$\mathrm{MTBF}={\mathrm{\&alpha;}}_0\mathrm{\&Gamma;}(1+\frac{1}{\widetilde{\mathrm{\&beta;}}})---\left(23\right)$

从尺度参数α角度考虑的MTBF的置信区间与从形状参数β角度考虑的MTBF的置信区间的并集即为二参数威布尔分布模型的最终MTBF的置信区间。

所述步骤四中威布尔混合模型参数可以采用传统的图估计法获得，也可以采用传统的数学分析方法，例如极大似然估计(MLE)方法获得，还可以采用本发明提出的图解法与遗传算法相结合的参数估计方法获得。

所述图解法与遗传算法相结合的参数估计方法如下：

首先，在WPP图中散点的两端作两条渐近线，然后利用图解法求出α₁、β₁、α₂和β₂的观测值α₁₀、β₁₀、α₂₀和β₂₀；

其次，给定宽松系数k(0≤k≤1)，定义遗传算法搜索区间为

$\left\{\begin{array}{c}{\mathrm{\&alpha;}}_1\&Element;[(1-k)\&CenterDot;{\mathrm{\&alpha;}}_{10},(1+k)\&CenterDot;{\mathrm{\&alpha;}}_{10}];\\ {\mathrm{\&beta;}}_1\&Element;[(1-k)\&CenterDot;{\mathrm{\&beta;}}_{10},(1+k)\&CenterDot;{\mathrm{\&beta;}}_{10}];\\ {\mathrm{\&alpha;}}_2\&Element;[(1-k)\&CenterDot;{\mathrm{\&alpha;}}_{20},(1+k)\&CenterDot;{\mathrm{\&alpha;}}_{20}];\\ {\mathrm{\&beta;}}_2\&Element;[(1-k)\&CenterDot;{\mathrm{\&beta;}}_{20},(1+k)\&CenterDot;{\mathrm{\&beta;}}_{20}].\end{array}\right.---\left(11\right)$

再次，运用遗传算法进行参数的优化，在搜索空间内产生初始种群，进行指定次数的选择、交叉和变异，并在循环结束后提取最优解，得到α₁、β₁、α₂和β₂的估计值；

最后，将α₁、β₁、α₂、β₂的估计值代入式（10)得到可靠度函数R(t)的拐点坐标（x,y）,解方程组 $\left\{\begin{array}{c}y=\mathrm{pexp}[-{\left(\frac{x}{{\mathrm{\&alpha;}}_1}\right)}^{{\mathrm{\&beta;}}_1}]+\mathrm{qexp}\left[{\left(\frac{x}{{\mathrm{\&alpha;}}_2}\right)}^{{\mathrm{\&beta;}}_2}\right]\\ p+q=1\end{array}\right.$ 得到p和q的估计值。

所述步骤五中在二重威布尔混合分布建模后还可以对得到的二重威布尔混合分布模型进行假设检验。

所述二重威布尔混合分布模型假设检验采用d检验法，具体如下：

将D_n与临界值D_n,α进行比较，若满足式(12)，则接受原假设；否则，拒绝原假设；

$D_n=\underset{-\&infin;<t<+\&infin;}{\sup }|F_n\left(t_i\right)-F_0\left(t_i\right)|=\max \left\{d_i\right\}<D_{n,\mathrm{\&alpha;}}---\left(12\right)$

式(12)中：

F₀(t)——原假设分布函数，F₀(t)=1-R(t)；

F_n(t)——样本大小为n的经验分布函数，采用中位秩法，有

$F_n\left(t\right)=\left\{\begin{array}{cc}0,& t<t_i\\ \frac{i-0.3}{n+0.4},& t_i\&le;t<t_{i+1}\\ 1,& t\&GreaterEqual;t_n\end{array}\right.---\left(13\right)$

式(12)、式(13)中，t_i是将n个实验数据由小到大的排列序列，i=1,2,3...(n-1),n为样本大小，即关键功能部件发生故障的总次数。

本发明利用改进的多重威布尔分布模型能够更好揭示功能部件产品寿命周期过程中故障分布规律，通过不同故障期的真实分布及其衔接关系，探寻现场运行数据内在的相互依赖关系，可以进而建立整机及其子系统的早期故障期、偶然故障期及其相互衔接的寿命分布模型（包括故障间隔时间概率密度模型、累积概率分布模型、可靠度模型和故障率模型），这是评估整机产品的可靠性水平和求解各种可靠性参数的理论依据，也是进一步分析产品故障规律、预测故障发展、研究其失效机理及制定可靠性提高策略的重要手段。本发明是基于运行状态的数控机床动力伺服刀架、刀库及机械手、高速主轴单元等功能部件的可靠性现场试验数据进行功能部件的可靠性建模和评估研究，采用多重威布尔分布模型揭示这些功能部件产品寿命周期过程中故障分布规律，并提出了综合参数估计方法，较好地解决了功能部件这类复杂系统的可靠性建模和评估问题。运用该种方法对数控机床功能部件进行系统的可靠性建模与评估研究在国内外公开发表的中英文文献中未见报道。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细说明。

图1本发明的数控机床三类功能部件的可靠性建模与可靠性评估方法的流程图。

图2具体实施实例中圆盘式刀库系统的WPP图。

图3具体实施实例中圆盘式刀库系统的WPP图的渐近线拟合图。

图4具体实施实例中圆盘式刀库间隔时间可靠度函数曲线图。

图5具体实施实例中圆盘式刀库间隔时间累积故障分布函数曲线图。

具体实施方式

数控机床功能部件可靠性分析的首要问题在于寻找确切反映功能部件产品失效机理并与故障或失效数据的分析结果相符合的寿命分布模型。为了确定数控机床功能部件这种由多个零部件组成、零部件之间具有故障相关性、零部件的故障分布不同或分布类型相同参数不同的复杂系统的可靠性模型，如果用传统的可靠性分析方法，单纯使用标准威布尔分布来描述,会有很大局限性。本发明主要解决的技术问题是提供一种基于运行状态的数控机床的动力伺服刀架、刀库及机械手、高速主轴单元等功能部件的可靠性建模和评估方法，采用多重威布尔分布模型揭示对整机可靠性有重要影响的功能部件产品寿命周期过程中故障分布规律，通过对来自运行现场的故障时间信息进行理论分析，并提出了综合参数估计方法，建立了比较符合实际的可靠性模型，由此使功能部件可靠性评估从观测值上升到理论高度，为评估整机产品的可靠性水平和求解各种可靠性参数提供理论依据。

本发明是提供一种基于运行状态的数控机床功能部件的可靠性建模和评估方法，具体实施流程图如图1所示，其步骤如下：

步骤1：进行数控机床三种功能部件的故障信息挖掘。产品可靠性强调在生产现场运转中的功能维持性。由于现代数控机床结构复杂，价格昂贵，且在使用条件下的不确定因素多，许多故障只有在用户使用中才能暴露和反映出来，所以通过机床现场实际运行考核机床功能部件的可靠性，挖掘故障信息，才能反映其在实际使用中发生故障的真实情况，因而本发明采用基于可靠性现场试验方法挖掘机床功能部件的故障信息。主要包括以下内容：

1)故障的定义、分类及计数准则。(a)故障定义在规定的条件下，规定的时间内，丧失规定功能或其性能指标超过规定界限。不能完成规定的功能，称为功能性故障，如刀库转位失调、撞刀等；一个或几个性能参数超出允许的变化范围，称为参数性故障，如加工精度超差等。

(b)故障分类：故障分为关联故障和非关联故障。关联故障是由于产品本身质量缺陷引起的，在解释试验或者计算可靠性量值时必须记入的故障。非关联故障是由于误用或维修不当以及其它外界因素所引起的，在解释试验或者计算可靠性量值时应予排出的故障。

(c)故障的计数准则。考核可靠性指标时，仅计入关联故障，但非关联故障也应记录以便分析和判断。关联故障的计数准则主要有：

●必须经更换元器件、零部件或附属设备才能排除的故障，如电源烧毁等。

●有寿命期要求的损耗件在其寿命期以内发生的故障，如密封圈在寿命期内的损坏。

●需要对电缆、接插件、印刷电路板等进行修整，以消除断路、短路和接触不良，方可排除的故障，如电气系统按钮线断。

●机床偶然出现失常或停运现象，但再次起动就能恢复正常工作；这种偶然事件如在72小时内累计达三次应计为一次关联性故障；不足三次则可作为非关联性故障处理。

●不是同一因素引起而同时发生的关联性故障，则应如数计入；如果是同一因素引起的，则只计一次。

●修复后累计工作不足24小时的，发生同一关联性故障，只计入一次。

2）可靠性信息的采集和故障记录规范。采集和积累故障数据是可靠性工作的重要组成部分。为确保原始信息本身的可靠性，与主机生产厂及其典型用户协同制订规范的《关键功能部件机床运行记录表》和《关键功能部件故障记录表》，如表1和表2所示。运行记录表记载考核机床个功能部件工作时间和停机时间，故障记录表对发生每一次故障都要进行详细记录。机床运行过程中，无论是否发生故障，每个班次下班时，操作工都要填写运行记录表。如果某个班次机床出现故障，经过维修，机床恢复运行后，维修工人按照故障记录表所列内容逐项进行故障记录。

表1关键功能部件机床运行记录表

表2关键功能部件故障记录表

修理工设备使用人

表中：工时是指台时与每台维修人数的乘积。

3）故障数据的分类和整理。定时结尾的现场试验结束后，将《运行记录表》和《故障记录表》收集起来，按照故障判据和故障类型进行故障数据的分类和整理，得到试验样本数据故障汇总表。根据制定的定时截尾实验方案对关键功能部件故障记录表进行处理。记录每台关键功能部件故障发生的时刻(即故障停机月日)，令考核起始时刻t=0，记为t₀，再根据故障发生的实际时间计算出关键功能部件每个故障对于起始时刻t₀的时刻t_i(i=1,2,…,n)。根据故障时刻得到的故障间隔时间，进行故障模型识别。

绘制故障间隔时间的WPP图(Weibull Probability Paper)，按增序排列数据，令t_i(i=1,2,…,n)标记这列有序数据。WPP图横坐标x_i及纵坐标y_i如下：

$\left\{\begin{array}{c}{x}_i=\ln \left(t_i\right)\\ y_i=\ln [-\ln (1-\frac{i-0.3}{n+0.4})]\end{array}\right.---\left(1\right)$

式（1）中，t_i(i=1,2,…,n)为功能部件第i次故障相对于考核起始时刻t₀的时刻；n是功能部件发生故障的总次数。

在普通坐标纸上绘制(x_i,y_i)，i=1,2,…,n。

步骤2：判断故障间隔时间WPP图是否具有明显的拐点，如果WPP图趋向于一条直线，则对其进行简单的二参数威布尔分布建模分析；如果具有明显的拐点，则对其采用威布尔混合模型建模分析。

步骤3：功能部件寿命分布模型的建立。主要按以下步骤进行：

1如果WPP图趋向于一条直线，则对其进行简单的二参数威布尔分布建模分析。简单的二参数威布尔分布建模、检验和评估均采用如下方法：

1.1威布尔分布

威布尔分布的概率密度函数为：

$f\left(t\right)=\left\{\begin{array}{cc}\frac{\mathrm{\&beta;}}{\mathrm{\&alpha;}}{\left(\frac{t-\mathrm{\&gamma;}}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}-1}\exp [-{\left(\frac{t-\mathrm{\&gamma;}}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}}],& t\&GreaterEqual;\mathrm{\&gamma;}\\ 0,& t<\mathrm{\&gamma;}\end{array}\right.---\left(2\right)$

累积故障分布函数为：

$F\left(t\right)=\left\{\begin{array}{cc}{\&Integral;}_0^{t}f\left(t\right)\mathrm{dt}=1-\exp [-{\left(\frac{t-\mathrm{\&gamma;}}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}}],& t\&GreaterEqual;\mathrm{\&gamma;}\\ 0,& t<\mathrm{\&gamma;}\end{array}\right.---\left(3\right)$

在式（2）、式（3）中，t为时间变量，β为形状参数，β>0；α为尺度参数，α>0；γ为位置参数，γ>0。在产品的故障分析中，β与产品的故障机理相联系，不同的β值伴随着不同的故障机理。当β<1时，呈早期故障期的寿命分布；当β=1时，呈偶然故障期的寿命分布；当β>1时，呈耗损故障期的寿命分布。α与工作条件的负载有关，负载大，则相应的α小；反之亦然。γ的变化影响概率密度曲线的平移位置，产品在t=γ之前不发生故障，在t=γ以后发生故障。

在实际应用中，往往假设在t=0时产品便发生故障。这样，式(2)、(3)便分别简化为：

$f\left(t\right)=\frac{\mathrm{\&beta;}}{\mathrm{\&alpha;}}{\left(\frac{t}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}-1}\exp [-{\left(\frac{t}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}}],t\&GreaterEqual;0$     (4)

$F\left(t\right)=1-\exp [-{\left(\frac{t}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}}],t\&GreaterEqual;0$     (5)

此处以二参数威布尔分布来研究无拐点的关键功能部件的故障间隔时间的分布规律。

1.2二参数威布尔分布的参数估计

威布尔分布参数估计常用的主要有极大似然估计、最小二乘法估计和图估计等几种。在这诸多的参数估计方法中，极大似然估计方法的精度相对较高，因此使用极大似然估计方法对参数进行估计。

公式(4)和(5)分别为两参数威布尔分布的密度函数和累积故障分布函数。其似然函数为：

$L(t;\mathrm{\&alpha;},\mathrm{\&beta;})=\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}\left(\frac{\mathrm{\&beta;}}{\mathrm{\&alpha;}}\right){\left(\frac{t_i}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}-1}\exp [-{\left(\frac{t_i}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}}]---\left(6\right)$

似然方程组经整理化简后为：

$\left\{\begin{array}{c}\frac{1}{\mathrm{\&beta;}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i^{\mathrm{\&beta;}}\ln t_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i^{\mathrm{\&beta;}}}-\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\ln t_i\\ {\mathrm{\&alpha;}}^{\mathrm{\&beta;}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i^{\mathrm{\&beta;}}\end{array}\right.---\left(7\right)$

式（6）、（7）中，n和t_i的含义与式（1）中相同。

通过迭代法即可求出α和β的估计值

和

将和

值代入式（4）、式（5）即可得到二参数威布尔分布模型（包括两参数威布尔分布的概率密度函数f(t)和累积故障分布函数F（t））。

1.3二参数威布尔分布拟合的假设检验

分布类型的判断原则上有理论法和统计法。理论法是根据失效机理制定的数学模型或者某种分布的性质推导出来的。统计法是根据大量试验统计求得的。很多同类性能在以往大量试验的基础上已经验证了其分布。对分布不明的情况，则应做大样本的试验判定分布类型，对已有经验参考的，则可做较小样本的试验，假设其分布类型，再进行相应的拟合性检验。常用的统计法有χ²检验法和d检验法。

χ²检验法一般只用于大样本，而且对于截尾样本容易犯第Ⅱ类错误，即接受了不正确的原假设。所以，在综合比较的基础上，本发明对以上所推导的关键功能部件故障间隔时间分布函数（即累积故障分布函数）进行d检验。

d检验法比χ²检验法精细，而且还适用于小样本的情况。d检验法是将n个试验数据按由小到大的次序排列，根据假设的分布，计算每个数据对应的F₀(t_i)，将其与经验分布函数F_n(t_i)进行比较，其中差值的最大绝对值即检验统计量D_n的观察值。将D_n与临界值D_n，α进行比较。满足下列条件，则接受原假设，否则拒绝原假设。

将D_n与临界值D_n,α进行比较，满足式（8），则接受原假设，否则拒绝原假设；

$D_n=\underset{-\&infin;<t<+\&infin;}{\sup }|F_n\left(t_i\right)-F_0\left(t_i\right)|=\max \left\{d_i\right\}<D_{n,\mathrm{\&alpha;}}---\left(8\right)$

式(8)中：

F₀(t_i)——原假设分布函数，即式（5）中的累积故障分布函数F(t)；

F_n(t)——样本大小为n的经验分布函数，采用中位秩法，有

$F_n\left(t\right)=\left\{\begin{array}{cc}0,& t<t_i\\ \frac{i-0.3}{n+0.4},& t_i\&le;t<t_{i+1}\\ 1,& t\&GreaterEqual;t_n\end{array}\right.---\left(9\right)$

D_n,α——临界值；

临界值根据显著水平，查下表可得。下表中α为显著水平。

用式（5）中的累积故障分布函数F(t)作为原假设分布函数F₀(x_i)即可对F(t)进行检验。

2如果WPP图具有明显的拐点，则对其采用二重威布尔混合模型建模分析。

2.1二重威布尔混合模型

二重威布尔混合模型的可靠度函数R(t)是：

$R\left(t\right)=p{R}_1\left(t\right)+q{R}_2\left(t\right)$

$=\mathrm{pexp}[-{\left(\frac{t}{{\mathrm{\&alpha;}}_1}\right)}^{{\mathrm{\&beta;}}_1}]+\mathrm{qexp}[-{\left(\frac{t}{{\mathrm{\&alpha;}}_2}\right)}^{{\mathrm{\&beta;}}_2}]---\left(10\right)$

式（10)中：

t——时间变量；

α₁和β₁——第一个子分布的形状参数和尺度参数；

α₂和β₂——第二个子分布的形状参数和尺度参数；

p和q——分别为第一个和第二个子分布在二重威布尔混合模型中的比重，并且0≤p≤1，0≤q≤1，p+q＝1；

对于整个关键功能部件，存在两类失效：突发性失效，这种失效模式的形状参数β<1；磨耗型失效，这种失效模式的形状参数β>1。当α₁=α₂，β₁=β₂时，该模型将成为一个简单的威布尔模型。

二重威布尔混合模型最常用的方法是图解法。传统参数的估计方法都是使用图估计法，将故障数据绘制成故障间隔时间WPP图（即威布尔概率图），然后通过数据曲线的渐近线，求得α₁、β₁、α₂、β₂之后得到函数的拐点坐标（x,y）,利用解方程组 $\left\{\begin{array}{c}y=\mathrm{pexp}[-{\left(\frac{x}{{\mathrm{\&alpha;}}_1}\right)}^{{\mathrm{\&beta;}}_1}]+\mathrm{qexp}\left[{\left(\frac{x}{{\mathrm{\&alpha;}}_2}\right)}^{{\mathrm{\&beta;}}_2}\right],\\ p+q=1\end{array}\right.$ 便可得到p和q的值。图解法在大多数情况下可以获得参数估计的可行解，但主观性很强，精度不高。对于式(10)而言，由于参数相对较多(5个)，传统的数学分析方法，例如极大似然估计(MLE)，通常比较复杂，而且很难获得最终解。本发明提出一种图解法和遗传算法相结合的参数估计方法。

遗传算法(Genetic Algorithm，GA)具有较强的全局优化能力，是一种自适应的、智能的搜索技术，其最成功的应用领域是复杂的非线性优化问题。遗传算法最关键的操作是适应度函数的选取，它影响到收敛的速度和解的精度。另外，合适的搜索空间的选择有助于更快的找到满意的解。本发明利用图解法为遗传算法提供初始解和搜索空间。下面给出参数估计的具体步骤。

首先作故障间隔时间的WPP图。利用本领域最常规的方法——图解法求出α₁、β₁、α₂和β₂的观测值α₁₀、β₁₀、α₂₀和β₂₀。具体方法是：在图中散点的两端作两条渐近线，并标记为L₁和L₂。根据L₁和L₂的表达式，然后由最小二乘法分别求得四个参数的观测值α₁₀、β₁₀、α₂₀和β₂₀。

给定宽松系数k(0≤k≤1)，定义遗传算法搜索区间为

$\left\{\begin{array}{c}{\mathrm{\&alpha;}}_1\&Element;[(1-k)\&CenterDot;{\mathrm{\&alpha;}}_{10},(1+k)\&CenterDot;{\mathrm{\&alpha;}}_{10}];\\ {\mathrm{\&beta;}}_1\&Element;[(1-k)\&CenterDot;{\mathrm{\&beta;}}_{10},(1+k)\&CenterDot;{\mathrm{\&beta;}}_{10}];\\ {\mathrm{\&alpha;}}_2\&Element;[(1-k)\&CenterDot;{\mathrm{\&alpha;}}_{20},(1+k)\&CenterDot;{\mathrm{\&alpha;}}_{20}];\\ {\mathrm{\&beta;}}_2\&Element;[(1-k)\&CenterDot;{\mathrm{\&beta;}}_{20},(1+k)\&CenterDot;{\mathrm{\&beta;}}_{20}].\end{array}\right.---\left(11\right)$

运用遗传算法进行参数的优化。在搜索空间内产生初始种群，进行指定次数的选择、交叉和变异，并在循环结束后提取最优解。

宽松系数k的值根据具体情况所得，一般取0.2至0.8，当然也可以根据最大的适应度来取。一般来说，也可以先去一个较小值，参与遗传算法，得到最终的最优适应度值，然后，逐渐增加，根据每次得到的最终的最优适应度值，确定k的值。

对于遗传算法的初始种群的产生、交叉和变异等具体操作，方法都是基本固定的，属于本领域的常规技术手段，因此不做过多说明。适应度函数是因具体情况而异的，现在d检验(柯尔莫哥洛夫-斯米尔诺夫检验，或k-s检验)的基础上进行适当的处理作为适应度函数。

d检验法是可靠性分析中常用的检验方法，它比χ²检验法精细，而且还适用于小样本的情况，d检验法是将n个数据按照由小到大的次序排列，根据假设的分布，计算每个数据对应的F₀(t_i)，将其与经验分布函数F_n(t_i)进行比较，其中差值的最大绝对值即检验统计量D_n的观察值。将D_n与临界值D_n,α进行比较。若满足式(12)，则接受原假设；否则，拒绝原假设。

将D_n与临界值D_n，α进行比较，若满足式(12)，则接受原假设；否则，拒绝原假设；

$D_n=\underset{-\&infin;<t<+\&infin;}{\sup }|F_n\left(t_i\right)-F_0\left(t_i\right)|=\max \left\{d_i\right\}<D_{n,\mathrm{\&alpha;}}---\left(8\right)$

式(12)中：

F₀(t)——原假设分布函数，这里，它是二重威布尔混合模型的分布函数，函数中的参数对应于种群中个体的染色体段。

F_n(t)——样本大小为n的经验分布函数，采用中位秩法，有

$F_n\left(t\right)=\left\{\begin{array}{cc}0,& t<t_i\\ \frac{i-0.3}{n+0.4},& t_i\&le;t<t_{i+1}\\ 1,& t\&GreaterEqual;t_n\end{array}\right.---\left(13\right)$

D_n，α——临界值，取显著性水平α=0.10，则由经验公式得：

本发明采用的适应度函数的表达式如下

Fitness(x)=D_n-D_n,α=max{|F_n(x_i)-F₀(x_i)|}-D_n,α    (14)

用F（t)=1-R(t)（R(t)即式（10）中的可靠度函数）作为原假设分布函数F₀(t_i)对R(t)进行检验。

适应度函数的取值越小，表明假设函数越贴近经验分布函数，即获得的假设分布的函数越精确”，在遗传迭代中被选择的几率越大。而且这个适应度函数还有一个重要的优势：如果对于某个体使Fitness(x)<0，则可以判定该个体中的参数估计值通过d检验，可以在遗传迭代中直接获得符合d检验的参数估计值。

步骤3：功能部件的可靠性评估。主要按以下步骤进行：

1二参数威布尔分布模型可靠性评估

1.1MTBF的观测值

故障间隔时间是指关键功能部件相邻两次故障间的工作时间。相应的可靠性指标为平均故障工作时间，用MTBF(Mean Time Between Failures)表示，它是故障间隔时间T的数学期望E(t)或记为

关键功能部件MTBF的观测值按下述公式进行计算：

$\mathrm{MTBF}=\frac{1}{N_0}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{r}_i}---\left(15\right)$

式(15)中：

N₀—在评定周期内关键功能部件累计故障频数；

n—关键功能部件抽样台数；

t_i—在评定周期内第i台关键功能部件的实际工作时间，单位是小时；

r_i在评定周期内第i台关键功能部件出现的故障频数。

1.2MTBF的点估计值

对于二参数威布尔分布模型，其MTBF的点估计值按下式计算：

$\mathrm{MTBF}={\&Integral;}_0^{\&infin;}\mathrm{tf}\left(t\right)\mathrm{dt}=\mathrm{\&alpha;}\&CenterDot;\mathrm{\&Gamma;}(1+\frac{1}{\mathrm{\&beta;}})---\left(16\right)$

式中：f(t)为故障间隔时间的概率密度函数。将式（7）通过迭代法计算出的α，β代入式（16）即可求出MTBF的点估计值。

所谓MTBF的观测值就是不考虑故障数据的分布规律，只是从观测或者测定的角度，人工测定得到的结果。通过测量或测定所得到的样本值。

MTBF的点估计是从故障数据的统计学意义考虑，所谓点估计就是由样本x₁，x₂，…x_n确定一个统计量，用它来估计总体的未知参数，称为总体参数的估计量。当具体的样本抽出后，可求出样本参数的值。用它做为总体参数的估计值，称做总体参数的点估计，实际上它就是总体未知参数的近似值。这样能够得到更多的信息，可以知道被测设备是处于何种时期，设备是处于何种状态。

MTBF的观测值和MTBF的点估计都能够从一个侧面反映MTBF的水平，也就是可靠性的水平，但是都不完整。MTBF的观测值反映的是试验期间设备可靠性水平的观测水平，而MTBF的点估计反映试验期间设备可靠性水平的真实水平。

如果MTBF的观测值与MTBF的点估计数值接近，一方面意味着本专有技术得到的二参数威布尔分布模型与实际分布函数比较相符，同时也证明是实验数据的真实可靠。

1.3区间估计

根据似然比（Likelihood Ratio）理论，待估计参数θ₁和θ₂的似然函数L(θ₁,θ₂)与极大似然函数

的比值的自然对数的负2倍服从置信水平为α，且自由度为1的χ²分布，即：

$-2\ln l=-2\ln \left(\frac{L({\mathrm{\&theta;}}_1,{\mathrm{\&theta;}}_2)}{L({\widehat{\mathrm{\&theta;}}}_1,{\widehat{\mathrm{\&theta;}}}_2)}\right)=-2\ln L({\mathrm{\&theta;}}_1,{\mathrm{\&theta;}}_2)+2\ln L({\widehat{\mathrm{\&theta;}}}_1,{\widehat{\mathrm{\&theta;}}}_2)~{\mathrm{\&chi;}}_{\mathrm{\&alpha;}}^2\left(1\right)$

通过该理论可以求得尺度参数和形状参数的置信区间。

二参数威布尔分布的对数似然函数为：

$\ln L(\mathrm{\&alpha;},\mathrm{\&beta;})=n\ln \mathrm{\&beta;}-\mathrm{n\&beta;}\ln \mathrm{\&alpha;}+(\mathrm{\&beta;}-1)\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\ln t_i-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\frac{t_i}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}}---\left(17\right)$

式（17）中，n和t_i的含义与式（1）中相同。

形状参数β的置信区间

原假设H₀:β=β₀，备择假设H:β≠β₀，则似然比为：

$l=\frac{L(\widetilde{\mathrm{\&alpha;}},{\mathrm{\&beta;}}_0)}{L(\widehat{\mathrm{\&alpha;}},\widehat{\mathrm{\&beta;}})}---\left(18\right)$

式（18）中：

和

是α和β的极大似然估计值；

是在假设H₀:β=β₀成立，参数α的条件极大似然估计，它可由对lnL(α,β₀)关于α施行极大化来求得：

$\widetilde{\mathrm{\&alpha;}}={\left(\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i^{\mathrm{\&beta;}}\right)}^{\frac{1}{\mathrm{\&beta;}}}---\left(19\right)$

同时

服从自由度为1的χ²分布。因此-2lnl是关于β₀的函数。设置信水平为γ，即其显著水平为1-γ，令X=-2lnl。因为X～χ²(1)（X～χ²(1)比表示，X服从自由度为1的卡方分布（又叫χ²分布））。

所以 ${\&Integral;}_0^{{\mathrm{\&chi;}}_{\mathrm{\&gamma;}}^2\left(1\right)}\frac{X^{-\frac{1}{2}}\exp (-\frac{X}{2})}{\sqrt{2}\mathrm{\&Gamma;}\left(\frac{1}{2}\right)}\mathrm{dX}=\mathrm{\&gamma;}$ 即 $X=-2\ln l{\&le;\mathrm{\&chi;}}_{\mathrm{\&gamma;}}^2\left(1\right)$ 成立时，接受原假设H₀。解此不等式就可以得到β₀值的集合，即为在置信水平为γ时，形状参数β的置信区间。再通过公式：

$\mathrm{MTBF}=\widetilde{\mathrm{\&alpha;}}\mathrm{\&Gamma;}(1+\frac{1}{{\mathrm{\&beta;}}_0})---\left(20\right)$

就可以得到通过β₀值计算出的MTBF的集合。

尺度参数α的置信区间

原假设H₀:α=α₀，备择假设H:α≠α₀，则似然比为：

$l=\frac{L({\mathrm{\&alpha;}}_0,\widetilde{\mathrm{\&beta;}})}{L(\widehat{\mathrm{\&alpha;}},\widehat{\mathrm{\&beta;}})}---\left(21\right)$

式中：

和

是α和β的极大似然估计值；

是在假设H₀:α=α₀成立，参数β的条件极大似然估计，它可由对lnL(α,β₀)关于β施行极大化来求得，令

得到方程：

$\frac{n}{\mathrm{\&beta;}}-n\ln {\mathrm{\&alpha;}}_0+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\ln t_i-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\frac{t_i}{{\mathrm{\&alpha;}}_0}\right)}^{\mathrm{\&beta;}}\ln \left(\frac{t_i}{{\mathrm{\&alpha;}}_0}\right)=0---\left(22\right)$

式（22）中，n和t_i的含义与式（1）中相同。

解方程(22)即可得到α₀对应的

同时

服从自由度为1的χ²分布。因此同样地，

成立时，接受原假设H₀。解此不等式就可以得到α₀值的集合，即为在置信水平为γ时，形状参数α的置信区间。再通过公式：

$\mathrm{MTBF}={\mathrm{\&alpha;}}_0\mathrm{\&Gamma;}(1+\frac{1}{\widetilde{\mathrm{\&beta;}}})---\left(23\right)$

就可以得到通过α₀值计算出的MTBF的集合。

MTBF的置信区间

上面由式(20)可以得到从尺度参数α的角度考虑得到MTBF的置信区间，由式(23)从形状参数β的角度考虑得到MTBF的置信区间，两者的并集即为二参数威布尔分布模型的最终MTBF的置信区间；例如，由式(20)得到MTBF的区间为[900,1000]同时由式(23)得到MTBF的区间为[950,1050]，则最终MTBF的置信区间为[900,1050]。

2二重威布尔混合模型可靠性评估

2.1MTBF的观测值

二重威布尔混合模型的MTBF的观测值与二参数威布尔分布的MTBF的观测值的公式相同，具体执行过程见1.1。

2.2MTBF的点估计值

对于二重威布尔混合模型，其MTBF的点估计值按下式计算：

$\mathrm{MTBF}=p\&CenterDot;{\mathrm{\&alpha;}}_1\&CenterDot;\mathrm{\&Gamma;}(1+\frac{1}{{\mathrm{\&beta;}}_1})+q\&CenterDot;{\mathrm{\&alpha;}}_2\&CenterDot;\mathrm{\&Gamma;}(1+\frac{1}{{\mathrm{\&beta;}}_2})---\left(24\right)$

式(24)中：

α₁和β₁——第一个子分布的形状参数和尺度参数；

α₂和β₂——第二个子分布的形状参数和尺度参数；

p和q——分别为第一个和第二个子分布在二重威布尔混合模型中的比重，并且0≤p≤1，0≤q≤1，p+q＝1

2.3MTBF的区间估计

MTBF的区间估计是根据故障时间数据求得可靠性特征量MTBF的一个置信区间，这个区间以一定的概率(即置信水平)包括未知参数MTBF的真值。由于寿命试验数据的随机性，试验总时间T^*可以看作是一个随机变量。根据国标GB5080.4—85，定时截尾的可靠性特征量区间估计，选择置信水平1-α=90%，即α=10%，可找到两个下侧分位数

和

使得：

$P[{\mathrm{\&chi;}}_{\frac{\mathrm{\&alpha;}}{2}}^2\left(2r\right)\&le;\frac{2T^{*}}{\mathrm{\&theta;}}\&le;{\mathrm{\&chi;}}_{1-\frac{\mathrm{\&alpha;}}{2}}^2(2r+2)]=1-\mathrm{\&alpha;}$

式中：r—发生故障的次数；

T^*—定时截尾试验总试验时间

θ—未知参数，此处代表MTBF

设用θ^*作为未知参数θ的估计值，误差不超过某一正数δ的概率P

P(θ^*-θ|<δ)=1-α    (25)

式（25）中：

α——显著性水平

1-α——参数θ位于随机区间(θ^*-δ,θ^*+δ)的概率。

(θ^*-δ,θ^*+δ)——置信区间，表示估计的精确性。

定时截尾寿命试验时平均寿命的置信区间

采用定时截尾的试验指标区间估计，置信水平1-α=90%，即α=10%，可找到两个分位数

和

使得

$P[{\mathrm{\&chi;}}_{\frac{\mathrm{\&alpha;}}{2}}^2\left(2r\right)\&le;\frac{2T^{*}}{\mathrm{\&theta;}}\&le;{\mathrm{\&chi;}}_{1-\frac{\mathrm{\&alpha;}}{2}}^2(2r+2)]=1-\mathrm{\&alpha;}---\left(26\right)$

即：

$P[\frac{2T^{*}}{{\mathrm{\&chi;}}_{1-\frac{\mathrm{\&alpha;}}{2}}^2(2r+2)}\&le;\mathrm{\&theta;}\&le;\frac{2T^{*}}{{\mathrm{\&chi;}}_{\frac{\mathrm{\&alpha;}}{2}}^2\left(2r\right)}]=1-\mathrm{\&alpha;}---\left(27\right)$

由此可得平均无故障工作时间MTBF的置信水平为90%的双侧置信区间为：”

$\frac{2T^{*}}{{\mathrm{\&chi;}}_{0.95}^2(2r+2)}\&le;\mathrm{\&theta;}\&le;\frac{2T^{*}}{{\mathrm{\&chi;}}_{0.05}^2\left(2r\right)}---\left(28\right)$

式中：

r——发生故障的次数

T^*——定时截尾试验总试验时间

（3）有益效果：本发明与现有技术相比具有如下特点。

一、针对高速主轴单元、刀库及机械手、动力伺服刀架是国产数控机床可靠性的薄弱环节，是影响整机可靠性水平的关键部件的实际，首次对三大功能部件进行系统的可靠性建模和评估研究。以往关于数控车床功能部件的可靠性研究主要集中在对某一功能部件的故障现象分析，以及功能部件的故障对数控机床主机可靠性的影响，缺乏对功能部件进行系统的可靠性建模和评估研究，特别是对高速主轴单元、刀库及机械手、动力伺服刀架三大功能部件进行全寿命周期的可靠性建模和评估，揭示三大功能部件产品在不同寿命周期过程中故障分布的规律，这在国内外公开发表的中英文文献中还未见报道。该研究成果是评估整机产品的可靠性水平和求解各种可靠性参数的理论依据，也是进一步分析产品故障规律、预测故障发展、研究其失效机理及制定可靠性提高策略的重要手段。

二、现有机床可靠性建模和评估技术大多是针对单个部件，即使是多部件系统也是假设各部件相互独立且故障分布相同的情况，建模常常是将部件的失效数据拟合成简单的二参数或三参数威布尔分布。由于数控机床功能部件是由多个零部件构成的复杂系统，零部件之间具有故障相关性，而且其故障可能是在多种失效机理共同作用下发生的,并且在不同的寿命阶段,不同的失效机理对系统的失效起到不同的主导作用。本发明利用改进的多重威布尔分布模型可以更好揭示功能部件产品寿命周期过程中故障分布规律，通过不同故障期的真实分布及其衔接关系，建立早期故障期、偶然故障期及其相互衔接的寿命分布模型，避免单纯使用标准威布尔分布来描述,会带来较大的误差的弊端，与现有技术相比更加符合实际工程情况。

三、由于产品可靠性强调在生产现场运转中的功能维持性。因此本发明主要通过机床现场实际运行考核机床功能部件的可靠性，在实际运行过程中挖掘故障信息，相比其他现有的实验室试验技术更能反映其在实际使用中发生故障的真实情况。本发明为确保原始信息本身的可靠性，与产品生产厂及其典型用户协同制订规范的《关键功能部件机床运行记录表》和《关键功能部件故障记录表》。制订详细的可靠性跟踪考核方案，包括故障的定义、分类及计数准则以及可靠性信息的采集和故障记录规范等大量文件，保障了三大功能部件可靠性建模和评估的真实可靠。

四、本发明对三大功能部件进行可靠性建模和评估的方法，不仅局限于数控机床功能部件产品的可靠性研究，也可应用于其他类型产品的可靠性建模和评估。

下面结合对数控机床的圆盘式刀库进行可靠性建模和评估为例对本发明的具体实施方式予以进一步说明。该例根据说明书中技术方案所制定的定时截尾实验方案对刀库系统故障进行记录。从2010年1月4日起到2010年12月31日截止共收集到圆盘式刀库55条故障数据。具体的建模和评估步骤如下：

步骤1：故障信息挖掘和故障模型识别。主要按如下方法进行：

1)根据圆盘式刀库故障记录表对每台受试机床进行跟踪记录。

2)由用户负责记录故障数据。一旦发生故障，立即根据故障判据和故障类型进行记录，恢复正常工作状态后继续观察。

3)进行中途检查。每隔一定时间，生产厂家或负责此项工作的有关人员到现场了解情况，并就具体问题进行指导。

4)对在考核试验中所获得的该圆盘式刀库的基本信息和故障信息进行记录。

5)本次试验属于有替换定时截尾试验。机床的工作方式为二班生产。在现场试验结束后，将《运行记录表》和《故障记录表》收集起来，按照故障判据和故障类型进行故障数据的分类和整理，得到试验样本数据故障表如表3所示。

表3试验样本故障数据表

序号	故障间隔时间(h)	序号	故障间隔时间(h)	序号	故障间隔时间(h)
						1	42.90	20	1973.59	39	2949.66
2	107.26	21	2037.95	40	2949.66
						3	235.97	22	2070.12	41	3003.29
4	311.05	23	2145.21	42	3003.29
						5	343.23	24	2177.38	43	3003.29
6	353.96	25	2349.00	44	3003.29
						7	525.58	26	2349.00	45	3056.92
8	922.44	27	2402.63	46	3228.53
						9	933.16	28	2434.81	47	3239.26
10	1008.25	29	2434.81	48	3260.71
						11	1104.78	30	2434.81	49	3271.44
12	1104.78	31	2584.97	50	3314.34
						13	1147.68	32	2606.42	51	3550.32
14	1222.77	33	2681.51	52	3614.67
						15	1233.49	34	2724.41	53	3775.56
16	1394.38	35	2788.77	54	3786.29
						17	1512.37	36	2788.77	55	3839.92
18	1598.18	37	2788.77
						19	1823.42	38	2863.85

6）根据表1绘制故障数据的WPP图进行故障模型识别，附图2是应用Matlab绘制的圆盘式刀库系统的WPP图，图中的散点并非在一条直线上分布，呈现一个具有拐点的曲线。因此，选用二重威布尔混合模型对圆盘式刀库系统的故障间隔时间建立可靠性模型，随后进行参数估计和假设检验。附图3给出了WPP图的两条渐近线，它们的表达式为

$\left[\begin{array}{c}{L}_1:y=0.7935x-5.6021\\ L_2:y=3.5224x-27.7585\end{array}\right]$

步骤2：建立圆盘式刀库寿命分布的二重威布尔混合模型。主要按以下步骤进行：

1参数估计

应用最小二乘法求得各个参数的观测值，取宽松系数k＝0.3，可求出各个参数的搜索范围。各个参数的观测值及其搜索范围如表4所示。通过计算可得各个参数的参数估计值，见表5。适应度函数小于0，因此接受原假设，圆盘式刀库系统的可靠度函数服从以表5中对应值为参数的二重威布尔混合模型。

表4参数的观测值及其搜索范围

参数	α1	β1	α2	β2	p
						观测值	1164.43	0.7935	2645.37	3.5224
最小值	815.10	0.5555	1851.76	2.4657	0
						最大值	1513.76	1.0316	3438.97	4.5791	1

表5参数的估计值及其适应度

参数	α1	β1	α2	β2	p	适应度
							估计值	1053.56	0.9815	2921.52	4.1137	0.2971	-0.1227

因此圆盘式刀库系统的可靠度函数可以表示为

$R\left(t\right)=0.2971\&times;\exp [-{\left(\frac{t}{1053.56}\right)}^{0.9815}]+0.7029\&times;\exp [-{\left(\frac{t}{2921.52}\right)}^{4.1137}]$

从而圆盘式刀库系统的累计故障分布函数可以表示为

$F\left(t\right)=1-R\left(t\right)$

$=1-0.2971\&times;\exp \left[{-\left(\frac{t}{1053.56}\right)}^{0.9815}\right]-0.7029\&times;\exp [-{\left(\frac{t}{2921.52}\right)}^{4.1137}]$

2分布拟合的假设检验

在技术方案中已经介绍了d检验方法，现假设故障间隔时间服从二重威布尔混合分布模型：

$F\left(t\right)=1-0.2971\&times;\exp [-{\left(\frac{t}{1053.56}\right)}^{0.9815}]-0.7029\&times;\exp [-{\left(\frac{t}{2921.52}\right)}^{4.1137}]$

由计算结果知，D_n的观察值为D_n=0.0418取显著性水平α=0.10，则由经验公式得：当n=55查表

由于D_n<D_n,α，因此接收原假设，认为该产品故障间隔时间服从威布尔分布。

$F\left(t\right)=1-0.2971\&times;\exp [-{\left(\frac{t}{1053.56}\right)}^{0.9815}]-0.7029\&times;\exp [-{\left(\frac{t}{2921.52}\right)}^{4.1137}]$

3可靠性模型

由前面分析得，该刀库故障间隔时间的可靠度函数R(t)和累积故障分布函数F(t)分别为：

$R\left(t\right)=0.2971\&times;\exp [-{\left(\frac{t}{1053.56}\right)}^{0.9815}]+0.7029\&times;\exp [-{\left(\frac{t}{2921.52}\right)}^{4.1137}]$

$F\left(t\right)=1-0.2971\&times;\exp [-{\left(\frac{t}{1053.56}\right)}^{0.9815}]-0.7029\&times;\exp [-{\left(\frac{t}{2921.52}\right)}^{4.1137}]$

该圆盘式刀库系统故障间隔时间的可靠度函数R(t)和累积故障分布函数F(t)如附图4和附图5所示。

步骤3：圆盘式刀库的可靠性评估。主要按以下步骤进行：

1MTBF的观测值

根据说明书技术方案公式(15)，该系列圆盘式刀库系统的MTBF的观测值为：

$\mathrm{MTBF}=\frac{1}{N_0}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{r}_i}=217.13\left(h\right)$

2MTBF的点估计值

根据说明书技术方案公式(24)，该系列圆盘式刀库系统的MTBF的点估计值为：

$\mathrm{MTBF}=p\&CenterDot;{\mathrm{\&alpha;}}_1\&CenterDot;\mathrm{\&Gamma;}(1+\frac{1}{{\mathrm{\&beta;}}_1})+q\&CenterDot;{\mathrm{\&alpha;}}_2\&CenterDot;\mathrm{\&Gamma;}(1+\frac{1}{{\mathrm{\&beta;}}_2})$

$=0.2971\&times;1053.56\&times;\mathrm{\&Gamma;}(1+\frac{1}{0.9815})+0.7029\&times;2921.52\&times;\mathrm{\&Gamma;}(1+\frac{1}{4.1137})$

$=2179.86\left(h\right)$

3MTBF的区间估计

根据公式(28)，该系列圆盘式刀库系统的平均无故障工作时间MTBF的置信水平为90%的双侧置信区间的置信区间为：

${\mathrm{\&theta;}}_{\min }=\frac{2T^{*}}{{\mathrm{\&chi;}}_{0.95}^2(2r+2)}=\frac{2\&times;119412.86}{{\mathrm{\&chi;}}_{0.95}^2(2\&times;55+2)}=1734.37\left(h\right)$

${\mathrm{\&theta;}}_{\max }=\frac{2T^{*}}{{\mathrm{\&chi;}}_{0.05}^2\left(2r\right)}=\frac{2\&times;119412.86}{{\mathrm{\&chi;}}_{0.05}^2(2\&times;55)}=2751.71\left(h\right)$

则双侧置信区间为：(1734.37，2751.71)

Claims (6)

1.一种数控机床三类功能部件的可靠性建模与可靠性评估方法，其特征在于包括下述步骤：

步骤一：

根据数控机床可靠性现场试验数据绘制故障间隔时间的WPP图，WPP图横坐标x_i及纵坐标y_i如下：

$\left\{\begin{array}{c}{x}_i=\ln \left(t_i\right)\\ y_i=\ln [-\ln (1-\frac{i-0.3}{n+0.4})]\end{array}\right.---\left(1\right)$

式(1)中，t_i(i＝1,2,…,)为关键功能部件第i次故障相对于考核起始时刻t₀的时刻；n是关键功能部件发生故障的总次数；

步骤二：

判断故障间隔时间WPP图是否趋向于一条直线；如果WPP图趋向于一条直线，则转步骤三；如果WPP图具有明显的拐点，则转步骤四；

步骤三：

对功能部件寿命进行二参数威布尔分布建模并对其可靠性进行评估；

步骤四：

采用二重威布尔混合分布建模并对其参数进行估计；所述二重威布尔混合模型的可靠度函数R(t)如下：

$R\left(t\right)=p{R}_1\left(t\right)+q{R}_2\left(t\right)$

$=\mathrm{pexp}[-{\left(\frac{t}{{\mathrm{\&alpha;}}_1}\right)}^{{\mathrm{\&beta;}}_1}]+\mathrm{qexp}[-{\left(\frac{t}{{\mathrm{\&alpha;}}_2}\right)}^{{\mathrm{\&beta;}}_2}]---\left(10\right)$

式(10)中：

t——时间变量；

α₁和β₁——第一个子分布的形状参数和尺度参数；

α₂和β₂——第二个子分布的形状参数和尺度参数；

p和q——分别为第一个和第二个子分布在二重威布尔混合模型中的比重，并且0≤p≤1，0≤q≤1，p+q＝1；

步骤五：

对二重威布尔混合模型进行可靠性评估，关键功能部件的平均无故障时间MTBF的观测值按式(15)进行计算：

$\mathrm{MTBF}=\frac{1}{N_0}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{r}_i}---\left(15\right)$

式(15)中：

N₀—在评定周期内关键功能部件累计故障频数；

n—关键功能部件抽样台数；

t_i—在评定周期内第i台关键功能部件的实际工作时间，单位是小时；

r_i—在评定周期内第i台关键功能部件出现的故障频数；

MTBF的点估计值按式(24)计算：

$\mathrm{MTBF}=p\&CenterDot;{\mathrm{\&alpha;}}_1\&CenterDot;\mathrm{\&Gamma;}(1+\frac{1}{{\mathrm{\&beta;}}_1})+q\&CenterDot;{\mathrm{\&alpha;}}_2\&CenterDot;\mathrm{\&Gamma;}(1+\frac{1}{{\mathrm{\&beta;}}_2})---\left(24\right)$

式(24)中：

α₁和β₁——第一个子分布的形状参数和尺度参数；

α₂和β₂——第二个子分布的形状参数和尺度参数；

p和q——分别为第一个和第二个子分布在二重威布尔混合模型中的比重，并且0≤p≤1，0≤q≤1，p+q＝1；

采用定时截尾的试验指标区间估计，获得MTBF置信水平为1-α的双侧置信区间为：

其中α为显著性水平，r为关键功能部件发生故障的次数，T^*为定时截尾试验总试验时间。

2.根据权利要求1所述的数控机床三类功能部件的可靠性建模与可靠性评估方法，其特征在于所述步骤三中二参数威布尔分布建模方法如下：

威布尔分布的概率密度函数为：

$f\left(t\right)=\frac{\mathrm{\&beta;}}{\mathrm{\&alpha;}}{\left(\frac{t}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}-1}\exp [-{\left(\frac{t}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}}],t\&GreaterEqual;0---\left(4\right)$

累积故障分布函数为：

$F\left(t\right)=1-\exp [-{\left(\frac{t}{\mathrm{\&alpha;}}\right)}^{\mathrm{\&beta;}}],t\&GreaterEqual;0---\left(5\right)$

式(4)、式(5)中，t为时间变量，β为形状参数，β＞0；α为尺度参数，α＞0；

利用式(7)通过迭代法求出α和β的估计值

和

将

和值代入式(4)、式(5)得到二参数威布尔分布的概率密度函数f(t)和累积故障分布函数F(t)；

$\left\{\begin{array}{c}\frac{1}{\mathrm{\&beta;}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i^{\mathrm{\&beta;}}\ln t_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i^{\mathrm{\&beta;}}}-\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\ln t_i\\ {\mathrm{\&alpha;}}^{\mathrm{\&beta;}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{t}_i^{\mathrm{\&beta;}}\end{array}\right.---\left(7\right)$

(7)中，t_i(i＝1,2,…,n)为关键功能部件第i次故障相对于考核起始时刻t₀的时刻；n是关键功能部件发生故障的总次数。

3.根据权利要求2所述的数控机床三类功能部件的可靠性建模与可靠性评估方法，其特征在于所述步骤三中在二参数威布尔分布建模后采用d检验法对二参数威布尔分布模型进行假设检验。

4.根据权利要求2所述的数控机床三类功能部件的可靠性建模与可靠性评估方法，其特征在于所述步骤三中，根据似然比理论对二参数威布尔分布模型的可靠性进行评估。

5.根据权利要求1所述的数控机床三类功能部件的可靠性建模与可靠性评估方法，其特征在于所述步骤四中威布尔混合模型参数采用图解法与遗传算法相结合的参数估计方法获得；图解法与遗传算法相结合的参数估计方法如下：

首先，在WPP图中散点的两端作两条渐近线，然后利用图解法求出α₁、β₁、α₂和β₂的观测值α₁₀、β₁₀、α₂₀和β₂₀；

其次，给定宽松系数k(0≤k≤1)，定义遗传算法搜索区间为

$\left\{\begin{array}{c}{\mathrm{\&alpha;}}_1\&Element;[(1-k)\&CenterDot;{\mathrm{\&alpha;}}_{10},(1+k)\&CenterDot;{\mathrm{\&alpha;}}_{10}];\\ {\mathrm{\&beta;}}_1\&Element;[(1-k)\&CenterDot;{\mathrm{\&beta;}}_{10},(1+k)\&CenterDot;{\mathrm{\&beta;}}_{10}];\\ {\mathrm{\&alpha;}}_2\&Element;[(1-k)\&CenterDot;{\mathrm{\&alpha;}}_{20},(1+k)\&CenterDot;{\mathrm{\&alpha;}}_{20}];\\ {\mathrm{\&beta;}}_2\&Element;[(1-k)\&CenterDot;{\mathrm{\&beta;}}_{20},(1+k)\&CenterDot;{\mathrm{\&beta;}}_{20}].\end{array}\right.---\left(11\right)$

再次，运用遗传算法进行参数的优化，在搜索空间内产生初始种群，进行指定次数的选择、交叉和变异，并在循环结束后提取最优解，得到α₁、β₁、α₂和β₂的估计值；

最后，将α₁、β₁、α₂、β₂的估计值代入式(10)得到可靠度函数R(t)的拐点坐标(x,y),解方程组 $\left\{\begin{array}{c}y=\mathrm{pexp}[-{\left(\frac{x}{{\mathrm{\&alpha;}}_1}\right)}^{{\mathrm{\&beta;}}_1}]+\mathrm{qexp}\left[{\left(\frac{x}{{\mathrm{\&alpha;}}_2}\right)}^{{\mathrm{\&beta;}}_2}\right]\\ p+q=1\end{array}\right.$ 得到p和q的估计值。

6.根据权利要求5所述的数控机床三类功能部件的可靠性建模与可靠性评估方法，其特征在于所述步骤五中在二重威布尔混合分布建模后采用d检验法对得到的二重威布尔混合分布模型进行假设检验。

Images