Note

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

Power to Gas with Heat Coupling

This is an example for power to gas with optional coupling to heat sector (via boiler OR Combined-Heat-and-Power (CHP))

A location has an electric, gas and heat bus. The primary source is wind power, which can be converted to gas. The gas can be stored to convert into electricity or heat (with either a boiler or a CHP).

[1]:
import pypsa
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pyomo.environ import Constraint

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

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

Combined-Heat-and-Power (CHP) parameterisation

This setup follows http://www.ea-energianalyse.dk/reports/student-reports/integration_of_50_percent_wind%20power.pdf pages 35-6 which follows http://www.sciencedirect.com/science/article/pii/030142159390282K

[2]:
# ratio between max heat output and max electric output
nom_r = 1.0

# backpressure limit
c_m = 0.75

# marginal loss for each additional generation of heat
c_v = 0.15

Graph for the case that max heat output equals max electric output

[3]:
fig, ax = plt.subplots(figsize=(9, 5))

t = 0.01
ph = np.arange(0, 1.0001, t)

ax.plot(ph, c_m * ph)
ax.set_xlabel("P_heat_out")
ax.set_ylabel("P_elec_out")
ax.grid(True)

ax.set_xlim([0, 1.1])
ax.set_ylim([0, 1.1])
ax.text(0.1, 0.7, "Allowed output", color="r")
ax.plot(ph, 1 - c_v * ph)

for i in range(1, 10):
    k = 0.1 * i
    x = np.arange(0, k / (c_m + c_v), t)
    ax.plot(x, k - c_v * x, color="g", alpha=0.5)

ax.text(0.05, 0.41, "iso-fuel-lines", color="g", rotation=-7)
ax.fill_between(ph, c_m * ph, 1 - c_v * ph, facecolor="r", alpha=0.5)

fig.tight_layout()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [3], line 1
----> 1 fig, ax = plt.subplots(figsize=(9, 5))
      3 t = 0.01
      4 ph = np.arange(0, 1.0001, t)

NameError: name 'plt' is not defined

Optimisation

[4]:
network = pypsa.Network()
network.set_snapshots(pd.date_range("2016-01-01 00:00", "2016-01-01 03:00", freq="H"))

network.add("Bus", "0", carrier="AC")
network.add("Bus", "0 gas", carrier="gas")

network.add("Carrier", "wind")
network.add("Carrier", "gas", co2_emissions=0.2)

network.add("GlobalConstraint", "co2_limit", sense="<=", constant=0.0)

network.add(
    "Generator",
    "wind turbine",
    bus="0",
    carrier="wind",
    p_nom_extendable=True,
    p_max_pu=[0.0, 0.2, 0.7, 0.4],
    capital_cost=1000,
)

network.add("Load", "load", bus="0", p_set=5.0)

network.add(
    "Link",
    "P2G",
    bus0="0",
    bus1="0 gas",
    efficiency=0.6,
    capital_cost=1000,
    p_nom_extendable=True,
)

network.add(
    "Link",
    "generator",
    bus0="0 gas",
    bus1="0",
    efficiency=0.468,
    capital_cost=400,
    p_nom_extendable=True,
)

network.add("Store", "gas depot", bus="0 gas", e_cyclic=True, e_nom_extendable=True)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [4], line 1
----> 1 network = pypsa.Network()
      2 network.set_snapshots(pd.date_range("2016-01-01 00:00", "2016-01-01 03:00", freq="H"))
      4 network.add("Bus", "0", carrier="AC")

NameError: name 'pypsa' is not defined

Add heat sector

[5]:
network.add("Bus", "0 heat", carrier="heat")

network.add("Carrier", "heat")

network.add("Load", "heat load", bus="0 heat", p_set=10.0)

network.add(
    "Link",
    "boiler",
    bus0="0 gas",
    bus1="0 heat",
    efficiency=0.9,
    capital_cost=300,
    p_nom_extendable=True,
)

network.add("Store", "water tank", bus="0 heat", e_cyclic=True, e_nom_extendable=True)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [5], line 1
----> 1 network.add("Bus", "0 heat", carrier="heat")
      3 network.add("Carrier", "heat")
      5 network.add("Load", "heat load", bus="0 heat", p_set=10.0)

NameError: name 'network' is not defined

Add CHP constraints

[6]:
# Guarantees ISO fuel lines, i.e. fuel consumption p_b0 + p_g0 = constant along p_g1 + c_v p_b1 = constant
network.links.at["boiler", "efficiency"] = (
    network.links.at["generator", "efficiency"] / c_v
)


def extra_functionality(network, snapshots):

    # Guarantees heat output and electric output nominal powers are proportional
    network.model.chp_nom = Constraint(
        rule=lambda model: network.links.at["generator", "efficiency"]
        * nom_r
        * model.link_p_nom["generator"]
        == network.links.at["boiler", "efficiency"] * model.link_p_nom["boiler"]
    )

    # Guarantees c_m p_b1  \leq p_g1
    def backpressure(model, snapshot):
        return (
            c_m
            * network.links.at["boiler", "efficiency"]
            * model.link_p["boiler", snapshot]
            <= network.links.at["generator", "efficiency"]
            * model.link_p["generator", snapshot]
        )

    network.model.backpressure = Constraint(list(snapshots), rule=backpressure)

    # Guarantees p_g1 +c_v p_b1 \leq p_g1_nom
    def top_iso_fuel_line(model, snapshot):
        return (
            model.link_p["boiler", snapshot] + model.link_p["generator", snapshot]
            <= model.link_p_nom["generator"]
        )

    network.model.top_iso_fuel_line = Constraint(
        list(snapshots), rule=top_iso_fuel_line
    )
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [6], line 3
      1 # Guarantees ISO fuel lines, i.e. fuel consumption p_b0 + p_g0 = constant along p_g1 + c_v p_b1 = constant
      2 network.links.at["boiler", "efficiency"] = (
----> 3     network.links.at["generator", "efficiency"] / c_v
      4 )
      7 def extra_functionality(network, snapshots):
      8
      9     # Guarantees heat output and electric output nominal powers are proportional
     10     network.model.chp_nom = Constraint(
     11         rule=lambda model: network.links.at["generator", "efficiency"]
     12         * nom_r
     13         * model.link_p_nom["generator"]
     14         == network.links.at["boiler", "efficiency"] * model.link_p_nom["boiler"]
     15     )

NameError: name 'network' is not defined
[7]:
network.lopf(network.snapshots, extra_functionality=extra_functionality)
network.objective
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [7], line 1
----> 1 network.lopf(network.snapshots, extra_functionality=extra_functionality)
      2 network.objective

NameError: name 'network' is not defined

Inspection

[8]:
network.loads_t.p
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [8], line 1
----> 1 network.loads_t.p

NameError: name 'network' is not defined
[9]:
network.links.p_nom_opt
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [9], line 1
----> 1 network.links.p_nom_opt

NameError: name 'network' is not defined
[10]:
# CHP is dimensioned by the heat demand met in three hours when no wind
4 * 10.0 / 3 / network.links.at["boiler", "efficiency"]
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [10], line 2
      1 # CHP is dimensioned by the heat demand met in three hours when no wind
----> 2 4 * 10.0 / 3 / network.links.at["boiler", "efficiency"]

NameError: name 'network' is not defined
[11]:
# elec is set by the heat demand
28.490028 * 0.15
[11]:
4.2735042
[12]:
network.links_t.p0
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [12], line 1
----> 1 network.links_t.p0

NameError: name 'network' is not defined
[13]:
network.links_t.p1
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [13], line 1
----> 1 network.links_t.p1

NameError: name 'network' is not defined
[14]:
pd.DataFrame({attr: network.stores_t[attr]["gas depot"] for attr in ["p", "e"]})
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [14], line 1
----> 1 pd.DataFrame({attr: network.stores_t[attr]["gas depot"] for attr in ["p", "e"]})

NameError: name 'pd' is not defined
[15]:
pd.DataFrame({attr: network.stores_t[attr]["water tank"] for attr in ["p", "e"]})
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [15], line 1
----> 1 pd.DataFrame({attr: network.stores_t[attr]["water tank"] for attr in ["p", "e"]})

NameError: name 'pd' is not defined
[16]:
pd.DataFrame({attr: network.links_t[attr]["boiler"] for attr in ["p0", "p1"]})
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [16], line 1
----> 1 pd.DataFrame({attr: network.links_t[attr]["boiler"] for attr in ["p0", "p1"]})

NameError: name 'pd' is not defined
[17]:
network.stores.loc["gas depot"]
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [17], line 1
----> 1 network.stores.loc["gas depot"]

NameError: name 'network' is not defined
[18]:
network.generators.loc["wind turbine"]
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [18], line 1
----> 1 network.generators.loc["wind turbine"]

NameError: name 'network' is not defined
[19]:
network.links.p_nom_opt
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [19], line 1
----> 1 network.links.p_nom_opt

NameError: name 'network' is not defined

Calculate the overall efficiency of the CHP

[20]:
eta_elec = network.links.at["generator", "efficiency"]

r = 1 / c_m

# P_h = r*P_e
(1 + r) / ((1 / eta_elec) * (1 + c_v * r))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [20], line 1
----> 1 eta_elec = network.links.at["generator", "efficiency"]
      3 r = 1 / c_m
      5 # P_h = r*P_e

NameError: name 'network' is not defined