I am trying to solve a large MILP, but it seems like dense problems can be very difficult for moderns solvers. I tried to solve the problem described below considering only constraints (1) and (2) that is the case where x is from the binary set B and continuous set C (so ignoring constraints (3), (4) and (5)). For 1000 variables (800 binary and 200 continuous) in total and 100k scenarios, it took 15 minutes to solve. As I increased scenarios to 200k, the problem was not solved after 7 hours.
I would like to use Benders Decomposition or another method to solve the problem, but I don't know how to formulate the master and subproblems.
The problem is described below:
$\textrm{Minimize }c^{T}x + \eta + \frac{1}{S(1-\alpha )}\sum_{s=1}^{S}[\xi ^{T}x-\eta]^{+}\\\\$
The objective function can be equivalently represented as following by introducing auxiliary variables $z_s\geq0$ along with constraint (2)
$\textrm{Minimize}_{x,y,z_s,\eta} \sum _{i}c_{i}x_{i} + \eta+\frac{1}{S(1-\alpha)}\sum_{s=1}^{S}z_{s}\\ \textrm{Subject to}\\ \begin{align} Ax \geq b \hspace{3cm}(1)\\ z_{s} \geq \sum_{i=1}^{n}\xi_{is}x_{i} - \eta \hspace{1cm} \forall s \hspace{1cm} (2)\\ l_{i}y_{i} \leq x_{i} \leq u_{i}y_{i} \hspace{1cm}\forall i \in N \hspace{0.25cm}and \hspace{0.25cm} i \in\mathbb{Z} \hspace{1cm}(3)\\ \sum_{i}^{}y_{i} \leq 10 \hspace{1cm}(4)\\ \sum_{i}^{}P_{i}x_{i} = 5000000 \hspace{1cm}\forall i \in N \,\,and \,\, i \in\mathbb{Z}\hspace{1cm}(5) \end{align}\\\\$
where
$\hspace{1cm}\bullet$ $x_i\in ${0,1}$ \hspace{0.25cm}\forall i \in B\hspace{0.5cm} \textrm{B is the index set of binary variables}$
$\hspace{1cm}\bullet$ $0\leq x_i\leq1\hspace{0.25cm}\forall i \in C\hspace{0.5cm} \textrm{C is the index set of continuous variables bounded between 0 and 1}$
$\hspace{1cm}\bullet$ $l_i\leq x_i\leq u_i\hspace{0.25cm}\forall i \in N\hspace{0.5cm} \textrm{N is the index set of continuous variables bounded between } l_i$ $\textrm{ and }$ $u_i$
$\hspace{1cm}\bullet$ $l_i\leq x_i\leq u_i\hspace{0.25cm}\forall i \in \mathbb{Z}\hspace{0.5cm} \mathbb{Z}\textrm{ is the index set of integer variables}$
$\hspace{1cm}\bullet$ $y_i\in ${0,1}$ \hspace{0.25cm}\forall i \in N \hspace{0.25cm} and \hspace{0.25cm}\forall i \in \mathbb{Z}\hspace{0.5cm} \textrm{variables y are introduced to enforce bounds on variables x and constraints (3), (4) and (5)}$
$\hspace{1cm}\bullet$ $\textrm{s=1,...S are set of scenarios}$
$\hspace{1cm}\bullet$ $\xi_s$ $\textrm{ are vector of uncertain coefficients or parameters associated with variable }x \textrm{ given for each scenario s}$
$\hspace{1cm}\bullet$ $\alpha = 0.01$
$\hspace{1cm}\bullet$ $P_i$ $\textrm{ are (certain) parameters }$
Could this problem be decomposed using Benders Decomposition, and what will be master and subproblems? If considering all constraints is too much to ask for, I would be happy with answer answer considering just constraints (1) and (2) only for which x is from set B and C (i.e. binary and continuous bounded between 0 and 1). Note: Constraint (1) does not involve any uncertain parameters and is not scenario dependent. It is just like a standard LP/MILP type constraint.
Does decomposition always require creating as many subproblems as the number of scenarios or could we create some subset? For example, the problem may be efficiently solved in extensive form for say 20k scenarios. I am not sure if a decomposition method could be applied if we solve extensive form(s) based on sets of 20k scenarios in each (i.e. breaking the problem in to 10 = 200k/20k large deterministic problems).
Update: I have provided R code for a simplified example to show extensive form of the model in R, considering the continuous variable x case. OMPR package does not seem to allow defining a variable to have different types (such as continuous, binary, etc) based on index.
suppressPackageStartupMessages(library(dplyr, quietly = TRUE))
suppressPackageStartupMessages(library(ROI))
library(ROI.plugin.glpk)
library(ompr)
library(ompr.roi)
extensive <- function(nvar,S,alpha) {
# Generate Random data
p <- matrix(runif(nvar, min=10, max=150), 1, nvar)
xi_temp <- matrix(runif(S*nvar,min=10,max=120), S, nvar)
xi <- sweep(-xi_temp,2,p,FUN="+")
result <- MIPModel() |>
# Variable type x is defined as continuous only here
add_variable(x[i],i=1:nvar, type = "continuous",lb=0, ub=1) |>
add_variable(z[i],i=1:S, type = "continuous",lb=0) |>
add_variable(eta,type = "continuous",lb=0) |>
add_constraint(sum_over(xi[j,i]*x[i], i=1:nvar)<= z[j] +eta, j =1:S) |> # Constraint (2)
add_constraint(sum_over(p[i]*x[i], i=1:nvar) >= 2000) |> # Constraint (1)
set_objective(eta+1/(1-alpha)*1/S*sum_over(z[j], j =1:S)) |>
solve_model(with_ROI(solver = "glpk"))
solution <- get_solution(result,x[i])
}
solution <- extensive(100,1000,0.95)
