Dispatch problems, Maps & Statistics

Recap: basic components¶

A more general form of an electricity dispatch problem¶

For an hourly electricity market dispatch simulation, PyPSA will solve an optimisation problem that looks like this

such that

Decision variables:

• $g_{i,s,t}$ is the generator dispatch at bus $i$, technology $s$, time step $t$,

• $f_{\ell,t}$ is the power flow in line $\ell$,

• $g_{i,r,t,\text{dis-/charge}}$ denotes the charge and discharge of storage unit $r$ at bus $i$ and time step $t$,

• $e_{i,r,t}$ is the state of charge of storage $r$ at bus $i$ and time step $t$.

Parameters:

• $o_{i,s}$ is the marginal generation cost of technology $s$ at bus $i$,

• $x_\ell$ is the reactance of transmission line $\ell$,

• $K_{i\ell}$ is the incidence matrix,

• $C_{\ell c}$ is the cycle matrix,

• $G_{i,s}$ is the nominal capacity of the generator of technology $s$ at bus $i$,

• $F_{\ell}$ is the rating of the transmission line $\ell$,

• $E_{i,r}$ is the energy capacity of storage $r$ at bus $i$,

• $\eta^{0/1/2}_{i,r,t}$ denote the standing (0), charging (1), and discharging (2) efficiencies.

South Africa & Mozambique system example¶

Compared to the previous example, we will consider a more complex system with more components (generators, transmission lines) and more buses. We also discuss basics of the plotting functionality built into PyPSA.

We have the following data:

• fuel costs in € / MWh$_{th}$

• efficiencies of thermal power plants in MWh$_{el}$ / MWh$_{th}$

• specific emissions in t$_{CO_2}$ / MWh$_{th}$

• power plant capacities in MW

• electrical load in MW

Building the network¶

By convention, PyPSA is imported without an alias:

First, we create a new network object which serves as the overall container for all components.

This time, we add two buses, one for each country, and provide geographic coordinates (x is longitude, y is latitude) so that we can later plot the network on a map:

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.

Let’s add a generator in Mozambique:

Add generators in South Africa (in a loop):

The complete n.generators DataFrame looks like this now:

Next, we’re going to add the electricity demand.

A positive value for p_set means consumption of power from the bus (in MW).

Finally, we add the connection between Mozambique and South Africa with a 500 MW line:

Optimisation¶

With all input data transferred into PyPSA’s data structure, we can now build and run the resulting optimisation problem. We can have a look at the optimisation problem that will be solved by PyPSA with the n.optimize.create_model() function. This function returns a linopy model object:

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 in your environment, to use them here!

Since we have already inspected the solver output in the previous notebook, we now suppress it with log_to_console=False to keep the notebook tidy:

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

or if you prefer it in relation to the generators nominal capacity

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

The p0 will tell you the flow from bus0 to bus1. p1 will tell you the flow from bus1 to bus0.

What about the shadow prices?

Model inspection¶

As we saw in the previous notebook, the underlying linopy optimisation model is accessible via n.model, e.g. to inspect variables, constraints — or, in this case, the shadow prices (duals) of the generation limits:

Basic network plotting¶

For plotting PyPSA network, we’re going to need the help of some friends:

PyPSA has a built-in plotting function based on matplotlib:

Since we have provided x and y coordinates for our buses, n.plot() will try to plot the network on a map by default. Of course, there’s an option to deactivate this behaviour:

The n.plot() function has a variety of styling arguments to tweak the appearance of the buses, the lines and the map in the background:

We can use the bus_sizes argument of n.plot() to display the regional distribution of load. First, we calculate the total load per bus:

The resulting pandas.Series we can pass to n.plot(bus_sizes=...):

The important point here is, that s needs to have entries for all buses, i.e. its index needs to match n.buses.index.

The bus_sizes argument of n.plot() can be even more powerful. It can produce pie charts, e.g. for the mix of electricity generation at each bus.

The dispatch of each generator, we can find at:

If we group this by the bus and carrier...

... we get a multi-indexed pandas.Series ...

... which we can pass to n.plot(bus_sizes=...):

How does this magic work? The plotting function will look up the colors specified in n.carriers for each carrier and match it with the second index-level of s.

Besides the static plots with n.plot(), PyPSA also has an interactive plotting function n.explore(), which renders the network on a zoomable map. We will see more of it in the SciGRID example below, but it already works for our small two-country network:

Modifying networks¶

Modifying data of components in an existing PyPSA network is as easy as modifying the entries of a pandas.DataFrame. For instance, if we want to reduce the cross-border transmission capacity between South Africa and Mozambique, we’d run:

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.

Global constraints for emission limits¶

In the example above, we happen to have some spare gas capacity with lower carbon intensity than the coal and oil generators. We could use this to lower the emissions of the system, but it will be more expensive. We can implement the limit of carbon dioxide emissions as a constraint.

This is achieved in PyPSA through Global Constraints which add constraints that apply to many components at once.

But first, we need to calculate the current level of emissions to set a sensible limit.

We can compute the emissions per generator (in tonnes of CO$_2$) in the following way.

where $\rho$ is the specific emissions (tonnes/MWh thermal) and $\eta$ is the conversion efficiency (MWh electric / MWh thermal) of the generator with dispatch $g$ (MWh electric):

Summed up, we get total emissions in tonnes:

So, let’s say we want to reduce emissions by 10%:

Let’s check how the new global constraint looks like in the optimisation problem:

The shadow price of the emission limit tells us by how much the system cost would decrease if we allowed one more tonne of CO$_2$:

Can we lower emissions even further? Say by another 5% points?

No! Without any additional capacities, we have exhausted our options to reduce emissions in that hour. The solver tells us that the problem is infeasible, i.e. there is no solution that satisfies all constraints. Better revert this change:

Data import and export¶

You may find yourself in a need to store PyPSA networks for later use. Or, maybe you want to import the genius PyPSA example that someone else uploaded to the web to explore.

Among other file formats, PyPSA networks can be stored as netCDF (.nc) files, as folders of CSV files, or as Excel (.xlsx) files.

• netCDF files have the advantage that they take up less space than CSV files and are faster to load.

• CSV and Excel might be easier to inspect by hand.

The approach is similar for all formats:

Example grid model with time resolution¶

Dispatch problem with German SciGRID network

SciGRID is a project that provides an open reference model of the European transmission network. The network comprises time series for loads and the availability of renewable generation at an hourly resolution for January 1, 2011 as well as approximate generation capacities in 2014. This dataset is a little out of date and only intended to demonstrate the capabilities of PyPSA.

We load it as a second network n2, so that our two-country network n stays available for the exercises at the end:

There are some infeasibilities without allowing extension. Moreover, to approximate so-called $N-1$ security, we don’t allow any line to be loaded above 70% of their thermal rating. $N-1$ security is a constraint that states that no single transmission line may be overloaded by the failure of another transmission line (e.g. through a tripped connection).

Because this network includes time-varying data, now is the time to look at another attribute of n: n.snapshots. Snapshots is the PyPSA terminology for time steps. In most cases, they represent a particular hour. They can be a pandas.DatetimeIndex or any other list-like attributes.

This index will match with any time-varying attributes of components:

We can use simple pandas syntax, to create an overview of the load time series...

... and the capacity factor time series of the renewable generators:

We can also inspect the total power plant capacities per technology:

Interactive network plotting¶

The network inputs and outputs can also be visualised on a map. The n.plot() function we already know uses matplotlib to create static plots, while n.explore() creates interactive plots based on pydeck. In the following, we will focus on the interactive plotting with n.explore().

The n.explore() function has a variety of styling arguments to tweak the appearance of the buses, the lines and the map in the background. For example, we can size the buses according to their load:

The bus_size argument of n.explore() can be even more powerful. It can show the nodal composition of power plant capacities as pie charts by grouping data by the bus and carrier attributes of the generators:

... for which we need to assign some colors to the technologies first:

So let’s solve the electricity market simulation for January 1, 2011. It’ll take a short moment.

Now, we can also plot model outputs, like the calculated power flows on the network map using the line_flow argument of n.explore():

Static plots remain useful for reports and publications. For instance, we can plot the average locational marginal prices (LMPs) with n.plot():

Statistics¶

The PyPSA n.statistics module provides a variety of pre-defined tables and plots to analyse optimisation results at different aggregation levels. For instance, the energy balance, the operational costs (OPEX), curtailment, prices and market value. We will use this module in more detail later in this workshop.

For example, we can plot the hourly dispatch grouped by carrier with a single line of code:

Or look at the energy balance as a table (in GWh):

Or the operational costs (OPEX) per technology:

There is much more to explore in PyPSA. If you are hooked, have a look at the documentation and the examples section.

Exercises¶

Modify some of the input data of the South Africa & Mozambique network n, first removing the emission limit global constraint:

Task 1: Model an outage of the transmission line by removing it. How does the model compensate for the lack of transmission?

Task 2: Double the wind capacity. How is the electricity price in South Africa affected? What generator is price-setting?

Task 3: South Africa 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?