我们在模拟退火（SA）算法求解旅行商（TSP）问题这篇推文讲解了SA求TSP问题的大致思路，今天为各位详细讲解一下这篇推文中的MATLAB代码（PS：我目前使用的是MATLAB R2019b）。公众号后台回复 SATSP 即可提取代码。


01 | SA_TSP主函数


%%      @作者：随心390%      @微信公众号：优化算法交流地%clc;clear;close all;%% 输入数据x=[38.24,39.57,40.56,36.26,33.48,37.56,38.42,37.52,41.23,41.17,36.08,38.47,38.15,37.51,35.49,39.36,38.09,36.09,40.44,40.33,40.37,37.57];y=[20.42,26.15,25.32,23.12,10.54,12.19,13.11,20.44,9.100,13.05,-5.210,15.13,15.35,15.17,14.32,19.56,24.36,23,13.57,14.15,14.23,22.56];vertexs=[x;y]';n=length(x);                    %城市数目h=pdist(vertexs);dist=squareform(h);             %距离矩阵%% 参数MaxOutIter=300;                     %外层循环最大迭代次数MaxInIter=15;                       %里层循环最大迭代次数T0=0.025;                           %初始温度alpha=0.99;                         %冷却因子pSwap=0.2;                          %选择交换结构的概率pReversion=0.5;                     %选择逆转结构的概率pInsertion=1-pSwap-pReversion;      %选择插入结构的概率%% 初始化currRoute=randperm(n);              %随机构造初始解currL=RouteLength(currRoute,dist);  %初始解总距离bestRoute=currRoute;                %初始将初始解赋值给全局最优解bestL=currL;                        %初始将初始解总距离赋值给全局最优解总距离BestLOutIter=zeros(MaxOutIter,1);   %记录外层循环的每一代最优解的总距离T=T0;                               %温度初始化%% 模拟退火for outIter=1:MaxOutIter    for inIter=1:MaxInIter        newRoute=Neighbor(currRoute,pSwap,pReversion,pInsertion);       %经过邻域结构后产生的新的路线        newL=RouteLength(newRoute,dist);                                %该路线的总距离        %如果新路线比当前路线更好，则更新当前路线，以及当前路线总距离        if newL<=currL             currRoute=newRoute;             currL=newL;        else             %如果新路线不如当前路线好，则采用退火准则，以一定概率接受新路线            delta=(newL-currL)/currL;           %计算新路线与当前路线总距离相差的百分比            P=exp(-delta/T);                    %计算接受新路线的概率            %如果0~1的随机数小于P，则接受新路线，并更新当前路线，以及当前路线总距离            if rand<=P                currRoute=newRoute;                 currL=newL;            end        end        %将当前路线与全局最优路线进行比较，如果当前路线更好，则更新全局最优路线，以及全局最优路线总距离        if currL<=bestL            bestRoute=currRoute;            bestL=currL;        end    end    %记录外层循环每次迭代的全局最优路线的总距离    BestLOutIter(outIter)=bestL;    %显示外层循环每次迭代的信全局最优路线的总距离    disp(['第' num2str(outIter) '次迭代：全局最优路线总距离 = ' num2str(BestLOutIter(outIter))]);    %更新当前温度    T=alpha*T;    %画出外层循环每次迭代的全局最优路线图    figure(1);    PlotRoute(bestRoute,x,y)    pause(0.01);end%% 打印外层循环每次迭代的全局最优路线的总距离变化趋势图figure;plot(BestLOutIter,'LineWidth',1);xlabel('迭代次数');ylabel('距离');

02 | RouteLength函数RouteLength函数是求一条路线的总距离。


%%      @作者：随心390%      @微信公众号：优化算法交流地%%输入route：               路线%输入dist：                距离矩阵%输出L：                   该路线总距离function L=RouteLength(route,dist)    n=length(route);    route=[route route(1)];    L=0;    for k=1:n         i=route(k);        j=route(k+1);         L=L+dist(i,j);     endend
03 | Neighbor函数

Neighbor函数是求出路线1经过邻域结构后的路线2。


%%      @作者：随心390%      @微信公众号：优化算法交流地%%输入route1：路线1%输入pSwap：选择交换结构的概率%输入pReversion：选择逆转结构的概率%输入pInsertion：选择插入结构的概率%输出route2：经过邻域变换后的路线2function route2=Neighbor(route1,pSwap,pReversion,pInsertion)index=Roulette(pSwap,pReversion,pInsertion);if index==1    %交换结构    route2=Swap(route1);elseif index==2    %逆转结构    route2=Reversion(route1);else    %插入结构    route2=Insertion(route1);endend
04 | Roulette函数

Roulette函数是轮盘赌选择，输出结果为序号。


%%      @作者：随心390%      @微信公众号：优化算法交流地%%输入pSwap：选择交换结构的概率%输入pReversion：选择逆转结构的概率%输入pInsertion：选择插入结构的概率%输出index：最终选择哪一个邻域结构，即序号：1 2 3function index=Roulette(pSwap,pReversion,pInsertion)p=[pSwap pReversion pInsertion];r=rand;c=cumsum(p);index=find(r<=c,1,'first');end
05 | Swap函数


%%      @作者：随心390%      @微信公众号：优化算法交流地%
%比如说有6个城市，当前解为123456，我们随机选择两个位置，然后将这两个位置上的元素进行交换。%比如说，交换2和5两个位置上的元素，则交换后的解为153426。
%输入route1：路线1%输出route2：经过交换结构变换后的路线2function route2=Swap(route1)n=length(route1);seq=randperm(n);I=seq(1:2);i1=I(1);i2=I(2);route2=route1;route2([i1 i2])=route1([i2 i1]);end
06 | Reversion函数


%%      @作者：随心390%      @微信公众号：优化算法交流地%
%有6个城市，当前解为123456，我们随机选择两个位置，然后将这两个位置之间的元素进行逆序排列。%比如说，逆转2和5之间的所有元素，则逆转后的解为154326。
%输入route1：路线1%输出route2：经过逆转结构变换后的路线2function route2=Reversion(route1)n=length(route1);seq=randperm(n);I=seq(1:2);i1=min(I);i2=max(I);route2=route1;route2(i1:i2)=route1(i2:-1:i1);end
07 | Insertion函数


%%      @作者：随心390%      @微信公众号：优化算法交流地%
%有6个城市，当前解为123456，我们随机选择两个位置，然后将这第一个位置上的元素插入到第二个元素后面。%比如说，第一个选择5这个位置，第二个选择2这个位置，则插入后的解为125346。
%输入route1：路线1%输出route2：经过插入结构变换后的路线2function route2=Insertion(route1)n=length(route1);seq=randperm(n);I=seq(1:2);i1=I(1);i2=I(2);if i1<i2    route2=route1([1:i1-1 i1+1:i2 i1 i2+1:end]);else    route2=route1([1:i2 i1 i2+1:i1-1 i1+1:end]);endend
08 | PlotRoute函数


%%      @作者：随心390%      @微信公众号：优化算法交流地%%输入route：路线%输入x,y：x,y坐标function PlotRoute(route,x,y)route=[route route(1)];plot(x(route),y(route),'k-o','MarkerSize',10,'MarkerFaceColor','w','LineWidth',1.5);xlabel('x');ylabel('y');end