Genin.一日一进步

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求解比较稳定，一般能取得较好的优化参数。