结构可靠度基于Matlab算法性能比较.docx
《结构可靠度基于Matlab算法性能比较.docx》由会员分享，可在线阅读，更多相关《结构可靠度基于Matlab算法性能比较.docx（19页珍藏版）》请在冰豆网上搜索。
结构可靠度基于Matlab算法性能比较
结构可靠度基于Matlab算法性能比较
摘要：
本文对只有一个失效模式的构件可靠度问题的各类算法进行比较，并对各种算法（JC法、Breitung法、一次渐进法、重要抽样蒙特卡洛法、蒙特卡洛法）发展历史和推导过程做出论述。
基于结构可靠度分析方法（一次二阶矩法，二次二阶矩法，蒙特卡洛法）的基本原理,采取不同的算例利用MATLAB软件编制了相应计算程序，并对计算结果进行了比较。
本文算例对比主要分为三大类：
1.线性功能函数与非线性功能函数的对比；2.相同的功能函数相同的分布不同的参数；3.相同的功能函数相同的参数不同的分布。
并通过对算例的比较得出每种方法的适用范围和精确程度，为以后工程中解决可靠度问题提供了理论和实例依据。
关键词：
结构可靠度；Matlab;JC法；Breitung法；一次渐进法；重要抽样蒙特卡罗法；蒙特卡罗法；
1绪论
1.1概述
为了保证结构、构件等物体的适用性和耐久性要求，工程人员要对构件物体的受力进行可靠度分析。
可靠度分析采用容许应力法即对材料强度进行一系列折减（即当材料所能承受的最大力为F,在工程应用中我们将F除以一个大于1的系数n，采用除数作为构件所能承受的最大的力进行受力分析。
但是这种方法有自身的缺陷性，例如对于一个受弯构建来讲，当力F增大2倍时，弯矩增大量却大于2倍。
因此，用容许应力法做结构可靠度设计有自身的缺陷，现在工程中已经较少采用。
随着概率论与数理统计的发展，分项系数法逐渐在可靠度分析中采用，分项系数法不像容许应力法那样只对应力进行一项折减，而是分别对不同的力不同的材料采取不同的分项系数进行可靠度计算，分项系数法现在在工程可靠度分析中广泛采用。
由于材料生产的不均匀性和结构作用荷载的随机性，导致抗力和荷载都属于符合某种情形的分布，通过概率论知识并依托Matlab数学分析软件我们可以计算出结构的可靠指标和失效概率，确定结构的安全程度。
1.2国内外研究现状
可靠度理论在国外开展的比较早，我国知道20世纪70年代才开始可靠度方面的研究，虽说我国开展的较晚但也做了大量的研究并提出了一些自己的理论，目前我国从事结构可靠度方面的研究仍然比较多。
可靠度理论在结构工程方面的应用研究开始得较早，目前已广泛应用于水利、煤炭、通信等多种行业。
早在1947年，前苏联的尔然尼钦【1】通过大量研究就提出了用一次二阶矩方法来估计结构的失效概率。
但那以后的一段时间里，这种计算失效概率的方法并没有得到发展仍然停留在古典可靠性理论，而且只适用于随机变量符合标准正太分布的情况。
其后，美国的C.A.Comell、A.H.s【2】和wH.Tang【3】在总结了大量的工程实践经验发展了工程技术中应用概率的概念和方法。
其中CA.Comell于1969年提出了作为衡量结构安全度的统一标准—与结构失效概率相联系的可靠度指标，并建立了结构安全度的二阶矩模式。
Lind【4】在总结前人经验的基础上提出了材料和构件分项系数的概念，将可靠度指标表达为设计人员习惯使用的分项系数形式。
随后美国伊利诺斯大学A.H.S在结构可靠度方面做出了较大的进步，提出了广义可靠度概率法，并对各种结构的分布情况做了分析。
1976年国际“结构安全度联合委员会”（Jcss）采用RackwitZ【5】和Fiessler【6】等提出的“当量正态”法以考虑随机变量实际分布的二阶矩模式，这种方法极大的提高了二阶距方法的精度。
在此后一段时间可靠度理论越来越完善，在结构设计中的应用也越来越广泛。
可靠度理论迎来的发展的黄金期。
20世纪70年代以来，各个国家都开始了基于概率论与数理统计理论的可靠度分析体制和相应的设计规范的建立。
加拿大【7】、联邦德国【8】、北欧五国、美国等国家先后建立了自己的可靠度设计规范，我国也在90年代建立了自己的可靠度设计规范。
目前我国土木工作者依托大量的工作实践提出了自己对于可靠度的理解，为世界范围内可靠度理论的发展做出了自己的贡献。
从可靠度理论开始应用于结构工程方面，世界范围内的科学家依托自己对于可靠度理论的理解并结合相应的工程实践提出了许多根据可靠度理论计算可靠指标和失效概率的方法，例如：
一次二阶距法（中心点法、JC法、实用分析法等）、二次二阶距法【9】、二次四阶矩法【10】、一次渐进法、蒙特卡洛法【11】（直接蒙特卡洛法、重要抽样蒙特卡洛法、方向抽样蒙特卡洛法等）等等方法。
每种方法都有自身的实用范围和优缺点，要总结相应方法的适用范围，在不同的情况下采用不同的计算方法，获取更加精确的结果。
2结构可靠度基本理论
2.1可靠性概念
从事工程结构设计的基本目的，是在一定的经济条件下，赋予结构以适当的可靠度，使结构在预定的使用期限内，能满足设计所预期的各种功能要求。
无论是房屋、桥梁、隧洞等结构，都必须满足以下四项基本功能要求【12】：
1.能承受在正常施工和正常使用时可能出现的各种作用；
2.在正常使用时具有良好的工作性能；
3.在正常维护下具有足够的耐久性能；
4.在偶然事件发生时（如地震、火灾等）及发生后，仍能保持必需的整体稳定性。
上述第1，4项为结构的安全性要求，第2项为结构的适用性要求，第3项是结构的耐久性要求。
结构若同时满足安全性，适用性和耐久性要求，则称该结构可靠，即结构的可靠性是结构安全性、适用性和耐久性的统称。
所谓结构的耐久性能，是指在一定的期限内结构材料和主体结构没有发生过大的老化或者变形影响适用者的舒适度和房屋失效概率增加。
所谓结构的安全性，是指结构在遇到偶然荷载（例如50年一遇的地震、大风等），结构不至于倒塌对人的生命财产造成一些损失。
可靠度的定义【13】：
结构在规定的时间内，在规定的条件下，完成预定功能的能力，称为结构可靠性。
对于特定情况而言，可具体解释成在规定的参照时期内，结构将不会达到某一特定极限状态的概率，即所谓的结构可靠度。
结构可靠度是结构可靠性的概率量度。
上述“规定的时间”，一般指结构设计基准期，目前世界上大多数国家普通结构的设计基准期均为50年。
由于荷载效应一般随设计基准期增长而增大（假如一个地区遇到7级地震的概率是20年一遇，那么设计基准期选为50年结构遇到7级地震的概率和次数就要比基准期为30年的结构要大，对结构造成的破坏也就更大），抵抗能力随着基准期增长降低，而影响结构抗力的材料性能指标则随设计基准期的增大而减小，因此结构可靠度与“规定的时间”有关，“规定的时间”越长，结构的可靠度越低。
一般情况下，总可以将影响结构可靠度的因素归纳为两个综合量，即结构或结构构件的荷载效应S和抗力R。
令
Z=g（R，S）=R—S
由于R和S都是随机变量，因此Z是一个随机变量，Z可能出现下列三种情况：
Z>0,结构可靠
Z<0,结构实效
Z=0,结构处于极限状态
结构可靠性设计既是求出Z大于0的概率即可靠度，Z小于0对应于结构的失效概率。
对于结构抗力R由于构件材料、尺寸、制造工艺的差别，造成R并不是一个固定不变的数而是一个关于材料、尺寸、制造差异的一个随机变量。
S也是一个随机变量（例如一个房屋所承受的雪荷载、风荷载等值并不是固定的一个数，我们可以通过对不同月份和年份的统计出雪荷载的大小，并根据数理统计原理统计出S的分布情况来）。
结构的极限状态方程如下：
Z=R—S
（1）
假设R和S为两个相互独立的正态随机变量，他们的均值和方差分别为μR，μS；σR，σS。
根据概率论知识可以得到：
μZ=μR-μS
（2）
σZ2=σR2+σS2（3）
令β=μZ/σZ，称β为结构的可靠指标。
β越小结构的失效概率越大；β越大，结构的失效概率越小，结构越安全。
2.2结构可靠度与失效概率
对于公式Z=R-S计算P（Z>0）,在概率论与数理统计中先求出极限状态方程的联合分布密度再求出极限状态方程的Z>0的区域，对联合分布密度函数在该区域内积分求出极限状态方程的概率。
对于Z<0即对联合分布密度在Z<0的区域内积分。
对于结构可靠度也采用同上面相应算法表示为Pr，fX（X1,X2,…Xn）为结构的联合密度函数，Z>0为联合分布函数的积分区域
（4）
相反，如果结构不能完成预定的功能，则称相应的概率为结构失效的概率，表示为Pf，即：
（5）
结构的可靠与失效为两个互不相容事件，因此，由概率论可知，结构的可靠概率Pr与失效概率Pf是互补的，即：
Pr+Pf=1（6）
根据概率论的只是我们知道对于比较复杂的分布情况直接求出功能函数联合密度函数的难度较大，因此为了简化可靠度计算的复杂程度提出了许多种近似计算方法。
下面章节将对本文所涉及的可靠度计算方法做出详细的介绍和对比。
3结构可靠度分析方法
3.1一次二阶矩法
利用高等数学中的泰勒公式将功能函数做泰勒展开并只保留展开后的一次项和二次项这种方法称为可靠度计算的一次二阶距法。
一次二阶距法早在1947年前苏联学者尼捏而然就提出了这种方法，后面有许多学者对这种方法做出了改进，这种方法精确度不是很高但是计算编程比较方便，在计算机不是很发达的时代采用较多，在这里只介绍本文用到的验算点法（JC法）。
3.1.1验算点法（JC法）
哈索弗尔（Hasofer）【14】和林德（Lind）、拉克维茨（Rackwitz）和菲斯莱（Fiessler）、帕洛赫摩（Paloheimo）和汉拉斯（Hannus【15】）等人提出了结构可靠度计算的验算点法。
验算点法对于非正态随机变量需将其近似为正态分布的线性随机变量，正态化的方法目前有三种：
（1）当量正态化，这是我国《建筑结构设计统一标准》和国际结构安全度联合会（（JCSS）推荐的方法（JC法）
（2）映射变换法，即采用数学变换的方法将非正态随机变量变换为正态随机变量。
（3）实用分析法，是近年山赵国藩【18】等提出的一种新的当量正态化方法，该方法是对JC法的改进。
验算点法是近十来年提出并在不断发展的结构可靠度分析方法，其基本思路与中心点法相仿，将非线性功能函数在验算点（在失效边界上）处作泰勒级数展开并保留至一次项，然后近似计算功能函数的平均值和标准差，再求可靠指标。
其特点是可以考虑非正态的随机变量，在计算量增加不多的情况下，可获得比中心点法更高的精度，且求得满足极限状态方程的“验算点”设计值，便于根据规范给出的标准值计算分项系数，以利于设计人员采用惯用
的多系数设计表达式。
设x*=（x*1x*2...，x*n）为极限状态面上一点，即满足Z=g（x*）=0，在x*处将极限状态方程Taylor展开并取至一次项得到ZL，并计算出ZL的均值LZμ和标准差σZL，结构的可靠度的表达式:
（7）
可以根据上式建立迭代公式，求解可靠度。
当基本变量X中含有非正态随机变量时，运用验算点法须事先处理非正态变量，这里用当量正态化法。
当量正态化条件要求在验算点x*i处Xi′和Xi的分布函数和概率密度函数分别对应相等，即：
（8）
（9）
验算点法的优点在于：
（1）它适用于随机变量为任意分布下结构可靠指标的求解，而且通俗易懂，计算速度快，计算精度又能满足工程的实际需要。
（2）它能给出一套固定的解题步骤，适合于编制计算程序和便于一般工程技术人员的应用[。
但其局限性在于：
（1）将极限状态方程在验算点处展为泰勒级数线性化极限状态方程，可能会带来显著性误差。
（2）由于将非正态变量等价正态化，也使计算带来误差。
（3）当在标准正态空间中的极限状态方程中有几个点到原点的距离取极值时，则问题的解将与初始迭代点有关，很可能得到的解是局部最优，而不是总体最优解。
3.2MonteCarlo抽样法
3.2.1直接蒙特卡罗法
Monte-Carlo法又称随机抽样技巧、概率模拟方法和统计试验法[15]。
其理论基础是概率论中的大数定律，具有模拟的收敛速度与基本随机向量的维数无关、极限状态函数的复杂程度与模拟过程无关、无需将状态函数线性化和随机变量当量正态化、能直接解决问题、数值模拟的误差可由模拟次数和精度较容易地加以确定的特点，因此其应用范围几乎没有什么限制。
利用电脑进行随机n（n比较大）次抽样的方法，将每次所抽取得结果带入到结构功能函数中，如果Z>0，g（x）=1，否则g（x）=0，这样就完成了一次计算，再产生下一随机数，重复上面的计算，直到完成预定的试验次数为止。
此时，失效概率为Pf=P（z<0）=
，式中，n是试验的总次数，k是试验中g（x）<0的次数，比值k/n是统计变量，对于低的失效概率或者n较小时，估算Pf值容易发生相当大的不确定性。
但是当模拟次数很大时，直至趋近于无穷大时，能够得出精确的Pf值；但是，当实际工程的结构破坏概率在10-8以下时，该法的模拟数日就会相当大，进而占用大量时间。
若基本随机变量相关时，利用条件概率密度，把多维问题化为一维问题来解决，其缺点就是它的计算中需耗费大量时间，获得的信息有限，只能得到失效概率或可靠指标，无法得到灵敏性系数，验算点信息。
因而在实际工程计算中较少采用，一般多用于检验一些新提出的计算方法的精度，或对几个方法进行比较。
3.2.2重要抽样蒙特卡洛法
蒙特卡洛法的模拟次数一般来说很可观，对于功能函数为隐式的结构系统，要进行成千上万次的有限元分析计算，效率很低。
重要抽样蒙特卡洛法其实于直接抽样额蒙特卡洛法没有本质上的区别，只是重要抽样蒙特卡洛法改变了抽样的“重心”，即将抽样的“重心”转移到对结构失效概率贡献最大的区域。
运用重要抽样蒙特卡洛法对随机变量进行抽样，在满足精度的条件下工作量减少提高了运算的效率。
设结构系统的基本变量矢量为X，其概率密度函数为
，则结构系统的失效概率：
其中，
，
表示第i个失效模式的安全裕量方程，m为系统所含失效模式的总数。
定义指示函数：
（10）
并引入重要抽样概率密度函数，V为重要抽样概率空间。
则
可以被写成
（11）
根据蒙特卡洛法的基本原理，失效概率的近似估计：
（12）
式中，N为仿真次数，由重要抽样概率密度函数抽样产生。
失效概率的方差：
（13）
由上式可见，失效概率估计值的方差主要依赖于重要抽样分布函数，也就是说，选择了合适的重要抽样密度函数，便可以有效减小方差，大大提高蒙特卡洛法的效率。
试验表明：
当结构精度为0.01时，重要抽样蒙特卡洛法要比直接抽样蒙特卡洛法抽样次数大大降低。
3.3一次渐近积分法
在一次二阶矩方法和二次二阶矩方法中，需要对基本随机变量进行当量正态化转换为正态分布的函数，而且对于非正态基本随机变量当量正态化造成一些误差，而且这种误差随着基本随机变量非线性程度的增加而提高。
因此学者们就提出了一种规避将基本随机变量当量正态化，直接计算结构可靠指标的方法，从而提高结构可靠指标的计算精度。
对失效概率的贡献主要是在结构失效最大可能点附近的积分，因此只要将积分局部化，集中在该点附近的失效区域内进行，就能够得到失效概率积分得近似结果。
失效概率的渐近积分是在失效最大可能点处，将基本变量概率密度函数的对数展开成Taylor级数并取至二次项，将功能函数也作Taylor级数展开，用所得超切平面或二次超曲面来逼近实际失效面，再利用一次二阶矩方法和二次二阶矩方法的成果即可完成失效概率的渐进积分。
在基本随机变量空间中用渐进积分方法计算结构的失效概率，无须变量空间的变换也不用到变量的累计分布函数，但要计算基本随机变量概率密度函数对数的一阶和二阶倒数，使处理问题的繁琐程度有所增加。
（1）求解
x*，根据解最优化问题的拉格朗日乘子法，求解x*。
（2）计算
。
（3）计算
。
（4）计算
。
一次渐进法避免将基本随机变量当量正态化，减小了由于基本随机变量当量正态化所带来的误差，使得计算精度有所提高，但是这导致的结果就是计算的复杂程度有所增加，目前用此方法进行结构可靠度计算还应用的较少。
3.4Breitung法
结构随机可靠度分析的一次二阶矩方法，概念比较简单，编程比较方便，因此在工程中有较广泛的应用，但是这种方法有因为只考虑到了功能函数的一阶倒数，因此，取得的关于功能函数的信息较少，得到的结果普遍差别较大，特别是当功能函数非线性程度较高时产生的误差更大。
功能函数的二阶倒数不仅考虑到了一阶倒数的信息，功能函数的二阶倒数还考可以运用功能函数在验算点附近的凹向、曲率等非线性性质，从而提高结构可靠度分析的精度。
目前在工程中应用较为广泛。
在功能函数极限状态曲面验算点附近计算结构可靠度，结构的非线性程度对于结构可靠度分析的影响很大，因此breitung法采用在验算点处将功能函数展开成泰勒级数并取至二次项，并以此二次函数代替原来的结构失效面；breitung法是以一次二阶距为方法为基础，先用一次二阶矩的方法取得可靠指标和失效概率，之后采用二次二阶矩方法对一次二阶矩方法取得的结构做出修正，取得的结果要比一次二阶矩方法取得的结果更精确一些。
Breitung法计算结构可靠指标分为两部分进行即：
首先采用一次二阶矩方法计算得到结构的可靠指标如本文上述的JC法的计算过程，之后运用二次二阶距方法计算结构的二阶倒数并计算功能函数的Hesse矩阵，并计算基本变量概率密度函数的导数。
本文给出breitung方法的计算步骤为：
计算β，采用一次二阶矩法；
计算各变量的偏导；
确定Hesse矩阵，采用正交规范化处理技术；
计算β，及失效概率；
Breitung法具有比一次二阶距方法更高的精度，在掌握了一次二阶距的编程方法之后，Breitung法先将一次二阶距方法程序写出来计算一次二阶距方法的可靠指标，再对Breitung方法的二次二阶倒数部分进行计算，从而获得二次二阶距方法的可靠指标。
4可靠度计算在Matlab环境下的实现与比较
算例1.1【16】一承载力为R的轴压短柱，承受荷载S作用。
已知R服从正态分布,；S服从对数正态分布,。
R,S相互独立。
试确定柱的受压承载能力的可靠指标。
柱的功能函数为Z=R-S。
利用MATLAB软件编写了JC法，一次渐进法，二次二阶矩法（Breitung），蒙特卡洛直接抽样法和蒙特卡洛重要法，具体程序在程序包中，计算结果如表1。
表1
项目
直接蒙特卡洛法
Breitung法
一次渐进法
JC法
重要抽样蒙特卡洛法
N=1e4
N=1e5
N=1e6
算例1.1
可靠指标（
）
1.8042
1.7872
1.6806
1.6666
1.8161
1.7995
1.8043
失效概率
0.0358
0.0359
0.0464
0.0478
0.0364
0.0359
0.0356
可靠指标相对于直接蒙特卡罗法误差
-0.94%
6.85%
7.63%
-0.66%
0.26%
0
算例1.2【16】一承载力为R的轴压短柱，承受荷载S作用。
已知R服从正态分布,；S服从对数正态分布,。
R,S相互独立。
试确定柱的受压承载能力的可靠指标。
柱的功能函数为Z=R-S。
利用MATLAB软件编写了JC法，一次渐进法，二次二阶矩法（Breitung），蒙特卡洛直接抽样法和蒙特卡洛重要法，具体程序在程序包中，计算结果如表2。
表2
项目
直接蒙特卡洛法
Breitung法
一次渐进法
JC法
重要抽样蒙特卡洛法
N=1e4
N=1e5
N=1e6
算例1.2
可靠指标（
）
3.6350
3.6430
3.6307
3.4057
3.7105
3.6455
3.6669
失效概率
（）
.3900
1.3069
1.4135
3.3001
1.0343
1.3341
1.2276
可靠指标相对于直接蒙特卡罗法误差
0
0.12%
6.31%
-2.08%
-0.29%
-0.88%
算例1.3【16】一承载力为R的轴压短柱，承受荷载S作用。
已知R服从正态分布,；S服从对数正态分布,。
R,S相互独立。
试确定柱的受压承载能力的可靠指标。
柱的功能函数为Z=R-S。
利用MATLAB软件编写了JC法，一次渐进法，二次二阶矩法（Breitung），蒙特卡洛直接抽样法和蒙特卡洛重要法，具体程序在程序包中，计算结果如表3。
表3
项目
直接蒙特卡洛法
Breitung法
一次渐进法
JC法
重要抽样蒙特卡洛法
N=1e4
N=1e5
N=1e6
算例1.3
可靠指标（
）
1.7669
1.7572
1.6504
1.6839
1.7403
1.7675
1.7663
失效概率
0.0386
0.0388
0.0494
0.0461
0.0409
0.0386
0.0387
可靠指标相对于直接蒙特卡罗法误差
0.55%
6.59%
4,.70%
1.51%
-0.03%
0.03%
算例1.4【15】一承载力为R的轴压短柱，承受荷载S作用。
已知R服从正态分布,；S服从对数正态分布,。
R,S相互独立。
试确定柱的受压承载能力的可靠指标。
柱的功能函数为Z=R-S。
利用MATLAB软件编写了JC法，一次渐进法，二次二阶矩法（Breitung），蒙特卡洛直接抽样法和蒙特卡洛重要法，具体程序在程序包中，计算结果如表4：
表4
项目
直接蒙特卡洛法
Breitung法
一次渐进法
JC法
重要抽样蒙特卡洛法
N=1e4
N=1e5
N=1e6
算例1.4
可靠指标（
）
3.5578
3.5594
3.5439
3.4144
3.4562
3.5704
3.5519
失效概率
（）
1.8700
1.8264
1.9712
3.1959
2.7390
1.7820
1.9124
可靠指标相对于直接蒙特卡罗法误差
-0.04%
0.39%
4.03%
2.86%
-0.35%
-0.35%
例题1设计为线性功能函数，设计方法采用功能函数和分布情况不发生变化，而参数发生变化。
分别通过均值不变标准差减小例1.2，均值增大标准差不变例1.3，均值增大标准差减小例1.4。
对四个例题分别进行对比讨论参数变化对于可靠指标精度的影响。
从四道例题总的来看JC法计算误差最大且结果偏大，重要抽样蒙特卡洛法相对误差最小同时重要抽样法抽样次数越多结果越精确，一次渐进法计算误差也较大，breitung法计算误差较小。
将例1.1同例1.2对比发现当均值不变标准差减小时breitung法和一次渐进法趋于精确，而重要抽样蒙特卡洛法精确度降低。
例1.2与1.3对比发现均值不变标准差减小与均值增大标准差不变所得结论相同。
将例1.1同例1.4对比发现当均值增大标准差减小时breitung法和一次渐进法精度较例1.1高。
由此可以得出结论，当功能函数和分布情况相同时，均值较大标准差较小的情况下采用breitung法和一次渐进法可以采用，重要抽样蒙特卡洛法可以采用抽样次数需大于1e6。
同时可以得到，重要抽样法所得结果随着抽样次数增加而精确。
算例2.1【16】已知结构的功能函数为，其中X1服从对数正态分布，X2服从极值1型分布，X3服从Weibull分布；其均值和标准差分别为;;。
试计算结构的可靠指标和失效概率。
利用MATLAB软件编写了JC法，一次渐进法，二次二阶矩法（Breitung），蒙特卡洛直接抽样法和蒙特卡洛重要法，具体程序在程序包中，计算结果如表5：
表5
项目
直接蒙特卡洛法
Breitung法
一次渐进法
JC法
重要抽样蒙特卡洛法
N=1e4
N
升级会员 