模拟退火算法解决TSP问题

模拟退火法

模拟退火法简单解释：为了找出地球上最高的山，一群兔子们开始想办法。首先兔子们用酒将自己灌醉了，它们随机地跳了很长的时间。在这期

间，它们可能走向高处，也可能踏入平地。但是，随着时间的流逝，它们渐渐清醒了并朝着最高的方向跳去。最后就达到了珠穆朗玛峰。

所谓的退火是指将固体加热到足够高的温度，使分子呈随机排列状态，然后逐步降温使之冷却，最后分子以低能状态排列，固体达到某种稳定状

态。

重要理解：

假设材料在状态i之下的能量为E(i),那么材料在温度T时从状态i进入状态j就遵循如下规律：

如果 E(j)<= E(i),接受该状态被转移

如果 E(j)> E(i)，则状态转移以如下的概率被接受：

     exp{[E(i)-E(j)]/KT}

 其中，K 是物理学中的波尔兹曼常数，T 是材料温度。

算法关键参数和操作的设定：

  从算法的流程上看，模拟退火算法包括三函数两准则，即状态产生函数、状态接受函数、温度更新函数、内循环终止准则和外循环终止准则，这

些环节的设计将决定模拟退火算法的优化性能。此外，初温的选择对模拟退火算法性能也有很大影响。

状态产生函数：

  原则：设计状态产生函数(领域函数)的出发点应该是尽可能保证产生的候选解遍布全部的解空间。通常，状态产生函数由两部分组成，即产生候

选解的方式和候选解产生的概率分布。

  方法：在当前状态的领域结构内以一定概率方式(均匀分布、正态分布、指数分布等)产生


状态接受函数：

  原则：函数一般以概率的方式给出，不同接受函数的差别主要在于接受概率的形式不同。设计状态接受概率，应该遵循以下原则：

  （1）在固定温度下，接受使目标函数下降的候选解的概率要大于使目标函数上升的候选解概率；

  （2）随温度的下降，接受使目标函数上升的解的概率要逐渐减少；

  （3）当温度趋于零时，只能接受目标函数下降的解。

  方法：状态接受函数的引入是模拟退火算法实现全局搜索的最关键的因素，模拟退火算法中通常用 min[1,exp(-△C/t)]作为状态接受函数。

初始温度：

  初始温度、温度更新函数、内循环终止准则和外循环终止准则通常被称为退火历程。

   原则：通过理论分析可以得到初温的解析式，但解决实际问题时难以得到精确的参数；实际应用时往往要让初温充分大。实验表明：初温越大，

获得高质量解的机率越大，但花费较多的计算时间。

   方法：

  （1）均匀抽样一组状态，以各状态目标值的方差为初温；

  （2）随机产生一组状态，确定两两状态间的最大目标差值，根据差值，利用一定的函数确定初温，譬如 t0 =-△/lnPr,其中Pr 为初始接受函数；

  （3）利用经验公式。

温度更新函数：

  温度更细函数，即温度的下降方式，用于在外循环中修改温度值。

内循环终止准则：

  常用的Metropolis 抽样准则：

   （1）检验目标函数的均值是否稳定；

   （2）连续若干步的目标值变化较小;

   （3）按一定的步数抽样。

外循环终止准则

（1）设置终止温度的阈值；

（2）设置外循环迭代次数；

（3）算法搜索到的最优值连续若干步保持不变；

（4）概率分析方法。

模拟退火的基本思想：

(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步。

以下是实例程序：

[plain] view plain copy print ?

1. main.m 文件  
2. zuobiao=[0.37 0.75 0.45 0.76 0.71 0.07 0.42 0.59 0.32 0.6 0.3 0.67 0.62 0.67 0.20 ...  
3.     0.35 0.27 0.94 0.82 0.37 0.61 0.42 0.6 0.39 0.53 0.4 0.63 0.5 0.98 0.68;  
4.     0.91 0.87 0.85 0.75 0.72 0.74 0.71 0.69 0.64 0.64 0.59 0.59 0.55 0.55 0.5...  
5.     0.45 0.43 0.42 0.38 0.27 0.26 0.25 0.23 0.19 0.19 0.13 0.08 0.04 0.02 0.85];  
6. plot(zuobiao(1,:),zuobiao(2,:),'g*'),hold on  
7. plot(zuobiao(1,:),zuobiao(2,:))  
8. length=max(size(zuobiao));  
9. %求初始距离..  
10.   
11. zhixu=randperm(length);  %随机生成一个路线经过点的顺序  
12. temp=zuobiao(1,:);  
13. newzuobiao(1,:)=temp(zhixu);  
14. temp=zuobiao(2,:);  
15. newzuobiao(2,:)=temp(zhixu);  
16. newzuobiao;  
17. f=juli(newzuobiao);  
18.   
19.   
20. %参数定义区--------------------------------------  
21. %初始温度为10000  
22. tmax=100;  
23. tmin=0.001;  
24. %温度下降速率  
25. down=0.95;  
26.   
27.   
28. %退火算法的函数..  
29. figure  
30. t=tmax;  
31. while t>tmin    
32.     for n=1:500  
33.         newzuobiao=newpath(zuobiao,length);  
34.         newf=juli(newzuobiao);  
35.         if newf<f  
36.            zuobiao=newzuobiao;  
37.            f=newf;  
38.         elseif   rand<exp(-(newf-f)/t)  
39.            zuobiao=newzuobiao;  
40.            f=newf;  
41.         end  
42.     end  
43.     huatu=[zuobiao,zuobiao(:,1)];  
44.     plot(huatu(1,:),huatu(2,:)),hold on  
45.     plot(huatu(1,:),huatu(2,:),'ro'),hold off  
46.     pause(0.00001);  
47.     t=t*down;  
48.     f;  
49. end  
50.   
51. newpath.m 文件(随机的交换两个点的坐标)  
52. function zuobiao=newpath(zuobiao,length)  
53. %随机交换两个点的坐标..  
54. a=ceil(rand(1,2)*length);  
55. qian=a(1);  
56. hou=a(2);  
57. temp=zuobiao(:,qian);  
58. zuobiao(:,qian)=zuobiao(:,hou);  
59. zuobiao(:,hou)=temp;  
60.   
61. juli.m 文件(计算按某种方式围成的圈的总路程)  
62. function lucheng=juli(zuobiao)  
63. length=max(size(zuobiao));  
64. s=0;  
65. for i=2:length  
66.     s=s+sqrt(sum((zuobiao(:,i)-zuobiao(:,i-1)).^2));  
67. end  
68. if length~=2  
69.     s=s+sqrt(sum((zuobiao(:,1)-zuobiao(:,length)).^2));  
70. end  
71.     lucheng=s;