Short-term power load prediction method based on modal aggregation and optimal integration

Abstract

The invention relates to a short-term power load prediction method and system based on modal aggregation and optimal integration. The method specifically comprises the following steps: S1, acquiring power load data in a power system; step S2, carrying out preliminary decomposition on the load sequence based on ICEEMDAN; s3, quantifying the complexity of the subsequences by adopting permutation entropy, and reconstructing the subsequences with high entropy similarity through K-medoids clustering to obtain an aggregation mode component; step S4, carrying out further noise reduction on the high-frequency noise component in the aggregation mode component by adopting a Savitzky-Golay filter; s5, adopting BiGRU as a predictor, and respectively training and predicting the aggregation mode components; s6, constructing a weight optimization model by taking prediction precision and stability as targets; step S7, adopting NSDBO to solve the optimization model, and obtaining a Pareto leading edge; s8, TOPSIS is used as a decision-making method, the relative close degree of each Pareto optimal solution is calculated, and an optimal weight scheme is determined; and S9, performing weighted stacking on the aggregation modal prediction sequence based on a decision scheme to obtain a final short-term power load prediction result.

Description

Technical Field

[0001] The present invention relates to the technical field of power load forecasting, and in particular to a short-term power load forecasting method based on modal aggregation and optimal integration. Background Art

[0002] The rapid development of new power systems and the large-scale adoption of clean energy sources such as wind and solar power have significantly increased the volatility and uncertainty of power systems. Accurate load forecasting provides crucial data support for power system operation planning and demand response. However, existing short-term load forecasting methods lack effective processing of raw load data and also lack consideration of model prediction stability.

[0003] Patent publication number CN110222877A proposes a load forecasting system and method based on a custom neural network. This method introduces a custom neural network, allowing users to input parameter changes, which the system then automatically executes, constructing a new neural network. This makes it more convenient and efficient to use. However, this technology exhibits strong nonlinear and non-stationary characteristics when dealing with raw power loads, and lacks effective load sequence preprocessing, resulting in forecasting accuracy that falls short of expectations.

[0004] Patent publication number CN117494889A proposes a method for establishing a short-term power load forecasting model. This method uses XGBoost as the forecasting model and optimizes model parameters through the OMA algorithm to improve the model prediction accuracy, but ignores the consideration of forecast stability.

[0005] In order to overcome the technical defects of the above technologies, the applicant proposed a short-term power load forecasting method based on modal aggregation and optimal integration. Summary of the Invention

[0006] The present invention solves the technical problems that the existing short-term power load forecasting technology cannot effectively process the original data set, does not fully utilize the differences of subsequences, and only considers the prediction accuracy during the prediction process, which makes the prediction stability poor. The present invention proposes a short-term power load forecasting method based on modal aggregation and optimal integration.

[0007] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0008] The short-term power load forecasting method based on modal aggregation and optimal integration includes the following steps:

[0009] Step S1: Obtain power load data D1 in the target power system;

[0010] Step S2: performing a preliminary decomposition on the acquired power load data D1 based on the improved complete set empirical mode decomposition to obtain a set of intrinsic mode functions IMFs;

[0011] Step S3: permutation entropy is used to quantify the complexity of the intrinsic mode function (IMF), and the intrinsic mode functions with high entropy similarity are reconstructed through K-medoids clustering to obtain a set of aggregated modal components (AIMF).

[0012] Step S4: Define the one with the largest entropy value in the aggregated modal component AIMF as the high-frequency aggregated modal component HAIMF, and the rest as the low-frequency aggregated modal component LAIMF. Use Savitzky-Golay filter to further reduce the noise of HAIMF to obtain the denoised HAIMF;

[0013] Step S5: Taking the low-frequency aggregate modal component LAIMF and the denoised HAIMF as input, respectively, and using BiGRU as the predictor, respectively, predictions are performed to obtain the prediction results of LAIMF and HAIMF;

[0014] Step S6: Taking prediction accuracy and stability as the objective function, the prediction results of LAIMF and HAIMF are used as input to construct an aggregate modal weight optimization model;

[0015] Step S7: using the NSDBO algorithm to solve the aggregate modal weight optimization model and obtain the Pareto frontier;

[0016] Step S8: Using TOPSIS as a decision-making method, calculate the relative closeness of each Pareto frontier optimal solution and determine the optimal weighting scheme;

[0017] Step S9: Based on the optimal weighting scheme, the prediction results of LAIMF and HAIMF are weighted and superimposed to obtain the final short-term power load forecast result.

[0018] In step S2, the specific steps of the improved complete set empirical mode decomposition are:

[0019] Step 1: Construct a new decomposition sequence:

[0020] x new =x+βE1(ω)

[0021] Step 2: x new Use envelope calculation to obtain the first residual;

[0022] r1= <M(x new )>

[0023] Step 3: Calculate the first modal component based on the results of step 2;

[0024] IMF1=x-r1

[0025] Step 4: further calculate the second residual and the second modal component;

[0026] r2=<M(r1+βE2(ω))>

[0027] IMF2=r1-r2

[0028] Step 5: Repeat the above steps to obtain the kth residual and modal component;

[0029] r k = <M(r k-1 +βE k (ω))>

[0030] IMF k =r k-1 -r k

[0031] Step 6: Combine IMF1 to IMF k , and obtain a set of intrinsic mode functions IMFi corresponding to the original load data;

[0032] Among them, x new is the new sequence constructed; rk is the kth residual component; <·> is the mean value calculation; M(·) is the envelope calculation; E k (·) is the kth modal component of EMD decomposition; x is the original sequence; ω is Gaussian white noise; β is the noise coefficient; IMF k is the kth eigenmode function.

[0033] In step S3, the complexity of each intrinsic mode function IMFi is quantified based on the permutation entropy, as follows:

[0034] Step 1. Construct the phase space matrix based on IMF1. The results are as follows:

[0035]

[0036] Step 2. Define each row in M ​​as a component, sort the elements in the component in ascending order, and obtain the symbol sequence S(j) consisting of the position column index j of each element:

[0037] S(j)=[j1,j2,…j m ]

[0038] Step 3: Take the ratio of the number of times each symbol sequence appears to the total number of times different symbol sequences appear as the probability of the symbol sequence, and calculate the probability of each symbol sequence [P1, P2, ..., P Z ], the permutation entropy is calculated as follows:

[0039]

[0040] Step 4. To facilitate analysis, normalize the entropy value to obtain the permutation entropy of IMF1:

[0041]

[0042] Step 5. Repeat Step 1 to Step 4 to calculate IMF1 to IMF2 respectively. k The permutation entropy of .

[0043] Where: M is the phase space matrix; τ represents the delay of the time series; m represents the embedding dimension; N represents the length of the time series, and H n represents the entropy normalization result, S(j) is the symbol sequence, j is the column index, P i is the symbol probability, H p is the permutation entropy, H IMF is the normalized value of permutation entropy.

[0044] In step S3, the intrinsic mode functions IMFi with high entropy similarity are reconstructed based on K-medoids to obtain the aggregated modal components AIMFi.

[0045] Step 1. Select the number of clusters and define the initial center H according to the number of clusters IMF ;

[0046] Step 2, calculate each H according to the Euclidean distance IMF With all centers H IMF The distance between IMF Assigned to the nearest center H IMF ;

[0047] Step 3, calculate all H in each cluster IMF To other H IMF The total distance of IMF As the new center point;

[0048] Step 4. Repeat the distribution of sample H IMF and Update Center H IMF , until the center H IMF No more changes occur, and the IMF reconstruction result AIMF is obtained.

[0049] In step S4, the entropy value of the aggregated modal component AIMFi is quantized again, and the high-frequency noise component HAIMF with the largest entropy value is further denoised using a Savitzky-Golay filter, specifically:

[0050] Step 1: Fix the size of the sliding window and perform polynomial fitting on the window data in HAIMF:

[0051]

[0052] Where: n is the fitting order, q a are polynomial coefficients; L c is the fitted value, and a is the degree of the polynomial.

[0053] Step 2. Calculate the least squares fitting residuals and solve the polynomial coefficients:

[0054]

[0055] Where: 2l+1 is the filter window size; L c is the fitting value; g c is the actual value; c is the index; ε is the least squares fitting residual;

[0056] Step 3, based on the polynomial coefficient q a Reconstruct the polynomial and obtain the denoised sequence HAIMF.

[0057] In step S5, BiGRU is used as the predictor, and GRU is specifically:

[0058] z t =σ(w z [h t-1 ,e t ])

[0059] r t =σ(w r [h t-1 ,e t ])

[0060] h t =tanh(w h [r t h t-1 ,e t ])

[0061] h t =z t h t +(1-z t )h t-1

[0062] Where r t is the reset gate; z t is the update gate; h t-1 is the hidden state at time t-1; h t is the hidden state at time t; is the candidate hidden state at time t; et is the input at time t; w z 、w r and w h Represent the weight coefficients of the update gate, reset gate and gate control unit respectively; σ and tanh represent the sigmoid function and hyperbolic tangent function respectively;

[0063] The specific principles of BiGRU are as follows:

[0064] h f ,t=GRU(e t ,h f,t-1 )

[0065] h b,t =GRU(e t ,h b,t-1 )

[0066] o t =f(w o [h f,t, h b,t ]+b o )

[0067] Where: h f,t and h b,t Represent the forward and backward hidden layer outputs at time t; h f,t-1 and h b,t-1 Represent the forward and backward hidden layer outputs at time t-1 respectively; o t is the neuron output at time t; w o represents the weight coefficient; b o represents bias; f represents activation function; e t is the input at time t; GRU stands for gated recurrent unit;

[0068] Based on BiGRU, the prediction results of LAIMF and HAIMF are obtained respectively.

[0069] In step S6, the aggregate modal weight optimization model takes prediction accuracy and stability as the goal, and the input is the prediction results of LAIMF and HAIMF;

[0070] The objective function of the aggregate modal weight optimization model is:

[0071]

[0072] The constraints of the aggregate modal weight optimization model are:

[0073]

[0074] Where: M APE is the mean absolute percentage error; S Dis the standard deviation; N is the length of the time series, Y i p is the load forecast value; Y i a is the actual load value; E is the mean load forecast error; w j is the weight coefficient; K is the modal aggregation number; is the prediction result of the j-th aggregated modal component at the i-th moment, specifically K-1 LAIMF prediction results and 1 HAIMF prediction result.

[0075] In step S7, the non-dominated sorting dung beetle optimization algorithm is used to solve the optimization model, specifically:

[0076] Step 1: Initialize the dung beetle population. Each dung beetle individual represents a weight coefficient {w1,w2,…,w K} combination, and calculate M APE and S D ;

[0077] Step 2, according to M APE and S D Perform non-dominated sorting, aggregate the solutions into different levels, and sort the individuals in the same level according to the crowding distance;

[0078] Step 3: Update the individual location information of the dung beetle according to its rolling, reproduction, foraging and stealing behaviors, and calculate the M of the new individual. APE and S D , perform non-dominated sorting and crowding distance calculation again;

[0079] Step 4: Repeat Step 3 until the maximum number of iterations is reached to obtain the Pareto optimal solution representing the weight coefficient combination.

[0080] In step S8, TOPSIS is used as a decision-making method to calculate the relative closeness of each Pareto optimal solution and determine the optimal weighting scheme, specifically:

[0081] Step 1: For each set of weight coefficients, construct the original judgment matrix G:

[0082]

[0083] Step 2: Normalize the matrix G to obtain the weighted standard decision matrix U

[0084] U=[u ij ] m×n I

[0085] in

[0086]

[0087] Step 3: Based on the weighted standard decision matrix U, obtain the positive and negative ideal solutions:

[0088]

[0089] Step 4: Calculate the closeness of each decision plan to the positive and negative ideal solutions:

[0090]

[0091] Step 5: Further calculate the relative closeness S i , and the relative closeness S i Sorting, the maximum value corresponds to the optimal solution, that is, the best weight coefficient, which is calculated as follows:

[0092]

[0093] Where: g ij represents the jth indicator of the i-th solution; G is the original judgment matrix; m is the number of solutions; n is the number of evaluation indicators; U is the weighted standard decision matrix; u ij is the weighted standard decision matrix element; I is the identity matrix; and are positive and negative ideal solutions; and are the closeness of positive and negative ideal solutions respectively; S i For relative closeness.

[0094] In step S9, the final load is calculated according to the weight coefficient corresponding to the maximum relative closeness, as follows:

[0095]

[0096] Where: Y load is the final load forecast result; w j is the weight coefficient; K is the modal aggregation number; is the prediction result of the j-th aggregated modal component, specifically K-1 LAIMF prediction results and 1 HAIMF prediction result.

[0097] Compared with the prior art, the present invention has the following technical effects:

[0098] 1) This paper combines secondary signal noise reduction technology, deep learning models, optimization algorithms, and decision-making techniques to propose a short-term load forecasting method. This method reduces the complexity of the original power load sequence through effective data preprocessing techniques, thereby improving load forecasting accuracy. Furthermore, by incorporating forecast stability as another performance evaluation objective, the generalizability of the resulting load forecasting method is greatly enhanced.

[0099] 2) The present invention reduces the complexity of the original sequence by using a data preprocessing method consisting of an improved complete set empirical mode decomposition, permutation entropy, K-medoids clustering and Savitzky-Golay filter, laying a solid foundation for effective short-term load forecasting;

[0100] 3) The present invention constructs a weight optimization model by considering the prediction stability, introduces the non-dominated sorting dung beetle optimization algorithm for super-weight optimization, and adopts TPOSIS as the decision method, so that the differences between modalities are fully utilized, which is conducive to improving the overall prediction effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0101] The present invention will be further described below with reference to the accompanying drawings and examples:

[0102] Figure 1 It is a prediction flow chart of the present invention. DETAILED DESCRIPTION

[0103] A short-term power load forecasting method based on modal aggregation and optimal integration includes the following steps:

[0104] Step S1: Obtain the power load data D1 in the target power system; for example, the power load data of the Belgian power system is selected as the analysis object, the time is from 00:00 on November 1, 2023 to 16:45 on December 31, 2023, the sampling frequency is 15 minutes, and the total number of samples is 5828;

[0105] Step S2: performing a preliminary decomposition on the acquired power load data D1 based on the improved complete set empirical mode decomposition to obtain a set of intrinsic mode functions IMFs;

[0106] Step S3: permutation entropy is used to quantify the complexity of the intrinsic mode function (IMF), and the intrinsic mode functions with high entropy similarity are reconstructed through K-medoids clustering to obtain a set of aggregated modal components (AIMF).

[0107] Step S4: Define the one with the largest entropy value in the aggregated modal component AIMF as the high-frequency aggregated modal component HAIMF, and the rest as the low-frequency aggregated modal component LAIMF. Use Savitzky-Golay filter to further reduce the noise of HAIMF to obtain the denoised HAIMF;

[0108] Step S5: Taking the low-frequency aggregate modal component LAIMF and the denoised HAIMF as input, respectively, and using BiGRU as the predictor, respectively, predictions are performed to obtain the prediction results of LAIMF and HAIMF;

[0109] Step S6: Taking prediction accuracy and stability as the objective function, the prediction results of LAIMF and HAIMF are used as input to construct an aggregate modal weight optimization model;

[0110] Step S7: using the NSDBO algorithm to solve the aggregate modal weight optimization model and obtain the Pareto frontier;

[0111] Step S8: Using TOPSIS as a decision-making method, calculate the relative closeness of each Pareto frontier optimal solution and determine the optimal weighting scheme;

[0112] Step S9: Based on the optimal weighting scheme, the prediction results of LAIMF and HAIMF are weighted and superimposed to obtain the final short-term power load forecast result.

[0113] In step S2, the specific steps of the improved complete set empirical mode decomposition are:

[0114] Step 1: Construct a new decomposition sequence:

[0115] x new =x+βE1(ω)

[0116] Step 2: x new Use envelope calculation to obtain the first residual;

[0117] r1= <M(x new )>

[0118] Step 3: Calculate the first modal component based on the results of step 2;

[0119] IMF1=x-r1

[0120] Step 4: further calculate the second residual and the second modal component;

[0121] r2=<M(r1+βE2(ω))>

[0122] IMF2=r1-r2

[0123] Step 5: Repeat the above steps to obtain the kth residual and modal component;

[0124] r k = <M(r k-1 +βE k (ω))>

[0125] IMF k =r k-1 -r k

[0126] Step 6: Combine IMF1 to IMF k, and obtain a set of intrinsic mode functions IMF corresponding to the original load data;

[0127] Among them, x new is the new sequence constructed; r k is the kth residual component; <·> is the mean value calculation; M(·) is the envelope calculation; E k (·) is the kth modal component of EMD decomposition; x is the original sequence; ω is Gaussian white noise; β is the noise coefficient; IMF k is the kth eigenmode function.

[0128] In step S3, the complexity of each intrinsic mode function IMF is quantified based on the permutation entropy, as follows:

[0129] Step 1. Construct the phase space matrix based on IMF1. The results are as follows:

[0130]

[0131] Step 2. Define each row in M ​​as a component, sort the elements in the component in ascending order, and obtain the symbol sequence S(j) consisting of the position column index j of each element:

[0132] S(j)=[j1,j2,…j m ]

[0133] Step 3: Take the ratio of the number of times each symbol sequence appears to the total number of times different symbol sequences appear as the probability of the symbol sequence, and calculate the probability of each symbol sequence [P1, P2, ..., P Z ], the permutation entropy is calculated as follows:

[0134]

[0135] Step 4. To facilitate analysis, normalize the entropy value to obtain the permutation entropy of IMF1:

[0136]

[0137] Step 5. Repeat Step 1 to Step 4 to calculate IMF1 to IMF2 respectively. k The permutation entropy of .

[0138] Where: M is the phase space matrix; τ represents the delay of the time series; m represents the embedding dimension; N represents the length of the time series, and H n represents the entropy normalization result, S(j) is the symbol sequence, j is the column index, P i is the symbol probability, H p is the permutation entropy, H IMF is the normalized value of permutation entropy.

[0139] In step S3, the intrinsic mode functions (IMFs) with high entropy similarity are reconstructed using K-medoids to obtain the aggregated modal components (AIMFs).

[0140] Step 1. Select the number of clusters and define the initial center H according to the number of clusters IMF ;

[0141] Step 2, calculate each H according to the Euclidean distance IMF With all centers H IMF The distance between IMF Assigned to the nearest center H IMF ;

[0142] Step 3, calculate all H in each cluster IMF To other H IMF The total distance of IMF As the new center point;

[0143] Step 4. Repeat the distribution of sample H IMF and Update Center H IMF , until the center H IMF No more changes occur, and the IMF reconstruction result AIMF is obtained.

[0144] In step S4, the entropy value of the aggregated modal component AIMF is quantized again, and the high-frequency noise component HAIMF with the largest entropy value is further denoised using a Savitzky-Golay filter, specifically:

[0145] Step 1: Fix the size of the sliding window and perform polynomial fitting on the window data in HAIMF:

[0146]

[0147] Where: 2l+1 is the filter window size; L c is the fitting value; g c is the actual value; c is the index; ε is the least squares fitting residual;

[0148] Step 2. Calculate the least squares fitting residuals and solve the polynomial coefficients:

[0149]

[0150] Where: 2l+1 is the filter window size; L c is the fitting value; g c is the actual value; c is the index;

[0151] Step 3, based on the polynomial coefficient q aReconstruct the polynomial and obtain the denoised sequence HAIMF.

[0152] In step S5, BiGRU is used as the predictor, and GRU is specifically:

[0153] z t =σ(w z [h t-1 ,e t ])

[0154] r t =σ(w r [h t-1 ,e t ])

[0155] h t =tanh(w h [r t h t-1 ,e t ])

[0156] h t =z t h t +(1-z t )h t-1

[0157] Where r t is the reset gate; z t is the update gate; h t-1 is the hidden state at time t-1; h t is the hidden state at time t; is the candidate hidden state at time t; e t is the input at time t; w z 、w r and w h Represent the weight coefficients of the update gate, reset gate and gate control unit respectively; σ and tanh represent the sigmoid function and hyperbolic tangent function respectively;

[0158] The specific principles of BiGRU are as follows:

[0159] h f,t =GRU(e t ,h f,t-1 )

[0160] h b,t =GRU(e t ,h b,t-1 )

[0161] o t =f(w o [h f,t, h b,t ]+bo )

[0162] Where: h f,t and h b,t Represent the forward and backward hidden layer outputs at time t; h f,t-1 and h b,t-1 Represent the forward and backward hidden layer outputs at time t-1 respectively; o t is the neuron output at time t; w o represents the weight coefficient; b o represents bias; f represents activation function; e t is the input at time t; GRU stands for gated recurrent unit;

[0163] Based on BiGRU, the prediction results of LAIMF and HAIMF are obtained respectively.

[0164] In step S6, the aggregate modal weight optimization model takes prediction accuracy and stability as the goal, and the input is the prediction results of LAIMF and HAIMF;

[0165] The objective function of the aggregate modal weight optimization model is:

[0166]

[0167] The constraints of the aggregate modal weight optimization model are:

[0168]

[0169] Where: M APE is the mean absolute percentage error; S D is the standard deviation; N is the length of the time series, Y i p is the load forecast value; Y i a is the actual load value; E is the mean load forecast error; w j is the weight coefficient; K is the modal aggregation number; is the prediction result of the j-th aggregated modal component at the i-th moment, specifically K-1 LAIMF prediction results and 1 HAIMF prediction result.

[0170] In step S7, the non-dominated sorting dung beetle optimization algorithm is used to solve the optimization model, specifically:

[0171] Step 1: Initialize the dung beetle population. Each dung beetle individual represents a weight coefficient {w1,w2,…,w K} combination, and calculate M APE and S D ;

[0172] Step 2, according to MAPE and S D Perform non-dominated sorting, aggregate the solutions into different levels, and sort the individuals in the same level according to the crowding distance;

[0173] Step 3: Update the individual location information of the dung beetle according to its rolling, reproduction, foraging and stealing behaviors, and calculate the M of the new individual. APE and S D , perform non-dominated sorting and crowding distance calculation again;

[0174] Step 4: Repeat Step 3 until the maximum number of iterations is reached to obtain the Pareto optimal solution representing the weight coefficient combination.

[0175] In step S8, TOPSIS is used as a decision-making method to calculate the relative closeness of each Pareto optimal solution and determine the optimal weighting scheme, specifically:

[0176] Step 1: For each set of weight coefficients, construct the original judgment matrix G:

[0177]

[0178] Step 2: Normalize the matrix G to obtain the weighted standard decision matrix U

[0179] U=[u ij ] m×n I

[0180] in

[0181]

[0182] Step 3: Based on the weighted standard decision matrix U, obtain the positive and negative ideal solutions:

[0183]

[0184] Step 4: Calculate the closeness of each decision plan to the positive and negative ideal solutions:

[0185]

[0186] Step 5: Further calculate the relative closeness S i , and the relative closeness S i Sorting, the maximum value corresponds to the optimal solution, that is, the best weight coefficient, which is calculated as follows:

[0187]

[0188] Where: g ijrepresents the jth indicator of the i-th solution; G is the original judgment matrix; m is the number of solutions; n is the number of evaluation indicators; U is the weighted standard decision matrix; u ij is the weighted standard decision matrix element; I is the identity matrix; and are positive and negative ideal solutions; and are the closeness of positive and negative ideal solutions respectively; S i For relative closeness.

[0189] In step S9, the final load is calculated according to the weight coefficient corresponding to the maximum relative closeness, as follows:

[0190]

[0191] Where: Y load is the final load forecast result; w j is the weight coefficient; K is the modal aggregation number; is the prediction result of the j-th aggregated modal component, specifically K-1 LAIMF prediction results and 1 HAIMF prediction result.

[0192] Table 1 is a table of prediction results of ablation experiments using the prediction model of the present invention;

[0193] Table 1

[0194]

[0195] As shown in Table 1, the unpredictability of the load sequence was initially reduced by the improved complete ensemble empirical mode decomposition, with its MAPE and SD reduced by 0.345% and 1.556MW respectively. Cluster analysis based on IB can effectively reduce the prediction error caused by multiple modeling of modal components, and the prediction stability is significantly enhanced. SGF is further used to reduce the unpredictability of high-frequency noise components, and the MAPE is relatively reduced by 4.743%. Compared with the average integration and optimal integration strategies, the former directly superimposes the prediction results of each aggregated mode with the same weight, fails to fully utilize the differences between the aggregated modes, resulting in a large prediction error, and its prediction result measurement index R 2 , RMSE, MAPE and SD are 0.990, 99.505MW, 1.406% and 96.735MW respectively. The optimal integration strategy explores the importance of the aggregated mode to the load forecast results and flexibly sets the weight coefficient, thereby effectively utilizing the differences in the aggregated modes. The measurement indicators R2, RMSE, MAPE and SD are 0.992, 85.744MW, 1.276% and 82.369 respectively, which are relatively improved by 0.202%, 13.829%, 9.246% and 14.851%.

[0196] Table 2 is a table showing the comparative experimental results of the prediction model of the present invention:

[0197] Table 2

[0198]

[0199] Comparing the two machine learning models, due to the lack of adaptability of SVM in the face of long time series data, its prediction accuracy and stability are lower than RF; comparing the two deep learning models, the GRU model is slightly better than LSTM in prediction accuracy, but has the problem of low prediction stability; comparing the machine learning model and the deep learning model, the latter has superior prediction accuracy and stability, and performs better in long time series data. Due to the strong non-stationary characteristics of power load, the load is not effectively processed and a single prediction model is used for direct prediction, resulting in a phenomenon of overall large prediction error and low prediction stability. Compared with the single prediction model, the prediction result R of the method proposed in this paper is better than that of the single prediction model. 2 The efficiency increased by 0.031, 0.081, 0.016 and 0.015 respectively; the RMSE decreased by 108.434MW, 206.766MW, 67.309MW and 61.728MW respectively; the MAPE decreased by 1.550%, 2.637%, 1.091% and 0.929% respectively; and the SD decreased by 96.68MW, 185.13MW, 41.384MW and 51.132MW respectively.

Claims

1. A short-term power load forecasting method based on modal aggregation and optimal integration, characterized in that: The following steps are involved: Step S1: Acquire power load data D1 in the target power system; Step S2: preliminarily decompose the acquired power load data D1 based on the improved complete set empirical mode decomposition to obtain a set of intrinsic mode functions IMF; Step S3: permutation entropy is used to quantify the complexity of the intrinsic mode function IMF, and the intrinsic mode functions with high entropy similarity are reconstructed through K-medoids clustering to obtain a set of aggregated modal components AIMF; Step S4: define the aggregate modal component AIMF with the largest entropy value as the high-frequency aggregate modal component HAIMF, and the rest as the low-frequency aggregate modal component LAIMF. Use Savitzky-Golay filter to further reduce the noise of HAIMF to obtain the denoised HAIMF; Step S5: Taking the low-frequency aggregate modal component LAIMF and the denoised HAIMF as inputs, respectively, and using BiGRU as a predictor, respectively, predictions are performed to obtain the prediction results of LAIMF and HAIMF; Step S6: Taking prediction accuracy and stability as the objective function, the prediction results of LAIMF and HAIMF are used as input to construct an aggregate modal weight optimization model; Step S7: using the NSDBO algorithm to solve the aggregation modal weight optimization model and obtain the Pareto frontier; Step S8: Using TOPSIS as a decision-making method, calculate the relative closeness of each Pareto frontier optimal solution and determine the optimal weighting scheme; Step S9: Based on the optimal weighting scheme, the prediction results of LAIMF and HAIMF are weighted and superimposed to obtain the final short-term power load forecast result.

2. The method according to claim 1, characterized in that In step S2, the specific steps of the improved complete set empirical mode decomposition are: Step 2-1: Construct a new decomposition sequence: x new =x+βE1(ω) Step 2-2: For x new Use envelope calculation to obtain the first residual; r1= <M(x new )> Step 2-3: Calculate the first modal component based on the result of step 2-2; IMF1=x-r1 Step 2-4: further calculate the second residual and the second modal component; r2=<M(r1+βE2(ω))> IMF2=r1-r2 Step 2-5: Repeat the above steps to obtain the kth residual and modal component; r k = <M(r k-1 +βE k (oh))> IMF k =r k-1 -r k Step 2-6: Combine IMF1 to IMF k , and obtain a set of intrinsic mode functions IMF corresponding to the original load data; Among them, x new is the new sequence constructed; r k is the kth residual component; <·> is the mean value calculation; M(·) is the envelope calculation; E k (·) is the kth modal component of EMD decomposition; x is the original sequence; ω is Gaussian white noise; β is the noise coefficient; IMF k is the kth eigenmode function.

3. The method according to claim 1, characterized in that In step S3, the complexity of each intrinsic mode function is quantified based on the permutation entropy, as follows: Step 3-1: Construct the phase space matrix based on IMF1. The result is as follows: Step 3-2: Define each row in M ​​as a component, arrange the elements in the component in ascending order, and obtain the symbol sequence S(j) consisting of the position column index j of each element: S(j)=[j1,j2,…j m ] Step 3-3: Take the ratio of the number of occurrences of each symbol sequence to the total number of occurrences of different symbol sequences as the probability of the symbol sequence, and calculate the probability of each symbol sequence [P1, P2, …, P Z ], the permutation entropy is calculated as follows: Step 3-4: For ease of analysis, normalize the entropy value to obtain the permutation entropy of IMF1: Step 3-5: Repeat steps 3-1 to 3-4 to calculate IMF1 to IMF k The permutation entropy of Where: M is the phase space matrix; τ represents the delay of the time series; m represents the embedding dimension; N represents the length of the time series, and H n represents the entropy normalization result, S(j) is the symbol sequence, j is the column index, P i is the symbol probability, H p is the permutation entropy, H IMF is the normalized value of the permutation entropy.

4. The method according to claim 3, characterized in that In step S3, the intrinsic mode function IMF with high entropy value similarity is reconstructed based on K-medoids, as follows: Step s3-1: Select the number of clusters and define the initial center H according to the number of clusters IMF ; Step s3-2: Calculate each H according to the Euclidean distance IMF With all centers H IMF The distance between IMF Assign to the nearest center H IMF ; Step s3-3: Calculate all H in each cluster IMF To other H IMF The total distance of IMF As a new center point; Step s3-4: Repeat the distribution of sample H IMF and Update Center H IMF , until the center H IMF No more changes occur, and the IMF reconstruction result AIMF is obtained.

5. The method according to claim 4, characterized in that In step S4, the entropy value of the aggregated modal component AIMF is quantized again, and the high-frequency noise component HAIMF with the largest entropy value is further denoised using a Savitzky-Golay filter, specifically: Step 4-1: Fix the size of the sliding window and perform polynomial fitting on the window data in HAIMF: Where: n is the fitting order, q a are polynomial coefficients; L c is the fitting value, a is the degree of the polynomial; Step 4-2: Calculate the least squares fitting residuals and solve for the polynomial coefficients: Where: 2l+1 is the filter window size; L c is the fitting value; g c is the actual value; c is the index; ε is the least squares fitting residual; Step 4-3: Based on the polynomial coefficient q a Reconstruct the polynomial and obtain the denoised sequence HAIMF.

6. The method according to any one of claims 1 to 5, characterized in that: In step S5, BiGRU is used as the predictor, and GRU is specifically: With t =σ(in z [h t-1 ,e t ]) r t =σ(w r [h t-1 ,e t ]) h t =tanh(w h [r t h t-1 ,e t ]) h t =z t h t +(1-z t )h t-1 In the formula, r t is the reset gate; t is the update gate; h t-1 is the hidden state at time t-1; h t is the hidden state at time t; is the candidate hidden state at time t; e t is the input at time t; w z 、w r and w h Represent the weight coefficients of the update gate, reset gate and gate control unit respectively; σ and tanh represent the sigmoid function and hyperbolic tangent function respectively; The specific principle of BiGRU is as follows: h f,t =GRU(e t ,h f,t-1 ) h b,t =GRU(e t ,h b,t-1 ) o t =f(w o [h f,t ,h b,t ]+b o ) Where: h f,t and h b,t Respectively represent the forward and backward hidden layer outputs at time t; h f,t-1 and h b,t-1 Respectively represent the forward and backward hidden layer outputs at time t-1; o t is the neuron output at time t; w o represents the weight coefficient; b o represents bias; f represents activation function; e t is the input at time t; GRU stands for Gated Recurrent Unit; Using BiGRU as the predictor, the prediction results of LAIMF and HAIMF are obtained respectively.

7. The method according to claim 6, characterized in that In step S6, the aggregate modal weight optimization model takes prediction accuracy and stability as the goal, and the input is the prediction results of LAIMF and HAIMF; The objective function of the aggregate modal weight optimization model is: The constraints of the aggregate modal weight optimization model are: Where: M APE is the mean absolute percentage error; S D is the standard deviation; N is the length of the time series, is the load forecast value; is the actual load value; E is the mean value of load prediction error; w j is the weight coefficient; K is the modal aggregation number; is the prediction result of the j-th aggregated modal component at the i-th moment, specifically K-1 LAIMF prediction results and 1 HAIMF prediction result.

8. The method according to claim 1 or 7, characterized in that: In step S7, the non-dominated sorting dung beetle optimization algorithm is used to solve the optimization model, specifically: Step 7-1: Initialize the dung beetle population. Each dung beetle individual represents a weight coefficient {w1,w2,…,w K } combination, and calculate M APE and S D ; Step 7-2: According to M APE and S D Perform non-dominated sorting, aggregate the solutions into different levels, and sort the individuals in the same level according to the crowding distance; Step 7-3: Update the individual position information of the dung beetle according to its rolling, reproduction, foraging and stealing behaviors, and calculate the M of the new individual. APE and S D , perform non-dominated sorting and crowding distance calculation again; Step 7-4: Repeat step 7-3 until the maximum number of iterations is reached to obtain the Pareto optimal solution representing the weight coefficient combination.

9. The method according to claim 1, characterized in that: In step S8, TOPSIS is used as a decision-making method to calculate the relative closeness of each Pareto optimal solution and determine the optimal weighting scheme, which is specifically: Step 8-1: For each set of weight coefficients, construct the original judgment matrix G: Step 8-2: Normalize the matrix G to obtain the weighted standard decision matrix U U=[u ij ] m×n AND in Step 8-3: Based on the weighted standard decision matrix U, obtain the positive and negative ideal solutions: Step 8-4: Calculate the closeness of each decision plan to the positive and negative ideal solutions: Step 8-5: Further calculate the relative closeness S i , and the relative closeness S i Sorting, the maximum value corresponds to the optimal solution, that is, the best weight coefficient, which is calculated as follows: Where: g ij represents the jth indicator of the ith solution; G is the original judgment matrix; m is the number of solutions; n is the number of evaluation indicators; U is the weighted standard decision matrix; u ij is the weighted standard decision matrix element; I is the identity matrix; and are positive and negative ideal solutions; and are the closeness of positive and negative ideal solutions respectively; S i For relative closeness.

10. The method according to claim 1, characterized in that In step S9, the final load is calculated according to the weight coefficient corresponding to the maximum relative closeness, as follows: Where: Y load is the final load forecast result; w j is the weight coefficient; K is the modal aggregation number; It is the prediction result of the j-th aggregated modal component, specifically K-1 LAIMF prediction results and 1 HAIMF prediction result.

Images