Note

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

Biomass, synthetic fuels and carbon management

In this example we show how to manage different biomass stocks with different potentials and costs, synthetic fuel production, direct air capture (DAC) and carbon capture and usage/sequestration/cycling (CCU/S/C).

Demand for electricity and diesel transport have to be met from various biomass sources, natural gas with possibility for carbon capture, electrolysis for hydrogen production, direct air capture of CO2, and diesel synthesis via Fischer-Tropsch.

The system has to reach a target of net negative emissions over the period.

All numbers/costs/efficiencies are fictitious to allow easy analysis.

[1]:
import pypsa
import numpy as np
import matplotlib.pyplot as plt
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [1], line 1
----> 1 import pypsa
      2 import numpy as np
      3 import matplotlib.pyplot as plt

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

First tell PyPSA that links can have multiple outputs by overriding the component_attrs. This can be done for as many buses as you need with format busi for i = 2,3,4,5,….

See https://pypsa.org/doc/components.html#link-with-multiple-outputs-or-inputs

[2]:
override_component_attrs = pypsa.descriptors.Dict(
    {k: v.copy() for k, v in pypsa.components.component_attrs.items()}
)
override_component_attrs["Link"].loc["bus2"] = [
    "string",
    np.nan,
    np.nan,
    "2nd bus",
    "Input (optional)",
]
override_component_attrs["Link"].loc["bus3"] = [
    "string",
    np.nan,
    np.nan,
    "3rd bus",
    "Input (optional)",
]
override_component_attrs["Link"].loc["efficiency2"] = [
    "static or series",
    "per unit",
    1.0,
    "2nd bus efficiency",
    "Input (optional)",
]
override_component_attrs["Link"].loc["efficiency3"] = [
    "static or series",
    "per unit",
    1.0,
    "3rd bus efficiency",
    "Input (optional)",
]
override_component_attrs["Link"].loc["p2"] = [
    "series",
    "MW",
    0.0,
    "2nd bus output",
    "Output",
]
override_component_attrs["Link"].loc["p3"] = [
    "series",
    "MW",
    0.0,
    "3rd bus output",
    "Output",
]
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [2], line 1
----> 1 override_component_attrs = pypsa.descriptors.Dict(
      2     {k: v.copy() for k, v in pypsa.components.component_attrs.items()}
      3 )
      4 override_component_attrs["Link"].loc["bus2"] = [
      5     "string",
      6     np.nan,
   (...)
      9     "Input (optional)",
     10 ]
     11 override_component_attrs["Link"].loc["bus3"] = [
     12     "string",
     13     np.nan,
   (...)
     16     "Input (optional)",
     17 ]

NameError: name 'pypsa' is not defined
[3]:
n = pypsa.Network(override_component_attrs=override_component_attrs)
n.set_snapshots(range(10))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [3], line 1
----> 1 n = pypsa.Network(override_component_attrs=override_component_attrs)
      2 n.set_snapshots(range(10))

NameError: name 'pypsa' is not defined

Add a constant electrical load

[4]:
n.add("Bus", "bus")
n.add("Load", "load", bus="bus", p_set=1.0)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [4], line 1
----> 1 n.add("Bus", "bus")
      2 n.add("Load", "load", bus="bus", p_set=1.0)

NameError: name 'n' is not defined

Add a constant demand for transport

[5]:
n.add("Bus", "transport")
n.add("Load", "transport", bus="transport", p_set=1.0)


n.add("Bus", "diesel")


n.add("Store", "diesel", bus="diesel", e_cyclic=True, e_nom=1000.0)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [5], line 1
----> 1 n.add("Bus", "transport")
      2 n.add("Load", "transport", bus="transport", p_set=1.0)
      5 n.add("Bus", "diesel")

NameError: name 'n' is not defined

Add a bus for Hydrogen storage.

[6]:
n.add("Bus", "hydrogen")

n.add("Store", "hydrogen", bus="hydrogen", e_cyclic=True, e_nom=1000.0)

# n.add("Load","hydrogen",
#      bus="hydrogen",
#      p_set=1.)

n.add("Link", "electrolysis", p_nom=2.0, efficiency=0.8, bus0="bus", bus1="hydrogen")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [6], line 1
----> 1 n.add("Bus", "hydrogen")
      3 n.add("Store", "hydrogen", bus="hydrogen", e_cyclic=True, e_nom=1000.0)
      5 # n.add("Load","hydrogen",
      6 #      bus="hydrogen",
      7 #      p_set=1.)

NameError: name 'n' is not defined

Allow production of diesel from H2 and CO2 using Fischer-Tropsch

[7]:
n.add(
    "Link",
    "FT",
    p_nom=4,
    bus0="hydrogen",
    bus1="diesel",
    bus2="co2 stored",
    efficiency=1.0,
    efficiency2=-1,
)

# minus sign because opposite to how fossil fuels used:
# CH4 burning puts CH4 down, atmosphere up
n.add("Carrier", "co2", co2_emissions=-1.0)

# this tracks CO2 in the atmosphere
n.add("Bus", "co2 atmosphere", carrier="co2")

# NB: can also be negative
n.add("Store", "co2 atmosphere", e_nom=1000, e_min_pu=-1, bus="co2 atmosphere")

# this tracks CO2 stored, e.g. underground
n.add("Bus", "co2 stored")

# NB: can also be negative
n.add("Store", "co2 stored", e_nom=1000, e_min_pu=-1, bus="co2 stored")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [7], line 1
----> 1 n.add(
      2     "Link",
      3     "FT",
      4     p_nom=4,
      5     bus0="hydrogen",
      6     bus1="diesel",
      7     bus2="co2 stored",
      8     efficiency=1.0,
      9     efficiency2=-1,
     10 )
     12 # minus sign because opposite to how fossil fuels used:
     13 # CH4 burning puts CH4 down, atmosphere up
     14 n.add("Carrier", "co2", co2_emissions=-1.0)

NameError: name 'n' is not defined

Direct air capture consumes electricity to take CO2 from the air to the underground store

[8]:
n.add(
    "Link",
    "DAC",
    bus0="bus",
    bus1="co2 stored",
    bus2="co2 atmosphere",
    efficiency=1,
    efficiency2=-1,
    p_nom=5.0,
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [8], line 1
----> 1 n.add(
      2     "Link",
      3     "DAC",
      4     bus0="bus",
      5     bus1="co2 stored",
      6     bus2="co2 atmosphere",
      7     efficiency=1,
      8     efficiency2=-1,
      9     p_nom=5.0,
     10 )

NameError: name 'n' is not defined

Meet transport with diesel

[9]:
n.add(
    "Link",
    "diesel car",
    bus0="diesel",
    bus1="transport",
    bus2="co2 atmosphere",
    efficiency=1.0,
    efficiency2=1.0,
    p_nom=2.0,
)

n.add("Bus", "gas")

n.add("Store", "gas", e_initial=50, e_nom=50, marginal_cost=20, bus="gas")

n.add(
    "Link",
    "OCGT",
    bus0="gas",
    bus1="bus",
    bus2="co2 atmosphere",
    p_nom_extendable=True,
    efficiency=0.5,
    efficiency2=1,
)


n.add(
    "Link",
    "OCGT+CCS",
    bus0="gas",
    bus1="bus",
    bus2="co2 stored",
    bus3="co2 atmosphere",
    p_nom_extendable=True,
    efficiency=0.4,
    efficiency2=0.9,
    efficiency3=0.1,
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [9], line 1
----> 1 n.add(
      2     "Link",
      3     "diesel car",
      4     bus0="diesel",
      5     bus1="transport",
      6     bus2="co2 atmosphere",
      7     efficiency=1.0,
      8     efficiency2=1.0,
      9     p_nom=2.0,
     10 )
     12 n.add("Bus", "gas")
     14 n.add("Store", "gas", e_initial=50, e_nom=50, marginal_cost=20, bus="gas")

NameError: name 'n' is not defined

Add a cheap and a expensive biomass generator.

[10]:
biomass_marginal_cost = [20.0, 50.0]
biomass_stored = [40.0, 15.0]

for i in range(2):
    n.add("Bus", "biomass" + str(i))

    n.add(
        "Store",
        "biomass" + str(i),
        bus="biomass" + str(i),
        e_nom_extendable=True,
        marginal_cost=biomass_marginal_cost[i],
        e_nom=biomass_stored[i],
        e_initial=biomass_stored[i],
    )

    # simultaneously empties and refills co2 atmosphere
    n.add(
        "Link",
        "biomass" + str(i),
        bus0="biomass" + str(i),
        bus1="bus",
        p_nom_extendable=True,
        efficiency=0.5,
    )

    n.add(
        "Link",
        "biomass+CCS" + str(i),
        bus0="biomass" + str(i),
        bus1="bus",
        bus2="co2 stored",
        bus3="co2 atmosphere",
        p_nom_extendable=True,
        efficiency=0.4,
        efficiency2=1.0,
        efficiency3=-1,
    )


# can go to -50, but at some point can't generate enough electricity for DAC and demand
target = -50
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [10], line 5
      2 biomass_stored = [40.0, 15.0]
      4 for i in range(2):
----> 5     n.add("Bus", "biomass" + str(i))
      7     n.add(
      8         "Store",
      9         "biomass" + str(i),
   (...)
     14         e_initial=biomass_stored[i],
     15     )
     17     # simultaneously empties and refills co2 atmosphere

NameError: name 'n' is not defined

Add a global CO\(_2\) constraint.

[11]:
n.add(
    "GlobalConstraint",
    "co2_limit",
    sense="<=",
    carrier_attribute="co2_emissions",
    constant=target,
)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [11], line 1
----> 1 n.add(
      2     "GlobalConstraint",
      3     "co2_limit",
      4     sense="<=",
      5     carrier_attribute="co2_emissions",
      6     constant=target,
      7 )

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

NameError: name 'n' is not defined

How do the different stores in the system behave?

[13]:
n.stores_t.e.plot(figsize=(9, 7), lw=3)
plt.tight_layout()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [13], line 1
----> 1 n.stores_t.e.plot(figsize=(9, 7), lw=3)
      2 plt.tight_layout()

NameError: name 'n' is not defined

Let’s have a look at the production

[14]:
n.links_t.p0[["biomass+CCS0", "biomass+CCS1", "OCGT+CCS", "DAC"]].plot(
    subplots=True, figsize=(9, 7)
)
plt.tight_layout()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [14], line 1
----> 1 n.links_t.p0[["biomass+CCS0", "biomass+CCS1", "OCGT+CCS", "DAC"]].plot(
      2     subplots=True, figsize=(9, 7)
      3 )
      4 plt.tight_layout()

NameError: name 'n' is not defined

At all times, the amount of carbon is constant!

[15]:
n.stores_t.e[["co2 stored", "co2 atmosphere", "gas", "diesel"]].sum(axis=1)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [15], line 1
----> 1 n.stores_t.e[["co2 stored", "co2 atmosphere", "gas", "diesel"]].sum(axis=1)

NameError: name 'n' is not defined