tsp.m

在大学生或者研究生数学建模竞赛中旅行商优化问题是一个很重要而又难解决的问题

function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=tsp
%% coryright:QBB
%% input the parameter
prompt={'输入你的名字:','蚂蚁的个数',...
    '最大迭代次数','表征信息素重要程度的参数alpha',...
    '表征启发式因子重要程度的参数beta',...
    '信息素蒸发系数rho','信息素增加强度系数Q'};
dlg_title = '输入参数对话框';
num_lines = [1 40];
def = {'guitar','31','200','1','5','0.1','100'};
ayf = inputdlg(prompt,dlg_title,num_lines,def,'on');
if ~strcmp(ayf{1},'guitar')
    errordlg('姓名不对','Error');
    return
end
temp=cellfun(@str2double,ayf(2:end));
temp=num2cell(temp);
[m,NC_max,Alpha,Beta,Rho,Q]=deal(temp{:});
[FileName,PathName] = uigetfile('*.dat','选择城市坐标的数据文件');
dataFile=strcat(PathName,FileName);
C=load(dataFile);
%% the first step:Initializing the variables
n=size(C,1);% n 表示城市的规模
D=zeros(n); % D 表示完全图的赋权邻接矩阵
for i=1:n
    for j=1:n
        if i~=j
            D(i,j)=sqrt(sumsqr(C(i,2:3)-C(j,2:3)));
        else
            D(i,j)=eps;
        end
    end
end
Eta=1./D;% Eta为启发因子
Tau=ones(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);%各代路线的平均长度
hh= waitbar(0,'Please wait...');
while NC<=NC_max%停止条件之一：达到最大迭代次数
    %% 第二步：将m只蚂蚁放到n个城市上
    Tabu(:,1)=randint(m,1,[1,n]);
    %% 第三步：m只蚂蚁按概率函数选择下一座城市，完成各自的周游
    for j=2:n
        for i=1:m
            visited=Tabu(i,1:(j-1));%已访问的城市
            J=setdiff(1:n,visited);%待访问的城市
            P=J;%待访问城市的选择概率分布
            %下面计算待选城市的概率分布
            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);
            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);
    for i=1:m
        R=Tabu(i,:);
        for j=1:n-1
            L(i)=L(i)+D(R(j),R(j+1));
        end
        L(i)=L(i)+D(R(1),R(n));
    end
    [L_best(NC),pos]=min(L);
    R_best(NC,:)=Tabu(pos,:);
    L_ave(NC)=mean(L);
    NC=NC+1;
    %% 第五步：更新信息素 (Ant-Cycle模型)
    Delta_Tau=zeros(n);
    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);
        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);
    mcc=NC/NC_max;
    waitbar(mcc,hh,['now....',num2str(round(100*mcc)),'%']);
end
close(hh);
%% 第七步：输出结果
[Shortest_Length,Pos]=min(L_best);
Shortest_Route=R_best(Pos,:);
% subplot(1,2,1)
figure(1);
DrawRoute(C,Shortest_Route);
% subplot(1,2,2)
axis equal
figure(2);
plot(L_best);
hold on
plot(L_ave);
%% 子函数
function DrawRoute(C,R)
%====================================================================
%  DrawRoute.m
%  画路线图的子函数
%--------------------------------------------------------------------
%  C    Coordinate        节点坐标，由一个N×2的矩阵存储
%  R    Route             路线
% it is better if use complication
%====================================================================
N=length(R);
scatter(C(:,2),C(:,3),5);
for i=1:size(C,1)
    text(C(i,2)+20,C(i,3)+20,num2str(i));
end
hold on
plot([C(R(1),2),C(R(N),2)],[C(R(1),3),C(R(N),3)])
for ii=2:N
    plot([C(R(ii-1),2),C(R(ii),2)],[C(R(ii-1),3),C(R(ii),3)],'r');
    pause(0.05);
end