{%hackmd @Hipp0/Hippotumuxthem %}

# 模擬退火算法

## 前言

最佳化問題實際上是一門非常古老的學科，它存在於各行各業中，是一門跨領域的學科。其中一
個最有名的例子是所謂的旅行商問題（Traveling Salesman Problem）。在這個例子中，該名旅行商需要跑遍 N（N > 1）個城市去推銷他的商品，而該些城市之間的距離都不一樣，這名推銷員需要從其中一個城市出發，而他老闆規定他必須把所有城市跑過一遍，請問這名旅行商應該如何繞才是最省時間（假定他的速度一直不變），也就是說，這名旅行商要找出一個最短距離的路徑。這個問題看起來簡單，實際上是一個非常複雜的數學問題。

譬如說，現在只有兩個城市，如果這名旅行商從其中一個城市出發，那這裡只有一個可供這名推銷員可選擇的路徑。假如現在有三個城市，路徑還是只有一個，因為他從其中一個城市出發繞一圈跟他反方向繞一圈的路徑距離是一樣的。意思是說假如該三個城市分別是 A，B，C。他從 A 出發，先到達 B 然後 C，再回到 A，跟他先到達 B，然後 C 再回到 A 所經過的距離是一樣的。假如現在有四個城市 A，B，C，D，那他就有三個選擇，他必須先把三個不同的路徑的距離算出來後再決定要選擇那一個。因此城市的數量愈多，可能的路徑也愈多，而且增加的速度是非線性的，十個城市的所有可能路徑就會有 1814400，假如他每一個路徑的距離都要先算出來後再作選擇，那他所需要的時間簡直就是天文數字，在實際的情況，他根本沒有那麼多的時間去做這件事，而設計一個有效的方法去找尋這個最短的路徑就是所謂最佳化問題的基本精神。
> (引用自: 物理期刊: 簡介導引模擬退火法及其應用)

另外，此問題為組合最佳化中的一個 NP-Hard 問題

### 複雜度

用動態規劃解 : $O(n^2 2^n)$
啟發式算法: 通常為多項式函數 ex: $O(n^3)$

但是啟發式算法就如同貪心一樣，不能確保為正確解。

## 啟發式算法

> 以下 "優化" 和 "最佳化" 是同一件事情。

### 啟發式算法

電腦科學的兩大基礎目標，就是發現可證明其執行效率良好且可得最佳解或次佳解的演算法。啟發式演算法則試圖一次提供一個或全部目標。

在數學優化和計算機科學中，「啟發式」（heuristic，源於希臘語 εὑρίσκω，意為「我找到、我發現」）是一種技術，當經典方法過於緩慢無法找到精確或近似解，或者經典方法在搜尋空間中無法找到任何精確解時，更快地解決問題。這種技術通過犧牲最優性(不確保最優解)、完備性、準確性或精確性來換取速度。某種程度上，啟發式方法可以視為一種捷徑。

不過有時候人們會發現在某些特殊情況下，啟發式演算法會得到很壞的答案或效率極差，然而造成那些特殊情況的資料結構，也許永遠不會在現實世界出現。因此現實世界中啟發式演算法很常用來解決問題 (因為通常都是普遍情況)。啟發式演算法處理許多實際問題時通常可以在合理時間內得到不錯的答案。

有一類的通用啟發式策略稱為元啟發演算法（metaheuristic），通常使用亂數搜尋技巧。他們可以應用在非常廣泛的問題上，但不能保證效率。

例子:(經典)
- $A^*$ 演算法
- 爬山演算法

### 元啟發式算法

> 簡單說就是想要一魚通吃，想找到所有問題都能使用某種方法解決。

(又稱萬能啟發式演算法、萬用啟發式演算法)

相較於傳統最佳化算法和迭代方法，元啟發式算法無法保證在某些類別的問題上找到全域最佳解。許多元啟發式算法實現了某種形式的隨機優化，因此找到的解依賴於生成的一組隨機變量。在組合優化中，有許多問題屬於NP完全問題，從相對低的複雜度開始就無法在可接受的時間內準確解決。此時，元啟發式算法通常能夠在比近似方法、迭代方法或簡單啟發式方法更少的計算量下提供良好的解。這一點在連續或混合整數優化領域同樣適用。因此，元啟發式算法是解決最佳化問題的有效方法。許多書籍和綜述論文都討論了此主題。根據元啟發式優化文獻回顧，Fred Glover 是提出「元啟發式」這一詞的創始人。

例子:
- 模擬退火法 (Simulated annealing algorithm, SA)
- 社會認知算法 (Social cognitive optimization, SCO)

### 無免費午餐（No-Free-Lunch）

無免費午餐（No-Free-Lunch）定理由 Wolpert 和 Macready 在1995年提出，證明了對於所有可能的優化問題，任何元啟發式算法在整體表現上都不可能優於其他算法。換句話說，沒有任何算法在所有優化問題上都能表現出色，這一定程度上限制了對元啟發式算法進行全面數學證明的可能性。

### 仿生元啟發式演算法

該類型演算法以生物的習性或群體生物行為作為靈感加以發展成為演算法。

例子
- 基因演算法 (Genetic algorithm, GA)
- 粒子群演算法 (Particle swarm optimization, PSO)
- 蟻群演算法 (Ant colony optimization, ACO)

> 還有很多奇奇怪怪的，但只要不是比較有名的通常都只是提出一個理論，然後沒有很精確的論證。

## 歷史

### 1952年：Robbins 和 Monro

貢獻：隨機優化方法
說明：提出了隨機優化的方法，為隨機搜尋在最佳化中的應用奠定基礎。

### 1954年：Barricelli

貢獻：進行了首個演化過程的模擬並應用於一般優化問題。
說明：這些模擬為後來的演化算法開辟了道路。

### 1963年：Rastrigin

貢獻：提出隨機搜尋 (random search)
說明：基於隨機性搜尋解空間，是一種早期的隨機最佳化方法。

### 1965年：Matyas

貢獻：隨機優化 (random optimization)
說明：提出了一種隨機優化策略，進一步探索了隨機性在優化問題中的應用。

### 1965年：Ingo Rechenberg

貢獻：首個演化策略算法
說明：引入了演化策略，為演化計算奠定基礎。

### 1975年：Holland

貢獻：基因演算法
說明：基因演算法的先驅，啟發了進化演算法的研究。

我寫過這篇筆記，不過只實作最簡單的 one-max 而已 -- [筆記連結](https://hackmd.io/@HIPP0/Genetic_Algorithm)

### 1983年：Kirkpatrick 等人

貢獻：模擬退火
說明：提出了模擬退火算法，用於避免局部最佳。

### 1992年：Dorigo

貢獻：蟻群最佳化 (螞蟻演算法)
說明：在博士論文中提出了蟻群最佳化，模仿蟻群尋路行為。

### 1995年：Wolpert 和 Macready

貢獻：無免費午餐定理 (no free lunch )
說明：證明了對於任意問題，沒有一種元啟發式算法能在所有情況下都優於其他算法。

## 爬山算法

爬山算法是一種局部優化的方法，是一種深度優先搜尋的改進，直白地講，就是當目前無法直接到達最優解，但是可以判斷兩個解哪個比較優的時候，根據一些回饋訊息產生一個新的可能解。

因此，爬山演算法每次在目前找到的最優方案 x 附近尋找一個新方案。如果這個新的解 $x'$ 比較優，那麼轉移到 $x'$，否則不變。

顯而易見，如果是單峰函數那肯定會收斂到全局最佳解。

> 那為啥都知道是單峰了還不三分搜尋呢?

答案是: 有些時候你不知道正解的寫法，或者本身的維度太大，很不容易寫三分搜尋。

但對於多數需要解的函數，爬山演算法很容易進入一個局部最優解，如下圖。


爬山演算法一般會引入溫度參數（類似模擬退火）。類比地說，爬山演算法就像是一隻兔子喝醉了在山上跳，它每次都會朝著它所認為的更高的地方（這往往只是個不準確的趨勢）跳，顯然它有可能一次跳到山頂，也可能跳過頭翻到對面去。不過沒關係，兔子翻過去後還會跳回來。顯然這個過程很沒有用，兔子永遠都找不到出路，所以在這個過程中兔子冷靜下來並在每次跳的時候更加謹慎，少跳一點，以到達合適的最優點。

兔子逐漸變得清醒的過程就是降溫過程，也就是溫度參數在爬山的時候會不斷減少。

關於降溫：降溫參數是略小於 1 的常數，一般在 [0.985, 0.999] 選取。
> 統計學

## 模擬退火 (SA)

> 不是 suffix array 哦 XD

模擬退火（英語：Simulated annealing，縮寫作SA）是一種逼近給定函數全局最佳的通用機率演算法，為一種元啟發算法。

「模擬退火」來自冶金學術語退火，是將材料加熱後再經特定速率冷卻的技術，目的是增大晶粒的體積，並且減少晶格中的缺陷，以改變材料的物理性質。材料中的原子原來會停留在使內能有局部最小值的位置，加熱使能量變大，原子會離開原來位置，而隨機在其他位置中移動。退火冷卻時速度較慢，使得原子有較多可能可以找到內能比原先更低的位置。

模擬退火的原理也和金屬退火的原理近似：我們將熱力學的理論套用到統計學上，將搜尋空間內每一點想像成空氣內的分子；分子的能量，就是它本身的動能；而搜尋空間內的每一點，也像空氣分子一樣帶有「能量」，以表示該點對命題的合適程度。演算法先以搜尋空間內一個任意點作起始：每一步先選擇一個「鄰居」，然後再計算從現有位置到達「鄰居」的機率。

可以證明，模擬退火算法所得解依機率收斂到全局最佳解。

### 能量函數 (目標函數)

在模擬退火算法中，優化問題的目標是找到目標函數 $f(x)$ 的極大極小值。這個目標函數在模擬退火中常被稱為能量函數 $E(x)$，目標找到 $E(x)$ 最小的解 $x$

****形式****

$$
\min\space f(x) = \min \space E(x), \space \space x \in S
$$

其中 $S$ 是解空間，包含所有可能的解。

### 能量變化

模擬退火之所以能有效搜索全域最優解，是因為它在一定概率下允許接受「能量更高」（即次優解）的解，從而跳出局部最優。這是模擬退火的重要機制，通過「能量變化」和「接受概率」來實現。

假設目前迭代解為 $x$ 下一步為 $x'$ 則能量變化為

$$
\Delta E = E(x') - E(x)
$$

如果

- $\Delta E \leq 0$ : 能量變小，接受新解 $x'$
- $\Delta E > 0$ : 新解比當前能量高，也就是新解比較差，但我們可以通過接受機率來選擇是否接受。

****接受機率****

$$
P(\Delta E) = e^{-\frac{\Delta E}{T}}
$$

這個公式來自波茲曼分佈，表示在物理系統中的粒子處於高能量狀態的概率。當溫度 T 高時，接受次優解的概率較高；當溫度下降時，系統趨向穩定，更傾向於接受能量較低的解。

### 降溫過程

在模擬退火中，溫度 T 的設定對結果影響很大。溫度通常隨著迭代次數逐漸降低，這一過程稱為降溫過程。降溫過程可以使用不同的降溫策略來實現，常見的降溫策略有以下幾種：

- 指數降溫
$$
T = \alpha^tT_0 \space \space \space 0 < \alpha < 1
$$
指數降溫能較快地降低溫度，但可能會導致過早收斂(不會收斂到全局最佳解)。

- 線性降溫
$$
T = T_0 - kt
$$
線性降溫使得溫度隨著時間緩慢降低，適合需要更平穩搜索的問題。 (同樣不確保收斂)

- 對數降溫
$$
T(t) = \frac{T_0}{log(1+t)}
$$
對數降溫較慢，理論上可以保證收斂到全域最優解，但實際應用中運行時間較長。

### 收斂性理論

根據 Geman 和 Geman 在 1984 證明了降溫速度

$$
T(t) = \frac{T_0}{log(1+t)}
$$

其中 $T_0$ 為常數。

則模擬退火在無限時間內（理論上）依機率收斂到全域最優解。這一證明依賴於馬可夫鏈的漸近性質，但在實際應用中，由於計算資源有限，無法達到無限降溫。

「依機率收斂」是機率論中的概念，表示隨著迭代次數增加，模擬退火找到全局最佳解的概率會逐漸增大，最終趨近於 1。這意味著，當我們增加迭代次數（即讓模擬退火算法的搜索過程運行得足夠長）時，找到全局最優解的可能性會愈來愈高。

數學上，我們可以用以下形式來描述依機率收斂：

$$
\lim_{t \rightarrow \infty} P(|x_t - x^*| \geq \epsilon) = 0
$$

其中:

- $x_t$ 是迭代到第 $t$ 次的解
- $x^*$ 是全局最佳解
- $P(x_t = x^*)$ 表示第 $t$ 次迭代找到全局最佳解 $x^*$ 的機率

模擬退火理論上的全局收斂性要求降溫過程非常緩慢。
- 算法需要進行非常多次迭代才能降到足夠低的溫度，因此計算時間非常長。
- 在實際問題中，這樣的慢速降溫策略會導致無法在合理的時間內完成，資源消耗過大，尤其對於大規模問題更是如此。

因此，雖然理論上保證了收斂，但在實際應用中，過慢的降溫速度使得模擬退火難以在有限時間內達到全局最優解，尤其在時間或計算資源受限的情況下。

反之如果降溫過程太快，算法可能會「凍結」在局部最優解附近，失去繼續探索其他解空間的機會，因此可能無法找到全局最優解。因此，依機率收斂到全局最佳解的證明需要溫度緩慢下降，以保證算法在更大的搜索空間中有足夠的機會找到全局最優解。

理論上，模擬退火在「無限時間」下有「概率 1」找到全局最優解，這意味著迭代次數越多，找到全局最優解的機率越高，但並不意味著有限次迭代一定能找到全局最優解。因此在實際應用中，我們通常只能在合理時間內得到「近似最優解」。

### 實際應用上的策略

設初始溫度 $T_0$、終止溫度 $T_k$、降溫係數 $d$

其中: $T_0$ 是一個比較大的數 ex: 100000，$d$ 接近 1 但小於 1 通常為 [0.985, 0.999]之間，而 $T_k$ 接近 0。

****降溫**** : 雖然對數降溫有理論在迭代次數夠多就會收斂，但通常那個太慢以至於達不到，所以較常選擇線性降溫 $T^{'} = d \cdot T$

****結束條件**** : $T < T_k$

維護答案上，並不是迭代完的那個解，而是過程中維護所有遇到的最佳解。 而有時候函數的峰太多了，可以分成很多段，每一段跑一次模擬退火，再取最佳解。

另外一種結束條件會用迭代次數 $N_{max}$ 好處是可以控制你想要的時間。

當然也可以兩種都看，如果再迭代次數前溫度就比 $T_k$ 小的話。

## 模擬退火解決 TSP 問題

TSP 的目標是找到一條路徑，使推銷員從一組城市中每個城市訪問一次後返回起點，並且總路徑最短。假設有 $n$ 個城市，定義每對城市之間的距離為 $d(i,j)$，則目標是找到一個排列 $\pi$，使得：

$$
\min \sum_{k=1}^n d(\pi(k), \pi(k+1))
$$

其中 $\pi(k)$ 表示第 $k$ 個訪問的城市，並且 $\pi(n+1) = \pi(1)$ 表示返回起點。

### 能量函數

在模擬退火解 TSP 的過程中，TSP 的「路徑總長度」被當作「能量函數」 $E(\pi)$

$$
E(\pi) = \sum_{k=1}^n d(\pi(k), \pi(k+1))
$$

### 狀態空間 (State Space)

在 TSP 問題中，狀態空間是所有可能解的集合。每個解代表一條遍歷所有城市的路徑，而在模擬退火中，每個「狀態」對應於一條具體的路徑。

假設有 $n$ 個城市，則所有可能的路徑數量為 $n!$ 排列組合數這就構成了一個巨大的狀態空間。例如，當 $n = 10$ 時，有 3628800 路徑，或者如果不考慮相反的話就是 1814400 條，在模擬退火中，我們無法遍歷所有可能路徑，因此選擇通過隨機搜索的方式在狀態空間中尋找最短路徑。

### 臨域解 (Neighborhood Solution)

在模擬退火中，每個狀態都有其對應的「鄰域解」，也就是從當前解通過一個小變動得到的新解。鄰域解決定了算法在每一步可以移動的範圍，影響算法的搜索範圍和搜索深度。

1. 交換法
- 隨機選擇兩個城市，交換它們在路徑中的位置。
- 例如: $A \rightarrow B \rightarrow C \rightarrow D$ 交換 $B$ $D$，得到 $A \rightarrow D \rightarrow C \rightarrow B$
2. 反轉法
- 隨機選擇路徑中的兩個城市，將這段城市之間的順序反轉。
- 例如: $A \rightarrow B \rightarrow C \rightarrow D \rightarrow E$ 反轉 $B$ $D$，得到 $A \rightarrow D \rightarrow C \rightarrow B \rightarrow E$
3. 移位法
- 隨機選擇一個城市，將它從當前位置移動到另一個位置。
例如: $A \rightarrow B \rightarrow C \rightarrow D$ 將 $B$ 移動到 $D$ 後面，得到 $A \rightarrow C \rightarrow D \rightarrow B$

### 能量變化

在模擬退火的每一步中，當生成了一個鄰域解後，算法會計算新解的「能量變化」，即路徑總長度的變化量：

$$
\Delta E = E(\pi') - E(\pi)
$$

### 溫度下降

雖然對數降溫有理論在迭代次數夠多就會收斂，但通常那個太慢以至於達不到，這邊選擇快速的指數降溫 $T^{'} = d \cdot T$

## 程式碼實作

採用策略:
降溫: 線性降溫
生成鄰域解: 交換法
結束條件: 迭代次數和最小溫度都看

### 城市位置

struct City {
    double x, y;
};

### 計算城市間距離

double distance(const City& a, const City& b) {
    return sqrt(pow(a.x - b.x, 2) + pow(a.y - b.y, 2));
}

### 計算路徑總長度

double totalDistance(const std::vector<int>& path, const std::vector<City>& cities) {
    double dist = 0.0;
    for (size_t i = 0; i < path.size() - 1; ++i) {
        dist += distance(cities[path[i]], cities[path[i + 1]]);
    }
    // 回到起點
    dist += distance(cities[path.back()], cities[path[0]]);
    return dist;
}

### 設定參數

採用迭代次數，你也可以改成判斷溫度

double initialTemp = 10000.0;
double toleranceTemp = 1e-8;
double coolingRate = 0.995;
int maxIter = 100000;

### 生成臨域解 (採用交換法)

std::vector<int> generateNeighbor(const std::vector<int>& path) {
    std::vector<int> newPath = path;
    int i = rand() % path.size();
    int j = rand() % path.size();
    std::swap(newPath[i], newPath[j]);
    return newPath;
}

### 模擬退火 (迭代部分)

for (int iter = 0; iter < maxIter; ++iter) {
    // 生成鄰域解
    std::vector<int> newPath = generateNeighbor(currentPath);
    double newEnergy = totalDistance(newPath, cities);
    double deltaEnergy = newEnergy - currentEnergy;

    // 決定是否接受新解
    if (deltaEnergy < 0 || exp(-deltaEnergy / temperature) > (rand() / (double)RAND_MAX)) {
        currentPath = newPath;
        currentEnergy = newEnergy;
        // 更新最佳解
        if (currentEnergy < bestEnergy) {
            bestPath = currentPath;
            bestEnergy = currentEnergy;
        }
    }

    // 降溫(指數降溫)
    temperature *= coolingRate;

    // 終止條件
    if (temperature < toleranceTemp) break;  // 溫度降到足夠低時停止
}

return bestPath;

### 完整程式碼

:::spoiler 完整程式碼

/*
 * 最佳路徑: 7 5 3 2 0 1 6 4 7
 * 長度: 25.8756 
*/

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <cstdlib>
#include <ctime>

// 城市位置結構
struct City {
    double x, y;
};

// 計算兩個城市之間的距離
double distance(const City& a, const City& b) {
    return sqrt(pow(a.x - b.x, 2) + pow(a.y - b.y, 2));
}

// 計算路徑的總距離
double totalDistance(const std::vector<int>& path, const std::vector<City>& cities) {
    double dist = 0.0;
    for (size_t i = 0; i < path.size() - 1; ++i) {
        dist += distance(cities[path[i]], cities[path[i + 1]]);
    }
    // 回到起點
    dist += distance(cities[path.back()], cities[path[0]]);
    return dist;
}

// 生成鄰域解：隨機交換兩個城市
std::vector<int> generateNeighbor(const std::vector<int>& path) {
    std::vector<int> newPath = path;
    int i = rand() % path.size();
    int j = rand() % path.size();
    std::swap(newPath[i], newPath[j]);
    return newPath;
}

// 模擬退火算法
std::vector<int> simulatedAnnealing(const std::vector<City>& cities, double initialTemp, double coolingRate, int maxIter, double toleranceTemp) {
    int n = cities.size();

    // 初始化路徑（按城市順序）
    std::vector<int> currentPath(n);
    for (int i = 0; i < n; ++i) currentPath[i] = i;
    
    // 隨機化初始路徑
    std::random_shuffle(currentPath.begin(), currentPath.end());

    double currentEnergy = totalDistance(currentPath, cities);
    std::vector<int> bestPath = currentPath;
    double bestEnergy = currentEnergy;
    double temperature = initialTemp;

    for (int iter = 0; iter < maxIter; ++iter) {
        // 生成鄰域解
        std::vector<int> newPath = generateNeighbor(currentPath);
        double newEnergy = totalDistance(newPath, cities);
        double deltaEnergy = newEnergy - currentEnergy;

        // 決定是否接受新解
        if (deltaEnergy < 0 || exp(-deltaEnergy / temperature) > (rand() / (double)RAND_MAX)) {
            currentPath = newPath;
            currentEnergy = newEnergy;
            // 更新最佳解
            if (currentEnergy < bestEnergy) {
                bestPath = currentPath;
                bestEnergy = currentEnergy;
            }
        }

        // 降溫
        temperature *= coolingRate;

        // 終止條件
        if (temperature < toleranceTemp) break; 
    }

    return bestPath;
}


int main() {
    srand(time(0));  // 設置隨機種子

    // 初始化城市位置
    std::vector<City> cities = {
        {0, 0}, {1, 5}, {2, 2}, {4, 3}, {7, 8}, {6, 2}, {3, 6}, {5, 5}
    };

    // 設定模擬退火參數
    double initialTemp = 1000.0;
    double toleranceTemp = 1e-8;
    double coolingRate = 0.995;
    int maxIter = 100000;

    // 運行模擬退火算法
    std::vector<int> bestPath = simulatedAnnealing(cities, initialTemp, coolingRate, maxIter, toleranceTemp);

    // 輸出最佳路徑和總距離
    std::cout << "最佳路徑: ";
    for (int city : bestPath) {
        std::cout << city << " -> ";
    }
    std::cout << bestPath[0] << std::endl;
    std::cout << "總距離: " << totalDistance(bestPath, cities) << std::endl;

    return 0;
}

:::

## Reference

- https://oi-wiki.org/misc/simulated-annealing/
- https://oi-wiki.org/misc/hill-climbing/
- https://ezcoding.bupt.edu.cn/ai/CODE_QUESTION/14745
- https://www.cnblogs.com/liuzhongkun/p/17286030.html
- 物理雙月刊（廿四卷二期）2002 年 4 月 簡介導引模擬退火法及其應用
- https://www.mit.edu/~dbertsim/papers/Optimization/Simulated%20annealing.pdf