美赛BOOM数学建模5-1模拟退火

注明：本文根据数学建模BOOM网课简单整理，自用

❑模型简介

❑ 模拟退火的“实例”

• 齐白石原是木匠，快30岁才正式学画。直到57岁以后画风开始转变，才真正有所成就

• 鲁迅原学医，25岁弃医从文，但37岁时才写出第一篇白话文小说《狂人日记》

• 项羽24岁巨鹿之战成名，26岁称王却走错了路子，三十出头就……

❑ 模拟退火的基本思想，就是走出舒适圈，多去“试一试”，万一成了呢？

• 舒适区：当前处于局部最优解（鲁迅那个年代能留学的都是佼佼者）

• 试一试：随机试探新解，有更好的解就直接选新解，没更好的则以一定概率选择新解

• 万一成了：找到了更优解甚至最优解，青史留名（载入语文/历史课本）

• 也可能没成：求的新解反而更差了（项羽没有去统一而是封王结果……）

• 一次没成，多试几次：继续随机试探（齐白石也不是转行后立马成名的）

❑贪心算法失效——陷入局部最优

❑ 从爬山说起

• 假如北海想要在日落前爬上一座山的最高峰（求最优解）

• 但山中云雾缭绕、视野受限，只能看到当前位置附近一定范围内的情况

• 贪心算法：在视野内选取一个点，如果更高，就过去；否则就不去

• 直到视野内没有更高的点，则输出当前解作为最终解（视频中有板书过程）

• 显然，当北海走到某个小山峰时，四周看不到更高的点，就不会走下山去、寻找更高的峰

• 从而陷入了局部最优解

❑ 模拟退火思想

• 那么在可见范围内，随机选择一点

• 1.如果该点比当前位置更高，就直接去该点（优化）

• 2.如果该点更低，那么就多掷几次硬币，结合该点与当前点的高度差决定去不去（一定概率）

• 在刚才的局部最优解的峰，会有一定概率走下了当前山峰，从而发现另一个山峰的上坡

• 从而就有可能走上新的更高峰

❑ 算法关键点

• 时间有限，需设置终止条件（“退火”过程中温度不断降低，剩余时间不多时就别乱跑）

• 且算法需要做到：若当前所处的山峰越高，前往低点的概率越低（概率的表达式）

• 这样才能使求得最优解的概率更大

• 关键在于，“概率”怎么求

❑适用赛题

❑ 可行解过多、 NP-hard问题

❑ NP-hard问题

• 旅行推销员问题（Travelling salesman problem, TSP）：给定一系列城市和每对城市之间的距离，求解访问每一座城市一次并回到起始城市的最短回路

• TSP问题是组合优化中的一个NP-hard问题，该问题的可行解是所有顶点的全排列，随着顶点数的增加，会产生组合爆炸，传统算法难以求解

• 此时就适合用模拟退火等启发式算法求近似解

❑ 适用赛题特点

• 有些规划类问题的可行解过多，传统算法运算时间过长

• 能够列出可行解的任一排列（例如题目要求最短回路，起码能先写出个回路作为初始化）

❑ 国赛题目：高等教育学费标准

• 假设高校所需1000万经费，涉及多种筹款路径，理论上的筹款组合方案过多，若暴力求解可能费时几十小时，此时就适合模拟退火求解

❑典型例题

❑ 遍历全国省会

• 已知全国34个省会城市（包括直辖市、自治区首府和港澳台）的经纬度坐标（第一个为北京）

• 现在需要从北京出发，到所有城市视察

• 要求每个城市只能到达一次，并最终回到北京

• 求视察路线方案，使得总路径最短

❑ 问题分析

• 若把所有可能的路线方案都列出来再求最小值，则有33！（33的阶乘）种方案

• 这是1036数量级，运算时间会很长

• 此时可考虑使用模拟退火

• 本期课程的原理讲解内容对于部分同学来说较难理解，应结合代码理解

❑原理与求解思路

❑ 所有可行解

• 本题的解空间就是固定起点和终点的所有可行路线的集合

• 题目要求以北京作为起点和终点，那么设起点北京为1号城市𝑥1

• 余下需要遍历的33个城市依次设为2，3，……，34号城市𝑥34

• 最后还需回到北京，设终点北京为35号城市𝑥35

• 因为起点与终点固定，则共有33！个可行解

❑ 目标函数：总路程最短

• 从1号𝑥1出发遍历34个城市最后回到1号𝑥35

• 第𝑖号城市𝑥𝑖与第𝑖 + 1号城市𝑥𝑖+1之间的路径长度设为𝑑𝑥𝑖𝑥𝑖+1

• 则目标函数：

❑ 1、初始化

• 设定一个初始温度𝑇0和问题的一个初始解𝑥(0)

• 有那么多种可行解，可以随便选一个作为初始解

• 例如本题，可以把原始数据的排列直接作为初始解，解为𝑥1𝑥2 𝑥3…… 𝑥35

• 或利用蒙特卡罗法求一个较好的初始解（课程第四部分代码求解会讲）

❑ 2、随机产生新解

• 在当前温度𝑇𝑖下随机求一个解𝑥′，制定产生新解的准则（准则不唯一，能确保随机即可）

• 假设上一步的解为：

• 随机选序号𝑢, 𝑣，将𝑢到𝑣的这部分转为逆序，得到解：

❑ 3、计算目标函数差值

• 随机得到的解𝑥′和当前解𝑥(𝑖)，两种方案总路径的差值记为𝛥𝑓 ：

❑ 4、是否接受新解

• 接受准则：

• 表达式中判断准则𝑒𝑥𝑝( − 𝛥𝑓/𝑇𝑖)是由热力学得到的灵感，后面补充讲解

• 题目求总路径最短，如果𝛥𝑓 < 0意味解𝑥′的总路径更短，则接受𝑥′作为新解， 𝑥(𝑖)= 𝑥′

• 𝛥𝑓 > 0意味解𝑥′总路径更长不如原解𝑥(𝑖)，但为防止陷入局部最优解，仍以一定概率接受𝑥′

• 在代码实现上，if 𝑒𝑥𝑝( − 𝛥𝑓/𝑇𝑖) ≥ rand 则接受𝑥′作为新解，否则保留解𝑥(𝑖)

• 显然， 𝛥𝑓越小，意味着p越大，那么接受𝑥′的概率就更大

❑ 5、马尔科夫过程

• 在当前温度T𝑖下，重复2、3、4步

• 温度Ti不变时，判断准则中的e𝑥𝑝( −𝛥𝑓/𝑇i) 由𝛥𝑓 决定，而𝛥𝑓 由随机解x′和当前解x(𝑖)决定

• 所以新解x（i+1）只与𝑥（i）有关，而与更早的x(i-1),x(i-2)...x(0)无关

• 显然这是个马尔科夫过程,x′在原解x(𝑖)的邻域中符合均匀分布*

❑ 6、退火过程

• 选定降温系数𝛼，求得新温度T𝑖+1 = 𝛼𝑇𝑖，此处设为0.999

• 再重复2、3、4、5步；然后继续降低温度……

❑ 7、结束条件

• 直到温度足够小，设终止温度为e = 10^−30，当𝑇 < 𝑒时终止迭代，输出最终解

• 退火过程足够慢、每个温度下寻找新解的次数足够多，则最终解是全局最优解的概率越大*

❑ *理论上模拟退火可以找到全局最优（该部分不要求掌握）

• 马尔科夫过程随机解𝑥′只与x(i)有关，所以x'在原解x(𝑖)的邻域中符合均匀分布

• 而x′被接受的状态由上页的概率表达式p 所决定

• 那么经过有限次的重复2、3、4步，在温度Ti下解的平衡态x𝑖满足玻尔兹曼分布：

• 平衡态Xi代表温度Ti下任意一种可能出现的解，S为所有可能出现的解的集合

• 式中f(i)代表该解下的总路径（目标函数值）

• 接着降低温度𝑇𝑖+1< T𝑖，在温度𝑇𝑖+1下有限次的重复2、3、4步

• 整个优化过程就是不断寻找新解和缓慢降温的交替过程，那么温度趋于0时得到最终解

• 重点在于，温度趋于0时得到的最终解，是不是接近于全局最优解？

• 上一页提到，温度𝑇𝑖 下解的平衡态𝑥𝑖 满足玻尔兹曼分布，根据玻尔兹曼分布和状态转移概率：(热力学统计物理，背景概念很多，此处只给结论)

• 从玻尔兹曼分布可知，温度T → 0时某个值大于0的状态（函数值）f(𝑥𝑖)的分布概率趋于0

• 也就意味着只有f(𝑥𝑖)尽可能趋于0的状态才会存在

• 对应现实问题，求得的最终解中只有fmin，也就是最优解（本题的最短路径）

• 从状态转移概率公式可知，当温度𝑇 → 0， 𝛥𝑓 ≥ 0时转移到更高函数值的概率趋近于0

• 也就意味着温度极低时，转移到更差解的概率越小、会以更大概率停留在当前解

• 结论：当寻找新解次数足够多、降低温度足够慢，最终会以接近于1的概率求得全局最优解

• 证明：根据玻尔兹曼分布和状态转移概率，求最终从一个极低的温度转移到𝑇 = 0后的概率

分布，结合两个公式可得：

❑代码求解

• 重点在于，如何在代码中实现课程第三部分原理讲解中的各种计算

• 求解结果每次都不一样

已知全国34个省会城市（包括直辖市、自治区首府和港澳台）的经纬度坐标（第一个为北京）；现在需要从北京出发，到所有城市视察，要求每个城市只能到达一次，并最终回到北京。求视察路线方案，使得总路径最短。

1、导入数据

因为题目提供的数据是经纬度，注意第一列是经度、第二列是纬度，写函数时不要写错了

求两个城市间距离，是求在地球上的弧线长度，也就是求一个扇形的弧长。

可使用distance函数求球体上两点间距离

d(i,j) = distance(lat1,lon1,lat2,lon2,radius);其中lat1、lon1是第一个点的纬度和经度，lat2,lon2是第二个点的纬度和经度，radius是球体的半径

clc, clear, close all
%数据在文件cities.xlsx中，注意要放在与本文件同一文件夹下
city = table2array(readtable('citys.xlsx','Range','B2:C35')); %table2array改变数据类型
n = size(city,1);      %城市距离初始化                                                                    
d = zeros(n,n+1);      %35号只会作为终点          
for i = 1:n
    for j = 1:n
            d(i,j)= distance(city(i,2),city(i,1),city(j,2),city(j,1),6371); % distance求圆心角的角度.第一个点的纬度、第一个点的经度
    end    
end
%各城市到35号即北京的距离作为第35列
for i=1:n
    d(i,35)= distance(city(i,2),city(i,1),city(1,2),city(1,1),6371);
end

2、蒙特卡罗法求初始解

之前讲过蒙特卡罗法。本题总共34个城市，从1号城市北京出发，遍历完所有34个城市后，返回北京，设返回的北京为第35号。

每次随机选一种路径方案，出北京外的其他33个城市，利用randperm(33)为33个从1到33的整数随机排列；那么1+randperm(33)就是2到34的整数随机排列。

path=[];
lenth=inf; %总路径及长度初始化，inf正无穷
for j=1:1000  %求较好的初始解，随机求1000种方案，挑出最好的作为初始方案
    temp_path=[1 1+randperm(33) 35];    % 当前解（方案）
    temp_lenth=0;   % 当前方案的总路径
    % 求该方案下总路径长度temp_lenth
    for i=1:34
        temp_lenth=temp_lenth+d(temp_path(i),temp_path(i+1));
    end
    % 如果该方案下总路径长度temp小于所记录的当前最短总路径长度long（初始为正无穷）
    if temp_lenth<lenth
        path=temp_path; lenth=temp_lenth;      % 将该路径方案temp_path记为最短路径方案path，将该方案的长度temp_lenth记为最短路径长度lenth
    end
    % 如此循环1000次，就从1000个随机方案里挑选出最优的方案作为初始方案，用于后面的模拟退火
end

3、模拟退火

while循环实现温度不断降低，温度初试为1，每次新的温度是之前的0.999倍；当温度低于10^(-30)时终止循环。

while循环内嵌套for循环，实现固定温度下的马尔科夫过程，即不断求新解、并判断是否接受。

判断是否接受新解，就是代码中的if和elseif两种情况

e=0.1^30;alpha=0.999;T=1000;markov=1000;   %这些参数都是可以改的
% for k=1:L  %退火过程
accept=0;rand_accept=0;refuse=0;
tic
while T>e
    for t=1:markov
        %新解随机选序号𝑢, 𝑣，将𝑢到𝑣的这部分转为逆序作为新解
        c=2+floor(33*rand(1,2));  %产生新解；floor向下取整，得到两个2到34的随机整数，（1，2）一行两列
        c=sort(c);  %随机选的两个点升序排序，用于接下来计算
        u=c(1);v=c(2);  %模型中的随机选的两个点u和v，u是序号小的那一个
        %计算目标函数值的增量
        df=d(path(u-1),path(v))+d(path(u),path(v+1))-...
            d(path(u-1),path(u))-d(path(v),path(v+1));
        if df<0 %接受准则
            path=[path(1:u-1),path(v:-1:u),path(v+1:35)]; %新路径u到v逆序
            lenth=lenth+df; 
            accept=accept+1;
        elseif exp(-df/T)>=rand   %rand产生0到1的随机数；不等号左边与温度T有关
            path=[path(1:u-1),path(v:-1:u),path(v+1:35)]; 
            lenth=lenth+df;
            rand_accept=rand_accept+1;
        else
            refuse=refuse+1;
        end
    end
toc
    T=T*alpha;
end
path(35)=1;     % 1号和35号都是北京，把35改成1，方便画图
plot(city(path ,1), city(path ,2),'o-');
    disp('最短路程：')
    disp(lenth)
    disp('直接接受新解次数：')
    disp(accept)
    disp('接受更差的随机解次数：')
    disp(rand_accept)
    disp('不接受随机解次数：')
    disp(refuse)    
for i = 1:n
    %对每个城市进行标号
    text(city(i,1),city(i,2),['   ' num2str(i)]); %number to string
end
xlabel('东经')
ylabel('北纬')