0

Hi there, I have written a Pyomo script to optimize a dispatch pattern with 3 components (thermal storage, boiler and heat pump). It seems to work until constraints are defined regarding the thermal storage. Could anyone help me understand why there is no optimal solution, and which constraint (its not the emission constraint) is defined wrong? Thank you for your help!

import pyomo.environ as pyo
from pyomo.environ import *
from pyomo.opt import SolverFactory
from pyomo.opt import TerminationCondition

# solver function

solver = SolverFactory("glpk")
model = ConcreteModel()

# creating a timeframe

number_of_period = 8760
model.time_set = pyo.Set(initialize=range(1, number_of_period+1), ordered=True)

# loading csv input files

input_assumptions = pd.read_csv('input_assumptions.csv')
input_price = pd.read_csv('prices.csv')
input_price_new = pd.DataFrame(input_price)
demand = input_assumptions['Total network demand (MWh)']
demand_new = pd.DataFrame(demand)

#parameters

model.Price = pyo.Param(model.time_set, initialize=dict(enumerate(input_price_new['Electricity prices'], 1)), within=pyo.Any)
model.Demand = pyo.Param(model.time_set, initialize=dict(enumerate(demand_new['Total network demand (MWh)'], 1)), within=pyo.Any)
model.HP_max = pyo.Param(model.time_set, initialize=dict(enumerate(input_assumptions['Max Heat pump'], 1)), within=pyo.Any)

#variables

model.b_min_cap = pyo.Param(initialize=0)
model.b_max_cap = pyo.Param(initialize=8.708)

model.b_soc = pyo.Var(model.time_set, domain = pyo.NonNegativeReals, bounds = (model.b_min_cap, model.b_max_cap))
model.b_discharge = pyo.Var(model.time_set, domain = pyo.NonNegativeReals, initialize = 0)
model.b_charge = pyo.Var(model.time_set, domain = pyo.NonNegativeReals, initialize = 0)

def HP_domain(model, i):
    upper_HP = input_assumptions['Max Heat pump']
    lower = input_assumptions['lower']
    return(0, upper_HP[i-1])

model.HP_var = pyo.Var(model.time_set, domain = pyo.NonNegativeReals, bounds=HP_domain, initialize = 0)

model.Boiler_var = pyo.Var(model.time_set, domain = pyo.NonNegativeReals, bounds=(0, 24), initialize = 0)

#definition of the objective

def minimise(model):
    return pyo.summation(model.Price, model.HP_var)

model.obj = pyo.Objective(rule=minimise, sense=pyo.minimize)

#constraints

def Soc_constraint(model, i):
    if i == model.time_set.first():
        return model.b_soc[i] == model.b_max_cap # starting at 100 % charge level 
    else: 
        return model.b_soc[i] == model.b_soc[i-1] - model.b_discharge[i-1] + model.b_charge[i-1]
model.soc_c = pyo.Constraint(model.time_set, rule=Soc_constraint)

def dispatch_rule(model, i):
    return sum(model.b_discharge[i] for i in range(1, 8761)) + sum(model.b_charge[i] for i in range(1, 8761)) == 0
model.Summ = pyo.Constraint(rule=dispatch_rule)

def discharge_constraint_two(model, i):
    return model.b_discharge[i]  <= model.b_soc[i] # Maximum discharge rate within a single hour
model.discharge_two = pyo.Constraint(model.time_set, rule=discharge_constraint_two)

def charge_constraint_one(model, i):
    return model.b_charge[i] <= model.HP_max[i] - model.Demand[i] # Maximum charge rate within a single hour
model.charge_one = pyo.Constraint(model.time_set, rule=charge_constraint_one)

def charge_constraint_two(model, i):
    return model.b_charge[i] <= 8.708 - model.b_soc[i] # Maximum charge rate within a single hour
model.charge_two = pyo.Constraint(model.time_set, rule=charge_constraint_two)

def Sum_Const(model, i):
    return model.HP_var[i] + model.Boiler_var[i] + model.b_discharge[i] - model.b_charge[i] == model.Demand[i]
model.SumConst = pyo.Constraint(model.time_set, rule=Sum_Const)

def Summation_Const(model, i):
    return sum(model.HP_var[i] for i in range(1, 8761)) + sum(model.Boiler_var[i] for i in range(1, 8761)) == sum(model.Demand[i] for i in range(1, 8761))
model.Summ = pyo.Constraint(rule=Summation_Const)

def Emission_Const(model, i):
    total_gas = sum(model.Boiler_var[i] for i in range(1, 8761))
    total_elec = sum(model.HP_var[i] for i in range(1, 8761))
    total_demand = sum(model.Demand[i] for i in range(1, 8761))
    return ((total_elec / 2.97 * 0.219) + (total_gas / 0.85 * 0.184)) / (total_demand * 0.984) == 0.1
model.Emission = pyo.Constraint(rule=Emission_Const)'''

GLPSOL--GLPK LP/MIP Solver 5.0
Parameter(s) specified in the command line:
 --mipgap 0.015 --write C:\Users\MATE~1.TOT\AppData\Local\Temp\tmpyopw_8dk.glpk.raw
 --wglp C:\Users\MATE~1.TOT\AppData\Local\Temp\tmp9aorspj7.glpk.glp --cpxlp
 C:\Users\MATE~1.TOT\AppData\Local\Temp\tmpmlk9d5az.pyomo.lp
Reading problem data from 'C:\Users\MATE~1.TOT\AppData\Local\Temp\tmpmlk9d5az.pyomo.lp'...
17523 rows, 43800 columns, 96358 non-zeros
201490 lines were read
Writing problem data to 'C:\Users\MATE~1.TOT\AppData\Local\Temp\tmp9aorspj7.glpk.glp'...
201484 lines were written
GLPK Simplex Optimizer 5.0
17523 rows, 43800 columns, 96358 non-zeros
Preprocessing...
17521 rows, 35038 columns, 87595 non-zeros
Scaling...
 A: min|aij| =  1.576e-06  max|aij| =  1.000e+00  ratio =  6.345e+05
GM: min|aij| =  7.640e-01  max|aij| =  1.309e+00  ratio =  1.713e+00
EQ: min|aij| =  5.836e-01  max|aij| =  1.000e+00  ratio =  1.713e+00
Constructing initial basis...
Size of triangular part is 17520
      0: obj =   0.000000000e+00 inf =   6.458e+04 (4)
Perturbing LP to avoid stalling [11026]...
  13932: obj =   2.831720898e+06 inf =   9.314e+02 (1) 84
  17298: obj =   2.831720898e+06 inf =   9.314e+02 (1)
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
Time used:   6.5 secs
Memory used: 37.8 Mb (39629941 bytes)
Writing basic solution to 'C:\Users\MATE~1.TOT\AppData\Local\Temp\tmpyopw_8dk.glpk.raw'...
61332 lines were written
Model unknown
Toth Mate
  • 1

1 Answers1

0

Welcome to the site. There are a couple things you should look at.

First, there are errors in your model, so I'm not sure what you are running that actually produces output. Specifically, in several of your constraints, you are omitting the index that needs to be passed to your function. Like this one:

def Summation_Const(model, i):
    return sum(model.HP_var[i] for i in range(1, 8761)) + sum(model.Boiler_var[i] for i in range(1, 8761)) == sum(model.Demand[i] for i in range(1, 8761))
model.Summ = pyo.Constraint(rule=Summation_Const)

You aren't passing in i.

Second, I'm concerned about your Emission constraint. Why is it an equality constraint rather than an upper limit? (<=)

Most importantly, you need to focus on the result statement and see that it reports the model is currently INFEASIBLE. So it is not just that an optimal solution can't be found, it is that you have some logic errors. You can comment out the constraints individually to see if you can find the culprit or do something as I suggest in this post:

Finding out reason of Pyomo model infeasibility

AirSquid
  • 10,214
  • 2
  • 7
  • 31