蒙特卡罗方法 在可靠性分析中的应用 （一） 此处看到的幻灯片标题 2011.12.28 何劼 2014.3

内容提要 蒙特卡罗方法概述 随机数的生成 随机变量的抽样 MC在PSA不确定性分析中的应用 MC在HRA不确定性分析中的应用 小结 2

蒙特卡罗方法概述 如何测量下面不规则图形的面积？ 3

蒙特卡罗方法概述 如何近似计算圆周率π？ 当N2=30000时，π的估计值与真实值只相差0.07% 4

蒙特卡罗方法概述(续) 蒙特卡罗(Monte Carlo, MC)方法又名随机模拟法或统计试验法。它是在第二次世界大战期间兴起和发展起来的。它的奠基人是冯.诺依曼(Von Neumann)。其主要思想是在计算机上模拟实际概率过程，然后加以统计处理。该方法与传统数学方法相比，具有直观性强、简便易行的优点，它能处理一些其它方法所不能的复杂问题，并且容易在计算机上实现。 可以采用蒙特卡罗方法解决的问题分为两类，一类是问题本身具有概率或随机事件背景，可以直接用蒙特卡罗随机试验模拟；另一类是问题本身是确定性问题，需要通过构造随机事件来逼近真实情景。 有记载的最早的蒙卡方法应用实例----蒲丰(Buffon)投针试验（1777） 5

蒙特卡罗方法概述(续) 蒙特卡罗方法的优点和缺点 优点 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。 受几何条件限制小。 收敛速度与问题的维数无关。 具有同时计算多个方案与多个未知量的能力。 误差容易确定。 程序结构简单，易于实现。 缺点 收敛速度慢。 误差具有概率性。 6

蒙特卡罗方法概述(续) 蒙特卡罗方法在PSA技术中的应用 PSA技术要素 蒙特卡罗方法应用范围 始发事件分析 始发事件频率估计及不确定性分析，特殊始发事件建模及系统失效仿真 系统分析 非能动系统可靠性分析，系统失效概率计算 仪控系统可靠性建模，Markov方法中方程组求解 数据分析 设备可靠性数据的参数估计，似然方程求解 人员可靠性分析 人员操作失误概率的不确定性传播 定量化 CDF/LRF的参数不确定性分析，PSA计算软件开发 外部事件 地震脆弱性分析 火灾起火频率、水淹始发事件频率不确定性分析 外部灾害发生频率估算（外部水淹、强风） PSA应用 在役检查管道劣化机理分析，结构可靠性分析(SRA)，管道失效频率定量化 7

蒙特卡罗方法概述 随机数的生成 随机变量的抽样 MC在PSA不确定性分析中的应用 MC在HRA不确定性分析中的应用 小结 8

随机数的生成 随机数与伪随机数 (pseudorandom number) 在用蒙特卡罗方法进行系统可靠性分析和不确定性分析时，采用随机数代替设备可靠性参数作为PSA模型的输入，随机数质量的好坏直接影响定量化结果。 通常所说的随机数是指［0，1］上均匀分布的连续型随机变量的子样。然而由于受到递推（迭代）公式和计算机字长的限制，在计算机上不可能生成“真正”的随机数，计算机生成的随机数只是近似［0，1］上均匀分布的随机数，称为（伪）随机数。 常用的生成伪随机数的方法是线性同余法(LCG)，例如有乘同余法、乘加同余法、加同余法等。 9

随机数的生成(续) RiskSpectrum中不确定性分析设置界面 Seed(种子)是什么东东？ 10

随机数的生成(续) 乘加同余法 M通常取计算机中可以表示的最大的数，例如对于32位的计算机，M=232=4,294,967,296 表示axi+c除以M之后的余数 乘加同余法 M通常取计算机中可以表示的最大的数，例如对于32位的计算机，M=232=4,294,967,296 a和c为（1，M-1）区间内的正整数，均为常数。 当i=0时，初始值x0称为种子，一般可任取一个正整数。 例如，当M=16，a=5，c=13，种子x0=12345时，计算过程如下： 11

随机数的生成(续) M=16，a=5，c=13，种子x0=12345 i xi ξi 12345 — 1 12345 — 1 (5*12345+13)/16=3858余10 10/16=0.625 2 (5*10+13)/16=3余15 15/16=0.9375 3 (5*15+13)/16=5余8 8/16=0.5 4 (5*8+13)/16=3余5 5/16=0.3125 5 (5*5+13)/16=2余6 6/16=0.375 得到随机数序列为：0.625，0.9375，0.5，0.3125，0.375，… 12

随机数的生成(续) M=16，a=5，c=13，种子x0=54321 i xi ξi 54321 — 1 54321 — 1 (5*54321+13)/16=16976余2 2/16=0.125 2 (5*2+13)/16=1余7 7/16=0.4375 3 (5*7+13)/16=3余0 0/16=0 4 (5*0+13)/16=0余13 13/16=0.8125 5 (5*13+13)/16=4余14 14/16=0.875 得到随机数序列为：0.125，0.4375，0，0.8125，0.875，… 种子唯一确定了随机数序列，当种子的取值给定之后整个随机数序列就不再“随机”了； 在进行蒙特卡罗模拟时，不应当刻意选定某个数值作为种子，否则就违背了“随机试验”的理念。任何蒙特卡罗试验的结论都不应当对种子的取值具有依赖性。 13

随机数的生成(续) 生成伪随机数时需要考虑的问题： 针对上面的问题，生成伪随机数的策略是： 伪随机数的质量（独立性、均匀性）； 伪随机数序列的周期和容量； 方法的经济性。 针对上面的问题，生成伪随机数的策略是： 选取合适的递推公式T，设置合适的参数值。例如对于乘加同余法，M的所有质数因子P，应当满足a≡1(mod P) 。当M=2s(s为正整数)时，取a≡1(mod 4)，c ≡1(mod 2)时乘加同余法生成的伪随机数质量最好； 用蒙特卡罗方法进行有限次试验时，尽量避免试验次数不要超过伪随机数序列的周期； 如果调用系统自带的随机数生成函数，如果不知道其后台算法，可以先随机生成种子再根据随机的种子生成随机数，避免在同一个随机数序列中多次取随机数。 14

随机数的生成(续) M=16，a=7，c=13，种子x0=12345 i xi ξi+1 12345 — 1 12345 — 1 (7*12345+13)/16=5401余12 12/16=0.75 2 (7*12+13)/16=6余1 1/16=0.0625 3 (7*1+13)/16=1余4 4/16=0.25 4 (7*4+13)/16=2余9 9/16=0.5625 5 (7*9+13)/16=4余12 10/16=0.75 循环 15

蒙特卡罗方法概述 随机数的生成 随机变量的抽样 MC在PSA不确定性分析中的应用 MC在HRA不确定性分析中的应用 小结 16

随机变量的抽样 由已知分布的随机抽样指的是由己知分布的随机变量总体中抽取简单子样。随机数序列是由单位均匀分布的总体中抽取的简单子样，属于一种特殊的由已知分布的随机抽样问题。 为方便起见，用ξn表示[0,1]区间上均匀分布的 随机数，用XF表示由己知分布F(x)中产生的简单 子样的个体。对于连续型分布，常用分布密度函 数f(x)表示总体的己知分布，用Xf表示由己知分 布密度函数f(x)产生的简单子样的个体。 常用的抽样方法有直接抽样法、舍选法等。目 前在PSA中常用的Gamma、Beta和Lognormal分布都 已有较成熟的抽样方法。 17

随机变量的抽样(续) 直接抽样法 对于任意给定的分布函数F(x)，直接抽样方法如下： 其中，ξ1，ξ2，…，ξn为随机数序列。为方便起见，将上式简化为： 直接抽样法的优点是抽样效率高、费用低，缺点是需要F(x)和F-1(x)都能写出解析的表达式。因此直接抽样法主要适用于离散型随机变量抽样。 使得F(t)≥ξ的所有t值当中最小的 18

随机变量的抽样(续) 直接抽样法（图示） F(t) ξ1 ξ2 ξn X1 X2 X3 19

随机变量的抽样(续) 舍选法（Acceptance-Rejection Method, AR） 舍选法是Von Neumann提出的，是应用最广泛也是研究成果最多的抽样方法。设密度函数f(x)只在有限区间上为正，且最大值有限，即当a≤x ≤ b时 独立地生成两个均匀随机数ξ，η~U(0,1)，令 Y=a+(b-a) ξ 如果 M η ≤f(ξ) 则取XF =ξ，否则重新生成ξ和η。 20

随机变量的抽样(续) 舍选法的几何意义 f(x) a b 21

随机变量的抽样(续) PSA常用分布的随机抽样 分布类 抽样方案 正态分布 U，V为[0,1]区间上的随机数 对数正态分布 将正态分布随机数取指数 Gamma分布 根据a>1和0<a<1的情况分别采用舍选法，可得到Г(a,1)，根据尺度特性可得到Г(a,b) Beta分布 由Gamma分布得到，若X1~Г(a,1)，X2~Г(b,1)，且二者独立，则X1/(X1+X2) ~Beta(a,b) 22

随机变量的抽样(续) 蒙特卡罗方法的收敛性 蒙特卡罗方法是由随机变量X的简单子样X1，X2，…，XN的算术平均值： 作为所求解的近似值。由大数定律可知，如X1，X2，…，XN独立同分布，且具有有限期望值（E(X)<∞），则 即随机变量X的简单子样的算术平均值，当子样数Ｎ充分大时，以概率1收敛于它的期望值E(X)。 23

随机变量的抽样(续) 蒙特卡罗方法的误差 根据中心极限定理，有如下近似式： 蒙特卡罗方法的误差ε定义为 上式中λα表示置信度为α的标准正态分布的分位点，σ表示随机变量X的标准差，N为样本容量。误差收敛速度的阶为 。显然，当给定置信度α后，误差ε由σ和N决定。要想减小ε，要么增大N，要么减小σ。在σ固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。所以，在做PSA不确定性分析时，如果将模拟次数由5000增加一倍到10000，误差并不会随之减小一倍。 24

蒙特卡罗方法概述 随机数的生成 随机变量的抽样 MC在PSA不确定性分析中的应用 MC在HRA不确定性分析中的应用 小结 25

MC在PSA不确定性分析中的应用 不确定性分析算法流程： 26

MC在PSA不确定性分析中的应用(续) 例1. 基本事件A、B的可靠性参数如下表，如何求顶事件不可用度的不确定性区间和概率密度函数（PDF）曲线？ 失效概率 分布类型 EF A 3E-4 LogN 5 B 1E-3 10 27

MC在PSA不确定性分析中的应用(续) 根据对数正态分布的性质，将基本事件A和B的 失效概率转换为分别服从LogN(-8.59，0.98)和 LogN(-7.89，1.40)分布。 生成5000个分别服从LogN(-8.59，0.98)和 LogN(-7.89，1.40)的随机数。 由于顶门是与门，因此将A和B两个基本事件的 5000个抽样数据分别相乘，得到顶事件不可用 度的5000个样点值。 28

MC在PSA不确定性分析中的应用(续) 将顶门的5000个样点值从小到大排序，分别取 出第50、100、150、…、5000个点，共100个点 值，对应于1%、2%、3%、…、100%分位点。 以这100个分位点值作为横坐标，以1%、2%、 3%、…、100%作为纵坐标，做出的曲线就是顶 门的不可用度的累积密度函数（cdf，分布函数） 图像。 29

MC在PSA不确定性分析中的应用(续) 顶事件不可用度的累计密度函数图像 30

MC在PSA不确定性分析中的应用(续) 由于概率密度函数是累计密度函数的导函数，对 累计密度函数求数值导数，即得到概率密度函数 （pdf）。 Pi ：百分数 PVi ：分位点值 以这100个分位点值作为横坐标，以cdf的数值导 数作为纵坐标，做出的曲线就是顶门的不可用度的 概率密度函数（pdf）图像，其中5%和95%对应的 函数值就是不确定性边界。 31

MC在PSA不确定性分析中的应用(续) 顶事件不可用度的概率密度函数图像 Mean：3.21E-07 5%分位点：4.11E-09 32

MC在PSA不确定性分析中的应用(续) 顶事件不可用度的概率密度函数图像（调整后） Mean：3.21E-07 5%分位点：4.11E-09 33

MC在PSA不确定性分析中的应用(续) RiskSpectrum(1.2.0.1)生成的pdf图像 Mean：3.43E-07 34

蒙特卡罗方法概述 随机数的生成 随机变量的抽样 MC在PSA不确定性分析中的应用 MC在HRA不确定性分析中的应用 小结 35

MC在HRA不确定性分析中的应用 THERP方法，NUREG/CR-1278 第7-18页: HCR/ORE方法，EPRI TR-100259 第7-3页: 36

MC在HRA不确定性分析中的应用(续) 例2. 如何求CIB-MAN01（SGTR事故后操纵员未能隔离故障SG）HEP的误差因子？ 序 号 序 号 CIB-MAN01子任务 描述 HEP 均值 压力 等级 信息源：THERP表格20-7或20-12 恢复因子 修正的 HEP 1 遗漏关闭MSIV的步 骤 3.75E-03 5 2 7.50E-02 1.41E-03 选错MSIV的控制器 1.33E-03 22 4.99E-04 3 遗漏关闭2个排污管 阀门的步骤（中相 关） ×0.15=5.63E-04 2.11E-04 4 选错排污管线上2个 阀门的控制器（中 相关） ×0.15=2.00E‑04 7.50E-05 遗漏关闭2个 MFW/SFW阀的步骤 （中相关） 6 选错主/启动给水管 线上2个阀门的控制 器（中相关） ×0.15=2.00E-04 PE=第1项+第2项+第3项+第4项+第5项+第6项 2.48E-03 37

MC在HRA不确定性分析中的应用(续) NUREG/CR-1278中HEP数据库 序号 HEP 描述 HEP (均值) THERP 表格 (序 号) EF 遗漏失误 当正确执行有校核规定的规程时 2 长规程>10条 3.75E-03 20-7 (2) 3 从一批外形相似的控制器中选错控制器 22 一个有清晰的模拟线路的盘台 1.33E-03 20-12 (4) 10 38

事件抽样VS参数抽样？ MC在HRA不确定性分析中的应用(续) 抽样过程类似于将人员操作事件CIB-MAN01看成一棵故障树，将6个操作步骤视为基本事件。 事件抽样VS参数抽样？ 39

RiskSpectrum中对参数相关性的处理 MC在HRA不确定性分析中的应用(续) RiskSpectrum中对参数相关性的处理 RiskSpectrum PSA程序在不确定性分析部分针对使用相同参数的基本事件提供了两种蒙特卡罗抽样类型：参数抽样和事件抽样。对于调用同一个参数的基本事件，分别应用参数相关和不相关的假设。 40

MC在HRA不确定性分析中的应用(续) CIB-MAN01的HEP不确定性分析结果 参数抽样 Mean：2.44E-03 误差因子：2.78 事件抽样 Mean：2.47E-03 5%分位点：1.00E-03 95%分位点：5.04E-03 误差因子：2.24 2.48E-03 41

MC在HRA不确定性分析中的应用(续) NUREG/CR-1278中对人员操作事件HEP误差因子的估计 42

蒙特卡罗方法概述 随机数的生成 随机变量的抽样 MC在PSA不确定性分析中的应用 MC在HRA不确定性分析中的应用 小结 43

小结 蒙特卡罗方法可以模拟各种随机变量的抽样。在 运用蒙特卡罗方法时，每个环节都要搞清楚背后的 技术点，尤其对于随机数的生成、随机变量的抽样 等算法不能出差错。 在PSA不确定性分析中，对于不确定的可靠性参 数，可以用蒙卡方法模拟服从给定分布的随机数， 通过多次模拟并经过数学运算，形成顶事件的参数 的若干个样值，最终得到不确定性区间和概率密度 函数图像。 在模拟人员失误概率时，如果需要考虑操纵员动 作之间的相关性，还需要对抽样过程进行适当调整。 44

45