R lpsolve如何定义约束旅行商

Question

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

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

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

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

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

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

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

我需要帮助来定义 A_(我的约束矩阵)。我应该如何将决策变量绑定(bind)到约束中？

 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)

Answer 1

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

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

下面是使用 LP 的旅行推销员的实现。解空间的大小为 n^3，其中 n 是距离矩阵中的行数。这表示 nxn 矩阵的 n 个连续副本，每个副本表示在时间 t for 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