Pulp Conditional constraint

Question

This is the problem I am trying out to resolve using Pulp :

The factory produces 2 objects x and y and sells them resp 10.5 Dollars and 8.5 Dollars; when the production of object x exceeds 20 units, then a maintenance cost of 10 euros is subtracted for our benefit. How to model this simple problem with solvers?

This is a IF conditional basical constraint. It seems to works using CPLEX :

import cplex
import docplex.mp
from docplex.mp.model import Model

maintenance_cost = 10
maintenance_trigger = 20

# Create model
model = Model(name='LP_example', log_output=True)

# Decisions variables
x = model.integer_var(name='x')
y = model.integer_var(name='y')
z = model.binary_var(name='z')

# Objective function
model.maximize(10.5 * x + 8.5 * y - z * maintenance_cost)

# Constraints
model.add_constraint(3 * x + 2 * y  <= 420)

# if then constraint
model.add_constraint(model.if_then(x >= maintenance_trigger, z == 1))


model.print_information() 
sol_model = model.solve()
model.print_solution()

but, I am bloqued, when using Pulp, how should I set the z variable correctly , how should I set the constraint ? I have no idea :

import pulp as p 
Lp_prob = p.LpProblem('Problem', p.LpMaximize)  

x = p.LpVariable("x", lowBound = 0, cat='Integer')   # Create a variable x >= 0 
y = p.LpVariable("y", lowBound = 0, cat='Integer')   # Create a variable y >= 0 
z = p.LpVariable("z", lowBound=0, cat='Binary')

Lp_prob +=  10.5 * x + 8.5 * y - maintenance_cost * z

Lp_prob += 3 * x + 2 * y  <= 420
 
# if then constraint
# The problematic part (Conditional constraint):
Lp_prob += z  == 1 # ???

status = Lp_prob.solve()   
print(p.value(x)," x produced objects")
print(p.value(y) , "y produced objects"  )
print(p.value(z) , "The value of the binary variable, used or not 1 or 0"  )
print(p.value(Lp_prob.objective) ,"our Profit" )

EDIT : I have found a solution there : How can I write an IF condition for my decision variable for Mixed Integer Linear Programming (MILP) using PuLP GLPK on Python?

So my pulp code now includes this "BigM" conditional constraint:

M1 = 1e6
Lp_prob += z >= (x - 20)/M1

I do not understand how M1 works, but the result seems to be ok now, any info appreciated.

Edit: I have removed both MIN production volume constraints, because my wish was only to understand if my BigM constraints was well written using Pulp. I don't know how to write "When the production of object x exceeds 20 units, then a maintenance cost of 10 dollars is subtracted for our benefit" and I'd like to add a second big M like this :"When the production of object x exceeds 50 units, then a maintenance cost of 20 dollars is subtracted for our benefit, not 10 dollars".

Answer 1

1e6 is overkill. Twice the maximum possible value of x is enough. You also need both lower and upper constraints using your binary variable. Finally, as in the comments, your binary variable is meaningless and will always be 1 because 100 > 20.

import pulp

maintenance_cost = 10
maintenance_trigger = 20

product_x = pulp.LpVariable(name='x', cat=pulp.LpInteger, lowBound=100)
product_y = pulp.LpVariable(name='y', cat=pulp.LpInteger, lowBound=40)
x_is_bulk = pulp.LpVariable(name='z', cat=pulp.LpBinary)

model = pulp.LpProblem(name='LP_example', sense=pulp.LpMaximize)
model.objective = 10.5*product_x + 8.5*product_y - x_is_bulk*maintenance_cost
model.addConstraint(3*product_x + 2*product_y <= 420)

# if x is bulk, x production must be at least 20
model.addConstraint(product_x >= maintenance_trigger*x_is_bulk)

# if x is not bulk, x production must be less than 20
x_max = (420 - 2*40)/3
model.addConstraint(product_x <= maintenance_trigger - 0.5 + x_is_bulk*2*x_max)

print(model)
model.solve()
assert model.status == pulp.LpStatusOptimal
print('x production:', product_x.value())
print('y production:', product_y.value())
print('x is bulk:', x_is_bulk.value() > 0.5)

LP_example:
MAXIMIZE
10.5*x + 8.5*y + -10*z + 0.0
SUBJECT TO
_C1: 3 x + 2 y <= 420

_C2: x - 20 z >= 0

_C3: x - 226.666666667 z <= 19.5

VARIABLES
100 <= x Integer
40 <= y Integer
0 <= z <= 1 Integer

Result - Optimal solution found

Objective value:                1550.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.01

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

x production: 100.0
y production: 60.0
x is bulk: True

Multi-step Maintenance Demonstration

I demonstrate how to do this with one integer covering the three maintenance (bulk) classes and one binary variable restricting the last (peak) class:

import pulp

x_min = 10
y_min = 40
x_coef = 3
y_coef = 2
production_max = 420
product_x = pulp.LpVariable(name='product_x', cat=pulp.LpInteger, lowBound=x_min - 0.5)
product_y = pulp.LpVariable(name='product_y', cat=pulp.LpInteger, lowBound=y_min - 0.5)
x_bulk_class = pulp.LpVariable(name='x_bulk_class', cat=pulp.LpInteger, lowBound=0, upBound=2)
x_peak_bulk = pulp.LpVariable(name='x_peak_bulk', cat=pulp.LpBinary)

maint_cost_increment = 10
maint_trigger_0 = 20
maint_trigger_1 = 50
maint_penalty = x_bulk_class*maint_cost_increment

model = pulp.LpProblem(name='LP_example', sense=pulp.LpMaximize)
model.objective = 12.95*product_x + 8.5*product_y - maint_penalty
model.addConstraint(x_coef*product_x + y_coef*product_y <= production_max, name='production')

# x production must always be more than the (non-peak-limited) penalty step function
# this is a step function with three integral values
maint_lobound = (maint_trigger_1 - maint_trigger_0)*(x_bulk_class - 1) + maint_trigger_0
model.addConstraint(product_x >= maint_lobound - 0.5, name='maint_lo')

# x production must be less than the upper penalty step function, until the peak
maint_upbound = maint_lobound + maint_trigger_1 - maint_trigger_0
x_max = (production_max - y_coef*y_min)/x_coef
relief = x_peak_bulk*2*x_max
model.addConstraint(product_x <= maint_upbound - 0.5 + relief, name='maint_up')

# the bulk peak segment occurs when bulk class >= 2
model.addConstraint(2*x_peak_bulk <= x_bulk_class, name='peak_lo')
model.addConstraint(1+x_peak_bulk >= x_bulk_class, name='peak_up')

print(model)
model.solve()
assert model.status == pulp.LpStatusOptimal
print('x production:', product_x.value())
print('y production:', product_y.value())
print('x bulk class:', x_bulk_class.value())
print('x peak bulk:', x_peak_bulk.value() > 0.5)

LP_example:
MAXIMIZE
12.95*product_x + 8.5*product_y + -10*x_bulk_class + 0.0
SUBJECT TO
production: 3 product_x + 2 product_y <= 420

maint_lo: product_x - 30 x_bulk_class >= -10.5

maint_up: product_x - 30 x_bulk_class - 226.666666667 x_peak_bulk <= 19.5

peak_lo: - x_bulk_class + 2 x_peak_bulk <= 0

peak_up: - x_bulk_class + x_peak_bulk >= -1

VARIABLES
9.5 <= product_x Integer
39.5 <= product_y Integer
0 <= x_bulk_class <= 2 Integer
0 <= x_peak_bulk <= 1 Integer


Result - Optimal solution found

Objective value:                1788.60000000
Enumerated nodes:               0
Total iterations:               2
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.01

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

x production: 18.0
y production: 183.0
x bulk class: 0.0
x peak bulk: False