上图，我的手机壁纸，嘻嘻嘻！！！
本文约6000字，预计阅读15分钟
今天的题目如下：
这是一道标准的tsp问题，我们可以用很多种方法去解决。
目标规划问题：
用lingo软件去写规划问题：
定义两个矩阵，d代表的是距离，s代表的是边是否存在
这个矩阵的输入有点麻烦，我们只需要保距离矩阵的下三角局矩阵就行。
if s(i,j)=1;代表从i走到j
if s(i,j)=0;代表这条边不存在
!TSP quesion;
MODEL:
SETS:
city/1..30/;
link(city,city)|&1#GT#&2:d,s;
ENDSETS
DATA:
d = 10.770
29.967 24.042
35.777 25.060 29.428
45.343 36.056 47.096 18.111
39.319 38.079 61.057 43.566 35.355
45.000 40.460 16.643 43.186 61.294 77.698
58.310 52.498 28.601 49.396 67.052 88.238 14.318
34.540 27.803 5.000 29.000 47.043 63.820 14.560 24.759
48.877 48.384 29.069 58.421 76.164 86.377 18.601 27.731 29.833
41.049 36.125 12.207 39.051 57.140 73.246 4.472 17.464 10.198 21.024
46.141 35.511 38.275 10.630 14.866 47.760 50.160 53.935 36.878 66.708 46.390
38.949 28.302 32.757 3.606 15.524 43.829 46.043 51.546 32.062 61.660 42.000 7.211
63.781 59.666 35.805 60.166 78.102 96.799 19.209 12.166 33.121 23.000 23.601 65.490 62.586
75.073 70.937 47.011 70.228 87.932 107.898 30.480 20.881 44.102 32.016 34.828 74.733 72.422 11.314
58.241 47.539 41.012 24.000 31.623 65.192 47.424 46.390 37.643 65.765 44.777 17.464 22.204 58.549 66.000
54.708 43.966 39.051 20.025 27.803 61.098 46.819 47.043 36.056 64.885 43.863 13.416 18.111 59.135 67.119 4.123
30.232 25.807 4.472 33.734 51.245 63.530 14.866 28.178 8.062 25.000 10.817 42.720 37.108 33.971 45.277 45.277 43.417
37.802 36.401 17.464 46.872 64.382 74.465 13.342 27.000 19.235 12.166 13.038 55.660 50.220 27.731 38.588 56.613 55.227 13.153
36.674 37.483 22.825 51.546 68.447 75.000 20.881 34.132 25.612 12.728 20.591 60.926 55.027 33.242 43.463 63.253 61.612 18.358 7.616
49.396 50.636 34.205 63.561 80.895 88.057 26.173 35.777 35.847 8.062 28.018 72.422 66.940 30.265 38.210 72.719 71.589 29.833 16.763 13.153
58.694 47.927 45.000 23.087 26.401 61.131 53.141 53.151 42.190 71.176 50.160 14.000 20.396 65.276 73.027 7.280 6.325 49.406 61.400 67.676 77.827
60.828 50.120 49.092 25.060 24.739 60.017 57.871 58.138 46.530 75.769 54.781 14.866 21.932 70.257 78.026 12.166 11.180 53.535 65.765 71.868 82.292 5.000
96.177 89.185 66.212 79.209 94.202 121.918 52.887 38.601 61.717 62.008 55.973 79.404 80.056 39.013 32.280 64.885 67.742 66.483 64.938 71.449 69.181 72.007 76.485
65.460 57.697 35.903 47.634 63.632 89.939 26.683 15.000 31.048 42.544 28.071 49.193 48.826 25.239 29.614 37.483 39.294 37.216 40.025 47.539 50.606 44.721 49.649 32.016
61.400 53.310 32.249 42.638 58.669 85.041 25.080 15.811 27.295 42.202 25.710 44.283 43.829 27.313 33.136 33.136 34.713 34.000 38.275 45.880 50.220 40.311 45.277 36.878 5.000
73.110 63.506 46.872 45.618 57.385 89.067 43.566 34.713 41.976 61.221 43.382 42.638 45.277 45.486 49.041 26.249 29.698 49.729 56.356 63.953 69.203 33.015 37.216 39.560 20.248 19.105
68.000 58.138 43.012 39.395 50.990 82.765 41.869 34.986 38.275 60.108 41.049 36.235 38.949 46.519 51.420 20.000 23.345 46.239 54.083 61.587 67.941 26.926 31.305 45.188 21.840 19.235 6.403
59.076 49.497 33.526 33.015 47.011 76.551 33.242 28.460 28.792 51.740 32.016 32.202 33.302 40.522 47.096 19.235 21.190 36.770 45.000 52.431 59.414 26.476 31.401 47.202 18.248 14.000 14.036 9.487
57.489 47.381 52.811 24.187 14.318 49.041 64.498 67.268 51.420 81.253 60.828 14.560 20.591 79.101 87.824 24.187 21.541 57.245 70.214 75.392 86.977 17.205 13.454 89.067 60.729 56.045 50.220 44.102 42.720;
ENDDATA
MIN=@SUM(link:d*s);
@SUM(city(j)|j#GT#1:S(j,1))=2;
!与第1个城市相连的有两个城市;
!与第i个城市相连有两个城市;
@FOR(city(i)|i#GT#1:
@SUM(city(j)|j#GT#i:s(j,i))+@SUM(city(k)|k#LT#i:s(i,k))=2);
@FOR(link:@BIN(s));
end
模拟退火算法
我们在过程当中，需要保存三种状态，一个是当前的状态，第二个是收到扰动生成的状态，第三个是最佳的状态。
我们给定测试的次数尽量大，用python跑一趟竟然用了50秒钟，而matlab的算法不超过5秒就能得到非常精确的解法
clc;
clear;
%% 得到距离矩阵
num=[41 94;37 84;54 67;25 62;
7 64;2 99;68 58;71 44;
54 62;83 69;64 60;18 54;
22 60;83 46;91 38;25 38;
24 42;58 69;71 71;74 78;
87 76;18 40;13 40;82 7;
62 32;58 35;45 21;41 26];
num(29,1)=44;
num(29,2)=35;
num(30,1)=4;
num(30,2)=50;
dis=zeros(30,30);
for i=1:29
for j=i+1:30
dis(i,j)=sqrt((num(i,1)-num(j,1))^2+(num(i,2)-num(j,2))^2);
dis(j,i)=dis(i,j);
end
end
alpha=0.99;
t=(1:100);
markovlen=10000;
solutionnew=(1:30);
index=30;
solutioncurrent=(1:30);
solutionbest=(1:30);
valuecurrent=99999999;
valuebest=valuecurrent;
res=[];
temp=t(100);
while temp >t(1)
for i=1:markovlen
% 扰动方式，二交换
if rand >0.5
while true
pos1=ceil(rand*30);
pos2=ceil(rand*30);
if pos1 ~= pos2
break;
end
end
solutiontemp=solutionnew(pos1);
solutionnew(pos1)=solutionnew(pos2);
solutionnew(pos2)=solutiontemp;
% 另一种扰动方式，三交换
else
while true
pos1=ceil(rand*30);
pos2=ceil(rand*30);
pos3=ceil(rand*30);
if (pos1~=pos2 && pos2~=pos3 && pos3~=pos1)
break;
end
end
if pos1>pos2
linshi=pos1;
pos1=pos2;
pos2=linshi;
end
if pos2 > pos3
linshi=pos3;
pos3=pos2;
pos2=linshi;
end
if pos1>pos2
linshi=pos1;
pos1=pos2;
pos2=linshi;
end
tempindex=solutionnew(pos1:pos2-1);
solutionnew(pos1:pos3-pos2+pos1)=solutionnew(pos2:pos3);
solutionnew(pos3-pos2+pos1+1:pos3)=tempindex;
% solutiontemp=solutionnew(pos1);
% solutionnew(pos1)=solutionnew(pos2);
% solutionnew(pos2)=solutionnew(pos3);
% solutionnew(pos3)=solutiontemp;
valuenew = 0;
for k=1:index-1
valuenew=valuenew+dis(solutionnew(k),solutionnew(k+1));
end
valuenew=valuenew+dis(solutionnew(1),solutionnew(30));
if valuenew < valuecurrent
valuecurrent = valuenew;
solutioncurrent = solutionnew;
if valuenew < valuebest
valuebest = valuenew;
solutionbest = solutionnew;
end
% 按照一定的概率去接受解
else
if rand < exp(-(valuenew-valuecurrent)/temp)
valuecurrent = valuenew;
solutioncurrent = solutionnew;
else
solutionnew=solutioncurrent;
end
end
end
end
temp=temp*alpha;
res=[res,valuebest];
temp
end
length=size(res);
x=(1:length(2));
plot(x,res);
xlabel('迭代次数');
ylabel('当前的最佳结果');
axis([0 1000 0 2000])
这是我暂时写出来的两种解法，当然还有状压DP的方式也可以去做，这块还没弄得太懂。
希望在大学毕业的时候，我们不会后悔过去四年在大学里虚度了自己的光阴。
无论是星星的闪烁还是快乐的生活都要主动伸出双手去掌握。
有疑问请后台留言：shirandexiaowo@foxmail.com
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
