现代优化方法第3章模拟退火算法

1、第3章 模拟退火算法(Simulated Annealing) Metropolis在1953年提出 Kirkpatrick在1983年应用在组合优化问题31模拟退火算法及模型1.退火是一种物理过程一种金属物体在加热至一定的温度z在状态空间D中后，它的所有自由运动z随着温度的下降，这些留在不同的状态逐渐停z在温度最低时， 结构排列重新以一定的2 由统计力学的研究表明，在温度T，停留在状态r满足波尔兹曼(Boltzmann)概率分布1- E(r)=PrE(E)rexp()Z(T)kT (31)BE(r)为状态r的能量，kB0为波尔兹曼常量，Z(T)数，为能量的一个随为概率分布的标准化因子- E(

2、 s )T)=sZ(exp()kT DBT, Z(T) |D|3 研究由(31)确定的函数随T变化的趋势z 选定两个能量E1E2，在同一个温度T，有Pr(=E1)Pr(=E2)0, T0(32)z 在同一个温度，(3.2)表示停留在能量小状态的概率比停留在能量大状态的概率要大z 当温度相当高时，(31)的概率分布使得每个状态的概率基本相同，接近平均值1|D|，|D|为状态空间D中状态的个数z 结合(32)，具有最低能量状态的波兹曼概率接近并超出平均值1|D|z Pr=E(rmin)关于温度T是单调下降的：当rmin为D中具有最低能量的状态时，有Pr=E(rmin)/ T0z 当T趋向于0时,P

3、r=E(rmin)1/|DO|, T0当温度趋向0时，(31)决定的概率渐近l|DO|，D0为具有最低能量的状态集合。由此可以得到，在温度趋向0时，留在最低能量状态的概率趋向1停z的概率变化趋势图 在能量最低状态 非能量最低状态温度最低时，只有能量最低的点的概率不为零4可以将组合优化问题同金属物体退火进行类比组合优化问题解 最优解费用函数金属物体状态能量最低的状态能量由以上的类比及(31)，组合优化的最优解可以类比为退火过程中能量的最低状态，也就是温度达到最低点时，(31)概率分布中具有最大概率的状态51简单的模拟退火算法为：x0任选一初始解；xi :=x0;k:=0;t0 :=tmax(初始

4、温度)2若在该温度达到内循环停止条件，转3；否则，从邻域N(xi)中随机选一xj，计算fij=f(xj)-f(xi);若fij0，则xi:=xj;否则，若exp(-fij/tk) random(0,1)时，则xi :=xj;重复23tk+1:=d(tk);k:=k+1;若满足停止条件， 终止计算；否则，回2接受较差解的概率：E=fj-fi0Pj / Pi = exp(-fij /tk) random(0,1)6 模拟退火的直观理解是：在一个给定的温度，搜索从一个状态随机地变化到另一个状态每一个状态到达的次数服从一个概率分布(从邻域中选一点)当温度很低时，以概率1停留在最优解T1T2T3 D0模

5、拟退火算法的每一次迭代都体现集中和扩散两个策略的平衡对遇到的下一个迭代解，如果这个解更好，则采用集中策略选择这个解为新解不然的话，采用扩散策略，以某一概率选择这个解为新解7 模拟退火算法的数学模型z在给定邻域结构后，模拟退火过程是从一个状态到另一个状态不断地随机游动z可用马尔可夫(Markov)链描述这一过程转移概率(transition probability)z当温度t为一确定值时，两个状态的转移概率定义为：pijj iGij (t ) Aij (t ),p (t ) = |D|1 - Gil (t ) Ail (t ),j = iijl =1,l i(3.1.5)产生概率(generat

6、ion probability)l Gij(t)称为从i到j的产生概率表示在状态i时，j状态被选取的概率如j为i的邻居，等概率选取，则j被选中的概率:(t ) = 1 / | N (i) |, j N (i)G0,ijj N (i)(3.1.6)接受概率(acceptance probability)l Aij(t)称为接受概率, 表示在状态i产生j 后，接受j的概率模拟退火算法中的接受概率：f (i) 1,f ( j)Aij (t ) = exp(-Df/ t ), f (i) 0；状态i和状态j相通:ijz从i到达j的首达时刻随量为：Tij =minn|X(0)=i,X(n)=j,n1,(

7、n)其概率 fij=Pr(Tij=n|X(0)=i)z迟早到达概率 fij = n=1 fij(n)z定理3.2.2fij0的充分必要条件是ijz当fii1时，状态i称为常返的，当fii1时称为非常返的z定理3.2.3 状态j是常返的，则以概率1，系统无穷次返回状态j状态j是非常返的，则以概率1，系统只有有限次返回状态ju = nfz定义( n)ii(3.2.7)in=1平均返回步数z正常返:当ui时为正常返，当ui 时为零常返z定理3.2.3表明常返是以概率1无穷次返回同一状态，(3.2.7)表明有些常返状态(正常返)的平均返回次数(步数)是有限的，而有些是无穷的z一个状态i称为周期(per

8、iodic) t (t正整数)的，如果pij(n)除nt，2t，3t，外均为0，且t是满足这个条件的最大整 数不存在周期的状态称为非周期状态z非周期的正常返状态称为遍历状态z一个集合C是闭集：iC，jC，有pij=0。即i不可达jz除整个状态空间外，没有别的闭集的马氏链称为不可约(irreducible)的马氏链闭集1闭集2z定理3.2.4不可约的有限状态时齐的马氏链都是正常返的z平稳分布(stationary distribution):一个概率分布vj | j1，2，即vj0，jl，2，vj=1，如果对马氏链的一步转移概率阵P=(pij)满足，v j= vipij(3.2.8)i =1则称

9、vj | j1，2，)为马氏链的平稳分布z定理3.2.5非周期不可约时齐马氏链是正常返的充分必要条件是存在唯一平稳分布vj | j1，2，满足，= v p=( n)vv pjiijiiji=1i=11= lim p(n)= 0vjijunj其中平均返回步数uj由(3.2.7)给出z定理3.2.6有限状态时齐的马氏链是不可约的一个充分条件为：任给两个状态i和j，存在n，使得 pij(n)0z定理3.2.7 若i和j相通，即ij，则它们或同为非周期的或同为周期的33 时齐算法的收敛性1 讨论(3.1.5)中A(t)和G(t)所满足的条件，使得定理3.2.5的平稳分布可以地表达2 定理3.3.3如式

10、(3.1.5)中A(t)和G(t)满足：(1) G(t)与t无关；(2) i,jD，都有GijGji且存在q1，l0，l1，lqD, l0i，lqj，使得Glk lk +1 ( t ) 0, k = 0,1,L , q - 1(3)i,j,kD，f(i)f(j)f(k),都有Aik(t)Aij (t)Ajk(t)(4)i,jD，f(i)f(j),都有Aij(t)1(5)i,jD，t0,f(i)f(j),都有0Aij(t)1;则模拟退火时齐算法的马氏链有平稳分布v=(v1，v2，v|D|)，满足Ai i (t )v (t ) =, i D(3.3.1)0 Ai(t )i0 jjD其中i0DOPT

11、最优状态集合3很容易验证:l (3.1.7)满足定理3.3.3的(3)(4)(5)l (2)要求任何两个状态或互为邻居或互不为邻居，当互为邻居时，它们的产生概率相同当任何两个状态或互为邻居或互不为邻居时，如下面的产生概率：满足（1）和（2）的Gij=Gji ,L为正数l （2）的其他条件需要邻域是连通的4定理3.3.4 在满足定理3.3.3的条件下，当满足i,jD,f(i)f(j)limt0Aij(t)=0时，则有5 定理3.3.3和定理3.3.4给出了时齐算法的全局收敛应满足的充分条件，但这些条件的限制是比较强的，如(3.1.6)的发生概率就不一定满足定理3.3.3的(2)6定理3.3.6

12、若Aij(t) ，Gij(t)满足定理3.3.3除（2）外的条件，且条件（2）改为任何两个状态i和j或互为邻居或互不为，且Gij满足：状态空间D对邻域连通，则平稳分布是7例3.3.1四个状态1，2，3，4的邻居分别为： N(1)=1,2,N(2)=1,2,3,4,N(3)=2,3,4, N(4)=2,3,4易证它们互为邻居或互不为邻居，但G12=1/2G21=1/4,说明定理3.3.6成立的条件比定理3.3.3要广泛34非时齐算法收敛性简介z 非时齐算法一步转移概率：pij(k-1,k)=Pr(X(k)=j|X(k-1)=i)j iGij (tk ) Aij (tk ),= |D|1 -Gil

13、 (tk ) Ail (tk ),j = il =1,l ilim tk = 0其中tk-1 tkkz多步转移概率：pij(m,k):m步在i，第k步在jz定义3.4.2若存在向量v=(v1,v2, ,v|D|)，满足且有则称马氏链为强遍历（strongly ergodic）(渐进稳定性，收敛于一个随机分布)z定理3.4.3在温度t时，一步转移概率（3.1.5）中的A(t) 按（3.1.7）定义，即G(t)为当存在k02，使得kk0，其中Dfmax=maxf(i),iD -minf(i),iD则马氏链X(tk)k=1,2为强遍历z本定理对应的非时齐马氏链收敛到全局最优解，即limkP(X(K)

14、DOPT)=1温度下降速度最快为3.5实现的技术问题z 从应用的角度讨论算法实现中的一些技术问题z 主要包括解的形式、邻域的构成、初始温度的选取、温度下降的规则、每一温度下马氏链的迭代步长和停止规则等z按理论要求达到平稳分布来应用模拟退火算法是不可能的，这是因为时齐的算法要在每一个温度迭代无穷步以达到平稳分布，而非时齐要求温度下降的迭代步数是指数次的从应用的角度来看，在可接受的时间里得到满意的解就可以了，因此本节介绍的技术问题无法保证模拟退火算法得到全局最优解z应用这些技术的模拟退火算法还是一种启发式算法351 解的形式和邻城结构解的表现形式直接决定于邻域的构造z第1个问题：f(x)=x2 ,

15、0x127,x 为整数。解可用:邻域:s=（1001011）7位进制表示N(S)=S|S是一个0-1码且|S-S|7= | s- s | k, k 1整数iii=1z第2个问题:TSP可用n个城市的一个排序表示问题的一个解。邻域采用2-optz以上两个问题的共同特征是：1) 每一个邻域所含的状态个数相同2) 每个邻居都是可行解3) 空间中的任何两个状态可达于是由时齐算法的理论，只要接受概率满足一定的条件，一定以概率1收敛到全局最优解z同样解的表达式和邻域结构，将第1个问题的方法应用到背包问题时，第2 个问题的方法应用到车间作业问题(第2 章242节)由于背包问题有包容积约束的限制和车间作业的死

16、锁现象，无法保证每个邻域中邻居都是可行解。因此，从理论上无法保证收敛到全局最优解。z处理这个问题有两类方法：用罚函数将不可行解可行化，这使得问题转化为上面讨论的两个问题，原有的方法可以继续使用采用罚函数的方法对不可行解进行处理时，一个主要缺陷是引起搜索区域的扩大，从而使计算时间增加为了克服这个不足，另一个方法是研究解和邻域结构，从理论上保证模拟退火算法以概率1收敛到全局最优解这一个方法要求对问题有相当深入的了解，因此，可以说其难度是比较大的z例3.5.1背包问题：3个物品，体积5，3，1，包容积为4，则(0,0,0)、(0,0,1)、(0,1,0)、(0,1,1)为可行解集合，采用罚值法，其余

17、4个解也为搜索对象。z罚值的方法可以很简单地保证理论要求的不可约和非周期性但如此计算一定会增加计算时间，而有时增加时间之多是不可以接受的如按例1.1.2的TSP的表示方法0,1n(n-1)，当n=10时，扩大31021倍z邻域的构造同计算时间紧密在一起，也就产生了对实际问题的理论分析同实际应用的352温度参数的控制温度参数是模拟退火时齐算法中一个最关键的参数之一主要包括起始温度的选取、温度的下降方法和停止温度的确定等1.起始温度的选取z从理论上来说，起始温度t0应保证平稳分布中每一状态的概率相等，也就是使Aij(t0)1z很容易估计一个值为的数，其中，t0K，K充分大maxf(j)| jDmi

18、nf(j)| jD实际计算中，可选K10，20，100等试验值z对一些问题，可简单估计,如TSP:maxf(j)| jD maxdij| ji,j=1,2,n minf(j)| jD mindij| ji,j=1,2,nz但有的时候，会出现比较难估计或很难精确地知道最大值和最小值而使上面的估计太粗z另外，K的选取过大会造成计算时间的增加，K过小则使算法过早陷入局部最优点z由此产生数值计算估计和统计推断估计方法z数值计算估计方法的基本思想是给出一个值0(0接近1，如0.9,0.8等)，对给定的初始温度t0，用以下的算法：z通常采用统计的方法估计费用函数的上下限（）统计的方法总是在一定的概率分布下

19、研究问题的均值、方差等如对大样本正态分布：样本均值统计量样本方差统计量目标值以9997%的概率落入该区间2时齐算法的温度下降方法1)tk+1=tk,k0, 012)tk=t0(K-k)/K，K为算法温度下降的总次数3)以理论中的平稳分布为依据，定义平稳分布时的期望、二阶矩(3.5.4)- (3.5.7)n个样本的均值、二阶矩方差统计量2（tk）:样本均值统计量（tk）：温度的变化为(3.5.12)4)其他3.时齐算法每一温度的迭代长度规则(1)固定长度：在每一温度，选代相同的步数，步数的选取同问题的大小有关，通常采用与邻域大小直接相关的规则如TSP:2-opt邻域有Cn2 =n(n-1)/2个

20、邻居,在同一温度,迭代步长可取n，n/2，n2，n2/2，nk等(2)由接受和拒绝的比率来控制迭代步数：z当温度很高时，每一个状态被接受的频率基本相同，而且几乎所有的状态都被接受此时，在同一温度应使迭代的步数尽量小z当温度渐渐变低时，越来越多的状态被拒绝如果在此温度的迭代太少，则可能造成过早地陷入局部最优状态z比较直观和有效的方法是随着温度的下降，将同一温度的选代步长增加实现的法是给定一个充分大的步长上限U和一个接受次数指标R，当接受次数等于R时，在此温度不再迭代而使温度下降，否则，一直迭代到上限步数实现的第二种方法是给定一个接受比率指标R、迭代步长上限U和下限L，每一温度至少迭代L步且同一温

21、度迭代的总次数和被接受的次数，当迭代超过L步时，若接受次数同总次数的比率不小于R时，在这一温度不再迭代而开始温度下降，否则，一直迭代到上限步数z同样可以用拒绝次数为指标(3)概率控：以概率1跳出任何一个局部最优解是概率控想的一个基本思实际应用中，(3)不适用 (1)和(2)比较容易实现4.算法的终止原则(1)零度法： 模拟退火的最终温度为零因而最为简单的原则是：给定一个比较小的正数，当温度tk时，算法停止(2)循环总数控：这一原则可分为两类z整个算法的总迭代步数为一固定数，它表示各温度时马氏链迭代数(内循环)的总和为一个给定的数这样很容易计算算法的复杂性，但需要合理分配内循环的长度和温度下降的

22、次数z内循环的次数由迭代长度规则的(2)和(3)决定，温度下降次数（外循环）为一个定值。这样的控对估计算法的复杂性有一定的解决的方法一是通过理论的估计，二则是类似迭代长度规则的(2)给出每一温度的迭代长度上限(3)基于不改进规则的控：在一个温度和给定的迭代次数内没有改进当前的局部最优解，则停止运算。模拟退火的一个基本思想是跳出局部最优解直观的结论是在较高的温度没能跳出局部最优解，则在低的温度跳出最优解的可能也比较小由此产生上面的停止原则(4)接受概率控： 该方法与(3)相同的思想给定一个指标f0是一个比较小的数，除当前局部最优解以外，其他状态的接受概率都小于f时，停止运算（5）邻域法：如采用定

23、理3.3.6的发生概率和(3.1.7)的接受概率，但当满足下式时次优邻域内的局部最优邻域内的局部最优到次优的接受概率除局部最优解以外状态的接受概率都小于邻域总点平均数，此时可以认为从局部最优解转移到其他状态的可能性很小，因此停止得终止温度(6) Lundy和Mees方法：从概率的意义给出一个判定方法(7) Aarts和Van Laarhoven法36应用案例下料问题（Cutting stock problem)z一般的套裁方式第二次下料时，重合部分归选上的那块或归面积较大的那块1第二次下料2第一次下料1K次下料后，有K1个不相交的矩形(边界可以相连)和两个有相同公共部分的矩形，如此共计K十1个

24、矩形z直观可知，套裁节约用料平面上的套裁问题是二维装箱问题，是一维装箱问题的推广，因此，这个优化问题是NPhardz如果将的下料顺序看成问题的一个解，现在必须给出如何从诸多个矩形中选择一个，及从何处下料当这些问题解决后， 就可以计算余料的大小，即目标值z优化问题就是确定一个最优的下料序用模拟退火求解时，其邻域的结构可以模仿TSP中的2opt来构造，交换两个的加工顺序于是问题解的形式和邻域结构已经得到z如何选择矩形原料下料的思想：在给定的一个下料序后，不求解这个下料序的最优下料方法(若求解，相当于对原问题求最优下料方法)，用一种启发式方法安排下料，在完全排完后，就可以计算余料的值，计为目标值z选择矩形原料的启发式规则是根据等待下料的，采用下列方案，(1)首选，的长或宽同原料的长或宽不超过一个给定的比较小的正数(2)其次,比较原料的长宽之和同长宽之和的差值的大小，从中选差值