模拟退火算法求解组合优化（附代码详解）
最开始接触模拟退火算法（SA）还是去年参加华为杯数学建模的时候，发现求解组合优化效果挺不错的，后续一些数模比赛涉及组合优化的也基本都是用他求解。靠着参加几个数模，8月份甚至天真的做了一份运筹优化的简历，投了顺丰，美团，阿里菜鸟的运筹优化岗位，真是too young, too simple. 菜鸟，美团直接一面挂，顺丰幸运的到了HR面，然而也没了消息，哈哈哈，这个岗位挺冷门的，虽然投的人比较少，但是可能要求也比较高吧，后面就老老实实投开发了。
SA算法模拟的是材料热处理的退火过程（虽然是机械的，但是真没学会），SA算法在某一初温下，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。SA是一种相对而言比较好实现的算法，但是因为需要较高的初始温度以及较低的降火速率和终止温度，因此目标函数调用次数较多，不适用于目标函数评估时间过长的模型（这类目标函数的优化可以采用代理模型辅助的昂贵优化算法，也叫贝叶斯优化或者EGO,三者其实说的都是一个东西，感兴趣的可以一起讨论啦，硕士期间主要就是研究这类算法）。
SA在实际应用中的主要流程如下（最小值问题）：
1：设置初始温度，终止温度，降火速率，迭代次数，生成初始解 x （算法迭代的起点）以及对应目标值f(x) .
2：以一定策略产生新解 \tilde{x} ，并计算目标函数增量 f(\tilde{x})
产生新解的策略与编码策略有关，比如（1）对于背包问题，直接进行 \left\{ 0,1,0,1,1,0\right\} 这种编码，1表示选择，0代表不选择，新解的产生可以通过随机交换两位置或者进行某位置的随机变异；而（2）对于TSP问题，比如进行 \left\{ 1,2,3,4,5 \right\} 这种编码，直接随机交换位置便可以实现新解的产生。待新解产生后，计算目标值的增量 \triangle=f(\tilde{x})-f(x) 。
3：判断新解是否被接受
这一步是SA算法的核心，一般采用Metropolis准则来判断解是否被接受，其接受概率 p的计算如下：
\begin{array}{c} p=\begin{cases} &1 {.......\bigtriangleup <0 } \\ &e^{-\frac{\bigtriangleup }{t}}....\bigtriangleup >0 \end{cases} \end{array}
当新解比原来的解好 (\Delta<0) 时，接受概率为1，即无条件接受新解，当新解比原来的解差时 (\Delta>0) ,这里并没有舍去它，而是按照一定概率接受它，防止陷入局部最优。此处 t 为当前的温度， 当前温度t 处于比较高的时候，接受概率 p 趋于1，因此当前更差的新解有更大概率被接受，进而可能跳出局部最优进行全局搜索。随着退火的进行，当前温度 t 减小到一个比较小的值时，接受概率趋近于0，不太能跳出当前搜索范围。因此整体而言，SA算法前期趋向于全局搜索，后期趋向于局部搜索，也就是论文里面常说的Exploration and Exploitation（探索与开发）。
4：新解替换当前解：代替变换部分，修正目标函数值即可。
整体的算法流程如图1下：
整体来说SA实现起来比较简单，在解决实际问题中，其实主要就是对优化问题的描述要清晰，编码方式选择好，然后套用框架即可，下面是一个关于求解背包问题（最大值）的Matlab程序，配合代码理解SA的原理会更好一些，有错误还请大家多多指出。
% 模拟退火实现背包问题（最大值问题）
%编制：mumusong 2020.06
clc;clear;close all
weight = [2,5,18,3,2,5,10,4,11,7,14,6]; %每个货物重量
price = [5,10,13,4,3,11,13,10,8,16,7,4]; %每个货物单价
weight_cons = 46; %背包重量约束
%初始化
alpha = 0.95; %退火系数
t_begin = 100; %最高温度
t_end = 2; %结束温度
t = t_begin; %当前退火温度
solution_new = ones(1,length(weight));
solution_current = zeros(1,length(weight));
value_current = 0;
value_best =0;
solution_best = solution_new;
counter = 0; %记录迭代次数
while t>t_end
counter = counter+1;
for i = 1:100 %马尔科夫长度
%生成新的解
index = randi([1,length(weight)],1,1);%随机生成1~10的一个整数
solution_new(1,index) =~solution_new(1,index);
% 判断是否满足约束
while sum(solution_new.*weight) > weight_cons
index = randi([1,length(weight)],1,1);
solution_new(1,index) =~solution_new(1,index);
end
% 判断是否接受当前解
value_new = sum(solution_new.*price);
probability = exp((value_new-value_current)/t); %当前概率
if (value_new>value_current)||probability>rand %新解更好或者劣解满足概率
value_current = value_new;
solution_current = solution_new;
else
solution_new = solution_current;
end
% 与最优解对比
if value_current>value_best
value_best = value_current;
solution_best = solution_current;
end
end
%保存每个过程的最优解
value_list(counter,:)= value_best;
solution_list(counter,:) = solution_best;
t = t*alpha;
end
%显示当前最优解
fprintf('最大价值：%f, 货物重量：%d\n',value_best,sum(solution_best.*weight));