An empiric-stochastic approach, based on normalization parameters, to simulate solar irradiance

TL;DR

This paper introduces an empiric-stochastic method using normalization parameters and kernel density estimation to accurately simulate solar irradiance, aiding in photovoltaic system design and testing.

Contribution

It presents a novel empiric-stochastic approach based on normalization parameters and probability density maps for solar irradiance estimation.

Findings

Good agreement between simulated and experimental solar power data

Method effectively estimates solar irradiance for system design

Supports pre-installation testing of photovoltaic systems

Abstract

The data acquisition of solar radiation in a locality is essential for the development of efficient designs of systems, whose operation is based on solar energy. This paper presents a methodology to estimate solar irradiance using an empiric-stochastic approach, which consists of the computation of normalization parameters from solar irradiance data. For this study, solar irradiance data was collected with a weather station during a year. Post-treatment included a trimmed moving average, to smooth the data, the performance a fitting procedure using a simple model, to recover normalization parameters, and the estimation of a probability density map by means of a kernel density estimation method. The normalization parameters and the probability density map allowed us to build an empiric-stochastic methodology that generates an estimate of the solar irradiance. In order to validate our method, simulated solar irradiance has been used to compute the theoretical generation of solar power, which in turn has been compared to experimental data, retrieved from a commercial photovoltaic system. Since the simulation results show a good agreement has been with the experimental data, this simple methodology can estimate the solar power production and may help consumers to design and test a photovoltaic system before installation.

Tables

Table 1. Table 1 : Astronomical seasons definition employed in this work.

Season	from day		to day
1	35	(February 4, 2017)	124	(May 4, 2017)
2	125	(May 5, 2017)	218	(August 6, 2017)
3	219	(August 7, 2017)	309	(November 5, 2017)
4	310	(November 6, 2017)	34	(February 3, 2018)

Table 2. Table 2 : Coefficients of the long-term variable trends, to be used in Eq.  14 .

$\tilde{Y}$	$y_0$	$y_1$	$y_2$
$\tilde{A}$	1.9790	-0.0017	0.0005
$\tilde{B}$	0.8990	0.0689	-0.0089
$\tilde{C}$	913.0363	103.6416	-54.6980

Table 3. Table 3 : Number of bins M x subscript 𝑀 𝑥 M_{x} and width Δ ​ x Δ 𝑥 \Delta x for the histograms of the residuals x 𝑥 x of each parameter A 𝐴 A , B 𝐵 B and C 𝐶 C with ther corresponding general trends A ~ ~ 𝐴 \tilde{A} , B ~ ~ 𝐵 \tilde{B} and C ~ ~ 𝐶 \tilde{C} , respectively, and the distribution parameters μ 𝜇 \mu and ν 𝜈 \nu of the distributions f A subscript 𝑓 𝐴 f_{A} , f B subscript 𝑓 𝐵 f_{B} and f C subscript 𝑓 𝐶 f_{C} of the corresponding residuals x A subscript 𝑥 𝐴 x_{A} , x B subscript 𝑥 𝐵 x_{B} and x C subscript 𝑥 𝐶 x_{C} .

Normalization	Trend	Residual	No. of bins	Width	Distribution
parameter		$x$	$M_x$	$\Delta x$	$f\left(x\right)$	$\mu$	$\nu$
$A$	$\tilde{A}$	$x_A=A-\tilde{A}$	15	0.0042	$f_A\left(x\right)$	0.0028	0.0064
$B$	$\tilde{B}$	$x_B=B-\tilde{B}$	15	0.0078	$f_B\left(x\right)$	0.0048	0.0111
$C$	$\tilde{C}$	$x_C=C-\tilde{C}$	15	20.7895	$f_C\left(x\right)$	26.0092	26.0947

Table 4. Table 4 : Nominal specifications of the Solartec™  S60PC-250 solar panel [ S ( 60 ) ]. STC - Standard Test Conditions; PTC - Photovoltaics for Utility Scale Applications Test Conditions; NOCT - Nominal Operating Cell Temperature

Parameter	Value
STC Power	250 W
PTC Power	226.47 W
NOCT	45∘
SCT Power per unit area	153.7 W$/$m2
Peak efficiency	15.39 %
Number of cells	60
Maximum power current (Imp)	8.17 A
Maximum power voltage (Vmp)	30.60 V
Short circuit current (Isc)	8.71 A
Open circuit voltage (Voc)	36.3 V

Equations

$$t=1440d+m,$$

$$\overline{E}\left(d,m\right)=\left(\frac{1}{2N+1-L}\right)\left\{\sum _{n=-N}^{N}E\left(d+n,m\right)-\sum _{l\in \mathcal{L}}E\left(l,m\right)\right\},$$

$$\overline{E}\left(d,m\right)=C\left(d\right)\left\{1-{\left[\frac{1}{B\left(d\right)}\right]}^2{\left[\frac{m}{m_c}-A\left(d\right)\right]}^2\right\},$$

$$m_{0\pm }=m_c\left(A\pm B\right).$$

$$E^{\ast }\left(m^{\ast }\right)$$

$$\mathcal{E}\left(m^{\ast }\right)=-{\left(m^{\ast }\right)}^2,$$

$$R^{\ast }\left(m^{\ast }\right)=E^{\ast }\left(m^{\ast }\right)-\mathcal{E}\left(m^{\ast }\right).$$

$$R^{\ast }\left(m^{\ast }\right)$$

$$m_j^{\ast }=-1+\frac{2j}{J}\text{with }j=0,1,2,\dots J\text{and }J=\lfloor \frac{m_c⟨B⟩}{5}\rfloor ,$$

$$\Delta r^{\ast }\left(m_j^{\ast }\right)$$

$$M_r$$

$${\hat{f}}_h\left(r^{\ast }\right)=\frac{1}{hS}\sum _{s=1}^{S}K\left(\frac{r^{\ast }-r_s^{\ast }}{h}\right)$$

$$K\left(u\right)=\frac{1}{\sqrt{2\pi }}\exp \left(\frac{-u^2}{2}\right).$$

$$h$$

$$\tilde{Y}(d)$$

$$\Delta x$$

$$f\left(x\right)=\frac{1}{\nu }\exp \left\{\frac{\left(x-\mu \right)}{\nu }-\exp \left[\frac{\left(x-\mu \right)}{\nu }\right]\right\},$$

$$I(d)=\frac{1}{60}{\int }_0^{1440}E\left(d,m\right)dm,$$

$$E\left(d,m\right)=C\left(d\right)\left[1+E^{\ast }\left(m^{\ast }\right)\right],$$

$$E\left(d,m\right)=C\left(d\right)\left[1-{\left(m^{\ast }\right)}^2+R^{\ast }\left(m^{\ast }\right)\right].$$

$$E\left(d,m\right)=\left[\tilde{C}\left(d\right)+x_C\right]\left[1-{\left(m^{\ast }\right)}^2+R^{\ast }\left(m^{\ast }\right)\right]$$

$$m^{\ast }=\frac{1}{\left[\tilde{B}\left(d\right)+x_B\right]}\left[\frac{m}{m_c}-\tilde{A}\left(d\right)-x_A\right].$$

$$I(d)$$

$${\int }_{-1}^1{R}^{\ast }\left(m^{\ast }\right)d{m}^{\ast }$$

$$\mathbf{E}\left[I(d)\right]$$

$${\int }_{-1}^1\mathbf{E}\left[R^{\ast }\left(m^{\ast }\right)\right]d{m}^{\ast }$$

$$\mathbf{E}\left[R^{\ast }\left(m_j^{\ast }\right)\right]={\left\{{\int }_{-\infty }^{\infty }{R}^{\ast }g\left[R^{\ast }\right]d{R}^{\ast }\right\}}_{m_j^{\ast }},$$

$$R^{\ast }\left(m_j^{\ast }\right)=R^{\ast }\left(m_{j-1}^{\ast }\right)+\left(\frac{2}{J}\right)r^{\ast }\left(m_j^{\ast }\right).$$

Full text

For this study, solar irradiance data was collected with a weather station during a year. Post-treatment included a trimmed moving average, to smooth the data, the performance a fitting procedure using a simple model, to recover normalization parameters, and the estimation of a probability density map by means of a kernel desnsity estimation method. The normalization parameters and the probability density map allowed us to build an empiric- stochastic methodology that generates an estimate of the solar irradiance. In order to validate our method, simulated solar irradiance has been used to compute the theoretical generation of solar power, which in turn has been compared to experimental data, retreived from a commercial photovoltaic system. Since the simulation results show a good agreement has been with the experimental data, this simple methodology can estimate the solar power production and may help consumers to design and test a photovoltaic system before installation.

††journal: Solar energy

1 Introduction

Solar radiation data provide information on how much of the energy provided by the sun is irradiated on the earth’s surface, during a particular time period. This data is essential for the design, cost and energy analysis of solar power systems, and the optimal exploitation of solar energy [Wittmann et al. (2008); Martín et al. (2010)]. In particular, an improvement of the solar irradiance forecast would be very useful in predicting the performance of photovoltaic systems [Li et al. (2016); Ibrahim and Anani (2017); Bora et al. (2018)]. Nevertheless, the measurement of global horizontal radiation (GHI) data is only available in specific regions, since the initial and maintainance cost of weather stations are high. Therefore, solar energy simulation techniques stand out among the available solutions to generate GHI data on any earth’s place, due to its relative low cost [Weiss and Hays (2004); Wan et al. (2015)]. This important advantage is extended by a significant attribute of simulations, i.e. they can provide data at any discrete time step, when the appropriate model is employed. For instance, data acquired during an hour can make a significant difference in the design, evaluation and improvement of solar energy systems [Yang et al. (2012); Jahani et al. (2017); Aoun and Bouchouicha (2017)].

Another situation arises when solar irradiance data is obtained with short time-steps. Meteorological or ambient conditions, such as sky blockages, relative humidity and ambient temperature, hinder the correct measurement of the full and clear sky solar irradiance [Gueymard (2012)]. However, these adverse conditions must be taken into account to provide accurate and more realistic predictions of the electrical charge or the thermal energy generated by solar power systems [Diagne et al. (2013)].

On the other hand, several techniques had been developed for the determination of the amount of solar irradiance incident on the Earth’s surface [Khatib et al. (2012); Wan-Nik et al. (2012)], with the purpose of improving the accuracy of the predictions. Among them, time series analysis [Martín et al. (2010); Yang et al. (2012)] and artificial neural networks [Inman et al. (2013); Kaushika et al. (2014); Gutierrez-Corea et al. (2016)] are the preferred approaches, when enough historical data is available. As well, empirical models that estimate solar radiation have been developed, offering low computational costs and great accuracy [Belcher and T. (2007); Nwokolo (2017)]. Additionally, the use of hybrid models has gained popularity, since they increase the accuracy of the predictions, by taking advantage of the strengths of the differents models [Inman et al. (2013); Gordon (2009); Ji and Chee (2011)].

This article is aimed to explain the implementation of an empiric-stochastic methodology to generate solar irradiance data. From experimental measurements, the method is based on the computation of normalization prarameters, average trends and probability distributions, which are used to perform simulations of the solar irradiance along a day. In Section 2, a detailed description of the methodology is provided. In Section 3, the results from the application of the methodology, departing from experimental measurements of solar irradiance, are presented and analysed. A description of the normalization parameters, average trends and probability distributions, computed for the measured data, is given. In Section 5, performed simulations of solar irradiance are compared with experimental data, in two different ways: 1) a direct comparison concerning solar irradiance and irradiant exposure; 2) an indirect comparison refering to solar-cell electric charge which requires the conversion of the simulated irradiance into electric charge, by means of a photovoltaic model. In Section 6, some final commments about the developed methodology and its employment are discussed.

2 Methodology: Data acquisition and post-treatment

The four calendar-based seasons are traditionally bounded by the solstices and equinoxes. For instance, the 2017-2018 seasons were: Spring from March 21st, 2017 (spring equinox) to June 21st, Summer from June 22st (summer solstice) to September 21st, Autumn from September 22st (autumn equinox) to December 20st, and Winter from December 21st (winter solstice) to March 20st, 2018. Herein we have decided to consider four non-traditional seasons, which are presented and described in Table 1. This convention has been chosen, since there is more likeness between the dates around a solstice or an equinox, from an astronomical view-point, than there is for the dates going from a solstice to an equinox, or vice versa.

The campus Chetumal of the Universidad de Quintana Roo (UQRoo-Chetumal) is located at the city of Chetumal, on the southern coast of the Yucatan Peninsula, in the mexican state of Quintana Roo. The weather station of the UQRoo-Chetumal is placed on the roof of building F, “Ing. Luis Felipe Leyva”, Engineering Workshop. Its geographic coordinates are 18*∘31’25” N and 88∘*16’7”, in latitude and longitude, respectively.

The GHI was measured with a thermopile pyranometer model 8-48, collecting data from February to August 2017 under all sky conditions. Hereafter, GHI will be called simply the irradiance $E(t)$, which is a function of time $t$, and the reference time $t_{0}=0$ is taken to be midnight on (March 21st) January the 1st, 2017. For simplicity, the temporal variable is decomposed as follows:

[TABLE]

where $t$ is the time in minutes with respect to $t_{0}$, $d$ indicates the day of the year (long-term variable) and $m$ denotes the elapsed time in minutes with respect to midnight on each day (short-term variable). Therefore, the irradiance can be now defined in terms of two time variables, i.e. $E(t)\rightarrow E(d,m)$, which accounts for the two different time scales previously described. Due to this decomposition, the irradiance data of all the days of each season can be plotted and overlapped in a single graphic, as it is shown in Figure 1 for the four seasons. Additionally, notice that the irradiance data has been obtained along the daytime, with a sampling frequency of 1$/$10 $min^{-1}$. The irradiance curves, together with the results of our methodology and analysis, are discussed in Section 3.

2.1 Trimmed moving average

In order to estimate the clear sky irradiance (CSI), a trimmed moving average (TMA) has been proposed. This approach consists in averaging the symmetric data, from $N$ days before to $N$ days after the date of interest $d$, excluding the $L$ lowest values or outliers. This TMA is given by:

[TABLE]

where $2N+1-L$ is the number of days employed to perform this averaging process, and the day $l\in\left[d-n,d+n\right]$ belongs to the set $\mathcal{L}$, formed by the $L$ days that present the lowest irradiance. The application of this numerical procedure smooths out short-term fluctuations, originating from strange meteorological conditions or instrument errors-deletions, and highlights long-term irradiance trends, the diffuse and global irradiance at the local position. The TMA has been applied to the data presented in Figure 1, with $N=5$ and $L=4$, yielding the smoothed curves displayed in Figure 2. It is important to emphasize that the objective of the TMA step in our methodology is not to replace the data, but to obtain more regular curves on which the application of fitting, normalizing and probability estimation procedures should be straightforward.

2.2 Normalization variables and master curve

As it is shown in Figure 1, the daily behavior of the irradiance $E$, as a function of the short-term variable $m$, is mostly represented by a bell-like shape. Therefore, the natural choice has been to fit the data with a parabola, in order to decouple the dependence on the long-term variable $d$ and the short-term variable $m$, and to obtain a simple expression to describe the irradiance $\overline{E}$ as a function of both variables. The parabolic fitting curve reads:

[TABLE]

where $m_{c}$ is a characteristic value of the short time-scale, which we define as one half-daytime at the equator $m_{c}=360$ minutes, and the parameters $A$, $B$ and $C$ depend on the long-term variable $d$. At the given latitude, $A\left(d\right)$ indicates the time, in terms of multiples of $m_{c}$, at which noon occurs, $B\left(d\right)$ denotes the half-daytime, also in terms of multiples of $m_{c}$, and $C\left(d\right)$ is the irradiance at noon, the power per unit area in W$/$m2. The evolution of the three fitting parameters, as functions of the long-term variable $d$, are shown in Figure 3. Since these parameters are closely related to the periodic motion of the earth around the sun, we expect them to show periodic behaviors with frequency $2\pi/365$ and to be in phase with the summer solstice. These ideas, together with a detailed description of the parameters trend, are discussed in Section 3.

The times of sunrise and nightfall at the given latitude can be approximated with the following expression, obtained as the time $m=m_{0\pm}$ of zero irradiance $\overline{E}=0$ from eq. (3):

[TABLE]

The negative subindex indicates sunrise, whereas the positive subindex indicaes nightfall. The analysis in this manuscript will be focused on the irradiance data occuring within the range $m\in\left[m_{0-},m_{0+}\right]$.

Now, using the information obtained from the fitting parameters, one may procede to normalize the original data (before the application of the TMA), with the purpose of finding a general trend for the global irradiance. The following normalization variables:

[TABLE]

make the irradiance/time data collapse into a general trend, which, in turn, may be described by the irradiance estimator $\mathcal{E}$ defined by the expression:

[TABLE]

which forecasts a general behavior of the normalized irradiance. As a consequence, the normalized daytime is $m^{\ast}\in\left[m_{0-}^{\ast},m_{0+}^{\ast}\right]$, where limits are $m_{0\pm}^{\ast}=\pm 1$, respectively. The normalization, described by eq. (5), has been applied to the irradiance data and a comparison with the irradiance estimator $\mathcal{E}$, given in eq. (6), is presented in Figure 4 as function of the normalized short-term variable $m^{\ast}$, for each season defined in Table 1. The estimator $\mathcal{E}$ seems to fairly describe the average trend of the normalized experimental data $E^{\ast}$, since the data density is higher around the parabolic shape of $\mathcal{E}$. More details about the normalization effects can be found in Section 3.

2.3 Residuals, discrete probability map and probability distribution map

Despite the affinity between the estimator $\mathcal{E}$ and the average trend of the normalized experimental data $E^{\ast}$, a significant dispersion cannot be neglected. To quantify this dispersion, the residuals between the normalized experimental data $E^{\ast}$ and the irradiance estimator $\mathcal{E}$ are computed as:

[TABLE]

where $R^{\ast}$ may be considered as a random variable, which clearly separates $E^{\ast}$ from its estimator $\mathcal{E}$. Nevertheless, in the course of the present analysis, we have realized that the evolution of the residuals is better described by a sequential variation. This is due to the fact that the evolution of the solar irradiance is correlated from one instant to another, which is noticeable in our data, despite the 10 minutes lapse between consecutive acquisition moments. Therefore, the residuals $R^{\ast}\left(m^{\ast}\right)$ may not be taken as random independent variables. We should consider a second order model, for which the rate of change of the residuals $r^{\ast}=d\,R^{\ast}/d\,m^{\ast}$ becomes the independent random variable. Following this statement, we define the normalized residual at a given value of the short-term variable with the expression:

[TABLE]

As it can be observed in Figure 4, the computation of the residuals only makes sense within the normalized daytime region $m^{\ast}\in\left[-1,1\right]$. Additionally, one can notice that it is difficult to identify the exact number of points concentrated by zones, which consequently may provoke the same effect on the residuals $R^{\ast}$ and its rate of change $r^{\ast}$. Therefore, a discretization $m^{\ast}_{j}$ of the dimensionless short-term variable $m^{\ast}$ is proposed, in order to enumerate the frequency of the residuals rate of change and yield a histogram for each discrete value of $m^{\ast}_{j}$ with $j=0,1,2,\dots J$. The values of $m^{\ast}_{j}$ are spaced by $2/J=0.0313$ dimensionless units, which corresponds approximately to 10.07 minutes, whereas its first and last values should be $-1$ and $1$, respectively. In other words, we propose the discretization given by:

[TABLE]

where $\langle B\rangle=0.8990$ is the average value of the half-daytime $B$ and the number of discrete values of $m^{\ast}_{j}$ is $J+1=65$. All the histogram are computed using the same number of bins $M_{r}$, which is the average of the required number of bins for each histogram $M_{r,j}$, computed with the Freedman Diaconis’ rule [Freedman and Diaconis (1981)] and with positions in the range $r^{\ast}\in[-J/2,J/2]$, which encloses properly the values of the rate of change of the residuals. Once $M_{r}$ is known, one may proceed to determinate the bin width $\Delta r^{\ast}$. In brief, the above-mentioned methodology is summarized by the following formulas:

[TABLE]

where $\ell\left(m^{\ast}_{j}\right)$ is the number of $r^{\ast}\left(m^{\ast}_{j}\right)$ samples and $\text{IQR}\left[\cdot\right]$ indicates the interquartile range of the data, at a given value of the short-term variable $m^{\ast}_{j}$ and in the trimmed range $r^{\ast}\in[-J/2,J/2]$. The equations in (10a) describe the Freedman Diaconis’ rule, whereas equations in (10b) indicate simple averaging and bin determination formulas. For this case, we have found the value of $\Delta r^{\ast}=0.0244$ and $M_{r}=85$. Finally, each histogram must be normalized to build up a discrete probability map, whose intensity (normalized frequency) describes the probability of a rate of change for the residuals $r^{\ast}$ at each value of the dimensionless short-term variable $m^{\ast}$. This procedure for obtaining a discrete probability map for the rate of change of the residuals $r^{\ast}$ was applied on each season, as depicted in Figure 5.

Additionally, in order to estimate a probability density map (PDM) from the data, a kernel density estimation (KDE) method [Silverman (1986)] was employed. This method, applied to the rate of change of the residuals $r^{\ast}$, consists in the calculation of the following density estimator:

[TABLE]

where $r^{\ast}_{s}$ are the discrete samples, with $1\leq s\leq S$ and $S$ being the total number of samples, $h$ is the bandwidth and $K(\cdot)$ is the kernel, which has been chosen to be a Gaussian function with mean $\mu=0$ and variance $\sigma^{2}=1$ reading:

[TABLE]

The bandwidth $h$ has been set to the maximum value of $h_{j}$, according to the discretization of $m^{\ast}_{j}$, with $j=0,1,2,\dots,J$. Using a normal distribution approximation, the following formulas are applied:

[TABLE]

where $\langle h_{j}\rangle$ is the average value of $h_{j}$ and $\hat{\sigma}\left(m^{\ast}_{j}\right)$ is the standard deviation of the $r^{\ast}\left(m^{\ast}_{j}\right)$ samples. This procedure for obtaining a PDM was applied on each of the four seasons described in Table 1, as shown in Figure 6. The value of $h$ was obtained taking into account the residuals of the 365 days of measurements, thus yielding a single value $h=0.0364$ employed to compute, homogeneously, the PDMs for all the seasons.

A complete discussion on the rate of change of the residuals and a detailed description of the discrete probability and probability density maps is given in Section 3.

3 Results

3.1 Short-term solar irradiance

The GHI data acquired for each season is presented in Figure 1. In the four cases, the global irradiance data shows the superposition of a bell shape and a noise signal. The bell shape suggestes an average behavior, whereas the noise is the consequence of any sky blockage, e.g. clouds and fleeting shadows, or meteorological condition, e.g. humidity, that decrease the GHI, and temporary reflecting surfaces or clear sky conditions that increase the GHI. The sunrise, between 5:30 AM (300 minutes) and 6:30 AM (390 minutes), and the sunset, between 5:30 PM (1050 minutes) and 6:30 AM (1110 minutes), are clearly observed and vary according to the season, as a consequence of earth’s orbital movement [Vieira et al. (2012)]. The sunrise and sunset are indicated by the change from or towards, respectively, a cero value of the solar irradiance. The maxima of the GHI, which show values around 1000 W$/$m2 for Seasons 1-3 and lower values around 800 for Season 4, occur in the range of the short-term variable $m\in[600,800]$ minutes. Some minutes apart from the sunrise and sunset instants, a parabolic shape is well defined for clear-sky days, where half-daytime occurs at $m=700$ minutes. For cloudy days, the irradiance lessens to values below 250 W$/$m2. A TMA was applied to the acquired data, resulting in smoothed curves $\overline{E}(m)$ which are shown in Figure 2. This procedure removes the effect of sky blockage and yields curves that highlight the aforementioned bell shape of the GHI curves. These smoothed results allow us to compute the normalization parameters for each day, introduced in eq. (3) and discussed in the following subsection. Using these parameters, we describe the GHI average trend by means of a normalized curve $E^{\ast}(m^{\ast})$ for each day, according to its corresponding season, with the definition of the dimensionless variables given by eq. (5). In Figure 4, the normalized irradiance data is compared to the master curve, a parabola described in eq.  (6). Even though most of the data points overlap in the surroundings of the parabolic trend, a significant dispersion is in sight. A detailed analysis of this dispersion is re-engaged after the following discussion about the normalization parameters and their behavior.

3.2 Long-term solar irradiance

As introduced in the previous subsection, the normalization parameters $A(d)$, $B(d)$ and $C(d)$, described by eq. (3), were obtained for each day. The evolution of the noon time $A(d)$, the half-daytime $B(d)$ and the peak irradiance $C(d)$ is presented in Figure  3, in terms of the long-term variable $d$. Additionally, the background color in the figure indicates the corresponding season. As it has been already mentioned, the values irradiance data depends on the meteorological conditions and temporary blocking/reflecting surfaces, which provokes a variation of the normalization parameters for each day. Nevertheless, due to the periodic motion of the earth within the solar system, the normalization parameters must be periodic functions with a frequency of $2\pi/365$ cycles per day, as discussed before in this article. It is also natural to consider that the peak irradiance occur at the summer solstice (day 172), and we expect that the long-term parameters $Y(d)=\left\{A(d),B(d),C(d)\right\}$ will show periodic behaviors and will be in phase with the summer solstice. Thus, the trends $\tilde{Y}(d)=\left\{\tilde{A}(d),\tilde{B}(d),\tilde{C}(d)\right\}$ of the three fitting parameters should be accurately described by Fourier series. We have noticed that a constant term and two cosine harmonic functions are enough, reading:

[TABLE]

where the coefficients $y_{l}$ with $l=0,1,2$, which were obtained by fitting each parameter $Y(d)$ with the trend function $\tilde{Y}(d)$, are given in Table 2. The trend functions $\tilde{Y}(d)=\left\{\tilde{A}(d),\tilde{B}(d),\tilde{C}(d)\right\}$ are also shown in Figure 3, as functions of the long-term variable $d$.

The residuals $x=\left\{A-\tilde{A},B-\tilde{B},C-\tilde{C}\right\}$ between the fitting parameters and the trend functions were computed, from which histograms were constructed to quantify the dispersion from each of the trends. Employing Rice’s rule [Rice (1953)] to determinate the number of bins $M_{x}$ and their width $\Delta x$.

[TABLE]

The values of $\Delta x$ and $M_{x}$, computed from the data presented in Figure 3, are given in Table 3 for each parameter. With these values, we obtain the histograms that are also shown at the bottom of Figure 3.

Additionally, probability distribution functions (PDFs) were found for the residuals $x=\left\{A-\tilde{A},B-\tilde{B},C-\tilde{C}\right\}$, by means of fitting procedures with the appropriated distribution shapes. Exponential-based distributions were suggested to fulfill the PDFs requirements, due to their likeness with the corresponding histograms. Gumbel distribution functions [Abramowitz and Stegun (1972)] had been employed :

[TABLE]

which expected value is $\mathbf{E}[x]=\mu-\nu\gamma$, being $\gamma\approx 0.5772$ the Euler-Mascheroni constant. Their corresponding fitting parameters are presented in Table 3, whereas their comparison and good agreement with the trends shown by the histograms are observed also at the bottom of Figure 3. These non-symmetric PDFs present their highest probability at positive values of the residuals $x$, relative to the complete range of $x$, and a negative skewness, i.e. the negative-resudials tail is longer. These observations are logical, since it is more natural and usual to undergo a decrease than an increase of GHI.

3.3 Probability maps

The information given by the data presented in Figure 4 is summarized in Figure 5, for each season. The discrete probability maps represent the residuals from the simple parabolic master curve, obtained as it has been described in subsection 2.3. The overlapping data presented in Figure 4 is now represented as a normalized frequency, indicating the probability at which a variation of the residuals happens at a given time of the day and season. Close to sunrise and nightfall, the most likely range of $r^{\ast}$, the rate of change of the residuals, is narrow. In contrast, the range of $r^{\ast}$ becomes wide around noon, since the probability distribution presents long tails for both positive and negative values of $r^{ast}$. A gradual increase of the tails starts from sunrise, reaches its maximum at noon, and gradually decreases until nightfall. It is important to notice that for any season and any time during a day, a value $r^{\ast}\approx 0$ is the most likely to occur than any other value of $r^{\ast}$. Nevertheless, small values $r^{\ast}\in(-0.1,0.1)$ present a signifitnact probability to happen.

The PDMs depicted in Figure 6 correspond to continuous versions of the discrete maps presented in Figure 5. These PDMs should be employed as follows:

For a given value of the short-term variable $m$ compute its dimensionless equivalent $m^{\ast}$. One must know beforehand the values of the normalization parameters $A(d)$ and $B(d)$, according to the day number. 2. 2.

According to the season, find the horizontal coordinate $m_{j}^{\ast}$ that is closest to the calculated value of $m^{\ast}$. 3. 3.

Along a vertical line, the vertical coordinate indicates the rate of change $r^{\ast}$ of the residuals $R^{\ast}$, recalling that $2/J=0.0313$ for the presented data. 4. 4.

Finally, recalling that herein $h=0.0364$, the color indicates the probability of $r^{\ast}$ occuring at $m^{\ast}$.

The PDMs were transformed into scaled-PDMs by multiply it by the bandwidth $h$, only to find a normalized map, which values remain in the range $[0,1]$, and compare with the corresponding discrete probability maps. Nevertheless, the objective of the PDMs is not only visual, since they are required as guidelines for the generation of the residuals (random numbers) for a given value of the short-term variable $m^{\ast}$ .

4 Stochastic model description

The daily radiant exposure $I$, in W$\cdot$h$/$m2, is obtained after integration of the irradiance $E(d,m)$ with respect to the short-term variable $m$:

[TABLE]

where the coefficient 1$/$60 is required to perform a units conversion from minutes to hours.

One may proceed to employ the irradiance model for $E(d,m)$ developed in this study, which reads:

[TABLE]

recalling that the normalized short-term variable $m^{\ast}$ is related to $m$ according to Eq.(5). The normalized irradiance $E^{\ast}\left(m^{\ast}\right)$ is proposed to be decomposed into a deterministic term, the irradiance estimator $\mathcal{E}\left(m^{\ast}\right)$, and a stochastic term, the residuals $R^{\ast}$, as proposed in Eq. (7). Also, the irradiance estimator has been introduced in Eq. (6), which leads to:

[TABLE]

Additionally, according to Table 3, the parameters $A$, $B$ and $C$ can be split, as well, into the deterministic terms $\tilde{A}$, $\tilde{B}$ and $\tilde{C}$, and the stochastic terms $x_{A}$, $x_{B}$ and $x_{C}$, respectively, which yields:

[TABLE]

This equation allows the simulation of irradiance data, by means of stochastic methods, whenever the probability distributions of the random variables $x_{A}$, $x_{B}$, $x_{C}$ and $r^{\ast}\left(m^{\ast}\right)$ are known.

The daily radiant exposure $I(d)$ can be obtained by performing the integral in eq. (17), applied to the expresion in eq. (19), with the change of variables $m\rightarrow m^{\ast}$ expressed in Eq.(5), and within the domain $m^{\ast}\in\left[-1,1\right]$, since $E^{\ast}=0$ is considered outside the corresponding boundaries. Therefore, remembering that $m_{c}=360$, we find the following expression:

[TABLE]

The first integral is equal to $4/3$, whereas the second one requires a more detailed analysis, which will be given later in this section. Additionally, according to Table 3, the parameters $B$ and $C$ can be split, as well, into the deterministic terms $\tilde{B}$ and $\tilde{C}$, and the stochastic terms $x_{B}$ and $x_{C}$, respectively, which yields:

[TABLE]

When only discrete values of $m^{\ast}_{j}$ are available, in order to estimate the value of the remaining integral in Eq. (22a), one may implement a numerical method. For instance, following the trapezoidal method, one finds:

[TABLE]

according to the definition of the discrete variable $m^{\ast}_{j}$, given in Eq. (9).

Finally, an average behavior is retrieved when the expected value of $I(d)$ is computed. Assuming that all the random variables are independent, one gets:

[TABLE]

The expected value of the residual $R^{\ast}$ at a given value of the discrete short-term variable $m^{\ast}_{j}$ should be obtained directly from the integral:

[TABLE]

where $g\left[R^{\ast}\right]$ is the probability density function of $R^{\ast}$. The generation of random values of $R^{\ast}$, recalling that $R^{\ast}$ depends on the independent random variable $r^{\ast}$, requires the following procedure. From eq. (8), and applying once more the trapezoidal integration method, a discrete version in a recursive form reads:

[TABLE]

Therefore, a random value of $R^{\ast}\left(m^{\ast}_{j}\right)$ is obtained from its forerunner $R^{\ast}\left(m^{\ast}_{j-1}\right)$ and a randomly generated value of $r^{\ast}\left(m^{\ast}_{j}\right)$. This process is performed ensuring that the generated value of $R^{\ast}$ does not trespass the boundaries fixed by the criterion $f_{h}[r^{\ast}]>0.001$ for the corresponding $m^{\ast}_{j}$, which indicate the maximum and minimum allowable values of $R^{\ast}\left(m^{\ast}_{j}\right)$.

We employ Eqs. (20b) to simulate the irradiance $E$ throughout the day, whereas Eq. (4) allows us to simulate the radiant exposure $I$ for an entire day. The random values of $x_{A}$, $x_{B}$, $x_{C}$ and $r^{\ast}(m^{\ast})$ had been generated by means of the inverse transform sampling method. Examples of the irradiance $E$ (measured data, expectation and simulation) as a function of the short-term variable $m$, for four dates, are presented in Fig. 7, whereas the daily radiant exposure $I$ (measured data, expectation and simulation) as a function of the long-term variable $d$, for a year, is presented in Fig. 8.

5 Validation tests

Once our empiric-stochastic model was completed, some validation tests were made to check its correct performance. The first validation test was performed to check the short-term reliability, simulations of solar irradiance over a day, whereas the second validation test was done to verify the long-term behavior, simulations of the radiant exposure over a year.

The firts test consisted in a comparisons between the measured, the expected and the simulated irradiance profiles $E(m)$. Examples are presented in Figure 7, for four different dates of the year (days 90, 180, 270 and 360). The measured data corresponds to the data depicted in Figure 1, the expected behaviors were computed using eqs. (3) and the corresponding normalization parameters, and the simulated profiles were constructed with eq. (20b) and the random generation of the corresponding parameters. All the procedures were described in the previous section. For the cases presented for all seasons, the expected behavior describes the trends of measurements and simulations, whereas a qualitative resemblance is observed between the measured and simulated data, which is an important objective of our model.

The second test consisted in a comparison between the measured, the expected and the simulated radiant exposures $I(d)$. An example is presented in Figure 8, for a year. The measured data corresponds to the integral, defined in eq. (17) and applied to each irradiance curve $E(m)$ in Figure 1, corresponding to a different day $d$. The expected behavior was obtained from eqs. (23) and the expected values of the corresponding probability distributions. The simulated data was constructed with eq. (4) and the random generation of the residuals. All the procedures were described in the previous section. A qualitative likeness between experimental and simulated curves is observed, whereas the expected curve clearly indicates the common trend.

Despite the good qualitative agreement and the good description of the average behaviors, a further validation test was performed. The comparison between the electric charge, generated with a photovoltaic system (PV) system throughout a month, and the statistics of simulated electric charge was done. On the one hand, the PV system is installed at the UQRoo-Chetumal campus, and consists of a series-connected array of eight solar panels, which characteristics are enlisted in Table 4. On the other hand, the statistics were obtained from a hundred simulations, which were performed for each day and for the aforementioned period of time, a month. Each simulation was employed as the input of a single-diode solar-cell model [Saloux et al. (2011); Ma et al. (2014); Chaibi et al. (2018)], to generate electric current data as the output. Subsequently, the electric current along a day was integrated over time to obtain the electric charge, for each of the hundred simulations of each day. Finally, some statistics were computed, i.e. median, quantiles and full range, for each day $d$. The comparison between the measured data and the simulated statistics can be observed in figure 9, where the quantiles and the full range of the simulated data are depicted as a box and whiskers, respectively. The measured data remains within the limits of the quantile boxes for most of the days, which indicates that the individual simulations should be qualitatively equivalent to the measured data and will provide a good prediction of the electric charge that can be generated along a month.

6 Conclusions

A hybrid approach to generate global horizontal irradiantion was introduced. For a given location, solar irradiance data can be employed to compute normalization parameters, average trends and probablity distributions, in order to gather the essential ingredients of the methodology. Once these information is known, we proceed to generate random numbers, according to the probability distributions, that are added to the average trends and subjected to a de-normalization procedure, to retrieve solar irradiance data, simulated for a single daytime or a complete year.

The proposed method was tested and validated with a direct and an undirect approach, by comparison with solar irradiance data, measured by a metheorological station, and electric charge, generated by a photovoltaic solar system. A good qualitative agreement between experiments and simulations was found.

For a given location, some experimental data is required to generate the parameters that the methodology employs. Nevertheless, the proposed methodology allows the implementation of an alternative inverse engineering sequence. An educated guess of the normalization parameters, average trends and probability distributions can be done, and afterwards, compare the results against a minimum amunt of experimental data.

The presented methodology is presented as an alternative and simple procedure to provide useful data , which together with an energy load profile, can be employed as the input of a succeeding analysis, for instance a simulation software for the control and surveillance of an intelligent energy network.