使用 lpsolve 在 R 中进行线性规划

Question

我正在尝试使用lpsolve 包解决R 中的线性规划问题。

问题来了：

这里是 R 中的可重现示例的示例：

    library("lpSolve")
a <- matrix(c(1,2,5,
              1/2,1,3,
              1/5,1/3,1),nrow=3,byrow=T)

#
f.obj <- c(1,0,0,0)

f.con <- matrix (c(
  1,1,-a[1,2],0, #Contraint 1 for a12
  1,-1,a[1,2],0, #Contraint 2 for a12
  1,1,0,-a[1,3], #Contraint 1 for a13
  1,-1,0,a[1,3], #Contraint 2 for a13
  1,0,1,-a[2,3], #Contraint 1 for a23
  1,0,-1,a[2,3], #Contraint 2 for a23
  0,1,1,1, #Contraint 3
  0,1,0,0, #Constraint 4
  0,0,1,0, #Constraint 4
  0,0,0,1 #Constraint 4

  ), nrow=10, byrow=TRUE)

f.dir <- c(rep("<=",6), "=",rep(">",3))

f.rhs <- c(rep(1,6),1,rep(0,3))

g <- lp ("max", f.obj, f.con, f.dir, f.rhs)
g$solution

如果我有一个7 X 7 或n x n 矩阵a，我可以手动解决这个小问题。我将如何指定约束 1 和 2，尤其是我正在努力定义与 a[i,j] 相关的约束？

a = matrix( 
  c(1,4,9,6,6,5,5,
    1/4,1,7,5,5,3,4,
    1/9,1/7,1,1/5,1/5,1/7,1/5,
    1/6,1/5,5,1,1,1/3,1/3,
    1/6,1/5,5,1,1,1/3,1/3,
    1/5,1/3,7,3,3,1,2,
    1/5,1/4,5,3,3,1/2,1
  ),nrow = 7,byrow =T)

上述矩阵的解决方案是0.986 0.501 0.160 0.043 0.060 0.060 0.1 0.075 任何帮助将不胜感激。

Answer 1

这是使用 for 循环的一种可能性。

正如我在评论中提到的，我认为您的条件 (4) 错误。这是我的建议。 我的想法是首先为约束（4）设置一个矩阵，然后为约束（3）设置一个矩阵 然后在循环中添加约束 (2) 和 (1)。请注意，在开始时，我不考虑与 \mu 对应的列。我将在最后添加此列。

n<- nrow(a)
f.cons<- diag(n)
f.cons<- rbind(f.cons, rep(1,n))

这会建立一个对应于约束 (4)（前 n 行）和约束 (3) 的矩阵。现在我向这个矩阵添加行，使用循环和命令 rbind.

for(i in 1:(n-1)){
  for(j in (i+1): n){
   x<- rep(0, n)
   x[i]<- 1  #x corresponds to (1)
   x[j]<- -a[i,j]
   y<- -x #y corresponds to (2)
   f.cons<- rbind(f.cons, rbind(x, y))
}

}

到目前为止，我已经忽略了第一列，它对应于 \mu。 我用这两条简单的线添加它：

 f.cons<- cbind(rep(1, nrow(f.cons)), f.cons)
 f.cons[1:(n+1), 1]=0

请注意，在我的矩阵 f.cond 中，前 n+1 行对应于约束 (3) 和 (4)！

Answer 2

已更新以包含修订后的约束 4，并进行了一些小的代码改进。

假设问题中的约束矩阵是正确的，这将使用combn 遍历所有 i x[1] 是i 的值，x[2] 是f 中j 的值。 make_cons 以与问题中所示相同的顺序返回约束矩阵，但如果可以使用这样的顺序，make_cons 中的 rbind 行可以简化为 rbind(cons1, cons2, cons3, cons4)。

make_cons <- function(a) {
   n <- nrow(a)
   f <- function(x) replace(numeric(n), x, c(1, -a[x[1], x[2]]))
   cons1 <- cbind(1, t(combn(1:n, 2, f)))
   cons2 <- cbind(1, -cons1[, -1])
   cons3 <- c(0, rep(1, n))
   cons4 <- cbind(0, diag(n))
   rbind(t(matrix(rbind(t(cons1), t(cons2)), ncol(cons1))), cons3, cons4)
}

# test

# a and f.con from question

a <- matrix(c(1, 0.5, 0.2, 2, 1, 0.333333333333333, 5, 3, 1), 3)
f.con <- matrix(c(1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 
  1, 1, 0, 0, -2, 2, 0, 0, 1, -1, 1, 0, 1, 0, 0, 0, -5, 5, -3, 
  3, 1, 0, 0, 1), 10)

all.equal(f.con, make_cons(a), check.attributes = FALSE)
## [1] TRUE