可靠性方法：改进一次二阶矩法的Matlab源代码

function AFOSM（NE）
%改进一次二阶矩法
%适用范围：随机变量为正态分布，变量间独立
%特点：线性功能函数获得失效概率精确解，非线性程度不高的功能函数获得失效概率较高精度的近似解
%特点：对于物理意义相同而数学表达是不同的问题具有统一的解
%缺点：不能反映功能函数的非线性对失效概率的影响
%缺点：在功能函数的非线性程度较大时，迭代算法受初始点影响较大
%缺点：对具有多个设计点问题，可能会陷入局部最优，甚至不收敛
%缺点：属基于功能函数梯度的方法，对极限状态方程表达式具有一定的依赖性，隐式函数梯度求解较困难
%缺点：设计点迭代可能不收敛，可能陷入局部最优，可使用遗传算法计算设计点
%输入参数：NE - 极限状态函数索引，整数
    global Prob
    tic

    if nargin<1 NE = 0; end                                 %如果未定义方程索引，默认使用索引0
    %1.定义符号变量
    syms x1  x2  x3  x4  x5  x6 x7  x8  x9  xa  xb  xc  xd  xe  xf  xg;
    X=[ x1  x2  x3  x4  x5  x6 x7  x8  x9  xa  xb  xc  xd  xe  xf  xg ];
    XX=[‘x1‘;  ‘x2‘;  ‘x3‘  ;‘x4‘ ; ‘x5‘;  ‘x6‘ ;‘x7‘;  ‘x8‘;  ‘x9‘;  ‘xa‘;  ‘xb‘;  ‘xc‘;  ‘xd‘;  ‘xe‘;  ‘xf‘; ‘xg‘];

    %2.获得方程相关公式、变量数目、变量分布、变量参数
    Prob = FunDist（ NE ）;

    %3.判断方法适用性
    for i=1:Prob.Nx
        if strcmp（Prob.Dist{i} ‘norm‘）==0 && strcmp（Prob.Dist{i} ‘Normal‘）==0
        	error（‘This method can only be used for normal distributed variables please use AFOSMJC method instead\n‘）;
    	end;
    end;
    
    %4.获得相应参数
    gFun=Prob.Fung;                                         %用符号变量表示极限状态函数
    for i=1:Prob.Nx                                         
        muX（i:）=Prob.Para{i}（1）;
        sigmaX（i:）=Prob.Para{i}（2）;
    end;

    %5.推导偏导函数，符号变量赋值
    for i=1:Prob.Nx
        dgFun（i:）=diff（gFunX（i））;                         %求极限状态函数偏导数表达式
    	eval（[XX（i:） ‘=muX（i:）;‘]）;                       %为符号变量赋值 - 变量均值
    end;
    
    %6.计算功能函数值及导数值
    gValue=eval（gFun）;                                      %求解极限状态函数在均值点处函数值
    for i=1:Prob.Nx                                        
        dgValue（i:）=eval（subs（dgFun（i）））;                  %求解极限状态函数在均值点处导数值
    end;
    
    %7.迭代求解MPP点
    x=muX;
    x0=repmat（epslength（muX）1）;
    IteratingStep=0;
    while norm（x-x0）/norm（x0） > 1e-6
        IteratingStep=IteratingStep+1;
        x0=x;
        %7.1符号变量赋值
        for i=1:Prob.Nx
            eval（[XX（i:） ‘=x（i）;‘]）;                       %为符号变量赋值 - xi迭代值
        end;
        %7.2计算功能函数值及导数值
        gValue=eval（gFun）;
        dgValue=eval（subs（dgFun））;                          %求在xi点的偏导数值dg/dxi|xi
        %7.3计算加权标准差
        gs=dgValue.*sigmaX;
        %如果gs值为零，需要偏移当前点位置，重新计算
        if（norm（gs）） == 0
            x=x+0.1;
            continue;
        end
        %7.4计算灵敏度系数
        alphaX=-gs/norm（gs）;
        %7.5计算可靠度指标
        beta=（gValue+dgValue.‘*（muX-x））/norm（gs）;
        %7.6计算新的x向量值
        x=muX+sigmaX.*alphaX*beta;
    end

    %8.求解函数值及失效概率
    MPP=eval（subs（x））;                                       %求解基本变量x的值
    Pf=normcdf（-beta）;                                       %失效概率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fprintf（‘××××××AFOSM 失效模式可靠性分析结果××××××\n‘）; 
    fprintf（‘设计点迭代步数为：  %d\n‘IteratingStep）; 
    fprintf（‘设计点坐标：‘）;
    for i=1:Prob.Nx
        fprintf（‘        %10.8f‘MPP（i））;
    end