基于matlab环境下蒙特卡罗法的实现

基于Matlab环境下蒙特卡罗法的实现针对应用蒙特卡罗对连续型分布采取直接抽样法解决结构可靠度所遇到的困难，提出利用MATLAB其强大数值计算功能来解决此类问题。利用MATLAB进行蒙特卡罗抽样模拟，在一定程度上减少了对连续型分布采用直接抽样时的困难，大大提高了计算效率。1.蒙特卡罗法蒙特卡洛方法是以数理统计原理为基础的，又称随机模拟方法，是随着电子计算机的发展而逐步发展起不来的一种独特的数值方法。用蒙特卡洛方法来研究事件的随机性是结构可靠度分析的一个重要方面。蒙特卡洛方法的优点是，它回避了结构可靠度分析中的数学困难，不需要考虑结构极限状态曲面的复杂性，只需要得到结构的响应即可；缺点是计算虽大，因此目前还不作为一种常规的结构可靠度分析的方法来使用，只适用于一些情况复杂的结构，由于其具有相对较高的精度，常用于结构可靠度各种近似方法计算精度的检验和计算结果的校核。直接抽样方法是蒙特卡洛分析最基本的一种方法，对于基本随机变量12(,,,)nXXXX，其概率密度函数为()fx，对应结构某一状态的功能函数为()Zgx。将随机样本值序列X代入功能函数()Zgx，若Z0，则模拟的结构失效一次。若总的模拟数为N，功能函数Z0的次数为fn，则结构失效概率fP的估计值ˆfP为：ˆffnPN(1.1)由伯努利大数定理：lim()1ffNnPPN(1.2)可得ˆfP以概率收敛于fP。失效概率的同样可以表达为：[()]()fPIgxfxdx(1.3)其中[()]Igx为()gx的示性函数，即：1()0[()]0()0gxIgxgx(1.4)则结构失效概率fP的估计值ˆfP为：11ˆ[()]NffiinPIgxNN(1.5)对于结构可靠度问题，其对应的结构失效概率的数量级通常为371010。由上式可以看出，这类小概率事件的结构失效问题，直接抽样法模拟一次()0Zgx的机会是非常不容易的。针对可靠性问题，提高蒙特卡洛方法的提高抽样效率的途径是增加Z0的次数，使得抽取的样本点有较多的机会落入失效域内，这时需要改变抽样的“重心”，这就是重要抽样的基本思想。2.最优化原理的蒙特卡罗法本文运用的最优化原理的蒙特卡罗法来实现程序的编写。最优化法的基本原理可以概括为：从一般情况出发，设有n个呈任意分布的独立随机变量影响结构的可靠度，极限状态函数为12(,,,)0nZXXXX。先将非正态变量当量正态化，其正态化原则可以用R-F法(拉科维茨-菲斯莱法)，得到可靠度指标，即2**2[()/]iixiXXu(2.1)可靠度指标β与失效概率ˆfP存在一一对应的关系，β越小，ˆfP越大，因此β是可以用来衡量结构稳定性的一个指标。而且考虑到各变量的中值与变异系数，β与K也存在线性关系，β越大，K越大。开始时设计验算点是未知的，只有通过β的优化求解才能确定，即把β看成是极限状态曲面上点12(,,,)nPXXX的函数，通过寻找β的最小值，可得β和设计验算点，即可建立优化数学模型2*2min[()/]ixixiXu(2.2)12.(,,,)0nstZgXXX(2.3)对于一般工程可靠度问题，由极限状态函数可将其中一个变量由其他变量表示1211(,,,,,,)iiinXgXXXXX，于是可转化成无约束优化模型2*2212111,min[()/](,,,,,,)nixixiiiniijXugXXXXX(2.4)在结构可靠度的数值模拟中，蒙特卡罗法模拟的收敛速度与随机变量的维数无关，极限状态函数的复杂程度与模拟过程无关，数值模拟的误差也可轻易地确定，从而确定模拟的次数和精度。但是，常用的蒙特卡罗法只能计算出结构失效概率，并以此换算出可靠度指标，无法得到验算点。最优化原理的蒙特卡罗法则可算出结构失效概率、可靠度指标和验算点。3.程序设计思路在最优化的蒙特卡洛法下，根据可靠度指标的几何意义，在标准正态坐标系中，可靠指标是原点到极限状态曲面的最短距离，而验算点为原点到极限状态曲面距离最短的点。根据此几何意义，可采用蒙特卡罗法计算验算点及可靠指标。求解这一最优化问题，可通过对已知分布的随机变量在极限状态曲面上进行抽样，计算各抽样点与原点距离，最小距离即为所求。计算步骤如下：（1）对每个随机变量1211,,,,,,iinXXXXX均产生M个具有已知分布的随机数；（2）由极限状态方程计算随机变量1211(,,,,,,)iiinXgXXXXX；（3）作映射变换*1[()]iiiYFX，其中*()iiFX为随机变量iX的概率分布系数；（4）计算21/21[],12,njjiiYjn，，。其中minj就是所要求的可靠度指标；（5）根据可靠度指标β与失效概率ˆfP的函数关系式，可以算出结构的可靠度为1()fP，()x为标准正态分布函数。4.计算案例已知非线性极限状态方程312312(,,)5670XZgXXXXXe，其中123,,XXX均服从正态分布，11222.18,0.0564;2.18,0.0564;XXXX332.18,0.0564XX。并且123,,XXX相互独立，求可靠度指标，可靠度fP及验算点***123,,XXX。（1）程序代码及说明（案例程序是基于Matlab7.0环境下编写）functionmengte(x)fprintf('抽样次数');n=input('n=');x1=normrnd(2.18,0.0654,1,n);x2=normrnd(0.6,0.078,1,n);x=(1134*(x1.*x2)).^0.5;x3=log(x);h1=normcdf(x1,2.18,0.0654);a=norminv(h1,0,1);h2=normcdf(x2,0.6,0.0786);b=norminv(h2,0,1);h3=normcdf(x3,3.49,0.03);c=norminv(h3,0,1);zb=sqrt(a.^2+b.^2+c.^2);B=min(zb);y=find(zb==B);k=1-normcdf(-B);fprintf('可靠度指标%f\n',B);fprintf('可靠度%f\n',k);fprintf('验算点为[%f,%f,%f]\n',x1(y),x2(y),x3(y));fprintf('-------------------\n');%定义带参变量的主函数%输入要进行抽样的次数%生成服从正态分布的随机变量x1%生成服从正态分布的随机变量x2%根据极限状态方程求出随机变量x3%将随机变量x1,x2,x3分别映射变换成a,b,c%进行高斯积分，以x1为分界，坐标轴左边的所有分布的累计概率%求满h1下的反函数值%计算n次抽样的可靠度指标%选择n次抽样中最小可靠度指标%定位于最小可靠度指标%计算结构的可靠度%输出求得的靠度指标%输出求得的靠度%输出求得的验算点（2）程序运行结果由下图可以看到，当抽样次数达到104时，其可靠度指标，可靠度及验算点已经趋于稳定，并且程序在抽样次数为106情况下运行时间不超过2秒钟。5.小结蒙特卡罗方法具有“精度随所抽样次数增大而增大”的优点，但是却存在“运算量随抽样次数增大而急剧增加”的缺点，而通过Matlab统计工具箱，利用Matlab强大的矩阵运算功能，实现了蒙特卡罗方法计算。并且Matlab编写程序所用语句比其它语言更少，效率更高。