故障诊断学与可靠性工程作业（吐血分享） - 图文

机械故障诊断学及可靠性

工程

MECHANICAL FAULT DIAGNOSTICS AND

RELIABILITY ENGINEERING

作业题目 可维护系统有效度仿真

作者姓名 黎 原 作者学号 S12080203030 学科专业 机械设计及理论 指导教师 吴月明教授

2013年6月

机械故障诊断学及可靠性工程

可维护系统有效度仿真

硕士研究生： 黎 原 硕士生学号： S12080202030 导

师： 吴月明教授

学科专业： 机械设计及理论 所 在 单 位： 机械工程学院

Mechanical fault diagnostics and Reliability Engineering

THE SIMULATION OF VALIDITY FOR

MAINTAINABLE SYSTEM

by Li Yuan

Supervisor: Professor Wu Yueming

Yanshan University

May, 2013

摘 要 摘 要

混联可维护性系统的有效性仿真是一种基于计算机计算，借助于一定的编程软件而进行的离线模拟仿真技术。这一仿真技术，对现实的实验研究以及产品的实际生产有着重要的预测和指导作用。

本文主要针对由三类零件所组成的混联系统进行有效度的模拟仿真，当然，数学模型也做了相应简化，忽略了零件维修时间。通过0-1分布函数产生的伪随机数，运用威布尔分布函数得出相应的产品使用寿命，由此，对故障零件的故障时间进行模拟仿真，由于计算机的计算效率很快，所以模拟的零件样本空间可以取得足够大，单个零件的采样截止时间也可以取得很大。根据蒙特卡罗算法，通过对每一类零件的样本空间进行取样，并统计一定时间间隔的故障个数，最终运用MATLAB进行编程进而模拟出整个混联系统的有效度。

关键词：混联系统；威布尔分布；蒙特卡罗算法；有效度仿真；

I

可维护系统有效度仿真 Abstract

The simulation of validity Maintainable hybrid system is a offline simulation technology, which is base on computer calculation by using certain programming software. This simulation technology plays an important role in the the prediction and guidance for real experimental study and the actual production of products.

This paper mainly gives the simulation of validity for the system which is composed by three kinds of parts. By the way, the corresponding mathematical model has been simplified, ignoring the repairing times of the parts. a serial of pseudo-random numbers are generated by 0-1 distribution function. By using the Weibull distribution function, corresponding product using lives are obtained, whereby we can simulate the fault time of the part. Due you the perfect computational efficiency, the simulated parts sample space can get as big as possible, at the same time, the sample deadline for individual parts can also get very large. According to Monte Carlo algorithm, by getting the number of parts in each type of sample space and the number of fault parts in certain time interval t and then eventually programming by using MATLAB, we can finally simulate the validity of the entire hybrid system.

Keywords: Hybrid systems; Weibull distribution; Monte Carlo algorithm; Simulation of validity;

Ⅱ

目 录 目 录

摘 要.............................................................................................................................. I ABSTRACT .................................................................................................................. II 第一章 可维护系统的组成.......................................................................................... 1 1.1 串联系统 ............................................................................................................. 1 1.2 并联系统 ............................................................................................................. 1 1.3 混联系统 ............................................................................................................. 2 第二章 随机数的产生.................................................................................................. 4 2.1 随机数的定义及性质 ......................................................................................... 4 2.2 随机数产生的方法 ............................................................................................. 4 2.2.1 0-1分布产生随机数.................................................................................. 4 2.2.2 随机数表产生随机数................................................................................. 5 2.2.3 物理方法产生随机数................................................................................. 5 2.3 伪随机数性质 ..................................................................................................... 5 第三章 威布尔分布应用.............................................................................................. 7 3.1 威布尔分布简介 ................................................................................................. 7 3.2 威布尔分布的用途 ............................................................................................. 7 3.3 威布尔分析方法 ................................................................................................. 8 3.4 分布类型的选择 ............................................................................................... 12 第四章 蒙特卡罗算法................................................................................................ 15 4.1蒙特卡罗算法简介 ............................................................................................ 15 4.2 算法的基本思想及方法 ................................................................................... 15 4.3 算法的收敛性及误差 ....................................................................................... 16 4.4 方法的特点及适用范围 ................................................................................... 18 4.4.1该方法的优点和缺点................................................................................ 18 4.4.2 方法主要的应用范围............................................................................... 20 4.5 MATLAB软件简介 ........................................................................................... 20 4.5.1 MATLAB产生的历史背景.......................................................................... 20

III

可维护系统有效度仿真 4.5.2 MATLAB的语言特点.................................................................................. 21 4.5.3 MATLAB的优势特点.................................................................................. 22 第五章 混联系统有效度仿真.................................................................................... 25 5.1 系统有效度仿真的总体思想 ........................................................................... 25 5.2 用蒙特卡罗模拟法求有效度的步骤 ............................................................... 25 5.3 MATLAB仿真的源程序 ................................................................................... 26 5.4 系统有效度仿真的MATLAB图形输出 ......................................................... 29 5.4.1 程序代码截图........................................................................................... 29 5.4.2 程序结果输出截图................................................................................... 30 致 谢............................................................................................................................ 32 参考文献...................................................................................................................... 33

Ⅳ

第一章 可维护系统的组成 第一章 可维护系统的组成

1.1 串联系统

假设一个系统由n 个子系统组成，当且仅当所有的子系统都能正常工作时，系统才能正常工作，这种系统称为串联系统。组成系统的所有单元中任一单元的故障就会导致整个系统故障的系统称串联系统。它属于非贮备可靠性模型，其逻辑框图如图所示。

1

2

3

……

涡轮 n

压气机 燃烧室 尾喷管

如果系统各个子系统的可靠性分别用R1，R2，…，Rn表示，则系统的可靠性：R=R1?R2?????Rn

如果系统的各个子系统的有效度分别用A1，A2，…，An来表示，则系统的有效度：A?A1?A2?????An。

1.2 并联系统

假如一个系统由n个子系统组成，只要有一个子系统能够正常工作，系统就能正常工作。组成系统的所有单元都故障时，系统才故障的系统叫并联系统，它属于工作贮备模型。其逻辑框图如图所示。

1 2 ????n

如果系统的各个子系统的可靠性分别用R1,R2,…,Rn来表示，则系统的可靠性：R=1?(1?R2)?(1?R2)?????(1?Rn)。

如果系统的各个子系统的有效度分别用A1,A2,…,An来表示，则系统的有效

1

可维护系统有效度仿真 度：A?1?(1?A1)?(1?A2)?????(1?An)。

A B C D 在并联系统中只有一个子系统是真正需要的，其余n-1个子系统都被称为冗余子系统。该系统随着冗余子系统数量的增加，其平均无故障时间也会增加。

串联就是一个有问题就会瘫痪，并联只要有一个能用就没有问题。

1.3 混联系统

如果一个系统既含有串联构成的子系统又含有并联构成的子系统，且两个子系统之间或者是串联关系或者是并联关系，则该系统就是一个混联的系统。可分为串并联系统和并串联系统。

串并联系统：总体串联，局部有并联，如图所示：

1 串并联系统

3 2 等效系统

3 12

如果系统的各个子系统的有效度分别用A1,A2,A3来表示，则统的有效度：

A??1?(1?A1)?(1?A2)??A3。

并串联系统：总体并联，局部有串联，如图所示：

2

第一章 可维护系统的组成 1 并串联系统

2 3 4 12 等效系统

34

如果系统的各个子系统的有效度分别用A1,A2,A3,A4来表示，则系统的有效度：A?1?(1?A1?A2)?(1?A3?A4)。

3

可维护系统有效度仿真 第二章 随机数的产生

2.1 随机数的定义及性质

由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的关系，属于一般和个别的关系，或者说属于共性和个性的关系。由具有已知分布的总体中产生简单子样，就是由简单子样中若干个性近似地反映总体的共性。

随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方法可以由它产生具有任意已知分布的简单子样。

由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号ξ表示。由随机数序列的定义可知，ξ1，ξ2，…是相互独立且具有相同单位均匀分布的随机数序列。也就是说，独立性、均匀性是随机数必备的两个特点。

随机数具有非常重要的性质：对于任意自然数s，由s个随机数组成的s维空间上的点（ξn+1，ξn+2，…ξn+s）在s维空间的单位立方体Gs上均匀分布，即对任意的ai，当0?ai?1，i?1,2,?,s时有如下等式成立：

P(?n?i?ai,i?1,?,s)??ai

s其中P(·)表示事件·发生的概率。反之，如果随机变量序列ξ1, ξ2…对于任意自然数s，由s个元素所组成的s维空间上的点（ξn+1，…ξn+s）在Gs上均匀分布，则它们是随机数序列。

由于随机数在蒙特卡罗方法中所处的特殊地位，它们虽然也属于由具有已知分布的总体中产生简单子样的问题，但就产生方法而言，却有着本质上的差别。

i?12.2 随机数产生的方法

2.2.1 0-1分布产生随机数

在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称，随机数序列，其中每一个体称为随机数。

单位均匀分布也称为[0,1]上的均匀分布，其分布密度函数为：

?1,0?x?1f(x)??

0,其他?其分布函数为 ：

4

第二章 随机数的产生 ?0,?F(x)??x,?1,?x?00?x?1x?1

Matlab中，生成（0,1）区间上均匀分布的随机变量。基本语法：rand([M,N,P ...])，其表示：生成排列成M*N*P... 多维向量的随机数。如果只写M，则生成M*M矩阵。

2.2.2 随机数表产生随机数

为了产生随机数，可以使用随机数表。随机数表是由0,1,…,9 十个数字组成，每个数字以0.1的等概率出现，数字之间相互独立。这些数字序列叫做随机数字序列。如果要得到n位有效数字的随机数，只需将表中每n个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为7634258910…，要想得到三位有效数字的随机数依次为0.763,0.425,0.891。

因为随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。

2.2.3 物理方法产生随机数

用物理方法产生随机数的基本原理是：利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。用来作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。

一般情况下，任意一个随机数在计算机内总是用二进制的数表示的：

???1?2?1??2?2?2????m?2?m

其中εi（i=1,2,…,m）或者为0，或者为1。

因此，利用物理方法在计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5。

用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适合在计算机上使用。

2.3 伪随机数性质

判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类：均匀性检验和独立性检验。

六十年代初，人们开始用定性的方法研究伪随机数序列的均匀性和独立性问题，简要叙述如下。

5

可维护系统有效度仿真 （1） 伪随机数的均匀性

这里只考虑伪随机数序列ξ1,ξ2…,ξn全体作为子样时的均匀性问题。其中n为伪随机数序列的最大容量。

对于任意的0≤x≤1，令Nn(x)表示伪随机数序列ξ1,ξ2…,ξn中适合不等式ξi< x，（其中 i=1,2,…,n）的个数，则：

?(n)?sup　|n0?x?1N(x)n?x|

标志伪随机数序列ξ1,ξ2…,ξn的均匀程度，称为均匀偏度。

?????n?，并令将伪随机数序列ξ1,ξ2…,ξn从小至大重新排列?1???2??1?1，则由δ(n)的定义，容易证明： ?0??0,?n?(n)?max?|?i??0?x?1??in|,|?i??1?i?|? n?很明显，对于固定的ｎ,δ(n)的值越小越好。它是描述伪随机数序列均匀程度的基本量。对于任意随机数序列，均有如下不等式成立：?(n)?时，所对应的伪随机数序列为最佳分布。

可以证明，伪随机数序列为最佳分布的充要条件是它取遍序列

34n2i?12n1n12n，当 ?(n)?12n的所有

值。对于计算机上使用的乘同余方法，按照前面介绍的方法选取a、x1时，所产生的伪随机数序列的均匀偏度?(n)?。对于乘加同余方法?(n)?，对于部

分伪随机数的均匀性问题通常用统计检验方法检验。 （2）伪随机数的独立性

对于任意0?x,y?1，令Nn(x,y)表示(ξ1,ξ2),(ξ2,ξ3),…, (ξn,ξn+1)中适合不等式?i?x,?i?1?y的个数，根据随机变量间相互独立的定义和频率近似概率的方

|n法，令?(n)?sup　0?x,y?1N(x,y)n?Nn(x)Nn(y)nn|，则ε(n)标志伪随机数序列ξ1,ξ2…，

ξn的独立程度，简称为独立偏度。对于固定的n, ε(n)的值越接近于零，伪随机数序列的独立性越好。

对于乘同余方法， ?(n)?1n?3n，对于乘加同余方法，?(n)?1n?4n，因

此，这两种方法的独立性都是很好的。

同伪随机数的均匀性问题一样，伪随机数序列的独立性问题也是对它的全体讨论的。若只考虑伪随机数的一部分，在通常情况下给出ε(i)是相当因难的。因此，伪随机数序列的独立性问题的统计检验方法同样是非常重要的。

6

第三章 威布尔分布应用 第三章 威布尔分布应用

3.1 威布尔分布简介

在所有可用的可靠性计算的分布当中，威布尔分布是唯一可用于工程领域的。在1937，Waloddi Weibull教授（1887－1979）创造性的提出了该种分布，它是用于失效数据分析分布中应用最广泛的分布之一，也用于寿命数据分析，因为系统或部件的寿命周期的测量也需要分析。

一位瑞典的工程师和一位数学家潜心研究冶金的失效，威布尔教授曾指出正态分布要求冶金的初始强度服从正态分布，而情况并非如此。他还指出对于功能需求可以包含各种分布，其中包括正态分布。

1951年他发表了代表作，“一个具有广泛适用性的统计分布函数”，威布尔教授声称寿命数据可以从威布尔分布族中选择最恰当的分布，然后用合适的参数进行合理准确的失效分析。他列举七种不同的情况来证明威布尔分布可顺利用于很多问题的分析。

对威布尔分布的最初反应是普遍诊断它太过完美以致于不真实。尽管如此，失效数据分析领域的先驱们还是开始应用并不断改进，直到1975年，美国空军才认可了它的优点并资助了威布尔教授的研究。

今天，威布尔分析涉及图表形式的概率分析以找出对于一个给定失效模式下最能代表一批寿命数据的分布。尽管威布尔分布在检测寿命数据以确定最合适的分布方面在世界范围内处于领先位置，但其它分布也会偶尔用于寿命数据分析包括指数分布，对数正态分布，正态分布，寿命数据有了对应的统计学分布，威布尔分析对预计产品寿命做了准备。这种具代表性的样本分布用来估计产品的重要寿命特征，如可靠性，某一时刻的失效率，产品的平均寿命及失效率。

3.2 威布尔分布的用途

威布尔分析广泛用于研究机械、化工、电气、电子、材料的失效，甚至人体疫病。威布尔分析最主要的优点在于它的功能：

（1） 提供比较准确的失效分析和小数据样本的失效预测，对出现的问题尽早的

制订解决方案。

（2） 为单个失效模式提供简单而有用的图表，使数据在不充足时，仍易于理解。 （3） 描述分布状态的形状可很好的选择相应的分布。 （4） 提供基于威布尔概率图的斜率的物理失效的线索。

虽然对数或对数正态分布的使用通常要至少20次失效或源于以往的经验，在只有2~3次失效时用威布尔分析非常好，在涉及安全性或极端费用时的失效结

7

可维护系统有效度仿真 果是很关键的。威布尔家族中的一员weibayes，在以往经验充足时甚至可用于无失效情况下。

威布尔分析一般用于以下方面失效数据的分析： （1） 研制、生产和服务 （2） 质量控制和设计缺陷 （3） 维修计划和替代方案 （4） 备用元件的预测 （5） 保障性分析

（6） 自然灾害（闪电袭击，暴风雪，强风，暴雪等）

威布尔分析新的应用包括医学研究，仪器校准，费用削减，材料性能和测量分析。

3.3 威布尔分析方法

威布尔分析研究的是通过在威布尔概率图上绘制单一失效模式的寿命数据来研究部件的寿命时间和它的可靠度之间的关系。威布尔分析最常用于描述元器件失效的时间，它们可以是电灯泡，滚珠轴承、电容、磁盘驱动器，打印机甚至是人。失效模式包括爆裂，折断，变形或由于腐蚀造成的疲劳，过应力，高温，佐佐木希北川景子板野友美观月雏乃初期致命失效，耗损等等。

当在威布尔概率图上绘制失效时间数据时，工程师们更愿意用median rank regression作为参数估计方法，median rank regression方法是通过用最小二乘法（曲线拟合），找到一条最佳拟合直线来将平方差减至最小，median rank regression被认为是标准参数估计方法，因为它通过大多数数据得出了正确结果。

典型的，水平刻度（x轴）度量部件的寿命，垂直刻度（Y轴）度量已知失效模式下的部件失效累积的百分数。

一个威布尔概率图沿着横坐标有一条线性/非线性的时间刻度，沿着纵坐标有另一条非线性的分布函数。这些非线性的刻度通过适当的数据模型选出。如果刻度与数据相匹配，图表就会呈现出一条直线。由于它们简单且有用，所以概率图表用于统计分析中已经很多年了。尽管如此，仍需注意的是用概率描绘的方法获得的分布参数是独立同分布的，这经常用于不可修的部件和系统，而对于可修系统的失效数据可能就不是这样。

在图中，威布尔概率图认为失效时间对应唯一的失效模型。当许多元器件在正常运转条件下被测试时，它们不会在同一时间因同一原因都失效。任一失效原因下的失效次数都会集中于平均值附近，次数过多或过少的情况都较少。由于寿命数据的分布如此，他们会服从某种分布。为了描述一种分布的形状，这种分布的形状取决于所要研究的内容，公式可由统计方法得出。如果已绘制的数据点落

8

第三章 威布尔分布应用 在直线附近，威布尔概率图便认为是合理的。

注意：Y轴上的值是从1%~99%的概率值，轴上各点之间的距离是不均匀的。威布尔概率图的X、Y轴上的点于点之间的距离是百分比的变化而不是点的变化。正如对数的刻度一样，1~2间的距离是100%的增加，与2~4间的距离相同，但那是另一个100%的增加。对数比例只为一些相似级数作铺垫。除了对问题有更深的洞察力，最直观的是对确认分布方法有帮助，该种方法可更好的将数据集构成一条直线。

如果用以前的数据表示发生的失效，将组件的失效寿命绘制成图是非常常见的。在这种情况下：

??1??Y轴通常为：ln?ln??? 1?F(t)????X轴为：ln?t?

Y轴的截距为：??ln???

双参数的威布尔分布目前在寿命数据分析中广泛应用：

?????t??R(t)?exp?????

???????其中：

t≥0,β＞0 且η>0 。这里,β和η分别是分布状态和比例参数。

因为双参数的威布尔分布有效地分析了初期致命失效，实用寿命的和耗损阶段的寿命数据，它也可用于失效率的增长，持续和递减。

定义了威布尔概率图的第一个参数是斜率β，它是形状参数，因为它确定了威布尔家族中哪一种分布相关性最好或可以描述数据。第二个参数是特征寿命，

9

可维护系统有效度仿真 伊塔（η）作为比例参数，因为它定义了分布状态的大部分。参数β和η可从寿命数据中估计，寿命数据总为正值。威布尔分析完成后，由图可看出威布尔概率的斜度和拟合度。

注意：三参数的威布尔分布应用也很广泛。第三个参数——位置，是一个常数，可从时间变量t中加上或减去。

威布尔危险函数或失效率依赖于β的值，因为β值说明了新或旧元件是否更有失效的可能，威布尔危险函数可以描绘出不同元件的浴盆曲线：

初期故障：在电子和制造业中，早期失效指在使用寿命的初期失效的概率极高，当β值于小1.0时，威布尔概率分布图表明较新的元件在正常使用时更有可能失效，被称为瞬时递减失效率。为中止电子和机械系统在早期故障的高失效率，制造商提供了产品接收测试，老练（burn－in）早期和环境应力筛选暂不先将系统交付客户。假如有部件在初期损失阶段没有失效，那么它的失效率应当是递减的，且它的可靠度增加。因此，旧元件被认为比新元件更好，因为新元件很可能在寿命的早期失效，而元件在早期失效阶段的检修是不合适的。

偶然故障：假设威布尔概率分布图以一个独立失效模型为基础，β为1.0说明失效率是常数或相对于时间独立。这意味着对于那些无故障运行至时间t的元件，在下一个单位时间内将不能保持恒定的百分比，称作恒定危险率或瞬时失效率。这使得威布尔概率图与指数分布一致。由于旧元件被认为与新元件一样好。检修通常是不适合的，唯一使系统或部件可靠度提高的方法是用随机失效进行重新设计。

早期损耗：在设计寿命时经常因为机械问题出现未预期的失效。当1.0<β<4.0时，大修或以低B－lives来替换元件会较经济。B-Lives指出给定总数的百分比失效时的时间。例如：B-1寿命是指总数的1%失效时的时间，而B-10寿命指总数的10%失效的时间。通过优化预防性维修计划，经历早期耗损的元件的可靠度和费用都会提高。

快速损耗：尽管一个元件设计寿命的β值大于4.0是需要引起重视，但多数斜率急剧升降的威布尔概率图在失效概率在忽视范围内有一个安全期，且发生失效的影响会超出设计寿命。斜率越大的直线，在失效时间内的变化越小且结果越可预知。对于有重大失效的元件，大修和检查会更经济，因为定时维护会较昂贵，所以当旧元件快损坏或失效时才会考虑，此时的失效称为瞬时增长失效率。

因为不同的斜率代表不同的失效类别，威布尔分布提供了可能引起失效的原因，下表列出了引起每一类失效的失效原因：

β值 β<1.0 类型 初期故障 斜述 当β<1.0,失效原因归结于： 10

第三章 威布尔分布应用 ? ? ? ? ? β=1.0 随机失效 不充足的burn-in或应力筛选 部件的质量问题 制造的质量问题 错误的安装,设置及使用 重做/刷新时出现的问题 维护中的人为错误 引发的失效而非固有的 意外事故和自然灾害(外来物体，闪电袭击，强风摧毁等) 当β=1.0，失效原因归结于： ? ? ? 1.0<β<4.0 早期损耗 当1.0<β<4.0，失效原因归结于： ? ? ? ? β>4.0 快速损耗 低循环疲劳 受力失效 腐蚀/侵蚀 制造过程 当β>4.0，除部件老化，还有以下原因引起失效： ? ? ? 材料的固有属性的缺陷(如陶瓷易碎) 制造过程中出现的严重问题 制造或材料上的细微变化 统计学家，数学家和工程师们已将统计分布简化为数学模型或描绘出某些行为。与其它统计分布相比,威布尔分布适于更广范围的寿命数据。威布尔概率密度函数是一个数学函数，用以描述与数据相适应的曲线。概率密度函数可用数学模型给出或用图形给出，其中图上X轴代表时间。威布尔家族中的不同成员有不同形状的概率密度函数。累积密度函数是概率密度函数曲线下的面积。威布尔分布的累积密度函数如下公式：

????t???F(t)?1?exp?????

???????其中:

η代表特征寿命(比例参数) β代表斜率(状态参数)

累积密度函数给出了时间t内的失效概率.参数η和β由失效时间进行估计，如果失效数据来自于威布尔分布，η和β的值代入累积密度函数的公式求出一定时间内元器件的失效预计。

11

可维护系统有效度仿真 特征寿命η和平均失效时间（MTTF）是相关的。特征寿命给出了系统或元器件寿命中的失效概率独立于失效分布参数的点。对所有威布尔分布来说，定义为63.2%的单元失效时的寿命。

虽然，威布尔教授最初提出用平均值作为MTTF值绘制在威布尔概率分布图的y轴上，现在是标准的工程方法用失效时间的中间值来划分寿命数据。表7-2展示了一个中间等级表（50%）作为10个数量的样本，由此形成莱奥纳多·杰克逊（Leonard Johnson）的等级公式。

因为在寿命数据中非均匀分布相当常见，所以中间值比均值更为准确些。一旦知道β和η，任意时间的失效概率都可轻易算出。

3.4 分布类型的选择

威布尔分布有不同形式的应用，包括单参数，双参数，三参数及混合威布尔分布，有时不属于威布尔分布的正态或对数正态分布也可用于寿命数据的分析，选择适合于特殊数据集的分布要以数据的数量和质量，以往经验以及良好的相关性测试为基础。下表描述了威布尔家族的各种分布。 双参数威布尔分布 双参数威布尔分布的所需参数是斜率和特征寿命。这种分布用小样本提供了正确合理的失效分析和失效预计。它尤其可以诊断出失效类型，例如初期损耗（尤其是电子产品），独立时间失效（意外事故和固有事件的发生）或耗损的构件（轴承、过滤器等），如果失效率递减（老练时期burn-in period）或递增（耗损阶段），或失效率保持恒定（随机失效阶段），推荐使用双参数威布尔分析。 指数分布 指数分布唯一的参数是失效率。指数分布可被视为威布尔分布的一种特例——β＝1。当部件的失效率恒定时，它的可靠度最好用威布尔分布或指数分布来描述。失效率恒定会产生无记忆属性，即一个使用过的部件的寿命与当前老化时间无关，因此可以说一个使用过的部件像一个新部件一样好（只有当β为1时，威布尔分布才是无记忆的）。因为指数分布假定没有初期故障或耗损阶段，所以区域内数据要经仔细测试以确保那些假设正确。对于指数分布，MTTF与失效率互为倒数。 瑞利分布 瑞利分布的唯一参数是特征寿命。当β＝2时，瑞利分布可被视为威布尔分布的一种特例。然而，它有其自身优点，是一个重要的分布，它不仅应用于可靠性问题也用于与通讯系统相关的噪声问题。作为与指数分布相似的一个单参数分布，瑞利分布可被用来描绘错误源的均方根值。如果失效率随时间线性增加时可推荐12

第三章 威布尔分布应用 使用瑞利分布。 Weibayes分布 WeiBays分布的唯一参数是特性寿命，也叫单参数威布尔分布。WeiBayes是威布尔分布的一种特例——斜率β定义如前所述。与Bayesian假定有关，当用传统的威布尔分析产生许多不确定因素时， Weibays便是一个有效的解决问题的方法。当样本小于10个失效数时，Weiboys分布比双参数威布尔分布更精确，而且它是在失效数为零时唯一可用的分布。例如，在对一个现存失效模式进行设计修改后，从测试中获取的有效数据可用于确定新设计中威布尔曲线的置信下限被成为Weibayes曲线。当元器件超出其设计寿命时，无失效的威布尔分析可延长其寿命。由于Weibayes分布可用于无失效测试要求的情况下，所以它的重要性全在于失效对安全性和极限成本的影响。 三参数威布尔分布 除斜率或特征寿命参数外，三参楼威布尔分布还包括一个位置参数t0，它定义了分布在时间上的位置。这个参数可转换时间刻度的原点，而且只有被双参数分布分析证明为是合适之后，才可被使用。（有关其它信息可参考172页“威布尔概率分布图的曲线数据”）。使用位置参数时，在生成威布尔概率分布图之前，t0的值可以从时间值中减去或加上。例如，在某段时间内如果失效率为零，那么时间刻度的原点应从0转换到t0处以反映此阶段为无失效保证阶段，修正量t0为一个正值，等于失效发生的最小时间。由于正式使用前会出现寿命（可靠性）损失，t0为一个负值。负修正对于仓库中闲置元件的腐蚀的情况是有用的。例如，橡胶部件，化学品和滚珠轴承，都会随贮藏时间的延长而腐蚀老化。当t0值用于数据修正时，结果可绘成一条直线。在没有经验的情况下，用三参数威布尔分布进行分析时通常需要至少20个失效数据。 Gumbel 在19世纪20年代，E.J.Gumbel第一个认真调查失效数据的极值，找到只有6个独立极值的分布。他的第III类最小极值分布与威布尔分布相同。Gumbel-(Lower)distribution也叫第I类极小值分布，是一种极小值分布。Gumbel-(upper)distribution也叫第一类极大值分布，是一种极大值分布。当失效数据为偶然性事件的结果并且失效数据取极值时，推荐用Gumbel分布。示例如自然灾害和最大载客量等。因为Gumbel分布（和正态分布）可用于预计高可靠度要求下无寿命数据情况下的负寿命值，所以在建立寿命数据模型时要谨慎使用。 13

可维护系统有效度仿真 尽管统计学家反对用极小样本，但安全性和明显的资金损失却决定了我们收集的数据的局限性。当仅有极少失效数据存在时，威布尔分析可提供有用的结果，因为：

（1） 耗损失效发生在最陈旧的单元上。多数失效结果被绘在威布尔概率分布图

的左下角B-0.1至B-1，这也正是在工程上最为关注的区域。

（2） 威布尔分析包括失效和中止。尽管中止没有失效严重，但会有上千的中止，

在B-0.1~B-1lives进行更准确的工程预计。

威布尔分布应用于失效机会倍增且第一次失效很重要的情况，也应用于线性衰退而不是加速衰退系统。当威布尔分布是非线性衰退而不是当前衰退的一个函数时，可使用对数正态分布。下表给出了正态和对数正态分布的描述。因为即使它们不属于威布尔分布但偶尔也用来做寿命数据的参数分析。多数威布尔软件可快速生成所有分布并自动为数据集选出一个最合适的分布。 正态分布或 分布 正态分布的两个参数是均值和标准差。正态分布是对称的，一态分布经常用于描述失效率随时间增加的设备。当失效时间可用某些随机变量的总和表示时才推荐使用正态分布。正态分布便于描述不同类型数据，它允许观测结果为负。由于时间t大于零后元件才会失效，所以寿命数据总是正的。因此正态分布并不能很好的描述寿命数据。分析人员并不为正态配置所困扰，因为利用正态分布的寿命数据也可绘制的威布尔概率分布图。 对数正态分布 对数正态分布的两个参数是均值和标准差。尽管对数正态分布与正态分布相似。它假定了随机变量的对数值是正态分布而不是随机变量本身服从正态分布。因而，所有值均为正，分布图像便不会向左倾斜。对数正态分布很可能是威布尔分布强有力的竞争对手。它多用于工程上的金属疲劳测试、维修数据（修复时间）、化学反应过程的仪器失效和维修，一些材料特性和非线性加速衰退。若有倍增的失效因素影响失效时间，推荐使用对数正态分布。例如，在衰退渐增的情况下，断裂形式是由于压力造成，且压力随裂缝增大而增加。对数正态分布的非工程应用包括私人收入分析，遗产继承和银行抵押。 高斯Gaussion般被称为贝尔曲线，该分布很重要且广泛用于概率统计中。正14

第四章 蒙特卡罗算法 第四章 蒙特卡罗算法

4.1蒙特卡罗算法简介

Monte Carlo方法，又名随机模拟法或统计实验法它是以概率统计理论为基础，解决一些很难直接用数学运算求解或用其他方法不能解决的复杂问题的一种近似计算法。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，Monte Carlo方法是一种有效的求出数值解的方法。

20世纪40年代以后，随着电子计算机的出现和发展，人们有可能用计算机来模拟这类实验和计算。计算机具有计算速度高和存储容量大的特点，采用数字模拟技术可以代替许多实际上非常庞大而复杂的实验，并迅速将实验结果进行运算处理，于是Monte Carlo 方法重新被提起，引起世人重视，应用日渐广泛。实际上，采用Monte Carlo 方法在计算机上建立模型来解决Buffon 问题是非常简单的。

蒙特卡罗算法是以概率和统计理论方法为基础的一种计算方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。由S.M.乌拉姆和J.冯·诺伊曼在20世纪40年代为研制核武器而首先提出 。为了求解数学、物理、工程技术以及管理等方面的问题 ，首先建立一个概率模型或随机过程，使它们的参数，如概率分布或数学期望等问题的解；然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征，并用算术平均值作为所求解的近似值。对于随机性问题，有时还可以根据实际物理背景的概率法则，用电子计算机直接进行抽样试验，从而求得问题的解答。

4.2 算法的基本思想及方法

传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。这也是我们采用该方法的原因。蒙特卡罗方法的基本原理及思想如下：当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。这就是蒙特卡罗方法的基本思想。

蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。它是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。可以把蒙特卡罗解题

15

可维护系统有效度仿真 归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。

（1）构造或描述概率过程

对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。 （2）实现从已知概率分布抽样

构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。 建立各种估计量： 一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。

4.3 算法的收敛性及误差

1.收敛性

蒙特卡罗方法是由随机变量X的简单子样X1，X2，?，XN的算术平均值：

XN?1N?Xi?1Ni

作为所求解的近似值。由大数定律可知，如果X1，X2，?，XN独立同分布，且具有有限期望值（E(X)<∞），则有：

??P?limXN?E(X)??1 ?N???16

第四章 蒙特卡罗算法 即随机变量X的简单子样的算术平均值XN，当子样数Ｎ充分大时，以概率1收敛于它的期望值E(X)。 2.误差

蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，?，XN独立同分布，且具有有限非零的方差σ2 ，即

0??2??(x?E(X))2f(x)dx??

其中，f(X)是X的分布密度函数。则有：

?N?P?XN?E(X)?x??lim???N????12??x?xe?t/2dt

2当N充分大时，有如下的近似式

????P?XN?E(X)????N??22????0e?t/2dt?1??

2其中α称为置信度，1－α称为置信水平。 这表明，不等式XN?E(X)?度的阶为O(N?1/2)。

通常，蒙特卡罗方法的误差ε定义为?????N

???N近似地以概率1－α成立，且误差收敛速

上式中??与置信度α是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出??。

下面给出几个常用的α与的数值：

α 0.5 0.05 0.003 λα 0.6745 1.96 3 关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差σ是未知的，必

??须使用其估计值?X?Ni?11N2i?(X?Ni?11Ni)2来代替，在计算所求量的同时，可计?。 算出?3.减小方差的方法

17

可维护系统有效度仿真 显然，当给定置信度α后，误差ε由σ和N决定。要减小ε，或者是增大N，或者是减小方差σ2。在σ固定的情况下，要把精度提高一个数量级，试验次数

N需增加两个数量级。因此，单纯增大N不是一个有效的办法。

另一方面，如能减小估计的均方差σ，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。因此降低方差的各种技巧，引起了人们的普遍注意。后面课程将会介绍一些降低方差的技巧。

一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中效率的概念。

22它定义为??c，其中c是观察一个子样的平均费用。显然??c越小，方法越有效。

4.4 方法的特点及适用范围

4.4.1该方法的优点和缺点

（1）能够比较逼真地描述具有随机性质的事物的特点及物理实验过程

从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。 （2）受几何条件限制小

在计算

Dss维空间中的任一区域Ds上的积分

g????g(x1,x2,?,xs)dx1dx2?dxs时，无论区域Ds的形状多么特殊，只要能给

(i),?,xs(i))，出描述Ds的几何特征的条件，就可以从Ds中均匀产生N个点(x1(i),x2得到积分的近似值。

gN?DsN?g(xi?1N(i)1(i),x2,?,xs(i))

其中Ds为区域Ds的体积。这是数值方法难以作到的。

另外，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。 （3）收敛速度与问题的维数无关

由误差定义可知，在给定置信水平情况下，蒙特卡罗方法的收敛速度为

O(N?1/2)，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算

时间的变化，不影响误差。也就是说，使用蒙特卡罗方法时，抽取的子样总数N

18

第四章 蒙特卡罗算法 与维数s无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。

（4）具有同时计算多个方案与多个未知量的能力

对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。

另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。 （5）误差容易确定

对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。

一般计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。 （6）程序结构简单，易于实现

在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。 （7）收敛速度慢

如前所述，蒙特卡罗方法的收敛速度为O(N?1/2)，一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问题，不如其他方法好。 （8）误差具有概率性

由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。

（9）在粒子输运问题中，计算结果与系统大小有关

经验表明，只有当系统的大小与粒子的平均自由程可以相比较时（一般在十个平均自由程左右），蒙特卡罗方法计算的结果较为满意。但对于大系统或小概率事件的计算问题，计算结果往往比真值偏低。而对于大系统，数值方法则是适用的。

19

可维护系统有效度仿真 因此，在使用蒙特卡罗方法时，可以考虑把蒙特卡罗方法与解析（或数值）方法相结合，取长补短，既能解决解析（或数值）方法难以解决的问题，也可以解决单纯使用蒙特卡罗方法难以解决的问题。这样，可以发挥蒙特卡罗方法的特长，使其应用范围更加广泛。

4.4.2 方法主要的应用范围

蒙特卡罗方法所特有的优点，使得它的应用范围越来越广。它的主要应用范围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。随着科学技术的发展，其应用范围将更加广泛。

蒙特卡罗方法在粒子输运问题中的应用范围主要包括：实验核物理，反应堆物理，高能物理等方面。

蒙特卡罗方法有很强的适应性，问题的几何形状的复杂性对它的影响不大。该方法的收敛性是指概率意义下的收敛，因此问题维数的增加不会影响它的收敛速度，而且存贮单元也很省，这些是用该方法处理大型复杂问题时的优势。因此，随着电子计算机的发展和科学技术问题的日趋复杂，蒙特卡罗方法的应用也越来越广泛。它不仅较好地解决了多重积分计算、微分方程求解、积分方程求解、特征值计算和非线性方程组求解等高难度和复杂的数学计算问题，而且在统计物理、核物理、真空技术、系统科学 、信息科学 、公用事业、地质、医学，可靠性及计算机科学等广泛的领域都得到成功的应用。

蒙特卡罗方法在实验核物理中的应用范围主要包括：通量及反应率，中子探测效率，光子探测效率，光子能量沉积谱及响应函数，气体正比计数管反冲质子谱，多次散射与通量衰减修正等方面。

4.5 MATLAB软件简介

4.5.1 MATLAB产生的历史背景

在20世纪70年代中期，Cleve Moler博士和其同事在美国国家科学基金的资助下开发了调用EISPACK和LINPACK的FORTRAN子程序库。EISPACK是特征值求解的FORTRAN程序库，LINPACK是解线性方程的程序库。在当时，这两个程序库代表矩阵运算的最高水平。

到20世纪70年代后期，身为美国New Mexico大学计算机系系主任的Clev e Moler，在给学生讲授线性代数课程时，想教学生使用EISPACK和LINPACK程序库，但他发现学生用FORTRAN编写接口程序很费时间，于是他开始自己动手，利用业余时间为学生编写EISPACK和LINPACK的接口程序。Cleve Moler给这个接口程序取名为MATLAB，该名为矩阵（matrix）和实验室（laboratory）两个

20

第四章 蒙特卡罗算法 英文单词的前三个字母的组合。在以后的数年里，MATLAB在多所大学里作为教学辅助软件使用，并作为面向大众的免费软件广为流传。

1983年春天，Cleve Moler到Stanford大学讲学，MATLAB深深地吸引了工程师John Little。John Little敏锐地觉察到MATLAB在工程领域的广阔前景。同年，他和 Cleve Moler、Sieve Bangert一起，用C语言开发了第二代专业版。这一代的MATLAB语言同时具备了数值计算和数据图示化的功能。

1984年，Cleve Moler和 John Lithe成立了MathWorks公司，正式把MATLAB推向市场，并继续进行MATLAB的研究和开发。

在当今30多个数学类科技应用软件中，就软件数学处理的原始内核而言，可分为两大类。一类是数值计算型软件，如 MATLAB、Xmath、Gauss等，这类软件长于数值计算，对处理大批数据效率高；另一类是数学分析型软件，如Mathematica、Maple等，这类软件以符号计算见长，能给出解析解和任意精度解，其缺点是处理大量数据时效率较低。MathWorks公司顺应多功能需求之潮流，在其卓越数值计算和图示能力的基础上，又率先在专业水平上开拓了其符号计算、文字处理、可视化建模和实时控制能力，开发了适合多学科、多部门要求的新一代科技应用软件MATLAB。经过多年的国际竞争，MATLAB 已经占据了数值型软件市场的主导地位。

4.5.2 MATLAB的语言特点

一种语言之所以能如此迅速地普及，显示出如此旺盛的生命力，是由于它有着不同于其他语言的特点。正如同FORTRAN和C等高级语言使人们摆脱了需要直接对计算机硬件资源进行操作一样，被称作为第四代计算机语言的MATLAB，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来。MATLAB的最突出的特点就是简洁。MATLAB用更直观的、符合人们思维习惯的代码，代替了C和FORTRAN语言的冗长代码。MATLAB给用户带来的是最直观、最简洁的程序开发环境。以下简单介绍一下MATLAB的主要特点。

语言简洁紧凑，使用方便灵活，库函数极其丰富。MATLAB程序书写形式自由，利用其丰富的库函数避开繁杂的子程序编程任务，压缩了一切不必要的编程工作。由于库函数都由本领域的专家编写，用户不必担心函数的可靠性。可以说，用MATLAB进行科技开发是站在专家的肩膀上。

具有FORTRAN和C等高级计算机语言知识的读者可能已经注意到，如果用FORTRAN或C语言去编写程序，尤其当涉及矩阵运算和画图时，编程会很麻烦。例如，如果用户想求解一个线性代数方程，就得编写一个程序块读入数据，然后再使用一种求解线性方程的算法（例如追赶法）编写一个程序块来求解方程，最后再输出计算结果。在求解过程中，最麻烦的要算第二部分。解线性方程的麻烦

21

可维护系统有效度仿真 在于要对矩阵的元素作循环，选择稳定的算法以及代码的调试都不容易。即使有部分源代码，用户也会感到麻烦，且不能保证运算的稳定性。解线性方程的程序用FORTRAN和C这样的高级语言编写至少需要好几十行。再如用双步QR方法求解矩阵特征值，如果用FORTRAN编写，至少需要四百多行，调试这种几百行的计算程序可以说很困难。

由于MATLAB是用C语言编写的，MATLAB提供了和C语言几乎一样多的运算符，灵活使用MATLAB的运算符将使程序变得极为简短。MATLAB既具有结构化的控制语句（如for循环、while循环、break语句和if语句），又有面向对象编程的特性。语法限制不严格，程序设计自由度大。例如，在MATLAB里，用户无需对矩阵预定义就可使用。程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行。MATLAB的图形功能强大。在FORTRAN和C语言里，绘图都很不容易，但在MATLAB里，数据的可视化非常简单。MATLAB还具有较强的编辑图形界面的能力。

MATLAB的缺点是，它和其他高级程序相比，程序的执行速度较慢。由于MATLAB的程序不用编译等预处理，也不生成可执行文件，程序为解释执行，所以速度较慢。

4.5.3 MATLAB的优势特点

（1）友好的工作平台和编程环境

MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。而且新版本的MATLAB提供了完整的联机查询、帮助系统，极大的方便了用户的使用。简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。 （2）简单易用的程序语言

Matlab一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M文件）后再一起运行。新版本的MATLAB语言是基于最为流行的C＋＋语言基础上的，因此语法特征与C＋＋语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种语言可移植性好、可拓展性极强，这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因。

22

第四章 蒙特卡罗算法 （3）强大的科学计算机数据处理能力

MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能。函数中所使用的算法都是科研和工程计算中的最新研究成果，而前经过了各种优化和容错处理。在通常情况下，可以用它来代替底层编程语言，如C和C++ 。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。MATLAB的这些函数集包括从最简单最基本的函数到诸如矩阵，特征向量、快速傅立叶变换的复杂函数。函数所能解决的问题其大致包括矩阵运算和线性方程组的求解、微分方程及偏微分方程的组的求解、符号运算、傅立叶变换和数据的统计分析、工程中的优化问题、稀疏矩阵运算、复数的各种运算、三角函数和其他初等数学运算、多维数组操作以及建模动态仿真等。 （4）出色的图形处理功能

MATLAB自产生之日起就具有方便的数据可视化功能，以将向量和矩阵用图形表现出来，并且可以对图形进行标注和打印。高层次的作图包括二维和三维的可视化、图象处理、动画和表达式作图。可用于科学计算和工程绘图。新版本的MATLAB对整个图形处理功能作了很大的改进和完善，使它不仅在一般数据可视化软件都具有的功能（例如二维曲线和三维曲面的绘制和处理等）方面更加完善，而且对于一些其他软件所没有的功能（例如图形的光照处理、色度处理以及四维数据的表现等），MATLAB同样表现了出色的处理能力。同时对一些特殊的可视化要求，例如图形对话等，MATLAB也有相应的功能函数，保证了用户不同层次的要求。另外新版本的MATLAB还着重在图形用户界面（GUI）的制作上作了很大的改善，对这方面有特殊要求的用户也可以得到满足。 （5）应用广泛的模块集合工具箱

MATLAB对许多专门的领域都开发了功能强大的模块集和工具箱。一般来说，它们都是由特定领域的专家开发的，用户可以直接使用工具箱学习、应用和评估不同的方法而不需要自己编写代码。目前，MATLAB已经把工具箱延伸到了科学研究和工程应用的诸多领域，诸如数据采集、数据库接口、概率统计、样条拟合、优化算法、偏微分方程求解、神经网络、小波分析、信号处理、图像处理、系统辨识、控制系统设计、LMI控制、鲁棒控制、模型预测、模糊逻辑、金融分析、地图工具、非线性控制设计、实时快速原型及半物理仿真、嵌入式系统开发、定点仿真、DSP与通讯、电力系统仿真等，都在工具箱（Toolbox）家族中有了自己的一席之地。

（6）实用的程序接口和发布平台

23

可维护系统有效度仿真 新版本的MATLAB可以利用MATLAB编译器和C/C++数学库和图形库，将自己的MATLAB程序自动转换为独立于MATLAB运行的C和C++代码。允许用户编写可以和MATLAB进行交互的C或C++语言程序。另外，MATLAB网页服务程序还容许在Web应用中使用自己的MATLAB数学和图形程序。MATLAB的一个重要特色就是具有一套程序扩展系统和一组称之为工具箱的特殊应用子程序。工具箱是MATLAB函数的子程序库，每一个工具箱都是为某一类学科专业和应用而定制的，主要包括信号处理、控制系统、神经网络、模糊逻辑、小波分析和系统仿真等方面的应用。

（7）应用软件开发（包括用户界面）

在开发环境中，使用户更方便地控制多个文件和图形窗口；在编程方面支持了函数嵌套，有条件中断等；在图形化方面，有了更强大的图形标注和处理功能，包括对性对起连接注释等；在输入输出方面，可以直接向Excel和HDF5进行连接。

24

第五章 混联系统有效度仿真 第五章 混联系统有效度仿真

5.1 系统有效度仿真的总体思想

蒙特卡罗法在系统有效度仿真理论中的应用，实际做法就是从0-1分布中随机地抽取一系列伪随机数，再通过威布尔分布计算出相应的随机数值作为零件的样本空间，并设定取样停止时间T，将所得到的寿命进行累加，如果累加的和大于取样停止时间T，则停止取样；反之，继续取样。每一次随机模拟相当于对一类零件抽取一个进行寿命试验，通过大量重复的随机抽样及比较，就可得到零件一系列的故障时间，从而可以求得零件的有效度。抽样次数愈多，则模拟精度愈高。要获得可靠的模拟计算结果，往往要进行至少千次以上甚至上万次的模拟。因此，随机模拟需由计算机完成。

5.2 用蒙特卡罗模拟法求有效度的步骤

(1) 确定零件的总数S=1000以及采样的终止时间T=900。 (2) 给威布尔的参数m，n，?赋值，其中默认?=0。

(3) 由0-1分布随机产生一系列伪随机数RN1、RN2、RN3、RN4、?? (4) 将上述随机带入威布尔分布的反函数求出一系列产品的寿命?T1、?T2、

?T3、?T4、??

(5) 计算每一个零件累计失效的时间T1、T2、T3、T4、?? 其中：T1=?T1

T2=?T1+?T2 T3=?T1+?T2+?T3 T4=?T1+?T2+?T3+?T4 ??

(6) 将T1、T2、T3、T4、??分别与T比较，如果前者小于T则继续进行取样，

直到所得值大于T时停止采样。

(7) 重复上述(3)(4)(5)(6)得到S=1000个零件的故障分布。 (8) 设定搜索的时间间隔V=(

160~

1200)T平均寿命，这里直接取V=0.01。

(9) 搜索统计第i个时间间隔内该零件的故障个数mi，计算该间隔内零件的有效

25

可维护系统有效度仿真 度Ai?S-miS，这样的坐标共N=

TV个，其中i=1,2,3,4，?,N。

(10) 重复上述过程，得到其他两类零件的有效度计算。

(11) 将前两种零件串联在与后一种零件并联，根据下述公式计算得出整个混联系

统的有效度：

计算第一类产品的有效度：A1(i,1)=(T-m(i,1))/T 计算第二类产品的有效度：A2(i,1)=(T-m(i,1))/T 计算第三类产品的有效度：A3(i,1)=(T-m(i,1))/T

计算第一和第二类产品串联所组成子系统的有效度：A12(i,1)=A1(i,1)*A2(i,1) 计算整个混联系统的有效度：A(i,1)=1-(1-A3(i,1))*(1-A12(i,1)) (12) 将上述混联系统有效度计算由MATLAB编程并输出二维图形。

5.3 MATLAB仿真的源程序

%以下程序为模拟由三类产品组成的可维护混联系统的有效度； %首先给出第一类和二类产品串联组成的子系统的有效度模拟； %然后给出上一子系统与第三类产品并联所组成总的系统的有效度模拟； T=900;%设定采样截止时间； s=1000;%采样的样本空间；

V=0.01;%设定有效度计算时间间隔； q(s,T)=0;%将方阵结尾为数据设为0

m1=1;%设定第一种威布尔函数的形状参数初始值 n1=2;%设定威布尔函数的尺度参数初始值 for i=1:s j=1;t=0; while t

r=rand;%获取0-1分布的随机数据 t=t+wblinv(r,m1,n1);%计算产品的寿命 q(i,j)=t;%给方阵元素赋值 j=j+1; end end

26

第五章 混联系统有效度仿真 N=T/V;m=zeros(N,1);A1=zeros(N,1);%设定变量的初值为0 for i=1:s for j=1:T

k=fix(q(i,j)/V);%采用取整的方法判断产品间隔时间内的故障数 m(k,1)=m(k,1)+1;%记录产品间隔时间内的故障数 end end for i=1:N

A1(i,1)=(T-m(i,1))/T;%计算第一类产品的有效度 end

q(s,T)=0;%将方阵结尾为数据设为0

m2=1;%设定第二种产品威布尔函数的形状参数初始值 n2=3;%设定第二种产品威布尔函数的尺寸参数初始值 for i=1:s j=1;t=0; while t

r=rand;%获取0-1分布的随机数据 t=t+wblinv(r,m2,n2);%计算产品的寿命 q(i,j)=t;%给方阵元素赋值 j=j+1; end end

N=T/V;m=zeros(N,1);A2=zeros(N,1);%设定变量的初值为0 for i=1:s for j=1:T

k=fix(q(i,j)/V);%采用取整的方法判断产品间隔时间内的故障数 m(k,1)=m(k,1)+1;%记录产品间隔时间内的故障数 end end for i=1:N

27

可维护系统有效度仿真 A2(i,1)=(T-m(i,1))/T;%计算第二类产品的有效度 end for i=1:N

A12(i,1)=A1(i,1)*A2(i,1);%第一和第二类产品串联所组成子系统的有效度 end

q(s,T)=0;%将方阵结尾为数据设为0

m3=1;%设定第三种产品威布尔函数的形状参数初始值 n3=5;%设定第三种产品威布尔函数的尺寸参数初始值 for i=1:s j=1;t=0; while t

r=rand;%获取0-1分布的随机数据 t=t+wblinv(r,m3,n3);%计算产品的寿命 q(i,j)=t;%给方阵元素赋值 j=j+1; end end

N=T/V;m=zeros(N,1);A3=zeros(N,1);%设定第三类产品变量的初值为0 for i=1:s for j=1:T

k=fix(q(i,j)/V);%采用取整的方法判断产品间隔时间内的故障数 m(k,1)=m(k,1)+1;%记录产品间隔时间内的故障数 end end for i=1:N

A3(i,1)=(T-m(i,1))/T;%计算第三类产品的有效度 end for i=1:N

A(i,1)=1-(1-A3(i,1))*(1-A12(i,1));%计算整个混联系统的有效度 end

28

第五章 混联系统有效度仿真 x(1:N,1)=0:V:T-V;%设定横向坐标变量

plot(x(1:T,1),A(1:T,1));%将并联系统的有效度以图表形式表示出来

5.4 系统有效度仿真的MATLAB图形输出

5.4.1 程序代码截图

29

可维护系统有效度仿真

5.4.2 程序结果输出截图

30

第五章 混联系统有效度仿真

31

可维护系统有效度仿真 致 谢

感谢吴老师这学期陪伴我们走过了一段美好的时光。吴老师给我的感觉总是那么的亲切、随和、自然，讲授的课程也能够循序渐进、由浅入深，本人在吴老师授课期间不仅学到学科内的知识，而且也收获了很多做人的道理。感谢吴老师这段时间以来对我们所做出的付出和牺牲，在这里，我代表所有喜欢和尊敬吴老师的同学，向吴老师致以最崇高的敬意以及最真挚的感谢。

同时，我也要感谢那些帮助过我的同学们。由于我以前没有接触过MATLAB软件，所以在使用的过程中遇到了很大的麻烦。幸好有同年级的同学以及高年级的学长帮助我，给我讲解编程要点，指导我如何使用软件。他们给了我很大的帮助，当然我也占用了他们很多时间，在这里，我也向他们表达我最诚恳的谢意。

32

参考文献 参考文献

[1] 肖刚, 李天柁著. 系统可靠性分析中的蒙特卡罗方法[M]. 科学出版社, 2003 [2] 孙荣恒. 应用数理统计[M]. 科学出版社, 2003

[3] 王少萍. 工程可靠性[M]. 北京航空航天大学出版社, 2000 [4] 刘惟信. 机械可靠性设计[M]. 清华大学出版社, 1996 [5] 王超. 机械可靠性工程[M]. 冶金工业出版社, 1992

[6] 夏长俊. 寿命服从威布尔分布的装备备件储备量确定[J]. 舰船科学技术. 2006(04) [7] 白恩健. 伪随机序列构造及其随机性分析研究[D]. 西安电子科技大学 2004 [8] 覃中平. 信息安全数学基础[M]. 清华大学出版社, 2006 [9] 肖国震. 伪随机序列及其应用[M]. 国防工业出版社, 1985

[10] 丁明, 李生虎. 基于序贯蒙特卡罗仿真的配电网可靠性评估模型[J]. 电网技术. 2004(03) [11] 吴义纯. 基于蒙特卡罗仿真的风力发电系统可靠性评价[J]. 电力自动化设备. 2004(12) [12] 丁明, 戴仁昶, 刘亚成, 宋云亭, 刘盛松. 概率稳定性的蒙特卡罗仿真[J]. 清华大学学报

(自然科学版). 1999(03)

[13] 肖刚,苏光辉,李天柁,贾斗南. 计算高可靠性系统失效概率的统计估计蒙特卡罗方法[J].

核科学与工程. 2000(01)

[14] Weibull W.A statistical distribution function of wide applicability. Journal of Applied

Mechanics . 1951

[15] Xie L Y.A knowledge-based multi-dimension discrete common cause failure model. Nuclear

Engineer The . 1998

[16] Xie,L.Y.Load-property order statistics interference models for system reliability evaluation.

Int. J. Performability Engineering . 2005

[17] Hagen E W.Common-mode/common-cause failure:a review. Nuclear Engineer The . 1980 [18] Epler E P.Common mode failure considerations in the design of system for protection and

control. Nuclear Safety . 1969

[19] Fleming K N.A reliability model for common mode failures in redundant safety system.

Nuclear Engineer The . 1975

[20] Hass A.The multiple prime random number generator. ACM Transactions on Mathematical

Software . 1987

[21] An H Z.A note on chaotic maps and time series. Athens Conference on Ap-plied Probability

and Time Series:Time Series Analysis in Memory of E.J.Hannan . 1996

[22] Tausworthe R C.Random numbers generated by linear re-currence modulo two. Mathematics

of Computation . 1965

33

参考文献 参考文献

[1] 肖刚, 李天柁著. 系统可靠性分析中的蒙特卡罗方法[M]. 科学出版社, 2003 [2] 孙荣恒. 应用数理统计[M]. 科学出版社, 2003

[3] 王少萍. 工程可靠性[M]. 北京航空航天大学出版社, 2000 [4] 刘惟信. 机械可靠性设计[M]. 清华大学出版社, 1996 [5] 王超. 机械可靠性工程[M]. 冶金工业出版社, 1992

[6] 夏长俊. 寿命服从威布尔分布的装备备件储备量确定[J]. 舰船科学技术. 2006(04) [7] 白恩健. 伪随机序列构造及其随机性分析研究[D]. 西安电子科技大学 2004 [8] 覃中平. 信息安全数学基础[M]. 清华大学出版社, 2006 [9] 肖国震. 伪随机序列及其应用[M]. 国防工业出版社, 1985

[10] 丁明, 李生虎. 基于序贯蒙特卡罗仿真的配电网可靠性评估模型[J]. 电网技术. 2004(03) [11] 吴义纯. 基于蒙特卡罗仿真的风力发电系统可靠性评价[J]. 电力自动化设备. 2004(12) [12] 丁明, 戴仁昶, 刘亚成, 宋云亭, 刘盛松. 概率稳定性的蒙特卡罗仿真[J]. 清华大学学报

(自然科学版). 1999(03)

[13] 肖刚,苏光辉,李天柁,贾斗南. 计算高可靠性系统失效概率的统计估计蒙特卡罗方法[J].

核科学与工程. 2000(01)

[14] Weibull W.A statistical distribution function of wide applicability. Journal of Applied

Mechanics . 1951

[15] Xie L Y.A knowledge-based multi-dimension discrete common cause failure model. Nuclear

Engineer The . 1998

[16] Xie,L.Y.Load-property order statistics interference models for system reliability evaluation.

Int. J. Performability Engineering . 2005

[17] Hagen E W.Common-mode/common-cause failure:a review. Nuclear Engineer The . 1980 [18] Epler E P.Common mode failure considerations in the design of system for protection and

control. Nuclear Safety . 1969

[19] Fleming K N.A reliability model for common mode failures in redundant safety system.

Nuclear Engineer The . 1975

[20] Hass A.The multiple prime random number generator. ACM Transactions on Mathematical

Software . 1987

[21] An H Z.A note on chaotic maps and time series. Athens Conference on Ap-plied Probability

and Time Series:Time Series Analysis in Memory of E.J.Hannan . 1996

[22] Tausworthe R C.Random numbers generated by linear re-currence modulo two. Mathematics

of Computation . 1965

33

本文来源：https://www.bwwdw.com/article/8azp.html