模拟退火算法（SA）

A，问题描述：已知100个城市的经纬度坐标如下（保存为city.txt）：

53.7121 15.3046 51.1758 0.032246.3253 28.2753 30.3313 6.9348

56.5432 21.4188 10.8198 16.252922.7891 23.1045 10.1584 12.4819

20.1050 15.4562 1.9451 0.205726.4951 22.1221 31.4847 8.9640

26.2418 18.1760 44.0356 13.540128.9836 25.9879 38.4722 20.1731

28.2694 29.0011 32.1910 5.869936.4863 29.7284 0.9718 28.1477

8.9586 24.6635 16.5618 23.614310.5597 15.1178 50.2111 10.2944

8.1519 9.5325 22.1075 18.55690.1215 18.8726 48.2077 16.8889

31.9499 17.6309 0.7732 0.465647.4134 23.7783 41.8671 3.5667

43.5474 3.9061 53.3524 26.725630.8165 13.4595 27.7133 5.0706

23.9222 7.6306 51.9612 22.851112.7938 15.7307 4.9568 8.3669

21.5051 24.0909 15.2548 27.21116.2070 5.1442 49.2430 16.7044

17.1168 20.0354 34.1688 22.75719.4402 3.9200 11.5812 14.5677

52.1181 0.4088 9.5559 11.421924.4509 6.5634 26.7213 28.5667

37.5848 16.8474 35.6619 9.933324.4654 3.1644 0.7775 6.9576

14.4703 13.6368 19.8660 15.12243.1616 4.2428 18.5245 14.3598

58.6849 27.1485 39.5168 16.937156.5089 13.7090 52.5211 15.7957

38.4300 8.4648 51.8181 23.01598.9983 23.6440 50.1156 23.7816

13.7909 1.9510 34.0574 23.396023.0624 8.4319 19.9857 5.7902

40.8801 14.2978 58.8289 14.522918.6635 6.7436 52.8423 27.2880

39.9494 29.5114 47.5099 24.066410.1121 27.2662 28.7812 27.6659

8.0831 27.6705 9.1556 14.130453.7989 0.2199 33.6490 0.3980

1.3496 16.8359 49.9816 6.082819.3635 17.6622 36.9545 23.0265

15.7320 19.5697 11.5118 17.388444.0398 16.2635 39.7139 28.4203

6.9909 23.1804 38.3392 19.995024.6543 19.6057 36.9980 24.3992

4.1591 3.1853 40.1400 20.303023.9876 9.4030 41.1084 27.7149

小明当前位置经纬度坐标为h(70,40)，小明欲前往这100个城市旅游。如何为小明设计一种方案使得小明从h点出发，游览完每一个城市，最后再回到h位置处后整个旅游路线最短。显然，这是一个旅行商问题（TSP），旅行商问题是一个经典的NP难问题。下面要介绍的模拟退火算法可以很好的解决类似这样的旅行商问题。

B，模拟退火算法：模拟退火算法的目标是求NP难组合优化问题的全局最优解。该方法 来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部 粒子 随 温升 变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到 平衡态 ，最后在 常温 时达到基态，内能减为最小。 模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合 概率 突跳特性在 解空间 中随机寻找 目标函数 的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。

C，算法处理过程：这是一个旅行商问题。我们依次给始发地点编号为1，旅游城市依次编号为2，3，…，101，最后我方始发地再重复编号为102（这样便于程序中计算）。则问题是求一个从点1 出发，走遍所有中间点，到达点102的一个最短路径。上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离（具体公式使用经纬度计算公式）。求解的模拟退火算法描述如下：

（1）解空间：解空间S 可表为{1,2,L,101,102 }的所有固定起点和终点的循环排列集合，本文中我们使用 MonteCarlo 方法求得一个较好的初始解。

（2）目标函数：此时的目标函数为侦察所有目标的路径长度或称代价函数：

（3）新解的产生

① 2 变换法：任选序号 u,v（u < v）交换u与v之间的顺序，此时的新路径如下，

② 3 变换法：任选序号u,v 和 w ，将u 和v 之间的路径插到 w之后，对应的新路径为（设u < v < w），

（4）代价函数

对于2变换法，路径差可以表示为

（5）接受准则

如果Δf < 0，则接受新的路径。否则，以概率exp(−Δf /T) 接受新的路径，即若exp(−Δf /T)大于 0到1之间的随机数则接受。

（6）降温

利用选定的降温系数α 进行降温即：T ←αT ，得到新的温度，这里我们取α = 0.999。

（7）结束条件

用选定的终止温度e = 10−30，判断退火过程是否结束。若T < e，算法结束，输出当前状态。

D，算法实现：

% 模拟退火算法

clc,clear;

load city.txt;

x = sj(:,1:2:8);x = x(:);

y = sj(:,2:2:8);y = y(:);

sj = [x,y];

d1 = [70,40];

sj = [d1;sj;d1];

sj = sj * pi /180;

%距离矩阵d

d = zeros(102);

for i = 1 : 101

    for j = i + 1 :102

       temp = cos(sj(i,1) - sj(j,1)) * cos(sj(i,2)) *cos(sj(j,2)) + sin(sj(i,2) * sin(sj(j,2)));

       d(i,j) = 6370 * acos(temp);

    end

end

d = d + d';

%产生较好的初始解

S0 = [];Sum = inf;

rand('state',sum(clock));

for j = 1 : 1000

    S = [11+randperm(100) 102];

    temp =0;

    for i = 1 :101

       temp = temp + d(S(i), S(i +1));

    end

    if temp <Sum

       S0 = S;Sum = temp;

    end

end

e = 0.1^30; L = 20000; at = 0.999; T = 1;

%下面模拟退火过程

while T > e

    for k = 1 :L

       c = 2 + floor(100 * rand(1,2));

       c = sort(c);

       c1 = c(1); c2 = c(2);

       %计算路径差

       df = (d(S0(c1 - 1),S0(c2)) + d(S0(c1),S0(c2+1)))- (d(S0(c1 - 1),S0(c1)) + d(S0(c2),S0(c2+1)));

       %接受准则

       if df < 0

           S0 =[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];

           Sum = Sum+ df;

       elseif exp(-df/T) > rand(1)

           S0 =[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];

           Sum = Sum+ df;

       end 

    end

    T = T *at;

end

E，运行结果

如下代码画出小明的旅游轨迹：

load city.txt;

x = sj(:,1:2:8);x = x(:);

y = sj(:,2:2:8);y = y(:);

sj = [x,y];

d1 = [70,40];

sj = [d1;sj;d1];

%下面path即使用模拟退火算法计算出来的最优路径 1和102表示始发点

path = [1     17    35     93     46    59     36     3    92     2     45    67     82     48    72     14     27    10     84     18    40     15     9    5     60     79    77 31     97     85    65     64     11    76     69     94    70     19     63    39     8     86    62     26     29    34     66     90    23     58     81    22     7     25 68    71     37     32    12     13     24    49     57     47    78     88     28    61     16     91    41     4     73    33     75     53    54     55     96    89     6 56   101     52     99    21   100    38     44     98    80     51     50    42     20     30    74     83     87    43     95   102];

n = 102;

  hold on;

    for i=1:n-1

        plot(sj(path(i),1),sj(path(i),2),'r*'); 

        line([sj(path(i),1)sj(path(i+1),1)],[sj(path(i),2) sj(path(i+1),2)]);  %只连线当前城市和下家城市      

    end

    plot(sj(path(n),1),sj(path(n),2),'r*');  

    line([sj(path(n),1)sj(path(1),1)],[sj(path(n),2) sj(path(1),2)]);     %最后一家城市连线第一家城市

hold off;

其中横坐标表示城市的经度，纵坐标表示城市的维度，每个标记的点代表一个城市。当小明按图中路线游览每一个城市时，所经历的路程较少。