模拟退火算法MATLAB实现

介绍

 算法试图随着控制参数T的降低，使目标函 数值f（内能E）也逐渐降低，直至趋于全局最 小值（退火中低温时的最低能量状态），算法 工作过程就像固体退火过程一样。

Metropolis准则——–以概率接受新状态

若在温度T，当前状态i (解1)→ 新状态 j(解2)

若E_j(目标函数f_2)<E_i(f_1)，则接受 j 为当前状态；

若E_j>E_i ，概率P=e^−E_j−E_i/KT大于[0,1)区间的 随机数，则仍接受状态 j 为当前状态； 若不成立，则保留状态 I 为当前状态。

算法其他参数的说明

退火过程由一组初始参数，即冷却进度表控制，它的核心是尽量使系统达到准平衡，以使算法在有限的时间内逼近最优解。冷却进度表包括：

1.控制参数的初值T_0:冷却开始的温度；
2.控制参数T的衰减函数：要将连续的降温过程，离散成降温过程中的一系列温度点，衰减函数即计算这一系列温度的表达式；
3. 控制参数T的终值T_f（停止准则）；
4. Markov链的长度L_k：任意温度T的迭代次数。

算法基本步骤

1、令T=T_0，随机生成一个初始解x_0，并计算相应的
      目标函数值E(x_0)；
2、令T等于冷却进度表中的下一个值T_i ；
3、根据当前解x_i进行扰动，产生一个新解x_j ，计相应
      的目标函数值E(x_j) ，得到 ∆e= E(x_j)−E(x_i);
4、 ∆e<0，则新解x_j被接受，作为新的当前解；
       ∆e>0，则新解x_j按概率P接受；
5、在温度T_i下，重复L_k次的扰动和接受过程，即执行
      步骤（3）、（4）；
6、判断T是否已达到T_f，是，则终止算法；否则转到
   （2）继续执行 

几点说明 

1、新解的产生    

  要求尽可能地遍及解空间的各个区域，这样，在某一 恒定温度下，不断产生新解时，就可能跳出局部最优解。

2、收敛的一般条件：

• 初始温度足够高；
• 热平衡时间足够长；
• 终止温度足够低；
• 降温过程足够缓慢；  

 举例

算法的MATLAB实现旅行商问题

模型：一名商人要到n 个不同的城市去推销商品，每2 个城市I 和j 之间的距离为d，如何选择一条路径使得商人每个城市走一遍后回到起点所走的路径最短。

例：有52座城市，已知每座城市的坐标，求每 个城市走一遍后回到起点，所走的路径最短。

 一、算法设计步骤

2.新解的产生（扰动）

（1）二变换法：任选序号u,v (设u<v<n)，交换u,v之间 的访问顺序。

（2）三变换法：任选序号u,v,w (设u≤v<w)， 将u,v之间 的路径插到w之后访问 

while t>=tf
 for r=1:Markov_length  
        if (rand < 0.5) 
     %随机产生0~1的数，若小于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;
            
 if (ind1 < ind2) && (ind2 < ind3)
            elseif (ind1 < ind3) && (ind3 < ind2)
                ind2 = tmp3;ind3 = tmp2;
            elseif (ind2 < ind1) && (ind1 < ind3)
                ind1 = tmp2;ind2 = tmp1;
            elseif (ind2 < ind3) && (ind3 < ind1)
                ind1 = tmp2;ind2 = tmp3; ind3 = tmp1;
            elseif (ind3 < ind1) && (ind1 < ind2)
                ind1 = tmp3;ind2 = tmp1; ind3 = tmp2;
            elseif (ind3 < ind2) && (ind2 < ind1)
                ind1 = tmp3;ind2 = tmp2; ind3 = tmp1;
            end      % ind1 < ind2 < ind3

            tmplist1 = sol_new((ind1+1):(ind2-1));  %u、v之间的城市
            sol_new((ind1+1):(ind1+ind3-ind2+1)) = ...
                sol_new((ind2):(ind3));   %将v到w的城市移到u后面
            sol_new((ind1+ind3-ind2+2):ind3) = ...
                tmplist1;     %u、v之间的城市移到w后面
        end

3.目标函数

访问所有城市的路径总长度：

模拟退火算法求出目标函数的最小值 

 % 计算目标函数即内能
        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
                % 冷却过程中最好的解保存下来´
                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