CN103116136B - 基于有限差分扩展卡尔曼算法的锂电池荷电状态估计方法 - Google Patents

Abstract

本发明公开了一种锂电池荷电状态的估计方法，该方法包括以下步骤：步骤一、给定初始值，对各协方差进行Cholesky分解；步骤二、状态的一步预测；步骤三、协方差的一步预测；步骤四、滤波增益；步骤五、更新状态的最优值；步骤六、滤波协方差更新。与现有技术相比，本发明的精度高于泰勒级数的一阶展开，而且，充分利用了由模型线性化产生的有效误差信息，对模型参数变动具有较强的鲁棒性。

Description

基于有限差分扩展卡尔曼算法的锂电池荷电状态估计方法

技术领域

本发明锂电池荷电状态预测技术，特别是涉及一种在锂电池实际应用过程中对其荷电状态进行估计的方法。

背景技术

电池作为电源已在通信、电力系统、军事装备、电动汽车等领域得到了广泛地应用。随着环保观念的日渐深入人心，越来越多的系统开始采用电池作为主要动力供给。在这些系统中，动力电池工作状态的好坏直接关系到整个系统的运行可靠性。为确保动力电池组的性能良好，延长电池组的使用寿命，须及时、准确地了解电池的运行状态，对电池进行合理有效地管理和控制。而电池荷电状态的精确估算是电池管理系统中最核心的技术。目前，对单体电池荷电状态估计的主要方法有开路电压法、内阻法、安时计量法、神经网络、卡尔曼滤波法和扩展卡尔曼滤波法等。

其中，开路电压法进行电池荷电状态估计时电池必须静置较长时间以达到稳定状态，而且只适用于电池电流非剧烈变化状态下的荷电状态估计，不能满足在线检测要求。内阻法需要精确测量电池的内阻，由于一般电池内阻在毫欧级，因此使用该方法会对测量仪器的要求非常高，难以在实际中加以应用。安时计量法易受到电流测量精度的影响，在高温或电流波动剧烈情况下，精度很差。神经网络法需要大量的样本数据进行训练才可以得到较好的精度。卡尔曼滤波法对电池模型要求较高，并且只适用于线性方程，这一问题的解决办法之一是采用扩展卡尔曼滤波法来进行荷电状态的估计。但是，扩展卡尔曼滤波法为了求取估计误差方差通过非线性函数的传播，需将非线性函数用泰勒级数展开并线性化，一方面增加了计算复杂度，另一方面，在该过程中忽略了高阶项，必定会引起误差。

发明内容

为克服现有技术的不足，本发明提出一种基于有限差分扩展卡尔曼算法的锂电池荷电状态估计方法，在锂电池的电池数学模型的基础上，利用有限差分卡尔曼滤波法进行剩余电量的估计，在对非线性方程进行线性化处理上，采用多项式近似技术和一阶中心差分法计算非线性函数的偏导数。

本发明提出了一种基于有限差分扩展卡尔曼算法的锂电池荷电状态估计方法，该方法包括以下步骤：

一种基于有限差分扩展卡尔曼算法的锂电池荷电状态估计方法，其特征在于，该方法包括以下步骤：

通过建立锂电池电量的数学模型，得到电池系统荷电状态方程x及荷电状态观测方程y如下：

x_k+1＝f(x_k,u_k，w_k)

y_k＝g(x_k,u_k，v_k)

其中，x_k∈Rⁿ和y_k∈R^m分别为k时刻的系统n维状态向量和m维观测向量，f:Rⁿ→Rⁿ和g:Rⁿ→R^m分别为系统状态函数和观测函数，w_k和v_k是互不相关的高斯白噪声，且假设噪声有如下统计学特性：

$\left\{\begin{array}{c}E\left[w_k\right]=w_k,\mathrm{Cov}[w_k,w_j^T]=Q_k{\mathrm{\δ}}_{\mathrm{kj}},\\ E\left[v_k\right]=v_k,\mathrm{Cov}[v_k,v_j^T]=R_k{\mathrm{\δ}}_{\mathrm{kj}},\\ \mathrm{Cov}[w_k,v_k]=0,\end{array}\right.$

其中：Q_k，R_k都为正定对称阵；δ_kj为Kronecker-δ函数；

步骤一、给定初始值，初始化x₀＝E[x₀],P₀＝E[(x₀-x₀)(x₀-x₀)^T]，对Q_k(过程噪声方差)，R_k(观测噪声方差)，（验前状态方差），P_k（验后状态方差）进行Cholesky分解，有

$Q=S_w{S}_w^T,$ $R=S_v{S}_v^T,$ $\stackrel{\‾}{P}={\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T,$ $P=S_x{S}_x^T,$

利用一阶中心差分计算非线性函数偏导数，即

F_x(k)＝(f(x_k+Δx_k,u_k,w_k)-f(x_k-Δx_k,u_k,w_k))/2Δx_k，

令Δx_k＝hS_x（h为步长调节系数），则

F_x(k)S_x＝S_xx＝{(f_i(x_k+hS_x，j,u_k,w_k)-f_i(x_k-hS_x，j,u_k,w_k))/2h，

其中，S_x，j为S_x的第j列。根据上式同样可导出下列各式。即

F_w(k)S_w＝S_xw＝{(f_i(x_k,u_k,w_k+hS_x,j)-f_i(x_k,u_k,w_k-hS_x,j))/2h}

$G_x\left(k\right){\stackrel{\‾}{S}}_x=S_{y\stackrel{\‾}{x}}=\{(g_i({\stackrel{\‾}{x}}_k+h{\stackrel{\‾}{S}}_{x,j},u_k,v_k)-g_i({\stackrel{\‾}{x}}_k-h{\stackrel{\‾}{S}}_{x,j},u_k,v_k))/2h\}$

$G_v\left(k\right)S_v=S_{\mathrm{yv}}=\{(g_i({\stackrel{\‾}{x}}_k,u_k,v_k+h{\stackrel{\‾}{S}}_{x,j})-g_i({\stackrel{\‾}{x}}_k{,u}_k,v_k-h{\stackrel{\‾}{S}}_{x,j}))/2h\},$

步骤二、得到基于锂电池系统上一状态预测出当前状态方程

${\stackrel{\‾}{x}}_{k+1}\≈f(x_k,u_k,w_k),$

步骤三、得到状态协方差的预测方程

$\stackrel{\‾}{P}\left(k\right)=F_x\left(k\right)P\left(k\right)F_x^T\left(k\right)+F_w\left(k\right)Q\left(k\right)F_w^T\left(k\right)$

$=F_x\left(k\right)S_x{S}_x^T{F}_x^T\left(k\right)+F_w\left(k\right)S_w{S}_w^T{F}_w^T\left(k\right)$

$=S_{\mathrm{xx}}{S}_{\mathrm{xx}}^T+S_{\mathrm{xw}}{S}_{\mathrm{xw}}^T,$

其中：是当前状态的验前状态协方差；

F_x(k)是利用一阶中心差分计算状态方程的偏导数；

P(k)是上一状态的验后状态协方差；

F_w(k)是利用一阶中心差分计算状态噪声的偏导数；

Q(k)是状态噪声协方差。

步骤四、得到滤波增益方程

$K_k=\stackrel{\‾}{P}\left(k\right)G_x^T\left(k\right){[G_x\left(k\right)\stackrel{\‾}{P}\left(k\right)G_x^T\left(k\right)+G_v\left(k\right)R\left(k\right)G_v^T\left(k\right)]}^{-1}$

$={\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T{\left(S_{y\stackrel{\‾}{x}}{\stackrel{\‾}{S}}_x^{-1}\right)}^T{[S_{y\stackrel{\‾}{x}}{S}_{y\stackrel{\‾}{x}}^T+S_{\mathrm{yv}}{S}_{\mathrm{yv}}^T]}^{-1}$

$={\stackrel{\‾}{S}}_x{S}_{y\stackrel{\‾}{x}}^T{[S_{y\stackrel{\‾}{x}}{S}_{y\stackrel{\‾}{x}}^T+S_{\mathrm{yv}}{S}_{\mathrm{yv}}^T]}^{-1},$

其中：K_k是卡尔曼滤波增益；

G_x(k)是利用一阶中心差分计算观测方程的偏导数；

G_v(k)是利用一阶中心差分计算观测噪声的偏导数；

R(k)是观测噪声协方差。

步骤五、根据当前状态的预测值，再收集当前状态的测量值。结合预测值和测量值，我们可以得到当前状态的最优值。得到更新状态的最优值

${\stackrel{\‾}{y}}_k=g({\stackrel{\‾}{x}}_k,u_k,v_k)$

$x_{k+1}={\stackrel{\‾}{x}}_{k+1}+K_k[y_k-{\stackrel{\‾}{y}}_k],$

其中，是当前状态的预测输出；

x_k+1是当前状态的最优值；

y_k是当前状态实际的测量值。

步骤六、得到滤波协方差更新

$P\left(k\right)=\stackrel{\‾}{P}\left(k\right)-K_k{G}_v\left(k\right)\stackrel{\‾}{P}\left(k\right)={\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T-K_k{G}_v\left(k\right){\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T$

$={\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T-{\stackrel{\‾}{S}}_x{S}_{y\stackrel{\‾}{x}}^T-K_k^T-K_k{S}_{y\stackrel{\‾}{x}}{\stackrel{\‾}{S}}_x^T+{\stackrel{\‾}{S}}_x{S}_{y\stackrel{\‾}{x}}^T{K}_k^T$

$={\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T-{\stackrel{\‾}{S}}_x{S}_{y\stackrel{\‾}{x}}^T{K}_k^T-K_k{S}_{y\stackrel{\‾}{x}}{\stackrel{\‾}{S}}_x^T+K_k{S}_{y\stackrel{\‾}{x}}{S}_{y\stackrel{\‾}{x}}^T{K}_k^T+K_k{S}_{\mathrm{yv}}{S}_{\mathrm{yv}}^T{K}_k^T$

$=[{\stackrel{\‾}{S}}_x-K_k{S}_{y\stackrel{\‾}{x}}{K}_k{S}_{\mathrm{yv}}]{[{\stackrel{\‾}{S}}_x-K_k{S}_{y\stackrel{-}{x}}{K}_k{S}_{\mathrm{yv}}]}^T,$

其中，P(k)是当前状态的验后状态协方差。

在实际电池运行过程中，实时采集锂电池的电压、电流信息，通过现场可编程门阵列/数字信号处理器运行电池荷电状态估计算法，将估算结果显示在液晶屏上。与此同时，对电池组充放电电路进行控制，保证电池运行状态的安全性和稳定性。

与现有技术相比，本方法的有益效果为，其精度高于泰勒级数的一阶展开，而且，充分利用了由模型线性化产生的有效误差信息，对模型参数变动具有较强的鲁棒性。

附图说明

图1是本发明的锂电池荷电状态的估计方法流程图；

图2是本发明的锂电池荷电状态的估计方法的具体实施系统的结构示意图。

具体实施方式

下面将结合附图对本发明的具体实施方式进行详细描述，这些实施方式若存在示例性的内容，不应解释成对本发明的限制。

锂电池的电池模型是个非线性系统，其中电池的开路电压和电池荷电状态呈非线性关系。扩展卡尔曼滤波算法在卡尔曼滤波算法的基础上，对非线性方程进行线性化处理，可以适用于非线性系统。当模型参数与过程参数精确匹配或基本匹配时，只要初始值选择适当，滤波过程可渐进收敛并得到状态的近似无偏估计值。然而，当模型参数不匹配或存在较大偏差时，扩展卡尔曼滤波算法的估计精度会大大下降甚至发散。

如图2所示，为本方法的实现系统举例，该实现系统包括现场可编程门阵列/数字信号处理器控制芯片、液晶显示模块、模数转换以及采集电路。在核心控制芯片中，实现电池荷电状态估计算法、充放电电路的控制、液晶显示时序控制、以及对电池各参数量的采集控制。锂电池的电压、电流、温度通过采集电路以及模数转换电路进入控制芯片，从而运行电池荷电状态估计算法，得到各单体电池的荷电状态值，并将电池实时状态显示在液晶屏上。从而实现了电池的荷电状态估计和实时监控。与此同时，通过实际要求控制充放电电路，优化电池组系统性能，保证系统工作的安全性和可靠性。

Claims (1)

1.一种基于有限差分扩展卡尔曼算法的锂电池荷电状态估计方法，其特征在于，该方法包括以下步骤：

通过建立锂电池电量的数学模型，得到电池系统荷电状态方程x及荷电状态观测方程y如下：

x_k+1＝f(x_k,u_k，w_k)

y_k＝g(x_k,u_k，v_k)，

其中，x_k∈Rⁿ和y_k∈R^m分别为k时刻的系统n维状态向量和m维观测向量，f:Rⁿ→Rⁿ和g:Rⁿ→R^m分别为系统状态函数和观测函数，w_k和v_k是互不相关的高斯白噪声，且假设噪声有如下统计学特性：

$\left\{\begin{array}{c}E\left[w_k\right]=w_k,\mathrm{Cov}[w_k,w_j^T]=Q_k{\mathrm{\δ}}_{\mathrm{kj}},\\ E\left[v_k\right]=v_k,\mathrm{Cov}[v_k,v_j^T]=R_k{\mathrm{\δ}}_{\mathrm{kj}},\\ \mathrm{Cov}[w_k,v_k]=0,\end{array}\right.$

其中：Q_k，R_k都为正定对称阵；δ_kj为Kronecker-δ函数；

步骤一、给定初始值，初始化对过程噪声方差Q_k、观测噪声方差R_k、验前状态方差验后状态方差进行Cholesky分解，有

$Q=S_w{S}_w^T,$ $R=S_v{S}_v^T,$ $\stackrel{\‾}{P}={\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T,$ $\hat{P}={\hat{S}}_x{\hat{S}}_x^T,$

S_w为过程噪声方差矩阵Q_kCholesky分解下的三角矩阵，S_v为观测噪声方差矩阵R_kCholesky分解的下三角矩阵，为验前状态方差矩阵Cholesky分解的下三角矩阵，为验后状态方差矩阵Cholesky分解的下三角矩阵；

利用一阶中心差分计算非线性函数偏导数，即

$F_x\left(k\right)=(f({\hat{x}}_k+\mathrm{\Δ}{\hat{x}}_k,u_k,w_k)-f({\hat{x}}_k-\mathrm{\Δ}{\hat{x}}_k,u_k,w_k))/2\mathrm{\Δ}{\hat{x}}_k,$

令h为步长调节系数，则

$F_x\left(k\right){\hat{S}}_x=S_{x\hat{x}}=\{(f_i({\hat{x}}_k+h{\hat{S}}_{x,j},u_k,w_k)-f_i({\hat{x}}_k-h{\hat{S}}_{x,j},u_k,w_k))/2h,$

其中，为的第j列，根据上式同样可导出下列各式，即

$F_w\left(k\right)S_w=S_{\mathrm{xw}}=\{(f_i({\hat{x}}_k,u_k,w_k+h{\hat{S}}_{x,j})-f_i({\hat{x}}_k,u_k,w_k-h{\hat{S}}_{x,j}))/2h\}$

$G_x\left(k\right){\stackrel{\‾}{S}}_x=S_{y\stackrel{\‾}{x}}=\{(g_i({\stackrel{\‾}{x}}_k+h{\stackrel{\‾}{S}}_{x,j},u_k,v_k)-g_i({\stackrel{\‾}{x}}_k-h{\stackrel{\‾}{S}}_{x,j},u_k,v_k))/2h\}$

$G_v\left(k\right)S_v=S_{\mathrm{yv}}=\{(g_i({\stackrel{\‾}{x}}_k,u_k,v_k+h{\stackrel{\‾}{S}}_{x,j})-g_i({\stackrel{\‾}{x}}_k,u_k,v_k-h{\stackrel{\‾}{S}}_{x,j}))/2h\},$

G_x(k)是利用一阶中心差分计算观测方程的偏导数；G_v(k)是利用一阶中心差分计算观测噪声的偏导数；

步骤二、得到基于锂电池系统上一状态预测出当前状态方程

${\stackrel{\‾}{x}}_{k+1}\≈f({\hat{x}}_k,u_k,w_k),$

步骤三、得到状态协方差的预测方程

$\begin{array}{c}\stackrel{\‾}{P}\left(k\right)=F_x\left(k\right)\hat{P}\left(k\right)F_x^T\left(k\right)+F_w\left(k\right)Q\left(k\right)F_w^T\left(k\right)\\ =F_x\left(k\right){\hat{S}}_x{\hat{S}}_x^T{F}_x^T\left(k\right)+F_w\left(k\right)S_w{S}_w^T{F}_w^T\left(k\right)\\ =S_{x\hat{x}}{S}_{x\hat{x}}^T+S_{\mathrm{xw}}{S}_{\mathrm{xw}}^T\end{array},$

其中：是当前状态的验前状态协方差；

F_x(k)是利用一阶中心差分计算状态方程的偏导数；

是上一状态的验后状态协方差；

F_w(k)是利用一阶中心差分计算状态噪声的偏导数；

Q(k)是状态噪声协方差；

步骤四、得到滤波增益方程

$\begin{array}{c}{K}_k=\stackrel{\‾}{P}\left(k\right)G_x^T\left(k\right){[G_x\left(k\right)\stackrel{\‾}{P}\left(k\right)G_x^T\left(k\right)+G_v\left(k\right)R\left(k\right)G_v^T\left(k\right)]}^{-1}\\ ={\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T{\left(S_{y\stackrel{\‾}{x}}{\stackrel{\‾}{S}}_x^{-1}\right)}^T{[S_{y\stackrel{\‾}{x}}{S}_{y\stackrel{\‾}{x}}^T+S_{\mathrm{yv}}{S}_{\mathrm{yv}}^T]}^{-1}\\ ={\stackrel{\‾}{S}}_x{S}_{y\stackrel{\‾}{x}}^T{[S_{y\stackrel{\‾}{x}}{S}_{y\stackrel{\‾}{x}}^T+S_{\mathrm{yv}}{S}_{\mathrm{yv}}^T]}^{-1}\end{array},$

其中：K_k是卡尔曼滤波增益；

G_x(k)是利用一阶中心差分计算观测方程的偏导数；

G_v(k)是利用一阶中心差分计算观测噪声的偏导数；

R(k)是观测噪声协方差；

步骤五、根据当前状态的预测值，再收集当前状态的测量值，结合预测值和测量值，我们可以得到当前状态的最优值，得到更新状态的最优值

${\stackrel{\‾}{y}}_k=g({\stackrel{\‾}{x}}_k,u_k,v_k)$

${\hat{x}}_{k+1}={\stackrel{\‾}{x}}_{k+1}+K_k[y_k-{\stackrel{\‾}{y}}_k],$

其中，是当前状态的预测输出；

是当前状态的最优值；

y_k是当前状态实际的测量值；

步骤六、得到滤波协方差更新

$\begin{array}{c}\hat{P}\left(k\right)=\stackrel{\‾}{P}\left(k\right)-K_k{G}_v\left(k\right)\stackrel{\‾}{P}\left(k\right)={\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T-K_k{G}_v\left(k\right){\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T\\ ={\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T-{\stackrel{\‾}{S}}_x{S}_{y\stackrel{\‾}{x}}^T{K}_k^T-K_k{S}_{y\stackrel{\‾}{x}}{\stackrel{\‾}{S}}_x^T+{\stackrel{\‾}{S}}_x{S}_{y\stackrel{\‾}{x}}^T{K}_k^T\\ ={\stackrel{\‾}{S}}_x{\stackrel{\‾}{S}}_x^T-{\stackrel{\‾}{S}}_x{S}_{y\stackrel{\‾}{x}}^T{K}_k^T-K_k{S}_{y\stackrel{\‾}{x}}{\stackrel{\‾}{S}}_x^T+K_k{S}_{y\stackrel{\‾}{x}}{S_{y\stackrel{\‾}{x}}^{T}K}_k^T+K_k{S}_{\mathrm{yv}}{S}_{\mathrm{yv}}^T{K}_k^T\\ =[{\stackrel{\‾}{S}}_x-K_k{S}_{y\stackrel{\‾}{x}}{K}_k{S}_{\mathrm{yv}}]{[{\stackrel{\‾}{S}}_x-K_k{S}_{y\stackrel{\‾}{x}}{K}_k{S}_{\mathrm{yv}}]}^T\end{array},$

其中，是当前状态的验后状态协方差。