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.

  • 1
    $\begingroup$ Is the uncertainty in the constraint RHSs, the objective function, the constraint matrix, or everywhere? $\endgroup$
    – Max
    Commented Jan 16, 2023 at 17:58
  • $\begingroup$ @Max the uncertainty is in the objective function only. $\endgroup$ Commented Jan 17, 2023 at 5:41

3 Answers 3


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

  • $\begingroup$ This makes perfect sense to me. Although my objective is maximization, I can at least see potential modifications. I would like to learn more about such formulations, do you have a resource that you can guide me to if possible? $\endgroup$ Commented Jan 17, 2023 at 17:32
  • 1
    $\begingroup$ @arvind rathoe My answer shows modification for maximization. $\endgroup$ Commented Jan 17, 2023 at 17:38
  • $\begingroup$ @arvindrathore I got quite a bit out of the Introduction to Stochastic Programming book by Birge and Louveaux, who deal with many of these types of problems. Since you're dealing with probability distributions, this might be a good fit. Additionally, perhaps Robust and Adaptive Optimization by Bertsimas and Den Hertog could be nice to get a very different perspective on 'optimisation under uncertainty'. Between those two books, many methods for dealing with uncertainty should be covered. $\endgroup$
    – Nelewout
    Commented Jan 20, 2023 at 12:50

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.

  • 1
    $\begingroup$ Whoops, didn't see @Nelewout 's answer until after I posted this. I'll leave my answer up for now, as it might be somewhat cleaner and clearer., especially given the availability of indicator functions in Gurobi. $\endgroup$ Commented Jan 17, 2023 at 13:16
  • $\begingroup$ Once I tried a chance constrained example like above and, to me, it worked. $\endgroup$ Commented Jan 17, 2023 at 14:18
  • $\begingroup$ @Sutanu This is not really chance-constrained. The "chance", meaning quantile, is in the objective function. However, my formulation moves the "chance" to the constraints as epigraph formulation for handle quantile in objective. $\endgroup$ Commented Jan 17, 2023 at 15:04
  • $\begingroup$ @MarkL.Stone just for my edification: do you know how Gurobi implements indicator constraints under the hood? I'd imagine that's something big-M-like as well, but would be interested to learn if they do something clever! $\endgroup$
    – Nelewout
    Commented Jan 21, 2023 at 11:31
  • 1
    $\begingroup$ Gurobi does not convert indicator constraints into Big M. You can read the entire thread or.stackexchange.com/questions/231/… for some insight as to how various solvers, including Gurobi, implement indicator constraints, and numerical differences vs. Big M. You may need a Gurobi employee to provide more detailed information. on Gurobi's implementation, and I don't now how forthcoming they would be in providing possibly proprietary details. $\endgroup$ Commented Jan 21, 2023 at 14:52

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

  • $\begingroup$ I think your answer would be more helpful, if you expanded it bit. For example, the magic constant "75" may seem odd without any explanation. $\endgroup$
    – Sune
    Commented Jan 16, 2023 at 16:29
  • $\begingroup$ The OP put an example of 75th percentile or probability 75% which in a normal dist is around 0.67 z-score. $\endgroup$ Commented Jan 16, 2023 at 16:55
  • $\begingroup$ I think I forgot to add an important detail. My data uncertainty is limited to only the objective function. In this case, is there a known way of solving any kind of solver? This is not a novel part of my research so I am hoping to reduce my workload on this. $\endgroup$ Commented Jan 17, 2023 at 5:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.