非齐次泊松过程课程设计.doc

课程名称： 《随机过程》

课程设计（论文）

题 目: 非齐次泊松过程

在数控机床可靠 性建模中的应用

学 院： 理学院 专 业： 数学与应用数学 班 级： 数学12-1班 学 生 姓 名： 王玲玲 学 生 学 号： 2012027149 指 导 教 师： 蔡吉花

2015 年 1月 3 日

随机过程课程设计

目 录

任务书………………………………………………………. …………. …. …….1 摘要………………………………………………………………. ……. …. …….1 前言…………………………………………………………. ……. …………. ….2

１ 非齐次泊松过程理论 ……………………………………………………………2

1.1 非齐次泊松过程的基本理论简介…………………………………………2 1.2 基于试验总时间法的趋势检验…………………………………………… 2 2 数控机床的非齐次泊松过程可靠性建模…… ………………………………3 2.1强度函数的建立………………………………… ………. ……. ………… 3 2.2 Ｋ台数控机床强度函数的参数估计… ……………… …. ………………4 2.3 非齐次泊松过程下的可靠性指标 …………………… … …. ……. ……5 3实例分析………………………………………………………………………5 4结束语…………………………………………………………………………7 5程序及结果……………………………………………………………………8 6参考文献 ……………………………………………………………………9 附录……………………………………………………………………………… 评阅书………………………………………… …………………………………

摘 要

基于试验总时间法对多样本随机截尾的数控机床现场数据进行趋势检验，在故障过程为浴盆曲线的趋势条件下构建了数控机床的非齐次泊松过程的可靠性模型。本文使用极大似然估计法对非齐次泊松过程的强度函数进行参数估计得到了该模型的可靠性指标，以6台加工中心的现场数据为例建立了非齐次泊松过程的可靠性模型。再通过matlab 曲线拟合，绘制出故障时间的曲线，通过曲线的拟合程度，可以确定非齐次泊松过程能够更恰当地表现故障的趋势。

关键词：数控机床 可靠性 非齐次泊松过程 浴盆曲线

前言

数控机床是由数目众多的零部件组成的复杂机电液可修系统。在其可靠性研究中，需要考虑维修活动对其可靠度的影响【1】。以往的数控机床可靠性建模方法，是将故障间隔时间视为独立同分布来分析其寿命分布，即假设维修活动是“修复

【2】

如新”而在实际生产中数控机床的维修活动是以调整或者更换一部分零部件和

元器件为主的，对于复杂的系统来说这种维修活动只能使产品恢复到正常功能维修前后可靠度并没有很大改变，因此将数控机床的维修活动视为“修复如旧”更加合理。非齐次泊松过程经常被用于建立“修复如旧”的维修策略且维修时间可忽略的可修系统可靠性模型用于模拟出现故障间隔时间的趋势【3-4】。结合数控机床的维修特点，使用非齐次泊松过程建立的可靠性模型更能贴近于复杂系统的生产实际。

本文提出了非齐次泊松过程的数控机床可靠性建模方法，并结合数控机床的失效特点，建立故障率为浴盆曲线的非齐次泊松过程可靠性模型。同时，结合具有随机截尾特点的多样本数控机床现场试验故障数据，对数控机床的可靠性进行了深入分析。

１．非齐次泊松过程理论

1.1 非齐次泊松过程的基本理论简介：

非齐次泊松过程是随机点过程的一种典型类型，当可修系统的相邻故障间隔呈现某种趋势时可以使用这种方法来描述。 非齐次泊松过程的重要参数【4】：ω(t ) 为强度函数，是非负函数；其 累 积 故 障 强 度 函 数W (t ) =⎰ω(u ) du ，表示在[0,t]中的平均故障数, 即E[N(t)]=W

0t

（t ），N[t] 代表在[0,t]出现的故障次数ｔ表示机床从观测开始后的运行时间。

βt

当强度函数为ω(t ) =() β-1时，成为威布尔过程。其中，α、β>0，α为

αα

尺度参数，β为形状参数。01, 表示不断恶化的(坏) 系统；β=1，表示系统服从指数分布。

1. ２基于试验总时间法的趋势检验：

本文采用基于试验总时间的方法，对具有多样本随机截尾现场数据的故障过程进行趋势检验【5】。

将在观测期间采集到的所有故障数据按照从大到小时间进行排序，得到t (i)

的时间序列。根据试验总时间的建模思想【6】，得到该序列的第ｉ个故障发生时的试验总时间：

T (t (i ) ) =⎰n (u ) du （1）

t (i )

随机过程课程设计

k

N =∑n i （2）

i =1

式中：n(u)表示在ｕ时刻观察到的数控机床数量，N 表示在观测期间的故障数；

，当时间序列t (i ) 的最后一个时间是n i 表示第ｉ台机床的故障数（共有Ｋ台机床）

故障数据时，n (i ) =n i -1；当时间序列t (i ) 的最后一个时间不是故障数据而是截尾数据时n (i ) =n i 。

在实际检验时同时使用Ｕ－检验、Ｊ－检验和Ｖ－检验等检验方法综合确定

。

有无趋势【7-8】其中Ｖ－检验如下：

Ｈ0：齐次泊松过程； Ｈ1：具有非单调趋势； V 1=

∑|T (t (i ) ) -

T (S ) T (S )

|-N ⨯V 2=

T (S ) 2T 2(S )

|T (t (i ) ) -|-N ⨯∑~N (0,1)

(3)

~N (0,1) (4)

V 3=

2∑log(T (S ) |2T (t (i ) ) -T (S ) |)

i =1

χ(2N )

2

~N (0,1) (5)

式中：Ｔ（Ｓ）表示总的观测时间。当｜Ｖ１｜＜ｚα／２，｜Ｖ２｜＜ｚα／２，Ｖ３＜χ２（２ｎ）时，接受Ｈ０。

一般情况下，Ｕ检验和Ｊ检验是检验具有单调趋势或齐次泊松过程和更新过程的故障数据的，而当故障数据具有非单调趋势时，则可以考虑Ｖ检验中的统计量，如表１所示。

表1 故障率和故障强度函数变化特性

２ 数控机床的非齐次泊松过程可靠性建模 2.1强度函数的建立

对于浴盆曲线趋势的故障过程，假设其故障强度函数由早期故障期和偶然故障期两部分组成，并且每一个阶段都是一个威布尔过程［９］，参数为尺度参数α

ｍ

和形状参数βｍ（ｍ＝１，２）。结合以上假设和多重威布尔分布模型的性质，则

随机过程课程设计

该数控机床故障强度函数为：

βt

ω(t ) =1() β1-1 （6）

α1α1

式中：α１，α２，β１，β２＞0

在（0，ｔ］内的平均故障个数为累积故障强度函数，即

由于该模型是由两重威布尔过程构成，其强度函数是具有非单调的浴盆曲线趋势，因此组成该模型的两个形状参数有（β１－１）（β２－１）＜0，则本文中假设β１＜１，β２＞１。

2.2 Ｋ台数控机床强度函数的参数估计

本文使用极大似然估计法对k 台样本的强度函数进行参数估计，k 台数控机床的故障数据是随机截尾的，第ｉ台的故障观测时间为［0，T ｉ］, 其中T i 为现场试验的截尾时间。t 0≡0，因此, 得到相应的似然函数为

n j K

T T β1t ij β1-1β2t ij β2-1

L K =∏{∏[() +() ]⨯exp[-(i ) β1-(i ) β2]} (8)

i =1

j =1

α1α1α2α2α1α2

似然函数的对数函数以及此对数函数对模型参数的偏导数为

n j K

T T β1t ij β1-1β2t ij β2-1

l K =∑{∑ln[() +() ]-[(i ) β1+(i ) β2]} (9)

i =1

j =1

α1α1

-(

α2α2

m

α1α2

m

βm 2t ij

) () β1-1

K n αα∂l K

=∑∑

∂αm i =1j =1β1t ij β-1β2t ij β

()() +()() α1α1α2α2

j

m 1

+

2-1

βm T i β

() (10) αm αj

由（10）可以得到

[1+βm ln(

n j

t ij

K

T T ∂l K

=∑∑-(i ) βm ln(i ) (11) ∂αm i =1j =1β1ij β1-1β2ij β2-1αm αm

()() +()()

αm αm

)](

t ij

αm

) βm -1

α1α1α2α2

m -1

∑T

i =1

K

K

βm

=∑∑

i =1j =1

K

n j

ββm t ij

(

β1β

)() β-1+(2)() βα1α1α2α2

1

ij ij

(12)

2-1

由累积故障函数可得

K K K

T i T i

∑W (T i ) =∑[() β1+(() β2)]=∑n i (13)

i =1

i =1

α1α2

i =1

以上式可导出：

随机过程课程设计

K

α2=∑

i =1

T i

[n i -(

T i

1

(14)

α1

)]β1

将式（9）转换为三参数的函数，即

n j K

βt βt

'=∑{∑ln[(1)(ij ) β1-1+(2)(ij ) β2-1]-n i (15) l K

i =1

j =1

α1α1α2α2

最终，似然函数的参数估计转化成以下的求最大化问题：

' max l K 约束条件：

一般情况，最大化问题都需要初始值。根据经验，在没有合适的初始值选择下，可以假设

2.3 非齐次泊松过程下的可靠性指标

（１）首次故障间隔时间的可靠度函数【10-11】从t ＝0开始直到第一个故障发生的时间T １,T 1的可靠度函数为

R 1(t ) =Pr(T i >t ) =e

-W (t )

=e ⎰0

-ω(u ) du

t

(16)

对于非齐次泊松过程模型的使用，如果能够估算出首次故障间隔时间的故障率函数，就能同时估计出产品整个寿命的强度函数。

（２）其他故障间隔时间的可靠度函数

在t 0时刻后的可靠度函数

R (t |t 0) =Pr(t +t 0>t 0) =e

-

⎰t 0

t 0+t

ω(u ) du

(17)

（３）平均故障间隔时间

瞬时平均故障间隔时间———故障强度函数ω（t ）表示单位时间发生故障的次数，则ω(t)的倒数表示一次故障所经过的时间，定义瞬时故障平均间隔时间为

T IMTBF =

1

(18) ω(t )

累积平均故障时间间隔———表示一段时间内的平均故障间隔时间，即累积故障强度函数的倒数，

t -t t -t

T CMTBF =21=t 221 (19)

W (t 1, t 2) ω(u ) du

⎰

t 1

３实例分析

以国内北京第一机床厂同一时期出厂的６台加工中心现场试验的故障数据

为例，其发生故障的时间如表2所示。首先，需要对这些数据进行趋势检验。根据2.2节中多样本的趋势检验方法，在显著性水平α=0.05，得到这批加工中心的统计量值如表3所示。

编号

1 2 3 4 5 6

故障时间/h

50.99 423.72 753.06 760.65 795.63 1005.80 1209.40 2509.2 3350.16 3801.90 3915.62 4011.10 5109.03 5197.12 5353.92 5845.90 5942.81 6106.49 6323.63 6474.60 6526.03 6827.10 7059.69 7460.86 8240.27 8745.00 9142.65

185.13 458.00 960.54 1005.87 3409.55 422.89 6061.44 6217.53 7479.45 7542.81 7775.96 7882.88 7994.25 8588.25 28.05 350.48 47.52 1560.23 1896.30 2541.10 3352.80 3915.12 4981.45 5112.97 5729.13 5812.46 5903.40 6109.13 6117.21 6275.28 6308.78 6348.21 6457.61 6620.46 6853.44 7005.85 7116.59 7249.74 7467.90 8088.96 8298.18 9509.28 131.09 785.61 287.51 870.56 2987.45 3500.75 4881.86 5136.51 5230.01 5376.53 5540.54 5746.57 6183.21 6505.13 6592.08 7125.03 7379.46 7703.03 7868.85 8275.74 8654.42 9032.10

148.17 578.80 1014.14 1952.18 2893.01 3287.36 3747.55 4279.01 4714.12 4839.79 5558.09 5600.66 6694.61 6855.49 7120.48 7368.47 7496.84 7659.20 8451.26 8638.805

51.98 359.4 956.72 1357.45 1549.56 2706.15 3417.46 4659.60 5150.64 5206.74 5483.61 5570.40 5651.25

表２ 加工中心故障数据

在V 检验中运行结果拒绝H 0可知故障数据具有非单调趋势，且由表１可知，故障发生过程呈浴盆曲线的趋势。

非齐次泊松过程是随机点过程的一种典型类型，当可修系统的相邻故障间隔呈现某种趋势时可以使用这种方法来描述。因此，建立扩展的非齐次泊松过程模型，并且对加工中心的故障强度函数参数进行估计，通过matlab 运行结果得到尺度参数及形状参数，代入公式得到：

ω(t ) =

0.62t 1.97t

() 0.62-1+() 1.97-1 239.59239.592641.902641.90

本文在故障数据基础上，根据典型威布尔过程对数控机床进行了参数估计，并与文中模型得到的结果进行了对比。由图1中可以明显看出，本文采用的扩展的非齐次泊松过程能够更恰当地表现故障的趋势。

拟合结果如图１所示，本文方法的拟合值与估计值相关系数

由于样本量较大，且相关系数接近1，所以拟合值与估计值之间线性相关，故障过程符合假设的强度函数为浴盆曲线的威布尔过程。其相关的可靠性指标如下。

（１）首次故障间隔时间的可靠度函数

图2表示加工中心从t=0时刻投入运行后，继续无故障工作的可靠度曲线。

（２）平均故障间隔时间

图3是瞬时平均故障间隔时间和从t=0时刻起的累计平均故障间隔时间的曲线。从图中可以看出，在观测时间的初始阶段，平均故障间隔时间较高，随着

观测时间的增大，平均故障间隔时间变小，对照原始故障数据，可以看出平均故障间隔时间的趋势与现实情况相符。

４结束语

基于试验总时间法的数控机床非齐次泊松过程可靠性建模方法，不仅可以解决随机截尾故障数据趋势检验的问题，同时“修复如旧”的前提假设相对于以往的可靠性方法更加适用于数控机床的维修。6台加工中心实例的研究结果表明，本文所建立的模型更贴近于数控机床的实际运行状态和可靠性水平。

5、程序及结果：

%程序1：V 检验求出该过程具有非单调趋势

X=[50.99 423.72 753.06 760.65 795.63 1005.80 1209.40 2509.2 3350.16 3801.90 3915.62 4011.10 5109.03 5197.12 5353.92 5845.90 5942.81 6106.49 6323.63 6474.60 6526.03 6827.10 7059.69 7460.86 8240.27 8745.00 9142.65 185.13 458.00 960.54 1005.87 3409.55 422.89 6061.44 6217.53 7479.45 7542.81 7775.96 7882.88 7994.25 8588.25 28.05 350.48 47.52 1560.23 1896.30 2541.10 3352.80 3915.12 4981.45 5112.97 5729.13 5812.46 5903.40 6109.13 6117.21 6275.28 6308.78 6348.21 6457.61 6620.46 6853.44 7005.85 7116.59 7249.74 7467.90 8088.96 8298.18 9509.28 131.09 785.61 287.51 870.56 2987.45 3500.75 4881.86 5136.51 5230.01 5376.53 5540.54 5746.57 6183.21 6505.13 6592.08 7125.03 7379.46 7703.03 7868.85 8275.74 8654.42 9032.10 148.17 578.80 1014.14 1952.18 2893.01 3287.36 3747.55 4279.01 4714.12 4839.79 5558.09 5600.66 6694.61 6855.49 7120.48 7368.47 7496.84 7659.20 8451.26 8638.805 51.98 359.4 956.72 1357.45 1549.56 2706.15 3417.46 4659.60 5150.64 5206.74 5483.61 5570.40 5651.25];

disp('X');

A=[1 2;1/2 1];

[n,n]=size(A);

x=ones(n,100);

y=ones(n,100);

m=zeros(1,100);

m(1)=max(x(:,1));

x(:,2)=A*y(:,1);

m(2)=max(x(:,2));

y(:,2)=x(:,2)/m(2);

p=0.0001;i=2;k=abs(m(2)-m(1));

while k>p

i=i+1;

x(:,i)=A*y(:,i-1);

m(i)=max(x(:,i));

y(:,i)=x(:,i)/m(i);

k=abs(m(i)-m(i-1));

end

a=sum(y(:,i));

w=y(:,i)/a;t=m(i);

disp(w);

disp(t); %以下是V 检验CI=(t-n)/(n-1);

CI=(t-n)/(n-1);

RI=[0 0 0.52 0.89 1.12 1.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];

CR=CI/RI(n);

if CR

disp('接受H0!');

disp('CI=');

disp(CI);

disp('CR=');

disp(CR);

else

disp('拒绝H0!');

end

拒绝H0

%程序2:绘制累计故障曲线

t=100:100:10000;

w=(t./239.59).^(0.62)+(t/2641.9).^(1.97);

plot(t,w,'-*')

%程序3：绘制加工中心的可靠度曲线

t=100:100:10000;

w=(0.62/239.59)*(t./239.59).^(-0.38)+(1.97/2641.90)*(t/2641.9).^(0.97);

-3

%程序4：求出尺度参数及形状参数

A %6台机器故障数据

size(A)

B=sort(A);

K=1:n; % n-故障数据样本容量

%绘制累计故障曲线

N=K./7;

Y=log(N);

X=log(B);

plot(X,Y,’+’);

% ANN输入向量及目标向量

P=X;

T=Y;

pause

clc

plot(P,T,’+’);

title(‘Training Vectors’);

xlabel(‘Input vector P’);

ylabel(‘Target Vector T’);

pause

clc

[R,Q]=size(P);

[S,Q]=size(T);

% 初始化权向量

[W,B]=solvelin(P,T);

lr=2*maxlinlr(P);

df=10;

me=100;

eg=0.2;

tP=[df me eg lr];

[W,B,epochs,error]=trainwh(W,B,P,T,TP);

pause

W %威布尔过程形状参数B

exp(B) %威布尔过程强度参数

W=239.59 2414.9

B=0.62 1.97

6 参考文献：

[1]易蓓玲．可靠性与维修性工程概论［Ｍ］．北京：清华大学出版社.2010.

[2]张波．应用随机过程．中国人民大学出版社．[书 号] 7300037690． 2005 年8月．

[3]Yazhou J,Molin W,Zhixin J.Probability distribution of machining center failures[J].Reliability Engineering and System Safety.1995.

[4] Emerging Tecnologyies and Factory Automation Proceedings.1995.

[5]王智明，杨建国，王国强. 多台数控机床最小维修的可靠性评估. 哈尔滨工业大学学报. 2011

[6]面向不完全维修的数控机床可靠性评估，张根保，李冬英，刘杰《机械工程学报》2014.

[7] 胡仁胜. 实时控制软件系统的可靠性分析. 华南师范大学学报.2001.

[8]张英芝，申桂香，薛玉霞. 随机截尾数控机床故障过程［Ｊ］．吉林大学学报：工学版.2007.

评 阅 书

课程名称： 《随机过程》

课程设计（论文）

题 目: 非齐次泊松过程

在数控机床可靠 性建模中的应用

学 院： 理学院 专 业： 数学与应用数学 班 级： 数学12-1班 学 生 姓 名： 王玲玲 学 生 学 号： 2012027149 指 导 教 师： 蔡吉花

2015 年 1月 3 日

随机过程课程设计

目 录

任务书………………………………………………………. …………. …. …….1 摘要………………………………………………………………. ……. …. …….1 前言…………………………………………………………. ……. …………. ….2

１ 非齐次泊松过程理论 ……………………………………………………………2

1.1 非齐次泊松过程的基本理论简介…………………………………………2 1.2 基于试验总时间法的趋势检验…………………………………………… 2 2 数控机床的非齐次泊松过程可靠性建模…… ………………………………3 2.1强度函数的建立………………………………… ………. ……. ………… 3 2.2 Ｋ台数控机床强度函数的参数估计… ……………… …. ………………4 2.3 非齐次泊松过程下的可靠性指标 …………………… … …. ……. ……5 3实例分析………………………………………………………………………5 4结束语…………………………………………………………………………7 5程序及结果……………………………………………………………………8 6参考文献 ……………………………………………………………………9 附录……………………………………………………………………………… 评阅书………………………………………… …………………………………

摘 要

基于试验总时间法对多样本随机截尾的数控机床现场数据进行趋势检验，在故障过程为浴盆曲线的趋势条件下构建了数控机床的非齐次泊松过程的可靠性模型。本文使用极大似然估计法对非齐次泊松过程的强度函数进行参数估计得到了该模型的可靠性指标，以6台加工中心的现场数据为例建立了非齐次泊松过程的可靠性模型。再通过matlab 曲线拟合，绘制出故障时间的曲线，通过曲线的拟合程度，可以确定非齐次泊松过程能够更恰当地表现故障的趋势。

关键词：数控机床 可靠性 非齐次泊松过程 浴盆曲线

前言

数控机床是由数目众多的零部件组成的复杂机电液可修系统。在其可靠性研究中，需要考虑维修活动对其可靠度的影响【1】。以往的数控机床可靠性建模方法，是将故障间隔时间视为独立同分布来分析其寿命分布，即假设维修活动是“修复

【2】

如新”而在实际生产中数控机床的维修活动是以调整或者更换一部分零部件和

元器件为主的，对于复杂的系统来说这种维修活动只能使产品恢复到正常功能维修前后可靠度并没有很大改变，因此将数控机床的维修活动视为“修复如旧”更加合理。非齐次泊松过程经常被用于建立“修复如旧”的维修策略且维修时间可忽略的可修系统可靠性模型用于模拟出现故障间隔时间的趋势【3-4】。结合数控机床的维修特点，使用非齐次泊松过程建立的可靠性模型更能贴近于复杂系统的生产实际。

本文提出了非齐次泊松过程的数控机床可靠性建模方法，并结合数控机床的失效特点，建立故障率为浴盆曲线的非齐次泊松过程可靠性模型。同时，结合具有随机截尾特点的多样本数控机床现场试验故障数据，对数控机床的可靠性进行了深入分析。

１．非齐次泊松过程理论

1.1 非齐次泊松过程的基本理论简介：

非齐次泊松过程是随机点过程的一种典型类型，当可修系统的相邻故障间隔呈现某种趋势时可以使用这种方法来描述。 非齐次泊松过程的重要参数【4】：ω(t ) 为强度函数，是非负函数；其 累 积 故 障 强 度 函 数W (t ) =⎰ω(u ) du ，表示在[0,t]中的平均故障数, 即E[N(t)]=W

0t

（t ），N[t] 代表在[0,t]出现的故障次数ｔ表示机床从观测开始后的运行时间。

βt

当强度函数为ω(t ) =() β-1时，成为威布尔过程。其中，α、β>0，α为

αα

尺度参数，β为形状参数。01, 表示不断恶化的(坏) 系统；β=1，表示系统服从指数分布。

1. ２基于试验总时间法的趋势检验：

本文采用基于试验总时间的方法，对具有多样本随机截尾现场数据的故障过程进行趋势检验【5】。

将在观测期间采集到的所有故障数据按照从大到小时间进行排序，得到t (i)

的时间序列。根据试验总时间的建模思想【6】，得到该序列的第ｉ个故障发生时的试验总时间：

T (t (i ) ) =⎰n (u ) du （1）

t (i )

随机过程课程设计

k

N =∑n i （2）

i =1

式中：n(u)表示在ｕ时刻观察到的数控机床数量，N 表示在观测期间的故障数；

，当时间序列t (i ) 的最后一个时间是n i 表示第ｉ台机床的故障数（共有Ｋ台机床）

故障数据时，n (i ) =n i -1；当时间序列t (i ) 的最后一个时间不是故障数据而是截尾数据时n (i ) =n i 。

在实际检验时同时使用Ｕ－检验、Ｊ－检验和Ｖ－检验等检验方法综合确定

。

有无趋势【7-8】其中Ｖ－检验如下：

Ｈ0：齐次泊松过程； Ｈ1：具有非单调趋势； V 1=

∑|T (t (i ) ) -

T (S ) T (S )

|-N ⨯V 2=

T (S ) 2T 2(S )

|T (t (i ) ) -|-N ⨯∑~N (0,1)

(3)

~N (0,1) (4)

V 3=

2∑log(T (S ) |2T (t (i ) ) -T (S ) |)

i =1

χ(2N )

2

~N (0,1) (5)

式中：Ｔ（Ｓ）表示总的观测时间。当｜Ｖ１｜＜ｚα／２，｜Ｖ２｜＜ｚα／２，Ｖ３＜χ２（２ｎ）时，接受Ｈ０。

一般情况下，Ｕ检验和Ｊ检验是检验具有单调趋势或齐次泊松过程和更新过程的故障数据的，而当故障数据具有非单调趋势时，则可以考虑Ｖ检验中的统计量，如表１所示。

表1 故障率和故障强度函数变化特性

２ 数控机床的非齐次泊松过程可靠性建模 2.1强度函数的建立

对于浴盆曲线趋势的故障过程，假设其故障强度函数由早期故障期和偶然故障期两部分组成，并且每一个阶段都是一个威布尔过程［９］，参数为尺度参数α

ｍ

和形状参数βｍ（ｍ＝１，２）。结合以上假设和多重威布尔分布模型的性质，则

随机过程课程设计

该数控机床故障强度函数为：

βt

ω(t ) =1() β1-1 （6）

α1α1

式中：α１，α２，β１，β２＞0

在（0，ｔ］内的平均故障个数为累积故障强度函数，即

由于该模型是由两重威布尔过程构成，其强度函数是具有非单调的浴盆曲线趋势，因此组成该模型的两个形状参数有（β１－１）（β２－１）＜0，则本文中假设β１＜１，β２＞１。

2.2 Ｋ台数控机床强度函数的参数估计

本文使用极大似然估计法对k 台样本的强度函数进行参数估计，k 台数控机床的故障数据是随机截尾的，第ｉ台的故障观测时间为［0，T ｉ］, 其中T i 为现场试验的截尾时间。t 0≡0，因此, 得到相应的似然函数为

n j K

T T β1t ij β1-1β2t ij β2-1

L K =∏{∏[() +() ]⨯exp[-(i ) β1-(i ) β2]} (8)

i =1

j =1

α1α1α2α2α1α2

似然函数的对数函数以及此对数函数对模型参数的偏导数为

n j K

T T β1t ij β1-1β2t ij β2-1

l K =∑{∑ln[() +() ]-[(i ) β1+(i ) β2]} (9)

i =1

j =1

α1α1

-(

α2α2

m

α1α2

m

βm 2t ij

) () β1-1

K n αα∂l K

=∑∑

∂αm i =1j =1β1t ij β-1β2t ij β

()() +()() α1α1α2α2

j

m 1

+

2-1

βm T i β

() (10) αm αj

由（10）可以得到

[1+βm ln(

n j

t ij

K

T T ∂l K

=∑∑-(i ) βm ln(i ) (11) ∂αm i =1j =1β1ij β1-1β2ij β2-1αm αm

()() +()()

αm αm

)](

t ij

αm

) βm -1

α1α1α2α2

m -1

∑T

i =1

K

K

βm

=∑∑

i =1j =1

K

n j

ββm t ij

(

β1β

)() β-1+(2)() βα1α1α2α2

1

ij ij

(12)

2-1

由累积故障函数可得

K K K

T i T i

∑W (T i ) =∑[() β1+(() β2)]=∑n i (13)

i =1

i =1

α1α2

i =1

以上式可导出：

随机过程课程设计

K

α2=∑

i =1

T i

[n i -(

T i

1

(14)

α1

)]β1

将式（9）转换为三参数的函数，即

n j K

βt βt

'=∑{∑ln[(1)(ij ) β1-1+(2)(ij ) β2-1]-n i (15) l K

i =1

j =1

α1α1α2α2

最终，似然函数的参数估计转化成以下的求最大化问题：

' max l K 约束条件：

一般情况，最大化问题都需要初始值。根据经验，在没有合适的初始值选择下，可以假设

2.3 非齐次泊松过程下的可靠性指标

（１）首次故障间隔时间的可靠度函数【10-11】从t ＝0开始直到第一个故障发生的时间T １,T 1的可靠度函数为

R 1(t ) =Pr(T i >t ) =e

-W (t )

=e ⎰0

-ω(u ) du

t

(16)

对于非齐次泊松过程模型的使用，如果能够估算出首次故障间隔时间的故障率函数，就能同时估计出产品整个寿命的强度函数。

（２）其他故障间隔时间的可靠度函数

在t 0时刻后的可靠度函数

R (t |t 0) =Pr(t +t 0>t 0) =e

-

⎰t 0

t 0+t

ω(u ) du

(17)

（３）平均故障间隔时间

瞬时平均故障间隔时间———故障强度函数ω（t ）表示单位时间发生故障的次数，则ω(t)的倒数表示一次故障所经过的时间，定义瞬时故障平均间隔时间为

T IMTBF =

1

(18) ω(t )

累积平均故障时间间隔———表示一段时间内的平均故障间隔时间，即累积故障强度函数的倒数，

t -t t -t

T CMTBF =21=t 221 (19)

W (t 1, t 2) ω(u ) du

⎰

t 1

３实例分析

以国内北京第一机床厂同一时期出厂的６台加工中心现场试验的故障数据

为例，其发生故障的时间如表2所示。首先，需要对这些数据进行趋势检验。根据2.2节中多样本的趋势检验方法，在显著性水平α=0.05，得到这批加工中心的统计量值如表3所示。

编号

1 2 3 4 5 6

故障时间/h

50.99 423.72 753.06 760.65 795.63 1005.80 1209.40 2509.2 3350.16 3801.90 3915.62 4011.10 5109.03 5197.12 5353.92 5845.90 5942.81 6106.49 6323.63 6474.60 6526.03 6827.10 7059.69 7460.86 8240.27 8745.00 9142.65

185.13 458.00 960.54 1005.87 3409.55 422.89 6061.44 6217.53 7479.45 7542.81 7775.96 7882.88 7994.25 8588.25 28.05 350.48 47.52 1560.23 1896.30 2541.10 3352.80 3915.12 4981.45 5112.97 5729.13 5812.46 5903.40 6109.13 6117.21 6275.28 6308.78 6348.21 6457.61 6620.46 6853.44 7005.85 7116.59 7249.74 7467.90 8088.96 8298.18 9509.28 131.09 785.61 287.51 870.56 2987.45 3500.75 4881.86 5136.51 5230.01 5376.53 5540.54 5746.57 6183.21 6505.13 6592.08 7125.03 7379.46 7703.03 7868.85 8275.74 8654.42 9032.10

148.17 578.80 1014.14 1952.18 2893.01 3287.36 3747.55 4279.01 4714.12 4839.79 5558.09 5600.66 6694.61 6855.49 7120.48 7368.47 7496.84 7659.20 8451.26 8638.805

51.98 359.4 956.72 1357.45 1549.56 2706.15 3417.46 4659.60 5150.64 5206.74 5483.61 5570.40 5651.25

表２ 加工中心故障数据

在V 检验中运行结果拒绝H 0可知故障数据具有非单调趋势，且由表１可知，故障发生过程呈浴盆曲线的趋势。

非齐次泊松过程是随机点过程的一种典型类型，当可修系统的相邻故障间隔呈现某种趋势时可以使用这种方法来描述。因此，建立扩展的非齐次泊松过程模型，并且对加工中心的故障强度函数参数进行估计，通过matlab 运行结果得到尺度参数及形状参数，代入公式得到：

ω(t ) =

0.62t 1.97t

() 0.62-1+() 1.97-1 239.59239.592641.902641.90

本文在故障数据基础上，根据典型威布尔过程对数控机床进行了参数估计，并与文中模型得到的结果进行了对比。由图1中可以明显看出，本文采用的扩展的非齐次泊松过程能够更恰当地表现故障的趋势。

拟合结果如图１所示，本文方法的拟合值与估计值相关系数

由于样本量较大，且相关系数接近1，所以拟合值与估计值之间线性相关，故障过程符合假设的强度函数为浴盆曲线的威布尔过程。其相关的可靠性指标如下。

（１）首次故障间隔时间的可靠度函数

图2表示加工中心从t=0时刻投入运行后，继续无故障工作的可靠度曲线。

（２）平均故障间隔时间

图3是瞬时平均故障间隔时间和从t=0时刻起的累计平均故障间隔时间的曲线。从图中可以看出，在观测时间的初始阶段，平均故障间隔时间较高，随着

观测时间的增大，平均故障间隔时间变小，对照原始故障数据，可以看出平均故障间隔时间的趋势与现实情况相符。

４结束语

基于试验总时间法的数控机床非齐次泊松过程可靠性建模方法，不仅可以解决随机截尾故障数据趋势检验的问题，同时“修复如旧”的前提假设相对于以往的可靠性方法更加适用于数控机床的维修。6台加工中心实例的研究结果表明，本文所建立的模型更贴近于数控机床的实际运行状态和可靠性水平。

5、程序及结果：

%程序1：V 检验求出该过程具有非单调趋势

X=[50.99 423.72 753.06 760.65 795.63 1005.80 1209.40 2509.2 3350.16 3801.90 3915.62 4011.10 5109.03 5197.12 5353.92 5845.90 5942.81 6106.49 6323.63 6474.60 6526.03 6827.10 7059.69 7460.86 8240.27 8745.00 9142.65 185.13 458.00 960.54 1005.87 3409.55 422.89 6061.44 6217.53 7479.45 7542.81 7775.96 7882.88 7994.25 8588.25 28.05 350.48 47.52 1560.23 1896.30 2541.10 3352.80 3915.12 4981.45 5112.97 5729.13 5812.46 5903.40 6109.13 6117.21 6275.28 6308.78 6348.21 6457.61 6620.46 6853.44 7005.85 7116.59 7249.74 7467.90 8088.96 8298.18 9509.28 131.09 785.61 287.51 870.56 2987.45 3500.75 4881.86 5136.51 5230.01 5376.53 5540.54 5746.57 6183.21 6505.13 6592.08 7125.03 7379.46 7703.03 7868.85 8275.74 8654.42 9032.10 148.17 578.80 1014.14 1952.18 2893.01 3287.36 3747.55 4279.01 4714.12 4839.79 5558.09 5600.66 6694.61 6855.49 7120.48 7368.47 7496.84 7659.20 8451.26 8638.805 51.98 359.4 956.72 1357.45 1549.56 2706.15 3417.46 4659.60 5150.64 5206.74 5483.61 5570.40 5651.25];

disp('X');

A=[1 2;1/2 1];

[n,n]=size(A);

x=ones(n,100);

y=ones(n,100);

m=zeros(1,100);

m(1)=max(x(:,1));

x(:,2)=A*y(:,1);

m(2)=max(x(:,2));

y(:,2)=x(:,2)/m(2);

p=0.0001;i=2;k=abs(m(2)-m(1));

while k>p

i=i+1;

x(:,i)=A*y(:,i-1);

m(i)=max(x(:,i));

y(:,i)=x(:,i)/m(i);

k=abs(m(i)-m(i-1));

end

a=sum(y(:,i));

w=y(:,i)/a;t=m(i);

disp(w);

disp(t); %以下是V 检验CI=(t-n)/(n-1);

CI=(t-n)/(n-1);

RI=[0 0 0.52 0.89 1.12 1.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];

CR=CI/RI(n);

if CR

disp('接受H0!');

disp('CI=');

disp(CI);

disp('CR=');

disp(CR);

else

disp('拒绝H0!');

end

拒绝H0

%程序2:绘制累计故障曲线

t=100:100:10000;

w=(t./239.59).^(0.62)+(t/2641.9).^(1.97);

plot(t,w,'-*')

%程序3：绘制加工中心的可靠度曲线

t=100:100:10000;

w=(0.62/239.59)*(t./239.59).^(-0.38)+(1.97/2641.90)*(t/2641.9).^(0.97);

-3

%程序4：求出尺度参数及形状参数

A %6台机器故障数据

size(A)

B=sort(A);

K=1:n; % n-故障数据样本容量

%绘制累计故障曲线

N=K./7;

Y=log(N);

X=log(B);

plot(X,Y,’+’);

% ANN输入向量及目标向量

P=X;

T=Y;

pause

clc

plot(P,T,’+’);

title(‘Training Vectors’);

xlabel(‘Input vector P’);

ylabel(‘Target Vector T’);

pause

clc

[R,Q]=size(P);

[S,Q]=size(T);

% 初始化权向量

[W,B]=solvelin(P,T);

lr=2*maxlinlr(P);

df=10;

me=100;

eg=0.2;

tP=[df me eg lr];

[W,B,epochs,error]=trainwh(W,B,P,T,TP);

pause

W %威布尔过程形状参数B

exp(B) %威布尔过程强度参数

W=239.59 2414.9

B=0.62 1.97

6 参考文献：

[1]易蓓玲．可靠性与维修性工程概论［Ｍ］．北京：清华大学出版社.2010.

[2]张波．应用随机过程．中国人民大学出版社．[书 号] 7300037690． 2005 年8月．

[3]Yazhou J,Molin W,Zhixin J.Probability distribution of machining center failures[J].Reliability Engineering and System Safety.1995.

[4] Emerging Tecnologyies and Factory Automation Proceedings.1995.

[5]王智明，杨建国，王国强. 多台数控机床最小维修的可靠性评估. 哈尔滨工业大学学报. 2011

[6]面向不完全维修的数控机床可靠性评估，张根保，李冬英，刘杰《机械工程学报》2014.

[7] 胡仁胜. 实时控制软件系统的可靠性分析. 华南师范大学学报.2001.

[8]张英芝，申桂香，薛玉霞. 随机截尾数控机床故障过程［Ｊ］．吉林大学学报：工学版.2007.

评 阅 书