导引模拟退火法.docx
《导引模拟退火法.docx》由会员分享，可在线阅读，更多相关《导引模拟退火法.docx（20页珍藏版）》请在冰豆网上搜索。
导引模拟退火法
简介导引模拟退火法及其应用
李世炳1邹忠毅2
中央研究院物理研究所
1e-mail:
spli@phys.sinica.edu.tw
2e-mail:
cichou@phys.sinica.edu.tw
1.前言：
最佳化问题实际上是一门非常古老的学科，它存在于各行各业中，是一门跨领域的学科。
其中一个最有名的例子是所谓的旅行商问题（TravelingSalesmanProblem）。
在这个例子中，该名旅行商需要跑遍N（N大于一）个城市去推销他的商品，而该些城市之间的距离都不一样，这名推销员需要从其中一个城市出发，而他老板规定他必须把所有城市跑过一遍，请问这名旅行商应该如何绕才是最省时间（假定他的速度一直不变），也就是说，这名旅行商要找出一个最短距离的路径。
这个问题看起来简单，实际上是一个非常复杂的数学问题。
譬如说，现在只有两个城市，如果这名旅行商从其中一个城市出发，那这里只有一个可供这名推销员可选择的路径。
假如现在有三个城市，路径还是只有一个，因为他从其中一个城市出发绕一圈跟他反方向绕一圈的路径距离是一样的。
意思是说假如该三个城市分别是A，B，C。
他从A出发，先到达B然后C，再回到A，跟他先到达B，然后C再回到A所经过的距离是一样的。
假如现在有四个城市A，B，C，D，那他就有三个选择，他必须先把三个不同的路径的距离算出来后再决定要选择那一个。
因此城市的数量愈多，可能的路径也愈多，而且增加的速度是非线性的，十个城市的所有可能路径就会有十八万一千四百四十个之多，假如他每一个路径的距离都要先算出来后再作选择，那他所需要的时间简直就是天文数字，在实际的情况，他根本没有那么多的时间去做这件事，而设计一个有效的方法去找寻这个最短的路径就是所谓最佳化问题的基本精神。
以上所介绍的旅行商问题（TravelingSalesmanProblem），是属于数学上所谓的NP问题。
这类问题之所以被称作NP问题，是因为随着问题中的变量的逐渐增加，它的复杂度却以惊人的速度在增加。
在过去人们对这类问题根本是束手无策，但随着计算机的出现，人们可以发展更有效的算法，并且在计算机上仿真，使得寻找这类问题的最佳化答案变成为可能。
最早期发展中的一个算法是人们称之为蒙地卡罗（MonteCarlo）算法，是Metropolis[1]等人在一九五三年提出。
虽然已经有了五十年的历史，但是因为它的简单性跟实用性，它还是被广泛运用在各种最佳化问题中。
它主要的概念是把所要研究的问题看成是一个统计系统，而一个统计系统的某一个温度是的状态分布是满足一个波兹曼（Boltzmann）概率分布函数。
因此在问题中寻找最佳化解答时，就是利用这个分布函数来选取答案。
譬如说在前面的旅行商问题中，我们先随意的把这N个城市连起来，如果N不是一个小数值，那我们刚刚连起来的路径几乎不可能是最短的路径。
在蒙地卡罗算法中，我们就从这个路径开始，来寻找最短的路径。
我们可以在这N个城市中随机的选取其中两个A跟B，并且把它们交换，意思是说我们的旅行商本来要先到达A城推销他的商品，然后经过若干城市后再到B城推销他的商品。
现在却是先到达B城推销他的商品，然后经过若干城市后再到A城推销他的商品。
经过这样的交换后，我们在把新的路径连起来，假如这个新的路径比旧的来得短的话，我们就把新的保留，但假如这个新的路径比旧的来得长的话，我们还是有一个概率把新的保留下来，而要不要保留这个新的路径就是用前述的波兹曼概率来决定。
模拟退火法（SimulatedAnnealing）是Kirkpatrick[2]等人在一九八三年提出并成功地应用在组合最佳化问题中，它是蒙地卡罗算法的推广。
退火是一种物理过程，一种金属物体再加热至一定的温度后，它的所有分子在状态空间中自由运动。
随着温度的下降，这些分子逐渐停留在不同的状态。
在温度最低时，分子重新以一定的结构排列，而分子的分布也就是以前面所述的以波兹曼（Boltzamnn）概率分布。
不同于上述所介绍的蒙地卡罗算法，模拟退火法中的温度是随着退火的时候有所改变，因此如何对温度作有效的调整就变成整个模拟退火法最重要的一环。
以上所介绍的算法都属于同一类型的算法，是所谓的局部搜索方法（localsearchmethod）。
这类问题虽然在局部搜索的能力很强，但是对于全局搜索（globalsearch）的能力却嫌不足，因此研究学者也开始开发具有全局搜索能力的算法，其中一个最早诞生且被广泛应用的算法就是所谓的遗传算法（GeneticAlgorithm---GA）。
遗传算法早期的研究工作始于二十世纪六十年代，在五十年代末六十年代初，一些生物学家开始利用计算机对遗传系统进行仿真。
在此期间，受到生物学家们模拟结果的启发，Holland和他的学生们首次应用模拟遗传算子来研究适应性中的人工问题。
之后，在六十年代中期，Holland[3]开发了一种数值模拟技术---遗传算法，其基本思想是利用类似于自然选择的方式来设计一套最佳化演算程序。
在随后的十年中，Holland致力于创造一种能表示任意计算机程序结构的遗传码，以拓展遗传算法的应用领域。
遗传算法是仿真前述生物进化过程的计算模型，它所处理的是染色体（chromosome），或者叫基因行个体（individuals）。
一定数量的个体组成了群体（population），也叫集团。
群体中个体的数目称为群体的大小（populationsize），也叫群体规模。
而各个体对环境的适应程度叫做适应度（fitness）。
遗传算法基本上是一种群体型操作，该操作以群体中的所有个体为对象。
选择（selection）、交叉（crossover）和突变（mutation）是遗传算法的三个主要操作算子，它们构成了所谓的遗传操作（geneticoperation），使算法具有了其它传统方法所没有的特性。
遗传算法中包含了五个基本要素：
（一）参数编码﹔
（二）初始群体的设定﹔（三）适应度函数的设计﹔（四）遗传操作设计﹔（五）控制参数设定（主要是指群体大小和使用遗传操作的概率等）。
这五个要素构成了遗传算法的核心内容。
经过了多年的研究，遗传算法的应用研究已经从初期的组合优化求解拓展到了许多更新、更多元化的应用方面。
并且与其它的进化规划（EvolutionProgramming，EP）以及进化策略（EvolutionStrategy，ES）等进化计算理论日益结合。
尽管遗传算法比其它传统搜索方法有更擅长全局搜索的能力，但是它的局部搜索的能力却嫌不足。
其它如禁忌搜索（TabuSearch）[4]，神经网络（NeuralNetworks）[5]，进化算法（EvoluationaryAlgorithms,EA）[6]等等，都是过去二三十年间研究人员开发出来的一些算法。
限于篇幅的关系，本文不作一一的介绍，有兴趣的读者可以从参考文献中得到相关的讯息。
2.向自然学习──两大类随机搜寻法的操作步骤
有趣的是，以上目前较成功的两大类型的随机搜寻法，都正好是由尝试模拟大自然中的现象而发展出来的。
第一个现象是物质世界中材料的冷却与结晶过程，例如一些高温融化的材料，在缓慢的降温冷却下，可能形成结构整齐的晶体。
另一种现象是生物世界的演化过程，在长时间及多次的世代交替后，能够在演化过程中成功存留的物种，也常常是较能适应环境的物种。
观察以上成功的两个真实世界中的最佳化过程，让人不得不佩服大自然的巧妙。
科学家们就向以上两个自然现象学习，发展出了两大类型的算法，就是模拟退火法（simulatedannealingmethod）与遗传算法（geneticalgorithm）。
以下我们就来介绍如何将这两种算法的基本想法化为计算机上的实际操作步骤。
2.1物质材料的冷却与结晶过程──模拟退火法
首先介绍模拟退火法，这个算法被希望能满足两个条件，第一，当温度够高时，系统的组态要能自由变化，也就是这系统的组态能在能量表面自由的移动或称为它在做无规行走（randomwalk）。
第二，当温度变小时，系统的组态在能量表面的移动将受到限制，并逐渐地向低能量的区域集中。
模拟退火法简单的操作方法如下：
（1）首先要针对问题选定目标函数，并且将这个目标函数视为一个广义的能量函数E（x）。
（2）引入一个温度T，这个温度T不一定要具有实际的物理意义。
并选择一个足够高的温度作为起始温度。
（3）在这一个温度作运算来仿真热平衡过程，通常使用metropolis方法来模拟热平衡过程，方法如下：
1.计算目前系统组态的能量E，随机的将系统组态改变一点，并计算变动后的系统组态的能量E’。
2.计算变动成功的机率P=min（1,exp（-（E’-E）/kT）），并以此为依据来决定系统组态是否改变。
（4）订定一个退火策略（或称为降温程序annealingschedule），由此决定在每一个温度所停留的时间，及降温的比例。
使用以上的演算步骤，科学家们已经在许多不同的问题上得到了成果。
但是对于较复杂的问题，要想成功的求得基态解，关键在于如何选好退火策略。
如果我们真的要模拟大自然成功的退火过程，往往要花上极大的计算时间。
对于一些多变量问题，我们目前所有的计算机计算能力，根本没有办法在短时间（一星期，一个月甚至一年）内完成。
所以有许多科学家们也正在尝试改进模拟退火法，来加强它的计算效率。
2.2生物系统的演化过程──遗传算法
接下来介绍遗传算法。
这个算法被希望能满足几个条件，首先生物（我们的各个系统组态）能透过交配与突变等种种方式把自己的一些特征遗传给下一代（新的系统组态），而在一代一代的遗传过程中透过环境（也就是我们的目标函数）的天择，使得后代适应环境的能力越来越强（也就是接近目标函数的最佳值）。
透过以下的方式，可以实践前面的想法：
（1）首先要针对问题选定目标函数（适应度函数），并且将这个目标函数作为一个适应环境的能力的指标，并设定控制参数。
（2）透过一定的编码技术，将系统组态编成一组基因码。
（3）先随机产生一群基因码，作为第一个亲代群体。
（4）透过以下三种遗传算子，完成世代交替及天择的演化过程：
1.交叉算子（crossover）──将亲代两个基因取出，透过交换部分基因码的方式，产生子代。
2.突变算子（mutation）──将亲代的部分基因码随机改变，产生子代基因。
3.选择算子（selection）──从子代基因中选择出较优秀的基因，使它们有机会参与下一代的繁殖过程。
用前面所提到的目标函数来决定基因的优秀程度
（5）反复进行上个迭代程序，直到迭代次数到达限定值或目标函数值收敛到某一个极值。
虽然遗传算法也已经应用在许多不同的问题上并且得到了一定的成果。
但是遗传算法的困难之处，在于以上所提到的各种算子，在使用时都有许的不同的控制参数。
对于较复杂的问题，我们很难选择出适当的控制参数，来有效地找出最佳解。
3.新的最佳化算法──导引模拟退火法
3.1就是没有好的模拟退火法的退火策略！
──新的算法反而由此产生
在介绍了以上两大类现有的较有效算法之后，接下来介绍一种由我们研究团队所发展出来的较新的最佳化算法──导引模拟退火法。
我们在2000年中以前，一直在尝试如何改良传统的模拟退火法，使它能更有效率地应用在多变量系统的求解基态问题上。
但是我们发现，要使模拟退火法一次就求出整体的最佳解是非常困难的。
在一次次的模拟计算中，绝大多数的试验结果总是得到能量较高的亚稳态解。
可是我们观察这些能量较基态为高的亚稳态解时，却发现到这一些亚稳态解之间有一些关联性。
所以我们便有了一些想法，并藉此发展出了一个最佳化算法。
我们的基本想法有三点：
第一，我们不再期望能由传统的模拟退火法一次就找到基态解,改为在较短时间内用它来找到许多亚稳态解。
第二，藉由分析这些亚稳态解的性质,找到有用的信息。
第三，想办法利用这些信息，来帮助下一代的搜寻过程。
这些想法使我们的算法不仅包含了传统模拟
退火法单纯地模拟物质世界中材料的冷却与结晶过程的表现，更具有了生物世界世代交替逐渐演化的精神。
接下来是要如何用数值方法来实践我们的基本想法？
首先我们要使用统计方法来分析亚稳态解的性质。
我们将这些亚稳态解做出统计分布图也就是所谓的直方图（histogram）。
再来观察亚稳态解的histogram可以知道这些解的可能出现区域，并由这些解的共有性质来去除不相关的自由度。
接下来我们在下一代的搜寻过程中增加histogram中，高可能性区域的搜寻机会，如此便可能找到更低能量的解。
我们用以下的示意图来说明，图
（一）中圆点代表在能量曲面中亚稳态解的位置，然后将亚稳态解做出统计分布图，藉由增加搜寻分布图高可能性区域，在第二代中找到基态。
当然这只是最简单的示意图，实际的问题通常不能简单地将高维度的能量曲面及分布图绘出来。
图
（一）导引模拟退火法的基本想法
由以上简单的示意图，我们也可以看出，这个新的算法的关键，就在于如何用上一代有用的信息来导引下一代的搜寻过程。
对于简单的问题，我们当然可以用个别变量位置的histogram来引导下一代的搜寻。
但是对于各种不同的较复杂的问题，上一代所共有的有用的信息，可能不仅仅是个别变量位置那么直观。
为了使我们的新算法，能更广泛的适应各种不同的问题，我们把这种用来引导下一代的搜寻过程的，上一代所共有的有用信息，以一种函数形式来表现。
我们称它为导引函数（guidingfunction），而这个新算法，我们就称它为导引模拟退火法（GuidedSimulatedAnnealing）。
3.2一个简单范例──Bergequation.
接下来用一个简单范例来介绍导引仿真退火法的操作方式。
考虑一个简单的能量函数[7]，
，
图
（二）简单的能量函数
其中，xi的范围在–1到1之间。
A=10，B=0.49，C=1.0E-5。
图
（二）表示xi与Ei的关系。
由图，我们可以看出，对于Ei有两个差距非常小的极小值（约1.4E-5），落在xi约等于–0.7与0.7处。
而整个能量函数E的整体最小能量也就落在所有的x都约等于–0.7处，并且能量函数E的所有的极小值一共有2的N次方个。
这个函数看起来简单，可是要用随机数值方法找出它的整体极小值，却不是一个简单的事情。
那么就来介绍如何使用导引模拟退火法来找出它的整体极小值。
首先我们决定相关参数，以100个变量为例（N=100，能量函数E一共有2100个极小值）。
决定每一代群体大小为20，退火过程每个温度有效移动100次，降温率0.85。
接着我们来执行导引模拟退火法的操作步骤：
第一步：
利用传统模拟退火法得到第一组亲代。
由于我们用了很快速的退火策略，所以第一组亲代群体，没有找到整体极小值。
第二步：
分析第一组群体的变量分布图作为导引函数。
对于这个简单的问题，我们可以直接用个别变量位置的histogram来作为导引函数（以（x）表示）。
第三步：
利用加入导引函数的退火过程得到子代。
有了导引函数后，我们就可以用它来帮助下一代群体的产生。
导引函数的使用方法如下，
（1）与传统模拟退火法不同，我们不再随机的选择变量并随机的改变变量。
在导引仿真退火法的操作中，变量的选择机会与变量所具有的导引函数成反比，也就是在许多变量中，越与上一代整体特征不符合的变量越容易被挑选到。
这个步骤也可以用下式表示
，
其中
表第i个变数被挑选的机率。
（2）挑选了变量之后，变量的改变也与该变量的导引函数成正比，也就是该变量较会被改变到导引函数值高的区域中。
用下式这个步骤表示
其中
表变量i被移动到位置xi的机率。
（3）接下来与传统模拟退火法相同，用波兹曼概率分布函数来决定移动成功与否。
第四步：
分析子代群体的变量分布作为新的导引函数。
然后并回到前一步，直到世代次数到达限定值或目标函数值收敛到某一个极值。
我们取其中的一个变量为例，图（三）可以看出第一代的解的分布集中在两个的极小值附近。
随着世代的增加，导引函数中解的分布越来越集中在两个的极小值附近。
同时这两个分布的大小，也随着世代的增加逐渐显出差距。
最后导引函数中只剩下一个尖锐的峰值，整体极小值也找到了。
另外我们也将导引模拟退火法的操作流程图展示如图（四）。
3.3导引模拟退火法的特点及限制
以上用了一个简单的例子来示范导引模拟退火法的操作过程。
接下来我们来看看这个方法的一些特色。
首先，我们引入了一个新的性质，也就是导引函数，它使搜寻过程中，系统组态较有可能在能量区面上，数个低能量区域间变换，使得系统组态也较不易陷入局部极小。
其次，导引函数可以有效地去掉不相关的自由度，使搜寻范围逐渐减少，进而
图（三）导引函数的变化
加速搜寻效率。
又因为我们所得到的每一个子代都是由各自独立的退火过程产生，使得这个算法可以容易地在分散或平行化计算机设施上执行，可大幅加强计算效率。
对于某些问题，我们可能可以用其它的方法预测解的大致范围。
也就是说，已经有了一个粗略的导引函数。
在这种情形下，导引模拟退火法很容易与其它方法配合使用，例如用其它数值或实验方法先取得粗略的导引函数，再用本法做进一步的计算，找出更精确的解。
另外这个方法在漏斗形的能量区面问题上特别有用。
这是因为这种能量区面，亚稳态解之间的关
图（四）导引模拟退火法的操作流程图
联性较强。
我们很容易就统计出亲代共有的特性，到到很好的导引函数来帮助下一代的搜寻。
在讨论了导引模拟退火法的一些特色后，我们再来看看和这个方法有相关性的两种方法，也就是模拟退火法与遗传算法与导引模拟退火法的异同处。
首先是模拟退火法，这两个方法表现出来的，都是使系统组态在能量曲面上做无规行走（randomwalk），并且透过波兹曼概率分布函数对能量的选择，使系统组态逐渐变化到能量最小值处。
但是导引模拟退火法多了一个方向性，使得本方法除了保留了原先模拟退火法良好的局部搜寻能力，更借着导引函数的帮助，加强了这个方法的全局搜寻能力。
接下来来讨论导引模拟退火法与遗传算法的异同处。
明显可见，这两个方法都是世代演化形式的演算，并不指望在第一代搜寻中便找到基态解。
而是透过世代演化的计算过程，逐渐地找到最好的解。
而且，在这两个方法的演算过程中，子代与亲代也都有一定的相似性。
可是对于导引模拟退火法的计算过程，子代不是仅由一或二个的亲代产生，它包含有上个世代整体特征。
或是说，亲代的整体特征透过导引函数遗传到了子代。
另外，由于在产生子代的过程中，使用了退火过程，也使得本方法比遗传算法具有更好的局部搜寻能力。
虽然导引模拟退火法有以上的一些优点，可是它也有以下的限制。
首先是每世代的群体规模大小。
由于我们必须统计分析每一代群体的共有特色，所以群体的规模不可以太小，可是过大的群体规模，又影响到计算时间。
其次，也是因为使用退火过程来产生子代。
所以我们还是必须要有一个退火策略，这个退火策略虽然不必像模拟退火法一样，需要极为缓慢的降温程序，但也还是不能太急促的降温。
最后，由于导引函数大大影响计算能力。
如果我们对于所面对的问题越了解，我们就越有可能写出较有效率地收集亲代群体共有特色的导引函数，进而加速计算。
4.导引模拟退火法的应用
在这一节，我们要介绍导引模拟退火法的实际应用。
以下的不同类型问题，是我们研究团队在目前所得到的一些结果。
我们希望透过用不同类型题目的测试，增进我们对这个新方法的了解，进而改进这个方法。
同时也希望这个方法能对这些问题所关联的研究领域有所贡献。
但是由于篇幅所限，我们只能够对各个问题做最简单的介绍，介绍的重点也偏重在导引模拟退火法的运用上，至于各个问题的深入讨论只好割爱了。
4.1与实验结果的配合──X光结晶学的运用
首先介绍导引模拟退火法如何与实验结果的配合，以一个X光结晶学的问题问题为例。
由于X光对晶体绕射实验一般只能观察到绕射光点的强度，所以要从绕射结果反推导出晶体的结构就变成了一个不容易的问题。
一般来说，科学家可以用所谓的directmethod等等方法[8][9]来解出晶体的结构。
可是我们也可以将这个问题用更直观的方式表现。
我们可以随意写出一个晶体结构，再将这个假设的晶体结构用计算机来计算出它对应的绕射光点的强度，如果它的绕射光点强度与实验所得到的强度相同，我们所假设的晶体结构可能就与真的实验所用的晶体结构相同[10][11]。
以上的想法再配合最小平方近似的技巧，可以用以下的能量函数表示，
，
其中
为第i个绕射光点的方向，
为假设的晶体结构的绕射光点强度，
为实验所观测到的绕射光点强度。
所以这个问题就变成，如何选择一个晶体结构使以上的能量函数变为最小。
图（五）Isoleucinomycin的结构
以下介绍导引仿真退火法运用在晶体结构问题的结果[12]。
以一个用传统模拟退火法不好解的大分子Isoleucinomycin（C60H102N6O18）[13]为例（图（五））。
我们使用原子在单位晶胞中空间分布作为这个问题的导引函数。
图（六）展现出各世代导引函数变化的结果，本图将立体的导引函数投影在二维平面上，其中深色代表高可能区域，表示建立导引函数时所用的取样宽度，并随着世代增加也逐渐缩小取样宽度。
从图中，一开始我们只能大概看出各分子的中心位置。
随着世代增加，可以看出导引函数逐渐收敛，终于解出所有原子位置。
4.2应用在原子丛集问题
其次介绍导引模拟退火法如何应用在原子丛集（atomiccluster）问题。
简单来说，假设知道两个原子之间的作用位能V（r1-r2），我们就可以算出两个原子的最低能量结构。
同样的，三个原子的最低能量结构，我们也许还可以用人工的方式算出。
随着原子数量增加，人工的方式越来越不可行。
那么，N个原子丛集的最低能量结构是什么[14][15]？
圖（六）導引函數的變化
以一个简单的原子之间的作用位能Lennard-Jonespotential为例，令它的能量函数表示为
其中rij为第i个原子与第j个原子的距离。
对于导引函数的取法，我们最直观的想法，也就是用原子的空间分布做为我们要的导引函数。
可是当我们在用导引模拟退火法计算原子丛集问题时，我们观察到在计算时各代群体的结构有一些对称性，这些群体的结构可以分为几个对称形状（如图（七））。
如果我们利用这些对称性，将导引函数再细分为几个不同对称性的大类，便可以大大加速计算过程。
我们计算了N不大于150个LJ原子的最低能量结构[16]，并得到与目前已知最好结构相同的结果[17]。
图（八）LJ27的结果
图（七）LJcluster常见的对称形状
图（八）表示27个LJ原子的各世代所找到的最低能量结构，由图可知在第五代时已找到类似最低能量结构的组合，但是一直到第九代才分辨出基态与次低的亚稳态的差别。
图（九）表示37个LJ原子的各世代平均能量与温度的关系，随着世代的增加，平均能量曲线逐渐下降。
图（九）LJ37的各世代平均能量与温度的关系
4.3离散系统的应用
以上两个问题，它们的变量都可以连续地变化。
接下来我们介绍另一大类问题，这类问题的变量不可以连续地变化，也就是所谓的离散系统。
首先是前言中所提到的旅行商问题，它的定义已在前面说明。
它的目标函数（或能量函数）也就是旅行商行走距离的路径。
对于这个问题，我们取导引函数为各城市连接频率的分布。
我们希望在用导引模拟退火法计算过程中，各新的世代在产生时能尽量保有前个世代中最常连接的一些城市的相对关系[18]。
用这个方法，我们得到了一些初步的结果，表
（一）表示我们的结果与世界最好纪录（由数据库TSPLIB中可得）的比较[19]，图表中的结果是我们在每一代中只用了一百五十个蒙地卡罗步所得到的，假如在每一代中增加蒙地卡罗步，我们将会得到更好的结果。
接下来介绍蛋白质折迭问题上的应用。
这个问题用最简单的说法，就是在我们知道一个蛋白质内胺
升级会员 