I am working on optimizing a heat dispatch model for a District Heating (DH) network. In this model, I have different heat assets that can produce heat at varying prices. The challenge arises from the fact that the price of heat produced by the Combined Heat and Power (CHP) plant depends on its heat load.

I want to include this dependency in the optimization process. Specifically, I aim to constrain the CHP heat load (model.CHP[t]) and its corresponding price ($model.price_CHP[t]$) to discrete points. At each time step, I have predefined arrays containing different heat load points ($chp.Line_PQ$) and their corresponding prices ($chp.Price_PQ[0]$).

Here are the relevant variables and parameters from my Pyomo model:

Variablesyour text

    model.CHP = pyomo.Var(model.T, domain=pyomo.NonNegativeReals)
    model.EB = pyomo.Var(model.T, domain=pyomo.NonNegativeReals)
    model.HP = pyomo.Var(model.T, domain=pyomo.NonNegativeReals)
    model.BO = pyomo.Var(model.T, domain=pyomo.NonNegativeReals)
    model.price_CHP = pyomo.Var(model.T, domain=pyomo.NonNegativeReals)```

    model.CHP = pyomo.Var(model.T, domain=pyomo.NonNegativeReals)
    model.EB = pyomo.Var(model.T, domain=pyomo.NonNegativeReals)
    model.HP = pyomo.Var(model.T, domain=pyomo.NonNegativeReals)
    model.BO = pyomo.Var(model.T, domain=pyomo.NonNegativeReals)
    model.price_CHP = pyomo.Var(model.T, domain=pyomo.NonNegativeReals)

Parameters

    model.price_EB = pyomo.Param(model.T, initialize=Day_price_EB)
    model.price_HP = pyomo.Param(model.T, initialize=Day_price_HP)
    model.price_BO = pyomo.Param(model.T, initialize=Day_price_BO)

And here's the objective function I'm using:

    def objective_func(model):
        return sum([model.CHP[t]*model.price_CHP[t] + model.HP[t]*model.price_HP[t] + model.EB[t]*model.price_EB[t] +      model.BO[t]*model.price_BO[t] for t in model.T])

For the constraints, I've attempted to implement them as follows:

    discrete_points_chp = chp.Line_PQ
    discrete_prices_chp = chp.Price_PQ[0]

    def constrain_chp_to_discrete_points(model, t):
        return model.CHP[t] in set(discrete_points_chp)

    model.constr_chp_discrete = pyomo.Constraint(model.T, rule=constrain_chp_to_discrete_points)

    def constrain_price_chp_to_discrete_points(model, t):
    return model.price_CHP[t] in set(discrete_prices_chp)

However, when I run the optimization, I encounter issues with these constraints. Could someone please provide guidance on how to correctly implement such constraints in Pyomo?

SO I'm trying to use a set to define the discrete options for my variables but this does not work....

The constraint

 x ∈ {p[0],p[1],...p[n-1]}   (p[i] are constants)

can be implemented as:

 x = sum(p[i]*y[i])
 sum(y[i]) = 1
 y[i] ∈ {0,1}  (binary variable)