Xwégnon Ghislain

Robin Girard

Senior Member, IEEE George Kariniotakis

Probabilistic Model for Spatio-Temporal Photovoltaic Power Forecasting

Keywords: Lasso, photovoltaic generation, probabilistic forecasts, quantile regression, reliability, sharpness, spatiotemporal

Photovoltaic (PV) power generation is characterized by significant variability. Accurate PV forecasts are a prerequisite to securely and economically operating electricity networks, especially in the case of large-scale penetration. In this paper, we propose a probabilistic spatio-temporal model for the PV power production that exploits production information from neighboring plants. The model provides the complete future probability density function of PV production for very short-term horizons (0-6 hours). The method is based on quantile regression and a L1 penalization technique for automatic selection of the input variables. The proposed modeling chain is simple, making the model fast and scalable to direct on-line application. The performance of the proposed approach is evaluated using a real-world test case, with a high number of geographically distributed PV installations and by comparison with state-ofthe-art probabilistic methods.

I. INTRODUCTION

R ENEWABLE Energy Sources (RES), including photo- voltaic (PV) production, are being developed in many countries as a response to the need for clean energy solutions. PV production is characterized by high variability and uncertainty due to its inherent dependence on changing meteorological conditions. Variability and uncertainty are a challenge for network operators especially in systems with large-scale RES integration. Short-term forecasting of PV production is a prerequisite for the economic and secure management of power systems, reduction of reserve costs [START_REF] Potter | Building a smarter smart grid through better renewable energy information[END_REF] and market participation of PV producers. It is also important to ensure the competitiveness of renewable energy technologies [START_REF] Pinson | Trading Wind Generation From Short-Term Probabilistic Forecasts of Wind Power[END_REF].

The literature features several methods to forecast PV production. Detailed reviews of the state-of-the-art are provided in [START_REF] Kariniotakis | Renewable Energy Forecasting: From Models to Applications[END_REF]- [START_REF] Sobri | Solar photovoltaic generation forecasting methods: A review[END_REF]. They can be classified according to the forecast horizon, the available data, and the type of approach, which may be based on statistics, physics or a hybrid combination [START_REF] Antonanzas | Review of photovoltaic power forecasting[END_REF], [START_REF] Zamo | A benchmark of statistical regression methods for short-term forecasting of photovoltaic electricity production, part I: Deterministic forecast of hourly production[END_REF]. Although early methods were deterministic, probabilistic approaches are increasingly popular since they provide additional information about the distribution of future production This work was carried out within the research project entitled "Improvement of PV power forecasting and predictive management including storage solutions", supported by the company Coruscant SA in the frame of its participation to a tender of the French Energy Regulator CRE for the development of PV plants above 250 kWc.

The authors are with MINES ParisTech, PSL -Research University, PERSEE -Centre for Processes, Renewable Energies and Energy Systems CS 10207 rue Claude Daunesse, 06904 Sophia Antipolis Cedex, France. (emails: xwegnon.agoua, robin.girard, georges.kariniotakis each with @minesparistech.fr) and uncertainty in the forecasts. Some of these probabilistic approaches are based on Numerical Weather Prediction (NWP) models or sky imaging and provide ensemble forecasts of the future PV generation [START_REF] Sperati | An application of the ecmwf ensemble prediction system for short-term solar power forecasting[END_REF]- [START_REF] Liu | Ensemble forecasting of solar irradiance by applying a mesoscale meteorological model[END_REF]. Analog ensembles [START_REF] Alessandrini | An analog ensemble for short-term probabilistic solar power forecast[END_REF], regression trees [START_REF] Zamo | A benchmark of statistical regression methods for short-term forecasting of photovoltaic electricity production. Part II: Probabilistic forecast of daily production[END_REF], [START_REF] Nagy | GEFCom2014: Probabilistic solar and wind power forecasting using a generalized additive tree ensemble approach[END_REF] and k-nearest neighbors (kNN) [START_REF] Huang | A semi-empirical approach using gradient boosting and -nearest neighbors regression for GEFCom2014 probabilistic solar power forecasting[END_REF]- [START_REF] Chu | Short-term probabilistic forecasts for Direct Normal Irradiance[END_REF] are also found in the related literature on probabilistic PV forecasting. A wide range of Artificial Neural Networks (ANN) based models also exist for short-term PV power production [START_REF] Yona | Application of Neural Network to One-Day-Ahead 24 hours Generating Power Forecasting for Photovoltaic System[END_REF]- [START_REF] Vagropoulos | ANN-based scenario generation methodology for stochastic variables of electric power systems[END_REF]. These models have evolved from simple neural networks, to radial neural networks (more suitable for time series prediction) and more recently to deep learning methods [START_REF] Zhang | Predictive Deep Boltzmann Machine for Multiperiod Wind Speed Forecasting[END_REF], [START_REF] Khodayar | Rough Deep Neural Architecture for Short-Term Wind Speed Forecasting[END_REF].

Recently the focus has shifted to improving short-term predictability of solar irradiance or PV production by considering off-site information as input to the models. Such models are commonly referred to as spatio-temporal models. Models for solar radiation forecasting are based on autoregressive models [START_REF] Glasbey | A spatiotemporal auto-regressive moving average model for solar radiation[END_REF]- [START_REF] Amaro E Silva | Impact of network layout and time resolution on spatio-temporal solar forecasting[END_REF], geostatistical models [START_REF] Journée | Improving the spatio-temporal distribution of surface solar radiation data by merging ground and satellite measurements[END_REF]- [START_REF] Quesada-Ruiz | Cloud-tracking methodology for intra-hour DNI forecasting[END_REF] or classification models [START_REF] Lazzaroni | Models for solar radiation prediction based on different measurement sites[END_REF]. Physical models combined with NWP forecasts [START_REF] Lorenz | Local and regional photovoltaic power prediction for large scale grid integration: Assessment of a new algorithm for snow detection[END_REF], [START_REF] Lorenz | Comparison of global horizontal irradiance forecasts based on numerical weather prediction models with different spatio-temporal resolutions[END_REF] and semi-parametric models [START_REF] Patrick | A semiparametric spatio-temporal model for solar irradiance data[END_REF] can also be used for spatio-temporal forecasting of solar radiation. ANN-based models are also used for spatio-temporal solar power forecasting [START_REF] Gutierrez-Corea | Forecasting short-term solar irradiance based on artificial neural networks and data from neighboring meteorological stations[END_REF]- [START_REF] Li | Machine learning for solar irradiance forecasting of photovoltaic system[END_REF]. For the case of PV power forecasting, most spatio-temporal models use off-site meteorological measurements and neighboring site information to improve the forecast for the site of interest. These models are autoregressive models [START_REF] Tascikaraoglu | Compressive Spatio-Temporal Forecasting of Meteorological Quantities and Photovoltaic Power[END_REF]- [START_REF] Yang | A novel ARX-based multi-scale spatio-temporal solar power forecast model[END_REF] or classification models [START_REF] Lazzaroni | Models for solar radiation prediction based on different measurement sites[END_REF]. The classification model proposed in [START_REF] Berdugo | Analog method for collaborative very-short-term forecasting of power generation from photovoltaic systems[END_REF] handles data acquisition and privacy issues.

The spatial-temporal models listed above provide deterministic forecasts. In contrast to day-ahead forecasting, very few spatial-temporal models for short-term probabilistic forecasting feature in the literature. These models are based on regression trees [START_REF] Nagy | GEFCom2014: Probabilistic solar and wind power forecasting using a generalized additive tree ensemble approach[END_REF], the kNN method [START_REF] Huang | A semi-empirical approach using gradient boosting and -nearest neighbors regression for GEFCom2014 probabilistic solar power forecasting[END_REF], the combination of a vectorial autoregressive model and gradient boosting [START_REF] Bessa | Probabilistic solar power forecasting in smart grids using distributed information[END_REF], multivariate predictive distributions, [START_REF] Golestaneh | Generation and evaluation of space-time trajectories of photovoltaic power[END_REF] and Gaussian random fields [START_REF] Zhang | Spatial-temporal solar power forecast through use of Gaussian Conditional Random Fields[END_REF].

One of the biggest challenges in spatio-temporal forecasting is to propose models that are able to handle the dimensionality issues raised by the large amount of explanatory variables. This is a prerequisite for the models so that they can generalize efficiently when they "see" new data. In this paper, we propose a short-term probabilistic forecasting methodology that exploits information available from a large number of PV installations and includes such a process for the automated selection of explanatory variables. The aim is to improve predictability in the time range of 0-6 hours.

In a previous paper [START_REF] Agoua | Short-term spatiotemporal forecasting of photovoltaic power production[END_REF], we proposed a deterministic spatio-temporal model for PV power forecasting. This paper makes several additional contributions, compared to both the previous paper and to the related literature on probabilistic PV forecasting. These contributions are: 1) the use of offsite data (neighboring production data) to generate probabilistic forecasts with improved properties and performance compared to the standard approach without off-site data 2) the introduction of an automated variable selection procedure for the PV power forecasts inside the model, which avoids the dimensionality problem 3) the integration of NWP forecasts as explanatory variables despite the short-term frame; this allows us to assess the effect of the NWP variables on the probabilistic forecasts. In probabilistic models in the literature [START_REF] He | Short-term power load probability density forecasting based on quantile regression neural network and triangle kernel function[END_REF], [START_REF] Bessa | Probabilistic solar power forecasting in smart grids using distributed information[END_REF], [START_REF] Tastu | Probabilistic Forecasts of Wind Power Generation Accounting for Geographically Dispersed Information[END_REF] the variable selection is done either empirically or with the variable importance criteria for decision trees; 4) two probabilistic models chosen from the most efficient ones in the literature are used as references for our model evaluation.

In contrast to methods proposed in the literature for this time range, our approach is not based on information from sky cameras or satellite images. It uses available data from off-site PV plants. The model uses the transfer of spatial and temporal information within the network of power plants to capture the meteorological perturbations. Unlike in [START_REF] Tascikaraoglu | A short-term spatiotemporal approach for Photovoltaic power forecasting[END_REF], these spatial and temporal correlations between power plants are useful for short-term horizons as they offer good spatial and temporal resolution.

Moreover, in the proposed model we use NWP forecasts as additional explanatory input, which is not standard practice in the literature for models developed especially for this time scale. NWPs are a standard input for models that address dayahead or longer horizons. Such models also cover horizons shorter than 6h. It is however widely recognized in renewable generation forecasting that, for the first hours, NWPs do not obtain performance improvements w.r.t. persistence. To achieve such improvements, it is necessary to include recent measurements, such as input. This observation, coupled with the need for frequent forecast updates (i.e. every 30 minutes) for applications like intraday trading, have motivated the development of purely statistical models based only on measurements. The standard spatiotemporal approach comes into this category. However, we have observed that including NWPs in this family of models can provide information on weather condition trends for horizons around 6 hours ahead. In cases like wind prediction, this inclusion brings improvements w.r.t. persistence, amounting to double the improvements without NWPs for 6 hours ahead [START_REF] Kariniotakis | Renewable Energy Forecasting: From Models to Applications[END_REF]. This inclusion of NWPs in the family of short-term models that operate with very frequent updates goes against the general tendency in the literature. However, this approach allows for better coupling between very short-term spatial-temporal forecasts and standard dayahead forecasts since NWPs reduce? the mismatch between two types of forecast that can be confusing for end-users. The most relevant NWP forecasts were selected by measuring their dependence on production data.

The generation of the probabilistic forecast was done using a quantile regression scheme. Therefore the variable selection process using L1-penalization that we propose in this paper is different from the one used in [START_REF] Agoua | Short-term spatiotemporal forecasting of photovoltaic power production[END_REF] as it is done in a L1-optimization model (quantile regression). We have also added the constraints of different parameters for each quantile and each time step.

This paper puts forward the non-stationary characteristic of PV power series. The spatio-temporal approach exploits correlations in the data from geographically dispersed sites. It is thus necessary to efficiently eliminate the deterministic effect that is related to the course of the sun. This issue has been treated by applying an advanced stationarization process to the production series as proposed in [START_REF] Agoua | Short-term spatiotemporal forecasting of photovoltaic power production[END_REF]. The evaluation is carried out on a real world test case of 185 PV installations in France. This test case is a good illustration of a system with high penetration of renewable energy and also allows us to demonstrate the full potential of a spatio-temporal approach for probabilistic forecasting when considering a highly dimensional problem.

The paper is structured as follows: the data and an analysis of the spatio-temporal correlations are presented in section II. The proposed spatio-temporal probabilistic models are presented in section III. The evaluation and analysis of the performances of the forecasts are treated in section IV. Finally, the conclusions of the study are discussed in section V.

II. THE INTEREST OF SPATIO-TEMPORAL MODELING

The aim of this section is to illustrate the interest and potential of spatio-temporal modeling in short-term PV power forecasting. A real-world case study is introduced followed by an analysis illustrating that useful information can be extracted from geographically distributed PV installations.

A. Data

The data set considered here, denoted as d, comprises the output of 905 PV power inverters in a mid-west region of France with peak power ranging from 3.2 kWc to 58 kWc. This number of inverters corresponds to 185 different PV power plants (set of inverters at the same location). The distance between them varies from 1 km to 230 km. The data relate to the period from November 2014 to March 2016 with a 15 min temporal resolution. The locations of the power plants are represented in Figure 1. After data cleaning, 136 power plants were retained. In the following, the power plants are labeled P i,1≤i≤136 . The production series have been stationarized employing the same procedure as that proposed in [START_REF] Agoua | Short-term spatiotemporal forecasting of photovoltaic power production[END_REF], which is used thereafter.

B. Spatial-Temporal Correlation

The existence of a spatial-temporal pattern is evaluated by an analysis of the correlation between the stationarized production series. In these series, the effects of the sun position on cross-correlations are absent.

Figure 2 presents the histograms of the lagged crosscorrelation values for three classes of distance. For a couple of stationarized series s 1 , s 2 for example, the value of the correlation is obtained by evaluating the correlation between s 1 and the lagged series lag(s 2 , k). The maximum lag used is the maximum horizon considered for the forecasts, here 6 hours. The same evaluation is made by switching s 1 and s 2 . The value retained is the maximum value. The correlation values are quite high, concentrated in the interval [0.4 -0.8] and, decrease with distance. Therefore, we can assume that these significant correlation values describe a spatial transfer of information between the power plants. This transfer, mainly due to cloud movements, can be used to anticipate meteorological perturbations and to improve the forecasts using production data only.

From a temporal point of view, Figure 3 presents the time lags for which the cross correlations described above reach their maximum. The idea behind this analysis is to determine the temporal limit of the propagation of correlation between the power plants. The figure shows that in the majority of cases, the maximum correlation is obtained for 15 minutes, but for some power plant couples it can be reached for 4 hours or more. The models proposed in the next section aim to exploit this spatial and temporal information.

III. PROBABILISTIC MODELS FOR SPATIAL-TEMPORAL PV FORECASTING

We consider here two models chosen from the most efficient probabilistic forecasting methods in the literature as reference Longitude difference (km) Time lag of maximum corr (x 15 min) Fig. 3. Time lags for which the correlation between each power plant couple is at a maximum. One point per couple. models: Kernel Density Estimation (KDE) and the Extreme Learning Machine (ELM) method proposed in [START_REF] Golestaneh | Very Short-Term Nonparametric Probabilistic Forecasting of Renewable Energy Generation #x2014; With Application to Solar Energy[END_REF]. The KDE model is a well-known method used to evaluate the future density of PV power [START_REF] Zhang | K-nearest neighbors and a kernel density estimator for GEFCom2014 probabilistic wind power forecasting[END_REF], [START_REF] Bacher | Online short-term solar power forecasting[END_REF]. An information-based criterion is proposed to select the input variables. These two models respectively use only "local" on-site information from the PV installation the power of which is predicted. These models can thus be used to evaluate the contribution of the proposed spatial-temporal probabilistic approach. The advanced method we propose is based on the quantile regression approach, which we adapted for high dimensionality problems. We describe below the principle of these methods, the conditions of their implementation and the modifications proposed to improve their efficiency.

A. Kernel Density Estimation

KDE is a non-parametric approach for density estimation [50]- [START_REF] Wand | Kernel Smoothing[END_REF]. This method reduces the level of estimation errors compared to parametric approaches as there is no hypothesis on the underlying distribution. The minimization problem consists in providing an estimation of the density f of a random variable X. The n-dimensional multivariate kernelestimator expression is:

f (x) = 1 N |H| N i=1 K H -1 (x -x i ) (1)
where x is the point where the evaluation is made, x i , i = 1 . . . N the data. H is a n × n matrix called the bandwidth matrix, which controls the smoothing and |H| its determinant. K is the kernel function. The kernel function K and the bandwidth matrix H are the two parameters which have to be determined to apply the KDE. The impact of the kernel function on the estimation quality is low and a standard Gaussian kernel implies a high computation cost [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF]. In this paper, we have chosen to use the Epanechnikov kernel function:

K(u) = 3 4 √ 5 (1 -u 2 /5) u ∈ [-1, 1]. (2)
The multivariate version of this kernel function is the product

K(u) = n j=1 K(u j ).
The smoothing matrix H is the most important parameter as it has a strong influence on the quality of the estimation [START_REF] Wand | Kernel Smoothing[END_REF]. Several approaches are possible for selecting the appropriate bandwidth matrix [START_REF] Mugdadi | A bandwidth selection for kernel density estimation of functions of random variables[END_REF], [START_REF] O'brien | A fast and objective multidimensional kernel density estimation method: fastKDE[END_REF]. In this paper, we use the smoothed cross validation technique [START_REF] Duong | Cross-validation Bandwidth Matrices for Multivariate Kernel Density Estimation[END_REF]. Moreover, as the PV production data are positive and bounded, we adapt the general method to our needs by using boundary correction as proposed in the literature [START_REF] Wand | Kernel Smoothing[END_REF], [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF].

Let Y ∈ R p be the random variable, the realizations of which are the production of a power plant and X ∈ R p the explanatory variables. The objective is to compute the probability density function of the conditioned random variable Y t+k |X t where t is the time the prediction is made and k the horizon. This density function is obtained by:

f Y t+k |Xt = f Y t+k ,Xt f Xt (3)
The forecast probability density function is :

fY t+k |Xt = 1 |H| N i=1 w(x, x i )K(H -1 (y -y i )) (4)
where

w(x, xi) = K H -1 (x -x i ) N j=1 K (H -1 (x -x j )) . (5)
The mutual information criterion [START_REF] Cover | Elements of information theory[END_REF] has been used for choosing the inputs of the KDE model. It is a measure of "distance" (Kullbac-Leiber) between two probability density functions. Here this criterion is used to evaluate the distance between the probability density functions of the PV production and those of the meteorological variables. This evaluation is carried out for each meteorological variable. The values of the information criterion are used to classify the meteorological variables according to their proximity to the probability density of the production, following which the most relevant inputs are selected for the KDE model. Let X be the random variable representing a meteorological variable and Y that of the PV production. If f X and f Y denote their respective probability density functions, the average mutual information between Y and X is:

I(X, Y ) = Y X f X,Y log f X,Y f X .f Y dxdy. (6)
A null value of this criterion indicates that the variables X and Y are independent.

B. ELM-Based Estimation

The Extreme Learning Machine method is a non-parametric density forecasting method presented and developed in [START_REF] Golestaneh | Very Short-Term Nonparametric Probabilistic Forecasting of Renewable Energy Generation #x2014; With Application to Solar Energy[END_REF]. It is based on feed forward neural networks but does not require tuning for the hidden layer parameters (the parameters are randomly selected making it fast). The proposed configuration of the ELM-based estimation model for this paper is the following:

• one model is tuned for each power plant;

• the predictors are the past solar observations of the site of interest with respective time lags and the NWP variables selected by the mutual information criterion presented in subsection III-A; • the results are generated as predicted quantiles for a 6-hour horizon with 15-min time-steps; the forecasting model is updated monthly • the non-crossing quantile constraint has been applied.

C. The proposed model: the QR-Lasso

The aim of quantile regression [START_REF] Koenker | Quantile Regression[END_REF] is to provide an estimation of the following cumulative distribution function:

F (y|X = x) = P(Y ≤ y|X = x) = E(1 Y ≤y |X = x). (7)
This estimation is made by providing information about a significant number of points describing the distribution known as the quantiles. For a continuous probability distribution function, the quantile of level α ∈ (0, 1), qα (Y |X) is

q α (Y |X) = F -1 (α) = inf{y, F(y|X = x) ≥ α}. (8)
This quantile is a solution of the minimization problem:

q α (Y ) = arg min g E[ρ α (Y -g(X))] (9)
where ρ α (u) = u (α -1 {u<0} ) is the loss function called the pinball loss. The problem ( 9) can be generalized to the conditional quantile of level α as :

q α (Y |X) = arg min g E[ρ α (Y -g(X))|X = x]. (10)
When assuming a linear relation between Y and X, the estimation can be made on the observation set using:

qα = arg min β n i=1 ρ α (Y i -X i β). (11)
The equation [START_REF] Zamo | A benchmark of statistical regression methods for short-term forecasting of photovoltaic electricity production. Part II: Probabilistic forecast of daily production[END_REF] does not have an explicit solution; the solution has been found numerically, despite the fact that the pinball loss is neither differentiable in 0 nor strictly convex. The problem is then transformed into a linear optimization problem solved by the simplex algorithm (for small samples) or interior point method (for large samples) [START_REF] Dantzig | Linear Programming and Extensions[END_REF]. The condition of non-crossing quantiles has been taken into account. In practice, this condition is implemented by adding some constraints to the problem that specify that each quantile value is positive, lower than the next value, and lower than the maximum.

D. Adaptation for the High Dimensional Input

Spatio-temporal forecasting, when considering numerous PV plants, raises the question of choosing the appropriate input variables for the models involved. Possible inputs include the production series of the predicted site and the neighboring sites, and the lagged version of these series. The NWP meteorological data can also be integrated into the models. Given that the PV plants may cover a large area of tens of kilometers, several points of the NWP grid can be considered. The above information represents a high amount of potential inputs for the forecasting models and raises the issue of their dimensionality and over-fitting. We propose here a variable selection process, which we have adapted to the quantile regression.

The method is based on a direct L1-penalization called the Least Absolute Shrinkage and Selection Operator (LASSO) [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF]. This penalization is applied directly to the estimation process of the quantile regression model. The estimator is defined for the quantile of level α as: βlasso = argmin

β {E[ρ α (Y -X β)] + λ β 1 } . ( 12
)
This penalization also helps to avoid over-fitting in the model. For values of the penalization parameter λ that are high enough, some coefficients β are set to zero, thus allowing a selection to be made between the input variables.

The problem ( 12) can be reformulated as

argmin β {E[ρ α (Y -X β)]} st β 1 ≤ Rλ . ( 13
)
The equations ( 12) and ( 13) are equivalent. The latter is the most common form of optimization problem. The Lasso is also more suitable for operational models as its run-time is lower than those of other methods, like L2-penalization, gradient boosting or the stepwise process. The L1 loss ensures a variable selection that makes the model perform well with a reasonable number of variables and is stricter than the L2 loss. More details about the regularization paths can be found in the literature [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF]. This modification of the standard quantile regression proposed here is essential in the spatial-temporal paradigm due to the high number of potential inputs. The new quantile regression model including this LASSO variable selection procedure is hereafter called QR-Lasso. The performance of this model is evaluated against the reference KDE model.

IV. EVALUATION OF THE FORECASTS

The performance of the proposed model is evaluated for horizons up to 6 hours ahead with a 15-min time step. It is considered that forecasts are updated every 15 minutes. This performance is compared to the reference model. The NWP data were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF). These forecasts are updated 4 times a day and are given with a 3-hourly temporal resolution. 15-min values were obtained through interpolation. The models were developed using the software R [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] with the packages quantreg, glmnet, elmNN and some specific functions developed for this study (matrice estimations for KDE, penalized quantile, etc). The probabilistic (distribution) forecasts will be characterized by their quantiles. The last third of the data set (around 5 months) was used as a testing set. The QR-Lasso model is updated for each time step for each decile. A single update step took approximately 5 seconds (on a 12 GB RAM computer).

A. Variables Selection Through Mutual Information

The variable selection procedure allows us to choose the relevant variables from numerous available inputs, measurements and NWPs, for the case where a large number of PV plants are involved in the spatio-temporal modeling. In Table I, we present the normalized mutual information between the PV production of four PV plants and the NWP variables. The mutual information has not been computed for the top net solar radiation variable since this variable is the main driver of the PV production process. The table shows that the two most important variables are the temperature and the humidity followed by the wind direction. The importance of the precipitation level is the lowest.

Using this process, we select the following input NWP variables for the reference KDE model, which will then be multivariate:

• the top net solar radiation (TSR),

• the temperature (T),

• the relative humidity (RH),

• the wind direction (WD).

B. Analysis of the Lasso Variable Selection

The potential input of the proposed spatio-temporal QR-Lasso model for each PV plant includes the measured power of the other 135 plants, the lagged values of these measurements and finally the NWP variables corresponding to the 4 surrounding NWP grid points. The subset of NWP variables selected by the mutual information criterion is considered. Then for each PV plant P i , i = 1, . . . , 136, the number of total potential inputs is 819 = 4 (NWP forecasts) + (136 × 6) -1 (6 lags). Figure 4 presents the number of nonnull coefficients related to PV power inputs selected by the QR-Lasso model for each decile. The maximum number of variables selected is 320, showing that the selection process is efficient. The performance of the variable selection process is also evaluated by examining the geographical position of the selected PV plants with respect to the position of the PV plant of interest. Figure 5 presents the positions of the selected PV plants when forecasting the median of the production of two PV plants for a 6-hour horizon. The numbers of PV plants selected from the 136 initial plants are 33 (left) and 43 (right). For the two power plants of interest presented, the choice of neighboring plants is quite homogeneous. Some power plants close to the plant of interest do not provide more knowledge to the model (in terms of error reduction) than the selected ones. To ensure the parsimony, these plants are thus not selected by the model.

C. The Performance of the QR-Lasso Spatial-Temporal Approach

The quality of the probabilistic forecasts produced with the selected variables can be evaluated with several criteria [START_REF] Pinson | Non-parametric probabilistic forecasts of wind power: required properties and evaluation[END_REF]. Here we limit our study to the two most common criteria, i.e. reliability and sharpness, and a combined criterion that represents overall skill.

Reliability is a criterion that describes the ability of the probabilistic forecast model to match the expected observation frequencies. For qα a predicted quantile of the level α, the reliability criterion verifies whether the associated proportion of data observed under this quantile is equal to the expected α. To evaluate the reliability of the full predictive density for each of our models, we evaluate the reliability for quantile q 10 to q 90 with an increasing step of 10. The quantiles can be interpreted as estimated coverage rates, which then allows us to evaluate the reliability through increasing coverage rate analysis. Formally, reliability can be defined for the quantile level (or coverage rate) α and horizon h as:

Rel α h = α - 1 N N t=1 1{y t+h|t ≤ q α t+h|t } ( 14
)
where t is the time when the prediction is made. Figure 6 presents for these quantiles (predictive intervals) the observed deviation from the nominal coverage for 3h hori- zons. Perfect reliability means that there is no deviation, and this is represented by the horizontal black line on the figure.

The deviation of the ELM model has also been plotted. Two deviations are plotted for the KDE representing two different smoothing matrices H1, H2. These smoothing matrices are respectively the results of the unbiased cross validation and the smooth cross-validation method (see [START_REF] Duong | Cross-validation Bandwidth Matrices for Multivariate Kernel Density Estimation[END_REF] for more details on the functions and minimization process). The figure shows that, depending on the values of the smoothing matrices, the predictive densities can be either over-estimated (red line for KDE/H1) or under-estimated (green line for KDE/H2). The values of the deviation for the QR-Lasso model do not exceed ±0.5% and are lower than those of the KDE model. Moreover, the deviations of the nominal reliability values of the QR-Lasso model are lower than those of the ELM model (purple line). Sharpness is a complementary analysis tool to reliability and can be used to evaluate the concentration of the predictive distributions. Sharpness is often evaluated by taking the average width of centered prediction intervals. After defining a set of prediction intervals, the distances between the two boundaries (length of interval) are computed and all of these distance values are averaged over the evaluation set [START_REF] Gneiting | Strictly Proper Scoring Rules, Prediction, and Estimation[END_REF]. The more concentrated the predictive distributions are, the sharper the forecasts, which is preferable. The sharpness is defined for a quantile level α for the horizon h as:

sharp α h = 1 N N t=1 q1-α 2 t+h|t - q α 2 t+h|t . (15)
The evaluation of sharpness (as inter-quantile interval lengths) is presented in Figure 7 for the same prediction intervals as on the reliability plot. The maximum observed production is used to normalize the interval lengths. The figure shows that the average interval length increases with the nominal quantile (coverage rate). The values range from 3% to 58% and are lower than those observed for the KDE model for quantiles up to q 60 . The inter-quantile interval lengths in the QR-Lasso

Nominal coverage rate (%) Average interval length (%Pmax)

D. Overall Probabilistic Forecasting Skill

In addition to reliability and sharpness, an overall evaluation criterion for probabilistic forecasts is the Continuous Rank Probability Score (CRPS). The CRPS evaluates the entire predictive distribution and can be seen as a criterion that combines reliability with sharpness. It is defined for a cumulative distribution function F and its observation y as

crps(F, y) = ∞ -∞ (F (x) -H(x -y)) 2 dx (16)
where H(x) is the Heaviside function, whose value is 0 for strictly negative x and 1 for positive x. The CRPS can be interpreted as a measure of distance between the observed cumulative distribution function and the forecasted one [START_REF] Hersbach | Decomposition of the Continuous Ranked Probability Score for Ensemble Prediction Systems[END_REF].

The average CRPS values over the 6-hour forecast horizon are presented in Table II for the reference KDE model with matrix H1 (the one with the best performances), the ELM model and the QR-Lasso model for five power plants in the test case. The CRPS values are lowest for the QR-Lasso model, confirming the improved forecasting performance compared to the KDE and ELM models. This analysis provides one more illustration of how well the QR-Lasso model performs. Moreover, the CRPS values are lower than those of the set of models presented in the review [START_REF] Zamo | A benchmark of statistical regression methods for short-term forecasting of photovoltaic electricity production. Part II: Probabilistic forecast of daily production[END_REF].

Figure 8 represents the predictive densities as a set of prediction intervals for six days of the testing set for the power plant P 1 . The QR-Lasso model performs quite well to anticipate the variations in production. The length of the quantile intervals is low compared to what can be observed in the literature [START_REF] Golestaneh | Very Short-Term Nonparametric Probabilistic Forecasting of Renewable Energy Generation #x2014; With Application to Solar Energy[END_REF]. Overall, most of the variations for days characterized by high production variability are predicted quite efficiently.

V. CONCLUSION

In this paper we proposed a probabilistic spatio-temporal model that exploits off-site information to improve short-term forecasting of photovoltaic production. The model is based on quantile regression that has been adapted to integrate a LASSO variable selection process. This makes it particularly suitable for situations where data from a high number of PV installations is available since it is able to handle directly the resulting high dimensionality and over-fitting issues that characterize the use of such large amounts of data.

A real world test case with a high number of PV plants was used to illustrate the full potential of the probabilistic spatiotemporal approach. The evaluation of the predicted densities was carried out employing evaluation criteria that are widely used for probabilistic forecasts. The model shows significant performance improvement compared to the reference kernel density model and to models in the literature. The probabilistic spatio-temporal model appears to be an interesting technique that performs well to predict the future densities of PV generation.

The perspectives of this work include the possibility of extending the forecast horizons to daily forecasts. This would emphasize the trade-off significance between NWP forecasts and historic data measurements, and help to define the horizon limit where the spatial-temporal approach does not permit improvement. Another improvement could be to consider a time-varying smoothing matrix to estimate the kernel density, which is a key parameter of this model.

Fig. 1 .Fig. 2 .

 12 Fig. 1. The power plants of the test case d.

Fig. 4 .

 4 Fig. 4. Number of non-null coefficients computed by the QR-Lasso model by decile level (the meteorological coefficients are not included). Four PV plants (P 1 -P 4 ) from the test case are represented and each line corresponds to one power plant. The horizon is 6 hours.

Fig. 5 .

 5 Fig. 5. The power plants selected by the QR-Lasso model (in blue) when forecasting the median of the future production for a 6-hour horizon for two PV plants (in red).

Fig. 6 .

 6 Fig. 6. Reliability of the predictive densities for 3-hour horizons. The red and green lines represent the KDE model with two cases of smoothing matrix. The purple and blue lines represent the ELM and QR-Lasso models. The values are given as average over the forecasting length on the testing set.

Fig. 7 .

 7 Fig. 7. Sharpness evaluation of the predictive densities for 3-hour horizons. The red and green lines are the interval lengths of the KDE model for two cases of smoothing matrix. The purple and blue lines are those of the ELM and QR-Lasso models.

TABLE I NORMALIZED

 I MUTUAL INFORMATION VALUES FOR 4 PV PLANTS OF THE DATA SET COMPUTED BETWEEN PV PRODUCTION AND NWP VARIABLES.

	NWP Variables	P 1	Mutual information (%) P 2 P 3 P 4
	Temperature (T)	72.38	70.74	80.79	78.22
	Relative Humidity (RH)	70.56	74.56	74.83	71.27
	Wind Direction (WD)	25.11	21.88	26.11	21.86
	Wind Speed (WS)	6.18		5.78	6.57	5.46
	Precipitation	0.18		0.48	0.26	0.27

TABLE II CRPS

 II OF THE SPATIO-TEMPORAL QR-LASSO MODEL AND THE REFERENCE KDE MODEL FOR FIVE POWER PLANTS.

	Power plant	KDE	CRPS (% Pmax) ELM QR-Lasso
	P1	7.20	7.00	5.20
	P2	8.32	7.42	4.14
	P3	7.65	7.57	5.23
	P4	8.67	7.68	5.14
	P5	8.32	7.08	4.38
	model are higher than those in the KDE model for upper
	quantiles.

  Fig. 8. Examples of predictive densities forecasted with the QR-Lasso model for 6 days in different months of the year. The quantiles are plotted from 10% to 90% with a 10% step. The time step is 15 min and the model is updated every 15 min.

									2015-11-15				2015-12-05
			90%				70	90%					90%
			80% 70%				60	80% 70%				80	80% 70%
			60%					60%					60%
	Power [% pmax]		50% meas 40% 30% 20% 10%			Power [% pmax]	20 30 40 50	50% meas 40% 30% 20% 10%			Power [% pmax]	40 60	50% meas 40% 30% 20% 10%
												20
							10
							0					0
							0	20	40	60	80	0	20	40	60	80
				Time (x 15min)				Time (x 15min)				Time (x 15min)
				2016-01-15				2016-02-10				2016-02-19
		80	90% 80%				80	90% 80%				40	90% 80%
			70%					70%					70%
	Power [% pmax]	40 60	60% 50% meas 40% 30% 20% 10%			Power [% pmax]	40 60	60% 50% meas 40% 30% 20% 10%			Power [% pmax]	30 20	60% 50% meas 40% 30% 20% 10%
		20					20					10
		0					0					0
		0	20	40	60	80	0	20	40	60	80	0	20	40	60	80
				Time (x 15min)				Time (x 15min)				Time (x 15min)

ACKNOWLEDGMENT

The authors would like to thank the French industrial Hespul for providing the PV data as well as the European Center for Medium-Range Weather Forecasts (ECMWF) for providing the NWP data.