Prévision d'ensemble par agrégation séquentielle appliquée à la prévision de production d'énergie photovoltaïque.

Résumé : Notre principal objectif est d'améliorer la qualité des prévisions de production d'énergie photovoltaïque (PV). Ces prévisions sont imparfaites à cause des incertitudes météorologiques et de l'imprécision des modèles statistiques convertissant les prévisions météorologiques en prévisions de production d'énergie. Grâce à une ou plusieurs prévisions météorologiques, nous générons de multiples prévisions de production PV et nous construisons une combinaison linéaire de ces prévisions de production. La minimisation du Continuous Ranked Probability Score (CRPS) permet de calibrer statistiquement la combinaison de ces prévisions, et délivre une prévision probabiliste sous la forme d'une fonction de répartition empirique pondérée. Dans ce contexte, nous proposons une étude du biais du CRPS et une étude des propriétés des scores propres pouvant se décomposer en somme de scores pondérés par seuil ou en somme de scores pondérés par quantile. Des techniques d'apprentissage séquentiel sont mises en oeuvre pour réaliser cette minimisation. Ces techniques fournissent des garanties théoriques de robustesse en termes de qualité de prévision, sous des hypothèses minimes. Ces méthodes sont appliquées à la prévision d'ensoleillement et à la prévision de production PV, fondée sur des prévisions météorologiques à haute résolution et sur des ensembles de prévisions classiques.

Ensemble forecasting using sequential aggregation

for photovoltaic power applications.

Abstract:

Our main objective is to improve the quality of photovoltaic power forecasts deriving from weather forecasts. Such forecasts are imperfect due to meteorological uncertainties and statistical modeling inaccuracies in the conversion of weather forecasts to power forecasts. First we gather several weather forecasts, secondly we generate multiple photovoltaic power forecasts, and finally we build linear combinations of the power forecasts. The minimization of the Continuous Ranked Probability Score (CRPS) allows to statistically calibrate the combination of these forecasts, and provides probabilistic forecasts under the form of a weighted empirical distribution function. We investigate the CRPS bias in this context and several properties of scoring rules which can be seen as a sum of quantile-weighted losses or a sum of threshold-weighted losses. The minimization procedure is achieved with online learning techniques. Such techniques come with theoretical guarantees of robustness on the predictive power of the combination of the forecasts. Essentially no assumptions are needed for the theoretical guarantees to hold. The proposed methods are applied to the forecast of solar radiation using satellite data, and the forecast of photovoltaic power based on high-resolution weather forecasts and standard ensembles of forecasts.

État de l'art et contributions

L'amélioration des prévisions de production d'énergie photovoltaïque (PV) contribue à une meilleure intégration de l'énergie photovoltaïque. Les prévisions météorologiques, les modèles de conversion météo-production et les techniques de post-traitement statistiques constituent trois axes d'amélioration. Pour notre part, le problème est abordé de la façon suivante : un prévisionniste, souhaitant fournir des prévisions probabilistes de production PV, récupère des prévisions météorologiques (potentiellement de sources variées). Dans ce cadre général, de multiples méthodes peuvent être testées et combinées. Les prévisions météorologiques constituent la base de nos méthodes de post-traitement qui impliquent des ensembles de prévisions. Un état de l'art des modèles statistiques utilisés pour le PV est proposé par Bacher et al. [START_REF] Bacher | Online short-term solar power forecasting[END_REF] et Inman et al. [START_REF] Inman | Solar forecasting methods for renewable energy integration[END_REF]. La récente revue de l'état de l'art [START_REF] Antonanzas | Review of photovoltaic power forecasting[END_REF] inclut une revue des techniques de prévision probabiliste et de prévision d'ensemble appliquées pour le PV. Un nombre restreint de publications mentionnent des techniques de post-traitement prenant en compte un ensemble de prévisions météorologiques pour le PV [Zam+14 ; Ale+15 ; SAM16]. Par ailleurs, Lorenz et al. [START_REF] Lorenz | Short term forecasting of solar irradiance by combining satellite data and numerical weather predictions[END_REF] combinent des prévisions météorologiques et des prévisions issues de données satellitaires.

Un prévisionniste disposant de multiples prévisions peut souhaiter les combiner de façon optimale, par exemple par agrégation séquentielle. Le prévisionniste donne alors un poids à chaque prévision et fournit la combinaison linéaire des prévisions. Ces poids sont déterminés par une règle de mise à jour qui ne prend en compte que l'information du passé disponible à chaque pas de temps. De plus, la prévision fournie dispose de garanties théoriques de performance valables pratiquement sans hypothèses (sur un a priori, un processus stochastique ou une distribution sous-jacente), cf. Cesa-Bianchi et Lugosi [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF], la revue de Shalev-Shwartz [START_REF] Shalev-Shwartz | Online learning and online convex optimization[END_REF] ou l'introduction de Stoltz [START_REF] Stoltz | Agrégation séquentielle de prédicteurs : méthodologie générale et applications à la prévision de la qualité de l'air et à celle de la consommation électrique[END_REF] en français. Les garanties théoriques sont essentiellement des garanties théoriques de robustesse (borne de regret), assurant que la combinaison des prévisions est sur le long terme au moins aussi performante que la meilleure prévision ou la meilleure combinaison de prévisions à poids fixes. Ces techniques ont déjà été testées sur des jeux de données variés : consommation d'électricité, concentration d'ozone, champs de vent et de pression [Sto10 ; MSM09 ; Mal10 ; Bau15 ; GGN16]. Nous montrons que l'agrégation séquentielle permet d'améliorer les prévisions d'ensoleillement au chapitre 2, d'après Thorey et al. [START_REF] Thorey | Ensemble forecast of solar radiation using TIGGE weather forecasts and HelioClim database[END_REF]. Dans ce travail, nous étudions les ensembles de prévisions de 6 centres de prévision météorologique : China Meteorological Administration (CMA), European Centre for Medium-Range Weather Forecasts (ECMWF), MetOffice (UKMO), Korea Meteorological Administration (KMA), Centro de Previsao Tempo e Estudos Climaticos (CPTEC), et Météo France (M.-F.), pour un total de 158 prévisions. Ces prévisions sont comparées aux observations d'origine satellitaire de la base de données HélioClim-3, sous forme d'un champ d'observations sur la France métropolitaine. Les observations satellitaires de notre jeu de données ont une résolution de 1/12 • , qui est plus fine que la résolution de 0.25 • à laquelle nous récupérons les prévisions. Nous trouvons que les ensembles de prévisions sont généralement sous-dispersés. Nous utilisons l'agrégation séquentielle pour combiner linéairement les 158 prévisions avec des poids pouvant varier en temps et en espace. À chaque point de grille i et pour chaque pas de temps t, nous utilisons une régression ridge escomptée pour trouver les poids u (i) t minimisant une certaine distance entre les prévisions x (i) m,t et les observations y (i) t passées. Les poids sont choisis de sorte qu'ils minimisent

J(u) = λ u -w ref 2 2 + t-1 t =1 1 + γ (t -t ) 2 (y (i) t -u • x (i) t ) 2 ,
où w ref est un vecteur de réference, λ est le paramètre de régularisation et γ est le paramètre d'escompte. Cette méthode diminue l'erreur de prévision de 20% par rapport à la prévision de référence HRES du ECMWF. De plus, cette méthode corrige les biais locaux spatiaux et produit des motifs spatiaux plus réalistes de l'ensoleillement.

Dans la littérature météorologique, la minimisation du CRPS est une stratégie commune pour calibrer les prévisions probabilistes [Gne+05 ; Sch14 ; JMA15]. Cependant, les techniques classiques n'offrent pas de garantie théorique de robustesse et ont recours à des hypothèses sur les distributions. Par exemple, le « Bayesian model averaging » (BMA) fournit un mélange de distributions paramétriques, une somme de gaussiennes [START_REF] Gneiting | Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation[END_REF] ou de distributions gamma [START_REF] Sloughter | Probabilistic wind speed forecasting using ensembles and Bayesian model averaging[END_REF][START_REF] Sloughter | Probabilistic quantitative precipitation forecasting using Bayesian model averaging[END_REF]. La régression non-homogène cale les paramètres d'une distribution paramétrique en fonction des attributs d'un ensemble de prévisions [Gne+05 ; Wil09 ; TG10]. Ainsi la moyenne et la variance d'une distribution gaussienne sont calées par un modèle linéaire en fonction de la moyenne et de la variance d'un ensemble de prévisions.

Pour fournir des prévisions probabilistes, nous proposons une approche innovante fournissant un mélange de distributions [START_REF] Thorey | Online learning with the CRPS for ensemble forecasting[END_REF], au chapitre 3. L'originalité de notre technique provient de l'utilisation de règles de mise à jour des poids issues de l'agrégation séquentielle pour minimiser le CRPS de la distribution empirique pondérée. Grâce à l'utilisation de l'agrégation séquentielle, nos prévisions bénéficient de garanties théoriques de robustesse. L'agrégation séquentielle a déjà été appliquée avec succès aux prévisions quantiles [GGN16 ; BP11], mais pas pour fournir directement une prévision en loi, jusqu'aux travaux de cette thèse [START_REF] Thorey | Online learning with the CRPS for ensemble forecasting[END_REF], et de Baudin [START_REF] Baudin | Prévision séquentielle par agrégation d'ensemble: application à des prévisions météorologiques assorties d'incertitudes[END_REF] et Zamo [START_REF] Zamo | Post-traitements statistiques de prévisions de vent déterministes et d'ensembles sur une grille[END_REF]. Les relations entre les prévisions quantiles et les prévisions en loi sont détaillées dans la Section 1.3.

Par ailleurs, notre étude améliore la connaissance du biais du CRPS pour les distributions empiriques. Le score S(G E , y) évalue la qualité de la prévision en loi G E de l'observation y. Supposons que G E dépende de variables aléatoires, alors le score S(G E , y) est également une variable aléatoire. Le biais du score S(G E , y) est la différence entre le score moyen E(S(G E , y)) et le score de la distribution moyenne S(E(G E ), y). Nous généralisons les résultats du biais du CRPS calculé pour un ensemble uniforme de prévisions [START_REF] Ferro | On the effect of ensemble size on the discrete and continuous ranked probability scores[END_REF], et nous proposons une minimisation juste du CRPS, sans hypothèse sur des distributions sous-jacentes. Cela est possible en groupant les prévisions par classes de membres échangeables, de façon assez similaire au travail de Fraley et al. [START_REF] Fraley | Calibrating multimodel forecast ensembles with exchangeable and missing members using Bayesian model averaging[END_REF] pour le BMA.

L'extension des résultats mentionnés ci-dessus est proposée au chapitre 4 pour des pertes, autres que le CRPS, pouvant se décomposer par seuil ou par quantile. De telles décompositions sont étudiées récemment par Gneiting et Ranjan [START_REF] Gneiting | Comparing density forecasts using thresholdand quantile-weighted scoring rules[END_REF] et Ehm et al. [START_REF] Ehm | Of quantiles and expectiles: consistent scoring functions, Choquet representations and forecast rankings[END_REF]. Le sujet n'est toutefois pas nouveau pour l'évaluation des prévisions d'évènement binaire [START_REF] Shuford | Admissible probability measurement procedures[END_REF]. La question de l'incertitude des observations est abordée à la section 4.2. En assimilation de données, un compromis est généralement trouvé entre de nouvelles observations et une ébauche en fonction de leurs niveaux de confiance respectifs. Des observations possédant un fort niveau de bruit ne seront que peu prises en compte dans la mise à jour du vecteur d'état. Le cas des observations bruitées est également étudié par Yang [START_REF] Yang | Combining forecasting procedures: some theoretical results[END_REF] en agrégation séquentielle. La perte quadratique est normalisée par le niveau de bruit des observations, et une borne de regret par rapport à l'espérance des observations est obtenue sur les pertes simples et sur les pertes normalisées. Notre cadre est différent puisqu'il concerne la prévision probabiliste. Notre objectif premier n'est pas d'obtenir des bornes de regret en espérance, mais de prévoir au mieux la réalité en délivrant une prévision probabiliste. Au lieu de recevoir une observation sans bruit (la réalité), nous considérons que le prévisionniste reçoit de multiples observations ou une distribution. Nous développons cette idée par une technique de moindres carrés généralisés appliquée au CRPS. Des liens sont établis avec des statistiques de test.

Nos techniques sont testées sur des cas d'étude de prévisions en loi pour le PV décrites aux chapitres 5, 6 et 7. Notre approche est résumée par la figure A. Tout d'abord, nous fournissons des prévisions pour 219 parcs de production PV situés en France à la résolution temporelle de 30 min, jusqu'à un horizon de 6 jours. Dans ce cas d'étude, nous utilisons la prévision déterministe HRES et la prévision d'ensemble ENS du ECMWF, Arpège (la prévision déterministe de Météo France) et PEARP (la prévision d'ensemble issue de Arpège). Nous étudions la prévision probabiliste issue des modèles statistiques opérationnels de conversion entre les données de production et les prévisions météorologiques. Nous comparons les performances des ensembles de prévisions et des prévisions quantiles issues de HRES et Arpège. De plus, nous montrons une amélioration de la qualité des prévisions grâce à l'agrégation séquentielle des prévisions de production.

Dans un second cas d'étude, nous étudions l'utilisation du système de prévision AROME (prévision Météo France à haute résolution). La haute résolution spatiotemporelle de AROME permet de générer de multiples prévisions. En plus des incertitudes météorologiques, plusieurs modèles statistiques sont construits pour tenir compte de la difficulté de convertir des prévisions météorologiques en prévisions de production. Grâce à l'agrégation séquentielle, nous montrons qu'il est possible de fournir une prévision de production calibrée issue des prévisions AROME.

Les systèmes de production photovoltaïque insulaires de La Réunion et de la Corse sont étudiés dans un troisième cas d'étude, pour des horizons de prévision courts entre 30 min et 4 h. Nous explorons alors la possibilité de mises à jour infra-journalières grâce à la construction de prévisions prenant en compte toute l'information disponible de la journée. Nous travaillons avec les prévisions météorologiques usuelles ainsi qu'avec des prévisions issues de données satellitaires estimant le mouvement des masses nuageuses. Ces prévisions satellitaires se montrent particulièrement intéressantes pour des horizons de prévision inférieurs à 2 h, alors que les prévisions météorologiques usuelles sont plus utiles pour les horizons de prévision plus longs. Dans tous nos cas d'étude, nous montrons que nos techniques améliorent la qualité des prévisions selon des scores d'évaluation déterministes et probabilistes.

State of the art and contributions

Improved photovoltaic (PV) power integration needs better power forecasts. Forecasters may pursue efforts to improve meteorological models, weather-based power models or statistical post-processing methods. For our part, we focus on the following case: a forecaster, willing to provide probabilistic PV power forecasts, retrieves meteorological forecasts (possibly from various sources). In this general setting, numerous state-ofthe-art methods can be tested and combined. Weather forecasts are the basis of our post-processing methods, involving ensembles of forecasts. State-of-the-art statistical models for PV are described in Bacher et al. [START_REF] Bacher | Online short-term solar power forecasting[END_REF] and Inman et al. [START_REF] Inman | Solar forecasting methods for renewable energy integration[END_REF]. Also, the recent state-of-the-art review [START_REF] Antonanzas | Review of photovoltaic power forecasting[END_REF] includes a review of probabilistic forecasting and ensemble forecasting for PV applications. Only few publications mention postprocessing methods with an ensemble of weather forecasts for PV [START_REF] Zamo | A benchmark of statistical regression methods for shortterm forecasting of photovoltaic electricity production. Part II: Probabilistic forecast of daily production[END_REF][START_REF] Alessandrini | An analog ensemble for short-term probabilistic solar power forecast[END_REF][START_REF] Sperati | An application of the ECMWF Ensemble Prediction System for short-term solar power forecasting[END_REF]. Interestingly, Lorenz et al. [START_REF] Lorenz | Short term forecasting of solar irradiance by combining satellite data and numerical weather predictions[END_REF] combines weather forecasts and shortterm prediction from satellite data.

A forecaster having multiple forecasts hopefully wishes to combine them in an optimal way. To do so with online learning, also called sequential aggregation, the forecaster gives a weight to each forecast and delivers the combination of the forecasts. The combination rules stemming from online learning depend only on the available past information of each forecast step and come with theoretical performance guarantees under essentially no assumptions (concerning prior weights, underlying stochastic processes or distributions), see the monograph Cesa-Bianchi and Lugosi [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF], or the gentle reviews Shalev-Shwartz [START_REF] Shalev-Shwartz | Online learning and online convex optimization[END_REF] or Stoltz [START_REF] Stoltz | Agrégation séquentielle de prédicteurs : méthodologie générale et applications à la prévision de la qualité de l'air et à celle de la consommation électrique[END_REF] in French. The theoretical performance guarantees are essentially robustness guarantees (regret bounds), ensuring that the combination of the forecasts performs on the long run at least as well as the best forecast or the best fixed combination of forecasts. These techniques have already been tested for several applications: electricity consumption, ozone concentration, wind and geopotential fields [Sto10; MSM09; Mal10; Bau15; GGN16].

We show that online learning may be used to improve solar irradiance forecasts in Chapter 2 based on Thorey et al. [START_REF] Thorey | Ensemble forecast of solar radiation using TIGGE weather forecasts and HelioClim database[END_REF]. In this paper, maps of surface solar irradiance are forecasted using ensembles of forecasts from the THORPEX Interactive Grand Global Ensemble (TIGGE) with a 6-h timestep. First, we study ensemble forecasts from 6 meteorological centers: China Meteorological Administration (CMA), European Centre for Medium-Range Weather Forecasts (ECMWF), MetOffice (UKMO), Korea Meteorological Administration (KMA), Centro de Previsao Tempo e Estudos Climaticos (CPTEC), and Météo France (M.-F.), providing a total of 158 forecasts. The forecasts are compared with observations derived from MeteoSat Second Generation (MSG) and provided by the HelioClim-3 database as gridded observations over metropolitan France. Satellite estimations of our data set have a resolution of 1/12 • which is much finer than the 0.25 • resolution at which the forecasts are retrieved. We find that the ensemble forecasts are generally under-dispersed, even when grouped together. Secondly, we use online learning techniques to linearly combine all 158 forecasts with weights that vary in space and time. For each grid point i and for each time t, we apply a discounted ridge regression to find the weights u (i) t minimizing a certain distance between the past forecasts x (i) m,t and the past observations y (i) t . The weights u (i) t are chosen as the minimizer of

J(u) = λ u -w ref 2 2 + t-1 t =1 1 + γ (t -t ) 2 (y (i) t -u • x (i) t ) 2 ,
where w ref is a reference vector, λ is the regularization parameter and γ is the discount parameter. This method decreases the forecast error by 20% compared to the reference forecast HRES from ECMWF. Besides, this method also corrects the spatial local biases and produces a more realistic spatial pattern of predicted irradiance.

In the meteorological literature, minimizing the CRPS is a common strategy to calibrate probabilistic forecasts [Gne+05; Sch14; JMA15]. However, standard techniques do not offer theoretical guarantees of robustness and usually resort to strong assumptions on the distributions. For example, Bayesian model averaging (BMA) techniques provide a mixture of parametric distributions, usually a Gaussian sum [START_REF] Gneiting | Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation[END_REF] or gamma distributions sum for wind and precipitation applications [START_REF] Sloughter | Probabilistic wind speed forecasting using ensembles and Bayesian model averaging[END_REF][START_REF] Sloughter | Probabilistic quantitative precipitation forecasting using Bayesian model averaging[END_REF]. Non-homogeneous regression fits the parameter of a parameterized distribution using characteristics of the ensemble of forecasts [START_REF] Gneiting | Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation[END_REF][START_REF] Wilks | Extending logistic regression to provide full-probabilitydistribution MOS forecasts[END_REF][START_REF] Thorarinsdottir | Probabilistic forecasts of wind speed: ensemble model output statistics by using heteroscedastic censored regression[END_REF]. For instance, a Gaussian distribution is fitted using a linear model between the mean of the distribution and the mean of the forecasts. Likewise the standard deviation of the Gaussian distribution is fitted according to the ensemble spread.

To provide probabilistic forecasts, we propose an innovative approach by combining multiple forecasts in a linear opinion pool [START_REF] Thorey | Online learning with the CRPS for ensemble forecasting[END_REF], in Chapter 3. The originality of our technique is to use combination rules deriving from online learning techniques in order to minimize the CRPS of the weighted empirical distribution function. Because we use online learning techniques, our forecasts come with theoretical guarantees of robustness. Online learning techniques have already been applied successfully to quantile prediction [GGN16; BP11], but not for directly delivering probabilistic forecasts, until this thesis [START_REF] Thorey | Online learning with the CRPS for ensemble forecasting[END_REF], and the works of Baudin [START_REF] Baudin | Prévision séquentielle par agrégation d'ensemble: application à des prévisions météorologiques assorties d'incertitudes[END_REF] and Zamo [START_REF] Zamo | Post-traitements statistiques de prévisions de vent déterministes et d'ensembles sur une grille[END_REF]. Relationships between quantile losses and probabilistic forecasts are detailed in Section 1.3.

Our study also improves the knowledge of the CRPS expectation for linear opinion pools. The score S(G E , y) evaluates the quality of the probabilistic forecast G E for the observation y. Assume G E relies on random variables, the score S(G E , y) is also a random variable. The bias of the score S(G E , y) is the difference between the average score E(S(G E , y)) and the score for the average distribution S(E(G E ), y). We generalize results on the bias of the CRPS computed with ensemble forecasts [START_REF] Ferro | On the effect of ensemble size on the discrete and continuous ranked probability scores[END_REF], and propose a new scheme to achieve fair CRPS minimization, without any assumption on the distributions. This is achieved by grouping forecasts in classes of exchangeable members, quite similarly to the work of Fraley et al. [START_REF] Fraley | Calibrating multimodel forecast ensembles with exchangeable and missing members using Bayesian model averaging[END_REF] for BMA. Extension of the above results is proposed in Chapter 4 to other losses than the CRPS, which admit a threshold decomposition or a quantile decomposition. Such decomposition was recently studied by Gneiting and Ranjan [START_REF] Gneiting | Comparing density forecasts using thresholdand quantile-weighted scoring rules[END_REF] and Ehm et al. [START_REF] Ehm | Of quantiles and expectiles: consistent scoring functions, Choquet representations and forecast rankings[END_REF]. This topic is not new for the evaluation of binary event forecasts [START_REF] Shuford | Admissible probability measurement procedures[END_REF]. The question of uncertainty in the observations is addressed in Section 4.2. In data assimilation, a trade-off is commonly achieved between observational information and the background state based on their respective level of confidence. Observations with a high level of noise will not have much importance for the state update. Observational noise is studied by Yang [START_REF] Yang | Combining forecasting procedures: some theoretical results[END_REF] for online learning. They normalize the square loss by the observation noise rate and obtain regret bounds on the normalized and unnormalized losses computed against the distribution mean. Our work is slightly different since we intend to provide probabilistic forecasts. Our first objective is not to obtain regret bounds in expectation, but to forecast as best as possible the reality by delivering a probabilistic forecast. Instead of receiving one single noiseless observation (the reality), we consider that the forecaster receives multiple observations or a distribution. We develop this idea using a generalized least-square method applied to the CRPS and provide connections to test-statistics.

Our new techniques are tested in several case studies of probabilistic PV power forecasts described in Chapters 5, 6 and 7. Our approach is summarized in Figure B. First, we provide forecasts for 219 photovoltaic (PV) power plants located in France with a 30-min timestep, up to 6 days of lead time. In this case study, we use ECMWF deterministic forecast HRES and ensemble ENS, Arpège (Météo France deterministic forecast) and PEARP (Météo France ensemble forecast). We investigate PV probabilistic forecasting with the operational statistical models between the production data and the weather forecasts. We compare the predictive performance of ensemble forecasts and quantile forecasts deriving from the deterministic HRES and Arpège. We also provide improvements on the calibration of the probabilistic forecasts by combining the production forecasts.

In a second case study, we investigate the use of AROME forecasts (Météo France high-resolution forecasts). The high spatio-temporal resolution of AROME enables to generate multiple forecasts accounting for various local scenarios. In addition to weather uncertainties, several statistical models are built to account for the difficulty of converting weather forecasts to production forecasts. Using online learning techniques, we show that we can provide calibrated production forecasts using the single highresolution forecasting system AROME.

The insular PV power production of Réunion and Corsica are investigated in a third case study focusing on short lead times from 30 min to 4 h. Here we explore the possibility of intraday updates with the generation of new members incorporating all the available information of the day. Namely, new forecasts are built using regular PV forecasts and available production data of the current day as inputs. We work with usual weather forecasts as well as predictions from satellite data derived from cloudmotion analysis. Prediction using satellite data are most useful for short lead times below 2 h, while regular weather forecasts deliver the most useful information for longer lead times. For all the case studies, we show that our technique provides forecast improvements for multiple evaluation tools both deterministic and probabilistic.

Introduction

This chapter aims at providing background knowledge on this thesis. First, we introduce photovoltaic power forecasts using weather forecasts, from the industrial context to statistical learning challenges. The theoretical background of sequential aggregation, or online learning with multiple experts, is then detailed. We also delineate our position on probabilistic forecasting and introduce probabilistic forecasts evaluation.

Forecasting photovoltaic power production with meteorological forecasts 1.Context

This thesis aims at improving forecasts of photovoltaic (PV) power production. To be more specific, the company "Electricité de France" (EDF) is particularly interested in forecasting PV power for several geographical areas, in metropolitan France or for insular electric systems. For example, we model an estimation of the total production of Metropolitan France, the total production of Corsica Island and the total production of Réunion. The installed capacity of renewable electricity is growing in France. The need for forecasts of renewable energy production grows with the amount of renewable energy delivered to the electrical grid. Over the past decade, the installed capacity of wind and solar power drastically increased. The current annual growth rate is close to 10%, see Table 1.1. Indeed, the installed capacity of solar power increased in France by 899 MW in year 2015 and 576 MW in 2016, to reach the amount of 6.8 GW. Over the same periods, the installed capacity of wind power increased in France by 1011 MW in year 2015 and 1345 MW in 2016, to reach the amount of 11.67 GW. This strong growth is significant compared to the total installed capacity of 129 GW at the beginning of 2016, with above 63 GW of nuclear capacity and above 25 GW of hydropower capacity. These numbers need to be put into perspective, because wind and solar are intermittent energy sources. One gigawatt of PV power is not equivalent to one gigawatt of nuclear power, at least in terms of load factor. The load factor is the ratio of the produced energy and the energy that would have been produced if the production level was maximal. In fact, the load factor of wind and solar power are quite low compared to the load factor of nuclear power plants (around 75%). The average load factor of solar power is of 15% with a strong seasonality (6% in December and 20% in summer). Besides, the load factor of wind power is around 24% (30-33% in winter and 15-20% in summer). We refer the interested reader to the periodic report "Overview of Renewable Electricity" * (in French "Panorama de l'électricité renouvelable"), and to the "Forecast assessment of electricity supply-demand balance" † (in French "Bilan prévisionnel de l'équilibre offredemande d'électricité en France"), both published by the French Transmission System Operator "Réseau de Transport d'électricité" (RTE).

What is at stake in renewable energy forecasting ?

The production forecasts are achieved at various geographical and temporal scales, depending on the production target. For example, the production can be local at the scale of a production unit, or at a larger geographical scale, i.e. a sum of local productions.

For operational purposes, Transmission System Operators and Distribution System Operators need production forecasts to anticipate constraints on the electrical grid, generated by supply-demand imbalances. For financial purposes, Balance Responsible Entities (BR) also need production forecasts. A BR must declare its injections and extractions of power on the grid, with constraints on its so-called balance perimeter.

Control parameters

Model Simulations Predictions

Assimilation A priori knowledge

Observations

Figure 1.1 -Illustration of data assimilation.

These constraints ensure the whole system to be balanced. When injections and extractions do not match the BR declaration, the BR suffers financial penalties. The BRs can trade electricity on the market to ensure the balance of its perimeter or take benefits from the market price. Therefore, improved production forecasts on their perimeter means better market positioning and avoided penalties for BRs. Historically in France, EDF was the main BR with solar and wind production inside its balance perimeter due to the mechanism of feed-in tariffs. The trends of the market evolution is that wind and solar energy producers will introduce their own product on the market, or sell it directly to a BR. Renewable energy forecasts can also be integrated to maintenance planning. The motivation behind is that the maintenance should be done -when the production is low for renewable energy producers, -without generating constraints on the grid for grid operators.

Weather forecasting

Wind and solar power productions are strongly dependent on the weather, hence wind and solar forecasts resort to numerical weather predictions (NWPs) for lead times between a few hours to several days.

How numerical weather predictions are generated ? NWPs are produced with data assimilation schemes, the most famous being the Kalman filter [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]. Data assimilation integrates several elements: a state vector, a model, a priori knowledge, observations, an observation function and control parameters, as represented in Fig- of data assimilation is the best evaluation of the state vector with the available information. After the reception of new observations, the state vector is updated to take into account this new information, given characteristics on observations errors, background errors on the state and model errors.

For numerical weather predictions:

-The state vector describes the atmospheric state and includes several fields such as horizontal wind speed and temperature.

-Atmospheric models are built from partial differential equations representing the atmospheric motions and many simplifications (or parameterizations), which tend to model the effect of a process rather than the process itself. The model is run from a given state to create predictions of the evolution of this state.

-The observation function relates the state space to the observation space. Indeed, the observations generally do not live in the same space as the state. For example, satellite images are used as observations, but do not always have direct counterparts in the model state.

-The size of the observation and state vectors are huge. They easily reaches 10 6 for the observations and 10 9 for the state vector, hence computational problems arise.

Numerical weather predictions used in this thesis.

In practice, several meteorological forecasting systems are used in this thesis. Their current characteristics are summarized in Table 1.2. Besides the model and the data assimilation scheme, the forecasts may differ in their covered geographical areas, their lead time, and their spatial and temporal resolutions. Ensemble forecasts are traditionally generated with the same model as a deterministic forecast, but at a coarser resolution with slight changes in the model physics, in the parameterizations, or with different perturbations of the initial conditions. Therefore ensemble forecasts integrate different sources of errors and quantify the forecast uncertainty. All along this thesis, we show the benefits of using several forecasts. This approach begins with the variety of weather forecasts that we use as inputs.

Météo France generates the global model Arpège, from which the 34 ensemble forecasts PEARP are derived. Météo France also provides mainly for France the local model AROME, designed to forecast severe weather events such as heavy rains in the south of France. The resolutions of AROME are quite fine (currently 1. the cost of limited geographical extension and limited lead time of two days. We also use atmospheric forecasts from ECMWF: the deterministic forecast HRES and the ensemble of forecasts ENS. For solar forecasting with a 6-h temporal resolution, TIGGE (THORPEX Interactive Grand Global Ensemble) ensembles from several meteorological centers were used. The detailed description of TIGGE ensembles is postponed to Chapter 2.

Solar radiation forecasting

For PV applications, the main meteorological fields are irradiance (or solar radiation) and 2-m temperature. Irradiance is directly related to the energy produced by the solar panels and high temperature at ground-level decreases the panel efficiency.

The position of the Sun compared to the position of the Earth controls directly the top-of-atmosphere irradiance. Along one year, the evolution of these positions generates the diurnal and the seasonal cycles. The diurnal cycle facilitates solar forecasting, since we know that every day begins with a null level of radiation. For a clear sky day, the main challenge is to estimate the bell-shape of solar radiation and the maximal reached value. The seasonal cycle is clearly depicted with monthly maximal value of solar radiation, see Figure 1.2. The ratio between the maximal value of January and June almost reach a factor 4. We highlight the fact that solar radiation data are timeaveraged, and usually indicate the average flux over the past period. For example, solar radiation at 12:00 with a 3-h temporal resolution is the average flux between 9:00 and 12:00.

From top-of-atmosphere to ground-level, solar radiation encounters many processes (absorption, scattering, reflection) mainly occurring in clouds. Hence irradiance maps are closely related to nebulosity maps. Cloud cover and nebulosity are often used as a proxy for solar radiation. Solar radiation maps of HelioClim (real-time estimations from satellite data), AROME and HRES forecasts (12-h lead time) are shown in Figure 1.3, for 2013-06-13 at 12:00 (UTC). Note that HelioClim and AROME maps are 1-h averages while the HRES map is a 3-h average. The spatial variations of irradiance are quite high since very cloudy areas may receive less than one third of the energy received by clear-sky areas. These variations are seen by the forecasts at the scale of a few tens of kilometers, or even less. The tremendous impacts of orographic effects in mountains areas with many scattered clouds can be seen in Figure 1.3. We see that spatial resolution is a key-issue in solar forecasting. At the scale of a gridpoint, one may ask whether the level of irradiance is correctly forecasted. However, at the scale of a regional map, one may ask whether the clouds are modeled at the correct time and location. Indeed, a cloud may be "seen" by the model, but located with an error of 50 km. As a conclusion, we note that:

-Observations may be local (ground-station measurements) or estimated with satellite data. Satellite data provide a better spatial representation, but may be more prone to calibration errors than local measurements.

-We need a validation criterion to assess the quality of a forecast. This validation criterion may depend on the spatial scale of the validation.

-The forecasted average over several grid points may be the best estimation for comparisons at the grid-point scale.

Using high-resolution forecasts such as AROME is a delicate task, see our contribution in Chapter 6.

PV power data and statistical modeling

As opposed to irradiance maps, production data of a PV power plant are localized. For lead times of several hours to several days, we usually work at a resolution of 30 min or 1 h. Averaging in time production data smooths the variability of the data, as shown in Figure 1.4, and facilitates the work of the forecaster. Indeed, the local production peaks are barely seen after time-averaging. The forecaster faces a dilemma when the production data show a large uncertainty. Should they attempt to forecast the average, q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 5 10 15 20 PV power q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 30min average 1h average 2h average or should they attempt to forecast local peaks although the prediction may not be at the correct time ? Weather forecasts are converted to PV forecasts with statistical or physical models. Besides weather variables, PV power production rely on the technology of the solar panels and on the incident angle of solar radiation on the solar panels. Hence physical models include the panels orientation and inclination. Solar radiation data are a good starting point to build a statistical model. Using solar radiation estimates from HelioClim (and not forecasts), we see in Figure 1.5 that linear models are a good approximation of the relationship between solar radiation and PV production, for a period of nearly 20 days. Here we simply applied a multiplicative factor to the solar radiation estimates. Advanced models often use clear-sky production profiles and clear-sky radiation profiles computed for each day of the year. In this setting, weather forecasts provide indications on possible production decreases along the day. The production P and solar forecasts I are converted respectively to the clear-sky indices τ P and τ I :

P = τ P P cc and I = τ I I cc ,
where clear sky production P cc and clear sky solar radiation I cc are the production and solar radiation in clear sky conditions. A statistical regression estimates the parameters a i between weather forecasts and the forecasts τ P of clear-sky index:

τ P = a 0 + a 1 τ I + a 2 τ 2 I + a 3 T cc .
In this example, the non-linear term τ 2 I and the total cloud cover T cc are added as explanatory variables. The statistical models of our PV forecasts are detailed in Chapter 5.

In practice, the forecaster encounters major difficulties:

-The forecaster cannot determine the best set of parameters for the current period, but only for a past period. This may result in drifting effects or inaccurate parameter estimates due to inter-year variability.

qqqqqqqqqqqq qq qq q q q qqq qq q q q q q q qq q q q q qqqqqqqqqqqqqqqqqqqqqq q q qq q q q qq q q q q q q q qqq q q qq q qqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q qqqq q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q qq q qq q qq qqqq qq qqqqqqqqqqqqqqqqqqqqqq q q q q q qq q q q q q q qq q q q q q q q qq qqqqqqqqqqqqqqqqqqqqqq q qqq q q qq qq q q q q q q q q q qq q qq qqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q qqqq q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqq q q qq q q q qqq q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q qq qqq q qq q q q q q q q q qqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q q qq q q q q qqq q q qqqqqqqqqqqqqqqqqqqqqq q q q q q q q qq qqqq q q q q q q q q q qq qqqqqqqqqqqqqqqqqqqqq q q q q q q q q qq q q q q q q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q qqqqq q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqq q q qq q q q q q qqq qqq q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqqq q q q q q q qq q qq qqq q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q qq qqq qq q q q q q q q q q qqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q qqqq q q q q q q q q q q q qqqqqqqqqqqqqqqqqqqqq q q q q q q q q q q q q q qqq qq q q qq q qqq qqqqqqqqqqq Figure 1.5 -Time-series of Helioclim (red line) and production data (black dots), with simple scaling of the data (for one plant in April 2013).

-Solar radiation forecasts may be quite far from the reality. Systematic biases such as seasonal biases may be corrected with statistical models taking into account the hour of the day, the time in the year and the lead time of the forecasts. However, specializing the statistical model increases the amount of model parameters, or reduces the available amount of data to fit several models.

-PV power observations are not always available, due to data corruption or lack of metering.

Interestingly, the forecaster can use statistical methods as downscaling methods. For example, statistical models are built between forecasts at a 3-h resolution and observations at a 30-min resolution in order to provide forecasts at the 30-min resolution. Thanks to the fine resolution of the observations, it is possible to go under the resolution of the forecasts. The same principle is applied by using HelioClim as solar radiation observations to improve the spatial resolution of the forecasts in Chapter 2.

Probabilistic forecasts of PV power with meteorological forecasts

Using multiple grid points gives a good insight into the spatial variability in the irradiance maps, see Figure 1.6. We show production data and solar forecasts from AROME, normalized by the daily maximum value. In a clear-sky situation (1.6(a)), the forecasts of nearby grid-points are almost identical. But in a situation with large uncertainties on the cloud cover (1.6(b)), the spatial variability of the solar forecasts is very high. For this particular day, the information of cloud cover variability is correct and noticeable in the production data. Our point here is to illustrate that using only one single forecast may not be appropriate, and that solar forecasts integrate complex spatio-temporal information.

Probabilistic forecasts acknowledge the inability of the forecaster to provide a priori a perfect estimate of the observation, that we consider here noiseless. Probabilistic forecasts are conveniently described by a Cumulative Distribution Function (CDF). They give the probability (according to the forecaster) that the observation is below a threshold. Confidence intervals are derived from the CDF of the forecaster. For example, the forecaster says that there is 5% chance that the production is below the number a and 95% chance that the observation is below the number b. The interval [a, b] is then a 90% confidence interval.

We illustrate probabilistic forecasts for PV applications with a 30-min temporal resolution in Figure 1.7. For two consecutive days, the forecaster delivers a probabilistic forecast and a deterministic forecast. The median of the probabilistic forecasts mainly Time (hours) PV power and solar forecasts q q q q q q q q q q q q q q q q q q q q q q q (a) 2013-04-15

5 10 15 20
Time (hours) PV power and solar forecasts q q q q q q q q q q q q q q q q q q q q q q q (b) 2013-04-08 indicates a change in the level of production between the two days. With the median only, one cannot know whether a large intraday variability is expected. This information is clearly given by the confidence intervals, which are much larger for the first day than for the second day. In this example, the observed production has indeed a larger intraday variability during the first day. The large confidence intervals reflect the inability of the forecaster to describe precisely the observation variations. Besides, the deterministic forecast provides additional information, such as production decreases during the first day. The integration of this information is possible thanks to sequential aggregation, as described in Section 1.2.

PV power

Time

Figure 1.7 -PV production (black), deterministic forecast (blue) and probabilistic forecasts (confidence intervals in shaded gray, median in white) for two consecutive days.

Sequential aggregation

Context

In this section, we introduce online learning in the expert setting. Just before time t > 0, a forecaster aims to deliver the best possible prediction for the observation y t ∈ Y, while an advice of M experts x m,t ∈ X is given to the forecaster. In other words, the forecaster receives several forecasts (or experts) and wishes to make an optimal use of this information. In real-world applications, the forecaster may also build experts and not only receive them. The key-point of sequential aggregation is nevertheless to ensure performance guarantees, given a set of arbitrary forecasts.

We work in the deterministic setting of individual sequences : no assumptions are made on the observations y t or on the experts x m,t . This setting is particularly interesting, because the performance guarantees described below are ensured no matter the received data (observations and forecasts). The methods providing performance guarantee with individual sequences are therefore intrinsically robust, see Cesa-Bianchi and Lugosi [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF] for an in-depth analysis.

In practice, the forecaster gives the weight u m,t to the expert x m,t , and provides the forecast s t (u t , x t ), where the function s t can be arbitrarily chosen by the forecaster. A simple combination of forecast s t (u t , x t ) = u t x t = M m=1 u m,t x m,t is often chosen. The general algorithm reads:

Initialization: u 1 ;
For each time index t = 1, 2, ..., T 1. get the vector of predictions data x t , 2. compute the forecaster's choice s t (u t , x t ), 3. get the verification y t and compute u t+1 , based on the update rule.

The initial weight vector u 1 is arbitrarily set, e.g., to [1/M, . . . , 1/M ] . Several examples of update rules are given in Section 1.2.3. We highlight here the fact that X and Y may be functional spaces. Hence it is possible to combine probabilistic forecasts (or CDFs) as described in Section 1.2.4.

Algorithm evaluation with regret bounds

The forecaster is willing to make the best possible forecast, where "best" refers to a real-valued loss function , measuring the distance between the observation and the forecast. The loss is negatively oriented, meaning that it should be as low as possible. The loss t (u) for time t is defined with implicit relationship with the experts x m,t and the observation y t . This notation is used because the weight vector u is the forecaster's choice while the experts and observations are arbitrarily given to the forecaster.

The expression "best possible forecast" introduces the notion of algorithm evaluation. The total loss of the forecaster may be divided into two terms:

T t=1 t (u t ) = inf u∈U T t=1 t (u) approximation error + T t=1 t (u t ) -inf u∈U T t=1 t (u)
estimation error , of approximation and estimation errors. The approximation error is the error of the best possible forecast with fixed weights in U. While the approximation error is considered as being unavoidable, the estimation error should be as low as possible. The estimation error is referred to as the cumulated regret of the forecaster, who competes against the best combination of experts. The cumulated regret of an algorithm A, generating the weights u t , is defined by

R T (A) = T t=1 t (u t ) -inf u∈U T t=1 t (u) .
Standard algorithms show a sublinearity of R T in T , expressed by lim T →∞ R T o(T ). This guarantees that the combination of the experts defined by the weights u t is asymptotically at least as good as the best fixed combination of experts. Consequently, the combination of experts is also asymptotically at least as good as the best experts and the best fixed subset of experts with uniform weights.

Online learning algorithms

In this section, we introduce several algorithms to find the weights u t . The update rules are often described by a regularized regression. The weight vector u t is then the minimizer of a cost function, with general formulation

u t = arg min u∈U        Ψ 0,t (u) + t-1 s=1 Ψ s,t (u, u s ) Regularization + t-1 s=1 s (u) Loss       
To obtain a general expression, we noted above s (u) instead of s (u) to indicate that the terms s (u) are eventually replaced by a linear or a quadratic approximation s (u) of s (u), for example around u s [KW97; HAK07; McM11]. The losses s (u) are to be minimized with a trade-off against the regularization terms Ψ s,t to control the variations of u t from the previous weights u s or from a reference weight vector u ref (or prior) implicitly defined in Ψ 0,t . For example, regularization terms can be of 2-norm

type u -u ref 2 2 , of Kullback-Leibler type m u m ln( um u ref m ) or of χ 2 type m (um-u ref m ) 2 u ref m [KW97]
. The χ 2 regularization may be seen as a second-order approximation around u ref of the Kullback-Leibler regularization.

Ridge Regression. A celebrated example is the ridge regression for the square loss ((u x t ) -y t ) 2 . The weight vector u t ∈ R M is chosen as the minimizer of

J(u) = λ u 2 2 + t-1 t =1 (y t -u x t ) 2 .
The parameter λ controls the trade-off between the losses and the 2-norm regularization.

The ridge regression has the following regret bound:

T t=1 (y t -u t x t ) 2 -min u∈B M T t=1 y t -u x t 2 O (ln T ) ,
under the assumption of bounded losses (y t -y t )

u m,t+1 = u m,t exp (-η m,t ) M m =1 u m ,t exp -η m ,t
.

The weight u m,t+1 is derived from u m,t with a multiplicative factor, accounting for the loss m,t . The idea is that the weight of a good expert (with a small loss compared to the other losses) increases with time, and the learning rate controls the speed of the variation of u t . Because the losses are linear, the best fixed combination of expert is the best expert: inf u∈U

T t=1 M m=1 u m m,t = min k T t=1 k,t .
In other words, the forecaster competes against the best expert. The algorithm EWA admits the following regret bound:

sup T t=1 M m=1 u m,t m,t -min k T t=1 k,t ln M η + η (b -a) 2 8 T .
The algorithm EWA may also be used with any convex function (u), thanks to Jensen's inequality ‡ . Using the bound for linear losses and Jensen's inequality, we have

T t=1 (u t ) -min k T t=1 k,t T t=1 M m=1 u m,t m,t -min k T t=1 k,t ln M η + η (b -a) 2 8 T .
We see that a limitation of EWA is that the forecaster competes against the best expert, but not again the best fixed combination of experts. The bound is minimized for the ‡. For a convex function ψ and a distribution D, we have ψ(Ez∼D(z)) Ez∼D(ψ(z)). Roughly speaking, the value of the expectation is inferior to the expectation of the values.

optimal learning rate η = (b -a) -1 8(ln M )/T , and the minimal value of the bounds, reached with η , is equal to (b -a) (ln M )T /2 = o(T ). A forecaster may use several learning rates according to the length of the experiment in order to circumvent the dependency in T of η . This is referred to as the doubling trick or its variants [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF][START_REF] Baudin | Prévision séquentielle par agrégation d'ensemble: application à des prévisions météorologiques assorties d'incertitudes[END_REF].

Exponentiated gradient (EG). Using EWA and the so-called gradient trick, the forecaster can also compete against the best convex combination. The resulting algorithm is Exponentiated Gradient, defined by the update rule

u m,t+1 = u m,t exp -η m,t M m =1 u m ,t exp -η m ,t .
Strictly speaking, EG is built by replacing m,t by the loss gradient

m,t = ∂ t ∂u m (u t ) ,
in the update rule of EWA.

The gradient trick uses the convexity of to obtain:

t (u t ) -t (u) (u t -u) ∇ t (u t ) = u t t -u t .
for any two vectors u t , u ∈ P M . Summing over time, we get the following regret bound inequalities:

T t=1 t (u t ) -inf u∈P T t=1 t (u) = sup u∈P T t=1 t (u t ) -t (u) sup u∈P T t=1 M m=1 u m,t m,t -u m m,t = T t=1 M m=1 u m,t m,t -min k T t=1 k,t ln M η + η ( b -a) 2 8 T .
The regret bound of EWA with the loss gradients gives the last inequality, with the assumptions of bounded gradients m,t ∈ [ a, b].

Polynomially weighted averages with multiple learning rates (ML-Poly).

The algorithm ML-Poly [START_REF] Gaillard | A Second-order Bound with Excess Losses[END_REF], described in Table 1.3, enjoys two desirable properties: it has no parameters and adapts to the difficulty of the data. The algorithm relies on the excess loss u t tm,t reflecting the regret of the forecaster to play the weighted combination u t instead of the mth expert. Designed for linear losses, the algorithm ML-Poly may be used with convex loss functions thanks to the gradient trick. The regret bound of ML-Poly is expressed against the best member for the linearized losses. For all sequences of losses m,t ∈ [0, 1], the cumulated loss of ML-Poly is bounded:

T t=1 u t t min 1 m M    T t=1 m,t + M (1 + ln(1 + T )) 1 + T t=1 (u t t -m,t ) 2    .
update the learning rate of each member η m,t = 1/ 1 + x -y

Loss

-1 0 1

Figure 1.8 -Deterministic losses: quadratic loss (black), quantile loss of level 0.9 (gray), log loss for y ∈ {0, 1} (dotted black).

As opposed to the aforementioned regret bounds, the bound of ML-Poly is of secondorder due to the term T t=1 (u t tm,t ) 2 . The interested reader is referred to Gaillard et al. [START_REF] Gaillard | A Second-order Bound with Excess Losses[END_REF] for a detailed analysis of second-order bounds, and to Koolen and Van Erven [START_REF] Koolen | Second-order Quantile Methods for Experts and Combinatorial Games[END_REF], Luo and Schapire [START_REF] Luo | Achieving All with No Parameters: AdaNor-malHedge[END_REF], and Wintenberger [START_REF] Wintenberger | Optimal learning with Bernstein online aggregation[END_REF] for further examples of algorithms showing second-order bounds. These bounds are adaptive in the sense that the algorithm performs well in the worst-case, but shows improved guarantees for easy data. The worst case scenario gives a bound O( √ M T ln T ), indicating that even in the worst case, the weighted forecast will perform at least as well as the best forecast. Besides, in the case of i.i.d. sequences of losses, with one expert significantly better than the others, the regret bound is practically constant.

Examples of loss functions

In the most common approach, the forecaster provides a point forecast (X , Y ⊂ R). Noting s(u, x) = y, the following losses are standard examples:

-Square loss: ( y -y) 2 .

-Quantile loss of level α: (H( y -y) -α) ( y -y), where H is the unit step function (H(x) = 1 if x > 0 and zero otherwise). Let y be a random variable described by the CDF F, the quantile loss of level α is minimized in average over the distribution of y by F -1 (α). The absolute loss |s(u, x) -y| equals twice the quantile loss of level 0.5.

-Log loss: -y ln( y) -(1 -y) ln(1 -y); often used in classification with X = [0, 1] and Y = [0, 1] or Y = {0, 1}.

See Figure 1.8 for a graphical representation of the losses.

In a probabilistic forecasting, the probability density function (PDF) g and the related CDF G are delivered, while the value y ∈ R is observed. For linear models G = m u m G m , the forecaster delivers a so-called linear opinion pool. In this setting, standard examples are -Quadratic score: -2g(y) + g 2 .

-Logarithm loss: -log(g(y)).

-Hyvärinen score: 2 g (y) g(y)

g (y)

g(y) 2 .
-Continuous ranked probability score (CRPS): (G(θ) -H(θ -y)) 2 dθ.

In our contributions, we use the square loss in Chapter 2 to improve deterministic forecasts of solar radiations, and the CRPS in Chapter 3 from a theoretical perspective in online learning and in Chapter 5 for practical applications with PV forecasts. The last section of this introduction is dedicated to a specific type of losses for probabilistic forecasting, which includes the CRPS.

Probabilistic forecasting with non-local strictly proper scoring rules

In this section, we introduce non local strictly proper scoring rules and explain why we resort to them.

Strict propriety. We are interested in a strictly proper scoring rule in the context of probabilistic forecasting, see Gneiting and Raftery [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF]. Say a forecaster provides a probabilistic forecast to predict y. We consider that y is a random variable described by the CDF F. The quality of the forecast is measured by a score S(G, y), where G is a CDF describing the probabilistic forecast. The scoring rule is said to be proper if E y [S(F, y)] E y [S(G, y)] and strictly proper if the strict inequality holds. In other words, with a strictly proper scoring rule, the average score is minimized only by the correct distribution F.

Locality. Let {A 1 , . . . , A K } be a partition of K possible events. By definition, two events A i and A j cannot occur simultaneously. The forecaster delivers p = (p 1 , . . . , p K ) reflecting an opinion on the probability of occurrence of each event. If the event A i occurred, a local scoring rule of p depends only on p i . In other words, a scoring rule is said to be local if the score only depends on the predicted probability of the event that actually occurred. Locality of the scoring rule is not always a desired property [START_REF] Bröcker | Scoring probabilistic forecasts: The importance of being proper[END_REF]. The logarithm score was shown to be a natural choice of strictly proper local scoring rule in the non-binary case (K > 2), see Bernardo [START_REF] Bernardo | Expected information as expected utility[END_REF] and Benedetti [START_REF] Benedetti | Scoring rules for forecast verification[END_REF] for example.

Why do we use non local strictly proper scoring rules ? Our objective is to provide probabilistic forecasts. To do so, we work with an ensemble of forecasts, which is conveniently described by a CDF. The corresponding probability density function is a combination of Dirac distributions, reaching 0 almost everywhere. The events "y is between θ k and θ k + dθ k " are therefore poorly predicted, and local scoring rules are hardly applicable in our case.

We introduce the case of binary events in Section 1.3.1, and probabilistic scores for continuous variables in Section 1.3.2. We follow here Shuford et al. [START_REF] Shuford | Admissible probability measurement procedures[END_REF], Schervish [START_REF] Schervish | A General Method for Comparing Probability Assessors[END_REF], and Buja et al. [START_REF] Buja | Loss functions for binary class probability estimation and classification: Structure and applications[END_REF] to give a simple explanation on the construction (but not the choice) of non-local strictly proper scoring rules. A discussion on the choice of strictly proper scoring rule may be found in Merkle and Steyvers [START_REF] Merkle | Choosing a strictly proper scoring rule[END_REF] and Lerch et al. [START_REF] Lerch | Forecaster's dilemma: Extreme events and forecast evaluation[END_REF]. A thorough analysis of non local strictly proper scoring rules and the decomposition with quantile and thresholds is given by Ehm et al. [START_REF] Ehm | Of quantiles and expectiles: consistent scoring functions, Choquet representations and forecast rankings[END_REF]. A simple example is given in Section 1.3.3, to provide a better understanding of the evaluation of a probabilistic forecast by the celebrated CRPS.

Binary case

We start from the binary case, where the forecaster wishes to predict the output of a single event and delivers the probability p of occurence of the event according to his opinion. A scoring rule for binary event is local and of the form

(p) = 1 ev S 1 (p) + (1 -1 ev )S 0 (1 -p) ,
where 1 ev = 1 if the event occurs and 0 otherwise. The forecaster suffers S 1 (p) when the event occurs and S 0 (1 -p) when the event does not occur.

The strict propriety imposes that the average score E( (p)) is minimized if (and only if) p is the true probability of occurrence f , and accordingly:

f S 1 (f ) -(1 -f )S 0 (1 -f ) = 0 , (1.1)
thanks to the stationnarity condition d dp (p = f ) = 0. Equivalence between the stationnarity condition and strict propriety is demonstrated by Shuford et al. [START_REF] Shuford | Admissible probability measurement procedures[END_REF].

The scoring rules may be taylored with respect to cost functions, and an infinite number of scoring rules are strictly proper. We restrict the study to the symmetrical case S 0 = S 1 , where we have

f S 1 (f ) = (1 -f )S 1 (1 -f ) , according to Equation 1.1. We define ζ(p) = pS 1 (p), with ζ being symmetrical with respect to 1/2: ζ(p) = ζ(1 -p). The strict propriety necessitates that ζ(p) is non-zero for p = f because for p ∈]0, 1[: d E( (p)) dp = f p ζ(p) - 1 -f 1 -p ζ(1 -p) = ζ(p) p(1 -p) (f -p) .
Losses are usually negatively oriented, meaning that (p) should be as low as possible.

Since the loss should decrease when p gets closer to f , we have ζ(p) < 0. This writing of the loss can also be interpreted with quantile losses. Let α be a decision threshold. We consider the decision-thresholded loss where a false negative costs 1 -α (for 1 ev = 1 and p α), and a false positive costs α (for 1 ev = 0 and α < p). The total loss (p) is a weighted sum of decision-thresholded losses, where the importance of the decision thresholds are given by ω(α):

(p) = 1 ev 1 p ω(α)(1 -α)dα + (1 -1 ev ) p 0 ω(α)αdα . ω(p) (p) (p(1 -p)) -3/2 1 ev 1-p p + (1 -1 ev ) p 1-p (p(1 -p)) -1/2 1 ev arcsin( √ 1 -p) + (1 -1 ev ) arcsin( √ p) -p(1 -p) p(1 -p) 1 ev ( (1-p) 3 3 -(1-p) 4 4 ) + (1 -1 ev )( p 3 3 -p 4 4 )
Table 1.4 -Example of scoring rules in the beta family ω(p) = p a-1 (1 -p) b-1 .

We identify

S 1 (p) = 1 p ω(α)(1 -α)dα and S 0 (1 -p) = p 0 ω(α)αdα. When ω(α) is symmetrical with respect to 1/2 (S 1 = S 0 ), we find ζ(p) = pS 1 (p) = -ω(p)p(1 -p).
We see that the function ζ is closely related to the weighting function ω of the thresholds. For any strictly proper score defined by S 1 , it is possible to define ω(p) = -S 1 (p)/(1-p), and consider the loss as a weighted sum of decision-thresholded losses.

Two classical examples are the quadratic (Brier) score and the logarithm score:

sq (p) = 1 ev (1 -p) 2 + (1 -1 ev )p 2 = (1 ev -p) 2 ln (p) = -(1 ev ln(p) + (1 -1 ev ) ln(1 -p))
respectively obtained for ω(p) = 2 and ω(p) = 1 p(1-p) . Other examples from Buja et al. [START_REF] Buja | Loss functions for binary class probability estimation and classification: Structure and applications[END_REF] are summarized in Table 1.4.

Ranked and continuous case

We now consider events of the type: "y is below the threshold θ k ". With K thresholdtype events, the binary case of Section 1.3.1 easily extends to

(p, y) = K k=1 φ k [H(θ k -y)S 1 (p k ) + H(y -θ k )S 0 (1 -p k )] ,
where p is the list of predicted probabilities p k for the events, the weights φ k > 0 define the importance of the thresholds θ k , and H is the Heaviside step function. We adopt the bold notation p to describe the list or vector with coordinates p k . The events "y is below the threshold θ k " are not a partition of events, consequently this scoring rule is not local. If the dimension of the weights φ k and the dimension of the observation y coincide, the loss (p, y) has the same dimension as y.

In the continous case, for each threshold θ, the weighting φ k is related to φ(θ)dθ. The following continous score is written with the CDF notation:

(G, y) = φ(θ) [H(θ -y)S 1 (G(θ)) + H(y -θ)S 0 (1 -G(θ))] dθ .
Using a uniform weighting scheme over the decision thresholds φ(x) = 1 with the quantile weighting ω(p) = 2 and ω(p) = 1 p(1-p) , we find the CRPS and the Continuous Ranked Ignorance score (CRIGN) defined by:

CRIGN = -H y ln G + (1 -H y ) ln(1 -G) CRPS = H y (1 -G) 2 + (1 -H y )G 2 = (H y -G) 2 ,

Examples with the CRPS

In this section, we provide simple examples with the CRPS, as an evaluation tool. The idea is that for a fixed value of the observation y, several CDFs may reach the same score. We illustrate how these CDFs are related to each other.

In Figure 1.9(a), we show the map of the CRPS for CDF of normal distribution N (µ, σ 2 ) with mean µ and variance σ 2 , for the observation y = 0. For a fixed value of the standard deviation σ, we see that the lowest value of the CRPS is reached for µ = y = 0. Besides, for a fixed value of the mean µ, the lowest value of the CRPS is not reached for the lowest spread σ = 0. Instead, the lowest CRPS value is reached when the spread σ is slightly higher than the error |µ -y| on the mean. The idea behind is that, for an incorrect forecast with an error on the mean, the spread of the distribution lowers the effect of the error |µ -y| > 0.

We compare CRPS values with the case N (1, 1) in Figure 1.9(b). The same CRPS values are obtained for g a = N (0.6022, 0), g b = N (1, 1) and g c = N (0, (2.577) 2 ), with respective CDFs G a , G b and G c . The distribution g a has no spread and a small error of 0.6 on the mean, while the distribution g c has a high spread but no error on the mean. We give an illustration of the CDFs G a , G b and G c , the PDFs g a , g b and g c and the errors (G a -H y ) 2 , (G b -H y ) 2 , (G c -H y ) 2 in Figure 1.10. The strong effect of the square in the CRPS is quite noticeable on the representation of (G b -H y ) 2 , and more generally on the domains where |G(θ) -H(θ -y)| 1. Consequently, the CRPS is not very sensitive to the distribution tails. On this subject, we address in Chapter 4 the question of other scoring rules than the CRPS.

g c g b g a y x 0 1 G c G b G a H y y x 0 1 (G a -H y ) 2 y x 0 1 (G b -H y ) 2 y x 0 1 (G c -H y ) 2

Thesis outline

In Chapter 2, we study multiple ensembles of forecasts of solar radiation, and use sequential aggregation to improve forecast maps of solar radiation.

In Chapter 3, we study the CRPS in the context of model mixture, specifically the bias of the CRPS of an ensemble of forecasts. We introduce sequential aggregation with the CRPS. In Chapter 4, we extend the results of Chapter 3 to other losses than the CRPS admitting a threshold or a quantile decomposition. We also investigate the question of noisy observations.

In Chapter 5 and 6, case studies of probabilistic forecasts of PV power are provided. We show that the methods introduced in Chapter 3 improve the quality of the forecasts for both deterministic and probabilistic evaluation tools. Chapter 7 focuses on PV power forecasts with intraday updates for insular systems such as Réunion and Corsica.

Publications

Introduction

Solar radiation forecasts and especially global horizontal irradiance (GHI) forecasts are needed for the integration of photovoltaic power (PV). The increasing installed capacity of photovoltaic power requires that solar radiation forecasts be always more accurate in terms of spatial and temporal resolutions.

Many meteorological centers provide solar radiation forecasts with two strategies: either a single deterministic forecast, or ensemble forecasts generally with coarser resolution. Deterministic forecasts have been extensively studied for solar and photovoltaic forecasts. Inman et al. [START_REF] Inman | Solar forecasting methods for renewable energy integration[END_REF] and Espinar et al. [START_REF] Espinar | Photovoltaic Forecasting: A state of the art[END_REF] review numerous modeling techniques to generate solar and PV forecasts from meteorological variables. Deterministic predictions from various meteorological centers are compared by Lorenz et al. [START_REF] Lorenz | Benchmarking of different approaches to forecast solar irradiance[END_REF] and Perez et al. [START_REF] Perez | Comparison of numerical weather prediction solar irradiance forecasts in the US, Canada and Europe[END_REF] in a broad range of sites. In the previously cited papers, some forecasting techniques resort to combinations of forecasts. These combinations are derived from a regression over a moving time-window and are applied to a few members.

Vernay et al. [V+13] listed several available maps of solar radiation deriving from cloud products of satellite observations. Even though solar radiation forecasts cover large areas, they are usually compared to ground measurements (PV or solar) or to satellite observations only at measurement sites [e.g., GDM80, among the firsts]. Indeed maps of satellite observations are not broadly used to assess the accuracy of solar radiation forecasts. Morcrette [START_REF] Morcrette | Radiation and cloud radiative properties in the European Centre for Medium Range Weather Forecasts forecasting system[END_REF] used satellite observations as reference to assess the performance of numerical weather predictions (NWP), but not predictions of solar irradiance (e.g., short-wave radiation). Perez et al. [START_REF] Perez | Comparing satellite remote sensing and ground network measurements for the production of site/time specific irradiance data[END_REF] studied the interactions between satellite observations and measurement sites with respect to the distances between sites. Due to the variety of error causes, Thelen and Edwards [START_REF] Thelen | Short-wave radiances: comparison between SEVIRI and the Unified Model[END_REF] restricted the comparison between NWP and satellite observations to reflectance for short-wave radiation. Dehghan et al. [START_REF] Dehghan | Evaluation and improvement of TAPM in estimating solar irradiance in Eastern Australia[END_REF] give emphasis on the spatial resolution of both NWP and satellite data at ground measurement sites.

Ensemble forecasts are classical in meteorology for any field with large uncertainty and for uncertainty quantification. However, no article using ensemble forecast for solar radiation was found in the literature, despite several conference presentations. Still, Yokohata et al. [START_REF] Yokohata | Reliability of multi-model and structurally different single-model ensembles[END_REF] studied climatological ensembles of top atmosphere radiation and radiation in cloud-free conditions.

In the framework of sequential aggregation, a single forecast is built as a linear combination of the ensemble members. The weights of the combination may depend on both time and space. The resulting aggregated forecast is hopefully more skillful than the ensemble members. Cesa-Bianchi and Lugosi [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF] detail the strong mathematical background of these methods, which is summarized and tested by Stoltz [START_REF] Stoltz | Agrégation séquentielle de prédicteurs : méthodologie générale et applications à la prévision de la qualité de l'air et à celle de la consommation électrique[END_REF] and Mallet et al. [START_REF] Mallet | Ozone ensemble forecast with machine learning algorithms[END_REF] and Mallet [START_REF] Mallet | Ensemble forecast of analyses: Coupling data assimilation and sequential aggregation[END_REF] on forecasts of respectively electricity consumption and ozone concentrations.

We propose here to compare the ensemble forecasts of solar radiation from TIGGE (THORPEX Interactive Grand Global Ensemble [START_REF] Bougeault | The THORPEX interactive grand global ensemble[END_REF]) and to combine them. The satellite observations from HelioClim are used in this article as reference observations. The use of both ensemble forecasts (from several sources) and satellite observations makes it possible to generate an aggregated forecast with local combinations on the spatial grid. The Integrated Forecast System (IFS) from ECMWF produces our reference forecast.

In Section 2.2, we describe the TIGGE data sets and we study the performance of the TIGGE ensembles. The HelioClim satellite observations are introduced in Section 2.2.4, along with a detailed comparison with TIGGE forecasts. Our sequential-aggregation strategy is introduced in Section 2.3. It is applied in Section 2.4, where the analysis of the results includes the composition of the ensemble, the spatial patterns in the forecast maps, the forecast time horizon and the sensitivity to the aggregation parameters.

Analysis of TIGGE solar radiation and HelioClim database 2.2.1 Description of TIGGE data

Several meteorological centers provide free-of-charge ensemble forecasts of solar radiation in TIGGE (Table 2.1). The data sets are available on TIGGE with a 2-day delay after the model ran. Detailed studies of the ensemble forecasts from TIGGE were achieved on geopotential height [START_REF] Buizza | A comparison of the ECMWF, MSC, and NCEP global ensemble prediction systems[END_REF], and 850-hPa and 2-m temperatures [START_REF] Hagedorn | Comparing TIGGE multimodel forecasts with reforecastcalibrated ECMWF ensemble forecasts[END_REF] for example. As far as we know, no such study exists for solar radiation.

The temporal resolution of TIGGE forecasts is 6 h with time horizon up to 15 days. Most meteorological centers provide at least two forecast sets per day. For the sake of clarity and because of the 6-h timestep of TIGGE forecasts, we focus on one forecast set per day for each meteorological center. With these constraints, 6 meteorological centers provide an ensemble: China Meteorological Administration (CMA), European Centre for Medium-Range Weather Forecasts (ECMWF), MetOffice (UKMO), Korea Meteorological Administration (KMA), Centro de Previsao Tempo e Estudos Climaticos (CPTEC), and Météo-France (M.-F.). We name "whole ensemble" the ensemble including all 158 members without consideration of their origins, as opposed to the 6 TIGGE center ensembles.

While each center ensemble has a native spatial resolution, our TIGGE data sets are obtained on a common regular grid with the spatial resolution of 0.25 • × 0.25 • , which is finer or close to the native resolutions. The resolution of the TIGGE data sets is coarser than the 0.125 • × 0.125 • resolution of the ECMWF deterministic forecast. Our study focuses on Metropolitan France and the surrounding areas for the 06:00-12:00 UTC accumulation period of day D with model runs starting at 18:00 or 24:00 UTC in day D-1. Thus we study the forecasts for either 12 h or 18 h of lead time. We do not study daily radiation but focus on the shortest timestep available in TIGGE. The flux values from TIGGE are averaged over the 6-h timestep so that the values of the forecasts for 12:00 UTC are expressed in W m -2 and correspond to the averaged flux between 06:00 and 12:00 UTC.

The nature of the solar radiation data sets in TIGGE is mostly net shortwave solar radiation (SSR) as defined for classical meteorological fields. Net shortwave solar radiation is the fraction of the downwards shortwave solar radiation (SSRD) absorbed by the ground on an horizontal plane. TIGGE data sets and are SSRD data. In the context of photovoltaic production, we are interested in SSRD. The well-known global horizontal irradiance (GHI) refers to SSRD.

The albedo coefficient α is the reflection coefficient of the ground. The albedo defines a relationship between SSR and SSRD, where the incident flux is divided between the absorbed flux and the reflected flux. Thus we deduce that:

SSRD = SSR 1 -α . (2.1)
Depending on the ground surface, the albedo coefficient can vary in space and time.

Analysis of the TIGGE ensembles of forecasts

Now we analyze and compare the ensembles of forecasts over the area spanning 41 • to 51.50 • in latitude and -5.50 • to 10 • in longitude. Two data sets are used in our study. The main data set consists of 350 consecutive days starting on 2012-01-02. The secondary data set is dated from 2011-06-09 to 2011-09-05 due to data availability and includes only 100 random locations. The secondary data set is only used as learning data set. An example of the annual average of the 6-hour forecasts at 12:00 UTC for one center (KMA) is provided in Figure 2.1, in order to show the large spatial variability of the average forecast.

We propose two sorting procedures so as to consistently number the forecasts in time. This step is required for the weights in our aggregation to be clearly associated with a given ensemble member and therefore to be meaningful. The two sorting procedures rely on the rank at each grid point. At any time and location, the first member of a sorted ensemble always provides the lowest value, and the last member gives the highest value. The first sorting procedure is applied to the 158 members of the whole ensemble, and the second sorting procedure is applied separately in each center ensemble.

Correlation coefficients between members are computed in order to quantify the similarities between the members (Figure 2.2). In the correlation matrices, the rectangles and squares materialize some separation between the different ensembles. Indeed the correlation matrices reveal that the members are especially close to one another within the same ensemble. In other words, one member from a given meteorological center ensemble is more correlated to another member of the same ensemble than to another 240 270 300 330 360 390 420 450

No sorting

Sorting by center Full sorting . Each row and column of the correlation matrices is dedicated to one given member. The rows and columns appear in the same order. Each matrix entry is the correlation coefficient between the merged data (all timesteps, 100 random locations) of a pair of members. On the left, the correlations are shown between raw members, as they are retrieved from TIGGE. In the middle, the correlations are computed after sorting within each center ensemble. On the right, the correlations are obtained after a full sorting of the forecasts.

C M A E C M W F U K M O K M A C P T E C M e t e o F r a n c e C M A E C M W F U K M O K M A C P T E C M e t e
member from a different ensemble. The CPTEC ensemble is very distinguishable due to the extremely high correlations between its own members and also due to the low resemblance between its members and the others.

Sorting has generally two effects on the ensembles. First, the sorted members with close ranks have higher correlations among them than non-sorted members. Second, the pairs of members with the lowest correlation coefficient are found between sorted members of extreme rank.

Reference performance measures

The strengths and weaknesses of the statistical indicators commonly used in solar forecasting are developed in Hoff et al. [START_REF] Hoff | Reporting of irradiance modeling relative prediction errors[END_REF]. The well-known root mean square error (RMSE) and mean absolute error (MAE) are classical performance indicators. The RMSE of the predictions y with respect to the observations y over the set S is given by

RM SE = 1 |S| s∈S ( y s -y s ) 2 , (2.2)
where |S| is the number of elements (cardinality) in S, and s indexes space or time or both. In case s describes all locations at one single timestep, the spatial RMSE is computed. The MAE is calculated in a similar way as

M AE = 1 |S| s∈S | y s -y s | . (2.
3)

It is noticeable that the errors y s -y s are computed independently for each s, hence errors due to geographical or temporal shifts are not detected as such. In case of missing predictions, the missing indices are excluded from the set S.

The average observed value is introduced to define relative indicators; e.g. for the relative RMSE:

rRM SE = RM SE 1 |S| s∈S y s , ( 2.4)
where rRMSE and RMSE are computed over the same set S. Except for the scores of Table 2.2, HelioClim is used as the observation y s . Because our main interest is to provide information over land and not over sea, all scores and regression coefficients are given for land locations only.

Comparison with HelioClim

Description of HelioClim data

Provided by Transvalor (Armines), HelioClim satellite observations are based on MSG satellite data (Meteosat Second Generation). The operational version HC3-v4 is used here. The HelioSat2 method [START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF] has been deployed to generate this He-lioClim database [START_REF] Blanc | The HelioClim Project: Surface Solar Irradiance Data for Climate Applications[END_REF]. The satellite estimations rely on instantaneous reflectance measurements. Data are acquired every 15 min, converted to an estimation of the incident radiation flux, and averaged over the 15-min timestep. The spatial resolution of HelioClim over France is natively between 3 and 5 kilometers; our HelioClim data were retrieved with a spatial resolution of 1/12 • , which is already much finer that the resolution of the forecasts. Our zone of interest includes five BSRN stations (Baseline Surface Radiation Network) [START_REF] Ohmura | Baseline Surface Radiation Network (BSRN/WCRP): New precision radiometry for climate research[END_REF] for the evaluation of HelioClim performance. The stations are located in Camborne (United Kingdom), Cener (Spain), Carpentras (South France), Palaiseau (North France), and Payerne (Switzerland), as exhibited in Figure 2.1. The evaluation results (Table 2.2) are computed over several years and show an average relative RMSE of 10.7% and an average relative MAE of 8%.

In Section 2.4, we wish to produce our own forecasts based on the previously described ensembles and we wish our forecasts to be at the finest available resolution, which is the spatial resolution of our HelioClim data. Consequently our study is carried out at the resolution of HelioClim. All the forecasts are interpolated by bilinear interpolation to reach the same resolution of 1/12 • , for a total amount of 127 × 187 grid points. On the high-resolution grid, the forecasts vary spatially slowly compared to the satellite observations.

TIGGE ensembles and HelioClim

The difference of nature between the HelioClim incident radiation (SSRD) and the ground-absorbed radiation from TIGGE (SSR) prevents a direct comparison. Therefore we test several conversion methods in Appendix 2.A. The method of linear conversion (lin) based on historical data shows the best RMSE. This method is used until the end of the section.

The monthly RMSEs are impacted by the seasonal variability (Figure 2.3). Solar radiation forecasting is more difficult between April and July, during the brightest days with large variability. It is worthy of notice that the ranking of the ensemble means is steady over time. The ensemble means of KMA and UKMO (after linear conversion) show better scores than the reference forecast, while their spatial resolution is natively poorer. An improved forecast, called TIGGE-mean, may be built as the mean of the TIGGE ensemble means (CPTEC excluded) with linear conversion. The score of the TIGGE-mean lies at 66.9 W m -2 , which is better than any TIGGE ensemble mean. The same analysis was undertaken with the MAE (not shown) and reveals similar trends: 54.4 W m -2 (reference forecast MAE), 50.2 W m -2 (UKMO ensemble mean with linear conversion MAE), 48.0 W m -2 (TIGGE-mean MAE).

We highlight the fact that no center ensemble is steadily the closest to the satellite observations, whatever the score discrepancy described above. If we track over the consecutive timesteps the origin of the best member at each location, we find that this origin changes from one timestep to the next with a frequency of 75%. Furthermore, the proportion of times and locations where the best member belongs to a given center is reported here: reference ECMWF 3%, CMA 10%, ensemble ECMWF 26%, UKMO 19%, KMA 18%, CPTEC 2%, Météo-France 23%. Considering the fact that the ECMWF reference forecast is one single member compared to 158 members, its frequency of being the closest member to the observation is rather high. On the opposite, the CPTEC ensemble and, to a lesser extent, the CMA ensemble show the lowest frequencies.

Rank histograms [And96; TVS99; HC97] evaluate the quality of the spread of an ensemble. Each observation is given a rank, which corresponds to the number of members with lower values than the considered observation. Then the distribution of rank frequencies can reveal the presence of under-dispersion (U-shaped histogram), overdispersion, and biases. The rank histograms of the center ensembles and the rank histogram of the whole 158-member ensemble are shown in Figure 2.4 and 2.5, for data converted by center (lin). The values of the bars of each histogram are normalized with respect to the total number of observations, so that the height of the bars of the ideally flat histogram is always 1. There are clear outliers on all histograms, which is a clear sign of general under-dispersion [START_REF] Hamill | Interpretation of rank histograms for verifying ensemble forecasts[END_REF]. The frequencies with which observations fall within an ensemble envelop are: CMA 34%, ECMWF 45%, UKMO 58%, KMA 56%, CPTEC 1%, Météo-France 59%, whole ensemble 89%. .5 -Rank histograms for the whole ensemble. In the ideal case of a flat histogram, the height of all the bars would be equal to 1. The gray scale indexes the number of ensembles whose spread contains the observation, from black (no ensemble) to white (all ensembles, but not observed because of the very low dispersion of CPTEC).

The whole ensemble is less under-dispersed than each center ensemble and exhibits a rank structure with peaks. The ranks of the major inner peaks (14, 28, 37, 78, 112, 121, 135 and 144) correspond to combinations of ensemble sizes, starting from either extreme peak. For example, the peak at rank 28 corresponds to the combination of the ensemble sizes of CPTEC and CMA (both 14). The peak at rank 135 corresponds to an ensemble of 23 members (UKMO or KMA) starting from the outer peak at rank 158. The peaks are generated by the 4% share of the observations (inner black bars) which do not fall within any of the center ensemble envelops. Indeed these observations are necessarily indexed at ranks related to combinations of ensemble sizes. We conclude that the peaks are due to the variety of the under-dispersed ensembles.

Ensemble forecast strategy: sequential aggregation

In this section we detail the principles of sequential aggregation. In particular, we explain the method of discounted ridge regression. The performance and robustness of discounted ridge regression have been tested for the case of air quality [START_REF] Mallet | Ozone ensemble forecast with machine learning algorithms[END_REF]. While time series of scalar fields were considered above, in this section we only consider scalar time series.

Notation

Let x m,t describe the m-th member of our forecast ensemble at time t, with m ∈ {1, . . . , M } indexing the members and t ∈ {1, . . . , T } indexing the forecast timesteps. The vector x t refers to the ensemble of forecasts [x 1,t , x 2,t , . . . , x M,t ] . At each timestep, the members should be conveniently combined with the weights w m,t to generate the forecast

y t = M m=1 w m,t x m,t (2.5)
of the observation y t .

Sequential aggregation: method

The aggregation weights w m,t are updated before the forecast step t, using only past observations y 1 , y 2 , . . . , y t-1 and past simulations x 1 , x 2 and x t-1 . For the discounted ridge regression with parameters (λ, γ), the weight vector

w t = [w 1,t , w 2,t , ..., w M,t ] is found through minimization of J(u) = λ u -w ref 2 2 + t-1 t =1 β γ (t -t ) × (y t -u • x t ) 2 (2.6) with β γ (t -t ) = 1 + γ (t -t ) 2 (2.7)
and w ref a reference weight vector chosen by the operator and constant in time. The parameter λ affects the distance between w and the reference vector w ref (usually set to zero or to [1/M, . . . , 1/M ] , following the ensemble mean). The function β γ gives higher importance to the most recent timesteps. When both λ and γ are set to zero, a simple recursive least-square regression is achieved.

The classical ridge regression without discount provides the theoretical guarantee that the final score of the aggregated forecast will be close to the final score of the best constant linear combination. Indeed we have

T t=1 1 T (y t -y t ) 2 -min u∈B M T t=1 1 T (y t -u • x t ) 2 O ln T T , (2.8)
under the assumption of bounded losses (y t -y t ) 2 . The condition u ∈ B M , where B M is a 2-norm ball in R M , means that u 2 2 is bounded. The best linear combination with constant weights (implicitly defined above by the second term in the left hand side) is named the oracle. The oracle is found by least-square regression over the whole set. By definition, the oracle shows better performance than any member in the ensemble. Consequently the aggregated forecast is more skillful in the long run than the best member. The discounted ridge regression provides asymptotically this guarantee, which is verified for each sequence of discounted regret.

Algorithm

Parameters: λ, γ, w ref ;

Initialization: w 1 ;
For each time index t = 1, 2, ..., T 1. get the predictions x t , 2. compute y t , with x t and w t , 3. get the observation y t and compute w t+1 .

The initial weight vector w 1 is arbitrarily set, e.g., to [1/M, . . . , 1/M ] .

Application

Experiment setup

The sequential aggregation with discounted ridge regression as described in 2.3.2 is applied independently at each location of the 127 × 187 grid. In a similar fashion to the study of Section 2.2, the forecast variable is the incident radiation flux integrated between 06:00 and 12:00 UTC. The TIGGE forecasts are available at 18:00 or 24:00 on day D-1 to forecast the quantity of interest for day D at 12:00, also named (D, 12:00).

The ensemble data are SSR data from TIGGE without any SSR-SSRD conversion because the aggregation does not depend on any multiplicative coefficient applied to the members. In other words, the weights of the aggregation produce multiplicative corrections and solve the issue of the nature of TIGGE data.

The aggregation parameters (λ, γ) are respectively set to 6 × 10 6 and 20 by default (see Section 2.4.2 for the assessment of the parameters). The values of the reference vector w ref are set to 1/M in order to drive the aggregated forecast towards the ensemble mean. Even though the ensemble mean may not be the most appropriate reference vector, the vector w ref mostly impacts the beginning of the aggregation so that the critical parameters are truly (λ, γ).

The aggregation may be achieved in a single step by choosing members from all of the center ensembles. Another approach involves two steps: a first aggregation within each center ensemble and a second aggregation with the resulting forecasts as members. Both procedures have been tested and only the procedure achieved in a single step is presented here since the second procedure did not lead to significantly different results.

In order to study the impact of the number of members M , the same amount of members are chosen from each meteorological center, but the center ensembles do not have the same size. Therefore, the members are chosen in a way that their ranks are linearly spaced and centered on the median of each center ensemble. The full sorting procedure is not impacted by this member selection because in this case the members are mixed up before sorting and rank selection. It is possible to realize aggregation with one single member, such as the reference forecast. In this case the weight plays the role of a local correction factor.

Missing data (only CPTEC ensemble) are replaced by the ensemble mean of the available members.

Results

Aggregation example with one center ensemble

In this section, the aggregation is first run with the ensemble KMA only (without sorting) and then with an additional deterministic member. In the first case, we indicate that the RMSE of the SSRD ensemble mean is 75.0 W m -2 , whereas the RMSE of the ensemble mean with linear conversion is 72.1 W m -2 . The RMSE of the aggregated forecast equals 70.0 W m -2 and may decrease when data from another source is included in the ensemble. For example, when the ECMWF deterministic forecast is added as a member, the RMSE reaches 67.5 W m -2 . The scores are to be compared to the 58.9 W m -2 RMSE of the oracle of the same ensemble and to the 65.1 W m -2 RMSE of the oracle with only KMA data.

One question is whether the improvement due to the deterministic member originates from its high spatial resolution. We therefore include the ECMWF deterministic forecast at the lower resolution of 0.25 • × 0.25 • which is obtained by averaging the fine 0.125 • × 0.125 • forecast. In that case, the RMSE is also equal to 67.5 W m -2 , which means that the resolution of the reference forecast is not the key factor. For comparison, the RMSE of the reference forecast corrected by discounted ridge regression is equal to 68.0 W m -2 .

Oracles with orthogonal members

We now want to quantify the potential improvements brought by sequential aggregation with all members, using a relevant oracle. The problem of overfitting can arise because of the large number of members ( 158 time sequence. In that case, the score of the oracle is artificially good, and the competition against the oracle does not carry any meaning. We present here oracles computed with the help of principal component analysis (PCA), so as to avoid overfitted oracles. The PCA generates orthogonal modes, and consequently orthogonal members. The 158 members are sorted within each center ensemble and then orthogonalized by PCA, independently at each grid point. We compute the oracle of the ensemble of size M , by selecting M orthogonal members based on the M modes explaining the largest variance share of the ensemble. In Figure 2.6 the RMSE of the oracle and the total amount of variance explained by its members are plotted against the number of orthogonal members. An indication of possible overfitting is shown in Figure 2.6, because the RMSE still decreases with the number of PCA modes while the share of unexplained variance is small. Indeed the 32 first orthogonal members explain 95% of the variance and generate an oracle with an RMSE of 54.4 W m -2 , whereas the RMSE of the oracle with all PCA modes (100% of explained variance) equals 33.9 W m -2 . Consequently, the 32-member oracle is considered as the relevant oracle; its score is a relevant evaluation of the best score that may be achieved by linear combination without overfitting.

Typical aggregation

In this paper, the "typical aggregation" refers to the aggregation with 30 members (5 per sorted center ensemble), with default parameters λ = 6 × 10 6 and γ = 20. The resulting forecast is simply referred to as typical aggregated forecast and is used below to illustrate spatial and temporal features resulting from sequential aggregation. Note that the typical aggregation is compared to all members in Figure 2.13 and to the ensemble means for each month in Figure 2.3. We highlight the fact that the typical aggregated forecast performs better than any of the individual members, steadily over time.

The maps of scores of RMSE, MAE and relative RMSE are shown in mountains severely degrades the forecast accuracy. Indeed the absolute largest errors are found in the region of the Alps and in the region of the Pyrenees. The coarse resolution of the members may explain the poor performance in the regions of complex terrain. Also, in the mountainous regions, the albedo shows higher temporal variations, which is difficult to catch for the aggregation. On the opposite, the best scores are steadily achieved in the inner lands of Spain and in the south-east of France. High relative errors are numerous in the northern area especially during the spring-summer period, because of the low values of the observations combined with large errors. For example in spring and summer period, the British area shows a relative RMSE higher than 20%, whereas the south-east of France reaches a relative RMSE below 12.5%, even though the two areas are associated with similar RMSEs over the same period.

We recall here that the relative RMSE of the satellite observations compared to BSRN stations is worth 10.7% on average.

Most of the members are computed with a low spatial resolution and do not show fine long-term spatial structures. On the contrary, the aggregated forecast is built independently at each location, which allows the procedure to adapt locally and to finally show fine structures, that are resolved by HelioClim. These structures are finer than any structures found in the ensemble members (due to their low resolutions), even in the ECMWF deterministic forecasts (Figure 2.8).

In order to compare predictions performance at each date, the spatial RMSE of various predictions are temporally sorted according to the performance of the ECMWF reference forecast (Figure 2.9). The overall trend is the same for the three forecasts. The aggregated forecast often shows the best performance, but it does not always perform better than the aggregation applied to the reference forecast only. One benefit of aggregation is that large improvements are achieved at the most difficult timesteps, even with one single member in the ensemble.

Members selection and sorting

The performance of the aggregated forecast is impacted by the number of members up to a certain extent (Figure 2.10). Beyond a total amount of roughly 60 members (10 members from each center ensemble), the scores reach a plateau. The full sorting procedure is the first to reach its own plateau, which may be caused by the merger process of the members with sorting. Besides, the sorting procedures impact the scores with limited consequences (less than 1 W m -2 ) for large enough ensembles. According to the observed performance, the sorting procedure by center should be preferred if more than 5 members per center are chosen. It is noticeable that the reference forecast brings a score improvement whatever the number of aggregated members. The best score achieved with default (λ, γ) equals 61.5 W m -2 for the RMSE and 43.6 W m -2 for the MAE. Compared to the scores of reference forecast, the improvements brought by the aggregation equal 21.2% for the RMSE and 19.8% for the MAE. Besides, compared to the scores of TIGGE-mean, these improvements are worth 8.1% for the RMSE and 9.1% for the MAE.

The aggregation ensemble may also be built with only 5 out of 6 center ensembles. We compare the aggregated forecast of 30 members, so that 6 members are chosen in each of 5 center ensembles. In this case, the reference aggregated forecast is the typical aggregated forecast with 5 members chosen in each of the 6 center ensembles. Either way, the members are sorted per ensemble. We found that the aggregated forecast with 5 centers performs better than the typical aggregated forecast, when the omitted members are from CMA or CPTEC. The largest improvement occurs when the CPTEC is excluded with a benefit of only 0.3 W m -2 . On the opposite, the score is the most severely degraded when the members from KMA are left out, generating the RMSE of 64.1 W m -2 , while the RMSE of the typical aggregated forecast equals 62.5 W m -2 . To conclude, not using all center ensembles to build the aggregation ensemble can generate both benefits and loss in terms of score. However, these benefits are lower than the loss caused by leaving out skillful forecasts. In practice, one may arguably keep all center ensembles.

Time horizon

We consider an ensemble forecast delivered at time t for the k next steps. The weights of the sequential aggregation are commonly used for the first step ahead (see the algorithm, in Section 2.3.3), but can also be used for the following timesteps t + k, where k > 1. In fact, the weights can even be used for timesteps that are not included in the sequential aggregation, such as the forecast for 18:00 UTC. This new framework allows us to generate the aggregated forecast for (D, 18:00), (D+1, 12:00) and (D+1, 18:00) with the weights initially computed for the prediction (D, 12:00). In terms of time horizon, the prediction for (D, 12:00) corresponds to the 12 h horizon, (D, 18:00) to 18 h, (D+1, 12:00) to 36 h and (D+1, 18:00) to 42 h. The time horizons of Météo-France forecasts are actually 6 h longer, because its predictions start at 18:00 and not 24:00. The night time predictions steps are skipped because of their very low values.

The RMSE of the ensemble means and the scores the typical aggregated forecast are shown in Figure 2.11, depending on the time horizon. The scores of CPTEC ensemble mean are not shown, because of their high values. The six presented forecasts share the same trends. The scores of predictions for 18:00 beat the scores for predictions 12:00 of the same day, by less than 4 W m -2 on average, which is consistent with the fact that the forecasts for 12:00 have higher values than the forecasts for 18:00. Furthermore, the predictions for day D+1 show an average degradation of more than 5 W m -2 compared to the prediction for day D. Additionally, we notice that the ensemble mean from ECMWF shows the steadiest performance over time and that the benefits of the aggregated forecast are retained for longer time horizons. Figure 2.12 -RMSE grid (W m -2 ) depending on the aggregation parameters, aggregating a 134-member ensemble using at most 30 sorted members from each center ensemble.

Influence of the aggregation parameters

We searched the best set of parameters (λ, γ) that provides the lowest RMSE for the data set of year 2011, already presented in Section 2.2.2. We found that at most 3 W m -2 were to be gained in terms of RMSE on a wide range of parameters (Figure 2.12(a)). Our default parameters (6.10 6 , 20) produce the score of 67.1 W m -2 , which is close to the best score achieved on the grid of Figure 2.12(a) (e. g. 66.7 W m -2 with λ = 6.10 6 and γ = 5). Furthermore, we found that the best set of parameters varies in terms of space and time. Indeed, if we exclude the 50 first timesteps in the RMSE, then the most appropriate parameters are λ = 10 7 and γ = 20. If we choose a posteriori the optimal (λ, γ) for each grid point, only 1 W m -2 is to be gained from default parameters. The minor variations of performance guarantee that we may choose our default parameters within one order of magnitude, and test them on another data set without significant loss of performance.

We produce the same analysis for the data set of year 2012 in Figure 2.12(b). In this case, the best parameters are λ = 5.10 6 and γ = 10, with the score of 61.5 W m -2 that is very similar to the score of 61.7 W m -2 obtained with default parameters. Therefore, the order of magnitude of the parameters for year 2012 can be deduced from the data set of year 2011. Besides we also found that if we choose the best (λ, γ) for each grid point, the gain does not exceed 1 W m -2 . As stated with the data set of year 2011, once the relevant order of magnitude for the parameters (λ, γ) is known, only little improvement is possible by adjusting them.

Conclusion

The ensemble forecasts from TIGGE provide a wide range of meteorological fields including net short-wave solar radiation, with the large timestep of 6 h. After conversion based on a constant albedo, the resulting ensembles are under-dispersed, even when grouped together. The performance of the ensemble means is assessed with RMSE and MAE and compared to the deterministic forecast from ECMWF. The reference forecast has higher resolution than the best TIGGE ensemble means, but produces similar scores. Sequential aggregation brings improvements to the TIGGE ensembles on several features, including RMSE and MAE, with theoretical guarantees. In this study, the aggregated forecast performs better than any member of the ensembles and any ensemble mean. On average, the aggregated forecast retrieves all spatial patterns, even at a much finer resolution than any of the members. Besides, the members combination proves to be consistent with the time horizon. Finally sequential aggregation is easyto-use for its parameters do not need accurate values.

Practical applications of the aggregation algorithm should investigate higher temporal resolutions, especially the hourly timestep. Next developments may focus on the study of uncertainty with sequential aggregation, which is possible using filtering [START_REF] Mallet | Minimax filtering for sequential aggregation: Application to ensemble forecast of ozone analyses[END_REF], but without the same theoretical robustness as in this paper. The introduction of multiple model runs per day with temporal interpolations may result in a robust framework for intraday forecasting. Furthermore, new procedures are to be developed to explore sequential aggregation with spatial dependencies between grid points, with the objective of better forecasting daily spatial patterns. There is also a need for quantifying the possible resolution improvements brought by the aggregation.

Appendix 2.A Conversion from SSR to SSRD and reference forecast 2.A.1 Methods

Several conversion methods allowing the estimation of the albedo are tested. Even though the KMA data set is already SSRD and does not need data conversion, the below methods are also tested on the KMA data set, for the sake of completeness.

First, we infer three conversions from our reference forecast. The SSR and SSRD forecasts from ECMWF are used to provide three SSR-SSRD empirical conversions (Table 2.3): with a constant 1.18 coefficient (glob), with local multiplicative coefficients (mult), and with local additive coefficients (add). The local multiplicative conversion is in fact the local estimation (in space and time) of the factor 1/(1 -α) based on ECMWF continuous forecast.

Although it may lead to local large errors due to seasonal changes in snow cover and vegetation, we model a constant albedo in space and time (glob). The linear relationship (Eq. (2.1)) between SSR and SSRD forecasts from ECMWF is found by linear regression on data from 100 random locations and for 350 days in year 2012. The resulting slope value (supposed to be equal to 1/(1 -α)) is 1.18 and the intercept value equals 13 W m -2 . The squared correlation coefficient R 2 SSR and SSRD forecasts are strongly correlated in practice. According to the ratio between ECMWF SSRD and SSR forecasts, the relative standard deviation (standard deviation divided by mean) of the coefficient 1/(1 -α) is equal to 8.8% on average over all the grid points. Consequently the local variations of the ratio 1/(1 -α) are rather small compared to its local mean value. Besides, more than 90% of the values of the coefficient 1/(1 -α) are comprised between 1.13 and 1.48.

Second, a constant albedo is computed for each center ensemble. This method is referred to as "lin" in Table 2.3 and is also referred to as the linear conversion below, because the method is based on linear regressions. Taking HelioClim observations as SSRD data and TIGGE ensemble means as SSR data, linear regressions are carried out on past data sets of year 2011 described in Section 2.2.2. The results of the regressions (Table 2.4) show two slopes lower than one (KMA and CPTEC), and intercepts ranging from -21 W m -2 to 41 W m -2 . Remind that KMA data is already SSRD data so that its slope is not related to the albedo. The diversity of slopes and intercepts values suggest that the albedo coefficient should be evaluated independently for each center ensemble. The abscissa is the member label and the members are grouped by origin. The dashed line is the score of the reference forecast. The dashed-dotted line is the score of a typical aggregated forecast (see Section 2.4.2).

2.A.2 Numerical results

The deterministic forecasts from ECMWF provide several scores. We build another deterministic forecast, called SSRD-lin, based on deterministic SSR and the linear conversion method. Compared to HelioClim, the RMSE of the ECMWF SSRD deterministic forecast equals 78.0 W m -2 while the RMSE of SSRD-lin equals 76.9 W m -2 . When the largest error percentiles of SSRD and SSRD-lin forecasts are compared, we see that the largest errors belong to SSRD-lin only beyond the 99 th percentile. On the one hand, the physical approach with local albedos (as opposed to SSRD-lin) does not generate the largest errors. On the other hand, local albedos do not provide the best score.

Furthermore, the resolution of the ECMWF deterministic forecast does not impact its performance. Indeed the ECMWF deterministic forecast at degraded resolution of 0.25 • × 0.25 • , by interpolation, shows the RMSE of 77.7 W m -2 when compared to HelioClim, which is below the RMSE of the reference forecast. The two predictions at different resolutions have similar performance, since the distance between them, measured with a root mean square discrepancy (in time and space), equals 3.0 W m -2 .

The RMSEs are shown in Figure 2.13 for all members without sorting procedures. The RMSEs between individual members and satellite observations are at least twice as big as the RMSEs between satellite observations and ground measurements (Table 2.2). It is noteworthy that members of a same center ensemble perform similarly. The scores strongly depend on the conversion methods and on the origin of each forecast. We highlight again the fact that the nature of KMA data is already SSRD.

The linear conversion (lin) is the best conversion method tested here except for the CPTEC ensemble. The multiplicative conversion method was supposed to provide accurate forecasts. However, here is confirmed the need for the conversion method to adapt to each center ensemble.

3 Online learning with the CRPS for ensemble forecasting

Introduction

The minimization of the CRPS is a common way to drive probabilistic forecasts [Gne+05; JMA15]. From diagnostic tools to modeling techniques, Gneiting and Katzfuss [START_REF] Gneiting | Probabilistic forecasting[END_REF] review the state of the art of probabilistic forecasting. Using several forecasts, based on various models and perturbed input data, is a common way to produce probabilistic forecasts [START_REF] Leutbecher | Ensemble forecasting[END_REF]. The roots of this framework known as ensemble forecasting is reviewed by Lewis [START_REF] Lewis | Roots of ensemble forecasting[END_REF]. Ensemble of forecasts is the raw material of the techniques proposed in this paper.

Sequential aggregation targets optimal combinations, as thoroughly introduced in the monograph Cesa-Bianchi and Lugosi [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF]. These techniques, also known under the scope of online learning, come with attractive theoretical guarantees of performance. Stoltz [START_REF] Stoltz | Agrégation séquentielle de prédicteurs : méthodologie générale et applications à la prévision de la qualité de l'air et à celle de la consommation électrique[END_REF] and Mallet et al. [START_REF] Mallet | Ozone ensemble forecast with machine learning algorithms[END_REF] and Mallet [START_REF] Mallet | Ensemble forecast of analyses: Coupling data assimilation and sequential aggregation[END_REF] summarized and tested these techniques on forecasts of respectively electricity consumption and ozone concentrations. Usually focused on scalar forecasting, sequential aggregation was applied to the Brier score and the quantile score by respectively Vovk and Zhdanov [START_REF] Vovk | Prediction With Expert Advice For The Brier Game[END_REF] and Biau and Patra [START_REF] Biau | Sequential quantile prediction of time series[END_REF]. In this paper, we use sequential aggregation in order to target the best CRPS, with theoretical guarantees that require essentially no assumption on the forecast or observation distributions. In this sense, our method is a non-parametric post-processing method. Our techniques generate weights for each ensemble member so as to produce a linear opinion pool, also known as model mixture [GM90; CW99]. Ranjan and Gneiting [START_REF] Ranjan | Combining probability forecasts[END_REF] provide mathematical grounds on these combinations. Our technique was first designed to work with an ensemble of scalar forecasts. Still, it can be applied when a parameterized distribution is associated to each forecast.

In Section 3.1, we describe the mathematical background on the CRPS. We provide contributions related to ensemble forecasting with discrete Cumulative Distribution Functions (CDFs). Our contributions are mainly generalizations of existing results to the case of combinations of forecasts with unequal weights, in a probabilistic framework. We also provide a framework to work with classes of members, compatible with fair probabilistic evaluations. In Section 3.2, we detail online learning techniques, with adaptation for probabilistic ensemble forecasting based on the CRPS. In Section 3.3, we illustrate the notions of Section 3.1 with numerical experiments, and we demonstrate our algorithms with numerical examples. We summarize several useful identities involving CDFs in Appendix 3.A.

Mathematical background

Bibliographical remarks

The evaluation of probabilistic forecasts is a long range discussion going on since Winkler and Murphy [START_REF] Winkler | Good" probability assessors[END_REF] and Savage [START_REF] Savage | Elicitation of personal probabilities and expectations[END_REF]; see Dawid [START_REF] Dawid | Comments on: Assessing probabilistic forecasts of multivariate quantities, with an application to ensemble predictions of surface winds[END_REF] for a detailed bibliographical analysis, and more recently Gneiting and Raftery [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF] and Candille and Talagrand [START_REF] Candille | Evaluation of probabilistic prediction systems for a scalar variable[END_REF] for detailed analyses. The Brier score was introduced by Brier [START_REF] Brier | Verification of Forecasts Expressed in Terms of Probability[END_REF] and Good [START_REF] Good | Rational decisions[END_REF] to evaluate probabilistic forecasts for a given threshold and a binary observation. The Continuous Ranked Probability Score (CRPS) can be viewed as a continuous version of the Brier score [Eps69; Mur71] for any threshold.

The Continuous Ranked Probability Score (CRPS)

We want to forecast a scalar quantity y called the verification and we suppose that y admits an underlying distribution that is described by the CDF F . The CRPS is considered as a realization of a random variable, and it is defined as

CRPS(G, y) = (G -H y ) 2 , (3.1)
where G is a CDF that is chosen by the forecaster in order to predict F , H is the unit (or Heaviside) step function, and H y (x) indicates a centered Heaviside function H(x -y). The CRPS is negatively oriented, meaning that lower scores imply better performance. Gneiting and Raftery [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF] show that the CRPS may also be written as

CRPS(G, y) = E(|X -y|) - 1 2 E(|X -X |) , (3.2)
where E is the expectation, and both X and X are two random variables drawn from G. A decomposition of the average CRPS was introduced by Hersbach [START_REF] Hersbach | Decomposition of the continuous ranked probability score for ensemble prediction systems[END_REF]. The decomposition of scores into divergence and uncertainty terms is explained in Bröcker [START_REF] Bröcker | Reliability, sufficiency, and the decomposition of proper scores[END_REF]. The average CRPS is decomposed as follows:

CRPS(G, y)dF(y) = (G -F) 2 + F(1 -F) , (3.3)
where y is integrated over the values defined by F (using Equation 3.28 of Appendix 3.A).

The CRPS is a strictly proper score, which means that it is minimized on average if and only if the forecaster's choice G is equal to F. This is a straightforward observation from Equation 3.3. The strict propriety of the CRPS can be compared to the non-strict propriety of the square loss [START_REF] Bröcker | Scoring probabilistic forecasts: The importance of being proper[END_REF], which reads

(y -E(X)) 2 = G -H y 2 , ( 3.4)
according to Equation 3.31 of Appendix 3.A. We see that minimizing the square loss and minimizing the CRPS (Equation 3.1) are rather different objectives, due to the location of the square function inside or outside the integral expression. The CRPS objective is more demanding, because in this case the integration is applied to a positive function.

The ensemble CRPS

In the case of ensemble forecasting, the forecaster relies on an ensemble of M members x m , m ∈ {1, . . . , M }, to construct a CDF. The empirical CDF G E of the ensemble is a step function with jumps of heights u m (called weights) at the members values x m . Thus we write G E (x) = M m=1 u m H(x -x m ). In order to satisfy the definition of a CDF, the weights u m should be nonnegative and sum to one, so that they produce a convex combination. Such weight vectors define the simplex P M of R M . The weights u m are to be optimized in order to minimize the CRPS.

The computation of the integral of Equation 3.1 is easy on step functions G E . When the CDF is a step function, we refer to the score as the ensemble CRPS:

CRPS(G E , y) = M m=1 u m |x m -y| - 1 2 M m,k=1 u m u k |x m -x k | . (3.5)
The derivation of Equation 3.1 is detailed in Appendix 3.B.

Without further information, the members are assumed to be i.i.d., thus the forecaster may arguably choose all weights equal to 1/M . By definition, a scoring rule depending on the verification y and i.i.d. members x m is fair if the average score is minimized when the members and the verification are sampled from the same distribution. Ferro et al. [START_REF] Ferro | On the effect of ensemble size on the discrete and continuous ranked probability scores[END_REF] show that the ensemble CRPS is unfair due to the finite size of the ensemble. In the next section, we generalize this result to the case of unequal weights, with non identically distributed members.

The bias of the score is an important topic in our optimization framework. Indeed if our objective function is intrinsically biased, then the resulting probabilistic forecast cannot be calibrated.

Bias of the ensemble CRPS with underlying mixture model

We consider that the members x m are independent samples from the CDFs G m , and that y is fixed. The purpose of this section is to compare the score obtained with the step function G E averaged according to the CDFs G m and the score obtained with the mixture model described by the average CDF G = u m G m .

Taking the expectation with respect to the members x m , we show that

E(CRPS(G E , y)) = H y -2 M m=1 u m G m H y + M m =k u m u k G m G k + M m=1 u 2 m G m , (3.6)
using Equation 3.28. The trick is that H 2 (x-x m ) = H(x-x m ), thus the average CRPS does not include G 2 m terms but G m terms instead. We conclude by introducing the terms

M m=1 u m G m and M m=1 u 2 m G 2
m in conjunction with Equation 3.29, 3.33 and 3.34:

E(CRPS(G E , y)) = H y + M m=1 u m G m -2 M m=1 u m G m H y + M m =k u m u k G m G k - M m=1 u m G m + M m=1 u 2 m G m = E(|X -y|) + M m =k u m u k G m G k - M m=1 u m G m + M m=1 u 2 m G m = E(|X -y|) + M m,k u m u k G m G k - M m=1 u m G m + M m=1 u 2 m G m - M m=1 u 2 m G 2 m = E(|X -y|) - 1 2 E(|X -X |) + 1 2 M m=1 u 2 m E(|X m -X m |) .
where X and X m are random variables with CDFs G and G m respectively. Therefore, we have

E(CRPS(G E , y)) = CRPS(G, y) + 1 2 M m=1 u 2 m E(|X m -X m |) . (3.7)
In the expectation of the ensemble CRPS, the diagonal terms u 2 m E |X m -X m | are missing, because the spread of each member is assumed to be null. The absence of the diagonal terms is the cause of the bias of the ensemble CRPS. One consequence of Equation 3.7 is that the expected ensemble score is never smaller than the score obtained with G. This fact also follows from the convexity of the CRPS, as pointed out by an anonymous reader.

As a consequence, the minimization of the ensemble CRPS should not be targeted because the solution of this optimization problem is not the underlying CDF of the verification. There is no contradiction with the strict propriety of the CRPS because, for the ensemble CRPS, the solution is only searched in a subspace made of step functions.

In the case of equal weights with i.i.d. members, Fricker et al. [START_REF] Fricker | Three recommendations for evaluating climate predictions[END_REF] detailed in their Appendix why minimizing the ensemble CRPS is misleading as stated above. Ferro et al. [START_REF] Ferro | On the effect of ensemble size on the discrete and continuous ranked probability scores[END_REF] exhibit a fair adjusted CRPS score, which includes correction terms to counteract the bias:

CRPS a (G E , y) = 1 M M m=1 |x m -y| - 1 2 M m,k=1 |x m -x k | M (M -1) (3.8) = CRPS(G E , y) - 1 2M M m,k=1 |x m -x k | M (M -1) . (3.9)
We see that rather than being a new score, the adjusted ensemble CRPS is a better estimation of the original CRPS, where the underlying distributions of the members are taken into account. In Equation 3.8, the dispersion of the ensemble E

(|X -X |) is estimated by M m,k=1 |x m -x k |/(M (M -1)).
In other terms, the bias terms

u 2 m E(|X m - X m |) of Equation 3.7 are taking into account in Equation 3.9 as E(|X -X |)/M 2 , by considering that E(|X m -X m |) = E(|X -X |).

Mixture model described by classes of members

We propose in this section a framework compatible with both ensemble forecasting and unbiased scores. In a standard model mixture design, a forecaster will assign weights to known parametric distributions [START_REF] Raftery | Using Bayesian model averaging to calibrate forecast ensembles[END_REF][START_REF] Grimit | The continuous ranked probability score for circular variables and its application to mesoscale forecast ensemble verification[END_REF]. We do not want to make assumptions on distributions, thus we use a standard ensemble forecasting framework, where the members are usually assumed to be sampled from unknown CDFs. The goal of this section is to show that despite the finite size of the ensemble, it is possible to use the CRPS by counteracting the discretization-induced bias. This framework is close to what is introduced in Fraley et al. [START_REF] Fraley | Calibrating multimodel forecast ensembles with exchangeable and missing members using Bayesian model averaging[END_REF], however this previous work focused on Bayesian Model Averaging (BMA), and did not include considerations on the CRPS.

We assume that ensemble members are grouped into classes within which members are i.i.d. In this new setting, a class C has a weight W C uniformly distributed among its members. The weight u m = W C /M C is assigned to the mth member of the ensemble, assuming that it belongs to class C and that class C has M C members. As an example, classes may be defined according to the rank of the members. Assuming that 10 members are available, two classes may be built by assigning the 5 members with the lowest values to the first class and the remaining members to the second class.

We introduce the CRPS using the classes. We call this score the class CRPS, and denote it

CRPS C (G C , y) = C∈C W C E(|X C -y|) - 1 2 C,D∈C W C W D E(|X C -X D |) . (3.10)
The terms of the class CRPS are detailed below.

For the class C, with M C members x C c associated to the random variables X C and X C , we have

E(|X C -y|) = M C c=1 |x C c -y|/M C , (3.11) E(|X C -X D |) = M C c=1 M D d=1 |x C c -x D d |/(M C M D ) , (3.12)
where class D is different from class C, and Then we sum inter-and intra-class dispersions to link them to inter-member differences

E(|X C -X C |) = M C c,c =1 |x C c -x C c |/(M C (M C -1)) . ( 3
|x m -x k |.
The key point is that inter-member differences for i.i.d. members comprise the intra-class dispersions. We note that

W 2 C M C c,c =1 |x C c -x C c | M C (M C -1) = M C M C -1 M C c,c =1 W C M C 2 |x C c -x C c | , (3.15)
and

M C /(M C -1) = 1 + 1/(M C -1) to obtain C,D∈C W C W D E(|X C -X D |) = C =D∈C W C W D E(|X C -X D |) + C∈C W 2 C E(|X C -X C |) = M m,k=1 u m u k |x m -x k | + C∈C 1 M C -1 M C c,c =1 W 2 C M 2 C |x C c -x C c | = M m,k=1 u m u k |x m -x k | + C∈C W 2 C M C E(|X C -X C |) = M m,k=1 u m u k |x m -x k | + M m=1 u 2 m E(|X Cm -X Cm |) ,
where C m is the class in which x m falls. To obtain the last equation, consider that E

(|X C -X C |) is counted M C times.
Compared to the ensemble CRPS, the class CRPS admits M additional terms corresponding to the dispersion of each member and resulting from the classes definition:

CRPS C (G C , y) = CRPS(G E , y) - 1 2 M m=1 u 2 m E(|X Cm -X Cm |) . (3.16)
In the case of a single class, the class CRPS is equal to the adjusted ensemble CRPS described in Section 3.1.4.

The i.i.d. assumption on the members can be seen as too strong. The exchangeability of the members is a relaxation of the i.i.d. assumption. By definition, the joint distribution function of exchangeable members is invariant under permutation of the arguments, thus the members are indistinguishable. We refer the reader to Ferro [START_REF] Ferro | Fair scores for ensemble forecasts[END_REF] for an analysis of fair scoring rules with the exchangeability assumption. In a few words, the user must investigate the (generally unknown) dependence structure and taylor the appropriate scoring rule accordingly. The simple case of pairwise uncorrelated members is however tractable. For the ensemble CRPS, the case of pairwise uncorrelated members is in practice equivalent to the case of i.i.d. members, because the terms |x m -x k | rely on pairwise correlations only. In the same way for the class CRPS, the assumption of pairwise uncorrelated members within each class and independent members between classes leads to similar results than i.i.d. members. Under the more general assumption of exchangeable members within each class, the definition of E(|X C -X C |) should take into account the dependence between members.

Also note that these assumptions are only needed to counter the bias in the ensemble CRPS. Our aggregation methods still remain applicable without such correction. The theoretical bounds described in the next section do not rely on any stochastic assumption on the prediction data and the verifications. The assumptions of i.i.d. members and the use of the class CRPS should only guide the choice of a loss function.

Online learning methods

Theoretical background

Up to this section, a single time t was considered. Now we introduce online learning techniques. In this setting, the forecaster receives prediction data D t and wishes to produce the best prediction of y t . In our case, prediction data are ensemble members and the algorithm gives a rule to compute the weights u m,t before each forecast time t. This rule takes into account only past information, and is therefore called the update rule. The goal of a given online learning algorithm is to provide the best possible weights according to a chosen loss function, for example the ensemble CRPS

CRP S E t (u) = M m=1 u m H m,t -H yt 2 ,
(3.17)

written above for time t. The notation t (u) emphasizes the importance of the weights, as opposed to the ensemble members and the verifications which are assumed to be given to the forecaster.

In practice, the algorithm reads Initialization: u 1 ;

For each time index t = 1, 2, ..., T 1. get prediction data D t , 2. compute the forecaster's choice with D t and u t , 3. get the verification y t and compute u t+1 , based on the update rule.

The initial weight vector u 1 is arbitrarily set, e.g., to [1/M, . . . , 1/M ] . The performance of an update rule comes with theoretical guarantee, where the forecaster's results are assessed against a reference, which is usually the best forecast with weights constant in time, called the oracle. An important aspect of these theoretical guarantees is that they come with essentially no stochastic assumption on the prediction data and the verifications. In this paper, the theoretical guarantees are regret bounds of the form

T t=1 t (u t ) -inf u∈P M T t=1 t (u) o(T ) , (3.18)
where t is assumed to be bounded. The bound of t can be arbitrarily small or large, so that this restriction is compatible with essentially all real world applications. Averaging the losses in time (i.e., dividing by T ) shows that an algorithm giving the weights u t is guaranteed to perform at least as well as any mixture model with weights constant in time and based on the same prediction data. This includes any individual forecast and any subset ensemble with uniform weights.

We now consider two algorithms: the online ridge regression and the exponentiated gradient method (EG). We introduce these methods in a general framework, and we show how the methods can be applied to the case of the CRPS. For the algorithm run with ensemble CRPS, a weight is explicitly given to each member. The quantities |x m,t -y t | and |x m,t -x k,t | are explicitly used in the minimization process. For the algorithm run with class CRPS, equal weights are given to all the members within a class. The weights W C,t are computed using the terms E(|X C,t -y t |) and E(|X C,t -X D,t |). Combining parameterized distributions is also possible with online learning techniques. It necessitates to compute the quantities E(|X m,t -y t |) and E(|X m,t -X k,t |). These quantities are tractable from the CDFs using Equation 3.29. They are computed in Grimit et al. [START_REF] Grimit | The continuous ranked probability score for circular variables and its application to mesoscale forecast ensemble verification[END_REF] for a Gaussian mixture distribution.

Ridge regression

The approach of the ridge regression can be directly expressed in terms of minimization. The update rule for time t + 1 and based on the loss is

u t+1 = argmin w∈R M λw w + t t =1 t (w) . (3.19)
The regularization term with parameter λ 0 controls the 2-norm of the weight vector. It is possible to add discount factors in the sum of the past losses, in order to give more importance to recent timesteps. At first sight, the ridge regression does not constrain the weights to be positive or sum to one. In practice, for the CRPS, we observed that these constraints are approximately satisfied after a spin up period.

Other regularization terms of the form λ(w -u 1 ) (wu 1 ) may also be used with arbitrary reference vector u 1 ∈ P M . The reader interested in recent advances in online regularized regression is addressed to Orabona et al. [START_REF] Orabona | A generalized online mirror descent with applications to classification and regression[END_REF].

For a given experiment length T , for any vector u ∈ P M , and if the CRPS losses t (u t ) are bounded, we have:

R T (u) = T t=1 t (u t ) - T t=1 t (u) O (ln T ) , (3.20)
so that the so-called regret R T (u) is sublinear. The Appendix 3.C details technical aspects, such as the proof for the bound 3.20, as well as guidelines to compute the weights. The ridge regression applied to the square loss (E(X) -y) 2 gives a similar regret bound in terms of square losses. Our proof for the CRPS is inspired from the proof of the regret bound for the square loss, concisely described by Cesa-Bianchi and Lugosi [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF]. We were helped by the work of Mallet et al. [START_REF] Mallet | Description of Sequential Aggregation Methods and their Performances for Ozone Ensemble Forecasting[END_REF], who demonstrated the case of multiple verification locations (also called stations) for the square loss. Our work transposes the results for the square loss with multiple locations to multiple Brier score with different thresholds, and to the CRPS.

Method Gradient loss

Ensemble CRPS

m,t = |x m,t -y t | - M k=1 u k,t |x m,t -x k,t | + y t - M k=1 u k,t x k,t Class CRPS C,t = E(|X C,t -y t |) - D∈C W D,t E(|X C,t -X D,t |) + y t -E(X t )
CRPS for general mixture models

m,t = E(|X m,t -y t |) - M k=1 u k,t E(|X m,t -X k,t |) + y t -E(X t )
Table 3.1 -Formulae of the loss gradients. Equations from Appendix 3.A are used for the simplifications. The terms of the form y t -E(X t ) do not impact the computation of the weights for EG, because they are independent of the member m or the class C.

Exponentiated gradient

Let the learning rate η be strictly positive, EG follows a multiplicative update rule of the form:

u m,t+1 = u m,t exp -η m,t M m =1 u m ,t exp -η m ,t , (3.21)
where

m,t = ∂ t ∂u m (u t ) . (3.22)
This update relates to Bayesian inference [START_REF] Catoni | Statistical Learning Theory and Stochastic Optimization, Lectures on Probability Theory and Statistics, Ecole d'été de Probabilités de Saint-Flour XXXI-2001[END_REF][START_REF] Audibert | Fast learning rates in statistical inference through aggregation[END_REF]. The algorithm EG admits a formulation in terms of cost function minimization, where the regularization function is the entropy function, also known as the Kullback-Leibler divergence [START_REF] Kivinen | Exponentiated Gradient versus Gradient Descent for Linear Predictors[END_REF]. The EG algorithm reads:

u t+1 = argmin w∈P M M m=1 w m ln( w m u m,t ) + ηw m m,t . (3.23)
Examples of loss gradients are provided in Table 3.1. The loss gradient m,t of the CRPS has two main terms: (i) E(|X m,t -y t |) accounting for the distance between the verification and the mth random variable X m,t , and (ii) the weighted sum of E(|X m,t -X k,t |) accounting for distances between X m,t and the X k,t . The first term controls a deviation from the median of the underlying distribution of the verifications, and the second term controls the dispersion of the mixture model. On average (on the verifications), the loss gradients are null if the verifications are correctly described by the forecaster's CDF.

The advantage of using the loss gradients is described (at least) in Devaine et al. [START_REF] Devaine | Forecasting electricity consumption by aggregating specialized experts[END_REF]. In a few words, using the loss gradients makes the algorithm competitive against the best convex combination with constant weights, whereas simply using the loss m,t = E(|X m,t -y t |) -0.5 E(|X m,t -X m,t |) would make the algorithm compete only against the best member. We insist on the fact that using the loss gradients provides the terms E(|X m,t -X k,t |) which are critical for the control of the ensemble spread. The theoretical guarantee for EG states that, if the loss function is convex with respect to u and admits a subgradient, and if the losses m,t are bounded within a constant interval [-a, a], then we have:

sup T t=1 t (u t ) -inf u∈P M T t=1 t (u) ln M η + η a 2 2 T , (3.24)
where the supremum is taken for all possible values of the members x m,t and the verifications y t , and η is the learning rate [START_REF] Devaine | Forecasting electricity consumption by aggregating specialized experts[END_REF]. For optimized values of η proportional to 1/ √ T , the regret is sublinear. The theoretical guarantee of Equation 3.24 is verified for the square loss and for the CRPS.

Numerical example

Simple model

We use the simple model described in Bröcker [START_REF] Bröcker | Evaluating raw ensembles with the continuous ranked probability score[END_REF]. The model is supposed to mimic local temperatures. We chose this model because the uncertainty terms are known, consequently we can easily draw conclusions from numerical tests.

We built the verifications y t from the exact time series

a t = (Asin(πω 1 t) + Bsin(πω 2 t)) 2 , (3.25)
combined with multiplicative and additive perturbation terms:

y t ∼ a t (1 + s 1 N (0, 1)) + s 2 N (0, 1) . (3.26)
Each term N (0, 1) represents an independent Gaussian noise with zero mean and a variance of one. The perturbation terms are sampled independently at each timestep.

The parameters are summarized in Table 3.2. The members are sampled as

x m,t ∼ a t (1 + s 1 N (0, d ens )) + s 2 N (0, d ens ) (3.27)
analogously to the verification distribution, but the standard deviation d ens describing the perturbations terms may differ from its optimal value (i.e., 1). The parameter d ens is also referred to as the dispersion parameter.

Experiments without online learning

In this first experiment, ensembles are built for different values of the dispersion parameter d ens . The members are drawn independently, and the weights of the members are taken constant and all equal to 1/M . As expected, the adjusted ensemble CRPS gets the lowest value when the ensemble shows the correct spread (i.e., for d ens = 1), see Figure 3.1. On the contrary, the best (non adjusted) ensemble CRPS is obtained for under-dispersed ensembles d ens < 1. The shift of the ensemble CRPS minimum from the ideal location d ens = 1 is larger for ensembles of small size, because the bias of the ensemble CRPS is proportional to 1/M . This is a direct illustration of the bias due to the limited size of the ensemble explained in Section 3.1.4.

Experiments with weight updates

Now we test online learning techniques and more specifically their ability to discriminate between members. We build an ensemble of M = 10 members, that is composed of two subensembles, or classes, of equal size. The first subensemble is defined by the same distribution than the verifications. The second subensemble follows a distribution controlled by d ens . If d ens = 1, then the whole ensemble is correctly dispersed. In other words, half of the members follow the correct distribution, while the second half can follow a different distribution.

An example of the temporal evolution of the weights is given in Figure 3.2. We used the algorithm EG (η = 0.05) with the gradients of the ensemble CRPS. At the middle of the experiment, we swap the dispersion parameters of the members. Correct members become incorrect members and conversely. The members with incorrect dispersion parameter (d ens = 1.5) see their weights decrease on average. After the swap, the weights of the newly incorrect members also decrease on average. The impact of the learning rate is shown in Figure 3.3, where a larger value η = 0.2 leads to a faster evolution of the weights. Note the difference of scales between Figures 3.2 The black crosses indicate that the CRPS of each member are used in EG (i). The light gray triangles indicate that the square loss gradients are used in EG (ii). This figure shows that not using the CRPS gradients favors the less dispersed members, even though they do not show the correct dispersion. The experiment of roughly ten years is repeated 200 times for each dispersion parameter. We used the learning parameters η = 0.05 and λ = 0.5.

for different learning algorithms. Here we did not include a change of the dispersion parameter at mid-experiment. The first subensemble therefore remains the correct one all the time. The discrimination procedure tests whether the algorithm makes a difference between the subensembles and whether the incorrect members are given lower weights than the correct members.

We show the importance of the CRPS gradients in EG for probabilistic forecasting. We show in Figure 3.4 the average weights of EG using: (i) m,t = |x m,t -y t |, using the CRPS without the gradients; or (ii) ˇ m,t = 2(u t x t -y t )x m,t , using the gradients of the square loss (u t x t -y t ) 2 , instead of the CRPS gradients m,t . We see that in either case, the members with the lowest dispersion parameter are the most weighted. The members with the correct distribution receive the highest weights when the incorrect members are over-dispersed (d ens > 1). Formulation (i) and (ii) do not tend to forecast the distribution of the verifications, but only the mean or the median of the distribution of the verifications. These formulations are therefore not suited for probabilistic forecasting, as opposed to the CRPS gradients (see below). Note that we can rewrite ˇ m,t = (x m,t -y t ) 2 -(x m,t -u t x t ) 2 plus terms independent of m. Thus using the gradients (or equivalently trying to get the best combination) is a diversification strategy compared to simply using (x m,t -y t ) 2 .

Using the same representation, the algorithms EG and the ridge regression are tested with the ensemble CRPS and the class CRPS, see Figure 3.5. The algorithms based on the class CRPS show correct discrimination: whatever the dispersion parameter of the wrongly dispersed members, the class with incorrect dispersion shows smaller weights on average. The sum of the weights attributed to the incorrect members stays below 0.5 (equal weights between the two subensembles). On the contrary, the algorithms based on the ensemble CRPS does not give a correct discrimination. When the dispersion parameter d ens is close to 0.70, the under-dispersed members receive larger weights than the correct members. We see that the minimization of the ensemble CRPS is misleading for an ensemble of small size. We interpret these results as direct consequences from the bias of the ensemble CRPS described in Section 3.1.4.

Conclusion

We introduced new tools for probabilistic forecasting using an ensemble of forecasts. Our algorithms use online learning techniques to produce forecast combinations that tend to minimize the CRPS. In the long run, they guarantee that the performance of the weighted ensemble is at least as good as the performance of the best weighted ensemble with weights constant in time. This theoretical guarantee holds with essentially no assumption on the distributions of the forecasts and verifications. In this sense, our method is a non-parametric post-processing method.

A new framework using classes of members is introduced in order to counteract the bias in the ensemble CRPS. With this framework and the proposed algorithms, numerical tests showed that our online learning techniques tend to give higher weights to the forecasts with the same distribution as the verifications.

The algorithms should now be tested against real data, in order to assess their potential in operational applications against Bayesian model averaging (BMA) or other post-processing techniques. The work of the forecaster is then to obtain numerous forecasts to combine. The methods do not require any assumptions on the forecasts to be applied (bias, spread, or any other stochastic or deterministic assumptions), except the loss boundedness. However, some good practices may be applied to improve the overall performance. For example, the forecasts can be altered before their inclusion in the ensemble, or additional forecasts may be derived from the raw ensemble. Also, it is recommended to draw ensembles with enough spread, so that they encompass the verifications. We argue that for most applications, the use of a multimodel ensemble combined with several post-processing techniques is an efficient way to obtain an ensemble to be calibrated with our algorithms. From a meteorological point of view, new members can be added to the ensemble by using nearby grid-points or time-shifted forecasts. This approach may be particularly efficient to account for the ability of a forecasting system to predict an event, but at the wrong time or location.

On theoretical side, a next step could be the inclusion of the uncertainty in the verifications. Also, other non local strictly proper scoring rules could serve as loss function.

Appendix 3.A Identities implying CDFs

Let the random variable Z be described by the probability density function K and the CDF K. We have for any real number x:

E(H(x -Z)) = K (Z)H(x -Z)dZ = K(x) , ( 3.28)
or equivalently E(H Z ) = K. The demonstration of the strict propriety of the CRPS uses this property for the integration over the CDF of the verifications. Let X and Z be two random variables described respectively by the CDFs G and K. We have:

E(|X -Z|) = G(1 -K) + K(1 -G) .
(3.29)

For G = K, the above quantity is the Gini mean difference, which is thoroughly introduced in the monograph of Yitzhaki and Schechtman [START_REF] Yitzhaki | The Gini Methodology: A primer on a statistical methodology[END_REF].

The product GK of CDFs is itself the CDF of the random variable max(X, Z). This can be used to explain simply Equation 3.29, using:

2 max(a, b) = |a -b| + a + b , (3.30)
for any (a, b) ∈ R 2 , and

E(Z) = +∞ -∞ (H(x) -K(x))dx . (3.31)
Let G = I i=1 u i G i and K = J j=1 w j K j be CDFs of mixture models with respectively I and J components, i.e., the G i and K j are CDFs, and the weight vectors u and w respectively belong to the simplexes P I and P J . Let X, X i , Z and Z j be random variables respectively following G, G i , K, and K j . We have

E(|X -Z|) = I i=1 J j=1 u i w j E(|X i -Z j |) , (3.32)
based on Equation 3.29. Indeed,

I i=1 J j=1 u i w j (G i (1 -K j ) + K j (1 -G i )) = I i=1 u i G i (1 - J j=1 w j K j ) + J j=1 w j K j (1 - I i=1 u i G i ) = E(|X -Z|) ,
because the weights w i and u j respectively sum to one, and using the linearity of integration.

It is straightforward to use Equation 3.32 to show to that

E (|X -y|) = M i=1 u i E (|X i -y|) , (3.33)
and that

E |X -X | = M i=1 u i E (|X i -X|) = M i,j=1 u i u j E |X i -X j | , (3.34)
with X and X j being random variables respectively described by G and G j .

because the weights u m sum to one. We highlight the fact that the diagonal terms u 2 m |x m -x m | are null, so that the double sum of Equation 3.35 is computed for m = k. The calculus of this section can also be written with expectations and random variables using the content of Appendix 3.A.

Appendix 3.C Regret bound of the ridge regression with the CRPS

This section is written for general model mixtures Gm, t and for general CDF F t for the verifications. For simplicity, we assume that the integrals of the CRPS can be computed on an interval [γ, Γ] of limited size. All the considered CDFs hit 0 at γ and 1 at Γ, which formalizes the assumption of bounded values for the members and the verifications. Thus the considered CDFs verify G m,t Γ -γ .

This appendix is structured as follows: (i) we exhibit an update rule between u t+1 and u t ; (ii) we bound the regret against a constant vector T t=1 t (u t ) -T t=1 t (u) by the regret against the best a posteriori vector T t=1 t (u t ) -T t=1 t (u t+1 ); (iii) we provide an interpretable regret bound by using the update rule and the convexity of t .

The CRPS has a quadratic form

t (u) = u G t G t u -2u F t G t + F 2 t , (3.36)
where F t (x) = H(x -y t ) and G t (x) is the vector of the CDFs G m,t (x). The cost function J t (u) = λu u + t t =1 t (u) is written in a quadratic matricial form with:

J t (u) = u A t+1 u -2u b t+1 + t t =1 F 2 t , (3.37)
where the vector b t is defined by:

b t = t-1 t =1 F t G t , ( 3.38)
and the matrix A t of size M × M is symmetrical positive-definite:

A t = λI M + t-1 t =1 G t G t , (3.39)
with I M the identity matrix. The matrix A t admits an inverse which is also symmetrical positive-definite. Note the trivial recurrence relation

J t+1 = t+1 + J t .
The weight u t+1 is by definition the minimizer of J t . Simple derivation gives the equality A t u t = b t . In practice, the weights are found via matrix inversion. Besides, a recurrence relation can be obtained. We successively deduce:

A t+1 u t+1 = b t+1 = b t + F t G t , = A t u t + F t G t , = A t+1 -G t G t u t + F t G t .
The recurrence relation holds for any quadratic definition of the loss , and is expressed as:

u t+1 -u t = A -1 t+1 F t -u t G t G t = - 1 2 A -1 t+1 ∇ t (u t ) . (3.40)

Demonstration of the regret bound

We iteratively use the fact that u t+1 is the minimizer of J t to get

J T (u) J T (u T +1 ) = T (u T +1 ) + J T -1 (u T +1 ) T (u T +1 ) + J T -1 (u T ) T t=1 t (u t+1 ) + λu 1 u 1 . (3.41)
The nonnegativity of λu 1 u 1 gives:

T t=1 t (u) T t=1 t (u t+1 ) -λu u . (3.42)
Thus the regret can be bounded:

R T (u) = T t=1 t (u t ) -t (u) λu u + T t=1 t (u t ) -t (u t+1 ) λu u + T t=1 (∇ t (u t )) (u t -u t+1 ) = λu u + 1 2 T t=1 (∇ t (u t )) A -1 t+1 ∇ t (u t ) , (3.43)
where we have used Equation 3.42, the convexity of the functions t and Equation 3.40. At this point of the demonstration, one may have the feeling that a logarithm bound can be obtained, because the matrix A t is a sum of t matrices, and because the logarithmic function is the primitive of the inverse function.

We define

Q t = A -1/2
t+1 G t and s t = (u t G t -F t ), so that the symmetry of

A -1/2 t+1 gives 1 2 (∇ t (u t )) A -1 t+1 ∇ t (u t ) = 2 s t Q t s t Q t .
The inequality of Cauchy-Schwartz gives

s t Q t s t Q t = M m=1 s t Q t m 2 M m=1 s 2 t [(Q t ) m ] 2 = s 2 t Q t Q t = t (u t ) G t A -1 t+1 G t . (3.44)
We continue with

G t A -1 t+1 G t = Tr A -1 t+1 G t G t = Tr A -1 t+1 G t G t = Tr I M -A -1 t+1 A t ln det A t+1 det A t . (3.45)
The first equality holds with the linearity of the integration and because z 1 Az 2 = Tr(Az 2 z 1 ) for any vectors z 1 , z 2 and matrix A. The inequality holds because A -1 t+1 A t is positive definite and 1 -1/x ln x for any x > 0.

At this step of the proof, we have shown that:

R T (u) λu u + 2 T t=1 t (u t ) ln det A t+1 det A t . (3.46)
We assume that the losses t (u t ) are bounded by a > 0. Then we easily reach:

R T (u) λu u + 2a ln det A T +1 λ M . (3.47)
The inequality of arithmetic and geometric means applied to the eigenvalues of A T +1 leads to the conclusion

det(A T +1 ) Tr A T +1 M M = M λ + T t=1 Tr G t G t M M M λ + M T (Γ -γ) M M , ( 3

Extension to threshold-weighted and quantile-weighted scoring rules

As introduced in Section 1.3, various ways exist to evaluate forecasts of binary events. A infinite amount of strictly proper scoring rules can be written as a sum of elementary quantile losses. Besides, a second summation on event thresholds allows to evaluate probabilistic forecasts of a scalar variable. In this section, the notation α refers to levels of quantile and the notation θ refers to thresholds. The interested reader is referred to the work of Buja et al. [START_REF] Buja | Loss functions for binary class probability estimation and classification: Structure and applications[END_REF], Gneiting and Ranjan [START_REF] Gneiting | Comparing density forecasts using thresholdand quantile-weighted scoring rules[END_REF], and Ehm et al. [START_REF] Ehm | Of quantiles and expectiles: consistent scoring functions, Choquet representations and forecast rankings[END_REF] for state-of-the-art articles on this subject, and Dawid [START_REF] Dawid | Comments on: Assessing probabilistic forecasts of multivariate quantities, with an application to ensemble predictions of surface winds[END_REF] for a detailed bibliographical analysis of assessments of probabilistic forecasts. General conditions ensuring strict propriety are not new since they were addressed by Shuford et al. [START_REF] Shuford | Admissible probability measurement procedures[END_REF], Savage [START_REF] Savage | Elicitation of personal probabilities and expectations[END_REF], and Schervish [START_REF] Schervish | A General Method for Comparing Probability Assessors[END_REF], among others. The choice between these scoring rules was also recently addressed by Merkle and Steyvers [MS13] and Lerch et al. [START_REF] Lerch | Forecaster's dilemma: Extreme events and forecast evaluation[END_REF].

In the remaining of this short introduction, we briefly review the score decomposition of Ehm et al. [START_REF] Ehm | Of quantiles and expectiles: consistent scoring functions, Choquet representations and forecast rankings[END_REF] for the celebrated CRPS. In Section 4.1.1, we generalize the results of Thorey et al. [START_REF] Thorey | Online learning with the CRPS for ensemble forecasting[END_REF] to threshold-weighted scoring rules. Quantile-weighted scoring rules are analyzed in the context of mixture models in Section 4.1.2.

Let G be a CDF delivered to forecast the observation y. The CRPS can be written

CRPS(G, y) = 2 S α,θ (G -1 (α), y)dαdθ , ( 4.1)
if G is an invertible function. The integration holds on the levels α ∈ [0, 1] of quantile and on the thresholds θ ∈ R. The function S α,θ (x, y) is defined by

S α,θ (x, y) = (H(x -y) -α)(H(x -θ) -H(y -θ)) (4.2) =      1 -α if y θ < x , α if x θ < y , 0 otherwise.
The forecast x of the observation y is evaluated with the score S α,θ (x, y) for the event y < θ. The quantity α is the cost of a false positive, and 1 -α is the cost a false negative.

Integrating over the levels α gives the standard formulation of the CRPS, indeed

2 S α,θ (G -1 (α), y)dα = 2 θ (1 -G(x))g(x)dx if y θ , 2 θ G(x)g(x)dx if θ < y , = (H(θ -y) -G(θ)) 2 ,
with the change of variable G(x) = α.

Integrating over the thresholds relates the CRPS to the quantile score QS α (x, y) (or pinball loss) of level α [START_REF] Koenker | Regression quantiles[END_REF]. The quantile decomposition of the CRPS was found by Laio and Tamea [START_REF] Laio | Verification tools for probabilistic forecasts of continuous hydrological variables[END_REF], and developed by Bröcker [START_REF] Bröcker | Evaluating raw ensembles with the continuous ranked probability score[END_REF] for ensemble forecasts. Using (H(x -θ) -H(y -θ))dθ = x -y, we have G = E(G E ), where the expectation is taken over the members, describes the random variable X. With Equation 4.5, we swap the members x m for Φ(x m ) and the verification y for Φ(y) in the expression of the weighted CRPS:

S α,θ (x, y)dθ = (H(x -y) -α)(x -y) = QS α (x, y) ,
wCRPS(G E , y) = M m=1 u m |Φ(x m ) -Φ(y)| - 1 2 M m,k=1 u m u k |Φ(x m ) -Φ(x k )| . (4.6)
Therefore, the bias of the threshold-weighted CRPS is explicitly found in

E(wCRPS(G E , y)) = E(|Φ(X) -Φ(y)|) - 1 2 E(|Φ(X) -Φ(X )|) + 1 2 M m=1 u 2 m E(|Φ(X m ) -Φ(X m )|) = wCRPS(G, y) + 1 2 M m=1 u 2 m E(|Φ(X m ) -Φ(X m )|) .
Also for the class CRPS, our results are extended to the class threshold-weighted CRPS by introducing the quantities of interest

M C c=1 |Φ(x C c )-Φ(y)|/M C , M C c=1 M D d=1 |Φ(x C c )- Φ(x D d )|/(M C M D ), and M C c,c =1 |Φ(x C c ) -Φ(x C c )|/(M C (M C -1)
). In order to give a better understanding on the effect of the Φ-transformation, we give the following example. Consider the case of a step function φ defined by: 4.1). We assume that all observations and members are in [0, 1]. Roughly speaking, in our example, errors in [0, z 0 ] cost more than errors in [z 0 , 1]. In this setting, we have

φ(x) =      1 if 0 x < z 0 κ if z 0 x 1 0 otherwise , with κ < 1, z 0 ∈ [0, 1] (see Figure
Φ(x) = x if 0 x < z 0 z 0 + κ(x -z 0 ) if z 0 < x 1 , and |Φ(x m ) -Φ(x k )| =      |x m -x k | if 0 x k , x m z 0 κ|x m -x k | if z 0 x k , x m 1 |z 0 -x m | + κ|x k -z 0 | if 0 x m z 0 x k 1 .
In this simple example, it is easy to see that the distance |Φ(x m ) -Φ(x k )| is high when x m and x k are in the domain of high value of φ. And conversely, the distance |Φ(x m ) -Φ(x k )| is low when x m and x k are in the domain of low value of φ. This is consistent with the intuition that the function φ defines domains with high or low importance.

Effect of quantile-weighting

We prove in this section that quantile weighting can also be understood as a data transformation, at least for the score gradients. We show that the data transformation due to the quantile weighting is related to both the weighting function ω and the forecaster's CDF G. We begin this section by rewriting quantile-weighted scores in a convenient way and studying the quantile-weighted score of an ensemble of forecasts.

Rewriting quantile-weighted scoring rules

Let the score S(G, y) be a quantile-weighted score defined by

S(G, y) = QS α (G -1 (α), y)ω(α)dα = QS G(X) (X, y)ω(G(X))dG(X) = E X QS G(X) (X, y)ω(G(X)) ,
where ω is strictly positive on open intervals in [0, 1]. Propriety is retained since the minimizer of E Y [S(G, Y)] verifies G(α) = F(α) for each 0 < α < 1, where Y is a random variable described by the CDF F.

From the definition of S α,θ in Equation 4.2, we have

S(G, y) = H(θ -y) 1 G(θ) (1 -α)ω(α)dα + H(y -θ) G(θ) 0 αω(α)dα dθ .
(4.7) We used the conditions y θ < x or equivalently G(y) G(θ) < α for the left term, and the conditions x θ < y or equivalently α G(θ) < G(y) for the right term.

The score S(G, y) can therefore be expressed as

S(G, y) = H y β 1 (G) + (1 -H y )β 0 (1 -G) , (4.8) by identifying β 1 (G) = 1 G (1 -α)ω(α)dα and β 0 (1 -G) = G 0 αω(α)dα .
For example with the CRPS, we have

β 1 (G) = (1 -G) 2 and β 0 (1 -G) = G 2 .
If ω is symmetrical with respect to 1/2, then β 1 = β 0 . Equation 4.7 proposes a formulation for quantile-weighted scoring rules emphasizing the interplay between ω and G. This formulation is particularly useful when analytical expressions of β 1 and β 0 are available. The linearity in H y is convenient to derive the uncertainty term of F and the divergence term between F and G in:

E Y [S(G, Y)] = Fβ 1 (G) + (1 -F)β 0 (1 -G) = F(β 1 (G) -β 1 (F)) + (1 -F)(β 0 (1 -G) -β 0 (1 -F)) divergence = d(F,G) + Fβ 1 (F) + (1 -F)β 0 (1 -F) uncertainty = e(F)
.

We highlight the fact that β 1 (G)-β 1 (F) may be rewritten under the form F G (1 -α)ω(α)dα and similarly for β 0 (1 -G) -β 1 (1 -F). This gives an explicit formulation of the divergence term. Besides, the uncertainty term is equal to E

Y [S(F, Y)].
By considering a distribution on the forecasts G, the above decomposition may be led further, as for any strictly proper scoring rule [START_REF] Bröcker | Reliability, sufficiency, and the decomposition of proper scores[END_REF]. To do so, we note F = E G [F] the climatological distribution, where the expectation is taken with respect to the frequency at which G is predicted (F is implicitly conditioned on G). The reliability E G [d(F, G)] measures the match between the forecasted probabilities G(θ) and the conditional frequencies F(θ). The resolution E G [d(F, F)] is the distance between the climatology F(θ) and the conditional frequencies F(θ). Such separation dates back to Murphy [START_REF] Murphy | A New Vector Partition of the Probability Score[END_REF] for the Brier score, Hersbach [START_REF] Hersbach | Decomposition of the continuous ranked probability score for ensemble prediction systems[END_REF] for the CRPS, and Bentzien and Friederichs [START_REF] Bentzien | Decomposition and graphical portrayal of the quantile score[END_REF] for quantile scores. The decomposition writes

E G,Y [S(G, Y)] = e(F) uncertainty + E G [d(F, G)] reliability -E G [d(F, F)] resolution .
For any quantile-weighted scoring rule, the explicit formulations of the uncertainty e and the divergence d give the score decomposition into uncertainty, reliability and resolution.

We finish this section with examples of quantile-weighted scoring rules. Several weighting functions ω of the beta family including the CRPS are defined in Table 4.1, and represented in Figure 4.2. The consequence of the choice of ω is shown in Figure 4.3, which illustrates the functions H y β 1 (G) + (1 -H y )β 0 (1 -G) in two situations: G (1) is the CDF of a Gaussian distribution with mean µ = 1 and variance σ 2 = 1, and G (2) is the CDF of a Gaussian distribution with mean µ = 0 and variance σ 2 = (2.577) 2 . In both cases, a null observation y = 0 is considered. Higher importance is given to the distribution tails when ω assigns higher weights to the extreme quantile levels near 0 and 1. Conversely, higher importance is given to the distribution center when ω assigns higher weights to the quantile levels close to 0.5. -G) for y = 0 and G the CDF of a Gaussian distribution with mean µ = 1 and variance σ 2 = 1 (left), and y = 0 and G the CDF of a Gaussian distribution with mean µ = 0 and variance σ 2 = (2.577) 2 (right). The weighting functions ω belong to the symmetrical beta family, ω(α) = α a-1 (1 -α) a-1 : a = 0 (violet), a = 1/2 (red), a = 1 (black), a = 2 (cyan). The functions are scaled to reach the same value at G -1 (0.5).

ω(α) β 1 (α) (α(1 -α)) -1 -ln(α) (CRIGN) (α(1 -α)) -1/2 arcsin( √ 1 -α) -α(1 -α) 2 (1 -α) 2 (CRPS) α(1 -α) (1-α) 3 3 -(1-α) 4 4 Table 4.1 -Examples of scoring rules in the symmetrical beta family ω(α) = α a-1 (1 -α) a-1 . a = 0 a = 1/2 a = 1 a = 2

Optimal quantiles locations for an ensemble of forecasts

In this section, we derive the optimal quantile positions of the M members x m with fixed weights u m when the combination G E of unit step functions is evaluated by a quantile-weighted scoring rule. Without lack of generality, we assume that the members are labeled according to their rank, i.e x m-1 < x m for any integer m ∈ [2, M ]. We note U m = k m u k with the convention U 0 = 0. The optimal location of the members is well-known for the CRPS [START_REF] Bröcker | Evaluating raw ensembles with the continuous ranked probability score[END_REF]. The CRPS expectation over the observations is minimized by members verifying F(x m ) = U m -um 2 . We show that the optimal locations for the CRPS are close to the optimal locations found for other quantileweighted scoring rules, except for the outer members x 1 and x M . For simplicity, the distribution of the observations is assumed to be bounded.

We recall that

E Y [S(G E , Y)] = Fβ 1 (G E ) + (1 -F)β 0 (1 -G E ) . (4.9)
Since G E is piecewise constant, the term depending on

x m in E Y [S(G E , Y)] is β 1 (U m-1 ) xm x m-1 F+β 1 (U m ) x m+1 xm F+β 0 (1-U m-1 ) xm x m-1 (1-F)+β 0 (1-U m ) x m+1 xm (1-F) ,
(4.10) for any integer m ∈ [2, M -1]. This expression is also valid for x 1 and x M if β 1 and β 0 are bounded because β 1 (0) = β 0 (0) = 0. In this case, we use the convention x 0 = sup F(x)=0 x and x M +1 = inf F(x)=1 x. If β 1 and β 0 are not bounded, the optimal location of the extreme members lies at the bounds of the observational distribution, i.e. x 1 = sup F(x)=0 x and x M = inf F(x)=1 x.

After differentiating Expression 4.10 with respect to x m , and setting the derivative to 0, the optimal location of the members are found with

F(x m ) = β 0 (1 -U m ) -β 0 (1 -U m-1 ) β 0 (1 -U m ) -β 0 (1 -U m-1 ) + β 1 (U m-1 ) -β 1 (U m ) . ( 4

.11)

A safety check for the CRPS indicates that Equation 4.11 is consistent with the results of Bröcker [START_REF] Bröcker | Evaluating raw ensembles with the continuous ranked probability score[END_REF]:

F(x m ) = U 2 m -(U m -u m ) 2 U 2 m -(U m -u m ) 2 + (1 -U m + u m ) 2 -(1 -U m ) 2 = u m u m × 2U m -u m 2U m -u m + 2 -2U m + u m = U m - u m 2 , with β 1 (α) = (1 -α) 2 and β 0 (1 -α) = α 2 .
The optimal quantiles of several scoring rules are shown in Figure 4.4. The location of the optimal members do not vary much depending on the scoring rule, except for the position of x 1 and x M . Larger variations are observed for a small amount of members M . The location of the members are barely distinguishable for M = 20. with weighting functions ω belonging to the symmetrical beta family, ω(α) = α a-1 (1 -α) a-1 : a = 0 (violet), a = 1/2 (red), a = 1 (black), a = 2 (cyan).

Model mixture and quantile-weighted scoring rules

Now we investigate the relationships between quantile-weighted scores and their gradients for model mixtures. Let G = m M u m G m be a model mixture. We compute the gradients of S(G, y) with Equation 4.7:

∂S(G, y) ∂u m = H y [-G m (1 -G)ω(G)] + (1 -H y )G m Gω(G) = G m (G -H y )ω(G) = 1/2 [(H y + G m -2H y G m ) -(G + G m -2GG m ) + G -H y ]ω(G) .
(4.12)

We check that we find the CRPS for ω = 2. In the Expression 4.12 of ∂S ∂um , we find a balance between how close G m is to H y against how close G m is to G. The distance is measured by terms of the form (K 1 + K 2 -2K 1 K 2 )ω(G), which can be seen as a generalization of the Gini mean difference with the additional ω(G) term (see Equation 4.5). Similarly to the CRPS, the distance between G m and G is a weighted distance between the mth members and the other members. Indeed we have

(G + G m -2GG m )ω(G) = M k=1 u k (G k + G m -2G k G m )ω(G) .
(4.13)

The quantile weighting can therefore be understood as a data transformation for the score gradients, recalling the results obtained for threshold-weighted scoring rules in Section 4.1.1. For quantile-weighted scoring rules, the weighting φ = ω(G) is applied in Equation 4.12 to the scores gradients. The data transformation related to the weighting ω(G) is its antiderivative x ω(G). Note that this data transformation depends on G.

To conclude this section, we study the effect of the quantile-weighting for ensemble forecasting, where G E and ω(G E ) are piecewise constant functions. This case is close to the example provided for threshold-weighted score in Section 4.1.1, with the piecewise constant weighting function φ. Let (z, z ) ∈ R 2 play the role of (x m , y) or (x m , x k ). We have

ω(α) ω(G E (x)) ω α, x 0 0.25 0.5 0.75 1 1 2 3 4 x 1 x 2 x 3 x 4 ω(α) ω(G E (x)) ω α, x 0 0.3 0.6 0.9 1 1 2 3 4 x 1 x 2 x 3 x 4
(H(x -z) + H(x -z ) -2H(x -z)H(x -z ))ω(G E (x))dx = | z z ω(G E )| , (4.14)
because the quantity (H(x -z) + H(x -z ) -2H(x -z)H(x -z )) = 1 if x is between z and z , and zero otherwise. Recalling that ω(G E ) is piecewise constant, we see that the gradient ∂S(G E ,y) ∂um is controlled by terms | z z ω(G E )| which are piecewise affine in z and z . The slope ω(G E (x)), with x ∈ [x m , x m+1 ], can be seen as the importance given to the interval [x m , x m+1 ].

In practice, the gradients may be computed as a sum of terms (x m+1 -x m )ω(U m ), accompanied by (x m * +1 -y)ω(U m * ) and (y -x m * )ω(U m * ), where x m * and x m * +1 are the closest members to the observation, i.e. x m * y x m * +1 .

We illustrate our example for M = 4 members, in Figure 4.5. We consider that x 1 < x 2 < x 3 < x 4 , without loss of generality. For this example, the function ω is convex. Higher importance is given to the levels of quantile close to 0 or 1 than those close to 0.5. The step function ω(G E ) gives a variable importance to the 3 intervals [x 1 , x 2 ], [x 2 , x 3 ] and [x 3 , x 4 ]. For uniformly distributed members u 1 = u 2 = u 3 = u 4 = 0.25, the importance ω(G E ) reaches the values ω(0.25) on the interval [x 1 , x 2 ], ω(0.5) on the interval [x 2 , x 3 ] and ω(0.75) on the interval [x 3 , x 4 ]. In our case where ω is symmetrical function (ω(α) = ω(1 -α)), the interval [x 2 , x 3 ] has the lowest importance and the importance of the intervals [x 1 , x 2 ] and [x 3 , x 4 ] are equal. We also show the case of non-uniform weights u 1 = 0.3, u 2 = 0.3, u 3 = 0.3, u 4 = 0.1, where the interval [x 1 , x 2 ] has a lower importance than the interval [x 3 , x 4 ]. Compared to the example of Section 4.1.1 with threshold-weighted score, the quantile-weighted score also gives a variable importance to domains of R. While the domains of variable importance are fixed for threshold-weighted score, the domains of variable importance rely on both u m and x m for quantile-weighted scores.

Illustrations of the loss gradients are proposed in Figure 4.6. We picked M = 101 linearly spaced members between -2 and 2 and show the loss gradients of the CRPS and the CRIGN for 3 observations y ∈ {0, 1.2, 1.9}. We see that the extreme levels of quantile receive relatively lower losses with the CRPS than with the CRIGN due to the higher concavity of the CRPS gradient loss (see the case y = 0). Besides, when the observation reaches a distribution tail (see the case y = 1.9), the relative magnitude of the CRIGN gradients is higher than the relative magnitude of the CRPS gradients. In other words, the weights u m learned by CRIGN minimization are more influenced by observations reaching the observations tails.

The study of the bias of quantile-weighted scores is left for future research. Using the notation G = E(G E ), the difficulty of the study of the bias in the general case is the comparison of E(β 1 (G E )) and E(β 0 (1 -G E )) to β 1 (G) and β 0 (1 -G).

Convexity of quantile-weighted scores

In order to use quantile-weighted scoring rules in sequential aggregation, we study the convexity of the score

S(G, y) = H y β 1 (G) + (1 -H y )β 0 (1 -G) ,
with respect to the weight vector u defining G = m M u m G m . We recall that convex loss functions are necessary to get regret bounds against the best fixed combination of experts. The convexity cannot be obtained by claiming a sum of convex functions, because the sum depends on the variable u which is concerned by the convexity.

From the mth gradient G m (G -H y )ω(G), we deduce the (m, k) component of the Hessian matrix A

A m,k = ∂ 2 S(G, y) ∂u m ∂u k = G m G k (ω(G) + (G -H y )ω (G)) , (4.15)
with ω being the derivative of ω. The score is convex if the Hessian matrix is positive definite, namely if for any vector w ∈ R M , w Aw 0.

The condition ω(G) + (G -H y )ω (G) 0 is sufficient to prove the convexity since we have

w Aw = M m,k=1 w m w k A m,k =   M m,k=1 w m w k G m G k   (ω(G) + (G -H y )ω (G)) = M m=1 w m G m 2 (ω(G) + (G -H y )ω (G)) .
We find the convexity of the CRPS in the case ω = 2. Besides, we prove the convexity of the CRIGN, obtained in the case ω(α) = 1/(α(1 -α)) = 1/α + 1/(1 -α). We recall the CRIGN definition: 2 , and

CRIGN = -H y ln G + (1 -H y ) ln(1 -G) . We have ω (α) = -1 α 2 + 1 (1 -α) 2 = 2α -1 α 2 (1 -α)
ω(G) + (G -H y )ω (G) = 1 G(1 -G) + (G -H y ) 2G -1 G 2 (1 -G) 2 = 1 G 2 (1 -G) 2 (G(1 -G) + (G -H y )(2G -1)) = 1 G 2 (1 -G) 2 (H y -G) 2 , 0 ,
which concludes the proof. The last equality holds since H y is either equal to 0 or 1.

The intermediate situation of ω(α) = 1/ α(1 -α) lies between the CRPS and the CRIGN. The related score is defined by

H y arcsin( √ 1 -G) + (1 -H y ) arcsin( √ G) -G(1 -G) . (4.16)
We also show convexity in this situation. Indeed we have

ω (α) = -(1 -2α)(α(1 -α)) -3/2 2 ,
from what we deduce that

ω(G) + (G -H y )ω (G) = (G(1 -G)) -3/2 (G(1 -G) - 1 2 (G -H y )(1 -2G) = (G(1 -G)) -3/2 (H y (1 -G)/2 + (1 -H y )G/2) 0 ,
using the same trick that H y is either equal to 0 or 1.

Conclusion

Threshold-weighted and quantile-weighted scores allow to focus on specific areas of interest on the real line or in the forecaster's distribution. We showed that the results of Thorey et al. [START_REF] Thorey | Online learning with the CRPS for ensemble forecasting[END_REF] can be generalized to threshold-weighted scoring rules. Indeed, the effect of the threshold-weighting is equivalent to a simple data transformation. The bias of the ensemble CRPS and the definition of the class CRPS therefore apply to the transformed data. Besides, we demonstrated several results for quantile-weighted scoring rules. After rewriting the score definition, the optimal locations for an ensemble of forecasts were derived. The loss gradients of a quantile-weighted score can be interpreted with a data transformation, but contrarily to the threshold-weighted scoring rules, the data transformation depends on the forecaster's distribution. However, the study of the bias of quantile-weighted scoring rules is left for future research.

Other areas of further work include the study of quantile-weighted scoring rules for extreme value verification, with a focus on the distribution tails. The CRIGN may not be a viable option because this score is not bounded even with bounded observational distribution. A similar study would of interest based on expectiles instead of quantiles. Besides, analytic expressions of the scores for parametric distributions and mixture of parametric distributions would facilitate the testing of these new scores. Evaluation and statistical learning with quantile-weighted scoring rules for real world data sets are shown in Appendix 7.A.

Probabilistic forecasting with observational noise

Bibliographical remarks

We investigate how a forecaster can take into account some knowledge about noise or perturbed observations in the setting of ensemble forecasting. A major issue in this field of investigation is that the noise distribution is generally unknown. The topic of noisy observations is widely studied in the data assimilation community, where a central question is to provide the best trade-off between new observations and background knowledge, both imperfect. When the observation errors are uncorrelated, the observation errors are often scaled by the standard deviation of the noise, such as for the Kalman filter. Besides, taking into account the noise by either perturbing the observation or perturbing the members at the corresponding level of noise can improve the forecasts verification [Sae+04; CT08]. In the case of independent additive noise on the observations, Bowler [START_REF] Bowler | Explicitly accounting for observation error in categorical verification of forecasts[END_REF] and Bowler [START_REF] Bowler | Accounting for the effect of observation errors on verification of MOGREPS[END_REF] demonstrate deconvolution methods to denoise the received observations. They rely on the fact that the distribution of a sum of independent random variables is the convolution of their PDFs.

We emphasize that the setting of noisy observations is seldom investigated in the framework of online learning with experts advice. The losses of each expert are usually considered i.i.d. or noiseless, but rarely with time-varying noise distributions. Still, Yang [START_REF] Yang | Combining forecasting procedures: some theoretical results[END_REF] shows regret bounds on the loss against the observation expectations (and not against the noisy observations) for the square loss. Interestingly, a normalization of the loss by the observation standard deviation is proposed, as in penalized least square regression. Cesa-Bianchi et al. [START_REF] Cesa-Bianchi | Online learning of noisy data[END_REF] focus mainly on expert noise and obtain unbiased estimates of the loss gradients for regret bounds in expectation. We discuss below the interest of the loss expectation in the context of probabilistic forecasting with noisy observations. In online classification with label noise, the work of Ben-David et al. [START_REF] Ben-David | Agnostic Online Learning[END_REF] shows regret bounds in expectation by considering the probability that the label shown to the forecaster is correct or not. Similar approaches are also studied in the batch setting [START_REF] Natarajan | Learning with noisy labels[END_REF], see Frénay and Verleysen [START_REF] Frénay | Classification in the presence of label noise: a survey[END_REF] for a review of classification with label noise.

Context

Say the forecaster receives a distribution or multiple observations instead of a single observation. A basic idea would be to switch from the CRPS to the average CRPS, where the expectation is taken according to the distribution of the observations. Let G be the forecasted CDF, F be the CDF of the observations received by the forecaster and Y be a random variable described by F. The average CRPS is equal to the squared difference between F and G, plus the uncertainty term of F:

E(CRPS(G, Y)) = (F -G) 2 + F(1 -F) .
Two elements should be considered:

(i) the reality, that we note y t , is not a random variable, but a fixed value. The superscript t here stands for "true". We emphasize that the distribution G should reflect the inability of the forecaster to provide perfect estimations of y t , whatever the level of observational noise. The assumed distribution of y t is in fact the most informative distribution that a forecaster could deliver. (ii) the observation y may be corrupted by some noise, from an inaccurate sensor for instance. If this noise is very large, the forecaster may not be willing to depict such a large uncertainty in its forecasts. In any case, the forecaster wishes to predict the reality y t and not the observation y. Consequently, the forecaster may look for a different objective than targeting the observational CDF F to compute the weights of the model mixture G = m M u m G m . We now give thought on strategies for the forecaster wishing both to apply ensemble post-processing methods and evaluate the resulting predictions. We include the aforementioned strategies in our list: -Empirical normalization: scale the members and the observations by an empirical value reflecting the level of noise, such as the standard deviation of the noise. Less weight is attributed to time steps with large inaccuracy. Besides Yang [START_REF] Yang | Combining forecasting procedures: some theoretical results[END_REF], we did not find theoretical guarantee and online learning algorithms taking these normalizing factors into account. For probabilistic forecasting, we suggest to normalize the data by the uncertainty term of the loss function, which is the Gini mean difference of the observations for the CRPS. -Forecasting noisy observations: playing the distribution of the noisy observations F, with the loss (F -G) 2 . This strategy is of limited interest because of the above point (ii). -Perturbing ensemble members: generate perturbed ensemble members according to the observational noise, and work with the perturbed members. However, if the perturbed distribution of the members match the observational distribution, it is not ensured that the distribution of the unperturbed members match the underlying distribution of the truth, as stated by Bröcker and Smith [START_REF] Bröcker | Scoring probabilistic forecasts: The importance of being proper[END_REF] in terms of scoring rules. -Noise deconvolution: generate denoised observations from deconvolution of the observational distribution. -Forecasting reality and not the observations: try to forecast samples from F, and not directly F. The motivation behind is that the reality y t is most likely located in domains with high observational density. We introduce a method to do so in Section 4.2.1.

Generalized least square with the CRPS

Let M forecasts be described by the CDFs G m , and Y CDF be the space of the CDFs. In common data assimilation formulations, the observation state contains several scalar observations, while our observation state describes the occurrence of threshold exceeding of only one observation. An observation y is described by the unit step function H y ∈ Y CDF centered on y. The CRPS between G = m M u m G m and y defined by (H y -G) 2 can be written with the inner product on functions:

(H y -G) (H y -G) = (H y -G) 2 .
(4.17)

We wish to obtain the best estimate of u * defined by the model

H y = m M u * m G m + e o , (4.18)
where e o is the observational error. We emphasize that in this model, the first objective of the forecaster is to be close to the CDF of a unit step function. In this sense, the forecaster targets samples from the distribution described by F and not directly F.

The assumption of model unbiasedness forces the equality m M u * m G m = F because E(H y ) = F. Consequently, e o = H y -F, and the covariance matrix of the observational errors R = E e o (e o ) is fully determined by F, whose knowledge is required.

We now choose a rule to find a good estimate of u * , in the form of a minimization solution. Let R -1 be the matrix inverse of R, then a good estimate of u is given by the minimizer of the generalized least squares (H y -G) R -1 (H y -G), as derived by [START_REF] Aitken | On Least Squares and Linear Combination of Observations[END_REF]. The idea is to shift the minimization problem in a space with uncorrelated, homoscedastic errors, and apply ordinary least square (OLS) in this space. Therefore, the theoretical guarantees of OLS apply in the transformed space. For example, the OLS estimator is unbiased and efficient. The solution of the OLS minimization has the minimum variance compared to other unbiased linear estimators. In data assimilation, the matrix R is commonly assumed to be diagonal, i.e. the errors are uncorrelated. This assumption does not hold in our case, see below the formula of R.

As we show below, this new method allows to target the distribution of the observations, while being quite different from minimizing the average CRPS. We relate the score expectation to the value of Pearson's χ 2 test statistic. We inform the reader that many open questions remain at the end of this chapter, which is present in this thesis because our preliminary theoretical results look promising.

Observational error covariance matrix R

Now we derive formulas for R and R -1 . In a continuous observation space, R is a symmetrical covariance operator such that

R(x, z) = E Y [(H Y (x) -F(x))(H Y (z) -F(z))] = E Y [H Y (x)H Y (z)] -F(x)F(z) = E Y [H Y (min(x, z))] -F(x)F(z) = F(min(x, z)) -F(x)F(z) = F(min(x, z))(1 -F(max(x, z))) , using E Y (H Y ) = F.
We note that the quantity R(x, z) reaches a maximum for F(x) = F(z) = 0.5 and that the diagonal terms R(x, x) are higher than the off-diagonal terms R(x, z) for any x and z = x.

Interestingly, the covariance operator is related to the Gini mean difference and to the variance of the random variable Y described by F:

Tr(R) = R(x, x)dx = F -F 2 = 1 2 E Y,Y (|Y -Y |) , (4.19) and R(x, z)dxdz = 1 2 E Y,Y ((Y -Y ) 2 ) = E Y∼F ((Y -Ȳ) 2 ) , (4.20)
where Ȳ = E(Y). Equation 4.20 was found in González Abril et al. [START_REF] González Abril | The similarity between the square of the coefficient of variation and the Gini index of a general random variable[END_REF].

For the sake of simplicity, we now consider a grid z 1 < z 2 < ... < z S discretizing the real line, with 0

< F(z 1 ) < • • • < F(z j ) < • • • < F(z S ) < 1.
Working in this simplified discrete case where Y CDF has a finite dimension allows us to derive expressions for both R and R -1 . A natural choice is to pick regularly spaced quantiles verifying F(z j+1 ) -F(z j ) = 1/(S + 1). We describe the symmetrical matrix R by its components:

R i,j = min(F i , F j )(1 -max(F i , F j )) , (4.21)
with the notation F j = F(z j ). The said inverse R -1 is tridiagonal with the following components on the ith line:

     0 . . . 0 -1 F i -F i-1 i-1 1 F i -F i-1 + 1 F i+1 -F i i -1 F i+1 -F i i+1 0 . . . 0     
with a slight difference for the first and the last line, respectively:

F 2 F 1 (F 2 -F 1 ) -1 F 2 -F 1 0 . . . 0 , and 0 . . . 0 -1 F S -F S-1 1 -F S-1 (1 -F S )(F S -F S-1 )
.

The general expression of the diagonal term is valid for all indices i with the notation, F 0 = 0, F S+1 = 1, z 0 = -∞, and z S+1 = +∞. Appendix 4.A exhibits a derivation of the general term of R -1 . The term F j+1 -F j is the probability of the event "Y ∈ B j ", where B j is the bin [z j , z j+1 ]. The presence of such terms in R -1 suggests that our penalized least square procedure is related to the occurrence of the events "Y ∈ B j ".

Examples of R and R -1 are shown in Figure 4.7 in the cases of uniform, piecewise uniform and Gaussian distributions for linearly spaced thresholds. We see that R is determined by F, while R -1 is determined by the slope of F. Steep gradients of F are related to low values in R -1 , while moderate gradients are related to high values in R -1 . The example of the Gaussian distribution shows the huge impact of very low gradients in the tails of the distribution. The case of the uniform distribution also illustrates the situation of linearly spaced quantiles (instead of thresholds), where equal probability is given to each bin (F(z j+1 ) -F(z j ) = 1/(S + 1)).

A new loss related in expectation to the χ 2 test

Now we compute the loss R -1 (u) = (e) R -1 (e), where e is a difference of discretized CDFs verifying e 0 = e S+1 = 0. We have

R -1 (u) = S j=1 e j e j -e j+1 F j+1 -F j + e j -e j-1 F j -F j-1 = S j=1 e j e j -e j+1 F j+1 -F j + S j=1
e j+1 e j+1 -e j F j+1 -

F j + e 2 1 F 1 = S j=1 (e j+1 -e j ) 2 F j+1 -F j + e 2 1 F 1 = S j=0
(e j+1 -e j ) 2 F j+1 -F j .

In the last line, the term j = 0 is

e 2 1
F 1 and the term j = S is e 2

S

(1-F S ) . We emphasize that this expression is valid for any e verifying e 0 = e S+1 = 0, hence the errors e i may be defined by either e i = H y (z i ) -G(z i ) or e i = F(z i ) -G(z i ). When F describes the climatological distribution and an observation y is received, the forecaster may be more willing to compare G with H y than with F.

At first look, the new loss could seem to be closely related to the CRPS, up to normalization factors F j+1 -F j . A normalized Ranked Probability Score (RPS) may be obtained by taking only the diagonal terms of R -1 , giving S j=1 (e j ) 2 F j (1-F j ) . However, the loss R -1 (u) is minimized when the difference e j+1 -e j of errors is equal to the difference F j+1 -F j for each j. We explain below why this is a major difference with the CRPS.

We restrict our study to the case of ensemble forecasting, where the CDF G E is a sum of unit step functions centered on the members x m . The difference of errors e j+1 -e j is then related to the presence of x m and y in the bin B j . Let 1 B j be the indicator function of B j (1 B j (z) = 1 if z ∈ B j and zero otherwise), and let u

(j) = G E (z j+1 ) -G E (z j ) = m M 1 B j (x m )
u m be the cumulated weights of the members in B j . For example, if x 4 , x 7 and y are in the bin B j and no other members are in this bin, then e j+1 -e j = 1 -u 4 -u 7 . If x 2 is the only member in B j and the observation is not in B j , then e j+1 -e j = -u 2 .

The loss may be rewritten as According to the location of the members and the location of the observation, the loss R -1 (u) simplifies to:

R -1 (u) = S j=0 1 B j (y) -u (j) 2 F j+1 -F j = S j=0 1 B j (y) 1 -u (j) 2 + (1 -1 B j (y)) u (j) 2
R -1 (u) = 1 -u (j * ) 2 F(z j * +1 ) -F(z j * ) + S j=0 j =j * u (j) 2 F(z j+1 ) -F(z j ) , (4.23)
where the bin B j * contains the observation. In terms of weight optimization, two different situations may occur:

-Case 1: the bin B j * contains the observation and at least one member. The optimal weights verify u (j * ) = 1 and u (j =j * ) = 0. In other words, the members outside of the correct bin should receive a null weight. -Case 2: the bin B j * contains the observation, but does not contain any member (u (j * ) = 0). The loss is minimized when u (j) = F(z j+1 ) -F(z j ) for all j such that the bin B j contains at least one member. Interestingly, the optimal weights do not depend on the location of the observation. To summarize the above cases, the members should try to be in the correct bin, but if no member is good enough to be in the correct bin, then the weights should be in agreement with F. We emphasize that the penalization of the members lying outside of the correct bin is not sensitive to the distance to the observation, which is in sharp contrast with the CRPS. Besides, the normalizing factors scale the loss related to u (j) with the probability that the observation is in B j . Interestingly, the tridiagonal structure [-1, 2, -1] of R -1 recalls two differentiations, one for each e in R -1 (u) one may say. We see that the generalized least square procedure shifts the problem from the CDF point of view to the PDF point of view.

We now check that the score is proper. Strict propriety is not obtained because of the discretization, but we show that the score is minimized by a CDF G verifying G(z j ) = F(z j ). Indeed, by taking the expectation over the observational distribution:

E (e) R -1 (e) = S j=0 (F j+1 -F j ) 1 -u (j) 2 + (1 -(F j+1 -F j )) u (j) 2 F j+1 -F j = S j=0 F j+1 -F j -u (j) 2 F j+1 -F j + (F j+1 -F j ) (1 -(F j+1 -F j )) F j+1 -F j = S j=0 F j+1 -F j -u (j) 2 F j+1 -F j + S .
The proof is complete since the optimal distribution G verifies G(z j+1 ) -G(z j ) = F(z j+1 ) -F(z j ) and G(z 0 ) = F(z 0 ) = 0. The uncertainty term of the score is equal to S, which depends only on the grid resolution, and not on the observational distribution. Besides, the divergence term of E (e) R -1 (e) is equal to Pearson's χ 2 test-statistic. Connections between χ 2 statistics and generalized least squares are not new [START_REF] Rao | Karl Pearson chi-square test the dawn of statistical inference[END_REF].

To the best of our knowledge, the approach of generalized least squares applied to discretized CDFs is innovative. The question of score bias is addressed in Appendix 4.A.

Discussion and further work

Assessing the flatness of the rank histogram is possible with the χ 2 statistic [And96], even though such verification tool may not be adequate because it treats all bins equally and forgets the rank structure [START_REF] Elmore | Alternatives to the chi-square test for evaluating rank histograms from ensemble forecasts[END_REF]. In any case, minimizing the new loss R -1 (u) may be seen as a rank histogram optimization procedure.

The χ 2 statistic is quite celebrated. Providing new insights on this statistic therefore seemed of interest to us. Further work could investigate connections with Cramer-von Mises criterion (G -F) 2 dF and Anderson-Darling statistic (G-F) 2 dF F(1-F) . Indeed, they appear to be related to simplifications of the χ 2 statistic by not taking R -1 into account or only partially with the diagonal only.

The distribution of the observations F is in general unknown. We emphasize that the lack of knowledge of F does not change the score non-strict propriety. Indeed, if the loss is built with the CDF F guess and the observations follow the CDF F, the scoring rule is still minimized in average if the forecaster's CDF follows F and not F guess . Using an incorrect F guess may not impair the scoring rule.

It would be interesting to test this new loss, for both learning (in online learning for example) and evaluating probabilistic forecasts. For real world data sets, the first task would be to determine the observational distribution. The climatological distribution of the observations could be a good starting point. Secondly, the forecaster should determine the appropriate thresholds z j or equivalently the bins B j . A rule of thumb for Pearson's test states that each bin should contain at least 5 observations. Photovoltaic power data sets may not be a good starting point, because observational uncertainty is seldom taken into account contrarily to other variables such as geopotential heights or temperature [START_REF] Saetra | Effects of Observation Errors on the Statistics for Ensemble Spread and Reliability[END_REF][START_REF] Candille | Impact of observational error on the validation of ensemble prediction systems[END_REF]. Further work may investigate the integration of such climatological knowledge in the loss function, and up to which level the forecaster may benefit from the integration of this knowledge.

We conjecture that theoretical results may possibly be obtained by using this new loss in sequential aggregation, for example showing a reduced variance of the weights due to the knowledge of observation uncertainty.

First we compute (RR -1 ) i,j for 1 < j < i S,

(RR -1 ) i,j = R i,j-1 R -1 j-1,j + R i,j R -1 j,j + R i,j+1 R -1 j+1,j = -F j-1 (1 -F i ) F j -F j-1 + F j (1 -F i )(F j+1 -F j-1 ) (F j+1 -F j )(F j -F j-1 ) - F j+1 (1 -F i ) F j+1 -F j = 1 -F i (F j+1 -F j )(F j -F j-1 ) × (-F j-1 (F j+1 -F j ) + F j (F j+1 -F j-1 ) -F j+1 (F j -F j-1 )) = 0 ,
and for 1 i < j < S,

(RR -1 ) i,j = R i,j-1 R -1 j-1,j + R i,j R -1 j,j + R i,j+1 R -1 j+1,j = -F i (1 -F j-1 ) F j -F j-1 + F i (1 -F j )(F j+1 -F j-1 ) (F j+1 -F j )(F j -F j-1 ) - F i (1 -F j+1 ) F j+1 -F j = F i (F j+1 -F j )(F j -F j-1 ) × (-(1 -F j-1 )(F j+1 -F j ) + (1 -F j )(F j+1 -F j-1 ) -(1 -F j+1 )(F j -F j-1 )) = 0 .
The above derivations are also valid for j = 1 with F 0 = 0 and for j = S with F S+1 = 1. Indeed

(RR -1 ) i,1 = R i,1 R -1 1,1 + R i,2 R -1 2,1 = F 1 (1 -F i )F 2 (F 2 -F 1 )F 1 - F 2 (1 -F i ) F 2 -F 1 = 0 , and (RR -1 ) i,S = R i,S-1 R -1 S-1,S + R i,S R -1 S,S = -F i (1 -F S-1 ) F S -F S-1 + F i (1 -F S )(1 -F S-1 ) (1 -F S )(F S -F S-1 ) = 0 .
Second we compute (RR -1 ) i,i for 1 < i < S. We have

(RR -1 ) i,i = R i,i-1 R -1 i-1,i + R i,i R -1 i,i + R i,i+1 R -1 i+1,i = -F i-1 (1 -F i ) F i -F i-1 + F i (1 -F i )(F i+1 -F i-1 ) (F i+1 -F i )(F i -F i-1 ) - F i (1 -F i+1 ) F i+1 -F i = 1 (F i+1 -F i )(F i -F i-1 ) × [-F i-1 (1 -F i )(F i+1 -F i ) + F i (1 -F i )(F i+1 -F i-1 ) -F i (1 -F i+1 )(F i -F i-1 )] = 1 , using F i (1 -F i-1 ) -F i-1 (1 -F i ) = F i -F i-1 .
The terms (RR -1 ) 1,1 and (RR -1 ) S,S are also equal to one:

(RR -1 ) 1,1 = R 1,1 R -1 1,1 + R 1,2 R -1 2,1 = F 1 (1 -F 1 )F 2 (F 2 -F 1 )F 1 - F 1 (1 -F 2 ) F 2 -F 1 = 1 , and (RR -1 ) S,S = R S,S-1 R -1 S-1,S + R S,S R -1 S,S = -F S-1 (1 -F S ) F S -F S-1 + F S (1 -F S )(1 -F S-1 ) (1 -F S )(F S -F S-1 ) = 1 , which completes the proof.
Bias of (e) R -1 (e).

Let the members X m be independent random samples described by the CDF G m , and let P (j) m be the probability that the mth member is in the bin B j . The derivations for the ensemble CRPS bias of Section 3.1.4 were achieved for the events X m θ, while in the current case, the events are of the type z j X m z j+1 . We show that similar derivations can be produced to show that (e) R -1 (e) is a biased score.

We recall Equation 4.22:

R -1 (u) = S j=0 1 B j (y) -u (j) 2 F j+1 -F j .
By taking the expectation over the members,

E u (j) = E m u m 1 B j (X m ) = m u m P (j) m = G(z j+1 ) -G(z j ) .
Besides,

E u (j) 2 = E   m,k u m u k 1 B j (X m )1 B j (X k )   = m =k u m u k P (j) m P (j) k + m u 2 m P (j) m = E u (j) 2 + m u 2 m P (j) m -P (j) m 2 .

Note that P (j)

m ∈ {0, 1} if the mth member lies always in the same bin. The lesser the variability of each member, the lesser the bias (like for the ensemble CRPS bias).

Consequently, each term E u (j) 2 generates an additional m u 2 m P (j) m -P

m -P (j) m 2 F j+1 -F j .
Unbiased estimates of the score can be delivered by counteracting the bias. An interpretation of the bias counteraction with classes of i.i.d. members is left for future research.

Application of online CRPS learning to probabilistic PV power forecasting

We provide probabilistic forecasts of photovoltaic (PV) production, for several PV plants located in France up to 6 days of lead time, with a 30-min timestep. First, we derive multiple forecasts from numerical weather predictions (ECMWF and Météo France), including ensemble forecasts. Second, our parameter-free online learning technique generates a weighted combination of the production forecasts for each PV plant. The weights are computed sequentially before each forecast using only past information. Our strategy is to minimize the Continuous Ranked Probability Score (CRPS). We show that our technique provides forecast improvements for both deterministic and probabilistic evaluation tools.

This chapter is a research paper written with Christophe Chaussin and Vivien Mallet, and submitted to International Journal of Forecasting.

Introduction

Improved photovoltaic power integration needs better power forecasts. Forecasters may pursue efforts to improve meteorological models, weather-based power models or statistical post-processing methods. For our part, we focus on the following case: a forecaster, willing to provide probabilistic PV power forecasts, retrieves multiple meteorological forecasts (possibly from various sources). In this general setting, numerous state-of-the-art methods can be tested and combined.

Meteorological forecasts can either be deterministic single forecasts or an ensemble of forecasts, usually at coarser resolution. Inman et al. [START_REF] Inman | Solar forecasting methods for renewable energy integration[END_REF] provide a review of PV forecasting methods with deterministic forecasts. Ensemble forecasting and more generally probabilistic forecasting has been widely covered in the meteorological community [START_REF] Gneiting | Probabilistic forecasting[END_REF]. Only recently, ensemble-based forecasting techniques are tested for PV [START_REF] Zamo | A benchmark of statistical regression methods for shortterm forecasting of photovoltaic electricity production. Part II: Probabilistic forecast of daily production[END_REF], while these techniques are more common for wind and wind power forecasting [START_REF] Ren | Ensemble methods for wind and solar power forecasting-A state-of-the-art review[END_REF].

A recent benchmark of deterministic and probabilistic PV forecasts is analyzed in Sperati et al. [START_REF] Sperati | The "Weather Intelligence for Renewable Energies" Benchmarking Exercise on Short-Term Forecasting of Wind and Solar Power Generation[END_REF], along with classical diagnostic tools. Probabilistic forecasts rely on the estimation of quantiles of the predicted probability density function (PDF). Quantile regression [START_REF] Almeida | PV power forecast using a nonparametric PV model[END_REF] and analogs [Ale+15; HP15] are amongst most popular techniques for quantile estimation in PV. These techniques do not require an ensemble of forecasts as they can rely only on the historical variability of the forecasts and production data.

A forecaster having multiple forecasts hopefully wishes to combine them in an optimal way. Online learning techniques provide rules for combining forecasts, see the monograph Cesa-Bianchi and Lugosi [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF]. The combination rules stemming from online learning depend only on the available past information at each forecast step and come with theoretical performance guarantee under essentially no assumptions (concerning prior weights, underlying stochastic process or distributions). These techniques have been tested for several applications: electricity consumption, ozone concentration, wind and geopotential fields, and solar irradiance [Sto10; MSM09; Mal10; Dev+13; Bau15; Tho+15]. This paper presents application results with our innovative approach [START_REF] Thorey | Online learning with the CRPS for ensemble forecasting[END_REF], whose purpose is to combine multiple forecasters in a linear opinion pool [START_REF] Genest | Allocating the weights in the linear opinion pool[END_REF][START_REF] Geweke | Optimal prediction pools[END_REF]. The originality of our technique is to use combination rules deriving from online learning techniques in order to minimize the CRPS of the weighted empirical distribution function. We stress here the fact that our method provides theoretical guarantee and that it does not rely on distribution assumptions. Besides, the algorithm has a low computational cost and is parameter-free. Our framework is inspired from the work of Gaillard et al. [START_REF] Gaillard | Additive models and robust aggregation for GEF-Com2014 probabilistic electric load and electricity price forecasting[END_REF], which focuses on quantile scoring functions.

Minimizing the CRPS is a common strategy in the meteorological literature to obtain calibrated probabilistic forecasts. However, standard techniques do not offer theoretical guarantees of robustness and usually resort to strong assumptions on the distributions. For example, Bayesian model averaging (BMA) techniques provide a mixture of parametric distributions, usually a Gaussian sum [START_REF] Gneiting | Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation[END_REF] the parameters of a parameterized distribution using characteristics of the ensemble of forecasts [Gne+05; Wil09; TG10]. For instance, a Gaussian distribution is fitted using a linear model between the mean of the distribution and the mean of the forecasts. Besides, likelihood maximization with the logarithm loss is not an appropriate tool in our setting since it fails to produce satisfactory scores for a mixture of Dirac distributions. A discussion on local scores such as the logarithm loss is adressed by Bröcker and Smith [START_REF] Bröcker | Scoring probabilistic forecasts: The importance of being proper[END_REF]. In Section 5.1, we introduce the production data sets and the forecasts from ECMWF and Météo France. We also detail our method to generate PV forecasts from meteorological data. The evaluation tools are described in Section 5.2. Our statistical postprocessing method is explained in Section 5.3. Numerical results and discussions are developed in Section 5.4.

Methods

Production and meteorological data

The production data cover 219 PV power plants in metropolitan France with 21 consecutive months (January 2012 to October 2013). The total power of the plants is referred to as France production. We wish to provide production forecasts for each power plant and for France production. The data are shown as load factor, i.e. scaled by the installed capacity. France production forecasts are the weighted sums of the plant forecasts w.r.t. the installed capacity of each plant.

Forecast data are summarized in Table 5.1 and 5.2. We use data from two meteorological centers (ECMWF and Météo France), both deterministic forecasts and ensembles of forecasts: HRES and ENS for ECMWF, and ARPEGE and PEARP for Météo France [Cou+91; Des+15; Pal+09], up to a lead time of 6 days. Note that the deterministic forecasts are not the unperturbed members of the ensembles of forecasts but different forecasts, with better resolution.

We are interested in predicting the 30-min average power output of the plants. We only show results for the following hours of the day 0600, 0900, 1200, 1500, and 1800 (where 0600 refers to 6:00 UTC) in order to save computation time and to avoid the issue of temporal interpolation of our forecast data. Ensemble forecasts are only solar irradiance forecasts while deterministic forecasts also include total cloud cover and 2-m temperature. The ensemble PEARP is available for longer lead times but only with a time step of 6 (and not 3) hours. Consequently, for our application we restricted the use of PEARP up to 2 days.

Conversion of meteorological forecasts to production forecasts

Our regression technique is inspired from Bacher et al. [START_REF] Bacher | Online short-term solar power forecasting[END_REF] and Lorenz et al. [START_REF] Lorenz | Irradiance forecasting for the power prediction of gridconnected photovoltaic systems[END_REF]. This regression technique has been successfully applied in the benchmark of Sperati et al. [START_REF] Sperati | The "Weather Intelligence for Renewable Energies" Benchmarking Exercise on Short-Term Forecasting of Wind and Solar Power Generation[END_REF], where the technique ranked first in the deterministic PV forecasting competition. For a given deterministic forecast, the following technique is applied for each time of the day and for each power plant independently. The training set ranges from early 2012 to February 2013 (nearly 400 days). The testing set with the remaining days of 2013 is about 240-day long.

First, clear sky indices τ P and τ I are generated from the production P and the solar forecasts I:

τ P = P P cc and τ I = I I cc , ( 5.1)
where clear sky production P cc and clear sky solar radiation I cc are the production and solar radiation in clear sky conditions. The clear sky profiles P cc and I cc are respectively estimated from the production P and the solar forecasts I thanks to quantile regression introduced in Section 5.1.3. The core of the statistical analysis is a linear regression between the production index τ P and the meteorological variables (the clear sky index τ I , the total cloud cover T cc, and the temperature T 2m ). Non-linear dependencies are taken into account by introducing several terms such as the squared clear sky index τ 2 I and cross terms between variables τ I (T 2m -T 2m ). The quantity T 2m -T 2m is the deviation of the temperature T 2m from its local average value T 2m . The linear regression estimates the coefficients a i to produce

τ P = a 0 + a 1 τ I + a 2 τ 2 I + a 3 T cc + a 4 τ I (T 2m -T 2m ) . ( 5.2)
A secondary statistical model is then fitted on the residuals τ P -τ P . The objective of this secondary model is to reduce the seasonal biases and other remaining errors of the first model. We use the elapsed time s from January 1 st of the current year, and the production forecasts τ P as inputs to build

τ P = b 0 τ P + λ b λ sin(λs) + b λ cos(λs) , ( 5.3)
a linear additive model on τ P and trigonometric polynomials of s, see Lorenz et al. [START_REF] Lorenz | Irradiance forecasting for the power prediction of gridconnected photovoltaic systems[END_REF].

The model parameters of the first two steps are set with the forecasts whose lead times are inferior to 24 h. The motivation behind is that the forecasts with short lead times are presumably the most accurate forecasts to fit the main parameters. Still, forecasting long lead times may necessitate slight corrections compared to short lead times, hence we introduce a third step, which takes into account the lead time of the forecasts. The third step is a multiplicative correction λ applied to τ P , such that λ τ P /τ P is equal to 1 on average.

The statistical regression scheme is slightly different for ensemble and deterministic forecasts. For ensemble forecasts, the input variable of the linear regression is simply the solar irradiance of the unperturbed member without other weather variables. The same conversion model is used for all the members of a given ensemble.

Quantile forecasts

For each deterministic production forecast, we build 19 quantile forecasts (of order 5 to 95) for a total of 38 additional forecasts. They are referred to as deterministic quantiles as opposed to the ensemble members. The idea is to train (on the training data set) a pick-up rule that gives production quantile forecasts according to the value of the deterministic production forecast, see [START_REF] Nielsen | Using quantile regression to extend an existing wind power forecasting system with probabilistic forecasts[END_REF]. We follow the idea that we should first precisely estimate the mean of the distribution and only then estimate the quantiles.

Quantile regression uses a piecewise linear asymmetric loss function QS α (x, y) called the quantile score (or pinball loss) of level α [KH01]:

QS α (x, y) = α(y -x) + + (1 -α)(x -y) + , ( 5.4)
where (•) + = max(•, 0). The expectation (over y) of QS α (x, y) is minimized if x correctly estimates the quantile of level α of y.

We apply quantile regressions on the residuals of the deterministic forecast obtained at the end of Section 5.1.2. The inputs are the deterministic forecast and trigonometric polynomials of s, similarly to the seasonal bias reduction of the deterministic forecasts (step 2 in Section 5.1.2). The quantile regressions are carried out independently for each lead time.

The clear sky profiles P cc and I cc are also based on quantile regressions but with only trigonometric polynomials of s as inputs. The chosen levels of quantiles to build clear sky profiles are close to 90%.

Concerning France production deterministic quantiles, they are not set to a weighted sum of quantiles of the plants, but they are determined from the deterministic forecast of France production. In other words, the deterministic forecasts of the plants are summed to generate France production forecast, and this forecast is used to generate quantile forecasts for France production.

At this point, the forecaster has a total of 50 + 34 + 2 × 19 + 2 = 124 forecasts up to the lead time of 48 hours, 90 forecasts up to the lead time of 96 hours, and 50 forecasts up to the lead time of 138 hours.

Linear opinion pools

Let the x m be M forecasts (or members). The unit Cumulative Distribution Function (CDF) H m (x) = H(x -x m ) equals 0 before x m and 1 otherwise. The forecaster's CDF G = m u m H m is designed as a weighted combination of unit step functions. The mth step of G is centered on x m and its height equals the weight u m . The weights u m are non-negative and sum to one (u in P M the simplex of R M ). This weighted CDF is also known as model mixture or linear opinion pool. Using a discrete CDF based on several forecasts allows us to model any CDF without distribution assumption.

The impact of the weights u m are illustrated in Figure 5.1 and 5.2. Two CDFs G and G using the same locations x m are shown in Figure 5.1. The CDF G is built with uniform weights u m = 1/M , while the weights u m of G are not uniform. We show in Figure 5.2 an illustration of probabilistic forecasts in two different cases: with equal weights for all members and with possibly different weights given by our online learning algorithm. A visual inspection indicates that the online learning algorithm provides a better estimation of the median and a larger spread of the distribution. We emphasize that methods involving weighted empirical distribution functions necessitate that the forecasts x m are sufficiently dispersed.

Evaluation

In the following we describe classical diagnostic tools used in Section 5.4, see for example the monograph of Jolliffe and Stephenson [START_REF] Jolliffe | Forecast Verification: A Practitioner's Guide in Atmospheric Science[END_REF] for further references.

We begin by describing the CRPS as it is at the heart of our learning method.

The CRPS

The CRPS is a classical scoring function in meteorology [Her00; CT05]. The CRPS is the generalization over all thresholds of the Brier score [START_REF] Brier | Verification of Forecasts Expressed in Terms of Probability[END_REF]. Let G be the cumulative distribution function of a forecaster describing the i.i.d. random variables X and X , and y be the observation revealed to the forecaster. The CRPS is defined as where H y (x) = H(x -y) is the CDF assigned to y, the unit step function H centered on y. The CRPS reduces to the absolute error for deterministic forecasts.

CRPS(G, y) = (G -H y ) 2 = E(|X -y|) - 1 2 E(|X -X |) , (5.5)
Assuming that y is a random variable, described by the CDF F, the averaged quantity E y (CRPS(G, y)) (on the observation) is minimized only for F = G. This property makes the CRPS a strictly proper scoring rule [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF], and as such it explains why the CRPS is a classical evaluation tool for probabilistic forecasts.

We highlight the fact the CRPS can also be written as a sum of quantile scores [START_REF] Gneiting | Comparing density forecasts using thresholdand quantile-weighted scoring rules[END_REF]:

CRPS(G, y) = 2 1 0 QS α (G -1 (α), y)dα .
(5.6)

The strategies of minimizing the CRPS or minimizing several quantile losses are therefore closely related. For a CDF step function, the corresponding CRPS is computed as:

CRPS M m=1 u m H m , y = M m=1 u m |x m -y| - 1 2 M m,k=1 u m u k |x m -x k | . (5.7)
which is also concisely noted (u) in Section 5.3.

Other diagnostic tools

The scores are all presented only for the test period, usually averaged over time. Besides the CRPS, we also show results for the celebrated RMSE and MAE for which our forecast is the weighted average M m=1 u m x m . The RMSE of the predictions y with respect to the observations y is given by

RMSE = 1 T T t=1 ( y t -y t ) 2 , ( 5.8)
and for the MAE:

MAE = 1 T T t=1 | y t -y t | .
(5.9)

We use daily scores to show the deterioration of the scores with the increasing lead time. To keep the range of the daily score consistent, the daily score is weighted by the average production of the related hour of the day ȳh . For a score S h depending on the lead time h, the daily score

S (d) = h S h × ȳh h ȳh (5.10)
is computed with summation over the available lead times h corresponding to the same daily lead time. Skill scores are useful to compare prediction performance. In this paper, the reference prediction chosen for skill scores is our forecast. Skill scores for a given score S are written

S skill pred = S ref -S pred S ref , ( 5.11)
so that our forecast shows better scores when the skill scores of the other forecasts are negative.

Online learning with the CRPS

Background

Our objective is to produce an optimal combination (of step functions), or more precisely, to minimize the regret due to unavoidable loss w.r.t. the best learning algorithm for a given class of algorithms. A learning algorithm, along with its weight update rule, is used to find weights for each time step using only available past information. In other words, we want to use an update rule that indicates the value of the weights u m,t and relies only on the values of the past forecasts and observations x m,t and y t with t < t.

Online learning algorithms come with a theoretical guarantee of long term performance. The guarantee is often expressed under the form of a regret bound:

sup T t=1 t (u t ) -inf u∈P M T t=1 t (u) o(T ) , (5.12)
where the supremum is taken over all possible values of x m,t and y t . The notation t refers to the CRPS as in Equation 5.7 with a highlighted dependency on the weight.

In the sense of the theoretical guarantee, our algorithm competes against the best combination with weights constant in time, which can be known only at the end of the experiment and is called the oracle. By definition, the oracle has a better score than individual forecasts (MAE of each forecast), and than any subset ensemble with uniform weights. It is common practice in online learning to use linearized losses, by computing the loss gradients w.r.t. the weights. For the CRPS, the loss gradient m,t of the mth forecaster can be written as

m,t = ∂ t ∂u m (u t ) = |x m,t -y t | - M k=1 u k,t |x m,t -x k,t | + y t - M k=1 u k,t x k,t .
(5.13)

The last two terms are identical for all forecasters and appear due to terms 1-M m=1 u m hidden in the expression of the CRPS, see Appendix B of Thorey et al. [START_REF] Thorey | Online learning with the CRPS for ensemble forecasting[END_REF]. The loss gradient is balanced between the distance of x m,t to y t and the weighted distance of x m,t to the ensemble members. A very good member is therefore close to the observation and far from the other members. A neutral member is equally distant to the observation and the other members.

The loss linearization shifts a regret against the best combination as Equation 5.12 to a regret against the best member only, as we now detail (see also Devaine et al. [START_REF] Devaine | Forecasting electricity consumption by aggregating specialized experts[END_REF]). The convexity and the differentiability of t gives

t (u t ) -t (u) (u t -u) ∇ t (u t ) = u t t -u t .
(5.14) for any two vectors u t , u ∈ P M . Summing over time, we get the following regret bound inequalities:

T t=1 t (u t )-inf u∈P T t=1 t (u) = sup u∈P T t=1 t (u t ) -t (u) (5.15) sup u∈P T t=1 u t t -u t (5.16) = T t=1 u t t -min expert k T t=1 k,t .
(5.17)

As a consequence, an algorithm formulated for linear losses M m=1 u m m and coming with theoretical guarantee (on the expression 5.17) may be used with any convex differentiable loss, by applying the algorithm on the gradient losses. A theoretical guarantee for the non linear losses t (u) is then obtained. In other word, knowing a regret bound for expression 5.17 provides a regret bound for expression 5.15.

Example of general algorithm

Initialization: u 1 ;
For each time index t = 1, 2, ..., T 1. get the vector of predictions data x t , update the learning rate of each member η m,t = 1/ 1 + t t =1 (u t tm,t ) 2 update the regret of each member

R m,t = R m,t-1 + u t t -m,t compute the weights u m,t+1 = η m,t (R m,t ) + / η t (R t ) +
Table 5.3 -ML-Poly algorithm, at time t after y t is given. The vectors η t and R t have M coordinates, respectively η m,t and R m,t . The functions (•) + applied to a vector are applied to all the vector's components.

2. compute the forecaster's choice G t with x t and u t , 3. get the verification y t and compute u t+1 , based on the update rule.

The initial weight vector u 1 is arbitrarily set, e.g., to [1/M, . . . , 1/M ] .

ML-Poly

In this article we use a learning algorithm from [GSE14] called ML-Poly for Polynomially weighted averages with multiple learning rates. The algorithm ML-Poly, described in Table 5.3, has no parameters. The algorithm relies on terms u t tm,t that compare the performance of each member to the performance of the weighted ensemble. The learning rate η m,t checks whether a forecaster's performance is in average close to the performance of the weighted forecast, and the regret R m,t quantifies the regret for not having given higher weights to a forecaster. The ideas behind ML-Poly are on the one hand an adaptation of the algorithm Prod of Cesa-Bianchi et al. [START_REF] Cesa-Bianchi | Improved second-order bounds for prediction with expert advice[END_REF] to multiple learning rates, and on the other hand the introduction of the polynomial potential described in Cesa-Bianchi and Lugosi [START_REF] Cesa-Bianchi | Potential-based algorithms in on-line prediction and game theory[END_REF] and giving the terms (R m,t ) + .

The regret bound of ML-Poly is expressed against the best member for the linearized losses. For all sequences of losses m,t ∈ [0, 1], the cumulated loss of ML-Poly is bounded:

T t=1 u t t min 1 m M    T t=1 m,t + M (1 + ln(1 + T )) 1 + T t=1 (u t t -m,t ) 2    .
(5.18) As opposed to the bound of Equation 5.12, the bound of ML-Poly is of second order due to the term T t=1 (u t tm,t ) 2 . The worst case scenario gives a bound O( √ M T ln T ), indicating that even in the worst case, the weighted forecast will perform at least as well as the best forecast. In the case of i.i.d. sequences of losses, the regret bound is practically constant. A detailed analysis of second-order bounds can be found in [START_REF] Gaillard | A Second-order Bound with Excess Losses[END_REF]. Besides, other algorithms showing second order bounds are described in Koolen and Van Erven [KV15], Luo and Schapire [START_REF] Luo | Achieving All with No Parameters: AdaNor-malHedge[END_REF], and Wintenberger [START_REF] Wintenberger | Optimal learning with Bernstein online aggregation[END_REF].

Application

Experiment setup

Local production data may be unfortunately unavailable for given days and plants. In such cases we removed the related data. However we did not modify France production capacity factor to account for local unavailability, because in our opinion, a challenging task for online learning technique is to reduce biases which may be caused by local null production.

The algorithm is run independently for each lead time and production site (including France production). We run the algorithm as if production data is available at the end of each day. For long lead times h where several observations arrive between the delivery of a forecast and the reception of the corresponding observation, we use shifted weights. At time t we compute u t+h to predict y t+h by using the weights u t+h-1 instead of u t in Table 5.3. For example with the shorthand notation u lead time, day and with a lead time of 36 h, the weights u 36 h,d were delivered at d -1 to forecast y 12:00,d . The weights u 36 h,d are updated to u 36 h,d+1 after y 12:00,d-1 is known at the end of d -1. The key point is that the weight update uses u 36 h,d instead of u 36 h,d-1 to check the combination performance against y 12:00,d-1 .

The production forecasts from PEARP and ENS are sorted by rank in order to associate clearly a weight with an ensemble member. As a result, all the members belong to one of the four sorted subensembles, except for the two deterministic forecasts.

We define a climatological reference for diagnostic purposes, called climatology forecast. For time t, we use 2 months of production data centered on t to estimate a so-called climatological mean and 19 quantiles of climatological production. The climatological mean is used for deterministic evaluations (bias, RMSE, MAE) and the quantiles are used for the CRPS. This method produces a rather 'skilled' reference because the climatology is not only evaluated on the training period but on a rolling period.

We define the raw forecast as the forecast with uniform weights. We use this forecast to assess the gain brought by our online learning algorithm.

The results are shown for PV production forecasts only, and not meteorological variables.

Results

In this section we only show the results for the individual plants. The results obtained for France production are quite similar to those obtained for the plants and are shown in Appendix 5.A.

Scores and skill scores

First we show the classical scores RMSE, MAE, CRPS and bias in Figure 5.3 on a daily average, see Equation 5.10. The confidence intervals indicate the variability of the scores obtained for the plants. The scores are shown for our weighted forecast as well as the raw forecast, the ECMWF deterministic forecast and the climatology forecast. Our weighted forecast gets the best scores up to a lead time of 4 days. Note that our forecast has a quite low bias. For days 5 and 6 the ENS members are the only members available in our study, hence a change of slope in the daily scores. Even for a lead time of 6 days, the climatological forecasts is the worst forecast. Therefore numerical weather predictions may be used to forecast plants productions for a lead time of several days at the 30-min timestep. It is noticeable that our regression scheme q q q q 1 2 3 4 5 6 0.05 0.07 0.09 Horizon CRPS q q q q q q 1 2 3 4 5 6 0.07 0.09 0.11 0.13 Horizon MAE q q q q q q 1 2 3 4 5 6 0.10 0.12 0.14 0.16 Horizon RMSE q q q q q q q Raw quantile det ECMWF Weighted

Figure 5.3 -RMSE, MAE, CRPS and bias for the daily scores, for all sites. The results are shown for 3 forecasts: our weighted forecast, the raw forecast (all members with uniform weights), the deterministic forecast of the ECMWF (and its quantiles for the CRPS). The climatology scores are the following : bias = -0.001, CRPS = 0.089, MAE = 0.139, RMSE = 0.167. The confidence intervals are derived from the scores of all sites. q q q q q q 1 2 3 4 5 6 -0.35 -0.25 -0.15 -0.05 Horizon CRPS Skill q q q q q q q Raw quantile det ECMWF quantile det Arpège ENS PEARP minimizes the CRPS and also achieves improvements on the other scores (RMSE, MAE and bias). We tried to identify situations where our algorithm provides a particular improvement over the raw ensemble, but we did not find discriminatory criteria. For example, the installed capacity of the plant or the CRPS of the raw ensemble are not explanatory statistics of the CRPS skill scores of the raw ensemble.

The CRPS skill scores of 5 ensembles (with uniform weights) are shown in Figure 5.4. The skill scores are assessed against our weighted forecast. The 5 ensembles are the 4 subensembles of our complete ensemble and the complete ensemble as well. Our weighted forecast performs better than any of the 5 ensembles. The best ensemble with uniform weights is (in average) the complete ensemble. This may be due to the variety of the forecasts in the complete ensemble. Although the quantile ensemble from HRES (quantile det ECMWF) performs well before 24 hours of lead time, it is beaten by the complete ensemble afterwards and its skill decreases with time. The skill of the ECMWF ensemble (ENS) increases notably with time, from the worst skill for day 1 to a satisfactory skill for day 4.

Diagnostic tools

Improvements are also shown for several other diagnostic tools but only for a lead time of 36 h (1200, D+1) for the sake of brevity. Better results are obtained for shorter lead times and conversely worst results are obtained for longer lead times. By better we mean improvement of our weighted forecast over the raw ensemble.

The spread-skill diagram checks whether the spread of an ensemble (binned into categories) is consistent with the error of the ensemble mean. The squared spread u m (x m -u x) 2 and the square error (u x -y) 2 are averaged in each bin and their rooted square are plotted against each other. The spread and the error should be ideally equal [For+14]. On the graph, the curves should match the first diagonal. The spread-skill diagram of the ensembles of our study is shown on Figure 5.5. We see that q q q q q q q q q q 0.00 0.05 0.10 0.15 0.20 0.25 0.00 0.05 0.10 0.15 0.20 0.25 Binned Spread RMSE q Raw quantile det ECMWF quantile det Arpège ENS PEARP Weighted our weighted forecasts are closer to the first diagonal than any other subensemble with uniform weights. Our weighted forecasts for plants are still under-dispersive, while the correction is better in the case of France production as shown in Appendix 5.A. The weights provided by the online learning algorithm are larger for the outer members of the ensemble, and especially the lowest members. Consequently the spread of the weighted ensemble is larger than the spread of the raw ensemble and the positive bias of the raw ensemble is mitigated. Besides, the ECMWF ensemble shows the lowest spread and the ensemble PEARP presents very large and very small spreads. However, when the ensemble PEARP shows a small spread, the error is quite larger than the spread.

The rank histogram [And96; TVS99; HC97] or Probability Integral Transform (PIT) is built with the values of the CDFs of the forecaster reached by the verifications along an experiment. The ideal rank histogram is flat. The rank histogram of our weighted forecast and the raw ensemble are shown in Figure 5.6. The rank histogram of our weighted forecast is closer to the ideal rank histogram than the rank diagram of the raw ensemble. The raw ensemble is under-dispersive, since it presents a U-shape. This 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Forecast probability Observed frequency q q q q q q q q q q Figure 5.7 -Reliability diagrams for lead times 36 hours for all sites; black circle: the raw ensemble; magenta triangles: our weighted forecast.

is consistent with the results shown on Figure 5.5.

For a given binary event, the reliability diagram checks whether the observed frequency and the forecasted frequency of the event match [START_REF] Atger | Estimation of the reliability of ensemble-based probabilistic forecasts[END_REF][START_REF] Bröcker | Increasing the reliability of reliability diagrams[END_REF]. The forecasted probabilities of the event are binned into categories. The observed frequency of the event for each category is the share of occurence of the event. The ideal reliability diagram shows a curve along the first diagonal. We use the following event "the production level is lower than the average production", where we use the climatological production defined above as local average production. We show the reliability diagrams of our weighted forecast and the raw ensemble in Figure 5.7. We see that our weighted forecast is very well calibrated for event with low probability, but tends to overpredict the occurence of the event when the event is highly likely.

Conclusion

We have applied the algorithm ML-Poly for the minimization of the CRPS, in order to provide probabilistic forecasts. The algorithm does not depend on any parameter or assumptions on distributions such as Gaussianity, and comes with theoretical guarantee of performance. The regret bound ensures our forecast to perform at least as well as the best forecast in the ensemble.

Our case study investigated the PV production of several power plants in France and the total production of the plants. We have shown that our weighted forecast improves on the raw ensemble, which is the best ensemble with uniform weights. Interestingly, we show that CRPS minimization brings improvement on classical scores for the ensemble mean and probabilistic diagnostic tools. Indeed, the forecasting capability measured by classical scores (RMSE, MAE, CRPS and bias) are improved by our online learning algorithm up to a lead time of 4 days. Besides, the online learning algorithm provides a spread correction as shown on the spread-skill diagrams and on the rank histograms. Horizon CRPS q q q q q q 1 2 3 4 5 6 0.035 0.045 0.055 Horizon MAE q q q q q q 1 2 3 4 5 6 0.045 0.055 0.065 0.075 Horizon RMSE q q q q q q q Raw quantile det ECMWF Weighted A visual inspection shows that the raw ensemble tends to underpredict the occurence of low production for a forecasted frequency between 0.3 and 0.7, when the event is likely to occur. Our weighted forecast does not show this tendency and is symmetrical with respect to the first diagonal although not perfectly aligned.

0.02 0.04 0.06 0.08 0.10 0.02 0.04 0.06 0.08 0.10 Binned Spread RMSE q q q q q q q q q Raw quantile det ECMWF quantile det Arpège ENS PEARP Weighted

Figure 5.11 -Spread skill for lead times 33, 36, 39 hours.

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Forecasted Frequency Observed Frequency q q q q q q q q q q 0.05 0.

PV probabilistic forecasts with the AROME high resolution forecasts

The high- Besides Arpège, Météo France also provides for France the local model AROME, designed to forecast severe weather events such as heavy rains in the south of France [START_REF] Seity | The AROME-France Convective-Scale Operational Model[END_REF]. AROME forecasts are operational since the end of 2008. Examples of solar radiation maps are provided in Figure 6.1 for satellite estimations HelioClim, AROME and the ECMWF deterministic forecast called HRES. We clearly see the differences in resolution in the three maps. We recall that HRES data have a 3-h resolution while AROME and HelioClim maps are shown here with a 1-h resolution. The weather forecasting system AROME provides high-resolution forecasts with a 2.5-km spatial resolution for our period of study in 2012-2013. In this chapter, we study how AROME forecasts can be used to provide probabilistic forecasts of PV power. This question is of interest for us for the following reasons. The high-resolution of AROME provides rich spatio-temporal information. How can we exploit it? AROME forecasts are designed for shorter lead times than those of HRES, ENS, Arpège and PEARP as used in Chapter 5. How does the predicting performance of AROME compare with other weather forecasts? Does using AROME forecasts as additional members improve the overall predicting performance? We recall that HRES is the ECMWF deterministic forecast, ENS is the ECMWF ensemble of forecasts, Arpège is a deterministic forecast from Météo France at a coarser spatial resolution than AROME, and PEARP is the ensemble of forecasts of Météo France based on Arpège.

For this study, we chose 10 power plants from the 219 power plants introduced in Chapter 5. The 10 sites were chosen to test the methods on plants with large capacities and to ensure a geographical coverage of metropolitan France. Compared to the study of Chapter 5, only little changes are to be noticed. The study period still runs from January 2012 to October 2013. Data are still shown after normalization by the plant installed capacity and are consequently dimensionless. Climatological estimations are still generated from a rolling period of 2 months. The statistical models providing conversion between weather forecasts and production forecasts are again trained during the period January 2012 to February 2013, and the remaining time of 2013 is our test period. AROME forecasts of the 00 UTC analysis are used. We work at the 30min temporal resolution with clear-sky interpolation of the solar forecasts. Statistical models are built for each lead time (11:00, 11:30, 12:00, ...), which coincides with the hour of the day in this chapter. The main differences with the previous case study are the limited number of 10 sites, the limited lead time of 24 h, and the larger amount of half-hourly models (not only 06:00, 09:00, 12:00, 15:00 and 18:00).

In the following, we introduce how we generate several forecasts thanks to AROME high-resolution. These forecasts are then combined with ML-Poly and CRPS minimization, which is the same online learning technique than previously.

Building an ensembles of forecasts from AROME forecasts

PV models are trained using almost the same model as in Section 5.1 with only solar radiation forecasts of AROME as input. A reference model is built with the following features: (i) clear-sky normalization of PV power and solar radiation, (ii) smoothing of solar radiation within a 100×100 km 2 square, (iii) linear regression between normalized PV production and solar radiation augmented with non-linear transformations such as square root and square, and (iv) multiplicative bias reduction. The main difference with Section 5.1 is that our reference model for AROME does not include seasonal bias reduction, which failed to provide satisfactory results with AROME.

Leveraging the high spatio-temporal resolution

Spatial smoothing. First, we compare our reference forecasts against models using the almost same features but smoothing areas of 50 × 50 km 2 and 25 × 25 km 2 . The predictive performance of the 3 forecasts are compared against each other in Figure 6.2 for all sites and all lead time. The MAEs * of the reference forecast beats the 50 × 50 km 2 forecast by less than 2%, and the 25 × 25 km 2 forecast by more than 4%. Consequently, the smoothed area of 25 × 25 km 2 is too narrow for AROME to provide accurate forecasts. We check in Figure 6.3 the difference between the forecasts time series for three consecutive days with high production variability. Graphically we see that all 3 forecasts are quite close from one another. Besides, for homogeneous forecasts during clear sky or very cloudy days, the 3 spatial averages are almost identical. Strategies are described below to obtain a wider variety in the forecasts.

Nearby data in time and space. Since AROME may forecast the presence of a cloud, but not at the correct time or location, an interesting way to leverage this spatiotemporal information is to somehow translate it. For instance, we use the (normalized) solar radiation forecasts for 11:30 as inputs in the 12:00 model to provide additional forecasts for 12:00. We emphasize that this is not equivalent to simply using 11:30 power forecasts for 12:00. Thrice more forecasts than before are hence available thanks to shifts of respectively +30 min and -30 min. Figure 6.4 shows the exact time and q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.04 0.08 0.12 0.00 0.04 0.08 0.12 mae ref 100 km mae ref 25 km q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.04 0.08 0.12 0.00 -Cross validation with training blocks of 5 days spaced with buffer blocks of 3 days for a total of 8 training subsets.

A total amount of 336 forecasts are generated with all possible configurations described above.

First sequential aggregation results with AROME meteorological experts

For each lead time, the 336 forecasts are combined with ML-Poly in a similar way as in Section 5.3.

A graphical check with time series. The raw ensemble (with uniform weights) and the weighted ensemble are shown with prediction intervals in Figure 6.6 for France production and also for one of the 10 sites. Four typical days are selected to illustrate a wide variety of situations. The algorithm ML-Poly seems to handle quite nicely the large amount of members. We note that France production appears much smoother than the power plant production, even though France production is here the sum of only 10 sites. For the single power plant production forecasts, strong variations may appear from one time step to another, even though time-shifted experts reduce the amplitude of these variations. Besides, the weighted ensemble seems at least to provide a better estimation of the median of the PDF than the raw ensemble. Both raw and weighted ensembles appear to be under-dispersed for these typical days.

Scores of the combination of AROME forecasts. The bias, MAE, RMSE and CRPS of several forecasts are summarized in Table 6.1 for France production. Once again, the bias of our weighted forecast is quite small, below 0.005. Compared to the raw ensemble (with uniform weights), we show score improvements of 11% for the MAE and the RMSE and 12% for the CRPS. Only little gain is achieved thanks to online learning for MAE and RMSE against AROME reference forecast (1% for the MAE and 3% for the RMSE). We emphasize that our method proceeds to CRPS minimization to generate probabilistic forecasts and not MAE or RMSE minimization for point-Pond Raw Climato AROME Ref Bias 0.002 0.013 -0.0 0.0 MAE 0.048 0.054 (-11%) 0.089 (-83%) 0.049 (-1%) RMSE 0.062 0.069 (-11%) 0.110 (-76%) 0.064 (-3%) CRPS 0.036 0.040 (-12%) 0.061 (-70%) Table 6.1 -Forecasts daily scores (bias, MAE, RMSE and CRPS) of France production for our weighted forecast (Pond), the raw ensemble (Raw), the climatological forecast (Climato) and AROME reference forecast (AROME Ref). Score skills against our weighted forecast for France PV power are indicated in parentheses. We emphasize that a skill of -10% means that the score is 10% worst than the corresponding score of our weighted forecast.

forecasts. Similar results are more finely observed with the bias, MAE and CRPS of each lead time for a single power plant and France production in Figures 6.7, 6.8, and 6.9. We clearly see that AROME reference forecast and our weighted forecast are amongst the best forecasts available. Interestingly, our weighted forecast provides larger improvements against the raw ensemble around noon when the production and level of errors are at their highest level.

Other probabilistic verification tools are shown below after additional members are included in the ensemble.

France production -0.005 0.010 Horizon (hours) bias q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 4 5 6 7 8 9 10.5 12 13.5 15 16.5 18 19.5 21

q Pond Raw Climato Arome Ref
Power plant production -0.03 -0.01 0.01 Horizon (hours) bias q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 4 5 6 7 8 9 10.5 12 13.5 15 16.5 18 19.5 21 France production 0.00 0.04 0.08 Horizon (hours) mae q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 4 5 6 7 8 9 10.5 12 13.5 15 16.5 18 19.5 21

q Pond Raw Climato Arome Ref
q Pond Raw Climato Arome Ref
Power plant production 0.00 0.10 0.20 Horizon (hours) mae q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 4 5 6 7 8 9 10.5 12 13.5 15 16.5 18 19.5 21 France production 0.00 0.04 Horizon (hours) crps q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 4 5 6 7 8 9 10.5 12 13.5 15 16.5 18 19.5 21

q Pond Raw Climato Arome Ref
q Pond Raw Climato
Power plant production 0.00 0.04 0.08 0.12 Horizon (hours) crps q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 4 5 6 7 8 9 10.5 12 13.5 15 16.5 18 19.5 21 q Pond Raw Climato Figure 6.9 -CRPS for each half-hour of our weighted forecast (Pond), the forecast with uniform weights (Raw), the climatological forecast (Climato) for one power plant and for France production.

Adding rolling quantiles experts

Until this point, the experts of our ensemble of forecasts are mainly built using spatio-temporal information of AROME and several parameterizations of our statistical model. The under-dispersion of the raw and the weighted ensemble encourages us to improve our set of experts, by including new experts. This is why we include 2 × 10 quantile experts in this section using quantile regression (10 experts) and quantile random forests (10 experts). These new experts are learned with rolling training periods of 90 days in order to take advantage of the recent available information. We recall that the previous models are static, trained on 13 months. We target the quantiles α ∈ {0.05, 0.15, . . . , 0.95} which are optimal for ensemble CRPS minimization [START_REF] Bröcker | Evaluating raw ensembles with the continuous ranked probability score[END_REF].

The new quantile experts are learned by using the best former experts as input variables. The former experts are ranked according to their average weight given by ML-Poly during a period of 20 days in January 2013, corresponding to the end of the training period. This criterion is chosen to define the added value of the expert for the ensemble, and to realize a model selection. The 5 experts with highest weights are the input variables to train quantile experts. The motivation behind is that we want to improve upon the already good experts, instead of going again through the statistical modeling process with the meteorological input variables. Also, the members with a very low weight are assumed not to help the weighted forecast. Unnecessary experts with very low weight are hence removed from the ensemble. This concerns about 40% of the experts whose average weight is below 0.001. The daily CRPS of the raw ensemble is improved by removing the least weighted members (6% gain for France production and 4% gain in average for all sites). Further work may investigate the use of online learning algorithms to build additional quantile experts as in Gaillard et al. [START_REF] Gaillard | Additive models and robust aggregation for GEF-Com2014 probabilistic electric load and electricity price forecasting[END_REF].

Lasso-penalized quantile regression. Ten new experts are trained with Lassopenalized quantile regression, inducing sparsity in the model thanks to the L1-norm regularization. In practice, a quantile expert of level α is obtained as a combination (a α t ) x best t of the best members selected above and nonlinear transformations of these experts (square, square root and inverse-logit). The parameters a α t are found through minimization of a cost function:

a α d,h = arg min b   d-1 d =d-90 QS α (b x best d ,h , y d ,h ) + λ||b|| 1   . (6.1)
The notation QS α (x, y) refers to the quantile score of level α between the prediction x and the observation y, see Equation 5.4. We adopt here the notation "(day=d, hour=h)" instead of the regular "time=t" to emphasize that quantile experts are learned independently for each half-hour of the day. The parameter λ of the regularization is taken with the default value 0.50. Quantile random forest. A major difference between random forests and quantile regression is that regression trees resort exclusively to analogue searches, and do not estimate linear or nonlinear effects of the input variables. Random forests [START_REF] Breiman | Random forests[END_REF] introduce randomization compared to regular binary regression trees. A large amount of regression trees are built using random training subsets (500 for example). Besides, only a random subset of the input variables are used at each node tree to find the splitting rule. Given new inputs, the output is commonly the average over the trees of the average value of the observations attached to the final leaf of each tree. Quantile random forest [START_REF] Meinshausen | Quantile regression forests[END_REF] enables to provide estimation of quantiles and not only the mean. To do so, the past observations of the final leaves are gathered to construct a CDF, from which a quantile is deduced.

Sequential aggregation results with AROME

statistically calibrated experts

Improvements with rolling quantile experts

We run the same experiment as described in Section 6.1.2, but with changes in the set of experts. The previous weighted forecast is referred to as Pond.0, while the new weighted forecast is referred to as Pond.1. The new ensemble comprises 3 subsets: the remaining forecasts from AROME, quantile regression experts and quantile random forest.

Half-hourly scores. The CRPS of each half-hour of the new ensemble and its subsets are shown in Figure 6.10. We see that Pond.1 gives better scores than its subsets, and quantile regression experts constitute the second best expert subset.

Average CRPS skill scores. The average CRPS skills against Pond.1 are summarized in Figure 6.11. All subsets of the ensemble taken with uniform weights are beaten by the weighted forecast Pond.1, by 5% for the quantile regression experts, 7% for the raw ensemble, and above 11% for the quantile random forest experts and remaining AROME experts. Besides, the score is improved by 7% compared to the previous France production 0.00 0.04 Horizon (hours) crps q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 4 5 6 7 8 9 10.5 12 13.5 15 16.5 18 19.5 21 q Pond Raw Climato Ens.arome Forest Quantile.reg

Power plant production 0.00 0.04 0.08 0.12 Horizon (hours) crps q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 4 5 6 7 8 9 10.5 12 13.5 15 16.5 18 19.5 21 q Pond Raw Climato Ens.arome Forest Quantile.reg Figure 6.10 -For the ensemble with rolling quantile experts, CRPS for each half-hour of our weighted forecast (Pond.1), the forecast with uniform weights (Raw.1), the climatological forecast (Climato), selected AROME forecasts (Ens.arome.1) quantile regression experts (Quantile.Reg), and quantile random forest experts (Forest) for one power plant and for France production.

q q q q q q q q q q q q q q q q q q q q q q q q q -0.15 -0.05 0.05 0.15 Hour Skill mae 4 6 8 11 14 17 20 q q q q q q q q q q q q q q q q -0.15 -0.05 0.05 0.15 Hour Skill rmse 4 6 8 11 14 17 20 q q q q q q q q q q q q q q q q q -0.15 France production for midday hours. q q q q q q q q 0.00 0.05 0.10 0.15 0.20 0.00 0.05 0.10 0.15 0.20 Binned Spread RMSE q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 6.15 -Spread skill of Pond.0 (red) and Pond.1 (black) for each site, computed for the 10 sites and France production for midday hours.

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Forecasted Frequency Observed Frequency (a) Reliability diagrams of Pond.0. q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Forecasted Frequency Observed Frequency q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q (b) Reliability diagrams of Pond.1. As a conclusion, we used the high-resolution of AROME forecasts to build a large amount of forecasts and deliver probabilistic forecasts. We demonstrated that adding rolling quantile members in the ensemble set improves the quality of the weighted forecasts, mainly for the CRPS thanks to adjustments of the ensemble spread. The online learning algorithm proved its efficiency for the tasks of member selection and improving the calibration of probabilistic forecasts. Further work may investigate the best ensemble of forecasts that one may build with AROME to provide perfectly calibrated forecasts.

Comparison of AROME with other forecasts from Météo France and ECMWF

We now combine the production forecasts introduced in the previous Chapters 5 and 6, from the weather forecasts Arpège, HRES, PEARP, ENS on the one hand and AROME on the other hand. We compare four probabilistic forecasts:

-Pond.smooth derived from Arpège, HRES, PEARP, ENS and quantile forecasts (at a coarser resolution than AROME) in Chapter 5.

-Pond.1 derived from AROME and quantile forecasts in the current Chapter 6.

-Pond.all whose set of experts is the union of the sets of expert of Pond.1 and Pond.all.

-Raw.all the uniformly weighted ensemble of Pond.all.

We emphasize that the algorithm ML-Poly with CRPS gradients is again used to generate the weights of Pond.all. The study is restricted to the 10 PV power plants selected above and the hour 12:00 UTC. Differences in previously shown score are due to the selected hour of verification. q q q q q q q q q q 2 4 6 8 10 -0.02 0.00 0.02 0.04 0.06 site label bias q q q q q q q q q q 2 4 6 8 10 0.08 0.10 0.12 0.14 site label mae q q q q q q q q q q 2 4 6 8 10 0.12 0.14 Comparison of deterministic forecasts. First, we compare the production forecasts derived from AROME, Arpège, HRES, Raw.all and Pond.all in Figure 6.17. We recall that AROME reference forecast is described at the beginning of Section 6.1. The bias of Pond.all is quite low, which is a frequent feature observed with weighted combinations. Besides, for each power plant, the lowest MAE and RMSE is almost always reached by Pond.all. The production forecast from Arpège shows the worst performance in average especially for difficult situations, i.e. for sites with high MAE and RMSE. The sites labeled 1 and 2 correspond to easy situations for Pond.all, especially since AROME forecasts show satisfactory results. The overall weights given to AROME forecasts is indeed higher for these 2 sites than for the other sites (not shown). Score skills against Pond.all are summarized in Table 6.2. The forecast Pond.all shows score improvements of MAE and RMSE above 7.5% compared to HRES, 9% compared to AROME.ref, and 15% compared to Arpège.

Comparison of probabilistic forecasts according to the seasonality and the forecasted level of production. Secondly, we compare the weighted forecasts Pond.1, Pond.smooth and Pond.all. For MAE, RMSE and CRPS, skill scores of Pond.1 and Pond.smooth are worst than -5% against Pond.all, see Table 6.4. The performance of Pond.1 and Pond.smooth vary largely from one power plant to another as shown in Figure 6.18. Indeed, the MAE, RMSE and CRPS of Pond.1 and Pond.smooth show a relative difference higher than 10% for the sites 2 and 4, but Pond.1 beats Pond.smooth at site 2 and Pond.smooth beats Pond.1 at site 4. The scores of Pond.all are almost always lower than those of Pond.1 and Pond.smooth. Hence merging all forecasts improves the weighted forecasts in a wide variety of situations. This statement is supported by the MAE depending on the seasonality and forecasted level of production in Table 6.3. The detailed comparison of Pond.1 and Pond.smooth highlights specific regimes. During the spring period, we see that it is much harder to forecast correctly for the sites 3-10 than sites 1-2. For the sites 1 and 2, a high level of production is more often forecasted and Pond.smooth is the worst forecast. However, a different regime is observed for the other sites 3-10, where Pond.1 provides better results than Pond.smooth for high level of forecasted production between February and May, but not for low level of forecasted production. Interestingly the converse is observed during summer where Pond.smooth provides much better results than Pond.1 for high level of forecasted production. Although Pond.all is seldom the best forecast with a large margin, Pond.all is very often amongst the best forecasts within a 0.003 MAE margin, which validates its robustness. The only exception concerns the sites 1 and 2 for high level of production during late winter and spring, where Pond.1 clearly beats Pond.all.

The online learning algorithm helps once more to improve the rank histogram of the weighted forecast against the raw forecast, as one can see in Figure 6.19. The added value of the online learning algorithm is also noticeable on the score gains above 4% for the RMSE, 5% for the MAE and 8% for the CRPS, when comparing Pond.all and Raw.all.

Best minimal ensemble subset. To conclude this analysis, we study the best minimal ensemble subset that compares favorably against Pond.all. The online learning experiment was run several times with subset ensembles. We excluded the ensembles PEARP and ENS from this analysis, because we want to show the achievable performance of a weighted forecast using only deterministic forecasts and their related quantile forecasts. The following four subsets were tested {AROME, HRES, ARPEGE}, {AROME, HRES}, {AROME, ARPEGE}, {HRES, ARPEGE}. In other words, adding PEARP and ENS members in the subset {AROME, HRES, ARPEGE} generates the full ensemble of the members that was used for Pond.all, or adding PEARP and ENS members in the subset {HRES, ARPEGE} generates the ensemble used for Pond.smooth. The weighted forecasts from the subsets {AROME, HRES, ARPEGE} and {AROME, HRES} provide very similar results than Pond.all. We found a relative difference inferior to 1% in terms of MAE, RMSE and CRPS between these three weighted forecasts. Consequently, the weather forecasts AROME and HRES seem sufficient to predict as accurately as possible PV power production. Indeed, including information from Arpège, PEARP and ENS does not improve the quality of the forecasts. We note again that adding unnecessary members do not degrade the forecast quality of the weighted ensemble. Other weighted forecasts generated without HRES or without AROME show lower prediction skills, but their performance are still quite close to those of Pond.all. Indeed, we find relative score differences below 8% in terms of MAE, RMSE and CRPS.

Discussion and perspectives

Probabilistic PV power forecasts are built with the high-resolution weather forecasts AROME. First, multiple forecasts were generated to take advantage of the rich spatio-temporal information of AROME. Secondly, we added rolling quantile forecasts in our ensemble thanks to quantile regression and quantile random forests. These forecasts were combined with online learning techniques so as to minimize the CRPS. We q q q q q q q q q q 2 4 6 8 10 0.08 0.09 0.10 0.11 0.12 0.13 site label mae q q q q q q q q q q 2 4 6 8 10 0.12 0.13 0.14 0.15 0.16 site label rmse q q q q q q q q q q 2 4 6 8 10 showed that statistically calibrated members provide a help in the calibration compared to purely meteorological members. Indeed, the rolling quantiles greatly help the calibration of the weighted forecast, even though little improvement is brought to the MAE and the RMSE. We also studied the predictive power of the weighted forecasts including Arpège, AROME, HRES, PEARP and ENS data. We showed that using multiple weather forecasts improves the weighted combination. A minimal subset ensemble of HRES and AROME forecasts proved to deliver accurate forecasts. The online learning algorithm proved to be a serious advantage in the forecast calibration. Indeed, the weighted forecast showed improved calibration compared to its raw ensemble in all cases shown above. Further work may investigate the calibration of the probabilistic PV forecast with only AROME in order to avoid forecast peaks. Providing an in-depth analysis of the predictive power of the weighted combination using only HRES and AROME data may be of independent interest. An ensemble built with specialized experts based on meteorological regimes could take benefits from the multiple weather forecasts. This raises the questions of which weather forecast and which quantile should be emphasized and under what circumstances. In an other setting, we also could combine solar radiation forecasts first with HelioClim maps or ground data as observation, and then convert the improved solar forecasts to production forecasts. This is left for future research.

PV probabilistic forecasts with intraday updates for insular systems

Operational forecasts of Corsica and Réunion include intraday updates of the forecasts with satellite-derived forecasts. These forecasts do not use intraday PV power observations, but only satellite information. We wish to go one step further to improve insular systems PV power forecasts by: (i) generating new predictions with intraday PV power observations, (ii) combining these newly built forecasts with day-to-day predictions and satellite predictions. We show that these new predictions greatly help to improve the quality of the forecasts, and we estimate the accuracy improvement brought by satellite data.

The insular electric systems of Corsica and Réunion are under the responsability of EDF-SEI. According to the 2016 yearly report on insular electricity production from EDF * , the install PV capacity in Réunion and Corsica are respectively of 187 MWc and 117 MWc, covering respectively 8.50% and 6.60% of the total production of each island. No interconnection is available for Réunion while Corsica is only connected to Italy. Consequently, the stability of these electric systems may largely suffer from inaccurate PV power forecasts. Réunion and Corsica are two different situations. It is easier to forecast for Corsica, which receives from many clear sky days, while Réunion weather is subject to numerous micro-climates and tropical atmospheric dynamics [START_REF] Lauret | Solar Forecasting in a Challenging Insular Context[END_REF][START_REF] Badosa | Reliability of day-ahead solar irradiance forecasts on Reunion Island depending on synoptic wind and humidity conditions[END_REF].

In this chapter, we intend to build seamless probabilistic forecasts using operational satellite forecasts and day-to-day forecasts for the total production of each island. Forecasting solar radiation with satellite images is not a new idea [START_REF] Hammer | Short-term forecasting of solar radiation: a statistical approach using satellite data[END_REF], for example tested with HelioClim data [START_REF] Dambreville | Very short term forecasting of the Global Horizontal Irradiance using a spatio-temporal autoregressive model[END_REF]. Besides, combinations of satellite and day-today forecasts were found in Lorenz et al. [START_REF] Lorenz | Short term forecasting of solar irradiance by combining satellite data and numerical weather predictions[END_REF], but not for probabilistic forecasting and not with an online learning algorithm. We focus on lead times between 30 min and several hours to bridge the gap between day-to-day forecasts and satellite forecasts, at a 30-min temporal resolution. The operational forecasts may be biased and not correctly calibrated, hence we generate new forecasts using real-time production data and the available forecasts. We emphasize that we rely on a stronger assumption than operational forecasts: the availability of real-time production data. We show that real-time updates make it possible to improve greatly the prediction accuracy compared to dayto-day forecasts. Online CRPS learning is used again to built forecasts combinations. For very short lead times, the challenge is to beat the persistence forecast, which is the projection at time t -∆t of the latest available production data y t-∆t to time t, by taking into account the clear sky evolution:

y persistence (t) = y(t -∆t) × P cc (t) P cc (t -∆t)
,

where P cc (t) is the clear sky production profile for time t.

Intraday PV updates experimental setup

Operational forecasts

Satellite forecasts of PV power are already operational, using data from the satellite Meteosat 10 for Corsica and Meteosat 7 for Réunion (soon replaced by Meteosat 8). The images are available every 15 min at a 3-km resolution for Corsica for both the visible and infrared spectrums. For Réunion, the resolutions of 2.50 km are available for visible images and 5 km for infrared images every 30 min. We now describe how these images are processed to deliver PV power forecasts. The procedure is applied at over 130 production sites in Réunion and over 20 production sites in Corsica. The total production of each island is the sum of the production of these sites.

Albedo images from the visible spectrum are transformed into cloud indices after scaling: cloud index = albedoalbedo min albedo max -albedo min ,

to clear the images from information on the ground state and on the hour of the day. The extreme albedo min and albedo max are found with the 3% and 97% quantiles of the 30 previous days. For morning forecasts up to the solar elevation angle of 5 • , infrared images of brightness temperature are used. Brightness temperature is the temperature of a black body with the same radiation intensity for this bandwidth. Infrared images are not processed to derive a cloud index.

The forecasting procedure resorts only to image processing and statistical learning but no physics. One statistical model for each site and for each month is trained with spline regression between production data and the pixel of the site in real-time images. In the model, the output production data are normalized by a clear sky profile. We emphasize that the same model with visible images is used for several hours of the day and all lead times. PV power forecasts are generated with these statistical models and cloud index forecasts that we now describe. First, a block-matching algorithm is applied to the last two cloud index images to estimate the cloud motion. Secondly, the cloud motion estimation defines a zone of interest in the last image. The average pixel value of the zone of interest gives the forecasted cloud index. The zone of interest is a triangle with one vertex at the site location and a 20 • angle. The triangle direction is the upstream cloud motion and its height is set by the lead time and the cloud motion speed.

We emphasize again that the operational satellite forecasts uses real-time data from satellites but no real-time PV power observations. Production data are not perfectly reliable, because data transmission from production meters to grid operators may fail.

Operational day-to-day production forecasts are also available. They rely on solar radiation, total cloud cover and 2-m temperature HRES forecasts from ECMWF of the 00 UTC analysis. Statistical models are built to deliver a deterministic production forecast using ECMWF data and quantile predictions associated with the deterministic forecast. These statistical models are similar to those of Chapter 5.

In the following, the operational satellite forecast is referred to as "prevsat" and the operational HRES forecast is referred to as the day-to-day forecast or "prevenir".

Building new forecasts with intraday updates

As in Chapter 6 with AROME forecasts, we generate rolling quantile predictions with Lasso-penalized quantile regression and quantile random forests. These quantile forecasts are trained with 5 major input variables:

-the last production forecast from either "prevenir" or "prevsat", -the last production data available, -their corresponding cumulated value over the current day to indicate the average level of production and forecast for this day, -the production data of the same hour of the previous day.

Name Description

Pond weighted forecast Raw raw ensemble with uniform weights prevenirQ quantile day-to-day forecasts qr_prevsat rolling quantile regression with satellite forecasts qrforest_prevsat rolling quantile random forests with satellite forecasts qr_prevenir rolling quantile regression with day-to-day forecasts qrforest_prevenir rolling quantile random forests with day-to-day forecasts Past day production data is especially interesting to predict very low level of production during sunrise and sunset. The statistical models use normalized data by a clear sky production profile, and the cumulated values are normalized by a cumulated clear sky profile. We generate models for each lead time and each half-hour of the day with a rolling training period of 3 months. For quantile regressions, nonlinear transformations (square, square root and inverse logit) of the 5 major input variables are included in the model to reach the total of 20 input variables.

A total of 60 new forecasts are generated with 15 level of quantiles, 2 statistical methods (quantile regression and random forests) and 2 input variable sets (day-to-day and satellite forecasts). Larger training sets are obtained with available data of nearby half-hour of the day. For instance, the (12:00 UTC, 1-h lead time) model is trained with the data initially gathered for the (11:30 UTC, 1-h lead time), (12:00 UTC, 1-h lead time) and (12:30 UTC, 1-h lead time) models.

In the full ensemble of experts, the following 82 forecasts are present:

-1 deterministic satellite forecast "prevsat", -1 deterministic day-to-day production forecast "prevenir" and 19 corresponding (non-rolling) quantiles "prevenirQ", -60 rolling quantiles introduced in this section, -1 persistence forecast "Persistence".

The forecasts are summarized in Table 7.1 with notation for each subensemble. Climatological productions (Climato) are also shown in the figures to illustrate the gain brought by forecasts.

Online learning experiment

The data sets range from May 2014 to June 2015 for Corsica and from September 2015 to September 2016 for Réunion. The online learning experiment starts at the beginning of these periods. The test periods begin 3 months later due to first rolling period. The satellite statistical models are trained in 2010-2011 for Réunion and 2013-2015 for Corsica. The day-to-day statistical models are trained in 2012-2014 for Réunion and from 2012 to June 2015 for Corsica. Our experiment is therefore not fully operational since the training and the testing periods of "prevsat" and "prevenir" overlap for Corsica. We show that the rolling quantile predictions are critical to improve the forecast accuracy. Since the rolling quantile predictions do not have any overlap between training and testing, we think that the conclusions of this chapter are rather robust.

We use the algorithm ML-Poly with CRPS gradients, as in Chapter 6 and 5. The weights were initially updated once per day for each lead time and half-hour of the day. As an example the weight of the mth member with 30-min for the 12:00 UTC PV production of day d that we note here u m,12:00,+30min,d was updated to u m,12:00,+30min,d+1 after y 12:00,d is received. Because the testing periods last less than a year, we use a trick for the weights to see more data. For a fixed lead time, the available data of nearby half-hour the day also produce a weight update. For instance, the observations and the forecasts of 1-h lead time of 11:30 UTC, 12:00 UTC and 12:30 UTC contribute to the update of the (12:00 UTC, 1-h lead time) weights. This setting is still operational because the weight updates can be achieved at the end of each day. The weighted combination of the forecasts is referred to as the weighted forecast.

Results

Time-series, spread and weights

Time-series of the weighted forecasts of Réunion and Corsica are shown in Figure 7.1 and 7.2. The day-to-day forecasts are also shown to illustrate their large spread and their smooth aspect. The error and the spread of the weighted forecasts increase with the lead time, hence more accurate forecasts are delivered for short lead times. Production sharp variations affect the weighted forecasts after the delay corresponding to the forecast lead time. The forecasts sharp variations are due to the rolling quantile forecasts and not to the weight updates, which occur at the end of the day.

The evolution of the probabilistic forecasts spreads in Figure 7.3 clearly indicates the level of uncertainty associated to each subensemble. For short lead times, the spread of rolling quantile forecasts is quite small thanks to the production data within the model inputs. For long lead times of a few hours, the spread of satellite-based forecasts are much larger than the spread of day-to-day forecasts, providing hints on the predictive power of satellite-derived forecasts for such long lead times. We also observe that quantile random forests deliver forecasts with larger spread than quantile regression while the same model inputs were used, and that probabilistic forecasts have a larger spread for Réunion than for Corsica. Besides, the raw ensemble shows quite a large spread, compared to those of our weighted forecast.

The average weight given to each member according to the lead time highlights the predictive power of the subensembles, see Figure 7.4 for Réunion and Figure 7.5 for Corsica. The small spread of the weighted ensemble at short lead times is due to the high weights given to forecasts with small spread, i.e. rolling quantile predictions. The weights given to the day-to-day forecasts increase with lead time, while the weights of the persistence forecast decrease with the lead time. This illustrates a trade-off between daily weather information and production information. The forecasts resulting from quantile random forests often receive lower weights for their extreme quantiles. This The day-to-day forecast (above) does not benefit from intraday updates. q q q q q q 0.5 1 1.5 2 2.5 3 0.03 0.04 0.05 0.06 0.07 0.08 Horizon (hours) Spread q q q q q q q q 0.5 1 1.5 2 2.5 3 3. may be due to the larger spread observed for quantile random forest subensembles. Moreover, satellite-derived forecasts receive higher weights for Corsica than for Réunion. To conclude this analysis, we note that quantile regression forecasts are more prone to receive lower weights than quantile random forest forecasts, while the former exhibit better performance than the latter. We also highlight the difficulty to interpret the average weights of one forecast, because they heavily rely on the behaviour of the other members. For example, 2 median forecasts with similar quality may either share the overall weight that would have been attributed if only one median forecast is present in the ensemble, or one median forecast may prevail over the other one, the latter receiving a very low weight compared to the former.

Probabilistic forecasts performance and calibration

Results for the 30-min lead time. We now focus on the scores obtained for the 30-min lead time, when satellite forecasts are possibly at their best, see Figure 7.6 and 7.7. Rolling quantiles, using real-time production, obviously predict much better than the initial forecasts, both satellite and day-to-day forecasts. This is verified for the rolling subensembles with the CRPS and for the rolling median quantiles for the MAE. We note the difficulty to beat the persistence forecast for this short lead time, while it is part of the rolling median inputs. The large MAE of the satellite forecast "prevsat" in the morning at Réunion (worst than climatology) are due to large biases in the forecasts. Moreover, the satellite forecasts at Corsica are not better than day-today forecasts for midday hours. The fact that a single model is used for all half-hours of the day may explain the poor performance of "prevsat", and supports our motivation to build models for each half-hour. Besides, the high weights given to satellite-derived quantile forecasts at Corsica encourages us not to drop satellite forecasts, but supports the need to improve the initial forecast "prevsat". Our weighted ensemble is amongst the best forecasts, but is beaten in MAE by the persistence forecast at Réunion. Score daily averages. Daily scores are average scores over all half-hours weighted 3-h lead time 0.00 0.06 0.12 Hours mae

                                                                   
3 4 5 6 7 8 9 10 11.5 13 14.5 16 17.5 19 by the mean observation as introduced in Equation 5.10. Our weighted forecast ranks first for almost all lead times and scores. As expected, the quality of the forecasts gets worse with increasing lead time, see Figures 7.8, 7.9, 7.10 and 7.11. The higher difficulty of forecasting for Réunion is seen with the scores of the weighted forecast. The MAE for Réunion lies between 0.029 (30 min) and 0.057 (3 h), while the MAE for Corsica lies between 0.021 (30 min) and 0.046 (4 h), producing a score relative difference above 20%. The satisfactory performance of the persistence forecast is tempered by its poor quality for lead times longer than 90 min. The satellite-derived quantile forecasts "qr_prevsat" are better suited for short lead times, especially for Corsica where they largely beat "qr_prevenir" and "qrforest_prevenir" quantile forecasts deriving from day-to-day forecasts (about 10% better). However, they show worst performance than day-to-day forecasts (with no updates) for lead times longer than 2 h in both Réunion and Corsica. In fact, day-to-day forecasts are amongst the best forecasts for long lead times, because little information is added from observations dating from several hours earlier. Interestingly, quantile regressions "qr_prevenir" and "qr_prevsat" provide respectively better results than the corresponding forecasts with quantile random forests, but they receive lower weights than their random forests counterparts.

The raw ensemble is always amongst best forecasts in terms of CRPS, but not for the MAE especially for short lead times. Consequently, including day-to-day forecasts in the ensemble set seems to degrade the distribution mean, but to improve the distribution spread for short lead times. In other words, the raw ensemble takes benefits from the large ensemble spread of day-to-day forecasts for the CRPS but not for the MAE, which depends only on the ensemble mean.

For Corsica at intermediate lead times of 120 min and 180 min, we observe noticeable effects of model blending since the raw ensemble and the weighted forecasts have better scores than other forecasts. These high gains may be due to the high diversity between the individual forecasts. Although it is difficult to show causality, we highlight the correlation between these two facts through the higher difference in Corsica between rolling quantile forecasts from day-to-day and satellite data. We checked the average absolute difference between the median quantiles of "qr_prevenir" and "qr_prevsat", and found that this average absolute difference is twice higher for Corsica than Réunion for the lead time of 30 min and 30% higher for the lead time of 3 h.

Added value of satellite information. In order to assess the utility of day-to-day and satellite forecasts, we run the online learning experiment on two ensemble subsets, either using satellite information or day-to-day forecasts. The persistence forecast is included in each subset. The first subset comprises "prevsat", "qr_prevsat", "qrfor-est_prevsat" and "Persistence", while the second subset comprises "prevenir", "pre-venirQ", "qr_prevenir", "qrforest_prevenir" and "Persistence". The weighted forecast without satellite information is referred to as "Pond.prevenir" and the weighted forecast without day-to-day forecasts is referred to as "Pond.prevsat". Results are shown in Figure 7.12. For Réunion, we find that satellite forecasts do not bring a significant amount of information since very similar forecasts are obtained with and without them. On the contrary, satellite forecasts are quite useful for Corsica. Indeed, they allow a 10% gain up to 2 h of lead time against the weighted forecast "Pond.prevenir" (without satellite q q q q q q 0.05 0.10 0.15 Horizon (hours) mae q q q q q q 0.5 1 1.5 2 2.5 3 q q q q q q -0.5 q q q q q q q q 0.02 0.06 0.10 0.14 Horizon (hours) mae q q q q q q q q 0.5 1 1.5 2 2.5 3 3.5 4 q q q q q q q q -0.5 -0.3 -0.1 0.1 Horizon (hours) mae skill q q q q q q q 0.5 1 1.5 2 2.5 3 3.5 4 information). Besides for Corsica, the weighted forecast without satellite information "Pond.prevenir" takes advantage of the persistence forecast for short lead times, which demonstrates once again the importance of real-time PV power observations for short lead times. Probabilistic forecasts calibration. The calibration of the weighted forecasts and the raw ensemble are now checked to validate the help brought by the online learning algorithm. Rank histograms are shown in Figure 7.13. We see that the weighted forecasts have much flatter rank histograms, compared to the over-dispersed raw ensemble. This statement is especially true for short lead times. Yet, our weighted forecasts show a light under-dispersion, clearly seen at the outer bars. Including a finer description of the rolling quantiles (61 quantiles instead of 15 quantiles) did not improve the weighted q q q q q q 0.025 0.035 0.045 0.055 Horizon (hours) crps 0.5 1 1.5 2 2.5 3 q q q q q q -0.5 q q q q q q q q 0.02 0.04 0.06 Horizon (hours) crps 0.5 1 1.5 2 2.5 3 3.5 4 q q q q q q q q -0.5 q q q q q q -0.5 -0.3 -0.1 0.1 Horizon (hours) crps skill 0.5 1 1.5 2 2.5 3 q q q q q q q q -0.5 forecast calibration (not shown).

The spread-skill diagrams validate many aforementioned statements from the rank histograms, from the evolution of the spread and from the evolution of the scores with the lead time. We see the strong error increase between the lead times of 30 min and 3 h for several level of errors in Figure 7.14. Here the over-dispersion of the raw ensemble is exposed, especially for short lead times. The light under-dispersion of our weighted forecast concerns small level of errors for Réunion, while an over-dispersion of our weighted forecasts is revealed for large levels of errors at Corsica. The improvements brought by online learning are clearly demonstrated, since the spread-skill diagrams of our weighted forecasts match better the first diagonal. The case of Réunion at long lead times mitigates this statement.

We conclude this analysis with reliability diagrams for the event "production is below the climatological production". They show satisfactory agreement with the first diagonal, especially for long lead times in Figure 7.15. The questionable reliability of our weighted forecast for short lead times may be due to the tighter spreads of the subensembles.

Conclusion and perspectives

Probabilistic PV power forecasts were built for Réunion and Corsica for lead times of 30 min to several hours. They relied on multiple forecasts based on day-to-day forecasts, satellite forecasts and real-time PV power observations. The multiple forecasts were Corsica Réunion q q q q q q q q 0.00 0.05 0.10 0.15 0.00 0.05 0.10 0.15 Binned Spread RMSE q q q q q q q q q q q q q q q q 0.00 0.05 0.10 0.15 0.00 0.05 0.10 0.15 Binned Spread RMSE q q q q q q q q Figure 7.14 -Spread-skill diagram for Corsica and Réunion for the lead times of 30 min (bold line) and 3 h (thin line) for the weighted forecast (black) and the raw ensemble (red). q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Forecasted Frequency Observed Frequency q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Forecasted Frequency Observed Frequency q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Forecasted Frequency Observed Frequency q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Forecasted Frequency Observed Frequency combined with an online learning algorithm to calibrate the weighted forecast, which outperformed almost always the individual forecasts. Consequently, we have shown that sequential aggregation enabled a seamless blending of the multiple sources of information for each lead time. For short lead times, the large improvement of the forecast quality was due to the rolling quantiles estimated with real-time production data. Satellite information proved to be very useful for Corsica, especially below 2 h of lead time. In contrast, we did not identify a gain brought by satellite information for Réunion.

Further work may focus on statistical modeling of satellite data. Indeed, the single statistical model used for all hours and lead times could be refined. Improved persistence forecast with ARMA (autoregressive moving average) models could also be of interest to quantify the gains brought by satellite data. Such study may include the identification of scenarios where satellite information demonstrates its value over time series analysis of production data. Furthermore, probabilistic forecasting of ramp events as already achieved for wind power in Bossavy et al. [START_REF] Bossavy | Forecasting ramps of wind power production with numerical weather prediction ensembles[END_REF] or in the review Gallego-Castillo et al. [START_REF] Gallego-Castillo | A review on the recent history of wind power ramp forecasting[END_REF] would give new indicators to grid operators, going beyond the confidence intervals of each half-hour of the day. Additional studies could also investigate other information sources such as AROME very short term forecasts (updated every hour) or sky cameras data. our weighted forecasts Pond are shown for several quantile-weighted scoring rules in Figure 7.16 for the Réunion data set. Only little variation is noticed in the skill scores of the raw ensemble, whatever the chosen evaluation metrics. This results is an indicator of the robustness of the results obtained with the CRPS. In the skill scores of the subensemble "qr_prevenir", we find a difference up to 10% between the evaluation with the cubic-CRPS and with the CRIGN, indicating that the subensemble provides a better description of the distribution central quantiles than the distribution tails.

A better discretization of the distribution tails of "qr_prevenir" mainly improves the CRIGN skill score. This is shown by increasing the subensemble size to 61 quantiles instead of the initial 15 quantiles. Increasing the subensemble size also decreases the difference in the skill scores between the CRIGN and the cubic-CRPS to 5%. As already mentioned, we did not find any advantage of running sequential aggregation algorithms with 61 rolling quantiles in each sub ensembles. Secondly, we run online learning experiment with gradients of quantile-weighted scores as loss function. We recall that these gradients are introduced in Section 4.1.2. We generate new weighted forecasts based on the optimization of the arcsin-CRPS and the CRIGN. The algorithm ML-Poly is run because it is parameter-free and its weight update scales well with the order of magnitude of the loss. For the CRIGN gradients, the terms 1/(α(1 -α)) are clipped between between 4 and 20. Consequently, quantiles of order below 5% and above 95% are equally weighted by the weighting function ω. Without surprise, we found differences in the weights of each member depending on the loss function, but without significant impact on the weighted forecast. Indeed, weighted forecasts learned with the CRPS, the arcsin-CRPS or the CRIGN show relative score differences below 1% in terms of CRPS, arcsin-CRPS and the CRIGN, as indicated in Table 7.3. Other evaluation tools such as the rank histogram or the reliability diagram did not help us to discriminate the weighted forecasts.

As a conclusion, the results of this Appendix are mainly robustness of the results in favor of the CRPS for both evaluation and learning. Indeed, the comparison of the raw ensemble against the weighted forecast (learned with the CRPS) indicates similar skill scores for various performance metrics. Besides, online learning with the CRPS, the arcsin-CRPS and the CRIGN provided very similar weighted distributions. Further study should investigate a wider range of data sets, concerning for example extreme-value verification with quantile-weighted scores focusing on the distribution tails. Score sensitivity to correct quantile estimation is another direction for future research as well as the impact of the ensemble size.

shows almost the same performance as the weighted forecast including all available power forecasts.

The study of intraday updates for insular systems is undertaken in Chapter 7 for Réunion and Corsica, where intraday satellite-derived forecasts and day-to-day weatherderived forecasts are blended. The major impact of integrating the latest available information is demonstrated. Besides, satellite information proves to enhance the forecast accuracy for Corsica but not for Réunion.

Online learning methods face practical limitations due to data unavailability. While it is always possible to replace missing data with moving averages, other ideas for missing data imputations are described below. In case of missing PV power observations, the average of nearby power plants productions may be used as replacing value. If the forecaster is only interested in the total production, adjusting the total installed capacity is an easy-to-implement solution. In case of missing forecasts, it may be possible to run the online learning experiment since its beginning by including only the available experts. Besides, the setting of sleeping (or specialized) experts is wellstudied, especially for the algorithm ML-Poly. In a few words, the weights are allowed to incorporate prior knowledge through the arbitrary confidences I m,t , which are also integrated in the regret bound. The weight u m,t becomes proportional to I m,t u m,t , and the loss m,t becomes I m,t m,t + (1 -I m,t ) k M u k,t k,t .

In our case studies, we encountered several open questions that we summarize below and leave for future research. Several of them concern the multidimensionality of the output variable.

1. In Chapter 2, it is shown that sequential aggregation improves solar forecasts. A new expert for PV forecasts can therefore be built using improved solar forecasts. What is the gain in PV power performance that is achievable thanks to this new expert, and for which spatio-temporal scales?

2. Can statistical errors from the conversion model and meteorological errors be separated? A first step would be to build a perfect statistical model with weather observations or satellite observations. However, observational noise and representativeness errors prevent from building this perfect model.

3. Better solar forecasts are obtained by running sequential aggregations independently at each grid point. Could spatial patterns be included in this setting? An idea would be to apply sequential aggregation over spatial averages covering specific regions. The definition of such spatial subsets and their evolution over time is however unclear.

4. We focused on the marginal distribution of each time step. Can we extend our framework to include forecast trajectories and probability of ramp events as in Bossavy et al. [START_REF] Bossavy | Forecasting ramps of wind power production with numerical weather prediction ensembles[END_REF] or in the review Gallego-Castillo et al. [START_REF] Gallego-Castillo | A review on the recent history of wind power ramp forecasting[END_REF] for wind power? This may require multivariate experts in the time dimension.

5. The CRPS generalizes to the energy score for multivariate probabilistic forecasts [START_REF] Szekely | E-Statistics: The energy of statistical samples[END_REF][START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF]. Probabilistic multivariate forecasting therefore seems compatible with online learning, provided that multivariate experts are available. We showed that quantile forecasts greatly improve the quality of the forecasts. Could multivariate quantile forecasts be generated? Possible applications include the challenging task of joint forecasting wind, PV and the electricity demand.

6. Say expert k performs significantly better than expert m in situation A k , and conversely say expert k performs significantly worst than expert m in situation A m . How should this prior knowledge be integrated? The forecaster may give prior weights to the experts m and k according to its opinion that situation A k or A m will occur, or the forecaster may build a new expert that is a combination of the experts m and k. This new expert should be closer to expert k in situation A k and closer to expert m in situation A m . Our opinion is that the latter setting allows for a greater flexibility.

Non-local strictly proper scoring rules are investigated in Chapter 3 for the CRPS and more generally in Chapter 4. Explanations for the CRPS bias and the definition of a fair CRPS with classes of members are proposed. Besides, we introduced improved formulations for threshold-weighted and quantile-weighted scoring rules, as well as better understanding of these scores for model mixtures. In the case study of Appendix 7.A, results obtained with the CRPS or other scoring rules are not significantly different. This is a robustness argument in favor of the CRPS. For the verification of extreme events, it might be of interest to compare strongly asymmetric quantile-weighted scoring rules and simple quantile scores.

We defined a generalized least-square alternative to the CRPS in Section 4.2, which includes the observational distribution. New perspectives are brought on the connections between Pearson's χ 2 , the Anderson-Darling test, and Cramer-von Mises test. This work appears as a starting point to further theoretical analysis. Case studies are needed to compare this new loss against common learning and verification procedures.

Figure

  Figure A -Description de notre approche des prévisions PV en loi.

Figure B -

 B Figure B -Description of our approach to PV power forecasts.

Contents 1 . 1

 11 Forecasting photovoltaic power production with meteorological forecasts . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Weather forecasting . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Solar radiation forecasting . . . . . . . . . . . . . . . . . . . . 1.1.4 PV power data and statistical modeling . . . . . . . . . . . . 1.1.5 Probabilistic forecasts of PV power with meteorological forecasts 1.2 Sequential aggregation . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Algorithm evaluation with regret bounds . . . . . . . . . . . 1.2.3 Online learning algorithms . . . . . . . . . . . . . . . . . . . 1.2.4 Examples of loss functions . . . . . . . . . . . . . . . . . . . . 1.3 Probabilistic forecasting with non-local strictly proper scoring rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Binary case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Ranked and continuous case . . . . . . . . . . . . . . . . . . . 1.3.3 Examples with the CRPS . . . . . . . . . . . . . . . . . . . .

ECMWFFigure 1

 1 Figure 1.3 -Maps of solar radiation for 2013-06-13 at 12:00 (UTC): HelioClim real-time estimation (left), AROME forecast (center), HRES forecast (right).

Figure 1 . 4 -

 14 Figure 1.4 -Production data of one plant, normalized by the installed capacity. The temporal averages for 30-min, 1-h and 2-h resolutions are shown.

Figure 1 . 6 -

 16 Figure1.6 -Production data (black line with white dots) and local solar radiation forecasts (5 × 5 nearby grid points from AROME). The data are normalized by the daily maximal value. For the clear day situation (left), the solar forecasts of the 5 × 5 nearby grid points are equal.

  Map of CRPS of normal distributions for the observation y = 0, depending on the mean and standard deviation of the distribution. Ratio of CRPS of normal distribution defined by its mean and standard deviation against the CRPS of N (1, 1).

Figure 1

 1 Figure 1.9 -CRPS for normal distribution defined by its mean and standard deviation, for the observation y = 0.

Figure 1

 1 Figure 1.10 -Representation of the PDFs g a = N (0.6022, 0), g b = N (1, 1) and g c = N (0, (2.577) 2 ) (top left), and their respective CDFs G a , G b and G c (top right), reaching the same CRPS of 0.602 at 0.001 precision (for y = 0). The squared difference between H y the Heaviside function centered on y and respectively G a , G b and G c are shown at the bottom.

Figure 2 .

 2 Figure 2.1 -Annual average of 6-hour forecasts at 12:00 UTC in W m -2 , for the KMA ensemble mean in year 2012. The white circles exhibit the locations of ground observation sites for HelioClim evaluation.

Figure 2 . 2 -

 22 Figure 2.2 -Matrices of correlation between members (R2 matrices). Each row and column of the correlation matrices is dedicated to one given member. The rows and columns appear in the same order. Each matrix entry is the correlation coefficient between the merged data (all timesteps, 100 random locations) of a pair of members. On the left, the correlations are shown between raw members, as they are retrieved from TIGGE. In the middle, the correlations are computed after sorting within each center ensemble. On the right, the correlations are obtained after a full sorting of the forecasts.

Figure 2

 2 Figure 2.3 -Monthly RMSEs of TIGGE ensemble means in W m -2 , after linear conversion (lin); (blue: CMA), (green: ECMWF), (red: UKMO), (cyan: KMA), (magenta: CPTEC) and (yellow: Météo-France). The dashed line is the score of the reference forecast. The dasheddotted line is the score of a typical aggregated forecast (see Section 2.4.2).

Figure 2 . 4 -Figure 2

 242 Figure 2.4 -Rank histograms of the center ensembles. In the ideal case of a flat histogram, the height of all the bars would be equal to 1.

Figure 2 .

 2 Figure 2.6 -RMSE (solid line with disks, left axis) and explained variance (dashed-dotted line, right axis) of oracles with orthogonal members plotted against number of orthogonal members.

Figure 2 .Figure 2

 22 Figure 2.7 -Maps of statistical scores for the typical aggregated forecast (unitless for rRMSE, in W m -2 for MAE and RMSE).

Figure 2 .Figure 2 . 9 -

 229 Figure 2.8 -Annual averaged estimation for 12:00 UTC in W m -2 . Although it is based on individual members with low resolution, the aggregated forecast shows fine structures that are comparable to those of HelioClim.

Figure 2 .

 2 Figure2.10 -RMSE in W m -2 against the maximal number of members chosen from each center ensemble; (blue: no sorting), (green: sorting full ensemble), (red: sorting by center), (magenta: sorting by center and reference forecast included in the ensemble).

Figure 2 .

 2 Figure 2.11 -Impact of the horizon up to 42 h on the RMSE of the aggregated forecast and ensemble mean (lin); (dash dot: aggregated with 60 members), (blue: CMA), (green: ECMWF), (red: UKMO), (cyan: KMA), and (yellow: Météo-France).

  E C M e t e o F r a n c e

Figure 2 .

 2 Figure2.13 -Annual RMSEs of all members in W m -2 according to SSR-SSRD conversion method; (blue: no conversion), (green: add), (red: mult), (cyan: glob) and (magenta: lin). The abscissa is the member label and the members are grouped by origin. The dashed line is the score of the reference forecast. The dashed-dotted line is the score of a typical aggregated forecast (see Section 2.4.2).

  .13) This last quantity can be seen as the dispersion associated to the i.i.d. members of class C. Note the bias correction of E(|X C -X C |) with the factor M C (M C -1). Now we show how the ensemble CRPS and the class CRPS are related. Summing among classes (which belong to the partition C of the set of the members), we have C∈C W C M C c=1 |x C c -y|/M C = M m=1 u m |x m -y| . (3.14)

Figure 3

 3 Figure 3.1 -Ensemble CRPS (dotted gray) and adjusted ensemble CRPS (solid black), for ensembles of various sizes (10, 20, 50) from light gray to dark gray. The dispersion parameter (x-axis) is d ens . The scores are averaged over nearly 200 years of data (73,000 timesteps). The (solid black) lines of the adjusted CRPS are approximately at the same location for all ensemble sizes.

Figure 3 . 2 -Figure 3 . 3 -Figure 3

 32333 Figure 3.2 -Temporal evolution of the weights u m,t , with learning rate η = 0.05. The weights of members with correct dispersion are in black, and the weights of members with the incorrect dispersion d ens = 1.5 are in gray.

Figure 3

 3 Figure 3.5 -Average cumulated weights of the members with (possibly) incorrect dispersion parameter. Both learning algorithms based on the CRPS are tested: EG (solid line) and ridge (dotted lined). The white circles indicate that the algorithm is run for class CRPS (equal weights within the class) and the black squares indicate that the weights are computed explicitly for each member.

Contents 4 . 1

 41 Extension to threshold-weighted and quantile-weighted scoring rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Effect of threshold-weighting . . . . . . . . . . . . . . . . . . 4.1.2 Effect of quantile-weighting . . . . . . . . . . . . . . . . . . . 4.2 Probabilistic forecasting with observational noise . . . . . . 4.2.1 Generalized least square with the CRPS . . . . . . . . . . . . 4.2.2 Discussion and further work . . . . . . . . . . . . . . . . . . . Appendix 4.A Supplementary material . . . . . . . . . . . . . . .

Figure 4

 4 Figure 4.1 -Illustration of the step function φ.

Figure 4 . 2 -

 42 Figure 4.2 -Scaled weighting function ω/ω(0.5) in the symmetrical beta family ω(α) = α a-1 (1α) a-1 .

Figure 4 . 3 -

 43 Figure 4.3 -Illustration of the error function H y β 1(G) + (1 -H y )β 0 (1 -G)for y = 0 and G the CDF of a Gaussian distribution with mean µ = 1 and variance σ 2 = 1 (left), and y = 0 and G the CDF of a Gaussian distribution with mean µ = 0 and variance σ 2 = (2.577) 2 (right). The weighting functions ω belong to the symmetrical beta family, ω(α) = α a-1 (1 -α) a-1 : a = 0 (violet), a = 1/2 (red), a = 1 (black), a = 2 (cyan). The functions are scaled to reach the same value at G -1 (0.5).

Figure 4

 4 Figure 4.4 -Optimal quantile levels F(x m ) of the members for M = 5 (left) and M = 20 (right)with weighting functions ω belonging to the symmetrical beta family, ω(α) = α a-1 (1 -α) a-1 : a = 0 (violet), a = 1/2 (red), a = 1 (black), a = 2 (cyan).

Figure 4 . 5 -

 45 Figure 4.5 -Illustration of the quantile-weighting function ω (black), for a convex function ω = 1 √

Figure 4 . 6 -

 46 Figure 4.6 -Centered loss gradients m -M k=1 u k k of the CRPS (top) and the CRIGN (bottom) for M=101 linearly spaced members in [-2, 2] with observations y = 0 (blue), y = 1.2 (green) and y = 1.9 (red).

F

  j+1 -F j (4.22) This score can be seen as a discretized quadratic score * (taking only the numerators) with normalizing factors. * . The quadratic score of the PDF g against the observation y is -2g(y) + g 2 .

(j) m 2 in

 2 Equation 4.22. They accumulate to the bias of the score, which is equal to m

Figure 5

 5 Figure 5.1 -Illustration of weighted CDFs. The CDFs G and G are built with the same locations x m . However the weight u m or u m given to a member is different for G and G.

Figure 5 . 2 -

 52 Figure 5.2 -Time series of France production forecasts scaled by the installed capacity (12 hours of lead time, for several consecutive days). Top: equal weights for all members, (b): our forecast with online learning of the weights. Real production is in red and the median of the forecasted distribution is in white.

Figure 5

 5 Figure 5.4 -CRPS skill scores for all sites based on Equation 5.11. The confidence intervals correspond to the scores of all sites.

Figure 5 Figure 5 . 6 -

 556 Figure 5.5 -Spread skill diagrams for 36 h of lead times for all sites.

Figure 5

 5 Figure 5.8 -RMSE, MAE, CRPS and bias for the daily scores of France production. The results are shown for 3 forecasts: our weighted forecast, the raw forecast (all members with uniform weights), the deterministic forecast of the ECMWF (and its quantiles for the CRPS). The climatology scores are the following : bias = -0.001, CRPS = 0.055, MAE = 0.081, RMSE = 0.101.

Figure 5

 5 Figure 5.9 -CRPS Skill score.

Figure 5

 5 Figure5.12 -Reliability diagrams for lead times 33, 36, 39 hours; black round: the raw ensemble; magenta triangles: our weighted forecast.

Figure 6 . 2 -Figure 6 . 3 -Figure 6

 62636 Figure 6.2 -Comparisons of the predictive performance in MAE of the 100 × 100, 50 × 50 and 25 × 25 km 2 models.

Figure 6 . 7 -

 67 Figure 6.7 -Bias for each half-hour of our weighted forecast (Pond), the forecast with uniform weights (Raw), the climatological forecast (Climato) and AROME reference forecast (AROME Ref) for one power plant and for France production. All 336 forecasts scores are shown in gray.

Figure 6 . 8 -

 68 Figure 6.8 -MAE for each half-hour of our weighted forecast (Pond), the forecast with uniform weights (Raw), the climatological forecast (Climato) and AROME reference forecast (AROME Ref) for one power plant and for France production. All 336 forecasts scores are shown in gray.

Figure 6 Figure 6 Figure 6

 666 Figure6.12 -Hourly skill of the weighted forecast with 336 forecasts against the weighted forecast computed rolling quantile experts for the MAE, RMSE and CRPS. Negative skills indicate that better performance are reached thanks to the quantile experts.

Figure 6 .

 6 Figure 6.16 -Comparison of reliability diagrams between the weighted forecasts, computed for midday hours for each site.

Figure 6

 6 Figure 6.19 -Rank histograms of Raw.all (a) and rank histogram of Pond.all (b).

Figure 7 . 2 -

 72 Figure7.2 -Time series of our weighted probabilistic forecast for 3 consecutive days for Corsica PV production (observation in red). The differences in the time series of intraday lead times of 3 hours, 1 hour and 30 minutes show the evolution of the forecast sharpness and precision. The day-to-day forecast (above) does not benefit from intraday updates.

Figure 7

 7 Figure 7.3 -Evolution of the average spread of probabilistic forecasts with the lead time for Réunion (left) and Corsica (right). Morning and evening data are not included in these spread averages.

Figure 7 . 5 -

 75 Figure 7.5 -Weight boxplot of each member according to the lead time for Corsica. The boxplot bounds are the first and third quartiles. Morning and evening data are not included in these weight averages.

Figure 7

 7 Figure 7.6 -MAE (top) and CRPS (bottom) half-hourly scores for 30-minute lead time at Réunion. The MAE of all individual forecasts are shown in gray.

Figure 7 .

 7 Figure 7.7 -MAE (top) and CRPS (bottom) half-hourly scores for 30-minute lead time at Corsica. The MAE of all individual forecasts are shown in gray.

Figure 7 . 8 -

 78 Figure 7.8 -Evolution of the MAE with lead time for all half-hours at Réunion, net values (left), skill scores against the weighted forecast (right).

Figure 7 . 9 -

 79 Figure 7.9 -Evolution of the MAE with lead time for all half-hours at Corsica, net values (left), skill scores against the weighted forecast (right).

Figure 7

 7 Figure 7.10 -Evolution of the CRPS with lead time for all half-hours at Réunion, net values (left), skill scores against the weighted forecast (right).

Figure 7

 7 Figure 7.11 -Evolution of the CRPS with lead time for all half-hours at Corsica, net values (left), skill scores against the weighted forecast (right).

Figure 7

 7 Figure7.12 -CRPS skills of the weighted forecasts with restricted subsets for Réunion (left) and Corsica (right) against Pond.all (the weighted forecast with all members).

Figure 7

 7 Figure 7.13 -Rank histograms of the raw ensemble (a) and weighted ensemble (b) for midday hours: 30-min lead time Réunion (top left), 3-h lead time Réunion (top right), 30-min lead time Corsica (bottom left), 3-h lead time Corsica (bottom right).

Figure 7 .

 7 Figure 7.15 -Reliability diagrams of the raw ensemble (red) and weighted ensemble (black) for midday hours: 30-min lead time Réunion (top left), 3-h lead time Réunion (top right), 30-min lead time Corsica (bottom left), 3-h lead time Corsica (bottom right).

Figure 7

 7 Figure7.16 -Skill scores against our weighted forecast Pond for the cubic-CRPS, the CRPS, the arcsin-CRPS and the CRIGN. The subensembles Raw (top), "qr_prevenir" with 15 quantiles (bottom left) and "qr_prevenir" with 61 quantiles (bottom right) are evaluated.

  Table1.1 -Growth and installed capacity (in GW) for solar, wind, nuclear and hydro power.

		Growth 2015 Growth 2016 Installed
				capacity
	Solar Power	+0.9	+0.6	6.8
	Wind Power	+1.0	+1.3	11.7
	Nuclear Power			63.1
	Hydropower			25.5

Table 1 .

 1 ure 1.1. These elements are not perfect and possibly corrupted by noise. The purpose 2 -Meteorological forecasts and their resolutions over Europe in 2016. The lead time is indicated for the forecast of the 00 UTC analysis. The temporal resolution of PEARP is only of 6-h after the lead time of 48 h, and the temporal resolution of ECMWF forecasts is of 6-h after the lead time of 150 h. The conversion from degrees to kilometers is roughly done with a factor 100.

		AROME Arpège PEARP HRES	ENS
	Spatial resolution	1.3 km	7.5 km 15.5 km	0.1 •	0.2 •
	Temporal resolution	1 h	1 h	3 h	3 h	3 h
	Lead time	2 days	4 days	4 days 10 days 15 days

  2 . The weight u is searched in a 2norm ball B M in R M . The assumption of bounded losses may not always hold, but this difficulty is overcome with an additional regularization term (u x t ) 2 [Vov01; AW01].

	We refer the reader to the introduction of Gerchinovitz [Ger11] for a concise analysis.
	Still, key points concerning the ridge regression are the logarithm regret bound, and
	its quadratic formulation. From an optimization perspective, the ridge regression is a
	very handy and powerful tool. A modified version of the ridge regression is used in
	Chapter 2.
	Exponentiated Weighted Average. A most standard algorithm, and among the
	oldest, is the Exponentiated Weighted Average (EWA) [LW94]. This algorithm uses
	linear losses M m=1 u m,t m,t with positive weights summing to one (defining the simplex
	P M ), and bounded losses m,t ∈ [a, b]. In EWA, the regularization is the Kullback-
	Leibler divergence of u t from u. The update rule of EWA with learning rate η reads
	u t+1 = argmin u∈P M	M m=1	u m ln	u m u m,t	+ ηu m m,t ,
	giving the update rule

Ensemble forecast of solar radiation using TIGGE weather forecasts and HelioClim database

	2 Medium-range forecasts (one day to two weeks) of solar radiation are com-
	monly assessed with a single forecast at a given location. In this paper,
	we forecast maps of surface solar irradiance, using ensembles of forecasts
	from the THORPEX Interactive Grand Global Ensemble (TIGGE) with a
	6-h timestep. We compare our forecasts with observations derived from Me-
	teoSat Second Generation (MSG) and provided by the HelioClim-3 database Chapter 2 is based on the paper Thorey et al. [Tho+15]: as gridded observations over metropolitan France. First, we study the en-Thorey, J., Mallet, V., Chaussin, C., Descamps, L. and Blanc, P. (2015). Ensemble sembles from six meteorological centers. Second, we use sequential aggrega-forecast of solar radiation using TIGGE weather forecasts and HelioClim database. tion to linearly combine all the forecasts with weights that vary in space and Solar Energy, 120, 232-243. time. Sequential aggregation updates the weights before any forecast, using
	available observations. We use the global numerical weather prediction from Chapter 3 is based on the paper Thorey et al. [TMB16]: the European Center for Medium-range Weather Forecasts (ECMWF) as a Thorey, J., Mallet, V. and Baudin, P. (2017), Online learning with the Continu-reference forecast. The issue of spatial resolution is discussed because the ous Ranked Probability Score for ensemble forecasting. Q.J.R. Meteorol. Soc., 143: low resolution forecasts from TIGGE are compared to high resolution irra-521-529. diance estimated from MSG data. We found that the TIGGE ensembles
	are under-dispersed but rather different from one to another. Aggregation Chapter 5 is submitted as a research article to the International Journal of Forecast-decreases the forecast error by 20%, and produces a more realistic spatial ing. pattern of predicted irradiance.
	This chapter was published as the paper: Thorey, J., Mallet, V., Chaussin, C.,
	Descamps, L., and Blanc, P. « Ensemble forecast of solar radiation using TIGGE
	weather forecasts and HelioClim database ». In: Solar Energy 120 (Oct. 2015), pp. 232-
	243.

Contents 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.2 Analysis of TIGGE solar radiation and HelioClim database 42

	2.2.1 Description of TIGGE data . . . . . . . . . . . . . . . . . . . 42
	2.2.2 Analysis of the TIGGE ensembles of forecasts . . . . . . . . . 43
	2.2.3 Reference performance measures . . . . . . . . . . . . . . . . 45
	2.2.4 Comparison with HelioClim . . . . . . . . . . . . . . . . . . . 45
	2

.3 Ensemble forecast strategy: sequential aggregation . . . . . 49

	2.3.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
	2.3.2 Sequential aggregation: method . . . . . . . . . . . . . . . . . 49
	2.3.3 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Table 2 .

 2 Note that KMA data are different from the other 1 -Overview of TIGGE ensembles available for solar radiation, with forecasts starting from 18:00 or 24:00 UTC for the following 108 to 360 hours (horizon). In total, there are 158 ensemble members.

	Center	Origin Number of members Run
	CMA	China	14	24:00
	ECMWF	UE	50	24:00
	UKMO	UK	23	24:00
	KMA	Korea	23	24:00
	CPTEC	Brazil	14	24:00
	Météo-France France	34	18:00

Table 2 .

 2 2 -Evaluation of HelioClim-3 estimation of SSRD compared to in-situ measurements (daily average of daytime data). The relative scores are given, followed by the absolute scores in brackets (in W m -2 ).

Table 2 .

 2 of 0.968 shows that 3 -Empirical conversion formula of TIGGE data x tigge from SSR to SSRD. The deterministic forecasts from ECMWF are named x ecmwf .

	Label Conversion formula
	glob	x tigge × 1.18
	mult x tigge ×	x ssrd ecmwf x ssr ecmwf
	add	x tigge -x ssr ecmwf -x ssrd ecmwf
	lin	x tigge × a center + b center
	Center	Slope a center Intercept b center R 2
			(W m -2 )
	CMA	1.11	18	0.72
	ECMWF	1.18	-21	0.78
	UKMO	1.10	19	0.81
	KMA	0.91	18	0.81
	CPTEC	0.82	41	0.48
	Météo-France	1.07	37	0.76

Table 2 .

 2

4 -Linear regression SSR-SSRD, based on data from year 2011.

3.1 Mathematical background . . . . . . . . . . . . . . . . . . . .

  Mallet, V., and Baudin, P. « Online learning with the CRPS for ensemble forecasting ». In: Quarterly Journal of the Royal Meteorological Society (2016).

	Ensemble forecasting resorts to multiple individual forecasts to produce a
	discrete probability distribution which accurately represents the uncertain-
	ties. Before every forecast, a weighted empirical distribution function is
	derived from the ensemble, so as to minimize the Continuous Ranked Prob-
	Contents
	3.3.1 Simple model . . . . . . . . . . . . . . . . . . . . . . . . . . .
	3.3.2 Experiments without online learning . . . . . . . . . . . . . .
	3.3.3 Experiments with weight updates . . . . . . . . . . . . . . . .
	Appendix 3.

ability Score (CRPS). We apply online learning techniques, which have previously been used for deterministic forecasting, and we adapt them for the minimization of the CRPS. The proposed method theoretically guarantees that the aggregated forecast competes, in terms of CRPS, against the best weighted empirical distribution function with weights constant in time. This is illustrated on synthetic data. Besides, our study improves the knowledge of the CRPS expectation for model mixtures. We generalize results on the bias of the CRPS computed with ensemble forecasts, and propose a new scheme to achieve fair CRPS minimization, with essentially no assumption on the distributions.

A slightly different version of this chapter was published as the paper: Thorey, J.,

3

.1.1 Bibliographical remarks . . . . . . . . . . . . . . . . . . . . . 3.1.2 The Continuous Ranked Probability Score (CRPS) . . . . . . 3.1.3 The ensemble CRPS . . . . . . . . . . . . . . . . . . . . . . . 3.1.4 Bias of the ensemble CRPS with underlying mixture model . 3.1.5 Mixture model described by classes of members . . . . . . . . 3.2 Online learning methods . . . . . . . . . . . . . . . . . . . . . 3.2.1 Theoretical background . . . . . . . . . . . . . . . . . . . . . 3.2.2 Ridge regression . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Exponentiated gradient . . . . . . . . . . . . . . . . . . . . . 3.3 Numerical example . . . . . . . . . . . . . . . . . . . . . . . .

A Identities implying CDFs . . . . . . . . . . . . . . Appendix 3.B Computation of the ensemble CRPS . . . . . . . Appendix 3.C Regret bound of the ridge regression with the CRPS

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 3 .

 3 2 -Parameters of the numerical experiment.

	s 1	s 2	A	B	ω 1	ω 2	T
	0.3 0.3 1.68 0.336 1/365.25 1/11 730

Table 5 .

 5 or gamma distributions sum for wind and precipitation applications[START_REF] Sloughter | Probabilistic wind speed forecasting using ensembles and Bayesian model averaging[END_REF][START_REF] Sloughter | Probabilistic quantitative precipitation forecasting using Bayesian model averaging[END_REF]. Non-homogeneous regression fits 2 -Forecast availability with lead time. PEARP is the Météo France ensemble, Det defines the deterministic forecasts Arpège and HRES, and ENS is the ECMWF ensemble.

	Label	Nature	Origin	Timestep Resolution Base time # of forecasts
	HRES	deterministic	ECMWF		3 h	0.13 •	0 h	1
	ARPEGE deterministic Météo France	1 h	0.10 •	0 h	1
	ENS	ensemble	ECMWF		3 h	0.25 •	0 h	50
	PEARP	ensemble	Météo France	3 h	0.20 •	18 h	34
	Table 5.1 -Forecast weather data. The indicated resolutions may change for further lead times
	than those of the present article.
		D	D + 1	D + 2 D + 3 D + 4 D + 5
		PEARP PEARP	x	x	x	x
		Det	Det	Det	Det	x	x
		ENS	ENS	ENS	ENS	ENS	ENS

.3 Discussion and perspectives . . . . . . . . . . . . . . . . . . . 146

  resolution forecasting system AROME delivers dense spatio-temporal information for short lead times. In this chapter, we generate numerous PV forecasts using AROME solar radiation forecasts by leveraging multiple statistical models and exploiting the available spatio-temporal information. These forecasts are then combined with online learning techniques. We study the calibration of the resulting forecasts and improve the calibration thanks to additional quantile predictions. Finally, AROME forecasts are combined with other forecasts from Météo France and ECMWF. The work was jointly carried out with Clément Dolou, whose internship at EDF R&D was supervised by the author of this PhD thesis and Christophe Chaussin. Improvements with rolling quantile experts . . . . . . . . . . 139 6.2.2 Comparison of AROME with other forecasts from Météo France and ECMWF . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 6

		Helioclim	AROME	ECMWF
		50
		48
	latitude	46
		44
		42
	Contents 6.1 Building an ensembles of forecasts from AROME forecasts 131 6.1.1 Leveraging the high spatio-temporal resolution . . . . . . . . 131 6.1.2 First sequential aggregation results with AROME meteorolog-ical experts . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.1.3 Adding rolling quantiles experts . . . . . . . . . . . . . . . . 138 6.2 Sequential aggregation results with AROME statistically calibrated experts . . . . . . . . . . . . . . . . . . . . . . . . . 139 -2 0 2 4 6 8 -4 -2 0 2 4 6 8 longitude -4 -2 0 2 4 6 Figure 6.1 -Maps of solar radiation (in W m -2 ) for 2013-06-13 at 12:00 (UTC): HelioClim 6.2.1 -4 real-time estimation (left), AROME forecast (center), ECMWF HRES forecast (right).	8

Table 6 .

 6

2 -Skill score against Pond.all of the three main deterministic production forecasts.

Table 6 .

 6 3 -MAE according to seasonality and forecasted level of production of Pond.all. The asterisk * and the star respectively indicate the worst and best score by line. These signs are attributed in case of a clear score difference (above 0.003) .

		MAE for sites 1-2
	Pond.1 Pond.smooth Pond.all Number of points
	Feb-May level>0.6 0.055	0.058	0.055	68
	Feb-May level<0.6 0.087	0.108 *	0.090	132
	May-Oct level>0.6 0.068	0.067	0.066	203
	May-Oct level<0.6 0.117	0.150 *	0.125	59
		MAE for sites 3-10
	Pond.1 Pond.smooth Pond.all Number of points
	Feb-May level>0.6 0.086	0.091 *	0.085	177
	Feb-May level<0.6 0.146 * 0.139	0.135	583
	May-Oct level>0.6 0.079 * 0.062	0.064	556
	May-Oct level<0.6 0.135	0.143 *	0.132	484

Table 6 .

 6 4 -Skill scores of the probabilistic forecasts against Pond.all.

Table 7 .

 7 1 -Summary of probabilistic forecasts in the set of experts.

  Figure 7.4 -Weight boxplot of each member according to the lead time for Réunion. The boxplot bounds are the first and third quartiles. Morning and evening data are not included in these weight averages.

		0.10		0.10	3-h lead time
		0.08		0.08
	Weight	0.04 0.06	Weight	0.04 0.06	prevenirQ qr_prevsat qrforest_prevsat qr_prevenir prevsat prevenir qr_prevenir prevsat prevenirQ qr_prevsat qrforest_prevsat prevenir
					qrforest_prevenir qrforest_prevenir
		0.02		0.02	Persistence Persistence
		0.00		0.00
			1 7 14 22 30 38 46 54 62 70 78 1 7 14 22 30 38 46 54 62 70 78
					Member index Member index
		0.10		0.10	1-h lead time 1-h lead time
		0.08		0.08
	Weight	0.04 0.06	Weight	0.04 0.06	prevenirQ qr_prevsat qrforest_prevsat qr_prevenir prevsat prevenir qr_prevenir prevsat prevenirQ qr_prevsat qrforest_prevsat prevenir
					qrforest_prevenir qrforest_prevenir
		0.02		0.02	Persistence Persistence
		0.00		0.00
			1 7 14 22 30 38 46 54 62 70 78 1 7 14 22 30 38 46 54 62 70 78
					Member index Member index
		0.10		0.10	30-min lead time 30-min lead time
		0.08		0.08
	Weight	0.04 0.06	Weight	0.04 0.06	prevenirQ qr_prevsat qrforest_prevsat qr_prevenir prevsat prevenir qr_prevenir prevsat prevenirQ qr_prevsat qrforest_prevsat prevenir
					qrforest_prevenir qrforest_prevenir
		0.02		0.02	Persistence Persistence
		0.00		0.00
			1 7 14 22 30 38 46 54 62 70 78 1 7 14 22 30 38 46 54 62 70 78
					Member index Member index

  Table7.3 -Daily averages of quantile-weighted scores according to the optimization criterion and the lead time.

	Evaluation		CRPS
	Optimization CRPS arcsin-CRPS CRIGN
	0.50 h	0.0213	0.0213	0.0214
	1.00 h	0.0287	0.0288	0.0290
	1.50 h	0.0332	0.0333	0.0337
	2.00 h	0.0365	0.0368	0.0372
	2.50 h	0.0390	0.0393	0.0397
	3.00 h	0.0409	0.0410	0.0415
	Evaluation		arcsin-CRPS
	Optimization CRPS arcsin-CRPS CRIGN
	0.50 h	0.0267	0.0266	0.0266
	1.00 h	0.0355	0.0355	0.0357
	1.50 h	0.0409	0.0410	0.0414
	2.00 h	0.0450	0.0452	0.0457
	2.50 h	0.0481	0.0483	0.0487
	3.00 h	0.0504	0.0504	0.0510
	Evaluation		CRIGN
	Optimization CRPS arcsin-CRPS CRIGN
	0.50 h	0.0733	0.0727	0.0726
	1.00 h	0.0965	0.0962	0.0967
	1.50 h	0.1115	0.1112	0.1122
	2.00 h	0.1226	0.1226	0.1236
	2.50 h	0.1313	0.1313	0.1322
	3.00 h	0.1381	0.1374	0.1387

Ces trois années ont passé à la fois très vite et très lentement, elles m'ont permis de rencontrer de nombreuses personnes comme Isabelle Herlin et tous les collègues de l'INRIA : Nathalie, Yann, Paul, Sylvain, Nicolas, Raphaël, Julien et Guillaume. Grâce à toi Isabelle, j'ai pu découvrir un peu les rouages d'une institution aussi riche que l'INRIA. Merci à tous pour votre enthousiasme et votre brin de folie, vos avis tranchés lors de débats passionnés autour du billard ou d'un tour de magie. Côté EDF, je souhaite remercier les groupes Météo appliquée et Environnement atmosphérique et plus particulièrement Stéphanie et Christophe pour leur gentillesse et leur attention, ainsi qu'Alain, l'homme au short qui ne sèche jamais, mais aussi Bénédicte, Sylvie, Augustin et tous les autres.Je tiens bien sûr à remercier ma famille pour son soutien indéfectible, les amis de toujours et les nouvelles personnes qui sont entrées dans ma vie. Et plus particulièrement ma femme Camille qui est toujours à mes côtés.

Remerciements

Je tiens à remercier tout d'abord mon superviseur de thèse Vivien Mallet, qui m'a fait confiance au cours de ces trois années en me proposant un sujet de thèse riche et une grande liberté pour avancer. Merci à toi Vivien pour ta grande disponibilité, ta patience et ton sens aigü du travail bien fait qui m'ont permis de progresser énormément.

Mes remerciements vont également à Jochen Broecker et Liliane Bel pour avoir accepté de rapporter ma thèse, ainsi qu'à Petra Friederichs, Olivier Mestre et Olivier Wintenberger en tant que membres du jury. Tout au long de ce travail, le comité de suivi de thèse composé notamment de Laurent Dubus, Luc Musson-Genon, Laurent Descamps, Yannig Goude et Philippe Blanc, m'a accompagné avec leurs précieux conseils. Merci également à Gilles Stoltz et Pierre Gaillard pour avoir répondu à mes questions théoriques.

Appendix 3.B Computation of the ensemble CRPS

We have

where Γ is the upper bound of the integral. Because the weights sum to one, we get the last simplification.

We rewrite the above expression using Equation 3.30:

thus CRPS(G, y) = 2 QS α (G -1 (α), y)dα .

In the following, we elaborate on the impact of adding weighting functions φ(θ) or ω(α) in Equation 4.1 instead of uniform weighting schemes.

Effect of threshold-weighting

Here, we show how our analysis on the bias of the ensemble CRPS and on the class CRPS is easily extended to the threshold-weighted CRPS. The results of Section 3.1 hold up to a change of variables, suggested by Ehm et al. [START_REF] Ehm | Of quantiles and expectiles: consistent scoring functions, Choquet representations and forecast rankings[END_REF].

The CRPS, seen as a sum of Brier scores can be extended to the threshold-weighted CRPS [START_REF] Matheson | Scoring rules for continuous probability distributions[END_REF], where different weights are used for each threshold. Using the strictly positive function φ, the threshold-weighted CRPS is defined by: wCRPS(G, y) = (G -H y ) 2 φ .

(4.3)

Let Φ be the antiderivative function of φ, the function Φ is strictly increasing and invertible because the function φ is strictly positive. We now get into further details by generalizing Equation 3.29 and Equation 3.31 of Appendix 3.A. Let X and Z be independent random variables respectively described by the CDFs G and K, and let H be the unit step function centered on 0. We have (H -K)φ = E(Φ(Z)) -Φ(0) , (4.4) using integration by parts, and, as proved below,

We see that weighting thresholds with φ is equivalent to a data transformation with Φ. We now demonstrate Equation 4.5:

We have successively used the fact that GK is the CDF of max(X, Z), Equation 4.4, the fact that for any x, z ∈ R, Φ(max(x, z)) = max(Φ(x), Φ(z)) because Φ is an increasing function, and 2 max(x, z) = |x -z| + x + z.

With the same notation as in Section 3.1.4, we show that the bias of the ensemble threshold-weighted CRPS with underlying mixture model is therefore similar to the unweighted case of Equation 3.7 up to the data transformation with Φ. Let G E be the weighted CDF described by the members x m and the weights u m . The members are assumed to be samples from the independent random variables X m . The CDF

Appendix 4.A Supplementary material

Proof of R inverse formula in the discretized case:

Here we demonstrate that the matrix noted R -1 is indeed the inverse of R. We make an intensive use of R i,j = F min(i,j) (1 -F max(i,j) ) and of the matrix product formula (AB) i,j = k S A i,k B k,j .

Appendix 5.A Results for France production

In this Appendix, we show the results for France production, while the results for the individual sites are shown in Section 5.4. The results of France production forecasts and plants forecasts are roughly similar. Our online learning algorithm provides improvements over the raw ensemble up to a lead time of a few days. Because it is easier to forecast the power output of the total production, the forecast quality is better than for individual sites. This statement is verified for all diagnostic tools shown below.

We show in Figure 5.8 the average bias, CRPS, MAE and RMSE for France production. We see the scores of the sites are more than twice as large as the score of France production, but for the bias. Our online learning algorithm provides improvement for bias, CRPS, MAE and RMSE up to a lead time of 4 days.

The CRPS skill scores are shown in Figure 5.9. Once again, the score trends are mostly equivalent to those obtained for the sites. The quantile ensemble from HRES (quantile det ECMWF) has good scores for short lead times and the raw ensemble is the best ensemble with uniform weights after 24 h of lead time. Our online learning algorithm provides an improvement of roughly 10% over the raw ensemble for the first 24 h of lead time. This improvement decreases with time quickly than for the sites. It is remarkable that the CRPS skill score of the quantile ensemble from Arpège ("quantile det Arpège") shows much better results for the plants than for France production. Indeed the skill score of "quantile det Arpège" is around -15% for the plants and is stable, while it is at least below -24% for France production. For days 5 and 6, the weights brought by our algorithm do not vary much from the uniform distribution. Consequently, the skill scores are close to one.

The following probabilistic diagnostic tools are only for 0900, 1200, and 1500 of day 2 (lead times 33, 36, 39 hours). In Figure 5.10, we compare the rank histogram for the raw ensemble and our weighted forecast. The raw ensemble is largely under-dispersive with a positive bias (over-estimation). Our online learning algorithm manages to reduce the under-dispersion of the raw ensemble. This statement is verified on the spread skill diagram in Figure 5.11. We see that the spread and errors of our weighted forecasts match approximately, while the other ensembles (with uniform weights) are underdispersive with respect to their errors.

A correction of the forecast reliability is illustrated in Figure 5.12. The event "the production level is lower than the average production" is used (same as for the sites). +30 min and -30 min forecasts for several days. We see that the time-shifted forecasts may make up for inaccurate forecasts of ramp events. This is especially true for the second day in Figure 6.4. Yet, the mean production level and the amplitude of the production variations are not correctly described at this point. In a similar manner, space-shifted forecasts are built using spatial information. Besides the spatial average of the region of averaging, the 0.25, 0.50 and 0.75 spatial quantiles of the region are also used as inputs of the statistical model built with the average solar radiation. In other words for the 100 × 100 km 2 region, we sort the 40 × 40 = 1600 solar radiation forecasts of each 2.5 km. Then the median, the 0.25-quantile and the 0.75-quantile solar radiation values are picked to generate three additional production forecasts. Another option would be to generate forecasts for all the grid-points of the region, but this option is obviously too costly. We show for several consecutive days the production forecasts corresponding to spatial quantiles in Figure 6.5. The mean forecast is flanked by the spatial quantile forecasts. We see that this method may generate forecast peaks, which illustrate the sharp spatial variations of AROME solar radiation.

Model parameterization. Our ensemble of forecasts already comprises 3 (smoothing areas) × 3 (time-shifted) × 4 (spatial estimates) = 36 forecasts, and yet only 3 statistical models were built with different smoothing radii. We kept the 3 smoothing areas because the temporal and the spatial variations seen over the 3 areas differ with the area extension. Depending on the power plant, forecasts focusing on either small-scale or large-scale variations may be of higher interest. The following configurations are also generated: -Statistical model without clear sky normalization.

-Seasonal bias reduction, described in Section 5.1.

-Clear-sky model with timeshift-optimistic calibration. After clear-sky normalization, we pick the normalized solar radiation in the 3 time-shifted estimates being the closest to the normalized production to train the model. The motivation behind is that such statistical model is presumably built with an artificially good training data set, less sensitive to weather forecasts inaccuracies. weighted forecast without quantile experts Pond.0. The poor performance of quantile random forest experts, especially for France production, may be due to an insufficient amount of learning data (90 points). These scores and the weight distribution clearly point out that the 20 quantile experts help the weighted forecast Pond.1. Indeed, the weights given to the new experts are more than twice higher than the uniform weight 1/M for both quantile regression and quantile random forest experts.

Comparison of Pond.0 with Pond.1. Half-hourly skills of Pond.0 against Pond.1 indicate that our new setting does not improve the MAE, but slightly improves the RMSE around 2%, and CRPS gains are much higher, up to 10%, see Figure 6.12. Hence a major difference between Pond.0 and Pond.1 is the improvement of the spread of the probabilistic forecasts, since the mean of the probabilistic forecast is only slightly improved. A graphical verification of time series supports this statement. We clearly see a spread correction between Pond.1 and Pond.0 in Figure 6.13. Now we check the calibration of Pond.0 and Pond.1. The new weighted forecast shows a rank histogram closer to a flat rank histogram, see Figure 6.14, mainly by reducing the height of the outer bars. This is consistent with the increase of the forecast spread. Also, spread skill diagrams of Pond.1 show a better agreement with the first diagonal than those of Pond.0, thanks to the wider spread of Pond.1 in Figure 6.15. Reliability diagrams for the event "the production is below the climatological median" of Pond.0 and Pond.1 appear to be quite similar in Figure 6.16. The forecast Pond.0 provides better reliability for predicting events with probability below 0.4, and both forecasts tend to over-predict the occurrence of low production with probability between 0.5 and 0.9.

We also investigated the effect of removing the unnecessary experts. The online learning algorithm was run with all 336 + 20 experts. Interestingly, this weighted combination gives very close results than the weighted combination where unnecessary experts are removed (not shown). Consequently, the online learning algorithm handles quite nicely a large amount of poor experts. .1 -Time series of our weighted probabilistic forecast for 3 consecutive days for Réunion PV production (observation in red). The differences in the time series of intraday lead times of 3 hours, 1 hour and 30 minutes show the evolution of the forecast sharpness and precision. The day-to-day forecast (above) does not benefit from intraday updates.

Appendix 7.A Empirical results of quantile-weighted scoring rules with real-world data

Quantile weighted scoring rules, introduced in Chapter 4, offer an alternative to the CRPS. They allow higher focus on either the distribution tails or the distribution center. Examples of quantile-weighted scores are shown in Table 7.2, with the notation of Chapter 4. In this appendix, we provide empirical results on quantile-weighted scoring rules with real-world data. We focus on the data set of Réunion, because similar results were found for Corsica and for the data set mixing ECMWF, AROME and other Météo France forecasts of Chapter 6.

First, we investigate whether similar results are obtained with the quantile-weighted scoring rules described in Table 7.2. In the remaining of this thesis, the CRPS is the main score for probabilistic forecasts evaluation. Hence we wish to know whether the same conclusions may be drawn with other evaluation tools. Skill scores against ω(α)

Score S(G, y)

Thesis conclusions

Our case studies demonstrate that combining forecasts derived from several meteorological centers or postprocessing techniques enables to improve the accuracy of solar radiation and PV power forecasts. To do so, we resorted to online learning techniques providing update rules of the combination weights. These techniques come with theoretical performance guarantees on the predictive power of the combination of the forecasts, under essentially no assumptions. These methods do not depend on the nature of the output variable, but their implementation on weather-related variables (solar radiation and PV power) are particularly interesting under at least two aspects: the large uncertainty of the forecasts due to the difficulty to forecast clouds at the correct time and location, and the spatial structure of the forecasts.

The following practical results are emphasized in our case studies:

-The weighted combination of forecasts more than often outperforms uniform averages of a subset of forecasts. -Including forecasts of poor quality in the ensemble does not degrade the accuracy of the weighted forecasts. -Online learning techniques minimizing the CRPS improve the calibration and the reliability of probabilistic forecasts, as verified for several evaluation tools (CRPS, rank histogram, spread skill diagram, and reliability diagram). We now briefly review more specific results obtained for each case study. Chapter 2 focuses on non-probabilistic forecasting of solar radiation. Several TIGGE ensembles of solar forecasts are compared against satellite-derived HelioClim maps. This work is among the first examples of sequential aggregation techniques applied to maps of forecasts with the work of Baudin [START_REF] Baudin | Prévision séquentielle par agrégation d'ensemble: application à des prévisions météorologiques assorties d'incertitudes[END_REF] and Zamo [START_REF] Zamo | Post-traitements statistiques de prévisions de vent déterministes et d'ensembles sur une grille[END_REF]. Our forecast favorably compares against HRES (ECMWF reference forecast) and spatial patterns are more finely described.

For PV power forecasting, statistical models between meteorological and production data are built, which brings an additional challenge. Both the meteorological and the conversion model inaccuracies should be taken into account when delivering probabilistic forecasts of PV power. Chapter 3 details a proof-of-concept case study with ECMWF and Météo France forecasts, where online learning techniques meet CRPS minimization and the calibration of probabilistic forecasts. This case study covers 219 PV sites, with the high 30-min temporal resolution. Even for the long lead time of 6 days, the superiority of weather forecasts over climatological averages is verified.

We include AROME forecasts in our study in Chapter 6. A large amount of members are generated to take into account the rich spatio-temporal information of AROME. Statistically calibrated forecasts from quantile regressions greatly improve the CRPS of the weighted forecasts. Besides, for short lead times below 24 h, the weighted combination of PV power forecasts generated with only AROME and HRES weather forecasts