如何在 PuLP 中编写条件约束
How to write a conditional constraint in PuLP
背景
我正在建立一个模型来优化使用 PuLP 的太阳能 + 电池 + electrolyser/fuel-cell 系统的调度。这个想法是通过在价格低时为电池充电并通过太阳能发电生产氢气，然后在价格高时对电池放电并燃烧氢气来最大化收入。
研究与故障排除
我是 PuLP 的新手，我花了一段时间来研究这个问题，使用在线指南和 SO 寻求帮助，但我的代码仍然返回垃圾答案。我花了一段时间浏览 SO 来寻找答案 (here, , and )。我想我知道可能是什么问题了。
潜在问题识别
我认为我的主要问题是我有一堆逻辑约束，我需要以不同的形式编写它们才能在 PuLP 中工作。我目前没有指标变量或 Big M 约束，我认为这是编写条件约束的有用方法。我试着自己写，但没有成功，因为我不知道如何用 PuLP 写它们（我相信它不允许你使用 IF 或 np.where，等等）
至于 Big M 约束，我尝试过以我认为合适的方式编写它们，但是 returns 错误 "Non-constant expressions cannot be multiplied".
问题示例
我在下面附上了我试图通过 PuLP 中的条件约束实现的代码示例。
在我的系统中，电解槽的运行容量最小 (10%)。因此，它在时间 t 的势能消耗需要类似于；
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
哪里
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)
请求帮助
有谁知道我如何在可用于 PuLP 约束的方法中编写上述逻辑参数？
让我给 bus example 添加一个逻辑约束，向您展示如何使用 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)
这给出了
Optimal
A = 0.0
B = 1.0
nbBus30 = 10.0
nbBus40 = 0.0
用 docplex API 你可以写同样的东西
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)
但您也可以编写简单的逻辑约束：
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)
背景
我正在建立一个模型来优化使用 PuLP 的太阳能 + 电池 + electrolyser/fuel-cell 系统的调度。这个想法是通过在价格低时为电池充电并通过太阳能发电生产氢气，然后在价格高时对电池放电并燃烧氢气来最大化收入。
研究与故障排除
我是 PuLP 的新手，我花了一段时间来研究这个问题，使用在线指南和 SO 寻求帮助，但我的代码仍然返回垃圾答案。我花了一段时间浏览 SO 来寻找答案 (here,
潜在问题识别
我认为我的主要问题是我有一堆逻辑约束，我需要以不同的形式编写它们才能在 PuLP 中工作。我目前没有指标变量或 Big M 约束，我认为这是编写条件约束的有用方法。我试着自己写，但没有成功，因为我不知道如何用 PuLP 写它们（我相信它不允许你使用 IF 或 np.where，等等）
至于 Big M 约束，我尝试过以我认为合适的方式编写它们，但是 returns 错误 "Non-constant expressions cannot be multiplied".
问题示例
我在下面附上了我试图通过 PuLP 中的条件约束实现的代码示例。
在我的系统中，电解槽的运行容量最小 (10%)。因此，它在时间 t 的势能消耗需要类似于；
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
哪里
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)
请求帮助
有谁知道我如何在可用于 PuLP 约束的方法中编写上述逻辑参数？
让我给 bus example 添加一个逻辑约束，向您展示如何使用 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)
这给出了
Optimal
A = 0.0
B = 1.0
nbBus30 = 10.0
nbBus40 = 0.0
用 docplex API 你可以写同样的东西
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)
但您也可以编写简单的逻辑约束：
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)