这是一篇关于模拟退火算法的总结博客，包括算法思想，算法步骤，Python实现，MATLAB实现，算法进阶等。

目录

1 算法思想2 算法步骤3 SA解函数最小值（python实现）4 SA解旅行商问题（MATLAB实现）5 算法进阶

1 算法思想

SA是基于Monte-Carlo 迭代求解策略的一种随机寻优算法；出发点是基于物理中固体物质的退火过程与组合最优化问题（NP完全问题）之间的相似性。GA其实也是一种Greedy算法，但它比Greedy改进的地方就在于，它的搜索过程用到了Metropolis准则，也就是以一定的概率来接受一个比当前解要差的解，因此有可能会跳出局部最优解，找到全局最优解。Monte-Carlo 基本思想：利用大量采样的方法来求解一些难以直接计算得到的积分。例如，假想你有一袋豆子，把豆子均匀地朝一个形状超级不规则的图形上撒，然后数这个图形之中有多少颗豆子，这个豆子的数目就是图形的面积。当你的豆子越小，撒的越多的时候，结果就越精确。借助计算机程序可以生成大量均匀分布坐标点，然后统计出图形内的点数，通过它们占总点数的比例和坐标点生成范围的面积就可以求出图形面积。（这里我没有明白SA中具体哪里用到了Monte-Carlo思想）物理中固体物质的退火过程：

首先物体刚开始处于非晶体状态（左图）。我们将固体加温至充分高，固体内部粒子随温升变为无序状，能量增大，可以自由运动和重新排列（中图）。再让其缓慢冷却，粒子渐趋有序，在每个温度都达到平衡态，最后完全冷却时能量减为最小，物体形成晶体形态，这就是退火（右图）。

2 算法步骤

look这个函数，我们要求它的最小值。

模拟退火算法是这么做的：

设定初始温度 T0T_0T0​，终止温度TfT_fTf​，降温速度0<rate<10<rate<10<rate<1，T0T_0T0​要取的很高，相当于加温到很高的温度。TfT_fTf​要取的很小，相当于基态，rateraterate越大，降温越慢。取一个起始点x，算出函数值 f(x)。当 T0>TfT_0>T_fT0​>Tf​时，随机让 x 移动，计算新的函数值 f(x+1)；根据Metropolis准则进行判断：

（1）如果新值 f(x + 1) < 当前值f(x) ，以概率1把当前值更新为新值；这很好理解，如果你发现移动到的新点求出来的值更小，肯定要更新当前值；

（2）如果f(x + 1) > f(x) ，计算概率P=exp(−f(x+1)−f(x)T0)P = exp(-\frac{f(x + 1)-f(x)}{T_0})P=exp(−T0​f(x+1)−f(x)​)。取一个随机数0<r<10<r<10<r<1，当r<Pr<Pr<P时把当前值更新为新值，否则不更新。可以发现，当前温度T0T_0T0​越大，指数函数括号里的值越大，接受这个新值的概率P越大。随着迭代次数的增加，T0T_0T0​会越来越小，接受新值的概率P也会越来越小，相当于渐渐冷却了下来，趋于稳定的状态。进行一次降温：T0=T0∗rateT_0=T_0*rateT0​=T0​∗rate，转到步骤2.

伪代码：

求函数f（x）的最小值T0: 系统的初始温度（高温状态）Tf:温度的下限rate: 控制降温的速率f(now)： 系统当前状态的评价函数值f(new): 系统新状态的评价函数值while (T0 > Tf){if f(new) - f(now) < 0f(now) = f(new)else if (exp(f(new) - f(now) / T0) > random(0,1))f(now) = f(new)else不更新f(now)T0 = T0 * rate;更新状态new }循环结束，得到模拟退火求得的最小值f(now)

3 SA解函数最小值（python实现）

求下列函数的最小值：

f(x)=x3−60x2−4x+6f(x)= x^3-60x^2-4x+6 f(x)=x3−60x2−4x+6

定义域为[0,100][0,100][0,100]

from __future__ import divisionimport numpy as npimport matplotlib.pyplot as pltimport math# define aim functiondef aimFunction(x):y = x ** 3 - 60 * x ** 2 - 4 * x + 6return yx = [i / 10 for i in range(1000)]y = [0 for i in range(1000)]for i in range(1000):y[i] = aimFunction(x[i])plt.plot(x, y)T = 1000 # initiate temperatureTmin = 10 # minimum value of terperaturex = np.random.uniform(low=0, high=100) # initiate xk = 50 # times of internal circulationy = 0 # initiate resultt = 0 # timewhile T >= Tmin:for i in range(k):# calculate yy = aimFunction(x)# generate a new x in the neighboorhood of x by transform functionxNew = x + np.random.uniform(low=-0.055, high=0.055) * Tif (0 <= xNew and xNew <= 100):yNew = aimFunction(xNew)if yNew - y < 0:x = xNewelse:# metropolis principlep = math.exp(-(yNew - y) / T)r = np.random.uniform(low=0, high=1)if r < p:x = xNewt += 1print(t)T = 1000 / (1 + t) #降温函数，也可使用T=0.9Tprint(x, aimFunction(x))

通过求导计算得到的精确最小值为 40.033，运行求得的最小值为40.072。

4 SA解旅行商问题（MATLAB实现）

代码来源

总函数：

%% I. 清空环境变量clear allclc%% II. 导入城市位置数据X = [16.4700 96.100016.4700 94.440020.0900 92.540022.3900 93.370025.2300 97.240022.0000 96.050020.4700 97.020017.2000 96.290016.3000 97.380014.0500 98.120016.5300 97.380021.5200 95.590019.4100 97.130020.0900 92.5500];%% III. 计算距离矩阵D = Distance(X); %计算距离矩阵n = size(D,1); %城市的个数%% IV. 初始化参数T0 = 1e10;% 初始温度Tf = 1e-30;% 终止温度q = 0.9; % 降温速率Time = ceil(double(solve([num2str(T0) '*(0.9)^x = ', num2str(Tf)]))); % 计算迭代的次数 T0 * (0.9)^x = Tfcount = 0;% 初始化迭代计数Obj = zeros(Time, 1);% 目标值矩阵初始化path = zeros(Time, n); % 每代的最优路线矩阵初始化%% V. 随机产生一个初始路线S1 = randperm(n);DrawPath(S1, X) % 画出初始路线disp('初始种群中的一个随机值:')OutputPath(S1); % 用箭头的形式表述初始路线rlength = PathLength(D, S1); % 计算初始路线的总长度disp(['总距离：', num2str(rlength)]);%% VI. 迭代优化while T0 > Tfcount = count + 1; % 更新迭代次数%temp = zeros(L,n+1);% 1. 产生新解S2 = NewAnswer(S1);% 2. Metropolis法则判断是否接受新解[S1, R] = Metropolis(S1, S2, D, T0);% Metropolis 抽样算法% 3. 记录每次迭代过程的最优路线if count == 1 || R < Obj(count-1) % 如果当前温度下的距离小于上个温度的距离，记录当前距离Obj(count) = R; elseObj(count) = Obj(count-1); endpath(count,:) = S1; % 记录每次迭代的路线T0 = q * T0;% 以q的速率降温end%% VII. 优化过程迭代图figureplot(1: count, Obj)xlabel('迭代次数')ylabel('距离')title('优化过程')grid on%% VIII. 绘制最优路径图DrawPath(path(end, :), X) % path矩阵的最后一行一定是最优路线%% IX. 输出最优解的路线和总距离disp('最优解:')S = path(end, :);p = OutputPath(S);disp(['总距离：', num2str(PathLength(D, S))]);

Metropolis准则判定函数:

function [S,R] = Metropolis(S1, S2, D, T)%% 输入% S1： 当前解% S2: 新解% D: 距离矩阵（两两城市的之间的距离）% T: 当前温度%% 输出% S： 更新后的当前解% R： 更新后的路线距离R1 = PathLength(D, S1); % 计算S1路线长度n = length(S1); % 城市的个数R2 = PathLength(D,S2);% 计算S2路线长度dC = R2 - R1; % 计算能力之差if dC < 0 % 如果能力降低 接受新路线S = S2; % 接受新解R = R2;elseif exp(-dC/T) >= rand % 以exp(-dC/T)概率接受新路线S = S2;R = R2;else % 不接受新路线S = S1;R = R1;endend

更新新解的函数:

function S2 = NewAnswer(S1)% 输入：当前解S1% 输出：新解S2n = length(S1);S2 = S1;a = round(rand(1, 2) * (n - 1) + 1); %产生两个随机位置，用于交换W = S2(a(1));S2(a(1)) = S1(a(2));S2(a(2)) = W;

计算两两城市之间的距离:

function D = Distance(citys)% 输入:各城市的位置坐标(citys)% 输出:两两城市之间的距离(D)n = size(citys, 1);D = zeros(n, n);for i = 1: nfor j = i + 1: nD(i, j) = sqrt((citys(i, 1) - citys(j, 1))^2 + (citys(i, 2) - citys(j, 2))^2);D(j, i) = D(i, j);endend

画路径函数：

function DrawPath(Route, citys)%% 画路径函数% 输入% Route 待画路径 % citys 各城市坐标位置%画出初始路线figureplot([citys(Route, 1); citys(Route(1), 1)], [citys(Route, 2); citys(Route(1), 2)], 'o-');grid on%给每个地点标上序号for i = 1: size(citys, 1)text(citys(i, 1), citys(i, 2), [' ' num2str(i)]);endtext(citys(Route(1), 1), citys(Route(1), 2), ' 起点');text(citys(Route(end), 1), citys(Route(end), 2), ' 终点');

用箭头的形式描绘出路径方向：

function p = OutputPath(Route)% 输入: 路径Route% 输出: 路径函数%给R末端添加起点即R(1)Route = [Route, Route(1)];n = length(Route);p = num2str(Route(1));for i = 2: np = [p, ' —> ', num2str(Route(i))];enddisp(p)

计算起点与终点之间的路径长度

function Length = PathLength(D, Route)%% 输入% D 两两城市之间的距离% Route 路线%% 输出起点与终点之间的路径长度LengthLength = 0;n = length(Route);for i = 1: (n - 1)Length = Length + D(Route(i), Route(i + 1));endLength = Length + D(Route(n), Route(1));

运行了上述代码后，得到的结果如下：

初始解:11 —> 10 —> 14 —> 9 —> 5 —> 3 —> 7 —> 4 —> 8 —> 1 —> 2 —> 13 —> 12 —> 6 —> 11总距离：62.7207SA求得的最优解:9 —> 11 —> 8 —> 13 —> 7 —> 12 —> 6 —> 5 —> 4 —> 3 —> 14 —> 2 —> 1 —> 10 —> 9总距离：29.3405

随着迭代次数的增加，SA找到的最短路径变化图如下：

可以看出，随着不断的迭代，距离是一直在不断下降的。我们不能保证最终的答案一定是最短的，迭代的次数越多，得到的结果就越接近正确答案。但是迭代次数过多就会需要很长的计算时间，所以要适当选择迭代次数。

5 算法进阶

1 算法特点：

如果rateraterate太小，降温就会太快，很可能最终得不到全局最优解。理论上来说，只要rateraterate足够大，降温过程足够慢，就一定可以找到全局最优解。但对于计算机来说，rateraterate太大会影响计算速度。所以要找到一个合适的方法来权衡解的性能和计算速度。可以为循环结束再加一个条件，当连续m次迭代都没有更新当前最优值时，可以认为已经找到了全局最优值，可以提前结束算法。找到一个合适的m值是一个问题。初始温度的选择和可行解的邻域也需要研究。模拟退火在求解规模较大的问题时，收敛速度很慢。

2 改进的可行方案：

设计合适的状态产生函数；设计高效的退火历程；避免状态的迂回搜索；采用并行搜索结构；避免陷入局部极小，改进对温度的控制方式；选择合适的初始状态；设计合适的算法终止准则。

3 改进的方式：

增加升温或重升温过程，避免陷入局部极小；增加记忆功能（记忆“Best so far”状态）；增加补充搜索过程（以最优结果为初始解）；对每一当前状态，采用多次搜索策略，以概率接受区域内的最优状态；结合其它搜索机制的算法；上述各方法的综合。

4 改进的思路：

记录“Best so far”状态，并即时更新；设置双阈值，使得在尽量保持最优性的前提下减少计算量，即在各温度下当前状态连续 m1 步保持不变则认为Metropolis抽样稳定，若连续 m2 次退温过程中所得最优解不变则认为算法收敛。

5 改进的退火过程

给定初温t0，随机产生初始状态s，令初始最优解s*=s，当前状态为s(0)=s，i=p=0；令t=ti，以t，s和s(i)调用改进的抽样过程，返回其所得最优解s’和当前状态s’(k)，令当前状态s(i)=s’(k)；判断C(s*)<C(s*’)? 若是，则令p=p+1；否则，令s*=s*’，p=0；退温ti+1=update(ti)，令i=i+1；判断p>m2? 若是，则转第(6)步；否则，返回第(2)步；以最优解s*作为最终解输出，停止算法。

6 改进的抽样过程

令k=0时的初始当前状态为s’(0)=s(i)，q=0；由状态s通过状态产生函数产生新状态s’，计算增量∆C’=C(s’)-C(s)；若∆C’<0，则接受s’作为当前解，并判断C(s*’)>C(s’)? 若是，则令s*’=s’，q=0；否则，令q=q+1。若∆C’>0，则以概率exp(-∆C’/t)接受s’作为下一当前状态；令k=k+1，判断q>m1? 若是，则转第(5)步；否则，返回第(2)步；将当前最优解s*’和当前状态s’(k)返回改进退火过程。

时间不变的噪声算法（TINA Time-invariant noise algorithm）

状态产生函数中扰动强度不随时间改变，而是和能量大小相关，能量大的扰动大，能量小的扰动小，能量为零，扰动也为零，算法停止。

单调升温(Monotonic temperature rising) SA

在算法退火后期，温度很低且陷入局部极小解的时，算法很难跳出。因此，可以适当重新提高温度，促使算法跳出。

记忆指导SA(Simulated Annealing with Memmory Guidance ,简记为SAMG)

增加一个记忆装置，存储算法计算过程产生的最好的解，以这个解为最终解。

自适应SA算法

根据邻域搜索进展的反馈信息, 自适应确定温度变化和邻域搜索强度特点：

退火过程中温度参数变化符合幅值递减的下降总趋势, 但不排除局部升温的可能, 以保证寻求到合适的温度序列, 避免陷入局部最优;

算法的终止条件依据退火温度和邻域搜索进展状态设计;

每一温度下算法的迭代次数随温度下降而递增, 邻域搜索强度依其对目标函数的贡献动态分配;

温度变化、邻域搜索和终止条件的控制机制由算法过程自动触发。

并行性

操作并行性：各个环节同时处理；

进程并行性：同时多个算法运行；

空间并行性：解空间分解分别处理，最终组合。

传统优化算法与全局优化算法的对比：

全局优化算法：

不依赖于初始条件；不与求解空间有紧密关系，对解域，无可微或连续的要求。求解稳健，但收敛速度慢。能获得全局最优。适合于求解空间不知的情况

传统优化算法：

依赖于初始条件。与求解空间有紧密关系，促使较快地收敛到局部解，但同时对解域有约束，如可微或连续。利用这些约束，收敛快。有些方法，如Davison-Fletcher-Powell直接依赖于至少一阶导数；共轭梯度法隐含地依赖于梯度。