Note

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

Screening curve analysis

Compute the long-term equilibrium power plant investment for a given load duration curve (1000-1000z for z \(\in\) [0,1]) and a given set of generator investment options.

[1]:
import pypsa
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

%matplotlib inline
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [1], line 1
----> 1 import pypsa
      2 import numpy as np
      3 import pandas as pd

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

Generator marginal (m) and capital (c) costs in EUR/MWh - numbers chosen for simple answer.

[2]:
generators = {
    "coal": {"m": 2, "c": 15},
    "gas": {"m": 12, "c": 10},
    "load-shedding": {"m": 1012, "c": 0},
}

The screening curve intersections are at 0.01 and 0.5.

[3]:
x = np.linspace(0, 1, 101)
df = pd.DataFrame(
    {key: pd.Series(item["c"] + x * item["m"], x) for key, item in generators.items()}
)
df.plot(ylim=[0, 50], title="Screening Curve", figsize=(9, 5))
plt.tight_layout()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [3], line 1
----> 1 x = np.linspace(0, 1, 101)
      2 df = pd.DataFrame(
      3     {key: pd.Series(item["c"] + x * item["m"], x) for key, item in generators.items()}
      4 )
      5 df.plot(ylim=[0, 50], title="Screening Curve", figsize=(9, 5))

NameError: name 'np' is not defined
[4]:
n = pypsa.Network()

num_snapshots = 1001
n.snapshots = np.linspace(0, 1, num_snapshots)
n.snapshot_weightings = n.snapshot_weightings / num_snapshots

n.add("Bus", name="bus")

n.add("Load", name="load", bus="bus", p_set=1000 - 1000 * n.snapshots.values)

for gen in generators:
    n.add(
        "Generator",
        name=gen,
        bus="bus",
        p_nom_extendable=True,
        marginal_cost=float(generators[gen]["m"]),
        capital_cost=float(generators[gen]["c"]),
    )
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [4], line 1
----> 1 n = pypsa.Network()
      3 num_snapshots = 1001
      4 n.snapshots = np.linspace(0, 1, num_snapshots)

NameError: name 'pypsa' is not defined
[5]:
n.loads_t.p_set.plot.area(title="Load Duration Curve", figsize=(9, 5), ylabel="MW")
plt.tight_layout()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [5], line 1
----> 1 n.loads_t.p_set.plot.area(title="Load Duration Curve", figsize=(9, 5), ylabel="MW")
      2 plt.tight_layout()

NameError: name 'n' is not defined
[6]:
n.lopf(solver_name="cbc")
n.objective
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [6], line 1
----> 1 n.lopf(solver_name="cbc")
      2 n.objective

NameError: name 'n' is not defined

The capacity is set by total electricity required.

NB: No load shedding since all prices are below 10 000.

[7]:
n.generators.p_nom_opt.round(2)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [7], line 1
----> 1 n.generators.p_nom_opt.round(2)

NameError: name 'n' is not defined
[8]:
n.buses_t.marginal_price.plot(title="Price Duration Curve", figsize=(9, 4))
plt.tight_layout()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [8], line 1
----> 1 n.buses_t.marginal_price.plot(title="Price Duration Curve", figsize=(9, 4))
      2 plt.tight_layout()

NameError: name 'n' is not defined

The prices correspond either to VOLL (1012) for first 0.01 or the marginal costs (12 for 0.49 and 2 for 0.5)

Except for (infinitesimally small) points at the screening curve intersections, which correspond to changing the load duration near the intersection, so that capacity changes. This explains 7 = (12+10 - 15) (replacing coal with gas) and 22 = (12+10) (replacing load-shedding with gas).

Note: What remains unclear is what is causing :nbsphinx-math:`l `= 0… it should be 2.

[9]:
n.buses_t.marginal_price.round(2).sum(axis=1).value_counts()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [9], line 1
----> 1 n.buses_t.marginal_price.round(2).sum(axis=1).value_counts()

NameError: name 'n' is not defined
[10]:
n.generators_t.p.plot(ylim=[0, 600], title="Generation Dispatch", figsize=(9, 5))
plt.tight_layout()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [10], line 1
----> 1 n.generators_t.p.plot(ylim=[0, 600], title="Generation Dispatch", figsize=(9, 5))
      2 plt.tight_layout()

NameError: name 'n' is not defined

Demonstrate zero-profit condition.

  1. The total cost is given by

[11]:
(
    n.generators.p_nom_opt * n.generators.capital_cost
    + n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0).sum()
    * n.generators.marginal_cost
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [11], line 2
      1 (
----> 2     n.generators.p_nom_opt * n.generators.capital_cost
      3     + n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0).sum()
      4     * n.generators.marginal_cost
      5 )

NameError: name 'n' is not defined
  1. The total revenue by

[12]:
(
    n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0)
    .multiply(n.buses_t.marginal_price["bus"], axis=0)
    .sum(0)
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [12], line 2
      1 (
----> 2     n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0)
      3     .multiply(n.buses_t.marginal_price["bus"], axis=0)
      4     .sum(0)
      5 )

NameError: name 'n' is not defined

Now, take the capacities from the above long-term equilibrium, then disallow expansion.

Show that the resulting market prices are identical.

This holds in this example, but does NOT necessarily hold and breaks down in some circumstances (for example, when there is a lot of storage and inter-temporal shifting).

[13]:
n.generators.p_nom_extendable = False
n.generators.p_nom = n.generators.p_nom_opt
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [13], line 1
----> 1 n.generators.p_nom_extendable = False
      2 n.generators.p_nom = n.generators.p_nom_opt

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

NameError: name 'n' is not defined
[15]:
n.buses_t.marginal_price.plot(title="Price Duration Curve", figsize=(9, 5))
plt.tight_layout()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [15], line 1
----> 1 n.buses_t.marginal_price.plot(title="Price Duration Curve", figsize=(9, 5))
      2 plt.tight_layout()

NameError: name 'n' is not defined
[16]:
n.buses_t.marginal_price.sum(axis=1).value_counts()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [16], line 1
----> 1 n.buses_t.marginal_price.sum(axis=1).value_counts()

NameError: name 'n' is not defined

Demonstrate zero-profit condition. Differences are due to singular times, see above, not a problem

  1. Total costs

[17]:
(
    n.generators.p_nom * n.generators.capital_cost
    + n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0).sum()
    * n.generators.marginal_cost
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [17], line 2
      1 (
----> 2     n.generators.p_nom * n.generators.capital_cost
      3     + n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0).sum()
      4     * n.generators.marginal_cost
      5 )

NameError: name 'n' is not defined
  1. Total revenue

[18]:
(
    n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0)
    .multiply(n.buses_t.marginal_price["bus"], axis=0)
    .sum()
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [18], line 2
      1 (
----> 2     n.generators_t.p.multiply(n.snapshot_weightings.generators, axis=0)
      3     .multiply(n.buses_t.marginal_price["bus"], axis=0)
      4     .sum()
      5 )

NameError: name 'n' is not defined