优化算法（三）—模拟退火算法（附MATLAB程序）

模拟退火算法（Simulated Annealing, SA）是一种基于概率的优化算法，旨在寻找全局最优解。该算法模拟金属退火过程中的物质冷却过程，逐渐降低系统的“温度”以达到全局优化的效果。它特别适用于解决复杂的组合优化问题。

一、模拟退火算法基本原理

模拟退火算法（Simulated Annealing, SA）是一种用于寻找复杂优化问题的全局最优解的随机优化算法。其基本原理可以通过以下几个核心概念和步骤来理解：

模拟退火过程：模拟退火算法受到金属退火过程的启发。金属在加热到高温后逐渐冷却，冷却过程中的原子在晶格中找到低能量的最稳定配置。类似地，模拟退火算法在解空间中进行搜索，从高温状态开始，逐渐降低温度，使解趋向于全局最优解。
接受准则：在搜索过程中，算法不仅接受当前解的邻域解，还允许接受较差的解，以避免陷入局部最优。接受较差解的概率与当前温度有关。这种接受准则可以通过以下公式描述：

其中， $\Delta E$ 是新解和当前解之间目标函数值的差（新解目标函数值 - 当前解目标函数值），T 是当前温度。

二、模拟退火算法基本公式推导

模拟退火算法的核心在于通过接受准则来决定是否接受一个新的解，即使这个新解在目标函数上可能更差。这个过程的数学基础可以追溯到统计物理中的Boltzmann分布。下面是模拟退火算法公式的详细推导过程。

2.1Boltzmann 分布

在物理学中，Boltzmann分布用于描述粒子在不同能量状态下的概率分布。对于一个能量为 E 的状态，在温度 T 下，其出现的概率为：

其中 k 是玻尔兹曼常数，Z是配分函数（用于归一化的常数），确保所有概率之和为1。

2.2目标函数与能量

在模拟退火算法中，目标函数 $f\left ( x \right )$ 可以看作是“能量”，即我们希望最小化的值。因此，假设当前解的目标函数值为 $E\left ( x \right )$ ，新解的目标函数值为 $E\left ( {x }'\right )$ ，则能量差为：

2.3接受准则

我们希望模拟退火算法能够在搜索过程中避免局部最优解，找到全局最优解。为此，我们需要决定是否接受一个新解 ${x}'$ ，即使它的目标函数值比当前解 $x$ 差。我们引入一个接受准则，使得在温度较高时，接受较差解的概率较大，而在温度降低时，接受较差解的概率减少。

根据Boltzmann分布的思想，我们可以用以下概率公式来决定是否接受一个新解 ${x}'$ ：

接受新解的概率（如果新解的目标函数值较差）：

其中 T 是当前温度， $\Delta E$ 是新解的目标函数值与当前解的目标函数值之差。这个公式是根据Boltzmann分布的形式推导出来的。

2.4温度更新

模拟退火算法中的温度 T 随时间逐渐降低，以模拟物理退火过程中的冷却。常用的温度更新公式为指数衰减：

温度更新公式：

其中 $0< \alpha < 1$ 是冷却率。这个公式表示温度在每一步迭代后乘以一个常数因子 $\alpha$ ，使得温度逐渐降低。

注意：

接受准则公式：其中，E(x)和 E(x′)分别为当前解和新解的目标函数值，T 是当前温度。
温度更新公式： $T_{next}=\alpha T$ 其中， $\alpha$ 是冷却率，通常 $0< \alpha < 1$ 。

模拟退火算法利用这些公式在解空间中进行搜索，平衡全局探索和局部优化，逐步找到全局最优解。

三、MATLAB程序仿真

下面是一个用 MATLAB 实现的模拟退火算法的示例程序。这个示例使用模拟退火算法来解决一个简单的优化问题。假设我们要最小化目标函数 $f\left ( x \right )$ ，你可以根据具体的问题调整目标函数和参数设置。

% 模拟退火算法 MATLAB 示例程序
% 目标：最小化目标函数 f(x)% 目标函数定义（这里以一个简单的二次函数为例）
objectiveFunction = @(x) x^2 - 4*x + 4; % f(x) = (x-2)^2% 参数设置
T0 = 100;       % 初始温度
Tmin = 1e-10;   % 最低温度
alpha = 0.95;   % 冷却率
maxIter = 1000; % 最大迭代次数% 初始化
x_current = randn(); % 初始解
f_current = objectiveFunction(x_current); % 初始目标函数值
T = T0; % 初始温度% 存储最优解
x_best = x_current;
f_best = f_current;% 模拟退火主循环
for iter = 1:maxIter% 生成邻域解x_new = x_current + randn(); % 在当前解附近随机生成新解f_new = objectiveFunction(x_new);% 计算目标函数值的差deltaF = f_new - f_current;% 决定是否接受新解if deltaF < 0 % 如果新解更优，直接接受x_current = x_new;f_current = f_new;else% 按概率接受较差解if rand < exp(-deltaF / T)x_current = x_new;f_current = f_new;endend% 更新最优解if f_current < f_bestx_best = x_current;f_best = f_current;end% 温度更新T = alpha * T;% 打印当前迭代信息fprintf('Iter: %d, x: %.4f, f(x): %.4f, T: %.4f\n', iter, x_current, f_current, T);% 检查是否达到停止温度if T < Tminbreak;end
end% 输出最优解
fprintf('最优解: x = %.4f, f(x) = %.4f\n', x_best, f_best);

代码说明

目标函数:
- objectiveFunction 是我们要最小化的函数。在这个示例中，我们使用了一个简单的二次函数 $f\left ( x \right )=\left ( x-2 \right )^{^{2}}$ 。你可以根据具体问题修改这个函数。
参数设置:
- T0 是初始温度。
- Tmin 是最低温度，算法在温度低于这个值时停止。
- alpha 是冷却率，控制温度的下降速度。
- maxIter 是最大迭代次数。
初始化:
- x_current 是初始解。
- f_current 是初始解的目标函数值。
- T 是当前温度。
模拟退火主循环:
- 在每次迭代中，生成一个新的解 x_new，计算其目标函数值 f_new。
- 根据目标函数值差 deltaF 和当前温度决定是否接受新解。
- 更新当前解为新解（如果接受），并更新最优解。
- 按冷却率 alpha 更新温度 T。
- 打印当前迭代的信息，包括当前解、目标函数值和温度。
输出:
- 最终输出找到的最优解和对应的目标函数值。