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.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `kirchhoff` formulation
---------------------------------------------------------------------------
RuntimeError 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.19.3/lib/python3.10/site-packages/pypsa/components.py:733, in Network.lopf(self, snapshots, pyomo, solver_name, solver_options, solver_logfile, formulation, keep_files, extra_functionality, multi_investment_periods, **kwargs)
726 logger.warning(
727 "You have defined one or more shunt impedances. "
728 "Shunt impedances are ignored by the linear optimal "
729 "power flow (LOPF)."
730 )
732 if pyomo:
--> 733 return network_lopf(self, **args)
734 else:
735 return network_lopf_lowmem(self, **args)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.19.3/lib/python3.10/site-packages/pypsa/opf.py:2435, 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)
2429 logger.warning(
2430 "Encountered nonzero ramp limits for links. These are ignored when running the optimization with `pyomo=True`."
2431 )
2433 snapshots = _as_snapshots(network, snapshots)
-> 2435 network_lopf_build_model(
2436 network,
2437 snapshots,
2438 skip_pre=skip_pre,
2439 formulation=formulation,
2440 ptdf_tolerance=ptdf_tolerance,
2441 )
2443 if extra_functionality is not None:
2444 extra_functionality(network, snapshots)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.19.3/lib/python3.10/site-packages/pypsa/opf.py:2202, in network_lopf_build_model(network, snapshots, skip_pre, formulation, ptdf_tolerance)
2199 elif formulation in ["ptdf", "cycles"]:
2200 define_sub_network_balance_constraints(network, snapshots)
-> 2202 define_global_constraints(network, snapshots)
2204 define_linear_objective(network, snapshots)
2206 # tidy up auxilliary expressions
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.19.3/lib/python3.10/site-packages/pypsa/opf.py:1745, in define_global_constraints(network, snapshots)
1738 c.lhs.constant += sum(
1739 attribute * network.stores.at[store, "e_initial"]
1740 for store in stores
1741 )
1743 global_constraints[gc] = c
-> 1745 l_constraint(
1746 network.model,
1747 "global_constraints",
1748 global_constraints,
1749 list(network.global_constraints.index),
1750 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.19.3/lib/python3.10/site-packages/pypsa/opt.py:227, in l_constraint(model, name, constraints, *args)
222 else:
223 raise KeyError(
224 '`sense` must be one of "==","<=",">=","><"; got: {}'.format(sense)
225 )
--> 227 v._data[i] = _GeneralConstraintData(constr_expr, v)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.19.3/lib/python3.10/site-packages/pyomo/core/base/constraint.py:294, in _GeneralConstraintData.__init__(self, expr, component)
292 self._expr = None
293 if expr is not None:
--> 294 self.set_value(expr)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.19.3/lib/python3.10/site-packages/pyomo/core/base/constraint.py:543, in _GeneralConstraintData.set_value(self, expr)
528 raise ValueError(
529 "Invalid constraint expression. The constraint "
530 "expression resolved to a trivial Boolean (%s) "
(...)
534 % (expr, "Feasible" if expr else "Infeasible",
535 expr, self.name))
537 else:
538 msg = ("Constraint '%s' does not have a proper "
539 "value. Found '%s'\nExpecting a tuple or "
540 "equation. Examples:"
541 "\n sum(model.costs) == model.income"
542 "\n (0, model.price[item], 50)"
--> 543 % (self.name, str(expr)))
544 raise ValueError(msg)
545 #
546 # Normalize the incoming expressions, if we can
547 #
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.19.3/lib/python3.10/site-packages/pyomo/core/base/component.py:277, in _ComponentBase.name(self)
274 @property
275 def name(self):
276 """Get the fully qualifed component name."""
--> 277 return self.getname(fully_qualified=True)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.19.3/lib/python3.10/site-packages/pyomo/core/base/component.py:900, in ComponentData.getname(self, fully_qualified, name_buffer, relative_to)
898 return base + index_repr(idx)
899 #
--> 900 raise RuntimeError("Fatal error: cannot find the component data in "
901 "the owning component's _data dictionary.")
RuntimeError: Fatal error: cannot find the component data in the owning component's _data dictionary.
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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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.19.3/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