整理一下数学建模会用到的算法，供比赛时候参考食用。
——————————————————————————————————————————
旅行商问题（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