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の行）で、空のモデルを作成し、私は自分の目標を最大化します。

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) 

私はonemore制約を作成したいと思います。これは、x個の変数のみを使用して解決します。だから私は2に設定すると、2つの変数で最適解を作り出すでしょう。私はset.type = "binary"を試しましたが、それは成功しませんでした。

Answer 1

@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) 
} 

しかし、2つのx（i）がゼロより大きい必要がある場合、問題は実行不可能になります。与えられた制約を満たすには、少なくとも3つが必要です。 （これはパラメータnによって行われる）。また、set.rowは、長期的にはadd.constraintよりも効率的です。

@ ErwinKalvelagenの第2のコメントを精緻化する別のアプローチは、可能な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]]) 
} 

なる両コードは同じソリューションを返す、しかし、最初のソリューションは、はるかに高速です：

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 

Answer 2

0以外の変数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はそれほど悪くはありません：N choose k = N*(N-1)/2可能な組み合わせ。