1

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.

josliber
  • 43,891
  • 12
  • 98
  • 133
rrbest
  • 1,619
  • 3
  • 14
  • 22

1 Answers1

3

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

josliber
  • 43,891
  • 12
  • 98
  • 133