Artificial intelligence enabled microgrid power generation prediction

3.1 System model

To have a better understanding of the PV electricity generation process, figuring out the PV system structure and setting up configurations of PV panels are very important.

Sandia National Laboratories, which operates under the U.S. Department of Energy has published one of the related mathematical models is: the Plane of Array (POA) Model [30], which figures out the mathematical relation between the solar energy that PV panels absorbed, and the solar radiation is necessary. POA represents the PV panel surface and the irradiance cast on POA can be calculated by equation (1) [30].

in which three main components of the POA irradiance, E b is the POA beam component, E g is the POA ground reflect component, E d is POA sky diffuse component [30]. The overall solar irradiance reflected on the solar panel is the summation of the irradiance from direct sun irradiance, irradiance reflected from the ground, and irradiance diffused from the sky.

Moreover, the POA beam component is decided by Direct Normal Irradiance and the angle between the sun rays and the PV panel, which is determined by not only the solar azimuth, and zenith angles, but also the tilt, azimuth angles of the PV panel [31], in which DNI denotes direct normal irradiance, AOI denotes angle of incidence.

POA formula’s second component ground reflected irradiance can be calculated as the following equation. The main affected variables include global horizontal irradiance GHI, the ground surface reflectivity, which is also called ground albedo ρ , and the tilt angle of the PV panel mentioned in the PV structure section θ T [32]. Here, ρ denotes albedo, which is highly dependent on the ground color and material. When the surface is dark, ρ ≈ 0 ; when the surface is bright white, ρ ≈ 1 [33].

Sky diffuse irradiance has several different theory models, like isotropic sky diffuse model [34], Hay and Davies Sky Diffuse Model [35], Reindl Sky Diffuse Model [36], and Perez Sky Diffuse Model [37]. This report used the simple isotropic model to demonstrate sky diffuse [34,38].

Also, the ground horizontal irradiance (GHI) could be calculated by the diffuse horizontal irradiance (DHI), direct normal irradiance (DNI), and solar zenith angles θ Z following the equation (5)  [39]. In [40], the author mentioned that the DHI value is around 10–20% of GHI value on a sunny day, however, when encountering a cloudy day, the DHI value is almost equal to the GHI value.

As a result, the irradiance reflected on the PV panel in total could be calculated by putting equations (2), (3), (4) into equation (1)

in which the DNI can be derived from an absolute cavity radiometer; POA can be obtained by a pyranometer; AOI is mainly determined by solar azimuth θ A , solar zenith θ Z , and surface tilt angle θ T . The albedo parameter ρ as mentioned before, denotes the reflectivity of the ground surface [31–33]. In practice, parameters θ A , and θ Z are related to the sun’s position which is changed by date (time-related features). The tilt angle of the surface of PV panel θ T will be dynamic as well if it is a solar tracking rack rather than a steady PV panel.

To provide an intuitive feeling of three components: beam, ground reflect, and sky diffuse, the numeric values were provided to understand each component’s contributions to the POA model. The daily average global, beam, and diffuse irradiance component measurements in Kimberly, Idaho are 413, 481, and 132 W/m 2 , respectively [41]. In another research [42], the author provided DNI, GHI, and DHI measurement values in Doha, Qatar, which are 200.4, 225.2, and 94.7 W/m 2 , respectively, among half-year records.

To sum up, despite solar irradiance related weather factors DNI and DHI, there are some time-related features ( θ A , θ Z ), some weather-related features (cloud cover, temperature), and other external environment features ( ρ , θ T ) that existed in the PV generation model.

3.2 Predictive framework

The PV generation forecasting process will follow the steps in the flow chart. After obtaining the original PV generation data and related weather data, two datasets will first be preprocessed to get rid of the missing and wrong data and aligned by their timestamps afterward. Then, the processed time, PV generation data, and weather data will be fed into the XGBoost model to decide on the most important and effective input features for the proposed resilient model. Finally, SARIMA and long short-term memory network (LSTM) algorithms are used predict PV generation in the short term, comparing models’ performance under different input features. The predictive framework flowchart is illustrated in Figure 1.

Figure 1

PV power generation prediction framework flowchart.

3.3 XGBoost model

CART tree algorithm, which is short for classical classification and regression tree, is the title for both trees: Classification tree (the prediction results are types) and regression tree (the prediction results are numeral). CART tree is a nonparametric decision tree algorithm [43]. CART algorithm first puts the input dataset in the root node, splits sub-nodes from the root node according to the attribute, and decides the best homogeneity for the threshold; the splitting process keeps going until a pure subset or meets the maximum node depth. The final node called leaf node is the one who holds the decision [8,43]. This whole iteration process will provide the relative best-fit model.

The eXtreme gradient boosting (XGBoost) is a supervised, improved gradient-boosted trees algorithm, which is integrated by many CART trees. The XGBoost’s regularized objective function L ( ϕ ) includes training loss and regularization term [44].

The differentiable convex training loss ( L ) measures the predictive ability of the model by comparing the difference between prediction value y i ˆ and measurement value y i . f k denotes independent tree structure and leaf weights.

Regularization term ( Ω ) refers to the part that controls the model complexity that prevents the model from overfitting, which can be demonstrated as follows [44,45]:

where Ω defines the complexity of the tree f ; λ determines the strength of the regularization; γ is a penalise nodes constant, when it is greater than 1; and T denotes the number of leaves in the tree.

The XGBoost algorithm optimizes the function by keep adding new trees to simulate the residuals from the last prediction rather than using methods in Euclidean space [43,44]. The purpose of the model is to find f p value to optimise the objective function by searching the smallest score [44]. In dataset D = ( x i , y i ) , ( x i denotes examples, and y i denotes features), the objective function L ( p ) can be presented as as follows [43,44]:

where y ˆ i ( p ) refers to the prediction value on p th iteration of i th instance [44].

To simplify the objective function, the constant terms ∑ i = 1 n l ( y i , y ˆ ( p − 1 ) ) could be removed [44]. By applying one of the common training loss functions mean squared error (MSE) and regularization term Ω ( f p ) , the objective function can be simplified as follows [43]:

After the process of second-order approximation, this objective equation (9) sums the first and second input gradient statistics g i and h i up together as G j and H j and then calculates the score according to the formula, which indicates the quality of the tree structure [43].

This article used the XGBoost model purely for time-series PV data predicting, which separates the various time factors (like months, weeks, days, and hours) from PV generation data and stores them in a tree. Finally, predict future PV generation based on the summarized time features relationship.

3.4 Seasonal ARIMA model

Autoregressive integrated moving average (ARIMA) is one of the effective univariate time series algorithms. As an extended version, seasonal autoregressive integrated moving average (SARIMA) algorithm supports both autoregressive and moving average functions [46], which means SARIMA would identify seasonal changing input data and make better predictions compared to ARIMA. Seasonal changing data refer to the training data value changes due to seasonal factors. Accordingly, we use the SARIMA model to train and predict the PV generation value. The SARIMA model can be formed as SARIMA ( ( p , d , q ) , ( P , D , Q , m ) ) in which ( p , d , q ) represents the nonseasonal feature from the data; ( P , D , Q , m ) represents the seasonal feature. In more detail, p represents trend auto-regression order; d represents trend difference order; q moving average order; P represents seasonal auto-regressive order; D represents seasonal difference order; Q represents seasonal moving average order; and m represents the number of time steps for a seasonal period [46]. The SARIMA model can be formulated as equation (11)

where y t denotes the value of the time series at time t , ϕ denotes the coefficients of the auto-regressive, θ denotes coefficients of the auto regressive forecast errors; ε denotes the moving average forecast error; η denotes coefficients of seasonal forecast errors; and ω denotes coefficients of seasonal auto regressive. The enumeration function was used to list all ( p , d , q ) and ( P , D , Q , m ) combinations to look for the best model fit, which is judged by the lowest Akaike information criterion (AIC) value [47].

3.5 Long short-term memory network (LSTM)

LSTM is one type of the RNN (recurrent neural network) especially overcome exploding and vanishing gradient problems during long-term dependencies, which has recurrent neurons to process the input data through the activation and formulate an output to the next neuron [48]. The LSTM is good at solving the sequence problems because of the feedforward network from the last training [49], which no longer suffers from Simple Recurrent Networks [50]. A typical vanilla LSTM model architecture can be demonstrated as multiple memory blocks shown in Figure 2 [48,51].

Figure 2

LSTM memory block demonstration.

In this special memory block, there are three essential multiplicative units: input gate, output gate, and forget gate, which always use sigmoid as their nonlinear activation function [48]. This memory block first takes previous memory c t − 1 , output h t − 1 , and input signal X t as input to feed through input gate [48,52–54].

Second, in forget gate, the memory block decides what information will be forgotten based on the input signal X t , memory c t − 1 , and output h t − 1 from the previous block [48,52–54].

Third step is generating output combining input signal X t , memory c t − 1 , and output h t − 1 from last iteration [48,52–54].

This memory block iteration will carry on when more input data has been fed through the LSTM model. LSTM networks there would be another loop in the model, which will provide the feedback value as an input vector from an output of the network to the input of the network [48].

4.1 Data preparation

The quality of the input dataset may affect the training model’s performance and accuracy [55], which means two data preprocessing steps: data cleaning and filtering missing data are necessary [56]. After the data preprocessing steps, the data normalized step is required to reduce the noise and normalize the dataset.

In this work, for the XGBoost, SARIMA, and LSTM algorithms, 1 year London area’s PV electricity generation record from Sheffield open-source PV live data (London area) and weather data were from MIDAS UK open weather data^[1]. The data were from the Heathrow station, and both of the datasets (PV and weather) start from 2020/01/01 to 2020/12/27, recording in every 60 min. Due to the synchronization failure and misoperations, there are some repeated data or vacant data. In the data preprocessing step, we removed the repeat data according to the date and time, and then left the vacant records (NaN) as noise. As a result, the PV generation dataset has in total of 8657 records after the preprocessing step.

The PV generation data itself are recorded in Megawatts, and hourly generation records can reach up to thousands of megawatts, sometimes even more than 5,000 MW. However, at night, due to lack of solar irradiance, the PV generation equals 0 mostly. The large periodical difference in the dataset will not benefit the learning process, which means that the data normalization process is necessary as well. MinMaxScaler function was used to scale all the PV generation data under [0,1] scale. The weather data air temperature are stored in the degree Celsius; the global solar irradiance mount is stored in kJ/m².

The Pearson correlation coefficient model was used to analyze the associations behind the feature data. It is a normalization evaluation method of the covariance of two features, which reviews the strength and direction of linear correlation between them [57]. On the one hand, the various weather data were studied including time-related features’ correlations with PV electricity generation, which can be illustrated in Figure 3.

Figure 3

Heatmap of PV generation features correlation.

By analyzing the correlation between features from the dataset, we can see that pure time features like hours, date, and month do not have much correlation with hourly PV generation because their correlation values with PV generation are close to 0. However, previous records of global irradiance and air temperature have a very high correlation to the next hour’s PV generation, which means they could be the main features that contribute to the prediction model.

On the other hand, the correlation between history PV generation (24 h before) and current PV generation was also studied. Because the POA model is very complex and contains various times, weather, and other features, we cannot easily have access to all of them. In contrast, the previous PV generation data are the accurate value simulated by all the factors. The correlation heatmap between history and current PV generation value is shown in Figure 4, where P_24 refers to PV generation from 24 h ago, and P_0 refers to the PV generation in the future 1 h, which needs to be predicted.

Figure 4

Heatmap of PV generation features correlation.

From Figure 4, P_1 to P_3 (previous 3 h) and P_21 to P_24 (yesterday 3 h after) PV generation have a high correlation (above 0.6) to the future 1 h PV generation prediction. The other hours of PV generation also have some associations with the future PV generation data.

To sum up, we used a total of 8752 records (from 2020/01/01 0:00 to 2020/12/21 23:00) to train the XGBoost model. The dataset is developed in six columns: month, date, time, air temperature, global irradiance, and PV generation. By applying different feature combinations to the algorithms to receive a high-performance model that requires minimum data. The testing dataset will be from 6 days of records (from 2020/12/22 0:00 to 2020/12/27 23:00).

4.2 XGBoost

To find the ideal prediction model, this work first used the XGBoost model to verify and analyze the hypothesis we concluded from the first Pearson correlation model that solar irradiance and air temperature are the core features to predict short-term PV generation. The XGBoost model predicts future PV generation based on the previous 1 h weather data and the time data without the data normalization process to keep the original data features. As a result, PV generation prediction uses megawatts as its measurement unit. PV generation data’s relationships in the XGBoost model among different features can be illustrated in Figure 5.

Figure 5

XGBoost time features analysis in PV generation data.

XGBoost feature importances function utilizes F score that indicates each feature’s percentage weight over the weights of all features. As Figure 5 shows, in XGBoost, solar irradiance has over 50% importance and air temperature has around 30% importance for this dataset.

For the hyperparameter optimization part, we still used the grid search method to automatically test the best fit of parameters for the XGBoost model. The XGBoost prediction model fit is shown in Figure 6, which is an XGBoost 1 h ahead prediction model using previous hours, solar irradiance, and air temperature as the training dataset.

Figure 6

XGBoost PV generation prediction.

From the figure, we can see that the XGBoost model only using previous one-hour solar irradiance, air temperature, and hours time data shows low prediction accuracy on the summit PV generation of the day. The average R 2 score is around 0.7; the root mean squared error (RMSE) value is around 467 (MW). This model’s performance is not very well for just using previous weather data and time data as its input. As a result, we decided to use previous PV generation data as input instead of using previous solar irradiance and air temperature. Because previous weather data are more related to previous PV generation instead of PV generation in the future. Also, previous PV generation is always accurate and more related to the predicting panel (containing more complex features) compared to the local weather data.

4.3 SARIMA

To confirm the correlation between history and current PV generation, we utilized the time-series algorithm SARIMA to predict the current generation value using the previous generation value analyzed by date and time. The auto-regression part of the SARIMA model measures the dependency by past observations, which will use the previous 8,700 data records as the training dataset. Also, from observation of the PV generation value, the value repeats in a similar tendency loop for 24 h of solar movement, which means the PV generation value can be set in a 24 h (24 records) seasonal period in a SARIMA model. The test data will be 6 days starting from 2020/12/22 to 2020/12/27.

Our goal is to find the best-fit value for SARIMA ( p , d , q ) ( P , D , Q , m ) to optimize the interest metric. We used the hyperparameter optimization method, Grid Search, to generate and test different combinations of variables p , d , q , and then evaluate the model accuracy by comparing their AIC value. AIC value is a common model selection theory, which will gradually minimize the mean squared error of prediction to find the best fit [58]. Normally, the lower AIC value is considered to be a better model fit.

After the grid search process, the best hyperparameters that fit the SARIMA model are ( 1 , 1 , 1 ) ( 1 , 1 , 1 , 24 ) with the lowest AIC score of 262.36. The SARIMA short-term PV prediction can be illustrated in Figure 7.

Figure 7

SARIMA PV generation prediction.

The SARIMA model’s coefficients are listed in Table 1, where ar.L1 and ma.L1 rows denote autoregressive (AR) and moving average (MA) coefficients for the nonseasonal component of the model, respectively. Similarly, “ar.S.L24” and “ma.S.L24” denote the AR and MA coefficients, respectively, the seasonal component of the model, where the seasonality is specified as 24 h. The sigma2 (sigma square) column denotes the variance of residual values; the coef column denotes weights of each feature; the std err column denotes standard error. The significance of each coefficient is assessed using the z -test, with the associated p -values provided in the “ P > ∣ z ∣ ” column. A p -value less than the conventional threshold of 0.05 indicates that the corresponding coefficient is statistically significant, suggesting that it has a significant association with the predicted PV generation.

As Figure 7 shows, the SARIMA model has relatively high accuracy and confidence (average mean absolute error [MAE], is around 0.0097, average RMSE is around 0.0196) on short-term PV generation predicting. The grey shadow in the figure is the confidence bounds of the prediction model, which refers to the uncertainty of each step of prediction. However, we can still see that there are imprecise predictions during the end of the day and the summit generation in a day.

4.4 LSTM

The LSTM model was first used to conduct one-step ahead PV generation prediction as a comparison group with the SARIMA model. We set the look-back window as 24 (last 1-day records) so that the LSTM algorithm will take previous 1-day records to predict the value. The LSTM model contains 50 neurons, mean square error as loss function, and the adam optimizer. LSTM model’s average score ( R 2 ) is around 0.96, and average RMSE is around 0.2. The best LSTM one-step ahead model we conduct can be illustrated in Figure 8.

Figure 8

LSTM PV generation prediction.

LSTM 1 h ahead short-term prediction using 24 h history records has a relatively good performance, even though sometimes it is still limited in predicting the summit PV generation value of the day and some other turning points. Then, to find a resilient and high-performance model, we also added the weather data (solar irradiance and air temperature) and tested other look-back window sizes on the LSTM model to reduce the input feature amount. Different input feature combinations’ requirements and performance are listed in Table 1.

4.5 Result comparison

We researched the time-series PV generation algorithms SARIMA and LSTM, which use previous PV generation data to predict future generations. We used four common evaluation methods R 2 score ( R 2 ), MAE, MSE, and RMSE.

R 2 score, which also means coefficient of determination, describes the accuracy of dependent variable changes according to the prediction of independent variables. In other words, it explains how well the data fit the model and R 2 = 1 represents the perfect fit. y i denotes the true value of the dataset, y ˆ i denotes the predicted value, and y ¯ i denotes the mean value.

MAE represents the average of absolute residuals between the prediction and true value, which can be calculated as follows:

MSE measures the average square errors between predicted and true values. Compared to MAE, it measures the variance of residuals rather than the average of residuals. It can be calculated as follows:

RMSE is the square root of the MSE value, which measures the standard deviation of residuals and can be calculated as follows:

In the algorithm results comparison table, * denotes the dataset with normalization; ‘denotes to the dataset with weather data (solar irradiance and air temperature). The number after LSTM denotes the look-back number. For example, LSTM*’ 5 refers to the LSTM model utilizing the normalization dataset and weather data, using the previous 5 h of PV generation data as a look-back window. For each of the models, we used an average score among 20 times model fitting (Table 2).

From the comparison, we can see that the LSTM model using 24 history PV generation records and weather data does have excellent performance (average R 2 score can reach 0.9563, and average RMSE is around 0.0195). However, various input features reduce the resilience and robustness of the model. When there are not that much PV generation data available or some of them are not accurate because of cyber attacks, the model will not perform or make an accurate prediction in that circumstance. Instead, we studied various LSTM models that require minimum input data to guarantee the model’s resilience. In the table, LSTM models only use the previous two records and three records of PV generation data as input and have relatively high performance (with an average R 2 score of around 0.90 and an average RMSE score is around 0.026).

However, when previous solar irradiance, air temperature, or other weather data are available combined with previous PV generation data, using multivariables LSTM algorithms will have a slightly better performance in short-term prediction (by comparing LSTM*’ models with LSTM* models).