模拟退火算法-旅行商问题-matlab实现
整理一下数学建模会用到的算法，供比赛时候参考食用。
——————————————————————————————————————————
旅行商问题（TSP）：
给定一系列城市和每对城市之间的距离，求解访问每一座城市一次并回到起始城市的最短回路。
它是组合优化中的一个NP困难问题，在运筹学和理论计算机科学中非常重要。
以下两个程序，在不同数据集合情况下表现有所差别，理论上第一个程序的解更为优化。
clear
clc
a = 0.99; %温度衰减函数的参数
t0 = 97; %初始温度
tf = 3; %终止温度
t = t0;
Markov_length = 10000; %Markov链长度
% load data.txt
% x = data(:, 1:2:8); x = x(:);
% y = data(:, 2:2:8); y = y(:);
% data = [70,40;x, y];
% coordinates = data;
coordinates = [
1 565.0 575.0; 2 25.0 185.0; 3 345.0 750.0;
4 945.0 685.0; 5 845.0 655.0; 6 880.0 660.0;
7 25.0 230.0; 8 525.0 1000.0; 9 580.0 1175.0;
10 650.0 1130.0; 11 1605.0 620.0; 12 1220.0 580.0;
13 1465.0 200.0; 14 1530.0 5.0; 15 845.0 680.0;
16 725.0 370.0; 17 145.0 665.0; 18 415.0 635.0;
19 510.0 875.0; 20 560.0 365.0; 21 300.0 465.0;
22 520.0 585.0; 23 480.0 415.0; 24 835.0 625.0;
25 975.0 580.0; 26 1215.0 245.0; 27 1320.0 315.0;
28 1250.0 400.0; 29 660.0 180.0; 30 410.0 250.0;
31 420.0 555.0; 32 575.0 665.0; 33 1150.0 1160.0;
34 700.0 580.0; 35 685.0 595.0; 36 685.0 610.0;
37 770.0 610.0; 38 795.0 645.0; 39 720.0 635.0;
40 760.0 650.0; 41 475.0 960.0; 42 95.0 260.0;
43 875.0 920.0; 44 700.0 500.0; 45 555.0 815.0;
46 830.0 485.0; 47 1170.0 65.0; 48 830.0 610.0;
49 605.0 625.0; 50 595.0 360.0; 51 1340.0 725.0;
52 1740.0 245.0;
];
coordinates(:,1) = [];
amount = size(coordinates,1); %城市的数目
%通过向量化的方法计算距离矩阵
dist_matrix = zeros(amount,amount);
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);
sol_new = 1:amount; %产生初始解，sol_new是每次产生的新解
sol_current = sol_new; %sol_current是当前解
sol_best = sol_new; %sol_best是冷却中的最好解
E_current = inf; %E_current是当前解对应的回路距离
E_best = inf; %E_best是最优解
p = 1;
rand('state', sum(clock));
for j = 1:10000
sol_current = [randperm(amount)];
E_current = 0;
for i=1:(amount-1)
E_current = E_current+dist_matrix(sol_current(i), sol_current(i+1));
end
if E_current<E_best
sol_best = sol_current;
E_best = E_current;
end
end
while t >= tf
for r = 1:Markov_length %Markov链长度
%产生随机扰动
if(rand < 0.5)
%两交换
ind1 = 0;
ind2 = 0;
while(ind1 == ind2)
ind1 = ceil(rand * amount);
ind2 = ceil(rand * amount);
end
tmp1 = sol_new(ind1);
sol_new(ind1) = sol_new(ind2);
sol_new(ind2) = tmp1;
else
%三交换
ind=ceil(amount*rand(1,3));
ind = sort(ind);
sol_new = sol_new(1, [1:ind(1)-1, ind(2)+1:ind(3),ind(1):ind(2),ind(3)+1:end]);
end
%检查是否满足约束
%计算目标函数值（即内能）
E_new = 0;
for i = 1:(amount - 1)
E_new = E_new + dist_matrix(sol_new(i),sol_new(i + 1));
end
%再算上从最后一个城市到第一个城市的距离
E_new = E_new + dist_matrix(sol_new(amount),sol_new(1));
if E_new < E_current
E_current = E_new;
sol_current = sol_new;
if E_new < E_best
E_best = E_new;
sol_best = sol_new;
end
else
%若新解的目标函数值大于当前解，
%则仅以一定概率接受新解
if rand < exp(-(E_new - E_current) / t)
E_current = E_new;
sol_current = sol_new;
else
sol_new = sol_current;
end
end
end
t = t * a; %控制参数t（温度）减少为原来的a倍
end
E_best = E_best+dist_matrix(sol_best(end), sol_best(1));
disp('最优解为:');
disp(sol_best);
disp('最短距离:');
disp(E_best);
data1 = zeros(length(sol_best),2 );
for i = 1:length(sol_best)
data1(i, :) = coordinates(sol_best(1,i), :);
end
data1 = [data1; coordinates(sol_best(1,1),:)];
figure
plot(coordinates(:,1)', coordinates(:,2)', '*k', data1(:,1)', data1(:, 2)', 'r');
title( [ '近似最短路径如下，路程为' , num2str( E_best ) ] ) ;另一种根据《数学建模算法与应用—司守奎》中算法改编：
clc;
clear;
close all;
coordinates = [
2 25.0 185.0; 3 345.0 750.0;
4 945.0 685.0; 5 845.0 655.0; 6 880.0 660.0;
7 25.0 230.0; 8 525.0 1000.0; 9 580.0 1175.0;
10 650.0 1130.0; 11 1605.0 620.0; 12 1220.0 580.0;
13 1465.0 200.0; 14 1530.0 5.0; 15 845.0 680.0;
16 725.0 370.0; 17 145.0 665.0; 18 415.0 635.0;
19 510.0 875.0; 20 560.0 365.0; 21 300.0 465.0;
22 520.0 585.0; 23 480.0 415.0; 24 835.0 625.0;
25 975.0 580.0; 26 1215.0 245.0; 27 1320.0 315.0;
28 1250.0 400.0; 29 660.0 180.0; 30 410.0 250.0;
31 420.0 555.0; 32 575.0 665.0; 33 1150.0 1160.0;
34 700.0 580.0; 35 685.0 595.0; 36 685.0 610.0;
37 770.0 610.0; 38 795.0 645.0; 39 720.0 635.0;
40 760.0 650.0; 41 475.0 960.0; 42 95.0 260.0;
43 875.0 920.0; 44 700.0 500.0; 45 555.0 815.0;
46 830.0 485.0; 47 1170.0 65.0; 48 830.0 610.0;
49 605.0 625.0; 50 595.0 360.0; 51 1340.0 725.0;
52 1740.0 245.0;
];
coordinates(:,1) = [];
data = coordinates;
% 读取数据
% load data.txt;
% x = data(:, 1:2:8); x = x(:);
% y = data(:, 2:2:8); y = y(:);
x = data(:, 1);
y = data(:, 2);
start = [565.0 575.0];
data = [start; data;start];
% data = [start; x, y;start];
% data = data*pi/180;
% 计算距离的邻接表
count = length(data(:, 1));
d = zeros(count);
for i = 1:count-1
for j = i+1:count
% temp = cos(data(i, 1)-data(j,1))*cos(data(i,2))*cos(data(j,2))...
% +sin(data(i,2))*sin(data(j,2));
d(i,j)=( sum( ( data( i , : ) - data( j , : ) ) .^ 2 ) ) ^ 0.5 ;
% d(i,j) = 6370*acos(temp);
end
end
d =d + d'; % 对称 i到j==j到i
S0=[]; % 存储初值
Sum=inf; % 存储总距离
rand('state', sum(clock));
% 求一个较为优化的解，作为初值
for j = 1:10000
S = [1 1+randperm(count-2), count];
temp = 0;
for i=1:count-1
temp = temp+d(S(i), S(i+1));
end
if temp<Sum
S0 = S;
Sum = temp;
end
end
e = 0.1^40; % 终止温度
L = 2000000; % 最大迭代次数
at = 0.999999; % 降温系数
T = 2; % 初温
% 退火过程
for k = 1:L
% 产生新解
c =1+floor((count-1)*rand(1,2));
c = sort(c);
c1 = c(1); c2 = c(2);
if c1==1
c1 = c1+1;
end
if c2==1
c2 = c2+1;
end
% 计算代价函数值
df = d(S0(c1-1), S0(c2))+d(S0(c1), S0(c2+1))-...
(d(S0(c1-1), S0(c1))+d(S0(c2), S0(c2+1)));
% 接受准则
if df<0
S0 = [S0(1: c1-1), S0(c2:-1:c1), S0(c2+1:count)];
Sum = Sum+df;
elseif exp(-df/T) > rand(1)
S0 = [S0(1: c1-1), S0(c2:-1:c1), S0(c2+1:count)];
Sum = Sum+df;
end
T = T*at;
if T<e
break;
end
end
data1 = zeros(2, count);
% y1 = [start; x, y; start];
for i =1:count
data1(:, i) = data(S0(1,i), :)';
end
figure
plot(x, y, 'o', data1(1, :), data1(2, :), 'r');
title( [ '近似最短路径如下，路程为' , num2str( Sum ) ] ) ;
disp(Sum);
S0