数学建模常用算法-模拟退火算法

一、模拟退火算法

模拟退火的灵感来源于物理中的 “退火过程”—— 将金属加热到高温后，缓慢冷却，金属原子会在热能作用下自由运动，逐渐形成能量最低的稳定结构。算法将这一过程抽象为数学模型：

“温度 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)=x²的最小值，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) = x² -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) = x² -4x +5的最优解：" << endl;cout << "最优x值：" << x << endl;cout << "最小函数值：" << f << endl;// 理论最优解：x=2，f(x)=1，实际运行结果会接近此值return 0;}