A short-term load forecasting method for a power system

Abstract

The present invention discloses a short-term load forecasting method for a power system, including: acquiring historical load data of the power system and preprocessing it; optimizing the decomposition mode number and penalty factor of the variational mode decomposition algorithm through the Tianying optimizer algorithm, and decomposing the load data into intrinsic mode function components with different central frequencies, including low-frequency intrinsic mode function components and high-frequency intrinsic mode function components; optimizing the LSSVM prediction models based on linear kernel functions and radial basis kernel functions respectively through the fireworks algorithm, predicting the low-frequency intrinsic mode function components and high-frequency intrinsic mode function components respectively to obtain the prediction values of the low-frequency intrinsic mode function components and the high-frequency intrinsic mode function components, and obtaining the final load forecasting result therefrom. The present invention uses the variational mode decomposition algorithm optimized by the Tianying optimizer algorithm to decompose the load data, and uses the LSSVM prediction model optimized by the fireworks algorithm to predict the load data, thereby improving the load forecasting accuracy.

Description

Technical Field

[0001] The present invention belongs to the technical field of power system load forecasting, and particularly relates to a short-term load forecasting method for power systems. Background Art

[0002] Short-term power system load forecasting is to predict the load power of the power system within one hour to one week, which plays a key role in the economic and safe operation strategies of the power system and in power trading. At present, the forecasting methods applicable to this forecasting modeling process can be divided into two types: statistical-based models and artificial intelligence-based models. Statistical-based models include methods such as multiple regression, ARMA, and ARIMA; artificial intelligence-based models include methods such as SVM, LSSVM, BP neural network, Deep Belief Network (DBN), and Recurrent Neural Network (RNN).

[0003] Since there will be noise in the load data during the acquisition process, which will affect the forecasting accuracy, decomposing the load sequence first, removing the noise, and then predicting each decomposed sequence one by one is a current research hotspot for load forecasting. Currently, the main decomposition methods include VMD decomposition, EMD decomposition, etc. However, when using VMD decomposition, if the decomposition mode number K and penalty factor α parameters in VMD decomposition are set according to empirical methods, it will affect the final load forecasting accuracy. Summary of the Invention

[0004] Object of the Invention: Aiming at the load forecasting accuracy problem caused by setting the K and α parameters in VMD decomposition according to empirical methods in the prior art, the present invention discloses a short-term load forecasting method for power systems, realizing short-term power system load forecasting and achieving the goal of reducing load forecasting errors.

[0005] Technical Solution: To achieve the above object of the invention, the present invention adopts the following technical solution:

[0006] A short-term load forecasting method for power systems, comprising the following steps:

[0007] S1. Obtain the historical load data of the power system and preprocess it;

[0008] S2. Optimize the decomposition mode number and penalty factor of the variational mode decomposition algorithm through the Tianying optimizer algorithm, and decompose the load data into intrinsic mode function components with different center frequencies through the optimized variational mode decomposition algorithm, wherein the intrinsic mode function components are divided into low-frequency intrinsic mode function components and high-frequency intrinsic mode function components;

[0009] S3. Optimize the penalty coefficient of the LSSVM prediction model based on the linear kernel function through the fireworks algorithm, and use the optimized LSSVM prediction model based on the linear kernel function to predict the low-frequency intrinsic mode function components to obtain the predicted values of the low-frequency intrinsic mode function components;

[0010] S4. Optimize the penalty coefficient and the kernel function width factor of the LSSVM prediction model based on the radial basis kernel function through the fireworks algorithm, and use the optimized LSSVM prediction model based on the radial basis kernel function to predict the high-frequency intrinsic mode function components to obtain the predicted values of the high-frequency intrinsic mode function components;

[0011] S5. Add the predicted values of the low-frequency intrinsic mode function components and the predicted values of the high-frequency intrinsic mode function components to obtain the final load prediction result.

[0012] Preferably, in step S2, the general Tianying optimizer algorithm optimizes the decomposition mode number and the penalty factor of the variational mode decomposition algorithm, including the following steps:

[0013] S21: Initialize the algorithm parameters of the Tianying optimizer algorithm, including the population size, the total number of iterations, the range of the decomposition mode number K, and the range of the penalty factor α;

[0014] S22: Initialize the Tianying positions, generate the initial population X, calculate the fitness values, and set the individual with the minimum fitness as the best individual;

[0015] X = (UB1 - LB1) × rand + LB1

[0016] where rand represents a random number taken in the interval [0, 1], UB1 is a vector composed of the upper limits of the ranges of the variable K and the variable α, and LB1 is a vector composed of the lower limits of the ranges of the variable K and the variable α;

[0017] S23: Start the iterative operation. The initial number of iterations is 0, and the number of iterations is incremented by 1 for each iteration. If the number of iterations is less than or equal to 2 / 3 of the total number of iterations, then a random number between [0, 1] is generated. If this number is less than 0.5, then update the positions according to Method 1, otherwise update the positions according to Method 2; if the number of iterations is greater than 2 / 3 of the total number of iterations, then another random number between [0, 1] is generated. If this number is less than 0.5, then update the positions according to Method 3, otherwise update the positions according to Method 4;

[0018] Calculate the fitness values of the population after updating the positions, compare the fitness of the current best individual and the best individual after iteration, and retain the optimal individual;

[0019] S24: Determine whether the number of iterations is less than the total number of iterations. If so, return to step S23 to continue the iteration. If the number of iterations is equal to the total number of iterations, stop the iteration and output the optimal solutions of the decomposition mode number and the penalty factor;

[0020] Among them, the average value of the signal difference is used as the fitness, and the calculation formula is as follows:

[0021]

[0022] Among them, is the sum of the decomposed IMF components, f(n) is the original sequence, and N is the number of discrete points in the sequence;

[0023] The methods for updating the position include:

[0024] Method 1: Vertical dive attack, which is mathematically expressed as:

[0025]

[0026] where t is the current number of iterations, X1(t + 1) is the position of the eagle using Method 1 at the (t + 1)-th iteration, X best (t) is the best position after the t-th iteration, T1 is the maximum number of iterations, X M (t) is the average position of the population after the t-th iteration, and rand1 represents a random number taken in the interval [0, 1];

[0027] Method 2: Contour flight and short glide attack, which is mathematically expressed as:

[0028] X2(t + 1) = X best (t) × Levy + X R (t) + (y - x) × rand2

[0029] where X2(t + 1) is the position of the eagle using Method 2 at the (t + 1)-th iteration, X R (t) is a random one of the N positions at the t-th iteration, N is the total number of eagles in the population, rand2 represents a random number taken in the interval [0, 1], Levy is the Levy flight distribution function, and y and x are in a spiral shape and satisfy:

[0030]

[0031] where r is the spiral radius and θ is the spiral angle, and there is:

[0032] r = r1 + 0.00565 × D1

[0033]

[0034] Among them, r1 is a value between [1, 20], and D1 is any integer from 1 to the length of the search space (Dim);

[0035] Method 3: Low-altitude flight and slow descent attack, which is mathematically expressed as:

[0036] X3(t + 1) = (X best (t) - X M (t)) × 0.1 - rand3 + ((UB1 - LB1) × rand4 + LB1) × 0.1

[0037] Among them, X3(t + 1) is the position of the eagle using Method 3 at the (t + 1)-th iteration, and rand3 and rand4 respectively represent a number randomly taken within the interval [0, 1];

[0038] Method 4: Land walking attack, which is mathematically expressed as:

[0039] X4(t + 1) = QF(t) × X best (t) - G1 × (1 - X(t)) × rand5 - G2 × Levy

[0040] Among them, X4(t + 1) is the position of the eagle using Method 4 at the (t + 1)-th iteration, X(t) is the position at the t-th iteration, rand5 represents a number randomly taken within the interval [0, 1], and QF(t) is the quality function of the balanced search strategy, and there is:

[0041]

[0042] Among them, rand6 represents a number randomly taken within the interval [0, 1];

[0043] G1 is various movements of the eagle during the prey's escape, and there is:

[0044] G1 = 2 × rand7 - 1

[0045] Among them, rand7 represents a number randomly taken within the interval [0, 1];

[0046] G2 is a decreasing value from 2 to 0, representing the flight slope of the eagle following the prey from the first position to the last position during the prey's escape, and there is:

[0047]

[0048] Among them, t is the current iteration number, and T1 is the total number of iterations.

[0049] Preferably, in step S2, the output intrinsic mode function components are divided according to the central frequency. Those with a central frequency less than or equal to the frequency threshold are low-frequency intrinsic mode function components, and those with a central frequency greater than the frequency threshold are high-frequency intrinsic mode function components.

[0050] Preferably, in step S3, the penalty coefficient of the LSSVM prediction model based on the linear kernel function is optimized by the fireworks algorithm, including the following steps:

[0051] S31: Initialize the parameters, including: the number of fireworks P, the maximum number of iterations T2, the number of mutant sparks b, the number of explosions q, the explosion radius h, the upper bound UB2 and the lower bound LB2 of the penalty coefficient;

[0052] S32: Generate the initial fireworks population X2, calculate the fitness of each firework, that is, the prediction error, and select the one with the minimum fitness as the best firework;

[0053] X2 = (UB2 - LB2) × rand + LB2

[0054] where rand is a random number between 0 and 1;

[0055] S33: Enter the iteration, the initial iteration number is 0, and calculate the explosion radius and the number of explosions of each firework according to the following formula:

[0056]

[0057] where f(i) is the fitness of the i-th firework individual, fr(i) is the explosion radius of the i-th firework, f min is the minimum fitness in the fireworks population, f max is the maximum fitness in the fireworks population, fn(i) is the number of sparks generated by the i-th firework, and ε is a machine epsilon;

[0058] S34: Select the d-th variable x id , 1 ≤ d ≤ D, D is the variable dimension, and generate the explosion spark ex id and the Gaussian mutation spark mx id ;

[0059] Check whether the newly generated sparks are out of bounds, and perform mapping processing on the sparks that exceed the boundary. If it exceeds the upper limit, take the upper limit value, and if it exceeds the lower limit, take the lower limit value;

[0060] ex id = x id + fr(i) × (2 × rand(1, E) - 1)

[0061] mx id = x id × e

[0062] Among them, rand(1, E) randomly selects a number within the interval (1, E), where E is the number of optimization variables, which is 1 here, and e is a Gaussian distribution random number;

[0063] S35: Calculate the fitness values of the explosion sparks and Gaussian mutation sparks, compare them with the fitness values of the fireworks, sort all individuals in ascending order of fitness values, select the top P individuals to enter the next iteration, and increment the iteration count by 1;

[0064] S36: Determine whether the current iteration count is less than the maximum iteration count. If so, return to step S33 to continue the iteration; otherwise, terminate the calculation, output the individual with the minimum fitness value and the fitness value, and the individual with the minimum fitness value is the optimal penalty coefficient.

[0065] Beneficial effects: Compared with the prior art, the present invention has the following remarkable beneficial effects:

[0066] 1. The present invention uses the AO algorithm to optimize the VMD algorithm to decompose the load data, obtaining the optimal combination of K and α in the VMD algorithm, laying a foundation for improving the load forecasting accuracy;

[0067] 2. For the IMF components with different frequencies obtained by decomposition, the present invention respectively establishes an LSSVM prediction model based on a linear kernel function and an LSSVM prediction model based on a radial basis kernel function, and optimizes them using the fireworks algorithm, greatly improving the load forecasting accuracy. Brief Description of the Drawings

[0068] Figure 1 is the step flow chart of the load forecasting method described in the present invention;

[0069] Figure 2 is the flow chart of the AO algorithm optimizing the AMD algorithm in the load forecasting method described in the present invention;

[0070] Figure 3 is the flow chart of the fireworks algorithm optimizing the LSSVM prediction model in the load forecasting method described in the present invention;

[0071] Figure 4 is the IMF component curve obtained by decomposing the AMD algorithm optimized by the AO algorithm;

[0072] Figure 5 is the result of load forecasting using the load forecasting method described in the present invention. Detailed Embodiments

[0073] The following further describes the present invention with reference to the accompanying drawings.

[0074] The present invention discloses a short-term load forecasting method for a power system, as Figure 1As shown, it includes the following steps:

[0075] Step S1: Obtain the historical load data of the power system and preprocess it. The preprocessing includes: using the mean imputation method to supplement the missing load data.

[0076] Step S2: Use the Tianying optimizer (AO) algorithm to optimize the variational mode decomposition (VMD) algorithm. Through the optimized VMD algorithm, decompose the load data into intrinsic mode function (IMF) components with different center frequencies.

[0077] The process of optimizing the VMD algorithm by the AO algorithm is as Figure 2 shown.

[0078] The VMD algorithm decomposes a sequence f(n) into K intrinsic mode function components u k , k = 1, 2,..., K. Taking the sum of the estimated bandwidths of all IMF components as the objective function and the sum of all IMF components being equal to f(n) as the constraint condition, and then through the Lagrange transformation, transform the constrained problem into an unconstrained problem. The objective function of the variational constraint model after the Lagrange transformation is as follows:

[0079]

[0080] where, {u k} = {u1, u2,..., u K} is the set of IMF components obtained by decomposition, {ω k} = {ω1, ω2,..., ω K} is the set of center frequencies of the IMF components, ω k is the center frequency of the IMF component u k , λ is the Lagrange multiplier, α is the quadratic penalty factor, K is the total number of IMF components, i.e., the decomposition mode number, δ n is the Dirac distribution, θ n is the partial derivative with respect to n, n represents the time node, and ||·||2 represents the 2-norm.

[0081] If the decomposition mode number K in the VMD algorithm is too large, it will lead to over-decomposition. If K is too small, it will lead to under-decomposition, so that some IMF components cannot be effectively identified; the quadratic penalty factor α will affect the bandwidth of the IMF components. If α is too large, modal aliasing is likely to occur. If α is too small, the information contained in the IMF components will be lacking. Therefore, the selection of these two parameters is crucial for the VMD algorithm. The general selection method is that K is subjectively set by humans, and α adopts the default value of 2000. However, this setting method ignores the interaction between parameters and can only obtain a relatively optimal solution.

[0082] The present invention optimizes the VMD algorithm using the AO algorithm, and uses the average value of the signal difference as the fitness of the AO algorithm for optimization to obtain the best combination of K and α. The calculation formula for the average value of the signal difference is as follows:

[0083]

[0084] where is the sum of the IMF components after decomposition, f(n) is the original sequence, and N is the number of discrete points in the sequence.

[0085] The AO algorithm is an intelligent optimization algorithm proposed based on the hunting method of the eagle. Most eagles will adopt different hunting methods according to different situations, generally including four methods:

[0086] Method 1: Vertical dive attack, which is mathematically expressed as:

[0087]

[0088] where t is the current iteration number, X1(t + 1) is the position of the eagle using Method 1 at the (t + 1)-th iteration, X best (t) is the best position after the t-th iteration, T1 is the maximum number of iterations, X M (t) is the average position of the population after the t-th iteration, and rand1 represents a random number taken in the interval [0, 1].

[0089] Method 2: Contour flight and short glide attack, which is mathematically expressed as:

[0090] X2(t + 1) = X best (t) × Levy + X R (t) + (y - x) × rand2

[0091] where X2(t + 1) is the position of the eagle using Method 2 at the (t + 1)-th iteration, X R (t) is a random one among the N positions after the t-th iteration, N is the total number of eagles in the population, rand2 represents a random number taken in the interval [0, 1], Levy is the Levy flight distribution function, and y and x are in a spiral shape, satisfying:

[0092]

[0093] where r is the spiral radius and θ is the spiral angle, and there is:

[0094] r = r1 + 0.00565 × D1

[0095]

[0096] Among them, r1 is a value between [1, 20], and D1 is any integer from 1 to the search space length (Dim).

[0097] Method 3: Low-altitude flight and slow descent attack, which is mathematically expressed as:

[0098] X3(t + 1) = (X best (t) - X M (t)) × 0.1 - rand3 + ((UB1 - LB1) × rand4 + LB1) × 0.1

[0099] Among them, X3(t + 1) is the position of the eagle using Method 3 at the (t + 1)-th iteration, rand3 and rand4 respectively represent a number randomly taken in the interval [0, 1], UB1 is the upper limit of the optimization variable, and LB1 is the lower limit of the optimization variable.

[0100] Method 4: Land walking attack, which is mathematically expressed as:

[0101] X4(t + 1) = QF(t) × X best (t) - G1 × (1 - X(t)) × rand5 - G2 × Levy

[0102] Among them, X4(t + 1) is the position of the eagle using Method 4 at the (t + 1)-th iteration, X(t) is the position at the t-th iteration, rand5 represents a number randomly taken in the interval [0, 1], QF(t) is the quality function for balancing the search strategy, and there is:

[0103]

[0104] Among them, rand6 represents a number randomly taken in the interval [0, 1];

[0105] G1 is various movements of the eagle during the prey's escape, and there is:

[0106] G1 = 2 × rand7 - 1

[0107] Among them, rand7 represents a number randomly taken in the interval [0, 1];

[0108] G2 is a decreasing value from 2 to 0, representing the flight slope of the eagle following the prey from the first position to the last position during the prey's escape, and there is:

[0109]

[0110] The specific steps for AO to optimize VMD decomposition are as follows:

[0111] Step S21: Initialize the AO algorithm parameters, including the number of AO populations, the total number of iterations, the range of variable K, and the range of variable α.

[0112] In an embodiment of the present invention, the number of AO populations can be set to 100, the total number of iterations to 50, the range of variable K to [2, 20], and the range of variable α to [0, 5000].

[0113] Step S22: Initialize the Tianying position, generate the initial population X, calculate the SDA fitness value, and set the individual with the minimum fitness as the best individual;

[0114] X = (UB1 - LBl) × rand + LBl

[0115] where rand represents a randomly selected number in the interval [0, 1], UB1 is a vector composed of the upper limits of the ranges of variable K and variable α, and LB1 is a vector composed of the lower limits of the ranges of variable K and variable α.

[0116] Step S23: Start the iterative operation. The initial iteration number t = 0, and the iteration number is incremented by 1 for each iteration. If the iteration number is less than or equal to 2 / 3 of the total number of iterations, then a number between [0, 1] is randomly generated. If this number is less than 0.5, the position is updated according to Method 1, otherwise it is updated according to Method 2; if the iteration number is greater than 2 / 3 of the total number of iterations, another number between [0, 1] is randomly generated. If this number is less than 0.5, the position is updated according to Method 3, otherwise it is updated according to Method 4. Calculate the fitness value of the population after updating the position, compare the fitness of the current best individual and the best individual after iteration, and retain the optimal individual.

[0117] Step S24: Determine whether the iteration number is less than the total number of iterations. If so, return to Step S23 to continue the iteration. If the iteration number is equal to the total number of iterations, stop the iteration and output the optimal solutions of K and α.

[0118] The AO algorithm optimizes and outputs the optimal solutions of the decomposition mode number K and the quadratic penalty factor α. The VMD algorithm uses the optimal solutions of the above decomposition mode number K and the quadratic penalty factor α to perform VMD decomposition on the load data and outputs several IMF components. The output several IMF components are divided according to the central frequency. Those with a central frequency less than or equal to the frequency threshold are low-frequency IMF components, and those with a central frequency greater than the frequency threshold are high-frequency IMF components.

[0119] As Figure 4 shown, in an embodiment of the present invention, the VMD algorithm optimized by the AO algorithm (i.e., the VMD algorithm using the optimal solutions of the decomposition mode number K and the quadratic penalty factor α optimized by the AO algorithm) for Figure 4The load data in [[]] was decomposed by VMD, and six IMF components, namely IMF1, IMF2, IMF3, IMF4, IMF5, and IMF6, were obtained. Among them, IMF6 with a divided center frequency less than or equal to 0.01 Hz is the low-frequency IMF component, and IMF1, IMF2, IMF3, IMF4, and IMF5 with a center frequency greater than 0.01 Hz are the low-high-frequency IMF components.

[0120] Step S3: For the decomposed low-frequency IMF component, a prediction model of LSSVM optimized by the fireworks algorithm based on the linear kernel function is used for prediction to obtain the predicted value of the low-frequency IMF component.

[0121] The basic principle of the least squares support vector machine (LSSVM) is as follows:

[0122] For the non-linear load prediction model:

[0123] f(x) = (ω, φ(x)) + b

[0124] where ω is the weight vector, φ(x) is the non-linear mapping from the input space to the high-dimensional feature space, and b is the bias.

[0125] Given a set of data point sets (x i , y i ), i = 1,..., l, x i ∈ R D is the historical load data, D is the dimension of the selected input variable, y i ∈ R is the expected value of the predicted quantity, and l is the total number of known data points. According to the principle of structural minimization, the LSSVM optimization objective can be expressed as

[0126]

[0127] s.t. ω T φ(x i ) + b + e i = y i , i = 1,..., l

[0128] where e i is the error, e = [e1, e2,..., e l ∈ R l×1 is the error vector, γ is the regularization parameter, which controls the degree of punishment for the error. Introduce the Lagrange multiplier λ i ∈ R l×1 , and the above formula can be transformed into:

[0129]

[0130] From the KKT conditions, we get

[0131]

[0132] Eliminating ω and e, the solution of the above equation is:

[0133]

[0134] Among them, K(x, x i ) is the kernel function, x is the D-dimensional input vector, x i is the center of the i-th kernel function, having the same dimension as x, σ is the normalization parameter, which determines the width of the function around the center point, ||x - x i || is the norm of the vector x - x i , representing the distance between x and x i .

[0135] Since the low-frequency IMF component has relatively small fluctuations, the present invention uses a simple and easy-to-operate linear kernel function LSSVM to predict the low-frequency IMF component. The formula of the linear kernel function is:

[0136]

[0137] Therefore, for the LSSVM prediction model based on the linear kernel function, only the penalty coefficient γ needs to be selected as the object optimized by the fireworks algorithm, and the minimum prediction error is taken as the objective function to establish the fireworks algorithm-optimized LSSVM prediction model based on the linear kernel function.

[0138] The objective function, that is, the prediction error formula, is:

[0139]

[0140] Among them, y(m) is the actual load data at the m-th time node of the selected prediction day, y′(m) is the predicted load data at the m-th time node, and M is the total number of time nodes.

[0141] The specific steps of the fireworks algorithm-optimized LSSVM prediction model based on the linear kernel function are as follows:

[0142] Step S31: Initialize the parameters, including: the number of fireworks, the maximum number of iterations, the number of mutated sparks, the number of explosions, the explosion radius, the upper limit of the variable, and the lower limit of the variable;

[0143] For example, set the number of fireworks P to 20, the maximum number of iterations T2 to 50, the number of mutated sparks b to 20, the number of explosions q to 10, the explosion radius h to 20, the upper limit UB2 of the penalty coefficient γ to 500, and the lower limit LB2 of the penalty coefficient γ to 0.1;

[0144] Step S32: Generate the initial fireworks population X2, calculate the fitness of each firework, that is, the prediction error, and select the one with the minimum fitness as the best firework;

[0145] X2 = (UB2 - LB2) × rand + LB2

[0146] Where rand is a random number between 0 and 1;

[0147] Step S33: Enter the iteration. The initial iteration number t = 0. Calculate the explosion radius and the number of explosions of each firework according to the following formulas respectively:

[0148]

[0149] Where f(i) is the fitness of the i-th firework individual, fr(i) is the explosion radius of the i-th firework, f min is the minimum fitness in the firework population, f max is the maximum fitness in the firework population, fn(i) is the number of sparks generated by the i-th firework, ε is a machine minimum used to avoid division by zero, and the minimum value generated each time is different;

[0150] Step S34: Select the d-th variable x id in the i-th firework individual x(i), 1 ≤ d ≤ D. Generate the explosion spark ex id and the Gaussian mutation spark mx id according to the following formula. Check whether the newly generated sparks exceed the limit. Perform mapping processing on the sparks that exceed the boundary. If it exceeds the upper limit, take the upper limit value; if it exceeds the lower limit, take the lower limit value.

[0151] ex id = x id + fr(i) × (2 × rand(1, E) - 1)

[0152] mx id = x id × e

[0153] Where rand(1, E) is a random number randomly selected within the interval (1, E), E is the number of optimization variables, which is 1 here, and e is a random number of Gaussian distribution;

[0154] Step S35: Calculate the fitness values of the explosion sparks and the Gaussian mutation sparks, compare them with the firework fitness values, sort all individuals in ascending order of fitness value, select the top P individuals to enter the next iteration, and increment the iteration number by 1;

[0155] Step S36: Determine whether the current iteration number t is less than the maximum iteration number T. If so, return to Step S33 to continue the iteration; otherwise, terminate the calculation and output the individual with the minimum fitness and the fitness value. The individual with the minimum fitness is the optimal penalty coefficient γ.

[0156] Step S4. For the decomposed high-frequency IMF components, since the data fluctuates greatly and is not easy to fit, a prediction model of LSSVM based on radial basis kernel function optimized by fireworks algorithm is used for prediction to obtain the predicted values of high-frequency IMF components.

[0157] The formula of the radial basis kernel function is:

[0158]

[0159] where σ is the kernel function width factor. For the LSSVM prediction model based on the radial basis kernel function, two parameters, γ and σ, need to be selected 2 as the optimization objects of the fireworks algorithm. Taking the prediction error minimum as the objective function, a prediction model of LSSVM based on the radial basis kernel function optimized by the fireworks algorithm is established. The specific steps are the same as those of the fireworks algorithm for optimizing the LSSVM prediction model based on the linear kernel function in Step S3, except that the upper and lower variable limits in Step S31 are adaptively modified. For example, the upper bound vector UB of the penalty coefficient and the kernel function width factor is set to [500, 10], and the lower bound vector LB of the penalty coefficient and the kernel function width factor is set to [0.1, 0.1].

[0160] Step S5. Add the predicted values of the low-frequency IMF components obtained in Step S3 and the predicted values of the high-frequency IMF components obtained in Step S4 to obtain the final load prediction result.

[0161] Taking the load data of a certain area as an example to verify the effectiveness of the load prediction method described in the present invention. The load of this area in a certain month (31 days) is collected, once every 15 minutes, with a total of 2976 data. The first 30 days are selected as the training set, and the 31st day is the test set. The input is the load of the three days before the prediction day. The final prediction result is as Figure 5 shown. The average absolute error percentage is 2.13%, indicating that the short-term load prediction model established by the present invention has a relatively high prediction accuracy.

[0162] The method described in the present invention can be used not only for short-term load prediction but also for medium- and long-term load prediction: short-term load prediction generally predicts the load of one day or one week for power grid scheduling, and not much training data is required; medium- and long-term load prediction generally predicts on a monthly and annual basis for infrastructure planning, and several years of historical data of this area are required for training, so the data scale is very large. The method described in the present invention needs to decompose the data and then predict one by one. If medium- and long-term load prediction is carried out, the training time is relatively long, and factors such as local economy and population need to be considered.

[0163] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention.

Claims

1. A short-term load forecasting method for a power system, characterized in that, It includes the following steps: S1. Obtain the historical load data of the power system and preprocess it; S2. Optimize the decomposition mode number and penalty factor of the variational mode decomposition algorithm through the Tianying optimizer algorithm. Decompose the load data into intrinsic mode function components with different central frequencies through the optimized variational mode decomposition algorithm. Among them, the intrinsic mode function components are divided into low-frequency intrinsic mode function components and high-frequency intrinsic mode function components; S3. Optimize the penalty coefficient of the LSSVM prediction model based on the linear kernel function through the fireworks algorithm. Predict the low-frequency intrinsic mode function components through the optimized LSSVM prediction model based on the linear kernel function to obtain the predicted values of the low-frequency intrinsic mode function components; S4. Optimize the penalty coefficient and kernel function width factor of the LSSVM prediction model based on the radial basis kernel function through the fireworks algorithm. Predict the high-frequency intrinsic mode function components through the optimized LSSVM prediction model based on the radial basis kernel function to obtain the predicted values of the high-frequency intrinsic mode function components; S5. Add the predicted values of the low-frequency intrinsic mode function components and the predicted values of the high-frequency intrinsic mode function components to obtain the final load prediction result; Among them, in step S3, optimizing the penalty coefficient of the LSSVM prediction model based on the linear kernel function through the fireworks algorithm includes the following steps: S31: Initialize the parameters, including: the number of fireworks P, the maximum number of iterations T2, the number of mutant sparks b, the number of explosions q, the explosion radius h, the upper limit UB2 and the lower limit LB2 of the penalty coefficient; S32: Generate the initial fireworks population X2, calculate the fitness of each firework, that is, the prediction error, and select the one with the minimum fitness as the best firework; X2 = (UB2 - LB2) × rand + LB2 where rand is a random number between 0 and 1; S33: Enter the iteration, the initial iteration number is 0, and calculate the explosion radius and the number of explosions of each firework according to the preset formula respectively; S34: Select the d-th variable xd in the i-th firework individual x(i), where 1 ≤ d ≤ D and D is the variable dimension, and generate the explosion spark ex according to the following formula id , and the Gaussian mutation spark mx id ; id ; Check whether the newly generated sparks exceed the limit, and perform mapping processing on the sparks that exceed the boundary. If it exceeds the upper limit, take the upper limit value, and if it exceeds the lower limit, take the lower limit value; ex id = x id + fr(i) × (2 × rand(1, E) - 1) mx id = x id × e where rand(1, E) is a random number drawn within the interval (1, E), E is the number of optimization variables, here it is 1, and e is a Gaussian distribution random number; S35: Calculate the fitness values of the explosion sparks and the Gaussian mutant sparks, and compare them with the firework fitness values. Arrange all individuals in ascending order of fitness values, select the first P individuals to enter the next iteration, and add 1 to the iteration number; S36: Judge whether the current iteration number is less than the maximum iteration number. If so, return to step S33 to continue the iteration, otherwise terminate the calculation, output the individual with the minimum fitness and the fitness value, and the individual with the minimum fitness is the optimal penalty coefficient.

2. The short-term load forecasting method for a power system according to claim 1, wherein, In step S33, calculate the explosion radius and the number of explosions of each firework according to the preset formula respectively. The preset formula includes: Among them, f(i) is the fitness of the i-th firework individual, fr(i) is the explosion radius of the i-th firework, f min is the minimum fitness in the firework population, f max is the maximum fitness in the firework population, fn(i) is the number of sparks generated by the i-th firework, and ε is a machine epsilon.

3. A short-term load forecasting method for a power system according to claim 1, characterized in that, In step S2, optimizing the decomposition mode number and penalty factor of the variational mode decomposition algorithm through the general Tianying optimizer algorithm includes the following steps: S21: Initialize the algorithm parameters of the Tianying optimizer algorithm, including the population size, the total number of iterations, the range of the decomposition mode number K, and the range of the penalty factor α; S22: Initialize the position of the eagle, generate the initial population X, calculate the fitness value, and set the individual with the minimum fitness as the best individual; X = (UB1 - LB1) × rand + LB1 where rand represents a randomly taken number in the interval [0, 1], UB1 is the vector composed of the upper limit of the range of variable K and the upper limit of the range of variable α, and LB1 is the vector composed of the lower limit of the range of variable K and the lower limit of the range of variable α; S23: Start the iterative operation. The initial number of iterations is 0, and the number of iterations is incremented by 1 for each iteration. If the number of iterations is less than or equal to 2 / 3 of the total number of iterations, then a number between [0, 1] is randomly generated. If this number is less than 0.5, update the position according to Method 1, otherwise update the position according to Method 2; if the number of iterations is greater than 2 / 3 of the total number of iterations, then a number between [0, 1] is randomly generated again. If this number is less than 0.5, update the position according to Method 3, otherwise update the position according to Method 4; Calculate the fitness value of the population after updating the position, compare the fitness of the current best individual and the best individual after iteration, and retain the optimal individual; S24: Determine whether the number of iterations is less than the total number of iterations. If so, return to step S23 to continue the iteration. If the number of iterations is equal to the total number of iterations, stop the iteration and output the optimal solutions of the decomposition mode number and the penalty factor; where the average value of the signal difference is used as the fitness, and the calculation formula is as follows: Among them, is the sum of the decomposed IMF components, f(n) is the original sequence, and N1 is the number of discrete points of the sequence; The methods for updating the position include: Method 1, vertical dive attack, which is mathematically expressed as: Among them, X1(t + 1) is the position of the t + 1-th iteration of the Tianying usage method 1, and X best (t) is the best position after the t-th iteration, T1 is the total number of iterations, and X M (t) is the average position of the population after the t-th iteration, and rand1 represents a number randomly taken in the interval [0, 1]; Method 2, contour flight and short glide attack, which is mathematically expressed as: X2(t + 1) = X best (t) × Levy + X R (t) + (y - x) × rand2 Among them, X2(t + 1) is the position of the t + 1-th iteration of the method 2 of the Tianying, and X R (t) is a random one among the N2 positions of the t-th iteration, N2 is the total number of Tianying in the population, rand2 represents a number randomly taken in the interval [0, 1], Levy is the Levy flight distribution function, and y and x are in a spiral shape, satisfying: where r is the spiral radius and θ is the spiral angle, and there is: r＝r1+0.00565×D1 where r1 is a value between [1, 20], and D1 is any integer from 1 to the length of the search space; Method 3, low-altitude flight and slow descent attack, which is mathematically expressed as: X3(t + 1) = (X best (t) - X M (t)) × 0.1 - rand3 + ((UB1 - LB1) × rand4 + LB1) × 0.1 where X3(t + 1) is the position of the eagle using Method 3 at the (t + 1)-th iteration, and rand3 and rand4 respectively represent a randomly taken number in the interval [0, 1]; Method 4, land walking attack, which is mathematically expressed as: X4(t + 1)= QF(t)×X best (t)- G1×(1 - X(t))×rand5 - G2×Levy where X4(t + 1) is the position of the eagle using Method 4 at the (t + 1)-th iteration, X(t) is the position at the t-th iteration, rand5 represents a randomly taken number in the interval [0, 1], and QF(t) is the quality function of the balanced search strategy, and there is: where rand6 represents a randomly taken number in the interval [0, 1]; G1 is various movements of the eagle during the prey's escape, and there is: G1 = 2 × rand7 - 1 where rand7 represents a randomly taken number in the interval [0, 1]; G2 is a decreasing value from 2 to 0, indicating the flight slope of the eagle following the prey from the first position to the last position during the prey's escape, and there is: where T1 is the total number of iterations.

4. A short-term load forecasting method for a power system according to claim 1, characterized in that, In step S2, the output intrinsic mode function components are divided according to the central frequency. Those with a central frequency less than or equal to the frequency threshold are low-frequency intrinsic mode function components, and those with a central frequency greater than the frequency threshold are high-frequency intrinsic mode function components.

Images