# 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, 2021 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, 2021 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.
– prubin
Oct 18, 2021 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, 2021 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, 2021 at 8:03

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