第14讲 智能计算模型 —模拟退火算法 今天我们来讲另外一种解决优化问题的算法-模拟退火算法。上一节课我们讲的遗传算法是通过仿照生物进化过程来设计算法，其关键是将生物进化过程中的自然选择-好的个体有更大可能遗传到下一代来保证算法使得目标函数朝着最优方向前进，而根据遗传中的交叉、变异来使得在向最优解前进过程中进行全局和局部的随机游走，从而避免陷入局部最优解，最终寻找全局最优解。遗传算法是一种群体行为。今天我们讲的模拟退火算法是仿照燃烧的物体退火过程中的行为来寻找最优解，其也是一种类比行为。

【主要内容】 介绍模拟退火算法的主要思想、算法流程以及在数学建模中的应用。 【主要目的】 掌握数学建模过程中模拟退火算法在求解优化问题中的应用。

固体退火原理: 将固体加温至充分高，再让其徐徐冷却。加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。

模拟退火算法基本原理： 模拟退火算法的基本原理是模拟固体退火过程中总是从能量高的状态向着能量最低的平衡态转换的思想来寻找最优解，其通过一个冷却温度表来控制这个过程，同时在每一个能量状态设计解的随机游走并以一定概率接收差的解，而且随着能量的降低，接收差的解的概率也显著降低，从而使得在高能量状态具有逃离局部最优解，在低能量状态具有收敛全局最优解的能力。 为什么能量高往能量低走可以转换到优化的思想中？ 相比遗传算法的思路，我们想想其最大的不同点在哪里？ 算法自然满足在温度高（差的解的时候）具有较大随机性，在温度低（好的解）的时候具有较小随机性，从而使得算法的收敛性天生就较好

模拟退火算法的一般描述： 模拟退火算法是一种模拟物理机理过程来确定优化方向的随机寻优算法。算法类比物理机理来设计最优性搜索方向，通过设计解的随机游动提高解的全局性，从而达到既具有全局搜索能力又具有收敛性保障，针对目标极值多、非线性程度高等一类优化问题有较好的效果。 关键问题： （1）退火温度控制； （2）随机搜索策略 这里模拟的物理机理，就是燃烧的物理退火过程

将固体退火过程与优化问题求解对应起来 高温 低温 T 能量高，粒子热运动活跃 目标函数大，解随机游动最大 固体退火 问题寻优 系统趋于平衡态，能量达到最小 解随机游动，趋向目标函数最小 这里有一个问题就是要注意，在固体退火中是有体现的，就是能量越低，粒子运动越弱，也就是说要设计使得在优化问题中温度越低，解随机游动最小 等温过程是算法的一个最关键的过程，如何实现？ 能量最低，粒子热运动最弱，系统稳定 目标函数最小，解随机游动最小

上述固体退火的第二个过程，即指定温度下系统达到平衡的过程可以用蒙特卡洛方法进行模拟。 (1)使用随机数发生器产生一个随机的分子构型； (2)对此分子构型的粒子坐标作无规则的改变，产生一个新的分子构型； (3)计算新的分子构型的能量； (4)比较新的分子构型与改变前的分子构型的能量变化，判断是否接受该构型。 物理上采用这样的仿真手段来模拟固体退火的第二个过程

若新的分子构型能量低于原分子构型的能量，则接受新的构型，使用这个构型重复再做下一次迭代。 若新的分子构型能量高于原分子构型的能量，则计算玻兹曼因子，并产生一个随机数。 若这个随机数大于所计算出的玻兹曼因子，则放弃这个构型，重新计算。 若这个随机数小于所计算出的玻兹曼因子，则接受这个构型，使用这个构型重复再做下一次迭代。 (5)如此进行迭代计算，直至最后搜索出低于所给能量条件的分子构型结束。 这里有个问题：你如何知道其一定会稳定？稳定就是分子不再动了，我们还需要从波兹曼概率分布来看， 我们把这个接收概率抽取出来，得到这样一个准则

Metropolis准则说明，若新状态的能量变小，则接收，若新状态能量变大，则以一定概率接收，即以一定概率接收差的状态。 假设固体中当前所有粒子的状态用状态i表示,该状态的能量为Ei，现给固体中粒子施加随机扰动得到一个新的状态j，该状态能量为Ej，则固体从状态i转变到状态j的接收概率为： 温度 波兹曼常数 我们下面再来分析一下接收概率的性质 Metropolis准则说明，若新状态的能量变小，则接收，若新状态能量变大，则以一定概率接收，即以一定概率接收差的状态。

若能量变化量相同，则温度越高，系统接收差的状态概率越大，反映出系统越活跃,反之，温度越低，系统越稳定。 能量改变量 温度 若能量变化量相同，则温度越高，系统接收差的状态概率越大，反映出系统越活跃,反之，温度越低，系统越稳定。 若温度相同，则能量改变量越大，系统接收差的状态概率越小。 接收概率与能量该变量以及温度有关 下面我们定性来看看接收概率与能量该变量以及温度的关系 下面我们来看看Aij曲线的变化 若温度无限大，则系统以近似1概率接收一切状态。若温度为0，则系统不接受任何差的状态。

这是从能量该变量大小来看温度不同是接收差的解的状态 能量改变不大时，不管是高温还是低温，都有可能接收差的状态（改变状态），其反映了在能量变化局部范围内高低温都可能改变，能量改变大（全局性）时低温不再接收状态变差。

温度低时，能量变化越小（局部）越可能状态改变，能量变化越大（全局)越不能改变，高温时，能量变化不同大小都有可能改变。 这是从温度变化角度来看能量变化 总结一下：在温度高时，系统状态变化概率大，但能量变化很大的概率还是小，温度低时，系统状态概率变化小，但小的能量变换（状态转化）还是有，但总的来说，系统是朝着最低能量变化 这就说明温度和能量，也对已对应到目标函数可以很好的配合来确定状态变换的方向，也就是目标寻优的方向 温度低时，能量变化越小（局部）越可能状态改变，能量变化越大（全局)越不能改变，高温时，能量变化不同大小都有可能改变。

将Metropolis准则与固体退火过程结合起来，就可以构造一种优化算法框架—模拟退火算法(Simulated Annealing Algirithm) 物理退火 解 粒子的状态 最优解 能量最低的粒子状态 目标或评价函数 能量 初温 退火初始状态 邻近解 转换状态 Metropolis采样过程 等温过程 控制参数下降 温度降低，冷却 注意这个地方是一个随机游动的构成，可能接收差的解 所以模拟退火算法实际上是有双重循环的，第一重是控制参数（也就是温度）在不断下降，第二重，是对于每一个控制参数，又存在一个随机游动

Metrapolis采样得到当前温度下目标最优解 模拟退火算法流程 选定初始温度和初始解 降温 Metrapolis采样得到当前温度下目标最优解 输出最优解 Y N 满足结束条件 或降到最低温 外循环： 降温过程 内循环： 等温过程 蒙特卡洛仿真

算法步骤 （1）初始化：初始温度T（充分大），初始解S，每个温度下的迭代次数为L （2）对k=1,2,...,L,执行以下过程 ① 产生新解S’; ② 计算目标增量ΔE=f(S’)-f(S)，其中f(S)为评价函数或目标函数 ③ 若ΔE<0，则以概率1接受S’作为当前新解，否则计算接收概率A

④ 产生一个随机数ξ，若ξ<A，则接受S’作为当前新解，否则放弃新解，当前解不变 ⑤ 如果循环终止或达到结束条件，则转下一步 （3）温度T按照某种规则减小，且T→0，以第（2）最终输出的解为当前状态重新转第（2）步执行。 注意到一点，模拟退火算法是单粒子行为，遗传算法是种群行为 下面我们从理论上来看看这个算法

理论描述 理论上可以用马尔科夫链来描述上述过程。在保温过程，系统是从一个状态通过随机游动转换到另一个状态，那么t温度下，两个状态之间的转移概率可以定义为： 其中|D|为状态集合D中所有状态的个数。也称pij为一步转移概率，其构成的矩阵称为一步转移矩阵。

Aij(t)称为接受概率，表示从当前在状态i时产生状态j接受状态j的概率，Metropolis准则中，接受概率定义为： Gij(t)称为从状态i到状态j的产生概率，表示当前在状态i时，状态j被作为下一个转换状态的概率，若状态j为状态i的邻居，则状态j被选中的概率为 N(i)为状态i的邻域。 Aij(t)称为接受概率，表示从当前在状态i时产生状态j接受状态j的概率，Metropolis准则中，接受概率定义为： 注意这里概率矩阵pij完全描述了各个状态之间的转换

上述模型中一步转移概率矩阵刻画了系统的粒子在各种状态之间的转换规律，可以证明，上述模型满足： （1）可达性：即任何一个状态都可以达到，这样就有了求全局最优解的可能； （2）渐进不依赖起点：即无论起点如何，都能达到全局最优解； （3）分布稳定性：包含两个内容，一是温度不变时，极限分布存在，二是当温度趋于0时，极限分布存在； （4）收敛到最优解：当温度趋于0时，最优状态的极限分布概率为1。 这就从理论上论证了模拟退火算法的有效性。但算法真正实现起来，还有几个问题需要解决

算法运行的三个关键技术问题： （1）温度如何下降； （2）新解如何产生； （3）参数如何设置。 有哪些参数？ 等温过程中，也就是内循环做多少次合适？初始温度多去多少？不管是内循环还是外循环什么时候停止？

温度下降控制 降温方式对算法性能影响很大，降温过快，可能丢失最优值点，而降温过慢导致算法收敛速度极大降低，同时为保证全局收敛性，一般要求温度下降趋于0。 （1）经典退火方式: T(t)=T0/ln(t+1)。降温过程非常慢，因此算法收敛性较慢。 （2）快速退火方式: T(t)=T0/(1+ɑt)。降温过程在高温区速度较快，低温区速度逐渐减慢。ɑ可用来改善降温曲线。 （3）指数衰减 : T(t+1)=ɑtT(t)，这里ɑt接近1，一般取0.8-0.99

新解产生：分为如下四个步骤： 第一步：由一个产生函数从当前解产生一个位于解空间的新解，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等；或者对当前解添加随机扰动得到新解。 第二步：计算与新解所对应的目标函数差。 第三步：判断新解是否被接受, 最常用的接受准则是Metropolis准则: 若ΔE<0则接受S′作为新的当前解S，否则以概率exp(-ΔE/T)接受S′作为新的当前解S。 第四步：当新解被确定接受时，用新解代替当前解，在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上进入下一轮迭代。

此时每个状态出现的概率相等，即为1/|D|。根据上式可以得到T0的一个估计 初温越高，获得全局最优解的可能性越大，但相应所需的计算时间将增加，效率将降低。从理论上来说，初始温度应使得下式尽量满足才能更好的保证全局最优性。 此时每个状态出现的概率相等，即为1/|D|。根据上式可以得到T0的一个估计 第一个公式说明的是初始情况下，从i状态跳到任何一种状态都是等可能的 实际不知道 K是充分大的数，可以通过试验得到。

（1）随机产生一组状态，确定两两状态之间的最大目标差Δmax，依据差值，再利用一定的函数确定初温，如T0= Δmax /Pr，其中Pr为初始接收概率； （2）采用统计推断得到最大最小目标函数差。随机采样一组样本（函数值），假设函数值服从正态分布，计算样本函数值的均值μ和方差σ ，按照3σ原则，目标值99.97%落在[μ - 3σ , μ +3σ]之间，所以函数值极差为 6σ ；

（3）数值估计 第一步：给定一个常量T；初始温度T0; R0（接近1，如0.9,0.8）,r=0,k=1； 第二步：对于一个给定迭代步长L，进行蒙特卡洛模拟，记录算法中被接受的状态和被拒绝的状态个数，计算接受状态与迭代步长L的比值rk； 第三步：若|rk-R0|<ɛ，停止计算，否则若rk-1<R0,rk<R0,则k=k+1,T0=T0+T，返回第二步，若rk-1和rk>=R0，则k=k-1,T0=T0-T，返回第二步，若rk-1>R0,rk<R0，则k=k+1,T0=T0+T/2，T=T/2，返回第二步，若rk-1<R0,rk>R0，则k=k-1,T0=T0-T/2，返回第二步 R0表示期望开始时候的接收概率，一般来说要大一点，下面我们来模拟查看接收概率的情况 通过这个过程得到合适的初始温度

参数设置-内循环（保温过程）迭代次数L设置 （1）固定长度：通常选取与邻域大小直接相关的量。有研究表明一般是邻域的平方级别。 （2）由接受比例和拒绝比例来控制迭代次数。 (1) 有研究表明一般是平方级 （2）当温度很高时，迭代次数应使得保持接受比例和拒绝比例基本相同条件下次数最少，而当温度较低时，为了避免过早陷入局部最优状态，迭代步长应适当增加。一种方法是可以设定一个迭代步长上限和下限以及接受比例R，每一温度至少迭代下限次，若没达到接受比例，则增加迭代步长。

内循环（保温过程）终止条件 （1）直接根据迭代次数终止内循环； （2）每个温度下只计算一次，但控制温度下降非常慢； （3）内循环中最优解持续不再发生变化。

参数设置-外循环（降温过程）终止条件 （1）通常可以使最终温度降到0，但这样往往使算法运行时间过长； （2）实际上温度一般不需要降到0，因为差解的接受概率已经非常接近温度是0的时候概率，所以，停止条件可以是 一个合适的较低的温度； 系统进入“稳态”，即既没有好解也没有坏解被接受或者称为接受概率控制法，若接受概率大于一定程度，如1/N 则停止 确定的迭代次数

算法比较 （1）局部搜索法 （2）爬山法 局部搜索 梯度下降法 模拟退火算法

比较 模拟退火 遗传算法 迭代变量 一个解 种群 新解产生 邻域随机搜索 交叉、变异 最优性 物理机理 Metropolis准则 选择算子 全局性 蒙特卡洛仿真 种群多样性 迭代控制 温度控制 进化代数 优点 优化准则随温度自适应 策略简单、易操作 不足 参数多，难控制 自适应不足

温度的控制管理和种蒙特卡洛仿真的有效性是模拟退火算法得以获得好效果的关键。 模拟退火算法小结 模拟退火算法通过对温度参数的控制基于局部随机仿真实现目标寻优，其蕴含深刻的思想：迭代寻优、局部随机仿真、解的稳态。其基于的Metropolis准则使得算法既反映出全局搜索能力又有局部搜索能力。 温度的控制管理和种蒙特卡洛仿真的有效性是模拟退火算法得以获得好效果的关键。 遗传算法中种群的引入是非常重要的，没有种群，全局搜索能力就没有保障

例1：求一元函数 在[-1,2]区间内的最大值 此函数在 [-1,2]有多个局部极值点，其全局最大点在约1.85处达到，最大值为3.85

线性温度下降

指数温度下降

对数温度下降

例2：求二元函数 在[0,5]区间内的最大值 此函数在[0,5]有多个局部极值点，其全局最大点在约（3,3）处达到，最大值为1

线性温度下降

指数温度下降

对数温度下降

例3：求二元函数 在[-5,5]区间内的最小值 此函数在[-5,5]局部极值点更多，其全局最大点在（0,0）处达到，最小值为0

线性温度下降

例4：十滴水游戏遗传算法求解 游戏规则：6*6的棋盘格中分别有0-4滴水，当前你有10滴水，可以加入到棋盘格中，若加入后棋盘格中水大于4滴，则将破灭并向上下左右四个方向各弹出一滴水，同样可能引起其他棋盘格中水破灭，若没破灭，则棋盘格中多一滴水。

算法设计 解产生：考虑到棋盘的特殊形式，棋盘位置只能为整数，因此此问题解不能直接使用连续实数变量表示，这就需要设计新解产生方式 设向量X=(x1,x2,…,x10)表示点击序列，xi表示一个整数对应棋盘的位置，如7对应的是棋盘第二行第一列，对于当前解，新解为X’=(x1’,x2’,…,x10’) 若产生新解的xi’,小于等于0或大于等于36，则 需要取整数

能量函数 该问题的能量函数不易直接表示出来，如目标希望是所用的水滴数最少，那么目标需要根据给出的序列模拟算出其清空棋盘所有水滴需要使用的最少水滴数。若10滴水全用完仍然未清空，则返回10。这个能量函数无法直接给出来，需要模拟程序运行。 温度下降 采用线性温度下降。 终止条件 采用最优解温度作为终止条件。

作业： 14.1 P369 第3题

谢 谢！