Wind power prediction method based on improved LSTM and FA-KELM

Abstract

The invention discloses a wind power prediction method based on improved LSTM and FA-KELM, and the method comprises the steps: obtaining the wind power data of a wind power plant, carrying out the normalization of the wind power data, and carrying out the decomposition of the wind power data through employing CEEMDAN; calculating the fuzzy entropy of the obtained subsequences, performing secondary decomposition on the subsequence with the maximum fuzzy entropy by using a PSO-VMD algorithm, and dividing other sequences into a high-entropy sequence and a low-entropy sequence according to the fuzzy entropy; the method comprises the following steps: improving an LSTM network, and replacing an input gate, a forgetting gate and an output gate in the LSTM network with a gate control structure with a cross coupling structure to obtain ConvsLSTM; mIC is adopted to screen out meteorological factors with large correlation with wind power in the meteorological factors of the wind power plant; and in combination with the obtained meteorological factor data, predicting the obtained high-entropy sequence by using ConvsLSTM, and predicting the low-entropy sequence by using FA-KELM. And superposing the predicted values of all the components, and carrying out reverse normalization to obtain a final ultra-short-term wind power prediction result. According to the method, the appropriate prediction model can be selected for the power subsequences with different complexities, and the overall wind power prediction precision is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of wind power prediction, and in particular to a wind power prediction method based on improved LSTM and FA-KELM. Background Art

[0002] As a leader in clean and renewable energy, wind power leads the transformation of energy structure with its pollution-free and renewable characteristics, but its intermittent and instability have also become new challenges for grid dispatching and stable operation. Accurately predicting wind power has become the key to improving the rationality, economy, safety and stability of grid planning. In the face of the nonlinear and non-stationary fluctuation characteristics of wind power, in-depth research on prediction technology is of great significance to promoting large-scale grid connection of wind power and ensuring the safety of power systems. In recent years, new prediction technologies such as artificial intelligence, machine learning, and hybrid models based on optimization algorithms have been widely used. These technologies can capture the short-term and long-term fluctuation trends of wind power more accurately through deep learning and analysis of a large amount of historical data, significantly improving the accuracy and reliability of wind power prediction.

[0003] Reference [1]: Sun Kun, Zhao Mengmeng, Shen Meina, et al. Random forest wind power prediction based on CEEMD and fuzzy entropy [J]. Smart Power, 2019, 47(10): 36-43. The power is decomposed into multiple subsequences using complete total empirical mode decomposition, and the fuzzy entropy value is used as the evaluation index for reorganization to reassemble subsequences with similar complexity into new sequences. Then, a random forest model is established for sequences with different fluctuation attributes and the model parameters are optimized for prediction. However, CEEMD is less effective in dealing with modal aliasing problems and cannot adaptively adjust the noise amplitude, resulting in low decomposition accuracy.

[0004] Reference [2]: Li Yu. Probabilistic prediction of wind power based on improved BNN-LSTM [J]. Microcomputer Applications, 2024, 40(03): 206-209. The TCNN network is used to process wind power historical data, extract the correlation features of time series data, and combine the mutual information entropy to find the factors with greater correlation with the meteorological data set of wind power, and use the Bayesian long short-term memory neural network (BNN-LSTM) for prediction. However, Bayesian inference introduces a complex calculation process, especially in large-scale data or deep networks, which significantly increases the training time, is difficult to train, and is prone to overfitting, limiting its effectiveness in practical applications.

[0005] Reference [3]: Ma Zhixia, Zhang Lin, Bayintana, et al. Ultra-short-term wind power forecasting based on adaptive quadratic decomposition and CNN-BiLSTM [J]. Acta Energiae Solaris Sinica, 2024, 45(06): 429-435. The wind power data is decomposed into multiple subsequences using the improved fully adaptive noise ensemble empirical mode decomposition (ICEEMDAN), and the high-frequency components are quadratically decomposed using the sparrow search algorithm optimized variational mode decomposition (SSA-VMD). The convolutional neural network (CNN) and the bidirectional long short-term memory network (BiLSTM) are used to predict the wind power in the ultra-short term. Although the single model prediction is applied to the decomposed subsequences, the modeling process is simplified, but the characteristics of each frequency band are not differentiated, and the complex dynamic relationship between different frequency bands is ignored, resulting in low prediction accuracy. Summary of the invention

[0006] The present invention provides a wind power prediction method and system based on improved LSTM and FA-KELM, which can more deeply represent the inherent change law of wind power sequence. The proposed combined prediction strategy can more effectively select appropriate prediction models for power subsequences of different complexities, thereby improving the overall wind power prediction accuracy.

[0007] The technical solution adopted by the present invention is:

[0008] The wind power prediction method based on improved LSTM and FA-KELM includes the following steps:

[0009] Step 1: Obtain wind power data of the wind farm, normalize all wind power data, and use fully adaptive noise ensemble empirical mode decomposition (CEEMDAN) to decompose the wind power data;

[0010] Step 2: Calculate the fuzzy entropy of the subsequences obtained in step 1, and use the PSO-VMD algorithm to perform secondary decomposition on the subsequence with the largest fuzzy entropy. The remaining sequences are divided into high entropy and low entropy sequences according to the fuzzy entropy.

[0011] Step 3: Improve the traditional LSTM network by introducing convolution operation in the input-state and state-state conversion process of the LSTM network, and replace the input gate, forget gate and output gate in the LSTM network with a new gating structure with a cross-coupling structure to obtain a convolutionally simplified long short-term memory network ConvsLSTM;

[0012] Step 4: Use the maximum correlation coefficient MIC to screen out the meteorological factors of the wind farm that have a greater correlation with wind power;

[0013] Step 5: Combined with the meteorological factor data obtained in step 4, the convolutional simplified long short-term memory network ConvsLSTM is used to predict the obtained high entropy sequence, and the firefly optimized kernel extreme learning machine FA-KELM is used to predict the low entropy sequence. Step 6: Superimpose the predicted values ​​of all components, denormalize, and obtain the final ultra-short-term wind power forecast result, using the mean absolute error MAE, mean square error MSE and determination coefficient R 2 These three evaluation indicators are used to compare and analyze the prediction results.

[0014] The step 1 comprises the following steps:

[0015] Step S1.1: First, normalize the wind power data to the interval [0,1]. The normalization formula is:

[0016]

[0017] Where, X is the original wind power data; X norm is the normalized data; X max and X min are the maximum and minimum values ​​of the original data respectively;

[0018] Step S1.2: Use CEEMDAN to decompose the normalized wind power data. The core idea is to add Gaussian white noise to the power sequence to be decomposed, and use empirical mode decomposition (EMD) to reduce modal aliasing. Through multiple iterations, the intrinsic mode function (IMF) and residual of the sequence are gradually decomposed. The specific steps include the following:

[0019] 1). Let the original power sequence be s(t), and add K times Gaussian white noise β with a mean of 0 k E k (w), where β k is the noise intensity, E k (w) is the k-th random noise, w is the noise sample;

[0020] 2). For each addition of noise, construct a new sequence s k (t) = s(t) + β k E k (w), where k = 1, 2, ..., K;

[0021] 3). For each new sequence s k (t) Perform EMD decomposition to obtain a series of IMFs and residuals. The nth subsequence obtained by the kth decomposition is called IMF. k,n (t);

[0022] 4). Calculate the mean of all subsequences of the first decomposition as the first IMF:

[0023]

[0024] In the formula, IMF k,1 (t) represents the first subsequence obtained after the k-th decomposition.

[0025] 5). Calculate the initial residual: r1(t)=s(t)-IMF1(t).

[0026] For n ≥ 2, the steps to calculate the nth IMF are as follows:

[0027] a). Construct a new sequence r n-1 (t)+β k E k (w), where r n-1 (t) is the residual after the n-1th iteration.

[0028] b). EMD decomposition: for the new sequence r n-1 (t)+β k E k (w) Perform EMD decomposition to obtain the first IMF k,1 (t);

[0029] c) Calculate IMF: By calculating all IMF k,1 (t), and subtract the impact of the previously calculated IMF (if any), to obtain the nth IMF:

[0030]

[0031] In the formula, IMF n (t) represents the nth subsequence obtained by the final decomposition, IMF k,1 (t) represents the first subsequence obtained by the k-th decomposition, IMF m (t) represents the previously calculated IMF, and m represents the sequence number of the mth subsequence calculated previously.

[0032] d). Update residual: r n (t) = r n-1 (t)-IMF n (t); r n (t) represents the residual of the nth subsequence, r n-1 (t) represents the residual of the n-1th subsequence, IMF n (t) represents the nth subsequence obtained by the final decomposition;

[0033] 6). Repeat step 5) until the stop condition is met.

[0034] The step 2 comprises the following steps:

[0035] Step 2.1: Calculating the fuzzy entropy of the wind power subsequence includes the following steps:

[0036] S2.1.1: Embed the wind power subsequence {u(i)} of length N into a space of dimension m, where u(i) represents the value of the subsequence at the i-th time point. Construct an embedding matrix, where each row X i = [u(i), u(i+1), …, u(i+m-1)] are all m-dimensional vectors, u(i+1) represents the value of the subsequence at the i+1th time point, u(i+m-1) represents the value of the subsequence at the i+m-1th time point i = 1, 2, …, N-m+1;

[0037] S2.1.2: Calculate any two rows in the embedding matrix X i and X j The distance d ij , which is usually defined as the absolute value of the maximum difference:

[0038]

[0039] In the formula, u(i+k) represents the value of the subsequence at the i+kth time point, u(j+k) represents the value of the subsequence at the j+kth time point, and k represents the number of time points k.

[0040] S2.1.3: Use a fuzzy membership function (such as a Gaussian function) to define the similarity μ(d ij ,r):

[0041]

[0042] In the formula, μ(d ij ,r) represents the similarity between the i-th sequence and the j-th sequence, and r is a predetermined threshold parameter used to control the decay rate of the similarity.

[0043] S2.1.4: Calculate the average similarity between each i and all other j:

[0044]

[0045] In the formula, φ m (r,i) is the average similarity between each i and all other j, where i and j represent the sequence numbers.

[0046] S2.1.5: Then calculate φ for all i m The average value of (r,i)Φ m (r):

[0047]

[0048] S2.1.6: Increase the embedding dimension to m+1 and repeat steps S2.1.1 to S2.1.5 to obtain Φ m+1 (r).

[0049] S2.1.7: Calculate the fuzzy entropy of the wind power subsequence:

[0050]

[0051] Where FE(m,r,N) represents the value of fuzzy entropy.

[0052] Step 2.2: Use VMD to decompose the wind power subsequence with the highest fuzzy entropy:

[0053] First, construct a variational problem, constraining all subsequences u k The sum of (t) is equal to the original sequence, which can be expressed as:

[0054]

[0055] Where f(t) represents the input wind power sequence, K represents the number of decomposed sequences, and u k (t) represents the kth sequence function, ω k represents the center frequency of the kth sequence function, t represents time, and δ(t) is the pulse function.

[0056] In order to solve the above variational problem, the Lagrangian multiplier λ(t) is introduced to obtain the augmented Lagrangian function L:

[0057]

[0058] In the formula, α represents the quadratic penalty parameter, which is used to balance the fidelity of the sequence and the compactness of the decomposition. k ,ω k ,λ) represents the augmented Lagrangian function; represents the bandwidth limitation term in VMD; f(t) represents the input power sequence; represents the regularization term.

[0059] The saddle point of the above augmented Lagrangian function is solved by the alternating direction multiplier method, that is, each power subsequence u is found k and the center frequency ω k For each iteration step, update u k and ω k The formula is as follows:

[0060]

[0061] In the formula, n represents the number of iterations; and Represent f(t), u i (t), λ(t), and Fourier transform of ; ω represents frequency, Represents the center frequency obtained by the n+1th iteration.

[0062] Through iterative updates, until the stopping condition is met, K subsequences u are obtained k and the corresponding center frequency ω k .

[0063] Step 2.3: Use the PSO algorithm to optimize the parameters of the VMD decomposition power subsequence, including the following steps:

[0064] S2.3.1: Randomly generate particle positions and velocities. The position X and velocity V of each particle represent K,α.

[0065] S2.3.2: For each particle, perform VMD decomposition using its position parameters and calculate the reconstruction error as the fitness value. The reconstruction error is defined as the difference between IMF1x(t) after the first decomposition and all the u k The difference between the sum of (t):

[0066]

[0067] Where Reconstruction Error represents the reconstruction error.

[0068] S2.3.3: Update the position and velocity of the particle according to the velocity update formula and position update formula.

[0069] Position update formula:

[0070]

[0071] Where, X i (t+1) represents the position information of the i-th particle in the t+1 generation, X i (t) represents the position information of the i-th particle in the t-th generation, V i (t+1) Represents the velocity information of the i-th particle in the t+1th generation.

[0072] Speed ​​update formula:

[0073]

[0074] Where σ is the inertia weight, c1 and c2 are acceleration factors, r1 and r2 are random numbers ranging from 0 to 1, and P i best is the optimal historical position of particle i itself, Gbest is the best historical position of all particles, V i (t+1) represents the velocity information of the i-th particle in the t+1 generation, V i (t) represents the velocity information of the i-th particle in the t-th generation, X i (t) Represents the position information of the i-th particle in the t-th generation.

[0075] S2.3.4: Update the historical optimal position P of each particle according to the fitness value i best and the global optimal position G among all particles best .

[0076] S2.3.5: Repeat the above steps until the stopping condition is met. The stopping condition is that the maximum number of iterations is reached or the error converges to a certain threshold.

[0077] In the step 3,

[0078] The LSTM network is a neural network that efficiently processes time series through memory units. Each unit contains an input gate, a forget gate, and an output gate. Its main function is to selectively forget information that does not need to be retained for a long time by tracking information from input to output and the information transmitted to the next unit, while updating valuable information and retaining it in the memory unit for a long time. In view of the problem that the LSTM network lacks the ability to process spatiotemporal data and has a long operation time, the present invention proposes a long short-term memory network ConvsLSTM combined with convolution simplification.

[0079] First, the convolution operation is introduced in the input-state and state-state conversion process of the LSTM network to enhance the network feature extraction capability. The forward propagation process is as follows:

[0080] i t =σ(W xi *x t +W hi *h t-1 +b i )

[0081] f t =σ(W xf *x t +W hf *h t-1 +b f )

[0082]

[0083] o t =σ(W xo *x t+W ho *h t-1 +b o )

[0084]

[0085] In the formula, * is the convolution operation, is the Hadamard product operation, W xi , W xf , W xc and W xo represents the input weight, W hf , W hi , W hc and W ho is the cycle weight, b i 、b f 、b c and b o is regarded as the input bias, σ(·) is the activation function, tanh(·) is the hyperbolic activation function, and f t 、i t 、c t , o t and h t represents the forget gate, input gate, memory cell state, output gate and network output of LSTM, h t-1 represents the network output at the previous time, c t-1 Indicates the state of the memory unit at the previous time, x t Represents the input of the power subsequence.

[0086] The built-in gating structure of the LSTM network involves multiple sets of parameters, which need to be iteratively updated during network training. In order to further reduce the number of LSTM parameters to be optimized and its training time, a new gating structure with a cross-coupling structure is used to replace the input gate, forget gate, and output gate in the LSTM network. Compared with the original LSTM network, the introduction of convolution operations improves the ability to process spatiotemporal data and reduces the model training time by sharing the parameters of each gating structure. At the same time, a peephole structure is introduced to consider the impact of storage units on network output. The structure of the proposed ConvsLSTM is shown in Figure 1. Figure 1 As shown, the forward propagation process is shown as follows:

[0087] i t =σ(W x *x t +W h *h t-1 +W c *c t-1 +b i )

[0088]

[0089] o t =σ(W x *x t +W h *h t-1 +W c *c t +b o )

[0090]

[0091] In the formula, * is the convolution operation, is the Hadamard product operation, W x is the input weight, W h is the cycle weight, W c is the peephole weight, b i , b c and b o is the input bias, σ(·) is the activation function, tanh(·) is the hyperbolic activation function, i t 、c t , o t and h t They represent the input gate, memory cell state, output gate and network output of LSTM respectively. t-1 represents the network output at the previous time, c t-1 Indicates the state of the memory unit at the previous time, x t Represents the input of the network.

[0092] In step 4,

[0093] MIC (maximum information coefficient) can capture the potentially complex relationship between wind farm meteorological factors and wind power by measuring the nonlinear correlation between two variables. In the analysis of wind farm meteorological factors, MIC can be used to screen meteorological factors that have a strong correlation with wind power, such as wind speed, wind direction, temperature, etc. The basic idea of ​​MIC is to grid the data, calculate the mutual information under each grid division, and then find the maximum value of the mutual information divided by the complexity of the grid division. The complexity of the grid division is usually measured by the number of grids or the degree of refinement of the grid division. The calculation formula of the maximum correlation coefficient MIC is expressed as:

[0094]

[0095] Where X and Y are two variables, a×b represents the grid division method, that is, X is divided into a intervals, Y is divided into b intervals, and the value of a×b is less than a given upper limit B to avoid overfitting, I(X; Y| G) represents the mutual information between X and Y under a given grid G ​​partitioning, log2min(a,b) is the complexity measure of the grid partitioning, which is used to normalize the mutual information.

[0096] In the present invention, by calculating the MIC value between different meteorological factors (such as wind speed, wind direction, temperature, etc.) and wind power, the impact of these factors on wind power can be quantified. Meteorological factors with higher MIC values ​​indicate that they have a stronger correlation with wind power, and therefore can be given priority for use in the wind power prediction model. Conversely, factors with lower MIC values ​​may have less impact on power, and perhaps their weight can be reduced or excluded in the model. In this way, MIC can screen out the most influential meteorological factors, thereby optimizing the prediction model and improving prediction accuracy and robustness. This process ensures that the wind power prediction model can make full use of the information of key meteorological factors and improve the operating efficiency and prediction ability of the wind farm.

[0097] The step 5 comprises the following steps:

[0098] Step S5.1: Extreme Learning Machine (ELM) is an efficient single hidden layer feedforward neural network that shortens the training time by randomly generating the weights and biases of the hidden layer. However, ELM performs poorly in dealing with nonlinear wind power prediction problems. To solve this problem, KELM introduces a kernel function based on ELM, replacing the random mapping in ELM with a stable kernel mapping, so that the model can more effectively process nonlinear data and improve the prediction accuracy and generalization ability. The algorithm steps are as follows:

[0099] The output function of the ELM network with L hidden neurons can be expressed as:

[0100]

[0101] In the formula, f elm (x) represents the output function, x is the input vector, x j represents the specific input value of each input vector, β i is the weight connecting the hidden neuron and the output neuron, g(i) is the activation function, ω i is the randomly assigned weight vector connecting the input nodes and hidden neurons, b i represents the bias of the random assignment of hidden layer neurons, j represents the ordinal number of the value in the input vector, and N represents the number of values ​​in the input vector.

[0102] The above formula can be redefined as:

[0103] Hβ=T

[0104]

[0105] In the formula, H represents the output matrix of the hidden layer, T represents the predicted expected output matrix, β is the weight connecting the hidden neurons and the output neurons; x1 represents the first number of the input vector, ω1 is the first weight vector, b1 represents the deviation of the random assignment of the first hidden layer neuron, ω L represents the Lth weight vector, b L represents the bias of the random assignment of the Lth hidden layer neuron, x N Represents the Nth number of the input vector.

[0106] In ELM, the above formula can be regarded as a linear function, and its least squares solution is:

[0107] β=H + T

[0108] In the formula, H + is the generalized inverse matrix of H.

[0109] H + =H T (HH T ) -1

[0110] In the formula, H T represents the transposed matrix of H, HH T represents a symmetric matrix.

[0111] The regularization coefficient and the identity matrix are introduced to solve the least squares solution of the hidden layer node output weight β:

[0112]

[0113] Where C is the regularization coefficient and I represents the N-order unit matrix.

[0114] KELM uses the kernel function matrix Ω ELM The random matrix HH that replaces ELM T , the present invention adopts RBF function as kernel function:

[0115]

[0116] Where σ represents the kernel parameter; K(x i ,x j ) represents the kernel function, x i and x j Represents the i-th and j-th data points in the dataset.

[0117] Combination of the above and Hβ=T, we can get the output expression of KELM:

[0118]

[0119] Where f(x) represents the output of KELM; K(x,x1) represents the kernel function used to measure the similarity between the new sample x and the first sample x1 in the training set; K(x,x N ) indicates that the kernel function is used to measure the new sample x and the Nth sample x in the training set. N The similarity between .

[0120] Step S5.2: In order to improve the accuracy of KELM in predicting wind power, the firefly algorithm (FA) is used to optimize the penalty coefficient C and kernel parameter σ of KELM;

[0121] First, initialize N fireflies, and the position of each firefly is x i =(σ i ,C i ), i=1,2,…,N,σ i represents the i-th kernel parameter, C i represents the i-th penalty coefficient, N represents the number of fireflies, x i Represents a solution of KELM, and defines the mean square error (MSE) of the KELM model on the validation set as the objective function:

[0122]

[0123] In the formula, f(x i ) represents the mean square error, y j is the actual value, is the predicted value of KELM, and n represents the length of the input sequence.

[0124] The brightness of each firefly I i =-f(x i ), the brighter fireflies attract the brighter fireflies, and the position of the brighter fireflies is x j , the position of the firefly with low brightness is x i , the formula for updating the position is:

[0125]

[0126] In the formula, r ij =x i -x j is the distance between fireflies i and j, β0 is the maximum attraction coefficient, γ is the light intensity absorption coefficient, rand is a random number between [0,1], and α is the step size coefficient.

[0127] After several iterations, find the firefly x with the highest brightness best =(σ best ,C best ), σbest represents the optimal kernel parameter, C best Represents the optimal penalty coefficient. When the maximum number of iterations is reached or the brightness change is less than the preset threshold, the optimization algorithm ends and the optimal parameters of KELM are obtained.

[0128] In step 6, the formula for denormalizing the prediction results is:

[0129] X=(X norm *(X max -X min ))+X min

[0130] Where X is the original wind farm power data, X norm is the normalized data; X max and X min are the maximum and minimum values ​​of the original data respectively.

[0131] The three evaluation indicators include:

[0132]

[0133]

[0134] Where n represents the length of the prediction sequence, y and Represent the actual value and predicted value of power respectively, is the average value of actual power.

[0135] On the other hand, the wind power prediction system based on improved LSTM and FA-KELM includes the following modules:

[0136] The data collection module collects the original power data set of the wind farm, wherein the original data set includes the wind speed, weather and power data of the wind farm with a sampling period of 15 minutes.

[0137] In the data preprocessing module, the original power data set is decomposed into multiple subsequences based on the combination of CEEMDAN and VMD, and the fuzzy entropy algorithm is used to perform fuzzy entropy classification on each subsequence.

[0138] In the model training module, the subsequences after fuzzy entropy classification are used as the input of ConvsLSTM and FA-KELM models for training to obtain the wind power prediction model.

[0139] The prediction module inputs the subsequence of the current wind farm to be predicted that has undergone data preprocessing into the trained wind power prediction model for prediction.

[0140] On the other hand, the present invention provides at least one memory and a processor; the memory is used to store computer programs and instructions, and the processor is used to execute the programs and instructions stored in the memory, so that the method described in the present invention can be implemented.

[0141] The present invention proposes a wind power prediction method and system based on improved LSTM and FA-KELM, and the technical effects are as follows:

[0142] 1) This paper verifies the effectiveness of the "decomposition-classification" strategy in wind power forecasting. Through the dual decomposition mechanism of CEEMDAN and VMD, the complex and volatile wind power series is decomposed into a series of relatively stable subsequences. This process significantly reduces the nonlinearity and non-stationarity of the original data, thereby reducing the burden of a single prediction model and achieving a significant improvement in prediction accuracy.

[0143] 2) Compared with the traditional prediction method based only on the single variable of wind power, the present invention innovatively introduces multi-source information such as weather conditions and wind speed changes as input features, further improving the accuracy of the prediction results.

[0144] 3) The subsequences obtained after decomposition were processed in detail according to their complexity. Specifically, subsequences of different complexities were respectively handed over to the ConvsLSTM and FA-KELM models for prediction. This fuzzy entropy-guided model customization strategy effectively utilizes the advantages of each algorithm within the specific complexity range of the sequence and further optimizes the prediction accuracy. 4) Experimental results show that the secondary decomposition strategy can further explore the inherent variation patterns of the original power data and reveal more potential prediction information, thereby achieving a significant improvement in prediction accuracy compared to a single decomposition. This discovery provides new ideas and methodological support for the prediction of complex time series. BRIEF DESCRIPTION OF THE DRAWINGS

[0145] Figure 1 It is a structural diagram of the ConvsLSTM unit of the present invention.

[0146] Figure 2 It is a prediction flow chart of the present invention.

[0147] Figure 3 It is the CEEMDAN decomposition result diagram of the present invention.

[0148] Figure 4 This is the result diagram of the secondary decomposition of IMF1 by PSO-VMD of the present invention.

[0149] Figure 5 It is the fuzzy entropy graph of the CEEMDAN subsequence of the present invention.

[0150] Figure 6It is the fuzzy entropy graph of the PSO-VMD secondary decomposition subsequence of the present invention.

[0151] Figure 7 It is a multi-factor correlation analysis diagram in the present invention.

[0152] Figure 8 is a comparison chart of the prediction results of the present invention and different models;

[0153] Fig. 9 for Figure 8 A partial enlarged view of point A. DETAILED DESCRIPTION

[0154] The present invention proposes an ultra-short-term wind power prediction model combining a convolutional simplified long short-term memory network (ConvsLSTM) and a firefly optimized kernel extreme learning machine (FA-KELM). First, the wind power sequence is secondary decomposed using fully adaptive noise ensemble empirical mode decomposition (CEEMDAN) and variational mode decomposition (VMD) optimized by a particle swarm algorithm; secondly, all subsequences are analyzed for complexity using fuzzy entropy to divide them into high entropy and low entropy sequences; then the maximum correlation coefficient (MIC) is used to screen out factors with a high correlation with wind power; further, a new gating structure of convolution operation and cross-coupling is introduced into the long short-term memory network (LSTM) to enhance the network feature extraction capability and improve the training efficiency, and the obtained ConvsLSTM model is used to predict the high entropy sequence, while the FA-KELM model is used to predict the low entropy sequence; finally, the predicted values ​​of each sequence are superimposed to obtain the final prediction result.

[0155] The wind power prediction method based on improved LSTM and FA-KELM first normalizes the historical data of the example of the present invention, then uses CEEMDAN to decompose the power sequence, calculates the fuzzy entropy of the subsequence, and then uses VMD to decompose the subsequence with the largest fuzzy entropy for a second time, and uses the PSO optimization algorithm to optimize the parameters. The subsequences decomposed twice use fuzzy entropy as an indicator of complexity and are divided into high entropy sequences and low entropy sequences. The high entropy sequence is predicted by using the LSTM improved by convolution, and the low entropy sequence is predicted by using the KELM optimized by FA. The predicted values ​​of different subsequences are superimposed and denormalized to obtain the final predicted value. The prediction process is as follows: Figure 2 shown.

[0156] Step a: Normalize the historical data of the wind farm and then use CEEMDAN to decompose the power series.

[0157] Step a1: First, normalize the historical data to the interval [0,1]. The normalization formula is:

[0158]

[0159] In the formula, X is the original wind power, X norm is the normalized data, X max and X min are the maximum and minimum values ​​of the original data respectively.

[0160] Step a2: Perform CEEMDAN decomposition on the normalized power data. The process is to add Gaussian white noise to the sequence to be decomposed and use EMD to reduce modal aliasing. It gradually decomposes the IMF and residual of the sequence through multiple iterations. The specific steps of the algorithm are as follows:

[0161] 1) Let the original power sequence be s(t), and add K times Gaussian white noise β with a mean of 0 k E k (w), where β k is the noise intensity, E k (w) is the k-th random noise, and w is the noise sample.

[0162] 2) For each addition of noise, construct a new sequence:

[0163] s k (t) = s(t) + β k E k (w)

[0164] In the formula, k = 1, 2,…, K.

[0165] 3) For each s k (t) Perform EMD decomposition to obtain a series of IMFs and residuals. The nth IMF obtained by the kth decomposition is called IMF k,n (t).

[0166] 4) Calculate the mean of all power subsequences of the first decomposition as the first IMF:

[0167]

[0168] 5) Calculate the initial residual r1(t)=s(t)-IMF1(t).

[0169] For n ≥ 2, the steps to calculate the nth IMF are as follows:

[0170] a) Construct a new sequence:

[0171] r n-1 (t)+β k E k (w)

[0172] In the formula, r n-1 (t) is the residual after the n-1th iteration.

[0173] b) EMD decomposition: Perform EMD decomposition on the new sequence to obtain the first IMF k,1 (t).

[0174] c) Calculate IMF: By calculating all IMF k,1 (t), and subtract the impact of the previously calculated IMF (if any), to obtain the nth IMF:

[0175]

[0176] d) Update the residual:

[0177] r n (t) = r n-1 (t)-IMF n (t)

[0178] 6) Repeat step 5 until the stop condition is met. The decomposition result is as follows: Figure 3 shown.

[0179] Step b: Calculate the fuzzy entropy of the subsequence obtained in step a. Figure 3 It can be seen that the IMF1 obtained after the first decomposition still has strong non-stationarity, and there will be large errors when directly predicted. PSO-VMD is used to perform a second sequence decomposition on the IMF1 with the largest fuzzy entropy:

[0180] Step b1: The steps of calculating the fuzzy entropy of the obtained power subsequence are as follows:

[0181] 1) Embed the power sequence {u(i)} of length N into a space of dimension m and construct an embedding matrix. Each row X of the matrix i =[u(i),u(i+1),…,u(i+m-1)] are all m-dimensional vectors, i=1,2,…,N-m+1.

[0182] 2) Calculate any two rows in the embedding matrix X i and X j The distance d ij :

[0183]

[0184] 3) Use fuzzy membership function (such as Gaussian function) to define similarity μ(d ij ,r):

[0185]

[0186] Where r is a predetermined threshold parameter that controls the decay rate of similarity.

[0187] 4) Calculate the average similarity between each i and all other j:

[0188]

[0189] 5) Then calculate φ for all i m The average value of (r,i)Φ m (r):

[0190]

[0191] 6) Increase the embedding dimension to m+1 and repeat steps 1) to 5) to get Φ m+1 (r).

[0192] 7) Calculate the fuzzy entropy of the power subsequence:

[0193]

[0194] Step b2: Use VMD to decompose the power subsequence with the highest fuzzy entropy. First, construct a variational problem to constrain all wind power subsequences u k The sum of (t) is equal to the original power sequence, which can be expressed as:

[0195]

[0196] In the formula, f(t) represents the input power sequence, K represents the number of decomposition sequences, and u k (t) represents the kth sequence, ω k represents the center frequency of the kth sequence, t represents time, and δ(t) is the pulse function.

[0197] In order to solve the above variational problem, the Lagrangian multiplier λ(t) is introduced to obtain the augmented Lagrangian function L:

[0198]

[0199] Where α represents the quadratic penalty parameter, which is used to balance the fidelity of the sequence and the compactness of the decomposition.

[0200] The saddle point of the above augmented Lagrangian function is solved by the alternating direction multiplier method, that is, each subsequence u is found k and the center frequency ω k For each iteration step, update u k and ω k The formula is as follows:

[0201]

[0202] In the formula, n represents the number of iterations, and Represent f(t), u i (t), λ(t), and The Fourier transform of .

[0203] Through iterative updates, until the stopping condition is met, K subsequences u are obtained k and the corresponding center frequency ω k .

[0204] Step b3: Use the PSO algorithm to optimize the sequence number K and penalty factor α in VMD. The process is as follows:

[0205] a) Initialization: Randomly generate the position and velocity of particles. The position X and velocity V of each particle represent K,α.

[0206] b) Calculate fitness: For each particle, perform VMD decomposition using its position and velocity parameters, and calculate the reconstruction error as the fitness value. The reconstruction error is defined as the difference between the original sequence x(t) and all subsequences u k The difference between the sum of (t):

[0207]

[0208] c) Update position and velocity: Update the position and velocity of the particle according to the velocity update formula and position update formula.

[0209] Position update formula:

[0210]

[0211] Where, X i (t) represents the position information of the i-th particle in the t-th generation, V i (t+1) Represents the velocity information of the i-th particle in the t+1th generation.

[0212] Speed ​​update formula:

[0213]

[0214] Where ω is the inertia weight, c1 and c2 are acceleration factors, r1 and r2 are random numbers (ranging between 0 and 1), and P i best is the optimal historical position of particle i itself, G best is the best historical position of all particles.

[0215] d) Update individual and global optimal solutions: Update the historical optimal position P of each particle according to the fitness value i best and the global optimal position G among all particles best .

[0216] e) Iteration: Repeat the above steps until the stopping condition is met (the maximum number of iterations is reached or the error converges to a certain threshold).

[0217] Other parameter settings of VMD: TAU is 1, TOL is 1e-7, DC is 0, and INIT is set to 1. The present invention uses the PSO particle swarm optimization algorithm to optimize the modal number K and penalty factor α of VMD, and obtains the optimal parameters K=5, α=1865. The decomposition results are as follows: Figure 4 As shown. Figure 4 It can be seen that the secondary decomposition can obtain more detailed frequency and time characteristics, reveal the complex structure that the first decomposition failed to completely separate, and more clearly identify the noise or fast fluctuation patterns.

[0218] Step c: Use fuzzy entropy to analyze the subsequences after CEEMDAN and VMD secondary decomposition. A high fuzzy entropy value indicates a high sequence complexity, and a low fuzzy entropy value indicates a low sequence complexity. Set the critical value λ0 = 0.25 to classify each sequence into a high entropy sequence or a low entropy sequence. The algorithm parameters are set to embedding dimension m = 2, and the threshold is 1.4 times the standard deviation. The fuzzy entropy values ​​and division results of all subsequences are shown in Figure 2. Figure 5 and Figure 6 shown.

[0219] Figure 5 and Figure 6 The classification results of each subsequence based on the fuzzy entropy value critical value of 0.25 are shown. Since high entropy subsequences often contain more irregularities and random fluctuations, more noise reduction or smoothing is required; while low entropy sequences are mostly trend parts in power, which are suitable for pattern recognition or feature extraction. Therefore, by classifying the entropy values ​​of subsequences, sequences with high complexity can be further preprocessed to ensure that the accuracy and stability of the model are improved in the subsequent analysis or prediction process.

[0220] Step d: Considering that the historical information of the wind farm contains meteorological information such as wind speed, temperature, air pressure, and historical power data, in order to avoid input redundancy, MIC is used to filter out the relevant meteorological data. The MIC values ​​of each factor and power are as follows: Figure 7 As shown by Figure 7 It can be seen that among the main factors affecting wind power, wind speed has the greatest correlation, so the wind speed at different heights of the wind tower and the wind speed of the generator hub are selected as multi-factor variables.

[0221] The calculation formula of MIC of wind farm meteorological information and wind power data is:

[0222]

[0223] Where X and Y are two variables, a×b represents the grid division method, that is, X is divided into a intervals, Y is divided into b intervals, and the value of a×b is less than a given upper limit B (to avoid overfitting), I(X; Y| G ) represents the mutual information between X and Y under a given grid G ​​partitioning, log2min(a,b) is the complexity measure of the grid partitioning, which is used to normalize the mutual information.

[0224] Step e: Combined with the meteorological factor data obtained in step d, the high entropy subsequence is predicted using ConvsLSTM, and the low entropy subsequence is predicted using FA-KELM. The optimal parameters of the FA-KELM prediction model for the low entropy subsequence are shown in Table 1.

[0225] Table 1 Optimal parameters of FA-KELM

[0226]

[0227] Step f: Denormalize the prediction results. The formula is:

[0228] X=(X norm *(X max -X min ))+X min

[0229] Where X is the power data after denormalization, X norm is the normalized data, X max and X min are the maximum and minimum values ​​of the original power data respectively.

[0230] The predicted values ​​of all subsequences are superimposed to obtain the final ultra-short-term wind power forecast result.

[0231] Step g: To verify the superiority of the model proposed in the present invention, multiple models such as CEEMDAN-(ConvsLSTM, FA-KELM), PSO-VMD-(ConvsLSTM, FA-KELM), CEEMDAN-ConvsLSTM, CEEMDAN-FA-KELM, PSO-VMD-ConvsLSTM, and PSO-VMD-FA-KELM are used for comparison. The prediction results of each model are shown in the figure. Figure 8 , Fig. 9 As shown in Table 2, it can be seen from the figure that the proposed model has a better effect in wind power prediction.

[0232] As shown in Table 2, compared with the single prediction model after one decomposition, the MAE of the model proposed in the present invention is reduced by 21.81%, 49.84%, and 25.41% respectively; the MSE is reduced by 31.71%, 71.4%, and 39.13% respectively; R 2 Compared with the model of combined prediction after one decomposition, MAE decreased by 7.49% and 61.55%; MSE decreased by 12.51% and 78.94%; R 2 The results show that the secondary decomposition can reveal the inherent variation law of wind power sequence more deeply, and the proposed combined forecasting strategy can more effectively select appropriate forecasting models for subsequences of different complexity, thereby improving the overall forecasting accuracy.

[0233] Table 2 Evaluation indicators of different models

[0234]

Claims

1. The wind power prediction method based on improved LSTM and FA-KELM is characterized by The following steps are involved: Step 1: Obtain wind power data of the wind farm, normalize the wind power data, and decompose the wind power data using fully adaptive noise ensemble empirical mode decomposition (CEEMDAN); Step 2: Calculate the fuzzy entropy of the subsequences obtained in step 1, and use the PSO-VMD algorithm to perform secondary decomposition on the subsequence with the largest fuzzy entropy. The remaining sequences are divided into high entropy and low entropy sequences according to the fuzzy entropy. Step 3: Improve the LSTM network by introducing convolution operations in the input-state and state-state conversion processes of the LSTM network, and replace the input gate, forget gate, and output gate in the LSTM network with a gate structure with a cross-coupling structure to obtain a convolutionally simplified long short-term memory network ConvsLSTM; Step 4: Use the maximum correlation coefficient MIC to screen out the meteorological factors of the wind farm that have a greater correlation with wind power; Step 5: Combined with the meteorological factor data obtained in step 4, the convolutional simplified long short-term memory network ConvsLSTM is used to predict the obtained high entropy sequence, and the firefly optimized kernel extreme learning machine FA-KELM is used to predict the low entropy sequence; Step 6: Superimpose the predicted values ​​of all components and denormalize them to obtain the final ultra-short-term wind power forecast result.

2. The wind power prediction method based on improved LSTM and FA-KELM according to claim 1 is characterized in that: The step 1 comprises the following steps: Step S1.1: First, normalize the wind power data to the interval [0,1]. The normalization formula is: Where, X is the original wind power data; X norm is the normalized data; X max and X min are the maximum and minimum values ​​of the original data respectively; Step S1.2: Use CEEMDAN to decompose the normalized wind power data, add Gaussian white noise to the power sequence to be decomposed, and use empirical mode decomposition (EMD) to reduce modal aliasing; through multiple iterations, gradually decompose the intrinsic mode function (IMF) and residual of the sequence.

3. The wind power prediction method based on improved LSTM and FA-KELM according to claim 2 is characterized in that: Step S1.2 specifically includes the following steps: 1). Let the original power sequence be s(t), and add K times Gaussian white noise β with a mean of 0 k E k (w), where β k is the noise intensity, E k (w) is the k-th random noise, w is the noise sample; 2). For each addition of noise, construct a new sequence s k (t) = s(t) + β k E k (w), where k = 1, 2, ..., K; 3). For each new sequence s k (t) Perform EMD decomposition to obtain a series of IMFs and residuals. The nth subsequence obtained by the kth decomposition is called IMF. k,n (t); 4). Calculate the mean of all subsequences of the first decomposition as the first IMF: In the formula, IMF k,1 (t) represents the first subsequence obtained by the k-th decomposition; 5). Calculate the initial residual: r1(t)=s(t)-IMF1(t); For n ≥ 2, the steps to calculate the nth IMF are as follows: a). Construct a new sequence r n-1 (t)+β k E k (w), where r n-1 (t) is the residual after the n-1th iteration; b). EMD decomposition: for the new sequence r n-1 (t)+β k E k (w) Perform EMD decomposition to obtain the first IMF k,1 (t); c) Calculate IMF: By calculating all IMF k,1 (t), and subtract the impact of the previously calculated IMF (if any), to obtain the nth IMF: In the formula, IMF n (t) represents the nth subsequence obtained by the final decomposition, IMF k,1 (t) represents the first subsequence obtained by the k-th decomposition, IMF m (t) represents the previously calculated IMF, and m represents the sequence number of the m-th subsequence that has been calculated previously; d). Update residual: r n (t) = r n-1 (t)-IMF n (t); r n (t) represents the residual of the nth subsequence, r n-1 (t) represents the residual of the n-1th subsequence, IMF n (t) represents the nth subsequence obtained by the final decomposition; 6). Repeat step 5) until the stop condition is met.

4. The wind power prediction method based on improved LSTM and FA-KELM according to claim 1 is characterized in that: The step 2 comprises the following steps: Step 2.1: Calculating the fuzzy entropy of the wind power subsequence includes the following steps: Step 2.2: Use VMD to decompose the subsequence with the highest fuzzy entropy; Step 2.3: Use PSO algorithm to optimize the parameters of VMD.

5. The wind power prediction method based on improved LSTM and FA-KELM according to claim 4 is characterized in that: Step 2.1 includes the following steps: S2.1.1: Embed the wind power subsequence {u(i)} of length N into a space of dimension m, where u(i) represents the value of the subsequence at the i-th time point; construct an embedding matrix, where each row X of the embedding matrix i = [u(i), u(i+1), …, u(i+m-1)] are all m-dimensional vectors, u(i+1) represents the value of the subsequence at the i+1th time point, u(i+m-1) represents the value of the subsequence at the i+m-1th time point i = 1, 2, …, N-m+1; S2.1.2: Calculate any two rows in the embedding matrix X i and X j The distance d ij , which is usually defined as the absolute value of the maximum difference: In the formula, u(i+k) represents the value of the subsequence at the i+kth time point, u(j+k) represents the value of the subsequence at the j+kth time point, and k represents the number of time points k; S2.1.3: Use fuzzy membership function to define similarity μ(d ij ,r): In the formula, μ(d ij ,r) represents the similarity between the i-th sequence and the j-th sequence, r is a predetermined threshold parameter used to control the decay speed of the similarity; S2.1.4: Calculate the average similarity between each i and all other j: In the formula, φ m (r,i) is the average similarity between each i and all other j, where i and j represent the sequence numbers; S2.1.5: Then calculate φ for all i m The average value of (r,i)Φ m (r): S2.1.6: Increase the embedding dimension to m+1 and repeat steps S2.1.1 to S2.1.5 to obtain Φ m+1 (r); S2.1.7: Calculate the fuzzy entropy of the wind power subsequence: Where FE(m,r,N) represents the value of fuzzy entropy.

6. The wind power prediction method based on improved LSTM and FA-KELM according to claim 4 is characterized by: In step 2.2, we first construct a variational problem, constraining all subsequences u k The sum of (t) is equal to the original sequence, expressed as: Where f(t) represents the input wind power sequence, K represents the number of decomposed sequences, and u k (t) represents the kth sequence function, ω k represents the center frequency of the kth sequence function, t represents time, and δ(t) is the pulse function; In order to solve the above variational problem, the Lagrangian multiplier λ(t) is introduced to obtain the augmented Lagrangian function L: Where α represents the quadratic penalty parameter, which is used to balance the fidelity of the sequence and the compactness of the decomposition. L(u k ,ω k ,λ) represents the augmented Lagrangian function; represents the bandwidth limitation term in VMD; f(t) represents the input power sequence; represents the regularization term; The saddle point of the above augmented Lagrangian function is solved by the alternating direction multiplier method, that is, each power subsequence u is found k and the center frequency ω k ; For each iteration step, update u k and ω k The formula is as follows: In the formula, n represents the number of iterations; and Represent f(t), u i (t), λ(t), and Fourier transform of ; ω represents frequency, represents the center frequency obtained in the n+1th iteration; Through iterative updates, until the stopping condition is met, K subsequences u are obtained k and the corresponding center frequency ω k .

7. The wind power prediction method based on improved LSTM and FA-KELM according to claim 4 is characterized by: Step 2.3 includes the following steps: S2.3.1: Randomly generate particle positions and velocities. The position X and velocity V of each particle represent K,α; S2.3.2: For each particle, use its position parameters to perform VMD decomposition and calculate the reconstruction error as the fitness value; the reconstruction error is defined as the sum of IMF1x(t) after the first decomposition and all the u k The difference between the sum of (t): Where, Reconstruction Error represents the reconstruction error; S2.3.3: Update the position and velocity of the particle according to the velocity update formula and the position update formula; Position update formula: Where, X i (t+1) represents the position information of the i-th particle in the t+1 generation, X i (t) represents the position information of the i-th particle in the t-th generation, V i (t+1) represents the velocity information of the i-th particle in the t+1 generation; Speed ​​update formula: Where σ is the inertia weight, c1 and c2 are acceleration factors, r1 and r2 are random numbers with values ​​between 0 and 1, and P i best is the optimal historical position of particle i itself, G best is the best historical position of all particles, V i (t+1) represents the velocity information of the i-th particle in the t+1 generation, V i (t) represents the velocity information of the i-th particle in the t-th generation, X i (t) represents the position information of the i-th particle in the t-th generation; S2.3.4: Update the historical optimal position P of each particle according to the fitness value i best and the global optimal position G among all particles best ; S2.3.5: Repeat the above steps until the stopping condition is met. The stopping condition is that the maximum number of iterations is reached or the error converges to a certain threshold.

8. The wind power prediction method based on improved LSTM and FA-KELM according to claim 1 is characterized by: In step 3, a long short-term memory network ConvsLSTM combined with convolution simplification is proposed: A gated structure with a cross-coupling structure is used to replace the input gate, forget gate, and output gate in the LSTM network. At the same time, a peephole structure is introduced to consider the impact of the storage unit on the network output. The forward propagation process of the ConvsLSTM structure is shown in the following formula: i t =σ(W x *x t +W h *h t-1 +W c *c t-1 +b i ) o t =σ(W x *x t +W h *h t-1 +W c *c t +b o ) In the formula, * is the convolution operation, is the Hadamard product operation, W x is the input weight, W h is the cycle weight, W c is the peephole weight, b i 、b c and b o is the input bias, σ(·) is the activation function, tanh(·) is the hyperbolic activation function, i t 、c t , o t and h t They represent the input gate, memory cell state, output gate and network output of LSTM respectively. t-1 represents the network output at the previous time, c t-1 Indicates the state of the memory unit at the previous time, x t Represents the input of the network.

9. The wind power prediction method based on improved LSTM and FA-KELM according to claim 1 is characterized by: In step 4, the calculation formula of the maximum correlation coefficient MIC is expressed as: Where X and Y are two variables, a×b represents the grid division method, that is, X is divided into a intervals, Y is divided into b intervals, and the value of a×b is less than a given upper limit B to avoid overfitting, I(X; Y| G ) represents the mutual information between X and Y under a given grid G ​​partitioning, log2min(a,b) is the complexity measure of the grid partitioning, which is used to normalize the mutual information.

10. The wind power prediction method based on improved LSTM and FA-KELM according to claim 1 is characterized by: The step 5 comprises the following steps: Step S5.1: KELM introduces a kernel function based on ELM and replaces the random mapping in ELM with a stable kernel mapping. The algorithm steps are as follows: The output function of the ELM network with L hidden neurons is expressed as: In the formula, f elm (x) represents the output function, x is the input vector, x j represents the specific input value of each input vector, β i is the weight connecting the hidden neuron and the output neuron, g(i) is the activation function, ω i is the randomly assigned weight vector connecting the input nodes and hidden neurons, b i represents the bias of the random assignment of hidden layer neurons, j represents the ordinal number of the value in the input vector, and N represents the number of values ​​in the input vector; The above formula can be redefined as: Hβ=T In the formula, H represents the output matrix of the hidden layer, T represents the predicted expected output matrix, β is the weight connecting the hidden neurons and the output neurons; x1 represents the first number of the input vector, ω1 is the first weight vector, b1 represents the deviation of the random assignment of the first hidden layer neuron, ω L represents the Lth weight vector, b L represents the bias of the random assignment of the Lth hidden layer neuron, x N Represents the Nth number of the input vector; In ELM, the above equation is regarded as a linear function, and its least squares solution is: β=H + T In the formula, H + is the generalized inverse matrix of H; H + =H T (HH T ) -1 In the formula, H T represents the transposed matrix of H, HH T represents a symmetric matrix; The regularization coefficient and the identity matrix are introduced to solve the least squares solution of the hidden layer node output weight β: In the formula, C is the regularization coefficient, I represents the N-order unit matrix; KELM uses the kernel function matrix Q ELM The random matrix HH that replaces ELM T , the present invention adopts RBF function as kernel function: Where σ represents the kernel parameter; K(x i ,x j ) represents the kernel function, x i and x j Represents the i-th and j-th data points in the data set; Combination of the above and Hβ=T, we get the output expression of KELM: Where f(x) represents the output of KELM; K(x,x1) represents the kernel function used to measure the similarity between the new sample x and the first sample x1 in the training set; K(x,x N ) indicates that the kernel function is used to measure the new sample x and the Nth sample x in the training set. N The similarity between Step S5.2: Using the firefly algorithm (FA) to optimize the penalty coefficient C and kernel parameter σ of KELM; First, initialize N fireflies, and the position of each firefly is x i =(σ i ,C i ), i=1,2,…,N,σ i represents the i-th kernel parameter, C i represents the i-th penalty coefficient, N represents the number of fireflies, x i Represents a solution of KELM, and defines the mean square error (MSE) of the KELM model on the validation set as the objective function: In the formula, f(x i ) represents the mean square error, y j is the actual value, is the predicted value of KELM, and n represents the length of the input sequence; The brightness of each firefly I i =-f(x i ), the brighter fireflies attract the brighter fireflies, and the position of the brighter fireflies is x j , the position of the firefly with low brightness is x i , the formula for updating the position is: In the formula, r ij =x i -x j is the distance between fireflies i and j, β0 is the maximum attraction coefficient, γ is the light intensity absorption coefficient, rand is a random number between [0,1], and α is the step size coefficient; After several iterations, find the firefly x with the highest brightness best =(σ best ,C best ), σ best represents the optimal kernel parameter, C best Represents the optimal penalty coefficient. When the maximum number of iterations is reached or the brightness change is less than the preset threshold, the optimization algorithm ends and the optimal parameters of KELM are obtained.

Images