This article is inspired by the Stack Exchange 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 using indicator constraints.
Because Gurobi's indicator constraints require a binary variable as the indicator variable, we model if \(x > y\) by enforcing \( x > y \rightarrow b=1 \) and \( x \leq y \rightarrow b=0 \). The binary variable \( b \) thus indicates if \( x > y \) is true \((b = 1 )\) or false \((b = 0)\). To model this logic, one can use the following big-\(M\) approach:
x &\geq y + \epsilon - M\cdot (1-b)\\
x &\leq y + M\cdot b\\
where \(\epsilon > 0\) represents a small positive tolerance to simulate a \(>\) constraint. Please note that Gurobi and other solvers do not support strict inequality constraints (\(<\) or \(>\)). 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. Setting \(M\) equal to the largest of the upper bounds of \(x\) and \(y\) (if they are present) often works well.
With the above formulation, we can use variable \(b\) to formulate the indicator constraints
b=1 &\rightarrow z = w_1\\
b=0 &\rightarrow z = w_2.
Python code modeling a conditional statement in Gurobi might look similar to the following:
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")
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")
# Build the rest of the model
Note that the above idea also works if some of the terms \(x\), \(y\), \(z\), \(w_1\), and \(w_2\) are constant values.