在 lpSolveAPI R 中添加计数约束

Question

我正在尝试深入研究优化，并且一直在摆弄 lpSolveAPI 包。示例将演示我的简单设置。

数据（每一行都有不同的变量）：

dput(test.df)
    structure(list(objective = c(-0.162235888601422, -0.168597233981057, 
    -0.165558234725657, -0.156096491294958, -0.15294764940114), constrain1 = c(1.045, 
    1.259, 1.792, 2.195, 2.802)), .Names = c("objective", "constrain1"
    ), row.names = c(NA, -5L), class = "data.frame")

library(lpSolveAPI)

创建空模型，有 5 个变量（test.df 中的行），我想最大化我的 objective。

test.lp <- make.lp(0, NROW(test.df))
set.objfn(test.lp, test.df$objective)
lp.control(test.lp,sense='max')

让我们添加一些约束。

add.constraint(test.lp, test.df$constrain1, "=", 2)
add.constraint(test.lp, rep(1, nrow(test.df)), "=", 1)
set.bounds(test.lp, upper = rep(0.5, nrow(test.df)))
set.bounds(test.lp, lower = rep(0, nrow(test.df)))
RowNames <- c("constraint1", "constraint2")
ColNames <- paste0("var", seq(1, nrow(test.df)))
dimnames(test.lp) <- list(RowNames, ColNames)

我想再创建一个约束，即在求解时仅使用 x 个变量。因此，如果我将其设置为 2，它将创建具有其中 2 个变量的最佳解决方案。我试过 set.type = "binary" 但没有成功。

Answer 1

我认为您想添加一个限制条件，将非零变量的数量 x(i) 限制为 2。在 LP 中无法真正完成计数，但可以将其表示为 MIP。

标准公式是引入二元变量 y(i)，其中：

x(i) <= y(i)*xup(i)      (implication: y(i)=0 => x(i)=0)
sum(i, y(i)) <= 2     
0 <= x(i) <= xup(i)      (bounds)
y(i) in {0,1}            (binary variables)

对于更大的问题，这比解决每个可能的组合更有效。尽管 N 中的 k=2 还不错：N choose k = N*(N-1)/2 种可能的组合。

Answer 2

这里有一些代码展示了如何将@ErwinKalvelagen 提到的技术应用于 R 中的 lpSolveAPI。要点是为问题添加额外的二进制变量。

library(lpSolveAPI)
fun1 <- function(n=3) {
    nx <- 5

    # set up lp
    lprec <- make.lp(0, 2*nx) # one var for the value x(i) and one var y(i) binary if x(i) > 0
    set.objfn(lprec, c(-0.162235888601422, -0.168597233981057, -0.165558234725657, -0.156096491294958, -0.15294764940114, rep(0, nx)))
    lp.control(lprec,sense='max')
    set.type(lprec, columns=seq(nx+1,2*nx), "binary") # y(i) are binary vars

    # add constraints as in the question
    set.bounds(lprec, upper=rep(0.5, nx), columns=seq(1,nx)) # lpsolve implicitly assumes x(i) >= 0, so no need to define bounds for that
    add.constraint(lprec, c(1.045, 1.259, 1.792, 2.195, 2.802, rep(0, nx)), "=", 2)
    add.constraint(lprec, c(rep(1, nx), rep(0, nx)), "=", 1)

    # additional constraints as suggested by @ErvinKarvelagen
    for (i in seq(1,nx)) add.constraint(lprec, xt=c(1, -0.5), type="<=", rhs=0, indices=c(i, nx+i)) # x(i)<=y(i)*0.5
    add.constraint(lprec, c(rep(0,nx), rep(1,nx)), "<=", n) # sum(y(i))<=2 (if set to 3, it finds a solution)

    # solve and print solution
    status <- solve(lprec)
    if(status!=0) stop("no solution found, error code=", status)
    sol <- get.variables(lprec)[seq(1, nx)]
    return(sol)
}

但是请注意，如果您要求只有两个 x(i) 大于零，则该问题变得不可行。您至少需要三个才能满足给定的限制。（这是由参数 n 完成的）。另请注意，从长远来看，set.row 比 add.constraint 更有效。

阐述@ErwinKalvelagen 的第二条评论，另一种方法是简单地从 5 个可能的变量组合中取出每个 n 并求解这 n 个变量。翻译成 R 代码，这将变成

fun2 <- function(n=3) {
    nx <- 5

    solve_one <- function(indices) {
        lprec <- make.lp(0, n) # only n variables
        lp.control(lprec,sense='max')
        set.objfn(lprec, c(-0.162235888601422, -0.168597233981057, -0.165558234725657, -0.156096491294958, -0.15294764940114)[indices])
        set.bounds(lprec, upper=rep(0.5, n))
        add.constraint(lprec, c(1.045, 1.259, 1.792, 2.195, 2.802)[indices],"=", 2)
        add.constraint(lprec, rep(1, n), "=", 1)
        status <- solve(lprec)
        if(status==0) 
            return(list(obj = get.objective(lprec), ind=indices, sol=get.variables(lprec)))
        else
            return(list(ind=indices, obj=-Inf))
    }

    sol <- combn(nx, n, FUN=solve_one, simplify=FALSE)
    best <- which.max(sapply(sol, function(x) x[["obj"]]))
    return(sol[[best]])
}

两个代码 return 相同的解决方案，但是第一个解决方案要快得多：

library(microbenchmark)
microbenchmark(fun1(), fun2(), times=1000, unit="ms")
#Unit: milliseconds
#   expr      min       lq      mean   median       uq       max neval
# fun1() 0.429826 0.482172 0.5817034 0.509234 0.563555  6.590409  1000
# fun2() 2.514169 2.812638 3.2253093 2.966711 3.202958 13.950398  1000

在 lpSolveAPI R 中添加计数约束

Adding count constraint in lpSolveAPI R

optimization

r

lpsolve