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

选取的算例为：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 ,  xopt0 , jacg , g );

if result == 1

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

   return

end

%beta迭代求解

g1 = gfunc(g,xopt1);

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 ,  xopt1 , jacg , g );

   temp = g2;

   g2 = gfunc(g,xopt2);

   g1 = temp;

   if result == 1

      break

   end

end

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

disp(beta2);

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

disp(xopt2);

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

disp(g2);

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