Matlab呼叫遺傳工具箱復現論文模型求解部分
原文轉載自：https://blog.csdn.net/robert_chen1988/article/details/52431594
論文來源：
https://www.sciencedirect.com/science/article/pii/S0045782599003898
function [xm,fv]=GAEOQ1
%%初始目標函式與約束條件
%求解變數t, k, 目標函式 f,約束條件 g1
syms t k z;
f=100/t+25*(25*t+k*10*sqrt(1+1.25*t))+100*10*sqrt(1+1.25*t)*int((z-k)*normpdf(z,0,1),z,k,inf)/t;
tmin=1e-2;
tmax=4;
kmin=0;
kmax=2;
g1=1-200*sqrt(2+0.75*t)*int((z-k)*normpdf(z,0,1),z,k,inf)/(50*t);%%將約束標準化，將右端變為1
NP=50; %進化多少代
%%初始化種群，種群長度40
size=20;
E=zeros(20,5); %前兩列為初始解，第三列為適應函式值，第四列記錄是否為可行解，第五列記錄違背約束條件的差值
E(:,1)=tmin+(tmax-tmin)*rand(size,1);
E(:,2)=kmin+(kmax-kmin)*rand(size,1);
fv=inf;%初始最優值為無窮大的值
D=zeros(NP,4);%用來記錄每代的最優解,平均值，最差解,最優解是否為可行解
%%計算適應函式罰函式值，判斷是否為可行解
for i=1:size
B=zeros(1,1);
B(1)=subs(g1,[t,k],E(i,(1:2)));
if B(1)>=0
E(i,4)=1;
E(i,3)=subs(f,[t,k],E(i,(1:2)));
else
E(i,4)=0;
E(i,3)=0;
end
if B(1)>=0
B(1)=0;
else
B(1)=abs(B(1));
end
E(i,5)=B(1);
end
fmax=max(E(:,3));
for i=1:size
if E(i,4)<1e-6
E(i,3)=fmax+E(i,5);
end
end
%%遺傳進化 %%到這步適應值還沒出錯
for g=1:NP %%原來錯誤在這裡，這個k跟前面的k重複了
%%競標賽選擇 %%小生態技術
M=zeros(size,2);%用來儲存優勝者的中間矩陣
for i=1:size
%A=randperm(size,6);
%dij1=sqrt(0.5*(E(A(1),1)-E(A(2),1))^2); %%小生態技術,只有單界時小生態技術沒法用
%dij2=sqrt(0.5*(E(A(1),1)-E(A(3),1))^2);
%dij3=sqrt(0.5*(E(A(1),1)-E(A(4),1))^2);
%dij4=sqrt(0.5*(E(A(1),1)-E(A(5),1))^2);
%dij5=sqrt(0.5*(E(A(1),1)-E(A(6),1))^2);
%if dij1<0.1
%if E(A(1),3)<=E(A(2),3)
%M(i,:)=E(A(1),(1:2));
%else
%M(i,:)=E(A(2),(1:2));
%end
%continue;
%elseif dij2<0.1
%if E(A(1),3)<=E(A(3),3)
%M(i,:)=E(A(1),(1:2));
%else
%M(i,:)=E(A(3),(1:2));
%end
%continue;
%elseif dij3<0.1
%if E(A(1),3)<=E(A(4),3)
%M(i,:)=E(A(1),(1:2));
%else
%M(i,:)=E(A(4),(1:2));
%end
%continue;
%elseif dij4<0.1
%if E(A(1),3)<=E(A(5),3)
%M(i,:)=E(A(1),(1:2));
%else
%M(i,:)=E(A(5),(1:2));
%end
%continue;
%elseif dij5<0.1
%if E(A(1),3)<=E(A(6),3)
%M(i,:)=E(A(1),(1:2));
%else
%M(i,:)=E(A(6),(1:2));
%end
%else
%M(i,:)=E(A(1),(1:2));
%end
%end
A=randperm(size,2);
if E(A(1),3)<=E(A(2),3)
M(i,:)=E(A(1),(1:2));
else
M(i,:)=E(A(2),(1:2));
end
end
%%模擬二進位制交叉生成後代
for j=1:size/2
if rand()>=0.5
A=randperm(size,2);
c=rand();
x2=max(M(A(1),1),M(A(2),1));
x1=min(M(A(1),1),M(A(2),1));
beita1_t=1+2*(x1-tmin)/(x2-x1);
rfa_t=2-beita1_t^(-2);
if c<=1/rfa_t
beita2_t=sqrt(rfa_t*c);
else
beita2_t=sqrt(1/(2-rfa_t*c));
end
E(2*j-1,1)=0.5*(x1+x2-beita2_t*(x2-x1));
E(2*j,1)=0.5*(x1+x2+beita2_t*(x2-x1));
end
%%只在可行解時出錯是怎麼回事？是不是變異的原因，已糾正
if rand()>0.5
c=rand();
x2=max(M(A(1),2),M(A(2),2));
x1=min(M(A(1),2),M(A(2),2));
beita1_t=1+2*(x1-kmin)/(x2-x1);
rfa_t=2-beita1_t^(-2);
if c<=1/rfa_t
beita2_t=sqrt(rfa_t*c);
else
beita2_t=sqrt(1/(2-rfa_t*c));
end
E(2*j-1,2)=0.5*(x1+x2-beita2_t*(x2-x1));
E(2*j,2)=0.5*(x1+x2+beita2_t*(x2-x1));
end
end
%%變異，變異會不會導致可行解不可行?單下界時不用變異,設定判斷條件防止過界
for i=1:size
nita=100+g;
pm=1/size+g*(1-1/size)/NP;
if rand()<pm
u=rand();
%x=E(i,1);
x=E(i,1);%%
deltamax=1;
if u<=0.5
delta_2=(2*u)^(1/(nita+1))-1;
else
delta_2=1-(2*(1-u))^(1/(nita+1));
end
if x+delta_2*deltamax>=tmin
E(i,1)=x+delta_2*deltamax;
end
end
if rand()<pm
u=rand();
x=E(i,2);
deltamax=1;
if u<=0.5
delta_1=(2*u)^(1/(nita+1))-1;
else
delta_1=1-(2*(1-u))^(1/(nita+1));
end
if x+delta_1*deltamax>=kmin
E(i,2)=x+delta_1*deltamax;
end
end
end
%%計運算元代罰函式值，判斷是否滿足可行解
for i=1:size
B(1)=subs(g1,[t,k],E(i,(1:2)));
if B(1)>=0
E(i,4)=1;
E(i,3)=subs(f,[t,k],E(i,(1:2)));%%跟直接算的結果不一樣，也跟EOQ得到的結果不一樣 eval出錯的原因
else
E(i,4)=0;
E(i,3)=0;
end
if B(1)>=0
B(1)=0;
else
B(1)=abs(B(1));
end
E(i,5)=B(1);
end
for i=1:size
if E(i,4)<1e-6
E(i,3)=fmax+E(i,5);
end
end
[Q,IX]=sort(E,1);
%Q=vpa(Q,4);%%
D(g,1)=Q(1,3);
D(g,2)=mean(Q(:,3));
D(g,3)=Q(size,3);
D(g,4)=E(IX(1,3),4);
if Q(1,3)<fv && D(g,4)==1
fv=Q(1,3);
xm=E(IX(1,3),(1:2));
xm=[xm,E(IX(1,3),4)];
end
end
%畫圖
k=1;
for i=1:NP
if D(i,4)==1
N(k,:)=D(i,:);
k=k+1;
end
end
plot(N(:,1),'r*');
hold on
plot(N(:,2),'b+');
hold on
plot(N(:,3),'ms');
legend('最優值','平均值','最差值');
hold off
end