Linear programming with conditional constraints in R

Question

I have a linear programming problem where I'm trying to select from a number of binary resources to optimize value, basically a knapsack problem. The issue I'm having is that the different resources have characteristics in common and I want to ensure that my final solution has either 0 or 2 of resources with a specific characteristic. Is there some way to accomplish this? I haven't been able to think of one or find one despite extensive searching. In my data, the decision variables are resources and the constraints are characteristics of those resources. Consider the following code:

library(lpSolve)
const_mat_so<-matrix(c(
  c(0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,1,0,0,1,0,1)
     ,c(0,  0,  0,  0,  0,  0,  0,  1,  0,  1,  0,1,1,0,0,1,1)
     ,c(0,  0,  0,  0,  0,  1,  1,  0,  0,  0,  1,0,1,0,1,0,0)
     ,c(1,  1,  0,  1,  1,  0,  0,  0,  1,  0,  0,0,0,0,0,0,0)
     ,c(8800,   8500,   7600,   8600,   8400,   7500, 7000, 8500,   8800,   7700,   6700,5500,1200,6700,9500,8700,6500)
     ,c(0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,0,0,0,0,0,0)
     ,c(0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,0,0,0,0,0,0)
     ,c(0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,0,0,0,0,0,0)
     ,c(0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  1,0,0,1,0,1,0)
     ,c(0,  0,  1,  0,  0,  0,  0,  0,  1,  0,  0,1,1,0,0,0,0)
     ,c(0,  0,  1,  0,  0,  0,  0,  1,  0,  0,  0,0,0,0,0,0,0)
     ,c(0,  0,  0,  1,  0,  0,  0,  0,  0,  1,  0,1,1,1,0,1,0)
     ,c(0,  1,  0,  0,  0,  0,  0,  1,  0,  0,  0,0,0,0,0,1,0)
     ,c(0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,0,0,0,1,0,0)
     ,c(0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  1,0,0,0,0,0,0)
     ),nrow=15,byrow = TRUE)

const_dir_so<-c("=","=","=","=","<=","<=","<=","<=","<=","<=","<=","<=","<=","<=","<=")

max_cost_so = 25000

objective_so = c(21.0, 19.3, 19.2, 18.8, 18.5, 16.6, 16.4, 16.4, 16.0, 16.0, 14.9, 14.6, 14.0, 13.9,12.0,5.5,24.6)

const_rhs_so<-c(1,1,1,1,25000,3,3,3,2,2,2,2,2,2,2)

x = lp ("max", objective_so, const_mat_so, const_dir_so, const_rhs_so,     all.bin=TRUE, all.int=TRUE 
    )

> x
Success: the objective function is 68.1

> x$solution
 [1] 1 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0

While the above produces a solution, it is not the solution I want because I actually want the last seven constraints to be >=2 or 0. I have no clue how to code this or whether it's possible. Any help would be appreciated. I'm not a linear programming whiz so please forgive any misconceptions regarding the approach.

Answer 1

My understanding is that each of the last 7 constraints are to be greater than 2 or equal to zero, i.e. not 1.

1) There are only 7 such constraints so there are 2^7 = 128 possibilities which is small enough that we can just run every one using the formulation given it he question without excessive runtime and then find the maximum of those.

dec2bin takes a base 10 (i.e. decimal) number and converts it to a binary vector of 0s and 1s. Running it on each number between 0 and 127 gives binary numbers such that the 1s correspond to constraints which are >= 2 (with the rest equal to 0).

dec2bin <- function(dec, digits = 7) {
   # see http://stackoverflow.com/questions/6614283/converting-decimal-to-binary-in-r
   tail(rev(as.integer(intToBits(dec))), digits)
}

runLP <- function(i) {
   bin <- dec2bin(i)
   n <- length(const_rhs_so)  # 15
   ix <- seq(to = n, length = length(bin))  # indexes of last 7 constraints, i.e. 9:15
   const_dir_so[ix] <- ifelse(bin, ">=", "=")
   const_rhs_so[ix] <- 2*bin
   lp("max", objective_so, const_mat_so, const_dir_so, const_rhs_so, all.bin = TRUE)
}

lpout <- lapply(0:127, runLP)
ixmax <- which.max(sapply(lpout, "[[", "objval"))
ans <- lpout[[ixmax]]
ans
ans$solution
tail(c(const_mat_so %*% ans$solution), 7)

giving:

> ans
Success: the objective function is 62 
> ans$solution
 [1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1
> tail(c(const_mat_so %*% ans$solution), 7) # last 7 constraint values
[1] 0 0 0 0 0 0 0

2) In @Erwin Kalvelagen's second alternative it refers to constraining variables but I think what was meant was that x in his answer is the value of the LHS of one of the last 7 constraints. That is, if C is the matrix of the original last 7 constraints then replace those original 7 constraints with these 14 constraints:

 Cx + D1 y <= 0
 Cx + D2 y >= 0

where D1 is a diagonal matrix whose diagonal elements are any sufficiently large negative number and D2 is a diagonal matrix whose diagonal elements are all -2. Here we are optimizing over x and y vectors of binary variables. The x variables are as in the question and there are 7 new y binary variables such that y[i] is 0 to constrain the ith of the last 7 original constraints to 0 or 1 to constrain it to 2 or more. The y variables are called bin in (1). The coefficients of the y variables in the objective are all zero.

In terms of lpSolve R code:

objective_so2 <- c(objective_so, numeric(7))
const_mat_so2 <- cbind(rbind(const_mat_so, const_mat_so[9:15, ]), 
                       rbind(matrix(0, 8, 7), diag(-100, 7), diag(-2, 7)))
const_dir_so2 <- c(const_dir_so, rep(">=", 7))
const_rhs_so2 <- c(const_rhs_so[1:8], numeric(14))
x2 = lp ("max", objective_so2, const_mat_so2, const_dir_so2, const_rhs_so2, all.bin = TRUE)

giving the same value of 62 as in (1). The y variables (last 7) are all 0 which also corresponds to (1). This also provides a double check as two methods have now given consistent answers.

> x2
Success: the objective function is 62 
> x2$solution
 [1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0

Answer 2

I believe LpSolve supports semi-continuous variables. A semi-continuous variable with lower bound L and upper bound U can assume the values 0 or between L and U. I am not sure the R package lpSolve supports this variable type.

However we can simulate this with an extra binary variable y and extra constraints. So you need to make your x variable continuous (or integer if you want only integer values) and add the constraints:

2*y <= x <= U*y

where U is an upper bound for x.

Answer 3

lpSolveAPI package provides more advanced interface to "lp_solve". As @Erwin Kalvelagen mentioned, "lp_solve" and lpSolveAPI supports Semi-continuous variable (Semi-continuous decision variables can take allowed values between their upper and lower bound as well as zero). And a constraint matrix enable you to transfer the outputs of 9-15th constraint formulas into 18-24th variables. For example (about 9th constraint), when x6 + x11 + x14 + x16 - x18 = 0, x6 + x11 + x14 + x16 = x18. So I think you can control x6 + x11 + x14 + x16 via semi-continuous variable, x18.

library(lpSolveAPI)

  ## add 18-24th cols to define the 18-24th variables
const_mat_so2 <- cbind(const_mat_so, rbind(matrix(0, nrow = 8, ncol = 7), diag(-1, 7)))

  ## [EDITED] make a model and set a constraint matrix and objective coefs
model <- make.lp(nrow(const_mat_so2), 0)

for(i in 1:ncol(const_mat_so2)) add.column(model, const_mat_so2[,i])
set.constr.type(model, c(const_dir_so[-c(9:15)], rep("=", 7)))
set.rhs(model, c(const_rhs_so[-c(9:15)], rep(0, 7)))   # each original output - 18-24th = 0

set.objfn(model, c(objective_so, rep(0, 7)))           # 18-24th are 0

  ## define semi-continuous and bounds.
set.semicont(model, col = 18:24)
set.bounds(model, lower = rep(1.9, 7), col = 18:24)   # default upper is Inf.

  ## define other things
set.type(model, col = 1:17, type = "binary")     # original variable
set.type(model, col = 18:24, type = "integer")   # outputs of original constraint formulas
lp.control(model, sense = "max")                 # do maximize

# write.lp(model, "filename.lp", "lp")   # if you want to watch the whole model

solve(model)
get.variables(model)
 # [1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 [18] 0 0 0 0 0 0 0
get.objective(model)
 # [1] 62
t(const_mat_so %*% res[1:17])
 #      [,1] [,2] [,3] [,4]  [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
 # [1,]    1    1    1    1 22300    1    0    0    0     0     0     0     0     0     0