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 non-linear 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 K_T.

1 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⁻²) 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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }F\left(z;a,b\right)=\underset{\mathrm{0}}{\overset{x}{\int }}f\left(\mathit{\xi };a,b\right)\mathrm{d}\mathit{\xi }=\mathrm{1}-{\left(\mathrm{1}-z^{\mathit{\alpha }}\right)}^{\mathit{\beta }}\iff \\ \text{(1)}& {\displaystyle }& {\displaystyle }Q\left(z\right)=F^{-\mathrm{1}}\left(z\right)=\mathit{\{}\mathrm{1}-(\mathrm{1}-z{)}^{\frac{\mathrm{1}}{a}}{\mathit{\}}}^{\frac{\mathrm{1}}{b}},\end{array}$$

where $z\in \left[\mathrm{0},\mathrm{1}\right]$ is standardized according to $z=\frac{z-z_{min}}{z_{max}-z_{min}}$, with z_min and z_max are minimum and maximum values determined from the empirical time series.

Figure 1 Double periodicity diagram of Solar radiation (i.e. averaged measured value for each hour) for Denver station.

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)

$$\begin{array}{}\text{(2)}& K_T={\displaystyle \frac{I}{I_{\mathrm{0}}}},K_T\in \left[\mathrm{0},\mathrm{1}\right]\end{array}$$

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₀ and so it accounts for the transparency of the atmosphere.

The denominator of Eq. (2) is a deterministic process, which can be determined from

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{I}_{\mathrm{0}}={\displaystyle \frac{\mathrm{12}}{\mathit{\pi }}}{G}_{\mathrm{sc}}d[{\displaystyle \frac{\mathrm{180}}{\mathit{\pi }}}\left({\mathit{\omega }}_{\mathrm{2}}-{\mathit{\omega }}_{\mathrm{1}}\right)\sin \left(\mathit{\delta }\right)\sin \left(\mathit{\phi }\right)\\ \text{(3)}& {\displaystyle }& {\displaystyle }+\cos \left(\mathit{\delta }\right)\cos \left(\mathit{\phi }\right)\left(\sin \left({\mathit{\omega }}_{\mathrm{2}}\right)-\sin \left({\mathit{\omega }}_{\mathrm{1}}\right)\right)],\end{array}$$

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

$$\begin{array}{}\text{(4)}& d=\mathrm{1}+\mathrm{0.0034}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }J}{\mathrm{365}}}-\mathrm{0.05}\right)\end{array}$$

where J=1 and 365 refer to 1 January and 31 December, respectively;

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\omega }={\displaystyle \frac{(\mathrm{St}-\mathrm{12})h}{\mathrm{24}\cdot h}}\mathrm{360}\iff \mathit{\omega }=\mathrm{15}(\mathrm{St}-\mathrm{12}),\\ \text{(5)}& {\displaystyle }& {\displaystyle }\mathit{\omega }\in [-\mathrm{180}{}^{\circ },\mathrm{180}{}^{\circ }]\end{array}$$

where St is the solar time, $\mathrm{St}\in \left[\mathrm{0},\mathrm{23}\right]$ so between successive hours ${\mathit{\omega }}_{\mathrm{1}}-{\mathit{\omega }}_{\mathrm{2}}=\mathrm{15}{}^{\circ }$;

$$\begin{array}{}\text{(6)}& \mathit{\delta }=-\mathrm{0.49}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }J}{\mathrm{365}}}+\mathrm{0.16}\right)\end{array}$$

Therefore, investigating the stochastic nature of clearness index K_T can lead to conclusions for the fluctuations in hourly solar radiation.

2 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

3.1 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 $\mathit{\lambda }\in [\mathrm{0},\mathrm{1}]$. The probability distribution and density functions of the proposed distribution are:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }F\left(x;\mathit{\lambda },a_{\mathrm{1}},b_{\mathrm{1}},a_{\mathrm{2}},b_{\mathrm{2}}\right)=\\ {\displaystyle }& {\displaystyle }\mathit{\lambda }F\left(x;a_{\mathrm{1},}{b}_{\mathrm{1}}\right)+\left(\mathrm{1}-\mathit{\lambda }\right)F\left(x;a_{\mathrm{2}},b_{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }F\left(x;\mathit{\lambda },a_{\mathrm{1}},b_{\mathrm{1}},a_{\mathrm{2}},b_{\mathrm{2}}\right)=\\ {\displaystyle }& {\displaystyle }\mathit{\lambda }\underset{\mathrm{0}}{\overset{\mathit{\chi }}{\int }}f\left(\mathit{\xi };{\mathit{\alpha }}_{\mathrm{1}},b_{\mathrm{1}}\right)\mathrm{d}\mathit{\xi }+\left(\mathrm{1}-\mathit{\lambda }\right)\underset{\mathrm{0}}{\overset{\mathit{\chi }}{\int }}f\left(\mathit{\xi };{\mathit{\alpha }}_{\mathrm{2}},b_{\mathrm{2}}\right)\mathrm{d}\mathit{\xi }\\ {\displaystyle }& {\displaystyle }F\left(x;\mathit{\lambda },a_{\mathrm{1}},b_{\mathrm{1}},a_{\mathrm{2}},b_{\mathrm{2}}\right)=\\ \text{(7)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }\left[\mathrm{1}-{\left(\mathrm{1}-x^{a_{\mathrm{1}}}\right)}^{b_{\mathrm{1}}}\right]+(\mathrm{1}-\mathit{\lambda })\left[\mathrm{1}-{\left(\mathrm{1}-x^{a_{\mathrm{2}}}\right)}^{b_{\mathrm{2}}}\right]{\displaystyle }& {\displaystyle }f=\mathit{\lambda }{f}_{\mathrm{1}}+\left(\mathrm{1}-\mathit{\lambda }\right)f_{\mathrm{2}}\iff f\left(x;\mathit{\lambda },a_{\mathrm{1}},b_{\mathrm{1}},a_{\mathrm{2}},b_{\mathrm{2}}\right)=\\ {\displaystyle }& {\displaystyle }\mathit{\lambda }f\left(x;a_{\mathrm{1},}{b}_{\mathrm{1}}\right)+\left(\mathrm{1}-\mathit{\lambda }\right)f\left(x;a_{\mathrm{2}},b_{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }f\left(x;\mathit{\lambda },a_{\mathrm{1}},b_{\mathrm{1}},a_{\mathrm{2}},b_{\mathrm{2}}\right)=\\ {\displaystyle }& {\displaystyle }\mathit{\lambda }\left[{\mathit{\alpha }}_{\mathrm{1}}{b}_{\mathrm{1}}{x}^{(a_{\mathrm{1}}-\mathrm{1})}{\left(\mathrm{1}-x^{a_{\mathrm{1}}}\right)}^{b_{\mathrm{1}}}\right]\\ \text{(8)}& {\displaystyle }& {\displaystyle }+(\mathrm{1}-\mathit{\lambda })\left[{\mathit{\alpha }}_{\mathrm{2}}{b}_{\mathrm{2}}{x}^{(a_{\mathrm{2}}-\mathrm{1})}{\left(\mathrm{1}-x^{a_{\mathrm{2}}}\right)}^{b_{\mathrm{2}}}\right]\end{array}$$

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).

Figure 2 Examples of comparison between empirical and theoretical distribution, where the latter is a linear combination of Kumaraswamy distributions. (a, b) Psitallia station, Greece (bimodal); (c, d) Barrow station, Alaska (unimodal).

3.2 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.

Figure 3 Climacogram of three investigated stations: Denver, H=0.78, Elizabeth city, H=0.75, Bluefield Virginia, H=0.64.

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.

4 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 diurnal-seasonal 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\left(x_{i,j};{\mathbf{p}}_{i,j}\right)$ to each one of the diurnal-seasonal processes (with the estimated, e.g. $i=\mathrm{12}\times j=\mathrm{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^{-\mathrm{1}}\left(F\right(x_{i,j};{\mathbf{p}}_{i,j});\mathrm{0},\mathrm{1})$. Note that in case the diurnal-seasonal processes are all Gaussian, the proposed homogenization is equivalent to the standard normalization scheme, where for each process the mean is subtracted and the residual is divided by the standard deviation.

In this implicit way, we manage to homogenize the timeseries $x_{i,j}\sim F(x_{i,j};{\mathbf{p}}_{i,j})$ to $y\sim N(\mathrm{0},\mathrm{1})$. In case where the marginal distribution is unknown or difficult to estimate, we may use non-linear transformation schemes based on the maximization of entropy (Koutsoyiannis et al., 2008; Dimitriadis and Koutsoyiannis, 2015b). It is noted that a more robust approach to reduce the 12×24 set of parameters would be to employ an analytical expression for the double solar periodicity (as done for the wind process in Deligiannis et al., 2016). This homogenization scheme has been applied to several processes such as wind (Deligiannis et al., 2016), solar radiation (Koudouris et al., 2017), wave height, wave period and wind for renewable energy production (Moschos et al., 2017), river discharge (Pizarro et al., 2018) and precipitation (Dimitriadis and Koutsoyiannis, 2018). However, it is noted that this scheme assumes stationary in the dependence structure rather cyclostationary (for such analyses see Koutsoyiannis et al., 2008, and references therein).

Figure 5 Comparison of distribution (a, c) and density (b, d) function of simulated and observed data, for hour 13:00: (a, b) May; (c, d) September.

Figure 6 Lag-1 autocorrelation of hourly time series for each month (a) empirical; (b) simulated.

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):

$$\begin{array}{}\text{(9)}& \mathit{\gamma }\left(k\right)=\mathrm{1}/(\mathrm{1}+k/q{)}^{\mathrm{2}-\mathrm{2}H}\end{array}$$

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).

Figure 7 Climacogram of observed and simulated K_T vs. white noise and a homogenized process without periodicities.

In this implicit manner, the marginal characteristics of each periodic part are exactly preserved (since the marginal distribution functions are implicitly handled through the proposed homogenization) and the expectation of the second-order 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 preserved 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-like dependence structure, which later it can be transformed back to the original distribution function $F(x_{i,j};{\mathbf{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 determined value of the hourly intensity of solar radiation at the top of the atmosphere.

5 Conclusions

The hourly marginal distribution of solar radiation is investigated via the clearness index K_T. Regarding the marginal distribution, after analysing a variety of stations with different regional climate data, it is concluded that the Kumaraswamy distribution cannot adequately describe the empirical data of hourly solar radiation. Therefore, a new bimodal distribution is constructed which is a weighted sum of two Kumaraswamy distributions This distribution is fully characterized by five parameters and can adequately fit the empirical hourly data of the clearness index in in the regions investigated in this study. Also, the dependence structure of solar radiation process is investigated via the climacogram tool. It is concluded that since the Hurst parameter is estimated higher than 0.5, the examined process exhibits long-term persistence and cannot be considered as a white noise nor a Markov process. Finally, a new preliminary stochastic model is proposed, capable 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.