Cyril Voyant
email: cyrilvoyant@gmail.com

Marc Muselli

Christophe Paoli

Marie-Laure Nivet

Numerical Weather Prediction (NWP) and hybrid ARMA/ANN model to predict global radiation

Keywords: Time Series forecasting, hybrid, Artificial Neural Networks, ARMA, Stationarity Periodic coefficients at time t

We propose in this paper an original technique to predict global radiation using a hybrid ARMA/ANN model and data issued from a numerical weather prediction model (ALADIN). We particularly look at the Multi-Layer Perceptron. After optimizing our architecture with ALADIN and endogenous data previously made stationary and using an innovative pre-input layer selection method, we combined it to an ARMA model from a rule based on the analysis of hourly data series. This model has been used to forecast the hourly global radiation for five places in Mediterranean area. Our technique outperforms classical models for all the places. The nRMSE for our hybrid model ANN/ARMA is 14.9% compared to 26.2% for the naïve persistence predictor. Note that in the stand alone ANN case the nRMSE is 18.4%. Finally, in order to discuss the reliability of the forecaster outputs, a complementary study concerning the confidence interval of each prediction is proposed

Introduction

Solar radiation is one of the principal energy sources, occupying a very important role in some engineering applications as production of electricity, heat and cold [START_REF] Muselli | Improved procedure for stand-alone photovoltaic systems sizing using METEOSAT satellite images[END_REF][START_REF] Cao | Study of forecasting solar irradiance using neural networks with preprocessing sample data by wavelet analysis[END_REF][START_REF] Kalogirou | Artificial neural networks in renewable energy systems applications: a review[END_REF][START_REF] Badescu | Modeling solar radiation at the earth's surface: recent advances[END_REF]. The process of converting sunlight to electricity without combustion creates power without pollution. It is certainly one of the most interesting themes in solar energy area. To use ideally this technology, it is necessary to understand and create efficient prediction models like done for example in [START_REF] Mueller | Rethinking satellite-based solar irradiance modelling: The SOLIS clear-sky module[END_REF][START_REF] Mellit | A simplified model for generating sequences of global solar radiation data for isolated sites: Using artificial neural network and a library of Markov transition matrices approach[END_REF] and Mubiru et al, 2008 [START_REF] Mueller | Rethinking satellite-based solar irradiance modelling: The SOLIS clear-sky module[END_REF][START_REF] Mubiru | Estimation of monthly average daily global solar irradiation using artificial neural networks[END_REF][START_REF] Mellit | A simplified model for generating sequences of global solar radiation data for isolated sites: Using artificial neural network and a library of Markov transition matrices approach[END_REF]. Insolation is defined as the solar radiation striking a surface at a certain time and place and is typically expressed in kilowatt hours per meter square (kWh/m²) [START_REF] Bird | Review, Evaluation, and Improvement of Direct Irradiance Models[END_REF][START_REF] Bird | A Simplified Clear Sky Model for Direct and Diffuse Insolation on Horizontal Surfaces[END_REF].

Many factors determine how much sunlight is available at a given location. We can mention the atmospheric conditions, the Earth's position in relation to the sun, and the site obstructions [START_REF] Badescu | Modeling solar radiation at the earth's surface: recent advances[END_REF][START_REF] Bird | A Simplified Clear Sky Model for Direct and Diffuse Insolation on Horizontal Surfaces[END_REF].

Atmospheric conditions that can also affect the amount of radiation received on the Earth's surface are the quantity of air molecules, water vapor, dust, ozone and carbon dioxide, the cloud cover, the air pollution, the dust storms, the volcanic eruptions, etc. [START_REF] Mueller | Rethinking satellite-based solar irradiance modelling: The SOLIS clear-sky module[END_REF][START_REF] Iqbal | An introduction to solar radiation[END_REF]. There is an interest to control the solar radiation prediction, as for example to identify the most optimal locations for developing solar power project or to maintain the grid stability in solar and conventional power management. Note that, in Europe, the White Paper on Energy (established in 1997) set a target (not yet achieved) of 12% of electricity production from renewable energies by 2010. The new European guidelines (in press) set a new target to 20% by 2020. The issues of the solar energy prediction are very important and mobilize a lot of research teams around the world and particularly in the Mediterranean area [START_REF] Kalogirou | Artificial neural networks in renewable energy systems applications: a review[END_REF][START_REF] Mellit | A simplified model for generating sequences of global solar radiation data for isolated sites: Using artificial neural network and a library of Markov transition matrices approach[END_REF][START_REF] Mellit | An adaptive model for predicting of global, direct and diffuse hourly solar irradiance[END_REF][START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF][START_REF] Christophe | Forecasting of preprocessed daily solar radiation time series using neural networks[END_REF]. In practice, the global radiation (or insolation) forecasting is the name given to the process used to predict the amount of solar energy available in the current and near terms. A lot of methods have been developed by experts around the world [START_REF] Mubiru | Estimation of monthly average daily global solar irradiation using artificial neural networks[END_REF][START_REF] Kaplanis | New methodologies to estimate the hourly global solar radiation; Comparisons with existing models[END_REF][START_REF] Elminir | Prediction of hourly and daily diffuse fraction using neural network, as compared to linear regression models[END_REF]. Often the Times Series (TS) mathematical formalism is necessary. It is described by sets of numbers that measures the status of some activity over time [START_REF] Gooijer | 25 years of time series forecasting[END_REF]. In primary studies [START_REF] Muselli | Improved procedure for stand-alone photovoltaic systems sizing using METEOSAT satellite images[END_REF][START_REF] Christophe | Forecasting of preprocessed daily solar radiation time series using neural networks[END_REF][START_REF] Voyant | Optimization of an artificial neural network dedicated to the multivariate forecasting of daily global radiation[END_REF] we have demonstrated that an optimized multi-layer perceptron (MLP) with endogenous input made stationary and exogenous inputs (meteorological data) can forecast the global solar radiation time series with acceptable errors. This prediction model has been compared to other prediction methods (ARMA, k-NN, Markov Chains, etc.) and we have concluded that ANN and ARMA were similar. Following these studies and in order to see if we can significantly improve our results, we decided to add weather forecast (instead of exogenous data previously used) as new inputs of our mode. We assumed that Numerical Weather Prediction (NWP) simulations tools compute data patterns essential for determining solar radiation. These weather forecast data present two advantages: firstly they are becoming more and more available through the internet and secondly the models provide spatially distributed data which are very relevant to the regional scale studies. Thanks to an agreement with Météo-France, which is the French meteorological organization (http://france.meteofrance.com), we had the opportunity to freely access to some of the forecasts of the French operational limited area model ALADIN [START_REF] Radnoti | The spectral limited area model ARPEGE/ALADIN[END_REF][START_REF] Yessad | Basics about ARPEGE/IFS, ALADIN and AROME in the cycle 37 of ARPEGE/IFS[END_REF][START_REF] François | The Météo-France NWP system: description, recent changes and plans[END_REF]. ALADIN forecasting results with confidence intervals in order to give more complete information to a final user, such as a power manager.

The paper is organized as follow. Section 2 describes the data we have used: radiation time series measured from meteorological stations and forecast data computed by the ALADIN numerical weather model. After recalling the principles of time series forecasting and the need to make stationary a time series, we present in section 3 the forecasting models (ARMA and ANN) that have allowed us to build our original hybrid method, and the variables selection approach used. The section 4 includes final results and experiences conducted during this study and showing that forecasting results can significantly be improved by selecting ANN or ARMA models according to their performances. The reliability of the predictions is also considered by computing the confidence intervals. The section 5 concludes and suggests perspectives.

Radiation time series and Numerical Weather Prediction forecast data

In this study we have used two types of data: radiation time series and meteorological forecasts from the ALADIN NWP model. In order to verify the robustness of our approach we chose to apply our methodology on five distinct stations located in Mediterranean coastal area. Figure 1 shows the location in Mediterranean area of the five weather stations studied. The second types of data we decided to use are the meteorological forecasts from the ALADIN NWP model. Météo-France proposed us a free access to some of the forecasts issued from their numerical weather prediction model called ALADIN-France. This model is a bi-spectral limited area, based on the assimilation of daily measurements, and driven using, for boundary data, the outputs of the ARPEGE (acronym of "Action de Recherche Petite Echelle Grande Echelle") global model providing also by French meteorological services. The model evolves in average every six months because the ALADIN code follows the ARPEGE one in its permanent evolution and we are actually on the 37 th cycle. For a better description of the model and its parameterization, the interested reader may refer e.g. to [START_REF] Radnoti | The spectral limited area model ARPEGE/ALADIN[END_REF][START_REF] Yessad | Basics about ARPEGE/IFS, ALADIN and AROME in the cycle 37 of ARPEGE/IFS[END_REF][START_REF] François | The Météo-France NWP system: description, recent changes and plans[END_REF]. The French NWP system is organized around the production of analyses at 00, 06, 12 and 18UTC, and the range of the forecast is 54 hours. The horizontal resolution of ALADIN-France is approximately 9.5 Km, with 60 levels vertically. These data are those used actually, that is to say in the 37 th cycles. The ALADIN model has more than twenty outputs available at a high temporal resolution of one hour but not the global radiation. These values are computed in all points of the computing grid with mesh size of 9.5 Km. Considering these facts we had three major choices to do.

First we had to select the forecast parameters to add as an input of our model. In a second time we had to choose the grid points for our five locations, and finally we add to select the analyses (between 00, 06, 12 and 18UTC) and ranges (1 to 54h) of the forecasts to take into account. For the ALADIN output, we based our choice on preceding works [START_REF] Voyant | Optimization of an artificial neural network dedicated to the multivariate forecasting of daily global radiation[END_REF] in which we analyzed the benefit of taken into account exogenous variables. In these studies we made some computations about the correlation between the global radiation and a lot of exogenous meteorological parameters. Among the 23

ALADIN possible outputs we chose those which seem to have a straight link to solar radiation [START_REF] Voyant | Optimization of an artificial neural network dedicated to the multivariate forecasting of daily global radiation[END_REF].

These data are pressure (P, Pa), nebulosity (N, Octas), rain precipitations (RP, mm) and temperature (T, °C). The second choice was on the selection of the grid points. We decided to simply use a proximity criterion and we chose the points of the ALADIN computational grid nearer to our five selected stations. We decided to take into account the analyses closer to the sunrise: we select the 06AM analyze with a 12 hours head forecast horizon. This couple of analyze and range allows us to cover, with a unique ALADIN extraction per day, the central hours we consider, that is to say 8:00AM to 4:00PM. However, we are aware that some of our choices can be discussed: first the periodic changes in the ALDIN model (more of ten cycles) may affect the training capability of the ANN, secondly an explicit help of a professional forecaster may comfort us in the choice of NWP model parameters, we will return to these restrictions in our conclusion.

Forecasting models

This section details the prerequisites that have enabled us to build our prediction model. After recalling the principles of time series forecasting and the need to make stationary a time series and how to validate it, we present the ARMA and ANN models. To conclude this section we present our final prediction model as a combination of ARMA and ANN models.

a. Time series analysis, need of stationarity and validation

There are different approaches to model Time Series (TS) [START_REF] Gooijer | 25 years of time series forecasting[END_REF]: ARMA [START_REF] Elminir | Prediction of hourly and daily diffuse fraction using neural network, as compared to linear regression models[END_REF][START_REF] Leisch | On the stationary of autoregessive neural network models[END_REF], ANN [START_REF] Sfetsos | Univariate and multivariate forecasting of hourly solar radiation with artificial intelligence techniques[END_REF][START_REF] Qi | Trend Time Series Modeling and Forecasting With Neural Networks[END_REF], are the ones we studied and which are the most effective based on our previous works. However the common base of all these models seems to be that a TS x t can be defined by a linear or non-linear model called f n (see Equation 1where t = n,n-1,…,p+1,p with n, the number of observations and p the number of parameters ; n p) [START_REF] Kim Tae Yoon | Artificial neural networks for nonstationary time series[END_REF][START_REF] Bourbonnais | Analyse des séries temporelles en économie[END_REF].

( ) (Eq 1)
To estimate the model, a stationarity hypothesis is often necessary. This condition usually implies a stable process [START_REF] Qi | Trend Time Series Modeling and Forecasting With Neural Networks[END_REF][START_REF] Kim Tae Yoon | Artificial neural networks for nonstationary time series[END_REF][START_REF] Hwang | Least squares estimation for critical random coefficient first-order autoregressive processes[END_REF][START_REF] Azadeh | An integrated fuzzy regression algorithm for energy consumption estimation with non-stationary data: A case study of Iran[END_REF]. This notion is directly linked to the fact that whether certain feature such as mean or variance change over time or remain constant. In fact, the time series is called weak-sense or weakly stationary if the first and the second moments are time invariant. In other word, if the first moment is constant and if the covariance is not time dependent like show on equation 2 [START_REF] Bourbonnais R Gis | Analyse des s ries temporelles : application l conomie et la gestion. 2 éd[END_REF][START_REF] Bourbonnais | Analyse des séries temporelles en économie[END_REF][START_REF] James | Time series analysis[END_REF].

[ ] ( )

and [ ] [( )( )] ( ) (Eq 2)
Note that an equivalent stationarity criterion must be fond with the simple correlation coefficient (corr). The relation linking the two parameters is

[ ] [ ] .
A stronger criterion is that the whole distribution (not only the mean and the variance) of the process does not depend on the time. The probability distribution F of the stochastic process x t is invariant under a shift in time. In this case the series is called strict stationary [START_REF] Azadeh | An integrated fuzzy regression algorithm for energy consumption estimation with non-stationary data: A case study of Iran[END_REF], the two moments shown in the Equation 2are stationary, but an other condition is also necessary (Equation 3).

( ) ( ) (Eq 3)
The stationarity hypothesis is an important tool in classic time series analysis. As it is primordial for the ARMA method, this rule stays also correct for neural network studies [START_REF] James | Time series analysis[END_REF]. In fact, all artificial networks are considered like functions approximation tools on a compact subset of . Moreover standard MLP (with at least 1 hidden layer) are asymptotically stationary, it converges to its stationary distribution, (i.e.

). Moreover, this kind of network can approximate any continuous and multivariate function. They cannot show "explosive" behavior or growing variance with time [START_REF] Leisch | On the stationary of autoregessive neural network models[END_REF][START_REF] Kurt | Multilayer feedforward networks are universal approximators[END_REF].

In practice a varying process may be considered to be close to stationary if it varies slowly and it is the modeling condition to use the MLP. Note that, the network can be trained to mimic a non-stationary process on a finite time interval. But the out-of-sample or prediction performance will be poor. Indeed, the network inherently cannot capture some important features of the process. Without pre-process, ANN and ARMA can be unappealing for many of the non-stationary problems encountered in practice [START_REF] Pollock | A handbook of time-series analysis, signal processing and dynamics[END_REF]. One way to overcome this problem is to transform a non-stationary series into a stationary (weakly or stronger if possible) one and then model the remainder by a stationary process.

In our case, we have developed a sophisticated method to make the global radiation stationary (X t , modeling of cycles). The original series has two periodicities (angular frequency and ) very difficult to overcome (Equation 4) because the amplitude a(t), b(t) and the cloud occurrence, have not simple expressions. Without the last term, this equation represents a clear sky estimation.

( ( ) ( ) ( ) ( )) (Eq 4)
The first periodicity is a classic yearly seasonality which can be erased with a ratio to trend: multiplicative scheme induced by the nature of the global radiation series on the Equation 5. In previous studies [START_REF] Christophe | Forecasting of preprocessed daily solar radiation time series using neural networks[END_REF], we have demonstrated that the clear sky index obtained with Solis model [START_REF] Mueller | Rethinking satellite-based solar irradiance modelling: The SOLIS clear-sky module[END_REF][START_REF] Pierre | A broadband simplified version of the Solis clear sky model[END_REF] is the more reliable for our locations. As the second daily seasonality is often not completely erased after this operation, then we use a method of seasonal correction (corrected for seasonal variance) based on the moving average [START_REF] Bourbonnais R Gis | Analyse des s ries temporelles : application l conomie et la gestion. 2 éd[END_REF][START_REF] Bourbonnais | Analyse des séries temporelles en économie[END_REF]. The chosen method is essentially interesting for the case of a deterministic nature of the series seasonality (true for the global radiation series) but not for the stochastic seasonality [START_REF] Zhang | Neural network forecasting for seasonal and trend time series[END_REF]. The steps we follow to make the series stationary are as following:

1-use of Solis model to establish the clear sky model of the considered location :

( ) ( ) ( ( ( )) ) ( ( )) (Eq 5)
2-calculate the ratio to trend to overcome the periodicity, the result is the clear sky index CSI :

( ) TS without trend as it is the case for the insolation. As matter of fact its annual average value is relatively constant at long term. In case of the periodic effect is significant, we can consider that the TS has a seasonal component. This test is proposed for the seasonality of the TS: daily and yearly phenomena. Before to perform this test it is necessary to construct a new scheme of the series, so , with (a and b are respectively related to the daily and yearly seasonality), where ( ) The test is constructed from the null hypothesis H0 meaning no periodic influence in the TS, and the alternative hypothesis H 1 meaning TS with periodic influence. While details can be found in [START_REF] Bourbonnais R Gis | Analyse des s ries temporelles : application l conomie et la gestion. 2 éd[END_REF], principle of the Fisher test is described by the following equations:

(
1-Compute of the "empirical" Fisher statistic ⁄ (V p the variance on period on the Equation 10 and V R the residue variance on the Equation 11). p is the number of measures per period and N the number of periods.

∑ ( 〈 〉 〈 〉 ) (Eq 10) ( )( ) ∑ ∑ ( 〈 〉 〈 〉 〈 〉 ) (Eq 11) 〈 〉 ( ) ∑ 〈 〉 ( ) ∑ (Eq 12) 〈 〉 ( ) ∑ ∑
For the a time series (daily effect) the parameter N=6 and p=3285, for the b series (yearly effect), N=6*365=2190 and p=9 (not 24 hours because only the hour between 8:00AM and 04:00PM are considered).

2-determination of the critical Fisher value for the degree of freedom v 1 = (p-1) et v 2 = (N-1).(p-1) and the = 0,05, so the value corresponding to:

(Eq 13)

3-If the H 0 hypothesis is rejected and H 1 is accepted. If The F c calculated from the data is greater than the critical value of F distribution for desired false-rejection probability (0.05), the TS is described as seasonal. Furthermore, we can estimate than, more the F c coefficient is important, and more the seasonality component is important.

In section dedicated to the experiments we will present the results obtained for the variation coefficient and the Fisher test on our CSI stationary processing with and without periodic coefficients for the five places studied.

b. ARMA

The ARMA method is certainly the most used with the prediction problems [START_REF] James | Time series analysis[END_REF][START_REF] Abrahart | Neural Network vs. ARMA Modelling: constructing benchmark case studies of river flow prediction[END_REF][START_REF] Bourbonnais R Gis | Analyse des s ries temporelles : application l conomie et la gestion. 2 éd[END_REF]. The ARIMA techniques are especially reference estimators in the prediction of global radiation field. It is a stochastic process coupling autoregressive component (AR) to a moving average component (MA).

This kind of model is commonly called ARMA (p, q) and is defined with p and q parameters (Equation 14).

( ∑ ) ( ∑ ) ( ) (Eq 14)
Where, x t is a time series, φ and θ are the parameters of the autoregressive and moving average part, L is the lag operator and is an error term distributed as a Gaussian white noise. The optimization of these parameters must be made depending on the type of the series studied. In the presented study, we chose to use Matlab© software and the Yule-Yalker fitting method [START_REF] James | Time series analysis[END_REF]. The criterion adopted to consider when an ARMA model 'fits' to the global radiation time series is the normalized root mean square error described by the Equation 15 [START_REF] Ahlburg | Error measures and the choice of a forecast method[END_REF].

√ [( ̂) ] √ [ ] (Eq 15)

The prediction error is generated by the prediction of two years of radiation not used during the ARMA parameters calculation step. Several experiments are needed to obtain the best model. Residual auto-correlogram tests have been computed to verify that the error term is a white noise. Before to use this method of forecasting, the global radiation time series is made stationary with the method described in the section 3-a (clear sky index with seasonal adjustment) and then, centered and reduced.

The models after optimization are very simple, the ARMA(1,0) for Ajaccio, Bastia, Montpellier and Nice or ARMA(2,0) for Marseille are successful. In the case of the ARMA(2,0) the prediction can be expressed by the Equation 16.

̂( ) ∑ ( ) ( ) (Eq 16)

c. Neural network and time series forecasting

Although a large range of different architecture of ANNs is available [START_REF] Hontoria | Generation of hourly irradiation synthetic series using the neural network multilayer perceptron[END_REF][START_REF] Yu | Handbook of neural network signal processing[END_REF], MultiLayer Perceptron (MLP) remains the most popular [START_REF] Balestrassi | Design of experiments on neural network s training for nonlinear time series forecasting[END_REF][START_REF] Mohamed | Radial Basis Function Network-based prediction of global solar radiation data: Application for sizing of a stand-alone photovoltaic system at Al-Madinah, Saudi Arabia[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 [START_REF] Zhang | Neural network forecasting for seasonal and trend time series[END_REF][START_REF] Zhang | Forecasting with artificial neural networks: The state of the art[END_REF][START_REF] Wei | Neural Network Earnings per Share Forecasting Models: A Comparative Analysis of Alternative Methods[END_REF] have validated this approach based on ANN for the non-linear modeling of time series. To forecast the 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] Costa | Improving generalization of MLPs with sliding mode control and the Levenberg-Marquardt algorithm[END_REF][START_REF] Faraday | Times Series Forecasting with Neural Networks: A Case Study[END_REF]. Considering the initial time series equation (Equation 1), we can transform this formula to the non-linear case of one hidden layer MPL with b related to the biases, f and g to the activation function of the output and hidden layer, and to the weights. The number of hidden nodes (H) and the number of the input node (In) allow to detail this transformation. (Equation 17):

̂ (∑ )

(Eq 17)

(∑ )

In the presented study, the MLP has been computed with the Matlab© software and its Neural Network toolbox. The characteristics chosen and related to previous work are the following: one hidden layer, the activation functions are the continuously and differentiable hyperbolic tangent (hidden) and linear (output), the Levenberg-Marquardt learning algorithm with a max fail parameter before stopping training equal to 5. This algorithm is an approximation to the Newton's method [START_REF] Zhang | Forecasting with artificial neural networks: The state of the art[END_REF] and is represented by the Equation 18(J() is the jacobian matrix, ( ) this transposed and e() the error between the N simulations and the N measures).

( ) ( ) [ ( ) ( ) ] ( ) ( ) ( ) ( ( ) ( ) ( ) ( )
) The input variables selection step is one of the key tasks of ANN optimization. We proposed to base the input selection on the use of a regression model. The methodology used by our team in the past was to compute all the coefficients of correlation between exogenous data at time t, t-1, etc. and the clear sky index at time t+1. We encounter problems related to the difficulty to find a significant limit to this coefficient. To work around the apparent classical student T-test permissiveness and to respect the parsimony principle, we fixed empirically a limit not really justified theoretically. In the presented paper, we proposed a pre-input layer selection method in order to choose from a pool of available data.

(Eq 18) ( ) ( ( ) ( ) ) ( ̂
We have used the parsimony principle and limited the approach to 10 lags for the endogenous case and 2 lags for the each ALADIN output parameters. These 18 data in total are grouped in the pre-input Layer and are shown on the scheme of the Figure 2. This method also uses the statistical Student Ttest to guide the choice but this time it is applied on the coefficients of a multiple linear regression.

Note that concerning the ALADIN data, we consider that the physical signification of the lags upper than one are not very relevant because they are predictions of anterior time steps and not measures. In first we generate a regression model for the 18 pre-inputs nodes. There are 10 endogenous nodes (In = 10), and 2 for each meteorological variable called the total number of exogenous nodes is 8). The Equation 19shows the simplest case of 10 endogenous nodes and Me nodes of exogenous parameter E [START_REF] Hwang | Least squares estimation for critical random coefficient first-order autoregressive processes[END_REF]:

( ) ̂ ( ) ( ) ∑ ∑ (Eq 19)
This equation can be expressed by matrix form like shown in the Equation 20 or in the developed form in the Equation 21. Y is a vector with the elements ( ), S' a matrix with concatenated endogenous and exogenous data and a vector with all the adjustable parameters of the regression.

For the estimation of this model the is N T data in the measurement history. The optimization corresponds to solve this equation to find correct values of linear regression weights.

The least square method is used. The classical estimator is defined by ̂ ( ) ( is the transpose of the matrix ). The dimension of the vector W LR is 19 (10 variables for the insolation, 8 for exogenous data, and 1 for the regression constant). The next step is to verify if some weights (related to each variable CSI*, N, RP, P and T called X i for 1<i<12) are not significantly different from zero. We use a Student T-test for all coefficients. The statistic used (t j ) corresponds to ( j is the standard deviation of the parameter j). It is often easier to interpret it with the notation √ ( ) ( is the error of prediction). The weight is not statistically equal to zero if t i is higher than the value 1.96 (large sample with 5% alpha level). Note that this test can be replaced by an equivalent one, based on the confidence interval of the weight considered ( [ ]). If the sign of the two limits are different, then we can assimilate the weight  j to zero. This is equivalent to search the sign of the produce of the two limits defined by (( ) (

)) (( ))). The value must be strictly positive in order to consider  j different to zero. The variables with an associate weight equal to zero are the lest correlated with the insolation at the lag t+1. It is the criterion chosen to select (or not) the variable in the input layer of the neural network. In fact our pre-input layer selection method can be resumed by the following rule:

. The application of this rule allows reducing the dimension of the input layer.

After this step, it is necessary to optimize the number of hidden nodes. The technique used is relatively standard; it consists to try several configurations by varying the number of nodes. In previous experiments [START_REF] Christophe | Forecasting of preprocessed daily solar radiation time series using neural networks[END_REF], we have seen that often the number of hidden nodes must be comparable to the number of input nodes. After optimization (choice of input and hidden nodes, activation function, etc.)

the output of the network can be expressed by the following expression (Equation 22) with ̂ ̂ ̂ ̂ and by the number of prediction of each meteorological data use (respectively p, n, rp and t). These predictions done with the model ALADIN depend on latitude, longitude, orography, temperature, humidity, etc. [START_REF] Radnoti | The spectral limited area model ARPEGE/ALADIN[END_REF].

̂ ( ) ∑ ( ∑ ( ) ∑ ̂( ) [ ] [ ] ∑ ̂( ) [ ] ∑ ̂( ) ∑ ̂( ) [ ] [ ] ) (Eq 22)
Note that some elements of the Matrix W 1 ( ) are equal to zero after the methodology related to the linear regression presented previously. The computing time is about 5 minutes with a Core i5 3.20GHz and 8Go of RAM. This computing time includes the clear sky generation, the seasonal adjustment, the statistical test, the ANN training and the average prediction computed for five ANNs with different initializing.

d. Hybrid methodology

As in primary studies we have demonstrated that performances of ARMA and MLP were similar (in univariate case), we decide to investigate the possible contribution of a hybrid methodology combining Solis model and seasonal adjustments, ARMA and ANN models with ALADIN forecast data. The idea of hybrid methodology was born from the observation of some cases, where the model ARMA, the model ANN and in generally the single predictors are not really efficient. It will be the case if an important accuracy is required or when the priorization (during a benchmark) of predictors is impossible. The solution to combining methods to understand and model the dynamic of the signal is then considered: it's the base of the data mining technics. Among the predictors hybridization, we can mention the ARMA-MLP [START_REF] Durdu | A hybrid neural network and ARIMA model for water quality time series prediction[END_REF][START_REF] Zhang | Time series forecasting using a hybrid ARIMA and neural network model[END_REF], the fuzzy inferences-MLP [START_REF] Sfetsos | Univariate and multivariate forecasting of hourly solar radiation with artificial intelligence techniques[END_REF], the wavelet-MLP [START_REF] Mellit | An adaptive wavelet-network model for forecasting daily total solar-radiation[END_REF],

the Markov-MLP [START_REF] Mellit | A simplified model for generating sequences of global solar radiation data for isolated sites: Using artificial neural network and a library of Markov transition matrices approach[END_REF] or Bayes-MLP [START_REF] Lauret | Bayesian neural network approach to short time load forecasting[END_REF]. Our choice is oriented to the ARMA-PMC (multivariate) method with time series made stationary. Concerning this model, we have used the classical approach shown previously but the predictor is now formed by two sub-predictors. There is a multitude of possible arrangement to construct it. We have opted to a model selection based on the transition linear/non-linear of the global radiation process. Linearity test exists (Lagrange multipliers), but are not really efficient when the two phenomena are concomitant in the same series. Our hypothesis is very simple, it consists to consider that when the nebulosity is low, the series is of linear character and the days with an important cloud occurrence, the series is of non-linear character. Thus, for the linear approach we use ARMA and for the non-linear approach, we select the PMC estimation.

The validation of the hypothesis is done in an in-press article not yet published.. For Zhang [START_REF] Zhang | Time series forecasting using a hybrid ARIMA and neural network model[END_REF] a hybrid model having both linear and nonlinear modeling abilities could be a good alternative for predicting time series data. By combining different models, different aspects of the underlying patterns may be captured. To construct the model selection we could only considered the day position in the year: during the winter months PMC would be selected and ARMA during the summer months. In this case, we would be considered a selection by rigid seasonality. This approach will not consider the sunny day in winter or the cloudy days in summer (although they are few). Consequently, we decided to design a methodology based on a not rigid, not repetitive and not well-marked seasonality. The new model is constructed as a stochastic model, which depends solely on the mistake made the previous hour. If the ARMA method was better at time t, (equivalent to sunny period, and indirectly to a linear phenomenon) then it will be still ARMA at t+1 else it will be the MLP. We can summarize this method with the following rule (Equation 23) where is the residue of the prediction:

| ( )| | ( )| ̂( ) ̂ ( ) ̂( ) ̂ ( ) (Eq 23)
The Figure 4 details the process followed by the hybrid methodology. The next section presents, results obtained using this hybrid model on the five places located in Mediterranean area presented in section 2.

Experiments, results and discussion

This section includes all the experiments and results conducted during this study. The first part of this section is dedicated to the stationary efficiency of our preprocessing methodology. The second part shows intermediate results on MLP trained with endogenous and ALADIN forecast data. Finally we presents results obtained with the hybrid model presented in the previous section.

a. Experiences about series makes Stationary

As we have seen previously in section 3.a, the different methods dedicated to the time series The VC X parameter is not the more interesting as a stationarity criterion but it is very simple to use.

The two other parameters show that the CSI * processing is the most efficient. There is only one case where the F c is not minimal with the CSI * (Bastia, F c(daily) =2.04). Indeed, with this method, 8 of 10 While the CSI * processing seems the most interesting, the CSI seems graphically give the first moment constant (average fixed around the year ~-0.4), but let the second central moment about the mean of the series with periodicity (variance more important during winter than summer). We can see that the curve related to the CSI * process seems the most non-seasonal, contrary to the CSI process. In winter the standard deviation has not increased. The results presented in the next section will argue if the stationarity of the series increases the prediction quality.

b. Results with MLP and ALADIN forecast data

This subpart proposes to analyze the impact of adding data from the ALADIN numerical model to a MLP. We have chosen to compare it with four other models (see Table 2). The first model which is considered as the reference is based on ARMA, the second is an endogenous MLP with our clear sky index preprocessing (CSI), the third is an endogenous MLP with CSI and seasonal correction (CSI * ), and the last is an MLP with ALADIN forecast data. In the following, we adopt a canonical form to present the MLP architecture obtained after optimization steps: (endo e , RP rp N n ,P p ,T If we analyze the annual error, the methodology V, based on the utilization of an optimized MLP with ALADIN forecasting data and seasonal adjustment with periodic coefficients (CSI * preprocessing) is the most relevant. If we analyze seasonally the errors, we find only four cases (Ajaccio, Bastia and Nice during spring and Montpellier during summer) which have a better forecaster: ARMA method.

For the other cases, the ARMA method seems equivalent to the ANN method without pretreatment (method II) but is worse than the ANN when the stationary of the series is the same (method III).

While the distinct adding of CSI * preprocessing (method III ; error is 13.0%) and ALADIN forecast data (method IV error is 13.2%) to the ANN model improves the prediction, their combined adding (method V) potentiates this effect. Because the sun light is abundant and there is low occurrence of clouds, the summer season is the most favorable to prediction.

c. Results with the hybrid methodology

This subpart proposes to present results of the hybrid method (see the Table 3) defined in section 3.d.

The ARMA and MLP models used are those presented in the previous section: models I and V. This methodology is based on the prediction of "ARMA" and "ANN" with ALADIN data, clear sky index and seasonal adjustments. The ratio ARMA/ANN represents the numbers of hourly ARMA simulations versus ANN simulations We can observe that the distribution of predictor use is 1/3 for ARMA and 2/3 for ANN. This observation confirms the results previously presented. In summer, the error is for the first time lower than 10% (Marseille). All the seasons and cities benefit of the use of this method coupling ARMA and ANN. The summer is the season where the gain is the less significant, certainly because it is the season where the stand-alone methodologies (ANN and ARMA) Ultimately, We can see that all the methodologies proposed in this paper have decreased the prediction error. In average on the five cities, the total nRMSE decrease of 11.3% against a naïve persistence predictor which has a average nRMSE equal to 26.2%. A step very interesting is the mixing between the ANN and ARMA, but do not forget that this result is only true with the type of ANN and ARMA considered. In the Figure 6, we can see the matching between measure and simulation (hybrid method presented previously) for all the cities. We can presume that there is a strong correlation between the two quantities. Except to Bastia, the cloudy periods (like the 71°-81° hour interval) seem correctly predicted. Concerning the case of Bastia, the model is not able to anticipate the nebulosity, and improvements are certainly necessary. [START_REF] Hwang | Least squares estimation for critical random coefficient first-order autoregressive processes[END_REF] This index is necessary to judge the relevance of the prediction. For example, the Figure 7 shows for Ajaccio the prediction and the confidence interval (average ± CI). Three models are evaluated: the ANN with only endogenous data, the ANN with ALADIN and the periodic coefficients and the hybrid model "ARMA" with "ANN+ ALADIN+PC". The period considered is one of the most complicated to forecast. It corresponds to the winter with a lot of cloudy days. It is interesting to see that the simple model endogenous ANN describes well the radiation, however, there are some outliers points (15 th or 41 th hour), with a confidence interval not significant. In the second curve related to the ANN+ALADIN+PC, the error is most regular, the atypical points seems non-existent. The third curve (both of ANN and ARMA) is visually the most interesting, despite that the confidence interval is yet incorrect for some points. In fact, with the hybrid method (ARMA+ANN), the ANN forecast is preponderant during the cloudy days; the ARMA predictions are used only during the sunny days. So the two last methods presented are relatively similar for the month considered.

5.

( ) | ( )| [ ] (Eq 24) ( ) ( ) ∑ ( ( ) )) [ ] (Eq
6.

Conclusion

We proposed in this paper an original technique to predict hourly global radiation time series using meteorological forecasts from a numerical weather prediction model. We optimized a Multi-Layer Perceptron (MLP) with ALADIN forecast data and endogenous data previously made stationary. We have used an innovative pre-input layer selection method and we have combined our optimized MLP to an Auto-Regressive and Moving Average (ARMA) model from a rule based on the analysis of hourly data series. This model has been successfully used to forecast the hourly global horizontal radiation for five places in Mediterranean area. In addition, this paper has allowed determining a stationarity process (method and control) for the global radiation time series. This step is primordial to correctly forecast the future values of insolation (annual average nRMSE gain over the five locations equals to 1.7%). The use of ALADIN forecasting data as an input of a MLP has shown a really great interest to improve the prediction (average nRMSE gain of 0.7%). These results could certainly be improved by a better comprehension of the complexity of the model and collaboration with a professional forecaster of Meteo-France who could help us in selection of data especially concerning the runs to consider. In last, the use of a hybrid method coupling ANN and ARMA predictors decrease much the prediction error (average nRMSE gain to 3.5%). If we compare the results with a standard model like persistence method (a typically naïve predictor) the nRMSE error is reduce by about 11.3%. The detail of the prediction error decrease following the different steps of the forecasting methodology is done on the Figure 8. The Last important point treated in this paper is the proposition of using confidence intervals in order to estimate the reliability of the prediction. The perceptive of this work are related the demonstration of the generalization of our model. We would show that the superior performance of this model is not likely to be a consequence of data mining (or data snooping). In fact, it should be sure that the model constructed in this way is not of limited practical value. Because of the intermittent nature of some energy sources (like PV or wind energy), they are included in a limited way in power systems. This limitation ensures the electrical grid stability. There are two methodologies to overcome this limitation problem: the storage of the overflow to redistribute it adequately and the prediction of the energy sources. This last solution allows to use at the right time other energy productions not dependent on weather, so to avoid a drop of the electrical current or maybe a black-out. The main problem of the management is the dispatch between these two energy types. This study deals with the one hour forecasting horizon. It is certainly the most important horizon for an electrical manager because it corresponds to the starting delay of some conventional energy sources like diesel engines or gas turbines for power generation. Actually the industrially prediction of the global radiation is often done with a simple persistence. Even if, the predictor is effective in sunny days, the present study shows that more successful forecasters exist. This result is especially interesting in the insular case (or generally to overcome the isolation problem), where the interconnection is limited and where the energy autonomy must be considered at medium or longterm.

In future studies, it would be very relevant to study an approach for solar irradiance forecasting 24 hours ahead using several MLP connected, to decrease the time step and to test the methodology on a real PV module. In addition, an important work will be to simplify the model while keeping an acceptable prediction. Indeed, even if our method is attractive, it could be complex and costly to implement for an electric power manager in grid stability context of power supply mixing renewable and conventional energy.

List of captions

List of tables

(

  http://www.cnrm.meteo.fr/aladin) is a hydrostatic model developed by Météo-France in collaboration with the European Centre for Medium Range Weather Forecasts (ECMWF). It is the result of a project launched in 1990 by Météo-France with the aim of developing a limited area model and today 15 countries are participating in the common work. We propose in this paper an original technique to predict hourly global radiation time series using meteorological forecasts from the ALADIN NWP model. After optimizing our MLP with ALADIN forecast data and endogenous data previously made stationary with an ad-hoc method, we combine it to an Auto-Regressive and Moving Average (ARMA) model from rules based on the analysis of hourly data series. Finally we present all

Figure 1 .

 1 Figure 1. The five studied stations marked in the Mediterranean sea: Ajaccio, Bastia, Montpellier, Marseille and Nice.

̂)

  In our case the parameter m takes the value 0.1 and 0.001 when the error, respectively, decreases or increases. Inputs are normalized on {-0.9,0.9} because the truncated interval gives better results than the full interval for the conditions of this study. Training, validation and testing data sets were respectively set to 80%, 20% and 0% (Matlab parameters). These phases concern the four first years and the global solar radiation test forecasting the two last years. The prediction methodology used in this paper is presented on the Figure3. The pre-input layer consists on endogenous data of radiation and meteorological forecasts (exogenous) from the ALADIN numerical weather prediction model. We can see the possibility to select only the interesting data in the input layer among endogenous data of radiation and exogenous data from ALADIN model.

Figure 3 :

 3 Figure 3 : Scheme of the prediction methodology based on MLP model and NWP model

Figure 4 :

 4 Figure 4: Scheme of the hybrid method forecasting (the L-box is the lag or backshift operator)

Figure 5 :

 5 Figure 5: Effects of the stationarity processing on the Bastia time series during 365 days (oct 2002oct 2003). The CSI * and the CSI are normalized.

Figure 6 :

 6 Figure 6: Comparison between measured and simulated global radiation done with the hybrid methodology "ARMA" and "ANN+ALADIN+PC+CSI".The 250 hours shown corresponding to half-February to half-March 2008)

Figure 7 :

 7 Figure 7: Prediction marked with confidence interval of the global radiation in Ajaccio during the 50 th first hours of January 2007. The line is related to the measure.

Figure 8 :

 8 Figure 8: Prediction error decrease following the different steps of the forecasting methodology

Figure 1 .

 1 Figure 1. The five studied stations marked in the Mediterranean sea: Ajaccio, Bastia, Montpellier,

Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 2345 Figure 2 : Scheme of the stationarity methodology

Figure 6 :Figure 7 :

 67 Figure 6: Comparison between measured and simulated global radiation done with the hybrid

Figure 8 :

 8 Figure 8: Prediction error decrease following the different steps of the forecasting methodology

Figure 2 : Scheme of the stationarity methodology

  To validate the stationary processing, we have chosen two tools. The first is the variation coefficient of

	H	) 1 (  t	X	t	1 	X	t	H gh	(t	)
	gh
			( R a t i o t o t r e n d ) -1							R a t i o t o t r e n d C le in d e x a r s k y
	3-calculate the moving average (MM(t)) considering that 2. Pre-treatment Post-treatment 1 t  CSI	Eq 6) is equal to the t CSI
	periodicity of the series. In the case of a 9 hours periodicity, 〈 〉 [ ( R a t i o t o m o v i n g a v e r a 1 ) -e g	corresponds to 4 hours R a t i o t o m S e a s o n a l a d ju o v i n g a v e r a g e s tm e n ts
	1  CP t									C	t
					*				*
			CSI	t		1	CSI	t
				Predictor
	the series (	√ [(	) ] [ ]

] (Eq 7) 4-operate a new ratio to trend with the moving average. The new coefficient are called the periodic coefficients C(t) : (Eq 8) 5-Compute of the average over one year E[ ] (done on 365x9=3285 values because only 9 hours per day have been considered) 6-Compute of the new stationary series ( ) [ ] (modulo 3285 hours) (Eq 9)

The Figure

2

shows the stationarity methodology concerning the global radiation time series prediction. ) ⁄ . Although this tool is easy to use, it is only dedicated to the cross comparison. It is not an absolute criterion allowing to distinguish the stationary and nonstationary process. To overcome this problem, there are a lot of available non-stationary tests: variance test of Fisher, unit root, Dickey-Fuller, KPSS, just to name a few. For more details, the reader can refer to Hamilton, 1994,

[START_REF] Pollock | A handbook of time-series analysis, signal processing and dynamics[END_REF] Bourbounais, 2008 [27,29,31]

. The quality of each method seems equivalent and some tests are not adapted to the seasonal case but to the extra-seasonality trends or differency stationary. We have chosen a very classic statistic test based on variance analysis: the Fisher test detailed by Bourbounais in 2008 (variance ratio to period by residue). This test involves a

Table 1 : Stationary efficiency of CSI and CSI * preprocessing for the five cities studied. Bold values signify that the stationarity is effective.

 1 forecasting have to be used with stationary time series. It should be emphasized that if the available data are nonstationary, they must be made stationary before applied to ARMA or ANN model. Based on this observation, it is necessary to study stationary methods and to see if it should be interesting in the case of the global radiation time series forecasting. The two tested processes, the clear sky index (CSI) and the CSI with use of seasonal correction and periodic coefficients (CSI * ) are compared to the original series without treatment. To know if the preprocessing makes the series stationary, we use criterions like the variation coefficient (VC X ) considered as a dispersion rate and the Fisher T-test related to the daily and yearly stationarity. The values given on the Table1are established after a normalization process between {-0.9,0.9}. According to the Fisher Table, we found that threshold of

	the F-Test is	in the
			Ajaccio	Bastia	Montpellier	Marseille	Nice
	Original series	VC F c(yearly) F c(daily)	-3.45 12.99 6.81	-2.48 12.44 1.01	-3.36 11.12 7.02	-3.62 15.07 6.53	-3.21 13.72 7.27
		VC	-0.40	-0.4	-0.37	-0.38	-0.42
	CSI	F c(yearly)	3.25	3.05	3.39	4.08	3.44
		F c(daily)	2.06	2.07	2.07	2.07	2.08
		VC	-0.09	-0.08	-1.42	-0.76	-1.2
	CSI *	F c(yearly)	0.83	0.87	0.72	0.70	0.81
		F c(daily)	2.05	2.04	1.41	1.80	1.52

  , p and t are the numbers of neurons activated in relation with endogenous data, precipitation, nebulosity, pressure and temperature, H and S are the number of hidden and output nodes. The exact description of the compared methods is: I. ARMA+PC : The best ARMA model with the CSI * preprocessing; II. ANN : The best endogenous MLP with the CSI preprocessing. Optimizations have been done with a standard method based on the interpretation of PACF partial autocorrelation factor[START_REF] Voyant | Optimization of an artificial neural network dedicated to the multivariate forecasting of daily global radiation[END_REF]; III. ANN+PC: The best MLP with CSI * preprocessing. Optimizations have been done with a standard method based on the autocorrelation factor and cross-correlation;

			III. ANN+PC (Endo 1,2,3 )x15x1	19.3	20.4	18.8	16.0	19.8
			IV. ANN +ALADIN (Endo 1:3,6:10 Pr 1,2 N 1,2 P 1,2 T 1,2 )	19.3	20.3	18.6	16.8	18.1
			V. ANN +ALADIN+PC (Endo 1,5,10 Pr 1,2 N 1,2 P 1,2 T 1,2 )	18.6	20.1	17.9	15.5	19.2
			Persistence	25.3	32.9	25.3	20.0	32.3
			I. ARMA(2,0)	18.9	23.9	19.0	11.8	21.4
	Marseille	II. ANN (Endo 1-10 )x10x1 III. ANN+PC (Endo 1,2,3,4 )x15x1 IV. ANN +ALADIN (Endo 1,2,6:10 PR 1,2 N 1,2 P 1,2 T 1 )	19.0 16.9 17.4	22.5 20.6 20.4	20.7 17.8 18.5	11.3 10.5 10.4	18.8 17.1 16.4
			V. ANN +ALADIN+PC (Endo 1,3,7 PR 1,2 N 1,2 P 1,2 T 1,2 )	16.3	19.6	16.6	10.3	16.4
			Persistence	26.4	32.1	24.5	21.1	37.1
			I. ARMA(1,0)	20.7	23.5	17.6	12.4	37.5
	Nice	II. ANN (Endo 1-10 )x10x1 III. ANN+PC (Endo 1,2,3 )x15x1	20.9 20.1	21.7 20.5	18.8 19.1	11.7 11.4	32.3 30.9
			IV. ANN +ALADIN (Endo 1,3,7:10 PR 1 N 1 P 1,2 T 1,2 )	20.1	20.4	18.1	11.4	32.3
			Models V. ANN +ALADIN+PC (Endo 1:5 PR 1,2 N 1,2 P 1,2 T 1,2 )	Annual 19.4	Winter 19.9	Spring 18.6	Summer 11.0	Autumn 30.1
			Persistence	25.1	34.7	25.2	21.4	33.9
			I. ARMA(1,0)	19.4	29.4	17.7	14.3	26.8
	Ajaccio	II. ANN (Endo 1-10 )x15x1 III. ANN+ PC (Endo 1,2,3,4 )x15x1	20.3 18.6	27.2 25.3	20.4 18.4	13.7 12.2	24.1 24.2
			IV. ANN +ALADIN (Endo 1,2,5,6,10 PR 1,2 N 1,2 T 1,2 )x15x1	19.0	26.8	19.1	12.3	23.0
			V. ANN +ALADIN+PC (Endo 1;5,9,10 PR 1,2 N 1 P 1,2 T 1,2 ) x15x1	17.8	24.9	18.1	11.7	22.0
			Persistence	27.1	35.0	27.1	22.6	34.4
			I. ARMA(1,0)	21.1	26.7	20.3	15.8	26.9
	Bastia	II. ANN (Endo 1-10 )x10x1 III. ANN+PC (Endo 1,2,3,4 )x15x1	22.8 20.8	27.3 24.9	23.4 21.4	16.1 14.9	25.7 24.9
			IV. ANN +ALADIN (Endo 1:3,8:10 PR 1 N 1,2 P 1,2 T 1,2 )	21.3	25.8	21.7	15.1	24.0
			V. ANN +ALADIN+PC (Endo 1:3,5 PR 1,2 N 1,2 P 1,2 T 1,2 )	19.9	24.4	20.5	14.2	23.3
	Montpellie	r	Persistence I. ARMA(1,0) II. ANN (Endo 1-10 )x10x1	26.9 20.1 20.8	32.6 23.5 22.4	25.9 18.7 20.2	24.6 15.5 17.9	33.2 21.9 19.3

t ) x H x S, where e, winter rp, nIV. ANN+ALADIN: An optimized MLP with the CSI preprocessing and pre-input layer selection method using ALADIN forecast data in input;

V. ANN+ALADIN+PC: Same as previous but with CSI * preprocessing.

Table 2 : Comparison between five forecasting methods for the five cities studied. The error is the nRMSE and best results are in bold.

 2

Table 3 : Annual and seasonal results for the hybrid ARMA/ANN model

 3

Table 1 :

 1 Stationary efficiency of CSI and CSI * preprocessing for the five cities studied. Bold values signify that the stationarity is effective.

Table 2 :

 2 Comparison between five forecasting methods for the five cities studied. The error is the nRMSE and best results are in bold.

Table 3 :

 3 Annual and seasonal results for the hybrid ARMA/ANN model

Acknowledgment

This work was partly supported by the Territorial Collectivity of Corsica. We thank the French National Meteorological Organization (Météo-France) and the CNRM (Centre National de Recherches Météorologiques) and particularly Ms. Delmas who has supervised the data collection from their data bank for the five synoptic stations and from the ALADIN model.


J Jin

S Leesirisan

M Song
email: m.song@lboro.ac.uk

Electrical Conductivity of Ion-doped Graphite/Polyethersulphone Composites

Keywords: Polyethersulphone, graphite, doping, ions, electrical conductivity 1. Introduction

The electrical conductivity of polyethersulphone (PES) insulating polymer was improved by incorporation of electrically conductive graphite and ions. An initial conducting pathway of the PES/graphite composites was formed at lower than 3 wt% of the filler content. LiCl was found to be an effective dopant for the improvement of the electrical conductivities of the PES/graphite composites. By doping with 0.06 wt% of LiCl the electrical conductivity was enhanced by 2 orders of magnitude. The enhancement resulted from intercalation of Li + ions into interlayer spaces of the graphite. Upon intercalation, an acceptor GIC, Li-GIC, was consequently formed. The stability of the improved electrical conductivities of the composites contributed with doped ions was assessed. The electrical conductivity of both undoped and doped graphite/PES composites slightly increased with increasing temperature and slightly decreased by physical ageing. The enhancement of the electrical conductivities by doping ions was stable at the high temperatures.

Conductive polymer composites based on insulating polymer matrices and formulated using electrically conductive fillers such as carbon black, carbon fibers, graphite, metal powder and aluminium flakes [1] have been widely investigated in academia and industry because of their outstanding multifunctional properties when compared to the pristine polymer used [START_REF] Bourdo | Graphite/Polyaniline (GP) composites: Synthesis and characterization[END_REF][START_REF] Kohjiya | Visualisation of carbon black networks in rubbery matrix by skeletonisation of 3D-TEM image[END_REF][START_REF] Yui | Morphology and electrical conductivity of injection-molded polypropylene/carbon black composites with addition of high-density polyethylene[END_REF][START_REF] Jdger | Scaling of the viscoelasticity of highly filled carbon black polyethylene composites above the melting point[END_REF][START_REF] Gojny | Evaluation and identification of electrical and thermal conduction mechanisms in carbon nanotube/epoxy composites[END_REF] They have been applied in various fields including light emitting diodes, batteries, electromagnetic shielding, antistatic coating or fibers, gas sensor, corrosion protection and activators [START_REF] Margolis | Conductive Polymers and Plastics[END_REF][START_REF] Rupprecht | Conductive Polymers and Plastics in Industrial Applications[END_REF][START_REF] Chung | Electromagnetic interference shielding effectiveness of carbon materials[END_REF].

Normally, the electrical conductivity of composite materials and as well their mechanical properties are influenced by dispersion, aspect ratio, shape, orientation of conductive fillers in insulating host matrices. In most cases, relatively large quantities of fillers are needed to reach the critical percolation value, as the filler particle size is at micrometer scales. Too high concentration of the conductive filler could lead to materials redundancy and detrimental mechanical properties [START_REF] Pinto | Conducting aluminum-filled nylon 6 composites[END_REF]. The electrical conductivity of composites, filled with micron sized carbon fibres, with the lowest aspect ratio [START_REF] Pinto | Conducting aluminum-filled nylon 6 composites[END_REF] loaded as high as 60 wt% was the similar as compared to 10wt% of the one with highest aspect ratio (300) [START_REF] Bourdo | Graphite/Polyaniline (GP) composites: Synthesis and characterization[END_REF][START_REF] Kohjiya | Visualisation of carbon black networks in rubbery matrix by skeletonisation of 3D-TEM image[END_REF][START_REF] Yui | Morphology and electrical conductivity of injection-molded polypropylene/carbon black composites with addition of high-density polyethylene[END_REF][START_REF] Jdger | Scaling of the viscoelasticity of highly filled carbon black polyethylene composites above the melting point[END_REF][START_REF] Gojny | Evaluation and identification of electrical and thermal conduction mechanisms in carbon nanotube/epoxy composites[END_REF]. Graphite is one type of carbon based filler and its electrical conductivity is 10 4 S/cm at ambient temperature as high as carbon fibers [START_REF] Ho | Handbook of Carbon, Graphite, Diamond and Fullerenes: Properties, Processing and Applications[END_REF][START_REF] Walsh | [END_REF]. Due to being layered structure of graphite, its layers are possibly expanded and exfoliated with treatment to form high aspect ratio graphite with thickness in nano-scale. In case of polymer/graphite composites, a percolation threshold depends on degree of expansion and exfoliation of graphite. The more surface area and the higher aspect ratio of graphite are, the lower amount is needed to form a continuous network. Recently, various useful methods have been applied to achieve these purposes. With thermal shock, expanded graphite can be obtained from the graphite intercalation compound [START_REF] Yasmin | Mechanical and thermal properties of graphite platelet/epoxy composites[END_REF][START_REF] Chen | PMMA/graphite nanosheets composite and its conducting properties[END_REF]. Sonication technique has been applied to this field for exfoliation of graphite [START_REF] Chen | PMMA/graphite nanosheets composite and its conducting properties[END_REF][START_REF] Wang | Preparation of polymer/oriented graphite nanosheet composite by electric fieldinducement[END_REF]. Enhancement of interfacial interactions between graphite and polymers has been successfully done through functionalisation of graphite surfaces [START_REF] Leesirisan | Polyethersulphone/Graphite conductive composites for coatings[END_REF][START_REF] Weng | Preparation and characterizations of nanoparticles from graphite via an electrochemically oxidizing method[END_REF]. Introduction of graphite into polymer matrices assisting by these methods induces volume conductivity in the range of 10 - 6 ~10 -4 S/cm at appropriate loading usually above 15-17 wt% prepared by melt blending processing [START_REF] Chen | Dispersion of Graphite Nanosheets in a Polymer Matrix and the Conducting Property of the nanocomposites[END_REF][START_REF] Dhakate | Expanded graphite-based electrically conductive composites as bipolar plate for PEM fuel cell[END_REF]. These improved conductivity properties are due to the synergistic effects of the dispersed graphite in the polymer [START_REF] Zhang | Flammability and thermal stability studies of styrene-butyl acrylate copolymer/graphite oxide nanocomposite[END_REF]. Those composites will be a benefit to coatings for electrostatic dissipation applications if its volume electrical conductivity goes up in to the range of about 10 -9 -10 -5 S/cm or possessing surface conductivity of 10 -10 -10 -6 S/cm. However, at these high contents of graphite the optimum mechanical properties of the composites were still unexpected [START_REF] Chen | Dispersion of Graphite Nanosheets in a Polymer Matrix and the Conducting Property of the nanocomposites[END_REF][START_REF] Dhakate | Expanded graphite-based electrically conductive composites as bipolar plate for PEM fuel cell[END_REF].

Ion-doping, which is a quite simple method, is often used for preparation of polymer electrolytes [START_REF] Armand | Second International Meeting on Solid Electrolytes[END_REF][START_REF] Aravindan | Investigations on Na + ion conducting polyvinylidenefluoride-co-hexafluoropropylene/poly ethylmethacrylate blend polymer electrolytes[END_REF][START_REF] Starkey | Plasticizer interactions with polymer and salt in propylene carbonate-poly(acrylonitrile)4ithium triflate[END_REF]. This method containing dissolved metal salts has been proposed as an alternative to liquid electrolyte in rechargeable batteries and enhancement of conductivity of polymer electrolyses [START_REF] Armand | Second International Meeting on Solid Electrolytes[END_REF][START_REF] Aravindan | Investigations on Na + ion conducting polyvinylidenefluoride-co-hexafluoropropylene/poly ethylmethacrylate blend polymer electrolytes[END_REF][START_REF] Starkey | Plasticizer interactions with polymer and salt in propylene carbonate-poly(acrylonitrile)4ithium triflate[END_REF]. Recently, only a few of articles were reported that the improvement of electrical conductivity for conductive filler/polymer composites was cooperated with ion-doping method [START_REF] Sang | Influences of doping approach on conductivity of composite alkaline solid polymer electrolyte PVA-HA-KOH-H2O[END_REF][START_REF] Pan | Electrical and structural analysis of conductive polyaniline/polyacrylonitrile composites[END_REF]. However, there was no any report regarding graphite/polymer systems. In graphite, each carbon has four valence electrons. Three of them are each in three sp 2 hybridised orbitals forming  bond to three nearest neighbour carbon atoms. The remaining electron delocalises in  orbital forming a  bond to a carbon atom in adjacent carbon layer. By bonding with intercalates, the graphite  bonds can gain electrons from or lose electrons to intercalates. In other words, intercalates can act as electron donors or acceptors in doping graphite. [START_REF] Dresselhaus | Lattice-dynamical model for alkali-metal-graphite intercalation compounds[END_REF]. If this method of adding ions takes into account for graphite/polymer systems, it could give rise to decrease in loading of filler and increase significantly the conductivity of the conductive composite. It may also aid in the improvement of the mechanical properties of conductive composites. In this communication we attempted to fabricate a series of electrical conductive graphite/polymer composites by combining doping ions. Polyethersulphone (PES), which is a very attractive engineering plastic polymer because of its superb properties [START_REF] Ulrich | Introduction to Industrial Polymer[END_REF], was selected as a model material. The electrical conductivity of PES composites was studied and their conduction mechanism was understood. For practical applications it needs to know whether the improved electrical conductivity of the composites contributed by doped-ions is stable when the materials are processed. So the effects of temperature and physical ageing on the electrical conductivities of the composites were also investigated.

Experimental

Materials

Expandable graphite flakes were supplied by Qing Dao Graphite Company (China).

Polyethersulphone (PES), molecular weight about 10,000 was obtained from Cytec Engineered Materials Limited (UK). Lithium chloride (LiCl), sodium chloride (NaCl) and potassium chloride (KCl) purchased from Aldrich Company (UK). N, N-Dimethylformamide (DMF), nitric acid and sulfuric acid were purchased from Fisher Scientific Company (UK).

Preparation of nitric treated graphite

First, expanded graphite was prepared by subjecting expandable graphite flake to a preheated furnace at a temperature of 850C for 2 mins. Then expanded graphite (1 g) was suspended in 40 ml of concentrated nitric acid. The mixture was heated and maintained the temperature at 60C using a magnetic stirrer hot plate for 4hrs. The resulting material was then filtered and washed with deionised water until pH 7 was reached and dried in a vacuum oven at a temperature of 100C for 2 days.

Preparation of PES/ graphite and ion-doped PES/graphite composites

PES/graphite and ion-doped PES/graphite composites were fabricated by solution method. PES and nitric treated or untreated graphite with different weight contents (1, 3, 5, 8 or 10 wt% related to PES) were dissolved in DMF (80 wt% related to PES), then, were heated and maintained at a temperature of 50C using a stirrer hot plate for 24 hrs. After that the solvent was separated by adding methanol into the solution and stirring for 15 mins at room temperature. The resulting materials were dried in a vacuum oven at 155C. Then, the composites were press moulded (240C, 5 Ton) between heated metal plates with lined aluminium foils into sheet form with a thickness of about 0.3 mm. For ion-doped PES/graphite composites with small amount of LiCl (0.02-0.06, 0.08 and 0.1 wt% related to graphite) in DMF was added into 5 wt% of untreated or nitric treated graphite before blending with PES polymer.

For comparison of conductivity the composites with the (NaCl) and (KCl) were also prepared, respectively. The contents of NaCl) and (KCl) were the same molar weight as the LiCl.

Electrical conductivity measurements

Four-Wired Digital Multimeter Model Solartron 7150 was used to measure electrical volume resistances of samples and silver paint as an electrode material. Press moulded samples (240 o C, 5 Ton) were cut into 1.9 x 1.9 cm 2 specimens with thickness about 0.3 mm, then, were applied silver paint all over top and bottom of the specimen surfaces and were dried in air. The measurements were conducted at room temperature or at elevated temperatures with 60 sec of electrificationa. The electrical volume resistances of physical aged samples were also measured. The samples were heated to 240 o C to eliminate thermal history, followed by aged at three different temperatures of T g -3 o C, T g -6 o C and T g -9 o C (T g is the glass transition temperature of the PES), each for three different aging times of 15, 40 and 60 min, respectively. The volume resistivity, v  , can be derived from the following Equation [START_REF][END_REF]:

  v v R t A   (1)
Where:

A = Effective area of the measuring electrode, cm 2 t = Average thickness of the specimen, cm v R = Measured volume resistance, k

The reciprocal of the volume resistivity is the volume conductivity.

Results and Discussion

Effect of addition of graphite on electrical conductivity of PES

Volume resistance, v R , of PES/untreated and PES/treated graphite composites with filler content 1-10 wt% prepared by the solution method was measured at room temperature by four-wired resistance measurement method. Their volume conductivities, v  , are plotted against filler content as shown in Figure 1. For very low filler contents (1-2.5 wt %), the volume conductivity of both systems are close to that of the PES matrix [29.30]. There is a sharp increase in the volume conductivity of PES/untreated graphite composites in the filler content ranging between 2.5 -3.0 wt%. Similarly, a sharp increase occurs at the same range of content for PES/treated graphite composites system. Electrical conductivity of both the PES/untreated and PES/treated graphite composites systems then level off with further increase in the amount of each filler. The former exhibited electrical conductivities about one order of magnitude higher than the latter in the filler content ranging 3-10 wt%. The electrical conductivity of the PES/untreated and PES/treated graphite composites at the filler content of 10 wt% are 12 and 11 orders of magnitude higher than that of the pure PES, respectively. PES consists of diphenylene sulphone groups and ether oxygen atoms. Bonds in the PES are strong covalent, its absence of mobile electrons makes it an electrical insulator, in contrast to electrical conductive feature of graphite [START_REF] Ram | Fundamentals of Polymer Engineering[END_REF].

The electrical conductivity of a polymer can improve with the addition of graphite.

A percolation threshold is able to achieve at lower amount of graphite, when graphite with high aspect ratio is employed. In this study, the percolation threshold of PES/untreated and PES/treated graphite composites is lower than 3 wt% which is much lower than that of composites prepared with conventional (micro-scale) graphite [START_REF] Pinto | Conducting aluminum-filled nylon 6 composites[END_REF]. The improvement in conductivity of the composites at low loading could be attributed to high aspect ratio of the graphite. The results also indicate that nitric treatment of graphite is not beneficial in electrical conductivity enhancement of the composites since the chemical modification of the graphite surfaces by oxidation leads to a reduction of the aspect ratio of the graphite [START_REF] Gojny | Evaluation and identification of electrical and thermal conduction mechanisms in carbon nanotube/epoxy composites[END_REF].

Improvement of electrical conductivity by ion doping

The addition of dilute LiCl content in the PES/treated graphite composites results in great improvement in the electrical conductivity as shown in Figure 2. The electrical conductivity gradually rose to 10 -5 S/cm and to 10 -4 S/cm with a LiCl content of 0.04% and 0.06%, respectively, from 10 -6 S/cm without the dopant. In contrast, the improvement in the electrical conductivity of the PES/untreated graphite composites was insignificant. With increasing LiCl content up to 0.06 wt%, the electrical conductivities were still in the same order of magnitude as the corresponding undoped composite.

Therefore, it is necessary to understand the effect of LiCl-doping on the electrical conductivities of the PES composites. An attempt to further increase the electrical conductivities of the PES/treated graphite composite was done by increasing the amount of LiCl. As expected, the volume conductivities of the composite changes with presence of LiCl. However, it was interesting to find out that the volume conductivity of the composite increased with increasing LiCl content, only when the content was very low (0.02-0.06 wt%). Beyond these concentrations, the volume conductivities show a downward trend (Figure 3). compound (Li-GIC) as follows [START_REF] Noel | Electrochemistry of graphite intercalation compounds[END_REF]:

  Li C x  Li C x   (2)
It is believed that there is no chloride ions intercalation due to their negative charges.

Lithium ions accept  electrons from the graphite, according to the direction of charge transfer, Li-GIC is considered as an acceptor GIC [START_REF] Ahuja | Calculated electronic and optical properties of a graphite intercalation compound: LiC6[END_REF][START_REF] Shioyama | The interactions of two chemical species in the interlayer spacing of graphite[END_REF]. Upon the intercalation, extra holes are created in  orbitals of carbon atoms in the graphite and Fermi energy un-treated treated is correspondingly lowered. The density of holes in Li-GIC could be higher than that in the pristine graphite. An electron from adjacent carbon atoms may hop into a hole leaving a new hole behind. Holes are thus mobile. In an electric field, holes are accelerated toward the electric field, which is in the opposite direction to that in which accelerated electrons are moving [START_REF] Ahuja | Calculated electronic and optical properties of a graphite intercalation compound: LiC6[END_REF][START_REF] Shioyama | The interactions of two chemical species in the interlayer spacing of graphite[END_REF]. The force reacting on holes in electric field is equal to eE , where e is the electric charge and E is the electric field. It is the high density and the high mobility of holes that lead to the higher conductivity of the doped composites.

Further increase in the dopant concentration introduces more holes on  orbitals of carbon atoms on graphite layers. Too high a density of holes has been reported to be accompanied by a decrease in their mobility due to a collision broadening [35]. The collision of electrons results in the change of their direction of motion by scattering.

Holes thus move more randomly toward the electric field. This suggests the existence of a dilute limit of lithium chloride content which can effectively improve the conductivities of the PES/treated graphite composite.

In order to show that the enhancements of the electrical conductivities of LiCl-doped composites is attributed to charge transfer between intercalate and graphite during intercalation, the effect of the size of intercalates and mixing time on the electrical conductivity were investigated.

Another two kinds of metal chloride, namely, NaCl and KCl were used as dopants for the PES/treated graphite composites. The basic idea is that an intercalation will occur easier for the intercalate having a smaller size, which will contribute to more charge transfer. The results indicate that the conductivities of each type of doped composites at the same molar content of dopants were in the same order of magnitude.

However, the conductivity of the LiCl-doped composites was found to be a little higher than that of the NaCl and KCl -doped composites. The conductivity of the composites containing KCl was found to be the lowest.

Na + and K + ions which are derived from NaCl and KCl behave as electron acceptors the same way as Li + ions do, i.e. they form Na-GIC and K-GIC and accept  electrons from the graphite upon intercalations. The sizes of neutral atoms of Li, Na and K are 0.145, 0.180 and 0.220 nm [36], respectively, and the original interlayer spacing of the graphite obtained from the XRD result is 0.343 nm [START_REF] Leesirisan | Polyethersulphone/Graphite conductive composites for coatings[END_REF]. There is a difference in the capability of the individual intercalates to diffuse into the original interlayer spaces of the graphite as reflected in the electrical conductivities of each doped composite system. Among the three intercalates, Li + ion is the smallest, and Li-doped PES/treated graphite composites correspondingly possesses the highest conductivity.

Mixing time of 0.08% of LiCl and the treated graphite before adding into the PES was varied from 3 hrs to 1, 3 and 7 days. The basic idea is that if intercalation processes occur, then, electrical conductivities of the doped composites will be correspondingly affected by mixing time. On the contrary, if there is no intercalation, the electrical conductivities of the doped composites should be independent on the mixing time employed. The results are presented in Figure 4.

The LiCl-doped composites prepared with 3 hrs and 1 day mixing time have the same level of electrical conductivity which are the two highest value among the four values. The volume conductivity decreases continuously with the prolonged times.

The difference between the lowest conductivity obtained from the doped composites prepared with 7 days mixing time and the highest conductivity values is as large as 3 orders of magnitude. The results imply that only a certain amount of Li + ions make the greatest improvement in the electrical conductivities and are consistent with

Stability of Electrical conductivity of doped graphite/PES composites

For material processing it is very important to know whether the electrical conductivities of the PES/graphite composites enhanced by doped ions are stable. The electrical conductivities of the LiCl-doped PES/TG5% composites with LiCl contents of 0.02, 0.06 and 0.1 wt% at room and elevated temperatures were measured and compared with undoped ones. The results are listed in Table 1. The electrical conductivity of the doped composites slightly increased with increasing temperature as same as the undoped ones. This means that the electrical conductivity enhanced by the dopant was stable at the high temperatures. The temperature dependence of the electrical conductivity of conducting composites having conductive particles dispersed in insulating matrices is normally expressed as [START_REF] Murugaraj | Electron transport in semiconducting nanoparticle and nanocluster carbon-polymer composites[END_REF]:

] ) exp[( ) ( 0 0    T T T   (3)
Where:

 = Conductivity, S/cm or S/m 0  = Proportionality constant T = Absolute temperature, K 0 T = a constant  = a constant
From Equation 3, a relationship:

  ) 1 ( ) ( ln T T 
is obtained. The electrical conductivity data shows a good fit with a value of  = 1/4 in the plot of

) ( ln T  versus 4 1  T
as shown in Figure 5. A value of  = 1/2 has been reported for many systems such as polyimide/carbon nanoparticle composite [START_REF] Long | Electrical properties of multi-walled carbon nanotube/polypyrrole nanocables: percolationdominated conductivity[END_REF] and muti-walled carbon nanotubes/polypyrrole nanocables [START_REF] Long | Electrical properties of multi-walled carbon nanotube/polypyrrole nanocables: percolationdominated conductivity[END_REF]; whereas  = 1/4 has been reported for carbon black filled polyimide nanocomposites [START_REF] Ne | Temperaturedependent transport properties in the semiconducting regime of nanoparticle carbon-polyimide composite films[END_REF].

1  T (K -1/4 )
The measured electrical conductivity values of composites at elevated temperature are mainly contributed by two factors which are coefficient of thermal expansion (CTE) of polymer matrix and conductive particles, and thermal activation of holes and electrons [START_REF] Murugaraj | Thermistor behaviour in a semiconducting polymer-nanoparticle composite film[END_REF]. The latter is the more prominent factor at elevated temperatures.

Ln((T))

PES/TG5% PES/TG5%/LiCl0.02% PES/TG5%/LiCl0.06% PES/TG5%/LiCl0.1%

1/T 1/4 (K -1/4 )
The difference between the CTE of the PES polymer matrix and the graphite contributes to a decrease in the electrical conductivity value. The CTE of the PES is 55  10 -6 /C [41] while that of graphite within the measured temperature range (20-190C) is about 22 10 -6 /C in the c-direction but is negative in the ab-direction [START_REF] Ho | Handbook of Carbon, Graphite, Diamond and Fullerenes: Properties, Processing and Applications[END_REF].

The coefficient of thermal expansion of materials is defined as the change in length per unit length for a change in temperature one degree Centigrade. Materials typically have a characteristic of positive expansion coefficient i.e. expand in response to heating but contract in response to cooling. The negative expansion coefficient in the ab-direction of the graphite means that the graphite will contract on heating [START_REF] Ho | Handbook of Carbon, Graphite, Diamond and Fullerenes: Properties, Processing and Applications[END_REF]. With increasing temperature, the PES expands by taking up more space the polymer chains are moving around more vigorously. In the mean time, the graphite sheets will shrink in ab-direction but expand in c-direction. The vigorous movement of the PES polymer chains may end up being a thick barrier in a conductive path or the shrinkage of graphite sheets may make the inter-particle distance wider.

Contrary to the CTE, the thermal activation of holes and electrons contributes to an increase in the electrical conductivities of the composites. With rising temperature, there are more thermally activated holes and electrons acquiring the necessary energy for jumping or hopping across the potential barriers. Generally, barrier heights in conducting particle filled composites can be significantly reduced when the filler loading is sufficiently high [START_REF] Long | Electrical properties of multi-walled carbon nanotube/polypyrrole nanocables: percolationdominated conductivity[END_REF]. The variation of barrier heights among the different composite systems filled with 5 wt% treated graphite is due to the difference in the effective amounts of Li + ions presenting in the gallery of the graphite. With the presence of the guest Li + intercalate, a hopping conduction can more easier occur along c-axis of graphite in which the barrier height has been lowered after the intercalation.

The electrical conductivities of the doped PES/treated graphite composites and undoped PES/treated graphite composites, which were aged at three different temperatures (below 3 o C, 6 o C or 10 o C of Tg.) for different times, were measured and compared with corresponding non-aged samples in Table 2. It can be seen that at a certain ageing temperature, the electrical conductivity of each aged doped composite sample only slightly decreased compared with that of the corresponding non-aged sample and their changes were similarly as the un-doped composites. It indicates that the electrical conductivities of the PES/treated graphite composites contributed by doping ions were not affected by the physical ageing.

As known, quenching a soften material through g T causes free volume to freeze in.

The free volume will be gradually lost from the material during ageing at a temperature well below g T by rotational movement of the PES polymer side chains.

Disappearance of the free volume causes densification of the polymer which may correspondingly result in the creation of gaps at the interfaces between the PES polymer and graphite sheets. This process will make a contribution to a decrease in the electrical conductivities of the composite.

During the densification, graphite sheets are probably forced further apart, in other words, are interrupted by barriers, whether space or PES matrix. If graphite sheets are still close to the nearest neighbouring graphite sheets or are interrupted by just thin barriers, a tunnelling process can occur directly between graphite sheets [START_REF] Celzard | Conduction mechanisms in some graphite -polymer composites: the effect of a direct-current electric field[END_REF]. If graphite sheets, on the other hand, are too far from the neighbouring graphite sheets, the electrical conductivities of composites will be dramatically decreased due to the absence of tunnelling process. Nonetheless, probability of tunnelling may be enhanced due to the presence of intermediate states such as impurities or dispersed defects [START_REF] Celzard | Conduction mechanisms in some graphite -polymer composites: the effect of a direct-current electric field[END_REF]. In the doped composite systems, electrons or holes from a graphite sheet may cross a thick barrier to another graphite sheet via one or several intermediate states i.e. Li + ions within the barrier. With passing ageing time, PES polymer side chains acquire more equilibrium conformation. Therefore, localised movement of the polymer side chains and also the variation of electrical conductivity of each aged samples become less with prolonged ageing time.

In doped PES/graphite composite systems, Li + ions can assist tunnelling across the insulating polymer if they are present at the vicinity of the graphite conductive paths.

However, the ions are actually incorporated at very small amount in the doped composite system and are randomly dispersed in the graphite gallery or the polymer matrix which may or may not in the vicinity of the network of conductive paths.

In addition, a lower ageing temperature could result in a slower of the localised rotational movement of PES polymer side chains [START_REF] Ferry | Viscoelastic Properies of Polymers[END_REF]. The overall effect of physical ageing on electrical conductivity of the composites aged at lower ageing temperatures can be discussed in the same way as those aged at higher ageing temperatures.

Conclusions

The electrical conductivity of PES polymer was improved by addition of graphite and doped-ions. Percolation threshold of PES/untreated and PES/treated graphite composites was lower than 3 wt% of the filler. It was found that nitric treatment of graphite was not beneficial in electrical conductivity enhancement of the composites since surface functionalisation by oxidation reaction leaded to a lower aspect ratio of graphite. LiCl was found to be an effective dopant for the improvement of the electrical conductivity of PES/treated graphite composites. Electrical conductivity of 5 wt% treated graphite/PES composite with the presence of a LiCl content of 0.06 wt% was 3.0  10 -4 S/cm, which was 2 orders of magnitude higher than that of the undoped composite and was in the same order of magnitude as that of PES filled with 8 wt% of treated graphite. The enhancement of electrical conductivity of the doped composites was resulted from intercalation of Li + ions into interlayer spaces of graphite. Upon intercalation, an acceptor GIC, Li-GIC, was consequently formed.

The electrical conductivity of both undoped and doped PES/treated graphite composites slightly increased with increasing temperature and slightly decreased by physical ageing. The increase conductivity resulted from increased thermally activated of holes and electrons. Physical ageing leaded to gap creation process at the interfaces between PES polymer and graphite sheets which accompanied free volume annihilation and densification of the polymer. In comparison between undoped and doped PES/graphite composite systems, it is believed that the improved electrical conductivities contributed by doping ions were stable at high temperatures.

Figure 1 :

 1 Figure 1: Volume conductivity of

Figure 2 :Figure 3 :

 23 Figure 2: Volume conductivity of

Figure 4

 4 Figure 4 which shows that 0.06 wt% of LiCl is the proper content. The exceeded

Figure 4 :

 4 Figure 4: Volume conductivity of

Figure 5

 5 Figure 5 Temperature dependence of electrical conductivity for LiCl-doped

4

 4

Table 1

 1 Volume conductivity of PES/TG5% and LiCl-doped PES/TG5% composites

	at various temperatures
			Volume conductivity (S/cm)
	Temperature
	(ºC)	PES/TG5%	PES/TG5%	PES/TG5%	PES/TG5%
		/LiCl 0.00%	/LiCl 0.02%	/LiCl 0.06%	/LiCl 0.1%
	20	4.2 10 -6	4.2 10 -6	3.1 10 -4	1.4 10 -4
	50	4.3 10 -6	8.3 10 -6	3.3 10 -4	1.4 10 -4
	75	4.8 10 -6	8.3 10 -6	3.4 10 -4	1.5 10 -4
	100	5.6 10 -6	9.1 10 -6	3.6 10 -4	1.7 10 -4
	125	6.7 10 -6	1.0 10 -5	3.7 10 -4	2.2 10 -4
	150	7.1 10 -6	1.0 10 -5	4.0 10 -4	2.4 10 -4
	190	7.7 10 -6	1.1 10 -5	4.3 10 -4	2.6 10 -4

Table 2

 2 The effect of aging time on electrical conductivity of PES/graphite compositesaged at three different temperatures

	Aging		Electrical Conductivity (S/cm)
	temperature o C	Sample	0 hr	1 hr	5 hrs	24 hrs
		PES/TG5%	8.4 10 -8	2.6 10 -8	1.8 10 -8	2.3 10 -8
	T g -3	PES/TG/LiCl 0.006%	5.6 10 -6	1.4 10 -6	1.3 10 -6	1.2 10 -6
		PES/TG5%	8.4 10 -8	2.1 10 -8	2.9 10 -8	3.0 10 -8
	T g -6	PES/TG5%/LiCl 0.006%	5.6 10 -6	1.2 10 -6	1.3 10 -6	1.2 10 -6
		PES/TG5%	8.4 10 -8	2.9 10 -8	3.6 10 -8	4.2 10 -8
	T g -10	PES/TG5%/LiCl 0.006%	5.6 10 -6	1.0 10 -6	1.0 10 -6	1.0 10 -6