roptimizationlpsolve

lpSolve issue with no feasible solution


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"

Solution

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

    max sum(x * m)
    sumRows(m*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 = 2) <= 1L)
    prob <- Problem(objective, constraints)
    res <- solve(prob)
    
    res$status
    ## [1] "optimal"
    
    res$value
    ## ...snip...
    
    r <- round(res$getValue(x))
    colSums(row(r) * r) # for each col show which row is 1
    ## [1] 1 9 3 4 5 7 8 6 1 1 2 3 1 1 5 4 1
    
    

    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))