故障诊断学与可靠性工程

机械故障诊断学及可靠性工程MECHANICALFAULTDIAGNOSTICSANDRELIABILITYENGINEERING作业题目可维护系统有效度仿真作者姓名作者学号学科专业机械设计及理论指导教师机械故障诊断学及可靠性工程可维护系统有效度仿真硕士研究生硕士生学号导师学科专业机械设计及理论所在单位机械工程学院MECHANICALFAULTDIAGNOSTICSANDRELIABILITYENGINEERINGTHESIMULATIONOFVALIDITYFORMAINTAINABLESYSTEM摘要混联可维护性系统的有效性仿真是一种基于计算机计算，借助于一定的编程软件而进行的离线模拟仿真技术。这一仿真技术，对现实的实验研究以及产品的实际生产有着重要的预测和指导作用。本文主要针对由三类零件所组成的混联系统进行有效度的模拟仿真，当然，数学模型也做了相应简化，忽略了零件维修时间。通过01分布函数产生的伪随机数，运用威布尔分布函数得出相应的产品使用寿命，由此，对故障零件的故障时间进行模拟仿真，由于计算机的计算效率很快，所以模拟的零件样本空间可以取得足够大，单个零件的采样截止时间也可以取得很大。根据蒙特卡罗算法，通过对每一类零件的样本空间进行取样，并统计一定时间间隔的故障个数，最终运用MATLAB进行编程进而模拟出整个混联系统的有效度。关键词混联系统；威布尔分布；蒙特卡罗算法；有效度仿真；ABSTRACTTHESIMULATIONOFVALIDITYMAINTAINABLEHYBRIDSYSTEMISAOFFLINESIMULATIONTECHNOLOGY,WHICHISBASEONCOMPUTERCALCULATIONBYUSINGCERTAINPROGRAMMINGSOFTWARETHISSIMULATIONTECHNOLOGYPLAYSANIMPORTANTROLEINTHETHEPREDICTIONANDGUIDANCEFORREALEXPERIMENTALSTUDYANDTHEACTUALPRODUCTIONOFPRODUCTSTHISPAPERMAINLYGIVESTHESIMULATIONOFVALIDITYFORTHESYSTEMWHICHISCOMPOSEDBYTHREEKINDSOFPARTSBYTHEWAY,THECORRESPONDINGMATHEMATICALMODELHASBEENSIMPLIFIED,IGNORINGTHEREPAIRINGTIMESOFTHEPARTSASERIALOFPSEUDORANDOMNUMBERSAREGENERATEDBY01DISTRIBUTIONFUNCTIONBYUSINGTHEWEIBULLDISTRIBUTIONFUNCTION,CORRESPONDINGPRODUCTUSINGLIVESAREOBTAINED,WHEREBYWECANSIMULATETHEFAULTTIMEOFTHEPARTDUEYOUTHEPERFECTCOMPUTATIONALEFFICIENCY,THESIMULATEDPARTSSAMPLESPACECANGETASBIGASPOSSIBLE,ATTHESAMETIME,THESAMPLEDEADLINEFORINDIVIDUALPARTSCANALSOGETVERYLARGEACCORDINGTOMONTECARLOALGORITHM,BYGETTINGTHENUMBEROFPARTSINEACHTYPEOFSAMPLESPACEANDTHENUMBEROFFAULTPARTSINCERTAINTIMEINTERVALTANDTHENEVENTUALLYPROGRAMMINGBYUSINGMATLAB,WECANFINALLYSIMULATETHEVALIDITYOFTHEENTIREHYBRIDSYSTEMKEYWORDSHYBRIDSYSTEMSWEIBULLDISTRIBUTIONMONTECARLOALGORITHMSIMULATIONOFVALIDITY目录摘要IABSTRACTII第一章可维护系统的组成111串联系统112并联系统113混联系统2第二章随机数的产生421随机数的定义及性质422随机数产生的方法422101分布产生随机数4222随机数表产生随机数5223物理方法产生随机数523伪随机数性质5第三章威布尔分布应用731威布尔分布简介732威布尔分布的用途733威布尔分析方法834分布类型的选择12第四章蒙特卡罗算法1541蒙特卡罗算法简介1542算法的基本思想及方法1543算法的收敛性及误差1644方法的特点及适用范围18441该方法的优点和缺点18442方法主要的应用范围2045MATLAB软件简介20451MATLAB产生的历史背景20452MATLAB的语言特点21453MATLAB的优势特点22第五章混联系统有效度仿真2551系统有效度仿真的总体思想2552用蒙特卡罗模拟法求有效度的步骤2553MATLAB仿真的源程序2654系统有效度仿真的MATLAB图形输出29541程序代码截图29542程序结果输出截图30致谢32参考文献33第一章可维护系统的组成11串联系统假设一个系统由N个子系统组成，当且仅当所有的子系统都能正常工作时，系统才能正常工作，这种系统称为串联系统。组成系统的所有单元中任一单元的故障就会导致整个系统故障的系统称串联系统。它属于非贮备可靠性模型，其逻辑框图如图所示。123N压气机燃烧室涡轮尾喷管如果系统各个子系统的可靠性分别用，表示，则系统的可1R2N靠性12NR如果系统的各个子系统的有效度分别用，来表示，则系统的1A2N有效度。12NAA12并联系统假如一个系统由N个子系统组成，只要有一个子系统能够正常工作，系统就能正常工作。组成系统的所有单元都故障时，系统才故障的系统叫并联系统，它属于工作贮备模型。其逻辑框图如图所示。12N如果系统的各个子系统的可靠性分别用,来表示，则系统的可靠1R2N性。2211NRR如果系统的各个子系统的有效度分别用,来表示，则系统的有效1A2N度。12NAAABCD在并联系统中只有一个子系统是真正需要的，其余N1个子系统都被称为冗余子系统。该系统随着冗余子系统数量的增加，其平均无故障时间也会增加。串联就是一个有问题就会瘫痪，并联只要有一个能用就没有问题。13混联系统如果一个系统既含有串联构成的子系统又含有并联构成的子系统，且两个子系统之间或者是串联关系或者是并联关系，则该系统就是一个混联的系统。可分为串并联系统和并串联系统。串并联系统总体串联，局部有并联，如图所示串并联系统等效系统312123如果系统的各个子系统的有效度分别用,来表示，则统的有效度1A23。123AA并串联系统总体并联，局部有串联，如图所示并串联系统等效系统12412334如果系统的各个子系统的有效度分别用,来表示，则系统的有效1A234度。1234AA第二章随机数的产生21随机数的定义及性质由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的关系，属于一般和个别的关系，或者说属于共性和个性的关系。由具有已知分布的总体中产生简单子样，就是由简单子样中若干个性近似地反映总体的共性。随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方法可以由它产生具有任意已知分布的简单子样。由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号表示。由随机数序列的定义可知，1，2，是相互独立且具有相同单位均匀分布的随机数序列。也就是说，独立性、均匀性是随机数必备的两个特点。随机数具有非常重要的性质对于任意自然数S，由S个随机数组成的S维空间上的点（N1，N2，NS）在S维空间的单位立方体GS上均匀分布，即对任意的AI，当时有如下等式成立01,2IA，1,SNIIIPAA其中P表示事件发生的概率。反之，如果随机变量序列1,2对于任意自然数S，由S个元素所组成的S维空间上的点（N1，NS）在GS上均匀分布，则它们是随机数序列。由于随机数在蒙特卡罗方法中所处的特殊地位，它们虽然也属于由具有已知分布的总体中产生简单子样的问题，但就产生方法而言，却有着本质上的差别。22随机数产生的方法22101分布产生随机数在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称，随机数序列，其中每一个体称为随机数。单位均匀分布也称为0,1上的均匀分布，其分布密度函数为1,0XFX其他其分布函数为0,11,XFXMATLAB中，生成（0,1）区间上均匀分布的随机变量。基本语法RANDM,N,P，其表示生成排列成MNP多维向量的随机数。如果只写M，则生成MM矩阵。222随机数表产生随机数为了产生随机数，可以使用随机数表。随机数表是由0,1,9十个数字组成，每个数字以01的等概率出现，数字之间相互独立。这些数字序列叫做随机数字序列。如果要得到N位有效数字的随机数，只需将表中每N个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为7634258910，要想得到三位有效数字的随机数依次为0763,0425,0891。因为随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。223物理方法产生随机数用物理方法产生随机数的基本原理是利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。用来作为随机数发生器的物理源主要有两种一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。一般情况下，任意一个随机数在计算机内总是用二进制的数表示的12M其中I（I1,2,M）或者为0，或者为1。因此，利用物理方法在计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为05。用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适合在计算机上使用。23伪随机数性质判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类均匀性检验和独立性检验。六十年代初，人们开始用定性的方法研究伪随机数序列的均匀性和独立性问题，简要叙述如下。（1）伪随机数的均匀性这里只考虑伪随机数序列1,2,N全体作为子样时的均匀性问题。其中N为伪随机数序列的最大容量。对于任意的0X1，令NNX表示伪随机数序列1,2,N中适合不等式I0。这里,和分别是分布状态和比例参数。因为双参数的威布尔分布有效地分析了初期致命失效，实用寿命的和耗损阶段的寿命数据，它也可用于失效率的增长，持续和递减。定义了威布尔概率图的第一个参数是斜率，它是形状参数，因为它确定了威布尔家族中哪一种分布相关性最好或可以描述数据。第二个参数是特征寿命，伊塔（）作为比例参数，因为它定义了分布状态的大部分。参数和可从寿命数据中估计，寿命数据总为正值。威布尔分析完成后，由图可看出威布尔概率的斜度和拟合度。注意三参数的威布尔分布应用也很广泛。第三个参数位置，是一个常数，可从时间变量T中加上或减去。威布尔危险函数或失效率依赖于的值，因为值说明了新或旧元件是否更有失效的可能，威布尔危险函数可以描绘出不同元件的浴盆曲线初期故障在电子和制造业中，早期失效指在使用寿命的初期失效的概率极高，当值于小10时，威布尔概率分布图表明较新的元件在正常使用时更有可能失效，被称为瞬时递减失效率。为中止电子和机械系统在早期故障的高失效率，制造商提供了产品接收测试，老练（BURNIN）早期和环境应力筛选暂不先将系统交付客户。假如有部件在初期损失阶段没有失效，那么它的失效率应当是递减的，且它的可靠度增加。因此，旧元件被认为比新元件更好，因为新元件很可能在寿命的早期失效，而元件在早期失效阶段的检修是不合适的。偶然故障假设威布尔概率分布图以一个独立失效模型为基础，为10说明失效率是常数或相对于时间独立。这意味着对于那些无故障运行至时间T的元件，在下一个单位时间内将不能保持恒定的百分比，称作恒定危险率或瞬时失效率。这使得威布尔概率图与指数分布一致。由于旧元件被认为与新元件一样好。检修通常是不适合的，唯一使系统或部件可靠度提高的方法是用随机失效进行重新设计。早期损耗在设计寿命时经常因为机械问题出现未预期的失效。当1040快速损耗当40，除部件老化，还有以下原因引起失效材料的固有属性的缺陷如陶瓷易碎制造过程中出现的严重问题制造或材料上的细微变化统计学家，数学家和工程师们已将统计分布简化为数学模型或描绘出某些行为。与其它统计分布相比,威布尔分布适于更广范围的寿命数据。威布尔概率密度函数是一个数学函数，用以描述与数据相适应的曲线。概率密度函数可用数学模型给出或用图形给出，其中图上X轴代表时间。威布尔家族中的不同成员有不同形状的概率密度函数。累积密度函数是概率密度函数曲线下的面积。威布尔分布的累积密度函数如下公式TTFEXP1其中代表特征寿命比例参数代表斜率状态参数累积密度函数给出了时间T内的失效概率参数和由失效时间进行估计，如果失效数据来自于威布尔分布，和的值代入累积密度函数的公式求出一定时间内元器件的失效预计。特征寿命和平均失效时间（MTTF）是相关的。特征寿命给出了系统或元器件寿命中的失效概率独立于失效分布参数的点。对所有威布尔分布来说，定义为632的单元失效时的寿命。虽然，威布尔教授最初提出用平均值作为MTTF值绘制在威布尔概率分布图的Y轴上，现在是标准的工程方法用失效时间的中间值来划分寿命数据。表72展示了一个中间等级表（50）作为10个数量的样本，由此形成莱奥纳多杰克逊（LEONARDJOHNSON）的等级公式。因为在寿命数据中非均匀分布相当常见，所以中间值比均值更为准确些。一旦知道和，任意时间的失效概率都可轻易算出。34分布类型的选择威布尔分布有不同形式的应用，包括单参数，双参数，三参数及混合威布尔分布，有时不属于威布尔分布的正态或对数正态分布也可用于寿命数据的分析，选择适合于特殊数据集的分布要以数据的数量和质量，以往经验以及良好的相关性测试为基础。下表描述了威布尔家族的各种分布。双参数威布尔分布双参数威布尔分布的所需参数是斜率和特征寿命。这种分布用小样本提供了正确合理的失效分析和失效预计。它尤其可以诊断出失效类型，例如初期损耗（尤其是电子产品），独立时间失效（意外事故和固有事件的发生）或耗损的构件（轴承、过滤器等），如果失效率递减（老练时期BURNINPERIOD）或递增（耗损阶段），或失效率保持恒定（随机失效阶段），推荐使用双参数威布尔分析。指数分布指数分布唯一的参数是失效率。指数分布可被视为威布尔分布的一种特例1。当部件的失效率恒定时，它的可靠度最好用威布尔分布或指数分布来描述。失效率恒定会产生无记忆属性，即一个使用过的部件的寿命与当前老化时间无关，因此可以说一个使用过的部件像一个新部件一样好（只有当为1时，威布尔分布才是无记忆的）。因为指数分布假定没有初期故障或耗损阶段，所以区域内数据要经仔细测试以确保那些假设正确。对于指数分布，MTTF与失效率互为倒数。瑞利分布瑞利分布的唯一参数是特征寿命。当2时，瑞利分布可被视为威布尔分布的一种特例。然而，它有其自身优点，是一个重要的分布，它不仅应用于可靠性问题也用于与通讯系统相关的噪声问题。作为与指数分布相似的一个单参数分布，瑞利分布可被用来描绘错误源的均方根值。如果失效率随时间线性增加时可推荐使用瑞利分布。WEIBAYES分布WEIBAYS分布的唯一参数是特性寿命，也叫单参数威布尔分布。WEIBAYES是威布尔分布的一种特例斜率定义如前所述。与BAYESIAN假定有关，当用传统的威布尔分析产生许多不确定因素时，WEIBAYS便是一个有效的解决问题的方法。当样本小于10个失效数时，WEIBOYS分布比双参数威布尔分布更精确，而且它是在失效数为零时唯一可用的分布。例如，在对一个现存失效模式进行设计修改后，从测试中获取的有效数据可用于确定新设计中威布尔曲线的置信下限被成为WEIBAYES曲线。当元器件超出其设计寿命时，无失效的威布尔分析可延长其寿命。由于WEIBAYES分布可用于无失效测试要求的情况下，所以它的重要性全在于失效对安全性和极限成本的影响。三参数威布尔分布除斜率或特征寿命参数外，三参楼威布尔分布还包括一个位置参数T0，它定义了分布在时间上的位置。这个参数可转换时间刻度的原点，而且只有被双参数分布分析证明为是合适之后，才可被使用。（有关其它信息可参考172页“威布尔概率分布图的曲线数据”）。使用位置参数时，在生成威布尔概率分布图之前，T0的值可以从时间值中减去或加上。例如，在某段时间内如果失效率为零，那么时间刻度的原点应从0转换到T0处以反映此阶段为无失效保证阶段，修正量T0为一个正值，等于失效发生的最小时间。由于正式使用前会出现寿命（可靠性）损失，T0为一个负值。负修正对于仓库中闲置元件的腐蚀的情况是有用的。例如，橡胶部件，化学品和滚珠轴承，都会随贮藏时间的延长而腐蚀老化。当T0值用于数据修正时，结果可绘成一条直线。在没有经验的情况下，用三参数威布尔分布进行分析时通常需要至少20个失效数据。GUMBEL在19世纪20年代，EJGUMBEL第一个认真调查失效数据的极值，找到只有6个独立极值的分布。他的第III类最小极值分布与威布尔分布相同。GUMBELLOWERDISTRIBUTION也叫第I类极小值分布，是一种极小值分布。GUMBELUPPERDISTRIBUTION也叫第一类极大值分布，是一种极大值分布。当失效数据为偶然性事件的结果并且失效数据取极值时，推荐用GUMBEL分布。示例如自然灾害和最大载客量等。因为GUMBEL分布（和正态分布）可用于预计高可靠度要求下无寿命数据情况下的负寿命值，所以在建立寿命数据模型时要谨慎使用。尽管统计学家反对用极小样本，但安全性和明显的资金损失却决定了我们收集的数据的局限性。当仅有极少失效数据存在时，威布尔分析可提供有用的结果，因为（1）耗损失效发生在最陈旧的单元上。多数失效结果被绘在威布尔概率分布图的左下角B01至B1，这也正是在工程上最为关注的区域。（2）威布尔分析包括失效和中止。尽管中止没有失效严重，但会有上千的中止，在B01B1LIVES进行更准确的工程预计。威布尔分布应用于失效机会倍增且第一次失效很重要的情况，也应用于线性衰退而不是加速衰退系统。当威布尔分布是非线性衰退而不是当前衰退的一个函数时，可使用对数正态分布。下表给出了正态和对数正态分布的描述。因为即使它们不属于威布尔分布但偶尔也用来做寿命数据的参数分析。多数威布尔软件可快速生成所有分布并自动为数据集选出一个最合适的分布。正态分布或高斯GAUSSION分布正态分布的两个参数是均值和标准差。正态分布是对称的，一般被称为贝尔曲线，该分布很重要且广泛用于概率统计中。正态分布经常用于描述失效率随时间增加的设备。当失效时间可用某些随机变量的总和表示时才推荐使用正态分布。正态分布便于描述不同类型数据，它允许观测结果为负。由于时间T大于零后元件才会失效，所以寿命数据总是正的。因此正态分布并不能很好的描述寿命数据。分析人员并不为正态配置所困扰，因为利用正态分布的寿命数据也可绘制的威布尔概率分布图。对数正态分布对数正态分布的两个参数是均值和标准差。尽管对数正态分布与正态分布相似。它假定了随机变量的对数值是正态分布而不是随机变量本身服从正态分布。因而，所有值均为正，分布图像便不会向左倾斜。对数正态分布很可能是威布尔分布强有力的竞争对手。它多用于工程上的金属疲劳测试、维修数据（修复时间）、化学反应过程的仪器失效和维修，一些材料特性和非线性加速衰退。若有倍增的失效因素影响失效时间，推荐使用对数正态分布。例如，在衰退渐增的情况下，断裂形式是由于压力造成，且压力随裂缝增大而增加。对数正态分布的非工程应用包括私人收入分析，遗产继承和银行抵押。第四章蒙特卡罗算法41蒙特卡罗算法简介MONTECARLO方法，又名随机模拟法或统计实验法它是以概率统计理论为基础，解决一些很难直接用数学运算求解或用其他方法不能解决的复杂问题的一种近似计算法。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，MONTECARLO方法是一种有效的求出数值解的方法。20世纪40年代以后，随着电子计算机的出现和发展，人们有可能用计算机来模拟这类实验和计算。计算机具有计算速度高和存储容量大的特点，采用数字模拟技术可以代替许多实际上非常庞大而复杂的实验，并迅速将实验结果进行运算处理，于是MONTECARLO方法重新被提起，引起世人重视，应用日渐广泛。实际上，采用MONTECARLO方法在计算机上建立模型来解决BUFFON问题是非常简单的。蒙特卡罗算法是以概率和统计理论方法为基础的一种计算方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。由SM乌拉姆和J冯诺伊曼在20世纪40年代为研制核武器而首先提出。为了求解数学、物理、工程技术以及管理等方面的问题，首先建立一个概率模型或随机过程，使它们的参数，如概率分布或数学期望等问题的解；然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征，并用算术平均值作为所求解的近似值。对于随机性问题，有时还可以根据实际物理背景的概率法则，用电子计算机直接进行抽样试验，从而求得问题的解答。42算法的基本思想及方法传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。这也是我们采用该方法的原因。蒙特卡罗方法的基本原理及思想如下当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。这就是蒙特卡罗方法的基本思想。蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。它是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。可以把蒙特卡罗解题归结为三个主要步骤构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。（1）构造或描述概率过程对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。（2）实现从已知概率分布抽样构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是0,1上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从0,1上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。建立各种估计量一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。43算法的收敛性及误差1收敛性蒙特卡罗方法是由随机变量X的简单子样X1，X2，XN的算术平均值1NII作为所求解的近似值。由大数定律可知，如果X1，X2，XN独立同分布，且具有有限期望值（EX），则有1LIMNPXE即随机变量X的简单子样的算术平均值，当子样数充分大时，以概率1收NX敛于它的期望值EX。2误差蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，XN独立同分布，且具有有限非零的方差2，即220XEXFXD其中，FX是X的分布密度函数。则有2/1LIMXTNNPE当N充分大时，有如下的近似式2/01TNXEED其中称为置信度，1称为置信水平。这表明，不等式近似地以概率1成立，且误差收敛N速度的阶为。1/2O通常，蒙特卡罗方法的误差定义为N上式中与置信度是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出。下面给出几个常用的与的数值050050003067451963关于蒙特卡罗方法的误差需说明两点第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差是未知的，必须使用其估计值来代替，在计算所求量的同时，2211NNIIIIX可计算出。3减小方差的方法显然，当给定置信度后，误差由和N决定。要减小，或者是增大N，或者是减小方差2。在固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。另一方面，如能减小估计的均方差，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。因此降低方差的各种技巧，引起了人们的普遍注意。后面课程将会介绍一些降低方差的技巧。一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中效率的概念。它定义为，其中C是观察一个子样的平均费用。显然越小，22C方法越有效。44方法的特点及适用范围441该方法的优点和缺点（1）能够比较逼真地描述具有随机性质的事物的特点及物理实验过程从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。（2）受几何条件限制小在计算S维空间中的任一区域DS上的积分时，无论区域DS的形状多么特殊，只要能给1212,SSSDGXDX出描述DS的几何特征的条件，就可以从DS中均匀产生N个点，12,IISX得到积分的近似值。12,NIISSIGX其中DS为区域DS的体积。这是数值方法难以作到的。另外，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。（3）收敛速度与问题的维数无关由误差定义可知，在给定置信水平情况下，蒙特卡罗方法的收敛速度为，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计1/2ON算时间的变化，不影响误差。也就是说，使用蒙特卡罗方法时，抽取的子样总数N与维数S无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。（4）具有同时计算多个方案与多个未知量的能力对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。（5）误差容易确定对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。一般计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。（6）程序结构简单，易于实现在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。（7）收敛速度慢如前所述，蒙特卡罗方法的收敛速度为，一般不容易得到精确度1/2ON较高的近似结果。对于维数少（三维以下）的问题，不如其他方法好。（8）误差具有概率性由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。（9）在粒子输运问题中，计算结果与系统大小有关经验表明，只有当系统的大小与粒子的平均自由程可以相比较时（一般在十个平均自由程左右），蒙特卡罗方法计算的结果较为满意。但对于大系统或小概率事件的计算问题，计算结果往往比真值偏低。而对于大系统，数值方法则是适用的。因此，在使用蒙特卡罗方法时，可以考虑把蒙特卡罗方法与解析（或数值）方法相结合，取长补短，既能解决解析（或数值）方法难以解决的问题，也可以解决单纯使用蒙特卡罗方法难以解决的问题。这样，可以发挥蒙特卡罗方法的特长，使其应用范围更加广泛。442方法主要的应用范围蒙特卡罗方法所特有的优点，使得它的应用范围越来越广。它的主要应用范围包括粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。随着科学技术的发展，其应用范围将更加广泛。蒙特卡罗方法在粒子输运问题中的应用范围主要包括实验核物理，反应堆物理，高能物理等方面。蒙特卡罗方法有很强的适应性，问题的几何形状的复杂性对它的影响不大。该方法的收敛性是指概率意义下的收敛，因此问题维数的增加不会影响它的收敛速度，而且存贮单元也很省，这些是用该方法处理大型复杂问题时的优势。因此，随着电子计算机的发展和科学技术问题的日趋复杂，蒙特卡罗方法的应用也越来越广泛。它不仅较好地解决了多重积分计算、微分方程求解、积分方程求解、特征值计算和非线性方程组求解等高难度和复杂的数学计算问题，而且在统计物理、核物理、真空技术、系统科学、信息科学、公用事业、地质、医学，可靠性及计算机科学等广泛的领域都得到成功的应用。蒙特卡罗方法在实验核物理中的应用范围主要包括通量及反应率，中子探测效率，光子探测效率，光子能量沉积谱及响应函数，气体正比计数管反冲质子谱，多次散射与通量衰减修正等方面。45MATLAB软件简介451MATLAB产生的历史背景在20世纪70年代中期，CLEVEMOLER博士和其同事在美国国家科学基金的资助下开发了调用EISPACK和LINPACK的FORTRAN子程序库。EISPACK是特征值求解的FORTRAN程序库，LINPACK是解线性方程的程序库。在当时，这两个程序库代表矩阵运算的最高水平。到20世纪70年代后期，身为美国NEWMEXICO大学计算机系系主任的CLEVEMOLER，在给学生讲授线性代数课程时，想教学生使用EISPACK和LINPACK程序库，但他发现学生用FORTRAN编写接口程序很费时间，于是他开始自己动手，利用业余时间为学生编写EISPACK和LINPACK的接口程序。CLEVEMOLER给这个接口程序取名为MATLAB，该名为矩阵（MATRIX）和实验室（LABORATORY）两个英文单词的前三个字母的组合。在以后的数年里，MATLAB在多所大学里作为教学辅助软件使用，并作为面向大众的免费软件广为流传。1983年春天，CLEVEMOLER到STANFORD大学讲学，MATLAB深深地吸引了工程师JOHNLITTLE。JOHNLITTLE敏锐地觉察到MATLAB在工程领域的广阔前景。同年，他和CLEVEMOLER、SIEVEBANGERT一起，用C语言开发了第二代专业版。这一代的MATLAB语言同时具备了数值计算和数据图示化的功能。1984年，CLEVEMOLER和JOHNLITHE成立了MATHWORKS公司，正式把MATLAB推向市场，并继续进行MATLAB的研究和开发。在当今30多个数学类科技应用软件中，就软件数学处理的原始内核而言，可分为两大类。一类是数值计算型软件，如MATLAB、XMATH、GAUSS等，这类软件长于数值计算，对处理大批数据效率高；另一类是数学分析型软件，如MATHEMATICA、MAPLE等，这类软件以符号计算见长，能给出解析解和任意精度解，其缺点是处理大量数据时效率较低。MATHWORKS公司顺应多功能需求之潮流，在其卓越数值计算和图示能力的基础上，又率先在专业水平上开拓了其符号计算、文字处理、可视化建模和实时控制能力，开发了适合多学科、多部门要求的新一代科技应用软件MATLAB。经过多年的国际竞争，MATLAB已经占据了数值型软件市场的主导地位。452MATLAB的语言特点一种语言之所以能如此迅速地普及，显示出如此旺盛的生命力，是由于它有着不同于其他语言的特点。正如同FORTRAN和C等高级语言使人们摆脱了需要直接对计算机硬件资源进行操作一样，被称作为第四代计算机语言的MATLAB，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来。MATLAB的最突出的特点就是简洁。MATLAB用更直观的、符合人们思维习惯的代码，代替了C和FORTRAN语言的冗长代码。MATLAB给用户带来的是最直观、最简洁的程序开发环境。以下简单介绍一下MATLAB的主要特点。语言简洁紧凑，使用方便灵活，库函数极其丰富。MATLAB程序书写形式自由，利用其丰富的库函数避开繁杂的子程序编程任务，压缩了一切不必要的编程工作。由于库函数都由本领域的专家编写，用户不必担心函数的可靠性。可以说，用MATLAB进行科技开发是站在专家的肩膀上。具有FORTRAN和C等高级计算机语言知识的读者可能已经注意到，如果用FORTRAN或C语言去编写程序，尤其当涉及矩阵运算和画图时，编程会很麻烦。例如，如果用户想求解一个线性代数方程，就得编写一个程序块读入数据，然后再使用一种求解线性方程的算法（例如追赶法）编写一个程序块来求解方程，最后再输出计算结果。在求解过程中，最麻烦的要算第二部分。解线性方程的麻烦在于要对矩阵的元素作循环，选择稳定的算法以及代码的调试都不容易。即使有部分源代码，用户也会感到麻烦，且不能保证运算的稳定性。解线性方程的程序用FORTRAN和C这样的高级语言编写至少需要好几十行。再如用双步QR方法求解矩阵特征值，如果用FORTRAN编写，至少需要四百多行，调试这种几百行的计算程序可以说很困难。由于MATLAB是用C语言编写的，MATLAB提供了和C语言几乎一样多的运算符，灵活使用MATLAB的运算符将使程序变得极为简短。MATLAB既具有结构化的控制语句（如FOR循环、WHILE循环、BREAK语句和IF语句），又有面向对象编程的特性。语法限制不严格，程序设计自由度大。例如，在MATLAB里，用户无需对矩阵预定义就可使用。程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行。MATLAB的图形功能强大。在FORTRAN和C语言里，绘图都很不容易，但在MATLAB里，数据的可视化非常简单。MATLAB还具有较强的编辑图形界面的能力。MATLAB的缺点是，它和其他高级程序相比，程序的执行速度较慢。由于MATLAB的程序不用编译等预处理，也不生成可执行文件，程序为解释执行，所以速度较慢。453MATLAB的优势特点（1）友好的工作平台和编程环境MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近WINDOWS的标准界面，人机交互性更强，操作更简单。而且新版本的MATLAB提供了完整的联机查询、帮助系统，极大的方便了用户的使用。简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。（2）简单易用的程序语言MATLAB一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M文件）后再一起运行。新版本的MATLAB语言是基于最为流行的C语言基础上的，因此语法特征与C语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种语言可移植性好、可拓展性极强，这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因。（3）强大的科学计算机数据处理能力MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能。函数中所使用的算法都是科研和工程计算中的最新研究成果，而前经过了各种优化和容错处理。在通常情况下，可以用它来代替底层编程语言，如C和C。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。MATLAB的这些函数集包括从最简单最基本的函数到诸如矩阵，特征向量、快速傅立叶变换的复杂函数。函数所能解决的问题其大致包括矩阵运算和线性方程组的求解、微分方程及偏微分方程的组的求解、符号运算、傅立叶变换和数据的统计分析、工程中的优化问题、稀疏矩阵运算、复数的各种运算、三角函数和其他初等数学运算、多维数组操作以及建模动态仿真等。（4）出色的图形处理功能MATLAB自产生之日起就具有方便的数据可视化功能，以将向量和矩阵用图形表现出来，并且可以对图形进行标注和打印。高层次的作图包括二维和三维的可视化、图象处理、动画和表达式作图。可用于科学计算和工程绘图。新版本的MATLAB对整个图形处理功能作了很大的改进和完善，使它不仅在一般数据可视化软件都具有的功能（例如二维曲线和三维曲面的绘制和处理等）方面更加完善，而且对于一些其他软件所没有的功能（例如图形的光照处理、色度处理以及四维数据的表现等），MATLAB同样表现了出色的处理能力。同时对一些特殊的可视化要求，例如图形对话等，MATLAB也有相应的功能函数，保证了用户不同层次的要求。另外新版本的MATLAB还着重在图形用户界面（GUI）的制作上作了很大的改善，对这方面有特殊要求的用户也可以得到满足。（5）应用广泛的模块集合工具箱MATLAB对许多专门的领域都开发了功能强大的模块集和工具箱。一般来说，它们都是由特定领域的专家开发的，用户可以直接使用工具箱学习、应用和评估不同的方法而不需要自己编写代码。目前，MATLAB已经把工具箱延伸到了科学研究和工程应用的诸多领域，诸如数据采集、数据库接口、概率统计、样条拟合、优化算法、偏微分方程求解、神经网络、小波分析、信号处理、图像处理、系统辨识、控制系统设计、LMI控制、鲁棒控制、模型预测、模糊逻辑、金融分析、地图工具、非线性控制设计、实时快速原型及半物理仿真、嵌入式系统开发、定点仿真、DSP与通讯、电力系统仿真等，都在工具箱（TOOLBOX）家族中有了自己的一席之地。（6）实用的程序接口和发布平台新版本的MATLAB可以利用MATLAB编译器和C/C数学库和图形库，将自己的MATLAB程序自动转换为独立于MATLAB运行的C和C代码。允许用户编写可以和MATLAB进行交互的C或C语言程序。另外，MATLAB网页服务程序还容许在WEB应用中使用自己的MATLAB数学和图形程序。MATLAB的一个重要特色就是具有一套程序扩展系统和一组称之为工具箱的特殊应用子程序。工具箱是MATLAB函数的子程序库，每一个工具箱都是为某一类学科专业和应用而定制的，主要包括信号处理、控制系统、神经网络、模糊逻辑、小波分析和系统仿真等方面的应用。（7）应用软件开发（包括用户界面）在开发环境中，使用户更方便地控制多个文件和图形窗口；在编程方面支持了函数嵌套，有条件中断等；在图形化方面，有了更强大的图形标注和处理功能，包括对性对起连接注释等；在输入输出方面，可以直接向EXCEL和HDF5进行连接。第五章混联系统有效度仿真51系统有效度仿真的总体思想蒙特卡罗法在系统有效度仿真理论中的应用，实际做法就是从01分布中随机地抽取一系列伪随机数，再通过威布尔分布计算出相应的随机数值作为零件的样本空间，并设定取样停止时间T，将所得到的寿命进行累加，如果累加的和大于取样停止时间T，则停止取样；反之，继续取样。每一次随机模拟相当于对一类零件抽取一个进行寿命试验，通过大量重复的随机抽样及比较，就可得到零件一系列的故障时间，从而可以求得零件的有效度。抽样次数愈多，则模拟精度愈高。要获得可靠的模拟计算结果，往往要进行至少千次以上甚至上万次的模拟。因此，随机模拟需由计算机完成。52用蒙特卡罗模拟法求有效度的步骤1确定零件的总数S1000以及采样的终止时间T900。2给威布尔的参数M，N，赋值，其中默认0。3由01分布随机产生一系列伪随机数RN1、RN2、RN3、RN4、4将上述随机带入威布尔分布的反函数求出一系列产品的寿命1、2、T3、4、T5计算每一个零件累计失效的时间T1、T2、T3、T4、其中T11T212T3123TT41234T6将T1、T2、T3、T4、分别与T比较，如果前者小于T则继续进行取样，直到所得值大于T时停止采样。7重复上述3456得到S1000个零件的故障分布。8设定搜索的时间间隔V，这里直接取V001。1602平均寿命9搜索统计第I个时间间隔内该零件的故障个数MI，计算该间隔内零件的有效度，这样的坐标共N个，其中I1,2,3,4，,N。IIMASTV10重复上述过程，得到其他两类零件的有效度计算。11将前两种零件串联在与后一种零件并联，根据下述公式计算得出整个混联系统的有效度计算第一类产品的有效度A1I,1TMI,1/T计算第二类产品的有效度A2I,1TMI,1/T计算第三类产品的有效度A3I,1TMI,1/T计算第一和第二类产品串联所组成子系统的有效度A12I,1A1I,1A2I,1计算整个混联系统的有效度AI,111A3I,11A12I,112将上述混