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
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()
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)
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)
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
)
[7]:
network.lopf(network.snapshots, extra_functionality=extra_functionality)
network.objective
INFO:pypsa.linopf:Prepare linear problem
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In [7], line 1
----> 1 network.lopf(network.snapshots, extra_functionality=extra_functionality)
2 network.objective
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pypsa/components.py:770, in Network.lopf(self, snapshots, pyomo, solver_name, solver_options, solver_logfile, formulation, keep_files, extra_functionality, multi_investment_periods, **kwargs)
768 return network_lopf(self, **args)
769 else:
--> 770 return network_lopf_lowmem(self, **args)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pypsa/linopf.py:1456, in network_lopf(n, snapshots, solver_name, solver_logfile, extra_functionality, multi_investment_periods, skip_objective, skip_pre, extra_postprocessing, formulation, keep_references, keep_files, keep_shadowprices, solver_options, warmstart, store_basis, solver_dir)
1453 n.determine_network_topology()
1455 logger.info("Prepare linear problem")
-> 1456 fdp, problem_fn = prepare_lopf(
1457 n, snapshots, keep_files, skip_objective, extra_functionality, solver_dir
1458 )
1459 fds, solution_fn = mkstemp(prefix="pypsa-solve", suffix=".sol", dir=solver_dir)
1461 if warmstart == True:
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pypsa/linopf.py:1125, in prepare_lopf(n, snapshots, keep_files, skip_objective, extra_functionality, solver_dir)
1122 define_objective(n, snapshots)
1124 if extra_functionality is not None:
-> 1125 extra_functionality(n, snapshots)
1127 n.binaries_f.write("end\n")
1129 # explicit closing with file descriptor is necessary for windows machines
Cell In [6], line 10, in extra_functionality(network, snapshots)
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 )
17 # Guarantees c_m p_b1 \leq p_g1
18 def backpressure(model, snapshot):
AttributeError: 'Network' object has no attribute 'model'
Inspection
[8]:
network.loads_t.p
[8]:
| Load |
|---|
| snapshot |
| 2016-01-01 00:00:00 |
| 2016-01-01 01:00:00 |
| 2016-01-01 02:00:00 |
| 2016-01-01 03:00:00 |
[9]:
network.links.p_nom_opt
[9]:
Link
P2G 0.0
generator 0.0
boiler 0.0
Name: p_nom_opt, dtype: float64
[10]:
# CHP is dimensioned by the heat demand met in three hours when no wind
4 * 10.0 / 3 / network.links.at["boiler", "efficiency"]
[10]:
4.273504273504273
[11]:
# elec is set by the heat demand
28.490028 * 0.15
[11]:
4.2735042
[12]:
network.links_t.p0
[12]:
| Link |
|---|
| snapshot |
| 2016-01-01 00:00:00 |
| 2016-01-01 01:00:00 |
| 2016-01-01 02:00:00 |
| 2016-01-01 03:00:00 |
[13]:
network.links_t.p1
[13]:
| Link |
|---|
| snapshot |
| 2016-01-01 00:00:00 |
| 2016-01-01 01:00:00 |
| 2016-01-01 02:00:00 |
| 2016-01-01 03:00:00 |
[14]:
pd.DataFrame({attr: network.stores_t[attr]["gas depot"] for attr in ["p", "e"]})
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/core/indexes/base.py:3803, in Index.get_loc(self, key, method, tolerance)
3802 try:
-> 3803 return self._engine.get_loc(casted_key)
3804 except KeyError as err:
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/_libs/index.pyx:138, in pandas._libs.index.IndexEngine.get_loc()
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/_libs/index.pyx:165, in pandas._libs.index.IndexEngine.get_loc()
File pandas/_libs/hashtable_class_helper.pxi:5745, in pandas._libs.hashtable.PyObjectHashTable.get_item()
File pandas/_libs/hashtable_class_helper.pxi:5753, in pandas._libs.hashtable.PyObjectHashTable.get_item()
KeyError: 'gas depot'
The above exception was the direct cause of the following exception:
KeyError 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"]})
Cell In [14], line 1, in <dictcomp>(.0)
----> 1 pd.DataFrame({attr: network.stores_t[attr]["gas depot"] for attr in ["p", "e"]})
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/core/frame.py:3804, in DataFrame.__getitem__(self, key)
3802 if self.columns.nlevels > 1:
3803 return self._getitem_multilevel(key)
-> 3804 indexer = self.columns.get_loc(key)
3805 if is_integer(indexer):
3806 indexer = [indexer]
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/core/indexes/base.py:3805, in Index.get_loc(self, key, method, tolerance)
3803 return self._engine.get_loc(casted_key)
3804 except KeyError as err:
-> 3805 raise KeyError(key) from err
3806 except TypeError:
3807 # If we have a listlike key, _check_indexing_error will raise
3808 # InvalidIndexError. Otherwise we fall through and re-raise
3809 # the TypeError.
3810 self._check_indexing_error(key)
KeyError: 'gas depot'
[15]:
pd.DataFrame({attr: network.stores_t[attr]["water tank"] for attr in ["p", "e"]})
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/core/indexes/base.py:3803, in Index.get_loc(self, key, method, tolerance)
3802 try:
-> 3803 return self._engine.get_loc(casted_key)
3804 except KeyError as err:
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/_libs/index.pyx:138, in pandas._libs.index.IndexEngine.get_loc()
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/_libs/index.pyx:165, in pandas._libs.index.IndexEngine.get_loc()
File pandas/_libs/hashtable_class_helper.pxi:5745, in pandas._libs.hashtable.PyObjectHashTable.get_item()
File pandas/_libs/hashtable_class_helper.pxi:5753, in pandas._libs.hashtable.PyObjectHashTable.get_item()
KeyError: 'water tank'
The above exception was the direct cause of the following exception:
KeyError 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"]})
Cell In [15], line 1, in <dictcomp>(.0)
----> 1 pd.DataFrame({attr: network.stores_t[attr]["water tank"] for attr in ["p", "e"]})
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/core/frame.py:3804, in DataFrame.__getitem__(self, key)
3802 if self.columns.nlevels > 1:
3803 return self._getitem_multilevel(key)
-> 3804 indexer = self.columns.get_loc(key)
3805 if is_integer(indexer):
3806 indexer = [indexer]
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/core/indexes/base.py:3805, in Index.get_loc(self, key, method, tolerance)
3803 return self._engine.get_loc(casted_key)
3804 except KeyError as err:
-> 3805 raise KeyError(key) from err
3806 except TypeError:
3807 # If we have a listlike key, _check_indexing_error will raise
3808 # InvalidIndexError. Otherwise we fall through and re-raise
3809 # the TypeError.
3810 self._check_indexing_error(key)
KeyError: 'water tank'
[16]:
pd.DataFrame({attr: network.links_t[attr]["boiler"] for attr in ["p0", "p1"]})
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/core/indexes/base.py:3803, in Index.get_loc(self, key, method, tolerance)
3802 try:
-> 3803 return self._engine.get_loc(casted_key)
3804 except KeyError as err:
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/_libs/index.pyx:138, in pandas._libs.index.IndexEngine.get_loc()
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/_libs/index.pyx:165, in pandas._libs.index.IndexEngine.get_loc()
File pandas/_libs/hashtable_class_helper.pxi:5745, in pandas._libs.hashtable.PyObjectHashTable.get_item()
File pandas/_libs/hashtable_class_helper.pxi:5753, in pandas._libs.hashtable.PyObjectHashTable.get_item()
KeyError: 'boiler'
The above exception was the direct cause of the following exception:
KeyError Traceback (most recent call last)
Cell In [16], line 1
----> 1 pd.DataFrame({attr: network.links_t[attr]["boiler"] for attr in ["p0", "p1"]})
Cell In [16], line 1, in <dictcomp>(.0)
----> 1 pd.DataFrame({attr: network.links_t[attr]["boiler"] for attr in ["p0", "p1"]})
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/core/frame.py:3804, in DataFrame.__getitem__(self, key)
3802 if self.columns.nlevels > 1:
3803 return self._getitem_multilevel(key)
-> 3804 indexer = self.columns.get_loc(key)
3805 if is_integer(indexer):
3806 indexer = [indexer]
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.1/lib/python3.10/site-packages/pandas/core/indexes/base.py:3805, in Index.get_loc(self, key, method, tolerance)
3803 return self._engine.get_loc(casted_key)
3804 except KeyError as err:
-> 3805 raise KeyError(key) from err
3806 except TypeError:
3807 # If we have a listlike key, _check_indexing_error will raise
3808 # InvalidIndexError. Otherwise we fall through and re-raise
3809 # the TypeError.
3810 self._check_indexing_error(key)
KeyError: 'boiler'
[17]:
network.stores.loc["gas depot"]
[17]:
attribute
bus 0 gas
type
carrier gas
e_nom 0.0
e_nom_extendable True
e_nom_min 0.0
e_nom_max inf
e_min_pu 0.0
e_max_pu 1.0
e_initial 0.0
e_initial_per_period False
e_cyclic True
e_cyclic_per_period True
p_set 0.0
q_set 0.0
sign 1.0
marginal_cost 0.0
capital_cost 0.0
standing_loss 0.0
build_year 0
lifetime inf
e_nom_opt 0.0
Name: gas depot, dtype: object
[18]:
network.generators.loc["wind turbine"]
[18]:
attribute
bus 0
control Slack
type
p_nom 0.0
p_nom_extendable True
p_nom_min 0.0
p_nom_max inf
p_min_pu 0.0
p_max_pu 1.0
p_set 0.0
q_set 0.0
sign 1.0
carrier wind
marginal_cost 0.0
build_year 0
lifetime inf
capital_cost 1000.0
efficiency 1.0
committable False
start_up_cost 0.0
shut_down_cost 0.0
min_up_time 0
min_down_time 0
up_time_before 1
down_time_before 0
ramp_limit_up NaN
ramp_limit_down NaN
ramp_limit_start_up 1.0
ramp_limit_shut_down 1.0
p_nom_opt 0.0
Name: wind turbine, dtype: object
[19]:
network.links.p_nom_opt
[19]:
Link
P2G 0.0
generator 0.0
boiler 0.0
Name: p_nom_opt, dtype: float64
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))
[20]:
0.9099999999999999