LOPF then non-linear power flow with SciGRID

This Jupyter Notebook is also available to download at: http://www.pypsa.org/examples/scigrid-lopf-then-pf.ipynb and can be viewed as an HTML page at: http://pypsa.org/examples/scigrid-lopf-then-pf.html.

In this example, the dispatch of generators is optimised using the linear OPF, then a non-linear power flow is run on the resulting dispatch.

The data files for this example are in the examples folder of the github repository: https://github.com/PyPSA/PyPSA.

Data sources

Grid: based on SciGRID Version 0.2 which is based on OpenStreetMap.

Load size and location: based on Landkreise (NUTS 3) GDP and population.

Load time series: from ENTSO-E hourly data, scaled up uniformly by factor 1.12 (a simplification of the methodology in Schumacher, Hirth (2015)).

Conventional power plant capacities and locations: BNetzA list.

Wind and solar capacities and locations: EEG Stammdaten, based on http://www.energymap.info/download.html, which represents capacities at the end of 2014. Units without PLZ are removed.

Wind and solar time series: REatlas, Andresen et al, "Validation of Danish wind time series from a new global renewable energy atlas for energy system analysis," Energy 93 (2015) 1074 - 1088.

NB:

All times in the dataset are UTC.

Where SciGRID nodes have been split into 220kV and 380kV substations, all load and generation is attached to the 220kV substation.

Warnings

This script and the data behind it are no longer supported. See https://github.com/PyPSA/pypsa-eur for a newer model that covers the whole of Europe.

This dataset is ONLY intended to demonstrate the capabilities of PyPSA and is NOT (yet) accurate enough to be used for research purposes.

Known problems include:

i) Rough approximations have been made for missing grid data, e.g. 220kV-380kV transformers and connections between close sub-stations missing from OSM.

ii) There appears to be some unexpected congestion in parts of the network, which may mean for example that the load attachment method (by Voronoi cell overlap with Landkreise) isn't working, particularly in regions with a high density of substations.

iii) Attaching power plants to the nearest high voltage substation may not reflect reality.

iv) There is no proper n-1 security in the calculations - this can either be simulated with a blanket e.g. 70% reduction in thermal limits (as done here) or a proper security constrained OPF (see e.g. http://www.pypsa.org/examples/scigrid-sclopf.ipynb).

v) The borders and neighbouring countries are not represented.

vi) Hydroelectric power stations are not modelled accurately.

viii) The marginal costs are illustrative, not accurate.

ix) Only the first day of 2011 is in the github dataset, which is not representative. The full year of 2011 can be downloaded at http://www.pypsa.org/examples/scigrid-with-load-gen-trafos-2011.zip.

x) The ENTSO-E total load for Germany may not be scaled correctly; it is scaled up uniformly by factor 1.12 (a simplification of the methodology in Schumacher, Hirth (2015), which suggests monthly factors).

xi) Biomass from the EEG Stammdaten are not read in at the moment.

xii) Power plant start up costs, ramping limits/costs, minimum loading rates are not considered.

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

import pypsa

import numpy as np

import pandas as pd

import os

import matplotlib.pyplot as plt

import cartopy.crs as ccrs

%matplotlib inline
In [2]:
#You may have to adjust this path to where 
#you downloaded the github repository
#https://github.com/PyPSA/PyPSA

csv_folder_name = os.path.dirname(pypsa.__file__) + "/../examples/scigrid-de/scigrid-with-load-gen-trafos/"

network = pypsa.Network(csv_folder_name=csv_folder_name)
The argument csv_folder_name for initiating Network() is deprecated, please use import_name instead.
WARNING:pypsa.io:
Importing PyPSA from older version of PyPSA than current version 0.15.0.
Please read the release notes at https://pypsa.org/doc/release_notes.html
carefully to prepare your network for import.

INFO:pypsa.io:Imported network  has buses, generators, lines, loads, storage_units, transformers

Plot the distribution of the load and of generating tech

In [3]:
fig,ax = plt.subplots(1,1,subplot_kw={"projection":ccrs.PlateCarree()})

fig.set_size_inches(6,6)

load_distribution = network.loads_t.p_set.loc[network.snapshots[0]].groupby(network.loads.bus).sum()

network.plot(bus_sizes=0.5*load_distribution,ax=ax,title="Load distribution")
Out[3]:
(<matplotlib.collections.PathCollection at 0x7fdc11b798d0>,
 <matplotlib.collections.LineCollection at 0x7fdc11b79e10>)
In [4]:
fig.tight_layout()
#fig.savefig('load-distribution.png')
In [5]:
network.generators.groupby("carrier")["p_nom"].sum()
Out[5]:
carrier
Brown Coal       20879.500000
Gas              23913.130000
Geothermal          31.700000
Hard Coal        25312.600000
Multiple           152.700000
Nuclear          12068.000000
Oil               2710.200000
Other             3027.800000
Run of River      3999.100000
Solar            37041.524779
Storage Hydro     1445.000000
Waste             1645.900000
Wind Offshore     2973.500000
Wind Onshore     37339.895329
Name: p_nom, dtype: float64
In [6]:
network.storage_units.groupby("carrier")["p_nom"].sum()
Out[6]:
carrier
Pumped Hydro    9179.5
Name: p_nom, dtype: float64
In [7]:
techs = ["Gas","Brown Coal","Hard Coal","Wind Offshore","Wind Onshore","Solar"]

n_graphs = len(techs)

n_cols = 3

if n_graphs % n_cols == 0:
    n_rows = n_graphs // n_cols
else:
    n_rows = n_graphs // n_cols + 1

    
fig, axes = plt.subplots(nrows=n_rows, ncols=n_cols, subplot_kw={"projection":ccrs.PlateCarree()})

size = 4

fig.set_size_inches(size*n_cols,size*n_rows)

for i,tech in enumerate(techs):
    i_row = i // n_cols
    i_col = i % n_cols
    
    ax = axes[i_row,i_col]
    
    gens = network.generators[network.generators.carrier == tech]
    
    gen_distribution = gens.groupby("bus").sum()["p_nom"].reindex(network.buses.index,fill_value=0.)
    
    network.plot(ax=ax,bus_sizes=0.2*gen_distribution)
    
    ax.set_title(tech)
    
    

Run Linear Optimal Power Flow on the first day of 2011

In [8]:
#to approximate n-1 security and allow room for reactive power flows,
#don't allow any line to be loaded above 70% of their thermal rating

contingency_factor = 0.7

network.lines.s_max_pu = contingency_factor
In [9]:
#There are some infeasibilities without small extensions                                                                                 
network.lines.loc[["316","527","602"],"s_nom"] = 1715


#the lines to extend to resolve infeasibilities can
#be found by
#uncommenting the lines below to allow the network to be extended

#network.lines["s_nom_original"] = network.lines.s_nom

#network.lines.s_nom_extendable = True
#network.lines.s_nom_min = network.lines.s_nom

#Assume 450 EUR/MVA/km
#network.lines.capital_cost = 450*network.lines.length
In [10]:
group_size = 4

solver_name = "cbc"

print("Performing linear OPF for one day, {} snapshots at a time:".format(group_size))

network.storage_units.state_of_charge_initial = 0.

for i in range(int(24/group_size)):
    #set the initial state of charge based on previous round
    if i>0:
        network.storage_units.state_of_charge_initial = network.storage_units_t.state_of_charge.loc[network.snapshots[group_size*i-1]]
    network.lopf(network.snapshots[group_size*i:group_size*i+group_size],
                 solver_name=solver_name,
                 keep_files=True)
    #network.lines.s_nom = network.lines.s_nom_opt
Performing linear OPF for one day, 4 snapshots at a time:
INFO:pypsa.pf:Slack bus for sub-network 0 is 1
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `angles` formulation
INFO:pypsa.opf:Solving model using cbc
Solver log file: '/tmp/tmpo_0d1pdp.cbc.log'
Solver solution file: '/tmp/tmp4qyatsk_.pyomo.soln'
Solver problem files: ('/tmp/tmp4qyatsk_.pyomo.lp',)
INFO:pypsa.opf:Optimization successful
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 1449687.25
  Upper bound: 1449687.25
  Number of objectives: 1
  Number of constraints: 14025
  Number of variables: 12281
  Number of nonzeros: 1780
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.91
  Wallclock time: 0.91
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: None
      Number of created subproblems: None
    Black box: 
      Number of iterations: 5335
  Error rc: 0
  Time: 0.9351990222930908
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
INFO:pypsa.pf:Slack bus for sub-network 0 is 1
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `angles` formulation
INFO:pypsa.opf:Solving model using cbc
Solver log file: '/tmp/tmpsxktaazk.cbc.log'
Solver solution file: '/tmp/tmpsmhk1j99.pyomo.soln'
Solver problem files: ('/tmp/tmpsmhk1j99.pyomo.lp',)
INFO:pypsa.opf:Optimization successful
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 873763.9023
  Upper bound: 873763.9023
  Number of objectives: 1
  Number of constraints: 14025
  Number of variables: 12281
  Number of nonzeros: 1772
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.93
  Wallclock time: 0.93
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: None
      Number of created subproblems: None
    Black box: 
      Number of iterations: 5752
  Error rc: 0
  Time: 0.9551925659179688
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
INFO:pypsa.pf:Slack bus for sub-network 0 is 1
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `angles` formulation
INFO:pypsa.opf:Solving model using cbc
Solver log file: '/tmp/tmpvrv2zxjt.cbc.log'
Solver solution file: '/tmp/tmpohkds5y3.pyomo.soln'
Solver problem files: ('/tmp/tmpohkds5y3.pyomo.lp',)
INFO:pypsa.opf:Optimization successful
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 790786.5949
  Upper bound: 790786.5949
  Number of objectives: 1
  Number of constraints: 14025
  Number of variables: 12281
  Number of nonzeros: 1769
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.89
  Wallclock time: 0.9
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: None
      Number of created subproblems: None
    Black box: 
      Number of iterations: 5476
  Error rc: 0
  Time: 0.9228682518005371
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
INFO:pypsa.pf:Slack bus for sub-network 0 is 1
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `angles` formulation
INFO:pypsa.opf:Solving model using cbc
Solver log file: '/tmp/tmpaxfuxuv0.cbc.log'
Solver solution file: '/tmp/tmpd0ybe96k.pyomo.soln'
Solver problem files: ('/tmp/tmpd0ybe96k.pyomo.lp',)
INFO:pypsa.opf:Optimization successful
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 1455232.678
  Upper bound: 1455232.678
  Number of objectives: 1
  Number of constraints: 14025
  Number of variables: 12281
  Number of nonzeros: 1767
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.77
  Wallclock time: 0.78
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: None
      Number of created subproblems: None
    Black box: 
      Number of iterations: 5278
  Error rc: 0
  Time: 0.8118035793304443
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
INFO:pypsa.pf:Slack bus for sub-network 0 is 1
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `angles` formulation
INFO:pypsa.opf:Solving model using cbc
Solver log file: '/tmp/tmpv_lsuq8o.cbc.log'
Solver solution file: '/tmp/tmpkdpd3_8z.pyomo.soln'
Solver problem files: ('/tmp/tmpkdpd3_8z.pyomo.lp',)
INFO:pypsa.opf:Optimization successful
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 2647879.912
  Upper bound: 2647879.912
  Number of objectives: 1
  Number of constraints: 14025
  Number of variables: 12281
  Number of nonzeros: 1761
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 1.13
  Wallclock time: 1.15
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: None
      Number of created subproblems: None
    Black box: 
      Number of iterations: 5972
  Error rc: 0
  Time: 1.1760950088500977
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
INFO:pypsa.pf:Slack bus for sub-network 0 is 1
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `angles` formulation
INFO:pypsa.opf:Solving model using cbc
Solver log file: '/tmp/tmpr3h3xeug.cbc.log'
Solver solution file: '/tmp/tmpq682caxf.pyomo.soln'
Solver problem files: ('/tmp/tmpq682caxf.pyomo.lp',)
INFO:pypsa.opf:Optimization successful
# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 2142594.539
  Upper bound: 2142594.539
  Number of objectives: 1
  Number of constraints: 14025
  Number of variables: 12281
  Number of nonzeros: 1769
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.89
  Wallclock time: 0.9
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: None
      Number of created subproblems: None
    Black box: 
      Number of iterations: 5542
  Error rc: 0
  Time: 0.9249532222747803
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
In [11]:
#if lines are extended, look at which ones are bigger
#network.lines[["s_nom_original","s_nom"]][abs(network.lines.s_nom - network.lines.s_nom_original) > 1]
In [12]:
p_by_carrier = network.generators_t.p.groupby(network.generators.carrier, axis=1).sum()
In [13]:
p_by_carrier.drop((p_by_carrier.max()[p_by_carrier.max() < 1700.]).index,axis=1,inplace=True)
In [14]:
p_by_carrier.columns
Out[14]:
Index(['Brown Coal', 'Gas', 'Hard Coal', 'Nuclear', 'Run of River', 'Solar',
       'Wind Offshore', 'Wind Onshore'],
      dtype='object', name='carrier')
In [15]:
colors = {"Brown Coal" : "brown",
          "Hard Coal" : "k",
          "Nuclear" : "r",
          "Run of River" : "green",
          "Wind Onshore" : "blue",
          "Solar" : "yellow",
          "Wind Offshore" : "cyan",
          "Waste" : "orange",
          "Gas" : "orange"}
#reorder
cols = ["Nuclear","Run of River","Brown Coal","Hard Coal","Gas","Wind Offshore","Wind Onshore","Solar"]
p_by_carrier = p_by_carrier[cols]
In [16]:
fig,ax = plt.subplots(1,1)

fig.set_size_inches(12,6)

(p_by_carrier/1e3).plot(kind="area",ax=ax,linewidth=4,colors=[colors[col] for col in p_by_carrier.columns])


ax.legend(ncol=4,loc="upper left")

ax.set_ylabel("GW")

ax.set_xlabel("")
/home/tom/miniconda3/envs/pypsa-eur/lib/python3.6/site-packages/pandas/plotting/_matplotlib/core.py:203: UserWarning: 'colors' is being deprecated. Please use 'color'instead of 'colors'
  "'colors' is being deprecated. Please use 'color'"
Out[16]:
Text(0.5, 0, '')
In [17]:
fig.tight_layout()
#fig.savefig("stacked-gen.png")
In [18]:
fig,ax = plt.subplots(1,1)
fig.set_size_inches(12,6)

p_storage = network.storage_units_t.p.sum(axis=1)
state_of_charge = network.storage_units_t.state_of_charge.sum(axis=1)
p_storage.plot(label="Pumped hydro dispatch",ax=ax,linewidth=3)
state_of_charge.plot(label="State of charge",ax=ax,linewidth=3)

ax.legend()
ax.grid()
ax.set_ylabel("MWh")
ax.set_xlabel("")
Out[18]:
Text(0.5, 0, '')
In [19]:
fig.tight_layout()
#fig.savefig("storage-scigrid.png")
In [20]:
now = network.snapshots[4]
In [21]:
print("With the linear load flow, there is the following per unit loading:")
loading = network.lines_t.p0.loc[now]/network.lines.s_nom
print(loading.describe())
With the linear load flow, there is the following per unit loading:
count    852.000000
mean      -0.002815
std        0.258588
min       -0.700000
25%       -0.131007
50%        0.003142
75%        0.122643
max        0.700000
dtype: float64
In [22]:
fig,ax = plt.subplots(1,1,subplot_kw={"projection":ccrs.PlateCarree()})
fig.set_size_inches(6,6)

network.plot(ax=ax,line_colors=abs(loading),line_cmap=plt.cm.jet,title="Line loading")
Out[22]:
(<matplotlib.collections.PathCollection at 0x7fdc0e8811d0>,
 <matplotlib.collections.LineCollection at 0x7fdc0e86d668>)
In [23]:
fig.tight_layout()
#fig.savefig("line-loading.png")
In [24]:
network.buses_t.marginal_price.loc[now].describe()
Out[24]:
count    585.000000
mean      15.737598
std       10.941995
min      -10.397824
25%        6.992120
50%       15.841190
75%       25.048186
max       52.150120
Name: 2011-01-01 04:00:00, dtype: float64
In [25]:
fig,ax = plt.subplots(1,1,subplot_kw={"projection":ccrs.PlateCarree()})
fig.set_size_inches(6,4)


network.plot(ax=ax,line_widths=pd.Series(0.5,network.lines.index))
plt.hexbin(network.buses.x, network.buses.y, 
           gridsize=20,
           C=network.buses_t.marginal_price.loc[now],
           cmap=plt.cm.jet)

#for some reason the colorbar only works with graphs plt.plot
#and must be attached plt.colorbar

cb = plt.colorbar()
cb.set_label('Locational Marginal Price (EUR/MWh)') 
In [26]:
fig.tight_layout()
#fig.savefig('lmp.png')

Look at variable curtailment

In [27]:
carrier = "Wind Onshore"

capacity = network.generators.groupby("carrier").sum().at[carrier,"p_nom"]
In [28]:
p_available = network.generators_t.p_max_pu.multiply(network.generators["p_nom"])
In [29]:
p_available_by_carrier =p_available.groupby(network.generators.carrier, axis=1).sum()
In [30]:
p_curtailed_by_carrier = p_available_by_carrier - p_by_carrier
In [31]:
p_df = pd.DataFrame({carrier + " available" : p_available_by_carrier[carrier],
                     carrier + " dispatched" : p_by_carrier[carrier],
                     carrier + " curtailed" : p_curtailed_by_carrier[carrier]})

p_df[carrier + " capacity"] = capacity
In [32]:
p_df["Wind Onshore curtailed"][p_df["Wind Onshore curtailed"] < 0.] = 0.
In [33]:
fig,ax = plt.subplots(1,1)
fig.set_size_inches(12,6)
p_df[[carrier + " dispatched",carrier + " curtailed"]].plot(kind="area",ax=ax,linewidth=3)
p_df[[carrier + " available",carrier + " capacity"]].plot(ax=ax,linewidth=3)

ax.set_xlabel("")
ax.set_ylabel("Power [MW]")
ax.set_ylim([0,40000])
ax.legend()
Out[33]:
<matplotlib.legend.Legend at 0x7fdc0f445fd0>