Park-level integrated energy system multi-element load prediction method and system

Abstract

The invention relates to a park-level integrated energy system multi-element load prediction method and system, and the method comprises the following steps: carrying out the intrinsic mode decomposition of an original electric, cold and thermal load sequence through employing ICEEMDAN, and obtaining a strong non-stationary component sequence; obtaining a stationary sequence component group and a non-stationary sequence component group of original electric, cold and heat loads by using VDM; performing data dimension reduction on the weather calendar rule feature data set by using LASSO; inputting the original power, cold and heat load non-stationary sequence components and the weather calendar rule principal component data set into a BiLSTM neural network to obtain predicted power, cold and heat load non-stationary sequence components; inputting the original power, cold and heat load stationary sequence components and the weather calendar rule principal component data set into an ARMA model to obtain predicted power, cold and heat load stationary sequence components; and adding the load stationary sequence component and the load non-stationary sequence component to obtain predicted power, cold and heat load values. Compared with the prior art, the prediction difficulty is reduced, and the prediction precision is improved.

Description

Technical Field

[0001] The present invention relates to the field of multi-element load forecasting for integrated energy systems, and in particular to a method and system for forecasting multi-element loads for integrated energy systems at a park level. Background Art

[0002] In order to solve the energy shortage and realize accurate energy management, the Integrated Energy System (IES) came into being. Multiple energy sources are supplied and converted in IES, which makes a certain coupling relationship between electricity, cooling and heating loads. With the development of the power market and the application of energy management systems, higher requirements are put forward for load forecasting.

[0003] Traditional load forecasting methods include statistical methods and machine learning methods, including Kalman filtering, multiple linear regression (MLR), support vector regression, random forest, etc. The statistical method model is simple and the prediction speed is fast, but it is only applicable to stationary sequences and it is difficult to accurately fit nonlinear sequences; machine learning has improved prediction performance when facing nonlinear sequences, but ignores the temporal relationship of samples.

[0004] In recent years, deep learning has been widely used in load forecasting. The most typical representative is the long short-term memory (LSTM) neural network and its variants. For example, patent application CN 114266374A discloses a multi-load forecasting method for an integrated energy system based on deep learning. The method aims at the diversified types of loads in the integrated energy system and the different load characteristics of each load, the insufficient information of the existing historical multi-load data mining, and the problem that the single prediction of multi-loads cannot take into account the characteristics of each load data. The method first analyzes the time series characteristics of the multi-load historical data, extracts the characteristic indicators based on the historical data, and uses them to deeply mine the load characteristics of the multi-load historical data in the integrated energy system. The input data set of the model is constructed based on the load characteristic indicators and the multi-load historical data. Then, based on the deep learning algorithm, the coupling information between the multi-loads is taken into account through the multi-task learning theory to realize the training of the multi-load short-term forecasting model, thereby completing the multi-load forecasting of the integrated energy system. In addition, deep belief networks and the like are also used. Although deep learning can mine the nonlinear relationship between sequences and learn their time series characteristics, the prediction accuracy needs to be improved when facing non-continuous sequences. Therefore, some scholars have proposed a combined prediction method to improve prediction accuracy by combining the respective advantages of different models.

[0005] However, most current research methods only predict power loads. Single load prediction methods ignore the coupling characteristics between multiple loads and cannot accurately predict IES multiple loads. Currently, there are few studies on IES multiple load prediction, and most of the current research is aimed at regional IES. For user-level IES, multiple loads are more random and volatile, and load prediction faces greater challenges. Summary of the invention

[0006] The purpose of the present invention is to overcome the defects of the above-mentioned prior art and to provide a multi-element load forecasting method and system for a park-level integrated energy system, thereby effectively improving the load forecasting accuracy.

[0007] The purpose of the present invention can be achieved by the following technical solutions:

[0008] A multi-element load forecasting method for a park-level integrated energy system comprises the following steps:

[0009] The original electric load sequence, the original cooling load sequence and the original heating load sequence are respectively subjected to intrinsic mode decomposition using the complete set empirical mode decomposition of improved adaptive noise to obtain a load decomposition sequence group, which is sorted according to the stationarity strength of each load decomposition sequence component to obtain a strong non-stationary component sequence through analysis;

[0010] The sequence component group obtained by applying variational mode decomposition to the strong non-stationary component sequence is subjected to hyperparameter optimization using the NRBO algorithm to obtain the stationary sequence component group and the non-stationary sequence component group of the original cooling load, the stationary sequence component group and the non-stationary sequence component group of the original heating load, and the stationary sequence component group and the non-stationary sequence component group of the original electric load;

[0011] LASSO is used to reduce the dimension of the weather calendar rule feature data set, remove redundant data, and obtain the weather calendar rule principal component data set;

[0012] The original electric load steady sequence component, the original cooling load steady sequence component, the original heating load steady sequence component and the weather calendar rule principal component data set are combined into an input data set, and a predicted electric load steady sequence component, a predicted cooling load steady sequence component and a predicted heating load steady sequence component are obtained through ARMA model prediction;

[0013] The original electric load non-stationary sequence component, the original cooling load non-stationary sequence component, the original heating load non-stationary sequence component and the weather calendar rule principal component data set are used to form a neural network input data set, and the predicted electric load non-stationary sequence component, the predicted cooling load non-stationary sequence component and the predicted heating load non-stationary sequence component are predicted by the BiLSTM neural network model;

[0014] All the predicted load steady sequence components predicted by the ARMA model and all the predicted load non-steady sequence components predicted by the BiLSTM neural network model are added together to obtain the final predicted electric load value, predicted cooling load value and predicted heating load value.

[0015] Furthermore, the steps of the complete set empirical mode decomposition of the improved adaptive noise are as follows:

[0016] Introduce the original load sequence;

[0017] A Gaussian white noise sequence with a mean of 0 and a normal distribution is added to the original load sequence, and residual components and IMF components of each order are obtained, and the final load decomposition sequence group is obtained through processing.

[0018] Furthermore, the residual component and the IMF component are as follows:

[0019] δ k = <M[δ k-1 +γ k-1 E k ω i (t)]>

[0020] f k =δ k-1 -δ k

[0021] In the formula, δ k is the k-th order residual component, f k is the k-th order IMF component, ω i (t) is the Gaussian white noise with a mean of 0 and a normal distribution added for the i-th time, γ k-1 is the adaptive coefficient for solving the k-1th order modal component, M(·) is the partial signal mean, E k (·) is the k-th modal component after VMD decomposition, and <·> is the overall mean.

[0022] Furthermore, the load decomposition sequence components in the final load decomposition sequence group are as follows:

[0023]

[0024] In the formula, load represents the sequence group of the original load sequence finally decomposed by ICEEMDAN, and F k represents the kth eigenmode component sequence in the final decomposition sequence group, F k (t) represents the load of the kth intrinsic mode component sequence at the tth moment in the final decomposition sequence group, U represents the number of sequence components after decomposition, and T represents the number of samples contained in a sequence component.

[0025] Furthermore, the strong non-stationary component sequence is subjected to hyperparameter optimization using the NRBO algorithm to obtain the following sequence group:

[0026]

[0027] In the formula, u r The component sequence group obtained by further decomposing the strongly non-stationary component sequence through VMD and then optimizing the hyperparameters using the NRBO algorithm is the rth eigenmode component sequence in the sequence component, u r (t) represents the load of the rth eigenmode component sequence at the tth moment in the sequence component after hyperparameter optimization, and R represents the number of sequence components obtained by decomposition.

[0028] Furthermore, the specific steps of using the NRBO algorithm for hyperparameter optimization are to use the NRM-based NRBO algorithm, define the search path by applying NRM, and use the vector set and operators NRSR and TAO to explore the search domain.

[0029] Furthermore, the NRBO algorithm searches for the best solution by generating an initial random population on the boundary of candidate solutions, each of which consists of decision vectors of different dimensions. The equation for generating a random population is as follows:

[0030]

[0031] in, It represents the load at the tth moment of the qth intrinsic mode component sequence in the non-stationary sequence components obtained after the load sequence is processed by ICEEMDAN and sorted by the stationarity strength. rand represents a random number between (0,1), and U1 represents the number of non-stationary sequence components.

[0032] Furthermore, the weather calendar rule feature data set is as follows:

[0033] data h =(month h ,day h ,week h ,hour h ,holiday h ,d h,1 ,d h,2 ,...,d h,Q )h∈[1,T]

[0034] In the formula, data h is the hth sample, h represents the hth hour, T is the number of samples, month h is the month serial number in the hth sample, day h is the day sequence number of the hth sample, weekh is the week number in the hth sample, hour h is the hour number in the hth sample, holiday h is the holiday serial number in the hth sample, For the qth * Types of weather factors, q * ∈[1,Q], Q is the number of weather factor types;

[0035] The principal component dataset of weather calendar rules after dimension reduction using the LASSO is as follows:

[0036]

[0037] In the formula, is the hth sample, h represents the hth hour, T is the number of samples, is the pth sample in the hth sample * elements, and H is the dimension of each sample.

[0038] Furthermore, the parameters of the ARMA model are estimated using a maximum likelihood estimation method or a least squares estimation method.

[0039] According to another aspect of the present invention, a multi-element load forecasting system for a park-level integrated energy system is provided, comprising:

[0040] An intrinsic mode decomposition module is used to perform intrinsic mode decomposition on the original electric load sequence, the original cooling load sequence and the original heating load sequence to obtain a load decomposition sequence group, sort the components of each load decomposition sequence according to the stationarity strength, and analyze and obtain a strong non-stationary component sequence;

[0041] A variational mode decomposition module is used to perform hyperparameter optimization on the sequence component group obtained by applying variational mode decomposition to the strong non-stationary component sequence using the NRBO algorithm to obtain a stationary sequence component group and a non-stationary sequence component group of the original cooling load, a stationary sequence component group and a non-stationary sequence component group of the original heating load, and a stationary sequence component group and a non-stationary sequence component group of the original electric load;

[0042] The weather data dimension reduction module is used to reduce the dimension of the weather calendar rule feature data set, remove redundant data, and obtain the weather calendar rule main component data set;

[0043] A load non-stationary sequence component prediction module is used to form a neural network input data set by combining the original electric load non-stationary sequence component, the original cooling load non-stationary sequence component, the original heating load non-stationary sequence component and the weather calendar rule principal component data set, and to obtain predicted electric load non-stationary sequence components, predicted cooling load non-stationary sequence components and predicted heating load non-stationary sequence components through a BiLSTM neural network model;

[0044] A load steady sequence component prediction module is used to form an input data set with the original electric load steady sequence component, the original cooling load steady sequence component, the original heating load steady sequence component and the weather calendar rule principal component data set, and obtain a predicted electric load steady sequence component, a predicted cooling load steady sequence component and a predicted heating load steady sequence component through ARMA model prediction;

[0045] The predicted value reconstruction module is used to add all the predicted load steady sequence components obtained by the ARMA model and all the predicted load non-steady sequence components obtained by the BiLSTM neural network model to obtain the final predicted electric load value, predicted cooling load value and predicted heating load value.

[0046] Compared with the prior art, the present invention has the following beneficial effects:

[0047] 1. The present invention uses ICEEMDAN to perform intrinsic mode decomposition on the original cooling load sequence, the original heating load sequence and the original electric load sequence respectively to obtain IMF components with different time scales, uses VMD to perform secondary decomposition on the strongly non-stationary variables obtained by decomposition, uses the NRBO algorithm to optimize the parameters of the load sequence obtained after the secondary decomposition, and then inputs it into the bidirectional long short-term memory network model, thereby improving the prediction accuracy.

[0048] 2. The present invention uses VMD to decompose multivariate loads into different intrinsic mode functions IMFs, and uses the compression estimation algorithm LASSO to reduce the dimension of the weather calendar rule feature set data, thereby reducing the difficulty of prediction.

[0049] 3. The present invention uses the BiLSTM neural network to deeply explore the coupling characteristics of multivariate loads, and predicts the decomposed non-stationary and stationary component sequences based on the BiLSTM neural network and the ARMA model respectively. Finally, the prediction results are reconstructed to obtain the final prediction results, which effectively improves the prediction accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1 A schematic diagram of the process of the multi-element load forecasting method for the park-level integrated energy system proposed by the present invention;

[0051] Figure 2This is a schematic diagram of LASSO path selection in the embodiment;

[0052] Figure 3 Schematic diagram of MAPE distribution of the prediction model proposed by the present invention in the embodiment. DETAILED DESCRIPTION

[0053] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. This embodiment is implemented based on the technical solution of the present invention, and provides a detailed implementation method and specific operation process, but the protection scope of the present invention is not limited to the following embodiments.

[0054] Abbreviations involved:

[0055] Integrated Energy System: Integrated Energy System, IES

[0056] Intrinsic Computing ExpressiveEmpirical Mode Decomposition With Adaptive Noise, ICEEMDAN

[0057] Variational Mode Decomposition: VMD

[0058] Newton-Raphson Method: Newton-Raphson Method, NRM

[0059] Newton-Raphson Optimizer: Newton-Raphson-Based Optimizer, NRBO

[0060] Compression estimation algorithm: Least Absolute Shrinkage and Selection Operator, LASSO

[0061] Auto-Regressive Moving Average, ARMR

[0062] Bidirectional Long Short-Term Memory Network: Bi-directional Long Short-Term Memory, BiLSTM

[0063] Example 1

[0064] This embodiment provides a multi-element load forecasting method for a park-level integrated energy system. Figure 1 As shown, the following steps are included:

[0065] This embodiment uses the cold, heat and electricity multi-load data of IES in Shandong Qingdao Sino-German Ecological Park as the experimental raw data. The input feature set includes two types: local weather and calendar rules. Meteorological data is collected from the official website of the National Renewable Energy Laboratory of the United States, including 8 categories of temperature, humidity, air pressure, wind speed, wind direction, cloud type, dew point and surface emissivity. At the same time, four types of calendar rules, namely time, day, month and holidays, are considered as input features, and the feature data set has a total of 12 dimensions.

[0066] S1. Use ICEEMDAN to perform intrinsic mode decomposition on the original electric, cold and hot load sequences respectively to obtain load decomposition sequence groups. Then sort them according to the stationarity strength of each load decomposition sequence component and analyze to obtain the strong non-stationary component sequence.

[0067] Data from January 1, 2020 to December 31, 2020 were selected. In order to fully cover the data for the whole year and fully consider the seasonal load changes, the data for each quarter were divided into training set and test set in a ratio of 8:2, and the training set was further divided by 5-fold cross-validation to construct a full-year data set. The time step was set to 1h.

[0068] In order to simply and quickly fill in the missing values ​​that may appear in the data, the median of the previous and next moments is used to complete the data. This processing can also reduce data noise to a certain extent. At the same time, considering that the units of electricity, cooling and heating loads are different, the three load data are normalized.

[0069] The original electrical load sequence is as follows:

[0070] x e (t),t∈[1,T]

[0071] In the formula, x e (t) is the original electric load at time t in the original electric load sequence, and T is the number of original electric loads;

[0072] The original cooling load sequence is as follows:

[0073] x c (t),t∈[1,T]

[0074] In the formula, x c (t) is the original cooling load at time t in the original cooling load sequence;

[0075] The original heat load sequence is as follows:

[0076] x h (t),t∈[1,T]

[0077] In the formula, x h (t) is the original heat load at time t in the original heat load sequence.

[0078] First, the ICEEMDAN modal decomposition is performed on the three loads of electricity, cooling and heating. The IMF1 generated after this is a strongly non-stationary variable, so VMD needs to be used to perform secondary modal decomposition on it. The IMF components generated by the VMD secondary modal decomposition are all stationary components, which can be predicted by the ARMA algorithm, and this process takes a short time. The IMF1 variables generated by the three load decompositions are decomposed into multiple IMF components again.

[0079] The steps of ICEEMDAN decomposition are as follows:

[0080] Introduce the original load sequence;

[0081] A Gaussian white noise sequence with a mean of 0 and normal distribution is added to the original load sequence, and the residual components and IMF components of each order are obtained, and the final load decomposition sequence group is obtained through processing.

[0082] The residual component and IMF component are as follows:

[0083] δ k = <M[δ k-1 +γ k-1 E k ω i (t)]>

[0084] f k =δ k-1 -δ k

[0085] In the formula, δ k is the k-th order residual component, f k is the k-th order IMF component, ω i (t) is the Gaussian white noise with a mean of 0 and a normal distribution added for the i-th time, γ k-1 is the adaptive coefficient for solving the k-1th order modal component, M(·) is the partial signal mean, E k (·) is the k-th modal component after VMD decomposition, and <·> is the overall mean.

[0086] The load decomposition sequence components in the final load decomposition sequence group are as follows:

[0087]

[0088] In the formula, load represents the sequence group of the original load sequence finally decomposed by ICEEMDAN, and F k represents the kth eigenmode component sequence in the final decomposition sequence group, F k(t) represents the load of the kth intrinsic mode component sequence at the tth moment in the final decomposition sequence group, U represents the number of sequence components after decomposition, and T represents the number of samples contained in a sequence component.

[0089] The stationarity strength of each sequence component in the final load decomposition sequence group is obtained by the approximate entropy method. According to the order of stationarity strength from left to right from weak to strong, each sequence component in the final load decomposition sequence group is sorted to obtain the sequence components after stationarity strength sorting.

[0090] The stationarity strength of each series component is as follows:

[0091] load stat =[ξ1,...,ξ k ,...,ξ u ],k∈[1,U]

[0092] In the formula, load stat is the stationarity strength set of each sequence component, ξ k It represents the stationarity strength of the kth sequence component, and U represents the number of sequence components after decomposition.

[0093] The sequence components after sorting by stationary strength are as follows:

[0094] load SE =[g1,...,g γ ,...,g U ],γ∈[1,U],g γ =[g γ (1),...,g γ (t),...,g γ (T)] T ,t∈[1,T]

[0095] In the formula, load SE represents the sequence component after the stationary strength is sorted, g γ represents the γth intrinsic mode component sequence in the sequence components after the stationary strength is sorted, g γ (t) represents the load of the γth intrinsic mode component sequence at the tth moment in the sequence components after the stationary strength is sorted, U represents the number of sequence components after decomposition, and T represents the number of samples contained in a sequence component.

[0096] The sequence component group load after stationarity strength sorting SE The corresponding stationarity strength of each series component is as follows:

[0097]

[0098] In the formula, is the stationary strength set of each sequence component in the sequence component group after the stationary strength is sorted, χ γ ,γ∈[1,U] represents the stationarity strength of the γth sequence component, and U represents the number of sequence components after decomposition;

[0099] Taking τ1 as the dividing line between stationary and non-stationary, and τ2 as the dividing line between non-stationary and strongly non-stationary, the stationary strength of each sequence is analyzed as follows:

[0100] For any γ∈[1,U], if χ γ <τ1, then the sequence component is classified as a stationary sequence component, if τ1≤χ γ ≤τ2, then the sequence component is classified as a non-stationary sequence component. γ >τ2, the sequence component is classified as a strongly non-stationary sequence component.

[0101] S2. The sequence component group obtained by applying VDM to the strongly non-stationary component sequence is optimized by using the NRBO algorithm to obtain the stationary sequence component group and non-stationary sequence component group of the original electricity, cooling and heating loads.

[0102] The strongly non-stationary component sequence is further decomposed by VMD to obtain a component sequence group, and then the NRBO algorithm is used to optimize the hyperparameters to obtain the following sequence group:

[0103]

[0104] In the formula, u r The component sequence group obtained by further decomposing the strongly non-stationary component sequence through VMD and then optimizing the hyperparameters using the NRBO algorithm is the rth eigenmode component sequence in the sequence component, u r (t) represents the load at the tth moment of the rth eigenmode component sequence in the sequence component after the component sequence group is further decomposed by VMD and then optimized by NRBO algorithm for hyperparameters. R represents the number of sequence components obtained by decomposition.

[0105] The specific steps of using the NRBO algorithm for hyperparameter optimization are to use the NRM-based NRBO algorithm, define the search path by applying NRM, and use vector sets and operators NRSR and TAO to explore the search domain.

[0106] The NRBO algorithm searches for the best solution by generating an initial random population on the boundary of candidate solutions. Each population consists of decision vectors of different dimensions. The equation for generating a random population is as follows:

[0107]

[0108] in, It represents the load at the tth moment of the qth intrinsic mode component sequence in the non-stationary sequence components obtained after the load sequence is processed by ICEEMDAN and sorted by the stationarity strength. rand represents a random number between (0,1), and U1 represents the number of non-stationary sequence components.

[0109] NRSR can enable the population to explore the feasible area more accurately and obtain a better position. TAO is introduced to improve the effectiveness and accuracy of the NRBO algorithm in dealing with real-world problems.

[0110] The stationary sequence group obtained by ICEEMDAN decomposition of the original load sequence and the strong stationary sequence group obtained by ICEEMDAN decomposition of the original load sequence are combined again by NRBO algorithm to obtain the strong stationary sequence group, and finally obtain the number of stationary sequence components. The stationary sequence group obtained by the secondary decomposition of the original load sequence is defined as Load s , the expression is as follows:

[0111] Load s =[A1,...,A α ,...,A I ],α∈[1,I],A α =[A α (1),...,A α (t),...,A α (T)] T ,t∈[1,T]

[0112] In the formula, A α A represents the αth intrinsic mode component sequence in the stable sequence component group obtained after the original load sequence is decomposed twice and sorted from left to right in order of stability strength from weak to strong; A α (t) represents the load at the tth moment of the αth intrinsic mode component sequence in the stable sequence component group obtained after the original load sequence is decomposed twice and sorted from left to right in order of stability strength from weak to strong; I represents the number of final stable sequence components after secondary decomposition optimization; T represents the number of samples contained in each sequence component.

[0113] The non-stationary sequence group obtained by decomposing the original load sequence by ICEEMDAN and the strongly stationary sequence group obtained by decomposing the original load sequence by ICEEMDAN are combined by decomposing and optimizing the non-stationary sequence group by NRBO algorithm, and finally J non-stationary sequence components are obtained. The non-stationary sequence group obtained by the secondary decomposition of the original load sequence is defined as Load ns , the expression is as follows:

[0114] Load ns=[B1,...,B β ,...,B J ],β∈[1,J],B β =[B β (1),...,B β (t),...,B β (T)] T ,t∈[1,T]

[0115] In the formula, B β It represents the βth intrinsic mode component sequence in the non-stationary sequence component group obtained after the original load sequence is decomposed twice and sorted from left to right in the order of stationarity strength from weak to strong; B β (t) represents the load at the tth moment of the βth intrinsic mode component sequence in the non-stationary sequence component group obtained after the original load sequence is decomposed twice and sorted from left to right in order of stationarity strength from weak to strong; J represents the number of final stationary sequence components after secondary decomposition optimization.

[0116] The original electric, cold and hot load sequences are respectively decomposed twice and sorted from left to right from weak to strong in terms of stationarity strength to obtain the stationary and non-stationary sequence component groups corresponding to the original electric, cold and hot load sequences.

[0117] Original electric load steady sequence component group as follows:.

[0118]

[0119] In the formula, represents the αth intrinsic mode sequence component in the stationary sequence component group; represents the load of the αth intrinsic mode sequence component at the tth moment in the stationary sequence component group; I represents the number of stationary sequence components; T represents the number of samples contained in each sequence component;

[0120] The original electric load non-stationary sequence component group as follows:

[0121]

[0122] In the formula, represents the βth intrinsic mode sequence component in the non-stationary sequence component group; represents the load of the βth intrinsic mode sequence component in the non-stationary sequence component group at the tth moment; J represents the number of non-stationary sequence components.

[0123] Original cooling load steady sequence component group as follows:

[0124]

[0125] In the formula, represents the αth intrinsic mode sequence component in the stationary sequence component group, represents the load of the αth intrinsic mode sequence component at the tth moment in the stationary sequence component group;

[0126] Original cooling load non-stationary sequence component group as follows:

[0127]

[0128] In the formula, represents the βth intrinsic mode sequence component in the non-stationary sequence component group, It represents the load of the βth intrinsic mode sequence component in the non-stationary sequence component group at the tth moment.

[0129] Original heat load steady sequence component group as follows:

[0130]

[0131] In the formula, represents the αth intrinsic mode sequence component in the stationary sequence component group, represents the load of the αth intrinsic mode sequence component at the tth moment in the stationary sequence component group;

[0132] Original heat load non-stationary series component group as follows:

[0133]

[0134] In the formula, represents the βth intrinsic mode sequence component in the non-stationary sequence component group, It represents the load of the βth intrinsic mode sequence component in the non-stationary sequence component group at the tth moment.

[0135] S3. Use LASSO to reduce the dimension of the weather calendar rule feature data set, remove redundant data, and obtain the weather calendar rule principal component data set.

[0136] The weather calendar rule feature dataset is as follows:

[0137] data h =(month h ,day h ,week h ,hour h ,holiday h ,d h,1 ,d h,2,...,d h,Q )h∈[1,T]

[0138] In the formula, data h is the hth sample, h represents the hth hour, T is the number of samples, month h is the month serial number in the hth sample, day h is the day sequence number of the hth sample, week h is the week number in the hth sample, hour h is the hour number in the hth sample, holiday h is the holiday serial number in the hth sample, For the qth * Types of weather factors, q * ∈[1,Q], Q is the number of weather factor types;

[0139] The weather calendar rule principal component dataset after dimensionality reduction is as follows:

[0140]

[0141] In the formula, is the hth sample, h represents the hth hour, T is the number of samples, is the pth sample in the hth sample * elements, and H is the dimension of each sample.

[0142] The factors that affect the load in this example include X1-air pressure, X2-wind speed, X3-holidays, X4-cloud type, X5-dew point, X6-humidity, X7-wind direction, X8-temperature, X9-surface emissivity, X1-air pressure, X2-wind speed, X3-holidays, X4-cloud type, X5-dew point, X6-humidity, X7-wind direction, X8-temperature, X9-surface e 10 —Month, X 11 —Time, X 12 On the first day, the data was normalized and the LASSO method was used to select variables. Figure 2 shown.

[0143] According to the AIC standard calculation, the AIC takes the minimum value when s = 0.386. After removing the factors with very small coefficients or 0, four main factors are obtained, namely X8-temperature, X6-humidity, X 10 —month, X3—holidays, as shown in Table 1.

[0144] Table 1 Results of principal component extraction by LASSO

[0145]

[0146]

[0147] S4. The original steady series components of electricity, cooling and heating loads and the principal component data set of weather calendar rules are combined into the input data set of the ARMA model to obtain the predicted steady series components of electricity, cooling and heating loads.

[0148] The original cooling load steady series component obtained by secondary decomposition, the original heating load steady series component obtained by secondary decomposition, the original electric load steady series component obtained by secondary decomposition and the weather calendar rule principal component data set are combined into an ARMA input data set, and prediction is performed using the ARMA algorithm model.

[0149] S5. The original steady sequence components of electricity, cooling and heating loads and the principal component data set of weather calendar rules are combined into the input data set of the BiLSTM neural network model to obtain the predicted non-steady sequence components of electricity, cooling and heating loads.

[0150] The non-stationary sequence components of the original electric load obtained by secondary decomposition, the non-stationary sequence components of the original cooling load obtained by secondary decomposition, the non-stationary sequence components of the original heat load obtained by secondary decomposition, and the principal component data set of weather calendar rules are used to form the neural network input data set, and the NRBO algorithm is used to optimize the hyperparameters of the data set obtained after secondary decomposition, and then it is input into the BiLSTM neural network model for prediction.

[0151] BiLSTM is composed of two independent LSTMs, which is a combination of LSTM neural network and bidirectional RNN, responsible for processing the input sequence from two directions (forward and reverse). This allows the model to obtain information before and after the current time step for training, effectively solving the long-term dependency problem, with stronger nonlinear expression ability, and can better handle uncertainty in data.

[0152] S6. Add all the load steady sequence components predicted by the ARMA model and all the load non-steady sequence components predicted by the BiLSTM neural network model to obtain the final predicted electric load value, predicted cooling load value and predicted heating load value.

[0153] In this embodiment, the root mean square error E RMSe , mean absolute percentage error E mAPE As an indicator for evaluating the results of the model predicting a single load, the expression is as follows:

[0154]

[0155] In the formula, n is the sample size of the test set, y i and They are respectively the predicted value and actual value of the cooling, heating and electricity loads at time i.

[0156] The comprehensive mean absolute percentage error I MAPEAs an indicator of the result of the comprehensive evaluation model predicting multivariate loads, the expression is as follows:

[0157]

[0158] In the formula, β m is the importance ratio of the mth type of load, is the MAPE of the mth type of load.

[0159] The BiLSTM network has two hidden layers, with 32 and 64 neurons respectively; the learning rate is 0.01; three output fully connected layers are set, with 1 layer each; the Adam algorithm is used as the training optimization algorithm; the time step parameter of LSTM is 24; and the number of iterations is 600.

[0160] In order to analyze the influence of different decomposition and reconstruction strategies on prediction accuracy, the secondary decomposition and reconstruction strategy proposed in the present invention is experimentally compared with no modal decomposition (direct prediction) and ICEEMDAN single modal decomposition. The error comparison is shown in Table 2 below.

[0161] Table 2 Error comparison under different decomposition and reconstruction strategies

[0162]

[0163] As can be seen from Table 2, the prediction accuracy of direct prediction without modal decomposition is the worst, with MAPE of electricity, cooling and heating loads being 2.46%, 4.64% and 4.03% respectively; the accuracy after primary modal decomposition is second best, with MAPE of the three loads being 1.48%, 2.36% and 2.46% respectively, which is nearly doubled compared with the direct prediction accuracy; the secondary modal decomposition accuracy using ICEEMDAN and VMD decomposition methods respectively is the highest among the three, with MAPE of electricity, cooling and heating loads being 0.46%, 1.03% and 0.93% respectively, which is nearly 4 times higher than the direct prediction accuracy, and the prediction accuracy of electricity load is more than 5 times higher, which is also doubled compared with the primary modal decomposition accuracy using only ICEEMDAN method.

[0164] In order to verify the superiority of the data processing method and prediction model proposed in this paper, the prediction model proposed in this invention (abbreviated as LVNA) is compared with several common combination models. The control models are: KPCA combined deep belief network (abbreviated as KPCA-DBN); LSTM-MLR model and MTL-BiLSTM model. The comparison of prediction result errors is shown in Table 3 below. The MAPE distribution of the model of the present invention is shown in Figure 3 shown.

[0165] Table 3 Error comparison of different prediction models

[0166]

[0167]

[0168] Comprehensive comparison shows that the prediction model proposed in this embodiment shows significant superiority in terms of fit with the actual load curve, especially at the mutation point of the load. Specifically, the MAPE errors of the model in predicting electricity, cooling and heating loads are 0.41%, 0.74% and 1.05%, respectively, and these error indicators are significantly lower than other reference models. In comparison, the prediction error of the KPCA-DBN model combining the deep belief network with KPCA is the most significant, and although the LSTM-MLR model and the improved BiLSTM model have improved prediction accuracy compared with the KPCA-DBN model, there is still a significant performance difference with the model proposed in this embodiment.

[0169] Further analysis shows that in the KPCA-DBN combination model, the deep belief network is composed of a bottom-level multi-layer restricted Boltzmann basis and a top-level back propagation neural network. The structure lacks a memory structure similar to the LSTM neural network, and its performance is slightly worse than that of the LSTM neural network. The LSTM-MLR combination model uses the LSTM neural network, which has greatly improved the prediction accuracy. However, since only historical data is trained in one direction, the global information of the historical data is not learned, and the correlation between the previous and next time is ignored, it is impossible to train an accurate prediction model. At the same time, due to the randomness and strong volatility of the park-level IES multivariate load, direct prediction using the LSTM-MLR combination model cannot obtain a high prediction accuracy.

[0170] In addition, to study the applicability of the proposed model for different prediction time scales, simulation experiments were conducted to observe and compare the performance of the proposed model under 15-minute and 1-hour prediction steps. For a more comprehensive evaluation, the simulation results of direct prediction using MTL-BiLSTM under 1-hour prediction step were also used for comparison. The error comparison is shown in Table 4 below.

[0171] Table 4 Prediction error under different prediction step sizes

[0172]

[0173] The comprehensive prediction results show that the prediction accuracy of the prediction model proposed in this embodiment decreases slightly when the prediction step is 1h. This may be because a larger prediction step increases the difficulty of the model in learning short-term fluctuations. Nevertheless, the simulation results show that although the prediction accuracy varies under different prediction step sizes, the proposed model is still effective at different time scales.

[0174] Example 2

[0175] This embodiment provides a multi-element load forecasting system for a park-level integrated energy system, including:

[0176] An intrinsic mode decomposition module is used to perform intrinsic mode decomposition on the original electric load sequence, the original cooling load sequence and the original heating load sequence to obtain a load decomposition sequence group, sort the components of each load decomposition sequence according to the stationarity strength, and analyze and obtain a strong non-stationary component sequence;

[0177] The variational mode decomposition module is used to optimize the hyperparameters of the sequence component groups obtained by applying variational mode decomposition to the strongly non-stationary component sequence using the NRBO algorithm, and obtain the stationary sequence component group and non-stationary sequence component group of the original cooling load, the stationary sequence component group and non-stationary sequence component group of the original heating load, and the stationary sequence component group and non-stationary sequence component group of the original electric load;

[0178] The weather data dimension reduction module is used to reduce the dimension of the weather calendar rule feature data set, remove redundant data, and obtain the weather calendar rule main component data set;

[0179] The load non-stationary sequence component prediction module is used to form a neural network input data set with the original electric load non-stationary sequence component, the original cooling load non-stationary sequence component, the original heating load non-stationary sequence component and the weather calendar rule principal component data set, and obtain the predicted electric load non-stationary sequence component, the predicted cooling load non-stationary sequence component and the predicted heating load non-stationary sequence component through the BiLSTM neural network model;

[0180] The load steady sequence component prediction module is used to form an input data set with the original electric load steady sequence component, the original cooling load steady sequence component, the original heating load steady sequence component and the weather calendar rule principal component data set, and obtain the predicted electric load steady sequence component, the predicted cooling load steady sequence component and the predicted heating load steady sequence component through the ARMA model prediction;

[0181] The predicted value reconstruction module is used to add all the predicted load steady sequence components obtained by the ARMA model and all the predicted load non-steady sequence components obtained by the BiLSTM neural network model to obtain the final predicted electric load value, predicted cooling load value and predicted heating load value.

[0182] The rest is the same as in Example 1.

[0183] The preferred specific embodiments of the present invention are described in detail above. It should be understood that a person skilled in the art can make many modifications and changes based on the concept of the present invention without creative work. Therefore, any technical solution that can be obtained by a person skilled in the art through logical analysis, reasoning or limited experiments based on the concept of the present invention on the basis of the prior art should be within the scope of protection determined by the claims.

Claims

1. A multi-element load forecasting method for a park-level integrated energy system, characterized in that: The following steps are involved: The original electric load sequence, the original cooling load sequence and the original heating load sequence are decomposed into intrinsic modes by using the complete set empirical mode decomposition of improved adaptive noise, and the load decomposition sequence group is obtained. The components of each load decomposition sequence are sorted according to the stationarity strength, and the strong non-stationary component sequence is obtained by analysis. The sequence component group obtained by applying variational mode decomposition to the strong non-stationary component sequence is subjected to hyperparameter optimization using the NRBO algorithm to obtain the stationary sequence component group and the non-stationary sequence component group of the original cooling load, the stationary sequence component group and the non-stationary sequence component group of the original heating load, and the stationary sequence component group and the non-stationary sequence component group of the original electric load; LASSO is used to reduce the dimension of the weather calendar rule feature data set, remove redundant data, and obtain the weather calendar rule principal component data set; The original electric load steady sequence component, the original cooling load steady sequence component, the original heating load steady sequence component and the weather calendar rule principal component data set are combined into an input data set, and a predicted electric load steady sequence component, a predicted cooling load steady sequence component and a predicted heating load steady sequence component are obtained through ARMA model prediction; The original electric load non-stationary sequence component, the original cooling load non-stationary sequence component, the original heating load non-stationary sequence component and the weather calendar rule principal component data set are used to form a neural network input data set, and the predicted electric load non-stationary sequence component, the predicted cooling load non-stationary sequence component and the predicted heating load non-stationary sequence component are predicted by the BiLSTM neural network model; All the predicted load steady sequence components predicted by the ARMA model and all the predicted load non-steady sequence components predicted by the BiLSTM neural network model are added together to obtain the final predicted electric load value, predicted cooling load value and predicted heating load value.

2. The multi-element load forecasting method for a park-level integrated energy system according to claim 1 is characterized in that: The steps of the complete ensemble empirical mode decomposition of the improved adaptive noise are as follows: Introduce the original load sequence; A Gaussian white noise sequence with a mean of 0 and a normal distribution is added to the original load sequence, and residual components and IMF components of each order are obtained, and the final load decomposition sequence group is obtained through processing.

3. The multi-element load forecasting method for a park-level integrated energy system according to claim 2 is characterized in that: The residual component and the IMF component are as follows: d k = <M[δ k-1 +g k-1 E k oh i (t)]) f k =d k-1 -d k In the formula, δ k is the k-th order residual component, f k is the k-th order IMF component, ω i (t) is the Gaussian white noise with a mean of 0 and a normal distribution added for the i-th time, γ k-1 is the adaptive coefficient for solving the k-1th order modal component, M(·) is the partial signal mean, E k (·) is the k-th modal component after VMD decomposition, and <·> is the overall mean.

4. The multi-element load forecasting method for a park-level integrated energy system according to claim 2 is characterized in that: The load decomposition sequence components in the final load decomposition sequence group are as follows: In the formula, load represents the sequence group of the original load sequence finally decomposed by ICEEMDAN, and F k represents the kth eigenmode component sequence in the final decomposition sequence group, F k (t) represents the load of the kth intrinsic mode component sequence at the tth moment in the final decomposition sequence group, U represents the number of sequence components after decomposition, and T represents the number of samples contained in a sequence component.

5. The multi-element load forecasting method for a park-level integrated energy system according to claim 1 is characterized in that: The sequence group obtained by performing hyperparameter optimization on the strong non-stationary component sequence using the NRBO algorithm is as follows: In the formula, u r The component sequence group obtained by further decomposing the strongly non-stationary component sequence through VMD and then optimizing the hyperparameters using the NRBO algorithm is the rth eigenmode component sequence in the sequence component, u r (t) represents the load of the rth eigenmode component sequence at the tth moment in the sequence components after hyperparameter optimization, and R represents the number of sequence components obtained by decomposition.

6. The multi-element load forecasting method for a park-level integrated energy system according to claim 1, characterized in that: The specific steps of using the NRBO algorithm for hyperparameter optimization are to use the NRM-based NRBO algorithm, define the search path by applying NRM, and use the vector set and operators NRSR and TAO to explore the search domain.

7. The multi-element load forecasting method for a park-level integrated energy system according to claim 6, characterized in that: The NRBO algorithm searches for the best solution by generating an initial random population on the boundary of candidate solutions. Each population consists of decision vectors of different dimensions. The equation for generating a random population is as follows: in, It represents the load at the tth moment of the qth intrinsic mode component sequence in the non-stationary sequence components obtained after the load sequence is processed by ICEEMDAN and sorted by the stationarity strength. rand represents a random number between (0,1), and U1 represents the number of non-stationary sequence components.

8. The multi-element load forecasting method for a park-level integrated energy system according to claim 1, characterized in that: The weather calendar rule feature data set is as follows: data h =(month h ,day h ,week h ,hour h holiday h ,d h,1 d h,2 ...,d h，Q )h∈[1,T] In the formula, data h is the hth sample, h represents the hth hour, T is the number of samples, month h is the month serial number in the hth sample, day h is the day sequence number of the hth sample, week h is the week number in the hth sample, hour h is the hour number in the hth sample, holiday h is the holiday serial number in the hth sample, For the qth * Types of weather factors, q * ∈[1,Q], Q is the number of weather factor types; The principal component dataset of weather calendar rules after dimension reduction using the LASSO is as follows: In the formula, is the hth sample, h represents the hth hour, T is the number of samples, is the pth sample in the hth sample * elements, and H is the dimension of each sample.

9. The multi-element load forecasting method for a park-level integrated energy system according to claim 1, characterized in that: The parameters of the ARMA model are estimated using the maximum likelihood estimation method or the least squares estimation method.

10. A multi-element load forecasting system for a park-level integrated energy system, characterized in that: include: An intrinsic mode decomposition module is used to perform intrinsic mode decomposition on the original electric load sequence, the original cooling load sequence and the original heating load sequence to obtain a load decomposition sequence group, sort the components of each load decomposition sequence according to the stationarity strength, and analyze and obtain a strong non-stationary component sequence; A variational mode decomposition module is used to perform hyperparameter optimization on the sequence component group obtained by applying variational mode decomposition to the strong non-stationary component sequence using the NRBO algorithm to obtain a stationary sequence component group and a non-stationary sequence component group of the original cooling load, a stationary sequence component group and a non-stationary sequence component group of the original heating load, and a stationary sequence component group and a non-stationary sequence component group of the original electric load; The weather data dimension reduction module is used to reduce the dimension of the weather calendar rule feature data set, remove redundant data, and obtain the weather calendar rule main component data set; A load non-stationary sequence component prediction module is used to form a neural network input data set by combining the original electric load non-stationary sequence component, the original cooling load non-stationary sequence component, the original heating load non-stationary sequence component and the weather calendar rule principal component data set, and to obtain predicted electric load non-stationary sequence components, predicted cooling load non-stationary sequence components and predicted heating load non-stationary sequence components through a BiLSTM neural network model; A load steady sequence component prediction module is used to form an input data set with the original electric load steady sequence component, the original cooling load steady sequence component, the original heating load steady sequence component and the weather calendar rule principal component data set, and obtain a predicted electric load steady sequence component, a predicted cooling load steady sequence component and a predicted heating load steady sequence component through ARMA model prediction; The predicted value reconstruction module is used to add all the predicted load steady sequence components obtained by the ARMA model and all the predicted load non-steady sequence components obtained by the BiLSTM neural network model to obtain the final predicted electric load value, predicted cooling load value and predicted heating load value.

Images