蚁群_VRP,蚁群算法,matlab

%%Q：信息更新参数 %%W:车辆载重量 clear clc close all tic %%初始化变量和系数 m=20; %m：蚂蚁数目 Alpha=1; %Alpha：重要度系数 Beta=1; %Beta：能见度系数 gama=5; Rho=0.15; %Beta：能见度系数 NC_max=100;%NC_max：最大迭代次数 Q=15; W=9; qq=0.05; C=[0 0; 0 -1; 0 3; -1 -2; -3 -3; 3 -1; -4 0; -4 -1; 1 -2; 1 -1; 1 3; 3 4; -3 0; 2 0; 1 -3; 2 -1; 2 1; 1 -4; -3 2; -1 -1]; %C：DC和工厂的坐标 t=[0 1.5 1.8 2 0.8 1.5 1.0 2.5 3.0 1.7 0.6 0.2 2.4 1.9 2.0 0.7 0.5 2.2 3.1 0.1]; %构造仓库和工厂等之间的距离矩阵 n=size(C,1);%n表示问题的规模（城市个数） D=Distanse(C); %构造节省量矩阵 U=zeros(n,n);%U表示工厂之间的连接和其与仓库之间连接能节省的距离 for i=1:n U(i,i)=eps; for j= 1:n U(i,j)=D(i,1)+D(j,1)-D(i,j); U(j,i)= U(i,j); end end U=U+1; load_w=0; Eta=1./D;%Eta为启发因子，这里设为距离的倒数 Tau=ones(n,n);%Tau为信息素矩阵 Tabu=zeros(m,n+4);%存储并记录路径的生成 NC=1;%迭代计数器 G_best_route=[NC_max,n+4];%各代最佳路线 G_best_length=inf.*ones(NC_max,1);%各代最佳路线的长度 length_ave=zeros(NC_max,1);%各代路线的平均长度 Obj=zeros(NC_max,1); %%第二步，把蚂蚁放到DC内 while NC<=NC_max%停止条件之一：达到最大迭代次数 Tabu(:,1)=ones(m,1); %%第三步，m只蚂蚁按照要求的方法选择工厂，并完成周游 for i=1:m for j=2:n+4 visited=Tabu(i, :); visited=visited(visited>0); to_visit=setdiff(1:n,visited); x=zeros(length(to_visit),1); if isempty(to_visit) % 按照规则选下一个工厂或者是回到仓库 Tabu(i,length(visited)+1)=1; load_w=0; continue; end for k=1:length(to_visit) x(k)=(Tau(visited(end),to_visit(k))^Alpha)*(Eta(visited(end),to_visit(k))^Beta);%*(U(visited(end),to_visit(k))^gama); end if rand<qq Select=find(max(x)); else x=x/(sum(x)); %按概率原则选取下一个城市 xcum=cumsum(x); Select=find(xcum>=rand,1,'first'); end load_w=load_w+t(to_visit(Select)); if load_w>W load_w=0; Tabu(i,length(visited)+1)=1; else Tabu(i,length(visited)+1)=to_visit(Select); end end end %% 第四步记录本代各种参数 L=zeros(m,1); for i=1:m MM=Tabu(i,:); R=MM(MM>0); for j=1:(length(R)-1) L(i)=L(i)+D(R(j),R(j+1)); end end [G_best_length(NC),pos]=min(L); G_best_route(NC,1:length(Tabu(pos(1),:)))=Tabu(pos(1),:); %%应用2-opt方法对最优解进行更新 select=find(G_best_route(NC,:)==1); FF=[];%记录新的路线 LM=0;%计算新的值 for a=1:(length(select)-1) y_G_best_route=G_best_route(NC,select(a):select(a+1));%选取第一条路线 al=length(y_G_best_route); T=zeros((length(select)-1),1);%记录各条路线的长度 %计算各条路线的长度 for d=1:(al-1) T(a)=T(a)+D(y_G_best_route(d),y_G_best_route(d+1)); end for b=2:(al-1) for c=(b+1):(al-1) DD=y_G_best_route; %交换位置 temp1=DD(b); temp2=DD(c); DD(b)=temp2; DD(c)=temp1; %记录距离 TT=zeros(1); for d=1:(al-1)%计算交换位置后的距离 TT=TT+D(DD(d),DD(d+1)); end if TT<T(a) T(a)=TT; y_G_best_route=DD; end end end if a>=2 y_G_best_route=y_G_best_route(2:al); end %更新路线和长度 FF=[FF,y_G_best_route]; LM=LM+T(a); end %更新最优路线和长度 G_best_length(NC)=LM; G_best_route(NC,1:length(FF))=FF; FF=[]; LM=0; %%2-opt优化完毕 length_ave(NC)=mean(L); %% 第五步更新信息素 Delta_Tau=zeros(n,n); for i=1:m MM=Tabu(i,:); R=MM(MM>0); for j=1:(length(R)-1) Delta_Tau(R(j),R(j+1))=Delta_Tau(R(j),R(j+1))+Q/L(i); end end Tau=(1-Rho).*Tau+Delta_Tau; %% 第六步：禁忌表清零 Tabu=zeros(m,n); load_w=0; if NC==1 || G_best_length(NC)<Obj(NC-1) Obj(NC)=G_best_length(NC); else Obj(NC)=Obj(NC-1); end NC=NC+1; end %%第七步：输出结果 [best_length,Pos]=min(G_best_length); best_route=G_best_route(Pos(1),:); best_route=best_route(best_route>0); disp(best_route); disp('最小行程为:') disp(best_length); toc %% 最优解的路径图 DrawPath(best_route(1:end-1),C) %% 优化过程迭代图 figure plot(1:NC_max,Obj) hold on plot(length_ave) xlabel('迭代次数') ylabel('距离') title('优化过程')