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
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_4538/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 Glpk solver
INFO:linopy.io: Writing time: 0.07s
INFO:linopy.solvers:GLPSOL--GLPK LP/MIP Solver 5.0
Parameter(s) specified in the command line:
 --lp /tmp/linopy-problem-i9hehm_7.lp --output /tmp/linopy-solve-olvsqdts.sol
Reading problem data from '/tmp/linopy-problem-i9hehm_7.lp'...
83 rows, 38 columns, 159 non-zeros
449 lines were read
GLPK Simplex Optimizer 5.0
83 rows, 38 columns, 159 non-zeros
Preprocessing...
52 rows, 37 columns, 127 non-zeros
Scaling...
 A: min|aij| =  2.000e-01  max|aij| =  3.120e+00  ratio =  1.560e+01
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part is 50
      0: obj =   0.000000000e+00 inf =   1.377e+02 (11)
     22: obj =   1.703538022e+05 inf =   0.000e+00 (0)
*    25: obj =   1.622539996e+05 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.1 Mb (104191 bytes)
Writing basic solution to '/tmp/linopy-solve-olvsqdts.sol'...

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

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/pypsa/optimization/optimize.py:357: FutureWarning: A value is trying to be set on a copy of a DataFrame or Series through chained assignment using an inplace method.
The behavior will change in pandas 3.0. This inplace method will never work because the intermediate object on which we are setting values always behaves as a copy.

For example, when doing 'df[col].method(value, inplace=True)', try using 'df.method({col: value}, inplace=True)' or df[col] = df[col].method(value) instead, to perform the operation inplace on the original object.


  n.df(c)[attr + "_opt"].update(df)
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.
[7]:
('ok', 'optimal')
[8]:
network.objective
[8]:
162253.9996

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.6489
generator    28.4900
boiler        4.2735
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.0000 21.3675 4.2735
2016-01-01 01:00:00 23.1854 21.3675 4.2735
2016-01-01 02:00:00 58.6489 0.0000 0.0000
2016-01-01 03:00:00 41.3708 21.3675 4.2735
[14]:
network.links_t.p1
[14]:
Link P2G generator boiler
snapshot
2016-01-01 00:00:00 -3.00000 -9.99999 -13.33332
2016-01-01 01:00:00 -13.91124 -9.99999 -13.33332
2016-01-01 02:00:00 -35.18934 -0.00000 -0.00000
2016-01-01 03:00:00 -24.82248 -9.99999 -13.33332
[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.64100 11.7298
2016-01-01 01:00:00 11.72980 0.0000
2016-01-01 02:00:00 -35.18930 35.1893
2016-01-01 03:00:00 0.81854 34.3708
[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.33333 6.66667
2016-01-01 01:00:00 -3.33333 10.00000
2016-01-01 02:00:00 10.00000 0.00000
2016-01-01 03:00:00 -3.33333 3.33333
[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.2735 -13.33332
2016-01-01 01:00:00 4.2735 -13.33332
2016-01-01 02:00:00 0.0000 -0.00000
2016-01-01 03:00:00 4.2735 -13.33332
[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.1893
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.927
Name: wind turbine, dtype: object
[20]:
network.links.p_nom_opt
[20]:
Link
P2G          58.6489
generator    28.4900
boiler        4.2735
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