Power load forecasting method in multi-energy coupling mode

Abstract

The present invention relates to a load forecasting method based on a multi-energy coupling scene, which comprises the following steps: step 1, establishing a multilevel indicator system of key influencing factors of load demand in a multi-energy mode; step 2, obtaining key influencing factors influencing the total load demand and the total supply of the multi-energy coupling load; step 3, normalizing the data of the key influencing factors extracted in step 2, and initializing population characteristic parameters of adaptive firework algorithm (AFWA); step 4, forecasting the regional total power demand and regional coupling power supply respectively by adopting LSSVM optimized by the AFWA algorithm; and step 5, forecasting the net power load demand on the regional coupling power supply. The present invention improves the calculation efficiency and model stability and also ensures the forecasting accuracy.

Claims

1. A load forecasting method based on a multi-energy coupling scene, comprising the following steps: Step 1, establishing a multilevel indicator system of key influencing factors of load demand in a multi-energy coupling mode, and collecting corresponding data sets according to key influencing factors in the indicator system; Step 2, extracting key influencing factors most correlated with a total load demand and a total supply of multi-energy coupling load in the indicator system established in step 1 respectively by using a maximum and minimum redundancy method, removing surplus influencing factors with redundancy information, and finally obtaining the key influencing factors influencing the total load demand and the total supply of the multi-energy coupling load; Step 3, normalizing the data of the key influencing factors extracted in step 2, and initializing population characteristic parameters of an AFWA adaptive firework algorithm; Step 4, based on the key influencing factor data normalized in step 3 and initialized population characteristic parameters of the AFWA adaptive firework algorithm, forecasting regional total power demand and regional coupling power supply respectively by using LSSVM optimized by AFWA algorithm; and Step 5, based on forecast results of the regional total power demand and the regional coupling power supply in step 4, forecasting a net power load demand of the regional coupling power supply.

2. The load forecasting method based on the multi-energy coupling scene according to claim 1, comprising the following steps: the multilevel indicator system of the key influencing factors of the load demand in the multi-energy coupling mode in step 1 comprises factors influencing the regional total load demand and factors influencing the multi-energy coupling supply; wherein the factors influencing the regional total load demand comprises economic development, industry structure and demand coupling benefit; (1) economic development comprises: three-level influencing factor indicators of gross regional product and price index; (2) industry structure comprises: three-level influencing factor indicators of contribution of regional primary industry, contribution of regional secondary industry, contribution of regional tertiary industry; (3) demand coupling benefit comprises three-level influencing factor indicators of investment in energy-saving systems; wherein the factors influencing the multi-energy coupling supply comprise: renewable energy power supply, combined supply of cooling, heating and electricity, other energy power supply and energy supply economy; (1) the renewable energy power supply comprises: three-level influencing factor indicators of wind power generating capacity, solar power capacity, hydroelectric power generating capacity and other renewable energy power generating capacity; (2) the combined supply of cooling, heating and electricity comprises: three-level influencing factor indicators of combined supply power; (3) other energy power supply comprises: three-level influencing factor indicators of garbage power generating capacity; (4) energy supply economy comprises: three-level influencing factor indicators of unit price of main energy and investment in energy base stations.

3. The load forecasting method based on the multi-energy coupling scene according to claim 1, wherein the step 2 specifically comprises: (1) obtaining correlation between each influencing factor and load value by calculating mutual information of each influencing factor: $\begin{array}{cc}I\left(x_i,y\right)=\int \int p\left(x_i,y\right)\log \frac{p\left(x_i,y\right)}{p\left(x_i\right)p\left(y\right)}{\mathrm{dx}}_i\mathrm{dy}& \left(1\right)\end{array}$ wherein I(x.sub.i, y) is positive correlation between x.sub.i and y, x.sub.i is the i.sup.th influencing factor, y is the load value, p(x.sub.i, y), p(x.sub.i) and p(y) represent combined probability density of x.sub.i and y, probability density of x.sub.i and probability density of y; (2) by calculating the maximum correlation of each influencing factor, selecting a set of m key influencing factors with maximum correlation with the load value: $\begin{array}{cc}\max D\left(S,y\right),D=\frac{1}{.Math.S.Math.}\underset{x_i\in S}{.Math.}I\left(x_i,y\right),& \left(2\right)\end{array}$ wherein S is the set of the influencing factors {x.sub.i}, |S|=m is the number of key influencing factors to be selected, and D is a calculation result of the correlation between the selected m influencing factors and the load value; (3) by calculating minimum redundancy of each influencing factor, selecting the set of m key influencing factors with minimum information redundancy with the influencing factors: $\begin{array}{cc}\min R\left(S\right),R=\frac{1}{{.Math.S.Math.}^2}\underset{x_i,x_j\in S}{.Math.}I\left(x_i,x_j\right),& \left(3\right)\end{array}$ wherein R is a calculation result of mutual redundancy of the selected m influencing factors; (4) based on mRMR criterion, determining key influencing factors by adopting a gain search method, wherein assuming that m−1 characteristic influencing factors exist in S, the selection basis of the m characteristic influencing factor is shown in formula (5): $\begin{array}{cc}\underset{x_j\in X-S_{m-1}}{\max }\left[I\left(x_i,y\right)-\frac{1}{m-1}\underset{x_i\in S_{m-1}}{.Math.}I\left(x_j,x_i\right)\right]& \left(4\right)\end{array}$ the mRMR criterion is as follows:
max ϕ(D,R),ϕ=D−R  (5) (5) obtaining the score ϕ of each key influencing factor through formulas (1)-(5), and further selecting top m key influencing factors.

4. The load forecasting method based on the multi-energy coupling scene according to claim 1, wherein the step 3 specifically comprises: (1) normalizing the original data, and converting a preprocessed data set of the original data into dimensionless data from 0 to 1; and selecting part of the data as a training set for subsequent AFWA-LSSVM model training; the min-max standardization method is recommended as the data normalization method in the present invention, and the calculation formula is shown in formula (6): $\begin{array}{cc}{x}^{*}=\frac{x-x_{\min }}{x_{\max }-x_{\min }}& \left(6\right)\end{array}$ wherein x* is a data normalization result, x is an actual value of the data set, x.sub.min is a minimum value in the data set, and x.sub.max is a maximum value in the data set; (2) initializing population characteristic parameters of the AFWA adaptive firework algorithm, comprising a firework position, a population size, number of fireworks, total number of sparks, upper and lower boundaries of sparks, explosion range, data dimension, function evaluation times and iteration times.

5. The load forecasting method based on the multi-energy coupling scene according to claim 1, wherein the step 4 specifically comprises: (1) first, using the AFWA algorithm to optimize penalty parameters and kernel function width parameters of the LSSVW algorithm; (2) then inputting key influencing factor load data after the redundancy processing in step 2 to perform the AFWA-LSSVM model training; (3) after completing the process of optimizing the LSSVM model training by using AFWA, inputting the key influencing factor data set at a to-be-forecast time point, and calculating final load forecast results of total load demand and multi-energy coupling supply respectively.

6. The load forecasting method based on the multi-energy coupling scene according to claim 5, wherein the step 4 (1) specifically comprises: {circle around (1)} the calculation formula for minimizing the calculation error through firework algorithm optimization is as follows: $\begin{array}{cc}\min f\left(X\right)s.t.g_t\left(X\right)\le 0\left(t=1,2,\mathrm{.Math.},m\right)& \left(7\right)\end{array}$ wherein ƒ(X) is a target function and is set as that an error between the LSSVM forecast value and the actual value is minimum; g.sub.t(X) is a constraint function, m is number of the constraint functions, and X is an n-dimensional optimization variable; {circle around (2)} initializing the data:
x.sub.ij(0)=x.sub.ij.sup.L+rand(0,1)(x.sub.ij.sup.U−x.sub.ij.sup.L)  (8) wherein x.sub.ij(0) is spatial positions of initial firework individuals, i.e., initialized penalty parameters and kernel function width parameters of LSSVM; x.sub.ij.sup.U and x.sub.ij.sup.L are upper boundary and lower boundary of j.sup.th dimension respectively, i.e., parameter boundaries; rand (0,1) represents random numbers i=1, 2, . . . , N and j=1, 2, . . . ,n generated in an azimuth greater than 0 and less than 1; {circle around (3)} calculating an explosion operator, wherein the explosion operator mainly comprises explosion intensity, explosion range and displacement variation, wherein the explosion intensity is represented by the number of sparks, and the calculation way is as follows: $\begin{array}{cc}{N}_i=\hat{N}.Math.\frac{Y_{\max }-f\left(X_i\right)+\mathrm{.Math.}}{\underset{i=1}{\overset{N}{.Math.}}\left(Y_{\max }-f\left(X_i\right)\right)+\mathrm{.Math.}}& \left(9\right)\end{array}$ wherein N.sub.t is number of sparks of i.sup.th firework and is a set value; {circumflex over (N)} is a constant of controlling the total number of the sparks; Y.sub.max=max(ƒ(X.sub.t)) is an adaption value of the individual with the worst adaptability and individual with the greatest error of the load forecast result; ƒ(X.sub.t) is an adaption value of individual X.sub.t, i.e. load forecast result, and ε is an extremely small constant for preventing denominator from being zero; meanwhile, to limit the number of the sparks, the following rules are set: $\begin{array}{cc}{N}_i=\{\begin{array}{c}\mathrm{round}\left(N_{\min }\right),N_i<N_{\min }\\ \mathrm{round}\left(N_{\max }\right),N_i>N_{\max }\\ \mathrm{round}\left(N_i\right),\mathrm{Others}\end{array}& \left(10\right)\end{array}$ wherein round( ) is a rounding function; $\begin{array}{cc}{A}_i=\hat{A}.Math.\frac{f\left(X_i\right)-Y_{\min }+\mathrm{.Math.}}{\underset{i=1}{\overset{N}{.Math.}}\left(f\left(X_i\right)-Y_{\min }\right)+\mathrm{.Math.}}& \left(11\right)\end{array}$ wherein A.sub.t an explosion range of the i.sup.th firework, Â is the constant of limiting the maximum explosion range, and Y.sub.min=min(ƒ(X.sub.t)) is the adaption value of the firework individual with the highest forecast accuracy; an explosion radius is a distance between a specific individual and an optimum individual: $\begin{array}{cc}\hat{s}=\underset{s_i}{\arg \min }\left(d\left(s_i,s^{*}\right)\right)& \left(12\right)\\ f\left(s_i\right)>f\left(X\right)& \left(13\right)\end{array}$ wherein d indicates the calculation of some distance, i.e. the distance between various sparks, and each spark represents the width parameter in the LSSVM function and the kernel function; s.sub.i indicates all sparks, s* indicates the individuals with best adaptability in all sparks and fireworks, and X indicates the fireworks; formula (12) and formula (13) indicate that the explosion radius is the distance satisfying two conditions: condition 1: selecting the distance as the shortest distance from the alternative individual to the optimum individual; condition 2: alternative individuals are individuals having worse adaptability than the fireworks of this generation; and then carrying out displacement operation for the fireworks from different dimensions:
Δx.sub.ij=x.sub.ij+rand(0,A.sub.i)  (14) {circle around (4)} Variation operator:
x.sub.ij=x.sub.ijg  (15) wherein g˜N(1,1) is Gaussian distribution with both mean and variance of 1; {circle around (5)} mapping rule: Mapping the sparks beyond the boundary to a limit range by using a modular arithmetic rule:
x.sub.ij=x.sub.ij.sup.L+|x.sub.ij.sup.L|%(x.sub.i.sup.U−x.sub.ij.sup.L)  (16) wherein % represents modular arithmetic, x.sub.ij represents the position of i.sup.th individual on j.sub.th dimension, and x.sub.ij.sup.U, x.sub.ij.sup.L indicate the upper boundary and lower boundary on the j.sub.th dimension respectively, {circle around (6)} selection operation: calculating and selecting the distance between two individuals by using Euclidean distance: $\begin{array}{cc}R\left(X_i\right)=\underset{q=1}{\overset{K}{.Math.}}d\left(X_i,X_q\right)=\underset{q=1}{\overset{K}{.Math.}}.Math.X_i-X_q.Math.& \left(17\right)\end{array}$ wherein R(X.sub.i), indicates the sum of the distance from X.sub.i to all other individuals, d(X.sub.i,X.sub.q) indicates the Euclidean distance between any two individuals X.sub.i, X.sub.q, and K are position sets of all sparks after the Gaussian compilation, q∈K; $\begin{array}{cc}p\left(X_i\right)=\frac{R\left(X_i\right)}{\underset{q\in K}{.Math.}R\left(X_i\right)}& \left(18\right)\end{array}$ wherein p(X.sub.i) indicates the selection probability of X.sub.i; and the fireworks or sparks X.sub.i selected after reaching the maximum iteration times is the optimum penalty parameter and width parameter factor vector.

7. The load forecasting method based on the multi-energy coupling scene according to claim 5, wherein the step 4 (2) specifically comprises: {circle around (1)} after the nonlinear mapping, the obtained final decision function is as follows:
ƒ(x)=w.sup.T.Math.φ(x)+b  (19) wherein w∈R.sup.k(k>d) is a high-dimensional characteristic weight vector and reflects how each key influencing factor affects the load size and the influencing degree of each key influencing factor on the load value; x.sub.i∈R.sup.4 is the input of d-dimensional training samples and is specifically the sample of each key influencing factor; y.sub.i∈R is the output of training samples and the load value corresponding to the corresponding key influencing factor; and b∈R is an offset; {circle around (2)} according to the minimization principle of the structural risk, the forecast error is minimized to find the optimum high-dimensional characteristic weight vector w and the offset b; and the optimization problem is obtained as follows: $\begin{array}{cc}\min \frac{1}{2}{w}^{T}w+r\underset{i=1}{\overset{n}{.Math.}}{\xi }_i^2& \left(20\right)\end{array}$ wherein r>0 is penalty parameter, and ξ.sub.i is a slack variable; {circle around (3)} compared with a standard support vector machine algorithm, differences exist in constraint conditions:
y[w.sup.T.Math.φ(x.sub.i)+b]=1−ξ.sub.i,i=1,2, . . . ,n  (21) {circle around (4)} Solving the optimization problem by using Lagrangian function: $\begin{array}{cc}L=\frac{1}{2}{w}^T.Math.w+r.Math.\frac{1}{2}\underset{i=1}{\overset{n}{.Math.}}{\xi }_i^2-\underset{i=1}{\overset{n}{.Math.}}{\alpha }_i\left\{y_i\left[w^T.Math.\phi \left(x_i\right)+b\right]-1+{\xi }_i\right\}s.t.\{\begin{array}{c}w=\underset{i=1}{\overset{n}{.Math.}}{\alpha }_i{y}_i\phi \left(x_i\right)\\ \underset{i=1}{\overset{n}{.Math.}}{\alpha }_i{y}_i=0\\ {\alpha }_i=r{\xi }_i\\ y_i\left[w^T.Math.\phi \left(x_i\right)+b\right]-1+{\xi }_i=0\end{array}& \left(22\right)\end{array}$ wherein α.sub.i is a Lagrangian multiplicator vector, and α.sub.i>0, i=1, 2, . . . , n; {circle around (5)} Obtaining the final forecasting function as shown in formula (23): $\begin{array}{cc}f\left(x\right)=\underset{i=1}{\overset{n}{.Math.}}{\alpha }_{i}K\left(x,x_i\right)+b& \left(23\right)\end{array}$ wherein K(x.sub.i,x.sub.j)=φ(x.sub.i).sup.Tφ(x.sub.j) is a kernel function satisfying Mercer conditions: a radial basis kernel function is used as the kernel function, as shown in formula (24): $\begin{array}{cc}K\left(x,x_i\right)=\exp \left(-\frac{{.Math.x-x_i.Math.}^2}{2g^2}\right)& \left(24\right)\end{array}$ wherein g is a width factor of the kernel function.

8. The load forecasting method based on the multi-energy coupling scene according to claim 5, wherein the step 4 (3) specifically comprises: {circle around (1)} inputting the data sets of total load key influencing factors at a required forecasting time period into the AFWA-LSSVM model completing the total load data training, and then running the algorithm to obtain the forecasting result of the total load, {circle around (2)} inputting the data sets of multi-energy coupling key influencing factors at a required forecasting time period into the AFWA-LSSVM model completing the multi-energy coupling data training, and then running the algorithm to obtain the forecasting result of the multi-energy coupling.

9. The load forecasting method based on the multi-energy coupling scene according to claim 1, wherein the step 5 specifically comprises: the net power load in step 5 is a difference between the total demand load and the regional coupling supply load, and the calculation method of the energy demand based on this scene is shown in formula (25):
L.sub.t=L.sub.n−L.sub.c  (25) wherein L.sub.t is the load demand to be provided by the power grid, i.e. the net power load demand of the regional coupling power supply, L.sub.n is the total demand of all loads, and L.sub.c is energy supply of a non-power-grid composition in the multi-energy coupling mode.

Description

DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a load forecasting flow chart in a multi-energy coupling mode of the present invention.

(2) FIG. 2 shows analysis results of multi-energy coupling key influencing factors of the present invention.

(3) FIGS. 3(A)-3(C) are a calculation schematic diagram of net loads of the present invention, in which FIG. 3(A) is the calculation schematic diagram of net loads for the multi-energy coupling (forecasting set), FIG. 3(B) is the calculation schematic diagram of net loads for the multi-energy coupling (training set), and FIG. 3(C) is the calculation schematic diagram of net loads for the total demand (forecasting set); and FIG. 3(D) is the calculation schematic diagram of net loads for the total demand (training set); and FIG. 3(E) is the calculation schematic diagram showing of the calculation result of power grid load demand (forecasting set); and FIG. 3(F) is the calculation schematic diagram showing of the calculation result of power grid load demand (training set).

(4) FIG. 4 is a diagram of load forecasting results of a power grid in district A of the present invention.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

(5) The present invention will be further described below in detail with reference to the accompanying drawings.

(6) At present, the core problem of the load forecasting research is key influencing factors and model selection. The key influencing factors are mainly analyzed in three aspects, i.e. economy, energy and environment. The model selection is to select the model applicable to the current scene mainly from novel technical means such as cloud computation, big data, artificial intelligence, etc. It can be concluded that the load forecasting technology has already formulated a lot of specific solutions. However, among the existing load forecasting research achievements, there is almost no power load forecasting research for the multi-energy coupling circumstance. Therefore, it is required to analyze the scene characteristics and key influencing factors according to the current development situation of the multi-energy coupling and to establish a suitable power load forecasting model in the multi-energy coupling circumstance according to the analysis result.

(7) Since the multi-energy coupling scene has developed for a relatively short time, the corresponding medium-long-term power load forecasting technology is immature and has the characteristics of less samples, more influencing factors, complex correlation, etc. In the designing process of the forecasting model and method, it is required to establish a complete and accurate indicator system based on these characteristics in the characteristic selection process of the multi-energy coupling scene on one hand; and on the other hand, it is required to ensure the full use of the information and accuracy and effectiveness of the forecasting result in the establishing process of the forecasting model of the multi-energy coupling scene.

(8) A power load forecasting method based on a multi-energy coupling scene, as shown in FIG. 1, includes the following steps:

(9) Step 1, establishing a multilevel indicator system of key influencing factors of load demand in a multi-energy coupling mode, and collecting corresponding data sets according to the key influencing factors in the indicator system.

(10) The multilevel indicator system of the key influencing factors of the load demand in the multi-energy coupling mode in step 1 includes factors influencing the regional load total demand and factors influencing the multi-energy coupling supply.

(11) The factors influencing the regional load total demand includes economic development, industry structure and demand coupling benefit;

(12) (1) the economic development includes: three-level influencing factor indicators of gross regional product and price index;

(13) (2) the industry structure includes: three-level influencing factor indicators of contribution of regional primary industry, contribution of regional secondary industry, contribution of regional tertiary industry;

(14) (3) the demand coupling benefit includes three-level influencing factor indicators of investment in energy-saving systems.

(15) The factors influencing the multi-energy coupling supply comprise: renewable energy power supply, combined supply of cooling, heating and electricity, other energy power supply and energy supply economy;

(16) (1) the renewable energy power supply includes: three-level influencing factor indicators of wind power generating capacity, solar power capacity, hydroelectric power generating capacity and other renewable energy power generating capacity;

(17) (2) the combined supply of cooling, heating and electricity includes: three-level influencing factor indicators of combined supply power;

(18) (3) other energy power supply includes: three-level influencing factor indicators of garbage power generating capacity;

(19) (4) energy supply economy includes: three-level influencing factor indicators of unit price of main energy and investment in energy base stations.

(20) In the present embodiment, FIG. 2 shows an analysis method of multi-energy coupling key influencing factors provided by the present invention: factors influencing the regional load total demand mainly consider the development and change rule of the load and social and economic influencing factors as well as overall coupling benefit brought by the overall investment and construction of the regional multi-energy coupling. The analysis on the multi-energy coupling supply mainly considers the direct contradiction effect of the electric power generated by the multi-energy self-supply ways of the regional energy on the energy demand as well as relevant influencing factors of the energy economic market. Based on this, a multilevel indicator system of load demand key influencing factors in a multi-energy coupling mode is obtained, and corresponding data sets are collected according to the analysis result of the key influencing factors.

(21) By considering the difficulty in obtaining the data and the data requirement for actual calculation, the key influencing factors are set, as shown in Table 1.

(22) TABLE-US-00001 TABLE 1 Key influencing factors of load demand in multi-energy coupling mode B. Key content A. Key content decomposition C. Influencing factor A1. Total demand of B11. Economic C11. Gross regional energy development product C12. Price index B12. Industry C21. Contribution of structure regional primary industry C22. Contribution of regional secondary industry C23. Contribution of regional tertiary industry B13. Demand C31. Investment in coupling benefit energy-saving systems A2. Multi-energy B21. Renewable C41. Wind power coupling supply (load) energy power supply generating capacity C42. Solar power supply C43. Hydroelectric power generating capacity C44. Other renewable energy power generating capacity B22. Combined C51. Combined supply of cooling, heating supply power and electricity B23. Other energy C61. Garbage power power supply generating capacity B24. Energy supply C71. Unit price of economy main energy (coal) C72. Investment on energy base station

(23) Table 1 mainly includes the three-level forecasting indicator system, wherein key content indicators of primary items include: A1, energy demand and A2, multi-energy coupling supply; A1, the energy total demand includes secondary key content decomposition indicators: B11, economic development, B12, industry structure, and B13, demand coupling benefit; A2, the multi-energy coupling supply (load) includes secondary key content decomposition indicators: B21, renewable energy power supply, B22, combined supply of cooling, heating and electricity, B23, other energy power supply, and B24, energy supply economy; B11, the economic development includes tertiary influencing factor indicators: C11, gross regional product, and C12, price index; B12, the industry structure includes tertiary influencing factor indicators: C21, contribution of regional primary industry, C22, contribution of regional secondary industry, and C23, contribution of regional tertiary industry; B13, the demand coupling benefit includes tertiary influencing factor indicators: C31, investment in energy-saving systems, B21, the renewable energy power supply includes tertiary influencing factor indicators: C41, wind power generating capacity, C42, solar power supply, C43, hydroelectric power generating capacity, and C44, other renewable energy power generating capacity; B22, the combined supply of cooling, heating and electricity includes tertiary influencing factor indicators: C51, combined supply power; B23, other energy power supply includes tertiary influencing factor indicators: C61, garbage power generating capacity; and B24, the energy supply economy includes tertiary influencing factor indicators: C71, unit price of main energy (coal), and C72, investment in energy base stations.

(24) Data sets of the corresponding indicator system are collected according to the key indicator system shown in table 1 to prepare for the model training and parameter optimization.

(25) In the present invention, the power grid load changes in the southwestern region of China from the first quarter of 2004 to the first quarter of 2019 are taken as an example for empirical analysis.

(26) Step 2, extracting key influencing factors most correlated with a total load demand and multi-energy coupling load total supply in the indicator system established in step 1 respectively by using a maximum and minimum redundancy method, removing surplus influencing factors with redundancy information, and finally obtaining the key influencing factors influencing the total load demand and the multi-energy coupling load total supply.

(27) Further, the step 2 specifically includes:

(28) Mutual information (MI) is a theoretical method for evaluating the correlation of variables. The minimal redundancy maximal relevance (mRMR) algorithm is a characteristic selection method for maximizing the correlation between the characteristic variables and target variables and minimizing the redundancy information based on the mutual information theory.

(29) (1) Obtaining the correlation between each influencing factor and the load value by calculating mutual information of each influencing factor:

(30) $\begin{array}{cc}I\left(x_i,y\right)=\int \int p\left(x_i,y\right)\log \frac{p\left(x_i,y\right)}{p\left(x_i\right)p\left(y\right)}{\mathrm{dx}}_i\mathrm{dy}& \left(1\right)\end{array}$

(31) wherein I(x.sub.i,y) is positive correlation between x.sub.i and y, is the i.sup.th influencing factor, y is the load value, p(x.sub.i,y), p(x.sub.i) and p(y) represent a combined probability density of x.sub.i and y, probability density of x.sub.i and probability density of y respectively.

(32) The correlation between each influencing factor and the load value can be specifically calculated by calculating mutual information, and the greater the correlation is, the bigger the I(x.sub.i,y) is.

(33) (2) Selecting a set of m key influencing factors with the greatest correlation with the load value by calculating the maximum correlation of each influencing factor:

(34) $\begin{array}{cc}\max D\left(S,y\right),D=\frac{1}{.Math.S.Math.}\underset{x_i\in S}{.Math.}I\left(x_i,y\right),& \left(2\right)\end{array}$

(35) wherein S is the set of the influencing factors {x.sub.i}, |S|=m is the number of key influencing factors to be selected, and D is a calculation result of the correlation between the selected m influencing factors and the load value.

(36) Through this process, a set of m key influencing factors with the greatest correlation with the load value is selected.

(37) (3) By calculating the minimum redundancy of each influencing factor, selecting a set of m key influencing factors with minimum information redundancy with the influencing factors:

(38) $\begin{array}{cc}\min R\left(S\right),R=\frac{1}{{.Math.S.Math.}^2}\underset{x_i,x_j\in S}{.Math.}I\left(x_i,x_j\right),& \left(3\right)\end{array}$

(39) wherein R is a calculation result of mutual redundancy of the selected m influencing factors.

(40) Through this process, a set of m key influencing factors with minimum information redundancy is selected from the influencing factors:

(41) (4) Based on the mRMR criterion, determining key influencing factors by adopting a gain search method, wherein assuming that m−1 characteristic influencing factors exist in S, the selection basis of the m.sup.th characteristic influencing factor is shown in formula (5):

(42) 0$\begin{array}{cc}\underset{x_j\in X-S_{m-1}}{\max }\left[I\left(x_i,y\right)-\frac{1}{m-1}\underset{x_i\in S_{m-1}}{.Math.}I\left(x_j,x_i\right)\right]& \left(4\right)\end{array}$

(43) The mRMR criterion is as follows:
max ϕ(D,R),ϕ=D−R  (5)

(44) The mutual information difference (MID) is selected here, that is, the selected influencing factor set shall have greater correlation D with the load value, and the redundancy R among the factors shall be small.

(45) (5) Obtaining the score ϕ of each key influencing factor (the m.sup.th influencing factor) through formulas (1)-(5), and further selecting top m key influencing factors.

(46) In present embodiment, the characteristic influencing factors are selected according to the formulas (1)-(5), and the selection results are shown in Table 2 and table 3. It can be seen that the local multi-energy coupling self-supply load supply is mainly influenced by the wind power and the hydroelectric power. The load total demand is apparently affected by the price index, gross domestic product and total investment in fixed assets.

(47) TABLE-US-00002 TABLE 2 Multi-energy coupling mRMR calculation results Cumulative Cumulative Cumulative value of value of Cumulative value of hydroelectric other Cumulative value of wind power power power value of solar natural gas generating generating generating Investment power output capacity capacity capacity amount generating (hundred (hundred (hundred (hundred (ten Characteristic capacity (ten million cubic million million million thousand Fuel price factor thousand kWhs) meters) kWhs) kWhs) kWhs) Yuan) index mRMR 0 −0.16 1.63 0.94 0.52 0.62 0.60 value Order 6 7 1 2 5 3 4

(48) TABLE-US-00003 TABLE 3 mRMR calculation results of total load demand Gross Total domestic Contribution Contribution Contribution investment product in rate of rate of rate of in fixed current primary secondary tertiary assets season Data industry in industry in industry in (hundred (hundred of Characteristic current season current season current season million million price factor (%) (%) (%) Yuan) Yuan) index mRMR 0 −0.066 0.73 0.96 0.94 1.30 value Order 5 6 4 2 3 1

(49) The multi-energy coupling characteristic influencing factors of top three scores are selected and inputted into the AFWA-SVM model to perform the subsequent calculation.

(50) Step 3, normalizing the data of the key influencing factors extracted in step 2, and initializing population characteristic parameters of an AFWA adaptive firework algorithm.

(51) The key influencing factors extracted in step 2 cannot be directly inputted into the model for forecast. In this step, the model input and model output are identified. The data is normalized for the model calculation. The parameters of the forecasting model are also determined.

(52) Further, the step 3 specifically includes:

(53) (1) first, normalizing the original data, and converting the preprocessed data set of the original data into dimensionless data from 0 to 1; and selecting part of the data as a training set for subsequent model training;

(54) The min-max standardization method is recommended as the data normalization method in the present invention. The calculation formula is shown in formula (6):

(55) $\begin{array}{cc}{x}^{*}=\frac{x-x_{\min }}{x_{\max }-x_{\min }}& \left(6\right)\end{array}$

(56) wherein x* is a data normalization result, x is an actual value of the data set, x.sub.min is a minimum value in the data set (the data set of key influencing factors and the load data set), and x.sub.max is a maximum value in the data set.

(57) (2) Then, initializing population characteristic parameters of the AFWA adaptive firework algorithm, including a fireworks location, a population size, number of fireworks, total number of sparks, upper and lower bound of sparks, explosion range, data dimension, function evaluation times and iteration times;

(58) Some specific values of the firework position, population size, number of fireworks, total number of sparks, upper and lower boundaries of sparks and explosion range have good effect according to the research of experts and are generally set and may also be adjusted. The data dimension is set according to the characteristics of the optimization model and is set as 2 here. The function evaluation times and the iteration times may be determined according to the operation efficiency and accuracy.

(59) In present embodiment, the parameters of the AFWA model are initialized. The initialized parameters are shown in Table 4.

(60) TABLE-US-00004 TABLE 4 Parameter setting of adaptive firework algorithm Population size 5 Number of fireworks of Gaussian explosion 5 Total number of sparks 200 Upper boundary parameters of sparks 100 Lower boundary parameters of sparks 0.1 Control parameters of the explosion range 100 Data dimension 2 Function evaluation times 10000 Iteration times 1

(61) Step 4, based on the data of the key influencing factors normalized in step 3 and initialized population characteristic parameters of the AFWA (adaptive firework algorithm), forecasting regional total power demand and regional coupling power supply respectively by using LSSVM (Least Squares Support Vector Machine) optimized by AFWA (adaptive firework algorithm).

(62) In Step 4, model optimization and calculation are carried out for the initialized parameters in step 3.

(63) Further, the step 4 specifically includes:

(64) (1) First, using the AFWA algorithm to optimize penalty parameters and kernel function width parameters of the LSSVW algorithm, thereby improving the forecasting accuracy and calculation efficiency of LSSVM;

(65) In present embodiment, in a way of minimizing the calculation error through the optimization of the firework algorithm, the AFWA algorithm is used to optimize the penalty parameters and kernel function width parameters of the LSSVM algorithm. The AFWA algorithm mainly includes four calculation parts, i.e. explosion operator, variation operation, mapping rule and selection operation. According to the algorithm, the firework position is first formed by initializing the data, i.e. original penalty parameters and width parameters. Through the explosion algorithm, the screening range is enlarged and the function is converged. Through the variation operation, novel fireworks (penalty parameter and width parameter positions) are formed to ensure the complete search. Overflow is prevented by the mapping rule. The optimum fireworks (penalty parameter and width parameter positions) are continuously selected by the selection operation to reduce the calculation error of LSSVM. The specific calculation way of the optimization problem is described as follows:

(66) Further, the step 4 (1) specifically includes:

(67) {circle around (1)} The calculation formula for minimizing the calculation error through the optimization of firework algorithm is as follows:

(68) $\begin{array}{cc}\min f\left(X\right)s.t.g_t\left(X\right)\le 0\left(t=1,2,\mathrm{.Math.},m\right)& \left(7\right)\end{array}$

(69) wherein ƒ(X) is a target function and is set as that the error between the LSSVM forecast value and the actual value is minimum (the LSSVM model is already inserted into the AFWA adaptive function). g.sub.t(X) is a constraint function (since the AFWA is used to optimize the two parameters of LSSVM, here is set as the parameter range), m is the number of the constraint functions, and X is n-dimensional optimization variable.

(70) The process of the AFWA algorithm is interpreted based on this.

(71) {circle around (2)} Initializing the Data:
x.sub.ij(0)=x.sub.ij.sup.L+rand(0,1)(x.sub.ij.sup.U−x.sub.ij.sup.L)  (8)

(72) wherein x.sub.ij(0) is a spatial position of initial firework individual, i.e. penalty parameter and kernel function width parameter of initialized LSSVM; x.sub.ij.sup.U and x.sub.ij.sup.L are upper boundary and lower boundary of j.sub.th dimension respectively, i.e. parameter boundary; rand (0,1) represents a random number i=1, 2, . . . , N and j=1, 2, . . . , n generated in an azimuth greater than 0 and less than 1;

(73) {circle around (3)} Calculating an explosion operator, wherein the explosion operator mainly includes explosion intensity, explosion range and displacement variation, the explosion intensity is represented by the number of sparks, and the calculation way is as follows:

(74) $\begin{array}{cc}{N}_i=\hat{N}.Math.\frac{Y_{\max }-f\left(X_i\right)+\mathrm{.Math.}}{\underset{i=1}{\overset{N}{.Math.}}\left(Y_{\max }-f\left(X_i\right)\right)+\mathrm{.Math.}}& \left(9\right)\end{array}$

(75) wherein N.sub.i is the number of sparks of i.sup.th firework and is a set value; {circumflex over (N)} is a constant for controlling the total number of the sparks; Y.sub.max=max(ƒ(X.sub.i)) is an adaptive value of the individual with the worst adaptability and the individual with the maximum error of the load forecasting result; ƒ(X.sub.i) is adaption value of the individual X.sub.i, i.e. load forecasting result, and ε is an extremely small constant for preventing denominator from being zero;

(76) Meanwhile, to limit the number of the sparks, the following rules are set:

(77) $\begin{array}{cc}{N}_i=\{\begin{array}{c}\mathrm{round}\left(N_{\min }\right),N_i<N_{\min }\\ \mathrm{round}\left(N_{\max }\right),N_i>N_{\max }\\ \mathrm{round}\left(N_i\right),\mathrm{Others}\end{array}& \left(10\right)\end{array}$

(78) wherein round( ) is a rounding function;

(79) $\begin{array}{cc}{A}_i=\hat{A}.Math.\frac{f\left(X_i\right)-Y_{\min }+\mathrm{.Math.}}{\underset{i=1}{\overset{N}{.Math.}}\left(f\left(X_i\right)-Y_{\min }\right)+\mathrm{.Math.}}& \left(11\right)\end{array}$

(80) wherein A.sub.i is an explosion range of the i.sup.th firework, Â is the constant of limiting the maximum explosion range, and Y.sub.min=min (ƒ(X.sub.i)) is the adaption value of the firework (width parameter and penalty parameter) individual with the highest forecasting accuracy;

(81) For the reasonability of the calculation method for the explosion radius of the traditional firework algorithm, the adaptive firework algorithm uses the generated sparks to calculate the optimum firework explosion radius so as to realize the adaptive adjustment of the explosion radius.

(82) An explosion radius is a distance between a specific individual and an optimum individual:

(83) $\begin{array}{cc}\hat{s}=\underset{s_i}{\arg \min }\left(d\left(s_i,s^{*}\right)\right)& \left(12\right)\\ f\left(s_i\right)>f\left(X\right)& \left(13\right)\end{array}$

(84) wherein d indicates the calculation of a distance, i.e. the distance among various sparks, and each spark represents the width parameter in the LSSVM function and the kernel function; s.sub.i indicates all sparks, s* indicates the individual with the best adaptability in all sparks and fireworks, and X indicates the fireworks;

(85) The formula (12) and formula (13) indicate that the explosion radius is the distance satisfying two conditions: a: the selected distance is the shortest distance from the alternative individual to the optimum individual; b: alternative individuals are individuals having worse adaptability than the fireworks of this generation;

(86) Then the displacement operation is carried out for the fireworks from different dimensions:
Δx.sub.ij=x.sub.ij+rand(0,A.sub.i)  (14)

(87) {circle around (4)} Variation Operator:

(88) Gaussian variation is mainly used in the present invention.
x.sub.ij=x.sub.ijg  (15)

(89) wherein g˜N(1,1) is Gaussian distribution with both mean and variance of 1.

(90) {circle around (5)} Mapping rule: the mapping rule is an algorithm that maps the sparks beyond the boundary to a limited range through a method, mainly including modular arithmetic rules, mirror reflection rules, random mapping rules, etc. Modular arithmetic rules are used here.
x.sub.ij=x.sub.ij.sup.L+|x.sub.ij.sup.L|%(x.sub.ij.sup.U−x.sub.ij.sup.L)  (16)

(91) wherein % represents modular arithmetic, x.sub.ij represents the position of i.sup.th individual on j.sub.th dimension, and x.sub.ij.sup.U, x.sub.ij.sup.L indicate the upper boundary and lower boundary on the j.sub.th dimension respectively.

(92) {circle around (6)} Selection operation: the selection operation adopts the distance-based selection and random selection. A distance between two individuals is calculated and selected by using Euclidean distance.

(93) $\begin{array}{cc}R\left(X_i\right)=\underset{q=1}{\overset{K}{.Math.}}d\left(X_i,X_q\right)=\underset{q=1}{\overset{K}{.Math.}}.Math.X_i-X_q.Math.& \left(17\right)\end{array}$

(94) wherein R(X.sub.i) indicates the sum of the distance from X.sub.i to all other individuals, d(X.sub.i,X.sub.q) indicates the Euclidean distance between any two individuals X.sub.iX.sub.q, and K is the position set of all sparks after the Gaussian compilation, q∈K;

(95) $\begin{array}{cc}p\left(X_i\right)=\frac{R\left(X_i\right)}{\underset{q\in K}{.Math.}R\left(X_i\right)}& \left(18\right)\end{array}$

(96) wherein p(x.sub.i) indicates the selected probability of X.sub.i; and the fireworks or sparks X.sub.i selected after reaching the maximum iteration times are the optimum penalty parameter and width parameter factor vectors.

(97) (2) Then, inputting the processed key influencing factor load data to perform the model training;

(98) The main model in the process of inputting the processed key influencing factor load data to perform the model training is LSSVM model. The optimum penalty parameter and width parameter value are obtained through the AFWA model in step 4(1), and the model forms a fitted forecasting function through the machine learning (continuous training) process, thereby realizing the forecasting process. Aiming at the minimum functional forecasting error in the training process, the optimum high-dimensional characteristic weight vector w and offset b are found, and the forecasting fitting function is finally formed to complete the training process.

(99) Further, the step 4 (2) specifically includes:

(100) {circle around (1)} After the nonlinear mapping ψ(x)=(φ(x.sub.i), φ(x.sub.2), . . . , φ(x.sub.n)), the final decision function is obtained as follows:
ƒ(x)=w.sup.T.Math.φ(x)+b  (19)

(101) wherein w∈R.sup.k(k>d) is a high-dimensional characteristic weight vector, i.e. the influencing degree of each key influencing factor on the load value; x.sub.i∈R.sup.d is the input of d-dimensional training samples and is specifically samples of each key influencing factor; y.sub.i∈R is the output of training samples and the load value corresponding to the corresponding key influencing factor; and b∈R is an offset.

(102) {circle around (2)} According to the minimization principle of the structure risk, the forecasting error is minimized to find the optimum high-dimensional characteristic weight vector w and the offset b. The optimization problem is obtained as follows:

(103) $\begin{array}{cc}\min \frac{1}{2}{w}^{T}w+r\underset{i=1}{\overset{n}{.Math.}}{\xi }_i^2& \left(20\right)\end{array}$

(104) wherein r>0 is the penalty parameter, and ξ.sub.i is a slack variable.

(105) {circle around (3)} Compared with a standard support vector machine algorithm, differences exist in constraint conditions:
y[w.sup.T.Math.φ(x)+b]=1−ξ.sub.i,i=1,2, . . . ,n  (21)

(106) {circle around (4)} Solving the optimization problem by using Lagrangian function:

(107) 0$\begin{array}{cc}L=\frac{1}{2}{w}^T.Math.w+r.Math.\frac{1}{2}\underset{i=1}{\overset{n}{.Math.}}{\xi }_i^2-\underset{i=1}{\overset{n}{.Math.}}{\alpha }_i\left\{y_i\left[w^T.Math.\phi \left(x_i\right)+b\right]-1+{\xi }_i\right\}s.t.\{\begin{array}{c}w=\underset{i=1}{\overset{n}{.Math.}}{\alpha }_i{y}_i\phi \left(x_i\right)\\ \underset{i=1}{\overset{n}{.Math.}}{\alpha }_i{y}_i=0\\ {\alpha }_i=r{\xi }_i\\ y_i\left[w^T.Math.\phi \left(x_i\right)+b\right]-1+{\xi }_i=0\end{array}& \left(22\right)\end{array}$

(108) wherein α.sub.i is a Lagrangian multiplicator vector, and α.sub.i>0, i=1, 2, . . . , n.

(109) {circle around (5)} Obtaining the final forecasting function as shown in formula (23):

(110) $\begin{array}{cc}f\left(x\right)=\underset{i=1}{\overset{n}{.Math.}}{\alpha }_{i}K\left(x,x_i\right)+b& \left(23\right)\end{array}$

(111) wherein κ(x.sub.i,x.sub.j)=φ(x.sub.i).sup.Tφ(x.sub.j) is a kernel function satisfying Mercer conditions:

(112) A radial basis kernel function is used as the kernel function, as shown in formula (24):

(113) $\begin{array}{cc}K\left(x,x_i\right)=\exp \left(-\frac{{.Math.x-x_i.Math.}^2}{2g^2}\right)& \left(24\right)\end{array}$

(114) wherein g is a width factor of the kernel function.

(115) (3) After completing the process of optimizing the LSSVM model training by using AFWA, inputting the data set of key influencing factors at a required forecasting time point, and calculating the final load forecasting results of total load demand and multi-energy coupling supply respectively.

(116) After the AFWA optimization and LSSVM model training process is completed, the required forecasting model can be obtained, and then as long as there is data input, the output of required forecasting results can be obtained. The data sets of key influencing factors of total load demand and multi-energy coupling supply at the required forecasting time point are inputted respectively, and the final load forecasting results of the total load demand and multi-energy coupling supply can be obtained by model calculation.

(117) In present embodiment, the penalty parameters and kernel function width parameters after the LSSVM is optimized by AFWA are shown in Table 5:

(118) TABLE-US-00005 TABLE 5 Parameter setting of optimized Lssvm Penalty parameter c Width parameter g Multi-energy coupling 30.9715 60.4294 Total load demand 26.4350 61.0462

(119) The optimized parameters are used to train the multi-energy coupling demand and total demand LSSVM model. After the training, the selected key influencing factor set is inputted for forecast. Through this step, the regional total demand forecasting result and the regional multi-energy coupling power supply forecasting result can be obtained. The results are further calculated to obtain the energy demand forecasting result in the multi-energy coupling mode.

(120) The final forecasting result data is shown in Table 6. In order to more intuitively reflect the effect of the present invention, the results are shown in FIGS. 3(A)-3(C). The graph of the corresponding training set in FIGS. 3(A)-3(C) shows the fitting degree of the existing data forecasting result. The multi-energy coupling (forecasting set) and total demand (forecasting set) in FIG. 3(A) and FIG. 3(C) are result sets of the actual forecast in present step, and the difference between the multi-energy coupling and the total demand is a calculation result (forecasting set) of the power grid load demand (subsequent calculation).

(121) TABLE-US-00006 TABLE 6 Multi-energy coupling forecasting result and total demand forecasting result Actual value Forecasting Actual value Forecasting (coupling self value (coupling (total value (total No. supply) self supply) demand) demand) 1 43.53 44.01408 1144.58 1157.58 2 23.02 23.62906 583.7 618.0737 3 20.87 21.49441 606.83 610.3525 4 68.93 69.64405 1458.72 1547.29 5 81.54 79.7132 1376.27 1357.008 6 42.79 44.34388 1182.28 1199.143 7 34.81 33.62083 885.21 765.9002 8 33.71 33.82143 864.45 927.1944 9 59.35 59.25783 1310.27 1312.341 10 74.12 76.71408 1501.11 1638.691 11 100.27 99.53899 1613.74 1533.311 12 52.79 53.38457 1235.36 1272.739 13 69.93 73.1467 1404.35 1469.58

(122) Step 5, forecasting a net power load demand of the regional coupling power supply;

(123) The net power load in step 5 is a difference between the total demand load and the regional coupling supply load.

(124) The load demand forecast of the present invention refers to the calculation of the load demand of the regional power user on the power grid under the multi-energy coupling situation.

(125) The net power load demand forecasting result of the regional coupling power supply obtained in step 4 is not the forecasting result of the regional demand on the power grid. The regional load demand on the power grid needs to be further calculated. The self supply of part of energy is realized in the multi-energy coupling mode, so that the regional load demand on the power grid is relatively reduced compared with the load demand without the multi-energy coupling mode. Therefore, the regional load demand should be the difference between the total power demand and the multi-energy coupling self supply. The calculation method of the energy demand based on this scene is shown in formula (25).
L.sub.t=L.sub.n−L.sub.c  (25)

(126) wherein L.sub.t is the load demand required to provide by the power grid, i.e. a deserved result of the forecasting model of the present invention, L.sub.n is the total demand of all loads, and L.sub.c is energy supply of a non-power-grid composition in the multi-energy coupling mode.

(127) In present embodiment, FIG. 4 is a structural schematic diagram of the load demand in the multi-energy coupling mode. The regional energy supply is satisfied by the power grid and the multi-energy coupling self supply of integrated energy. In order to ensure the power supply, an accurate forecasting result of regional load demand on the power grid is needed. To solve this problem, the difference between the total load forecasting result and the multi-energy coupling supply forecasting result is used as the load forecasting result of the power grid.

(128) The total demand forecasting result and the multi-energy coupling forecasting result are calculated according to the formula (25) to obtain the forecasting result of the regional power demand on the power grid in the multi-energy coupling circumstance. In present embodiment, the calculation result is shown in Table 7 and intuitively shown as the forecasting result (forecasting set) of the load demand on the power grid in FIGS. 3(A)-3(C).

(129) TABLE-US-00007 TABLE 7 Forecasting result of regional power demand on power grid Actual Forecasting value value (demand (demand Actual Forecasting Actual of load of load value value value Forecasting users on users on Sequence (coupling (coupling (total value (total the power the power number self-supply) self-supply) demand) demand) grid) grid) 1 43.53 44.01408 1144.58 1157.58 1101.05 1113.566 2 23.02 23.62906 583.7 618.0737 560.68 594.4447 3 20.87 21.49441 606.83 610.3525 585.96 588.8581 4 68.93 69.64405 1458.72 1547.29 1389.79 1477.646 5 81.54 79.7132 1376.27 1357.008 1294.73 1277.295 6 42.79 44.34388 1182.28 1199.143 1139.49 1154.799 7 34.81 33.62083 885.21 765.9002 850.4 732.2794 8 33.71 33.82143 864.45 927.1944 830.74 893.373 9 59.35 59.25783 1310.27 1312.341 1250.92 1253.083 10 74.12 76.71408 1501.11 1638.691 1426.99 1561.977 11 100.27 99.53899 1613.74 1533.311 1513.47 1433.772 12 52.79 53.38457 1235.36 1272.739 1182.57 1219.355 13 69.93 73.1467 1404.35 1469.58 1334.42 1396.434

(130) To verify the advantages of the multi-energy forecasting model proposed in the present invention, results shown in Table 8 can be obtained by comparing the calculation efficiency of other load forecasting models.

(131) TABLE-US-00008 TABLE 8 Model comparison evaluation t MSE RMSE MAE MAPE % R{circumflex over ( )}2 Coupling mRMR-AFWA-LSSVM 21.31 2.05 1.431 1.10 2.08 0.9961 AFWA-LSSVM 43.97 0.66 0.81 0.71 2.01 0.9976 LSSVM 10.55 31.99 5.65 4.30 7.46 0.9395 AFWA-SVM 27.60 12.6 3.55 3.41 7.07 0.9761 BP 11.20 105.32 10.26 29.20 8.12 0.9237 Total mRMR-AFWA-LSSVM 26.99 4545.26 67.41 52.33 4.56 0.9561 demand AFWA-LSSVM 51.33 4455.12 66.74 51.79 4.51 0.9570 LSSVM 8.94 5877.98 76.66 70.22 6.75 0.9432 AFWA-SVM 33.59 7298.54 85.43 71.40 6.13 0.9295 BP 12.30 130532.40 361.29 1033.31 7.01 0.9031

(132) It can be discovered by analyzing the calculation results that the forecasting accuracy for the five indicators, i.e. MSE, RMSE, MAE, MAPE and R{circumflex over ( )}2 are as follows:

(133) AFWA−LSSVM>mRMR−AFWA−LSSVM>>LSSVM−AFWA−SVM>BP. In the coupling forecast, the MAPE difference is 0.07%, 5.39%, −0.90% and 1.03% respectively, and in the total demand forecast, the MAPE difference is 0.05%, 2.19%, 0.58% and 0.88% respectively. However, for the total calculation time, LSSVM<mRMR−AFWA−LSSVM<BP<AFWA−LSSVM<AFWA−SVM, and the total calculation time is 19.49, 23.5, 48.3, 61. 19 and 95.3 respectively.

(134) After the mRMR selection, the AFWA-LSSVM forecasting accuracy decreases slightly, while the forecasting time is reduced. Since the indicators selected by examples of the present invention are indicators screened by experience, the accuracy of mRMR decreases slightly in this process. However, the mRMR can greatly reduce the data processing workload in the actual load forecasting process with a large indicator system and rich data and may not have great impact on the forecasting accuracy. Therefore, mRMR-AFWA-LSSVM is generally an optimum forecasting model of the actual load.

(135) mRMR-AFWA-LSSVM has good forecasting accuracy and relatively high calculation performance. The present invention proposes the calculation method for forecasting the regional load demand on the power grid in the multi-energy coupling mode for the first time, which has significant reference value for power planning and economic and reliable operation of the power grid.

(136) It shall be emphasized that the embodiments of the present invention are illustrative rather than restrictive. Therefore, the present invention includes, but is not limited to the embodiments described in the detailed description. Other implementation modes obtained by those skilled in the art according to the technical solution of the present invention also belong to the protection scope of the present invention.