lpSolve issue with no feasible solution

Question

I'm trying to select 5 values in the first row, 1 value in rows 2-9, and maximize the objective function. The only thing I can't figure out how to embed in the problem is that I can only choose a value once - which is to say, if value 49.27 is selected for row 1, I can't also choose value 49.27 for row 9.

obj <- c(
  49.2723892,
  31.4629786,
  25.3265757,
  23.3915431,
  19.5677222,
  19.2621303,
  18.2912147,
  16.7899278,
  14.5498144,
  12.7440869,
  11.6331623,
  11.0012557,
  7.9473880,
  7.8453557,
  3.7182273,
  0.2084929,
  0.0100000
)

mat <- as.matrix(c(1, 0, 0, 0, 0, 0, 0, 0, 1)) %>%
  cbind( c(0, 1, 1, 0, 0, 0, 0, 0, 1),
         c(1, 0, 0, 0, 0, 0, 0, 0, 1),
         c(0, 0, 1, 1, 0, 0, 0, 0, 1),
         c(0, 0, 0, 0, 1, 1, 0, 0, 1),
         c(0, 0, 0, 0, 0, 1, 1, 0, 1),
         c(0, 0, 0, 0, 0, 0, 1, 1, 1),
         c(0, 0, 0, 0, 1, 1, 0, 0, 1),
         c(1, 0, 0, 0, 0, 0, 0, 0, 1),
         c(1, 0, 0, 0, 0, 0, 0, 0, 1),
         c(0, 1, 0, 1, 0, 0, 0, 0, 1),
         c(0, 1, 1, 0, 0, 0, 0, 0, 1),
         c(1, 0, 0, 0, 0, 0, 0, 0, 1),
         c(1, 0, 0, 0, 0, 0, 0, 0, 1),
         c(1, 0, 0, 0, 0, 0, 0, 0, 1),
         c(0, 0, 0, 0, 0, 0, 0, 1, 1),
         c(0, 0, 0, 0, 0, 0, 0, 1, 1)
  )

const.dir <- c("==", "==", "==", "==","==","==","==", "==", "==")
const.rhs <- c(5, 1, 1, 1, 1, 1, 1, 1, 1)

sol <- lp(direction = "max", obj,
          mat, const.dir, const.rhs)

Am I missing something?

Why is it returning the following:

"Error: no feasible solution found"

Answer 1

Let the decision variables be a 0-1 matrix x with the same dimensions as mat and reformulate as:

max sum(x %*% obj)
sumRows(mat * x) == c(5, 1, ..., 1)
sumRows(x) == c(5, 1, ..., 1)
sumCols(x) <= 1

It is a bit easier to specify this with the CVXR package. Note that sum_entries(..., axis = 1) is the row sums and sum_entries(..., axis = 2) is the column sums when using this package.

library(CVXR)

nr <- nrow(mat)
nc <- ncol(mat)

x <- Variable(nr, nc, boolean = TRUE)
objective <- Maximize(sum(x %*% obj))
constraints <- list(
  sum_entries(mat * x, axis = 1) == const.rhs, 
  sum_entries(x, axis = 1) == const.rhs, 
  sum_entries(x, axis = 2) <= 1L)
prob <- Problem(objective, constraints)
res <- solve(prob)
 
res$status
## [1] "optimal"
 
res$value
## [1] 261.2402

r <- round(res$getValue(x))
res$getValue(x)
## ...snip...

colSums(row(r) * r) # for each col show which row is 1
##  [1] 1 2 1 4 5 7 8 6 1 1 9 3 1 0 0 0 0

Note

Inputs used

const.rhs <- c(5, 1, 1, 1, 1, 1, 1, 1, 1)
mat <- structure(c(1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 
1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 
0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 
0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 
1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 
1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 
0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 
1, 0, 0, 0, 0, 0, 0, 0, 1, 1), dim = c(9L, 17L))