An Integrated Framework Based on an Improved Gaussian Process Regression and Decomposition Technique for Hourly Solar Radiation Forecasting

Abstract

The precise forecast of solar radiation is exceptionally imperative for the steady operation and logical administration of a photovoltaic control plant. This study proposes a hybrid framework (CBP) based on complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN), an enhanced Gaussian process regression with a newly designed physical-based combined kernel function (PGPR), and the backtracking search optimization algorithm (BSA) for solar radiation forecasting. In the CEEMDAN-BSA-PGPR (CBP) model, (1) the CEEMDAN is executed to divide the raw solar radiation into a few sub-modes; (2) PACF (partial autocorrelation coefficient function) is carried out to pick the appropriate input variables; (3) PGPR is constructed to predict each subcomponent, respectively, with hyperparameters optimized by BSA; (4) the final forecasting result is produced by combining the forecasted sub-modes. Four hourly solar radiation datasets of Australia are introduced for comprehensive analysis and several models available in the literature are established for multi-step ahead prediction to demonstrate the superiority of the CBP model. Comprehensive comparisons with the other nine models reveal the efficacy of the CBP model and the superb impact of CEEMDAN blended with the BSA, respectively. The CBP model can produce more precise results compared with the involved models for all cases using different datasets and prediction horizons. Moreover, the CBP model is less complicated to set up and affords extra decision-making information regarding forecasting uncertainty.

1. Introduction

1.1. Background and Motivation

The current global energy crisis has added new urgency to speed up clean energy (wind and solar PV) transitions and, once again, emphasized the vital role of renewable energy. Notably, clean energy power generation, such as wind and solar PV, has great potential in reducing the dependence of the power sector on many fossil fuels including coal, oil, and natural gas. As stated in the latest Renewables Global Status Report in 2021, the new global renewable power capacity has broken a new record, reaching nearly 315 GW. Together, solar PV and wind sources accounted for nearly 90% of all new renewable power additions. Solar PV represented more than half of the additions (around 175 GW) and wind power another 102 GW [1]. Solar electricity is becoming one of the most appealing energy-saving sources due to its tremendous potential. However, solar energy exhibits burdensome characteristics such as intermittency, fluctuation, and stochasticity due to the effects of global climate change and human intervention [2]. Therefore, integrating PV energy sources into a power grid becomes more difficult [3]. As a result, accurate and dependable forecasting information about solar radiation is extremely important for choosing the most cost-effective power grid strategy and power storage capacity [4].

1.2. Literature Review

Recently, a variety of successful applications have been applied to the forecasting of solar energy. These methodologies can be broadly classified into the three most important approaches: physical models (PMs), traditional time series (TS) models, and machine learning (ML) models [2]. The PMs rely heavily on Numerical Weather Prediction (NWP). Based on high-performance computers and the actual situation of the atmosphere, NWP solves a number of thermo-hydrodynamic equations with weather evolution to predict the future solar radiation intensity. In its calculation, multiple data sources are required, such as climate data and satellite cloud observations. Therefore, the prediction capability is determined by the consistency of weather conditions [5]. Thereby, these processes are intricate and time-consuming. Traditional TS models are primarily based on time series and mainly rely on historical measured data to predict solar radiation value in future, such as the autoregressive (AR) model [6], and variants of the AR model such as ARMA [7,8] and ARIMA [9], as well as the Grey model. In early solar radiation forecasting, the classical TS models were commonly utilized; however, their performance is relatively poor due to the fixed parameters and inability to adjust to shifting time sequences. Many researchers made efforts to strengthen the model performance by introducing other techniques. Ding et al. [10] proposed an improved Grey model with real-time updated parameters to perform photovoltaic power forecasting. Unfortunately, the Grey model cannot deal with a massive amount of stochastic information effectively, which leads to the degradation of prediction performance.

In contrast, ML models are more commonly used today, thanks to the development of high-performance computers, big data mining techniques, and artificial intelligence theory. Most widely used are artificial neural networks (ANN) [11], random forests (RF), support vector machines (SVM) [12], extreme learning machines (ELM) [13], and heuristic intelligent optimization algorithms based on the above-mentioned methods [14]. Yagli et al. [15] investigated the forecasting abilities of 68 off-the-shelf ML methods for solar irradiance forecasting in seven locations that had different meteorological characteristics. According to their findings, there is no one-size-fits-all model, although there are some preferred options for each sky and climate. Based on ANN, Moreira et al. [16] created an ANN ensemble to forecast medium-term photovoltaic power. They provided a versatile method that enhanced the forecasting results by allowing for the modification of various experimental elements, such as the base forecasting module and the desired forecast horizon. They also suggested using the feature selection method to eliminate redundant information. Sun et al. [17] proposed a progressive back-propagation neural network (BP) model for PV module temperature. The step-by-step BP model produced more precise prediction results relative to the direct one. Benali et al. [18] investigated the performance of two famous ML methods, ANN and RF, in hourly solar radiation forecasting. Experimental analysis suggested that the RF technique outperformed ANN in forecasting solar data. Cervone et al. [19] presented a new model by integrating ANN and Analog Ensemble (AnEn) together. The results indicated that the combined model outperforms each method running individually.

More recently, the hybrid model with decomposition techniques and/or optimization algorithms has been shown to exhibit a more desirable forecasting performance and its use has grown rapidly. Zhang et al. [14] provided a comprehensive comparison of the single BP, BP with optimum parameters searched by particle swarm optimization (PSO), and several widely used statistical models for daily solar radiation prediction. They found that PSO-BP was more effective than all involved benchmark models, confirming that the compound model surpasses the standalone model. Feng et al. [20] developed a new hybrid model in which PSO is used to find the suitable parameters of ELM. The PSO-ELM has better forecasting ability than the five other ML models in forecasting daily solar radiation. Undoubtedly, the hybrid models integrated with several complementary techniques are more promising.

From the aspect of real operation, the uncertainty of forecasting results is necessary to include as additional information for decision-making. Unfortunately, most existing research seems to ignore that factor. A novel ML approach, namely GPR, has the advantages of traditional ML in dealing with non-linear problems, but also inherits the flexible inductive reasoning ability of Bayesian methods. The capability of nonlinear fitting and uncertainty quantification of GPR has been successfully tested in a variety of applications, including wind energy forecasting [21] and monthly streamflow forecasting [22]. Additionally, some researchers made efforts to enhance the ability of GPR by utilizing combined kernel functions [23] or by bringing in emerging optimization algorithms such as the PSO algorithm and Bayesian Optimization [24]. To the best of our knowledge, the advantage of GPR has not been fully investigated in the solar energy domain.

Due to the influence of environmental and meteorological variables such as irradiance and cloudiness, observed solar radiation series show some obvious features such as strong stochasticity, complex nonlinearity, and being non-stationary. The historical sequence of solar radiation is formed by many components, including tendency, periodicity, and randomness. To improve the precision of model predictions, it is vital to employ suitable series decomposition techniques to dig the concealed information in the observed series. The widely used decomposition approaches contain variational mode decomposition (VMD) [25], discrete wavelet transform (DWT) [26], empirical mode decomposition (EMD) [27], and their variants including empirical wavelet transform (EWT) [28], ensemble EMD (EEMD) [29], and complementary EEMD (CEEMD) [5]. For example, Zhang and Wei [30] utilized the WT to process the original daily solar radiation, and then used principal component analysis (PCA) to decrease input dimension; after that, ELM optimized by bat algorithm (BA) was established to make predictions. To address the mode mixing of EMD, the high computational cost of EEMD, and to make the reconstructed time series noise-free, another improved version of EMD was proposed [31]. The name of the new version is complete EEMD with adaptive noise (CEEMDAN), and it is appropriate for decomposing a non-stationary signal into several more stationary sub-modes. Therefore, this study adopts the CEEMDAN to disassemble the initial solar series into several simpler sub-series to improve forecasting accuracy.

1.3. Contribution and Organization

Through the above systematic analysis of the existing studies, some deficiencies and problems are found. The standard GPR used a single kernel function, which is incompetent for accurate real solar radiation series forecasting due to the fact the periodicity, tendency, and randomness features of the solar radiation are intertwined with each other. When utilizing the forecasting model in real operation, the uncertainty of forecasting results is not provided.

Therefore, to provide more accurate forecasting results, a compound framework, namely CEEMDAN-BSA-PGPR (CBP), is put forward by integrating the CEEMDAN, BSA, and GPR models with a newly designed physical-based kernel function. The CBP is established to adequately perform from one-hour to three-hour solar radiation forecasting. Contributions of this study mainly include the following:

(1)

To reduce the volatility of the original solar data, the non-stationary solar radiation series is divided into several stable sub-series with various scales via the CEEMDAN method. These stable components are simpler than the original series for forecasting, which is beneficial to enhance forecasting accuracy.

(2)

According to the inherent features, including periodicity, tendency, and randomness, of the original solar radiation series, a new physical information-based kernel function is designed by using several commonly used kernel functions. A GPR model based on the new physical kernel (PGPR) is constructed. The model PGPR can realize both deterministic and directly probabilistic solar radiation prediction, which is more in line with the actual demand during decision-making.

(3)

To handle local optimality and initialization dependence disadvantages of the conjugate gradient (CG) algorithm, widely used for searching the optimal hyper-parameters, BSA is executed to search for the optimum hyperparameters of PGPR.

The following is the layout of the remainder of the paper: Section 2 presents the theoretical foundation of this study, including CEEMDAN, PGPR, BSA, and the CBP; a comprehensive case study is given in Section 3; Experimental results of all involved models using four datasets and for all horizons are presented in Section 4; Section 5 analyzes and testifies to the positive effects of the CEEMDAN and the BSA algorithm; Section 6 summarizes the main findings of this paper. The nomenclature for abbreviations is shown in Appendix A.

2. Methodology

2.1. Solar Radiation Prediction Based on ML

Solar radiation forecasting based on ML can be stated as:

$$s\left(t+L\right)=\phi \left(s\left(t-d_1+1\right),Other\left(t-d_2+1\right)\right)$$

(1)

where $s\left(t+L\right)$ represents the computed solar radiation value at t + L in the future, L is the corresponding forecasting horizon, $s\left(t-d_1+1\right)$ stands for the measured solar radiation value up to $t-d_1+1$ time steps, $Other\left(t-d_2+1\right)$ is the other variable having great relevance with the solar radiation at time t + L, such as precipitation, surface pressure, temperature, etc., and $\phi \left(·\right)$ is the forecasting function which can be replaced by one of the three types of methods mentioned in the introduction.

PMs can sometimes be replaced by ML models because of issues such as a lack of modeling data, a poor understanding of physical mechanisms, or an inaccurate forecast of long-term weather conditions. In addition, ML models are simple to construct and can be easily expanded by utilizing additional optimization algorithms or series decomposition methods to create new hybrid models, resulting in highly accurate forecasting.

2.2. Complete Ensemble Empirical Mode Decomposition with Adaptive Noise

The CEEMDAN [31] retains the advantages of the EEMD method to effectively eliminate the mode-mixing hassle of EMD via repeatedly feeding white noise to the preliminary series. Simultaneously, it can likewise use less computation sources and improve the decomposition efficiency. The primary steps of CEEMDAN are as follows:

(1)

Create a series $S\left(t\right)+{\epsilon }_0{n}^i\left(t\right)$ by adding additional white noise $n^i\left(t\right)$ to the original series $S\left(t\right)$, where ${\epsilon }_0$ is the noise amplitude;

(2)

Execute EMD N times to generate the first components $IM{F}_1^i\left(t\right)\left(i=1,2,\cdots ,N\right)$. The final first component $IM{F}_1^{}\left(t\right)$ is calculated by taking the mean of $IM{F}_1^i\left(t\right)$. The corresponding first residual series $r_1^{}\left(t\right)$ is obtained by $r_1\left(t\right)=S\left(t\right)-IM{F}_1^{}\left(t\right)$;

(3)

Decompose the set of noise-added residual series $r_1\left(t\right)+{\epsilon }_1{E}_1\left(n^i\left(t\right)\right)$ to produce the second $IM{F}_2^i\left(t\right)\left(i=1,2,\cdots ,N\right)$; $E_1\left(n^i\left(t\right)\right)$ is the first sub-mode of the $n^i\left(t\right)$. Thus, the final second $IM{F}_2^{}\left(t\right)$ is calculated by taking the average of $IM{F}_2^i\left(t\right)$;

(4)

The subsequent processes of decomposition are identical to Steps 2–3. After calculating the k-th residual series, the (k + 1)-th mode can be obtained. Continue doing so until the residual series exhibits monotony. The remaining modes of CEEMDAN can be represented mathematically as:

$$r_k\left(t\right)=r_{k-1}\left(t\right)-IM{F}_k\left(t\right)=r_{k-1}\left(t\right)-1/N{\displaystyle {\sum }_{i=1}^{N}IM{F}_k^i\left(t\right)}$$

(2)

The following can be used to compute the final residual R(t):

$$R\left(t\right)=S\left(t\right)-{\displaystyle {\sum }_{k=1}^{K}IM{F}_k\left(t\right)}$$

(3)

Therefore, the original solar data can be written as:

$$S\left(t\right)=R\left(t\right)+{\displaystyle {\sum }_{k=1}^{K}IM{F}_k\left(t\right)}$$

(4)

2.3. Gaussian Process Regression with Physical-Based Kernel Function (PGPR)

2.3.1. Basic Principles of Gaussian Process Regression

The GPR, a famous ML method, is a non-parametric Bayesian method [32]. It was proposed by integrating Bayesian and statistical theory. Therefore, the GPR possesses both advantages of the Bayesian method and ML. These characteristics grant the GPR model superiority in handling high-dimensional complicated nonlinear issues; therefore it is widely used in all walks of life [32].

The training input sample of GPR can be denoted as $S_n=\left\{\left({\mathit{X}}_i,y_i\right)\right\},i=1,\cdots ,n,\hspace{1em}\mathit{X}\in R^m, y\in R$, where X is the input data with n observations for the corresponding target output y with m-dimensional factors. The stochastic process state set of the X obeys a joint Gaussian distribution with n dimensions. According to the definition, the state set g of the stochastic process is a Gaussian process. Thereby, its unique probability density function, denoted GP, can be obtained according to its mean function $E\left(X\right)$ and covariance function matrix $K\left(X,X\right)$.

$$g\left(\mathit{X}\right)\sim GP\left(E\left(\mathit{X}\right),K\left(\mathit{X},\mathit{X}\right)\right)$$

(5)

Then, in accordance with the features of GP, the y of the training dataset $S_n$ and the $y_{test}$ of the testing dataset $S_n^{test}=\left\{\left(x_{test},y_{test}\right)|x_{test}\in R^m,y_{test}\in R\right\}$ obey the following multivariate Gaussian distribution:

$$\left[\begin{array}{c}\mathit{y}\\ y_{test}\end{array}\right]\sim N\left(0,\left[\begin{array}{cc}K(\mathit{X},\mathit{X})+{\sigma }_n^2{I}_n& K(\mathit{X},x_{test})\\ K(x_{test},\mathit{X})& k(x_{test},x_{test})\end{array}\right]\right)$$

(6)

where $K(x_{test},\mathit{X})$ equals $K{(\mathit{X},x_{test})}^T$, and $k(x_{test},x_{test})$ is the covariance for $x_{test}$ itself.

After that, according to the Bayesian principle, the posterior distribution can be inferred as:

$$\begin{array}{l}p(y_{test}|\mathit{X},\mathit{y},x_{test})\sim N(E\left(y_{test}\right),\mathrm{cov}(y_{test}))\\ E\left(y_{test}\right)=K(x_{test},\mathit{X}){\left[K\left(\mathit{X},\mathit{X}\right)+{\sigma }_n^2{I}_n\right]}^{-1}y\\ \mathrm{cov}(y_{test})=k(x_{test},x_{test})-K(x_{test},\mathit{X}){\left(K(\mathit{X},\mathit{X})+{\sigma }_n^2{I}_n\right)}^{-1}K(\mathit{X},x_{test})\end{array}$$

(7)

where $E\left(y_{test}\right)$ is the expected value, and $\mathrm{cov}(y_{test})$ represents the posterior variance for $y_{test}$. The $\mathrm{cov}(y_{test})$ is used to measure the forecasting uncertainty.

2.3.2. Covariance Functions and Their Hyperparameters Determination

As stated above, GP can be obtained by using its $E\left(X\right)$ and $K\left(X,X\right)$. Based on the statistical principles, the GP can be easily transformed into a new Gaussian distribution with $E\left(X\right)=0$. As a result, figuring out the GPR model’s $K\left(X,X\right)$ is the most important part of solving it. The $K\left(X,X\right)$ of GPR are numerous. Three commonly used single $K\left(X,X\right)$ are presented:

(1)

Squared-Exponential (SE) kernel:

$$k_{SE}={\sigma }_f^2\exp \left(-r/2l^2\right)$$

(8)

(2)

Matern (Ma) 5/2 kernel:

$$k_{Ma}={\sigma }_f^2\left(1+\sqrt{5}\frac{r}{l}+\frac{5r^2}{3l^2}\right)\exp \left(-\sqrt{5}\frac{r}{l}\right)$$

(9)

(3)

Rational Quadratic (RQ) kernel:

$$k_{RQ}={\sigma }_f^2{\left(1+r^2/2\alpha l^2\right)}^{-\alpha }$$

(10)

(4)

Linear kernel:

$$k_{Lin}={\sigma }_f^2\left(x-c\right)\left(x’-c\right)$$

(11)

(5)

Random noise (RN) kernel:

$$k_{RN}={\sigma }_f^2{\delta }_{x{x}^{’}}$$

(12)

where $r=‖x-x^{\prime }‖$, l is the relevance determination; ${\sigma }_f^2$ is the variance and is used to control the degree of local relevance; ${\delta }_{x{x}^{’}}$ is the Kronecker’s delta function for x and $x^{\prime }$; the $\theta =\{l,{\sigma }_f^2\}$ is called the hyperparameter set of the kernel function.

Usually, hyperparameters $\theta$ have a significant impact on the training process and forecasting outcomes. It can be iteratively searched by the usage of the maximum likelihood method to build a log marginal likelihood function $L(\theta )$ of a training dataset. The $L(\theta )$ is calculated by:

$$L(\theta )=0.5{\mathit{y}}^T{C}^{-1}\mathit{y}+0.5\log \left|C\right|+0.5 n\log 2\pi$$

(13)

After that, the partial derivative of the hyperparameters can be obtained as:

$$\begin{array}{l}\frac{\partial L(\theta )}{\partial {\theta }_i}=0.5tr\left(((C^{-1}\mathit{y}){(C^{-1}\mathit{y})}^T-C^{-1})\frac{\partial C}{\partial {\theta }_i}\right)\\ C=K_n+{\sigma }_n^2{I}_n\end{array}$$

(14)

At last, the optimum $\theta$ is sought by minimizing Equation (14) via the conjugate gradient (CG) optimization. Using the best $\theta$ values, the mean value $E\left(y_{test}\right)$ and the posterior variance $\mathrm{cov}(y_{test})$ can be obtained using Equation (7).

Based on the sigma rule, the predicted interval under a given confidence level $1-\alpha$ is:

$$(E\left(y_{test}\right)-\frac{\mathrm{cov}(y_{test})}{\sqrt{n}}{z}_{\alpha /2},E\left(y_{test}\right)+\frac{\mathrm{cov}(y_{test})}{\sqrt{n}}{z}_{\alpha /2})$$

(15)

2.3.3. Physical-Based Kernel Function

A physical-based kernel function using several common kernels is designed for the GPR to model solar radiation. Generally, solar radiation data are periodic; therefore a periodic kernel is adopted. Additionally, to capture the amplitude of data changes in a day, the SE kernel is combined with the periodic function. Moreover, a linear and RN kernel is also considered to describe linearity and noise in solar radiation. Therefore, the new physical-based kernel function is defined as:

$$k_{CK}={\theta }_1^2\exp \left(-r/2{\theta }_2^2-2\sin \left(\frac{\pi r}{{\theta }_4^{}}\right)/{\theta }_3^2\right)+\hspace{0.33em}{\theta }_5^2\left(x-c\right)\left(x’-c\right)+{\theta }_6^2{\delta }_{x{x}^{’}}$$

(16)

To avoid too strong a dependence on the optimization result on the initial value, and easily plunging into the local optima of the standard CG method, BSA is adopted to find the optimum hyperparameters, ${\theta }_{i}s$, of the combined covariance function. The principle of BSA is introduced in the next subsection.

2.4. Backtracking Search Optimization Algorithm

BSA [33] is an advanced optimization method to solve complicated numerical issues. BSA utilizes well-known operators, such as selection, mutation, and crossover, of the genetic algorithms (GAs) in a new structure. In addition to these operators, BSA incorporates a number of novel strategies, such as a memory that stores a population from a historical generation selected at random. Civicioglu’s research [33] reveals that there are five key steps involved in putting BSA into action. The following is a simple introduction of these five stages:

(a) Initialization

In this stage, individuals of two populations are randomly initialized in the search space by a uniform distribution U, as shown below:

$$P_e^{i,j}=U(lo{w}_j,u{p}_j),\hspace{0.33em}\hspace{0.33em}{P}_{old}^{i,j}=U(lo{w}_j,u{p}_j),\hspace{0.33em}\hspace{0.33em}i=1,2,\cdots ,N_{pop};j=1,2,\cdots ,D$$

(17)

where $P_{old}$ is the historical population, $P_e$ presents the evolution population, $N_{pop}$ and D are the size of the population and the dimension of variables, respectively, and $lo{w}_j$ and $u{p}_j$ are the predefined lower and upper limitations of variables to be optimized.

(b) Selection-I

Prior to each iteration, an option is used to update the $P_{old}$ using the following "if-then" rule:

$$\begin{array}{l}if\hspace{0.33em}{R}_1<R_2\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}then\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{P}_{old}^{i,j}=P_e^{i,j},R_1,R_2\in U(0,1)\\ P_{old}^{}:=permuting(P_{old}^{})\end{array}$$

(18)

where $R_1$ and $R_2\hspace{0.33em}$ are two random numbers subject to uniform distribution in the range [0,1], and they are used to judge whether $P_{old}$ should be replaced by $P_e$ in the current generation; $:=$ is an operator for update; $permuting\left(\cdot \right)$ is a function that alters the order of individuals in $P_{old}$.

(c) Mutation

The initial format of the trial population $P_t$ is generated by:

$$\begin{array}{l}{P}_t=P_e+F\cdot (P_{old}-P_e)\\ F=3\cdot rndn,rndn\sim N(0,1)\end{array}$$

(19)

where F is a random number and controls the magnitude of the search direction matrix $(P_{old}-P_e)$.

BSA can take some of the experiences from antecedent generations, due to the utilization of the $P_{old}$ in the mutation process.

(d) Crossover

In this stage, the final format of $P_t$ can be obtained. The crossover operator is expressed by:

$$P_t{}^{i,j}=\left\{\begin{array}{c}{P}_e^{i,j},\hspace{0.33em}\hspace{0.33em}if ma{p}_{ij}=1\\ P_t{}^{i,j},otherwise\end{array}\right.\hspace{0.33em}\hspace{0.33em}$$

(20)

where $map$ is a binary matrix with size $N_{pop}\times D$, and determines the directions of crossover operator.

(e) Selection-II

A greedy selection strategy is utilized to create a new population of the next generation. This means that the individuals in the $P_t$ who have higher fitness values will be chosen to update the individuals in the $P_e$ that correspond to them.

$$P_e^{i,j}=\left\{\begin{array}{c}{P}_t{}^{i,j},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}if fitness(P_t{}^{i,j})\hspace{0.33em}<fitness(P_e^{i,j})\\ P_e^{i,j},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}otherwise\end{array}\right.\hspace{0.33em}\hspace{0.33em}$$

(21)

2.5. Flowchart of CEEMDAN-BSA-PGPR (CBP) Model

The implementation of the developed CBP model is given in Figure 1 and summarized as follows:

• Step 1: For a nonlinear and nonstationary solar radiation series $S\left(t\right)$, employ CEEMDAN to decompose it into a residue, $R\left(t\right)$, and several intrinsic mode functions, (IMFs) $IM{F}_k\left(k=1,2,\cdots ,K\right)$;
• Step 2: For each component ($IM{F}_k$ or $R\left(t\right)$), apply the PACF to determine suitable input variables;
• Step 3: Based on the selected input, establish PGPR for each component ($IM{F}_k$ or $R\left(t\right)$) and apply the BSA to obtain the best hyperparameters of the PGPR model;
• Step 4: The well-tuned models are utilized separately to do multi-hour-ahead forecasting using the testing dataset for every IMF and the residue;
• Step 5: The final forecasting values of the original solar radiation are obtained by aggregating all prediction results of all the IMFs and the residue.

3. Case Studies

3.1. Data Collection

Solar radiation data, from 2017 to 2018 and with a 1 h interval, are collected from Plataforma Solar de Almeria (PSA), which is located at 37.094416° N, 2.359850° W 490.7 m in Spain, and are used. The data are obtained from the SolarGIS database (https://solargis.com/cn/products/time-series-and-tmy-data/useful-resources, accessed on 15 October 2022), which is obtained by monitoring information from MTSAT and HIMAWARI satellites. Four datasets representing different seasons are utilized for experimental analysis. Figure 2 presents hourly GHI in W/m² (4663h in total) from different seasons, where Datasets 1–4 represent Spring, Summer, Autumn, and Winter. In Figure 2, the training and testing datasets are represented by the blue and red curves, respectively. Data with a 0 value at night have been deleted.

Table 1. Statistics for all datasets.

Datasets	Ave. (W/m2)	Std. (W/m2)	Cv.	Cs.	Max. (W/m2)	Min. (W/m2)	Sample Size	Data Period
Spring	526.4	340.0	0.65	−0.13	1029	1	1327	2017/3/20~2017/6/20
Summer	483.0	326.7	0.68	−0.03	1033	1	1347	2017/6/21~2017/9/21
Autumn	337.4	224.3	0.66	0.08	804	1	998	2017/9/22~2017/12/20
Winter	300.2	216.9	0.72	0.36	817	1	991	2017/12/21~2018/3/19

summarizes the statistics for the four datasets.

3.2. Parameter Setting

To comprehensively analyze the efficacy of the proposed CBP model, many commonly used models are constructed to do comprehensive comparisons. These models include BP, GRNN, RBF, ELM, SVR, GPR-SE, PGPR, BSA-PGPR, and DT-PGPR. The number of neurons in the hidden layers for BP and ELM is searched by the grid search (GS) in [1,100] with an interval of 1 [34]. The parameters of SVM, including kernel parameter λ and error penalty ρ, are also determined by the GS [34]. According to Reference [21], the parameters of GRNN and RBF are determined by a five-fold cross-validation method (CV) in (0,1) and [0.5, 1.5]. Besides the above five standalone models, the GPR with a SE kernel function has been established to attest to the effectiveness of the physical-based kernel function applied to the GPR model.

Considering the inherent non-stationarity and nonlinearity of the original solar radiation series, for the proposed CBP model, the CEEMDAN is first executed to divide the initial data into several IMFs and a residue. The related parameters of the CEEMDAN are as shown below: the realization number N is 20; the white noise amplitude is 0.2; the maximum iterations is 500 [35]. Figure 3 depicts the 11 subseries for the Spring dataset.

3.3. Determination of Input Variables

Input feature selection is a vital aspect of ML modelling. In general, if the input features are insufficient, the forecasting results may not be satisfactory [5]. On the contrary, too many inputs may lead to an additional computational burden. The antecedent solar radiation will be utilized as an input vector for modelling and forecasting. The optimal lag time of an input which has a significant impact on the target solar data needs to be determined. The PACF method, a widely used method, is adopted, and its principle can be found in Reference [36]. The PACF values for 24 lags with a 95% confidence interval (CI) for all datasets are calculated. As for the Spring dataset, the correlograms are presented in Figure 4.

3.4. Performance Criteria

In order to conduct a synthetical evaluation of the proposed model’s forecasting performance, four widely used evaluation criteria, Pearson’s correlation coefficient (R), root mean square error (RMSE), mean absolute error (MAE), and mean average percentage error (MAPE) [22], are utilized. These indices are calculated by

$$R=\left(\frac{{\displaystyle {\sum }_{i=1}^N\left(S_{obs,i}-\overline{S}{}_{obs}\right)\left(S_{fore,i}-\overline{S}{}_{fore}\right)}}{\sqrt{{\displaystyle {\sum }_{i=1}^N{\left(S_{obs,i}-\overline{S}{}_{obs}\right)}^2}}\sqrt{{\displaystyle {\sum }_{i=1}^N{\left(S_{fore,i}-\overline{S}{}_{fore}\right)}^2}}}\right),\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}-1<R<1$$

(22)

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{\left(S_{obs,i}-S_{fore,i}\right)}^2}}$$

(23)

$$\mathrm{MAE}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N\left|S_{obs,i}-S_{fore,i}\right|}$$

(24)

$$\mathrm{MAPE}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N\left|\frac{S_{obs,i}-S_{fore,i}}{S_{obs,i}}\right|}$$

(25)

where $S_{obs,i}$ and $S_{fore,i}$ are the i-th measured and forecasted solar radiation, respectively, $\overline{S}{}_{obs}$ and $\overline{S}{}_{fore}$ are the mean of the observed and computed solar radiation data, respectively, and N is the whole number of samples.

Improved percentage indices of these four metrics between model 1 and model 2 are calculated by [25]:

$$P_{index}=\left(inde{x}_1-inde{x}_2\right)/inde{x}_1$$

(26)

4. Results

The 10 models, including BP, GRNN, RBF, ELM, SVR, GPR-SE, PGPR, BSA-PGPR, CEEMDAN-PGPR (CEN-PGPR), and CBP, established in this study are utilized on the four datasets for from one-hour- to three-hour-in-advance solar radiation forecasting. These models are evaluated via four indices. Additionally, other results analysis approaches, such as histograms, scatter, radar graphs, and area graphs are also adopted to evaluate model performances more visually. In this section, the developed model, CBP, is compared with the other benchmark models, including seven standalone models (BP, GRNN, RBF, ELM, SVR, GPR-SE, and PGPR) and two hybrid models (BSA-PGPR and CEN-PGPR).

4.1. Results of 1-Hour Ahead Forecasting

The four metrics provided by the proposed model, CBP, and nine contrast models on the testing dataset for all datasets are presented in

Table 2. Evaluation indices for one-hour-ahead forecasting.

Datasets	Metrics	BP	GRNN	RBF	ELM	SVR	GPR-SE	PGPR	BSA-PGPR	CEN-PGPR	CBP
Spring	R	0.9879	0.9735	0.9747	0.9764	0.9814	0.9880	0.9831	0.9881	0.9965	0.9973
	RMSE (W/m2)	54	85	81	79	67	54	63	54	29	26
	MAE (W/m2)	31	66	60	59	42	31	34	31	17	17
	MAPE (%)	1.42	3.14	3.05	2.93	1.97	2.15	2.32	2.18	1.02	0.95
Summer	R	0.9790	0.8785	0.9464	0.9518	0.9759	0.8822	0.9728	0.9762	0.9960	0.9969
	RMSE (W/m2)	64	177	109	103	69	159	73	71	28	25
	MAE (W/m2)	41	115	77	75	44	92	49	46	19	17
	MAPE (%)	1.43	5.30	3.66	4.11	1.25	5.59	1.52	1.66	0.83	0.79
Fall	R	0.9420	0.8491	0.8633	0.9027	0.9272	0.8922	0.9300	0.9405	0.9829	0.9857
	RMSE (W/m2)	60	111	92	81	71	84	67	62	33	31
	MAE (W/m2)	43	92	68	64	54	59	49	44	24	23
	MAPE (%)	0.52	1.00	1.08	0.71	0.60	0.85	0.60	0.64	0.31	0.30
Winter	R	0.9202	0.8663	0.8532	0.8486	0.8928	0.9235	0.9246	0.9255	0.9796	0.9817
	RMSE (W/m2)	100	135	144	134	113	97	96	95	51	50
	MAE (W/m2)	77	108	118	109	85	73	73	71	39	37
	MAPE (%)	3.41	4.18	7.47	5.93	3.42	4.35	4.01	3.73	2.00	1.92

. Additionally, the indices values of the 10 models are visually presented in Figure 5. M1 to M10 are BP, GRNN, RBF, ELM, SVR, GPR-SE, PGPR, BSA-PGPR, CEN-PGPR and CBP, respectively.

According to Table 2 and Figure 5, the forecasting performance of the Winter dataset is the poorest for almost all cases, then the Summer, followed by the Spring, and, lastly, the Autumn dataset. This may be due to the gradual increase in the value of the testing dataset in the Winter dataset compared to the training dataset (Figure 2). In addition, the Summer solar radiation testing dataset has obvious volatility which is different from that of the training data, which increases the difficulty of prediction (Figure 2).

Based on a thorough analysis of Table 2 and Figure 5, a few conclusions can be drawn: (a) Among seven single models, the PGPR model produces the highest R-value and the lowest RMSE, MAE, and MAPE values compared to the other six standalone models for all datasets, except for the Spring dataset. For the Spring dataset, the evaluation indices of the PGPR are close to those of the best single model, indicating that the PGPR model enhances the interpretability without sacrificing precision. In contrast, the GRNN performs the worst. For the Spring dataset, the R, RMSE, MAE, and MAPE of the PGPR model are 0.9831, 63 W/m², 34 W/m², and 2.32%, respectively. Indices of the GRNN are 0.9735, 85 W/m², 66 W/m², and 3.14%, respectively. This result demonstrates the effective nonlinear fitting ability of the PGPR model. (b) Among all hybrid and single models, the proposed CBP model generated the lowest RMSE, MAE, and MAPE values and the highest R-value. For the Spring dataset, the R, RMSE, MAE, and MAPE of the CBP model are 0.9973, 26 W/m², 17 W/m², and 0.95%, respectively, indicating an excellent forecasting performance. (c) The four indices of the BSA-PGPR are better than those of the PGPR for all datasets. The four metric values of the CBP are always superior to those of the CEN-PGPR for all four datasets. From the comparison of the PGPR (CEN-PGPR) and BSA-PGPR (CBP) models, it is shown that the BSA can be a useful tool to search the optimum hyperparameters of the PGPR model. (d) The comparison of the CBP and BSA-PGPR models indicates that the CEEMDAN can be adopted as a capable alternative in improving the accuracy of solar radiation forecasting.

To intuitively present the forecasting performance and further testify to the efficiency of the CBP model, the single SVR, the best standalone model PGPR, and the hybrid model BSA-PGPR are used for further analysis. The one-hour-ahead forecasting results of these four models for the four datasets are portrayed in Figure 6, Figure 7, Figure 8 and Figure 9. A line chart at the bottom of these graphs displays forecasting errors between actual and predicted values of solar radiation. At the top of these graphs are forecasting and real curves. When compared to the forecasting curves provided by SVR, PGPR, and BSA-PGPR, it is evident that the forecasting curves of the CBP model are the closest to the observed curves. Additionally, the forecasting errors of the proposed CBP are around zero, with a narrow range of variation, significantly lower than those of other contrast models, especially for extreme situations. As for the Spring dataset, the residual value of the proposed CBP is relatively small, that of the PGPR has an obvious underestimation, and that of the BSA-PGPR has a significant overestimation between 350 hours and 375 hours (marked in Figure 6 as a blue ellipse). According to the scatter plots (Figure 10), it is shown that the forecasted solar radiation of the CBP model is uniformly scattered around the 45° line with a narrower range in the description of solar radiation values with different amplitudes than other contrast models, which graphically illustrates the superiority of the CBP model.

The prediction intervals (PIs) with two confidence levels obtained by the proposed CBP are visually depicted in Figure 11. There are a few extreme fluctuations in a small number of samples that fall outside of 90% PIs for all datasets. Furthermore, the fact that the PIs for these two confidence levels are sufficiently narrow draws the conclusion that the proposed CBP model is capable of producing the applicable PI.

4.2. Results of Multi-Hour Ahead Forecasting

In practical application, multi-hour-ahead forecasting of solar radiation is also required. The multi-hour-in-advance forecasting information is of great help to photovoltaic grid connection and power dispatching. Thereby, the CBP model is also executed to do two- and three-hour-ahead forecasting. The forecasting results of the developed CBP model are compared with the other nine benchmark models. Two-hour- and three-hour-ahead forecasting results are given in

Table 3. Evaluation indices for two-hour-ahead forecasting.

Datasets	Metrics	BP	GRNN	RBF	ELM	SVR	GPR-SE	PGPR	BSA-PGPR	CEN-PGPR	CBP
Spring	R	0.9718	0.9478	0.9579	0.9486	0.9436	0.9750	0.9723	0.9751	0.9915	0.9927
	RMSE (W/m2)	86	117	105	114	116	84	86	81	46	43
	MAE (W/m2)	57	95	83	92	84	59	56	54	31	29
	MAPE (%)	3.70	3.44	3.69	6.06	6.34	3.56	3.62	3.26	1.90	1.83
Summer	R	0.9480	0.8111	0.8940	0.8998	0.9292	0.8488	0.9336	0.9435	0.9864	0.9927
	RMSE (W/m2)	109	219	150	149	127	183	116	110	51	37
	MAE (W/m2)	73	141	107	111	91	122	83	74	38	27
	MAPE (%)	3.76	6.04	7.04	5.59	3.71	7.43	4.04	3.62	1.39	1.36
Fall	R	0.8367	0.7049	0.7769	0.7868	0.7458	0.8015	0.7988	0.8457	0.9633	0.9697
	RMSE (W/m2)	98	152	116	119	134	110	112	97	49	45
	MAE (W/m2)	73	115	86	96	107	82	84	73	36	35
	MAPE (%)	0.96	1.44	1.29	1.12	1.14	1.27	1.07	1.05	0.52	0.52
Winter	R	0.8254	0.7612	0.7267	0.7245	0.7369	0.8299	0.8346	0.8408	0.9622	0.9654
	RMSE (W/m2)	145	170	176	177	170	142	139	137	70	68
	MAE (W/m2)	116	141	147	146	139	115	112	110	53	53
	MAPE (%)	5.03	4.91	8.50	5.97	7.55	5.59	5.06	4.88	3.10	2.61

and

Table 4. Evaluation indices for three-hour-ahead forecasting.

Datasets	Metrics	BP	GRNN	RBF	ELM	SVR	GPR-SE	PGPR	BSA-PGPR	CEN-PGPR	CBP
Spring	R	0.9497	0.9242	0.9383	0.9173	0.9040	0.9571	0.9587	0.9585	0.9805	0.9813
	RMSE (W/m2)	112	139	132	139	147	112	110	106	68	67
	MAE (W/m2)	78	116	108	117	117	85	84	78	48	47
	MAPE (%)	2.27	3.58	4.33	5.89	6.31	2.97	3.00	3.08	2.46	1.95
Summer	R	0.9046	0.7885	0.8486	0.8460	0.8589	0.8152	0.8806	0.9037	0.9719	0.9845
	RMSE (W/m2)	144	236	184	182	179	207	160	147	73	54
	MAE (W/m2)	100	158	135	141	138	145	117	105	56	40
	MAPE (%)	6.13	7.29	8.49	7.54	7.52	8.07	6.77	6.52	2.30	1.90
Fall	R	0.7693	0.6666	0.7617	0.6781	0.5326	0.7781	0.6881	0.7924	0.9314	0.9437
	RMSE (W/m2)	115	152	116	141	175	113	137	110	67	61
	MAE (W/m2)	88	119	91	118	148	85	107	83	50	48
	MAPE (%)	0.88	1.40	1.21	1.37	1.67	1.19	1.37	1.12	0.76	0.70
Winter	R	0.7159	0.6372	0.6747	0.6224	0.5339	0.7258	0.7341	0.7406	0.9268	0.9346
	RMSE (W/m2)	177	198	191	202	219	173	171	169	97	93
	MAE (W/m2)	144	168	163	167	180	144	141	139	79	76
	MAPE (%)	6.08	6.61	8.13	6.03	9.08	6.11	5.27	4.69	5.03	3.53

, respectively. Moreover, the corresponding heatmap for R, vertical histograms for RMSE, radar plots for MAE, and area graphs for MAPE, for two-hour- and three-hour-ahead forecasting are visually presented in Figure 12 and Figure 13, respectively. From these results, listed in Table 3 and Table 4 and Figure 12 and Figure 13, it is reasonable to draw the same conclusion as Section 4.1. Therefore, the effectiveness and superiority of the proposed CBP have been attested.

5. Discussion

In light of the preceding comparative analysis of one-hour- to three-hour-ahead solar radiation forecasting results, the impacts of the CEEMDAN on smoothing raw data and the BSA searched the optimal hyperparameters on forecasting performance are further discussed. The improvements of R, RMSE, MAE, and MAPE, provided by the model with decomposition technique and without decomposition technique, are further compared and listed in

Table 5. Improvements of models with decomposition technique relative to that without decomposition technique.

Datasets	Improvements (%)	CEN-PGPR vs. PGPR			CBP vs. BSA-PGPR
		1-hour	2-hour	3-hour	1-hour	2-hour	3-hour
Spring	PR	1.37	1.97	2.27	0.93	1.80	2.37
	PRMSE	54.09	47.24	37.84	52.49	47.82	37.04
	PMAE	47.85	44.63	42.30	44.75	45.69	39.26
	PMAPE	56.06	47.44	17.87	56.32	43.84	36.61
Summer	PR	2.39	5.66	10.37	2.13	5.22	8.94
	PRMSE	62.41	56.42	54.76	65.00	66.38	62.96
	PMAE	61.17	54.03	52.31	62.37	63.96	62.10
	PMAPE	45.72	65.58	66.00	52.80	62.57	70.86
Fall	PR	5.69	20.60	35.35	4.81	14.67	19.10
	PRMSE	51.32	56.45	51.32	49.32	53.54	44.63
	PMAE	51.85	56.81	53.36	48.52	51.28	42.56
	PMAPE	47.99	51.74	44.65	53.53	50.87	37.16
Winter	PR	5.96	15.30	26.25	6.08	14.82	26.20
	PRMSE	46.91	50.00	43.13	47.93	49.94	45.00
	PMAE	46.67	52.41	44.30	47.30	52.15	45.04
	PMAPE	50.16	38.81	4.59	48.68	46.46	24.76

. All the improved percentages are over zero, which indicates that the CEEMDAN has a significant positive impact on strengthening the forecasting performance of solar radiation. Compared with the BSA-PGPR method, the P_RMSE of the proposed CBP model for four datasets are 52.49%, 65%, 49.32%, and 47.93%, respectively, for one-hour-ahead forecasting; 37.04%, 62.96%, 44.63%, and 45%, respectively, for three-hour-ahead forecasting. By embedding CEEMDAN, the prediction accuracy of the CBP is more attractive relative to the BSA-PGPR, and the average promotion rate in terms of RMSE is approximately 50%.

The four improvement percentage metrics, including P_R, P_RMSE, P_MAE, and P_MAPE, obtained by the model with the optimization algorithm and without the optimization algorithm are listed in

Table 6. Improvements of models with BSA optimization relative to those without BSA optimization.

Datasets	Improvements (%)	BSA-PGPR vs. PGPR			CBP vs. CEN-PGPR
		1-hour	2-hour	3-hour	1-hour	2-hour	3-hour
Spring	PR	0.51	0.29	−0.02	0.07	0.12	0.08
	PRMSE	13.58	5.54	3.17	11.8	7.05	1.95
	PMAE	8.02	4.02	7.45	2.63	6.22	2.65
	PMAPE	5.93	9.89	−2.68	6.95	3.86	26.2
Summer	PR	0.35	1.06	2.62	0.09	0.64	1.28
	PRMSE	3.81	5.06	8.56	11.66	36.52	33.59
	PMAE	7.62	10.04	9.77	11.71	41.78	39.45
	PMAPE	−9.43	10.40	3.81	5.1	2.65	21.31
Fall	PR	1.13	5.87	15.15	0.29	0.66	1.3
	PRMSE	8.60	13.15	19.81	5.09	7.95	9.63
	PMAE	10.33	13.29	22.64	4.3	2.24	4.95
	PMAPE	−8.02	1.67	18.82	3.61	0.11	8.51
Winter	PR	0.10	0.74	0.88	0.21	0.33	0.83
	PRMSE	1.04	1.84	0.92	3.01	1.73	4.35
	PMAE	2.13	1.86	1.78	3.39	1.34	3.18
	PMAPE	6.94	3.66	11.07	4.36	18.63	42.6

. It indicates that the metric values of BSA-PGPR are slightly poorer than those of PGPR in a handful of the cases, such as three-step forecasting for Spring, and one-step forecasting for Winter and Spring. In other cases, the BSA-PGPR performs well, while the R, RMSE, MAE, and MAPE of CBP are superior to those of CEN-PGPR in all cases. This phenomenon indicates that the positive effect of BSA can be amplified by introducing the decomposition technique and highlights the advantages of combining CEEMDAN and BSA.

6. Conclusions

This study develops and attests the forecasting ability of a compound framework based on an improved GPR and CEEMDAN (called CBP) for multi-step hourly solar radiation prediction. In the improved GPR, a newly physical-based kernel is designed to enhance forecasting performance, and BSA is applied to optimize hypermeters of the physical-based GPR (PGPR). To attest the superiority of the proposed model, nine models, namely BP, GRNN, RBF, ELM, SVR, GPR-SE, PGPR, BSA-PGPR, and CEN-PGPR are established. The following is a summary of the study’s main findings:

(1)

The proposed CBP model provided the best results, with the biggest R and lowest RMSE, MAE, and MAPE, which attested to the excellent forecasting ability of the CBP relative to the other comparative models;

(2)

The forecasting results produced by the PGPR model are close to or superior to those of the GPR model, indicating that the newly designed kernel function enhances the interpretability of the GPR model while maintaining accuracy;

(3)

The proposed PGPR model outperformed many widely used standalone BP, GRNN, RBF, ELM, and SVR models, which revealed the effectiveness of the CBP model relative to other comparative models;

(4)

The BSA, applied in searching appropriate hyperparameters of the PGPR, along with the use of CEEMDAN to extract detailed information, can enhance the forecasting ability of CBP. Using BSA alone to optimize PGPR parameters does not always achieve an ideal forecasting performance.

In a word, the proposed CBP model can be a viable option for improving solar radiation forecasting. In the future, the prediction ability of the CBP can be further validated by employing other solar radiation data with various time scales and for other locations. The CBP model will also be used with other kinds of data, such as wind energy data.