模拟退火算法--北海数学建模课程笔记

xzzzzq

于 2024-08-12 14:27:06 发布

阅读量508

点赞数 18

文章标签：数学建模模拟退火算法笔记

本文链接：https://blog.csdn.net/weixin_61963212/article/details/141124497

版权

本文为北海的数模课程学习笔记，课程出自微信公众号：数学建模BOOM

模拟退火算法简介

模拟退火（Simulated Annealing, SA）是一种概率型优化算法，它受到金属退火过程的启发。在金属加工中，退火是将金属加热到一定温度然后缓慢冷却的过程，这样可以减少金属内部的缺陷，提高其整体性能。模拟退火算法将这一物理过程抽象为一种优化策略，用于解决数学上的优化问题，特别是那些具有复杂搜索空间的问题。

模拟退火算法的基本原理

初始状态：算法从一个随机解或者某个初始解开始
温度参数：设定一个初始温度 T，这个温度参数在算法过程中会逐渐降低。温度逐渐降低的过程即为退火的过程
领域搜索：在当前解的邻域内随机选择一个新的解。邻域是指与当前解相近的解的集合，例如在旅行商问题中，邻域可以是交换两个城市的位置。
接收准则：新解被接受为新的当前解的条件有两个：
(1) 如果新解比当前解更优，则总是接受新解。
(2) 如果新解比当前解差，以一定的概率接受新解，这个概率由当前温度决定，通常使用Boltzmann公式.
冷却过程：根据某个冷却方案降低温度，然后重复上述过程。这个可以自行设定，降温越缓慢，迭代次数越多。

特点

概率接受：相较于贪心算法而言，模拟退火允许算法以一定的概率接受较差的解，这有助于算法跳出局部最优解，探索更广阔的解空间。
温度控制：温度参数的设置和降低策略对算法性能有很大影响。
适用性广：适用于连续和离散的优化问题，特别是对于非线性、多峰值、多维数的问题。

优缺点

优点：简单易实现，不需要了解问题的具体结构，对初始解不敏感，能够找到全局最优解的概率较高。
缺点：缺点：收敛速度可能较慢，需要精心设计冷却方案，对于某些问题可能效率不高。

基于matlab软件内实现模拟退火算法解决TSP问题

问题重述

假如你是一位公司员工，你们公司在全国34个省会城市（包括直辖市，自治区首府，港澳台）都有一个自己的分公司，现在，你们老板想要去每一个城市考察一下当地子公司的业绩水平，于是将这个任务交给了你，现在从北京出发，如何设计出行路线，使得出行距离最短，并且每个城市只去一次，最终回到北京，已知每一个城市的经纬度坐标。

问题分析

很明显这是一个经典的旅行商（TSP）问题，在传统情况下，如果采取枚举法，那边有33!(33的阶乘)种可能的情况，用计算机求解太过于复杂，此时呢可以考虑使用模拟退火方法。

原理和思路

由问题可知，本题的所有可行解是一个很庞大的解空间，数量级在10的36次方。
本题的解空间就是固定起点和终点的所有可行路线的集合。

以北京为起点城市，记为x₁
遍历剩下的所有城市序号为x₂,x₃,…,直到34号城市x₃₄
最终回到北京，记为x₃₅。

我们的目标是在不重复的经过每一座城市，使得出行的总路程最短
从1号城市x₁出发，经过剩下的33个城市，最终回到35号城市北京x₃₅。
在这个过程种，我们可以将第x_i和第x_i+1号城市之间的距离设置为d_xixi+1。
那我们的目标函数便可设置如下：
目标函数表达式

求解思路

初始化
设定一个初始温度T₀和问题的一个初始解x(0)
因为这个问题存在很多解，所以随便选取一个作为初始解即可
本题可以直接选取城市序号排列1-34的顺序作为初始解，但是为了更好的应对其他问题，本文采取蒙特卡罗法进行初始解求解（在代码中具体实现）
随机产生一个新解
在当前温度T_i下随机产生一个新解x^’，产生新解的方式不唯一，确保随机即可
假设在上一步中，我们得到的解为

这里的u,v为随机选取的序号，我们将u-v这段进行逆序，即得到新解如下
计算目标函数差值
得到的随机解x^’和我们的当前解x(0)之间存在路程差，这一部分的计算公式为
是否接受新解
接受的准则为
马尔可夫过程
在当前温度下，重复进行第2，3，4步，选取出当前温度下的最优解
退火过程
选定降温系数a，求得降温后的新温度为T_i+1=aT_i，这里的a=0.999，降温较为缓慢
在新的温度下，继续重复2，3，4，5步，寻找新温度下的最优解
结束条件
当温度足够小，达到终结温度时，停止迭代，输出结果。
当退火速度足够慢，在每个温度下寻找最优解的次数足够多的时候，最终得到的全局最优解的概率越大。

代码实现

下面结合matlab代码，对该问题进行求解
第一步，导入数据，并且对所有城市之间的距离进行求解

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

第二步，利用蒙特卡洛法，求初始解

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

第三步，模拟退火的主循环过程

e=0.01^30;  % 终结温度
alpha=0.999;  % 降温系数
T=1000;  % 初始温度
markov=100;  % 马尔可夫过程的循环次数

accept=0;  % 接受新解的次数
rand_accept=0;  % 在概率允许下接收次解的次数
refuse=0;  % 拒绝当前解的次数
n2=0;  % 初始化降温次数
num2=[];  % 初始化每一个温度下的最短距离的值
tic;  % 启动计时器
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;  % 进行每一次的降温
    n2=n2+1;  % 统计所进行降温的次数
    num2(n2)=lenth;  % 对每个温度下的最优解做以统计
end
toc;  % 显示计时器时间

第四步结果可视化

num1=1:n2;  % 作为x轴
plot(num1,num2);
xlabel('次数')
ylabel('最短距离')
title('最短距离变化图');
path(35)=1;     % 1号和35号都是北京，把35改成1，方便画图
plot(city(path ,1), city(path ,2),'o-');   
xlabel('东经');
ylabel('北纬');
title('最短距离连线图');
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)]); 
end