CN105093121A - 似然度函数粒子滤波的动力电池电荷状态估计方法和系统 - Google Patents

Abstract

本发明为一种似然度函数粒子滤波的动力电池电荷状态估计方法和系统，由电池Thevenin模型，得到状态和量测方程，参数初始化后，进行状态预测计算状态预测值的均值和协方差，再重新采样、重构采样分布函数。计算电池端电压预测值，计算粒子权值，权值归一化和计算有效粒子数。有效粒子数N_eff与有效粒子数阈值N_thr比较，当N_eff小于N_thr，采用拉普拉斯分布作为似然度函数，并引入方差调节因子和工况适应因子，以自适应地修改似然度函数的方差，适应动力电池不同工况。最终得到更新的SOC估计值和协方差。本系统微控制器连接电压和电流传感器，微控制器内有各程序执行模块。本发明加大有效粒子数；有效避免了方差的过修正；估算精度优。

Description

似然度函数粒子滤波的动力电池电荷状态估计方法和系统

技术领域

本发明涉及电动汽车动力电池电荷状态估计技术领域，具体为利用具有方差调整因子的拉普拉斯分布似然度函数粒子滤波的动力电池电荷状态估计方法和系统。

背景技术

近年来，由于锂离子电池的优异性能，使其被广泛应用于电动汽车、便携式电子设备及航空航天等领域。动力电池电荷状态(StateofCharge，SOC)是电池管理系统(BatteryManagementSystem，BMS)描述电池状态的一个重要参数，对电荷状态准确的估算，可防止电池的过充和过放，有效延长电池的使用寿命。

目前研究者已给出了诸多实时估算电荷状态SOC的方法。

安时法通过精确测量电流实现电荷状态估计，是电荷状态估计的最基本方法，也是目前SOC测量的标准方法。

基于电压的SOC估计方法是通过SOC-OCV(opencircuitvoltage,开路电压)关系表,查表获得SOC，但OCV的测量需要电池长时间的静置。

卡尔曼滤波(KF，KalmanFilter)也是电荷状态估计的常用方法。考虑到电池的非线性特性，此SOC估计方法又进一步改进为有扩展卡尔曼滤波(EKF，ExtentedKalmanFilter)，无迹卡尔曼滤波(UKF，UnscentedKalmanFilter)，自适应无迹卡尔曼滤波，粒子滤波(PF，ParticleFilter)，无迹粒子滤波(UPF，UnscentedParticleFilter)等。

高斯粒子滤波(GPF)算法是粒子滤波(PF)算法的一种变形，高斯滤波和粒子滤波结合，称为高斯粒子滤波(GaussianParticleFilter，GPF)。它通过高斯分布来近似未知变量的后验分布,实时性上要优于粒子滤波算法,性能上要优于扩展卡尔曼滤波(EKF)、无味卡尔曼滤波(UKF)等算法。高斯粒子滤波比其它的高斯滤波方法适用性更强，能处理更多非线性动态系统问题；而与一般的粒子滤波相比，因为高斯粒子滤波用高斯分布近似后验分布，所以只要所用的高斯分布是正确的，就不会产生粒子退化问题，就不需要对粒子进行重采样，从而使算法的计算量降低，复杂度也降低。

高斯粒子滤波以预测概率密度函数N(x_k-1|k-1,P_k|k-1)，作为重要性概率密度函数，但没有考虑最新的量测信息，重要性采样粒子分布与后验概率分布产生的粒子存在较大偏差，归一化重要性权值只集中在部分粒子上，导致其他粒子权重几乎接近于0，使得权值方差增大，故有效粒子数减小，滤波性能下降。

发明内容

本发明的目的是设计一种似然度函数粒子滤波的动力电池电荷状态估计方法，以端电压预测值与量测值之差作为似然度函数的输入，在有效粒子数较少时,用拉普拉斯分布的似然度函数替换高斯分布的似然度函数，更新计算权值，拉普拉斯分布的似然度函数在峰值之后逐渐减小，下降速度慢，似然度函数曲线与横轴之间的高度降低缓慢，且可自适应调节，增加有效粒子数，拓宽了下一时刻抽样范围，提高了滤波稳定性；引入方差调节因子γ，自适应地修改似然度函数的方差，有效地避免了方差的过修正问题。

本发明的另一目的是设计一种实现上述似然度函数粒子滤波的动力电池电荷状态估计方法的似然度函数粒子滤波的动力电池电荷状态估计系统。

本发明设计的一种似然度函数粒子滤波的动力电池电荷状态估计方法，采用目前最广泛使用的一种电池等效模型Thevenin模型，描述电池的静态和动态性能。电池的极化电阻R_p与电池的极化电容C_p并联构成一阶RC结构，表示电池的极化反应，RC两端电压为U_p(t)；串接欧姆电阻R₀和Uoc，Uoc为电池的开路电压OCV，采样得到电池端电压U(t)和流过欧姆内阻R₀的电流i(t)。

电池Thevenin模型数学表达式如下：

$\{\begin{array}{c}\frac{{\mathrm{du}}_p\left(t\right)}{dt}=-\frac{U_p\left(t\right)}{R_p{C}_p}+\frac{i\left(t\right)}{C_p}\\ U\left(t\right)=U_{OC}\left(t\right)-R_{0}i\left(t\right)-U_p\left(t\right)\end{array}---\left(1\right)$

用后向差分法对Thevenin电池模型离散化，整理后得到对应的后向差分离散模型：

U(k)＝aU(k-1)+bI(k)+cI(k-1)+(1-a)U_OC(k)(2)

式中：U_OC(k)表示k时刻的开路电压；U(k)为k时刻的电池输出端电压；I(k)为k时刻的回路电流；a,b,c为离散模型参数。

采用含遗忘因子的递推最小二乘法(FFRLS，Forgetting-factorrecursiveleastsquaremethod)，进行电池后向差分模型的参数辨识，求得a、b、c、U_OC(k)的值。根据方程(2)，对应得到k时刻的电池端电压的预测值U(k)。

定义SOC如下：

$SOC\left(t\right)={\mathrm{SOC}}_0-1/Q_N{\∫}_0^t{\mathrm{\η}}_{0}i\left(t\right)dt---\left(3\right)$

式中：SOC₀是SOC的初始值；Q_N为电池的额定容量；η₀为电池的充放电库伦效率，本发明η₀＝1；i(t)为放电电流。

结合电池的离散模型，对电池特性描述如下，

状态方程 ${\mathrm{SOC}}_k={\mathrm{SOC}}_{k-1}-\frac{i_{k-1}T}{Q_N}---\left(4\right)$

量测方程U_K＝U_OC(SOC_k)-R₀i_k-U_p,k(5)

式中：U_p,k为k时刻极化电容两端电压；T为采样周期；U_OC(SOC_k)表示电池的开路电压和SOC的非线性关系。

U_OC(SOC_k)与SOC关系如下：

U_oc(SOC_k)＝a₁SOC_k ⁸+a₂SOC_k ⁷+a₃SOC_k ⁶+a₄SOC_k ⁵+

(6)

a₅SOC_k ⁴+a₆SOC_k ³+a₇SOC_k ²+a₈SOC_k+a₉

根据实验得到的SOC和在线辨识得到的开路电压U_OC，运用含遗忘因子的递推最小二乘法(FFRLS)，求出式(6)中的a₁～a₉。

根据式(5)，由SOC的预测值SOC_k得到端电压的预测值U_k。

本方法主要包括如下步骤：

步骤Ⅰ、参数初始化

设定SOC的初始值和滤波协方差初始值，分别记作x_0|0和P_0|0，设定采样粒子数为N，N取值为30～100，较佳方案N为45～55，有效粒子数阈值N_thr。

根据经验选定N_thr＝(0.5～0.85)N，较佳方案为N_thr＝(0.73～0.81)N。

步骤Ⅱ、重要性采样

从正态分布，即高斯分布N(x_k-1|k-1,P_k-1|k-1)中随机采样N个粒子记作为从N(x_k-1|k-1,P_k-1|k-1)中采样得到的大小为N的以x_k-1|k-1为中心值、P_k-1|k-1为方差的高斯分布；x_k-1|k-1和P_k-1|k-1分别为k-1时刻滤波均值和协方差，即SOC_k-1及其协方差。

Ⅱ-1、状态预测

采用步骤Ⅰ中的式(4)

状态方程 ${\mathrm{SOC}}_k={\mathrm{SOC}}_{k-1}-\frac{i_{k-1}T}{Q_N}$

进行状态预测，即求得：

其中，为k-1时刻的对SOC的第i个采样点；为k时刻从预测粒子集采样的第i个采样点；i_k-1为k-1时刻的电流，放电时电流值为正，充电时电流值为负，；T为采样周期，Q_N为电池总容量。

Ⅱ-2、计算状态预测值的均值x_k|k-1和协方差P_k|k-1

$x_{k|k-1}=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}{\hat{x}}_k^i---\left(7\right)$

$P_{k|k-1}=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}({\hat{x}}_k^i-x_{k|k-1}){({\hat{x}}_k^i-x_{k|k-1})}^T---\left(8\right)$

Ⅱ-3、重新采样、重构采样分布函数

以状态预测值的均值x_k-1|k-1和协方差P_k-1|k-1重构高斯分布得N(x_k|k-1,P_k|k-1)，并重新采样为从高斯分布N(x_k|k-1,P_k|k-1)中采样的N个粒子。

步骤Ⅲ、计算粒子权值并归一化

Ⅲ-1、计算电池端电压预测值

将粒子代入量测方程(5)中计算电池端电压预测值 表示k时刻粒子对应的端电压预测值集合。

即： $y_{1,k}^i=U_{oc}\left({\tilde{x}}_k^i\right)-R_0{i}_k-U_{p,k},$

由式(6)计算得到。

Ⅲ-2、计算粒子权值

k时刻第i个粒子的权值由高斯分布的似然度函数计算得到，即 ${\mathrm{\ω}}_k^i=1/\sqrt{2\mathrm{\π}R}\exp (-1/(2R\left){\left(y_{1,k}-y_{1,k}^i\right)}^2\right),$ R为高斯似然度函数的方差，根据经验赋值。y_1,k为k时刻的量测值，表示k时刻第i个粒子对应的端电压预测值。

Ⅲ-3、权值归一化和计算有效粒子数N_eff

${\mathrm{\ω}}_k^i=\frac{{\widehat{\mathrm{\ω}}}_k^i}{\underset{i=1}{\overset{N}{\Σ}}{\widehat{\mathrm{\ω}}}_k^i}---\left(9\right)$

$N_{eff}=\frac{1}{\underset{i=1}{\overset{N}{\Σ}}{\left({\mathrm{\ω}}_k^i\right)}^2}---\left(10\right)$

为归一化后的权值；N_eff为相对于采样粒子数N而言的有效粒子数。 $N_{eff}=N/(1+Var({\mathrm{\ω}}_k^i\left)\right),$ 为的方差。

$x_{k|k}=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_k^i{\tilde{x}}_k^i---\left(11\right)$

$P_{k|k}=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_k^i({\tilde{x}}_k^i-x_{k|k}){({\tilde{x}}_k^i-x_{k|k})}^T---\left(12\right)$

x_k|k为状态估计值；P_k|k为协方差估计值。

步骤Ⅳ、有效粒子数N_eff与有效粒子数阈值N_thr的比较

高斯粒子滤波的预测粒子分布即先验概率分布π(x_k|k-1)，代表粒子数分布；后验概率分布即粒子实际概率分布，代表粒子的权值分布；先验概率分布和后验概率分布存在较大偏差。当高斯后验概率分布与预测粒子分布的重叠越少，预测粒子分布中对应的有效粒子数越少。权值均集中于部分粒子，有效粒子数少，由上述式(12)可知协方差P_k|k-1小，下一时刻按照N(x_k|k-1,P_k|k-1)采样的范围小，容易丢失目标，滤波性能差。

主要表现有两种情况：

1)高斯后验概率分布与预测粒子分布中心点重合，后验概率分布曲线峰值A大，而预测粒子分布曲线虽然存在峰值B，但相比较来讲尖峰平坦，峰值B较小；A大于B的三倍或更多；后验概率分布峰值两侧下降越快，高斯后验概率分布与预测粒子分布重叠越少，使预测粒子分布中对应的有效粒子数愈少。

2)高斯后验概率分布与预测粒子分布中心点不重合，后验概率分布曲线位于预测粒子分布曲线的尾部。当高斯后验概率分布与预测粒子分布中心点相距越大，高斯后验概率分布与预测粒子分布重叠越少，使预测粒子分布中对应的有效粒子数愈少。

拉普拉斯分布与高斯后验概率分布比较，峰值相同情况的两种分布曲线，在峰值两侧幅值均逐渐降低。将曲线与横轴间的高度定义为拖尾厚度。在接近横轴时，拉普拉斯分布的幅值下降比高斯后验概率分布的幅值下降更缓慢，即拉普拉斯分布的拖尾厚度大于高斯后验概率分布的拖尾厚度。

与高斯后验概率分布峰值相同的拉普拉斯分布，无论是拉普拉斯分布与预测粒子分布最大峰值中心点重合或不重合，拉普拉斯分布与预测粒子分布的重叠都比高斯后验概率分布多，即，预测粒子分布中有效粒子数增大。增大了滤波的协方差P_k|k-1，下一时刻按照N(x_k|k-1,P_k|k-1)进行抽样，采样范围更大，增强了滤波的稳定性。

因此根据步骤Ⅰ设定的有效粒子数阈值N_thr，确定是否更换似然度函数。

当N_eff小于N_thr，进入步骤Ⅴ；否则进入步骤Ⅵ。

步骤Ⅴ、更换似然度函数

采用拉普拉斯分布作为似然度函数。

拉普拉斯分布作为似然度函数增加了与预测粒子分布的重叠部分的面积，但对应增加了离量测值较远的粒子的权值。而离量测值较远的粒子的权值越大，使得当前时刻滤波的值偏差越大。

本发明引入方差调节因子γ，γ＝1+θ；有效粒子数退化率θ，θ＝(N_thr-N_eff)/N_thr。N_thr为有效粒子数阈值；N_eff表示有效粒子数，由公式(10)计算得到。

拉普拉斯分布的概率密度函数，或称为似然度函数如下：

$f\left(x\right|\mathrm{\μ},R)=\frac{1}{2R}\exp (-\frac{|x-\mathrm{\μ}|}{R})---\left(13\right)$

式中：x为预测方程所得状态预测值，或者x为量测方程所得量测值对应物理量的预测值，μ和R分别是该拉普拉斯分布的中心值和方差。

引入γ后，方程(13)改写为：

$\begin{array}{c}f\left(x\right|\mathrm{\μ},R\mathrm{\δ})=\frac{1}{2\left(R\mathrm{\γ}\right)}\exp (-\frac{|x-\mathrm{\μ}|}{R\mathrm{\γ}})\\ =\frac{1}{2\left(R\right(1+\left(N_{thr}-N_{eff}\right)/N_{thr}\left)\right)}\exp (-\frac{|x-\mathrm{\μ}|}{R(1+(N_{thr}-N_{eff})/N_{thr})})\end{array}---\left(14\right)$

当γ增大时拉普拉斯分布的拖尾厚度增加。根据γ和θ的计算公式，当N_eff越小，θ越大，γ越大，拖尾厚度越大，从而实现自适应地修改似然度函数的方差，有效地避免了方差的过修正问题。

动力电池在不同工况下，SOC变化快慢不一样，预测粒子分布和似然度函数的相对位置及形状不一样，为此，在γ的定义式中引入工况适应因子ε，γ＝1+εθ，ε～(0,1]。

代入式(14)得：

$\begin{array}{c}f\left(x\right|\mathrm{\μ},R\mathrm{\γ})=f\left(x\right|\mathrm{\μ},R\mathrm{\ϵ})\\ =\frac{1}{2\left(R\right(1+\mathrm{\ϵ}\left(N_{thr}-N_{eff}\right)/N_{thr}\left)\right)}\exp \left(-\frac{|x-\mathrm{\μ}|}{R(1+\mathrm{\ϵ}(N_{thr}-N_{eff})/N_{thr})}\right)\end{array}$

动力电池工况负载变化不剧烈时，SOC变化相对较慢，即前后两个时刻的SOC真实值相近，以N(x_k|k-1,P_k|k-1)为重要性概率密度函数的预测概率密度函数的中心，与似然度函数中心较为接近，出现似然度函数位于状态转移概率密度函数尾部的概率小，似然度函数的拖尾厚度越小ε取较小值。

动力电池工况剧烈时，SOC相对变化较快，即前后两个时刻的SOC真实值相差较大，预测概率密度函数的中心与似然度函数中心相隔较远的可能性大，出现似然度函数位于状态转移概率密度函数尾部的概率大，可使似然度函数的拖尾厚度大，ε取较大值。本发明根据DST工况，取ε＝0.4～0.6。

借鉴高斯滤波中的方法，用代替拉普拉斯分布的概率密度函数(13)中的x-μ，得到修改后的拉普拉斯分布似然度函数

${\mathrm{\ω}}_k^i=1/\left(2\mathrm{\δ}R\right)\exp (-|y_{1,k}-y_{1,k}^i|/(\mathrm{\δ}R\left)\right)$

其中x、μ、R的赋值分别参照高斯分布的似然度函数中的y_1,k、R，然后进行归一化。

步骤Ⅵ、状态更新和协方差更新

即执行步骤Ⅲ-3的式(11)(12)，得到更新的SOC估计值和更新的协方差。

步骤Ⅶ、判断滤波时间大于设定运行时间，则结束，否则，k＝k+1，返回步骤Ⅱ。

本发明设计的似然度函数粒子滤波的动力电池电荷状态估计系统包括微控制器、电压传感器和电流传感器，微控制器还连接显示器。微控制器配有通用接口，电压传感器和电流传感器接入模数转换电路，再经通用接口连接微控制器。

微控制器的中心处理器连接有程序存储器和数据存储器，程序存储器含有电压电流数据采集模块、根据参数辨识的U_OC(SOC_k)与SOC关系模块，端电压估计及端电压预测偏差计算模块，有效粒子数N_eff与有效粒子数阈值N_thr的比较模块，方差调节因子、有效粒子数退化率及修改后的拉普拉斯分布似然度函数计算模块，状态更新和协方差更新模块以及SOC显示模块；数据存储器存放初始化参数。中心处理器连接CAN总线接口、PC104总线接口和以太网接口中的一种或多种，以供功能扩展使用。

与现有技术相比，本发明似然度函数粒子滤波的动力电池电荷状态估计方法和系统的优点为：1、针对高斯粒子滤波的动力电池电荷状态估计方法的重要性概率密度函数未考虑最新的量测信息，导致有效粒子数减少，滤波性能下降的问题，在有效粒子数小于有效粒子数阈值时，以拉普拉斯分布作为似然度函数，加大有效粒子数；2、引入方差调节因子增加了离量测值较远的粒子的权值，实现自适应地修改似然度函数的方差，有效地避免了方差的过修正问题；3、在方差调节因子的定义式中引入工况适应因子，根据动力电池工况负载变化，调节预测粒子分布和似然度函数的相对位置及形状；4、本方法的估算精度要优于高斯粒子滤波的动力电池电荷状态估计方法，且当SOC初始值偏差高达80％时，本发明方法依然迅速收敛，保持较高精度。

附图说明

图1为本似然度函数粒子滤波的动力电池电荷状态估计方法实施例的流程图；

图2为本似然度函数粒子滤波的动力电池电荷状态估计方法实施例所采用的动力电池等效模型的原理图；

图3为本似然度函数粒子滤波的动力电池电荷状态估计方法实施例步骤Ⅳ中高斯后验概率分布、拉普拉斯后验概率分布和预测粒子分布中心点重合，且二前者相对于第三者呈尖峰状态；

图4为本似然度函数粒子滤波的动力电池电荷状态估计方法实施例步骤Ⅳ中高斯后验概率分布和拉普拉斯后验概率分布位于预测粒子分布的尾部；

图5为本似然度函数粒子滤波的动力电池电荷状态估计方法实施例步骤Ⅳ中不同的γ条件下的拉普拉斯分布的拖尾厚度对比；

图6为本似然度函数粒子滤波的动力电池电荷状态估计方法实施例(PF-LFAF)工况适应因子ε取0.5所得SOC与高斯粒子滤波(GPF)所得SOC结果对比图，本图所示为SOC初值SOC₀＝0.8、采样粒子数N＝50；

图7为本似然度函数粒子滤波的动力电池电荷状态估计方法实施例(PF-LFAF)工况适应因子ε取0.5所得SOC与高斯粒子滤波(GPF)所得SOC结果对比图，本图所示为SOC初值SOC₀＝0.2、采样粒子数N＝50；

图8为本似然度函数粒子滤波的动力电池电荷状态估计方法实施例(PF-LFAF)工况适应因子ε取0.5所得SOC与高斯粒子滤波(GPF)所得SOC结果对比图，本图所示为SOC初值SOC₀＝0.8、采样粒子数N＝100；

图9为本似然度函数粒子滤波的动力电池电荷状态估计方法实施例(PF-LFAF)工况适应因子ε取0.5所得SOC与高斯粒子滤波(GPF)所得SOC结果对比图，本图所示为SOC初值SOC₀＝0.2、采样粒子数N＝100。

图10为本似然度函数粒子滤波的动力电池电荷状态估计系统实施例结构框图。

具体实施方式

似然度函数粒子滤波的动力电池电荷状态估计方法实施例

本似然度函数粒子滤波的动力电池电荷状态估计方法实施例，其流程如图1所示。采用目前最广泛使用的一种电池等效模型Thevenin模型，描述电池的静态和动态性能。如图2所示，电池的极化电阻R_p与电池的极化电容C_p并联构成一阶RC结构，表示电池的极化反应，RC两端电压为U_p(t)；串接欧姆电阻R₀和Uoc，Uoc为电池的开路电压OCV，采样得到电池端电压U(t)和流过欧姆内阻R₀的电流i(t)。

电池Thevenin模型数学表达式如下：

$\{\begin{array}{c}\frac{{\mathrm{dU}}_p\left(t\right)}{dt}=-\frac{U_p\left(t\right)}{R_p{C}_p}+\frac{i\left(t\right)}{C_p}\\ U\left(t\right)=U_{OC}\left(t\right)-R_{0}i\left(t\right)-U_p\left(t\right)\end{array}---\left(1\right)$

用后向差分法对Thevenin电池模型离散化，整理后得到对应的后向差分离散模型：

U(k)＝aU(k-1)+bI(k)+cI(k-1)+(1-a)U_OC(k)(2)

式中：U_OC(k)表示k时刻的开路电压；U(k)为k时刻的电池输出端电压；I(k)为k时刻的回路电流；a,b,c为离散模型参数。

采用含遗忘因子的递推最小二乘法(FFRLS)，进行电池后向差分模型的参数辨识，求得a、b、c、U_OC(k)的值。根据方程(2)，对应得到k时刻的电池端电压的预测值U(k)。

定义SOC如下：

$SOC\left(t\right)={\mathrm{SOC}}_0-1/Q_N{\∫}_0^t{\mathrm{\η}}_{0}i\left(t\right)dt---\left(3\right)$

式中：SOC₀是SOC的初始值；Q_N为电池的额定容量；η₀为电池的充放电库伦效率，本例η₀＝1；i(t)为放电电流。

结合电池的离散模型，对电池特性描述如下，

状态方程 ${\mathrm{SOC}}_k={\mathrm{SOC}}_{k-1}-\frac{i_{k-1}T}{Q_N}---\left(4\right)$

量测方程U_K＝U_OC(SOC_k)-R₀i_k-U_p,k(5)

式中：U_p,k为k时刻极化电容两端电压；T为采样周期；U_OC(SOC_k)表示电池的开路电压和SOC的非线性关系。

U_OC(SOC_k)与SOC关系如下：

U_oc(SOC_k)＝a₁SOC_k ⁸+a₂SOC_k ⁷+a₃SOC_k ⁶+a₄SOC_k ⁵+

a₅SOC_k ⁴+a₆SOC_k ³+a₇SOC_k ²+a₈SOC_k+a₉(6)

根据实验得到的SOC和在线辨识得到的开路电压U_OC，运用含遗忘因子的递推最小二乘法(FFRLS)，求出式(6)中的a₁～a₉。

根据式(5)，由SOC的预测值SOC_k得到端电压的预测值U_k。

本方法主要包括如下步骤：

步骤Ⅰ、参数初始化

设定SOC的初始值和滤波协方差初始值，分别记作x_0|0和P_0|0，设定采样粒子数为N，本例N＝50，有效粒子数阈值N_thr＝0.78N＝39。

步骤Ⅱ、重要性采样

从高斯分布N(x_k-1|k-1,P_k-1|k-1)中随机采样N个粒子记作 为从N(x_k-1|k-1,P_k-1|k-1)中采样得到的大小为N的以x_k-1|k-1为中心值、P_k-1|k-1为方差的高斯分布；x_k-1|k-1和P_k-1|k-1分别为k-1时刻滤波均值和协方差，即SOC_k-1及其协方差。

Ⅱ-1、状态预测

采用步骤Ⅰ中的式(4)

状态方程 ${\mathrm{SOC}}_k={\mathrm{SOC}}_{k-1}-\frac{i_{k-1}T}{Q_N}$

进行状态预测，即求得：其中，为k-1时刻对SOC的第i个采样点；为k时刻从预测粒子集采样的第i个采样点；i_k-1为k-1时刻的电流，放电时电流值为正，充电时电流值为负；T为采样周期，Q_N为电池总容量。

Ⅱ-2、计算状态预测值的均值x_k|k-1和协方差P_k|k-1

$x_{k|k-1}=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}{\hat{x}}_k^i---\left(7\right)$

$P_{k|k-1}=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}({\hat{x}}_k^i-x_{k|k-1}){({\hat{x}}_k^i-x_{k|k-1})}^T---\left(8\right)$

Ⅱ-3、重新采样、重构采样分布函数

以状态预测值的均值x_k-1|k-1和协方差P_k-1|k-1重构高斯分布得N(x_k|k-1,P_k|k-1)，并重新采样 ${\left\{{\tilde{x}}_k^i\right\}}_{i=1}^N~N(x_{k|k-1},P_{k|k-1}),$

为从高斯分布N(x_k|k-1,P_k|k-1)中采样的N个粒子。

步骤Ⅲ、计算粒子权值并归一化。

Ⅲ-1、计算电池端电压预测值

将粒子代入量测方程(5)中计算电池端电压预测值 表示k时刻粒子对应的端电压预测值集合。

即： $y_{1,k}^i=U_{oc}\left({\tilde{x}}_k^i\right)-R_0{i}_k-U_{p,k},$

由式(6)计算得到。

Ⅲ-2、计算粒子权值

k时刻第i个粒子的权值由高斯分布的似然度函数计算得到，即 ${\mathrm{\ω}}_k^i=1/\sqrt{2\mathrm{\π}R}\exp (-1/(2R\left){\left(y_{1,k}-y_{1,k}^i\right)}^2\right),$ R为高斯似然度函数的方差，根据经验赋值。y_1,k为k时刻的量测值，表示k时刻第i个粒子对应的端电压预测值。

Ⅲ-3、权值归一化和计算有效粒子数N_eff

${\mathrm{\ω}}_k^i=\frac{{\widehat{\mathrm{\ω}}}_k^i}{\underset{i=1}{\overset{N}{\Σ}}{\widehat{\mathrm{\ω}}}_k^i}---\left(9\right)$

$N_{eff}=\frac{1}{\underset{i=1}{\overset{N}{\Σ}}{\left({\mathrm{\ω}}_k^i\right)}^2}---\left(10\right)$

为归一化后的权值；N_eff为相对于采样粒子数N而言的有效粒子数。 $N_{eff}=N/(1+Var({\mathrm{\ω}}_k^i\left)\right),$ 为的方差。

$x_{k|k}=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_k^i{\tilde{x}}_k^i---\left(11\right)$

$P_{k|k}=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_k^i({\tilde{x}}_k^i-x_{k|k}){({\tilde{x}}_k^i-x_{k|k})}^T---\left(12\right)$

x_k|k为状态估计值；P_k|k为协方差估计值。

步骤Ⅳ、有效粒子数N_eff与有效粒子数阈值N_thr的比较

高斯粒子滤波的预测粒子分布即先验概率分布π(x_k|k-1)，代表粒子数分布；后验概率分布即粒子实际概率分布，代表粒子的权值分布；先验概率分布和后验概率分布存在较大偏差。

拉普拉斯分布与高斯后验概率分布比较，峰值相同情况的两种分布曲线，在峰值两侧幅值均逐渐降低。将曲线与横轴间的高度定义为拖尾厚度。在接近横轴时，拉普拉斯分布的幅值下降比高斯后验概率分布的幅值下降更缓慢，即拉普拉斯分布的拖尾厚度大于高斯后验概率分布的拖尾厚度。

与高斯后验概率分布峰值相同的拉普拉斯分布，无论是拉普拉斯分布与预测粒子分布最大峰值中心点重合或不重合，拉普拉斯分布与预测粒子分布的重叠都比高斯后验概率分布多，，如图3和4所示，即预测粒子分布中有效粒子数增大。增大了滤波的协方差P_k|k-1，下一时刻按照N(x_k|k-1,P_k|k-1)进行抽样，采样范围更大，增强了滤波的稳定性。

当有效粒子数N_eff小于有效粒子数阈值N_thr，进入步骤Ⅴ；否则进入步骤Ⅵ。

步骤Ⅴ、更换似然度函数

采用拉普拉斯分布作为似然度函数。

步骤Ⅴ-1、计算有效粒子数退化率θ＝(N_thr-N_eff)/N_thr。N_thr为有效粒子数阈值；N_eff表示有效粒子数；本例N_thr＝0.75N，

步骤Ⅴ-2、计算方差调节因子γ＝1+εθ；

当γ增大时拉普拉斯分布的拖尾厚度增加，如图5所示。ε表示工况适应因子，本例ε＝0.5；

步骤Ⅴ-3、用代替拉普拉斯分布的概率密度函数中的x-μ，得到修改后的拉普拉斯分布似然度函数

${\mathrm{\ω}}_k^i=1/\left(2\mathrm{\δ}R\right)\exp (-|y_{1,k}-y_{1,k}^i|/(\mathrm{\δ}R\left)\right)$

其中x、μ、R的赋值分别参照高斯分布的似然度函数中的y_1,k、R，然后进行归一化。

步骤Ⅵ、状态更新和协方差更新

即执行步骤Ⅲ-3的式(11)(12)，得到更新的SOC估计值和更新的协方差。

步骤Ⅶ、判断滤波时间大于设定的运行时间，则结束，否则，k＝k+1，返回步骤Ⅱ。

似然度函数粒子滤波的动力电池电荷状态估计系统实施例

本似然度函数粒子滤波的动力电池电荷状态估计系统实施例如图10所示，本例的微控制器为基于x86系统的嵌入式微控制器，如PC104板。本例微控制器的中心处理器连接程序存储器和数据存储器，程序存储器含有电压电流数据采集模块、根据参数辨识的U_OC(SOC_k)与SOC关系模块，端电压估计及端电压预测偏差计算模块，有效粒子数N_eff与有效粒子数阈值N_thr的比较模块，方差调节因子、有效粒子数退化率及修改后的拉普拉斯分布似然度函数计算模块，状态更新和协方差更新模块以及SOC显示模块；数据存储器存放初始化参数。

电压传感器和电流传感器接入模数转换电路，再经通用输入输出接口(通用I/O接口)连接微控制器的中心处理器。微控制器的中心处理器连接显示器。

中心微处理器连接CAN总线接口、PC104总线接口和以太网接口，以便功能扩展使用。

在不同SOC初值和不同粒子数N的情况下，本似然度函数粒子滤波的动力电池电荷状态估计方法PF-LFAF和高斯粒子滤波GPF的动力电池电荷状态估计方法的所得SOC结果以及SOC实际值的对比如图6至图9所示。

传统实验所得的该型号动力电池的SOC作为SOC实际值。

图6至9的纵坐标为SOC值，横坐标为时间，大图单位为10⁴秒，局部放大的右上角小图单位为秒。其中实线表示SOC实际值，短线虚线表示本似然度函数粒子滤波的动力电池电荷状态估计方法PF-LFAF所得SOC，点虚线表示高斯粒子滤波GPF的动力电池电荷状态估计方法的所得SOC。

图6SOC初值SOC₀＝0.8、粒子数N＝50；

图7SOC初值SOC₀＝0.2、粒子数N＝50；

图8SOC初值SOC₀＝0.8、粒子数N＝100；

图9SOC初值SOC₀＝0.2、粒子数N＝100。

图7和9中可清楚看到本实施例所得SOC曲线短线虚线与实际SOC曲线实线最为接近。各图的右上角小图中可看到当本实施例电荷状态初值与SOC实际值相差较大的情况下，本实施例PF-LFAF法的收敛速度要快于GPF法。

为了进一步说明本实施例方法相对于高斯粒子滤波GPF的动力电池电荷状态估计方法具有良好的适应性，列出不同SOC初始值相同粒子数下的实验数据结果，见表1和表2。

表1粒子数为50时两种方法估算误差对照表

N	SOC0	PF-LFAF平均误差	GPF平均误差
				50	1	0.0512％	0.1469％
50	0.8	0.0527％	0.1574％
				50	0.6	0.1374％	发散
50	0.4	0.1427％	发散
				50	0.2	0.2797％	发散

表2粒子数为100时两种方法估算误差对照表

N	SOC0	PF-LFAF平均误差	GPF平均误差
				100	1	0.0476％	0.0702％
100	0.8	0.0491％	0.1442％
				100	0.6	0.0636％	0.1608％
100	0.4	0.0888％	发散
				100	0.2	0.2045％	发散

由表1可知，粒子数为50时本例PF-LFAF法的估算精度要优于GPF法，且当SOC初始值偏差分别为40％、60％、80％时PF-LFAF法依然收敛，保持较高精度。

上述实施例，仅为对本发明的目的、技术方案和有益效果进一步详细说明的具体个例，本发明并非限定于此。凡在本发明的公开的范围之内所做的任何修改、等同替换、改进等，均包含在本发明的保护范围之内。

Claims (8)

1.一种似然度函数粒子滤波的动力电池电荷状态估计方法，采用Thevenin模型描述电池的静态和动态性能；电池的极化电阻R_p与电池的极化电容C_p并联构成一阶RC结构，RC两端电压为U_p(t)；串接欧姆电阻R₀和Uoc，Uoc为电池的开路电压OCV，采样得到电池端电压U(t)和流过欧姆内阻R₀的电流i(t)；

电池Thevenin模型数学表达式如下：

$\left\{\begin{array}{c}\frac{{\mathrm{dU}}_p\left(t\right)}{dt}=-\frac{U_p\left(t\right)}{R_p{C}_p}+\frac{i\left(t\right)}{C_p}\\ U\left(t\right)=U_{OC}\left(t\right)-R_{0}i\left(t\right)-U_p\left(t\right)\end{array}\right.---\left(1\right)$

用后向差分法对Thevenin电池模型离散化，整理后得到对应的后向差分离散模型：

U(k)＝aU(k-1)+bI(k)+cI(k-1)+(1-a)U_OC(k)(2)

式中：U_OC(k)表示k时刻的开路电压；U(k)为k时刻的电池输出端电压；I(k)为k时刻的回路电流；a,b,c为离散模型参数；

采用含遗忘因子的递推最小二乘法，进行电池后向差分模型的参数辨识，求得a、b、c、U_OC(k)的值，根据方程(2)，对应得到k时刻的电池端电压的预测值U(k)；

定义SOC如下：

$SOC\left(t\right)={\mathrm{SOC}}_0-1/Q_N{\∫}_0^t{\mathrm{\η}}_{0}i\left(t\right)dt---\left(3\right)$

式中：SOC₀是SOC的初始值；Q_N为电池的额定容量；η₀为电池的充放电库伦效率，本发明η₀＝1；i(t)为放电电流；

结合电池的离散模型，对电池特性描述如下，

状态方程 ${\mathrm{SOC}}_k={\mathrm{SOC}}_{k-1}-\frac{i_{k-1}T}{Q_N}---\left(4\right)$

量测方程U_K＝U_OC(SOC_k)-R₀i_k-U_p,k(5)

式中：U_p,k为k时刻极化电容两端电压；T为采样周期；U_OC(SOC_k)表示电池的开路电压和SOC的非线性关系；

U_OC(SOC_k)与SOC关系如下：

U_oc(SOC_k)＝a₁SOC_k ⁸+a₂SOC_k ⁷+a₃SOC_k ⁶+a₄SOC_k ⁵+

a₅SOC_k ⁴+a₆SOC_k ³+a₇SOC_k ²+a₈SOC_k+a₉(6)

根据实验得到的SOC和在线辨识得到的开路电压U_OC，运用含遗忘因子的递推最小二乘法(FFRLS)，求出式(6)中的a₁～a₉；

根据式(5)，由SOC的预测值SOC_k得到端电压的预测值U_k；

其特征在于主要包括如下步骤：

步骤Ⅰ、参数初始化

设定SOC的初始值和滤波协方差初始值，分别记作x_0|0和P_0|0，设定采样粒子数为N，N取值为30～100；有效粒子数阈值为N_thr，设定N_thr＝(0.5～0.85)N；

步骤Ⅱ、重要性采样

从正态分布，即高斯分布N(x_k-1|k-1,P_k-1|k-1)中随机采样N个粒子记作为从N(x_k-1|k-1,P_k-1|k-1)中采样得到的大小为N的以x_k-1|k-1为中心值、P_k-1|k-1为方差的高斯分布；x_k-1|k-1和P_k-1|k-1分别为k-1时刻滤波均值和协方差，即SOC_k-1及其协方差；

Ⅱ-1、状态预测

采用步骤Ⅰ中的式(4)

状态方程 ${\mathrm{SOC}}_k={\mathrm{SOC}}_{k-1}-\frac{i_{k-1}T}{Q_N}$

进行状态预测，即求得： ${\hat{x}}_k^i=x_{k-1}^i-i_{k-1}T/Q_N;$

其中，为k-1时刻的对SOC的第i个采样点；为k时刻从预测粒子集采样的第i个采样点；i_k-1为k-1时刻的电流，放电时电流值为正，充电时电流值为负；T为采样周期，Q_N为电池总容量；

Ⅱ-2、计算状态预测值的均值x_k|k-1和协方差P_k|k-1

$x_{k|k-1}=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}{\hat{x}}_k^i---\left(7\right)$

$P_{k|k-1}=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}({\hat{x}}_k^i-x_{k|k-1}){({\hat{x}}_k^i-x_{k|k-1})}^T---\left(8\right)$

Ⅱ-3、重新采样、重构采样分布函数

以状态预测值的均值x_k-1|k-1和协方差P_k-1|k-1重构高斯分布得N(x_k|k-1,P_k|k-1)，并重新采样

为从高斯分布N(x_k|k-1,P_k|k-1)中采样的N个粒子；

步骤Ⅲ、计算粒子权值并归一化

Ⅲ-1、计算电池端电压预测值

将粒子代入量测方程(5)中计算电池端电压预测值 表示k时刻粒子对应的端电压预测值集合；

即： $y_{1,k}^i=U_{oc}\left({\tilde{x}}_k^i\right)-R_0{i}_k-U_{p,k},$

由式(6)计算得到；

Ⅲ-2、计算粒子权值

k时刻第i个粒子的权值由高斯分布的似然度函数计算得到，即 ${\mathrm{\ω}}_k^i=1/\sqrt{2\mathrm{\π}R}\exp (-1/(2R\left){\left(y_{1,k}-y_{1,k}^i\right)}^2\right),$ R为高斯似然度函数的方差，根据经验赋值；y_1,k为k时刻的量测值，表示k时刻第i个粒子对应的端电压预测值；

Ⅲ-3、权值归一化和计算有效粒子数N_eff

${\mathrm{\ω}}_k^i=\frac{{\widehat{\mathrm{\ω}}}_k^i}{\underset{i=1}{\overset{N}{\Σ}}{\widehat{\mathrm{\ω}}}_k^i}---\left(9\right)$

$N_{eff}=\frac{1}{\underset{i=1}{\overset{N}{\Σ}}{\left({\mathrm{\ω}}_k^i\right)}^2}---\left(10\right)$

为归一化后的权值；N_eff为相对于采样粒子数N而言的有效粒子数； $N_{eff}=N/(1+Var({\mathrm{\ω}}_k^i\left)\right)$ 为的方差；

$x_{k|k}=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_k^i{\tilde{x}}_k^i---\left(11\right)$

$P_{k|k}=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_k^i({\tilde{x}}_k^i-x_{k|k}){({\tilde{x}}_k^i-x_{k|k})}^T---\left(12\right)$

x_k|k为状态估计值；P_k|k为协方差估计值；

步骤Ⅳ、有效粒子数N_eff与有效粒子数阈值N_thr的比较

高斯粒子滤波的预测粒子分布即先验概率分布π(x_k|k-1)，代表粒子数分布；后验概率分布即粒子实际概率分布，代表粒子的权值分布；

当N_eff小于N_thr，进入步骤Ⅴ；否则进入步骤Ⅵ；

步骤Ⅴ、更换似然度函数

采用拉普拉斯分布作为似然度函数；

拉普拉斯分布的概率密度函数，或称为似然度函数如下：

$f\left(x\right|\mathrm{\μ},R)=\frac{1}{2R}\exp (-\frac{\left|x-\mathrm{\μ}\right|}{R})---\left(13\right)$

式中：x为预测方程所得状态预测值，或者x为量测方程所得量测值对应物理量的预测值，μ和R分别是该拉普拉斯分布的中心值和方差；

借鉴高斯滤波的方法，用代替拉普拉斯分布的概率密度函数(13)中的x-μ，得到修改后的拉普拉斯分布似然度函数

${\mathrm{\ω}}_k^i=1/\left(2\mathrm{\δ}R\right)\exp (-|y_{1,k}-y_{1,k}^i|/(\mathrm{\δ}R\left)\right);$

其中x、μ、R的赋值分别参照高斯分布的似然度函数中的y_1,k、R，然后进行归一化；

步骤Ⅵ、状态更新和协方差更新

即执行步骤Ⅲ-3的式(11)(12)，得到更新的SOC估计值和更新的协方差；

步骤Ⅶ、判断滤波时间大于设定运行时间，则滤波结束；否则，k＝k+1，返回步骤Ⅱ。

2.根据权利要求1所述的似然度函数粒子滤波的动力电池电荷状态估计方法，其特征在于：

所述步骤Ⅰ中设定采样粒子数N取值为45～55。

3.根据权利要求1所述的似然度函数粒子滤波的动力电池电荷状态估计方法，其特征在于：

所述步骤Ⅰ中设定有效粒子数阈值N_thrr＝(0.73～0.81)N。

4.根据权利要求1所述的似然度函数粒子滤波的动力电池电荷状态估计方法，其特征在于：

所述步骤Ⅴ更换似然度函数时引入方差调节因子γ，γ＝1+θ；θ为有效粒子数退化率，θ＝(N_thr-N_eff)/N_thr；

引入γ后，方程(13)改写为：

$\begin{array}{c}f\left(x\right|\mathrm{\μ},R\mathrm{\δ})=\frac{1}{2\left(R\mathrm{\γ}\right)}\exp (-\frac{\left|x-\mathrm{\μ}\right|}{R\mathrm{\γ}})\\ =\frac{1}{2\left(R\right(1+\left(N_{thr}-N_{eff}\right)/N_{thr}\left)\right)}\exp (-\frac{\left|x-\mathrm{\μ}\right|}{R\left(1+\left(N_{thr}-N_{eff}\right)/N_{thr}\right)})\end{array}---\left(14\right)$

当γ增大时拉普拉斯分布的拖尾厚度增加。

5.根据权利要求4所述的似然度函数粒子滤波的动力电池电荷状态估计方法，其特征在于：

所述方差调节因子γ的定义式中引入工况适应因子ε，γ＝1+εθ，ε～(0,1]，代入式(14)得：

$\begin{array}{c}f\left(x\right|\mathrm{\μ},R\mathrm{\γ})=f\left(x\right|\mathrm{\μ},R\mathrm{\ϵ})\\ =\frac{1}{2\left(R\right(1+\mathrm{\ϵ}\left(N_{thr}-N_{eff}\right)/N_{thr}\left)\right)}\exp (-\frac{\left|x-\mathrm{\μ}\right|}{R\left(1+\mathrm{\ϵ}\left(N_{thr}-N_{eff}\right)/N_{thr}\right)})\end{array}---\left(15\right)$

动力电池工况负载变化不剧烈时，SOC变化相对较慢，似然度函数的拖尾厚度小，ε取较小值；

动力电池工况剧烈时，SOC相对变化较快，似然度函数的拖尾厚度大，ε取较大值。

6.根据权利要求5所述的似然度函数粒子滤波的动力电池电荷状态估计方法，其特征在于：

所述工况适应因子ε＝0.4～0.6。

7.根据权利要求1至6中任一项所述的似然度函数粒子滤波的动力电池电荷状态估计方法所设计的似然度函数粒子滤波的动力电池电荷状态估计系统；包括微控制器、电压传感器和电流传感器，微控制器还连接显示器；微控制器配有通用输入输出接口，电压传感器和电流传感器接入模数转换电路，再经通用接口连接微控制器；其特征在于：

所述微控制器的中心处理器连接有程序存储器和数据存储器，程序存储器含有电压电流数据采集模块、根据参数辨识的U_OC(SOC_k)与SOC关系模块，端电压估计及端电压预测偏差计算模块，有效粒子数N_eff与有效粒子数阈值N_thr的比较模块，方差调节因子、有效粒子数退化率及修改后的拉普拉斯分布似然度函数计算模块，状态更新和协方差更新模块以及SOC显示模块；数据存储器存放初始化参数。

8.根据权利要求7所述的似然度函数粒子滤波的动力电池电荷状态估计系统；其特征在于：

所述微控制器的中心处理器连接CAN总线接口、PC104总线接口和以太网接口中的一种或多种。