A load forecasting method integrating multi-source feature enhancement and chaotic optimization

Abstract

This invention discloses a load forecasting method integrating multi-source feature enhancement and chaotic optimization, comprising the following steps: preprocessing the data; selecting a core feature set from the preprocessed data; performing wavelet decomposition on the time series in the core feature set to obtain multiple frequency components; generating abrupt change features and integrating them into the corresponding frequency components, reconstructing them to form an enhanced feature set; calculating the comprehensive similarity between the target prediction date and historical dates, selecting a similar date sample set based on the comprehensive similarity, and loading the corresponding historical optimal model parameters to initialize the current prediction model; constructing an integrated prediction model; optimizing the integration weights of each branch model; and performing inverse normalization on the fusion result to obtain the final load forecast value. The advantages of this invention are: it achieves multi-dimensional breakthroughs in feature quality, abrupt change capture, optimization efficiency, and scenario adaptation, providing a high-precision and high-stability technical solution for power loads under complex operating conditions.

Description

Technical Field

[0001] This invention relates to the technical field of power system load forecasting, and in particular to a load forecasting method that integrates multi-source feature enhancement and chaotic optimization. Background Technology

[0002] Electricity load forecasting is a core supporting technology for ensuring the safe and stable operation of the power system and improving the efficiency of energy resource allocation. With the increasing penetration of renewable energy, the diversification of user electricity consumption behavior, the frequent occurrence of extreme weather, and the frequent start-ups and shutdowns of industrial equipment (such as port cranes and industrial park refrigeration units), load sequences exhibit five typical characteristics: strong non-stationarity, multi-periodicity, weather sensitivity, sudden fluctuations, and feature redundancy. Existing forecasting methods still have many shortcomings that urgently need to be addressed: First, the static historical window is not sufficiently adapted. Traditional methods often use fixed-length historical data segments to extract features, which makes it difficult to dynamically match the characteristics of load caused by seasonal changes and equipment start-ups and shutdowns. They are also susceptible to interference from outdated data, resulting in low utilization of effective samples. Second, the fusion of multi-source features is insufficient. Existing technologies often simply splice together weather and load correlation features, failing to accurately select core effective features and neglecting the differentiated mining of linear and nonlinear correlations. Redundant information interferes with model training. Third, the capture of load mutations is lacking. The sudden load changes caused by equipment start-ups and shutdowns and extreme weather are not fully considered, and single frequency decomposition cannot meet the modeling needs of sudden fluctuations and regular trends. Fourth, the integrated optimization strategy is rigid. Many methods use fixed weights to fuse the prediction results of multiple models. Hyperparameter adjustment relies on human experience, which easily leads to local optima and makes it difficult to balance prediction accuracy and stability. In summary, there is an urgent need for a prediction method that can accurately screen multi-source features, capture sudden load changes, dynamically adapt to scenarios, and optimize integration effects in order to improve the reliability of predictions under complex working conditions. Summary of the Invention

[0003] The purpose of this invention is to address the shortcomings of existing technologies by providing a load forecasting method that integrates multi-source feature enhancement and chaotic optimization. This method integrates multi-source data through a triple feature screening mechanism to accurately identify linear and nonlinear core features and eliminate redundant information; it integrates multi-frequency features from wavelet decomposition with iterative cumulative square sum mutation point detection to achieve collaborative modeling of load trends, periodic fluctuations, and sudden changes; it uses a tent chaos-improved multi-objective Jaya optimization algorithm to achieve adaptive adjustment of integrated weights and hyperparameters, avoiding manual intervention and local optima; and it achieves dynamic adaptation based on comprehensive similarity day selection to improve the prediction generalization ability under multiple operating conditions.

[0004] The objective of this invention is achieved through the following technical solutions: A load forecasting method integrating multi-source feature enhancement and chaotic optimization is characterized by the following steps: S1: Collect hourly load data, meteorological data, and equipment operating status data for the target area, and preprocess the data; S2: Using Pearson correlation coefficient, MIC and LASSO The triple feature selection mechanism selects the core feature set from the preprocessed data; S3: Perform wavelet decomposition on the time series in the core feature set to obtain multiple frequency components; use ICSS The algorithm detects mutation points in the sequence, generates mutation point features, and integrates the mutation point features into the corresponding frequency components. After reconstruction, an enhanced feature set is formed. S4: Calculate the comprehensive similarity between the target prediction date and historical dates, select a sample set of similar dates based on the comprehensive similarity, and load the corresponding historical best model parameters to initialize the current prediction model; S5: Construct an integrated prediction model, which includes features for processing full-dimensional features. Transformer Branches, used to capture local features CNN Branching, a simplification for handling linear trends SVR Branching and improvements for strengthening temporal dependencies Mogrifier-LSTM Branching; the integrated weights of each branch model are optimized using a multi-objective Jaya optimization algorithm based on tent chaotic mapping. S6: Based on the optimized fusion weights, the prediction results of each branch model are fused, and the fusion results are inversely normalized to obtain the final load prediction value.

[0005] Step S1 specifically includes the following steps: Based on hourly load data, meteorological data, and equipment operating status data of the target area, the time field is parsed to generate hourly and date dimension information, and the load and meteorological fields are extracted to form the original dataset. X raw For the original dataset X raw Outliers in the data are represented by 3 σ The criteria are used for identification and correction, and the mean of each data column is calculated. μ and standard deviation σ This will exceed [ μ- 3 σ, μ+ 3 σ Outliers in the interval are replaced with the mean of data from two adjacent time points. The corrected formula is as follows: = ( ∉[ μ- 3 σ, μ+ 3σ ]); In the formula, For the original dataset X raw middle Abnormal data values ​​at any given time. , The original datasets are respectively X raw middle Normal data values ​​for the time before and the time after the current moment; For the original dataset X raw Missing values ​​in the data are filled using a combination of forward and backward padding to ensure the continuity of the time series data. The modified and padded dataset is normalized to the [0,1] interval using min-max normalization. The normalization formula is as follows: = ; In the formula, The normalized value. , These are the minimum and maximum values ​​of the feature columns in the corrected and filled dataset, respectively. The data values ​​of the feature columns in the corrected and filled dataset; the minimum and maximum parameters of the feature normalizer and the label normalizer are saved synchronously; Derived features are constructed based on the time dimension, forming a preprocessed dataset containing both original and derived features. X pre .

[0006] Step S2 specifically includes the following steps: Computational preprocessing of datasets X pre Pearson correlation coefficient between each feature and the predicted label r The predicted label is the predicted load, and the calculation formula is as follows: ; In the formula, x t , y t Preprocessed datasets X pre Feature sequences and label sequences in t Data values ​​at any given time , Preprocessed datasets X pre The mean of the feature sequence and the label sequence, nThis represents the number of data samples. Set a filtering threshold so that | r |≥0.3, retain features with significant linear correlations to form a preliminary feature set. X step1 ; use MINE The toolkit calculates the initial feature set. X step1 Features and tags MIC The details are as follows: Calculate the initial screening feature set X step1 Feature sequences X With label sequence Y joint probability distribution P ( X , Y ), characteristic sequence X Marginal probability distribution P ( X ), tag sequence Y Marginal probability distribution P ( Y Mutual information value MI ( X , Y The calculation formula is: ; In the formula, x , y They are the characteristic sequences. X With label sequence Y Data values, P ( x , y )for x , y The joint probability, P ( x )for x The marginal probability, P ( y )for y The marginal probability; Mutual information value MI ( X , Y The maximum information coefficient is obtained by normalization. MIC ( X , Y The formula is: ; In the formula, H ( X ) is a characteristic sequence X Information entropy ,H ( Y ) is a label sequence Y Information entropy ; Set a filtering threshold to make MIC ≥ 0.2, retaining features with significant nonlinear correlations to form a secondary selection feature set. X step2 ; Secondary screening feature set X step2 As input features, construct LASSO Regularization models use the absolute value of penalty coefficients to compress and remove redundant features. LASSO The objective function of the regularization model is as follows: ; In the formula, α For regularization parameters, β 0 represents the intercept term. β i For secondary screening of feature sets X step2 The i The regression coefficients of each feature sequence. p For secondary screening of feature sets X step2 Feature dimensions, y t For secondary screening of feature sets X step2 Tag sequence in t Data values ​​at any given time x ti For secondary screening of feature sets X step2 The i The feature sequences are in t Data values ​​at any given time; Features with absolute regression coefficients greater than 1e-5 are retained, forming a core feature set of 10-15 dimensions. X core The training and test sets are randomly divided in an 8:2 ratio for model training and parameter optimization.

[0007] Step S3 specifically includes the following steps: For core feature set X core The load sequence and meteorological sequence in the data are adopted. db Four wavelet bases are used for 8-level decomposition, based on Mallat The algorithm decomposes the original sequence into a set of low-frequency approximation coefficients. cA Groups 8 and 7 high-frequency detail coefficients cD 1 ~cD 7. The decomposition process can be represented as: ; In the formula, x The time series sequence to be decomposed is... coeffs The set of coefficients obtained from the decomposition. pywt.wavedec ( · ) is the core function in the PyWavelets library used to perform one-dimensional multi-scale wavelet decomposition; Based on frequency characteristics, the decomposition coefficients are classified into three types of frequency components: long-term trend term. cA 8. Medium-term cycle items cD 4 ~ cD 6. Short-term fluctuations cD 1 ~cD 2; Calculate the iterative residual of the sequence a t With cumulative sum of squares C t The specific calculation formula is as follows: ; in, For the front t The mean of the time series, For the front t Variance of the time series; Define mutation point detection statistics ICSS ,when ICSS> At 1.96, this moment is determined to be a breakout point. The statistical formula is as follows: ; In the formula, m The total length of the sequence. C m This is the cumulative sum of squares for the entire sequence from the start point to the end point; Based on the detected mutation point locations, a mutation point feature vector is generated. change future The abrupt change features are incorporated into the decomposition coefficients of the corresponding frequency components at a ratio of 0.1 to 0.3 times the standard deviation. An enhanced feature set fusing multiple frequencies and abrupt change information is then reconstructed using inverse wavelet transform. X enhance The inverse transform process can be represented as: ; In the formula, coeffs fusion This is the set of wavelet coefficients after incorporating the features of abrupt change points. pywt.waverec ( · ) is a multi-level one-dimensional discrete wavelet inverse transform function in the PyWavelets library.

[0008] Step S4 specifically includes the following steps: Meteorological characteristics and load characteristics are selected as the core indicators for the evaluation of similar days. Among them, meteorological characteristics reflect the similarity of the external environment, and load characteristics reflect the similarity of internal time series characteristics. A comprehensive similarity calculation is performed, and the Euclidean distance of meteorological characteristics is calculated between the target prediction date and each historical date. d 1 and load characteristics European distance d 2. The calculation formula is as follows: ; In the formula, f target，i , f hist，i The target date and the historical date are respectively the first i Meteorological characteristics, g target，j , g hist，j The target date and the historical date are respectively the first j Item load characteristics; After normalizing the distance index, the overall similarity is calculated using an equal-weight fusion method. S The formula is as follows: ; In the formula, max( d 1) max() d 2) These are the maximum values ​​of the corresponding distance indicators; Parameter adaptation was performed on the similar day dataset, based on comprehensive similarity. S Sort the samples in descending order and select the top-50 similar days to build an adaptive training dataset. Load the historical best model parameters corresponding to this similar day dataset to initialize the current prediction model and ensure that the model quickly converges to the optimal state in dynamically changing scenarios.

[0009] Step S5 specifically includes the following steps: Build Transformer , CNN ,simplify SVR ,improve Mogrifier-LSTM The four-branch model architecture divides the modeling work according to frequency feature adaptability: Transformer Branch: Through Linear layer, LayerNorm Layers and GELU Activation function, integrating and enhancing feature sets X enhance The full-dimensional cross-frequency correlation information is used to output the sequence mean as the branch prediction result. out trans The calculation process can be expressed as: out trans =mean ( GELU ( LayerNorm ( Linear ( X enhance )))) In the formula, mean ( · This indicates the operation of averaging along the sequence dimension; CNN Branch: Uses convolutional kernels of size 3 and number 64. Conv 1 d layer, BatchNorm 1 d Layers and GELU The activation function captures the local nonlinear characteristics of short-term fluctuations, and the branch prediction results are output after global average pooling. out cnn The calculation process can be expressed as: out cnn =GlobalAvgPool ( GELU ( BatchNorm 1 d ( Conv 1 d ( X enhance ,3,64))))). In the formula, GlobalAvgPool Indicates global average pooling; simplify SVR Branch: Through Linear Layers and LayerNorm After layer preprocessing of features, a support vector regression model is used to adapt the linear correlation features of the long-term trend term, and the branch prediction results are output. out svr The calculation process can be expressed as: out svr =SVR ( LayerNorm ( Linear(flatten ( X enhance )))) In the formula, flatten Indicates the flattening operation; improve Mogrifier-LSTM Branching: Employs 3 rounds of input-hidden state interaction. Mogrifier The mechanism enhances the ability to capture long-term dependencies through bidirectional means. LSTM Layered modeling of the temporal correlation of mid-term periodic terms, with the output sequence mean used as the branch prediction result. out LSTM The calculation process can be expressed as: out LSTM =mean ( BiLSTM ( Mogrifier ( X enhance , hidden ,3)))). In the formula, mean ( · () indicates taking the mean along the sequence dimension. BiLSTM bidirectional LSTM layer, hidden To hide the state dimension; The model branches are jointly trained based on the training set and the similar day adaptation dataset until the validation set loss converges. The ensemble weights of the four-branch model are optimized using a multi-objective Jaya optimization algorithm based on tent chaotic mapping, as detailed below: Population initialization: An initial population is generated through tent chaotic mapping. The individual dimensions of the population are 4, corresponding to the ensemble weights of the four models. ω 1. ω 2. ω 3. ω 4. The chaotic mapping formula is as follows: x i =low + α ·( high - low )·(1-2| r -0.5|); In the formula, x i For the process i The chaotic variable values ​​after the next iteration r A random number in the interval [0,1]. α Here are the chaos coefficients. α = 0.4, low , high These are the lower and upper bounds of the population boundary, respectively. low = -1, high = 1; Multi-objective optimization: Using the mean absolute error of the validation set and the prediction variance as optimization objectives, the population update mechanism of the multi-objective Jayer optimization algorithm is used to iteratively update individuals, select the Pareto optimal solution, and ensure that the ensemble weights satisfy the following: ω 1+ ω 2+ ω 3+ ω 4 = 1; Ensemble prediction: Based on the optimized ensemble weights, the weighted sum of the prediction results of the four-branch model is calculated to obtain the final ensemble prediction result. outensemble The calculation formula is as follows: out ensemble =ω 1· out trans + ω 2. out cnn + ω 3· out svr + ω 4. out LSTM .

[0010] Step S6 specifically includes the following steps: The parameters of the label normalizer saved in the preprocessing stage are called to standardize the prediction results. out ensemble Inverse transformation to load forecast values ​​with physical units y pred The formula is: y pred =out ensemble ·( y max - y min )+ y min ; In the formula, y max , y min These are the maximum and minimum values ​​of the label sequence, respectively.

[0011] The advantages of this invention are: 1. In feature selection, "Pearson+" MIC + LASSO "The triple mechanism increases the utilization rate of core features by 50%, the rate of redundant information removal by 40%, and the efficiency of model training by more than 35%." 2. In mutation capture, "wavelet decomposition + ICSS "The feature enhancement system improves the accuracy of capturing load surges caused by extreme weather and equipment start-up and shutdown by 30% to 35%, and reduces prediction errors by 25% to 30% under extreme operating conditions." 3. In terms of optimization efficiency, the chaotic multi-objective Jaya optimization algorithm reduces the time spent on integrated weight and hyperparameter optimization by 40%, improves the balance between prediction accuracy and stability by 45%, and avoids the subjectivity and local optima of manual parameter tuning. 4. Regarding scenario adaptability, the comprehensive similarity daily dynamic adaptation mechanism ensures that the model's average absolute percentage error (ASE) varies across multiple scenarios, including normal operating conditions, extreme weather, and intensive equipment start-ups and shutdowns.MAPE The accuracy is stable within 3% to 5%, while traditional single models often have errors exceeding 10% in extreme scenarios; 5. In terms of engineering practicality, based on routine monitoring data, no additional hardware is required, single-sample prediction time is < 80ms, and it supports... CPU / GPU Automatic adaptation, can be directly integrated into Python Environment or encapsulation as API The model parameters are automatically optimized through data, reducing operation and maintenance costs and application thresholds. It also has wide functional compatibility and can flexibly adapt to multiple tasks such as load forecasting and energy consumption forecasting, and has a wide range of engineering application scenarios. Attached Figure Description

[0012] Figure 1 This is the overall flowchart of the load prediction method integrating multi-source feature enhancement and chaotic optimization of the present invention; Figure 2 This is a block diagram of the load prediction method integrating multi-source feature enhancement and chaotic optimization of the present invention. Figure 3 This diagram illustrates the working principle of each module used in the load prediction method that integrates multi-source feature enhancement and chaotic optimization in this invention. Detailed Implementation

[0013] The features and other related features of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments, so as to facilitate understanding by those skilled in the art: Example: Figures 1 - 3 As shown, this embodiment relates to a load forecasting method that integrates multi-source feature enhancement and chaotic optimization. The method mainly includes the following steps: S1: Collect hourly load data, meteorological data, and equipment operation status data for the target area (all data time granularity is 1 hour, data time span is not less than 365 days, to ensure the integrity of time series characteristics), and preprocess the data.

[0014] In this embodiment, step S1 is implemented through a data preprocessing module. Step S1 specifically includes the following steps: Based on hourly load data (obtainable through the power grid dispatch platform), meteorological data (including temperature, humidity, irradiance, and wind speed, obtainable through the meteorological monitoring platform), and equipment operating status data (including power factor and whether the equipment is out of service) for the target area, the time field is parsed to generate hourly and date-dimensional information, and the load and meteorological fields are extracted to form the original dataset. X raw For the original dataset X raw Outliers in the data are represented by 3 σ The criteria are used for identification and correction, and the mean of each data column is calculated. μ and standard deviationσ This will exceed [ μ- 3 σ, μ+ 3 σ Outliers in the interval are replaced with the mean of data from two adjacent time points. The corrected formula is as follows: = ( ∉[ μ- 3 σ, μ+ 3 σ ]); In the formula, For the original dataset X raw middle t Abnormal data values ​​at any given time. , The original datasets are respectively X raw middle t Normal data values ​​for the time before and the time after the current moment.

[0015] For the original dataset X raw Missing values ​​are filled using a combination of forward and backward padding to ensure the continuity of time series data.

[0016] To eliminate the influence of different units on model training, min-max normalization is adopted. MinMaxScaler The corrected and padded dataset is normalized to the [0,1] interval using the following normalization formula: = ; In the formula, The normalized value. , These are the minimum and maximum values ​​of the feature columns in the corrected and filled dataset, respectively. The data values ​​of the feature columns in the dataset after correction and filling; the feature normalizer is saved synchronously. scaler X ) and label normalizer ( scaler y The minimum and maximum parameters of ) are used for the inverse normalization of subsequent prediction results; Based on the time dimension, derived features (including hourly periodic features, daily average load, and daily load peak-to-valley difference) are constructed to form a preprocessed dataset containing both original and derived features. X pre .

[0017] S2: Using Pearson correlation coefficient, MIC ( Maximal Information Coefficient (maximum information coefficient) and LASSO ( Least Absolute Shrinkage and Selection Operator The triple feature selection mechanism (minimum absolute shrinkage and selection operator) selects the core feature set from the preprocessed data.

[0018] In this embodiment, step S2 is implemented through a multi-dimensional feature filtering module. Step S2 specifically includes the following steps: 1. Linear screening using Pearson correlation coefficient: Computational preprocessing of datasets X pre Pearson correlation coefficient between each feature and the predicted label r The predicted label is the predicted load, and the calculation formula is as follows: ; In the formula, x t , y t Preprocessed datasets X pre Feature sequences and label sequences in t Data values ​​at any given time , Preprocessed datasets X pre The mean of the feature sequence and the label sequence, n This represents the number of data samples. Set a filtering threshold so that | r |≥0.3, retain features with significant linear correlations to form a preliminary feature set. X step1 .

[0019] 2. MIC Non-linear filtering: use MINE The (Maximal Information-based Nonparametric Exploration) toolkit calculates the initial feature set. X step1 Features and tags MIC This indicator effectively captures nonlinear associations by quantifying the strength of the maximum information correlation between two variables, and is not limited by the type of variable distribution. Its core calculation formula is based on the definition of information entropy, as follows: Calculate the initial screening feature set X step1 Feature sequences X With label sequence Yjoint probability distribution P ( X , Y ), characteristic sequence X Marginal probability distribution P ( X ), tag sequence Y Marginal probability distribution P ( Y Mutual information value MI ( X , Y The calculation formula is: ; In the formula, x , y They are the characteristic sequences. X With label sequence Y Data values, P ( x , y )for x , y The joint probability, P ( x )for x The marginal probability, P ( y )for y The marginal probability; Mutual information value MI ( X , Y The maximum information coefficient is obtained by normalization. MIC ( X , Y The formula is: ; In the formula, H ( X ) is a characteristic sequence X Information entropy , H ( Y ) is a label sequence Y Information entropy ; Set a filtering threshold to make MIC ≥ 0.2, retaining features with significant nonlinear correlations to form a secondary selection feature set. X step2 .

[0020] 3. LASSO Regularized redundancy compression: Secondary screening feature set X step2 As input features, construct LASSORegularization models use the absolute value of penalty coefficients to compress and remove redundant features. LASSO The objective function of the regularization model is as follows: ; In the formula, α For regularization parameters, α= 0.01, β 0 represents the intercept term. β i For secondary screening of feature sets X step2 The i The regression coefficients of each feature sequence. p For secondary screening of feature sets X step2 Feature dimensions, y t For secondary screening of feature sets X step2 Tag sequence in t Data values ​​at any given time x ti For secondary screening of feature sets X step2 The i The feature sequences are in t Data values ​​at any given time; Features with absolute regression coefficients greater than 1e-5 are retained, forming a core feature set of 10-15 dimensions. X core The training and test sets are randomly divided in an 8:2 ratio for model training and parameter optimization.

[0021] S3: Perform wavelet decomposition on the time series in the core feature set to obtain multiple frequency components; use ICSS ( Iterated Cumulative Sums of Squares The iterative cumulative sum of squares (ITS) algorithm detects mutation points in the sequence, generates mutation point features, and integrates the mutation point features into the corresponding frequency components. After reconstruction, an enhanced feature set is formed.

[0022] In this embodiment, step S3 is implemented through a multi-dimensional feature enhancement module. Step S3 specifically includes the following steps: 1. Wavelet decomposition for multi-frequency component extraction: For core feature set X core The load sequence and meteorological sequence in the data are adopted. db Four wavelet bases are used for 8-level decomposition, based on Mallat The algorithm (a multi-resolution signal processing method based on wavelet transform) decomposes the original sequence into a set of low-frequency approximation coefficients. cA Groups 8 and 7 high-frequency detail coefficients cD 1 ~cD7. The decomposition process can be represented as: ; In the formula, x The time series sequence to be decomposed is... coeffs The set of coefficients obtained from the decomposition. pywt.wavedec ( · )for PyWavelets Library (a type of library) Python The core function for performing one-dimensional multi-scale wavelet decomposition in a free and open-source library for discrete wavelet transform and correlation analysis; Based on frequency characteristics, the decomposition coefficients are classified into three types of frequency components: long-term trend term. cA 8 (corresponding to low-frequency, slow-changing characteristics), medium-term periodic term cD 4 ~cD 6 (corresponding to mesoscale cyclical fluctuations), short-term fluctuation term cD 1 ~cD 2 (corresponding to high-frequency sudden fluctuations).

[0023] 2. ICSS Mutation point detection and feature fusion: Calculate the iterative residual of the sequence a t With cumulative sum of squares C t The specific calculation formula is as follows: ; in, For the front t The mean of the time series, For the front t Variance of the time series; Define mutation point detection statistics ICSS ,when ICSS> At 1.96, this moment is determined to be a breakout point. The statistical formula is as follows: ; In the formula, m The total length of the sequence. C m This is the cumulative sum of squares for the entire sequence from the start point to the end point; Based on the detected mutation point locations, a mutation point feature vector is generated. change future (The mutation point is marked as 1, and all other times are marked as 0). The mutation point features are incorporated into the decomposition coefficients of the corresponding frequency components at a ratio of 0.2 times the standard deviation. The enhanced feature set, which integrates multi-frequency and mutation point information, is then reconstructed through inverse wavelet transform. X enhance The inverse transform process can be represented as: ; In the formula, coeffs fusion This is the set of wavelet coefficients after incorporating the features of abrupt change points. pywt.waverec ( · )for PyWavelets The library contains multi-level one-dimensional discrete wavelet inverse transform functions.

[0024] S4: Calculate the comprehensive similarity between the target prediction date and historical dates, select a sample set of similar dates based on the comprehensive similarity, and load the corresponding historical best model parameters to initialize the current prediction model.

[0025] In this embodiment, step S4 is implemented through the integrated similar day dynamic adaptation module. Step S4 specifically includes the following steps: Meteorological characteristics (including temperature, humidity, irradiance, and wind speed) and load characteristics (including hourly cycle characteristics, daily average load, and daily peak-valley difference in load) are selected as the core indicators for evaluating similar days. Among them, meteorological characteristics reflect the similarity of the external environment, while load characteristics reflect the similarity of internal temporal characteristics. A comprehensive similarity calculation is performed, and the Euclidean distance of meteorological characteristics is calculated between the target prediction date and each historical date. d 1 and load characteristics European distance d 2. The calculation formula is as follows: ; In the formula, f target，i , f hist，i The target date and the historical date are respectively the first i Meteorological characteristics, g target，j , g hist，j The target date and the historical date are respectively the first j Item load characteristics; After normalizing the distance index, the overall similarity is calculated using an equal-weight fusion method. S The formula is as follows: ; In the formula, max( d 1) max() d 2) These are the maximum values ​​of the corresponding distance indicators; Parameter adaptation was performed on the similar day dataset, based on comprehensive similarity. S Sort the samples in descending order and select the top-50 similar days to construct an adaptation training dataset (including a similar day feature dataset). X sim Similar date tag datasets Y simThe system loads the historical best model parameters (including initial weights and learning rates of multiple models) corresponding to the similar day dataset to initialize the current prediction model, ensuring that the model quickly converges to the optimal state in dynamically changing scenarios.

[0026] S5: Construct an ensemble prediction model, which includes features for handling full-dimensional features. Transformer (Transformer) Branch, used to capture local features CNN ( Convolutional Neural Network (Convolutional Neural Network) branch, a simplified approach for handling linear trends SVR ( Support Vector Regression Support vector regression (SVR) branch and improvements for strengthening temporal dependencies Mogrifier-LSTM (integrated) Mogrifier Mechanism LSTM ) branch; adopting a multi-objective Jaya (based on tent chaotic mapping) improved version. MOJAYA Optimize the algorithm and optimize the ensemble weights of each branch model.

[0027] In this embodiment, step S5 is implemented through a chaos optimization multi-model ensemble module. Step S5 specifically includes the following steps: Build Transformer , CNN ,simplify SVR ,improve Mogrifier-LSTM The four-branch model architecture divides the modeling work according to frequency feature adaptability: 1. Transformer Branch: Through Linear Layer (linear projection layer) LayerNorm Layer (normalized layer) and GELU (Gaussian error linear unit) activation function, integrating and enhancing feature sets. X enhance The full-dimensional cross-frequency correlation information is used to output the sequence mean as the branch prediction result. out trans The calculation process can be expressed as: out trans =mean ( GELU ( LayerNorm ( Linear ( X enhance )))) In the formula, mean ( · ) indicates the operation of averaging along the sequence dimension.

[0028] 2. CNN Branch: Uses convolutional kernels of size 3 and number 64. Conv 1 dLayer (one-dimensional convolutional layer) BatchNorm 1 d Layer (one-dimensional batch normalization layer) and GELU The activation function captures the local nonlinear characteristics of short-term fluctuations, and the branch prediction results are output after global average pooling. out cnn The calculation process can be expressed as: out cnn =GlobalAvgPool ( GELU ( BatchNorm 1 d ( Conv 1 d ( X enhance ,3,64))))). In the formula, GlobalAvgPool This indicates global average pooling.

[0029] 3. Simplify SVR Branch: Through Linear Layers and LayerNorm After layer preprocessing of features, a support vector regression model is used to adapt the linear correlation features of the long-term trend term, and the branch prediction results are output. out svr The calculation process can be expressed as: out svr =SVR ( LayerNorm ( Linear(flatten ( X enhance )))) In the formula, flatten This indicates a flattening operation.

[0030] 4. Improvement Mogrifier-LSTM Branching: Employs 3 rounds of input-hidden state interaction. Mogrifier Mechanism (to) LSTM / RNN Enhancement techniques) enhance the ability to capture long-term time-series dependencies through bidirectional... LSTM (Long Short-Term Memory Network) layers model the temporal correlation of medium-term periodic terms, and the output sequence mean is used as the branch prediction result. out LSTM The calculation process can be expressed as: out LSTM =mean ( BiLSTM ( Mogrifier ( X enhance , hidden ,3)))). In the formula, mean ( ·() indicates taking the mean along the sequence dimension. BiLSTM bidirectional LSTM layer, hidden To hide the state dimension.

[0031] The model branches are jointly trained based on the training set and the similar day adaptation dataset until the validation set loss converges. The ensemble weights of the four-branch model are optimized using a multi-objective Jaya optimization algorithm based on tent chaotic mapping, as detailed below: 1. Population Initialization: An initial population is generated through tent chaotic mapping. The individual dimensions of the population are 4, corresponding to the ensemble weights of the four models. ω 1. ω 2. ω 3. ω 4. The chaotic mapping formula is as follows: x i =low+α ·( high-low )·(1-2| r -0.5|); In the formula, x i For the process i The chaotic variable values ​​after the next iteration r A random number in the interval [0,1]. α Here are the chaos coefficients. α = 0.4, low , high These are the lower and upper bounds of the population boundary, respectively. low= -1, high= 1.

[0032] 2. Multi-objective optimization: using the validation set mean absolute error (MAE) MAE Using the prediction variance as the optimization objective, individuals are iteratively updated through the population update mechanism (crossover-perturbation strategy) of the multi-objective Jaya optimization algorithm to select the Pareto optimal solution and ensure that the ensemble weights satisfy: ω 1+ ω 2+ ω 3+ ω 4 = 1.

[0033] 3. Ensemble Prediction: Based on the optimized ensemble weights, the weighted sum of the prediction results of the four-branch model is calculated to obtain the final ensemble prediction result. out ensemble The calculation formula is as follows: out ensemble =ω 1· out trans + ω 2. out cnn + ω 3· out svr + ω 4. out LSTM .

[0034] S6: Based on the optimized fusion weights, the prediction results of each branch model are fused, and the fusion results are inversely normalized to obtain the final load prediction value.

[0035] In this embodiment, step S6 is implemented through the prediction output module. Step S6 specifically includes the following steps: Call the label normalizer saved in the preprocessing stage ( scaler y The parameters of ) will standardize the prediction results. out ensemble Inverse transformation to load forecast values ​​with physical units y pred The formula is: y pred =out ensemble ·( y max - y min )+ y min ; In the formula, y max , y min These are the maximum and minimum values ​​of the label sequence, respectively.

[0036] The beneficial technical effects of this embodiment are as follows: 1. In feature selection, "Pearson+" MIC + LASSO "The triple mechanism increases the utilization rate of core features by 50%, the rate of redundant information removal by 40%, and the efficiency of model training by more than 35%." 2. In mutation capture, "wavelet decomposition + ICSS "The feature enhancement system improves the accuracy of capturing load surges caused by extreme weather and equipment start-up and shutdown by 30% to 35%, and reduces prediction errors by 25% to 30% under extreme operating conditions." 3. In terms of optimization efficiency, the chaotic multi-objective Jaya optimization algorithm reduces the time spent on integrated weight and hyperparameter optimization by 40%, improves the balance between prediction accuracy and stability by 45%, and avoids the subjectivity and local optima of manual parameter tuning. 4. Regarding scenario adaptability, the comprehensive similarity daily dynamic adaptation mechanism ensures that the model's average absolute percentage error (ASE) varies across multiple scenarios, including normal operating conditions, extreme weather, and intensive equipment start-ups and shutdowns. MAPE The accuracy is stable within 3% to 5%, while traditional single models often have errors exceeding 10% in extreme scenarios; 5. In terms of engineering practicality, based on routine monitoring data, no additional hardware is required, single-sample prediction time is < 80ms, and it supports... CPU / GPU Automatic adaptation, can be directly integrated into Python Environment or encapsulation as API The model parameters are automatically optimized through data, reducing operation and maintenance costs and application thresholds. It also has wide functional compatibility and can flexibly adapt to multiple tasks such as load forecasting and energy consumption forecasting, and has a wide range of engineering application scenarios.

Claims

1. A load forecasting method integrating multi-source feature enhancement and chaotic optimization, characterized in that... The method includes the following steps: S1: Collect hourly load data, meteorological data, and equipment operating status data for the target area, and preprocess the data; S2: Using Pearson correlation coefficient, MIC and LASSO The triple feature selection mechanism selects the core feature set from the preprocessed data; S3: Perform wavelet decomposition on the time series in the core feature set to obtain multiple frequency components; use ICSS The algorithm detects mutation points in the sequence, generates mutation point features, and integrates the mutation point features into the corresponding frequency components. After reconstruction, an enhanced feature set is formed. S4: Calculate the comprehensive similarity between the target prediction date and historical dates, select a sample set of similar dates based on the comprehensive similarity, and load the corresponding historical best model parameters to initialize the current prediction model; S5: Construct an integrated prediction model, which includes features for processing full-dimensional features. Transformer Branches, used to capture local features CNN Branching, a simplification for handling linear trends SVR Branching and improvements for strengthening temporal dependencies Mogrifier-LSTM Branching; the integrated weights of each branch model are optimized using a multi-objective Jaya optimization algorithm based on tent chaotic mapping. S6: Based on the optimized fusion weights, the prediction results of each branch model are fused, and the fusion results are inversely normalized to obtain the final load prediction value.

2. The load forecasting method integrating multi-source feature enhancement and chaotic optimization as described in claim 1, characterized in that... Step S1 specifically includes the following steps: Based on hourly load data, meteorological data, and equipment operating status data of the target area, the time field is parsed to generate hourly and date dimension information, and the load and meteorological fields are extracted to form the original dataset. X raw For the original dataset X raw Outliers in the data are represented by 3 σ The criteria are used for identification and correction, and the mean of each data column is calculated. μ and standard deviation σ This will exceed [ μ- 3 σ, μ+ 3 σ Outliers in the interval are replaced with the mean of data from two adjacent time points. The corrected formula is as follows: = （ ∉[ μ- 3 σ, μ+ 3 σ ]）； In the formula, For the original dataset X raw middle Abnormal data values ​​at any given time. , The original datasets are respectively X raw middle Normal data values ​​for the time before and the time after the current moment; For the original dataset X raw Missing values ​​in the data are filled using a combination of forward and backward padding to ensure the continuity of the time series data. The modified and padded dataset is normalized to the [0,1] interval using min-max normalization. The normalization formula is as follows: = ； In the formula, The normalized value. , These are the minimum and maximum values ​​of the feature columns in the corrected and filled dataset, respectively. The data values ​​of the feature columns in the corrected and filled dataset; the minimum and maximum parameters of the feature normalizer and the label normalizer are saved synchronously; Derived features are constructed based on the time dimension, forming a preprocessed dataset containing both original and derived features. X pre .

3. The load forecasting method integrating multi-source feature enhancement and chaotic optimization as described in claim 2, characterized in that... Step S2 specifically includes the following steps: Computational preprocessing of datasets X pre Pearson correlation coefficient between each feature and the predicted label r The predicted label is the predicted load, and the calculation formula is as follows: ； In the formula, x t , y t Preprocessed datasets X pre Feature sequences and label sequences in t Data values ​​at any given time , Preprocessed datasets X pre The mean of the feature sequence and the label sequence, n This represents the number of data samples. Set a filtering threshold so that | r |≥0.3, retain features with significant linear correlations to form a preliminary feature set. X step1 ; use MINE The toolkit calculates the initial feature set. X step1 Features and tags MIC The details are as follows: Calculate the initial screening feature set X step1 Feature sequences X With label sequence Y joint probability distribution P ( X , Y ), characteristic sequence X Marginal probability distribution P ( X ), tag sequence Y Marginal probability distribution P ( Y Mutual information value MI ( X , Y The calculation formula is: ； In the formula, x , y They are the characteristic sequences. X With label sequence Y Data values, P ( x , y )for x , y The joint probability, P ( x )for x The marginal probability, P ( y )for y The marginal probability; Mutual information value MI ( X , Y The maximum information coefficient is obtained by normalization. MIC ( X , Y The formula is: ； In the formula, H ( X ) is a characteristic sequence X Information entropy , H ( Y ) is a label sequence Y Information entropy ; Set a filtering threshold to make MIC ≥ 0.2, retaining features with significant nonlinear correlations to form a secondary selection feature set. X step2 ; Secondary screening feature set X step2 As input features, construct LASSO Regularization models use the absolute value of penalty coefficients to compress and remove redundant features. LASSO The objective function of the regularization model is as follows: ； In the formula, α For regularization parameters, β 0 represents the intercept term. β i For secondary screening of feature sets X step2 The i The regression coefficients of each feature sequence. p For secondary screening of feature sets X step2 Feature dimensions, y t For secondary screening of feature sets X step2 Tag sequence in t Data values ​​at any given time x ti For secondary screening of feature sets X step2 The i The feature sequences are in t Data values ​​at any given time; Features with absolute regression coefficients greater than 1e-5 are retained, forming a core feature set of 10-15 dimensions. X core The training and test sets are randomly divided in an 8:2 ratio for model training and parameter optimization.

4. The load forecasting method integrating multi-source feature enhancement and chaotic optimization as described in claim 3, characterized in that... Step S3 specifically includes the following steps: For core feature set X core The load sequence and meteorological sequence in the data are adopted. db Four wavelet bases are used for 8-level decomposition, based on Mallat The algorithm decomposes the original sequence into a set of low-frequency approximation coefficients. cA Groups 8 and 7 high-frequency detail coefficients cD 1 ~cD 7. The decomposition process can be represented as: ； In the formula, x The time series sequence to be decomposed is... coeffs The set of coefficients obtained from the decomposition. pywt.wavedec ( · ) is the core function in the PyWavelets library used to perform one-dimensional multi-scale wavelet decomposition; Based on frequency characteristics, the decomposition coefficients are classified into three types of frequency components: long-term trend term. cA 8. Medium-term cycle items cD 4 ~cD 6. Short-term fluctuations cD 1 ~cD 2; Calculate the iterative residual of the sequence a t With cumulative sum of squares C t The specific calculation formula is as follows: ； in, For the front t The mean of the time series, For the front t Variance of the time series; Define mutation point detection statistics ICSS ,when ICSS> At 1.96, this moment is determined to be a breakout point. The statistical formula is as follows: ； In the formula, m The total length of the sequence. C m This is the cumulative sum of squares for the entire sequence from the start point to the end point; Based on the detected mutation point locations, a mutation point feature vector is generated. change future The abrupt change features are incorporated into the decomposition coefficients of the corresponding frequency components at a ratio of 0.1 to 0.3 times the standard deviation. An enhanced feature set fusing multiple frequencies and abrupt change information is then reconstructed using inverse wavelet transform. X enhance The inverse transform process can be represented as: ； In the formula, coeffs fusion This is the set of wavelet coefficients after incorporating the features of abrupt change points. pywt.waverec ( · ) is a multi-level one-dimensional discrete wavelet inverse transform function in the PyWavelets library.

5. The load forecasting method integrating multi-source feature enhancement and chaotic optimization as described in claim 4, characterized in that... Step S4 specifically includes the following steps: Meteorological characteristics and load characteristics are selected as the core indicators for the evaluation of similar days. Among them, meteorological characteristics reflect the similarity of the external environment, and load characteristics reflect the similarity of internal time series characteristics. A comprehensive similarity calculation is performed, and the Euclidean distance of meteorological characteristics is calculated between the target prediction date and each historical date. d 1 and load characteristics European distance d 2. The calculation formula is as follows: ； In the formula, f target，i , f hist，i The target date and the historical date are respectively the first i Meteorological characteristics, g target，j , g hist，j The target date and the historical date are respectively the first j Item load characteristics; After normalizing the distance index, the overall similarity is calculated using an equal-weight fusion method. S The formula is as follows: ； In the formula, max( d 1) max() d 2) These are the maximum values ​​of the corresponding distance indicators; Parameter adaptation was performed on the similar day dataset, based on comprehensive similarity. S Sort the samples in descending order and select the top-50 similar days to build an adaptive training dataset. Load the historical best model parameters corresponding to this similar day dataset to initialize the current prediction model and ensure that the model quickly converges to the optimal state in dynamically changing scenarios.

6. The load forecasting method integrating multi-source feature enhancement and chaotic optimization as described in claim 5, characterized in that... Step S5 specifically includes the following steps: Build Transformer , CNN ,simplify SVR ,improve Mogrifier-LSTM The four-branch model architecture divides the modeling work according to frequency feature adaptability: Transformer Branch: Through Linear layer, LayerNorm Layers and GELU Activation function, integrating and enhancing feature sets X enhance The full-dimensional cross-frequency correlation information is used to output the sequence mean as the branch prediction result. out trans The calculation process can be expressed as: out trans =mean （ GELU （ LayerNorm （ Linear （ X enhance ））））； In the formula, mean ( · This indicates the operation of averaging along the sequence dimension; CNN Branch: Uses convolutional kernels of size 3 and number 64. Conv 1 d layer, BatchNorm 1 d Layers and GELU The activation function captures the local nonlinear characteristics of short-term fluctuations, and the branch prediction results are output after global average pooling. out cnn The calculation process can be expressed as: out cnn =GlobalAvgPool （ GELU （ BatchNorm 1 d （ Conv 1 d （ X enhance ，3,64））））； In the formula, GlobalAvgPool Indicates global average pooling; simplify SVR Branch: Through Linear Layers and LayerNorm After layer preprocessing of features, a support vector regression model is used to adapt the linear correlation features of the long-term trend term, and the branch prediction results are output. out svr The calculation process can be expressed as: out svr =SVR （ LayerNorm （ Linear (flatten) （ X enhance ））））； In the formula, flatten Indicates the flattening operation; improve Mogrifier-LSTM Branching: Employs 3 rounds of input-hidden state interaction. Mogrifier The mechanism enhances the ability to capture long-term dependencies through bidirectional means. LSTM Layered modeling of the temporal correlation of mid-term periodic terms, with the output sequence mean used as the branch prediction result. out LSTM The calculation process can be expressed as: out LSTM =mean （ BiLSTM （ Mogrifier （ X enhance ， hidden ,3）））； In the formula, mean ( · () indicates taking the mean along the sequence dimension. BiLSTM bidirectional LSTM layer, hidden To hide the state dimension; The model branches are jointly trained based on the training set and the similar day adaptation dataset until the validation set loss converges. The ensemble weights of the four-branch model are optimized using a multi-objective Jaya optimization algorithm based on tent chaotic mapping, as detailed below: Population initialization: An initial population is generated through tent chaotic mapping. The individual dimensions of the population are 4, corresponding to the ensemble weights of the four models. ω 1. ω 2. ω 3. ω 4. The chaotic mapping formula is as follows: x i =low+α ·（ high-low ）·（1-2| r -0.5|）； In the formula, x i For the process i The chaotic variable values ​​after the next iteration r A random number in the interval [0,1]. α Here are the chaos coefficients. α= 0.4, low , high These are the lower and upper bounds of the population boundary, respectively. low= -1, high= 1; Multi-objective optimization: Using the mean absolute error of the validation set and the prediction variance as optimization objectives, the population update mechanism of the multi-objective Jayer optimization algorithm is used to iteratively update individuals, select the Pareto optimal solution, and ensure that the ensemble weights satisfy the following: ω 1+ ω 2+ ω 3+ ω 4 = 1； Ensemble prediction: Based on the optimized ensemble weights, the weighted sum of the prediction results of the four-branch model is calculated to obtain the final ensemble prediction result. out ensemble The calculation formula is as follows: out ensemble =ω 1· out trans + ω 2· out cnn + ω 3· out svr + ω 4· out LSTM 。 7. The load forecasting method integrating multi-source feature enhancement and chaotic optimization as described in claim 6, characterized in that... Step S6 specifically includes the following steps: The parameters of the label normalizer saved in the preprocessing stage are called to standardize the prediction results. out ensemble Inverse transformation to load forecast values ​​with physical units y pred The formula is: y pred =out ensemble ·（ y max - y min ）+ y min ； In the formula, y max , y min These are the maximum and minimum values ​​of the label sequence, respectively.

Images