I'm currently writing a Python/Pyomo script that optimizes the combined usage of a battery storage system and a solar power plant.
I want to constrain the combined energy output of the solar power plant and the battery so that it is at every timestep lower than a specific constant value (i.e., the maximum allowed power generation of the solar power plant).
The error that I get with this constraint is that the optimization finishes but Pyomo vars seemingly did not get initialized at all.
My Pyomo code works fine with all the other constrains that I added, but when I add the following constraint ("maxSingleDischarge") the solver (glpk) seems to crash (but finishes?) or it doesn't do anything at all. Here's a short extract from the code, with the relevant pieces:
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)
As I mentioned earlier, adding this constraint causes vars like Ein or Eout to not be initialized at all. When I want to look at the optimization results, I get errors similar to:
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]
Other constraints work fine, as for example the constraint below, where I limit the sum of discharged energy from the battery within the last 24 hours to a specific value (model.Dmax). Interestingly, this constraint is in my opinion defined in a very similar way, but even a bit more complex than the one I struggle with.
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)
My question is, where have I made a mistake in implementing the constraint "maxSingleDischarge"? Is there a better way to define that constraint? I hope I added enough information so that you might have an idea what the problem could be. I'm thankful for any suggestion how to properly add this constraint.
Edit:
Because it was requested, I add a shortened version of the code with some data with which the error is reproducable:
A main file including the data
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)
And a file "stackoverflow_bat" handling the battery optimization:
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
When the constraint maxSingleDischarge is removed
#model.maxSingleDischarge = Constraint(model.T, rule=maxSingleDischarge)
the optimization works perfectly fine. Otherwise I get aforementioned error.
It is not the most elegant code ever written and I'm sure there are lots of things to improve, but I hope it helps with identifying the problem in my code. Thank you!
Here's what's happening: When you add that constraint in, you model is infeasible. The solver can sometimes deduce that without going through any machinations of putting values in your variables, and in that case, when you go to inspect/use the results, it barfs because the variable (and likely many others) are never used.
You are committing one of the cardinal sins of optimization: not checking the solver status before proceeding. There are a couple ways to do this. If it is a stand-alone solve, you can print the captured results. You could also (if interested) select tee=True in the solve to see the log on screen. Or, if you intend to grab the result and do something with it as you are, you can use pyomo's built in function to check the result as shown below. Note, if you substitute this code into yours at the solve statement, it will print out the constraint you are having trouble with and you'll see why it is infeasible rather quickly.
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