Pyomo优化问题(蓄热+锅炉+热泵)

Question

你好，我已经编写了一个Pyomo脚本来优化一个包含3个组件(蓄热、锅炉和热泵)的调度模式。它似乎起作用，直到定义了有关蓄热的约束。有谁能帮助我理解为什么没有最优解，以及哪些约束(不是排放约束)被定义错了？谢谢你的帮助！

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

Answer 1

欢迎来到现场。有几件事你应该看看。

首先，您的模型中有错误，所以我不确定您正在运行的是什么，这实际上会产生输出。具体来说，在几个约束中，您忽略了需要传递给函数的索引。就像这个：

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)

你不能通过i。

第二，我担心你的排放限制。为什么它是一个等式约束而不是上限？(<=)

最重要的是，您需要关注结果语句，并确保它报告模型当前是不可行的。所以，不只是找不到最优解，而是有一些逻辑错误。你可以逐个评论这些限制，看看你是否能找到罪魁祸首，或者像我在这篇文章中建议的那样做些什么：

Finding out reason of Pyomo model infeasibility