The Effect of Using Aerosol Variables on the Performance of Deep Learning-based GHI Forecasting Models

1. Introduction

The growth in solar energy adoption worldwide is motivated by the world’s commitment to achieving the sustainable development goals of the United Nations in producing clean energy and the declining cost of solar power plants and their technology. As per expectations, solar will become the primary energy source by 2050 [1]. Saudi Arabia, as one of the leading oil producers, has planned to shift the country’s dependence on oil and adopt renewable energy sources. Since the country is in the six place of the world in potential for producing solar energy [2], it has launched a series of projects that aims to produce 9.5 GW of renewable energy by the end of 2023 [3]. Saudi Arabia’s first utility-scale solar PV project is the Sakaka 300-megawatt solar power station, which was linked to the national grid in November 2019.

Weather change challenges solar energy production and makes PV power generation instable. Managing smart grids and energy markets efficiently requires short-term accurate forecasting of the power output of PV systems. This has encouraged researchers to work on innovative methods to forecast solar radiation. Currently, three types of methods are used for this task: numerical and simulation weather prediction models, statistical models, and artificial intelligence (AI) based models [4,5]. Lately, DL methods have shown superior performance compared to other types of forecasting methods [6,7,8,9].

In earlier work [10], the performance of several DL-based forecasting models in next-hour GHI forecasting was compared using lagged features of meteorological variables and solar radiation measurements. The focus was on comparing the forecasters’ performances in multiple locations with various climates. In this work, the performance of five new DL models is compared carrying out the same task but using different sets of features, including weather and aerosol variables along with solar radiation components. The goal is to understand the effect of using various combinations of features on forecasting accuracy. Also, two data sources are used for the same location to validate the results. The contributions of this paper are summarized as follows.

• Studying the effect of using aerosol variables on the performance of five new DL-based models for a next-hour GHI forecasting task using data from a location with a hot desert climate
• Using two different data sources, ground-based and satellite-based to validate the forecasting results.
• Presenting the forecasting results using visualization and several evaluation metrics, including root mean square error (RMSE), mean absolute error (MAE), mean absolute percentage error (MAPE), and forecast skills (FS).

The rest of this paper is structured as follows. Section 2 discusses the related work and identifies the research gap. Then, Section 3 explains the methodology, including the data preprocessing, models’ development, implementation, and evaluation metrics. Next, Section 4 discusses GHI forecasting results using different sets of features, which include lagged GHI values, weather variables, and aerosol measurements. Finally, Section 5 concludes the work.

2. Related Work

Alkhayat and Mehmood conducted an widespread literature review on DL-based solar energy forecasting methods in [11]. In most of these studies, lagged GHI values and weather variables, such as air temperature (AT), wind speed (WS), and wind direction (WD), are used as features for prediction. The effect of using exogenous inputs, such as weather data, in addition to endogenous inputs, such as historical GHI values, on the performance of statistical solar irradiance forecasting models has been widely studied, as shown in the literature [12,13,14,15]. This has also been studied in relation to machine learning models (ML) [16,17,18,19,20], and a few studies have compared DL-based models’ performance with and without exogenous inputs, such as [21,22,23,24]. For example, Lee et al. [21] compared the performance of an LSTM-AE model and a hybrid model of convolutional neural network (CNN) and long short-term memory (LSTM) with and without weather inputs. Their results show improvement in forecasting using weather inputs. Similarly, Castangia et al. [22] compared five ML models’ performances, including the CNN model and LSTM models, for short-term GHI forecasting with exogenous and endogenous inputs. Their results show improvement in the models’ performances in relation to the following eight features: UV index, cloud cover (CC), AT, relative humidity (RH), dew point (DP), wind bearing, hour of the day, and sunshine duration. Likewise, Omar et al. [23] compared the performance of an LSTM model and a radial basis function neural network (RBFNN) with exogenous and endogenous inputs. They found that the LSTM model performed better without exogenous variables, such as pressure (P), zenith angle (ZA), AT, and RH. In [24], Omar et al. used a novel feature selection method called weather recursive feature elimination along with an LSTM model for GHI forecasting. They found that the model achieved the best performance with GHI, RH, direct normal irradiance (DNI), and diffuse horizontal irradiance (DHI).

Aerosol measurements, as one of the exogenous inputs, have been used in the literature to improve physical and statistical GHI forecasting models [25,26,27,28,29,30]. However, few studies have included aerosol measurements in the features to develop ML models for the GHI prediction [31,32,33,34,35]. For example, Alfadda et al. [31] tested four ML methods, namely, multilayer perceptron (MLP), decision tree regression (DTR), support vector regression (SVR), and k-nearest neighbors (kNN), for hourly GHI prediction using features that include two aerosol measurements: aerosol optical depth (AOD) and the angstrom exponent (AE). Since the data used in their work was gathered from Riyadh in Saudi Arabia, which is a high turbidity location, they found that the inclusion of aerosol measurements in the features improved GHI forecasting. Similarly, Kumar et al. [32] used AOD to develop their artificial neural network (ANN) model to predict the next 3-hour GHI in Delhi, India, whereas Zuo et al. [33] used it along with other meteorological parameters to develop an LSTM model to forecast the next 10-minute (min) GHI in China. In addition, Si et al. [34] developed a hybrid model for solar irradiance prediction that utilizes satellite images and meteorological data, including AOD. Their model combines a CNN model to extract features from sky images and an MLP model to forecast GHI in Shandong province, China. Also, Zhu et al. [35] used an ensemble model to predict the next 10-min GHI, which consists of a multiple regression model, an SVR, and an MLP. The effect of AOD and another nine meteorological variables on the model’s performance was studied in this work. 

Table 1 summarizes the related studies by identifying the forecasting method, features used as inputs, data source type, and the main results.

Research Gap

Two research gaps are recognized in this section. The first is the lack of comparative studies on DL-based GHI forecasting models’ performance with and without exogenous inputs, as what has been carried out in the literature with statistical models [12,13,14,15] and traditional ML models [16,17,18,19,20]. The second is the dearth of research on the effect of aerosol measurements, as one of the exogenous inputs, on ML-based GHI forecasting models’ performance, as what has been done with physical and statistical GHI forecasting models [25,26,27,28,29,30]. Thus, there is a need to perform more studies on the performance of DL-based models with various sets of features, including weather and aerosol measurement inputs, which is the aim of athis work. To the best of our knowledge, this is the first paper that studies the effect of using aerosol variables on the performance of DL-based models for a next-hour GHI forecasting task using data from a location with a hot desert climate.

3. Methodology

First, Section 3.1 describes the data preprocessing steps, including data collection, data cleaning, feature extraction, and data normalization and dividing. Next, Section 3.2 explains the development process for five DL models, which are LSTM, gated recurrent unit (GRU), bidirectional LSTM (BiLSTM), bidirectional GRU (BiGRU), and LSTM-AE. Section 3.3 provides the implementation details of all the models. Finally, Section 3.4 clarifies the performance evaluation metrics.

3.1. Data Preprocessing

Three data preprocessing steps are explained in this section: data collection, data cleaning and feature extraction, and data normalization and dividing.

3.1.1. Data Collection

Four datasets were used, all of which were collected at the same location, which is King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia (22.305° N 39.103° E). The location of KAUST can be seen on the map of Saudi Arabia shown in Figure 1. The Köppen climate classification of this location is hot desert climate (BWh). Figure 2 shows the monthly averages of temperature and precipitation for the years 1990-2019 according to ClimateCharts.net [36].

The first dataset was gathered by King Abdullah City for Atomic and Renewable Energy (K.A.CARE)[37] at the mid-range (tier 2) station. This contains a rotating shadow band radiometer (RSR), which provides fundamental solar resource data, plus basic meteorological instruments. These ground-based measurements were taken with a resolution of 1 min and nominal uncertainty of +/−5% (sub-hourly) [38]. 1-min readings were averaged into hourly readings. The dataset covers the period from 1 January 2016 to 31 March 2021; however, observations on some days are missing because of maintenance scheduling or devices failure. The total number of missing days is 290.

The second dataset was collected by the aerosol robotic network (AERONET) [39] with level 2.0 quality, which means the data is cloud-screened and quality-controlled. This dataset contains hourly ground-based measurements of AOD and AE. AOD measures the amount of direct sunlight that is blocked from reaching the ground by aerosol particles, such as dust, smoke, and pollution, which could absorb sunlight or cause it to scatter. A low value of AOD corresponds to a clean atmosphere, whereas a high value corresponds to hazy or dusty conditions [40]. AE quantifies the particle size of atmospheric aerosols or clouds, and it is inversely related to the average size of the particles in the aerosol. Consequently, low AE values suggest a strong presence of coarse aerosols relating to dust events [41]. This dataset covers the period from 1 January 2016 to 31 March 2021. However, 902 days are missing from the dataset.

The third dataset was collected by the National Solar Radiation Database (NSRDB), which is available on the website of the National Renewable Energy Laboratory [42]. Since this data is satellite-derived measurements, there are no missing values. The time resolution is 1 hour, and the spatial resolution is 4 km. This dataset spans from 1 January 2017 to 31 December 2019.

The fourth dataset was collected by NASA and accessed through the Goddard interactive online visualization and analysis infrastructure (GIOVANNI) website [43]. This data is the product of the modern era retrospective-analysis for research and applications, Version 2 (MERRA-2), which assimilates space-based observation of aerosols. The time resolution is 1 hour, and the spatial resolution is 50 km. This dataset spans from 1 January 2017 to 31 December 2019 with no missing values.

3.1.2. Data Cleaning and Feature Extraction

In this section, the data cleaning and feature extraction process of four datasets are detailed.

• K.A.CARE dataset

This dataset contains hourly values of nine attributes, as follows:

• Output: GHI as watt-hour per square meter (Wh/m²)
• DHI as Wh/m²
• DNI as Wh/m²
• ZA as degree °
• AT as Celsius (° C)
• WS taken at 3m as a meter per second (m/s)
• WD taken at 3 m as m/s
• Barometric pressure (BP) as Pascal (Pa)
• RH as a percentage (%)

The correlation between GHI and the GHI values of earlier hours on the same day was calculated. It is found that the correlation is high with only the values from the last three hours. GHI correlation with GHI_lag3 is 0.57 and with GHI_lag4 is 0.32. Therefore, the previous three hours’ records were used as inputs to predict the next hour’s GHI. Lagged features can be created with the shift method in the Pandas library or by reshaping inputs with timestep equals three. Figure 3 shows the correlation matrix of the K.A.CARE dataset after using the last three hours’ values of the nine attributes listed before as features. Correlations between GHI and DHI_lag1, DHI_lag2, DHI_lag3, DNI_lag1, DNI_lag2, DNI_lag3, GHI_lag1, GHI_lag2, and GHI_lag3 are strong positive correlations. On the other hand, correlations between GHI and ZA_lag1, and ZA_lag2 are strong negative correlations.

The GHI values for the same hour on previous days might also be important in predicting GHI. For example, today's GHI value at 12 p.m. could be highly correlated with the GHI value at 12 p.m. yesterday or the day before. Therefore, the correlation between GHI and its value for the same hour on the previous seven days was calculated, then on the previous 15 days, and lastly for the last months up to the previous year. Figure 4 shows the correlation matrix of the K.A.CARE dataset after using the previous day up to the previous year's same-hour GHI values as features. Surprisingly, the correlation of GHI and the previous year’s same-hour GHI value (GHI_360D) is equal to that of GHI and the previous hour value (GHI_lag1), which is 0.94. Also, all the previous year’s GHI values have a higher correlation with GHI than the correlation with the GHI value of the previous 2 hours (GHI_lag2).

The final number of features is 50, as listed in Table 2. As mentioned before, there are many missing values. Records where GHI is null or less than 1 (night hours) were eliminated. The interpolation method was used to fill some null values of GHI in the previous n days. For example, if the GHI value of the same hour on the previous day (GHI_1D) was missing, it was filled by the GHI value of the same hour on the previous 2 or 3 days (GHI_2D) or (GHI_3D). If both values were unavailable, then GHI_1D was filled by the average of GHI_4D, GHI_5D, and GHI_6D. The same method was used to fill all the missing values of GHI_1D up to GHI_360D. Even after using interpolation, many GHI previous n days values remain missing because there are no consecutive records to use for filling. For example, the first day in the dataset is January 1, 2016. Therefore, the GHI_1D feature for all that day’s records will remain missing because the previous days’ records are unavailable. These remaining null records were deleted after interpolation was used.

• AERONET dataset

This dataset contains the hourly values of four attributes, as follows:

• AOD at 500 nm (AOD_500)
• AOD at 551 nm (AOD_551)
• AE for the wavelength range from 440 to 675 nm (440-675_AE)
• Optical Air Mass (OAM)

Since the number of null values in AOD_551 is larger than the number of null values in AOD_500, it was decided to use AOD_500 as a feature and use the former only to fill null values in the latter. Then, the K.A.CARE dataset was merged with this dataset using date and time columns. Figure 5 shows the correlation matrix of the merged datasets. It is seen that the GHI correlation with AOD_500 and 440-675 AE is very weak, whereas the GHI correlation with OAM is moderate (coefficient= 0.43).

Merging the K.A.CARE dataset with the AERONET dataset resulted in 5731 hourly records. These records were further analyzed to see the effect of the existence or absence of dust on GHI and its components DNI and DHI. First, these records were filtered to meet the following two conditions: 440-675_AE is less or equal to 0.2 and AOD_500 is greater or equal to 0.9. A small value of AE indicates the existence of large coarse particles, such as dust while the opposite is true for AOD [41]. Records for forty-one hours met these thresholds for dusty weather. Out of these, thirty-three hours have a DHI over GHI value greater than 60%. This means that there is an 80% chance that when the weather is dusty, more than 60% of the GHI amount is scattered or diffused. Second, the records were filtered to meet the following two conditions: 440-675_AE is greater or equal to 1.4 and AOD_500 is smaller or equal to 0.2. A high value of AE and a low value of AOD indicate very clear weather [44]. Records for sixty-nine hours met these thresholds for clear weather. Out of these, fifty-six hours have DNI over GHI value greater than 70%. This means that there is an 81% chance that when the weather is very clear and there are no large particles to cause sunlight to scatter, more than 70% of the GHI amount is received in a direct line. Figure 6 (a) shows an example of a clear day in which DNI is greater than GHI for all hours. On the other hand, Figure 6 (b) shows an example of a dusty day in which DHI is almost equal to or a big part of the GHI value because dust causes sunlight to diffuse.

The final number of features is nineteen, as listed in Table 3. Compared to Table 2, thirty-one features were eliminated because of the large number of records missing after merging. Therefore, all lag 1 features were kept, and lag 2 and 3 features were removed. Also, regarding previous n days’ GHI values, only 4 days' GHI values were used to avoid having many null records.

Records where GHI is null or less than 1 (night hours) were eliminated. Then, the interpolation method was used to fill some null values of GHI in the last n days, as explained previously. The null records remaining after using interpolation were deleted.

Table 3. K.A.CARE & AERONET dataset features.

Time t features	Time t-1 features		Tim t features last n days
GHI (output)	GHI_lag1	WD_lag1	GHI_1D
	DNI_lag1	RH_lag1	GHI_2D
Hour	DHI_lag1	BP_lag1	GHI_3D
Day	ZA_lag1	AOD_500_lag1	GHI_4D
Month	AT_lag1	440-675_AE_lag1
	WS_lag1	OAM_lag1

Table 3. K.A.CARE & AERONET dataset features.

Time t features	Time t-1 features		Tim t features last n days
GHI (output)	GHI_lag1	WD_lag1	GHI_1D
	DNI_lag1	RH_lag1	GHI_2D
Hour	DHI_lag1	BP_lag1	GHI_3D
Day	ZA_lag1	AOD_500_lag1	GHI_4D
Month	AT_lag1	440-675_AE_lag1
	WS_lag1	OAM_lag1

• NSRDB dataset

This dataset contains hourly values of nine attributes, as follows:

• Output: GHI as w/m²
• DHI as w/m²
• DNI as w/m²
• ZA as degree °
• AT as Celsius (° C)
• WS as a meter per second (m/s)
• WD as m/s
• BP as Millibar
• RH as a percentage (%)

Figure 7 shows the correlation matrix of the NSRDB dataset after using the previous three hours’ records of the nine attributes listed before as features. The correlations between GHI and AT_lag1, DHI_lag1, DHI_lag2, DHI_lag3, DNI_lag1, DNI_lag2, DNI_lag3, GHI_lag1, GHI_lag2, and GHI_lag3 are strong positive correlations. On the other hand, the correlations between GHI and ZA_lag1, ZA_lag2, ZA_lag3, RH_lag1, and RH_lag2 are strong negative correlations. Compared to Figure 3, the correlations between GHI and AT or RH are not significant in the K.A.CARE dataset. This might be attributed to the incompleteness of this dataset in addition to the differences between the ground-based and satellite-based measurements.

Figure 8 shows the correlation matrix of the NSRDB dataset after using the previous day up to last year’s same-hour GHI values as features. As noted in Figure 4, all last year’s same-hour GHI values have a higher correlation coefficient than the GHI of 2 hours earlier on the same day.

The final number of features in the NSRDB dataset is 50, and they are the same features as those of the K.A.CARE dataset listed in Table 2.

• GIOVANNI NNI dataset

This dataset contains hourly values of ten attributes as follows:

• Dust extinction aerosol optical thickness 550 nm (DUEXTTAU)
• Dust extinction aerosol optical thickness 550 nm - PM 2.5 (DUEXTT25)
• Total aerosol extinction aerosol optical thickness 550 nm (TOTEXTTAU)
• Dust column mass density (DUCMASS) as kg m^-2
• Dust column mass density - PM 2.5 (DUCMASS25) as kg m^-2
• Dust surface mass concentration (DUSMASS) as kg m^-3
• Dust surface mass concentration - PM 2.5 (DUSMASS25) as kg m^-3
• Dust scattering aerosol optical thickness 550 nm - PM 1.0 (DUSCATFM)
• Total aerosol scattering aerosol optical thickness 550 nm (TOTSCATAU)
• Total Aerosol Angstrom parameter 470-870 nm (TOTANGSTR)

First, the NSRDB dataset was merged with this dataset using date and time columns.

Figure 9 shows the correlation matrix of the merged datasets. It is noticeable that GHI correlations with all ten aerosol attributes are very weak.

The total number of records for daylight hours after merging is 13169. These records were further analyzed to see the effect of the existence or absence of dust on GHI and its components DNI and DHI. First, these records were filtered to meet two conditions: TOTANGSTR is less or equal to 0.1, and TOTEXTTAU and TOTSCATAU are greater than or equal to 0.7. A small TOTANGSTR value indicates the existence of large coarse particles, such as dust, while the opposite is true for TOTEXTTAU and TOTSCATAU [45,46]. Records for seventy-one hours meet these thresholds for dusty weather. Out of these, fifty-seven hours have a DHI over GHI value greater than 60%. This means that there is an 80% chance that wahen the weather is dusty, more than 60% of the GHI amount is scattered or diffused. Second, the records were filtered to meet the two conditions: first, TOTANGSTR is greater or equal to 1.0 and second, TOTEXTTAU and TOTSCATAU are smaller than or equal to 0.2. A high value of TOTANGSTR and a low value of TOTEXTTAU and TOTSCATAU indicate very clear weather [44,47]. Records for ninety-seven hours met these thresholds for very clear weather. Out of these, eighty-five hours have a DNI over GHI value greater than 70%. This means that when the weather is very clear and there are no large particles to cause sunlight to scatter, there is an 88% chance that more than 70% of the GHI amount is received in a direct line. Figure 10 (a) shows an example of a clear day, in which DNI is greater than GHI for all hours. On the other hand, Figure 10 (b) shows an example of a dusty day, in which DHI is almost equal to or a big part of the GHI value because dust causes sunlight to diffuse.

The final number of features is 39, as listed in Table 4. Compared to Table 3, this table includes 10 aerosol variables instead of the 3 variables available in the AERONET dataset. Also, it includes highly correlated lag 2 and lag 3 features with GHI in addition to the previous seven days’ GHI values. These additional features were included because this merged dataset does not have a lot of missing records as the AREONET and K.A.CARE merged dataset.

Records, where GHI is less than 1 (night hour), were eliminated. Also, null records, which resulted from the unavailability of previous days’ GHI values, were deleted.

3.1.3. Data Normalization and Dividing

Min-max scaler was used to normalize data features to the range of [0,1]. Denormalization to the normal range was applied after the training was finished. Data was divided into 70% for training, and 30% for validation and testing. Table 5 shows each dataset used in this work. It clarifies the period covered, the number of missing days, and the total number of hourly records used for training, validation, and testing using the dividing percentages specified earlier. The table also shows the GHI mean, standard deviation (SD), and variance (var) of each data portion in addition to the percentage for three weather conditions (sunny, partly clear, unclear) outa of the total records. Sunny weather here means the DNI value is 80% of the GHI value or higher [48], whereas unclear weather means the DHI value is 90% of GHI or larger due to clouds, haze, or dust [48,49]. Partly clear weather does not meet the aforementioned conditions.

3.2. Models’ Development

In this section, five DL-based models used for GHI forecasting are described. These models are LSTM, GRU, BiLSTM, BiGRU, and LSTM-AE.

3.2.1. LSTM

An LSTM model was developed for next-hour GHI forecasting, as clarified in Figure 11, which comprises three LSTM layers for feature extraction and one dense layer for GHI prediction. Refer to [50,51] to know more about LSTM.

3.2.2. GRU

A GRU model was developed for next-hour GHI forecasting, as clarified in Figure 12, which comprises three GRU layers for feature extraction and one dense layer for GHI prediction. Refer to [52] to know more about GRU.

3.2.3. BiLSTM

A BiLSTM model was developed for next-hour GHI forecasting, as clarified in Figure 13, which comprises two LSTM layers and one BiLSTM layer for feature extraction, followed by one dense layer to make GHI predictions. Refer to [50,53] to know more about BiLSTM.

3.2.4. BiGRU

A BiGRU model was developed for next-hour GHI forecasting, as clarified in Figure 14, which comprises two GRU layers and one BiGRU layer for feature extraction, followed by one dense layer to make GHI predictions. Refer to [54] to know more about BiGRU.

3.2.5. LSTM-AE

An LSTM-AE model was developed for next-hour GHI forecasting, as clarified in Figure 15, in which both the encoder and decoder consist of two LSTM layers, followed by a dense layer to make GHI predictions. Refer to [55,56] to know more about LSTM-AE.

3.3. Implementation

Keras library was used to create DL models. The experiments were performed on a laptop with NVIDIA GeForce RTX 3070 GPU and 16 GB memory. The hyperparameters used are 200 epochs, a batch size of 256, and a learning rate equal to 0.001. For models’ optimization, Adam algorithm was used in addition to a dropout layer with a value of 0.2 and weight decay equal to 0.000001. ReLU activation function was used for all layers and the loss function was the mean squared error (MSE).

3.4. Evaluation Metrics

In this work, four well-known performance evaluation metrics are utilized to evaluate the forecasting models: MAE [57], RMSE [57], MAPE [58], and FS [59]. The equations used to calculate these metrics are given below:

$$\begin{array}{c}MAE=\frac{1}{N}{\displaystyle \sum _{\mathrm{i}=1}^N}\left|{actual GHI}_i-{forecast GHI}_i\right|\end{array}$$

(1)

$$\begin{array}{c}RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum _{\mathrm{i}=1}^N}{({actual GHI}_i-{forecast GHI}_i)}^2}\end{array}$$

(2)

$$\begin{array}{c}MAPE=\frac{1}{N} {\displaystyle \sum _{\mathrm{i}=1}^N}\left|\frac{actual {GHI}_i-{forecast GHI}_{\mathrm{i}}}{actual {GHI}_i}\right|\times 100\%\end{array}$$

(3)

$$\begin{array}{c}FS=1-\frac{{Metric}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{e}\mathrm{d}}}{{Metric }_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}} \ast 100 \%\end{array}$$

(4)

4. Results and Discussion

This section illustrates and discusses the results of the forecasting performance of five DL-based models using three different types of inputs: lagged GHI features (Section 4.1), weather and solar radiation components features (Section 4.2), and aerosol features (Section 4.3). Finally, all the results of experiments using the FS of all the models are presented (Section 4.4).

4.1. Effect of Using Lagged GHI Features on Forecasting

To study this effect, all five forecasting models were trained and tested three times with the same records. In the first experiment, 26 features were used in training, including the GHI values for the previous 3 hours to the forecasting hour and the GHI values for the same forecasting hour on previous days over the last year. The second experiment included 16 features after eliminating 10 highlighted features from the first experiment, which represent the GHI value of the same forecasting hour 90 days up to 1 year ago. In the last experiment, only GHI values for the previous 3 hours to the forecasting hour were used. Table 6 shows the features used in all three experiments.

Figure 16 displays the average RMSE results of 30 runs of five forecasting models for the K.A.CARE dataset on the left and the NSRDB dataset on the right. It is noted that using 26 features (K.A.CARE-1) did not improve the RMSE results compared to using 16 features (K.A.CARE-2). In fact, it increased RMSE slightly by a maximum of 5 points. On the other hand, using 16 features (K.A.CARE-2) compared to using only 6 features (K.A.CARE-3) delivered contradictory results. With the LSTM, BiGRU, and LSTM-AE models, using 6 features enlarged RMSE results, whereas using 16 features instead of 6 improved the RMSE results of the GRU and BiLSTM models by less than 1 point, which could be related to the stochastic nature of DL models training. Looking at the NSRDB dataset results, using 26 features (NSRDB-1) made the RMSE results worse compared to using 16 features (NSRDB-2). Also, the results of experiments 2 and 3 (NSRDB-2 & NSRDB-3) are contradictory. Using 16 features instead of 6 slightly improved the RMSE results of the LSTM and BiGRU models, whereas no improvement was achieved with the remaining three models.

Tt is deduced from the RMSE results of both datasets that using the GHI values for the previous 3 hours is satisfactory and that including the GHI values for the same hour on previous days does not improve forecasting. Additionally, the LSTM-AE model achieves the best RMSE result with 6 features only, which is equal to 50.69 for the K.A.CARE dataset and 63.51 for the NSRDB dataset.

The average MAE results of 30 runs of five forecasting models are shown in Figure 17. The K.A.CARE dataset results on the left and the NSRDB dataset results are on the right. It is noted that using 26 features (K.A.CARE-1) did not improve the MAE results compared to using 16 features (K.A.CARE-2). In fact, it increased MAE slightly by a maximum of 4 points. Also, using 16 features (K.A.CARE-2) compared to using only 6 features (K.A.CARE-3) did not improve the MAE results, except for the BiLSTM model with less than 1 point difference, which could be related to the stochastic nature of the DL models training. On the other hand, the NSRDB dataset results for experiments 1 and 2 (NSRDB-1 vs. NSRDB-2) show that using 26 features made the MAE results worse compared to using 16 features and using 16 features compared to 6 also did not help (NSRDB-2 vs. NSRDB-3). An improvement of less than 1 point in the MAE results of the LSTM and BiLSTM models is mostly related to testing stochastic error.

It is inferred from the MAE results of both datasets that using GHI values for the previous 3 hours is satisfactory and that including GHI values for the same hour on previous days does not improve forecasting. Additionally, the LSTM model achieves the best MAE result with only 6 features, which is equal to 25.8 for the K.A.CARE dataset and 42.04 for the NSRDB dataset.

The average MAPE results of 30 runs of five forecasting models are shown in Figure 18. The K.A.CARE dataset results appear on the left and the NSRDB dataset on the right. It is noted that using 26 features (K.A.CARE-1) did not improve the MAPE results compared to using 16 features (K.A.CARE-2), except for the BiGRU model with 10 points improvement. Unlike the RMSE and MAE results, using 16 features (K.A.CARE-2) compared to using only 6 features (K.A.CARE-3) enhanced the MAPE results significantly, excluding for the BiGRU and LSTM-AE models. MAPE of the LSTM improved from 33.26 to 20.97 and GRU went from 50.51 to 21.69. The best improvement is with the BiLSTM model, which went from 71.16 to 38.03, and the least improvement is with the LSTM-AE model, which went from 21.33 to 19.4. Regarding the NSRDB dataset results, it is seen that using 26 features (NSRDB-1) enhanced the MAPE results slightly for the BiLSTM and LSTM-AE models, whereas the MAPE result improved by almost 7 points for the BiGRU model, going from 52.39 to 45.6. When the results of the NSRDB-2 and NSRDB-3 are compared, it is observed that using 16 features enhanced the MAPE results for all models, except for the BiGRU model. However, since the MAPE results of the BiGRU model improved with 26 features, this can be considered as a stochastic error.

From Figure 18, it is concluded that using more features improves MAPE results, unlike the situation with the RMSE and MAE results. Additionally, the LSTM-AE model achieves the best MAPE result with 26 features, which is equal to 20.21 for the K.A.CARE dataset and 29.51 for the NSRDB dataset.

4.2. Effect of Using Weather and Solar Radiation Components Features on Forecasting

To study this effect, all five models were trained and tested twice with the same records. First, training was conducted with 50 features (experiment 1), including GHI with other weather-lagged features, and then again with GHI-lagged features only (experiment 2) without the highlighted features, as clarified in Table 7. Both experiments were conducted twice using the K.A.CARE dataset and the NSRDB dataset.

Figure 19 illustrates the average RMSE results of 30 runs of five forecasting models for the K.A.CARE dataset on the left and the NSRDB dataset on the right. It is noted that using 50 features, including the results of weather and solar radiation components’ lagged features (K.A.CARE-1), is slightly better than using GHI-lagged features only (K.A.CARE-2). The difference ranges from 0.77 with BiLSTM to almost 6 points with LSTM-AE. It is also noted that LSTM-AE achieves the best RMSE result, equal to 52.94 with 50 features. On the other hand, using weather and solar radiation components’ lagged features with the second dataset (NSRDB-1) slightly improved the RMSE results of three models only: BiLSTM, BiGRU, and LSTM_AE, whereas using 26 features (NSRDB-2) delivered better RMSE results for the LSTM and GRU models. Again, the LSTM-AE model achieved the best RMSE result equal to 66.06 with 50 features.

From Figure 19, it is concluded that using weather in addition to solar radiation components’ lagged features improved the RMSE results slightly. However, this slight improvement might not be worth the loss in efficiency due to the increase in the number of parameters. Besides, the LSTM-AE model achieves the best RMSE results with 50 features.

Figure 20 presents the average MAE results of 30 runs of five forecasting models for the K.A.CARE dataset on the left and the NSRDB dataset on the right. It is noted that using 50 features, including weather and solar radiation components’ lagged features (K.A.CARE-1), is slightly better than using GHI-lagged features only (K.A.CARE-2). The difference ranges from 1 with the BiLSTM model to almost 5 points with the GRU model. It is also noted that the GRU model achieves the best MAE result, equal to 28.23 with 50 features. On the other hand, using weather and solar radiation components’ lagged features with the second dataset (NSRDB-1) slightly improves the MAE results by less than 2 points. It is also noted that the LSTM model achieves the best MAE result, equal to almost 44 with 26 features (NSRDB-2).

It is concluded from the MAE results of both datasets that using weather in addition to solar radiation components’ lagged features improved the MAE results slightly. However, this slight improvement might not be worth the loss in efficiency due to the increase in the number of parameters.

The average MAPE results of 30 runs of five forecasting models are shown in Figure 21. The K.A.CARE dataset results appear on the left and the NSRDB dataset on the right. From the figure, it is seen that using 50 features, including weather and solar radiation components’ lagged features (K.A.CARE-1), improved the MAPE results of (K.A.CARE-2) slightly for the LSTM, BiGRU, and LSTM-AE models. However, the GRU result did not improve, and the BiLSTM model result became worse by almost 10 points. It is also noted that the LSTM-AE model achieves the best MAPE result of around 20 in both experiments. On the other hand, using weather and solar radiation components’ lagged features with the second dataset (NSRDB-1) did not improve the MAPE results. Also, the LSTM-AE model achieves the best MAPE result of around 29 in both experiments.

It is concluded from the MAPE results of both datasets that using weather in addition to solar radiation components’ lagged features might improve the MAPE results slightly. However, this slight improvement might not be worth the loss in efficiency due to the increase in the number of parameters.

4.3. Effect of Using Aerosol Features on Forecasting

To examine this effect, all five models were trained and tested twice with the same records using the K.A.CARE and AERONET merged dataset, once with 35 features (experiment 1), including 3 aerosol features, and once without the highlighted aerosol features (experiment 2), as clarified in Table 8.

Another two experiments were conducted using the NSRDB and GIOVANNI merged dataset, once with 39 features (experiment 1), including 10 aerosol features, and once without the highlighted aerosol features (experiment 2), as clarified in Table 9.

Figure 22 displays the average RMSE results of 30 runs of five forecasting models for the K.A.CARE and AERONET merged dataset on the left and the NSRDB and GIOVANNI merged dataset on the right. It is noted that using aerosol features (K.A.CARE & AERONET-1) slightly enhanced the RMSE results for the BiLSTM and LSTM-AE models only and worsened the results for other models. On the other hand, using aerosol features with the second dataset (NSRDB & GIOVANNI-1) delivered results very similar to those of experiment 2 in which aerosol features were eliminated (NSRDB & GIOVANNI-2).

From both datasets, it is concluded that using aerosol features might not improve RMSE at all and make a slight improvement at best. Also, the LSTM-AE model achieves the best RMSE with aerosol features, which is equal to 46.19 with the first dataset and 62.46 with the second dataset.

The average MAE results of 30 runs of five forecasting models are shown in Figure 23. The results of K.A.CARE and AERONET merged dataset appear on the left and the results of NSRDB and GIOVANNI merged dataset on the right. It is noted that using aerosol features (K.A.CARE & AERONET-1) slightly enhanced the MAE results for the GRU, BiLSTM, and LSTM-AE models. Also, the LSTM-AE model achieves the best MAE with aerosol features, which is 25.69. On the other hand, using aerosol features with the second dataset (NSRDB & GIOVANNI-1) slightly improved the MAE results for all models except the BiLSTM. Also, the LSTM and GRU models achieve the best MAE value of around 37. From both datasets, it is concluded that using aerosol features might slightly improve the MAE.

Figure 24 displays the average MAPE results of 30 runs of five forecasting models for the K.A.CARE and AERONET merged dataset on the left and the NSRDB and GIOVANNI merged dataset on the right. It is noted that using aerosol features (K.A.CARE & AERONET-1) slightly enhanced the MAPE results for the GRU, BiLSTM, and LSTM-AE models. Also, the GRU and LSTM-AE models achieved the best MAPE with aerosol features, which is around 8. On the other hand, using aerosol features with the second dataset (NSRDB & GIOVANNI-1) slightly enhanced the MAPE results for the LSTM, BiLSTM, and LSTM-AE models. Also, the LSTM model achieved the best MAPE value of around 24. From both datasets, it is concluded that using aerosol features might slightly improve the MAPE, but this improvement might not pay off due to the loss in efficiency.

4.4. FS of all Models

In this section, FS results of all the models are presented. FS measures the enhancement in forecasting compared to the persistence method. This metric helps to evaluate the models’ performance in comparison to other models developed using different datasets.

Figure 25 shows the FS of the LSTM model for three metrics: RMSE, MAE, and MAPE. The results for all the experiments discussed in Section 4.1, 4.2, and 4.3 are presented in this figure. The name of each dataset along with the number of features used are clarified, for example, K.A.CARE-50 means the K.A.CARE dataset with 50 features. From Figure 25, it is noted that the FS _RMSE of the LSTM model ranges from 55% with NSRDB-50 to 66% with K.A.CARE & AERONET-16. Also, the FS _MAE ranges from 67% with NSRDB-50 to 80% with K.A.CARE-6, whereas the FS _MAPE ranges from 72% with K.A.CARE and AERONET-19 to 84% with K.A.CARE-50 and K.A.CARE-16. The most significant improvement is seen in the MAPE metric first, followed by the MAE metric.

Figure 26 displays the FS of the GRU model for three metrics: RMSE, MAE, and MAPE. The results for all experiments discussed in Section 4.1, 4.2, and 4.3 are presented in this figure. The name of each dataset along with the number of features used is clarified. For example, K.A.CARE-16 means the K.A.CARE dataset with 16 features. From Figure 26, the FS _RMSE of the GRU model ranges from 55% with NSRDB-50 to 65% with K.A.CARE-16, K.A.CARE-6, and K.A.CARE and AERONET-16. Also, the FS _MAE ranges from 65% with NSRDB -50, NSRDB-26, and NSRDB -16 to 79% with K.A.CARE-6, whereas the FS _MAPE ranges from 61% with K.A.CARE-6 to 83% with K.A.CARE-50, K.A.CARE-26, and K.A.CARE-16. The most significant improvement is in the MAPE metric first, followed by the MAE metric.

Figure 27 shows the FS of the BiLSTM model for three metrics: RMSE, MAE, and MAPE. The results for all experiments discussed in Section 4.1, 4.2, and 4.3 are presented in this figure. The name of each dataset along with the number of features used is clarified. For example, NSRDB-50 means NSRDB dataset with 50 features. From Figure 27, the FS _RMSE of the BiLSTM model ranges from 58% with NSRDB-26 to 65% with K.A.CARE and AERONET-19. Also, the FS _MAE ranges from 66% with NSRDB-50 and NSRDB-26 to 76% with K.A.CARE and AERONET-19, whereas the FS _MAPE ranges from 44% with K.A.CARE-6 to 71% with K.A.CARE and AERONET-19. The most significant improvement is in the MAE metric first, followed by the MAPE metric.

Figure 28 shows the FS of the BiGRU model for three metrics: RMSE, MAE, and MAPE. The results for all the experiments discussed in Section 4.1, 4.2, and 4.3 are presented in this figure. The name of each dataset along with the number of features used are clarified. For example, NSRDB-6 means the NSRDB dataset with 6 features. From Figure 28, the FS _RMSE of the BiGRU model ranges from 55% with NSRDB-26 to 65% with K.A.CARE and AERONET-16. Also, the FS _MAE ranges from 61% with NSRDB-26 to 75% with K.A.CARE-6, and K.A.CARE and AERONET-16, whereas the FS _MAPE ranges from 48% with K.A.CARE-16 to 68% with K.A.CARE and AERONET-16. The most significant improvement is in the MAE metric first, followed by the MAPE metric.

Figure 29 shows the FS of the LSTM-AE model for three metrics: RMSE, MAE, and MAPE. The results for all experiments discussed in Section 4.1, 4.2, and 4.3 are presented in this figure. The name of each dataset along with the number of features used are clarified, for example, K.A.CARE and AERONE-16 means the K.A.CARE and AERONET merged dataset with 16 features. From Figure 29, the FS _RMSE of the LSTM-AE model ranges from 59% with NSRDB-26 and NSRDB-16 to 67% with K.A.CARE and AERONET-19. Also, the FS _MAE ranges from 67% with NSRDB-26 to 79% with K.A.CARE and AERONET-19, whereas the FS _MAPE ranges from 71% with K.A.CARE and AERONET-16 to 85% with K.A.CARE-50 and K.A.CARE-16.The most significant improvement is seen in the MAPE metric first, followed by the MAE metric.

From all the figures in this section, it is concluded that the FS results of the five models are better than the persistence method by at least 44% and at most 85%. The best FS results are realized by the LSTM-AE model. The most improved metric is the FS _MAPE for the LSTM, GRU, and LSTM-AE models, whereas it is the FS_MAE for the BiGRU and BiLSTM models. The least improved metric is the FS_RMSE for the LSTM, GRU, and LSTM-AE models, whereas it is FS_MAPE for the BiGRU and BiLSTM models.

5. Conclusion

This paper aims to evaluate the ability of five DL-based models to provide accurate forecasting of next-hour GHI using different sets of features, including weather variables, aerosol variables, along solar radiation components. Fourteen experiments were conducted in which LSTM, GRU, BiLSTM, BiGRU, and LSTM-AE models were tested using different feature sets to investigate their effect on forecasting results. The findings are summarized as follows:

• Although the GHI values of the same forecasting hour on previous days have a stronger or equal correlation with the output than the GHI values of the previous three hours on the same day, using the latter in forecasting provides better accuracy, especially if measured by RMSE or MAE. However, only MAPE was improved when the GHI values of the same forecasting hour on previous days were used for prediction. Therefore, the decision about the inclusion of more GHI-lagged features depends on the performance metric of interest and the size of the dataset.
• Using weather, aerosol, and solar radiation components’ lagged features improves RMSE, MAE, and MAPE results slightly. However, this slight improvement might not be worth the loss in efficiency due to the increase in the number of parameters. Therefore, the decision about the inclusion of these features depends on a tradeoff between performance and efficiency.
• The LSTM-AE model provides the best forecasting results with all feature sets, followed by the LSTM and GRU models, whereas the BiLSTM and BiGRU models provide the worst.
• The best forecast skills results are achieved by the LSTM-AE model, which reaches 85%.
• FS _MAPE is the most improved metric for the LSTM, GRU, and LSTM-AE models, whereas it is FS_MAE for the BiGRU and BiLSTM models.
• The best RMSE, MAE, and MAPE results are 46.19, 25.69, and 8.18 achieved by the LSTM-AE model with the K.A.CARE and AERONET merged dataset with 19 features.
• Regarding datasets, all results associated with the NSRDB dataset are worse than results associated with the K.A.CARE dataset. Ground-based measurements are more accurate than satellite-based observations and thus provide better forecasting. However, ground-based data suffer from a huge number of missing values due to device malfunction or maintenance scheduling. It is safe to use satellite data for model development purpose and assume that results would be better with ground-based data.