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

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

Info

Publication number
CN114779088B
CN114779088B CN202210415333.5A CN202210415333A CN114779088B CN 114779088 B CN114779088 B CN 114779088B CN 202210415333 A CN202210415333 A CN 202210415333A CN 114779088 B CN114779088 B CN 114779088B
Authority
CN
China
Prior art keywords
particle
working process
battery
kth
state
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210415333.5A
Other languages
Chinese (zh)
Other versions
CN114779088A (en
Inventor
张九思
李翔
罗浩
蒋宇辰
吴诗梦
田纪伦
尹珅
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology Shenzhen
Original Assignee
Harbin Institute of Technology Shenzhen
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology Shenzhen filed Critical Harbin Institute of Technology Shenzhen
Priority to CN202210415333.5A priority Critical patent/CN114779088B/en
Publication of CN114779088A publication Critical patent/CN114779088A/en
Application granted granted Critical
Publication of CN114779088B publication Critical patent/CN114779088B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R31/00—Arrangements for testing electric properties; Arrangements for locating electric faults; Arrangements for electrical testing characterised by what is being tested not provided for elsewhere
    • G01R31/36—Arrangements for testing, measuring or monitoring the electrical condition of accumulators or electric batteries, e.g. capacity or state of charge [SoC]
    • G01R31/367—Software therefor, e.g. for battery testing using modelling or look-up tables
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R31/00—Arrangements for testing electric properties; Arrangements for locating electric faults; Arrangements for electrical testing characterised by what is being tested not provided for elsewhere
    • G01R31/36—Arrangements for testing, measuring or monitoring the electrical condition of accumulators or electric batteries, e.g. capacity or state of charge [SoC]
    • G01R31/392—Determining battery ageing or deterioration, e.g. state of health

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Secondary Cells (AREA)
  • Tests Of Electric Status Of Batteries (AREA)

Abstract

基于最大期望‑无迹粒子滤波的电池剩余使用寿命预测方法，本发明涉及电池剩余使用寿命预测方法。本发明的目的是为了解决现有电池容量突然的增加将对电池的剩余使用寿命预测产生很大的误差的问题。过程为：1、提取第k次工作过程中的电池容量数据；2、构建基于无迹粒子滤波的动态电池退化模型；3、采用基于最大期望算法自适应地估计电池退化模型中的过程噪声和测量噪声；过程噪声和测量噪声用于第k+1次工作过程中，构建基于无迹粒子滤波的动态电池退化模型；4、判断第k次工作过程是否发生容量再生现象；5、求解电池剩余使用寿命和电池容量的置信区间；本发明用于电池的剩余使用寿命预测领域。

Figure 202210415333

A method for predicting the remaining service life of a battery based on maximum expectation-unscented particle filtering, and the invention relates to a method for predicting the remaining service life of a battery. The purpose of the present invention is to solve the problem that a sudden increase in the capacity of the existing battery will cause a large error in the prediction of the remaining service life of the battery. The process is: 1. Extract the battery capacity data in the kth working process; 2. Construct a dynamic battery degradation model based on unscented particle filter; 3. Use the maximum expectation algorithm to adaptively estimate the process noise and Measurement noise; process noise and measurement noise are used in the k+1th working process to construct a dynamic battery degradation model based on unscented particle filter; 4. Determine whether capacity regeneration occurs in the kth working process; 5. Solve the battery remaining Confidence intervals of service life and battery capacity; the invention is used in the field of remaining service life prediction of batteries.

Figure 202210415333

Description

基于最大期望-无迹粒子滤波的电池剩余使用寿命预测方法A method for predicting remaining useful life of batteries based on maximum expectation-unscented particle filtering

技术领域Technical Field

本发明涉及状态空间理论、统计分析和锂离子电池的剩余使用寿命预测相结合的学科交叉领域，具体涉及电池剩余使用寿命预测方法。The present invention relates to an interdisciplinary field combining state space theory, statistical analysis and prediction of the remaining service life of lithium-ion batteries, and in particular to a method for predicting the remaining service life of batteries.

背景技术Background Art

作为重要的储能设备，锂离子电池以其能量密度高、工作温度范围宽、循环寿命长等优势成为电动汽车的重要部件。随着锂离子电池充放电次数的增加，电池的容量将会逐渐下降，这将大大降低系统或者设备的可靠性，甚至在严重时容易发生灾难性事故。因此，准确地预测锂离子电池的剩余使用寿命(Remaining useful Life，RUL)，换句话说，准确地估计电池从当前时刻到完全失效时所经历的时间对于系统的安全性和可靠性具有至关重要的作用。其意义在于能够减少用电设备的维修成本，降低故障发生的概率，实现视情维修。As an important energy storage device, lithium-ion batteries have become an important component of electric vehicles due to their high energy density, wide operating temperature range, and long cycle life. As the number of times lithium-ion batteries are charged and discharged increases, the capacity of the battery will gradually decrease, which will greatly reduce the reliability of the system or equipment, and even in serious cases, catastrophic accidents may occur. Therefore, accurately predicting the remaining useful life (RUL) of lithium-ion batteries, in other words, accurately estimating the time from the current moment to the complete failure of the battery, plays a vital role in the safety and reliability of the system. Its significance lies in reducing the maintenance cost of electrical equipment, reducing the probability of failure, and realizing condition-based maintenance.

目前锂离子电池的剩余使用寿命预测方法从总体上可以分为基于物理失效模型的方法、数据驱动的方法和基于混合经验模型的方法。由于精确的电池机理建模存在着很大的困难，因此基于物理失效模型的方法存在着一定的局限性。另一方面，数据驱动的方法无需事先了解电池的物理化学知识，通过电池工作过程中产生的数据便可以预测剩余使用寿命。然而，数据驱动的学习类方法要求在线数据还被要求遵循与历史数据相似的分布，否则学习类的方法很难保证预测的准确性。基于混合经验模型的方法具备很强的解释性和可靠性，通过将先验知识融入到模型中，构建电池健康因子和电流、电压以及循环次数等参数之间的算法模型，然后采用扩展卡尔曼滤波(Extended Kalman Filter，EKF)、无迹卡尔曼滤波(Unscented Kalman Filter，UKF)以及粒子滤波(Particle Filter，PF)等算法进行参数辨识，从而完成剩余使用寿命预测任务。At present, the remaining service life prediction methods of lithium-ion batteries can be generally divided into methods based on physical failure models, data-driven methods, and methods based on hybrid empirical models. Since accurate battery mechanism modeling is very difficult, the methods based on physical failure models have certain limitations. On the other hand, data-driven methods do not require prior knowledge of the physical and chemical knowledge of the battery, and the remaining service life can be predicted through the data generated during the battery operation process. However, data-driven learning methods require that online data also follow a distribution similar to that of historical data, otherwise it is difficult for learning methods to ensure the accuracy of prediction. The method based on hybrid empirical models has strong interpretability and reliability. By incorporating prior knowledge into the model, an algorithm model between battery health factors and parameters such as current, voltage, and number of cycles is constructed, and then algorithms such as extended Kalman filter (Extended Kalman Filter, EKF), unscented Kalman filter (Unscented Kalman Filter, UKF) and particle filter (Particle Filter, PF) are used for parameter identification to complete the remaining service life prediction task.

虽然数据驱动的方法在电池的剩余使用寿命预测领域得以广泛的运用，但是这些学习类的方法在本质上是黑箱模型，需要大量的历史数据训练算法模型，并且不能量化电池在退化过程中的不确定性。值得注意的是，从总体上而言，锂离子电池的容量随着充放电循环次数的增加是逐渐降低的。然而，当电池处于静置状态时，电池中的电化学反应产物会减少，电化学性能将发生短暂的恢复，这便是电池容量的再生现象。电池容量突然的增加将对电池的剩余使用寿命预测产生很大的误差。因此，准确判断电池容量的再生点至关重要。与此同时，现有的一些基于混合经验模型的方法中模型的过程噪声和测量噪声都需要人为设定，如何采用自适应的方式对噪声进行估计需要深入的研究。Although data-driven methods are widely used in the field of battery remaining service life prediction, these learning methods are essentially black box models that require a large amount of historical data to train the algorithm model and cannot quantify the uncertainty of the battery during degradation. It is worth noting that, in general, the capacity of lithium-ion batteries gradually decreases with the increase in the number of charge and discharge cycles. However, when the battery is in a static state, the electrochemical reaction products in the battery will decrease, and the electrochemical performance will recover briefly, which is the regeneration phenomenon of the battery capacity. The sudden increase in battery capacity will cause a large error in the prediction of the remaining service life of the battery. Therefore, it is crucial to accurately determine the regeneration point of the battery capacity. At the same time, in some existing methods based on hybrid empirical models, the process noise and measurement noise of the model need to be set manually, and how to estimate the noise in an adaptive way requires in-depth research.

发明内容Summary of the invention

本发明的目的是为了解决现有电池容量突然的增加将对电池的剩余使用寿命预测产生很大的误差的问题，而提出基于最大期望-无迹粒子滤波的电池剩余使用寿命预测方法。The purpose of the present invention is to solve the problem that a sudden increase in the existing battery capacity will cause a large error in the prediction of the remaining service life of the battery, and to propose a battery remaining service life prediction method based on maximum expectation-unscented particle filtering.

最大期望-无迹粒子滤波的电池剩余使用寿命预测方法具体过程为：The specific process of the battery remaining service life prediction method based on maximum expectation-unscented particle filtering is as follows:

步骤1、提取第k次工作过程中的电池容量数据；Step 1, extracting the battery capacity data during the kth working process;

步骤2、根据步骤1中提取的第k次工作过程中的电池容量数据，构建基于无迹粒子滤波的动态电池退化模型；Step 2: According to the battery capacity data of the kth working process extracted in step 1, a dynamic battery degradation model based on unscented particle filtering is constructed;

步骤3、采用基于最大期望算法自适应地估计电池退化模型中的过程噪声和测量噪声；Step 3, adaptively estimating process noise and measurement noise in the battery degradation model using a maximum expectation algorithm;

过程噪声和测量噪声用于第k+1次工作过程中，构建基于无迹粒子滤波的动态电池退化模型；The process noise and measurement noise are used to construct a dynamic battery degradation model based on unscented particle filtering during the k+1th working process;

步骤4、判断第k次工作过程是否发生容量再生现象；Step 4, determining whether capacity regeneration occurs during the kth working process;

步骤5、根据步骤2中构建的电池退化模型和当前的电池容量值，通过趋势外推求解电池剩余使用寿命；Step 5: According to the battery degradation model constructed in step 2 and the current battery capacity value, the remaining battery life is solved by trend extrapolation;

根据步骤2中构建的电池退化模型和当前的电池容量值，计算电池容量的置信区间。Based on the battery degradation model constructed in step 2 and the current battery capacity value, calculate the confidence interval of the battery capacity.

本发明的有益效果为：The beneficial effects of the present invention are:

本发明针对单个电池提出了一种新颖的基于无迹粒子滤波-最大期望-Wilcoxon秩和检验的方法预测电池的剩余使用寿命。在无需大量历史数据的前提下，仅仅依据电池容量数据构建基于无迹粒子滤波的动态退化模型描述电池的退化过程。采用基于最大期望算法的自适应地估计动态退化模型中的过程噪声和测量噪声，并通过Wilcoxon秩和检验的方法准确检测电池容量的再生点，进而预测电池的剩余使用寿命，本发明解决了现有电池容量突然的增加将对电池的剩余使用寿命预测产生很大的误差的问题。The present invention proposes a novel method based on untraceable particle filtering-maximum expectation-Wilcoxon rank sum test to predict the remaining service life of a single battery. Without the need for a large amount of historical data, a dynamic degradation model based on untraceable particle filtering is constructed based only on the battery capacity data to describe the degradation process of the battery. The process noise and measurement noise in the dynamic degradation model are adaptively estimated based on the maximum expectation algorithm, and the regeneration point of the battery capacity is accurately detected through the Wilcoxon rank sum test method, thereby predicting the remaining service life of the battery. The present invention solves the problem that a sudden increase in the existing battery capacity will cause a large error in the prediction of the remaining service life of the battery.

1、传统的数据驱动方法需要大量的历史数据训练算法模型，本发明无需有标签的历史数据对于单个电池提出了一种基于无迹粒子滤波的动态退化模型，并提出最大期望算法利用电池工作过程中的容量数据自适应地更新模型过程噪声和测量噪声参数。1. The traditional data-driven method requires a large amount of historical data to train the algorithm model. The present invention does not require labeled historical data and proposes a dynamic degradation model based on unscented particle filtering for a single battery. It also proposes a maximum expectation algorithm that uses the capacity data during the battery operation process to adaptively update the model process noise and measurement noise parameters.

2、对于电池容量再生点，本发明提出了一种基于Wilcoxon秩和检验的方法，对粒子的先验概率密度分布和后验概率密度分布之间的差异进行衡量，由此准确检测电池容量的再生点，减少由于电池容量突增而导致的剩余使用寿命预测误差。2. For the battery capacity regeneration point, the present invention proposes a method based on the Wilcoxon rank sum test to measure the difference between the prior probability density distribution and the posterior probability density distribution of the particles, thereby accurately detecting the regeneration point of the battery capacity and reducing the remaining service life prediction error caused by the sudden increase in battery capacity.

3、考虑到电池在退化过程中的不确定性很难得到量化，本发明提出的预测方法能够通过无迹粒子滤波算法中的粒子的概率分布来确定电池容量的置信区间，由此描述电池退化过程中的不确定性。3. Considering that the uncertainty of the battery in the degradation process is difficult to quantify, the prediction method proposed in the present invention can determine the confidence interval of the battery capacity through the probability distribution of particles in the unscented particle filter algorithm, thereby describing the uncertainty in the battery degradation process.

附图说明BRIEF DESCRIPTION OF THE DRAWINGS

图1是本发明的工作流程图；Fig. 1 is a workflow diagram of the present invention;

图2是NASA Ames预测中心电提供的数据集中4个电池的电池容量衰减曲线示意图；FIG2 is a schematic diagram of battery capacity decay curves of four batteries in a data set provided by NASA Ames Forecast Center;

图3a是本发明中剩余使用寿命预测方法对于B0005电池中退化模型中

Figure BDA0003605634160000031
的自适应估计结果图；FIG. 3a is a diagram of the remaining service life prediction method of the present invention for the degradation model of the B0005 battery.
Figure BDA0003605634160000031
Adaptive estimation result diagram of ;

图3b是本发明中剩余使用寿命预测方法对于B0005电池中退化模型中

Figure BDA0003605634160000032
的自适应估计结果图；FIG. 3b is a diagram of the remaining service life prediction method of the present invention for the degradation model of the B0005 battery.
Figure BDA0003605634160000032
Adaptive estimation result diagram of ;

图3c是本发明中剩余使用寿命预测方法对于B0005电池中退化模型中

Figure BDA0003605634160000033
的自适应估计结果图；FIG. 3c is a diagram of the remaining service life prediction method of the present invention for the degradation model of the B0005 battery.
Figure BDA0003605634160000033
Adaptive estimation result diagram of ;

图3d是本发明中剩余使用寿命预测方法对于B0005电池中退化模型中

Figure BDA0003605634160000034
的自适应估计结果图；FIG. 3d is a diagram of the remaining service life prediction method of the present invention for the degradation model of the B0005 battery.
Figure BDA0003605634160000034
Adaptive estimation result diagram of ;

图3e是本发明中剩余使用寿命预测方法对于B0005电池中退化模型中

Figure BDA0003605634160000035
的自适应估计结果图；FIG. 3e is a diagram of the remaining service life prediction method of the present invention for the degradation model of the B0005 battery.
Figure BDA0003605634160000035
Adaptive estimation result diagram of ;

图4a是采用Wilcoxon秩和检验的方法对于编号为B0005的电池的容量再生点检验的置信度结果图；FIG4a is a graph showing the confidence results of the capacity regeneration point test of the battery numbered B0005 using the Wilcoxon rank sum test method;

图4b是编号为B0005的电池容量再生点的展示图；FIG4 b is a diagram showing a battery capacity regeneration point numbered B0005;

图5a是本发明中剩余使用寿命预测方法对于B0005电池在不同循环下剩余使用寿命预测值和真实值的结果图，其中EM-UPF代表最大期望-无迹粒子滤波算法，EM-UPF-W代表最大期望-无迹粒子滤波-Wilcoxon秩和检验算法，RUL Ground Truth代表剩余使用寿命的真实值；FIG5a is a result diagram of the predicted value and true value of the remaining service life of a B0005 battery under different cycles by the remaining service life prediction method of the present invention, wherein EM-UPF represents the maximum expectation-unscented particle filter algorithm, EM-UPF-W represents the maximum expectation-unscented particle filter-Wilcoxon rank sum test algorithm, and RUL Ground Truth represents the true value of the remaining service life;

图5b是本发明中剩余使用寿命预测方法对于B0006电池在不同循环下剩余使用寿命预测值和真实值的结果图，其中EM-UPF代表最大期望-无迹粒子滤波算法，EM-UPF-W代表最大期望-无迹粒子滤波-Wilcoxon秩和检验算法，RUL Ground Truth代表剩余使用寿命的真实值；FIG5 b is a result diagram of the predicted value and true value of the remaining service life of a B0006 battery under different cycles by the remaining service life prediction method of the present invention, wherein EM-UPF represents the maximum expectation-unscented particle filter algorithm, EM-UPF-W represents the maximum expectation-unscented particle filter-Wilcoxon rank sum test algorithm, and RUL Ground Truth represents the true value of the remaining service life;

图5c是本发明中剩余使用寿命预测方法对于B0007电池在不同循环下剩余使用寿命预测值和真实值的结果图，其中EM-UPF代表最大期望-无迹粒子滤波算法，EM-UPF-W代表最大期望-无迹粒子滤波-Wilcoxon秩和检验算法，RUL Ground Truth代表剩余使用寿命的真实值；FIG5c is a result diagram of the predicted value and true value of the remaining service life of the B0007 battery under different cycles by the remaining service life prediction method of the present invention, wherein EM-UPF represents the maximum expectation-unscented particle filter algorithm, EM-UPF-W represents the maximum expectation-unscented particle filter-Wilcoxon rank sum test algorithm, and RUL Ground Truth represents the true value of the remaining service life;

图5d是本发明中剩余使用寿命预测方法对于B0018电池在不同循环下剩余使用寿命预测值和真实值的结果图，其中EM-UPF代表最大期望-无迹粒子滤波算法，EM-UPF-W代表最大期望-无迹粒子滤波-Wilcoxon秩和检验算法，RUL Ground Truth代表剩余使用寿命的真实值；FIG5 d is a result diagram of the predicted value and true value of the remaining service life of the B0018 battery under different cycles by the remaining service life prediction method of the present invention, wherein EM-UPF represents the maximum expectation-unscented particle filter algorithm, EM-UPF-W represents the maximum expectation-unscented particle filter-Wilcoxon rank sum test algorithm, and RUL Ground Truth represents the true value of the remaining service life;

图6a是本发明中剩余使用寿命预测方法对于B0005在第60个工作循环下的电池容量的99％置信区间图；FIG6 a is a 99% confidence interval diagram of the battery capacity of B0005 at the 60th working cycle according to the remaining service life prediction method of the present invention;

图6b是本发明中剩余使用寿命预测方法对于B0006在第60个工作循环下的电池容量的99％置信区间图；FIG6 b is a 99% confidence interval diagram of the battery capacity of B0006 at the 60th working cycle according to the remaining service life prediction method of the present invention;

图6c是本发明中剩余使用寿命预测方法对于B0007在第60个工作循环下的电池容量的99％置信区间图；FIG6 c is a 99% confidence interval diagram of the battery capacity of B0007 at the 60th working cycle according to the remaining service life prediction method of the present invention;

图6d是本发明中剩余使用寿命预测方法对于B0018在第60个工作循环下的电池容量的99％置信区间图。FIG6 d is a 99% confidence interval diagram of the battery capacity of B0018 at the 60th working cycle according to the remaining service life prediction method of the present invention.

具体实施方式DETAILED DESCRIPTION

具体实施方式一：本实施方式基于最大期望-无迹粒子滤波的电池剩余使用寿命预测方法具体过程为：Specific implementation method 1: The specific process of the battery remaining service life prediction method based on maximum expectation-unscented particle filtering in this implementation method is as follows:

步骤1、提取第k次工作过程中的电池容量数据；Step 1, extracting the battery capacity data during the kth working process;

提取电池在第k次实际工作循环中的容量数据，即可以在电池每个工作循环后进行充放电测试，获得相应循环下的电池容量数据，作为电池剩余使用寿命预测方法的输入。The capacity data of the battery in the kth actual working cycle is extracted, that is, a charge and discharge test can be performed after each working cycle of the battery to obtain the battery capacity data under the corresponding cycle as the input of the battery remaining service life prediction method.

步骤2、根据步骤1中提取的第k次工作过程中的电池容量数据，构建基于无迹粒子滤波的动态电池退化模型，用于描述电池的退化过程；Step 2: According to the battery capacity data of the kth working process extracted in step 1, a dynamic battery degradation model based on unscented particle filtering is constructed to describe the degradation process of the battery;

步骤3、采用基于最大期望算法自适应地估计电池退化模型中的过程噪声和测量噪声；Step 3, adaptively estimating process noise and measurement noise in the battery degradation model using a maximum expectation algorithm;

过程噪声和测量噪声用于第k+1次工作过程中，构建基于无迹粒子滤波的动态电池退化模型；The process noise and measurement noise are used to construct a dynamic battery degradation model based on unscented particle filtering during the k+1th working process;

具体而言，针对基于无迹粒子滤波的动态退化模型中的过程噪声和测量噪声，采用最大期望算法进行自适应估计，利用电池从初始状态到当前时刻的所有容量数据对过程噪声和测量噪声进行递归更新；Specifically, for the process noise and measurement noise in the dynamic degradation model based on unscented particle filtering, the maximum expectation algorithm is used for adaptive estimation, and all the capacity data of the battery from the initial state to the current moment are used to recursively update the process noise and measurement noise;

步骤4、判断第k次工作过程是否发生容量再生现象；Step 4, determining whether capacity regeneration occurs during the kth working process;

本发明提出了一种基于Wilcoxon秩和检验的方法，对无迹粒子滤波(UnscentedKalman Filter，UPF)中粒子的先验概率密度分布和后验概率密度分布之间的差异进行衡量，由此准确检测电池容量的再生点，减少由于电池容量突增而导致的剩余使用寿命预测误差；The present invention proposes a method based on Wilcoxon rank sum test to measure the difference between the prior probability density distribution and the posterior probability density distribution of particles in the unscented Kalman Filter (UPF), thereby accurately detecting the regeneration point of the battery capacity and reducing the remaining service life prediction error caused by the sudden increase of the battery capacity;

步骤5、根据步骤2中构建的电池退化模型和当前的电池容量值，通过趋势外推求解电池剩余使用寿命；Step 5: According to the battery degradation model constructed in step 2 and the current battery capacity value, the remaining battery life is solved by trend extrapolation;

根据步骤2中构建的电池退化模型和当前的电池容量值，计算电池容量的置信区间；评估剩余使用寿命的在线预测效果：采用均方根误差(Root Mean Square Error，RMSE)和绝对平均误差(Mean Absolute Error，MAE)这两个指标评估本发明所提出的电池剩余使用寿命预测方法的效果。According to the battery degradation model constructed in step 2 and the current battery capacity value, the confidence interval of the battery capacity is calculated; the online prediction effect of the remaining service life is evaluated: the root mean square error (RMSE) and the mean absolute error (MAE) are used to evaluate the effect of the battery remaining service life prediction method proposed in the present invention.

具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤2中根据步骤1中提取的第k次工作过程中的电池容量数据，构建基于无迹粒子滤波的动态电池退化模型，用于描述电池的退化过程；具体过程为：Specific implementation method 2: This implementation method is different from specific implementation method 1 in that in 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 step 1 to describe the battery degradation process; the specific process is:

相对于卡尔曼滤波算法和扩展卡尔曼滤波算法，粒子滤波的优势在于能够处理非线性和非高斯系统的参数辨识问题。需要注意的是，传统的粒子滤波算法直接将先验概率密度函数作为重要性函数。然而，从先验分布中产生的粒子不一定能够反映真实的粒子分布，并且忽略了系统当前时刻的测量信息。为此，本发明考虑采用无迹卡尔曼滤波算法为粒子滤波产生重要性函数。具体而言，即将无迹卡尔曼滤波算法的后验概率密度分布作为粒子滤波的先验概率密度分布。Compared with the Kalman filter algorithm and the extended Kalman filter algorithm, the advantage of the particle filter is that it can handle the parameter identification problem of nonlinear and non-Gaussian systems. It should be noted that the traditional particle filter algorithm directly uses the prior probability density function as the importance function. However, the particles generated from the prior distribution may not necessarily reflect the true particle distribution, and ignore the measurement information of the system at the current moment. For this reason, the present invention considers using the unscented Kalman filter algorithm to generate the importance function for the particle filter. Specifically, the posterior probability density distribution of the unscented Kalman filter algorithm is used as the prior probability density distribution of the particle filter.

基于无迹粒子滤波的电池退化模型可以采用状态空间的形式表示成公式(1)所示的形式：The battery degradation model based on unscented particle filtering can be expressed in the form of state space as shown in formula (1):

Figure BDA0003605634160000051
Figure BDA0003605634160000051

其中，xk为第k个工作循环下的状态变量，xk＝[ak,bk,ck,dk]T，ak为第k个工作循环下的状态变量的第一分量，bk为第k个工作循环下的状态变量的第二分量，ck为为第k个工作循环下的状态变量的第三分量，dk为第k个工作循环下的状态变量的第四分量，T为转置，k表示电池工作循环次数，xk-1为第k-1个工作循环下的状态变量，yk为电池容量，uk为过程噪声项，uk＝[ua,ub,uc,ud]T，ua为第k个工作循环下状态变量的第一分量的噪声项，ub为第k个工作循环下状态变量的第二分量的噪声项，uc为第k个工作循环下状态变量的第三分量的噪声项，ud为第k个工作循环下状态变量的第四分量的噪声项，vk为测量噪声项，vk∈R1×1，1×1为长宽为1的矩阵，f()代表状态转移方程的映射关系，h()代表测量方程的映射关系；Wherein, xk is the state variable under the kth working cycle, xk = [ ak , bk , ck , dk ] T , ak is the first component of the state variable under the kth working cycle, bk is the second component of the state variable under the kth working cycle, ck is the third component of the state variable under the kth working cycle, dk is the fourth component of the state variable under the kth working cycle, T is the transpose, k represents the number of battery working cycles, xk -1 is the state variable under the k-1th working cycle, yk is the battery capacity, uk is the process noise term, uk = [ ua , ub , uc , ud ] T , ua is the noise term of the first component of the state variable under the kth working cycle, ub is the noise term of the second component of the state variable under the kth working cycle, uc is the noise term of the third component of the state variable under the kth working cycle, ud is the noise term of the fourth component of the state variable under the kth working cycle, vk is the measurement noise term, vk∈R 1×1 , 1×1 is a matrix with a length and width of 1, f() represents the mapping relationship of the state transfer equation, and h() represents the mapping relationship of the measurement equation;

公式(1)中的状态转移方程xk＝f(xk-1)+uk可以写成公式(2)中的形式：The state transfer equation x k =f(x k-1 )+ uk in formula (1) can be written in the form of formula (2):

Figure BDA0003605634160000061
Figure BDA0003605634160000061

其中，

Figure BDA0003605634160000062
为过程变量噪声项ua的方差，
Figure BDA0003605634160000063
为过程变量噪声项ub的方差，
Figure BDA0003605634160000064
为过程变量噪声项uc的方差，
Figure BDA0003605634160000065
为过程变量噪声项ud的方差；ak-1为第k-1个工作循环下的状态变量的第一分量，bk-1为第k-1个工作循环下的状态变量的第二分量，ck-1为第k-1个工作循环下的状态变量的第三分量，dk-1为第k-1个工作循环下的状态变量的第四分量；
Figure BDA0003605634160000066
表示ua服从均值为0，方差为
Figure BDA0003605634160000067
的正态分布，
Figure BDA0003605634160000068
为均值为0，方差为
Figure BDA0003605634160000069
的正态分布，
Figure BDA00036056341600000610
为均值为0，方差为
Figure BDA00036056341600000611
的正态分布，
Figure BDA00036056341600000612
为为均值为0，方差为
Figure BDA00036056341600000613
的正态分布，
Figure BDA00036056341600000614
为均值为0，方差为
Figure BDA00036056341600000615
的正态分布；in,
Figure BDA0003605634160000062
is the variance of the process variable noise term u a ,
Figure BDA0003605634160000063
is the variance of the process variable noise term u b ,
Figure BDA0003605634160000064
is the variance of the process variable noise term u c ,
Figure BDA0003605634160000065
is the variance of the process variable noise term ud ; a k-1 is the first component of the state variable under the k-1th working cycle, b k-1 is the second component of the state variable under the k-1th working cycle, c k-1 is the third component of the state variable under the k-1th working cycle, and d k-1 is the fourth component of the state variable under the k-1th working cycle;
Figure BDA0003605634160000066
It means that u a has a mean of 0 and a variance of
Figure BDA0003605634160000067
The normal distribution of
Figure BDA0003605634160000068
The mean is 0 and the variance is
Figure BDA0003605634160000069
The normal distribution of
Figure BDA00036056341600000610
The mean is 0 and the variance is
Figure BDA00036056341600000611
The normal distribution of
Figure BDA00036056341600000612
The mean is 0 and the variance is
Figure BDA00036056341600000613
The normal distribution of
Figure BDA00036056341600000614
The mean is 0 and the variance is
Figure BDA00036056341600000615
Normal distribution of

相应地，测量方程yk＝h(xk)+vk可以写成公式(3)中的形式：Accordingly, the measurement equation y k =h(x k )+v k can be written in the form of formula (3):

Figure BDA00036056341600000616
Figure BDA00036056341600000616

其中，

Figure BDA00036056341600000617
为测量噪声的方差；
Figure BDA00036056341600000618
为均值为0，方差为
Figure BDA00036056341600000619
的正态分布；in,
Figure BDA00036056341600000617
is the variance of the measurement noise;
Figure BDA00036056341600000618
The mean is 0 and the variance is
Figure BDA00036056341600000619
Normal distribution of

整个无迹粒子滤波算法的实施步骤可以总结如下：The implementation steps of the entire unscented particle filter algorithm can be summarized as follows:

(1)粒子初始化，过程为：(1) Particle initialization, the process is:

根据初始状态变量x0，在初始状态变量x0的正态分布p(x0)中通过蒙特卡洛方法，也就是随机抽样的方式可以产生粒子集合

Figure BDA00036056341600000620
Ns为粒子的总数，且每个粒子的初始化权重值为
Figure BDA00036056341600000621
According to the initial state variable x 0 , a particle set can be generated in the normal distribution p(x 0 ) of the initial state variable x 0 by the Monte Carlo method, that is, random sampling.
Figure BDA00036056341600000620
N s is the total number of particles, and the initialization weight value of each particle is
Figure BDA00036056341600000621

初始化粒子的均值

Figure BDA00036056341600000622
以及初始化粒子的方差S0可以写成公式(4)—(5)所示的形式：Initialize the mean of the particles
Figure BDA00036056341600000622
And the variance S 0 of the initialized particle can be written as shown in formulas (4)-(5):

Figure BDA0003605634160000071
Figure BDA0003605634160000071

Figure BDA0003605634160000072
Figure BDA0003605634160000072

其中，i为粒子的编号，

Figure BDA0003605634160000073
为初始化粒子的均值，S0为初始化粒子的方差，E()表示期望运算符，T表示矩阵转置运算符，r0为初始化粒子；Where i is the number of the particle,
Figure BDA0003605634160000073
is the mean of the initialized particle, S 0 is the variance of the initialized particle, E() represents the expectation operator, T represents the matrix transposition operator, and r 0 is the initialized particle;

(2)通过UKF算法计算粒子所表征的状态变量rk的后验期望

Figure BDA0003605634160000074
和协方差
Figure BDA0003605634160000075
过程为：(2) Calculate the posterior expectation of the state variable r k represented by the particle through the UKF algorithm
Figure BDA0003605634160000074
and covariance
Figure BDA0003605634160000075
The process is:

在完成k-1次迭代过程后，构建各粒子的Sigma采样点集

Figure BDA0003605634160000076
如公式(6)所示：After completing k-1 iterations, construct the Sigma sampling point set for each particle
Figure BDA0003605634160000076
As shown in formula (6):

Figure BDA0003605634160000077
Figure BDA0003605634160000077

其中，

Figure BDA0003605634160000078
为第i个粒子在k-1次工作过程中的均值向量，
Figure BDA0003605634160000079
为第i个粒子在k-1次工作过程中的协方差矩阵，参数α为UKF中的系数值；in,
Figure BDA0003605634160000078
is the mean vector of the i-th particle in the k-1 working process,
Figure BDA0003605634160000079
is the covariance matrix of the ith particle in the k-1 working process, and the parameter α is the coefficient value in UKF;

各粒子的Sigma点集的样本点个数为2n+1，n为状态变量的数目，各粒子的Sigma点集的样本点的权重表达式分别如公式(7)—(9)所示：The number of sample points of the Sigma point set of each particle is 2n+1, where n is the number of state variables. The weight expressions of the sample points of the Sigma point set of each particle are shown in formulas (7) to (9):

Figure BDA00036056341600000710
Figure BDA00036056341600000710

Figure BDA00036056341600000711
Figure BDA00036056341600000711

Figure BDA00036056341600000712
Figure BDA00036056341600000712

其中，

Figure BDA00036056341600000713
为第i个粒子第1个Sigma点集的样本点所代表的过程变量权重，
Figure BDA00036056341600000714
为第i个粒子第1个Sigma点集的样本点所代表的测量变量权重，
Figure BDA00036056341600000715
为第i个粒子第j个Sigma点集的样本点所代表的过程变量权重，
Figure BDA00036056341600000716
为第i个粒子第j个Sigma点集的样本点所代表的测量变量权重，参数β,λ为UKF中的系数值，n为状态变量的数目；in,
Figure BDA00036056341600000713
is the weight of the process variable represented by the sample point of the first Sigma point set of the i-th particle,
Figure BDA00036056341600000714
is the weight of the measured variable represented by the sample point of the first Sigma point set of the i-th particle,
Figure BDA00036056341600000715
is the weight of the process variable represented by the sample point of the jth Sigma point set of the i-th particle,
Figure BDA00036056341600000716
is the weight of the measured variable represented by the sample point of the jth Sigma point set of the i-th particle, the parameters β and λ are the coefficient values in UKF, and n is the number of state variables;

在此基础上，UKF的状态估计表达式可以写成公式(10)—(11)所示的形式：On this basis, the state estimation expression of UKF can be written as shown in formulas (10)-(11):

Figure BDA00036056341600000717
Figure BDA00036056341600000717

Figure BDA0003605634160000081
Figure BDA0003605634160000081

其中，

Figure BDA0003605634160000082
为第k次工作过程第i个粒子中第j个Sigma点所代表的先验状态向量，f()为状态转移方程的映射关系，
Figure BDA0003605634160000083
为第k-1次工作过程第i个粒子中第j个Sigma点所代表的后验状态向量，
Figure BDA0003605634160000084
为第k-1次工作过程第i个粒子所代表的噪声变量，
Figure BDA0003605634160000085
为第k次工作过程第i个粒子的加权先验状态向量；in,
Figure BDA0003605634160000082
is the prior state vector represented by the jth Sigma point in the i-th particle in the k-th working process, f() is the mapping relationship of the state transfer equation,
Figure BDA0003605634160000083
is the posterior state vector represented by the jth Sigma point in the i-th particle in the k-1th working process,
Figure BDA0003605634160000084
is the noise variable represented by the i-th particle in the k-1th working process,
Figure BDA0003605634160000085
is the weighted prior state vector of the ith particle in the kth working process;

协方差估计的表达式则可以写成公式(12)所示的形式：The expression of covariance estimation can be written as shown in formula (12):

Figure BDA0003605634160000086
Figure BDA0003605634160000086

其中，

Figure BDA0003605634160000087
为第k次工作过程第i个粒子中状态变量的先验协方差，
Figure BDA0003605634160000088
为第k次工作过程第i个粒子中噪声变量的协方差矩阵，
Figure BDA0003605634160000089
diag表示对角矩阵，
Figure BDA00036056341600000810
为第k次工作过程过程变量噪声项ua的方差，
Figure BDA00036056341600000811
为第k次工作过程过程变量噪声项ub的方差，
Figure BDA00036056341600000812
为第k次工作过程的过程变量噪声项uc的方差，
Figure BDA00036056341600000813
为第k次工作过程过程变量噪声项ud的方差；in,
Figure BDA0003605634160000087
is the prior covariance of the state variables in the ith particle of the kth working process,
Figure BDA0003605634160000088
is the covariance matrix of the noise variable in the i-th particle of the k-th working process,
Figure BDA0003605634160000089
diag represents a diagonal matrix,
Figure BDA00036056341600000810
is the variance of the process variable noise term u a of the kth working process,
Figure BDA00036056341600000811
is the variance of the process variable noise term u b of the kth working process,
Figure BDA00036056341600000812
is the variance of the process variable noise term u c of the kth working process,
Figure BDA00036056341600000813
is the variance of the process variable noise term ud in the kth working process;

相应的，对于前向的一步预测，有公式(13)—(14)所示的形式：Correspondingly, for the one-step forward prediction, we have the form shown in formulas (13)-(14):

Figure BDA00036056341600000814
Figure BDA00036056341600000814

Figure BDA00036056341600000815
Figure BDA00036056341600000815

其中，

Figure BDA00036056341600000816
为第i个粒子中第j个Sigma点所代表的状态向量，
Figure BDA00036056341600000817
为第i个粒子中第j个Sigma点所代表的UKF前向一步预测值，
Figure BDA00036056341600000818
为第i个粒子所代表的UKF预测值，
Figure BDA00036056341600000819
为第k次工作过程第i个粒子的测量噪声项，h()代表测量方程的映射关系；in,
Figure BDA00036056341600000816
is the state vector represented by the jth Sigma point in the ith particle,
Figure BDA00036056341600000817
is the UKF one-step-forward prediction value represented by the j-th Sigma point in the i-th particle,
Figure BDA00036056341600000818
is the UKF prediction value represented by the i-th particle,
Figure BDA00036056341600000819
is the measurement noise term of the ith particle in the kth working process, and h() represents the mapping relationship of the measurement equation;

在此基础上，可以计算自协方差矩阵

Figure BDA00036056341600000820
和互协方差矩阵
Figure BDA00036056341600000821
如公式(15)—(16)所示：On this basis, the autocovariance matrix can be calculated
Figure BDA00036056341600000820
and the cross-covariance matrix
Figure BDA00036056341600000821
As shown in formulas (15)-(16):

Figure BDA00036056341600000822
Figure BDA00036056341600000822

Figure BDA0003605634160000091
Figure BDA0003605634160000091

其中，

Figure BDA0003605634160000092
为第k次工作过程第i个粒子的噪声方差；in,
Figure BDA0003605634160000092
is the noise variance of the i-th particle in the k-th working process;

由此可以获得卡尔曼滤波增益向量Kk如公式(17)所示：Thus, the Kalman filter gain vector K k can be obtained as shown in formula (17):

Figure BDA0003605634160000093
Figure BDA0003605634160000093

其中，

Figure BDA0003605634160000094
为第i个粒子卡尔曼滤波增益向量，
Figure BDA0003605634160000095
为自协方差矩阵，
Figure BDA0003605634160000096
为互协方差矩阵，随着第k次迭代过程中新的电池容量数据的加入，UKF的状态更新和协方差更新可以表达成公式(18)—(19)所示的形式：in,
Figure BDA0003605634160000094
is the Kalman filter gain vector of the i-th particle,
Figure BDA0003605634160000095
is the autocovariance matrix,
Figure BDA0003605634160000096
is the mutual covariance matrix. With the addition of new battery capacity data in the kth iteration, the state update and covariance update of UKF can be expressed as shown in formulas (18)-(19):

Figure BDA0003605634160000097
Figure BDA0003605634160000097

Figure BDA0003605634160000098
Figure BDA0003605634160000098

其中，

Figure BDA0003605634160000099
为通过UKF算法计算粒子的状态变量rk的后验期望(UKF算法状态变量的更新值)，
Figure BDA00036056341600000910
为通过UKF算法计算粒子的状态变量rk的协方差(UKF算法协方差矩阵的更新值)；in,
Figure BDA0003605634160000099
is the posterior expectation of the particle state variable r k calculated by the UKF algorithm (the updated value of the UKF algorithm state variable),
Figure BDA00036056341600000910
is the covariance of the particle state variable r k calculated by the UKF algorithm (the updated value of the UKF algorithm covariance matrix);

根据公式(18)—(19)中的结果，对

Figure BDA00036056341600000911
求取均值获得
Figure BDA00036056341600000912
对
Figure BDA00036056341600000913
求取均值获得
Figure BDA00036056341600000914
According to the results in formulas (18)-(19),
Figure BDA00036056341600000911
Get the mean
Figure BDA00036056341600000912
right
Figure BDA00036056341600000913
Get the mean
Figure BDA00036056341600000914

基于

Figure BDA00036056341600000915
和
Figure BDA00036056341600000916
获得在UKF算法下粒子的后验概率密度分布
Figure BDA00036056341600000917
based on
Figure BDA00036056341600000915
and
Figure BDA00036056341600000916
Obtain the posterior probability density distribution of particles under the UKF algorithm
Figure BDA00036056341600000917

(3)粒子权重更新，过程为：(3) Particle weight update, the process is:

根据(2)获得粒子的重要性密度函数可以表示为：

Figure BDA00036056341600000918
According to (2), the importance density function of the particle can be expressed as:
Figure BDA00036056341600000918

其中，

Figure BDA00036056341600000919
表示第k次工作过程第i个粒子的重要性密度函数，
Figure BDA00036056341600000920
为第i个粒子初始状态到第k-1次工作过程所代表的状态变量，y1:k表示第1次工作过程到第k次工作过程的电池容量数据，
Figure BDA00036056341600000921
为第i个粒子初始状态到第k-1次工作过程所代表的状态变量和第1次工作过程到第k次工作过程的电池容量数据先验条件下
Figure BDA00036056341600000922
的概率密度函数，
Figure BDA00036056341600000923
为均值为
Figure BDA00036056341600000924
方差为
Figure BDA00036056341600000925
的正态分布，～为服从某分布的含义；in,
Figure BDA00036056341600000919
represents the importance density function of the i-th particle in the k-th working process,
Figure BDA00036056341600000920
is the state variable represented by the initial state of the i-th particle to the k-1th working process, y 1:k represents the battery capacity data from the 1st working process to the kth working process,
Figure BDA00036056341600000921
The state variables represented by the initial state of the i-th particle to the k-1th working process and the battery capacity data from the 1st working process to the kth working process under the prior conditions
Figure BDA00036056341600000922
The probability density function of
Figure BDA00036056341600000923
The mean is
Figure BDA00036056341600000924
The variance is
Figure BDA00036056341600000925
Normal distribution, ~ means obeying a certain distribution;

在此基础上，粒子的权重更新表达式可以写成公式(20)中的形式：On this basis, the particle weight update expression can be written as in formula (20):

Figure BDA00036056341600000926
Figure BDA00036056341600000926

其中，

Figure BDA0003605634160000101
为第i个粒子在k-1个工作循环下的权重，
Figure BDA0003605634160000102
为第i个粒子在k个工作循环下的权重，
Figure BDA0003605634160000103
为yk在
Figure BDA0003605634160000104
条件下的概率密度函数，
Figure BDA0003605634160000105
为
Figure BDA0003605634160000106
在
Figure BDA0003605634160000107
条件下的概率密度函数；in,
Figure BDA0003605634160000101
is the weight of the ith particle under k-1 working cycles,
Figure BDA0003605634160000102
is the weight of the ith particle under k working cycles,
Figure BDA0003605634160000103
For y k
Figure BDA0003605634160000104
The probability density function under the condition,
Figure BDA0003605634160000105
for
Figure BDA0003605634160000106
exist
Figure BDA0003605634160000107
Probability density function under the conditions;

Figure BDA0003605634160000108
在经过粒子的权重归一化，可以获得归一化的粒子权重
Figure BDA0003605634160000109
表达式如下所示：
Figure BDA0003605634160000108
After the particle weight is normalized, the normalized particle weight can be obtained
Figure BDA0003605634160000109
The expression is as follows:

Figure BDA00036056341600001010
Figure BDA00036056341600001010

(4)粒子重采样，过程为：(4) Particle resampling, the process is:

为了改善粒子滤波算法的状态估计性能，本发明采用重采样的方式保留粒子集合中权重较大的粒子，减少权重较小的粒子。假设有效粒子数目为Nef，可以定义成公式(22)中的形式：In order to improve the state estimation performance of the particle filter algorithm, the present invention adopts a resampling method to retain particles with larger weights in the particle set and reduce particles with smaller weights. Assuming that the number of effective particles is Nef , it can be defined as in formula (22):

Figure BDA00036056341600001011
Figure BDA00036056341600001011

设定阀值Nth，如果有效粒子数目Nef小于Nth，则需要进行粒子的重采样；Set the threshold N th . If the number of valid particles Nef is less than N th , resampling of particles is required.

由于随机重采样实现过程简单且具备较高的采样效率，故本发明选用随机重采样的方式实现粒子滤波的重采样过程；Since random resampling is simple to implement and has high sampling efficiency, the present invention uses random resampling to implement the resampling process of particle filtering;

如果有效粒子数目Nef大于等于Nth，则不需要进行粒子的重采样；If the number of effective particles Nef is greater than or equal to Nth , there is no need to resample the particles;

在完成(4)后，每一个粒子的权重为

Figure BDA00036056341600001012
(需要重采样和不需要重采样的每一个粒子的权重皆为
Figure BDA00036056341600001013
)；After completing (4), the weight of each particle is
Figure BDA00036056341600001012
(The weights of each particle that needs to be resampled and does not need to be resampled are
Figure BDA00036056341600001013
);

(5)状态更新，过程为：(5) Status update: The process is as follows:

在重采样获得的粒子权重的基础上，即可通过公式(23)—(24)获得无迹粒子滤波算法(UPF)后验条件下的状态变量期望

Figure BDA00036056341600001014
和协方差
Figure BDA00036056341600001015
Based on the particle weights obtained by resampling, the expected state variables under the posterior condition of the unscented particle filter (UPF) can be obtained by formulas (23)-(24):
Figure BDA00036056341600001014
and covariance
Figure BDA00036056341600001015

Figure BDA00036056341600001016
Figure BDA00036056341600001016

Figure BDA00036056341600001017
Figure BDA00036056341600001017

其中，

Figure BDA00036056341600001018
为UPF算法后验条件下的状态变量期望，
Figure BDA00036056341600001019
为UPF算法后验条件下的协方差。in,
Figure BDA00036056341600001018
is the expectation of the state variables under the posterior condition of the UPF algorithm,
Figure BDA00036056341600001019
is the covariance under the posterior condition of the UPF algorithm.

其它步骤及参数与具体实施方式一相同。The other steps and parameters are the same as those in the first embodiment.

具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述步骤3中采用基于最大期望算法自适应地估计电池退化模型中的过程噪声和测量噪声；Specific implementation method three: This implementation method is different from specific implementation method one or two in that the process noise and measurement noise in the battery degradation model are adaptively estimated based on the maximum expectation algorithm in step 3;

具体而言，针对基于无迹粒子滤波的动态退化模型中的过程噪声和测量噪声，采用最大期望算法进行自适应估计，利用电池从初始状态到当前时刻的所有容量数据对过程噪声和测量噪声进行递归更新；Specifically, for the process noise and measurement noise in the dynamic degradation model based on unscented particle filtering, the maximum expectation algorithm is used for adaptive estimation, and all the capacity data of the battery from the initial state to the current moment are used to recursively update the process noise and measurement noise;

具体过程为：The specific process is:

过程噪声和测量噪声所构成的未知参数向量可以表示为：

Figure BDA0003605634160000111
The unknown parameter vector composed of process noise and measurement noise can be expressed as:
Figure BDA0003605634160000111

EM算法的本质在于通过在后验概率的前提下最大化联合似然函数p(x0:k,y1:k|Ξ)来实现未知参数Ξ的极大似然估计。在已知自初始到第k个工作循环的电池容量数据y1:k＝[y1,y2,…,yk]的前提下，可以构建联合对数似然函数Lk(Ξ)如公式(25)所示：The essence of the EM algorithm is to achieve the maximum likelihood estimation of the unknown parameter Ξ by maximizing the joint likelihood function p(x 0:k ,y 1:k |Ξ) under the premise of the posterior probability. Under the premise that the battery capacity data y 1:k = [y 1 ,y 2,…,y k ] from the initial to the kth working cycle is known, the joint log-likelihood function L k (Ξ) can be constructed as shown in formula (25):

Figure BDA0003605634160000112
Figure BDA0003605634160000112

其中，

Figure BDA0003605634160000113
为第t次工作过程第i个粒子的状态变量，
Figure BDA0003605634160000114
Figure BDA0003605634160000115
为第i个粒子的第t次工作过程中的过程变量值，Ns为粒子的总数，且每个粒子的初始化权重值为
Figure BDA0003605634160000116
i为粒子的编号，
Figure BDA0003605634160000117
为第1次工作过程到第k次工作过程第i个粒子所表征的状态变量，
Figure BDA0003605634160000118
为未知参数条件下的第1次工作过程到第k次工作过程第i个粒子所表征的状态变量与第1次工作过程到第k次工作过程的电池容量数据所构成的联合概率密度函数，
Figure BDA0003605634160000119
为未知参数条件下第1次工作过程到第k次工作过程第i个粒子所表征的状态变量的条件概率密度函数，
Figure BDA00036056341600001110
为未知参数条件和第1次工作过程到第k次工作过程第i个粒子所表征的状态变量条件下第1次工作过程到第k次工作过程的电池容量数据的条件概率密度函数，
Figure BDA0003605634160000121
为未知参数条件和第t-1次工作过程第i个粒子状态变量条件下第t次工作过程第i个粒子状态变量的条件概率密度函数，
Figure BDA0003605634160000122
为未知参数条件和第t次工作过程第i个粒子状态变量条件下第t次工作过程电池容量数据的条件概率密度函数，
Figure BDA0003605634160000123
第t-1次工作过程第i个粒子状态变量，
Figure BDA0003605634160000124
为未知参数条件下的第1次工作过程到第k次工作过程第i个粒子所表征的状态变量与第1次工作过程到第k次工作过程的电池容量数据的关系，
Figure BDA0003605634160000125
为未知参数条件下第1次工作过程到第k次工作过程第i个粒子所表征的状态变量的关系，
Figure BDA0003605634160000126
为未知参数条件和第1次工作过程到第k次工作过程第i个粒子所表征的状态变量条件下第1次工作过程到第k次工作过程的电池容量数据的关系，
Figure BDA0003605634160000127
未知参数条件和第t-1次工作过程第i个粒子状态变量条件下第t次工作过程第i个粒子状态变量的关系，
Figure BDA0003605634160000128
为未知参数条件和第t次工作过程第i个粒子状态变量条件下第t次工作过程电池容量数据的关系；in,
Figure BDA0003605634160000113
is the state variable of the ith particle in the tth working process,
Figure BDA0003605634160000114
Figure BDA0003605634160000115
is the process variable value of the t-th working process of the i-th particle, N s is the total number of particles, and the initialization weight value of each particle is
Figure BDA0003605634160000116
i is the number of the particle,
Figure BDA0003605634160000117
is the state variable represented by the i-th particle from the 1st working process to the k-th working process,
Figure BDA0003605634160000118
is the joint probability density function composed of the state variables represented by the i-th particle from the 1st working process to the k-th working process under unknown parameter conditions and the battery capacity data from the 1st working process to the k-th working process,
Figure BDA0003605634160000119
is the conditional probability density function of the state variable represented by the i-th particle from the 1st working process to the k-th working process under unknown parameter conditions,
Figure BDA00036056341600001110
is the conditional probability density function of the battery capacity data from the 1st working process to the kth working process under the condition of unknown parameters and the state variables represented by the i-th particle from the 1st working process to the kth working process,
Figure BDA0003605634160000121
is the conditional probability density function of the state variable of the ith particle in the t-th working process under the condition of unknown parameters and the state variable of the ith particle in the t-1th working process,
Figure BDA0003605634160000122
is the conditional probability density function of the battery capacity data of the t-th working process under the condition of unknown parameters and the i-th particle state variable of the t-th working process,
Figure BDA0003605634160000123
The state variable of the i-th particle in the t-1th working process,
Figure BDA0003605634160000124
is the relationship between the state variable represented by the i-th particle from the 1st working process to the k-th working process under unknown parameter conditions and the battery capacity data from the 1st working process to the k-th working process,
Figure BDA0003605634160000125
is the relationship between the state variables represented by the i-th particle from the 1st working process to the k-th working process under unknown parameter conditions,
Figure BDA0003605634160000126
is the relationship between the battery capacity data from the 1st working process to the kth working process under the unknown parameter conditions and the state variable conditions represented by the i-th particle from the 1st working process to the kth working process,
Figure BDA0003605634160000127
The relationship between the unknown parameter conditions and the state variables of the i-th particle in the t-th working process under the condition of the i-th particle state variables in the t-1th working process,
Figure BDA0003605634160000128
is the relationship 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;

基于步骤2中的UPF的自适应状态更新模型公式(1)—(3)中状态转移方程和测量方程中的内容，可以获得公式(26)—(27)中的关系表达式：Based on the contents of the state transfer equation and measurement equation in the adaptive state update model formula (1)-(3) of the UPF in step 2, the relational expressions in formulas (26)-(27) can be obtained:

Figure BDA0003605634160000129
Figure BDA0003605634160000129

Figure BDA00036056341600001210
Figure BDA00036056341600001210

其中，

Figure BDA00036056341600001211
为均值为
Figure BDA00036056341600001212
方差为
Figure BDA00036056341600001213
的正态分布，t为工作过程的编号，
Figure BDA00036056341600001214
为第t次工作过程第i个粒子的噪声变量方差，
Figure BDA00036056341600001215
第t次工作过程第i个粒子的过程变量协方差矩阵，
Figure BDA00036056341600001216
diag表示对角矩阵，
Figure BDA00036056341600001217
为分别第t次工作过程过程变量噪声项ua，ub，uc，ud的方差；in,
Figure BDA00036056341600001211
The mean is
Figure BDA00036056341600001212
The variance is
Figure BDA00036056341600001213
The normal distribution of , t is the number of the work process,
Figure BDA00036056341600001214
is the noise variable variance of the ith particle in the tth working process,
Figure BDA00036056341600001215
The process variable covariance matrix of the ith particle in the tth working process is:
Figure BDA00036056341600001216
diag represents a diagonal matrix,
Figure BDA00036056341600001217
are the variances of the process variable noise terms ua , ub , uc , ud of the t-th working process respectively;

将公式(26)—(27)代入到公式(25)中，且忽略常数项，则联合对数似然函数可以进一步表示为公式(28)所示的形式：Substituting formulas (26)-(27) into formula (25) and ignoring the constant term, the joint log-likelihood function can be further expressed as shown in formula (28):

Figure BDA0003605634160000131
Figure BDA0003605634160000131

其中，∝表示正比于的含义；Among them, ∝ means proportional to;

在已知自初始到第k个工作循环的电池容量数据y1:k＝[y1,y2,…,yk]的前提下，对于在EM算法的第m次迭代过程，自适应退化模型的未知参数估计向量可以表示为：Under the premise that the battery capacity data y 1:k = [y 1 ,y 2 ,…,y k ] from the initial to the kth working cycle is known, for the mth iteration process of the EM algorithm, the unknown parameter estimation vector of the adaptive degradation model can be expressed as:

Figure BDA0003605634160000132
Figure BDA0003605634160000132

其中，

Figure BDA0003605634160000133
为第m次迭代过程中第k次工作过程噪声变量ua的方差，
Figure BDA0003605634160000134
为第m次迭代过程中第k次工作过程噪声变量ub的方差，
Figure BDA0003605634160000135
为第m次迭代过程中第k次工作过程噪声变量uc的方差，
Figure BDA0003605634160000136
为第m次迭代过程中第k次工作过程噪声变量ud的方差，
Figure BDA0003605634160000137
为第m次迭代过程中第k次工作测量噪声变量的方差；in,
Figure BDA0003605634160000133
is the variance of the noise variable u a in the kth working process during the mth iteration,
Figure BDA0003605634160000134
is the variance of the noise variable u b in the kth working process during the mth iteration,
Figure BDA0003605634160000135
is the variance of the noise variable u c in the kth working process during the mth iteration,
Figure BDA0003605634160000136
is the variance of the noise variable ud in the kth working process during the mth iteration,
Figure BDA0003605634160000137
The variance of the noise variable measured for the kth job during the mth iteration;

根据公式(25)—(28)的推导，以及EM算法的基本原理；According to the derivation of formulas (25)-(28) and the basic principle of EM algorithm;

对于第m+1次迭代过程，可以分为E步和M步，可以表示成公式(29)—(30)所示：For the m+1th iteration process, it can be divided into the E step and the M step, which can be expressed as shown in formulas (29)-(30):

(1)E步骤：计算(1) Step E: Calculation

Figure BDA0003605634160000138
Figure BDA0003605634160000138

(2)M步骤：计算(2) Step M: Calculation

Figure BDA0003605634160000139
Figure BDA0003605634160000139

其中，

Figure BDA00036056341600001310
为Ξ在
Figure BDA00036056341600001311
条件下的联合对数似然函数，
Figure BDA00036056341600001312
为r1:k,y1:k在
Figure BDA00036056341600001313
条件下的期望运算；
Figure BDA00036056341600001314
为第m+1次迭代过程对第k个工作过程的未知参数估计向量，通过迭代E步骤和M步骤直到满足某一个收敛判据为止，从而获得退化模型的未知参数值；in,
Figure BDA00036056341600001310
For
Figure BDA00036056341600001311
The joint log-likelihood function under the condition,
Figure BDA00036056341600001312
For r 1:k ,y 1:k
Figure BDA00036056341600001313
Expected operation under conditions;
Figure BDA00036056341600001314
The unknown parameter estimation vector of the kth working process in the m+1th iteration process is obtained by iterating the E step and the M step until a certain convergence criterion is met, thereby obtaining the unknown parameter value of the degradation model;

根据公式(28)，公式(29)可以扩展地写成公式(31)所示的形式：According to formula (28), formula (29) can be expanded to the form shown in formula (31):

Figure BDA0003605634160000141
Figure BDA0003605634160000141

在公式(31)中，

Figure BDA0003605634160000142
Figure BDA0003605634160000143
和
Figure BDA0003605634160000144
是第i个粒子基于自初始到第k个工作循环的电池容量数据y1:k＝[y1,y2,…,yk]的隐含状态变量的条件期望；In formula (31),
Figure BDA0003605634160000142
Figure BDA0003605634160000143
and
Figure BDA0003605634160000144
is the conditional expectation of the implicit state variables of the ith particle based on the battery capacity data y 1:k = [y 1 , y 2 , …, y k ] from the initial to the kth working cycle;

为了计算在后验条件下的状态变量期望，本发明采用Rauch-Tung-Striebel(RTS)最优平滑算法进行后向平滑实现。基于UPF算法的前向迭代，可以获得粒子所代表的状态变量的期望和协方差的估计如公式(32)和公式(33)所示：In order to calculate the expectation of the state variables under the posterior condition, the present invention adopts the Rauch-Tung-Striebel (RTS) optimal smoothing algorithm for backward smoothing. Based on the forward iteration of the UPF algorithm, the expectation and covariance estimation of the state variables represented by the particles can be obtained as shown in formula (32) and formula (33):

Figure BDA0003605634160000145
Figure BDA0003605634160000145

Figure BDA0003605634160000146
Figure BDA0003605634160000146

其中，

Figure BDA0003605634160000151
为第i个粒子的状态变量在RTS算法中的期望，
Figure BDA0003605634160000152
为第i个粒子的状态变量在RTS算法中的协方差，
Figure BDA0003605634160000153
为无迹粒子滤波算法第i个粒子于第k次工作过程的后验状态变量期望，
Figure BDA0003605634160000154
为无迹粒子滤波算法第i个粒子于第k次工作过程的后验协方差矩阵；in,
Figure BDA0003605634160000151
is the expectation of the state variable of the ith particle in the RTS algorithm,
Figure BDA0003605634160000152
is the covariance of the state variable of the ith particle in the RTS algorithm,
Figure BDA0003605634160000153
is the posterior state variable expectation of the i-th particle in the k-th working process of the unscented particle filter algorithm,
Figure BDA0003605634160000154
is the posterior covariance matrix of the i-th particle in the k-th working process of the unscented particle filter algorithm;

相应地，k-1和k工作循环之间的状态变量之间的协方差可以表示成公式(34)所示的形式：Accordingly, the covariance between the state variables between k-1 and k working cycles can be expressed as shown in formula (34):

Figure BDA0003605634160000155
Figure BDA0003605634160000155

其中，

Figure BDA0003605634160000156
为第i个粒子所代表的k-1和k工作循环之间的状态变量之间的协方差，
Figure BDA0003605634160000157
为第i个粒子的UKF增益，
Figure BDA0003605634160000158
为k-1和k工作循环之间的状态变量之间的协方差，I4×4为四维单位阵，
Figure BDA0003605634160000159
为无迹粒子滤波算法中第i个粒子所代表的k-1次工作过程的状态变量之间的后验协方差矩阵；in,
Figure BDA0003605634160000156
is the covariance between the state variables between the k-1 and k working cycles represented by the ith particle,
Figure BDA0003605634160000157
is the UKF gain of the ith particle,
Figure BDA0003605634160000158
is the covariance between the state variables between k-1 and k working cycles, I 4×4 is the four-dimensional unit matrix,
Figure BDA0003605634160000159
is the posterior covariance matrix between the state variables of the k-1 working processes represented by the i-th particle in the unscented particle filter algorithm;

根据公式(32)—(34)中的初始化条件，即可以进行RTS平滑运算；RTS平滑增益的表达式如公式(35)所示：According to the initialization conditions in formulas (32)-(34), RTS smoothing operation can be performed; the expression of RTS smoothing gain is shown in formula (35):

Figure BDA00036056341600001510
Figure BDA00036056341600001510

其中，

Figure BDA00036056341600001511
为第i个粒子在第t个工作过程的RTS平滑增益，
Figure BDA00036056341600001512
为第i个粒子在第t个工作过程的UKF后验协方差，
Figure BDA00036056341600001513
为第i个粒子在第t个工作过程的UKF前向一步预测协方差；in,
Figure BDA00036056341600001511
is the RTS smoothing gain of the ith particle in the tth working process,
Figure BDA00036056341600001512
is the UKF posterior covariance of the ith particle in the tth working process,
Figure BDA00036056341600001513
is the UKF forward one-step prediction covariance of the ith particle in the tth working process;

相应地，后向迭代过程中的状态变量

Figure BDA00036056341600001514
和协方差
Figure BDA00036056341600001515
的更新如公式(36)—(37)所示：Accordingly, the state variables in the backward iteration process are
Figure BDA00036056341600001514
and covariance
Figure BDA00036056341600001515
The update of is shown in formulas (36)-(37):

Figure BDA00036056341600001516
Figure BDA00036056341600001516

Figure BDA00036056341600001517
Figure BDA00036056341600001517

其中，

Figure BDA0003605634160000161
为第i个粒子在RTS后向平滑中对于第t+1个工作过程的后验状态向量，
Figure BDA0003605634160000162
为第i个粒子在UKF前向预测中对于第t+1个工作过程状态向量的预测值，
Figure BDA0003605634160000163
为第i个粒子在RTS后向平滑中对于第t个工作过程的先验状态向量，
Figure BDA0003605634160000164
为第i个粒子在RTS后向平滑中对于第t个工作过程的后验状态向量，
Figure BDA0003605634160000165
为第i个粒子在RTS后向平滑中对于第t+1个工作过程的协方差，
Figure BDA0003605634160000166
为第i个粒子为UKF前向预测中对于第t+1个工作过程协方差的预测值，
Figure BDA0003605634160000167
为第i个粒子在UKF前向运算中对于第t个工作过程的后验协方差，
Figure BDA0003605634160000168
为第i个粒子在RTS后向平滑中对于第t个工作过程的协方差；in,
Figure BDA0003605634160000161
is the posterior state vector of the ith particle in the RTS backward smoothing for the t+1th working process,
Figure BDA0003605634160000162
is the predicted value of the state vector of the t+1th working process of the ith particle in the UKF forward prediction,
Figure BDA0003605634160000163
is the prior state vector of the ith particle in the RTS backward smoothing for the tth working process,
Figure BDA0003605634160000164
is the posterior state vector of the ith particle in RTS backward smoothing for the tth working process,
Figure BDA0003605634160000165
is the covariance of the ith particle in the RTS backward smoothing for the t+1th working process,
Figure BDA0003605634160000166
is the predicted value of the covariance of the t+1th working process in the UKF forward prediction for the ith particle,
Figure BDA0003605634160000167
is the posterior covariance of the ith particle in the UKF forward operation for the tth working process,
Figure BDA0003605634160000168
is the covariance of the ith particle in the RTS backward smoothing for the tth working process;

除此之外，相邻两个循环之间状态矩阵之间的协方差可以表示成公式(38)所示的形式：In addition, the covariance between the state matrices between two adjacent cycles can be expressed as shown in formula (38):

Figure BDA0003605634160000169
Figure BDA0003605634160000169

其中，

Figure BDA00036056341600001610
为RTS后向平滑中第t-1个工作过程和第t个工作过程状态变量的协方差，
Figure BDA00036056341600001611
为RTS后向平滑中第t个工作过程和第t+1个工作过程状态变量的协方差；in,
Figure BDA00036056341600001610
is the covariance of the state variables of the t-1th working process and the tth working process in RTS backward smoothing,
Figure BDA00036056341600001611
is the covariance of the state variables of the t-th working process and the t+1-th working process in RTS backward smoothing;

基于公式(32)—(38)，可以求解条件期望表达式

Figure BDA00036056341600001612
Figure BDA00036056341600001613
和
Figure BDA00036056341600001614
如式(39)—(42)所示：Based on formulas (32)-(38), the conditional expectation expression can be solved
Figure BDA00036056341600001612
Figure BDA00036056341600001613
and
Figure BDA00036056341600001614
As shown in formulas (39)-(42):

Figure BDA00036056341600001615
Figure BDA00036056341600001615

Figure BDA00036056341600001616
Figure BDA00036056341600001616

Figure BDA00036056341600001617
Figure BDA00036056341600001617

Figure BDA00036056341600001618
Figure BDA00036056341600001618

其中，

Figure BDA0003605634160000171
为RTS后向平滑中第t-1个工作过程第i个粒子所代表的状态变量；in,
Figure BDA0003605634160000171
is the state variable represented by the i-th particle in the t-1-th working process in RTS backward smoothing;

将公式(39)—(42)中的条件期望代入到公式(31)中，可以写成公式(43)所示的形式Substituting the conditional expectations in formulas (39)-(42) into formula (31), it can be written as shown in formula (43):

Figure BDA0003605634160000172
Figure BDA0003605634160000172

基于公式(31)—(43)中E步骤的计算结果，根据公式(30)计算EM算法的M步骤，针对退化模型的未知参数估计向量

Figure BDA0003605634160000173
中每一个参数的偏导数，令其值为0，可以求解出公式(44)—(45)的结果。Based on the calculation results of step E in formulas (31)-(43), the M step of the EM algorithm is calculated according to formula (30), and the unknown parameter estimation vector of the degradation model is
Figure BDA0003605634160000173
By taking the partial derivative of each parameter in and setting its value to 0, we can solve the results of formulas (44)-(45).

Figure BDA0003605634160000174
Figure BDA0003605634160000174

Figure BDA0003605634160000175
Figure BDA0003605634160000175

通过不断进行E步骤和M步骤的迭代，直到满足

Figure BDA0003605634160000176
或者达到最大迭代次数10次为止停止迭代，由此实现关于退化模型中未知参数的自适应估计；By continuously iterating the E and M steps until the
Figure BDA0003605634160000176
Or the iteration is stopped when the maximum number of iterations reaches 10, thereby realizing adaptive estimation of unknown parameters in the degradation model;

其中，ε为人为设定的参数。Among them, ε is a parameter set artificially.

其它步骤及参数与具体实施方式一或二相同。The other steps and parameters are the same as those in the first or second embodiment.

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述步骤4中判断第k次工作过程是否发生容量再生现象；Specific implementation method 4: This implementation method is different from any one of specific implementation methods 1 to 3 in that in step 4, it is determined whether the capacity regeneration phenomenon occurs during the kth working process;

本发明提出了一种基于Wilcoxon秩和检验的方法，对无迹粒子滤波(UnscentedKalman Filter，UPF)中粒子的先验概率密度分布和后验概率密度分布之间的差异进行衡量，由此准确检测电池容量的再生点，减少由于电池容量突增而导致的剩余使用寿命预测误差；The present invention proposes a method based on Wilcoxon rank sum test to measure the difference between the prior probability density distribution and the posterior probability density distribution of particles in the unscented Kalman Filter (UPF), thereby accurately detecting the regeneration point of the battery capacity and reducing the remaining service life prediction error caused by the sudden increase of the battery capacity;

具体过程为：The specific process is:

在步骤2中基于UPF的自适应状态更新模型中，通过公式(1)—(19)可以获得粒子在UKF算法中的后验概率密度分布，由此作为粒子滤波的先验概率密度分布。接下来，通过公式(20)—(24)可以获得经过粒子滤波的后验概率密度分布。为了方便表达，在经过UKF算法更新后的电池容量，详细地说，经过公式(18)求解的状态变量代入到公式(3)获得电池容量的集合为

Figure BDA0003605634160000181
在经过UPF算法更新后的电池容量，详细地说，经过公式(23)求解的状态变量代入到公式(3)获得电池容量的集合为
Figure BDA0003605634160000182
In the adaptive state update model based on UPF in step 2, the posterior probability density distribution of the particle in the UKF algorithm can be obtained by formulas (1)-(19), which is used as the prior probability density distribution of the particle filter. Next, the posterior probability density distribution after the particle filter can be obtained by formulas (20)-(24). For the convenience of expression, the battery capacity after the UKF algorithm is updated. In detail, the state variables solved by formula (18) are substituted into formula (3) to obtain the set of battery capacities:
Figure BDA0003605634160000181
In detail, the battery capacity after the UPF algorithm updates is obtained by substituting the state variables solved by formula (23) into formula (3) to obtain the set of battery capacities:
Figure BDA0003605634160000182

由于电池在静置一段时间后将发生容量再生的现象，这将导致

Figure BDA0003605634160000183
和
Figure BDA0003605634160000184
之间的分布产生一定的分布差异。Wilcoxon秩和检验的作用在于衡量两组数据样本之间的分布差异。其优势是无需对数据样本的分布作任何的特殊假设，从而能够适用于一些复杂的数据分布情况。为此，本发明将采用Wilcoxon秩和检验的方式衡量
Figure BDA0003605634160000185
和
Figure BDA0003605634160000186
之间的分布差异，提出假设H0如下所示：Since the battery will regenerate capacity after being left idle for a period of time, this will result in
Figure BDA0003605634160000183
and
Figure BDA0003605634160000184
The distribution between the two groups produces a certain distribution difference. The function of the Wilcoxon rank sum test is to measure the distribution difference between two groups of data samples. Its advantage is that it does not require any special assumptions about the distribution of the data samples, so it can be applied to some complex data distribution situations. To this end, the present invention will use the Wilcoxon rank sum test to measure
Figure BDA0003605634160000185
and
Figure BDA0003605634160000186
The distribution difference between them makes the hypothesis H 0 as follows:

假设H0：两个总体

Figure BDA0003605634160000187
和
Figure BDA0003605634160000188
在显著性水平α下具备相同的概率密度分布；Hypothesis H 0 : Two populations
Figure BDA0003605634160000187
and
Figure BDA0003605634160000188
Have the same probability density distribution at the significance level α;

为了检验假设H0，Wilcoxon秩和检验可以总结为以下的步骤：To test the hypothesis H 0 , the Wilcoxon rank sum test can be summarized as follows:

(1)将两个总体

Figure BDA0003605634160000189
和
Figure BDA00036056341600001810
混合获得总体Qk，依照电池容量的数值从小到大进行统一编号，定义每个电池容量样本对应的序数为秩；(1) The two populations
Figure BDA0003605634160000189
and
Figure BDA00036056341600001810
The population Q k is obtained by mixing, and the battery capacity values are uniformly numbered from small to large, and the ordinal number corresponding to each battery capacity sample is defined as the rank;

(2)计算样本总体

Figure BDA00036056341600001811
的样本所对应的秩之和T；(2) Calculate the sample population
Figure BDA00036056341600001811
The sum of the ranks corresponding to the samples of T;

(3)根据粒子总数Ns和显著性水平γ查阅秩和检验表，可以获得秩和下限T1、秩和上限T2以及相应的置信度水平p；(3) According to the total number of particles Ns and the significance level γ, the rank sum test table can be consulted to obtain the rank sum lower limit T1 , the rank sum upper limit T2 and the corresponding confidence level p;

(4)设置区间[T1,T2]，如果T在区间内部，也就是p>γ，则接受假设H0，认为两个总体

Figure BDA0003605634160000191
和
Figure BDA0003605634160000192
在水平α下无显著性差异；否则拒绝假设H0，认为两个总体
Figure BDA0003605634160000193
和
Figure BDA0003605634160000194
在水平α下存在显著性差异；(4) Set the interval [T 1 ,T 2 ]. If T is within the interval, that is, p>γ, then accept the hypothesis H 0 and assume that the two populations are
Figure BDA0003605634160000191
and
Figure BDA0003605634160000192
There is no significant difference at level α; otherwise, the hypothesis H 0 is rejected and the two populations are considered
Figure BDA0003605634160000193
and
Figure BDA0003605634160000194
There is a significant difference at level α;

总结而言，一旦假设H0被拒绝，则说明

Figure BDA0003605634160000195
和
Figure BDA0003605634160000196
的分布存在较大的差异，这也就意味着在第k个工作循环中电池产生了容量再生的现象。In summary, once the hypothesis H 0 is rejected, it means
Figure BDA0003605634160000195
and
Figure BDA0003605634160000196
There is a big difference in the distribution of , which means that the battery has capacity regeneration in the kth working cycle.

其它步骤及参数与具体实施方式一至三之一相同。The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述步骤5中根据步骤2中构建的电池退化模型和当前的电池容量值，通过趋势外推求解电池剩余使用寿命；Specific implementation method 5: This implementation method is different from any one of specific implementation methods 1 to 4 in that in step 5, the remaining service life of the battery is solved by trend extrapolation based on the battery degradation model constructed in step 2 and the current battery capacity value;

根据步骤2中构建的电池退化模型和当前的电池容量值，计算电池容量的置信区间；Calculate the confidence interval of the battery capacity based on the battery degradation model constructed in step 2 and the current battery capacity value;

具体过程为：The specific process is:

电池的剩余使用寿命预测值即为当容量首次退化到失效阈值时所经过的循环数目；The predicted value of the remaining useful life of the battery is the number of cycles when the capacity first degrades to the failure threshold;

根据步骤2中的电池退化模型的状态变量值和当前的电池容量值yk可以通过趋势外推的方式求解电池的剩余使用寿命，具体过程为：According to the state variable value of the battery degradation model in step 2 and the current battery capacity value yk, the remaining service life of the battery can be solved by trend extrapolation. The specific process is:

1)、令预测的起始点v＝k，将公式(23)状态变量期望(期望即为均值)带入公式

Figure BDA0003605634160000197
计算电池容量；1) Let the starting point of the prediction be v = k, and substitute the expected value of the state variable (the expected value is the mean value) in formula (23) into the formula
Figure BDA0003605634160000197
Calculate battery capacity;

若检测到第k个工作循环中电池产生了容量再生的现象，则需将预测的起始点v＝k替代为v＝k-1；If it is detected that the battery has regenerated capacity in the kth working cycle, the predicted starting point v=k needs to be replaced by v=k-1;

2)、如果计算出的电池容量yv小于等于失效阈值(失效阈值设置的)，则执行3)；2) If the calculated battery capacity y v is less than or equal to the failure threshold (the failure threshold is set), execute 3);

如果计算出的电池容量yv大于失效阈值，则令v＝v+1，通过公式

Figure BDA0003605634160000198
计算电池容量，重复执行2)，直到解算出来的yv小于等于电池容量的失效阈值；If the calculated battery capacity y v is greater than the failure threshold, let v = v + 1, and use the formula
Figure BDA0003605634160000198
Calculate the battery capacity and repeat step 2) until the calculated y v is less than or equal to the failure threshold of the battery capacity;

3)、通过公式(46)获得电池剩余使用寿命：3) Obtain the remaining battery life through formula (46):

RULk＝v–k (46)RUL k = v–k (46)

4)、如果计算出的电池容量yv大于失效阈值，获取Ns个粒子的均值和标准差，均值加减正负三倍标准差作为电池容量置信区间上下线。4) If the calculated battery capacity y v is greater than the failure threshold, obtain the mean and standard deviation of N s particles, and use the mean plus or minus three times the standard deviation as the upper and lower limits of the battery capacity confidence interval.

假设获得了电池在第k次过程的电池容量数据yk，然后通过上述的UPF，EM等算法步骤，最后便是开启了自第k个点开始的电池容量预测，此时第k个工作循环便是预测的起始点。如果检测到第k个工作循环中电池产生了容量再生的现象，则需要将预测的起始点yk替代为yk-1，以减少预测误差。Assuming that the battery capacity data y k of the battery in the kth process is obtained, then through the above-mentioned UPF, EM and other algorithm steps, the battery capacity prediction starting from the kth point is finally started. At this time, the kth working cycle is the starting point of the prediction. If it is detected that the battery has a capacity regeneration phenomenon in the kth working cycle, the starting point of the prediction y k needs to be replaced by y k-1 to reduce the prediction error.

其它步骤及参数与具体实施方式一至四之一相同。The other steps and parameters are the same as those in Specific Embodiments 1 to 4.

评估剩余使用寿命的在线预测效果：采用均方根误差(Root Mean Square Error，RMSE)和绝对平均误差(Mean Absolute Error，MAE)这两个指标评估本发明所提出的电池剩余使用寿命预测方法的效果。均方根误差和绝对平均误差和的表达式如式(47)—(48)所示：Evaluation of the online prediction effect of the remaining service life: The root mean square error (RMSE) and the mean absolute error (MAE) are used to evaluate the effect of the battery remaining service life prediction method proposed in the present invention. The expressions of the root mean square error and the mean absolute error are shown in equations (47)-(48):

Figure BDA0003605634160000201
Figure BDA0003605634160000201

Figure BDA0003605634160000202
Figure BDA0003605634160000202

其中，N为测试数据的样本数量，a为样本的序号，RULpa和RULta分别为第a个样本剩余使用寿命的预测值和真实值。均方根误差和绝对平均误差的值越小，则说明本发明提出的剩余使用寿命预测方法效果越好。Wherein, N is the number of samples of the test data, a is the serial number of the sample, RUL pa and RUL ta are the predicted value and true value of the remaining useful life of the ath sample, respectively. The smaller the values of the root mean square error and the absolute mean error, the better the effect of the remaining useful life prediction method proposed in the present invention.

采用以下实施例验证本发明的有益效果：The following examples are used to verify the beneficial effects of the present invention:

实施例一：Embodiment 1:

本发明采用NASA Ames预测中心电提供的电池数据集验证所提出的基于EM-UPF-W的剩余使用寿命预测方法。该数据集包含了18650型号的4个编号为B0005、B0006、B0007和B0018的锂离子电池在24℃的条件下工作所产生的数据。这4个额定容量为2Ah的电池工作模式包含了充电模式和放电模式。对于充电过程，电池以1.5A的恒定电流进行充电，直到电池电压上升到4.2V，然后转变成恒压充电，直到充电电流下降到20mA。对于放电模式，电池以恒定的2A电流放电，直到4个电池两端的电压分别下降到不同的阈值。在此基础上，NASA电池数据集提供了每一个充放电循环结束时的电池容量，因此可以直接用于电池剩余使用寿命的动态预测。考虑到不同电池的工作条件不同，这4个电池的退化阈值分别设定为1.4Ah,1.22Ah,1.6Ah和1.4Ah。这4个电池的电池容量衰减曲线如图2所示。The present invention uses the battery data set provided by NASA Ames Prediction Center to verify the proposed remaining service life prediction method based on EM-UPF-W. The data set contains data generated by four 18650 lithium-ion batteries numbered B0005, B0006, B0007 and B0018 working under 24°C. The working modes of these four batteries with a rated capacity of 2Ah include charging mode and discharging mode. For the charging process, the battery is charged at a constant current of 1.5A until the battery voltage rises to 4.2V, and then converted to constant voltage charging until the charging current drops to 20mA. For the discharge mode, the battery is discharged at a constant current of 2A until the voltages at both ends of the four batteries drop to different thresholds. On this basis, the NASA battery data set provides the battery capacity at the end of each charge and discharge cycle, so it can be directly used for the dynamic prediction of the remaining service life of the battery. Considering the different working conditions of different batteries, the degradation thresholds of these four batteries are set to 1.4Ah, 1.22Ah, 1.6Ah and 1.4Ah respectively. The battery capacity decay curves of these four batteries are shown in Figure 2.

步骤1，提取在第k次工作过程中的电池容量数据：NASA电池数据集提供了每一个充放电循环结束时的电池容量，由此可以获得相应循环下的电池容量数据，作为电池剩余使用寿命预测方法的输入。Step 1, extract the battery capacity data during the kth working process: The NASA battery data set provides the battery capacity at the end of each charge and discharge cycle, from which the battery capacity data under the corresponding cycle can be obtained as the input of the battery remaining service life prediction method.

步骤2，构建电池的退化模型：仅仅依据步骤1中的电池容量数据，构建基于无迹粒子滤波的动态退化模型描述电池的退化过程。通过对于粒子初始化、UKF状态变量和协方差的估计、粒子权重更新以及重采样等步骤获得隐含的参数变量。Step 2, build a battery degradation model: Based only on the battery capacity data in step 1, build a dynamic degradation model based on unscented particle filtering to describe the battery degradation process. Obtain implicit parameter variables through steps such as particle initialization, UKF state variable and covariance estimation, particle weight update, and resampling.

步骤3，自适应估计退化模型中的过程噪声和测量噪声：采用基于最大期望算法自适应估计退化模型中的过程噪声和测量噪声。具体而言，针对基于无迹粒子滤波的动态退化模型中的过程噪声和测量噪声，采用最大期望算法进行自适应估计，利用电池从初始状态到当前时刻的所有容量数据对过程噪声和测量噪声进行递归更新。本发明以编号为B0005的电池为例展示算法对于过程噪声和测量噪声的自适应参数估计结果，如图3a、3b、3c、3d、3e所示。从图3a、3b、3c、3d、3e中可以看出，所有的模型参数在20个工作循环以内便可以实现收敛，因此，本发明所提出的基于最大期望的的参数估计方法有着很好的收敛效果。Step 3, adaptively estimate the process noise and measurement noise in the degradation model: adaptively estimate the process noise and measurement noise in the degradation model based on the maximum expectation algorithm. Specifically, for the process noise and measurement noise in the dynamic degradation model based on the untraceable particle filter, the maximum expectation algorithm is used for adaptive estimation, and the process noise and measurement noise are recursively updated using all the capacity data of the battery from the initial state to the current moment. The present invention takes the battery numbered B0005 as an example to demonstrate the adaptive parameter estimation results of the algorithm for process noise and measurement noise, as shown in Figures 3a, 3b, 3c, 3d, and 3e. It can be seen from Figures 3a, 3b, 3c, 3d, and 3e that all model parameters can converge within 20 working cycles. Therefore, the parameter estimation method based on maximum expectation proposed in the present invention has a good convergence effect.

步骤4，判断第k次工作过程是否发生容量再生现象：采用基于Wilcoxon秩和检验的方法对于无迹粒子滤波中粒子的先验概率密度分布和后验概率密度分布之间的差异进行衡量，由此准确检测电池容量的再生点，减少由于电池容量突增而导致的剩余使用寿命预测误差。本发明图4a、4b以编号为B0005的电池为例展示Wilcoxon秩和检验的方法对于电池容量再生点的检验结果。Step 4, determine whether capacity regeneration occurs during the kth working process: use a method based on the Wilcoxon rank sum test to measure the difference between the prior probability density distribution and the posterior probability density distribution of particles in the unscented particle filter, thereby accurately detecting the regeneration point of the battery capacity and reducing the remaining service life prediction error caused by the sudden increase in battery capacity. Figures 4a and 4b of the present invention take the battery numbered B0005 as an example to show the test results of the Wilcoxon rank sum test method for the battery capacity regeneration point.

步骤5，预测电池剩余使用寿命：根据步骤2中的电池退化模型和当前的电池容量值，通过趋势外推求解电池剩余使用寿命。需要注意的是，如果检测到第k个工作循环中电池产生了容量再生的现象，则需要将预测的起始点yk替代为yk-1，以减少预测误差。最后，分别将每个粒子分别代入模型便可以获得电池容量的置信区间。Step 5, predict the remaining battery life: According to the battery degradation model in step 2 and the current battery capacity value, the remaining battery life is solved by trend extrapolation. It should be noted that if the battery capacity regeneration phenomenon is detected in the kth working cycle, the starting point y k of the prediction needs to be replaced by y k-1 to reduce the prediction error. Finally, each particle is substituted into the model to obtain the confidence interval of the battery capacity.

图5a、5b、5c、5d展示了本发明提出的方法在NASA电池数据集上4个电池的剩余使用寿命真实值和预测值。可以看出，本发明提出的方法中的剩余使用寿命预测值和真实值十分贴近，由此可知本发明所提出的方法有着很好的预测效果。Figures 5a, 5b, 5c, and 5d show the true and predicted values of the remaining service life of four batteries on the NASA battery data set using the method proposed in the present invention. It can be seen that the predicted value of the remaining service life in the method proposed in the present invention is very close to the true value, which shows that the method proposed in the present invention has a good prediction effect.

图6a、6b、6c、6d展示了电池容量99％的置信区间，可以看出电池容量的真实值基本都落在置信区间以内，由此可知本发明提出的方法能够很好地描述电池在退化过程中的不确定性。Figures 6a, 6b, 6c, and 6d show the 99% confidence interval of the battery capacity. It can be seen that the true value of the battery capacity basically falls within the confidence interval. Therefore, it can be seen that the method proposed in the present invention can well describe the uncertainty of the battery in the degradation process.

步骤6，评估剩余使用寿命的在线预测效果：采用均方根误差(Root Mean SquareError，RMSE)和绝对平均误差(Mean Absolute Error，MAE)这两个指标评估本发明所提出的电池剩余使用寿命预测方法的效果。Step 6, evaluating the online prediction effect of the remaining service life: using the two indicators of root mean square error (RMSE) and mean absolute error (MAE) to evaluate the effect of the battery remaining service life prediction method proposed in the present invention.

表1剩余使用寿命的预测结果Table 1 Prediction results of remaining service life

Figure BDA0003605634160000211
Figure BDA0003605634160000211

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。The present invention may also have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims (2)

1.基于最大期望-无迹粒子滤波的电池剩余使用寿命预测方法，其特征在于：所述方法具体过程为：1. A method for predicting the remaining service life of a battery based on maximum expectation-unscented particle filtering, characterized in that: the specific process of the method is as follows: 步骤1、提取第k次工作过程中的电池容量数据；Step 1, extracting the battery capacity data during the kth working process; 步骤2、根据步骤1中提取的第k次工作过程中的电池容量数据，构建基于无迹粒子滤波的动态电池退化模型；Step 2: According to the battery capacity data of the kth working process extracted in step 1, a dynamic battery degradation model based on unscented particle filtering is constructed; 步骤3、采用基于最大期望算法自适应地估计电池退化模型中的过程噪声和测量噪声；Step 3, adaptively estimating process noise and measurement noise in the battery degradation model using a maximum expectation algorithm; 过程噪声和测量噪声用于第k+1次工作过程中，构建基于无迹粒子滤波的动态电池退化模型；The process noise and measurement noise are used to construct a dynamic battery degradation model based on unscented particle filtering during the k+1th working process; 步骤4、判断第k次工作过程是否发生容量再生现象；Step 4, determining whether capacity regeneration occurs during the kth working process; 步骤5、根据步骤2中构建的电池退化模型和当前的电池容量值，通过趋势外推求解电池剩余使用寿命；Step 5: According to the battery degradation model constructed in step 2 and the current battery capacity value, the remaining battery life is solved by trend extrapolation; 根据步骤2中构建的电池退化模型和当前的电池容量值，计算电池容量的置信区间；Calculate the confidence interval of the battery capacity based on the battery degradation model constructed in step 2 and the current battery capacity value; 所述步骤2中根据步骤1中提取的第k次工作过程中的电池容量数据，构建基于无迹粒子滤波的动态电池退化模型；具体过程为：In 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 step 1; the specific process is: 基于无迹粒子滤波的电池退化模型可以采用状态空间的形式表示成公式(1)所示的形式：The battery degradation model based on unscented particle filtering can be expressed in the form of state space as shown in formula (1):
Figure FDA0004171951290000011
Figure FDA0004171951290000011
 其中，xk为第k个工作循环下的状态变量，xk＝[ak,bk,ck,dk]T，ak为第k个工作循环下的状态变量的第一分量，bk为第k个工作循环下的状态变量的第二分量，ck为为第k个工作循环下的状态变量的第三分量，dk为第k个工作循环下的状态变量的第四分量，T为转置，k表示电池工作循环次数，xk-1为第k-1个工作循环下的状态变量，yk为电池容量，uk为过程噪声项，uk＝[ua,ub,uc,ud]T，ua为第k个工作循环下状态变量的第一分量的噪声项，ub为第k个工作循环下状态变量的第二分量的噪声项，uc为第k个工作循环下状态变量的第三分量的噪声项，ud为第k个工作循环下状态变量的第四分量的噪声项，vk为测量噪声项，vk∈R1×1，1×1为长宽为1的矩阵，f()代表状态转移方程的映射关系，h()代表测量方程的映射关系；Wherein, xk is the state variable under the kth working cycle, xk = [ ak , bk , ck , dk ] T , ak is the first component of the state variable under the kth working cycle, bk is the second component of the state variable under the kth working cycle, ck is the third component of the state variable under the kth working cycle, dk is the fourth component of the state variable under the kth working cycle, T is the transpose, k represents the number of battery working cycles, xk -1 is the state variable under the k-1th working cycle, yk is the battery capacity, uk is the process noise term, uk = [ ua , ub , uc , ud ] T , ua is the noise term of the first component of the state variable under the kth working cycle, ub is the noise term of the second component of the state variable under the kth working cycle, uc is the noise term of the third component of the state variable under the kth working cycle, ud is the noise term of the fourth component of the state variable under the kth working cycle, vk is the measurement noise term, vk∈R 1×1 , 1×1 is a matrix with a length and width of 1, f() represents the mapping relationship of the state transfer equation, and h() represents the mapping relationship of the measurement equation; 公式(1)中的状态转移方程xk＝f(xk-1)+uk可以写成公式(2)中的形式：The state transfer equation x k =f(x k-1 )+ uk in formula (1) can be written in the form of formula (2):
Figure FDA0004171951290000021
Figure FDA0004171951290000021
 其中，
Figure FDA0004171951290000022
为过程变量噪声项ua的方差，
Figure FDA0004171951290000023
为过程变量噪声项ub的方差，
Figure FDA0004171951290000024
为过程变量噪声项uc的方差，
Figure FDA0004171951290000025
为过程变量噪声项ud的方差；ak-1为第k-1个工作循环下的状态变量的第一分量，bk-1为第k-1个工作循环下的状态变量的第二分量，ck-1为第k-1个工作循环下的状态变量的第三分量，dk-1为第k-1个工作循环下的状态变量的第四分量；
Figure FDA0004171951290000026
表示ua服从均值为0，方差为
Figure FDA0004171951290000027
的正态分布，
Figure FDA0004171951290000028
为均值为0，方差为
Figure FDA0004171951290000029
的正态分布，
Figure FDA00041719512900000210
为均值为0，方差为
Figure FDA00041719512900000211
的正态分布，
Figure FDA00041719512900000212
为为均值为0，方差为
Figure FDA00041719512900000213
的正态分布，
Figure FDA00041719512900000214
为均值为0，方差为
Figure FDA00041719512900000215
的正态分布；in,
Figure FDA0004171951290000022
is the variance of the process variable noise term u a ,
Figure FDA0004171951290000023
is the variance of the process variable noise term u b ,
Figure FDA0004171951290000024
is the variance of the process variable noise term u c ,
Figure FDA0004171951290000025
is the variance of the process variable noise term ud ; a k-1 is the first component of the state variable under the k-1th working cycle, b k-1 is the second component of the state variable under the k-1th working cycle, c k-1 is the third component of the state variable under the k-1th working cycle, and d k-1 is the fourth component of the state variable under the k-1th working cycle;
Figure FDA0004171951290000026
It means that u a has a mean of 0 and a variance of
Figure FDA0004171951290000027
The normal distribution of
Figure FDA0004171951290000028
The mean is 0 and the variance is
Figure FDA0004171951290000029
The normal distribution of
Figure FDA00041719512900000210
The mean is 0 and the variance is
Figure FDA00041719512900000211
The normal distribution of
Figure FDA00041719512900000212
The mean is 0 and the variance is
Figure FDA00041719512900000213
The normal distribution of
Figure FDA00041719512900000214
The mean is 0 and the variance is
Figure FDA00041719512900000215
Normal distribution of 测量方程yk＝h(xk)+vk可以写成公式(3)中的形式：The measurement equation y k =h(x k )+v k can be written in the form of formula (3):
Figure FDA00041719512900000216
Figure FDA00041719512900000216
 其中，
Figure FDA00041719512900000217
为测量噪声的方差；
Figure FDA00041719512900000218
为均值为0，方差为
Figure FDA00041719512900000219
的正态分布；in,
Figure FDA00041719512900000217
is the variance of the measurement noise;
Figure FDA00041719512900000218
The mean is 0 and the variance is
Figure FDA00041719512900000219
Normal distribution of (1)粒子初始化，过程为：(1) Particle initialization, the process is: 根据初始状态变量x0，在初始状态变量x0的正态分布p(x0)中通过蒙特卡洛方法产生粒子集合
Figure FDA00041719512900000220
Ns为粒子的总数，且每个粒子的初始化权重值为
Figure FDA00041719512900000221
According to the initial state variable x 0 , a particle set is generated by the Monte Carlo method in the normal distribution p(x 0 ) of the initial state variable x 0
Figure FDA00041719512900000220
N s is the total number of particles, and the initialization weight value of each particle is
Figure FDA00041719512900000221
 初始化粒子的均值
Figure FDA00041719512900000222
以及初始化粒子的方差S0可以写成公式(4)—(5)所示的形式：Initialize the mean of the particles
Figure FDA00041719512900000222
And the variance S 0 of the initialized particle can be written as shown in formulas (4)-(5):
Figure FDA00041719512900000223
Figure FDA00041719512900000223
Figure FDA00041719512900000224
Figure FDA00041719512900000224
 其中，i为粒子的编号，
Figure FDA00041719512900000225
为初始化粒子的均值，S0为初始化粒子的方差，E()表示期望运算符，T表示矩阵转置运算符，r0为初始化粒子；Where i is the number of the particle,
Figure FDA00041719512900000225
is the mean of the initialized particle, S 0 is the variance of the initialized particle, E() represents the expectation operator, T represents the matrix transposition operator, and r 0 is the initialized particle; (2)通过UKF算法计算粒子所表征的状态变量rk的后验期望
Figure FDA0004171951290000031
和协方差
Figure FDA0004171951290000032
过程为：(2) Calculate the posterior expectation of the state variable r k represented by the particle through the UKF algorithm
Figure FDA0004171951290000031
and covariance
Figure FDA0004171951290000032
The process is: 在完成k-1次迭代过程后，构建各粒子的Sigma采样点集
Figure FDA0004171951290000033
如公式(6)所示：After completing k-1 iterations, construct the Sigma sampling point set for each particle
Figure FDA0004171951290000033
As shown in formula (6):
Figure FDA0004171951290000034
Figure FDA0004171951290000034
 其中，
Figure FDA0004171951290000035
为第i个粒子在k-1次工作过程中的均值向量，
Figure FDA0004171951290000036
为第i个粒子在k-1次工作过程中的协方差矩阵，参数α为UKF中的系数值；in,
Figure FDA0004171951290000035
is the mean vector of the i-th particle in the k-1 working process,
Figure FDA0004171951290000036
is the covariance matrix of the ith particle in the k-1 working process, and the parameter α is the coefficient value in UKF; 各粒子的Sigma点集的样本点个数为2n+1，n为状态变量的数目，各粒子的Sigma点集的样本点的权重表达式分别如公式(7)—(9)所示：The number of sample points of the Sigma point set of each particle is 2n+1, where n is the number of state variables. The weight expressions of the sample points of the Sigma point set of each particle are shown in formulas (7) to (9):
Figure FDA0004171951290000037
Figure FDA0004171951290000037
Figure FDA0004171951290000038
Figure FDA0004171951290000038
Figure FDA0004171951290000039
Figure FDA0004171951290000039
 其中，
Figure FDA00041719512900000310
为第i个粒子第1个Sigma点集的样本点所代表的过程变量权重，
Figure FDA00041719512900000311
为第i个粒子第1个Sigma点集的样本点所代表的测量变量权重，
Figure FDA00041719512900000312
为第i个粒子第j个Sigma点集的样本点所代表的过程变量权重，
Figure FDA00041719512900000313
为第i个粒子第j个Sigma点集的样本点所代表的测量变量权重，参数β,λ为UKF中的系数值，n为状态变量的数目；in,
Figure FDA00041719512900000310
is the weight of the process variable represented by the sample point of the first Sigma point set of the i-th particle,
Figure FDA00041719512900000311
is the weight of the measured variable represented by the sample point of the first Sigma point set of the i-th particle,
Figure FDA00041719512900000312
is the weight of the process variable represented by the sample point of the jth Sigma point set of the i-th particle,
Figure FDA00041719512900000313
is the weight of the measured variable represented by the sample point of the jth Sigma point set of the i-th particle, the parameters β and λ are the coefficient values in UKF, and n is the number of state variables; 在此基础上，UKF的状态估计表达式可以写成公式(10)—(11)所示的形式：On this basis, the state estimation expression of UKF can be written as shown in formulas (10)-(11):
Figure FDA00041719512900000314
Figure FDA00041719512900000314
Figure FDA00041719512900000315
Figure FDA00041719512900000315
 其中，
Figure FDA00041719512900000316
为第k次工作过程第i个粒子中第j个Sigma点所代表的先验状态向量，f()为状态转移方程的映射关系，
Figure FDA00041719512900000317
为第k-1次工作过程第i个粒子中第j个Sigma点所代表的后验状态向量，
Figure FDA00041719512900000318
为第k-1次工作过程第i个粒子所代表的噪声变量，
Figure FDA00041719512900000319
为第k次工作过程第i个粒子的加权先验状态向量；in,
Figure FDA00041719512900000316
is the prior state vector represented by the jth Sigma point in the i-th particle in the k-th working process, f() is the mapping relationship of the state transfer equation,
Figure FDA00041719512900000317
is the posterior state vector represented by the jth Sigma point in the i-th particle in the k-1th working process,
Figure FDA00041719512900000318
is the noise variable represented by the i-th particle in the k-1th working process,
Figure FDA00041719512900000319
is the weighted prior state vector of the ith particle in the kth working process; 协方差估计的表达式则可以写成公式(12)所示的形式：The expression of covariance estimation can be written as shown in formula (12):
Figure FDA0004171951290000041
Figure FDA0004171951290000041
 其中，
Figure FDA0004171951290000042
为第k次工作过程第i个粒子中状态变量的先验协方差，
Figure FDA0004171951290000043
为第i个粒子第k次工作过程第i个粒子中噪声变量的协方差矩阵，
Figure FDA0004171951290000044
diag表示对角矩阵，
Figure FDA0004171951290000045
为第k次工作过程过程变量噪声项ua的方差，
Figure FDA0004171951290000046
为第k次工作过程过程变量噪声项ub的方差，
Figure FDA0004171951290000047
为第k次工作过程的过程变量噪声项uc的方差，
Figure FDA0004171951290000048
为第k次工作过程过程变量噪声项ud的方差；in,
Figure FDA0004171951290000042
is the prior covariance of the state variables in the ith particle of the kth working process,
Figure FDA0004171951290000043
is the covariance matrix of the noise variable in the ith particle in the kth working process of the ith particle,
Figure FDA0004171951290000044
diag represents a diagonal matrix,
Figure FDA0004171951290000045
is the variance of the process variable noise term u a of the kth working process,
Figure FDA0004171951290000046
is the variance of the process variable noise term u b of the kth working process,
Figure FDA0004171951290000047
is the variance of the process variable noise term u c of the kth working process,
Figure FDA0004171951290000048
is the variance of the process variable noise term ud in the kth working process; 对于前向的一步预测，有公式(13)—(14)所示的形式：For the one-step forward prediction, we have the form shown in formulas (13)-(14):
Figure FDA0004171951290000049
Figure FDA0004171951290000049
Figure FDA00041719512900000410
Figure FDA00041719512900000410
 其中，
Figure FDA00041719512900000411
为第i个粒子中第j个Sigma点所代表的状态向量，
Figure FDA00041719512900000412
为第i个粒子中第j个Sigma点所代表的UKF前向一步预测值，
Figure FDA00041719512900000413
为第i个粒子所代表的UKF预测值，
Figure FDA00041719512900000414
为第k次工作过程第i个粒子的测量噪声项，h()代表测量方程的映射关系；in,
Figure FDA00041719512900000411
is the state vector represented by the jth Sigma point in the ith particle,
Figure FDA00041719512900000412
is the UKF one-step-forward prediction value represented by the j-th Sigma point in the i-th particle,
Figure FDA00041719512900000413
is the UKF prediction value represented by the i-th particle,
Figure FDA00041719512900000414
is the measurement noise term of the ith particle in the kth working process, and h() represents the mapping relationship of the measurement equation; 计算自协方差矩阵
Figure FDA00041719512900000415
和互协方差矩阵
Figure FDA00041719512900000416
如公式(15)—(16)所示：Compute the autocovariance matrix
Figure FDA00041719512900000415
and the cross-covariance matrix
Figure FDA00041719512900000416
As shown in formulas (15)-(16):
Figure FDA00041719512900000417
Figure FDA00041719512900000417
Figure FDA00041719512900000418
Figure FDA00041719512900000418
 其中，
Figure FDA00041719512900000419
为第k次工作过程第i个粒子的噪声方差；in,
Figure FDA00041719512900000419
is the noise variance of the i-th particle in the k-th working process; 由此可以获得卡尔曼滤波增益向量Kk如公式(17)所示：Thus, the Kalman filter gain vector K k can be obtained as shown in formula (17):
Figure FDA00041719512900000420
Figure FDA00041719512900000420
 其中，
Figure FDA00041719512900000421
为第i个粒子卡尔曼滤波增益向量，
Figure FDA00041719512900000422
为自协方差矩阵，
Figure FDA00041719512900000423
为互协方差矩阵，随着第k次迭代过程中新的电池容量数据的加入，UKF的状态更新和协方差更新可以表达成公式(18)—(19)所示的形式：in,
Figure FDA00041719512900000421
is the Kalman filter gain vector of the i-th particle,
Figure FDA00041719512900000422
is the autocovariance matrix,
Figure FDA00041719512900000423
is the mutual covariance matrix. With the addition of new battery capacity data in the kth iteration, the state update and covariance update of UKF can be expressed as shown in formulas (18)-(19):
Figure FDA0004171951290000051
Figure FDA0004171951290000051
Figure FDA0004171951290000052
Figure FDA0004171951290000052
 其中，
Figure FDA0004171951290000053
为通过UKF算法计算粒子的状态变量rk的后验期望，
Figure FDA0004171951290000054
为通过UKF算法计算粒子的状态变量rk的协方差；in,
Figure FDA0004171951290000053
To calculate the posterior expectation of the particle state variable r k by the UKF algorithm,
Figure FDA0004171951290000054
is the covariance of the particle state variable r k calculated by the UKF algorithm; 根据公式(18)—(19)中的结果，对
Figure FDA0004171951290000055
求取均值获得
Figure FDA0004171951290000056
对
Figure FDA0004171951290000057
求取均值获得
Figure FDA0004171951290000058
基于
Figure FDA0004171951290000059
和
Figure FDA00041719512900000510
获得在UKF算法下粒子的后验概率密度分布
Figure FDA00041719512900000511
According to the results in formulas (18)-(19),
Figure FDA0004171951290000055
Get the mean
Figure FDA0004171951290000056
right
Figure FDA0004171951290000057
Get the mean
Figure FDA0004171951290000058
based on
Figure FDA0004171951290000059
and
Figure FDA00041719512900000510
Obtain the posterior probability density distribution of particles under the UKF algorithm
Figure FDA00041719512900000511
 (3)粒子权重更新，过程为：(3) Particle weight update, the process is: 根据(2)获得粒子的重要性密度函数可以表示为：
Figure FDA00041719512900000512
According to (2), the importance density function of the particle can be expressed as:
Figure FDA00041719512900000512
 其中，
Figure FDA00041719512900000513
表示第k次工作过程第i个粒子的重要性密度函数，
Figure FDA00041719512900000514
为第i个粒子初始状态到第k-1次工作过程所代表的状态变量，y1:k表示第1次工作过程到第k次工作过程的电池容量数据，
Figure FDA00041719512900000515
为第i个粒子初始状态到第k-1次工作过程所代表的状态变量和第1次工作过程到第k次工作过程的电池容量数据先验条件下
Figure FDA00041719512900000516
的概率密度函数，
Figure FDA00041719512900000517
为均值为
Figure FDA00041719512900000518
方差为
Figure FDA00041719512900000519
的正态分布，～为服从某分布的含义；in,
Figure FDA00041719512900000513
represents the importance density function of the i-th particle in the k-th working process,
Figure FDA00041719512900000514
is the state variable represented by the initial state of the i-th particle to the k-1th working process, y 1:k represents the battery capacity data from the 1st working process to the kth working process,
Figure FDA00041719512900000515
The state variables represented by the initial state of the i-th particle to the k-1th working process and the battery capacity data from the 1st working process to the kth working process under the prior conditions
Figure FDA00041719512900000516
The probability density function of
Figure FDA00041719512900000517
The mean is
Figure FDA00041719512900000518
The variance is
Figure FDA00041719512900000519
Normal distribution, ~ means obeying a certain distribution; 在此基础上，粒子的权重更新表达式可以写成公式(20)中的形式：On this basis, the particle weight update expression can be written as in formula (20):
Figure FDA00041719512900000520
Figure FDA00041719512900000520
 其中，
Figure FDA00041719512900000521
为第i个粒子在k-1个工作循环下的权重，
Figure FDA00041719512900000522
为第i个粒子在k个工作循环下的权重，
Figure FDA00041719512900000523
为yk在
Figure FDA00041719512900000524
条件下的概率密度函数，
Figure FDA00041719512900000525
为
Figure FDA00041719512900000526
在
Figure FDA00041719512900000527
条件下的概率密度函数；in,
Figure FDA00041719512900000521
is the weight of the ith particle under k-1 working cycles,
Figure FDA00041719512900000522
is the weight of the ith particle under k working cycles,
Figure FDA00041719512900000523
For y k
Figure FDA00041719512900000524
The probability density function under the condition,
Figure FDA00041719512900000525
for
Figure FDA00041719512900000526
exist
Figure FDA00041719512900000527
Probability density function under the conditions;
Figure FDA00041719512900000528
在经过粒子的权重归一化，可以获得归一化的粒子权重
Figure FDA00041719512900000529
表达式如下所示：
Figure FDA00041719512900000528
After the particle weight is normalized, the normalized particle weight can be obtained
Figure FDA00041719512900000529
The expression is as follows:
Figure FDA00041719512900000530
Figure FDA00041719512900000530
 (4)粒子重采样，过程为：(4) Particle resampling, the process is: 假设有效粒子数目为Nef，可以定义成公式(22)中的形式：Assuming that the effective number of particles is Nef , it can be defined as in formula (22):
Figure FDA0004171951290000061
Figure FDA0004171951290000061
 设定阀值Nth，如果有效粒子数目Nef小于Nth，则需要进行粒子的重采样；Set the threshold N th . If the number of valid particles Nef is less than N th , resampling of particles is required. 选用随机重采样的方式实现粒子滤波的重采样过程；The random resampling method is used to implement the resampling process of particle filtering; 如果有效粒子数目Nef大于等于Nth，则不需要进行粒子的重采样；If the number of effective particles Nef is greater than or equal to Nth , there is no need to resample the particles; 在完成(4)后，每一个粒子的权重为
Figure FDA0004171951290000062
After completing (4), the weight of each particle is
Figure FDA0004171951290000062
 (5)状态更新，过程为：(5) Status update: The process is as follows: 在重采样获得的粒子权重的基础上，即可通过公式(23)—(24)获得无迹粒子滤波算法后验条件下的状态变量期望
Figure FDA0004171951290000063
和协方差
Figure FDA0004171951290000064
Based on the particle weights obtained by resampling, the expected state variables of the unscented particle filter algorithm under the posterior condition can be obtained through formulas (23)-(24):
Figure FDA0004171951290000063
and covariance
Figure FDA0004171951290000064
Figure FDA0004171951290000065
Figure FDA0004171951290000065
Figure FDA0004171951290000066
Figure FDA0004171951290000066
 其中，
Figure FDA0004171951290000067
为UPF算法后验条件下的状态变量期望，
Figure FDA0004171951290000068
为UPF算法后验条件下的协方差；in,
Figure FDA0004171951290000067
is the expectation of the state variables under the posterior condition of the UPF algorithm,
Figure FDA0004171951290000068
is the covariance under the posterior condition of the UPF algorithm; 所述步骤3中采用基于最大期望算法自适应地估计电池退化模型中的过程噪声和测量噪声；具体过程为：In step 3, the process noise and measurement noise in the battery degradation model are adaptively estimated based on the maximum expectation algorithm; the specific process is: 过程噪声和测量噪声所构成的未知参数向量可以表示为：
Figure FDA0004171951290000069
The unknown parameter vector composed of process noise and measurement noise can be expressed as:
Figure FDA0004171951290000069
 在已知自初始到第k个工作循环的电池容量数据y1:k＝[y1,y2,…,yk]的前提下，可以构建联合对数似然函数Lk(Ξ)如公式(25)所示：Under the premise that the battery capacity data y 1:k =[y 1 ,y 2 ,…,y k ] from the initial to the kth working cycle is known, the joint log-likelihood function L k (Ξ) can be constructed as shown in formula (25):
Figure FDA00041719512900000610
Figure FDA00041719512900000610
 其中，rt (i)为第t次工作过程第i个粒子的状态变量，
Figure FDA0004171951290000071
Figure FDA0004171951290000072
为第i个粒子的第t次工作过程中的过程变量值，Ns为粒子的总数，且每个粒子的初始化权重值为
Figure FDA0004171951290000073
i为粒子的编号，
Figure FDA0004171951290000074
为第1次工作过程到第k次工作过程第i个粒子所表征的状态变量，
Figure FDA0004171951290000075
为未知参数条件下的第1次工作过程到第k次工作过程第i个粒子所表征的状态变量与第1次工作过程到第k次工作过程的电池容量数据所构成的联合概率密度函数，
Figure FDA0004171951290000076
为未知参数条件下第1次工作过程到第k次工作过程第i个粒子所表征的状态变量的条件概率密度函数，
Figure FDA0004171951290000077
为未知参数条件和第1次工作过程到第k次工作过程第i个粒子所表征的状态变量条件下第1次工作过程到第k次工作过程的电池容量数据的条件概率密度函数，
Figure FDA0004171951290000078
为未知参数条件和第t-1次工作过程第i个粒子状态变量条件下第t次工作过程第i个粒子状态变量的条件概率密度函数，
Figure FDA0004171951290000079
为未知参数条件和第t次工作过程第i个粒子状态变量条件下第t次工作过程电池容量数据的条件概率密度函数，
Figure FDA00041719512900000710
第t-1次工作过程第i个粒子状态变量，
Figure FDA00041719512900000711
为未知参数条件下的第1次工作过程到第k次工作过程第i个粒子所表征的状态变量与第1次工作过程到第k次工作过程的电池容量数据的关系，
Figure FDA00041719512900000712
为未知参数条件下第1次工作过程到第k次工作过程第i个粒子所表征的状态变量的关系，
Figure FDA00041719512900000713
为未知参数条件和第1次工作过程到第k次工作过程第i个粒子所表征的状态变量条件下第1次工作过程到第k次工作过程的电池容量数据的关系，
Figure FDA00041719512900000714
未知参数条件和第t-1次工作过程第i个粒子状态变量条件下第t次工作过程第i个粒子状态变量的关系，yt|rt (i),Ξ为未知参数条件和第t次工作过程第i个粒子状态变量条件下第t次工作过程电池容量数据的关系；Among them, r t (i) is the state variable of the i-th particle in the t-th working process,
Figure FDA0004171951290000071
Figure FDA0004171951290000072
is the process variable value of the t-th working process of the i-th particle, N s is the total number of particles, and the initialization weight value of each particle is
Figure FDA0004171951290000073
i is the number of the particle,
Figure FDA0004171951290000074
is the state variable represented by the i-th particle from the 1st working process to the k-th working process,
Figure FDA0004171951290000075
is the joint probability density function composed of the state variables represented by the i-th particle from the 1st working process to the k-th working process under unknown parameter conditions and the battery capacity data from the 1st working process to the k-th working process,
Figure FDA0004171951290000076
is the conditional probability density function of the state variable represented by the i-th particle from the 1st working process to the k-th working process under unknown parameter conditions,
Figure FDA0004171951290000077
is the conditional probability density function of the battery capacity data from the 1st working process to the kth working process under the condition of unknown parameters and the state variables represented by the i-th particle from the 1st working process to the kth working process,
Figure FDA0004171951290000078
is the conditional probability density function of the state variable of the ith particle in the t-th working process under the condition of unknown parameters and the state variable of the ith particle in the t-1th working process,
Figure FDA0004171951290000079
is the conditional probability density function of the battery capacity data of the t-th working process under the condition of unknown parameters and the i-th particle state variable of the t-th working process,
Figure FDA00041719512900000710
The state variable of the i-th particle in the t-1th working process,
Figure FDA00041719512900000711
is the relationship between the state variable represented by the i-th particle from the 1st working process to the k-th working process under unknown parameter conditions and the battery capacity data from the 1st working process to the k-th working process,
Figure FDA00041719512900000712
is the relationship between the state variables represented by the i-th particle from the 1st working process to the k-th working process under unknown parameter conditions,
Figure FDA00041719512900000713
is the relationship between the battery capacity data from the 1st working process to the kth working process under the unknown parameter conditions and the state variable conditions represented by the i-th particle from the 1st working process to the kth working process,
Figure FDA00041719512900000714
The relationship between the unknown parameter condition and the state variable of the i-th particle in the t-1-th working process, y t |r t (i) ,Ξis the relationship 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 in the t-th working process; 基于公式(1)—(3)中状态转移方程和测量方程中的内容，可以获得公式(26)—(27)中的关系表达式：Based on the contents of the state transfer equation and the measurement equation in formulas (1)-(3), the relational expressions in formulas (26)-(27) can be obtained:
Figure FDA00041719512900000715
Figure FDA00041719512900000715
Figure FDA0004171951290000081
Figure FDA0004171951290000081
 其中，
Figure FDA0004171951290000082
为均值为
Figure FDA0004171951290000083
方差为
Figure FDA0004171951290000084
的正态分布，t为工作过程的编号，
Figure FDA0004171951290000085
为第t次工作过程第i个粒子的噪声变量方差，
Figure FDA0004171951290000086
第t次工作过程第i个粒子的过程变量协方差矩阵，
Figure FDA0004171951290000087
diag表示对角矩阵，
Figure FDA0004171951290000088
为分别第t次工作过程过程变量噪声项ua，ub，uc，ud的方差；in,
Figure FDA0004171951290000082
The mean is
Figure FDA0004171951290000083
The variance is
Figure FDA0004171951290000084
The normal distribution of , t is the number of the work process,
Figure FDA0004171951290000085
is the noise variable variance of the ith particle in the tth working process,
Figure FDA0004171951290000086
The process variable covariance matrix of the ith particle in the tth working process is:
Figure FDA0004171951290000087
diag represents a diagonal matrix,
Figure FDA0004171951290000088
are the variances of the process variable noise terms ua , ub , uc , ud of the t-th working process respectively; 将公式(26)—(27)代入到公式(25)中，且忽略常数项，则联合对数似然函数可以进一步表示为公式(28)所示的形式：Substituting formulas (26)-(27) into formula (25) and ignoring the constant term, the joint log-likelihood function can be further expressed as shown in formula (28):
Figure FDA0004171951290000089
Figure FDA0004171951290000089
 其中，∝表示正比于的含义；Among them, ∝ means proportional to; 在已知自初始到第k个工作循环的电池容量数据y1:k＝[y1,y2,…,yk]的前提下，对于在EM算法的第m次迭代过程，自适应退化模型的未知参数估计向量可以表示为：
Figure FDA00041719512900000810
Under the premise that the battery capacity data y 1:k = [y 1 ,y 2 ,…,y k ] from the initial to the kth working cycle is known, for the mth iteration process of the EM algorithm, the unknown parameter estimation vector of the adaptive degradation model can be expressed as:
Figure FDA00041719512900000810
 其中，
Figure FDA00041719512900000811
为第m次迭代过程中第k次工作过程噪声变量ua的方差，
Figure FDA00041719512900000812
为第m次迭代过程中第k次工作过程噪声变量ub的方差，
Figure FDA00041719512900000813
为第m次迭代过程中第k次工作过程噪声变量uc的方差，
Figure FDA00041719512900000814
为第m次迭代过程中第k次工作过程噪声变量ud的方差，
Figure FDA00041719512900000815
为第m次迭代过程中第k次工作测量噪声变量的方差；in,
Figure FDA00041719512900000811
is the variance of the noise variable u a in the kth working process during the mth iteration,
Figure FDA00041719512900000812
is the variance of the noise variable u b in the kth working process during the mth iteration,
Figure FDA00041719512900000813
is the variance of the noise variable u c in the kth working process during the mth iteration,
Figure FDA00041719512900000814
is the variance of the noise variable ud in the kth working process during the mth iteration,
Figure FDA00041719512900000815
The variance of the noise variable measured for the kth job during the mth iteration; 对于第m+1次迭代过程，可以分为E步和M步，可以表示成公式(29)—(30)所示：For the m+1th iteration process, it can be divided into the E step and the M step, which can be expressed as shown in formulas (29)-(30): (1)E步骤：计算(1) Step E: Calculation
Figure FDA0004171951290000091
Figure FDA0004171951290000091
 (2)M步骤：计算(2) Step M: Calculation
Figure FDA0004171951290000092
Figure FDA0004171951290000092
 其中，
Figure FDA0004171951290000093
为Ξ在
Figure FDA0004171951290000094
条件下的联合对数似然函数，
Figure FDA0004171951290000095
为r1:k,y1:k在
Figure FDA0004171951290000096
条件下的期望运算；
Figure FDA0004171951290000097
为第m+1次迭代过程对第k个工作过程的未知参数估计向量；in,
Figure FDA0004171951290000093
For
Figure FDA0004171951290000094
The joint log-likelihood function under the condition,
Figure FDA0004171951290000095
For r 1:k ,y 1:k
Figure FDA0004171951290000096
Expected operation under conditions;
Figure FDA0004171951290000097
is the unknown parameter estimation vector of the kth working process in the m+1th iteration process; 根据公式(28)，公式(29)可以扩展地写成公式(31)所示的形式：According to formula (28), formula (29) can be expanded to the form shown in formula (31):
Figure FDA0004171951290000098
Figure FDA0004171951290000098
 在公式(31)中，
Figure FDA0004171951290000099
Figure FDA00041719512900000910
和
Figure FDA00041719512900000911
是第i个粒子基于自初始到第k个工作循环的电池容量数据y1:k＝[y1,y2,…,yk]的隐含状态变量的条件期望；In formula (31),
Figure FDA0004171951290000099
Figure FDA00041719512900000910
and
Figure FDA00041719512900000911
is the conditional expectation of the implicit state variables of the ith particle based on the battery capacity data y 1:k = [y 1 , y 2 , …, y k ] from the initial to the kth working cycle; 基于UPF算法的前向迭代，可以获得粒子所代表的状态变量的期望和协方差的估计如公式(32)和公式(33)所示：Based on the forward iteration of the UPF algorithm, the expectation and covariance estimates of the state variables represented by the particles can be obtained as shown in formulas (32) and (33):
Figure FDA0004171951290000101
Figure FDA0004171951290000101
Figure FDA0004171951290000102
Figure FDA0004171951290000102
 其中，
Figure FDA0004171951290000103
为第i个粒子的状态变量在RTS算法中的期望，
Figure FDA0004171951290000104
为第i个粒子的状态变量在RTS算法中的协方差，
Figure FDA0004171951290000105
为无迹粒子滤波算法第i个粒子于第k次工作过程的后验状态变量期望，
Figure FDA0004171951290000106
为无迹粒子滤波算法第i个粒子于第k次工作过程的后验协方差矩阵；in,
Figure FDA0004171951290000103
is the expectation of the state variable of the ith particle in the RTS algorithm,
Figure FDA0004171951290000104
is the covariance of the state variable of the ith particle in the RTS algorithm,
Figure FDA0004171951290000105
is the posterior state variable expectation of the i-th particle in the k-th working process of the unscented particle filter algorithm,
Figure FDA0004171951290000106
is the posterior covariance matrix of the i-th particle in the k-th working process of the unscented particle filter algorithm; 相应地，k-1和k工作循环之间的状态变量之间的协方差可以表示成公式(34)所示的形式：Accordingly, the covariance between the state variables between k-1 and k working cycles can be expressed as shown in formula (34):
Figure FDA0004171951290000107
Figure FDA0004171951290000107
 其中，
Figure FDA0004171951290000108
为第i个粒子所代表的k-1和k工作循环之间的状态变量之间的协方差，
Figure FDA0004171951290000109
为第i个粒子的UKF增益，
Figure FDA00041719512900001010
为k-1和k工作循环之间的状态变量之间的协方差，I4×4为四维单位阵，
Figure FDA00041719512900001011
为无迹粒子滤波算法中第i个粒子所代表的k-1次工作过程的状态变量之间的后验协方差矩阵；in,
Figure FDA0004171951290000108
is the covariance between the state variables between the k-1 and k working cycles represented by the ith particle,
Figure FDA0004171951290000109
is the UKF gain of the ith particle,
Figure FDA00041719512900001010
is the covariance between the state variables between k-1 and k working cycles, I 4×4 is the four-dimensional unit matrix,
Figure FDA00041719512900001011
is the posterior covariance matrix between the state variables of the k-1 working processes represented by the i-th particle in the unscented particle filter algorithm; 根据公式(32)—(34)中的初始化条件，即可以进行RTS平滑运算；RTS平滑增益的表达式如公式(35)所示：According to the initialization conditions in formulas (32)-(34), RTS smoothing operation can be performed; the expression of RTS smoothing gain is shown in formula (35):
Figure FDA00041719512900001012
Figure FDA00041719512900001012
 其中，
Figure FDA00041719512900001013
为第i个粒子在第t个工作过程的RTS平滑增益，
Figure FDA00041719512900001014
为第i个粒子在第t个工作过程的UKF后验协方差，
Figure FDA00041719512900001015
为第i个粒子在第t个工作过程的UKF前向一步预测协方差；in,
Figure FDA00041719512900001013
is the RTS smoothing gain of the ith particle in the tth working process,
Figure FDA00041719512900001014
is the UKF posterior covariance of the ith particle in the tth working process,
Figure FDA00041719512900001015
is the UKF forward one-step prediction covariance of the ith particle in the tth working process; 相应地，后向迭代过程中的状态变量
Figure FDA0004171951290000111
和协方差
Figure FDA0004171951290000112
的更新如公式(36)—(37)所示：Accordingly, the state variables in the backward iteration process are
Figure FDA0004171951290000111
and covariance
Figure FDA0004171951290000112
The update of is shown in formulas (36)-(37):
Figure FDA0004171951290000113
Figure FDA0004171951290000113
Figure FDA0004171951290000114
Figure FDA0004171951290000114
 其中，
Figure FDA0004171951290000115
为第i个粒子在RTS后向平滑中对于第t+1个工作过程的后验状态向量，
Figure FDA0004171951290000116
为第i个粒子在UKF前向预测中对于第t+1个工作过程状态向量的预测值，
Figure FDA0004171951290000117
为第i个粒子在RTS后向平滑中对于第t个工作过程的先验状态向量，
Figure FDA0004171951290000118
为第i个粒子在RTS后向平滑中对于第t个工作过程的后验状态向量，
Figure FDA0004171951290000119
为第i个粒子在RTS后向平滑中对于第t+1个工作过程的协方差，
Figure FDA00041719512900001110
为第i个粒子为UKF前向预测中对于第t+1个工作过程协方差的预测值，
Figure FDA00041719512900001111
为第i个粒子在UKF前向运算中对于第t个工作过程的后验协方差，
Figure FDA00041719512900001112
为第i个粒子在RTS后向平滑中对于第t个工作过程的协方差；in,
Figure FDA0004171951290000115
is the posterior state vector of the ith particle in the RTS backward smoothing for the t+1th working process,
Figure FDA0004171951290000116
is the predicted value of the state vector of the t+1th working process of the ith particle in the UKF forward prediction,
Figure FDA0004171951290000117
is the prior state vector of the ith particle in the RTS backward smoothing for the tth working process,
Figure FDA0004171951290000118
is the posterior state vector of the ith particle in RTS backward smoothing for the tth working process,
Figure FDA0004171951290000119
is the covariance of the ith particle in the RTS backward smoothing for the t+1th working process,
Figure FDA00041719512900001110
is the predicted value of the covariance of the t+1th working process in the UKF forward prediction for the ith particle,
Figure FDA00041719512900001111
is the posterior covariance of the ith particle in the UKF forward operation for the tth working process,
Figure FDA00041719512900001112
is the covariance of the ith particle in the RTS backward smoothing for the tth working process; 除此之外，相邻两个循环之间状态矩阵之间的协方差可以表示成公式(38)所示的形式：In addition, the covariance between the state matrices between two adjacent cycles can be expressed as shown in formula (38):
Figure FDA00041719512900001113
Figure FDA00041719512900001113
 其中，
Figure FDA00041719512900001114
为RTS后向平滑中第t-1个工作过程和第t个工作过程状态变量的协方差，
Figure FDA00041719512900001115
为RTS后向平滑中第t个工作过程和第t+1个工作过程状态变量的协方差；in,
Figure FDA00041719512900001114
is the covariance of the state variables of the t-1th working process and the tth working process in RTS backward smoothing,
Figure FDA00041719512900001115
is the covariance of the state variables of the t-th working process and the t+1-th working process in RTS backward smoothing; 基于公式(32)—(38)，可以求解条件期望表达式
Figure FDA00041719512900001116
Figure FDA00041719512900001117
和
Figure FDA00041719512900001118
如式(39)—(42)所示：Based on formulas (32)-(38), the conditional expectation expression can be solved
Figure FDA00041719512900001116
Figure FDA00041719512900001117
and
Figure FDA00041719512900001118
As shown in formulas (39)-(42):
Figure FDA0004171951290000121
Figure FDA0004171951290000121
Figure FDA0004171951290000122
Figure FDA0004171951290000122
Figure FDA0004171951290000123
Figure FDA0004171951290000123
Figure FDA0004171951290000124
Figure FDA0004171951290000124
 其中，
Figure FDA0004171951290000125
为RTS后向平滑中第t-1个工作过程第i个粒子所代表的状态变量；in,
Figure FDA0004171951290000125
is the state variable represented by the i-th particle in the t-1-th working process in RTS backward smoothing; 将公式(39)—(42)中的条件期望代入到公式(31)中，可以写成公式(43)所示的形式Substituting the conditional expectations in formulas (39)-(42) into formula (31), it can be written as shown in formula (43):
Figure FDA0004171951290000126
Figure FDA0004171951290000126
 基于公式(31)—(43)中E步骤的计算结果，根据公式(30)计算EM算法的M步骤，针对退化模型的未知参数估计向量
Figure FDA0004171951290000127
中每一个参数的偏导数，令其值为0，可以求解出公式(44)—(45)的结果；Based on the calculation results of step E in formulas (31)-(43), the M step of the EM algorithm is calculated according to formula (30), and the unknown parameter estimation vector of the degradation model is
Figure FDA0004171951290000127
Taking the partial derivative of each parameter in and setting its value to 0, we can solve the results of formulas (44)-(45);
Figure FDA0004171951290000131
Figure FDA0004171951290000131
Figure FDA0004171951290000132
Figure FDA0004171951290000132
 通过不断进行E步骤和M步骤的迭代，直到满足
Figure FDA0004171951290000133
或者达到最大迭代次数10次为止停止迭代，由此实现关于退化模型中未知参数的自适应估计；By continuously iterating the E and M steps until the
Figure FDA0004171951290000133
Or the iteration is stopped when the maximum number of iterations reaches 10, thereby realizing adaptive estimation of unknown parameters in the degradation model; 其中，ε为人为设定的参数；Among them, ε is a parameter set artificially; 所述步骤4中判断第k次工作过程是否发生容量再生现象；具体过程为：In step 4, it is determined whether capacity regeneration occurs during the kth working process; the specific process is as follows: 在经过UKF算法更新后的电池容量的集合为
Figure FDA0004171951290000134
在经过UPF算法更新后的电池容量的集合为
Figure FDA0004171951290000135
The set of battery capacities after being updated by the UKF algorithm is
Figure FDA0004171951290000134
The set of battery capacities after being updated by the UPF algorithm is
Figure FDA0004171951290000135
 假设H0：两个总体
Figure FDA0004171951290000136
和
Figure FDA0004171951290000137
在显著性水平α下具备相同的概率密度分布；Hypothesis H 0 : Two populations
Figure FDA0004171951290000136
and
Figure FDA0004171951290000137
Have the same probability density distribution at the significance level α; 为了检验假设H0，Wilcoxon秩和检验可以总结为以下的步骤：To test the hypothesis H 0 , the Wilcoxon rank sum test can be summarized as follows: (1)将两个总体
Figure FDA0004171951290000138
和
Figure FDA0004171951290000139
混合获得总体Qk，依照电池容量的数值从小到大进行统一编号，定义每个电池容量样本对应的序数为秩；(1) The two populations
Figure FDA0004171951290000138
and
Figure FDA0004171951290000139
The population Q k is obtained by mixing, and the battery capacity values are uniformly numbered from small to large, and the ordinal number corresponding to each battery capacity sample is defined as the rank; (2)计算样本总体
Figure FDA00041719512900001310
的样本所对应的秩之和T；(2) Calculate the sample population
Figure FDA00041719512900001310
The sum of the ranks corresponding to the samples of T; (3)根据粒子总数Ns和显著性水平γ查阅秩和检验表，可以获得秩和下限T1、秩和上限T2以及相应的置信度水平p；(3) According to the total number of particles Ns and the significance level γ, the rank sum test table can be consulted to obtain the rank sum lower limit T1 , the rank sum upper limit T2 and the corresponding confidence level p; (4)设置区间[T1,T2]，如果T在区间内部，也就是p>γ，则接受假设H0，认为两个总体
Figure FDA00041719512900001311
和
Figure FDA00041719512900001312
在水平α下无显著性差异；否则拒绝假设H0，认为两个总体
Figure FDA00041719512900001313
和
Figure FDA00041719512900001314
在水平α下存在显著性差异；(4) Set the interval [T 1 ,T 2 ]. If T is within the interval, that is, p>γ, then accept the hypothesis H 0 and assume that the two populations are
Figure FDA00041719512900001311
and
Figure FDA00041719512900001312
There is no significant difference at level α; otherwise, the hypothesis H 0 is rejected and the two populations are considered
Figure FDA00041719512900001313
and
Figure FDA00041719512900001314
There is a significant difference at level α; 也就意味着在第k个工作循环中电池产生了容量再生的现象。This means that the battery has capacity regeneration in the kth working cycle. 2.根据权利要求1所述基于最大期望-无迹粒子滤波的电池剩余使用寿命预测方法，其特征在于：所述步骤5中根据步骤2中构建的电池退化模型和当前的电池容量值，通过趋势外推求解电池剩余使用寿命；2. The method for predicting the remaining useful life of a battery based on maximum expectation-unscented particle filtering according to claim 1, characterized in that: in step 5, the remaining useful life of the battery is solved by trend extrapolation according to the battery degradation model constructed in step 2 and the current battery capacity value; 根据步骤2中构建的电池退化模型和当前的电池容量值，计算电池容量的置信区间；Calculate the confidence interval of the battery capacity based on the battery degradation model constructed in step 2 and the current battery capacity value; 具体过程为：The specific process is: 1)、令预测的起始点v＝k，将公式(23)状态变量期望带入公式
Figure FDA0004171951290000141
计算电池容量；1) Let the starting point of the prediction v = k, and substitute the state variable expectation of formula (23) into the formula
Figure FDA0004171951290000141
Calculate battery capacity; 若检测到第k个工作循环中电池产生了容量再生的现象，则需将预测的起始点v＝k替代为v＝k-1；If it is detected that the battery has regenerated capacity in the kth working cycle, the predicted starting point v=k needs to be replaced by v=k-1; 2)、如果计算出的电池容量yv小于等于失效阈值，则执行3)；2) If the calculated battery capacity y v is less than or equal to the failure threshold, execute 3); 如果计算出的电池容量yv大于失效阈值，则令v＝v+1，通过公式
Figure FDA0004171951290000142
计算电池容量，重复执行2)，直到解算出来的yv小于等于电池容量的失效阈值；If the calculated battery capacity y v is greater than the failure threshold, let v = v + 1, and use the formula
Figure FDA0004171951290000142
Calculate the battery capacity and repeat step 2) until the calculated y v is less than or equal to the failure threshold of the battery capacity; 3)、通过公式(46)获得电池剩余使用寿命：3) Obtain the remaining battery life through formula (46): RULk＝v–k (46)RUL k = v–k (46) 4)、如果计算出的电池容量yv大于失效阈值，获取Ns个粒子的均值和标准差，均值加减正负三倍标准差作为电池容量置信区间上下线。4) If the calculated battery capacity y v is greater than the failure threshold, obtain the mean and standard deviation of N s particles, and use the mean plus or minus three times the standard deviation as the upper and lower limits of the battery capacity confidence interval.
CN202210415333.5A 2022-04-20 2022-04-20 Battery remaining service life prediction method based on maximum expected-unscented particle filtering Active CN114779088B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210415333.5A CN114779088B (en) 2022-04-20 2022-04-20 Battery remaining service life prediction method based on maximum expected-unscented particle filtering

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210415333.5A CN114779088B (en) 2022-04-20 2022-04-20 Battery remaining service life prediction method based on maximum expected-unscented particle filtering

Publications (2)

Publication Number Publication Date
CN114779088A CN114779088A (en) 2022-07-22
CN114779088B true CN114779088B (en) 2023-05-12

Family

ID=82431009

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210415333.5A Active CN114779088B (en) 2022-04-20 2022-04-20 Battery remaining service life prediction method based on maximum expected-unscented particle filtering

Country Status (1)

Country Link
CN (1) CN114779088B (en)

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150349385A1 (en) * 2014-04-01 2015-12-03 Medtronic, Inc. Method and System for Predicting Useful Life of a Rechargeable Battery
CN105445671A (en) * 2015-12-29 2016-03-30 北京航天测控技术有限公司 Lithium ion battery service life prediction method based on traceless particle filtering
CN106055775B (en) * 2016-05-27 2018-12-11 哈尔滨工业大学 A kind of service life of secondary cell prediction technique that particle filter is combined with mechanism model
CN109932656A (en) * 2019-04-17 2019-06-25 合肥工业大学 A Lithium Battery Life Estimation Method Based on IMM-UPF
CN112415414A (en) * 2020-10-09 2021-02-26 杭州电子科技大学 Method for predicting remaining service life of lithium ion battery
CN112487702B (en) * 2020-10-26 2024-04-05 湖州师范学院 Method for predicting residual service life of lithium ion battery
CN112949026B (en) * 2021-01-19 2023-05-23 中国人民解放军火箭军工程大学 Age and state dependence considered degradation equipment residual life prediction method
CN112800616B (en) * 2021-02-05 2023-07-18 中国人民解放军空军工程大学 Adaptive Prediction Method of Remaining Life of Equipment Based on Proportional Accelerated Degradation Modeling
CN113093020B (en) * 2021-04-02 2022-07-12 中国矿业大学 A method for predicting the remaining service life of lithium-ion batteries based on LSTM neural network
CN114252797B (en) * 2021-12-17 2023-03-10 华中科技大学 A method for predicting the remaining service life of lithium batteries based on uncertainty estimation

Also Published As

Publication number Publication date
CN114779088A (en) 2022-07-22

Similar Documents

Publication Publication Date Title
CN111443294B (en) Method and device for indirectly predicting remaining life of lithium ion battery
CN113030744B (en) Battery health condition prediction method, system and medium based on health factor extraction
CN107957562B (en) An online prediction method for the remaining life of lithium-ion batteries
CN104505894B (en) Power management system and state estimation method based on mining lithium ion batteries
CN110457789B (en) Lithium ion battery residual life prediction method
CN113917337A (en) Battery state of health estimation method based on charging data and LSTM neural network
CN112782591B (en) Lithium battery SOH long-term prediction method based on multi-battery data fusion
CN113253116A (en) Lithium ion battery state of charge estimation method and storage medium
CN109061506A (en) Lithium-ion-power cell SOC estimation method based on Neural Network Optimization EKF
CN114325450A (en) Lithium ion battery health state prediction method based on CNN-BilSTM-AT hybrid model
CN110058160B (en) Lithium battery state of health prediction method based on square root extended Kalman filter
CN108872866A (en) A kind of charge states of lithium ion battery dynamic evaluation and long-acting prediction fusion method
CN114636932A (en) Method and system for predicting remaining service life of battery
CN111220920B (en) State of charge calculation method for decommissioned lithium-ion batteries based on H∞ unscented Kalman filter algorithm
CN114397577A (en) New energy automobile lithium battery health state assessment method based on ASTUKF-GRA-LSTM model
CN114839538A (en) Method for extracting degradation characteristics of lithium ion battery for estimating residual life
CN114487890A (en) An improved long short-term memory neural network for lithium battery state of health estimation method
CN114660497A (en) Lithium ion battery service life prediction method aiming at capacity regeneration phenomenon
Chang et al. Application of radial basis function neural network, to estimate the state of health for LFP battery
CN112014757A (en) Battery SOH estimation method integrating capacity increment analysis and genetic wavelet neural network
Yang et al. Lithium-ion battery life cycle prediction with deep learning regression model
CN115542182A (en) SOH estimation method for series battery pack of mobile energy storage system
CN114779088B (en) Battery remaining service life prediction method based on maximum expected-unscented particle filtering
CN114415041A (en) Battery remaining capacity estimation method, device and terminal equipment
CN117826000A (en) A method for estimating the health status of lithium-ion batteries

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant