How to write a conditional constraint in PuLP

Question

Background

I'm setting up a model to optimize dispatch of a solar + battery + electrolyser/fuel-cell system using PuLP. The idea is to maximize revenue by charging the Battery and producing Hydrogen from solar generation when prices are low, then discharging the Battery and burning Hydrogen when prices are high.

research & troubleshooting

I'm new to PuLP, I've spent a while working on this, using online guides and SO for assistance but my code is still returning rubbish answers. I've spent a while browsing SO searching for answers (here, here, here and here). and I think I know what the issue might be.

Potential Problem Identification

I think my main problem is that I have a bunch of logical constraints which I need to write in a different form to work in PuLP. I currently have no indicator variables or Big M constraints which I believe are useful way of writing conditional constraints. I tried to write them myself, but with no luck because I didn't know how to write them in PuLP (which I believe doesn't allow you to use IF or np.where, etc)

As for Big M constraints, I've had a go at writing them in what I believe to be an appropriate way, but it returns the error "Non-constant expressions cannot be multiplied".

Problem Example

I've attached sample of code below of what I am trying to achieve via a conditional constraint in PuLP.

In my system the Electroylyser has a minimal operational capacity (10%). Therefore its potential energy draw at time t needs to be something like;

if (SE[t] + BE[t]) < E_MinCons_kwh:
   (SE[t] + BE[t]) == 0
elif (HZ_Size_nm3 - HSC[t]) < E_MinProd_nm3:
   (SE[t] + BE[t]) == 0
else:
   (SE[t] + BE[t]) <= e_kwh
   (SE[t] + BE[t]) <= (HZ_Size_nm3 - HSC[t]) * nm3_to_kwh)
   (SE[t] + BE[t]) >= E_MinCons_kwh

Where

SE[t] = electricity input from Solar powering the Electrolyser at time t

BE[t] = electricity input from Battery powering the Electrolyser at time t

E_MinCons_kwh = The minimal power draw of the Electrolyser (~18 kWh)

e_kwh = the maximum power draw of the Electrolyser (~180 KWh)

HZ_Size_nm3 = The size of the hydrogen Storage tank (~5,000 nm3)

HSC[t] = the current storage level of the hydrogen tank at time t (in nm3)

HSC[t+1] = HSC[t] + SH[t] + BH[t] - HF[t] - HM[t]

SH[t] = volume of Hydrogen produced from solar electricity at time t (in nm3)

BH[t] = volume of Hydrogen produced from Battery Discharge at time t (in nm3)

HF[t] = volume of Hydrogen consumed by the Fuel-Cell at time t (in nm3)

HM[t] = volume of Hydrogen sold to market at time t (in nm3)

E_MinProd_nm3 = Hydrogen produced when Electroylser operating at min capacity (~3.6 nm3)

nm3_to_kwh = Conversion factor, kwh required to produce 1 nm3 of hydrogen (~5.4 kwh)

Requested Assistance

Does anyone know How I can write the above logical argument, in a method which can be used in PuLP constraints?

I invite you to visit the dedicated Operations Research Stack Exchange site, where you can access a community focused on these topics. Could cross post there. https://or.stackexchange.com/ — dhasson, Jan 22 '20 at 02:56

score 2 · Accepted Answer · answered Jan 22 '20 at 12:05

let me add a logical constraint to the bus example to show you how to use big M with pulp.

import pulp
import cplex

bus_problem = pulp.LpProblem("bus", pulp.LpMinimize)

nbBus40 = pulp.LpVariable('nbBus40', lowBound=0, cat='Integer')
nbBus30 = pulp.LpVariable('nbBus30', lowBound=0, cat='Integer')

# Objective function
bus_problem += 600 * nbBus40 + 480 * nbBus30, "cost"

# Constraints
bus_problem += 40 * nbBus40 + 30 * nbBus30 >= 300

bus_problem.solve(pulp.CPLEX())
print(pulp.LpStatus[bus_problem.status])

for variable in bus_problem.variables():
    print ("{} = {}".format(variable.name, variable.varValue))

print()
print("with if nb buses 40 more than 3  then nbBuses30 more than 7")

M=100

#if then constraint

A = pulp.LpVariable('A', lowBound=0, cat='Binary')
B = pulp.LpVariable('B', lowBound=0, cat='Binary')

bus_problem += A<=B


bus_problem += nbBus40-2<=M*(A)
bus_problem +=nbBus40-3>=-M*(1-A)
bus_problem +=nbBus30-6<=M*(B)
bus_problem +=nbBus30-7>=-M*(1-B)

which gives

Optimal
A = 0.0
B = 1.0
nbBus30 = 10.0
nbBus40 = 0.0

with the docplex API you could write the same

from docplex.mp.model import Model

mdl = Model(name='buses')
nbbus40 = mdl.integer_var(name='nbBus40')
nbbus30 = mdl.integer_var(name='nbBus30')

A = mdl.binary_var(name='A')
B = mdl.binary_var(name='B')


mdl.add_constraint(nbbus40*40 + nbbus30*30 >= 300, 'kids')
mdl.minimize(nbbus40*500 + nbbus30*400)

mdl.solve()

for v in mdl.iter_integer_vars():
   print(v," = ",v.solution_value)

print()
print("with if nb buses 40 more than 3  then nbBuses30 more than 7")

M=100

#if then constraint



mdl.add_constraint(A<=B)


mdl.add_constraint(nbbus40-2<=M*(A))
mdl.add_constraint(nbbus40-3>=-M*(1-A))
mdl.add_constraint(nbbus30-6<=M*(B))
mdl.add_constraint(nbbus30-7>=-M*(1-B))

mdl.minimize(nbbus40*500 + nbbus30*400)

mdl.solve()

for v in mdl.iter_binary_vars():
    print(v," = ",v.solution_value)

for v in mdl.iter_integer_vars():
    print(v," = ",v.solution_value)

but you could also write simply logical constraints:

from docplex.mp.model import Model

mdl = Model(name='buses')
nbbus40 = mdl.integer_var(name='nbBus40')
nbbus30 = mdl.integer_var(name='nbBus30')
mdl.add_constraint(nbbus40*40 + nbbus30*30 >= 300, 'kids')
mdl.minimize(nbbus40*500 + nbbus30*400)

mdl.solve()

for v in mdl.iter_integer_vars():
   print(v," = ",v.solution_value)

print()
print("with if nb buses 40 more than 3  then nbBuses30 more than 7")

#if then constraint
mdl.add_constraint(mdl.if_then(nbbus40>=3,nbbus30>=7))

mdl.minimize(nbbus40*500 + nbbus30*400)

mdl.solve()



for v in mdl.iter_integer_vars():
    print(v," = ",v.solution_value)

Thanks Alex, I hadn't heard of Docplex before, I'll rewrite my query using it, and see how I go! — Bobby Heyer, Jan 22 '20 at 23:25

How to write a conditional constraint in PuLP

1 Answers1