这两天宅在家里闲着是在是没什么意思，就总结一下以前数模时自己学到的一些算法吧，都是一些非常非常基础的算法，就作为学习笔记总结一下学到的最最干货的东西吧，以后学习科研都能用得上，也希望对各位学弟学妹们有所帮助吧，刚好2020MCM/ICM马上快开始了（今年没有参赛），希望对各位参加建模的同学有所帮助。关于数模的一些其他的干货（建模流程、写作、时间安排等等）我都放到这两个视频里面了，感兴趣的同学可以看一下：av73775318(美赛)、av66218067(国赛/省赛/校赛)。
模拟退火算法（SA）
说算法之前先介绍主要用模拟退火解决的问题吧
TSP问题（旅行商问题）：如下图所示，每个圈代表一个城市，各城市之间的坐标固定，居民可以通过两城市之间的连线由一个城市去另一个城市。有一个商人从城市1出发，他要走遍其他城市且只经过这些城市一次，最后回到城市1,问他怎么规划路线，可以使他走过的路线最短。
小偷背包问题(0-1背包问题）：有一个小偷行窃时发现共10件物品，质量分别为：10,8,5,9,7,4,3,2,5,6（kg），对应价值分别为15,6,10,14,8,9,6,8,3,8（价值单位），但是该小偷只能通过背包带走物品，该背包能承受的最大质量为40kg，问他需要带走那些物品，可以使他的利益最大化。
下面简单总结一下模拟退火的干货（详细介绍请参考教科书与相关论坛）：
模拟退火算法模拟的是退火过程（将金属试件加热到一定温度后在另一较低温度下冷却，如果冷却速度足够低，则冷却过程中任一温度下该试件都处于热平衡状态，试件最终会达到热稳定状态，此时试件的内能最低）(应用时将控制参数视为温度，将目标函数值视为内能，控制参数足够大且在循环过程中不断变小，目标函数值不断朝向“最小”的方向变化，直至最终达到稳定为止)。模拟退火算法是一种概率算法（表现在“内能”未达到“最低”而“稳定”时对“新内能”的选取），用以在大空间内寻找最优解，能够求出全局最优解（而不是陷入局部最优，这是梯度下降法求大空间最优解的缺点）。
模拟退火算法的基本步骤简述如下：
1.规定退火开始时的初始温度，随机生成一个初始解（解1）（该初始解是全局的初始解），计算对应的目标函数值（值1）
2.令温度降低，对上一步的全局初始解施加扰动（注意：如何施加扰动是模拟退火算法的关键），得到新解（解2）及新目标函数值（值2）
3.若值2较值1小，则接受解2，作为新的解1，否则以一定概率接受解1，作为新的解1
4.不断迭代上述过程，直至温度不断降到最低
下面给出TSP问题的模拟退火MATLAB程序(题目是2019MCM/ICM的D题，将卢浮宫的逃生点或出入口化成34个点，用优化算法得到各点之间的相对二维坐标，再将负责人遍历各点的问题转化为tsp问题)：
已知各点及其二维坐标:(点序号 x坐标 y坐标；)
1 -12.0 -180.0; 2 8.0 -180.0; 3 18.0 -180.0; 4 28.0 -180.0; 5 8.0 -180.0; 6 -12.0 -80.0; 7 8.0 -80.0; 8 18.0 -80.0; 9 28.0 -80.0; 10 8.0 -90.0; 11 -32.0 28.0; 12 -12.0 20.0; 13 -2.0 28.0; 14 8.0 20.0; 15 18.0 20.0; 16 28.0 20.0; 17 28.0 10.0; 18 8.0 0.0; 19 18.0 0.0; 20 28.0 0.0; 21 48.0 20.0; 22 48.0 30.0; 23 -32.0 128.0; 24 -12.0 120.0; 25 8.0 120.0; 26 18.0 120.0; 27 8.0 100.0; 28 28.0 100.0; 29 48.0 130.0; 30 38.0 120.0; 31 8.0 220.0; 32 18.0 220.0; 33 8.0 200.0; 34 38.0 228.0;
代码如下（参考卓金武编 MATLAB在数学建模中的应用 北航出版社）：
clc clear a = 0.99; t0 = 97; tf = 3; t = t0; Markov_length = 10000; coordinates = [ 1 -12.0 -180.0; 2 8.0 -180.0; 3 18.0 -180.0; 4 28.0 -180.0; 5 8.0 -180.0; 6 -12.0 -80.0; 7 8.0 -80.0; 8 18.0 -80.0; 9 28.0 -80.0; 10 8.0 -90.0; 11 -32.0 28.0; 12 -12.0 20.0; 13 -2.0 28.0; 14 8.0 20.0; 15 18.0 20.0; 16 28.0 20.0; 17 28.0 10.0; 18 8.0 0.0; 19 18.0 0.0; 20 28.0 0.0; 21 48.0 20.0; 22 48.0 30.0; 23 -32.0 128.0; 24 -12.0 120.0; 25 8.0 120.0; 26 18.0 120.0; 27 8.0 100.0; 28 28.0 100.0; 29 48.0 130.0; 30 38.0 120.0; 31 8.0 220.0; 32 18.0 220.0; 33 8.0 200.0; 34 38.0 228.0; ]; coordinates(:,1) = []; amount = size(coordinates,1); dist_matrix = zeros(amount, amount); coor_x_tmp1 = coordinates(:,1) * ones(1,amount); coor_x_tmp2 = coor_x_tmp1'; coor_y_tmp1 = coordinates(:,2) * ones(1,amount); coor_y_tmp2 = coor_y_tmp1'; dist_matrix = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + (coor_y_tmp1-coor_y_tmp2).^2); sol_new = 1:amount; % sol_new是每次产生的新解；sol_current是当前解；sol_best是冷却中的最好解； E_current = inf;E_best = inf; % E_current是当前解对应的回路距离； % E_new是新解的回路距离； % E_best是最优解的 sol_current = sol_new; sol_best = sol_new; p = 1; while t>=tf for r=1:Markov_length% 产生随机扰动 if (rand < 0.5) % 随机决定是进行两交换还是三交换 % 两交换 ind1 = 0; ind2 = 0; while (ind1 == ind2) ind1 = ceil(rand.*amount); ind2 = ceil(rand.*amount); end tmp1 = sol_new(ind1); sol_new(ind1) = sol_new(ind2); sol_new(ind2) = tmp1; else % 三交换 ind1 = 0; ind2 = 0; ind3 = 0; while (ind1 == ind2) || (ind1 == ind3) || (ind2 == ind3) || (abs(ind1-ind2) == 1) ind1 = ceil(rand.*amount); ind2 = ceil(rand.*amount); ind3 = ceil(rand.*amount); end tmp1 = ind1;tmp2 = ind2;tmp3 = ind3; % 确保ind1 < ind2 < ind3 if (ind1 < ind2) && (ind2 < ind3); elseif (ind1 < ind3) && (ind3 < ind2) ind2 = tmp3;ind3 = tmp2; elseif (ind2 < ind1) && (ind1 < ind3) ind1 = tmp2;ind2 = tmp1; elseif (ind2 < ind3) && (ind3 < ind1) ind1 = tmp2;ind2 = tmp3; ind3 = tmp1; elseif (ind3 < ind1) && (ind1 < ind2) ind1 = tmp3;ind2 = tmp1; ind3 = tmp2; elseif (ind3 < ind2) && (ind2 < ind1) ind1 = tmp3;ind2 = tmp2; ind3 = tmp1; end tmplist1 = sol_new((ind1+1):(ind2-1)); sol_new((ind1+1):(ind1+ind3-ind2+1)) = sol_new((ind2):(ind3)); sol_new((ind1+ind3-ind2+2):ind3) = tmplist1; end %检查是否满足约束 % 计算目标函数值（即内能） 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) E_current = E_new; sol_current = sol_new; else sol_new = sol_current; end end end t=t.*a; end disp('最优解为：') disp(sol_best) disp('最短距离：') disp(E_best)
路线结果：
16 15 14 13 12 11 27 26 25 24 23 33 31 32 34 29 30 28 22 21 9 8 7 10 4 3 2 5 1 6 18 19 20 17
这篇先写到这里，有什么意见或建议可以留言，谢谢大家
参考：卓金武等.MATLAB在数学建模中的应用.北京航空航天大学出版社
