Example of linear optimal power flow 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.

Available as a Jupyter notebook at http://www.pypsa.org/examples/lopf-with-heating.ipynb.

In [1]:
# make the code as Python 3 compatible as possible 
from __future__ import print_function, division

import pypsa

import numpy as np

import pandas as pd
In [2]:
network = pypsa.Network()
In [3]:
#add three buses of AC and heat carrier each
for i in range(3):
    network.add("Bus","electric bus {}".format(i),v_nom=20.)
    network.add("Bus","heat bus {}".format(i),carrier="heat")
In [4]:
print(network.buses)
attribute       v_nom type    x    y carrier  v_mag_pu_set  v_mag_pu_min  \
electric bus 0   20.0       0.0  0.0      AC           1.0           0.0   
heat bus 0        1.0       0.0  0.0    heat           1.0           0.0   
electric bus 1   20.0       0.0  0.0      AC           1.0           0.0   
heat bus 1        1.0       0.0  0.0    heat           1.0           0.0   
electric bus 2   20.0       0.0  0.0      AC           1.0           0.0   
heat bus 2        1.0       0.0  0.0    heat           1.0           0.0   

attribute       v_mag_pu_max control sub_network  
electric bus 0           inf      PQ              
heat bus 0               inf      PQ              
electric bus 1           inf      PQ              
heat bus 1               inf      PQ              
electric bus 2           inf      PQ              
heat bus 2               inf      PQ              
In [5]:
network.buses["carrier"].value_counts()
Out[5]:
heat    3
AC      3
Name: carrier, dtype: int64
In [6]:
#add three lines in a ring
for i in range(3):
    network.add("Line","line {}".format(i),
                bus0="electric bus {}".format(i),
                bus1="electric bus {}".format((i+1)%3),
                x=0.1,
                s_nom=1000)
In [7]:
print(network.lines)
attribute            bus0            bus1 type    x    r    g    b   s_nom  \
line 0     electric bus 0  electric bus 1       0.1  0.0  0.0  0.0  1000.0   
line 1     electric bus 1  electric bus 2       0.1  0.0  0.0  0.0  1000.0   
line 2     electric bus 2  electric bus 0       0.1  0.0  0.0  0.0  1000.0   

attribute  s_nom_extendable  s_nom_min    ...      terrain_factor  \
line 0                False        0.0    ...                 1.0   
line 1                False        0.0    ...                 1.0   
line 2                False        0.0    ...                 1.0   

attribute  num_parallel  v_ang_min  v_ang_max  sub_network  x_pu  r_pu g_pu  \
line 0              1.0       -inf        inf                0.0   0.0  0.0   
line 1              1.0       -inf        inf                0.0   0.0  0.0   
line 2              1.0       -inf        inf                0.0   0.0  0.0   

attribute  b_pu  s_nom_opt  
line 0      0.0        0.0  
line 1      0.0        0.0  
line 2      0.0        0.0  

[3 rows x 23 columns]
In [8]:
#connect the electric to the heat buses with heat pumps with COP 3.
for i in range(3):
    network.add("Link",
                "heat pump {}".format(i),
                bus0="electric bus {}".format(i),
                bus1="heat bus {}".format(i),
                p_nom=100,
                efficiency=3.)
In [9]:
print(network.links)
attribute              bus0        bus1 type  efficiency  p_nom  \
heat pump 0  electric bus 0  heat bus 0              3.0  100.0   
heat pump 1  electric bus 1  heat bus 1              3.0  100.0   
heat pump 2  electric bus 2  heat bus 2              3.0  100.0   

attribute    p_nom_extendable  p_nom_min  p_nom_max  p_set  p_min_pu  \
heat pump 0             False        0.0        inf    0.0       0.0   
heat pump 1             False        0.0        inf    0.0       0.0   
heat pump 2             False        0.0        inf    0.0       0.0   

attribute    p_max_pu  capital_cost  marginal_cost  length  terrain_factor  \
heat pump 0       1.0           0.0            0.0     0.0             1.0   
heat pump 1       1.0           0.0            0.0     0.0             1.0   
heat pump 2       1.0           0.0            0.0     0.0             1.0   

attribute    p_nom_opt  
heat pump 0        0.0  
heat pump 1        0.0  
heat pump 2        0.0  
In [10]:
#add carriers
network.add("Carrier","gas",
           co2_emissions=0.27)
network.add("Carrier","biomass",
           co2_emissions=0.)
In [11]:
print(network.carriers)
attribute  co2_emissions
gas                 0.27
biomass             0.00
In [12]:
#add a gas generator at bus 0
network.add("Generator","gas generator",
            bus="electric bus 0",
            p_nom=100,
            marginal_cost=50,
            carrier="gas",
            efficiency=0.3)

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

#add a boiler at all heat buses
for i in range(3):
    network.add("Generator","boiler {}".format(i),
            bus="heat bus {}".format(i),
            p_nom=1000,
            efficiency=0.9,
            marginal_cost=20.,
            carrier="gas")
In [13]:
print(network.generators)
attribute                     bus control type   p_nom  p_nom_extendable  \
gas generator      electric bus 0      PQ        100.0             False   
biomass generator  electric bus 1      PQ        100.0             False   
boiler 0               heat bus 0      PQ       1000.0             False   
boiler 1               heat bus 1      PQ       1000.0             False   
boiler 2               heat bus 2      PQ       1000.0             False   

attribute          p_nom_min  p_nom_max  p_min_pu  p_max_pu  p_set    ...      \
gas generator            0.0        inf       0.0       1.0    0.0    ...       
biomass generator        0.0        inf       0.0       1.0    0.0    ...       
boiler 0                 0.0        inf       0.0       1.0    0.0    ...       
boiler 1                 0.0        inf       0.0       1.0    0.0    ...       
boiler 2                 0.0        inf       0.0       1.0    0.0    ...       

attribute          start_up_cost  shut_down_cost min_up_time  min_down_time  \
gas generator                0.0             0.0           0              0   
biomass generator            0.0             0.0           0              0   
boiler 0                     0.0             0.0           0              0   
boiler 1                     0.0             0.0           0              0   
boiler 2                     0.0             0.0           0              0   

attribute          initial_status  ramp_limit_up  ramp_limit_down  \
gas generator                   1            NaN              NaN   
biomass generator               1            NaN              NaN   
boiler 0                        1            NaN              NaN   
boiler 1                        1            NaN              NaN   
boiler 2                        1            NaN              NaN   

attribute          ramp_limit_start_up  ramp_limit_shut_down  p_nom_opt  
gas generator                      1.0                   1.0        0.0  
biomass generator                  1.0                   1.0        0.0  
boiler 0                           1.0                   1.0        0.0  
boiler 1                           1.0                   1.0        0.0  
boiler 2                           1.0                   1.0        0.0  

[5 rows x 27 columns]
In [14]:
#add electric loads
for i in range(3):
    network.add("Load","electric load {}".format(i),
                bus="electric bus {}".format(i),
                p_set=i*10)
In [15]:
#add heat loads
for i in range(3):
    network.add("Load","heat load {}".format(i),
                bus="heat bus {}".format(i),
                p_set=(3-i)*10)
In [16]:
print(network.loads)
attribute                   bus type  p_set  q_set  sign
electric load 0  electric bus 0         0.0    0.0  -1.0
electric load 1  electric bus 1        10.0    0.0  -1.0
electric load 2  electric bus 2        20.0    0.0  -1.0
heat load 0          heat bus 0        30.0    0.0  -1.0
heat load 1          heat bus 1        20.0    0.0  -1.0
heat load 2          heat bus 2        10.0    0.0  -1.0
In [17]:
print(network.loads.p_set)
electric load 0     0.0
electric load 1    10.0
electric load 2    20.0
heat load 0        30.0
heat load 1        20.0
heat load 2        10.0
Name: p_set, dtype: float64
In [18]:
#function for the LOPF

def run_lopf():
    network.lopf(keep_files=True)
    print("Objective:",network.objective)
    print(network.generators_t.p)
    print(network.links_t.p0)
    print(network.loads_t.p)
In [19]:
run_lopf()
INFO:pypsa.pf:Slack bus for sub-network 0 is electric bus 0
INFO:pypsa.pf:Slack bus for sub-network 1 is heat bus 0
INFO:pypsa.pf:Slack bus for sub-network 2 is heat bus 1
INFO:pypsa.pf:Slack bus for sub-network 3 is heat bus 2
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `angles` formulation
INFO:pypsa.opf:Solving model using glpk
Solver log file: '/tmp/tmplq7c4fy5.glpk.log'
Solver solution file: '/tmp/tmpx6nsweo0.glpk.raw'
Solver problem files: ('/tmp/tmp9k8pbjk9.pyomo.lp',)
INFO:pypsa.opf:Optimization successful
Objective: 2500.0
     gas generator  biomass generator  boiler 0  boiler 1  boiler 2
now           50.0                0.0       0.0       0.0       0.0
     heat pump 0  heat pump 1  heat pump 2
now         10.0     6.666667     3.333333
     electric load 0  electric load 1  electric load 2  heat load 0  \
now              0.0             10.0             20.0         30.0   

     heat load 1  heat load 2  
now         20.0         10.0  
In [20]:
#rerun with marginal costs for the heat pump operation

network.links.marginal_cost = 10
run_lopf()
INFO:pypsa.pf:Slack bus for sub-network 0 is electric bus 0
INFO:pypsa.pf:Slack bus for sub-network 1 is heat bus 0
INFO:pypsa.pf:Slack bus for sub-network 2 is heat bus 1
INFO:pypsa.pf:Slack bus for sub-network 3 is heat bus 2
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `angles` formulation
INFO:pypsa.opf:Solving model using glpk
Solver log file: '/tmp/tmp5ic7mgaa.glpk.log'
Solver solution file: '/tmp/tmppp9gux18.glpk.raw'
Solver problem files: ('/tmp/tmpq9gz83rp.pyomo.lp',)
INFO:pypsa.opf:Optimization successful
Objective: 2700.00000000001
     gas generator  biomass generator  boiler 0  boiler 1  boiler 2
now           30.0                0.0      30.0      20.0      10.0
     heat pump 0  heat pump 1  heat pump 2
now          0.0          0.0          0.0
     electric load 0  electric load 1  electric load 2  heat load 0  \
now              0.0             10.0             20.0         30.0   

     heat load 1  heat load 2  
now         20.0         10.0  
In [21]:
#rerun with no CO2 emissions

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

run_lopf()
INFO:pypsa.pf:Slack bus for sub-network 0 is electric bus 0
INFO:pypsa.pf:Slack bus for sub-network 1 is heat bus 0
INFO:pypsa.pf:Slack bus for sub-network 2 is heat bus 1
INFO:pypsa.pf:Slack bus for sub-network 3 is heat bus 2
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `angles` formulation
INFO:pypsa.opf:Solving model using glpk
Solver log file: '/tmp/tmpdlz6emja.glpk.log'
Solver solution file: '/tmp/tmpx7ep6s8i.glpk.raw'
Solver problem files: ('/tmp/tmp1__z1jsk.pyomo.lp',)
INFO:pypsa.opf:Optimization successful
Objective: 5200.0
     gas generator  biomass generator  boiler 0  boiler 1  boiler 2
now            0.0               50.0       0.0       0.0       0.0
     heat pump 0  heat pump 1  heat pump 2
now         10.0     6.666667     3.333333
     electric load 0  electric load 1  electric load 2  heat load 0  \
now              0.0             10.0             20.0         30.0   

     heat load 1  heat load 2  
now         20.0         10.0