matlab模拟退火算法（SA）详解（二）旅行商问题（TSP）详解

 

旅行商问题(Traveling Salesman Problem,TSP)代表一类组合优化问题,在物流配送、计算机网络、电子地图、交通疏导、电气布线等方面都有重要的工程和理论价值,引起了许多学者的关注 。TSP 简单描述为:一名商人要到n 个不同的城市去推销商品,每2个城市i和j之间的乐离为d,如何选择一条路径使得商人每个城市走一遍后回到起点,所走的路径最短。 TSP是典型的组合优化问题,并且是一个 NP难题。 TSP 描述起来很简单,早期的研容者使用精确算法求解该问题,常用的方法包括分枝定界法、线性规划法和动态规划法等，但是可能的路径总数随城市数目n是呈指数型增长的,所以当城市数目在100个以上时,一般很难结确地求出其全局最优解。随着人工智能的发展,出现了许多独立于问题的智能优化算法,如蚁群算法、遗传算法、模拟退火、禁忌搜索、神经网络、粒子群优化算法、免疫算法等,通过模拟或解释某些自然现象或过程而得以发展。模拟退火算法具有高效、鲁棒、通用、灵活的优点。将模拟退火算法引入TSP求解,可以避免在求解过程中陷入TSP的局部最优。
算法设计步骤:
(1)TSP问题的解空间和初始解 
 TSP的解空间S是遍访每个城市恰好一次的所有回路,是所有城市排列的集合。TSP问
题的解空间 S 可表示为{1,2,....n}的所有排列的集合,即
S={(C1,C2,..,Cn.) | (C1,C2,...，C)}
其中每一个排列S,表示遍访n个城市的一个路径,Ci=j表示第i次访问城市j。模拟退火算法的最优解和初始状态没有强的依赖关系,故初始解为随机函数生成一个{1,2,...,n}的随机排列作为S。.
(2)目标函数
TSP 问题的目标函数即为访问所有城市的路径总长度,也可称为代价函数,现在 TSP 问题的求解就是通过模拟退火算法求出目标函数C(c1,C2,...,C.)的最小值,相应地,S=(ci,c2,,c:)即为 TSP问题的最优解。 
(3)新解产生
新解的产生对问题的求解非常重要。新解可通过分别或者交替使用以下两种方法来
产生:
二变换法:任选序号u、v(设u<v<n),交换u和v之间的访问顺序。
三变换法:任选序号u、v(设u<v<n),u.v.w(设u<<w),将u和之间的路径插到w之后访问。
(4)目标函数差 
计算变换前的解和变换后目标函数的差值:
△C'=C(s0)-C(s)
(5)Metropolis接受准则 
以新解与当前解的目标函数差定义接受概率,即 
P = (1)1 , △C<0 向好的方向发展

       (2)exp(-△C/T),△C>0 向坏的方向发展

                                                   

关于TSP的数据可以访问这个网站下载https://wwwproxy.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/ 

或者去我的CSDN资源下载，打赏一点积分 https://download.csdn.net/download/viafcccy/11295980

我们使用的是berlin52数据集进行测试，其中包含柏林的真实的52座城市坐标位置。

clear %所有变量全部删除
clc %清除命令行窗口中的数据
a = 0.99; %温度衰减参数
t0 = 97; tf = 3; t = t0; %初始温度 终值
Markov_length = 10000; %Markov链长度
coordinates = [
1 565.0 575.0;
2 25.0 185.0;
3 345.0 750.0;
4 945.0 685.0;
5 845.0 655.0;
6 880.0 660.0;
7 25.0 230.0;
8 525.0 1000.0;
9 580.0 1175.0;
10 650.0 1130.0;
11 1605.0 620.0;
12 1220.0 580.0;
13 1465.0 200.0;
14 1530.0 5.0;
15 845.0 680.0;
16 725.0 370.0;
17 145.0 665.0;
18 415.0 635.0;
19 510.0 875.0;
20 560.0 365.0;
21 300.0 465.0;
22 520.0 585.0;
23 480.0 415.0;
24 835.0 625.0;
25 975.0 580.0;
26 1215.0 245.0;
27 1320.0 315.0;
28 1250.0 400.0;
29 660.0 180.0;
30 410.0 250.0;
31 420.0 555.0;
32 575.0 665.0;
33 1150.0 1160.0;
34 700.0 580.0;
35 685.0 595.0;
36 685.0 610.0;
37 770.0 610.0;
38 795.0 645.0;
39 720.0 635.0;
40 760.0 650.0;
41 475.0 960.0;
42 95.0 260.0;
43 875.0 920.0;
44 700.0 500.0;
45 555.0 815.0;
46 830.0 485.0;
47 1170.0 65.0;
48 830.0 610.0;
49 605.0 625.0;
50 595.0 360.0;
51 1340.0 725.0;
52 1740.0 245.0;
];
coordinates(:,1) = []; %去除第一行的城市编号
amount = size(coordinates,1); %计算城市个数

%通过向量化的方法计算距离矩阵（保存所有两点之间的距离）
coor_x_tmp1 = coordinates(:,1) * ones(1,amount);
coor_x_tmp2 = coor_x_tmp1';
coor_y_tmp1 = coordinates(:,2) * ones(1,amount);
coor_y_tmp2 = coor_y_tmp1';
dist_matrix = sqrt((coor_x_tmp1 - coor_x_tmp2).^2 + (coor_y_tmp1 - coor_y_tmp2).^2);

%产生初始解
sol_new = 1:amount;
%sol代表当前路线 new代表每次产生的新解 current代表当前解 best代表冷却过程中的最优解
E_current = inf;
E_best = inf;
%E代表目标函数的值 也就是距离和 new代表新解的遍历距离 best最优解
sol_current = sol_new;
sol_best = sol_new;
p = 1; %初始化P

figure;
hold on;
box on;
xlim([tf,t0]);
title('优化过程')
xlabel('当前温度')
ylabel('当前最优解')

while t >= tf %当温度大于停止温度时
    for r = 1:Markov_length %Markov链长度
        %产生随机扰动
        if(rand < 0.5) %随机决定而交换还是三交换
            %二交换
            ind1 = 0; ind2 = 0;
            while(ind1 == ind2)
                ind1 = ceil(rand.*amount); %朝正无穷大四舍五入
                ind2 = ceil(rand.*amount);
            end
            tmp1 = sol_new(ind1);
            sol_new(ind1) = sol_new(ind2);
            sol_new(ind2) = tmp1;
        else
            %三交换
            ind1 = 0; ind2 = 0; ind3 = 0;
            while(ind1 == ind2) || (ind1 == ind3) || (ind2 == ind3) || (abs(ind1 - ind2) == 1) %三变换两个点之间没有路径
                ind1 = ceil(rand.*amount);
                ind2 = ceil(rand.*amount);
                ind3 = ceil(rand.*amount);
            end
                %tmp1 = ind1;tmp2 = ind2;tmp3 = ind3;
                %排序ind1<ind2<ind3
                ind = [ind1 ind2 ind3];
                ind_sort = sort(ind);
                ind1 = ind_sort(1);
                ind2 = ind_sort(2);
                ind3 = ind_sort(3);
                %进行交换
                
                tmplist1 = sol_new((ind1 + 1):(ind2 - 1));
                sol_new((ind1 + 1):(ind1 + ind3 - ind2 + 1)) = sol_new((ind2) : (ind3));
                sol_new((ind1 + ind3 - ind2 + 2):ind3) = tmplist1;
        end
        %检查是否满足约束条件
        %计算目标函数函数值（类比于内能）
        E_new = 0;
        for i = 1 : (amount - 1)
            E_new = E_new + dist_matrix(sol_new(i),sol_new(i+1));
        end
        %加上最后一个到第一个的距离
        E_new = E_new + dist_matrix(sol_new(amount),sol_new(1)); 
        if E_new < E_current %当出现更优解时 更新为更优解的路线和目标函数值
            E_current = E_new;
            sol_current = sol_new;
            if E_new < E_best
                pre_E_best = E_best;
                E_best = E_new;
                sol_best = sol_new;
            end
        else
        %一定范围内接受非优解
            if rand < exp(-(E_new - E_current)./t)
                E_current = E_new;
                sol_current = sol_new;
            else
                sol_new = sol_current;
            end
        end
    end
    pre_t = t;
    t = t.* a;%更新控制参数t为原来的a倍
    line([pre_t,t],[pre_E_best,E_best]);
    pause(0.0001);

end

disp('最优解为:')
for i = 1 : amount
fprintf('->%d', sol_best(i));
end
fprintf('\n');
disp('最短距离为:')
disp(E_best);

最优解为:
->26->47->13->14->52->11->51->33->43->10->9->8->41->19->45->32->49->1->22->31->18->3->17->21->42->7->2->30->23->20->50->29->16->46->44->34->35->36->39->40->37->38->48->24->5->15->6->4->25->12->28->27
最短距离为:
   7.5444e+03、

                    

最后，我们使用了模拟退火和遗传算法完成了对于TSP问题的简单求解，这两种算法都是针对全局优化产生的，做出以下小结：

 

全局优化问题主要用于求解组合类或者不适合大规模遍历的问题,其中遗传算法和模拟退火算法最为常用,这两种算法能解决80%以上的全局优化问题。对于具体的全局优化问 题,也不妨两种方法都尝试一下,以便确定哪种算法更优。对MATLAB全局优化的各个求解器,一般选择默认的参数,这些默认的参数适应能力相对较强,也可以调整一下参数,以便研究参数对最优结果的影响。