Matlab利用模拟退火算法求解旅行商(TSP)问题

简介

先引入一个例题:

旅行商问题（TSP）

假设有一个旅行商人要拜访n个城市，已知这n个城市的坐标，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。

这种问题被归类为NP-hard问题，NP是指非确定性多项式（non-deterministic polynomial，缩写NP）。所谓的非确定性是指，可用一定数量的运算去解决多项式时间内可解决的问题。即有限的时间内无法精确解决的问题。

旅行商问题（TSP）就是一个NP-hard问题

如果采用普通搜索算法，即穷举法，在城市数量比较少（10个及以下）的时候还可以解决，一旦城市超过15个，我们就不可能在短时间内求解，数量更多时，求解时间可能要上亿年，所以穷举法求解旅行商问题是不可行的。

我们今天采用启发式算法之一——模拟退火算法

模拟退火算法求解的结果通常不是全局最优解，但是这个算法可以在有限时间内求解NP-hard问题的近似解，近似解和全局最优解通常差距极小

相关概念

温度：每次冷却后会降温，温度越高，当前的解越不稳定，温度较低时，当前的解很稳定.

候选解：符合题目条件的解，又称可行解

新解：已知当前的解，我们人为规定一个法则，生成一个新的候选解

Metropolis准则：先假设我们的目标是最小化问题，求解最小值问题。如果生成的新解的目标值小于当前解的目标值，则直接接受新解；如果生成的新解的目标值大于当前解的目标函数值，则概率接受新解，接受新解的概率取决于两解的目标值之差以及当前温度，当前的温度越高，接受新解的概率越大，两解的目标值之差越大，接受新解的概率越小

用公式表达就是

                               

 

链长：在每个温度下要进行L次迭代，这个L就是链长

冷却：主要是指数冷却，线性冷却

指数冷却：

线性冷却:

基本思想

算法流程：

初始化各个参数以及初始解

当前温度不符合终止准则时：

        迭代次数k<L:

                生成新解

                根据Metropolis准则判断是否接受新解

                k=k+1

        end

        迭代次数到达时，冷却（线性冷却，指数冷却）

end

输出最终解

算法流程图

实例：模拟退火算法求解TSP（旅行商）问题 

数据集下载(tsp文件可以把后缀名改成txt):

链接：https://pan.baidu.com/s/1_8RMaoOwoYhFdxeyzl8-0Q?pwd=1234 
提取码：1234 
--来自百度网盘超级会员V4的分享

注意:下面的函数同时在当前文件夹下,运行程序TSP.m,才可以运行

主程序TSP.m 

%初始参数：T0(初始温度) T_end(结束温度) L(链长) alpha(降温速率) S1(初始解) 
% 1.计算任意两点间的距离矩阵
% 2.生成新解
% 3.根据Metropolis准则判断是否接受新解
% 4.继续迭代，直到达到了迭代次数
% 5.降温
% 6.重复上面的2~5，直到温度降到最低
% 定义的函数：
% 1.根据每个城市的坐标计算任意距离矩阵函数Distence.m
% 2.绘制路线图Drawpath.m
% 3.计算一个解的距离的函数CalDist.m
% 4.生成新解的函数NewSolution.m
% 5.Metropolis准则函数,传入参数为(df,T_now),返回0或1。1代表接受新解,Metropolis.m
% 6.输出解的函数Disp.m

% 总体思路：给定初始温度，温度大于最低温度的时候，代表计算还在继续。每迭代一次，温度以指数规律减少。每次迭代时，指定链长L是操作的最大次数，我们每次迭代一共生成L次新解

%bayg29是一个包括29个点的城市坐标数据集，每行第一个元素是横坐标，第二个元素是纵坐标
bayg29=[1150	1760
630	1660
40	2090
750	1100
750	2030
1030	2070
1650	650
1490	1630
790	2260
710	1310
840	550
1170	2300
970	1340
510	700
750	900
1280	1200
230	590
460	860
1040	950
590	1390
830	1770
490	500
1840	1240
1260	1500
1280	790
490	2130
1460	142
1260	1910
360	1980];
%citydataset代表当前的(N*2)数据集
citydataset=bayg29;
%模型的初始参数可以随时修改，L是链长，代表同一个温度下迭代的次数
T0=50000;T_end=1e-10;L=400;alpha=0.99;
[Num,~]=size(citydataset);%Num是城市的个数
S1=randperm(Num);%S1代表目前的解
DistMatrix=Distance(citydataset);%距离矩阵，第i行第j列是城市i到城市j的距离
disp(['初始路线的长度是：',num2str(CalDist(S1,DistMatrix))])
T_now=T0;
Drawpath(S1,citydataset)
%Tset=T_now;
Tset=1;p=1;
Calset=CalDist(S1,DistMatrix);
while T_now>T_end
    k=1;
    while k<=L
        S_new=NewSolution(S1);
        df=CalDist(S_new,DistMatrix)-CalDist(S1,DistMatrix);
        if Metropolis(df,T_now)==1 %接受新解
            S1=S_new;
        end
        k=k+1;
    end
    T_now=T_now*alpha;
    %Tset=[Tset,T_now];
    Tset=[Tset,p];p=p+1;
    cal=CalDist(S1,DistMatrix);
    Calset=[Calset,cal];
end
figure
plot(Tset,Calset)
xlabel('迭代次数')
ylabel('当前的总路线长度')
Drawpath(S1,citydataset)
disp('最终的路线是：')
Disp([S1])
disp(['最终路线的长度是：',num2str(CalDist(S1,DistMatrix))])

Diatance.m (生成距离矩阵,距离矩阵的第i行第j列代表第i个城市与第j个城市之间的距离)

function y=Distance(Label)
% Label是城市1~N的坐标
[Num,~]=size(Label);
y=zeros(Num,Num);%坐标矩阵初始化
for i=1:Num
    for j=1:Num
        y(i,j)=sqrt((Label(i,1)-Label(j,1))^2+(Label(i,2)-Label(j,2))^2);
    end
end

Metropolis.m (Metropolis准则)

function L=Metropolis(df,T)
if df<0 %新解比当前解小的时候，一定接受新解
    L=1;
else
    r=rand(1);
    if r<exp(-df/T) %新解大于当前解的时候，概率接受新解，接受概率取决于当前温度和新解与旧解之差
        L=1;
    else
        L=0;
    end
end

CalDist.m (根据当前解和城市的距离矩阵,生成这条路线的距离总和)

function sum=CalDist(s,Matrix)% s是当前解，Matrix是距离矩阵
[~,Num]=size(s);
sum=0;
for i=1:Num-1
    sum=sum+Matrix(s(i),s(i+1));
end
sum=sum+Matrix(Num,1);

NewSolution.m (用随机交换两个城市的方法生成新解)

function s1=NewSolution(s)
[~,Num]=size(s);
R=round(rand(1,2)*(Num-1)+1);
s1=s;
s1(R(1))=s(R(2));
s1(R(2))=s(R(1));

Drawpath.m(已知序列和各个城市的坐标矩阵,绘制路线图) 

function Drawpath(S,coordinate)
[~,Num]=size(S);
ChormMatrix=zeros(2,Num+1);
for i=1:Num
    ChormMatrix(1,i)=coordinate(S(i),1);
    ChormMatrix(2,i)=coordinate(S(i),2);
end
ChormMatrix(1,Num+1)=ChormMatrix(1,1);
ChormMatrix(2,Num+1)=ChormMatrix(2,1);
figure
hold on
plot(ChormMatrix(1,:),ChormMatrix(2,:))

Disp.m (已知一个序列,输出路线)

function Disp(S)
[~,Num]=size(S);
p=num2str(S(1));
for i=2:Num-1
    p=[p,'->',num2str(S(i))];
end
disp(p)

求解结果:

初始解:

 

最终解 

 

>> TSP

初始路线的长度是：27863.5694

最终的路线是：

27->7->25->23->8->24->13->16->19->11->22->17->14->18->15->4->10->20->2->21->1->28->6->12->9->5->26->29

最终路线的长度是：9665.4694

注意:模拟退火算法每次运行的结果都不一定相同,所以可以多试几次,得到一个比较小的解