模拟退火算法

1、模拟退火算法汇报人： 许炯楼2014xxxxxx 李娜2014200909冒亚婷2014200922李园园20142009231 1 31.1 模拟退火算法的来源及基本原理w算法的提出算法的提出 模拟退火算法最早的思想由模拟退火算法最早的思想由Metropolis等（等（1953）提出，）提出，1983年年Kirkpatrick等将其应用于组合优化。等将其应用于组合优化。w算法的目的算法的目的 解决解决NP复杂性复杂性问题；问题； 克服优化过程陷入局部极小；克服优化过程陷入局部极小； 克服初值依赖性。克服初值依赖性。 41.1 模拟退火算法的来源及基本原理n物理退火过程：物理退火过程：退火是指

2、将固体加热到足够高的温度，使分子呈随机排列状态，退火是指将固体加热到足够高的温度，使分子呈随机排列状态，然后逐步降温使之冷却，最后分子以低能状态排列，固体达到某然后逐步降温使之冷却，最后分子以低能状态排列，固体达到某种稳定状态。种稳定状态。 加温过程加温过程增强粒子的热运动，消除系统原先可能存在的非均匀态；等温过程等温过程对于与环境换热而温度不变的封闭系统，系统状态的自发变化总是 朝自由能减少的方向进行，当自由能达到最小时，系统达到平衡态；冷却过程冷却过程使粒子热运动减弱并渐趋有序，系统能量逐渐下降，从而得到低能的 晶体结构。 51.1 模拟退火算法的来源及基本原理模仿自然界退火現象而得，利用

3、了物理中固体物质的退火过程与一般优化问题的相似性； 从某一初始温度开始，伴随温度的不断下降，结合概率突跳特性在解空间中随机寻找全局最优解。 61.1 模拟退火算法的来源及基本原理w数学表述数学表述 在温度在温度T，分子停留在状态，分子停留在状态r满足满足Boltzmann概率分布概率分布DsBBBTksETZTZkrrEETkrETZrEEP)(exp)()(Boltzmann0)()(exp)(1)(子：为概率分布的标准化因常数。为的能量，表示状态机变量，表示分子能量的一个随 71.1 模拟退火算法的来源及基本思想w数学表述数学表述 在在同一个温度同一个温度T，选定两个能量，选定两个能量E1

4、01 81.1 模拟退火算法的来源及基本思想BoltzmanBoltzman概率分布概率分布告诉我们：告诉我们： （1）在同一个温度，分子停留在能量小状态的概率大于停留在能量大状态）在同一个温度，分子停留在能量小状态的概率大于停留在能量大状态的概率的概率 （2）温度越高，不同能量状态对应的概率相差越小；温度足够高时，各状）温度越高，不同能量状态对应的概率相差越小；温度足够高时，各状态对应概率基本相同。态对应概率基本相同。 （3）随着温度的下降，能量最低状态对应概率越来越大；温度趋于）随着温度的下降，能量最低状态对应概率越来越大；温度趋于0时，时，其状态趋于其状态趋于1 91.1 模拟退火算法的

5、来源及基本原理w数学表述数学表述 若若|D|为状态空间为状态空间D中状态的个数，中状态的个数，D0是具有最低能量的状态集合：是具有最低能量的状态集合： 当温度很高时，每个状态概率基本相同，接近平均值当温度很高时，每个状态概率基本相同，接近平均值1/|D|； 状态空间存在超过两个不同能量时，具有最低能量状态的概率超出平均值状态空间存在超过两个不同能量时，具有最低能量状态的概率超出平均值1/|D| ； 当温度趋于当温度趋于0时，分子停留在最低能量状态的概率趋于时，分子停留在最低能量状态的概率趋于1。能量最低状态能量最低状态 非能量最低状态非能量最低状态 101.1 模拟退火算法的来源及基本原理wM

6、etropolis准则（准则（1953）以概率接受新状态以概率接受新状态 固体在恒定温度下达到热平衡的过程可以用固体在恒定温度下达到热平衡的过程可以用Monte Carlo方法（计算机方法（计算机随机模拟方法）加以模拟，虽然该方法简单，但必须大量采样才能得到比较精随机模拟方法）加以模拟，虽然该方法简单，但必须大量采样才能得到比较精确的结果，计算量很大。确的结果，计算量很大。 若在温度若在温度T，当前状态，当前状态i 新状态新状态j; 若若Ej=randrom0,1 s=sj; Until 抽样稳定准则满足；抽样稳定准则满足； 退温退温tk+1=update(tk)并令并令k=k+1； Unti

7、l 算法终止准则满足；算法终止准则满足； 输出算法搜索结果。输出算法搜索结果。u基本步骤 121.1模拟退火算法的基本思想和步骤模拟退火算法的基本思想和步骤 给定初温给定初温t=t0，随机产生初始状态，随机产生初始状态s=s0，令，令k=0； Repeat Repeat 产生新状态产生新状态sj=Genete(s)； if min1,exp-(C(sj)-C(s)/tk=randrom0,1 s=sj; Until 抽样稳定准则满足；抽样稳定准则满足； 退温退温tk+1=update(tk)并令并令k=k+1； Until 算法终止准则满足；算法终止准则满足； 输出算法搜索结果。输出算法搜索结

8、果。u影响优化结果的主要因素三函数两准则三函数两准则初始温度初始温度 131.1模拟退火算法的基本思想和步骤模拟退火算法的基本思想和步骤u模拟退火算法的步骤Step1 设定初始温度设定初始温度t = tmax, 任选初始解任选初始解r = r0Step2 内循环内循环 Step2.1 从从r的邻域中随机选一个解的邻域中随机选一个解rt, 计算计算r和和rt对应目标函对应目标函 数值数值, 如如rt对对 应目标函数值较小，则令应目标函数值较小，则令r = rt; 否则若否则若 exp(E(rt)E(r)/t)random(0,1), 则令则令r=rt. Step2.2 不满足内循环停止条件时，重

9、复不满足内循环停止条件时，重复Step2.1Step3 外循环外循环 Step3.1 降温降温t = decrease(t) Step3.2 如不满足外循环停止条件，则转如不满足外循环停止条件，则转Step2；否则算法结束；否则算法结束1. 达到终止温度2. 达到迭代次数3. 最优值连续若干步 保持不变1. 目标函数均值稳定2. 连续若干步的目标值变化较小3. 固定的抽样步数模拟退火算法的马氏链描述模拟退火算法的马氏链描述2 2 152.1马尔科夫链u定义) 1()(Pr) 1(,) 1 (,) 0()(Pr )( )(,1021inXjnXinXiXiXjnXZnkXkkXss，满足称为马尔

10、可夫链，若随机序列时刻状态变量的取值。为间，为所有状态构成的解空令一步转移概率：一步转移概率：n步转移概率：步转移概率：) 1()(Pr) 1(,inXjnXnpji若解空间有限，称马尔可夫链为若解空间有限，称马尔可夫链为有限状态有限状态；)0()(Pr)(,iXjnXpnji若若 ，称马尔可夫链为，称马尔可夫链为时齐的时齐的。) 1()(,npnpZnjiji 162.2 模拟退火算法与马尔科夫链模拟退火算法与马尔科夫链w模拟退火算法对应了一个马尔可夫链模拟退火算法对应了一个马尔可夫链 模拟退火算法：新状态接受概率仅依赖于新状态和当前状态，并由温度加以控制。 若固定每一温度，算法均计算马氏链

11、的变化直至平稳分布，然后下降温度，则称为时齐算法； 若无需各温度下算法均达到平稳分布，但温度需按一定速率下降，则称为非时齐算法。w分析收敛性分析收敛性模拟退火算法关键参数和操作的设计3 3 * *3 模拟退火算法关键参数和操作的设计原则原则 产生的候选解应遍布全部解空间方法方法 在当前状态的邻域结构内以一定概率方式（均匀分布、正态分布、指数分布等）产生 状态产生函数状态产生函数 状态接受函数状态接受函数原则原则 (1)在固定温度下，接受使目标函数下降的候选解的概率要大于使目标函数上升的候选解概率； (2)随温度的下降，接受使目标函数上升的解的概率要逐渐减小； (3)当温度趋于零时，只能接受目标

12、函数下降的解。方法方法 具体形式对算法影响不大 一般采用一般采用min1,exp(-C/t) * *3 模拟退火算法关键参数和操作的设计收敛性分析收敛性分析 通过理论分析可以得到初温的解析式，但解决实际问题时难以得到精确的参数； 初温应充分大；实验表明实验表明 初温越大，获得高质量解的机率越大，但花费较多的计算时间；方法方法 （1）均匀抽样一组状态，以各状态目标值得方差为初温； （2）随机产生一组状态，确定两两状态间的最大目标值差，根据差值，利用一定的函数确定初温； （3）利用经验公式。 初温初温 * *3 模拟退火算法关键参数和操作的设计时齐算法的温度下降函数时齐算法的温度下降函数 （1）

13、，越接近1温度下降越慢，且其大小可以不断变化； （2） ，其中t0为起始温度，K为算法温度下降的总次数。 温度更新函数温度更新函数10 , 0 ,1kttkk0tKkKtk * *3 模拟退火算法关键参数和操作的设计非时齐模拟退火算法非时齐模拟退火算法 每个温度下只产生一个或少量候选解时齐算法时齐算法常用的常用的Metropolis抽样稳定准则抽样稳定准则 （1）检验目标函数的均值是否稳定； （2）连续若干步的目标值变化较小； （3）按一定的步数抽样。 内循环终止准则内循环终止准则 外循环终止准则外循环终止准则常用方法常用方法 （1）设置终止温度的阈值； （2）设置外循环迭代次数； （3）算法

14、搜索到的最优值连续若干步保持不变； （4）概率分析方法。 实例计算4 4成都市三环线-绕城高速西北区域公交路线的设计 * *4.1 问题描述发展成都市经济文化的发展导致三环路以外居民不断普及，绕城高速区聚集了越来越多的企业、学校和商业网点。123经济的发展凸显了公共交通网络的不完善，成都市现有的公交线路覆盖范围主要涉及主城区，但郊区的公交线路和数量都很少，无法满足人们的出行要求。鉴于此，对未开通公交车的地区进行线路的优化是十分有必要的。优化目标：明确选择成都三环线绕城高速的西北区域范围（必须包括大丰镇、安靖镇、犀浦镇和高新西区等重点区域），利用数学模型，在考虑人口密度、交通路况等基础上，对该区

15、域的公交路线进行合理的规划设计，以满足居民在该区域内或区域外工作和生活的需要； * *4.2 调查分析 在成都三环线绕城高速的西北区域范围（包括大丰镇、安靖镇、犀浦镇和高新西区等重点区域）内： 一共收集到了通过该区域的公交车辆3838辆（线路7676条，往返线路可以不一致），在该区域范围内的公交站点，以及每条公交线路的总长、在区域范围内的实际长度和空间直线长度，同时也刻画了公交站点在区域范围内的相对位置。注：在进行线路优化时，总的公交线路条数是不变的（7676条），不添 加额外的公交线路。 时间代价矩阵 * *4.3 模型建立（1）数据准备客流矩阵OD距离矩阵 公交网络矩阵21路况矩阵345（

16、2）公交空间网络解空间M的建立记原来区域内的路线组成的交通网络为S0，根据S0的每一条线路在区域上的起讫点，运用 前n条最短路算法，确定S0的每一条路线的n条备用路线，建立区域内公交网络解空间M。记原来区域内的路线组成的交通网络为S0，显然 S0 属于M，所设计的新公交路线便从M中选取。（生成代码见附录2）Dijkstra其中： * *4.3 模型建立（3）目标函数的确定对于每一个交通网络 根据相应的数据库得出的数据所确定时间成本函数 ，我们的目标即在解空间内找到最优解S使得 最小，即求： 表示从站点i到j站点的第m条直达公交线路上的人流量； 表示从站点i到j站点的第m条需换乘公交线路上的人流

17、量； 表示第m条直达公交线路上从站点i到站点j所需行车时间； 表示第m条需换乘公交线路上从站点i到站点j所需行车时间。( )mmmmmijmijrrtrtrijijijiji N j N rDRi N j N trTRMinT Sdtdt mrijdmtrijdmrijtmtrijt,1,1( ,1)( , )mmmmrrrijk kk kk kri jtDv ,1,1( ,1)( , )r( , )mmlmtrtrrijk kk klmk ktri jtDvtr i j +5K （v ：第m条公交线路上相邻站点k与k+1之间的平均速度；v ：第m条公交线路上相邻站点k与k+1之间的距离；,1

18、mrk kv,1mrk kD * *4.3 模型建立v ：第m条需换乘公交线路上相邻站点k与k+1之间的距离；v ：第m条需换乘公交线路上换乘的次数。,1mtrk kDK转化为解空间的表达z( )( ,)( )( ,)( )mmmmmijmijrrtrtrijijijiji N j N rDRi N j N trTRMinSdS odtSdS odtS 客户评价内容S：解空间产生的一个公交路线网络M：公交路线网络解空间 Pih：第i条路线的原起始点 li：解空间中第i条线路的一条备选路线 Pit：第i条路线的原始终点Di ：第i条路的原始距离Cover1：交通网络的1级覆盖率Cover2：交通

19、网络的2级覆盖率. .1( )0.92( )0.5SMstCover SCoverS511.51.4ihiitiiiiiPlPlll D . .st * *4 利用算法求解第一步 确定目标函数第二步 利用模拟退火算法进行最优解计算根据路况矩阵，以及odod矩阵和相关参数确定目标函数T(S)T(S) 设置初始温度T，和截止温度T （min），以及截止指数k（初始为0）退火指数r20时输出S，否则返回（2）（4）)() (STSTt0t SSS TkzbetP )( S)( tPRSS 第三步 列出所得最优解，进行合理性测试。 * *4.4 算法程序function S ,L3new = tuih

20、uoNew(Sc,OD, condition , L,L3,L4 ) % L3 原始解 L4 L4 解空间S=Sc;S=Sc;S2=Sc;S2=Sc;k=0;T=T=100;00;k_ _max=21T_min=20;T_min=20;r=0.95;r=0.95;while Twhile TT T_minmin & k0if de0 S=S2; S=S2; k= =0;elseelse if ( exp( de/T ) rand ( 1 ) ) if ( exp( de/T ) rand ( 1 ) ) S=S2; S=S2; k= =0; end endendendT=r*T;k=k+1;e

21、ndendL3new=L3;L3new=L3;endend模拟退火算法function A = buildA( S )function A = buildA( S )%UNTITLED2 Summary of this function goes here%UNTITLED2 Summary of this function goes here% Detailed explanation goes here% Detailed explanation goes here%到了第三问这里有buildA for count version two.buildA for count version

22、two.A=zeros(244,244);A=zeros(244,244);for i=1:244for i=1:244for j=1:244for j=1:244if length(Si,j)=1if length(Si,j)=1 temp=(Si,j); temp=(Si,j); aaa=mean(temp,1); aaa=mean(temp,1); A(i,j)=aaa(2); A(i,j)=aaa(2);endendendendendendfor i=1:244for i=1:244 for j=1:244 for j=1:244 if A(i,j)=0&i=j if A(i,j)=0

23、&i=j A(i,j)=inf; A(i,j)=inf; end end end endendend代价矩阵A代码 * *4.4 算法程序function pathaim,dist=dijkstra2(D,s,aim)function pathaim,dist=dijkstra2(D,s,aim)%Dijkstra最短路算法MatlabMatlab程序用于求从起始点s s到其它各点的最短路%D为赋权邻接矩阵%d为s s到其它各点最短路径的长度%DD记载了最短路径生成树tic;tic;m,n=size(D);m,n=size(D);d=inf.d=inf.* *ones(1,m);ones(1,

24、m);d(1,s)=0;d(1,s)=0;dd=zeros(1,m);dd=zeros(1,m);dd(1,s)=1;dd(1,s)=1;y=s;y=s;for iii=1:mfor iii=1:m pathiii,1(1)=iii; pathiii,1(1)=iii;endendwhile length(find(dd=1)mwhile length(find(dd=1)m for i=1:m for i=1:m if dd(i)=0 if dd(i)=0 d(i)=min(d(i),d(y)+D(y,i); d(i)=min(d(i),d(y)+D(y,i); if d(i)=d(y)+D(y,i) if d(i)=d(y)+D(y,i) pathi=pathy;pathi(length(pathy,1)+1)=i; pathi=pathy;pathi(length(pathy,1)+1)=i; end end end end end end ddd=inf; d