A stochastic model for the hourly solar radiation process for application in renewable resources management

Since the beginning of the 21st century, the scientific community has made huge leaps to exploit renewable energy sources, with solar radiation being one of the most important. However, the variability of solar radiation has a significant impact on solar energy conversion systems, such as in photovoltaic systems, characterized by a fast and nonlinear response to incident solar radiation. The performance prediction of these systems is typically based on hourly or daily data because those are usually available at these time scales. The aim of this work is to investigate the stochastic nature and time evolution of the solar radiation process for daily and hourly scale, with the ultimate goal of creating a new cyclostationary stochastic model capable of reproducing the dependence structure and the marginal distribution of hourly solar radiation via the clearness index KT .


Introduction
Human activities are either explicitly or implicitly linked with the dynamic behaviour of the solar radiation process.As a result, during the last two decades extensive research (Ettoumi et al., 2002;Tovar-Pescador, 2008;Reno et al., 2012;Tsekouras and Koutsoyiannis, 2014) has been carried out into the stochastic nature of the solar radiation process (e.g.marginal distribution, dependence structure etc.).Although many popular distributions, used in geophysics such as (Gamma, Pareto, Lognormal, Logistic, mixture of two Normal etc.) are suggested in the literature (Ayodele and Ogunjuyigbe, 2015;Jurado et al., 1995;Aguiar and Collares-Pereira, 1992;Hollands and Huget, 1983) and may exhibit a good fit, they cannot adequately fit the right tail of the distribution.This can be explained, considering that the right boundary of the process varies at a seasonal scale.Also, the maximum value of solar radiation that can be measured at the earth surface is the solar constant (i.e.G sc = 1367 W m −2 ) and therefore, distributions which are not right bounded should not be applied for solar radiation.Koudouris et al. (2017) introduce a new marginal distribution (i.e.Kumaraswamy distribution) for daily solar radiation process which is verified by three goodness of fit tests.The Kumaraswamy distribution may not be a very popular distribution, but it originates from the Beta family of distributions (Jones, 2009) and exhibits some technical advantages in modelling, e.g.invertible closed form of the cumulative distribution function: where z ∈ [0, 1] is standardized according to z = z−z min z max −z min , with z min and z max are minimum and maximum values determined from the empirical time series.
In this research framework a more extensive analysis is being conducted for hourly solar radiation.The marginal distribution for hourly scale and the dependence structure of the examined process are being investigated with the ultimate goal to synthesize a preliminary cyclostationary stochastic model capable of generating synthetic series calibrated from regional climate data.Furthermore, one the most common characteristics of hydrometeorological processes is the double periodicity (i.e. the diurnal and seasonal variation of the process); therefore, solar radiation exhibits same behaviour (e.g.Fig. 1).The seasonality occurs due to the deterministic movement of the earth in orbit around the sun and around its axis of rotation.
In order to proceed to the hourly scale data analysis, the double periodicity must be deduced from the measured data.To achieve that, the clearness index K T is introduced Eq. ( 2) The index can be considered as a ratio of both deterministic and stochastic elements.The clearness index K T , describes all the meteorological stochastic influences (e.g.cloudiness, dew point, temperature, atmospheric aerosols) as it is the ratio of the actual solar radiation measured on the ground I , to that available at the top of the atmosphere I 0 and so it accounts for the transparency of the atmosphere.The denominator of Eq. ( 2) is a deterministic process, which can be determined from where d is the solar eccentricity (Eq.4), G sc is the solar constant, ω is the hour angle (Eq.5), ϕ is the latitude and δ is the declination (Eq.6); these quantities are given as where J = 1 and 365 refer to 1 January and 31 December, respectively; where St is the solar time, St ∈ [0, 23] so between successive hours Therefore, investigating the stochastic nature of clearness index K T can lead to conclusions for the fluctuations in hourly solar radiation.

Experimental data
The analysis of hourly solar radiation is conducted via the clearness index K T .In order to examine the process behaviour on a world-wide scale, data from both United States of America and Greece are examined.Data for Greece were retrieved from the Hydrological Observatory of Athens (HOA).The network consists of more than 10 stations located in Attica region measuring environmental variables of hydrometeorological interest.Each station is equipped with a data logger which records with 10 min interval.These measurements were aggregated to mean hourly data through the open software "Hydrognomon" (http://hydrognomon. org, last access: 9 August 2018).For the Attica region, the K T timeseries were generated utilizing Eq. ( 3) by producing hourly solar top of atmosphere intensity at the local station.For USA, the data base of NRLE (National Renewable Energy Laboratory)-NSRDB (National Solar Radiation Database) was used, which contains more than 1500 stations with hourly solar radiation, but only 40 of them include time series with more than 14 years of measurements of hourly global solar radiation measured on horizontal surface.
3 Hourly solar radiation stochastic investigation

Marginal distribution
According to previous research (Koudouris et al., 2017), daily solar radiation can be modelled by the Kumaraswamy distribution Eq. ( 1).A new investigation for the hourly marginal distribution of hourly solar radiation is being conducted through the clearness index K T .Firstly, for every examined station, the hourly empirical data K T is divided into 288 times series (i.e.24 hourly time series for each month which are constructed with approximately 30 observations during days for a certain number of years).This technique is usually found to be sufficient allowing the construction of a stationary model for the solar radiation time series.The Kumaraswamy distribution is fitted to the empirical data of the clearness index, to evaluate the statistical properties of the solar data under study.Furthermore, three goodness of fit tests (Marsaglia and Marsaglia, 2004;Csorgo and Farway, 1996;Burnham and Anderson, 2003) were applied (i.e.Kolmogorov-Smirnov, Cramer von Misses and the Anderson Darling) by setting the significance level at 5 %.For the Psyttalia station (Greece), from the 288 times series only 170 are considered with a mean value of solar radiation much larger than zero (0, 1367], as a result of the absence of the sunshine during the night period.We calculated that for only 44 of them the Kumaraswamy distribution is appropriate.
The empirical probability distributions that were constructed for these results, seem to indicate that the clearness index and therefore the solar radiation exhibit a bimodal behaviour.As a result, from the analysis of all HOA network stations, the Kumaraswamy distribution cannot describe the empirical data sufficiently well and thus it is an insufficient distribution for modelling hourly solar radiation in the Attica region.Nevertheless, the empirical probability distributions that were constructed from NSRDB stations indicate that hourly solar radiation does not always exhibit a bimodal behaviour at any geographic location.This is due to the fact that hourly solar radiation is extremely correlated with the cloudiness process.Consequently, in regions where the cloudiness process does not exhibit fluctuations in small time scales, the Kumaraswamy distribution seem to adequately fit the hourly empirical data.This assumption is confirmed after investigating the marginal distribution of the clearness index at the Barrow station in Alaska where cloudiness does not exhibit a highly varying behaviour in contrast to the Attica region.For a better representation of the multivariate analysis scenarios of clearness index, a linear combination of Kumaraswamy distributions is proposed.A new distribution is constructed considering the sum of two Kumaraswamy distributions conventionally chosen and weighted by a parameter λ ∈ [0, 1].
The probability distribution and density functions of the proposed distribution are: The parameters are calculated via the least square estimator method or the maximum likelihood estimation.The resulting f (x) and F (x) exhibit good agreement with the empirical ones (e.g.Fig. 2),notably when solar radiation exhibits a bimodal behaviour (e.g.Fig. 2c, b).

Dependence structure
It is well known from the literature that hydrometeorological processes show a large variability (often linked to the maximum entropy; Koutsoyiannis, 2011) at different time scales, exhibiting so-called long-term persistence (LTP) or else the Hurst-Kolmogorov dynamics (Koutsoyiannis, 2002).
In order to investigate if solar radiation exhibits an LTP behaviour, the Hurst parameter (representing the dependence structure of the process) is estimated via the climacogram tool (i.e.double logarithmic plot of variance of the averaged process versus averaging time scale k; Koutsoyiannis, 2010).The Hurst parameter is identified only at large scales.Thus, in Fig. 3, it is identified for time scales > 1 year (8760 h), where the seasonal variation is averaged out and equals the half slope of the climacogram, as scale tends to infinity, plus 1.The climacogram has some important statistical advantages if compared to the autocovariance and power spectrum (Dimitriadis and Koutsoyiannis, 2015a;Harrouni et al., 2005).Exploiting the Hurst parameter, the persistence behaviour of the process is quantified and examined.The analysis of both NSRDB and HOA stations estimates the Hurst parameter larger than 0.5 (where the latter value corresponds to a white-noise behaviour); consequently, the examined process indicates long-term persistence and cannot be considered as white noise nor a Markov process.

Methodology and application of the model
In this section, we describe a simple methodology to produce synthetic hourly solar radiation.After we cautiously select a marginal distribution model (e.g., Eq. 8, Fig. 5), we estimate the distribution parameters p i,j for each diurnalseasonal process x i,j (e.g., 12 × 24 different set of parameters; Fig. 4).Then we homogenize the timeseries (Dimitriadis, 2017) by applying the distribution function F x i,j ; p i,j to each one of the diurnal-seasonal processes (with the estimated, e.g.i = 12 × j = 24, set of parameters p i,j ) and then, by employing the standard Gaussian (or any other) distribution function to each diurnal-seasonal process, i.e. y = N −1 (F (x i,j ; p i,j ); 0, 1).Note that in case the processes the homogenization standard scheme, for the is In this implicit timeseries ∼ F (x i,j ; p i,j ) to y ∼ N (0, 1).In case where distribution difficult the maximization et al., 2008; noted approach 12 × 24 of employ analytical expression the (as the Deligiannis 2016).scheme been applied to such wind (Deligiannis lar (Koudouris riod and wind for and and However, noted that this scheme rather (for such analyses see Koutet al., 2008, and references therein).The above homogenization enables the estimation and modelling of the dependence structure after the effect of the double periodicity has been approximately removed.This homogenization scheme also enables approximating the correlations among the diurnal processes for the same scale (Fig. 6).After the estimation of the N (0, 1) homogenized timeseries, we estimate the dependence structure through the second-order climacogram and we fit a generalized-Hurst Kolmogorov (GHK) model (Dimitriadis and Koutsoyiannis, 2018): where γ is the standardized variance, k the scale (h), q a scale-parameter and H the Hurst parameter (for the examined process we estimate q = 2 and H = 0.83, considering the bias effect; Fig. 7).
In this implicit manner, the marginal characteristics of each periodic part are exactly preserved (since the marginal distribution functions are implicitly handled through the pro-posed homogenization) and the expectation of the secondorder dependence structure (e.g.correlation function) is also exactly preserved after properly adjusting for bias through the mode or expected value of the estimator (Dimitriadis, 2017).It is noted that higher-order moments of processes with HK behaviour cannot be adequately in an implicit manner (see an illustrative example in Dimitriadis and Koutsoyiannis, 2018, their Appendix D) and thus, for a more accurate preservation of the dependence structure an explicit algorithm is necessary (Dimitriadis and Koutsoyiannis, 2018).We may use a simple generation scheme, such as the sum of AR(1) models (SAR; Dimitriadis and Koutsoyiannis, 2015b), that can synthesize any N (0, 1) autoregressive-G.Koudouris   like dependence structure, which later it can be transformed back to the original distribution function F (x i,j ; p i,j ) and so in this way produce a double-periodic process with the desired marginal distribution for each diurnal-seasonal cycle as well as the desired dependence structure.Finally, we multiply each value of the synthetic K T with the deterministically intensity radiation the top atmosphere.hourly distribution of solar radiation is invesvia the clearness index K T .Regarding marginal a of ferent it the distribution new distribution which is distribution characterized parameters and adequately fit the hourly index in in this study.the of solar radiation is that the parameter mated than 0.5, the process exhibits persistence noise a Markov Finally, preliminary proposed, of reproducing the clearness sky index K T and so the hourly solar radiation.The model can maintain and preserve the probability density function by means of the first four central moments and also the Hurst parameter which represents the correlation and the persistent behaviour.

Figure 4 .Figure 5 .
Figure 4. Results of the simulation model for the Denver Station: (a) 2-year synthetic timeseries of hourly K T ; (b) yearly average of hourly solar radiation observations vs. synthetic values of 15 years simulation.
et al.: A stochastic model for the hourly solar radiation process