R Lpsolve Binary Find All Possible Solutions

How to get lpsolveAPI to return all possible solutions?

Discovered that this was dupe with:
Get multiple solutions for 0/1-Knapsack MILP with lpSolveAPI

Here is the code I used to solve this adapted from the link:

# find first solution
status<-solve(lpmodel) 
sols<-list() # create list for more solutions
obj0<-get.objective(lpmodel) # Find values of best solution (in this case four)
counter <- 0 #construct a counter so you wont get more than 100 solutions

# find more solutions
while(counter < 100) {
  sol <- get.variables(lpmodel)
  sols <- rbind(sols,sol)
  add.constraint(lpmodel,2*sol-1,"<=", sum(sol)-1) 
  rc<-solve(lpmodel)
  if (status!=0) break;
  if (get.objective(lpmodel)<obj0) break;
  counter <- counter + 1
}
sols

This finds all six solutions:

> sols
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
sol 1    0    0    1    0    1    0    1   
sol 1    0    0    1    0    1    1    0   
sol 0    1    0    1    0    1    0    1   
sol 0    1    0    1    0    1    1    0   
sol 0    1    0    0    1    1    1    0   
sol 0    1    0    0    1    1    0    1   

Sorry for the dupe everyone. If anyone knows another way built into the 'lpSolveAPI' package I'd still love to know.

All feasible solutions using lpSolveAPI in R

There are two different approaches to enumerate all feasible integer solutions.

The first is to add a cut to the constraints that forbids the previously found solution. I.e.

 1. solve problem
 2. if infeasible: done
 3. record optimal solution x[i] into a[i]
 4. add cut of the form
      sum(i, (2*a[i]-1)*x[i]) <= sum(i,a[i])-1
 5. go to step 1. 

I assumed here the binary variables are x[i].

The second approach is to use the Cplex solution pool. This is much faster if there are many solutions.

PS: LpSolveApi documents things like get.solutioncount and select.solution. I don't think that actually works.

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

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

---

Related Topics