【老生谈算法】matlab实现模拟退火算法——模拟退火算法

模拟退火算法 matlab程序

---

1、文档下载：

本算法已整理成文档如下。有需要的朋友可以点击下载

序号	文档(点击下载)
本项目文档	模拟退火算法-matlab2.doc

---

2.算法详解:

模拟退火算法来源于固体退火原理，将固体加热到足够高，然后慢慢冷却。加热时，固体内的颗粒随着温度的升高而变得无序，内部可以增加。然而，当颗粒慢慢冷却时，颗粒逐渐变得有序，在每个温度下达到平衡状态，最终在室温下达到基态，最大限度地减少内部能量。

根据Metropolis温度T时粒子平衡的概率为e-ΔE/(kT)，E是温度T时的内能，ΔE为其变量，k为Boltzmann常数。将内能E模拟为目标函数值f，温度T演变成控制参数t，即解决组合优化问题的模拟退火算法:从初始解i和控制参数初始值t开始，对当前解重复产生新的解决方案→计算目标函数差→接受或放弃的迭代，并逐渐衰减t值，算法终止时的当前解是近似最优解，这是一个启发性的随机搜索过程，基于蒙特卡罗迭代求解。

冷却进度表退火过程(Cooling Schedule)控制，包括控制参数的初始值及其衰减因子Δt、每个t值的迭代次数L和停止条件S。

模拟退火算法可分为三空间、目标函数和初始解三部分。

模拟退火的基本思想： (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->然后转第二步。

新解模拟退火算法的产生和接受可分为以下四个步骤： 第一步是从当前的解决方案中产生一个新的解决方案；为了便于后续计算和接受，减少算法耗时，通常选择当前的新解决方案，如构成新解决方案的全部或部分元素的替换和交换，并注意到新解决方案的转换方法决定了当前却进度表的选择有一定的影响。 第二步是计算与新解决方案对应的目标函数差。由于目标函数差仅由变换部分产生，因此最好根据增量计算目标函数差。事实表明，这是大多数应用程序计算目标函数差的最快方法。 第三步是判断新解决方案是否被接受。判断的基础是接受标准。最常用的接受标准是Metropo1is准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。 第四步是当新解决方案被确定和接受时，用新解决方案代替当前解决方案，这只需要实现当前解决方案对应于产生新解决方案时的变换部分，并修改目标函数值。此时，当前的解决方案已经实现了迭代。下一轮测试可以在此基础上开始。当新解决方案被判定为放弃时，下一轮测试将继续基于当前解决方案。

模拟退火算法与初始值无关。S(算法迭代的起点)无关；模拟退火算法具有渐近收敛性，理论上已被证明是概率l 模拟退火算法是并行的。

模拟退火算法流程图：

一、glory.m文件

T_max=input('please input the start temprature'); T_min=input('please input the end temprature'); iter_max=input('please input the most interp steps on the fit temp'); s_max=input('please input the most steady steps ont the fit temp'); T=T_max; load d:\address.txt; order1=randperm(size(address,1))';%产生初始解。 plot(address(order1,1),address(order1,2),'*r-') totaldis1=distance(address,order1);
while T>=T_min
    iter_num=1;
    s_num=1;
      plot(T,totaldis1)
    hold on
    while iter_num<iter_max&s_num<s_max;
        order2=exhgpath(order1);
        totaldis2=distance(address,order2);
        R=rand;
        DeltaDis=totaldis2-totaldis1;
        if DeltaDis<0;
            order1=order2;
             totaldis1=totaldis2;
        elseif (exp((totaldis1-totaldis2)/(T))>R)
            order1=order2;
            totaldis1=totaldis2;
        else s_num=s_num+1;
        end
        iter_num=iter_num+1;
    end
    T=T*0.99;
end
order1
totaldis1
figure(2)
plot(address(order1,1),address(order1,2),'*r-')

二、distance.m文件

function y=distance(address,order)
nmb=size(address,1);
y=0;
for i=1:nmb-1
    y=y+sqrt((address(order(i+1),1)-address(order(i),1))^2+(address(order(i+1),2)-address(order(i),2))^2);
end
y=y+sqrt((address(order(i+1),1)-address(order(1),1))^2+(address(order(i+1),2)-address(order(1),2))^2);

三、exhgpath.m文件

function y=exhgpath(order)
   while 1
       b=size(order,1);
       r=unidrnd(b,1,2);
       if r(1)-r(2)~=0
           break
       end
   end
   b=order(r(2));
   order(r(2))=order(r(1));
   order(r(1))=b;
   y=order;