Electric vehicle charging load prediction method and system

Abstract

The invention discloses an electric vehicle charging load prediction method and system, and the method comprises the steps: collecting and obtaining the historical charging load time series data of an electric vehicle, and employing an ICEEMDAN model to carry out the decomposition of the charging load time series data, and obtaining a plurality of components; performing sample entropy calculation on the components, and performing signal classification by adopting K-medoids clustering according to a sample entropy result to obtain a high-frequency component F1, an intermediate-frequency component F2 and a low-frequency component F3; carrying out secondary decomposition on the high-frequency component F1 by adopting a characteristic mode decomposition (FMD) model to obtain a plurality of components; the MTS-Mixers model, the Crossform model and the DeepESN model are optimized by using an improved artemisinin AO optimization algorithm, and a high-frequency component F1, an intermediate-frequency component F2 and a low-frequency component F3 are predicted; the prediction results of the AO-MTSMixers model, the AO-Crossform model and the AO-DeepESN model are reconstructed, and the prediction result of the charging load is obtained; compared with a traditional method, the charging load prediction accuracy can be improved, and powerful support is provided for charging use and management of the electric vehicle.

Description

Technical Field

[0001] The present invention relates to a charging load prediction method and system, and in particular to an electric vehicle charging load prediction method and system, and belongs to the technical field of charging load prediction. Background Art

[0002] With the booming development of the electric vehicle industry and the increasing prominence of environmental issues, electric vehicles are gradually becoming an important choice for future transportation. However, the popularity of electric vehicles has also brought new challenges, one of which is the prediction of charging load. Accurate charging load prediction is of great significance for power grid dispatching, charging facility planning, and electric vehicle charging management. It can help power grid companies to reasonably arrange power resources and avoid overload or underload of the power system. It can also provide scientific planning and construction basis for charging facility operators to ensure that electric vehicle users can obtain charging services conveniently and quickly.

[0003] At present, the commonly used charging load prediction methods mainly include: prediction methods based on statistical models, machine learning, and deep learning. The prediction method based on statistical models mainly uses statistical analysis to derive the changing rules and trends of the charging load. This type of method is simple and easy to implement, but the prediction accuracy is low and it is difficult to capture the complex changing characteristics of the charging load. The prediction method based on machine learning mainly uses machine learning algorithms to analyze and model the charging load data. This type of method has high prediction accuracy, but requires a large amount of training data and computing resources, and the model generalization ability is limited. The prediction method based on deep learning mainly uses deep learning algorithms to automatically learn the complex features in the charging load data. This type of method has high prediction accuracy, but the model is complex, the training is difficult, and a large amount of training data is required.

[0004] Electric vehicle charging load data is affected by many factors, such as vehicle type, driving habits, distribution of charging facilities, weather conditions, etc., which makes the charging load data present complex nonlinear characteristics and difficult to predict directly. In addition, existing prediction models often ignore the characteristics of different frequency components in the charging load data, resulting in insufficient prediction accuracy. In order to overcome the above defects, it is necessary to further study and develop new electric vehicle charging load prediction methods, such as hybrid prediction methods that combine multiple models to improve prediction accuracy and adaptability to meet practical application needs. Summary of the invention

[0005] Purpose of the invention: The purpose of the present invention is to provide a method and system for predicting electric vehicle charging load that can improve prediction accuracy.

[0006] Technical solution: The electric vehicle charging load prediction method described in the present invention comprises:

[0007] (1) Collecting and obtaining historical charging load time series data of electric vehicles, and using the ICEEMDAN model to decompose the charging load time series data into several components;

[0008] (2) performing sample entropy calculation on the components, and performing signal classification using K-medoids clustering according to the sample entropy results to obtain a high-frequency component F1, a mid-frequency component F2, and a low-frequency component F3;

[0009] (3) The high-frequency component F1 is decomposed twice using the characteristic mode decomposition (FMD) model to obtain several components;

[0010] (4) The improved artemisinin AO optimization algorithm is used to optimize the time step μ, number of channels ν and prediction step ρ in the MTS-Mixers model to obtain the AO-MTS_Mixers model, and the AO-MTS_Mixers model is used to batch predict the component signals of the high-frequency component F1; the improved artemisinin AO optimization algorithm is used to optimize the data partition length L in the Crossformer model. seg , the number of layers l and the total dimension D are optimized to obtain the AO-Crossformer model, and the AO-Crossformer model is used to predict the intermediate frequency component F2; the improved artemisinin AO optimization algorithm is used to optimize the storage layer leakage rate γ and the storage layer size P in the DeepESN model to obtain the AO-DeepESN model, and the AO-DeepESN model is used to predict the low frequency component F3;

[0011] (5) The prediction results of the AO-MTS_Mixers model, AO-Crossformer model and AO-DeepESN model are reconstructed to obtain the prediction results of the charging load.

[0012] Furthermore, in step (1), the ICEEMDAN model is used to decompose the charging load time series data to obtain several components, including:

[0013] Assume that the original charging load time series data signal is X = {x1, x2, ... x I}, where I represents the signal length; add a set of independent Gaussian white noise E1(ω (ie) ), the formula is as follows:

[0014]

[0015] Among them, ω (ie) Indicates the i-th e White noise, i e=1,2,…,I; E1(·) represents the first intrinsic mode function IMF1 component of a signal ICEEMDAN decomposition; β0 represents the weight coefficient when adding the noise component;

[0016] For each signal after adding noise Introduce the local mean operator M(·) to calculate its corresponding k-th order IMF k Components, by subtracting the local mean from the original signal, we get the k-th order IMF k Serving size:

[0017] Repeat the decomposition steps until the total number of components K is reached, k = 2, 3, ..., K, and find the kth IMF k Serving size:

[0018]

[0019] Where k = 2, 3, ..., K; r k represents the k-order residual; r k-1 represents the k-1 order residual; c k represents the kth IMF k β k-1 Represents the k-1th order weight coefficient.

[0020] Furthermore, the step (2) comprises:

[0021] Let each IMF component be represented by IMF n(t) , n∈1,2,…N, each IMF component IMF n(t) Partition into contiguous subsequences of m dimensions:

[0022] IMF n(t) ={x1,x2,…,x N},

[0023]

[0024] in, For the i IMF subsequences, N represents IMF n(t) Length;

[0025] Defining a subsequence and subsequence The maximum distance between corresponding elements is Calculate the Euclidean distance between each subsequence and other subsequences. The specific formula is as follows:

[0026]

[0027] in, Represents a subsequence Middle Elements and subsequences Middle The absolute value of the difference between the elements;

[0028] Calculate the i IMF The similarity between subsequences Given a threshold statistics The number of vectors is compared with the total number of vectors Nm, and the formula is as follows:

[0029]

[0030] in, express the number of

[0031] Calculate the sample entropy of each subsequence, The results are averaged and logarithmically transformed to obtain the sample entropy The formula is as follows:

[0032]

[0033] in, Indicates the i IMF The similarity of subsequences; Represents the average similarity of m-dimensional vectors; represents the sample entropy of the m-dimensional i-th subsequence;

[0034] Increase the dimension, add 1 to the dimension m, and get the similarity for m+1 vectors And the sample entropy

[0035] When the sample entropy change value fluctuates in the range of [0.1, 0.25], the iteration is stopped and the sample entropies SE1, SE2, …, SE of the K IMF components are calculated. K Composed a sample entropy vector S is , S is ={SE1,SE2,…,SE K};

[0036] From S is Randomly select three objects o1, o2 and o3, representing the initial category centers of high frequency, medium frequency and low frequency respectively. and Assign data objects to clusters, and calculate the distance between the remaining K-3 objects and the center point of each initial category. The calculation formula is as follows:

[0037]

[0038] Among them, o is represents the is-th sample entropy vector; represents the js-th center point, where js ∈ [1, 3]; represents the dimension of the data object; and each object is assigned to the cluster C represented by the nearest initial class center js , js = 1, 2, 3;

[0039] Introduce the formula J to measure the clustering effect:

[0040]

[0041] Among them, represents the is and spatial Euclidean distance; C js represents the js-th cluster; K represents the total number of sample entropies;

[0042] For each cluster C js , consider all non-center points within the cluster Calculate the total distance from all points within the cluster to the new center point when one of the non-center points is replaced by the center point, and denote it as J';

[0043] If J' < J, the total distance decreases after replacement, then make the replacement Otherwise, keep the original center point unchanged;

[0044] Repeat the steps of assigning data objects to clusters and updating cluster center points until the center points no longer change or reach the maximum number of iterations;

[0045] Finally, output 3 center points and their corresponding clusters C js , js = 1, 2, 3, representing the high-frequency, medium-frequency, and low-frequency components F1, F2, and F3 respectively.

[0046] Furthermore, the step (3) includes:

[0047] (31) Load the original high-frequency component F1, input the number of frequency band segments H and the number of decomposition modes ψ; Initialize the FIR filter bank with Z filters and start iteration, with the initial iteration ie = 1;

[0048] Evenly divide the frequency band of the original signal into Z segments, and the specific formulas for its upper and lower cut-off frequencies are as follows:

[0049]

[0050] Among them, f l represents the lowest cut-off frequency; f uIndicates the highest cutoff frequency; f s Indicates the sampling frequency of the original high-frequency signal;

[0051] (32) Passing Iterate to obtain the filtered modal component signal; where F1 represents the original high-frequency signal; represents the zth modal component signal of the ieth iteration; represents the response of the zth filter at the ieth iteration, z = 1, 2, ..., Z-1; * represents the convolution operation;

[0052] (33) Using the original high-frequency signal F1, decompose the modal component signal and the estimated period Update the filter coefficients; where, The autocorrelation spectrum reaches a local maximum after crossing the zero point The time is estimated; after completing one iteration, set ie=ie+1;

[0053] Determine whether the number of iterations has reached the preset number, if not, continue to update the filter and estimate the cycle;

[0054] (34) Calculate the correlation coefficient between each two modes and construct a correlation coefficient matrix CC (p×q) ; Select from the two modal components with the largest correlation coefficients and eliminate the mode with relatively small correlation kurtosis. At this time, Z-1 modes remain;

[0055] (35) discard the filter corresponding to the eliminated mode, and iteratively update the correlation coefficients of the filters corresponding to other modes again; repeat steps (33) to (34) until the number of remaining modal components is consistent with the preset ψ value, and when (Z-1) = ψ, stop the iteration;

[0056] If the number of remaining modal components is inconsistent with the preset ψ value, the process returns to continue = initialization.

[0057] Furthermore, the step (4) of using the improved artemisinin AO optimization algorithm to optimize the time step μ, the number of channels v and the prediction step ρ in the MTS-Mixers model includes:

[0058] A group of candidate solutions are generated by random initialization, each individual solution represents a parameter to be optimized in the AO-MTS_Mixers model, including: time step μ, number of channels ν, and prediction step ρ;

[0059] Assume that the entire population is initialized, denoted as A, the AO population size is Nb, the maximum number of iterations is Maxt, and the variable dimension is D. The formula is as follows:

[0060]

[0061] Among them, Hb and Lb represent the upper and lower bounds of the solution space range; Rb represents a set of random number sequences, and the value range is [0, 1]. The initial solution is generated using the random number sequence.

[0062] The Tent mapping is introduced to generate a chaotic sequence, and the population is initialized, updated, and optimized to make the initial solution as evenly distributed as possible within the solution space. The formula is as follows:

[0063]

[0064] Among them, Rb′ represents the chaotic sequence introduced by the Tent mapping; τ represents the control parameter, and τ ∈ rand(0, 1);

[0065] The population is updated to:

[0066] A Nb,D = Lb + Rb′ × (Hb - Lb)

[0067] The fit(·) function is introduced to represent the fitness value of an individual, and the individual with the best fitness value is selected as the optimal solution in the initial population.

[0068] The AO position is updated, and the comprehensive elimination stage is entered. The probability factor Q and the concentration factor θ are introduced as the guides for exploring the global solution space. When θ1 < Q, the position update formula for an individual is as follows:

[0069]

[0070] Among them, represents the position of the ib-th individual in the jb-th dimension at the (iter + 1)-th iteration; represents the position of the ib-th individual in the jb-th dimension at the iter-th iteration; θ represents the concentration factor; iter represents the current iteration number; represents the position of the optimal individual in the jb-th dimension at the iter-th iteration; θ1 represents a random number with a range of [0, 1];

[0071] Among them, the formula for the concentration factor θ is as follows:

[0072]

[0073] The formula for the probability coefficient Q is as follows:

[0074]

[0075] Among them, iter represents the current iteration number; Maxt represents the maximum iteration number;

[0076] The piecewise nonlinear convergence factor θ′ is introduced to slow down the factor attenuation rate based on the sine function and optimize the global search capability of the AO algorithm. The formula is as follows:

[0077]

[0078] AO performs local search and enters the local clearing stage; introduces the step size coefficient σ, continues to update the individual position for local search, and performs precise search in the searched solution space area to find the optimal solution. The formula is as follows:

[0079]

[0080] b1,b2,b3~U(1,N),b1≠b2≠b3

[0081] in, Indicates the position of the ibth individual in the iter+1th iteration; Represents an individual position randomly selected in the iterth iteration; and Indicates the positions of the other two individuals randomly selected in the iterth iteration, and b1, b2 and b3 are random and different indexes; Fit norm (ib) represents the normalized fitness value of the ith individual; fit(ib) represents the fitness value of the ibth individual; min(fit) and max(fit) represent the minimum and maximum fitness values; σ represents the step size coefficient;

[0082] The step coefficient σ of the AO algorithm is optimized by introducing the adaptive inertia weight. The formula is as follows:

[0083]

[0084] Among them, σ(iter) represents the step size coefficient in the iter-th iteration; b represents a parameter used to control the upper and lower limits of the weight; φ represents a parameter that controls the search smoothness;

[0085] Randomly choose to keep the current individual position or jump to the global optimal position to increase the algorithm's escape ability:

[0086]

[0087] Among them, best ib,jb represents the subvector of the current optimal solution of the jb-th dimension;

[0088] When the preset number of iterations is met, the iteration stops, otherwise the cycle continues until the optimal solution of time step μ, number of channels v and prediction step ρ is obtained.

[0089] Furthermore, in step (4), the AO-MTS_Mixers model is used to predict the component signals of the high-frequency component F1 in batches, including:

[0090] Assume that the input high-frequency component multivariate time series is χ h , The output sequence is

[0091] The input time series is standardized as follows:

[0092]

[0093] in, represents the data after standardization; χ h represents the input time series; norm(·) represents data normalization;

[0094] The obtained standardized data is encoded to obtain the processed input sequence. The formula is as follows:

[0095]

[0096] in, represents the input sequence after encoding operation; PE(·) represents sinusoidal position encoding;

[0097] The attention mechanism is used to learn the information features of the time dimension of the input data sequence to capture the characteristics of the time series. The formula is as follows:

[0098]

[0099] in, represents the temporal interaction feature sequence extracted by multi-head attention; Temporal(·) represents the extraction operation of information features in the time dimension; Attn(·) represents the feature extraction operation of the attention mechanism;

[0100] The extracted standardized data and the data after learning the time features are used to learn the information of the feature dimension. A feedforward neural network is used to learn the information of the feature dimension so that the model considers both the time series data and the time series features. The formula is as follows:

[0101]

[0102] in, Represents the sequence of extracted feature dimension information; Channel(·) represents the extraction operation of feature dimension information; FFN(·) represents the feedforward neural network operation;

[0103] A linear projection layer is used to perform information fusion and output mapping. The formula is as follows:

[0104]

[0105] in, represents the predicted output sequence of the ih-th high-frequency component, ih = 1, 2, 3, …, ψ; Linear(·) represents the output operation after information fusion and linear mapping of the two sequences.

[0106] Furthermore, in step (4), the data in the Crossformer model is divided into a length L by using the improved artemisinin AO optimization algorithm. seg , the number of layers l and the total dimension D are optimized to obtain the AO-Crossformer model, and the AO-Crossformer model is used to predict the intermediate frequency component F2, including:

[0107] Set the time series data set F2(t) of the intermediate frequency component F2 as the input data matrix The predicted output is Among them, T ime represents a past time point, τ time Indicates a point in time in the future. represents a matrix;

[0108] Use DSW technology to divide the data in each dimension into L seg The segmentation formula is as follows:

[0109]

[0110] Among them, X 1:T Represents the original historical data; X ts,d Represents the data at the tsth time in the d-dimensional original historical data; Indicates that the length in d dimension is L seg The icth data segment;

[0111] Then, linear projection is used to map the divided original data into the embedding space, and a time position code is given to each small segment:

[0112]

[0113] Among them, h i,d represents a time series segment of a data vector; represents a learnable projection matrix; Indicates the embedding position encoding of the icth data segment in the d dimension;

[0114] The linear projection of the segmented data and the position encoding of the segmented data are added to obtain the two-dimensional vector array H, and the formula is as follows:

[0115] H={h ic,d |1≤ic≤T ime / L seg ,1≤d≤D}

[0116] A two-stage attention mechanism is used to capture the cross-temporal and cross-dimensional dependencies of multivariate time series. The two-dimensional matrix Z is used to represent the input of the TSA layer, and MSA operations are performed on different dimensions:

[0117]

[0118] in, Represents the result of the MSA layer operation with dimension d; Z d Represents the data of all moments with dimension d; LayerNorm(·) represents the normalization layer operation;

[0119] Through multi-layer perceptron nonlinear transformation,

[0120]

[0121] Among them, Z time Represents the nonlinear transformation result of multilayer perceptron; Represents the result of multi-layer perceptron operation;

[0122] By Z time Capturing the cross-time dependencies of the same dimension, Z time As the input of the cross-dimensional stage, the dependencies between different dimensions are obtained;

[0123] Y=Z dim =TSA(Z)

[0124] Among them, Z dim represents the output of a multilayer perceptron; Z represents the input matrix of the TSA layer, and Y represents the output matrix of the TSA layer;

[0125] In each layer of the encoder except the first layer, two adjacent vectors in the time domain are merged to obtain a new layer, and TSA is used to capture the dependency of the layer; after obtaining Ni+1 output results of the encoder, the decoder uses the Ni+1 layer for prediction, and the formula is as follows:

[0126]

[0127] Where, l = 0,…,Ni; represents the icth data segment in the prediction result of the lth layer and dth dimension; represents a learnable matrix; represents the output of the decoder; Represents the prediction result of the lth layer; Indicates the final prediction result.

[0128] Furthermore, in step (4), the improved artemisinin AO optimization algorithm is used to optimize the storage layer leakage rate γ and the storage layer size P in the DeepESN model to obtain the AO-DeepESN model, and the low-frequency component F3 is predicted using the AO-DeepESN model, including:

[0129] Assume that DeepP-ESN consists of P storage layers and P-1 encoder layers to form a P-layer network; suppose there are N(iN) neurons in the iNth storage layer and M(jN) neurons in the jNth encoder layer, where 1≤iN≤P, 1≤jN≤P-1; the time length is set to T, the input signal matrix of the data is represented by X=[x(1),x(2),…,x(T)], the training signal matrix is ​​represented by d=[d(1),d(2),…,d(T)], and in time step t, x(t)∈R D and d(t)∈R L , the specific formula is as follows:

[0130]

[0131] Among them, z (iN) (t) represents the internal state of the iNth storage layer at time t; f(·) represents the nonlinear activation function tanh in the storage layer; represents the input of the iNth storage layer at time t+1; represents the activation value of the iNth storage layer at time t; represents the activation value of the iNth storage layer at time t+1, 1≤t≤T; W in(iN) represents the input weight of the iNth storage layer; W res(iN) represents the recurrent weight of the iNth storage layer;

[0132] The input of the iNth storage layer is the output of the iN-1th encoder layer. The function Represents a composite function, including high-dimensional projection and updating the state of the storage layer. The formula is updated as follows:

[0133]

[0134] Given the state of the previous pool, an unsupervised dimensionality reduction technique is used Encode and generate encoding features; the encoding process of the jNth encoder is expressed as:

[0135]

[0136] in, represents the activation value of the jNth encoder layer at time t, where 1≤t≤T; Represents a dimensionality reduction function, which converts the high-dimensional state vector Mapping to low-dimensional space; f enc (·) represents the activation function of the encoder; W enc(jN) represents the input weight of the jNth encoder layer;

[0137] According to the status of the previous storage layer, the status of the last storage layer of layer P is obtained. The specific formula is as follows:

[0138]

[0139] in, represents the activation value of the Pth storage layer at time t+1; symbol Represents a combination operator; Represents the high-dimensional projection of the jNth storage layer and the state composite function of the updated storage layer; Represents the encoding process of the jNth encoding layer; represents the coding function used by the jNth coding layer; u(t+1) represents the input signal at time t+1;

[0140] The output y(t+1) at time t+1 is as follows:

[0141]

[0142] Among them, W out represents the output weight; Represents the input feature matrix at time t+1;

[0143] Through the matrix W out With the input feature matrix The connection formula is as follows:

[0144]

[0145] in, Represents the final output matrix.

[0146] Furthermore, the step (5) comprises:

[0147] Prediction results for the AO-MTS-Mixers model Prediction results of AO-Crossformer Prediction results of AO-DeepESN model The weighted sum composite prediction model is fused, and the formula is:

[0148]

[0149] Where Y represents the predicted value of the charging load by the composite prediction model; α1, α2 and α3 represent random fusion coefficients in the range [0,1].

[0150] Based on the same inventive concept, the present invention also provides an electric vehicle charging load prediction system, comprising:

[0151] An initialization module is used to collect and obtain the historical charging load time series data of electric vehicles, and use the ICEEMDAN model to decompose the charging load time series data to obtain several components;

[0152] A classification module is used to calculate sample entropy of the components, and classify the signals using K-medoids clustering according to the sample entropy results to obtain a high-frequency component F1, a medium-frequency component F2, and a low-frequency component F3;

[0153] A secondary decomposition module is used to perform secondary decomposition on the high-frequency component F1 using a characteristic mode decomposition FMD model to obtain several components;

[0154] The preliminary prediction module is used to optimize the time step μ, number of channels ν and prediction step ρ in the MTS-Mixers model using the improved artemisinin AO optimization algorithm to obtain the AO-MTS_Mixers model, and use the AO-MTS_Mixers model to predict the component signals of the high-frequency component F1 in batches; use the improved artemisinin AO optimization algorithm to divide the data into the length L in the Crossformer model seg , the number of layers l and the total dimension D are optimized to obtain the AO-Crossformer model, and the AO-Crossformer model is used to predict the intermediate frequency component F2; the improved artemisinin AO optimization algorithm is used to optimize the storage layer leakage rate γ and the storage layer size P in the DeepESN model to obtain the AO-DeepESN model, and the AO-DeepESN model is used to predict the low frequency component F3;

[0155] The fusion prediction module is used to reconstruct the prediction results of the AO-MTS_Mixers model, the AO-Crossformer model and the AO-DeepESN model to obtain the prediction results of the charging load.

[0156] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages: 1. The present invention uses ICEEMDAN and FMD model algorithms to perform multi-level decomposition on the charging load data of electric vehicles, so as to more accurately extract the key signal components reflecting the charging load of electric vehicles; this refined signal processing method helps to capture the slight changes in the charging of electric vehicles in different use stages, thereby improving the accuracy of charging load prediction; 2. By combining K-medoids clustering and sample entropy calculation, the adaptability to signals of different frequencies is enhanced; the present invention can effectively identify and separate high-frequency, medium-frequency and low-frequency signal components; this frequency division processing method enables the prediction model to process signals of different characteristics more specifically, thereby improving the adaptability and overall prediction performance of the prediction model; 3. The selection and application of the prediction model are optimized. The present invention adopts three prediction models, AO-MTS-Mixers, AO-Crossformer and AO-DeepESN, which are respectively suitable for processing signal components of different frequencies; this sub-model prediction method can make full use of the prediction advantages of different models and improve the efficiency and accuracy of the prediction. BRIEF DESCRIPTION OF THE DRAWINGS

[0157] Figure 1 is a flow chart of a method according to an embodiment of the present invention;

[0158] Figure 2 It is a flow chart of the AO-MTS-Mixers prediction model of an embodiment of the present invention;

[0159] Figure 3 It is a flow chart of the AO-Crossformer prediction model of an embodiment of the present invention;

[0160] Figure 4 This is a flow chart of the AO-DeepESN prediction model of an embodiment of the present invention. DETAILED DESCRIPTION

[0161] In order to enable those skilled in the art to better understand the solution of the present application, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application.

[0162] As attached Figure 1 As shown, the electric vehicle charging load prediction method of this embodiment includes:

[0163] (1) Collecting and obtaining historical charging load time series data of electric vehicles, and using the ICEEMDAN model to decompose the charging load time series data into several components;

[0164] (2) performing sample entropy calculation on the components, and performing signal classification using K-medoids clustering according to the sample entropy results to obtain a high-frequency component F1, a mid-frequency component F2, and a low-frequency component F3;

[0165] (3) The high-frequency component F1 is decomposed twice using the characteristic mode decomposition (FMD) model to obtain several components;

[0166] (4) The improved artemisinin AO optimization algorithm is used to optimize the time step μ, number of channels ν and prediction step ρ in the MTS-Mixers model to obtain the AO-MTS_Mixers model, and the AO-MTS_Mixers model is used to batch predict the component signals of the high-frequency component F1; the improved artemisinin AO optimization algorithm is used to optimize the data partition length L in the Crossformer model. seg , the number of layers l and the total dimension D are optimized to obtain the AO-Crossformer model, and the AO-Crossformer model is used to predict the intermediate frequency component F2; the improved artemisinin AO optimization algorithm is used to optimize the storage layer leakage rate γ and the storage layer size P in the DeepESN model to obtain the AO-DeepESN model, and the AO-DeepESN model is used to predict the low frequency component F3;

[0167] (5) The prediction results of the AO-MTS_Mixers model, AO-Crossformer model and AO-DeepESN model are reconstructed to obtain the prediction results of the charging load.

[0168] Specifically, in order to improve the accuracy of electric vehicle charging load prediction, this embodiment proposes an ICEEMDAN-FMD secondary decomposition signal algorithm. First, the ICEEMDAN model is used to perform a preliminary decomposition of the charging load time series data to obtain 10 different components; then, the sample entropy of these components is calculated, and based on the calculation results, the K-medoids clustering method is used to classify the signals to distinguish three signal categories: high frequency, medium frequency and low frequency; then, for the high-frequency signal components, the FMD model is further used for secondary decomposition to obtain more refined 4 component signals; then, in the prediction stage, the MTS-Mixers prediction model is used to predict the 4 high-frequency component signals in batches, the Crossformer prediction model is used to predict the medium-frequency signals, and the DeepESN prediction model is used for the low-frequency signal components, and the improved AO optimization algorithm is used to optimize the parameters of the MTS-Mixers model (time step μ, number of channels ν, prediction step ρ), the parameters of the Crossformer model (data partition length L seg, number of layers l, total dimension D) and parameters of the DeepESN model (storage layer leakage rate γ, storage layer size P) are optimized. Finally, the prediction results of the three models are integrated and reconstructed to obtain the final electric vehicle charging load prediction results.

[0169] In step (1), the preprocessed charging load time series data is used as the input signal, and the original charging load time series data signal is set to X = {x1, x2, ... x I}, where I represents the signal length; add a set of independent Gaussian white noise to the original signal X The specific implementation process is as follows:

[0170]

[0171] Among them, ω (ie) Indicates the i-th e White noise, i e =1,2,…,I; E1(·) represents the first intrinsic mode function IMF1 component of a signal ICEEMDAN decomposition; β0 represents the weight coefficient when adding the noise component.

[0172] Furthermore, for each signal after adding noise The local mean operator M(·) is introduced to calculate the corresponding first-order IMF1 component. The first-order IMF1 component is obtained by subtracting the local mean from the original signal:

[0173]

[0174] Where r1 represents the first-order residual; c1 represents the first IMF1 value; and M(·) represents the local mean operator.

[0175] Continue with similar processing for the second IMF2 component:

[0176]

[0177] Among them, r2 represents the second-order residual; c2 represents the second IMF2 value; β1 represents the first-order weight coefficient.

[0178] Repeat the above steps until the predetermined total decomposition order K is reached, k = 2, 3, K, K (K represents the number of IMFs, K = 10), and find the kth IMF k Serving size:

[0179]

[0180] Where k = 2, 3, ..., K; r k represents the k-order residual; r k-1 represents the k-1 order residual; ck represents the kth IMF k β k-1 Represents the k-1th order weight coefficient.

[0181] From the above process, after multiple iterations, the original signal X is decomposed into K IMF components, that is, 10 IMF i Components, the original signal is successfully decomposed into the sum of 10 different frequency components.

[0182] In step (2), for the 10 signal components decomposed by the ICEEMDAN algorithm, the sample entropy of each component IMF is calculated to measure the signal complexity, and K-medoids clustering is used for signal classification. The specific implementation process is as follows:

[0183] First determine the subsequence length, assuming that each IMF component is represented by IMF n(t) , n∈1,2,L 10, each IMF component IMF n(t) Divide into m-dimensional continuous subsequences,

[0184] IMF n(t) ={x1,x2,…,x N},

[0185]

[0186] in, For the i IMF subsequences, N represents IMF n(t) Length;

[0187] For each subsequence Calculate other subsequences (where i≠j) the Euclidean distance between Defining a sequence with sequence The maximum distance between corresponding elements is Calculate the Euclidean distance between each subsequence and other subsequences. The specific formula is as follows:

[0188]

[0189] in, Represents a subsequence and The maximum distance between corresponding elements; Represents a subsequence Middle Elements and subsequences Middle The absolute value of the difference between the elements.

[0190] Furthermore, calculate the i-th IMF The similarity between subsequences Given a threshold statistics The number of vectors is compared with the total number of vectors Nm. The specific formula is as follows:

[0191]

[0192] in, express the number of

[0193] Further, the sample entropy of each subsequence is calculated. The results are averaged and logarithmically transformed to obtain the sample entropy The specific formula is as follows:

[0194]

[0195] in, Indicates the i IMF The similarity of subsequences; Represents the average similarity of m-dimensional vectors; represents the sample entropy of the m-dimensional i-th subsequence;

[0196] Further, we increase the dimension and add 1 to the dimension m, that is, we get the similarity for m+1 vectors. And the sample entropy

[0197] When m reaches a certain preset value or the sample entropy change value fluctuates in the range of [0.1, 0.25], the iteration stops. Finally, the sample entropy SE1, SE2, L, SE 10 Composed a sample entropy vector S is , S is ={SE1,SE2,…,SE 10}.

[0198] According to the sample entropy vector S calculated above is ={SE1,SE2,…,SE 10}, randomly select 3 objects o1, o2, o3, representing the initial category centers of high frequency, medium frequency and low frequency respectively Assign data objects to clusters, calculate the distance between the remaining 7 objects and the center point of each initial category, and assign each object to the cluster C represented by the nearest initial category center. js ,js=1,2,3.

[0199] The distance calculation formula is as follows:

[0200]

[0201] Among them, is Represents the is-th sample entropy vector; represents the jsth center point, js∈[1,3]; Represents the dimension of the data object;

[0202] Furthermore, formula J is introduced to measure the clustering effect:

[0203]

[0204] in, Indicates o is and The spatial Euclidean distance of C js represents the jsth cluster; K represents the total number of sample entropy.

[0205] Furthermore, for each cluster C js , considering all non-central points in the cluster Calculate one of the non-central points When replaced by the center point, the total distance from all points in the cluster to the new center point is recorded as J'.

[0206] If the total distance decreases after replacement, that is, J'<J, then replace Otherwise, keep the original center point unchanged.

[0207] Repeat the steps of assigning data objects to clusters and updating cluster centers until the centers no longer change or the maximum number of iterations is reached.

[0208] Finally, 3 center points are output and its corresponding cluster C js ,js=1,2,3, representing the three major components of high frequency, medium frequency and low frequency F1, F2, F3 respectively.

[0209] In step (3), the high-frequency signal F1 obtained by the K-medoids cluster analysis in step (2) is further decomposed by a secondary method using the FMD algorithm to further extract the detailed features of the high-frequency signal. The specific implementation process is as follows:

[0210] (31) Load the original high-frequency signal F1 and input the preset parameters, namely the number of frequency band divisions H and the number of decomposition modes ψ, ψ = 4. Initialize the FIR filter bank through Z filters and start iteration, with the initial iteration being ie = 1.

[0211] Furthermore, in order to realize the Hanning window initialization FIR filter bank, the frequency band of the original signal is evenly divided into Z segments, and the specific formulas of the upper and lower cutoff frequencies are as follows:

[0212]

[0213] Among them, f l Indicates the minimum cutoff frequency; f u Indicates the highest cutoff frequency; f s Indicates the sampling frequency of the original high-frequency signal.

[0214] (32) Further, filter update and period estimation are performed. Iterate to obtain the filtered modal component signal. Where F1 represents the original high-frequency signal; represents the zth modal component signal of the ieth iteration; represents the response of the zth filter at the ieth iteration, z = 1, 2, …, Z-1; “*” represents the convolution operation.

[0215] (33) Further, the original high-frequency signal F1 and the decomposed modal component signal are used and the estimated period Update the filter coefficients; where, The autocorrelation spectrum reaches a local maximum after crossing the zero point The moment is estimated; after completing one iteration, set ie=ie+1.

[0216] Determine whether the number of iterations reaches the pre-iteration number. If not, continue to update the filter and estimate the cycle; otherwise, continue with the following operations.

[0217] (34) Calculate the correlation coefficient between each two modes and construct a correlation coefficient matrix CC (p×q) . Select from the two modal components with the largest correlation coefficients, eliminate the mode with relatively small correlation kurtosis, and then Z-1 modes remain.

[0218]

[0219] Among them, CC pq represents the correlation coefficient between the two modes; u p (n) and u q (n) represent any two selected modes respectively; and Represent the means of the two modes respectively;

[0220] (35) Finally, discard the filter corresponding to the eliminated mode, and iterate again to update the correlation coefficients of the filters corresponding to other modes. Repeat steps (33) to (34) until the number of remaining modal components is consistent with the preset ψ value, that is, (Z-1) = ψ, and stop the iteration. If the condition is not met, return to continue initialization.

[0221] In step (4), according to the decomposition result of step (3),

[0222] (41) The AO-MTS_Mixers prediction model is used to predict the four component signals f1, f2, f3, and f4 obtained by decomposing the high-frequency component F1 in batches. The component signals f1, f2, f3, and f4 are respectively used as input sequences for prediction, and their data sequences are set as χ h .like Figure 2 The specific implementation process is as follows:

[0223] Step 1: Initialization of the improved artemisinin (AO) optimization algorithm. A group of candidate solutions are generated by random initialization. Each individual solution represents a parameter that needs to be optimized in the AO-MTS_Mixers model, including: time step μ, number of channels ν, and prediction step ρ. Assume that the entire population is initialized, denoted as A. The AO population size is Nb, the maximum number of iterations is Maxt, the variable dimension is D, and the specific formula is as follows:

[0224]

[0225] Among them, Hb and Lb represent the upper and lower bounds of the solution space; Rb represents a set of random number sequences whose value range is [0,1], and the random number sequence is used to generate the initial solution.

[0226] Furthermore, the Tent mapping is introduced to generate chaotic sequences, and the population is initialized and updated to optimize the initial solution so that the initial solution is distributed as evenly as possible in the solution space. The specific formula is as follows:

[0227]

[0228] Where Rb' represents the chaotic sequence introduced with Tent mapping; τ represents the control parameter, τ∈rand(0,1).

[0229] Furthermore, the population is updated as:

[0230] A Nb,D =Lb+Rb′×(Hb-Lb)

[0231] In addition, the fit(·) function is introduced to represent the fitness value of the individual, and the individual with the best fitness value is selected as the optimal solution in the initial population.

[0232] The second step is to update the AO position and enter the comprehensive elimination stage. The probability factor Q and the concentration factor θ are introduced as guides to explore the global solution space. When θ1<Q, the individual position update formula is as follows:

[0233]

[0234] in, Indicates the position of the ib-th individual in the jb-th dimension at the iter+1-th iteration; It represents the position of the ib-th individual in the jb-th dimension in the iter-th iteration; θ represents the concentration factor; iter represents the current number of iterations; represents the position of the optimal individual in the jbth dimension at the iterth iteration; θ1 represents a random number in the range [0,1];

[0235] Among them, the specific formula of concentration factor θ is as follows:

[0236]

[0237] The formula for the probability coefficient Q is as follows:

[0238]

[0239] Among them, iter represents the current number of iterations; Maxt represents the maximum number of iterations;

[0240] Furthermore, since the concentration factor θ′ exhibits nonlinear attenuation, the present invention introduces a piecewise nonlinear convergence factor, which slows down the factor attenuation rate based on a sine function and optimizes the global search capability of the AO algorithm. The specific formula is as follows:

[0241]

[0242] In the third step, AO performs local search and enters the local clearing stage. The step size coefficient σ is introduced to continue updating the position of the individual for local search, and perform precise search in the searched solution space area to find the optimal solution. The specific formula is as follows:

[0243]

[0244] b1,b2,b3~U(1,N),b1≠b2≠b3

[0245] in, Indicates the position of the ibth individual in the iter+1th iteration; Represents an individual position randomly selected in the iterth iteration; and Indicates the positions of the other two individuals randomly selected in the iterth iteration, and b1, b2 and b3 are random and different indexes; Fit norm (ib) represents the normalized fitness value of the ith individual; fit(ib) represents the fitness value of the ibth individual; min(fit) and max(fit) represent the minimum and maximum fitness values; σ represents the step size coefficient;

[0246] The step length coefficient σ is used to control the search step length of an individual in a local area. In order to perform a more detailed search in a local area, the present invention introduces an adaptive inertia weight to optimize the step length coefficient σ of the AO algorithm. The specific formula is as follows:

[0247]

[0248] Among them, σ(iter) represents the step size coefficient in the iter-th iteration; b represents a parameter used to control the upper and lower limits of the weight; φ represents a parameter that controls the search smoothness;

[0249] Step 4: Post-consolidation stage, in order to prevent the AO algorithm from falling into the local optimal solution too early, the algorithm's escape ability is increased by randomly selecting to maintain the current individual position or jump to the global optimal position.

[0250]

[0251] Among them, best ib,jb represents the subvector of the current optimal solution of the jb-th dimension;

[0252] Step 5: Determine whether the iterative optimization meets the preset number of iterations. If it does, the iteration stops. Otherwise, the cycle continues until the optimal solution is obtained, which is the three hyperparameters of the MTS-Mixers prediction model: time step μ, number of channels ν, and prediction step ρ.

[0253] Step 6: Substitute the three hyperparameters μ, ν, and ρ of the MTS-Mixers prediction model obtained by the improved AO optimization algorithm into the MTS-Mixers prediction model to predict the high-frequency components f1, f2, f3, and f4 in batches.

[0254] Assume that the input high-frequency component multivariate time series is χ h , the high-frequency component dataset is divided into training set and test set in a ratio of 7:3. Is a μ×ν historical time series matrix, where μ represents the sequence length and ν represents the number of channels in the time series. h Contains the historical data used to train the model. The goal of the prediction task is to predict all values ​​of the ν channels in the next μ time steps. represents a ρ×ν forecast time series data matrix, Where ρ represents the prediction step size. Contains the future time steps predicted by the model.

[0255] Step 7: Data standardization of the MTS-Mixers prediction model. First, the input time series is standardized. The overall prediction structure of MTS-Mixers can be summarized as learning the input χ h and output The model simplifies the convolution layer and attention filters into two modules for extracting features and a linear projection module to output the final measured result. The specific formula is as follows:

[0256]

[0257] in, represents the data after standardization; χ h represents the input time series; norm(·) represents data normalization.

[0258] Step 8: Information feature learning operation of the time dimension of the MTS-Mixers prediction model. The obtained standardized data is further encoded to obtain the processed input sequence. The specific formula is as follows:

[0259]

[0260] in, represents the input sequence after encoding operation; PE(·) represents sinusoidal position encoding.

[0261] After obtaining the standardized and encoded data, the attention mechanism is used to learn the information features of the time dimension of the input data sequence to capture the characteristics of the time series. The specific formula is as follows:

[0262]

[0263] in, Represents the temporal interaction feature sequence extracted by multi-head attention; Temporal(·) represents the extraction operation of information features in the time dimension; Attn(·) represents the feature extraction operation of the attention mechanism.

[0264] Step 9: Information learning operation of the feature dimension of the MTS-Mixers prediction model. Further, the information learning of the feature dimension is performed on the extracted and standardized data and the data after learning the time features, and the feedforward neural network is used to learn the information of the feature dimension so that the model can consider both time series data and time series features. The specific formula is as follows:

[0265]

[0266] in, Represents the sequence of extracted feature dimension information; Channel(·) represents the extraction operation of feature dimension information; FFN(·) represents the feedforward neural network operation.

[0267] Step 10: Linear projection operation of the MTS-Mixers prediction model. Finally, a linear projection layer performs information fusion and output mapping. The specific formula is as follows:

[0268]

[0269] in, represents the predicted output sequence of the ih-th high-frequency component, ih = 1, 2, 3, 4; Linear(·) represents the output operation after information fusion and linear mapping of the two sequences.

[0270] (42) The AO-Crossformer prediction model is used to predict the intermediate frequency component F2, such as Figure 3 As shown, the specific implementation process of the AO-Crossformer prediction model is:

[0271] Step 1: The AO optimization algorithm optimizes the parameters of the Crossformer prediction model to find the optimal value. These parameters include the data partition length L seg , number of layers l and total dimension D, and substitute them into the Crossformer model to predict the intermediate frequency component F2(t).

[0272] Step 2: Split the time series data set F2(t) of the intermediate frequency component F2 into a training set and a test set in a ratio of 7:3, and set it as the input data matrix The predicted output is X Time+1 , Among them, T ime represents a past time point, τ time represents a future time point, D represents the total dimension, Represents a matrix.

[0273] Step 3: Dimensional Segmentation Embedding Layer (DSW)

[0274] Use DSW technology to divide the data in each dimension into L seg The specific formula is as follows:

[0275]

[0276] Among them, X 1:T Represents the original historical data; X ts,d Represents the data at the tsth time in the d-dimensional original historical data; Indicates that the length in d dimension is L segThe icth data segment;

[0277] Then, linear projection is used to map the divided original data into the embedding space, and a time position code is given to each small segment:

[0278]

[0279] Among them, h i,d represents a time series segment of a data vector; represents a learnable projection matrix; Indicates the embedding position encoding of the icth data segment in the d dimension;

[0280] Furthermore, the linear projection of the segmented data is added to its position code to obtain a two-dimensional vector array H, the specific formula is as follows:

[0281] H={h ic,d |1≤ic≤T ime / L seg ,1≤d≤D}

[0282] Step 4: Two-Stage Attention (TSA)

[0283] A two-stage attention mechanism (TSA) is used to capture the cross-temporal and cross-dimensional dependencies of multivariate time series. The two-dimensional matrix Z is used to represent the input of the TSA layer, and MSA operations are performed on different dimensions:

[0284]

[0285] in, Represents the result of the MSA layer operation with dimension d; Z d Represents the data of all moments with dimension d; LayerNorm(·) represents the normalization layer operation;

[0286] Furthermore, a multilayer perceptron (MLP) is used to perform nonlinear transformation.

[0287]

[0288] Among them, Z time Represents the nonlinear transformation result of multilayer perceptron; Represents the result of multi-layer perceptron operation.

[0289] Furthermore, through Z time It captures the cross-time dependencies of the same dimension, and then Z time As the input of the cross-dimensional stage, it can obtain the dependencies between different dimensions.

[0290] Y=Z dim =TSA(Z)

[0291] Among them, Z dim represents the output of a multilayer perceptron; Z represents the input matrix of the TSA layer, and Y represents the output matrix of the TSA layer.

[0292] Step 5: Hierarchical Encoder-Decoder (HED)

[0293] Furthermore, in each layer of the encoder (except the first layer), two adjacent vectors in the time domain are merged to obtain a new layer, and TSA is used to capture the dependency of the layer. After obtaining Ni+1 output results of the encoder, the decoder uses the Ni+1 layer for prediction. The specific formula is as follows:

[0294]

[0295] Where, l = 0,…,Ni; represents the icth data segment in the prediction result of the lth layer and dth dimension; represents a learnable matrix; represents the output of the decoder; Represents the prediction result of the lth layer; Indicates the final prediction result.

[0296] (43) The AO-DeepESN prediction model is used to predict the low-frequency component F3, as shown in Figure 4 As shown, the specific implementation process of the AO-DeepESN prediction model is as follows:

[0297] Step 1: The AO optimization algorithm optimizes the parameters of the DeepESN prediction model and solves the optimal values. These parameters include the storage layer leakage rate γ and the storage layer size P, and are substituted into the DeepESN model to predict the low-frequency component F3.

[0298] Step 2: Split the low-frequency component F3 time series data set F3(t) into training and test sets in a ratio of 7:3, and set the DeeP-ESN to consist of P storage layers and P-1 encoder layers to form a P-layer network. Assume that there are N(iN) neurons in the iNth storage layer and M(jN) neurons in the jNth encoder layer, where 1≤iN≤P, 1≤jN≤P-1; the time length is set to T, the input signal matrix of the data is represented by X=[x(1),x(2),…,x(T)], the training signal matrix is ​​represented by d=[d(1),d(2),…,d(T)], and in time step t, x(t)∈R D and d(t)∈R L , the specific formula is as follows:

[0299]

[0300]

[0301] Among them, z (iN) (t) represents the internal state of the iNth storage layer at time t; f(·) represents the nonlinear activation function tanh in the storage layer; represents the input of the iNth storage layer at time t+1; represents the activation value of the iNth storage layer at time t; represents the activation value of the iNth storage layer at time t+1, 1≤t≤T; W in(iN) represents the input weight of the iNth storage layer; W res(iN) represents the cyclic weight of the iNth storage layer; γ represents the leakage rate of the storage layer.

[0302] The input of the iNth storage layer is the output of the iN-1th encoder layer. The function Represents a composite function, including high-dimensional projection and updating the state of the storage layer. The formula is updated as follows:

[0303]

[0304] in, Represents a composite function.

[0305] Furthermore, given the state of the previous pool, an unsupervised dimensionality reduction technique is used Encode and generate encoding features; the encoding process of the jNth encoder is expressed as:

[0306]

[0307] in, represents the activation value of the jNth encoder layer at time t, where 1≤t≤T; Represents a dimensionality reduction function, which converts the high-dimensional state vector Mapping to low-dimensional space; f enc (·) represents the activation function of the encoder; W enc(jN) represents the input weights of the jNth encoder layer.

[0308] Furthermore, according to the state of the previous storage layer, the state of the last P-th storage layer is obtained. The specific formula is as follows:

[0309]

[0310] in, represents the activation value of the Pth storage layer at time t+1; symbol Represents a combination operator; Represents the high-dimensional projection of the jNth storage layer and the state composite function of the updated storage layer; Represents the encoding process of the jNth encoding layer; represents the coding function used by the jNth coding layer; u(t+1) represents the input signal at time t+1;

[0311] Furthermore, the output y(t+1) at time t+1 is expressed as:

[0312]

[0313] Among them, W out represents the output weight; Represents the input feature matrix at time t+1;

[0314] Through the matrix W out With the input feature matrix The connection formula is as follows:

[0315]

[0316] in, Represents the final output matrix.

[0317] In step (5), the prediction result of the high-frequency component F1(t) obtained by the above AO-MTS-Mixers prediction model is The prediction results of the intermediate frequency component F2(t) obtained using the AO-Crossformer prediction model Using AO-DeepESN low frequency component The prediction results of the charging load are further integrated by using a weighted sum composite prediction model. The specific formula is as follows:

[0318]

[0319] Among them, Y represents the predicted value of the charging load by the composite prediction model; α1, α2 and α3 represent random fusion coefficients in the range [0,1], generally set α1∈[0.9,1], α2∈[0.6,0.9], α3∈[0,0.6].

[0320] This embodiment may also include step (6): In order to evaluate the training effect of the prediction model, error analysis is required. The following error parameter evaluation indicators are used: relative error RE, mean square relative error MRE, absolute error AE, and mean absolute error MAE are used to verify the performance of the model. The calculation formulas of the corresponding indicators are as follows:

[0321] Relative error RE:

[0322]

[0323] Mean square relative error MRE:

[0324]

[0325] Absolute error AE:

[0326]

[0327] Mean absolute error (MAE):

[0328]

[0329] Where N represents the number of sample data points. Represents the model's predicted value for the i-th data point, Y i Represents the actual value of the i-th data point.

[0330] Based on the same inventive concept, this embodiment also provides an electric vehicle charging load prediction system, including:

[0331] An initialization module is used to collect and obtain the historical charging load time series data of electric vehicles, and use the ICEEMDAN model to decompose the charging load time series data to obtain several components;

[0332] A classification module is used to calculate sample entropy of the components, and classify the signals using K-medoids clustering according to the sample entropy results to obtain a high-frequency component F1, a medium-frequency component F2, and a low-frequency component F3;

[0333] A secondary decomposition module is used to perform secondary decomposition on the high-frequency component F1 using a characteristic mode decomposition FMD model to obtain several components;

[0334] The preliminary prediction module is used to use the improved artemisinin AO optimization algorithm to optimize the time step μ, the number of channels v and the prediction step ρ in the MTS-Mixers model to obtain the AO-MTS_Mixers model, and use the AO-MTS_Mixers model to predict the component signals of the high-frequency component F1 in batches; use the improved artemisinin AO optimization algorithm to divide the data into the length L in the Crossformer modelseg , the number of layers l and the total dimension D are optimized to obtain the AO-Crossformer model, and the AO-Crossformer model is used to predict the intermediate frequency component F2; the improved artemisinin AO optimization algorithm is used to optimize the storage layer leakage rate γ and the storage layer size P in the DeepESN model to obtain the AO-DeepESN model, and the AO-DeepESN model is used to predict the low frequency component F3;

[0335] The fusion prediction module is used to reconstruct the prediction results of the AO-MTS_Mixers model, the AO-Crossformer model and the AO-DeepESN model to obtain the prediction results of the charging load.

Claims

1. A method for predicting electric vehicle charging load, characterized in that: include: (1) Collecting and obtaining historical charging load time series data of electric vehicles, and using the ICEEMDAN model to decompose the charging load time series data into several components; (2) performing sample entropy calculation on the components, and performing signal classification using K-medoids clustering according to the sample entropy results to obtain a high-frequency component F1, a mid-frequency component F2, and a low-frequency component F3; (3) The high-frequency component F1 is decomposed twice using the characteristic mode decomposition (FMD) model to obtain several components; (4) The improved artemisinin AO optimization algorithm is used to optimize the time step μ, number of channels ν and prediction step ρ in the MTS-Mixers model to obtain the AO-MTS_Mixers model, and the AO-MTS_Mixers model is used to batch predict the component signals of the high-frequency component F1; the improved artemisinin AO optimization algorithm is used to optimize the data partition length L in the Crossformer model. seg , the number of layers l and the total dimension D are optimized to obtain the AO-Crossformer model, and the AO-Crossformer model is used to predict the intermediate frequency component F2; the improved artemisinin AO optimization algorithm is used to optimize the storage layer leakage rate γ and the storage layer size P in the DeepESN model to obtain the AO-DeepESN model, and the AO-DeepESN model is used to predict the low frequency component F3; (5) The prediction results of the AO-MTS_Mixers model, AO-Crossformer model and AO-DeepESN model are reconstructed to obtain the prediction results of the charging load.

2. The electric vehicle charging load prediction method according to claim 1, characterized in that: In step (1), the ICEEMDAN model is used to decompose the charging load time series data to obtain several components, including: Assume that the original charging load time series data signal is X = {x1, x2, ... x I }, where I represents the signal length; add a set of independent Gaussian white noise to the original signal X The formula is as follows: in, Indicates the i-th e White noise, i e =1,2,…,I; E1(·) represents the first intrinsic mode function IMF1 component of a signal ICEEMDAN decomposition; β0 represents the weight coefficient when adding the noise component; For each signal after adding noise Introduce the local mean operator M(·) to calculate its corresponding k-th order IMF k Components, by subtracting the local mean from the original signal, we get the k-th order IMF k Serving size: Repeat the decomposition steps until the total number of components K is reached, k = 2, 3, ..., K, and find the kth IMF k Serving size: Where k = 2, 3, ..., K; r k represents the k-order residual; r k-1 represents the k-1 order residual; c k represents the kth IMF k β k-1 Represents the k-1th order weight coefficient.

3. The electric vehicle charging load prediction method according to claim 1, characterized in that: The step (2) comprises: Let each IMF component be represented by IMF n(t) , n∈1,2,…N, each IMF component IMF n(t) Partition into contiguous subsequences of m dimensions: IMF n(t) ={x1,x2,…,x N }, in, For the i IMF subsequences, N represents IMF n(t) Length; Defining a subsequence and subsequence The maximum distance between corresponding elements is Calculate the Euclidean distance between each subsequence and other subsequences. The specific formula is as follows: in, Represents a subsequence Middle Elements and subsequences Middle The absolute value of the difference between the elements; Calculate the i IMF The similarity between subsequences Given a threshold statistics The number of vectors is compared with the total number of vectors Nm, and the formula is as follows: in, express the number of Calculate the sample entropy of each subsequence, The results are averaged and logarithmically transformed to obtain the sample entropy The formula is as follows: in, Indicates the i IMF The similarity of subsequences; Represents the average similarity of m-dimensional vectors; represents the sample entropy of the m-dimensional i-th subsequence; Increase the dimension, add 1 to the dimension m, and get the similarity for m+1 vectors And the sample entropy When the sample entropy change value fluctuates in the range of [0.1, 0.25], the iteration is stopped and the sample entropies SE1, SE2, …, SE of the K IMF components are calculated. K Composed a sample entropy vector S is , S is ={SE1,SE2,…,SE K }; From S is Randomly select three objects o1, o2 and o3, representing the initial category centers of high frequency, medium frequency and low frequency respectively. and Assign data objects to clusters, and calculate the distance between the remaining K-3 objects and the center point of each initial category. The calculation formula is as follows: Among them, o is represents the is-th sample entropy vector; represents the js-th center point, where js ∈ [1, 3]; represents the dimension of the data object; and each object is assigned to the cluster C represented by the nearest initial class center js , where js = 1, 2, 3; Formula J is introduced to measure the clustering effect: in, Indicates o is and The spatial Euclidean distance of js represents the jsth cluster; K represents the total number of sample entropy; For each cluster C js , considering all non-central points in the cluster Calculate one of the non-central points When replaced by the center point, the total distance from all points in the cluster to the new center point is recorded as J′; If J′ < J and the total distance decreases after replacement, then perform replacement Otherwise, keep the original center point unchanged; Repeat the steps of assigning data objects to clusters and updating cluster centers until the center no longer changes or the maximum number of iterations is reached; Finally, 3 center points are output and its corresponding cluster C js ,js=1,2,3, representing the three major components of high frequency, medium frequency and low frequency F1, F2 and F3 respectively.

4. The electric vehicle charging load prediction method according to claim 1, characterized in that: The step (3) comprises: (31) Load the original high-frequency component F1, input the frequency band division number H and the decomposition mode number ψ; initialize the FIR filter bank through Z filters and start iteration, the initial iteration is ie = 1; The frequency band of the original signal is evenly divided into Z segments, and the specific formulas for the upper and lower cutoff frequencies are as follows: Among them, f l Indicates the minimum cutoff frequency; f u Indicates the highest cutoff frequency; f s Indicates the sampling frequency of the original high-frequency signal; (32) Passing Iterate to obtain the filtered modal component signal; where F1 represents the original high-frequency signal; represents the zth modal component signal of the ieth iteration; represents the response of the zth filter at the ieth iteration, z = 1, 2, ..., Z-1; * represents the convolution operation; (33) Using the original high-frequency signal F1, decompose the modal component signal and the estimated period Update the filter coefficients; where, The autocorrelation spectrum reaches a local maximum after crossing the zero point The time is estimated; after completing one iteration, set ie=ie+1; Determine whether the number of iterations has reached the preset number, if not, continue to update the filter and estimate the cycle; (34) Calculate the correlation coefficient between each two modes and construct a correlation coefficient matrix CC (p×q) ; Select from the two modal components with the largest correlation coefficients and eliminate the mode with relatively small correlation kurtosis. At this time, Z-1 modes remain; (35) discard the filter corresponding to the eliminated mode, and iteratively update the correlation coefficients of the filters corresponding to other modes again; repeat steps (33) to (34) until the number of remaining modal components is consistent with the preset ψ value, and when (Z-1) = ψ, stop the iteration; If the number of remaining modal components is inconsistent with the preset ψ value, return to continue initialization.

5. The electric vehicle charging load prediction method according to claim 1, characterized in that: The step (4) of using the improved artemisinin AO optimization algorithm to optimize the time step μ, the number of channels v and the prediction step ρ in the MTS-Mixers model includes: A group of candidate solutions are generated by random initialization, each individual solution represents a parameter to be optimized in the AO-MTS_Mixers model, including: time step μ, number of channels ν, and prediction step ρ; Assume that the entire population is initialized, denoted as A, the AO population size is Nb, the maximum number of iterations is Maxt, and the variable dimension is D. The formula is as follows: Among them, Hb and Lb represent the upper and lower bounds of the solution space; Rb represents a set of random number sequences with a value range of [0,1], and the initial solution is generated using the random number sequence; The Tent map is introduced to generate a chaotic sequence, and the population is initialized and updated to optimize so that the initial solution is distributed as evenly as possible in the solution space. The formula is as follows: Where Rb′ represents the chaotic sequence introduced with Tent mapping; τ represents the control parameter, τ∈rand(0,1); Update the population to: A Nb,D =Lb+Rb′×(Hb-Lb) The fit(·) function is introduced to represent the fitness value of an individual, and the individual with the best fitness value is selected as the optimal solution in the initial population. The AO position is updated, entering the comprehensive elimination stage. The probability factor Q and the concentration factor θ are introduced as the guides to explore the global solution space. When θ1 < Q, the individual position update formula is as follows: in, Indicates the position of the ib-th individual in the jb-th dimension at the iter+1-th iteration; It represents the position of the ib-th individual in the jb-th dimension in the iter-th iteration; θ represents the concentration factor; iter represents the current number of iterations; represents the position of the optimal individual in the jbth dimension at the iterth iteration; θ1 represents a random number in the range [0,1]; Among them, the formula for the concentration factor θ is as follows: The formula for the probability coefficient Q is as follows: Among them, iter represents the current iteration number; Maxt represents the maximum iteration number; The piecewise non-linear convergence factor θ′ is introduced to slow down the factor decay rate based on the sine function and optimize the global search ability of the AO algorithm. The formula is as follows: AO performs local search and enters the local elimination stage; the step size coefficient σ is introduced to continue updating the individual position for local search, and precise search is carried out within the searched solution space area to find the optimal solution. The formula is as follows: b1, b2, b3 ∼ U(1, N), b1 ≠ b2 ≠ b3 in, Indicates the position of the ibth individual in the iter+1th iteration; Represents an individual position randomly selected in the iterth iteration; and Indicates the positions of the other two individuals randomly selected in the iterth iteration, and b1, b2 and b3 are random and different indexes; Fit norm (ib) represents the normalized fitness value of the ith individual; fit(ib) represents the fitness value of the ibth individual; min(fit) and max(fit) represent the minimum and maximum fitness values; σ represents the step size coefficient; The adaptive inertia weight is introduced to optimize the step size coefficient σ of the AO algorithm. The formula is as follows: Among them, σ(iter) represents the step size coefficient in the iter-th iteration; b represents the parameter used to control the upper and lower limits of the weight; φ represents the parameter for controlling the search smoothness; Randomly select to keep the current individual position or jump to the global optimal position to increase the escape ability of the algorithm: Among them, best ib,jb represents the subvector of the current optimal solution of the jb-th dimension; When the preset iteration number is satisfied, the iteration stops; otherwise, continue the loop until the optimal solutions of the time step μ, the number of channels v, and the prediction step ρ are obtained.

6. The electric vehicle charging load prediction method according to claim 1, characterized in that: The batch prediction of each component signal of the high-frequency component F1 using the AO-MTS_Mixers model described in step (4) includes: Assume that the input high-frequency component multivariate time series is χ h , The output sequence is The input time series is standardized. The formula is as follows: in, represents the data after standardization; χ h represents the input time series; norm(·) represents data normalization; The obtained standardized data is encoded to obtain the processed input sequence. The formula is as follows: in, represents the input sequence after encoding operation; PE(·) represents sinusoidal position encoding; The attention mechanism is adopted to learn the information features in the time dimension of the input data sequence and capture the features of the time series. The formula is as follows: in, represents the temporal interaction feature sequence extracted by multi-head attention; Temporal(·) represents the extraction operation of information features in the time dimension; Attn(·) represents the feature extraction operation of the attention mechanism; The information learning in the feature dimension is performed on the standardized data and the data after learning the time features. The feedforward neural network is used for the information learning in the feature dimension, enabling the model to consider both the time series data and the time series features. The formula is as follows: in, Represents the sequence of extracted feature dimension information; Channel(·) represents the extraction operation of feature dimension information; FFN(·) represents the feedforward neural network operation; An information fusion and output mapping are performed by a linear projection layer. The formula is as follows: in, represents the predicted output sequence of the ih-th high-frequency component, ih = 1, 2, 3, …, ψ; Linear(·) represents the output operation after information fusion and linear mapping of the two sequences.

7. The electric vehicle charging load prediction method according to claim 1, characterized in that: Step (4) uses the improved artemisinin AO optimization algorithm to divide the data in the Crossformer model into a length L seg , the number of layers l and the total dimension D are optimized, and the intermediate frequency component F2 is predicted using the optimized Crossformer prediction model, including: Set the time series data set F2(t) of the intermediate frequency component F2 as the input data matrix The predicted output is Among them, T ime represents a past time point, τ time Indicates a point in time in the future. represents a matrix; Use DSW technology to divide the data in each dimension into L seg The segmentation formula is as follows: Among them, X 1:T Represents the original historical data; X ts,d Represents the data at the tsth time in the d-dimensional original historical data; Indicates that the length in d dimension is L seg The icth data segment; Then, the original data after division is mapped to the embedding space by linear projection, and a time position encoding is given to each small segment: Among them, h i,d represents a time series segment of a data vector; represents a learnable projection matrix; Indicates the embedding position encoding of the icth data segment in the d dimension; The linear projection of the segmented data is added to the position encoding of the segmented data to obtain a two-dimensional vector array H. The formula is as follows: H={h ic,d |1≤ic≤T ime / L seg ,1≤d≤D} The two-stage attention mechanism is adopted to capture the cross-time and cross-dimension dependencies of the multivariate time series. The input of the TSA layer is represented by a two-dimensional matrix Z, and the MSA operation is performed on different dimensions respectively: in, Represents the result of the MSA layer operation with dimension d; Z d Represents the data of all moments with dimension d; LayerNorm(·) represents the normalization layer operation; Non-linear transformation is performed through a multi-layer perceptron, Among them, Z time Represents the nonlinear transformation result of multilayer perceptron; Represents the result of multi-layer perceptron operation; By Z time Capturing the cross-time dependencies of the same dimension, Z time As the input of the cross-dimensional stage, the dependencies between different dimensions are obtained; Y=Z dim =TSA(Z) Among them, Z dim represents the output of a multilayer perceptron; Z represents the input matrix of the TSA layer, and Y represents the output matrix of the TSA layer; In each layer of the encoder except the first layer, the two adjacent vectors in the time domain are combined to obtain a new layer, and the dependencies of this layer are captured using TSA; after obtaining the Ni+1 output results of the encoder, the decoder uses Ni+1 layers for prediction. The formula is as follows: in, represents the icth data segment in the prediction result of the lth layer and dth dimension; represents a learnable matrix; represents the output of the decoder; Represents the prediction result of the lth layer; Indicates the final prediction result.

8. The electric vehicle charging load prediction method according to claim 1, characterized in that: Step (4) uses the improved artemisinin AO optimization algorithm to optimize the storage layer leakage rate γ and the storage layer size P in the DeepESN model, and uses the optimized DeepESN prediction model to predict the low-frequency component F3, including: Assume that DeepP-ESN consists of P storage layers and P-1 encoder layers to form a P-layer network; suppose there are N(iN) neurons in the iNth storage layer and M(jN) neurons in the jNth encoder layer, where 1≤iN≤P, 1≤jN≤P-1; the time length is set to T, the input signal matrix of the data is represented by X=[x(1),x(2),…,x(T)], the training signal matrix is ​​represented by d=[d(1),d(2),…,d(T)], and in time step t, x(t)∈R D and d(t)∈R L , the specific formula is as follows: Among them, z (iN) (t) represents the internal state of the iNth storage layer at time t; f(·) represents the nonlinear activation function tanh in the storage layer; represents the input of the iNth storage layer at time t+1; represents the activation value of the iNth storage layer at time t; represents the activation value of the iNth storage layer at time t+1, 1≤t≤T; W in(iN) represents the input weight of the iNth storage layer; W res(iN) represents the cyclic weight of the iNth storage layer; The input of the iNth storage layer is the output of the iN-1th encoder layer. The function Represents a composite function, including high-dimensional projection and updating the state of the storage layer. The formula is updated as follows: Given the state of the previous pool, an unsupervised dimensionality reduction technique is used Encode and generate encoding features; the encoding process of the jNth encoder is expressed as: in, represents the activation value of the jNth encoder layer at time t, where 1≤t≤T; Represents a dimensionality reduction function, which converts the high-dimensional state vector Mapping to low-dimensional space; f enc (·) represents the activation function of the encoder; W enc(jN) represents the input weight of the jNth encoder layer; According to the status of the previous storage layer, the status of the last storage layer of layer P is obtained. The specific formula is as follows: in, represents the activation value of the Pth storage layer at time t+1; symbol Represents a combination operator; Represents the high-dimensional projection of the jNth storage layer and the state composite function of the updated storage layer; Represents the encoding process of the jNth encoding layer; represents the coding function used by the jNth coding layer; u(t+1) represents the input signal at time t+1; The output y(t+1) at time t+1 is as follows: Among them, W out represents the output weight; Represents the input feature matrix at time t+1; Through the matrix W out With the input feature matrix The connection formula is as follows: in, Represents the final output matrix.

9. The electric vehicle charging load prediction method according to claim 1, characterized in that: The step (5) comprises: Prediction results for the AO-MTS-Mixers model Prediction results of AO-Crossformer Prediction results of AO-DeepESN model The weighted sum composite prediction model is fused, and the formula is: Where Y represents the predicted value of the charging load by the composite prediction model; α1, α2 and α3 represent random fusion coefficients in the range [0,1].

10. An electric vehicle charging load prediction system, characterized in that: include: An initialization module is used to collect and obtain the historical charging load time series data of electric vehicles, and use the ICEEMDAN model to decompose the charging load time series data to obtain several components; A classification module is used to calculate sample entropy of the components, and classify the signals using K-medoids clustering according to the sample entropy results to obtain a high-frequency component F1, a medium-frequency component F2, and a low-frequency component F3; A secondary decomposition module is used to perform secondary decomposition on the high-frequency component F1 using a characteristic mode decomposition FMD model to obtain several components; The preliminary prediction module is used to optimize the time step μ, number of channels ν and prediction step ρ in the MTS-Mixers model using the improved artemisinin AO optimization algorithm to obtain the AO-MTS_Mixers model, and use the AO-MTS_Mixers model to predict the component signals of the high-frequency component F1 in batches; use the improved artemisinin AO optimization algorithm to divide the data into the length L in the Crossformer model seg , the number of layers l and the total dimension D are optimized to obtain the AO-Crossformer model, and the AO-Crossformer model is used to predict the intermediate frequency component F2; the improved artemisinin AO optimization algorithm is used to optimize the storage layer leakage rate γ and the storage layer size P in the DeepESN model to obtain the AO-DeepESN model, and the AO-DeepESN model is used to predict the low frequency component F3; The fusion prediction module is used to reconstruct the prediction results of the AO-MTS_Mixers model, the AO-Crossformer model and the AO-DeepESN model to obtain the prediction results of the charging load.

Images