This article is inspired by the stackoverflow post How to write if else statement in Linear programming?

Conditional statements such as `if` \(x > y\) `then` \(z = w_1\) `else` \(z = w_2\), where \(x,y,z,w_1,w_2\) are all optimization variables, can be implemented in Gurobi by introducing an auxiliary binary variable \(b\) and then using Gurobi's indicator constraints.

Because Gurobi's indicator constraints require a binary variable as indicator variable, we model the term `if` \(x > y\) by enforcing that \(b=1\) if \(x > y\) and \(b=0\) if \(x \leq y\). To achieve this, one can use a big-\(M\) approach as follows:

\[\begin{align}

x &\geq y + \epsilon - M\cdot (1-b)\\

x &\leq y + M\cdot b\\

b&\in \{0,1\}

\end{align}\]

where \(\epsilon\) represents a small positive (\(>0\)) tolerance to simulate a \(>\) constraint. Please note that Gurobi does not support strict inequality constraints \(<,>\). In many applications, the \(\epsilon\) tolerance can be dropped, but you should evaluate for your specific use case whether setting \(\epsilon=0\) is feasible. The big-\(M\) value should be chosen as tight as possible to avoid numerical issues. Often setting \(M\) to the largest of the upper bounds of \(x\) and \(y\) (if they are present) works well.

With the above formulation, we can use variable \(b\) to formulate the indicator constraints

\[\begin{align}

b=1 &\rightarrow z = w_1\\

b=0 &\rightarrow z = w_2

\end{align}\]

A Python code modeling a conditional statement in Gurobi might look similar to

import gurobipy as gp

from gurobipy import GRB

# Create a new model

m = gp.Model("test")

# Create variables

x = m.addVar(ub=10, vtype=GRB.CONTINUOUS, name="x")

y = m.addVar(ub=5, vtype=GRB.CONTINUOUS, name="y")

z = m.addVar(vtype=GRB.CONTINUOUS, name="z")

w1 = m.addVar(vtype=GRB.CONTINUOUS, name="w1")

w2 = m.addVar(vtype=GRB.CONTINUOUS, name="w2")

b = m.addVar(vtype=GRB.BINARY, name="b")

# Constants

eps = 0.0001

M = 10 + eps # smallest possible given bounds on x and y

# Model if x>y then b = 1, otherwise b = 0

m.addConstr(x >= y + eps - M * (1 - b), name="bigM_constr1")

m.addConstr(x <= y + M * b, name="bigM_constr2")

# Add indicator constraints

m.addConstr((b == 1) >> (z == w1), name="indicator_constr1")

m.addConstr((b == 0) >> (z == w2), name="indicator_constr2")

# Rest of model

# [...]

Note that the above idea also works if some of the terms \(x,y,z,w_1,w_2\) are constant values.