多种结构可靠度计算方法的快速实现徐　港1,3　王　青2　王永明3(1.华中科技大学土木与力学学院,武汉430074;2.广西大学土木建筑工程学院,南宁530004;3.三峡大学土木水电学院,宜昌440332)[摘　要]　本文在总结多种结构可靠度计算方法的基础上,提出了应用Matlab 快速实现这些算法的设想,并对常用的一次二阶矩法、蒙特卡罗法以实例的形式介绍了计算过程。

[关键词]　结构可靠度;一次二阶矩法;Matlab ;蒙特卡罗法[中图分类号]　T U31112 [文献标识码]　A [文章编号]　10012523X (2004)0620007203FAST REALIZATION OF SEVERAL CALCU LATION METH ODS OFSTRUCTURAL RE LIABI LITYXu G ang Qing Wang Y ong 2ming[Abstract ]　Summing up several calculation method of structural reliability ,the thesis presents the assumption that we can realize itfleetly on Matlab ,and the fast realization of s ome usually method such as first 2order second 2m oment method and M onte Carlo method.[K eyw ords ]　S tructural reliability ;First 2order second 2m oment method ;Matlab ;M onte Carlo method 收稿日期:2004-02-28作者简介:徐　港(19742),男,内蒙古包头市人,毕业于武汉水利电力大学,现为华中科技大学在读硕士生。

1　概述可靠度的计算方法从研究的对象来说可分为点可靠度计算方法和体系可靠度计算方法。

由于可靠度研究本身的复杂性,目前对结构体系可靠度的研究还很不成熟,仍处于探索阶段。

而结构点可靠度的计算方法已较成熟,主要有:一次二阶矩法、高次高阶矩法、响应面法、蒙特卡罗法及随机限元法等[1]。

但这些方法在研究或应用中存在的一个共同难点,就是涉及到大量的数学运算。

通常的做法是利用计算机高级语言编程求解,但这样一来无疑增大了这些计算方法应用的难度。

因为它不仅要求人们要有较好的编程能力,同时还应熟练掌握各种数学算法。

那么,是否有一种能快速、准确地实现多种结构可靠度计算方法的好办法呢?经笔者实践,认为充分利用Matlab 的强大数值计算功能,便可很好地实现这一设想。

2　Matlab 简介Matlab 是由Mathw orks 公司开发的,它不仅是一个强大的集数值计算、符号运算及图形处理等功能于一体的可跨操作系统平台的科学计算软件,同时又是一种更高级,更自由的计算机语言,几乎能满足所有的计算需求。

Matlab 有20多个工具箱,如:统计工具箱、偏微分工具箱、优化工具箱、神经网络工具箱、模糊逻辑工具箱等等,汇集了大量数学、统计、科学和工程所需的函数[2]。

其中与可靠度分析最直接相关的便是统计工具箱,包含了20多种随机变量分布类型的概率分布、参数估计与假设检验、线性模型与非线性模型分析、多元统计分析、试验设计以及统计工序管理的相关函数。

下面以点可靠度分析计算中最常用的一次二阶矩法和蒙特卡罗法为例来阐述本文的观点。

3　一次二阶矩法一次二阶矩法是实际工程中最主要的计算结构可靠度的方法,按计算精度及简化条件的不同又可分为:均值一次二阶矩法、改进一次二阶矩法、JC 法及几何法等。

而其中较常用的是改进一次二阶矩法和JC 法。

改进一次二阶矩法适用于结构功能函数所含基本随机变量为独立、正态变量情况。

其主要计算难点就是解方程组困难,传统的做法无论是手算还是机算都要迭代求解,故绝大多数情况也只能求得近似解,且求解过程繁杂。

但在Matlab 中则可利用其强大的符号计算功能快速的求得精确解,如以下算例:例:已知极限状态方程为Z =g (f ,w )=fw -1140=0,且f 、w 均服从正态分布,方差μ,变异系数δ分别为:μf =38,δf =0110;μw =54,δw =0105。

求可靠指标β。

对本题详细求解过程见参考文献[3],代入相关数据运算便可得出如下方程组:cos θf =- 3.8w3(2.7f 3)2+(3.8w 3)27第31卷第6期2004年6月建　筑　技　术　开　发Building T echnique DevelopmentV ol.31,N o.6Jun.2004cosθw=- 2.7f3(2.7f3)2+(3.8w3)2f3=μf+βσf cosθf=38+3.8βcosθfw3=μw+βσw cosθw=54+2.7βcosθwz=f3w3-1140=0式中:f3,w3为验算点座标。

在Matlab中求解以上5个方程组,只须通过如下短短8句(%后为注释语句)便直接求得β的精确解为412614。

%22222222222222222222syms csl cs2w fb%这一句是定义基本变量;%22222222以下五句就是按matlab语法把上面的五个方程列出222222r=-3.8・wΠ[(2.7・f)∧2+(3.8・w)∧2]∧0.5-cs1;s=-2.7・fΠ[(2.7・f)∧2+(3.8・w)∧2]∧0.5-cs2;t=38+3.8・b・cs1-f;u=54+2.7・b・cs2-w;y=f・w-1140;%22222222222222222222result=s olve(r,s,t,u,y)　%本句为求解这五个方程组。

b=subs(result.b)　%β值。

%2222222222222222222语法简单明了,运行求解也只须4秒钟。

一般而言,结构构件功能函数不会全部由正态基本随机变量构成。

所以,实际工程中更为常用的是JC法。

JC法的基本思路是:对非正态基本随机变量作当量正态化处理,将其转化为等效正态随机变量,然后即利用改进的一次二阶矩法求结构可靠指标[4]。

我国《建筑结构设计统一标准》、《铁路工程结构设计统一标准》中即采用此法。

所以,JC法与改进一次二阶矩法的求解过程基本相同,仍然可用Matlab的符号计算方法来求解,只是要增加对非正态随机变量当量化的过程。

另外,值得说明的是,即使该法用数值迭代方式求解,用Matlab来实现比一般语言也要快捷的多,例如在求等效正态变量的均值与方差时涉及到求标准正态分布函数的反函数,则只需输入:X=NORMI NV(P,M U,SIG M A)便可得到均值为M U,方差为SIG M A,概率为P的正态分布的分位点值,根本不用去查表或编程。

再如解线性方程组时,设系数矩阵为A,常数项矩阵为D则输入:result=A\D便求出方程组解,诸如此类的例子还很多。

4　蒙特卡罗法蒙特卡罗法是最直观、精确、获取信息最多,对高次非线性问题最有效的结构可靠度计算方法。

其基本原理是对各随机变量进行大量抽样,结构失效次数占抽样数的频率即为失效概率。

该法的主要难点在于:一是随机数的生成方法:二是抽样数大小的确定。

对于问题二,主要与计算机的硬件水平有关,当然抽样技术的改进也能起到一定的作用,有关这一问题详见参考文献[3]等有关文献,而对于问题一,在Matlab中除极值型分布外各种概率分布的随机数均可由相应的随机数发生器直接得到。

例如: R=RAND(M,N)产生服从(0～1)均匀分布的m行n列的随机变量数组R。

R=NORMRND(M U,SIG M A,M,N)产生服从正态分布,N(M U,SIG M A2)的m行n列的随机变量数组R。

R=LOG NRND(M U,SIG M A,M,N)产生ln R服从正态分布,N(M U,SIG M A2)的m行n列的随机变量数组R。

而对于极值型分布的随机变量数组,在Matlab中也只需通过对(0～1)均匀分布得到的随机变量数组做简单变转化,如通过以下语句: r=rand(1,n); R=M U-0.453M U3dlt-0.77973M U3dlt3log [-log(r)];产生R服从均值为M U变异系数为dlt的极值Ⅰ型分布的n个随机变量。

有了相应分布的随机变量,则可利用Matlab的点运算功能,直接将相关数组代入功能函数进行运算,而不必像一般高级语言进行编程循环计算。

最后,统计出功能函数值中不大于0的值的个数再除以总的抽样个数,便求得失效概率。

举例如下:例:设极限状态方程Z=X1+X2-X3-X4=0,各变量的均值和变异系数X1=(2234132,011),对数正态分布;X2= (949159,011),对数正态分布;X3=(152119,01109),正态分布;X4=(49611,01292),极值Ⅰ型分布。

Matlab求解的全部语句(%部分为注释语句)如下:%222222以下一句用来提示输入抽样个数2222222n=input(′请输入抽样数N:′);%222222以下三句用来生成服从对数正态分布的随机变量X1的n个随机数222222M U=log[2234.32Πsqrt(1+0.1∧2)];sg=sqrt[log(1+0.1∧2)];x1=log nrnd(M U,sg,1,n);%222222以下三句用来生成服从对数正态分布的随机变量X2的n个随机数222222M U=log[949.59Πsqrt(1+0.1∧2)];sg=sqrt[log(1+0.1∧2)];x2=log nrnd(M U,sg,1,n);%222222下一句用来生成服从正态分布的随机变量X3的n个随机数222222x3=normrnd(1521.9,1521.9・0.109,1,n);%222222下二句用来生成服从极值Ⅰ型分布的随机变量X4的n个随机数222222tem=rand(1,n);(下转第12页)8第6期徐　港等:多种结构可靠度计算方法的快速实现第31卷图1　双层十杆衍架示意图2　多目标EP 的一次演化过程由于子目标权重向量(ω1,ω2,ω3)和构造级别不劣于关系所需的阈值cm s ,dm s ,im s ,nm s 体现了决策者的偏好序和愿意承担的风险,因此有必要在权重和阈值的不同组合下对结果进行灵敏度分析。

灵敏度分析结果如表1所示。

表1　灵敏度分析成果组合序列阈值cm sdm sim snm s权重比例偏好解Δ1Δ2W10.70.40.60.52∶3∶56831134120.70.40.60.53∶4∶35529169330.40.50.50.72∶3∶55932145440.40.50.50.73∶4∶374401265 由表1可以看出,组合1、2是在较强的阈值关系下对不同的权重比例的计算,即当决策者愿意承担较大的风险时,对权重越大的目标,决策者愿意牺牲其它子目标的利益来补偿此目标,其结果越优:但对组合3、4,在较弱的阈值关系下,决策者不愿意承担较大的风险,权重越大的目标越容易 受其它子目标竞争的包围,相反结果越劣。