I'm currently trying to model an energy system model with V2G storage as my first go with pyomo. The aim of the model is to minimize costs while satisfying the demand for electricity. This can either be done through intermittent generation from solar, onshore wind and offshore wind or by discharging electricity from electric vehicles. I have attached a picture of the equations of the model. Model equations However, I'm experiencing some problems with my code and get the error code:
289093 lines were written
GLPK Simplex Optimizer 5.0
52562 rows, 35044 columns, 140160 non-zeros
Preprocessing...
26279 rows, 26280 columns, 105113 non-zeros
Scaling...
A: min|aij| = 2.860e-05 max|aij| = 1.107e+00 ratio = 3.872e+04
GM: min|aij| = 7.727e-03 max|aij| = 1.294e+02 ratio = 1.675e+04
EQ: min|aij| = 5.971e-05 max|aij| = 1.000e+00 ratio = 1.675e+04
Constructing initial basis...
Size of triangular part is 26279
0: obj = 8.061601191e+09 inf = 1.219e+07 (23830)
7587: obj = 5.758397626e+09 inf = 6.047e+06 (16900) 69
14917: obj = 5.758397626e+09 inf = 3.688e+06 (11839) 29
LP HAS NO PRIMAL FEASIBLE SOLUTION
glp_simplex: unable to recover undefined or non-optimal solution
If you need actual output for non-optimal solution, use --nopresol
...
N_of = 0
ERROR: evaluating object as numeric value: V_discharge[0]
(object: <class 'pyomo.core.base.var._GeneralVarData'>)
No value for uninitialized NumericValue object V_discharge[0]
Output exceeds the size limit. Open the full output data in a text editor---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_15212\450560137.py in
80 print("N_on =", model.N_on())
81 print("N_of =", model.N_of())
---> 82 print("Z =", model.Z())
83
c:\Users\emok\Anaconda3\lib\site-packages\pyomo\core\base\expression.py in __call__(self, exception)
59 if self.expr is None:
60 return None
---> 61 return self.expr(exception=exception)
62
63 def is_named_expression_type(self):
c:\Users\emok\Anaconda3\lib\site-packages\pyomo\core\expr\base.py in __call__(self, exception)
113
114 """
--> 115 return evaluate_expression(self, exception)
116
117 def __str__(self):
c:\Users\emok\Anaconda3\lib\site-packages\pyomo\core\expr\visitor.py in evaluate_expression(exp, exception, constant)
1240
1241 try:
...
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 V_discharge[0]
When I initialize V_discharge I just get 0 for all my results. I have a feeling the trouble has to do with the dynamic aspect of my model as if it doesn't update the variable V (energy storage in the EVs) properly. Moreover, I have successfully solved the static model with no storage option.
Therefore, I have the following questions:
All help would greatly be appreciated.
import pandas as pd
from pyomo.environ import \*
from pyomo.opt import SolverFactory
from pyomo.core import Constraint
# Load data into a Pandas DataFrame
data = pd.read_csv(r'D:\data\data.csv')
# Create a Pyomo model
model = ConcreteModel()
# Define the sets
model.hours = Set(initialize=data.index, doc='Set of hours')
# Define the parameters
model.k_ipv = Param(default=40)
model.k_ion = Param(default=30)
model.k_iof = Param(default=46)
model.k_vv2g = Param(default=226.56)
model.eta_charge = Param (default=0.924)
model.eta_discharge = Param(default=0.903)
model.E_s = Param(default=160425)
model.alpha = Param(default=0.5)
#Define data
model.CF_pv = Param(model.hours, initialize=data['CF_pv'].to_dict())
model.CF_on = Param(model.hours, initialize=data['CF_on'].to_dict())
model.CF_of = Param(model.hours, initialize=data['CF_of'].to_dict())
model.pv2g = Param(model.hours, initialize=data['pv2g'].to_dict())
model.d_ev = Param(model.hours, initialize=data['d_ev'].to_dict())
# Define the variables
model.N_pv = Var(domain=NonNegativeReals, initialize=0)
model.N_on = Var(domain=NonNegativeReals, initialize=0)
model.N_of = Var(domain=NonNegativeReals, initialize=0)
model.V = Var(model.hours, domain=NonNegativeReals)
model.V_charge = Var(model.hours, domain=NonNegativeReals)
model.V_discharge = Var(model.hours, domain=NonNegativeReals)
model.C = Var(model.hours, domain=NonNegativeReals)
# Market balance constraints
model.d1 = Constraint(model.hours, rule=lambda model, h: model.CF_pv[h] * model.N_pv +
model.CF_on[h] * model.N_on +
model.CF_of[h] * model.N_of + model.V_discharge[h] == data['d'][h])
model.d2 = Constraint(model.hours, rule=lambda model, h: model.CF_pv[h] * model.N_pv +
model.CF_on[h] * model.N_on +
model.CF_of[h] * model.N_of == data['d'][h] + model.V_charge[h] + model.C[h])
#State of charge constraint
def s1_constraint(model,h):
if h == model.hours.first():
return model.V[h] == model.alpha*model.E_s / 2
else:
return model.V[h] == model.V[h-1] + model.eta_charge * model.V_charge[h] - model.V_discharge[h]/model.eta_discharge - data['d_ev'][h]
model.s1 = Constraint(model.hours, rule = s1_constraint)
"""No free lunch"""
model.s2 = Constraint(rule=lambda model: model.V[model.hours.first()] == model.V[model.hours.last()])
"""Maximum energy storage must be less or equal to capacity"""
model.s3 = Constraint(model.hours, rule=lambda model, h: model.V[h] <= model.alpha*model.E_s)
# Charging and discharging constraints
""" Maximum discharge rate within a single hour """
model.disc1 = Constraint(model.hours, rule=lambda model, h: model.V_discharge[h] <= model.alpha*data['pv2g'][h])
""" Maximum charge rate within a single hour """
model.char1 = Constraint(model.hours, rule=lambda model, h: model.V_charge[h] <= model.alpha*data['pv2g'][h])
# Define objective function
def obj_expression(model):
return model.k_ipv * model.N_pv + model.k_ion * model.N_on + model.k_iof * model.N_of + model.k_vv2g * sum(model.V_discharge[h] for h in model.hours)
model.Z = Objective(expr=obj_expression(model), sense=minimize) # Use the defined objective expression with the 'h' variable
# Solve the optimization problem
solver = SolverFactory('glpk')
#solver.solve(model)
results = solver.solve(model,tee=True)
# Access the solution
print("N_pv =", model.N_pv())
print("N_on =", model.N_on())
print("N_of =", model.N_of())
print("Z =", model.Z())
The error you are getting shows "LP HAS NO PRIMAL FEASIBLE SOLUTION" in the report, which is the key issue. Your model as written is infeasible. I think your formulation is suspect in several areas. A few things to consider:
in constraints d1, d2... If you subtract one of those equations from the other you are left with:
v_discharge = -v_charge - c
which can only be true if they are all zero as you have them defined as non-negative. So there is a problem with your sign convention for charging & discharging.
I'm not sure why you need d2. Typically you just need to cover the demand side of flow balance, I think:
pv + solar + discharge - charge + outside_source >= demand
Your d1 constraint is problematic. Typically this should be a greater than or equals constraint, not an equality. What happens if demand in any particular hour is zero? You've then forced all of your other variables to be zero'ed...
It isn't clear what model.C is. You should comment all of your variables...
If you are still stuck after working on the formulation, you should edit your post with the re-formulation, and cut & paste about 10 rows or so of your .csv file so your error can be reproduced.