CN104181470A - 一种基于非线性预测扩展卡尔曼滤波的电池soc估计方法 - Google Patents

Abstract

本发明公开了一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法，它具有运算简单与精度高的优点；它包括步骤：1)将整个方法所需要的时间平均分成N个时间段，每一个时间段代表一步，即第k个时间段为第k步；其中k≤N，k与N均为正整数；2)根据基尔霍夫电压和电流定理，建立电池系统的模型；得到电池系统模型的误差矩阵d(k)与电池系统的模型误差分配矩阵G(k)；3)对步骤2)所得的先验状态估计方程进行补偿；4)求解第k+1步电池系统的后验状态估计方程，得到SOC值；5)根据步骤4)的第k+1步电池系统的后验状态估计结果，与电池的SOC真值进行比较，验证非线性预测卡尔曼滤波NPEKF算法的有效性，k＝n+1，并转到步骤3)，直至第N步执行完。

Description

一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法

技术领域

本发明涉及电池的电荷状态SOC估计方法，尤其涉及一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法。

背景技术

车载动力电池作为电动汽车的关键部件，其性能对整车的动力性、经济性和安全性至关重要，是制约电动汽车规模发展的关键因素。为最大限度地发挥动力电池的性能并延长电池的使用寿命，对电池进行有效管理至关重要，而准确获得电池的荷电状态(state of charge，SOC)是电池管理最核心的技术。电池SOC估计是判断电池是否过充过放，是否需要均衡或更换某一单体电池的重要依据。因此，提高电池SOC估计精度，对于提高电池使用效率、延长电池循环寿命以及保障电池安全可靠，有着重要的意义。

目前常用的电池SOC估计方法可以分为以下四类：①基于安时计量的估计方法，主要有放电实验法和Ah计量法；②基于电池表征参数测量值的估计方法，主要有开路电压法和内阻法；③基于经验方程和数学模型的估计方法，主要有线性模型法、神经网络法、模糊逻辑法、支持向量机法和相关向量机法；④基于电池等效电路模型的估计方法，主要有卡尔曼滤波(KF)法、观测器法和粒子滤波法。近几年新兴的基于电池等效电路模型的估计方法因其模型物理意义明确，易于参数辨识，已成为电池研究的热点。其中扩展卡尔曼滤波(EKF)算法备受关注，它适应于电流波动比较剧烈的SOC估计，非常适合应用在电动汽车上。但是EKF有效应用的前提条件是：

1)要求系统模型准确已知；

2)过程噪声被假设为协方差已知的零均值高斯白噪声过程。

在基于EKF的电池SOC估计中，电池动态模型中存在的误差被作为过程噪声来处理，而且一般假定为零均值、协方差已知的高斯白噪声。然而，这种假设往往缺少理论根据，不能正确反映电池模型误差的真实特性，从而导致SOC估计精度下降，甚至使滤波器发散。

发明内容

本发明的目的就是为解决上述问题，提出了一种基于非线性预测扩展卡尔曼滤波(NPEKF)的电池SOC估计方法，能够在电池模型存在较大模型误差的情况下准确估计电池SOC，有效地解决存在显著模型误差的非线性滤波问题，与传统基于EKF的SOC估计算法相比较，NPEKF算法具有以下优点：1)过程噪声的统计特性无需事先准确知道(不必满足高斯分布)，对于未知的电池模型误差，作为滤波器解的一部分可以实时地得到它的一步预测估计值；2)电池模型误差的结构形式不受限制，可以是线性的，也可以是非线性的；3)省去了由于协方差的求解而带来的运算负担，而且不必计算非线性状态函数与测量函数的雅可比矩阵。

为了实现上述目的，本发明所采用的技术方案为：

一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法，包括步骤：

1)将整个方法所需要的时间平均分成N个时间段，每一个时间段代表一步，即第k个时间段为第k步；其中k≤N，k与N均为正整数；

2)根据基尔霍夫电压和电流定理，建立电池系统的模型；以电池的电流为输入变量，以SOC值与支路电压为状态变量，得到电池系统的先验状态估计方程；以电池两端的电压为输出变量，以电池的电流为输入变量，以SOC值与支路电压为状态变量，得到电池系统的测量方程；根据电池系统模型的测量方程，得到电池系统模型的误差矩阵d(k)与电池系统的模型误差分配矩阵G(k)；

3)对步骤2)所得的先验状态估计方程进行补偿：利用电池系统的模型误差矩阵d(k)与电池系统的模型误差分配矩阵G(k)对步骤2)所得的先验状态估计方程进行补偿，得到精准的先验状态估计方程，k＝n，n为正整数，n≤N；

4)求解第k+1步电池系统的后验状态估计方程，得到SOC值；

5)根据步骤4)得到的第k+1步电池系统的后验状态估计结果，与电池的SOC真值进行比较，验证非线性预测卡尔曼滤波NPEKF算法的有效性，k＝n+1，并转到步骤3)，直至第N步执行完。

所述步骤3)中对先验状态估计方程进行补偿的具体步骤为：

31)对电池系统的模型误差分配矩阵G(k)进行初始化；

32)求解电池系统的模型误差矩阵d(k)；

33)对步骤2)所得的先验状态估计方程进行补偿：

${\hat{x}}_{k+1}^{-}=f({\hat{x}}_k,u_k)+G\left(k\right)d\left(k\right);$

其中为第k+1步的先验状态估计，为关于的函数，x_k为状态变量；u_k为输入变量。

所述步骤31)还包括对NPEKF滤波器的初始化：

给定状态初始估计值给定电池系统的模型误差加权矩阵W、过程激励噪声协方差矩阵Q、测量噪声方差矩阵R和误差协方差矩阵P₀的初始值。

所述步骤3)中电池系统的模型误差矩阵d(k)的具体求解过程为：

321)求解参数Λ(Δt)：

Λ(Δt)是1×1维的对角阵，其值为：

Λ(Δt)＝Δt²/2；

322)求解参数

的计算公式如下：

$Z[\hat{x}\left(t\right),\mathrm{\Δt}]=\mathrm{\Δt}{\×L}_f^1+({\mathrm{\Δt}}^2/2)\×L_f^2;$

式中，分别是关于的第1阶和第2阶李导数，为输出变量估计关于的函数，为第k+1步的状态估计关于的函数，输入变量u_k＝i_k，i_k为第k个时刻点的电池端电流；

323)求解参数

是1×2维的矩阵，其行元素如下：

$S\left[{\hat{x}}_k\right]=\left\{L_{g_1}\right[L_f^1\left(h({\hat{x}}_k,u_k)\right)],L_{g_2}[L_f^1\left(h({\hat{x}}_k,u_k)\right)\left]\right\}$

式中，g₁和个g₂分别是G(k)中的第1和第2列元素；是先沿后沿g₁的李导数；是先沿后沿g₂的李导数。

根据上述公式可得到：

$S\left[{\hat{x}}_k\right]=\left\{L_{g_1}\right[L_f^1\left(h({\hat{x}}_k,u_k)\right)],L_{g_2}[L_f^1\left(h({\hat{x}}_k,u_k)\right)\left]\right\}$

式中，U_oc,_k为第k步电池开路电压，SOC为电池的电荷状态；

324)求解模型误差矩阵d(k)：

$\begin{array}{c}d\left(k\right)=-{\{{\left[\mathrm{\Λ}\left(\mathrm{\Δt}\right)S\left({\hat{x}}_k\right)\right]}^T{R}^{-1}\mathrm{\Λ}\left(\mathrm{\Δt}\right)S\left({\hat{x}}_k\right)+W\}}^{-1}\\ \×{\left[\mathrm{\Λ}\left(\mathrm{\Δt}\right)S\left({\hat{x}}_k\right)\right]}^T{R}^{-1}\left[Z({\hat{x}}_k,\mathrm{\Δt}){+{\hat{y}}_k-y_{k+1}}_{}\right]\end{array}$

其中，R为测量噪声协方差矩阵，W为加权矩阵，Δt＝t_k+1-t_k为采样时间间隔，为第k步的状态估计，为第k步的y_k的估计，y_k+1为第k+1步的测量真值，y_k+1＝U_bat,k+1，U_bat,k+1第k+1步的电池端电压。

所述步骤4)求解电池系统的后验状态估计方程的具体步骤为：

41)计算第k+1步的增益矩阵K_k+1；

42)求解后验估计

${\hat{y}}_{k+1}=h({\hat{x}}_{k+1}^{-},u_k)+v_{k+1};$

${\hat{x}}_{k+1}={\hat{x}}_{k+1}^{-}+K_{k+1}[y_{k+1}-{\hat{y}}_{k+1}];$

其中，为第k+1步测量值估计值，为，y_k+1为第k+1步的测量真值，具体指代为第k+1步电池的端电压U_bat,k+1，u_k+1＝i_k+1，i_k+1为第k+1步的电池电流，为第k+1步状态的先验估计。

所述步骤41)中第k+1步的增益矩阵K_k+1的具体求解过程为：

411)计算第k+1步误差协方差矩阵的先验估计

${\hat{P}}_{k+1}^{-}=A_{k+1}{\hat{p}}_k{A_{k+1}}^T+w_{k+1}{Q}_k{w}_{k+1}^T;$

其中，A_k+1为第k+1步的系统矩阵，为第k步的误差协方差矩阵的后验估计，Q_k为第k步的过程激励噪声协方差矩阵，w_k+1表示过程激励噪声，为正态分布的白色噪声：w_k+1～N(0,R_k+1)；

412)计算第k+1步的增益矩阵K_k+1：

$K_{k+1}={\hat{P}}_{k+1}^{-}{H}_{k+1}^T{[H_{k+1}{\hat{P}}_{k+1}^{-}{H}_{k+1}^T+v_{k+1}{R}_{k+1}{v}_{k+1}^T]}^{-1};$

其中，为第k+1步的误差协方差矩阵的先验估计，H_k+1为第k+1步的测量矩阵，v_k+1为过程观测噪声，为正态分布v_k+1～N(0,R_k+1)。

所述步骤411)第k步的误差协方差矩阵的后验估计的具体求解为：

${\hat{P}}_k=(1-K_k{H}_k){\hat{P}}_{k+1}^{-};$

其中K_k为第k步的增益矩阵，H_k为第k步的测量矩阵，为第k步的误差协方差矩阵的先验估计。

本发明的有益效果为：

1、对于未知的模型误差，该算法能够将其作为滤波器解的一部分，可以实时地得到它的一步预测估计值，并用于修正电池的动态模型，提高了SOC估计精度；

2、模型误差的结构形式不受限制；

3、算法运算简单；

4、有效地解决存在显著动态模型误差情况下的非线性电池系统SOC估计问题，其估计误差在4％以内；

5、适用于电池的各种复杂工况。

附图说明

图1为本发明的二阶RC电池等效电路模型图，其中下标d为放电方向的参数；下标c为充电方向的参数。

图2为本发明的系统流程图；

图3为本发明具体流程图；

图4为本发明的二阶RC模型在脉冲放电下的端电压及其误差曲线图。

图5为本发明的二阶RC模型在脉冲充电下的端电压及其误差曲线图。

图6为本发明的NPEKF滤波算法的流程图。

具体实施方式

为了更好的了解本发明的技术方案，下面结合附图对本发明作进一步说明。

1、二阶RC模型

利用非线性预测扩展卡尔曼滤波器估计电池SOC，需要精确建立电池模型。搭建电池模型是指应用数学理论尽量全面地去描述实际电池的响应特性和内部特性。所谓响应特性是指电池的端电压与负载电流的对应关系；内部特性是指电池的内部变量欧姆内阻、极化内阻和极化电压与SOC、温度间的关系。

如图1所示为本发明的二阶RC等效电路模型，包括运行时间电路和I-V特性电路。

所述运行时间电路包括电池的自放电电阻R_d、电容C_Q和电流控制电流源电路，电阻R_d与电容C_Q并联在电流控制电流源的受控源两端，独立电源的一端接地。所述I-V特性电路包括两路支路，每个支路包括三组一个电容与一个电阻并联组成的RC回路，且每个支路的两组RC回路并联有切换开关。所述I-V特性电路包括欧姆内阻R₀、电化学极化内阻R₁、电化学极化电容C₁、浓差极化内阻R₂、浓差极化电容C₂和电流控制电流源、电压控制电压源电路，其中：电压控制电压源电路的受控源的正极连接两路，一路连接二极管D_d后连接电阻R_1d、电阻R_2d和电阻R_od后连接电池的正极，一路反接二极管D_c后连接电阻R_1c、电阻R_2c和电阻R_oc后连接电池的正极，电容C_1d并联在电阻R_1d两端，电容C_2d并联在电阻R_2d两端，电容C_1c并联在电阻R_1c两端，电容C_2c并联在电阻R_2c两端，电压控制电压源电路的受控源正、负极之间的电压为电池开路电压OCV。

所述运行时间电路和I-V特性电路通过一个流控电流源和一个压控电压源建立联系，当对电池进行充放电时，负载电流i_bat通过流控电流源对电容C_Q进行充放电，从而改变C_Q存储的电量，从而表征电池SOC的变化，C_Q两端电压OCV也随之变化，进而控制I-V特性电路的压控电压源OCV＝f(SOC)随SOC的变化而变化。负载电流i_bat对电容C_Q进行放电时，当C_Q两端电压OCV达到放电终止电压时即可获得电池总的运行时间。

所述电容C_Q表示电池的可用容量，C_Q＝3600·C_Ah·f₁·f₂，CAh为用安时为单位的电池容量，f₁和f₂分别是电池循环寿命和温度的修正因子。

所述电流控制电流源的受控源的电流为电池的端电流i_bat，当电池进行充放电时负载电流i_bat通过电流控制电流源对电容C_Q进行充放电，改变电容C_Q中存储的电量，从而表征电池SOC的变化。

所述电流控制电流源的受控源两端的电压为电池开路电压OCV(或记为U_oc)，所述开路电压OCV与SOC存在较强的非线性关系，采用一个压控电压源来模拟(如图1所示)，具体关系式为：

$U_{\mathrm{OC}}=\mathrm{OCV}=a_1+a_2\ln \mathrm{SOC}+a_3\ln (1-\mathrm{SOC})+\frac{a_4}{\mathrm{SOC}}+a_5\mathrm{SOC}---\left(1\right)$

式中，a₁-a₅为常数，由实验数据基于最小二乘法辨识得到。

图2-3为本发明的流程图。

如图4-5所示分别为本发明的二阶RC模型在脉冲放电和充电下的端电压及其误差曲线图。从图中可以看出，采用本文提出的变阶数RC等效电路模型能很好地反应电池的脉冲充放电过程，这说明该模型是准确的，可以利用的。在恒流充放电的过程中所产生的误差要比在静止的阶段产生的误差要大一点，这是因为拟合的是零输入零状态响应。总体来说，静置阶段电池的误差几乎为零，充放电阶段电池的误差小于0.02V。

2、二阶RC电池模型的数学方程

以电池放电为例，在二阶RC模型的基础上，结合安时法，根据基尔霍夫电压和电流原理，建立了如下数学模型：

U_bat(t)＝U_OC(SOC)-i_bat(t)R₀-U₁(t)-U₂(t)   (2)

$\left\{\begin{array}{c}{U}_1\left(t\right)=i\left(t\right)R_1(1-e^{t/R_1{C}_1})\\ U_2\left(t\right)=i\left(t\right)R_2(1-e^{t/R_2{C}_2})\end{array}\right.---\left(3\right)$

式中，U_bat为电池的端电压；R₀为欧姆内阻；R₁为电化学极化内阻，C₁为电化学极化电容；R₂为浓差极化电阻，C₂为浓差极化电容；U_oc为电池的开路电压；i_bat(t)为电池的端电流，以下简记为i。

要利用NPEK计算动力电池的SOC，需要将SOC作为NPEK的状态变量之一。利用Ah法对SOC进行估算的离散化表达式为：

${\mathrm{SO}\hat{C}}_{k+1}^{-}={\mathrm{SO}\hat{C}}_k+\frac{\mathrm{\η}*\mathrm{\Δt}}{C}{i}_k+w_k---\left(4\right)$

式中，η为库仑系数，可以通过电池充放电试验得到，电池充电时，η＝1，放电时，η<1；i_k为采样时刻点k处的电流，具有符号性，充电时i_k>0，放电时i_k<0；C为电池标称容量，单位是Ah，Δt＝t_k+1-t_k是采样间隔。

根据图1所示的动力电池模型，选取RC网络中电容的电压和SOC作为系统状态变量，系统的状态方程为：

$\left[\begin{array}{c}{\mathrm{SOC}}_{k+1}^{-}\\ U_{1,k+1}^{-}\\ U_{2,k+1}^{-}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& \exp (-\mathrm{\&Delta;t}/R_1{C}_1)& 0\\ 0& 0& \exp (-\mathrm{\&Delta;t}/R_2{C}_2)\end{array}\right]*\left[\begin{array}{c}{\mathrm{SOC}}_k\\ U_{1,k}\\ U_{2,k}\end{array}\right]+\left[\begin{array}{c}-\mathrm{\&eta;}*\mathrm{\&Delta;t}/C\\ R_1*(1-\mathrm{esp}(-\mathrm{\&Delta;t}/R_1{C}_1))\\ R_2*(1-\exp (-\mathrm{\&Delta;t}/R_2{C}_2))\end{array}\right]*i_k---\left(5\right)$

上述方程可记为：

${\hat{x}}_{k+1}^{-}=A_k{x}_k+B_k{u}_k=f(x_k,u_k)---\left(6\right)$

其中，状态变量为：

$x_k=\left[\begin{array}{c}{\mathrm{SOC}}_k\\ U_{1,k}\\ U_{2,k}\end{array}\right]---\left(7\right)$

输入变量为：

u_k＝i_k   (8)

系统矩阵为：

$A_k=\frac{\&PartialD;f(x_k,u_k)}{\&PartialD;x}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{esp}(-\mathrm{\&Delta;t}/R_1{C}_1)& 0\\ 0& 0& \mathrm{esp}(-\mathrm{\&Delta;t}/R_2{C}_2)\end{array}\right]---\left(9\right)$

输入控制矩阵为：

$B_k=\frac{\&PartialD;f(x_k,u_k)}{\&PartialD;u}=\left[\begin{array}{c}-\mathrm{\&eta;}*\mathrm{\&Delta;t}/C\\ R_1*(1-\exp (-\mathrm{\&Delta;t}/R_1{C}_1))\\ R_2*(1-\exp (-\mathrm{\&Delta;t}/R_2{C}_2))\end{array}\right]---\left(10\right)$

根据图1所示的动力电池模型，可得到电池系统的测量方程的离散形式为：

${\hat{U}}_{\mathrm{bat},k+1}={\left[\begin{array}{c}0\\ -1\\ -1\end{array}\right]}^T*\left[\begin{array}{c}{\mathrm{SOC}}_{k+1}^{-}\\ U_{1,k+1}^{-}\\ U_{2,k+1}^{-}\end{array}\right]-R_0*i_{k+1}+U_{\mathrm{OC},k}---\left(11\right)$

上述(11)式可记为：

${\hat{y}}_{k+1}=H_{k+1}{\hat{x}}_{k+1}^{-}+D_{k+1}{u}_{k+1}+v_{k+1}=h({\hat{x}}_{k+1}^{-},u_{k+1})+v_{k+1}---\left(12\right)$

式中，随机信号v_k+1表示过程观测噪声，假设为正态分布的白色噪声，且v_k+1～N(0,R_k+1)，R_k+1为第k+1步的观测噪声协方差。

测量矩阵为：

$H_{k+1}=\frac{\&PartialD;h({\hat{x}}_{k+1}^{-},)}{\&PartialD;x}=\left[\begin{array}{ccc}\frac{{\mathrm{dU}}_{\mathrm{OC}}}{\mathrm{dSOC}}{|}_{{\mathrm{SOC}}_{k+1}={\mathrm{SO}\hat{C}}_{k+1}}& -1& -1\end{array}\right]---\left(13\right)$

$D_{k+1}=\frac{\&PartialD;h({\hat{x}}_{k+1}^{-},u_{k+1})}{\&PartialD;u}=\left[R_0\right]---\left(14\right)$

输出变量为：

y_k+1＝U_bat,k+1   (15)

3、输出估计的泰勒级数近似

对方程(5)中的进行泰勒级数展开，可得到

${\hat{U}}_{\mathrm{bat},k+1}={\hat{U}}_{\mathrm{bat},k}+\mathrm{\&Delta;t}*{\stackrel{\stackrel{\&CenterDot;}{^}}{U}}_{\mathrm{bat},k}+\frac{{\mathrm{\&Delta;t}}^2}{2}*{\stackrel{\stackrel{\&CenterDot;\&CenterDot;}{^}}{U}}_{\mathrm{bat},k}+\&CenterDot;\&CenterDot;\&CenterDot;---\left(16\right)$

整理上式，得到

${\hat{U}}_{\mathrm{bat},k+1}={\hat{U}}_{\mathrm{bat},k}+Z[{\hat{x}}_k,\mathrm{\&Delta;t}]+\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left[{\hat{x}}_k\right]d\left(k\right)+\&CenterDot;\&CenterDot;\&CenterDot;---\left(17\right)$

式中，Λ(Δt)是1×1维的对角阵，其值为：

Λ(Δt)＝Δt²/2   (18)

的计算公式如下：

$Z[\hat{x}\left(t\right),\mathrm{\&Delta;t}]=\mathrm{\&Delta;t}\&times;L_f^1+({\mathrm{\&Delta;t}}^2/2)\&times;L_f^2---\left(19\right)$

式中，是关于的m阶李导数，定义为：

$\left\{\begin{array}{c}{L}_f^m\left(h({\hat{x}}_k,u_k)\right)=h({\hat{x}}_k,u_k),m=0\\ L_f^m\left(h({\hat{x}}_k,u_k)\right)=\frac{{\&PartialD;L}_f^{m-1}\left(h({\hat{x}}_k,u_k)\right)}{\&PartialD;\hat{x}}f({\hat{x}}_k,u_k),m\&GreaterEqual;1\end{array}\right.---\left(20\right)$

根据上述定义可以得到：

$L_f^1=\frac{\&PartialD;h({\hat{x}}_k,u_k)}{\&PartialD;\hat{x}}f\left(\hat{x}\right)=\frac{\&PartialD;h({\hat{x}}_k,u_k)}{\&PartialD;\hat{x}}f({\hat{x}}_k,u_k)=(\frac{{\&PartialD;U}_{\mathrm{OC},k}}{\&PartialD;\mathrm{SOC}}-1-1)\&times;f({\hat{x}}_k,u_k)---\left(21\right)$

$\begin{array}{c}{L}_f^2=\frac{{\&PartialD;L}_f^1}{\&PartialD;\hat{x}}f({\hat{x}}_k,u_k)=\&PartialD;\frac{(\frac{{\&PartialD;U}_{\mathrm{OC},k}}{\&PartialD;\mathrm{SOC}}-1-1)\&times;f({\hat{x}}_k,u_k)}{\&PartialD;\hat{x}}\&times;f({\hat{x}}_k,u_k)\\ =(\frac{{\&PartialD;}^2{U}_{\mathrm{OC},k}}{{\&PartialD;\mathrm{SOC}}^2}f({\hat{x}}_k,u_k)+(\frac{{\&PartialD;U}_{\mathrm{OC},k}}{\&PartialD;\mathrm{SOC}}-1-1)\frac{\&PartialD;f({\hat{x}}_k,u_k)}{\&PartialD;\hat{x}})\&times;f({\hat{x}}_k,u_k)\end{array}---\left(22\right)$

式中，是1×2维的矩阵，其行元素给出如下：

$S\left[{\hat{x}}_k\right]=\left\{L_{g_1}\right[L_f^1\left(h({\hat{x}}_k,u_k)\right)],L_{g_2}[L_f^1\left(h({\hat{x}}_k,u_k)\right)\left]\right\}---\left(23\right)$

式中，g₁和个g₂分别是G(k)中的第1和第2列元素；是先沿后沿g₁的李导数；是先沿后沿g₂的李导数。

根据上述公式可得到：

$S\left[{\hat{x}}_k\right]=\frac{{\&PartialD;L}_f^1}{\&PartialD;\hat{x}}G\left(k\right)=(\frac{{\&PartialD;}^2{U}_{\mathrm{OC},k}}{{\&PartialD;\mathrm{SOC}}^2}f({\hat{x}}_k,u_k)+(\frac{{\&PartialD;U}_{\mathrm{OC},k}}{\&PartialD;\mathrm{SOC}}-1-1)\frac{\&PartialD;f({\hat{x}}_k,u_k)}{\&PartialD;\hat{x}})*G\left(k\right)---\left(24\right)$

5、基于NPEKF的SOC估计

一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法，如图2所示，它的步骤总结为：

第一步，建立电池数学模型。

以二阶RC电池模型为例，选取RC网络中电容的电压和SOC作为系统状态变量，系统的状态方程为：

$\left[\begin{array}{c}{\mathrm{SOC}}_{k+1}^{-}\\ U_{1,k+1}^{-}\\ U_{2,k+1}^{-}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{esp}(-\mathrm{\&Delta;t}/R_1{C}_1)& 0\\ 0& 0& \mathrm{esp}(-\mathrm{\&Delta;t}/R_2{C}_2)\end{array}\right]**\left[\begin{array}{c}{\mathrm{SOC}}_k\\ U_{1,k}\\ U_{2,k}\end{array}\right]+\left[\begin{array}{c}-\mathrm{\&eta;}*\mathrm{\&Delta;t}/C\\ R_1*(1-\exp (-\mathrm{\&Delta;t}/R_1{C}_1))\\ R_2*(1-\exp (-\mathrm{\&Delta;t}/R_2{C}_2))\end{array}\right]*i_k---\left(25\right)$

上述方程可记为：

${\hat{x}}_{k+1}^{-}=A_k{x}_k+B_k{u}_k=f(x_k,u_k)---\left(26\right)$

其中，状态变量为：

$x_k=\left[\begin{array}{c}{\mathrm{SOC}}_k\\ U_{1,k}\\ U_{2,k}\end{array}\right]---\left(27\right)$

输入变量为：

u_k＝i_k   (28)

系统矩阵为：

$A_k=\frac{\&PartialD;f(x_k,u_k)}{\&PartialD;x}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \exp (-\mathrm{\&Delta;t}/R_1{C}_1)& 0\\ 0& 0& \exp (-\mathrm{\&Delta;t}/R_2{C}_2)\end{array}\right]---\left(29\right)$

输入控制矩阵为：

$B_k=\frac{\&PartialD;f(x_k,u_k)}{\&PartialD;u}=\left[\begin{array}{c}-\mathrm{\&eta;}*\mathrm{\&Delta;t}/C\\ R_1*(1-\exp (-\mathrm{\&Delta;t}/R_1{C}_1))\\ R_2*(1-\exp (-\mathrm{\&Delta;t}/R_2{C}_2))\end{array}\right]---\left(30\right)$

根据图1所示的动力电池模型，可得到电池系统的测量方程的离散形式为：

${\hat{U}}_{\mathrm{bat},k+1}={\left[\begin{array}{c}0\\ -1\\ -1\end{array}\right]}^T*\left[\begin{array}{c}{\mathrm{XOC}}_{k+1}^{-}\\ U_{1,k+1}^{-}\\ U_{2,k+1}^{-}\end{array}\right]-R_0*i_{k+1}+U_{\mathrm{OC},k}---\left(31\right)$

上述(27)式可记为：

${\hat{y}}_{k+1}=H_{k+1}{\hat{x}}_{k+1}^{-}+D_{k+1}{u}_{k+1}+v_{k+1}=h({\hat{x}}_{k+1}^{-},u_{k+1})+v_{k+1}---\left(32\right)$

式中，随机信号v_k+1表示过程观测噪声，假设为正态分布的白色噪声，且v_k+1～N(0,R_k+1)，R_k+1为第k+1步的观测噪声协方差。

测量矩阵为：

$H_{k+1}=\frac{\&PartialD;h({\hat{x}}_{k+1}^{-},)}{\&PartialD;x}=\left[\begin{array}{ccc}\frac{{\mathrm{dU}}_{\mathrm{OC}}}{\mathrm{dSOC}}{|}_{{\mathrm{SOC}}_{k+1}={\mathrm{SO}\hat{C}}_{k+1}}& -1& -1\end{array}\right]---\left(33\right)$

$D_{k+1}=\frac{\&PartialD;h({\hat{x}}_{k+1}^{-},u_{k+1})}{\&PartialD;u}=\left[R_0\right]---\left(34\right)$

输出变量为：

y_k+1＝U_bat,k+1   (35)

上式中，U_bat为电池的端电压；R₀为欧姆内阻；R₁为电化学极化内阻，C₁为电化学极化电容；R₂为浓差极化电阻，C₂为浓差极化电容；U₁，U₂分别是电容C₁和C₂的端电压；U_oc为电池的开路电压；i_k为k时刻的电池端电流。

第二步，NPEKF滤波器的初始化。

给定状态初始估计值加权矩阵W以及测量噪声方差矩阵R和误差协方差矩阵P₀。

状态估计值初始化：

$x_0=\left[\begin{array}{c}{\mathrm{SOC}}_0\\ U_{\mathrm{1,0}}\\ U_{\mathrm{2,0}}\end{array}\right]=\left[\begin{array}{c}{\mathrm{SOC}}_0\\ 0\\ 0\end{array}\right]---\left(36\right)$

式中，SOC₀为估计初值。

误差协方差初始化：

$P_0=\left[\begin{array}{ccc}10& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]---\left(37\right)$

过程激励噪声协方差：

$Q=\left[\begin{array}{ccc}0.05& 0& 0\\ 0& 0.1& 0\\ 0& 0& 0.1\end{array}\right]---\left(38\right)$

测量噪声协方差：

$R=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0.1& 0\\ 0& 0& 0.1\end{array}\right]---\left(39\right)$

选取G(k)为：

$G\left(k\right)=\left[\begin{array}{cc}0& 0\\ 1& 0\\ 0& 1\end{array}\right]---\left(40\right)$

结合先验知识选取加权矩阵W为：

$W=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]---\left(41\right)$

其中，W为模型误差加权矩阵，可根据经验事先给定。值得说明的是，模型误差加权矩阵W的取值对估计效果的影响很大，但没有固定的选取规则。通过仿真发现，若模型误差是加性的或者初始估计误差过大时，W的值不应过大；若模型误差为乘性的，W的值可以大一些；测量误差一般对W的取值不构成影响。

第三步，计算参量

设已获得了k时刻的状态估计值利用计算出中间参数矩阵 在接到k+1时刻的测量信息y(k+1)后，根据公式(26)预测[k，k+1]间隔内的模型误差项d(k)，k＝n，n为自然数。

Λ(Δt)是1×1维的对角阵，其值为：

Λ(Δt)＝Δt²/2   (42)

的计算公式如下：

$Z[\hat{x}\left(t\right),\mathrm{\&Delta;t}]=\mathrm{\&Delta;t}{\&times;L}_f^1+({\mathrm{\&Delta;t}}^2/2)\&times;L_f^2---\left(43\right)$

式中，是关于的m阶李导数，定义为：

$\left\{\begin{array}{c}{L}_f^m\left(h({\hat{x}}_k,u_k)\right)=h({\hat{x}}_k,u_k),m=0\\ L_f^m\left(h({\hat{x}}_k,u_k)\right)=\frac{{\&PartialD;L}_f^{m-1}\left(h({\hat{x}}_k,u_k)\right)}{\&PartialD;\hat{x}}f({\hat{x}}_k,u_k),m\&GreaterEqual;1\end{array}\right.---\left(44\right)$

根据上述定义可以得到：

$L_f^1=\frac{\&PartialD;h({\hat{x}}_k,u_k)}{\&PartialD;\hat{x}}f\left(\hat{x}\right)=\frac{\&PartialD;h({\hat{x}}_k,u_k)}{\&PartialD;\hat{x}}f({\hat{x}}_k,u_k)=(\frac{{\&PartialD;U}_{\mathrm{OC},k}}{\&PartialD;\mathrm{SOC}}-1-1)\&times;f({\hat{x}}_k,u_k)---\left(45\right)$

$\begin{array}{c}{L}_f^2=\frac{{\&PartialD;L}_f^1}{\&PartialD;\hat{x}}f({\hat{x}}_k,u_k)=\&PartialD;\frac{(\frac{{\&PartialD;U}_{\mathrm{OC},k}}{\&PartialD;\mathrm{SOC}}-1-1)\&times;f({\hat{x}}_k,u_k)}{\&PartialD;\hat{x}}\&times;f({\hat{x}}_k,u_k)\\ =(\frac{{\&PartialD;}^2{U}_{\mathrm{OC},k}}{{\&PartialD;\mathrm{SOC}}^2}f({\hat{x}}_k,u_k)+(\frac{{\&PartialD;U}_{\mathrm{OC},k}}{\&PartialD;\mathrm{SOC}}-1-1)\frac{\&PartialD;f({\hat{x}}_k,u_k)}{\&PartialD;\hat{x}})\&times;f({\hat{x}}_k,u_k)\end{array}---\left(46\right)$

式中，是1×2维的矩阵，其行元素给出如下：

$S\left[{\hat{x}}_k\right]=\left\{L_{g_1}\right[L_f^1\left(h({\hat{x}}_k,u_k)\right)],L_{g_2}[L_f^1\left(h({\hat{x}}_k,u_k)\right)\left]\right\}---\left(47\right)$

式中，g₁和个g₂分别是G(k)中的第1和第2列元素；是先沿后沿g₁的李导数；是先沿后沿g₂的李导数。

根据上述公式可得到：

$S\left[{\hat{x}}_k\right]=\frac{{\&PartialD;L}_f^1}{\&PartialD;\hat{x}}G\left(k\right)=(\frac{{\&PartialD;}^2{U}_{\mathrm{OC},k}}{{\&PartialD;\mathrm{SOC}}^2}f({\hat{x}}_k,u_k)+(\frac{{\&PartialD;U}_{\mathrm{OC},k}}{\&PartialD;\mathrm{SOC}}-1-1)\frac{\&PartialD;f({\hat{x}}_k,u_k)}{\&PartialD;\hat{x}})*G\left(k\right)---\left(48\right)$

第四步，预测模型误差d(k)。

$\begin{array}{c}d\left(k\right)={\{{\left[\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left({\hat{x}}_k\right)\right]}^T{R}^{-1}\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left({\hat{x}}_k\right)+W\}}^{-1}\\ \&times;{\left[\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left({\hat{x}}_k\right)\right]}^T{R}^{-1}[Z({\hat{x}}_k,\mathrm{\&Delta;t})+{\hat{y}}_k-y_{k+1}]\end{array}---\left(49\right)$

把d(k)的预测值代入状态估计方程，将状态估计值传播到k+1时刻，即得到

${\hat{x}}_{k+1}^{-}=f({\hat{x}}_k,u_k)+G\left(k\right)d\left(k\right)---\left(50\right)$

式中，d(k)为2×1维的模型误差矩阵；G(k)为3×2维的模型误差分配矩阵。-代表先验，^代表估计，为在已知第k+1步以前状态情况下第k+1步的先验状态估计；为已知测量真值y_k+1时，第k+1步的后验状态估计。

第六步，误差协方差矩阵预测

${\hat{P}}_{k+1}^{-}=A_{k+1}{\hat{P}}_k{A_{k+1}}^T+w_{k+1}{Q}_k{w}_{k+1}^T---\left(51\right)$

式中，Q_k+1为第k+1步的过程噪声协方差，随机信号w_k+1表示过程激励噪声，假设为正态分布的白色噪声：w_k～N(0,Q_k)。

第七步，计算增益矩阵K_k+1：

$K_{k+1}={\hat{P}}_{k+1}^{-}{H}_{k+1}^T{[H_{k+1}{\hat{P}}_{k+1}^{-}{H}_{k+1}^t+v_{k+1}{R}_{k+1}{v}_{k+1}^T]}^{-1}---\left(52\right)$

式中，R_k+1为第k+1步的观测噪声协方差，随机信号v_k+1表示过程观测噪声，假设为正态分布的白色噪声，且v_k+1～N(0,R_k+1)。

第六步，状态修正

${\hat{y}}_{k+1}=h({\hat{x}}_{k+1}^{-},u_k)+v_{k+1}---\left(53\right)$

${\hat{x}}_{k+1}={\hat{x}}_{k+1}^{-}+K_{k+1}[y_{k+1}-{\hat{y}}_{k+1}]---\left(54\right)$

式中，y_k+1为第k+1步的测量真值，本发明中具体指代电池的端电压U_bat,k+1；第k+1步的测量估计值，即电池的端电压的估计值；u_k+1为第k+1步的系统输入值，本发明中具体指代电池的端电流i_bat,k+1。

第九步，误差协方差矩阵修正P_k+1：

${\hat{P}}_{k+1}=(1-K_{k+1}{H}_{k+1}){\hat{P}}_{k+1}^{-}---\left(55\right)$

第十步，根据方程(54)得到第k+1步的估计结果，并与SOC真值进行比较，验证NPEKF算法的有效性，k＝n+1，并转到第三步。

图6给出了在UDDS循环工况下基于NPEKF的SOC估计结果。NPEKF的初始参数设置为：SOC真值初始值SOC_t0＝0.9，SOC估计初始值SOC_e0＝0.8。从实验结果可以看出，SOC估计在初始SOC有较大误差的情况下能迅速向实际值逼近，之后SOC估计值在SOC真值附近波动，最大估计误差小于4％，证明了本发明提出的NPEKF算法的有效性。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (7)

1.一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法，其特征是，包括步骤：

1)将整个方法所需要的时间平均分成N个时间段，每一个时间段代表一步，即第k个时间段为第k步；其中k≤N，k与N均为正整数；

2)根据基尔霍夫电压和电流定理，建立电池系统的模型；以电池的电流为输入变量，以SOC值与支路电压为状态变量，得到电池系统的先验状态估计方程；以电池两端的电压为输出变量，以电池的电流为输入变量，以SOC值与支路电压为状态变量，得到电池系统的测量方程；根据电池系统的模型的测量方程，得到电池系统模型的误差矩阵d(k)与电池系统的模型误差分配矩阵G(k)；

3)对步骤2)所得的先验状态估计方程进行补偿：利用电池系统的模型误差矩阵d(k)与电池系统的模型误差分配矩阵G(k)对步骤2)所得的先验状态估计方程进行补偿，得到精准的先验状态估计方程，k＝n，n为正整数，n≤N；

4)求解第k+1步电池系统的后验状态估计方程，得到SOC值；

5)根据步骤4)得到的第k+1步电池系统的后验状态估计结果，与电池的SOC真值进行比较，验证非线性预测卡尔曼滤波NPEKF算法的有效性，k＝n+1，并转到步骤3)，直至第N步执行完。

2.如权利要求1所述的一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法，其特征是，所述步骤3)中对先验状态估计方程进行补偿的具体步骤为：

31)对电池系统的模型误差分配矩阵G(k)进行初始化；

32)求解电池系统的模型误差矩阵d(k)；

33)对步骤2)所得的先验状态估计方程进行补偿：

${\hat{x}}_{k+1}^{-}=f({\hat{x}}_k,u_k)+G\left(k\right)d\left(k\right);$

其中为第k+1步的先验状态估计，为关于的函数，x_k为状态变量；u_k为输入变量。

3.如权利要求2所述的一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法，其特征是，所述步骤31)还包括对NPEKF滤波器的初始化：

给定状态初始估计值给定电池系统的模型误差加权矩阵W、过程激励噪声协方差矩阵Q、测量噪声方差矩阵R和误差协方差矩阵P₀的初始值。

4.如权利要求书1-2任一所述的一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法，其特征是，所述步骤3)中电池系统的模型误差矩阵d(k)的具体求解过程为：

321)求解参数Λ(Δt)：

Λ(Δt)是1×1维的对角阵，其值为：

Λ(Δt)＝Δt²/2；

322)求解参数

的计算公式如下：

$Z[\hat{x}\left(t\right),\mathrm{\&Delta;t}]=\mathrm{\&Delta;t}{\&times;L}_f^1+({\mathrm{\&Delta;t}}^2/2)\&times;L_f^2;$

式中，分别是关于的第1阶和第2阶李导数，为输出变量估计关于的函数，为第k+1步的状态估计关于的函数，输入变量u_k＝i_k，i_k为第k个时刻点的电池端电流；

323)求解参数

是1×2维的矩阵，其行元素如下：

$S\left[{\hat{x}}_k\right]=\left\{L_{g_1}\right[L_f^1\left(h({\hat{x}}_k,u_k)\right)],L_{g_2}[L_f^1\left(h({\hat{x}}_k,u_k)\right)\left]\right\}$

式中，g₁和个g₂分别是G(k)中的第1和第2列元素；是先沿后沿g₁的李导数；是先沿后沿g₂的李导数；

根据上述公式可得到：

$S\left[{\hat{x}}_k\right]=\frac{{\&PartialD;L}_f^1}{{\&PartialD;\hat{x}}_k}G\left(k\right)=(\frac{{\&PartialD;}^2{U}_{\mathrm{OC},k}}{{\&PartialD;}^2\mathrm{SOC}}f({\hat{x}}_k,u_k)+(\frac{{\&PartialD;U}_{\mathrm{OC},k}}{\&PartialD;\mathrm{SOC}}-1-1)\frac{\&PartialD;f({\hat{x}}_k,u_k)}{{\&PartialD;\hat{x}}_k})*G\left(k\right)$

式中，U_oc,k为第k步电池开路电压，SOC为电池的电荷状态；

324)求解模型误差矩阵d(k)：

$\begin{array}{c}d\left(k\right)=-{\{{\left[\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left({\hat{x}}_k\right)\right]}^T{R}^{-1}\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left({\hat{x}}_k\right)+W\}}^{-1}\\ \&times;{\left[\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)S\left({\hat{x}}_k\right)\right]}^T{R}^{-1}\left[Z({\hat{x}}_k,\mathrm{\&Delta;t}){+{\hat{y}}_k-y_{k+1}}_{}\right]\end{array}$

其中，R为测量噪声协方差矩阵，W为加权矩阵，Δt＝t_k+1-t_k为采样时间间隔，为第k步的状态估计，为第k步的y_k的估计，y_k+1为第k+1步的测量真值，y_k+1＝U_bat,k+1，U_bat,k+1第k+1步的电池端电压。

5.如权利要求1所述的一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法，其特征是，所述步骤4)求解电池系统的后验状态估计方程的具体步骤为：

41)计算第k+1步的增益矩阵K_k+1；

42)求解后验估计

${\hat{y}}_{k+1}=h({\hat{x}}_{k+1}^{-},u_k)+v_{k+1};$

${\hat{x}}_{k+1}={\hat{x}}_{k+1}^{-}+K_{k+1}[y_{k+1}-{\hat{y}}_{k+1}];$

其中，为第k+1步测量值估计值，y_k+1为第k+1步的测量真值，具体指代为第k+1步电池的端电压U_bat,k+1，u_k+1＝i_k+1，i_k+1为第k+1步的电池电流，为第k+1步状态的先验估计。

6.如权利要求5所述的一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法，其特征是，所述步骤41)中第k+1步的增益矩阵K_k+1的具体求解过程为：

411)计算第k+1步误差协方差矩阵的先验估计

${\hat{P}}_{k+1}^{-}=A_{k+1}{\hat{p}}_k{A_{k+1}}^T+w_{k+1}{Q}_k{w}_{k+1}^T;$

其中，A_k+1为第k+1步的系统矩阵，为第k步的误差协方差矩阵的后验估计，Q_k为第k步的过程激励噪声协方差矩阵，w_k+1表示过程激励噪声，为正态分布的白色噪声：w_k+1～N(0,R_k+1)；

412)计算第k+1步的增益矩阵K_k+1：

$K_{k+1}={\hat{P}}_{k+1}^{-}{H}_{k+1}^T{[H_{k+1}{\hat{P}}_{k+1}^{-}{H}_{k+1}^T+v_{k+1}{R}_{k+1}{v}_{k+1}^T]}^{-1};$

其中，为第k+1步的误差协方差矩阵的先验估计，H_k+1为第k+1步的测量矩阵，v_k+1为过程观测噪声，为正态分布v_k+1～N(0,R_k+1)。

7.如权利要求6所述的一种基于非线性预测扩展卡尔曼滤波的电池SOC估计方法，其特征是，所述步骤411)第k步的误差协方差矩阵的后验估计的具体求解为：

${\hat{P}}_k=(1-K_k{H}_k){\hat{P}}_{k+1}^{-};$

其中K_k为第k步的增益矩阵，H_k为第k步的测量矩阵，为第k步的误差协方差矩阵的先验估计。