如何在PuLP中编写条件约束

Question

背景

我正在建立一个使用PuLP的太阳能+电池+电解质/燃料电池系统的优化调度模型。这样做的目的是在价格较低的情况下，通过给电池充电和利用太阳能发电生产氢气，然后在价格高的时候排放电池和燃烧氢气，从而实现收入最大化。

研究与故障排除

我是PuLP新手，我花了一段时间在这方面工作，使用在线指南等帮助，但我的代码仍在返回垃圾答案。我花了一段时间浏览搜索答案(这里、这里、这里和这里)。我想我知道问题可能是什么。

潜在问题识别

我认为我的主要问题是，我有许多逻辑约束，需要用不同的形式编写才能在PuLP中工作。我目前没有指标变量或大M约束，我认为这是写条件约束的有用方法。我试着自己写，但没有运气，因为我不知道如何用PuLP (我认为它不允许您使用IF或np.where等等)

至于大M约束，我试着用我认为合适的方式编写它们，但它返回了“非常数表达式不能乘以”的错误。

问题示例

下面是我试图通过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约束中使用的方法来编写上述逻辑参数吗？

Answer 1

让我在总线实例中添加一个逻辑约束，向您展示如何使用带纸浆的大M。

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

使用双工 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)