Ограничить сумму параметров и var на любом временном шаге
В настоящее время я пишу сценарий Python/Pyomo, который оптимизирует совместное использование системы хранения аккумуляторов и солнечной электростанции.
Я хочу ограничить совокупную выработку энергии солнечной электростанции и батареи так, чтобы на каждом временном шаге она была ниже определенного постоянного значения (т. Е. Максимально допустимой выработки электроэнергии солнечной электростанции).
Ошибка, которую я получаю с этим ограничением, заключается в том, что оптимизация завершается, но переменные Pyomo, похоже, вообще не инициализируются.
Мой код Pyomo отлично работает со всеми другими ограничениями, которые я добавил, но когда я добавляю следующее ограничение («maxSingleDischarge»), решатель (glpk), похоже, дает сбой (но завершает работу?) или вообще ничего не делает. Вот короткая выдержка из кода с соответствующими фрагментами:
model = ConcreteModel()
model.T = Set(doc='hour of year', initialize=df.hour.tolist(), ordered=True) # time index
model.PV = Param(model.T,initialize=PV_curve) # time series of PV production (constant, initialized with a list of int)
model.Ein = Var(model.T, domain=NonNegativeReals) # Charged energy of the battery (optimized)
model.Eout = Var(model.T, domain=NonNegativeReals) # Discharged energy of the battery (optimized)
model.DSys = Param(initialize=maxDischarge_System) #Parameter setting the maximum allowed combined discharge (initialized with a constant int)
def maxSingleDischarge(model, t):
return (model.Eout[t] + model.PV[t]) <= model.DSys
model.maxSingleDischarge = Constraint(model.T, rule=maxSingleDischarge)
Как я упоминал ранее, добавление этого ограничения приводит к тому, что такие переменные, как Ein или Eout, вообще не инициализируются. Когда я хочу посмотреть на результаты оптимизации, я получаю ошибки, подобные:
Traceback (most recent call last):
File "<input>", line 2, in <module>
File "<input>", line 2, in <listcomp>
File "pyomo\core\expr\numvalue.pyx", line 153, in pyomo.core.expr.numvalue.value
File "pyomo\core\expr\numvalue.pyx", line 140, in pyomo.core.expr.numvalue.value
ValueError: No value for uninitialized NumericValue object Ein[0]
Другие ограничения работают нормально, как, например, ограничение ниже, где я ограничиваю сумму разряженной энергии от батареи за последние 24 часа определенным значением (model.Dmax). Интересно, что это ограничение, по моему мнению, определяется очень похожим образом, но даже немного более сложным, чем то, с которым я борюсь.
def timed_discharge_limit(model, t):
"Limit on discharge within a 24 hour period"
max_t = model.T.last()
if t < max_t - 24:
return sum(model.Eout[i] for i in range(t, t + 24)) <= model.Dmax
else:
return Constraint.Skip
if model.Dmax != 0:
model.limit_out = Constraint(model.T, rule=timed_discharge_limit)
Мой вопрос: где я допустил ошибку при реализации ограничения «maxSingleDischarge»? Есть ли лучший способ определить это ограничение? Надеюсь, я добавил достаточно информации, чтобы вы могли понять, в чем может быть проблема. Я благодарен за любые предложения, как правильно добавить это ограничение.
Редактировать:
Поскольку это было запрошено, я добавляю сокращенную версию кода с некоторыми данными, с которыми воспроизводится ошибка:
Основной файл, включающий данные
from stackoverflow_bat import *
time_index = np.arange(0,24*7)
#PV production [kW]
PV_curve =[5410.47, 5530.605000000001, 5620.59, 5638.1849999999995, 5631.915,
5625.72, 5625.48, 5625.42, 5625.3150000000005, 5625.225, 5619.195,
5625.225, 5631.285, 5625.465, 5631.599999999999, 5631.735, 5631.6900000000005,
5631.96, 5638.304999999999, 5740.65, 5824.485, 5818.35, 5818.32, 5824.215,
5818.23, 5818.26, 5818.320000000001, 5866.32, 5860.17, 5866.125, 5866.005,
5836.514999999999, 5859.27, 5859.09, 5859.0, 5853.03, 5859.15, 5853.195,
5847.27, 5853.254999999999, 5853.375, 5871.599999999999, 5955.629999999999,
5979.51, 6009.525, 6009.419999999999, 6009.299999999999, 6009.195, 6009.24,
6009.3, 6009.375, 6021.525, 6051.4349999999995, 6051.375, 6051.254999999999,
6051.15, 6056.985, 6056.73, 6056.549999999999, 6050.370000000001, 6050.175,
6050.04, 6049.875, 6013.665, 5959.65, 5953.74, 5959.92, 5953.9349999999995,
5960.084999999999, 5960.175, 5954.160000000001, 5960.085, 5960.085000000001,
5953.949999999999, 5953.875, 5953.77, 5941.665, 5881.65, 5869.74, 5869.845,
5869.9349999999995, 5875.769999999999, 5869.709999999999, 5875.665,
5851.634999999999, 5821.679999999999, 5821.844999999999, 5821.9349999999995,
5822.055, 5822.325, 5828.730000000001, 5835.165, 5907.525, 5967.464999999999,
5985.135, 5972.775, 5972.52, 5966.13, 5953.725, 5857.44, 5803.5, 5779.59,
5779.74, 5785.799999999999, 5785.814999999999, 5779.635, 5779.545, 5767.515,
5701.500000000001, 5695.785, 5695.9349999999995, 5696.070000000001,
5696.1900000000005, 5702.264999999999, 5702.295, 5702.355, 5696.475, 5702.46,
5702.370000000001, 5696.219999999999, 5696.070000000001, 5695.949999999999,
5695.799999999999, 5695.575, 5683.41, 5593.320000000001, 5575.44, 5575.635,
5575.679999999999, 5569.530000000001, 5575.485, 5569.619999999999, 5581.755,
5569.77, 5581.844999999999, 5575.785, 5575.71, 5569.665, 5575.605000000001,
5569.485, 5551.365, 5521.41, 5521.5, 5521.41, 5509.35, 5473.35, 5449.365000000001,
5443.469999999999, 5443.56, 5443.5, 5443.455, 5437.469999999999, 5425.349999999999,
5371.245, 5353.425, 5359.56, 5365.695, 5359.785, 5365.875, 5365.95, 5365.98,
5359.995, 5366.115, 5360.13, 5366.175, 5360.16, 5366.115000000001, 5359.965]
# energy price [€/MWh]
price = [ 1.2060e+01, -1.0000e-01, -6.6000e-01, -1.9900e+00, -2.4200e+00,
-1.8600e+00, -3.3500e+00, 0.0000e+00, -2.2000e-01, 2.9000e-01,
7.3000e-01, 1.1400e+00, 1.3400e+00, 1.5300e+00, 5.7800e+00,
2.3850e+01, 3.6540e+01, 4.6030e+01, 5.2710e+01, 5.4950e+01,
4.9050e+01, 4.4140e+01, 4.5960e+01, 3.5000e+01, 8.2060e+01,
7.3030e+01, 5.2930e+01, 4.6120e+01, 5.0080e+01, 7.8100e+01,
1.0508e+02, 1.4064e+02, 1.4708e+02, 1.4923e+02, 1.4561e+02,
1.4335e+02, 1.4438e+02, 1.4376e+02, 1.4820e+02, 1.5534e+02,
1.6300e+02, 1.7132e+02, 1.7474e+02, 1.6490e+02, 1.5300e+02,
1.4167e+02, 1.3491e+02, 1.2422e+02, 1.3001e+02, 1.2000e+02,
1.1876e+02, 1.1500e+02, 1.1363e+02, 1.1627e+02, 1.3892e+02,
1.5992e+02, 1.6917e+02, 1.7118e+02, 1.6857e+02, 1.6949e+02,
1.6450e+02, 1.6091e+02, 1.6306e+02, 1.7142e+02, 1.7499e+02,
1.7127e+02, 1.6902e+02, 1.6099e+02, 1.4962e+02, 1.3469e+02,
1.1944e+02, 1.0397e+02, 1.0510e+02, 9.0850e+01, 8.3100e+01,
8.3060e+01, 9.3940e+01, 1.1508e+02, 1.5995e+02, 1.6001e+02,
1.5504e+02, 1.5792e+02, 1.5792e+02, 1.7413e+02, 1.6745e+02,
1.5745e+02, 1.5606e+02, 1.4974e+02, 1.5799e+02, 1.7856e+02,
1.6736e+02, 1.7756e+02, 1.7589e+02, 1.4997e+02, 1.4797e+02,
1.0382e+02, 1.0801e+02, 9.7780e+01, 9.5090e+01, 5.3850e+01,
5.6350e+01, 7.5610e+01, 1.0498e+02, 1.1155e+02, 1.2394e+02,
1.4213e+02, 1.4003e+02, 1.4580e+02, 1.4235e+02, 1.4281e+02,
1.5103e+02, 1.5948e+02, 1.7273e+02, 1.9202e+02, 1.9467e+02,
1.7798e+02, 1.6399e+02, 1.5200e+02, 1.4464e+02, 1.2700e+02,
1.0721e+02, 1.0110e+02, 9.6730e+01, 8.6760e+01, 9.9090e+01,
9.9970e+01, 1.0504e+02, 1.4203e+02, 1.4413e+02, 1.4190e+02,
1.3894e+02, 1.2344e+02, 1.2222e+02, 1.1609e+02, 1.3942e+02,
1.5591e+02, 1.6391e+02, 1.5593e+02, 1.6350e+02, 1.4674e+02,
1.3939e+02, 1.3190e+02, 1.2994e+02, 1.2151e+02, 1.0808e+02,
1.0505e+02, 1.0103e+02, 1.0502e+02, 1.0806e+02, 1.3931e+02,
1.5592e+02, 1.5594e+02, 1.5251e+02, 1.6390e+02, 1.5991e+02,
1.6221e+02, 1.6393e+02, 1.6394e+02, 1.6481e+02, 1.8211e+02,
1.6867e+02, 1.7498e+02, 1.7840e+02, 1.6690e+02, 1.6002e+02,
1.5596e+02, 1.5528e+02, 1.0808e+02]
qt = 1 # factor for adjusting time index spacing relative to 1 hour;
# qt = 1 means the time index represents hours
maxStorage = 2000 # maximum capacity of the battery [kWh]
maxFlow = 2000 # maximum charging/discharging power of the battery [kW]
maxDischargePerDay = 4000 # maximum amount of energy allowed to be discharged within a day [kW]
efficiency = 0.9 # round trip efficiency of the battery
minSoC = 0 # minimum allowed State of Charge of the battery [kWh]
maxSoC = 2000 # maximum allowed State of Charge of the battery [kWh]
df = pd.DataFrame(data = {"time_index":time_index, "price": price})
results = optimize_year(df,PV_curve, last_model_time_index=df["time_index"].values[-1],maxFlow=maxFlow,
maxStorage=maxStorage,maxDischargePerDay=maxDischargePerDay,
efficiency=efficiency, minSoC=minSoC, maxSoC=maxSoC,qt=qt,
maxDischarge_System=13000,allowed_load=2000)
И файл stackoverflow_bat, отвечающий за оптимизацию батареи:
import numpy as np
import pandas as pd
from pyomo.environ import *
def model_to_df(model, first_time_index, last_time_index):
time_index = range(model.T[first_time_index + 1], model.T[last_time_index + 1] + 1)
Ein = [value(model.Ein[i]) for i in time_index]
Eout = [value(model.Eout[i]) for i in time_index]
price = [model.P.extract_values()[None][i] for i in time_index]
SoC = [value(model.SoC[i]) for i in time_index]
df = {"time_index":time_index,
"Ein":Ein,
"Eout":Eout,
"price":price,
"SoC":SoC}
return pd.DataFrame(data=df)
def optimize_year(df,PV_curve, first_model_time_index=0, last_model_time_index=8759, maxFlow = 4000,
maxStorage=16000 ,maxDischargePerDay=1000, efficiency=0.9, minSoC = 0, maxSoC = 1, qt=1,
maxDischarge_System=None, allowed_load=0):
# reduce data if necessary
df = df.loc[first_model_time_index:last_model_time_index, :]
model = ConcreteModel()
# Model parameters
model.T = Set(doc='time index', initialize=df.time_index.tolist(), ordered=True)
model.Rmax = Param(initialize=maxFlow, doc='max charging or discharging power (kW)')
model.Smax = Param(initialize=maxStorage, doc='Max capacity (kWh)')
model.Dmax = Param(initialize=maxDischargePerDay, doc='Max amount of discharge in a day (24*qt time_index periode)')
model.P = Param(initialize=df.price.tolist(), doc='prices')
if not maxDischarge_System is None:
model.DSys = Param(initialize=maxFlow, doc='Max allowed discharge by whole system: battery + PV (kW)')
else:
model.DSys = Param(initialize=maxDischarge_System, doc='Max allowed discharge by whole system: battery + PV (kW)')
model.PV = Param(model.T,initialize=PV_curve[0:len(df)], doc='PV production')
model.maxCharge = Param(model.T,initialize=np.array(PV_curve[0:len(df)])+allowed_load, doc='Max possible charge by whole system (kW)')
# Charge, discharge, and SoC to be optimized
model.Ein = Var(model.T, domain=NonNegativeReals)
model.Eout = Var(model.T, domain=NonNegativeReals)
model.SoC = Var(model.T, bounds=(0, model.Smax))
# Constraints
def SoC(model, t):
# Battery State of Charge, including efficiency losses for charging and discharging
# Set state of charge at first time index to half of max SoC
if t == model.T.first():
return model.SoC[t] == model.Smax/2
else:
return (model.SoC[t] == (model.SoC[t - 1]
+ (model.Ein[t - 1] * np.sqrt(efficiency))
- (model.Eout[t - 1] / np.sqrt(efficiency))))
model.charge_state = Constraint(model.T, rule=SoC)
def model_end_state(model, t):
# Set state of charge at last time index to half of max SoC
if t == model.T.last():
return model.SoC[t] == model.Smax / 2
else:
return Constraint.Skip
model.charge_end_state = Constraint(model.T, rule=model_end_state)
def minimum_SoC(model, t):
# minimum allowed SoC
return model.SoC[t] >= minSoC
model.minimum_SoC = Constraint(model.T, rule=minimum_SoC)
def maximum_SoC(model, t):
# maximum allowed SoC
return model.SoC[t] <= maxSoC
model.maximum_SoC = Constraint(model.T, rule=maximum_SoC)
def max_discharge_constraint(model, t):
# maximum allowed discharge by the battery alone
return model.Eout[t] <= model.Rmax
model.max_discharge = Constraint(model.T,rule=max_discharge_constraint)
def charge_constraint(model, t):
#Maximum charge within a single time_index: Limiting factor are PV production and the amount the battery is
#allowed to charge
if model.Rmax.value <= model.maxCharge[t]:
constraint = model.Rmax.value
else:
constraint = model.maxCharge[t]
return model.Ein[t] <= constraint
model.charge = Constraint(model.T, rule=charge_constraint)
def discharge_constraint(model, t):
# Limiting the amount of energy discharged by the battery to the amount of energy stored in the battery
# (including losses)
return model.Eout[t] <= model.SoC[t] / np.sqrt(efficiency)
model.discharge_constraint = Constraint(model.T, rule=discharge_constraint)
def neg_price_constraint(model, t):
# Prohibit discharging of the battery when prices are negative
if df.loc[t, 'price'] < 0:
return model.Eout[t] <= 0
else:
return Constraint.Skip
model.neg_price_constraint = Constraint(model.T, rule=neg_price_constraint)
def discharge_limit(model, t):
# Limit the amount of energy discharged within a specific timeframe (24*qt)
max_t = model.T.last()
if t < max_t - 24 * qt:
return sum(model.Eout[i] for i in range(t, t + 24*qt)) <= model.Dmax
else:
return Constraint.Skip
if model.Dmax > 0:
model.discharge_limit = Constraint(model.T, rule=discharge_limit)
def maxSingleDischarge(model, t):
# Limit on discharge of the whole system within a timestep
return (model.Eout[t] + model.PV[t]) <= model.DSys
model.maxSingleDischarge = Constraint(model.T, rule=maxSingleDischarge)
# Profit function to be optimized
profit = sum((model.PV[t] - model.Ein[t] + model.Eout[t])/1000*df.loc[t, 'price'] for t in model.T)
model.objective = Objective(expr=profit, sense=maximize)
solver = SolverFactory('glpk')
solver.solve(model)
results = model_to_df(model, first_time_index=first_model_time_index,
last_time_index=last_model_time_index)
return results
Когда ограничение maxSingleDischarge удалено
#model.maxSingleDischarge = Constraint(model.T, rule=maxSingleDischarge)
оптимизация работает отлично. В противном случае я получаю вышеупомянутую ошибку.
Это не самый элегантный код из когда-либо написанных, и я уверен, что есть что улучшить, но надеюсь, что он поможет выявить проблему в моем коде. Спасибо!
Спасибо @AirSquid. Я воздержался от добавления слишком большого количества кода, потому что беспокоился, что люди не захотят работать с большим количеством кода. Я добавил сокращенную версию, включающую некоторые данные, воспроизводящие ошибку.
Ответы 2
Возможно ли, что ваша модель неосуществима и поэтому не имеет значения, связанного с Ein[0]?
Если model.PV[t] больше model.DSys, это приведет к тому, что модель будет неосуществимой, поскольку при любом значении model.Eout[t] ограничение не может быть удовлетворено.
Спасибо за ваш ответ. Это хороший момент, но я проверил модель.PV, и ее максимальное значение в любой момент времени составляет около 12,5 МВт. Если для model.DSys установлено значение 13 МВт, model.PV всегда должен быть меньше, а model.Eout должен быть свободным для оптимизации в рамках всего, что доступно между model.PV и model.DSys. Или я здесь что-то не понимаю?
@Voric, не могли бы вы также предоставить результаты решателя glpk для ответа на вопрос, возможно, это поможет лучше понять проблему
Вот что происходит: когда вы добавляете это ограничение, ваша модель становится неосуществимой. Решатель иногда может прийти к такому выводу, не прибегая к каким-либо махинациям с вводом значений в ваши переменные, и в этом случае, когда вы собираетесь проверить/использовать результаты, он выводит из строя, потому что переменная (и, вероятно, многие другие) никогда не используются.
Вы совершаете один из главных грехов оптимизации: не проверяете состояние решателя перед тем, как продолжить. Есть несколько способов сделать это. Если это отдельное решение, вы можете распечатать полученные результаты. Вы также можете (если интересно) выбрать tee=True в решении, чтобы увидеть журнал на экране. Или, если вы намерены получить результат и сделать с ним что-нибудь, как есть, вы можете использовать встроенную функцию pyomo, чтобы проверить результат, как показано ниже. Обратите внимание: если вы подставите этот код в свой в операторе решения, он распечатает ограничение, с которым у вас возникли проблемы, и вы довольно быстро поймете, почему это невозможно.
solver = SolverFactory('glpk')
solve_res = solver.solve(model)
# print it: good idea!
print(solve_res)
# check it:
if check_optimal_termination(solve_res): # <--- built in pyomo function
results = model_to_df(model, first_time_index=first_model_time_index,
last_time_index=last_model_time_index)
return results
else:
# non-optimal termination.... do something
model.maxSingleDischarge.pprint()
return None
Это очень помогло в поиске проблемы. Большое спасибо, я не знал об этих инструментах для проверки результатов решателя, но с этого момента я буду строго их использовать!
Здесь происходит что-то еще. Это поможет вам получить лучший ответ, если вы опубликуете достаточно кода, чтобы воспроизвести ошибку. Я бы посоветовал вам отредактировать (или просто добавить) свой пост и включить достаточно вашей модели (и некоторых игрушечных данных), чтобы воспроизвести ошибку.