资源简介
可靠性算法,改进一次二阶矩法的Matlab源代码,里面包含部分测试例子,可直接在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
评论
共有 条评论