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

可靠性算法，改进一次二阶矩法的Matlab源代码，结合Rackwitz-Fiessler方法，能够考虑任意分布的随机变量，里面包含部分测试例子，可直接在Matlab软件中调用执行，文件中包含详细的注释。

function AFOSMJC（NE）
%AFOSM法+JC法（Rackwitz-Fiessler方法）
%适用范围：随机变量为任意分布，变量间独立
%特点：保证等价前后失效概率具有较高的近似精度，特征点取MPP点
%缺点：设计点迭代可能不收敛，可能陷入局部最优，可使用遗传算法计算设计点
%输入参数：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.获得相应参数
    gFun=Prob.Fung;                                         %用符号变量表示极限状态函数
    for i=1:Prob.Nx   
        [muX（i:）sigmaX（i:）]=RandomVariableStat（Probi）;
    end; 
    
    %4.推导偏导函数，符号变量赋值
    for i=1:Prob.Nx
        dgFun（i:）=diff（gFunX（i））;                         %求极限状态函数偏导数表达式
    	eval（[XX（i:） ‘=muX（i:）;‘]）;                       %为符号变量赋值 - 变量均值
    end;
    
    %5.计算功能函数值及导数值
    gValue=eval（gFun）;                                      %求解极限状态函数在均值点处函数值
    for i=1:Prob.Nx                                        
        dgValue（i:）=eval（subs（dgFun（i）））;                  %求解极限状态函数在均值点处导数值
    end;

    %6.迭代求解MPP点
    x=muX;
    x0=repmat（epslength（muX）1）;
    IteratingStep=0;
    while norm（x-x0）/norm（x0） > 1e-6
        IteratingStep=IteratingStep+1;
        x0=x;
        %6.1符号变量赋值
        for i=1:Prob.Nx
            eval（[XX（i:） ‘=x（i）;‘]）;                       %为符号变量赋值 - xi迭代值
        end;
        %6.2计算功能函数值及导数值
        gValue=eval（gFun）;
        dgValue=eval（subs（dgFun））;                          %求在xi点的偏导数值dg/dxi|xi
        %%%%%%%%%%%%%%%%%%%%
        %6.3 JC法对非正态分布变量进行转换处理
        [muXsigmaX]=JCMethod（Prob x）;
        %%%%%%%%%%%%%%%%%%%%
        %6.4计算加权标准差
        gs=dgValue.*sigmaX;
        %如果gs值为零，需要偏移当前点位置，重新计算
        if（norm（gs）） == 0
            x=x+0.1;
            continue;
        end
        %6.5计算灵敏度系数
        alphaX=-gs/norm（gs）;
        %6.6计算可靠度指标
        beta=（gValue+dgValue.‘*（muX-x））/norm（gs）;
        %6.7计算新的x向量值
        x=muX+sigmaX.*alphaX*beta;
    end
    
    %7.求解函数值及失效概率
    MPP=eval（subs（x））;                                       %求解基本变量x的值
    Pf=normcdf（-beta）;                                       %失效概率
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fprintf（‘××××××AFOSMJC 失效模式可靠性分析结果××××××\n‘）; 
    fprintf（‘设计点迭代步数：    %d\n‘IteratingStep）; 
    fprintf（‘设计点坐标：‘）;
    for i=1:Prob.Nx
        fprintf（‘        %10.8f‘MPP（i））;
    end    
    fprintf（‘\n‘）; 
    fprintf（‘设计点g值：         g（x）=%10.10f\n‘gValue）;
    fprintf（‘可靠度指标：        β=%10.10f\n‘beta）;
    fprintf（‘失效概率：          Pf=%10.10f\n‘Pf）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    toc
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%