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)
Input In [7], in <cell line: 1>()
----> 1 network.lopf(network.snapshots, extra_functionality=extra_functionality)
2 network.objective
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pypsa/components.py:750, in Network.lopf(self, snapshots, pyomo, solver_name, solver_options, solver_logfile, formulation, keep_files, extra_functionality, multi_investment_periods, **kwargs)
744 if pyomo:
745 logger.warning(
746 "Solving optimisation problem with pyomo."
747 "In PyPSA version 0.21 the default will change to ``n.lopf(pyomo=False)``."
748 "Explicitly set ``n.lopf(pyomo=True)`` to retain current behaviour."
749 )
--> 750 return network_lopf(self, **args)
751 else:
752 return network_lopf_lowmem(self, **args)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pypsa/opf.py:2428, 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)
2422 logger.warning(
2423 "Encountered nonzero ramp limits for links. These are ignored when running the optimization with `pyomo=True`."
2424 )
2426 snapshots = _as_snapshots(network, snapshots)
-> 2428 network_lopf_build_model(
2429 network,
2430 snapshots,
2431 skip_pre=skip_pre,
2432 formulation=formulation,
2433 ptdf_tolerance=ptdf_tolerance,
2434 )
2436 if extra_functionality is not None:
2437 extra_functionality(network, snapshots)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pypsa/opf.py:2195, in network_lopf_build_model(network, snapshots, skip_pre, formulation, ptdf_tolerance)
2192 elif formulation in ["ptdf", "cycles"]:
2193 define_sub_network_balance_constraints(network, snapshots)
-> 2195 define_global_constraints(network, snapshots)
2197 define_linear_objective(network, snapshots)
2199 # tidy up auxilliary expressions
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pypsa/opf.py:1737, in define_global_constraints(network, snapshots)
1730 c.lhs.constant += sum(
1731 attribute * network.stores.at[store, "e_initial"]
1732 for store in stores
1733 )
1735 global_constraints[gc] = c
-> 1737 l_constraint(
1738 network.model,
1739 "global_constraints",
1740 global_constraints,
1741 list(network.global_constraints.index),
1742 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.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.20.0/lib/python3.10/site-packages/pyomo/core/base/constraint.py:297, in _GeneralConstraintData.__init__(self, expr, component)
295 self._expr = None
296 if expr is not None:
--> 297 self.set_value(expr)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pyomo/core/base/constraint.py:546, in _GeneralConstraintData.set_value(self, expr)
531 raise ValueError(
532 "Invalid constraint expression. The constraint "
533 "expression resolved to a trivial Boolean (%s) "
(...)
537 % (expr, "Feasible" if expr else "Infeasible",
538 expr, self.name))
540 else:
541 msg = ("Constraint '%s' does not have a proper "
542 "value. Found '%s'\nExpecting a tuple or "
543 "equation. Examples:"
544 "\n sum(model.costs) == model.income"
545 "\n (0, model.price[item], 50)"
--> 546 % (self.name, str(expr)))
547 raise ValueError(msg)
548 #
549 # Normalize the incoming expressions, if we can
550 #
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pyomo/core/base/component.py:281, in _ComponentBase.name(self)
278 @property
279 def name(self):
280 """Get the fully qualifed component name."""
--> 281 return self.getname(fully_qualified=True)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pyomo/core/base/component.py:930, in ComponentData.getname(self, fully_qualified, name_buffer, relative_to)
924 return name_buffer[id(self)]
925 else:
926 #
927 # No buffer, we can do what we are going to do all the time after we
928 # deprecate the buffer.
929 #
--> 930 return base + index_repr(self.index())
931 #
932 raise RuntimeError("Fatal error: cannot find the component data in "
933 "the owning component's _data dictionary.")
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pyomo/core/base/component.py:856, in ComponentData.index(self)
847 """
848 Returns the index of this ComponentData instance relative
849 to the parent component index set. None is returned if
(...)
852 to the parent component's index set.
853 """
854 parent = self.parent_component()
855 if ( parent is not None and
--> 856 self._index is not NOTSET and
857 parent[self._index] is not self ):
858 # This error message is a bit goofy, but we can't call self.name
859 # here--it's an infinite loop!
860 raise DeveloperError(
861 "The '_data' dictionary and '_index' attribute are out of "
862 "sync for indexed %s '%s': The %s entry in the '_data' "
863 "dictionary does not map back to this component data object."
864 % (parent.ctype.__name__, parent.name, self._index))
865 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.20.0/lib/python3.10/site-packages/pandas/core/indexes/base.py:3621, in Index.get_loc(self, key, method, tolerance)
3620 try:
-> 3621 return self._engine.get_loc(casted_key)
3622 except KeyError as err:
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pandas/_libs/index.pyx:136, in pandas._libs.index.IndexEngine.get_loc()
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pandas/_libs/index.pyx:163, in pandas._libs.index.IndexEngine.get_loc()
File pandas/_libs/hashtable_class_helper.pxi:5198, in pandas._libs.hashtable.PyObjectHashTable.get_item()
File pandas/_libs/hashtable_class_helper.pxi:5206, 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)
Input In [14], in <cell line: 1>()
----> 1 pd.DataFrame({attr: network.stores_t[attr]["gas depot"] for attr in ["p", "e"]})
Input In [14], 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.20.0/lib/python3.10/site-packages/pandas/core/frame.py:3505, in DataFrame.__getitem__(self, key)
3503 if self.columns.nlevels > 1:
3504 return self._getitem_multilevel(key)
-> 3505 indexer = self.columns.get_loc(key)
3506 if is_integer(indexer):
3507 indexer = [indexer]
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pandas/core/indexes/base.py:3623, in Index.get_loc(self, key, method, tolerance)
3621 return self._engine.get_loc(casted_key)
3622 except KeyError as err:
-> 3623 raise KeyError(key) from err
3624 except TypeError:
3625 # If we have a listlike key, _check_indexing_error will raise
3626 # InvalidIndexError. Otherwise we fall through and re-raise
3627 # the TypeError.
3628 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.20.0/lib/python3.10/site-packages/pandas/core/indexes/base.py:3621, in Index.get_loc(self, key, method, tolerance)
3620 try:
-> 3621 return self._engine.get_loc(casted_key)
3622 except KeyError as err:
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pandas/_libs/index.pyx:136, in pandas._libs.index.IndexEngine.get_loc()
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pandas/_libs/index.pyx:163, in pandas._libs.index.IndexEngine.get_loc()
File pandas/_libs/hashtable_class_helper.pxi:5198, in pandas._libs.hashtable.PyObjectHashTable.get_item()
File pandas/_libs/hashtable_class_helper.pxi:5206, 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)
Input In [15], in <cell line: 1>()
----> 1 pd.DataFrame({attr: network.stores_t[attr]["water tank"] for attr in ["p", "e"]})
Input In [15], 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.20.0/lib/python3.10/site-packages/pandas/core/frame.py:3505, in DataFrame.__getitem__(self, key)
3503 if self.columns.nlevels > 1:
3504 return self._getitem_multilevel(key)
-> 3505 indexer = self.columns.get_loc(key)
3506 if is_integer(indexer):
3507 indexer = [indexer]
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pandas/core/indexes/base.py:3623, in Index.get_loc(self, key, method, tolerance)
3621 return self._engine.get_loc(casted_key)
3622 except KeyError as err:
-> 3623 raise KeyError(key) from err
3624 except TypeError:
3625 # If we have a listlike key, _check_indexing_error will raise
3626 # InvalidIndexError. Otherwise we fall through and re-raise
3627 # the TypeError.
3628 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.20.0/lib/python3.10/site-packages/pandas/core/indexes/base.py:3621, in Index.get_loc(self, key, method, tolerance)
3620 try:
-> 3621 return self._engine.get_loc(casted_key)
3622 except KeyError as err:
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pandas/_libs/index.pyx:136, in pandas._libs.index.IndexEngine.get_loc()
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pandas/_libs/index.pyx:163, in pandas._libs.index.IndexEngine.get_loc()
File pandas/_libs/hashtable_class_helper.pxi:5198, in pandas._libs.hashtable.PyObjectHashTable.get_item()
File pandas/_libs/hashtable_class_helper.pxi:5206, 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)
Input In [16], in <cell line: 1>()
----> 1 pd.DataFrame({attr: network.links_t[attr]["boiler"] for attr in ["p0", "p1"]})
Input In [16], 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.20.0/lib/python3.10/site-packages/pandas/core/frame.py:3505, in DataFrame.__getitem__(self, key)
3503 if self.columns.nlevels > 1:
3504 return self._getitem_multilevel(key)
-> 3505 indexer = self.columns.get_loc(key)
3506 if is_integer(indexer):
3507 indexer = [indexer]
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.20.0/lib/python3.10/site-packages/pandas/core/indexes/base.py:3623, in Index.get_loc(self, key, method, tolerance)
3621 return self._engine.get_loc(casted_key)
3622 except KeyError as err:
-> 3623 raise KeyError(key) from err
3624 except TypeError:
3625 # If we have a listlike key, _check_indexing_error will raise
3626 # InvalidIndexError. Otherwise we fall through and re-raise
3627 # the TypeError.
3628 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