Method for predicting residual life of lithium battery

Abstract

The invention discloses a method for predicting the residual life of a lithium battery. The method belongs to the field of lithium battery health management, and comprises the following operation steps: collecting operation data of a lithium battery, and performing primary decomposition on a lithium battery capacity sequence by using a TVFEMD model; performing secondary decomposition on high-frequency and residual components in the IMF signal by using an EWT model; performing sample entropy calculation on the medium-frequency and low-frequency components after primary decomposition and the high-frequency, medium-frequency and low-frequency components after secondary decomposition, and performing clustering analysis on a classified component set by adopting a K-means clustering method; predicting a high-frequency component after K-means clustering by using a BKA-GLIP model, and predicting an intermediate-frequency component and a low-frequency component by using a BKA-ESN model; and reconstructing prediction results of the two models to obtain a prediction result of the remaining service life of the lithium battery. According to the method, the accuracy of predicting the remaining service life of the lithium battery is improved, and powerful support is provided for use and management of the lithium battery.

Description

Technical Field

[0001] The invention belongs to the technical field of lithium battery remaining life prediction and relates to a method for lithium battery remaining life prediction. Background Art

[0002] As the development and utilization of new energy sources have attracted more and more attention; in the field of new energy, lithium-ion batteries are widely used in various electronic devices and electric vehicles due to their high energy density, long life and environmental protection; however, during the use of lithium batteries, their performance will gradually decline until they can no longer meet the use needs. Therefore, predicting the remaining life of lithium batteries is of great significance to ensure the safe operation of equipment and improve the utilization rate of batteries.

[0003] At present, the commonly used methods for predicting the remaining life of lithium batteries mainly include model-based prediction methods and data-based prediction methods; the model-based prediction method is to predict the remaining life of the battery based on the battery operating data and model parameters by establishing a mathematical model of the battery; however, due to the complexity and uncertainty of the internal process of the battery, it is difficult to establish an accurate battery model; the data-based prediction method is to collect the operating data of the battery, such as voltage, current, charge and discharge time, temperature and capacity, and then use a data-driven method to predict the remaining life of the battery; this method does not require the establishment of a complex battery model and therefore has better applicability; but since the operating data of the battery usually has the characteristics of nonlinearity, non-stationarity and uncertainty, how to effectively extract and process these data to improve the accuracy of the prediction is still a challenging problem; at present, traditional data processing methods are difficult to effectively extract the feature information of the data; and most of the commonly used prediction methods are based on a single prediction model and cannot process signals of different frequencies at the same time, resulting in the accuracy of the prediction results being affected. Summary of the invention

[0004] In view of the above problems, the present invention aims to provide a method for predicting the remaining service life of lithium batteries, which can improve the accuracy of predicting the remaining service life of lithium batteries and provide strong support for the use and management of lithium batteries.

[0005] The technical solution of the present invention is: a method for predicting the remaining life of a lithium battery according to the present invention, the operating steps are as follows:

[0006] Step (1): Collecting the operating data of the lithium battery; using the TVFEMD model to decompose the lithium battery capacity sequence, which is divided into high frequency, medium frequency, low frequency and residual components from high to low frequency;

[0007] The operation data includes information on discharge voltage, discharge current, discharge time and capacity sequence;

[0008] Step (2): Use the EWT model to perform secondary decomposition on the high-frequency and residual components in the IMF signal. The IMF components after secondary decomposition are high-frequency, medium-frequency, and low-frequency components according to the frequency from high to low.

[0009] Step (3): Calculate sample entropy for the intermediate frequency and low frequency components after the first decomposition and the high frequency, intermediate frequency and low frequency components after the second decomposition. Classify the components after the second decomposition into high frequency, intermediate frequency and low frequency component sets according to the sample entropy evaluation results. Perform cluster analysis on the classified component sets using the K-means clustering method to obtain three signal components of high frequency, intermediate frequency and low frequency.

[0010] Step (4): The high-frequency components after K-means clustering are predicted using the GLIP prediction model, and the medium-frequency and low-frequency components are predicted using the ESN network;

[0011] Step (5): In order to improve the prediction performance of the two prediction models, the BKA optimization algorithm is used to optimize the hyperparameters α, β, and γ in the GLIP prediction model and the hyperparameters of the reserve pool spectrum radius SR, the reserve pool sparsity SD, and the input unit scale IS in the ESN network;

[0012] Step (6): Finally, the predictions of the two models are added and reconstructed to obtain the remaining service life prediction result of the lithium battery.

[0013] Furthermore, in step (1), in the step of using the TVFEMD model to decompose the lithium battery capacity sequence, the implementation process of the TVFEMD model to decompose the lithium battery capacity sequence is as follows:

[0014] Step (1.1): Set the original battery capacity signal to x(t), and use Hilbert transform to obtain the signal x'(t) after the original battery capacity signal x(t) is transformed. The formula is as follows:

[0015]

[0016] Wherein, X(t) represents the instantaneous amplitude of the battery capacity signal x(t); φ(t) represents the instantaneous phase of the battery capacity signal x(t); x'(t) represents the Hilbert transform signal of the battery capacity signal x(t);

[0017] The formula of the analytical signal h(t) corresponding to the battery capacity signal x(t) is expressed as:

[0018] h(t)=X(t)e jφ(t) =α1e jφ(t) +α2e jφ(t)

[0019] Where α1 and α2 represent the instantaneous amplitude of the component signal. The phase difference represents the instantaneous amplitude;

[0020] The formula for instantaneous amplitude X(t) is:

[0021] X 2 (t) = α1 2 (t)+α2 2 (t)+2α1(t)α2(t)cos[φ1(t)-φ2(t)]

[0022] Where φ1(t) and φ2(t) represent the instantaneous phase of the component signal, and α1(t) and α2(t) represent the instantaneous amplitude of the signal component at the tth moment;

[0023] Step (1.2): Determine the local maximum value X(t) of the instantaneous amplitude X(t) max ) sequence and the local minimum sequence X(t min ), and interpolate the local maximum curve β1(t) and the local minimum curve β2(t) respectively, and then calculate the amplitudes α1(t) and α2(t), the instantaneous frequency components φ1'(t) and φ2'(t), the formula is as follows:

[0024]

[0025] Where β1(t) and β2(t) represent the local maximum and local minimum; α1(t) and α2(t) represent the new amplitude;

[0026] The formula for the instantaneous frequency component φ'(t) of the battery capacity signal x(t) is:

[0027]

[0028] X 2 (t max )φ'(t max ) and X 2 (t min )φ'(t min ) is interpolated to obtain η1(t) and η2(t), then the formulas for the instantaneous frequency components φ1'(t) and φ2'(t) are:

[0029]

[0030] In the formula, η1(t) and η2(t) represent the 2 (t max )φ'(t max ) and X 2 (t min )φ'(t min) is obtained by interpolation;

[0031] Calculate the local cutoff frequency φ of the signal b '(t), the formula is as follows:

[0032]

[0033] Step (1.3): Filter the uncalculated signal using a time-varying filter to obtain the local mean, and calculate the signal z(t); as shown in the following formula:

[0034] z(t)=cos[∫φ' b (t)dt]

[0035] In the formula, φ b '(t) is the local cutoff frequency of the signal; when z(t) is at an extreme value, the time point t is extracted max and t min As node m, a B-spline approximation time-varying filter is constructed, and the signal x(t) is subjected to B-spline approximation time-varying filtering. The approximation result is recorded as m(t), which represents the local mean;

[0036] Continuously update the signal x(t) and define the stopping standard θ(t). When the signal satisfies θ(t)≤ε, ε is the given bandwidth threshold, the signal component at this time is recorded as IMF i component signal; otherwise, let x(t) = x(t) - m(t), and repeat the above steps (1.1)-(1.3) to continue decomposing the signal until the stopping criterion is met and the decomposition is completed;

[0037] Wherein, the formula of the stopping criterion θ(t) is as follows:

[0038]

[0039] Where B(t) represents the instantaneous bandwidth of the signal; φ avg (t) represents the weighted mean instantaneous frequency;

[0040] According to steps (1.1) to (1.3), the TVFEMD model decomposes the lithium battery capacity sequence and extracts a series of IMF signals, which are high frequency, medium frequency, and low frequency three signal components IMF1, IMF2, and IMF3 respectively according to the frequency from high to low;

[0041] Step (1.4): When all IMF signals are extracted, the remaining signal is the residual component R(t); as shown in the following formula:

[0042]

[0043] Furthermore, in step (2), the implementation process of using the EWT model to perform secondary decomposition on the high frequency and residual components in the IMF signal obtained after the primary decomposition is:

[0044] Step (2.1): Perform Fourier transform on the signal IMF1 retained in step (1) and the residual signal R(t);

[0045] Step (2.2): Through the empirical scaling function With the empirical wavelet function ψ n (ω) realize the adaptive properties of H-1 bandpass filters;

[0046] Among them, the empirical scaling function With the empirical wavelet function ψ n The calculation formula of (ω) is as follows:

[0047]

[0048]

[0049] In the formula, represents the empirical scaling function; ψ n (ω) represents the empirical wavelet function; ω represents the frequency, ω n represents the frequency of the nth boundary; β(·) represents the transformation function; τ n represents the coefficient function, τ n =γω n , (0<γ<1);

[0050] Step (2.3): Based on the empirical scaling function The inner product of the empirical wavelet function and the original signal is used to obtain the scaling coefficient, and the inner product of the empirical wavelet function and the original signal is used to obtain the empirical wavelet coefficient. The formula is as follows:

[0051]

[0052] In the formula, t represents the time point; and represents the empirical scaling function and the empirical wavelet function; and represents the Fourier transform of the empirical scaling function and the empirical wavelet function; F -1 [·] denotes inverse Fourier transform;

[0053] Step (2.4): The signal IMF1 and the residual signal R(t) are decomposed and reconstructed by the EWT algorithm as follows:

[0054]

[0055] In the formula, k represents the number of decomposed components, and its value is k=3.

[0056] Furthermore, in step (3), the implementation process of cluster analysis using the K-means clustering method is as follows:

[0057] Step (3.1): Let each IMF component be IMF n(t) , n∈1,2,…8, each IMF component IMF n(t) Divide into m-dimensional continuous subsequences as follows:

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

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

[0060] Where N represents IMF i(t) Length;

[0061] For each subsequence x i , calculate other subsequences x j The Euclidean distance d(x i ,x j ); define the sequence x i With the sequence x j The maximum distance between corresponding elements is d max (x i ,x j ); its formula is as follows:

[0062]

[0063] Step (3.2): Given a threshold r, count d[x i ,x j ]<r and is ratioed to the total number of vectors Nm, recorded as The formula is as follows:

[0064]

[0065] In the formula, num{d[x i ,x j ]<r} means d max (x i ,x j )<the number of r;

[0066] right The average of the results is calculated using the following formula:

[0067]

[0068] SE m (i,r)=-ln[(B m (r))]

[0069] In the formula, represents the similarity of the i-th subsequence; B m (r) represents the average similarity of m-dimensional vectors; SE m (i, r) represents the sample entropy of the m-dimensional i-th subsequence, and then the dimension m is added by 1, that is, the similarity B is obtained for m+1 vectors m+1 (r) and sample entropy SE m+1 (i,r);

[0070] Step (3.3): The sample entropies of the eight IMF components SE1, SE2, …, SE8 are combined into a sample entropy vector S i , S i ={SE1, SE2, …, SE8}; select 3 initial center cluster numbers C1, C2, C3, representing high-frequency, medium-frequency and low-frequency signals respectively; for each sample entropy S i , calculate its relationship with each center C j The distance is usually calculated using the Euclidean distance, which is as follows:

[0071]

[0072] In the formula, S i represents the i-th sample entropy vector; C j represents the number of the jth cluster, j∈[1,3];

[0073] For each sample entropy value S i , find the center C that is closest to it j , and assign it to this center, the formula is as follows:

[0074] C j* =argmin j d(S i ,C j )

[0075] In the formula, C j* Indicates distance S i nearest center;

[0076] Step (3.4): Recalculate the center of each cluster and take the average value of the entropy of all samples assigned to the cluster as the new cluster center. The formula is as follows:

[0077]

[0078] In the formula, |C j | indicates assignment to cluster C j The number of samples;

[0079] Repeat iterative steps (3.3) to (3.4) until the new cluster center is basically unchanged. Finally, each IMF component will be assigned to the corresponding cluster, and the clustering results are output: high-frequency component F1(t), medium-frequency component F2(t) and low-frequency component F3(t).

[0080] Furthermore, in step (4), the implementation process of using the GLIP prediction model to predict the high frequency component is:

[0081] Step (4.1): Initialization of the BKA optimization algorithm: Randomly initialize to generate a group of candidate solutions, each individual solution represents one of the parameters of the GLIP model, including hyperparameters α, β, and γ;

[0082] Among them, α and β represent the global and local weight coefficients respectively; γ represents the acceleration of weight change; the population size of the BKA optimization algorithm is set to N, the maximum number of iterations is T, and the variable dimension is j; at the same time, the f(·) function is introduced to represent the fitness value of the individual. During the initialization process, BKA selects the individual with the best fitness value as the optimal solution X in the initial population. L , that is, setting the optimal solution as the hyperparameter of the GLIP prediction model;

[0083] Step (4.2): Update the BKA position and iteratively optimize the hyperparameters of the GLIP prediction model; in each iteration, update the position of the candidate solution according to the fitness function to find the hyperparameter settings of the GLIP prediction model;

[0084] The formula for updating the BKA position is as follows:

[0085]

[0086] In the formula, and Respectively represent the position solutions of the ith individual in the jth dimension at the tth iteration and the t+1th iteration; r is a random number between [0,1]; p is a constant; t is the number of iterations completed so far; T is the total number of iterations; n is the adhesion coefficient;

[0087] The formula for BKA to continue position updating is as follows:

[0088]

[0089] m=2×sin(r+π / 2)

[0090] In the formula, represents the optimal solution of the j-th dimension individual in the t-th iteration; F i represents the current position of any individual in the jth dimension at the tth iteration; F ri represents the fitness value of any individual at the random position of the jth dimension at the tth iteration; C(0,1) represents the Cauchy mutation;

[0091] Step (4.3): According to the obtained best individual fitness value, BKA selects the best individual and saves it, that is, updates the hyperparameters of the GLIP prediction model;

[0092]

[0093] In the formula, X best represents the individual optimal solution; F best Represents the optimal fitness value of an individual; represents the optimal solution of the j-th dimension individual in the t-th iteration; and They represent the position solutions of the ith individual in the jth dimension at the tth iteration and the t+1th iteration respectively;

[0094] Step (4.4): Determine whether the iterative optimization meets the maximum iteration condition. If it does, the iteration stops. Otherwise, the cycle continues until the optimal solution is obtained, i.e., the three hyperparameters α, β, and γ of the GLIP prediction model.

[0095] Step (4.5): Substitute the three hyperparameters α, β, and γ of the GLIP prediction model obtained by the BKA optimization algorithm into the GLIP prediction model to predict the high-frequency component F1(t); input the high-frequency component F1(t), set its data sequence to X, the input length to L, the output length to O, and the hyperparameters α, β, and γ obtained after the above optimization, and split the high-frequency component data set into the training set X in a ratio of 7:1:2. train , validation set X val and the test set X test ;

[0096] Step (4.6): GLIP global prediction: for the entire training set X train Perform a discrete Fourier transform to obtain the amplitude A and frequency f; then, select the high amplitude frequency f from the amplitude * , then its corresponding period T * =1 / f * , let T * =[T1, T2, ...]; then, using the period T * To construct the global basis function Θ g , the formula is as follows:

[0097] Θ g=[sin(C1t),cos(C1t),1]

[0098] Where, C1 = 2π / T * , t represents the length of the training set, T * Indicates high frequency period;

[0099] The sparse identification process is obtained according to the global basis function, and the global parameter Ξ is solved through sparse identification. g , to achieve global prediction, the optimization formula of sparse recognition is as follows:

[0100]

[0101] Where X represents the training set data sequence; ||X-Θ g Ξ g || 2 Represents the global basis function for the training set X train The fitting error of g || indicates the global parameter Ξ g l1 norm regularization;

[0102] By using the coordinate descent method similar to LASSO to solve the above optimization problem, we get Θ g ,make To predict the validation set X val and the test set X test The data in the prediction result is recorded as X gp ;

[0103] Using the global prediction result X gp , construct a new basis function The formula is as follows:

[0104]

[0105] In the formula, X g =[X train ,X gp ], indicating the univariate prediction result X gp With training set X train splicing;

[0106] The new basis function Alternative Θ g , re-solve the optimization problem of sparse identification and obtain the global prediction result considering the network coupling relationship between variables

[0107] Step (4.7): Verify the recognition and determine whether to perform local rolling prediction: The GLIP prediction model is used for verification, to evaluate the effectiveness of global recognition and prediction, to observe whether the use of local features will perform better than global prediction, and to prepare for subsequent local rolling prediction;

[0108] Assume that the global prediction value for the validation set data is Calculate the forecast error for each variable i and compare it with the local rolling forecast The error is compared; the implementation formula is as follows:

[0109]

[0110] Where N b Indicates the number of validation sets; Represents the sequence data set used for verification, that is, the true value; Representation sequence The global prediction value of Representation sequence The local rolling prediction value; k1 represents a hyperparameter; i represents a variable; L input Indicates the length of the input sequence;

[0111] If all variables satisfy the above formula, it means that the global forecast is effective enough and there is no need for local rolling forecast. If some variables do not satisfy the above formula, it means that the global forecast performs poorly on these variables and local rolling forecast is needed.

[0112] Step (4.8): Local prediction: Determine whether to adopt local rolling prediction based on the verification recognition result of global prediction; before local prediction, perform a simple preprocessing process on the input data, perform weighted processing on each window of the test set, and construct a local recognition curve x * , which is used to capture the global trend and local fluctuations of the data; its formula is as follows:

[0113]

[0114] In the formula, x represents the original input data; y represents the corresponding result obtained by global prediction; * Indicates the construction of a new input curve; ω represents the weight vector, which is used to control the proportion of the original input data x and the global prediction result y in the weighted average. Represents element-wise multiplication;

[0115] ω=(linspace(α 1 / γ ,β 1 / γ ,L batch )) γ

[0116] In the formula, the weight ω uses the linear interpolation function linspace to generate a linearly changing weight vector whose length is equal to L batch ; α and β represent the global and local weight coefficients respectively; γ represents the acceleration of weight change;

[0117] Weighted input curve x * Perform Fourier transform to extract high-frequency components and convert them into potential cycles Next, use the cycle To construct the local basis function Θ l , the formula is as follows:

[0118] Θ l =[sin(C3t),cos(C3t)]

[0119] In the formula, Indicates potential cycle The coefficient of

[0120] The global basis function Θ g , store basis function Θ s and the local basis function Θ l Fusion is performed to obtain the final basis function library Θ pred , the formula is as follows:

[0121] Θ pred =[Θ g ,Θ s ,Θ l ];

[0122] Solving the local prediction parameter Ξ through sparse identification pred Implement local prediction and its sparse recognition optimization problem; set the prediction time to t = L input +1,L input +2,…, using Ξ pred and θ pred Perform local rolling prediction to obtain the final local prediction result

[0123] Step (4.9): Output the final prediction result; when all variables meet the verification and identification conditions of the global prediction, the final prediction result is the global prediction result When some variables do not meet the verification and identification conditions of the global prediction, the final prediction result is the local prediction result. Finally, the prediction result of the high-frequency component F1(t) is obtained as

[0124] Furthermore, in step (5), the implementation process of the BKA optimization algorithm to optimize the prediction of the intermediate frequency and low frequency components of the ESN prediction model is as follows:

[0125] Step (5.1): Use the BKA optimization algorithm to optimize the parameters of the ESN prediction model and find the optimal value;

[0126] The parameters include the reserve pool spectrum radius SR, the reserve pool sparsity degree SD and the input unit scale IS, which are substituted into the ESN prediction model to predict the intermediate frequency component F2(t) and the low frequency component F3(t) respectively;

[0127] Step (5.2): Assume that the ESN prediction model has K input nodes, N reserve pool nodes, and L output nodes; split the time series data set of the intermediate frequency component F2(t) or the low frequency component F3(t) into a training set and a test set in a ratio of 7:3, and set them as the input data u(t). Then, the update formulas of the reserve pool state x(t) and the output result y(t) are as follows:

[0128] x(t+1)=(1-α)x(t)+ατ[W in u(t+1)+Wx(t)+W back y(t)]

[0129] y(t+1)=σ[W out (u(t+1),x(t+1),y(t))]

[0130] Where x(t) and x(t+1) represent the state of the reserve pool at time t and time t+1; y(t) and y(t+1) represent the output state of the reserve pool at time t and time t+1; u(t+1) represents the input data accepted by the ESN model at time t+1; α represents the leakage rate, that is, the speed of the reserve pool update; W in , W, W back represents the input weight matrix, the reserve pool weight matrix, and the output feedback weight matrix; W out represents the output weight matrix; τ and σ represent the activation function;

[0131] The output weight W is solved by regularized ridge regression method. out , the formula is as follows:

[0132] Y target =W out X

[0133] W out =Y target W T (XX T +βI) -1

[0134] Where Y target represents the target value; X represents the internal state matrix of the reserve pool, X Trepresents the transpose of the internal state matrix of the reservoir; β represents the regularization coefficient; I represents the identity matrix; W T represents the transpose of the weight matrix;

[0135] Step (5.3): Input the test set data into the trained ESN prediction model, through the reserve pool state x(t) and output weight W out The prediction result is calculated as follows:

[0136] y'(t)=tanh(W out x(t)

[0137] Where y'(t) represents the predicted output;

[0138] Finally, the prediction results of the mid-frequency component F2(t) and the low-frequency component F3(t) are obtained. and

[0139] Furthermore, in step (6), the implementation process of superimposing and reconstructing the prediction results is as follows:

[0140] Prediction results for the obtained high-frequency component F1(t) Prediction results of mid-frequency component F2(t) and low-frequency component F3(t) and The weighted sum composite prediction model is used to integrate the prediction results of lithium batteries, and the formula is:

[0141]

[0142] Where Y represents the predicted value of the composite prediction model for the RUL of lithium batteries; α1, α2, and α3 represent random fusion coefficients in the range [0, 1]. Generally, α1∈[0.9, 1], α2∈[0.6, 0.9], and α3∈[0, 0.5] are set.

[0143] The beneficial effects of the present invention are as follows: 1. The present invention uses the TVFEMD model and the EWT model to perform multi-level decomposition on the battery capacity data, so as to more accurately extract the key signal components reflecting the battery performance degradation; this refined signal processing method helps to capture the slight changes of the battery in different use stages, thereby improving the accuracy of the remaining service life prediction; 2. By combining K-means 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 two prediction models, BKA-GLIP and BKA-ESN, which are respectively suitable for processing signal components of different frequencies; this model-based 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

[0144] Figure 1 is the overall system flow chart of the present invention;

[0145] Figure 2 is a flow chart of the BKA-GLIP prediction model of the present invention;

[0146] Figure 3 It is a flow chart of the BKA-ESN prediction model of the present invention. DETAILED DESCRIPTION

[0147] The specific technical scheme of the present invention is further described in detail below with reference to specific examples.

[0148] As shown in the figure, the method for predicting the remaining life of a lithium battery according to the present invention has the following operating steps:

[0149] Step 1, collect the operating data of the lithium battery, including discharge voltage, discharge current, discharge time and capacity sequence information; use the TVFEMD model to decompose the lithium battery capacity sequence, which is divided into high frequency, medium frequency, low frequency and residual components from high to low frequency;

[0150] Step 2: Use the EWT model to perform secondary decomposition on the high-frequency and residual components in the IMF signal. The IMF components after secondary decomposition are high-frequency, medium-frequency, and low-frequency components according to the frequency from high to low.

[0151] Step 3: Calculate the sample entropy of the intermediate frequency and low frequency components after the first decomposition and the high frequency, intermediate frequency and low frequency components after the second decomposition. According to the sample entropy evaluation results, classify the components after the second decomposition into high frequency, intermediate frequency and low frequency component sets. Use the K-means clustering method to perform cluster analysis on the classified component sets to obtain three signal components of high frequency, intermediate frequency and low frequency.

[0152] Step 4: The high-frequency components after K-means clustering are predicted using the global local information model (GLIP), and the mid-frequency and low-frequency components are predicted using the echo state network (ESN). In order to improve the prediction performance of the two prediction models, the black kite optimization algorithm (BKA) is used to optimize the hyperparameters α, β, and γ in the GLIP model and the hyperparameters SR, SD, and IS in the ESN model;

[0153] Step 5: Reconstruct the prediction results of the two models to obtain the prediction results of the remaining service life of the lithium battery. The specific flow chart is as follows: Figure 1 shown.

[0154] The specific implementation process is:

[0155] In step (1), the specific implementation process of using the TVFEMD model to decompose the lithium battery capacity sequence is:

[0156] S1.1. In the present invention, the pre-processed capacity data is used as the input signal, the original battery capacity signal is set to x(t), and the signal after the Hilbert transform is obtained from the original battery capacity signal x(t) is x'(t). The specific formula is as follows:

[0157]

[0158]

[0159] Wherein, X(t) represents the instantaneous amplitude of the battery capacity signal x(t); φ(t) represents the instantaneous phase of the battery capacity signal x(t); x'(t) represents the Hilbert transform signal of the battery capacity signal x(t);

[0160] The analytical signal corresponding to the battery capacity signal x(t) can be expressed as:

[0161] h(t)=X(t)e jφ(t) =α1e jφ(t) +α2e jφ(t)

[0162] Where h(t) represents the analytical signal of signal x(t); α1 and α2 represent the instantaneous amplitude of the component signal. The phase difference represents the instantaneous amplitude;

[0163] Furthermore, the instantaneous amplitude X(t) can be expressed as:

[0164] X 2 (t) = α1 2 (t)+α2 2 (t)+2α1(t)α2(t)cos[φ1(t)-φ2(t)]

[0165] Where φ1(t) and φ2(t) represent the instantaneous phase of the component signal, and α1(t) and α2(t) represent the instantaneous amplitude of the signal component at the tth moment;

[0166] S1.2. Determine the local maximum value X(t) of the instantaneous amplitude X(t) max ) sequence and the local minimum sequence X(t min ), and interpolate the local maximum curve β1(t) and the local minimum curve β2(t) respectively, and then calculate the amplitudes α1(t) and α2(t), the instantaneous frequency components φ1'(t) and φ2'(t). The specific formula is as follows:

[0167]

[0168] Where β1(t) and β2(t) represent the local maximum and local minimum; α1(t) and α2(t) represent the new amplitude;

[0169] The instantaneous frequency component φ'(t) of the battery capacity signal x(t) is:

[0170]

[0171] X 2 (t max )φ'(t max ) and X 2 (t min )φ'(t min ) is interpolated to obtain η1(t) and η2(t), then the instantaneous frequency components φ1'(t) and φ2'(t) are:

[0172]

[0173] In the formula, η1(t) and η2(t) represent the 2 (t max )φ'(t max ) and X 2 (t min )φ'(t min ) is obtained by interpolation;

[0174] Furthermore, the local cutoff frequency φ of the signal is calculatedb '(t), the specific formula is as follows:

[0175]

[0176] S1.3. Filter the uncalculated signal using a time-varying filter to obtain the local mean and calculate the signal z(t); as follows:

[0177] z(t)=cos[∫φ' b (t)dt]

[0178] In the formula, φ b '(t) is the local cutoff frequency of the signal; when z(t) is at an extreme value, the time point t is extracted max and t min As node m, a B-spline approximation time-varying filter can be constructed to perform B-spline approximation time-varying filtering on the signal x(t). The approximation result is recorded as m(t), which represents the local mean.

[0179] Continuously update the signal x(t) and define the stopping standard θ(t). When the signal satisfies θ(t)≤ε, ε is the given bandwidth threshold, the signal component at this time is recorded as IMF i component signal; otherwise, let x(t) = x(t) - m(t), and repeat the above process to continue decomposing the signal until the stopping criterion is met and the decomposition is completed;

[0180] The specific formula of the stopping criterion θ(t) is as follows:

[0181]

[0182] Where B(t) represents the instantaneous bandwidth of the signal; φ avg (t) represents the weighted mean instantaneous frequency;

[0183] From the above process, the TVFEMD model decomposes the lithium battery capacity sequence once, and the original signal x(t) is decomposed into a series of IMF signals, which are high frequency, medium frequency, and low frequency three signal components IMF1, IMF2, and IMF3 according to the frequency from high to low;

[0184] S1.4. When all IMF signals are extracted, the remaining signal is the residual component R(t); that is, the following formula:

[0185]

[0186] In step (2), the high-frequency signal and residual signal obtained by the TVFEMD algorithm in S1 are further decomposed by the EWT algorithm to further extract the detailed features of the high-frequency signal and the residual signal. The specific implementation process of the EWT algorithm is as follows:

[0187] S2.1. Perform Fourier Transform (FT) on the signal IMF1 retained in S1 and the residual signal R(t) to convert the time domain signal into a frequency domain signal. The specific formula of Fourier Transform is as follows:

[0188]

[0189] S2.2, construct the scaling function and wavelet function required for the EWT algorithm; through the empirical scaling function With the empirical wavelet function ψ n (ω) realizes the adaptive properties of (H-1) bandpass filters, where the empirical scaling function With the empirical wavelet function ψ n The calculation formula of (ω) is as follows:

[0190]

[0191]

[0192] In the formula, represents the empirical scaling function; ψ n (ω) represents the empirical wavelet function; n represents the number of decomposition levels; ω represents the frequency, ω n represents the frequency of the nth boundary; β(·) represents the transformation function; τ n represents the coefficient function, τ n =γω n , (0<γ<1);

[0193] S2.3. Based on the empirical scaling function The inner product of the empirical wavelet function and the original signal is used to obtain the scale coefficient W(0,t), and the inner product of the empirical wavelet function and the original signal is used to obtain the empirical wavelet coefficient W(k,t). This coefficient contains the characteristic information of the high-frequency signal IMF1 and the residual signal IMF4 in different frequency bands. The specific formula is as follows:

[0194]

[0195] In the formula, t represents the time point; and represents the empirical scaling function and the empirical wavelet function; and represents the Fourier transform of the empirical scaling function and the empirical wavelet function; F -1[·] denotes inverse Fourier transform;

[0196] S2.4, reconstruct the signal IMF1 and the residual signal R(t) according to the scale coefficient W(0,t) and the empirical wavelet coefficient W(k,t);

[0197] The signal IMF1 and the residual signal R(t) are decomposed and reconstructed by the EWT algorithm as follows:

[0198]

[0199] In the formula, k represents the number of decomposed components, which is set to k = 3 here;

[0200] Furthermore, through the EWT algorithm, the signal IMF1 and the residual signal R(t) are decomposed into multiple components, each component corresponds to a different frequency band; IMf4, IMf5, and IMf6 correspond to the high-frequency signal decomposition results of IMF1, and IMf7, IMf8, and IMf9 correspond to the decomposition results of R(t).

[0201] In step (3), for the eight signal components decomposed by the TVFEMD-EWT algorithm, the sample entropy of each component IMF is calculated to measure the signal complexity, and K-means clustering is used for signal classification. The specific implementation process is as follows:

[0202] S3.1. First determine the subsequence length, assuming that each IMF component is represented by IMF n(t) , n∈1,2,…8, each IMF component IMF n(t) Divide into m-dimensional continuous subsequences,

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

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

[0205] Where N represents IMF i(t) Length;

[0206] For each subsequence x i , calculate other subsequences x j (where k≠j) the Euclidean distance d(x i ,x j ); define the sequence x i With 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:

[0207]

[0208] S3.2. Calculate the similarity B between subsequences i m (r); given a threshold r, statistics d[x i ,x j ]<r and compared with the total number of vectors Nm, the specific formula is as follows:

[0209]

[0210] In the formula, num{d[x i ,x j ]<r} means d max (x i ,x j )<the number of r;

[0211] Further, the sample entropy of each subsequence is calculated, and the i m The results of (r) are averaged and logarithmically transformed to obtain the sample entropy SE(i,r). The specific formula is as follows:

[0212]

[0213] SE m (i,r)=-ln[(B m (r))]

[0214] In the formula, represents the similarity of the i-th subsequence; B m (r) represents the average similarity of m-dimensional vectors; SE m (i,r) represents the sample entropy of the i-th subsequence m-dimensional;

[0215] Further, we increase the dimension and add 1 to the dimension m, that is, we get the similarity B for m+1 vectors. m+1 (r) and sample entropy SE m+1 (i,r);

[0216] S3.3, when m reaches a certain preset value or the sample entropy does not change much, stop the iteration; finally, the sample entropies SE1, SE2, ..., SE8 of the calculated 8 IMF components are combined into a sample entropy vector S i , S i ={SE1,SE2,…,SE8};

[0217] According to the sample entropy vector S calculated above i , select three initial center cluster numbers C1, C2, C3, representing high-frequency, medium-frequency and low-frequency signals respectively; for each sample entropy S i , calculate its relationship with each center C j The distance is usually calculated using the Euclidean distance, and the specific formula is as follows:

[0218]

[0219] In the formula, S i represents the i-th sample entropy vector; C j represents the number of the jth cluster, j∈[1,3];

[0220] For each sample entropy value S i , find the center C that is closest to it j , and assign it to this center. The specific formula is as follows:

[0221] C j* =argmin j d(S i ,C j )

[0222] In the formula, C j* Indicates distance S i nearest center;

[0223] S3.4. Recalculate the center of each cluster and take the average value of the entropy of all samples assigned to the cluster as the new cluster center. The specific formula is as follows:

[0224]

[0225] In the formula, |C j | indicates assignment to cluster C j The number of samples;

[0226] The iterative clustering process is repeated until the new cluster center is basically unchanged. Finally, each IMF component will be assigned to the corresponding cluster, and the clustering results are output: high-frequency component F1(t), medium-frequency component F2(t) and low-frequency component F3(t).

[0227] In step (4): according to the signal classification result of the above K-means clustering, the high-frequency component F1(t) is predicted using the BKA-GLIP prediction model; the high-frequency component F1(t) is used as the input sequence, and its data sequence is set to X; the specific implementation process is as follows:

[0228] S4.1. Initialization of the Black Kite (BKA) optimization algorithm: A group of candidate solutions are generated by random initialization. Each individual solution represents one of the parameters of the GLIP model, including hyperparameters α, β, and γ. α and β represent the global and local weight coefficients, respectively. γ represents the acceleration of weight change.

[0229] Assuming that the population size of the BKA optimization algorithm is N, the maximum number of iterations is T, and the variable dimension is j, the specific formula is as follows:

[0230] X i =BK lb +rand(BK ub -BK lb )

[0231] In the formula, X i represents individual solution, i∈[1,N]; BK lb , BK ub They represent the upper and lower bounds of the range of the j-th dimension individual solution respectively; rand represents a random number between [0,1]; at the same time, the f(·) function is introduced to represent the fitness value of the individual;

[0232] Furthermore, during the initialization process, BKA selects the individual with the best fitness value as the optimal solution X in the initial population. L , that is, the optimal solution is set as the hyperparameter of the GLIP prediction model; the specific formula is as follows:

[0233] f best =Best(f(X i ))

[0234] X L =X(find(f best = = f(X i )))

[0235] In the formula, f best represents the optimal individual fitness value; X L represents the optimal solution in the initial population;

[0236] S4.2, BKA position update, iteratively optimize the hyperparameters of the GLIP prediction model; in each iteration, the position of the candidate solution is updated according to the fitness function to find a better hyperparameter setting for GLIP;

[0237] BKA position update, the specific formula is as follows:

[0238]

[0239] In the formula, and Respectively represent the position solutions of the ith individual in the jth dimension at the tth iteration and the t+1th iteration; r is a random number between [0,1]; p is a constant; t is the number of iterations completed so far; T is the total number of iterations; n is the adhesion coefficient;

[0240] BKA continues to update the position, the specific formula is as follows:

[0241]

[0242] m=2×sin(r+π / 2)

[0243] In the formula, represents the optimal solution of the j-th dimension individual in the t-th iteration; F i represents the current position of any individual in the jth dimension at the tth iteration; F ri represents the fitness value of any individual at the random position of the jth dimension at the tth iteration; C(0,1) represents the Cauchy mutation;

[0244] S4.3. Based on the obtained best individual fitness value, BKA selects the best individual and saves it, that is, updates the hyperparameters of the GLIP prediction model;

[0245]

[0246] In the formula, X best represents the individual optimal solution; F best Represents the optimal fitness value of an individual; represents the optimal solution of the j-th dimension individual in the t-th iteration; and They represent the position solutions of the ith individual in the jth dimension at the tth iteration and the t+1th iteration respectively;

[0247] S4.4, determine whether the iterative optimization meets the maximum iteration condition, if it does, the iteration stops, otherwise it continues to loop until the optimal solution is obtained, that is, the three hyperparameters α, β, and γ of the GLIP prediction model;

[0248] S4.5, the three hyperparameters α, β, and γ of the GLIP prediction model obtained by the BKA optimization algorithm are substituted into the GLIP prediction model to predict the high-frequency component F1(t);

[0249] Furthermore, the input high-frequency component F1(t) is determined, its data sequence is set to X, the input length is L, the output length is O, and the hyperparameters α, β, and γ obtained after the above optimization are set; the high-frequency component data set is divided into the training set X and the training set X in a ratio of 7:1:2. train , validation set X val and the test set X test ;

[0250] S4.6, GLIP global prediction:

[0251] First, for the entire training set X train Perform a discrete Fourier transform (DFT) to obtain the amplitude A and frequency f; then, select the high amplitude frequency f from the amplitude * , then its corresponding period T * =1 / f * , let T * =[T1, T2, ...]; then, using the period T * To construct the global basis function Θ g , the specific formula is as follows:

[0252] Θ g =[sin(C1t),cos(C1t),1]

[0253] Where, C1 = 2π / T * , t represents the length of the training set, T * Indicates high frequency period;

[0254] Furthermore, the sparse identification process can be obtained according to the global basis function, and the global parameter Ξ can be solved by sparse identification. g , to achieve global prediction, the specific formula for the optimization problem of sparse recognition is as follows:

[0255]

[0256] Where X represents the training set data sequence; ||X-Θ g Ξ g || 2 Represents the global basis function for the training set X train The fitting error of g || indicates the global parameter Ξ g l1 norm regularization of is used to achieve sparsity;

[0257] Furthermore, by using the coordinate descent method similar to LASSO to solve the above optimization problem, we get Θ g , let t = L train +1,L train +2,… to predict the validation set X val and the test set X test The data in the prediction result is recorded as X gp , that is, to achieve global prediction;

[0258] Furthermore, using the global prediction result X gp , construct a new basis function The specific formula is as follows:

[0259]

[0260] Among them, X g =[X train ,X gp ], indicating the univariate prediction result X gp With training set X train splicing;

[0261] The new basis function Alternative Θ g , re-solve the optimization problem of sparse identification to obtain the global prediction result considering the network coupling relationship between variables

[0262] S4.7. Construct storage basis functions: Divide the test set into multiple batches, each batch contains input sequences and output sequences of a certain length; perform discrete Fourier transform (DFT) on each variable in each batch to obtain its frequency domain features, extract high-frequency components from them, and store them in a set; convert the high-frequency components extracted from all batches into corresponding period T s * , using the extracted periodic set, construct the storage basis function Θ s ;

[0263] The specific formula is as follows:

[0264] Θ s =[sin(C2t),cos(C2t)]

[0265] Where, C2 = 2π / T s * , T s * Indicates the period corresponding to the high-frequency components extracted from the test set;

[0266] S4.8, verify and identify whether to perform local rolling prediction;

[0267] Furthermore, the GLIP model is used for verification and recognition to evaluate the effectiveness of global recognition and prediction, to observe whether the use of local features will perform better than global prediction, and to prepare for subsequent local rolling predictions;

[0268] Assume that the global prediction value for the validation set data is Calculate the forecast error for each variable i and compare it with the local rolling forecast The specific implementation formula is as follows:

[0269]

[0270] Where N bIndicates the number of validation sets; Represents the sequence data set used for verification, that is, the true value; Representation sequence The global prediction value of Representation sequence The local rolling prediction value; k1 represents a hyperparameter; i represents a variable; L input Indicates the length of the input sequence;

[0271] If all variables satisfy the above formula, it means that the global forecast is effective enough and there is no need for local rolling forecast. If some variables do not satisfy the above formula, it means that the global forecast performs poorly on these variables and local rolling forecast is needed.

[0272] S4.9, local prediction: Determine whether to adopt local rolling prediction based on the verification recognition result of global prediction; before local prediction, perform a simple preprocessing process on the input data, perform weighted processing on each window of the test set, and construct a local recognition curve x * , which is used to capture the global trend and local fluctuation of the data; the specific formula is as follows:

[0273]

[0274] In the formula, x represents the original input data; y represents the corresponding result obtained by global prediction; * Indicates the construction of a new input curve for local recognition and prediction; ω represents the weight vector, which is used to control the proportion of the original input data x and the global prediction result y in the weighted average. Represents element-wise multiplication;

[0275] ω=(linspace(α 1 / γ ,β 1 / γ ,L batch )) γ

[0276] In the formula, the weight ω uses the linear interpolation function linspace to generate a linearly changing weight vector whose length is equal to L batch ; Among them, α and β represent the global and local weight coefficients respectively; γ represents the acceleration of weight change;

[0277] Furthermore, for the weighted input curve x * Perform Fourier transform (DFT) to extract the high-frequency components and convert them into potential cycles Next, use the cycle To construct the local basis function Θ l , the specific formula is as follows:

[0278] Θ l=[sin(C3t),cos(C3t)]

[0279] In the formula, Indicates potential cycle The coefficient of

[0280] Furthermore, the global basis function θ g , store basis function Θ s and the local basis function Θ l Fusion is performed to obtain the final basis function library Θ pred , the specific formula is as follows:

[0281] Θ pred =[Θ g ,Θ s ,Θ l ]

[0282] Furthermore, the local prediction parameter Ξ is solved by sparse identification pred , realize local prediction, and its sparse recognition optimization problem; set the prediction time to t = L input +1,L input +2,…, using Ξ pred and θ pred Perform local rolling prediction to obtain the final local prediction result

[0283] S4.10, output the final prediction result:

[0284] When all variables meet the verification and identification conditions of the global prediction, the final prediction result is the global prediction result. When some variables do not meet the verification and identification conditions of the global prediction, the final prediction result is the local prediction result. Finally, the prediction result of the high-frequency component F1(t) is obtained as

[0285] In addition, the BKA-ESN prediction model is used to predict the intermediate frequency component F2(t) and the low frequency component F3(t), wherein the specific implementation process of the BKA-ESN prediction model is as follows:

[0286] Step 1: Optimize the parameters of the ESN prediction model according to the BKA optimization algorithm in S4 above to find the optimal values. These parameters include the reserve pool spectrum radius SR, the reserve pool sparsity SD and the input unit scale IS. Substitute them into the ESN model to predict the mid-frequency component F2(t) and the low-frequency component F3(t) respectively.

[0287] Step 2: Set the ESN to have K input nodes, N reserve nodes, and L output nodes; split the time series data set of the intermediate frequency component F2(t) (or low frequency component F3(t)) into a training set and a test set in a ratio of 7:3, and set them as the input data u(t). The specific formulas for updating the reserve state x(t) and the output result y(t) are as follows:

[0288] x(t+1)=(1-α)x(t)+ατ[W in u(t+1)+Wx(t)+W back y(t)]

[0289] y(t+1)=σ[W out (u(t+1),x(t+1),y(t))]

[0290] Where x(t) and x(t+1) represent the state of the reserve pool at time t and time t+1; y(t) and y(t+1) represent the output state of the reserve pool at time t and time t+1; u(t+1) represents the input data accepted by the ESN model at time t+1; α represents the leakage rate, that is, the speed of the reserve pool update; W in , W, W back represents the input weight matrix, the reserve pool weight matrix, and the output feedback weight matrix; W out represents the output weight matrix; τ and σ represent the activation function;

[0291] The output weight W is solved by regularized ridge regression method. out , the specific formula is as follows:

[0292] Y target =W out X

[0293] W out =Y target W T (XX T +βI) -1

[0294] Where Y target represents the target value; X represents the internal state matrix of the reserve pool, X T represents the transpose of the internal state matrix of the reservoir; β represents the regularization coefficient; I represents the identity matrix; W T represents the transpose of the weight matrix;

[0295] Furthermore, the test set data is input into the trained ESN prediction model, and the pool state x(t) and output weight W are used to calculate the prediction result. out Calculate the prediction results. The specific formula is as follows:

[0296] y'(t)=tanh(W out x(t)

[0297] Where y'(t) represents the predicted output;

[0298] Finally, the prediction results of the mid-frequency component F2(t) and the low-frequency component F3(t) are obtained. and

[0299] In step (6), the prediction result of the high frequency component F1(t) obtained by the BKA-GLIP prediction model is Prediction results of mid-frequency component F2(t) and low-frequency component F3(t) obtained using the BKA-ESN prediction model and Furthermore, a weighted sum composite prediction model is used to integrate the prediction results of the lithium battery obtained above, and its specific formula is:

[0300]

[0301] Where Y represents the predicted value of the composite prediction model for the RUL of lithium batteries; α1, α2, and α3 represent random fusion coefficients in the range [0, 1]. Generally, α1∈[0.9, 1], α2∈[0.6, 0.9], and α3∈[0, 0.6] are set.

[0302] In order to evaluate the training effect of the prediction model, error analysis is required; the following error parameter evaluation indicators are used: root mean square error RMSE, mean absolute error MAE, and mean absolute percentage error MAPE are used to verify the performance of the model; the calculation formulas of the corresponding indicators are as follows:

[0303] Root mean square error formula:

[0304]

[0305] Mean absolute error formula:

[0306]

[0307] Mean absolute percentage error formula:

[0308]

[0309] In the formula, represents the predicted value; y i represents the real value; n represents the quantity;

[0310] Furthermore, in order to ensure the accuracy of the RUL prediction of lithium batteries, the remaining operating data of the collected lithium batteries, including voltage, current, charge and discharge time, temperature and other data information, are brought into the BKA-GLIP prediction model for further training and prediction to improve the RUL prediction results of lithium batteries.

Claims

1. A method for predicting the remaining life of a lithium battery, characterized in that: The operation steps are as follows: Step (1): Collecting the operating data of the lithium battery; using the TVFEMD model to decompose the lithium battery capacity sequence, which is divided into high frequency, medium frequency, low frequency and residual components from high to low frequency; Step (2): Use the EWT model to perform secondary decomposition on the high-frequency and residual components in the IMF signal. The IMF components after secondary decomposition are high-frequency, medium-frequency, and low-frequency components according to the frequency from high to low. Step (3): Calculate sample entropy for the intermediate frequency and low frequency components after the first decomposition and the high frequency, intermediate frequency and low frequency components after the second decomposition. Classify the components after the second decomposition into high frequency, intermediate frequency and low frequency component sets according to the sample entropy evaluation results. Perform cluster analysis on the classified component sets using the K-means clustering method to obtain three signal components of high frequency, intermediate frequency and low frequency. Step (4): The high-frequency components after K-means clustering are predicted using the GLIP prediction model, and the medium-frequency and low-frequency components are predicted using the ESN network; Step (5): In order to improve the prediction performance of the two prediction models, the BKA optimization algorithm is used to optimize the hyperparameters α, β, and γ in the GLIP prediction model and the hyperparameters of the reserve pool spectrum radius SR, the reserve pool sparsity SD, and the input unit scale IS in the ESN network; Step (6): Finally, the prediction results of the two models are reconstructed to obtain the remaining service life prediction results of the lithium battery.

2. The method for predicting the remaining life of a lithium battery according to claim 1, characterized in that: In step (1), the operating data includes information on discharge voltage, discharge current, discharge time and capacity sequence.

3. The method for predicting the remaining life of a lithium battery according to claim 1, characterized in that: In step (1), the TVFEMD model is used to decompose the lithium battery capacity sequence. The implementation process of the TVFEMD model to decompose the lithium battery capacity sequence is as follows: Step (1.1): Set the original battery capacity signal to x(t), and use Hilbert transform to obtain the signal x'(t) after the original battery capacity signal x(t) is transformed. The formula is as follows: Wherein, X(t) represents the instantaneous amplitude of the battery capacity signal x(t); φ(t) represents the instantaneous phase of the battery capacity signal x(t); x'(t) represents the Hilbert transform signal of the battery capacity signal x(t); The formula of the analytical signal h(t) corresponding to the battery capacity signal x(t) is expressed as: h(t)=X(t)e jφ(t) =α1e jφ(t) +α2e jφ(t) In the formula, h(t) represents the analytical signal of signal x(t), α1 and α2 represent the instantaneous amplitude of the component signal, The phase difference represents the instantaneous amplitude; The formula for instantaneous amplitude X(t) is expressed as: X 2 (t)=α1 2 (t)+α2 2 (t)+2α1(t)α2(t)cos[φ1(t)-φ2(t)] Where φ1(t) and φ2(t) represent the instantaneous phase of the component signal, and α1(t) and α2(t) represent the instantaneous amplitude of the signal component at the tth moment; Step (1.2): Determine the local maximum value X(t) of the instantaneous amplitude X(t) max ) sequence and the local minimum sequence X(t min ), and interpolate the local maximum curve β1(t) and the local minimum curve β2(t) respectively, and then calculate the amplitudes α1(t) and α2(t), the instantaneous frequency components φ1'(t) and φ2'(t), the formula is as follows: Where β1(t) and β2(t) represent the local maximum and local minimum; α1(t) and α2(t) represent the new amplitude; The formula for the instantaneous frequency component φ'(t) of the battery capacity signal x(t) is: X 2 (t max )φ'(t max ) and X 2 (t min )φ'(t min ) is interpolated to obtain η1(t) and η2(t), then the formulas for the instantaneous frequency components φ1'(t) and φ2'(t) are: In the formula, η1(t) and η2(t) represent the 2 (t max )φ'(t max ) and X 2 (t min )φ'(t min ) is obtained by interpolation; Calculate the local cutoff frequency φ of the signal b '(t), the formula is as follows: Step (1.3): Filter the uncalculated signal using a time-varying filter to obtain the local mean, and calculate the signal z(t); as shown in the following formula: z(t)=cos[∫φ' b (t)dt] In the formula, φ b '(t) is the local cutoff frequency of the signal; when z(t) is at an extreme value, the time point t is extracted max and t min As node m, a B-spline approximation time-varying filter is constructed, and the signal x(t) is subjected to B-spline approximation time-varying filtering. The approximation result is recorded as m(t), which represents the local mean; Step (1.4): When all IMF signals are extracted, the remaining signal is the residual component R(t); as shown in the following formula:

4. The method for predicting the remaining life of a lithium battery according to claim 3, characterized in that: In step (1.3), the signal x(t) is continuously updated, and the stopping criterion θ(t) is defined. When the signal satisfies θ(t)≤ε, ε is a given bandwidth threshold, the signal component at this time is recorded as IMF i component signal; otherwise, let x(t) = x(t) - m(t), and repeat the above steps (1.1)-(1.3) to continue decomposing the signal until the stopping criterion is met and the decomposition is completed; Wherein, the formula of the stopping criterion θ(t) is as follows: Where B(t) represents the instantaneous bandwidth of the signal; φ avg (t) represents the weighted mean instantaneous frequency; According to steps (1.1) to (1.3), the TVFEMD model decomposes the lithium battery capacity sequence and extracts a series of IMF signals, which are high frequency, medium frequency, and low frequency signal components IMF1, IMF2, and IMF3, respectively, from high to low frequency.

5. The method for predicting the remaining life of a lithium battery according to claim 1, characterized in that: In step (2), the implementation process of using the EWT model to perform secondary decomposition on the high frequency and residual components in the IMF signal obtained after the primary decomposition is: Step (2.1): Perform Fourier transform on the signal IMF1 retained in step (1) and the residual signal R(t); Step (2.2): Through the empirical scaling function With the empirical wavelet function ψ n (ω) realize the adaptive properties of H-1 bandpass filters; Among them, the empirical scaling function With the empirical wavelet function ψ n The calculation formula of (ω) is as follows: In the formula, represents the empirical scaling function; ψ n (ω) represents the empirical wavelet function; ω represents the frequency, ω n represents the frequency of the nth boundary; β(·) represents the transformation function; τ n represents the coefficient function, τ n =γω n , (0<γ<1); Step (2.3): Based on the empirical scaling function The inner product of the empirical wavelet function and the original signal is used to obtain the scaling coefficient, and the inner product of the empirical wavelet function and the original signal is used to obtain the empirical wavelet coefficient. The formula is as follows: In the formula, t represents the time point; and represents the empirical scaling function and the empirical wavelet function; and represents the Fourier transform of the empirical scaling function and the empirical wavelet function; F -1 [·] denotes inverse Fourier transform; Step (2.4): The signal IMF1 and the residual signal R(t) are decomposed and reconstructed by the EWT algorithm as follows: In the formula, k represents the number of decomposed components, and its value is k=3.

6. The method for predicting the remaining life of a lithium battery according to claim 1, characterized in that: In step (3), the implementation process of cluster analysis using the K-means clustering method is as follows: Step (3.1): Let each IMF component be IMF n(t) , n∈1,2,…8, each IMF component IMF n(t) Divide into m-dimensional continuous subsequences as follows: 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 i(t) Length; For each subsequence x i , calculate other subsequences x j The Euclidean distance d(x i ,x j ); define the sequence x i With the sequence x j The maximum distance between corresponding elements is d max (x i ,x j ); its formula is as follows: Step (3.2): Given a threshold r, count d[x i ,x j ]<r and is ratioed to the total number of vectors Nm, recorded as The formula is as follows: In the formula, num{d[x i ,x j ]<r} means d max (x i ,x j )<the number of r; right The average of the results is calculated using the following formula: THIS m (i,r)n-ln[(B m (r))] In the formula, represents the similarity of the i-th subsequence; B m (r) represents the average similarity of m-dimensional vectors; SE m (i, r) represents the sample entropy of the m-dimensional i-th subsequence, and then the dimension m is added by 1, that is, the similarity B is obtained for m+1 vectors m+1 (r) and sample entropy SE m+1 (i,r); Step (3.3): The sample entropies of the eight IMF components SE1, SE2, …, SE8 are combined into a sample entropy vector S i , S i ={SE1, SE2, …, SE8}; select 3 initial center cluster numbers C1, C2, C3, representing high-frequency, medium-frequency and low-frequency signals respectively; for each sample entropy S i , calculate its relationship with each center C j The distance is usually calculated using the Euclidean distance, which is as follows: In the formula, S i represents the i-th sample entropy vector; C j represents the number of the jth cluster, j∈[1,3]; For each sample entropy value S i , find the center C that is closest to it j , and assign it to this center, the formula is as follows: In the formula, Indicates distance S i nearest center; Step (3.4): Recalculate the center of each cluster and take the average value of the entropy of all samples assigned to the cluster as the new cluster center. The formula is as follows: In the formula, |C j | indicates assignment to cluster C j The number of samples; Repeat iterative steps (3.3) to (3.4) until the new cluster center is basically unchanged. Finally, each IMF component will be assigned to the corresponding cluster, and the clustering results are output: high-frequency component F1(t), medium-frequency component F2(t) and low-frequency component F3(t).

7. The method for predicting the remaining life of a lithium battery according to claim 1, characterized in that: In step (4), the implementation process of using the GLIP prediction model to predict the high frequency component is: Step (4.1): Initialization of the BKA optimization algorithm: Randomly initialize to generate a group of candidate solutions, each individual solution represents one of the parameters of the GLIP model, including hyperparameters α, β, and γ; Among them, α and β represent the global and local weight coefficients respectively; γ represents the acceleration of weight change; the population size of the BKA optimization algorithm is set to N, the maximum number of iterations is T, and the variable dimension is j; at the same time, the f(·) function is introduced to represent the fitness value of the individual. During the initialization process, BKA selects the individual with the best fitness value as the optimal solution X in the initial population. L , that is, setting the optimal solution as the hyperparameter of the GLIP prediction model; Step (4.2): Update the BKA position and iteratively optimize the hyperparameters of the GLIP prediction model; in each iteration, update the position of the candidate solution according to the fitness function to find the hyperparameter settings of the GLIP prediction model; The formula for updating the BKA position is as follows: In the formula, and Respectively represent the position solutions of the ith individual in the jth dimension at the tth iteration and the t+1th iteration; r is a random number between [0,1]; p is a constant; t is the number of iterations completed so far; T is the total number of iterations; n is the adhesion coefficient; The formula for BKA to continue position updating is as follows: m=2×sin(r+π / 2) In the formula, represents the optimal solution of the j-th dimension individual in the t-th iteration; F i represents the current position of any individual in the jth dimension at the tth iteration; F ri represents the fitness value of any individual at the random position of the jth dimension at the tth iteration; C(0,1) represents the Cauchy mutation; Step (4.3): According to the obtained best individual fitness value, BKA selects the best individual and saves it, that is, updates the hyperparameters of the GLIP prediction model; Where, X best represents the individual optimal solution; F best Represents the optimal fitness value of an individual; represents the optimal solution of the j-th dimension individual in the t-th iteration; and They represent the position solutions of the ith individual in the jth dimension at the tth iteration and the t+1th iteration respectively; Step (4.4): Determine whether the iterative optimization meets the maximum iteration condition. If it does, the iteration stops. Otherwise, the cycle continues until the optimal solution is obtained, i.e., the three hyperparameters α, β, and γ of the GLIP prediction model. Step (4.5): Substitute the three hyperparameters α, β, and γ of the GLIP prediction model obtained by the BKA optimization algorithm into the GLIP prediction model to predict the high-frequency component F1(t); input the high-frequency component F1(t), set its data sequence to X, the input length to L, the output length to O, and the hyperparameters α, β, and γ obtained after the above optimization, and split the high-frequency component data set into the training set X in a ratio of 7:1:

2. train , validation set X val and the test set X test ; Step (4.6): GLIP global prediction: for the entire training set X train Perform a discrete Fourier transform to obtain the amplitude A and frequency f; then, select the high amplitude frequency f from the amplitude * , then its corresponding period T * =1 / f * , let T * =[T1, T2, ...]; then, using the period T * To construct the global basis function Θg, the formula is as follows: E g =[sin(C1t),cos(C1t),1] Where, C1 = 2π / T * , t represents the length of the training set, T * Indicates high frequency period; The sparse identification process is obtained according to the global basis function, and the global parameter Ξ is solved through sparse identification. g , to achieve global prediction, the optimization formula of sparse recognition is as follows: Where X represents the training set data sequence; ||X-Θ g Ξ g || 2 Represents the global basis function for the training set X train The fitting error of g || indicates the global parameter Ξ g l1 norm regularization; By using the coordinate descent method similar to LASSO to solve the above optimization problem, we get Θ g ,make To predict the validation set X val and the test set X test The data in the prediction result is recorded as X gp ; Using the global prediction result X gp , construct a new basis function The formula is as follows: Where, X g =[X train ,X gp ], indicating the single variable prediction result X gp With training set X train splicing; The new basis function Alternative Θ g , re-solve the optimization problem of sparse identification and obtain the global prediction result considering the network coupling relationship between variables Step (4.7): Verify the recognition and determine whether to perform local rolling prediction: The GLIP prediction model is used for verification, to evaluate the effectiveness of global recognition and prediction, to observe whether the use of local features will perform better than global prediction, and to prepare for subsequent local rolling prediction; Assume that the global prediction value for the validation set data is Calculate the forecast error for each variable i and compare it with the local rolling forecast The error is compared; the implementation formula is as follows: Where N b Indicates the number of validation sets; Represents the sequence data set used for verification, that is, the true value; Representation sequence The global prediction value of Representation sequence The local rolling prediction value; k1 represents a hyperparameter; i represents a variable; L input Indicates the length of the input sequence; If all variables satisfy the above formula, it means that the global forecast is effective enough and there is no need for local rolling forecast. If some variables do not satisfy the above formula, it means that the global forecast performs poorly on these variables and local rolling forecast is needed. Step (4.8): Local prediction: Determine whether to adopt local rolling prediction based on the verification recognition result of global prediction; before local prediction, perform a simple preprocessing process on the input data, perform weighted processing on each window of the test set, and construct a local recognition curve x * , which is used to capture the global trend and local fluctuations of the data; its formula is as follows: In the formula, x represents the original input data; y represents the corresponding result obtained by global prediction; * Indicates the construction of a new input curve; ω represents the weight vector, which is used to control the proportion of the original input data x and the global prediction result y in the weighted average. Represents element-wise multiplication; ω=(linspace(α 1 / γ ,b 1 / γ ,L batch )) γ In the formula, the weight ω uses the linear interpolation function linspace to generate a linearly changing weight vector whose length is equal to L batch ; α and β represent the global and local weight coefficients respectively; γ represents the acceleration of weight change; Weighted input curve x * Perform Fourier transform to extract high-frequency components and convert them into potential cycles Next, use the cycle To construct the local basis function Θ l , the formula is as follows: I l =[sin(C3t),cos(C3t)] In the formula, Indicates potential cycle The coefficient of The global basis function Θ g , store basis function Θ s and the local basis function Θ l Fusion is performed to obtain the final basis function library Θ pred , the formula is as follows: I pred =[Θ g ,I s ,I l ]; Solving the local prediction parameter Ξ through sparse identification pred Implement local prediction and its sparse recognition optimization problem; set the prediction time to t = L input +1,L input +2,…,using Ξ pred and θ pred Perform local rolling prediction to obtain the final local prediction result Step (4.9): Output the final prediction result; when all variables meet the verification and identification conditions of the global prediction, the final prediction result is the global prediction result When some variables do not meet the verification and identification conditions of the global prediction, the final prediction result is the local prediction result. Finally, the prediction result of the high-frequency component F1(t) is obtained as 8. The method for predicting the remaining life of a lithium battery according to claim 1, characterized in that: In step (5), the implementation process of the BKA optimization algorithm to optimize the ESN prediction model for medium-frequency and low-frequency component prediction is as follows: Step (5.1): Use the BKA optimization algorithm to optimize the parameters of the ESN prediction model and find the optimal value; Step (5.2): Assume that the ESN prediction model has K input nodes, N reserve pool nodes, and L output nodes; split the time series data set of the intermediate frequency component F2(t) or the low frequency component F3(t) into a training set and a test set in a ratio of 7:3, and set them as the input data u(t). Then, the update formulas of the reserve pool state x(t) and the output result y(t) are as follows: x(t+1)=(1-α)x(t)+ατ[W in u(t+1)+Wx(t)+W back y(t)] y(t+1)=σ[W out (u(t+1),x(t+1),y(t))] Where x(t) and x(t+1) represent the state of the reserve pool at time t and time t+1; y(t) and y(t+1) represent the output state of the reserve pool at time t and time t+1; u(t+1) represents the input data accepted by the ESN model at time t+1; α represents the leakage rate, that is, the speed of the reserve pool update; W in , W, W back represents the input weight matrix, the reserve pool weight matrix, and the output feedback weight matrix; W out represents the output weight matrix; τ and σ represent the activation function; The output weight W is solved by regularized ridge regression method. out , the formula is as follows: Y target =W out X W out =Y target W T (XX T +βI) -1 Where Y target represents the target value; X represents the internal state matrix of the reserve pool, X T represents the transpose of the internal state matrix of the reservoir; β represents the regularization coefficient; I represents the identity matrix; W T represents the transpose of the weight matrix; Step (5.3): Input the test set data into the trained ESN prediction model, through the reserve pool state x(t) and output weight W out The prediction result is calculated as follows: y'(t)=tanh(W out x(t)) Where y'(t) represents the predicted output; Finally, the prediction results of the mid-frequency component F2(t) and the low-frequency component F3(t) are obtained. and 9. The method for predicting the remaining life of a lithium battery according to claim 8, characterized in that: In step (5.1), the parameters include the reserve pool spectrum radius SR, the reserve pool sparsity SD and the input unit scale IS, which are substituted into the ESN prediction model to predict the mid-frequency component F2(t) and the low-frequency component F3(t) respectively.

10. The method for predicting the remaining life of a lithium battery according to claim 1, characterized in that: In step (6), the implementation process of superimposing and reconstructing the prediction results is as follows: Prediction results for the obtained high-frequency component F1(t) Prediction results of mid-frequency component F2(t) and low-frequency component F3(t) and The weighted sum composite prediction model is used to integrate the prediction results of lithium batteries, and the formula is: Where Y represents the predicted value of the composite prediction model for the RUL of lithium batteries; α1, α2, and α3 represent random fusion coefficients in the range [0, 1], and α1∈[0.9, 1], α2∈[0.6, 0.9], and α3∈[0, 0.5].

Images