Introduction to pypsa and linopy: a simple dispatch problem

Introduction to `pypsa` and `linopy`: a simple dispatch problem#

PyPSA stands for Python for Power System Analysis.

PyPSA is an open source Python package for simulating and optimising modern energy systems that include features such as

conventional generators with unit commitment (ramp-up, ramp-down, start-up, shut-down),
time-varying wind and solar generation,
energy storage with efficiency losses and inflow/spillage for hydroelectricity
coupling to other energy sectors (electricity, transport, heat, industry),
conversion between energy carriers (e.g. electricity to hydrogen),
transmission networks (AC, DC, other fuels)

PyPSA can be used for a variety of problem types (e.g. electricity market modelling, long-term investment planning, transmission network expansion planning), and is designed to scale well with large networks and long time series.

Compared to building power system by hand in linopy, PyPSA does the following things for you:

manage data inputs
build optimisation problem
communicate with the solver
retrieve and process optimisation results
manage data outputs

Dependencies#

pandas for storing data about network components and time series
numpy and scipy for linear algebra and sparse matrix calculations
matplotlib and cartopy for plotting on a map
networkx for network calculations
linopy for handling optimisation problems

Note

Documentation for this package is available at https://pypsa.readthedocs.io.

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 pypsa highspy

Basic components#

Component	Description
Network	Container for all components.
Bus	Node where components attach.
Carrier	Energy carrier or technology (e.g. electricity, hydrogen, gas, coal, oil, biomass, on-/offshore wind, solar). Can track properties such as specific carbon dioxide emissions or nice names and colors for plots.
Load	Energy consumer (e.g. electricity demand).
Generator	Generator (e.g. power plant, wind turbine, PV panel).
Line	Power distribution and transmission lines (overhead and cables).
Link	Links connect two buses with controllable energy flow, direction-control and losses. They can be used to model: HVDC links HVAC lines (neglecting KVL, only net transfer capacities (NTCs)) conversion between carriers (e.g. electricity to hydrogen in electrolysis)
StorageUnit	Storage with fixed nominal energy-to-power ratio.
Store	Storage with separately extendable energy capacity.
GlobalConstraint	Constraints affecting many components at once, such as emission limits.

Note

Links in the table lead to documentation for each component.

Getting started#

Simple electricity market problem#

generator 1: “gas” – marginal cost 70 EUR/MWh – capacity 50 MW

generator 2: “nuclear” – marginal cost 10 EUR/MWh – capacity 100 MW

load: “Consumer” – demand 120 MW

single time step (“now”)

single node (“Springfield”)

Building a basic network#

# By convention, PyPSA is imported without an alias:
import pypsa

# First, we create a network object which serves as the container for all components
n = pypsa.Network()

Empty PyPSA Network
Components: none
Snapshots: 1

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 (i.e. Kirchhoff’s Current Law).

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", "Springfield", v_nom=380, carrier="AC")

Index(['Springfield'], dtype='object')

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	v_mag_pu_set	v_mag_pu_min	v_mag_pu_max	control	generator	sub_network
Bus
Springfield	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.readthedocs.io/en/latest/components.html#bus, and analogous for all other components.

The n.add() function lets you add any component to the network object n:

n.add(
    "Generator",
    "gas",
    bus="Springfield",
    p_nom_extendable=False,
    marginal_cost=70,  # €/MWh
    p_nom=50,  # MW
)
n.add(
    "Generator",
    "nuclear",
    bus="Springfield",
    p_nom_extendable=False,
    marginal_cost=10,  # €/MWh
    p_nom=100,  # MW
)

Index(['nuclear'], dtype='object')

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.

As a result, the n.generators DataFrame looks like this:

n.generators

	bus	control	p_nom	p_nom_mod	p_nom_extendable	p_nom_min	p_nom_max	p_min_pu	p_max_pu	...	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	weight	p_nom_opt
Generator
gas	Springfield	PQ	50.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0
nuclear	Springfield	PQ	100.0	0.0	False	0.0	inf	0.0	1.0	...	0	0	1	0	NaN	NaN	1.0	1.0	1.0	0.0

2 rows × 37 columns

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",
    "Small town",
    bus="Springfield",
    p_set=120,  # MW
)

Index(['Small town'], dtype='object')

n.loads

	bus	carrier	type	p_set	q_set	sign	active
Load
Small town	Springfield			120.0	0.0	-1.0	True

Optimisation#

The design principle of PyPSA is that basically each component is associated with a set of variables and constraints that will be added to the optimisation model based on the input data stored for the components.

For this dispatch problem, PyPSA will solve an optimisation problem that looks like this

(1)#\[\begin{equation} \min_{g_{s,t};} \sum_{t,s} o_{s} g_{s,t} \end{equation}\]

such that

(2)#\[\begin{align} 0 & \leq g_{s,t} \leq G_{s} & \text{generation limits : generator} \\ D_t &= \sum_s g_{s,t} & \text{market clearing : bus} \\ \end{align}\]

Decision variables:

\(g_{s,t}\) is the generator dispatch of technology \(s\) at time \(t\)

Parameters:

\(o_{s}\) is the marginal generation cost of technology \(s\)
\(G_{s}\) is the nominal capacity of technology \(s\)
\(D_t\) is the power demand in Springfiled at time \(t\)

With all input data transferred into the PyPSA’s data structure (network), 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 adhering to the constraints defined in the network.

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(solver_name="highs")

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['Springfield'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['Springfield'], dtype='object', name='Bus')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io: Writing time: 0.01s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 2 primals, 5 duals
Objective: 2.40e+03
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper were not assigned to the network.

Running HiGHS 1.9.0 (git hash: fa40bdf): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e+01, 7e+01]
  Bound  [0e+00, 0e+00]
  RHS    [5e+01, 1e+02]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-5); columns 0(-2); elements 0(-6) - Reduced to empty
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-dl4iu_fg
Model status        : Optimal
Objective value     :  2.4000000000e+03
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
Writing the solution to /tmp/linopy-solve-wfdkvo0h.sol

('ok', 'optimal')

Let’s have a look at the results. The network object n contains now the objective value and the results for the decision variables.

n.objective

2400.0

Since the power flow and dispatch are generally time-varying quantities, these are stored in a different locations 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	gas	nuclear
snapshot
now	20.0	100.0

n.buses_t.marginal_price

Bus	Springfield
snapshot
now	70.0

n.generators_t.mu_upper

Generator
snapshot
now

Explore pypsa model#

n.model

Linopy LP model
===============

Variables:
----------
 * Generator-p (snapshot, Generator)

Constraints:
------------
 * Generator-fix-p-lower (snapshot, Generator-fix)
 * Generator-fix-p-upper (snapshot, Generator-fix)
 * Bus-nodal_balance (Bus, snapshot)

Status:
-------
ok

n.model.constraints

linopy.model.Constraints
------------------------
 * Generator-fix-p-lower (snapshot, Generator-fix)
 * Generator-fix-p-upper (snapshot, Generator-fix)
 * Bus-nodal_balance (Bus, snapshot)

n.model.constraints["Generator-fix-p-upper"]

Constraint `Generator-fix-p-upper` [snapshot: 1, Generator-fix: 2]:
-------------------------------------------------------------------
[now, gas]: +1 Generator-p[now, gas]         ≤ 50.0
[now, nuclear]: +1 Generator-p[now, nuclear] ≤ 100.0

n.model.constraints["Bus-nodal_balance"]

Constraint `Bus-nodal_balance` [Bus: 1, snapshot: 1]:
-----------------------------------------------------
[Springfield, now]: +1 Generator-p[now, gas] + 1 Generator-p[now, nuclear] = 120.0

n.model.objective

Objective:
----------
LinearExpression: +70 Generator-p[now, gas] + 10 Generator-p[now, nuclear]
Sense: min
Value: 2400.0

(If time allows) Let’s write optimization problem manually amd reproduce the PyPSA model#

PyPSA optimisation module is based on Linopy — an open-source framework for formulating, solving, and analyzing optimization problems with Python.

With Linopy, you can create optimization models within Python that consist of decision variables, constraints, and optimization objectives. You can then solve these instances using a variety of commercial and open-source solvers (specialised software).

Linopy supports a wide range of problem types, including:

Linear programming
Integer programming
Mixed-integer programming
Quadratic programming

Note

Documentation for this package is available at https://linopy.readthedocs.io.

import linopy

# remove all constraints
for name in list(n.model.constraints):
    n.model.remove_constraints(name)

n.model.constraints

linopy.model.Constraints
------------------------
<empty>

# remove objective
n.model.objective.expression = linopy.LinearExpression(None, model=n.model)

n.model.objective.expression

LinearExpression
----------------
+0

# set objective expression to the problem
n.model.objective.expression = (
    n.model["Generator-p"].loc[:, "gas"] * n.generators.marginal_cost.loc["gas"]
    + n.model["Generator-p"].loc[:, "nuclear"]
    * n.generators.marginal_cost.loc["nuclear"]
)

# check it is there
n.model.objective

Objective:
----------
LinearExpression: +70 Generator-p[now, gas] + 10 Generator-p[now, nuclear]
Sense: min
Value: 2400.0

n.model.add_constraints(
    n.model["Generator-p"].sum(dim="Generator") == n.loads.p_set["Small town"],
    name="nodal_balance",
)

Constraint `nodal_balance` [snapshot: 1]:
-----------------------------------------
[now]: +1 Generator-p[now, gas] + 1 Generator-p[now, nuclear] = 120.0

n.model.add_constraints(n.model["Generator-p"].loc[:, "gas"] >= 0, name="p_lower_gas")

Constraint `p_lower_gas` [snapshot: 1]:
---------------------------------------
[now]: +1 Generator-p[now, gas] ≥ -0.0

n.model.add_constraints(
    n.model["Generator-p"].loc[:, "nuclear"] >= 0, name="p_lower_coal"
)

Constraint `p_lower_coal` [snapshot: 1]:
----------------------------------------
[now]: +1 Generator-p[now, nuclear] ≥ -0.0

n.model.add_constraints(
    n.model["Generator-p"].loc[:, "gas"] <= n.generators.p_nom.loc["gas"],
    name="p_upper_gas",
)

Constraint `p_upper_gas` [snapshot: 1]:
---------------------------------------
[now]: +1 Generator-p[now, gas] ≤ 50.0

n.model.add_constraints(
    n.model["Generator-p"].loc[:, "nuclear"] <= n.generators.p_nom.loc["nuclear"],
    name="p_upper_nuclear",
)

Constraint `p_upper_nuclear` [snapshot: 1]:
-------------------------------------------
[now]: +1 Generator-p[now, nuclear] ≤ 100.0

# check that we did a good job
n.model.constraints

linopy.model.Constraints
------------------------
 * nodal_balance (snapshot)
 * p_lower_gas (Generator, snapshot)
 * p_lower_coal (Generator, snapshot)
 * p_upper_gas (Generator, snapshot)
 * p_upper_nuclear (Generator, snapshot)

Let’s ensure that we get the same results as when using the n.optimize() function

n.model.solve(solver_name="highs")

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io: Writing time: 0.02s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 2 primals, 5 duals
Objective: 2.40e+03
Solver model: available
Solver message: optimal

Running HiGHS 1.9.0 (git hash: fa40bdf): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e+01, 7e+01]
  Bound  [0e+00, 0e+00]
  RHS    [5e+01, 1e+02]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-5); columns 0(-2); elements 0(-6) - Reduced to empty
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-6jauvs59
Model status        : Optimal
Objective value     :  2.4000000000e+03
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
Writing the solution to /tmp/linopy-solve-3cwnm_q6.sol

('ok', 'optimal')

n.generators_t.p

Generator	gas	nuclear
snapshot
now	20.0	100.0

Introduction to pypsa and linopy: a simple dispatch problem

Contents

Introduction to pypsa and linopy: a simple dispatch problem#

Dependencies#

Basic components#

Getting started#

Simple electricity market problem#

Building a basic network#

Optimisation#

Explore pypsa model#

(If time allows) Let’s write optimization problem manually amd reproduce the PyPSA model#

Introduction to `pypsa` and `linopy`: a simple dispatch problem#