博文
旅行推销商问题TSP的动态规划解法
|
2016-03-17 20:35:36
利用动态规划方法是可以精确求解旅行推销商问题(Traveling Salesman Problem)的, 虽然这种方法只适用于求解小规模的问题. 这个算法我一直没有弄清楚, 最近有个问题需要使用类似的算法来解决, 所以我就仔细研究了一下这个算法. 下面是网上几篇资料的总结.
我们先考虑一个小规模的问题, 只有4个城市, 城市间的距离由下面的矩阵 C 决定.
| ij | 0 | 1 | 2 | 3 |
|---|---|---|---|---|
| 0 | 0 | 3 | 6 | 7 |
| 1 | 5 | 0 | 2 | 3 |
| 2 | 6 | 4 | 0 | 2 |
| 3 | 3 | 7 | 5 | 0 |
值得注意的是, 这个矩阵是不对称的, 也就是说 <span class="MathJax" id="MathJax-Element-2-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">C<mrow class="MJX-TeXAtom-ORD">ij≠C<mrow class="MJX-TeXAtom-ORD">ji" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">Cij≠Cji. 如果这个矩阵是对称的, 算法还可以简化一些.
假定我们从城市0出发, 城市1, 2, 3每个经过一次, 最后回到城市0, 那么求解的递归树可以表示如下:
d(0,{1,2,3})d[0,7]=3+7=10d(1,{2,3})d[1,6]=2+5=73d(2,{1,3})d[2,5]=4+6=106d(3,{1,2})d[3,3]=5+9=147d(1,{})d[3,0]=5d(2,{})d[2,0]=6d(3,{})d[3,0]=3d(2,{3})d[2,4]=2+3=52d(3,{2})d[3,2]=5+6=113d(1,{3})d[1,4]=3+3=64d(3,{1})d[3,1]=7+5=122d(1,{2})d[1,2]=2+6=87d(2,{1})d[2,1]=4+5=9525d(3,{})d[3,0]=337d(2,{})d[2,0]=62d(1,{})d[1,0]=54Fig.1其中的 <span class="MathJax" id="MathJax-Element-3-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">d<mo stretchy="false">(i,<mo fence="false" stretchy="false">{j,k,l<mo fence="false" stretchy="false">}<mo stretchy="false">)" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">d(i,{j,k,l}) 表示由城市 <span class="MathJax" id="MathJax-Element-4-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">i" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">i 出发, 集合 <span class="MathJax" id="MathJax-Element-5-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML"><mo fence="false" stretchy="false">{j,k,l<mo fence="false" stretchy="false">}" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">{j,k,l} 中的城市 <span class="MathJax" id="MathJax-Element-6-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">j,k,l" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">j,k,l 每个经过一次需要的最小路程, 箭头上的值表示两个城市之间的距离. 很显然, <span class="MathJax" id="MathJax-Element-7-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">d<mo stretchy="false">(0,<mo fence="false" stretchy="false">{1,2,3<mo fence="false" stretchy="false">}<mo stretchy="false">)" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">d(0,{1,2,3}) 就是我们最终要求的值. 这个值可以一步一步分解下去, 最终分解为每个城市到城市0的距离, <span class="MathJax" id="MathJax-Element-8-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">d<mo stretchy="false">(i,<mo fence="false" stretchy="false">{<mo fence="false" stretchy="false">}<mo stretchy="false">)" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">d(i,{}). 将过程反过来, 向上递推, 就可以得到需要的值了.
为了方便求解并记录路径, 我们可以使用二进制数表示城市集合. 一个 <span class="MathJax" id="MathJax-Element-9-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">n" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">n 位的二进制数可以表示 <span class="MathJax" id="MathJax-Element-10-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">n" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">n 个城市的集合. 当某位为1时, 表示这个位所代表的城市出现在集合中. 我们最多有3个城市需要经历, 所以需要 <span class="MathJax" id="MathJax-Element-11-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">23=8" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">23=8 个集合. 每个集合对应的数字如下:
| 编号 | 二进制 | 集合 |
|---|---|---|
| 0 | 000 | {} |
| 1 | 001 | {1} |
| 2 | 010 | {2} |
| 3 | 011 | {1,2} |
| 4 | 100 | {3} |
| 5 | 101 | {1,3} |
| 6 | 110 | {2,3} |
| 7 | 111 | {1,2,3} |
在上面图中, <span class="MathJax" id="MathJax-Element-12-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">d<mo stretchy="false">[i,j<mo stretchy="false">]" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">d[i,j] 中的 <span class="MathJax" id="MathJax-Element-13-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">j" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">j 使用的就是表中对应的集合. 根据递推关系
其中的 <span class="MathJax" id="MathJax-Element-15-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">j<mi class="MJX-variant" mathvariant="normal">∖<mo fence="false" stretchy="false">{k<mo fence="false" stretchy="false">}" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">j∖{k} 表示从集合 <span class="MathJax" id="MathJax-Element-16-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">j" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">j 中去除集合 <span class="MathJax" id="MathJax-Element-17-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML"><mo fence="false" stretchy="false">{k<mo fence="false" stretchy="false">}" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">{k}. 集合 <span class="MathJax" id="MathJax-Element-18-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">j" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">j 去除 城市 <span class="MathJax" id="MathJax-Element-19-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">k" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">k 后所形成的集合, 其对应的编号为 <span class="MathJax" id="MathJax-Element-20-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">j−2<mrow class="MJX-TeXAtom-ORD">k−1" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">j−2k−1. 判断城市 <span class="MathJax" id="MathJax-Element-21-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">k" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">k 是否属于集合 <span class="MathJax" id="MathJax-Element-22-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">j" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">j 的方法是使用二进制与操作, 若j & 2^(k-1) == 0, <span class="MathJax" id="MathJax-Element-23-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">k" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">k 不属于集合 <span class="MathJax" id="MathJax-Element-24-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">j" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">j.
在写代码实现时, 我们可以使用一个二维表格(数组), 其大小为 <span class="MathJax" id="MathJax-Element-25-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">n2<mrow class="MJX-TeXAtom-ORD">n−1" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">n2n−1, 根据上面的关系式逐次填充计算值, 最终得到结果.
| ij | 0 {} | 1 {1} | 2 {2} | 3 {1,2} | 4 {3} | 5 {1,3} | 6 {2,3} | 7 {1,2,3} |
|---|---|---|---|---|---|---|---|---|
| 0 | x | 待求 | ||||||
| 1 | 已知 | x | C1 | x | B1 | x | A | x |
| 2 | 已知 | C2 | x | x | A1 | B | x | x |
| 3 | 已知 | B2 | A2 | C | x | x | x | x |
上表中的x表示无须计算, 其他字母表示依赖关系, 待求值依赖A, B, C, 而A依赖A1, A2, B, C类似. 因此, 填表时我们需要按 <span class="MathJax" id="MathJax-Element-26-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">j" role="presentation" style="margin:0px;padding:0px;display:inline;line-height:normal;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;position:relative;">j 的顺序来填, 这样才不会差生依赖问题.
为了输出路径, 我们还需要另一个同样大小的二维数组, 记录每次决定的路径.
下面的Fortran代码是根据资料中的C代码改写的, 放在这里供参考.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71 | Program TSP
integer:: Linp=1integer N, M, i, j, k
real*8:: RMAX=1.d10, Rmin, Rtmp
integer, allocatable:: P(:,:)
real*8, allocatable:: C(:,:), D(:,:)
open(Linp, file='C.dat')
read(Linp, *) N ! 读入节点个数
N = N-1! 邻接矩阵, 下标从0开始allocate(C(0:N, 0:N))
read(Linp, *) ( (C(i, j), j=0,N), i=0,N ) ! 读入数据
M=2**N-1! 确定集合大小, 分配数组. D存放阶段最优值, P存放最优策略allocate(D(0:N, 0:M), P(0:N, 0:M))
P =-1! 初值, 设为不可能的值
D =-1.d0
D(1:N, 0) = C(1:N, 0) ! 0列赋初值! 填表do j=1, M-1! 最后一列不在循环中计算, 节省时间do i=1, N
if(Iand(j, 2**(i-1))==0) then! i不属于集合j
Rmin=RMAX
do k=1, N
if(Iand(j, 2**(k-1))/=0) then! k属于集合j
Rtmp = C(i, k)+D(k, j-2**(k-1)) ! 从j中去掉k即将k对应的二进制位置0if(Rtmp<Rmin) thenRmin=Rtmp
D(i,j) = Rmin ! 阶段最优值
P(i,j) = k ! 最优决策end ifend ifend doend ifend doend do! 计算总体最优值
Rmin=RMAX
do k=1, N
Rtmp = C(0, k)+D(k, M-2**(k-1))
if(Rtmp<Rmin) thenRmin=Rtmp
D(0,M) = Rmin ! 整体最优值
P(0,M) = k ! 最优决策end ifend do! 输出最优路径
i=0; j=M
do while(j>0)
i = P(i, j) ! 下一步去往哪个结点
j = j-2**(i-1) ! 从j中去掉i结点print*, i
end do! 输出矩阵write(*, '(A, 1000I3)') 'ij', (j, j=0, M)
do i=0, N
write(*, '(1000I3)') i, (int(D(i,j)), j=0, M)
end dowrite(*, '(/, A, 1000I3)') 'ij', (j, j=0, M)
do i=0, N
write(*, '(1000I3)') i, (P(i,j), j=0, M)
end do |
◆本文地址: http://jerkwin.github.io/2016/03/17/旅行推销商问题TSP的动态规划解法/, 转载请注明◆
https://blog.sciencenet.cn/blog-548663-974691.html
上一篇：GROMACS分析教程：使用g_select计算平均滞留时间
下一篇：使用GROMACS和VMD绘制空间分布函数SDF的步骤
