模拟退火法含程序.docx
《模拟退火法含程序.docx》由会员分享，可在线阅读，更多相关《模拟退火法含程序.docx（23页珍藏版）》请在冰点文库上搜索。
模拟退火法含程序
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
#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[1..ncity]中。
数组iorder[1..ncity]表示途径城市的顺序。
在输出项中，iorder中的元素将被置为数字l到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;ii1=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;
iseed=111
for（j=1;j<=100;j++）{//试验100个温度值
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;
trnspt（iorder,ncity,n）;//输送工作结束
}
}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
#defineALEN（a,b,c,d）sqrt（（（b）-（a））*（（b）-（a））+（（d）-（c））*（（d）-（c））））
floatrevcst（floatx[]，floaty[]，intiorder[]，intncity，intn[]）
/*
该子程序返回的是反转某给定路径所需的代价函数值。
参数中ncity为城市数；数组x[1..ncity]，y[1..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[1]和n[2]分别为将要被反转的路径的起点和终点城市，n[3]和n[4]则分别为紧随n[1]之前和紧随n[2]之后的两个城市的标号，其中n[3]和n[4]由函数revcst给出。
在输出端，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
#defineALEN（a,b,c,d）sqrt（（（b）-（a））*（（b）-（a））+（（d）-（c））*（（d）-（c）））
floattrncat（floatx[],floaty[],intiorder[],intncity,intn[]）
/*
该子程序返回的是输送某段给定路径所需的代价函数值。
输入参数中ncity为城市的个数；x[1..ncity]和y[1..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[1..ncity]给出当前的路径顺序；数组n共有6个元素，其意义分别为：
n[1]和n[2]分别代表将要被输送的路径段的起点城市和终点城市；n[3]和n[4]为两个相邻的城市，n[1]至n[2]路径段即将放入它们中间;n[5]和n[6]分别为n[1]之前和n[2]之后的两个城市。
在这六个元索中n[4],n[5]和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[6]和n[3]间的城市个数
nn=1;
for（j=1;j<=m1;j++）{
jj=1+（（j+n[1]-2）%ncity）;//复制所选路径段
jorder[nn++]=iorder[jj];
}
if（m2>0）{
for（j=1;j<=m2;j++）{//复制n[4]到n[5]间的路径段
jj=1+（（j+n[4]-2）%ncity）;
jorder[nn++]=iorder[jj];
}
}
if（m3>0）{
for（j=1;j<=m3;j++）{//最后，复制n[6]到n[3]间的路径段
jj=1+（（j+n[6]-2）%ncity）;
jorder[nn++]=iorder[jj];
}
}
for（j=1;j<=ncity;j++）//将jorder拷贝回iorder
iorder[j]=jorder[j];
free_ivector（jorder;1,ncity）;
}
#include
intmetrop（float,de,floatt）
/*
该子程序为米特罗波利斯算法的程序实现。
metrop返回的是一个布尔型变量，由该变量决定是否接受一个使目标函数e产生改变量de的重构形。
如果de<0则metrop=1（真），而当de>0时，metrop为真的概率是exp（一de/t），这里才是一个由退火进程决定的温度值。
*/
{
floatran3（long*idum）;
staticlonggljdum=1;
returnde<0.0||ran3（&gljdum）}
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成比例）添加到存贮的函数值中（该函数值与单纯形的每个顶点都有关联），再从每个被当作替代的新点的函数值中减去一个类似的随机变量。
和普通的米特罗波利斯方法一样，这种方法总是接受一个真正的下降步骤。
但有时也接受一次上升步骤。
在极限过程T→0中，该算法恰好变成了下降的单纯形法，并收敛到一个局部极小点。
在T的某有限值处，单纯形将扩展到一定的规模，其大小接近于在这个温度值所能达到的区域；然后单纯形在这个区域内部做随机的滚动翻转式布朗运动，并在该过程中抽取一些新的、近似随机的点作为样本。
一个区域被利用的效率与其狭窄度（对椭球状山谷，指它的主轴比率）及方向性均无关。
如果温度降低得足够缓慢，那么单纯形将极有可能收缩到那个区域内，而那个区域内包含已遇到的最低相对极小。
由于在所有模拟退火法的应用场合中，“足够缓慢”一词的含意根据问题的不同可以有相当大的细微区别，因而成功与否往往取决于退火进程选择。
下面几种可能性我们认为值一试：
升级会员 