Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Redispatch Example with SciGRID network

In this example, we compare a 2-stage market with an initial market clearing in two bidding zones with flow-based market coupling and a subsequent redispatch market (incl. curtailment) to an idealised nodal pricing scheme.

[1]:
import pypsa
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from pypsa.descriptors import get_switchable_as_dense as as_dense
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [1], line 1
----> 1 import pypsa
      2 import matplotlib.pyplot as plt
      3 import cartopy.crs as ccrs

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/__init__.py:10
      1 # -*- coding: utf-8 -*-
      4 """
      5 Python for Power Systems Analysis (PyPSA)
      6
      7 Grid calculation library.
      8 """
---> 10 from pypsa import (
     11     components,
     12     contingency,
     13     descriptors,
     14     examples,
     15     geo,
     16     io,
     17     linopf,
     18     linopt,
     19     networkclustering,
     20     opf,
     21     opt,
     22     optimization,
     23     pf,
     24     plot,
     25 )
     26 from pypsa.components import Network, SubNetwork
     28 __version__ = "0.21.2"

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/components.py:50
     37 from pypsa.io import (
     38     export_to_csv_folder,
     39     export_to_hdf5,
   (...)
     47     import_series_from_dataframe,
     48 )
     49 from pypsa.opf import network_lopf, network_opf
---> 50 from pypsa.optimization.optimize import OptimizationAccessor
     51 from pypsa.pf import (
     52     calculate_B_H,
     53     calculate_dependent_values,
   (...)
     62     sub_network_pf,
     63 )
     64 from pypsa.plot import iplot, plot

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/optimization/__init__.py:7
      1 #!/usr/bin/env python3
      2 # -*- coding: utf-8 -*-
      3 """
      4 Build optimisation problems from PyPSA networks with Linopy.
      5 """
----> 7 from pypsa.optimization import abstract, constraints, optimize, variables
      8 from pypsa.optimization.optimize import create_model

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/optimization/constraints.py:9
      6 import logging
      8 import pandas as pd
----> 9 from linopy.expressions import LinearExpression, merge
     10 from numpy import arange, cumsum, inf, nan, roll
     11 from scipy import sparse

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/__init__.py:9
      1 #!/usr/bin/env python3
      2 # -*- coding: utf-8 -*-
      3 """
      4 Created on Wed Mar 10 11:03:06 2021.
      5
      6 @author: fabulous
      7 """
----> 9 from linopy import model, remote
     10 from linopy.expressions import merge
     11 from linopy.io import read_netcdf

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/model.py:22
     20 from linopy import solvers
     21 from linopy.common import best_int, replace_by_map
---> 22 from linopy.constraints import (
     23     AnonymousConstraint,
     24     AnonymousScalarConstraint,
     25     Constraints,
     26 )
     27 from linopy.eval import Expr
     28 from linopy.expressions import LinearExpression, ScalarLinearExpression

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/constraints.py:21
     18 from scipy.sparse import coo_matrix
     19 from xarray import DataArray, Dataset
---> 21 from linopy import expressions, variables
     22 from linopy.common import (
     23     _merge_inplace,
     24     has_assigned_model,
   (...)
     27     replace_by_map,
     28 )
     31 class Constraint(DataArray):

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/expressions.py:23
     20 from xarray.core.dataarray import DataArrayCoordinates
     21 from xarray.core.groupby import _maybe_reorder, peek_at
---> 23 from linopy import constraints, variables
     24 from linopy.common import as_dataarray
     27 def exprwrap(method, *default_args, **new_default_kwargs):

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/variables.py:398
    393     roll = varwrap(DataArray.roll)
    395     rolling = varwrap(DataArray.rolling)
--> 398 @dataclass(repr=False)
    399 class Variables:
    400     """
    401     A variables container used for storing multiple variable arrays.
    402     """
    404     labels: Dataset = Dataset()

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:1211, in dataclass.<locals>.wrap(cls)
   1210 def wrap(cls):
-> 1211     return _process_class(cls, init, repr, eq, order, unsafe_hash,
   1212                           frozen, match_args, kw_only, slots,
   1213                           weakref_slot)

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:959, in _process_class(cls, init, repr, eq, order, unsafe_hash, frozen, match_args, kw_only, slots, weakref_slot)
    956         kw_only = True
    957     else:
    958         # Otherwise it's a field of some type.
--> 959         cls_fields.append(_get_field(cls, name, type, kw_only))
    961 for f in cls_fields:
    962     fields[f.name] = f

File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:816, in _get_field(cls, a_name, a_type, default_kw_only)
    812 # For real fields, disallow mutable defaults.  Use unhashable as a proxy
    813 # indicator for mutability.  Read the __hash__ attribute from the class,
    814 # not the instance.
    815 if f._field_type is _FIELD and f.default.__class__.__hash__ is None:
--> 816     raise ValueError(f'mutable default {type(f.default)} for field '
    817                      f'{f.name} is not allowed: use default_factory')
    819 return f

ValueError: mutable default <class 'xarray.core.dataset.Dataset'> for field labels is not allowed: use default_factory
[2]:
solver = "cbc"

Load example network

[3]:
o = pypsa.examples.scigrid_de(from_master=True)
o.lines.s_max_pu = 0.7
o.lines.loc[["316", "527", "602"], "s_nom"] = 1715
o.set_snapshots([o.snapshots[12]])
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [3], line 1
----> 1 o = pypsa.examples.scigrid_de(from_master=True)
      2 o.lines.s_max_pu = 0.7
      3 o.lines.loc[["316", "527", "602"], "s_nom"] = 1715

NameError: name 'pypsa' is not defined
[4]:
n = o.copy()  # for redispatch model
m = o.copy()  # for market model
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [4], line 1
----> 1 n = o.copy()  # for redispatch model
      2 m = o.copy()  # for market model

NameError: name 'o' is not defined
[5]:
o.plot();
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [5], line 1
----> 1 o.plot();

NameError: name 'o' is not defined

Solve original nodal market model o

First, let us solve a nodal market using the original model o:

[6]:
o.lopf(solver_name=solver, pyomo=False)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [6], line 1
----> 1 o.lopf(solver_name=solver, pyomo=False)

NameError: name 'o' is not defined

Costs are 301 k€.

Build market model m with two bidding zones

For this example, we split the German transmission network into two market zones at latitude 51 degrees.

You can build any other market zones by providing an alternative mapping from bus to zone.

[7]:
zones = (n.buses.y > 51).map(lambda x: "North" if x else "South")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [7], line 1
----> 1 zones = (n.buses.y > 51).map(lambda x: "North" if x else "South")

NameError: name 'n' is not defined

Next, we assign this mapping to the market model m.

We re-assign the buses of all generators and loads, and remove all transmission lines within each bidding zone.

Here, we assume that the bidding zones are coupled through the transmission lines that connect them.

[8]:
for c in m.iterate_components(m.one_port_components):
    c.df.bus = c.df.bus.map(zones)

for c in m.iterate_components(m.branch_components):
    c.df.bus0 = c.df.bus0.map(zones)
    c.df.bus1 = c.df.bus1.map(zones)
    internal = c.df.bus0 == c.df.bus1
    m.mremove(c.name, c.df.loc[internal].index)

m.mremove("Bus", m.buses.index)
m.madd("Bus", ["North", "South"]);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [8], line 1
----> 1 for c in m.iterate_components(m.one_port_components):
      2     c.df.bus = c.df.bus.map(zones)
      4 for c in m.iterate_components(m.branch_components):

NameError: name 'm' is not defined

Now, we can solve the coupled market with two bidding zones.

[9]:
m.lopf(solver_name=solver, pyomo=False)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [9], line 1
----> 1 m.lopf(solver_name=solver, pyomo=False)

NameError: name 'm' is not defined

Costs are 214 k€, which is much lower than the 301 k€ of the nodal market.

This is because network restrictions apart from the North/South division are not taken into account yet.

We can look at the market clearing prices of each zone:

[10]:
m.buses_t.marginal_price
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [10], line 1
----> 1 m.buses_t.marginal_price

NameError: name 'm' is not defined

Build redispatch model n

Next, based on the market outcome with two bidding zones m, we build a secondary redispatch market n that rectifies transmission constraints through curtailment and ramping up/down thermal generators.

First, we fix the dispatch of generators to the results from the market simulation. (For simplicity, this example disregards storage units.)

[11]:
p = m.generators_t.p / m.generators.p_nom
n.generators_t.p_min_pu = p
n.generators_t.p_max_pu = p
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [11], line 1
----> 1 p = m.generators_t.p / m.generators.p_nom
      2 n.generators_t.p_min_pu = p
      3 n.generators_t.p_max_pu = p

NameError: name 'm' is not defined

Then, we add generators bidding into redispatch market using the following assumptions:

  • All generators can reduce their dispatch to zero. This includes also curtailment of renewables.

  • All generators can increase their dispatch to their available/nominal capacity.

  • No changes to the marginal costs, i.e. reducing dispatch lowers costs.

With these settings, the 2-stage market should result in the same cost as the nodal market.

[12]:
g_up = n.generators.copy()
g_down = n.generators.copy()

g_up.index = g_up.index.map(lambda x: x + " ramp up")
g_down.index = g_down.index.map(lambda x: x + " ramp down")

up = (
    as_dense(m, "Generator", "p_max_pu") * m.generators.p_nom - m.generators_t.p
).clip(0) / m.generators.p_nom
down = -m.generators_t.p / m.generators.p_nom

up.columns = up.columns.map(lambda x: x + " ramp up")
down.columns = down.columns.map(lambda x: x + " ramp down")

n.madd("Generator", g_up.index, p_max_pu=up, **g_up.drop("p_max_pu", axis=1))

n.madd(
    "Generator",
    g_down.index,
    p_min_pu=down,
    p_max_pu=0,
    **g_down.drop(["p_max_pu", "p_min_pu"], axis=1)
);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [12], line 1
----> 1 g_up = n.generators.copy()
      2 g_down = n.generators.copy()
      4 g_up.index = g_up.index.map(lambda x: x + " ramp up")

NameError: name 'n' is not defined

Now, let’s solve the redispatch market:

[13]:
n.lopf(solver_name=solver, pyomo=False)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [13], line 1
----> 1 n.lopf(solver_name=solver, pyomo=False)

NameError: name 'n' is not defined

And, as expected, the costs are the same as for the nodal market: 301 k€.

Now, we can plot both the market results of the 2 bidding zone market and the redispatch results:

[14]:
fig, axs = plt.subplots(
    1, 3, figsize=(20, 10), subplot_kw={"projection": ccrs.AlbersEqualArea()}
)

market = (
    n.generators_t.p[m.generators.index]
    .T.squeeze()
    .groupby(n.generators.bus)
    .sum()
    .div(2e4)
)
n.plot(ax=axs[0], bus_sizes=market, title="2 bidding zones market simulation")

redispatch_up = (
    n.generators_t.p.filter(like="ramp up")
    .T.squeeze()
    .groupby(n.generators.bus)
    .sum()
    .div(2e4)
)
n.plot(
    ax=axs[1], bus_sizes=redispatch_up, bus_colors="blue", title="Redispatch: ramp up"
)

redispatch_down = (
    n.generators_t.p.filter(like="ramp down")
    .T.squeeze()
    .groupby(n.generators.bus)
    .sum()
    .div(-2e4)
)
n.plot(
    ax=axs[2],
    bus_sizes=redispatch_down,
    bus_colors="red",
    title="Redispatch: ramp down / curtail",
);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [14], line 1
----> 1 fig, axs = plt.subplots(
      2     1, 3, figsize=(20, 10), subplot_kw={"projection": ccrs.AlbersEqualArea()}
      3 )
      5 market = (
      6     n.generators_t.p[m.generators.index]
      7     .T.squeeze()
   (...)
     10     .div(2e4)
     11 )
     12 n.plot(ax=axs[0], bus_sizes=market, title="2 bidding zones market simulation")

NameError: name 'plt' is not defined

We can also read out the final dispatch of each generator:

[15]:
grouper = n.generators.index.str.split(" ramp", expand=True).get_level_values(0)

n.generators_t.p.groupby(grouper, axis=1).sum().squeeze()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [15], line 1
----> 1 grouper = n.generators.index.str.split(" ramp", expand=True).get_level_values(0)
      3 n.generators_t.p.groupby(grouper, axis=1).sum().squeeze()

NameError: name 'n' is not defined

Changing bidding strategies in redispatch market

We can also formulate other bidding strategies or compensation mechanisms for the redispatch market.

For example, that ramping up a generator is twice as expensive.

[16]:
n.generators.loc[n.generators.index.str.contains("ramp up"), "marginal_cost"] *= 2
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [16], line 1
----> 1 n.generators.loc[n.generators.index.str.contains("ramp up"), "marginal_cost"] *= 2

NameError: name 'n' is not defined

Or that generators need to be compensated for curtailing them or ramping them down at 50% of their marginal cost.

[17]:
n.generators.loc[n.generators.index.str.contains("ramp down"), "marginal_cost"] *= -0.5
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [17], line 1
----> 1 n.generators.loc[n.generators.index.str.contains("ramp down"), "marginal_cost"] *= -0.5

NameError: name 'n' is not defined

In this way, the outcome should be more expensive than the ideal nodal market:

[18]:
n.lopf(solver_name=solver, pyomo=False)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [18], line 1
----> 1 n.lopf(solver_name=solver, pyomo=False)

NameError: name 'n' is not defined

Costs are now 559 k€ compared to 301 k€.

[ ]:

[ ]: