模拟退火法含程序.doc

10.9 模拟退火法 模拟退火法(参见1,2)作为一种适合于求解大规模的优化问题的技术，近来已引起极大的关注。特别是当优化问题有很多局部极值而全局极值又很难求出时，模拟退火法尤其有效。在实用上，它有效地“解决了”著名的旅行推梢员问题，即在必须依次访问每一个城市(共有N个城市)的前提下，为旅行推销员设计一条能够返回起点的最短旅程。模拟退火方法还被成功地用于设计复杂的集成电路，也就是说如何最佳地安排几十万个电路元件，使它们全部集成在一个很小的硅片上，而相互连接的线路之间的缠绕能够达到最小(参见3,4).尽管模拟退火法的功效非凡，但它的算法实现却相对地简单，这一点似乎有些不可思议。 请注意，我们上面提到的两个例子都属于组合极小化问题。现本类问题通常也有一个目标函数，但是函数的定义域并不是简单地由N个连续参变量组成的N维空间，而是一个离散的巨大空间，例如，由所有可能的城市旅行路线组成的集合，或者硅片电路元件的所有可能的分配方式的集合。构形空间中元素的数量相当巨大，根本不可能穷举，而且因为集合是离散的，我们也不可能“沿合适的方向连续下降”。因此在构形空间中，“方向”概念就没有什么意义了。 后面我们还将介绍如何在其有连续控制参数的空间中利用模拟退火法。这种应用实际上要比组合问题复杂一些，因为其中又要出现“狭长山谷”的情况。正如在下文中我们将看到的，模拟退火法的试探步骤是“随机”的。但在一个狭窄且漫长的等高线山谷中，几乎所有的随机步骤都呈向上的趋势，因此，算法中需要增加一些技巧。 模拟退火的核心思想与热力学的原理颇为相似，而且尤其类似于液体流动和结晶以及金属冷却和退火方式。在高温下，一种液体的大最分子彼此之间进行着相对自由的移动。如果该流体慢慢地冷却下来，热能可动性便会消失。大量原子常常能够自行排列成行，形成一个纯净的晶体，该晶体在各个方向上都被完全有序地排列在几百万倍于单个原子大小的距离之内。对于这个系统来说，晶体状态是能量最低状态；而所有缓慢冷却的系统都可以自然达到这个能量最低状态，这可以说是一个令人惊奇的事实。实际上，如果某种液体金属被迅速冷却或被“猝熄”，那么它不会达到这一状态，而只能达到一种具有较高能量的多晶状态或非结晶状态。 因此，这一过程的本质在于缓缓地致冷，以争取充足的时间，让大量原子在丧失可动性之前进行重新分布。这就是所谓退火在技术上的定义，同时也是确保达到低能量状态所必需的条件。尽管我们的比喻并不算贴切，但是迄今为止本身所讨论的所有极小化算法，确实与快速冷却猝熄有某种关联之处。以往我们处理问题的方式都是:从初始点开始，立即沿下降方向前进，走得越远越好，似乎这样才能迅速求得问题的解。但是，正如前面常常提到的，这种方法往往只能求得局部极小点，却求不到整体最小点。自然界本身的极小化算法则基于一种截然不同的方式，所谓的玻尔兹曼(Boltzmann)概率分布 （10.9.1）表达了这样一种思想，即:一个处于热平衡状态且具有温度T的系统，其能量按照概率，分布于所有不同能量状态之中。即使在很低的温度下，系统也有可能（虽然这种可能性很小）处于一个较高的能量状态。因此，相应地，系统也能够获得摆脱局部能量极小点的机会。并找到一个更好的、更接近于整体的极小点。式(10.9.1)中的参数k(称为玻尔兹曼常数)是一个自然常数，它的作用是将温度与能量联系起来。换句话说，在有些情况下系统的能量可上升，也可下降，但是温度越低，显著上升的可能性就越小。 1953年，米特罗波利斯(Metropolis)及其合作者们首次将这种原理渗透到数值计算中。他们对一个模拟热力学系统提供了一系列选择项，并假设：系统构形从能量变化到能量的概率为。很显然，如果，将大于1；在这类情况下，对构形的能量变化任意指定一个概率值，也就是说，该系统总是取这个选择项。这种格式总是采取下降过程，偶尔采取上升步骤。目前已被公认为米特罗波利斯算法。 为了将米特波利斯算法应用于热力学以外的系统，必须提供以下几项基本要素：1. 对可能的系统构形的一种描述。2. 一个有关构形内部随机变化的生成函数，这些变化将作为“选择项”提交给该系统。3. 一个目标函数 (类似于能量)，求解的极小值，即为算法所要完成的工作。4. 一个控制参数T(类似于温度)和一个退火进程，该进程用来说明系统是如何从高值向低值降低的，例如在温度T时每次下降步骤中要经过多少次随机的构形变化以及该步长是多大等等。应说明的是，这里“高”和“低”的含义，还有进程表的确定，都需要一定的物理知识和/或一些摸索的实验。10.9.1组合极小化：旅行推销员问题下面是我们用“旅行推销员问题”为具体实例说明模拟退火法的应用。假设一个推销员，要去N个分别位于的城市进行推销，并于最后返回他原来所在的城市，要求每个城市只能去一次，而且所经过的路径要尽可能地短。这个问题属于一类所谓“NP-完全问题”。这类问题求出一个精确解所需的计算时间是随N的增加以指数exp(常数N)增长的。当N不断增大时，运行时间将迅速增加，进而导致费用高到令人难以接受的程度。旅行推销员问题也是众多极小化问题中的一种，它的目标函数具有多个局部极小值。在实际应用当中，常常有足够多的条件可以从多个极小值中选出一个最小的，这个最小值即使不是绝对最小，也相当接近绝对最小了。退火法的目的就是要获得这个最小值，同时又要将计算量限制在N的低阶次的数量级上。旅行推销员问题也是按模拟退火问题的方式进行处理的。具体如下：1.构形 将N个城市分别标记为中的数，其中每个城市具有坐标。 一个构形就是数字的一个排列，可以解释为推销员途径的城市的顺序。2.调整 林(Lin)曾经提出过一种所谓“转移有效集”，这里的“转移”包括两种类型：(a)移走路径的某一段，然后对这段路径上的城市用相反次序重新进行排列，并用后者来代替前者。(b)移走某段路径，并用位于城市间的随机选取的另一段路径来取代被移走的路径。3.目标函数 在旅行推销员问题的一种最简单的形式中，目标函数被定义为旅途的总长度 （10.9.2）这里认为第N+1个点与第1个点是重合的。但是为了表明模拟退火法的灵活性，我们还要用到下面的技巧：假设推销员无端端地害怕飞越密西西比河，在这种情况下，我们对每个城市给定一个参数，如果该城市位于密西西比河以东，取1；若在密西西比河以西则取-1。对于目标函数，我们将其改写为： （10.9.3）由于每过一次河都将以作为惩罚，因而现在我们设计的算法的目标，就变成了寻找尽可能回避过河的最短路径。路径长度对过河次数的相对重要性将由我们选择的来确定。图10. 9. 1表明了所得的结果。显然，这种技巧可以推广到包含许多相互冲突的目的要求的极小化问题当中。 (A) (B)(C)图(A)表示的是从四个随机分布的城市中间找到的一条(接近)最短路径，中间竖虚线标识的是一条河流，但这是对过河没有附加您罚项的情况。在图(B)中对过河施加的惩罚项很大，而图中所示的解本身的过河次数也相应地只有少得不能再少的两次。在图(C)中惩罚项为负，这就是说，推销员实际上成了恣意偷渡的走私者！图10.9. 1用模拟退火法解决旅行推销员问题 4.退火进程 这一步需要借助试验来确定。首先要进行一些随机调整，然后利用它们来确定从转移到转移过程中将会遇到的值之范围。对参数T取一初始值(这个初始值要远远大于通常所能遇到的的最大值)，并以倍增的步长下减，每次使T总共减少10%。我们拿每个新的常数T值去试各种100N重构形，或10N成功的重构形，无论哪个在前出现就取哪个。当实在不能再进一步减小时，则停止。下面的旅行推销员程序利用了米特罗波利斯(Metropolis)算法，并展示了模拟退火技术应用于组合问题的几个主要方面。# include # include # define TFACTR 0.9 /退火进程：每步中t的下降值由该因子决定# define ALEN(a,b,c,d) aqrt (b)-(a)*(b)-(a)+(d)-(c)*(d)-(c)void anneal(float x, float y, int iorder, int ncity)/* 本算法用于求解在ncity个城市之间作往返旅行的最短路径，其中这ncity个城市的位置坐标存贮在数组xl.ncity和y1.ncity中。数组iorder1.ncity表示途径城市的顺序。在输出项中，iorder中的元素将被置为数字l到ncity的某排对，本程序将返回它所能求出的最佳选择路径。*/int irbit1(unsigned long *iseed);int metrop(float de, float t);float ran3(long *idum);float revcst(float x,float y,int iorder,int ncity,int n);void reverse(int iorder,int ncity,int n);float trncst(float x,float y,int iorder,int ncity,int n);void trnspt(int iorder,int ncity, int n);int ans,nover,nlimit,il,i2;int i,j,k,nsucc,nn,idec;static int n7;long idum;unsigned long iseed;float path,de,t;nover=100*ncity; /在任何温度下，所试路径的最大次数nlimit=10*ncity; /在继续进行之前，成功的路径改变的最大次数path=0.0;t=0.5;for(i=1;incity;i+) /计算初始路径的长度 i1=iorderi; i2=iorderi+1: path+=ALEN(xil,xi2,yi1,yi2);i1=iorderncity; /将路径头尾相连并结束循环i2=iorder1path+=ALEN(xi1,xi2,yi1,yi2);idum=-1;iseed=111for (j=1;j=100;j+) /试验100个温度值 nsucc=0; for (k=1;k=n1) +n2; nn=1+(n1-n2+ncity-1) % ncity); /nn为不位于当前段上的城市数 while (nn=nlimit) break; /如果成功的转移超过次数则提前结束printf(n %s %10.6f %s &12.6f n,T=,t, Path Length=,path);printf(Successful Moves: %6dn，nsucc);t*=TFACTR; /退火进程if (nsucc=0) return; /如果不成功则返回#include #define ALEN(a,b,c,d) sqrt(b)-(a)*(b)-(a)+(d)-(c)*(d)-(c)float revcst(float x，float y，intiorder，int ncity，int n)/*该子程序返回的是反转某给定路径所需的代价函数值。参数中ncity为城市数；数组x1.ncity，y1.ncity为这些城市的位置坐标；iordern. ncity为当前路线；数组n的头两个元素n1和n2分别代表将要被反转的路径段的起点城市和终点城市。子程序的输出项de为反转所需的代价值，但真正的反转过程并不是由该程序执行。*/ float xx5,yy5,de; int i,ii; n3=1+ (n1+ncity-2) % ncity); /找出n1以前的城市 n4=1+(n2 % ncity); /.找出n2以后的城市 for (j=1;j=4;j+) ii=iordernj; /求四个所涉及到的城市的坐标 xxj=xii; yyj=yii; de=-ALEN(xx1,xx3,yy1,yy3); /计算断开路径段两端所需的代价 de-=ALEN(xx2,xx4,yy2,yy4); /以及用相反顺序重新连接所需的代价 de+=ALEN(xx1,xx4,yy1,yy4); de+=ALEN(xx2,xx3,yy2,yy3); return de;void reverse (int iorder，int ncity，int n)/*该子程序的作用是执行一段路径的反转过程。输入参数iorderl.ncity给出当前的路线顺序；向量n的前四个元素中，n1和n2分别为将要被反转的路径的起点和终点城市，n3和n4则分别为紧随n1之前和紧随n2之后的两个城市的标号，其中n3和n4由函数revcst给出。在输出端，iorder又将被作为返回值，它的定义是n1 到n2路段被反转后的行程路线。*/ int nn,j,k,l,itmp; nn= (1 + (n2-n1+ncity) % ncity) )/2; /为实现反转操作，必须交换这么多城市 for (j=1;j=nn;j+) k=1 + (n1+j-2) % ncity); /从段的端部开始，顺序交换城市对， l=1 + (n2-j+ncity) % ncity); /直至到达路径的中心 itmp=iorderk; iorderk =iorder l; iorderl=itmp; # include # define ALEN(a,b,c,d) sqrt(b)-(a)*(b)-(a)+(d)-(c)*(d)-(c)float trncat (float x, float y, int iorder, int ncity, int n)/*该子程序返回的是输送某段给定路径所需的代价函数值。输入参数中ncity为城市的个数；x1.ncity和y1.ncity 为这些城市的位置坐标，数组n的前三个元索分别为：将要被输送的路径段的起、止点城市以及这段路径将要被插入处的标志城市(插在该城市之后)，该子程序的输出项de为计算出来的输送代价值，但实际的输送操作井不由该子程序执行。*/ float xx7,yy7,de; int j,ii; n4=1+(n3 % ncity); /找出位于n3之后的城市. n5=1+(n1+ncity-2) % ncity); / .位于n1之前的一个城市. n6=1+(n2 % ncity); / .位于n2之后的一个城市. for (j=1;j=6;j+) ii = iordernj; /求出六个城市的有关坐标 xxj=xii; yyj=yii; de = -ALEN(xx2,xx6,yy2,yy6); /计算下列操作所需代价，断开n1到n2de -= ALEN(xx1,xx5,yyl,yy5); /间的路径。打开n3和n4之间的空间;de -= ALEN(xx3,xx4,yy3,yy4); /连接这个空间中的路径段;以及连接n5、n6de += ALEN(xx1,xx3,yy1,yy3);de += ALEN(xx2,xx4,yy2,yy4);de += ALEN(xx5,xx6,yy5,yy6);return de;# include nrutil.hvoid trnspt(int iorder, int ncity, int n)/*该子程序的作用是执行真正的段输送操作，输入参数iorder1.ncity给出当前的路径顺序；数组n共有6个元素，其意义分别为：n1和n2分别代表将要被输送的路径段的起点城市和终点城市；n3和n4为两个相邻的城市，n1至n2路径段即将放入它们中间;n5和n6分别为n1之前和n2之后的两个城市。在这六个元索中n 4,n5和n6由函数trncst给出。在输出端，iorder将根据路径段的移动和变化作出相应的修改。*/ int ml,m2,m3,nn,j,jj,*jorder; jorder=ivector(1,ncity) ; m1=1+(n2-n1+ncity) % ncity); /找出位于n1和n2间城市个数 m2=1+(n5-n4+ncity) % ncity); /n4到n5间的城市个数 m3=1+(n3-n6+ncity) % ncity); /n6和n3间的城市个数nn=1;for(j=1; j0) for(j=1;j0) for(j=1;j=m3;j+) /最后，复制n6到n3间的路径段 jj=1+(j+n6-2) % ncity); jordernn+=iorderjj; for (j=1;j=ncity;j+) /将jorder拷贝回iorder iorderj=jorderj;free_ivector(jorder;1,ncity);# include int metrop (float, de, float t )/*该子程序为米特罗波利斯算法的程序实现。metrop返回的是一个布尔型变量，由该变量决定是否接受一个使目标函数e产生改变量de的重构形。如果de0时，metrop为真的概率是exp(一de/t)，这里才是一个由退火进程决定的温度值。*/ float ran3(long *idum); static long gljdum =1 ; return de0.0 | ran3(&gljdum) exp(-de/t);10.9.2 模拟遇火法在连续极小化问题中的应用 模拟退火法的基本思想也可以应用于具有连续N维控制参数的极小化问题当中，例如，在某个函数(这里x为一个N维向量)有许多局部极小点的情况下，求解它的(理想的全局的)极小值。这时米特罗波利斯算法所需的四要素可以具体化为：1)目标函数即为的值；2)系统状态描态即为N维空间中的点x；3)类似于温度的控制参数T以及一个使T逐渐降低的退火进程仍为原先的定义；4)描述构形内部随机变化的发生器即为一个从x到x+x采取随机步骤的方法。在上述四要素中问题最大的是最后一条。目前已发表的文献中7-10，介绍了几种不同的选择x办法。但我们认为，这些方法都不算成功。问题在于“效率”二字：当局部的向下运动存在时，如果某个随机变化发生器几乎总是作出向上运动的决策，那么它的效率就很低。我们认为，一个好的发生器在等高线的“窄谷”中仍应保持高效性；当算法在接近收敛到极小点处时，它的效率也不应越变越低。在上面我们提到的几篇文献中，除7中介绍的方法之外，其他所有方法都表现不同程度的低效性。 下面我们将要介绍的这种方法，利用了下降单纯形法(见第10.4节)的一种修改后的形式。首先我们将米特罗波利斯算法四要素中的系统状态描述，由单个点x改为一个其有N+1个点的单纯形；单纯形的“操作”和第10.4节中介绍的相同，分为反射、扩张和收缩三种。然后我们将一个正的、呈对数分布的随机变量(与温度T成比例)添加到存贮的函数值中 (该函数值与单纯形的每个顶点都有关联)，再从每个被当作替代的新点的函数值中减去一个类似的随机变量。和普通的米特罗波利斯方法一样，这种方法总是接受一个真正的下降步骤。但有时也接受一次上升步骤。在极限过程T0中，该算法恰好变成了下降的单纯形法，并收敛到一个局部极小点。 在T的某有限值处，单纯形将扩展到一定的规模，其大小接近于在这个温度值所能达到的区域；然后单纯形在这个区域内部做随机的滚动翻转式布朗运动，并在该过程中抽取一些新的、近似随机的点作为样本。一个区域被利用的效率与其狭窄度(对椭球状山谷，指它的主轴比率)及方向性均无关。如果温度降低得足够缓慢，那么单纯形将极有可能收缩到那个区域内，而那个区域内包含已遇到的最低相对极小。 由于在所有模拟退火法的应用场合中，“足够缓慢”一词的含意根据问题的不同可以有相当大的细微区别，因而成功与否往往取决于退火进程选择。下面几种可能性我们认为值一试：每经过m步移动之后，将T减到(1)T，这里/m的具体值要通过实验确定。设总的移动步数为K，每经过m步移动之后将T减到 ，其中k是到目前为止所经过的步数的累加值。为一常数，可取为1，2: 4等。的最佳值取决于各种深度的相对于极小的统计分布，稍大一些的值在较低温度时，需化费的迭代次数将更多；每经过m步骤移动之后，置T为乘，其中为一个阶数为1的常数，其具体值由试验确定；为目前单纯形中最小的函数值；为曾经遇到的最佳函数值。但应注意的是，T的降低幅度一次不要超过某个分数值。另一个策略方法上的问题是，当单纯形的某个顶点被放弃并让位于“永远的最佳点”时，是否需要采取重新开始的步骤(但我们必须保证在进行这项工作时，这个“永远的最佳点”当前并不在单纯形中!)。对于有些问题重新开始(例如，只要温度降低了因子3即执行重新开始步骤)的效果极佳，而对于另外一些问题，重新开始不仅没有任何效果反而会产生负作用。上述两种截然相反的情况，我们都找到了例子可以作为例证。将下面的程序amebsa同第10.4节中与之相应的amoeba进行比较，你会发现，参数iter的使用方式在两个程序中略有不同。# include # include nrutil.h# define GET_PSUM for (n=1;n=ndim;n+) for (sum=0.0,m=1l;myhi) ihi=1; ilo=2; ynhi=yhi; yhi=ylo; ylo=ynhi; for (i=3;i=mpts;i+) /对单纯形中的点进行循环 yt=yi+tt * log (ran1(&idum); /更多的热起伏运动 if (yt yhi) ynhi=yhi; ihi=i; yhi=yt; else if (yt ynhi) ynhi=yt; rtol=2.0 * fabs(yhi-ylo)/(fabs(yhi)+fabs(ylo); /计算从最高点到最低点的范圈，若合乎要求，/则返回 if (rtolftol | *iter0) /若返回，将最佳点和最佳位放入槽1中 swap=y1; y1=yilo; yilo=swap; for (n=1;n=ndim;n+) swap=p1n; p1n=pilon; pilon=swap; break; *iter -= 2; /开始新的一轮迭代，首先从高点通过单纯形的表面，/以因子-1做外插，即从离点对单纯形进行反射。ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi. -1.0);if (ytry = ynhi ) /反射后的点不如次高点，因此需要找一个中间的较低点, /即做一次一维收缩 ysave=yhi; ytry=amotsa(p,y,psum,ndim,Pb,yb,funk, ihi,&yhi,0.5); if (ytry = ysave) /似乎无法摆脱高点，最好围绕最低点也即最佳点进行收缩 for (i=1;i=mpts;i+) if (i!=ilo) for (j=1;jndim;j+) psumj=0.5 * (pij+piloj); pij=psumj; yi= (* funk) (psum); *iter -= ndim;GET_PSUM /重新计算psum else +(*iter); /纠正计数器free_vector(psum,1,ndim);# include # include nrutil.hextern long idum; /在主程序中定义和初始赋值extern float tt; /在amobsa中定义float amotsa(float *p，float y,float psum，int ndim，float pb，float *yb，float(*funk)(float), int ihi，float *yhi, float fac)/从高点通过单纯形的表面，以因子fac作外推，如果得到的新点较好，则用它取代高点。 float ran1(long *idum); int j; float fac1,fac2,yflu,ytry, *ptry; ptry=vector(1, ndim); facl=(1.0 - fac)/ndim; fac2=fac1-fac; for (j=1;j=ndim;j+) ptryj =psumj * fac1 - pihij * fac2; ytry= (* funk) (ptry); if (ytry = *yb) /存储至令最好的 for(j=1;j=ndim;j+) pbj=ptryj; *yb = ytry; yflu=ytry-tt*log(ran1(&idum); /我们曾经对所有的当前顶点添加热起伏运动， if (yflu*yhi) /但这里我们要减少热起伏，目的是给单纯形一 yihi=ytry; /个热布朗运动，就象接受任何提议的变化一样 *yhi=yflu; for (j=1;j=ndim;j+) psumj += ptryj-pihij; pihij=ptryj; free_vector(ptry,1,ndim); return yflu; 模拟退火法在优化方法中将处于什么样的地位、扮演什么样的角色？对于这个问题目前还没有足够的实践经验来明确回答它。但是这种方法确有几个极为惹人注目的特点，与其他优化技术相比显得尤其突出： 首先，这种方法不象“急功近利的小人”那样表现“贪婪”，也就是说，它能迅速求出某些局部极小值，但那些实际上是不利的极小值并不会让它轻易上当。假如给定一些足够一般的重构形，这种方法会在深度小于T的局部极小值之间自由地游动。当T值降低，这些具备频繁访问条件的极小值的数目也将逐渐减少。 第二点，这种方法中有关构形的判定是趋向于以一种合乎逻辑的顺序进行的。当控制参数T很大时，那些引起最大能量差异的诸多变化将被筛去。T减小时，这些判定就变得更恒定，相应地，算法将注意力就更多地转向解的较精细的改进方面。例如，在旅行推销员问题中，如果密西比河迂回于旅途之中，而且值很大，那么，