大家好，又见面了，我是你们的朋友全栈君。
题目：编制改进一次二阶矩法计算可靠指标的程序，并给出算例，要求提供源程序，算法语言不限。
选取的算例为：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_r(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_r(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
发布者：全栈程序员-用户IM，转载请注明出处：https://javaforall.cn/142942.html原文链接：https://javaforall.cn
【正版授权，激活自己账号】： Jetbrains全家桶Ide使用，1年售后保障，每天仅需1毛
【官方授权 正版激活】： 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...
![[转]组合数取模【转自AekdyCoin的组合数取模】](http://javaforall.cn/wp-content/themes/justnews/themer/assets/images/lazy.png)