A method for predicting the remaining life of lithium batteries

Abstract

This invention discloses a method for predicting the remaining lifespan of lithium batteries, belonging to the field of lithium battery health management. The operational steps are as follows: Collecting operational data of the lithium battery; performing a primary decomposition of the lithium battery capacity sequence using the TVFEMD model; performing a secondary decomposition of the high-frequency and residual components in the IMF signal using the EWT model; calculating the sample entropy of the mid- and low-frequency components after the primary decomposition, and the high-, mid-, and low-frequency components after the secondary decomposition; performing cluster analysis on the classified component sets using the K-means clustering method; predicting the high-frequency components after K-means clustering using the BKA-GLIP model, and predicting the mid- and low-frequency components using the BKA-ESN model; reconstructing the prediction results of the two models to obtain the predicted remaining lifespan of the lithium battery. This invention helps improve the accuracy of predicting the remaining lifespan of lithium batteries, providing strong support for the use and management of lithium batteries.

Description

Technical Field

[0001] This invention belongs to the field of lithium battery remaining life prediction technology, and relates to a method for predicting the remaining life of lithium batteries. Background Technology

[0002] With increasing attention being paid to the development and utilization of new energy sources, lithium-ion batteries are widely used in various electronic devices and electric vehicles due to their advantages such as high energy density, long lifespan, and environmental friendliness. However, the performance of lithium batteries gradually declines during use until they can no longer meet the usage requirements. Therefore, predicting the remaining lifespan of lithium batteries is of great significance for ensuring the safe operation of equipment and improving battery utilization.

[0003] Currently, commonly used methods for predicting the remaining life of lithium batteries mainly include model-based prediction methods and data-based prediction methods. Model-based prediction methods predict the remaining life of the battery by establishing a mathematical model of the battery and using the battery's operating data and model parameters. However, due to the complexity and uncertainty of the battery's internal processes, establishing an accurate battery model is quite difficult. Data-based prediction methods collect battery operating data, such as voltage, current, charge / discharge time, temperature, and capacity, and then use a data-driven approach to predict the remaining life of the battery. This method does not require the establishment of a complex battery model, thus having better applicability. However, since battery operating data usually has nonlinear, non-stationary, and uncertain characteristics, how to effectively extract and process this data to improve the accuracy of prediction remains a challenging problem. Currently, traditional data processing methods are difficult to effectively extract the feature information of the data; and most commonly used prediction methods are based on a single prediction model, which cannot process signals of different frequencies simultaneously, thus affecting the accuracy of the prediction results. Summary of the Invention

[0004] To address the aforementioned problems, the present invention aims to provide a method for predicting the remaining lifespan of lithium batteries, which can improve the accuracy of predicting the remaining lifespan of lithium batteries and provide strong support for the use and management of lithium batteries.

[0005] The technical solution of the present invention is as follows: A method for predicting the remaining life of lithium batteries, the operation steps of which are as follows:

[0006] Step (1): Collect the operating data of the lithium battery; use the TVFEMD model to decompose the lithium battery capacity sequence into high frequency, medium frequency, low frequency and residual components according to the frequency from high to low;

[0007] The operational 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, mid-frequency and low-frequency components according to their frequency from high to low.

[0009] Step (3): Calculate the sample entropy for the mid-frequency and low-frequency components after the first decomposition and the high-frequency, mid-frequency and low-frequency components after the second decomposition. Based on the sample entropy evaluation results, classify the components of the two decompositions into high-frequency, mid-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: high-frequency, mid-frequency and low-frequency.

[0010] Step (4): Use the GLIP prediction model to predict the high-frequency components after K-means clustering, and use the ESN network to predict the mid-frequency and low-frequency components;

[0011] Step (5): To improve the prediction performance of the two prediction models, the BKA optimization algorithm is used to optimize the hyperparameters α, β, γ in the GLIP prediction model and the hyperparameters SR, SD, and IS in the ESN network.

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

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

[0014] Step (1.1): Let the original battery capacity signal be x(t), and use Hilbert transform to obtain the transformed signal x'(t), which is as follows:

[0015]

[0016] In the formula, 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); and x'(t) represents the Hilbert transform signal of the battery capacity signal x(t).

[0017] The formula for the analytic 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] In the formula, α1 and α2 represent the instantaneous amplitudes of the component signals. The phase difference represents the instantaneous amplitude;

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

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

[0022] In the formula, φ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 time t.

[0023] Step (1.2): Determine the local maxima X(t) of the instantaneous amplitude X(t). max Sequence and Local Minimum Sequence X(t) min The local maxima curve β1(t) and local minima curve β2(t) are obtained by interpolation. Then, the amplitudes α1(t) and α2(t), and the instantaneous frequency components φ1'(t) and φ2'(t) are calculated, and the formulas are as follows:

[0024]

[0025] In the formula, β1(t) and β2(t) represent local maxima and local minima, respectively; α1(t) and α2(t) represent the new amplitudes.

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

[0027]

[0028] For X 2 (t max )φ'(t max ) and X 2 (t min )φ'(t min After interpolation, η1(t) and η2(t) are obtained. 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 values ​​of X. 2 (t max )φ'(t max ) and X 2 (t min )φ'(t minThe interpolation is obtained by )

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

[0032]

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

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

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

[0036] The signal x(t) is continuously updated, and a stopping criterion θ(t) is defined. When the signal satisfies θ(t) ≤ ε, where ε is a given bandwidth threshold, the signal component at this point is denoted 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] The formula for the stopping criterion θ(t) is as follows:

[0038]

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

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

[0041] Step (1.4): After all IMF signals have been 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 performing a secondary decomposition of the high-frequency and residual components in the IMF signal obtained after the first decomposition using the EWT model is as follows:

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

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

[0046] Among them, the empirical scaling function With empirical wavelet function ψ n The formula for calculating (ω) 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 scaling coefficients are obtained by multiplying the empirical wavelet function with the original signal, and the empirical wavelet coefficients are obtained by multiplying the empirical wavelet function with the original signal. The formulas are as follows:

[0051]

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

[0053] Step (2.4): After the signal IMF1 and the residual signal R(t) are decomposed and reconstructed by the EWT algorithm, the following results are obtained:

[0054]

[0055] In the formula, k represents the number of components to be decomposed, and the 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 an IMF. n(t) For n∈1,2,…8, each IMF component is... n(t) Divide into m-dimensional continuous subsequences, as shown in the following equation:

[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] In the formula, N represents IMF i(t) Length;

[0061] For each subsequence x i Calculate other subsequences x j The Euclidean distance d(x) between them i ,x j Define sequence x i With 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, calculate d[x] i ,x j The number of vectors less than r is compared with the total number of vectors Nm, denoted as [[r]]. Its formula is as follows:

[0064]

[0065] In the formula, num{d[x i ,x j ]<r} represents 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, B represents the similarity of the i-th subsequence; m (r) represents the average similarity of m-dimensional vectors; SE m (i,r) represents the sample entropy of the i-th subsequence in m dimensions. Increasing the dimension m by 1 gives the similarity B for m+1 vectors. m+1 (r) and sample entropy SE m+1 (i,r);

[0070] Step (3.3): Combine the calculated sample entropies SE1, SE2, ..., SE8 of the eight IMF components into a sample entropy vector S. i S i ={SE1,SE2,…,SE8}; Select 3 initial cluster centers C1, C2, C3, representing high-frequency, mid-frequency, and low-frequency signals respectively; For each sample entropy S i Calculate its relationship with each center C j The distance is calculated and assigned to the nearest center; the distance is usually calculated using Euclidean distance, with the following formula:

[0071]

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

[0073] For each sample entropy value S i Find the nearest center C. j And it is allocated to this center, according to the following formula:

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

[0075] In the formula, C j* Distance S i The nearest center;

[0076] Step (3.4): Recalculate the center of each cluster, taking the average entropy of all samples assigned to that cluster as the new cluster center, using the following formula:

[0077]

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

[0079] Repeat steps (3.3) to (3.4) until the new cluster centers remain essentially unchanged. Finally, each IMF component will be assigned to the corresponding cluster, and the clustering results will be output: high-frequency component F1(t), mid-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 high-frequency components is as follows:

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

[0082] Where α and β represent the global and local weighting coefficients, respectively; γ represents the acceleration of weight change; let the population size of the BKA optimization algorithm be N, the maximum number of iterations be T, and the variable dimension be j; simultaneously, the f(·) function is introduced to represent the fitness value of an 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): BKA position update, 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 Let represent the positional solutions of the i-th individual in the j-th dimension at the t-th and t+1-th iterations, respectively; r is a random number between [0,1]; p represents a constant; t represents the number of iterations completed so far; T represents the total number of iterations; and n represents the adhesion coefficient.

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

[0088]

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

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

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

[0092]

[0093] In the formula, X best F represents the optimal solution for an individual. best This represents the optimal fitness value of an individual; Let represent the optimal solution for the j-th dimension individual in the t-th iteration; and Let these represent the positional solutions of the i-th individual in the j-th dimension at the t-th and t+1-th iterations, respectively.

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

[0095] Step (4.5): Based on the three hyperparameters α, β, and γ of the GLIP prediction model obtained by the BKA optimization algorithm, substitute them into the GLIP prediction model to predict the high-frequency component F1(t); input the high-frequency component F1(t), set its data sequence as X, input length as L, output length as O, and the hyperparameters α, β, and γ obtained after the above optimization, and divide the high-frequency component dataset into training set X in a ratio of 7:1:2. train Validation set X val and 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 higher amplitude frequency f from the amplitude. * Then its corresponding period T * =1 / f * 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] In the formula, C1 = 2π / T * t represents the duration of the training set, T * Indicates a high-frequency period;

[0099] The sparse identification process is derived from the global basis function, and the global parameter Ξ is solved through sparse identification. g To achieve global prediction, the optimization problem for sparse identification is formulated as follows:

[0100]

[0101] In the formula, X represents the training set data sequence; ||X-Θ g Ξ g || 2 This represents the global basis function on the training set X. train The fitting error; λ||Ξ g || indicates the global parameter Ξ g l1 norm regularization;

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

[0103] Using the global prediction result X gp Construct new basis functions The formula is as follows:

[0104]

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

[0106] New basis functions Replace Θ g The optimization problem of sparse identification is solved again to obtain a global prediction result that considers the network coupling relationship between variables.

[0107] Step (4.7): Verify and identify whether to perform local rolling prediction: The GLIP prediction model is verified and identified to evaluate the effectiveness of global identification and prediction, observe whether using local features will perform better than global prediction, and prepare for subsequent local rolling prediction.

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

[0109]

[0110] In the formula, N b Indicates the number of validation sets; This represents the sequence dataset used for verification, i.e., the true values; Represents a sequence The global predicted value; Represents a sequence The local rolling prediction value; k1 represents the hyperparameter; i represents the variable; L input Indicates the length of the input sequence;

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

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

[0113]

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

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

[0116] In the formula, the weight ω is generated using the linear interpolation function linspace, which produces a linearly varying weight vector with a length equal to L. batch α and β represent the global and local weighting coefficients, respectively; γ represents the acceleration of the weight change.

[0117] For the weighted input curve x * Perform a Fourier transform to extract the high-frequency components and convert them into latent periods. Next, utilize the cycle To construct local basis functions Θ l The formula is as follows:

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

[0119] In the formula, Indicates potential period The coefficient;

[0120] global basis function Θ g Storage basis function Θ s and local basis functions Θ l The base function library Θ is obtained by merging the base functions. pred The formula is as follows:

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

[0122] Solving the local prediction parameter Ξ using sparse identification pred The optimization problem of sparse identification in achieving local prediction; the prediction time is set to t = L. input +1,L input +2,…, using Ξ pred and Θ pred The final local prediction result is obtained by performing local rolling prediction.

[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 for global prediction, the final prediction result is a local prediction result. Finally, the prediction result of the high-frequency component F1(t) is obtained as follows:

[0124] Furthermore, in step (5), the implementation process of the BKA optimization algorithm to optimize the prediction of mid-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 solve for the optimal value;

[0126] The parameters include the reservoir spectral radius SR, the reservoir sparsity SD, and the input cell 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.

[0127] Step (5.2): Assume the ESN prediction model has K input nodes, N reservoir nodes, and L output nodes; divide the time series dataset of the mid-frequency component F2(t) or the low-frequency component F3(t) into a training set and a test set in a 7:3 ratio, and set them as the input data u(t). Then the update formulas for the reservoir 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] In the formula, x(t) and x(t+1) represent the states of the reservoir at time t and time t+1, respectively; y(t) and y(t+1) represent the output states of the reservoir at time t and time t+1, respectively; u(t+1) represents the input data received by the ESN model at time t+1; α represents the leakage rate, i.e., the rate at which the reservoir is updated; W in W, W back W represents the input weight matrix, the reservoir weight matrix, and the output feedback weight matrix. out The output weight matrix is ​​represented by τ and σ; τ and σ represent the activation functions.

[0131] The output weights W are solved using a 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] In the formula, Y target X represents the target value; X represents the internal state matrix of the reservoir. TThe transpose of the internal state matrix of the reservoir is represented by β; β represents the regularization coefficient; I represents the identity matrix; W T This represents the transpose of the weight matrix;

[0135] Step (5.3): Input the test set data into the trained ESN prediction model, using the reservoir state x(t) and output weights W. out The prediction result is calculated using the following formula:

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

[0137] In the formula, y'(t) represents the predicted output;

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

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

[0140] Based on the obtained high-frequency component F1(t) prediction results Prediction results of mid-frequency component F2(t) and low-frequency component F3(t) and The prediction results for lithium batteries obtained by using a weighted sum composite prediction model are fused, and the formula is as follows:

[0141]

[0142] In the formula, Y represents the predicted value of lithium battery RUL by the composite prediction model; α1, α2, and α3 represent random fusion coefficients in the range [0,1], generally set α1∈[0.9,1], α2∈[0.6,0.9], and α3∈[0,0.5].

[0143] The beneficial effects of this invention are as follows: 1. By employing the TVFEMD and EWT models to perform multi-level decomposition of battery capacity data, this invention can more accurately extract key signal components reflecting battery performance degradation. This refined signal processing method helps to capture subtle changes in the battery at different stages of use, thereby improving the accuracy of remaining lifespan prediction. 2. By combining K-means clustering and sample entropy calculation, the adaptability to signals of different frequencies is enhanced. This invention can effectively identify and separate high-frequency, mid-frequency, and low-frequency signal components. This frequency division processing method enables the prediction model to process signals with different characteristics more specifically, thereby improving the adaptability and overall prediction performance of the prediction model. 3. Optimizing the selection and application of prediction models, this invention adopts two prediction models, BKA-GLIP and BKA-ESN, which are suitable for processing signal components of different frequencies. This model-based prediction method can fully utilize the prediction advantages of different models, improving prediction efficiency and accuracy. Attached Figure Description

[0144] Figure 1 This is the overall system flowchart of the present invention;

[0145] Figure 2 This is a flowchart of the BKA-GLIP prediction model in this invention;

[0146] Figure 3 This is a flowchart of the BKA-ESN prediction model of the present invention. Detailed Implementation

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

[0148] As shown in the figure, the operation steps of the method for predicting the remaining life of lithium batteries according to the present invention are as follows:

[0149] Step 1: Collect the operating data of the lithium battery, including discharge voltage, discharge current, discharge time, and capacity sequence; use the TVFEMD model to decompose the lithium battery capacity sequence into high frequency, medium frequency, low frequency, and residual components according to frequency from high to low.

[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, mid-frequency, and low-frequency components according to their frequency from high to low.

[0151] Step 3: Calculate the sample entropy for the mid-frequency and low-frequency components after the first decomposition, and the high-frequency, mid-frequency and low-frequency components after the second decomposition. Based on the sample entropy evaluation results, classify the components of the two decompositions into high-frequency, mid-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: high-frequency, mid-frequency and low-frequency.

[0152] Step 4: Use the Global Local Information (GLIP) model to predict the high-frequency components after K-means clustering, and use the Echo State Network (ESN) model to predict the mid-frequency and low-frequency components. In order to improve the prediction performance of the two prediction models, use the Black-winged Kite Optimization Algorithm (BKA) to optimize the hyperparameters α, β, γ in the GLIP model and the hyperparameters SR, SD, IS in the ESN model.

[0153] Step 5: Reconstruct the prediction results of the two models to obtain the predicted remaining lifespan of the lithium battery. The specific flowchart is as follows: Figure 1 As shown.

[0154] The specific implementation process is as follows:

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

[0156] S1.1 In this invention, the preprocessed capacity data is used as the input signal. The original battery capacity signal is set as x(t), and the Hilbert transform is used to obtain the transformed signal x'(t). The specific formula is as follows:

[0157]

[0158]

[0159] In the formula, 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); and 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] In the formula, h(t) represents the analytic signal of x(t); α1 and α2 represent the instantaneous amplitudes of the component signals. 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] In the formula, φ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 time t.

[0166] S1.2 Determine the local maxima X(t) of the instantaneous amplitude X(t). max Sequence and Local Minimum Sequence X(t) min The local maxima curve β1(t) and local minima curve β2(t) are obtained by interpolation. Then, the amplitudes α1(t) and α2(t), and the instantaneous frequency components φ1'(t) and φ2'(t) are calculated. The specific formulas are as follows:

[0167]

[0168] In the formula, β1(t) and β2(t) represent local maxima and local minima, respectively; α1(t) and α2(t) represent the new amplitudes.

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

[0170]

[0171] For X 2 (t max )φ'(t max ) and X 2 (t min )φ'(t min After interpolation, η1(t) and η2(t) are obtained. Then, the instantaneous frequency components φ1'(t) and φ2'(t) are:

[0172]

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

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

[0175]

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

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

[0178] In the formula, φ b z(t) is the local cutoff frequency of the signal; when z(t) is at an extreme value, the extraction time point t is... max and t min As a 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 denoted as m(t), which represents the local mean.

[0179] The signal x(t) is continuously updated, and a stopping criterion θ(t) is defined. When the signal satisfies θ(t) ≤ ε, where ε is a given bandwidth threshold, the signal component at this point is denoted 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 for the stopping criterion θ(t) is as follows:

[0181]

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

[0183] From the above process, the TVFEMD model decomposes the lithium battery capacity sequence once. 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 After all IMF signals have been extracted, the remaining signal is the residual component R(t); that is, the following equation:

[0185]

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

[0187] S2.1. Perform a 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 for the Fourier Transform is as follows:

[0188]

[0189] S2.2 Construct the scaling function and wavelet function required for the EWT algorithm; using the empirical scaling function... With empirical wavelet function ψ n (ω) realizes the adaptive property of (H-1) bandpass filters, where the empirical scaling function is... With empirical wavelet function ψ n The formula for calculating (ω) 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 yields the scaling coefficient W(0,t), and the inner product of the empirical wavelet function and the original signal yields 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 a point in time; and Represents the empirical scaling function and the empirical wavelet function; and F represents the Fourier transform of the empirical scaling function and the empirical wavelet function; -1[·] indicates the inverse Fourier transform;

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

[0197] After decomposition and reconstruction of the signal IMF1 and the residual signal R(t) using the EWT algorithm, the following results are obtained:

[0198]

[0199] In the formula, k represents the number of components to be decomposed, 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 corresponding 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 obtained by the TVFEMD-EWT algorithm, the sample entropy of the IMF of each component 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 each IMF component is represented as an IMF. n(t) For n∈1,2,…8, each IMF component is... 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] In the formula, N represents IMF i(t) Length;

[0206] For each subsequence x i Calculate other subsequences x j The Euclidean distance d(x) between (where k≠j) i ,x j Define sequence x i With sequence x j The maximum distance between corresponding elements is d. max (x i,x j ); Calculate the Euclidean distance between each subsequence and all other subsequences, using the following formula:

[0207]

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

[0209]

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

[0211] Furthermore, the sample entropy of each subsequence is calculated for B. i m The average of the results (r) is then calculated, and a logarithmic transformation is performed 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, B represents the similarity of the i-th subsequence; m (r) represents the average similarity of m-dimensional vectors; SE m (i,r) represents the m-dimensional sample entropy of the i-th subsequence;

[0215] Furthermore, we increase the dimension by adding 1 to the dimension m, that is, we obtain 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 changes little, stop the iteration; finally, combine the calculated sample entropies SE1, SE2, ..., SE8 of the 8 IMF components into a sample entropy vector S. i S i ={SE1,SE2,…,SE8};

[0217] Based on the above calculation, the sample entropy vector S is obtained. i Three initial cluster centers, C1, C2, and C3, are selected to represent high-frequency, mid-frequency, and low-frequency signals, respectively; for each sample entropy S... i Calculate its relationship with each center C j The distance is calculated and assigned to the nearest center; the distance is usually calculated using Euclidean distance, as shown in the following formula:

[0218]

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

[0220] For each sample entropy value S i Find the nearest center C. j And it is allocated to this center, according to the following formula:

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

[0222] In the formula, C j* Distance S i The nearest center;

[0223] S3.4 Recalculate the center of each cluster, using the average entropy of all samples assigned to that 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] Repeat the iterative clustering process until the new cluster centers remain essentially unchanged. Finally, each IMF component will be assigned to its corresponding cluster, and the clustering results will be output as follows: high-frequency component F1(t), mid-frequency component F2(t), and low-frequency component F3(t).

[0227] In step (4): Based on the signal classification results of the K-means clustering above, 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 as X; the specific implementation process is as follows:

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

[0229] Assuming 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 Denotes an individual solution, i∈[1,N]; BK lb BK ub represents the upper and lower bounds of the solution range of the j-th individual, 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 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 X represents the optimal individual fitness value; L This 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, update the position of the candidate solution according to the fitness function to find better hyperparameter settings for GLIP;

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

[0238]

[0239] In the formula, and Let represent the positional solutions of the i-th individual in the j-th dimension at the t-th and t+1-th iterations, respectively; r is a random number between [0,1]; p represents a constant; t represents the number of iterations completed so far; T represents the total number of iterations; and n represents the adhesion coefficient.

[0240] BKA continues with position updates, using the following formula:

[0241]

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

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

[0244] S4.3 Based on the obtained optimal 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 F represents the optimal solution for an individual. best This represents the optimal fitness value of an individual; Let represent the optimal solution for the j-th dimension individual in the t-th iteration; and Let these represent the positional solutions of the i-th individual in the j-th dimension at the t-th and t+1-th iterations, respectively.

[0247] S4.4 Determine whether the iterative optimization meets the maximum iteration condition. If it does, the iteration stops; otherwise, continue the loop until the optimal solution is obtained, i.e., the three hyperparameters α, β, and γ of the GLIP prediction model.

[0248] S4.5. Based on the three hyperparameters α, β, and γ of the GLIP prediction model obtained by the BKA optimization algorithm, substitute them 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 as X, the input length is L, the output length is O, and the hyperparameters α, β, and γ obtained after the above optimization are defined; the high-frequency component dataset is divided into a training set X in a ratio of 7:1:2. train Validation set X val and 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 higher amplitude frequency f from the amplitude. * Then its corresponding period T * =1 / f * 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] In the formula, C1 = 2π / T * t represents the duration of the training set, T * Indicates a high-frequency period;

[0254] Furthermore, the sparse identification process can be derived from the global basis function, and the global parameter Ξ can be solved through sparse identification. g To achieve global prediction, the specific formula for optimizing sparse identification is as follows:

[0255]

[0256] In the formula, X represents the training set data sequence; ||X-Θ g Ξ g || 2 This represents the global basis function on the training set X. train The fitting error; λ||Ξ g || indicates the global parameter Ξ g l1 norm regularization is used to achieve sparsity;

[0257] Furthermore, by using a coordinate descent method similar to LASSO to solve the above optimization problem, Θ is obtained. g Let t = L train +1,L train +2,... to predict the validation set X val and test set X test The data in the dataset, and the prediction result is denoted as X. gp This means achieving global prediction;

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

[0259]

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

[0261] New basis functions Replace Θ g The optimization problem of sparse identification is solved again to obtain a global prediction result that takes into account the network coupling relationship between variables.

[0262] S4.7 Constructing Storage Basis Functions: Divide the test set into multiple batches, each batch containing an input sequence and an output sequence of a certain length; perform a 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; transform the high-frequency components extracted from all batches into their corresponding periods T. s * Using the extracted periodic set, a storage basis function Θ is constructed. s ;

[0263] The specific formula is as follows:

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

[0265] In the formula, C2 = 2π / T s * T s * This 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 undergoes validation to evaluate the effectiveness of global recognition and prediction, observe whether using local features performs better than global prediction, and prepare for subsequent local rolling predictions.

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

[0269]

[0270] In the formula, N bIndicates the number of validation sets; This represents the sequence dataset used for verification, i.e., the true values; Represents a sequence The global predicted value; Represents a sequence The local rolling prediction value; k1 represents the hyperparameter; i represents the variable; L input Indicates the length of the input sequence;

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

[0272] S4.9 Local Prediction: Based on the verification and recognition results of the global prediction, determine whether to adopt local rolling prediction; before local prediction, perform a simple preprocessing process on the input data, weight each window of the test set, and construct a local recognition curve x. * This is used to capture the global trend and local fluctuations of data; the specific formula is as follows:

[0273]

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

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

[0276] In the formula, the weight ω is generated using the linear interpolation function linspace, which produces a linearly varying weight vector with a length equal to L. batch Where α and β represent the global and local weighting coefficients, respectively; γ represents the acceleration of the weight change.

[0277] Furthermore, for the weighted input curve x * Perform a Fourier transform (DFT) to extract the high-frequency components and convert them into latent periods. Next, utilize the cycle To construct local basis functions Θ l The specific formula is as follows:

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

[0279] In the formula, Indicates potential period The coefficient;

[0280] Furthermore, the global basis function Θ g Storage basis function Θ s and local basis functions Θ l The base function library Θ is obtained by merging the base functions. pred The specific formula is as follows:

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

[0282] Furthermore, the local prediction parameter Ξ is solved through sparse identification. pred To achieve local prediction, the optimization problem of sparse identification is solved; the prediction time is set 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 validation and identification conditions for global prediction, the final prediction result is the global prediction result. When some variables do not meet the verification and identification conditions for global prediction, the final prediction result is a local prediction result. Finally, the prediction result of the high-frequency component F1(t) is obtained as follows:

[0285] In addition, the BKA-ESN prediction model is used to predict the mid-frequency component F2(t) and the low-frequency component F3(t). 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, and solve for the optimal values. These parameters include the reservoir spectral radius SR, the reservoir sparsity SD, and the input cell 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 up the ESN with K input nodes, N reservoir nodes, and L output nodes; divide the mid-frequency component F2(t) (or low-frequency component F3(t)) time series dataset into training and test sets in a 7:3 ratio, and set it as the input data u(t). The specific formulas for updating the reservoir 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] In the formula, x(t) and x(t+1) represent the states of the reservoir at time t and time t+1, respectively; y(t) and y(t+1) represent the output states of the reservoir at time t and time t+1, respectively; u(t+1) represents the input data received by the ESN model at time t+1; α represents the leakage rate, i.e., the rate at which the reservoir is updated; W in W, W back W represents the input weight matrix, the reservoir weight matrix, and the output feedback weight matrix. out The output weight matrix is ​​represented by τ and σ; τ and σ represent the activation functions.

[0291] The output weights W are solved using a 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] In the formula, Y target X represents the target value; X represents the internal state matrix of the reservoir. T The transpose of the internal state matrix of the reservoir is represented by β; β represents the regularization coefficient; I represents the identity matrix; W T This represents the transpose of the weight matrix;

[0295] Furthermore, the test set data is input into the trained ESN prediction model, using the reservoir state x(t) and output weights W. out The prediction result is calculated using the following formula:

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

[0297] In the formula, y'(t) represents the predicted output;

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

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

[0300]

[0301] In the formula, Y represents the predicted RUL of the lithium battery by the composite prediction model; α1, α2, and α3 represent random fusion coefficients in the range [0,1], generally set as α1∈[0.9,1], α2∈[0.6,0.9], and α3∈[0,0.6].

[0302] To evaluate the training performance of the prediction model, error analysis is required. The following error parameters are used as evaluation metrics: Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE) to validate the model's performance. The calculation formulas for these metrics 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, Indicates the predicted value; y i Indicates the actual value; n represents the quantity;

[0310] Furthermore, to ensure the accuracy of lithium battery RUL prediction, the collected remaining operating data of the lithium battery, including voltage, current, charge / discharge time, temperature, and other data, are input into the BKA-GLIP prediction model for further training and prediction, thereby improving the RUL prediction results of the lithium battery.

Claims

1. A method for predicting the remaining life of lithium batteries, characterized in that, The steps are as follows: S1: Collects operating data of the lithium battery; using The model decomposes the lithium battery capacity sequence into high-frequency, mid-frequency, low-frequency and residual components according to the frequency from high to low. S2: Use Model pair The high-frequency and residual components in the signal are decomposed into high-frequency, mid-frequency and low-frequency components according to their frequency from high to low. S3: Calculate the sample entropy for the mid-frequency and low-frequency components after the first decomposition, and the high-frequency, mid-frequency, and low-frequency components after the second decomposition. Based on the sample entropy evaluation results, classify the components from the two decompositions into sets of high-frequency, mid-frequency, and low-frequency components. Clustering methods perform cluster analysis on the classified component sets to obtain three signal components: high frequency, mid frequency, and low frequency. S4: Yes High-frequency components after clustering The prediction model makes predictions, using mid-frequency and low-frequency components. Networks make predictions; S5: To improve the prediction performance of the two prediction models, use Optimization algorithm Hyperparameters in prediction models , Hyperparameter reservoir spectral radius in the network 、 Sparseness of the reservoir 、 Input cell scale Optimize; The Optimization algorithm The prediction model performs the following process for mid-frequency and low-frequency component prediction: S51: Adopted Optimization algorithm The parameters of the prediction model are optimized, and the optimal value is obtained. S52: Let Predictive models have One input node, Each reserve pool node One output node; the intermediate frequency component or low-frequency components The time series dataset was split into training and test sets in a 7:3 ratio, and these sets were designated as the input data. Then the state of the reserve pool and output results The update formula is: In the formula, and Indicates that the reserve pool is Time and The state at any given moment; and Indicates that the reserve pool is Time and Output status at any given moment; Represents in time The input data accepted by the model; This indicates the leakage rate, which is the rate at which the reservoir is replenished. This represents the input weight matrix, the reserve pool weight matrix, and the output feedback weight matrix; This represents the output weight matrix; and Indicates the activation function; The output weights are solved using regularized ridge regression. Its formula is: In the formula, Indicates the target value; This represents the internal state matrix of the reservoir. This represents the transpose of the internal state matrix of the reservoir; Represents the regularization coefficient; Represents the identity matrix; This represents the transpose of the weight matrix; S53: Input the test set data into the trained... In the prediction model, the state of the reserve pool is used as a reference. and output weights The prediction result is calculated using the following formula: In the formula, Indicates the predicted output; Finally, the intermediate frequency component is obtained. and low frequency components Prediction results and ; S6: Finally, the prediction results of the two models are reconstructed to obtain the prediction results of the remaining service life of the lithium battery.

2. The method for predicting the remaining life of lithium batteries according to claim 1, characterized in that, The operating data described in S1 includes information on discharge voltage, discharge current, discharge time, and capacity sequence.

3. The method for predicting the remaining life of lithium batteries according to claim 1, characterized in that, S1 as described The model decomposes the lithium battery capacity sequence as follows: S11: Set the original battery capacity signal to The original battery capacity signal was obtained by using Hilbert transform. The transformed signal is Its formula is: In the formula, Indicates battery capacity signal The instantaneous amplitude; Indicates battery capacity signal The instantaneous phase; Indicates battery capacity signal Hilbert transform signal; Battery capacity signal Corresponding analytical signal The formula is: In the formula, Indicates signal The analytical signal, This represents the instantaneous amplitude of the component signal. The phase difference represents the instantaneous amplitude; Instantaneous amplitude The formula is: In the formula, This represents the instantaneous phase of the component signal. Indicates the first The instantaneous amplitude of the signal component at a given moment; S12: Determine the instantaneous amplitude Local maxima Sequence and Local Minimum Sequence And interpolate them separately to obtain local maximum curves. With local minimum curve Next, calculate the amplitude. and Instantaneous frequency components and Its formula is: In the formula, Represents local maxima and local minima; Indicates the new amplitude; Battery capacity signal instantaneous frequency components The formula is: ; right and Perform interpolation to obtain and Then the instantaneous frequency component and The formula is: In the formula, Indicates passage and The interpolation is used to obtain the result; Calculate the local cutoff frequency of the signal The formula is as follows: ; S13: Filter the uncalculated signal using a time-varying filter to obtain a local mean, and then calculate the signal. ; as shown in the following formula: In the formula, The signal's local cutoff frequency; when When at an extreme value, extract the time point. and As a node That is, to construct a B-spline approximation time-varying filter for the signal. Perform time-varying filtering using B-spline approximation, and denot the approximation result as follows. , representing the local mean; S14: When all After all the signals have been extracted, the remaining signals are the residual components. ; as shown in the following formula: .

4. The method for predicting the remaining life of lithium batteries according to claim 3, characterized in that, In S13, the signal is continuously updated. and define the stopping criteria therein. When the signal satisfies hour, Given a bandwidth threshold, the signal component at this point is denoted as... Component signal; otherwise, let Repeat steps S11-S13 to continue decomposing the signal until the stopping criterion is met, thus completing the decomposition. The stopping criteria The formula is: In the formula, Indicates the instantaneous bandwidth of the signal; This represents the instantaneous frequency of the weighted mean; According to S11-S13, The model performs a decomposition of the lithium battery capacity sequence, extracting a series of... The signal is divided into three components according to frequency, from high to low: high frequency, medium frequency, and low frequency. .

5. The method for predicting the remaining life of lithium batteries according to claim 1, characterized in that, The implementation process for the secondary decomposition described in S2 is as follows: S21: Transfer the signal held in S1 With residual signal Perform a Fourier transform; S22: Through empirical scaling function With empirical wavelet function Implement the adaptive properties of H-1 bandpass filters; The empirical scaling function With empirical wavelet function The calculation formula is: In the formula, Represents the empirical scaling function; Represents the empirical wavelet function; Indicates frequency, Indicates the first The frequency of each boundary; Represents the transformation function; Represents a coefficient function. ; S23: Based on the empirical scaling function The scaling coefficients are obtained by multiplying the empirical wavelet function with the original signal, and the empirical wavelet coefficients are obtained by multiplying the empirical wavelet function with the original signal. The formulas are as follows: In the formula, Indicates a point in time; 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; Indicates the inverse Fourier transform; S24: Signal With residual signal go through After decomposition and reconstruction by the algorithm, the result is: In the formula, This represents the number of components in the decomposition, and its value is... .

6. The method for predicting the remaining life of a lithium battery according to claim 1, characterized in that, The use described in S3 The implementation process of cluster analysis using clustering methods is as follows: S31: Let each Components , , each Quantity Divided into A continuous subsequence of dimension is given by the following formula: In the formula, express Length; For each subsequence Calculate other subsequences Euclidean distance between Define sequence with sequence The maximum distance between corresponding elements is ; Its formula is as follows: ; S32: Given a threshold ,statistics The number of vectors and the total number of vectors. The ratio is denoted as . Its formula is: In the formula, express The number; right The average of the results is calculated using the following formula: In the formula, Indicates the first The similarity of each subsequence; express The average similarity of dimensional vectors; Indicates the first Subsequences The sample entropy of dimension, and then the dimension Add 1, that is, for Similarity is obtained from the vectors. and sample entropy ; S33: Calculate the 8... Sample entropy of components Form a sample entropy vector , Choose 3 initial cluster centers. , representing high-frequency, mid-frequency, and low-frequency signals respectively; for each sample entropy Calculate its relationship with each center The distance is calculated and assigned to the nearest center; the distance is usually calculated using Euclidean distance, with the formula: In the formula, Indicates the first A sample entropy vector; Indicates the first Number of clusters, ; For each sample entropy value Find the center that is closest to it. And it is allocated to this center, using the following formula: In the formula, Indicates distance The nearest center; S34: Recalculate the center of each cluster, using the average entropy of all samples assigned to that cluster as the new cluster center. The formula is as follows: In the formula, Indicates assignment to cluster The number of samples; Repeat steps S33-S34 until the new cluster centers remain essentially unchanged. Finally, each... The components will be assigned to the corresponding clusters, and the clustering results will be output: high-frequency components. Intermediate frequency components with low-frequency components .

7. The method for predicting the remaining life of lithium batteries according to claim 1, characterized in that, S4 adopts The process by which the prediction model predicts high-frequency components is as follows: S41: Optimization algorithm initialization: Randomly generate a set of candidate solutions, each individual solution representing... One of the parameters of the model is the hyperparameter. ; in, These represent the global and local weighting coefficients, respectively. Let the acceleration be a change in weight; The population size of the optimization algorithm is The maximum number of iterations is Variable dimension is At the same time, introduce The function represents the fitness value of an individual. During initialization, The individual with the best fitness value is selected as the optimal solution in the initial population. The optimal solution is set as follows: Hyperparameters of the prediction model; S42: Location update, iterative optimization The hyperparameters of the prediction model are determined; in each iteration, the positions of candidate solutions are updated according to the fitness function to find... Hyperparameter settings for the prediction model; The The formula for position update is: In the formula, and They represent the first The individual in the first Victor The second iteration and the first The positional solution in the next iteration; for Random numbers between; Represents a constant; This indicates the number of iterations that have been completed so far; Indicates the total number of iterations; Indicates the adhesion coefficient; The formula for continuing position updates is: In the formula, Indicates the first The iteration of the ... The optimal solution for a given individual; It means that any individual in the first... The iteration of the ... The current position of the dimension; It means that any individual in the first... The iteration of the ... Fitness values ​​at random locations in a dimension; Indicates Cauchy mutation; S43: Based on the obtained optimal individual fitness value Select the best individual and save it; that is, update. Hyperparameters of the prediction model; In the formula, Represents the optimal solution for an individual; This represents the optimal fitness value of an individual; Indicates the first The iteration of the ... The optimal solution for a given individual; and They represent the first The individual in the first Victor The second iteration and the first The positional solution in the next iteration; S44: Determine if the iterative optimization satisfies the maximum iteration condition. If it does, stop the iteration; otherwise, continue the loop until the optimal solution is obtained. Three hyperparameters of the prediction model ; S45: According to Optimization algorithm obtained Three hyperparameters of the prediction model Substitute it into High-frequency components in the prediction model Prediction; input high-frequency components Set its data sequence as The input length is The output length is and the hyperparameters obtained after the above optimization. The high-frequency component dataset was split into training sets in a 7:1:2 ratio. Validation set and test set ; S46: Global prediction: For the entire training set Perform a discrete Fourier transform to obtain the amplitude. and frequency Then, select high-amplitude frequencies from the amplitude range. Then its corresponding period ,make Next, using the cycle To construct global base functions Its formula is: In the formula, , Indicates the duration of the training set. Indicates a high-frequency period; The sparse identification process is derived from the global basis function, and the global parameters are solved through sparse identification. To achieve global prediction, the optimization problem for sparse identification is formulated as follows: In the formula, This represents a sequence of training set data. Represents the global basis function on the training set The fitting error; Indicates global parameters of Norm regularization; By using a coordinate descent method similar to LASSO to solve the above optimization problem, we obtain... ,make To predict the validation set and test set The data in the image, the prediction result is denoted as ; Using global prediction results Construct new basis functions Its formula is: In the formula, Indicates univariate prediction results With training set splicing; New basis functions Alternative The optimization problem of sparse identification is solved again to obtain a global prediction result that considers the network coupling relationship between variables. ; S47: Verification and identification to determine whether to perform local rolling prediction: The GLIP prediction model is verified to evaluate the effectiveness of global identification and prediction, observe whether using local features will outperform global prediction, and prepare for subsequent local rolling prediction. Suppose the global prediction for the validation set data is... Calculate the prediction error for each variable and compare it with the local rolling prediction. The error is compared; the formula is as follows: In the formula, Indicates the number of validation sets; This represents the sequence dataset used for verification, i.e., the true values; Represents a sequence The global predicted value; Represents a sequence Local rolling prediction value; Indicates hyperparameters; Represents variables; Indicates the length of the input sequence; If all variables satisfy the above formula, it means that the global prediction is effective enough and there is no need to perform local rolling prediction; if some variables do not satisfy the above formula, it means that the global prediction performs poorly on these variables and local rolling prediction is needed. S48: Local Prediction: Based on the verification and recognition results of the global prediction, determine whether to adopt local rolling prediction; before local prediction, perform a simple preprocessing process on the input data, weight each window of the test set, and construct a local recognition curve. It is used to capture the global trend and local fluctuations of data; its formula is: In the formula, This represents the original input data; This represents the corresponding result obtained from the global prediction; This indicates the construction of a new input curve; This represents a weight vector used to control the original input data. and global prediction results The proportion in the weighted average Represents element multiplication; In the formula, the weights Using linear interpolation functions Generate a linearly varying weight vector with a length equal to... ; These represent the global and local weighting coefficients, respectively. This indicates the acceleration due to the change in weight. For weighted input curves Perform a Fourier transform to extract the high-frequency components and convert them into latent periods. Next, using the cycle To construct local basis functions Its formula is: In the formula, Indicates potential period The coefficient; global basis functions Storage basis functions and local basis functions The base function library is then merged to obtain the final base function library. Its formula is: ; Local prediction parameters are solved by sparse identification. Achieving local prediction involves optimizing sparse identification; the prediction time is set to... ,use and The final local prediction result is obtained by performing local rolling prediction. ; S49: Output the final prediction result; when all variables meet the verification and identification conditions for global prediction, the final prediction result is the global prediction result. When some variables do not meet the verification and identification conditions for global prediction, the final prediction result is a local prediction result. Finally, the high-frequency components are obtained. The prediction result is .

8. The method for predicting the remaining life of a lithium battery according to claim 1, characterized in that, The parameters described in S51 include the reservoir spectral radius. Sparseness of the reserve pool and input unit scale Substitute Predictive model for mid-frequency components and low frequency components Make separate predictions.

9. The method for predicting the remaining life of a lithium battery according to claim 1, characterized in that, The process of superimposing and reconstructing the various prediction results described in S6 is as follows: For the obtained high frequency components Prediction results Mid-frequency components and low frequency components Prediction results and The prediction results for lithium batteries are obtained by fusing the results using a weighted sum composite prediction model, and the formula is as follows: In the formula, This indicates the composite prediction model for lithium batteries. The predicted value; Indicates within the range The random fusion coefficient between them is set. .

Images