Forecasting Intraday Power Output by a Set of PV Systems using Recurrent
  Neural Networks and Physical Covariates
Accurate intraday forecasts of the power output by PhotoVoltaic (PV) systems
are critical to improve the operation of energy distribution grids. We describe
a neural autoregressive model that aims to perform such intraday forecasts. We
build upon a physical, deterministic PV performance model, the output of which
is used as covariates in the context of the neural model. In addition, our
application data relates to a geographically distributed set of PV systems. We
address all PV sites with a single neural model, which embeds the information
about the PV site in specific covariates. We use a scale-free approach which
relies on the explicit modeling of seasonal effects. Our proposal repurposes a
model initially used in the retail sector and discloses a novel truncated
Gaussian output distribution. An ablation study and a comparison to alternative
architectures from the literature shows that the components in the best
performing proposed model variant work synergistically to reach a skill score
of 15.72% with respect to the physical model, used as a baseline.
Pierrick Bruneau, David Fiorelli, Christian Braun, Daniel Koster
Hybrid-Physical Probabilistic Forecasting for a Set of Photovoltaic Systems using Recurrent Neural Networks

###### Abstract

Accurate intra-day forecasts of the power output by PhotoVoltaic (PV) systems are critical to improve the operation of energy distribution grids. We describe a hybrid-physical model, which aims at improving deterministic intra-day forecasts, issued by a PV performance model fed by Numerical Weather Predictions (NWP), by using them as covariates in the context of an autoregressive recurrent neural model. Our proposal repurposes a neural model initially used in the retail sector, and discloses a novel truncated Gaussian output distribution. We experimentally compare many model variants to alternatives from the literature, and an ablation study shows that the components in the best performing variant work synergistically to reach a skill score of 7.54% with respect to the NWP-driven PV performance model baseline.

## 1 Introduction

Grids of PV systems have become an inevitable component in the modern and future energy distribution systems. However, due to weather conditions, the magnitude of PV power production is fluctuating, while the supply to consumers requires to be adapted to the demand at each point in time. Distribution system operators (DSOs) have increasing and specific requirements for PV power forecasts. Indeed, fluctuating renewables could cause operational issues (e.g. grid congestion), which call for active grid operation. In this context, _intra-day_ forecasts (i.e., for the whole day to come at a fixed time in day) of PV power are critical in view to facilitate operations. Also, many forecasting models issue point forecasts, but hardly characterize the uncertainty attached to their forecasts, when such information can be critical for a DSO in order to quantify and mitigate risks in an optimal way.

In [10], PV power production is forecasted using a deterministic PV performance model, which involves regional solar irradiance forecasts issued by a Numerical Weather Prediction (NWP) service as inputs. The underlying hypothesis is that solar irradiance is fairly smooth over limited regional areas, and the production curve specific to a PV system will be mainly influenced by how it converts this solar energy to PV power according to its specifications. In [11], authors of the present paper introduced a model which performs intra-day probabilistic forecasts of PV power production. It combines the PV
performance model above with a model based on Long-Short-Term Memory (LSTM) cells [14]. This kind of combination of a model based on a set of physical equations to a statistical model is coined as _hybrid-physical_ approaches [1]. For training and evaluation, it uses real data provided by Electrois, a DSO in Luxembourg. Results show that this new model improves the baseline performance, while coping with local effects such as PV system shading. The former paper rather targets solar energy specialists, with little details unvealed about how the Machine Learning model acting as cornerstone to the approach has been designed and trained. The present paper aims at filling this gap, by providing entirely new material focusing on this complementary view. Specifically, the purpose of the present paper is to focus on neural time series forecasting aspects in the context of this application.

The specific contributions of the present work mainly focus on the design of a model architecture and a training procedure which meets the operational needs of a DSO. Our proposal is based on an existing LSTM-based model [15], which we present in a concise and effective way in order to make this work self-contained. In addition, we design a novel truncated Gaussian output component, which we plug in to the LSTM-based model. In Section 2, we give a structured survey of the related work which positions the problem faced by a DSO, and motivates which existing work could be reused or repurposed to suit our needs. After describing our model proposal in Section 3, we provide a thorough experimental evaluation in Section 4. Several variants of our model proposal are compared to alternative models from the literature, and an ablation study allows to emphasize the specific contribution of each of its components. Finally we recall some qualitative results to underline how local effects, tainting the PV performance model, are mitigated using our approach.

## 2 Related Work

In Section 2.1, we review seminal PV power forecasting methods. Then in Section 2.2, we survey time series forecasting as addressed in the literature on Machine Learning (ML), from the perspective of their repurposing potential to the PV power forecasting application. Section 2.3 reviews existing hybrid-physical models which relate the most closely to our proposal. Finally, Section 2.4 focuses on the peculiarities in terms of forecasting structure and validation which come with ML approaches to time series forecasting.

### PV Power Forecasting

Most approaches in PV power forecasting model the conversion chain of solar irradiance to electrical power in a PV system. Thus they follow a two-step approach: first, forecasting the solar irradiance, then converting this irradiance to PV power forecasts [1, 1]. The most common way to forecast solar irradiance relies on NWP systems such as the European Centre for Medium-Range Weather Forecasts (ECWMF) Ensemble Prediction System [2]. Everyday, it issues hourly regional forecasts for a range of meteorological variables (including solar irradiance) on the 10 days to come, thus including the intra-day range. Intra-day may be defined as up to 6h forecast horizons in the literature [1]. In this paper, we deviate from this defini
tion by considering intra-day as the next 24h, starting at midnight of the same day. In view to improve NWP forecasts, or to avoid having to rely on such systems, solar irradiance can also be forecasted using statistical methods and ML models. The simplest include persistence models, which are often adjusted using clear sky models [11]. [20] also review various ML techniques which have been employed to this purpose, e.g., AutoRegressive (AR) models, Feed-Forward Networks (FFN) and k-Nearest Neighbors. In this range of contributions, [12] address the intra-day hourly prediction of solar irradiance using an ensemble of FFN models. Specifically, they implement rolling forecasts by specializing each model to a current time and a prediction horizon.

PV power forecasting is reviewed in detail by [1]. In this landscape, several approaches aim at modelling directly the series of PV power values, without having to rely on solar irradiance forecasts. [10] propose short-term forecasts (\(<\)30mn) which exploit cross-correlation of PV measurements in a grid of PV systems. They hypothesize that clouds casting over a given PV system have lagged influence on other systems downwind. They optimize associated time lag and cloud motion vector. [1] also consider a spatially distributed set of PV panels. They directly use PV power values, without converting proxy information such as solar irradiance or cloud density. Similarly to [10], they focus on correlations among stations to help accounting for intermittency due to clouds. They report intra-day forecasts are useful for an energy trading strategy, while hour-ahead and intra-hour serve for managing demand response.

[14] present AR approaches to PV power forecasting. They focus on forecasting one and two hours ahead, where NWP models tend to under-perform. Several models are compared, among which are persistence, linear models such as AutoRegressive Integrated Moving Average (ARIMA) [15], and FFN. They found out that FFN perform the best, with improvements brought by the optimization of FFN parameters, input selection and structure using a Genetic Algorithm. They conjecture that binning data according to the associated season, and learning per-bin models should improve forecasting ability overall, even if this approach has been recently criticized [16].

### ML approaches for time series forecasting

Among other related work, Section 2.1 surveyed some contributions which involved ML methods in view to forecast solar irradiance and PV power production. In this section, we generalize this view, by surveying recent work in time series forecasting at large. Methods in this section were generally not applied to the application context considered in the present paper, but could be repurposed _a priori_. Besides neural and ARIMA models, seminal ways to forecast time series include the Croston method, which is an exponential smoothing model dedicated to intermittent demand forecasting [17]. It is notably used as baseline method in [16], along with the Innovation State-Space Model (ISSM) [18].

Modern, so-called _deep_ neural network architectures exploit the structure of data - sequential, in the case of the Long-Short-Term Memory (LSTM) [19] and the Gated Recurrent Unit (GRU) [12], 2D or 3D, in the case of convolutional networks. Even though recurrent models such as the LSTM would appear as outdated, they are still popular in the recent literature thanks to improvements to training and inference procedures carried by modern toolboxes
such as Tensorflow [1] or MXNet [15], as well as the encoder-decoder mechanism, in which a context sequence is encoded and conditions the prediction of the target sequence. It has been initially codified for the GRU, and transferred to the LSTM, leading to continued usage in recent contributions [16, 17]. These models contrast with the seminal Feed-Forward Network (FFN), in which all layers are fully connected to previous and next layers, up to the activation layer [11].

Salinas et al. propose DeepAR, which implements flexible forecasting for univariate time series [17]. Formally, it defines a _context_ interval (a chunk of past values) and a _prediction_ interval (the set of values to be forecasted), and the model is optimized end-to-end w.r.t. prediction interval forecasts. This contrasts with linear models such as ARIMA, which are optimized for one time step ahead forecasts. Also, instead of point forecasts, DeepAR predicts model parameters, which can be used to compute sample paths, and empirical quantiles, which can be highly valuable in our context. In this case, a family of probability distributions has to be chosen so as to fit the time series at hand. The model was initially aimed at retail business applications, but it can be adapted to other types of data just by changing the family of the output probability distribution. It is based on the LSTM model architecture. The model supports the adjunction of covariates, i.e., time series which are available for both context and prediction interval at forecast time. By repeating the same value for the whole intervals, static covariates are also supported. Such flexibility meets all the requirements of our hybrid-physical approach (i.e., NWP covariates and PV system descriptive features).

In retail applications, input data may have highly variable magnitude (e.g., depending on item popularity or time in the year). The authors observe an approximate power law between item magnitude and frequency (i.e., values are less likely as they get large, and reciprocally). They claim that grouping items to learn group-specific models or performing group-specific normalizations, as previously done in the solar and PV power forecasting literature [14, 15] are not good strategies in such case. They propose a simple alternative scheme, where samples are scaled by an item dependent factor computed using the context interval.

[18] propose a time series classification model inspired by Inception-v4 [16]. Basically, they are transferring the multi-scale pattern extraction capability of convolutional neural networks to 1D data such as time series. A convolutional encoder is tested by Wen et al. in the context of their multi-horizon quantile forecaster [21]. Instead of forecasting probabilistic model parameters, this model directly forecasts quantiles in a non-parametric fashion. However, it suffers from the quantile crossing problem: forecasted values may have ranks inconsistent with the quantile they are attached too.

[11] is another alternative to DeepAR. Similarly to [21], it does not rely on probability distribution outputs, and implements conditional quantile functions using regression splines instead. Spline parameters are fit using a neural network directly minimizing the Continuous Ranked Probability Score (CRPS), which is then used as a loss function. This results in a more flexible output distribution, and an alternative to other flexible schemes (e.g. mixture of distributions in the context of [17]). However, it currently1
lacks a publicly available implementation.

Multivariate forecasting consists in modelling and forecasting multiple time series simultaneously, by contrast to univariate forecasting. The seminal way to achieve this, is with the Vector AutoRegression (VAR) model, which is an extension of the linear autoregressive model to multiple variables [14]. As this model has hard time dealing with many variables (e.g., items in the retail domain), neural network-based models such as DeepVAR were designed as an alternative [10]. It can be thought of as a multivariate extension to [11]. DeepVAR models interaction between time series, e.g., as resulting from causality or combined effects. It uses a Copula model, which models interactions between time series, and elegantly copes with time series of varying magnitude, alleviating the need for an explicit scaling mechanism. A single multivariate state variable underlying a LSTM model is used for all time series. Empirical cumulative distributions serve to compute sample paths and quantiles. Only static covariates (i.e., constant in the context and prediction intervals) were considered in this paper. They define a low-rank parametrization, which opens the possibility to deal with a very large number of time series.

Some prior work involved deep learning in the context of PV power forecasting. For example, [1] consider one hour ahead forecasts (instead of intra-day as aimed at in this paper). They used a single LSTM layer without the encoder-decoder mechanism. Also, they consider point forecasts. Data for two PV sites is used for the experiments, with roughly the same power magnitude for both sites (approx. 3.5kW). Models are trained for each site separately. Alternatively, in this paper we address all sites with a single model in a scale free approach, addressing an arbitrary number of sites with little to no model size overhead. Finally, the locations associated to the datasets are subject to a dry climate, which is simpler to predict [2]. Our application testbed is a temperate area, subject to frequent and abrupt changes on a daily basis, therefore much more challenging to predict.

[14] also address PV power forecasting, by decorrelating scale free forecasts using LSTM, from seasonal effects modelled separately using time correlation features and partial daily pattern prediction. However, they focus on forecasts aggregated at a daily scale, when we consider hourly data in this paper. In addition, our approach is end-to-end, with no independent modelling of seasonal effects.

### Hybrid-physical approaches

[1] focus on wind speed and power forecasting. They consider hourly forecasts for 72 hours ahead, using Numerical Weather Prediction (NWP) forecasts as additional inputs, thus proposing an early combination of observation data with NWP covariates. The recurrent model used then (diagonal recurrent neural networks [12]) was superseded by alternative models such as LSTM [10] and GRU [15] in the recent literature as seen in the previous section. Another early work is proposed by [13], who combine NWP covariates and neural networks for hourly and daily solar irradiance forecasting.

[17] present a regional PV power forecasting system, with hourly resolution up to 72h ahead. Their approach combines clustering and numerical optimization, and it is compared to regression methods such as ElasticNet [18], SARIMAX [19], or Random Forests [19]. The CRPS metric is used
for evaluation. Their approach is not autoregressive, rather they directly predict future PV power from solar irradiance and temperature forecasts obtained from a proprietary system which refines NWP forecasts according to local conditions. Alternatively, our approach tries to combine the benefits of using NWP forecasts with an autoregressive model of the PV power observations.

### Forecast structure and validation

Figure 1 distinguishes _regular_ forecasts from _rolling_ forecasts, which are the two main strategies to consider for extracting fixed sized blocks from time series. For simplicity, the figure considers hourly forecasts and 24-hour context and prediction intervals, but the definition is straightforward to generalize. In brief, two consecutive regular forecasts are offset by the size of the prediction interval, when rolling forecasts are offset by the frequency of the time series (hourly, in Figure 1). In other words, with 24-hour prediction and context intervals, regular forecasts happen on the same time every day, whereas rolling forecasts are issued every hour for the whole prediction interval. In this process, forecasts beyond the next hour are refreshed every hour.

Works such as [11] consider regular forecasts, as the forecast time is tied to availability of NWP data. As we use similar NWP covariates, regular forecast are a requirement in our work too. Alternatively, [12] address rolling forecasts by having a distinct model for each possible starting time in the day. Let us note that some models (e.g., [13]) allow to encode seasonal information on the predicted time steps (e.g. hour in day, day in week) as covariates. Therefore, they can be used indistinctively with regular and rolling forecasts, provided an adapted training set is available.

Time series forecasting models are typically evaluated using a variant of the Root-Mean-Square Error (RMSE) metric. When quantiles can be computed, the Continuous Ranked Probability Score (CRPS) rates the compatibility of an observation with a set of quantiles [10]. This metric is generally used for evaluating models which output forecasting quantiles [10, 11, 12].

[12] discuss the problem of cross-validation, and more generally validation, in the context of time series forecasting. Original formulations of cross-validation methods often assume that data items are independent. They cannot be used out of the box with time series, as the sequential structure of the latter invalidates the underlying hypotheses. They recognize that the

Figure 1: Distinction between regular and rolling forecasts. \(t_{0}\) denotes the present time step in the context of a given data sample.
seminal way to validate models with time series is to train a model using the first \(t_{0}\) values, and validate using the last \(T-t_{0}\) values. However, this way is hardly compatible with cross-validation schemes, and yields weak test error estimations. In virtue of the bias-variance tradeoff [1], this issue has moderate impact on models with strong bias such as linear models. However, over-parametrized models, such as most neural models presented in Section 2.2 (e.g. [12, 12, 13, 14]) can be significantly affected, and exhibit a strong tendency to overfit (even if recent theory shows that with careful consideration this problem can be correctly addressed [15]). For mitigation, [1] recommend blocked cross-validation, in which time series segments of size \(\tau<<T\) are used for independent training and validation for model selection and test error computation. As we also use deep learning as a building block in our approach, we carefully consider these recommendations in our experimental design (see Section 4).

## 3 Model Description

The survey of related works in Section 2 led us to choose DeepAR [12] as a framework to develop our hybrid-physical implementation. We adapted the official implementation of the model [1] to suit our needs. Another model which offers the relevant flexibility as well as a public implementation is the model by [17]. We will compare to this model in our experimental section.

Alternatively, PV sites could have been considered as dimensions in a multivariate forecasting problem, thus possibly forecasting all sites at a given time at once using DeepVAR [12]. However, its limitation to static covariates prevents us from implementing the projected hybrid-physical approach. Also, it is unclear how new PV sites, as well as missing values, which are pretty common as PV sites may witness independent breakdowns and interruptions of measurements or data communication, can be handled with such multivariate modelling. Alternatively, we choose to model the PV site using covariates.

### DeepAR model

In the remainder of the paper, for clarity of the derivations, scalar variables are represented in normal font, vector variables in bold font, and matrix variables in capital bold font. This section 3.1 essentially paraphrases [12], but ensures the present paper is self-contained, while introducing the necessary formalism.

Let us assume we have a data set of \(N\) univariate time series, each with fixed size \(T\). Each observed time series in \(\{\mathbf{z}_{n}\}_{n\in 1,\ldots,N}\) may relate to an item in store in the retail context, or to a distinct PV system in the context addressed in this paper. \(t_{0}\in[1\ldots T]\) denotes the present time, i.e. the latest time point for which we assume \(z_{n,t}\) is known when issuing the forecast. \([1\ldots t_{0}]\) is then the _condition_ interval, and \([t_{0}+1\ldots T]\) is the _prediction_ interval. The goal of the model is to forecast \(\mathbf{z}_{n,t_{0}+1:T}=[z_{n,t_{0}+1},\ldots,z_{n,T}]\) with the knowledge of \(\mathbf{z}_{n,0:t_{0}}=[z_{n,0},\ldots,z_{n,t_{0}}]\). We also consider a set of covariates \(\mathbf{X}_{n,1:T}=[\mathbf{x}_{n,1},\ldots,\mathbf{x}_{n,T}]\) which are known for \(t\in 1,\ldots,T\) at time \(t_{0}\).

In this context, the model is defined by the following product of likelihood factors, also summarized in Figure 2:
\[Q_{\Theta} =\prod_{n=1}^{N}\prod_{t=t_{0}+1}^{T}q_{\Theta}(z_{n,t}|\mathbf{z}_{n,1:t-1},\mathbf{X}_{n,1:T})\] \[=\prod_{n=1}^{N}\prod_{t=t_{0}+1}^{T}p(z_{n,t}|\theta(\mathbf{h}_{ n,t},\Theta)) \tag{1}\]

The model is both autoregressive and recurrent, as state variable:

\[\mathbf{h}_{n,t}=\Theta(\mathbf{h}_{n,t-1},z_{n,t-1},\mathbf{x}_{n,t}) \tag{2}\]

is obtained from LSTM model \(\Theta\) in which the state variable and observation of the previous time step are both reinjected. The model also depends on parametrized function \(\theta\), which learns the mapping between the state variable \(\mathbf{h}\) and parameters of the probability distribution \(p\).

In effect, as seen in Figure 2, during training time observations are injected as \(z_{n,t-1}\) in Equation (2). However at test time actual observations are not available for the prediction interval, so we sample \(\tilde{z}_{n,t}\sim p\), and inject them as proxy observations. Doing so yields sample paths, which can be repeated and serve to compute empirical quantiles in the prediction interval, instead of simple point estimates. In this paper, when point estimates \(\hat{z}_{n,t}\) are needed, we take them as the empirical median of a set of sample paths. In Figure 2, we can see that the LSTM _encodes_ the context interval into \(\mathbf{h}\), which is then _decoded_ for the prediction interval. The same LSTM model \(\Theta\) is used for encoding and decoding.

The negative log of expression (1) is used as the loss function for training all parameters in the model in an end-to-end fashion. The form of function \(\theta\) depends on the probabilistic model in expression (1): for example, if \(p\) is chosen as a Gaussian, appropriate functions would be:

Figure 2: Illustration of the DeepAR model. Observed variables are represented as shaded boxes, and latent variables as blank boxes. For the context interval, \(z\) variables are always known. For the prediction interval, the model behaves differently at training and test time. At test time, \(\tilde{z}\) variables are sampled according to \(p\), forming sample paths. Plain lines represent dependencies between random variables, and the dashed line highlights the reinjected sample.
\[\theta_{\mu}(\mathbf{h}_{n,t}) =\mathbf{w}_{\mu}\mathbf{h}_{n,t}+b_{\mu}\] \[\theta_{\sigma}(\mathbf{h}_{n,t}) =\log(1+\exp(\mathbf{w}_{\sigma}\mathbf{h}_{n,t}+b_{\sigma})) \tag{3}\]

We note that the softplus function in (3) ensures \(\sigma\) is mapped as a positive real number. Among possible probabilistic models and mapping functions, the official DeepAR implementation [1], used in the experiments for this paper, features Gaussian, Student, negative binomial, and mixture distributions. The mixture distribution composes several distributions from the same nature using mixture weights, which have their dedicated \(\theta\) function.

### Positive Gaussian likelihood model

As PV power measurements are bound to be non-negative real numbers, a contribution of this paper is to allow for the Gaussian distribution to be truncated from below at 0, referred to as the _positive Gaussian_ distribution in the remainder of this paper. Formally this yields:

\[p(z_{n,t}|\theta_{\mu},\theta_{\sigma})=\frac{1}{\sigma\sqrt{2\pi}}\frac{\exp \left(-\frac{1}{2}\frac{(z_{n,t}-\theta_{\mu})^{2}}{\theta_{\sigma}{}^{2}} \right)}{1-\Phi(-\frac{\theta_{\mu}}{\theta_{\sigma}})} \tag{4}\]

With \(\Phi\) the cumulative distribution function of the standard Gaussian (i.e., with mean 0 and standard deviation 1). Besides adapting the loss function (see Equation (1)) to this new probability distribution function, the same \(\theta_{\sigma}\) function as the Gaussian distribution can be used. To make sure the range of \(\theta_{\mu}\) is also positive, for the positive Gaussian we use:

\[\theta_{\mu}(\mathbf{h}_{n,t})=\log(1+\exp(\mathbf{w}_{\mu}\mathbf{h}_{n,t}+b_ {\mu}))\]

From an application of the Smirnov transformation [1] to the case at hand, samples from a positive Gaussian distribution can be obtained as:

\[\tilde{z}=\Phi^{-1}\Bigg{(}\Phi\big{(}-\frac{\theta_{\mu}}{\theta_{\sigma}} \big{)}+\tilde{u}\bigg{(}1-\Phi\big{(}-\frac{\theta_{\mu}}{\theta_{\sigma}} \big{)}\bigg{)}\Bigg{)}\theta_{\sigma}+\theta_{\mu} \tag{5}\]

where \(\tilde{u}\) is a uniform sample in \([0,1]\).

## 4 Experiments

### Data

Section 3 presented the forecasting model underlying our experiments in general terms, but here we recall that we focus on a specific application and its peculiarities: forecasting the power output of a set of PV systems.

The variable to forecast (\(z\) in Section 3) is the average power output of a PV system during the last hour in Watts. As hypothesized in Section 2.4, it is thus a hourly time series. For our experiments, we used data recorded by 119 PV systems located in Luxembourg between 01/01/2020 and 31/12/2021.
They are dispatched in a relatively small (4 \(\times\) 4 km) area. These PV systems are managed by Electris, a DSO in Luxembourg which collaborated with the authors of this paper in the context of a funded research project. Besides PV power measurements, each time step is associated to intra-day, day-ahead, and 2 days ahead predictions by the physical PV performance model described in [13], which uses ECMWF NWP solar irradiance forecasts as input (referred to as _24h_, _48h_, and _72h_ NWP forecasts, respectively, in the remainder of the paper). We use these tier forecasts as covariates (**X** in Section 3), as they will be available beforehand for the prediction interval.

The model also supports the adjunction of _static_ covariates, which are constant for a given time series. Relating to Section 3, we note that this simply amounts to set associated \(\mathbf{x}_{n,1:T}\) to a constant. In the context of the present work, we consider a _system ID_ categorical feature, which is simply the system ID converted to a categorical feature with 119 modalities. We also consider _system description_ continuous features. Among the set of descriptors provided by system vendors and characteristics of their setup, we retain the following features as they are expected to influence PV power curves and magnitude: the _exposition_ of the system (in degrees), its _inclination_ (in degrees), its nominal _power_ (in Watts) and its _calibration factor_ (unitless, tied to the system on-site setup). As the DeepAR implementation expects normally distributed features, we standardize features so that they have zero mean and unit standard deviation.

As the nomimal power of our PV systems varies over a large range (from 1.4kW to 247kW), a scaling scheme is necessary to properly handle measured values. As mentioned in Section 2, we address this as implemented in DeepAR by dividing all measurements in a given sample by \(\frac{1}{t_{0}}\sum_{1}^{t_{0}}|z_{t}|\). Also, as NWP covariates are expected to be distributed similarly to their associated measurements, they are normalized likewise.

### Loss function and metrics

The sum of negative log-likelihoods of observations \(z_{n,t}\) (on top of Figure 2) is used as a loss function to fit all model parameters \(\{\mathbf{w}_{\mu},\mathbf{w}_{\sigma},\Theta\}\) in an end-to-end fashion. As commonly done in the literature (see Section 2.2), we use a fixed size for the context and prediction intervals in our experiments. As we are interested in intra-day regular forecasts (see Section 2.4), with hourly data this means that the prediction interval has size 24. The midnight run of the ECMWF NWP forecasts for the 3 days to come are broadcasted each day early in the morning, but before the sunrise, so while PV power is still obviously zero. We assume also the associated covariates are instantly available then. For simplicity, we thus choose midnight as the reference time in day for the regular forecasts (i.e., \(t_{0}\) in Section 3). In this context, 24h, 48h and 72h NWP covariates associated to predicted time step \(t_{0}+h\) will have been issued at time steps \(t_{0}\), \(t_{0}-24\) and \(t_{0}-48\), respectively.

As the maximal horizon of the collected NWP forecasts is 72h, and we previously set the intra-day prediction interval as 24h, as a rule of thumb we used 48h as the context interval so that a training sample covers 72h. In practice, preliminary tests showed that using a larger context interval would not bring visible improvements, and using a multiple of the prediction interval size facilitates the creation of train and test data sets.
To measure model performance, RMSE-based metrics are common in energy utility companies, notably as they penalize large errors [16]. We used the normalized Root-Mean-Square Error (nRMSE) defined as:

\[\text{nRMSE}(\hat{\mathbf{Z}},\mathbf{Z})=\sqrt{\sum_{n=1}^{N}\frac{1}{N}\frac{ \frac{1}{T-t_{0}}\sum_{t_{0}+1}^{T}(\hat{z}_{nt}-z_{nt})^{2}}{P_{n}^{2}}} \tag{6}\]

with \(P_{n}\) the nominal power of PV system \(n\), \(\hat{\mathbf{Z}}=\{\hat{z}_{nt}\}\) the estimated point forecast, and \(\mathbf{Z}=\{z_{nt}\}\) the observed power. This nRMSE allows to measure the performance of a point estimate forecast in such way that PV systems with larger nominal power do not dominate the error metric. This is a field requirement, as PV systems have private owners, who have to be treated equally, irrespective of the nominal power of their system. In practice, nRMSE can be interpreted as a percentage of the PV system nominal power. To evaluate the performance of a proposed system w.r.t. a reference, the _skill score_ is derived from RMSE metrics as:

\[\text{Skill score}=1-\frac{\text{nRMSE}(\hat{\mathbf{Z}},\mathbf{Z})}{\text {nRMSE}(\hat{\mathbf{Z}}_{\text{ref}},\mathbf{Z})} \tag{7}\]

As presented in Section 3, the models trained in the context of this work output prediction quantiles. We use the median as the point estimate forecast; in addition, we compute the CRPS metric, commonly used in related work [17, 18, 19], which rates the quality of prediction quantiles as a whole:

\[\text{CRPS}(F^{-1},\mathbf{Z}) =\frac{1}{N(T-t_{0})}\sum_{n=1}^{N}\sum_{t=t_{0}+1}^{T}\int_{0}^{ 1}2\Lambda_{\alpha}(F^{-1}(\alpha),z_{nt})d\alpha \tag{8}\] \[\Lambda_{\alpha}(F^{-1}(\alpha),z_{nt}) =(\alpha-\mathcal{I}_{[z_{nt}<F^{-1}(\alpha)]})(z_{nt}-F^{-1}( \alpha))\]

with \(F^{-1}\) the quantile function of the predictor (which returns the quantile level in Watts associated to a probability \(\alpha\in]0,1[\)), \(\Lambda_{\alpha}\) the _quantile loss_, and \(\mathcal{I}_{[c]}\) the indicator function associated to logical clause \(c\). As discussed in Section 3, the quantile function is estimated empirically using a set of sample paths \(\{\hat{\mathbf{z}}_{n}\}\). Intuitively, the quantile loss gets larger when observations are far from the median as measured by distribution quantiles. This also allows to penalize models which are excessively confident, and somehow reward models which are able to better estimate the expected accuracy of their point forecast. In our experiments, we use 100 paths per sample, which can be used as input to return empirical quantiles.

PV power is naturally zero at night time. Therefore, including these time steps in metric computation is likely to bias nRMSE and CRPS towards 0. To prevent this, we exclude nightly time steps by limiting terms accounted for in Equation (6) and (8) to time steps where the value averaged over the whole data set is significantly different from zero.

### Validation scheme

Following recommendations by [17], we create training samples by cutting the data set in fixed size 72h segments, with \(t_{0}\) in each segment being midnight
24h before the end of the segment. Assuming we extracted the segment at the beginning of the time series, we then shift the offset 24h forward, so that the next segment includes the previous prediction interval in its context interval (see, e.g., top of Figure 1). As we treat all PV systems as independent time series, this results in \(O(CD)\) series, with \(C\) the number of PV systems and \(D\) the number of values \(t_{0}\) can take in the original time series.

Our number of PV systems and temporal collection bounds would yield 86989 samples. However, PV systems may exhibit missing or erroneous measurements due to several reasons (e.g., power outage, bad manipulation, faulty sensor). Figure 3 summarizes how missing values are distributed in the data set. The l.h.s. of Figure 3 shows that missing values are not uniformly distributed across PV systems. Approximately one third has no missing value, another third has a bit less than 20% of missing values, and the last third between 25% and 50%. The under-representation of this last third can be problematic. The r.h.s. of Figure 3 shows that these missing values are not evenly distributed in time: this indicates that a group of system may have been offline for a contiguous time frame during late winter and spring. Actually, most missing values are linked to systems started later than the others in year 2020. The associated periods are therefore under-represented, but we note that any month has at most 30% missing data. In the remainder, we consider that this bias remains in a range which makes uniform sampling w.r.t. time acceptable for building training batches. In order to facilitate processing, and as samples cuts are aligned with day frames, we detect and exclude days matching one of the following patterns: more than two consecutive missing values, measurements blocked to a constant value, visual inspection for aberrant values. This results in 67666 valid day frames.

The PV systems are distributed in a relatively small area in Luxembourg: therefore, it is expected that prediction intervals for different systems but same absolute time attached to \(t_{0}\) will be highly correlated. In order to validate this intuition, we computed all \(D\) intra-day correlation matrices between systems in our data set. Specifically, we defined intra-day time steps (excluding nightly time steps) as observations and PV systems as variables, resulting in \(D\frac{C(C-1)}{2}\) distinct correlation values. We observe that the median of the distribution of these correlation values is 0.95, which confirms a very high correlation between systems for a given day. As a consequence, sampling uniformly training, validation and test sets in the \(O(CD)\) series would result in _data leakage_, i.e., the model will be able to overfit without harming test error as identical samples (up to scale) will be scattered in training, validation and test sets. Let

Figure 3: _l.h.s._: Proportion of missing values per system. _r.h.s._: Proportion of missing data per associated month.
us note an unexpected positive benefit of this strong intra-day correlation: the under-representation of some systems is then much less problematic. The only remaining issue would pertain to estimating the parameters associated to the static categorical modalities of these systems, if using system ID static covariates. We hypothesise that at least 50% of represented day frames is sufficient to perform this estimation.

To prevent the data leakage problem, we first group the time series by the absolute time attached to their respective \(t_{0}\) (hence grouping samples for all PV systems for a given \(t_{0}\)), and sample 60% of the \(D\) time steps as the training set. We use a static temporal pattern, in order to ensure that each month and season is represented fairly evenly. Validation and test sets are uniformly sampled as half of the remaining 40% irrespective of the PV system, which is not an issue as the goal of validation error is to be an estimate of the test error, provided parameters are not explicitly fitted on the former. The validation set is used to implement early stopping, and select the model before it starts to overfit. The test set serves to compute the metrics described in Section 4.2. To choose the cut between validation and test, and ensure validation error is a fair proxy of the test error, we resample cuts until the RMSE and nRMSE between ground truth and intra-day NWP forecasts (which are considered as constant and known in advance, and are the most relevant baseline to compare to) are equal up to a small threshold. Using the cutting procedure defined so far, we obtain 40670 training, 13498 validation and 13498 test samples.

### Hyper-parameters

The models were trained using the Adam optimizer with learning rate \(10^{-3}\), batch size 64, for 200 epochs. Samples are reshuffled at the beginning of each epoch. In the end, we implement a form of early stopping, by selecting the model with best validation error. DeepAR uses 2 LSTM layers by default, we stick to this parametrization. Two free parameters remain then: the LSTM hidden layer size, and the number of components when a mixture distribution output is used. Figure 4 shows the results of a hyper-parameter search of these parameters using the Gaussian distribution as mixture components. Median results from 6 independently trained models are displayed in the graphs, along with \(\pm 1\) standard deviation error bars. First we determine the optimal LSTM hidden layer size using a single component mixture on the l.h.s. of Figure 4. We see that beyond 100 units, the validation performance is only marginally improved: we therefore retain 100 as hidden layer size for all our experiments. The r.h.s. of Figure 4 shows the validation performance w.r.t. the number of mixture components using this layer size. Using 2 components yields the best results, so we retain this mixture size for the remainder of the experiments.

### Model comparison strategy

Besides these hyper-parameters, a large search space of models results from using a given distribution output (Gaussian, Student, Positive Gaussian), one or two mixture components, whether or not using NWP, system ID, or system description covariates. Instead of comparing all possible combinations, we opted for a stepwise approach. We first compared 2-component mixtures of Gaussian, Student, and positive Gaussian components, with NWP covariates,
but without system static covariates (models 3 to 5 in Table 1). We chose this configuration as a middle ground, as using the NWP features implements the hybrid-physical approach advocated in this paper. This first comparison highlighted that the Student distribution output yields significantly inferior results in terms of nRMSE. We thus exclude the Student distribution from further tested combinations. This middle ground can be compared to the alternative quantile regression architecture proposed by [23], with access to equivalent input information (model 13).

Then we performed an ablation study of our middle ground. First we tested the consequence of using a single component instead of mixture (models 1 and 2). We also evaluated performance without NWP covariates, i.e. implementing purely autoregressive models (models 10 and 11). The latter are compared to a FFN model (model 12), which forecasts \(\mathbf{z}_{t_{0}+1:T}\) as a function of \(\mathbf{z}_{1:t_{0}}\) without any form of autoregression or covariates support, so that the value of using LSTM cells is evaluated _per se_.

Finally, the improvement brought by using the system ID (models 6 and 8) and system description features (models 7 and 9) to the middle ground configuration is evaluated. System ID is also used for model 11, as a way to estimate the best performance that can be reached (i.e. using a mixture of the best distribution output _a posteriori_) without NWP features. This can be of practical interest, as access to ECMWF solar irradiance forecasts is free for research purpose, but requires a subscription for industrial applications.

### Results and interpretation

Results are given in Table 1. In [13], skill scores are computed using a 24h persistence model, adjusted according to the clear sky PV power for the day under consideration. This is a common baseline in the solar energy domain [10]. In this paper, we rather consider that the NWP covariates are the baseline against which our results have to be evaluated. So we use 24h NWP covariates as \(\hat{\mathbf{Z}}_{\text{ref}}\) in Equation (7) for computing skill scores presented in Table 1. This is a stronger baseline for skill score computation, as it was shown to significantly outperform persistence forecasts [15]. In this experimental section, this means that a model with a skill score lower than 0 is not able to beat the covariates it is given among its inputs.

First focusing on the nRMSE metric through skill scores, we see that middle ground models with Gaussian and positive Gaussian yield some improvement,

Figure 4: Validation (orange) and test (blue) nRMSE curves for variable LSTM layer sizes (_l.h.s._) and numbers of mixture components (_r.h.s._).
with skill scores 1.04% and 2.23%, respectively. On the other hand, using the Student distribution yields skill score -0.53%. Using the Student distribution is therefore not even able to do better than just copying part of its inputs. This motivated pruning this option beyond the middle ground. We also note that the quantile regression method of [17] (model 13) obtains a skill score of -7.17%, which justifies our choice of DeepAR as the framework for our proposal.

Using a single Gaussian does also not match the baseline. On the other hand, a single positive Gaussian component yields a skill score of 0.50%. Using a mixture of distributions therefore contributes to improving the performance. Then, combining the system ID covariate to the mixture of Gaussian and positive Gaussian components yields skill scores 2.16% and 7.54%, respectively. We see that an important performance gap results from the joint usage of the system ID covariate and the positive Gaussian component. Using the system description covariates, these scores drop to 1.89% and 5.12%. Using the system ID is therefore the best option, but the system description features have the advantage to generalize to new systems without having to fully retrain the model, which can be useful in production conditions. We note that the improvement brought by using the system ID as a covariate provides anecdotal evidence supporting the hypothesis formulated in Section 4.3 regarding the imbalance of systems representation in the dataset. Using the system description covariates is beneficial (e.g., causes skill score increase by 2.89 points with the mixture of positive Gaussian components), but to a lesser extent than using the system ID. We can relate this to the fact that they do not fully reflect some local effects which already taint the PV performance model (see Section 1).

If ignoring NWP covariates, the performance of models is very significantly degraded. The mixtures of Gaussian and positive Gaussian components get 26.3% and 23.0% negative skill scores, respectively. Both models are still doing better than the alternative FFN architecture (28.4%), but the fallback consisting of not using NWP features comes at a high cost in terms of performance. We note that the gap between model 10 and 11 (2.7%) is not as large as the gap between their counterparts using NWP features (model 3 and 8, 6.5%). This con

\begin{table}
\begin{tabular}{c l l l l l l} \hline \multicolumn{2}{c}{**ID**} & \multicolumn{1}{c}{**Output**} & \multicolumn{1}{c}{**NWP**} & \multicolumn{1}{c}{**Static**} & \multicolumn{1}{c}{**nRMSE (\%)**} & \multicolumn{1}{c}{**Skill (\%)**} & \multicolumn{1}{c}{**CRPS (-)**} \\ \hline \multicolumn{2}{c}{} & \multicolumn{5}{c}{_Baseline_} & \multicolumn{1}{c}{} \\ \hline \multicolumn{2}{c}{-} & PV perf. model & Yes & - & 9.651 & - & - \\ \hline \multicolumn{2}{c}{_Single-component DeepAR_} & \multicolumn{1}{c}{} \\ \hline
1 & Gaussian & Yes & - & 9.697(\(\pm\)0.035) & -0.48 & 0.639(\(\pm\)0.004) \\
2 & Positive & Yes & - & 9.603(\(\pm\)0.041) & 0.50 & 0.650(\(\pm\)0.015) \\ \hline \multicolumn{2}{c}{_Mixture DeepAR_} & \multicolumn{1}{c}{} \\ \hline
3 & Gaussian & Yes & - & 9.551(\(\pm\)0.089) & 1.04 & 0.626(\(\pm\)0.002) \\
4 & Student & Yes & - & 9.702(\(\pm\)0.023) & -0.53 & 0.629(\(\pm\)0.004) \\
5 & Positive & Yes & - & 9.436(\(\pm\)0.049) & 2.23 & 0.618(\(\pm\)0.005) \\
6 & Gaussian & Yes & System ID & 9.443(\(\pm\)0.090) & 2.16 & 0.605(\(\pm\)0.004) \\
7 & Gaussian & Yes & System descr. & 9.469(\(\pm\)0.070) & 1.89 & 0.620(\(\pm\)0.004) \\
8 & Positive & Yes & System ID & **8.923**(\(\pm\)**0.044) & **7.54** & **0.579(\(\pm\)0.007)** \\
9 & Positive & Yes & System descr. & 9.157(\(\pm\)0.041) & 5.12 & 0.596(\(\pm\)0.004) \\
10 & Gaussian & No & - & 12.185(\(\pm\)0.055) & -26.3 & 0.831(\(\pm\)0.008) \\
11 & Positive & No & System ID & 11.873(\(\pm\)0.036) & -23.0 & 0.811(\(\pm\)0.006) \\ \hline \multicolumn{2}{c}{_Alternative models_} & \multicolumn{1}{c}{} \\ \hline
12 & FFN & No & - & 12.392(\(\pm\)0.089) & -28.4 & 0.879(\(\pm\)0.017) \\
13 & Wen et al. & Yes & - & 10.343(\(\pm\)0.070) & -7.17 & 0.686(\(\pm\)0.009) \\ \hline \end{tabular}
\end{table}
Table 1: nRMSE and CRPS test metrics for the range of compared models. The results of the best performing model are bold-faced. Median results and standard deviations are estimated from 6 models trained independently.
firms a synergy between the elements composing of the best-performing model (i.e. mixture of positive Gaussian components, NWP covariates, and system ID).

Figure 5 illustrates the relationship between nRMSE and CRPS metrics. We see that they are strongly correlated, so in the context of our experiments, models with the best skill scores generally yield the best prediction intervals. Models significantly below the linear regression fit will tend to provide better prediction intervals. This is the case with model 4 (mixture of Student components). However, its CRPS is still not as good as that of the other middle ground models (3 and 5). On the other hand, model 2 (single positive Gaussian component) has significantly degraded CRPS. All other models are quite close to the regression line. This includes models based on the mixture of positive Gaussian components: using a mixture seems to fix the discrepancy observed with model 2.

As alternative designs, we considered using unscaled NWP features (i.e. not scaling them along with \(\mathbf{z}\) values as described in Section 4.1), and using the weighted sampling scheme described in [13], which samples training examples with frequency proportional to the magnitude of \(\mathbf{z}\) values. We did not report results with these alternative designs as they brought systematic degradation to the performance. We note that the poor performance of weighted sampling w.r.t. the nRMSE metric is expected, as the latter balances the leverage of systems with large nominal power. We also tried to use both static covariates (i.e., system ID and description) simultaneously, but this led to a slight degradation compared to the respective model using only the system ID covariate. This is expected, as the system ID alone already encodes system di

Figure 5: Illustration of the relationship between nRMSE and CRPS metrics. The line is the result of a linear regression of the points in the graph. Glyph shapes and colors recall characteristics of the respective models. For better legibility, outlying models (i.e. those not using ECWF covariates) are excluded, even though they were also used for fitting the linear regression.
versity. In addition, as we already saw above, system description features may reflect the system setup in a biased way - even though using them is better than using no covariates at all.

### Discussion and qualitative examples

In the previous section we evaluated our proposed models using global metrics. In this section, we aim at providing more detailed insight into our results by analyzing model performance at the system level. The displayed examples were obtained using the best performing model identified in the previous section (i.e. model 8). First we compute per-sample nRMSE metrics for the test set, group them according to their associated system ID, and rank the systems according to the difference between model 8 and baseline nRMSE. In other words, the higher a system is in this ranking, the more DeepAR outperforms the PV performance model for this system. For all but 7 systems over 118, model 8 performs better than the baseline. As a worst case example, we first consider system 115 which comes last in this ranking.

On the l.h.s. of Figure 6, we display the sample of this system associated to

Figure 6: Two test samples associated to system 115. Confidence bounds are displayed for the prediction interval using green shades. The observations and NWP forecasts of the context interval are prepended.

Figure 7: Two test samples associated to system 44. Confidence bounds are displayed for the prediction interval using green shades. The observations and NWP forecasts of the context interval are prepended.
the lowest nRMSE with the DeepAR model. This is a typical clear day example, where the prediction is fairly easy for the neural model. We note that for this instance, forecasts stick more closely to the observations curve than the baseline. The confidence intervals are naturally tight, reflecting the high confidence of the neural model in this case. On the r.h.s. of Figure 6, for the same system, we display a sample for which the difference between the two models is among the largest. In this case, DeepAR is not able to keep up with the sudden peak of PV power. The 24h NWP covariates somehow informed on this peak, but this information was not used by model 8, which acted conservatively regarding the observations in the context interval.

In Figure 7, we consider samples from system 44. This system is the second best of our ranking, and has been identified as problematic for the PV performance model because of a double-pitched roof, not reflected by the system description features [13]. On both sides of Figure 7, the systematic shift of the PV performance model is clearly visible. We also see that model 8 is able to completely ignore this shift, and come up with sensible forecasts. The figure also shows how the confidence interval is tighter when the PV production curve was straightforward to forecast (l.h.s.), and broader when the day ahead is harder to forecast (r.h.s.).

## 5 Conclusion

Eventually, we are able to improve power forecasts obtained from an already strong PV performance model. By comparing many model variants, our experiments highlight the best working configuration, which uses the PV performance model forecasts as covariates, a mixture of positive Gaussians as output distribution, and a static categorical covariate reflecting the associated system ID. The positive Gaussian output allows to deal effectively with the _bell-shaped_ data profile typical of solar energy applications, and the system ID feature allows to model local effects which went previously unnoticed with the PV performance model alone.

In future work, we plan to refine and explore novel neural model designs. For example, quantile regression methods more recent than [14] will be explored. Also, we will further investigate how to deal with novel systems being added to the grid without having to retrain the full model. We saw that using system description features is an effective fallback, but these features do not account for local effects such as a double-pitched roof, so they remain suboptimal. We will also consider longer prediction intervals (e.g. day-ahead and 2 days ahead).

Despite visible success, models trained for this work tended to overfit, and relied critically on early stopping. This is mostly due to our measures taken to prevent data leakage: when segmenting 2 years of data at the day scale, despite all our measures, training and test sets are unlikely to be identically distributed. We addressed this problem in the most straightforward and conservative way, but it seems related to the domain shift problem characterized by the domain adaptation literature [15]. Adapting contributions from this area to the peculiarities of our application is left for future work.
Acknowledgements

This work was supported by the Luxembourg National Research Fund (FNR) in the framework of the FNR BRIDGES Project _CombiCast_

(BRIDGES18/IS/12705349/Combi-Cast). Furthermore, the authors would like to thank our partner Electirs (a brand of Hoffmann Freres Energie et Bois s.a r.l.), for their trust, the very supportive partnership throughout the whole project duration, and their contribution to the common project, financially as well as in terms of manpower and data.
