Problem 7.4#
Integrated Energy Grids
Problem 7.4
In this problem, we will model a node that contains an electricity, a gas and a heat bus. There is a demand of 50 MWh of electricity and 40 MWh of heat. Electricity can be produced using an Open Cycle Gas Turbine (OCGT) with an efficiency of 0.35 or using a Combined Heat and Power (CHP) unit. Heat can be produced with the CHP unit or with a gas boiler. The gas boiler has an efficiency of 0.9. The OCGT, gas boiler and CHP unit have a marginal cost due to the fuel cost of 20 EUR/MWh\(_{th}\)
(a) In this section, we assume that the CHP has a fixed efficiency of 0.3 when producing electricity and 0.3 when producing heat. Model the OCGT and gas boiler, using a link and the CHP unit using and multilink element in PyPSA. Add a gas store to the gas bus that represents an unlimited supply of gas. Optimize the system and calculate which technologies are supplying the electricity and heat demand.
(b) In this section, we assume that the CHP unit can be operated either in condensing mode or in back-pressure mode. In practice, this means the feasible space for operating the CHP unit is determined by the iso-fuel lines (with constant \(c_v\)=-0.15) and the back-pressure line (with constant \(c_m\)=0.75). Optimize the system and calculate which technologies are supplying the electricity and heat demand.
Note: This problem is based on the PyPSA examples Multi-Link usage: CHP with fixed heat-power ratio and Power to Gas with Heat Coupling
Note
If you have not yet set up Python on your computer, you can execute this tutorial in your browser via Google Colab. Click on the rocket in the top right corner and launch “Colab”. If that doesn’t work download the .ipynb file and import it in Google Colab.
Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.
!pip install numpy pypsa
We start by creating the network and adding the buses and loads.
import pypsa
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
network = pypsa.Network()
network.add("Bus", "electricity")
network.add("Load", "electricity load", bus="electricity", p_set=50)
network.add("Bus", "heat")
network.add("Load", "heat load", bus="heat", p_set=40)
network.add("Bus", "gas")
# We add a gas store with energy capacity and an initial filling level much higher than the required gas consumption,
# this way gas supply is unlimited
network.add("Store", "gas", e_initial=1e6, e_nom=1e6, bus="gas")
Index(['gas'], dtype='object')
network.buses
| v_nom | type | x | y | carrier | unit | location | v_mag_pu_set | v_mag_pu_min | v_mag_pu_max | control | generator | sub_network | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Bus | |||||||||||||
| electricity | 1.0 | 0.0 | 0.0 | AC | 1.0 | 0.0 | inf | PQ | |||||
| heat | 1.0 | 0.0 | 0.0 | AC | 1.0 | 0.0 | inf | PQ | |||||
| gas | 1.0 | 0.0 | 0.0 | AC | 1.0 | 0.0 | inf | PQ |
We add the links representing the OCGT and the CHP. For the later, we can use a multilink, that is a link that outputs to two buses (electricity and heating).
network.add(
"Link",
"OCGT",
bus0="gas",
bus1="electricity",
p_nom=1000,
marginal_cost=20,
efficiency=0.35,
)
network.add(
"Link",
"gas boiler",
bus0="gas",
bus1="heat",
marginal_cost=20,
p_nom=1000,
efficiency=0.9,
)
network.add(
"Link",
"CHP",
bus0="gas",
bus1="electricity",
bus2="heat",
p_nom=1000,
marginal_cost=20,
efficiency=0.3,
efficiency2=0.3,
)
network.links
| bus0 | bus1 | type | carrier | efficiency | active | build_year | lifetime | p_nom | p_nom_mod | ... | min_down_time | up_time_before | down_time_before | ramp_limit_up | ramp_limit_down | ramp_limit_start_up | ramp_limit_shut_down | p_nom_opt | bus2 | efficiency2 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Link | |||||||||||||||||||||
| OCGT | gas | electricity | 0.35 | True | 0 | inf | 1000.0 | 0.0 | ... | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | NaN | |||
| gas boiler | gas | heat | 0.90 | True | 0 | inf | 1000.0 | 0.0 | ... | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | NaN | |||
| CHP | gas | electricity | 0.30 | True | 0 | inf | 1000.0 | 0.0 | ... | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | heat | 0.3 |
3 rows × 36 columns
network.optimize();
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['OCGT', 'gas boiler', 'CHP'], dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['gas'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'heat', 'gas'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io: Writing time: 0.02s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 5 primals, 12 duals
Objective: 3.24e+03
Solver model: available
Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Link-fix-p-lower, Link-fix-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-energy_balance were not assigned to the network.
Running HiGHS 1.10.0 (git hash: fd86653): Copyright (c) 2025 HiGHS under MIT licence terms
LP linopy-problem-xq1kvhqg has 12 rows; 5 cols; 18 nonzeros
Coefficient ranges:
Matrix [3e-01, 1e+00]
Cost [2e+01, 2e+01]
Bound [0e+00, 0e+00]
RHS [4e+01, 1e+06]
Presolving model
0 rows, 0 cols, 0 nonzeros 0s
0 rows, 0 cols, 0 nonzeros 0s
Presolve : Reductions: rows 0(-12); columns 0(-5); elements 0(-18) - Reduced to empty
Solving the original LP from the solution after postsolve
Model name : linopy-problem-xq1kvhqg
Model status : Optimal
Objective value : 3.2380952381e+03
Relative P-D gap : 0.0000000000e+00
HiGHS run time : 0.00
Writing the solution to /tmp/linopy-solve-p4_9gjgc.sol
network.loads_t.p
| Load | electricity load | heat load |
|---|---|---|
| snapshot | ||
| now | 50.0 | 40.0 |
network.links_t.p0
| Link | OCGT | gas boiler | CHP |
|---|---|---|---|
| snapshot | |||
| now | 28.571429 | -0.0 | 133.333333 |
network.links_t.p1
| Link | OCGT | gas boiler | CHP |
|---|---|---|---|
| snapshot | |||
| now | -10.0 | 0.0 | -40.0 |
network.links_t.p2
| Link | OCGT | gas boiler | CHP |
|---|---|---|---|
| snapshot | |||
| now | 0.0 | 0.0 | -40.0 |
(b) In this section, we assume that the CHP unit can be operated either in condensing mode or in back-pressure mode. In practice, this means the feasible space for operating the CHP unit is determined by the iso-fuel lines (with constant \(c_v\)=-0.15) and the back-pressure line (with constant \(c_m\)=0.75). Optimize the system and calculate which technologies are supplying the electricity and heat demand.
Combined-Heat-and-Power (CHP) parameterisation#
The parameters describing the operation of the CPH unit follow https://doi.org/10.1016/0301-4215(93)90282-K
# backpressure limit
c_m = 0.75
# marginal loss of electricity generation for each additional generation of heat
c_v = 0.15
We can plot a figure that shows the operating space of the CHP units
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("Heat generation")
ax.set_ylabel("Electricity generation")
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()
We can build the network and add the buses, loads and links representing the OCGT and CHP unit. We use two idependent link to represent the generation of electricity and heat in the CHP unit and later we will impose constraints limiting the feasible space (as shown in the figure above).
network = pypsa.Network()
network.add("Bus", "electricity")
network.add("Load", "electricity load", bus="electricity", p_set=50)
network.add("Bus", "heat")
network.add("Load", "heat load", bus="heat", p_set=40)
network.add("Bus", "gas")
# We add a gas store with energy capacity and an initial filling level much higher than the required gas consumption,
# this way gas supply is unlimited
network.add("Store", "gas", e_initial=1e6, e_nom=1e6, bus="gas")
network.add(
"Link",
"OCGT",
bus0="gas",
bus1="electricity",
p_nom=1000,
marginal_cost=20,
efficiency=0.35,
)
network.add(
"Link",
"CHP generator",
bus0="gas",
bus1="electricity",
efficiency=0.3,
p_nom=1000,
marginal_cost=20,
)
network.add(
"Link",
"CHP boiler",
bus0="gas",
bus1="heat",
p_nom=1000,
)
network.add(
"Link",
"gas boiler",
bus0="gas",
bus1="heat",
marginal_cost=20,
p_nom=1000,
efficiency=0.9,
)
Index(['gas boiler'], dtype='object')
network.loads
| bus | carrier | type | p_set | q_set | sign | active | |
|---|---|---|---|---|---|---|---|
| Load | |||||||
| electricity load | electricity | 50.0 | 0.0 | -1.0 | True | ||
| heat load | heat | 40.0 | 0.0 | -1.0 | True |
We add the CHP constraints to limite the feasible operation space
# Guarantees isofuel lines, i.e. increase in heat generation is proportional to decrease in electricity generation
network.links.at["CHP boiler", "efficiency"] = (network.links.at["CHP generator", "efficiency"] / c_v)
model = network.optimize.create_model()
link_p = model.variables["Link-p"]
# Guarantees back-pressure line
model.add_constraints(
c_m * network.links.at["CHP boiler", "efficiency"] * link_p.sel(Link="CHP boiler")
- network.links.at["CHP generator", "efficiency"] * link_p.sel(Link="CHP generator")
<= 0,
name="backpressure",)
# Guarantees top iso fuel line
model.add_constraints(
link_p.sel(Link="CHP boiler") + link_p.sel(Link="CHP generator")
- network.links.at["CHP generator", "p_nom"]
<= 0,
name="top_iso_fuel_line",
)
network.optimize.solve_model()
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['OCGT', 'CHP generator', 'CHP boiler', 'gas boiler'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['gas'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'heat', 'gas'], dtype='object', name='Bus')
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[14], line 10
6 link_p = model.variables["Link-p"]
8 # Guarantees back-pressure line
9 model.add_constraints(
---> 10 c_m * network.links.at["CHP boiler", "efficiency"] * link_p.sel(Link="CHP boiler")
11 - network.links.at["CHP generator", "efficiency"] * link_p.sel(Link="CHP generator")
12 <= 0,
13 name="backpressure",)
16 # Guarantees top iso fuel line
17 model.add_constraints(
18 link_p.sel(Link="CHP boiler") + link_p.sel(Link="CHP generator")
19 - network.links.at["CHP generator", "p_nom"]
20 <= 0,
21 name="top_iso_fuel_line",
22 )
File /opt/hostedtoolcache/Python/3.11.12/x64/lib/python3.11/site-packages/linopy/variables.py:414, in Variable.__rmul__(self, other)
410 """
411 Right-multiply variables with a coefficient.
412 """
413 try:
--> 414 return self.to_linexpr(other)
415 except TypeError:
416 return NotImplemented
File /opt/hostedtoolcache/Python/3.11.12/x64/lib/python3.11/site-packages/linopy/variables.py:316, in Variable.to_linexpr(self, coefficient)
293 def to_linexpr(
294 self,
295 coefficient: int
(...) 300 | DataArray = 1,
301 ) -> expressions.LinearExpression:
302 """
303 Create a linear expression from the variables.
304
(...) 314 Linear expression with the variables and coefficients.
315 """
--> 316 coefficient = as_dataarray(coefficient, coords=self.coords, dims=self.dims)
317 ds = Dataset({"coeffs": coefficient, "vars": self.labels}).expand_dims(
318 TERM_DIM, -1
319 )
320 return expressions.LinearExpression(ds, self.model)
File /opt/hostedtoolcache/Python/3.11.12/x64/lib/python3.11/site-packages/linopy/common.py:258, in as_dataarray(arr, coords, dims, **kwargs)
256 arr = numpy_to_dataarray(arr, coords=coords, dims=dims, **kwargs)
257 elif isinstance(arr, (np.number, int, float, str, bool, list)):
--> 258 arr = DataArray(arr, coords=coords, dims=dims, **kwargs)
260 elif not isinstance(arr, DataArray):
261 supported_types = [
262 np.number,
263 str,
(...) 269 DataArray,
270 ]
File /opt/hostedtoolcache/Python/3.11.12/x64/lib/python3.11/site-packages/xarray/core/dataarray.py:479, in DataArray.__init__(self, data, coords, dims, name, attrs, indexes, fastpath)
477 data = _check_data_shape(data, coords, dims)
478 data = as_compatible_data(data)
--> 479 coords, dims = _infer_coords_and_dims(data.shape, coords, dims)
480 variable = Variable(dims, data, attrs, fastpath=True)
482 if not isinstance(coords, Coordinates):
File /opt/hostedtoolcache/Python/3.11.12/x64/lib/python3.11/site-packages/xarray/core/dataarray.py:184, in _infer_coords_and_dims(shape, coords, dims)
182 dims_tuple = tuple(dims)
183 if len(dims_tuple) != len(shape):
--> 184 raise ValueError(
185 "different number of dimensions on data "
186 f"and dims: {len(shape)} vs {len(dims_tuple)}"
187 )
188 for d in dims_tuple:
189 if not hashable(d):
ValueError: different number of dimensions on data and dims: 0 vs 1
network.links_t.p1
| Link | OCGT | CHP generator | CHP boiler | gas boiler |
|---|---|---|---|---|
| snapshot | ||||
| now | -20.0 | -30.0 | -40.0 | 0.0 |