This paper applies standard regression techniques to forecast the country-wide consumption of electrical energy. All considered machine learning algorithms are available as a part of the Scikit-learn library. Besides the fine-tuning of regression hyperparameters, several data preparation techniques are employed to improve the forecasting accuracy. It is demonstrated that forecasting for 24 hours ahead is possible with good accuracy and has practical significance.
• • • the volume of the market decreases, the price of the electricity increases, the volatility of the volume increases.
The laws of physics apply to electrical circuits regardless of the scale. There are some electricity losses associated with resistance, but usually, they are negligible. If the amount of electrical energy traded and transmitted is measured on substations, we can conclude that the amount of produced electricity is exactly equal to the amount of consumed electricity. That is, for the purpose of this paper we can use the following terms interchangeably: energy production, energy consumption, and market volume. When the country is considered an open system, the following equation applies.
production + import = consumption + export. (1)
The dataset used in this research represents the time range from July 1st, 2020, to December 31st,
2021. For historical reasons, the time range from July 1st, 2019, to June 30th, 2020, did not contain
bilateral market data [
It is often the case the modeled process is affected by other external factors not represented via input parameters from the original dataset. The outside temperature influences the consumption of electricity as more energy is needed in winter for heating and in summer for air-conditioning. Two columns with hourly data were added to the dataset representing the temperature for Ukraine and its capital, see the dependencies below in Figure 3. The location representing the country was selected as its linear geographic center with decimal GPS coordinates 48.379433N 31.165580E.
Another important factor is the periodicity in the consumption of electrical energy. For example, at night people need less electricity than in the daytime. Similarly, on weekends the electricity consumption is lower than on weekdays. This paper considers four cycle types: daily, weekly, monthly, and yearly. One of the next sections will analyze whether these additions are helpful.
The problem is how to feed time representation to the machine learning algorithm in a way that
similar moments in time would be interpreted as close by the algorithm. As shown in Figure 4, hour
values 23 and 0 are close on the timescale, but they are distant in real-valued representation. One of
the possible solutions to this problem is to calculate the sine and cosine of the cycle phase [
composite primary key. In addition to the original 4 attribute columns with market volume data now we have 10 more columns. The temperature data were downloaded from the site https://openweathermap.org, the periodic columns were calculated using an algorithm written in
The usage of additional input parameters typically provides better regression results. If we need to forecast market volumes for 24 hours ahead then it makes sense to take into account the available data for the last 24 hours (at least). The machine learning algorithms and library functions expect that both input and output parameters are represented as one record. So, as a data preparation step, the data displayed in Figure 6 were resampled into the following columns, where M1 suffix means the parameter was taken one hour ago, P1 suffix means the parameter was taken one hour later, etc.
resampling were not fully qualified. The dataset was split into training and testing parts using the
standard library function train_test_split from sklearn.model_selection namespace [
To measure the influence of input parameters, we used the nearest neighbors regression model represented by class KNeighborsRegressor from sklearn.neighbors namespace. This machine learning algorithm provides quite competitive results and has a small number of hyperparameters to optimize.
The Python code snippets that implement this functionality are provided in Appendix 1. The complexity of the algorithm is hidden behind fit and prediction methods. Other regression and classification algorithms also reuse these methods, so the substitution of one algorithm instead of another is relatively simple.
The metrics used to measure the discrepancy between the test set and forecasted data are given in Table 1. Here yi is the output value from the i-th record in the testing dataset, fi is the predicted value for the i-th record, y is the average output value over the test dataset. These formulas are considered in the context of one selected output column representing the market volume.
Now we need to evaluate the effect of additional parameters and history length on prediction accuracy. Table 2 shows the accuracy improvements after adding temperature and periodic
R2 = 1 − MAPE = MAE = ∑i ( yi − fi )2 ∑i ( yi − y )2 1 n
∑ n i=1 yi − fi
yi 1 n
∑ yi − fi n i=1 (2) (3) (4) parameters. It appears all additional parameters are useful, but the overall effect is rather minor. Here are the parameters for the starting model.
Input columns: Bilateral, DayAhead, Intraday, Balancing Output columns: BilateralP24, DayAheadP24, IntradayP24, BalancingP24
Input columns: Bilateral, DayAhead, Intraday, Balancing, BilateralM1, DayAheadM1, IntradayM1, BalancingM1, ..., BilateralM23, DayAheadM23, IntradayM23, BalancingM23, TempUkr, TempKiev, SinDay, CosDay, SinWeek, CosWeek, SinMonth, CosMonth, SinYear, CosYear
The high dimensionality of input space is typically considered a problem, especially with noisy data. On the other hand, not all input parameters explored so far have equal contribution to the quality of results. So, it would be helpful to try removing the parameters that provide less useful information than others.
It appears this is not complex with the class SelectFromModel from sklearn.feature_selection
namespace [
So far, all the results were related to 24-hour forecasting. Figures 8 and 9 below show the R2 score and mean absolute percentage error for the range from 1 and up to 24 hours. The one-hour forecasting provides the best results. It is also worth noting that bilateral and day-ahead markets have much better predictability than balancing markets. As for the intraday market, it has a low mean absolute error just because the size of this market is small.
The 24-hour prediction error for all four markets can be measured on the test set, which represents 20 % of the original dataset. For convenience in representation and analysis, the test set was sorted by real market volume. The predicted values are shown in Figures 10–17 with dots. The probability distribution of error is shown using histograms. An interesting finding is that forecasting error is not always Gaussian. In particular, this is the case for bilateral and intraday market volumes.
The curve representing balancing market volume in Figure 16 crosses the zero line, also it has more negative values than positive. This can be interpreted as that market players tend to overbuy electricity in other markets, so they need to sell more on average at the last moment. Let’s note that this inefficiency can be mitigated by the usage of forecasting models.
So far, all the results were obtained with the nearest neighbor regressor. And it makes sense to
explore the performance of other algorithms on the same column configuration that is represented
in Appendix 3. The output parameters were selected for 24-hour forecasting. The results shown in
Table 4 and Table 5 include the comparison with classic instruments like multi-layer perceptron [
It is worth noting that some algorithms do not natively support multi-output configuration, so it was needed to use the class MultiOutputRegressor to overcome this problem and cover four electrical energy markets with one machine learning model.
Tables 4–6 represent the following characteristics obtained for different machine learning models: R2 score, mean absolute percentage error, and mean absolute error. It appears, that for this specific task, the ensemble methods are much better than others, and the winning algorithm Histogram Gradient Boosting is one of them. Also, it is one of the fastest and it can flawlessly handle datasets with missing values. On the current dataset, the training phase takes about 20 seconds.
Two different training approaches were used for multi-layer perceptron: quasi-Newton optimizer (QNO) and stochastic gradient descent (SGD). The first algorithm uses analytic solution to weight optimization problem, while the second algorithm employs an iterative process to find the minimum of error function. In both cases the architecture was the same, the perceptron had four layers (that is two hidden layers). Empirically, this architecture was more successful than three or five-layer perceptrons. Meanwhile, all these configurations are universal approximators.
It makes sense to explain the high mean absolute percentage error for the balancing column in Table 5. This is not the error, also a couple of zero values in the input dataset were replaced to improve MAPE figures. As shown in Figure 16, most of the values for this column are located close to zero. So, the calculations according to formula (3) involve the division by a small value. In other words, the MAPE metric is just not that adequate for the balancing column.
This paper demonstrates that proper forecasting model selection is a multi-stage process that may involve data selection, data preprocessing, data augmentation, selection of machine learning algorithm, optimization of hyperparameters, etc. While all computations for this work were done on a regular 8-core machine, the creation of the MLOps pipeline may require much more powerful computation resources.
The pre-trained model can be saved into a file for subsequent reuse in the production
environment. There are two formats popular among Python developers: .joblib and .pickle. In
addition, there is .onnx format that can be loaded not just in Python, but also in faster .NET or
Javabased applications [
There are dedicated solutions that can host the machine learning models using a microservice
approach like Seldon Core [
The forecasting accuracy that was obtained for electrical energy markets is different. Nevertheless, the 1 % error for 24-hour forecasting of the bilateral market looks impressive. Such a forecast can be useful at the country scale to ensure required fuel supply, plan import/export operations, reduce electricity costs, etc. Similar research can be done for more specific datasets from commercial and state energy enterprises.
The function that evaluates the input parameter set using nearest neighbors regressor: def evaluate_input_features \ (training_inputs, testing_inputs, training_outputs, testing_outputs): evaluation_regressor.fit(training_inputs, training_outputs) predicted_training_outputs = evaluation_regressor.predict(training_inputs) predicted_testing_outputs = evaluation_regressor.predict(testing_inputs) print_evaluation_metrics(training_outputs, testing_outputs, \ predicted_training_outputs, predicted_testing_outputs) def print_evaluation_metrics(training_outputs, testing_outputs, \ predicted_training_outputs, predicted_testing_outputs): training_score = r2_score(training_outputs, predicted_training_outputs, multioutput = 'raw_values', force_finite = True) print(F"Training dataset R2 score(s): {training_score}") training_percentage_error = mean_absolute_percentage_error \ (training_outputs, predicted_training_outputs, multioutput = 'raw_values') print(F"Training percentage error(s): {training_percentage_error}") train_standard_deviation = mean_squared_error \ (training_outputs, predicted_training_outputs, \ multioutput = "raw_values", squared = False) print("Training standard deviation(s):", train_standard_deviation) train_absolute_error = mean_absolute_error \
(training_outputs, predicted_training_outputs, multioutput = 'raw_values') print("Training mean absolute error(s):", train_absolute_error) test_score = r2_score(testing_outputs, predicted_testing_outputs, multioutput = 'raw_values', force_finite = True) print(F"Testing dataset R2 score(s): {test_score}") test_percentage_error = mean_absolute_percentage_error \
(testing_outputs, predicted_testing_outputs, multioutput = 'raw_values') print(F"Testing percentage error(s): {test_percentage_error}") test_standard_deviation = mean_squared_error \ (testing_outputs, predicted_testing_outputs, \ multioutput = "raw_values", squared = False) print("Testing standard deviation(s):", test_standard_deviation) test_absolute_error = mean_absolute_error \
(testing_outputs, predicted_testing_outputs, multioutput = 'raw_values') print("Testing mean absolute error(s):", test_absolute_error)
Identifying the features that provide higher information entropy: random_forest = RandomForestRegressor \
(n_estimators = 100, criterion = 'squared_error', ccp_alpha = 0.0) random_forest.fit(dataset_inputs, dataset_outputs) random_forest.feature_importances_ optimized_model = SelectFromModel(random_forest, \
threshold = "0.12 * mean", prefit = True) optimized_inputs = optimized_model.transform(dataset_inputs) Appendix 3
Most informative input parameters selected with a random forest model: selected_names = \ ['Bilateral', 'DayAhead', 'Intraday', 'Balancing', 'BilateralM1', 'DayAheadM1', 'IntradayM1', 'BalancingM1', 'DayAheadM2', 'IntradayM2', 'BalancingM2', 'IntradayM6', 'BilateralM7', 'DayAheadM7', 'IntradayM7', 'BalancingM7', 'BilateralM8', 'DayAheadM8', 'IntradayM8', 'BilateralM9', 'DayAheadM9', 'IntradayM9', 'BilateralM10', 'DayAheadM10', 'IntradayM10', 'BilateralM13', 'BilateralM14', 'IntradayM14', 'BalancingM14', 'BilateralM15', 'DayAheadM15', 'IntradayM15', 'BalancingM15', 'BilateralM16', 'DayAheadM16', 'IntradayM16', 'BilateralM17', 'DayAheadM17', 'IntradayM17', 'BilateralM18', 'DayAheadM18', 'IntradayM18', 'IntradayM19', 'BilateralM21', 'BilateralM22', 'DayAheadM22', 'IntradayM22', 'BilateralM23', 'DayAheadM23', 'IntradayM23', 'BalancingM23', 'TempUkr', 'TempKiev', 'CosDay', 'SinWeek', 'CosWeek', 'SinMonth', 'CosMonth', 'SinYear', 'CosYear']
MultiOutputRegressor(HistGradientBoostingRegressor (loss = 'squared_error', learning_rate = 0.20, max_iter = 300, early_stopping = False, scoring = 'loss', random_state = 1)) MultiOutputRegressor(AdaBoostRegressor (estimator = DecisionTreeRegressor(criterion = 'squared_error', splitter = 'best', max_depth = None, min_samples_split = 2, max_features = None, random_state = 1), n_estimators = 10, learning_rate = 1.0, loss = 'square', random_state = 1)) MultiOutputRegressor(GradientBoostingRegressor (loss = 'squared_error', learning_rate = 0.39, n_estimators = 100, subsample = 1.0, criterion = 'squared_error', max_depth = 6, random_state = 1, max_leaf_nodes = None, ccp_alpha = 0.0)) ExtraTreesRegressor(n_estimators = 100, criterion = 'squared_error', max_depth = None, max_features = 1.0, bootstrap = False, n_jobs = 8, random_state = 1, ccp_alpha = 0.0) KNeighborsRegressor(n_neighbors = 3, weights = 'distance',
algorithm = 'auto', p = 1, metric = 'minkowski', n_jobs = 8) RandomForestRegressor(n_estimators = 100, criterion = 'squared_error', max_features = 1.0, bootstrap = True, ccp_alpha = 0.0, n_jobs = 8, random_state = 1) MultiOutputRegressor(NuSVR(nu = 0.4, C = 1000000.0,
kernel = 'rbf', gamma = 'scale', shrinking = True, max_iter = -1)) MLPRegressor(hidden_layer_sizes = (200, 200,), activation = 'relu', solver = 'lbfgs', alpha = 0.0000, max_iter = 5000, random_state = 1) MLPRegressor(hidden_layer_sizes = (200, 200,), activation = 'relu', solver = 'adam', alpha = 0.0002, max_iter = 200, batch_size = min(50, training_set_size), shuffle = True, random_state = 1, early_stopping = False) ElasticNet(alpha = 1.0, l1_ratio = 1.0, fit_intercept = True, max_iter = 1000, positive = False, random_state = 1, selection = 'cyclic') LinearRegression(fit_intercept = True, n_jobs = 8) MultiOutputRegressor(BayesianRidge(max_iter = 300, tol = 0.001, alpha_init = None, lambda_init = 1.0, fit_intercept = True))