CODEWITHSUNDEEP
linear-programmingmixed-integer-programming
bobby
bobby

Reputation: 113

reification of a constraint in linear programming using open source solvers

How can we reify a constraint (and enforce it) in linear programming:

delta == 1 ==> a*x <= b
delta == 1 ==> a*x >= b
delta == 1 ==> a*x == b

where delta is a boolean variable

I am aware that through commercial solvers like gurobi, cplex we have already a built in capability of indicator constraints, but while using an open source linear solver how can we express constraint reification.

Upvotes: 2

Views: 260

Answers (3)

Sebastian Wozny
Sebastian Wozny

Reputation: 17526

With the python library PuLP which ships bundling the open source Coin-OR suite we can realize it as a big-M constraint:

import pulp
problem = pulp.LpProblem("BigMConstraintProblem", pulp.LpMinimize)

a = 2
b = 5
M = 100  # Big M value

x = pulp.LpVariable("x", lowBound=0, cat='Continuous')
delta = pulp.LpVariable("delta", lowBound=0, upBound=1, cat='Binary')

problem += a * x

# Add the Big M constraint
problem += a * x <= b + (1 - delta) * M

problem.solve()

PuLP ships with both a linear and integer solver.

To install PuLP:

python -m pip install pulp

Upvotes: 0

Erwin Kalvelagen
Erwin Kalvelagen

Reputation: 16782

(1) Using big-M constraints. E.g. for the first:

     a*x <= b + (1-δ)M

where M is a large enough constant. This needs to be chosen with some care: the smaller, the better (but it should be large enough to make the constraint non-binding when δ=0).

(2) Using SOS1 variables. This does not require a big-M constant. Not all open-source solvers/modeling tools support SOS1 variables.

(3) Using indicator constraints. Again no big-M constants are needed. Often not available in open-source solvers and modeling tools. It would be good if there is more support for this as it is a very useful modeling concept.

Upvotes: 1

Bhartendu Awasthi
Bhartendu Awasthi

Reputation: 1029

Would advise you to read Model building in Mathematical programming by HP Williams, especially chapter 9

I have used Google-ortool's linear solver (python API) to express your indicator constraints. You can install ortools via pip install ortools. Solvers such cbc, SCIP, etc gets shipped with ortools

from ortools.linear_solver import pywraplp

# CASE 1 : 
    
# delta == 1 ==> a*x <= b
# concrete example : 
# delta == 1 ==> 3x + 2y <= 8

solver = pywraplp.Solver("", pywraplp.Solver.SCIP_MIXED_INTEGER_PROGRAMMING)

delta = solver.BoolVar("") 
x = solver.IntVar(0, 10, name= "")
y = solver.IntVar(0, 10, name= "")

big_M = 8  # upper bound for the expression a*x - b

# a*x - b <= big_M * (1 - delta)
solver.Add(3*x + 2*y - 8 <= big_M * (1-delta))    

# check : fix delta == 1
delta.SetBounds(lb = 1, ub = 1)
solver.Maximize(3*x + 2*y)

sol = solver.Solve()

solver.Objective().Value() # results in 8
x.SolutionValue() # x = 2
y.SolutionValue() # y = 1
delta.SolutionValue() # delta = 1


# CASE 2 :

# delta == 1 ==> a*x <= b
# concrete example : 
# delta == 1 => 3*x + 2*y >= 8

solver = pywraplp.Solver("", pywraplp.Solver.SCIP_MIXED_INTEGER_PROGRAMMING)

delta = solver.BoolVar("") 
x = solver.IntVar(0, 10, name= "")
y = solver.IntVar(0, 10, name= "")

small_M = -8 # lower bound for the expression a*x - b

# a*x >= small_M * (1 - delta) + b
solver.Add(3*x + 2*y >= -small_M * (1 - delta) + 8)

# check : fix delta == 1
delta.SetBounds(lb = 1, ub = 1)
solver.Minimize(3*x + 2*y)

sol = solver.Solve()

solver.Objective().Value() # results in 8
x.SolutionValue() # x = 0
y.SolutionValue() # y = 4
delta.SolutionValue() # delta = 1


# CASE 3 :

# delta == 1 ==> a*x <= b
# concrete example : 
# delta == 1 => 3*x + 2*y == 13

solver = pywraplp.Solver("", pywraplp.Solver.SCIP_MIXED_INTEGER_PROGRAMMING)

delta = solver.BoolVar("") 
x = solver.IntVar(0, 10, name= "")
y = solver.IntVar(0, 10, name= "")

big_M = 13  # upper bound for the expression a*x - b
small_M = -13 # lower bound for the expression a*x - b

# a*x - b <= big_M * (1 - delta)
# a*x >= small_M * (1 - delta) + b

solver.Add(3*x + 2*y - 13 <= big_M * (1-delta))
solver.Add(3*x + 2*y >= -small_M * (1 - delta) + 8)

# check : fix delta == 1
delta.SetBounds(lb = 1, ub = 1)
solver.Maximize(3*x + 2*y)

sol = solver.Solve()

solver.Objective().Value() # results in 8
x.SolutionValue() # x = 0
y.SolutionValue() # y = 4
delta.SolutionValue() # delta = 1

Upvotes: 2

Related Questions