anttsp.m

一个用matlab编写基本的蚁群算法解决TSP问题程序

%file:  anttsp.m
%保留每次迭代的最优解
%以max(t^a*d^(-b))为依据找最优路径,与保留的最优路径比较
%n表示城市数目
%r表示轨迹持久性
%a表示轨迹相对重要性
%b表示能见度相对重要性
%m表示蚂蚁数目
clear
x=[41 37 54 25 7 2 68 71 54 83 64 18 22 83 91 ...
25 24 58 71 74 87 18 13 82 62 58 45 41 44 4];
y=[94 84 67 62 64 99 58 44 62 69 60 54 60 46 38 ...
38 42 69 71 78 76 40 40 7 32 35 21 26 35 50];
n=30;%n表示城市数目
c=100;%c表示初始信息浓度
q=1000000;
NC=100;
r=0.9;%r表示轨迹持久性
a=1.5;%a表示轨迹相对重要性
b=2;%b表示能见度相对重要性
m=50;%m表示蚂蚁数目
for  i=1:n 
    for j=1:n
dij(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);%dij保存任意两个城市之间的距离
    end
end
for i=1:n
dij(i,i)=0.01;%自己和自己的距离定义为0.01
end
min10=100000;%
t=ones(n)*c;%t表示信息浓度
for nc=1:NC;%NC为迭代次数
tabu=ones(m,n);%tabu表示蚂蚁访问过的城市集合，即禁忌表
tabu(:,1)=0;%表示第一个城市已经过，不允许再去该城市，所以第一列全为0
path=ones(m,n);%表示每只蚂蚁所走的路径
for k=1:m%k表示第k只蚂蚁
    for step=1:n-1%表示第k只蚂蚁第step次移动
        ta=t.^a;
        tb=dij.^(-b);
        td=ta.*tb;
        pd=tabu(k,:).*td(path(k,step),:);%此处蚂蚁下一步允许的城市，是由tabu(k,:)决定，path(k,step)表示蚂蚁所在的城市
        pk=pd/sum(pd);%pk表示转移概率，
        rk=rand;%产生一个0到1之间的随机数
        spk=0;%？？
        j=1;
        while j<=n;%此处的循环只要产生的随机数小于转移概率之和（前面的）就跳出
            if rk<spk+pk(j)
                break;
            else
                spk=spk+pk(j);
                j=j+1;
            end
        end
        tabu(k,j)=0;%改变禁忌表，使第k行第j个元素变0？？
        path(k,step+1)=j;%修改path，使第k行第step+1个元素变j？？
    end
end                        %%到此处为止，m只蚂蚁所走的所有路径都可以在path中找到
dt=zeros(n);%dt表示信息增量
for i=1:m
    ltsp(i)=ca_tsp(n,path(i,:),dij);%每一只蚂蚁所走的路径长度
    for k=1:n-1
        dt(path(i,k),path(i,k+1))=dt(path(i,k),path(i,k+1))+q/ltsp(i);
        dt(path(i,k+1),path(i,k))=dt(path(i,k),path(i,k+1));
    end
    dt(path(i,n),path(i,1))=dt(path(i,n),path(i,1))+q/ltsp(i);
    dt(path(i,1),path(i,n))=dt(path(i,n),path(i,1));
end
[min1 i]=min(ltsp);
if min1<min10
    min10=min1;
    c0=path(i,:);
end
t=r*t+dt;
end
ta=t.^a;
tb=dij.^(-b);
td=ta.*tb;
k=3;
ts(1)=1;
td(:,1)=0;
[my,i]=max(td(1,:));
ts(2)=i;
td(:,i)=0;
while k<=n
    [my,i]=max(td(i,:));
    ts(k)=i;
    td(:,i)=0;
    k=k+1;
end
ts;
ltsp0=ca_tsp(n,ts,dij);
if min10<ltsp0
    ts=c0;
    ltsp0=min10;
end
k=1;
while k<=n
    x1(k)=x(ts(k));
    y1(k)=y(ts(k));
    k=k+1;
end
x1(n+1)=x1(1);
y1(n+1)=y1(1);
line(x1,y1)
hold on
plot(x,y,'o');
ltsp0