PYOMO优化问题（热存储＆＃x2B;锅炉＆＃x2B;热泵）

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。

其次，我担心您的排放限制。为什么它是平等约束而不是上限？ （＆lt; =）

最重要的是，您需要关注结果语句，并看到它报告了该模型当前是不可行的。因此，不仅找不到最佳解决方案，还因为您有一些逻辑错误。您可以单独评论约束，以查看您是否可以找到罪魁祸首或做我在这篇文章中所建议的事情：

找出PYOMO模型的原因