Method applied to electric vehicle charging load prediction

Abstract

The invention discloses a method for predicting the charging load of an electric vehicle, and the method comprises the steps: collecting the historical charging load data of the electric vehicle, and carrying out the one-time decomposition of the time series data of the charging load through employing an ICEEMDAN model, and obtaining 10 components; performing sample entropy calculation on components obtained by primary decomposition, performing signal classification by adopting K-medoids clustering according to a sample entropy result to obtain a high-frequency component, an intermediate-frequency component and a low-frequency component respectively, optimizing parameters in an MTS-Mixers model, a Crossform model and a DeepESN model by using an improved AO optimization algorithm, and performing signal classification by adopting K-medoids clustering to obtain a high-frequency component, an intermediate-frequency component and a low-frequency component; and predicting the high-frequency component, the intermediate-frequency component and the low-frequency component by using an MTS-Mixers model, a Crossform model and a DeepESN model, and finally reconstructing prediction results of the three models to obtain a prediction result of the charging load.

Description

Technical Field

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

[0002] With the booming electric vehicle industry and increasingly prominent environmental issues, electric vehicles are becoming an important option for future transportation. However, the widespread adoption of electric vehicles also brings new challenges, one of which is charging load forecasting. Accurate charging load forecasting is crucial for grid scheduling, charging facility planning, and electric vehicle charging management. It helps grid companies rationally allocate power resources and avoid overload or underload on the power system. It also provides charging facility operators with a scientific basis for planning and construction, ensuring that electric vehicle users have convenient and efficient access to charging services.

[0003] Currently, commonly used charging load forecasting methods include those based on statistical models, machine learning, and deep learning. Statistical model-based forecasting methods primarily use statistical analysis to derive patterns and trends in charging load changes. While simple and easy to implement, these methods suffer from low prediction accuracy and struggle to capture the complex variations in charging load. Machine learning-based forecasting methods primarily utilize machine learning algorithms to analyze and model charging load data. While these methods offer high prediction accuracy, they require extensive training data and computing resources, and their model generalization capabilities are limited. Deep learning-based forecasting methods primarily utilize deep learning algorithms to automatically learn the complex features of charging load data. These methods offer high prediction accuracy, but their models are complex, training is challenging, and they require extensive training data.

[0004] Electric vehicle charging load data is influenced by numerous factors, such as vehicle type, driving habits, charging facility distribution, and weather conditions. This results in complex, nonlinear characteristics, making direct prediction difficult. Furthermore, existing prediction models often overlook the characteristics of different frequency components in charging load data, resulting in insufficient prediction accuracy. To overcome these limitations, further research and development of new EV charging load prediction methods is needed, 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 for predicting the charging load of electric vehicles, which can improve the accuracy of charging load prediction and provide strong support for the charging use and management of electric vehicles.

[0006] Technical solution: The method for predicting charging load of electric vehicles described in the present invention is characterized by comprising the following steps:

[0007] S1: Collect historical charging load data of electric vehicles. The data characteristics of historical charging load data include the starting state of charge, the charging end state of charge, the charging time distribution, and the vehicle charging interval. Use the ICEEMDAN model to decompose the historical charging load data to obtain several IMF components;

[0008] S2: Calculate the sample entropy of several IMF components obtained by the primary decomposition, and use K-medoids clustering to classify the signal based on the sample entropy results to obtain three signal components: high-frequency component F1, medium-frequency component F2, and low-frequency component F3;

[0009] S3: The high-frequency component F1 is decomposed twice using the FMD eigenmode decomposition model to obtain four component signals f1, f2, f3 and f4;

[0010] S4: The MTS-Mixers prediction model is used to predict the four component signals f1, f2, f3 and f4 of the high-frequency component F1 in batches, the Crossformer prediction model is used to predict the mid-frequency component F2, and the DeepESN prediction model is used to predict the low-frequency component F3;

[0011] S5: Reconstruct the prediction results of the three models to obtain the prediction result of the charging load.

[0012] Preferably, the specific process of the ICEEMDAN model in S1 decomposing the historical charging load data to obtain a number of IMF components includes the following steps:

[0013] S1.1: Set the original charging load time series data signal to X = {x1, x2, ...x I}, where I represents the signal length, and a set of independent Gaussian white noise E1(ω (i) ), the specific formula is as follows:

[0014] X (i) =X+β0E1(ω (i) )

[0015] Among them, ω (i) represents the i-th white noise added, i = 1, 2, 3, ..., I, E(·) represents the first intrinsic mode function IMF1 component of the ICEEMDAN decomposition of a signal; β0 represents the weight coefficient when the noise component is added;

[0016] S1.2: For each signal X after adding noise (i) , introduce the local mean operator M(·) to calculate its corresponding k-th order IMF kComponent, the k-th order IMF is obtained by subtracting the local mean from the original signal k Quantity;

[0017] S1.3: Repeat the decomposition steps until the predetermined total decomposition order K is reached, k = 2, 3...K, where K represents the number of IMFs. Find the kth IMF k Quantity;

[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-1 order weight coefficient;

[0020] S1.4: Repeat the above process for multiple iterations. The original signal X is decomposed into K IMF components. That is, the original signal X is decomposed into the sum of several different frequency components.

[0021] Preferably, the specific implementation process of sample entropy calculation and K-medoids clustering in S2 includes the following steps:

[0022] S2.1: Assume that each IMF component is represented as IMF n(t) , n∈1,2,...,N, each IMF component IMF n(t) Divide into continuous subsequences of m dimensions:

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

[0024] x i ={x i ,x i+1 ,…,x i+m-1}, i=1~Nm-1

[0025] Where N represents IMF n(t) length;

[0026] S2.2: Define the sequence x i and the sequence x j The maximum distance between corresponding elements is d max (x i ,x j ), calculate the Euclidean distance between each subsequence and other subsequences. The specific formula is as follows:

[0027]

[0028] Among them, d max (x i ,x j ) represents the subsequence x i and x j The maximum distance between corresponding elements, Represents a subsequence x i Middle elements and subsequence x j Middle The absolute value of the difference between the elements;

[0029] S2.3: Calculate the similarity between subsequences statistics The number of vectors is compared with the total number of vectors Nm. The specific formula is as follows:

[0030]

[0031] in express the number of

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

[0033]

[0034] in, represents the similarity of the i-th subsequence; Represents the average similarity of m-dimensional vectors; represents the sample entropy of the m-dimensional i-th subsequence;

[0035] S2.5: Increase the dimension and then add 1 to the dimension m, that is, get the similarity for m+1 vectors And sample entropy

[0036] S2.6: When m reaches the preset value or the sample entropy stops changing, the iteration stops and the sample entropies SE1, SE2, ..., SE n Form a sample entropy vector S i , S i ={SE1,SE2,…,SE n};

[0037] S2.7: From the sample entropy vector S i Randomly select three objects o1, o2 and o3, representing the initial category centers of high frequency, medium frequency and low frequency respectively and Calculate the distance between the remaining objects and the center point of each initial category, and assign each object to the cluster C represented by the nearest initial category center. j ,j=1,2,3, the distance calculation formula is as follows:

[0038]

[0039] Among them i represents the i-th sample entropy vector; represents the jth center point, j∈[1,3]; Represents the dimension of the data object;

[0040] S2.8: Introduce formula J to measure the clustering effect:

[0041]

[0042] in, Indicates the oth i and The spatial Euclidean distance of C j represents the jth cluster; K represents the total number of sample entropy;

[0043] S2.9: For each cluster C j , considering all non-central points in the cluster Calculate one of the non-center points When the center point is replaced, the total distance from all points in the cluster to the new center point is recorded as J′;

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

[0045] S2.10: 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, and finally output 3 center points and its corresponding cluster C j ,j=1,2,3, representing the three categories of high-frequency, medium-frequency and low-frequency signals F1, F2 and F3 respectively.

[0046] Preferably, the S3 specifically includes the following steps:

[0047] S3.1: Load the original high-frequency signal F1 and enter the preset parameters: filter length L, number of frequency band divisions H, and number of hierarchical modes ψ, where ψ = 4. Initialize the FIR filter using Z filters and start iteration, with the initial iteration i = 1.

[0048] S3.2: Divide the frequency band of the original signal into Z segments evenly. The upper and lower cutoff frequency formulas are as follows:

[0049]

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

[0051] S3.3: Pass-through Iterate to obtain the filtered modal component signal, where represents the zth modal component signal of the i-th iteration; represents the response of the zth filter at the i-th iteration, z = 1, 2, ..., Z-1, and "*" represents the convolution operation;

[0052] S3.4: Use the original high-frequency signal F1 to 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 moment of estimation is obtained, one iteration is completed, and i=i+1 is set;

[0053] S3.5: Determine whether the number of iterations reaches the pre-iteration number. If not, continue to update the filter and estimate the cycle; otherwise, proceed to the next step.

[0054] S3.6: 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 coefficient, eliminate the mode with relatively small correlation kurtosis, and at this time there are Z-1 modes remaining;

[0055] S3.7: Discard the filter corresponding to the eliminated mode and iteratively update the correlation coefficients of the filters corresponding to other modes. Repeat the above steps until the number of remaining modal components is consistent with the preset ψ, that is, (Z-1) = ψ, and stop the iteration;

[0056] S3.8: If the conditions are not met, return to continue the Hanning window initialization FIR filter bank and use the remaining Z-1 modal components, that is, ψ modal components, as the final component signals of the characteristic mode decomposition, that is, 4 component signals f1, f2, f3 and f4.

[0057] As a preferred method, in order to improve the prediction performance of the prediction model in S4, the improved artemisinin AO optimization algorithm is used to adjust the time step μ, the number of channels υ and the prediction step ρ in the MTS-Mixers model, and the data partition length L in the Crossformer model. seg, the number of layers l and the total dimension D, the storage layer leakage rate γ and the storage layer size P in the DeepESN model are optimized respectively.

[0058] As a preference, the improved artemisinin AO optimization algorithm is used to optimize the three hyperparameters of the MTS-Mixers model: time step μ, number of channels υ and prediction step ρ, and the data partition length L in the Crossformer model. seg , the number of layers l and the total dimension D, as well as the specific steps for optimizing the storage layer leakage rate γ and storage layer size P in the DeepESN model are:

[0059] S4.11: Randomly initialize and generate a group of candidate solutions. 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, and the total dimension is D. The specific formula is:

[0060]

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

[0062] S4.12: Introduce the Tent map to generate chaotic sequences and perform initialization and update optimization on the population so that the initial solutions are distributed as evenly as possible in the solution space. The specific formula is as follows:

[0063]

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

[0065] S4.13: 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 the individual, and the individual with the best fitness value is selected as the optimal solution in the initial group.

[0068] AO position is updated and enters the comprehensive elimination stage, introducing probability factor Q and concentration factor As a guide to explore the global solution space, when When , the individual position update formula is as follows:

[0069]

[0070] in, represents the position of the i-th individual in the j-th dimension at the t+1-th iteration; Indicates the position of the i-th individual in the j-th dimension in the t-th iteration; t represents the current number of iterations; represents the position of the optimal individual in the jth dimension at the tth iteration; Represents a random number in the range [0,1];

[0071] Concentration factor The specific formula is:

[0072]

[0073] Where t represents the current number of iterations; Maxt is the maximum number of iterations;

[0074] S4.14: Introduce a piecewise nonlinear convergence factor, based on the sine function to slow down the factor attenuation rate and optimize the global search capability of the AO algorithm. The specific formula is as follows:

[0075]

[0076] S4.15: AO performs local search and enters the local clearance phase. It introduces the step size coefficient σ and continues to update the individual position to perform local search. It performs precise search in the searched solution space area to find the optimal solution. The specific formula is as follows:

[0077]

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

[0079] in, represents the position of the i-th individual in the t+1th iteration; represents the position of an individual randomly selected in the tth iteration; and Indicates the positions of the other two individuals randomly selected in the tth iteration, and b1, b2 and b3 are random and different indices; Fit norm (i) represents the normalized fitness value of the i-th individual; fit(i) represents the fitness value of the i-th individual; min(fit) and max(fit) represent the minimum and maximum fitness values, respectively; the step coefficient σ is a random value between [0.1, 0.6];

[0080] S4.16: Introduce adaptive inertia weight to optimize the step size coefficient σ of the AO algorithm. The specific formula is as follows:

[0081]

[0082] Where σ(t) represents the step size coefficient in the tth iteration; b represents the parameter that controls the upper and lower limits of the weight; φ represents the parameter that controls the search smoothness;

[0083] S4.17: To prevent the AO algorithm from falling into the local optimal solution too early, the algorithm's escape ability is increased by randomly choosing to maintain the current individual position or jump to the global optimal position. The specific formula is:

[0084]

[0085] Among them, best i,j Represents the subvector of the current optimal solution of the j-th dimension;

[0086] S4.18: Determine whether the iterative optimization meets the maximum iteration condition. If it does, the iteration stops. Otherwise, the loop continues until the optimal solution is obtained. These are the three hyperparameters of the MTS-Mixers prediction model: time step μ, number of channels υ, and prediction step ρ.

[0087] S4.19: Repeat S4.11-S4.18 to obtain the data partition length L in the Crossformer model seg , the number of layers l and the total dimension D, as well as the storage layer leakage rate γ and storage layer size P in the DeepESN model.

[0088] As a preference, the three hyperparameters of the MTS-Mixers prediction model obtained by the improved AO optimization algorithm, namely, the time step μ, the number of channels υ, and the prediction step ρ, are substituted into the MTS-Mixers prediction model to predict the four component signals f1, f2, f3, and f4 of the high-frequency component F1 in batches, specifically including the following steps:

[0089] S4.21: Set the input high-frequency component multivariate time series to χ h , Where μ represents the sequence length, v represents the number of channels of the time series, and the output sequence is Where ρ represents the prediction step size;

[0090] S4.22: Normalize the input time series. The specific formula is as follows:

[0091]

[0092] in, represents the standardized data; χ h represents the high-frequency component of the input multivariate time series; norm(·) represents data normalization;

[0093] S4.23: Perform encoding on the obtained normalized data to obtain the processed input sequence. The specific formula is as follows:

[0094]

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

[0096] S5.24: 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:

[0097]

[0098] 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;

[0099] S4.25: The extracted and 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 can consider both time series data and time series features. The specific formula is as follows:

[0100]

[0101] 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;

[0102] S4.26: Finally, a linear projection layer performs information fusion and output mapping. The specific formula is as follows:

[0103]

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

[0105] As a preference, the data partition length L in the Crossformer model obtained by the improved AO optimization algorithm is seg The specific implementation process of substituting the number of layers l and the total dimension D into the Crossformer prediction model to predict the intermediate frequency component F2 is as follows:

[0106] S4.31: Set the intermediate frequency component F2(t) time series data set as the input data matrix The predicted output is Where T ime represents a past time point, τ time represents a future time point, D represents the total dimension, represents a matrix;

[0107] S4.32: Use DSW technology to divide the data in each dimension into L seg , the specific formula is as follows:

[0108]

[0109] where X 1:T Indicates original historical data; X t,d Represents the data at time t in dimension d in the original historical data; Indicates that the length in d dimension is L seg The i-th data segment;

[0110] S4.33: Use linear projection to map the divided raw data into the embedding space, and give each small segment a time position code:

[0111]

[0112] where h i,d represents a time series segment of a data vector; represents the learnable projection matrix; Indicates the position encoding of the i-th data segment in the d-th dimension;

[0113] S4.34: The linear projection of the segmented data and its positional encoding are added to obtain the two-dimensional vector array H. The specific formula is as follows:

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

[0115] S4.35: Two-Stage Attention 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 the MSA operation is performed on different dimensions separately:

[0116]

[0117] in Represents the MSA layer operation result of dimension d; Z dRepresents the data of dimension d at all times; LayerNorm(·) represents the normalization layer operation;

[0118] S4.36: Perform nonlinear transformations using a multilayer perceptron. The specific formula is as follows:

[0119]

[0120] where Z time Represents the nonlinear change result of multilayer perceptron; Represents the result of multi-layer perceptron operation;

[0121] S4.37: Through Z time Capture the cross-time dependencies of the same dimension, and then Z time As input to the cross-dimensional stage to obtain the dependencies between different dimensions:

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

[0123] where Z dim represents the output of the multilayer perceptron; Z, Z represents the input matrix of the TSA layer, and Y represents the output matrix of the TSA layer;

[0124] S4.38: In each layer of the encoder except the first layer, two adjacent vectors in the time domain are merged to form a new layer, and the dependencies of this layer are captured using TSA. After obtaining the N+1 output results of the encoder, the decoder uses the N+1 layers (0, 1, ..., N) for prediction. The specific formula is as follows:

[0125]

[0126] Where l represents the number of layers, Represents the i-th data segment in the prediction result of the l-th layer d dimension; represents a learnable matrix; represents the output of the decoder; Represents the prediction result of the lth layer; Indicates the final prediction result.

[0127] As a preference, the specific implementation process of substituting the storage layer leakage rate γ and storage layer size P in the DeepESN model obtained by the improved AO optimization algorithm into the DeepESN model prediction model to predict the low-frequency component F3 is as follows:

[0128] S4.41: Assume that DeepPESN consists of P storage layers and P-1 encoder layers to form a P-layer network. Assume that there are N(i) neurons in the i-th storage layer and M(j) neurons in the j-th encoder layer, where 1≤i≤P, 1≤j≤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:

[0129]

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

[0131] S4.42: The input of the i-th storage layer is the output of the i-1-th encoder layer, introducing the function To express the conforming function, including high-dimensional projection and updating the state of the storage layer, the formula is updated as follows:

[0132]

[0133] in represents a composite function;

[0134] S4.43: Given the previous state of the reservoir, use unsupervised dimensionality reduction techniques Encode and generate encoding features. The encoding process of the j-th encoder is expressed as follows:

[0135]

[0136] in, represents the activation value of the j-th encoder layer at time t, 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(i)represents the input weight of the j-th encoder layer;

[0137] S4.44: Based on the states of the storage layers from the 1st layer to the P-1th layer, the state of the last storage layer at the Pth layer is obtained. The specific formula is as follows:

[0138]

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

[0140] S4.45: Output y(t+1) at time t+1, specifically:

[0141]

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

[0143] S4.46: Through the matrix W out With the input feature matrix The specific formula is as follows:

[0144]

[0145] in Represents the final output matrix.

[0146] As a preferred embodiment, the specific implementation process of S5 is:

[0147] The prediction results of the high-frequency component F1(t) obtained by using the AO-MTS-Mixers prediction model The prediction results of the intermediate frequency component F2(t) obtained using the AO-Crossformer prediction model And the prediction results of AO-DeepESN low-frequency component F3(t) The prediction result of charging load is obtained by using a weighted sum composite prediction model to combine the results. The specific formula is:

[0148]

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

[0150] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages:

[0151] (1) The present invention uses the ICEEMDAN and FMD model algorithms to perform multi-level decomposition on the charging load data of electric vehicles, which can more accurately extract the key signal components reflecting the charging load of electric vehicles. This refined signal processing method helps to capture the subtle changes in electric vehicle charging at different usage stages, thereby improving the accuracy of charging load prediction.

[0152] (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 of the prediction model and the overall prediction performance;

[0153] (3) Optimize the selection and application of prediction models. 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 prediction. BRIEF DESCRIPTION OF THE DRAWINGS

[0154] Figure 1 It is the overall system flow chart of the present invention.

[0155] Figure 2 Flowchart of the AO-MTS-Mixers prediction model in the present invention.

[0156] Figure 3 Flowchart of the AO-Crossformer prediction model in the present invention.

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

[0158] The technical solution of the present invention will be further described below with reference to the accompanying drawings.

[0159] See attached Figures 1 to 4 FIG. 1 shows a method for predicting charging load of electric vehicles according to the present invention, which is characterized by comprising the following steps:

[0160] S1: Collect historical charging load data of electric vehicles. The data characteristics of historical charging load data include the starting state of charge, the charging end state of charge, the charging time distribution, and the vehicle charging interval. The ICEEMDAN model is used to decompose the historical charging load data to obtain 10 IMF components.

[0161] The specific steps include:

[0162] S1.1: Set the original charging load time series data signal to X = {x1, x2, ...x I}, where I represents the signal length, and a set of independent Gaussian white noise E1(ω (i) ), the specific formula is as follows:

[0163] X (i) =X+β0E1(ω (i) )

[0164] Among them, ω (i) represents the i-th white noise added, i = 1, 2, 3, ..., I, E(·) represents the first intrinsic mode function IMF1 component of the ICEEMDAN decomposition of a signal; β0 represents the weight coefficient when the noise component is added;

[0165] S1.2: For each signal X after adding noise (i) , introduce the local mean operator M(·) to calculate its corresponding k-th order IMF k Component, the k-th order IMF is obtained by subtracting the local mean from the original signal k Quantity;

[0166] S1.3: Repeat the decomposition steps until the predetermined total decomposition order K is reached, k = 2, 3...K, K represents the number of IMFs, K = 10, and find the kth IMF k Quantity;

[0167]

[0168] 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-1 order weight coefficient;

[0169] S1.4: Repeat the above process for multiple iterations, and the original signal X is decomposed into K IMF components, i.e. 10 IMF i Components, that is, the original signal X is decomposed into the sum of 10 different frequency components.

[0170] S2: Calculate the sample entropy of the 10 components obtained from the first decomposition, and use K-medoids clustering to classify the signal based on the sample entropy results to obtain three signal components: high-frequency component F1, medium-frequency component F2, and low-frequency component F3;

[0171] The specific implementation process of sample entropy calculation and K-medoids clustering includes the following steps:

[0172] S2.1: Assume that each IMF component is represented as IMF n(t) , n∈1,2,...,10, each IMF component IMF n(t) Divide into continuous subsequences of m dimensions:

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

[0174] x i ={x i ,x i+1 ,…,x i+m-1}, i=1~Nm-1

[0175] Where N represents IMF n(t) length;

[0176] S2.2: Define the sequence x i and the sequence x j The maximum distance between corresponding elements is d max (x i ,x j ), calculate the Euclidean distance between each subsequence and other subsequences. The specific formula is as follows:

[0177]

[0178] Among them, d max (x i ,x j ) represents the subsequence x i and x j The maximum distance between corresponding elements, Represents a subsequence x i Middle elements and subsequence x j Middle The absolute value of the difference between the elements;

[0179] S2.3: Calculate the similarity between subsequences statistics The number of vectors is compared with the total number of vectors Nm. The specific formula is as follows:

[0180]

[0181] in express the number of

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

[0183]

[0184] in, represents the similarity of the i-th subsequence; Represents the average similarity of m-dimensional vectors; represents the sample entropy of the m-dimensional i-th subsequence;

[0185] S2.5: Increase the dimension and then add 1 to the dimension m, that is, get the similarity for m+1 vectors And sample entropy

[0186] S2.6: When m reaches the preset value or the sample entropy stops changing, the iteration stops and the sample entropies SE1, SE2, …, SE 10 Composed of a sample entropy vector S i , S i ={SE1,SE2,…,SE 10};

[0187] S2.7: From the sample entropy vector S i Randomly select three objects o1, o2 and o3, representing the initial category centers of high frequency, medium frequency and low frequency respectively and 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. j ,j=1,2,3, the distance calculation formula is as follows:

[0188]

[0189] Among them i represents the i-th sample entropy vector; represents the jth center point, j∈[1,3]; Represents the dimension of the data object;

[0190] S2.8: Introduce formula J to measure the clustering effect:

[0191]

[0192] in, Indicates the oth i and The spatial Euclidean distance of C j represents the jth cluster; K represents the total number of sample entropy;

[0193] S2.9: For each cluster C j , considering all non-central points in the cluster Calculate one of the non-center points When the center point is replaced, the total distance from all points in the cluster to the new center point is recorded as J′;

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

[0195] S2.10: 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, and finally output 3 center points and its corresponding cluster C j ,j=1,2,3, representing the three categories of high-frequency, medium-frequency and low-frequency signals F1, F2 and F3 respectively.

[0196] S3: The high-frequency component F1 is decomposed again using the FMD eigenmode decomposition model to further extract the detailed features of the high-frequency signal and obtain four component signals f1, f2, f3 and f4. The specific steps include:

[0197] S3.1: Load the original high-frequency signal F1 and enter the preset parameters: filter length L, number of frequency band divisions H, and number of hierarchical modes ψ, where ψ = 4. Initialize the FIR filter using Z filters and start iteration, with the initial iteration i = 1.

[0198] S3.2: Divide the frequency band of the original signal into Z segments evenly. The upper and lower cutoff frequency formulas are as follows:

[0199]

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

[0201] S3.3: Pass-through Iterate to obtain the filtered modal component signal, where represents the zth modal component signal of the i-th iteration; represents the response of the zth filter at the i-th iteration, z = 1, 2, ..., Z-1, and "*" represents the convolution operation;

[0202] S3.4: Use the original high-frequency signal F1 to 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 moment of estimation is obtained, one iteration is completed, and i=i+1 is set;

[0203] S3.5: Determine whether the number of iterations reaches the pre-iteration number. If not, continue to update the filter and estimate the cycle; otherwise, proceed to the next step.

[0204] S3.6: 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 coefficient, eliminate the mode with relatively small correlation kurtosis, and at this time there are Z-1 modes remaining;

[0205]

[0206] Among them CC (p×q) represents the correlation coefficient of the two modes; u p (n) and u q (n) represents any two selected modes respectively; and Represent the means of the two modes respectively;

[0207] S3.7: Discard the filter corresponding to the eliminated mode and iteratively update the correlation coefficients of the filters corresponding to other modes. Repeat the above steps until the number of remaining modal components is consistent with the preset ψ, that is, (Z-1) = ψ, and stop the iteration;

[0208] S3.8: If the conditions are not met, return to continue the Hanning window initialization FIR filter bank and use the remaining Z-1 modal components, that is, ψ modal components, as the final component signals of the characteristic mode decomposition, that is, 4 component signals f1, f2, f3 and f4.

[0209] S4: The four component signals f1, f2, f3 and f4 of the high-frequency component F1 are predicted in batches using the MTS-Mixers prediction model, the intermediate-frequency component F2 is predicted using the Crossformer prediction model, and the low-frequency component F3 is predicted using the DeepESN prediction model. In order to improve the prediction performance of the prediction model, the improved artemisinin AO optimization algorithm is used to adjust the time step μ, the number of channels υ and the prediction step ρ in the MTS-Mixers model, and the data partition length L in the Crossformer model. seg , the number of layers l and the total dimension D, the storage layer leakage rate γ and the storage layer size P in the DeepESN model are optimized respectively.

[0210] The improved artemisinin AO optimization algorithm is used to optimize the three hyperparameters of the MTS-Mixers model, namely, the time step μ, the number of channels υ, and the prediction step ρ, and the data partition length L in the Crossformer model. seg , the number of layers l and the total dimension D, as well as the specific steps for optimizing the storage layer leakage rate γ and storage layer size P in the DeepESN model are:

[0211] S4.11: Randomly initialize and generate a group of candidate solutions. 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, and the total dimension is D. The specific formula is:

[0212]

[0213] Where Hb and Lb represent the upper and lower bounds of the solution space respectively; Rb represents a set of random number sequences whose value range is [0,1], and the initial solution is generated using the random number sequence;

[0214] S4.12: Introduce the Tent map to generate chaotic sequences and perform initialization and update optimization on the population so that the initial solutions are distributed as evenly as possible in the solution space. The specific formula is as follows:

[0215]

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

[0217] S4.13: The population is updated to:

[0218] A Nb,D =Lb+Rb'×(Hb-Lb)

[0219] 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 group.

[0220] AO position is updated and enters the comprehensive elimination stage, introducing probability factor Q and concentration factor As a guide to explore the global solution space, when When , the individual position update formula is as follows:

[0221]

[0222] in, represents the position of the i-th individual in the j-th dimension at the t+1-th iteration; Indicates the position of the i-th individual in the j-th dimension in the t-th iteration; t represents the current number of iterations; represents the position of the optimal individual in the jth dimension at the tth iteration; Represents a random number in the range [0,1];

[0223] Concentration factor The specific formula is:

[0224]

[0225] Where t represents the current number of iterations; Maxt is the maximum number of iterations;

[0226] S4.14: Concentration Factor It presents nonlinear attenuation, introduces a piecewise nonlinear convergence factor, slows down the factor attenuation rate based on the sine function, and optimizes the global search capability of the AO algorithm. The specific formula is as follows:

[0227]

[0228] S4.15: AO performs local search and enters the local clearance phase. It introduces the step size coefficient σ and continues to update the individual position to perform local search. It performs precise search in the searched solution space area to find the optimal solution. The specific formula is as follows:

[0229]

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

[0231] in, represents the position of the i-th individual in the t+1th iteration; represents the position of an individual randomly selected in the tth iteration; and Indicates the positions of the other two individuals randomly selected in the tth iteration, and b1, b2 and b3 are random and different indices; Fit norm (i) represents the normalized fitness value of the i-th individual; fit(i) represents the fitness value of the i-th individual; min(fit) and max(fit) represent the minimum and maximum fitness values, respectively; the step coefficient σ is a random value between [0.1, 0.6];

[0232] S4.16: To perform a more detailed search in a local area, an adaptive inertia weight is introduced to optimize the step size coefficient σ of the AO algorithm. The specific formula is as follows:

[0233]

[0234] Where σ(t) represents the step size coefficient in the tth iteration; b represents the parameter that controls the upper and lower limits of the weight; φ represents the parameter that controls the search smoothness;

[0235] S4.17: In the post-consolidation phase, to prevent the AO algorithm from falling into the local optimal solution too early, the algorithm's escape ability is increased by randomly choosing to maintain the current individual position or jump to the global optimal position. The specific formula is:

[0236]

[0237] Among them, best i,j Represents the subvector of the current optimal solution of the j-th dimension;

[0238] S4.18: Determine whether the iterative optimization meets the maximum iteration condition. If it does, the iteration stops. Otherwise, the loop continues until the optimal solution is obtained. These are the three hyperparameters of the MTS-Mixers prediction model: time step μ, number of channels υ, and prediction step ρ.

[0239] S4.19: Repeat S4.11-S4.18 to obtain the data partition length L in the Crossformer model seg , the number of layers l and the total dimension D, as well as the storage layer leakage rate γ and storage layer size P in the DeepESN model.

[0240] The three hyperparameters of the MTS-Mixers prediction model obtained by the improved AO optimization algorithm, namely, time step μ, number of channels υ, and prediction step ρ, are substituted into the MTS-Mixers prediction model to predict the four component signals f1, f2, f3, and f4 of the high-frequency component F1 in batches. Specifically, the following steps are included:

[0241] S4.21: Set the input high-frequency component multivariate time series to χ h, the high-frequency component dataset is divided into training set and test set in a ratio of 7:3, Where μ represents the sequence length, ν represents the number of channels in the time series, and the output sequence is Where ρ represents the prediction step size;

[0242] S4.22: Normalize the input time series. The specific formula is as follows:

[0243]

[0244] in, represents the standardized data; χ h represents the high-frequency component of the input multivariate time series; norm(·) represents data normalization;

[0245] S4.23: Perform encoding on the obtained normalized data to obtain the processed input sequence. The specific formula is as follows:

[0246]

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

[0248] S5.24: 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:

[0249]

[0250] 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;

[0251] S4.25: The extracted and 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 can consider both time series data and time series features. The specific formula is as follows:

[0252]

[0253] 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;

[0254] S4.26: Finally, a linear projection layer performs information fusion and output mapping. The specific formula is as follows:

[0255]

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

[0257] The data in the Crossformer model obtained by the improved AO optimization algorithm is divided into lengths L seg The specific implementation process of substituting the number of layers l and the total dimension D into the Crossformer prediction model to predict the intermediate frequency component F2 is as follows:

[0258] S4.31: Set the intermediate frequency component F2(t) time series data set as the input data matrix The predicted output is Where T ime represents a past time point, τ time represents a future time point, D represents the total dimension, represents a matrix;

[0259] S4.32: Use DSW technology to divide the data in each dimension into L seg , the specific formula is as follows:

[0260]

[0261] where X 1:T Indicates original historical data; X t,d Represents the data at time t in dimension d in the original historical data; Indicates that the length in d dimension is L seg The i-th data segment;

[0262] S4.33: Use linear projection to map the divided raw data into the embedding space, and give each small segment a time position code:

[0263]

[0264] where h i,d represents a time series segment of a data vector; represents the learnable projection matrix; Indicates the position encoding of the i-th data segment in the d-th dimension;

[0265] S4.34: The linear projection of the segmented data and its positional encoding are added to obtain the two-dimensional vector array H. The specific formula is as follows:

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

[0267] S4.35: Two-Stage Attention 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 the MSA operation is performed on different dimensions separately:

[0268]

[0269] in Represents the MSA layer operation result of dimension d; Z d Represents the data of dimension d at all times; LayerNorm(·) represents the normalization layer operation;

[0270] S4.36: Perform nonlinear transformations using a multilayer perceptron. The specific formula is as follows:

[0271]

[0272] where Z time Represents the nonlinear change result of multilayer perceptron; Represents the result of multi-layer perceptron operation;

[0273] S4.37: Through Z time Capture the cross-time dependencies of the same dimension, and then Z time As input to the cross-dimensional stage to obtain the dependencies between different dimensions:

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

[0275] where Z dim represents the output of the multilayer perceptron; Z, Z represents the input matrix of the TSA layer, and Y represents the output matrix of the TSA layer;

[0276] S4.38: In each layer of the encoder except the first layer, two adjacent vectors in the time domain are merged to form a new layer, and the dependencies of this layer are captured using TSA. After obtaining the N+1 output results of the encoder, the decoder uses the N+1 layers (0, 1, ..., N) for prediction. The specific formula is as follows:

[0277]

[0278] Where l represents the number of layers, Represents the i-th data segment in the prediction result of the l-th layer d dimension; represents a learnable matrix; represents the output of the decoder; Represents the prediction result of the lth layer; Indicates the final prediction result.

[0279] Substituting the storage layer leakage rate γ and storage layer size P in the DeepESN model obtained by the improved AO optimization algorithm into the DeepESN model prediction model to predict the low-frequency component F3 is as follows:

[0280] S4.41: Assume that DeepPESN consists of P storage layers and P-1 encoder layers to form a P-layer network. Assume that there are N(i) neurons in the i-th storage layer and M(j) neurons in the j-th encoder layer, where 1≤i≤P, 1≤j≤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:

[0281]

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

[0283] S4.42: The input of the i-th storage layer is the output of the i-1-th encoder layer, introducing the function To express the conforming function, including high-dimensional projection and updating the state of the storage layer, the formula is updated as follows:

[0284]

[0285] in represents a composite function;

[0286] S4.43: Given the previous state of the reservoir, use unsupervised dimensionality reduction techniques Encode and generate encoding features. The encoding process of the j-th encoder is expressed as follows:

[0287]

[0288] in, represents the activation value of the j-th encoder layer at time t, 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(i) represents the input weight of the j-th encoder layer;

[0289] S4.44: Based on the states of the storage layers from layer 1 to layer P-1, the state of the last storage layer at layer P is obtained. The specific formula is as follows:

[0290]

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

[0292] S4.45: Output y(t+1) at time t+1, specifically:

[0293]

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

[0295] S4.46: Through the matrix W out With the input feature matrix The specific formula is as follows:

[0296]

[0297] in Represents the final output matrix.

[0298] S5: Reconstruct the prediction results of the three models to obtain the prediction result of the charging load, which specifically includes the following steps:

[0299] The prediction results of the high-frequency component F1(t) obtained by using the AO-MTS-Mixers prediction model The prediction results of the intermediate frequency component F2(t) obtained using the AO-Crossformer prediction model And the prediction results of AO-DeepESN low-frequency component F3(t) The prediction result of charging load is obtained by using a weighted sum composite prediction model to combine the results. The specific formula is:

[0300]

[0301] Wherein, 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]. In this embodiment, α1∈[0.9, 1], α2∈[0.6, 0.9], and α3∈[0, 0.6].

[0302] S6: To evaluate the training effect of the prediction model, error analysis is required. The following error parameter evaluation indicators are used to verify the performance of the model: relative error RE, mean square relative error MRE, absolute error AE, and mean absolute error MAE. The calculation formulas for the corresponding indicators are as follows:

[0303] Relative error RE:

[0304]

[0305] Mean square relative error MRE:

[0306]

[0307] Absolute error AE:

[0308]

[0309] Mean Absolute Error (MAE):

[0310]

[0311] 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.

Claims

1. A method for predicting charging load of electric vehicles, characterized by: The specific steps include: S1: Collect historical charging load data of electric vehicles. The data characteristics of historical charging load data include the starting state of charge, the charging end state of charge, the charging time distribution, and the vehicle charging interval. Use the ICEEMDAN model to decompose the historical charging load data to obtain several IMF components; S2: Calculate the sample entropy of several IMF components obtained by the primary decomposition, and use K-medoids clustering to classify the signal based on the sample entropy results to obtain three signal components: high-frequency component F1, medium-frequency component F2, and low-frequency component F3; S3: The high-frequency component F1 is decomposed twice using the FMD eigenmode decomposition model to obtain four component signals f1, f2, f3 and f4; S4: The MTS-Mixers prediction model is used to predict the four component signals f1, f2, f3 and f4 of the high-frequency component F1 in batches, the Crossformer prediction model is used to predict the mid-frequency component F2, and the DeepESN prediction model is used to predict the low-frequency component F3; S5: Reconstruct the prediction results of the three models to obtain the prediction result of the charging load.

2. The method for predicting charging load of electric vehicles according to claim 1, characterized in that: The specific process of the ICEEMDAN model in S1 decomposing the historical charging load data to obtain several IMF components includes the following steps: S1.1: Set the original charging load time series data signal to X = {x1, x2, ...x I }, where I represents the signal length, and a set of independent Gaussian white noise E1(ω (i) ), the specific formula is as follows: X (i) =X+β0E1(ω (i) ) Among them, ω (i) represents the i-th white noise added, i = 1, 2, 3, ..., I, E(·) represents the first intrinsic mode function IMF1 component of the ICEEMDAN decomposition of a signal; β0 represents the weight coefficient when the noise component is added; S1.2: For each signal X after adding noise (i) , introduce the local mean operator M(·) to calculate its corresponding k-th order IMF k Component, the k-th order IMF is obtained by subtracting the local mean from the original signal k Quantity; S1.3: Repeat the decomposition steps until the predetermined total decomposition order K is reached, k = 2, 3…K, where K represents the number of IMFs. Find the kth IMF k Quantity; 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-1 order weight coefficient; S1.4: Repeat the above process for multiple iterations, and the original signal X is decomposed into K IMF components, i.e. 10 IMF i Components, that is, the original signal X is decomposed into the sum of several different frequency components.

3. The method for predicting charging load of electric vehicles according to claim 1, characterized in that: The specific implementation process of sample entropy calculation and K-medoids clustering in S2 includes the following steps: S2.1: Assume that each IMF component is represented as IMF n(t) , n∈1,2,...,N, each IMF component IMF n(t) Divide into continuous subsequences of m dimensions: IMF n(t) ={x1,x2,…,x N } x i ={x i ,x i+1 ,…,x i+m-1 },i=1~N-m-1 Where N represents IMF n(t) length; S2.2: Define the sequence x i and the sequence x j The maximum distance between corresponding elements is d max (x i ,x j ), calculate the Euclidean distance between each subsequence and other subsequences. The specific formula is as follows: Among them, d max (x i ,x j ) represents the subsequence x i and x j The maximum distance between corresponding elements, Represents a subsequence x i Middle elements and subsequence x j Middle The absolute value of the difference between the elements; S2.3: Calculate the similarity between subsequences statistics The number of vectors is compared with the total number of vectors Nm. The specific formula is as follows: in express the number of S2.4: Calculate the sample entropy of each subsequence, The results are averaged and logarithmically transformed to obtain the sample entropy The specific formula is as follows: in, represents the similarity of the i-th subsequence; Represents the average similarity of m-dimensional vectors; represents the sample entropy of the m-dimensional i-th subsequence; S2.5: Increase the dimension and then add 1 to the dimension m, that is, get the similarity for m+1 vectors And sample entropy S2.6: When m reaches the preset value or the sample entropy stops changing, the iteration stops and the sample entropies SE1, SE2, ..., SE of the calculated IMF components are converted to n Composed of a sample entropy vector S i , S i ={SE1,SE2,…,SE n }; S2.7: From the sample entropy vector S i Randomly select three objects o1, o2 and o3, representing the initial category centers of high frequency, medium frequency and low frequency respectively and Calculate the distance between the remaining objects and the center point of each initial category, and assign each object to the cluster C represented by the nearest initial category center. j ,j=1,2,3, the distance calculation formula is as follows: Among them i represents the i-th sample entropy vector; represents the jth center point, j∈[1,3]; Represents the dimension of the data object; S2.8: Introduce formula J to measure the clustering effect: in, Indicates the oth i and The spatial Euclidean distance of C j represents the jth cluster; K represents the total number of sample entropy; S2.9: For each cluster C j , considering all non-central points in the cluster Calculate one of the non-center points When the center point is replaced, the total distance from all points in the cluster to the new center point is recorded as J′; If the total distance decreases after replacement, that is, J'<J, then replace Otherwise, keep the original center point unchanged; S2.10: 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, and finally output 3 center points and its corresponding cluster C j ,j=1,2,3, representing the three categories of high-frequency, medium-frequency and low-frequency signals F1, F2 and F3 respectively.

4. The method for predicting charging load of electric vehicles according to claim 1, wherein: The S3 specifically includes the following steps: S3.1: Load the original high-frequency signal F1 and enter the preset parameters: filter length L, number of frequency band divisions H, and number of hierarchical modes ψ, where ψ = 4. Initialize the FIR filter using Z filters and start iteration, with the initial iteration i = 1. S3.2: Divide the frequency band of the original signal into Z segments evenly. The upper and lower cutoff frequency formulas are as follows: Among them, f l Indicates the lowest cutoff frequency; f u Indicates the highest cutoff frequency; f s Indicates the sampling frequency of the original high-frequency signal; S3.3: Pass-through Iterate to obtain the filtered modal component signal, where represents the zth modal component signal of the i-th iteration; represents the response of the zth filter at the i-th iteration, z = 1, 2, ..., Z-1, and "*" represents the convolution operation; S3.4: Use the original high-frequency signal F1 to 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 moment of estimation is obtained, one iteration is completed, and i=i+1 is set; S3.5: Determine whether the number of iterations reaches the pre-iteration number. If not, continue to update the filter and estimate the cycle; otherwise, proceed to the next step. S3.6: 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 coefficient, eliminate the mode with relatively small correlation kurtosis, and at this time there are Z-1 modes remaining; S3.7: Discard the filter corresponding to the eliminated mode and iteratively update the correlation coefficients of the filters corresponding to other modes. Repeat the above steps until the number of remaining modal components is consistent with the preset ψ, that is, (Z-1) = ψ, and stop the iteration; S3.8: If the conditions are not met, return to continue the Hanning window initialization FIR filter bank and use the remaining Z-1 modal components, that is, ψ modal components, as the final component signals of the characteristic mode decomposition, that is, 4 component signals f1, f2, f3 and f4.

5. The method for predicting charging load of electric vehicles according to claim 1, characterized in that: In order to improve the prediction performance of the prediction model in S4, the improved artemisinin AO optimization algorithm was used to adjust the time step μ, channel number υ and prediction step ρ in the MTS-Mixers model, and the data partition length L in the Crossformer model. seg , the number of layers l and the total dimension D, the storage layer leakage rate γ and the storage layer size P in the DeepESN model are optimized respectively.

6. The method for predicting charging load of electric vehicles according to claim 5, characterized in that: The improved artemisinin AO optimization algorithm is used to optimize the three hyperparameters of the MTS-Mixers model: time step μ, number of channels υ and prediction step ρ, and the data partition length L in the Crossformer model. seg , the number of layers l and the total dimension D, as well as the specific steps for optimizing the storage layer leakage rate γ and storage layer size P in the DeepESN model are: S4.11: Randomly initialize and generate a group of candidate solutions. 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, and the total dimension is D. The specific formula is: Where Hb and Lb represent the upper and lower bounds of the solution space range respectively; Rb represents a sequence of random numbers, whose value range is [0, 1], and the initial solution is generated using the random number sequence. S4.12: Introduce the Tent mapping to generate a chaotic sequence, initialize and update the population for optimization, so that the initial solution is distributed as evenly as possible in the solution space. The specific formula is as follows: Where Rb' represents the chaotic sequence introduced by the Tent mapping; τ represents the control parameter, τ ∈ rand(0, 1); S4.13: The population is updated as follows: A Nb,D =Lb+Rb'×(Hb-Lb) 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 group. The position of AO is updated, entering the comprehensive elimination stage. Introduce the probability factor Q and the concentration factor θ as guides to explore the global solution space. When θ1 < Q, the individual position update formula is as follows: in, represents the position of the i-th individual in the j-th dimension at the t+1-th iteration; Indicates the position of the i-th individual in the j-th dimension in the t-th iteration; t represents the current number of iterations; represents the position of the optimal individual in the jth dimension at the tth iteration; θ1 represents a random number in the range [0,1]; The specific formula for the concentration factor θ is: Where t represents the current iteration number; Maxt is the maximum iteration number; S4.14: Introduce a piecewise non-linear convergence factor, based on the sine function to slow down the decay rate of the factor, and optimize the global search ability of the AO algorithm. The specific formula is as follows: S4.15: AO performs local search, entering the local elimination stage. Introduce the step size coefficient σ, and continue to update the individual position for local search, and perform precise search within the searched solution space area to find the optimal solution. The specific formula is as follows: in, represents the position of the i-th individual in the t+1th iteration; represents the position of an individual randomly selected in the tth iteration; and Indicates the positions of the other two individuals randomly selected in the tth iteration, and b1, b2 and b3 are random and different indices; Fit norm (i) represents the normalized fitness value of the i-th individual; fit(i) represents the fitness value of the i-th individual; min(fit) and max(fit) represent the minimum and maximum fitness values, respectively; the step coefficient σ is a random value between [0.1, 0.6]; b1, b2, b3 ∼ U(1, N), b1 ≠ b2 ≠ b3 S4.16: Introduce an adaptive inertia weight to optimize the step size coefficient σ of the AO algorithm. The specific formula is as follows: Where, σ(t) represents the step size coefficient in the t-th iteration; b represents the parameter that controls the upper and lower limits of the weight; φ represents the parameter that controls the search smoothness; Among them, best i,j Represents the subvector of the current optimal solution of the j-th dimension; S4.17: To prevent the AO algorithm from falling into the local optimal solution prematurely, randomly select to keep the current individual position or jump to the global optimal position to increase the escape ability of the algorithm. The specific formula is: S4.19: Repeat S4.11-S4.18 to obtain the data partition length L in the Crossformer model seg , the number of layers l and the total dimension D, as well as the storage layer leakage rate γ and storage layer size P in the DeepESN model.

7. The method for predicting charging load of electric vehicles according to claim 6, characterized in that: S4.18: Judge whether the iterative optimization meets the maximum iteration condition. If it is met, the iteration stops; otherwise, continue to loop until the optimal solution is obtained, that is, the three hyperparameters of the MTS-Mixers prediction model: the time step μ, the number of channels υ, and the prediction step ρ; S4.21: Set the input high-frequency component multivariate time series to χ h , Where μ represents the sequence length, v represents the number of channels of the time series, and the output sequence is Where ρ represents the prediction step size; Substitute the three hyperparameters of the MTS-Mixers prediction model, the time step μ, the number of channels υ, and the prediction step ρ obtained by the improved AO optimization algorithm into the MTS-Mixers prediction model to predict the four component signals f1, f2, f3, and f4 of the high-frequency component F1 in batches. The specific steps are as follows: in, represents the standardized data; χ h represents the high-frequency component of the input multivariate time series; norm(·) represents data normalization; S4.22: Standardize the input time series. The specific formula is as follows: in, represents the input sequence after encoding operation; PE(·) represents sinusoidal position encoding; S4.23: Perform an encoding operation on the obtained standardized data to obtain the processed input sequence. The specific formula is as follows: S5.24: After obtaining the standardized and encoded data, use the attention mechanism to perform information feature learning on the input data sequence in the time dimension to capture the features of the time series. The specific 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; S4.25: The extracted and 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 can consider both time series data and time series features. The specific 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; S4.26: Finally, a linear projection layer performs information fusion and output mapping. The specific formula is as follows: in, represents the predicted output sequence of the i-th high-frequency component, i = 1, 2, 3, 4, and Linear(·) represents the output operation after information fusion and linear mapping of the two sequences.

8. The method for predicting charging load of electric vehicles according to claim 6, characterized in that: The data in the Crossformer model obtained by the improved AO optimization algorithm is divided into lengths L seg The specific implementation process of substituting the number of layers l and the total dimension D into the Crossformer prediction model to predict the intermediate frequency component F2 is as follows: S4.31: Set the intermediate frequency component F2(t) time series data set as the input data matrix The predicted output is Where T ime represents a past time point, τ time represents a future time point, D represents the total dimension, represents a matrix; S4.32: Use DSW technology to divide the data in each dimension into L seg , the specific formula is as follows: where X 1:T Indicates original historical data; X t,d Represents the data at time t in dimension d in the original historical data; Indicates that the length in d dimension is L seg The i-th data segment; S4.33: Use linear projection to map the divided raw data into the embedding space, and give each small segment a time position code: where h i,d represents a time series segment of a data vector; represents the learnable projection matrix; Indicates the position encoding of the i-th data segment in the d-th dimension; S4.34: The linear projection of the segmented data and its positional encoding are added to obtain the two-dimensional vector array H. The specific formula is as follows: H={h i,d |1≤i≤T ime / L seg ,1≤d≤D} S4.35: Two-Stage Attention 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 the MSA operation is performed on different dimensions separately: in Represents the MSA layer operation result of dimension d; Z d Represents the data of dimension d at all times; LayerNorm(·) represents the normalization layer operation; S4.36: Perform nonlinear transformations using a multilayer perceptron. The specific formula is as follows: where Z time Represents the nonlinear change result of multilayer perceptron; Represents the result of multi-layer perceptron operation; S4.37: Through Z time Capture the cross-time dependencies of the same dimension, and then Z time As input to the cross-dimensional stage to obtain the dependencies between different dimensions: Y=Z dim =TSA(Z) where 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; S4.38: In each layer of the encoder except the first layer, two adjacent vectors in the time domain are merged to form a new layer, and the dependencies of this layer are captured using TSA. After obtaining the N+1 output results of the encoder, the decoder uses the N+1 layers (0, 1, ..., N) for prediction. The specific formula is as follows: Where l represents the number of layers, l = 0,...,N; Represents the i-th data segment in the prediction result of the l-th layer d dimension; represents a learnable matrix; represents the output of the decoder; Represents the prediction result of the lth layer; Indicates the final prediction result.

9. The method for predicting charging load of electric vehicles according to claim 6, characterized in that: Substituting the storage layer leakage rate γ and storage layer size P in the DeepESN model obtained by the improved AO optimization algorithm into the DeepESN model prediction model to predict the low-frequency component F3 is as follows: S4.41: Assume that DeepPESN consists of P storage layers and P-1 encoder layers to form a P-layer network. Assume that there are N(i) neurons in the i-th storage layer and M(j) neurons in the j-th encoder layer, where 1≤i≤P, 1≤j≤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 (i) (t) represents the internal state of the i-th storage layer at time t; f(·) represents the nonlinear activation function in the storage layer; represents the input of the i-th storage layer at time t+1; represents the activation value of the i-th storage layer at time t; W represents the activation value of the i-th storage layer at time t+1, 1≤t≤T; in(i) represents the input weight of the i-th storage layer; W res(i) represents the cyclic weight of the i-th storage layer; γ represents the leakage rate of the storage layer; S4.42: The input of the i-th storage layer is the output of the i-1-th encoder layer, introducing the function To express the conforming function, including high-dimensional projection and updating the state of the storage layer, the formula is updated as follows: in represents a composite function; S4.43: Given the previous state of the reservoir, use unsupervised dimensionality reduction techniques Encode and generate encoding features. The encoding process of the j-th encoder is expressed as follows: in, represents the activation value of the j-th encoder layer at time t, 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(i) represents the input weight of the j-th encoder layer; S4.44: Based on the states of the storage layers from layer 1 to layer P-1, the state of the last storage layer at layer P is obtained. The specific formula is as follows: in represents the activation value of the Pth storage layer at time t+1; Chinese symbol Represents a combination operator; Represents the high-dimensional projection of the j-th storage layer and the state composite function of the updated storage layer; Represents the encoding process of the jth encoding layer; represents the coding function used by the jth coding layer; u(t+1) represents the input signal at time t+1; S4.45: Output y(t+1) at time t+1, the specific formula is: Among them, W out represents the output weight; Represents the input feature matrix at time t+1; S4.46: Through the matrix W out With the input feature matrix The specific formula is as follows: in Represents the final output matrix.

10. The method for predicting charging load of electric vehicles according to claim 1, characterized in that: The specific implementation process of S5 is as follows: The prediction results of the high-frequency component F1(t) obtained by using the AO-MTS-Mixers prediction model The prediction results of the intermediate frequency component F2(t) obtained using the AO-Crossformer prediction model And the prediction results of AO-DeepESN low-frequency component F3(t) The prediction result of charging load is obtained by using a weighted sum composite prediction model to combine the results. The specific formula is: Where Y represents the predicted value of the charging load by the composite prediction model, and α1, α2, and α3 represent random fusion coefficients in the range [0, 1].

Images