仓库源文，站点原文

layout: post title: 旅行推销商问题TSP的动态规划解法 categories:

科 tags:
数理
fortran viz: true math: true

2016-03-17 20:35:36

利用动态规划方法是可以精确求解旅行推销商问题(Traveling Salesman Problem)的, 虽然这种方法只适用于求解小规模的问题. 这个算法我一直没有弄清楚, 最近有个问题需要使用类似的算法来解决, 所以我就仔细研究了一下这个算法. 下面是网上几篇资料的总结.

我们先考虑一个小规模的问题, 只有4个城市, 城市间的距离由下面的矩阵 $C$ 决定.

i\j	0	1	2	3
0	0	3	6	7
1	5	0	2	3
2	6	4	0	2
3	3	7	5	0

值得注意的是, 这个矩阵是不对称的, 也就是说 $C{ij} \neq C{ji}$. 如果这个矩阵是对称的, 算法还可以简化一些.

假定我们从城市0出发, 城市1, 2, 3每个经过一次, 最后回到城市0, 那么求解的递归树可以表示如下:

<pre style="display:none"><textarea name="viz" id="viz-1"> digraph TSP4 { rankdir=TB; N0 [label="d(0,{1,2,3}) \n d[0,7]=3+7=10"] N10 [label="d(1,{}) \n d[3,0]=5"] N20 [label="d(2,{}) \n d[2,0]=6"] N30 [label="d(3,{}) \n d[3,0]=3"] N11 [label="d(1,{2,3}) \n d[1,6]=2+5=7"] N12 [label="d(2,{1,3}) \n d[2,5]=4+6=10"] N13 [label="d(3,{1,2}) \n d[3,3]=5+9=14"] N21 [label="d(2,{3}) \n d[2,4]=2+3=5"] N22 [label="d(3,{2}) \n d[3,2]=5+6=11"] N23 [label="d(1,{3}) \n d[1,4]=3+3=6"] N24 [label="d(3,{1}) \n d[3,1]=7+5=12"] N25 [label="d(1,{2}) \n d[1,2]=2+6=8"] N26 [label="d(2,{1}) \n d[2,1]=4+5=9"] N0 -> N11 [label="3"] N0 -> N12 [label="6"] N0 -> N13 [label="7"] N11 -> N21 [label="2"] N11 -> N22 [label="3"] N12 -> N23 [label="4"] N12 -> N24 [label="2"] N13 -> N25 [label="7"] N13 -> N26 [label="5"] N21 -> N30 [label="2"] N22 -> N20 [label="5"] N23 -> "d(3,{}) \n d[3,0]=3" [label="3"] N24 -> N10 [label="7"] N25 -> "d(2,{}) \n d[2,0]=6" [label="2"] N26 -> "d(1,{}) \n d[1,0]=5" [label="4"] } </textarea></pre><figure id="fig-viz-1"><figurecaption>Fig.1</figurecaption></figure>

其中的 $d(i,\{j,k,l\})$ 表示由城市 $i$ 出发, 集合 $\{j,k,l\}$ 中的城市 $j,k,l$ 每个经过一次需要的最小路程, 箭头上的值表示两个城市之间的距离. 很显然, $d(0, \{1,2,3\})$ 就是我们最终要求的值. 这个值可以一步一步分解下去, 最终分解为每个城市到城市0的距离, $d(i,\{\})$. 将过程反过来, 向上递推, 就可以得到需要的值了.

为了方便求解并记录路径, 我们可以使用二进制数表示城市集合. 一个 $n$ 位的二进制数可以表示 $n$ 个城市的集合. 当某位为1时, 表示这个位所代表的城市出现在集合中. 我们最多有3个城市需要经历, 所以需要 $2^3=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}

在上面图中, $d[i,j]$ 中的 $j$ 使用的就是表中对应的集合. 根据递推关系

$$d[i,j]=\begin{cases} C{i0}, & j=0 \ \ \underset { {k }\in j} \min { C{ik} +d[k,j \backslash {k}] }=\underset {{k}\in j} \min { C_{ik} +d[k,j-2^{k-1}] }, & j \neq 0 \end{cases}$$

其中的 $j \backslash \{k\}$ 表示从集合 $j$ 中去除集合 $\{k\}$. 集合 $j$ 去除城市 $k$ 后所形成的集合, 其对应的编号为 $j-2^{k-1}$. 判断城市 $k$ 是否属于集合 $j$ 的方法是使用二进制与操作, 若j & 2^(k-1) == 0, $k$ 不属于集合 $j$.

在写代码实现时, 我们可以使用一个二维表格(数组), 其大小为 $n 2^{n-1}$, 根据上面的关系式逐次填充计算值, 最终得到结果.

i\j	0<br/>{}	1<br/>{1}	2<br/>{2}	3<br/>{1,2}	4<br/>{3}	5<br/>{1,3}	6<br/>{2,3}	7<br/>{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类似. 因此, 填表时我们需要按 $j$ 的顺序来填, 这样才不会差生依赖问题.

为了输出路径, 我们还需要另一个同样大小的二维数组, 记录每次决定的路径.

下面的Fortran代码是根据资料中的C代码改写的, 放在这里供参考.

<table class="highlighttable"><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 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</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">Program TSP integer:: Linp=1 integer 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对应的二进制位置0 if(Rtmp<Rmin) then Rmin=Rtmp D(i,j) = Rmin ! 阶段最优值 P(i,j) = k ! 最优决策 end if end if end do end if end do end do ! 计算总体最优值 Rmin=RMAX do k=1, N Rtmp = C(0, k)+D(k, M-2**(k-1)) if(Rtmp<Rmin) then Rmin=Rtmp D(0,M) = Rmin ! 整体最优值 P(0,M) = k ! 最优决策 end if end 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)') 'i\j', (j, j=0, M) do i=0, N write(*, '(1000I3)') i, (int(D(i,j)), j=0, M) end do write(*, '(/, A, 1000I3)') 'i\j', (j, j=0, M) do i=0, N write(*, '(1000I3)') i, (P(i,j), j=0, M) end do </pre></div> </td></tr></table>