I am writing an optimization for the operation of a battery with a solar pv panel in a domestic home using Pyomo and using the glpk solver. When there is not excess solar the optimisation works perfectly but when there is too much solar I get an uninitialized error. Here's my code:
import numpy as np
import pandas as pd
from pyomo.environ import *
from pulp import *
power = pd.read_csv('Battery 2 Hourly Power flow for the month with time.csv')
solar_forecast = list(power['Solar_(W)'])
demand = list(power['Household_(W)']*(-1.0))
power['Time'] = np.arange(720)
Here is a sample of the data:
Hour Grid_(W) Household_(W) Solar_(W) Storage_(W) Time
00:00 392.095000 -379.960833 0.0 -12.134167 0
01:00 205.441667 -193.115833 0.0 -12.325000 1
02:00 142.570000 -130.807500 0.0 -11.763333 2
03:00 96.456667 -84.829167 0.0 -11.632500 3
04:00 100.940833 -89.330833 0.0 -11.610833 4
05:00 132.357500 -120.559167 0.0 -11.797500 5
06:00 182.733333 -171.605833 0.000000 -11.126667 6
07:00 96.900833 -127.610000 42.106667 -11.397500 7
08:00 35.052500 -228.003333 304.839167 -111.890000 8
09:00 1.385833 -308.922500 489.586667 -182.050833 9
10:00 2086.825833 -2596.776667 492.475000 17.478333 10
11:00 959.740000 -1696.840833 846.346667 -109.245833 11
12:00 2648.243333 -3006.825000 374.156667 -15.579167 12
13:00 1826.446667 -2055.055833 240.062500 -11.453333 13
14:00 1410.041667 -1551.870000 153.488333 -11.658333 14
15:00 536.535000 -572.946667 48.553333 -12.143333 15
16:00 1079.993333 -1069.868333 0.890000 -11.016667 16
17:00 2269.711667 -2260.661667 0.000000 -9.050833 17
Optimization
model = ConcreteModel()
# Define model parameters
model.T = Set(initialize=power.Time.tolist(), ordered=True, doc='hour of week')
model.Ein = Param(initialize=2400, doc='Max rate of power flow (W) in')
model.Eout = Param(initialize=-2400, doc='Max rate of power flow (W) out')
model.Smax = Param(initialize=4800, doc='Max storage (Wh)')
model.Smin = Param(initialize=480, doc='Min storage (Wh)')
model.price = Param(initialize = 0.000124, doc = 'Hourly Price of Electricity')
#Decision Variables
model.E = Var(model.T, domain=Reals) # Battery Power Flow
model.S = Var(model.T, domain=NonNegativeReals) # Battery State of Charge
model.excess_pv = Var(model.T, domain=Reals) # Excess PV Forecasted
## Set Constraints
def storage_state(model,t):
'Storage changes with flows in/out'
if t == model.T.first():
return model.S[t] == (model.Smax + model.Smin)/2
else:
return model.S[t] == (model.S[t-1] - model.E[t-1])
model.charge_state = Constraint(model.T, rule = storage_state)
def max_charge(model,t):
'Max hourly charge of the battery'
return model.E[t] <= model.Ein
model.pos_charge = Constraint(model.T, rule = max_charge)
def max_discharge(model,t):
'Max hourly discharge of the battery'
return model.E[t] >= model.Eout
model.discharge = Constraint(model.T, rule = max_discharge)
def max_stateofcharge(model,t):
'Max storage in the battery'
return model.S[t] <= model.Smax
model.max_soc = Constraint(model.T, rule = max_stateofcharge)
def min_stateofcharge(model, t):
'Min storage in the battery'
return model.S[t] >= model.Smin
model.min_soc = Constraint(model.T, rule = min_stateofcharge)
def stop_final_discharge(model, t):
'Prevents battery discharging in final hour'
return model.E[t] <= model.S[t]
model.final = Constraint(model.T, rule = stop_final_discharge)
def excess_solar(model, t):
'Calculating the excess solar if the battery is full'
if model.S[t] == model.Smax:
return model.excess_pv[t] == demand[t] - solar_forecast[t]
else:
return model.excess_pv[t] == 0
model.solar_constraint = Constraint(model.T, rule = excess_solar)
def demand_cons(model, t):
'Demand must be postive'
if model.excess_pv[t] <= 0:
return demand[t] - solar_forecast[t] - model.E[t] >= 0
else:
return demand[t] + model.excess_pv[t] - solar_forecast[t] - model.E[t] >= 0
model.demand_constraint = Constraint(model.T, rule = demand_cons)
# objective
cost = sum(model.price*(demand[t] - solar_forecast[t] - model.E[t]) for t in model.T)
model.cost = Objective(expr = cost, sense=minimize)
# solve
solvername='glpk'
solverpath_folder='C:\\w64'
solverpath_exe='C:\\w64\\glpsol'
opt = SolverFactory(solvername,executable=solverpath_exe)
results = opt.solve(model)
print('Cost = €', round(value(model.cost()),2))
If I don't have an excess PV I don't need the variable model.excess_pv = Var(model.T, domain=Reals) and the demand constraint looks like this:
def demand_cons(model, t):
'Demand must be postive'
return demand[t] - solar_forecast[t] - model.E[t] >= 0
model.demand_constraint = Constraint(model.T, rule = demand_cons)
This is the error I am getting when trying to run the code with the data that has excess pv
ERROR: evaluating object as numeric value: S[0]
(object: <class 'pyomo.core.base.var._GeneralVarData'>)
No value for uninitialized NumericValue object S[0]
ERROR: Rule failed when generating expression for constraint solar_constraint
with index 0: ValueError: No value for uninitialized NumericValue object
S[0]
ERROR: Constructing component 'solar_constraint' from data=None failed:
ValueError: No value for uninitialized NumericValue object S[0]
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-4-cd4b3e07af68> in <module>
63 return Constraint.Skip
64
---> 65 model.solar_constraint = Constraint(model.T, rule = excess_solar)
66
67 def demand_cons(model, t):
c:\users\stephen\python\lib\site-packages\pyomo\core\base\block.py in __setattr__(self, name, val)
542 # Pyomo components are added with the add_component method.
543 #
--> 544 self.add_component(name, val)
545 else:
546 #
c:\users\stephen\python\lib\site-packages\pyomo\core\base\block.py in add_component(self, name, val)
1087 _blockName, str(data))
1088 try:
-> 1089 val.construct(data)
1090 except:
1091 err = sys.exc_info()[1]
c:\users\stephen\python\lib\site-packages\pyomo\core\base\constraint.py in construct(self, data)
834 for index in self.index_set():
835 self._setitem_when_not_present(
--> 836 index, self.rule(block, index)
837 )
838 except Exception:
c:\users\stephen\python\lib\site-packages\pyomo\core\base\util.py in __call__(self, parent, idx)
302 return self._fcn(parent, *idx)
303 else:
--> 304 return self._fcn(parent, idx)
305
306
<ipython-input-4-cd4b3e07af68> in excess_solar(model, t)
58 def excess_solar(model, t):
59 'Calculating the excess solar if the battery is full'
---> 60 if model.S[t] == model.Smax:
61 return model.excess_pv[t] == demand[t] - solar_forecast[t]
62 else:
pyomo\core\expr\logical_expr.pyx in pyomo.core.expr.logical_expr.EqualityExpression.__nonzero__()
pyomo\core\expr\numeric_expr.pyx in pyomo.core.expr.numeric_expr.ExpressionBase.__call__()
c:\users\stephen\python\lib\site-packages\pyomo\core\expr\visitor.py in evaluate_expression(exp, exception, constant)
1052 visitor = _EvaluationVisitor()
1053 try:
-> 1054 return visitor.dfs_postorder_stack(exp)
1055
1056 except ( TemplateExpressionError, ValueError, TypeError,
c:\users\stephen\python\lib\site-packages\pyomo\core\expr\visitor.py in dfs_postorder_stack(self, node)
582 _sub = _argList[_idx]
583 _idx += 1
--> 584 flag, value = self.visiting_potential_leaf(_sub)
585 if flag:
586 _result.append( value )
c:\users\stephen\python\lib\site-packages\pyomo\core\expr\visitor.py in visiting_potential_leaf(self, node)
960
961 if node.is_numeric_type():
--> 962 return True, value(node)
963 elif node.is_logical_type():
964 return True, value(node)
pyomo\core\expr\numvalue.pyx in pyomo.core.expr.numvalue.value()
pyomo\core\expr\numvalue.pyx in pyomo.core.expr.numvalue.value()
ValueError: No value for uninitialized NumericValue object S[0]
Does anyone know how I should formulate the constraint so that it resolves this issue? Any help will be greatly appreciated and thanks in advance for all the help.
As @Erwin has mentioned you can't use a pyomo variable in an if statement. Because at the time of executing the if statement the value of the variable is unknown. See the code comment below outlining the lines with errors.
def excess_solar(model, t):
'Calculating the excess solar if the battery is full'
if model.S[t] == model.Smax: #here is the line with error
return model.excess_pv[t] == demand[t] - solar_forecast[t]
else:
return model.excess_pv[t] == 0
model.solar_constraint = Constraint(model.T, rule = excess_solar)
def demand_cons(model, t):
'Demand must be postive'
if model.excess_pv[t] <= 0: #error here too
return demand[t] - solar_forecast[t] - model.E[t] >= 0
else:
return demand[t] + model.excess_pv[t] - solar_forecast[t] - model.E[t] >= 0
model.demand_constraint = Constraint(model.T, rule = demand_cons)
I'm thinking that if you are using variables to determine how the constraint is going to be formulated there may be an issue in the overall logic. I haven't played around with your example but I am thinking that model.S should be part of the excess_solar constraint rather than used as an if statement. Wouldn't the optimal battery state be impacted by excess solar?
Something similar may be the case for demand constraint.
Maybe something like this for the excess solar constraint:
def excess_solar(model, t):
'Calculating the excess solar if the battery is full'
return model.excess_pv[t] + (model.S[t] - model.Smax) + (demand[t] - solar_forecast[t]) <= 0
If you want to make sure model.excess_pv is 0 when the battery is full you can set a bound on the variable which makes it greater than 0. This will just result in slack on the constraint.