一、模拟退火法
模拟退火法(参见[1,2])作为一种适合于求解大规模的优化问题的技术，近来已引起极大的关注。特别是当优化问题有很多局部极值而全局极值又很难求出时，模拟退火法尤其有效。在实用上，它有效地“解决了”著名的旅行推梢员问题，即在必须依次访问每一个城市(共有N个城市)的前提下，为旅行推销员设计一条能够返回起点的最短旅程。模拟退火方法还被成功地用于设计复杂的集成电路，也就是说如何最佳地安排几十万个电路元件，使它们全部集成在一个很小的硅片上，而相互连接的线路之间的缠绕能够达到最小(参见[3,4]).尽管模拟退火法的功效非凡，但它的算法实现却相对地简单，这一点似乎有些不可思议。
请注意，我们上面提到的两个例子都属于组合极小化问题。现本类问题通常也有一个目标函数，但是函数的定义域并不是简单地由N个连续参变量组成的N维空间，而是一个离散的巨大空间，例如，由所有可能的城市旅行路线组成的集合，或者硅片电路元件的所有可能的分配方式的集合。构形空间中元素的数量相当巨大，根本不可能穷举，而且因为集合是离散的，我们也不可能“沿合适的方向连续下降”。因此在构形空间中，“方向”概念就没有什么意义了。
后面我们还将介绍如何在其有连续控制参数的空间中利用模拟退火法。这种应用实际上要比组合问题复杂一些，因为其中又要出现“狭长山谷”的情况。正如在下文中我们将看到的，模拟退火法的试探步骤是“随机”的。但在一个狭窄且漫长的等高线山谷中，几乎所有的随机步骤都呈向上的趋势，因此，算法中需要增加一些技巧。
模拟退火的核心思想与热力学的原理颇为相似，而且尤其类似于液体流动和结晶以及金属冷却和退火方式。在高温下，一种液体的大最分子彼此之间进
行着相对自由的移动。如果该流体慢慢地冷却下来，热能可动性便会消失。大量原子常常能够自行排列成行，形成一个纯净的晶体，该晶体在各个方向上都被完全有序地排列在几百万倍于单个原子大小的距离之内。对于这个系统来说，晶体状态是能量最低状态；而所有缓慢冷却的系统都可以自然达到这个能量最低状态，这可以说是一个令人惊奇的事实。实际上，如果某种液体金属被迅速冷却或被“猝熄”，那么它不会达到这一状态，而只能达到一种具有较高能量的多晶状态或非结晶状态。
因此，这一过程的本质在于缓缓地致冷，以争取充足的时间，让大量原子在丧失可动性之前进行重新分布。这就是所谓退火在技术上的定义，同时也是确保达到低能量状态所必需的条件。
尽管我们的比喻并不算贴切，但是迄今为止本身所讨论的所有极小化算法，确实与快速冷却猝熄有某种关联之处。以往我们处理问题的方式都是:从初始点开始，立即沿下降方向前进，走得越远越好，似乎这样才能迅速求得问题的解。但是，正如前面常常提到的，这种方法往往只能求得局部极小点，却求不到整体最小点。自然界本身的极小化算法则基于一种截然不同的方式，所谓的玻尔兹曼(Boltzmann)概率分布
()()
ob/
Pr-（1）
E
~
kT
E
Exp
表达了这样一种思想，即:一个处于热平衡状态且具有温度T的系统，其能量按照概率，分布于所有不同能量状态E之中。即使在很低的温度下，系统也有可能（虽然这种可能性很小）处于一个较高的能量状态。因此，相应地，系统也能够获得摆脱局部能量极小点的机会。并找到一个更好的、更接近于整体的极小点。式(1)中的参数k(称为玻尔兹曼常数)是一个自然常数，它的作用是将温度与能量联系起来。换句话说，在有些情况下系统的能量可上升，也可下降，但
是温度越低，显著上升的可能性就越小。
1953年，米特罗波利斯(Metropolis)及其合作者们首次将这种原理渗透到数值计算中。他们对一个模拟热力学系统提供了一系列选择项，并假设：系统构形从能量1E 变化到能量2E 的概率为()[]kT E E Exp p 12--=。很显然，如果12E E <，p 将大于1；在这类情况下，对构形的能量变化任意指定一个概率值1=p ，也就是说，该系统总是取这个选择项。这种格式总是采取下降过程，偶尔采取上升步骤。目前已被公认为米特罗波利斯算法。
为了将米特波利斯算法应用于热力学以外的系统，必须提供以下几项基本要素：
1. 对可能的系统构形的一种描述。
2. 一个有关构形内部随机变化的生成函数，这些变化将作为“选择项”提交给该系统。
3. 一个目标函数E (类似于能量)，求解E 的极小值，即为算法所要完成的工作。
4. 一个控制参数T(类似于温度)和一个退火进程，该进程用来说明系统是如何从高值向低值降低的，例如在温度T 时每次下降步骤中要经过多少次随机的构形变化以及该步长是多大等等。应说明的是，这里“高”和“低”的含义，还有进程表的确定，都需要一定的物理知识和/或一些摸索的实验。
模拟退火算法的模型
模拟退火算法可以分解为解空间、目标函数和初始解三部分。
模拟退火的基本思想:
(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)， 每个T 值的迭代次数L
(2) 对k=1，……，L做第(3)至第6步：
(3) 产生新解S′
(4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数
(5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)接受S′作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解，结束程序。
终止条件通常取为连续若干个新解都没有被接受时终止算法。
(7) T逐渐减少，且T->0，然后转第2步。
模拟退火算法新解的产生和接受可分为如下四个步骤：
第一步，是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
第二步，是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
第三步，是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropo1is准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。
第四步，是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃
时，则在原当前解的基础上继续下一轮试验。
模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率1 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。
模拟退火算法的简单应用
作为模拟退火算法应用，讨论货郎担问题(Travelling Salesman Problem ，简记为TSP)：设有n 个城市，用数码1,…,n 代表。城市i 和城市j 之间的距离为d(i ，j) i, j=1,…,n ．TSP 问题是要找遍访每个域市恰好一次的一条回路，且其路径总长度为最短.。
求解TSP 的模拟退火算法模型可描述如下：
解空间：解空间S 是遍访每个城市恰好一次的所有路经，是{1，……，
n}的所有循环排列的集合，S 中的成员记为{}n w w w ,,21，
n w w w ,,21是1,2,……,n 的一个排列，表明从1w 城市出发，依次经过n w w w ,,32城市，再返回1w 城市。初始解可选为(1,……, n) ;
目标函数：目标函数为访问所有城市的路径总长度或称为代价函数; 我们要求的最优路径为目标函数(代价函数)为最小值时对应的路径。 新路径的产生：随机产生1和n 之间的两相异数k 和m ，不妨假设k ()n m m k k w w w w w w w ,1121,,,,,++ 变为新路径： ()n m k k m w w w w w w w ,1121,,,,,++。 这种变换方法就是将k 和m 对应的两个城市在路径序列中交换位置,称为2-opt 映射。 根据上述描述，模拟退火算法求解TSP问题的流程框图如下 图2 模拟退火算法的流程框图 也可以采用其他的变换方法，有些变换有独特的优越性，有时也将它们交替使用，得到一种更好方法。例如： 随机产生1和n 之间的两相异数k 和m ，不妨假设k ()n m m k k w w w w w w w ,1121,,,,,++ 变为新路径： ()n k k k m m w w w w w w w w 2,1121,,,,,,++- 上述变换方法就是将k 和m 之间对应的多个城市在路径序列颠倒位置, 可简单说成是“逆转中间或者逆转两端”。 根据上述分析，可写出用模拟退火算法求解TSP 问题的伪程序： Procedure TSPSA: begin init-of-T; // T 为初始温度 S={1，……，n}; //S 为初始值 termination=false; while termination=false begin for i=1 to L do begin generate(S′form S); //从当前回路S 产生新回路S′ Δt:=f(S′))-f(S); //f(S)为路径总长 IF(Δt<0) OR (EXP(-Δt/T)>Random -of-[0,1]) S=S′; IF the-halt-condition-is-TRUE THEN termination=true; End; T_lower; End; End 模拟退火算法的应用很广泛，可以较高的效率求解最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)、调度问题(Scheduling Problem)等等。 模拟退火算法的参数控制问题 模拟退火算法的应用很广泛，可以求解NP完全问题，但其参数难以控制，其主要问题有以下三点： (1) 温度T的初始值设置问题。 温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。实际应用过程中，初始温度一般需要依据实验结果进行若干次调整。 (2) 退火速度问题。 模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索(退火)是相当必要的，但这需要计算时间。实际应用中，要针对具体问题的性质和特征设置合理的退火平衡条件。 (3) 温度管理问题。 温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中，由于必须考虑计算复杂度的切实可行性等问题，常采用如下所示的降温方式：T(t+1)＝k×T(t) 式中k为正的略小于1.00的常数，t为降温的次数 二、有记忆的模拟退火法 传统的模拟退火算法思路清晰、原理简单、可用于解决许多实际问题。但存在很多弊病，国内外学者并对其作相应的改进，得到了改进的模拟退火算法，包括：加温退火法、有记忆的模拟退火算法、带返回搜索的模拟退火算法、多次寻优法、回火退火算法等。针对特征选择的特点，本文选用了有记忆的模拟退火算法作为搜索策略。 在传统的模拟退火过程中，算法终止于一个预先规定的停止准则S ，如：控制参数t 的值小于某个充分小的正数；相继的若干个Markov 链中解未得到任何改善；两个相继Markov 链所得解的差的绝对值小于某个充分小的正数等。但由于模拟退火算法的搜索过程是随机的，且当t 值较大时可以接受部分恶化解，而随着t 值的衰减，恶化解被接受的概率逐渐减小直至为0。另一方面，某些当前解要达到最优解时必须经过暂时恶化的“山脊”，因此，上述这些停止准则无法保证最终解正好是整个搜索过程中曾经达到的最优解。 因此，在传统的算法中增加一个记忆器，使之能够记住搜索过程中遇到过的最好解，当退火结束时，将所得最终解与记忆器中的解比较并取较优者作为最后结果。 三、组合极小化：旅行推销员问题 下面是我们用“旅行推销员问题”为具体实例说明模拟退火法的应用。假设一个推销员，要去N 个分别位于()i i y x ,的城市进行推销，并于最后返回他原来 所在的城市，要求每个城市只能去一次，而且所经过的路径要尽可能地短。这个间题属于一类所谓“NP-完全问题”。这类问题求出一个精确解所需的计算时间是随N 的增加以指数exp(常数×N)增长的。当N 不断增大时，运行时间将迅速增加，进而导致费用高到令人难以接受的程度。旅行推销员问题也是众多极小化问题中的一种，它的目标函数E 具有多个局部极小值。在实际应用当中，常常有足够多的条件可以从多个极小值中选出一个最小的，这个最小值即使不是绝对最小，也相当接近绝对最小了。退火法的目的就是要获得这个最小值，同时又要将计算量限制在N 的低阶次的数量级上。 旅行推销员问题也是按模拟退火问题的方式进行处理的。具体如下： 1.构形 将N 个城市分别标记为N i ,2,1=中的数，其中每个城市具有坐标 ()i i y x ,。 一个构形就是数字N 1的一个排列，可以解释为推销员途径的城市的顺序。 2.调整 林(Lin)曾经提出过一种所谓“转移有效集”，这里的“转移”包括两种类型：(a)移走路径的某一段，然后对这段路径上的城市用相反次序重新进行排列，并用后者来代替前者。(b)移走某段路径，并用位于城市间的随机选取的另一段路径来取代被移走的路径。 3.目标函数 在旅行推销员问题的一种最简单的形式中，目标函数E 被定义为旅途的总长度 ()()∑=++-+-≡=N i i i i i y y x x L E 12121 （10.9.2） 这里认为第N+1个点与第1个点是重合的。但是为了表明模拟退火法的灵活性，我们还要用到下面的技巧：假设推销员无端端地害怕飞越密西西比河，在这种情况下，我们对每个城市给定一个参数i μ，如果该城市位于密西西比河 以东，i μ取1；若在密西西比河以西则取-1。对于目标函数E ，我们将其改写为： ()()()∑=+++?? ????-+-+-=N i i i i i i i y y x x E 1212121μμλ （10.9.3） 由于每过一次河都将以λ4作为惩罚，因而现在我们设计的算法的目标，就变成了寻找尽可能回避过河的最短路径。路径长度对过河次数的相对重要性将由我们选择的λ来确定。图10. 9. 1表明了所得的结果。显然，这种技巧可以推广到包含许多相互冲突的目的要求的极小化问题当中。 (A) (B) (C) 图(A)表示的是从四个随机分布的城市中间找到的一条(接近)最短路径，中间竖虚线标识的是一条河流，但这是对过河没有附加您罚项的情况。在图(B)中对过河施加的惩罚项很大，而图中所示的解本身的过河次数也相应地只有少得不能再少的两次。在图(C)中惩罚项为负，这就是说，推销员实际上成了恣意偷渡的走私者！ 图1 用模拟退火法解决旅行推销员问题 4.退火进程这一步需要借助试验来确定。首先要进行一些随机调整，然后利用它们来确定从转移到转移过程中将会遇到的E ?值之范围。对参数T取一初始值(这个初始值要远远大于通常所能遇到的E ?的最大值)，并以倍增的步长下减，每次使T总共减少10%。我们拿每个新的常数T值去试各种1OON重构形，或1ON成功的重构形，无论哪个在前出现就取哪个。当实在不能再进一步减小E 时，则停止。 下面的旅行推销员程序利用了米特罗波利斯(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个城市的位置坐标存贮在数组x[l..ncity]和y[1..ncity]中。数组iorder[1..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 n[7]; long idum; unsigned long iseed; float path,de,t; nover=100*ncity;//在任何温度下，所试路径的最大次数 nlimit=10*ncity;//在继续进行之前，成功的路径改变的最大次数 path=0.0; t=0.5; for(i=1;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; 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," Path Length=",path); printf("Successful Moves: %6d\n"，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为城市数；数组x[1..ncity]，y[1..ncity]为这些城市的位置坐标；iorder[n.. ncity]为当前路线；数组n的头两个元素n[1]和n[2]分别代表将要被反转的路径段的起点城市和终点城市。子程序的输出项de为反转所需的代价值，但真正的反转过程并不是由该程序执行。 */ { float xx[5],yy[5],de; int i,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]); return de; } void reverse (int iorder[]，int ncity，int n[]) /* 该子程序的作用是执行一段路径的反转过程。输入参数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]路段被反转后的行程路线。 */ { int nn,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 # 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为城市的个数；x[1..ncity]和y[1..ncity] 为这些城市的位置坐标，数组n的前三个元索分别为：将要被输送的路径段的起、止点城市以及这段路径将要被插入处的标志城市(插在该城市之后)，该子程序的输出项de为计算出来的输送代价值，但实际的输送操作井不由该子程序执行。 */ { float xx[7],yy[7],de; int j,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]); return de; } # include "nrutil.h" void trnspt(int iorder[], int ncity, int n[]) /* 该子程序的作用是执行真正的段输送操作，输入参数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将根据路径段的移动和变化作出相应的修改。 */ { int ml,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=l; 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 int metrop (float, de, float t ) /* 该子程序为米特罗波利斯算法的程序实现。metrop返回的是一个布尔型变量，由该变量决定是否接受一个使目标函数e产生改变量de的重构形。如果de<0则metrop=1(真)，而当de>O时，metrop为真的概率是exp(一de/t)，这里才是一个由退火进程决定的温度值。 */ { float ran3(long *idum); static long gljdum =1 ; return de<0.0 || ran3(&gljdum) < exp(-de/t); } 四、模拟退火法在连续极小化问题中的应用 模拟退火法的基本思想也可以应用于具有连续N维控制参数的极小化问题当中，例如，在某个函数()x f(这里x为一个N维向量)有许多局部极小点的情况下，求解它的(理想的全局的)极小值。这时米特罗波利斯算法所需的四要素可以具体化为： 1)目标函数即为()x f的值； 2)系统状态描态即为N维空间中的点x； 3)类似于温度的控制参数T以及一个使T逐渐降低的退火进程仍为原先的定义； 4)描述构形内部随机变化的发生器即为一个从x到x+△x采取随机步骤的方法。 在上述四要素中问题最大的是最后一条。目前已发表的文献中[7-10]，介绍了几种不同的选择△x办法。但我们认为，这些方法都不算成功。问题在于“效率”二字：当局部的向下运动存在时，如果某个随机变化发生器几乎总是作出向上运动的决策，那么它的效率就很低。我们认为，一个好的发生器在等高线的“窄谷”中仍应保持高效性；当算法在接近收敛到极小点处时，它的效率也不应越变越低。在上面我们提到的几篇文献中，除[7]中介绍的方法之外，其他所有方法都表现不同程度的低效性。 下面我们将要介绍的这种方法，利用了下降单纯形法的一种修改后的形式。首先我们将米特罗波利斯算法四要素中的系统状态描述，由单个点x改为一个其有N+1个点的单纯形；单纯形的“操作”分为反射、扩张和收缩三种。然后我们将一个正的、呈对数分布的随机变量(与温度T成比例)添加到存贮的函数值中(该函数值与单纯形的每个顶点都有关联)，再从每个被当作替代的新点的函数值中减去一个类似的随机变量。和普通的米特罗波利斯方法一样，这种方法总是接受一个真正的下降步骤。但有时也接受一次上升步骤。在极限过程T→0中，该算法恰好变成了下降的单纯形法，并收敛到一个局部极小点。 在T的某有限值处，单纯形将扩展到一定的规模，其大小接近于在这个温度值所能达到的区域；然后单纯形在这个区域内部做随机的滚动翻转式布朗运动，并在该过程中抽取一些新的、近似随机的点作为样本。一个区域被利用的效率与其狭窄度(对椭球状山谷，指它的主轴比率)及方向性均无关。如果温度降 低得足够缓慢，那么单纯形将极有可能收缩到那个区域内，而那个区域内包含已遇到的最低相对极小。 由于在所有模拟退火法的应用场合中，“足够缓慢”一词的含意根据问题的不同可以有相当大的细微区别，因而成功与否往往取决于退火进程选择。下面几种可能性我们认为值一试： ·每经过m 步移动之后，将T 减到(1ε-)T ，这里。ε/m 的具体值要通过实验确定。 ·设总的移动步数为K ，每经过m 步移动之后将T 减到()αT k T T -=10 ，其中k 是到目前为止所经过的步数的累加值。α为一常数，可取为1，2: 4等。α的最佳值取决于各种深度的相对于极小的统计分布，稍大一些的α值在较低温度时，需化费的迭代次数将更多； ·每经过m 步骤移动之后，置T 为β乘b f f -1，其中β为一个阶数为1的常 数，其具体值由试验确定；1f 为目前单纯形中最小的函数值；b f 为曾经遇到的最 佳函数值。但应注意的是，T 的降低幅度一次不要超过某个分数值γ。 另一个策略方法上的问题是，当单纯形的某个顶点被放弃并让位于“永远的最佳点”时，是否需要采取重新开始的步骤(但我们必须保证在进行这项工作时，这个“永远的最佳点”当前并不在单纯形中!)。对于有些问题重新开始(例如，只要温度降低了因子3即执行重新开始步骤)的效果极佳，而对于另外一些问题，重新开始不仅没有任何效果反而会产生负作用。上述两种截然相反的情况，我们都找到了例子可以作为例证。 将下面的程序amebsa 采用改进下降单纯形法。 # include # include "nrutil.h" # define GET_PSUM \ for (n=1;n<=ndim;n++) {\ for (sum=0.0,m=1l;m<=mpts;m++) sum += p[m][n];\ psum[n]=sum; } extern long idum; //由主程序定义和初始值 float tt ;//与amotsa进行通信 void amebsa (float **p,float y[],int ndim, float pb[],float *yb,float ftol, float (*funk)(float []),int *iter,float temptr) /* 用摸拟退火与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]将依次返回已遇到的最佳函数值和最佳点(这个最佳点可以水远不是单纯形中的点)。 */ { float amotsa(float **p, float y[], float psum[], int ndim, float pb[], float *yb, float (*funk) (float []), int ihi, float *yhi, float fac); float ran1 (long *idum); int i,ihi,ilo,j,m,n,mpts=nidm+1; float rtol,sum,swap,yhi,ylo,ynhi,ysave,yt,ytry, *psum; psum=vector (1,ndim); tt = -temptr; GET_PSUM for (;;) { ilo=1;//确定最高点(即差点)、次高点和最低点(即最佳点) ihi=2; ynhi=ylo=y[1]+tt*log(ran1(&idum));//我们所“看到”的顶点总是处于 //随机的热起伏状态 yhi=y[2]+tt * log(ran1(&idum)); 模拟退火算法模拟退火算法是一种通用的随机搜索算法，是局部搜索算法的扩展。它的思想是再1953 年由metropolis 提出来的，到1983 年由kirkpatrick 等人成功地应用在组合优化问题中。 模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis 准则，粒子在温度T 时趋于平衡的概率为e- △ E/(kT)，其中E为温度T时的内能，AE为其改变量，k 为Boltzmann 常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f ,温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解-计算目标函数差T接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooli ng Schedule)控制，包括控制参数的初值t 及其衰减因子△ t、每个t值时的迭代次数L和停止条件S。 模拟退火算法新解的产生和接受可分为如下四个步骤：第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。 第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。 第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropo1is准则：若厶t‘ <0 则接受S'作为新的当前解S,否则以概率exp(- △ t‘ /T) 接受S'作为新的当前解S。 第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。 可在此基础上开始下一轮试验。而当新解被判定为舍弃时， 实验用例： 用模拟退火算法解决如下10个城市的TSP 问题，该问题最优解为691.2 opt f 。 表1 10个城市的坐标 城市 X 坐标 Y 坐标 城市 X 坐标 Y 坐标 3 0.4000 0.4439 8 0.8732 0.6536 编程实现 用MATLAB 实现模拟退火算法时，共编制了5个m 文件，分别如下 1、swap.m function [ newpath , position ] = swap( oldpath , number ) % 对 oldpath 进 行 互 换 操 作 % number 为 产 生 的 新 路 径 的 个 数 % position 为 对 应 newpath 互 换 的 位 置 m = length( oldpath ) ; % 城 市 的 个 数 newpath = zeros( number , m ) ; position = sort( randi( m , number , 2 ) , 2 ); % 随 机 产 生 交 换 的 位 置 for i = 1 : number newpath( i , : ) = oldpath ; % 交 换 路 径 中 选 中 的 城 市 newpath( i , position( i , 1 ) ) = oldpath( position( i , 2 ) ) ; newpath( i , position( i , 2 ) ) = oldpath( position( i , 1 ) ) ; end 2、pathfare.m function [ objval ] = pathfare( fare , path ) % 计 算 路 径 path 的 代 价 objval % path 为 1 到 n 的 排 列 ，代 表 城 市 的 访 问 顺 序 ； % fare 为 代 价 矩 阵 ， 且 为 方 阵 。 [ m , n ] = size( path ) ; objval = zeros( 1 , m ) ; for i = 1 : m for j = 2 : n objval( i ) = objval( i ) + fare( path( i , j - 1 ) , path( i , j ) ) ; end objval( i ) = objval( i ) + fare( path( i , n ) , path( i , 1 ) ) ; end 2014中南大学数学建模模拟竞赛第一轮 承诺书 我们仔细阅读了中国大学生数学建模竞赛的竞赛规则. 我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。 我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。 我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。 我们参赛选择的题号是（从A/B/C/D中选择一项填写）： B 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）：中南大学 参赛队员(打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人(打印并签名)： 日期： 2014 年 8 月 11 日赛区评阅编号（由赛区组委会评阅前进行编号）： 2014高教社杯全国大学生数学建模竞赛 编号专用页 赛区评阅编号（由赛区组委会评阅前进行编号）： 全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）： 商品期货交易策略 摘要 我国的期货发展历史已有十多年，吸引了大量交易者的参与，如何从中获取相对稳定的收益成为交易者非常关注的问题。本文旨在为交易者谋得最大盈利，通过数据分析，找到影响价格因素，对价格波动进行分类并预测，从而建立交易模型。 本文通过对数据抽样，拟合检验，建立主成分分析模型（模型1），找到影响价格因素指标，回归分析检验结果；再建立聚类分析模型（模型2），对波动方式进行分类，并建立小波神经网络预测模型（模型3）对价格趋势作出预测，最后建立期货获利交易模型（模型4），使交易者获得最大盈利。 模型1：主成分分析模型 由于对价格有影响的因素众多，而由SPSS得到的散点图和相关系数表可发现，成交价与B1价、S1价和日期有极其显著的关系，但许多变量之间可能存在信息上的重叠。故选用了主成分分析模型，进行贡献率的判定。利用SPSS软件，将数据标准化（数据见附件1），并获得相关系数表和特征方程，提取特征值大于1的前4 个主成分，通过计算可得到每个主成分前的系数，即特征向量。计算可得出主成分表达式。最后可由主成分综合模型中根据每个因素的贡献率判定对价格的影响因素。最后利用MatlabR2012a 软件进行回归分析检验。 模型2：聚类分析模型 为找到不同波动方式的类型，先利用MatlabR2012a软件绘出时间-盈利走势图，在此基础上选择盈利最大周期，3个交易日；然后选择R性聚类分析，对变量进行相似性度量，对相似性大的变量进行聚类。利用SPSS软件，将10个相关变量进行组内链接，皮尔逊相关测量区间的相关性方法作出聚类图，共分为8组（表2），最后给出分析得到的交易量、持仓量和价格的关系。 模型3：小波神经网络预测模型 为了对价格的后期走势作出预测，按交易者的投资来看必然是短期预测，故采用精确度较高的小波神经网络进行预测。利用MatlabR2012a软件，选取3个输入节点，6个隐含层节点和1个输出节点，对9天的数据进行训练，修正，另外10天的数据进行预测，分别反复训练200次和500次，得到预测结果与实际结果高精确度吻合（见图4-5），说明该预测模型合理。 模型4：期货获利交易模型 根据前两问得出价格相关因素和价格的预测，为使交易者盈利最大，建立期货获利交易模型，在原先盈利函数上扣除手续费、保证金，利用线性规划方法，设立约束条件，目标函数为最大盈利，最后利用MatlabR2012a软件进行求解得到月所有日最大收益为515700元。 关键字：主成分分析聚类分析小波神经网络预测期货获利价格波动最大盈利 例已知敌方100个目标的经度、纬度如下： 我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。 这是一个旅行商问题。我们依次给基地编号为1，敌方目标依次编号为2，3，…，101，最后我方基地再重复编号为102（这样便于程序中计算）。距离矩阵102102)(?=ij d D ，其中ij d 表示表示j i ,两点的距离，102,,2,1, =j i ，这里D 为实对称矩阵。则问题是求一个从点1出发，走遍所有中间点，到达点102的一个最短路径。 上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设B A ,两点的地理坐标分别为),(11y x ，),(22y x ，过B A ,两点的大圆的劣弧长即为两点的实际距离。以地心为坐标原点O ，以赤道平面为XOY 平面，以0度经线圈所在的平面为XOZ 平面建立三维直角坐标系。则B A ,两点的直角坐标分别为： )s i n ,c o s s i n ,c o s c o s (11111y R y x R y x R A )s i n ,c o s s i n ,c o s c o s (22222y R y x R y x R B 其中6370=R 为地球半径。 B A ,两点的实际距离 ?? =R d arccos ， 化简得 ]s i n s i n c o s c o s )(a r c c o s [c o s 212121y y y y x x R d +-=。 求解的模拟退火算法描述如下： （1）解空间 解空间S 可表为{102,101,,2,1 }的所有固定起点和终点的循环排列集合，即 }102,}101,,3,2{),,(,1|),,{(102101211021===ππππππ的循环排列为 S 其中每一个循环排列表示侦察100个目标的一个回路，j i =π表示在第i 次侦察j 点，初始解可选为)102,,2,1( ，本文中我们使用Monte Carlo 方法求得一个较好的初始解。 （2）目标函数 此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求 ∑=+= 101 1 211 102),,,(min i i i d f π ππππ 一.爬山算法( Hill Climbing ) 介绍模拟退火前，先介绍爬山算法。爬山算法是一种简单的贪心搜索算 法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到 达到一个局部最优解。 爬山算法实现很简单，其主要缺点是会陷入局部最优解，而不一定能搜 索到全局最优解。如图1所示：假设C点为当前解，爬山算法搜索到A点 这个局部最优解就会停止搜索，因为在A点无论向那个方向小幅度移动都不 能得到更优的解。 二.模拟退火(SA,Simulated Annealing)思想（跟人一样找不 到最优解就最产生疑惑，我到底需不需要坚持，随着时间的推移，逐渐的慢慢的放弃去追寻最优解的念头） 爬山法是完完全全的贪心法，每次都鼠目寸光的选择一个当前最优解，因此只能搜索到局部的最优值。模拟退火其实也是一种贪心算法，但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。 以图1为例，模拟退火算法在搜索到局部最优解A后，会以一定的概率接受到E的移动。也许经过几次这样的不是局部最优的移动后会到达D点，于是就跳出了局部最大值A。 若J( Y(i+1) )>= J( Y(i) ) (即移动后得到更优解)，则总是接受该移动 若J( Y(i+1) )< J( Y(i) ) (即移动后的解比当前解要差)，则以一定的概率接受移动，而且这个概率随着时间推移逐渐降低（逐渐降低才能趋向稳定） 这里的“一定的概率”的计算参考了金属冶炼的退火过程，这也是模拟退火算法名称的由来。 根据热力学的原理，在温度为T时，出现能量差为dE的降温的概率为P(dE)，表示为： P(dE) = exp( dE/(kT) ) 其中k是一个常数，exp表示自然指数，且dE<0。这条公式说白了就是：温度越高，出现一次能量差为dE的降温的概率就越大；温度越低，则出现降温的概率就越小。又由于dE总是小于0（否则就不叫退火了），因此dE/kT < 0 ，所以P(dE)的函数取值范围是(0,1) 。 随着温度T的降低，P(dE)会逐渐降低。 我们将一次向较差解的移动看做一次温度跳变过程，我们以概率P(dE)来接受这样的移动。 关于爬山算法与模拟退火，有一个有趣的比喻：（有点意思） 2012年数学建模竞赛试题 注意事项（请参赛队员详细阅读！） 1. 凯里学院校内数学建模竞赛丁2012年6月29日8: 00至7月 1日20 : 00举行。 2. 参赛队可在A、B两题中任选其中一题，可以使用各种图书资料、网络信息、计算机和软件以及各种实验手段。 3. 答卷论文请提交WORD文档方式的A4纸电子稿。并按下列要求制作。 论文用白色A4纸单面打印；上下左右各留出至少 2.5厘米的贞边距； 从左侧装订。 封面：只需填上所选论文题目（注明A或B）及参赛队序号，其他一律不要。 首页：论文题目、摘要（含模型的主要特点、建模方法和主要结果）。 正文：问题提出、问题分析、模型假设、符号说明、模型建立、模型求 解、计算方法设计和软件实现、模型结果分析和检验、模型优缺点分析等。 4. 论文从第三页开始编写贞码，贞码必须位丁每贞贞脚中部，用阿拉伯数字从“ 1”开始连续编号。 论文题目用三号黑体字、一级标题用四号黑体字，并居中；二级、三 级标题用小四号黑体字，左端对齐（不居中）。论文中其他汉字一律采用 小四号宋体字，行距用单倍行距，打印时应尽量避免彩色打印。 提请大家注意：摘要应该是一份简明扼要的详细摘要（包括关键词）， 在整篇论文评阅中占有重要权重，请认真书写（注意篇幅不能超过一页，且无需译成英文）。评阅时将首先根据摘要和论文整体结构及概貌对论文优劣进行初步筛选引用别人的成果或其他公开的资料（包括网上查到的资料）必须按照规定的参考文献的表述方式在正文引用处和参考文献中均明确列出。正文引用处用方括号标示参考文献的编号，如［1］［3］等；引用书籍还必须指出贞码。参考文献按正文中的引用次序列出，其中书籍的表述方式为： ［编号］作者，书名，出版地：出版社，出版年。 参考文献中期刊杂志论文的表述方式为： ［编号］作者，论文名，杂志名，卷期号：起止贞码，出版年。 参考文献中网上资源的表述方式为： ［编号］作者，资源标题，网址，访问时间（年月日）。 5. 竞赛评奖以模型假设的合理性、建模的创造性、结果的正确性、文字表述的活晰程度为主要标准。 6. 答卷（电子稿）务必丁2012年7月1日20:00 —22:00交到凯里学院数学实验室潘东云或雷学红老师处。 凯里学院数学建模领导小组 2012年06月28日 数学模型的分类 按模型的数学方法分： 几何模型、图论模型、微分方程模型、概率模型、最优控制模型、规划论模型、马氏链模型等 按模型的特征分： 静态模型和动态模型，确定性模型和随机模型，离散模型和连续性模型，线性模型和非线性模型等 按模型的应用领域分： 人口模型、交通模型、经济模型、生态模型、资源模型、环境模型等。 按建模的目的分： 预测模型、优化模型、决策模型、控制模型等 一般研究数学建模论文的时候，是按照建模的目的去分类的，并且是算法往往也和建模的目的对应 按对模型结构的了解程度分： 有白箱模型、灰箱模型、黑箱模型等 比赛尽量避免使用，黑箱模型、灰箱模型，以及一些主观性模型。 按比赛命题方向分： 国赛一般是离散模型和连续模型各一个，2016美赛六个题目（离散、连续、运筹学/复杂网络、大数据、环境科学、政策） 数学建模十大算法 1、蒙特卡罗算法 （该算法又称随机性模拟算法，是通过计算机仿真来解决问题的算法，同时可以通过模拟可以来检验自己模型的正确性，比较好用的算法） 2、数据拟合、参数估计、插值等数据处理算法 （比赛中通常会遇到大量的数据需要处理，而处理数据的关键就在于这些算法，通常使用Matlab作为工具） 3、线性规划、整数规划、多元规划、二次规划等规划类问题 （建模竞赛大多数问题属于最优化问题，很多时候这些问题可以用数学规划算法来描述，通常使用Lindo、Lingo软件实现） 4、图论算法 （这类算法可以分为很多种，包括最短路、网络流、二分图等算法，涉及到图论的问题可以用这些方法解决，需要认真准备） 5、动态规划、回溯搜索、分治算法、分支定界等计算机算法 （这些算法是算法设计中比较常用的方法，很多场合可以用到竞赛中） 6、最优化理论的三大非经典算法：模拟退火法、神经网络、遗传算法 （这些问题是用来解决一些较困难的最优化问题的算法，对于有些问题非常有帮助，但是算法的实现比较困难，需慎重使用） 7、网格算法和穷举法 （当重点讨论模型本身而轻视算法的时候，可以使用这种暴力方案，最好使用一些高级语言作为编程工具） 8、一些连续离散化方法 （很多问题都是从实际来的，数据可以是连续的，而计算机只认的是离散的数据，因此将其离散化后进行差分代替微分、求和代替积分等思想是非常重要的） 9、数值分析算法 （如果在比赛中采用高级语言进行编程的话，那一些数值分析中常用的算法比如方程组求解、矩阵运算、函数积分等算法就需要额外编写库函数进行调用） 10、图象处理算法 （赛题中有一类问题与图形有关，即使与图形无关，论文中也应该要不乏图片的这些图形如何展示，以及如何处理就是需要解决的问题，通常使用Matlab进行处理） 算法简介 1、灰色预测模型（必掌握） 解决预测类型题目。由于属于灰箱模型，一般比赛期间不优先使用。 满足两个条件可用： ①数据样本点个数少，6-15个 ②数据呈现指数或曲线的形式 2、微分方程预测（高大上、备用） 微分方程预测是方程类模型中最常见的一种算法。近几年比赛都有体现，但其中的要求，不言而喻。学习过程中 无法直接找到原始数据之间的关系，但可以找到原始数据变化速度之间的关系，通过公式推导转化为原始数据的关系。 3、回归分析预测（必掌握） 求一个因变量与若干自变量之间的关系，若自变量变化后，求因变量如何变化； 样本点的个数有要求： ①自变量之间协方差比较小，最好趋近于0，自变量间的相关性小； ②样本点的个数n＞3k+1，k为自变量的个数； 数学建模常用的十大算法==转 (2011-07-24 16:13:14) 转载▼ 1. 蒙特卡罗算法。该算法又称随机性模拟算法，是通过计算机仿真来解决问题的算法，同时可以通过模拟来检验自己模型的正确性，几乎是比赛时必用的方法。 2. 数据拟合、参数估计、插值等数据处理算法。比赛中通常会遇到大量的数据需要处理，而处理数据的关键就在于这些算法，通常使用MA TLAB 作为工具。 3. 线性规划、整数规划、多元规划、二次规划等规划类算法。建模竞赛大多数问题属于最优化问题，很多时候这些问题可以用数学规划算法来描述，通常使用Lindo、Lingo 软件求解。 4. 图论算法。这类算法可以分为很多种，包括最短路、网络流、二分图等算法，涉及到图论的问题可以用这些方法解决，需要认真准备。 5. 动态规划、回溯搜索、分治算法、分支定界等计算机算法。这些算法是算法设计中比较常用的方法，竞赛中很多场合会用到。 6. 最优化理论的三大非经典算法：模拟退火算法、神经网络算法、遗传算法。这些问题是用来解决一些较困难的最优化问题的，对于有些问题非常有帮助，但是算法的实现比较困难，需慎重使用。 7. 网格算法和穷举法。两者都是暴力搜索最优点的算法，在很多竞赛题中有应用，当重点讨论模型本身而轻视算法的时候，可以使用这种暴力方案，最好使用一些高级语言作为编程工具。 8. 一些连续数据离散化方法。很多问题都是实际来的，数据可以是连续的，而计算机只能处理离散的数据，因此将其离散化后进行差分代替微分、求和代替积分等思想是非常重要的。 9. 数值分析算法。如果在比赛中采用高级语言进行编程的话，那些数值分析中常用的算法比如方程组求解、矩阵运算、函数积分等算法就需要额外编写库函数进行调用。 10. 图象处理算法。赛题中有一类问题与图形有关，即使问题与图形无关，论文中也会需要图片来说明问题，这些图形如何展示以及如何处理就是需要解决的问题，通常使用MA TLAB 进行处理。 以下将结合历年的竞赛题，对这十类算法进行详细地说明。 以下将结合历年的竞赛题，对这十类算法进行详细地说明。 2 十类算法的详细说明 2.1 蒙特卡罗算法 大多数建模赛题中都离不开计算机仿真，随机性模拟是非常常见的算法之一。 举个例子就是97 年的A 题，每个零件都有自己的标定值，也都有自己的容差等级，而求解最优的组合方案将要面对着的是一个极其复杂的公式和108 种容差选取方案，根本不可能去求解析解，那如何去找到最优的方案呢？随机性模拟搜索最优方案就是其中的一种方法，在每个零件可行的区间中按照正态分布随机的选取一个标定值和选取一个容差值作为一种方案，然后通过蒙特卡罗算法仿真出大量的方案，从中选取一个最佳的。另一个例子就是去年的彩票第二问，要求设计一种更好的方案，首先方案的优劣取决于很多复杂的因素，同样不可能刻画出一个模型进行求解，只能靠随机仿真模拟。 2.2 数据拟合、参数估计、插值等算法 数据拟合在很多赛题中有应用，与图形处理有关的问题很多与拟合有关系，一个例子就是98 年美国赛A 题，生物组织切片的三维插值处理，94 年A 题逢山开路，山体海拔高度的插值计算，还有吵的沸沸扬扬可能会考的“非典”问题也要用到数据拟合算法，观察数据的 在数学建模中常用的方法：类比法、二分法、量纲分析法、差分法、变分法、图论法、层次分析法、数据拟合法、回归分析法、数学规划（线性规划，非线性规划，整数规划，动态规划，目标规划）、机理分析、排队方法、对策方法、决策方法、模糊评判方法、时间序列方法、灰色理论方法、现代优化算法（禁忌搜索算法，模拟退火算法，遗传算法，神经网络）。 用这些方法可以解下列一些模型：优化模型、微分方程模型、统计模型、概率模型、图论模型、决策模型。拟合与插值方法（给出一批数据点，确定满足特定要求的曲线或者曲面，从而反映对象整体的变化趋势）：matlab可以实现一元函数，包括多项式和非线性函数的拟合以及多元函数的拟合，即回归分析，从而确定函数；同时也可以用matlab实现分段线性、多项式、样条以及多维插值。 在优化方法中，决策变量、目标函数（尽量简单、光滑）、约束条件、求解方法是四个关键因素。其中包括无约束规则（用fminserch、fminbnd实现）线性规则（用linprog实现）非线性规则、（用fmincon实现）多目标规划（有目标加权、效用函数）动态规划（倒向和正向）整数规划。 回归分析：对具有相关关系的现象，根据其关系形态，选择一个合适的数学模型，用来近似地表示变量间的平均变化关系的一种统计方法（一元线性回归、多元线性回归、非线性回归），回归分析在一组数据的基础上研究这样几个问题：建立因变量与自变量之间的回归模型（经验公式）；对回归模型的可信度进行检验；判断每个自变量对因变量的影响是否显著；判断回归模型是否适合这组数据；利用回归模型对进行预报或控制。相对应的有线性回归、多元二项式回归、非线性回归。 逐步回归分析：从一个自变量开始，视自变量作用的显著程度，从大到地依次逐个引入回归方程：当引入的自变量由于后面变量的引入而变得不显著时，要将其剔除掉；引入一个自变量或从回归方程中剔除一个自变量，为逐步回归的一步；对于每一步都要进行值检验，以确保每次引入新的显著性变量前回归方程中只包含对作用显著的变量；这个过程反复进行，直至既无不显著的变量从回归方程中剔除，又无显著变量可引入回归方程时为止。（主要用SAS来实现，也可以用matlab软件来实现）。 聚类分析：所研究的样本或者变量之间存在程度不同的相似性，要求设法找出一些能够度量它们之间相似程度的统计量作为分类的依据，再利用这些量将样本或者变量进行分类。 系统聚类分析—将n个样本或者n个指标看成n类，一类包括一个样本或者指标，然后将性质最接近的两类合并成为一个新类，依此类推。最终可以按照需要来决定分多少类，每类有多少样本（指标）。 系统聚类方法步骤： 1.计算n个样本两两之间的距离 2.构成n个类，每类只包含一个样品 3.合并距离最近的两类为一个新类 4.计算新类与当前各类的距离（新类与当前类的距离等于当前类与组合类中包含的类的距离最小值）， 若类的个数等于1，转5，否则转3 5.画聚类图 6.决定类的个数和类。 判别分析：在已知研究对象分成若干类型，并已取得各种类型的一批已知样品的观测数据，在此基础上根据某些准则建立判别式，然后对未知类型的样品进行判别分类。 距离判别法—首先根据已知分类的数据，分别计算各类的重心，计算新个体到每类的距离，确定最短的距离（欧氏距离、马氏距离） Fisher判别法—利用已知类别个体的指标构造判别式（同类差别较小、不同类差别较大），按照判别式的值判断新个体的类别 Bayes判别法—计算新给样品属于各总体的条件概率，比较概率的大小，然后将新样品判归为来自概率最大的总体 模糊数学：研究和处理模糊性现象的数学（概念与其对立面之间没有一条明确的分界线）与模糊数学相关的问题：模糊分类问题—已知若干个相互之间不分明的模糊概念，需要判断某个确定事物用哪一个模糊概念来反映更合理准确；模糊相似选择—按某种性质对一组事物或对象排序是一类常见的问题，但是用来比 解析模拟退火算法 一.爬山算法(Hill Climbing) 介绍模拟退火前，先介绍爬山算法。爬山算法是一种简单的贪心搜索算法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到达到一个局部最优解。 爬山算法实现很简单，其主要缺点是会陷入局部最优解，而不一定能搜索到全局最优解。如图1所示：假设C点为当前解，爬山算法搜索到A点这个局部最优解就会停止搜索，因为在A点无论向那个方向小幅度移动都不能得到更优的解。 二.模拟退火(SA,Simulated Annealing)思想 爬山法是完完全全的贪心法，每次都鼠目寸光的选择一个当前最优解，因此只能搜索到局部的最优值。模拟退火其实也是一种贪心算法，但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以图1为例，模拟退火算法在搜索到局部最优解A后，会以一定的概率接受到E的移动。也许经过几次这样的不是局部最优的移动后会到达D点，于是就跳出了局部最大值A。 模拟退火算法描述： 若J(Y(i+1))>=J(Y(i))(即移动后得到更优解)，则总是接受该移动 若J(Y(i+1)) SARS传播的数学模型 (轩辕杨杰整理) 摘要 本文分析了题目所提供的早期SARS传播模型的合理性与实用性，认为该模型可以预测疫情发展的大致趋势，但是存在一定的不足.第一，混淆了累计患病人数与累计确诊人数的概念；第二，借助其他地区数据进行预测，后期预测结果不够准确；第三，模型的参数L、K的设定缺乏依据，具有一定的主观性. 针对早期模型的不足，在系统分析了SARS的传播机理后，把SARS的传播过程划分为：征兆期，爆发期，高峰期和衰退期4个阶段.将每个阶段影响SARS 传播的因素参数化，在传染病SIR模型的基础上，改进得到SARS传播模型.采用离散化的方法对本模型求数值解得到：北京SARS疫情的预测持续时间为106天，预测SARS患者累计2514人，与实际情况比较吻合. 应用SARS传播模型，对隔离时间及隔离措施强度的效果进行分析，得出结论：“早发现，早隔离”能有效减少累计患病人数；“严格隔离”能有效缩短疫情持续时间. 在建立模型的过程中发现，需要认清SARS传播机理，获得真实有效的数据.而题目所提供的累计确诊人数并不等于同期累计患病人数，这给模型的建立带来不小的困难. 本文分析了海外来京旅游人数受SARS的影响，建立时间序列半参数回归模型进行了预测，估算出SARS会对北京入境旅游业造成23.22亿元人民币损失，并预计北京海外旅游人数在10月以前能恢复正常. 最后给当地报刊写了一篇短文，介绍了建立传染病数学模型的重要性. 1．问题的重述 SARS （严重急性呼吸道综合症，俗称：非典型肺炎）的爆发和蔓延使我们认识到，定量地研究传染病的传播规律，为预测和控制传染病蔓延创造条件，具有很高的重要性.现需要做以下工作： （1） 对题目提供的一个早期模型，评价其合理性和实用性. （2） 建立自己的模型，说明优于早期模型的原因；说明怎样才能建立一个真正能够预测以及能为预防和控制提供可靠、足够信息的模型，并指出这样做的困难；评价卫生部门采取的措施，如：提前和延后5天采取严格的隔离措施，估计对疫情传播的影响. （3） 根据题目提供的数据建立相应的数学模型，预测SARS 对社会经济的影响. （4） 给当地报刊写一篇通俗短文，说明建立传染病数学模型的重要性. 2．早期模型的分析与评价 题目要求建立SARS 的传播模型，整个工作的关键是建立真正能够预测以及能为预防和控制提供可靠、足够的信息的模型.如何结合可靠、足够这两个要求评价一个模型的合理性和实用性，首先需要明确： 合理性定义 要求模型的建立有根据，预测结果切合实际. 实用性定义 要求模型能全面模拟真实情况，以量化指标指导实际. 所以合理的模型能为预防和控制提供可靠的信息；实用的模型能为预防和控制提供足够的信息. 2.1早期模型简述 早期模型是一个SARS 疫情分析及疫情走势预测的模型， 该模型假定初始时刻的病例数为0N ， 平均每病人每天可传染K 个人（K 一般为小数），K 代表某种社会环境下一个病人传染他人的平均概率，与全社会的警觉程度、政府和公众采取的各种措施有关.整个模型的K 值从开始到高峰期间保持不变，高峰期后 10天的范围内K 值逐步被调整到比较小的值，然后又保持不变. 平均每个病人可以直接感染他人的时间为L 天.整个模型的L 一直被定为20.则在L 天之内，病例数目的增长随时间t (单位天)的关系是： t k N t N )1()(0+?= 考虑传染期限L 的作用后，变化将显著偏离指数律，增长速度会放慢.采用半模拟循环计算的办法，把到达L 天的病例从可以引发直接传染的基数中去掉. 2.2早期模型合理性评价 根据早期模型对北京疫情的分析与预测，其先将北京的病例起点定在3月1日，经过大约59天在4月29日左右达到高峰，然后通过拟合起点和4月20日以后的数据定出高峰期以前的K =0.13913.高峰期后的K 值按香港情况变化，即10天范围内K 值逐步被调整到0.0273.L 恒为20.由此画出北京3月1日至5月7日疫情发展趋势拟合图像以及5月7日以后的疫情发展趋势预测图像，如图1. 模拟退火算法 模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis准则，粒子在温度T时趋于平衡的概率为e-ΔE/(kT)，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。 模拟退火算法的模型 模拟退火算法可以分解为解空间、目标函数和初始解三部分。 模拟退火的基本思想: (1)初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起 点)，每个T值的迭代次数L (2) 对k=1，……，L做第(3)至第6步： (3) 产生新解S′ (4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数 (5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T) 接受S′作为新的当前解. (6) 如果满足终止条件则输出当前解作为最优解，结束程序。终止条件通常取为连续若干个新解都没有被接受时终止算法。 (7) T逐渐减少，且T->0，然后转第2步。 算法对应动态演示图： 模拟退火算法新解的产生和接受可分为如下四个步骤： 第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。 第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。 第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropo1is准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。 第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则 Problem A Warmer Days or Sour Grapes ? The high quality of wines（葡萄酒）produced in the Finger Lakes Region（五指湖区）of upstate （北部）New York is widely known. Proximity（接近）to lakes tempers the climate and makes it more suitable for growing several varieties of premium（独特）grapes: R iesling（雷司令）, G ewürztraminer（琼瑶浆）, C hardonnay（霞多丽）, M erlot（梅洛）, P inot Noir（黑比诺）, and Cabernet F ranc（品丽珠）. (There are many more, but we will restrict（限制）the discussion to these six to simplify（简化）the modeling.) Each variety has its own preferred “average temperature” range but is also different in its susceptibility（感受性）to diseases and ability to withstand（抵抗）short periods of unusually cold temperature. As our local climate changes, the relative suitability of these varieties will be changing as well. A forward-looking winery（酒厂）has hired your team to help with the long-term planning. You will need to recommend a) the proportion（比例）of the total vineyard（葡萄园）to be used for growing each of the above six varieties; b) and when should these changes be implemented （实施）(based on observed temperatures and/or current market prices for each type of wine). Naturally, the winery is interested in maximizing its annual profit. But since the latter （后者）is weather-dependent, it might vary a lot year-to-year. You are also asked to evaluate the trade-offs （权衡）between optimizing the expected/average case versus the worst(-realistic-)scenario（情景）. Things to keep in mind: Climate modeling is complicated（复杂）and predicting the rate of “global warming” is a hotly debated area. For the purposes of this problem, assume that the annual average temperature in Ithaca（伊萨卡）, NY will increase by no more than 4°C by the end of this century. It is not all about the average temperature – a short snap（临时）of sub- zero（零度）temperature in late Ferburay or early March (after the vines already started getting used to warmer weather) is far more damaging than the same low temperature would be in the middle of the winter. It takes at least 3 years for a newly planted vine to start producing grapes suitable for winemaking. Problem B Outlook of Car-to-Car Tech SAN FRANCISCO -- After more than a decade of research into car-to-car communications, U.S. auto safety regulators took a step forward today by unveiling their plan for requiring cars to have wireless gear that will enable them to warn drivers of danger. 模拟退火算法简介与实例 2010-07-10 12:30:55| 分类：algorithms | 标签：|字号大中小订阅 摘要 模拟退火算法是S. Kirkpatrick, C. D. Gelatt和M. P. Vecchi在1983年所发明。是一种典型的概率模拟算法（Monte Carlo算法），其基本想想与冶金上的退火有相似之处，在一个相当大的空间内搜索最优解，而每次只搜索与自己临近的状态。此算法被证明以接近概率1接近最优解。其中有较好的物理思想，是模拟类算法中的典范。模拟退火算法由于要计算相临状态，这与Ising模拟的计算模拟有相似之处，因此本文也将对Ising做一个介绍。本文介绍算法的基本思想并做一个例子求解TSP问题（旅行商问题），重在介绍算法思想，具体算法的优化与改进不是本文涵盖范围。 1. Ising模型 Ising模型描述的是物体的铁磁性质，在铁和镍这类金属中，当温度低于居里温度时，原子的自旋自发地倾向某个方向，而产生宏观磁矩。温度高于居里温度时,自旋的取向非常紊乱,因而不产生净磁矩。当温度从大于或小于两边趋于居里温度时，金属的比热容趋于无限大。这是物质在铁磁性状态和非铁磁性状态之间的相变。伊辛模型就是模拟铁磁性物质的结构，解释这类相变现象的一种粗略的模型。它的优点在于，用统计物理方法，对二维情形求得了数学上严格的解。这就使得铁磁性物质相变的大致特征，获得了理论上的描述。 1.1模型描述 这个模型所研究的系统是由N个阵点排列成n维周期性点阵，这里n=2。点阵的几何构形可以是立方的或六角形的,每个阵点上都赋予一个取值+1或-1的自旋变量i,如果i＝+1,即第N个阵点的自旋向上；如i＝-1，即第个N阵点的自旋向下并且认为只是最近邻的自旋之间有相互作用。点阵的位形用一组自旋变量（这里i＝2）来确定，如下图所示 图1，模型图示图2，最近临磁子 1.2模型计算 1)两个相临磁子趋向平行能量最低，即两个磁子的自旋方向非平行与平行。能量相差ΔE。 2)每个磁子的磁矩为m,总的磁矩为每个磁子的磁矩和。模拟退火算法原理及matlab源代码
模拟退火算法(MATLAB实现)
商品期货交易数学建模
数学建模 模拟退火
爬山算法、模拟退火算法、遗传算法
数学建模模拟试题
数学建模常用算法模型
数学建模中常见的十大模型
数学建模方法归类(很全很有用)
模拟退火算法介绍
SARS传播的数学模型 数学建模全国赛优秀论文
模拟退火算法算法的简介及程序
2015年美国数学建模竞赛第二次模拟赛题
模拟退火算法简介与实例
相关主题
文本预览