模拟退火——和随机梯度下降SGD有异曲同工之妙（附旅行推销员问题MATLAB求解讲解）

MikingG

于 2023-09-02 21:42:56 发布

阅读量470

点赞数 1

分类专栏：从零开始学数学建模学习笔记文章标签：算法数学建模

本文链接：https://blog.csdn.net/weixin_64123373/article/details/132461647

版权

从零开始学数学建模学习笔记专栏收录该内容

24 篇文章 5 订阅

订阅专栏

模拟退火（Simulated Annealing, SA）是一种启发式的随机搜索算法，其灵感来源于固体物理中的退火过程。在固体物理学中，物质通过加热然后慢慢冷却来减少其晶体结构的缺陷，达到一个能量较低的稳定状态。模拟退火算法借鉴了这一概念，旨在找到一个给定问题的近似全局最优解。

模拟退火和随机梯度下降SGD有异曲同工之妙

随机性：SGD 和模拟退火都包含了一定程度的随机性。在SGD中，每次更新模型参数时都是基于一小批次（mini-batch）的随机样本，而不是整个数据集。在模拟退火中，接受新解的概率也受到随机性的影响。这种随机性可以帮助算法跳出局部最小值，从而更有可能找到全局最优解。
逐步减小步长：在SGD中，学习率（或步长）逐渐减小，以便在接近最优解时更加稳定。而在模拟退火中，温度参数也逐渐降低，从而减小了接受新解的概率，使算法趋于稳定。
调整参数：SGD中的学习率和模拟退火中的温度参数都需要精心选择和调整，以确保算法能够有效地找到最优解。这种参数调整过程需要一定的经验和实验。
局部搜索和全局搜索：SGD通常用于训练神经网络等深度学习模型，用于局部搜索参数空间，找到模型的局部最优解。而模拟退火通常用于更通用的优化问题，可以在参数空间中进行更广泛的搜索，以找到全局最优解。

虽然SGD和模拟退火有一些相似之处，但它们也有明显的区别。SGD主要用于优化目标函数的导数，通常应用于机器学习和深度学习中，而模拟退火可以用于更广泛的优化问题，不一定需要目标函数的导数信息。

模型简介

基本思想：就是走出舒适圈（处于 局部最优 解），多去“试一试”（随机试探新解），万一成了呢？

也就是说，在可见范围内，随机选择一点：

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

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

关键在于

时间有限：“退火”过程中温度不断降低，温度可以看成是时间。
合理概率：当前的高度越高，前往低点的概率就越低。

这其实是个马尔科夫过程，因为新解只与当前解有关。

适用赛题

可行解过多、NP难问题。

典型的就是旅行推销员问题（ Travelling salesman problem, TSP ）。

还有多种组合的方案。比如要筹1000万经费，涉及多种筹款路径。

典型例题

题目描述

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

原理与求解思路

目标函数 $\min f(x_1,x_2,...,x_{35})= \sum_{i=1}^{34} d_{x_i x_{x+1}}$

1.初始化

初始温度 $T_0$ 和初始解 $x(0)$

初始解求法：原始数据的某种排列作为初始解/蒙特卡洛法求一个较好的初始解

2.随机产生新解

准则不唯一，能确保随机即可

3.计算目标函数差值

这种部分逆序的中间路程不会变，只有头尾会变。

差值就是：新头尾-旧头尾。

4.是否接受新解

如果路程更短，也就是差值为负数，肯定接受。如果路程更长，在一开始（温度高）时接受的概率大，后来（温度低）时接受的概率小。

5.退火过程

$T_{i+1}= \alpha T_i$ ,降温系数 $\alpha$ 可自己设置，比如0.999

6.结束条件

温度低于终止温度e。

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

代码求解

1、导入数据

因为题目提供的数据是经纬度，注意第一列是经度、第二列是纬度，写函数时不要写错了
求两个城市间距离，是求在地球上的弧线长度，也就是求一个扇形的弧长。

可使用distance函数求球体上两点间距离，使用distance前需要安装Mapping Toolbox

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'));
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; %总路径及长度初始化
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=1;markov=1;   %这些参数都是可以改的

% for k=1:L  %退火过程
accept=0;rand_accept=0;refuse=0;
while T>e
    for t=1:markov
        %新解随机选序号𝑢, 𝑣，将𝑢到𝑣的这部分转为逆序作为新解
        c=2+floor(33*rand(1,2));  %产生新解；floor向下取整，得到两个2到34的随机整数
        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
    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('北纬')