模拟退火算法

一、模拟退火算法简介

      

模拟退火算法（Simulated  Annealing, SA）是一种通用的随机搜索算法，是对局部搜索算法的扩展。与一般局部搜索算法不同，SA以一定的概率选择邻域中目标值相对较小的状态，是一种理论上的全局最优算法。

模拟退火算法是源于对热力学中退火过程的模拟，在某一给定初温下，通过缓慢下降温度参数，使算法能够在多项式时间内给出一个近似最优解。

二、模拟退火算法与金属退火过程对比

   

物理退火过程：

加温过程----固态溶解为液态的过程，分子的分布从有序的晶体态转化为无序的液态，消除原有的非均匀态。

等温过程----在退火中，需要保证系统在每一个温度下都达到充分的平衡状态。

冷却过程----液态凝固成固态晶体的过程。系统能量逐渐下降，从而得到低能有序的晶体结构。

模拟退火过程：

设定初始高温，相当于物理退火的加温过程。初始温度要足够高，在实际应用中，要根据以往的经验，通过反复实验来确定T0的值。

热平衡达到，相当于物理退火的等温过程。是指在一个给定温度下，SA用特殊的抽样策略进行随机搜索，最终达到平衡状态的过程。这是SA算法的内循环过程。

降温函数，相当于物理退火的冷却过程。用来控制温度的下降方式，这是SA算法的外循环过程。常用的降温函数有Tk+1=Tk-DT，Tk+1=Tk*r，其中r∈(0.95,0.99)。

  1953年，Metropolis等提出了一种重要性采样法(Metropolis准则)，即以概率接受新状态。

   Metropolis准则

   在温度t，由当前状态i产生新状态j，两者的能量分别为Ei 和 Ej ,若Ei >Ej ，则接受新状态j为当前状态；否则，以一定概率：

         

 来接受状态j，其中k为Boltzmann常数。当这种过程多次重复后，系统将趋于能量较低的平衡态。

  三、模拟退火算法流程图

  

  四、模拟退火算法求解TSP问题

   

 

  参数的设定

 •初始解:随机选定一条路径，如[A,B,C,D,E,A]

 •初始温度T0=100

 •终止温度Tf=20

 •降温函数Tk+1=Tk-DT，DT=40

 •热平衡达到:内循环次数n=3

代码实现如下:

%======================================================================================

   %%%一个旅行商人要拜访全国31个省会城市，需要选择最短的路径%%%%


%%%模拟退火算法解决TSP问题%%%%%%%

clear all; %清除所有变量
close all; %清图
clc ;      %清屏
C=[
1304 2312;
3639 1315;
4177 2244;
3712 1399;
3488 1535;
3326 1556;
3238 1229;
4196 1004;
4312 790;
4386 570;
3007 1970;
2562 1756;
2788 1491;
2381 1676;
1332 695;
3715 1678;
3918 2179;
4061 2370;
3780 2212;
3676 2578;
4029 2838;
4263 2931;
3429 1908;
3507 2367;
3394 2643;
3439 3201;
2935 3240;
3140 3550;
2545 2357;
2778 2826;
2370 2975
];%31个省会坐标


n=size(C,1); %TSP问题的规模，即城市数目
T=100*n;     %初始温度
L=100;       %马尔科夫链的长度
K=0.99;      %衰减参数


%%%城市坐标结构体%%%%%%%
city=struct([]);


for i=1:n
    city(i).x=C(i,1);
    city(i).y=C(i,2);
end


l=1;        %统计迭代次数
len(l)=func5(city,n); %每次迭代后路线的长度


figure(1);

while T>0.001
    %%%%%%%%%%%%%%%多次迭代扰动，温度降低前多次试验%%%%%%%%
    for i=1:L
    %%%%%%%%%%%%%%%计算原路线总距离%%%%%%%%%
    len1=func5(city,n);
    %%%%%%%%%%%%%%%产生随机扰动%%%%%%%%%
    %%%%%%%%%%%%%%%随机置换两个不同的城市的坐标%%%%%%%%%
    p1=floor(1+n*rand);
    p2=floor(1+n*rand);
    while p1==p2
        p1=floor(1+n*rand);
        p2=floor(1+n*rand);
    end
    tmp_city=city;
    %%交换元素
    tmp=tmp_city(p1);
    tmp_city(p1)=tmp_city(p2);
    tmp_city(p2)=tmp;
    %%%%%%%%%%%%%%%计算新路线总距离%%%%%%%%%
    len2=func5(tmp_city,n);
    %%%%%%%%%%%%%%%新老距离的差值，相当于能量%%%%%%%%%
    delta_e=len2-len1;
    %%%%%%%%%%%%%%%新路线好于旧路线，用新路线替代旧路线%%%%%%%%%
    if delta_e<0
        city=tmp_city;
    else
        %%%%%%%%%%%%%%%以一定概率选择是否接受%%%%%%%%%
        if exp(-delta_e/T)>rand()
            city=tmp_city;
        end
    end
    end
     l=l+1;          %迭代次数加1
     
     %%%%%%%%%%%%%%%计算新路线的距离%%%%%%%%%
     len(l)=func5(city,n);
     %%%%%%%%%%%%%%%温度不断下降%%%%%%%%%
     T=T*K;
     for i=1:n-1
         plot([city(i).x,city(i+1).x],[city(i).y,city(i+1).y],'bo-');
         hold on;
     end
      plot([city(n).x,city(1).x],[city(n).y,city(1).y],'ro-');
      
      title(['优化最短距离:',num2str(len(l))]);
         hold off;
         pause(0.005);
end
  
  figure(2)
  plot(len);
  xlabel('迭代次数');
  ylabel('目标函数值');
  title('适应度进化曲线');

%计算距离的函数
function len=func5(city,n)
 len=0;
 for i=1:n-1
    len=len+sqrt((city(i).x-city(i+1).x)^2+(city(i).y-city(i+1).y)^2);
 end
   len=len+sqrt((city(n).x-city(1).x)^2+(city(n).y-city(1).y)^2);
end

%=========================================================================================