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
ERROR 1: PROJ: proj_create_from_database: Open of /home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/share/proj failed
[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()
../_images/examples_power-to-gas-boiler-chp_6_0.png

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_4550/1732597205.py:2: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
  network.set_snapshots(pd.date_range("2016-01-01 00:00", "2016-01-01 03:00", freq="H"))

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)

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
)

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(
    network.links.at["generator", "efficiency"] * nom_r * link_p_nom["generator"]
    - network.links.at["boiler", "efficiency"] * link_p_nom["boiler"]
    == 0,
    name="heat-power output proportionality",
)

# Guarantees c_m p_b1  \leq p_g1
model.add_constraints(
    c_m * network.links.at["boiler", "efficiency"] * link_p.sel(Link="boiler")
    - network.links.at["generator", "efficiency"] * link_p.sel(Link="generator")
    <= 0,
    name="backpressure",
)

# Guarantees p_g1 +c_v p_b1 \leq p_g1_nom
model.add_constraints(
    link_p.sel(Link="boiler")
    + link_p.sel(Link="generator")
    - link_p_nom.sel({"Link-ext": "generator"})
    <= 0,
    name="top_iso_fuel_line",
)

network.optimize.solve_model()
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io: Writing time: 0.07s
INFO:linopy.solvers:Log file at /tmp/highs.log.
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.7.0 (git hash: 27ccfaa): Copyright (c) 2024 HiGHS under MIT licence terms
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
35 rows, 26 cols, 98 nonzeros  0s
Presolve : Reductions: rows 35(-48); columns 26(-12); elements 98(-61)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0    -7.0546311750e-05 Pr: 16(126.838); Du: 0(2.16492e-12) 0s
         29     1.6225399956e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 29
Objective value     :  1.6225399956e+05
HiGHS run time      :          0.00
[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]:
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]:
attribute
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
capital_cost                     0.0
standing_loss                    0.0
build_year                         0
lifetime                         inf
e_nom_opt                  35.189349
Name: gas depot, dtype: object
[19]:
network.generators.loc["wind turbine"]
[19]:
attribute
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
q_set                            0.0
sign                             1.0
carrier                         wind
marginal_cost                    0.0
marginal_cost_quadratic          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
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]:
0.9099999999999999