Background - I'm trying to optimize the operation of a battery energy storage system (BESS). I've used the This Github Project as the basis.
Problem - Most BESS systems have a depth of discharge limit. I.e. for a 400 MWh battery, you can only use about 320 MWhs of it. The last 20% is essentially reserved/unusable. I've been trying to add this as an additional constraint but I've had no luck
What I've tried
I've tried putting in bounds to my Battery Storage Param, I've tried adding a constraint (as below) but both have resulted in the error 'No value for uninitialized NumericValue object Ein[0]' (Ein being the param for Energy In)
Attempt 1
model.S = Var(model.T, bounds=(model.Smin, model.Smax))
Attempt 2
def dod_limit(model, t):
"Depth of Discharge limit"
return model.Smin[t] <= model.S[t]
model.dod = Constraint(model.T, rule=dod_limit)
Help Requested
Does anyone know how I can write the constraint for the depth of discharge limit?
Full code for completeness
# ---------------------
# Main Functions
# ---------------------
def model_to_df(model, first_period, last_period):
"""
Create a dataframe with hourly charge, discharge, state of charge, and
price columns from a pyomo model. Only uses data from between the first
(inclusive) and last (exclusive) hours.
Parameters
----------
model : pyomo model
Model that has been solved
first_period : int
First hour of data to keep
last_period: int
The final hour of data to keep
Returns
-------
dataframe
"""
# Need to increase the first & last hour by 1 because of pyomo indexing
# and add 1 to the value of last model hour because of the range
# second model hour by 1
periods = range(model.T[first_period + 1], model.T[last_period + 1] + 1)
Ein = [value(model.Ein[i]) for i in periods]
Eout = [value(model.Eout[i]) for i in periods]
rrp = [model.P.extract_values()[None][i] for i in periods]
charge_state = [value(model.S[i]) for i in periods]
df_dict = dict(
period=periods,
Ein=Ein,
Eout=Eout,
rrp=rrp,
charge_state=charge_state
)
df = pd.DataFrame(df_dict)
return df
def optimize_year(df, mwh, mw, ef, soc, first_model_period, last_model_period):
"""
Optimize the charge/discharge behavior of a battery storage unit over a
full year. Assume perfect foresight of electricity prices. The battery
has a discharge constraint equal to its storage capacity and round-trip
efficiency of 85%.
Parameters
----------
df : dataframe
dataframe with columns of hourly LBMP and the hour of the year
mwh : int
size of bess
mw : int
storage capacity of bess
ef : double
round trip efficiency of bess
soc : double
State of Charge
first_model_period : int, optional
Set the first hour of the year to be considered in the optimization
(the default is 0)
last_model_period : int, optional
Set the last hour of the year to be considered in the optimization (the
default is 8759)
Returns
-------
dataframe
hourly state of charge, charge/discharge behavior, lbmp, and time stamp
"""
# Filter the data
df = df.loc[first_model_period:last_model_period, :]
model = ConcreteModel()
# Define model parameters
model.T = Set(initialize=df.period.tolist(), doc='periods of year', ordered=True)
model.Rmax = Param(initialize=mw / 2, doc='Max rate of power flow (MW) in or out', within=Any)
model.Smax = Param(initialize=mwh, doc='Max storage (MWh)', within=Any)
model.Smin = Param(initialize=mwh * 0.2, doc='Min storage (MWh)', within=Any)
model.Dmax = Param(initialize=mwh * np.sqrt(ef) * 0.8, doc='Max discharge in 24 hour period', within=Any)
model.P = Param(initialize=df.rrp.tolist(), doc='price for each period', within=Any)
eta = ef # Round trip storage efficiency
# Charge, discharge, and state of charge
# Could use bounds for the first 2 instead of constraints
model.Ein = Var(model.T, domain=NonNegativeReals)
model.Eout = Var(model.T, domain=NonNegativeReals)
model.S = Var(model.T, bounds=(model.Smin, model.Smax))
# Set all constraints
def storage_state(model, t):
'Storage changes with flows in/out and efficiency losses'
# Set first period state of charge to 0
if t == model.T.first():
return model.S[t] == soc
else:
return (model.S[t] == (model.S[t - 1]
+ (model.Ein[t - 1] * np.sqrt(eta))
- (model.Eout[t - 1] / np.sqrt(eta))))
model.charge_state = Constraint(model.T, rule=storage_state)
def dod_limit(model, t):
"Depth of Discharge limit"
return model.Smin[t] <= model.S[t]
model.dod = Constraint(model.T, rule=dod_limit)
def discharge_constraint(model, t):
"Maximum dischage within a single period"
return model.Eout[t] <= model.Rmax
model.discharge = Constraint(model.T, rule=discharge_constraint)
def charge_constraint(model, t):
"Maximum charge within a single period"
return model.Ein[t] <= model.Rmax
model.charge = Constraint(model.T, rule=charge_constraint)
# Without a constraint the model would discharge in the final period
# even when SOC was 0.
def positive_charge(model, t):
'Limit discharge to the amount of charge in battery, including losses'
return model.Eout[t] <= model.S[t] * np.sqrt(eta)
model.positive_charge = Constraint(model.T, rule=positive_charge)
def discharge_limit(model, t):
"Limit on discharge within a 24 hour period"
max_t = model.T.last()
# Check all t until the last 24 hours
if t < max_t-24:
return sum(model.Eout[i] for i in range(t, t+24)) <= model.Dmax
else:
return Constraint.Skip
model.limit_out = Constraint(model.T, rule=discharge_limit)
# Define the battery income, expenses, and profit
income = sum(df.loc[t, 'rrp'] * model.Eout[t] for t in model.T)
expenses = sum(df.loc[t, 'rrp'] * model.Ein[t] for t in model.T)
profit = income - expenses
model.objective = Objective(expr=profit, sense=maximize)
# Solve the model
solver = SolverFactory('glpk')
solver.solve(model)
results_df = model_to_df(model, first_period=first_model_period, last_period=last_model_period)
results_df['local_time'] = df.loc[:, 'local_time']
return results_df
