Extraction-Condensing CHP¶

This example demonstrates how to model an extraction-condensing Combined Heat and Power (CHP) plant with flexible heat-power ratios. For an example of a CHP plant with a fixed heat-power ratio, see the backpressure CHP example.

In this example, a location has an electric, gas, and heat bus. The primary energy source is wind power, which can be converted to gas. The gas can be stored to convert into electricity and/or heat when needed, using either a boiler or a CHP unit.

In [2]:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import pypsa

import matplotlib.pyplot as plt import numpy as np import pandas as pd import pypsa

CHP parametrisation¶

The setup roughly follows Grohnheit (1993). We define the ratio between maximum heat and power output of the CHP unit as nom_r, the backpressure limit as c_m, and the marginal loss for each additional generation of heat as c_v.

In [3]:

nom_r = 1.0
c_m = 0.75
c_v = 0.15

nom_r = 1.0 c_m = 0.75 c_v = 0.15

The feasible operational space of the CHP unit is shown in the graph below:

In [4]:

fig, ax = plt.subplots()

t = 0.01
ph = np.arange(0, 1.0001, t)

ax.plot(ph, c_m * ph, color="k")
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="gray")
ax.plot(ph, 1 - c_v * ph, color="k")

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="orange", linestyle="--", linewidth=0.75)

ax.text(0.05, 0.41, "iso-fuel-lines", color="orange", rotation=-7)
ax.fill_between(ph, c_m * ph, 1 - c_v * ph, facecolor="#ddd")

fig, ax = plt.subplots() t = 0.01 ph = np.arange(0, 1.0001, t) ax.plot(ph, c_m * ph, color="k") 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="gray") ax.plot(ph, 1 - c_v * ph, color="k") 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="orange", linestyle="--", linewidth=0.75) ax.text(0.05, 0.41, "iso-fuel-lines", color="orange", rotation=-7) ax.fill_between(ph, c_m * ph, 1 - c_v * ph, facecolor="#ddd")

Out[4]:

<matplotlib.collections.FillBetweenPolyCollection at 0x7ecb115cb8c0>

Optimisation¶

First, we add the power and gas sectors with a wind generator, a power-to-gas unit, a gas storage unit, and the turbine of the CHP unit:

In [5]:

n = pypsa.Network()
n.set_snapshots(pd.date_range("2025-01-01 00:00", "2025-01-01 03:00", freq="h"))

n.add("Bus", "0 power", carrier="AC")
n.add("Bus", "0 gas", carrier="gas")

n.add("Carrier", ["wind", "gas"])

n.add(
    "Generator",
    "wind turbine",
    bus="0 power",
    carrier="wind",
    p_nom_extendable=True,
    p_max_pu=[0.0, 0.2, 0.7, 0.4],
    capital_cost=1000,
)

n.add("Load", "load", bus="0 power", p_set=5)

n.add(
    "Link",
    "power-to-gas",
    bus0="0 power",
    bus1="0 gas",
    efficiency=0.6,
    capital_cost=1000,
    p_nom_extendable=True,
)

n.add(
    "Link",
    "generator",
    bus0="0 gas",
    bus1="0 power",
    efficiency=0.468,
    capital_cost=400,
    p_nom_extendable=True,
)

n.add("Store", "gas depot", bus="0 gas", e_cyclic=True, e_nom=1000);

n = pypsa.Network() n.set_snapshots(pd.date_range("2025-01-01 00:00", "2025-01-01 03:00", freq="h")) n.add("Bus", "0 power", carrier="AC") n.add("Bus", "0 gas", carrier="gas") n.add("Carrier", ["wind", "gas"]) n.add( "Generator", "wind turbine", bus="0 power", carrier="wind", p_nom_extendable=True, p_max_pu=[0.0, 0.2, 0.7, 0.4], capital_cost=1000, ) n.add("Load", "load", bus="0 power", p_set=5) n.add( "Link", "power-to-gas", bus0="0 power", bus1="0 gas", efficiency=0.6, capital_cost=1000, p_nom_extendable=True, ) n.add( "Link", "generator", bus0="0 gas", bus1="0 power", efficiency=0.468, capital_cost=400, p_nom_extendable=True, ) n.add("Store", "gas depot", bus="0 gas", e_cyclic=True, e_nom=1000);

Next, we add the heat sector with the boiler of the CHP unit.

In [6]:

n.add("Bus", "0 heat", carrier="heat")

n.add("Carrier", "heat")

n.add("Load", "heat load", bus="0 heat", p_set=10)

n.add(
    "Link",
    "boiler",
    bus0="0 gas",
    bus1="0 heat",
    efficiency=0.9,
    capital_cost=300,
    p_nom_extendable=True,
)

n.add("Store", "water tank", bus="0 heat", e_cyclic=True, e_nom_extendable=True)

n.add("Bus", "0 heat", carrier="heat") n.add("Carrier", "heat") n.add("Load", "heat load", bus="0 heat", p_set=10) n.add( "Link", "boiler", bus0="0 gas", bus1="0 heat", efficiency=0.9, capital_cost=300, p_nom_extendable=True, ) n.add("Store", "water tank", bus="0 heat", e_cyclic=True, e_nom_extendable=True)

Next, we need to add some constraints that ensure that combination of turbine and boiler output of the CHP is operationally feasible. This is done by implementing the operational space shown in the graph above with a set of linear constraints.

In [7]:

# Guarantees ISO fuel lines, i.e. fuel consumption p_b0 + p_g0 = constant along p_g1 + c_v p_b1 = constant (b=boiler, g=generator)
n.links.at["boiler", "efficiency"] = n.links.at["generator", "efficiency"] / c_v
boiler_eff = float(n.links.at["boiler", "efficiency"])
generator_eff = float(n.links.at["generator", "efficiency"])

m = n.optimize.create_model()

p = m.variables["Link-p"]
p_nom = m.variables["Link-p_nom"]

# Guarantees heat output and electric output nominal powers are proportional
m.add_constraints(
    generator_eff * nom_r * p_nom.loc["generator"] - boiler_eff * p_nom.loc["boiler"]
    == 0,
    name="heat-power output proportionality",
)

# Guarantees c_m p_b1 <= p_g1
m.add_constraints(
    p.loc[:, "boiler"] * c_m * boiler_eff - p.loc[:, "generator"] * generator_eff <= 0,
    name="backpressure",
)

# Guarantees p_g1 +c_v p_b1 <= p_g1_nom
m.add_constraints(
    p.loc[:, "boiler"] + p.loc[:, "generator"] - p_nom.loc["generator"] <= 0,
    name="top_iso_fuel_line",
)

n.optimize.solve_model(log_to_console=False)

# Guarantees ISO fuel lines, i.e. fuel consumption p_b0 + p_g0 = constant along p_g1 + c_v p_b1 = constant (b=boiler, g=generator) n.links.at["boiler", "efficiency"] = n.links.at["generator", "efficiency"] / c_v boiler_eff = float(n.links.at["boiler", "efficiency"]) generator_eff = float(n.links.at["generator", "efficiency"]) m = n.optimize.create_model() p = m.variables["Link-p"] p_nom = m.variables["Link-p_nom"] # Guarantees heat output and electric output nominal powers are proportional m.add_constraints( generator_eff * nom_r * p_nom.loc["generator"] - boiler_eff * p_nom.loc["boiler"] == 0, name="heat-power output proportionality", ) # Guarantees c_m p_b1 <= p_g1 m.add_constraints( p.loc[:, "boiler"] * c_m * boiler_eff - p.loc[:, "generator"] * generator_eff <= 0, name="backpressure", ) # Guarantees p_g1 +c_v p_b1 <= p_g1_nom m.add_constraints( p.loc[:, "boiler"] + p.loc[:, "generator"] - p_nom.loc["generator"] <= 0, name="top_iso_fuel_line", ) n.optimize.solve_model(log_to_console=False)

WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['0 power'], dtype='object', name='name')

WARNING:pypsa.consistency:The following links have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['power-to-gas'], dtype='object', name='name')

/tmp/ipykernel_7700/1779618540.py:6: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  m = n.optimize.create_model()

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.08s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 37 primals, 82 duals
Objective: 1.62e+05
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-ext-e-lower, Store-ext-e-upper, Store-energy_balance, heat-power output proportionality, backpressure, top_iso_fuel_line were not assigned to the network.

Out[7]:

('ok', 'optimal')

In [8]:

n.objective

n.objective

Out[8]:

162253.99956169186

Results¶

Let's start by inspecting the optimised conversion capacities:

In [9]:

n.links.p_nom_opt

n.links.p_nom_opt

Out[9]:

name
power-to-gas    58.648915
generator       28.490028
boiler           4.273504
Name: p_nom_opt, dtype: float64

The CHP boiler is dimensioned by the heat demand met in three hours when there is no wind supply. The CHP generator is set by the heat demand.

In [10]:

display(4 * 10 / 3 / float(n.links.at["boiler", "efficiency"]))
display(28.490028 * 0.15)

display(4 * 10 / 3 / float(n.links.at["boiler", "efficiency"])) display(28.490028 * 0.15)

4.273504273504273

4.2735042

In [11]:

n.links_t.p0.round(2)

n.links_t.p0.round(2)

Out[11]:

name	power-to-gas	generator	boiler
snapshot
2025-01-01 00:00:00	5.00	21.37	4.27
2025-01-01 01:00:00	23.19	21.37	4.27
2025-01-01 02:00:00	58.65	-0.00	-0.00
2025-01-01 03:00:00	41.37	21.37	4.27

In [12]:

n.links_t.p1.round(2)

n.links_t.p1.round(2)

Out[12]:

name	power-to-gas	generator	boiler
snapshot
2025-01-01 00:00:00	-3.00	-10.0	-13.33
2025-01-01 01:00:00	-13.91	-10.0	-13.33
2025-01-01 02:00:00	-35.19	0.0	0.00
2025-01-01 03:00:00	-24.82	-10.0	-13.33

In [13]:

pd.DataFrame({attr: n.stores_t[attr]["gas depot"] for attr in ["p", "e"]}).round(2)

pd.DataFrame({attr: n.stores_t[attr]["gas depot"] for attr in ["p", "e"]}).round(2)

Out[13]:

	p	e
snapshot
2025-01-01 00:00:00	22.64	11.73
2025-01-01 01:00:00	11.73	-0.00
2025-01-01 02:00:00	-35.19	35.19
2025-01-01 03:00:00	0.82	34.37

In [14]:

pd.DataFrame({attr: n.stores_t[attr]["water tank"] for attr in ["p", "e"]}).round(2)

pd.DataFrame({attr: n.stores_t[attr]["water tank"] for attr in ["p", "e"]}).round(2)

Out[14]:

	p	e
snapshot
2025-01-01 00:00:00	-3.33	6.67
2025-01-01 01:00:00	-3.33	10.00
2025-01-01 02:00:00	10.00	-0.00
2025-01-01 03:00:00	-3.33	3.33

In [15]:

pd.DataFrame({attr: n.links_t[attr]["boiler"] for attr in ["p0", "p1"]}).round(2)

pd.DataFrame({attr: n.links_t[attr]["boiler"] for attr in ["p0", "p1"]}).round(2)

Out[15]:

	p0	p1
snapshot
2025-01-01 00:00:00	4.27	-13.33
2025-01-01 01:00:00	4.27	-13.33
2025-01-01 02:00:00	-0.00	0.00
2025-01-01 03:00:00	4.27	-13.33

Finally, let's calculate the overall efficiency of the CHP system.

In [16]:

eta_elec = n.links.at["generator", "efficiency"]

r = 1 / c_m

# P_h = r*P_e
(1 + r) / ((1 / eta_elec) * (1 + c_v * r))

eta_elec = n.links.at["generator", "efficiency"] r = 1 / c_m # P_h = r*P_e (1 + r) / ((1 / eta_elec) * (1 + c_v * r))

Out[16]:

np.float64(0.9099999999999999)