3

我想在 R 中编写旅行商问题。我将首先从 3 个城市开始,然后我将扩展到更多城市。下面的距离矩阵给出了 3 个城市之间的距离。目标(如果有人不知道)是推销员将从一个城市开始,并将访问其他 2 个城市,这样他必须走最短距离。

在以下情况下,他应该从纽约或洛杉矶出发,然后前往芝加哥,然后前往其余城市。我需要帮助来定义 A_ (我的约束矩阵)。

我的决策变量将与距离矩阵具有相同的维度。它将是一个 1,0 矩阵,其中 1 表示从等于行名的城市到等于列名的城市的旅行。例如,如果一个推销员从纽约旅行到芝加哥,第 1 行中的第二个元素将是 1。我的列名和行名是 ny、chicago 和 LA

通过查看问题的解决方案,我得出结论,我的约束将是:

行总和必须小于 1,因为他不能从同一个城市离开两次

列总和必须小于 1,因为他不能两次进入同一个城市

矩阵元素的总和必须为 2,因为销售员将访问 2 个城市并从 2 个城市离开。

我需要帮助来定义 A_ (我的约束矩阵)。我应该如何将我的决策变量与约束联系起来?

 ny=c(999,9,20)
 chicago=c(9,999,11)
 LA=c(20,11,999)
 distances=cbind(ny,chicago,LA)


 dv=matrix(c("a11","a12","a13","a21","a22","a23","a31","a32","a33"),nrow=3,ncol=3)

 c_=c(distances[1,],distances[2,],distances[3,])
 signs = c((rep('<=', 7)))
 b=c(1,1,1,1,1,1,2)
 res = lpSolve::lp('min', c_, A_, signs, b,  all.bin = TRUE)
4

2 回答 2

2

您的解决方案存在一些问题。首先是您考虑的限制并不能保证所有城市都会被访问——例如,路径可能只是从 NY 到 LA,然后再返回。这可以很容易地解决,例如,要求每一行和每一列的总和恰好为 1,而不是最多为 1(尽管在这种情况下,您会找到一个旅行推销员之旅,而不仅仅是一条路径)。

更大的问题是,即使我们解决了这个问题,您的约束也不能保证所选顶点实际上在图中形成一个循环,而不是多个较小的循环。而且我认为您对问题的陈述不能解决这个问题。

这是使用 LP 的 Traveling Salesman 的实现。解空间的大小为 n^3,其中 n 是距离矩阵中的行数。这表示 nxn 矩阵的 n 个连续副本,每个副本表示在时间t为遍历的边1<=t<=n。约束保证

  1. 每一步最多遍历一条边
  2. 曾经的顶点只被访问过一次
  3. 遍历的第 i 条边的起点与第 i-1 条边的终点相同

这避免了多个小循环的问题。例如,如果有四个顶点,则该序列(12)(21)(34)(43)不是有效的解,因为第二条边的端点与第三条边(21)的起点不匹配(34)

tspsolve<-function(x){
   diag(x)<-1e10
   ## define some basic constants
   nx<-nrow(x)
   lx<-length(x)
   objective<-matrix(x,lx,nx)
   rowNum<-rep(row(x),nx)
   colNum<-rep(col(x),nx)
   stepNum<-rep(1:nx,each=lx)

   ## these constraints ensure that at most one edge is traversed each step
   onePerStep.con<-do.call(cbind,lapply(1:nx,function(i) 1*(stepNum==i)))
   onePerRow.rhs<-rep(1,nx)

   ## these constraints ensure that each vertex is visited exactly once
   onceEach.con<-do.call(cbind,lapply(1:nx,function(i) 1*(rowNum==i)))
   onceEach.rhs<-rep(1,nx)

   ## these constraints ensure that the start point of the i'th edge
   ## is equal to the endpoint of the (i-1)'st edge
   edge.con<-c()
   for(s in 1:nx){
     s1<-(s %% nx)+1    
     stepMask<-(stepNum == s)*1
     nextStepMask<- -(stepNum== s1)
     for(i in 1:nx){        
       edge.con<-cbind(edge.con,stepMask * (colNum==i) + nextStepMask*(rowNum==i))
     }
   }
   edge.rhs<-rep(0,ncol(edge.con))

   ## now bind all the constraints together, along with right-hand sides, and signs
   constraints<-cbind(onePerStep.con,onceEach.con,edge.con)
   rhs<-c(onePerRow.rhs,onceEach.rhs,edge.rhs)
   signs<-rep("==",length(rhs))
   list(constraints,rhs)

   ## call the lp solver
   res<-lp("min",objective,constraints,signs,rhs,transpose=F,all.bin=T)

   ## print the output of lp
   print(res)

   ## return the results as a sequence of vertices, and the score = total cycle length
   list(cycle=colNum[res$solution==1],score=res$objval)
}

这是一个例子:

set.seed(123)
x<-matrix(runif(16),c(4,4))
x
##           [,1]      [,2]      [,3]      [,4]
## [1,] 0.2875775 0.9404673 0.5514350 0.6775706
## [2,] 0.7883051 0.0455565 0.4566147 0.5726334
## [3,] 0.4089769 0.5281055 0.9568333 0.1029247
## [4,] 0.8830174 0.8924190 0.4533342 0.8998250
tspsolve(x)
## Success: the objective function is 2.335084 
## $cycle
## [1] 1 3 4 2
## 
## $score
## [1] 2.335084

我们可以使用原始的蛮力搜索来检查这个答案的正确性:

tspscore<-function(x,solution){
    sum(sapply(1:nrow(x), function(i) x[solution[i],solution[(i%%nrow(x))+1]])) 
}

tspbrute<-function(x,trials){
  score<-Inf
  cycle<-c()
  nx<-nrow(x)
  for(i in 1:trials){
    temp<-sample(nx)
    tempscore<-tspscore(x,temp)
    if(tempscore<score){
      score<-tempscore
      cycle<-temp
    }
  }
  list(cycle=cycle,score=score)
}

tspbrute(x,100)
## $cycle
## [1] 3 4 2 1
## 
## $score
## [1] 2.335084

请注意,即使这些答案名义上不同,它们也代表相同的周期。

但是,对于较大的图表,蛮力方法不起作用:

> set.seed(123)
> x<-matrix(runif(100),10,10)
> tspsolve(x)
Success: the objective function is 1.296656 
$cycle
 [1]  1 10  3  9  5  4  8  2  7  6

$score
[1] 1.296656

> tspbrute(x,1000)
$cycle
 [1]  1  5  4  8 10  9  2  7  6  3

$score
[1] 2.104487

这种实现对于小矩阵非常有效,但正如预期的那样,随着它们变大,它开始严重恶化。在大约 15x15 时,它开始减速很多:

timetsp<-function(x,seed=123){
    set.seed(seed)
    m<-matrix(runif(x*x),x,x)   
    gc()
    system.time(tspsolve(m))[3]
}

sapply(6:16,timetsp)
## elapsed elapsed elapsed elapsed elapsed elapsed elapsed elapsed elapsed elapsed 
## 0.011   0.010   0.018   0.153   0.058   0.252   0.984   0.404   1.984  20.003 
## elapsed 
## 5.565
于 2013-10-05T00:04:04.697 回答
1

您可以使用该gaoptim包来解决排列/实值问题 - 它是纯 R,所以它不是那么快:

欧元旅游问题(见?optim)

 eurodistmat = as.matrix(eurodist)

 # Fitness function (we'll perform a maximization, so invert it)
 distance = function(sq)
 {
   sq = c(sq, sq[1])
   sq2 <- embed(sq, 2)
   1/sum(eurodistmat[cbind(sq2[,2], sq2[,1])])
 }

 loc = -cmdscale(eurodist, add = TRUE)$points
 x = loc[, 1]
 y = loc[, 2]
 n = nrow(eurodistmat)

 set.seed(1)

 # solving code
 require(gaoptim)
 ga2 = GAPerm(distance, n, popSize = 100, mutRate = 0.3)
 ga2$evolve(200)
 best = ga2$bestIndividual()
 # solving code

 # just transform and plot the results
 best = c(best, best[1])
 best.dist = 1/max(ga2$bestFit())
 res = loc[best, ]
 i = 1:n

 plot(x, y, type = 'n', axes = FALSE, ylab = '', xlab = '')
 title ('Euro tour: TSP with 21 cities')
 mtext(paste('Best distance found:', best.dist))
 arrows(res[i, 1], res[i, 2], res[i + 1, 1], res[i + 1, 2], col = 'red', angle = 10)
 text(x, y, labels(eurodist), cex = 0.8, col = 'gray20')
于 2013-10-09T14:20:06.877 回答