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

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.,"2nd bus efficiency","Input (optional)"]
override_component_attrs["Link"].loc["efficiency3"] = ["static or series","per unit",1.,"3rd bus efficiency","Input (optional)"]
override_component_attrs["Link"].loc["p2"] = ["series","MW",0.,"2nd bus output","Output"]
override_component_attrs["Link"].loc["p3"] = ["series","MW",0.,"3rd bus output","Output"]
[3]:
n = pypsa.Network(override_component_attrs=override_component_attrs)
n.set_snapshots(range(10))

Add a constant electrical load

[4]:
n.add("Bus","bus")
n.add("Load","load",bus="bus",
      p_set=1.)

Add a constant demand for transport

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


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


n.add("Store","diesel",bus="diesel",
      e_cyclic=True,
      e_nom=1000.)

Add a bus for Hydrogen storage.

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

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

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

n.add("Link","electrolysis",
      p_nom=2.,
      efficiency=0.8,
      bus0="bus",
      bus1="hydrogen")

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.,
      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.)

#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")

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.)

Meet transport with diesel

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

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)

Add a cheap and a expensive biomass generator.

[10]:
biomass_marginal_cost = [20.,50.]
biomass_stored = [40.,15.]

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.,
          efficiency3 = -1)


#can go to -50, but at some point can't generate enough electricity for DAC and demand
target = -50

Add a global CO\(_2\) constraint.

[11]:
n.add("GlobalConstraint","co2_limit",
      sense="<=",
      carrier_attribute="co2_emissions",
      constant=target)
[12]:
n.lopf();
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `kirchhoff` formulation
INFO:pypsa.opf:Solving model using glpk
WARNING:pyomo.solvers:Could not locate the 'glpsol' executable, which is required for solver 'glpk'
---------------------------------------------------------------------------
ApplicationError                          Traceback (most recent call last)
/tmp/ipykernel_240/1812968345.py in <module>
----> 1 n.lopf();

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pypsa/components.py in lopf(self, snapshots, pyomo, solver_name, solver_options, solver_logfile, formulation, keep_files, extra_functionality, multi_investment_periods, **kwargs)
    646
    647         if pyomo:
--> 648             return network_lopf(self, **args)
    649         else:
    650             return network_lopf_lowmem(self, **args)

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pypsa/opf.py in network_lopf(network, snapshots, solver_name, solver_io, skip_pre, extra_functionality, multi_investment_periods, solver_logfile, solver_options, keep_files, formulation, ptdf_tolerance, free_memory, extra_postprocessing)
   1663                               solver_logfile=solver_logfile, solver_options=solver_options,
   1664                               keep_files=keep_files, free_memory=free_memory,
-> 1665                               extra_postprocessing=extra_postprocessing)

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pypsa/opf.py in network_lopf_solve(network, snapshots, formulation, solver_options, solver_logfile, keep_files, free_memory, extra_postprocessing)
   1563             network.results = network.opt.solve(*args, suffixes=["dual"], keepfiles=keep_files, logfile=solver_logfile, options=solver_options)
   1564     else:
-> 1565         network.results = network.opt.solve(*args, suffixes=["dual"], keepfiles=keep_files, logfile=solver_logfile, options=solver_options)
   1566
   1567     if logger.isEnabledFor(logging.INFO):

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pyomo/opt/base/solvers.py in solve(self, *args, **kwds)
    510         """ Solve the problem """
    511
--> 512         self.available(exception_flag=True)
    513         #
    514         # If the inputs are models, then validate that they have been

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pyomo/opt/solver/shellcmd.py in available(self, exception_flag)
    123             if exception_flag:
    124                 msg = "No executable found for solver '%s'"
--> 125                 raise ApplicationError(msg % self.name)
    126             return False
    127         return True

ApplicationError: No executable found for solver 'glpk'

How do the different stores in the system behave?

[13]:
n.stores_t.e.plot(figsize = (9,7), lw=3)
plt.tight_layout()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/tmp/ipykernel_240/1715769439.py in <module>
----> 1 n.stores_t.e.plot(figsize = (9,7), lw=3)
      2 plt.tight_layout()

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pandas/plotting/_core.py in __call__(self, *args, **kwargs)
    970                     data.columns = label_name
    971
--> 972         return plot_backend.plot(data, kind=kind, **kwargs)
    973
    974     __call__.__doc__ = __doc__

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pandas/plotting/_matplotlib/__init__.py in plot(data, kind, **kwargs)
     69             kwargs["ax"] = getattr(ax, "left_ax", ax)
     70     plot_obj = PLOT_CLASSES[kind](data, **kwargs)
---> 71     plot_obj.generate()
     72     plot_obj.draw()
     73     return plot_obj.result

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pandas/plotting/_matplotlib/core.py in generate(self)
    284     def generate(self):
    285         self._args_adjust()
--> 286         self._compute_plot_data()
    287         self._setup_subplots()
    288         self._make_plot()

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pandas/plotting/_matplotlib/core.py in _compute_plot_data(self)
    451         # no non-numeric frames or series allowed
    452         if is_empty:
--> 453             raise TypeError("no numeric data to plot")
    454
    455         self.data = numeric_data.apply(self._convert_to_ndarray)

TypeError: no numeric data to plot

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()
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
/tmp/ipykernel_240/447846780.py in <module>
----> 1 n.links_t.p0[["biomass+CCS0", "biomass+CCS1", "OCGT+CCS", "DAC"]].plot(subplots=True, figsize = (9,7))
      2 plt.tight_layout()

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pandas/core/frame.py in __getitem__(self, key)
   3462             if is_iterator(key):
   3463                 key = list(key)
-> 3464             indexer = self.loc._get_listlike_indexer(key, axis=1)[1]
   3465
   3466         # take() does not accept boolean indexers

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pandas/core/indexing.py in _get_listlike_indexer(self, key, axis)
   1312             keyarr, indexer, new_indexer = ax._reindex_non_unique(keyarr)
   1313
-> 1314         self._validate_read_indexer(keyarr, indexer, axis)
   1315
   1316         if needs_i8_conversion(ax.dtype) or isinstance(

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pandas/core/indexing.py in _validate_read_indexer(self, key, indexer, axis)
   1372                 if use_interval_msg:
   1373                     key = list(key)
-> 1374                 raise KeyError(f"None of [{key}] are in the [{axis_name}]")
   1375
   1376             not_found = list(ensure_index(key)[missing_mask.nonzero()[0]].unique())

KeyError: "None of [Index(['biomass+CCS0', 'biomass+CCS1', 'OCGT+CCS', 'DAC'], dtype='object')] are in the [columns]"

At all times, the amount of carbon is constant!

[15]:
n.stores_t.e[["co2 stored","co2 atmosphere","gas","diesel"]].sum(axis=1)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
/tmp/ipykernel_240/3838290976.py in <module>
----> 1 n.stores_t.e[["co2 stored","co2 atmosphere","gas","diesel"]].sum(axis=1)

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pandas/core/frame.py in __getitem__(self, key)
   3462             if is_iterator(key):
   3463                 key = list(key)
-> 3464             indexer = self.loc._get_listlike_indexer(key, axis=1)[1]
   3465
   3466         # take() does not accept boolean indexers

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pandas/core/indexing.py in _get_listlike_indexer(self, key, axis)
   1312             keyarr, indexer, new_indexer = ax._reindex_non_unique(keyarr)
   1313
-> 1314         self._validate_read_indexer(keyarr, indexer, axis)
   1315
   1316         if needs_i8_conversion(ax.dtype) or isinstance(

~/checkouts/readthedocs.org/user_builds/pypsa-docs-staging/envs/latest/lib/python3.7/site-packages/pandas/core/indexing.py in _validate_read_indexer(self, key, indexer, axis)
   1372                 if use_interval_msg:
   1373                     key = list(key)
-> 1374                 raise KeyError(f"None of [{key}] are in the [{axis_name}]")
   1375
   1376             not_found = list(ensure_index(key)[missing_mask.nonzero()[0]].unique())

KeyError: "None of [Index(['co2 stored', 'co2 atmosphere', 'gas', 'diesel'], dtype='object')] are in the [columns]"