How can I write an IF condition for my decision variable for Mixed Integer Linear Programming (MILP) using PuLP GLPK on Python?

Question

I am trying to solve an optimization problem using mixed integer linear programming on PuLP with GLPK solver on Python. So far I have been successful solving basic optimization problems with constraints, such as:

prob = LpProblem("MILP", LpMinimize)
x1 = LpVariable("x1",lowBound=0, cat = 'Binary')
x2 = LpVariable("x2", cat = 'Continuous')
prob += 4*x1 + x2, "Objective Function"
prob += x2 - 4*x1 <= 0
prob += x2 - 2*x1 >= 0
status = prob.solve()
LpStatus[status]
value(x1), value(x2), value(prob.objective)

This gives an optimum result where x1 = 1.0, x2 = 3.0 and Objective Function = 7.0

What I'm trying to figure out is how can I solve an optimization problem with an if condition in, for example, the following constraint:

x1 > 0 IF x2 > 2

or something like:

x1 > 0 IF x2 == 3

Bascially, how can I integrate an if conditional statement into the MILP constraints.

Answer 1

The search term you are looking for is "indicator variable", or "big-M constraint".

As far as I know PULP doesn't directly support indicator variables, so a big-M constraint is the way to go.

A Simple Example: x1 <= 0 IF x2 > 2

from pulp import *

prob = LpProblem("MILP", LpMaximize)
x1 = LpVariable("x1", lowBound=0, upBound=10, cat = 'Continuous')
x2 = LpVariable("x2", lowBound=0, upBound=10, cat = 'Continuous')

prob += 0.5*x1 + x2, "Objective Function"

b1 = LpVariable("b1", cat='Binary')

M1 = 1e6
prob += b1 >= (x1 - 2)/M1

M2 = 1e3
prob += x2 <= M2*(1 - b1)

status = prob.solve()
print(LpStatus[status])
print(x1.varValue, x2.varValue, b1.varValue, pulp.value(prob.objective))

We want a constraint x1 <= 0 to exist when x2 > 2. When x2 <= 2 no such constraint exists (x1 can be either positive or negative).

First we create a binary variable:

b1 = LpVariable("b1", cat='Binary')

Choose this to represent the condition x2 > 2. The easiest way to achieve this adding a constraint:

M1 = 1e6
prob += b1 >= (x2 - 2)/M1

Here M1 is the big-M value. It needs to be chosen such that for the largest possible value of x2 the expression (x2-2)/M is <=1. It should be as small as possible to avoid numerical/scaling issues. Here a value of 10 would work (x2 has upper bound of 10).

To understand how this contraint works, think of the cases, for x2<=2 the right-hand-side is at-most 0, and so has no effect (lower bound of a binary variable already set to 0). However if x2>2 the right hand side will force b1 to be more than 0 - and as a binary variable it will be forced to be 1.

Finally, we need to build the required constraint:

M2 = 1e3
prob += x1 <= M2*(b1 - 1)

Again to understand how this constraint works, consider the cases, if b1 is true (1) the constraint is active and becomes: x1 <= 0. If b1 is false ('0') the constraint becomes x1 <= M2, provided M2 is large enough this will have no effect (here it could be as small as 10 as x1 already has an upper bound of 10.

In the full code above if you vary the coefficient of x1 in the objective function you should notice that b1 is activated/de-activated and the additional constraint applied to x1 as expected.