一次二阶矩法计算可靠指标的程序

题目：编制改进一次二阶矩法计算可靠指标的程序，并给出算例，要求提供源程序，算法语言不限。

选取的算例为：z=g(x,y)=x*y-1140，其中x,y服从正态分布，μx=38,Vx=0.1, μy=38,Vy=0.05 本程序采用Matlab编写。 选取β1=3.0，β2=2.5

计算结果为：可靠指标β＝4.2672，最终验算点为：(22.8430 , 49.9060)，在验算点处功能函数值为：1.2354e-004

%保存为strRlbt.m，在Matlab命令窗口中输入strRlbt执行即可 N = 2;%变量个数 miu = [38 54];%均值 v = [0.1 0.05];%变异系数 sgma = miu .* v;%方差 syms x y

g = sym('x * y - 1140');%功能函数 jacg = jacobian( g ,[x;y]);%计算雅可比矩阵 initvalue = [miu;v;sgma];%用作函数参数 %选取beta，定义x0=miu beta1 = 3.0; xopt0 = [38 54];

alpha0 = zeros(1,2);

[ alpha1 , xopt1 , result ] = calforbeta( initvalue , beta1 , alpha0 , xopt0 , jacg , g ); if result == 1

disp('第一次假定的饧次煽恐副辏'); return end

%再次假定beta beta2 = 2.5;

xopt0 = miu - beta2 * alpha1 .* sgma; gvalue = jacgfunc(jacg,xopt1);

alpha0 = (sgma .* gvalue) / sqrt(sum((sgma .* gvalue).^2));

[ alpha2 , xopt2 , result ] = calforbeta( initvalue , beta2 , alpha0 , if result == 1

disp('第二次假定的饧次煽恐副辏'); return end

?ta迭代求解 g1 = gfunc(g,xopt1);

xopt0 , jacg , g ); g2 = gfunc(g,xopt2); eps = 0.1; %精度 while abs(g2) > eps temp = beta2;

beta2 = beta2 - (beta2 - beta1)/(g2 - g1) * g2; beta1 = temp;

[ alpha2 , xopt2 , result ] = calforbeta( initvalue , beta2 , alpha1 , temp = g2; g2 = gfunc(g,xopt2); g1 = temp; if result == 1 break end end

disp('可靠指标为：'); disp(beta2);

disp('最终验算点为：'); disp(xopt2);

disp('在验算点处功能函数值为：'); disp(g2);

xopt1 , jacg , g );

function g_out = gfunc( g , x_in ) %功能函数值计算 x = x_in(1); y = x_in(2);

g_out = eval(g);%函数值 %将以上内容保存为gfunc.m

function g_out = jacgfunc( jacg , x_in ) %功能函数偏导数计算，即雅可比矩阵计算 x = x_in(1); y = x_in(2); for i = 1:2

g_out(i) = eval(jacg(i));%1为对x的导数,2为对y的导数 end

%将以上内容保存为jacgfunc.m

function [ alpha1 , xopt1 ,result ] = calforbeta( initvalue , beta0 , alpha0 , xopt0 , jacg , g) %对选取的beta进行计算 result = 0;

N = length(xopt0); alpha = alpha0; xopt = xopt0; %initvalue为初始值

miu = initvalue(1,:);%第一行为均值 v = initvalue(2,:);%第二行为变异系数 sgma = initvalue(3,:);%第三行为方差 eps = 0.1; while 1

%功能函数达到精度则退出循环,result=1表示计算出可靠指标 if abs(gfunc(g,xopt0)) < eps alpha1 = alpha0; xopt1 = xopt0; result = 1; break; end

%计算alpha和新的验算点xopt gvalue = jacgfunc(jacg,xopt);

sgmaz = sqrt(sum((sgma .* gvalue).^2)); alpha0 = sgma .* gvalue / sgmaz;

xopt0 = miu - beta0 * alpha0 .* sgma; sum1 = sum((alpha - alpha0).^2); sum2 = sum((xopt - xopt0).^2); alpha = alpha0; xopt = xopt0;

%醝和验算点xi达到精度则退出循环 if sum1 < 0.001 | sum2 < 0.001 alpha1 = alpha0; xopt1 = xopt0; break; end end

%将以上内容保存为calforbeta.m