Issue getting decision variable to optimize properly in Pyomo

Question

I have a simple multi-period optimization problem I am working on using pyomo. The goal of the model is to determine which hours a power plant should be on or off based on the Spark Spread (Power Price - Gas Price * Heat Rate + Variable Costs) for that hour. The Spark spread can be negative, which would indicate that the plant should be off, or positive which means the plant should be running.

Currently the results are a showing that the plant just gets turned on and runs despite the spark spread being negative.

How can I get the plant to turn on and off at each time step, given the spark spread for that hour?

I am sure this is a fairly simple solution, but I am very new to pyomo and optimization problems, so any guidance and help would be much appreciated.

gas_price = [2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81]
power_price = [26.24,23.8,21.94,20.4,21.2,19.98,19.34,18.83,19.19,18.48,21,21.77,23.45,26.53,29.85,31.8,28.7]

priceDict = dict(enumerate(power_price))
gasDict = dict(enumerate(gas_price))

m = en.ConcreteModel()

m.Time = en.RangeSet(0, len(power_price)-1)
m.powerPrice = en.Param(m.Time, initialize=priceDict)
m.gasPrice = en.Param(m.Time, initialize=gasDict)

m.generation = en.Var(m.Time, bounds=(0,800),
initialize=0)
m.spark = en.Var(m.Time,initialize=0)
m.heatRate = en.Var(m.Time,initialize=7)
m.vom = en.Var(m.Time,initialize=2)

m.max_gen = en.Param(initialize=800)


def Obj_fn(m):
    return sum((m.spark[i]*m.generation[i]) for i in m.Time)
m.total_cost = en.Objective(rule=Obj_fn,sense=en.maximize)

# 7 is the heat rate of the plant
def spark_rule(m,i):
    return (m.spark[i] == m.powerPrice[i]-(m.gasPrice[i]*7+m.vom[i]))
m.hourly_spark = en.Constraint(m.Time,rule=spark_rule)

def generation_rule(m,i):
    return (0<=m.generation[i]<=m.max_gen)
m.t_generation_rule = en.Constraint(m.Time, rule=generation_rule)


opt = SolverFactory("clp",executable='C:\\clp.exe')
results = opt.solve(m)

The output of the model is currently:

Time Generation     Spark Spread
1   0               6.57
2   800             4.13
3   800             2.27
4   800             0.73
5   800             1.53
6   800             0.31
7   800            -0.33
8   800            -0.84
9   800            -0.48
10  800            -1.19
11  800             1.33
12  800             2.1
13  800             3.78
14  800             6.86
15  800            10.18
16  800            12.13
17  800             9.03

Answer 1

I might be wrong here, but I think you actually meant to define heatRate and vom as parameters, rather than variables.

This leads to a weird problem because the plant will naturally use its max power whenever the "spark" price is positive, and go to 0 whenever the spark price is negative. I guess you'll add more constraints later on.

If heatRate and vom are fixed, then you can redefine the problem in the following way:

from pyomo import environ as pe

gas_price = [2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81,2.81]
power_price = [26.24,23.8,21.94,20.4,21.2,19.98,19.34,18.83,19.19,18.48,21,21.77,23.45,26.53,29.85,31.8,28.7]

priceDict = dict(enumerate(power_price))
gasDict = dict(enumerate(gas_price))

m = pe.ConcreteModel()

m.Time = pe.RangeSet(0, len(power_price)-1)

# this are all input parameters
m.powerPrice = pe.Param(m.Time, initialize=priceDict)
m.gasPrice = pe.Param(m.Time, initialize=gasDict)
m.vom = pe.Param(default=7)
m.heatRate = pe.Param(default=2)
m.maxGen = pe.Param(default=800)

# this is a "dependent" parameter
m.spark = pe.Param(m.Time,
    initialize = lambda m,t: m.powerPrice[t]-(m.gasPrice[t]*7+m.vom)
)

# this is the only variable
m.generation = pe.Var(m.Time, 
    initialize=0,
    bounds = (0, m.maxGen)
)

def Obj_fn(m):
    return sum((m.spark[t]*m.generation[t]) for t in m.Time)

m.total_cost = pe.Objective(rule=Obj_fn,sense=pe.maximize)