Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Pyomo is an open-source framework for formulating, solving and analysing optimisation problems with Python.

One can embed within Python an optimization model consisting of decision variables, constraints, and an optimisation objective, and solve these instances using commercial and open-source solvers (specialised software).

Pyomo supports a wide range of problem types, including:

  • Linear programming

  • Quadratic programming

  • Nonlinear programming

  • Mixed-integer linear programming

  • Mixed-integer quadratic programming

  • Mixed-integer nonlinear programming

  • Stochastic programming

  • Generalized disjunctive programming

  • Differential algebraic equations

  • Bilevel programming

  • Mathematical programs with equilibrium constraints

Basic Components

Almost all Pyomo objects exist within the pyomo.environ namespace:

import pyomo.environ as pe

A model instance is created by calling the pe.ConcreteModel() constructor and and assigninig it to a local variable. It’s good practice to choose a short variable name to reduce the verbosity of the code.

m = pe.ConcreteModel()
type(m)
pyomo.core.base.PyomoModel.ConcreteModel

Decision Variables

Variables are the unknowns of an optimisation problems and are intended to be given values by solving an optimisation problem.

In pyomo, decision variables are created with the pe.Var() constructor.

In the code below, we create two decision variables:

x,y∈R+x, y \in \mathbb{R}^+
(1)

  • The name you assign the object to (e.g. m.x) must be unique in any given model.

  • The pe.Var() constructor can take a variety of keyword arguments.

  • For instance, within=... would set the variable domain.

  • There are several pre-defined domains, like pe.NonNegativeReals.

  • But you can also manually define the bounds=... with a tuple representing lower and upper bounds.

m.x = pe.Var(within=pe.NonNegativeReals)
m.y = pe.Var(bounds=(0, None))

Objectives

An objective is a function of variables that returns a value that an optimization package attempts to maximize or minimize.

In pyomo, objectives are created with the pe.Objective() constructor.

In the code below, we create a simple objective function:

max⁡x,y4x+3y\max_{x,y} 4x + 3y
(2)

  • Like for variables, the object returned by pe.Objective() is assigned to the model under a unique name.

  • The keyword argument expr=... can be a mathematical expression or a function that returns an expression.

  • The keyword argument sense=... determines whether the objective should be maximised or minimised. You can conveniently access the optimisation sense through pe.minimize and pe.maximize. The default is pe.minimise.

m.obj = pe.Objective(expr=4 * m.x + 3 * m.y, sense=pe.maximize)

You can inspect almost all pyomo components with the .pprint() function:

m.obj.pprint()
obj : Size=1, Index=None, Active=True
    Key  : Active : Sense    : Expression
    None :   True : maximize : 4*x + 3*y

Constraints

Constraints are equality or inequality expressions that constrain the feasible space of the decision variables.

In pyomo, constraints are created with the pe.Constraint() constructor.

In the code below, we create three simple constraints in different ways:

x≤42x+y≤10x+y≤10\begin{align} x & \leq 4 \\ 2x + y & \leq 10 \\ x + y & \leq 10 \\ \end{align}
(3)

  • Again, the object return by pe.Constraint() is assigned to the model under a unique name.

  • Like for objectives, constraints take a mathematical expression under expr=....

  • Expressions can also be represented by a 3-tuple consisting of lower bound, body, upper bound.

m.c1 = pe.Constraint(expr=m.x <= 4)
m.c2 = pe.Constraint(expr=2 * m.x + m.y <= 10)
m.c3 = pe.Constraint(expr=(None, m.x + m.y, 10))

We can also call .pprint() on the whole model m to get a full representation of the optimisation problem built:

m.pprint()
2 Var Declarations
    x : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  None :  None : False :  True : NonNegativeReals
    y : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  None :  None : False :  True :  Reals

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : 4*x + 3*y

3 Constraint Declarations
    c1 : Size=1, Index=None, Active=True
        Key  : Lower : Body : Upper : Active
        None :  -Inf :    x :   4.0 :   True
    c2 : Size=1, Index=None, Active=True
        Key  : Lower : Body    : Upper : Active
        None :  -Inf : 2*x + y :  10.0 :   True
    c3 : Size=1, Index=None, Active=True
        Key  : Lower : Body  : Upper : Active
        None :  -Inf : x + y :  10.0 :   True

6 Declarations: x y obj c1 c2 c3

Solving optimisation problems

Solver installation

Solvers are needed to compute solutions to the optimization models developed using Pyomo.

There exists a large variety of solvers. In many cases, they specialise in certain problem types or solving algorithms, e.g. linear or nonlinear problems.

The open-source solvers are sufficient to handle meaningful Pyomo models with hundreds to several thousand variables and constraints. However, as applications get large or more complex, there may be a need to turn to a commercial solvers (which often provide free academic licenses).

Pyomo has support for a wide variety of solvers, but they need to be installed separately!

For this course, we use HiGHS, which is already in the course environment.

The open-source solvers are sufficient to handle meaningful Pyomo models with hundreds to several thousand variables and constraints. However, as applications get large or more complex, there may be a need to turn to a commercial solvers (with free academic licenses).

Calling solvers from pyomo

First, we create an instance of a solver (aka optimiser):

opt = pe.SolverFactory("appsi_highs")  # try also 'cbc', 'glpk', 'gurobi'

We can tell this optimiser to solve the model m:

log = opt.solve(m, tee=True)

The keyword argument tee=True is optional causes a more verbose display of solver log outputs.

And then display the solver outputs:

log.write()
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Lower bound: 30.0
  Upper bound: 30.0
  Number of objectives: 1
  Number of constraints: 0
  Number of variables: 0
  Sense: maximize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Termination message: TerminationCondition.optimal
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

While we can read from the message above that our problem was solved successfully, we can also formally check by accessing the reported status and comparing it to a pre-defined range of status codes provided in pe.SolverStatus:

assert log.solver.status == pe.SolverStatus.ok

Other states of solved optimisation models are:

pe.SolverStatus._member_names_
['ok', 'warning', 'error', 'aborted', 'unknown']

Retrieving optimisation results

Primal Values

Objective value:

m.obj()
30.0

Decision variables

m.x()
-0.0
m.y()
10.0

Dual values (aka shadow prices)

Access to dual values in scripts is similar to accessing primal variable values, except that dual values are not captured by default so additional directives are needed before optimization to signal that duals are desired.

To signal that duals are desired, declare a Suffix (another pyomo component) with the name “dual” on the model.

m.dual = pe.Suffix(direction=pe.Suffix.IMPORT)

Resolve:

pe.SolverFactory("appsi_highs").solve(m).write()
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Lower bound: 30.0
  Upper bound: 30.0
  Number of objectives: 1
  Number of constraints: 0
  Number of variables: 0
  Sense: maximize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Termination message: TerminationCondition.optimal
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

Then, we can access the shadow prices of individual constraints:

m.dual[m.c1]
0.0
m.dual[m.c2]
1.0

Indexed Components and Sets

Let’s first create a fresh model instance:

m = pe.ConcreteModel()

Almost all Pyomo components can be indexed.

While constructing them, all non-keyword arguments are assumed to represent indices.

The indices can be any iterable object.

Sets and Range Sets

In many cases, it will be useful to make use of another basic pyomo component: the set, which is constructed by pe.Set() or pe.RangeSet(). These objects are handled very similarly to variables, constraints and objectives.

The pe.RangeSet will create a numbered index, starting from the value 1.

m.A = pe.RangeSet(2)
m.A.pprint()
A : Dimen=1, Size=2, Bounds=(1, 2)
    Key  : Finite : Members
    None :   True :   [1:2]

With pe.Set() we can also pass more general indices:

B = ["wind", "solar"]
m.B = pe.Set(initialize=B)
m.B.pprint()
B : Size=1, Index=None, Ordered=Insertion
    Key  : Dimen : Domain : Size : Members
    None :     1 :    Any :    2 : {'wind', 'solar'}

Indexed Variables

We can then use these sets to create pyomo components, like variables:

m.x = pe.Var(m.A)
m.x.pprint()
x : Size=2, Index=A
    Key : Lower : Value : Upper : Fixed : Stale : Domain
      1 :  None :  None :  None : False :  True :  Reals
      2 :  None :  None :  None : False :  True :  Reals

There are no restrictions as to how many dimensions an indexed variable can have:

m.y = pe.Var(m.A, m.B)
m.y.pprint()
y : Size=4, Index=A*B
    Key          : Lower : Value : Upper : Fixed : Stale : Domain
    (1, 'solar') :  None :  None :  None : False :  True :  Reals
     (1, 'wind') :  None :  None :  None : False :  True :  Reals
    (2, 'solar') :  None :  None :  None : False :  True :  Reals
     (2, 'wind') :  None :  None :  None : False :  True :  Reals

It is not strictly necessary to use set objects to create indexed variables. We can use any iterable object.

m.z = pe.Var(B)

List comprehension on indexed variables

List comprehensions are quite useful to process indexed variables in a compact format:

m.c1 = pe.Constraint(expr=sum(m.z[b] for b in m.B) <= 10)
m.c1.pprint()
c1 : Size=1, Index=None, Active=True
    Key  : Lower : Body               : Upper : Active
    None :  -Inf : z[wind] + z[solar] :  10.0 :   True

Indexed Constraints

Pyomo used the concept of construction rules to specify indexed constraints, rather than expr=....

When building indexed constraints, particular indices are passed as argument to a rule (function) that returns an expression for each index.

def construction_rule(m, a):
    return sum(m.y[a, b] for b in m.B) == 1
m.c2 = pe.Constraint(m.A, rule=construction_rule)
m.c2.pprint()
c2 : Size=2, Index=A, Active=True
    Key : Lower : Body                   : Upper : Active
      1 :   1.0 : y[1,wind] + y[1,solar] :   1.0 :   True
      2 :   1.0 : y[2,wind] + y[2,solar] :   1.0 :   True

You can access individual constraints of a set of indexed constraints like this:

m.c2[2].pprint()
{Member of c2} : Size=2, Index=A, Active=True
    Key : Lower : Body                   : Upper : Active
      2 :   1.0 : y[2,wind] + y[2,solar] :   1.0 :   True

It is also possible to use an alternative decorator notation to reduce redundancy in rule definitions.

@m.Constraint(m.A)
def c3(m, a):
    return sum(m.y[a, b] for b in B) >= 1

Electricity Market Examples

Single bidding zone, single period

We want to minimise operational cost of an example electricity system representing South Africa subject to generator limits and meeting the load:

min⁡gs∑sosgs\min_{g_s} \sum_s o_s g_s
(4)

such that

gs≤Gsgs≥0∑sgs=d\begin{align} g_s &\leq G_s \\ g_s &\geq 0 \\ \sum_s g_s &= d \end{align}
(5)

We are given the following information on the South African electricity system:

Marginal costs in EUR/MWh

marginal_costs = {
    "Wind": 0,
    "Coal": 30,
    "Gas": 60,
    "Oil": 80,
}

Power plant capacities in MW

capacities = {"Coal": 35000, "Wind": 3000, "Gas": 8000, "Oil": 2000}

Inelastic demand in MW

load = 42000

First step is to initialise a new model, and since we want to know about the shadow prices, we additionally pass the respective Suffix.

m = pe.ConcreteModel()
m.dual = pe.Suffix(direction=pe.Suffix.IMPORT)

We create a set SS for the technologies:

m.S = pe.Set(initialize=capacities.keys())

And the decision variables for the generator dispatch:

m.g = pe.Var(m.S, within=pe.NonNegativeReals)

Our objective is to minimise total operational costs:

m.cost = pe.Objective(expr=sum(marginal_costs[s] * m.g[s] for s in m.S))
m.cost.pprint()
cost : Size=1, Index=None, Active=True
    Key  : Active : Sense    : Expression
    None :   True : minimize : 30*g[Coal] + 0*g[Wind] + 60*g[Gas] + 80*g[Oil]

subject to the capacity limits of the generators

@m.Constraint(m.S)
def generator_limit(m, s):
    return m.g[s] <= capacities[s]
m.generator_limit.pprint()
generator_limit : Size=4, Index=S, Active=True
    Key  : Lower : Body    : Upper   : Active
    Coal :  -Inf : g[Coal] : 35000.0 :   True
     Gas :  -Inf :  g[Gas] :  8000.0 :   True
     Oil :  -Inf :  g[Oil] :  2000.0 :   True
    Wind :  -Inf : g[Wind] :  3000.0 :   True

and meeting energy demand

m.energy_balance = pe.Constraint(expr=sum(m.g[s] for s in m.S) == load)
m.energy_balance.pprint()
energy_balance : Size=1, Index=None, Active=True
    Key  : Lower   : Body                                : Upper   : Active
    None : 42000.0 : g[Coal] + g[Wind] + g[Gas] + g[Oil] : 42000.0 :   True

Then, we solve the problem:

pe.SolverFactory("appsi_highs").solve(m).write()
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Lower bound: 1290000.0
  Upper bound: 1290000.0
  Number of objectives: 1
  Number of constraints: 0
  Number of variables: 0
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Termination message: TerminationCondition.optimal
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

Hopefully, the optimisation was successful (check!).

We can use pandas to retrieve and process the optimisation results.

import pandas as pd

This is the optimal generator dispatch (in MW):

pd.Series(m.g.get_values())
Coal 35000.0 Wind 3000.0 Gas 4000.0 Oil 0.0 dtype: float64

And the market clearing price we can read from the shadow price of the energy balance constraint (i.e. the added cost of increasing electricity demand by one unit):

market_price = m.dual[m.energy_balance]
market_price
60.0

There are further interesting shadow prices. The dual values of the generator limits tell us by how much the objective would be reduced if the capacity of the respective generator would be increased. It is either zero (cf. complementary slackness) for the generators for which their capacity limit is not binding (i.e. oil and gas), or, if the constraint is binding, it denotes the difference between the market clearing price and the marginal cost of the respective generator.

Retrieve all shadow prices at once:

pd.Series(m.dual.values(), m.dual.keys())
generator_limit[Coal] -30.0 generator_limit[Wind] -60.0 generator_limit[Gas] -0.0 generator_limit[Oil] -0.0 [None] 60.0 dtype: float64

A more targeted alternative for the generator limits:

pd.Series({s: m.dual[m.generator_limit[s]] for s in m.S})
Coal -30.0 Wind -60.0 Gas -0.0 Oil -0.0 dtype: float64

Two bidding zones with transmission

Let’s add a spatial dimension, such that the optimisation problem is expanded to

min⁡gi,s,fℓ∑soi,sgi,s\min_{g_{i,s}, f_\ell} \sum_s o_{i,s} g_{i,s}
(6)

such that

gi,s≤Gi,sgi,s≥0∑sgi,s−∑ℓKiℓfℓ=diKCL∣fℓ∣≤Fℓline limits∑ℓCℓcxℓfℓ=0KVL\begin{align} g_{i,s} &\leq G_{i,s} \\ g_{i,s} &\geq 0 \\ \sum_s g_{i,s} - \sum_\ell K_{i\ell} f_\ell &= d_i & \text{KCL} \\ |f_\ell| &\leq F_\ell & \text{line limits} \\ \sum_\ell C_{\ell c} x_\ell f_\ell &= 0 & \text{KVL} \end{align}
(7)

In this example, we connect the previous South African electricity system with a hydro generation unit in Mozambique through a single transmission line. Note that because a single transmission line will not result in any cycles, we can neglect KVL in this case.

We are given the following data (all in MW):

capacities = {
    "South Africa": {"Coal": 35000, "Wind": 3000, "Gas": 8000, "Oil": 2000},
    "Mozambique": {"Hydro": 1200},
}
transmission = 500
loads = {"South Africa": 42000, "Mozambique": 650}
marginal_costs["Hydro"] = 0  # MWh

First, let’s create a new model instance:

m = pe.ConcreteModel()
m.dual = pe.Suffix(direction=pe.Suffix.IMPORT)

We now have two sets: one for the countries, one for the technologies, based on which we can create the generator dispatch variables.

m.countries = pe.Set(initialize=loads.keys())
technologies = list(capacities["South Africa"].keys() | capacities["Mozambique"])
m.technologies = pe.Set(initialize=technologies)
m.g = pe.Var(m.countries, m.technologies, within=pe.NonNegativeReals)

We als need an additional variable for the flow.

m.f = pe.Var()

The objective remains unchanged:

m.cost = pe.Objective(
    expr=sum(marginal_costs[s] * m.g[c, s] for s in m.technologies for c in m.countries)
)

Generator limits also remain quite similar albeit for the additional index. For missing entries, we set the upper limit to zero (e.g. this country-technology combination can not produce electricity).

@m.Constraint(m.countries, m.technologies)
def generator_limit(m, c, s):
    return m.g[c, s] <= capacities[c].get(s, 0)
m.generator_limit.pprint()
generator_limit : Size=10, Index=countries*technologies, Active=True
    Key                       : Lower : Body                  : Upper   : Active
       ('Mozambique', 'Coal') :  -Inf :    g[Mozambique,Coal] :     0.0 :   True
        ('Mozambique', 'Gas') :  -Inf :     g[Mozambique,Gas] :     0.0 :   True
      ('Mozambique', 'Hydro') :  -Inf :   g[Mozambique,Hydro] :  1200.0 :   True
        ('Mozambique', 'Oil') :  -Inf :     g[Mozambique,Oil] :     0.0 :   True
       ('Mozambique', 'Wind') :  -Inf :    g[Mozambique,Wind] :     0.0 :   True
     ('South Africa', 'Coal') :  -Inf :  g[South Africa,Coal] : 35000.0 :   True
      ('South Africa', 'Gas') :  -Inf :   g[South Africa,Gas] :  8000.0 :   True
    ('South Africa', 'Hydro') :  -Inf : g[South Africa,Hydro] :     0.0 :   True
      ('South Africa', 'Oil') :  -Inf :   g[South Africa,Oil] :  2000.0 :   True
     ('South Africa', 'Wind') :  -Inf :  g[South Africa,Wind] :  3000.0 :   True

The energy balance constraint is replaced by KCL, where we take into account local generation as well as incoming or outgoing flows.

We also need the incidence matrix of this network (here it’s very simple!) and assume some direction for the flow variable. Here, we picked the orientation from South Africa to Mozambique. This means that if the values for the flow variable m.f are positive South Africa exports to Mozambique and vice versa if the variable takes negative values.

@m.Constraint(m.countries)
def kcl(m, c):
    sign = -1 if c == "Mozambique" else 1  # minimal incidence matrix
    return sum(m.g[c, s] for s in m.technologies) - sign * m.f == loads[c]
m.kcl.pprint()
kcl : Size=2, Index=countries, Active=True
    Key          : Lower   : Body                                                                                                                : Upper   : Active
      Mozambique :   650.0 :           g[Mozambique,Hydro] + g[Mozambique,Coal] + g[Mozambique,Gas] + g[Mozambique,Wind] + g[Mozambique,Oil] + f :   650.0 :   True
    South Africa : 42000.0 : g[South Africa,Hydro] + g[South Africa,Coal] + g[South Africa,Gas] + g[South Africa,Wind] + g[South Africa,Oil] - f : 42000.0 :   True

We also need to constrain the transmission flows to the line’s rated capacity. Here, we use a 3-tuple to specify upper and lower bounds in one go.

m.line_limit = pe.Constraint(expr=(-transmission, m.f, transmission))
m.line_limit.pprint()
line_limit : Size=1, Index=None, Active=True
    Key  : Lower  : Body : Upper : Active
    None : -500.0 :    f : 500.0 :   True

Then, we can solve and inspect results:

pe.SolverFactory("appsi_highs").solve(m).write()
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Lower bound: 1260000.0
  Upper bound: 1260000.0
  Number of objectives: 1
  Number of constraints: 0
  Number of variables: 0
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Termination message: TerminationCondition.optimal
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

m.cost()
1260000.0
pd.Series(m.g.get_values()).unstack()
Loading...
m.f()
-500.0
pd.Series(m.dual.values(), m.dual.keys())
generator_limit[South Africa,Hydro] -60.0 generator_limit[South Africa,Coal] -30.0 generator_limit[South Africa,Gas] -0.0 generator_limit[South Africa,Wind] -60.0 generator_limit[South Africa,Oil] -0.0 generator_limit[Mozambique,Hydro] -0.0 generator_limit[Mozambique,Coal] -0.0 generator_limit[Mozambique,Gas] -0.0 generator_limit[Mozambique,Wind] -0.0 generator_limit[Mozambique,Oil] -0.0 kcl[South Africa] 60.0 kcl[Mozambique] -0.0 [None] 60.0 dtype: float64

In pyomo it is also possible to deactivate (and activate) constraints, and fix variables to a pre-defined value, which alters the results

m.line_limit.deactivate()
m.g["South Africa", "Coal"].fix(34000)
pe.SolverFactory("appsi_highs").solve(m).write()
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Lower bound: 1287000.0
  Upper bound: 1287000.0
  Number of objectives: 1
  Number of constraints: 0
  Number of variables: 0
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Termination message: TerminationCondition.optimal
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

pd.Series(m.g.get_values()).unstack()
Loading...
m.cost()
1287000.0

Single bidding zone with several periods

In this example, we consider multiple time periods (labelled [1,2,3,4]) to represent variable wind generation and changing load.

min⁡gs,t∑sosgs,t\min_{g_{s,t}} \sum_s o_{s} g_{s,t}
(8)

such that

gs,t≤g^s,tGi,sgs,t≥0∑sgs,t=dt\begin{align} g_{s,t} &\leq \hat{g}_{s,t} G_{i,s} \\ g_{s,t} &\geq 0 \\ \sum_s g_{s,t} &= d_t \end{align}
(9)
wind_ts = [0.3, 0.6, 0.4, 0.5]
load_ts = [42000, 43000, 45000, 46000]
m = pe.ConcreteModel()
m.dual = pe.Suffix(direction=pe.Suffix.IMPORT)
techs = capacities["South Africa"].keys()

m.S = pe.Set(initialize=techs)
m.T = pe.RangeSet(4)
m.g = pe.Var(m.S, m.T, within=pe.NonNegativeReals)
m.cost = pe.Objective(expr=sum(marginal_costs[s] * m.g[s, t] for s in m.S for t in m.T))
@m.Constraint(m.S, m.T)
def generator_limit(m, s, t):
    cf = wind_ts[t - 1] if s == "Wind" else 1
    return m.g[s, t] <= cf * capacities["South Africa"][s]
m.generator_limit.pprint()
generator_limit : Size=16, Index=S*T, Active=True
    Key         : Lower : Body      : Upper   : Active
    ('Coal', 1) :  -Inf : g[Coal,1] : 35000.0 :   True
    ('Coal', 2) :  -Inf : g[Coal,2] : 35000.0 :   True
    ('Coal', 3) :  -Inf : g[Coal,3] : 35000.0 :   True
    ('Coal', 4) :  -Inf : g[Coal,4] : 35000.0 :   True
     ('Gas', 1) :  -Inf :  g[Gas,1] :  8000.0 :   True
     ('Gas', 2) :  -Inf :  g[Gas,2] :  8000.0 :   True
     ('Gas', 3) :  -Inf :  g[Gas,3] :  8000.0 :   True
     ('Gas', 4) :  -Inf :  g[Gas,4] :  8000.0 :   True
     ('Oil', 1) :  -Inf :  g[Oil,1] :  2000.0 :   True
     ('Oil', 2) :  -Inf :  g[Oil,2] :  2000.0 :   True
     ('Oil', 3) :  -Inf :  g[Oil,3] :  2000.0 :   True
     ('Oil', 4) :  -Inf :  g[Oil,4] :  2000.0 :   True
    ('Wind', 1) :  -Inf : g[Wind,1] :   900.0 :   True
    ('Wind', 2) :  -Inf : g[Wind,2] :  1800.0 :   True
    ('Wind', 3) :  -Inf : g[Wind,3] :  1200.0 :   True
    ('Wind', 4) :  -Inf : g[Wind,4] :  1500.0 :   True

@m.Constraint(m.T)
def energy_balance(m, t):
    return sum(m.g[s, t] for s in m.S) == load_ts[t - 1]
m.energy_balance.pprint()
energy_balance : Size=4, Index=T, Active=True
    Key : Lower   : Body                                        : Upper   : Active
      1 : 42000.0 : g[Coal,1] + g[Wind,1] + g[Gas,1] + g[Oil,1] : 42000.0 :   True
      2 : 43000.0 : g[Coal,2] + g[Wind,2] + g[Gas,2] + g[Oil,2] : 43000.0 :   True
      3 : 45000.0 : g[Coal,3] + g[Wind,3] + g[Gas,3] + g[Oil,3] : 45000.0 :   True
      4 : 46000.0 : g[Coal,4] + g[Wind,4] + g[Gas,4] + g[Oil,4] : 46000.0 :   True

Now, we can solve. Let’s also time how long it takes to solve this problem with a small IPython tool.

%%timeit -n 1 -r 1
pe.SolverFactory("appsi_highs").solve(m).write()
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Lower bound: 6082000.0
  Upper bound: 6082000.0
  Number of objectives: 1
  Number of constraints: 0
  Number of variables: 0
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Termination message: TerminationCondition.optimal
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
7.42 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)

m.cost()
6082000.0
pd.Series(m.dual.values(), m.dual.keys())
generator_limit[Coal,1] -30.0 generator_limit[Coal,2] -30.0 generator_limit[Coal,3] -50.0 generator_limit[Coal,4] -50.0 generator_limit[Wind,1] -60.0 generator_limit[Wind,2] -60.0 generator_limit[Wind,3] -80.0 generator_limit[Wind,4] -80.0 generator_limit[Gas,1] -0.0 generator_limit[Gas,2] -0.0 generator_limit[Gas,3] -20.0 generator_limit[Gas,4] -20.0 generator_limit[Oil,1] -0.0 generator_limit[Oil,2] -0.0 generator_limit[Oil,3] -0.0 generator_limit[Oil,4] -0.0 energy_balance[1] 60.0 energy_balance[2] 60.0 energy_balance[3] 80.0 energy_balance[4] 80.0 dtype: float64
pd.Series(m.g.get_values()).unstack(0)
Loading...

Single bidding zone with several periods and storage

Now, we want to expand the optimisation model with a storage unit to do price arbitrage to reduce oil consumption.

We have been given the following characteristics of the storage:

storage_energy = 6000  # MWh
storage_power = 1000  # MW
efficiency = 0.9  # discharge = charge
standing_loss = 0.00001  # per hour

To model a storage unit, we need three additional variables for the discharging and charging of the storage unit and for its state of charge (energy filling level). We can directly define the bounds of these variables in the variable definition:

m.discharge = pe.Var(m.T, bounds=(0, storage_power))
m.charge = pe.Var(m.T, bounds=(0, storage_power))
m.soc = pe.Var(m.T, bounds=(0, storage_energy))

Then, we implement the storage consistency equations,

et=(1−standing loss)⋅et−1+η⋅gcharge,t−1η⋅gdischarge,te_{t} = (1-\text{standing loss}) \cdot e_{t-1} + \eta \cdot g_{charge, t} - \frac{1}{\eta} \cdot g_{discharge, t}
(10)

For the initial period, we set the state of charge to zero.

e0=0e_{0} = 0
(11)
@m.Constraint(m.T)
def storage_consistency(m, t):
    if t == 1:
        return m.soc[t] == 0
    return (
        m.soc[t]
        == (1 - standing_loss) * m.soc[t - 1]
        + efficiency * m.charge[t]
        - 1 / efficiency * m.discharge[t]
    )
m.storage_consistency.pprint()
storage_consistency : Size=4, Index=T, Active=True
    Key : Lower : Body                                                                        : Upper : Active
      1 :   0.0 :                                                                      soc[1] :   0.0 :   True
      2 :   0.0 : soc[2] - (0.99999*soc[1] + 0.9*charge[2] - 1.1111111111111112*discharge[2]) :   0.0 :   True
      3 :   0.0 : soc[3] - (0.99999*soc[2] + 0.9*charge[3] - 1.1111111111111112*discharge[3]) :   0.0 :   True
      4 :   0.0 : soc[4] - (0.99999*soc[3] + 0.9*charge[4] - 1.1111111111111112*discharge[4]) :   0.0 :   True

And we also need to modify the energy balance to include the contributions of storage discharging and charging.

For that, we should first remove the existing energy balance constraint, which we seek to overwrite.

m.del_component(m.energy_balance)
@m.Constraint(m.T)
def energy_balance(m, t):
    return sum(m.g[s, t] for s in m.S) + m.discharge[t] - m.charge[t] == load_ts[t - 1]
pe.SolverFactory("appsi_highs").solve(m).write()
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Lower bound: 6017200.65599352
  Upper bound: 6017200.65599352
  Number of objectives: 1
  Number of constraints: 0
  Number of variables: 0
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Termination message: TerminationCondition.optimal
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

m.cost()
6017200.65599352
pd.Series(m.g.get_values()).unstack(0)
Loading...
pd.Series(m.charge.get_values())
1 0.0 2 1000.0 3 0.0 4 0.0 dtype: float64
pd.Series(m.discharge.get_values())
1 1000.0000 2 0.0000 3 800.0000 4 9.9918 dtype: float64
pd.Series(m.soc.get_values())
1 -0.000000 2 900.000000 3 11.102111 4 0.000000 dtype: float64
pd.Series(m.dual.values(), m.dual.keys())
generator_limit[Coal,1] -30.00000 generator_limit[Coal,2] -30.00000 generator_limit[Coal,3] -49.99920 generator_limit[Coal,4] -50.00000 generator_limit[Wind,1] -60.00000 generator_limit[Wind,2] -60.00000 generator_limit[Wind,3] -79.99920 generator_limit[Wind,4] -80.00000 generator_limit[Gas,1] -0.00000 generator_limit[Gas,2] -0.00000 generator_limit[Gas,3] -19.99920 generator_limit[Gas,4] -20.00000 generator_limit[Oil,1] -0.00000 generator_limit[Oil,2] -0.00000 generator_limit[Oil,3] -0.00000 generator_limit[Oil,4] -0.00000 [Unattached ConstraintData] 60.00000 [Unattached ConstraintData] 60.00000 [Unattached ConstraintData] 80.00000 [Unattached ConstraintData] 80.00000 storage_consistency[1] -71.99784 storage_consistency[2] -71.99856 storage_consistency[3] -71.99928 storage_consistency[4] -72.00000 energy_balance[1] 60.00000 energy_balance[2] 60.00000 energy_balance[3] 79.99920 energy_balance[4] 80.00000 dtype: float64

Exercise

  • Using the conversion efficiencies and specific emissions from the lecture slides, add a constraint that limits the total emissions in the four periods to 50% of the unconstrained optimal solution. How does the optimal objective value and the generator dispatch change?

  • Reimplement the storage consistency constraint such that the initial state of charge is not zero but corresponds to the state of charge in the final period of the optimisation horizon.

  • What parameters of the storage unit would have to be changed to reduce the objective? What’s the sensitivity?

References
  1. Bynum, M. L., Hackebeil, G. A., Hart, W. E., Laird, C. D., Nicholson, B. L., Siirola, J. D., Watson, J.-P., & Woodruff, D. L. (2021). Pyomo — Optimization Modeling in Python. In Springer Optimization and Its Applications. Springer International Publishing. 10.1007/978-3-030-68928-5
Data Science for Energy System Modelling
Plotly and Streamlit
Solutions
Solutions Python