I'm trying to create a model optimization for a energy storage system using pyomo. Using the demand in kWh from an household and the electricity prices, I would like to minimize the cost charging and discharging the battery at the right time. I already have a working model for 1 year of data and the model is able to find an optimal solution (see code below). However, when I try find a model for only three months (let's say from October to December), pyomo returns with a termination condition "unbound" but I can't figure out why.
Model for 1 year of data:
#Battery parameters
battery_capacity = 252
total_energy = 13.1
usable_energy = 12.4
nominal_voltage = 51.8
ratio = total_energy/usable_energy
power = usable_energy / ratio
nominal_current = usable_energy / nominal_voltage * 1000
recharging_hours = battery_capacity/nominal_current
battery_level_threshold = 0
model = pyo.ConcreteModel()
#Set time period
model.T = pyo.Set(initialize=pyo.RangeSet(len(df_2021)),ordered=True)
#PARAMETERS:
model.b_efficiency = pyo.Param(initialize=0.9)
model.b_min_cap = pyo.Param(initialize=0)
model.b_max_cap = pyo.Param(initialize=12.4)
model.b_charge_power = pyo.Param(initialize=power)
model.b_discharge_power = pyo.Param(initialize=power)
model.spot_prices = pyo.Param(model.T,initialize=dict(enumerate(df_2021["Price/kWh"],1)),within=pyo.Any)
model.demand = pyo.Param(model.T, initialize=dict(enumerate(df_2021["Demand"],1)),within=pyo.Any)
#Variables : also the variable has to be indexed with the time T
model.b_soc = pyo.Var(model.T, domain = pyo.NonNegativeReals, bounds = (model.b_min_cap, model.b_max_cap))
model.b_discharge = pyo.Var(model.T, domain = pyo.NonNegativeReals)
model.b_charge = pyo.Var(model.T, domain = pyo.NonNegativeReals)
model.elect_purchased = pyo.Var(model.T, domain = pyo.NonNegativeReals)
#CONSTRAINTS
#Purchase constraint
def purchase_constraint(model,t):
return model.elect_purchased[t] >= model.demand[t] - model.b_discharge[t] + model.b_charge[t]
#State of charge constraint
def soc_constraint(model,t):
if t == model.T.first():
return model.b_soc[t] == model.b_max_cap / 2 #- model.b_discharge[t] + model.b_charge[t]
else:
return model.b_soc[t] == model.b_soc[t-1] - model.b_discharge[t-1] + model.b_charge[t-1]
#Discharge and charge constraints
def discharge_constraint_1(model,t):
""" Maximum discharge rate within a single hour """
return model.b_discharge[t] <= model.b_discharge_power
def discharge_constraint_2(model,t):<br/>
""" Sets the maximum energy available to be discharged as the SOC - minimum SOC """
return model.b_discharge[t] <= model.b_soc[t] - model.b_min_cap
def charge_constraint_1(model,t):
""" Maximum charge rate within a single hour """
return model.b_charge[t] <= model.b_charge_power
def charge_constraint_2(model,t):<br/>
""" Sets the maximum energy available to be cahrge as the SOC max """
return model.b_charge[t] <= model.b_max_cap - model.b_soc[t]
model.purchase_c = pyo.Constraint(model.T, rule = purchase_constraint)
model.soc_c = pyo.Constraint(model.T, rule = soc_constraint)<br/>
model.discharge_c1 = pyo.Constraint(model.T,rule = discharge_constraint_1)
model.discharge_c2 = pyo.Constraint(model.T,rule = discharge_constraint_2)
model.charge_c1 = pyo.Constraint(model.T,rule = charge_constraint_1)
model.charge_c2 = pyo.Constraint(model.T,rule = charge_constraint_2)
#OBJECTIVE
expr = sum(model.elect_purchased[t] * model.spot_prices[t] for t in model.T)
model.objective = pyo.Objective(rule = expr, sense = pyo.minimize)
opt = pyo.SolverFactory('cbc',executable='/usr/bin/cbc')
results = opt.solve(model)
results.write()
Result 1 year data
The optimal solution is found
However, when I change the dataset, using the SAME model structure and constraints, pyomo doesn't find a solution.
Model for 3 months:
#Battery parameters
battery_capacity = 252
total_energy = 13.1
usable_energy = 12.4
nominal_voltage = 51.8
ratio = total_energy/usable_energy
power = usable_energy / ratio
nominal_current = usable_energy / nominal_voltage * 1000
recharging_hours = battery_capacity/nominal_current
battery_level_threshold = 0
model = pyo.ConcreteModel()
#Set time period
model.T = pyo.Set(initialize=pyo.RangeSet(len(df_2021)),ordered=True)
#PARAMETERS:
model.b_efficiency = pyo.Param(initialize=0.9)
model.b_min_cap = pyo.Param(initialize=0)
model.b_max_cap = pyo.Param(initialize=12.4)
model.b_charge_power = pyo.Param(initialize=power)
model.b_discharge_power = pyo.Param(initialize=power)
model.spot_prices = pyo.Param(model.T,initialize=dict(enumerate(df_2021["Price/kWh"],1)),within=pyo.Any)
model.demand = pyo.Param(model.T, initialize=dict(enumerate(df_2021["Demand"],1)),within=pyo.Any)
#Variables : also the variable has to be indexed with the time T
model.b_soc = pyo.Var(model.T, domain = pyo.NonNegativeReals, bounds = (model.b_min_cap, model.b_max_cap))
model.b_discharge = pyo.Var(model.T, domain = pyo.NonNegativeReals)
model.b_charge = pyo.Var(model.T, domain = pyo.NonNegativeReals)
model.elect_purchased = pyo.Var(model.T, domain = pyo.NonNegativeReals)
#CONSTRAINTS
#Purchase constraint
def purchase_constraint(model,t):
return model.elect_purchased[t] >= model.demand[t] - model.b_discharge[t] + model.b_charge[t]
#State of charge constraint
def soc_constraint(model,t):
if t == model.T.first():
return model.b_soc[t] == model.b_max_cap / 2 #- model.b_discharge[t] + model.b_charge[t]
else:
return model.b_soc[t] == model.b_soc[t-1] - model.b_discharge[t-1] + model.b_charge[t-1]
#Discharge and charge constraints
def discharge_constraint_1(model,t):
""" Maximum discharge rate within a single hour """
return model.b_discharge[t] <= model.b_discharge_power
def discharge_constraint_2(model,t):<br/>
""" Sets the maximum energy available to be discharged as the SOC - minimum SOC """
return model.b_discharge[t] <= model.b_soc[t] - model.b_min_cap
def charge_constraint_1(model,t):
""" Maximum charge rate within a single hour """
return model.b_charge[t] <= model.b_charge_power
def charge_constraint_2(model,t):<br/>
""" Sets the maximum energy available to be cahrge as the SOC max """
return model.b_charge[t] <= model.b_max_cap - model.b_soc[t]
model.purchase_c = pyo.Constraint(model.T, rule = purchase_constraint)
model.soc_c = pyo.Constraint(model.T, rule = soc_constraint)<br/>
model.discharge_c1 = pyo.Constraint(model.T,rule = discharge_constraint_1)
model.discharge_c2 = pyo.Constraint(model.T,rule = discharge_constraint_2)
model.charge_c1 = pyo.Constraint(model.T,rule = charge_constraint_1)
model.charge_c2 = pyo.Constraint(model.T,rule = charge_constraint_2)
#OBJECTIVE
expr = sum(model.elect_purchased[t] * model.spot_prices[t] for t in model.T)
model.objective = pyo.Objective(rule = expr, sense = pyo.minimize)
opt = pyo.SolverFactory('cbc',executable='/usr/bin/cbc')
results = opt.solv
Pyomo returns:
WARNING: Loading a SolverResults object with a warning status into
model.name="unknown";
- termination condition: unbounded
- message from solver: <undefined>
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: None
Upper bound: inf
Number of objectives: 1
Number of constraints: 16993
Number of variables: 11329
Number of nonzeros: 2832
Sense: minimize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: warning
User time: -1.0
System time: 0.45
Wallclock time: 0.58
Termination condition: unbounded
Termination message: Model was proven to be unbounded.
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Black box:
Number of iterations: 0
Error rc: 0
Time: 0.6151924133300781
Since the only change between the two runs is the length of the dataframe, I don't really know where to look for the error.
Given that the model works on some data, but not on an alternate data source, we can obviously focus a bit on the data set (which isn't shown).
We have a huge clue in the error report that the problem is unbounded. This means that there is nothing to prevent the objective function from running away to infinity, or negative infinity in the case of a minimization problem. So, let's look at your objective function. You are:
sum(demand[t] * price[t] for t in T)
The domain of your demand variable is set to non-negative real numbers (good) and price[t] is a parameter read in from data, and you are minimizing. So the only way this could run away to negative infinity is if there are one (or more) negative prices in your data, which the solver would then act on and demand would be unbounded.
So, comb your data set. if you are using pandas, just use a logical search for rows where price is < 0, you'll likely find at least one. Then you'll have to decide if that is a typo or if it is realistic to have a negative price (which could happen in some systems), and if it is legit, you will have to impose some other constraint to limit the model in that situation.