clc
clear
plane=[2 16 610 1; %plane performance matrix 飞机性能
2 21 815 1;
2 15 540 1
2 9 750 1];
xp=-1000; yp=-1000; %air port location 机场位置
x0=0; y0=0; %search boundery 搜索边界
x1=1000; y1=500;
dx=20; dy=-20; %ocean current 洋流
p=2; %the density of grid(km per grid) 网格密度
Q=100; %quantity of population
gen=240; %genetic algebra 遗传代数
mode=1; %概率平均分布模式
width=50; %the width of scan area
M=0.3; %mutation rate;
%initial 初始化
plane_num=sum(plane(:,1));
[r,c]=size(plane);
grad(1)=plane(1,1);
for i=2:r
grad(i)=grad(i-1)+plane(i,1);
end
population=zeros(Q,2*plane_num);
population(:,1:2:end)=unidrnd((x1-x0)/p+1,[Q,plane_num])-1+x0/p;
population(:,2:2:end)=unidrnd(((y1-y0)/p+1),[Q,plane_num])-1+y0/p;
if mode==1
P=zeros((x1-x0)/p,(y1-y0)/p+1)+(1/(((x1-x0)/p*(y1-y0)/p+1)));
end
[P_r,P_c]=size(P);
for j=1:Q
p_temp(:,:,j)=P;
fitness(j)=0;
for i=1:plane_num
flag=sum(i>grad)+1;%标记具体型号飞机
x=population(j,2*i-1);
y=population(j,2*i);
d_left(i)=plane(flag,2)*plane(flag,3)-p*2*sqrt((x-xp/p)^2+(y-yp/p)^2);
area=d_left(i)*plane(flag,4);
len=area/width;
width_net=width/p;
len_net=floor(len/p);
%constrct the plane scan area
for m=1:len_net
for n=1:width_net
if x+m>P_r||y+n>P_c
continue;
end
fitness(j)=fitness(j)+p_temp(x+m,y+n,j);
p_temp(x+m,y+n,j)=-0.001;
end
end
end
end
for t=1:gen
hold on
[Q,W]=size(population);
% st=3;
% %单点交叉
% for i=st:2:Q
% if rand<0.5
% loc=ceil(rand*W);
% temp=population(i,1:loc);
% population(i,1:loc)=population(i+1,1:loc);
% population(i+1,1:loc)=temp;
% end
% end
t
%choose
popu_table=[population fitness'];
popu_table=sortrows(popu_table,-(W+1));
population=[];
population(Q,W)=0;
population(1,:)=popu_table(1,1:end-1);
fitness=sort(fitness);
fitness(end);
deno=sum(1:Q-1);
p_(1)=1/deno;
for i=2:Q-1
p_(i)=p_(i-1)+i/deno;
end
for i=2:Q
fate=rand;
[r,c]=find(p_>fate);
loc=length(c)+1;
population(i,:)=popu_table(loc,1:end-1);
end
clear popu_table loc deno p_ r c
%fitness
for j=1:Q
% if j==1
% % fill([0 x1/p x1/p 0 ],[0 0 y1/p y1/p],'w')
% drawnow
% end
p_temp(:,:,j)=P;
fitness(j)=0;
for i=1:plane_num
flag=sum(i>grad)+1;%标记具体型号飞机
x=population(j,2*i-1);
y=population(j,2*i);
d_left(i)=plane(flag,2)*plane(flag,3)-p*2*sqrt((x-xp/p)^2+(y-yp/p)^2);
area=d_left(i)*plane(flag,4);
len=area/width;
width_net=width/p;
len_net=floor(len/p);
%constrct the plane scan area
% for m=1:len_net
% for n=1:width_net
% if x+m>P_r||y+n>P_c
% continue;
% end
% fitness(j)=fitness(j)+p_temp(x+m,y+n,j);
% p_temp(x+m,y+n,j)=-0.0001;
% end
% end
if x+len_net>=P_r && y+width_net>=P_c %%%%右上出界
fitness(j)=fitness(j)+sum(sum(p_temp(x:end,y:end,j)));
p_temp(x:end,y:end,j)=-0.0001;
if flag==1 &&j==1
fill([x x1/p x1/p x ],[y y y1/p y1/p],'r','EraseMode','xor')
elseif flag==2 &&j==1
fill([x x1/p x1/p x ],[y y y1/p y1/p],'g','EraseMode','xor')
elseif flag==3 &&j==1
fill([x x1/p x1/p x ],[y y y1/p y1/p],'b','EraseMode','xor')
elseif flag==4 &&j==1
fill([x x1/p x1/p x ],[y y y1/p y1/p],'y','EraseMode','xor')
end
axis([0 x1/p 0 y1/p]);
drawnow
break;
elseif x+len_net>=P_r %&& y+width_net>P_c 右出界
fitness(j)=fitness(j)+sum(sum(p_temp(x+1:end,y+1:y+width_net,j)));
p_temp(x:end,y+width_net-1,j)=-0.0001;
if flag==1 &&j==1
fill([x x1/p x1/p x ],[y y y+width_net y+width_net],'r')
elseif flag==2 &&j==1
fill([x x1/p x1/p x ],[y y y+width_net y+width_net],'g')
elseif flag==3 &&j==1
fill([x x1/p x1/p x ],[y y y+width_net y+width_net],'b')
elseif flag==4 &&j==1
fill([x x1/p x1/p x ],[y y y+width_net y+width_net],'y')
end
axis([0 x1/p 0 y1/p]);
drawnow
break;
elseif y+width_net>=P_c %上出界
fitness(j)=fitness(j)+sum(sum(p_temp(x+1:x+len_net,y+1:end,j)));
p_temp(x:x+len_net-1,y:end,j)=-0.0001;
if flag==1 &&j==1
fill([x x+len_net x+len_net x ],[y y y1/p y1/p],'r')
elseif flag==2 &&j==1
fill([x x+len_net x+len_net x ],[y y y1/p y1/p],'g')
elseif flag==3 &&j==1
fill([x x+len_net x+len_net x ],[y y y1/p y1/p],'b')
elseif flag==4 &&j==1
fill([x x+len_net x+len_net x ],[y y y1/p y1/p],'y')
end
axis([0 x1/p 0 y1/p]);
drawnow
break;
else
if x+len_net-1>P_r ||y+width_net-1>P_c
break;
end
check1=x+len_net-1;
check2=y+width_net-1;
tmp1=(sum(p_temp(x+1:x+len_net,y+1:y+width_net,j)));
tmp2=sum(tmp1);
fitness(j)=fitness(j)+tmp2;
p_temp(x+1:x+len_net,y+1:y+width_net,j)=-0.0001;
if flag==1 &&j==1
fill([x x+len_net x+len_net x ],[y y y+width_net y+width_net],'r')
elseif flag==2 &&j==1
fill([x x+len_net x+len_net x ],[y y y+width_net y+width_net],'g')
elseif flag==3 &&j==1
fill([x x+len_net x+len_net x ],[y y y+width_net y+width_net],'b')
elseif flag==4 &&j==1
fill([x x+len_net x+len_net x ],[y y y+width_net y+width_net],'y')
end
axis([0 x1/p 0 y1/p]);
drawnow
end
end
end
hold off
% [x__,y__]=find(p_temp(:,:,1)==-0.0001);
% h=plot(x__,y__,'r.','EraseMode','xor');
%
% drawnow
%mutation
% population(2:end,:)=population(2:end,:)+population(2:end,:).*(rand-0.5)*M;
if t<gen/4*1
population(2:end,:)=population(2:end,:)+20*(rand(Q-1,plane_num*2)-0.5);
elseif t<gen/4*2
population(2:end,:)=population(2:end,:)+10*(rand(Q-1,plane_num*2)-0.5);
elseif t<gen/4*3
population(2:end,:)=population(2:end,:)+5*(rand(Q-1,plane_num*2)-0.5);
else
population(2:end,:)=population(2:end,:)+1*(rand(Q-1,plane_num*2)-0.5);
end
population=floor(population);
population(population<=0)=1;
if t~=gen
clf reset
end
xlabel('Length');
ylabel('Width');
title('Danamic Result Of Solution ')
legend('P-3C','P-8A','C-130','IL-76');
elite(t)=max(fitness);
max(fitness)
elite(t)=max(elite);
end
figure
plot(1:gen,elite,'rs','LineWidth',1,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',3)
xlabel('Genneration');
ylabel('Fitness');
title('The curve of optimization')
