用MATLAB实现蒙特卡罗法计算结构可靠度

1、用 MATLAB 实现蒙特卡罗法计算结构可靠度冯晓波 杨 桦(武汉大学水利水电学院 430072摘要 针对应用蒙特卡罗直接抽样法解决结构可靠度 所遇到的困 难 , 提 出利用 MATLAB 的 强大数 值计算 功能解 决 此类 问题 。 利 用 MATLAB 进行蒙 特卡罗抽样模拟 , 在一定程度上减少了对连续型分布 采用直接抽样时的困难 , 大大提 高 了计算效率 。关键词 蒙特卡罗法 MATLAB 可靠度 随机数1蒙特卡罗法在结构可靠性分析中运用蒙特卡罗法 1方法 , 首先考虑各基本变量相 互独立的情况。 设 基本变 量 x 1, x 2, , , x n 分别有 分布函数 F x 1 (

2、x 1 , F x2(x 2 , , , F xn(x n , 因为 F xi(x i 为 0, 1区间上的一个数 , 可以将 其与由蒙特 卡罗法产 生的随机数 r j 对应。 这样 , 便可得到 x i =F -1xi(r j , i =1, 2, , , n 。 对于每一个 r j 值 , 可以得到一 组对应 的基本变 量 x 1, x 2, , , x n 。 将这组值 代入 功 能函数 g(x 1, x 2, , , x n , 便得到一个值 , 该值 若小于 等于 0, 则 在程序中记录一次功能函数的实现 , 大于 0则不记入。 再对另一 随机数重复进行这些计算 , 直 到完 成预定

3、 的循 环次 数。 假定 所 进行的循环次数为 K 次 , g(x 1, x 2, , , x n 0的次数为 m 次 , 则可得失效概率为 P f =m/K (要求 K 足够大 。若基本 变量相 关 , 可利 用条件 概率密度 , 把多维 问题化 为 一维问题来解决 , 具 体见文 献 2。 因 此从理论 上说 , 该 方法 的 应用几乎没有什么限制。但是在 实际问 题中 , 连 续型分 布是很复 杂的。 有 的只能 给 分布函数的解析表达式 , 但给 不出 其反函 数的 解析 表达式 , 如 著名的 B 分布 ; 有的则连分 布函数 的解 析表达 式都 给不出 , 如 正态 分布。所以通常

4、情况下 对连续型 分布采 用直接 抽样是 有 一定困难的 。本文利用 MATLAB 的强大数值计算功能 , 实现了 在 MATLAB 中采用蒙特卡罗直 接抽样计算结构可 靠度 , 较好 地 解决了上述问题。2在 MATLAB 中实现蒙特卡罗法在结构可靠 度中的计算2. 1MATLAB 简介MATLAB 是一种功能极其强 大的科学和工程计算数学软件 系统。它汇 集了大量数学、 统计、 科学和工 程所需的函 数 , 除 了 具有类似于其他计算机编程语言的编程特 性外 , 对计算数学 领 域的特定数学问题 ,MATLAB 都给出了该问题的各种高效算法。 范 围 的 系 , 例如常用的矩阵代数运算、

5、数组运算等 , 使用 户可以以多 种 形式快速地操作数据 集。与 Basic 、 Fortran 、 C 等编程 语言 相比 , MATLAB 具有编程简单直 观、 用 户界面 友善、 开放性 强等 特点 , 大大提高了编程效 率。将 MATLAB 用于 蒙特 卡罗 法的 一个 显 著优点是它拥有功能强大的随机数发生器指令。2. 2随机数的产生通常在用蒙特卡罗直接 抽样法 时 , 必须 解 x =F -1(r 以 求得服从相应分布类型的 随机变 量 x , 再代 入功能 函数 求解。 而 MATLAB(6. 0版 提供了 23种随机变量分布类型 的随机数发 生器 , 如正态分布、 对数正态分布

6、、 泊松 分布、 威布尔分布 等 (这 些分布基本包括了工程 实际中 出现 的变量 分布 情况 , 可直 接 产生变量 x 以 代入功能函数 , 省去了可能会带来 很大麻烦的 求 分布函数反函数这一步 , 极大地提高 了效率。常 用指令举例 如 下 :r =r and (m, n产生 m 行 n 列的 (01 间均布随机数 组 r 。r =normrnd (Mu, Sigma, m, n产生服从 N (L , R 2 分布的 m 行 n 列随机变量数组 r 。r =lg nrnd (M u, Sigma, m, n产生 ln R 服从 N (L , R 2 分布的 m 行 n 列随机变量数组

7、r 。 MATLAB 目前版本 尚无极值 型分布的随 机数发生 器 , 故 对 极值型分布 , 仍须求解相应反函数。2. 3编程求解由于 MATLAB 能直接产生 服从 各相 应概 率分 布函 数的 随 机变量数组 , 从而使编程过程大为简化。设功能函数为 g (x 1, x 2, , , x n , x i (i =1, 2, , , n 为 服 从各相应概率分布函数的随机变量。 首 先用相应 随机数发生 器 指令产生 m c n c 的随机变量数组 x i (k , l , (i =1, 2, , , n ; k = 1, 2, , , m c ; l =1, 2, , , n c 。 然

8、后 , 将各数组中各元 素一一对应 代入功能函数。 由于 MATLAB 提供 了数组 运算 指令 , 避 免了 使 用循环语句 , 从而执行更快。 最后得到功能 函数结果数 组 , 统 计 0j , :50中国农村水利水电 #2002年第 8期P f =j m c n c如果其中某一随机变量为极值型分布 , 则先用均匀分布随 机数发生器指令 r =r and(m c , n c , 产生一个 m c 行 n c 列的均布 随机变量数组 r , 然后用式 x =F -1(r 解出服从极 值分布的随 机变量数组 , 再与用前述方法产生的其 余变量一 起代入功能 函 数求解即可。 程序框图如图 1

9、:图 1 用 MATLAB 实现蒙特卡罗法程序框图3 算 例(1 算例 1。 失效函数 g =2. 5-N C/(1+E H lg (P +$P /P 有 6个相互独立的正态变量 , 而且包含对数运算 , 具体 数值见表 1。表 1 数 据变量 N C E H P $P 均 值 1. 0000. 3961. 190168. 0003. 7200. 350分布系数0. 100. 250. 150. 050. 050. 20根据前面所述 , 直接用正态分布随 机数发生 器指令产生 服 从各相应统计特征的变 量数组元素 , 按 MATLAB 提供的数组 计 算指令代入功能函数 , 再统计小于 等于

10、0的功能函数值个数 即 可 , 非常简单直观。最后得结果为 P f =6. 4010-3, 抽样次数 16万次 , 仅 需 2、 3s 。用 JC 法计算 结果为 P f =7. 3610-3, 计 算较为麻烦。(2 算例 2。设一轴心受压 构件受 到恒定 荷载和临 时性 荷 载 2种压力 , 恒载产生的轴向力 N 1为正态分布 , L N 1=1159. 1kN , R N 1=81. 1kN ; 临时性 荷载产生的轴向力为极值 I 型分布 , L N 2=765. 5kN , R N 2=222kN ; 截面承载能力 R 为对数正态 分 布 , L R =4560kN , R R =729

11、. 6kN 。极限状态方程为 :Z =g(R, N 1, N 2 =R -N 1-N 2=0此例有一极值型分布 , 故先用均布 随机数发 生器指令产 生 (, 值分布的随机 变量。在 MATLAB 中 只要 2条 指令 即可 完成 上 述过程 , 远比 Basic 、 Fortran 、 C 等语 言简单 高效。 最后此 例结 果 为 P f =4. 2710-5。4 结 语(1 蒙特卡罗法在结构可靠度的数 值模拟中 具有收敛速 度 与基本变量的维数无关 , 极限状态函数 的复杂程 度与模拟过 程 无关且无需 将状态函数线 性化和 随机 变量 /当 量正 态化 0的 特 点 , 具有直接解决问

12、题的能力。它的普 遍适用性 使其成为当 今 可靠度理论研究的热点 , 必将在结构可 靠性分析 中发挥越来 越 大的作用。(2 在蒙特卡罗法计算中 , 当 结构失效概 率很小时 , 抽样 模 拟次数很大 , 会大大增 加计 算量 , 尤其 当功 能函 数的计 算需 要 进行结构整 体 有 限 元分 析 时 更 是 如 此。但 随 着 各 种 改 进 方 法 45和高性能计算机的出现 , 这些问题将得到较大的解决。(3 MATLAB 的强大功能为结构 可靠度计算提供了便利 , 研 究人员可迅速编出科学高效的计算程序 , 大大 提高了效率。 在 西方 ,MATLAB 的应用已遍及现 代科学 界和工

13、程界 , 相比 之下 , 我国工程界人员对其了 解和应 用就少很 多 , MATLAB 的强大 功 能等待去发掘。t参 考 文 献1 徐钟济 . 蒙特卡罗法 M. 上海 :上海科学技术出版社 , 19852 何水清 , 王善 . 结构可靠性分析与设计 M. 北京 :国防工业出版社 ,19933 张志涌 . 精通 MA TLAB5M1第 3版 1北京 :北京航空航天大学出版社 , 20004 李亚东 . 重要抽样模 拟及其在 结构可靠 度计算中 的应用 J. 西 南交通大学学报 , 1992(45 程耿东 , 蔡文学 . 结构 可靠度 计 算的 近似 重要 抽样 模拟 及其 应 用J. 工程力学 , 1997, 14(2作者简介 冯晓波 , 男 , 29岁 , 助教 , 博士研究生(收稿日期 :2002-04-02#信 息 #欢 迎 订 阅 下 列 书 刊5中国农村水利水电 61994年 , 1995年精装合订本 , 每本各 48. 00元 ; 1996年 、 1997年 、 1998年 、 1999年 精装合 订本 , 每本 各 96. 00元 ; 2000年 、 2001年精 装 合 订 本