Как написать условное ограничение в PuLP - PullRequest
Новая
       5

Как написать условное ограничение в PuLP

0 голосов
/ 22 января 2020

Справочная информация

Я настраиваю модель для оптимизации отправки системы "солнечная батарея + батарея + электролизер / топливный элемент" с использованием PuLP. Идея состоит в том, чтобы максимизировать доход за счет зарядки батареи и производства водорода из солнечной генерации, когда цены низкие, затем разрядки батареи и сжигания водорода, когда цены высокие.

исследование и устранение неисправностей

Я новичок в PuLP, я потратил некоторое время, работая над этим, используя интерактивные руководства и SO для помощи, но мой код все еще возвращает ответы rubbi sh. Я провел некоторое время, просматривая SO в поисках ответов ( здесь , здесь , здесь и здесь ). и я думаю, что знаю, в чем может быть проблема.

Идентификация потенциальной проблемы

Я думаю, что моя главная проблема в том, что у меня есть куча логических ограничений, которые мне нужно пишите в другой форме для работы в ПУЛП. В настоящее время у меня нет индикаторных переменных или ограничений Big M, которые, на мой взгляд, являются полезным способом записи условных ограничений. Я пытался написать их сам, но безуспешно, потому что я не знал, как написать их в PuLP (что, я считаю, не позволяет использовать IF или np.where, et c)

Что касается ограничений Big M, у меня было go при их написании, что я считаю правильным способом, но оно возвращает ошибку "Неконстантные выражения не могут быть умножены".

Пример задачи

Ниже приведен пример кода, который я пытаюсь достичь с помощью условного ограничения в PuLP.

В моей системе Electroylyser имеет минимальную работоспособность (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?

1 Ответ

1 голос
/ 22 января 2020

позвольте мне добавить логическое ограничение к примеру шины , чтобы показать вам, как использовать большие значения 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

с помощью 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)
...