最近报了数学建模比赛，学习一下模拟退火算法，网上搜了各种各样的学习视频和网站，索性按照自己的理解做个能看懂的总结。
1 算法思想
- SA是基于Monte-Carlo 迭代求解策略的一种随机寻优算法；出发点是基于物理中固体物质的退火过程与组合最优化问题（NP完全问题）之间的相似性。GA其实也是一种Greedy算法，但它比Greedy改进的地方就在于，它的搜索过程用到了Metropolis准则，也就是以一定的概率来接受一个比当前解要差的解，因此有可能会跳出局部最优解，找到全局最优解。
- Monte-Carlo 基本思想：利用大量采样的方法来求解一些难以直接计算得到的积分。例如，假想你有一袋豆子，把豆子均匀地朝一个形状超级不规则的图形上撒，然后数这个图形之中有多少颗豆子，这个豆子的数目就是图形的面积。当你的豆子越小，撒的越多的时候，结果就越精确。借助计算机程序可以生成大量均匀分布坐标点，然后统计出图形内的点数，通过它们占总点数的比例和坐标点生成范围的面积就可以求出图形面积。（这里我没有明白SA中具体哪里用到了Monte-Carlo思想）
- 物理中固体物质的退火过程：
首先物体刚开始处于非晶体状态（左图）。我们将固体加温至充分高，固体内部粒子随温升变为无序状，能量增大，可以自由运动和重新排列（中图）。再让其缓慢冷却，粒子渐趋有序，在每个温度都达到平衡态，最后完全冷却时能量减为最小，物体形成晶体形态，这就是退火（右图）。
2 算法步骤
look这个函数，我们要求它的最小值。
模拟退火算法是这么做的：
- 设定初始温度
，终止温度
，降温速度
，
越大，降温越慢。取一个起始点x，算出函数值 f(x)。
- 当
时，随机让 x 移动，计算新的函数值 f(x+1)；
- 根据Metropolis准则进行判断：
（1）如果新值 f(x + 1) < 当前值f(x) ，以概率1把当前值更新为新值；这很好理解，如果你发现移动到的新点求出来的值更小，肯定要更新当前值；
（2）如果f(x + 1) > f(x) ，计算概率。取一个随机数
，当
时把当前值更新为新值，否则不更新。可以发现，当前温度
- 进行一次降温：
，转到步骤2.
伪代码：
求函数f（x）的最小值
T0: 系统的初始温度（高温状态）
Tf: 温度的下限
rate: 控制降温的速率
f(now)： 系统当前状态的评价函数值
f(new): 系统新状态的评价函数值
while (T0 > Tf)
{
if f(new) - f(now) < 0
f(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解旅行商问题（MATLAB）
- 总函数：
%% I. 清空环境变量
clear all
clc
%% II. 导入城市位置数据
X = [16.4700 96.1000
16.4700 94.4400
20.0900 92.5400
22.3900 93.3700
25.2300 97.2400
22.0000 96.0500
20.4700 97.0200
17.2000 96.2900
16.3000 97.3800
14.0500 98.1200
16.5300 97.3800
21.5200 95.5900
19.4100 97.1300
20.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 = Tf
count = 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 > Tf
count = 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;
else
Obj(count) = Obj(count-1);
end
path(count,:) = S1; % 记录每次迭代的路线
T0 = q * T0; % 以q的速率降温
end
%% VII. 优化过程迭代图
figure
plot(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;
end
end
- 更新新解的函数:
function S2 = NewAnswer(S1)
% 输入：当前解S1
% 输出：新解S2
n = 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: n
for j = i + 1: n
D(i, j) = sqrt((citys(i, 1) - citys(j, 1))^2 + (citys(i, 2) - citys(j, 2))^2);
D(j, i) = D(i, j);
end
end
- 画路径函数：
function DrawPath(Route, citys)
%% 画路径函数
% 输入
% Route 待画路径
% citys 各城市坐标位置
%画出初始路线
figure
plot([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)]);
end
text(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: n
p = [p, ' —> ', num2str(Route(i))];
end
disp(p)
- 计算起点与终点之间的路径长度
function Length = PathLength(D, Route)
%% 输入
% D 两两城市之间的距离
% Route 路线
%% 输出
起点与终点之间的路径长度Length
Length = 0;
n = length(Route);
for i = 1: (n - 1)
Length = Length + D(Route(i), Route(i + 1));
end
Length = Length + D(Route(n), Route(1));
运行了上述代码后，得到的结果如下：
初始解:
11 —> 10 —> 14 —> 9 —> 5 —> 3 —> 7 —> 4 —> 8 —> 1 —> 2 —> 13 —> 12 —> 6 —> 11
总距离：62.7207
SA求得的最优解:
9 —> 11 —> 8 —> 13 —> 7 —> 12 —> 6 —> 5 —> 4 —> 3 —> 14 —> 2 —> 1 —> 10 —> 9
总距离：29.3405
随着迭代次数的增加，SA找到的最短路径变化图如下：
4 SA解旅行商问题（python，待补充）
5 对SA算法的思考
- 如果
- 理论上来说，只要
- 可以为循环结束再加一个条件，当连续m次迭代都没有更新当前最优值时，可以认为已经找到了全局最优值，可以提前结束算法。找到一个合适的m值是一个问题。
- 初始温度的选择和可行解的邻域也需要研究。
- 其他问题待补充，欢迎讨论。
文章出处登录后可见！
已经登录？立即刷新
