Category Archives: EuroSun2008-13

Nonlinear Analysis of Daily Global Solar Radiation Time Series

G. Lopez1*, F. J. Batlles2 and J. L. Bosch2

1 Departamento de Ingenieria Electrica y Termica, Escuela Politecnica Superior, Universidad de Huelva,

21819, Huelva, Spain

2 Departamento de Fisica Aplicada, Universidad de Almeria, Almeria, Spain
Corresponding Author, gabriel. lopez@die. uhu. es


In this work we analyse a daily global irradiance time series into the framework of chaotic dynamic systems in order to examine a possible underlying nonlinear behaviour. Employed methods are based on a phase space reconstruction from the measured data and are devoted to the calculation of the properties of an underlying attractor, such as the Lyapunov exponents. Researches on these dynamical system invariants will point out the presence of chaos. We also use local lineal models as a test for nonlinearity.

The global solar radiation data were measured at the radiometric station of the University of Almeria (Spain) during eight years. Results have shown the non-existence of any attractor in the phase space for the global irradiance time series. Negative Lyapunov exponents exclude a chaotic behaviour that might allow a better short term prediction than autoregressive models, and the idea of the existence of a nonlinear differential equation system. These results match with those obtained from applying local linear models for prediction, of which estimations suggest that the data are best described by a linear stochastic process.

Keywords: solar radiation, forecasting, chaos, nonlinear time series.

1. Introduction

Information on the availability of solar radiation is needed in many applications dealing with the exploitation of solar energy. Particularly, global solar radiation is one of the most important input parameters for any solar energy system and different techniques have been developed to model and forecast it. Once the solar energy system (like concentrated solar power plants) is running, prediction of power load, normally done on an hourly basis with a prediction horizon between 1 and 24 hours, is instrumental for planning and operation of the total power system, e. g. for buying or selling power or for solving the unit commitment and dispatch problems.

Studies about solar radiation time series and other meteorological variables with autoregressive or stochastic models [1-3] have experienced an important growth in the last years, since synthetic sequences statistically indistinguishable from the original ones are needed to design the solar devices properly. However, using this type of model, prediction from past values is limited due to the random character exhibited by these time series. The main factor of randomness affecting global solar radiation data is due to variations in the sky cloud cover, which make the radiation on cloudy days difficult to predict.

Appearing of several nonlinear dynamical models developed by Lorenz [4], with a fully irregular and complex behaviour (generically named chaos), allows to study the nature of fluctuations in solar radiation data from a new methodology. Time series from many system evolutions, apparently evolving in a random way, can now be predicted with higher precision than with traditional ARMA models, at least for short term predictions. Detection of chaos in a time series would thus allow a better modelling of time series against statistically based models.

In this paper, we analyse a daily global solar irradiance time series into the framework of chaotic dynamic systems in order to examine a possible underlying nonlinear behaviour. Existence of chaos should help to improve the short term predictions performed by means of the traditional
statistical techniques. Employed methods are based on a phase space reconstruction from the measured data and are devoted to the calculation of the properties of an underlying attractor, such as the Lyapunov exponents. We also use local lineal models as a test for nonlinearity.

2. Experimental Data

The horizontal global solar radiation time series (symbolized by {xt}, being t the time step in days) consisted of 2945 daily values. They were obtained by integration from experimental values averaged every ten minutes during the years 1990-1992 and every five minutes during the years 1993-1998. The measurements were recorded at the radiometric station of the University of Almeria (36.8° N, 2.44° W) in south-eastern Spain by means of a Kipp and Zonen CM-11 pyranometer. Measurements from other two pyranometers (LI-COR, model LI-200) were also available. They were used to detect and replace any anomalous value measured by the CM-11 pyranometer. The calibration constants of the pyranometers were checked periodically by our research group. Measurements by the LI-200 pyranometers were also corrected following [5].

3. Methodology and results

Residual kriging model evaluation

Подпись: Fig.4. June and October H maps based on the residual kriging method on a 1 km x 1 km grid.

Again, as for the ordinary kriging, a seasonal pattern is present, with the summer months showing the lowest RMSE values and the winter months the highest ones. Particularly, lowest RMSE relative value are found for February 0.54 MJ m-2day-1 (4.5%) and the highest for October 1.31 MJ m-2day-1 (10.2 %). MAE ranges from 0.22 MJ m-2day-1 (2.2%) in January to 0.70 MJ m-2day-1 (5.4%) in October. Regarding the ME, the model tends to slightly overestimate the solar radiation, for all the months, but ME values are fairly low. Figure 4 shows the estimated H maps based on the residual kriging procedure.

. Relation between cloud index and clear sky index

The clear sky index is an appropriate and thus an established quantity for the description of the transmission properties of the Earth’s atmosphere.

It is defined by

k* = Ig / Idear (3)

where the current global radiation Ig is related to the expected clear sky radiation Idear. In addition a model is needed, which describes the atmospheric conditions, a so called clear sky model, see sec. 2.2.

[Beyer et al. 96] introduced the relation

k* = 1 — n (4)

into the Heliosat method. This relation can be used without a previous calibration with ground data, but it is necessary that the clear sky model is optimally adapted to the regarded area. By Fontoynont et al. [2] relation (4) was improved in order to detail the case of complete cloudiness and to limit the clear sky index to meaningful values.

. Ground albedo

The influence of the albedo on the SSI varies with the solar zenith angle and the atmospheric turbidity. This influence is assessed by computing the difference between the SSIs obtained with the ground albedo equal to respectively 0 and 0.2 (Fig. 5).

The influence of the albedo on the SSI decreases as the wavelength increase because the intensity of diffusion decreases. It increases with increasing atmospheric turbidity and so, decreases also with increasing visibility. This influence decreases when the solar zenith angle is growing because of the decreasing of the SSI; this is also true for the relative deviation. Presence of clouds increases the relative influence of the albedo: the larger tc, the larger the relative deviation.

Solar Radiation Dataset

The dataset used in this study is based on pyranometric global solar irradiance ground measurements which belongs to stations of the Spanish National Radiometric Network from Agencia Espanola de Meteorologia (AEMet) sited in Madrid (40.45°N, 3.72°W, 680 m.). The period of data acquired goes from 1st of January 1979 up to 31st of December 2003 having 17376 of half daily values.

The quality of pyranometric available data and temporal acquisition of measurements are two crucial factors related to adjustment of proposed models. In this study an exhaustive quality analysis has been done on both important factors which influence high quality measurements (based on BSRN recommendations). Manufacturer of pyranometers used is the Dutch firm Kipp&Zonen and belong to CM11 series.

2. Methodology

The mean daily solar radiation presents a seasonal variation due principally to the variation of the position of the earth to the sun and a stochastic aspect introduced by atmospheric components (aerosols, H2O, O3,…).

Half daily ground solar radiation shows a seasonal trend with fluctuations from day to day due to cloudiness, changing air mass and aerosols. The seasonal trend can be separated in several different ways: (1) “clearness sky index” based on extracting the influence of latitude, dividing between measured ground solar radiation and extraterrestrial solar radiation values, (2) subtraction of the annual harmonic (first harmonic of Fourier analysis), (3) “lost component” which is based on the subtraction of each value of the extraterrestrial radiation the corresponding value of ground solar radiation. In general, approach (1) leads to a procedure much simpler from the point of view of the user, but it is not suitable because statistical forecasting methods needs predictors with stationary behaviour and Gaussian distribution.

Review to solar radiation forecasting methods

Three different methods are currently proposed to forecast solar radiation. The first is based on time series models, which use a series of the daily average of the measured solar radiation as input data. It has been shown in different papers that a RMSE of 5.1 % [7], 8.3 % [8], and 8.4 % [9] can be obtained, for local predictions using ANN — wavelet methods. Whereas in [7] and [8] the Continuous Wavelet Transform (CTW) is used, [9] employs the Discrete Wavelet Transform (DWT). The author of [8] uses the day number of the year and defuzzficated cloud cover information from the weather forecast service, as auxiliary information. These models are applicable for sites where the solar radiation was measured during one year [7], two years [8], or longer time intervals [7], [9]. The second method is able to forecast the motion of clouds using satellite imaging over the earth surface. It can forecast the solar radiation for any site or area, but the uncertainty related to the utilized models increases substantially over 22% to 30%, for forecast horizons larger than six hours, as reported by Lorenz [10]. The author obtained the former value for low and the latter for high cloud variability. The third method, used in the present article, is based on the NWP with its statistical correction MOS.

Study area and observations

This study has been conducted within the Natural Park of Sierra Magina, in Southern Spain (37.6° N, 3.4° W) (Fig. 1). The area has a rough extension of 30 km x 48 km and shows a relatively complex topography, with a minimum height above sea level of 600 m and a maximum of 2150 m, with the 50% of the grids points above 800 m. The maximum surface slope is 56.94° and the 25% of the grids have a slope higher than 16.41°. The predominant surface orientation is around North. This aspects distribution benefits the existence of terrain patches with no direct solar beam along the year. The local climate is that of a typical middle and high Mediterranean mountainous environment, with dry hot summers, cold winters and relatively high precipitation during autumn and spring.

Located within the park, our research group have deployed a total of 11 radiometric stations (Figure 2). The radiometric stations are equipped, among other instrumentation, with HOBO (Onset Corporation) data-loggers that record the global irradiance registered by Licor 200-SZ radiometers. An annual calibration and inter-comparison of the sensors is carried out and the instrumental error is estimated to be less than 5%, typically 3%. The data set consists on data recorded on a 10 minutes basis and, for the purposes of this work, they were averaged to a hourly basis. In order to avoid the cosine error of the sensors, those data corresponding to a solar height below 7° were discarded. The location of the stations cover a range of elevations, slopes and aspects as wide as was possible. Table 1 shows the topographic features of the stations.


Fig. 1. Localization of the study area.







Подпись: 9Mr

Fig. 2. Location of the stations. Results at stations 4, 5 and 11 are evaluated in this study.



Elevation(m. a.s. l.)

Slope (°)

Aspect (°)







Mancha Real

















































Table 1. Topographic features of the stations. Aspect values are presented in decimal degrees measured from North, increasing clockwise. Results presented in this study correspond to stations 4, 5 and 11 (underlined).

Comparison to maps derived from satellite data

Differences between global solar radiation interpolated from ground level measurements (Fig. 2 values) and published in [7] were shown on Fig. 6. As it is observed, a non-homogeneous distribution of significant relative differences (from +25 to -60 %) indicates the high dependence of solar radiation from the local conditions in every location. Therefore, it should be recommended a calibration of satellite data from a number of ground level measurements (ideally, from all of the available stations) that can represent these local effects.

2.2 Self-learning experience

From the 34 students registered in the Solar Radiation subject, 32 started their work selecting a meteorological station to process. 31 students passed the subject in the first attempt (including a written exam by 60 % of the total qualification); just one student required to repeat the written exam, in order to pass the subject. The other two students abandon the subject.

The interest of the students in this work was higher than the theoretical issues and practical problems of the subject. Specially, processing meteorological datasets from real stations selected by the students was an extra incentive. And, in addition, processed datasets result very useful to get a solar resource evaluation in the region over the studied period, as only minor revisions of students’ calculations were required.

3. Conclusions

In this work an estimation of the annual solar resource over Galicia was developed, in the framework of a self-learning experience with the collaboration of students of the MSc in Renewable Energies and Energetic Sustainability. This educational experience was very successful, as the interest of the students was very high and their results were very useful to achieve the solar resource maps.

From the solar global irradiation distribution, solar resource in this region is extremely connected to the local conditions (local cloudiness, rain, fog), that are very variable because of the changeable weather and complex topography of this region. Although the validity of these results is limited to an annual period, irradiation distribution shows the necessity to consider local conditions in the application of any methodology to estimate solar resources at Galicia.

As the main conclusion, satellite data processing has to be completed with a high number of ground level measurements. Other techniques, as Digital Terrain Models (DTM) [17], allow to estimating ground level irradiation considering topographic effects; these technique can help to cover complex terrain zones with lack of ground level data.


Meteorological dataset provided by MeteoGalicia (Xunta de Galicia) from its web page is acknowledged. Dataset processing was partially developed by the students of the MSc in Renewable Energies and Energetic Sustainability of the University of Santiago de Compostela (2007-08 academic year).

This work was partially funded by Galician regional R&D Programme (Xunta de Galicia) under project 07REM02CT.


[1] Vera Mella N. Atlas climatico de irradiation solar a partir de imagenes del satelite NOAA. Aplicacion a la peninsula iberica. PhD Thesis. Universitat Politecnica de Catalunya, Barcelona, Spain; 2005.

[2] Sharmer k. Towards a new atlas of solar radiation in Europe. International Journal of Solar Energy, 1994; 15: 81-87

[3] Baldasano J. M., Soriano C., Flores H. Atlas de radiacio solar a Catalunya. Institut Catala d’Energia. Barcelona, Spain; 2001.

[4] Alnaser W. E., Eliagoubi B., Al-Kalak A., Trabelsi H., Al-Maalej M., El-Sayed H. M., Alloush M. First solar radiation atlas of the Arab world. Renewable Energy 2004; 29: 1085-1107

[5] Font Tullot I. Atlas de la radiation solar en Espana. Instituto Nacional de Meteorologia, Ministerio de Transportes, Turismo y Comunicaciones. Madrid, Spain; 1984.

[6] Nrnez, M., Reyes, J. J., Marroquin, A., Ramiro, A. Mapas de valores medios mensuales de irradiacion solar estimados para Extremadura a partir de otros datos meteorologicos. XXVIII Jornadas Cientlficas de la Asociacion Meteorologica Espanola. Badajoz, Spain, 2004.

[7] Vazquez Vazquez M., Santos Navarro J. M., Prado Cerqueira M. T., Vazquez Rios D., Rodrigues Fernandes F. M. Atlas de radiacion solar de Galicia. Universidad de Vigo. Vigo, Spain; 2005.

[8] Rigollier, C., Lefevre, M. and Wald, L. The method Heliosat-2 for deriving shortwave solar radiation from satellite images. Solar Energy 2004; 77: 159-169.

[9] Pettazzi A., Souto J. A. Analysis of the incoming solar radiation and other significant parameters to estimate the solar resource at eight sites in Galicia (NW Spain), Proceedings of EUROSUN 2008, 1st International Conference on Solar Heating, Cooling and Buildings., Lisbon, Portugal, 2008

[10] Salson S., Souto J. A. Automatic weather stations network of the department of environment of Galicia: data acquisition, validation and quality control, Proceedings of the 3rd international conference on experiences with automatic weather stations, Torremolinos, Spain; 2003.

[11] Pettazzi A., Souto J. A., Salson S. EOAS, a shared joint atmospheric observation site of MeteoGalicia. Proceedings of 4th ICEAWS — International Conference on Experiences with Automatic Weather Stations, Lisbon, Portugal; 2006.

[12] Davies J. A. Validation of models for estimating solar radiation on horizontal surfaces. Report available from the IEA, Downsview, Ontario, Canada; 1988.

[13] Iqbal M. An introduction to solar radiation, Academic Press, San Diego, CA; 1983.

[14] Instituto Nacional de Meteorologia. Guia resumida del clima en Espana 1971-2000. Instituto Nacional de Meteorologia, D25.3, Ministerio de Medio Ambiente. Madrid, Spain; 2001.

[15] Batlles F. J., Martinez-Durban M., Miralles I., Ortega R., Barbero F. J., Tovar-Pescador J., Pozo — Vazquez D., Lopez G. Evaluation de los recursos energeticos solares en zonas de topografia compleja. XII Congreso Iberico y VII Congreso Ibero Americano de Energia Solar. Vigo, Spain; 2004.

[16] Martinez Cortizas A., Perez Alberti A. Atlas Climatico de Galicia, Xunta de Galicia, Santiago de Compostela, Spain; 1999.

[17] Tovar-Pescador J., Pozo-Vazquez D., Molina A., Batlles F. J., Lopez G. Mejora en la estimation de la irradiancia solar en zonas de topografia compleja mediante modelos digitales del terreno. XII Congreso Iberico y VII Congreso Ibero Americano de Energia Solar. Vigo, Spain; 2004.

Phase space reconstruction


The first step in any analysis of chaotic data is to reconstruct an attractor (defined as set of points in phase space visited by a signal trajectory after transients are gone) from the data. Taken [6] suggested a method of phase space reconstruction that is known as the method of delays. This method consists of embedding the time series into a d-dimensional space, ^d, which is equivalent to the original unknown phase space composed of all the dynamical variables. Figure 1a shows the lag space plot for the global solar radiation time series using a lag of one. Points are distributed around a 45° straight line denoting a periodic nonchaotic time series. This is also evidenced observing the correlogram showed in Fig. 1b. An annual cycle of 365 days is presented. This is an expected result as the daily values of solar radiation are mainly affected by the earth motion.

xt (MJ/m2) Lags (days)

Fig. 1. a) Phase space structure of the global solar radiation time series using lags of one; b) Correlogram for

the global solar radiation time series.

In order to remove periodicity and trend, a transformation of the time series is undertaken by differencing. The new data are arithmetic differences between pairs of observations using a lag of one day (Axt = xt+I — xt). The differenced time series is denoted as {Axt}. Figure 2a shows the phase space plot for the differenced time series. The periodic pattern shown in the above figure seems to be removed. This result is corroborated analysing the correlogram displayed in Figure 2b. Autocorrelation coefficients fall to within the random-like zone after the first lag and, for the most part, remain in that zone thereafter. It is thus assumed that periodicity is removed.

Although point distribution in Fig. 2a seems to be due to some non-random complex underlying behaviour, the presence of a chaotic attractor in that phase space is not clear. To search the
optimum time delay г we have used the method of mutual information [7]. Mutual information, like autocorrelation, tries to measure the extent to which values Axt+m are related to values of Axt, at a given lag. However, mutual information has the advantage of using probabilities, rather than a linear basis, to asses the correlation and thus, nonlinear correlations are taken into accounts. The software implementation of that algorithm (and of those used hereafter) is from the TISEAN package [8].


Подпись: Fig. 3. Mutual information as a function of lag for the differenced time series.

As a prescription, a good candidate for the time delay г will be where the first marked minimum of mutual information occurs. From Fig. 3 it is noted that the first minimum is for г = 3. Once the time delay is chosen, the next step is to select the embedding dimension. Kennel et al. [9] introduce the method called the false nearest neighbors, which calculates the minimal embedding dimension. Two points in a d-dimensional phase space are false neighbors, in the sense of some distance function, when the distance, ||xjd — x/||, is small but the distance ||xi+1d+1 — Xj+1d+1|| in a d+1- dimensional phase space is not. Given a distance function, the Euclidean distance in our case, and some threshold size, the percentage of false near neighbors becomes zero as the dimension of phase space goes to the minimal embedding dimension. Figure 4 shows the percentage of false nearest neighbors as a function of the embedding dimension d. The main characteristic is that the fraction of false nearest neighbors does not fall to zero as d increases. This implies that the differenced time series has residual ‘noise’ in it, and thus the existence of low dimensional chaos is to be discarded.

Ordinary kriging vs residual kriging

The comparison of the ordinary and the residual kriging models results, in terms of the validation dataset, shows several interesting features. Residual kriging provides better estimates for all the months, with relative improvement in RMSE ranging from 5% to 20%. The maximum

improvement is found for January, with (RMSE value of 0.46 versus 0.60, 23% for relative improvement). In June, relative improvement is about 5% and in October about 9%. Similar results are found for the MAE values, with relative improvements ranging from 5% to 20%. Particularly, the maximum improvement is found during December (MAE 0.31 versus 0.4, that is 22% of relative improvement). A slightly improvement is also found for the ME values, changing from a small underestimation using the ordinary kriging to almost negligible error values using the residual kriging. Finally, R2 values show an overall improvement, being particularly high for January, February and March.

Additionally, mean values are fairly reproduced by both ordinary and residual kriging methodologies, with a maximum error (overestimation) of about 0.5 MJ m-2day-1 (2.6%) during September. On the other hand, standard deviation estimation is considerable better reproduced by the residual kriging method than by the ordinary kriging. Maximum improvements are found for the winter months For instance, for January, the observed standard deviation is 0.99 MJ m-2day-1, while the value provide by the ordinary kriging is 0.62 MJ m-2day-1 (37% error) and the estimated using the residual kriging is 0.94 MJ m-2day-1 (5% error). For the summer months, residual kriging also provides better standard deviation estimates, but improvements are lower.

Regarding the maximum values, similar estimates are provided by both kriging methods during the summer months. On the other hand, residual kriging estimates are considerable better for winter months. As far as minimum values is concerned, residual kriging perform considerable better than the ordinary kriging method for all the months. Maximum differences are found for the summer months.


Overall, the ordinary kriging method is able to provide fair estimates of the solar resources in the area of the study, with RMSE values ranges from 0.63 MJ m-2day-1 (6.2%) in June to around 1.44 MJ m-2day-1 (11.2%) in October. Nevertheless, by the inclusion of external information in the interpolation procedure, the residual kriging estimates shows considerable lower errors.

Particularly, the inclusion as external explanatory variable of the semi-sky-view factor (which accounts for topographic shadows cast and is able to explain between 15% and 45% of the spatial variability) give rise to relative improvement in RMSE values ranging from 5% in the summer month to more than 20% in the autumn and winter months. Particularly, RMSE values of the residual kriging estimates ranges from 1.44 MJ m-2day-1 (5.5%) in the June to around 1.31 MJ m — 2day-1 (10.2 %) in October. Explained variance also shows a considerable improvement compared to the ordinary kriging method, with all month showing R2 values above 0.92.


[1] J. Mubiru, K. Karume, M. Majaliwa, E. J.Banda, T. Otiti, (2006). Interpolating methods for solar radiation in Uganda. Theor. Appl. Climatol. 88, 259-263.

[2] I. O.Odeh, A. B.McBratney, D. J.Chittleborough, (1995). Further results on prediction of soil properties from terrain attributes. Geoderma, 67, 215-226.

[3] S. Ahmed, G. deMarsily, (1987). Comparison of geostatistical methods for estimating transmissivity.

Water Resour. Res. 23, 1717-1737.

[4] D. L.Phillips, J. Dolph, D. Marks, (1992). A comparison of geostatistical procedures for spatial analysis for precipitation in mountainous terrain. Agric. Forest Meteorol. 58, 1031-1047

[5] J. P.Delhomme, (1978). Kriging in the hydrosciences. Adva. Water Resour. 1, 251-266.

[6] Clark Labs, (2006). IDRISI Software, Clark University, MA, USA.