A Hybrid Model for Forecasting Sunspots Time Series Based on Variational Mode Decomposition and Backpropagation Neural Network Improved by Firefly Algorithm

The change of the number of sunspots has a great impact on the Earth's climate, agriculture, communications, natural disasters, and other aspects, so it is very important to predict the number of sunspots. Aiming at the chaotic characteristics of monthly mean of sunspots, a novel hybrid model for forecasting sunspots time-series based on variational mode decomposition (VMD) and backpropagation (BP) neural network improved by firefly algorithm (FA) is proposed. Firstly, a set of intrinsic mode functions (IMFs) are obtained by VMD decomposition of the monthly mean time series of the sunspots. Secondly, the firefly algorithm is introduced to initialize the weights and thresholds of the BP neural network, and a prediction model is established for each IMF. Finally, the predicted values of these components are calculated to obtain the final predict results. Comparing BP model, FA-BP model, EMD-BP model, and VMD-BP model, the simulation results show that the proposed algorithm has higher prediction accuracy and can be used to forecast the time series of sunspots.


Introduction
Sunspot is the most basic and obvious activity in the solar activity, and the sunspot numbers are an indicator of the total solar activity level. Time-series data of sunspot are considered to be nonlinear, nonstationary, and chaotic, which are widely used to evaluate the effectiveness of nonlinear time-series models [1]. Sunspot has a very significant impact on the Earth's magnetic field, climate, geology, animals, and plants and also endangers the living environment of human beings. erefore, the observation and prediction of the sunspot numbers are of great significance. At present, with a large amount of long-term historical data, researchers have proposed a variety of prediction methods to forecast years, monthly mean, and peak of sunspots, such as chaos theory [2,3], various neural networks, and their optimization algorithms. Among them, Ding et al. [4] used BP neural network to predict monthly mean of smoothed sunspot area. Zhao et al. [5] used RBF neural network to predict the smoothed monthly mean sunspot numbers, but the prediction error was gradually enlarged with the prolongation of forecasting time. Artificial neural networks inherently slowed down the learning speed and have poor fault tolerance and easily fell into local minima. ese defects have not been well improved, thus the prediction accuracy cannot be further improved. e rise of intelligent optimization algorithms has greatly improved these shortcomings of artificial neural networks. Guan et al. [6] used quantum-behaved particle swarm optimization to optimize the weights and thresholds of BP neural networks to predict sunspot numbers. Li et al. [7,8] used genetic algorithms to optimize BP neural networks to predict shortterm traffic flow chaotic time series. e experimental results show that the optimized models are better.
Empirical mode decomposition (EMD), wavelet decomposition, and other methods provide an idea for processing nonlinear signals. In references [9,10], EMD has been applied to the decomposition prediction of nonlinear time series, but the EMD has the effect of mode mixing and point effects [11]. e introduction of ensemble empirical mode decomposition (EEMD) [12] and variational mode decomposition (VMD) [13] has improved the mode mixing to some extent. Especially VMD transforms signal decomposition into nonrecursive variational mode decomposition mode. Its number of components is also less than EMD and EEMD, and it shows better noise robustness.
In recent years, hybrid models of EMD, EEMD, and VMD combined with artificial neural networks have been widely used in the field of prediction [9,10,14,15]. e hybrid forecasting models formed by combining VMD with other algorithms have been successfully applied in many fields, such as financial time series [16], stock price evaluation [17], wind speed [18], and so on. Lahmiri [16] developed a new model of VMD and general regression neural network (GRNN) to analyze and predict economic and financial time series. In reference [17], the VMD and BP neural network were combined to predict the stock price well. In reference [18], the VMD was used to predict the wind speed sequence and then greatly reduced the data complexity and improved the prediction accuracy. erefore, this paper proposes a novel hybrid forecasting model based on VMD and firefly algorithm (FA) to optimize BP neural network (FA-BP) and applies the model to the forecast of sunspot numbers.

Variational Mode Decomposition.
VMD is an effective signal decomposition method. Its overall framework is a variational problem [13,19], which is different from the EMD of circular filtering. Its each mode is assumed to be a finite bandwidth with a different center frequency and the goal is to minimize the sum of the estimated bandwidths for each mode. e algorithm can be divided into the structure and solution of variational problem. e detailed description is as follows.

2.1.1.
e Structure of Variational Problem. Assume the original signal f is decomposed into K modal functions u k (t)(k � 1, 2, . . . , K), such that the variational problem can be described as a constrained variational problem. e problem is targeted at minimizing the sum of the estimated bandwidth for each mode, and the constraint is that the sum of modes is equal to the input signal f, as follows:

e Solution of Variational Problem
(a) e above constrained variational problem can be changed into a nonbinding variational problem by introducing a quadratic penalty factor C and Lagrange multipliers θ(t). Where C guarantees the reconstruction accuracy of the signal and θ(t) maintains the rigor of the constraint. e augmented Lagrange is denoted as (b) e alternate direction multiplier method (ADMM) is used to solve the saddle points of the above variational problem, and then the u n+1 k , ω n+1 k , and θ n+1 are updated alternately (n represents the number of iterations), which is given by (c) Given a discriminant accuracy e > 0, the convergence condition of the stop iteration is as follows: e specific process of the VMD algorithm is summarized as follows: (3) Judge whether or not u k meets the convergence condition (6), and repeat the above steps to update parameters until the convergence stop condition is satisfied. (4) e corresponding modal subsequences are obtained according to the given modal number.  [20,21]. e network prediction model mainly realizes the training process through the forward propagation of the input signal and the backpropagation of the error signal. It can parallelize large-scale data and has a certain degree of robustness and fault tolerance. e typical structure of three-layer BP neural network is shown in Figure 1.

FA-BP
In Figure 1, . , Y n is output value of BP neural network, and Y k is the expected output. e input vectors are propagated layer by layer from input layer, hidden layer, and output layer. e number of input layer nodes and output layer nodes are determined by the features and dimensions of training samples, respectively. e number of hidden layer nodes needs to be selected artificially according to the situation. Setting the transfer function between the layers as used sigmoid function, the output of the jth node of the hidden layer is where ω ij is the connection weight of the ith node of the input layer to the jth node of the hidden layer, a j is the threshold of the jth node of the hidden layer, and l is the number of hidden layer nodes. e hidden layer output H j is further passed back to get the kth node of the output layer as where ω jk is the connection weight of the jth node of the hidden layer to the kth node of the output layer, b k is the threshold of the kth node of the output layer, and m is the number of output layer nodes. According to the network prediction output Y k and the expected output Y h , calculate the network prediction error E: BP neural network corrects the weights and thresholds of all layers along the direction of error reduction and repeatedly modifies the weights and thresholds until the algorithm converges to obtain satisfactory error precision.

Firefly Algorithm.
Firefly algorithm (FA) was developed by Cambridge scholar Yang [22]. It is a swarm intelligence optimization algorithm based on the idealized behavior of the flashing characteristics of fireflies. Comparing with other optimization algorithms, the main advantages of firefly algorithm are simple structure, less adjustment parameters, and better spatial search capability [23]. erefore, the application field of the algorithm is also quite extensive such as function optimization, path planning, resource management, and so on [24][25][26][27].
In this algorithm, the brightness of a firefly and mutual attractiveness are the key factors that determine the firefly's movement.
e fitness value of the firefly's position determines its brightness. e brighter the brightness, the better its position. Attractiveness and brightness are related. e higher brightness of the fireflies has a stronger attraction, so other fireflies of lower brightness approach in this direction.
rough the interaction of brightness and attractiveness, fireflies will eventually gather in the highest brightness position to achieve the optimization of objective function. In order to facilitate the establishment of the algorithm model, the following definition is given: Definition 1. e relative brightness of a firefly: where I 0 is the original light intensity which is related to the objective function value, c is the light absorption coefficient, generally defined as a constant 1, and r ij is the Cartesian distance between firefly i and firefly j. e expression is as follows: Definition 2. e attractiveness of a firefly is as follows: where β 0 is the maximum attractiveness which represents the degree of attraction at the location of the maximum brightness. Equation (12) describes the property that the attractiveness of fluorescence emitted by fireflies to other individuals decreases with distance and the absorption coefficient of the propagation medium. β 0 can be taken as 1 for most application problems. e value of c has a great influence on the performance of the algorithm. eoretically, we can take c ∈ [0, ∞), but in practice we take usually c ∈ [0.1, 10].
Definition 3. Updated formula for the position of firefly i moved by firefly j is as follows: where t is the iteration number of the algorithm, x i , x j are the spatial positions of the firefly i and j, α is the step factor, and taking the constant at [0, 1], rand is the random factor that obeys the uniform distribution on [0, 1]. e process of firefly optimization algorithm is as follows: (1) Initialize the basic parameters of FA and set the number of firefly population m, maximum attractiveness β 0 , light absorption coefficient c, step size α, the maximum number of iterations max T, and search accuracy ε. (2) Randomly initialize the location of fireflies and calculate the target value of fireflies which is used as their maximum brightness I 0 . (3) Calculate the relative brightness I and attractiveness β of the fireflies in the group by (10) and (12) and Computational Intelligence and Neuroscience determine the direction of the firefly's movement according to the relative brightness. (4) According to Equation (13), to update the spatial location of fireflies, the firefly at the best position is randomly perturbed by α(rand − 0.5) to prevent premature convergence and fall into the local optimal solution. (5) According to the location of the updated firefly, recalculate the brightness of fireflies. (6) Turn on (7) when the search accuracy or the maximum number of searches is satisfied. Otherwise, the number of searches will be increased by 1 and turn (3) to the next search. (7) Output the global extreme point and the optimal individual value.

FA-BP Prediction Model.
In the BP neural network, the two kinds of training parameters of weight matrices ω ij , ω jk and thresholds a j , b k have significant influences on the prediction accuracy. e basic principle of optimizing FA-BP neural network is to use the better global search ability of FA to optimize the topology, connection weights, and thresholds of neural network. Considering the network parameters of BP neural network as firefly individuals, the fitness function of firefly individuals and the training error output of the corresponding network model form a linear relationship in the process of network model training: According to the above equation, it shows that the smaller the training error, the greater the individual fitness value. In the firefly algorithm, the greater the fitness value, the greater the brightness, indicating that firefly individuals perform better on the optimization problem. By searching for the best individual, which is the optimal network parameter, the firefly position is constantly updated, and then the weights and thresholds are substituted into the BP neural network to predict the sample. Specific steps are as follows: (1) Determine the neural network structure according to the input sample characteristics and output requirements.

Hybrid VMD-FA-BP Forecasting Model.
In this section, the proposed VMD-FA-BP model is established for the monthly mean of sunspots forecasting. e basic structure of the hybrid forecasting method is as follows: (

Data Decomposition Preprocessing.
e time series of monthly mean sunspots data comes from the official website of the solar influence data analysis center (SIDC). e experiment selects the observation data from January 1917 to December 2016 as the experimental sample, a total of 1200 data. e time-series diagram is shown in Figure 4. Taking into the randomness of the monthly mean of sunspots, direct prediction there will be a greater error. In order to improve the accuracy of prediction, the data complexity needs to be reduced. e original sequence is decomposed to generate multiple subsequences by VMD. Set the number of subsequences K before proceeding with VMD decomposition. However, we can see from the previous test that for the average monthly sequence of sunspots, the following subsequences tend to be similar when K > 6, so we choose K � 6 in this paper. e original data VMD decomposition is shown in Figure 5. e same data sample is also decomposed by EMD as shown in Figure 6.    Figure 3: e flowchart of VMD-FA-BP model.

Computational Intelligence and Neuroscience
As shown in Figures 5 and 6, the time series of sunspots are decomposed into 6 IMFs by VMD, and 6 IMFs and 1 residual component are obtained by EMD decomposition. From the decomposition results, it can be seen that the IMF3 component decomposed by the VMD method is most similar to the original signal waveform. e waveform distortion is small, and the trend of the high-frequency component variation is relatively stable. So, it is good for prediction. Although the low-frequency component fluctuates greatly, the prediction error is limited. Because the final VMD forecast result is an accumulation of each IMF's forecast result, the characteristics of the VMD decomposition result help to improve the prediction accuracy. However, the IMFs discomposed by EMD have different characteristics. End-point effects occur during the decomposition process. Large fluctuations occur at both ends of the IMF component and affect the entire component sequence continuously. e error would be more serious and the decomposition results would be seriously distorted, especially for low-frequency IMF components. Moreover, the drastic fluctuations of the IMF will lead to large errors in the final forecast based on the EMD. erefore, IMFs decomposed by VMD are more suitable for the establishment of hybrid forecasting model than IMFs decomposed by EMD [28][29][30].

Performance Standards of Prediction Accuracy.
is study uses the following three error indicators to verify the validity and practicability of the proposed prediction model: mean absolute error (MAE), root mean squared error (RMSE), and mean absolute percentage error (MAPE). e error of prediction value is quantified by using the performance indexes of MAE, RMSE, and MAPE. e smaller the value, the better the prediction accuracy. ese formulas are as follows: where x(t) is forecast data and x(t) is original data. In order to avoid the errors caused by the difference and randomness of each prediction result, this paper averages the results of 10 predictions, that is, the error value of the final prediction.  0  200  400  600  800  1000  1200  -10  0  10  IMF1   IMF2   0  200  400  600  800  1000  1200  -50  0  50   IMF3   0  200  400  600  800  1000  1200  -50  0  50   IMF4   0  200  400  600  800  1000   Computational Intelligence and Neuroscience used as the prediction samples.

Data Prediction Processing and Results
e BP neural network prediction model adopts the network structure of 5-8-1 and adopts the same structure for the FA-BP network prediction model, that is, the number of nodes in the input layer is 5, the number of nodes in the input layer is 8, and the number of nodes in the input layer is 1. ere are 5 × 8 + 8 × 1 � 48 weights, 8 + 1 � 9 thresholds. e individual code length of the firefly algorithm is 48 + 9 � 57. e parameters are set as follows: the number of fireflies m � 30, the maximum attraction β 0 � 1, the light absorption coefficient c � 1, the step factor α � 0.2, and the maximum number of iterations max T � 100. e above parameter settings and sample data are applied to all model tests in this article to ensure fair and efficient comparisons between different forecast models.
According to the above settings, each IMF after VMD decomposition is predicted and reconstructed to obtain the final prediction value. e prediction result is shown in Figure 7. It can be seen from the figure that the blue line represents the actual monthly mean of sunspots and the red line represents the predicted value of the hybrid forecasting model. It can be seen that the VMD-FA-BP model proposed in this paper is good for fitting the original data and can predict the monthly mean of sunspots very well. In order to facilitate comparison, BP neural network, FA-BP, EMD-BP, and VMD-BP prediction models are used to predict the same sunspot monthly mean time series, as shown in Figure 8. e partial enlargement is shown in Figure 9. While the forecasting performance evaluation results are displayed in Table 1 and Figure 10, respectively.
It can be seen from Table 1 and Figure 9 that the three error indicators of the VMD-FA-BP prediction model proposed in this paper are less than other prediction models, where MAE � 1.2208, RMSE � 1.7117, and MAPE � 0.0435, and its forecasting accuracy is improved by one order of magnitude compared with other models.
is obvious change is chiefly due to the following reasons. First, as a new adaptive multiresolution technique, the VMD is more robust to analyze noisy signals such as sunspot data than EMD and reduces the mode mixing existed in the EMD, making each IMF used for the following predictions better. Secondly, the processing time required by VMD, EMD, and BP neural network is significantly short. When processing monthly mean sunspot time series, the VMD and EMD computational time is 10s30 and 12s26, respectively. e BP neural network takes on average 0.95 s for learning and testing, compared with the traditional autoregressive moving average model in [14,28], which takes more time for modeling and forecasting (more than 5 min) and saves a lot of time and improves the prediction efficiency and accuracy. en, the firefly algorithm is a feasible and effective group intelligent optimization algorithm with high optimization precision and convergence speed. e optimal initial weights and thresholds of the BP neural network are obtained by using the firefly optimization algorithm, and the FA-BP model is applied to predict each IMF, and the obtained effect is better than a single BP model. Finally, the predicted values of each IMF are accumulated to obtain the final predicted values, thus the best prediction accuracy is shown in all the models of Table 1.
Hybrid model is more beneficial for modeling and improving prediction accuracy by using the advantages of each single model. It is proved that the variational mode decomposition combined with BP neural network prediction model using firefly algorithm in this paper can predict the variation trend of the monthly mean of sunspots time series well, which is a good prediction model.

Conclusions
In this paper, aiming at the problem of forecasting the monthly mean time series of sunspots, a hybrid forecasting model based on variational mode decomposition and firefly algorithm to optimize BP neural network is proposed. Firstly, a set of intrinsic mode functions (IMFs) are obtained Computational Intelligence and Neuroscience by VMD decomposition of the monthly mean of the sunspots. Secondly, a prediction model is established for each IMF to improve the prediction accuracy. In order to overcome the defects of the slow learning speed and easy to fall into local minimum in BP neural network, firefly algorithm is used to optimize the BP neural network. e optimal weights and thresholds are obtained. e training samples of the network are trained, and then the test samples are predicted. Finally, the predicted values of these components are summed to obtain the final prediction result, and the predicted value is compared with the real value to calculate the error. e experimental results show that the combination method proposed in this paper can be used to predict the monthly mean of sunspots and improve the prediction accuracy and reduce the error compared with other models. erefore, the model proposed in this paper has a good reference value for forecasting actual problems accurately such as short-time traffic flow and other time series.
Data Availability e smoothed monthly mean sunspot number time series in this paper comes from the data analysis center of the Sun Observatory of Belgium (http://sidc.oma.be/silso/datafiles). It records the sunspot data from 1749 to the present. In particular, the data analysis center significantly revised the number of sunspots on July 1, 2015. e data in this paper were derived from the revised dataset.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.    Computational Intelligence and Neuroscience