Cyril Voyant
email: cyril.voyant@ch-castelluccio.fr

Pierrick Haurant

Marc Muselli

Christophe Paoli

Marie-Laure Laure Nivet

Time series modeling and large scale global solar radiation forecasting from geostationary satellites data

Keywords: Time series, Artificial neural networks, Irradiance, Prediction, Gamma index. Prediction error (measurement-prediction) and Distance (also called DTA) and intensity Tolerances

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The production and use of non-renewable resources based on fossil fuels combustion are responsible of real public health problem and raise environmental concerns. There are lots of alternatives such as photovoltaic and wind energy sources, which one of the main advantages are the renewable and inexhaustible aspects and the main disadvantage is related to their intermittencies [START_REF] Hocaoglu | Stochastic approach for daily solar radiation modeling[END_REF][START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF]. These non-continuities can cause a demand/production unbalance involving irrelevant wind or solar systems uses. To overcome this problem, it is necessary to predict the resource and to manage the transition between different energies sources [START_REF] Bouhouras | Installation of PV systems in Greece-Reliability improvement in the transmission and distribution system[END_REF][START_REF] Darras | PV output power fluctuations smoothing: The MYRTE platform experience[END_REF]Muselli et al., 1998b). Considering the grid manager's point of view [START_REF] Köpken | Assimilation of Meteosat radiance data within the 4D-Var system at ECMWF: Assimilation experiments and forecast impact[END_REF], needs in terms of prediction can be distinguished according to the considered horizon: following days, next day by hourly step, next hour and next few minutes. We choose in this paper to focus only on h+1 horizon prediction of global radiation as a first step (one hour in advance).

Of course, we are aware concerning the importance of other horizons (Voyant et al., 2013b). Note that it is appreciable to know the eventual fluctuations at least 30 minutes ahead (ignition delay of turbines)

for an ideal electrical grid management [START_REF] Troccoli | Seasonal climate forecasting[END_REF]. With efficient prediction tools dedicated to grid managers, the PV part in the mix energy would be increased; actually in France, the intermittent energy contribution is limited to 30%.

Several prediction methods have been developed by experts and can be divided in three main groups: methods using mathematical formalism of times series (TS) [START_REF] De Gooijer | 25 years of time series forecasting[END_REF][START_REF] Elminir | Prediction of hourly and daily diffuse fraction using neural network, as compared to linear regression models[END_REF][START_REF] Hamilton | Disaggregation of satellite derived irradiance maps: Evaluation of the process and application to Corsica[END_REF], numerical weather predictions (NWP) and models based on clouds detection [START_REF] Inness | NWP Modelsthe Basic Principles[END_REF][START_REF] Perez | Chapter 2 -Semi-Empirical Satellite Models[END_REF]. In this study, we have chosen to study prediction methods of the first group and we will study if this methodology can be an alternative to the NWP models. Not that in the literature, the NWP models are compared against ground measurements and the error established is approximately 30-40% but depends on the orography and micro-climate studied [START_REF] Paulescu | Weather modeling and forecasting of PV systems operation[END_REF]. The time series formalism (TS) and modelling is often used in 1 dimensional

(1-D) global radiation predictions, i.e. related to one measurement system at ground level [START_REF] Pons | Mapping a topographic global solar radiation model implemented in a GIS and refined with ground data[END_REF]. Persistence, autoregressive models, multilayer perceptron (MLP) and more widely artificial neural network (ANN) often applied to this aim [START_REF] Hocaoglu | Stochastic approach for daily solar radiation modeling[END_REF][START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF][START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF].

In this paper, we will complete the first prediction results exposed in [START_REF] Haurant | Hourly global radiation prediction from geostationary satellite data[END_REF] and we will show that the TS formalism applied to surface solar irradiances (SSI) estimations reduced to one point (1-D approach) can be generalized in the 2-D case, even if no ground detector is present.

Alternative approaches are available in [START_REF] Loyola | Applications of neural network methods to the processing of earth observation satellite data[END_REF][START_REF] Rahimikhoob | Comparative study of statistical and artificial neural network's methodologies for deriving global solar radiation from NOAA satellite images[END_REF]. From this point of view, satellite derived SSI maps extracted from the HelioClim-3 (HC-3) [START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF] database and centered on Corsica are used as hourly 2-D data generator. Each data series is processed with stochastic estimators in order to generate 1158 predictions per hour (1158 pixels per map separately treated). In fact, for an overall year of prediction, it is necessary to generate more than 10 million of hourly predictions (24x365x1158). The purpose of this paper is to generate one hour in advance predicted global radiation maps of a specific area from HC-3 SSI maps. But as data used are available for the entire surface, the method can be easily generalized. The geographical effects are taken into account using clear sky index in addition to temporal or seasonal phenomena [START_REF] Allan | Combining satellite data and models to estimate cloud radiative effect at the surface and in the atmosphere[END_REF].

The uncertainty of the used satellite derived SSI maps is about 16-23% (http://www.soda-is.com/eng/helioclim/helioclim3_uncertainty_eng.html). Moreover, in this paper, inputs of stochastic models are previously measured values, however, in a previous paper we have shown that exogenous data (parameters such as temperature, air pressure etc.) improve the prediction efficiency (Voyant et al., 2013b). As it was our first experiment in using geostationary satellites data, we have preferred to start without the multivariate case.

In the next section, the Meteosat images acquisitions, the SSI computation by Heliosat-II model described in [START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF][START_REF] Gueymard | Clear-sky irradiance predictions for solar resource mapping and large-scale applications: Improved validation methodology and detailed performance analysis of 18 broadband radiative models[END_REF] and clear sky index computing methodologies [START_REF] Maini | Weather Satellites[END_REF] will be first explained. Then we will detail the methodologies of prediction we have tested, taking care to explain first the TS formalism dedicated to the 2-D global solar radiation modeling. Then we will expose results of evaluation between modeling and HC-3 SSI over 2 years allowing the cross-comparison of models. Finally we will close the paper with a conclusion.

Materials and methods

Before to expose the time series modeling in the 2-D case, two brief sections are dedicated to the global radiation estimations from satellite acquisitions, the determination of clouds and the clear sky index described in [START_REF] Şenkal | Modeling of solar radiation using remote sensing and artificial neural network in Turkey[END_REF][START_REF] Singh | A comparison of the performance of Kalpana and Meteosat-7 WV radiances in the ECMWF NWP model[END_REF][START_REF] Thies | Satellite based remote sensing of weather and climate: recent achievements and future perspectives[END_REF]. The clear sky index based on the European Solar Radiation Atlas (ESRA) model [START_REF] Rigollier | On the clear sky model of the ESRA -European Solar Radiation Atlaswith respect to the heliosat method[END_REF] has the characteristic to define a TS made stationary [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF][START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF] and therefore directly usable by most stochastic models (as MLP for example).

2.1.

Satellite-derived hourly radiation

Many methods allow a precise knowledge of spatial distribution and temporal behavior of SSI over a territory. Among these methods, remote sensing allows processing efficiency satellite images to estimate the cloud cover related to the studied area and so, indirectly the SSI. The interest of SSI estimated from satellite data for a great area has been demonstrated in (Muselli et al., 1998a;[START_REF] Perez | Chapter 2 -Semi-Empirical Satellite Models[END_REF]. In this study, SSI was extracted from the Helioclim-3 database [START_REF] Vernay | Characterizing measurements campaigns for an innovative calibration approach of the global horizontal irradiation estimated by HelioClim-3[END_REF] centered on Corsica. This database provides hourly SSI maps with a nadir spatial resolution of 2.5 km obtained from Meteosat Second Generation (MSG) images treated by the Heliosat-2 model [START_REF] Guemene Dountio | On the reliability of HELIOSAT method: A comparison with experimental data[END_REF]. This model uses calibrated satellite images that are converted into irradiance.

They result from radiation-atmosphere interactions: a fluctuation of the signal measured by the sensors is interpreted as a variation of the cloud cover. Heliosat-2 computes the SSI from the clear sky radiation multiplied by a cloud index which expresses the clearness of the atmosphere for the pixel and the time t (equation 1).

(1)

In this equation, is the reflectance (or apparent albedo) measured by the sensor, is the apparent albedo of the brightest clouds and is the apparent albedo of the ground under clear sky. [START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF] have explicated the relation between the parameter and the clear sky index noted as described in the equation 2.

(2) [START_REF] Rigollier | On the clear sky model of the ESRA -European Solar Radiation Atlaswith respect to the heliosat method[END_REF]. In this model, the global horizontal irradiance under clear sky is split into direct and diffuse components.

Their computations are based on the Linke turbidity factor [START_REF] Remund | Worldwide Linke turbidity information[END_REF] that summarizes the turbidity of the atmosphere, and hence the attenuation of the direct beam and the generation of the diffuse fraction [START_REF] Kasten | The linke turbidity factor based on improved values of the integral Rayleigh optical thickness[END_REF].

Studied zone and data

Corsica is a Mediterranean island located in the golf of Genoa. It stretches out on 183 km long between latitudes 41° 20' 02" N and 43° 01' 31" N and 83.5 km wide between longitudes 8° 32' 30" E and 9° 33' 38" E extending to an area of 8722 km 2 . The island's seaside globally benefits from a Mediterranean climate characterized by hot dry summers and mild wet winters. Corsica presents one of the greatest solar potential of metropolitan France. The island's mountainous inner has a climate of height elevations or even alpine climate. The summers are hot and sunny but the weather is more instable on winters and the precipitations are abundant. Thus, the geographical distribution of solar potential and the dynamics of meteorology are particularly heterogeneous in such a territory. The public meteorological network is composed of only six stations. Only three of them provide hourly radiations measurements while the others provide 10 days integrated data (Haurant et al., 2012). In this context, predictions from 1-D TS modeling can be realized for three punctual locations in Corsica, not enough to allow a global management of electrical network and photovoltaic plants. Thus 2D prediction map from satellite SSI maps can be a good alternative related to the lack of measured data.

In the present study, 1158 points are considered, defined by the HelioClim-3 meshgrid centered on Corsica (not necessarily coincident with the native MSG pixels). Figure 1

Prediction models

There are a lot of different approaches to model time series and comprehensive reviews are available [START_REF] Hamilton | Disaggregation of satellite derived irradiance maps: Evaluation of the process and application to Corsica[END_REF][START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF]Şahin et al., 2012.;Voyant et al., 2013). The common base of all these models is that an element of TS ( ) can be defined by a linear or non-linear model called f n (see Equation 3where t = n,n-1,…,p+1,p with n, the number of observations and p the 184 km number of parameters with n p,and the error term of the prediction à time t+1 and for the pixel ) [START_REF] Crone | Stepwise Selection of Artificial Neural Networks Models for Time Series Prediction[END_REF].

(3)

In the extended territory case, the previous equation is applicable for each pixel of the image as described in the Figure 2. The following part is dedicated to particular families of functions , especially related to the persistence and its scaled form, artificial neural networks and to the cloudless approximation and modeling. Note that the ARMA approach is not tested here, because in a previous study (Voyant et al., 2013b) we show that MLP and ARMA give the same results for the global radiation prediction in Corsica Island in the univariate case. i th row j th line x y SSI at t-1 SSI at t SSI at t-2 . .

Persistence and scaled persistence

The first of studied prediction methods is the persistence method; the simplest way of producing a forecast. The persistence method assumes that the atmospheric conditions are stationary as shown in the equation 4 (Voyant et al., 2013a).

(4)

To take into account the fact that the apparent position of the sun is not identical between t and t+1, it is possible to correct the previous form with a clear sky ratio term (see equation 5) [START_REF] Sfetsos | Univariate and multivariate forecasting of hourly solar radiation with artificial intelligence techniques[END_REF] and to generate the scaled persistence.

(5)

Clear sky approximation

Another way to estimate the global radiation the next hour is to consider that we are dealing exclusively with sunny days (without cloud cover). In this case, the prediction is done only with determinist model as clear sky model presented previously (see equation 6) [START_REF] Rigollier | On the clear sky model of the ESRA -European Solar Radiation Atlaswith respect to the heliosat method[END_REF].

(6)

Artificial neural networks

To estimate the global radiation with stochastic models, a stationary hypothesis is often necessary. This result, initially shown for autoregressive and moving average processes [START_REF] Hamilton | Disaggregation of satellite derived irradiance maps: Evaluation of the process and application to Corsica[END_REF]) can be also applicable for the study and prediction from artificial neural networks (ANN) [START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF]. In fact, before to use ANN, time series must be made stationary. In previous work we have demonstrated that the clear sky index is a stationary computed series that can be directly operated by neural networks [START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF].

Although a large range of different architecture of ANN is available, MultiLayer Perceptron (MLP) remains the most popular [START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF]. In particular, feed-forward MLP networks with two layers (one hidden layer and one output layer) are often used for modeling and forecasting time series. Several studies have validated this approach based on ANN for the non-linear modeling of time series. To forecast the next value of a time series, a fixed number p of past values are set as inputs of the MLP, the output is the prediction of the future value [START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF]. Considering the initial time series equation (equation 3), we can transform this formula in the non-linear case of one hidden layer MLP with b related to the biases, f and g to the activation functions of the output and hidden layer, and to the weights (see equation 7). The number of hidden nodes (H) and input nodes (In) allow detailing this transformation. Note that as MLP is a stationary estimator, it is the clear sky index which is used during the modeling which defines a time series made stationary. To build a predictor of the global radiation, it is possible to construct and train a MLP with the clear sky estimation at time t+1 and the clear sky index of previous lags.

With

(7)

In In order to optimize the different MLP, it is essential to construct a tool able to determine the number and the nature of the MLP inputs. We propose to use the automutual information tool (MI(X t , L  X t ) for > 1 in equation 8) which is a quantity measuring the automutual dependence of a variables X t (L and  the lag operator and associated order). In fact, this formalism replaces and generalizes the autocorrelation which allows to measures only the linear relationship. Mutual information is more general

(x i ,y j ) (x i ,y j ) (x i ,y j ) (x i ,y j )
and measures the reduction of uncertainty in L  X t after observing X t [START_REF] Parviz | Mutual Information Based Input Variable Selection Algorithm and Wavelet Neural Network for Time Series Prediction[END_REF]. So MI [START_REF] Jiang | Mutual information algorithms[END_REF] can measure non-monotonic and other more complicated relationships. It can be expressed as a combination of marginal and conditional entropies (respectively and ) as described in the equation 8 [START_REF] Huang | Effective feature selection scheme using mutual information[END_REF].

) X L )-H(X )=H(X X ,L MI(X t t t t t           
Entropy corresponds to a measure of unpredictability or information content and can be written by the expression detailed in the equation 9 (entropy of a discrete random variable X t with possible values x={x 1 ,..x n }).

        
One may also define the conditional entropy of two events X t and L  X t (Equation 10). This quantity should be understood as the amount of randomness of the random variable X t given that you know the value of Y t =L  X t (possible values y={y 1 ,..y n }).

(10)

The definition of the joint probability distribution function (p(x,y)) and marginal probabilities (p(x) and p(y)), allows to define a new form of the mutual information as described in the equation 11.

(11)

We consider that the maximum lag to take as input of the MLP corresponds to the first minimum of the automutual information graph ( ). In example presented in the Figure 4 The first minimum of auto-mutual information concerning all the points of the meshing ( 1158different operations) has a median value equal to 7 time lags (min=5, max=10, mean value=7.63 and standard deviation=1.08). Thereby, for the overall territory, the number of inputs (In) is takes equal to 7 (i.e. 7 hours of passed measurements). Related to previous studies (Voyant et al., 2013a[START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF], we consider that the optimal configuration of the MLP corresponds to one hidden layer provided with 7 hidden neurons (size of inputs and hidden nodes should be equivalent).

Results

Before exposing prediction results, it is important to define error metrics. During this study, two parameters are dedicated to the standard TS analysis (see Equations 12 and 13) and one has been constructed to take into account the 2-D aspect of predictions and measurements (Equation 14 and 15).

The nRMSE metrics is expressed with the expected value (E) which is in this case, the temporal mean computed on the overall test sampling.

(12)

(13)

In the present study, the clear sky index is used in all the predictive methodologies and is related to some geographical (altitude, longitude or latitude) and modeling parameters (Linke turbidity factor or relative air mass). Because of the sampling process of the HC-3 SSI, it is possible that for some pixels, these parameters could not be adapted (especially altitude, if it is in a mountainous region). In this context, we need a tool giving information about the fact that neighbor pixels could have a better forecasted value. In physics and particularly in medical physics, the gamma index is a criterion for comparing data from two matrixes. This interpretation is based on the assumption that sometime the discrepancy between measures and calculi from the model may not be due to an error of the map calculation algorithm or delivery hardware, but simply due to experimental error (different offsets or spatial resolutions). In the high gradient region, a small offset could induce a significant error between real pixel and forecasted one while the modeling is good. More abstractly, the gamma index takes two (p), the gamma index ( ) is defined from the distance (noted , Equation 14) between points in reference (or measurement;

) and evaluation distributions (or prediction ; ), note that polar coordinates are used but a Cartesian approach is also possible: Each pixel of the image is associated with a gamma index value. The higher is this index, the less the prediction is good so the determination of γ throughout the predicted irradiance distribution provides a presentation that quantitatively indicates the calculation accuracy. Moreover, regions where  > 1 correspond to locations where the calculation does not meet the acceptance criteria. Note that, it possible to construct the same type of test with the hypothesis that Tol r +∞, in this case the physical distance is not taken into account, only the pixels intensity are considered ( ).

With this index, it is possible to construct a statistical test based on the null hypothesis H0: the prediction related to the pixel is relevant. If > 1, H0 is rejected else H0 is accepted. This test is represented in the Figure 5 for each pixel in a space composed of pixel intensity and spatial coordinates. The acceptance criteria form an ellipsoid surface, the major axis scales of which are determined by individual acceptance criteria and the center of which is located at the estimated point in question. When the calculated intensity distribution surface passes through the ellipsoid, the calculation passes the acceptance test for the measurement point. The measure of acceptability is the multidimensional distance between the measurement and calculation points in both the pixel intensity and the physical distance, scaled as a fraction of the acceptance criteria. For an overall image, it is possible to synthetize the test computing the gamma passing rate (%GP; for each map it is the percentage of pixels passing the test), higher is this percentage and better is the prediction.

Performances of HC-3 SSI estimations

The data used during this study are related to global solar estimation from the Helioclim-3 database. The performances of the satellite derivate SSI estimations have been evaluated in the case of Corsica. In that view, they have been compared to hourly radiometric measurements provided by three meteorological stations in Corsica (Fig 1 ,table 1) that cover the period from April 1 st , 2004 to December 17 th , 2006.

The table 1 summarizes seasonal performances of SSI estimations for the three stations. It is shown these performances when irradiations under 0.1 Wh/m² and 10 Wh/m² are eliminated. In the first case all the night values are removed (~6:00-22:00 in summer; period corresponding to daytime)

and concerning the second case, one or two additional hours are removed. Imposing such a threshold allows limiting uncertainties due to the air optical air mass artifacts at sunrise and sunset and shading effects. Firstly, we can see that the 10 Wh/m² threshold allows decreasing the errors of 0.5 to 1.5 points demonstrating that the sunrise and sunset artefacts induce important errors, especially in Corte (the only mountainous site) but also in Bastia. Don't forget that the nRMSE is also linked with errors in the reference data: the location of radiometric stations is not necessarily relative to the center of each pixel. We can also observe that the cloudy periods (winter -autumn) present the most important error (nRMSE between 20 and 30 %) showing that SSI under cloudy sky are particularly difficult to be calculated by the model as we could have expected. However, during cloudy months the increasing of observed nRMSE can be also or partly explained by a decreasing of the reference irradiation (see equation 13). In fact, the HelioClim3 database (due to the pixel size) is not very sensitive to small clouds formations while they can have a big impact on the ground measures. Finally we can note that the estimations are correlated with measurements since yearly and seasonal cross correlation coefficients are over 0.96. Yearly nRMSE are between 18.0 and 27.2 %, comparable but lightly over the nRMSE of Meteosat-8 estimations for the same stations and the same periods (Haurant et al., 2012). However, these data have been selected for this study because HelioClim-3 is one of the databases that offer SSI maps (or data series) on a long period but others models exist (CM-SAF, Land-SAF, etc.) as described in [START_REF] Ineichen | Satellite Application Facilities irradiance products: hourly time step comparison and validation over Europe[END_REF]. Also, the HelioClim-3 extracted maps have higher resolution than Meteosat-8 images. HelioClim-3 allows also the possibility of near real-time SSI acquisitions: the day after in standard mode and on specific requests, 15 minutes MSG acquisitions plus some minutes on a selected region.

Periods Threshold Ajaccio

Irradiation prediction validation

Now the sensibility study is done concerning the validation of the HC-3 data, in this part we will expose the performance of the four methods of prediction described in the section 2. During all the experimentations only hours between 8:00 and 18:00 are used (UTC), outside this interval the choice to predict the global radiation is not really justified. The table 2 describes the comparison between the four predictors. The nRMSE, the gamma index and the gamma passing rate are given for the four seasons. The gamma index is related to a = 1 pixel (2.5 km) and = 10%. We choose intuitively these values because this parameter is never used in global radiation forecasting. If gamma is less than 1 it corresponds to a prediction error lower than 10% with a precision better than 2.5 km. The two better predictors are the scaled persistence and the MLP, their nRMSE is less than 17%, the average gamma index is lower than 1 (except 1.21 for MLP in spring) and the gamma passing rate (%GP) is upper than 65% for the scaled persistence and upper than 45% for MLP.

quality of the prediction is to generate the algebraic error ( between prediction (for typical days) and measurement as shown in figure 7. If in winter the error is significant (a lot of red pixels indicating that absolute error is upper than 100Wh/m²), in summer it is close to 0Wh/m² (a lots of light blue pixels). In this example, there is not significant relationship between the altitude and the error of prediction (mountain ridge in the center of the island longitudinally). To appreciate visually the quality of prediction, it is essential to draw profiles of HC-3 SSI map and predictions as in the Figures 91011. In these plots we see that during the summer period the prediction is effective unlike the other seasons.

Conclusions

With this article, we add a new original approach for global radiation prediction. A methodology usually applied for 1-D time series modeling has been modified and used to generate a 2-D predictions map for an overall territory (Corsica Island in our case). Classical time series analyze tools were used, however a lot of new methods are developed:

-usually, the TS formalism is applied for a specific location and not covering an overall spatial territory (resolution of 2.5km); -the mutual information as MLP optimization parameter is rarely used; -the scaled persistence is very rarely used in time series prediction; -the gamma index is never used in global radiation forecasting.

The methodology applied admits two sources of uncertainty, the first one is linked to the Helioclim-3 database and the global radiation computing from satellite acquisitions, the second one is linked to the prediction method used. The value of nRMSE including the two error components is close to 37% (respectively 20.1% and 16.5%). Concerning the gamma passing rate, our study shows that for winter and spring the scaled persistence gives the best results (respectively 67.7% and 86%), for autumn, it is the simple persistence (95.3%) and for summer it is the MLP (99.8%). We think that this new methodology based on the 2-D time series modeling with scaled persistence (or MLP) has to be compared against the complex model NWP (with different pixels sizes and origins offsets). The next step will be to compare in Corsica, using the gamma index and nRMSE, the prediction of NWP reference methodologies (NOAA and AROME) and the time series approaches (MLP and scaled persistence). Moreover, if in the 1-D case the MLP gives very good result, in the 2-D case the interest of this methodology is not shown in this paper, new adaptations are necessary (specialized MLP by geographical clustering, use of adjacent pixels values, round off the optimization of each MLP, etc.).

Finally, with MLP there is the operational possibility to integrate to the forecast schema, existing realtime measurements of irradiation and/or PV production at different locations in the region of interest. This approach which includes a possible calibration with historical existing databases and autocorrection in the real-time basis is certainly a good way to dramatically decrease the nRMSE of the forecast (Kalman filter, recursive least squares filter, etc.). Moreover, in this paper, inputs of models are previously measured values, the multivariate case (with exogenous data as inputs: wind speed, temperature, humidity, etc.) could improve the result and it will be a good perspective of future work.

.

List of captions

Table 1. Comparison between helioclim-3 satellite SSI and ground measurements for the 3 studied stations in Corsica (nRMSE in %)

Table 2. Comparison of the 4 predictors, the best values are in bold

  x={ x 1 … x n } Y t / y={ y 1 … y n } Time series X t and Y t, and possible entropies and the mutual information (bit) Albedo of the ground under clear sky (unitless) p(x), p(y) and P(x,y) Marginal and joint probabilities distribution function of X t , and Y t Clear sky index (unitless) Normalized root mean square error (

  shows this meshgrid and three of the principal cities of the island where there are weather stations. In fact, each SSI pixel is a generator of one global radiation time series. Eight years of HC-3 hourly global radiations are available, from 2005 to 2012 for these three places. The two last years (2011-2012) are used for testing the different predictors and the six first years (2005-2010) are dedicated to the learning phase concerning the MLP approach.

Figure 1 .

 1 Figure 1. Points of the HelioClim3 meshgrid and location of the 3 meteorological stations providing hourly solar radiation measurements, circle locates Ajaccio (41°55'N and 8°44'E, elev. 0-787 m), square locates Corte (42°18'N and 9°09'E, elev. 300-2626 m) and triangle locates Bastia (42°42'N and 9°27'E, elev. 0-963 m)

Figure 2 .

 2 Figure 2. TS definition from successive satellite acquisitions. In dark grey the pixel at

  the presented study, the MLP has been computed with the Matlab© software and its Neural Network toolbox. The characteristics chosen and related to previous works(Voyant et al., 2013b): one hidden layer, the activation functions are the continuously and differentiable hyperbolic tangent (hidden) and the linear function (output), the Levenberg-Marquardt (approximation to the Newton's method) learning algorithm with a max fail parameter before stopping training equal to 3(Voyant et al., 2013b). The detail of the MLP architecture used is shown in the Figure3. In fact a similar MLP is defined for all the pixels of the radiation map: 1158 different MLP without any hypothesis on the spatial coherency.

Figure 3 .

 3 Figure 3. MLP architecture for the pixel (x i ,y i )

  (the MI versus the time lag for a pixel near Ajaccio), there is a first minimum corresponding to the 10 th time lag, and thereby inducing a MLP with 10 inputs ( .

Figure 4 .

 4 Figure 4. Example of automutual information of the clear sky index in the Ajaccio site

  2D arrays and compares each element of array_1 with spatially-nearby elements of array_2. The gamma test was first introduced by[START_REF] Low | A technique for the quantitative evaluation of dose distributions[END_REF] as a single metric that combined features of both pixel intensity difference (PID; as describe in equations 12 or 13) and distance to agreement (DTA; e.g. distance between iso-intensity curves), allowing a robustly performing in the regions prone to failure. It is the minimum Euclidean distance in a normalized intensity pixel-distance space. In this way, both PID and DTA are taken into account for every point in the evaluated distribution[START_REF] Thomas | A Comparison of Four Indices for Combining Distance and Dose Differences[END_REF] generating a dimensionless metric. For a measured image (m) and predicted image

  intensity difference at point in the predicted image and point in the measured image ( , so the gamma index ( is defined for two tolerance terms ( by the equation 15 (Tol r is in fact a fixed DTA and Tol I , a fixed PID). The minimum radial distance between the measurement point and the calculation points (expressed as a surface in the intensity-distance space of the Figure 5) induced the γ index metrics.

Figure 5 .

 5 Figure 5. Schematic representation of the gamma index test. Inside the ellipse, the test is accepted (red point) and outside, it is rejected (black point).

Figure 7 .

 7 Figure 7. Forecasted SSI with scaled persistence, HC-3 SSI and algebraic error (Color map in Wh/m²)

ForecastsFigure 8 .

 8 Figure 8. Distribution of the annual nRMSE (a. spatial repartition and b. histogram)

Figure 9 .

 9 Figure 9. Profiles of HC-3 SSI data and scaled persistence estimation of the global radiation (respectively line and marks) in a pixel localized near Ajaccio

Figure 1 .

 1 Figure 1. Points of the HelioClim3 meshgrid and location of the 3 meteorological stations providing hourly solar radiation measurements, circle locates Ajaccio (41°55'N and 8°44'E, elev. 0-787 m), square locates Corte (42°18'N and 9°09'E, elev. 300-2626 m) and triangle locates Bastia (42°42'N and 9°27'E, elev. 0-963 m) Figure 2. TS definition from successive satellite acquisitions. In dark grey the pixel at time t defining the global radiation notedThe intensity of successive pixels defines each TS ( ). Figure3. MLP architecture for the pixel (x i ,y i ) Figure4. Example of automutual information of the clear sky index in the Ajaccio site Figure5. Schematic representation of the gamma index test. Inside the ellipse, the test is accepted (red point) and outside, it is rejected (black point).Figure6. Gamma test for 4 typical days with scaled persistence. Black pixel if the test is rejected and white pixel if it is passed Figure 7. Forecasted SSI with scaled persistence, HC-3 SSI and algebraic error (Color map in Wh/m²) Figure 8. Distribution of the annual nRMSE (a. spatial repartition and b. histogram) Figure 9. Profiles of HC-3 SSI data and scaled persistence estimation of the global radiation (respectively line and marks) in a pixel localized near Ajaccio Figure 10. Profiles of HC-3 SSI data and scaled persistence estimation of the global radiation (respectively line and marks) in pixel localized near Corte Figure 11. Profiles of HC-3 SSI data and scaled persistence estimation of the global radiation (respectively line and maks) in pixel localized near Bastia

Figure 6 .

 6 Figure 2. TS definition from successive satellite acquisitions. In dark grey the pixel at time t defining the global radiation notedThe intensity of successive pixels defines each TS ( ). Figure3. MLP architecture for the pixel (x i ,y i ) Figure4. Example of automutual information of the clear sky index in the Ajaccio site Figure5. Schematic representation of the gamma index test. Inside the ellipse, the test is accepted (red point) and outside, it is rejected (black point).Figure6. Gamma test for 4 typical days with scaled persistence. Black pixel if the test is rejected and white pixel if it is passed Figure 7. Forecasted SSI with scaled persistence, HC-3 SSI and algebraic error (Color map in Wh/m²) Figure 8. Distribution of the annual nRMSE (a. spatial repartition and b. histogram) Figure 9. Profiles of HC-3 SSI data and scaled persistence estimation of the global radiation (respectively line and marks) in a pixel localized near Ajaccio Figure 10. Profiles of HC-3 SSI data and scaled persistence estimation of the global radiation (respectively line and marks) in pixel localized near Corte Figure 11. Profiles of HC-3 SSI data and scaled persistence estimation of the global radiation (respectively line and maks) in pixel localized near Bastia

Table 2 .

 2 Comparison of the 4 predictors, the best values are in bold

	*gamma passing rate
	Criteria	Scaled persistence	MLP	Persistence	Clear sky model
	nRMSE (%)	16.61	16.54	34.11	45.04
	 wi  sp  su  au	0.84 0.36 0.04 0.27	1.21 0.94 0.27 0.58	1.21 1.3 0.28 0.51	4.9 1.01 0.12 0.23
	%GP wi *	67.7	45.6	47.6	0.01
	%GP sp	86.0	58.4	12.4	78.4
	%GP su	99.5	99.8	99.7	98.6
	%GP au	93.8	94.5	95.3	94.5

Concerning the clear sky model, the gamma passing rate is very good in summer (more than 98%; related to few cloud occurrence) and very bad in winter (less than 1%; related to a lot of cloud occurrences). Note that the variability of the global radiation from satellite acquisitions is less important than in the case of ground-based measurements. If in the first case the nRMSE of the prediction is close to 20% and in the second it is less than 17%. In fact, the "small clouds" are not well considered in the Helioclim-3 data (because of the pixels size), there is sometime overestimation and then underestimation regarding ground measurements. As scaled persistence and MLP predictions are equivalent, the simplicity of the first one, led us to think it is the most interesting. In the Figure 6 the gamma test is presented for 4 typical days related to each season with the scaled persistence. Black pixels if the test is rejected and white if it is passed. During summer and autumn, there are a lot of white pixels, almost all the predictions are in the ellipse 2.5km/10%. For the other seasons, results are less uniform: in spring the number of black and white pixels is equivalent but in winter, there are more black than white pixels. Another way to see the

Acknowledgement

Thanks to Armines and to the SoDa community for making the HelioClim data available (http://www.helioclim.org/index.html).


Maurice Robert Kibler
email: kibler@ipnl.in2p3.fr

On Two Ways to Look for Mutually Unbiased Bases

Keywords: finite-dimensional quantum mechanics, quantum information, MUBs, SIC POVMs, equiangular lines, equiangular vectors

Two equivalent ways of looking for mutually unbiased bases are discussed in this note. The passage from the search for d+1 mutually unbiased bases in C d to the search for d(d+1) vectors in C d 2 satisfying constraint relations is clarified. Symmetric informationally complete positive-operatorvalued measures are briefly discussed in a similar vein.

Introduction

The concept of mutually unbiased bases (MUBs) plays an important role in finite-dimensional quantum mechanics and quantum information (for more details, see [START_REF] Vourdas | Quantum systems with finite Hilbert space[END_REF][START_REF] Tolar | Feynman's path integral and mutually unbiased bases[END_REF][START_REF] Kibler | An angular momentum approach to quadratic Fourier transform, Hadamard matrices, Gauss sums, mutually unbiased bases, unitary group and Pauli group[END_REF][START_REF] Durt | On mutually unbiased bases[END_REF] and references therein). Let us recall that two orthonormal bases {|aα : α = 0, 1, . . . , d -1} and {|bβ : β = 0, 1, . . . , d-1} in the d-dimensional Hilbert space C d (endowed with an inner product denoted as | ) are said to be unbiased if the modulus of the inner product aα|bβ of any vector |bβ with any vector |aα is equal to 1/ √ d. It is known that the maximum number of MUBs in C d is d + 1 and that this number is reached when d is a power of a prime integer. In the case where d is not a prime integer, it is not known if one can construct d + 1 MUBs (see [START_REF] Durt | On mutually unbiased bases[END_REF] for a review).

In a recent paper [START_REF] Kibler | Equiangular vectors approach to mutually unbiased bases[END_REF], it was discussed how the search for d + 1 mutually unbiased bases in C d can be approached via the search for d(d + 1) vectors in C d 2 satisfying constraint relations. It is the main aim of this note to make the results in [START_REF] Kibler | Equiangular vectors approach to mutually unbiased bases[END_REF] more precise and to show that the two approaches (looking for d + 1 MUBs in C d or for d(d + 1) vectors in C d 2 ) are entirely equivalent. The central results are presented in Sections 2 and 3. In Section 4, parallel developments for the search of a symmetric informationally complete positive-operator-valued measure (SIC POVM) are considered in the framework of similar approaches. Some concluding remarks are given in the last section.

The two approaches

It was shown in [START_REF] Kibler | Equiangular vectors approach to mutually unbiased bases[END_REF] how the problem of finding d + 1 MUBs in C d , i.e., d + 1 bases

B a = {|aα : α = 0, 1, . . . , d -1} (1)
satisfying

| aα|bβ | = δ α,β δ a,b + 1 √ d (1 -δ a,b ) (2)
can be transformed in the problem of finding d(d + 1) vectors w(aα) in C d 2 , of components w pq (aα), satisfying

w pq (aα) = w qp (aα), p, q ∈ Z/dZ (3) d-1 p=0 w pp (aα) = 1 (4)
and

d-1 p=0 d-1 q=0 w pq (aα)w pq (bβ) = δ α,β δ a,b + 1 d (1 -δ a,b ) (5)
with a, b = 0, 1, . . . , d and α, β = 0, 1, . . . , d -1 in (1)-( 5). (In this paper, the bar denotes complex conjugation.) This result was described by Proposition 1 in [START_REF] Kibler | Equiangular vectors approach to mutually unbiased bases[END_REF]. In fact, the equivalence of the two approaches (in C d and C d 2 ) requires that each component w pq (aα) be factorized as

w pq (aα) = ω p (aα)ω q (aα) (6)
for a = 0, 1, . . . , d and α = 0, 1, . . . , d -1, a condition satisfied by the example given in [START_REF] Kibler | Equiangular vectors approach to mutually unbiased bases[END_REF]. The factorization of w pq (aα) follows from the fact that the operator Π aα defined in [START_REF] Kibler | Equiangular vectors approach to mutually unbiased bases[END_REF] is a projection operator. The introduction of ( 6) in ( 3), ( 4) and ( 5) leads to some simplifications. First, (6) implies the hermiticity condition (3). Second, by introducing ( 6) into ( 4) and ( 5), we obtain

d-1 p=0 |ω p (aα)| 2 = 1 (7)
and

d-1 p=0 ω p (aα)ω p (bβ) 2 = δ α,β δ a,b + 1 d (1 -δ a,b ) (8)
respectively. It is clear that [START_REF] Appleby | Group theoretic, Lie algebraic and Jordan algebraic formulations of the SIC existence problem[END_REF] follows from ( 8) with a = b and α = β. Therefore, (3) and ( 7) are redundant in view of ( 5) and ( 6). As a consequence, Maurice R. Kibler

Acta Polytechnica

PREPRINT Proposition 1 in [START_REF] Kibler | Equiangular vectors approach to mutually unbiased bases[END_REF] can be precised and reformulated in the following way.

Proposition 1. For d ≥ 2, finding d + 1 MUBs in C d (if they exist) is equivalent to finding d(d + 1) vectors w(aα) in C d 2 , of components w pq (aα) such that d-1 p=0 d-1 q=0 w pq (aα)w pq (bβ) = δ α,β δ a,b + 1 d (1 -δ a,b ) (9)
and w pq (aα) = ω p (aα)ω q (aα), p, q ∈ Z/dZ [START_REF] Brierley | Constructing mutually unbiased bases in dimension six[END_REF] where a, b = 0, 1, . . . , d and α, β = 0, 1, . .

(M aα ) pq = ω p (aα)ω q (aα), p, q ∈ Z/dZ (11)
and satisfying the trace relations

Tr (M aα M bβ ) = δ α,β δ a,b + 1 d (1 -δ a,b ) (12)
where a, b = 0, 1, . . . , d and α, β = 0, 1, . . . , d -1.

Equivalence

Suppose that we have a complete set {B a : a = 0, 1, . . . , d} of d + 1 MUBs in C d , i.e., d(d + 1) vectors |aα satisfying (2), then we can find d(d + 1) vectors w(aα) in C d 2 , of components w pq (aα), satisfying [START_REF] Butterley | Numerical evidence for the maximum number of mutually unbiased bases in dimension six[END_REF] and [START_REF] Brierley | Constructing mutually unbiased bases in dimension six[END_REF]. This can be achieved by introducing the projection operators

Π aα = |aα aα| ( 13
)
where a = 0, 1, . . . , d and α = 0, 1, . . . , d -1. In fact, it is sufficient to develop Π aα in terms of the E pq generators of the GL(d, C) complex Lie group; the coefficients of the development are nothing but the w pq (aα) complex numbers satisfying ( 9) and [START_REF] Brierley | Constructing mutually unbiased bases in dimension six[END_REF], see [START_REF] Kibler | Equiangular vectors approach to mutually unbiased bases[END_REF] for more precisions.

Reciprocally, should we find d(d + 1) vectors w(aα) in C d 2 , of components w pq (aα), satisfying ( 9) and ( 10), then we could construct d(d + 1) vectors |aα satisfying (2). This can be done by means of a diagonalization procedure of the matrices

M aα = d-1 p=0 d-1 q=0 w pq (aα)E pq ( 14
)
where a = 0, 1, . . . , d and α = 0, 1, . . . , d -1. An alternative and more simple way to obtain the |aα vectors from the w(aα) vectors is as follows. Equation (8) leads to

d-1 p=0 ω p (aα)ω p (bβ) = δ α,β δ a,b + 1 √ d (1 -δ a,b ) (15)
to be compared with [START_REF] Tolar | Feynman's path integral and mutually unbiased bases[END_REF]. Then, the |aα vectors can be constructed once the w(aα) vectors are known. The solution, in matrix form, is

|aα =      ω 0 (aα) ω 1 (aα) . . . ω d-1 (aα)      (16) a = 0, 1, . . . , d α = 0, 1, . . . , d -1 (17)
Therefore, we can construct a complete set {B a : a = 0, 1, . . . , d} of d+1 MUBs from the knowledge of d(d+ 1) vectors w(aα). Note that, for fixed a and α, the |aα vector is an eigenvector of the M aα matrix with the eigenvalue 1. This establishes a link with the above-mentioned diagonalization procedure.

A parallel problem

The present work takes its origin in [START_REF] Albouy | A unified approach to SIC-POVMs and MUBs[END_REF] where some similar developments were achieved for the search of a SIC POVM. Symmetric informationally complete positive-operator-valued measures play an important role in quantum information. Their existence in arbitrary dimension is still the object of numerous studies (see for instance [START_REF] Appleby | Group theoretic, Lie algebraic and Jordan algebraic formulations of the SIC existence problem[END_REF]). A SIC POVM in dimension d can be defined as a set of d 2 nonnegative operators P x = |Φ x Φ x | acting on C d and satisfying 1 d

d 2 x=1 P x = I (18)
and

Tr (P x P y ) = dδ x,y + 1 d + 1 ( 19
)
where I is the identity operator. The search for such a SIC POVM amounts to find

d 2 vectors |Φ x in C d satisfying 1 d d 2 x=1 |Φ x Φ x | = I (20)
and

| Φ x |Φ y | = dδ x,y + 1 d + 1 (21)
with x, y = 1, 2, . . . , d 2 . The P x operator can be developed as

P x = d-1 p=0 d-1 q=0 v pq (x)E pq (22)
so that the determination of d 2 operators P x (or

d 2 vectors |Φ x ) is equivalent to the determination of d 2 vectors v(x), of components v pq (x), in C d 2 .
In the spirit of the preceding sections, we have the following result.

PREPRINT vol. no. / On Two Ways to Look for Mutually Unbiased Bases

Proposition 3. For d ≥ 2, finding a SIC POVM in C d (if it exists) is equivalent to finding d 2 vectors v(x) in C d 2 , of components v pq (x) such that 1 d d 2 x=1 v pq (x) = δ p,q , p, q ∈ Z/dZ (23) d-1 p=0 d-1 q=0 v pq (x)v pq (y) = dδ x,y + 1 d + 1 ( 24
)
and v pq (x) = ν p (x)ν q (x), p, q ∈ Z/dZ (25)

where x, y = 1, 2, . . . , d 2 .

Concluding remarks

The equivalence discussed in this work of the two ways of looking at MUBs amounts in some sense to the equivalence between the search for equiangular lines in C d and for equiangular vectors in C d 2 (cf. [START_REF] Godsil | Equiangular lines, mutually unbiased bases, and spin models[END_REF]). Equiangular lines in C d correspond to Observe that the modulus disappears and the 1/ √ d factor is replaced by 1/d when passing from (26) to (27). It was questioned in [START_REF] Kibler | Equiangular vectors approach to mutually unbiased bases[END_REF] if the equiangular vectors approach can shed light on the still unsolved question to know if one can find d+1 MUBs when d is not a (strictly positive) power of a prime integer. In the case where d is not a power of a prime, the impossibility of finding d(d + 1) vectors w(aα) or d(d + 1) matrices M aα satisfying the conditions in Propositions 1 and 2 would mean that d + 1 MUBs do not exist in C d . However, it is hard to know if one approach is better than the other. It is the hope of the author that the equiangular vectors approach be tested in the d = 6 case for which one knows only three MUBs instead of d + 1 = 7 in spite of numerous numerical studies (see [START_REF] Butterley | Numerical evidence for the maximum number of mutually unbiased bases in dimension six[END_REF][START_REF] Brierley | Constructing mutually unbiased bases in dimension six[END_REF][START_REF] Mcnulty | On the impossibility to extend triples of mutually unbiased product bases in dimension six[END_REF] and references therein for an extensive list of related works).

| aα|bβ | = 1 √ d for a = b ( 26
Similar remarks apply to SIC POVMs. The existence problem of SIC POVMs in arbitrary dimension is still unsolved although SIC POVMs have been constructed in every dimension d ≤ 67 (see [START_REF] Appleby | Group theoretic, Lie algebraic and Jordan algebraic formulations of the SIC existence problem[END_REF] The parallel between MUBs and SIC POVM characterized by the couples of equations ( 26)-( 29), ( 27)-( 30) and ( 28)-(31) should be noted. These matters shall be the subject of a future work.

w

  ) while equiangular vectors in C d 2 correspond to w(aα) • w(bβ) = 1 d for a = b (27)where the w(aα) • w(bβ) inner product in C d 2 is depq (aα)w pq (bβ) (28)

  and references therein). For SIC POVMs, the equiangular lines in C d correspond to| Φ x |Φ y | = 1 √ d + 1 for x = y (29)and the equiangular vectors inC d 2 to v(x) • v(y) = 1 d + 1 for x = y (30)where the v(x) • v(y) inner product in C d 2 is defined asv(x) • v(y) =

  . , d -1. This result can be transcribed in matrix form. Therefore, we have the following proposition. Proposition 2. For d ≥ 2, finding d + 1 MUBs in C d (if they exist) is equivalent to finding d(d + 1) matrices M aα of dimension d, with elements

Acknowledgements

The material contained in the present note was planned to be presented at the eleventh edition of the workshop Analytic and Algebraic Methods in Physics (AAMP XI). Unfortunately, the author was unable to participate in AAMP XI. He is greatly indebted to Miloslav Znojil for suggesting him to submit this work to Acta Polytechnica.