Solar Forecasting in a Challenging Insular Context

This paper aims at assessing the accuracy of different solar forecasting methods in the case of an insular context. Two sites of La Reunion Island, Le Tampon and Saint-Pierre, are chosen to do the benchmarking exercise. Reunion Island is a tropical island with a complex orography where cloud processes are mainly governed by local dynamics. As a consequence, Reunion Island exhibits numerous micro-climates. The two aforementioned sites are quite representative of the challenging character of solar forecasting in the case of a tropical island with complex orography. Hence, although distant from only 10 km, these two sites exhibit very different sky conditions. This work focuses on day-ahead and intra-day solar forecasting. Day-ahead solar forecasts are provided by the European Center for Medium-Range Weather Forecast (ECMWF). This organization maintains and runs the Numerical Weather Prediction (NWP) model named Integrated Forecast System (IFS). In this work, post-processing techniques are applied to refine the output of the IFS model for day-ahead forecasting. Statistical models like a recursive linear model or a nonlinear model such as an artificial neural network are used to produce the intra-day solar forecasts. It is shown that a combination of the IFS model and the neural network model further improves the accuracy of the forecasts.


Introduction
Accurate solar forecasts at various time steps are needed in order to increase the integration of renewables into electricity grids [1]. This statement is reinforced in the case of insular grids. The fluctuating character of the solar energy together with the fact that the island electricity grid is isolated may put in danger the stability of the grid and consequently the balance between supply and demand.
In addition, solar forecasting may be very challenging in an insular context as islands (like Réunion Island for instance) may usually experience a high spatial and temporal variability of the solar resource [2]. Previous works regarding irradiance forecasting for the power PV prediction of grid-connected systems were mainly done for large-scale continental grids [3][4][5]. Due to the small scale of the climatic phenomena, forecasting the solar irradiance in insular territories addresses new issues.
In order to cope with specific plant operations, solar forecasts must be provided with different granularities and horizons [6]. In addition, the recent development of grid-connected storages associated with intermittent renewables (solar, wind and wave) also requires power PV forecasts in order to optimize their operational management [7,8]. PV power output is directly related to the level of received global solar irradiance. Thus, forecasting the global solar irradiance is the key feature for power PV forecasting. In this work, we will focus on day-ahead and intra-day solar irradiance forecasts with an hourly granularity. cloud-motion vectors (CMVs) from satellite images are the most suitable [10,11]. In addition, finally, for very short-term horizons, from few minutes up to six hours, the literature is dominated by statistical approaches based on linear or nonlinear models [12][13][14][15][16].
NWP models may exhibit (large) biases [9]. Post-processing methods like Model Output Statistics (MOS) are frequently applied to reduce these biases or to refine the output of NWP models, as detailed local weather features are generally not resolved by NWP predictions [3]. In addition, spatial averaging can contribute to reduce forecast errors, particularly in situations with variable clouds, that are difficult to predict ( [4,17]). To our best knowledge, these post-processing techniques have been only tested in a continental context ( [3][4][5]) and it appears interesting to assess in this survey the performances of these techniques in the case of an insular environment. In this work, day-ahead forecasts are provided by the European Center for Medium-Range Weather Forecast (ECMWF). This organization maintains and runs the global Numerical Weather Prediction (NWP) model Integrated Forecast System (IFS). Two post-processing techniques will be used to refine the IFS forecasts for the insular sites of Saint-Pierre and Le Tampon in Réunion Island.
Regarding the intra-day solar forecasting, two statistical models namely an Artificial Neural Network (ANN) and a recursive linear technique (ARMA.RLS) will be evaluated. Furthermore, some day-ahead parameters produced by ECMWF will be used as exogenous inputs of the ANN model. It will be shown that the combination of the IFS model and the ANN model further improves the accuracy of the forecasts.
The remainder of this paper is organized as follows: Section 2 is devoted to the sites' description and data preprocessing. Section 3 sets out an in-depth analysis of the sky conditions experienced by each site while Section 4 will detail the errors metrics used to assess the forecasting accuracy of the different methods. Section 5 describes the day-ahead solar forecasts with a special emphasis on the MOS post-processing techniques. Section 6 presents the results of the intra-day solar forecasting techniques. Finally, Section 7 gives some concluding remarks.

Sites Description
Réunion Island exhibits a particular meteorological context dominated by a large diversity of micro-climates. Two main regimes of cloudiness are superposed: the clouds driven by synoptic conditions over the Indian Ocean and the orographic cloud layer generated by the local reliefs [2]. The data used to build the models are Global Horizontal Irradiance (GHI) measured by the PIMENT laboratory at the meteorological stations of Saint-Pierre and Le Tampon located in the southern part of Réunion Island (see Figure 1). Saint-Pierre is a coastal site while Le Tampon is an inland site. Table 1 gives the sites' characteristics.

Solar Radiation Data
Measurements are available on an hourly basis and two years of data (2012 and 2013) are used respectively for the building and appraisal of the models. The stations measure the GHI every six seconds and the 1-min averages are recorded. The hourly data correspond to the average of the previous 60 min of measurements. The solar irradiance is measured with a secondary standard pyranometer (CMP 11 from Kipp and Zonen). The uncertainty of the pyranometers is˘3.0% for the daily sum of GHI. Measurement quality is an essential asset in any solar resource forecasting study. The sites of Saint-Pierre and Le Tampon are well maintained and have followed the radiometric techniques regarding calibration, maintenance and quality control. Each data point has been processed with SERI-QC quality control procedure [18].
In this work, for comparison purposes, we also include two continental US sites namely Desert Rock (36.6N; 116.0W; 1007 m¨a.s.l) and Fort Peck (48.3N; 105.1W; 634 m¨a.s.l). These two stations are part of the NOAA's SURFRAD network [19].

Clear Sky Index
It corresponds to the ratio of the measured GHI to the theoretical GHI observed under clear sky conditions pGH I Clear q. With this methodology, the models designed in this work are dedicated to the forecasting of the stochastic part of the global radiation due to the cloud cover, leaving the geometric and the deterministic part to be modeled by the clear sky model.
The clear sky irradiance is generated with the BIRD model [20]. This simple model calculates estimates of GH I Clear with acceptable accuracy and with only few inputs [21]. For this study, the Aerosol Optical Depths (AODs) and the water vapor column parameters were retrieved from the AERONET site REUNION_ST_DENIS (20˝52'S; 55˝28') [22] located in the north of the island. 4.5 years were used to compute the climatological means of these atmospheric input parameters. The latter were used as inputs of the BIRD model. Similarly, the ozone atmospheric content was retrieved from the World Ozone Monitoring Mapping provided by the Canadian government [23]. The values of the input parameters of the BIRD model used for the two studied sites are given in Table 1.
The forecasting accuracy of the models proposed in this work depends on the accuracy of the clear sky model used to derive the clear sky index. In order to evaluate this error induced by the BIRD model, only the clear sky periods are considered. The clear sky hours were detected using the Ineichen method [24] applied to the two years of hourly-measured global irradiance. The last three lines of Table 1 give the number of clear sky hours detected, the relative Root Mean Square Error (rRMSE) and the relative Mean Bias Error (rMBE)-see the definition of these metrics below of the corresponding clear sky irradiance produced by the BIRD model. The performances of the BIRD model in this work are consistent with previous results [25]. Notice also that the site of Le Tampon experiences fewer occurrences of clear sky hours (see Table 1).

Filtering Methodology
Low solar elevations induce complex reflection phenomena that are not properly taken into account by the pyranometer. As a consequence, the values of the measured GHI are often not reliable. Furthermore, the amount of solar energy received at ground level in this condition is very small. As a consequence, data corresponding to a solar zenith angle (Θ Z ) larger than to 85˝are removed. This filtering removes less than 1% of the total annual sum of solar energy. Put in other words, night times and low solar elevations were not taken into account for the building and the test of the models. In addition, this filtering process allows to discard data with less precision as measurement uncertainties associated to pyranometers are typically much higher than˘3.0% for Θ Z > 85˝.

Sites Analysis
This section aims at analyzing the sky conditions experienced by each site. Figure 2 plots the distribution of the clear sky index for each site. In Table 2, based on the previous work done by [26], we define thresholds to distinguish different sky conditions. We also include in Figure 2, the kt˚distribution of the two continental sites namely Desert Rock and Fort Peck. The site of Desert Rock is representative of an arid climate dominated by clear skies. One may also notice that the kt˚distributions of the sites of Fort Peck and Saint-Pierre are quite similar. Indeed, as mentioned above, Saint-Pierre is a coastal site with rather high occurrences of clear skies. Conversely, the site of Le Tampon exhibits higher occurrences of cloudy/intermediate skies.
The characterization of the variability regime of solar irradiation at a particular site is a key factor to understand the level of error of a solar forecasting method. Put in other words, the variability experienced by a site has an impact on the forecasting accuracy of the different forecasting methods. Table 1 gives a metric that characterizes the annual solar variability of a site. This metric proposed by [27] is the standard deviation of the change in the clear sky index defined by: Var The values of solar variability in Table 1 have been calculated for a time scale Δt = 1 hour (as we deal with hourly GHI records) and for a time span that corresponds to two years of hourly data. A site with variability above 0.2 is considered as experiencing very unstable conditions [27,28].
As seen, the variability of the site of Le Tampon is above this threshold. Conversely, the variability of Desert Rock and Fort Peck are respectively 0.14 and 0.18. One may also notice that the variability of Saint-Pierre is close to that of Fort-Peck. Sub-Section 6.5.2 below will discuss the impact of the sky conditions on the forecasting accuracy of the different methods.

Error Metrics
In the realm of the solar forecasting community, the commonly used error metrics for point forecasts are: the Root Mean Square Error (RMSE), Mean Absolute Error (MAE) and Mean Bias Error (MBE). These error metrics are defined by the following equations: where N is the number of points in the dataset for the considered period. Relative values of these metrics (rRMSE, rMAE and rMBE) are obtained by normalization to the mean ground measured irradiance of the considered period. As noted by Lorenz in [1], the MAE is appropriate for applications with linear cost functions, that is, where the costs that are caused by a wrong forecast are proportional to the forecast error. The RMSE is more sensitive to large forecast errors, and hence is suitable for applications where small errors are more tolerable and larger errors cause disproportionately high costs. In the following, we provide the standard set of error metrics but with a special emphasis on the rRMSE.
In this work, we also include an additional metric: the forecast skill parameter s. The latter proposed by [29] is given by (%) = 1 − × 100 where stands for the RMSE The values of solar variability in Table 1 have been calculated for a time scale ∆t " 1 h (as we deal with hourly GHI records) and for a time span that corresponds to two years of hourly data. A site with variability above 0.2 is considered as experiencing very unstable conditions [27,28].
As seen, the variability of the site of Le Tampon is above this threshold. Conversely, the variability of Desert Rock and Fort Peck are respectively 0.14 and 0.18. One may also notice that the variability of Saint-Pierre is close to that of Fort-Peck. Sub-Section 6.5.2 below will discuss the impact of the sky conditions on the forecasting accuracy of the different methods.

Error Metrics
In the realm of the solar forecasting community, the commonly used error metrics for point forecasts are: the Root Mean Square Error (RMSE), Mean Absolute Error (MAE) and Mean Bias Error (MBE). These error metrics are defined by the following equations: RMSE " where N is the number of points in the dataset for the considered period. Relative values of these metrics (rRMSE, rMAE and rMBE) are obtained by normalization to the mean ground measured irradiance of the considered period. As noted by Lorenz in [1], the MAE is appropriate for applications with linear cost functions, that is, where the costs that are caused by a wrong forecast are proportional to the forecast error. The RMSE is more sensitive to large forecast errors, and hence is suitable for applications where small errors are more tolerable and larger errors cause disproportionately high costs. In the following, we provide the standard set of error metrics but with a special emphasis on the rRMSE.
In this work, we also include an additional metric: the forecast skill parameter s. The latter proposed by [29] is given by s p%q "ˆ1´R MSE method RMSE pers˙ˆ1 00 where RMSE method stands for the Further, the higher the s-skill score is, the better the improvement is.

Initial ECMWF Forecast
The European Center for Medium-Range Weather Forecast (ECMWF) is an intergovernmental organization that provides operational forecasts. It maintains and runs the numerical weather prediction (NWP) model Integrated Forecast System (IFS). NWP models outperform forecasts based on satellite data or ground measurements (only) for horizons longer than 5 h [1,11]. In order to supply forecasts that can be used by the grid operator for day-ahead scheduling, we retrieved the data generated by the IFS at 12 h 00 UTC (16 h 00 in Réunion Island). The forecasts correspond to hourly data with a spatial resolution of 0.125˝ˆ0.125˝(approximately 14 kmˆ14 km in Réunion-see Figure 3).
Atmosphere 2016, 7, 18 6 of each new proposed forecasting method and is the RMSE of the persistence model. With this definition, the persistence model has a forecast skill s = 0% and a value of s = 100% denotes a perfect forecast. Negative values of s indicate that the forecasting model fails to outperform the persistence model while positive values of s means that the forecasting method improves on persistence. Further, the higher the s-skill score is, the better the improvement is.

Initial ECMWF Forecast
The European Center for Medium-Range Weather Forecast (ECMWF) is an intergovernmental organization that provides operational forecasts. It maintains and runs the numerical weather prediction (NWP) model Integrated Forecast System (IFS). NWP models outperform forecasts based on satellite data or ground measurements (only) for horizons longer than 5 hours [1,11]. In order to supply forecasts that can be used by the grid operator for day-ahead scheduling, we retrieved the data generated by the IFS at 12 h 00 UTC (16 h 00 in Réunion Island). The forecasts correspond to hourly data with a spatial resolution of 0.125° × 0.125° (approximately 14 km x 14 km in Réunion-see Figure 3).

Post-Processing the IFS Model
In order to obtain an optimized local prediction, post-processing techniques are used to refine the output of NWP models. More precisely, post-processing the NWP model output consists in using a statistical method (also called Model Output Statics (MOS)) to correct systematic deviations due to either model errors or local influences not resolved by the NWP model. Indeed, although spatial resolution has increased during the last years, the NWP models still do not resolve the local weather details [1]. Therefore, there is room for improvement as regards the original NWP forecasts.
For the case study of Germany, Lorenz et al. [3] showed that the ECMWF forecasts could be refined with a MOS technique that consisted of a bias correction. In this work, a prior step to this bias correction was to apply a spatial averaging procedure on the original ECWMF forecasts.

Spatial Averaging
Spatial averaging of irradiance forecasts can lead to improved forecasts (notably in terms of RMSE) by smoothing the variations in variable sky conditions (due to changing cloud cover) [3,17]. For the site of Saint-Pierre and Le Tampon, we applied the spatial averaging technique (i.e., arithmetic average of surrounding pixels) to potentially improve the accuracy of the day-ahead forecasts. As shown by Figure 4, the spatial averaging technique did not improve the forecast accuracy. Unlike [3] which obtained the lowest RMSE by averaging over 4 × 4 pixels, here, the lowest RMSE is obtained

Post-Processing the IFS Model
In order to obtain an optimized local prediction, post-processing techniques are used to refine the output of NWP models. More precisely, post-processing the NWP model output consists in using a statistical method (also called Model Output Statics (MOS)) to correct systematic deviations due to either model errors or local influences not resolved by the NWP model. Indeed, although spatial resolution has increased during the last years, the NWP models still do not resolve the local weather details [1]. Therefore, there is room for improvement as regards the original NWP forecasts.
For the case study of Germany, Lorenz et al. [3] showed that the ECMWF forecasts could be refined with a MOS technique that consisted of a bias correction. In this work, a prior step to this bias correction was to apply a spatial averaging procedure on the original ECWMF forecasts.

Spatial Averaging
Spatial averaging of irradiance forecasts can lead to improved forecasts (notably in terms of RMSE) by smoothing the variations in variable sky conditions (due to changing cloud cover) [3,17]. For the site of Saint-Pierre and Le Tampon, we applied the spatial averaging technique (i.e., arithmetic average of surrounding pixels) to potentially improve the accuracy of the day-ahead forecasts. As shown by Figure 4, the spatial averaging technique did not improve the forecast accuracy. Unlike [3] Atmosphere 2016, 7, 18 7 of 17 which obtained the lowest RMSE by averaging over 4ˆ4 pixels, here, the lowest RMSE is obtained for the nearest pixel. This result may come from the small scale of the island or also because the average predictions for land and surrounding ocean are mixed. Therefore, we can conjecture that the spatial averaging of the global IFS forecasts may be irrelevant in the case of small insular territories.
Atmosphere 2016, 7, 18 7 for the nearest pixel. This result may come from the small scale of the island or also because the average predictions for land and surrounding ocean are mixed. Therefore, we can conjecture that the spatial averaging of the global IFS forecasts may be irrelevant in the case of small insular territories. In order to further illustrate the challenging character of solar forecasting in the case of a small tropical island with complex orography, Figure 5a plots the ground observations of Saint-Pierre vs. the ground observations of Le Tampon. As seen, although there is a small distance between the two sites, the discrepancy between the two ground measurements is relatively important. A correlation coefficient of 0.64 calculated between the two * time series of each site reinforces this statement. As shown by Figure 5b, it is not the case for the ECMWF forecasts. A correlation coefficient of 0.96 confirms that the ECMWF forecasts are rather similar for the two sites as a consequence of the coarse spatial resolution of ECMWF for a territory like Réunion Island.

Bias Correction with An Artificial Neural Network
Consistent error patterns allow for MOS to be used to produce a bias reduction function for future forecasts. In the case of Germany, as the original ECMWF forecasts showed a considerable In order to further illustrate the challenging character of solar forecasting in the case of a small tropical island with complex orography, Figure 5a plots the ground observations of Saint-Pierre vs. the ground observations of Le Tampon. As seen, although there is a small distance between the two sites, the discrepancy between the two ground measurements is relatively important. A correlation coefficient of 0.64 calculated between the two kt˚time series of each site reinforces this statement. As shown by Figure 5b, it is not the case for the ECMWF forecasts. A correlation coefficient of 0.96 confirms that the ECMWF forecasts are rather similar for the two sites as a consequence of the coarse spatial resolution of ECMWF for a territory like Réunion Island.
Atmosphere 2016, 7, 18 7 for the nearest pixel. This result may come from the small scale of the island or also because the average predictions for land and surrounding ocean are mixed. Therefore, we can conjecture that the spatial averaging of the global IFS forecasts may be irrelevant in the case of small insular territories. In order to further illustrate the challenging character of solar forecasting in the case of a small tropical island with complex orography, Figure 5a plots the ground observations of Saint-Pierre vs. the ground observations of Le Tampon. As seen, although there is a small distance between the two sites, the discrepancy between the two ground measurements is relatively important. A correlation coefficient of 0.64 calculated between the two * time series of each site reinforces this statement. As shown by Figure 5b, it is not the case for the ECMWF forecasts. A correlation coefficient of 0.96 confirms that the ECMWF forecasts are rather similar for the two sites as a consequence of the coarse spatial resolution of ECMWF for a territory like Réunion Island.

Bias Correction with An Artificial Neural Network
Consistent error patterns allow for MOS to be used to produce a bias reduction function for future forecasts. In the case of Germany, as the original ECMWF forecasts showed a considerable

Bias Correction with An Artificial Neural Network
Consistent error patterns allow for MOS to be used to produce a bias reduction function for future forecasts. In the case of Germany, as the original ECMWF forecasts showed a considerable kt (predicted ECMWF clear sky index) and the cosine of the solar zenith angle (cosΘ Z ). The corrected forecast is then obtained by subtracting the modeled bias from the original predicted values. This approach eliminated bias and reduced root mean square error (RMSE) of hourly forecasts by 5% for 24 h forecasts. Similarly to [3], we implement a bias correction, but based on an Artificial Neural Network (ANN) model instead of a polynomial model, which allows considering additional input variables ANNs are data driven approaches capable of performing a non-linear mapping between sets of input and output variables. An ANN with d inputs, m hidden neurons and a single linear output unit defines a non-linear parameterized mapping from an input vector x " px 1 , x 2 ,¨¨¨, x d q to an output y given by the following relationship: Each of the m hidden units are related to the tangent hyperbolic function f pxq " e x´e´x˘{`ex`e´x˘. The parameter vector w "` w j ( , w ji ( , b 1 , b 2˘g overns the non-linear mapping and is estimated during a phase called the training or learning phase. During this phase, the ANN is trained by the scaled conjugate gradient algorithm using a dataset (called training set) of N input and output examples. The second phase, called the generalization phase, consists of evaluating the ability of the ANN to generalize, that is to say, to give correct outputs when it is confronted with examples that were not seen during the training phase. Careful attention must be put on the building of the model, as a too complex ANN will easily overfit the training data. Several techniques like pruning [30], Bayesian regularization [31] or multi-objective genetic algorithm [32] can be employed to control the ANN complexity. In this work, we used the Bayesian Technique in order to automatically control the ANN complexity and therefore improve the generalization capability of the model [31]. The Bayesian approach offers significant advantages over the classical ANN learning process. Among others, one can cite the optimization of the ANN model using only the training dataset (i.e., no independent or validation dataset is needed). In addition, the Bayesian method offers a means to select the optimal number of hidden nodes by performing a model comparison. The interested reader is referred to [31] for details regarding the Bayesian approach.
In our application, we used a feed-forward ANN with only one hidden layer as it was theoretically proven that only one layer is sufficient to approximate any continuous function [33]. The ANN output y (i.e., the modeled bias) is related to three input variables namely the predicted clear sky index ( x kt ), the cosine of the zenith angle (cosΘ Z ) and the predicted total cloud cover ( z TCC). The MOS-corrected ECWMF forecasts are then obtained by subtracting the modeled bias from the original ECMWF forecasts.

Results of Day-Ahead Solar Forecasting and Discussion
An operational experimental set-up was used to implement the bias correction over the year 2012. More precisely, a sliding window of 90 days was used to train the ANN model and the bias correction was tested on the next day following the sliding window. In addition, let us recall that this MOS correction is not applied for Θ Z > 85˝. We compare our approach with the MOS technique employed by [3]. Let us recall that, in this previous work, the bias correction function was a polynomial of fourth order in the variables clear sky index ( x kt ) and the cosine of the zenith angle (cosΘ Z ). Table 3 lists the values of the error metrics.
As shown by Table 3, post-processing of the day-ahead forecasts produced by ECMWF allows improving their quality. Indeed, the bias of the corrected ECMWF forecasts is efficiently reduced. However, a relatively small improvement is observed for the quadratic error (rRMSE), which is in agreement with [3]. Overall, the ANN technique performs slightly better than the polynomial model. This improvement may come from the adding of the Total Cloud Cover (TCC) as a third input variable. Further, the MOS technique based on the 4th order polynomial model is a linear parameterized model while the ANN model is able to reproduce more complex non-linear relationships between the inputs and the output. In addition, we analyze the capability of the different models to distinguish between different sky conditions: clear sky, cloudy sky and overcast sky, defined by the thresholds given in Table 2. A forecast correctly predicting the sky conditions, i.e., falling in the same range of the measured clear sky index according to Table 2 is labeled as "good" forecast in the following. In other words, a "good" forecast corresponds to a correctly predicted event or state. Figure 6 shows the percentage of "good" forecasts for the initial ECMWF model and for the two post-processing methods. In this figure, the percent value is the number of correctly predicted states over the number of states as measured. As shown by Figure 6, according to the site under study, the behavior of the post-processing techniques is different. For the site of Saint-Pierre, the MOS techniques increase the rate of "good" forecasts in case of clear sky. However, the MOS techniques also reduced the percentage of "good" forecasts in case of cloudy or overcast skies. For Le Tampon, the bias correction method increases the rate of "good" forecasts in case of cloudy or overcast skies but decreases this rate in case of clear skies. Figure 6 gives also the number of correct forecasts for all-sky conditions. For Saint-Pierre, the overall share of correctly predicted events by the post-techniques is increasing. However, for le Tampon, no improvement is observed. Overall, regarding the two post-processing techniques, it appears again that ANN performs slightly better than the polynomial model.
Finally, while the MOS correction method led to a significant bias reduction, it appears that more work must be done in order to improve the results. As a bias correction can only correct systematic deviations, some original "good" forecasts have been unnecessarily corrected. Future work may consist in carrying an in-depth study of the error analysis that could pinpoint in detail the origin of the errors. This detailed bias analysis may lead (in addition to the clear sky index, TCC and cosΘ Z ) to the identification of another relevant input variables for the ANN.
Atmosphere 2016, 7, 18 9 parameterized model while the ANN model is able to reproduce more complex non-linear relationships between the inputs and the output. In addition, we analyze the capability of the different models to distinguish between different sky conditions: clear sky, cloudy sky and overcast sky, defined by the thresholds given in Table 2. A forecast correctly predicting the sky conditions, i.e., falling in the same range of the measured clear sky index according to Table 2 is labeled as "good" forecast in the following. In other words, a "good" forecast corresponds to a correctly predicted event or state. Figure 6 shows the percentage of "good" forecasts for the initial ECMWF model and for the two post-processing methods. In this figure, the percent value is the number of correctly predicted states over the number of states as measured. As shown by Figure 6, according to the site under study, the behavior of the post-processing techniques is different. For the site of Saint-Pierre, the MOS techniques increase the rate of "good" forecasts in case of clear sky. However, the MOS techniques also reduced the percentage of "good" forecasts in case of cloudy or overcast skies. For Le Tampon, the bias correction method increases the rate of "good" forecasts in case of cloudy or overcast skies but decreases this rate in case of clear skies. Figure 6 gives also the number of correct forecasts for allsky conditions. For Saint-Pierre, the overall share of correctly predicted events by the post-techniques is increasing. However, for le Tampon, no improvement is observed. Overall, regarding the two postprocessing techniques, it appears again that ANN performs slightly better than the polynomial model.
Finally, while the MOS correction method led to a significant bias reduction, it appears that more work must be done in order to improve the results. As a bias correction can only correct systematic deviations, some original "good" forecasts have been unnecessarily corrected. Future work may consist in carrying an in-depth study of the error analysis that could pinpoint in detail the origin of the errors. This detailed bias analysis may lead (in addition to the clear sky index, TCC and cos Θ ) to the identification of another relevant input variables for the ANN.

Reference Models
Following Perez et al. [11], we propose to test our forecasting methods against reference models like persistence, smart persistence and climatology.
The persistence (Pers) model is expressed as follows: x kt˚pt`hq " kt˚ptq Notice that in all the following equations describing the forecasting models, we used the following notations: a variable with a hat (ˆ) is a forecasted variable, h = 1,2,3, . . . 6 is the forecasting time horizon also called the lead time and t denotes the time when the forecasts are generated.
This model assumes that the clear sky index for each time horizon h depends only on the previous measured value kt˚ptq, which means that the sky conditions remain invariant between time 't' and time 't + h'. The next model represents an easy way to improve the persistence model and it is called Smart Persistence (Smart-Pers). It consists in forecasting the clear sky index for each time horizon 'h' as the mean of the 'h' previous measured clear sky values. Smart-persistence model is defined by the following equation: Finally, we propose the climatological mean model, which is independent of the forecasting time horizon. More precisely, this model performs a constant forecast of the clear sky index that corresponds to its mean historical value: x kt˚pt`hq " mean pkt˚q Unlike persistence, it presents a limit in terms of accuracy for longer horizons. In our case, we used the average clear sky index of the year 2012 in order to forecast the clear sky index of the year 2013.

Intra-Day Solar Forecasts with A Linear Recursive Model (ARMA.RLS Model)
The Auto Regressive Moving Average (ARMA) model is a popular linear technique in the realm of solar forecasting. In particular, it has been extensively studied in renewable energy forecasting and, owing to its parsimony, it has turned out to be a very tough competitor to beat. Applications include, among others, forecasting of wind power generation [34], online power forecasting [35] and wave energy flux [36].
A general formulation of an ARMA (p,q) model with p autoregressive (AR) terms and q moving average (MA) terms is given by [37]: ptq is an independent and identically distributed random variable with a zero mean. The vector Θ "`θ 0 .θ 1 .¨¨¨.θ p .Φ 1. .¨¨¨.Φ q˘c ontains the set of parameters to be estimated. A classic setting based on a least-squares method and a training data set can be used to estimate the set of parameters [37].
However, in this work, to estimate the model's parameters, we chose a variation of the least squares method, namely the Recursive Least Square (RLS) method (see [38] for details of implementation). This method offers the advantage of reducing the computational cost for estimating the model's parameters. In addition, the parameters are updated in real-time as new data become available. This contrasts with more intensive estimation methods operating on a sliding window where the estimation process is being carried out at each time step. Hence, the RLS method is particularly well suited in an operational context where forecasts have to be timely delivered. In addition, it must be stressed that contrary to non-linear models like ANNs, no training set is necessary here.
Regarding the structure of the ARMA model, the use of the classical ACF (Auto Correlation Function) and PACF (Partial Auto Correlation Function) techniques [37] led to the selection of the following orders p " 6 and q " 2.

Intra-Day Solar Forecasts with An ANN Model
In this work, two versions of the ANN approach are proposed. The first one consists in designing an ANN to predict next values of solar irradiance from only past measured values of the irradiance kt˚pt´iq i.e., no exogenous variables are used. This first model is given by the following equation: The number of past input values p is calculated with the auto mutual information factor (see [14] for details of implementation of the auto mutual information factor).

ECWMF Forecasts as A Means to Improve Intra-Day Solar Forecasting
In addition to the past ground measurements kt˚pt´iq , the second ANN takes additional exogenous inputs provided by the day-head ECMWF forecasts i.e., the forecasted clear sky index ( { ktE CMWF ), the total cloud cover ( { TCC ECMWF ) and the cosine of the solar zenith angle (cosθ z q for the considered forecast horizon h. This second model, denoted by the acronym ANN+ECMWF, has the following form: (12) Again, in both cases (Equations (11) and (12)), we used an ANN with only one hidden layer whose complexity was automatically controlled through the use of a Bayesian regularization technique [31] and year 2012 was chosen as the training dataset.
Finally, regarding the implementation of the ARMA.RLS and ANN models, we chose to create one model per forecasting time horizon in order to avoid the propagation of error that is usually observed while running the same model iteratively.

Models with Only Past Input Measurements
This section evaluates the accuracy of the different methods in the case of only past on-site measurements data as inputs for the models. Tables 4-6 give (on the one-year testing period i.e., year 2013) the rMBE, rRMSE and rMAE values of the different methods for forecasting time horizon ranging from 1h to 6h. Mean GHI is given for each site from which one can infer the absolute values from the relative values.
First, Table 4 shows also a slight improvement of the rMBE brought by the ANN and ARMA.RLS models compared to the Pers and Smart-Pers models.
Second, regarding the rRMSE and rMAE metrics, Tables 5 and 6 and Figure 7 show that the ANN and ARMA.RLS methods perform better than the reference models (Pers, Smart-pers and Climatology) for forecasting time horizons greater than one hour. Further, the gain in rRMSE or rMAE increases with the forecasting horizon.
Interestingly, Table 5 and Figure 7 show that the simple linear recursive ARMA.RLS model performs equally well than the nonlinear ANN model in terms of rRMSE. It even produces slightly better forecasts in the case of Saint-Pierre. However, in terms of rMAE, it appears that the ANN produces slightly better forecasts. One may also notice (see Figure 7) that the performances of the linear ARMA.RLS and nonlinear ANN models tend towards that of the climatological mean. This behavior is consistent, as these methods tend to asymptotically model the mean of the data.

13
One may also notice (see Figure 7) that the performances of the linear ARMA.RLS and nonlinear ANN models tend towards that of the climatological mean. This behavior is consistent, as these methods tend to asymptotically model the mean of the data.  Table 7 lists the s-skill scores of the different forecasting techniques for each forecasting time horizon. As shown by Table 7, the s-skill scores of the methods increase with the forecasting time horizon. The proposed models (ARMA.RLS, ANN) exhibit positive scores and therefore perform better than the persistence model.   Table 7 lists the s-skill scores of the different forecasting techniques for each forecasting time horizon. As shown by Table 7, the s-skill scores of the methods increase with the forecasting time horizon. The proposed models (ARMA.RLS, ANN) exhibit positive scores and therefore perform better than the persistence model. This section aims at assessing the impact of the sky conditions on the forecasting performance of the different methods. Figure 7 plots the performance of the different models (except the ANN + ECMWF model for the sites of Desert Rock and Fort Peck). Let us recall that the site of Desert Rock experiences high occurrences of clear skies and possesses a low annual variability. Fort Peck exhibits a somewhat less but still high occurrence of clear sky conditions with an annual variability similar to that of Saint-Pierre.
The performance of the models is obviously the best for the site of Desert Rock. At the opposite, more variable sky conditions lead to worse forecasting performances (case of Le Tampon).
As mentioned above, the s-skill scores of the ARMA.RLS and ANN methods demonstrate that these techniques outperform the references models (Persistence and Climatology) whatever the site under study. In other words, we can argue that the ARMA.RLS and ANN models offer a significant improvement over Persistence and Climatology by an amount that is fairly independent of the sky conditions experienced by a site.
It is interesting to note that the sites of Saint-Pierre and Fort Peck, which exhibit almost the same variability and kt˚distributions, display almost the same rRMSE values. At a first glance, one can conjecture a link between variability and forecastability. Some recent works [39,40] have shown that the type of climate and the associate cloud patterns have a profound impact in the forecasting performance. The authors of these works ([39,40]) also went further by proposing some metrics, which were related to solar variability to a priori assess the accuracy of solar forecasting models (i.e., to estimate the forecasting performance before any forecasts are produced). In addition, [17] showed increasing RMSE values with increasing variability.

Models with Exogenous Inputs Provided by ECMWF (ANN+ECMWF Model)
Tables 4-6 and Figure 7 clearly demonstrate the improvement (at each forecast horizon) of the forecasting accuracy obtained from the combination of day-ahead ECMWF forecasts with on-site measured irradiance. One may also notice that higher values of s-skill scores (Table 7) are obtained with the ANN+ECMWF model. Figure 7 also confirms that the rRMSE values of the ANN+ECMWF model tend towards that of the initial day-ahead ECMWF forecasts for Le Tampon and Saint-Pierre when the forecast horizon increases. The proposed combination of different type of forecasting models through the use of an ANN may be a viable alternative to the standard linear combination of forecasting models currently implemented by the solar forecasting community.

Main Conclusions
This work proposed to evaluate the accuracy of day-ahead and intra-day forecasts in a particular insular context. Two sites of Réunion Island, which are geographically close, but with very different sky conditions, were chosen for the benchmarking exercise.
Day-ahead forecasts were provided by the ECMWF. Model output statistics (MOS) methods were applied in an attempt to refine the ECMWF forecasts. It was shown that, due to the small scale of the island, a technique like spatial averaging was inefficient for the case of insular sites like Saint-Pierre and Le Tampon.
As the original day-ahead forecasts were biased, a bias correction method based on two techniques was tested. While the MOS correction method led to a significant bias reduction, it appears that more work must be done in order to improve the results. For instance, future work may consist in finding systematic dependencies of the bias in dependence on new parameters that could better describe the meteorological situation.
Intraday solar forecasts (from 1 h to 6 h ahead) were produced by two statistical models: a simple linear technique with recursive estimation of the model's parameters and a non-linear artificial neural network. A first experimental set-up that consists in using only past on-site irradiance measurements was proposed. It was shown that the proposed models clearly outperformed the reference models like persistence or smart persistence for forecasting horizon greater that 1h. Interestingly, it was found that a linear recursive technique like ARMA.RLS performed equally well or slightly better than a non-linear method such as an artificial neural network. This finding if confirmed may favor in the future the use of such linear recursive model for short-term solar forecasting. In addition, the recursive estimation of the model's parameters makes the method very well suited to online forecasting.
The second experimental set-up proposed the combination of day-ahead ECMWF forecasts with past ground measurements. It was demonstrated that the incorporation of additional exogenous inputs into the ANN model clearly improved the accuracy of the intra-day forecasts.
A future and interesting step could be the comparison of the proposed methods with satellite-based forecasts produced by the cloud motion vector (CMV) technique in the case of an insular context. Finally, preliminary results showed an impact of the sky conditions on the forecasting accuracy of the different methods. Insular sites like Le Tampon, which exhibits higher solar variability due to clouds that are formed locally, are prone to have the worse forecasting performance than less variable (continental or insular) sites.