模拟退火及旅行商算法

小白刚加入数模编程大军，涉及到一些算法，特此归纳，并加入一些自己的想法，类比此题进行另一题的程序修改。

模拟退火算法的模型 
模拟退火算法可以分解为解空间、目标函数和初始解三部分。 
模拟退火的基本思想: 
(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)， 每个T值的迭代次数L 
(2) 对k=1，……，L做第(3)至第6步： 
(3) 产生新解S′ 
(4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数 
(5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)接受S′作为新的当前解. 
(6) 如果满足终止条件则输出当前解作为最优解，结束程序。 
     终止条件通常取为连续若干个新解都没有被接受时终止算法。 

(7) T逐渐减少，且T->0，然后转第2步。

例  已知敌方100个目标的经度、纬度如下：

经度      纬度

经度      纬度

经度      纬度

经度      纬度

53.7121   15.3046

51.1758    0.0322

46.3253   28.2753

30.3313    6.9348

56.5432   21.4188

10.8198   16.2529

22.7891   23.1045

10.1584   12.4819

20.1050   15.4562

1.9451    0.2057

26.4951   22.1221

31.4847    8.9640

26.2418   18.1760

44.0356   13.5401

28.9836   25.9879

38.4722   20.1731

28.2694   29.0011

32.1910    5.8699

36.4863   29.7284

0.9718   28.1477

8.9586   24.6635

16.5618   23.6143

10.5597   15.1178

50.2111   10.2944

8.1519    9.5325

22.1075   18.5569

0.1215   18.8726

48.2077   16.8889

31.9499   17.6309

0.7732    0.4656

47.4134   23.7783

41.8671    3.5667

43.5474    3.9061

53.3524   26.7256

30.8165   13.4595

27.7133    5.0706

23.9222    7.6306

51.9612   22.8511

12.7938   15.7307

4.9568    8.3669

21.5051   24.0909

15.2548   27.2111

6.2070    5.1442

49.2430   16.7044

17.1168   20.0354

34.1688   22.7571

9.4402    3.9200

11.5812   14.5677

52.1181    0.4088

9.5559   11.4219

24.4509    6.5634

26.7213   28.5667

37.5848   16.8474

35.6619    9.9333

24.4654    3.1644

0.7775    6.9576

14.4703   13.6368

19.8660   15.1224

3.1616    4.2428

18.5245   14.3598

58.6849   27.1485

39.5168   16.9371

56.5089   13.7090

52.5211   15.7957

38.4300    8.4648

51.8181   23.0159

8.9983   23.6440

50.1156   23.7816

13.7909    1.9510

34.0574   23.3960

23.0624    8.4319

19.9857    5.7902

40.8801   14.2978

58.8289   14.5229

18.6635    6.7436

52.8423   27.2880

39.9494   29.5114

47.5099   24.0664

10.1121   27.2662

28.7812   27.6659

8.0831   27.6705

9.1556   14.1304

53.7989    0.2199

33.6490    0.3980

1.3496   16.8359

49.9816    6.0828

19.3635   17.6622

36.9545   23.0265

15.7320   19.5697

11.5118   17.3884

44.0398   16.2635

39.7139   28.4203

6.9909   23.1804

38.3392   19.9950

24.6543   19.6057

36.9980   24.3992

4.1591    3.1853

40.1400   20.3030

23.9876    9.4030

41.1084   27.7149

 我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。

load ij.txt;%加载目标数据ij.txt存放矩阵x=ij(:,1:2:8);%把sj中属于1,3,5,7的列所有元素赋给x,此时,x的行的大小为sj的行的大小,列的大小为4,x=x(:);%x=x(:)的意思是把x拉成一维列向量(应该是列向量)y=ij(:,2:2:8);y=y(:);sj=[x y];%将xy合并成一个矩阵d1=[70,40];%起始点sj=[d1;sj;d1];%首尾添加上出发点和终点sj=sj*pi/180;%角度化为弧度%计算距离d=zeros(102);%先建一个102x102的空的距离矩阵，再把对应距离填进去%A(R cos x1 cos y1 , R sin x1 cos y1 , R sin y1 )%B(R cos x2 cos y2 , R sin x2 cos y2 , R sin y2 )%其中 R = 6370 为地球半径。%A, B 两点的实际距离:%d = R arccos[cos(x1 − x2 ) cos y1 cos y2 + sin y1 sin y2 ]for i=1:101    for j=i+1:102        d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))...            +sin(sj(i,2))*sin(sj(j,2)));    endend%上述的d为上三角矩阵，根据对称性，化为对称矩阵,点击打开链接d矩阵包含各点间的距离信息d=d+d';S0=[];Sum=inf;rand('state',sum(clock));%作用是定义一个随时间变化的初值%matlab里面的随机生成函数基本都是以rand为基函数通过函数关系式得到，%比如 normrnd，unidrnd等，你每次重启matlab后运行已编好的含随机数生成的函数你将得到相同的结果，%比如我的电脑上重启matlab运行unidrnd(100),每次的值都是82，这是因为rand函数的初值都一样，%所以为了避免上述问题经常在程序前运行或加命令rand('state',sum(clock))，这样重启matlab，运行随机数生成值就不同了。for i=1:1000%循环1000次，尽可能让所有排序情况都出现一遍    s=[1,1+randperm(100),102];%s为中间2到101的排列组合    temp=0;    for j=1:101        temp=temp+d(s(j),s(j+1));    end    if temp<Sum       Sum=temp;       S0=s;%保存当前最小值的排序情况    endende=0.1^30;%选定的终止温度 e = 10 −30 ，判断退火过程是否结束at=0.999;%选定的降温系数α进行降温即：T =αTL=20000;T=1;%退火过程for k=1:Lc=2+floor(100*rand(1,2));%floor函数其功能是“向下取整”，或者说“向下舍入”，即取不大于x的最大整数%ceil函数向上取整c=sort(c);% [a,b]=sort(X)是按列从小到大排序，而 [a,b]=sort(X,2)是按行，b为排序情况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)%以概率exp(−df/T)接受新的路径%注意时elseif而不是else if    S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];    Sum=Sum+df;end;T=T*at;%退火if T<e%达到终止温度    break;end;end;%输出路径和路径长度S0,Sum%画出其中一个巡航路径plot(sj(S0,1)/pi*180,sj(S0,2)/pi*180);hold onplot(sj(S0,1)/pi*180,sj(S0,2)/pi*180,'rx');axis([-5,75,-5,45]);


附ij.txt
53.7121   15.304651.1758    0.032246.3253   28.275330.3313    6.934856.5432   21.418810.8198   16.252922.7891   23.104510.1584   12.481920.1050   15.45621.9451    0.205726.4951   22.122131.4847    8.964026.2418   18.176044.0356   13.540128.9836   25.987938.4722   20.173128.2694   29.001132.1910    5.869936.4863   29.72840.9718   28.14778.9586   24.663516.5618   23.614310.5597   15.117850.2111   10.29448.1519    9.532522.1075   18.55690.1215   18.872648.2077   16.888931.9499   17.63090.7732    0.465647.4134   23.778341.8671    3.566743.5474    3.906153.3524   26.725630.8165   13.459527.7133    5.070623.9222    7.630651.9612   22.851112.7938   15.73074.9568    8.366921.5051   24.090915.2548   27.21116.2070    5.144249.2430   16.704417.1168   20.035434.1688   22.75719.4402    3.920011.5812   14.567752.1181    0.40889.5559   11.421924.4509    6.563426.7213   28.566737.5848   16.847435.6619    9.933324.4654    3.16440.7775    6.957614.4703   13.636819.8660   15.12243.1616    4.242818.5245   14.359858.6849   27.148539.5168   16.937156.5089   13.709052.5211   15.795738.4300    8.464851.8181   23.01598.9983   23.644050.1156   23.781613.7909    1.951034.0574   23.396023.0624    8.431919.9857    5.790240.8801   14.297858.8289   14.522918.6635    6.743652.8423   27.288039.9494   29.511447.5099   24.066410.1121   27.266228.7812   27.66598.0831   27.67059.1556   14.130453.7989    0.219933.6490    0.39801.3496   16.835949.9816    6.082819.3635   17.662236.9545   23.026515.7320   19.569711.5118   17.388444.0398   16.263539.7139   28.42036.9909   23.180438.3392   19.995024.6543   19.605736.9980   24.39924.1591    3.185340.1400   20.303023.9876    9.403041.1084   27.7149
输出结果
S0 =1    17    35    93    43    36     3    45    67     2    92    87    83    74    20    30    82    48    72    14    27    10

    84    97    85    77    31    79    60     9    40    18    15    42    50    80    51    98    38     5    75    33    73    13

    24    16    91    41     4    76    69    11    64    65    94    70    19    63    62    26    29    34    66    90    86     8

    39    78    47    57    88    61    49    28    23    58    81    22    25    68     7    71    37    32    12    53    54    89

     6    96    55    44   100    56    21    99   101    52    46    59    95   102

Sum =  4.2699e+004


数据修改小结
1.由于数据矩阵的不同，x，y的取值也有所不同，依据新数据修改此处，y同理。
x=ij(:,1:2:8);
2.其余起点位置，数据数量等据情况修改，此处计算的是飞机绕地球飞行的弧线距离，若用于TSP问题的直线距离，则要修改程序公式
for i=1:101    for j=i+1:102        d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))...            +sin(sj(i,2))*sin(sj(j,2)));    endend
3.博主所用新数据当时大约30个，于是修改了下100，原理也不太懂，有了解的欢迎骚扰呀哈哈
c=2+floor(100*rand(1,2));