Introduction to Energy Modelling - A ML Approach

API Calls

The code below uses the new RDP Search service which is an excellent tool for finding content, in our case Physical Assets such as windfarms - but there are many other search views to help you out such as instruments, people, entities and many more (reference guide). In the code below we provide a filter which contains the query syntax we need to narrow down the Physical Asset universe to the ones we want - wind farms located in Germany.

3519 rows × 5 columns

Next we provide the list of physical asset RICs from the previous search to our get_data API call, along with a few fields we are interested in.

3519 rows × 6 columns

Visualization Workflow

We will be creating some visualisations using geopandas and shape files - to do this we need to create a geometry point type (a coordinate of longitude and latitude) so we can work with polygonal maps. Here we zip our longitude and latitude together to form this new Point datatype. We store this in a new geometry column.

3519 rows × 8 columns

Open source wind data

Next we wish to download some wind data. There is plenty of open source data available - the data used here is sourced from Deutscher Wetterdienst (see further resources section for more detail) - this is from various weather monitoring stations in Germany - which record wind speed at a height of 10 meters. Maybe you have store of such or similar information in your own company. Here we are just reading it into a dataframe, wind. Data basis: Deutscher Wetterdienst, aggregated.

5 rows × 52 columns

We also want to get total wind output data for Germany - this I got from our Point Connect data service.

Now we add two additional features Total Wind Speed (measuring across all weather stations) and Mean Wind Speed. Data basis: Deutscher Wetterdienst, average and total added.

5 rows × 53 columns

We now join our wind output dataframe to our main wind features dataframe using an inner join on Date.

8 rows × 54 columns

Locating Weather Monitoring Stations

For more detailed modelling I have included this section of code, where we can download the geographic locations of each of our wind monitoring sites. This will allow us to calculate the closest wind monitoring station to each of the windfarms and perhaps get a more accurate indication of local weather conditions on output. First we load the locations of our wind monitoring sites and then create a new geometry type for them. Data basis: Deutscher Wetterdienst, gridded data reproduced graphically

Calculate nearest weather station for each windfarm

Here we will calculate the closest weather monitoring station for each of our windfarms using thee sciipy.spatial.distance cdist routine which calculates the distance from one point to a range of coordinates - then we just select the minimum (closest) one and store it to our dataframe.

Visualizations using GeoPandas, Shapefiles and Matplotlib

We will now do some visualisations working with GeoPandas and shapefiles. GeoPandas extends the datatypes used by pandas to allow spatial operations on geometric types. Shapefiles are geospatial vector data formats that allow for storing of geographic type features such as lines, points and polygons (areas). There are many shapefiles available on the internet. I am using one I found for Germany and have included this is the the github repo containing this notebook so you can replicate. Here we just read in the file (I have also printed it so you can take a look at it) and then plot it using matplotlib.

GID_0 NAME_0 GID_1 NAME_1 GID_2 \
0 DEU Germany DEU.1_1 Baden-Württemberg DEU.1.1_1
1 DEU Germany DEU.1_1 Baden-Württemberg DEU.1.1_1
2 DEU Germany DEU.1_1 Baden-Württemberg DEU.1.1_1
3 DEU Germany DEU.1_1 Baden-Württemberg DEU.1.1_1
4 DEU Germany DEU.1_1 Baden-Württemberg DEU.1.1_1
... ... ... ... ... ...
11297 DEU Germany DEU.16_1 Thüringen DEU.16.22_1
11298 DEU Germany DEU.16_1 Thüringen DEU.16.22_1
11299 DEU Germany DEU.16_1 Thüringen DEU.16.22_1
11300 DEU Germany DEU.16_1 Thüringen DEU.16.22_1
11301 DEU Germany DEU.16_1 Thüringen DEU.16.22_1

NAME_2 GID_3 NAME_3 GID_4 \
0 Alb-Donau-Kreis DEU.1.1.1_1 Allmendingen DEU.1.1.1.1_1
1 Alb-Donau-Kreis DEU.1.1.1_1 Allmendingen DEU.1.1.1.2_1
2 Alb-Donau-Kreis DEU.1.1.2_1 Blaubeuren DEU.1.1.2.1_1
3 Alb-Donau-Kreis DEU.1.1.2_1 Blaubeuren DEU.1.1.2.2_1
4 Alb-Donau-Kreis DEU.1.1.3_1 Blaustein DEU.1.1.3.1_1
... ... ... ... ...
11297 Weimarer Land DEU.16.22.9_1 Nordkreis Weimar DEU.16.22.9.12_1
11298 Weimarer Land DEU.16.22.9_1 Nordkreis Weimar DEU.16.22.9.13_1
11299 Weimarer Land DEU.16.22.9_1 Nordkreis Weimar DEU.16.22.9.14_1
11300 Weimarer Land DEU.16.22.9_1 Nordkreis Weimar DEU.16.22.9.15_1
11301 Weimarer Land DEU.16.22.9_1 Nordkreis Weimar DEU.16.22.9.16_1

NAME_4 VARNAME_4 TYPE_4 ENGTYPE_4 CC_4 \
0 Allmendingen None Gemeinde Town 084255001002
1 Altheim None Gemeinde Town 084255001004
2 Berghülen None Gemeinde Town 084255002017
3 Blaubeuren None Stadt Town 084255002020
4 Blaustein None Stadt Town 084250141141
... ... ... ... ... ...
11297 Rohrbach None Gemeinde Town 160715013081
11298 Sachsenhausen None Gemeinde Town 160715013082
11299 Schwerstedt None Gemeinde Town 160715013085
11300 Vippachedelhausen None Gemeinde Town 160715013092
11301 Wohlsborn None Gemeinde Town 160715013097

geometry
0 (POLYGON ((9.747226715087947 48.32448959350592...
1 (POLYGON ((9.76099872589117 48.32106399536138,...
2 POLYGON ((9.82752799987793 48.45702743530279, ...
3 POLYGON ((9.775817871093693 48.36040115356457,...
4 POLYGON ((9.82752799987793 48.45702743530279, ...
... ...
11297 POLYGON ((11.37983036041265 51.07011032104498,...
11298 POLYGON ((11.39456176757824 51.03350448608404,...
11299 POLYGON ((11.30850887298595 51.05867004394531,...
11300 POLYGON ((11.19609355926508 51.10009002685547,...
11301 POLYGON ((11.37327384948736 51.02006912231445,...

[11302 rows x 15 columns]

Next we just convert our df dataframe containing all our windfarm geometry points and our location dataframe that contains the location of our weather monitoring stations to GeoDataFrames and plot both the map, and overlay our windfarms as green dots and our weather monitoring stations as blue dots. Data basis: Deutscher Wetterdienst, gridded data reproduced graphically.

Get temperature data

Here we are using hourly air temperature. This could be a demand side feature as colder temperatures might lead to more heating use or warmer temperatures might lead to more airconditioning etc. Note we are only taking the data from one station in middle of Germany, but for our exploratory purposes this will give us some indication of temperature conditions. Data basis: Deutscher Wetterdienst. More information is available in the Further Resources section.

5 rows × 55 columns

Data on Solar energy

Here we are using hourly sum of solar incoming radiation data. At night the sun does not shine - so we should expect zero readings and from sunrise it should start to increase till sunset. Note, this is for only one area in the centre of Germany but for our exploratory purposes it will at least give us some indication on sunny conditions. Data basis: Deutscher Wetterdienst More information is available in the Further Resources section.

5 rows × 56 columns

Getting Prices Data (our Labels)

Thus far we have been building our feature matrix (our wind, temperature and solar features) and now we wish to build our label dataset - in our case we wish to look at hourly EPEX Spot energy prices. All price data on Refinitiv Platforms is in GMT - whereas all the ther data is in local (CET) - therefore we need to shift our price data forward by 1hr to align it correctly. The source is our Point Connect repository - again have included this file in the github repo so you can replicate.

Sometimes ML routines can work better with smoothed data so I have included a moving average, I don't use it here but include it for reference.

Creating z-scores

We also wish to create some z-scores to normalise variation across the dataset.

41617 rows × 4 columns

We also need to normalize our feature matrix before handing it to our ML routine.

30977 rows × 116 columns

Temporal Considerations

We want to compare the same hourly data for each day - so we need to create a daily series for each hour in the day as it is these series we want to model. The reason is that demand/supply profiles are different each hour of the day so - its not correct to compare hours in a continuous series ie hour 01 to hour 02 etc. The correct set up should be day 1 hour 1, day 2 hour 1, day 3 hour 1 etc. In this way we build up 24 discrete hourly profiles - as it is these we wish to model. Below we just extract the hour of our timestamp.

Implement Machine Learning Model

For our modelling we will use an xgboost model - which a gradient-boosted decision tree class of algorithm - which is quite popular in the literature. We will apply the simplest vanilla implementation here to see if we can get any usable results. Obviously the ML area is quite technical and requires some mastery but you can use this codebase to extend and experiement as appropriate.

Split our dataset in to train, evaluation and test sets

We need to split our dataset into a training set, an evaluation set and a test set. We will hand the the ML routine only the training and evaluation datasets to generate its models. We will then test predictions made by that model against the test dataset to see if we can find anything useful. We split the data roughly at the halfway mark using 10% of our data as evaluation set.

array(['00', '01', '02', '03', '04', '05', '06', '07', '08', '09', '10',
'11', '12', '13', '14', '15', '16', '17', '18', '19', '20', '21',
'22', '23'], dtype=object)

00 (606, 117) (60, 117) (609, 117)
01 (611, 117) (61, 117) (634, 117)
02 (610, 117) (61, 117) (637, 117)
03 (602, 117) (61, 117) (624, 117)
04 (611, 117) (61, 117) (626, 117)
05 (610, 117) (61, 117) (628, 117)
06 (594, 117) (58, 117) (572, 117)
07 (611, 117) (60, 117) (632, 117)
08 (610, 117) (61, 117) (633, 117)
09 (598, 117) (61, 117) (589, 117)
10 (612, 117) (61, 117) (635, 117)
11 (613, 117) (61, 117) (631, 117)
12 (591, 117) (60, 117) (612, 117)
13 (613, 117) (61, 117) (634, 117)
14 (613, 117) (61, 117) (637, 117)
15 (598, 117) (60, 117) (614, 117)
16 (609, 117) (61, 117) (629, 117)
17 (608, 117) (61, 117) (630, 117)
18 (588, 117) (59, 117) (610, 117)
19 (607, 117) (61, 117) (631, 117)
20 (611, 117) (61, 117) (637, 117)
21 (599, 117) (61, 117) (619, 117)
22 (612, 117) (61, 117) (637, 117)
23 (612, 117) (61, 117) (633, 117)

Create feature and label sets

We now need to seperate our training (pri_train), evaluation (pri_eval) and test (pri_test) dataframes into their feature and label set components. Here we create a small function to do this for us. The training, evaluation and test datasets contain both the raw and the normalised (_zscore) data. We only wish to pass the normalised feature data to the ML routine so we simply select the relevant zscore features. The function also takes in which target label parameter, which we use to select the column we want to model for.

Define our XGBoost model - hand it our train, eval and test sets

So we are now ready to hand each of our 3 datasets for each hour from 0 to 23 to the XGBoost routine to generate a model for each hour. We define our XGB model using the Regressor class and a squared error objective function. We then pass the training feature and label sets to the fit function along with our evalution set. Once the modeel is generated, we will then store the feature importance graph and then generate a prediction for each hour and then calculate the Mean Squared Error (MSE) for each prediction versus the observed label.

Generate Visual Output

We have just created 24 hourly xgboost models, we now want to plot our feature importance (the higher the F-score the more important the feature) for each hourly series, as well as the accuracy of our ML generated predictions - using in our case the mean squared error (MSE) calculation we created above. We can easily do this is 2 steps.

Initial Observations

From the line charts above, we can see that the generated models are able to predict some of the downward spikes in spot energy prices. The MSE varies from 0.73 (hour 01) to 1.56 (hour 6). There does seem to be quite some variation in accuracy even with this very simple and limited model. However, none of these are particularly accurate, MSE of 0 being perfect and MSE of 1 indicating 1 standard deviation.

We can see for example that on some days prices went much lower than our prediction - this could be due for example to a day which was both windy and sunny (remember these energy sources are not stored but dumped to the grid - leading to the negative prices - the negative prices discourage use of thermal sources). We did have one feature included for solar energy - but it was for only one location in the centre of the country, as opposed to wind which we had much greater coverage for. Conversely, there are other days where the prediction was much lower than the observed, so perhaps it wasn't that sunny at all locations or thermal capacity had been dialled back etc. So we could expand our model to include new more complete features such more solar locations, thermal features and such.

Our model does not really tell us anything about upward spikes as we have only included one demand side type of feature - air temperature - and only for one station in the centre of Germany - in our feature matrix. Again we could expand our model with features that might help us with this demand side modelling for example.

Another next step might be to use only a subset of the highest F-score features and re-run the models to see if that can improve things, or indeed look at a suite of different classes of algorithm - we have only used xgboost but there are many other that we could implement within the main iterative routine I have given above.

Further, we have been looking only at observed values. Remember, one of the features we are using is actual Wind Output - we wouldn't know this ahead of time. So clearly we need to know about forecast data. We know that energy markets are based around forecast data. So another next step would be to incorporate all of that data but that is beyond the scope of this introductory article.

Summary

A lot of ground was covered in this article so let's summarize. We retrieved a set of physical asset (windfarms in our case) RICs from our excellent new RDP Search API. We used these to request some reference-type data for these windfarms using our Eikon Data API. We then created a new geometry 'Point' datatype where we stored locational co-ordinates for our windfarms. We followed this by loading some open source wind data, wind output data from our Point Connect repository and joined these together, creating 2 additional features along the way. We also loaded some open source solar energy and air temperature data. We then retrived the locations of our weather monitoring stations and as before, created a new geometry datatype for those. We were then able to calculate the closest weather monitoring station for each of our windfarms using the scipy.spatial cdist routine - though we didn't use it in ths introductory model other than for visualisation - I included it as it might prove useful were you to develop more complex models. We then downloaded a shape file for Germany, converted our pandas dataframes to GeoPandas dataframes and ploted overlays for both windfarms and wind monitoring stations. Next we read in some price data (EPEX Spot) to use for our label data, shifted it forward one hour to align with our local time (CET) features and we then started to normalise our datasets using z-scores. After joining our feature and label dataframes, we then added a feature called hour which tells us which hour the timestamp refers to. We then split our data into training, evaluation and test sets and defined our feature and label sets for each of the 24 hours in the day. Next we defined 24 different XGBoost ML models (one for each hour of the day), handing it the appropriate training and evaluation sets and ploted our feature importance F-scores. Finally, we used our model to generate a label prediction based on our hold-out or test feature matrix and plotted that against the observed label values.

Whilst a simple and far from complete example, I hope this article can provide you with an idea of how to approach this area and a practical codebase to experiment and explore further.

Further Resources

Further Resources for Eikon Data API

• Overview
• Quick Start
• Documentation
• Tutorials
• Q&A Forums

For Content Navigation in Eikon - please use the Data Item Browser Application: Type 'DIB' into Eikon Search Bar.

Further Resources for Refinitv Data Plaform (RDP) APIs

• Overview
• Quick Start
• Documentation
• Downloads
• Tutorials
• Q&A Forums

For Content Navigation in RDP - please go here

For more information on the Point Connect service please go here

Further Resources for Open Source German Weather Data kindly made public by DWD

Hourly Wind Data

Hourly Air Temperature - I used Station 01266, and field TT_TU

Hourly Solar Readings - I used Station 00183, and field FG_LBERG

Source of geospatial base data:

Surveying authorities of the Länder and Federal Agency for Cartography and Geodesy (https://www.bkg.bund.de).

Source of satellite data: EUMETSAT (https://www.eumetsat.int) NOAA (https://www.noaa.gov)

Further Resources for Shapefiles

German SHP Map