模擬退火法含程序(數(shù)模講稿)

一、模擬退火法模擬退火法（參見［1,2］）作為一種適合于求解大規(guī)模的優(yōu)化問題的技術(shù)，近來已引起極大的關(guān)注。特別是當(dāng)優(yōu)化問題有很多局部極值而全局極值又很難求出時，模擬退火法尤其有效。在實用上，它有效地“解決了”著名的旅行推梢員問題，即在必須依次訪問每一個城市（共有N個城市）的前提下，為旅行推銷員設(shè)計一條能夠返回起點的最短旅程。模擬退火方法還被成功地用于設(shè)計復(fù)雜的集成電路，也就是說如何最佳地安排幾十萬個電路元件，使它們?nèi)考稍谝粋€很小的硅片上，而相互連接的線路之間的纏繞能夠達到最?。▍⒁姡?,4］）.盡管模擬退火法的功效非凡，但它的算法實現(xiàn)卻相對地簡單，這一點似乎有些不可思議。請注意，我們上面提到的兩個例子都屬于組合極小化問題?，F(xiàn)本類問題通常也有一個目標(biāo)函數(shù)，但是函數(shù)的定義域并不是簡單地由N個連續(xù)參變量組成的N維空間，而是一個離散的巨大空間，例如，由所有可能的城市旅行路線組成的集合，或者硅片電路元件的所有可能的分配方式的集合。構(gòu)形空間中元素的數(shù)量相當(dāng)巨大，根本不可能窮舉，而且因為集合是離散的，我們也不可能“沿合適的方向連續(xù)下降”。因此在構(gòu)形空間中，“方向”概念就沒有什么意義了。后面我們還將介紹如何在其有連續(xù)控制參數(shù)的空間中利用模擬退火法。這種應(yīng)用實際上要比組合問題復(fù)雜一些，因為其中又要出現(xiàn)“狹長山谷”的情況。正如在下文中我們將看到的，模擬退火法的試探步驟是“隨機”的。但在一個狹窄且漫長的等高線山谷中，幾乎所有的隨機步驟都呈向上的趨勢，因此，算法中需要增加一些技巧。模擬退火的核心思想與熱力學(xué)的原理頗為相似，而且尤其類似于液體流動和結(jié)晶以及金屬冷卻和退火方式。在高溫下，一種液體的大最分子彼此之間進行著相對自由的移動。如果該流體慢慢地冷卻下來，熱能可動性便會消失。大量原子常常能夠自行排列成行，形成一個純凈的晶體，該晶體在各個方向上都被完全有序地排列在幾百萬倍于單個原子大小的距離之內(nèi)。對于這個系統(tǒng)來說晶體狀態(tài)是能量最低狀態(tài)；而所有緩慢冷卻的系統(tǒng)都可以自然達到這個能量最低狀態(tài)，這可以說是一個令人驚奇的事實。實際上，如果某種液體金屬被迅速冷卻或被“猝熄”，那么它不會達到這一狀態(tài)，而只能達到一種具有較高能量的多晶狀態(tài)或非結(jié)晶狀態(tài)。因此，這一過程的本質(zhì)在于緩緩地致冷，以爭取充足的時間，讓大量原子在喪失可動性之前進行重新分布。這就是所謂退火在技術(shù)上的定義，同時也是確保達到低能量狀態(tài)所必需的條件。盡管我們的比喻并不算貼切，但是迄今為止本身所討論的所有極小化算法確實與快速冷卻猝熄有某種關(guān)聯(lián)之處。以往我們處理問題的方式都是:從初始點開始，立即沿下降方向前進，走得越遠越好，似乎這樣才能迅速求得問題的解但是，正如前面常常提到的，這種方法往往只能求得局部極小點，卻求不到整體最小點。自然界本身的極小化算法則基于一種截然不同的方式，所謂的玻爾茲曼（Boltzmann）概率分布ProbE~ExpE/kT （1）表達了這樣一種思想，即:一個處于熱平衡狀態(tài)且具有溫度T的系統(tǒng)，其能量按照概率，分布于所有不同能量狀態(tài)e之中。即使在很低的溫度下，系統(tǒng)也有可能（雖然這種可能性很?。┨幱谝粋€較高的能量狀態(tài)。因此，相應(yīng)地，系統(tǒng)也能夠獲得擺脫局部能量極小點的機會。并找到一個更好的、更接近于整體的極小點。式（1）中的參數(shù)k（稱為玻爾茲曼常數(shù)）是一個自然常數(shù)，它的作用是將溫度與能量聯(lián)系起來。換句話說，在有些情況下系統(tǒng)的能量可上升，也可下降，但是溫度越低，顯著上升的可能性就越小。1953年，米特羅波利斯（Metropolis）及其合作者們首次將這種原理滲透到數(shù)值計算中。他們對一個模擬熱力學(xué)系統(tǒng)提供了一系列選擇項，并假設(shè)：系統(tǒng)構(gòu)形從能量e變化到能量E的概率為p=ExpL（E-EkT]o很顯然，如果E<E，1221’21p將大于1;在這類情況下，對構(gòu)形的能量變化任意指定一個概率值p=1，也就是說，該系統(tǒng)總是取這個選擇項。這種格式總是采取下降過程，偶爾采取上升步驟。目前已被公認(rèn)為米特羅波利斯算法。為了將米特波利斯算法應(yīng)用于熱力學(xué)以外的系統(tǒng)，必須提供以下幾項基本要素：對可能的系統(tǒng)構(gòu)形的一種描述。一個有關(guān)構(gòu)形內(nèi)部隨機變化的生成函數(shù)，這些變化將作為“選擇項”提交給該系統(tǒng)。一個目標(biāo)函數(shù)E（類似于能量），求解E的極小值，即為算法所要完成的工作。一個控制參數(shù)T（類似于溫度）和一個退火進程，該進程用來說明系統(tǒng)是如何從高值向低值降低的，例如在溫度T時每次下降步驟中要經(jīng)過多少次隨機的構(gòu)形變化以及該步長是多大等等。應(yīng)說明的是，這里“高”和“低”的含義，還有進程表的確定，都需要一定的物理知識和/或一些摸索的實驗。模擬退火算法的模型模擬退火算法可以分解為解空間、目標(biāo)函數(shù)和初始解三部分。模擬退火的基本思想:（1）初始化：初始溫度T（充分大），初始解狀態(tài)S（是算法迭代的起點），每個T值的迭代次數(shù)L對k=l,……,L做第(3)至第6步：產(chǎn)生新解9計算增量氐?。?')-。?),其中C(S)為評價函數(shù)若MvO則接受S'作為新的當(dāng)前解，否則以概率exp(-M/T)接受S，作為新的當(dāng)前解.如果滿足終止條件則輸出當(dāng)前解作為最優(yōu)解,結(jié)束程序。終止條件通常取為連續(xù)若干個新解都沒有被接受時終止算法。T逐漸減少，且T->0,然后轉(zhuǎn)第2步。模擬退火算法新解的產(chǎn)生和接受可分為如下四個步驟：第一步，是由一個產(chǎn)生函數(shù)從當(dāng)前解產(chǎn)生一個位于解空間的新解；為便于后續(xù)的計算和接受，減少算法耗時，通常選擇由當(dāng)前新解經(jīng)過簡單地變換即可產(chǎn)生新解的方法，如對構(gòu)成新解的全部或部分元素進行置換、互換等，注意到產(chǎn)生新解的變換方法決定了當(dāng)前新解的鄰域結(jié)構(gòu)，因而對冷卻進度表的選取有一定的影響。第二步，是計算與新解所對應(yīng)的目標(biāo)函數(shù)差。因為目標(biāo)函數(shù)差僅由變換部分產(chǎn)生，所以目標(biāo)函數(shù)差的計算最好按增量計算。事實表明，對大多數(shù)應(yīng)用而言，這是計算目標(biāo)函數(shù)差的最快方法。第三步，是判斷新解是否被接受,判斷的依據(jù)是一個接受準(zhǔn)則，最常用的接受準(zhǔn)則是Metropolis準(zhǔn)則：若At'vO則接受S'作為新的當(dāng)前解S,否則以概率exp(-At'/T)接受S'作為新的當(dāng)前解So第四步,是當(dāng)新解被確定接受時,用新解代替當(dāng)前解,這只需將當(dāng)前解中對應(yīng)于產(chǎn)生新解時的變換部分予以實現(xiàn)，同時修正目標(biāo)函數(shù)值即可。此時，當(dāng)前解實現(xiàn)了一次迭代?？稍诖嘶A(chǔ)上開始下一輪試驗。而當(dāng)新解被判定為舍棄

時，則在原當(dāng)前解的基礎(chǔ)上繼續(xù)下一輪試驗。模擬退火算法與初始值無關(guān)，算法求得的解與初始解狀態(tài)S（是算法迭代的起點）無關(guān)；模擬退火算法具有漸近收斂性，已在理論上被證明是一種以概率1收斂于全局最優(yōu)解的全局優(yōu)化算法；模擬退火算法具有并行性。模擬退火算法的簡單應(yīng)用作為模擬退火算法應(yīng)用，討論貨郎擔(dān)問題（TravellingSalesmanProblem，簡記為TSP）:設(shè)有n個城市，用數(shù)碼l,...,n代表。城市i和城市j之間的距離為d（i，j）i,j=l,...,n.TSP問題是要找遍訪每個域市恰好一次的一條回路，且其路徑總長度為最短.。求解TSP的模擬退火算法模型可描述如下:解空間:解空間S是遍訪每個城市恰好一次的所有路經(jīng)， 是｛1，……n｝的所有循環(huán)排列的集合,S中的成員記為h,w,…w｝，w,w,…w是1,2,……,n12n12n的一個排列，表明從w城市出發(fā)，依次經(jīng)過w,w,…w城市，再返回w城123n1市。初始解可選為（1,……,n）;目標(biāo)函數(shù):目標(biāo)函數(shù)為訪問所有城市的路徑總長度 或稱為代價函數(shù);我們要求的最優(yōu)路徑為目標(biāo)函數(shù)（代價函數(shù)）為最小值時對應(yīng)的路徑。新路徑的產(chǎn)生:隨機產(chǎn)生1新路徑的產(chǎn)生:隨機產(chǎn)生1和n之間的兩相異數(shù)k和m，不妨假設(shè)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對應(yīng)的兩個城市在路徑序列中交換位置稱為2-opt映射。

根據(jù)上述描述，模擬退火算法求解TSP問題的流程框圖如下圖2模擬退火算法的流程框圖也可以采用其他的變換方法，有些變換有獨特的優(yōu)越性，有時也將它們交替使用，得到一種更好方法。例如：隨機產(chǎn)生1和n之間的兩相異數(shù)k和m，不妨假設(shè)kvm，則將原路徑Cw,w,…w,w,…w,w…w1 2 k k+1 m m+1, n變?yōu)樾侣窂剑篊w,w,…w,w,…w,w,w…w1 2 m m-1 k+1 k, k+2 n上述變換方法就是將k和m之間對應(yīng)的多個城市在路徑序列顛倒位置,可簡單說成是“逆轉(zhuǎn)中間或者逆轉(zhuǎn)兩端”。根據(jù)上述分析，可寫出用模擬退火算法求解TSP問題的偽程序：ProcedureTSPSA:begininit-of-T; //T為初始溫度S={1，……，n}; 〃S為初始值termination=false;whiletermination=falsebeginfori=1toLdobegingenerate(SformS); 〃從當(dāng)前回路S產(chǎn)生新回路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模擬退火算法的應(yīng)用很廣泛，可以較高的效率求解最大截問題(MaxCutProblem)、0-1背包問題(ZeroOneKnapsackProblem)、圖著色問題(GraphColouringProblem)、調(diào)度問題(SchedulingProblem)等等。模擬退火算法的參數(shù)控制問題模擬退火算法的應(yīng)用很廣泛，可以求解NP完全問題，但其參數(shù)難以控制，其主要問題有以下三點：溫度T的初始值設(shè)置問題。溫度T的初始值設(shè)置是影響模擬退火算法全局搜索性能的重要因素之一、初始溫度高，則搜索到全局最優(yōu)解的可能性大，但因此要花費大量的計算時間；反之，則可節(jié)約計算時間，但全局搜索性能可能受到影響。實際應(yīng)用過程中，初始溫度一般需要依據(jù)實驗結(jié)果進行若干次調(diào)整。退火速度問題。模擬退火算法的全局搜索性能也與退火速度密切相關(guān)。一般來說，同一溫度下的“充分”搜索(退火)是相當(dāng)必要的，但這需要計算時間。實際應(yīng)用中，要針對具體問題的性質(zhì)和特征設(shè)置合理的退火平衡條件。溫度管理問題。溫度管理問題也是模擬退火算法難以處理的問題之一。實際應(yīng)用中，由于必須考慮計算復(fù)雜度的切實可行性等問題，常采用如下所示的降溫方式：T(t+1)=kxT(t)式中k為正的略小于1.00的常數(shù)，t為降溫的次數(shù)二、有記憶的模擬退火法傳統(tǒng)的模擬退火算法思路清晰、原理簡單、可用于解決許多實際問題。但存在很多弊病，國內(nèi)外學(xué)者并對其作相應(yīng)的改進，得到了改進的模擬退火算法，包括：加溫退火法、有記憶的模擬退火算法、帶返回搜索的模擬退火算法、多次尋優(yōu)法、回火退火算法等。針對特征選擇的特點，本文選用了有記憶的模擬退火算法作為搜索策略。在傳統(tǒng)的模擬退火過程中，算法終止于一個預(yù)先規(guī)定的停止準(zhǔn)則S,如：控制參數(shù)t的值小于某個充分小的正數(shù)；相繼的若干個Markov鏈中解未得到任何改善；兩個相繼Markov鏈所得解的差的絕對值小于某個充分小的正數(shù)等。但由于模擬退火算法的搜索過程是隨機的，且當(dāng)t值較大時可以接受部分惡化解，而隨著t值的衰減，惡化解被接受的概率逐漸減小直至為0。另一方面，某些當(dāng)前解要達到最優(yōu)解時必須經(jīng)過暫時惡化的“山脊”，因此，上述這些停止準(zhǔn)則無法保證最終解正好是整個搜索過程中曾經(jīng)達到的最優(yōu)解。因此，在傳統(tǒng)的算法中增加一個記憶器，使之能夠記住搜索過程中遇到過的最好解，當(dāng)退火結(jié)束時，將所得最終解與記憶器中的解比較并取較優(yōu)者作為最后結(jié)果。三、組合極小化：旅行推銷員問題

下面是我們用“旅行推銷員問題”為具體實例說明模擬退火法的應(yīng)用。假

設(shè)一個推銷員，要去N個分別位于(x,y)的城市進行推銷，并于最后返回他原來

ii所在的城市，要求每個城市只能去一次，而且所經(jīng)過的路徑要盡可能地短。這個間題屬于一類所謂“NP-完全問題”這類問題求出一個精確解所需的計算時間是隨N的增加以指數(shù)exp(常數(shù)XN)增長的。當(dāng)N不斷增大時，運行時間將迅速增加，進而導(dǎo)致費用高到令人難以接受的程度。旅行推銷員問題也是眾多極小化問題中的一種，它的目標(biāo)函數(shù)e具有多個局部極小值。在實際應(yīng)用當(dāng)中，常常有足夠多的條件可以從多個極小值中選出一個最小的，這個最小值即使不是絕對最小，也相當(dāng)接近絕對最小了。退火法的目的就是要獲得這個最小值，同時又要將計算量限制在N的低階次的數(shù)量級上。旅行推銷員問題也是按模擬退火問題的方式進行處理的。具體如下：1?構(gòu)形將N個城市分別標(biāo)記為i=1,2,N中的數(shù)，其中每個城市具有坐標(biāo)一個構(gòu)形就是數(shù)字1N的一個排列，可以解釋為推銷員途徑的城市的順序。2?調(diào)整林(Lin)曾經(jīng)提出過一種所謂“轉(zhuǎn)移有效集”這里的“轉(zhuǎn)移”包括兩種類型：(a)移走路徑的某一段，然后對這段路徑上的城市用相反次序重新進行排列，并用后者來代替前者。(b)移走某段路徑，并用位于城市間的隨機選取的另一段路徑來取代被移走的路徑。3.目標(biāo)函數(shù)在旅行推銷員問題的一種最簡單的形式中，目標(biāo)函數(shù)e被定義為旅途的總長度E二L三迓* —》+(y_y—F (10.9.2)Xii+1 ii+1i=1這里認(rèn)為第N+1個點與第1個點是重合的。但是為了表明模擬退火法的靈活性，我們還要用到下面的技巧：假設(shè)推銷員無端端地害怕飛越密西西比河，在這種情況下，我們對每個城市給定一個參數(shù)卩，如果該城市位于密西西比河i以東，卩取1；若在密西西比河以西則取-1。對于目標(biāo)函數(shù)E，我們將其改寫為:iE=為打 _J2+(y_y_X+九Q-卩》 (10.9.3)L'ii+1 ii+1 ii+1_i=1由于每過一次河都將以4九作為懲罰，因而現(xiàn)在我們設(shè)計的算法的目標(biāo)，就變成了尋找盡可能回避過河的最短路徑。路徑長度對過河次數(shù)的相對重要性將由我們選擇的九來確定。圖10.9.1表明了所得的結(jié)果。顯然，這種技巧可以推廣到包含許多相互沖突的目的要求的極小化問題當(dāng)中。

0 .5(C)圖（A）表示的是從四個隨機分布的城市中間找到的一條（接近）最短路徑，中間豎虛線標(biāo)識的是一條河流，但這是對過河沒有附加您罰項的情況。在圖（B）中對過河施加的懲罰項很大，而圖中所示的解本身的過河次數(shù)也相應(yīng)地只有少得不能再少的兩次。在圖（C）中懲罰項為負(fù)，這就是說，推銷員實際上成了恣意偷渡的走私者！圖1用模擬退火法解決旅行推銷員問題4.退火進程這一步需要借助試驗來確定。首先要進行一些隨機調(diào)整，然后利用它們來確定從轉(zhuǎn)移到轉(zhuǎn)移過程中將會遇到的AE值之范圍。對參數(shù)T取一初始值（這個初始值要遠遠大于通常所能遇到的AE的最大值），并以倍增的步長下減，每次使T總共減少10%。我們拿每個新的常數(shù)T值去試各種100N重構(gòu)形,或1ON成功的重構(gòu)形，無論哪個在前出現(xiàn)就取哪個。當(dāng)實在不能再進一步減小E時，則停止。F面的旅行推銷員程序利用了米特羅波利斯(Metropolis)算法，并展示了模擬退火技術(shù)應(yīng)用于組合問題的幾個主要方面。#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個城市的位置坐標(biāo)存貯在數(shù)組x[l..ncity]和y[l..ncity]中。數(shù)組iorder[l..ncity]表示途徑城市的順序。在輸出項中，iorder中的元素將被置為數(shù)字1到ncity的某排對，本程序?qū)⒎祷厮芮蟪龅淖罴堰x擇路徑。*/{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; //在任何溫度下，所試路徑的最大次數(shù)nlimit=10*ncity; //在繼續(xù)進行之前，成功的路徑改變的最大次數(shù)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]);

//將路徑頭尾相連并結(jié)束循環(huán)i1=iorder[ncity]; ////將路徑頭尾相連并結(jié)束循環(huán)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));//?段的結(jié)尾if(n[2]>=n[1])++n[2];nn=1+((n[1]-n[2]+ncity-1)%ncity);//nn為不位于當(dāng)前段上的城市數(shù)}while(nn<3);idec=irbit1(&iseed);//確定是否做段反轉(zhuǎn)或段輸送if(idec==0){ //做輸送n[3]=n[2]+(int)(abs(nn-2)*ran3(&idum))+1;n[3]=1+((n[3]-1)%ncity);//輸送到一個不在當(dāng)前路徑上的某處de=trncst(x,y,iorder,ncity,n); //計算代價ans=metrop(de,t); //做預(yù)測if(ans){++nsucc;path+=de;//輸送工作結(jié)束//做段反轉(zhuǎn)////輸送工作結(jié)束//做段反轉(zhuǎn)//計算代價//作預(yù)測//完成段反轉(zhuǎ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; //如果成功的轉(zhuǎn)移超過次數(shù)則提前結(jié)束}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[])/*該子程序返回的是反轉(zhuǎn)某給定路徑所需的代價函數(shù)值。參數(shù)中ncity為城市數(shù)；數(shù)組x[l..ncity],y[l..ncity]為這些城市的位置坐標(biāo)；iorder[n..ncity]為當(dāng)前路線；數(shù)組n的頭兩個元素n[1]和n[2]分別代表將要被反轉(zhuǎn)的路徑段的起點城市和終點城市。子程序的輸出項de為反轉(zhuǎn)所需的代價值,但真正的反轉(zhuǎn)過程并不是由該程序執(zhí)行。*/{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]]; //求四個所涉及到的城市的坐標(biāo)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[])/*該子程序的作用是執(zhí)行一段路徑的反轉(zhuǎn)過程。輸入?yún)?shù)iorder[l..ncity]給出當(dāng)前的路線順序；向量n的前四個元素中，n⑴和n[2]分別為將要被反轉(zhuǎn)的路徑的起點和終點城市，n⑶和n[4]則分別為緊隨n[1]之前和緊隨n[2]之后的兩個城市的標(biāo)號，其中n[3]和n[4]由函數(shù)revest給出。在輸出端，iorder又將被作為返回值，它的定義是n[1]到n[2]路段被反轉(zhuǎn)后的行程路線。*/{intnn,j,k,l,itmp;nn=(1+((n[2]-n[1]+ncity)%ncity))/2;//為實現(xiàn)反轉(zhuǎn)操作，必須交換這么多城市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[])/*該子程序返回的是輸送某段給定路徑所需的代價函數(shù)值。輸入?yún)?shù)中ncity為城市的個數(shù)；x[l..ncity]和y[l..ncity]為這些城市的位置坐標(biāo)，數(shù)組n的前三個元索分別為：將要被輸送的路徑段的起、止點城市以及這段路徑將要被插入處的標(biāo)志城市(插在該城市之后)，該子程序的輸出項de為計算出來的輸送代價值，但實際的輸送操作井不由該子程序執(zhí)行。*/{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]]; //求出六個城市的有關(guān)坐標(biāo)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[])/*該子程序的作用是執(zhí)行真正的段輸送操作，輸入?yún)?shù)iorder[l..ncity]給出當(dāng)前的路徑順序；數(shù)組n共有6個元素，其意義分別為：n⑴和n[2]分別代表將要被輸送的路徑段的起點城市和終點城市;n[3]和n[4]為兩個相鄰的城市,n[1]至n[2]路徑段即將放入它們中間;n⑸和n⑹分別為n[1]之前和n[2]之后的兩個城市。在這六個元索中n[4],n⑸和n[6]由函數(shù)trncst給出。在輸出端，iorder將根據(jù)路徑段的移動和變化作出相應(yīng)的修改。*/{intml,m2,m3,nn,j,jj,*jorder;jorder=ivector(1,ncity);m1=1+((n[2]-n[1]+ncity)%ncity);//找出位于n[1]和n[2]間城市個數(shù)m2=1+((n[5]-n[4]+ncity)%ncity);〃?n[4]到n[5]間的城市個數(shù)m3=1+((n[3]-n[6]+ncity)%ncity);//??n⑹和n[3]間的城市個數(shù)nn=l;for(j=1;j<=m1;j++){jj=1+((j+n[1]-2)%ncity); //復(fù)制所選路徑段jorder[nn++]=iorder[jj];}if(m2>0){for(j=1;jv=m2;j++){ 〃復(fù)制n[4倒n⑸間的路徑段jj=1+((j+n[4]-2)%ncity);jorder[nn++]=iorder[jj];}if(m3>0){for(j=1;jv=m3;j++){ 〃最后，復(fù)制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)/*該子程序為米特羅波利斯算法的程序?qū)崿F(xiàn)。metrop返回的是一個布爾型變量，由該變量決定是否接受一個使目標(biāo)函數(shù)e產(chǎn)生改變量de的重構(gòu)形。如果de<0則metrop=l(真)，而當(dāng)de>O時，metrop為真的概率是exp(—de/t)，這里才是一個由退火進程決定的溫度值。*/{floatran3(long*idum);staticlonggljdum=1;returndev0.0||ran3(&gljdum)vexp(-de/t);}模擬退火法的基本思想也可以應(yīng)用于具有連續(xù)N維控制參數(shù)的極小化問題當(dāng)中，例如，在某個函數(shù)f6）（這里x為一個N維向量）有許多局部極小點的情況下，求解它的（理想的全局的）極小值。這時米特羅波利斯算法所需的四要素可以具體化為：1） 目標(biāo)函數(shù)即為f（x）的值；2） 系統(tǒng)狀態(tài)描態(tài)即為N維空間中的點x；3）類似于溫度的控制參數(shù)T以及一個使T逐漸降低的退火進程仍為原先的定義；4）描述構(gòu)形內(nèi)部隨機變化的發(fā)生器即為一個從x到x+Ax采取隨機步驟的方法。在上述四要素中問題最大的是最后一條。目前已發(fā)表的文獻中［7-10］，介紹了幾種不同的選擇辦法。但我們認(rèn)為，這些方法都不算成功。問題在于“效率”二字：當(dāng)局部的向下運動存在時，如果某個隨機變化發(fā)生器幾乎總是作出向上運動的決策，那么它的效率就很低。我們認(rèn)為，一個好的發(fā)生器在等高線的“窄谷”中仍應(yīng)保持高效性；當(dāng)算法在接近收斂到極小點處時，它的效率也不應(yīng)越變越低。在上面我們提到的幾篇文獻中，除［7］中介紹的方法之外，其他所有方法都表現(xiàn)不同程度的低效性。下面我們將要介紹的這種方法，利用了下降單純形法的一種修改后的形式。首先我們將米特羅波利斯算法四要素中的系統(tǒng)狀態(tài)描述，由單個點x改為一個其有N+1個點的單純形；單純形的“操作”分為反射、擴張和收縮三種。然后我們將一個正的、呈對數(shù)分布的隨機變量（與溫度T成比例）添加到存貯的函數(shù)值中（該函數(shù)值與單純形的每個頂點都有關(guān)聯(lián)），再從每個被當(dāng)作替代的新點的函數(shù)值中減去一個類似的隨機變量。和普通的米特羅波利斯方法一樣，這種方法總是接受一個真正的下降步驟。但有時也接受一次上升步驟。在極限過程T-0中，該算法恰好變成了下降的單純形法，并收斂到一個局部極小點。在T的某有限值處，單純形將擴展到一定的規(guī)模，其大小接近于在這個溫度值所能達到的區(qū)域；然后單純形在這個區(qū)域內(nèi)部做隨機的滾動翻轉(zhuǎn)式布朗運動，并在該過程中抽取一些新的、近似隨機的點作為樣本。一個區(qū)域被利用的效率與其狹窄度（對橢球狀山谷，指它的主軸比率）及方向性均無關(guān)。如果溫度降低得足夠緩慢，那么單純形將極有可能收縮到那個區(qū)域內(nèi)，而那個區(qū)域內(nèi)包含已遇到的最低相對極小。由于在所有模擬退火法的應(yīng)用場合中，“足夠緩慢”一詞的含意根據(jù)問題的不同可以有相當(dāng)大的細(xì)微區(qū)別，因而成功與否往往取決于退火進程選擇。下面幾種可能性我們認(rèn)為值一試：?每經(jīng)過m步移動之后，將T減到（1_e）T,這里。*/m的具體值要通過實驗確定。?設(shè)總的移動步數(shù)為K,每經(jīng)過m步移動之后將T減到t二TG-kT、，其0'中k是到目前為止所經(jīng)過的步數(shù)的累加值。d為一常數(shù)，可取為1，2:4等。（X的最佳值取決于各種深度的相對于極小的統(tǒng)計分布，稍大一些的X值在較低溫度時，需化費的迭代次數(shù)將更多；?每經(jīng)過m步驟移動之后，置T為0乘f-f，其中0為一個階數(shù)為1的常1b數(shù)，其具體值由試驗確定；f為目前單純形中最小的函數(shù)值；f為曾經(jīng)遇到的最1b佳函數(shù)值。但應(yīng)注意的是，T的降低幅度一次不要超過某個分?jǐn)?shù)值丫。另一個策略方法上的問題是，當(dāng)單純形的某個頂點被放棄并讓位于“永遠的最佳點”時，是否需要采取重新開始的步驟（但我們必須保證在進行這項工作時，這個“永遠的最佳點”當(dāng)前并不在單純形中!）。對于有些問題重新開始（例如，只要溫度降低了因子3即執(zhí)行重新開始步驟）的效果極佳，而對于另外一些問題，重新開始不僅沒有任何效果反而會產(chǎn)生負(fù)作用。上述兩種截然相反的情況，我們都找到了例子可以作為例證。將下面的程序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的下降單純形法相結(jié)合的方法求多元函數(shù)funk(x)的極小值。其中x[1..ndim]為一個ndim維向量。作為輸入?yún)?shù)的矩陣p[1..ndim][1..ndim〕共有ndim+1行，每行均為一個ndim維向量。分別代表初始單純形的各個頂點。amebsa的輸入項還包括矢量y[1..ndim+1]、浮點效ftol以及iter和temptr。其中y的各個元素將被初始化為函數(shù)funk在P的ndim+1個頂點(即行)的值；ftol為函數(shù)值所要達到的相對收斂容許限，一旦滿足要求，程序?qū)⒈M早返回;iter和temptr的意義分別是函數(shù)值的計算次數(shù)和溫度。程序在某退火溫度temptr處進行iter次函數(shù)值計算后返回。而接下來要做的事就是根據(jù)退火進程調(diào)低溫度temptr；重置iter并再次調(diào)用該程序(每次調(diào)用時其他參數(shù)保持不變)。如果iter以正值返回，則說明程序正常收斂。如果第一次調(diào)用時yb被初始化為一個很大的數(shù)，則yb和pb[1..ndim]將依次返回已遇到的最佳函數(shù)值和最佳點(這個最佳點可以水遠不是單純形中的點)。*/{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));//我們所“看到”的頂點總是處于//隨機的熱起伏狀態(tài)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++){ //對單純形中的點進行循環(huán)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){//得到了比最佳點還要好的結(jié)果，因此用因子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//糾正計數(shù)器//重新計算psum//糾正計數(shù)器//在主程序中定義和初始賦值//在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));//我們曾經(jīng)對所有的當(dāng)前頂點添加熱起伏運動，if(yflu<*yhi){