6
$\begingroup$

I am using the rDEA package's linear optimization function multi_glpk_solve_LP(which is built on top of the Rglpk package's linear optimization function Rglpk_solve_LP) to optimise the following problem:

## minimize: 20 x_1 + 34 x_2 + 22 x_3 + 27 x_3
## subject to:
## x_1 + x_2 >= 0
## x_1 + x_2 <= 4650
## x_1 + x_2 + x_3 >= 0
## x_1 + x_2 + x_3 <= 4650
## x_1 + x_2 + x_3 + x_4 == 0
## x_1, x_2, x_3, x_4 are real numbers

## Bounds
## 0 <= x_1 <= 4650
## -4650 <= x_2 <= 4650
## -4650 <= x_3 <= 4650
## -4650 <= x_4 <= 4650

Until now there is no issue as I wrote the problem as follows and the results are correct:

max_wg <- 18600/4
max <- F
prices <- c(20, 34, 22, 27)
d_const <- matrix(c(1,1,1,1,1,1,0,1,1,0,0,1), nrow = length(prices) - 1)
v <- c(rep(2, nrow(d_const)-1),1)
d_const <- d_const[rep(1:nrow(d_const), times = v), ]

dir <- c(rep(c(">=", "<="), (nrow(d_const)-1)/2), "==")
rhs <- c(rep(c(0, max_wg), (nrow(d_const)-1)/2), 0)
bounds <- list(lower = list(ind = c(1L, 2L, 3L, 4L), val = c(0, -max_wg, -max_wg, -max_wg )),
               upper = list(ind = c(1L, 2L, 3L, 4L), val = c(max_wg, max_wg, max_wg, max_wg)))

multi_glpk_solve_LP(obj = prices, mat = d_const, dir = dir, rhs = rhs, bounds = bounds, max = max)

However, I would like to optimize the same problem now, but where bounds are a function of the sum of x_i so I should have :

## minimize: 20 x_1 + 34 x_2 + 22 x_3 + 27 x_3
## subject to:
## x_1 + x_2 >= 0
## x_1 + x_2 <= 4650
## x_1 + x_2 + x_3 >= 0
## x_1 + x_2 + x_3 <= 4650
## x_1 + x_2 + x_3 + x_4 == 0
## x_1, x_2, x_3, x_4 are real numbers

## Bounds
## 0 <= x_1 <= 4650

## -(60% + 1/4 * x_1/max_wg) * 4650 <= x_2 <= 4650 if x_1 <= max_wg * 60%
## -(60% + 1/4 * x_1/max_wg) * 4650 <= x_2 <= (3/2 * (100 - x_1/max_wg) + 40%) * 4650 if x_1 > max_wg * 60% and x_1 <=  max_wg * 80%
## -((20% * x_1/max_wg ) + 90%) * 4650 <= x_2 <= (3/2 * (100 - x_1/max_wg) + 40%) * 4650 if x_1 >  max_wg * 80%

## -(60% + 1/4 * (x_1+x_2)/max_wg ) * 4650 <= x_3 <= 4650 if (x_1+x_2) <= max_wg * 60%
## -(60% + 1/4 * (x_1+x_2)/max_wg ) * 4650 <= x_3 <= (3/2 * (100 - (x_1+x_2)/max_wg ) + 40%) * 4650 if (x_1+x_2)> max_wg * 60% and (x_1+x_2)<=  max_wg * 80%
## -((20% * (x_1+x_2)/max_wg ) + 90%) * 4650 <= x_3 <= (3/2 * (100 - (x_1+x_2)/max_wg ) + 40%) * 4650 if (x_1+x_2)>  max_wg * 80%

## -(60% + 1/4 * (x_1+x_2+x_3)/max_wg ) * 4650 <= x_4 <= 4650 if (x_1+x_2+x_3) <= max_wg * 60%
## -(60% + 1/4 * (x_1+x_2+x_3)/max_wg ) * 4650 <= x_4 <= (3/2 * (100 - (x_1+x_2+x_3)/max_wg ) + 40%) * 4650 if (x_1+x_2+x_3)> max_wg * 60% and (x_1+x_2+x_3)<=  max_wg * 80%
## -((20% * (x_1+x_2+x_3)/max_wg ) + 90%) * 4650 <= x_4 <= (3/2 * (100 - (x_1+x_2+x_3)/max_wg ) + 40%) * 4650 if (x_1+x_2+x_3)>  max_wg * 80%

If there is no package in R that can do this I could eventually do it in Python. Basically the bounds are a function of the cumulative sums of the values of the optimization.

EDIT: Made the bounds condition clearer. The answer that I was given does not explain how to make the y function a function of xi

$\endgroup$
5
  • $\begingroup$ I migrated this from Quant.SE, I would appreciate some feedback whether regular users here feel this question fits this site well. I believe it’s operations research but is this considered too basic or too easy or something else? $\endgroup$
    – Bob Jansen
    Oct 18 at 15:15
  • 1
    $\begingroup$ This is a valid question here. It would have been better with some ideas on how you think you could do it, and it would be even better if it were language agnostic (work with a mathematical model instead directly, as this is a pure modeling question). $\endgroup$
    – Kuifje
    Oct 18 at 15:16
  • 1
    $\begingroup$ The contingent bounds stated in the second part of the question don't make sense. For instance, the second bound say "... if x_1 <= 50% * 4650*4/2", but that expression reduces to 4650, which is the stated upper bound for x_1 ... so the "if" is always true. $\endgroup$
    – prubin
    Oct 18 at 15:47
  • 1
    $\begingroup$ Just to explicitly say what Kuifje's answer implicitly assumes: A "bound" that depends on other variables is not a bound at all, it is a constraint. $\endgroup$
    – Ben Voigt
    Oct 19 at 17:53
  • $\begingroup$ @BenVoigt This means that I need to add them as consttraints in the matrix to optimize ok thank you. But then do I have to specify bounds on something other than x1 ? $\endgroup$ Oct 20 at 8:03
6
$\begingroup$

You want to model constraints of the form $$ x_1 \le a_1 \quad \Longrightarrow \quad x_2 \in [b_2,c_2] $$ Define a binary variable $y \in \{0,1\}$ and use constraints \begin{align*} 0 &\le x_1 \le a_1 + M_1(1-y) \tag{1}\\ b_2 - M_3(1-y) &\le x_2 \le c_2 +M_2 (1-y) \tag{2} \end{align*}

Constraint $(1)$ enforces $y=1 \; \Longrightarrow \; x_1 \le a_1$, constraint $(2)$ enforces $y=1 \; \Longrightarrow \; x_1 \in [b_2,c_2]$. $M_1,M_2,M_3$ are large constants, and as @Mark L. Stone mentions in the comment section, these constants should be the smallest largest constants you can think of, that is, they should be "big enough". For example, you could use $M_1 := L_1 -a_1$, where $L_1$ is your upper bound on $x_1$.

Define binary variables for the other intervals, adapt the equations, and make sure only one interval is used per variable with $\sum_i y_i = 1 $.

$\endgroup$
4
  • 1
    $\begingroup$ I don't know this R package that is being used in the question, but if it supports indicator constraints then those might be better to use than big-M-constraints. $\endgroup$ Oct 18 at 15:16
  • $\begingroup$ Yes you can try indicator constraints. Checkout this post : or.stackexchange.com/questions/231/… and $\endgroup$
    – Kuifje
    Oct 18 at 15:20
  • 2
    $\begingroup$ I suggest you give more guidance on the Big M values, otherwise "disaster" may ensue.. As Johan Lofberg says, it should be called just big enough M. $\endgroup$ Oct 18 at 15:24
  • $\begingroup$ Ok I get it so y = 1 would mean that the condition is met so sum(xi) > the threshhold and so the boudaries have to be adjusted. But in your answer yi is not a function of the xj with j < i $\endgroup$ Oct 18 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.