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 matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pypsa
%matplotlib inline
[2]:
import logging
logging.basicConfig(level="INFO")
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
[3]:
# 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
[4]:
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#
[5]:
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)
/tmp/ipykernel_3198/1732597205.py:2: FutureWarning:
'H' is deprecated and will be removed in a future version, please use 'h' instead.
[5]:
Index(['gas depot'], dtype='object')
Add heat sector
[6]:
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)
[6]:
Index(['water tank'], dtype='object')
Add CHP constraints
[7]:
# 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
)
boiler_eff = float(network.links.at["boiler", "efficiency"])
generator_eff = float(network.links.at["generator", "efficiency"])
model = network.optimize.create_model()
link_p = model.variables["Link-p"]
link_p_nom = model.variables["Link-p_nom"]
# Guarantees heat output and electric output nominal powers are proportional
model.add_constraints(
generator_eff * nom_r * link_p_nom.loc["generator"]
- boiler_eff * link_p_nom.loc["boiler"]
== 0,
name="heat-power output proportionality",
)
# Guarantees c_m p_b1 \leq p_g1
model.add_constraints(
link_p.loc[:, "boiler"] * c_m * boiler_eff
- link_p.loc[:, "generator"] * generator_eff
<= 0,
name="backpressure",
)
# Guarantees p_g1 +c_v p_b1 \leq p_g1_nom
model.add_constraints(
link_p.loc[:, "boiler"] + link_p.loc[:, "generator"] - link_p_nom.loc["generator"]
<= 0,
name="top_iso_fuel_line",
)
network.optimize.solve_model()
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['P2G'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['0'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io: Writing time: 0.05s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 38 primals, 83 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-ext-e-lower, Store-ext-e-upper, Store-energy_balance, backpressure, top_iso_fuel_line 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-1wzk9x6q has 83 rows; 38 cols; 159 nonzeros
Coefficient ranges:
Matrix [2e-01, 3e+00]
Cost [3e+02, 1e+03]
Bound [0e+00, 0e+00]
RHS [5e+00, 1e+01]
Presolving model
38 rows, 29 cols, 98 nonzeros 0s
31 rows, 22 cols, 90 nonzeros 0s
Dependent equations search running on 7 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
30 rows, 21 cols, 88 nonzeros 0s
Presolve : Reductions: rows 30(-53); columns 21(-17); elements 88(-71)
Solving the presolved LP
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
0 -8.8568328443e-05 Pr: 15(126.838) 0s
24 1.6225399956e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model name : linopy-problem-1wzk9x6q
Model status : Optimal
Simplex iterations: 24
Objective value : 1.6225399956e+05
Relative P-D gap : 3.5874407454e-16
HiGHS run time : 0.00
Writing the solution to /tmp/linopy-solve-f1qvtlyu.sol
[7]:
('ok', 'optimal')
[8]:
network.objective
[8]:
162253.99956169186
Inspection#
[9]:
network.loads_t.p
[9]:
| Load | load | heat load |
|---|---|---|
| snapshot | ||
| 2016-01-01 00:00:00 | 5.0 | 10.0 |
| 2016-01-01 01:00:00 | 5.0 | 10.0 |
| 2016-01-01 02:00:00 | 5.0 | 10.0 |
| 2016-01-01 03:00:00 | 5.0 | 10.0 |
[10]:
network.links.p_nom_opt
[10]:
Link
P2G 58.648915
generator 28.490028
boiler 4.273504
Name: p_nom_opt, dtype: float64
[11]:
# CHP is dimensioned by the heat demand met in three hours when no wind
4 * 10.0 / 3 / network.links.at["boiler", "efficiency"]
[11]:
np.float64(4.273504273504273)
[12]:
# elec is set by the heat demand
28.490028 * 0.15
[12]:
4.2735042
[13]:
network.links_t.p0
[13]:
| Link | P2G | generator | boiler |
|---|---|---|---|
| snapshot | |||
| 2016-01-01 00:00:00 | 5.000000 | 21.367521 | 4.273504 |
| 2016-01-01 01:00:00 | 23.185404 | 21.367521 | 4.273504 |
| 2016-01-01 02:00:00 | 58.648915 | -0.000000 | -0.000000 |
| 2016-01-01 03:00:00 | 41.370809 | 21.367521 | 4.273504 |
[14]:
network.links_t.p1
[14]:
| Link | P2G | generator | boiler |
|---|---|---|---|
| snapshot | |||
| 2016-01-01 00:00:00 | -3.000000 | -10.0 | -13.333333 |
| 2016-01-01 01:00:00 | -13.911243 | -10.0 | -13.333333 |
| 2016-01-01 02:00:00 | -35.189349 | 0.0 | 0.000000 |
| 2016-01-01 03:00:00 | -24.822485 | -10.0 | -13.333333 |
[15]:
pd.DataFrame({attr: network.stores_t[attr]["gas depot"] for attr in ["p", "e"]})
[15]:
| p | e | |
|---|---|---|
| snapshot | ||
| 2016-01-01 00:00:00 | 22.641026 | 11.729783 |
| 2016-01-01 01:00:00 | 11.729783 | -0.000000 |
| 2016-01-01 02:00:00 | -35.189349 | 35.189349 |
| 2016-01-01 03:00:00 | 0.818540 | 34.370809 |
[16]:
pd.DataFrame({attr: network.stores_t[attr]["water tank"] for attr in ["p", "e"]})
[16]:
| p | e | |
|---|---|---|
| snapshot | ||
| 2016-01-01 00:00:00 | -3.333333 | 6.666667 |
| 2016-01-01 01:00:00 | -3.333333 | 10.000000 |
| 2016-01-01 02:00:00 | 10.000000 | -0.000000 |
| 2016-01-01 03:00:00 | -3.333333 | 3.333333 |
[17]:
pd.DataFrame({attr: network.links_t[attr]["boiler"] for attr in ["p0", "p1"]})
[17]:
| p0 | p1 | |
|---|---|---|
| snapshot | ||
| 2016-01-01 00:00:00 | 4.273504 | -13.333333 |
| 2016-01-01 01:00:00 | 4.273504 | -13.333333 |
| 2016-01-01 02:00:00 | -0.000000 | 0.000000 |
| 2016-01-01 03:00:00 | 4.273504 | -13.333333 |
[18]:
network.stores.loc["gas depot"]
[18]:
bus 0 gas
type
carrier gas
e_nom 0.0
e_nom_mod 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
marginal_cost_quadratic 0.0
marginal_cost_storage 0.0
capital_cost 0.0
standing_loss 0.0
active True
build_year 0
lifetime inf
e_nom_opt 35.189349
Name: gas depot, dtype: object
[19]:
network.generators.loc["wind turbine"]
[19]:
bus 0
control PQ
type
p_nom 0.0
p_nom_mod 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
e_sum_min -inf
e_sum_max inf
q_set 0.0
sign 1.0
carrier wind
marginal_cost 0.0
marginal_cost_quadratic 0.0
active True
build_year 0
lifetime inf
capital_cost 1000.0
efficiency 1.0
committable False
start_up_cost 0.0
shut_down_cost 0.0
stand_by_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
weight 1.0
p_nom_opt 90.927022
Name: wind turbine, dtype: object
[20]:
network.links.p_nom_opt
[20]:
Link
P2G 58.648915
generator 28.490028
boiler 4.273504
Name: p_nom_opt, dtype: float64
Calculate the overall efficiency of the CHP
[21]:
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))
[21]:
np.float64(0.9099999999999999)
[ ]: