R中带条件约束的线性规划

Question

我有一个线性编程问题,我试图从一些二进制资源中选择优化值,基本上是一个背包问题.我遇到的问题是不同的资源具有共同的特征,我想确保我的最终解决方案具有0或2个具有特定特征的资源.有没有办法实现这个目标？尽管进行了广泛的搜索,我仍然无法想到一个或找到一个.在我的数据中,决策变量是资源,约束是这些资源的特征.请考虑以下代码:

library(lpSolve)
const_mat_so<-matrix(c(
  c(0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,1,0,0,1,0,1)
     ,c(0,  0,  0,  0,  0,  0,  0,  1,  0,  1,  0,1,1,0,0,1,1)
     ,c(0,  0,  0,  0,  0,  1,  1,  0,  0,  0,  1,0,1,0,1,0,0)
     ,c(1,  1,  0,  1,  1,  0,  0,  0,  1,  0,  0,0,0,0,0,0,0)
     ,c(8800,   8500,   7600,   8600,   8400,   7500, 7000, 8500,   8800,   7700,   6700,5500,1200,6700,9500,8700,6500)
     ,c(0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,0,0,0,0,0,0)
     ,c(0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,0,0,0,0,0,0)
     ,c(0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,0,0,0,0,0,0)
     ,c(0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  1,0,0,1,0,1,0)
     ,c(0,  0,  1,  0,  0,  0,  0,  0,  1,  0,  0,1,1,0,0,0,0)
     ,c(0,  0,  1,  0,  0,  0,  0,  1,  0,  0,  0,0,0,0,0,0,0)
     ,c(0,  0,  0,  1,  0,  0,  0,  0,  0,  1,  0,1,1,1,0,1,0)
     ,c(0,  1,  0,  0,  0,  0,  0,  1,  0,  0,  0,0,0,0,0,1,0)
     ,c(0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,0,0,0,1,0,0)
     ,c(0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  1,0,0,0,0,0,0)
     ),nrow=15,byrow = TRUE)

const_dir_so<-c("=","=","=","=","<=","<=","<=","<=","<=","<=","<=","<=","<=","<=","<=")

max_cost_so = 25000

objective_so = c(21.0, 19.3, 19.2, 18.8, 18.5, 16.6, 16.4, 16.4, 16.0, 16.0, 14.9, 14.6, 14.0, 13.9,12.0,5.5,24.6)

const_rhs_so<-c(1,1,1,1,25000,3,3,3,2,2,2,2,2,2,2)

x = lp ("max", objective_so, const_mat_so, const_dir_so, const_rhs_so,     all.bin=TRUE, all.int=TRUE 
    )

> x
Success: the objective function is 68.1

> x$solution
 [1] 1 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0

虽然上面提出了一个解决方案,但它不是我想要的解决方案,因为我实际上希望最后七个约束是> = 2或0.我不知道如何编码或者是否可能.任何帮助,将不胜感激.我不是一个线性编程高手,所以请原谅任何关于这种方法的误解.

Answer 1

我的理解是最后 7 个约束中的每一个都大于 2 或等于 0，即不是 1。

1) 这样的约束只有 7 个，因此有 2^7 = 128 种可能性，这种可能性足够小，我们可以使用给定问题的公式来运行每一种可能性，而无需过多的运行时间，然后找到其中的最大值。

dec2bin接受一个以 10 为基数（即十进制）的数字，并将其转换为由 0 和 1 组成的二进制向量。对 0 到 127 之间的每个数字运行它会给出二进制数，其中 1 对应于 >= 2 的约束（其余等于 0）。

dec2bin <- function(dec, digits = 7) {
   # see http://stackoverflow.com/questions/6614283/converting-decimal-to-binary-in-r
   tail(rev(as.integer(intToBits(dec))), digits)
}

runLP <- function(i) {
   bin <- dec2bin(i)
   n <- length(const_rhs_so)  # 15
   ix <- seq(to = n, length = length(bin))  # indexes of last 7 constraints, i.e. 9:15
   const_dir_so[ix] <- ifelse(bin, ">=", "=")
   const_rhs_so[ix] <- 2*bin
   lp("max", objective_so, const_mat_so, const_dir_so, const_rhs_so, all.bin = TRUE)
}

lpout <- lapply(0:127, runLP)
ixmax <- which.max(sapply(lpout, "[[", "objval"))
ans <- lpout[[ixmax]]
ans
ans$solution
tail(c(const_mat_so %*% ans$solution), 7)

给予：

> ans
Success: the objective function is 62 
> ans$solution
 [1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1
> tail(c(const_mat_so %*% ans$solution), 7) # last 7 constraint values
[1] 0 0 0 0 0 0 0

2） 在@Erwin Kalvelagen的第二个替代方案中，它指的是约束变量，但我认为他的答案的意思x是最后7个约束之一的LHS值。也就是说，如果C是原来最后 7 个约束的矩阵，那么用这 14 个约束替换原来的 7 个约束：

 Cx + D1 y <= 0
 Cx + D2 y >= 0

其中D1是对角矩阵，其对角元素是任意足够大的负数，D2是对角矩阵，其对角元素全部为-2。在这里，我们正在优化二元变量的x向量y。变量x如问题中所示，有 7 个新的y二元变量，使得 y[i] 为 0，以将最后 7 个原始约束的第 i 个约束为 0 或 1，以将其约束为 2 或更多。这些变量在（1）中y被调用。目标中变量bin的系数均为零。y

就 lpSolve R 代码而言：

objective_so2 <- c(objective_so, numeric(7))
const_mat_so2 <- cbind(rbind(const_mat_so, const_mat_so[9:15, ]), 
                       rbind(matrix(0, 8, 7), diag(-100, 7), diag(-2, 7)))
const_dir_so2 <- c(const_dir_so, rep(">=", 7))
const_rhs_so2 <- c(const_rhs_so[1:8], numeric(14))
x2 = lp ("max", objective_so2, const_mat_so2, const_dir_so2, const_rhs_so2, all.bin = TRUE)

给出与 (1) 中相同的值 62。变量y（最后 7 个）均为 0，也对应于 (1)。这也提供了双重检查，因为两种方法现在给出了一致的答案。

> x2
Success: the objective function is 62 
> x2$solution
 [1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0