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
