aco.txt

一共四个蚁群算法 程序源码

%    the procedure of ant colony algorithm for VRP
%
%    %    %    %    %    %    %    %    %    %    %

%initialize the parameters of ant colony algorithms
load data.txt;
d=data(:,2:3);
g=data(:,4);

m=31; % 蚂蚁数
alpha=1;
belta=4;% 决定tao和miu重要性的参数
lmda=0;
rou=0.9;　%衰减系数
q0=0.95;
% 概率
tao0=1/(31*841.04);%初始信息素
Q=1;%　蚂蚁循环一周所释放的信息素
defined_phrm=15.0;   % initial pheromone level value 
QV=100;  % 车辆容量
vehicle_best=round(sum(g)/QV)+1;　%所完成任务所需的最少车数
V=40;
% 计算两点的距离 
for i=1:32;
    for j=1:32;
       dist(i,j)=sqrt((d(i,1)-d(j,1))^2+(d(i,2)-d(j,2))^2);
    end;
end;
%给tao miu赋初值
for i=1:32;
       for j=1:32;
           if i~=j;
               %s(i,j)=dist(i,1)+dist(1,j)-dist(i,j);
               tao(i,j)=defined_phrm;
               miu(i,j)=1/dist(i,j); 
           end;                  
       end;
end;

for k=1:32;
     for k=1:32;
         deltao(i,j)=0;
     end;
end;            
best_cost=10000;       
for n_gen=1:50;
   print_head(n_gen); 
  for i=1:m;
     %best_solution=[];
     print_head2(i);
     sumload=0;
     cur_pos(i)=1;
     rn=randperm(32);
     n=1;
     nn=1;
     part_sol(nn)=1;
     %cost(n_gen,i)=0.0;
     n_sol=0;   % 由蚂蚁产生的路径数量
     M_vehicle=500;
     t=0;  %最佳路径数组的元素数为0
          
     while sumload<=QV;
               
        for k=1:length(rn);
            if sumload+g(rn(k))<=QV;
                gama(cur_pos(i),rn(k))=(sumload+g(rn(k)))/QV;
                A(n)=rn(k);
                n=n+1;
            end;
        end;
       fid=fopen('out_customer.txt','a+');
        fprintf(fid,'%s  %i\t','the current position is:',cur_pos(i));  
        fprintf(fid,'\n%s','the possible customer set is:')
        fprintf(fid,'\t%i\n',A);
        fprintf(fid,'------------------------------\n');
        fclose(fid);
         p=compute_prob(A,cur_pos(i),tao,miu,alpha,belta,gama,lmda,i);
        maxp=1e-8;
        na=length(A);
        for j=1:na;
               if p(j)>maxp
                   maxp=p(j);
                   index_max=j;
               end;
        end;
           
        old_pos=cur_pos(i);
        if rand(1)<q0
            cur_pos(i)=A(index_max);
        else 
            krnd=randperm(na);
            cur_pos(i)=A(krnd(1));       
            bbb=[old_pos cur_pos(i)];
            ccc=[1 1];
            if bbb==ccc;
                cur_pos(i)=A(krnd(2)); 
            end;
        end;
        
        tao(old_pos,cur_pos(i))=taolocalupdate(tao(old_pos,cur_pos(i)),rou,tao0);%对所经弧进行局部更新
        
        sumload=sumload+g(cur_pos(i));

        nn=nn+1;
        part_sol(nn)=cur_pos(i);
        temp_load=sumload;
                           
        if cur_pos(i)~=1;
            rn=setdiff(rn,cur_pos(i));
            n=1;
            A=[];
        end;
        
        if cur_pos(i)==1;  % 如果当前点为车场,将当前路径中的已访问用户去掉后,开始产生新路径
           if setdiff(part_sol,1)~=[];
                n_sol=n_sol+1;  % 表示产生的路径数,n_sol=1,2,3,..5,6...,超过5条对其费用加上车辆的派遣费用
                fid=fopen('out_solution.txt','a+');
                fprintf(fid,'%s%i%s','NO.',n_sol,'条路径是:');
                fprintf(fid,'%i  ',part_sol);
                fprintf(fid,'\n');
                fprintf(fid,'%s','当前的用户需求量是:');
                fprintf(fid,'%i\n',temp_load);
                fprintf(fid,'------------------------------\n');
                fclose(fid);       
                
                % 对所得路径进行路径内3-opt优化
                final_sol=exchange(part_sol);
                             
                for nt=1:length(final_sol); % 将所有产生的路径传给一个数组
                    temp(t+nt)=final_sol(nt);
                end;
                t=t+length(final_sol)-1;
                
                sumload=0;
                final_sol=setdiff(final_sol,1);
                rn=setdiff(rn,final_sol);
                part_sol=[];
                final_sol=[];
                nn=1;
                part_sol(nn)=cur_pos(i);
                A=[];
                n=1;
                
            end;   
        end;
                
        if setdiff(rn,1)==[];% 产生最后一条终点不为1的路径
            n_sol=n_sol+1;
            nl=length(part_sol);
            part_sol(nl+1)=1;%将路径的最后1位补1
            
            % 对所得路径进行路径内3-opt优化
            final_sol=exchange(part_sol);           
           
            for nt=1:length(final_sol); % 将所有产生的路径传给一个数组
                temp(t+nt)=final_sol(nt);                
            end;
            
            cost(n_gen,i)=cost_sol(temp,dist)+M_vehicle*(n_sol-vehicle_best);   %计算由蚂蚁i产生的路径总长度
            
            for ki=1:length(temp)-1;
                deltao(temp(ki),temp(ki+1))=deltao(temp(ki),temp(ki+1))+Q/cost(n_gen,i);
            end; 
           
            if cost(n_gen,i)<best_cost;
                best_cost=cost(n_gen,i);
                old_cost=best_cost;
                best_gen=n_gen;  % 产生最小费用的代数
                best_ant=i; %产生最小费用的蚂蚁
                best_solution=temp;
            end;
                                  
            if i==m;　　%如果所有蚂蚁均完成一次循环，，则用最佳费用所对应的路径对弧进行整体更新
                for ii=1:32;
                    for jj=1:32;
                        tao(ii,jj)=(1-rou)*tao(ii,jj);
                    end;
                end;
                
                for kk=1:length(best_solution)-1;
                    tao(best_solution(kk),best_solution(kk+1))=tao(best_solution(kk),best_solution(kk+1))+deltao(best_solution(kk),best_solution(kk+1));
                end; 
            end;     
                      
            fid=fopen('out_solution.txt','a+');
            fprintf(fid,'%s%i%s','NO.',n_sol,'路径是:');
            fprintf(fid,'%i ',part_sol);
            fprintf(fid,'\n');
            fprintf(fid,'%s %i\n','当前的用户需求量是:',temp_load);
            fprintf(fid,'%s %f\n','总费用是:',cost(n_gen,i));
            fprintf(fid,'------------------------------\n');
            fprintf(fid,'%s\n','最终路径是:');
            fprintf(fid,'%i-',temp);
            fprintf(fid,'\n');
            fclose(fid); 
            temp=[];
            break;
        end;
    end;
    
  end;
end;

 

 




第二个：


function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta,Rho,Q)

%%-------------------------------------------------------------------------

%% 主要符号说明

%% C n个城市的坐标，n×2的矩阵

%% NC_max 最大迭代次数

%% m 蚂蚁个数

%% Alpha 表征信息素重要程度的参数

%% Beta 表征启发式因子重要程度的参数

%% Rho 信息素蒸发系数

%% Q 信息素增加强度系数

%% R_best 各代最佳路线

%% L_best 各代最佳路线的长度

%%=========================================================================


%%第一步：变量初始化

n=size(C,1);%n表示问题的规模（城市个数）

D=zeros(n,n);%D表示完全图的赋权邻接矩阵

for i=1:n

for j=1:n

if i~=j

D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;

else

D(i,j)=eps;      %i=j时不计算，应该为0，但后面的启发因子要取倒数，用eps（浮点相对精度）表示

end

D(j,i)=D(i,j);   %对称矩阵

end

end

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);   %各代最佳路线的长度

L_ave=zeros(NC_max,1);        %各代路线的平均长度


while NC<=NC_max        %停止条件之一：达到最大迭代次数，停止

%%第二步：将m只蚂蚁放到n个城市上

Randpos=[];   %随即存取

for i=1:(ceil(m/n))

Randpos=[Randpos,randperm(n)];

end

Tabu(:,1)=(Randpos(1,1:m))';    %此句不太理解？


%%第三步：m只蚂蚁按概率函数选择下一座城市，完成各自的周游

for j=2:n     %所在城市不计算

for i=1:m     

visited=Tabu(i,1:(j-1)); %记录已访问的城市，避免重复访问

J=zeros(1,(n-j+1));       %待访问的城市

P=J;                      %待访问城市的选择概率分布

Jc=1;

for k=1:n

if length(find(visited==k))==0   %开始时置0

J(Jc)=k;

Jc=Jc+1;                         %访问的城市个数自加1

end

end

%下面计算待选城市的概率分布

for k=1:length(J)

P(k)=(Tau(visited(end),J(k))^Alpha)*(Eta(visited(end),J(k))^Beta);

end

P=P/(sum(P));

%按概率原则选取下一个城市

Pcum=cumsum(P);     %cumsum，元素累加即求和

Select=find(Pcum>=rand); %若计算的概率大于原来的就选择这条路线

to_visit=J(Select(1));

Tabu(i,j)=to_visit;

end

end

if NC>=2

Tabu(1,:)=R_best(NC-1,:);

end


%%第四步：记录本次迭代最佳路线

L=zeros(m,1);     %开始距离为0，m*1的列向量

for i=1:m

R=Tabu(i,:);

for j=1:(n-1)

L(i)=L(i)+D(R(j),R(j+1));    %原距离加上第j个城市到第j+1个城市的距离

end

L(i)=L(i)+D(R(1),R(n));      %一轮下来后走过的距离

end

L_best(NC)=min(L);           %最佳距离取最小

pos=find(L==L_best(NC));

R_best(NC,:)=Tabu(pos(1),:); %此轮迭代后的最佳路线

L_ave(NC)=mean(L);           %此轮迭代后的平均距离

NC=NC+1                      %迭代继续



%%第五步：更新信息素

Delta_Tau=zeros(n,n);        %开始时信息素为n*n的0矩阵

for i=1:m

for j=1:(n-1)

Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i);           

%此次循环在路径（i，j）上的信息素增量

end

Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);

%此次循环在整个路径上的信息素增量

end

Tau=(1-Rho).*Tau+Delta_Tau; %考虑信息素挥发，更新后的信息素

%%第六步：禁忌表清零

Tabu=zeros(m,n);             %%直到最大迭代次数

end

%%第七步：输出结果

Pos=find(L_best==min(L_best)); %找到最佳路径（非0为真）

Shortest_Route=R_best(Pos(1),:) %最大迭代次数后最佳路径

Shortest_Length=L_best(Pos(1)) %最大迭代次数后最短距离

subplot(1,2,1)                  %绘制第一个子图形

DrawRoute(C,Shortest_Route)     %画路线图的子函数

subplot(1,2,2)                  %绘制第二个子图形

plot(L_best)

hold on                         %保持图形

plot(L_ave,'r')

title('平均距离和最短距离')     %标题

function DrawRoute(C,R)

%%=========================================================================

%% DrawRoute.m

%% 画路线图的子函数

%%-------------------------------------------------------------------------

%% C Coordinate 节点坐标，由一个N×2的矩阵存储

%% R Route 路线

%%=========================================================================


N=length(R);

scatter(C(:,1),C(:,2));

hold on

plot([C(R(1),1),C(R(N),1)],[C(R(1),2),C(R(N),2)],'g')

hold on

for ii=2:N

plot([C(R(ii-1),1),C(R(ii),1)],[C(R(ii-1),2),C(R(ii),2)],'g')

hold on

end

title('旅行商问题优化结果 ')




这是一个例子：



1.function Q = FINDabc(A) %给定B1，B2两点坐标计算平面参数a，b，c 
v = 1e-3*[1 1 1 1];r=1.3;l=0;
K=6;
a = RAND;
a = RESHAP_MAX(A,a,K);
while(sum(STAND(A,a(K,:))>v)~=0)
l=0;
[Q,w]=REF(a,r,K);
while((~D_ZONE1(Q))|(COMP(A,Q)>=COMP(A,a(K,:))))
[Q,w]=REF(a,w,K); 
if (w<=eps^12)

a(K,:)=[];
K=K-1; 
l = 1;
end
if(l==1)
break;
end
end
if(~l)
a(K,:)=Q; 
end 
if(K<=1)
a=FINDabc(A);
break;
end 
a = RESHAP_MAX(A,a,K);
end
if(K<=1)
Q = a(1,:);
else Q=sum(a)./K;
end

11
2.function P=TEXT(i) %对i组[a,b,c]检验，返回[a,b,c,a',b',c']
j = 1;
while(j<=i)
A(1,1:2)= CREAT_RAND12;
A(1,3)=CREAT_RAND3;
if(D_ZONE(A(1,:))==1)
b(j,:)=A(1,:);j=j+1;
end
end
for j=1:i
c(j,:)=Fx(b(j,:));
d(j,:)=FINDabc(c(j,:));
e(j,:)=(d(j,:)-b(j,:))./b(j,:);
end
P=[b,d,e];
plot(e(:,1),'b--o');hold on;
plot(e(:,2),'g-.x');
plot(e(:,3),'r:+');


3.function p = recognize(i) %分辨率的影响
j = 1;
while(j<=i)
a(1,1:2) = CREAT_RAND12 ();
a(1,3) = CREAT_RAND3();
if(D_ZONE(a(1,:)))
b(j,:) = a(1,:);
j = j+1;
end
end
for j = 1:i
c(j,:) = Fx(b(j,:));
e(j,:)=1e-4*round(1e4*c(j,:));
end
for j = 1:i 
f(j,1:3) = FINDabc(e(j,:));
end 
hl = (b-f)./b;
for j = 1:3
h(:,j)=sum(abs(hl(:,j)))/i*ones(i,1);
end
plot(h(:,1),'--r');hold on;plot(h(:,2),'-y');plot(h(:,3),':k');
p = h;
12
4.function Q_C12 =CREAT_RAND12 () 
if(y_C12<=0.5)
Q_C12(1) =0.0017 *rand(1);
else
Q_C12(1) = -0.0017*rand(1);
end
y_C12 = rand(1);
if(y_C12<=0.5)
Q_C12(2) = 0.0017*rand(1);
else
Q_C12(2) = -0.0017*rand(1);
End


5.function Q_C = CREAT_RAND3() 
y_C = rand(1);
if(y_C<=0.5)
Q_C = 3*1e-6*rand(1);
else
Q_C = -3*1e-6*rand(1);
end


6.function y = D_ZONE(Q_D) 
p = pi/180/10;