模擬退火算法-旅行商問題-matlab實現
阿新 • • 發佈：2019-01-20
是否 -i ren isp close 交換 coo 準則 matrix
整理一下數學建模會用到的算法，供比賽時候參考食用。
——————————————————————————————————————————
旅行商問題（TSP）：
給定一系列城市和每對城市之間的距離，求解訪問每一座城市一次並回到起始城市的最短回路。
它是組合優化中的一個NP困難問題，在運籌學和理論計算機科學中非常重要。
以下兩個程序，在不同數據集合情況下表現有所差別，理論上第一個程序的解更為優化。
1 clear 2 clc 3 a = 0.99; %溫度衰減函數的參數 4 t0 = 97; %初始溫度 5 tf = 3; %終止溫度6 t = t0; 7 Markov_length = 10000; %Markov鏈長度 8 9 % load data.txt 10 % x = data(:, 1:2:8); x = x(:); 11 % y = data(:, 2:2:8); y = y(:); 12 % data = [70,40;x, y]; 13 % coordinates = data; 14 coordinates = [ 15 1 565.0 575.0; 2 25.0 185.0; 3 345.0 750.0; 16 4 945.0 685.0; 5845.0 655.0; 6 880.0 660.0; 17 7 25.0 230.0; 8 525.0 1000.0; 9 580.0 1175.0; 18 10 650.0 1130.0; 11 1605.0 620.0; 12 1220.0 580.0; 19 13 1465.0 200.0; 14 1530.0 5.0; 15 845.0 680.0; 20 16 725.0 370.0; 17 145.0 665.0; 18 415.0 635.0; 21 19 510.0 875.0; 20 560.0 365.0; 21 300.0 465.0; 22 22 520.0 585.0; 23 480.0 415.0; 24 835.0 625.0; 23 25 975.0 580.0; 26 1215.0 245.0; 27 1320.0 315.0; 24 28 1250.0 400.0; 29 660.0 180.0; 30 410.0 250.0; 25 31 420.0 555.0; 32 575.0 665.0; 33 1150.0 1160.0; 26 34 700.0 580.0; 35 685.0 595.0; 36 685.0 610.0; 27 37 770.0 610.0; 38 795.0 645.0; 39 720.0 635.0; 28 40 760.0 650.0; 41 475.0 960.0; 42 95.0 260.0; 29 43 875.0 920.0; 44 700.0 500.0; 45 555.0 815.0; 30 46 830.0 485.0; 47 1170.0 65.0; 48 830.0 610.0; 31 49 605.0 625.0; 50 595.0 360.0; 51 1340.0 725.0; 32 52 1740.0 245.0; 33 ]; 34 coordinates(:,1) = []; 35 amount = size(coordinates,1); %城市的數目 36 %通過向量化的方法計算距離矩陣 37 dist_matrix = zeros(amount,amount); 38 coor_x_tmp1 = coordinates(:,1) * ones(1,amount); 39 coor_x_tmp2 = coor_x_tmp1‘; 40 coor_y_tmp1 = coordinates(:,2) * ones(1,amount); 41 coor_y_tmp2 = coor_y_tmp1‘; 42 dist_matrix = sqrt((coor_x_tmp1 - coor_x_tmp2).^2 + (coor_y_tmp1 - coor_y_tmp2).^2); 43 44 45 sol_new = 1:amount; %產生初始解，sol_new是每次產生的新解 46 sol_current = sol_new; %sol_current是當前解 47 sol_best = sol_new; %sol_best是冷卻中的最好解 48 E_current = inf; %E_current是當前解對應的回路距離 49 E_best = inf; %E_best是最優解 50 p = 1; 51 52 53 rand(‘state‘, sum(clock)); 54 55 for j = 1:10000 56 sol_current = [randperm(amount)]; 57 E_current = 0; 58 for i=1:(amount-1) 59 E_current = E_current+dist_matrix(sol_current(i), sol_current(i+1)); 60 end 61 if E_current<E_best 62 sol_best = sol_current; 63 E_best = E_current; 64 end 65 end 66 67 68 69 70 while t >= tf 71 for r = 1:Markov_length %Markov鏈長度 72 %產生隨機擾動 73 if(rand < 0.5) 74 %兩交換 75 ind1 = 0; 76 ind2 = 0; 77 while(ind1 == ind2) 78 ind1 = ceil(rand * amount); 79 ind2 = ceil(rand * amount); 80 end 81 tmp1 = sol_new(ind1); 82 sol_new(ind1) = sol_new(ind2); 83 sol_new(ind2) = tmp1; 84 else 85 %三交換 86 ind1 = 0; 87 ind2 = 0; 88 ind3 = 0; 89 while( (ind1 == ind2) || (ind1 == ind3) || (ind2 == ind3) || (abs(ind1 -ind2) == 1) ) 90 ind1 = ceil(rand * amount); 91 ind2 = ceil(rand * amount); 92 ind3 = ceil(rand * amount); 93 end 94 tmp1 = ind1; 95 tmp2 = ind2; 96 tmp3 = ind3; 97 %確保 ind1 < ind2 < ind3 98 if(ind1 < ind2) && (ind2 < ind3); 99 elseif(ind1 < ind3) && (ind3 < ind2) 100 ind1 = tmp1; ind2 = tmp3; ind3 = tmp2; 101 elseif(ind2 < ind1) && (ind1 < ind3) 102 ind1 = tmp2; ind2 = tmp1; ind3 = tmp3; 103 elseif(ind2 < ind3) && (ind3 < ind1) 104 ind1 = tmp2; ind2 = tmp3; ind3 = tmp1; 105 elseif(ind3 < ind1) && (ind1 < ind2) 106 ind1 = tmp3; ind2 = tmp1; ind3 = tmp2; 107 elseif(ind3 < ind2) && (ind2 < ind1) 108 ind1 = tmp3; ind2 = tmp2; ind3 = tmp1; 109 end 110 111 tmplist1 = sol_new((ind1 + 1):(ind2 - 1)); 112 sol_new((ind1 + 1):(ind1 + (ind3 - ind2 + 1) )) = sol_new((ind2):(ind3)); 113 sol_new((ind1 + (ind3 - ind2 + 1) + 1):(ind3)) = tmplist1; 114 end 115 116 %檢查是否滿足約束 117 118 %計算目標函數值（即內能） 119 E_new = 0; 120 for i = 1:(amount - 1) 121 E_new = E_new + dist_matrix(sol_new(i),sol_new(i + 1)); 122 end 123 %再算上從最後一個城市到第一個城市的距離 124 E_new = E_new + dist_matrix(sol_new(amount),sol_new(1)); 125 126 if E_new < E_current 127 E_current = E_new; 128 sol_current = sol_new; 129 if E_new < E_best 130 E_best = E_new; 131 sol_best = sol_new; 132 end 133 else 134 %若新解的目標函數值大於當前解， 135 %則僅以一定概率接受新解 136 if rand < exp(-(E_new - E_current) / t) 137 E_current = E_new; 138 sol_current = sol_new; 139 else 140 sol_new = sol_current; 141 end 142 143 end 144 end 145 146 t = t * a; %控制參數t（溫度）減少為原來的a倍 147 end 148 149 E_best = E_best+dist_matrix(sol_best(end), sol_best(1)); 150 151 disp(‘最優解為:‘); 152 disp(sol_best); 153 disp(‘最短距離:‘); 154 disp(E_best); 155 156 data1 = zeros(length(sol_best),2 ); 157 for i = 1:length(sol_best) 158 data1(i, :) = coordinates(sol_best(1,i), :); 159 end 160 161 data1 = [data1; coordinates(sol_best(1,1),:)]; 162 163 164 figure 165 plot(coordinates(:,1)‘, coordinates(:,2)‘, ‘*k‘, data1(:,1)‘, data1(:, 2)‘, ‘r‘); 166 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
模擬退火算法-旅行商問題-matlab實現