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
WARNING:pypsa.components:Solving optimisation problem with pyomo.In PyPSA version 0.21 the default will change to ``n.lopf(pyomo=False)``.Explicitly set ``n.lopf(pyomo=True)`` to retain current behaviour.
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `kirchhoff` formulation
---------------------------------------------------------------------------
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.0/lib/python3.10/site-packages/pypsa/components.py:773, in Network.lopf(self, snapshots, pyomo, solver_name, solver_options, solver_logfile, formulation, keep_files, extra_functionality, multi_investment_periods, **kwargs)
767 if pyomo:
768 logger.warning(
769 "Solving optimisation problem with pyomo."
770 "In PyPSA version 0.21 the default will change to ``n.lopf(pyomo=False)``."
771 "Explicitly set ``n.lopf(pyomo=True)`` to retain current behaviour."
772 )
--> 773 return network_lopf(self, **args)
774 else:
775 return network_lopf_lowmem(self, **args)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.0/lib/python3.10/site-packages/pypsa/opf.py:2416, 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)
2410 logger.warning(
2411 "Encountered nonzero ramp limits for links. These are ignored when running the optimization with `pyomo=True`."
2412 )
2414 snapshots = _as_snapshots(network, snapshots)
-> 2416 network_lopf_build_model(
2417 network,
2418 snapshots,
2419 skip_pre=skip_pre,
2420 formulation=formulation,
2421 ptdf_tolerance=ptdf_tolerance,
2422 )
2424 if extra_functionality is not None:
2425 extra_functionality(network, snapshots)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.0/lib/python3.10/site-packages/pypsa/opf.py:2182, in network_lopf_build_model(network, snapshots, skip_pre, formulation, ptdf_tolerance)
2179 elif formulation in ["ptdf", "cycles"]:
2180 define_sub_network_balance_constraints(network, snapshots)
-> 2182 define_global_constraints(network, snapshots)
2184 define_linear_objective(network, snapshots)
2186 # tidy up auxilliary expressions
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.0/lib/python3.10/site-packages/pypsa/opf.py:1725, in define_global_constraints(network, snapshots)
1718 c.lhs.constant += sum(
1719 attribute * network.stores.at[store, "e_initial"]
1720 for store in stores
1721 )
1723 global_constraints[gc] = c
-> 1725 l_constraint(
1726 network.model,
1727 "global_constraints",
1728 global_constraints,
1729 list(network.global_constraints.index),
1730 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.0/lib/python3.10/site-packages/pypsa/opt.py:219, in l_constraint(model, name, constraints, *args)
214 else:
215 raise KeyError(
216 '`sense` must be one of "==","<=",">=","><"; got: {}'.format(sense)
217 )
--> 219 v._data[i] = _GeneralConstraintData(constr_expr, v)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.0/lib/python3.10/site-packages/pyomo/core/base/constraint.py:299, in _GeneralConstraintData.__init__(self, expr, component)
297 self._expr = None
298 if expr is not None:
--> 299 self.set_value(expr)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.0/lib/python3.10/site-packages/pyomo/core/base/constraint.py:556, in _GeneralConstraintData.set_value(self, expr)
541 raise ValueError(
542 "Invalid constraint expression. The constraint "
543 "expression resolved to a trivial Boolean (%s) "
(...)
547 % (expr, "Feasible" if expr else "Infeasible",
548 expr, self.name))
550 else:
551 msg = ("Constraint '%s' does not have a proper "
552 "value. Found '%s'\nExpecting a tuple or "
553 "equation. Examples:"
554 "\n sum(model.costs) == model.income"
555 "\n (0, model.price[item], 50)"
--> 556 % (self.name, str(expr)))
557 raise ValueError(msg)
558 #
559 # Normalize the incoming expressions, if we can
560 #
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.0/lib/python3.10/site-packages/pyomo/core/base/component.py:316, in _ComponentBase.name(self)
313 @property
314 def name(self):
315 """Get the fully qualifed component name."""
--> 316 return self.getname(fully_qualified=True)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.0/lib/python3.10/site-packages/pyomo/core/base/component.py:965, in ComponentData.getname(self, fully_qualified, name_buffer, relative_to)
959 return name_buffer[id(self)]
960 else:
961 #
962 # No buffer, we can do what we are going to do all the time after we
963 # deprecate the buffer.
964 #
--> 965 return base + index_repr(self.index())
966 #
967 raise RuntimeError("Fatal error: cannot find the component data in "
968 "the owning component's _data dictionary.")
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.0/lib/python3.10/site-packages/pyomo/core/base/component.py:891, in ComponentData.index(self)
882 """
883 Returns the index of this ComponentData instance relative
884 to the parent component index set. None is returned if
(...)
887 to the parent component's index set.
888 """
889 parent = self.parent_component()
890 if ( parent is not None and
--> 891 self._index is not NOTSET and
892 parent[self._index] is not self ):
893 # This error message is a bit goofy, but we can't call self.name
894 # here--it's an infinite loop!
895 raise DeveloperError(
896 "The '_data' dictionary and '_index' attribute are out of "
897 "sync for indexed %s '%s': The %s entry in the '_data' "
898 "dictionary does not map back to this component data object."
899 % (parent.ctype.__name__, parent.name, self._index))
900 return self._index
AttributeError: '_GeneralConstraintData' object has no attribute '_index'
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.0/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.0/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.0/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.0/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.0/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.0/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.0/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.0/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.0/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.0/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.0/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.0/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.0/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.0/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.0/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