Prediction of Intraday Electricity Supply Curves

Abstract

The electricity market in Spain, as in many European countries, is organized into daily, intraday, and reserve markets. This project aims to predict the supply curves in the Spanish intraday market that have six sessions with different horizons of application, using information from the market itself. To achieve this, we approximate these curves using a non-uniform grid of points and evaluate the quality of these approximations with a weighted distance, both based on empirical market data. We employ neural network models, including multilayer perceptrons (MLPs), convolutional neural networks (CNNs), long short-term memory (LSTM), bidirectional LSTM (BiLSTM), and a Transformer network alongside a naive model for benchmarking. The MLP and CNN models demonstrated significant improvements in predicting these supply curves for the six market sessions.

1. Introduction

Having reliable predictions of future electricity transaction volumes and prices is valuable information for both generators and consumers, as well as for the system operator. For the generators and consumers, this data provides an advantage when planning and optimizing their processes and finances. For system operators, it allows them to anticipate issues related to the physical reality of electricity (such as potential overloads) and prevent them. One potentially accurate approach is to forecast these prices indirectly by intersecting supply and demand curves. As we will see in the next subsection, the electricity market is organized into different submarkets that occur at different times of the day, the offers of which have different application horizons. In each of these submarkets, the supply curves are obtained from the production proposals (energy and price) made by the participating companies.

The objective of this paper is to obtain a reliable prediction of the supply curves of the intraday market, which is a less studied market than the daily market that has a single session. The intraday market has six sessions during the period analyzed. In the next section, we will describe the Spanish electricity market structure, and we will show the motivations for the approach that we followed in predicting the supply curves for the intraday market.

Spanish Electricity Market

The enactment of Law 54/1997 on 27 November 1997 initiated a progressive liberalization of the Spanish electricity sector [1]. This reform involved allowing third-party access to networks, creating a structured energy trading market, and minimizing government control over system management. Currently [2], electricity generation and commercialization are liberalized activities in the market, while high-voltage transmission (which remains a monopoly controlled by the government) and low-voltage distribution are regulated activities. Figure 1 shows a representation of the structure of the Spanish electricity sector (see [3,4]) (In Figure 1: Audax Energía, Comisión Nacional de los Mercados y de la Competencia, Endesa, EDP HC Energía, Iberdrola, Gas Natural Fenosa, Nexus, Red Eléctrica de España and Viesgo are official names of companies and regulatory institutions).

The Spanish electricity market is actually a set of interconnected markets designed to meet the changing needs of electricity supply and demand. In addition to electricity transactions planned for weeks, months, or years ahead, many trades are conducted with only hours or even minutes of notice before the actual generation of the energy. Given that electrical energy is difficult to store, production at each moment must match the demand. Supply and demand have a volatile nature due to weather conditions, socioeconomic factors, and other variables, highlighting the need to constantly review negotiations in line with current forecasts. Additionally, the energy transmission network must be physically capable of executing what has been agreed upon.

Following the 2007 integration of the Spanish electricity market with the Portuguese one in the so-called MIBEL (Iberian Electricity Market), the resulting market was coupled with Europe progressively from 2013, establishing a single price in the connected regions when the system allows it.

In the daily market [6], which is organized by OMIE (Iberian Energy Market Operator), energy transactions between generators and consumers are conducted for each of the 24 h of the following day (which is referred to as the ‘dispatch day’ and is denoted by D). Once the bids from producers for each hour of the next day are submitted, they are arranged in ascending order, forming the market supply curve for each hour. The bids from buyers are ordered in descending order of price, forming the market demand curve. Finally, the market price for the following day is set by the intersection point of the supply and demand curves. It follows a pay-as-clear system, meaning that all participants, both buyers and sellers, pay the market price established for that period.

Following the daily market, participants have the opportunity to buy and sell electricity again in the intraday market, held a few hours before real-time operations, or in the continuous intraday market, which operates under a European crossborder framework [7]. Intraday markets, also organized by OMIE, provide a crucial mechanism for market players to fine-tune their schedules from the daily market to better align with real-time energy demands.

Lastly, there are reserve markets through which the system operator, Red Eléctrica de España, ensures the physical feasibility of the previously negotiated generation and distribution of energy, making adjustments to the agreements.

In this paper, we focus our attention on the intraday markets, with our goal being to forecast the supply curves of these markets, specifically using information from the market itself and its predecessor. In the period of interest for this project (see Section 3.1), these markets were organized into six auction sessions, which are carried out sequentially, each with its own application horizon.

Table 1. Schedules and horizons of intraday markets. Information available at: https://www.omies.es [7], accessed on 6 September 2024.

	1st Session	2nd Session	3rd Session	4th Session	5th Session	6th Session
Market opening	14:00 D-1	17:00 D-1	21:00 D-1	1:00 D	4:00 D	9:00 D
Market close	15:00 D-1	17:50 D-1	21:50 D-1	1:50 D	4:50 D	9:50 D
Matching	15:00 D-1	17:50 D-1	21:50 D-1	1:50 D	4:50 D	9:50 D
PIBCA program publication	15:07 D-1	17:57 D-1	21:57 D-1	1:57 D	4:57 D	9:57 D
PHF publication of the OSs	16:20 D-1	18:20 D-1	22:20 D-1	2:20 D	5:20 D	10:20 D
Market horizon	24 h	28 h	24 h	20 h	17 h	12 h
(Time period)	(1–24 D)	(21–24 D-1 & 1–24 D)	(1–24 D)	(5–24 D)	(8–24 D)	(13–24 D)

shows the schedules and market horizons for these sessions.

In Table 1, we can see that the times at which the programs and results of a session are published allow that information to be used to predict the curve of the following session. For example, the results and trading programs of the first session are published at 16:20 at the latest on day D-1, and the second session occurs between 17:00 and 17:50 on the same day. Similarly, the results of the daily market are known before the first session of the intraday market, and, therefore, their information can be used to predict the supply curve of the first session. This is the predictive modeling proposal that we will use in this project, using data from the previous market session to predict the next session. We, therefore, pose the following research questions that will be addressed in this project:

• Q1: Is it possible to predict supply curves more accurately than using the curves from the previous day?
• Q2: Is it useful to incorporate information from the previous market to predict curves from intraday market sessions?

The prediction of supply curves in the electricity market is not a new topic in the literature, but, as we will see in the next section, previous research has focused on the prediction of the day-ahead market or the reserve market, which are markets that occur once a day. Our research aims to close the gap in terms of a type of market structure in several daily sessions.

The rest of the paper is organized into four sections. Section 2 presents the state of the art regarding the prediction of supply curves. Section 3 develops the proposed methodology and describes the considered models. The results are presented in Section 4. Finally, the conclusions and possible extensions are discussed in Section 5.

2. State of the Art for Supply Curve Prediction

Most of the published works in this field focus on directly predicting the energy price and demand or other scalar variables. Three recent reviews on the techniques used for such purposes can be found in [8,9,10].

Our objective is to predict an offer curve, which is a more complex task. Nonetheless, obtaining this prediction is, indirectly, a prediction of the energy price, which is achieved by obtaining the intersection point of the forecasts of the supply and demand curves, as proposed by [11]. The two main challenges in predicting a supply curve are that it is not a prediction of a scalar quantity and it is not a continuous function.

An approach that addresses the first challenge is the procedure proposed by [12], which uses functional data nonparametric techniques to model residual demand curves. Since they use nonparametric kernel-type estimators, they are limited by the problem known as the curse of dimension and, therefore, only consider a single functional regressor. These authors assume that the curves are sufficiently smooth with up to two derivatives, which contradicts the second challenge.

In [11], LASSO regression is used to forecast a short number of bid classes obtained from the supply and demand curves in the EPEX market. This paper used the term “X-model” to denote the process of obtaining the price as a result of crossing the supply and demand curves. This is the first work that addresses both challenges since it takes into account the fact that the functions to be predicted are stepwise. In [13], a modification of the X-model is proposed, where transformed versions of the curves are used by assuming a perfectly inelastic demand. Both methods use a small number (only 16) of bid classes, which makes their approximation of the curves imprecise. In fact, as the motivation of [11,13] is price prediction, the authors do not study the behavior of the procedure as a curve predictor.

Parametric and nonparametric functional autoregressive models were used in [14] to model the supply and demand curves of the Italian electricity market. Again, the techniques used assume smoothness conditions that are not satisfied by step functions.

Another approach is to approximate the curves using others from a parametric family and predict the relevant parameters. In the literature, various procedures of this kind can be found, such as using radial basis functions in [15] or using polynomials and sigmoid functions [16]. In [15], univariate autoregressive models were used for each of the parameters, and weekly seasonality was modeled using dummy variables. The authors of [16] used vector autoregressive models with seasonal harmonic component periods of 24 and 168. This approach has the attraction that the number of parameters to be predicted is small, but it does not solve the challenge regarding the non-continuity of the curves.

In [17], a continuous piecewise version of the supply curve is used to approximate it with great precision. This approach has the disadvantage that it increases the number of nodes to approximate the curves.

The authors of [18] followed the approach of using a grid of points to generate stepwise functions for the supply curves at the Spanish secondary market, which is one of the reserve markets managed by Red Eléctrica de España. The authors consider random forest, gradient bootsting, dense neural networks, and long short-term memory models to predict supply curves assumed to be observed on a fixed grid with 50 nodes.

In [19], a distance learning procedure is proposed to forecast the day-ahead curves for the daily market. This approach does not need to approximate the curves because it relies on the distances between them. It proposes the prediction of the distances between the curves on day D+1 and the curves in the training set based on the distances between the curves on day D and those curves in the training set. The procedure is highly computationally intensive and it is not clear how other scalar covariates will be incorporated.

All the papers mentioned focus on daily markets, with the exception of [18], which studies a secondary reserve market. None of them have studied the intraday market with its respective sessions. It should be noted that this structure of the electricity market, which is subdivided into daily markets, reserve markets, and intraday markets, is common in many European countries (see [20,21]). Prediction in the intraday market has the added difficulty of several sessions being carried out on the same day and different prediction horizons, as shown in Table 1. In this project, we turn this difficulty into an advantage because we can use the information from the previous session to predict the session that follows. As we can see in Table 1, the publication times of the session results make this approach feasible.

3. Methodology

We now present the methodology followed in this paper in detail, dividing it into four parts: the source of the raw data and its format, the preprocessing performed once the data was obtained, the error metrics and the subsequent construction of the models, and the generation of the predictions. The Python 3.1 code implemented for this article is available upon request to the authors.

3.1. Raw Data

The raw data were obtained from the OMIE website, https://www.omie.es, accessed on 25 April 2024. The information is organized in files with a .1 extension (which can be read as .csv). Each file corresponds to the data of a specific market on a given day. The naming convention for these files works as follows:

• For the daily market, the names follow the format curva_pbc_uof_yyyymmdd.1, where yyyy refers to the year, mm to the month, and dd to the day the file refers to (to the date being negotiated, not the date the negotiation took place).
• For the intraday markets, the structure is similar, curva_pibc_uof_yyyymmddxx.1, where xx refers to the number of the market, going from 01 to 06. For instance, the file curva_pibc_uof_2021032704.1 contains data related to the fourth intraday market on 27 March 2021.

Each file starts with a block of metadata, followed by the relevant information. The information is presented in

Table 2. Example of raw data format. Source: OMIE, accessed on 25 April 2024.

Hour	Date	Country	Unit	Offer Type	Buy/Sell Energy	Buy/Sell Price	Offered/ Matched
1	01/01/2020	MI	ENDEC04	C	22.0	180.30	O
1	01/01/2020	MI	ENDE01	C	1237.7	180.30	O
1	01/01/2020	MI	EE21C01	C	1434.7	180.30	O

.

• Hour. This refers to the hour for which offers and demands are made. The number represents the position of that hour in the day, regardless of what the clock shows. For example, on a day when the clock is moved forward, the period from 03:00 to 04:00 corresponds to the third hour of the day, whereas, on a normal day, it would correspond to the fourth hour of the day.
• Date. This refers to the day for which offers and demands are made. It is given in the dd/mm/yyyy format.
• Country. This refers to the market to which the offer to buy or sell is presented. The codes are as follows: MI (Iberian Market: both Portugal and Spain), ES (Spain), or PT (Portugal).
• Unit. An alphanumeric code that is used to identify the current offer or demand.
• Offer type. This refers to whether the offer is for buying (C) or selling (V).
• Buy/Sell Energy. The amount of energy offered or demanded, expressed in MWh and in the Spanish number format.
• Buy/Sell Price. The bid or ask price for the respective energy, expressed in EUR/MWh and in the Spanish number format.
• Offered/Matched. This indicates whether the energy–price pair corresponds to the offers curve (O) or the resulting matching program (C).

• Choice of Date Range

For the development of this work, we used only the information from 2022 and 2023. The reason for using only these two years (despite having data from a longer period) is due to changes in market regulations. In particular, the price ranges of the markets have varied substantially over the years. Although we normalized the prices to standardize the limits, the density of their distribution varies depending on the original price range, and this affects the models’ performance. The two most recent changes are

• 20 May 2021 [22]. Maximum and minimum offer price limits were set at −500 EUR/MWh and 3000 EUR/MWh for the daily market and −9999 EUR/MWh and 9999 EUR/MWh for the intraday market.
• 11 May 2022 [23]. The price range from −500 EUR/MWh to 3000 EUR/MWh of the daily market was changed to the new one from −500 EUR/MWh to 4000 EUR/MWh. Despite this change, there are few offers above the previous limit of 3000 EUR/MWh. Therefore, we decided, instead of normalizing, to divide by the previous legal range (which remains, in practice, as the current range).

3.2. Preprocessing

3.2.1. Train-Test Split

The data were split into training and testing sets, with the training set covering the dates from 1 January 2022 to 31 December 2022 (both dates inclusive) and the testing set covering the remaining dates. Regarding the partition of the training set into train-train and train-validation, the first one covers from 1 January 2022 to 31 October 2022 (both dates inclusive), while the second one covers the rest of 2022.

3.2.2. Data Preprocessing

The first step taken before preprocessing was to verify that all the files in our interest existed and that these files were in the correct format (for example, not empty); we did not find any errors.

Once that was carried out, the following process was carried out for each of the files. Once the file was read as a .csv file, we filtered the data to retain only the sell offers from the original offer curve (not the resulting matching program). Then, the numerical values were changed from the Spanish number format to the English one. Moreover, the numerical data were formatted into Python’s float type, and the names of some variables were shortened.

Then, the offers were grouped by day and hour, grouping them again into larger 24-h blocks (day by day). However, this poses a problem. In Spain, the clocks are set forward 1 h in the spring and set back 1 h in the fall [24]. Therefore, the affected dates have a different number of hours, which may conflict with our desired homogeneous data format. The specifics are as follows:

• The markets with a 24 h application horizon (the daily market and the first and third intraday markets) are affected by these time changes. Some days have 23 h, while others have 25. Since these represent a tiny percentage of the total dates, we chose to remove the last hour on 25 h days and duplicate the last hour on 23 h days, ensuring we always had 24 h.
• The second daily market, with its 28 h application period spread over 2 days, required a more detailed examination. Its application horizon covers an entire day and the last 4 h of the previous day. This market is only affected by time changes when they impact the 24 h day it applies to, resulting in markets with either 27 or 29 sessions. We handled this full day in the same way as the 24 h horizon markets.
• The offers and demands in the intraday markets with a horizon of less than 24 h (the fourth, fifth, and sixth) pertain to the respective last hours of the day, and therefore, these markets are not affected by the time changes. However, on 10 October 2022 (the clocks were set back on this date [24]), the fourth intraday market conducted 21 h instead of 20. Given the minor nature of the issue, we addressed it by removing the first hour and keeping the last 20, which is in line with the market structure and the format of our data.

Another potential issue is that, due to anomalous situations, technical complications, or regulatory changes, a day might have an abnormal number of sessions. Fortunately, there were no such problems in our date range.

Once the issues related to temporal homogeneity had been assessed, each hourly block was sorted in ascending order by offer price. Next, we normalized the price values and summed the energies corresponding to these prices cumulatively, replacing the original energy value and finally obtaining the offers curve. We merged the hourly blocks into a larger dataframe of 24 h and attached this block to all the previous ones, creating a single dataframe for all prices and another for all energies (logically, this process is carried out once for each market).

It should be noticed that the dimensions of these dataframes are not homogeneous since they vary depending on how many different prices are present in the corresponding session. Structural homogeneity in our data is addressed in Section 3.2.4.

3.2.3. Representation of Curves

Regarding the energy-price relationship, we can represent it using either the P-Q representation (with prices on the x-axis and energy on the y-axis) or the Q-P representation (where the axes are reversed). We opted for the first option because, in this way, the curves always have the same domain (the range of prices is the same). This approach provided a more uniform framework.

3.2.4. Approximation Procedures for Curves

Since our goal was to predict a curve, we needed a method to synthesize its information and work with this synthesis as input and output for our models. Two straightforward approaches are to approximate it with curves from a parametric family (and work with these parameters) or to approximate it using a grid of points (and work with them), interpolating the intermediate values. We chose the second option since the use of a parametric family usually requires assumptions about the smoothness of the functions, and we know that the supply functions are step functions.

There are three issues to consider:

• Interpolation technique. Since we were working with nondecreasing step functions, it was reasonable to interpolate consistently with this property, using right-side constant interpolation (i.e., the interpolated values are equal to the largest grid point value smaller than them).
• Uniformity of the grid. Despite the broad price ranges presented by the markets, in reality, most offers occur within a more limited range. Given our interest in accurately representing this subinterval, it made more sense to use a non-uniform grid, allocating more density to these empirically observed regions. In Figure 2 and Figure 3, we observe that, although both interpolations generally fit the original curve well, the non-uniform interpolation captures the section of the curve with a higher density of transactions much better.
• Grid size. There is a trade-off between precise representation and the informational weight of the grid (which we wanted to reduce compared to the raw data). We opted for a grid of 150 points.

3.3. Error Metrics

To define an approximation error, we needed to establish when two curves are considered similar. We took the Euclidean distance, $d_2$, which is commonly used in functional data analysis, which is defined as

$$d_2(f,g)=\sqrt{\int {(f-g)}^2},$$

(1)

where f and g are curves. The $d_2$ distance in (1) measures how that separation varies across the entire domain. Although it is generally a good tool, in our case, it presented two problems.

• Integrals are a continuous operation, whereas we only had a grid of points. The accuracy of our approximation to the true value of the integral would be better if we used a more dense grid of points. Therefore, we calculate the approximation

  $$\sqrt{\sum _x{\left(f\left(x\right)-g\left(x\right)\right)}^2\Delta x}$$

  (2)

  by not using the non-uniform grid previously discussed but by reconstructing a more dense curve. The value of $\Delta x$ in (2) was set to 0.0001. Note that by scaling the prices we have, the integral is performed on the interval [0, 1].
• As mentioned, almost all offers and demands occur within a price range that does not cover everything considered by the regulations. See, for instance, Figure 2 and Figure 3, where the scaled prices are concentrated in the interval [0.4, 0.6]. When using the interpolation we just mentioned, the non-uniformity introduced in the grid is lost. To bring it back, we modified the $d_2$ distance formula. The idea is to use a density factor, K, that gives more weight to the more typical price ranges:

  $$d_K(f,g)=\sqrt{\int {(f-g)}^{2}K}.$$

  (3)
  The K factor in (3) is obtained as an estimate of the density of the bid prices for each session. This means that if, in a session, the majority of bids have their (scaled) prices in the interval [0.4, 0.6], then the differences in that interval have a greater weight than outside that interval. Each session has or can have a different K factor since it depends on the prices of its offers.

Given that we had a huge number of predictions, we calculated statistical metrics to summarize the errors made: the mean, median, standard deviation (SD), and mean absolute deviation (MAD).

3.4. Models and Predictions

3.4.1. Input and Output Format

We opted for a direct multi-output approach. Particularly, for predicting a specific market, we incorporated all data from the preceding market on the same day as well as all data from the previous day of the target market. For the first intraday market, we used the daily market as its predecessor. Therefore, the inputs had a dimension of $(H_s+H_p)\times 150$, where $H_s$ is the number of horizon hours of the market we were trying to predict, and $H_p$ is the length of the horizon for the preceding market. Consequently, the output had dimensions of $H_s\times 150$.

3.4.2. Correction of Predictions

The neural networks we used to predict the curves, which are reconstructed from the predicted grid values using the interpolation discussed in Section 3.2.4, successfully captured their shape. However, they failed to generate curves satisfying the assumption of being nondecreasing. To address this, we opted for the following procedure to produce a corrected curve:

$$t\left(x\right)=\left\{\begin{array}{cc}c\left(x\right)\hfill & \mathrm{if}c\left(x\right)>t\left(y\right)\mathrm{for}\mathrm{all}y<x,\hfill \\ t\left(z\right)\hfill & \mathrm{if}\mathrm{there}\mathrm{exists}z<x\mathrm{such}\mathrm{that}t\left(z\right)\ge t\left(y\right)\mathrm{for}\mathrm{all}y<x\mathrm{and}t\left(z\right)\ge c\left(x\right),\hfill \end{array}\right.$$

(4)

where c is the predicted curve and t is the corrected curve. An example of the subtle but crucial effect of the transformation defined in (4) can be observed in Figure 4.

3.4.3. Selection of Hyperparameters

All the models (except for the naive one) have hyperparameters to be tuned. To do this, we used the partitioning of the training set into the train-train and train-validation sets mentioned earlier. The hyperparameter selection that performed best in this partition was used for the final prediction in the train-test partition.

Given the scale of our models, testing a large number of hyperparameters entailed an unmanageable computational cost. Therefore, we opted for the following approach: for models with 300 or fewer hyperparameter combinations, we tested all of them. If the number of possibilities exceeded this value, we randomly selected 300 combinations (always using the same seed). In all models, this value never accounted for less than 40% of the total possible combinations, providing greater confidence that the selected set is close to the optimal in terms of performance.

3.4.4. Model Comparison

Finally, we needed a tool that allowed us to determine whether one model was statistically significantly better than another. We used the Diebold-Mariano test [25] (see also [26]), which compares whether the two forecasts ${\left\{a_i\right\}}_{i=1}^n$ and ${\left\{b_i\right\}}_{i=1}^n$ of a series of real numbers, ${\left\{x_i\right\}}_{i=1}^n$, are significantly different. The null and alternative hypotheses of this test are

• $H_0$. Both forecasts have the same accuracy, meaning that the observed differences are likely due to randomness;
• $H_1$. One forecast has better accuracy than the other, meaning that the predictive difference between the models is unlikely to be due to randomness.

To calculate the statistics, we obtain the series of the differences in losses.

$$d_i={(x_i-a_i)}^2-{(x_i-b_i)}^2,$$

(5)

as well as their mean and autocovariance, which are, respectively,

$$d={\displaystyle \frac{1}{n}}\sum _{i=1}^n{d}_i\mathrm{and}{\gamma }_k={\displaystyle \frac{1}{n}}\sum _{i=k+1}^n(d_i-d)(d_{i-k}-d).$$

(6)

Finally, the value of the statistic is obtained from the expression

$$\mathrm{DM}={\displaystyle \frac{d}{\sqrt{\frac{1}{n}({\gamma }_0+2{\sum }_{k=1}^{h-1}{\gamma }_k})}},$$

(7)

where it is generally adequate to use the value $h=\sqrt[3]{n}+1$. Under the null hypothesis, the DM statistic distributes as an $\mathcal{N}(0,1)$. In the case of rejecting the null hypothesis, it follows that the first model has better predictions if DM has a negative value; if DM has a positive value, then the second model is the one with the best performance.

The problem we faced was that, as we have seen, the Diebold-Mariano test is designed to compare two time series of real numbers as predictions of another real number series, but the objects in our time series are curves. We addressed this issue by considering not the time series of curves but the sequence of errors (which are real numbers), meaning we substitute $x_i-a_i$ and $x_i-b_i$ in (5) with the prediction errors of the two curves’ prediction methods. That is, instead of $x_i-a_i$ and $x_i-b_i$, we use $d_k(x_i,a_i)$ and $d_k(x_i,b_i)$, respectively. We have used the db.test function included in the forecast package in R 4.2.2 [27].

3.4.5. Naive Model

The naive model assumes that the offer curve for a given day is exactly the same as that of the previous day. Obviously, this model is computationally inexpensive and served as a benchmark to determine whether the other models were significantly better.

3.4.6. Multilayer Perceptron Model

The multilayer perceptron [28] (MLP) is the simplest neural network. For an application of MLP in energy forecasting, see [29]. It consists of several layers of neurons (computational unit that processes inputs to produce an output, contributing to the network’s overall task), with each neuron being fully connected to all neurons in the previous and subsequent layers. By adjusting the weights and biases (parameters) using a back-propagation algorithm and a cost function, the network learns the patterns in the training data. The considered hyperparameters were the following:

• epoch. The number of complete passes of the training dataset (epochs). The assessed values were 5, 10, and 20;
• n. The number of layers besides the input and output layers. We considered one, two, three, and four layers;
• neuron. The number of neurons in each layer. The tested values were 5, 10, 20, 50, 100, 200, and 400;
• act. The activation function, which is a mathematical function applied to each node, introduces nonlinearity into the network and allows it to learn complex patterns. We tried two different functions: ReLU and tanh;
• batch. The number of training examples used in one forward and backward pass of the training process. The tested values were 3, 10, and 30.

The selected hyperparameters for the MLP model are available in

Table A1. Selected hyperparameters for the MLP model.

Market	best_mean	best_epoch	best_n	best_neuron	best_act	best_batch
1st	153,721.7	20	3	200	relu	3
2nd	121,656.5	20	3	400	relu	3
3rd	107,365.2	20	3	200	relu	3
4th	97,758.4	20	1	100	relu	3
5th	66,961.5	20	4	100	relu	3
6th	67,873.2	20	3	100	relu	3

in Appendix A.

3.4.7. CNN Model

Convolutional neural network (CNN) [30] models are primarily used for visual data analysis (for example, image classification or object detection). For an application of a CNN in energy demand forecasting, see [31]. The main features of these networks are convolutional layers (where convolutional filters are applied, producing feature maps), pooling layers (which reduce the spatial dimensions of the feature maps), and fully connected layers (used to make final predictions). The considered hyperparameters were the following:

• epoch. The number of complete passes of the training dataset (epochs). The assessed values are 50 and 100;
• n_dense. The number of fully connected layers besides the output layer. We studied one, two, and three layers;
• n_conv. The number of convolutional and pooling layers (we always introduce a pooling layer after each convolutional layer). We considered one, two, and three layers;
• neuron. The number of neurons in each fully connected layer. The tested values are 5, 10, 20, 50, and 100;
• fil. The number of filters in each convolutional layer. We considered eight and 16 filter cases;
• kernel. The dimensions of the kernel matrix in each convolutional layer. The values studied are $3\times 3$, $5\times 5$, and $7\times 7$.

The selected hyperparameters for the CNN model are available in

Table A2. Selected hyperparameters for the CNN model.

Market	best_mean	best_epoch	best_n_dense	best_n_conv	best_fil	best_kernel	best_neuron
1st	159,051	50	2	1	16	7	50
2nd	117,692	50	2	1	8	5	100
3rd	95,489.6	100	1	3	8	5	100
4th	89,170.1	100	1	3	8	3	100
5th	58,594.2	100	2	1	8	5	100
6th	65,483.3	50	3	3	8	7	50

in Appendix A.

3.4.8. LSTM Model

Long short-term memory (LSTM) models [32] are recurrent neural networks used for time-series data and natural language processing. For an application of LSTM in energy consumption forecasting, see [33]. This model successfully addresses the vanishing or exploding gradient problem that other recurrent networks face, thus understanding long-range dependencies in sequence data. Instead of regular nodes, such as in MLP, LSTM networks have cells (arranged in layers), which are more complex structures that can update the information they store by using mechanisms such as the forget gate or the cell state. In our model, we also added fully connected layers after the LSTM layers. The considered hyperparameters were the following:

• n_LSTM. The number of LSTM layers. We tested one and two layers;
• n_dense. The number of fully connected layers besides the output layer. The values studied are one, two, three, and four;
• unit. The number of cells per LSTM layer. The assessed values are 5, 10, 15, 20, and 50;
• neuron. The number of neurons in each fully connected layer. We tried 10, 40, 50, and 100 neurons;
• dropout. The random deactivation of some input units to reduce overfitting. The studied values are 0, 0.1, and 0.2.

The selected hyperparameters for the LSTM model are available in

Table A3. Selected hyperparameters for the LSTM model.

Market	best_mean	best_n_LSTM	best_n_dense	best_unit	best_neuron	best_dropout
1st	163,648	1	4	20	40	0.2
2nd	147,153	1	2	10	100	0.2
3rd	133,267	1	1	20	50	0
4th	125,865	1	3	10	100	0
5th	102,703	1	1	100	100	0
6th	101,212	1	2	100	40	0

in Appendix A.

3.4.9. BiLSTM Model

Bidirectional long short-term memory (BiLSTM) [32] networks consist of two LSTM networks, processing sequential data both forwards and backward. This allows the model to capture information from both past and future contexts in a sequence. For an application of BiLSTM in power consumption prediction, see [34]. Once again, this structure is followed by fully connected layers in our model design. The considered hyperparameters were the following:

• n_LSTM. The number of LSTM layers. We tested one and two layers;
• n_dense. The number of fully connected layers besides the output layer. The considered values are one, two, three, and four;
• unit. The number of cells per LSTM layer. We tried 10, 20, 50 and 100 units;
• neuron. The number of neurons in each fully connected layer. The assessed values are 10, 40, 50, and 100;
• dropout. The random deactivation of some input units to reduce overfitting. The studied values are 0, 0.1, and 0.2.

The selected hyperparameters for the BiLSTM model are available in

Table A4. Selected hyperparameters for the BiLSTM model.

Market	best_mean	best_n_LSTM	best_n_dense	best_unit	best_neuron	best_dropout
1st	168,036	1	4	50	100	0
2nd	137,046	1	3	10	40	0.1
3rd	134,913	1	3	20	40	0
4th	120,819	1	2	10	100	0.1
5th	93,953.1	1	4	50	100	0
6th	93,598.4	1	4	10	100	0.1

in Appendix A.

3.4.10. CNN+LSTM Model

This model consists of CNN layers, followed by LSTM layers, and finally, fully connected layers. The considered hyperparameters were the following:

• n_CNN. The number of convolutional and pooling layers (we always introduce a pooling layer after each convolutional one). We considered one and two layers;
• n_LSTM. The number of LSTM layers. The considered values are one and two;
• n_dense. The Number of fully connected layers besides the output layer. We tried one, two, three, and four layers;
• filter. The number of filters in each convolutional layer. We considered eight and 16 filters;
• kernel. The dimensions of the kernel matrix in each convolutional layer. The assessed values are $3\times 3$, $5\times 5$, and $7\times 7$;
• unit. Number of cells per LSTM layer. We tried 10, 20, 50, and 100 units;
• dropout. The random deactivation of some input units to reduce overfitting. The studied values are 0, 0.1, and 0.2 (ranging from 0 to 1);
• neurons. The number of neurons in each fully connected layer: 5, 10, 20, 50, and 100 neurons were considered.

The selected hyperparameters for the CNN+LSTM model are available in

Table A5. Selected hyperparameters for the CNN+LSTM model.

Market	best_mean	best_n_CNN	best_n_LSTM	best_n_dense	best_filter	best_kernel	best_unit	best_dropout	best_neuron
1st	214,795	2	1	2	8	3	20	0.2	50
2nd	194,554	2	2	4	8	5	100	0	5
3rd	187,181	2	1	4	8	7	20	0.1	10
4th	164,056	2	2	4	8	3	20	0.2	10
5th	138,093	1	1	3	16	5	20	0.2	20
6th	125,495	1	1	4	16	3	10	0	10

in Appendix A.

3.4.11. Transformer Model

Transformers are a relatively new neural network architecture [35] and have revolutionized natural language processing. The Transformer model was used for electricity load forecasting in [36]. Unlike previous frameworks (such as recurrent neural networks), Transformers do not process information sequentially. Instead, they use a self-attention mechanism that allows them to capture the relationships between words in a sequence regardless of their position. The considered hyperparameters were the following:

• epoch. The number of complete passes of the training dataset (epochs). The assessed values are 20 and 50;
• n_dense. The number of fully connected layers besides the output layer. One, two, and three layers were considered;
• n_trans. The number of Transformer blocks. We tried one, two, and three layers;
• head. The number of parallel attention mechanisms. The tested values are eight and 16;
• dim. The size of the query, key, and value vectors used in the attention mechanism. The assessed values are 10, 20, 50, and 75;
• dropout. The random deactivation of some input units to reduce overfitting. The studied values are 0, 0.1, and 0.2;
• neuron. The number of neurons in each fully connected layer. We considered 20, 50, 100, and 200 neurons.

The selected hyperparameters for this model are available in

Table A6. Selected hyperparameters for the Transformer model.

Market	best_mean	best_epoch	best_n_dense	best_n_trans	best_head	best_dim	best_dropout	best_neuron
1st	210,978	50	3	1	8	10	0.1	200
2nd	196,384	50	3	3	8	20	0	200
3rd	185,167	50	3	3	8	10	0.2	200
4th	163,819	50	3	2	8	50	0	200
5th	139,121	50	3	3	8	50	0.1	200
6th	126,295	50	3	1	8	10	0	200

in Appendix A.

4. Results

4.1. Error Metrics

We now present the error metrics for the seven considered models.

Table 3. Summary statistics of the prediction error with the naive model.

Market	Mean	Median	SD	MAD
1st	161,035	119,650	144,772	57,635
2nd	201,648	153,984	170,172	82,586
3rd	137,397	101,724	128,090	50,949
4th	190,072	149,545	155,449	77,801
5th	198,397	157,803	158,881	84,706
6th	117,832	88,306	101,669	44,601

,

Table 4. Summary statistics of the prediction error with the MLP model.

Market	Mean	Median	SD	MAD
1st	182,371	142,978	140,956	67,747
2nd	126,912	101,131	91,371	42,882
3rd	106,281	85,497	74,003	31,410
4th	99,136	79,452	69,457	30,375
5th	71,913	57,084	54,572	22,910
6th	75,030	61,623	50,305	24,885

,

Table 5. Summary statistics of the prediction error with the CNN model.

Market	Mean	Median	SD	MAD
1st	157,291	121,796	120,665	53,156
2nd	119,827	97,033	88,847	39,610
3rd	88,953	71,184	66,079	26,554
4th	85,012	68,360	61,345	27,868
5th	60,332	46,474	48,765	17,705
6th	77,024	64,683	50,501	25,960

,

Table 6. Summary statistics of the prediction error with the LSTM model.

Market	Mean	Median	SD	MAD
1st	807,905	631,225	525,548	241,955
2nd	454,070	438,301	207,141	115,938
3rd	493,333	407,004	299,729	120,577
4th	305,885	284,560	142,275	77,029
5th	621,313	507,990	408,766	177,791
6th	369,563	325,991	185,801	86,222

,

Table 7. Summary statistics of the prediction error with the BiLSTM model.

Market	Mean	Median	SD	MAD
1st	9,949,320	7,113,820	8,181,740	3,309,050
2nd	1,335,310	509,853	1,676,470	243,729
3rd	396,317	368,620	163,949	80,360
4th	363,259	346,863	189,974	95,969
5th	318,561	298,485	156,222	78,313
6th	4,114,860	3,600,920	2,997,620	1,386,870

,

Table 8. Summary statistics of the prediction error with the CNN+LSTM model.

Market	Mean	Median	SD	MAD
1st	237,219	192,560	165,774	73,009
2nd	270,741	238,881	190,208	105,736
3rd	224,195	187,946	153,898	74,907
4th	199,779	164,432	142,618	67,457
5th	178,912	144,792	132,740	61,505
6th	180,031	158,938	129,836	74,575

and

Table 9. Summary statistics of the prediction error with the Transformer Model.

Market	Mean	Median	SD	MAD
1st	202,487	157,074	146,240	47,723
2nd	236,436	192,165	168,999	74,396
3rd	203,799	164,773	146,345	55,942
4th	182,406	146,881	132,737	55,306
5th	161,450	129,051	116,126	43,136
6th	143,626	120,655	98,667	41,839

show the summary statistics of the prediction error using the $d_K$ measure. We see that the MLP and CNN models (Table 4 and Table 5, respectively) generally obtain good results when compared to the naive model (Table 3). This is not only observed in the centrality measures (mean and median) but also in the dispersion measures (SD and MAD). This is relevant because it implies that the errors with MLP and CNN are smaller and also more concentrated. LSTM-type models, i.e., LSTM, BiLSTM, and CNN+LSTM, do not perform well since their mean and median values are generally much higher than those obtained using the naive procedure. The Transformer model (Table 9) has mixed behavior.

We now group the mean prediction errors using the $d_K$ distance to facilitate the comparison of the models. In

Table 10. Comparison of the mean prediction errors using the $d_K$ measure.

Market	Naive	MLP	CNN	LSTM	BiLSTM	CNN+LSTM	Transformer
1st	161,035	182,371	157,291	807,905	9.94932 × 106	237,219	202,487
2nd	201,648	126,912	119,827	454,070	1.33531 × 106	270,741	236,436
3rd	137,397	106,281	88,953.2	493,333	396,317	224,195	203,799
4th	190,072	99,136.6	85,012.3	305,885	363,259	199,779	182,406
5th	198,397	71,913.8	60,331.9	621,313	318,561	178,912	161,450
6th	117,832	75,030.5	77,023.6	369,563	4.11486 × 106	180,031	143,626

The cells highlighted in red, green, and blue denote the three best-performing models.

, the top three models for each market are colored red, green, and blue. We see that the CNN model obtains the best results in the first five sessions and is the second best in the sixth. In the sixth session, the best model is MLP, although the difference with CNN is small. We also see that the Naive model obtains good results, which indicates that, despite its simplicity, it is a good benchmark. It can be seen that the naive method is the second best in the first session market, only surpassed by the CNN procedure. The Transformer model improves the naive model in the fourth and fifth sessions.

Clearly, the best results are obtained for the second to sixth sessions using the MLP and CNN models. The percentage reduction in the error means with the MLP model compared to the naive procedure is between 22.6% in the third session and 63.8% in the fifth session. The reductions with the CNN model are between 35.2% and 69.6% in those sessions. In the fourth and fifth sessions, where the Transformer model is better than the naive procedure, the reductions are 4.0% and 18.6%, respectively. Therefore, they are lower than with MLP and CNN.

4.2. Model Comparison

The pairwise results of the Diebold–Mariano test between the top three models (which are mainly the naive model, the MLP model, and the CNN model, based on the previous results) are presented in

Table 11. The results of the Diebold–Mariano test for the comparison of the naive and MLP models.

Market	DM	p-Value
1st	−4.9711	≪0.005
2nd	31.9851	≪0.005
3rd	17.5485	≪0.005
4th	32.7746	≪0.005
5th	46.5074	≪0.005
6th	19.0330	≪0.005

,

Table 12. The results of the Diebold–Mariano test for the comparison of the naive and CNN models.

Market	DM	p-Value
1st	6.1157	≪0.005
2nd	34.4745	≪0.005
3rd	21.8768	≪0.005
4th	34.0464	≪0.005
5th	47.9174	≪0.005
6th	18.4540	≪0.005

and

Table 13. The results of the Diebold–Mariano test for the comparison of the MLP and CNN models.

Market	DM	p-Value
1st	34.2805	≪0.005
2nd	6.7702	≪0.005
3rd	12.7954	≪0.005
4th	16.5296	≪0.005
5th	22.0421	≪0.005
6th	−1.9850	0.0239

. It should be noticed that the obtained p-values are, with the exception of only one, much less than $\alpha =0.005$, which is the value recommended in [37] in order to provide solid evidence, to reduce the rate of false positives, and improve the reproducibility of the research.

Table 11 and Table 12 present the results of the Diebold–Mariano test (7) for the comparison of the naive model versus our proposals, MLP, and the CNN models. It is observed that the MLP model significantly improves the benchmark in five of the intraday market sessions but is outperformed in the first session. For its part, the CNN model significantly improves the naive model in all sessions. The results in these two tables answer affirmatively to the two research questions (Q1 and Q2) that we posed in the introduction since we achieved better prediction results than the naive procedure using information from the market and from the session preceding the session to be predicted.

Finally, Table 13 shows the Diebold-Mariano test for the comparison of the MLP and CNN models. We find that the CNN model significantly improves the MLP in the first five market sessions. In the sixth session, although it is true that the mean error of the MLP is lower, the differences are not significant using $\alpha =0.005$.

In Figure 5, boxplots of the prediction errors using the naive, MLP, and CNN models are shown. The improvements of MLP and CNN are evident in the second to sixth markets, where the boxes of the proposed models are below the median of the naive model. In the case of the first market, the naive model is better than MLP and slightly (but significantly) worse than CNN. We also see that the interquartile range for the MLP and CNN models is smaller than for the naive procedure. Therefore, their predictions have less dispersion. This is particularly clear in the second to sixth sessions.

5. Conclusions and Extensions

We have concluded that the MLP and CNN models are significantly better than the naive model when it comes to predicting energy curves, broadly reducing the error incurred in some cases. These models use information from the same session from the previous day and information from the previous session on the same day. By improving the naive procedure in this way, we have responded positively to the two research questions posed. Between these two models, at the $\alpha =0.005$ significance level, it is concluded that CNN performs better in the first five markets, while in the last one, the performances are similar. On the other hand, we observe that the LSTM and BiLSTM models give, by far, the worst results. The CNN+LSTM composite model, again, improves the results thanks to the inclusion of the convolutional layers. The Transformer model does not manage to improve on the naive model, except in a couple of markets.

Non-uniform point grids as an approximation to supply curves have proven effective in obtaining accurate predictions, thereby reducing the size of the objects handled as well as the computational time required. However, this aspect is one of the limitations of our study in the sense that the size of the grid implies a high dimension of the input tensor of the neural network models. Therefore, it would be interesting to study procedures that allow for reducing this dimension.

Another limitation that also relates to the dimensionality of the input features is that we have not taken into account meteorological variables that influence the availability of electricity produced through solar or wind technologies. Therefore, a possible extension of this work (maintaining the techniques used) could be to increase the amount of market information used for predictions, for example, meteorological variables such as temperature, solar radiation, precipitation, wind, etc., which affect both generators and consumers, could also be included.

Although this project has focused on predicting supply curves, both the preprocessing and predictive models can be used for demand curves as well. By doing so, an estimate of the energy price would be obtained by intersecting both curves. This will be the subject of future research.