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

Question

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

Answer 1

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