问题:柏林有52座城市,每座城市的坐标数据在coordinates中,求TSP的解答。
a = 0.99;
T0 = 97; Tf = 3; t = T0;
MarkovLength = 10000;
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);
%使用矩阵方法计算巧妙得到距离矩阵
distMatrix = zeros(amount, amount);
coorXTemp1 = coordinates(:,1) * ones(1,amount);
coorXTemp2 = coorXTemp1';
coorYTemp1 = coordinates(:,2) * ones(1,amount);
coorYTemp2 = coorYTemp1';
distMatrix = sqrt((coorXTemp1-coorXTemp2).^2 + (coorYTemp1-coorYTemp2).^2);
newSol =1:amount;
curResult = inf; bestResult = inf;
curSol = newSol; bestSol = newSol;
p=1;
while t >= Tf
for r = 1:MarkovLength
%产生随机扰动
if rand<0.5
%两交换
index1 = 0; index2 = 0;
while index1 == index2
index1 = ceil(rand*amount);
index2 = ceil(rand*amount);
end
%三段式交换
tmp = newSol(index1);
newSol(index1) = newSol(index2);
newSol(index2) = tmp;
else
%三交换
index1 = 0; index2 = 0; index3 = 0;
while (index1 == index2)||(index1 == index3)||(index2 == index3)||(abs(index1-index2) == 1)
index1 = ceil(rand*amount);
index2 = ceil(rand*amount);
index3 = ceil(rand*amount);
end
tmp1 = index1; tmp2 = index2; tmp3 = index3;
%确保index1 < index2 < index3
if(index1 < index2)&&(index2<index3);
elseif (index1 < index3)&&(index3<index2)
index2 = tmp3; index3 = tmp2;
elseif (index2 < index1)&&(index1<index3)
index1 = tmp2; index2 = tmp1;
elseif (index2 < index3)&&(index3<index1)
index1 = tmp2; index2 = tmp3; index3 = tmp1;
elseif (index3 < index1)&&(index1<index2)
index1 = tmp3; index2 = tmp1; index3 = tmp2;
elseif (index3 < index2)&&(index2<index1);
index1 = tmp3; index3 = tmp1;
end
tmpRoute = newSol((index1+1):(index2-1));
newSol((index1+1):(index1 + index3-index2 + 1)) = newSol(index2:index3);
newSol((index1 + index3 - index2 + 2):index3) = tmpRoute;
end
newResult = 0;
for i = 1:(amount-1)
newResult = newResult + distMatrix(newSol(i),newSol(i+1));
end
newResult = newResult + distMatrix(newSol(amount), newSol(1));
if newResult < curResult
curResult = newResult;
curSol = newSol;
if newResult < bestResult
bestResult = newResult;
bestSol = newSol;
end
else
if rand < exp(-(newResult - curResult)/t)
curResult = newResult;
curSol = newSol;
else
newSol = curSol;
end
end
end
t = t * a;
end
disp('最优解为:');
disp(bestSol);
disp('最短距离为:');
disp(bestResult);
运行结果
最优解为:
1 至 23 列19 41 8 9 10 43 33 51 11 52 14 13 47 26 27 28 12 25 4 6 15 5 24
24 至 46 列
48 38 37 40 39 36 35 34 44 46 16 29 50 20 23 30 2 7 42 21 17 3 18
47 至 52 列
31 22 1 49 32 45
最短距离为:
7.5444e+03