Benjamin Dummitt - 6 months ago 43
R Question

Linear programming with conditional constraints in R

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.

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