# How can I solve a linear optimization problem with bounds that are a function of the decision

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

• 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? Oct 18 '21 at 15:15
• 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). Oct 18 '21 at 15:16
• 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. Oct 18 '21 at 15:47
• 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. Oct 19 '21 at 17:53
• @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 ? Oct 20 '21 at 8:03

## 1 Answer

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

• 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. Oct 18 '21 at 15:16
• Yes you can try indicator constraints. Checkout this post : or.stackexchange.com/questions/231/… and Oct 18 '21 at 15:20
• 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. Oct 18 '21 at 15:24
• 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 Oct 18 '21 at 16:05