运用matlab编程实现蚁群算法的低压电力线设计1.zip

clear all; clc; %%第一步：变量初始化 m=40;Alpha=1;Beta=1;Rho=0.1;NC_max=75;Q=100;q0=0.2;s=2;d=23; C =[ 37 129 94 424 361 5 91 14 380 373 200 288 163 217 21 211 155 389 430 439 104 40 298 0 489 275 488 461 25 282 316 426 451 258 319 222 127 217 289 313 31 380 135 212 424 266 297 156 421 156 33 432 330 389 351 124 352 253 325 142 221 239 26 423 282 326 43 30 386 150 186 156 374 259 377 270 84 437 305 212 371 83 499 277 487 368 52 266 432 243 277 410 77 327 180 308 222 151 371 354 221 490 274 11 324 454 70 341 66 305 150 314 309 372 197 74 112 350 471 29 ]; %scatter(C(:,1),C(:,2)); %hold on; n=size(C,1);%n表示问题的规模（城市个数）结果：n=40 得到第一列维数。即40个城市。 D=ones(n,n);%D表示完全图的赋权邻接矩阵 31*31矩阵 %%%%%%%%求距离 for i=1:n for j=1:n if i<j D(i,j)=sqrt((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2); end D(j,i)=D(i,j); end end %D %%%%%%%分布时延 delay=zeros(n,n); for i=1:n for j=1:n if D(i,j)==1 delay(i,j)=0; elseif D(i,j)>=1&D(i,j)<100 delay(i,j)=0.002; elseif D(i,j)>=100 & D(i,j)<200 delay(i,j)=0.004; elseif D(i,j)>=200 & D(i,j)<300 delay(i,j)=0.006; elseif D(i,j)>=300 & D(i,j)<400 delay(i,j)=0.008; elseif D(i,j)>=400 & D(i,j)<500 delay(i,j)=0.01; elseif D(i,j)>=500 & D(i,j)<600 delay(i,j)=0.012; end end end %delay Eta=1./D;%Eta为启发因子，这里设为距离的倒数 Tau=ones(n,n);%Tau为信息素矩阵 Tabu=zeros(m,n);%存储并记录路径的生成 NC=1;%迭代计数器 R_best=zeros(NC_max,n);%各代最佳路线 L_best=inf.*ones(NC_max,1);%各次最佳路线的时延 inf表示正无穷大 while NC<=NC_max%停止条件之一：达到最大迭代次数 %%%第二步：放在初始位置 Tabu(:,1)=s; %%%%第三步：建立路由表 R=zeros(n,n); for i=1:n for j=1:n if D(i,j)<350 R(i,j)=1; else R(i,j)=0; end end end R(i,j)=R(i,j); %%%%第四步：蚂蚁按概率计算下一跳 for i=1:m%蚂蚁依次循环 for j=2:n%城市 visited=Tabu(i,1:(j-1));%已访问的城市，是个1行j-1列的向量，纪录了第i只蚂蚁走过的城市号 避免环路 J=zeros(1,(n-j+1));%待访问的城市，是个1行n-j+1列的向量 P=J;%待访问城市的选择概率分布 Jc=1;%Jc在这里是什么意思呢？没什么意思，只是为Jc赋初值 for k=1:n %寻找第i只蚂蚁未走过的城市编号 if(length(find(visited==k))==0)&(R(visited(end),k)==1) J(Jc)=k;%将未走过的城市编号依次赋给J中第Jc个元素 Jc=Jc+1;%通过循环，进行下一个未走过城市的查找 end % end %这个循环又是什么作用？ end if J(1)==0 break; end %%%%%下面计算待选城市的概率分布 for k=1:length(J) if J(k)~=0 P(k)=(Tau(visited(end),J(k))^Alpha)*(Eta(visited(end),J(k))^Beta);%Tau是个n行n列的矩阵，存储着任一城市和其它城市之间的信息素浓度，visited(end)表示现在蚂蚁的位置，J(k)中存储着可选择的各城市编号；Eta同理 end end P=P/(sum(P));%计算各个城市的选择概率 %%%%%%按概率原则选取下一个城市 q=rand; if q<=q0 Select=find(P==max(P)); to_visit=J(Select(1)); Tabu(i,j)=to_visit; elseif q>q0 Pcum=cumsum(P);%cumsum(A)--计算各维元素的和；例:P = 0.0169 0.0059 0.9772 ,则 Pcum =0.0169 0.0228 1.0000 Select=find(Pcum>=rand);%此句是什么意思呢? rand本身是产生一个从0到1随机分布的小数,为什么将Pcum>=rand作为选择下一个城市的判决条件?--算法的规定! to_visit=J(Select(1));%将向量J中第Select(1)个值作为下一个要访问的城市编号 Tabu(i,j)=to_visit;%将本次访问的城市编号纪录到禁忌矩阵中 end if Tabu(i,j)==d break; end end end Tabu(i,j)=Tabu(i,j); %%第五步：记录本次迭代最佳路线 L=zeros(m,1);%l长为m的列向量,纪录了m只蚂蚁各自周游的总长度 for i=1:m %第m只蚂蚁 V=Tabu(i,:);%V为一行向量,纪录了第i只蚂蚁依次周游的各个城市 for j=1:(n-1) if V(j+1)~=0 L(i)=L(i)+delay(V(j),V(j+1)); else break; end end end L_best(NC)=min(L); pos=find(L==L_best(NC));%找出是第几只蚂蚁周游的距离最短，其中执行find命令后返回的是元素的序列号，而不是符合条件的元素本身 R_best(NC,:)=Tabu(pos(1),:);%将第NC次产生的最优路径纪录到R_best矩阵中 因为走出最短路径的有好多只，所以用pos(1),将第一只的路径放入r_best。 L_ave(NC)=mean(L); NC=NC+1; %%%%第六步：更新信息素 Delta_Tau=zeros(n,n); for i=1:m for j=1:(n-1) if Tabu(i,j+1)~=0 Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i); else break; end end end Tau=(1-Rho).*Tau+Delta_Tau; %%%%%第七步：禁忌表清零 Tabu=zeros(m,n); end L_best L_ave pos=find(L_best==min(L_best)); shortest_Route=R_best(pos(1),:) shortest_length=L_best(pos(1)) % %%%%第八步：目的节点处理 % %%%%%1.排序 % R_best; % B=zeros(1,n); % for j=1:(NC_max-1) % for i=1:(NC_max-j) % if L_best(i)>L_best(i+1) % B=R_best(i,:); % R_best(i,:)=R_best((i+1),:); % R_best((i+1),:)=B; % A=L_best(i); % L_best(i)=L_best(i+1); % L_best(i+1)=A; % end % end % end % %L_best % r_best=R_best % % M=zeros(1,n) % % for i=1:m %第m只蚂蚁 % % V=r_best(i,:);%V为一行向量,纪录了第i只蚂蚁依次周游的各个城市 % % for j=1:(n-1) % % if V(j+1)~=0 % % M(i)=M(i)+delay(V(j),V(j+1)); % % else % % break; % % end % % end % % end % % %%%%%2.取出主路径上的邻节点，放入F(1,n)中 % Nbr=zeros(length(shortest_Route),n); % A=shortest_Route(shortest_Route~=0); % for i=1:length(A) % k=1; % for j=1:n % if R(A(i),j)==1 % Nbr(i,k)=j; % k=k+1; % end % end % end % Nbr; % F=zeros(1,n); % F=unique(Nbr); % % %%%%%3.判断是否是简单相关 % % for i=1:NC_max % for j=1:n % flag=0; % for k=1:length(F) % if R_best(i,j)==F(k) % flag=1; % break; % else % flag=0; % end % end % if flag==0 % R_best(i,:)=0;%如果不是简单相关路径则把这行赋0. % break; % end % end % end %%%%第九步 画路线图 % R_best % R_best(~any(R_best,2),:)=[];%删除全0行 % % unique(R_best,'rows');%删除重复 % subplot(1,2,1) % for i=1:3 % R_best(i,:) % DrawRoute(C,R_best(i,:)); % hold on % end subplot(1,2,1) plot(L_best) hold on plot(L_ave)