Express logical constraints : one constraint implies another in linear programming (open source solver)

How can I express the following in linear programming:

2*x + 3*y >= 15 ==> m + n <= 5

x, y, m and n are decision variables (continuous or integer) with relevant upper and lower bounds

I know with CPLEX (OPL and with python API), we can write if-then constraints like :

2*x + 3*y >= 15 => m + n <= 5 (OPL)
model.add_if_then(2*x + 3*y >= 15, m + n <= 5) [python API, docplex]

But how to express above in an open source solver ?

Using indicator constraints, you can model the logic as follows:

2*x + 3*y >= 15 ==> w == 1
w == 1 ==> m + n <= 5

w is a binary variable

indicator are generally (at least in Gurobi) of the form : if binary variable (here w) is true, then that would force the linear constraint to be satisfied. So in Gurobi, you would have to express the following :
2*x + 3*y >= 15 ==> w == 1 with its logical contra-positive i.e.
w == 0 ==> 2*x + 3*y < 15

Strict inequalities aren't allowed in a math programming framework, but we can closely approximate this constraint with 2*x + 3*y <= 15 - epsilon for some small epsilon > 0. Note that we can use epsilon = 1 if all variables are integers.

So finally:

w == 0 ==> 2*x + 3*y <= 15 - epsilon
w == 1 ==> m + n <= 5

In Gurobi :

model.addGenConstrIndicator(w, 0, 2*x + 3*y <= 15 - epsilon)
model.addGenConstrIndicator(w, 1, m + n <= 5)

I am not sure whether indicator constraints are supported in open source solvers.

For open source solvers, however, you can use big M to express following:

2*x + 3*y >= 15 ==> w == 1

implementation:
a*x - (M + e) * w <= b - e
a, x represent coefficients and variables
e (read: epsilon) is some small tolerance value beyond which we will 
   regard the constraint as having been broken. 
   incase variables are integers, we can take e = 1
M (big_M) is the upper bound on 2*x + 3*y - 15
b is the RHS of the constraint

w == 1 ==> m + n <= 5

implementation:
a*x - b <= M * (1 - w)
M (big_M) is the upper bound on : m + n - 5

Below is the implementation in Google OR-Tools, linear solver (python API)

from ortools.linear_solver import pywraplp

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

e = 10e-5 #epsilon

x = s.NumVar(lb = 0, ub = 5, name = "x")
y = s.NumVar(lb = 0, ub = 5, name = "y")

m = s.NumVar(lb = 0, ub = 10, name = "x")
n = s.NumVar(lb = 0, ub = 10, name = "x")

# 2*x + 3*y >= 15 ==> m + n <= 5

# introduce w as a binary variable, then, enforce following
# 2*x + 3*y >= 15 ==> w == 1
# w == 1 ==> m + n <= 5

# to test the above, lets fix x = 3 and y = 4
# inequality : 2*x + 3*y >= 15 now will be true,
# and m + n <= 5 should hold true post the solve
x.SetLb(3)
x.SetUb(3)

y.SetLb(4)
y.SetUb(4)

# add w as a binary variable:
w = s.BoolVar("")

# 2*x + 3*y >= 15 ==> w == 1

# a*x - (M + e) * w <= b - e
# a, x represent coefficients and variables
# e (read: epsilon) is some small tolerance value beyond which we will 
#   regard the constraint as having been broken 
#   is variables are integers, we can take e = 1
# M (big_M) is the upper bound on 2*x + 3*y - 15
# b is the RHS of the constraint

big_M_1 = 2*x.ub() + 3*y.ub() - 15
s.Add(2*x + 3*y - (big_M_1 + e) * w <= 15 - e)

# now we express : w == 1 ==> m + n <= 5

# a*x - b <= M * (1 - w)
# M (big_M) is the upper bound on : m + n - 5

big_M_2 = m.ub() + n.ub() - 5
s.Add(m + n - 5 <= big_M_2 * (1 - w))

s.Maximize(x + y + m + n)

s.Solve()

s.Objective().Value() # objective function value

print("x == " + str(x.SolutionValue())) # x = 3
print("y == " + str(y.SolutionValue())) # y = 4
print("m == " + str(m.SolutionValue())) # m = 5
print("n == " + str(n.SolutionValue())) # n = 0

print("w == " + str(w.SolutionValue())) # w = 1

as you see since we fixed x, y to be 3 and 4 respectively, 2x + 3y >= 15 is satisfied, which forces m + n <= 5 (m = 5 and n = 0 in optimal solution). This is despite of the fact that we are maximising x + y + m + n (without the constraint, m, n would take their upper bound values.

The proper solution is strongly dependent on the possible ranges of x, y, a, b and their integrality. The following is a demonstration only, and you'll need to tune the parameters for whatever it is you're actually doing, in particular the values of M and p for the big-M constraints.

import pulp

m = pulp.LpVariable('m', lowBound=0, upBound=40)
n = pulp.LpVariable('n', lowBound=0, upBound=40)
x = pulp.LpVariable('x', lowBound=2.4)
y = pulp.LpVariable('y', lowBound=3.2)  # Increase lower bound to see constraint changes
b = pulp.LpVariable('b', cat=pulp.LpBinary)
prob = pulp.LpProblem(name='conditional_demo', sense=pulp.LpMinimize)
prob.objective = x + y - m - n

# If 2x + 3y < 15, b == 0.
# If 2x + 3y >= 15, b == 1.
# Epsilon must be used for the first expression, otherwise b is unconstrained at 2x + 3y == 15
# p must be a coefficient sufficiently large that it accounts for the highest possible values of x,y
epsilon = 1e-3
p = 3
prob.addConstraint(b >= ((2*x + 3*y)/15 - 1)/p + epsilon)
prob.addConstraint(b <= (2*x + 3*y)/15)

# If b == 0, then m + n is unconstrained
# If b == 1, then m + n <= 5
# M should be sufficiently large that it will always exceed realistic values of m + n - 5
M = 100
prob.addConstraint(m + n <= 5 + (1-b)*M)

print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal
print('m =', m.value())
print('n =', n.value())
print('x =', x.value())
print('y =', y.value())
print('b =', b.value(), 'because 2x+3y =', 2*x.value() + 3*y.value())

m = 5.0
n = 0.0
x = 2.4
y = 5.2
b = 1.0 because 2x+3y = 20.400000000000002

Note: since constraints on m,n are one-sided, b may also be one-sided; that is, you can loosen it from a binary variable to a singly-bounded integral variable:

import pulp

m = pulp.LpVariable('m', lowBound=0, upBound=40)
n = pulp.LpVariable('n', lowBound=0, upBound=40)
x = pulp.LpVariable('x', lowBound=6)
y = pulp.LpVariable('y', lowBound=6)  # Change lower bound to see constraint changes
b = pulp.LpVariable('b', cat=pulp.LpInteger, upBound=1)
prob = pulp.LpProblem(name='conditional_demo', sense=pulp.LpMinimize)
prob.objective = x + y - m - n

# If 2x + 3y < 15, b == 0.
# If 2x + 3y >= 15, b == 1.
# Epsilon must be used for the first expression, otherwise b is unconstrained at 2x + 3y == 15
# p must be a coefficient sufficiently large that it accounts for the highest possible values of x,y
epsilon = 1e-3
p = 3
prob.addConstraint(b >= ((2*x + 3*y)/15 - 1)/p + epsilon)
prob.addConstraint(b <= (2*x + 3*y)/15)

# If b == 0, then m + n is unconstrained
# If b == 1, then m + n <= 5
# M should be sufficiently large that it will always exceed realistic values of m + n - 5
M = 100
prob.addConstraint(m + n <= 5 + (1-b)*M)

print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal
print('m =', m.value())
print('n =', n.value())
print('x =', x.value())
print('y =', y.value())
print('b =', b.value(), 'because 2x+3y =', 2*x.value() + 3*y.value())