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

随机过程课程设计课程名称： 随机过程 课程设计（论文）题 目: 非齐次泊松过程 在数控机床可靠 性建模中的应用 学 院： 理学院 专 业： 数学与应用数学 班 级： 数学12-1班 学 生 姓 名： 王玲玲 学 生 学 号： 2012027149 指 导 教 师： 蔡吉花 2015 年 1月 3 日目 录任务书.1摘要.1前言.2 非齐次泊松过程理论 21.1 非齐次泊松过程的基本理论简介21.2 基于试验总时间法的趋势检验 22 数控机床的非齐次泊松过程可靠性建模 32.1强度函数的建立 . . 32.2 台数控机床强度函数的参数估计 . 42.3 非齐次泊松过程下的可靠性指标 . . 53实例分析54结束语75程序及结果86参考文献 9附录 评阅书 随机过程课程设计 随机过程 课程设计任务书姓名王玲玲学号19指导教师蔡吉花设计题目非齐次泊松过程在数控机床可靠性建模中的应用理论要点使用极大似然估计法对非齐次泊松过程的强度函数进行参数估计并且得到该模型的可靠性指标设计目标结合数控机床的维修特点，使用非齐次泊松过程建立可靠性模型贴近于复杂系统的生产实际，同时，结合具有随机截尾特点的多样本数控机床现场试验故障数据对数控机床的可靠性进行深入分析。研究方法步骤1 用观察，调查，统计，抽样等方法对统计的数据进行趋势检测；2 根据动态数据作出相关图形，进行相关分析；3. 然后做出模型，把平均故障间隔时间的趋势与现实情况相比较。预期结果利用非齐次泊松建立的模型更贴近于数控机床的实际运行状态和可靠性水平。计划与进步的安排第一步（1-2天）分析题目，查找资料；第二步（3-4天）针对性学习相关知识，整理思路；第三步（5-6天）编写程序；第四步（7天）用程序计算，写出结论。参考资料1易蓓玲可靠性与维修性工程概论北京：清华大学出版社.2010.2张波应用随机过程中国人民大学出版社书 号 7300037690 2005 年8月3Yazhou J,Molin W,Zhixin J.Probability distribution of machining center failuresJ.Reliability Engineering and System Safety.1995. 5王智明，杨建国，王国强. 多台数控机床最小维修的可靠性评估.哈尔滨工业大学学报. 2011填写时间2013.1.3 摘 要 基于试验总时间法对多样本随机截尾的数控机床现场数据进行趋势检验，在故障过程为浴盆曲线的趋势条件下构建了数控机床的非齐次泊松过程的可靠性模型。本文使用极大似然估计法对非齐次泊松过程的强度函数进行参数估计得到了该模型的可靠性指标，以6台加工中心的现场数据为例建立了非齐次泊松过程的可靠性模型。再通过matlab曲线拟合，绘制出故障时间的曲线，通过曲线的拟合程度，可以确定非齐次泊松过程能够更恰当地表现故障的趋势。关键词：数控机床 可靠性 非齐次泊松过程 浴盆曲线前言数控机床是由数目众多的零部件组成的复杂机电液可修系统。在其可靠性研究中，需要考虑维修活动对其可靠度的影响【1】。以往的数控机床可靠性建模方法，是将故障间隔时间视为独立同分布来分析其寿命分布，即假设维修活动是“修复如新”【2】而在实际生产中数控机床的维修活动是以调整或者更换一部分零部件和元器件为主的，对于复杂的系统来说这种维修活动只能使产品恢复到正常功能维修前后可靠度并没有很大改变，因此将数控机床的维修活动视为“修复如旧”更加合理。非齐次泊松过程经常被用于建立“修复如旧”的维修策略且维修时间可忽略的可修系统可靠性模型用于模拟出现故障间隔时间的趋势【3-4】。结合数控机床的维修特点，使用非齐次泊松过程建立的可靠性模型更能贴近于复杂系统的生产实际。本文提出了非齐次泊松过程的数控机床可靠性建模方法，并结合数控机床的失效特点，建立故障率为浴盆曲线的非齐次泊松过程可靠性模型。同时，结合具有随机截尾特点的多样本数控机床现场试验故障数据，对数控机床的可靠性进行了深入分析。非齐次泊松过程理论1.1 非齐次泊松过程的基本理论简介：非齐次泊松过程是随机点过程的一种典型类型，当可修系统的相邻故障间隔呈现某种趋势时可以使用这种方法来描述。非齐次泊松过程的重要参数【4】：为强度函数，是非负函数；其 累 积 故 障 强 度 函 数，表示在0,t中的平均故障数,即EN(t)=W（t），Nt 代表在0,t出现的故障次数表示机床从观测开始后的运行时间。当强度函数为时，成为威布尔过程。其中，、0，为尺度参数，为形状参数。01, 表示不断恶化的(坏)系统；=1，表示系统服从指数分布。1.基于试验总时间法的趋势检验：本文采用基于试验总时间的方法，对具有多样本随机截尾现场数据的故障过程进行趋势检验【5】。将在观测期间采集到的所有故障数据按照从大到小时间进行排序，得到t(i)的时间序列。根据试验总时间的建模思想【6】，得到该序列的第个故障发生时的试验总时间： （1） （2）式中：n(u)表示在时刻观察到的数控机床数量， 表示在观测期间的故障数；表示第台机床的故障数（共有台机床），当时间序列的最后一个时间是故障数据时，=-1；当时间序列的最后一个时间不是故障数据而是截尾数据时=。在实际检验时同时使用检验、检验和检验等检验方法综合确定有无趋势【7-8】。其中检验如下：0：齐次泊松过程；1：具有非单调趋势； (3) (4) (5)式中：（）表示总的观测时间。当，（）时，接受。一般情况下，检验和检验是检验具有单调趋势或齐次泊松过程和更新过程的故障数据的，而当故障数据具有非单调趋势时，则可以考虑检验中的统计量，如表所示。表1 故障率和故障强度函数变化特性 数控机床的非齐次泊松过程可靠性建模2.1强度函数的建立对于浴盆曲线趋势的故障过程，假设其故障强度函数由早期故障期和偶然故障期两部分组成，并且每一个阶段都是一个威布尔过程，参数为尺度参数和形状参数（，）。结合以上假设和多重威布尔分布模型的性质，则该数控机床故障强度函数为： （6）式中：，0在（0，内的平均故障个数为累积故障强度函数，即由于该模型是由两重威布尔过程构成，其强度函数是具有非单调的浴盆曲线趋势，因此组成该模型的两个形状参数有（）（）0，则本文中假设，。2.2 台数控机床强度函数的参数估计本文使用极大似然估计法对k台样本的强度函数进行参数估计，k台数控机床的故障数据是随机截尾的，第台的故障观测时间为0，T,其中Ti为现场试验的截尾时间。t00，因此,得到相应的似然函数为 (8)似然函数的对数函数以及此对数函数对模型参数的偏导数为 (9) (10)由（10）可以得到 (11) (12)由累积故障函数可得 (13)以上式可导出： (14)将式（9）转换为三参数的函数，即 (15)最终，似然函数的参数估计转化成以下的求最大化问题：约束条件：一般情况，最大化问题都需要初始值。根据经验，在没有合适的初始值选择下，可以假设2.3 非齐次泊松过程下的可靠性指标（）首次故障间隔时间的可靠度函数【10-11】从t0开始直到第一个故障发生的时间T,T1的可靠度函数为 (16)对于非齐次泊松过程模型的使用，如果能够估算出首次故障间隔时间的故障率函数，就能同时估计出产品整个寿命的强度函数。（）其他故障间隔时间的可靠度函数在t0时刻后的可靠度函数 (17)（）平均故障间隔时间瞬时平均故障间隔时间故障强度函数（t）表示单位时间发生故障的次数，则(t)的倒数表示一次故障所经过的时间，定义瞬时故障平均间隔时间为 (18)累积平均故障时间间隔表示一段时间内的平均故障间隔时间，即累积故障强度函数的倒数， (19)实例分析以国内北京第一机床厂同一时期出厂的台加工中心现场试验的故障数据为例，其发生故障的时间如表2所示。首先，需要对这些数据进行趋势检验。根据2.2节中多样本的趋势检验方法，在显著性水平=0.05，得到这批加工中心的统计量值如表3所示。编号故障时间/h12345650.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.2528.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.28131.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检验中运行结果拒绝H0可知故障数据具有非单调趋势，且由表可知，故障发生过程呈浴盆曲线的趋势。非齐次泊松过程是随机点过程的一种典型类型，当可修系统的相邻故障间隔呈现某种趋势时可以使用这种方法来描述。因此，建立扩展的非齐次泊松过程模型，并且对加工中心的故障强度函数参数进行估计，通过matlab运行结果得到尺度参数及形状参数，代入公式得到：本文在故障数据基础上，根据典型威布尔过程对数控机床进行了参数估计，并与文中模型得到的结果进行了对比。由图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);y(:,1)=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 kp 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);enda=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 CR0.10 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);plot(t,w,r-*)%程序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;pauseclcplot(P,T,+);title(Training Vectors);xlabel(Input vector P);ylabel(Target Ve