tspsa.m

matlab 编写的模拟退火算法用于求解TSP问题

function [f,T]=TSPSA
%% copyright:QBB
%% TSP问题（货郎担问题，旅行商问题）的模拟退火算法通用malab源程序
% f目标最优值,T最优路线,d距离矩阵,t0初始温度,tf结束温度

%% input the parameter
prompt={'输入你的名字:','初时温度','结束温度'};
dlg_title = '输入参数对话框';
num_lines = [1 40];
def = {'guitar','80','20'};
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);
[t0,tf]=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
%% 计算初时路径长度
L=100*n;
t=t0;
pi0=1:n;
min_f=conputeDist(pi0,D);
p_min=pi0;
%% 开始退火
while t>tf
    for k=1:L;
        [d_f,pi_1]=exchange_2(pi0,D);
        r_r=rand;
        if d_f<0
            pi0=pi_1;
        elseif exp(-d_f/t)>r_r
            pi0=pi_1;
        end
    end
    f_temp=conputeDist(pi0,D);
    if min_f>f_temp
        min_f=f_temp;
        p_min=pi0;
    end
    t=0.87*t;
end
f=min_f;
T=p_min;
%% 输出结果
DrawRoute(T,C);
axis equal
%% 子函数
function DrawRoute(R,C)
%====================================================================
%  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
%% 要调用的子程序,用于产生新解
function [d_f,pi1]=exchange_2(pi0,d)
%%
newData=randperm(length(pi0));
u=newData(1);
v=newData(2);
L1=conputeDist(pi0,d);
pi1=pi0;
pi1([u,v])=pi1([v,u]);
L2=conputeDist(pi1,d);
d_f=diff([L1 L2]);
%% 计算距离的子程序
function dist=conputeDist(x,d)
dist=0;x=cat(2,x,x(1));
for i=1:length(x)-1
  dist=dist+d(x(i),x(i+1));
end