Introduction to pyomo#
Note
This material is in part adapted from the following resources:
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
Note
Documentation for this package is available at http://www.pyomo.org/.
Note
If you have not yet set up Python on your computer, you can execute this tutorial in your browser via Google Colab. Click on the rocket in the top right corner and launch “Colab”. If that doesn’t work download the .ipynb file and import it in Google Colab.
Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.
!pip install pandas pyomo highspy
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
Note
Pyomo distinguishes between concrete and abstract models. The difference is that concrete models have real-valued input parameters in the building stage, whereas abstract models are initially build with placeholders for the input parameters.
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:
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:
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 throughpe.minimizeandpe.maximize. The default ispe.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:
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 can be installed easily via pip and is already in the course environment:
pip install highspy
If you want to try the other solvers, follow the links to their documentation for installation instructions.
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).
Note
For most solvers, it is possible to pass additional parameters before running the solver to configure things like solver tolerances, number of threads, type of solving algorithm (e.g. barrier or simplex). Read more about it in the Pyomo Documentation.
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)
Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
WARNING: No semi-integer/integer variables in model with non-empty integrality
Presolving model
2 rows, 2 cols, 4 nonzeros
2 rows, 2 cols, 4 nonzeros
Presolve : Reductions: rows 2(-1); columns 2(-0); elements 4(-1)
Solving the presolved LP
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
0 -2.9999888912e+00 Ph1: 2(2); Du: 1(2.99999) 0s
2 -3.0000000000e+01 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model status : Optimal
Simplex iterations: 2
Objective value : 3.0000000000e+01
HiGHS run time : 0.00
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: -1
# ----------------------------------------------------------
# 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()
Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
# ==========================================================
# = 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: -1
# ----------------------------------------------------------
# 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=y_index
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:
such that
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 \(S\) 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()
Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
# ==========================================================
# = 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: 1
# ----------------------------------------------------------
# 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
such that
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)
Note
Note that we could have been more efficient by only creating variables for the combinations of country and technology that do exist, rather than creating variables for all combinations.
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=generator_limit_index, 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,Oil] + g[Mozambique,Wind] + g[Mozambique,Coal] + g[Mozambique,Gas] + g[Mozambique,Hydro] + f : 650.0 : True
South Africa : 42000.0 : g[South Africa,Oil] + g[South Africa,Wind] + g[South Africa,Coal] + g[South Africa,Gas] + g[South Africa,Hydro] - 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()
Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
# ==========================================================
# = 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: 1
# ----------------------------------------------------------
# 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()
| Coal | Gas | Hydro | Oil | Wind | |
|---|---|---|---|---|---|
| Mozambique | 0.0 | 0.0 | 1150.0 | 0.0 | 0.0 |
| South Africa | 35000.0 | 3500.0 | -0.0 | 0.0 | 3000.0 |
m.f()
-500.0
pd.Series(m.dual.values(), m.dual.keys())
generator_limit[South Africa,Oil] -0.0
generator_limit[South Africa,Wind] -60.0
generator_limit[South Africa,Coal] -30.0
generator_limit[South Africa,Gas] -0.0
generator_limit[South Africa,Hydro] -60.0
generator_limit[Mozambique,Oil] -0.0
generator_limit[Mozambique,Wind] -0.0
generator_limit[Mozambique,Coal] -0.0
generator_limit[Mozambique,Gas] -0.0
generator_limit[Mozambique,Hydro] -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()
Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
# ==========================================================
# = 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: 1
# ----------------------------------------------------------
# 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()
| Coal | Gas | Hydro | Oil | Wind | |
|---|---|---|---|---|---|
| Mozambique | -0.0 | 0.0 | 1200.0 | 0.0 | -0.0 |
| South Africa | 34000.0 | 4450.0 | -0.0 | 0.0 | 3000.0 |
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.
such that
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=generator_limit_index, 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()
Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
# ==========================================================
# = 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: 1
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Termination condition: optimal
Termination message: TerminationCondition.optimal
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
3.8 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)
| Coal | Gas | Oil | Wind | |
|---|---|---|---|---|
| 1 | 35000.0 | 6100.0 | 0.0 | 900.0 |
| 2 | 35000.0 | 6200.0 | 0.0 | 1800.0 |
| 3 | 35000.0 | 8000.0 | 800.0 | 1200.0 |
| 4 | 35000.0 | 8000.0 | 1500.0 | 1500.0 |
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,
For the initial period, we set the state of charge to zero.
@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()
Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]
Copyright (c) 2023 HiGHS under MIT licence terms
# ==========================================================
# = 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: 1
# ----------------------------------------------------------
# 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)
| Coal | Gas | Oil | Wind | |
|---|---|---|---|---|
| 1 | 35000.0 | 5100.0 | 0.0000 | 900.0 |
| 2 | 35000.0 | 7200.0 | 0.0000 | 1800.0 |
| 3 | 35000.0 | 8000.0 | 0.0000 | 1200.0 |
| 4 | 35000.0 | 8000.0 | 1490.0082 | 1500.0 |
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 _GeneralConstraintData] 60.00000
[Unattached _GeneralConstraintData] 60.00000
[Unattached _GeneralConstraintData] 80.00000
[Unattached _GeneralConstraintData] 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?