###### Abstract

Trading on electricity markets requires accurate information about the
realization of electricity prices and the uncertainty attached to the
predictions. We present a probabilistic forecasting approach for day-ahead
electricity prices using the fully data-driven deep generative model called
normalizing flows. Our modeling approach generates full-day scenarios of day-
ahead electricity prices based on conditional features such as residual load
forecasts. Furthermore, we propose extended feature sets of prior realizations
and a periodic retraining scheme that allows the normalizing flow to adapt to
the changing conditions of modern electricity markets. In particular, we
investigate the impact of the energy crisis ensuing from the Russian invasion
of Ukraine. Our results highlight that the normalizing flow generates high-
quality scenarios that reproduce the true price distribution and yield highly
accurate forecasts. Additionally, our analysis highlights how our improvements
towards adaptations in changing regimes allow the normalizing flow to adapt to
changing market conditions and enables continued sampling of high-quality day-
ahead price scenarios.

###### Abstract

Trading on electricity markets requires accurate information about the
realization of electricity prices and the uncertainty attached to the
predictions. We present a probabilistic forecasting approach for day-ahead
electricity prices using the fully data-driven deep generative model called
normalizing flows. Our modeling approach generates full-day scenarios of day-
ahead electricity prices based on conditional features such as residual load
forecasts. Furthermore, we propose extended feature sets of prior realizations
and a periodic retraining scheme that allows the normalizing flow to adapt to
the changing conditions of modern electricity markets. In particular, we
investigate the impact of the energy crisis ensuing from the Russian invasion
of Ukraine. Our results highlight that the normalizing flow generates high-
quality scenarios that reproduce the true price distribution and yield highly
accurate forecasts. Additionally, our analysis highlights how our improvements
towards adaptations in changing regimes allow the normalizing flow to adapt to
changing market conditions and enables continued sampling of high-quality day-
ahead price scenarios.

Multivariate Scenario Generation of Day-Ahead Electricity Prices using
Normalizing Flows

Hannes Hilgera,b, Dirk Witthauta,b, Manuel Dahmena, Leonardo Rydin Gorjãoc,
Julius Trebbiena,b, Eike Cramerd,a,∗

  * a

Forschungszentrum Jülich GmbH, Institute of Energy and Climate Research,
Energy Systems Engineering (IEK-10), Jülich 52428, Germany

  * b

Institute for Theoretical Physics, University of Cologne, 50937 Köln, Germany

  * c

Faculty of Science and Technology, Norwegian University of Life Sciences, 1432
Ås, Norway

  * d

RWTH Aachen University, Process Systems Engineering (AVT.SVT), Aachen 52074,
Germany

Keywords: Electricity Prices; Time Series Forecasting; Normalizing Flows;
Scenario Generation

.

Multivariate Scenario Generation of Day-Ahead Electricity Prices using
Normalizing Flows

Hannes Hilgera,b, Dirk Witthauta,b, Manuel Dahmena, Leonardo Rydin Gorjãoc,
Julius Trebbiena,b, Eike Cramerd,a,∗

  * a

Forschungszentrum Jülich GmbH, Institute of Energy and Climate Research,
Energy Systems Engineering (IEK-10), Jülich 52428, Germany

  * b

Institute for Theoretical Physics, University of Cologne, 50937 Köln, Germany

  * c

Faculty of Science and Technology, Norwegian University of Life Sciences, 1432
Ås, Norway

  * d

RWTH Aachen University, Process Systems Engineering (AVT.SVT), Aachen 52074,
Germany

Keywords: Electricity Prices; Time Series Forecasting; Normalizing Flows;
Scenario Generation

##  1 Introduction

The transition to a sustainable energy system challenges the stable operation
of electric power systems and markets [1]. Renewable power generation is
intrinsically uncertain and fluctuates on various time scales from minutes to
seasons [2, 3]. Furthermore, electricity markets evolve in time due to changes
in the generation mix or the regulatory framework. The European electricity
markets underwent a fundamental change in late 2021 caused by the energy
crisis induced by the Russian invasion of Ukraine leading to exploding prices
and substantial changes in the statistical behavior of the electricity prices
[4, 5]. In general, the distribution of prices is heavy-tailed with
occasionally negative values, and price changes are strongly correlated over
several hours [6]. We argue that any model of electricity prices should
reflect these properties and that a model should be probabilistic to capture
the intrinsic uncertainty of the market. Furthermore, electricity price
forecasting models need to be able to adapt to changes in the behavior of the
markets. For instance, changes such as the increasing prices resulting from
the Russian war against Ukraine would render a non-adaptable model useless.

In this article, we present a data-driven, adaptable, and probabilistic
forecasting model to generate scenarios of day-ahead electricity prices. Our
model learns the conditional distribution of day-ahead electricity prices
based on forecasts of external factors such as wind and solar power generation
and load. We model all 24 hours of day-ahead prices for a given day as a
multivariate joint probability distribution. This multivariate probabilistic
forecasting approach reflects the fundamental structure of the day-ahead
electricity markets, where all 24 hourly prices are set simultaneously [7]. To
learn the conditional probability distribution, we use conditional normalizing
flows [8, 9], which we previously used for wind power scenario generation [10]
and prediction of intraday electricity prices [11]. The conditional
normalizing flow is a deep generative model [12] based on invertible neural
networks [13]. The ensemble approach provides several advantages over simpler
methods like point forecasting or forecasting of mean and standard deviation
[14]. First, scenario forecasts can produce potentially complicated, non-
Gaussian forecast distributions. Second, each scenario is intrinsically
consistent, i.e., correlations between the time steps are considered and
reproduced. Finally, the generated scenarios enable the formulation and
solution of stochastic optimization of problems to plan ahead under
uncertainty [10, 15].

We design our model architecture to be robust to changes in the overall market
behavior such as the price increase resulting from the energy crisis in 2021.
The model inherits price, demand, and renewable power generation data from the
previous day as conditional inputs. Thus, the model can rapidly detect changes
and adapt accordingly. Furthermore, we propose a periodic model update through
regular retraining steps. The retraining allows the model to compensate for
fundamental changes in market structure and behavior such as regulatory
changes or the increasing share of renewables.

The model is trained and tested using data from the German-Luxembourg day-
ahead electricity market and power system. We evaluate the model performance
and provide a detailed statistical analysis, comparing predictions and the
actual price time series. The results show that the model reproduces the
intricate statistical properties of the price time series, including the
heavy-tailed distribution as well as conditional distributions, temporal
correlations, and the impact of the European energy crisis.

The article is organized as follows: We first provide some background on the
European electricity markets and review the state of the art in electricity
price forecasting in Section 2. Then, we describe the concept and
implementation of the normalizing flow in Section 3. Our results on the model
performance and the statistical properties of prices and scenarios are given
in Section 4. Finally, we summarize and discuss our results in Section 5.

##  2 Background

This Section reviews the structural setup of the European electricity markets
including the day-ahead bidding markets. In the second part of the Section, we
review the state-of-the-art in electricity price forecasting.

###  2.1 European electricity markets

Stable operation of an electric power system requires that power generation
and demand are balanced at all times [16]. In the European system, power
generation is mainly coordinated through trading on electricity markets on
different time scales, e.g., in hourly or quarter-hourly intervals. Each
market participant has to align the physical net amount of electrical energy
that is produced or consumed in a given time window to the “virtual” amount of
electrical energy that is bought or sold on the electricity markets in that
particular time window [17]. For instance, a wind farm operator is required to
market the exact amount of electricity produced in any given quarter-hourly
time window. This process ensures a physical balance between power generation
and demand on the system level. Residual imbalances lead to deviations of the
grid frequency from its set value of 50 Hz and are corrected in real time via
the load-frequency control systems [16, 17]. Generally, the daily and weekly
patterns of buy and sell decisions lead to complex fluctuations of electricity
prices [6].

Market participants may buy and sell electricity either via direct power
purchase agreements, which may be agreed on months or years in advance, or via
trading on an electricity exchange [17]. On the exchanges, electricity is
traded on the futures markets and the spot markets. Power futures are long-
term contracts that regard delivery dates months or years in advance. On the
spot markets, electricity is traded with delivery dates on the next day (day-
ahead) [18] or the same day (intraday) [19].

Trading is organized in bidding zones and we will focus on the Germany-
Luxembourg bidding zone (Germany-Austria-Luxembourg until October 1, 2018).
For this article, we will restrict our analysis to the European Power Exchange
EPEX Spot [7], which has the highest trading volume for the Germany-Luxembourg
bidding zone. Furthermore, we focus on the day-ahead market, the most
important spot market in terms of trading volume [18]. At EPEX Spot,
electricity is traded in hourly windows for the 24 hours of the following day.
Market participants place buy and sell orders until 12:00 noon. Then, the
hourly prices are determined according to the market clearing principle: The
highest price that finds a buyer in each hour is determined as the market
clearing price for that hour. Every unit of electricity is traded at the
market clearing price in each respective hour. This is commonly referred to as
“pay-as-cleared”. Predicting this market clearing price is the central
objective of this article.

The European electricity markets are strongly affected by the global energy
crisis of 2021 and 2022. Energy prices soared in many regions of the world
following the Russian invasion of Ukraine [4]. Europe was particularly
strongly affected due to the dependence on Russian fossil fuels of various
countries.

Figure 1: Time series of day-ahead mean prices of each day from April 20, 2016
to December 31, 2022. We consider October 1, 2021 as the beginning of the
2021/22 energy crisis (shaded period). Data from EPEX Spot, taken from the
ENTSO-E transparency platform [20].

Figure 1 shows the daily average day-ahead prices in the Germany-Luxembourg
bidding zone from April 20, 2016, to December 31, 2022. The average price
level soared from around 30 EUR/MWh before the crisis to around 200 EUR/MWh
during the crisis, with peaks up to 800 EUR/MWh. Notably, the energy crisis
began well before the actual invasion (late February 2022) due to rising
political tensions in the preceding months. As the beginning is not clearly
defined, we use October 1, 2021, as a reference date during our analysis.

###  2.2 Electricity price forecasting

The field of electricity price forecasting is well established and receives
contributions from economics and technical fields like engineering, computer
science, and physics [21]. Most works are concerned with day-ahead electricity
prices [22], but there is a growing number of works considering intraday
electricity prices as well [23, 24], e.g., our previous work on normalizing
flows [11].

Traditionally, electricity price forecasting relied on statistical time series
models such as autoregressive (ARIMA, LASSO) models [21]. However, with the
increase in computing power and research on neural network regression, deep
learning became one of the drivers for continuous development in electricity
price forecasting [25]. Here, artificial neural networks and time series
neural networks like Long-Short Term Memory models (LSTMs) are the workhorse
methods [26, 27]. Recent advances in machine learning have benefited both
model development as well as feature selection for forecasting. For instance,
our previous work uses SHapley Additive exPlanations(SHAP) values to dissect
the functional relationship between electricity prices and relevant features
beyond the merit order principle [28]. In a similar work, Tschora et al. [29]
use SHAP values to identify correlations between bidding zones to improve
their forecasting performance.

Despite the increased understanding of modern electricity markets, the
realization of day-ahead electricity prices remains a stochastic process.
Thus, measures of uncertainty such as probabilistic forecasts can greatly
improve the reliability of the predictions [30]. Other approaches to quantify
the uncertainty include ensemble forecasts [31] or generation of prediction
intervals for neural network forecasts [32].

Most of the published approaches to forecasting day-ahead electricity prices
rely on a step-by-step forecasting approach, e.g., in autoregressive models
[22]. Notably, such an approach contrasts the actual procedure of settling the
day-ahead bidding markets, where all 24 hourly price intervals are set
simultaneously (c.f. Section 2.1). Instead, multivariate forecasting matches
the fundamental structure of the day-ahead market. Ziel and Weron [33] compare
univariate and multivariate forecasting and report improved performance for
the multivariate case. Other works combine multivariate forecasting with
Schaake shuffles to obtain probabilistic methods [34].

There are a few works considering adaptations towards the changing prices
resulting from the 2021/2022 European energy crisis. Examples include adaptive
preprocessing [35] and our previous work on probabilistic forecasting using
LSTMs [27].

##  3 Methods and Data

###  3.1 Fundamentals of normalizing flows

Normalizing flows are a class of deep generative models using invertible
transformations. The concept of normalizing flows was first introduced by
Tabak and Vanden-Eijnden [36] and Tabak and Turner [37] about ten years ago. A
generative model describes the probability distribution of a given data set
and can generate new samples from that distribution. Notably, other generative
models like variational autoencoders [38] and generative adversarial networks
[39] give an implicit representation of the probability distribution, i.e.,
they only allow for sampling. Normalizing flows, however, provide an explicit
representation of the probability distribution, i.e., the probability density
function (PDF), which enables mathematically consistent and efficient training
via likelihood maximization. We refer to Papamakarios et al. [13] for a
comprehensive review of normalizing flows.

The target data, in our case the day-ahead electricity prices, is represented
by a random vector $X\in\mathbb{R}^{D}$. The model learns a diffeomorphism
[13], i.e., a differentiable invertible transformation

$\displaystyle f:\;$ $\displaystyle\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$
$\displaystyle x\mapsto f(x)$

that maps $X$ to another random variable $Z$ following a well-known base
distribution. The most common choice for the base distribution is a
multivariate standard normal Gaussian distribution, i.e., $Z\sim
N(0,\mathbf{I})$. Using the diffeomorphism, normalizing flows provide an
explicit representation of the PDF of the target variable $X$ via a change of
variables [13], i.e.,

$\displaystyle p_{X}(x)=p_{Z}(f(x))\left|\det J_{f}(x)\right|^{-1},$ (1)

where $J_{f}(x)$ denotes the Jacobian of the function $f$ at the point $x$.
This direct representation allows for sampling according to $p_{X}(x)$ by
first sampling $z$ from the Gaussian $p_{Z}(z)$ and then transforming it
through the inverse transformation, i.e., computing $x=f^{-1}(z)$.

Using the explicit PDF in Equation (1), a normalizing flow is trained via
likelihood maximization [13]. Let $x_{1},x_{2},\ldots,x_{N}$ denote the data
points from the respective training set. Then, the function $f$ is chosen such
that it minimizes the negative log-likelihood

$\displaystyle
NLL=-\sum_{i=1}^{N}\log\left[p_{Z}\left(f(x_{i})\right)\left|\det
J_{f}(x_{i})\right|^{-1}\right].$ (2)

In practice, $f$ is chosen as an invertible neural network with a finite set
of parameters $\theta$.

The baseline normalizing flow can be extended to conditional statistics [9,
10], where the probability distribution depends on another variable
$y\in\mathbb{R}^{L}$. This conditional input is taken into account by
generalizing the flow to

$\displaystyle f:\;$
$\displaystyle\mathbb{R}^{D}\times\mathbb{R}^{L}\rightarrow\mathbb{R}^{D}$
$\displaystyle x,y\mapsto f(x,y).$

For every fixed value of $y$, the restricted function $x\mapsto f(x,y)$ must
be differentiable and invertible [8] w.r.t. $x$. Then, the conditional PDF is
given by

$\displaystyle p_{X|Y}(x|y)=p_{Z}\left(f(x,y)\right)\left|\det
J_{f}(x,y)\right|^{-1},$ (3)

where $J_{f}(x,y)$ denotes the Jacobian with respect to the variable $x$.

###  3.2 Model architecture and training

We implement the conditional normalizing flow using the real non-volume
preserving transformation (RealNVP) [40] with an extension to include
conditional features [8, 10]. RealNVP uses affine coupling layers that
construct highly flexible transformations that guarantee the invertibility of
the overall transformation. The coupling layers are built on so-called
conditioner models that introduce nonlinearity into the transformation. For
more details on normalizing flows and their implementation, we refer to the
review article by Papamakarios et al. [13], the original work on RealNVP by
Dinh et al. [40], and our previous work [14, 10].

As conditional inputs, we use the concatenation of seven 24-dimensional
forecast profiles, which amounts to a 168-dimensional conditional input vector
$y$ that is passed to the conditional RealNVP layers. The conditional inputs
include the day-ahead forecasts of wind and solar generation and load for
every hour of the following day as these features show the highest influence
on the realization of the day-ahead prices [28]. Furthermore, the conditional
input also includes the wind, solar, and load forecasts and the day-ahead
price realization of the previous day. The latter information allows the model
to scale the predicted day-ahead prices. Recall that the same wind, solar, and
load vectors result in very different day-ahead prices before and during the
energy crisis [27]. Therefore, any model that is trained prior and deployed
during the energy crisis is likely to perform poorly. Including information
from the previous day solves this problem for two reasons. First, it provides
a typical price level for the respective period. Second, the model can learn
that a certain set of wind, solar, and load profiles resulted in a certain
day-ahead price profile on the previous day. Including this additional
information enables the model to extrapolate what the wind, solar, and load
forecasts for the next day might result in. The robustness of the model
performance is assessed in detail in Section 4.

We scale all power data, i.e., the wind, solar, and demand data, by a factor
of 1.1 times their historical maximum to obtain features between 0 and 1. All
price data is scaled by a constant factor of 100. Note that normalizing flows
are not restricted to any specific interval, but scaling the data typically
improves its performance [13].

In the final stage, the model contains a decoding step that reduces the
dimensionality of the day-ahead electricity price data. By this step, we
mitigate a problem that repeatedly occurs in energy time series forecasting:
The strong correlation of time steps means that the target data $X$ lies on a
lower-dimensional manifold in the target space $\mathbb{R}^{D}$ [14]. In such
a case, normalizing flows typically learn smeared-out distributions [41] and
generate noisy scenarios [14]. We mitigate this problem by dimensionality
reduction to a lower dimensional space using principal component analysis
(PCA) [12, 14]. That is, we encode an original data point $x$ according to
$x^{\prime}:=U^{\top}(x-\bar{x})$, where $U$ is the matrix of principal
components, $\bar{x}$ is the sample mean, and $\top$ denotes the transpose.
The normalizing flow is trained on the encoded data
$x^{\prime}\in\mathbb{R}^{D^{\prime}}$, and scenarios are decoded using the
inverse of $U^{\top}$, i.e., $x:=\bar{x}+U\,x^{\prime}$. In practice, we use
an encoding into $D^{\prime}=14$ dimensions, which explains $99.5\%$ of the
variance of the original data.

To test the performance of the normalizing flow, we do not use a fixed train-
test-split but implement a retraining scheme: Every 90 days, the normalizing
flow is newly trained on all available data until that point. For instance,
the normalizing flow might be newly trained at the end of 2018 with all data
available until then (Jan 2016 - Dec 2018). This retraining also includes
adjustments of the scaling factors for preprocessing, if necessary. The newly
trained normalizing flow is then used for scenario generation for the
following 90 days. For instance, at the beginning of April 2019, the
normalizing flow is then retrained again with all available data (Jan 2016 -
Mar 2019). It is this retraining scheme that allows the model to take into
account changing market conditions, as the normalizing flow regularly gains
new training samples that might exhibit novel market behavior.

###  3.3 Implementation and hyperparameter optimization

The normalizing flow is implemented in Python 3.9.13 using TensorFlow 2.10.0
[42]. The code for the normalizing flow is based on our prior studies [10] and
open source libraries from [43]. The PCA calculations are done using scikit-
learn 1.1.2 for Python [44].

Table 1: Hyperparameter optimization of the normalizing flow. We test different combinations of hyperparameters in two steps and evaluate the performance in terms of the mean absolute error (MAE). hyperparameter | values (1st step) | values (2nd step)
---|---|---
coupling layers | 2, 3, 4, 5 | 3, 4, 5
network depth | 2, 3, 4, 5 | 2, 3, 4
network width | 14, 21, 28 | 14, 21
epochs | 500, 750, 1000, 1500 | 1000

We use fully connected neural networks to implement the conditioner models.
Thus, the model contains four hyperparameters: the number of coupling layers,
the depth of each network describing the conditioner models, the number of
nodes in each hidden layer of the conditioner models, and the number of
training epochs. These hyperparameters are optimized in two steps using on the
JURECA DC supercomputer at Forschungszentrum Jülich [45].

For hyperparameter optimization, we train the normalizing flow on all
available data except for the first 90 days, which are set aside as the test
set. First, we train one model instance for the 192 different hyperparameter
combinations as listed in the center of Table 1. In each case, we evaluate the
mean absolute error (MAE) of the scenario mean and discard hyperparameter
values that lead to high MAE values. We find that normalizing flows with just
two coupling layers tend to underfit the data and thus discard this
configuration. In a second step, we train each parameter combination eight
times and evaluate the mean and the standard deviation of the MAE to avoid an
influence from stochastic effects in the training. Therefore, we reduced the
number of parameter combinations according to the results of the first step,
keeping only 18 combinations. We list the six best-performing hyperparameter
combinations in Table 2.

Table 2: Results of the second step of hyperparameter optimization. We only
show the six hyperparameter combinations with the lowest averaged MAE. For
each hyperparameter combination, we train eight models and report the mean and
the standard deviation over the eight runs.

| coupling
---
layers
| hidden
---
layers
| hidden
---
nodes
| mean absolute error
---
[EUR/MWh]
5 | 2 | 21 | $\textbf{11.11}\pm\textbf{0.56}$
4 | 2 | 14 | $11.17\pm 0.24$
3 | 2 | 21 | $11.21\pm 0.34$
3 | 3 | 21 | $11.29\pm 0.20$
4 | 3 | 21 | $11.31\pm 0.28$
4 | 2 | 21 | $11.32\pm 0.21$

We find that the differences in performance between the different models are
small and therefore the choice of hyperparameters appears to only play a minor
role in the examined ranges. In the following, we choose the best-performing
hyperparameter combination w.r.t. the MAE (coupling layers: 5, number of
hidden layers: 2, number of nodes: 21, epochs: 1000).

###  3.4 Benchmark models

To assess the performance of the normalizing flow, we consider two benchmark
models for scenario generation. First, an uninformed historical model
generates samples by randomly drawing from the pool of past price
realizations. Hence, the model ignores all forecasts but captures typical
daily profiles. We include this model to represent a valid reference point and
lower bound for the model performance examination. Second, an informed
historical model generates samples using a k-nearest-neighbors approach. It
generates scenarios by drawing the historical price realizations of days with
the closest conditional inputs. In other words, the generated scenarios
consist of price profile realizations of the historical days with the most
similar conditions. The conditional vectors are a 96-dimensional concatenation
of wind, solar, and load forecasts and the price realization of the previous
day. The k-nearest-neighbors model is implemented using the NearestNeighbors
function from scikit-learn 1.1.2 in Python [44].

###  3.5 Data Sources

We use data from the ENTSO-E transparency platform [20] from January 2016 to
December 2022, which were retrieved via the restful API provided by ENTSO-E
[7] using the entsoe-py open-source implementation for Python. The day-ahead
price is the price of the EPEX Spot day-ahead auction, for the Germany-
Luxembourg bidding zone (Germany-Austria-Luxembourg prior to October 1st
2018). The day-ahead load forecast is the expected hourly load in the Germany-
Luxembourg bidding zone. The day-ahead solar and wind forecasts are the
expected hourly production of each generation type in the Germany-Luxembourg
bidding zone. We use the ENTSO-E forecasts because they provide a coherent
publicly available reference data source, although market participants
typically use a variety of different forecasting products (cf. the discussion
in [46]).

##  4 Results

This Section applies the RealNVP-based normalizing flow to generate scenarios
of the day-ahead electricity price. We start by showing examples of generated
scenarios in reference to the benchmark models. Next, we perform a statistical
analysis of the generated scenarios and, finally, evaluate the scenarios using
ensemble forecasting scores such as the energy score and the variogram score.

###  4.1 Initial examples

This Section provides a qualitative overview of the capabilities of the
conditional normalizing flow. Figure 2 shows three selected examples of
ensemble forecasts and the associated conditional inputs. The first row of
Figure 2 shows a typical day in May 2017. The load profile of that day is on a
low level and, thus, the production from solar and wind is also low. This
translates into a typical price profile with two peaks, one in the morning and
another in the afternoon. Prices are lower at noon due to stronger solar power
generation and during the night due to a lower load. Overall, the shape and
price level of the realization are well predicted by the scenarios from the
normalizing flow.

Figure 2: Example forecasts for May 7, 2017 (top), November 28, 2017 (center)
and August 22, 2020 (bottom). Left column shows solar generation forecast
(yellow), wind generation forecast (blue) and load forecast (red) for each
day. Right column show 50 generated scenarios (blue) according to the
conditions forecasts and respective price realization (black) for comparison.

The second row in Figure 2 shows a day, where the expected wind power
production is high in the morning hours but decreases throughout the day. This
is well reflected in the day-ahead price profile scenarios, where the price
peak in the afternoon is higher due to a higher residual load compared to that
in the morning hours. Again, the generated scenarios tightly mirror the actual
realization.

The third row in Figure 2 shows a day where the expected load is low (a
typical Saturday), and the solar and wind production are both expected to be
quite high, especially during the noon hours. This results in a day-ahead
price realization with a deep price dip to almost 0 EUR/MWh around noon. The
model predicts this price dip, with some scenarios even going into the
negative price range. At this hour, the distribution of the predicted prices
becomes strongly non-Gaussian with a clear negative skew.

###  4.2 Statistical verification of normalizing flow-generated scenarios

Electricity price time series have intricate statistical properties [6], e.g.,
heavy-tailed PDFs. In this Section, we analyze whether the normalizing flow is
able to reproduce the statistics of the actual time series. To this end, we
compare the histograms of hourly prices in the realizations and scenarios as
well as the leading statistical moments. In the scenario histograms, we scale
the number of occurrences by the number of samples in order to match the
realization histograms.

Figure 3: Histogram of prices of all generated scenarios compared to the
histogram of the actual day-ahead price time series (“realizations”). Dotted
line is a Gaussian fit onto the realizations histogram. The value $D$ gives
the Kullback–Leibler divergence between scenario and realization histogram.
Time series ranges from April 20, 2016 to December 31, 2022.

Figure 3 shows the histogram for the entire period of analysis from April 20,
2016 to December 31, 2022. As motivated in our previous work [47], Figure 3
shows the histogram in logarithmic scaling to allow for an analysis of the
tails of the distribution. Overall, the scenario histogram matches the
realizations histogram very well, which is also reflected in the similarity of
the statistical moments listed in Table 3. The scenarios slightly
underestimate the likelihood of high prices. We attributed this discrepancy to
the fact that the model needs a certain amount of training data from the new
regime to adjust to changing market conditions. At the onset of the energy
crisis, the normalizing flow systematically underestimates the electricity
prices as the price increase is not yet included in the training data.
However, this period is rather short due to our retraining scheme such that we
observe a very good overall agreement.

Table 3: Mean $\mu$, standard deviation $\sigma$, skewness $s$, and kurtosis $k$ of the day-ahead price time series (“Realizations”) and of the scenarios generated by the normalizing flow. We provide the normalized central moments for the entire time period under observation as well as separately for the time before and during the energy crisis. |  | $\mu$ | $\sigma$ | $s$ | $k$
---|---|---|---|---|---
20.04.16 - | realizations | 72.87 | 90.75 | 2.94 | 10.11
31.12.22 | scenarios | 70.19 | 84.98 | 2.75 | 9.16
20.04.16 - | realizations | 40.38 | 23.63 | 1.29 | 7.17
30.09.21 | scenarios | 39.17 | 22.83 | 1.54 | 8.57
01.10.21 - | realizations | 218.43 | 129.25 | 0.91 | 1.05
31.12.22 | scenarios | 209.74 | 115.97 | 0.74 | 1.72

We emphasize that the histograms have an unusual shape, which differs
considerably from the histograms for the period 2015 to 2019 analyzed in [6].
This is a direct result of the overlay of distributions from different market
regimes, i.e., before and during the energy crisis (Figure 1). For a more
detailed analysis, we show separate histograms for the two market periods in
Figure 4. The distribution of prices during the energy crisis vastly differs
from the distribution before the crisis. Notably, the scenarios show a good
overall match to the realizations, demonstrating the normalizing flow’s
capability to learn and sample from complex non-Gaussian distributions.

Figure 4: Histograms of prices of generated scenarios compared to histograms
of the actual day-ahead price time series (“Realizations”). The normalizing
flow is trained on all available data at the given time. The left side shows
histograms for time series before October 1, 2021. The right side shows
histograms for time series after October 1, 2021. Note the different scales on
the x-axis. Dotted lines present Gaussian fits onto the realizations
histograms. The value $D$ gives the Kullback–Leibler divergences between
scenario and realization histograms.

The histograms show that negative electricity prices almost never occur after
the onset of the energy crisis. However, the normalizing flow generates
scenarios that include negative prices. The absence of negative electricity
prices has both economic and regulatory reasons [48, 49]. In the German
market, wind turbines and large solar PV installations receive subsidies
(“Marktprämie”) that are given by the difference between a fixed reference
value (“Anzulegender Wert”) and the average market price level. In a high-
price market regime, the average market price level exceeds the reference
value and the subsidies drop to zero. In such a case, wind and solar plants
curtail generation to avoid negative prices. Hence, the price frequently
decreases to zero or small positive values but rarely to negative values.
Figure 4 shows a small peak around zero but very few values below zero. As the
proposed normalizing flow scheme is fully data-driven, this regulatory
mechanism cannot be enforced explicitly. The scenarios fail to represent the
respective features of the distribution and the scenario histogram is
smoothed, missing the peak at zero and the sharp decrease below. We argue that
this mismatch results from the change in regime being ahead of the adoption of
training data. Moreover, the training data for the normalizing flow after the
onset of the energy crisis still includes data from previous years containing
negative prices.

Beyond the full distributions, we emphasize that the normalizing flow also
reproduces the _marginal_ distributions. Figure 5 shows the histograms for two
hourly windows starting at 6:00 and 12:00. The probability for high prices is
higher at 6:00, while the probability for low or negative prices is higher at
12:00. Again, this is well explainable through a typical solar profile and the
merit order effect. The generated scenarios reflect this behavior and produce
different distributions for different hours of the day.

Figure 5: Marginal price histogram of generated scenarios vs. realizations at
6:00 (left) and 12:00 (right).

###  4.3 Forecasting performance

We provide a quantitative assessment of the performance of the normalizing
flow in reference to the two benchmark scenario generation methods. First, we
show the absolute error and its mean (MAE), i.e., the MAE of the hourly mean
values of the generated scenarios. We emphasize that the MAE is designed to
evaluate point forecasts. Thus, our MAE analysis is limited to the mean of the
generated scenarios. Results are provided in Figure 6 in comparison to the two
benchmark models introduced in Section 3.4. We find that the normalizing flow
strongly outperforms the two benchmark models in terms of the MAE. In
particular, under shifting market conditions such as in 2022, the normalizing
flow approach with retraining holds up well. For the period from 2019/01/30 to
2020/02/08, we find an MAE of $3.88$ EUR/MWh for the mean value of the
normalizing flow scenarios. This value is comparable to our recent results
using LSTMs [27], reporting state-of-the-art performance with an MAE of $3.73$
EUR/MWh for the year 2019. For the period between the year 2019-2022, our
previous work [27] finds an MAE of $11.92$ EUR/MWh. Again, the normalizing
flow yields competitive results with an MAE of $11.11\pm 0.56$ EUR/MWh over
the entire period of 2016–2022 (Table 2).

Considering the different time periods, the results from [27] are slightly
better. The time before the energy crisis which generally leads to a lower
MAE, due to the lower absolute price values, is more strongly represented in
the full data set used for the normalizing flow training. Nevertheless, we
find that the normalizing flow is generally competitive even in terms of the
MAE despite not being designed to produce point forecasts.

Figure 6:  Performance of the normalizing flow in comparison to the benchmark
models. The upper plots compare the mean absolute error for 2019 (left) and
2022 (right). The center plots compare the energy score for 2019 (left) and
2022 (right). The bottom plots compare the variogram score for 2019 (left) and
2022 (right). The black vertical bar indicates the sample median. The boxes
indicates the ranges between 75% and 25%, and the whiskers indicate the range
between 97.5% and 2.5%.

Opposed to the MAE, there are metrics to specifically evaluate the quality of
probabilistic forecasts. As in our previous work on intraday price forecasting
[11], we use the energy score (ES) and the variogram score (VS). The energy
score [50, 51] is defined as

$\displaystyle\text{ES}=\frac{1}{N}\sum_{s=1}^{N}\left\|\mathbf{\lambda}-\mathbf{\hat{\lambda}_{s}}\right\|_{2}-\frac{1}{2N^{2}}\sum_{s=1}^{N}\sum_{s^{\prime}=1}^{N}\left\|\mathbf{\hat{\lambda}_{s}}-\mathbf{\hat{\lambda}_{s}^{\prime}}\right\|_{2}.$
(4)

Here, $\lambda$ is the 24-dimensional realized price profile for a given day
and $\hat{\lambda}_{s}$ is the price profile per scenario $s$. The operator
$\left\|\cdot\right\|_{2}$ denotes the Euclidean norm and $N$ is the number of
scenarios used to compute the energy score. The first term on the right side
of Equation (4) measures the distance between the scenarios and the
realization. The second term measures the diversity of the samples. The VS
[52] quantifies whether the forecasts correctly describe the correlations
between the individual time steps. It is defined as

$\displaystyle\text{VS}=\frac{1}{N}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}\left(|\lambda_{t}-\lambda_{t^{\prime}}|^{\gamma}-\frac{1}{N}\sum_{s=1}^{N}|\hat{\lambda}_{t,s}-\hat{\lambda}_{t^{\prime},s}|^{\gamma}\right)^{2}.$
(5)

The parameter $\gamma$ is referred to as variogram order and is typically set
to $\gamma=0.5$ [52]. Both the energy score and the variogram score are
negatively oriented, i.e., a lower score indicates a better result. Again, $N$
is the number of scenarios and $T$ is the number of time steps within each
scenario.

Figure 6 shows box plots of the MAE, the energy score, and the variogram score
distributions in comparison to the two benchmark models, before and after the
energy crisis. The results show that the normalizing flow yields substantially
lower values for both scores, thus indicating a much better agreement with the
realizations. Furthermore, the normalizing flow consistently outperforms the
benchmark methods for both periods and even increases its advantage in the
period of 2022. Notably, the absolute values of the MAE, energy score, and
variogram score increase by about a factor of ten with the onset of the energy
crisis. This increase is expected as the absolute prices increase by a similar
factor. We made the same observation in our previous work on price forecasting
using LSTMs [46].

In summary, the analysis in this Section shows how conditional normalizing
flows can generate high-quality scenarios of day-ahead electricity prices. The
normalizing flow-generated scenarios consistently outperform the benchmark
scenarios in all evaluations. Notably, the normalizing flow achieves state-of-
the-art forecasting performance for deterministic forecasting and realistic
scenario generation. The generated scenarios allow for quantification of the
prediction uncertainty. The adaptive retraining of the normalizing flow
produces high-quality results throughout the transition of market regimes.

###  4.4 Correlations

Figure 7:  Reproduction of correlations in the price time series. We
investigate the joint probability distribution for two points in time (left:
$t_{1}$ = 3:00 and $t_{2}$=5:00, right: $t_{1}$ = 6:00 to $t_{2}$=8:00). We
show the histograms of the realizations (top panels) and scenarios (center
panels). The lower panel shows statistics of the increments
$\Delta=\mbox{price}(t_{2})-\mbox{price}(t_{1})$.

An important advantage of multivariate scenario forecasting is that each
scenario is intrinsically consistent, i.e., every generated scenario reflects
correlations present in the actual price time series. Mathematically speaking,
the normalizing flow learns the distribution of a random vector $X$ describing
the prices of an entire day instead of an individual hour.

We test the capability of the model by fixing two points in time, $t_{1}$ and
$t_{2}$, and comparing the _joint_ probability distribution of the respective
prices. Figure 7 shows the histograms of the occurrences for the two
respective times. In the early morning, prices increase from $t_{1}$ = 3:00 to
$t_{2}$=5:00 in a characteristic way (cf. Figure 2). Hence, the joint PDF is
concentrated above the bisector. Later, between $t_{1}$ = 6:00 to
$t_{2}$=8:00, the prices mostly decrease and the joint PDF is concentrated
slightly below the bisector. In both cases, Figure 7 shows that the
normalizing flow reproduces the joint PDF aptly. Hence, the normalizing flow
successfully learns the correlations between the different points in time.

For a more detailed analysis, we consider the price increments
$\Delta=\mbox{price}(t_{2})-\mbox{price}(t_{1})$ and compute their histogram.
Overall, we find a good agreement of the scenarios and realizations in terms
of the increment statistics. However, in both examples, the actual
realizations show a sharp peak at an increment of $\Delta\approx 0$, which is
not reproduced. We believe this peak results from complex regulatory aspects
of the market. For instance, we have already discussed in Section 4.2 that the
regulation of renewable subsidies leads to an increased likeliness of a price
of 0 EUR/MWh or slightly above. Hence, there is an increased likeliness that
the price stays at one value for several hours leading to an increment of
$\Delta\approx 0$ EUR/MWh. The normalizing flow does not learn this
characteristic such that the increment distribution is smoothed compared to
the actual data. Again, we expect this model behavior to change with the
inclusion of more training data from the later periods and the exclusion of
data from earlier periods where negative prices were more common.

###  4.5 Errors and uncertainties

Quantification of forecasting uncertainty is of high importance in many
applications. We consequently study whether the normalizing flow can provide a
measure of confidence for its forecasts. In particular, we examine the
following question: If the model’s scenario forecast is off for a particular
hour, did the model express uncertainty about the outcome? Figure 8 compares
the standard deviation of the hourly forecast distribution to the absolute
error of the mean value of the generated scenarios. The scatter plot reveals
that higher values for the absolute error of the expected forecast usually
correlate with a higher forecast standard deviation, as events with a high
absolute error of the expected forecast and low forecast standard deviation
occur only rarely. In other words, the normalizing flow rarely makes bad
predictions with high confidence. Notably, there is a lower bound of the
standard deviation in the forecasts, which increases over time with the change
of the market regime.

Similar to Figure 6, the change in behavior over time in Figure 8 shows
increasing absolute errors and standard deviations for later periods. Again,
this behavior is consistent with our expectation of adjusting towards the
high-price regime with higher variance after the onset of the energy crisis.

Figure 8: Standard deviation of hourly forecast distribution against absolute
error of the expected forecast. Each dot represents one hour. Color represents
the date according to the colorbar.

##  5 Conclusion and Outlook

We present a multivariate probabilistic forecasting approach for day-ahead
electricity prices based on normalizing flows. Our normalizing flow
implementation incorporates relevant feature information to learn the
conditional multivariate probability distribution of the vector of day-ahead
electricity prices. We train our model via direct log-likelihood maximization
to achieve mathematically consistent and efficient training. The trained model
allows for sampling day-specific scenarios of electricity price time series
that are intrinsically consistent and match the fundamental market structure
of the day-ahead bidding market of the EPEX Spot markets by generating full
24-hour scenarios.

Our analysis shows that the normalizing flow yields high-quality scenarios
with a good representation of the actual price realization and informative
uncertainty quantification that indicates the reliability of the forecasts in
a quantitative way. The conditional normalizing flow significantly outperforms
two benchmark scenario generation approaches based on historic sampling.

We propose a periodic retraining scheme to continuously adapt the normalizing
flow to the changes in market regimes such as the onset of the energy crisis,
in 2021. With brief delays, the normalizing flow adapts to the changing
markets and generates high-quality scenarios. This retraining scheme could
prove useful for analyzing and modeling other strongly non-stationary time
series as well.

## Acknowledgments

The authors gratefully acknowledge the computing time granted through JARA on
the supercomputer JURECA [45] at Forschungszentrum Jülich. E.C. gratefully
acknowledges the financial support of the Kopernikus project SynErgie 3 by the
Federal Ministry of Education and Research (BMBF) and the project supervision
by the project management organization Projektträger Jülich. D.W., H.H., J.T.,
and M.D. received funding from the Helmholtz Association of German Research
Centers.

## Bibliography

  * [1] Dirk Witthaut, Frank Hellmann, Jürgen Kurths, Stefan Kettemann, Hildegard Meyer-Ortmanns, and Marc Timme.  Collective nonlinear dynamics and self-organization in decentralized power grids.  Reviews of Modern Physics, 94(1):015005, 2022.
  * [2] Mehrnaz Anvari, Gerald Lohmann, Matthias Wächter, Patrick Milan, Elke Lorenz, Detlev Heinemann, M Reza Rahimi Tabar, and Joachim Peinke.  Short term fluctuations of wind and solar power systems.  New Journal of Physics, 18(6):063027, 2016.
  * [3] Iain Staffell and Stefan Pfenninger.  The increasing impact of weather on electricity supply and demand.  Energy, 145:65–78, 2018.
  * [4] Andreas Goldthau and Simone Tagliapietra.  Energy crisis: five questions that must be answered in 2023.  Nature, 612(7941):627–630, 2022.
  * [5] Philipp C Böttcher, Leonardo Rydin Gorjão, Christian Beck, Richard Jumar, Heiko Maass, Veit Hagenmeyer, Dirk Witthaut, and Benjamin Schäfer.  Initial analysis of the impact of the ukrainian power grid synchronization with continental europe.  Energy Advances, 2(1):91–97, 2023.
  * [6] Chengyuan Han, Hannes Hilger, Eva Mix, Philipp C Böttcher, Mark Reyers, Christian Beck, Dirk Witthaut, and Leonardo Rydin Gorjão.  Complexity and persistence of price time series of the european electricity spot market.  PRX Energy, 1(1):013002, 2022.
  * [7] European Network of Transmission System Operators for Electricity.  Transparency platform RESTful API – user guide.  https://transparency.entsoe.eu/content/static_content/Staticcontent/webapi/Guide.html, 2023.
  * [8] Christina Winkler, Daniel Worrall, Emiel Hoogeboom, and Max Welling.  Learning likelihoods with conditional normalizing flows, 2020.
  * [9] Kashif Rasul, Abdul-Saboor Sheikh, Ingmar Schuster, Urs M Bergmann, and Roland Vollgraf.  Multivariate probabilistic time series forecasting via conditioned normalizing flows.  In International Conference on Learning Representations, 2021.
  * [10] Eike Cramer, Leonard Paeleke, Alexander Mitsos, and Manuel Dahmen.  Normalizing flow-based day-ahead wind power scenario generation for profitable and reliable delivery commitments by wind farm operators.  Computers & Chemical Engineering, 166:107923, 2022.
  * [11] Eike Cramer, Dirk Witthaut, Alexander Mitsos, and Manuel Dahmen.  Multivariate probabilistic forecasting of intraday electricity prices using normalizing flows.  Applied Energy, 346:121370, 2023.
  * [12] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep learning.  MIT press, 2016.
  * [13] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan.  Normalizing flows for probabilistic modeling and inference.  The Journal of Machine Learning Research, 22(1):2617–2680, 2021.
  * [14] Eike Cramer, Alexander Mitsos, Raúl Tempone, and Manuel Dahmen.  Principal component density estimation for scenario generation using normalizing flows.  Data-Centric Engineering, 3:e7, 2022.
  * [15] Mario Beykirch, Tim Janke, and Florian Steinke.  Bidding and scheduling in energy markets: Which probabilistic forecast do we need?  In 2022 17th International Conference on Probabilistic Methods Applied to Power Systems (PMAPS), pages 1–6. IEEE, 2022.
  * [16] Prabha Kundur.  Power system stability.  Power system stability and control, 10:7–1, 2007.
  * [17] Bundesnetzagentur.  Definitionen der marktakteuere und deren daten.  https://www.bundesnetzagentur.de/SharedDocs/Downloads/DE/Sachgebiete/Energie/Unternehmen_Institutionen/DatenaustauschUndMonitoring/MaStR/DefinitionenMarktakteuere.pdf, 2023.
  * [18] Ronald Huisman, Christian Huurman, and Ronald Mahieu.  Hourly electricity prices in day-ahead markets.  Energy Economics, 29(2):240–248, 2007.
  * [19] Priyanka Shinde and Mikael Amelin.  A literature review of intraday electricity markets and prices.  In 2019 IEEE Milan PowerTech, pages 1–6, 2019.
  * [20] European Network of Transmission System Operators for Electricity.  ENTSO-E transparency platform.  https://transparency.entsoe.eu, 2023.
  * [21] Rafał Weron.  Electricity price forecasting: A review of the state-of-the-art with a look into the future.  International Journal of Forecasting, 30(4):1030–1081, 2014.
  * [22] Jesus Lago, Grzegorz Marcjasz, Bart De Schutter, and Rafał Weron.  Forecasting day-ahead electricity prices: A review of state-of-the-art algorithms, best practices and an open-access benchmark.  Applied Energy, 293:116983, 2021.
  * [23] Georg Wolff and Stefan Feuerriegel.  Short-term dynamics of day-ahead and intraday electricity prices.  International Journal of Energy Sector Management, 11:557–573, 2017.
  * [24] Marcel Kremer, Rüdiger Kiesel, and Florentina Paraschiv.  An econometric model for intraday electricity trading.  Philosophical Transactions of the Royal Society A, 379(2202):20190624, 2021.
  * [25] Arkadiusz Jedrzejewski, Jesus Lago, Grzegorz Marcjasz, and Rafał Weron.  Electricity price forecasting: The dawn of machine learning.  IEEE Power and Energy Magazine, 20(3):24–31, 2022.
  * [26] Gaurav Kapoor and Nuttanan Wichitaksorn.  Electricity price forecasting in new zealand: A comparative analysis of statistical and machine learning models with feature selection.  Applied Energy, 347:121446, 2023.
  * [27] Julius Trebbien, Sebastian Pütz, Benjamin Schäfer, Heidi S Nygård, Leonardo Rydin Gorjão, and Dirk Witthaut.  Probabilistic forecasting of day-ahead electricity prices and their volatility with LSTMs.  arXiv preprint arXiv:2310.03339, 2023.
  * [28] Julius Trebbien, Leonardo Rydin Gorjão, Aaron Praktiknjo, Benjamin Schäfer, and Dirk Witthaut.  Understanding electricity prices beyond the merit order principle using explainable AI.  Energy and AI, 13:100250, 2023.
  * [29] Léonard Tschora, Erwan Pierre, Marc Plantevit, and Céline Robardet.  Electricity price forecasting on the day-ahead market using machine learning.  Applied Energy, 313:118752, 2022.
  * [30] Jakub Nowotarski and Rafał Weron.  Recent advances in electricity price forecasting: A review of probabilistic forecasting.  Renewable and Sustainable Energy Reviews, 81:1548–1568, 2018.
  * [31] Michał Narajewski and Florian Ziel.  Ensemble forecasting for intraday electricity prices: Simulating trajectories.  Applied Energy, 279:115801, 2020.
  * [32] Abbas Khosravi, Saeid Nahavandi, and Doug Creighton.  Quantifying uncertainties of neural network-based electricity price forecasts.  Applied Energy, 112:120–129, 2013.
  * [33] Florian Ziel and Rafał Weron.  Day-ahead electricity price forecasting with high-dimensional structures: Univariate vs. multivariate modeling frameworks.  Energy Economics, 70:396–420, 2018.
  * [34] Oliver Grothe, Fabian Kächele, and Fabian Krüger.  From point forecasts to multivariate probabilistic forecasts: The schaake shuffle for day-ahead electricity price forecasting.  Energy Economics, 120:106602, 2023.
  * [35] Carlos Sebastián, Carlos E González-Guillén, and Jesús Juan.  An adaptive standardisation model for day-ahead electricity price forecasting.  arXiv preprint arXiv:2311.02610, 2023.
  * [36] Esteban G Tabak and Eric Vanden-Eijnden.  Density estimation by dual ascent of the log-likelihood.  Communications in Mathematical Sciences, 8(1):217–233, 2010.
  * [37] Esteban G Tabak and Cristina V Turner.  A family of nonparametric density estimation algorithms.  Communications on Pure and Applied Mathematics, 66(2):145–164, 2013.
  * [38] Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  In Yoshua Bengio and Yann LeCun, editors, 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014.
  * [39] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In Advances in Neural Information Processing Systems, pages 2672–2680. MIT Press, Cambridge, MA, USA, 2014.
  * [40] Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using realNVP.  arXiv preprint arXiv:1605.08803, 2016.
  * [41] Johann Brehmer and Kyle Cranmer.  Flows for simultaneous manifold learning and density estimation.  Advances in Neural Information Processing Systems, 33:442–453, 2020.
  * [42] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, et al.  Tensorflow: Large-scale machine learning on heterogeneous distributed systems.  arXiv preprint arXiv:1603.04467, 2016.
  * [43] Keras Team.  Keras documentation: Density estimation using real NVP.  https://keras.io/examples/generative/real_nvp, 2023.
  * [44] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  The Journal of machine Learning research, 12:2825–2830, 2011.
  * [45] Philipp Thörnig.  JURECA: Data Centric and Booster Modules implementing the Modular Supercomputing Architecture at Jülich Supercomputing Centre.  Journal of Large-Scale Research Facilities, 7:A182, 2021.
  * [46] Julius Trebbien.  Explainable artificial intelligence and deep learning for analysis and forecasting of complex time series: Applications to electricity prices.  https://kups.ub.uni-koeln.de/70766/, 2023.
  * [47] Eike Cramer, Leonardo Rydin Gorjão, Alexander Mitsos, Benjamin Schäfer, Dirk Witthaut, and Manuel Dahmen.  Validation methods for energy time series scenarios from deep generative models.  IEEE Access, 10:8194–8207, 2022.
  * [48] Amani Joas.  Energy price cap in germany & curtailment of renewables.  https://flex-power.energy/energyblog/energy-price-cap-tax-renewables-curtailment, 2022.
  * [49] Amani Joas.  The economics of curtailing renewables like wind & solar.  https://flex-power.energy/energyblog/the-economics-of-curtailing-renewables, 2023.
  * [50] Tilmann Gneiting and Adrian E Raftery.  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association, 102, 2007.
  * [51] P. Pinson and R. Girard.  Evaluating the quality of scenarios of short-term wind power generation.  Applied Energy, 2012.  Publisher: Elsevier.
  * [52] Michael Scheuerer and Thomas Hamill.  Variogram-based proper scoring rules for probabilistic forecasts of multivariate quantities.  Monthly Weather Review, 2015.


#  $f(R,T)$ Analogue Gravity in $(2+1)$ D graphene sheet

M. Lapola<EMAIL_ADDRESS>São Paulo State University (Unesp),
Institute of Geosciences and Exact Sciences, Physics Department  L. A.
Barreiro  São Paulo State University (Unesp),
Institute of Geosciences and Exact Sciences, Physics Department

(August 28, 2024)

###### Abstract

We examine the analogue gravity model within the context of f(R,T) gravity
applied to graphene. The derivation of the Lagrangian density in two
dimensions (2D) is undertaken, accounting for the altered gravitational
effects as characterized by the function f(R,T). The Lagrangian encompasses
the quasiparticle field $\psi(x)$, its adjoint $\overline{\psi}$, the
effective metric tensor $g^{\mu\nu}$, and the gauge field $A_{\nu}$. The
equations of motion are established through variational principles applied to
the Lagrangian, resulting in modified Dirac equations. We discuss the
interpretation of the additional terms in the equations of motion and their
significance in capturing the modified gravitational dynamics in the graphene
system. Our findings contribute to the understanding of analogue gravity
models and their applications in condensed matter systems.

Analogue gravity; Graphene, $f(R,T)$ theory; Lagrangian density; Modified
Dirac equations

###### pacs:

33.15.Ta

##  I Introduction

The study of analogue gravity extends beyond its intrinsic theoretical
elegance and experimental feasibility. The insights gained from analogue
gravity models contribute to our understanding of fundamental physics,
including the nature of spacetime and the interplay between gravity and
quantum mechanics (Barcelo, ) . Furthermore, the potential applications of
analogue gravity in the field of condensed matter physics may lead to
technological advancements and innovative materials with tailored properties
(Castro, ). In this investigation, our primary focus is on the application of
the analogue gravity framework to graphene, a two-dimensional material
distinguished by its unique electronic (Galeratti, ). Specifically, we adopt
the $f(R,T)$ theory, wherein the Lagrangian density incorporates a function
$f(R,T)$ to describe the modified gravitational effects, as described in the
seminal article (Harko, ). As one see, $f$ is a function dependent on the
Ricci scalar curvature $R$ and the trace of the energy-momentum tensor.

The goal of this study is to derive the Lagrangian and equations of motion for
the quasiparticle field in graphene (Horava, ; Kibis, ) within the $f(R,T)$
formalism and analyze the implications of the additional terms in the
equations of motion. Subsequently, we meticulously scrutinize the implications
arising from the additional terms introduced in the equations of motion under
the $f(R,T)$ framework. This endeavor seeks to advance our understanding of
the intricacies involved in the application of $f(R,T)$ theory to graphene,
shedding light on the nuanced interplay between modified gravitational
dynamics and the electronic characteristics of this two-dimensional material.

We start by defining the $(2+1)D$ Schwarzschild-like metric and calculating
the relevant quantities such as the Ricci scalar $R$ and the trace of the
energy-momentum tensor $T$. We then construct the Lagrangian density by
combining the appropriate terms involving the quasiparticle field, the
effective metric tensor, and the gauge field (Landau, ; Vozmediano, ). The
form of the function $f(R,T)$ is chosen to capture the modified gravitational
dynamics in the system. Finally, we vary the Lagrangian to obtain the
equations of motion (Peskin, ), which provide insights into the dynamics of
the quasiparticles and their interactions with the effective metricand gauge
field.

##  II $f(R,T)$ gravity and the energy-momentum tensor

In the context of the $f(R,T)$ gravitational theory, as described in (Harko,
), and in the general form, the action reads

$S=\int\left[\frac{f(R,T)}{16\pi}+\mathcal{L}_{m}\right]\sqrt{-g}d^{4}x,$ (1)

with $g$ being the determinant of the metric, $f(R,T)$ a function of the Ricci
curvature scalar $R$ and the trace of the energy-momentum tensor $T$. In our
case, the matter lagrangian are $\mathcal{L}_{m}=\rho$, with $\rho$ being the
density.

Varying this action with respect to $g_{\mu\nu}$, we obtain the following
field equations

$f_{R}(R,T)R_{\mu\nu}-\frac{1}{2}f(R,T)g_{\mu\nu}+\left(g_{\mu\nu}\boxempty-\nabla_{\mu}\nabla_{\nu}\right)f_{R}(R,T)$
$=8\pi T_{\mu\nu}+f_{T}(R,T)\left(T_{\mu\nu}-\rho g_{\mu\nu}\right),$ (2)

where $f_{R}(R,T)=\nicefrac{{\partial f}}{{\partial R}}$ and
$f_{T}(R,T)=\nicefrac{{\partial f}}{{\partial T}}$. The general expression for
the covariant derivative of the energy-momentum tensor, as expressed in Harko
is given by

$\displaystyle\nabla^{\mu}T_{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{f_{T}(R,T)}{8\pi-f_{T}(R,T)}\left[\left(T_{\mu\nu}+\rho
g_{\mu\nu}\right)\nabla^{\mu}lnf_{T}(R,T)\right.$ (3)
$\displaystyle\hphantom{\frac{f_{T}(R,T)}{8\pi-
f_{T}(R,T)}}+\left.\nabla^{\mu}\rho
g_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\nabla^{\mu}T\right].$

The covariant derivative of the energy-momentum tensor, unlike Einstein’s
general relativity, are not null. This outcome might be viewed as the
formation or annihilation of matter, implying a lack of energy conservation.
However, if we impose the conservation, i.e., $\nabla^{\mu}T_{\mu\nu}=0$, it
is possible to discover a particular version of the function $f(R,T)$ that
meets this requirement.

##  III Dirac equation in $f(R,T)$ spacetime

Graphene flat sheets can be described as two-dimensional analogs of
relativistic systems for massless fermions. In this context, the Fermi
velocity ($v_{F}$) determines the maximum attainable speed within the graphene
system, analogous to the limiting speed observed in relativistic systems. In
(1+2) dimensional curved spacetime, taking into account the $f(R,T)$ gravity,
the Dirac equations of motion can be obtained from the Lagrangian

$\mathcal{L}=i\hbar v_{F}\int
d^{3}x\sqrt{g}f(R,T)\left(\overline{\psi}\gamma^{\mu}\mathcal{D\psi}\right),$
(4)

where $\mathcal{D}$ represents the covariant derivative and
$\sqrt{g}=\sqrt{(-det(g^{\mu\nu}))}$. The metric tensor $g^{\mu\nu}$ encodes
the effective metric of the system, which depends on the specific graphene
configuration or lattice structure. The Dirac matrices $\gamma_{\mu}$ satisfy
the Clifford Algebra

$\left\\{\gamma^{\mu},\gamma^{\nu}\right\\}=2g^{\mu\nu}\mathds{1}$ (5)

The Euler Lagrange equations result in the following equations of motion for
$\psi$ and $\overline{\psi}$

$i\hbar v_{F}\gamma^{\mu}\mathcal{D}\psi-\partial_{\mu}\left(i\hbar
v_{F}\sqrt{g}f(R,T)\gamma^{\mu}\psi\right)=0$ (6)

$i\hbar v_{F}\gamma^{\mu}\mathcal{D}\overline{\psi}-\partial_{\mu}\left(i\hbar
v_{F}\sqrt{g}f(R,T)\gamma^{\mu}\overline{\psi}\right)=0$ (7)

These equations describe the dynamics of the field $\psi$ and its adjoint
$\overline{\psi}$. The field $\psi(x)$ represents the quasiparticle field in
graphene. It is important to note that extra terms are determined by the exact
form of the covariant derivative $\mathcal{D}$ and the gamma matrix
representation used. The covariant derivative
$\mathcal{D}=\partial_{\nu}+eA_{\nu}$ includes the coupling of the
quasiparticles to the effective gauge field $A_{\nu}$, which describes the
modified gravitational effects in the $f(R,T)$ theory. Note that if $f(R,T)=0$
we recover the equations of motion for the massless Dirac fermions.

Overall, this Lagrangian captures the dynamics of the quasiparticles in the
graphene system within the analogue gravity framework, accounting for the
modified gravitational effects described by the function $f(R,T)$ and their
interaction with the effective metric and gauge field. The curvature effects
can be characterized for the extra terms in Eqs. (11) and (12), coupling the
Dirac spinors to the metric $g_{\mu\nu}$ and the form of $f(R,T)$ function,
which can be can be treated as ansatz.

For a more simple model, we can neglect some possible topological defects in
the lattice and consider an ideal homogeneous distribution. The extra terms in
the equations of motion beyond the standard Dirac equation can be interpreted
as the additional effects arising from the modified gravitational dynamics
described by the function $f(R,T)$ in the analogue gravity model.

Let’s examine the additional terms in the equations of motion: For the
quasiparticle field $\psi(x)$ equation:

$\partial_{\mu}\left(i\hbar v_{F}\sqrt{g}f(R,T)\gamma^{\mu}\psi\right)=0,$ (8)

and for the adjoint quasiparticle field $\overline{\psi}(x)$ equation:

$\partial_{\mu}\left(i\hbar
v_{F}\sqrt{g}f(R,T)\gamma^{\mu}\overline{\psi}\right)=0$ (9)

These additional terms can have various interpretations depending on the
specific form of $f(R,T)$ and the chosen gauge field. They could represent
modifications to the dispersion relations, energy levels, or coupling
strengths of the quasiparticles due to the presence of the modified gravity
sector. They might also introduce new interactions or affect the transport
properties of the quasiparticles in the graphene system.

The precise interpretation of these extra terms would require a detailed
analysis and consideration of the specific properties of the chosen $f(R,T)$
function, the gauge field, and their impact on the graphene system within the
analogue gravity framework.

##  IV Massless fermions in $f(R,T)$ spacetime

Now, in another way to investigate we consider that the electrons, i.e., de
Dirac masslesss fermions exists in a (2+1) dimensional perfect fluid. This
perfect fluid is the homogeneous and isotropic honeycomb structure of a
graphene sheet.

In this case, the action is given by

$S=\int\sqrt{-g}d^{3}x\left[\frac{f(R,T)}{2\kappa^{2}}+\mathcal{L}_{m}\right],$
(10)

where $\kappa^{2}=8\pi G$, and again $f(R,T)$ is the function of Ricci scalar
curvature $R$ and $T$ is the trace of the energy-momentum tensor of a (2+1)D
perfect fluid, given by

$T=T_{\mu}^{\mu}=(-\rho+2p)$ (11)

with $\rho$and $p$ being the energy density and the pressure repsctively.

Assuming that $\mathcal{L}_{m}$ is given by the massless Dirac fermions
Lagrangian, i.e.

$\mathcal{L}_{m}=i\overline{\psi}\gamma^{\mu}\partial_{\mu}\psi.$ (12)

Varying this action (13) in relation to the metric one obtains the general
form of (2+1) dimensional field equations given by the Eq.(2). Then, let´s
assume the perfect fluid describes the graphene sheet in wich we have the
electrons represented by the matter Lagrangian showed in Eq. (15). Also
assuming the linear form for $f(R,T)$ function as

$f(R,T)=R+\alpha T,$ (13)

with $\alpha$ being a positive constant. With these above assumptions, the
field equation becomes

$R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}=\kappa^{2}T_{\mu\nu}+\alpha
T_{\mu\nu}-\frac{\alpha
T}{2}g_{\mu\nu}-\alpha\left(i\overline{\psi}\gamma^{\mu}\partial_{\mu}\psi\right)g_{\mu\nu},$
(14)

where, on the left side of this equation we have the Einstein tensor in (2+1)
dimensions, and in the right hand side we have the terms of correction in the
energy-momentum tensor and the massless Dirac fermions term coupled with the
metric tensor $g_{\mu\nu}$. Note that if $\alpha=0$, we recover the general
relativity field equations.

##  V Torsion effects

Torsion is a geometric property that describes the twisting of a surface. In a
$(2+1)D$ sheet, torsion can arise due to curvature or strain. Torsion has a
significant impact on the electronic properties of $(2+1)D$ sheets. In the
presence of torsion, the Dirac equation is modified to include an additional
term:

$\frac{i}{4}\gamma^{\mu\nu}T_{\mu\nu}\psi$ (15)

where $T_{\mu\nu}$ is the torsion tensor and $\gamma^{\mu\nu}$ are the Dirac
matrices.

Torsion is a geometric property related to the amount of twisting or warping
that occurs in a material. In the context of a graphene-like material such as
Graphenylene, torsion can affect the electronic and mechanical properties of
the material. In particular, the torsion of a Graphene like sheet can lead to
a modification of its band structure, which determines the electronic
properties of the material. Torsion can also affect the mechanical behavior of
the nanosatructure, by altering its stiffness and strength.

Furthermore, the presence of torsion can lead to the emergence of new physical
phenomena, such as the generation of gauge fields that can affect the
electronic properties of the material. Therefore, understanding the effects of
torsion on nanostructures is important for both fundamental research and
practical applications.

To derive the equations of motion for the fermions and then the equation
expressing the electronic density in the graphene sheet considering the
effects of torsion, we start with a $(2+1)D$ metric as follows

$ds^{2}=-g(r)dt^{2}+h(r)dr^{2}+r^{2}d\theta^{2},$ (16)

where $g(r)$ and $h(r)$ are radial functions, as we can find in wormholes and
black holes formalism111For ordinary Schwarzschild metric in four dimensions
we have $g(r)=\left(1-\frac{2GM}{r}\right)$ and
$h(r)=\left(1-\frac{2GM}{r}\right)^{-1}$. Here $G$ is the four-dimensional
Newton constant and $M$ is the mass of the black hole.. The Euler-Lagrange
equations of motion, including the torsion term, result in

$\displaystyle i\hbar v_{F}\gamma^{\mu}\mathcal{D}_{\mu}\psi$
$\displaystyle-\partial_{\mu}(i\hbar v_{F}\sqrt{g}f(R,T)\gamma^{\mu}\psi)$
(17) $\displaystyle-i\hbar
v_{F}\sqrt{g}\partial_{\mu}f(R,T)\gamma^{\mu}\psi=0$

$\displaystyle i\hbar v_{F}\gamma^{\mu}\mathcal{D}_{\mu}\overline{\psi}$
$\displaystyle-\partial_{\mu}(i\hbar
v_{F}\sqrt{g}f(R,T)\gamma^{\mu}\overline{\psi})$ (18) $\displaystyle-i\hbar
v_{F}\sqrt{g}\partial_{\mu}f(R,T)\gamma^{\mu}\overline{\psi}=0$

These equations describe the dynamics of the fermions in the graphene sheet
within the framework of $f(R,T)$ gravity and considering the given metric.

Now, to express the electronic density in the graphene sheet, we need to
consider the conservation equation. The electronic density $n$ can be obtained
from the equation of continuity, which in curved spacetime takes the form:

$\nabla_{\mu}(\rho u^{\mu})=0$ (19)

where $\rho$ is the charge density, $u^{\mu}$ is the four-velocity of the
charge carriers, and $\nabla_{\mu}$ represents the covariant derivative. In
the case of graphene, we can assume that the charge carriers move with
negligible velocity compared to the speed of light. Therefore, we can consider
$u^{\mu}=(1,0,0,0)$ in suitable coordinates. Using the definition of the
charge density $n=\rho/e$, where $e$ is the elementary charge, and taking into
account that $u^{\mu}=(1,0,0,0)$, the equation of continuity simplifies to:

$\partial_{\mu}(\sqrt{-g}n)=0$ (20)

Expanding the covariant derivative, this becomes:

$\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}n^{\mu})=0$ (21)

$\partial_{\mu}n^{\mu}+n^{\mu}\Gamma_{\mu\sigma}^{\sigma}=0$ (22)

where $\Gamma_{\mu\sigma}^{\sigma}$ are the Christoffel symbols associated
with the metric.

This equation expresses the conservation of charge in the graphene sheet,
accounting for the effects of curvature and torsion. The exact form of
$n^{\mu}$ would depend on the specific properties of the charge carriers in
graphene and the details of the underlying theory.

##  VI Results

To visualize the effects of torsion on a $(2+1)D$ sheet, we generate a 3D
image using Python and Mayavi. We first generate a $(2+1)D$ sheet with torsion
effects using the code we developed earlier. We then use Mayavi to plot the
sheet and apply a color map to represent the electronic density as showed in
the figure 1.

Figure 1: The torsion effect in a $(2+1)D$ sheet, showing the electron density
in specific regions. This simulation was generated using a Phyton code in
which the parameters can be changed.

In this article, we investigate the influence of torsion on the electronic
properties of graphenylene sheets using the Palatini formalism of general
relativity. Torsion is a measure of the twist in the material, which can
result from the deformation of the honeycomb lattice structure. We find that
the torsion can significantly alter the electronic band structure of
graphenylene sheets, leading to the emergence of new electronic states.

To illustrate the effects of torsion on the electronic properties of
graphenylene sheets, we generate 3D images of the sheets using a Python code
that incorporates torsion. The code uses the Palatini formalism to derive the
field equations for a curved 2D space and then applies the torsion metric to
create a graphenylene sheet with torsion effects.

Our analysis shows that the torsion-induced modifications to the electronic
properties of graphenylene sheets can be understood in terms of an effective
gauge field. The gauge field is proportional to the torsion, and its presence
leads to a topological phase transition in the electronic band structure. This
effect can be understood as an analogue of the Aharonov-Bohm effect in
condensed matter systems.

In resume, the most important electronic property that changes in a
Graphenylene sheet when torsion is introduced is the band gap. In a pristine
Graphenylene sheet, the electronic band structure is gapless, with the valence
and conduction bands touching at the Dirac point, resulting in high electron
mobility and high conductivity. However, when torsion is introduced, it breaks
the sublattice symmetry and leads to the opening of a band gap. The magnitude
of the band gap depends on the magnitude and direction of the torsion angle,
as well as the size and shape of the Graphenylene sheet.

Therefore, by applying torsion to a Graphenylene sheet, we can effectively
tune its electronic properties and create a semiconducting material with a
tunable band gap. This makes Graphenylene a promising material for a wide
range of electronic applications, including transistors, solar cells, and
sensors.

Our results demonstrate the importance of considering torsion effects in the
study of graphenylene sheets and provide a new avenue for exploring the
intriguing electronic properties of these materials. We hope that our work
will inspire further research in the field of analogue gravity in condensed
matter systems.

## References

  * (1) Barcelo, Carlos, Stefano Liberati, and Matt Visser. "Analogue gravity." Living reviews in relativity 14 (2011): 1-159.
  * (2) Castro Neto, A. H., Guinea, F., Peres, N. M. R., Novoselov, K. S., Geim, A. K. (2009). The electronic properties of graphene. Reviews of modern physics, 81(1), 109.
  * (3) Galeratti, Antonio. Graphene, Dirac equation and analogue gravity. Physica Scripta, v. 97, n. 6, p. 064005, 2022.
  * (4) Harko, T., Lobo, F. S., Nojiri, S. I., and Odintsov, S. D. (2011). f (R, T) gravity. Physical Review D, 84(2), 024020.
  * (5) F. D. M. (1988). Model for a quantum Hall effect without Landau levels: Condensed-matter realization of the "parity anomaly". Physical Review Letters, 61(18), 2015.
  * (6) Horava, P. (2009). Quantum gravity at a Lifshitz point. Physical review D, 79(8), 084008.
  * (7) Kibis, O. V., Lozovik, Y. E. (2015). Graphene and graphene-like materials as platforms for the study of condensed-matter analogues of cosmological phenomena. Physics-Uspekhi, 58(10), 981.
  * (8) Landau, L. D., Lifshitz, E. M. (1975). The classical theory of fields (Vol. 2). Butterworth-Heinemann.
  * (9) Vozmediano, M. A. H., Katsnelson, M. I., Guinea, F. (2010). Gauge fields in graphene. Physics Reports, 496(4-5), 109-148.
  * (10) Peskin, Michael E. An introduction to quantum field theory. CRC press, 2018.