matlab实现模拟退火解决TSP问题

学习内容：

使用MATLAB利用模拟退火算法实现TSP问题的解决。TSP问题，已知城市的坐标或者其相距的距离，找出一条能够从起点出发，走遍所有城市，并回到起点的路程最短的路径。可以输入城市坐标或者距离矩阵进行求解。

学习产出：

源代码以及注释如下，可以根据注释看懂代码编写思路以及作用。

一、模拟退火算法主函数

clear all;
clc;
close all;
%%导入数据，初始化种群
%导入城市坐标
X = [16.47 96.10
     16.47 94.44
     20.09 92.54
     22.39 93.37
     25.23 97.24
     22.00 96.05
     20.47 97.02
     17.20 96.29
     16.30 97.38
     14.05 98.12
     16.53 97.38
     21.52 95.59
     19.41 97.13
     20.09 92.55];
T0 = 1000;                       %设置初始温度
Tend = 1e-3;                     %设置终止温度
L = 200;                         %在每一个温度下的链数
q = 0.9;                         %降温因数
D = Distance(X);                 %距离矩阵
N = size(D,1);                   %城市个数
S1 = randperm(N);                %初始解
DrawPath(S1,X);                  %画出初始解的路线图
disp('初始化种群的一个随机解:')
OutputPath(S1);                    
Rlength = PathLength(D,S1);
disp(['总距离',num2str(Rlength)]);
%计算需要降温的总次数
Time = Timecounter(T0,Tend,L);
count = 0;
Obj = zeros(Time,1);               %初始化目标值记录矩阵
track = zeros(Time,N);             %初始化最优值路径记录矩阵
%% 开始迭代
while T0 > Tend
    count = count + 1;
    temp = zeros(L,N + 1);
    for k = 1:L
        S2 = NewAnswer(S1);             %产生新的解
        [S1,R] = Metropolis(S1,S2,D,T0);
        temp(k,:) = [S1 R];             %记录产生的路径以及最优值
    end
    [d0,index] = min(temp(:,end));      %开始更新
    if count == 1 || d0 < Obj(count - 1)
        Obj(count) = d0;
    else
        Obj(count) = Obj(count - 1);
    end
    track(count,:) = temp(index,1:end - 1); %记录最优值的路径
    T0 = q*T0;
end
%% 绘图
figure
plot(1:count,Obj);
xlabel('迭代次数');
ylabel('距离');
title('优化过程');
DrawPath(track(end,:),X)

disp('最优值')
S = track(end,:);
OutputPath(S);
disp(['总距离：' num2str(PathLength(D,S))]);

二、距离矩阵计算函数

function D = Distance(a)
%%求距离矩阵
%输入：
%a       点的位置
%输出：
%D       距离矩阵
row = size(a,1);
D = zeros(row,row);
for i = 1:row
    for j = 1:row
        D(i,j) = ((a(i,1) - a(j,1))^2 + (a(i,2) - a(j,2))^2)^0.5;
        D(j,i) = D(i,j);
    end
end
end

三、路线图显示函数

function DrawPath(Chrom,X)
%%画出路线图
%输入：
%Chrom      选择的个体
%X          各个城市的坐标
R = [Chrom(1,:) Chrom(1,1)];
figure;
hold on;
plot(X(:,1),X(:,2),'o','color',[0.5 0.5 0.5]);                %画出所有的城市
plot(X(Chrom(1,1),1),X(Chrom(1,1),2),'rv','MarkerSize',20);   %画出起点
for i = 1:size(X,1)
    text(X(i,1) + 0.05,X(i,2) + 0.05,num2str(i),'color',[1,0,0]);
end
%%画箭头
A = X(R,:);
row = size(A,1);
for i = 2:row
    [arrowx,arrowy] = dsxy2figxy(gca,A(i - 1:i,1),A(i - 1:i,2));
    annotation('textarrow',arrowx,arrowy,'HeadWidth',8,'color',[0,0,1]);
end
hold off
%添加说明
xlabel('横坐标')
ylabel('纵坐标')
title('轨迹图');
box on;
end

四、箭头显示函数

function varargout = dsxy2figxy(varargin)
%%本函数用于画箭头 
    if length(varargin{1}) == 1 && ishandle(varargin{1}) ...
                                && strcmp(get(varargin{1},'type'),'axes')   
        hAx = varargin{1};
        varargin = varargin(2:end);
    else
        hAx = gca;
    end
    if length(varargin) == 1
        pos = varargin{1};
    else
        [x,y] = deal(varargin{:});
    end
    axun = get(hAx,'Units');
    set(hAx,'Units','normalized'); 
    axpos = get(hAx,'Position');
    axlim = axis(hAx);
    axwidth = diff(axlim(1:2));
    axheight = diff(axlim(3:4));
    if exist('x','var')
        varargout{1} = (x - axlim(1)) * axpos(3) / axwidth + axpos(1);
        varargout{2} = (y - axlim(3)) * axpos(4) / axheight + axpos(2);
    else
        pos(1) = (pos(1) - axlim(1)) / axwidth * axpos(3) + axpos(1);
        pos(2) = (pos(2) - axlim(3)) / axheight * axpos(4) + axpos(2);
        pos(3) = pos(3) * axpos(3) / axwidth;
        pos(4) = pos(4) * axpos(4 )/ axheight;
        varargout{1} = pos;
    end
    set(hAx,'Units',axun)
end

五、逆转排列函数

function S = intercross(S2)
%%个体选中某一段进行反向排序
%例如32311变为31321，既选择了中间三个数，将其反向排列
%输入：
%S2       刚输入的序列
%输出：
%S2       新产生的序列
L = length(S2);
r1 = randsrc(1,1,[1:L]);
r2 = randsrc(1,1,[1:L]);
while r1 == r2
    r1 = randsrc(1,1,[1:L]);
    r2 = randsrc(1,1,[1:L]);
end
s = min([r1 r2]);
e = max([r1 r2]);
    for i = s:e
       S2(:,s:e) = S2(:,e:-1:s);
       S = S2;      
    end
end

六、模拟退火算法选择函数

function [S,R] = Metropolis(S1,S2,D,T)
%该函数用于判读是否选择新的路径
%输入：
%S1              原来的路径
%S2              新的路径
%D               城市的距离矩阵
%T               当前的温度
%输出：
%S               选择后的路径
%R               选择的路径的长度
R1 = PathLength(D,S1);
R2 = PathLength(D,S2);
dc = R2 - R1;
if dc < 0
    S = S2;
    R = R2;
elseif exp(-dc/T) >= rand
    S = S2;
    R = R2;
else
    S = S1;
    R = R1;
end
end

七、新路径生成函数

function S2 = NewAnswer(S1)
%该函数用于交换两个位置产生新的的路径
%输入：
%S1        原来的路径
%输出：
%S2        新的路径
N = length(S1);
S2 = S1;
a = round(rand(1,2)*(N - 1) + 1);
W = S2(a(1));
S2(a(1)) = S1(a(2));
S2(a(2)) = W;
S2 = intercross(S2);   %将对S2进行逆转交换
end

八、命令行输出路径函数

function p = OutputPath(R)
%%在命令窗口给出路线图
R = [R,R(1)];                %从起点回到起点
N = length(R);
p = num2str(R(1));
for i = 2:N
    p = [p,'->',num2str(R(i))];
end
disp(p)

九、计算路径长度函数

function len  = PathLength(D,Chrom)
%%计算每个个体路线的总长度
%输入：
%D          距离矩阵
%Chrom      个体
%输出：
%len        长度
[row,col] = size(D);
NIND = size(Chrom,1);
len = zeros(NIND,1);
for i = 1:NIND
    p = [Chrom(i,:) Chrom(i,1)];
%     i1 = p(1:end - 1);                     %将坐标进行转换，算出距离长度
%     i2 = p(2:end);
%     len(i,1) = sum(D((i1 - 1)*col + i2));
    lentemp = 0;
    for j = 1:length(p) - 1
        lentemp = lentemp + D(p(j),p(j+1));
    end
    len(i,1) = lentemp;
end
end

十、迭代次数计算函数

function Time = Timecounter(T0,Tend,L)
%该函数用于求解在特定条件下需要迭代的次数
%输入：
%T0        初始温度
%Tend      最低温度
%L         温度降低因子
%输出：
%Time      最多需要迭代的次数
syms x
eqns = T0.*(L).^x == Tend;
Time = solve(eqns,x);
Time = ceil(double(Time));
end