Improving wind & solar nowcasting using meteorological
satellite image stream
submitted in partial fulfillment for the degree of master of science
Tam´
as Elekes
master information studies
data science
faculty of science
university of amsterdam
2019-07-03
Internal Supervisor External Supervisor Title, Name Ana Lucic, MSc Fons de Leeuw, MSc Affiliation University of Amsterdam Dexter Energy Services Email a.lucic@uva.nl fons@dexterenergy.nl
1
Introduction
The electric grids around the world are facing a major transformation in the last decade. Rising electricity consumption — along with global awareness of cli-mate change — is leading energy companies to invest in renewable energy resources. Although renewable resources like solar and wind energy are cleaner and a better energy solution from a sustainability point of view, connecting these power sources to electric grids have raised concerns of utility companies re-cently. This is due to the possible negative impacts of power fluctuations generated by these systems on network operations. These fluctuations can also lead to unstable operations of an electric grid prior to fault conditions, high power swings in the feeders, and un-acceptable power outages. Thus, energy suppliers need to plan how much electricity they need to have available at a given time. When the electricity con-sumption does not match the generation a so called imbalance price is calculated. Energy farms want to make their imbalance costs less. Because renewable energy supply greatly depends on weather conditions as well on private consumption, solar and wind fore-casts are an important tool for planning electricity usage and generation. Dexter Energy Services pro-vides these forecasts to energy farm operators. The aim of this thesis is to gather what are consid-ered as best practices about wind and solar forecast-ing systems. We deployed a baseline model based on state of the art (SOTA) techniques using Numerical Weather Predictions (NWPs).
NWPs are mathematical models that predict the weather based on current weather conditions. Mete-orological institutes run it on regional or global scale, with different accuracy. The result is a large database of quarterly weather metrics like temperature, solar radiation, wind speed and air humidity. These pre-dictions are available for 1-3 days ahead.
The baselines are set up with simple machine learn-ing techniques such as feature engineerlearn-ing, regression and power curve fitting. Our main contribution is ex-amining possible extensions of these SOTA methods by utilizing the satellite imagery stream from EU-METSAT. This satellite stream consists of images recorded from a geostationary viewpoint above
Eu-rope captured at different wavelengths.
Nowcasting is a technique for very short-range fore-casting that maps the current weather, then uses an estimate of its speed and direction of movement to forecast the weather a short period ahead. The satel-lite imagery is used for experiments with optical flow and deep neural network based methods. The objec-tive of these techniques is to create accurate nowcasts up to 6 hours and use them as an input for a short term forecasts.
We evaluate in terms of both accuracy of the forecasts as well as monetary savings.
1.1 Research questions
1.1.1 Main research objective
• The goal of the research is to determine whether image nowcasting techniques such as optical flow can offer better input features for solar energy yield forecasts than those based on NWPs in the time window of intraday electricity trades, 0-6 hours.
1.1.2 Research Questions
• RQ1: What is the performance of current state of the art models for day-ahead solar and wind forecasts?
• RQ2: How can geostationary satellite images be preprocessed and pipelined to be able to use them in machine learning applications?
• RQ3: How does the widely used optical flow method to predict satellite image stream per-forms?
• RQ4: What features can be extracted from now-casted satellite images for solar forecasting? • RQ5: Does including the features from RQ4
along with features from NWPs improves short term solar forecasting? Can using only the cre-ated nowcasted features result in similar accu-racy?
2
Preliminaries
2.1 Problem description
Forecasting renewable energy production has some similarities to other forecasting domains. For exam-ple, seasonal effects similarly as in financial markets are present, but in renewable energy generation it is due to icing, maintenance, amortization of assets and of course the seasonality of the weather.
Electricity farms have monetary interest in more ac-curate forecasts, since more accuracy means less im-balance costs for operating the farm. Imim-balance costs are a penalty cost determined by a central institute to compensate for the over or under delivery of elec-tricity for the grid. For Netherlands and a large part of Germany this grid operator institute is TenneT. Imbalance costs are centrally regulated market that aims to maintain the available electricity on the de-sired level. See 2.2.1 for details.
2.1.1 Photovoltaic power prediction
A good metric on solar energy production is the amount of solar irradiance exposed to the photo-voltaic panel. There are three types of solar irra-diance:
• Direct Normal Irradiance (DNI) is measured at the surface of the unit perpendicular to the sun. It excludes diffuse radiations from nearby sur-faces, only includes the direct beam from the sun. Atmospheric losses are included, but they are dependent on cloud cover, air pollution, moisture and others.
• Diffuse Horizontal Irradiance (DHI) is measured as the light scattered by molecules and particles in the atmosphere, so not directly coming from the sun.
• Global Horizontal Irradiance (GHI) is the total irradiance from the sun on a horizontal surface on Earth. It accounts for the zenith angle of the sun.
GHI = DHI + DNI · cos(z), (1) where z is the zenith angle of the sun.
More details on how this is used for setting up the solar baseline model can be found in 5.2.1.
2.1.2 Wind energy and power curves Any wind forecasting system has to account for im-portant dimensions of the optimized windfarm, such as • Location • Power curve • Rotor diameter • Hub height • Availability / uptime • Weather models • Surface roughness
The biggest factor in wind power generation is wind speed, but converting it to actual wind power is bet-ter for modelling the actual power that can be har-vested for it. Calculating wind energy from wind speed:
Pwind=
1 2%Av
3, (2)
where % is air density, A is rotor area, v is wind speed The theoretical maximum power that can be ex-tracted from the wind, is approximately 59.3% and it’s called Betz’s law. Find the proof in Betz (1966). It’s based on the idea, that if 100% of the kinetic en-ergy would be used the air would have to stop after the turbine. It would prevent further air passing the blades and it would stop spinning.
Wake effect is another factor influencing the power predictions. It happens when several wind turbines are set up in a row. When the wind blows parallel to this row, the first turbine can take from the full power of the wind, but the subsequent turbines will generate less due to the wind being slowed down. It is not directly accounted for in our model, but since the wind prediction model is recalculated for different wind directions, the effect of this force can be seen in the predictions.
Wind turbines have a cut in speed, which is the mini-mum wind speed needed to start the turbines to gen-erate power. For safety reasons there is also a cut out speed which mechanically slows down or stops the
turbine if the wind speed is over this threshold. See 5.2.2 for details of the wind baseline model.
2.2 Improving short term forecasts with nowcasted images
The forecasts relying on NWPs have an inherent error due to the inaccuracies of the spatio-temporal shifts and errors in peak weather values. Also some localized types of storms are hard to predict. Thus, there is a potential to improve forecasts using a real time and accurate weather data source. In this work, experiments will be conducted to test if good quality short term weather forecast is achievable. These are then converted to input features to electricity yield forecast systems and it will be investigated whether these modifications lead to better performance. Therefore the pipeline (what partly answers RQ2) of our system is:
Satellite image stream → Preprocessing → Nowcasting images → Extract features → Use features in forecasting systems
2.2.1 Imbalance costs
Power cannot be stored efficiently, hence supply and demand always have to be in balance. To match supply and demand, market players in the electric-ity system have to plan ahead the buying and selling of electricity. Issuing imbalance costs are a regulatory method to maintain the balance between the in flow-ing electricity and the withdrawal of electricity from the electricity grid. The following explanation is only accurate for a specific company’s (TenneT) imbalance cost calculations. TenneT works in Netherlands and some parts of Germany. Other regions have other methods of calculating imbalance price, but in every case it reflects the price of the electricity that has to be added or removed from the grid.
There are three roles in the balancing system: • Transmission System Operator (TSO):
The grid operator that is responsible for the sta-ble frequency over the grid. They are respon-sible for monitoring, maintaining and restoring the balance between supply and demand of elec-trical power in it’s area. TenneT is this operator in this case.
• Balance Responsible Party (BRP): All par-ties that connect to the electricity grid must be an accredited BRP. Each BRP obliged by law to send a schedule to the TSO for each day. They are financially responsible for it’s imbal-ance, that is the deviance from this trade sched-ule, and pays or receives the imbalance price for each imbalance settlement period (ISP). The ISP is fixed at 15 minutes.
• Balancing Service Provider (BSP): A mar-ket party from which the TSO can ask frequency restoration reserves for its balancing task. Their main purpose is to ensure operational security on the grid. The used resources are called Fre-quency Restoration Reserves (FRR). They are contracted to be able to supply or to be able to store excess production of electricity.
See how the determining day-ahead and imbalance price at Section 8.
Figure 1: Supply demand curve that is used for merit order calculating day ahead electricity prices
3
Related Work
3.1 Wind forecasting
Pinson (2006) has an extensive overview of many as-pects of wind power forecasting. The author pro-posed SOTA techniques that are using NWPs and historical output values, the 15 minute time series data from the farm, containing energy yield values over that quarter. The author states that the most common reference method for both solar and wind forecasting is persistence. This is a naive predictor, stating that weather based production will have sim-ilar outputs as the current:
pt+1= pt (3)
Due to the fact that weather is unlikely to drastically change in the short term, it can produce good results for up to 3-6 hours. It’s only worthwhile to develop a solar or wind forecasting method, if it’s able to beat this reference. We used this reference for the nowcasted satellite images.
In general forecasting is considered in three different timescales:
• From milliseconds up to a few minutes (used for active control of the turbine)
• For the following 48-72 hours (for management of energy trading and imbalance costs)
• For longer time scales up to 5-7 days ahead (for planning of the cost effective maintenance of the farms)
In this thesis we only forecast for up to 6 hours, to be used in the intra-day electricity market.
The physical approaches are discussed, where the whole environment of the turbine is modelled. Physi-cal attributes such as surface roughness, surrounding objects, whether if the park is onshore/offshore. The wind speed is converted to wind power and theoretical wind power curves are used. Recently empirically derived power curves are used more often. The statistical approaches are based on one or more models that couple historical values of power and wind power measures. The models create a mapping from weather values such as wind
speed, temperature, air pressure to historical power values of the farm. The physical properties are not decomposed and accounted for, but expertise on the problem is needed, in order to be able to choose the right explanatory variables. A train set is used for optimizing the and a test set is used to evaluate models. The models that uses the analyst’s expertise on the task are called structural models, such as linear regressions or decision trees. On the other hand, when little to no domain expertise are built into the model are called black-box models, like neural networks or support vector machines.
Marˇciukaitis (2017) proposes a technique to estimate the actual power curve of a turbine or a park. It is done by defining four parametric formulas and it is preceded by data cleaning based on a method by Tukey (1977). Each formula was scored based on the actuals, and then a comparison was made between the resulting derived power curves. They also propose partitioning the datasets by wind direction. Techniques from Tukey (1977) were used to visualize the available data for better understanding were used to visualize the data and get rid of the outliers. The wind forcasting system relies greatly on this work, however further data processing and techniques were used. See Section 5.2.2 for details.
3.2 Solar forecasting
See Emil Isaksson (2018) for a research on different photovoltaic panel output predictions. They used NWPs as input features and benchmarked five types of machine learning techniques. They predicted for up to six hours, due to their NWP source was pro-duced four times a day, only for six hours. In the fea-ture engineering they account for the sun zenith an-gle and they address the inaccuracies in the weather grid by introducing lagged variables. In our baseline model seen in Section 5.2.1 we also addressed these but with more care. Not only using the Suns zenith angle but calculating GHI from it, and the shifts in weather data with the moving windows technique. The models they tried:
• Autoregressive integrated moving average (ARIMA)
• Lasso regression • K-Nearest Neighbors
• Gradient Boosting Regression Trees (GBRT) • Artificial Neural Nets (ANN)
The baseline models were Persistence and Climatol-ogy was used. ClimatolClimatol-ogy is based on using the aver-age of historical data as the forecasted value. GBRT and ANN performed the best, but GBRT takes much less time to run. Based on this work we choose XG-Boost (to see accuracy) and multivariate regression (for getting feature importance) in Section 6.3. 3.3 Nowcasting satellite images
Optical flow technique was first proposed in Bruce D. Lucas (1981) and since there have been various usages in generating next steps in moving imagery and it is the current state of the art technique to do so. It is also a known technique for using its results in weather nowcasting. Georgy Ayzel and Winterrath (2018a) reviews some common optical flow algorithms and establish a set of benchmark procedures. See Section 5.3 for details on the used algorithms and the experiment with the satellite imagery. For benchmarking the software reposi-tory Georgy Ayzel and Winterrath (2018b) was used. Most of the research regarding meteorological satel-lite image nowcasting seem to focus on precipitation since that is the most important for weather alerts and safety systems. However using them to obtain more diverse weather features or as am input for other forecasting systems are a novel approach.
4
Data
4.1 NWPs
• KNMI HARMONIE, NWP day ahead predic-tions (KNMI (2013))
• European Centre for Medium-Range Weather Forecasts (ECMWF), NWP day ahead predic-tions (ECMWF (2019))
NWPs coupled with historical data of energy out-puts of renewable energy parks as a time series dataset. Since the weather values are from a day ahead weather simulation there is an error built in
the data that is very hard to filter, because it can be shifted in time, have inaccuracies in peak weather values. One solution for this, is that different window sizes were used around the examined farm and calcu-lated mean/min/max values over that window. Using those could account for spatio-temporal inaccuracies in the weather data.
As previously mentioned, NWPs are mathematical models that include current weather conditions, and calculate different weather parameters for a 2D grid over a specified area, in a time-series manner. Two NWP datsets were used with success, the Dutch KNMI and the European ECMWF. ICON-EU (Wet-terdienst (2015)) was also tried, it’s accuracy was not satisfactory. Every morning they run their models and make it available for further processing. For the purpose of solar and wind predictions a post process-ing is done, with calculatprocess-ing mean/min/max proper-ties over different window sizes of the local weather. See Figure 8 and Figure 7 for solar radiation value and Figure 9 for wind speed count at a farm loca-tion. It’s easy to spot some physical properties from plotting and exploring the data such as:
• Shorter days in the winter
• Less solar radiation in the winter due to more clouds and the angle of the sun
• Noise from the sea mask because less cloud over the ocean causing bigger window have higher val-ues of radiation if containing segments of the sea • General wind direction
• No output of solar panels in the night
4.2 Dexter datasets
Dexter provides actual load values (which are the tar-get value in forcasting) from solar and wind farms indexed with the same time-series as NWPs. This is an ideal training set for machine learning models for forecasting load values.
4.3 Satellite imagery
EUMETSAT’s “Rapid Scan High Rate SEVIRI Level 1.5 Image Data - MSG” is an image stream provided by geostationary satellites orbiting earth. The cap-tured raw data is processed and constitutes one of the main products of the Meteosat Second
Genera-tion (MSG) system. The designaGenera-tion ‘Level 1.5’ cor-responds to image data that has been corrected for all unwanted radiometric and geometric effects, has been geolocated using a standardised projection and it has been calibrated and radiance-linearised. The Level 1.5 data is suitable for the derivation of me-teorological products and further meme-teorological pro-cessing. (EUMETSAT (2009)) The SEVIRI instru-ment is designed to produce the image of the Earth from a spinning geostationary satellite.
Figure 2: The principle of SEVIRI Earth imaging.
It is a time series if of images every 5 minutes. The images are provided in twelve channels, correspond-ing to spectral bands. They vary in spectral char-acteristics, dynamic range, operating temperature of the detectors and number of detectors simultaneously acquiring image information. The images are clipped to only include Europe. Since they are geostation-ary photos of a geoid shape they require a computa-tionally heavy preprocessing step in order to convert them to a near equidistant version, the transverse mercator projection. This is important because if the neural net is being fed with the original distorted image, the learned filters of a convolutional neural net would not work on any part of the image, just where it was learned. Also useful for optical flow based transformations to have pixel to real distances same across the image. See details of the preprocess-ing in Section 5.1
To analyze solar and wind properties, the presence
Figure 3: The transverse mercator projection. Source: http://richarddmorey.org/map/
and movement of water vapor moisture in the up-per and middle levels of the atmosphere can be used. There is two channels with water vapor. Further re-search can be identifying useful features in different spectral bands. See Table 11 for details on spectral bands. Dataset description can be found in EUMET-SAT (2017). See description of spectral bands in Sec-tion 8.
The aim with these images is to use them as input for image nowcasting methods and forecast the following 6 hours worth of images with them.
5
Methodology
The hypothesis is that it’s possible to improve the baselines for short forecasts as 0-6 hours by creating features that are based on satellite image nowcasts. 5.1 Preprocessing the images
This section is answering RQ2. First step with the EUMETSAT dataset is to convert from its own geo-stationary projection with an offset to the east by 9.5°to an equidistant projection. The perfect equidis-tant projection is unfortunately impossible, as the earth is a sphere and the target projection is a plane. The solution is to take a projection where the target area has very low distance error on average for the examined area. The target projection used is trans-verse mercator. See Figure 3. Europe is in the middle line of it where there is the least distortion.
view-point. Since a pixel always represents the same loca-tion a mapping for each pixel to a latitude longitude values were created. PROJ.4 projections library was used to convert coordinates to meters.
The projections were defined with projections strings that PROJ.4 projection library can use. See PROJ contributors (2019) for the documentation of this li-brary, and the references to build a projection string. The projection string is a key/value list that can be used to initialize a PROJ.4 projection. To define the image itself, the original projection was defined like:
+proj = geos + lon 0 = 9.5 + h = 35785831.0 (4) where geos is geostationary viewpoint, lon is degrees to shift to the east, and h is height of the view point above the earth in meters.
Transverse mercator, the target projection:
+proj = tmerc+lat 0 = 0+lon 0 = 9.5+ellps = intl + units = m + no defs (5) where tmerc projection is transverse mercator, cen-tered at lat langitude and lon longitude (same place as the satellite is), the ellipse of the earth is defined as International 1909 (Hayford), units are in meters. The +no defs parameter means to ignore any other parameters from default settings files.
Using a pixel to meters ratio from the data descrip-tion EUMETSAT (2017) the pixel values for each co-ordinate is calculated then placed the original image on this grid, with linearly interpolating between the pixels of the original image. In the end of the prepro-cessing, the transformed images and a lookup table is available that makes two way conversion possible, from pixel coordinates to geolocation and back. This is useful at the feature extraction phase.
5.2 Baseline methods
To answer RQ1 two separate systems were set up for day-ahead forecasting. They are using NWPs and historical output values and based on current indus-try best practices.
5.2.1 Solar-specific Baseline
The approach taken for solar prediction was to cal-culate solar radiance values, including correction of the angle of the panel. For these calculations pvlib was utilized, see William F. Holmgren and Mikofski (2018) for details. pvlib can be queried for meteoro-logical information such as the position of the sun, zenith and azimuth positions and calculate solar ra-diation values that were explained in Section 2.1.1. There are many factors that have an effect on the efficiency of the solar panels; cooler panels work bet-ter (temperature); obstructing mabet-terials also impede performance (pressure, humidity, objects in the en-vironment). Following this logic cold air and wind cools panels, slightly improving their energy yield. Humidity can collect as water droplets on the panel and not only degrading their performance, but but also amortizing them. Wind can also dry panels, so it can reverse the effects of humidity.
Since time of the day, and even season has effect on the efficiency of the panels, categorical features were created. Namely a day got divided to 96 quarters and used as a categorical feature. Similarly using the number of the month the data point is used as fea-ture. Then a regression model was created from the available features and their higher polynomials. Or-dinary least squares (OLS) and automatic relevance determination (ARD) regressions were used.
5.2.2 Wind-Specific Baseline
The wind forecasting system made during this thesis relies on empirically derived power curves, using non linear curve fitting. The main idea is the estimation of the power curve by fitting a parametric formula based on the findings in Marˇciukaitis (2017) on the cleaned historical power output (measured in watts). If there is no historical output (for example in the case of a newly set up park) a general power curve, the power curve supplied by the manufacturer could be used. As recordings of the actual electricity output is gathered there is a forecast blending period where the actual prediction is a weighted mix of the general forecast and the forecast based on limited historical data.
Figure 4: Comparison over time of the actual output (red) reference model (green) and the baseline model (blue) created
Different data enriching methods such as scaling the wind speed to hub height, factoring with air density were implemented. Categorizing wind direction and keeping track of a different power curve per direc-tion. The number of directions is arbitrarily chosen however further development is possible for a dynami-cally chosen, overlapping wind direction intervals can achieve better accuracy.
To incorporate air density just by accounting it in the wind power (as in 2) didn’t result in satisfactory improvement so different power curves would have to be estimated for ranges in air density, but when the data gets fragmented into parts we’ll end up with in-complete datasets (e.g. only observations where wind speed was high and temperature was high or other way round).
The power curve is modeled based on Marˇciukaitis (2017). The parametric formula used for fitting it:
p(v) = pmax(1 + (
a v)
b)c, (6)
where a, b, c is coefficients of the curve fitting, v is the wind speed after corrections.
In Figure 10 a general curve from a turbines specifi-cation can be seen and a fitted one. The fitted one captures the environmental effects in the curve, and even shows that the maximum peak output is less than the specified (this can be because the turbine itself is using electricity or amortization).
5.3 Optical flow algorithm
To answer RQ3 an experiment was designed to gen-erate nowcasted images with optical flow, and mea-sure accuracy by comparing them to the observed images.
Figure 5: Extracting different areas from nowcasted images. Color scale is normalized to the min/max values in each crops.
Optical flow method assumes the motion of observed objects is dependent between time step t and t+1. For the experiments with Optical Flow algorithm and their evaluation was made based on the work Georgy Ayzel and Winterrath (2018a). There are two underlying methods implemented: (1) Sparse and (2) Dense. These methods have 1-1 variations, Sparse, SparseSD, Dense and DenseRotation. In the Sparse group the key idea is to identify corners in t-k steps, where k is the number of previous ob-servations and track these corners to t step, then ex-trapolate each pixel for all n nowcasted future steps. In the Dense group methods utilize the Dense Inverse Search algorithm (see Till Kroeger (2016) for the al-gorithm) for estimating the velocity of each image pixel based on two consecutive satellite images. It is optimized in regards for time complexity, so it is useful for real time applications. Sill, Dense methods take more time than the Sparse ones but also over performs them.
As a baseline model Persistence was used, which is a trivial, but surprisingly well performing model. See equation and concept at Equation 3. See them in de-tail and benchmarks in Georgy Ayzel and Winterrath (2018a). The used source code can be found in this repository, Georgy Ayzel and Winterrath (2018b). 5.4 Exporting nowcasted features
The concept in creating input for ML algorithms from the nowcasted images is to generate similar features as the preprocessed NWPs have. Over different win-dow sizes mean/min/max values of the gray-scale pixel values were saved along the raw pixel values, see Figure 5. RQ4 was assessed with this feature extraction.
5.5 Experiment to use nowcasted features To evaluate the how effective the features exported from nowcasted images were an experiment was set up. Fifteen solar farms location and historical output was used with 349 days of training data. The most important features from NWPs were determined so-lar predictions. Two models were trained three times, with different feature sets, with NWP+nowcasts and separately NWPs and nowcast features. This answers RQ5.
6
Results
6.1 Results of baseline methods
The results of the baseline models answer RQ1. Wind forecasts were compared SOTA methods cur-rently used in the industry, for solar only our accu-racy is available, but it is used in real life solar farm optimization.
6.1.1 Wind results
Figure 4 shows an example time period of the baseline wind forecasts and an available SOTA system. In the first half of 2019 our system performed 2.2% better in terms of normalized MAE. Comparing imbalance costs for the same period, and same wind farms it generates about 3300€ less costs, about 3% better in monetary standpoint.
6.1.2 Solar results
A baseline model was created based on the meth-ods listed in Section 5.2.1. Figure 12 shows actuals and predictions plotted over the hours of a day. It shows that it generally over predicts in the morning, but fairly good other times. Overall the final model achieves around 24% normalized MAE. In monetary terms in the first half of 2019 about 4000€ was payed per Gigawatt hours for the predicted farms. Given that there are no other major costs in generating a unit of solar electricity this is still a cheaper option than non-renewable types of electricity.
6.2 Nowcasted satellite images
Six hours of subsequent satellite images were created with all methods, and the used spectral bands, over 8 months, 37 chosen days. The training data was 2 images, with 15 minute difference at 8 in the morn-ing, except for Sparse, that has 12 images, 3 hours. All images were gray-scale with one channel. Each evaluation metric is the average over these results. Spectral bands with wavelength where water vapour is reflecting the most (WV62 & WV73) have very low error, what is a good indicator if clouds are above an area. Three evaluation metrics were calculated from the nowcasted images at 1h/3h/6h times, using the actual images for the same period as test dataset. See partial results in Table 1 and for channels and evalu-ation metrics see Table 4, Table 6 and Table 5. Dense models perform the best and DenseRotation over per-forms slightly Dense. The nowcasted images took less then a minute to produce with these methods. Per-sistence also produced similar accuracies, and over-performed methods in the Sparse group.
6.3 Using features from nowcasted images in solar panel yield prediction
An XGBoost (XGB) (XGBoost developers (2016)) model was set up, and after training the impact of each feature could be determined as Figure 13 shows. ”pix” prefix is nowcasted and ”harm” means nwp fea-tures. As the image shows nowcasted features defini-tively have an impact on the models performance. A regression model Automatic Relevance Determina-tion (ARD) Regression Pedregosa et al. (2011) was also tried. This model performs worse than XGBoost, but the impact of the nowcasted features can be seen at Figure 6. This plot shows the features ordered by impact on the model and the color scale corresponds to the value of the datapoints.
With both models while nowcasted features have ex-planatory value, and can improve a model however using them alone didn’t result nearly as good accu-racies as NWPs do.
1/3/6 hours Dense DenseRotation Sparse SparseSD Persistence IR97 .068/.119/.162 .068/.12/.164 .112/.173/.235 .111/.187/.259 .086/.14/.182 VIS8 .149/.271/.324 .149/.279/.326 .169/.274/.475 .201/.346/.405 .17/.292/.332 WV73 .045/.079/.117 .045/.079/.115 .087/.147/.208 .088/.151/.214 .064/.102/.137
Table 1: Normalized MAE averages from a water vapour, visible and infra red channel
xgb nwp xgb nwp now xgb now ard nwp ard nwp now ard now IR97 0.270933 0.264292 0.487953 0.369093 0.366500 0.557702 WV73 0.270264 0.267837 0.502214 0.369023 0.366347 0.557565 VIS8 0.270387 0.265134 0.387808 0.369099 0.365606 0.482962
Table 2: XGB and ARD model accuracies with the different features from NWPs and nowcasts. A water vapour, visible and infra red channel
Figure 6: Feature importance in an ARD regression based model.
7
Future work
Making more accurate nowcasted images would definitively make better forecasts. Using neural net-work based models such as TrajectoryGRU (Shi et al. (2017)) could generate more accurate images there-fore better extracted input features. There are works already converting pixel values to solar radiance from
EUMETSAT imagery. It might act as input to the original solar baseline model, or just be a more exact representation than the raw pixel values. See details in Kada Bouchouicha and Aoun (2016). Psychical properties of spectral bands could be utilized.
8
Conclusion
Using realtime data and nowcasting methods have a great potential to help energy transmission to renew-able sources. RQ1 was answered by creating baseline models that matched industry standards. The pre-processing and the resulting pipeline that was defined in Section 2.2 gives answer to RQ2. Optical flow based nowcasting was tried to answer RQ3 however it definitely could help to try other more advanced techniques to generate the nowcasted images. Fea-ture extraction was also tried (RQ4), creating values that is a suitable extension for the processed NWP datasets already used, however there is still research to be done to utilize the characteristics of the spec-tral bands. The final experiment that referred RQ5 tried the created features and did result in minor im-provement in short term solar predictions, however to use only satellite images in solar predictions needs further research.
Acknowledgements I want to thank the help of my supervisors, Ana Lucic, Fons de Leeuw, my teachers in the University of Amsterdam and my colleagues at Dexter Energy Services.
References
Albert Betz. Introduction to the theory of flow ma-chines. Oxford: Pergamon Press, 1966.
Pierre Pinson. Estimation of the uncertainty in wind power forecasting. PhD thesis, CEP -Centre ´Energ´etique et Proc´ed´es, 2006. URL https://pastel.archives-ouvertes.fr/ pastel-00002187.
Linas Martiˇsauskas Benas Jokˇsas Giedrius Geceviˇcius Athanasios Sfetsos Marˇciukaitis, Inga ˇZutautait˙e. Non-linear regression model for wind turbine power curve, renewable energy. 2017. URL https://www.sciencedirect.com/science/ article/pii/S096014811730544X.
J.W. Tukey. Exploratory data analysis. 1977. URL https://www.stat.berkeley.edu/~brill/ Papers/EDASage.pdf.
Mikael Karpe Conde Emil Isaksson. Solar Power Forecasting with Machine Learning Techniques. PhD thesis, KTH ROYAL IN-STITUTE OF TECHNOLOGY SCHOOL OF ENGINEERING SCIENCES, 2018. URL https://www.math.kth.se/matstat/ seminarier/reports/M-exjobb18/180601f.pdf. Takeo Kanade Bruce D. Lucas. An itera-tive image registration technique with an application to stereo vision. 1981. URL https://pdfs.semanticscholar.org/51fe/ a461cf3724123c888cb9184474e176c12e61.pdf. Maik Heistermann Georgy Ayzel and Tanja
Winterrath. Optical flow models as an open benchmark for radar-based precipitation now-casting (rainymotion v0.1). 09 2018a. URL https://www.geosci-model-dev.net/12/1387/ 2019/gmd-12-1387-2019.pdf.
Maik Heistermann Georgy Ayzel and Tanja Win-terrath. rainymotion: Python library for radar-based precipitation nowcasting radar-based on optical flow techniques, 2018b. URL https://github. com/hydrogo/rainymotion.
KNMI. Harmonie modelgegevens, 2013. URL http: //projects.knmi.nl/datacentrum/catalogus/ catalogus/content/nl-nwp-harm-grid-p1. htm.
ECMWF. Nwp data, 2019. URL https://www. ecmwf.int/en/forecasts.
Deutscher Wetterdienst. Icon-eu nwp forecast data, 2015. URL https://www.dwd.de/ EN/ourservices/nwp_forecast_data/nwp_ forecast_data.html.
EUMETSAT. Rapid Scan High Rate SEVIRI Level 1.5 Image Data - MSG, 03 2009. URL https://navigator.eumetsat.int/product/ EO:EUM:DAT:MSG:HRSEVIRI.
EUMETSAT. MSG Level 1.5 Image Data Format Description, 09 2017. URL http://www.eumetsat.int/website/
wcm/idc/idcplg?IdcService=GET_FILE& dDocName=PDF_TEN_05105_MSG_IMG_DATA& RevisionSelectionMethod=LatestReleased& Rendition=Web.
PROJ contributors. PROJ coordinate transformation software library. Open Source Geospatial Founda-tion, 2019. URL https://proj.org/.
Clifford W. Hansen William F. Holmgren and Mark A. Mikofski. pvlib python: a python pack-age for modeling solar energy systems. Journal of Open Source Software, 3(29), 884, 2018. URL https://doi.org/10.21105/joss.00884. Dengxin Dai Luc Van Gool Till Kroeger, Radu
Tim-ofte. Fast optical flow using dense inverse search. 2016. URL https://arxiv.org/pdf/ 1603.03590.pdf.
XGBoost developers. XGBoost, an optimized dis-tributed gradient boosting library, 2016. URL https://xgboost.readthedocs.io/en/latest/. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Pas-sos, D. Cournapeau, M. Brucher, M. Perrot, and
E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825–2830, 2011.
Xingjian Shi, Zhihan Gao, Leonard Lausen, Hao Wang, Dit-Yan Yeung, Wai-kin Wong, and Wang-chun Woo. Deep learning for precipitation nowcast-ing: a benchmark and a new model. In Advances in Neural Information Processing Systems, 2017. URL https://arxiv.org/pdf/1706.03458.pdf. Nour El Islam Bachari Kada Bouchouicha,
Ab-delhak Razagui and Nouar Aoun. Estima-tion of hourly global solar radiation using msg-hrv images. International Journal of Applied Environmental Sciences, 2016. URL https://www.ripublication.com/ijaes16/ ijaesv11n2_01.pdf.
TENNET. Imbalance Pricing System, 03 2019. URL https://www.tennet.eu/fileadmin/ user_upload/SO_NL/ALG_imbalance_pricing_ system.doc.pdf.
STEPHEN TJEMKES DIETER JUST JOCHEN KERKMANN SERGIO ROTA ALAIN RATIER JOHANNES SCHMETZ, PAOLO PILI. An introduction to meteosat second generation (msg). 2002. URL https://eumetsat.int/ website/wcm/idc/idcplg?IdcService=GET_ FILE&dDocName=PDF_SCI_BAMS0702_INTRO-MSG& RevisionSelectionMethod=LatestReleased& Rendition=Web.
Determining day-ahead and imbalance price
The general principle BRPs aim to reduce imbalance, and instead of relying on the TSOs regulations, it’s in their best interest to actively keep the balance. Prior to the delivery day, each BRP submits a schedule for the next day about their expected input to the grid per ISP. The TSO checks whether the expected electricity transactions add up to zero, so that the supply and demand of electricity is balanced for every ISP of the delivery day.
During the delivery day BRPs should deliver electricity regarding their submitted schedules. If a power imbalance occurs at any point, TenneT will take measures to restore this within an ISP by determining imbalance price. It’s in the BRPs best interest to solve the imbalance, so it might worth for them to tune up, or down the electricity generation in exchange the imbalance price.
For power imbalances BRPs place bids — of the amount of electricity they want to sell or buy with a price for it — that results in a supply-demand curve called the merit order. At the equilibrium point the day-ahead electricity price is determined for the delivery day. See Figure 1 for an example. The delivery day when the state of the power imbalance happens the TSO is solving it by every minute calculate a new merit order. These orders include the deviations of the BRPs schedules, that materialize next day as settlement prices and the reserves of BSPs.
Table 3 shows who pays whom settlement price in different imbalance scenarios. Further details of imbalance cost can be found in TENNET (2019).
EUMETSAT spectral bands
Short description and general usage of the 11 spectral channels of EUMETSAT satellite imagery that were used in this thesis:
• VIS0.6, VIS0.8: These channels are in the human visible wavelengths. Used for cloud detection, cloud tracking, location identification, aerosol, land surface and vegetation monitoring.
• (N)IR1.6: Near Infra-red channel. Discriminates between snow, cloud, ice and water clouds, provides aerosol information.
• IR3.9: Primarily for low cloud and fog detection. Also supports measurement of land and sea surface temperature at night and increases the low level wind coverage from cloud tracking.
• WV6.2, WV7.3: These two channels are the water vapor channels for observing water vapor and winds. • IR8.7: The channel provides quantitative information on thin cirrus clouds and supports the
discrimi-nation between ice and water clouds.
• IR9.7: This channel contains ozone radiances. It’s possible to use it as input to an NWP. Wind motion in the lower stratosphere and the state of ozone field can be observed.
• IR10.8, IR12.0: Cirrus clouds, volcanic ash clouds can be detected in this band and temperatures could be inferred.
• IR13.4: The CO2 absorption channel. In cloud-free areas, it will contribute to temperature information from the lower troposphere.
HRV channel was not used due to the images being much higher resolution and causing errors with mem-ory constraints in the preprocessing phase. Further details on spectral channels can be found in JO-HANNES SCHMETZ (2002).
Figure 14 shows how the categorical variable, quarter of the day in a dependece plot. The x values are possible values (that translate from 8am to 4pm in 15 min intervals) and the y values show how the feature
Regulation/ Price Positive (+) Negative (-) Upward TSO→ BSP BSP→ TSO Downward BSP→ TSO TSO→ BSP
Table 3: The direction of imbalance settlement payments
1/3/6 hours Dense DenseRotation Sparse SparseSD Persistence IR108 .073/.129/.176 .073/.131/.179 .111/.172/.228 .112/.188/.252 .097/.156/.204 IR120 .07/.124/.169 .071/.126/.171 .107/.166/.221 .109/.182/.244 .094/.151/.197 IR134 .055/.097/.134 .055/.098/.134 .092/.145/.197 .093/.157/.213 .075/.12/.155 IR16 .157/.272/.313 .157/.278/.315 .172/.262/.42 .201/.334/.379 .178/.292/.323 IR39 .145/.245/.334 .145/.246/.338 .228/.322/.429 .221/.344/.452 .166/.266/.357 IR87 .077/.138/.188 .078/.14/.191 .119/.184/.244 .119/.2/.267 .101/.165/.216 IR97 .068/.119/.162 .068/.12/.164 .112/.173/.235 .111/.187/.259 .086/.14/.182 VIS6 .179/.317/.398 .18/.324/.401 .217/.339/.54 .235/.392/.473 .214/.354/.426 VIS8 .149/.271/.324 .149/.279/.326 .169/.274/.475 .201/.346/.405 .17/.292/.332 WV62 .046/.08/.124 .046/.079/.121 .097/.171/.242 .1/.176/.25 .062/.104/.148 WV73 .045/.079/.117 .045/.079/.115 .087/.147/.208 .088/.151/.214 .064/.102/.137
Table 4: Normalized MAE averages
1/3/6 hours Dense DenseRotation Sparse SparseSD Persistence IR108 .913/.755/.571 .912/.748/.556 .761/.484/.157 .74/.415/.059 .835/.635/.426 IR120 .914/.756/.572 .914/.751/.56 .754/.466/.132 .732/.396/.028 .835/.633/.424 IR134 .932/.8/.624 .932/.797/.622 .736/.415/.037 .7/.33/-.089 .858/.689/.504 IR16 .816/.511/.32 .814/.491/.313 .707/.406/-.106 .685/.275/.033 .754/.44/.272 IR39 .792/.504/.186 .791/.498/.165 .456/.111/-.345 .448/-.021/-.489 .726/.414/.086 IR87 .916/.765/.591 .916/.759/.575 .776/.515/.2 .76/.461/.136 .847/.659/.462 IR97 .949/.856/.739 .949/.852/.731 .828/.62/.342 .822/.587/.285 .91/.795/.668 VIS6 .859/.61/.356 .858/.59/.348 .741/.436/-.081 .759/.424/.135 .794/.517/.267 VIS8 .851/.588/.362 .85/.571/.36 .732/.411/-.169 .722/.345/.059 .797/.528/.328 WV62 .973/.907/.763 .973/.908/.774 .717/.324/-.143 .727/.357/-.127 .94/.832/.666 WV73 .956/.861/.708 .956/.862/.717 .714/.348/-.077 .708/.324/-.144 .901/.774/.616 Table 5: R2 averages
1/3/6 hours Dense DenseRotation Sparse SparseSD Persistence IR108 .057/.095/.125 .057/.096/.127 .091/.134/.175 .098/.146/.184 .078/.116/.145 IR120 .056/.094/.125 .057/.096/.127 .092/.136/.178 .099/.148/.187 .079/.116/.145 IR134 .05/.085/.115 .05/.086/.116 .094/.14/.183 .103/.152/.193 .072/.106/.133 IR16 .074/.119/.155 .074/.12/.155 .101/.143/.192 .098/.148/.185 .087/.129/.16 IR39 .062/.099/.125 .062/.1/.127 .096/.128/.163 .092/.136/.165 .073/.109/.133 IR87 .056/.093/.123 .056/.095/.125 .089/.13/.171 .095/.141/.178 .076/.113/.141 IR97 .044/.074/.098 .045/.075/.099 .077/.115/.154 .083/.124/.16 .059/.088/.11 VIS6 .065/.106/.139 .065/.107/.14 .091/.131/.175 .086/.131/.163 .08/.119/.149 VIS8 .063/.102/.132 .063/.103/.133 .089/.129/.179 .088/.133/.164 .075/.111/.137 WV62 .025/.045/.072 .025/.045/.07 .075/.118/.154 .077/.119/.157 .037/.061/.085 WV73 .036/.064/.093 .036/.064/.091 .089/.135/.176 .091/.138/.18 .055/.082/.107
Table 6: Normalized RMSE averages
xgb nwp xgb nwp now xgb now ard nwp ard nwp now ard now VIS6 0.270344 0.262484 0.417479 0.369097 0.365463 0.496132 VIS8 0.270387 0.265134 0.387808 0.369099 0.365606 0.482962 WV73 0.270264 0.267837 0.502214 0.369023 0.366347 0.557565 IR120 0.272767 0.264167 0.467942 0.369155 0.367498 0.553722 IR97 0.270933 0.264292 0.487953 0.369093 0.366500 0.557702 IR87 0.271907 0.264892 0.466992 0.369107 0.367269 0.550218 IR108 0.272408 0.269292 0.472425 0.371879 0.369931 0.554962 IR39 0.269482 0.268590 0.479673 0.367524 0.366682 0.554574 IR134 0.270188 0.264078 0.486242 0.368861 0.367299 0.555977 WV62 0.270344 0.266690 0.510012 0.369097 0.367486 0.555448 IR16 0.270344 0.270084 0.406220 0.369097 0.366255 0.481871
Table 7: XGBoost and Automatic Relevance Determination Regression model results with the different features from NWPs and nowcasts.
Figure 7: Sample solar radiation data over four days
importance changes with its value. The color codes are the importance of a nowcasted feature. This figure shows that in early afternoon quarter of the day has higher impact on the result, and in the morning when quarter of the day has low impact pixel from the nowcasts have more explanatory value. Later in the day this phenom disappears.
Figure 8: Sample solar radiation data over one and a half years
Figure 10: Power curve fitting, default is the manufacturer’s provided power curve, fitted is based on observations
