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