夜深人靜寫演算法（十二）- 模擬退火

阿新 • • 發佈：2018-12-05

一、引例 1、函式最值 函式最值分為函式最大值和函式最小值，最小值即定義域內函式的最小值，最大值即定義域內函式的最大值。函式最大（小）值的幾何意義為函式影象的最高（低）點的縱座標。那麼，讓我們來看幾種簡單的情況： 1) 一次函式

圖一-1-1 畫出函式影象如下：

圖一-1-2 在定義域[x1,x2]內，函式的最小值和最大值分別取在兩端點上。 2) 二次函式

圖一-1-3

圖一-1-4 由於該函式存在極值點，所以定義域的不同，函式最值也不盡相同。當然也可以通過函式影象分析出來，分三種情況討論：

a)當定義域[x1,x2]滿足x2<=x0時，最大值取在x1，最小值取在x2； b)當定義域[x1,x2]滿足x1>=x0時，最大值取在x2，最小值取在x1； c)當定義域[x1,x2]滿足x1<x0<x2時，最大值為兩端點的小者，最小值取在x0； 3) 周期函式

圖一-1-5

圖一1-6 這個函式為標準的正弦函式，週期為2π，所以計算最值的時候可以將定義域平移到[0,2π)區間，然後分情況討論即可（參考二次函式的情況）。接下來，讓我們看下更加複雜的情況下的函式最值。 4) 複合函式

圖一-1-7 這種函式乍看之下，很難想象它的函式影象，那麼我們可以利用描點法繪製出它的函式影象，點越多越精確。得到了一個令人絕望的函式影象。

圖一-1-8 這種情況下，如何求它在給定的定義域下的函式最值呢？這就是本文今天要介紹的演算法---模擬退火。模擬退火演算法是一種通用的優化演算法，理論上演算法具有概率的全域性優化效能，目前已在工程中得到了廣泛應用，諸如VLSI、生產排程、控制工程、機器學習、神經網路、訊號處理等領域。 二、袋鼠跳 1、從天而降的袋鼠 在介紹模擬退火之前，我們先把上面提到的問題用一種有趣的方式描述出來，這樣有助於理解模擬退火的演算法核心。我們將上文的函式定義為y=f(x)，影象想象成一座山，上面有無數個連綿起伏的山峰和山谷，現在我們要想辦法找到這座山的最低谷，如圖二-1-1所示。

圖二-1-1 首先，我們在x座標上隨機採點，從天而降一些袋鼠。每隻袋鼠有自己的思維，它們都可以各自選擇是往左跳還是往右跳。

圖二-1-2 巨集觀來說，它們的行為是並行的（實際上，只需要用一個列舉來模擬每隻袋鼠的行為，並非真正的CPU平行計算），所以我們其實只需要考慮某一隻袋鼠的行為。 2、選擇性跳躍 對於每隻袋鼠，我們首先定義一個步長S，如果當前袋鼠的位置為x，那麼在x-S和x+S兩個位置中挑選一個相對較低的位置（即min{f(x-S),f(x+S)}），比較這個位置和當前袋鼠所在的位置。如果比當前點更低，那麼毫不猶豫的跳過去，否則以一定概率選擇性的跳躍。

圖二-1-3 3、退火（減少步長） 然後，等比例減少步長S（一般是乘上一個係數，比如0.99），繼續迭代這樣的操作，直到步長S小於某個精度，那時候袋鼠也跳不動了。這樣做，只要引數控制的好，袋鼠有很大機率落在最低點，但是也不一定。所以我們選了很多袋鼠，然後並行迭代後取最後落在最低點的那隻即可。好了，如果能夠聽懂，說明你已經無形之中領悟了模擬退火的精髓。 三、模擬退火 1、初衷 模擬退火演算法來源於固體退火原理，是基於蒙特卡羅迭代求解策略的一種隨機概率尋優演算法，其出發點是基於物理中固體物質的退火過程與一般組合優化問題之間的相似性。模擬退火演算法從某一較高初溫出發，隨著溫度引數的不斷下降，結合概率的跳躍特性，在解空間中隨機尋找目標函式的全域性最優解。即在區域性最優解時能夠概率性的跳出，並且最終趨於全域性最優。根據Metropolis準則，粒子在溫度T時趨於平衡的概率為：

（E為溫度T時的內能，ΔE為其改變數，k為Boltzmann常數）模擬退火演算法更加直觀的一種解釋：想象成一個坑坑窪窪的土地，從天上丟下來一些彈性非常好的球，由於能量損失，這些球每次的彈跳高度都會下降，直到一直彈不起來。那麼一開始它總能彈到一些本來不屬於初始下降區域的地方，即跳出區域性最優解。慢慢逼近全域性最優解。 2、演算法過程 (0) 初始化：初始溫度T(充分大)，初始解狀態X(是演算法迭代的起點)，每個溫度T下的迭代次數L； (1) for i = 1 to L 執行(2)(3)(4)(5) (2) 產生新解X′（以某種策略生成新解）； (3) 計算溫度增量ΔT=h(X′)-h(X)，其中h(X)為估價函式； (4) 若ΔT<0則接受X′作為新的當前解，否則以概率e(-ΔT/T)接受X′作為新的當前解； (5) 如果滿足終止條件則輸出當前解作為最優解，結束程式。終止條件通常取為連續若干個新解都沒有被接受時終止演算法，否則當溫度到達某個最小精度時結束。 (6) 否則，逐漸減少T，然後轉(1)。 3、演算法剖析 這裡的溫度T就是袋鼠的跳躍步長，會隨著時間的推移逐漸減少（等比下降，即退火過程）；初始解X是隨機生成的，可以認為是上述問題中函式影象的x座標，L就是袋鼠的個數。步驟(0)即初始化；步驟(1)則是平行計算每隻袋鼠的跳躍情況；步驟(2)產生左右兩個解，並且取更優解作為候選解X'。步驟(3)這一步出現了一個新的名詞，估價函式。估價函式即對當前解的一個評估，從而決定是否要從當前解轉向候選解。很顯然，在求函式最小值的問題中，估價函式正好是函式值本身，即h(X)=f(X)。求最大值時，估價函式是函式值的相反數，即h(X)=-f(X)（原因是在該問題的前提下，估價函式一定是越小越優）。步驟(4)對應了袋鼠的選擇性跳躍中的“選擇”二字。如果能夠確認更優，則直接替換當前解，否則按照上述概率進行篩選。逐漸迭代溫度T，直接溫度滿足某個精度為止，這個很好理解，因為溫度對應步長，當溫度很小時，左右移動意義也不大了。 4、演算法實現 1）類定義

class simulatedAnnealing {
    static const double minTemperature;      // 穩態（最低）溫度
    static const double deltaTemperature;    // 溫度下降率
    static const int solutionCount;          // 並行解個數
    static const int candidateCount;         // 每個解的迭代次數
private:
    Point3D bound;
    Point3D x[MAXC];
    Point3DSet pointSet;
    double temperature;

    bool valid(const Point3D& pt);
    double randIn01();
    Point3D getRandomPoint();
    Vector3D getRandomDirection();
    Point3D getNext(const Point3D& now);
public:
    void start(double temper, Point3D ptBound, Point3DSet& ptSet);
    double evaluateFunc(const Point3D& pt);
    Point3D getSolution();
    static simulatedAnnealing& Instance();
};
const double simulatedAnnealing::minTemperature = 1e-4;
const double simulatedAnnealing::deltaTemperature = 0.95;
const int simulatedAnnealing::solutionCount = 10;
const int simulatedAnnealing::candidateCount = 30;

2）退火主流程
三個引數temper代表退火的初始溫度、ptBound代表最優解的範圍、ptSet代表資料集。 temper：初始化溫度，即初溫temper，直接賦值給類變數temperature即可。 ptBound：從ptBound範圍內隨機挑選solutionCount個點作為初始點集，ptBound表示的是一個點，所以實際的初始點的範圍就是從原點到ptBound隨機取值，然而對於最優解在負數座標的情況，需要將整個座標系平移，即平移到所有解都位於原點之上。這裡的點是個抽象的概念，可以是一維的（上文提到的x-y函式），也可以是二維的（平面座標系上的點），也可以是三維的（空間座標系中的點），當然也可以是更高維度的。我寫的是三維的模板，如果要轉換成二維只要將第三維座標z置0即可。 ptSet：輔助資料，用於求特定的問題，例如問題描述為求一個任意多邊形的最大內接圓，那麼這裡任意的多邊形就用這個ptSet點集來表示。退火主流程分成三步走： a）初始化引數和初始解集生成（程式碼中的步驟0和1）；初始溫度選擇很重要，選的越大迭代時間越久越精確，會以時間作為代價。所以需要權衡兩者之間的平衡。 b）對每個初始解進行最優化篩選和替換（程式碼中的步驟2、3、4、5、6、7） minTemperature為最低溫度控制，即當temperature小於等於這個溫度時演算法終止，一般問題如果求解精度在0.01，那麼最低溫度最好控制在0.001，即儘量比求解精度小1個數量級，這是為了提高準確性。 solutionCount作為初始解集數量，為事先定義的常量，在情況比較特殊（比如單調性很明顯的單峰單谷問題）時可以取1，取得越大效率越慢，但是更加精確；相反，效率高但是精度有可能會下降。 candidateCount為當前某個解的下一候選解的數量，即上文提到的袋鼠選擇往左還是往右，對於一維的情況，即函式求極值，只需要將candidateCount設為2即可（因為只有左右兩個方向），對於更高維的情況，則需要隨機方向，所以取值可能是30、50、100，具體問題具體分析，通常是實驗結果。 getNext(x[i])返回的是當前解x[i]的下一個候選解，並且用evaluateFunc計算估價函式後，取估價值最小的作為實際候選解。這個候選解需要迭代candidateCount次，選取最優，然後和當前解x[i]進行估價函式的比較。如果候選解的估價值更小，則直接替換當前解；否則，以一定概率接收候選解。注意：這裡的一定概率也可能是零概率，也就是第二個分支永遠跑不到，比如求單谷函式的最小值，完全沒有必要接收估價更大的解，這種情況下，模擬退火演算法退化成“梯度下降”演算法。 c）退火（程式碼中的步驟8） temperature 直接乘上一個常量deltaTemperature，一般取0.9左右，根據具體情況而定，越大精度越高，時間越慢；反之，則精度越低，時間越快。
來看下演算法實現：

void simulatedAnnealing::start(double temper, Point3D ptBound, Point3DSet& ptSet) {
    // 0.初始化溫度
    temperature = temper;
    bound = ptBound;
    pointSet = ptSet;
    int i, j;

    // 1.隨機生成solutionCount個初始解
    for(i = 0; i < solutionCount; ++i) {
        x[i] = getRandomPoint();
    }

    while (temperature > minTemperature) {
        // 2.對每個當前解進行最優化選擇
        for(i = 0; i < solutionCount; ++i) {
            double nextEval = INF;
            Point3D nextOpt;
            // 3.對於每個當前解，隨機選取附近的candidateCount個點，並且將最優的那個解保留
            for(j = 0; j < candidateCount; ++j) {
                Point3D next = getNext(x[i]);
                if(!valid(next)) {
                    continue;
                }
                double Eval = evaluateFunc(next);
                if(Eval < nextEval) {
                    nextEval = Eval;
                    nextOpt = next;
                }
            }

            // 4.沒有生成可行解
            if(nextEval >= INF)
                continue;

            // 5.計算生成的最優解和原來的解進行比較
            double deltaEval = evaluateFunc(nextOpt) - evaluateFunc(x[i]);
            if(deltaEval < 0) {
                // 6.比原來的解更優，直接替換
                x[i] = nextOpt;
            }else {
                // 7.沒有原來的解優，則以一定概率進行接收
                // 這個概率上限會越來越小，直到最後趨近於0
                // 理論上，這個分支也可能不考慮
                if( randIn01() < exp(-deltaEval/temperature) ) {
                    x[i] = nextOpt;
                }    
            }
        }
        // 8.退火
        temperature *= deltaTemperature;
    }
}

       3）估計函式
模擬退火中的估價函式，類似A*，即通過某個解估算出一個價值（整數或者浮點數皆可），我們定義為估價函式越小，解越優。不同問題，估價函式的實現也不盡相同，這個在下一章節中會詳細介紹。      4）其它的一些小函式就不介紹了，完整C++程式碼，請參見：模擬退火演算法C++實現
        5、演算法複雜度
       最後，分析一下實際退火過程的時間複雜度。首先最外層迭代，從最高溫度到最低溫度經過一個等比式的迭代，令初溫為T，最低溫度TMin，溫度下降率DT，下降次數為N，那麼有：

從而得出N的最大值為：

然後來看每次迭代，解集數量定義為L，每個解的候選解數量為C，每次計算估價函式的時間複雜度為H，則得到演算法的時間複雜度為：

      得到演算法複雜度後，我們就可以根據實際情況，設定L和C的值，權衡時間和精度進行演算法求解了。

      四、模擬退火演算法舉例
      最後，我們對這個演算法展開一些問題的探索，加深對演算法本身的理解。       1、函式最值
      以函式最小值為例，回到一開始提到的那個複合函式：

【例題1】給定x1,x2，求該函式在區間[x1, x2]上的的最小值，精確到小數點後2位。直接套用模擬退火，初始溫度可以取函式定義域的區間長度x2-x1，最小溫度0.001，溫度下降率0.95，初始解集10個，候選解集2個，估價函式直接取函式值，實現如下：

double simulatedAnnealing::evaluateFunc(double x) {
    double sx = sin(x);        // sin(x)計算比較費時，預先計算出來
    return sx*sx - 2*sx + x*sx;
}

2、曲面最值
【例題2】給定一個二次曲面方程，保證這個曲面有一個最低點，求這個最低點的座標。

圖4-2-1 這個問題其實就是求一個座標，這個座標在曲面上，且座標的縱座標分量最小（即y座標）。因為有曲面方程，所以如果x和z的值確定，那麼求解一個一元二次方程就能算出y，所以起初我們可以任意隨機一個曲面上的點，然後通過模擬退火不斷將點逼近到最低點。估價函式的實現比較簡單：

double simulatedAnnealing::evaluateFunc(const Point3D& pt) {
    return pt.y;
}

3、最 遠的最近點【例題3】給定一個矩形區域W,H(W,H<=10000) 和 N(N<=1000)個危險點，求矩形內一個點，這個點離最近的那個危險點最遠（要求精確到小數點後1位）。如圖紅色點代表危險點，藍色點代表離最近的危險點最遠的那個最優解所在的位置。

圖四-3-1 初始化溫度T=max(W,H),隨機初始化一個矩形內的點作為初始解，將初始解作為當前解進行迭代計算。每次以當前解(圖中的黑點)為圓心,溫度T為半徑生成一些方向隨機的候選點集(圖中天藍和深藍的點)，然後選擇一個離所有危險點(圖中紅色點)最小距離最大的點(圖中深藍色點)作為當前解的候選替代者（暫時不進行替換）。

圖四-3-2
模擬退火演算法中，估價函式的值越小越優，而這題求的是最大距離，也就是離最近的那個危險點距離更大的點更優，所以估價函式需要在求出最小距離後加上一個相反數。將候選最優解和當前解的估價函式作差，如果估價差值△E<0，說明候選最優解比當前解更優，則用候選最優解替換當前解；否則，隨機一個0-1的隨機數p，如果p<e(-△E/T)，則用候選最優解替換當前解，否則不進行替換。每次計算完畢，則將溫度下降一定比率，直到溫度退化到一定精度，則演算法結束。這個演算法中，溫度T為每次選擇候選點集的偏移半徑；每個候選點到最近危險點的距離的相反數作為估價函式；估價函式差值作為退火時的能量差值作為新解取捨的判斷依據。當然，初始化生成的位置可能不夠好，導致最後陷入區域性最優解，所以，我們可以隨機生成一個初始解集合進行“並行”（並非真正的並行）計算。引數設定：最小溫度0.001，溫度下降率0.95，初始解集5個，候選解集30個，估價函式實現如下：

double simulatedAnnealing::evaluateFunc(const Point3D& pt) {
    double minDist = INF;
    for(int i = 0; i < pointSet.n; ++i) {
        double dist = (pointSet.p[i] - pt).len();
        if(dist < minDist) 
            minDist = dist;
    }
    return - minDist;
}

      其中pointSet就是上文我們提到的輔助資料，用來記錄題目中那些“危險點”。
       4、最近的最遠點
      【例題4】給定一個矩形區域W,H(W,H<=10000) 和 N(N<=1000)個點，求矩形內一個點，這個點離最遠的那個點最近（要求精確到小數點後1位）。這個問題和例題3，正好相反，同樣可以利用模擬退火的性質求解，唯一不同的就是估價函式不需要取相反數，並且選擇點的時候是選取距離最大的那個點。
       5、任意多邊形的最大內接圓
      【例題5】給定一個任意簡單多邊形（可能是凹多邊形）以及一個半徑為r的圓，問這個圓能否放入這個多邊形內。

圖四-5-1 作為一個判定性問題，同樣可以用模擬退火來進行判定。首先，初始化一些初始點集，對於每個點，以它到所有多邊形線段的距離的最小值的相反數作為估價函式。因為我們想要的是讓這個點到多邊形所有線段的最小距離最大，所以估價函式需要取相反數。仔細一想，這個問題和“最遠的最近點”是一樣的，唯一不同的是它有跳出條件，即當這個最小距離已經大於等於r時，可以跳出迭代。
需要注意的是，圓的圓心需要在多邊形內，否則可能出現以圖四-5-2的非法情況。