MATLAB 模拟退火算法
模拟退火算法原理
模拟退火算法是一种通用概率算法，用来在一个大的搜寻空间内寻找问题最优解。
其思想源于固体的退火过程：将固体加热至足够高的温度，缓慢冷却，其内能就由很大缓慢趋于内能最小。
足够高的温度对应随机的解，缓慢冷却对应对解进行一次随机扰动，内能对应目标函数。即将随机解不断进行扰动，根据目标函数变化以一定概率接受解，不断重复，解就会趋近于最优解。
接受新解概率——Metropolis准则：
设从当前状态 \(i\) 生成新状态 \(j\) ，若 \(E_j<E_i\)，则新状态 \(j\) 作为当前状态，否则，以概率 \(exp[\frac{-(E_j-E_i)}{\beta×T}]\) 接受当前状态。\(\beta\) 为常数，通常取1；\(t\) 是此时温度。
内能对应由于目标函数（求最小值，如果求最大值可加负号求最小），如果新解优于旧解，则接受，否则以一定概率接受。
就是因为模拟退火能够有一定概率接受较差的解，所以才可能跳出局部最优解，寻找到全局最优解。模拟退火的效果和 \(\beta\) 有非常大的关系，最后我们详细分析。
模拟退火算法步骤
参数名 | 说明 | 取值要求 |
|---|---|---|
T0 | 初始温度 | T0应足够大，但也应考虑计算量 |
T | 当前温度 | |
Tf | 终止温度 | 应足够小，比如0.01~5 |
α | 温度的衰减参数，T{k+1}=αT | 可取0.5~0.99；较大的α搜索范围更大 |
Lk | Markov链长度，即相同温度下循环次数 | 当 α 较大时，为减小循环次数，Lk可较小 |
算法步骤
令T = T0，随机生成初始解 x0，计算目标函数 E(x0)
根据实际情况对当前解进行扰动，产生新解，再次计算目标函数。
比较目标函数，根据Metropolis准则接受新解。
重复扰动 Lk 次，达到次数后再按系数α降低温度。
判断T是否达到Tf ，达到则终止算法，否则继续扰动。
注： 一般会把退火过程中碰到的最优解也保存下来，与最后一个解一同输出，因为有可能在以一定概率接受较差解的时候接受了较差解。
扰动方法
对应实际问题，都有不同的扰动方法以及不同的目标函数设计方法。
解 x 一般是能够表示一种实际方法的量，比如一个数组，根据x可以得出目标函数E的值。
根据实际情况如何变动，设计对解 x 进行改变，也就实现了扰动。
产生的随机扰动要保证可能取遍所有可能情况。
实例——旅行商问题
问题描述
一名商人要到n个不同的城市去推销商品，已知两两城市之间的距离，选择一条路径使商人经过所有城市并回到起点，使路程总长度最短。
搭建退火算法框架
我们先事先准备好实际问题所需要的数据，比如两两城市之间的距离矩阵等。
除此之外还要初始化退火算法的各种初始值，搭建温度循环和Markov链循环两个循环框架。对于解x和最优结果E，要设置三个参数new、current、best，分别用来存储当前值，上一轮的结果，历史最优解。
clear
alpha = 0.9;
beta = 0.01;
T0 = 97;
Tf = 3;
T = T0;
Markov_length = 10000;
coolingTime = ceil((log(Tf) - log(T0))/log(alpha));
coordinates_x = [565 25 345 945 845 880 25 525 580 650 1605 1220 1465 ...
1530 845 725 145 415 510 560 300 520 480 835 975 1215 1320 1250 660 ...
410 420 575 1150 700 685 685 770 795 720 760 475 95 875 700 555 830 ...
1170 830 605 595 1340 1740];
coordinates_y = [575 185 750 685 655 660 230 1000 1175 1130 620 580 200 5 ...
680 370 665 635 875 365 465 585 415 625 580 245 315 400 180 250 555 ...
665 1160 580 595 610 610 645 635 650 960 260 920 500 815 485 65 610 ...
625 360 725 245];
amount = size(coordinates_x,2);
dist_matrix = zeros(amount);
for i = 1:amount
for j = 1:amount
dist_matrix(i,j) = sqrt((coordinates_x(i)-coordinates_x(j)).^2 +...
(coordinates_y(i)-coordinates_y(j)).^2);
end
end
sol_new = 1:amount;
sol_current = sol_new;
sol_best = sol_new;
E_new = inf;
E_current = inf;
E_best = inf;
E_bestHistory = zeros(1, coolingTime);
E_currentHistory = zeros(1, coolingTime);
loopTimes = 0;
while T >= Tf
loopTimes = loopTimes + 1;
for r = 1:Markov_length
% 扰动程序
end
T = T.*alpha;
E_bestHistory(1, loopTimes) = E_best;
E_currentHistory(1, loopTimes) = E_current;
end
扰动的程序
对实际问题，核心是如何扰动。本题开始时初始化城市序号数列sol，按照默认顺序初始化可以满足初始值的E很大，然后扰动可以设计为任意交换两个城市的次序。根据sol我们可以转换成实际中的路线。
ind1 = 0; ind2 = 0;
while(ind1 == ind2)
ind1 = randi(amount);
ind2 = randi(amount);
end
tmp1 = sol_new(ind1);
sol_new(ind1) = sol_new(ind2);
sol_new(ind2) = tmp1;
E_new = 0;
for i = 1:(amount - 1)
E_new = E_new + dist_matrix(sol_new(i),sol_new(i+1));
end
E_new = E_new + dist_matrix(sol_new(amount),sol_new(1));
if E_new < E_current
E_current = E_new;
sol_current = sol_new;
if E_new < E_best
E_best = E_new;
sol_best = sol_new;
end
else
if rand <exp(- (E_new - E_current) ./T ./ beta)
E_current = E_new;
sol_current = sol_new;
else
sol_new = sol_current;
end
end
完整程序
clear
alpha = 0.9;
beta = 0.01;
T0 = 97;
Tf = 3;
T = T0;
Markov_length = 10000;
coolingTime = ceil((log(Tf) - log(T0))/log(alpha));
coordinates_x = [565 25 345 945 845 880 25 525 580 650 1605 1220 1465 ...
1530 845 725 145 415 510 560 300 520 480 835 975 1215 1320 1250 660 ...
410 420 575 1150 700 685 685 770 795 720 760 475 95 875 700 555 830 ...
1170 830 605 595 1340 1740];
coordinates_y = [575 185 750 685 655 660 230 1000 1175 1130 620 580 200 5 ...
680 370 665 635 875 365 465 585 415 625 580 245 315 400 180 250 555 ...
665 1160 580 595 610 610 645 635 650 960 260 920 500 815 485 65 610 ...
625 360 725 245];
amount = size(coordinates_x,2);
dist_matrix = zeros(amount);
for i = 1:amount
for j = 1:amount
dist_matrix(i,j) = sqrt((coordinates_x(i)-coordinates_x(j)).^2 +...
(coordinates_y(i)-coordinates_y(j)).^2);
end
end
sol_new = 1:amount;
sol_current = sol_new;
sol_best = sol_new;
E_new = inf;
E_current = inf;
E_best = inf;
E_bestHistory = zeros(1, coolingTime);
E_currentHistory = zeros(1, coolingTime);
loopTimes = 0;
while T >= Tf
loopTimes = loopTimes + 1;
for r = 1:Markov_length
ind1 = 0; ind2 = 0;
while(ind1 == ind2)
ind1 = randi(amount);
ind2 = randi(amount);
end
tmp1 = sol_new(ind1);
sol_new(ind1) = sol_new(ind2);
sol_new(ind2) = tmp1;
E_new = 0;
for i = 1:(amount - 1)
E_new = E_new + dist_matrix(sol_new(i),sol_new(i+1));
end
E_new = E_new + dist_matrix(sol_new(amount),sol_new(1));
if E_new < E_current
E_current = E_new;
sol_current = sol_new;
if E_new < E_best
E_best = E_new;
sol_best = sol_new;
end
else
if rand <exp(- (E_new - E_current) ./T ./ beta)
E_current = E_new;
sol_current = sol_new;
else
sol_new = sol_current;
end
end
end
T = T.*alpha;
E_bestHistory(1, loopTimes) = E_best;
E_currentHistory(1, loopTimes) = E_current;
end
讨论退火参数对结果的影响
在你理解了上述关于此算法的基本操作以后，我们来分析此算法要获得好的效果的参数条件。
总结来说，这个算法就是两个循环，外循环会改变温度，内循环是普通的重复。算法的精髓就在于这个温度的改变，温度所起的作用就是控制差解的接受概率。如果一直不接受差解，那么就和普通的循环完全没有区别，如果概率一直很高，那么解就会一直乱跑，同样和普通循环一样。我们希望的是，在温度比较高的时候，接受差解的概率高，而在温度低的时候，接受差解的概率低，这样才会使解稳定。
那么，我们首先要保证的是，降温系数足够大，初始温度足够大，结束温度足够低，同一温度下循环次数足够大，这是为了保证有足够多的温度能够进行算法。否则只经历两三个温度，肯定是不行的。
接下来，我们详细分析接受概率的函数，函数表达式如下：
\[ exp[\frac{-(E_j-E_i)}{\beta×T}] \]
温度 T 在分母上，这就是为了保证上面说的，在 T 高的时候概率大，在 T 低的时候概率小，T 的初值和终值的范围要足够大才行。
分子是目标函数结果相减，说明离当前接受解越接近，则概率越大，这样可以避免一个很差的解被接受。
最关键的就是 \(\beta\) ，我看有的书上写的是通常取1，这是个很可怕的事情，这个值一定是要根据实际情况来选的，直接决定这个算法是否有效。很简单的道理，每个实际问题，E 的范围是不定的，可能只有0-10，也可能是几万以上，而这个 \(\beta\) 就相当于避免这种影响，E 如果很大，\(\beta\) 就要大，E 很小 \(beta\) 就要小。至于如何直观选取呢？最简单的办法就是作出随降温次数变化的当前接受解和当前最优解，根据曲线很好判别。
以上面的例子为例，我们已经记录了每个文档下的 E_current 和 E_best ，画曲线：
figure()
plot(E_currentHistory)
hold on
plot(E_bestHistory)
% FigureHandle2D 是为了方便自定义的类
fh = FigureHandle2D();
fh.plotDecoration('降温次数', '目标函数', '随降温次数变化情况', ...
{'当前接受解', '当前最优解'})
fh.setFontSize(16)
当 \(\beta\) 取 1，结果比较理想，如下：
可以看到，一开始红线和蓝线离得比较开，说明接受值不是出现过的最优值，有利于在全局寻找，但是在温度降下来以后，就基本重合，此时就比较稳定，寻找这个范围内的局部最优。
我们在确定 \(\beta\) 以后，就可以提高 \(\alpha\) 让降温次数变多，来跑程序。
如果 \(\beta\) 设置的比较大，比如 100 ，就会变成这样：
很显然，红线和蓝线一直离得很开，说明接受概率一直很大，这样就会不停乱接受，结果也会不好。
如果 \(\beta\) 很小，比如 0.001，就会变成这样：
很明显，红线和蓝线一直重合，说明根本没有接受过比较差的解，这样就和一直循环是一样的，结果也会不好。
总结
退火算法其实很简单，框架是固定的，关键就在于如何去扰动，只要保证扰动的改变很小，而且基本等概率会取到所有可能的值就可以了。
参数的合理选择才能保证算法有效性，否则算法毫无意义，尤其掌握 \(\beta\) 的值。
本文章使用limfx的vsocde插件快速发布