Adding count constraint in lpSolveAPI R

I'm trying to dive into optimization and have been messing around with lpSolveAPI package. Example will demonstrate my simple setup.

data (every row holds different variable):

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)

create empty model, with 5 variables (rows in test.df) and I want to maximize my objective.

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

Lets add few constraints.

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)

I would like to create onemore constraint, which would be that use only x number of variables in solve. So if I set it up to 2, it would create optimal solution with 2 of those variables. I have tried set.type = "binary" but that wasn't successful.

Solution

Here some code that shows how to apply the technique mentioned by @ErwinKalvelagen to lpSolveAPI in R. The main point is to add additional binary variables to the problems.

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

Note, however, that the problem becomes infeasible if you require that only two x(i) are greater zero. You need at least three to meet the constraints given. (this is done by the parameter n). Also note that set.row is more efficient than add.constraint in the long term.

Elaborating @ErwinKalvelagen's second comment, another approach is to simply take every n out of the 5 possible variable combination and solve for these n variables. Translated to R code, this would become

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

Both codes return the same solution, however the first solution is much faster:

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