Get multiple solutions for 0/1-Knapsack MILP with lpSolveAPI

Question

Reproducable Example:

I described a simple 0/1-Knapsack problem with lpSolveAPI in R, which should return 2 solutions:

library(lpSolveAPI)

lp_model= make.lp(0, 3)
set.objfn(lp_model, c(100, 100, 200))
add.constraint(lp_model, c(100,100,200), "<=", 350)
lp.control(lp_model, sense= "max")
set.type(lp_model, 1:3, "binary")

lp_model
solve(lp_model)
get.variables(lp_model)
get.objective(lp_model)
get.constr.value((lp_model))
get.total.iter(lp_model)
get.solutioncount(lp_model)

Problem:

But get.solutioncount(lp_model) shows that there's just 1 solution found:

> lp_model
Model name: 
           C1   C2   C3         
Maximize  100  100  200         
R1        100  100  200  <=  350
Kind      Std  Std  Std         
Type      Int  Int  Int         
Upper       1    1    1         
Lower       0    0    0         
> solve(lp_model)
[1] 0
> get.variables(lp_model)
[1] 1 0 1
> get.objective(lp_model)
[1] 300
> get.constr.value((lp_model))
[1] 350
> get.total.iter(lp_model)
[1] 6
> get.solutioncount(lp_model)
[1] 1

I would expect that there are 2 solutions: 1 0 1 and 0 1 1.

I tried to pass the num.bin.solns argument of lpSolve with solve(lp_model, num.bin.solns=2), but the number of solutions remained 1.

Question:

How can I get the two correct solutions? I prefer using lpSolveAPI as the API is really nice. If possible I'd like to avoid to use lpSolve directly.

Answer 1

Looks like that is broken. Here is a DIY approach for your specific model:

# first problem
rc<-solve(lp_model)
sols<-list()
obj0<-get.objective(lp_model)
# find more solutions
while(TRUE) {
   sol <- round(get.variables(lp_model))
   sols <- c(sols,list(sol))
   add.constraint(lp_model,2*sol-1,"<=", sum(sol)-1)
   rc<-solve(lp_model)
   if (rc!=0) break;
   if (get.objective(lp_model)<obj0-1e-6) break;
}
sols

The idea is to cut off the current integer solution by adding a constraint. Then resolve. Stop when no longer optimal or when the objective starts to deteriorate. Here is some math background.

You should see now:

> sols
[[1]]
[1] 1 0 1

[[2]]
[1] 0 1 1

Update

Below in the comments it was asked why coefficients of the cut have the form 2*sol-1. Again have a look at the derivation. Here is a counter example:

           C1   C2        
Maximize    0   10        
R1          1    1  <=  10
Kind      Std  Std        
Type      Int  Int        
Upper       1    1        
Lower       0    0       

With "my" cuts this will yield:

> sols
[[1]]
[1] 0 1

[[2]]
[1] 1 1

while using the suggested "wrong" cuts will give just:

> sols
[[1]]
[1] 0 1