CN114779088B - Battery remaining service life prediction method based on maximum expected-unscented particle filtering - Google Patents

Abstract

The invention relates to a battery remaining service life prediction method based on maximum expected-unscented particle filtering. The invention aims to solve the problem that the sudden increase of the capacity of the existing battery causes great error in the prediction of the residual service life of the battery. The process is as follows: 1. extracting battery capacity data in the kth working process; 2. constructing a dynamic battery degradation model based on unscented particle filtering; 3. adaptively estimating process noise and measurement noise in the battery degradation model by using a maximum expectation-based algorithm; the process noise and the measurement noise are used for constructing a dynamic battery degradation model based on unscented particle filtering in the k+1th working process; 4. judging whether capacity regeneration occurs in the kth working process; 5. solving a confidence interval of the residual service life and the battery capacity of the battery; the method is used for the field of prediction of the residual service life of the battery.

Description

Battery remaining service life prediction method based on maximum expected-unscented particle filtering

Technical Field

The invention relates to the discipline intersection field combining state space theory, statistical analysis and residual service life prediction of a lithium ion battery, in particular to a method for predicting the residual service life of the battery.

Background

As an important energy storage device, a lithium ion battery has the advantages of high energy density, wide working temperature range, long cycle life and the like, and becomes an important part of an electric automobile. As the number of charging and discharging times of the lithium ion battery increases, the capacity of the battery gradually decreases, which greatly reduces the reliability of the system or the device, and even in serious cases, catastrophic accidents are liable to occur. Therefore, accurately predicting the remaining life (Remaining useful Life, RUL) of a lithium ion battery, in other words, accurately estimating the time that the battery has elapsed from the current time to the time of complete failure has a critical role in the safety and reliability of the system. The method has the advantages of reducing the maintenance cost of the electric equipment, reducing the probability of fault occurrence and realizing the maintenance according to conditions.

The current residual service life prediction method of the lithium ion battery can be generally divided into a method based on a physical failure model, a data driving method and a method based on a mixed experience model. Because of the great difficulty in modeling the exact battery mechanism, there are certain limitations to methods based on physical failure models. On the other hand, the data driving method can predict the residual service life through the data generated in the working process of the battery without knowing the physical and chemical knowledge of the battery in advance. However, the data-driven learning-class approach requires that the online data also be required to follow a distribution similar to the historical data, otherwise it is difficult for the learning-class approach to guarantee the accuracy of the predictions. The method based on the mixed experience model has strong interpretation and reliability, an algorithm model between parameters such as battery health factors, current, voltage and cycle times is built by integrating priori knowledge into the model, and then algorithms such as extended Kalman filtering (Extended Kalman Filter, EKF), unscented Kalman filtering (Unscented Kalman Filter, UKF) and Particle Filter (PF) are adopted for parameter identification, so that the residual service life prediction task is completed.

While data-driven methods are widely used in the field of battery life prediction, these learning-based methods are black-box models in nature, require extensive historical data training algorithm models, and do not quantify the uncertainty of the battery in the degradation process. It is noted that, in general, the capacity of a lithium ion battery gradually decreases with an increase in the number of charge and discharge cycles. However, when the battery is in a stationary state, electrochemical reaction products in the battery are reduced, and electrochemical performance is briefly restored, which is a regeneration phenomenon of the battery capacity. A sudden increase in battery capacity will create a significant error in the prediction of the remaining useful life of the battery. Therefore, it is important to accurately judge the regeneration point of the battery capacity. Meanwhile, in some existing methods based on mixed experience models, both process noise and measurement noise of the models need to be set manually, and how to estimate the noise in an adaptive manner needs to be studied intensively.

Disclosure of Invention

The invention aims to solve the problem that the sudden increase of the existing battery capacity will cause great error to the residual service life prediction of the battery, and provides a battery residual service life prediction method based on maximum expected-unscented particle filtering.

The method for predicting the remaining service life of the battery through maximum expected-unscented particle filtering specifically comprises the following steps:

step 1, extracting battery capacity data in the kth working process;

step 2, constructing a dynamic battery degradation model based on unscented particle filtering according to the battery capacity data in the kth working process extracted in the step 1;

step 3, adaptively estimating process noise and measurement noise in the battery degradation model by adopting a maximum expectation-based algorithm;

the process noise and the measurement noise are used for constructing a dynamic battery degradation model based on unscented particle filtering in the k+1th working process;

step 4, judging whether capacity regeneration occurs in the kth working process;

step 5, solving the residual service life of the battery through trend extrapolation according to the battery degradation model constructed in the step 2 and the current battery capacity value;

and (3) calculating a confidence interval of the battery capacity according to the battery degradation model constructed in the step (2) and the current battery capacity value.

The beneficial effects of the invention are as follows:

the invention provides a novel method for predicting the residual service life of a battery based on unscented particle filtering-maximum expectation-Wilcoxon rank sum test aiming at a single battery. On the premise of not needing a large amount of historical data, a dynamic degradation model based on unscented particle filtering is constructed only according to battery capacity data to describe the degradation process of the battery. The invention adopts a method of adaptively estimating process noise and measurement noise in a dynamic degradation model based on a maximum expected algorithm, and accurately detecting the regeneration point of the battery capacity by a Wilcoxon rank sum test method so as to further predict the residual service life of the battery.

1. The traditional data driving method requires a large amount of historical data to train an algorithm model, the invention provides a dynamic degradation model based on unscented particle filtering for a single battery without the need of labeled historical data, and provides a maximum expected algorithm to adaptively update model process noise and measurement noise parameters by using capacity data in the battery working process.

2. For the battery capacity regeneration point, the invention provides a method based on Wilcoxon rank sum test, which measures the difference between the prior probability density distribution and the posterior probability density distribution of particles, thereby accurately detecting the battery capacity regeneration point and reducing the residual service life prediction error caused by the sudden increase of the battery capacity.

3. Considering that the uncertainty of the battery in the degradation process is difficult to quantify, the prediction method provided by the invention can determine the confidence interval of the battery capacity through the probability distribution of particles in the unscented particle filtering algorithm, thereby describing the uncertainty in the degradation process of the battery.

Drawings

FIG. 1 is a workflow diagram of the present invention;

FIG. 2 is a schematic diagram of a battery capacity fade curve for 4 batteries in a data set provided by the NASA Ames prediction center electricity;

FIG. 3a is a graph of the residual life prediction method of the present invention for a degradation model in a B0005 battery

Is a self-adaptive estimation result graph of (1);

FIG. 3B shows the residual life prediction method of the present invention for the degradation model in a B0005 battery

Adaptive estimation result map of (2)；

FIG. 3c shows the residual life prediction method of the present invention for the degradation model of the B0005 battery

Is a self-adaptive estimation result graph of (1); />

FIG. 3d is a graph of the residual life prediction method of the present invention for a degradation model in a B0005 battery

Is a self-adaptive estimation result graph of (1);

FIG. 3e shows the residual life prediction method of the present invention for the degradation model of the B0005 battery

Is a self-adaptive estimation result graph of (1);

FIG. 4a is a graph of confidence results for a capacity regeneration point test for battery number B0005 using the Wilcoxon rank sum test method;

fig. 4B is a display of a battery capacity regeneration point numbered B0005;

FIG. 5a is a graph showing the results of the residual life prediction method of the present invention for a residual life prediction value and a true value of a B0005 battery under various cycles, wherein EM-UPF represents a maximum expected-unscented particle filter algorithm, EM-UPF-W represents a maximum expected-unscented particle filter-Wilcoxon rank sum test algorithm, and RUL group Truth represents a true value of the residual life;

FIG. 5B is a graph showing the results of the residual life prediction method of the present invention for a residual life prediction value and a true value of a B0006 battery under different cycles, wherein EM-UPF represents a maximum expected-unscented particle filter algorithm, EM-UPF-W represents a maximum expected-unscented particle filter-Wilcoxon rank sum test algorithm, and RUL group Truth represents a true value of the residual life;

FIG. 5c is a graph showing the results of the residual life prediction method of the present invention for a residual life prediction value and a true value of a B0007 battery under different cycles, wherein EM-UPF represents a maximum expected-unscented particle filter algorithm, EM-UPF-W represents a maximum expected-unscented particle filter-Wilcoxon rank sum test algorithm, and RUL group Truth represents a true value of the residual life;

FIG. 5d is a graph showing the results of the residual life prediction method of the present invention for the residual life prediction value and the true value of the B0018 battery under different cycles, wherein EM-UPF represents the maximum expected-unscented particle filter algorithm, EM-UPF-W represents the maximum expected-unscented particle filter-Wilcoxon rank sum test algorithm, and RUL group Truth represents the true value of the residual life;

FIG. 6a is a graph of a 99% confidence interval for battery capacity at the 60 th duty cycle for B0005 for the remaining useful life prediction method of the present invention;

FIG. 6B is a graph of a 99% confidence interval for battery capacity of B0006 at the 60 th duty cycle for the remaining useful life prediction method of the present invention;

FIG. 6c is a graph of a 99% confidence interval for battery capacity of B0007 at the 60 th duty cycle for the remaining useful life prediction method of the present invention;

fig. 6d is a graph of the remaining useful life prediction method of the present invention with respect to the 99% confidence interval of battery capacity of B0018 at the 60 th duty cycle.

Detailed Description

The first embodiment is as follows: the battery remaining service life prediction method based on maximum expected-unscented particle filtering in the embodiment specifically comprises the following steps:

step 1, extracting battery capacity data in the kth working process;

and extracting capacity data of the battery in the kth practical working cycle, namely performing charge and discharge test after each working cycle of the battery to obtain battery capacity data under the corresponding cycle, and taking the battery capacity data as input of a battery residual service life prediction method.

Step 2, constructing a dynamic battery degradation model based on unscented particle filtering according to the battery capacity data in the kth working process extracted in the step 1, and describing a battery degradation process;

step 3, adaptively estimating process noise and measurement noise in the battery degradation model by adopting a maximum expectation-based algorithm;

The process noise and the measurement noise are used for constructing a dynamic battery degradation model based on unscented particle filtering in the k+1th working process;

specifically, for process noise and measurement noise in a dynamic degradation model based on unscented particle filtering, performing adaptive estimation by adopting a maximum expected algorithm, and recursively updating the process noise and the measurement noise by using all capacity data of a battery from an initial state to a current moment;

step 4, judging whether capacity regeneration occurs in the kth working process;

the invention provides a method based on Wilcoxon rank sum test, which is used for measuring the difference between the prior probability density distribution and the posterior probability density distribution of particles in unscented particle filter (Unscented Kalman Filter, UPF), so that the regeneration point of the battery capacity is accurately detected, and the residual service life prediction error caused by the sudden increase of the battery capacity is reduced;

step 5, solving the residual service life of the battery through trend extrapolation according to the battery degradation model constructed in the step 2 and the current battery capacity value;

calculating a confidence interval of the battery capacity according to the battery degradation model constructed in the step 2 and the current battery capacity value; evaluating the on-line prediction effect of the remaining service life: the effect of the battery remaining service life prediction method proposed by the invention is evaluated by adopting two indexes of root mean square error (Root Mean Square Error, RMSE) and absolute average error (Mean Absolute Error, MAE).

The second embodiment is as follows: the first difference between the present embodiment and the specific embodiment is that in the step 2, a dynamic battery degradation model based on unscented particle filtering is constructed according to the battery capacity data in the kth working process extracted in the step 1, so as to describe the degradation process of the battery; the specific process is as follows:

the particle filtering has the advantage of being capable of handling the parameter identification problems of nonlinear and non-gaussian systems relative to the kalman filtering algorithm and the extended kalman filtering algorithm. It should be noted that the conventional particle filtering algorithm directly uses the prior probability density function as the importance function. However, the particles generated from the a priori distribution are not necessarily capable of reflecting the actual particle distribution and the measurement information at the current moment of the system is ignored. For this reason, the present invention contemplates the use of unscented Kalman filtering algorithms to generate importance functions for particle filtering. Specifically, the posterior probability density distribution of the unscented kalman filter algorithm is regarded as the prior probability density distribution of the particle filter.

The battery degradation model based on unscented particle filtering can be expressed in the form of a state space as shown in formula (1):

wherein ,x_k Is the state variable at the kth working cycle, x _k ＝[a _k ,b _k ,c _k ,d _k ] ^T ，a _k A first component, b, being a state variable at the kth duty cycle _k A second component, c, being a state variable at the kth duty cycle _k D, being the third component of the state variable at the kth duty cycle _k A fourth component of the state variable at the kth duty cycle, T being the transpose, k representing the number of battery duty cycles, x _k-1 Is the state variable at the k-1 th working cycle, y _k For battery capacity, u _k U is a process noise term _k ＝[u _a ,u _b ,u _c ,u _d ] ^T ，u _a Noise term, u, being the first component of the state variable at the kth duty cycle _b Noise term, u, being the second component of the state variable at the kth duty cycle _c A noise term which is the third component of the state variable at the kth duty cycle, u _d A noise term, v, which is the fourth component of the state variable at the kth duty cycle _k To measure noise term, v _k ∈R ^1×1 1×1 is a matrix with a length and width of 1, f () represents a mapping relation of a state transition equation, and h () represents a mapping relation of a measurement equation;

state transition equation x in equation (1) _k ＝f(x _k-1 )+u _k Can be written in the form of equation (2):

wherein ,

for the process variable noise term u _a Variance of->

For the process variable noise term u _b Variance of->

For the process variable noise term u _c Variance of->

For the process variable noise term u _d Is a variance of (2); a, a _k-1 A first component, b, being a state variable at the kth-1 duty cycle _k-1 A second component, c, being a state variable at the kth-1 duty cycle _k-1 A third component, d, being a state variable at the (k-1) th duty cycle _k-1 A fourth component that is a state variable at the kth-1 duty cycle; />

Represents u _a Obeying the mean value to be 0, variance to be +.>

Normal distribution of->

Mean value 0, variance +.>

Normal division of (2)The cloth is made of a cloth material,

mean value 0, variance +.>

Normal distribution of->

Is 0 as mean and 0 as variance->

Is used for the normal distribution of the (c),

mean value 0, variance +.>

Is a normal distribution of (2);

accordingly, measurement equation y _k ＝h(x _k )+v _k Can be written in the form of equation (3):

wherein ,

for measuring the variance of the noise; />

Mean value 0, variance +.>

Is a normal distribution of (2);

the implementation steps of the whole unscented particle filter algorithm can be summarized as follows:

(1) Particle initialization, the process is:

according to the initial state variable x ₀ In the initial state variable x ₀ Normal distribution p (x) ₀ ) The particle set can be generated by the Monte Carlo method, namely a random sampling mode

N _s Is the total number of particles, and the initialization weight value of each particle is +.>

Initializing the mean of particles

Initializing the variance S of particles ₀ Can be written in the form shown in the formulas (4) - (5):

Wherein i is the number of the particle,

to initialize the mean value of the particles, S ₀ To initialize the variance of the particles, E () represents the desired operator, T represents the matrix transpose operator, r ₀ To initialize the particles;

(2) Calculating state variable r characterized by particles through UKF algorithm _k Posterior expectation of (2)

Sum of covariance->

The process is as follows:

after completing the k-1 iterative process, a Sigma sampling point set of each particle is constructed

As shown in formula (6):

/>

wherein ,

for the mean vector of the ith particle during k-1 work>

The covariance matrix of the ith particle in the k-1 working process is adopted, and the parameter alpha is a coefficient value in UKF;

the number of sample points of the Sigma point set of each particle is 2n+1, n is the number of state variables, and the weight expressions of the sample points of the Sigma point set of each particle are shown in formulas (7) - (9):

wherein ,

process variable weights represented by sample points of the ith particle, sigma Point set 1, for example>

Weight of the measurement variable represented by the sample point of the 1 st Sigma point set for the ith particle, for>

Process variable weights represented by sample points of the jth Sigma point set for the ith particle, +.>

The weight of a measured variable represented by a sample point of the jth Sigma point set of the ith particle, the parameter beta, lambda are coefficient values in UKF, and n is the number of state variables;

On this basis, the state estimation expression of the UKF can be written in the form shown in formulas (10) - (11):

wherein ,

for the prior state vector represented by the j Sigma point in the ith particle of the kth working process, f () is the mapping relation of the state transition equation, +.>

For the posterior state vector represented by the j-th Sigma point in the ith particle of the kth-1 th course of operation,/->

For the noise variable represented by the ith particle of the kth-1 working process +.>

A weighted prior state vector for the ith particle of the kth course of operation;

the expression for covariance estimation can then be written in the form shown in equation (12):

wherein ,

a priori covariance of the state variables in the ith particle for the kth course of operation, +.>

Covariance matrix of noise variable in ith particle for kth working process, +.>

diag represents a diagonal matrix, < >>

For the kth process variable noise term u _a Variance of->

For the kth process variable noise term u _b Variance of->

Process variable noise term u for the kth course of operation _c Variance of->

For the kth process variable noise term u _d Is a variance of (2);

accordingly, for forward one-step prediction, there are forms shown in formulas (13) - (14):

wherein ,

for the state vector represented by the jth Sigma point in the ith particle, +.>

For UKF forward one-step prediction represented by the jth Sigma point in the ith particle,/I>

Predicted value of UKF represented by the ith particle, < >>

For the measurement noise item of the ith particle in the kth working process, h () represents the mapping relation of the measurement equation;

based on this, an autocovariance matrix can be calculated

And cross covariance matrix->

As shown in the formulas (15) - (16):

wherein ,

noise variance of the ith particle for the kth course of operation;

thereby, a Kalman filtering gain vector K can be obtained _k Such as formula(17) The following is shown:

wherein ,

for the ith particle Kalman filter gain vector, -/->

Is an autocovariance matrix>

As for the cross covariance matrix, with the addition of new battery capacity data in the kth iteration, the state update and covariance update of the UKF can be expressed as the forms shown in formulas (18) - (19):

wherein ,

to calculate the state variable r of a particle by the UKF algorithm _k Posterior expectation (updated value of UKF algorithm state variable),>

to calculate the state variable r of a particle by the UKF algorithm _k Covariance matrix (updated values of the UKF algorithm covariance matrix);

according to the results in the formulas (18) - (19), for

Obtaining the mean value->

For->

Obtaining the mean value->

Based on

and />

Obtaining posterior probability Density distribution of particles under UKF Algorithm->

(3) The particle weight is updated by the following steps:

the importance density function of the particles obtained according to (2) can be expressed as:

wherein ,

importance density function representing the ith particle of the kth course of operation,/for>

For the state variables represented by the ith particle initial state through the kth-1 working process, y _1:k Battery capacity data indicating the 1 st to kth operation, +.>

State variables represented for the ith particle initial state through the kth-1 th operation and the 1 st to kth operationUnder the prior condition of pool capacity data +.>

Is a function of the probability density of (c) in the (c),

is the mean value of +.>

Variance is->

Normal distribution of (a) -is the meaning obeying a certain distribution;

on this basis, the weight update expression of the particle can be written as in formula (20):

/>

wherein ,

for the weight of the ith particle under k-1 duty cycles,/for the particle>

Weight of the ith particle under k duty cycles, +.>

Is y _k At->

Probability density function under conditions->

Is->

At->

Probability density function under conditions;

after the weight normalization of the particles, normalized particle weights can be obtained >

The expression is as follows:

(4) The particle resampling process comprises the following steps:

in order to improve the state estimation performance of the particle filtering algorithm, the invention adopts a resampling mode to reserve particles with larger weight in a particle set and reduce the particles with smaller weight. Let the effective particle number be N _ef Can be defined as in equation (22):

setting a threshold value N _th If the effective particle number N _ef Less than N _th Resampling of the particles is required;

the random resampling implementation process is simple and has higher sampling efficiency, so the invention adopts a random resampling mode to realize the resampling process of particle filtering;

if the effective particle number N _ef N is greater than or equal to _th Then no resampling of the particles is required;

after completion of (4), each particle has a weight of

(the weight of each particle requiring resampling and not requiring resampling is +.>

)；

(5) The state updating process comprises the following steps:

based on the weight of the particles obtained by resampling, the expected state variable under the posterior condition of the unscented particle filter algorithm (UPF) can be obtained through formulas (23) - (24)

Sum of covariance->

wherein ,

for state variables under the posterior condition of the UPF algorithm, it is desirable to have, < +.>

Is covariance under the posterior condition of the UPF algorithm.

Other steps and parameters are the same as in the first embodiment.

And a third specific embodiment: the difference between the present embodiment and the first or second embodiment is that in the step 3, a process noise and a measurement noise in the battery degradation model are adaptively estimated based on a maximum expectation algorithm;

specifically, for process noise and measurement noise in a dynamic degradation model based on unscented particle filtering, performing adaptive estimation by adopting a maximum expected algorithm, and recursively updating the process noise and the measurement noise by using all capacity data of a battery from an initial state to a current moment;

the specific process is as follows:

the unknown parameter vector composed of the process noise and the measurement noise can be expressed as:

/>

the essence of the EM algorithm is that by maximizing the joint likelihood function p (x _0:k ,y _1:k I) to achieve maximum likelihood estimation of the unknown parameter. Battery capacity data y from initial to kth cycle is known _1:k ＝[y ₁ ,y2,…,y _k ]On the premise of (1) can construct a joint log-likelihood function L _k (xi) is as shown in formula (25):

wherein ,

for the state variable of the ith particle of the t-th working process,/for the first time>

N, which is the value of the process variable in the t-th working process of the ith particle _s Is the total number of particles, and the initialization weight value of each particle is +.>

i is the number of the particle, ">

State variables characterized by the ith particle for the 1 st to kth course of operation, +.>

A joint probability density function composed of state variables represented by the ith particles from the 1 st working procedure to the kth working procedure under the unknown parameter condition and battery capacity data from the 1 st working procedure to the kth working procedure,>

a conditional probability density function for a state variable characterized by the ith particle of the 1 st to kth operation under unknown parameter conditions, +.>

Conditional probability density function of battery capacity data for 1 st to kth operation under conditions of unknown parameters and state variables characterized by the 1 st to kth operation's ith particles, +.>

A conditional probability density function of the ith particle state variable of the t-1 th working process under the unknown parameter condition and the ith particle state variable of the t-1 th working process,>

conditional probability density function of battery capacity data of the t-th working process under unknown parameter conditions and the state variable of the i-th particle of the t-th working process, < +.>

Ith particle state variable of t-1 th working procedure,/- >

State variables characterized by the ith particle for the 1 st to kth operation under unknown parameter conditions are related to the battery capacity data for the 1 st to kth operation,/h>

For the relation of the state variables characterized by the ith particle of the 1 st to kth working process under the condition of unknown parameters, +.>

For the relation between the unknown parameter conditions and the battery capacity data from the 1 st to the kth operation under the condition of the state variable characterized by the i particles from the 1 st to the kth operation, ">

Relation between unknown parameter conditions and ith particle state variable of the t-1 th working process under the condition of the ith particle state variable of the t-1 th working process, +.>

The relation between the unknown parameter condition and the battery capacity data of the ith working process under the condition of the state variable of the ith particle of the ith working process;

based on the contents of the state transition equations and the measurement equations in the adaptive state update model equations (1) to (3) of the UPF in step 2, the relational expression in equations (26) to (27) can be obtained:

wherein ,

is the mean value of +.>

Variance is->

T is the number of the working process, < ->

Noise variance of the ith particle for the t-th working process, +. >

Process variable covariance matrix of ith particle of the t-th working process,/th working process>

diag represents a diagonal matrix, < >>

For the respective t-th working process variable noise term u _a ，u _b ，u _c ，u _d Is a variance of (2);

substituting equation (26) - (27) into equation (25) and ignoring the constant term, the joint log likelihood function may be further expressed as shown in equation (28):

wherein, oc represents a proportional meaning;

battery capacity data y from initial to kth cycle is known _1:k ＝[y ₁ ,y ₂ ,…,y _k ]On the premise that for the mth iteration process of the EM algorithm, the unknown parameter estimation vector of the adaptive degradation model can be expressed as:

wherein ,

is the noise variable u of the kth working process in the mth iteration process _a Variance of->

Is the noise variable u of the kth working process in the mth iteration process _b Variance of->

Is the noise variable u of the kth working process in the mth iteration process _c Variance of->

Is the noise variable u of the kth working process in the mth iteration process _d Variance of->

Measuring the variance of the noise variable for the kth work in the mth iteration process;

derivation according to formulas (25) - (28), and the basic principle of EM algorithm;

for the m+1th iteration process, it can be divided into E step and M step, and can be expressed as shown in formulas (29) — (30):

(1) E, step E: calculation of

(2) M step: calculation of

wherein ,

is of xi->

Joint log likelihood function under conditions, +.>

R is _1:k ,y _1:k At->

Expected operation under conditions; />

Estimating a vector for an unknown parameter of the kth working process for the (m+1) th iteration process, and obtaining an unknown parameter value of the degradation model by iterating the step E and the step M until a certain convergence criterion is met;

according to equation (28), equation (29) can be written expansively in the form shown in equation (31):

in the formula (31) of the present invention,

and

is the ith particle based on battery capacity data y from initial to kth duty cycle _1:k ＝[y ₁ ,y ₂ ,…,y _k ]Conditional expectation of implicit state variables of (2);

in order to calculate the state variable expectation under posterior conditions, the invention adopts a Rauch-Tung-Striebel (RTS) optimal smoothing algorithm to carry out backward smoothing. Based on the forward iteration of the UPF algorithm, an estimate of the expected and covariance of the state variable represented by the particle can be obtained as shown in equations (32) and (33):

wherein ,

for the state variable of the ith particle, expected in RTS algorithm, +.>

Covariance of state variables for the ith particle in RTS algorithm, +.>

For the posterior state variable expectation,/-for the ith particle of the unscented particle filter algorithm during the kth operation >

The posterior covariance matrix of the ith particle in the kth working process of the unscented particle filtering algorithm is obtained;

accordingly, the covariance between the state variables between k-1 and k duty cycles can be expressed as shown in equation (34):

wherein ,

covariance between state variables between k-1 and k duty cycle represented by the ith particle, +.>

UKF gain for the ith particle, < ->

For covariance between state variables between k-1 and k duty cycles，I _4×4 Is a four-dimensional unit array->

A posterior covariance matrix between state variables of k-1 working processes represented by an ith particle in the unscented particle filtering algorithm;

according to the initialization condition in the formula (32) - (34), RTS smoothing operation can be performed; the expression of the RTS smoothing gain is shown in formula (35):

wherein ,

for the i-th particle RTS smoothing gain during the t-th operation,/>

UKF posterior covariance for ith particle during the t-th working process,/>

Predicting covariance for UKF of the ith particle in the t working process in one step;

accordingly, state variables in a backward iterative process

Sum of covariance->

The update of (a) is shown in the following formulas (36) - (37):

wherein ,

for the posterior state vector of the ith particle for the t +1 th course of action in RTS backward smoothing,

Predicted value for the (t+1) th working process state vector in UKF forward prediction for the (i) th particle,>

for the a priori state vector of the ith particle for the t-th working process in RTS backward smoothing, +.>

For the posterior state vector of the ith particle for the t-th working process in RTS backward smoothing, +.>

Covariance for the (t+1) th working process in RTS backward smoothing for the (i) th particle,/>

For the ith particle as predicted value for the t+1th working process covariance in UKF forward prediction, +.>

For the posterior covariance of the ith particle for the t-th working process in UKF forward operation, +.>

Covariance for the ith particle for the t-th course of work in RTS backward smoothing;

in addition, the covariance between the state matrices between two adjacent loops can be expressed as shown in equation (38):

wherein ,

covariance of t-1 th working procedure and t-th working procedure state variable in RTS backward smoothing, +.>

Covariance of state variables of the t-th working process and the t+1th working process in RTS backward smoothing;

based on formulas (32) - (38), the conditional desired expression can be solved

and />

As shown in the formulas (39) - (42):

wherein ,

the state variable represented by the ith particle of the t-1 th working process in RTS backward smoothing;

Substituting the conditional expectation in the equation (39) - (42) into the equation (31) can be written in the form shown in the equation (43)

Based on the calculation result of step E in formulas (31) - (43), step M of EM algorithm is calculated according to formula (30), and vector is estimated for unknown parameters of degradation model

The partial derivative of each parameter is set to 0, and the results of formulas (44) - (45) can be solved. />

By continuously iterating the E step and the M step until the conditions are satisfied

Or stopping iteration until the maximum iteration number is 10 times, thereby realizing self-adaptive estimation on unknown parameters in the degradation model;

wherein epsilon is a parameter set manually.

Other steps and parameters are the same as in the first or second embodiment.

The specific embodiment IV is as follows: the difference between the present embodiment and one to three embodiments is that in the step 4, it is determined whether the capacity regeneration phenomenon occurs in the kth working process;

the invention provides a method based on Wilcoxon rank sum test, which is used for measuring the difference between the prior probability density distribution and the posterior probability density distribution of particles in unscented particle filter (Unscented Kalman Filter, UPF), so that the regeneration point of the battery capacity is accurately detected, and the residual service life prediction error caused by the sudden increase of the battery capacity is reduced;

The specific process is as follows:

in the adaptive state update model based on UPF in step 2, the posterior probability density distribution of particles in the UKF algorithm can be obtained through formulas (1) - (19), thereby serving as the prior probability density distribution of particle filtering. Next, a posterior probability density distribution of the particle filtered can be obtained by the formulas (20) - (24). For convenience of expression, in the battery capacity updated by the UKF algorithm, specifically, the state variable solved by the formula (18) is substituted into the formula (3) to obtain the battery capacity set as

In detail, the state variable solved by the formula (23) is substituted into the formula (3) to obtain the battery capacity set of +.>

Since the battery will undergo capacity regeneration after standing for a period of time, this will result in

and />

The distribution between them produces a certain distribution difference. The effect of the Wilcoxon rank sum test is to measure the difference in distribution between two sets of data samples. This has the advantage that no special assumptions need be made about the distribution of the data samples, so that it can be adapted to some complex data distribution situations. For this purpose, the invention will measure +. >

and />

The distribution difference between them, let us assume H ₀ The following is shown:

suppose H ₀ : two totalities

and />

The same probability density distribution is provided at the significance level alpha;

to test hypothesis H ₀ The Wilcoxon rank sum test can be summarized as the following steps:

(1) Combining two populations

and />

Mixing to obtain the total Q _k Uniformly numbering according to the values of the battery capacities from small to large, and defining the ordinal number corresponding to each battery capacity sample as a rank;

(2) Calculating sample population

A sum T of ranks corresponding to the samples of (2);

(3) According to the total number N of particles _s And significance level gamma consults a rank sum check table, a rank sum lower limit T can be obtained ₁ Rank and upper limit T ₂ And a corresponding confidence level p;

(4) Set interval [ T ] ₁ ,T ₂ ]If T is within the interval, i.e. p>Gamma, then accept hypothesis H ₀ Consider two general populations

and />

No significant difference at level α; otherwise reject hypothesis H ₀ Consider two general->

and />

There is a significant difference at level α;

summarizing, once H is assumed ₀ Is rejected, then explain

and />

There is a large difference in the distribution of (a), which means that the battery generates a phenomenon of capacity regeneration in the kth cycle.

Other steps and parameters are the same as in one to three embodiments.

Fifth embodiment: the difference between the embodiment and the specific embodiment is that in the step 5, the remaining service life of the battery is solved through trend extrapolation according to the battery degradation model constructed in the step 2 and the current battery capacity value;

calculating a confidence interval of the battery capacity according to the battery degradation model constructed in the step 2 and the current battery capacity value;

the specific process is as follows:

the predicted value of the residual service life of the battery is the number of cycles passed when the capacity is degraded to the failure threshold value for the first time;

according to the state variable value of the battery degradation model in step 2 and the current battery capacity value y _k The residual service life of the battery can be solved by a trend extrapolation mode, and the specific process is as follows:

1) Let the predicted starting point v=k, bring the equation (23) state variable expectation (expectation being mean) into the equationA kind of electronic device with high-pressure air-conditioning system

Calculating the battery capacity;

if the phenomenon that the capacity of the battery is regenerated in the kth working cycle is detected, replacing the predicted starting point v=k with v=k-1;

2) If the calculated battery capacity y _v Less than or equal to the failure threshold (failure threshold set), then execution 3);

if the calculated battery capacity y _v Greater than the failure threshold, let v=v+1, by the formula

Calculating battery capacity, repeating 2) until y is calculated _v A failure threshold value of the battery capacity or less;

3) Obtaining the remaining service life of the battery through a formula (46):

RUL _k ＝v–k (46)

4) If the calculated battery capacity y _v Greater than failure threshold, obtain N _s And the average value and the standard deviation of the individual particles are added and subtracted by plus and minus three times of the standard deviation of the average value to be used as the upper line and the lower line of the confidence interval of the battery capacity.

Assume that battery capacity data y of the battery in the kth process is obtained _k Then through the above algorithm steps of UPF, EM, etc., the battery capacity prediction from the kth point is started, and the kth working cycle is the starting point of the prediction. If it is detected that the capacity regeneration phenomenon of the battery occurs in the kth duty cycle, a predicted starting point y is required _k Is replaced by y _k-1 To reduce prediction errors.

Other steps and parameters are the same as in one to four embodiments.

Evaluating the on-line prediction effect of the remaining service life: the effect of the battery remaining service life prediction method proposed by the invention is evaluated by adopting two indexes of root mean square error (Root Mean Square Error, RMSE) and absolute average error (Mean Absolute Error, MAE). The expression of the root mean square error and the absolute average error sum is shown in the formulas (47) - (48):

Wherein N is the number of samples of the test data, a is the serial number of the samples, RUL _pa and RUL_ta Respectively, a predicted value and a true value of the residual service life of the a sample. The smaller the root mean square error and the absolute average error are, the better the residual service life prediction method provided by the invention is.

The following examples are used to verify the benefits of the present invention:

embodiment one:

the invention adopts a battery data set provided by a NASA Ames prediction center electricity to verify the proposed residual service life prediction method based on the EM-UPF-W. The data set contains data generated by operating model 18650 of 4 lithium ion batteries, numbered B0005, B0006, B0007, and B0018, at 24 ℃. The 4 battery operation modes with rated capacity of 2Ah include a charge mode and a discharge mode. For the charging process, the battery was charged at a constant current of 1.5A until the battery voltage rose to 4.2V, and then turned into constant voltage charging until the charging current dropped to 20mA. For the discharge mode, the battery is discharged at a constant 2A current until the voltage across the 4 batteries drops to different thresholds, respectively. On this basis, the NASA battery dataset provides battery capacity at the end of each charge-discharge cycle and thus can be used directly for dynamic prediction of the remaining useful life of the battery. The degradation thresholds of these 4 cells were set to 1.4Ah,1.22Ah,1.6Ah and 1.4Ah, respectively, considering the different operating conditions of the different cells. The battery capacity fade curves for these 4 cells are shown in fig. 2.

Step 1, extracting battery capacity data in the kth working process: the NASA battery data set provides the battery capacity at the end of each charge-discharge cycle, whereby battery capacity data at the corresponding cycle can be obtained as input to the battery remaining life prediction method.

Step 2, constructing a degradation model of the battery: and (3) constructing a dynamic degradation model based on unscented particle filtering to describe the degradation process of the battery only according to the battery capacity data in the step (1). Implicit parameter variables are obtained by steps of particle initialization, estimation of UKF state variables and covariance, particle weight updating, resampling, and the like.

Step 3, adaptively estimating process noise and measurement noise in the degradation model: process noise and measurement noise in the degradation model are adaptively estimated using a maximum expectation-based algorithm. Specifically, for process noise and measurement noise in a dynamic degradation model based on unscented particle filtering, a maximum expected algorithm is adopted to perform adaptive estimation, and all capacity data of a battery from an initial state to a current moment are used for recursively updating the process noise and the measurement noise. The present invention takes the battery numbered B0005 as an example to show the results of adaptive parameter estimation of the algorithm for process noise and measurement noise, as shown in fig. 3a, 3B, 3c, 3d, 3 e. As can be seen from fig. 3a, 3b, 3c, 3d, 3e, all model parameters can achieve convergence within 20 working cycles, so the parameter estimation method based on the maximum expectation provided by the invention has a good convergence effect.

Step 4, judging whether capacity regeneration phenomenon occurs in the kth working process: the difference between the prior probability density distribution and the posterior probability density distribution of particles in unscented particle filtering is measured by adopting a Wilcoxon rank sum test-based method, so that the regeneration point of the battery capacity is accurately detected, and the residual service life prediction error caused by sudden increase of the battery capacity is reduced. Fig. 4a, 4B of the present invention illustrate the results of the Wilcoxon rank sum test method on battery capacity regeneration points, taking the battery numbered B0005 as an example.

Step 5, predicting the remaining service life of the battery: based on the battery degradation model and the current battery capacity value in step 2, switching onAnd (5) extrapolating the over-trend to solve the residual service life of the battery. It should be noted that if the phenomenon of capacity regeneration of the battery in the kth duty cycle is detected, it is necessary to predict the starting point y _k Is replaced by y _k-1 To reduce prediction errors. Finally, each particle is respectively substituted into the model, so that the confidence interval of the battery capacity can be obtained.

Fig. 5a, 5b, 5c, 5d show the actual and predicted values of the remaining life of 4 cells on the NASA cell dataset for the proposed method. It can be seen that the predicted value and the true value of the residual service life in the method provided by the invention are very close, so that the method provided by the invention has a very good prediction effect.

Fig. 6a, 6b, 6c, 6d show confidence intervals of 99% of the battery capacity, and it can be seen that the true values of the battery capacity are basically within the confidence intervals, so that the method proposed by the invention can well describe the uncertainty of the battery in the degradation process.

Step 6, evaluating the online prediction effect of the residual service life: the effect of the battery remaining service life prediction method proposed by the invention is evaluated by adopting two indexes of root mean square error (Root Mean Square Error, RMSE) and absolute average error (Mean Absolute Error, MAE).

TABLE 1 prediction of remaining useful life

The present invention is capable of other and further embodiments and its several details are capable of modification and variation in light of the present invention, as will be apparent to those skilled in the art, without departing from the spirit and scope of the invention as defined in the appended claims.

Claims (2)

1. The battery remaining service life prediction method based on maximum expected-unscented particle filtering is characterized by comprising the following steps of: the method comprises the following specific processes:

step 1, extracting battery capacity data in the kth working process;

step 2, constructing a dynamic battery degradation model based on unscented particle filtering according to the battery capacity data in the kth working process extracted in the step 1;

Step 3, adaptively estimating process noise and measurement noise in the battery degradation model by adopting a maximum expectation-based algorithm;

the process noise and the measurement noise are used for constructing a dynamic battery degradation model based on unscented particle filtering in the k+1th working process;

step 4, judging whether capacity regeneration occurs in the kth working process;

step 5, solving the residual service life of the battery through trend extrapolation according to the battery degradation model constructed in the step 2 and the current battery capacity value;

calculating a confidence interval of the battery capacity according to the battery degradation model constructed in the step 2 and the current battery capacity value;

in the step 2, a dynamic battery degradation model based on unscented particle filtering is constructed according to the battery capacity data in the kth working process extracted in the step 1; the specific process is as follows:

the battery degradation model based on unscented particle filtering can be expressed in the form of a state space as shown in formula (1):

wherein ,x_k Is the state variable at the kth working cycle, x _k ＝[a _k ,b _k ,c _k ,d _k ] ^T ，a _k A first component, b, being a state variable at the kth duty cycle _k A second component, c, being a state variable at the kth duty cycle _k D, being the third component of the state variable at the kth duty cycle _k A fourth component of the state variable at the kth duty cycle, T being the transpose, k representing the number of battery duty cycles, x _k-1 Is the state variable at the k-1 th working cycle, y _k For battery capacity, u _k U is a process noise term _k ＝[u _a ,u _b ,u _c ,u _d ] ^T ，u _a Noise term, u, being the first component of the state variable at the kth duty cycle _b Noise term, u, being the second component of the state variable at the kth duty cycle _c A noise term which is the third component of the state variable at the kth duty cycle, u _d A noise term, v, which is the fourth component of the state variable at the kth duty cycle _k To measure noise term, v _k ∈R ^1×1 1×1 is a matrix with a length and width of 1, f () represents a mapping relation of a state transition equation, and h () represents a mapping relation of a measurement equation;

state transition equation x in equation (1) _k ＝f(x _k-1 )+u _k Can be written in the form of equation (2):

wherein ,

for the process variable noise term u _a Variance of->

For the process variable noise term u _b Variance of->

For the process variable noise term u _c Variance of->

For the process variable noise term u _d Is a variance of (2); a, a _k-1 A first component, b, being a state variable at the kth-1 duty cycle _k-1 A second component, c, being a state variable at the kth-1 duty cycle _k-1 A third component, d, being a state variable at the (k-1) th duty cycle _k-1 A fourth component that is a state variable at the kth-1 duty cycle; />

Represents u _a Obeying the mean value to be 0, variance to be +.>

Normal distribution of->

Mean value 0, variance +.>

Normal distribution of->

Mean value 0, variance +.>

Normal distribution of->

Is 0 as mean and 0 as variance->

Normal distribution of->

Mean value 0, variance +.>

Is a normal distribution of (2);

measurement equation y _k ＝h(x _k )+v _k Can be written in the form of equation (3):

wherein ,

for measuring the variance of the noise; />

Mean value 0, variance +.>

Is a normal distribution of (2);

(1) Particle initialization, the process is:

according to the initial state variable x ₀ In the initial state variable x ₀ Normal distribution p (x) ₀ ) In generating particle sets by the Monte Carlo method

N _s Is the total number of particles, and the initialization weight value of each particle is +.>

Initializing the mean of particles

Initializing the variance S of particles ₀ Can be written in the form shown in the formulas (4) - (5):

wherein i is the number of the particle,

to initialize the mean value of the particles, S ₀ To initialize the variance of the particles, E () represents the desired operator, T represents the matrix transpose operator, r ₀ To initialize the particles;

(2) Calculating state variable r characterized by particles through UKF algorithm _k Posterior expectation of (2)

Sum of covariance- >

The process is as follows:

after completing the k-1 iterative process, a Sigma sampling point set of each particle is constructed

As shown in formula (6):

wherein ,

for the mean vector of the ith particle during k-1 work>

The covariance matrix of the ith particle in the k-1 working process is adopted, and the parameter alpha is a coefficient value in UKF;

the number of sample points of the Sigma point set of each particle is 2n+1, n is the number of state variables, and the weight expressions of the sample points of the Sigma point set of each particle are shown in formulas (7) - (9):

wherein ,

process variable weights represented by sample points of the ith particle, sigma Point set 1, for example>

Weight of the measurement variable represented by the sample point of the 1 st Sigma point set for the ith particle, for>

Process variable weights represented by sample points of the jth Sigma point set for the ith particle, +.>

The weight of a measured variable represented by a sample point of the jth Sigma point set of the ith particle, the parameter beta, lambda are coefficient values in UKF, and n is the number of state variables;

on this basis, the state estimation expression of the UKF can be written in the form shown in formulas (10) - (11):

/>

wherein ,

for the prior state vector represented by the j Sigma point in the ith particle of the kth working process, f () is the mapping relation of the state transition equation, +. >

For the posterior state vector represented by the j-th Sigma point in the ith particle of the kth-1 th course of operation,/->

For the noise variable represented by the ith particle of the kth-1 working process +.>

A weighted prior state vector for the ith particle of the kth course of operation;

the expression for covariance estimation can then be written in the form shown in equation (12):

wherein ,

a priori covariance of the state variables in the ith particle for the kth course of operation, +.>

Covariance matrix of noise variable in ith particle for kth working process of ith particle, +.>

diag represents diagonal momentArray (S)>

For the kth process variable noise term u _a Variance of->

For the kth process variable noise term u _b Variance of->

Process variable noise term u for the kth course of operation _c Variance of->

For the kth process variable noise term u _d Is a variance of (2);

for forward one-step prediction, there are forms shown in formulas (13) - (14):

wherein ,

for the state vector represented by the jth Sigma point in the ith particle, +.>

For UKF forward one-step prediction represented by the jth Sigma point in the ith particle,/I>

U represented by the ith particleKF predictive value,/-A>

For the measurement noise item of the ith particle in the kth working process, h () represents the mapping relation of the measurement equation;

Calculating an autocovariance matrix

And cross covariance matrix->

As shown in the formulas (15) - (16):

wherein ,

noise variance of the ith particle for the kth course of operation;

thereby, a Kalman filtering gain vector K can be obtained _k As shown in formula (17):

wherein ,

for the ith particle Kalman filter gain vector, -/->

As auto-covariance momentArray (S)>

As for the cross covariance matrix, with the addition of new battery capacity data in the kth iteration, the state update and covariance update of the UKF can be expressed as the forms shown in formulas (18) - (19):

/>

wherein ,

to calculate the state variable r of a particle by the UKF algorithm _k Posterior expectation of->

To calculate the state variable r of a particle by the UKF algorithm _k Is a covariance of (2);

according to the results in the formulas (18) - (19), for

Obtaining the mean value->

For->

Obtaining the mean value->

Based on->

and />

Obtaining posterior probability Density distribution of particles under UKF Algorithm->

(3) The particle weight is updated by the following steps:

the importance density function of the particles obtained according to (2) can be expressed as:

wherein ,

importance density function representing the ith particle of the kth course of operation,/for>

For the state variables represented by the ith particle initial state through the kth-1 working process, y _1:k Battery capacity data indicating the 1 st to kth operation, +.>

State variables represented for the ith particle initial state to the kth-1 operation and battery capacity data prior conditions for the 1 st to kth operation +.>

Probability density function of>

Is the mean value of +.>

Variance is->

Normal distribution of (a) -is the meaning obeying a certain distribution;

on this basis, the weight update expression of the particle can be written as in formula (20):

wherein ,

for the weight of the ith particle under k-1 duty cycles,/for the particle>

Weight of the ith particle under k duty cycles, +.>

Is y _k At->

Probability density function under conditions->

Is->

At->

Probability density function under conditions;

after the weight normalization of the particles, normalized particle weights can be obtained>

The expression is as follows:

(4) The particle resampling process comprises the following steps:

let the effective particle number be N _ef Can be defined as in equation (22):

setting a threshold value N _th If the effective particle number N _ef Less than N _th Resampling of the particles is required;

a resampling process of particle filtering is realized by selecting a random resampling mode;

if the effective particle number N _ef N is greater than or equal to _th Then no resampling of the particles is required;

after completion of (4), each particle has a weight of

(5) The state updating process comprises the following steps:

based on the weight of the particles obtained by resampling, the expected state variable under the posterior condition of the unscented particle filter algorithm can be obtained through formulas (23) - (24)

Sum of covariance->

wherein ,

for state variables under the posterior condition of the UPF algorithm, it is desirable to have, < +.>

Covariance under the posterior condition of UPF algorithm;

in the step 3, adaptively estimating process noise and measurement noise in a battery degradation model based on a maximum expected algorithm; the specific process is as follows:

the unknown parameter vector composed of the process noise and the measurement noise can be expressed as:

battery capacity data y from initial to kth cycle is known _1:k ＝[y ₁ ,y ₂ ,…,y _k ]On the premise of (1) can construct a joint log-likelihood function L _k (xi) is as shown in formula (25):

wherein ,r_t ⁽ⁱ⁾ As a state variable of the ith particle of the t-th working process,

n, which is the value of the process variable in the t-th working process of the ith particle _s Is the total number of particles, and the initialization weight value of each particle is +.>

i is the number of the particle, ">

State variables characterized by the ith particle for the 1 st to kth course of operation, +. >

A joint probability density function composed of state variables represented by the ith particles from the 1 st working procedure to the kth working procedure under the unknown parameter condition and battery capacity data from the 1 st working procedure to the kth working procedure,>

a conditional probability density function for a state variable characterized by the ith particle of the 1 st to kth operation under unknown parameter conditions, +.>

Conditional probability density function of battery capacity data for 1 st to kth operation under conditions of unknown parameters and state variables characterized by the 1 st to kth operation's ith particles, +.>

A conditional probability density function of the ith particle state variable of the t-1 th working process under the unknown parameter condition and the ith particle state variable of the t-1 th working process,>

a conditional probability density function for the battery capacity data of the ith operating process under the unknown parameter conditions and the state variable of the ith particle of the ith operating process,

ith particle state variable of t-1 th working procedure,/->

State variables characterized by the ith particle for the 1 st to kth operation under unknown parameter conditions are related to the battery capacity data for the 1 st to kth operation,/h >

For the relation of the state variables characterized by the ith particle of the 1 st to kth working process under the condition of unknown parameters, +.>

For the relationship of the unknown parameter conditions and the battery capacity data from the 1 st to the kth operation under the condition of the state variable characterized by the i-th particles from the 1 st to the kth operation,

relationship between unknown parameter condition and ith particle state variable of the t-1 th working process under the condition of the ith particle state variable of the t-1 th working process, y _t |r _t ⁽ⁱ⁾ The xi is the relation between the unknown parameter condition and the battery capacity data of the t-th working process under the condition of the i-th particle state variable of the t-th working process;

based on the contents of the state transition equation and the measurement equation in the formulas (1) — (3), the relational expression in the formulas (26) — (27) can be obtained:

wherein ,

is the mean value of +.>

Variance is->

T is the number of the working process, < ->

Noise variance of the ith particle for the t-th working process, +.>

Process variable covariance matrix of ith particle of the t-th working process,/th working process>

diag represents a diagonal matrix, < >>

For the respective t-th working process variable noise term u _a ，u _b ，u _c ，u _d Is a variance of (2);

substituting equation (26) - (27) into equation (25) and ignoring the constant term, the joint log likelihood function may be further expressed as shown in equation (28):

Wherein, oc represents a proportional meaning;

battery capacity data y from initial to kth cycle is known _1:k ＝[y ₁ ,y ₂ ,…,y _k ]On the premise that for the mth iteration process of the EM algorithm, the unknown parameter estimation vector of the adaptive degradation model can be expressed as:

wherein ,

is the noise variable u of the kth working process in the mth iteration process _a Variance of->

Is the noise variable u of the kth working process in the mth iteration process _b Variance of->

Is the noise variable u of the kth working process in the mth iteration process _c Variance of->

Is the noise variable u of the kth working process in the mth iteration process _d Variance of->

Measuring the variance of the noise variable for the kth work in the mth iteration process;

for the m+1th iteration process, it can be divided into E step and M step, and can be expressed as shown in formulas (29) — (30):

(1) E, step E: calculation of

(2) M step: calculation of

wherein ,

is of xi->

Joint log likelihood function under conditions, +.>

R is _1:k ,y _1:k At the position of

Expected operation under conditions; />

Estimating a vector for unknown parameters of the kth working process for the (m+1) th iteration process;

according to equation (28), equation (29) can be written expansively in the form shown in equation (31):

in the formula (31) of the present invention,

and

is the ith particle based on battery capacity data y from initial to kth duty cycle _1:k ＝[y ₁ ,y ₂ ,…,y _k ]Conditional expectation of implicit state variables of (2);

based on the forward iteration of the UPF algorithm, an estimate of the expected and covariance of the state variable represented by the particle can be obtained as shown in equations (32) and (33):

/>

wherein ,

for the state variable of the ith particle, expected in RTS algorithm, +.>

Covariance of state variables for the ith particle in RTS algorithm, +.>

For the posterior state variable expectation,/-for the ith particle of the unscented particle filter algorithm during the kth operation>

The posterior covariance matrix of the ith particle in the kth working process of the unscented particle filtering algorithm is obtained;

accordingly, the covariance between the state variables between k-1 and k duty cycles can be expressed as shown in equation (34):

wherein ,

for the covariance between the state variables between k-1 and k duty cycles represented by the ith particle,

UKF gain for the ith particle, < ->

For covariance between state variables between k-1 and k duty cycles, I _4×4 Is a four-dimensional unit array->

A posterior covariance matrix between state variables of k-1 working processes represented by an ith particle in the unscented particle filtering algorithm;

according to the initialization condition in the formula (32) - (34), RTS smoothing operation can be performed; the expression of the RTS smoothing gain is shown in formula (35):

wherein ,

for the i-th particle RTS smoothing gain during the t-th operation,/>

UKF posterior covariance for ith particle during the t-th working process,/>

Predicting covariance for UKF of the ith particle in the t working process in one step;

accordingly, state variables in a backward iterative process

Sum of covariance->

The update of (a) is shown in the following formulas (36) - (37):

wherein ,

for the posterior state vector of the ith particle for the t+1th working process in RTS backward smoothing, +.>

Predicted value for the (t+1) th working process state vector in UKF forward prediction for the (i) th particle,>

for the a priori state vector of the ith particle for the t-th working process in RTS backward smoothing, +.>

For the posterior state vector of the ith particle for the t-th working process in RTS backward smoothing, +.>

Covariance for the (t+1) th working process in RTS backward smoothing for the (i) th particle,/>

For the ith particle as predicted value for the t+1th working process covariance in UKF forward prediction, +.>

For the posterior covariance of the ith particle for the t-th working process in UKF forward operation, +.>

Covariance for the ith particle for the t-th course of work in RTS backward smoothing;

in addition, the covariance between the state matrices between two adjacent loops can be expressed as shown in equation (38):

wherein ,

for covariance of the t-1 st working process and the t-th working process state variable in RTS backward smoothing,

covariance of state variables of the t-th working process and the t+1th working process in RTS backward smoothing;

based on formulas (32) - (38), the conditional desired expression can be solved

and />

As shown in the formulas (39) - (42):

wherein ,

the state variable represented by the ith particle of the t-1 th working process in RTS backward smoothing;

substituting the conditional expectation in the equation (39) - (42) into the equation (31) can be written in the form shown in the equation (43)

Based on the calculation result of step E in formulas (31) - (43), step M of EM algorithm is calculated according to formula (30), and vector is estimated for unknown parameters of degradation model

The partial derivative of each parameter of the formula (44) to (45) can be solved by making the value of the partial derivative 0Fruit;

by continuously iterating the E step and the M step until the conditions are satisfied

Or stopping iteration until the maximum iteration number is 10 times, thereby realizing self-adaptive estimation on unknown parameters in the degradation model;

wherein epsilon is a manually set parameter;

judging whether capacity regeneration occurs in the kth working process in the step 4; the specific process is as follows:

The battery capacity set updated by UKF algorithm is that

The set of battery capacities after update by UPF algorithm is +.>

Suppose H ₀ : two totalities

and />

The same probability density distribution is provided at the significance level alpha;

to test hypothesis H ₀ The Wilcoxon rank sum test can be summarized as the following steps:

(1) Combining two populations

and />

Mixing to obtain the total Q _k Uniformly numbering according to the values of the battery capacities from small to large, and defining the ordinal number corresponding to each battery capacity sample as a rank;

(2) Calculating sample population

A sum T of ranks corresponding to the samples of (2);

(3) According to the total number N of particles _s And significance level gamma consults a rank sum check table, a rank sum lower limit T can be obtained ₁ Rank and upper limit T ₂ And a corresponding confidence level p;

(4) Set interval [ T ] ₁ ,T ₂ ]If T is within the interval, i.e. p>Gamma, then accept hypothesis H ₀ Consider two general populations

and />

No significant difference at level α; otherwise reject hypothesis H ₀ Consider two general->

and />

There is a significant difference at level α;

meaning that the battery has a capacity regeneration phenomenon in the kth cycle.

2. The method for predicting remaining life of a battery based on maximum expected-unscented particle filtering of claim 1, wherein: in the step 5, according to the battery degradation model constructed in the step 2 and the current battery capacity value, solving the residual service life of the battery through trend extrapolation;

Calculating a confidence interval of the battery capacity according to the battery degradation model constructed in the step 2 and the current battery capacity value;

the specific process is as follows:

1) Let the predicted starting point v=k, take the equation (23) state variable into the equation as desired

Calculating the battery capacity; />

If the phenomenon that the capacity of the battery is regenerated in the kth working cycle is detected, replacing the predicted starting point v=k with v=k-1;

2) If the calculated battery capacity y _v And (3) if the failure threshold value is less than or equal to the failure threshold value, executing the step (3);

if the calculated battery capacity y _v Greater than the failure threshold, let v=v+1, by the formula

Calculating battery capacity, repeating 2) until y is calculated _v A failure threshold value of the battery capacity or less;

3) Obtaining the remaining service life of the battery through a formula (46):

RUL _k ＝v–k (46)

4) If the calculated battery capacity y _v Greater than failure threshold, obtain N _s And the average value and the standard deviation of the individual particles are added and subtracted by plus and minus three times of the standard deviation of the average value to be used as the upper line and the lower line of the confidence interval of the battery capacity.

Images