模拟退火算法求解TSP问题

本章会介绍在matlab上运行模拟退火算法求解TSP问题

目录

• 一、问题提出
• 二、导入数据与生成距离矩阵
  • 1.导入数据与变量初始化
  • 2.生成距离矩阵
• 三、进行模拟退火
  • 1.参数设置
  • 2.进入循环
    • （1）最外层循环
    • （2）生成随机数，替换解序列
    • （3）计算目标函数值的增量
    • （4）进入判断语句--是否接受结果

一、问题提出

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

二、导入数据与生成距离矩阵

1.导入数据与变量初始化

readtable可以读取excel文件，但读取后是一个table结构，需要用table2array将其转换成double结构
这里我们读取的是34个城市的经纬度坐标，是一个34行，2列的矩阵，但最后我们要返回北京，因此把北京定义为第35个城市。
又由于我们不会从第35个城市出发，因此在距离矩阵中只定义了第i个城市到第35个城市的距离，没有定义从第35个城市到第i个城市的距离。也就是没有35行但有35列

clc, clear, close all
%数据在文件cities.xlsx中，注意要放在与本文件同一文件夹下
city = table2array(readtable('citys.xlsx','Range','B2:C35'));
n = size(city,1);      %城市距离初始化                                                                    
d = zeros(n,n+1);

2.生成距离矩阵

distance函数可以通过输入两点的经纬度和圆的半径求两点的距离

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

三、进行模拟退火

在进行模拟退火之前，我们可以用蒙特卡洛模拟求一个解作为初始解，这里我们略去这一步骤

1.参数设置

退火温度为e；降温速度为alpha；T为初始温度；markov为在某一温度下的循环次数
e=0.1^30;alpha=0.999;T=1;markov=10;

2.进入循环

我们来分开讨论这个循环的代码，先附上总代码

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

（1）最外层循环

给定退出循环的条件：当温度过低，低于我们所给的推出阈值
在每经过一轮迭代后，温度都会降低

while T>e
    T=T*alpha;
end

（2）生成随机数，替换解序列

在该温度下，循环markov次，这里的参数可以改，但是不建议改成很大，会极大地增加运行时间
由于第一个城市和最后一个城市都是给定的，因此我们只需要生成两个从2-34的整数即可，33*rand是一个定义在0-33的有理数，向下取整后就得到了一个定义在0-32的整数（因为rand是开区间），此时再加上2就得到了一个从2-34的整数，再将其进行升序排序，就得到了替换序列的节点

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是序号小的那一个

（3）计算目标函数值的增量

按公式计算目标函数的增量
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));

（4）进入判断语句--是否接受结果

在接受准则后，通过向量连接的方式更新路径，再在原先所求的距离上加上我们所得的df
不断地迭代该步骤，得到最后的全局最优解

 if df<0 %接受准则
            path=[path(1:u-1),path(v:-1:u),path(v+1:35)]; %新路径u到v逆序
            lenth=lenth+df; 
        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;
        end