Matlab學習手記——旅行商TSP問題：模拟退火法

• 原理

• 結果 

• 坐标檔案

5.6500000e+02   5.7500000e+02
   2.5000000e+01   1.8500000e+02
   3.4500000e+02   7.5000000e+02
   9.4500000e+02   6.8500000e+02
   8.4500000e+02   6.5500000e+02
   8.8000000e+02   6.6000000e+02
   2.5000000e+01   2.3000000e+02
   5.2500000e+02   1.0000000e+03
   5.8000000e+02   1.1750000e+03
   6.5000000e+02   1.1300000e+03
   1.6050000e+03   6.2000000e+02
   1.2200000e+03   5.8000000e+02
   1.4650000e+03   2.0000000e+02
   1.5300000e+03   5.0000000e+00
   8.4500000e+02   6.8000000e+02
   7.2500000e+02   3.7000000e+02
   1.4500000e+02   6.6500000e+02
   4.1500000e+02   6.3500000e+02
   5.1000000e+02   8.7500000e+02
   5.6000000e+02   3.6500000e+02
   3.0000000e+02   4.6500000e+02
   5.2000000e+02   5.8500000e+02
   4.8000000e+02   4.1500000e+02
   8.3500000e+02   6.2500000e+02
   9.7500000e+02   5.8000000e+02
   1.2150000e+03   2.4500000e+02
   1.3200000e+03   3.1500000e+02
   1.2500000e+03   4.0000000e+02
   6.6000000e+02   1.8000000e+02
   4.1000000e+02   2.5000000e+02
   4.2000000e+02   5.5500000e+02
   5.7500000e+02   6.6500000e+02
   1.1500000e+03   1.1600000e+03
   7.0000000e+02   5.8000000e+02
   6.8500000e+02   5.9500000e+02
   6.8500000e+02   6.1000000e+02
   7.7000000e+02   6.1000000e+02
   7.9500000e+02   6.4500000e+02
   7.2000000e+02   6.3500000e+02
   7.6000000e+02   6.5000000e+02
   4.7500000e+02   9.6000000e+02
   9.5000000e+01   2.6000000e+02
   8.7500000e+02   9.2000000e+02
   7.0000000e+02   5.0000000e+02
   5.5500000e+02   8.1500000e+02
   8.3000000e+02   4.8500000e+02
   1.1700000e+03   6.5000000e+01
   8.3000000e+02   6.1000000e+02
   6.0500000e+02   6.2500000e+02
   5.9500000e+02   3.6000000e+02
   1.3400000e+03   7.2500000e+02
   1.7400000e+03   2.4500000e+02
           

• 源碼

function [sol_best, E_best] = SAA_TSP(coordinates, Markov_length, ratio, t_start, t_end)
%{
函數功能：模拟退火求解旅行商問題；
輸入：
  coordinates：城市坐标，兩列；
  Markov_length：Markov 鍊長度；
  ratio：溫度衰減系數；
  t_start：初始溫度；
  t_end：結束溫度；
輸出：
  sol_best：最優路線；
  E_best：最短距離；
示例：
clear; clc;
ratio = 0.99;    
t_start = 100; 
t_end = 1;        
Markov_length = 5000;
coordinates = load('CityCoordinates.txt');
[sol_best, E_best] = SAA_TSP(coordinates, Markov_length, ratio, t_start, t_end);
disp(['最優解為：', num2str(sol_best)]);
disp(['最短距離：', num2str(E_best)]);
x_path = [coordinates(sol_best, 1); coordinates(sol_best(1), 1)];
y_path = [coordinates(sol_best, 2); coordinates(sol_best(1), 2)];
plot(x_path, y_path, 'LineWidth', 2)
%} 
% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 
if nargin < 5
   t_end = 0.01 * t_start;
end
if nargin < 4
   t_start = 100; 
end
if nargin < 3
   ratio = 0.99;
end
if nargin < 2
   Markov_length = 5000;
end
if nargin < 1
   error('輸入參數不足！');
end
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);
sol_new = 1 : amount;         % 産生初始解
% sol_new是每次産生的新解；sol_current是目前解；sol_best是冷卻中的最好解；
% E_current是目前解對應的回路距離；E_new是新解的回路距離；E_best是最優解的
E_current = inf;        
E_best = inf;
sol_current = sol_new; 
sol_best = sol_new;
t = t_start;          
while t >= t_end
    % 目前溫度下循環
    for r = 1 : Markov_length
        % 産生随機擾動
        if (rand < 0.5) % 随機決定是進行兩交換還是三交換
            % 兩交換
            ind = randperm(amount, 2);
            ind1 = ind(1);
            ind2 = ind(2);            
            tmp1 = sol_new(ind1);
            sol_new(ind1) = sol_new(ind2);
            sol_new(ind2) = tmp1;
        else
            % 三交換
            ind = randperm(amount, 3);
            % 確定ind1 < ind2 < ind3
            ind = sort(ind);
            ind1 = ind(1);
            ind2 = ind(2);
            ind3 = ind(3);
            tmplist1 = sol_new(ind1 + 1 : ind2 - 1);
            sol_new(ind1 + 1 : ind1 + ind3 - ind2 + 1) = sol_new(ind2 : ind3);
            sol_new(ind1 + ind3 - ind2 + 2 : ind3) = tmplist1;
        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 * ratio; % 降溫
end