Script to add load, generators, missing lines and transformers to SciGRID

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

This script does some post-processing on the original SciGRID dataset version 0.2 and then adds load, generation, transformers and missing lines to the SciGRID dataset.

The intention is to create a model of the German electricity system that is transparent in the sense that all steps from openly-available raw data to the final model can be followed. The model is NOT validated and may contain errors.

Some of the libraries used for attaching the load and generation are not on github, but can be downloaded at

http://fias.uni-frankfurt.de/~hoersch/

The intention is to release these as free software soon. We cannot guarantee to support you when using these libraries.

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.

Warning

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 pandas as pd

import numpy as np

from six import iteritems
from six.moves import range

import os

import matplotlib.pyplot as plt

%matplotlib inline

Read in the raw SciGRID data

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

folder_prefix = os.path.dirname(pypsa.__file__) + "/../examples/scigrid-de/"
In [3]:
#note that some columns have 'quotes because of fields containing commas'
vertices = pd.read_csv(folder_prefix+"scigrid-151109/vertices_de_power_151109.csvdata",sep=",",quotechar="'",index_col=0)

vertices.rename(columns={"lon":"x","lat":"y","name":"osm_name"},inplace=True)
In [4]:
print(vertices["voltage"].value_counts(dropna=False))
220000;110000                 129
380000;110000                 101
380000;220000;110000           64
220000                         62
380000                         58
NaN                            34
380000;220000                  24
220000;65000                    3
110000                          3
400000;220000                   2
400000;220000;110000            2
400000                          2
150000                          2
380000;110000;20000             1
400000;220000;150000            1
400000;130000                   1
400000;225000                   1
380000;150000;110000            1
400000;225000;90000             1
220000;150000;60000             1
220000;155000                   1
380000;220000;110000;20000      1
Name: voltage, dtype: int64
In [5]:
links = pd.read_csv(folder_prefix+"scigrid-151109/links_de_power_151109.csvdata",sep=",",quotechar="'",index_col=0)
links.rename(columns={"v_id_1":"bus0","v_id_2":"bus1","name":"osm_name"},inplace=True)

links["cables"].fillna(3,inplace=True)
links["wires"].fillna(2,inplace=True)

links["length"] = links["length_m"]/1000.
In [6]:
print(links["voltage"].value_counts(dropna=False))
380000    422
220000    399
300000      2
400000      1
450000      1
Name: voltage, dtype: int64
In [7]:
## Drop the DC lines

for voltage in [300000,400000,450000]:
    links.drop(links[links.voltage == voltage].index,inplace=True)
In [8]:
## Build the network

network = pypsa.Network()

pypsa.io.import_components_from_dataframe(network,vertices,"Bus")

pypsa.io.import_components_from_dataframe(network,links,"Line")

Add specific missing AC lines

In [11]:
# Add AC lines known to be missing in SciGRID                                                                                                             
# E.g. lines missing because of OSM mapping errors.                                                                                                       
# This is no systematic list, just what we noticed;                                                                                                       
# please tell SciGRID and/or Tom Brown (brown@fias.uni-frankfurt.de)                                                                                      
# if you know of more examples                                                                                                                            

columns = ["bus0","bus1","wires","cables","voltage"]

data = [["100","255",2,6,220000], # Niederstedem to Wengerohr                                                                                             
        ["384","351",4,6,380000], # Raitersaich to Ingolstadt                                                                                             
        ["351","353",4,6,380000], # Ingolstadt to Irsching                                                                                                
        ]

last_scigrid_line = int(network.lines.index[-1])

index = [str(i) for i in range(last_scigrid_line+1,last_scigrid_line+1 + len(data))]

missing_lines = pd.DataFrame(data,index,columns)

#On average, SciGRID lines are 25% longer than the direct distance
length_factor = 1.25

missing_lines["length"] = [length_factor*pypsa.geo.haversine(network.buses.loc[r.bus0,["x","y"]],network.buses.loc[r.bus1,["x","y"]])[0,0] for i,r in missing_lines.iterrows()]
In [12]:
pypsa.io.import_components_from_dataframe(network,missing_lines,"Line")
In [13]:
network.lines.tail()
Out[13]:
b b_pu bus0 bus1 c_nfkm cables capital_cost frequency from_relation g ... terrain_factor type v_ang_max v_ang_min voltage wires wkt_srid_4326 x x_ohmkm x_pu
824 0.0 0.0 333 495 13.7 6.0 0.0 50.0 5502026.0 0.0 ... 1.0 inf -inf 380000.0 2.0 SRID=4326;LINESTRING(9.19698235871826 54.85109... 0.0 0.25 0.0
825 0.0 0.0 398 495 13.7 6.0 0.0 50.0 5502026.0 0.0 ... 1.0 inf -inf 380000.0 2.0 SRID=4326;LINESTRING(9.72460611042234 54.28827... 0.0 0.25 0.0
826 0.0 0.0 100 255 NaN 6.0 0.0 NaN NaN 0.0 ... 1.0 inf -inf 220000.0 2.0 NaN 0.0 NaN 0.0
827 0.0 0.0 384 351 NaN 6.0 0.0 NaN NaN 0.0 ... 1.0 inf -inf 380000.0 4.0 NaN 0.0 NaN 0.0
828 0.0 0.0 351 353 NaN 6.0 0.0 NaN NaN 0.0 ... 1.0 inf -inf 380000.0 4.0 NaN 0.0 NaN 0.0

5 rows × 37 columns

Determine the voltage of the buses by the lines which end there

In [14]:
network.lines.voltage.value_counts()
Out[14]:
380000.0    424
220000.0    400
Name: voltage, dtype: int64
In [15]:
buses_by_voltage = {}

for voltage in network.lines.voltage.value_counts().index:
    buses_by_voltage[voltage] = set(network.lines[network.lines.voltage == voltage].bus0)\
                                | set(network.lines[network.lines.voltage == voltage].bus1)
In [16]:
# give priority to 380 kV
network.buses["v_nom"] = 380
network.buses.loc[buses_by_voltage[220000],"v_nom"] = 220
network.buses.loc[buses_by_voltage[380000],"v_nom"] = 380
In [17]:
network.buses.v_nom.value_counts(dropna=False)
Out[17]:
380    277
220    218
Name: v_nom, dtype: int64

Connect buses which are < 850m apart

There are pairs of buses less than 850m apart which are not connected in SciGRID, but clearly connected in OpenStreetMap (OSM).

The reason is that the relations for connections between close substations do not appear in OSM.

Here they are connected with 2 circuits of the appropriate voltage level (an asumption).

850m is chosen as a limit based on manually looking through the examples.

The example 46-48 (Marzahn) at 892 m apart is the first example of close substations which are not connected in reality.

In [18]:
# Compute the distances for unique pairs

pairs = pd.Series()

for i,u in enumerate(network.buses.index):
    vs = network.buses[["x","y"]].iloc[i+1:]
    distance_km = pypsa.geo.haversine(vs,network.buses.loc[u,["x","y"]])

    to_add = pd.Series(data=distance_km[:,0],index=[(u,v) for v in vs.index])
    
    pairs = pd.concat((pairs,to_add))
In [19]:
pairs.sort_values().head()
Out[19]:
(443, 444)    0.024582
(465, 471)    0.038207
(483, 484)    0.075587
(250, 251)    0.089066
(419, 420)    0.100549
dtype: float64
In [20]:
# determine topology so we can look what's actually connected
network.determine_network_topology()
In [22]:
# Example all substations which are close to                                                                                                              
# each other geographically by not connected in network.adj                                                                                               

start = 0  #km                                                                                                                                            
stop = 1 #km                                                                                                                                              

for (u,v),dist in pairs.sort_values().iteritems():

    if dist < start:
        continue

    #only go up to pairs stop km apart                                                                                                                    
    if dist > stop:
        break

    #ignore if they're already connected                                                                                                                  
    if u in network.graph().adj[v]:
        continue


    print(u,v,dist)

    u_x = network.buses.at[u,"x"]
    u_y = network.buses.at[u,"y"]
    v_x = network.buses.at[v,"x"]
    v_y = network.buses.at[v,"y"]

    #have a look what's going on in OSM                                                                                                                   
    print("https://www.openstreetmap.org/#map=18/{}/{}".format(u_y,u_x))
    print("https://www.openstreetmap.org/#map=18/{}/{}".format(v_y,v_x))
465 471 0.0382066887159
https://www.openstreetmap.org/#map=18/50.8441529722/6.92909269081
https://www.openstreetmap.org/#map=18/50.8438169469/6.92897905392
250 251 0.0890664669456
https://www.openstreetmap.org/#map=18/51.6736430745/7.7155434195
https://www.openstreetmap.org/#map=18/51.6730996381/7.7145945335
419 420 0.100549382242
https://www.openstreetmap.org/#map=18/54.3540286882/6.02502872786
https://www.openstreetmap.org/#map=18/54.3549329471/6.02503294994
251 252 0.139286625328
https://www.openstreetmap.org/#map=18/51.6730996381/7.7145945335
https://www.openstreetmap.org/#map=18/51.6726636622/7.71270093431
194 195 0.162288711745
https://www.openstreetmap.org/#map=18/51.5771047777/6.6830656826
https://www.openstreetmap.org/#map=18/51.5756455059/6.68310697664
249 250 0.201052985128
https://www.openstreetmap.org/#map=18/51.6743422007/7.71823231963
https://www.openstreetmap.org/#map=18/51.6736430745/7.7155434195
2 296 0.217467718393
https://www.openstreetmap.org/#map=18/52.5438533224/9.11321007473
https://www.openstreetmap.org/#map=18/52.5457924749/9.11362795987
250 252 0.22423259083
https://www.openstreetmap.org/#map=18/51.6736430745/7.7155434195
https://www.openstreetmap.org/#map=18/51.6726636622/7.71270093431
54 367 0.237790544199
https://www.openstreetmap.org/#map=18/49.2524597289/8.43846751212
https://www.openstreetmap.org/#map=18/49.251402647/8.4413154852
131 267 0.255219620964
https://www.openstreetmap.org/#map=18/51.1225390519/6.80109333957
https://www.openstreetmap.org/#map=18/51.1247880606/6.8018236571
91 92 0.267503917553
https://www.openstreetmap.org/#map=18/52.4801491896/7.30561448538
https://www.openstreetmap.org/#map=18/52.4824259958/7.30433875539
249 251 0.286382404109
https://www.openstreetmap.org/#map=18/51.6743422007/7.71823231963
https://www.openstreetmap.org/#map=18/51.6730996381/7.7145945335
186 191 0.291957501887
https://www.openstreetmap.org/#map=18/51.5261875491/6.71114527446
https://www.openstreetmap.org/#map=18/51.5267908666/6.71525259692
129 173 0.29597806878
https://www.openstreetmap.org/#map=18/51.2606098249/6.62582790177
https://www.openstreetmap.org/#map=18/51.2587133681/6.62284323547
96 202 0.332431884295
https://www.openstreetmap.org/#map=18/52.4691114903/7.31784994378
https://www.openstreetmap.org/#map=18/52.4714572265/7.32089257542
53 322 0.349721512651
https://www.openstreetmap.org/#map=18/48.0993557006/7.7534954654
https://www.openstreetmap.org/#map=18/48.0964057189/7.75186249821
24 398 0.383041086878
https://www.openstreetmap.org/#map=18/54.2914209868/9.72699196993
https://www.openstreetmap.org/#map=18/54.2882702503/9.72460611042
133 156 0.409565029791
https://www.openstreetmap.org/#map=18/50.8635385009/6.84265066151
https://www.openstreetmap.org/#map=18/50.8598647244/6.84223115148
448 458 0.414353807553
https://www.openstreetmap.org/#map=18/51.4592583943/7.42474524551
https://www.openstreetmap.org/#map=18/51.4575100936/7.41946378046
451 458 0.415017456637
https://www.openstreetmap.org/#map=18/51.4607823587/7.41658261385
https://www.openstreetmap.org/#map=18/51.4575100936/7.41946378046
249 252 0.424643273603
https://www.openstreetmap.org/#map=18/51.6743422007/7.71823231963
https://www.openstreetmap.org/#map=18/51.6726636622/7.71270093431
187 195 0.439508350931
https://www.openstreetmap.org/#map=18/51.5764228882/6.68934277039
https://www.openstreetmap.org/#map=18/51.5756455059/6.68310697664
216 236 0.465178778054
https://www.openstreetmap.org/#map=18/48.9144168738/9.19078743692
https://www.openstreetmap.org/#map=18/48.912963438/9.19675652645
246 391 0.477174425603
https://www.openstreetmap.org/#map=18/50.0826078105/8.95870856177
https://www.openstreetmap.org/#map=18/50.0860157333/8.95464412134
396 410 0.520731065901
https://www.openstreetmap.org/#map=18/48.5099297353/10.4015868558
https://www.openstreetmap.org/#map=18/48.5146084081/10.4012813869
448 451 0.590359937034
https://www.openstreetmap.org/#map=18/51.4592583943/7.42474524551
https://www.openstreetmap.org/#map=18/51.4607823587/7.41658261385
448 450 0.60850191077
https://www.openstreetmap.org/#map=18/51.4592583943/7.42474524551
https://www.openstreetmap.org/#map=18/51.4570521583/7.41670791196
198 478 0.774034099174
https://www.openstreetmap.org/#map=18/51.6768446668/7.97065689371
https://www.openstreetmap.org/#map=18/51.6766515665/7.9818783186
470 473 0.781865929275
https://www.openstreetmap.org/#map=18/51.7500675543/7.82599043341
https://www.openstreetmap.org/#map=18/51.7550403465/7.83402074323
198 253 0.816599958644
https://www.openstreetmap.org/#map=18/51.6768446668/7.97065689371
https://www.openstreetmap.org/#map=18/51.6735262619/7.98122159643
112 143 0.836332808631
https://www.openstreetmap.org/#map=18/51.5758105021/7.3531118694
https://www.openstreetmap.org/#map=18/51.5683686081/7.35135784055
46 48 0.892118502762
https://www.openstreetmap.org/#map=18/52.5356342713/13.524198731
https://www.openstreetmap.org/#map=18/52.5374565238/13.5370442127
In [23]:
# From examining the map, it seems that all cases where substations                                                                                       
# are less than 850m apart are connected in reality                                                                                                       
# The first one to fail is 46-48 (Marzahn) at 892 m                                                                                                       

# Connect these substations                                                                                                                               

limit = 0.85

for (u,v),dist in pairs.sort_values().iteritems():

    #only go up to pairs stop km apart                                                                                                                    
    if dist > limit:
        break

    #ignore if they're already connected                                                                                                                  
    if u in network.graph().adj[v]:
        continue


    kv_u = network.buses.at[u,"v_nom"]
    kv_v = network.buses.at[v,"v_nom"]

    print(u,v,dist,kv_u,kv_v)
    
    last_scigrid_line = int(network.lines.index[-1])
    
    voltage = max(kv_u,kv_v)*1000
    
    wires = {220000 : 2, 380000 : 4}[voltage]
    
    cables = 6
    
    df = pd.DataFrame([[u,v,length_factor*dist,wires,cables,voltage]],columns=["bus0","bus1","length","wires","cables","voltage"],index=[str(last_scigrid_line+1)])

    pypsa.io.import_components_from_dataframe(network,df,"Line")
465 471 0.0382066887159 220 220
250 251 0.0890664669456 220 220
419 420 0.100549382242 380 380
251 252 0.139286625328 220 220
194 195 0.162288711745 380 220
249 250 0.201052985128 380 220
2 296 0.217467718393 380 380
250 252 0.22423259083 220 220
54 367 0.237790544199 380 380
131 267 0.255219620964 380 220
91 92 0.267503917553 380 380
249 251 0.286382404109 380 220
186 191 0.291957501887 380 220
129 173 0.29597806878 220 380
96 202 0.332431884295 220 380
53 322 0.349721512651 380 220
24 398 0.383041086878 220 380
133 156 0.409565029791 380 220
448 458 0.414353807553 380 220
451 458 0.415017456637 380 220
249 252 0.424643273603 380 220
187 195 0.439508350931 380 220
216 236 0.465178778054 380 220
246 391 0.477174425603 380 220
396 410 0.520731065901 380 380
448 451 0.590359937034 380 380
448 450 0.60850191077 380 380
198 478 0.774034099174 380 380
470 473 0.781865929275 380 380
198 253 0.816599958644 380 380
112 143 0.836332808631 380 220

Split buses with more than one voltage; add trafos between

This code splits the buses where you have 220 and 380 kV lines landing.

In [24]:
network.lines.voltage.value_counts()
Out[24]:
380000.0    451
220000.0    404
Name: voltage, dtype: int64
In [25]:
buses_by_voltage = {}

for voltage in network.lines.voltage.value_counts().index:
    buses_by_voltage[voltage] = set(network.lines[network.lines.voltage == voltage].bus0)\
                                | set(network.lines[network.lines.voltage == voltage].bus1)
In [26]:
network.buses.v_nom=380
network.buses.loc[buses_by_voltage[220000],"v_nom"] = 220
network.buses.loc[buses_by_voltage[380000],"v_nom"] = 380
In [27]:
overlap = buses_by_voltage[220000] & buses_by_voltage[380000]
len(overlap)
Out[27]:
96
In [28]:
## build up new buses and transformers to import


buses_to_split = [str(i) for i in sorted([int(item) for item in overlap])]
buses_to_split_df = network.buses.loc[buses_to_split]

buses_to_split_df.v_nom = 220

buses_to_split_220kV = [name + "_220kV" for name in buses_to_split_df.index]

buses_to_split_df.index = buses_to_split_220kV

trafos_df = pd.DataFrame(index=buses_to_split)
trafos_df["bus0"] = buses_to_split
trafos_df["bus1"] = buses_to_split_220kV
trafos_df["x"] = 0.1
#This high a nominal power is required for feasibility in LOPF
trafos_df["s_nom"] = 2000
In [29]:
pypsa.io.import_components_from_dataframe(network,buses_to_split_df,"Bus")
pypsa.io.import_components_from_dataframe(network,trafos_df,"Transformer")
In [30]:
##reconnect lines to the correct voltage bus

for line in network.lines.index:
    bus0 = network.lines.at[line,"bus0"]
    bus1 = network.lines.at[line,"bus1"]
    v0 = network.buses.at[bus0,"v_nom"]
    v1 = network.buses.at[bus1,"v_nom"]
    v = network.lines.at[line,"voltage"]
    if v0 != v/1000.:
        print(line,v0,v)
        network.lines.at[line,"bus0"] = bus0+"_220kV"
    if v1 != v/1000.:
        network.lines.at[line,"bus1"] = bus1+"_220kV"
3 380.0 220000.0
6 380.0 220000.0
8 380.0 220000.0
9 380.0 220000.0
16 380.0 220000.0
19 380.0 220000.0
32 380.0 220000.0
40 380.0 220000.0
41 380.0 220000.0
52 380.0 220000.0
53 380.0 220000.0
58 380.0 220000.0
78 380.0 220000.0
99 380.0 220000.0
127 380.0 220000.0
133 380.0 220000.0
140 380.0 220000.0
143 380.0 220000.0
145 380.0 220000.0
156 380.0 220000.0
157 380.0 220000.0
160 380.0 220000.0
161 380.0 220000.0
165 380.0 220000.0
166 380.0 220000.0
167 380.0 220000.0
170 380.0 220000.0
171 380.0 220000.0
172 380.0 220000.0
173 380.0 220000.0
176 380.0 220000.0
180 380.0 220000.0
193 380.0 220000.0
195 380.0 220000.0
196 380.0 220000.0
198 380.0 220000.0
199 380.0 220000.0
203 380.0 220000.0
206 380.0 220000.0
208 380.0 220000.0
209 380.0 220000.0
223 380.0 220000.0
224 380.0 220000.0
228 380.0 220000.0
229 380.0 220000.0
234 380.0 220000.0
235 380.0 220000.0
236 380.0 220000.0
237 380.0 220000.0
240 380.0 220000.0
243 380.0 220000.0
244 380.0 220000.0
262 380.0 220000.0
275 380.0 220000.0
278 380.0 220000.0
283 380.0 220000.0
306 380.0 220000.0
308 380.0 220000.0
314 380.0 220000.0
322 380.0 220000.0
323 380.0 220000.0
327 380.0 220000.0
333 380.0 220000.0
334 380.0 220000.0
335 380.0 220000.0
337 380.0 220000.0
338 380.0 220000.0
341 380.0 220000.0
342 380.0 220000.0
345 380.0 220000.0
350 380.0 220000.0
351 380.0 220000.0
363 380.0 220000.0
364 380.0 220000.0
365 380.0 220000.0
372 380.0 220000.0
373 380.0 220000.0
374 380.0 220000.0
375 380.0 220000.0
380 380.0 220000.0
381 380.0 220000.0
382 380.0 220000.0
393 380.0 220000.0
396 380.0 220000.0
413 380.0 220000.0
414 380.0 220000.0
418 380.0 220000.0
422 380.0 220000.0
424 380.0 220000.0
425 380.0 220000.0
427 380.0 220000.0
428 380.0 220000.0
434 380.0 220000.0
442 380.0 220000.0
463 380.0 220000.0
467 380.0 220000.0
473 380.0 220000.0
485 380.0 220000.0
486 380.0 220000.0
497 380.0 220000.0
503 380.0 220000.0
506 380.0 220000.0
507 380.0 220000.0
508 380.0 220000.0
509 380.0 220000.0
518 380.0 220000.0
520 380.0 220000.0
528 380.0 220000.0
530 380.0 220000.0
536 380.0 220000.0
539 380.0 220000.0
555 380.0 220000.0
564 380.0 220000.0
565 380.0 220000.0
566 380.0 220000.0
567 380.0 220000.0
575 380.0 220000.0
577 380.0 220000.0
583 380.0 220000.0
584 380.0 220000.0
601 380.0 220000.0
603 380.0 220000.0
604 380.0 220000.0
608 380.0 220000.0
609 380.0 220000.0
610 380.0 220000.0
612 380.0 220000.0
624 380.0 220000.0
626 380.0 220000.0
627 380.0 220000.0
629 380.0 220000.0
632 380.0 220000.0
634 380.0 220000.0
646 380.0 220000.0
648 380.0 220000.0
650 380.0 220000.0
654 380.0 220000.0
659 380.0 220000.0
660 380.0 220000.0
661 380.0 220000.0
669 380.0 220000.0
670 380.0 220000.0
673 380.0 220000.0
674 380.0 220000.0
678 380.0 220000.0
679 380.0 220000.0
680 380.0 220000.0
699 380.0 220000.0
705 380.0 220000.0
708 380.0 220000.0
716 380.0 220000.0
717 380.0 220000.0
720 380.0 220000.0
744 380.0 220000.0
745 380.0 220000.0
747 380.0 220000.0
748 380.0 220000.0
755 380.0 220000.0
757 380.0 220000.0
760 380.0 220000.0
766 380.0 220000.0
768 380.0 220000.0
774 380.0 220000.0
779 380.0 220000.0
780 380.0 220000.0
783 380.0 220000.0
803 380.0 220000.0
814 380.0 220000.0
819 380.0 220000.0
826 380.0 220000.0
830 380.0 220000.0
832 380.0 220000.0
836 380.0 220000.0
In [31]:
#determine the connected components

network.determine_network_topology()


#remove small isolated networks
for sn in network.sub_networks.obj:
    buses = sn.buses().index
    branches = sn.branches().index
    
    if len(buses) < 5:
        print("Dropping Sub-Network {} because it only has {} buses".format(sn,len(buses)))
        #print(buses.index)
        #print(len(branches),branches.index)
        for bus in buses:
            network.remove("Bus",bus)
        for branch in branches:
            network.remove("Line",branch[1])
    else:
        print("Keeping Sub-Network {} because it has {} buses".format(sn,len(buses)))

#rebuild topology

network.determine_network_topology()
Keeping Sub-Network SubNetwork 0 because it has 585 buses
Dropping Sub-Network SubNetwork 1 because it only has 2 buses
Dropping Sub-Network SubNetwork 2 because it only has 1 buses
Dropping Sub-Network SubNetwork 3 because it only has 1 buses
Dropping Sub-Network SubNetwork 4 because it only has 2 buses
In [29]:
colors = network.lines.voltage.map(lambda v: "g" if v == 220000 else "r" if v == 380000 else "c")

network.plot(line_colors=colors)
Out[29]:
(<matplotlib.collections.PathCollection at 0x7f7789af5d10>,
 <matplotlib.collections.LineCollection at 0x7f7789b199d0>)

Recalculate all electrical properties

In [30]:
network.lines["type"] = network.lines.voltage.map({220000 : "Al/St 240/40 2-bundle 220.0",
                                                   380000 : "Al/St 240/40 4-bundle 380.0"})

network.lines["num_parallel"] = network.lines.cables/3.*network.lines.wires/network.lines.voltage.map({220000 : 2., 380000 : 4.})

network.lines["s_nom"] = 3.**0.5*network.lines.voltage/1000.*network.lines.num_parallel*network.lines.voltage.map({220000 : 2., 380000 : 4.})*0.65

Attach the load

In [31]:
#import FIAS libraries for attaching data

#this script uses old versions of the FIAS libraries and 
#has not yet been updated to the new versions

#the latest versions are available at
#https://github.com/FRESNA/vresutils

#if you get it working with the new versions, please
#tell us! It shouldn't be too hard...

try:
    import vresutils, load
except:
    print("Oh dear! You don't have FIAS libraries, so you cannot add load :-(")
In [32]:
import load
In [33]:
from vresutils import graph as vgraph
from vresutils import shapes as vshapes
from vresutils import grid as vgrid
from vresutils import dispatch as vdispatch
from shapely.geometry import Polygon
from load import germany as DEload
import networkx as nx
In [34]:
#bounding poly for Germany for the Voronoi - necessary
#because some SciGRID points lie outside border vshapes.germany()
poly = Polygon([[5.8,47.],[5.8,55.5],[15.2,55.5],[15.2,47.]])
In [35]:
def generate_dummy_graph(network):
    """Generate a dummy graph to feed to the FIAS libraries.
    It adds the "pos" attribute and removes the 380 kV duplicate
    buses when the buses have been split, so that all load and generation
    is attached to the 220kV bus."""
    
    graph = pypsa.descriptors.OrderedGraph()
    
    graph.add_nodes_from([bus for bus in network.buses.index if bus not in buses_to_split])
    
    #add positions to graph for voronoi cell computation
    for node in graph.nodes():
        graph.node[node]["pos"] = np.array(network.buses.loc[node,["x","y"]],dtype=float)
    
    return graph
In [36]:
graph = generate_dummy_graph(network)
In [37]:
graph.name = "scigrid_v2"
In [38]:
def voronoi_partition(G, outline):
    """                                                                                                                                                   
    For 2D-embedded graph `G`, within the boundary given by the shapely polygon                                                                           
    `outline`, returns `G` with the Voronoi cell region as an additional node                                                                             
    attribute.                                                                                                                                            
    """
    #following line from vresutils.graph caused a bug
    #G = polygon_subgraph(G, outline, copy=False)
    points = list(vresutils.graph.get_node_attributes(G, 'pos').values())
    regions = vresutils.graph.voronoi_partition_pts(points, outline, no_multipolygons=True)
    nx.set_node_attributes(G, 'region', dict(zip(G.nodes(), regions)))

    return G
In [39]:
voronoi_partition(graph, poly)
Out[39]:
<pypsa.descriptors.OrderedGraph at 0x7fda831b3790>
In [40]:
#NB: starts at midnight CET, 23:00 UTC
load = DEload.timeseries(graph, years=[2011, 2012, 2013, 2014])
Serving call to read_entsoe_consumption from file load.europe.read_entsoe_consumption_2011_2012_2013_2014_countries.Germany.pickle of cache: 1.7 msec
Serving call to read_entsoe_consumption from file load.europe.read_entsoe_consumption__.pickle of cache: 15.6 msec
Serving call to landkreise from file load.germany.landkreise_ver5__.pickle of cache: 628.9 usec
Serving call to landkreise from file vresutils.shapes.landkreise__.pickle of cache: 4.0 msec
/home/vres/data/tom/lib/load/europe.py:35: FutureWarning: how in .resample() is deprecated
the new syntax is .resample(...).sum()
  load = load.resample('AS', how="sum")
In [41]:
#Kill the Timezone information to avoid pandas bugs
load.index = load.index.values
In [42]:
#Take the first year (in UTC time - we don't set time zone because of a Pandas bug)
network.set_snapshots(pd.date_range("2011-01-01 00:00","2011-12-31 23:00",freq="H"))

print(network.snapshots)
DatetimeIndex(['2011-01-01 00:00:00', '2011-01-01 01:00:00',
               '2011-01-01 02:00:00', '2011-01-01 03:00:00',
               '2011-01-01 04:00:00', '2011-01-01 05:00:00',
               '2011-01-01 06:00:00', '2011-01-01 07:00:00',
               '2011-01-01 08:00:00', '2011-01-01 09:00:00',
               ...
               '2011-12-31 14:00:00', '2011-12-31 15:00:00',
               '2011-12-31 16:00:00', '2011-12-31 17:00:00',
               '2011-12-31 18:00:00', '2011-12-31 19:00:00',
               '2011-12-31 20:00:00', '2011-12-31 21:00:00',
               '2011-12-31 22:00:00', '2011-12-31 23:00:00'],
              dtype='datetime64[ns]', length=8760, freq='H')
In [43]:
#temporary load scaling factor for Germany load in relation to ENTSO-E hourly load
#based roughly on Schumacher & Hirth (2015)
#http://www.feem.it/userfiles/attach/20151191122284NDL2015-088.pdf
#In principle rescaling should happen on a monthly basis

load_factor = 1.12

for bus in graph.nodes():
    network.add("Load",bus,bus=bus,
                p_set = pd.Series(data=load_factor*1000*load.loc[network.snapshots,bus],index=network.snapshots))
In [44]:
%matplotlib inline
In [45]:
pd.DataFrame(load.sum(axis=1)).plot()
Out[45]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fda80f8f810>
In [46]:
load_distribution = network.loads_t.p_set.loc[network.snapshots[0]].groupby(network.loads.bus).sum()
In [47]:
network.plot(bus_sizes=load_distribution)
Out[47]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fda808ac510>
In [48]:
total_load = load.sum(axis=1)
In [49]:
monthly_load = total_load.resample("M").sum()
In [50]:
monthly_load.plot(grid=True)
Out[50]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fda7ba50090>

Attach conventional generators from BNetzA list

In [51]:
from vresutils import shapes as vshapes

def read_kraftwerksliste(with_latlon=True):                                                                              
                                                                                                              
    kraftwerke = pd.read_csv('../../lib/vresutils/data/Kraftwerksliste_CSV_deCP850ed.csv',                                         
                             delimiter=';', encoding='utf-8', thousands='.', decimal=',')                                
    def sanitize_names(x):                                                                                               
        try:                                                                                                             
            x = x[:x.index('(')]                                                                                         
        except ValueError:                                                                                               
            pass                                                                                                         
        return x.replace(u'\n', u' ').strip()
    kraftwerke.columns = kraftwerke.columns.map(sanitize_names)
    
    def sanitize_plz(x):
        try:
            x = x.strip()
            if len(x) > 5:
                x = x[:5]
            return float(x)
        except (ValueError, AttributeError):
            return np.NAN
    kraftwerke.PLZ = kraftwerke.PLZ.apply(sanitize_plz)
    if with_latlon:
        postcodes = {pc: sh.centroid
                     for pc, sh in iteritems(vshapes.postcodeareas())
                     if sh is not None}
        kraftwerke['lon'] = kraftwerke.PLZ.map({pc: c.x for pc, c in iteritems(postcodes)})
        kraftwerke['lat'] = kraftwerke.PLZ.map({pc: c.y for pc, c in iteritems(postcodes)})
        #kraftwerke.dropna(subset=('lon','lat'), inplace=True)                                                           

    kraftwerke[u'Type'] = kraftwerke[u"Auswertung Energieträger"].map({
        u'Erdgas': u'Gas',
        u'Grubengas': u'Gas',
        u'Laufwasser': u'Run of River',
        u'Pumpspeicher': u'Pumped Hydro',
        u'Speicherwasser (ohne Pumpspeicher)': u'Storage Hydro',
        u'Mineralölprodukte': u'Oil',
        u'Steinkohle': u'Hard Coal',
        u'Braunkohle': u'Brown Coal',
        u'Abfall': u'Waste',
        u'Kernenergie': u'Nuclear',
        u'Sonstige Energieträger\n(nicht erneuerbar) ': u'Other',
        u'Mehrere Energieträger\n(nicht erneuerbar)': u'Multiple',
        u'Biomasse' : u'Biomass',
        u'Deponiegas' : u'Gas',
        u'Klärgas' : u'Gas',
        u'Geothermie' : u'Geothermal',
        u'Windenergie (Onshore-Anlage)' : u'Wind Onshore',
        u'Windenergie (Offshore-Anlage)' : u'Wind Offshore',
        u'Solare Strahlungsenergie' : u'Solar',
        u'Unbekannter Energieträger\n(nicht erneuerbar)' : u'Other'
    })

    return kraftwerke
In [52]:
power_plants = read_kraftwerksliste()
Serving call to postcodeareas from file vresutils.shapes.postcodeareas_ver2__.pickle of cache: 85.7 msec
In [53]:
power_plants[power_plants[u"Unternehmen"] == "EEG-Anlagen < 10 MW"].groupby(u"Type").sum()
Out[53]:
PLZ Netto-Nennleistung lon lat
Type
Biomass NaN 5444.6 NaN NaN
Gas NaN 497.0 NaN NaN
Geothermal NaN 31.7 NaN NaN
Run of River NaN 1438.4 NaN NaN
Solar NaN 37008.7 NaN NaN
Wind Onshore NaN 24394.9 NaN NaN
In [54]:
power_plants.groupby(u"Type").sum()
Out[54]:
PLZ Netto-Nennleistung lon lat
Type
Biomass 3016585.0 6678.72 639.913564 3213.997738
Brown Coal 2485734.0 22969.50 655.576325 3787.688000
Gas 14304542.0 29356.03 2579.020814 13424.719920
Geothermal NaN 31.70 NaN NaN
Hard Coal 5979186.0 31474.60 868.212204 5166.446722
Multiple 45772.0 152.70 7.118728 51.686962
Nuclear 618554.0 12068.00 77.599904 400.252062
Oil 2616997.0 4165.20 563.075179 2803.594992
Other 1516594.0 3037.80 303.710411 1750.250968
Pumped Hydro 1991089.0 9239.50 474.335182 2210.915360
Run of River 6890713.0 3999.10 899.795675 4087.882170
Solar 3484421.0 39190.00 1730.533616 7279.165506
Storage Hydro 262091.0 1445.00 39.067293 197.592367
Waste 4111475.0 1703.20 756.186790 4066.710387
Wind Offshore NaN 1587.40 NaN NaN
Wind Onshore 22458563.0 37315.30 6484.788929 33525.165765
In [55]:
import random

#NB: bnetza extracted from BNetzA using

#./Kraftwerksdaten.ipynb


def backup_capacity_german_grid(G):   

    from shapely.geometry import Point

    plants = power_plants
    plants = plants[plants["Kraftwerksstatus"] == u"in Betrieb"]
    
    #remove EEG-receiving power plants - except biomass, these will be added later
    
    #it's necessary to remove biomass because we don't have coordinates for it
    
    for tech in ["Solar","Wind Onshore","Wind Offshore","Biomass"]:
        plants = plants[plants['Type'] != tech]
    
    cells = {n: d["region"]
             for n, d in G.nodes_iter(data=True)}

    def nodeofaplant(x):
        if np.isnan(x["lon"]) or np.isnan(x["lat"]):
            return random.choice(list(cells.keys()))
        p = Point(x["lon"], x["lat"])
        for n, cell in iteritems(cells):
            if cell.contains(p):
                return n
        else:
            return min(cells, key=lambda n: cells[n].distance(p))
    nodes = plants.apply(nodeofaplant, axis=1)

    capacity = plants['Netto-Nennleistung'].groupby((nodes, plants[u'Type'])).sum() / 1e3
    capacity.name = 'Capacity'

    return capacity
In [56]:
cap = backup_capacity_german_grid(graph)
In [57]:
cap.describe(),cap.sum(),type(cap)
Out[57]:
(count    477.000000
 mean       0.218376
 std        0.442888
 min        0.000100
 25%        0.023100
 50%        0.062700
 75%        0.193700
 max        4.730000
 Name: Capacity, dtype: float64, 104.16513, pandas.core.series.Series)
In [58]:
print(cap[pd.isnull(cap)])
     Type        
141  Pumped Hydro   NaN
275  Other          NaN
Name: Capacity, dtype: float64
In [59]:
cap.fillna(0.1,inplace=True)
In [60]:
cap.index.levels[1]
Out[60]:
Index([u'Brown Coal', u'Gas', u'Geothermal', u'Hard Coal', u'Multiple',
       u'Nuclear', u'Oil', u'Other', u'Pumped Hydro', u'Run of River',
       u'Storage Hydro', u'Waste'],
      dtype='object', name=u'Type')
In [61]:
m_costs = {"Gas" : 50.,
           "Brown Coal" : 10.,
           "Hard Coal" : 25.,
           "Oil" : 100.,
           "Nuclear" : 8.,
           "Pumped Hydro" : 3.,
           "Storage Hydro" : 3.,
           "Run of River" : 3.,
           "Geothermal" : 26.,
           "Waste" : 6.,
           "Multiple" : 28.,
           "Other" : 32.}

default_cost = 10.
In [62]:
for (bus_name,tech_name) in cap.index:
    print(bus_name,tech_name,cap[(bus_name,tech_name)])
    if tech_name == "Pumped Hydro":
        network.add("StorageUnit",bus_name + " " + tech_name,
                bus=bus_name,p_nom=1000*cap[(bus_name,tech_name)],
                marginal_cost=m_costs.get(tech_name,default_cost),
                carrier=tech_name,
                max_hours = 6,
                efficiency_store=0.95,
                efficiency_dispatch=0.95)
    else:
        network.add("Generator",bus_name + " " + tech_name,
                bus=bus_name,p_nom=1000*cap[(bus_name,tech_name)],
                marginal_cost=m_costs.get(tech_name,default_cost),
                carrier=tech_name)   
1 Gas 0.121
1 Hard Coal 0.272
100_220kV Pumped Hydro 0.1445
102 Gas 0.0679
108 Run of River 0.0631
108 Waste 0.038
111 Gas 0.0135
112 Gas 0.0135
112 Run of River 0.0111
114 Hard Coal 0.019
114 Pumped Hydro 0.138
115 Brown Coal 4.73
116 Gas 0.5433
118 Gas 0.095
118 Waste 0.012
119 Brown Coal 0.0093
119 Gas 0.0169
119 Hard Coal 0.0615
121 Brown Coal 0.0246
121 Pumped Hydro 0.238
123_220kV Brown Coal 3.455
123_220kV Gas 0.0965
126 Gas 0.002
126 Oil 0.06
127 Brown Coal 1.8
127 Waste 0.024
128_220kV Run of River 0.0364
129_220kV Hard Coal 0.0615
12_220kV Gas 0.15
12_220kV Waste 0.03
130 Hard Coal 0.085
130 Waste 0.03
131 Gas 0.5863
132 Storage Hydro 0.029
133 Brown Coal 0.137
133 Waste 0.0334
134 Brown Coal 0.118
136 Gas 0.838
136 Hard Coal 0.765
137 Gas 0.0114
139 Hard Coal 0.76
14 Run of River 0.027
140 Gas 0.0111
140 Pumped Hydro 0.1
141 Gas 0.0647
141 Pumped Hydro 0.1
142 Gas 0.0098
142 Hard Coal 0.449
145 Brown Coal 0.0753
145 Gas 0.5365
145 Other 0.0798
145 Waste 0.0451
146 Gas 0.084
147 Gas 0.23
148 Gas 0.0102
148 Oil 0.0862
148 Waste 0.0537
149 Run of River 0.1287
150 Gas 0.0158
153 Waste 0.012
155 Waste 0.0126
156_220kV Brown Coal 0.0289
157_220kV Brown Coal 0.857
157_220kV Gas 0.022
158 Pumped Hydro 0.1
16 Gas 0.0114
160_220kV Hard Coal 0.103
160_220kV Waste 0.014
161 Gas 0.05858
166 Gas 0.3336
166 Hard Coal 1.3245
166 Pumped Hydro 0.0015
167 Gas 0.026
167 Run of River 0.0469
169 Gas 0.0337
170 Gas 0.0231
170 Waste 0.0615
174 Gas 0.0245
174 Other 0.038
175 Hard Coal 0.095
175 Other 0.021
176 Gas 0.0033
177 Multiple 0.0327
178 Gas 0.0092
178 Storage Hydro 0.015
180 Gas 0.0254
180 Oil 0.3335
180 Waste 0.0289
183 Gas 0.0225
183 Hard Coal 0.766
183 Oil 0.038
183 Waste 0.024
184_220kV Gas 0.0223
185 Waste 0.0331
188 Gas 0.274
18_220kV Gas 0.0528
19 Gas 0.0106
19 Hard Coal 0.37
190 Gas 0.0134
194 Hard Coal 1.017
196 Hard Coal 0.695
196 Other 0.7
197 Run of River 0.0401
198 Gas 0.015
200 Brown Coal 0.0275
200 Gas 0.0074
201_220kV Gas 0.0581
202 Gas 1.8701
202 Nuclear 1.329
203 Gas 0.0205
204 Hard Coal 0.31
205 Gas 0.417
205 Pumped Hydro 0.153
209 Hard Coal 0.026
21 Gas 0.127
21 Hard Coal 0.194
21 Waste 0.0164
210 Waste 0.0147
211 Hard Coal 0.721
211 Waste 0.0116
212 Run of River 0.024
213 Gas 0.1046
214 Gas 0.027
215 Gas 0.0472
215 Waste 0.0734
217 Gas 0.108
217 Hard Coal 0.508
217 Waste 0.017
219 Gas 0.2545
219 Waste 0.034
221 Geothermal 0.0078
222 Gas 0.099
223 Gas 0.087
224 Other 0.1761
225 Gas 0.0116
225 Hard Coal 0.0134
225 Waste 0.015
226 Gas 0.02
226 Run of River 0.0001
227 Run of River 0.011
228 Pumped Hydro 0.1
229 Gas 0.0245
229 Hard Coal 1.958
22_220kV Oil 0.2985
233_220kV Gas 0.004
233_220kV Hard Coal 0.0207
235 Other 0.0197
235 Pumped Hydro 0.47
238 Run of River 0.0121
239 Gas 0.8
240 Gas 0.126
245 Hard Coal 0.0615
246 Hard Coal 0.51
247_220kV Oil 0.0345
247_220kV Run of River 0.6332
248_220kV Hard Coal 0.895
249 Waste 0.0146
25 Pumped Hydro 0.196
250_220kV Other 0.3966
254 Gas 0.0107
255 Run of River 0.0848
258 Gas 0.032
258 Other 0.028
260 Hard Coal 0.1385
266 Gas 0.0599
266 Oil 0.0928
266 Pumped Hydro 0.0988
266 Run of River 0.214
266 Waste 0.009
268 Pumped Hydro 0.1
26_220kV Nuclear 1.41
270 Gas 0.05
275 Gas 0.042
275 Other 0.1
275 Waste 0.0125
276 Hard Coal 0.11
276 Run of River 0.016
279_220kV Brown Coal 0.0908
279_220kV Pumped Hydro 1.0452
28 Hard Coal 0.26
28 Oil 0.101
28 Run of River 0.3118
281_220kV Gas 0.074
282_220kV Oil 0.017
285 Gas 0.0181
287 Gas 0.064
288 Hard Coal 0.0185
289 Gas 0.0334
289 Hard Coal 0.794
290 Waste 0.0157
291 Brown Coal 1.105
292_220kV Brown Coal 0.0203
292_220kV Gas 0.216
292_220kV Waste 0.059
293 Pumped Hydro 0.3776
294 Gas 0.0135
294 Waste 0.0112
295 Gas 0.119
295 Hard Coal 0.0433
298 Gas 0.0375
299 Hard Coal 0.69
3 Run of River 0.0188
301 Gas 0.013
302 Brown Coal 0.0335
302 Gas 0.128
302 Other 0.056
302 Storage Hydro 0.079
302 Waste 0.0147
304 Gas 0.0342
304 Oil 0.0248
305 Nuclear 1.31
305 Oil 0.136
308_220kV Gas 0.0274
308_220kV Nuclear 1.275
30_220kV Run of River 0.003
31 Gas 0.1954
310 Hard Coal 0.778
313 Brown Coal 0.875
315_220kV Run of River 0.0774
316 Gas 0.0758
316 Waste 0.01
318 Run of River 0.0627
320 Gas 0.0098
320 Pumped Hydro 0.09
321 Gas 0.0002
321 Waste 0.0136
322_220kV Gas 0.0851
323 Gas 0.0129
326 Gas 0.3198
327 Gas 0.0269
328 Run of River 0.0365
329 Run of River 0.05
329 Storage Hydro 0.5
32_220kV Gas 0.1091
32_220kV Pumped Hydro 0.22
33 Gas 0.0318
330 Nuclear 1.41
330 Run of River 0.0358
331 Gas 0.227
331 Run of River 0.2175
331 Waste 0.0125
332_220kV Gas 0.0979
332_220kV Run of River 0.0252
333 Pumped Hydro 0.525
335 Gas 0.054
336_220kV Gas 0.0528
336_220kV Run of River 0.089
337 Gas 0.0134
337 Run of River 0.0801
338_220kV Waste 0.054
34 Hard Coal 0.757
34 Oil 0.056
340 Gas 0.2996
340 Storage Hydro 0.0825
341 Brown Coal 0.134
341 Gas 0.159
341 Oil 0.111
341 Other 0.0427
341 Waste 0.058
342 Gas 0.0195
342 Pumped Hydro 0.0797
343 Pumped Hydro 1.052
344 Gas 0.06995
347 Gas 0.211
347 Hard Coal 0.089
347 Oil 0.327
348 Gas 0.4589
350 Gas 0.0186
351_220kV Other 0.0233
351_220kV Run of River 0.0637
351_220kV Waste 0.0255
352 Hard Coal 0.472
352 Oil 0.046
352 Run of River 0.0733
353_220kV Gas 1.391
353_220kV Run of River 0.0233
354 Gas 0.5405
354 Hard Coal 0.333
354 Waste 0.04
358 Hard Coal 0.318
36 Gas 0.01106
36 Hard Coal 0.28
360_220kV Gas 0.042
360_220kV Hard Coal 0.779
360_220kV Waste 0.016
362_220kV Pumped Hydro 0.1
362_220kV Run of River 0.0795
363 Other 0.016
364 Waste 0.0116
365 Gas 0.0818
365 Hard Coal 0.688
365 Waste 0.036
367 Nuclear 1.402
368 Pumped Hydro 0.91
370 Brown Coal 0.0568
370 Other 0.085
373 Run of River 0.025
374 Gas 0.1002
374 Hard Coal 0.725
374 Pumped Hydro 0.1
376 Geothermal 0.0006
376 Hard Coal 0.875
377 Waste 0.011
378 Gas 0.26044
378 Hard Coal 0.769
379 Pumped Hydro 0.0398
37_220kV Gas 0.052
37_220kV Nuclear 1.36
37_220kV Waste 0.0584
380 Gas 0.25
380 Waste 0.0225
382_220kV Gas 0.122
382_220kV Waste 0.0207
383 Waste 0.0244
385_220kV Gas 0.0584
385_220kV Hard Coal 0.0248
385_220kV Oil 0.0005
386 Gas 0.023
389 Pumped Hydro 0.164
39 Gas 2.2993
390 Gas 1.0013
390 Hard Coal 0.0174
390 Waste 0.018
394 Gas 0.0265
395 Gas 0.49
396 Nuclear 2.572
397 Pumped Hydro 0.1445
398 Oil 0.087
399 Gas 0.0351
399 Waste 0.0095
40 Run of River 0.048
401 Gas 0.031
402 Oil 0.088
402 Run of River 0.06
402 Storage Hydro 0.158
403_220kV Run of River 0.019
404_220kV Run of River 0.024
405 Gas 0.056
405 Run of River 0.0846
409 Gas 0.0085
409 Hard Coal 0.1126
409 Oil 0.0699
409 Waste 0.0195
410 Run of River 0.0201
414 Gas 0.1937
417 Run of River 0.1373
418 Run of River 0.066
421_220kV Waste 0.0673
422 Gas 0.0135
423 Pumped Hydro 0.36
423 Run of River 0.0368
424 Oil 0.07
424 Pumped Hydro 0.295
425 Other 0.2885
426 Gas 0.0324
428 Brown Coal 0.074
429 Gas 0.1121
429 Hard Coal 0.012
43 Gas 0.196
43 Waste 0.017
430 Brown Coal 0.0385
433 Gas 0.0156
434 Run of River 0.0128
434 Storage Hydro 0.124
435 Brown Coal 0.164
438 Run of River 0.146
439 Gas 0.0117
44 Waste 0.026
440 Other 0.017
441 Gas 0.182
441 Hard Coal 0.0615
442 Gas 0.0265
444 Gas 0.43
445 Run of River 0.2171
446 Gas 0.0147
447 Gas 0.075
45 Gas 0.296
45 Geothermal 0.0231
451 Run of River 0.1044
453 Run of River 0.03
455 Gas 0.0169
456 Run of River 0.0001
457 Waste 0.0296
463 Gas 0.0258
464 Gas 0.0519
464 Oil 0.0663
464 Run of River 0.0254
465 Gas 0.035
465 Oil 0.08
468 Gas 0.0175
468 Hard Coal 0.079
468 Oil 0.0048
468 Waste 0.016
47 Gas 0.444
47 Run of River 0.0317
474 Gas 0.39
475 Oil 0.2506
478 Gas 0.095
479 Pumped Hydro 0.1
481 Pumped Hydro 0.1
484 Gas 0.065
484 Other 0.006
484 Run of River 0.0667
484 Storage Hydro 0.0455
485 Run of River 0.0265
486 Hard Coal 0.1385
486 Oil 0.0114
487 Brown Coal 0.875
487 Oil 0.0114
49 Brown Coal 0.352
49 Hard Coal 0.066
49 Waste 0.0375
492 Hard Coal 0.054
492 Waste 0.0569
493 Pumped Hydro 0.1
495 Hard Coal 0.147
51 Pumped Hydro 0.623
51 Storage Hydro 0.02
52_220kV Gas 0.39
52_220kV Hard Coal 1.347
56_220kV Gas 0.0117
58 Gas 0.0525
58 Hard Coal 0.0269
59 Hard Coal 0.324
5_220kV Run of River 0.0046
61_220kV Gas 0.334
61_220kV Hard Coal 0.35
62 Hard Coal 0.119
62 Run of River 0.0099
62 Waste 0.044
63 Brown Coal 0.465
63 Pumped Hydro 0.043
63 Run of River 0.022
64 Run of River 0.0112
65 Run of River 0.1645
66 Pumped Hydro 0.127
66 Run of River 0.0051
67 Gas 0.0148
67 Hard Coal 0.3
67 Oil 0.088
67 Other 0.15
67 Storage Hydro 0.392
67 Waste 0.033
68_220kV Run of River 0.0042
69_220kV Gas 0.0115
69_220kV Waste 0.0188
70 Gas 0.423
70 Waste 0.0156
71 Brown Coal 0.049
71 Gas 0.1611
71 Waste 0.0099
72 Brown Coal 0.9
72 Gas 0.151
72 Oil 0.12
72 Pumped Hydro 0.198
74_220kV Gas 0.122
74_220kV Waste 0.0233
75 Brown Coal 0.049
77 Gas 0.0053
78_220kV Pumped Hydro 0.0798
79 Brown Coal 0.75
79 Waste 0.0187
8 Pumped Hydro 0.1191
81 Brown Coal 2.835
81 Other 0.1245
83 Brown Coal 0.75
83 Other 0.0536
84 Hard Coal 0.323
87_220kV Gas 0.1387
87_220kV Hard Coal 0.1759
87_220kV Multiple 0.12
88 Gas 0.0398
89 Other 0.606
92 Gas 0.066
94_220kV Geothermal 0.0002
95_220kV Gas 0.0243
98 Pumped Hydro 0.247
99_220kV Gas 0.1136
99_220kV Hard Coal 0.6556

Add renewables

In [63]:
import generation.germany as DEgen

reload(DEgen)

generation = DEgen.timeseries_eeg(graph)
Caching call to timeseries_eeg in generation.germany.timeseries_eeg_scigrid_v2_.pickle: 
.. Get eeg onshore wind timeseries from REatlas: 
.. .. Caching call to eeg_windonlayouts_per_class in generation.germany.eeg_windonlayouts_per_class_scigrid_v2_Cutout_beckerEurope_2011_2014_.pickle: 
.. .. .. Caching call to capacity_layouts_wind in generation.germany.capacity_layouts_wind__cutout.Cutout_beckerEurope_2011_2014_normed.False_turbines.e325f1f75c97f4e5a4146db18354d997c44e9185.pickle: 
.. .. .. .. Serving call to postcodes2gridcells from file generation.germany.postcodes2gridcells_Cutout_beckerEurope_2011_2014_.pickle of cache: 121.2 msec
.. .. .. .. 4.5 sec
.. .. .. 317.9 sec
.. .. 890.8 sec
.. Get eeg offshore wind timeseries from REatlas: 
.. .. Caching call to eeg_windofflayouts_per_class in generation.germany.eeg_windofflayouts_per_class_scigrid_v2_Cutout_beckerEurope_2011_2014_.pickle: 49.0 msec
.. .. 424.1 sec
.. Get eeg solar timeseries from REatlas: 
.. .. Caching call to eeg_solarlayouts in generation.germany.eeg_solarlayouts_scigrid_v2_Cutout_beckerEurope_2011_2014_.pickle: 
.. .. .. Caching call to capacity_layout_solar in generation.germany.capacity_layout_solar__cutout.Cutout_beckerEurope_2011_2014_normed.False.pickle: 
.. .. .. .. Serving call to postcodes2gridcells from file generation.germany.postcodes2gridcells_Cutout_beckerEurope_2011_2014_.pickle of cache: 49.5 msec
.. .. .. .. 2.4 sec
.. .. .. 42.1 sec
.. .. 118.8 sec
.. 1435.0 sec
In [64]:
generation.items
Out[64]:
Index([u'solar', u'wind', u'windoff', u'windon'], dtype='object')
In [65]:
#Kill the Timezone information to avoid pandas bugs
generation.major_axis = generation.major_axis.values
In [66]:
generation.loc[["wind","solar"],network.snapshots,:].sum(axis=2).plot()
Out[66]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fda7b853c90>
In [67]:
solar = generation.loc["solar",network.snapshots,:].sum(axis=1)
solar.describe()
Out[67]:
count    8760.000000
mean        5.978042
std         8.963907
min         0.000000
25%         0.000000
50%         0.000000
75%        11.014865
max        34.577076
dtype: float64
In [68]:
#make sure the ordering of the minor axis is correc
generation.minor_axis = graph.nodes()

Get the capacities correct

In [69]:
cutout = vresutils.reatlas.Cutout(cutoutname="Europe_2011_2014", username="becker")
In [70]:
def panel_capacity(panel):
    """
    Returns the panel capacity in MW.
    
    Parameters
    ----------
    panel : string
        Panel name, e.g. "Sunpower"
    
    Returns
    -------
    capacity : float
        In MW
    """
    c = vresutils.reatlas.solarpanelconf_to_solar_panel_config_object(panel)
    return c['A'] + c['B'] * 1000 + c['C'] * np.log(1000)
In [71]:
solar_layouts = DEgen.eeg_solarlayouts(graph,cutout)
Serving call to eeg_solarlayouts from file generation.germany.eeg_solarlayouts_scigrid_v2_Cutout_beckerEurope_2011_2014_.pickle of cache: 137.9 msec
In [72]:
panel_cap = panel_capacity(solar_layouts[0]["panel"])
solar_caps = pd.Series(solar_layouts[1].sum(axis=(1,2))*panel_cap,
                       graph.nodes())
In [73]:
solar_caps.describe(),solar_caps.sum()
Out[73]:
(count    489.000000
 mean       0.075750
 std        0.093728
 min        0.000028
 25%        0.012453
 50%        0.042313
 75%        0.100619
 max        0.849577
 dtype: float64, 37.041524779087617)
In [74]:
(generation.solar.max()/solar_caps).describe()
Out[74]:
count    489.000000
mean       0.983755
std        0.018969
min        0.948734
25%        0.970006
50%        0.984081
75%        0.996938
max        1.043963
dtype: float64
In [75]:
windon_layouts = DEgen.eeg_windonlayouts_per_class(graph,cutout)
Serving call to eeg_windonlayouts_per_class from file generation.germany.eeg_windonlayouts_per_class_scigrid_v2_Cutout_beckerEurope_2011_2014_.pickle of cache: 842.1 msec
In [76]:
windon_capacities = pd.DataFrame(index=graph.nodes())
for turbine_items in windon_layouts:
    name = turbine_items[0]["onshore"]
    turbine_cap = np.array(vresutils.reatlas.turbineconf_to_powercurve_object(name)["POW"]).max()
    print(name,turbine_cap)
    windon_capacities[name] = turbine_items[1].sum(axis=(1,2))*turbine_cap/1000.
Vestas_V25_200kW 0.2001
Vestas_V47_660kW 0.66
Bonus_B1000_1000kW 1.0
Suzlon_S82_1.5_MW 1.5
Vestas_V66_1750kW 1.75
Vestas_V80_2MW_gridstreamer 2.0
Siemens_SWT_2300kW 2.3
Vestas_V90_3MW 3.0
In [77]:
windon_caps = windon_capacities.sum(axis=1)
windon_caps.describe(),windon_caps.sum()
Out[77]:
(count    489.000000
 mean       0.076360
 std        0.146817
 min        0.000000
 25%        0.002402
 50%        0.020759
 75%        0.085740
 max        1.297764
 dtype: float64, 37.339895329191734)
In [78]:
(generation.windon.max()/windon_caps).describe()
Out[78]:
count    488.000000
mean       0.999618
std        0.004390
min        0.909028
25%        0.999927
50%        0.999961
75%        0.999985
max        1.000000
dtype: float64
In [79]:
windoff_layouts = DEgen.eeg_windofflayouts_per_class(graph,cutout)
Serving call to eeg_windofflayouts_per_class from file generation.germany.eeg_windofflayouts_per_class_scigrid_v2_Cutout_beckerEurope_2011_2014_.pickle of cache: 585.6 msec
In [80]:
windoff_capacities = pd.DataFrame(index=graph.nodes())
for i,turbine_items in enumerate(windoff_layouts):
    name = turbine_items[0]["offshore"]
    turbine_cap = np.array(vresutils.reatlas.turbineconf_to_powercurve_object(name)["POW"]).max()
    print(name,turbine_cap)
    #add an index to name to avoid duplication of names
    windoff_capacities[name+"-" + str(i)] = turbine_items[1].sum(axis=(1,2))*turbine_cap/1000.
Siemens_SWT_2300kW_hornsrev2 2.3
NREL_ReferenceTurbine_5MW_offshore 5.0
NREL_ReferenceTurbine_5MW_offshore 5.0
Siemens_SWT_107_3600kW 3.6
Siemens_SWT_107_3600kW 3.6
NREL_ReferenceTurbine_5MW_offshore 5.0
In [81]:
windoff_capacities.sum()
Out[81]:
Siemens_SWT_2300kW_hornsrev2-0          0.0483
NREL_ReferenceTurbine_5MW_offshore-1    0.0600
NREL_ReferenceTurbine_5MW_offshore-2    0.4000
Siemens_SWT_107_3600kW-3                1.4472
Siemens_SWT_107_3600kW-4                0.1080
NREL_ReferenceTurbine_5MW_offshore-5    0.9100
dtype: float64
In [82]:
windoff_caps = windoff_capacities.sum(axis=1)
windoff_caps.describe(),windoff_caps.sum()
Out[82]:
(count    489.000000
 mean       0.006081
 std        0.081966
 min        0.000000
 25%        0.000000
 50%        0.000000
 75%        0.000000
 max        1.447200
 dtype: float64, 2.9735)
In [83]:
(generation.windoff.max()/windoff_caps).describe()
Out[83]:
count    5.000000e+00
mean     1.000000e+00
std      1.263743e-08
min      1.000000e+00
25%      1.000000e+00
50%      1.000000e+00
75%      1.000000e+00
max      1.000000e+00
dtype: float64
In [84]:
network.plot(bus_sizes=1000*windoff_caps)
Out[84]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fda7b7e2710>
In [85]:
network.plot(bus_sizes=1000*windon_caps)
Out[85]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fda7a39b690>
In [86]:
network.plot(bus_sizes=1000*solar_caps)
Out[86]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fda79235a50>
In [87]:
d = {"windoff" : {"full_name" : "Wind Offshore", "caps" : windoff_caps},
    "windon" : {"full_name" : "Wind Onshore", "caps" : windon_caps},
    "solar" : {"full_name" : "Solar", "caps" : solar_caps},
     }

for tech in ["windoff",'windon',"solar"]:
    caps = d[tech]["caps"]
    caps = caps[caps != 0]
    
    for i in caps.index:
        network.add("Generator","{} {}".format(i,d[tech]["full_name"]),
                    p_nom=caps[i]*1000.,dispatch="variable",
                    bus=i,carrier=d[tech]["full_name"],
                    p_max_pu=generation[tech].loc[network.snapshots,i]/caps[i])
In [88]:
csv_folder_name = "../../lib/data/de_model/scigrid-with-load-gen-trafos"

network.export_to_csv_folder(csv_folder_name)
In [89]:
network.set_snapshots(network.snapshots[:24])

csv_folder_name = "../../lib/pypsa/examples/scigrid-de/scigrid-with-load-gen-trafos"


network.export_to_csv_folder(csv_folder_name)
In [ ]: