旅行推销商问题TSP的动态规划解法

类别:    标签: 数理 fortran   阅读次数:   版权: (CC) BY-NC-SA

利用动态规划方法是可以精确求解旅行推销商问题(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, 那么求解的递归树可以表示如下:

</br>Fig.1

其中的 $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$ 使用的就是表中对应的集合. 根据递推关系

其中的 $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
随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: GROMACS分析教程:使用g_select计算平均滞留时间
后一篇: 使用GROMACS和VMD绘制空间分布函数SDF的步骤

访问人次(2015年7月 9日起): | 最后更新: 2017-12-13 23:46:45 UTC | 版权所有 © 2008 - 2017 Jerkwin