模拟退火法含程序(数模讲稿)

一、模拟退火法模拟退火法（参见［1,2］）作为一种适合于求解大规模的优化问题的技术，近来已引起极大的关注。特别是当优化问题有很多局部极值而全局极值又很难求出时，模拟退火法尤其有效。在实用上，它有效地“解决了”著名的旅行推梢员问题，即在必须依次访问每一个城市（共有N个城市）的前提下，为旅行推销员设计一条能够返回起点的最短旅程。模拟退火方法还被成功地用于设计复杂的集成电路，也就是说如何最佳地安排几十万个电路元件，使它们全部集成在一个很小的硅片上，而相互连接的线路之间的缠绕能够达到最小（参见［3,4］）.尽管模拟退火法的功效非凡，但它的算法实现却相对地简单，这一点似乎有些不可思议。请注意，我们上面提到的两个例子都属于组合极小化问题。现本类问题通常也有一个目标函数，但是函数的定义域并不是简单地由N个连续参变量组成的N维空间，而是一个离散的巨大空间，例如，由所有可能的城市旅行路线组成的集合，或者硅片电路元件的所有可能的分配方式的集合。构形空间中元素的数量相当巨大，根本不可能穷举，而且因为集合是离散的，我们也不可能“沿合适的方向连续下降”。因此在构形空间中，“方向”概念就没有什么意义了。后面我们还将介绍如何在其有连续控制参数的空间中利用模拟退火法。这种应用实际上要比组合问题复杂一些，因为其中又要出现“狭长山谷”的情况。正如在下文中我们将看到的，模拟退火法的试探步骤是“随机”的。但在一个狭窄且漫长的等高线山谷中，几乎所有的随机步骤都呈向上的趋势，因此，算法中需要增加一些技巧。模拟退火的核心思想与热力学的原理颇为相似，而且尤其类似于液体流动和结晶以及金属冷却和退火方式。在高温下，一种液体的大最分子彼此之间进行着相对自由的移动。如果该流体慢慢地冷却下来，热能可动性便会消失。大量原子常常能够自行排列成行，形成一个纯净的晶体，该晶体在各个方向上都被完全有序地排列在几百万倍于单个原子大小的距离之内。对于这个系统来说晶体状态是能量最低状态；而所有缓慢冷却的系统都可以自然达到这个能量最低状态，这可以说是一个令人惊奇的事实。实际上，如果某种液体金属被迅速冷却或被“猝熄”，那么它不会达到这一状态，而只能达到一种具有较高能量的多晶状态或非结晶状态。因此，这一过程的本质在于缓缓地致冷，以争取充足的时间，让大量原子在丧失可动性之前进行重新分布。这就是所谓退火在技术上的定义，同时也是确保达到低能量状态所必需的条件。尽管我们的比喻并不算贴切，但是迄今为止本身所讨论的所有极小化算法确实与快速冷却猝熄有某种关联之处。以往我们处理问题的方式都是:从初始点开始，立即沿下降方向前进，走得越远越好，似乎这样才能迅速求得问题的解但是，正如前面常常提到的，这种方法往往只能求得局部极小点，却求不到整体最小点。自然界本身的极小化算法则基于一种截然不同的方式，所谓的玻尔兹曼（Boltzmann）概率分布ProbE~ExpE/kT （1）表达了这样一种思想，即:一个处于热平衡状态且具有温度T的系统，其能量按照概率，分布于所有不同能量状态e之中。即使在很低的温度下，系统也有可能（虽然这种可能性很小）处于一个较高的能量状态。因此，相应地，系统也能够获得摆脱局部能量极小点的机会。并找到一个更好的、更接近于整体的极小点。式（1）中的参数k（称为玻尔兹曼常数）是一个自然常数，它的作用是将温度与能量联系起来。换句话说，在有些情况下系统的能量可上升，也可下降，但是温度越低，显著上升的可能性就越小。1953年，米特罗波利斯（Metropolis）及其合作者们首次将这种原理渗透到数值计算中。他们对一个模拟热力学系统提供了一系列选择项，并假设：系统构形从能量e变化到能量E的概率为p=ExpL（E-EkT]o很显然，如果E<E，1221’21p将大于1;在这类情况下，对构形的能量变化任意指定一个概率值p=1，也就是说，该系统总是取这个选择项。这种格式总是采取下降过程，偶尔采取上升步骤。目前已被公认为米特罗波利斯算法。为了将米特波利斯算法应用于热力学以外的系统，必须提供以下几项基本要素：对可能的系统构形的一种描述。一个有关构形内部随机变化的生成函数，这些变化将作为“选择项”提交给该系统。一个目标函数E（类似于能量），求解E的极小值，即为算法所要完成的工作。一个控制参数T（类似于温度）和一个退火进程，该进程用来说明系统是如何从高值向低值降低的，例如在温度T时每次下降步骤中要经过多少次随机的构形变化以及该步长是多大等等。应说明的是，这里“高”和“低”的含义，还有进程表的确定，都需要一定的物理知识和/或一些摸索的实验。模拟退火算法的模型模拟退火算法可以分解为解空间、目标函数和初始解三部分。模拟退火的基本思想:（1）初始化：初始温度T（充分大），初始解状态S（是算法迭代的起点），每个T值的迭代次数L对k=l,……,L做第(3)至第6步：产生新解9计算增量氐〜。©')-。©),其中C(S)为评价函数若MvO则接受S'作为新的当前解，否则以概率exp(-M/T)接受S，作为新的当前解.如果满足终止条件则输出当前解作为最优解,结束程序。终止条件通常取为连续若干个新解都没有被接受时终止算法。T逐渐减少，且T->0,然后转第2步。模拟退火算法新解的产生和接受可分为如下四个步骤：第一步，是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。第二步，是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。第三步，是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则：若At'vO则接受S'作为新的当前解S,否则以概率exp(-At'/T)接受S'作为新的当前解So第四步,是当新解被确定接受时,用新解代替当前解,这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃

时，则在原当前解的基础上继续下一轮试验。模拟退火算法与初始值无关，算法求得的解与初始解状态S（是算法迭代的起点）无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率1收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。模拟退火算法的简单应用作为模拟退火算法应用，讨论货郎担问题（TravellingSalesmanProblem，简记为TSP）:设有n个城市，用数码l,...,n代表。城市i和城市j之间的距离为d（i，j）i,j=l,...,n.TSP问题是要找遍访每个域市恰好一次的一条回路，且其路径总长度为最短.。求解TSP的模拟退火算法模型可描述如下:解空间:解空间S是遍访每个城市恰好一次的所有路经， 是｛1，……n｝的所有循环排列的集合,S中的成员记为h,w,…w｝，w,w,…w是1,2,……,n12n12n的一个排列，表明从w城市出发，依次经过w,w,…w城市，再返回w城123n1市。初始解可选为（1,……,n）;目标函数:目标函数为访问所有城市的路径总长度 或称为代价函数;我们要求的最优路径为目标函数（代价函数）为最小值时对应的路径。新路径的产生:随机产生1新路径的产生:随机产生1和n之间的两相异数k和m，不妨假设kvm，则将原路径1,w,・・・w,w,・・・w,w•…w1 2 k k+1 m m+1, 」）。）。w,w,・・・w,w,・・・w,w•…w1 2 m k+1 k m+1,这种变换方法就是将k和m对应的两个城市在路径序列中交换位置称为2-opt映射。

根据上述描述，模拟退火算法求解TSP问题的流程框图如下图2模拟退火算法的流程框图也可以采用其他的变换方法，有些变换有独特的优越性，有时也将它们交替使用，得到一种更好方法。例如：随机产生1和n之间的两相异数k和m，不妨假设kvm，则将原路径Cw,w,…w,w,…w,w…w1 2 k k+1 m m+1, n变为新路径：Cw,w,…w,w,…w,w,w…w1 2 m m-1 k+1 k, k+2 n上述变换方法就是将k和m之间对应的多个城市在路径序列颠倒位置,可简单说成是“逆转中间或者逆转两端”。根据上述分析，可写出用模拟退火算法求解TSP问题的伪程序：ProcedureTSPSA:begininit-of-T; //T为初始温度S={1，……，n}; 〃S为初始值termination=false;whiletermination=falsebeginfori=1toLdobegingenerate(SformS); 〃从当前回路S产生新回路S，△t:=f(S'))-f(S); 〃f(S)为路径总长IF(AtvO)OR(EXP(-At/T)>Random-of-[0,1])S=SZ;IFthe-halt-condition-is-TRUETHENtermination=true;End;T_lower;End;End模拟退火算法的应用很广泛，可以较高的效率求解最大截问题(MaxCutProblem)、0-1背包问题(ZeroOneKnapsackProblem)、图着色问题(GraphColouringProblem)、调度问题(SchedulingProblem)等等。模拟退火算法的参数控制问题模拟退火算法的应用很广泛，可以求解NP完全问题，但其参数难以控制，其主要问题有以下三点：温度T的初始值设置问题。温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。实际应用过程中，初始温度一般需要依据实验结果进行若干次调整。退火速度问题。模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索(退火)是相当必要的，但这需要计算时间。实际应用中，要针对具体问题的性质和特征设置合理的退火平衡条件。温度管理问题。温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中，由于必须考虑计算复杂度的切实可行性等问题，常采用如下所示的降温方式：T(t+1)=kxT(t)式中k为正的略小于1.00的常数，t为降温的次数二、有记忆的模拟退火法传统的模拟退火算法思路清晰、原理简单、可用于解决许多实际问题。但存在很多弊病，国内外学者并对其作相应的改进，得到了改进的模拟退火算法，包括：加温退火法、有记忆的模拟退火算法、带返回搜索的模拟退火算法、多次寻优法、回火退火算法等。针对特征选择的特点，本文选用了有记忆的模拟退火算法作为搜索策略。在传统的模拟退火过程中，算法终止于一个预先规定的停止准则S,如：控制参数t的值小于某个充分小的正数；相继的若干个Markov链中解未得到任何改善；两个相继Markov链所得解的差的绝对值小于某个充分小的正数等。但由于模拟退火算法的搜索过程是随机的，且当t值较大时可以接受部分恶化解，而随着t值的衰减，恶化解被接受的概率逐渐减小直至为0。另一方面，某些当前解要达到最优解时必须经过暂时恶化的“山脊”，因此，上述这些停止准则无法保证最终解正好是整个搜索过程中曾经达到的最优解。因此，在传统的算法中增加一个记忆器，使之能够记住搜索过程中遇到过的最好解，当退火结束时，将所得最终解与记忆器中的解比较并取较优者作为最后结果。三、组合极小化：旅行推销员问题

下面是我们用“旅行推销员问题”为具体实例说明模拟退火法的应用。假

设一个推销员，要去N个分别位于(x,y)的城市进行推销，并于最后返回他原来

ii所在的城市，要求每个城市只能去一次，而且所经过的路径要尽可能地短。这个间题属于一类所谓“NP-完全问题”这类问题求出一个精确解所需的计算时间是随N的增加以指数exp(常数XN)增长的。当N不断增大时，运行时间将迅速增加，进而导致费用高到令人难以接受的程度。旅行推销员问题也是众多极小化问题中的一种，它的目标函数e具有多个局部极小值。在实际应用当中，常常有足够多的条件可以从多个极小值中选出一个最小的，这个最小值即使不是绝对最小，也相当接近绝对最小了。退火法的目的就是要获得这个最小值，同时又要将计算量限制在N的低阶次的数量级上。旅行推销员问题也是按模拟退火问题的方式进行处理的。具体如下：1•构形将N个城市分别标记为i=1,2,N中的数，其中每个城市具有坐标一个构形就是数字1N的一个排列，可以解释为推销员途径的城市的顺序。2•调整林(Lin)曾经提出过一种所谓“转移有效集”这里的“转移”包括两种类型：(a)移走路径的某一段，然后对这段路径上的城市用相反次序重新进行排列，并用后者来代替前者。(b)移走某段路径，并用位于城市间的随机选取的另一段路径来取代被移走的路径。3.目标函数在旅行推销员问题的一种最简单的形式中，目标函数e被定义为旅途的总长度E二L三迓* —》+(y_y—F (10.9.2)Xii+1 ii+1i=1这里认为第N+1个点与第1个点是重合的。但是为了表明模拟退火法的灵活性，我们还要用到下面的技巧：假设推销员无端端地害怕飞越密西西比河，在这种情况下，我们对每个城市给定一个参数卩，如果该城市位于密西西比河i以东，卩取1；若在密西西比河以西则取-1。对于目标函数E，我们将其改写为:iE=为打 _J2+(y_y_X+九Q-卩》 (10.9.3)L'ii+1 ii+1 ii+1_i=1由于每过一次河都将以4九作为惩罚，因而现在我们设计的算法的目标，就变成了寻找尽可能回避过河的最短路径。路径长度对过河次数的相对重要性将由我们选择的九来确定。图10.9.1表明了所得的结果。显然，这种技巧可以推广到包含许多相互冲突的目的要求的极小化问题当中。

0 .5(C)图（A）表示的是从四个随机分布的城市中间找到的一条（接近）最短路径，中间竖虚线标识的是一条河流，但这是对过河没有附加您罚项的情况。在图（B）中对过河施加的惩罚项很大，而图中所示的解本身的过河次数也相应地只有少得不能再少的两次。在图（C）中惩罚项为负，这就是说，推销员实际上成了恣意偷渡的走私者！图1用模拟退火法解决旅行推销员问题4.退火进程这一步需要借助试验来确定。首先要进行一些随机调整，然后利用它们来确定从转移到转移过程中将会遇到的AE值之范围。对参数T取一初始值（这个初始值要远远大于通常所能遇到的AE的最大值），并以倍增的步长下减，每次使T总共减少10%。我们拿每个新的常数T值去试各种100N重构形,或1ON成功的重构形，无论哪个在前出现就取哪个。当实在不能再进一步减小E时，则停止。F面的旅行推销员程序利用了米特罗波利斯(Metropolis)算法，并展示了模拟退火技术应用于组合问题的几个主要方面。#include<stdio.h>#include<math.h>#defineTFACTR0.9 //退火进程：每步中t的下降值由该因子决定#defineALEN(a,b,c,d)aqrt(((b)-(a))*((b)-(a))+((d)-(c))*((d)-(c)))voidanneal(floatx[],floaty[],intiorder[],intncity)/*本算法用于求解在ncity个城市之间作往返旅行的最短路径，其中这ncity个城市的位置坐标存贮在数组x[l..ncity]和y[l..ncity]中。数组iorder[l..ncity]表示途径城市的顺序。在输出项中，iorder中的元素将被置为数字1到ncity的某排对，本程序将返回它所能求出的最佳选择路径。*/{intirbit1(unsignedlong*iseed);intmetrop(floatde,floatt);floatran3(long*idum);floatrevcst(floatx[],floaty[],intiorder[],intncity,intn[]);voidreverse(intiorder[],intncity,intn[]);floattrncst(floatx[],floaty[],intiorder[],intncity,intn[]);voidtrnspt(intiorder[],intncity,intn[]);intans,nover,nlimit,il,i2;inti,j,k,nsucc,nn,idec;staticintn[7];longidum;unsignedlongiseed;floatpath,de,t;nover=100*ncity; //在任何温度下，所试路径的最大次数nlimit=10*ncity; //在继续进行之前，成功的路径改变的最大次数path=0.0;t=0.5;for(i=1;i<ncity;i++){ //计算初始路径的长度i1=iorder[i];i2=iorder[i+1]:'path+=ALEN(x[il],x[i2],y[i1],y[i2]);

//将路径头尾相连并结束循环i1=iorder[ncity]; ////将路径头尾相连并结束循环i2=iorder[1]path+=ALEN(x[i1],x[i2],y[i1],y[i2]);idum=-1;//试验//试验100个温度值for(j=1;j<=100;j++){nsucc=0;for(k=1;k<=nover;k++){do{n[1]=1+(int)(ncity*ran3(&idum));//选择段的起始点„n[2]=1+(int)((ncity-1)*ran3(&idum));//„段的结尾if(n[2]>=n[1])++n[2];nn=1+((n[1]-n[2]+ncity-1)%ncity);//nn为不位于当前段上的城市数}while(nn<3);idec=irbit1(&iseed);//确定是否做段反转或段输送if(idec==0){ //做输送n[3]=n[2]+(int)(abs(nn-2)*ran3(&idum))+1;n[3]=1+((n[3]-1)%ncity);//输送到一个不在当前路径上的某处de=trncst(x,y,iorder,ncity,n); //计算代价ans=metrop(de,t); //做预测if(ans){++nsucc;path+=de;//输送工作结束//做段反转////输送工作结束//做段反转//计算代价//作预测//完成段反转}else{de=revcst(x,y,iorder,ncity,n);ans=metrop(de,t);if(ans){++nsucc;path+=de;reverse(iorder,ncity,n);if(nsucc>=nlimit)break; //如果成功的转移超过次数则提前结束}printf("\n%s%10.6f%s&12.6f\n","T=",t,"PathLength=",path);printf("SuccessfulMoves:%6d\n"，nsucc);t*=TFACTR; //退火进程if(nsucc==0)return; //如果不成功则返回}}#include<math.h>#defineALEN(a,b,c,d)sqrt(((b)-(a))*((b)-(a))+((d)-(c))*((d)-(c))))floatrevcst(floatx[]，floaty[]，intiorder[]，intncity，intn[])/*该子程序返回的是反转某给定路径所需的代价函数值。参数中ncity为城市数；数组x[l..ncity],y[l..ncity]为这些城市的位置坐标；iorder[n..ncity]为当前路线；数组n的头两个元素n[1]和n[2]分别代表将要被反转的路径段的起点城市和终点城市。子程序的输出项de为反转所需的代价值,但真正的反转过程并不是由该程序执行。*/{floatxx[5],yy[5],de;inti,ii;n[3]=1+((n[1]+ncity-2)%ncity); 〃找出n[1]以前的城市n[4]=1+(n[2]%ncity); //..找出n[2]以后的城市for(j=1;j<=4;j++){ii=iorder[n[j]]; //求四个所涉及到的城市的坐标xx[j]=x[ii];yy[j]=y[ii];}de=-ALEN(xx[1],xx[3],yy[1],yy[3]);//计算断开路径段两端所需的代价de-=ALEN(xx[2],xx[4],yy[2],yy[4]);//以及用相反顺序重新连接所需的代价de+=ALEN(xx[1],xx[4],yy[1],yy[4]);de+=ALEN(xx[2],xx[3],yy[2],yy[3]);returnde;}voidreverse(intiorder[],intncity,intn[])/*该子程序的作用是执行一段路径的反转过程。输入参数iorder[l..ncity]给出当前的路线顺序；向量n的前四个元素中，n⑴和n[2]分别为将要被反转的路径的起点和终点城市，n⑶和n[4]则分别为紧随n[1]之前和紧随n[2]之后的两个城市的标号，其中n[3]和n[4]由函数revest给出。在输出端，iorder又将被作为返回值，它的定义是n[1]到n[2]路段被反转后的行程路线。*/{intnn,j,k,l,itmp;nn=(1+((n[2]-n[1]+ncity)%ncity))/2;//为实现反转操作，必须交换这么多城市for(j=1;j<=nn;j++){k=1+((n[1]+j-2)%ncity);//从段的端部开始，顺序交换城市对，l=1+((n[2]-j+ncity)%ncity); //直至到达路径的中心itmp=iorder[k];iorder[k]=iorder[l];iorder[l]=itmp;}}#include<math.h>#defineALEN(a,b,c,d)sqrt(((b)-(a))*((b)-(a))+((d)-(c))*((d)-(c)))floattrncat(floatx[],floaty[],intiorder[],intncity,intn[])/*该子程序返回的是输送某段给定路径所需的代价函数值。输入参数中ncity为城市的个数；x[l..ncity]和y[l..ncity]为这些城市的位置坐标，数组n的前三个元索分别为：将要被输送的路径段的起、止点城市以及这段路径将要被插入处的标志城市(插在该城市之后)，该子程序的输出项de为计算出来的输送代价值，但实际的输送操作井不由该子程序执行。*/{floatxx[7],yy[7],de;intj,ii;n[4]=1+(n[3]%ncity); 〃找出位于n[3]之后的城市..n[5]=1+((n[1]+ncity-2)%ncity);//..位于n[1]之前的一个城市..n[6]=1+(n[2]%ncity); //..位于n[2]之后的一个城市..for(j=1;j<=6;j++){ii=iorder[n[j]]; //求出六个城市的有关坐标xx[j]=x[ii];yy[j]=y[ii];}de=-ALEN(xx[2],xx[6],yy[2],yy[6]); 〃计算下列操作所需代价，断开n[1]到n[2]de-=ALEN(xx[1],xx[5],yy[l],yy[5]); 〃间的路径。打开n[3]和n[4]之间的空间;de-=ALEN(xx[3],xx[4],yy[3],yy[4]); //连接这个空间中的路径段;以及连接n[5]、n[6]de+=ALEN(xx[1],xx[3],yy[1],yy[3]);de+=ALEN(xx[2],xx[4],yy[2],yy[4]);de+=ALEN(xx[5],xx[6],yy[5],yy[6]);returnde;}#include"nrutil.h"voidtrnspt(intiorder[],intncity,intn[])/*该子程序的作用是执行真正的段输送操作，输入参数iorder[l..ncity]给出当前的路径顺序；数组n共有6个元素，其意义分别为：n⑴和n[2]分别代表将要被输送的路径段的起点城市和终点城市;n[3]和n[4]为两个相邻的城市,n[1]至n[2]路径段即将放入它们中间;n⑸和n⑹分别为n[1]之前和n[2]之后的两个城市。在这六个元索中n[4],n⑸和n[6]由函数trncst给出。在输出端，iorder将根据路径段的移动和变化作出相应的修改。*/{intml,m2,m3,nn,j,jj,*jorder;jorder=ivector(1,ncity);m1=1+((n[2]-n[1]+ncity)%ncity);//找出位于n[1]和n[2]间城市个数m2=1+((n[5]-n[4]+ncity)%ncity);〃„n[4]到n[5]间的城市个数m3=1+((n[3]-n[6]+ncity)%ncity);//•„n⑹和n[3]间的城市个数nn=l;for(j=1;j<=m1;j++){jj=1+((j+n[1]-2)%ncity); //复制所选路径段jorder[nn++]=iorder[jj];}if(m2>0){for(j=1;jv=m2;j++){ 〃复制n[4倒n⑸间的路径段jj=1+((j+n[4]-2)%ncity);jorder[nn++]=iorder[jj];}if(m3>0){for(j=1;jv=m3;j++){ 〃最后，复制n[6]到n[3]间的路径段jj=1+((j+n[6]-2)%ncity);jorder[nn++]=iorder[jj];}}for(j=1；j<=ncity；j++) 〃将jorder拷贝回iorderiorder[j]=jorder[j];free_ivector(jorder;1,ncity);}#includevmath.h>intmetrop(float,de,floatt)/*该子程序为米特罗波利斯算法的程序实现。metrop返回的是一个布尔型变量，由该变量决定是否接受一个使目标函数e产生改变量de的重构形。如果de<0则metrop=l(真)，而当de>O时，metrop为真的概率是exp(—de/t)，这里才是一个由退火进程决定的温度值。*/{floatran3(long*idum);staticlonggljdum=1;returndev0.0||ran3(&gljdum)vexp(-de/t);}模拟退火法的基本思想也可以应用于具有连续N维控制参数的极小化问题当中，例如，在某个函数f6）（这里x为一个N维向量）有许多局部极小点的情况下，求解它的（理想的全局的）极小值。这时米特罗波利斯算法所需的四要素可以具体化为：1） 目标函数即为f（x）的值；2） 系统状态描态即为N维空间中的点x；3）类似于温度的控制参数T以及一个使T逐渐降低的退火进程仍为原先的定义；4）描述构形内部随机变化的发生器即为一个从x到x+Ax采取随机步骤的方法。在上述四要素中问题最大的是最后一条。目前已发表的文献中［7-10］，介绍了几种不同的选择办法。但我们认为，这些方法都不算成功。问题在于“效率”二字：当局部的向下运动存在时，如果某个随机变化发生器几乎总是作出向上运动的决策，那么它的效率就很低。我们认为，一个好的发生器在等高线的“窄谷”中仍应保持高效性；当算法在接近收敛到极小点处时，它的效率也不应越变越低。在上面我们提到的几篇文献中，除［7］中介绍的方法之外，其他所有方法都表现不同程度的低效性。下面我们将要介绍的这种方法，利用了下降单纯形法的一种修改后的形式。首先我们将米特罗波利斯算法四要素中的系统状态描述，由单个点x改为一个其有N+1个点的单纯形；单纯形的“操作”分为反射、扩张和收缩三种。然后我们将一个正的、呈对数分布的随机变量（与温度T成比例）添加到存贮的函数值中（该函数值与单纯形的每个顶点都有关联），再从每个被当作替代的新点的函数值中减去一个类似的随机变量。和普通的米特罗波利斯方法一样，这种方法总是接受一个真正的下降步骤。但有时也接受一次上升步骤。在极限过程T-0中，该算法恰好变成了下降的单纯形法，并收敛到一个局部极小点。在T的某有限值处，单纯形将扩展到一定的规模，其大小接近于在这个温度值所能达到的区域；然后单纯形在这个区域内部做随机的滚动翻转式布朗运动，并在该过程中抽取一些新的、近似随机的点作为样本。一个区域被利用的效率与其狭窄度（对椭球状山谷，指它的主轴比率）及方向性均无关。如果温度降低得足够缓慢，那么单纯形将极有可能收缩到那个区域内，而那个区域内包含已遇到的最低相对极小。由于在所有模拟退火法的应用场合中，“足够缓慢”一词的含意根据问题的不同可以有相当大的细微区别，因而成功与否往往取决于退火进程选择。下面几种可能性我们认为值一试：•每经过m步移动之后，将T减到（1_e）T,这里。*/m的具体值要通过实验确定。•设总的移动步数为K,每经过m步移动之后将T减到t二TG-kT、，其0'中k是到目前为止所经过的步数的累加值。d为一常数，可取为1，2:4等。（X的最佳值取决于各种深度的相对于极小的统计分布，稍大一些的X值在较低温度时，需化费的迭代次数将更多；•每经过m步骤移动之后，置T为0乘f-f，其中0为一个阶数为1的常1b数，其具体值由试验确定；f为目前单纯形中最小的函数值；f为曾经遇到的最1b佳函数值。但应注意的是，T的降低幅度一次不要超过某个分数值丫。另一个策略方法上的问题是，当单纯形的某个顶点被放弃并让位于“永远的最佳点”时，是否需要采取重新开始的步骤（但我们必须保证在进行这项工作时，这个“永远的最佳点”当前并不在单纯形中!）。对于有些问题重新开始（例如，只要温度降低了因子3即执行重新开始步骤）的效果极佳，而对于另外一些问题，重新开始不仅没有任何效果反而会产生负作用。上述两种截然相反的情况，我们都找到了例子可以作为例证。将下面的程序amebsa采用改进下降单纯形法。#include<math.h>#include"nrutil.h"#defineGET_PSUM\for(n=1;n<=ndim;n++){\for(sum=0.0,m=1l;m<=mpts;m++)sum+=p[m][n];\psum[n]=sum;}externlongidum; //由主程序定义和初始值floattt; //与amotsa进行通信voidamebsa(float**p,floaty[],intndim,floatpb[],float*yb,floatftol,float(*funk)(float[]),int*iter,floattemptr)/*用摸拟退火与Nelder、Mead的下降单纯形法相结合的方法求多元函数funk(x)的极小值。其中x[1..ndim]为一个ndim维向量。作为输入参数的矩阵p[1..ndim][1..ndim〕共有ndim+1行，每行均为一个ndim维向量。分别代表初始单纯形的各个顶点。amebsa的输入项还包括矢量y[1..ndim+1]、浮点效ftol以及iter和temptr。其中y的各个元素将被初始化为函数funk在P的ndim+1个顶点(即行)的值；ftol为函数值所要达到的相对收敛容许限，一旦满足要求，程序将尽早返回;iter和temptr的意义分别是函数值的计算次数和温度。程序在某退火温度temptr处进行iter次函数值计算后返回。而接下来要做的事就是根据退火进程调低温度temptr；重置iter并再次调用该程序(每次调用时其他参数保持不变)。如果iter以正值返回，则说明程序正常收敛。如果第一次调用时yb被初始化为一个很大的数，则yb和pb[1..ndim]将依次返回已遇到的最佳函数值和最佳点(这个最佳点可以水远不是单纯形中的点)。*/{floatamotsa(float**p,floaty[],floatpsum[],intndim,floatpb[],float*yb,float(*funk)(float[]),intihi,float*yhi,floatfac);floatran1(long*idum);inti,ihi,ilo,j,m,n,mpts=nidm+1;floatrtol,sum,swap,yhi,ylo,ynhi,ysave,yt,ytry,*psum;psum=vector(1,ndim);tt=-temptr;GET_PSUMfor(;;){ilo=1; //确定最高点(即差点)、次高点和最低点(即最佳点)ihi=2;ynhi=ylo=y[1]+tt*log(ran1(&idum));//我们所“看到”的顶点总是处于//随机的热起伏状态yhi=y[2]+tt*log(ran1(&idum));if(ylo>yhi){ihi=1;ilo=2;ynhi=yhi;yhi=ylo;ylo=ynhi;}for(i=3;i<=mpts;i++){ //对单纯形中的点进行循环yt=y[i]+tt*log(ran1(&idum)); //更多的热起伏运动if(yt<=ylo){ilo=i;ylo=yt;}if(yt>yhi){ynhi=yhi;ihi=i;yhi=yt;}elseif(yt>ynhi){ynhi=yt;}}rtol=2.0*fabs(yhi-ylo)/(fabs(yhi)+fabs(ylo));//计算从最高点到最低点的范//围，若合乎要求，则返回if(rtol<ftol||*iter<0){ //若返回，将最佳点和最佳位放入槽1中swap=y[1];y[1]=y[ilo];y[ilo]=swap;for(n=1;n<=ndim;n++){swap=p[1][n];p[1][n]=p[ilo][n];p[ilo][n]=swap;}break;}*iter-=2; //开始新的一轮迭代，首先从高点通过单纯形的表面，//以因子-1做外插，即从离点对单纯形进行反射。ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi.-1.0);

if(ytry<=ylo){//得到了比最佳点还要好的结果，因此用因子2再做一次外推ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi,2.0);}elseif(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;j<ndim;j++){psum[j]=0.5*(p[i][j]+p[ilo][j]);p[i][j]=psum[j];}y[i]=(*funk)(psum);}}//重新计算psum//纠正计数器//重新计算psum//纠正计数器//在主程序中定义和初始赋值//在amobsa中定义}}else++(*iter);}free_vector(psum,1,ndim);}#include<math.h>#include"nrutil.h"externlongidum;externfloattt;floatamotsa(float**p，floaty[],floatpsum[]，intndim，floatpb[]，float*yb，float(*funk)(float[]),intihi，float*yhi,floatfac)//从高点通过单纯形的表面，以因子fac作外推，如果得到的新点较好，则用它取代高点。{floatran1(long*idum);intj;floatfac1,fac2,yflu,ytry,*ptry;ptry=vector(1,ndim);'facl=(1.0-fac)/ndim;fac2=fac1-fac;for(j=1;j<=ndim;j++)ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;ytry=(*funk)(ptry);if(ytry<=*yb){ //存储至令最好的for(j=1;j<=ndim;j++)pb[j]=ptry[j];*yb=ytry;}yflu=ytry-tt*log(ran1(&idum));//我们曾经对所有的当前顶点添加热起伏运动，if(yflu<*yhi){