LOPF with coupling to heating sector

Note

You can download this example as a Jupyter notebook or start it in interactive mode.

LOPF with coupling to heating sector#

In this example three locations are optimised, each with an electric bus and a heating bus and corresponding loads. At each location the electric and heating buses are connected with heat pumps; heat can also be supplied to the heat bus with a boiler. The electric buses are connected with transmission lines and there are electrical generators at two of the nodes.

[1]:

%pip install seaborn

Requirement already satisfied: seaborn in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (0.13.2)
Requirement already satisfied: numpy!=1.24.0,>=1.20 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from seaborn) (2.2.5)
Requirement already satisfied: pandas>=1.2 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from seaborn) (2.2.3)
Requirement already satisfied: matplotlib!=3.6.1,>=3.4 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from seaborn) (3.10.1)
Requirement already satisfied: contourpy>=1.0.1 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (1.3.2)
Requirement already satisfied: cycler>=0.10 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (4.57.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (1.4.8)
Requirement already satisfied: packaging>=20.0 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (25.0)
Requirement already satisfied: pillow>=8 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (11.2.1)
Requirement already satisfied: pyparsing>=2.3.1 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (3.2.3)
Requirement already satisfied: python-dateutil>=2.7 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (2.9.0.post0)
Requirement already satisfied: pytz>=2020.1 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from pandas>=1.2->seaborn) (2025.2)
Requirement already satisfied: tzdata>=2022.7 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from pandas>=1.2->seaborn) (2025.2)
Requirement already satisfied: six>=1.5 in /home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages (from python-dateutil>=2.7->matplotlib!=3.6.1,>=3.4->seaborn) (1.17.0)
Note: you may need to restart the kernel to use updated packages.

[2]:

import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

import pypsa

sns.set(rc={"figure.figsize": (9, 5)})

[3]:

network = pypsa.Network()

Add three buses of AC and heat carrier each

[4]:

for i in range(3):
    network.add("Bus", f"electric bus {i}", v_nom=20.0)
    network.add("Bus", f"heat bus {i}", carrier="heat")
network.buses

[4]:
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
electric bus 0 20.0 0.0 0.0 AC 1.0 0.0 inf PQ
heat bus 0 1.0 0.0 0.0 heat 1.0 0.0 inf PQ
electric bus 1 20.0 0.0 0.0 AC 1.0 0.0 inf PQ
heat bus 1 1.0 0.0 0.0 heat 1.0 0.0 inf PQ
electric bus 2 20.0 0.0 0.0 AC 1.0 0.0 inf PQ
heat bus 2 1.0 0.0 0.0 heat 1.0 0.0 inf PQ 
[5]:

network.buses["carrier"].value_counts()

[5]:

carrier
AC      3
heat    3
Name: count, dtype: int64

Add three lines in a ring

[6]:

for i in range(3):
    network.add(
        "Line",
        f"line {i}",
        bus0=f"electric bus {i}",
        bus1=f"electric bus {(i + 1) % 3}",
        x=0.1,
        s_nom=1000,
    )
network.lines

[6]:
bus0 bus1 type x r g b s_nom s_nom_mod s_nom_extendable ... v_ang_min v_ang_max sub_network x_pu r_pu g_pu b_pu x_pu_eff r_pu_eff s_nom_opt
Line
line 0 electric bus 0 electric bus 1 0.1 0.0 0.0 0.0 1000.0 0.0 False ... -inf inf 0.0 0.0 0.0 0.0 0.0 0.0 0.0
line 1 electric bus 1 electric bus 2 0.1 0.0 0.0 0.0 1000.0 0.0 False ... -inf inf 0.0 0.0 0.0 0.0 0.0 0.0 0.0
line 2 electric bus 2 electric bus 0 0.1 0.0 0.0 0.0 1000.0 0.0 False ... -inf inf 0.0 0.0 0.0 0.0 0.0 0.0 0.0

3 rows × 31 columns

Connect the electric to the heat buses with heat pumps with COP 3

[7]:

for i in range(3):
    network.add(
        "Link",
        f"heat pump {i}",
        bus0=f"electric bus {i}",
        bus1=f"heat bus {i}",
        p_nom=100,
        efficiency=3.0,
    )
network.links

[7]:
bus0 bus1 type carrier efficiency active build_year lifetime p_nom p_nom_mod ... shut_down_cost min_up_time 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
Link
heat pump 0 electric bus 0 heat bus 0 3.0 True 0 inf 100.0 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
heat pump 1 electric bus 1 heat bus 1 3.0 True 0 inf 100.0 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0
heat pump 2 electric bus 2 heat bus 2 3.0 True 0 inf 100.0 0.0 ... 0.0 0 0 1 0 NaN NaN 1.0 1.0 0.0

3 rows × 34 columns

Add carriers

[8]:

network.add("Carrier", "gas", co2_emissions=0.27)
network.add("Carrier", "biomass", co2_emissions=0.0)
network.carriers

[8]:
co2_emissions color nice_name max_growth max_relative_growth
Carrier
gas 0.27 inf 0.0
biomass 0.00 inf 0.0

Add a gas generator at bus 0, a biomass generator at bus 1 and a boiler at all heat buses

[9]:

network.add(
    "Generator",
    "gas generator",
    bus="electric bus 0",
    p_nom=100,
    marginal_cost=50,
    carrier="gas",
    efficiency=0.3,
)

network.add(
    "Generator",
    "biomass generator",
    bus="electric bus 1",
    p_nom=100,
    marginal_cost=100,
    efficiency=0.3,
    carrier="biomass",
)

for i in range(3):
    network.add(
        "Generator",
        f"boiler {i}",
        bus=f"heat bus {i}",
        p_nom=1000,
        efficiency=0.9,
        marginal_cost=20.0,
        carrier="gas",
    )

network.generators

[9]:
bus control type p_nom p_nom_mod p_nom_extendable p_nom_min p_nom_max p_min_pu p_max_pu ... min_up_time min_down_time up_time_before down_time_before ramp_limit_up ramp_limit_down ramp_limit_start_up ramp_limit_shut_down weight p_nom_opt
Generator
gas generator electric bus 0 PQ 100.0 0.0 False 0.0 inf 0.0 1.0 ... 0 0 1 0 NaN NaN 1.0 1.0 1.0 0.0
biomass generator electric bus 1 PQ 100.0 0.0 False 0.0 inf 0.0 1.0 ... 0 0 1 0 NaN NaN 1.0 1.0 1.0 0.0
boiler 0 heat bus 0 PQ 1000.0 0.0 False 0.0 inf 0.0 1.0 ... 0 0 1 0 NaN NaN 1.0 1.0 1.0 0.0
boiler 1 heat bus 1 PQ 1000.0 0.0 False 0.0 inf 0.0 1.0 ... 0 0 1 0 NaN NaN 1.0 1.0 1.0 0.0
boiler 2 heat bus 2 PQ 1000.0 0.0 False 0.0 inf 0.0 1.0 ... 0 0 1 0 NaN NaN 1.0 1.0 1.0 0.0

5 rows × 37 columns

Add electric loads and heat loads.

[10]:

for i in range(3):
    network.add(
        "Load",
        f"electric load {i}",
        bus=f"electric bus {i}",
        p_set=i * 10,
    )

for i in range(3):
    network.add(
        "Load",
        f"heat load {i}",
        bus=f"heat bus {i}",
        p_set=(3 - i) * 10,
    )

network.loads

[10]:
bus carrier type p_set q_set sign active
Load
electric load 0 electric bus 0 0.0 0.0 -1.0 True
electric load 1 electric bus 1 10.0 0.0 -1.0 True
electric load 2 electric bus 2 20.0 0.0 -1.0 True
heat load 0 heat bus 0 30.0 0.0 -1.0 True
heat load 1 heat bus 1 20.0 0.0 -1.0 True
heat load 2 heat bus 2 10.0 0.0 -1.0 True

We define a function for the LOPF

[11]:

def run_lopf():
    network.optimize()
    df = pd.concat(
        [
            network.generators_t.p.loc["now"],
            network.links_t.p0.loc["now"],
            network.loads_t.p.loc["now"],
        ],
        keys=["Generators", "Links", "Line"],
        names=["Component", "index"],
    ).reset_index(name="Production")

    sns.barplot(data=df, x="index", y="Production", hue="Component")
    plt.title(f"Objective: {network.objective}")
    plt.xticks(rotation=90)
    plt.tight_layout()

[12]:

run_lopf()

WARNING:pypsa.consistency:The following lines have carriers which are not defined:
Index(['line 0', 'line 1', 'line 2'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow:
Index(['line 0', 'line 1', 'line 2'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['heat pump 0', 'heat pump 1', 'heat pump 2'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electric bus 0', 'heat bus 0', 'electric bus 1', 'heat bus 1',
       'electric bus 2', 'heat bus 2'],
      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: 11 primals, 29 duals
Objective: 2.50e+03
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Link-fix-p-lower, Link-fix-p-upper, Kirchhoff-Voltage-Law 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-aczwx38h has 29 rows; 11 cols; 42 nonzeros
Coefficient ranges:
  Matrix [1e+00, 2e+01]
  Cost   [2e+01, 1e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e+01, 1e+03]
Presolving model
4 rows, 8 cols, 14 nonzeros  0s
Dependent equations search running on 4 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
4 rows, 8 cols, 14 nonzeros  0s
Presolve : Reductions: rows 4(-25); columns 8(-3); elements 14(-28)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Ph1: 0(0) 0s
          4     2.5000000000e+03 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-aczwx38h
Model status        : Optimal
Simplex   iterations: 4
Objective value     :  2.5000000000e+03
Relative P-D gap    :  1.8189894035e-16
HiGHS run time      :          0.00
Writing the solution to /tmp/linopy-solve-aj1rolpv.sol
../_images/examples_lopf-with-heating_19_2.png

Now, rerun with marginal costs for the heat pump operation.

[13]:

network.links.marginal_cost = 10
run_lopf()

WARNING:pypsa.consistency:The following lines have carriers which are not defined:
Index(['line 0', 'line 1', 'line 2'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow:
Index(['line 0', 'line 1', 'line 2'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['heat pump 0', 'heat pump 1', 'heat pump 2'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following sub_networks have carriers which are not defined:
Index(['0'], dtype='object', name='SubNetwork')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electric bus 0', 'heat bus 0', 'electric bus 1', 'heat bus 1',
       'electric bus 2', 'heat bus 2'],
      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: 11 primals, 29 duals
Objective: 2.70e+03
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Link-fix-p-lower, Link-fix-p-upper, Kirchhoff-Voltage-Law 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-oxquyte9 has 29 rows; 11 cols; 42 nonzeros
Coefficient ranges:
  Matrix [1e+00, 2e+01]
  Cost   [1e+01, 1e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e+01, 1e+03]
Presolving model
4 rows, 8 cols, 14 nonzeros  0s
Dependent equations search running on 4 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
4 rows, 7 cols, 13 nonzeros  0s
Presolve : Reductions: rows 4(-25); columns 7(-4); elements 13(-29)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Ph1: 0(0) 0s
          4     2.7000000000e+03 Pr: 0(0); Du: 0(3.55271e-15) 0s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-oxquyte9
Model status        : Optimal
Simplex   iterations: 4
Objective value     :  2.7000000000e+03
Relative P-D gap    :  1.6842494477e-16
HiGHS run time      :          0.00
Writing the solution to /tmp/linopy-solve-fi6g857w.sol
../_images/examples_lopf-with-heating_21_2.png

Finally, rerun with no CO2 emissions.

[14]:

network.add("GlobalConstraint", "co2_limit", sense="<=", constant=0.0)

run_lopf()

WARNING:pypsa.consistency:The following lines have carriers which are not defined:
Index(['line 0', 'line 1', 'line 2'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow:
Index(['line 0', 'line 1', 'line 2'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['heat pump 0', 'heat pump 1', 'heat pump 2'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following sub_networks have carriers which are not defined:
Index(['0'], dtype='object', name='SubNetwork')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electric bus 0', 'heat bus 0', 'electric bus 1', 'heat bus 1',
       'electric bus 2', 'heat bus 2'],
      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: 11 primals, 30 duals
Objective: 5.20e+03
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Link-fix-p-lower, Link-fix-p-upper, Kirchhoff-Voltage-Law 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-us6cb32k has 30 rows; 11 cols; 46 nonzeros
Coefficient ranges:
  Matrix [3e-01, 2e+01]
  Cost   [1e+01, 1e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e+01, 1e+03]
Presolving model
4 rows, 4 cols, 10 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-30); columns 0(-11); elements 0(-46) - Reduced to empty
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-us6cb32k
Model status        : Optimal
Objective value     :  5.2000000000e+03
Relative P-D gap    :  1.7490282726e-16
HiGHS run time      :          0.00
Writing the solution to /tmp/linopy-solve-056g0xfa.sol
../_images/examples_lopf-with-heating_23_2.png