- 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, 那么求解的递归树可以表示如下:
其中的 $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 {} |
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类似. 因此, 填表时我们需要按 $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=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
|
