模拟退火算法应用——求解TSP问题
仅作自己学习使用
一、问题
旅行商问题(TSP) 是要求从一个城市出发，依次访问研究区所有的城市，并且只访问一次不能走回头路，最后回到起点，求一个使得总的周游路径最短的城市访问顺序。
采用模拟退火算法求解TSP问题，很自然的想到退火的目标函数（优化函数）应该就是总的周游距离。那么在算法中如何体现呢？那就是把城市的坐标放在一个n×2的矩阵中，矩阵中存放城市的顺序就是依次周游城市的路径，所以在求解过程中会不断的产生新的更优解（周游顺序，在算法中体现就是城市坐标的存放顺序），有了这个关键的思路就很好解决了。
二、Matlab代码
clear
clc
T1 = cputime;
C = [
% 各个城市坐标
39.91, 116.39; % 北京
31.22, 121.48; % 上海
23.13, 113.27; % 广州
22.54, 114.06; % 深圳
30.67, 104.06; % 成都
34.27, 108.93; % 西安
31.98, 118.75; % 南京
39.92, 116.36; % 天津
28.71, 115.83; % 南昌
45.75, 126.63; % 哈尔滨
36.07, 120.38; % 青岛
38.04, 114.48; % 石家庄
29.59, 106.54; % 重庆
26.08, 119.30; % 福州
30.25, 120.16; % 杭州
28.19, 112.97; % 长沙
25.03, 102.73; % 昆明
35.68, 139.76; % 东京
37.56, 126.97; % 首尔
1.35, 103.82; % 新加坡
13.41, 103.86; % 金边
21.03, 105.85; % 河内
3.14, 101.69; % 吉隆坡
39.90, 32.85; % 安卡拉
37.97, 23.73; % 雅典
38.71, -9.14; % 里斯本
41.89, 12.50; % 罗马
52.52, 13.41; % 柏林
55.75, 37.62; % 莫斯科
48.86, 2.35; % 巴黎
];
n = length(C); % 获取城市的个数
T = 100 * n; % 初始温度
L = 10; % 马尔可夫链长度
K = 0.986; % 降温系数
%% 构建城市坐标结构体
city = struct([]);
for i = 1:n
city(i).x = C(i,1); % 经度
city(i).y = C(i,2); % 纬度
end
%% 开始退火
% 统计迭代次数
count = 1;
% 计算每次迭代后的总距离（第一次就是初始时，按照坐标的顺序计算的距离）
Dist(count) = GetDist(city,n);
figure(1)
% 当温度无限趋于0度时停止迭代
while T > 0.01
% 每次降温 均进行多次迭代
for i = 1:L
% 计算原路线周游距离
len1 = GetDist(city,n);
% 产生随机扰动(随机交换两个城市的坐标)
p1 = floor(1 + n * rand()); % rand函数产生一个0，1之间均匀分布的实数，包含0但不包含1
p2 = floor(1 + n * rand()); % 因此这个表达式可以产生一个从1到n的随机数
while (p1 == p2)
p1 = floor(1 + n * rand());
p2 = floor(1 + n * rand());
end
temp_city = city;
% 交换第P1个城市和第P2个城市的坐标
temp = temp_city(p1);
temp_city(p1) = temp_city(p2);
temp_city(p2) = temp;
% 计算新路线的周游距离
len2 = GetDist(temp_city,n);
% 新、老路线的差值（相当于能量）
delta = len2 - len1;
if(delta<0)
% 新路线的评估函数更小（记住，模拟退火算法相当于是一个求函数极小值的算法）
city = temp_city; % 更新原路线（变量里存放城市的顺序也就是访问城市的顺序）
else
% Metropolis接受准则（概率选择更差的解）
if exp((len1-len2)/T) > rand()
% 记住这个概率的公式，指数部分一定是要个负数，概率的值不可能超过1
city = temp_city;
end
end
end
% 本次迭代结束，统计迭代次数加1
count = count + 1;
% 将本次迭代的最优解放在len中
Dist(count) = GetDist(city,n);
%% 本次退火结束，降温
T = T * K;
% 按照新的城市的顺序，把这些城市画出来
for i = 1: n-1
plot([city(i).x,city(i+1).x],[city(i).y,city(i+1).y],'bo-');
hold on;
end
plot([city(n).x,city(1).x],[city(n).y,city(1).y],'ro-');
title(['优化最短距离：', num2str(Dist(count))]);
hold off
pause(0.005); % 动态显示出每次的搜索结果
end
T2 = cputime;
figure(2)
plot(Dist,LineWidth=2)
xlabel("迭代次数")
ylabel("目标函数值")
title("适应度进化曲线","搜索时间："+(T2-T1)+" s")
%% 评估函数
function result = GetDist(city,n)
% 计算总的周游路径长度(评估函数)
% city是各个城市的坐标
result = 0;
for i = 1:n-1
result = result + sqrt((city(i).x - city(i+1).x)^2 + (city(i).y - city(i+1).y)^2);
end
result = result + sqrt((city(n).x - city(1).x)^2 + (city(n).y - city(1).y)^2);
end
三、效果
四、问题
大家可以试一试更多的城市，当有很多城市的坐标相差不大时，在最后的搜索结果中，会出现一个非常奇怪的问题，就是在周游图中，有些城市消失了，检查存放城市的city结构体，是存放着这些坐标的，这里如果有知道的朋友还请多多批评指教，我将及时改正。