一、粒子群演算法簡介

1 引言
模擬退火演算法(Simulated Annealing，SA)的思想最早由Metropolis等人於1953年提出：Kirkpatrick於1983年第一次使用模擬退火演算法求解組合最優化問題[1] 。模擬退火演算法是一種基於MonteCarlo迭代求解策略的隨機尋優演算法， 其出發點是基於物理中固體物質的退火過程與一般組合優化問題之間的相似性。其目的在於為具有NP(Non-deterministic Polynomial) 複雜性的問題提供有效的近似求解演算法，它克服了其他優化過程容易陷入區域性極小的缺陷和對初值的依賴性。
模擬退火演算法是一種通用的優化演算法，是區域性搜尋演算法的擴充套件。它不同於區域性搜尋演算法之處是以一定的概率選擇鄰域中目標值大的劣質解。從理論上說，它是一種全域性最優演算法。模擬退火演算法以優化問
題的求解與物理系統退火過程的相似性為基礎， 它利用Metropolis演算法並適當地控制溫度的下降過程來實現模擬退火，從而達到求解全域性優化問題的目的[2].
模擬退火演算法是一種能應用到求最小值問題的優化過程。在此過程中，每一步更新過程的長度都與相應的引數成正比，這些引數扮演著溫度的角色。與金屬退火原理相類似，在開始階段為了更快地最小
化，溫度被升得很高，然後才慢慢降溫以求穩定。
目前，模擬退火演算法迎來了興盛時期，無論是理論研究還是應用研究都成了十分熱門的課題[3-7]。尤其是它的應用研究顯得格外活躍，已在工程中得到了廣泛應用，諸如生產排程、控制工程、機器學習、神經網路、模式識別、影象處理、離散/連續變數的結構優化問題等領域。它能有效地求解常規優化方法難以解決的組合優化問題和複雜函式優化問題，適用範圍極廣。
模擬退火演算法具有十分強大的全域性搜尋效能，這是因為比起普通的優化搜方法，它採用了許多獨特的方法和技術：在模擬退火演算法中，基本不用搜索空間的知識或者其他的輔助資訊，而只是定義鄰域
結構，在其鄰域結構內選取相鄰解，再利用目標函式進行評估：模擬退火演算法不是採用確定性規則，而是採用概率的變遷來指導它的搜尋方向，它所採用的概率僅僅是作為一種工具來引導其搜尋過程朝著更優化解的區域移動。因此，雖然看起來它是一種盲目的搜尋方法，但實際上有著明確的搜尋方向。

2 模擬退火演算法理論
模擬退火演算法以優化問題求解過程與物理退火過程之間的相似性為基礎，優化的目標函式相當於金屬的內能，優化問題的自變數組合狀態空間相當於金屬的內能狀態空間，問題的求解過程就是找一個組
合狀態， 使目標函式值最小。利用Me to polis演算法並適當地控制溫度的下降過程實現模擬退火，從而達到求解全域性優化問題的目的[8].
2.1物理退火過程
模擬退火演算法的核心思想與熱力學的原理極為類似。在高溫下，液體的大量分子彼此之間進行著相對自由移動。如果該液體慢慢地冷卻，熱能原子可動性就會消失。大量原子常常能夠自行排列成行，形成一個純淨的晶體，該晶體在各個方向上都被完全有序地排列在幾百萬倍於單個原子的距離之內。對於這個系統來說，晶體狀態是能量最低狀態，而所有緩慢冷卻的系統都可以自然達到這個最低能量狀態。實際上，如果液態金屬被迅速冷卻，則它不會達到這一狀態，而只能達到一種只有較高能量的多晶體狀態或非結晶狀態。因此，這一過程的本質在於緩慢地冷卻，以爭取足夠的時間，讓大量原子在喪失可動性之前進行重新分佈，這是確保能量達到低能量狀態所必需的條件。簡單而言，物理退火過程由以下幾部分組成：加溫過程、等溫過程和冷卻過程。

加溫過程
其目的是增強粒子的熱運動，使其偏離平衡位置。當溫度足夠高時，固體將熔解為液體，從而消除系統原先可能存在的非均勻態，使隨後進行的冷卻過程以某一平衡態為起點。熔解過程與系統的能量增
大過程相聯絡，系統能量也隨溫度的升高而增大。

等溫過程
通過物理學的知識得知，對於與周圍環境交換熱量而溫度不變的封閉系統，系統狀態的自發變化總是朝著自由能減小的方向進行：當自由能達到最小時，系統達到平衡態。

冷卻過程
其目的是使粒子的熱運動減弱並逐漸趨於有序，系統能量逐漸下降，從而得到低能量的晶體結構。

2.2 模擬退火原理
模擬退火演算法來源於固體退火原理，將固體加溫至充分高，再讓其徐徐冷卻。加溫時，固體內部粒子隨溫升變為無序狀，內能增大：而徐徐冷卻時粒子漸趨有序，在每個溫度上都達到平衡態，最後在常
溫時達到基態，內能減為最小。模擬退火演算法與金屬退火過程的相似關係如表7.1所示。根Metropolis準則， 粒子在溫度7時趨於平衡的概率為exp(-▲E/T) ， 其中E為溫度7時的內能， ▲E為其改變數。用
固體退火模擬組合優化問題，將內能E模擬為目標函式值，溫度7演化成控制引數，即得到解組合優化問題的模擬退火演算法：由初始解%和控制引數初值7開始，對當前解重複“產生新解→計算目標函式差一接受或捨棄”的迭代，並逐步減小T值，演算法終止時的當前解即為所得近似最優解， 這是基MonteCarlo迭代求解法的一種啟發式隨機搜尋過程。退火過程由冷卻進度表控制，包括控制引數的初值7及其衰
減因子K、每個7值時的迭代次數I和停止條件。

2.3 模擬退火演算法思想

模擬退火的主要思想是：在搜尋區間隨機遊走(即隨機選擇點)，再利用Metropolis抽樣準則， 使隨機遊走逐漸收斂於區域性最優解。而溫度是Metropolis演算法中的一個重要控制引數，可以認為這個引數的大小控制了隨機過程向區域性或全域性最優解移動的快慢。Metropolis是一種有效的重點抽樣法， 其演算法為：系統從一個能量狀態變化到另一個狀態時，相應的能量從E變化到E，其概率為

這樣經過一定次數的迭代，系統會逐漸趨於一個穩定的分佈狀態。
重點抽樣時，新狀態下如果向下，則接受(區域性最優)；若向上(全域性搜尋)，則以一定概率接受。模擬退火演算法從某個初始解出發，經過大量解的變換後，可以求得給定控制引數值時組合優化問題的相對最優解。然後減小控制引數7的值， 重複執行Metropolis演算法，就可以在控制引數T趨於零時，最終求得組合優化問題的整體最優解。控制引數的值必須緩慢衰減。溫度是Metropolis演算法的一個重要控制引數， 模擬退火可視為遞減控制引數7時Metro pl is演算法的迭代。開始時7值大， 可以接受較差的惡化解；隨著7的減小，只能接受較好的惡化解；最後在7趨於0時，就不再接受任何惡化解了。
在無限高溫時，系統立即均勻分佈，接受所有提出的變換。T的衰減越小， 7到達終點的時間越長； 但可使馬爾可夫(Markov) 鏈減小，以使到達準平衡分佈的時間變短。

2.4 模擬退火演算法的特點
模擬退火演算法適用範圍廣，求得全域性最優解的可靠性高，演算法簡單，便於實現；該演算法的搜尋策略有利於避免搜尋過程因陷入區域性最優解的缺陷，有利於提高求得全域性最優解的可靠性。模擬退火演算法具
有十分強的魯棒性，這是因為比起普通的優化搜尋方法，它採用許多獨特的方法和技術。主要有以下幾個方面：
(1)以一定的概率接受惡化解。
模擬退火演算法在搜尋策略上不僅引入了適當的隨機因素，而且還引入了物理系統退火過程的自然機理。這種自然機理的引入，使模擬退火演算法在迭代過程中不僅接受使目標函式值變“好”的點，而且還能夠以一定的概率接受使目標函式值變“差”的點。迭代過程中出現的狀態是隨機產生的，並且不強求後一狀態一定優於前一狀態，接受概率隨著溫度的下降而逐漸減小。很多傳統的優化演算法往往是確定性
的，從一個搜尋點到另一個搜尋點的轉移有確定的轉移方法和轉移關係，這種確定性往往可能使得搜尋點遠達不到最優點，因而限制了演算法的應用範圍。而模擬退火演算法以一種概率的方式來進行搜尋，增加了搜尋過程的靈活性。
(2)引進演算法控制引數。
引進類似於退火溫度的演算法控制引數，它將優化過程分成若干階段， 並決定各個階段下隨機狀態的取捨標準， 接受函式由Metropolis演算法給出一個簡單的數學模型。模擬退火演算法有兩個重要的步驟：一
是在每個控制引數下，由前迭代點出發，產生鄰近的隨機狀態，由控制引數確定的接受準則決定此新狀態的取捨，並由此形成一定長度的隨機Markov鏈； 二是緩慢降低控制引數， 提高接受準則， 直至控制引數趨於零，狀態鏈穩定於優化問題的最優狀態，從而提高模擬退火演算法全域性最優解的可靠性。
(3)對目標函式要求少。
傳統搜尋演算法不僅需要利用目標函式值，而且往往需要目標函式的導數值等其他一些輔助資訊才能確定搜尋方向：當這些資訊不存在時，演算法就失效了。而模擬退火演算法不需要其他的輔助資訊，而只是
定義鄰域結構，在其鄰域結構內選取相鄰解，再用目標函式進行評估。

2.5模擬退火演算法的改進方向
在確保一定要求的優化質量基礎上，提高模擬退火演算法的搜尋效率，是對模擬退火演算法改進的主要內容[9-10]。有如下可行的方案：選擇合適的初始狀態；設計合適的狀態產生函式，使其根據搜尋程序的
需要表現出狀態的全空間分散性或區域性區域性：設計高效的退火過程；改進對溫度的控制方式：採用並行搜尋結構；設計合適的演算法終止準則：等等。
此外，對模擬退火演算法的改進，也可通過增加某些環節來實現。
主要的改進方式有：
(1)增加記憶功能。為避免搜尋過程中由於執行概率接受環節而遺失當前遇到的最優解，可通過增加儲存環節，將到目前為止的最好狀態儲存下來。
(2)增加升溫或重升溫過程。在演算法程序的適當時機，將溫度適當提高，從而可啟用各狀態的接受概率，以調整搜尋程序中的當前狀態，避兔演算法在區域性極小解處停滯不前。
(3)對每一當前狀態，採用多次搜尋策略，以概率接受區域內的最優狀態，而不是標準模擬退火演算法的單次比較方式。
(4)與其他搜尋機制的演算法(如遺傳演算法、免疫演算法等)相結合。可以綜合其他演算法的優點，提高執行效率和求解質量。

3 模擬退火演算法流程
模擬退火演算法新解的產生和接受可分為如下三個步驟：
(1)由一個產生函式從當前解產生一個位於解空間的新解；為便於後續的計算和接受，減少演算法耗時，通常選擇由當前解經過簡單變換即可產生新解的方法。注意，產生新解的變換方法決定了當前新解
的鄰域結構，因而對冷卻進度表的選取有一定的影響。
(2)判斷新解是否被接受，判斷的依據是一個接受準則，最常用的接受準則是Metropolis準則：若AK 0， 則接受X作為新的當前解否則， 以概率exp(-▲E/7) 接受X”作為新的當前解X.
(3)當新解被確定接受時，用新解代替當前解，這隻需將當前解中對應於產生新解時的變換部分予以實現，同時修正目標函式值即可。此時，當前解實現了一次迭代，可在此基礎上開始下一輪試驗。若當新解被判定為捨棄，則在原當前解的基礎上繼續下一輪試驗。模擬退火演算法求得的解與初始解狀態(演算法迭代的起點)無關，具有漸近收斂性，已在理論上被證明是一種以概率1收斂於全域性最優解的優化演算法。模擬退火演算法可以分解為解空間、目標函式和初始解三部分。該演算法具體流程如下[8]：
(1)初始化：設定初始溫度7(充分大)、初始解狀態%(是算
法迭代的起點)、每個7值的迭代次數L：
(2)對k=1，…，L做第(3)至第(6)步；
(3)產生新解X；
(4) 計算增量A BE(X) -E(X) ， 其中E() ) 為評價函式：
(5)若AK0，則接受X作為新的當前解，否則以概率
exp(-▲E/7) 接受X作為新的當前解；
(6)如果滿足終止條件，則輸出當前解作為最優解，結束程式：
(7)7逐漸減小，且T→0，然後轉第(2)步。
模擬退火演算法流程如圖7.1所示。

4 關鍵引數說明
模擬退火演算法的效能質量高，它比較通用，而且容易實現。不過，為了得到最優解，該演算法通常要求較高的初溫以及足夠多次的抽樣，這使演算法的優化時間往往過長。從演算法結構知，新的狀態產生函
數、初溫、退溫函式、Markov鏈長度和演算法停止準則， 是直接影響演算法優化結果的主要環節。
狀態產生函式
設計狀態產生函式應該考慮到儘可能地保證所產生的候選解遍佈全部解空間。一般情況下狀態產生函式由兩部分組成，即產生候選解的方式和產生候選解的概率分佈。候選解的產生方式由問題的性質決
定，通常在當前狀態的鄰域結構內以一定概率產生。
初溫
溫度7在演算法中具有決定性的作用，它直接控制著退火的走向。由隨機移動的接受準則可知：初溫越大，獲得高質量解的機率就越大，且Metropolis的接收率約為1。然而， 初溫過高會使計算時間增加。
為此，可以均勻抽樣一組狀態，以各狀態目標值的方差為初溫。
退溫函式
退溫函式即溫度更新函式，用於在外迴圈中修改溫度值。目前，最常用的溫度更新函式為指數退溫函式， 即T(n+1) =KXT(n) ， 其中0<K1是一個非常接近於1的常數。
Markov鏈長度L的選取
Markov鏈長度是在等溫條件下進行迭代優化的次數， 其選取原則是在衰減引數7的衰減函式己選定的前提下，L應選得在控制引數的每一取值上都能恢復準平衡，一般L取100~1000.
演算法停止準則
演算法停止準則用於決定演算法何時結束。可以簡單地設定溫度終值T，當F=T，時演算法終止。然而，模擬火演算法的收斂性理論中要求T趨向於零，這其實是不實際的。常用的停止準則包：設定終止溫度的閾
值，設定迭代次數閾值，或者當搜尋到的最優值連續保持不變時停止搜尋。

二、案例及完整原始碼

1 案例

2 完整程式碼

%%%%%%%%%%%%%%%%%%%%%%模擬退火演算法解決函式極值%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                      %清除所有變數
close all;                      %清圖
clc;                            %清屏
D=10;                           %變數維數 
Xs=20;                          %上限                                
Xx=-20;                         %下限
%%%%%%%%%%%%%%%%%%%%%%%%%%%冷卻表引數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L = 200;                        %馬可夫鏈長度
K = 0.998;                      %衰減引數
S = 0.01;                       %步長因子
T=100;                          %初始溫度
YZ = 1e-8;                      %容差
P = 0;                          %Metropolis過程中總接受點
%%%%%%%%%%%%%%%%%%%%%%%%%%隨機選點 初值設定%%%%%%%%%%%%%%%%%%%%%%%%%
PreX = rand(D,1)*(Xs-Xx)+Xx;
PreBestX = PreX;
PreX =  rand(D,1)*(Xs-Xx)+Xx;
BestX = PreX;
%%%%%%%%%%%每迭代一次退火一次(降溫), 直到滿足迭代條件為止%%%%%%%%%%%%
deta=abs( func1( BestX)-func1(PreBestX));
while (deta > YZ) && (T>0.001)
    T=K*T; 
    %%%%%%%%%%%%%%%%%%%%%在當前溫度T下迭代次數%%%%%%%%%%%%%%%%%%%%%%
    for i=1:L  
        %%%%%%%%%%%%%%%%%在此點附近隨機選下一點%%%%%%%%%%%%%%%%%%%%%
            NextX = PreX + S* (rand(D,1) *(Xs-Xx)+Xx);
            %%%%%%%%%%%%%%%%%邊界條件處理%%%%%%%%%%%%%%%%%%%%%%%%%%
            for ii=1:D
                if NextX(ii)>Xs | NextX(ii)<Xx
                    NextX(ii)=PreX(ii) + S* (rand *(Xs-Xx)+Xx);
                end
            end            
        %%%%%%%%%%%%%%%%%%%%%%%是否全域性最優解%%%%%%%%%%%%%%%%%%%%%%
        if (func1(BestX) > func1(NextX))
            %%%%%%%%%%%%%%%%%%保留上一個最優解%%%%%%%%%%%%%%%%%%%%%
            PreBestX = BestX;
            %%%%%%%%%%%%%%%%%%%此為新的最優解%%%%%%%%%%%%%%%%%%%%%%
            BestX=NextX;
        end
        %%%%%%%%%%%%%%%%%%%%%%%% Metropolis過程%%%%%%%%%%%%%%%%%%%
        if( func1(PreX) - func1(NextX) > 0 )
            %%%%%%%%%%%%%%%%%%%%%%%接受新解%%%%%%%%%%%%%%%%%%%%%%%%
            PreX=NextX;
            P=P+1;
        else
            changer = -1*(func1(NextX)-func1(PreX))/ T ;
            p1=exp(changer);
            %%%%%%%%%%%%%%%%%%%%%%%%接受較差的解%%%%%%%%%%%%%%%%%%%%
            if p1 > rand        
                PreX=NextX;
                P=P+1;         
            end
        end
    trace(P+1)=func1( BestX);    
    end
    deta=abs( func1( BestX)-func1 (PreBestX)); 
end
disp('最小值在點:');
BestX
disp( '最小值為:');
func1(BestX)
figure
plot(trace(2:end))
xlabel('迭代次數')
ylabel('目標函式值')
title('適應度進化曲線')

%%%%%%%%%%%%%%%%%%%%%%%適應度函式%%%%%%%%%%%%%%%%%%%%%%%%
function result=func1(x)
summ=sum(x.^2);
result=summ;

三、解題過程及執行結果

1 解題過程

2 執行結果

四、matlab版本及參考文獻

1 matlab版本
2014a

2 參考文獻
[1] 包子陽,餘繼周,楊杉.智慧優化演算法及其MATLAB例項（第2版）[M].電子工業出版社，2016.
[2]張巖,吳水根.MATLAB優化演算法原始碼[M].清華大學出版社，2017.