L Mazorra Aguiar

B Pereira

P Lauret

F Díaz

M David

Combining solar irradiance measurements, satellite-derived data and a numerical weather prediction model to improve intra-day solar forecasting

Isolated power systems need to generate all the electricity demand with their own renewable resources. Among the latter, solar energy may account for a large share. However, solar energy is a fluctuating source and the island power grid could present an unstable behavior with a high solar penetration. Global Horizontal Solar Irradiance (GHI) forecasting is an important issue to increase solar energy pro duction into electric power system. This study is focused in hourly GHI forecasting from 1 to 6 h ahead. Several statistical models have been successfully tested in GHI forecasting, such us autoregressive (AR), autoregressive moving average (ARMA) and Artificial Neural Networks (ANN). In this paper, ANN models are designed to produce intra day solar forecasts using ground and exogenous data. Ground data were obtained from two measurement stations in Gran Canaria Island. In order to improve the results obtained with ground data, satellite GHI data (from Helioclim 3) as well as solar radiation and Total Cloud Cover forecasts provided by the European Centre for Medium Range Weather Forecasts (ECMWF) are used as additional inputs of the ANN model. It is shown that combining exogenous data (satellite and ECMWF forecasts) with ground data further improves the accuracy of the intra day forecasts.

Introduction

According to the trend shown for renewable energy generation, the proportion of these kinds of energies in the power system will increase in the next years. Renewable energies show a fluctuating generation profile because of their dependence on meteorological conditions. Since grid operators need to keep under control these variations, in order to accommodate the input/output balance of the system, forecasting methods are necessary to improve the coupling of renewable sources. This becomes an important issue when facing the fact that photovoltaic and solar thermo electrical plants are widely spread in most of the present power systems.

The output power of a PV plant is mainly correlated with the global solar irradiance received on the plan of the modules. The calculation of this received irradiance requires the knowledge of the direct and diffuse components of the solar irradiance. But in many cases as in our study, no measurement of these two com ponents exists. Thus, these two components are commonly derived from the global horizontal irradiance (GHI) with decomposition models [1e3]. When no assumption is done on the inclination of the PV modules, the key factor becomes the GHI. So, in this work, we will only focus on the forecasting of the GHI. Indeed, a reliable Global Horizontal Irradiance (GHI) forecasting model is considered as an important tool to avoid unstable behavior in the electrical grid, maintaining the balance between demand and supply [START_REF] Heinemann | Forecasting of Solar Radiation. Solar Energy Resource Managmente for Electricity Generation from Local Level to Global Scale[END_REF][START_REF] Wittmann | Case studies on the use of solar irradiance forecast for optimized operation strategies of solar thermal power plants[END_REF].

Isolated power systems, such as the one in Gran Canaria Island, need to generate all the electricity demand solely with their own resources. The island power grid could present an unstable behavior if there is a large scale solar energy introduction. This is aggravated when a great variability in solar radiation conditions are shown in the Canary Islands, caused by the diversity of the climatic areas that emerge as a by product of the differences in its orog raphy [START_REF] Díaz | An adaptive solar radiation numerical model[END_REF][START_REF] Diagne | Post-processing of solar irradiance forecasts from WRF model at Reunion Island[END_REF].

GHI forecasting has been developed in the last years using a wide range of methods. The most common forecasting models used in this field are:

Statistical models, models based on GHI time series prediction: most common statistical models used in GHI forecasting are linear models, such us autoregressive (AR) and autoregressive moving average (ARMA), and machine learning techniques, such as Artificial Neural Networks (ANNs) [START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF]. Forecasting using Ground based sky images: models based on sky images obtained with 180 vision angle cameras. Sky images lead to know cloud cover conditions for a few minutes ahead. Satellite images models. Geostationary satellites get atmosphere images all around the Earth with temporal resolution less than an hour. The high development occurred in the recent years in satellite data retrieving makes this technique a very useful tool to improve GHI forecasting [START_REF] Perez | Validation of short and medium therm operational solar radiation forecasts in the US[END_REF][START_REF] Hammer | Short-term forecasting of solar radiation: a statistical approach using satellite data[END_REF]. Numerical Weather Prediction models (NWP) based on physical models for estimating atmosphere conditions, including clouds formation and dissolution. Physical models are described with differential equations solved with numerical methods. NWPs models offer time horizons forecasting from few hours to 15 days ahead [START_REF] Heinemann | Forecasting of Solar Radiation. Solar Energy Resource Managmente for Electricity Generation from Local Level to Global Scale[END_REF][START_REF] Perez | Comparison of numerical weather prediction solar irradiance forecasts in the US, Canada and Europe[END_REF].

Depending on the forecasting time horizon and granularity, forecasting results are used for different purposes and are based on different input parameters [START_REF] Kostylev | Solar power forecasting performance toward industry standards[END_REF][START_REF] Sengupta | Best Practices Handbook for the Collection and Use of Solar Resource Data for Solar Energy Applications[END_REF]. Ground based sky images models show high precision information about cloud cover condi tions and movement for intra hourly forecasting horizons [START_REF] Lorenz | Prediction of solar irradiance and photovoltaic power[END_REF]. On the other hand, forecasting models based on cloud motion vectors derived from geostationary satellite images offer accurate results for hourly prediction range up to 6 h ahead [START_REF] Perez | Validation of short and medium therm operational solar radiation forecasts in the US[END_REF][START_REF] Hammer | Short-term forecasting of solar radiation: a statistical approach using satellite data[END_REF]. NWPs forecasting precision varies depending on the selected time scale and the geographic zone of the study. Heinemann [START_REF] Heinemann | Forecasting of Solar Radiation. Solar Energy Resource Managmente for Electricity Generation from Local Level to Global Scale[END_REF] showed that acquiring the radiation values for clear skies without deviation is possible. NWPs models associated with a post processing method using hourly ground measurements showed accurate results for predictions from 6 h onwards [START_REF] Diagne | Review of solar irradiance forecasting methods and a proposition for small-scale insular grids[END_REF].

The statistical method proposed in this paper is based on Arti ficial Neural Network (ANN) for GHI hourly forecasting from 1 to 6 h ahead. ANNs are considered a very attractive method for GHI forecasting because of their capacity to establish relationships be tween an input and output datasets [START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF]. ANNs have been used for solar radiation prediction using past ground measurements. Many papers have described ANN accurate results for GHI predictions with time horizons from few hours ahead, such as explained in Ref.

[16e19], to 24 h ahead daily irradiation forecasting [START_REF] Rehman | Artificial neural network estimation of global solar radiation using air temperature and relative humidity[END_REF][START_REF] Bosch | Daily solar irradiation estimation over a mountainous area using artificial neural networks[END_REF]. Statistical models based on ANN offer the possibility to combine historical GHI ground data with other meteorological data. Rehman [START_REF] Rehman | Artificial neural network estimation of global solar radiation using air temperature and relative humidity[END_REF] used air temperature and ground relative humidity for daily irradiation forecasting. Kemmoku [START_REF] Kemmoku | Daily insolation forecasting using a multi-stage neural network[END_REF] designed an ANN to forecast daily irradiation using atmospheric pressure and several meteo rological parameters. Sfetsos & Coonick [START_REF] Sfetsos | Univariate and multivariate forecasting of hourly solar radiation with artificial intelligence techniques[END_REF] used GHI ground data, air temperature, wind velocity and pressure for solar irradiance hourly forecasting. Moreover, several papers studied the influence of sunshine duration, air temperature, relative humidity, latitude and longitude to obtain GHI hourly and daily forecasting [START_REF] Mellit | A 24-h forecast of solar irradiance using artificial neural network: application for performance prediction of a grid-connected PV plant at Trieste, Italy[END_REF][START_REF] Ghanbarzadeh | Solar radiation forecasting based on meteorological data using artificial neural networks[END_REF][START_REF] Mohandes | Estimation of global solar radiation using artificial neural networks[END_REF].

ANNs may even lead to more accurate GHI forecasts by combining ground measurement time series analysis with exoge nous data, such as ground irradiance derived from satellite images and different NWPs data. In the last years, several studies used different techniques such as autoregressive models (AR) [26e28], artificial neural networks [START_REF] Marquez | Hybrid solar forecasting method uses satellite imaging and ground telemetry as inputs to ANNs[END_REF][START_REF] Aguiar | Use of satellite data to improve solar radiation forecasting with Bayesian artificial neural networks[END_REF] and genetic algorithms for choosing the relevant information from satellite derived data [START_REF] Zagouras | On the role of lagged exogenous variables and spatio-temporal correlations in improving the accuracy of solar forecasting methods[END_REF].

In this paper satellite derived and ground data are combined with solar radiation and total cloud cover predicted by the European Centre for Medium Range Weather Forecasts (ECMWF). The satellite derived data is selected as in Ref. [START_REF] Aguiar | Use of satellite data to improve solar radiation forecasting with Bayesian artificial neural networks[END_REF], using correla tion analysis results. In this work, we used the Bayesian framework [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF] which improves ANNs learning process. In addition, the Bayesian approach provides techniques to optimize the ANN ar chitecture and to select the most relevant inputs [START_REF] Lauret | Bayesian neural network approach to short time load forecasting[END_REF][START_REF] Saldana | Modelizaci on del recurso solar utilizando redes neuronales artificiales y su aplicaci on a la generaci on de mapas topogr aficos de radiaci on[END_REF].

The organization of this paper is as follows. In section 2, we present the ground data, satellite derived data and ECMWF data used in the present paper. The clear sky model is described in Section 3 and we analized Ground data distribution in Section 4. In Section 5, we detail the forecasting methods theory, in Section 6 we explain the ANN architecture optimization and Section 7 is devoted to the selection of the most relevant satellite inputs. In Section s8, we show the error metrics used in this paper and in Section 9 we detail the results of the different models and discuss the influence of the different exogenous data to improve GHI forecasting. Finally, Section 10 gives main remarks and conclusions.

Available data

The datasets used in the present study to build the forecasting model include ground, satellite and numerical predicted data. All the data will be converted into an hourly basis in order to compare and contrast all the information in a coherent manner. The tem poral availability for all type of data was restricted to the year 2005 because this is the only year in which ground and satellite data overlaps: ground data have been processed up to the year 2005, meanwhile the satellite data recollection started on the same year.

The application of the ANN models requires the division in training and testing datasets. Therefore, due to the restriction of only one year, we took the decision to select three weeks of a month for training and one week for testing purposes since this study only possess one year available. In this manner we can assure that training and testing datasets share similarities regarding meteo rological behaviors, while keeping a ratio of 75% training and 25% testing [START_REF] Zagouras | On the role of lagged exogenous variables and spatio-temporal correlations in improving the accuracy of solar forecasting methods[END_REF].

Ground data set

Ground data are obtained from two measurement stations located in Gran Canaria Island (Spain) and managed by the Canary Islands Technology Institute (Instituto Tecnol ogico de Canarias). These stations are part of a network of 23 stations in the Canary Islands, validated by the elaboration of a solar map of the zone [START_REF] Díaz | An adaptive solar radiation numerical model[END_REF][START_REF]Mapasolarcanarias[END_REF]. The selected stations are representative for the climatic di vision present on this island [START_REF] Díaz | An adaptive solar radiation numerical model[END_REF][START_REF] Mazorra | Typical meteorological year (TMY) evaluation for power generation in gran Canaria island, Spain[END_REF]. This climatic contrast between the northern region, with higher cloudiness, and the southern one, with clearer days and less variability, is well known and produced by the strong influence of the trade winds of the area.

The equipment used to acquire global solar horizontal irradi ance (GHI) is a secondary standard pyranometer, the CMP 11 of Kipp & Zonen, with 3% accuracy for daily sum of GHI. Data are recorded for every 5 s, 1 min average and later assembled into an hourly basis. In order to maintain enough quality of the data, all the measurement data were treated with the SERI QC control software [START_REF] Maxwell | Users Manual for SERI QC Software, Assessing the Quality of Solar Radiation Data[END_REF][START_REF] Younes | Clear-sky classification procedures and models using a world-wide data-basee[END_REF]. This treatment consists on filtering values of radiation that are negative or higher than the top of the atmosphere.

Furthermore, in global solar radiation forecasting, it is usual to remove night data and anomalous readings from the equipment at dawn and sunset, or other shadowing effects. Therefore, a filter based on the solar zenith angle (SZA) was introduced to establish a frontier of 80 for all data. This threshold is validated in the bibli ography [START_REF] Blanc | Atlas du Potentiel Solaire Photovoltaique et Thermodynamique en region PACA, Centre Energetique et Procedes[END_REF], and also concurs with the careful observations of the shadow effect for each station.

Satellite derived data set

Satellite information used is this work was retrieved from the Helioclim 3 database version 5 (HC3v5). All this information has been processed by the Heliosat 2 method, where the images taken from the Meteosat geostationary satellite network are converted into solar radiation in the area of study [START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF][START_REF] Blanc | The HelioClim project: surface solar irradiance data for climate applications[END_REF].

The selected area contains the entire island of Gran Canaria as well as a significant portion of sea at the north east, motivated by the knowledge and influence of the trade winds in the Canary Islands (see Fig. 1). This area is defined, in decimal degrees, by the coordinates latitude [þ28.7500 to þ27.2500], and longitude [ 16.0000 to 14.5000], resulting in a grid of 61 Â 55 pixels of information, where each pixel possesses a spatial resolution of 3 Â 3 km 2 . The information provided consists of global horizontal irradiation and irradiation on top of the atmosphere (ToA). All of these parameters are handled in Wh/m 2 with a temporal resolution of 15 min, assembling then into an hourly manner in order to be in accordance with the hourly ground data.

An initial assessment of the satellite information shows that the relative Root Mean Square Error between these data and the ground data is coherent. These estimations have an average 12.2% rRMSE at C0_Pozo Izquierdo and 27.8% rRMSE at C1 Las Palmas. Further assessment on the satellite data values indicates a consis tency with the climatic behavior of the island, following a tendency very similar to the ground data values [START_REF] Aguiar | Use of satellite data to improve solar radiation forecasting with Bayesian artificial neural networks[END_REF].

ECMWF data set

In this work, we used the predicted data provided by the ECMWF (European Centre for Medium Range Weather Forecast), granted by the Laboratoire de Physique et Ing enierie Math ematique pour l'Energie et l'environnement (PIMENT) from the Universit e de La R eunion.

This information was selected for the island of Gran Canaria for an area located from 27.5 to 28.5 of latitude north, and 15 to 16 longitude west Fig. 2. For each pixel of the selected grid there is a spatial resolution of 16 Â 16 km 2 and 21vertical levels of resolution. On the order hand, all the data concerning the ECMWF for the year 2005 comes within 3 h intervals, therefore an interpolation of the values into an hourly basis was necessary to compare data sets. This was achieved using a method that takes into account the conser vation of the solar energy [START_REF] Espinar | Report on the Harmonization and Qualification of Meteorological Data: Project ENDORSE, Energy Downstream Service, Providing Energy Components for GMES Grant Agreement[END_REF].

This ECMWF extraction provides a lot of information for a great number of meteorological variables, but in this study we only extracted the following variables described by Latitude, Longitude and Time:

Total Cloud Cover (TCC), with values between 0 and 1 using a cloud index, Surface Solar Radiation Downwards (SSRD), for accumulative values of J/m 2 within two instants.

Data transformation: use of the clear sky index

In order to work with statistical models, all datasets have to be transformed into stationary time series, which mean that auto correlation structures are constant over time [START_REF] Chatfield | Time Series Analysis, an Introduction[END_REF]. The approach used in this work was the Clear Sky index (k*), extensively applied in the bibliography [START_REF] Sengupta | Best Practices Handbook for the Collection and Use of Solar Resource Data for Solar Energy Applications[END_REF][START_REF] Diagne | Review of solar irradiance forecasting methods and a proposition for small-scale insular grids[END_REF], in order to remove the seasonal and daily trend from the hourly time series. This index is obtained by applying the following formula:

k à I g I cs (1)
where I g is the hourly measured GHI and I cs is the hourly irradiance computed by a clear sky model for a determined location and time. Many Clear Sky models have been proposed depending on the different climatic parameters used as inputs [START_REF] Reno | Global Horizontal Irradiance Clear Sky Models: Implementation and Analysis[END_REF][START_REF] Younes | Quality control of solar radiation data: present status and proposed new approaches[END_REF]. This survey is based on the Bird model [START_REF] Bird | Simplified the Clear Sky Moel for Direct and Diffuse Insolation on Horizontal Surfaces[END_REF], which is well known to provide ac curate results with only a few meteorological inputs [START_REF] Badescu | Accuracy analysis for fifty-four clear-sky solar radiation models using routine hourly global irradiance measurements in Romania[END_REF]: Aerosol Optical Depths (AOD500 mm and AOD380 mm), water vapour and Ozone atmospheric content. AODs and water vapour column were obtained from AERONET network [START_REF] Holben | AERONET a federated instrument network and data archive for aerosol characterization[END_REF]. We obtained data from 2008 to 2014 in order to calculate a climatological annual mean of the atmospheric parameters for all the years. On the order hand, the Ozone was retrieved from the World Ozone Monitoring Mapping provided by the Canadian Government [START_REF]Canada's Environment, World Ozone Monitoring Mapping[END_REF] and also aggregated in an annual mean.

An evaluation of the Bird Clear Sky model was made for two measurement stations in Gran Canaria Island: C0 Pozo Izquierdo [START_REF] Ineichen | Comparison of eight clear sky broadband models against 16 independent data banks[END_REF], allowing the detection of clear sky hours from the ground dataset in order to compare them with the values given by the clear sky model. This accuracy is quantified by the relative Root Mean Square Error (% rRMSE). The results of the Bird method for the selected locations are listed in Table 1.

The model accuracy for the number of clear sky hours found for each site is around 4% rRMSE. These errors are in the same order as those obtained through the McClear model, which uses AOD, Ozone and water vapour from MACC project [START_REF] Lefevre | McClear: a new model estimating downwelling solar radiation at ground level in clear-sky conditions[END_REF]. This means that the variability of these atmospheric input parameters can be consid ered as negligible for the case of Gran Canaria and climatological annual mean values are enough.

This result does not indicate any difference between southern and northern stations, while the number of clear sky hours is clearly higher in Pozo Izquierdo (southern), which is consistent with the behavior in the climate of the island. On the other hand, the variability of each site is given by the method proposed by Hoff and Perez [START_REF] Hoff | Quantifying PV power output variability[END_REF], with the standard deviation of the change in the clear sky index, were southern station also takes a lower value compared with northern one, maintaining the coherence with the known climate patterns in the area of study.

Ground data analysis

In order to characterize ground data of each location we decided to classify days in all the available data by mean of intraday clear sky index average, Eq (2), and intraday variability, Eq (3) and ( 4), as explained in Ref. [START_REF] Dambreville | Prevision du rayonnement solaire global par teledection pur la gestion de la production d'energie photovoltaique[END_REF][START_REF] Lauret | Characterization of the intraday variability regime of solar irradiation of climatically distinct locations[END_REF]. With these two parameters we could clearly classify the days in nine different frames, see Fig. 3. We divided intraday means into A, B and C days, where A is the heavy clouded type and C data type is the clearer one. In the other hand, the classification by roman numbers I, II and III follow the vari ability, where a I day stands for a stable variability and a type III describes a great variability through the day. By this distribution an example of a stable and clear day will fit into a type CIII.

k à P H d t 1 k * ðtÞ H d (2) Dk à ðtÞ k à ðt þ 1Þ k à ðtÞ (3) stdðDk Ã Þ P H d t 1 Dk à ðtÞ Dk à 2 H d v u u u u t ( 4
)
All the data represented in Fig. 3 and on Table 2 show great coherence with the empirical observations of the weather of Gran Canaria Island. There are clearer days in the southern part (C0), with 166 days in type C, than the northern one (C1), with 234 falling in the type B class.

Forecasting models description

The solar radiation parameter used by the statistical models proposed in this work is the clear sky index. Indeed, the general function used to relate input and output datasets is as follows:

c k * ðt þ hÞ F k * G ðtÞ
Where c k * ðt þ hÞ is the clear sky index predicted for time horizon h, from h 1 to h 6 h with hourly granularity, k * G ðt iÞ is clear sky index from the ground dataset at the location of interest for 'i' past values and k * Exn ðt jÞ corresponds to the different exogenous clear sky indices (inferred from satellite or ECMWF data) for past values. One the most important decision to take into account by the modeler is the number of inputs for the model. The model may obtain less accurate results if the ANN architecture is unnecessarily complex because of irrelevant inputs.

The dataset for the training process consists in pairs of input output vectors with 'p' number of patterns, P fX i ; y i g p i . The gen eral function F is established during the training process. In this work, X i matrix contains ground radiation data, satellite derived radiation data and ECMWF predicted parameters, while y i vector contains ground solar irradiance for the forecasting time horizon. Clear sky index is the parameter used during the training and testing process, however, results discussed in section 6 will be in terms of global solar irradiances, obtained from Eq. ( 1) as b I g c k à $I cs .

Naive models

In the literature, it is common practice to compare results ob tained by the proposed models against some naïve or reference models in order to discuss the improvement. In this study, we have worked with one naïve model and a climatological one. The naïve model is the so called Smart Persistence (Smart Pers) model [START_REF] Perez | A new operational solar resource forecast model service for PV fleet simulation[END_REF]. The latter only works with ground historical radiation data from the location of interest, using clear sky index series. It consists in forecasting clear sky index for each time horizon 'h' as the mean of the 'h' previous clear sky values [START_REF] Hoff | Modeling PV fleet output variability[END_REF]. Smart persistence model is defined by Eq. ( 6).

c k * À t þ h Á mean k * ðtÞ; …; k * À t h Á (6)

Climatological model

The climatological model [START_REF] Lorenz | Prediction of solar irradiance and photovoltaic power[END_REF], which is based on the historical mean (see Eq. ( 7)), provides a constant forecasting whatever the horizon is:

c k * ðt þ hÞ mean À k * ðtÞ Á historical mean (7)

Forecasting ANN based model

Artificial Neural Networks (ANNs) are statistical methods capable to establish a relation between an input and output data sets. ANNs are composed by several individual units, called neu rons, and connected by weights. These neurons received information from an external parameter or another neuron. During the training process the associated weights are modified to obtain the optimal weight distribution. Usually, each neuron applies an activation function to calculate the correspondent output. In this paper the function applied is the hyperbolic tangent function:

f ðxÞ tanhðxÞ e x e Àx e x þ e Àx (8)
One of the most common ANN model used in engineering ap plications is based on the Multilayer Perceptron (MLP). The MLP consists in one layer of inputs, one hidden layer made of several non linear neurons and one output layer with no feedback or lateral connections between them. The output layer is made of one neuron with an identity activation function (see Eq. ( 9)).

c k * ðt þ hÞ a 0 þ X h j 1 a j $f " b 0;j þ X p i 1 b i;j $ k * ðt iÞ # (9)
The optimal weights are computed optimizing a cost function. The most common function is the average of square difference between the network output c k * ðt þ hÞ and the desired output, target set 'k*' (see Eq. ( 10)). Firstly, a random weight set is estab lished and then, during the training process, are optimized by minimizing a cost function by backpropagation algorithm. The backpropagation training algorithm used to minimize E(w) is the scale conjugate gradient.

EðwÞ

1 2

X N i 1 c k * ðt þ hÞ k * i ! 2 (10)
To improve the accuracy of the model, one of the most impor tant issues is the network architecture. Too complex ANNs may lead to a good approximation of the training dataset, but poor results when a new test dataset is used (this problem is called overfitting). It is generally recommended to control the model complexity in order to achieve good generalization results [START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF][START_REF] Lauret | Bayesian neural network approach to short time load forecasting[END_REF]. Several model complexity control techniques exist and in this paper a Bayesian regularization framework is proposed [START_REF] Lauret | Bayesian and Sensitivity analysis approaches to modeling the direct solar irradiance[END_REF][START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF][START_REF] Lauret | A node pruning algorithm based on the Fourier amplitude sensitivity test method[END_REF]. The Bayesian approach considers a probability density function over the weight space and introduces new hyperparameters a and b to control model complexity. The optimal values correspond to the maximum of the probability density function optimizing the following cost function, Eq. [START_REF] Perez | Comparison of numerical weather prediction solar irradiance forecasts in the US, Canada and Europe[END_REF].

SðwÞ b 2 X N i 1 c k * ðt þ hÞ k * i ! 2 þ a 2 X m i 1 w 2 i (11)
As described in Lauret, 2006 [START_REF] Lauret | Bayesian and Sensitivity analysis approaches to modeling the direct solar irradiance[END_REF], to find the optimal hyper parameters a and b and the optimal weight vector w MP we followed an iterative procedure. At the beginning, weights are randomly initialized using a Gaussian distribution and the hyperparameters use small values. The optimal weights sets are computed using the standard scale conjugant gradient algorithm, minimizing S(w).

Once we obtained the first optimized weights set, hyperparameters are obtained using the Bayesian Framework numerical imple mentation. These steps are repeated until the regularized error is equal to half the number of data points, S(w) N/2.

ANN architecture optimization

In this study, the NETLAB toolbox developed in Matlab was used for ANNs computation [START_REF] Nabney | NETLAB: Algorithms for Pattern Recognition[END_REF]. We focused on improving ANNs fore casting results combining ground data and exogenous data as inputs.

ANN model complexity decision is one of the most important problems to study in ANNs forecasting. In this paper, we focus on number of inputs and hidden units. Techniques such as Master optimization [START_REF] Cornaro | Master optimization process based on neural networks ensemble for 24-h solar irradiance forecast[END_REF] or genetic algorithm [START_REF] Marquez | Hybrid solar forecasting method uses satellite imaging and ground telemetry as inputs to ANNs[END_REF] can be used to optimize the ANN models. However, in this work, based on our experience [START_REF] Lauret | Bayesian neural network approach to short time load forecasting[END_REF], we chose to use the Bayesian technique.

The Bayesian framework also offers different possibilities to study the model architecture. Automatic Relevance Determination (ARD) technique assigns a different regularization coefficient to each input and allows determining the most relevant inputs during the training process [START_REF] Lauret | Bayesian neural network approach to short time load forecasting[END_REF][START_REF] Saldana | Modelizaci on del recurso solar utilizando redes neuronales artificiales y su aplicaci on a la generaci on de mapas topogr aficos de radiaci on[END_REF]. ARD information of six inputs was studied in order to choose the most relevant for GHI forecasting. The best result for all time horizons and both stations was obtained using 6 inputs of ground measurements inputs. Fig. 4 shows an example with a high relevance for all six inputs in C0 station.

Once the number of ground lags inputs is fixed, Bayesian framework is used for choosing the number of hidden units by computing the probability of each model (evidence of the model) [START_REF] Aguiar | Use of satellite data to improve solar radiation forecasting with Bayesian artificial neural networks[END_REF]. The highest evidence of the model is used to select the best option for each station and time horizon [START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF][START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF][START_REF] Penny | Bayesian Neural Networks for classification: how useful is the evidence framework?[END_REF]. Most of results show low number units as we can see in Fig. 5 for station C0 and time horizon 3 h. However the difference between %rRMSE and % rMAE for the testing dataset in each experiment were not so pronounced.

To avoid overfitting we used %rRMSE and %rMAE obtained on test datasets to decide the final result. Bayesian framework has been widely used and tested for forecasting [START_REF] Lauret | Bayesian neural network approach to short time load forecasting[END_REF] obtaining very good results to avoid overfitting. In this case, Fig. 6 shows the final result obtained for training and testing datasets forecasting using classic ANN training and Fig. 7 shows the results obtained with Bayesian Framework, where overfitting problem is not present.

Selection of the satellite data inputs

ANNs have been trained using different combination of input data. Every ANNs forecast hourly GHI with six past ground clear sky data as inputs. The exogenous data added to the ground clear sky data are solar radiation data derived from Helioclim 3 and solar radiation and total cloud cover (TCC) predictions from ECMWF. Finally, in addition to the Smart persistence and climatology models, the ANN models computed and compared in this paper are the following: ANN model with only past ground data as inputs (denoted herein NN). ANN model with past ground data and satellite data as inputs (denoted herein NN þ SAT). ANN model with ground data and ECMWF radiation data as inputs (denoted herein NN þ ECMWF). ANN model with ground data, satellite data and ECMWF radia tion data as inputs (denoted herein NN þ ECMWF þ SAT).

One of the most important decisions is to select the best satellite information to improve GHI hourly forecasting. We considered a maximum number of 30 satellite derived radiation data based on the calculation explained by Mazorra et al. [START_REF] Aguiar | Use of satellite data to improve solar radiation forecasting with Bayesian artificial neural networks[END_REF]. In order to study weather conditions in Gran Canaria Island and establish the best satellite pixels between whole gridded data at the different sta tions, we computed the Pearson correlation between satellite and ground data [START_REF] Dambreville | Very short term forecasting of the Global Horizontal Irradiance using a spatio-temporal autoregressive model[END_REF][START_REF] Zagouras | On the role of lagged exogenous variables and spatio-temporal correlations in improving the accuracy of solar forecasting methods[END_REF]. The huge amount of satellite derived data makes the computation difficult, so we applied a median filter for each 3 Â 3 satellite pixels. Consequently, a Super pixel was created computing GHI median value of every 3 Â 3 group of pixels, Fig. 8.

In each station the correlation between ground data at t 0 and each satellite derived pixel with time lags from the same temporal period to a maximum of 3 h earlier, Eq. ( 12), was computed. In that way we establish a relation between gridded pixels for different time lags and the ground data for the present moment. Therefore, we have information about incoming events from the surroundings. C k * ði; jÞ h corr k à ground ðtÞ; k à satellite ðt hÞ for h 0; 1; 2& 3 [START_REF] Kostylev | Solar power forecasting performance toward industry standards[END_REF] Thus, ANNs inputs represent the surrounding pixels with highest relate with ground station data for different time lags. As explained in Ref. [START_REF] Aguiar | Use of satellite data to improve solar radiation forecasting with Bayesian artificial neural networks[END_REF], the distribution of pixels in the different time lagged images is according to the number of correlation values over 0.5 in each station [START_REF] Zagouras | On the role of lagged exogenous variables and spatio-temporal correlations in improving the accuracy of solar forecasting methods[END_REF]. To decide the best selection of pixel, we computed six different tests during the training process. As explained in Ref. [START_REF] Aguiar | Use of satellite data to improve solar radiation forecasting with Bayesian artificial neural networks[END_REF], the distribution of pixels in the different time lagged images is according to the number of correlation values over 0.5 in each station [START_REF] Zagouras | On the role of lagged exogenous variables and spatio-temporal correlations in improving the accuracy of solar forecasting methods[END_REF]. For each station, four tests were calculated using the annual correlation results and two more test used a different distribution every quarterly group of images. The best was selected using the lowest %rRMSE and %rMAE for the testing dataset.

Performance error metrics

In order to evaluate the performance of each method, we used two standard error metrics widely used in the solar forecasting community: the Root Mean Square Error (RMSE) and the Mean Absolute Error (MAE) [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF][START_REF] Willmott | Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance[END_REF]. Dividing both absolute error by the average of the hourly GHI data we compute their relative metrics (% rRMSE and %rMAE). In this paper, to study the quality of the fore casting methods, we provide the relative errors. To compare the relative improvement of the different models with respect to the persistence simple model, we calculate the forecast skill parameter [START_REF] Coimbra | Overview of solar-forecasting methods and a metric for accuracy evaluation[END_REF].

RMSE

1 N

X N i 1 GHI forecast;i GHI measured;i 2 v u u t (13) MAE 1 N X N i 1 GHI forecast;i GHI measured;i (14)
SKILLð%Þ 1 RMSE method RMSE persistence ! Â 100 (15)

Forecasting results

Once we decided the optimal satellite derived pixels, the different ANN performance results were compared. Tables 3e8 present the forecast performance for both station in terms of % rRMSE, RMSE and %SKILL from time horizon 1 to 6 using annual testing data set. Each row corresponds to the forecasting methods used in this paper. In all cases, forecasting models show worse % rRMSE error as time horizon increases. The Smart Persistence in creases from 15%rRMSE in 'h 1' to 26%rRMSE in 'h 6' at C0 station, while at C1 station error results oscillate between 27% and 42%. On the other hand, all models based on ANNs architecture lead to error results around 22% at C0 and 34% at C1 for time horizon 'h 6'. Different error results, almost 12% difference for 'h 6', have been observed between northern and southern stations. This result is consistence with cloud formation processes explained and observed in the north. Furthermore, it is also observed a dichotomy between the improvements of non linear methods depending on the station and sky conditions. At station C0, ANN models show an improvement between 1.5% and 4% with Smart Persistence model for 'h 1' and 'h 6'. While at C1, ANN models lead to a gain from 3% to 6%. Models with exogenous data improve on ANN models based only on ground data (NN) for all cases and at both stations. Figs. 9 and 10 show the influence of satellite radiation data and ECMWF radiation and TCC data in terms of %rRMSE for both stations.

At station C0, NN þ SAT obtain the best results for the first 3 h ahead compared to NN þ ECMWF (improvement over 1%). This results are consistent with the fact that satellite data comes from time lags from 'h 0' to 'h 3' hours. In the last 3 h, NN þ ECMWF is similar to the satellite data (NN þ SAT). The combination of both exogenous data, NN þ ECMWF þ SAT, leads to the best error results in terms of %rRMSE in general.

At station C1, as well as C0, satellite data achieve better results for the first 3 h, while ECMWF improve almost 2% sat ellite data for 'h 6'. The combination of both exogenous data (Satellite and ECMWF) leads to the best error results for all time horizons.

Tables 5 and8 present the %SKILL for all time horizons in both stations. Even if in southern stations, the results with ANN þ ECMWF þ SAT are much better than in the north, the SKILL parameter shows very similar results for both areas. Furthermore, the SKILL forecast increases with time horizon, which means that the more far ahead in time, the better results we get with ANN þ ECMWF þ SAT method compared with persistence model. Moreover, Figs. 11 and12 show the results for both stations in terms of %SKILL. It is easily observed that in station C0 the satellite data offer the best results for the first 3 h, while both models obtained very similar results for the last 3 h. In station C1, NN þ SAT model obtain the best results from 'h 1' to 'h 3' while NN þ ECMWF is the best model from 'h 4' to 'h 6'. For both stations, overall, the As a consequence, the model recommended for the whole Gran Canaria Island is NN þ ECMWF þ SAT. This model improves NN model (only with ground data) from 1% for 'h 1' to almost 2% for 'h 6' at station C0 and around 1.5% for 'h 1' to almost 3.5% for 'h 6' at station C1.

Conclusions

This work proposes to combine ground measurements with exogenous inputs provided by satellite and NWP data in order to improve intra day solar forecasting.

In order to explain the influence of each exogenous data on forecasts, satellite data and data from ECMWF were first studied separately and in combination. At Station C0, for the first 3 h forecast, inputting satellite data (NN þ SAT) obtained better results compared to inputting data from ECMWF (NN þ ECMWF). How ever, in the last 3 h, similar results were obtained by including ECMWF data. In this case, the best results were obtained by combining both models (NN þ ECMWF þ SAT) as compared with NN. In Station C1, this same effect is more clearly observable. In the first 3 h, the best model is the one that uses satellite data (NN þ SAT) and in the last hours the best model is the ECMWF data model (NN þ ECMWF). Similarly, the best results were obtained by combining these models (NN þ ECMWF þ SAT). The conclusion is that, in general combining exogenous satellite data and ECMWF data provides the best forecast results for both stations.

We would recommend working with the NN þ ECMWF þ SAT model on the entire island. With this model, at the C0 Pozo Izquierdo station, %rRMSE increases with the forecast horizon from 15.47% for 'h 1' to 22.17% for 'h 6', and for C1 Las Palmas, it fluctuated between 24.15% and 34.09%.

More work has to be done in order to cope with all meteoro logical situations. For instance, with more training data or by implementing a specific procedure (random sampling of the training cases), one can create a model ensemble that may be able to cope with the different meteorological conditions. Lastly, future lines of research should carry out a more effective selection of satellite pixels. Similarly, wind and humidity would give us very important information about solar radiation behavior. We propose carrying out a detailed study of wind and humidity in the various layers of the atmosphere in order to determine which layer offers the best results and relates these information with satellite pixels.

Fig. 1 .

 1 Fig. 1. Geographic distribution of satellite-derived data obtained from HC3v5 for Gran Canaria Island. Courtesy of MINES ParisTech/ARMINES.

Fig. 2 .

 2 Fig. 2. ECMWF grid for the domain.

[ 27 .

 27 8175 ºN, 15.4244 ºE, 47 m] and C1 Las Palmas [28.1108 ºN, 15.4269 ºE, 17 m]. The assessment of the accuracy of the Bird clear sky model was calculated with the Ineichen method

Fig. 3 .

 3 Fig. 3. Distribution of days in the plane k à À stdðDk Ã Þ of ground dataset for station C1Las Palmas (up) and station C0-Pozo Izquierdo (down). Intraday mean is divided in A (0 < k à < 0.5), B (0.5 < k à < 0.9) and C (0.9 < k à ) days. Daytype variability is divided in I (0< stdðDk à Þ<0.05) great variability, II (0.05< stdðDk à Þ<0.15) and III (0.15< stdðDk à Þ) a steady one.

Fig. 4 .

 4 Fig. 4. ARD results using 6 h ground measurement lag data as inputs in ANNs at station C0 Pozo Izquierdo for time horizon h 3 h.

Fig. 5 .

 5 Fig. 5. Logev information using 6 h ground measurement lag data and satellite data as inputs in ANNS at station C0 Pozo Izquierdo for time horizon h 3 h.

Fig. 6 .

 6 Fig. 6. Measured data vs Forecasted data for the training (left) and testing (right) datasets using classic ANN for station C0.

Fig. 7 .

 7 Fig. 7. Measured data vs Forecasted data for the training (left) and testing (right) datasets using Bayesian Framework for station C0.

Fig. 8 .

 8 Fig. 8. Superpixel (3 Â 3) selection at station C1-Las Palmas using Test-1 for time-lagged correlation images, t 0, 1, 2 & 3 h. Black area shows selected superpixels. Each pixel represents the intercorrelation with k * between ground measurement and each satellite pixel.

Fig. 9 .

 9 Fig. 9. %rRMSE for testing dataset at station C0 Pozo Izquierdo for different time horizons using ANN models with satellite and ECMWF radiation data.

Fig. 10 .

 10 Fig. 10. %rRMSE for testing dataset at station C1-Las Palmas for different time horizons using ANN models with satellite and ECMWF radiation data.

Fig. 11 .

 11 Fig. 11. %SKILL for testing dataset at station C0 Pozo Izquierdo for different time horizons using ANN models with satellite and ECMWF radiation data.

Fig. 12 .

 12 Fig. 12. %SKILL for testing dataset at station C1-Las Palmas for different time horizons using ANN models with satellite and ECMWF radiation data.

Table 1

 1 Measurement stations characteristics, Bird model inputs and performance.

	C0 Pozo Izquierdo	C1 Las Palmas

Table 2

 2 Distribution of day type categories in Gran Canaria.

		C0 Pozo Izquierdo	C1 Las Palmas
	Total A (0 < k * < 0.5)	27	7.78%	58	16.86%
	A I (0< stdðDk * Þ<0.05)	3		2
	A II (0.05< stdðDk * Þ<0.15)	9		2 6
	A III (0.15< stdðDk * Þ)	15		30
	Total B (0.5 < k * < 0.9)	154	44.38%	234	68.02%
	B I (0< stdðDk * Þ<0.05)	3		1
	B II (0.05< stdðDk * Þ<0.15)	62		78
	B III (0.15< stdðDk * Þ)	89		155
	Total C (0.9 < s)	166	47.84%	52	15.12%
	C I (0< stdðDk * Þ<0.05)	71		13
	C II (0.05< stdðDk * Þ<0.15)	83		32
	C III (0.15< stdðDk

* Þ) 12

Table 3

 3 RMSE for testing dataset at station C0 Pozo Izquierdo for different time horizons (in bold the best result).

	Model	1 h	2 h	3 h	4 h	5 h	6 h
	SMART-PERS	92.47 124.64 140.10 144.44 141.50 138.69
	CLI	145.13 145.17 145.25 145.33 145.42 145.50
	NN	88.20 113.39 125.13 127.12 131.68 130.03
	NN þ ECMWF	86.40 112.93 117.70 116.94 120.01 120.88
	NN þ SAT	83.58 104.58 113.51 116.15 117.31 120.60
	NN þ ECMWF þ SAT	84.00 106.17 110.51 114.93 118.89 120.43

Acknowledgments

This work has been supported by the "Convocatoria 2014 e Proyectos I þ D þ I e Programa Estatal de Investigaci on, Desarrollo e Innovaci on orientada a los retos de la sociedad" of the Spanish Government, grant contracts CTM2014 55014 C3 1 R, as part of the project "Integraci on de nuevas metodologías en simulaci on de campos de vientos, radiaci on solar y calidad del aire" and by the "C atedra Endesa Red de la Universidad de Las Palmas de Gran Canaria en su Convocatoria del curso 2014 2015". The autors are very grateful to Canary Islands Technology Institute (Instituto Tecnol ogico de Canarias, I.T.C.) for providing ground measurements databases and to Dr. Philippe Blanc, Responsable des activit es de recherche sur l' evaluation des ressources energ etiques renouvelables of MINES Par isTech/ARMINES, for his help with satellite databases.