Operations Research Stack Exchange is a question and answer site for operations research and analytics professionals, educators, and students. It only takes a minute to sign up.
Sign up to join this community
I have a research problem where my Mixed Integer Linear Program has data that follow probability distributions. I am approaching this by creating some m instances through realizations of these random variables (data).
How do I approach this problem in Gurobi?
For a better explanation, my data will be simple numerical data once the random variables are realized. But I will have created 1000 such instances and I want my model to perform optimally on say the 75th quantile case.
One way to do this is to rewrite the objective somewhat.
I'm going to start from an objective of the following form: $$ \min_x c(\zeta)^\top x, $$ with $c(\zeta)$ a cost vector depending on a random vector $\zeta$ with known probability distribution, and $x$ a decision vector, both of suitable dimension.
Now, the problem using $\zeta$ directly is that the support of its probability distribution could be very large. Instead, we can sample a representative set of $n \in \mathbb{N}$ realisations $\zeta_1, \zeta_2, \ldots, \zeta_n$ from the probability distribution of $\zeta$. If these samples adequately capture the randomness in $\zeta$, the resulting sampling solution should also be close to the optimal solution had we used $\zeta$ directly.
We can restate the sampling objective using an epigraph formulation, as $$ \min_{x, t} t $$ with $t \in \mathbb{R}$, subject to $$ t \ge c(\zeta_i)^\top x, \qquad i = 1, 2, \ldots, n. $$
Using binary variables $z_i \in \{0, 1\}$ for $i = 1,\ldots, n$, we can rewrite this to add the condition that a fraction $\alpha \in [0,1]$ of these constraints may be violated: \begin{align} \textstyle \sum_{i = 1}^n z_i &\le \alpha n, \\ t &\ge c(\zeta_i)^\top x - Mz_i \qquad i = 1, 2, \ldots, n, \end{align} for some $M$ big enough. A concrete value for $M$ could be the objective value you get when optimising with $M = \alpha = 0$.
For the $0.75$ quantile, set $\alpha = 0.25$.
The following formulation is more or less the same as the formulation used in How to represent a constraint on the kth-smallest function?. So thanks to @RobPratt for providing a comment which improved that formulation, and hence this one.
Le $N$ be the number off scenarios. In your case, 1000.
Let $q$ be the applicable quantile. In your case, 0.75
Let $f(i)$ be the non-robustified objective value for scenario $i$.
Declare an optimization variable $t$, which will be the objective function to be minimized in the robustified problem.
Declare $b$ as a vector of $N$ binary optimization variables.
Add the indicator constraints:
$b(i) = 1 \rightarrow f(i) \le t$ for $i=1,...,N$
(If maximizing instead of minimizing, instead use $b(i) = 1 \rightarrow f(i) \ge t$ for $i=1,...,N$)
Add the constraint $\Sigma_{i=1}^N b(i) \ge q N$
If this were implemented in a modeling system (and solver other than Gurobi) which did not have indicator constraints, the Big M equivalent could be used instead.
Robust optimization ($+/-$) is still a research area and I am not sure if Gurobi has an automated way unless you want to follow this example from Wolfram. In case it's stochastic programming one way is to apply sample approximation average (SAA) if continuous distribution and parameters are known. You can generate sample data at cdf $\ge 75$% then take expected value of the sample. You can repeat this process itself. If you are using chance constraints then you may need to determine the z-score at $75$% and then use appropriate relation.
You can check LINDO example and GAMS example
