【優化求解】基於matlab模擬退火算法求解函數極值問題【含Matlab源碼 1203期】

1、粒子羣算法簡介

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趨向於零，這實際上是不實際的。經常使用的中止準則包：設置終止溫度的閾
值，設置迭代次數閾值，或者當搜索到的最優值連續保持不變時中止搜索。

2、案例及完整源代碼

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;

3、解題過程及運行結果

1 解題過程

2 運行結果

4、matlab版本及參考文獻

1 matlab版本
2014a

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