故障诊断学与可靠性工程.doc

此文档收集于网络 如有侵权 请联系网站删除 精品文档 机械故障诊断学及可靠性机械故障诊断学及可靠性 工程工程 MECHANICAL FAULT DIAGNOSTICS AND RELIABILITY ENGINEERING 作业题目 可维护系统有效度仿真 作者姓名作者姓名 作者学号作者学号 学科专业学科专业 机械设计及理论 指导教师指导教师 此文档收集于网络 如有侵权 请联系网站删除 精品文档 机械故障诊断学及可靠性工程机械故障诊断学及可靠性工程 可维护系统有效度仿真 硕士研究生 硕士生学号 导师 学科专业 机械设计及理论 所 在 单 位 机械工程学院 此文档收集于网络 如有侵权 请联系网站删除 精品文档 Mechanical fault diagnostics and Reliability Engineering THE SIMULATION OF VALIDITY FOR MAINTAINABLE SYSTEM 此文档收集于网络 如有侵权 请联系网站删除 精品文档 摘摘 要要 混联可维护性系统的有效性仿真是一种基于计算机计算 借助于一定的编 程软件而进行的离线模拟仿真技术 这一仿真技术 对现实的实验研究以及产 品的实际生产有着重要的预测和指导作用 本文主要针对由三类零件所组成的混联系统进行有效度的模拟仿真 当然 数学模型也做了相应简化 忽略了零件维修时间 通过 0 1 分布函数产生的伪 随机数 运用威布尔分布函数得出相应的产品使用寿命 由此 对故障零件的 故障时间进行模拟仿真 由于计算机的计算效率很快 所以模拟的零件样本空 间可以取得足够大 单个零件的采样截止时间也可以取得很大 根据蒙特卡罗 算法 通过对每一类零件的样本空间进行取样 并统计一定时间间隔的故障个 数 最终运用 MATLAB 进行编程进而模拟出整个混联系统的有效度 关键词 关键词 混联系统 威布尔分布 蒙特卡罗算法 有效度仿真 此文档收集于网络 如有侵权 请联系网站删除 精品文档 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 此文档收集于网络 如有侵权 请联系网站删除 精品文档 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 个子系统组成 当且仅当所有的子系统都能正常工作时 系统才能正常工作 这种系统称为串联系统 组成系统的所有单元中任一单元 的故障就会导致整个系统故障的系统称串联系统 它属于非贮备可靠性模型 其逻辑框图如图所示 123n 压气机压气机 燃烧室燃烧室涡轮涡轮 尾喷管尾喷管 如果系统各个子系统的可靠性分别用 表示 则系统的可 1 R 2 R n R 靠性 12 n R RRR 如果系统的各个子系统的有效度分别用 来表示 则系统的 1 A 2 A n A 有效度 12n AAAA 1 2 并联系统 假如一个系统由 n 个子系统组成 只要有一个子系统能够正常工作 系统 就能正常工作 组成系统的所有单元都故障时 系统才故障的系统叫并联系统 它属于工作贮备模型 其逻辑框图如图所示 1 1 2 2 n n 如果系统的各个子系统的可靠性分别用 来表示 则系统的可靠 1 R 2 R n R 性 22 1 1 1 1 n RRRR 如果系统的各个子系统的有效度分别用 来表示 则系统的有效 1 A 2 A n A 度 12 1 1 1 1 n AAAA 此文档收集于网络 如有侵权 请联系网站删除 精品文档 A A B C C D D 在并联系统中只有一个子系统是真正需要的 其余 n 1个子系统都被称为冗 余子系统 该系统随着冗余子系统数量的增加 其平均无故障时间也会增加 串联就是一个有问题就会瘫痪 并联只要有一个能用就没有问题 1 3 混联系统 如果一个系统既含有串联构成的子系统又含有并联构成的子系统 且两个 子系统之间或者是串联关系或者是并联关系 则该系统就是一个混联的系统 可分为串并联系统和并串联系统 串并联系统 总体串联 局部有并联 如图所示 串并联系统 等效系统 3 1 2 123 如果系统的各个子系统的有效度分别用 来表示 则统的有效度 1 A 2 A 3 A 123 1 1 1 AAAA 并串联系统 总体并联 局部有串联 如图所示 此文档收集于网络 如有侵权 请联系网站删除 精品文档 并串联系统 等效系统 12 4 12 3 34 如果系统的各个子系统的有效度分别用 来表示 则系统的有效 1 A 2 A 3 A 4 A 度 1234 1 1 1 AAAAA 此文档收集于网络 如有侵权 请联系网站删除 精品文档 第二章第二章 随机数的产生随机数的产生 2 1 随机数的定义及性质 由具有已知分布的总体中抽取简单子样 在蒙特卡罗方法中占有非常重要 的地位 总体和子样的关系 属于一般和个别的关系 或者说属于共性和个性 的关系 由具有已知分布的总体中产生简单子样 就是由简单子样中若干个性 近似地反映总体的共性 随机数是实现由已知分布抽样的基本量 在由已知分布的抽样过程中 将 随机数作为已知量 用适当的数学方法可以由它产生具有任意已知分布的简单 子样 由于随机数在蒙特卡罗方法中占有极其重要的位置 我们用专门的符号 表示 由随机数序列的定义可知 1 2 是相互独立且具有相同单位均匀 分布的随机数序列 也就是说 独立性 均匀性是随机数必备的两个特点 随机数具有非常重要的性质 对于任意自然数 s 由 s 个随机数组成的 s 维 空间上的点 n 1 n 2 n s 在 s 维空间的单位立方体 Gs 上均匀分布 即对任意的 ai 当时有如下等式成立 011 2 i ais 1 1 s n iii i Paisa 其中 P 表示事件 发生的概率 反之 如果随机变量序列 1 2 对于任 意自然数 s 由 s 个元素所组成的 s 维空间上的点 n 1 n s 在 Gs 上均 匀分布 则它们是随机数序列 由于随机数在蒙特卡罗方法中所处的特殊地位 它们虽然也属于由具有已 知分布的总体中产生简单子样的问题 但就产生方法而言 却有着本质上的差 别 2 2 随机数产生的方法 2 2 1 0 1 分布产生随机数 在连续型随机变量的分布中 最简单而且最基本的分布是单位均匀分布 由该分布抽取的简单子样称 随机数序列 其中每一个体称为随机数 单位均匀分布也称为 0 1 上的均匀分布 其分布密度函数为 1 01 0 x f x 其他 其分布函数为 此文档收集于网络 如有侵权 请联系网站删除 精品文档 0 0 01 1 1 x F xxx x 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 物理方法产生随机数 用物理方法产生随机数的基本原理是 利用某些物理现象 在计算机上增 加些特殊设备 可以在计算机上直接产生随机数 这些特殊设备称为随机数发 生器 用来作为随机数发生器的物理源主要有两种 一种是根据放射性物质的 放射性 另一种是利用计算机的固有噪声 一般情况下 任意一个随机数在计算机内总是用二进制的数表示的 12 12 222 m m 其中 i i 1 2 m 或者为0 或者为1 因此 利用物理方法在计算机上产生随机数 就是要产生只取0或1的随机 数字序列 数字之间相互独立 每个数字取0或1的概率均为0 5 用物理方法产生的随机数序列无法重复实现 不能进行程序复算 给验证 结果带来很大困难 而且 需要增加随机数发生器和电路联系等附加设备 费 用昂贵 因此 该方法也不适合在计算机上使用 2 3 伪随机数性质 判断伪随机数序列是否满足均匀和相互独立的要求 要靠统计检验的方法 实现 对于伪随机数的统计检验 一般包括两大类 均匀性检验和独立性检验 此文档收集于网络 如有侵权 请联系网站删除 精品文档 六十年代初 人们开始用定性的方法研究伪随机数序列的均匀性和独立性 问题 简要叙述如下 1 伪随机数的均匀性 这里只考虑伪随机数序列 1 2 n 全体作为子样时的均匀性问题 其中 n 为伪随机数序列的最大容量 对于任意的0 x 1 令 Nn x 表示伪随机数序列 1 2 n 中适合不等式 i0 这里 和 分别是分布状态和比例参数 因为双参数的威布尔分布有效地分析了初期致命失效 实用寿命的和耗损 此文档收集于网络 如有侵权 请联系网站删除 精品文档 阶段的寿命数据 它也可用于失效率的增长 持续和递减 定义了威布尔概率图的第一个参数是斜率 它是形状参数 因为它确定 了威布尔家族中哪一种分布相关性最好或可以描述数据 第二个参数是特征寿 命 伊塔 作为比例参数 因为它定义了分布状态的大部分 参数 和 可 从寿命数据中估计 寿命数据总为正值 威布尔分析完成后 由图可看出威布 尔概率的斜度和拟合度 注意 三参数的威布尔分布应用也很广泛 第三个参数 位置 是一个 常数 可从时间变量 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 失效原因归结于 不充足的 burn in 或应力筛选 部件的质量问题 制造的质量问题 错误的安装 设置及使用 重做 刷新时出现的问题 1 0随机失效当 1 0 失效原因归结于 维护中的人为错误 引发的失效而非固有的 意外事故和自然灾害 外来物体 闪电袭 击 强风摧毁等 1 0 4 0早期损耗当 1 0 4 0快速损耗当 4 0 除部件老化 还有以下原因引起 失效 材料的固有属性的缺陷 如陶瓷易碎 制造过程中出现的严重问题 制造或材料上的细微变化 统计学家 数学家和工程师们已将统计分布简化为数学模型或描绘出某些 行为 与其它统计分布相比 威布尔分布适于更广范围的寿命数据 威布尔概率 密度函数是一个数学函数 用以描述与数据相适应的曲线 概率密度函数可用 数学模型给出或用图形给出 其中图上 X 轴代表时间 威布尔家族中的不同成 员有不同形状的概率密度函数 累积密度函数是概率密度函数曲线下的面积 威布尔分布的累积密度函数如下公式 t tFexp1 此文档收集于网络 如有侵权 请联系网站删除 精品文档 其中 代表特征寿命 比例参数 代表斜率 状态参数 累积密度函数给出了时间 t 内的失效概率 参数 和 由失效时间进行估 计 如果失效数据来自于威布尔分布 和 的值代入累积密度函数的公式 求出一定时间内元器件的失效预计 特征寿命 和平均失效时间 MTTF 是相关的 特征寿命给出了系统或 元器件寿命中的失效概率独立于失效分布参数的点 对所有威布尔分布来说 定义为 63 2 的单元失效时的寿命 虽然 威布尔教授最初提出用平均值作为 MTTF 值绘制在威布尔概率分布图 的 y 轴上 现在是标准的工程方法用失效时间的中间值来划分寿命数据 表 7 2 展示了一个中间等级表 50 作为 10 个数量的样本 由此形成莱奥纳 多 杰克逊 Leonard Johnson 的等级公式 因为在寿命数据中非均匀分布相当常见 所以中间值比均值更为准确些 一旦知道 和 任意时间的失效概率都可轻易算出 3 4 分布类型的选择 威布尔分布有不同形式的应用 包括单参数 双参数 三参数及混合威布 尔分布 有时不属于威布尔分布的正态或对数正态分布也可用于寿命数据的分 析 选择适合于特殊数据集的分布要以数据的数量和质量 以往经验以及良好 的相关性测试为基础 下表描述了威布尔家族的各种分布 双参数威布尔 分布 双参数威布尔分布的所需参数是斜率和特征寿命 这种分布用小 样本提供了正确合理的失效分析和失效预计 它尤其可以诊断出 失效类型 例如初期损耗 尤其是电子产品 独立时间失效 意外事故和固有事件的发生 或耗损的构件 轴承 过滤器等 如果失效率递减 老练时期 burn in period 或递增 耗损阶 段 或失效率保持恒定 随机失效阶段 推荐使用双参数威布 尔分析 指数分布指数分布唯一的参数是失效率 指数分布可被视为威布尔分布的 一种特例 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 年代 E J Gumbel 第一个认真调查失效数据的极值 此文档收集于网络 如有侵权 请联系网站删除 精品文档 找到只有 6 个独立极值的分布 他的第 III 类最小极值分布与威 布尔分布相同 Gumbel Lower distribution 也叫第 I 类极小值 分布 是一种极小值分布 Gumbel upper distribution 也叫第 一类极大值分布 是一种极大值分布 当失效数据为偶然性事件 的结果并且失效数据取极值时 推荐用 Gumbel 分布 示例如自 然灾害和最大载客量等 因为 Gumbel 分布 和正态分布 可用 于预计高可靠度要求下无寿命数据情况下的负寿命值 所以在建 立寿命数据模型时要谨慎使用 尽管统计学家反对用极小样本 但安全性和明显的资金损失却决定了我们 收集的数据的局限性 当仅有极少失效数据存在时 威布尔分析可提供有用的 结果 因为 1 耗损失效发生在最陈旧的单元上 多数失效结果被绘在威布尔概率分布 图的左下角 B 0 1 至 B 1 这也正是在工程上最为关注的区域 2 威布尔分析包括失效和中止 尽管中止没有失效严重 但会有上千的中 止 在 B 0 1 B 1lives 进行更准确的工程预计 威布尔分布应用于失效机会倍增且第一次失效很重要的情况 也应用于线 性衰退而不是加速衰退系统 当威布尔分布是非线性衰退而不是当前衰退的一 个函数时 可使用对数正态分布 下表给出了正态和对数正态分布的描述 因 为即使它们不属于威布尔分布但偶尔也用来做寿命数据的参数分析 多数威布 尔软件可快速生成所有分布并自动为数据集选出一个最合适的分布 正态分布或 高斯 Gaussion 分布 正态分布的两个参数是均值和标准差 正态分布是对称的 一 般被称为贝尔曲线 该分布很重要且广泛用于概率统计中 正 态分布经常用于描述失效率随时间增加的设备 当失效时间可 用某些随机变量的总和表示时才推荐使用正态分布 正态分布 便于描述不同类型数据 它允许观测结果为负 由于时间 t 大 于零后元件才会失效 所以寿命数据总是正的 因此正态分布 并不能很好的描述寿命数据 分析人员并不为正态配置所困扰 因为利用正态分布的寿命数据也可绘制的威布尔概率分布图 对数正态分布对数正态分布的两个参数是均值和标准差 尽管对数正态分布 与正态分布相似 它假定了随机变量的对数值是正态分布而不 是随机变量本身服从正态分布 因而 所有值均为正 分布图 像便不会向左倾斜 对数正态分布很可能是威布尔分布强有力 的竞争对手 它多用于工程上的金属疲劳测试 维修数据 修 此文档收集于网络 如有侵权 请联系网站删除 精品文档 复时间 化学反应过程的仪器失效和维修 一些材料特性和 非线性加速衰退 若有倍增的失效因素影响失效时间 推荐使 用对数正态分布 例如 在衰退渐增的情况下 断裂形式是由 于压力造成 且压力随裂缝增大而增加 对数正态分布的非工 程应用包括私人收入分析 遗产继承和银行抵押 此文档收集于网络 如有侵权 请联系网站删除 精品文档 第四章第四章 蒙特卡罗算法蒙特卡罗算法 4 1 蒙特卡罗算法简介 Monte Carlo 方法 又名随机模拟法或统计实验法它是以概率统计理论为 基础 解决一些很难直接用数学运算求解或用其他方法不能解决的复杂问题的 一种近似计算法 对于那些由于计算过于复杂而难以得到解析解或者根本没有 解析解的问题 Monte Carlo 方法是一种有效的求出数值解的方法 20 世纪 40 年代以后 随着电子计算机的出现和发展 人们有可能用计算 机来模拟这类实验和计算 计算机具有计算速度高和存储容量大的特点 采用 数字模拟技术可以代替许多实际上非常庞大而复杂的实验 并迅速将实验结果 进行运算处理 于是 Monte Carlo 方法重新被提起 引起世人重视 应用日渐 广泛 实际上 采用 Monte Carlo 方法在计算机上建立模型来解决 Buffon 问 题是非常简单的 蒙特卡罗算法是以概率和统计理论方法为基础的一种计算方法 将所求解 的问题同一定的概率模型相联系 用电子计算机实现统计模拟或抽样 以获得 问题的近似解 为象征性地表明这一方法的概率统计特征 故借用赌城蒙特卡 罗命名 由 S M 乌拉姆和 J 冯 诺伊曼在 20 世纪 40 年代为研制核武器而首先 提出 为了求解数学 物理 工程技术以及管理等方面的问题 首先建立一 个概率模型或随机过程 使它们的参数 如概率分布或数学期望等问题的解 然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征 并用算 术平均值作为所求解的近似值 对于随机性问题 有时还可以根据实际物理背 景的概率法则 用电子计算机直接进行抽样试验 从而求得问题的解答 4 2 算法的基本思想及方法 传统的经验方法由于不能逼近真实的物理过程 很难得到满意的结果 而 蒙特卡罗方法由于能够真实地模拟实际物理过程 故解决问题与实际非常符合 可以得到很圆满的结果 这也是我们采用该方法的原因 蒙特卡罗方法的基本 原理及思想如下 当所要求解的问题是某种事件出现的概率 或者是某个随机 变量的期望值时 它们可以通过某种 试验 的方法 得到这种事件出现的频 率 或者这个随机变数的平均值 并用它们作为问题的解 这就是蒙特卡罗方 法的基本思想 蒙特卡罗方法通过抓住事物运动的几何数量和几何特征 利用数学方法来 加以模拟 即进行一种数字模拟实验 它是以一个概率模型为基础 按照这个 模型所描绘的过程 通过模拟实验的结果 作为问题的近似解 可以把蒙特卡 此文档收集于网络 如有侵权 请联系网站删除 精品文档 罗解题归结为三个主要步骤 构造或描述概率过程 实现从已知概率分布抽样 建立各种估计量 1 构造或描述概率过程 对于本身就具有随机性质的问题 如粒子输运问题 主要是正确描述和模 拟这个概率过程 对于本来不是随机性质的确定性问题 比如计算定积分 就 必须事先构造一个人为的概率过程 它的某些参量正好是所要求问题的解 即 要将不具有随机性质的问题转化为随机性质的问题 2 实现从已知概率分布抽样 构造了概率模型以后 由于各种概率模型都可以看作是由各种各样的概率 分布构成的 因此产生已知概率分布的随机变量 或随机向量 就成为实现蒙 特卡罗方法模拟实验的基本手段 这也是蒙特卡罗方法被称为随机抽样的原因 最简单 最基本 最重要的一个概率分布是 0 1 上的均匀分布 或称矩形分布 随机数就是具有这种均匀分布的随机变量 随机数序列就是具有这种分布的 总体的一个简单子样 也就是一个具有这种分布的相互独立的随机变数序列 产生随机数的问题 就是从这个分布的抽样问题 在计算机上 可以用物理方 法产生随机数 但价格昂贵 不能重复 使用不便 另一种方法是用数学递推 公式产生 这样产生的序列 与真正的随机数序列不同 所以称为伪随机数 或伪随机数序列 不过 经过多种统计检验表明 它与真正的随机数 或随机 数序列具有相近的性质 因此可把它作为真正的随机数来使用 由已知分布随 机抽样有各种方法 与从 0 1 上均匀分布抽样不同 这些方法都是借助于随机 序列来实现的 也就是说 都是以产生随机数为前提的 由此可见 随机数是 我们实现蒙特卡罗模拟的基本工具 建立各种估计量 一般说来 构造了概 率模型并能从中抽样后 即实现模拟实验后 我们就要确定一个随机变量 作 为所要求的问题的解 我们称它为无偏估计 4 3 算法的收敛性及误差 1 收敛性 蒙特卡罗方法是由随机变量 X 的简单子样 X1 X2 XN 的算术平均值 1 1 N Ni i XX N 作为所求解的近似值 由大数定律可知 如果 X1 X2 XN 独立同分布 且具有有限期望值 E X 则有 1 lim N N PXE X 此文档收集于网络 如有侵权 请联系网站删除 精品文档 即随机变量 X 的简单子样的算术平均值 当子样数 充分大时 以概率 1 收 N X 敛于它的期望值 E X 2 误差 蒙特卡罗方法的近似值与真值的误差问题 概率论的中心极限定理给出了 答案 该定理指出 如果随机变量序列 X1 X2 XN 独立同分布 且具有有 限非零的方差 2 即 22 0 xE Xf x dx 其中 f X 是 X 的分布密度函数 则有 2 21 2 lim x t N x N N PXE Xxedt 当N充分大时 有如下的近似式 2 2 0 2 1 2 t N PXE Xedt N 其中 称为置信度 1 称为置信水平 这表明 不等式近似地以概率 1 成立 且误差收敛 N XE X N 速度的阶为 1 2 O N 通常 蒙特卡罗方法的误差 定义为 N 上式中与置信度 是一一对应的 根据问题的要求确定出置信水平后 查标准正态分布表 就可以确定出 下面给出几个常用的 与的数值 0 50 050 003 0 67451 963 关于蒙特卡罗方法的误差需说明两点 第一 蒙特卡罗方法的误差为概率 误差 这与其他数值计算方法是有区别的 第二 误差中的均方差 是未知的 必须使用其估计值来代替 在计算所求量的同时 22 11 11 NN ii ii XX NN 可计算出 3 减小方差的方法 此文档收集于网络 如有侵权 请联系网站删除 精品文档 显然 当给定置信度 后 误差 由 和N决定 要减小 或者是 增大N 或者是减小方差 2 在 固定的情况下 要把精度提高一个数量级 试验次数N需增加两个数量级 因此 单纯增大N不是一个有效的办法 另一方面 如能减小估计的均方差 比如降低一半 那误差就减小一半 这相当于N增大四倍的效果 因此降低方差的各种技巧 引起了人们的普遍注 意 后面课程将会介绍一些降低方差的技巧 一般来说 降低方差的技巧 往往会使观察一个子样的时间增加 在固定 时间内 使观察的样本数减少 所以 一种方法的优劣 需要由方差和观察一 个子样的费用 使用计算机的时间 两者来衡量 这就是蒙特卡罗方法中效率 的概念 它定义为 其中 c 是观察一个子样的平均费用 显然越小 2 c 2 c 方法越有效 4 4 方法的特点及适用范围 4 4 1 该方法的优点和缺点 1 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程 从这个意义上讲 蒙特卡罗方法可以部分代替物理实验 甚至可以得到物 理实验难以得到的结果 用蒙特卡罗方法解决实际问题 可以直接从实际问题 本身出发 而不从方程或数学表达式出发 它有直观 形象的特点 2 受几何条件限制小 在计算 s 维空间中的任一区域 Ds 上的积分 时 无论区域 Ds 的形状多么特殊 只要能给 1212 s ss D gg x xx dx dxdx 出描述 Ds 的几何特征的条件 就可以从 Ds 中均匀产生 N 个点 12 iii s xxx 得到积分的近似值 12 1 N iii s Ns i D gg xxx N 其中 Ds 为区域 Ds 的体积 这是数值方法难以作到的 另外 在具有随机性质的问题中 如考虑的系统形状很复杂 难以用一般 数值方法求解 而使用蒙特卡罗方法 不会有原则上的困难 3 收敛速度与问题的维数无关 由误差定义可知 在给定置信水平情况下 蒙特卡罗方法的收敛速度为 与问题本身的维数无关 维数的变化 只引起抽样时间及估计量计 1 2 O N 算时间的变化 不影响误差 也就是说 使用蒙特卡罗方法时 抽取的子样总 此文档收集于网络 如有侵权 请联系网站删除 精品文档 数 N 与维数 s 无关 维数的增加 除了增加相应的计算量外 不影响问题的误 差 这一特点 决定了蒙特卡罗方法对多维问题的适应性 而一般数值方法 比如计算定积分时 计算时间随维数的幂次方而增加 而且 由于分点数与维 数的幂次方成正比 需占用相当数量的计算机内存 这些都是一般数值方法计 算高维积分时难以克服的问题 4 具有同时计算多个方案与多个未知量的能力 对于那些需要计算多个方案的问题 使用蒙特卡罗方法有时不需要像常规 方法那样逐个计算 而可以同时计算所有的方案 其全部计算量几乎与计算一 个方案的计算量相当 例如 对于屏蔽层为均匀介质的平板几何 要计算若干 种厚度的穿透概率时 只需计算最厚的一种情况 其他厚度的穿透概率在计算 最厚一种情况时稍加处理便可同时得到 另外 使用蒙特卡罗方法还可以同时得到若干个所求量 例如 在模拟粒 子过程中 可以同时得到不同区域的通量 能谱 角分布等 而不像常规方法 那样 需要逐一计算所求量 5 误差容易确定 对于一般计算方法 要给出计算结果与真值的误差并不是一件容易的事情 而蒙特卡罗方法则不然 根据蒙特卡罗方法的误差公式 可以在计算所求量的 同时计算出误差 对干很复杂的蒙特卡罗方法计算问题 也是容易确定的 一般计算方法常存在着有效位数损失问题 而要解决这一问题有时相当困 难 蒙特卡罗方法则不存在这一问题 6 程序结构简单 易于实现 在计算机上进行蒙特卡罗方法计算时 程序结构简单 分块性强 易于实 现 7 收敛速度慢 如前所述 蒙特卡罗方法的收敛速度为 一般不容易得到精确度 1 2 O N 较高的近似结果 对于维数少 三维以下 的问题 不如其他方法好 8 误差具有概率性 由于蒙特卡罗方法的误差是在一定置信水平下估计的 所以它的误差具有 概率性 而不是一般意义下的误差 9 在粒子输运问题中 计算结果与系统大小有关 经验表明 只有当系统的大小与粒子的平均自由程可以相比较时 一般在 十个平均自由程左右 蒙特卡罗方法计算的结果较为满意 但对于大系统或小 概率事件的计算问题 计算结果往往比真值偏低 而对于大系统 数值方法则 是适用的 此文档收集于网络 如有侵权 请联系网站删除 精品文档 因此 在使用蒙特卡罗方法时 可以考虑把蒙特卡罗方法与解析 或数值 方法相结合 取长补短 既能解决解析 或数值 方法难以解决的问题 也可 以解决单纯使用蒙特卡罗方法难以解决的问题 这样 可以发挥蒙特卡罗方法 的特长 使其应用范围更加广泛 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 两个英文单词的前三个字母的组合 在以后的数年里 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 语言去编写程序 尤其当涉及矩阵运算和画图时 编程会很麻烦 例如 如果用户想求解一个线性代数方程 就得编写一个程序块读入数据 然 后再使用一种求解线性方程的算法 例如追赶法 编写一个程序块来求解方程 最后再输出计算结果 在求解过程中 最麻烦的要算第二部分 解线性方程的 此文档收集于网络 如有侵权 请联系网站删除 精品文档 麻烦在于要对矩阵的元素作循环 选择稳定的算法以及代码的调试都不容易 即使有部分源代码 用户也会感到麻烦 且不能保证运算的稳定性 解线性方 程的程序用 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 的商业化以及软件本身的