數學建模常用算法-模擬退火算法

一、模擬退火算法

模擬退火的靈感來源于物理中的 “退火過程”—— 將金屬加熱到高溫后，緩慢冷卻，金屬原子會在熱能作用下自由運動，逐漸形成能量最低的穩定結構。算法將這一過程抽象為數學模型：

“溫度 T”：對應物理中的溫度，控制算法探索新解的 “自由度”—— 溫度越高，越容易接受較差的解，幫助跳出局部最優；溫度越低，越傾向于接受更優解，逐漸收斂到穩定解。

“能量 E”：對應優化問題中的目標函數值（如函數值、路徑長度等），算法的目標是找到 “能量最低” 的狀態（即目標函數的最優解）。

“冷卻策略”：溫度隨迭代逐步降低的規則，決定了算法的探索效率和收斂速度。

二、核心原理：Metropolis 接受準則

模擬退火算法的核心是Metropolis 接受準則—— 它決定了算法是否接受一個新生成的解，即使這個解比當前解更差。具體規則如下：

設當前解為\(x\)，對應的目標函數值為\(f(x)\)；新生成的解為\(x_1\)，對應的函數值為\(f(x_1)\)。

若新解更優（如求最小值時\(f(x_1) < f(x)\)），則直接接受新解，更新當前解為\(x_1\)。

若新解較差（如求最小值時\(f(x_1) > f(x)\)），則以概率\(P\)接受新解：

\(P = \exp\left( \frac{f(x) - f(x_1)}{T} \right)\)

- 溫度\(T\)越高，\(P\)越大，越容易接受較差解；

- 溫度\(T\)越低，\(P\)越小，越難接受較差解；

- 當\(T \to 0\)時，\(P \to 0\)，算法僅接受更優解，收斂到穩定狀態。

三、代碼逐行解析：從框架到細節

結合你提供的模擬退火算法框架，我們逐部分拆解代碼邏輯，理解每一步的作用和設計思路。

1. 算法參數初始化

double T = 2000; // 初始溫度：控制初始探索范圍double dT = 0.99; // 溫度衰減系數：控制冷卻速度double eps = 1e-14; // 溫度下限：算法終止條件

初始溫度\(T\)：通常設為較大值（如 2000、1000），保證初始階段有足夠的 “自由度” 探索解空間，避免過早陷入局部最優。

衰減系數\(dT\)：一般取 0.95~0.99 之間的小數，系數越接近 1，溫度下降越慢，探索越充分，但迭代次數越多；系數越小，冷卻越快，可能導致探索不充分。

溫度下限\(eps\)：當溫度低于該值時，算法基本不再接受較差解，此時可終止迭代，輸出當前最優解。

2. 目標函數定義

// 用自變量計算函數值，支持單變量或多變量（如f(x,y,z)）double func(int x, ... ) {// 此處根據具體問題實現目標函數計算double ans = ....... // 例：若求f(x)=x2的最小值，ans = x*x;return ans;}

這是算法的 “核心輸入”，需根據具體建模問題定義。例如：

- 若求解函數\(f(x) = x^2 - 4x + 5\)的最小值，func需返回\(x^2 - 4x + 5\)；

- 若解決 TSP 問題，func需計算當前路徑的總長度（輸入為城市序列，輸出為路徑和）。

支持多變量優化（如\(f(x,y)\)），只需在參數列表中添加變量（如func(int x, int y, ...)），后續新解生成和更新時同步處理即可。

3. 初始解生成

// 生成初始解（隨機值）double x = rand(); // 初始自變量x0（單變量情況）double f = func(x,...); // 初始解對應的目標函數值f(x0)

初始解通常隨機生成（利用rand()函數），保證解的隨機性，避免初始位置對結果的影響。

若問題有明確的解空間限制（如\(x \in [0, 100]\)），需在初始生成時添加約束（如x = rand() % 101），或后續通過定義域檢查修正。

4. 核心迭代：退火過程

while(T > eps) { // 溫度未降到下限，持續迭代// 步驟1：生成新解x1（幅度與溫度正相關）double dx = (2*rand() - RAND_MAX) * T;// 解釋：(2*rand() - RAND_MAX)生成[-RAND_MAX, RAND_MAX]的隨機數，// 除以RAND_MAX后范圍變為[-1, 1]，再乘以T，使新解的探索幅度隨溫度降低而減小。// 步驟2：定義域檢查（可選，根據問題需求添加）while(x + dx > 上限 || x + dx < 下限) { // 若新解超出可行域dx = (2*rand() - RAND_MAX) * T; // 重新生成新解，直到符合要求}double x1 = x + dx; // 最終合法的新解// 步驟3：計算新解的目標函數值double f1 = func(x1, ...); // f1 = f(x1)// 步驟4：Metropolis接受準則——判斷是否接受新解// 情況1：新解更優（此處假設求最小值，即f1 < f）if(f1 < f) {x = x1; // 更新自變量為新解f = f1; // 更新目標函數值為新值// 若為多變量（如x,y），需同步更新：y = y1;}// 情況2：新解較差，按概率接受else {// 計算接受概率P = exp((f - f1)/T)（求最小值時，f - f1為負數，P∈(0,1)）double prob = exp((f - f1) / T);// 生成[0,1)的隨機數，若小于P則接受新解if(prob * RAND_MAX > rand()) {x = x1;f = f1;}}// 步驟5：冷卻溫度——按衰減系數降低溫度T = T * dT;}

新解生成邏輯：新解的探索幅度與溫度\(T\)正相關，這是模擬退火的關鍵設計 —— 高溫時大步探索，低溫時精細調整，既保證全局搜索能力，又兼顧局部收斂效率。

定義域檢查：若問題對解有明確范圍限制（如資源分配問題中 “人數不能為負”），必須添加此步驟，否則可能生成無效解，導致算法出錯。

接受準則細節：需根據 “求最大值” 或 “求最小值” 調整公式：

- 求最小值：較差解的接受概率為\(P = \exp((f - f1)/T)\)（\(f1 > f\)，分子為負，\(P < 1\)）；

- 求最大值：較差解的接受概率為\(P = \exp((f1 - f)/T)\)（\(f1 < f\)，分子為負，\(P < 1\)）。

5. 輸出最優解

cout << "最優自變量x：" << x << endl;

cout << "最優目標函數值f(x)：" << f << endl;

迭代結束后，輸出當前找到的最優解（自變量\(x\)）和對應的目標函數值\(f(x)\)。

由于算法存在隨機性（初始解、新解生成均為隨機），建議多次運行算法，取多次結果中的最優值，提高結果的可靠性。

四、參數調優：讓算法更高效

模擬退火算法的性能很大程度上依賴參數設置，以下是常見的調優技巧：

初始溫度\(T\)：

- 若\(T\)過小：初始探索不充分，易陷入局部最優；

- 若\(T\)過大：迭代次數過多，算法效率低；

- 調優建議：可先設較大值（如 1000、2000），觀察迭代次數，再根據時間限制調整。

衰減系數\(dT\)：

- 若\(dT\)過小（如 0.8）：溫度下降快，探索不充分；

- 若\(dT\)過大（如 0.995）：迭代次數多，耗時久；

- 調優建議：多數問題取 0.98~0.99，平衡探索效率和收斂速度。

溫度下限\(eps\)：

一般設為\(1e-10 \sim 1e-15\)，無需頻繁調整 —— 當溫度低于此值時，算法已基本收斂，繼續迭代意義不大。

迭代次數：

若問題復雜（如多變量、解空間大），可增加 “固定迭代次數” 作為額外終止條件（如while(T > eps && iter < 1e6)），避免算法超時。

五、應用場景：模擬退火能解決哪些問題？

模擬退火算法的適用范圍極廣，尤其適合數學建模中 “復雜、多局部最優” 的優化問題，常見場景包括：

函數極值求解：單變量或多變量函數的全局最大值 / 最小值（如\(f(x,y) = x^2 + y^2 - 2x\)的最小值）。
組合優化問題：

旅行商問題（TSP）：尋找訪問所有城市且路徑最短的路線；
背包問題：在重量限制下，選擇物品使總價值最大；
圖著色問題：用最少顏色給圖的頂點著色，相鄰頂點顏色不同。

工程優化問題：資源分配、生產調度、參數優化（如機器學習模型的超參數調優）。

六、實戰案例：用模擬退火求解函數極值

以 “求解\(f(x) = x^2 - 4x + 5\)的最小值” 為例，我們完善代碼并運行：

#include <iostream>#include <cmath>#include <cstdlib>#include <ctime>using namespace std;// 目標函數：f(x) = x2 -4x +5double func(double x) {return x*x - 4*x + 5;}int main() {srand((unsigned int)time(NULL)); // 初始化隨機種子，保證每次運行結果不同// 1. 參數初始化double T = 2000.0; // 初始溫度double dT = 0.99; // 衰減系數double eps = 1e-14; // 溫度下限double x = rand() % 100; // 初始解x0（隨機取0~99的整數）double f = func(x); // 初始函數值// 2. 退火迭代while(T > eps) {// 生成新解x1double dx = (2.0 * rand() - RAND_MAX) / RAND_MAX * T; // 范圍[-T, T]double x1 = x + dx;// 計算新解的函數值double f1 = func(x1);// Metropolis準則if(f1 < f) { // 新解更優，接受x = x1;f = f1;} else { // 按概率接受較差解double prob = exp((f - f1)/T);if(prob * RAND_MAX > rand()) {x = x1;f = f1;}}// 冷卻溫度T *= dT;}// 3. 輸出結果cout << "函數f(x) = x2 -4x +5的最優解：" << endl;cout << "最優x值：" << x << endl;cout << "最小函數值：" << f << endl;// 理論最優解：x=2，f(x)=1，實際運行結果會接近此值return 0;}