I am interested in building what is essentially an or statement in a linear program. Currently I am using lpSolveAPI in R, but I also hope to get answers on the best practices in linear programming for this type of problem.
Here is some working code using the lpSolveAPI package:
library(lpSolveAPI)
# Defines the 'score' which we seek to maximize
mat <- data.frame(score = c(0.04, 0.04, -.11, -.03, 0.04, 0.04, -0.06, 0.01))
# Side specifies whether something is a 'sell' or a 'buy'
mat[, 'side'] <- c(1, 1, -1, -1, 1, 1, -1, -1)
# 'a' and 'b' are descriptive factors defining what something is thing a or thing b.
mat[, 'a'] <- c(1, 1, 1, 1, 0, 0, 0, 0)
mat[, 'b'] <- c(0, 0, 0, 0, 1, 1, 1, 1)
# the lower bound for all sides = -1 have a lower bound of some negative value and an upper bound of zero
mat[, "lb"] <- c(0, 0, -0.2, -0.4, 0, 0, -.1, -.2)
# the upper bound for all sides = 1 have a lower bound of zero and an upper bound of 1
mat[, 'ub'] <- c(1, 1, 0, 0, 1, 1, 0, 0)
# this is just initializing our variables field which will be populated later
mat[, 'x'] <- rep(0, 8)
# using the lpSolveAPI package, create a program with the number of rows in the matrix 'mat'
LP <- make.lp(0, nrow(mat))
set.objfn(LP, -mat[, 'score'])
set.bounds(LP, lower = mat[, 'lb'])
set.bounds(LP, upper = mat[, 'ub'])
# This constraint requires that the sum of all of x must be equal to zero. In other words, the amount of sells equals the amount of buys
add.constraint(LP, rep(1, nrow(mat)), "=", 0)
solve(LP)
get.objective(LP)
get.variables(LP)
mat$x <- get.variables(LP)
When you run this code and look at the results you'll see that x6 = 0.9, x7 = -0.1, and x8 = -0.2. The interpretation of this is that 90% of b is purchased and 30% of b is sold.
I am looking to build a constraint that requires that if a or b is sold in any spot, it cannot also be bought. (and vice versa)
In this particular case, I would expect that the optimal solution would be that 20% of a is sold (x3) and 20% of b is bought (x5 or x6). x = (0, 0, -0.2, 0, 0.2, 0, 0, 0)
Said another way, if you choose to have a value other than zero in x1 and/or x2, you must have 0 in x3 and x4. Similarly, if you choose to have a value other than zero in x3 and/or x4, you must have a value of zero in x1 and x2.
In general, constraints of the form "all of these variables are either non-negative or non-positive" require the introduction of a binary variable that indicates if the positive or negative case has been selected.
Let's take the case of the first four variables; we will introduce a binary variable p1 that indicates whether were are in the non-negative regime (p1=1) or the non-positive regime (p1=0). Right now we have the following constraints:
0 <= x1 <= 1
0 <= x2 <= 1
-0.2 <= x3 <= 0
-0.4 <= x4 <= 0
If p1=1, then the lower bounds must be at least 0, and if p1=0 then the upper bounds must be no more than 0. We can rewrite our constraints as:
0 <= x1 <= 1 * p1
0 <= x2 <= 1 * p1
-0.2 * (1-p1) <= x3 <= 0
-0.4 * (1-p1) <= x4 <= 0
Basically, we need to multiply all the positive upper bounds by p1, which has no impact if we're in the non-negative space and sets the upper bound to 0 if we're in the non-positive case. Similarly, we need to multiple all the negative lower bounds by (1-p1), which has no impact if we're in the non-positive space and sets the lower bound to 0 if we're in the non-negative case.
We similarly define a new binary variable p2 for the second set of variables and update our constraints as follows:
0 <= x5 <= 1 * p2
0 <= x6 <= 1 * p2
-0.1 * (1-p2) <= x7 <= 0
-0.2 * (1-p2) <= x8 <= 0
In code, this would look like:
n <- nrow(mat)
LP <- make.lp(0, n+2)
set.objfn(LP, c(-mat[, 'score'], 0, 0))
set.type(LP, n+1, "binary")
set.type(LP, n+2, "binary")
set.bounds(LP, lower = c(mat[, 'lb'], 0, 0))
set.bounds(LP, upper = c(mat[, 'ub'], 1, 1))
for (i in 1:n) {
add.constraint(LP, c(i == (1:n), (mat[i,'a'] == 1) * mat[i,'lb'], (mat[i,'b'] == 1) * mat[i,'lb']), ">=", mat[i,'lb'])
add.constraint(LP, c(i == (1:n), -(mat[i,'a'] == 1) * mat[i,'ub'], -(mat[i,'b'] == 1) * mat[i,'ub']), "<=", 0)
}
add.constraint(LP, c(rep(1, n), 0, 0), "=", 0)
solve(LP)
get.objective(LP)
# [1] -0.058
round(head(get.variables(LP), -2), 3)
# [1] 0.0 0.0 -0.2 -0.4 0.0 0.6 0.0 0.0
The optimal solution found, with objective value 0.058, actually outperforms the one found by hand by the OP (which only has objective value 0.03).