整理一下數學建模會用到的算法，供比賽時候參考食用。
——————————————————————————————————————————
旅行商問題（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