I am using LpSolveAPI package to solve a linear programming problem in R. This gives me one solution, but I would like to know if there are other solutions. Here is an example. Suppose I want to solve the following:
library(lpSolveAPI)
MyMatrix <- c(1,0,0,0,0,
0,0,0,1,0,
0,1,1,0,1,
0,1,0,1,1,
1,0,0,0,0,
0,0,1,0,1)
MyMatrix <- matrix(MyMatrix,nrow=6,ncol=5,byrow=TRUE)
ObjFun=c(50,30,100,100,200)
lprec <- make.lp(0, 5)
set.objfn(lprec, ObjFun)
lp.control(lprec, sense="max")
set.type(lprec,c(1:5),"binary")
#not sure if this can be done easier (how to use the whole matrix MyMatrix?)
for (i in 1:6){
add.constraint(lprec, MyMatrix[i, ], "<=", 1)
}
This gives the following problem:
We can solve this:
solve(lprec)
get.objective(lprec)
MySolution=get.variables(lprec)
And the solution will be (1 0 0 0 1). However, there is also another solution, namely (1,0,1,1,0). We can easily check this (first line checks the constraints, while the second line calculates the value of the objective function):
MyMatrix %*% MySolution
ObjFun %*% MySolution
AnotherSolution=c(1,0,1,1,0)
MyMatrix %*% AnotherSolution
ObjFun %*% AnotherSolution
Is there a way to get all solutions? Right now it only gives MySolution, while I would also like to get AnotherSolution. Is there a way to do that?
I know that in Matlab there is a function lcon2vert which can find vertices in a set defined by constraints - then in a "standard" LP problem one could evaluate the objective function in each vertex (since we know that we will have a corner solution). But here we also have a mixed integer solution.
Alternatively, if this is possible in other packages, I am ok if it is not LPsolveAPI. For example, the same problem can defined as:
library(magrittr)
library(ROI)
library(ompr)
library(ROI.plugin.glpk)
library(ompr.roi)
model <- MIPModel() %>%
add_variable(solution[i], i=1:ncol(MyMatrix), type = 'binary') %>%
set_objective(sum_over(ObjFun[i] * solution[i], i=1:ncol(MyMatrix), sense = 'max')) %>%
add_constraint(sum_over(solution[i] * MyMatrix[j, i], i=1:ncol(MyMatrix)) <= 1, j=1:nrow(MyMatrix))
result <- solve_model(model, with_ROI("glpk", verbose = TRUE))
MyResult=result$solution
Here MyResult is also (1,0,0,0,1). Thank you in advance!