Добавление ограничения количества в lpSolveAPI R

Вот некоторый код, который показывает, как применить технику, упомянутую @ErwinKalvelagen, к lpSolveAPI в R. Суть в том, чтобы добавить к задачам дополнительные бинарные переменные.

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, другой подход состоит в том, чтобы просто взять каждое n из 5 возможных комбинаций переменных и решить для этих 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

Я думаю, вы хотите добавить ограничение, которое ограничивает количество ненулевых переменных 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) 

Для более крупных задач это может быть намного эффективнее, чем решение всех возможных комбинаций. Хотя k=2 из N не так уж и плохо: N choose k = N*(N-1)/2 возможных комбинаций.