Hybrid Approach of Fractional Generalized Pareto Motion and Cosine Similarity Hidden Markov Model for Solar Radiation Forecasting

Abstract

Power from solar energy is not reliable, due to weather-related factors, which diminishes the power system’s reliability. Therefore, this study suggests a way to predict the intensity of solar irradiance using various statistical algorithms and artificial intelligence. In particular, we suggest the use of a hybrid predictive model, combining statistical properties and historical data training. In order to evaluate the maximum prediction steps of solar irradiance, the maximum Lyapunov exponent was applied. Then, we used the cosine similarity algorithm in the hidden Markov model for the initial prediction. The combination of the Hurst exponent and tail parameter revealed the self-similarity and long-range dependence of the fractional generalized Pareto motion, which enabled us to consider the iterative predictive model. The initial prediction was substituted into a stochastic differential equation to achieve the final prediction, which prevents error propagation. The effectiveness of the hybrid model was demonstrated in the case study.

1. Introduction

1.1. Research Background

As the demand for green energy increases, the consumption of solar energy constantly increases. Since photovoltaic power mostly depends on solar irradiation, its power generation is not reliable. This concern can partially be addressed by the development of the corresponding predictive algorithms [1]. With the prediction of solar irradiance, an estimation of effective daily surface radiation can be achieved to instruct the configuration of the power system.

Investigation of the effective daily surface radiation shows its highly chaotic behavior [2]. The forecasting process for surface radiation is susceptible to small disturbances in the initial conditions [3]. As we know, the maximum Lyapunov exponent is often applied to measure the degree of chaos [4]. If the Lyapunov exponent is positive, then the time series exhibits chaotic behavior. Higher values of the maximum Lyapunov exponent increase the computational complexity of the forecasting process.

It was proven in a previously published work that solar irradiance exhibits long-range dependence and possesses self-similarity properties [5,6]. Long-range dependence means that the autocorrelation between different states of the signal is strong [7]. A process is said to be long-range-dependent if its autocorrelation decays to zero according to the power law. The definition of self-similarity is that the partial segment of the process resembles the entirety of the process [8]. The Hurst exponent is known as a self-similarity parameter, i.e., a measure of fractal property. The time series exhibits self-similarity if the value of the Hurst exponent ranges between 0 and 1. If the Hurst exponent of the time series ranges within (0.5, 1), then the series has long-range dependence.

Meteorological quantities are stable in a site with photovoltaic panels. Therefore, training the model from the historical effective daily surface radiation has significant importance for the forecasting of solar irradiation [9].

1.2. Literature Review

There are two data-driven predictive methods for solar irradiance: the artificial intelligence model and the statistical model [10,11]. Artificial intelligence methods require large amounts of training data [12]. Statistical predictive algorithms use statistical features from the data [13].

A neural network ensemble was used to predict solar irradiance in Nigeria [14]. In a previous study [15], a wavelet recurrent neural network was proposed to exploit the correlation between solar irradiance and weather-related factors. Elsewhere, a convolutional neural network model was constructed for prediction with automatically adjusted hyperparameters [16].

The combination of a support vector machine with signal decomposition has been used to provide the forecasting horizon for several hours [17]. In another study [18], a hidden Markov model was used for the construction of the solar irradiance predictive model. The forward–backward and Baum–Welch algorithms can also be efficiently implemented in the hidden Markov model of solar irradiance [19,20].

In comparison with the conventional Gaussian process, a multi-task Gaussian process has better adaptability and accuracy [21]. In another study [22], the probability distribution of a clearness index is proposed for the prediction of global solar radiation. Parametric probabilistic forecasting models are also constructed, based on Beta distribution and power distribution, in [23].

1.3. Research Highlights

In this study, we derive the fractional generalized Pareto motion (fGPm) with respect to the Riemann–Liouville integral and show its huge variance [24]. The fGPm is a self-similar and long-range-dependent process; thus, it can be used to predict solar radiation. In order to integrate more closely with the historical information on solar irradiance, the feature vector is constructed based on cosine similarity in the hidden Markov model (CS-HMM) [25]. The maximum prediction steps (MPS) for chaotic time series are the reciprocal of the maximum Lyapunov exponent, which is obtained from a small-data algorithm [26].

The CS-HMM algorithm does not refer to the long-range dependence and self-similarity of solar irradiance. The error in the difference iterative prediction can be propagated because the increment prediction is conducted based on the prediction value. In this paper, the CS-HMM algorithm and fGPm are combined to create a hybrid forecasting model for solar radiation. In this hybrid model, both the statistical properties and the historical data are fully employed. The initial prediction is carried out using the CS-HMM algorithm. The increment prediction in the fGPm prediction model is performed, based on the initial prediction; therefore, the error propagation is solved. We highlight the fact that the hybrid algorithm is beneficial for ensuring the reliability of the power system.

1.4. Organization of the Article

In Section 2, the definition of the fGPm model is proposed. Self-similarity and long-range dependence are proven for the fGPm model, and the difference iterative prediction model is derived. The CS-HMM forecasting model is proposed in Section 3. In Section 4, a hybrid algorithm is proposed, combining statistical properties and historical data training. The MPS in the hybrid algorithm is calculated with the maximum Lyapunov exponent. Solar radiation forecasting is performed in the case study for the purposes of validation. The work conducted in the paper id summarized in the conclusion.

2. Difference Iterative Prediction Model, Based on fGPm

2.1. Definition and Simulation of the fGPm

In several studies [27,28], a novel form of generalized Pareto distribution is introduced, which is referred to as a generalized double Pareto distribution:

$$f_D(x|\mu ,\delta ,\alpha )=\frac{1}{2\delta }{[1+\frac{|x-\mu |}{\alpha \delta }]}^{-1-\alpha },$$

(1)

where $\delta$ is the scaling parameter, $\mu$ is the location parameter, which gives information about the mean, and $\alpha$ is the tail parameter, bounded by the interval (0, 2).

Generalized Pareto motion (GPm) represents the corresponding random process. The fGPm can be constructed by using the GPm to drive the Riemann–Liouville integral.

$$P_{H,\alpha }(t)=\frac{1}{\mathrm{\Gamma }(H+0.5)}{\displaystyle \underset{0}{\overset{t}{\mathrm{\int }}}{(t-u)}^{H-0.5}dGPm}(u)$$

(2)

Note that $\frac{t^{H-0.5}}{\Gamma (H+0.5)}$ as $\mathrm{\Phi }(t)$. Thus, fGPm can be considered as the convolution of the differential of GPm and $\mathrm{\Phi }(t)$.

$$\begin{gathered} P_{H,\alpha }(t)={\displaystyle \underset{0}{\overset{t}{\int }}\frac{{(t-u)}^{H-0.5}}{\mathrm{\Gamma }(H+0.5)}}\frac{dGPm(u)}{du}du \\ =\frac{dGPm(t)}{dt}*\frac{t^{H-0.5}}{\mathrm{\Gamma }(H+0.5)} \\ =\frac{dGPm(t)}{dt}*\mathrm{\Phi }(t) \end{gathered}$$

(3)

The flow chart of the generation algorithm is depicted in Figure 1.

We performed 200 iterations of the fGPm (see Figure 2). The variance shown represents the measure of the dispersion of the distribution. A previous study [24] proved that the fGPm has a huge degree of variance, which causes the process to be heavy-tailed. The heavy-tailed characteristics of the fGPm indicate that the simulation path contains frequent long jumps.

A three-dimensional fGPm random walk with 200 iterations is depicted in Figure 3. The particle is moving in the three-dimensional space, with the initial position as its origin. In each iteration of the random walk, the moving direction for the particle is random and the jump length is subject to fGpm. For the heavy-tailed stochastic process, the jump lengths are usually long, and the area of movement is vast.

2.2. Long-Range Dependence and Self-Similarity of the fGPm

The parameters in the generalized Pareto distribution have a similar meaning as the Lévy distribution [29]. When the Hurst exponent, H, of the fractional Lévy stable motion is greater than $\frac{1}{\alpha }$, the previous values and successive values are positively related, which ensures the prediction accuracy. Therefore, the fractional Lévy stable motion has long-range dependence for $H\in (\frac{1}{\alpha },1)$. Similarly, the long-range dependence of the fGPm also requires that $H\in (\frac{1}{\alpha },1)$.

In the stochastic process, $X(t)$ is self-similar if a constant $k$ exists to make the following condition holds true for any constant $c$:

$$X(ct)\triangleq c^{k}X(t),$$

(4)

where $\triangleq$ denotes equality in the distribution.

Similar to the Lévy stable motion, the GPm is ${\alpha }^{-1}$ and is self-similar:

$$GPm(ct)\triangleq c^{\frac{1}{\alpha }}GPm(t).$$

(5)

Note that the fGPm is a self-similar process with the parameter of $H-\frac{1}{2}+\frac{1}{\alpha }$, and the derivation is as follows:

$$\begin{gathered} P_{H,\alpha }(ct)=\frac{1}{\mathrm{\Gamma }(H+0.5)}{\displaystyle \underset{0}{\overset{ct}{\int }}{(ct-u)}^{H-0.5}dGPm(u)} \\ =\frac{c^{H-0.5}}{\mathrm{\Gamma }(H+0.5)}{\displaystyle \underset{0}{\overset{ct}{\int }}{(t-\frac{u}{c})}^{H-0.5}dGPm(u)} \\ =\frac{c^{H-0.5}}{\mathrm{\Gamma }(H+0.5)}{\displaystyle \underset{0}{\overset{ct}{\int }}{(t-\frac{u}{c})}^{H-0.5}d(GPm(\frac{u}{c})}{c}^{\frac{1}{\alpha }}) \\ \stackrel{\frac{u}{c}\to u}{=}\frac{c^{H-\frac{1}{2}+\frac{1}{\alpha }}}{\Gamma (H+0.5)}{\displaystyle \underset{0}{\overset{t}{\int }}{(t-u)}^{H-0.5}dGPm(u)} \\ =c^{H-\frac{1}{2}+\frac{1}{\alpha }}{P}_{H,\alpha }(t) \end{gathered}$$

(6)

2.3. The Difference and Iteration in the fGPm Predictive Model

The difference iterative prediction model is based on the Langevin-type stochastic differential equation (SDE), driven by the fGPm [30]:

$$d{X}_{H,\alpha }(t)=\eta X_{H,\alpha }(t)d(t)+\sigma X_{H,\alpha }(t)dfGPm(t),$$

(7)

where $\eta$ and $\sigma$ denote the global trend and volatility rate, respectively.

Note that the increment of the fGPm is not analytical; therefore, we apply the Maruyama notation for the derivation of the predictive model [31]:

$$dfGPm(t){=\mathrm{w}}_{H,\alpha }\left(\mathrm{t}\right){(dt)}^{H-\frac{1}{2}+\frac{1}{\alpha }}$$

(8)

where $w_{H,\alpha }(t)$ is the white noise conforming to a generalized double Pareto distribution [24]. The location parameter is set to be zero for further derivation and the other parameters are estimated from the actual incremental data.

This provides the facility to combine the discretized forms of Equations (7) and (8):

$$\Delta X_{H,\alpha }(t)=\eta X_{H,\alpha }(t)\Delta t+\sigma X_{H,\alpha }(t)w_{H,\alpha }(t){(\Delta t)}^{H-\frac{1}{2}+\frac{1}{\alpha }},$$

(9)

Therefore, we can conclude the prediction formula for the future values:

$$\begin{gathered} X_{H,\alpha }(t+\Delta t)=X_{H,\alpha }(t)+\Delta X_{H,\alpha }(t) \\ =X_{H,\alpha }(t)+\eta X_{H,\alpha }(t)\Delta t+\sigma X_{H,\alpha }(t)w_{H,\alpha }(t){(\Delta t)}^{H-\frac{1}{2}+\frac{1}{\alpha }}. \end{gathered}$$

(10)

Let $\Delta t$ be a predictive time frame:

$$\begin{gathered} X_{H,\alpha }(i+1)=X_{H,\alpha }(i)+\Delta X_{H,\alpha }(i) \\ =X_{H,\alpha }(i)+\eta X_{H,\alpha }(i)+\sigma X_{H,\alpha }(i)w_{H,\alpha }(i), \end{gathered}$$

(11)

where $i$ is the prediction time point.

This concludes the derivation for the difference iterative prediction model, based on the fGPm. The forecasting process is an additive iteration of the differences between two adjacent time points.

2.4. Parameter Estimation

In order to estimate the parameters for the predictive model, the increments of the dataset are calculated.

The estimation of location parameters is the mean of the incremental data:

$$\widehat{\mu }=\frac{1}{n}{\displaystyle \sum _{i=1}^n\Delta X_{H,\alpha }(i)}.$$

(12)

In order to estimate the scaling and tail parameters, we use the maximum likelihood estimation, which reads as follows:

$$\begin{gathered} L(x|\widehat{\mu },\delta ,\alpha )={\displaystyle \prod _{i=1}^n\frac{1}{2\delta }{[1+\frac{|\Delta x_i-\widehat{\mu }|}{\alpha \delta }]}^{-1-\alpha }} \\ =\frac{1}{2^n}{\alpha }^n{(\alpha \delta )}^{n\alpha }{\displaystyle \prod _{i=1}^n\frac{1}{{(\alpha \delta +|-\widehat{\mu }+\Delta x_i|)}^{1+\alpha }}}. \end{gathered}$$

(13)

Then we take the logarithm to both sides of Equation (13), and we get the following:

$$\ln L(x|\widehat{\mu },\delta ,\alpha )=-n\ln 2+n\ln \alpha +n\alpha \ln \alpha +n\alpha \ln \delta -(1+\alpha ){\displaystyle \sum _{i=1}^n\ln (\alpha \delta +|-\widehat{\mu }+\Delta x_i}|).$$

(14)

Furthermore, the partial derivatives of Equation (14) are taken with respect to the scaling parameter and tail parameter. After elementary transformations, we obtain the following identities:

$$\frac{\partial \ln L(x|\widehat{\mu },\delta ,\alpha )}{\partial \delta }=\frac{n\alpha }{\delta }-\alpha (1+\alpha ){\displaystyle \sum _{i=1}^n\frac{1}{\alpha \delta +|-\widehat{\mu }+\Delta x_i|}}=0$$

(15)

$$\frac{\partial \ln L(x|\widehat{\mu },\delta ,\alpha )}{\partial \alpha }=\frac{n}{\alpha }+n(\ln \alpha \delta +1)-\delta (1+\alpha ){\displaystyle \sum _{i=1}^n\frac{1}{\alpha \delta +|-\widehat{\mu }+x_i|}}-{\displaystyle \sum _{i=1}^n\ln (}\alpha \delta +|-\widehat{\mu }+\Delta x_i|)=0.$$

(16)

To estimate parameter $\eta$, we define a new quantity, $X_{\eta }$.

$$X_{\eta }(i)=\frac{\Delta X_{H,\alpha }(i)}{X_{H,\alpha }(i)}=\eta +\sigma w_{H,\alpha }(i)$$

(17)

The mean of the white noise is zero, therefore Equation (18) can be established.

$$\begin{gathered} E(X_{\eta })=E(\eta +\sigma w_{H,\alpha }) \\ =E(\eta )+E(\sigma w_{H,\alpha }) \\ =E(\eta ). \end{gathered}$$

(18)

Therefore, the estimation of $\eta$ reads as follows:

$$\begin{gathered} \widehat{\eta }=E(\eta )=E(X_{\eta }) \\ =\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_{\eta }}(i)=\frac{1}{n}{\displaystyle \sum _{i=1}^n\frac{\Delta X_{H,\alpha }(i)}{X_{H,\alpha }(i)}} \end{gathered}$$

(19)

In order to calculate parameter $\sigma$, we introduce a new variable, $X_{\sigma }$, which reads as follows (see Equation (20)):

$$X_{\sigma }(i)=X_{\eta }(i)-\widehat{\eta }=\sigma w_{H,\alpha }(i).$$

(20)

Then, one can estimate the value of $\sigma$:

$$\begin{gathered} \widehat{\sigma }=E(\sigma )=E(\frac{X_{\sigma }}{w_{H,\alpha }}) \\ =\frac{1}{n}{\displaystyle \sum _{i=1}^n\frac{X_{\sigma }(i)}{w_{H,\alpha }(i)}}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\frac{X_{\eta }(i)-\widehat{\eta }}{w_{H,\alpha }(i)}} \end{gathered}$$

(21)

3. Combination of the Hidden Markov Model and Cosine-Similarity

The core of the hidden Markov model lies in the directed graph. Any state variable that is not observed directly is called a hidden variable. As time passes, the system switches among different states. Each state variable has an observed variable. In the prediction of meteorological quantities, extracting information from historical data is important because the climate in the area is stable. Therefore, cosine similarity is employed for similarity analysis of the solar irradiance within the same time frames (one year). Our CS-HMM model represents a hidden Markov model, in which the feature vector is calculated by using a cosine similarity measure.

Let V and W be n-dimensional vectors; one can visualize the similarity of the two vectors as the cosine of the angle between V and W. From a geometric standpoint, the angle between V and W is given by:

$$COS[V,W]=\frac{VW}{\left|\right|V\left|\right|\left|\right|W\left|\right|}=\frac{{\displaystyle \sum _{i=1}^n{v}_i{w}_i}}{({\displaystyle \sum _{i=1}^n{v}_i{}^2}{\displaystyle \sum _{i=1}^n{w}_i{}^2{)}^{\frac{1}{2}}}},$$

(22)

where $v_i$ and $w_i$ are the elements of vectors $V$ and $W$.

In order to calculate the feature vector, we need to constitute a radiation vector. The values in the vector are the solar radiation for the same day on five consecutive years. Then, the variation of the solar irradiance between two adjacent values in the radiation vector is calculated to form the variation vector. The three-dimensional vectors, named $V_1$ and $V_2$, are defined as follows:

$$V_1=[v_1,v_2,v_3],$$

(23)

$$V_2=[v_2,v_3,v_4],$$

(24)

where $v_1,v_2,v_3,v_4$ are the elements of the variation vector.

The cosine similarity between $V_1$, $V_2$, and a common criterion vector, $W$, are calculated.

The elements in the three-dimensional criterion vector are randomly chosen as 0 and 1.

$$c_1=COS[V_1,W]$$

(25)

$$c_2=COS[V_2,W]$$

(26)

Equation (27) defines the feature vector, $F$:

$$F=[c_1,c_2]$$

(27)

For the predictive analysis, we select a value for the training set to be the candidate for prediction. Then, we calculate the feature vector between the candidate and the four last values of the training set. The vector with the maximum probability is preferred, and the corresponding candidate is the predicted value. If we set a prediction length, n, and repeat the process for n times, then the solar radiation prediction, based on CS-HMM, is completed.

The nonhomogeneous Markov jump system is suitable for describing the dynamical systems when subject to random changes [32]. In solar irradiation prediction, the disturbance in the system is subject to the historical data and meteorological conditions. Therefore, it is more appropriate to employ a cosine similarity to extract historical information. The Takagi–Sugeno fuzzy model and singular system theory are not applied in this research, due to the excess computational complexity involved [33,34].

4. Hybrid Algorithm, Based on fGPm and CS-HMM

4.1. Calculation for the Number of Prediction Steps

To begin the hybrid prediction, the MPS must be determined. Within the MPS, the accuracy of the prediction can be secured. The MPS is the reciprocal of the maximum Lyapunov exponent, which is obtained via the small-data method [26]:

$$L=\frac{1}{\max Ly},$$

(28)

where max Ly is the maximum Lyapunov exponent.

4.2. Hybrid Predictive Algorithm

In this paper, the predictions based on statistical properties and artificial intelligence are combined to form a hybrid prediction model. Both the statistical characteristics and the training of historical data are utilized. The hybrid algorithm, which is based on fGPm and CS-HMM, can be divided into six steps:

(1)

The MPS is obtained by calculating the maximum Lyapunov exponent.

(2)

The long-range dependence condition for the fGPm difference iterative prediction model is tested.

(3)

The initial prediction of the solar radiation is carried out, based on the CS-HMM algorithm.

(4)

The initial prediction results are fed to the SDE, driven by fGPm to predict future increments.

(5)

The Maruyama notation is employed in the discretized SDE.

(6)

The prediction is finalized with the difference and iteration.

The flow chart of the prediction is in Figure 4. The essence of the hybrid algorithm isa two-step prediction. First, the CS-HMM algorithm is utilized to reach the initial prediction results. Second, the initial prediction results are substituted into the SDE for the prediction of increments. The final prediction results are the additive iteration of present values and increments.

5. Case Study

5.1. Dataset of Solar Irradiation

In order to validate our investigation, solar radiation forecasting is performed in the case study. The solar radiation is collected by the meteorological facility in Beijing, which can be found in the Supplementary Materials. The data for this is collected once a day. The radiation data from 1 July 2014 to 8 October 2014 are used as the training data. The training data are depicted in Figure 5, in which we can visualize the fluctuation of solar radiation.

As we see from Figure 5, the fluctuation of the solar irradiance is stationary and periodic. The augmented Dickey–Fuller test is performed using MATLAB, which confirms the stationary nature of the dataset [35].

5.2. The Number of Prediction Steps for Solar Irradiation Prediction

The logarithmic regression is plotted in Figure 6 for the shortest distance in the reconstructed phase space. As we can see from Figure 6, the slope of the straight line is the maximum Lyapunov exponent, which is 0.0199. First, we take the integral part of the reciprocal, then we can set the MPS to be 50 time steps.

5.3. Self-Similarity and the Long-Range Dependence of Solar Radiation Prediction

We should recall that the condition for long-range dependence is the following: $\alpha H>1$. We must validate the estimated Hurst and tail parameters to ensure that they satisfy the criterion.

The Hurst exponent was calculated as the slope of the straight line (see Figure 7). An estimation of the Hurst exponent, evaluated by the R/S method, gives 0.7198. The value of the Hurst exponent verifies the long-range dependence and self-similarity of the data.

The calculated tail parameter, $\alpha$, takes the value of 1.9708. The product between the Hurst exponent and the tail parameter is greater than 1, which shows that the model has long-range dependence.

5.4. Evaluation of the Hybrid Predictive Model

Prediction of solar irradiance was carried out with the hybrid algorithm, the fGPm difference iterative prediction model, and the CS-HMM algorithm. The prediction results are summarized in Figure 8. The evaluation criteria were calculated for three predictive models, as shown in

Table 1. Evaluation criteria for the predictive methods.

Predictive Method	MAE	RMSE	Relative Error
			Max (%)	Mean (%)
Hybrid	0.8984	1.9864	2.97	1.13
fGPm	0.9265	2.0752	3.06	1.22
CS-HMM	0.9574	2.1224	3.14	1.31

. The mean absolute error (MAE) and the root mean squared error (RMSE) are given by the following formulas:

$$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n|s_i}-k_i|,$$

(29)

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{(s_i-k_i)}^2}},$$

(30)

where $s_i$ is the ith value of the prediction and $k_i$ is the ith value of the actual data.

5.5. Discussion and Future Work

It can be inferred from Table 1 that the hybrid predictive algorithm shows better performance than the other two methods. The CS-HMM does not take into account the long-range dependence and self-similarity of the data. The disadvantage of the fGPm predictive method is as follows: the current increment is predicted on the basis of the present prediction. The prediction error can be propagated and enlarged during the prediction process. In the hybrid algorithm, long-range dependence and the self-similarity of solar radiation are considered. The error propagation is stopped by utilizing the initial prediction of CS-HMM to form the increment prediction.

Our further investigations will be aimed at an extension of the predictive model. The power load series is chaotic; therefore, the MPS needs to be calculated [36]. The power load series is long-range-dependent and self-similar [37,38]. As reported in an earlier study [39], the historical data can be used for the prediction of power load. Therefore, the hybrid predictive model can be extended to power load prediction.

6. Conclusions

In this study, we have suggested a hybrid predictive method that is based on the fGPm and CS-HMM models. Self-similarity and the long-range dependence of the fGPm are proven. The CS-HMM employs cosine similarity to extract information from the historical variation. In the hybrid algorithm, the MPS is determined with the maximum Lyapunov exponent. The initial prediction is performed by the CS-HMM and is then substituted into the fGPm difference iterative model to predict the successive values. The hybrid algorithm was then verified, using the dataset of solar irradiance.

The hybrid method utilizes the statistical properties of the dataset and avoids error propagation. The method’s predictive accuracy can be improved if the initial predictive algorithm is optimized.