Electricity Markets#
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 -q pypsa highspy
What is PyPSA?#
PyPSA stands for Python for Power System Analysis. The documentation describes it as follows:
PyPSA is an open-source Python framework for optimising and simulating modern power and energy systems that include features such as conventional generators with unit commitment, variable wind and solar generation, hydro-electricity, inter-temporal storage, coupling to other energy sectors, elastic demands, and linearised power flow with loss approximations in DC and AC networks. PyPSA is designed to scale well with large networks and long time series. It is made for researchers, planners and utilities with basic coding aptitude who need a fast, easy-to-use and transparent tool for power and energy system analysis.
For us, it’s a tool to model electricity markets and investment planning as optimisation problems as discussed in the slides.
PyPSA is a way to describe energy systems in Python code, interface with solver software, and analyse and visualise the results.
Note
For a full reference to the optimisation problem description, see the PyPSA User Guide
Input Data#
Let’s dive straight in and get acquainted with PyPSA by building simple electricity market models.
We have the following data:
fuel costs in € / MWh\(_{th}\)
fuel_cost = dict(
coal=8,
gas=100,
oil=48,
hydro=1,
wind=0,
)
efficiencies of thermal power plants in MWh\(_{el}\) / MWh\(_{th}\)
efficiency = dict(
coal=0.33,
gas=0.58,
oil=0.35,
hydro=1,
wind=1,
)
specific emissions in t\(_{CO_2}\) / MWh\(_{th}\)
# t/MWh thermal
emissions = dict(
coal=0.34,
gas=0.2,
oil=0.26,
hydro=0,
wind=0,
)
power plant capacities in MW
power_plants = {
"A": {"coal": 35000, "wind": 3000, "gas": 8000, "oil": 2000},
"B": {"hydro": 1200},
}
fixed electrical load in MW
loads = {
"A": 42000,
"B": 650,
}
Building the Network#
By convention, PyPSA is imported without an alias:
import pypsa
First, we create a new network object which serves as the overall container for all components.
n = pypsa.Network()
The second component we need are buses. Buses are the fundamental nodes of the network, to which all other components like loads, generators and transmission lines attach. They enforce energy conservation for all elements feeding in and out of it.
Components can be added to the network n using the n.add() function. It takes the component name as a first argument, the name of the component as a second argument and possibly further parameters as keyword arguments. Let’s use this function, to add buses for each country to our network:
n.add("Bus", "A", v_nom=380, carrier="AC")
n.add("Bus", "B", v_nom=380, carrier="AC");
For each class of components, the data describing the components is stored in a pandas.DataFrame. For example, all static data for buses is stored in n.buses
n.buses
| v_nom | type | x | y | carrier | unit | location | v_mag_pu_set | v_mag_pu_min | v_mag_pu_max | control | generator | sub_network | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Bus | |||||||||||||
| A | 380.0 | 0.0 | 0.0 | AC | 1.0 | 0.0 | inf | PQ | |||||
| B | 380.0 | 0.0 | 0.0 | AC | 1.0 | 0.0 | inf | PQ |
You see there are many more attributes than we specified while adding the buses; many of them are filled with default parameters which were added. You can look up the field description, defaults and status (required input, optional input, output) for buses here https://pypsa–1250.org.readthedocs.build/en/1250/user-guide/components/buses/, and analogous for all other components.
The method n.add() also allows you to add multiple components at once. For instance, multiple carriers for the fuels with information on specific carbon dioxide emissions, a nice name, and colors for plotting. For this, the function takes the component name as the first argument and then a list of component names and then optional arguments for the parameters. Here, scalar values, lists, dictionary or pandas.Series are allowed. The latter two needs keys or indices with the component names.
n.add(
"Carrier",
["coal", "gas", "oil", "hydro", "wind"],
co2_emissions=emissions,
nice_name=["Coal", "Gas", "Oil", "Hydro", "Onshore Wind"],
color=["grey", "indianred", "black", "aquamarine", "dodgerblue"],
);
The n.add() function is very general. It lets you add any component to the network object n. For instance, in the next step we add generators for all the different power plants.
In Region B:
n.add(
"Generator",
"B hydro",
bus="B",
carrier="hydro",
p_nom=1200, # MW
);
In Region A (in a loop):
for tech, p_nom in power_plants["A"].items():
n.add(
"Generator",
f"A {tech}",
bus="A",
carrier=tech,
efficiency=efficiency.get(tech, 1),
p_nom=p_nom,
marginal_cost=fuel_cost.get(tech, 0) / efficiency.get(tech, 1),
)
As a result, the n.generators DataFrame looks like this:
n.generators["marginal_cost"]
Generator
B hydro 0.000000
A coal 24.242424
A wind 0.000000
A gas 172.413793
A oil 137.142857
Name: marginal_cost, dtype: float64
Next, we’re going to add the electricity demand.
A positive value for p_set means consumption of power from the bus.
n.add(
"Load",
"A electricity demand",
bus="A",
p_set=loads["A"],
);
n.add(
"Load",
"B electricity demand",
bus="B",
p_set=loads["B"],
);
n.loads
| bus | carrier | type | p_set | q_set | sign | active | |
|---|---|---|---|---|---|---|---|
| Load | |||||||
| A electricity demand | A | 42000.0 | 0.0 | -1.0 | True | ||
| B electricity demand | B | 650.0 | 0.0 | -1.0 | True |
Finally, we add the connection between Region B and Region A with a 500 MW transmission line, modelled as a link:
n.add(
"Link",
"A-B",
bus0="A",
bus1="B",
p_min_pu=-1,
p_nom=500,
);
n.links
| bus0 | bus1 | type | carrier | efficiency | active | build_year | lifetime | p_nom | p_nom_mod | ... | shut_down_cost | min_up_time | min_down_time | up_time_before | down_time_before | ramp_limit_up | ramp_limit_down | ramp_limit_start_up | ramp_limit_shut_down | p_nom_opt | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Link | |||||||||||||||||||||
| A-B | A | B | 1.0 | True | 0 | inf | 500.0 | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 |
1 rows × 34 columns
Optimisation#
With all input data transferred into PyPSA’s data structure, we can now build and run the resulting optimisation problem. In PyPSA, building, solving and retrieving results from the optimisation model is contained in a single function call n.optimize(). This function optimizes dispatch and investment decisions for least cost.
The n.optimize() function can take a variety of arguments. The most relevant for the moment is the choice of the solver (e.g. “highs” and “gurobi”). They need to be installed on your computer, to use them here!
n.optimize(log_to_console=False)
Let’s have a look at the results.
Since the power flow and dispatch are generally time-varying quantities, these are stored in a different location than e.g. n.generators. They are stored in n.generators_t. Thus, to find out the dispatch of the generators, run
n.generators_t.p
| Generator | B hydro | A coal | A wind | A gas | A oil |
|---|---|---|---|---|---|
| snapshot | |||||
| now | 1150.0 | 35000.0 | 3000.0 | 1500.0 | 2000.0 |
or if you prefer it in relation to the generators nominal capacity
n.generators_t.p / n.generators.p_nom
| Generator | B hydro | A coal | A wind | A gas | A oil |
|---|---|---|---|---|---|
| snapshot | |||||
| now | 0.958333 | 1.0 | 1.0 | 0.1875 | 1.0 |
You see that the time index has the value ‘now’. This is the default index when no time series data has been specified and the network only covers a single state (e.g. a particular hour).
Similarly you will find the power flow in transmission lines at
n.lines_t.p0
| Line |
|---|
| snapshot |
| now |
n.lines_t.p1
| Line |
|---|
| snapshot |
| now |
The p0 will tell you the flow from bus0 to bus1. p1 will tell you the flow from bus1 to bus0.
What about the market clearing prices in the electricity market?
n.buses_t.marginal_price
| Bus | A | B |
|---|---|---|
| snapshot | ||
| now | 172.413793 | -0.0 |
Modifying networks#
Modifying data of components in an existing PyPSA network can be done by modifying the entries of a pandas.DataFrame. For instance, if we want to reduce the cross-border transmission capacity between Region A and Region B, we’d run:
n.links.loc["A-B", "p_nom"] = 400
n.links
| bus0 | bus1 | type | carrier | efficiency | active | build_year | lifetime | p_nom | p_nom_mod | ... | shut_down_cost | min_up_time | min_down_time | up_time_before | down_time_before | ramp_limit_up | ramp_limit_down | ramp_limit_start_up | ramp_limit_shut_down | p_nom_opt | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Link | |||||||||||||||||||||
| A-B | A | B | AC | 1.0 | True | 0 | inf | 400.0 | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 500.0 |
1 rows × 34 columns
n.optimize(solver_name="highs")
You can see that the production of the hydro power plant was reduced and that of the gas power plant increased owing to the reduced transmission capacity.
n.generators_t.p
| Generator | B hydro | A coal | A wind | A gas | A oil |
|---|---|---|---|---|---|
| snapshot | |||||
| now | 1050.0 | 35000.0 | 3000.0 | 1600.0 | 2000.0 |
Exploring the optimization model#
n.model
Linopy LP model
===============
Variables:
----------
* Generator-p (snapshot, Generator)
* Link-p (snapshot, Link)
Constraints:
------------
* Generator-fix-p-lower (snapshot, Generator-fix)
* Generator-fix-p-upper (snapshot, Generator-fix)
* Link-fix-p-lower (snapshot, Link-fix)
* Link-fix-p-upper (snapshot, Link-fix)
* Bus-nodal_balance (Bus, snapshot)
Status:
-------
ok
n.model.constraints["Generator-fix-p-upper"]
Constraint `Generator-fix-p-upper` [snapshot: 1, Generator-fix: 5]:
-------------------------------------------------------------------
[now, B hydro]: +1 Generator-p[now, B hydro] ≤ 1200.0
[now, A coal]: +1 Generator-p[now, A coal] ≤ 35000.0
[now, A wind]: +1 Generator-p[now, A wind] ≤ 3000.0
[now, A gas]: +1 Generator-p[now, A gas] ≤ 8000.0
[now, A oil]: +1 Generator-p[now, A oil] ≤ 2000.0
n.model.constraints["Bus-nodal_balance"]
Constraint `Bus-nodal_balance` [Bus: 2, snapshot: 1]:
-----------------------------------------------------
[A, now]: +1 Generator-p[now, A coal] + 1 Generator-p[now, A wind] + 1 Generator-p[now, A gas] + 1 Generator-p[now, A oil] - 1 Link-p[now, A-B] = 42000.0
[B, now]: +1 Generator-p[now, B hydro] + 1 Link-p[now, A-B] = 650.0
n.model.objective
Objective:
----------
LinearExpression: +24.24 Generator-p[now, A coal] + 172.4 Generator-p[now, A gas] + 137.1 Generator-p[now, A oil]
Sense: min
Value: 1398632.6317348
Exercises#
Task 1: Model an outage of the transmission line by setting it’s capacity to zero. How does the model compensate for the lack of transmission?
Task 2: Double the wind capacity. How is the electricity price affected? What generator is price-setting?
Task 3: Region A brings a new 1000 MW nuclear power plant into operation with an estimated marginal electricity generation cost of 20 €/MWh. Will that affect the electricity price?