Matlab中求解非连续型函数最优解
问题描述:
有一组装件由七个零件组成,这七个零件参数服从正态分布
x1~N(0.1,(0.005/3) ^2);
x2~N(0.3,0.0052^2);
x3~N(0.1,(0.005/3)^2);
x4~N(0.1,0.0052^2);
x5~N(1.5,(0.225/3)^2);
x6~N(16,(0.8/3)^2);
x7~N(0.75,(0.0375/3)^2);
该组装件的性能函数:
目标值是y0 = 1.5,偏离±0.1时,产品为次品,损失1000元;偏离±0.3时为废品,损失9000元。
(1)计算生产一件该产品的损失期望值;
(2)以各零件参数的均值u为变量,波动范围为20%,方差为固定值,试确定各零件的均值来优化产品性能。
1*第一问求解:
用normrnd创建七个服从正态分布的随机变量,对y进行大量实验,统计其落在某区间的概率,以求得损失期望值。
syms k j n m
k=0;
n=0;
m=0;
for i=1:1000
x1 = normrnd(0.1,(0.005/3));
x2 = normrnd(0.3,0.005);
x3 = normrnd(0.1,(0.005/3));
x4 = normrnd(0.1,0.005);
x5 = normrnd(1.5,(0.225/3));
x6 = normrnd(16,(0.8/3));
x7 = normrnd(0.75,(0.0375/3));
y = 174.42*(x1/x5)*(x3/(x2-x1))^0.85*sqrt((1-2.62*(1-0.36*(x4/x2)^(-0.56))^(3/2)*(x4/x2)^1.16)/(x6*x7));
if(1.4<y&&y<1.6)
k=k+1; %统计合格品个数
elseif(((1.2<y&&y<1.4)||y==1.4)||((1.6<y&&y<1.8)||y==1.6))
n=n+1; %统计次品个数
elseif(y>1.8||y==1.8||y<1.2||y==1.2)
m=m+1; %统计废品个数
end
end
k
n
m
j = (n*1000+m*9000)/1000 %损失
2*第二问求解
分别用fmincon和ga遗传算法求解,二者均可以设置所求变量的上下限和其满足的约束函数;二者均需另外建立方法的代码文件,再在另一代码文件中调用;前者要给出变量初值,后者不用。
用fmincon求解
#fun.m#代码文件
function y0 = fun(x)
y1 = 'abs(174.42*(normrnd(x(1),(0.005/3))/normrnd(x(5),(0.225/3)))*(normrnd(x(3),(0.005/3))/(normrnd(x(2),0.005)-normrnd(x(1),(0.005/3))))^0.85*sqrt((1-2.62*(1-0.36*(normrnd(x(4),0.005)/normrnd(x(2),0.005))^(-0.56))^(3/2)*(normrnd(x(4),0.005)/normrnd(x(2),0.005))^1.16)/(normrnd(x(6),(0.8/3))*normrnd(x(7),(0.0375/3))))-1.5)';
y2 = '+abs(174.42*(normrnd(x(1),(0.005/3))/normrnd(x(5),(0.225/3)))*(normrnd(x(3),(0.005/3))/(normrnd(x(2),0.005)-normrnd(x(1),(0.005/3))))^0.85*sqrt((1-2.62*(1-0.36*(normrnd(x(4),0.005)/normrnd(x(2),0.005))^(-0.56))^(3/2)*(normrnd(x(4),0.005)/normrnd(x(2),0.005))^1.16)/(normrnd(x(6),(0.8/3))*normrnd(x(7),(0.0375/3))))-1.5)';
for i = 1:100
y1 = strcat(y1,y2); %前面y1,y2定义为字符串,用strcat将y1,y2拼接起来,目的是使其求出的解更接近最优解
end
y0 = eval(y1); %eval将字符串转化为命令
end
#fmincon_method#代码文件
syms k j n m
k=0;
n=0;
m=0;
x0 = [0.1;0.3;0.1;0.1;1.5;16;0.75]; %初值
lb = [0.08;0.24;0.08;0.08;1.2;12.8;0.6]; %变量下限
ub = [0.12;0.36;0.12;0.12;1.8;19.2;0.9]; %变量上限
fval = 150;
while (fval>100)
[x,fval] = fmincon(@fun,x0,[],[],[],[],lb,ub); %其中四个[]表示等式约束和不等式约束均为空,x是所求变量的数列,fval是#fun.m#中所求函数的“最优值”(可能只是局部最优,甚至因为变量本身有随机性而产生的错误结果)
end
x
fval
for i=1:1000
x1 = normrnd(x(1),0.005/3);
x2 = normrnd(x(2),0.005);
x3 = normrnd(x(3),0.005/3);
x4 = normrnd(x(4),0.005);
x5 = normrnd(x(5),0.225/3);
x6 = normrnd(x(6),0.8/3);
x7 = normrnd(x(7),0.0375/3);
y = 174.42*(x1/x5)*(x3/(x2-x1))^0.85*sqrt((1-2.62*(1-0.36*(x4/x2)^(-0.56))^(3/2)*(x4/x2)^1.16)/(x6*x7));
if(1.4<y&y<1.6)
k=k+1; %统计合格品个数
elseif(((1.2<y&y<1.4)|y==1.4)|((1.6<y&y<1.8)|y==1.6))
n=n+1; %统计次品个数
elseif(y>1.8|y==1.8|y<1.2|y==1.2)
m=m+1; %统计废品个数
end
end
k
n
m
j = (n*1000+m*9000)/1000 %损失
用ga遗传算法求解
#ga_fun.m#代码文件
function f = ga_fun(x)
y1 = 'abs(174.42*(normrnd(x(1),(0.005/3))/normrnd(x(5),(0.225/3)))*(normrnd(x(3),(0.005/3))/(normrnd(x(2),0.005)-normrnd(x(1),(0.005/3))))^0.85*sqrt((1-2.62*(1-0.36*(normrnd(x(4),0.005)/normrnd(x(2),0.005))^(-0.56))^(3/2)*(normrnd(x(4),0.005)/normrnd(x(2),0.005))^1.16)/(normrnd(x(6),(0.8/3))*normrnd(x(7),(0.0375/3))))-1.5)';
y2 = '+abs(174.42*(normrnd(x(1),(0.005/3))/normrnd(x(5),(0.225/3)))*(normrnd(x(3),(0.005/3))/(normrnd(x(2),0.005)-normrnd(x(1),(0.005/3))))^0.85*sqrt((1-2.62*(1-0.36*(normrnd(x(4),0.005)/normrnd(x(2),0.005))^(-0.56))^(3/2)*(normrnd(x(4),0.005)/normrnd(x(2),0.005))^1.16)/(normrnd(x(6),(0.8/3))*normrnd(x(7),(0.0375/3))))-1.5)';
for i = 1:13
y1 = strcat(y1,y2);
end
f = eval(y1);
end
#ga_method#代码文件
syms k j n m
k=0;
n=0;
m=0;
lb = [0.08;0.24;0.08;0.08;1.2;12.8;0.6];
ub = [0.12;0.36;0.12;0.12;1.8;19.2;0.9];
[x,fval,exitflag] = ga(@ga_fun,7,[],[],[],[],lb,ub); %7表示有七个变量
x
fval
for i=1:1000
x1 = normrnd(x(1),0.005/3);
x2 = normrnd(x(2),0.005);
x3 = normrnd(x(3),0.005/3);
x4 = normrnd(x(4),0.005);
x5 = normrnd(x(5),0.225/3);
x6 = normrnd(x(6),0.8/3);
x7 = normrnd(x(7),0.0375/3);
y = 174.42*(x1/x5)*(x3/(x2-x1))^0.85*sqrt((1-2.62*(1-0.36*(x4/x2)^(-0.56))^(3/2)*(x4/x2)^1.16)/(x6*x7));
if(1.4<y&y<1.6)
k=k+1; %统计合格品个数
elseif(((1.2<y&y<1.4)|y==1.4)|((1.6<y&y<1.8)|y==1.6))
n=n+1; %统计次品个数
elseif(y>1.8|y==1.8|y<1.2|y==1.2)
m=m+1; %统计废品个数
end
end
k
n
m
j = (n*1000+m*9000)/1000 %损失
总结,用fmincon求解有一定随机性,想获得良好的优化参数有一定的运气成分,用ga求解比较稳定,一般能取得较好的优化参数。