%% 本程序旨在基于模拟退火算法计算TSP问题
% 符号与变量说明
%
% a：温度衰减函数的参数
%
% t0：初始温度
%
% tf：最终温度
%
% t：当前温度
%
% coordinates：城市的x，y坐标
%
% Markov_length：Markov链长度
%
% dist_matrix：每个城市到另一个城市的距离矩阵（对称矩阵）
%
% tmp：临时变量
%
% E：路程总距离
%
% sol：路程选择
%
% sol_new是每次产生的新解；sol_current是当前解；sol_best是冷却中的最好解；
%
% p 为绘制动画演示而保留的参数
%
% total_sol 为绘制动画演示而保留的参数
%
% E_total 为绘制收敛图而保留的参数
clear,clc,close all
%%
% 载入数据，计算距离矩阵，绘图
coordinates = xlsread('城市坐标.xlsx',1,'B2:C53');
amount = size(coordinates,1); % 城市的数目
coor_x_tmp1 = coordinates(:,1) * ones(1,amount);
coor_x_tmp2 = coor_x_tmp1';
coor_y_tmp1 = coordinates(:,2) * ones(1,amount);
coor_y_tmp2 = coor_y_tmp1';
dist_matrix = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + (coor_y_tmp1-coor_y_tmp2).^2);
%%
% 初始化算法
a = 0.99;t0 = 97; tf = 3; t = t0;
Markov_length = 10000;
sol_new = 1:amount; sol_current = sol_new; sol_best = sol_new; %设置初始解
E_new = total_dist(dist_matrix,sol_new);E_best = E_new;E_current=E_new; %计算初始解的长度
plot(coordinates(:,1),coordinates(:,2),'ro');
title('路线图')
xlabel('X坐标')
ylabel('Y坐标')
hold on
title('城市分布')
total_sol=zeros(10000,amount);
total_sol(1,:)=sol_current;
%%
% 开始迭代
tic;n=1;j=1;
while t>=tf
for r=1:Markov_length % Markov链长度
% 产生随机扰动
sol_new=get_sol(sol_current,amount);
% 计算新解的长度，判断是否保留
E_new = total_dist(dist_matrix,sol_new);
if E_new < E_current
sol_current = sol_new;
E_current = E_new;
p=1;
elseif E_new == E_current
p=0;
else
if rand < exp(-(E_new-E_current)/t)
sol_current = sol_new;
E_current = E_new;
p=1;
else
p=0;
end
end
% 保存最优解
if E_current < E_best
E_best = E_current;
sol_best = sol_new;
end
% 保留数据更改便于后期绘画
if p
n=n+1;
total_sol(n,:)=sol_current;
end
E_total(j)=E_current;
j=j+1;
end
t=t.*a; % 控制参数t（温度）减少为原来的a倍
end
total_sol(n+1:end,:)=[];
%%
% 数据输出：
disp(['最优解为：',num2str(sol_best)])
disp(['最短距离：',num2str(E_best)])
disp(['程序耗时',num2str(toc),'s'])
%%
% 数据可视化：
figure
subplot(1,2,1)
plot(E_total)
title('收敛图')
xlabel('迭代次数')
ylabel('当前路线长度')
ylim([7000,1.6e4])
subplot(1,2,2)
x = coordinates(:,1);y = coordinates(:,2);
x = x(sol_best);y = y(sol_best);
plot(x,y,'b-')
hold on
plot([x(end),x(1)],[y(end),y(1)],'b-')
hold off
title('路线图')
xlabel('X坐标')
ylabel('Y坐标')
%%
% 动画
figure
title('路线图')
xlabel('X坐标')
ylabel('Y坐标')
h=plot(coordinates(:,1),coordinates(:,2),'r*');
hold on
ind=total_sol(1,:);
x = coordinates(:,1);y = coordinates(:,2);
x = x(ind);y = y(ind);
h1=plot(x,y,'b-');
h2=plot([x(end),x(1)],[y(end),y(1)],'b-');
pause
for i = 2:n
delete(h1),delete(h2);
ind=total_sol(i,:);
x = coordinates(:,1);y = coordinates(:,2);
x = x(ind);y = y(ind);
h1=plot(x,y,'b-');
h2=plot([x(end),x(1)],[y(end),y(1)],'b-');
drawnow
end
%% 函数部分
%%
% 距离计算函数：
function f=total_dist(dist_matrix,sol)
amount=length(sol);
f = 0;
for i = 1 : (amount-1)
f= f + dist_matrix(sol(i),sol(i+1));
end
f = f +dist_matrix(sol(amount),sol(1));
end
%%
% 获取新解函数：
function sol_new=get_sol(sol_current,amount)
if (rand < 0.5)
% 二交换 ind1与ind2交换顺序
ind1 = 0; ind2 = 0;
while (ind1 == ind2)
ind1 = ceil(rand*amount);
ind2 = ceil(rand*amount);
end
sol_new=sol_current;
sol_new(ind1)=sol_current(ind2);
sol_new(ind2)=sol_current(ind1);
else
% 三交换 把 ind1 和 ind2 中的序列置于 ind3 后面
ind1 = 0; ind2 = 0; ind3 = 0;
while (ind1 == ind2) || (ind1 == ind3) || (ind2 == ind3) || (abs(ind1-ind2) == 1)
ind1 = ceil(rand.*amount);
ind2 = ceil(rand.*amount);
ind3 = ceil(rand.*amount);
end
% 确保ind1 < ind2 < ind3
tmp = sort([ind1,ind2,ind3]);
ind1 = tmp(1);
ind2 = tmp(2);
ind3 = tmp(3);
sol_new = sol_current;
sol_new((ind1+1):(ind1+ind3-ind2+1)) = sol_current((ind2):(ind3));
sol_new((ind1+ind3-ind2+2):ind3) = sol_current((ind1+1):(ind2-1));
end
end
数据文件 城市坐标.xlsx 请关注微信公众号 我不是wc 回复 016 获取
2021年9月6日20:26:15
|