Fruit fly is a kind of flying insect, which is very sensitive to the external environment because of its superior olfactory and vision. Firstly, the olfactory organ is used to obtain the odor floating in the air. Then, it will distinguish the general direction of the food source and fly to the source of food. Finally, the fruit fly can discover the position of food by its keen vision, and then fly to the position. The process of searching food for the fruit flies can be simulated as follows (Mitić et al., 2015):

1) Randomly initialize the population size, maximal number of iterations, and position coordinates ( x , y ) of the group in a set interval.

In the ordinary FOA, the individual fruit fly seeks the food source with the pre-set step size. Obviously, if the step size is too small, the search space will be limited, and it will cause the problem of local optimum. On the contrary, if the step size is too large, its local search ability will become weaker, and the convergence rate will slow down. To deal with these issues, the setting of step size should adhere to the following principles. In the initial phase of iterations, the step size should be large to ensure global optimization performance. On the contrary, in the later stage, the step size should be small to ensure local search performance.

Therefore, there are a few successful algorithms for the improvement of step size of fruit flies, such as the decreasing step fruit fly optimization algorithm (DSFOA) (Hu et al., 2017), self-adaptive step fruit fly optimization algorithm (FFOA) (Yu et al., 2016), and improved fruit fly optimization algorithm (IFFO) (Pan et al., 2014). In the DSFOA and IFFO, the step size decreased quickly in the initial phase of iteration, which cannot guarantee the global optimization performance of the algorithm. In this paper, we propose the multivariate adaptive step size, which can be demonstrated as follows: L i = L 0 ⋅ exp [ − N ( G i G max ) α ] , (8) where L 0 is the initial step size, G i is the current number of iterations, G max is the maximum number of iterations, N is a positive integer, and the exponential factor α is a constant within ( 0,10 ) . The positive integer N and the exponential factor α control the decreasing rate of step size and realize the better local search performance. In order to choose proper values of N and α , the convergence ability of algorithm under different parameter values is compared. The initial step size L 0 is set to 20, and the maximum number of iterations G max is set to 100. Figure 1 gives the variations of step size in Eq. 8 corresponding to different values of α when N = 10 , N = 15 , N = 20 , and N = 25 , respectively.

As shown in Figure 1, the step size decreases gradually from 20 to 0 with the increasing iteration number and different values of α correspond to different step size change trends. In the initial stage of iterations, the algorithm has the largest step size, which can guarantee the global optimum. As the iteration number increases, the capability of local search is gradually enhanced to find the local optimum value, which can be seen from the rapid decline in the curves. Therefore, the dynamic step size can realize the balance of global search capability and local optimization ability.

Besides, from the subfigures in Figure 1, the step size changes relatively symmetrical when α = 3 , α = 5 , and α = 7 . When α = 1 , the curve drops sharply from the beginning, which means the step size will become small even before achieving the global optimum, and the step size cannot achieve 20 at the beginning of iteration. The moment when step size begins to decline is a bit later when α = 9 . There seems to be no difference in convergence performance when N takes different values. So in Empirical Analysis, we will test the performance of the proposed model with different values of N and α to search for the optimal value, and we will substitute the optimal N and α into the model for short-term load forecasting.

The GRNN is a kind of neural network using the radial basis function and has been very popular in applications in recent years. It can establish the implicit mapping relationship according to the sample data, so that the output can converge the optimal regression surface. Once the sample is determined, the only goal is the determination of smoothing parameter in the kernel function (Ozturk and Turan, 2012; Kumar and Malik, 2016).

Assuming that f ( x , y ) is the joint probability density function of random variable X and variable Y , the observed value of X is x 0 , and the regression of Y with respect to X is Y ^ ( x 0 ) = ∫ − ∞ ∞ y f ( x , y ) d y ∫ − ∞ ∞ f ( x 0 , y ) d y . (9) Based on the Parzen non-parametric estimation, the density function f ( x 0 , y ) can be estimated by the sampled dataset { x i , y i } i = 1 n : f ( x 0 , y ) = 1 n ( 2 π ) p + 1 2 σ 1 σ 2 … σ p σ y ∑ i = 1 n e − d ( x 0 , x i ) e − d ( y , y i ) , (10) d ( x 0 , x i ) = ∑ j = 1 n [ ( x 0 j − x i j ) / σ j ] 2 , (11) d ( y , y i ) = ( y − y i ) 2 , (12) where n is the sample size, p is the dimension of random variable X , and σ is the width coefficient of the Gaussian function, which is called the smoothing parameter.

Substituting Eq. 10 into Eq. 9 yields Y ^ ( x 0 ) = ∑ i = 1 n ( e − d ( x 0 , x i ) ∫ − ∞ ∞ y e − d ( y , y 0 ) d y ) ∑ i = 1 n ( e − d ( x 0 , x i ) ∫ − ∞ ∞ e − d ( y , y 0 ) d y ) . (13) Note that ∫ − ∞ ∞ z e − z 2 d z = 0 , (Eq. 13) can be simplified as follows: Y ^ ( x 0 ) = ∑ i = 1 n y i e − d ( x 0 , x i ) ∑ i = 1 n e − d ( x 0 , x i ) . (14) The predicted value in Eq. 14 is the weighted sum of the observations of the dependent variable, and the weights are e − d ( x 0 , x i ) . The GRNN is composed of input layers, pattern layers, summation layers, and output layers. Once the learning samples are determined, the structure of neural network and the connection weights between neurons are completely determined. Therefore, the GRNN does not need to adjust the connection weight values between neurons, but to adjust the transfer function of each unit by changing the smoothing parameter to obtain the best regression result, which is different from the traditional error backward propagation algorithm. Thus, a key step in the GRNN is to determine the value of the smoothing parameter.

In this paper, the MAFOA is applied to optimize the smoothing parameters in the GRNN. The MAFOA-GRNN takes the root mean square error (RMSE) of GRNN as the fitness function of MAFOA, so as to calculate the smell concentration in each iteration. Part of the training data are used in the MAFOA to select the best parameters for the GRNN. When the algorithm reaches the maximum number of iterations, the location of the fruit fly with best smell concentration is obtained. Then, these optimal parameters will be used in the GRNN to get the optimal prediction model. The flowchart of the MAFOA-GRNN model is shown in Figure 2.

The power load data used in this paper are hourly and obtained from a power grid in Wuhan with 2,880 observations ranging from January 1, 2014, to April 30, 2014, which are shown in Figure 3. In this section, we predict the power load of the last day of each month. The in-sample data are power load data of each month except the last day, and the out-of-sample data are the power load data of the last day of each month.

As shown in Figure 3, the short-term power load has obvious periodicity. Therefore, historical load data are an important reference for forecasting. In order to accurately predict the power load, the factors influencing the power load should be considered as much as possible. The factors related to load forecasting include date classification (weekday, weekend, holiday), daily temperature (maximum, minimum, average temperature), and weather condition.

Combining the influence factor, the improved GRNN adopts a three-layer network structure. The input variables of the GRNN are shown in Table 1, and the corresponding output vector is the power load value at t o’clock on day d .

The original load data are normalized to eliminate the impact of the dimensions between indicators. In addition, the input variables of GRNN in Table 1 should be numerical data, so we quantify the above weather factors and date type factors. Meanwhile, we propose an efficient interval partition technique to handle temperature and weather types:

1) Normalization of load data. All load data are normalized by using the linear transformation method, given by

This paper uses the normalized root mean square error (NRMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE) as the evaluation criteria, given by N R M S E = 100 y ¯ 1 N ∑ i = 1 N ( y i − y ^ i ) 2 , (16) M A E = 1 N ∑ i = 1 N | y i − y ^ i | , (17) M A P E = 1 N ∑ i = 1 N | y i − y ^ i y i | , (18) where y ¯ is the mean of value, y ^ i is the predicted value, y i is the observation value, and N is the number of data.

Although the NRMSE, MAE, and MAPE can be used as criteria to obtain model predicted loss values, it cannot be verified whether the comparison result is statistically significant. To solve this problem, Diebold and Mariano (1994) proposed the Diebold–Mariano (DM) test to test the statistical significance of different prediction models. Assume that model B and model T do the forecasting task in period t at the same time, and we wonder if there are significant differences in the performance between the two models. The original hypothesis is that the forecast accuracy for two models is the same, which is equivalent to the mean value of relative loss function of 0. The DM statistics is defined as follows: D M = d ¯ 2 π f d ( 0 ) T , (19) where d ¯ = 1 T ∑ t = 1 T d t (20) is the sample mean loss differential, in which d t = L o s s T − L o s s B (21) is the relative loss function, where L o s s T and L o s s B are the loss function of predicted errors of test model T and benchmark method B at time t, respectively.

Note that, in this paper, the mean-squared prediction error (MSPE) is used as the loss function: L o s s i = 1 N ∑ i = 1 N ( y t − y ^ i t ) 2 , (22) where y ^ i t is the predicted value of model i at time t. f d ( 0 ) = 1 2 π ∑ τ = − ∞ ∞ γ d ( τ ) (23) is the spectral density of relative loss function at frequency zero. γ d ( τ ) = E [ ( d t − μ ) ( d t − τ − μ ) ] (24) is the autocovariance of d t at displacement τ , where μ is the population mean loss differential.

If the p -value corresponding to D M is less than the significant level, which normally is 0.01 or 0.05, the original hypothesis is rejected; otherwise, it cannot be rejected.

To determine the values of parameters α and N , we apply α = 1,3,5,7,9 and N = 10,15,20,25 into the model to test the performance of the model. In this section, the data from January 1, 2014, to January 30, 2014, in Wuhan are used as training data, and the load data on January 31, 2014, are regarded as test data. Finally, the anti-normalization processing is carried out, and the NRMSE, MAE, and MAPE are calculated.

Table 3 shows the prediction errors for different α and N . We can see that whatever value N takes, three types of errors are obviously higher than others when α = 1 . When α = 3 and N = 25 , the NRMSE and MAE are small. The NRMSE, MAE, and MAPE are small at the same time when α = 5 and N = 20 , which is consistent with that reported in the FFOA proposed by Yu et al. (2016). When α = 9 and N = 20 , the NRMSE and MAPE are both small. But the smallest values of NRMSE, MAE, and MAPE are obtained when α = 7 and N = 15 , which means the forecasting performance is the best at this moment. So, we can initially claim that when α = 7 and N = 15 , the step size can realize the balance of global search capability and local optimization ability. Figure 4 shows the fitness curves of MAFOA when α = 7 and N = 15 .

Figure 4 provides the fitness variation with the increase of iteration number. Figure 4A shows the convergence situation of different N when α = 7 , and Figure 4B shows the convergence performance of different α when N = 15 . It can be seen that, under the condition of α = 7 , when we choose N = 15 , the algorithm has the minimal fitness value and arrives at its optimal value much more quickly than other conditions; under the condition of N = 15 , when we choose α = 7 , the algorithm has the same performance. Accordingly, α = 7 and N = 15 are perceived as an ideal choice in the step size formula.

After choosing α = 7 and N = 15 , we train the power load data of each month except the last day to predict the load of the last day of each month. Table 4 shows the prediction results obtained by the MAFOA-GRNN algorithm from January to April 2014, respectively. The relative errors are basically within 2%, and the accuracy is high.

In order to test the forecasting performance of the proposed model, the backpropagation (BP) neural network, support vector machine (SVM), GRNN, PSO-GRNN, FOA-GRNN, and DSFOA-GRNN are regarded as benchmark models to be compared with MAFOA-GRNN in short-term power load forecasting. The PSO was proposed by Kennedy and Eberhart in 1995, which was inspired by the swarm behavior of birds. The FOA proposed by Pan in 2012 was also used in this work. Since PSO and FOA are both classical optimization algorithms that have been widely utilized in research, we have chosen PSO-GRNN and FOA-GRNN as benchmark models. The DSFOA proposed by Hu et al. in 2017 is an improvement algorithm of FOA. With the decreasing step size in mind, the DSFOA performed well in optimizing the spread parameter of GRNN. The flight distance is updated referring to the sigmoid function. So, DSFOA-GRNN has also been compared with our proposed model. Besides, some other basic prediction models are also taken into account, such as the BP neural network and SVM. Figure 5 shows the relative error curves of the single models on January 31. Figure 6 shows the relative error curves of the hybrid models on January 31.

It can be seen from Figure 5 that, in the commonly applied forecasting methods, the GRNN has the best prediction ability. Figure 6 shows that the proposed method can accurately predict the overall trend of power load, and the fitting effect is very good. From the relative error curves, it can be seen that MAFOA-GRNN can offer a better predicting performance and higher precision than DSFOA-GRNN, FOA-GRNN, and GRNN. In addition, the relative errors of MAFOA-GRNN are more stable, and the majority are below 0.02, which demonstrates that the improved FOA is perceived as an ideal method in optimizing model parameters during GRNN training.

Then, the anti-normalization processing is carried out, and the comparison results of NRMSE, MAE, and MAPE evaluation criteria are shown in Figures 7–9. Table 5 shows the error analysis of the training set and test set.

Obviously, MAFOA-GRNN has the smallest NRMSE, MAE, and MAPE, followed by FOA-GRNN, but the BP neural network has the worst performance. Besides, the prediction error of the training set and test set has no obvious difference, which indicates that MAFOA-GRNN has high generalization performance. According to the comparison results, it can be concluded that MAFOA-GRNN outperforms other models in both accuracy and stability. Table 5 demonstrates the same conclusions as above.

Although the NRMSE, MAE, and MAPE can be used as criteria to obtain model-predicted loss values, it cannot be verified whether the comparison result is statistically significant. To statistically compare the differences between the prediction accuracy of different models, the DM statistics test is carried out in this paper, and the results are shown in Table 6. For all the benchmark models, the values of the MAFOA-GRNN model proposed in this paper are below 0.05, which indicates that the predictive ability of the MAFOA-GRNN model is better than that of DSFOA-GRNN, DSFOA-GRNN, GRNN, SVM, and BP neural network under the confidence interval of 95%.

According to the above comparisons, the following three main conclusions can be summarized:

1) The proposed MAFOA-GRNN outperforms the GRNN, which indicates that the MAFOA can optimize the smoothing parameter of GRNN effectively.