案例详解MATLAB实现模拟退火法

模拟退火法作为一种适合于求解大规模的优化问题的技术，近年来已引起极大的关注。特别是当优化问题有很多局部极值而全局极值又很难求出时，模拟退火法尤其有效。在实用上，它有效地“解决了”著名的旅行推销员问题，即在必须依次访问每一个城市（共有N个城市）的前提下，为旅行推销员设计一条能够返回起点的最短旅程。模拟退火方法还被成功地用于设计复杂的集成电路，也就是说如何最佳地安排几十万个电路元件，使它们全部集成在一个很小的硅片上，而相互连接的线路之间的缠绕能够达到最小。尽管模拟退火法的功效非凡，但它的算法实现却相对地简单，这一点似乎有些不可思议。

注意：上面提到的两个问题都属于组合极小化问题。本类问题一般也有一个目标函数，但是函数的定义域并不是简单地由N个连续参变量组成的N维空间，而是一个离散的巨大空间，例如，由所有可能的城市旅行路线组成的集合，或者硅片电路元件所有可能的分配方式的集合。构型空间中元素的数量相当巨大，根本不可能穷举，而且因为集合是离散的，也不可能“沿合适的方向连续下降”。因此在构型空间中，“方向”概念就没有什么意义了。

在连续控制参数的空间中利用模拟退火法，这种应用实际上要比组合问题复杂一些，因为其中又要出现“狭长山谷”的情况。正如在下文中将看到的，模拟退火法的试探步骤是“随机”的。

但在一个狭窄且漫长的等高线山谷中，几乎所有的随机步骤都呈向上的趋势，因此，算法中需要增加一些技巧。

模拟退火的核心思想与热力学的原理颇为相似，而且尤其类似于液体流动和结晶以及金属冷却和退火方式。在高温下，一种液体的最大分子彼此之间进行着相对自由的移动。如果该流体慢慢地冷却下来，热能可动性便会消失。大量原子常常能够自行排列成行，形成一个纯净的晶体，该晶体在各个方向上都被完全有序地排列在几百万倍于单个原子大小的距离之内。对于这个系统来说，晶体状态是能量最低状态；而所有缓慢冷却的系统都可以自然达到这个能量最低状态，这可以说是一个令人惊奇的事实。实际上，如果某种液体金属被迅速冷却或被“猝熄”，那么它不会达到这一状态，而只能达到一种具有较高能量的多晶状态或非结晶状态。

因此，这一过程的本质在于缓缓地致冷，以争取充足的时间，让大量原子在丧失可动性之前进行重新分布。这就是所谓退火在技术上的定义，同时也是确保达到低能量状态所必需的条件。

1953年，米特罗波利斯（Metropolis）及其合作者们首次将上述原理渗透到数值计算中。他们对一个模拟热力学系统提供了一系列选择项，并假设：系统构形从能量E1变化到能量E2的概率为p=Exp[-（E2-E1）/kT]。很显然，如果E2<E1，p将大于1；在这类情况下，对构形的能量变化任意指定一个概率值p=1，也就是说，该系统总是取这个选择项。这种格式总是采取下降过程，偶尔采取上升步骤。目前已被公认为米特罗波利斯算法。

为了将米特罗波利斯算法应用于热力学以外的系统，必须提供以下几项基本要素：

（1）对可能的系统构型的一种描述。

（2）一个有关构型内部随机变化的生成函数，这些变化将作为“选择项”提交给该系统。

（3）一个目标函数E（类似于能量），求解E的极小值，即为算法所要完成的工作。

（4）一个控制参数T（类似于温度）和一个退火进程，该进程用来说明系统是如何从高值向低值降低的，例如在温度T时每次下降步骤中要经过多少次随机的构型变化以及该步长是多大等等。应说明的是，这里“高”和“低”的含义，还有进程表的确定，都需要一定的物理知识和一些摸索的实验。

栗子：　组合极小化：旅行推销员问题

下面以“旅行推销员问题”为具体实例说明模拟退火法的应用。假设一个推销员，要去N个分别位于（xi，yi）的城市进行推销，并于最后返回他原来所在的城市，要求每个城市只能去一次，而且所经过的路径要尽可能地短。这个问题属于一类所谓“NP-完全问题”。这类问题求出一个精确解所需的计算时间是随N的增加以指数exp（常数×N）增长的。当N不断增大时，运行时间将迅速增加，进而导致费用高到令人难以接受的程度。旅行推销员问题也是众多极小化问题中的一种，它的目标函数E具有多个局部极小值。在实际应用当中，常常有足够多的条件可以从多个极小值中选出一个最小的，这个最小值即使不是绝对最小，也相当接近绝对最小值了。退火法的目的就是要获得这个最小值，同时又要将计算量限制在N的低阶次的数量级上。

旅行推销员问题也是按模拟退火问题的方式进行处理的。具体如下：

（1）构型　将N个城市分别标记为i=1，2，…，N中的数，其中每个城市具有坐标（xi，yi），见表11-2-1。一个构型就是数字1，2，…，N的一个排列，可以解释为推销员途径的城市的顺序。

（2）调整　林（Lin）曾经提出过一种所谓“转移有效集”，这里的“转移”包括两种类型：（a）移走路径的某一段，然后对这段路径上的城市用相反次序重新进行排列，并用后者来代替前者。（b）移走某段路径，并用位于城市间的随机选取的另一段路径来取代被移走的路径。

（3）目标函数　在旅行推销员问题的一种最简单的形式中，目标函数E被定义为旅途的总长度

表11-2-1　52座城市的坐标

这里认为第N+1个点与第1个点是重合的。

11.2.2　MATLAB求解

1.TSP算法流程

根据以上对TSP的算法描述，可以写出用模拟退火算法解TSP问题的流程图11-2-1所示：

2.加载数据文件

52座城市的坐标数据文件li11_2_1.m

图11-2-1　TSP的模拟退火流程

当调用数据文件函数时，包含城市坐标信息的矩阵载入到inputcities的变量中。该数据文件的第一列通常是城市数量，之后两列为城市的坐标。初始化函数，根据直角坐标生成城市距离矩阵。

3.计算总距离的函数

这是一个城市间计算距离的函数，根据给定路径计算该路径对应总路程。

TSP问题的成本函数是城市之间的距离。调用此函数将计算n个城市之间的距离。

4.绘制路线的函数

这是一个绘制路线的函数，根据给定的城市坐标生成的城市矩阵绘制城市的坐标，根据给定的路线绘制线路。

5.交换城市的函数

这是一个用于城市交换的函数，它从某路径的邻域中随机的选择一个新的路径。

6.执行模拟退火的函数

7.用模拟退火算法解TSP问题的操作步骤

①将2～6的5个程序按对应程序名存盘；

②编写li11_2_1zhu.m程序代码存盘；

③在命令窗口执行li11_2_1zhu.m程序。

52个城市的TSP运行结果如图11-2-2。

图11-2-2

该优化路径的总路程近似为7.5946e+003。这是由52个城市构成的一个对称TSP实例，利用上述算法对该实例进行模拟退火求解，设定初始温度T0=99，冷却速率为0.999，互换城市数为2。经过20次仿真，得到的最优解与已知最优解（7.5444e+003）非常接近，所需时间也令人满意。