锂离子电池建模与状态估计研究

Lithium-Ion Battery Modeling and State Estimation Research

摘要: 通过混合储能技术，提升配电网的动态增容能力与灵活调节能力，增强电力系统源网荷储的协调控制水平。首先对储能电池进行数学建模，再基于数学模型，用数学公式精准模拟真实电池特性，进而实现对电池系统的仿真和算法开发以及电池系统的SOC估计和健康估计。

Abstract: Through hybrid energy storage technology, the dynamic capacity increases, the capacity and flexible adjustment ability of the distribution network are improved, and the coordination and control level of source-grid-load-storage of the power system is enhanced. Firstly, the energy storage battery is mathematically modeled, and then based on the mathematical model, the real battery characteristics are accurately simulated with mathematical formulas, and then the simulation and algorithm development of the battery system and the SOC estimation and health estimation of the battery system are realized.

1. 引言

以锂离子电池为核心，运用数学建模与等效电路建模相融合的方法，构建单体锂电池的电化学储能电池数字模型。研究过程采用数字镜像建模、模型理论分析以及实际测试数据验证相结合的方式，在获得时变参数耦合的储能电池模型后，依据状态估计理论，对储能电池荷电状态(state of charge，SOC)估计算法展开深入探究，旨在实现对锂电池荷电状态的科学、精准估计。

2. 锂电池的建模

锂离子电池的一阶Thevenin等效电路模型的数学建模如下图1所示，模型的输出电压公式可以表示为：

$U_{\text{oc}}=U_{\text{L}}+R_{\text{o}}{I}_{\text{L}}+U_{C_{\text{th}}}$

$I_{\text{L}}=\frac{U_{C_{\text{th}}}}{R_{\text{th}}}+C_{\text{th}}\frac{\text{d}{U}_{C_{\text{th}}}}{\text{d}t}$

其中，U_oc表示开路电压；R_o表示电池串联内阻；R_th表示电池极化电阻；C_th表示电池极化电容；I_L表示电池输出电流；U_L表示电池输出电压。

Figure 1. First-order Thevenin equivalent circuit model

图1. 一阶Thevenin等效电路模型

上面的表达式在通过拉普拉斯变换后可以得到如下的输入–输出的差分方程，令$\tau =R_{\text{th}}{C}_{\text{th}}$ ，则

$\tau U_{\text{oc}}\left(s\right)s+U_{\text{oc}}\left(s\right)=R_{\text{o}}\tau I_{\text{L}}\left(s\right)s+\cdots +\left(R_{\text{o}}+R_{\text{th}}\right)I_{\text{L}}\left(s\right)+\tau U_{\text{L}}\left(s\right)s+U_{\text{L}}\left(s\right)$

令$s=\frac{x\left(k\right)-x\left(k-1\right)}{\text{Δ}t}$ , $a=\tau$ , $b=R_0\tau$ , $c=R_0+R_{\text{th}}$ ，其中$\text{Δ}t$ 表示采样时间

$\begin{array}{c}{U}_{\alpha }\left(k\right)-U_{\text{L}}\left(k\right)=\frac{a}{\Delta t+a}\left[U_{\alpha }\left(k-1\right)-U_{\text{L}}\left(k-1\right)\right]+\frac{c\Delta t+b}{\Delta t+a}{I}_{\text{L}}\left(k\right)+\frac{-b}{\Delta t+a}{I}_{\text{L}}\left(k-1\right)\end{array}$

令$k_1=\frac{a}{\text{Δ}t+a}$ , $k_2=\frac{c\text{Δ}t+b}{\text{Δ}t+a}$ , $k_3=\frac{-b}{\text{Δ}t+a}$ ，则上式表示为

$U_{\alpha }\left(k\right)-U_L\left(k\right)=k_1\left[U_{\alpha }\left(k-1\right)-U_L\left(k-1\right)\right]+k_2{I}_L\left(k\right)+k_3{I}_L\left(k-1\right)$

从上述公式可知，锂离子电池在k时刻的输出电压，是$k-1$ 时刻输出电压与k、$k-1$ 时刻输出电流的线性组合形式，这表明其具有自回归结构。因此，能够构建出通过在线数据辨识来计算锂离子电池输出电压值的数学模型。

而自回归线性模型是在线数据建模的一种常用方法，其大体上可以表示为如下形式

$y\left(k\right)={\displaystyle \sum }_{i=1}^N{\alpha }_{i}y\left(k-i\right)+{\displaystyle \sum }_{j=0}^M{\beta }_{j}u\left(k-j\right)$

用最小二乘法为例来说明系数的整定过程：比如给定一批数据$\left\{y_i\in R,x_i\in R^n,i=1,\cdots ,N\right\}$ ，使其满足$y_i=x_i^{\text{T}}\theta$ ，此时辨识$n\times 1$ 维参数向量$\theta =\left[{\theta }_1,\cdots ,{\theta }_n\right]$ ，写成矩阵形式为$Y_N=X_N\theta$ ，其中的$Y_N\in R^N$ ，$X_N$ 是个$N\times n$ 的矩阵，矩阵的每一行对应一个$x_i^{\text{T}}$ 。估计的参数$\widehat{\theta }$ 越准，则模型的预测值$X_N\widehat{\theta }$ 的结果越接近$Y_N$ ，因此我们将误差最小化作为研究参数的目标

$\min {\Vert Y_N-X_N\widehat{\theta }\Vert }_2$

对$\widehat{\theta }$ 求导可得$X_N^{\text{T}}\left(Y_N-X_N\widehat{\theta }\right)=0$，则

$\widehat{\theta }={\left(X_N^{\text{T}}{X}_N\right)}^{-1}{X}_N^{\text{T}}{Y}_N$

由上述公式可见，若在t时刻获取N组数据，便可用该公式对模型参数$\widehat{\theta }$ 进行辨识。然而，若针对每个时刻的新数据都对全部历史数据进行求解，会因计算量巨大而耗费大量时间与存储资源。所以，可采用最小递推二乘法进行计算。此算法可理解为，每一步更新模型参数时，都是以上一时刻得到的模型参数为基础，并依据当前时刻的数据予以修正[1]，从而有效减少计算负担，提升计算效率。即

${\widehat{\theta }}_t={\widehat{\theta }}_{t-1}+\text{Δ}{\widehat{\theta }}_t$

根据这一思想，获得具体计算过程如下：定义变量$P_N^{-1}=\left(X_N^{\text{T}}{X}_N\right)$ ，即$P_N^{-1}=\left[x_1,x_2,\cdots ,x_N\right]\left[x_1^{\text{T}};x_2^{\text{T}};\cdots ;x_N^{\text{T}}\right]={\displaystyle {\sum }_{i=1}^{N-1}{x}_i{x}_i^{\text{T}}}+x_N{x}_N^{\text{T}}=P_{N-1}^{-1}+x_N{x}_N^{\text{T}}$ 同理，$X_N^{\text{T}}{Y}_N=X_{N-1}^{\text{T}}{Y}_{N-1}+x_N{y}_N$ 。根据定义的变量对上述公式重新简化后可得

$\begin{array}{c}{\widehat{\theta }}_t=P_N\left(P_{N-1}^{-1}{\widehat{\theta }}_{t-1}+x_N{y}_N\right)\\ =P_N\left(P_{N-1}^{-1}-x_N{x}_N^{\text{T}}\right){\widehat{\theta }}_{t-1}+P_N{x}_N{y}_N\\ ={\widehat{\theta }}_{t-1}+P_N{x}_N{y}_N-P_N{x}_N{x}_N^{\text{T}}{\widehat{\theta }}_{t-1}\\ ={\widehat{\theta }}_{t-1}+P_N\left(x_N{y}_N-x_N{x}_N^{\text{T}}{\widehat{\theta }}_{t-1}\right)\end{array}$

因为$P_N^{-1}=P_{N-1}^{-1}+x_N{x}_N^{\text{T}}$ ，并且根据矩阵逆引理

${\left[A+BCD\right]}^{-1}=A^{-1}-A^{-1}B{\left[C^{-1}+D{A}^{-1}B\right]}^{-1}D{A}^{-1}$

可以获得递归最小二乘的解为

$\{\begin{array}{l}{P}_N=P_{N-1}-P_{N-1}{x}_N{\left[I+x_N^{\text{T}}{P}_{N-1}{x}_N\right]}^{-1}{x}_N^{\text{T}}{P}_{N-1}\\ {\widehat{\theta }}_t={\widehat{\theta }}_{t-1}+P_N{x}_N\left(y_N-x_N^{\text{T}}{\widehat{\theta }}_{t-1}\right)\end{array}$

其中，$y_N-x_N^{\text{T}}{\widehat{\theta }}_{t-1}$相当于用$t-1$ 时刻的参数来进行t时刻参数预测时所带来的偏差，$P_N{x}_N$ 则相当于修正系数。

在锂离子电池一阶等效电路模型中，令$y\left(k\right)=U_{\text{oc}}\left(k\right)-U_{\text{L}}\left(k\right)$ ，${\Phi }^{\text{T}}\left(\begin{array}{c}k\end{array}\right)=\left[\begin{array}{ccc}y\left(k-1\right)& I_{\text{L}}\left(\begin{array}{c}k\end{array}\right)& I_{\text{L}}\left(\begin{array}{c}k-1\end{array}\right)\end{array}\right]$ ，$\theta ={\left[\begin{array}{ccc}{k}_1& k_2& k_3\end{array}\right]}^{\text{T}}$ ，可以采用根据最小二乘法改进得到的递推最小二乘法来完成对模型参数$\widehat{\theta }\left(k\right)$ 估计，得到的公式如下：

$\{\begin{array}{l}\widehat{\theta }\left(k\right)=\widehat{\theta }\left(k-1\right)+K\left(k\right)\left(y\left(k\right)-{\Phi }^{\text{T}}\left(k\right)\widehat{\theta }\left(k-1\right)\right)\\ K\left(k\right)=P\left(k\right)\Phi \left(k\right)=\frac{P\left(k-1\right)\Phi \left(k\right)}{\left(I+\Phi \left(k\right)P\left(k-1\right){\Phi }^{\text{T}}\left(k\right)\right)}\\ P\left(k\right)=\left(I-K\left(k\right){\Phi }^{\text{T}}\left(k\right)\right)P\left(k-1\right)\end{array}$

在对上述模型参数$\theta ={\left[\begin{array}{ccc}{k}_1& k_2& k_3\end{array}\right]}^{\text{T}}$ 进行整定后，可以进一步将当前时刻该等效电路的物性参数进行计算，其计算过程如下，令$k_0=\text{Δ}t+a$ ，则

$\{\begin{array}{l}\tau =k_1\Delta t/\left(1-k_1\right)\\ R_0=-k_3\left(\Delta t+\tau \right)/\tau \\ R_{th}=\left(k_2+k_3\right)\left(\Delta t+\tau \right)/\Delta t-R_0\\ C_{th}=\tau /R_{th}\end{array}$

这种方法具备两方面显著优势：其一，降低了对电池离线测试数据的依赖程度，让数字镜像建模方法更契合已建成的锂离子电池储能系统，显著增强了项目方案的适用性；其二，通过该方法建立起从电池外特性电压到电池内特性电阻、电容等参数的估计关系，为储能电池未来运行状态的模拟提供了便利条件[2]。

在运用递归最小二乘法完成对电池物性参数的估计后，基于电池数字镜像系统在未来模拟场景的应用中，为提高模拟精度，需要深入探究电池SOC对电池物性参数变化的影响过程。因此，在第一阶段工作完成电池物性参数估算后，需继续借助最小二乘法拟合出电池SOC与不同物性参数之间的映射函数，进而形成本项目所提出的两阶段动态建模流程，其具体实现步骤如下图2所示。

3. 电池SOC估计方法研究

美国先进电池联合会(USABC)对SOC给出了一种明确的定义方式：在特定的温度条件下，当电池处于放电过程时，将电池所剩余的电量与电池的额定电量的比值，界定为该电池的SOC。其定义公式可表示为：

$C_{\text{soc}}=\frac{Q_{\text{Re}}}{Q_a}\times 100\text{\%}$

根据电池荷电状态的定义就可以看出来，荷电状态的变化主要表现为电池电流随着时间的电荷累积作用。所以，由安时积分法所得的电池荷电状态计算公式如下所示：所以，基于安时积分过程的SOC计算方式，适用于在线检测的粗放型锂离子电池场景。传统安时积分法的公式为：

$SOC\left(t\right)=SO{C}_0-\frac{1}{C_{\text{nom}}}{\displaystyle {\int }_0^t\eta I\left(\tau \right)\text{d}\tau }$

Figure 2. The two-stage dynamic modeling process of the equivalent circuit model of lithium-ion battery based on online data

图2. 基于在线数据的锂离子电池等效电路模型两阶段动态建模过程

该方法下，当前时刻电池的电荷状态，取决于前一时刻电池的荷电状况以及此段时间内电池荷电状态的变动量。其中，荷电状态的变动量可通过电流在时间上的累积效果，即电流与时间乘积的总和得到确定。基于上述原理，安时积分法所推导出的电池荷电状态计算公式如下：

$C_{\text{soc}}\left(k+1\right)=C_{\text{soc}}\left(k+1\right)+\frac{\eta {\displaystyle {\int }_k^{k+1}{I}_a\left(k\right)\text{d}t}}{Q_a}$

然而，该方法存在明显缺陷，在其积分累积进程中，测量所产生的误差会持续影响SOC的估计结果，使得计算误差不断增大。在线性定常系统中，卡尔曼估计算法是一种最优化递推估计方法，常被用于对于测量值的校正和估计，在电池荷电状态的估计中卡尔曼估计的改进方法成为了电池荷电状态估计中的重要方法。其工作原理为，通过利用电池电压测量值与电池模型输出电压的计算值之间的偏差，来校正由安时积分法计算出的电池荷电状态估计值，即可以实现对于电池荷电状态的实时估计，其具体的计算过程如下。

基于系统状态空间描述$\{\begin{array}{l}X\left(k+1\right)=AX\left(k\right)+B{I}_0\left(k\right)\\ U\left(k+1\right)=HX\left(k\right)\end{array}$ 的卡尔曼滤波基本公式为：

$\widehat{X}{\left(k+1\right)}^{-}=A\widehat{X}\left(k\right)+B{I}_0\left(k\right)$

$P_k^{-}=A{P}_{k-1}{A}^{\text{T}}+Q$

$K_k=P_k^{-}{H}^{\text{T}}{\left(H{P}_k^{-}{H}^{\text{T}}+R\right)}^{-1}$

$\widehat{X}\left(k\right)=\widehat{X}{\left(k\right)}^{-}+K_k\left(U_0\left(k\right)-H\widehat{X}{\left(k\right)}^{-}\right)$

$P_k=P_k^{-}-K_{k}H{P}_k^{-}$

其中，噪声协方差、过程噪声及误差协方差、修正矩阵、验前和验后的估计误差协方差矩阵分别表示为$R$ 、$Q$ 、$K_k$ 、$P_k$ 和$P_k^{-}$ ，k时刻系统状态变量估计值可以表示为$\widehat{X}\left(k\right)={\left[\begin{array}{cc}{C}_{\text{soc}}\left(k\right)& U_{\text{p}}\left(k\right)\end{array}\right]}^{\text{T}}$。

鉴于传统的卡尔曼估计方法受限于线性、定常特性，难以契合非线性、时变特性显著的电池系统。于是，基于线性化理念的扩展卡尔曼估计方法应运而生，并在相关领域得以广泛应用[3]。

拓展卡尔曼滤波的状态方程为：

$x_k=f\left(x_{k-1},u_k\right)+w_k$

$Z_k=g\left(X_k,U_k\right)+V_k$

观测方程为：

$y_k=h\left(x_k\right)+v_k$

时间更新方程为：

${\widehat{X}}_k^{-}=f\left({\widehat{X}}_{k-1},U_{k-1}\right)$

$P_k^{-}=A_{k-1}{P}_{k-1}{A}_{k-1}^{\text{T}}+Q_{k-1}$

状态变量为：

$x={\left[SOC{V}_{\text{rc}1}{V}_{\text{rc}2}\right]}^{\text{T}}$

状态更新方程为：

$K_k=P_k^{-}{H}_k^{\text{T}}{\left(H_k{P}_k^{-}{H}_k^{\text{T}}+R_k\right)}^{-1}$

${\widehat{X}}_k={\widehat{X}}_k+K_k\left(Z_k-g\left({\widehat{X}}_k,U_k\right)\right)$

$P_k=\left(I-K_k{H}_k\right)P_k^{-}$

将该拓展卡尔曼算法应用于电池二阶等效模型后可得：

$\left[\begin{array}{c}{V}_1\left(k\right)\\ V_2\left(k\right)\\ SO{C}_k\end{array}\right]=A_k\ast \left[\begin{array}{c}{V}_1\left(k-1\right)\\ V_2\left(k-1\right)\\ SO{C}_{k-1}\end{array}\right]+B_k\ast i_{k-1}+W_{k-1}$

$V_B\left(k\right)=\left[\begin{array}{ccc}-1& -1& 0\end{array}\right]\left[\begin{array}{c}{V}_1\left(k\right)\\ V_2\left(k\right)\\ SO{C}_k\end{array}\right]-i_k\ast R\left(k\right)+U_{oc\nu }+V_k$

考虑电池数字镜像模型存在物性参数时变特点，可以基于扩展卡尔曼方式，利用电池的测量电压和数字镜像模型输出电压，修正安时积分对于电池SOC的计算值。为检验所运用的卡尔曼估计方法的实际效果，依据LFP电芯的测试数据展开相应计算，关于电压以及荷电状态(SOC)的误差结果如下图3所示。

从下述图中结果对比可以看出，基于电池测量电压修正电池SOC的Kalman估计方法可以有效地降低安时积分法造成的误差累计，降低了电池SOC估计误差，同时降低了电池在低SOC条件下电池镜像系统的模型电压与测量电压的误差，有效地保障了电池镜像系统的精度。

4. SOH的估计方法研究

在锂离子数字建模和电池SOC估计的基础上，对锂电池的电池健康状态估计进行研究，从电池的全生命周期放电量入手，提出一种新的基于累计放电量的SOH计算方法，用电池累计放电量一定程度上代替等效循环次数，从而减少等效循环时因复杂反应导致的误差。

锂离子电池老化过程中，剩余循环次数、电池内阻以及电池容量是三个不可或缺的关键因素，如下图4所示。

(a)

(b)

Figure 3. Comparison of the definition of state of health of lithium-ion batteries. (a) The results of the ampere-hour integral method; (b) Results based on the Kalman filter

图3. 锂离子电池健康状态定义方式对比图。(a) 安时积分法结果；(b) 基于卡尔曼滤波法的结果

从图4中可以看出，对电池SOH (健康状态)计算指标的选择对于准确地评估电池性能起到至关重要的作用。使用电池容量和电阻作为计算SOH指标的方法能够直接反映电池的储能水平和内阻情况。但在实际的运用场景里，保证电池进行完全的充放电的过程很难做到；同时，电阻的测量也往往受到电池工作状态、环境温度、电池内部复杂电化学过程等多种因素的影响，而使用循环使用寿命作为电池SOH的计算指标则更加贴近电池的实际使用情况。电池在使用的过程当中，随着充放电次数的不断增加，其性能也会逐渐地下降[4]。因此，通过计算电池实际应用过程的等效循环次数，可以更准确地评估电池的健康状态，其计算公式为：

$N_{\text{eq}}=\frac{{\displaystyle \sum DO{D}_i\cdot Q_{\text{discharge},i}}}{C_{\text{nom}}}$

其中的$DO{D}_i$ 为第i次的循环放电深度(0%~100%)，$Q_{\text{discharge,}i}$ 为第i次的放电容量(Ah)，$C_{\text{nom}}$ 为电池的标称容量(Ah)。尽管该方法的计算流程相对繁琐，不过它能够更为实际、精准地展现出电池性能的衰减状况。

Figure 4. Comparison of the definition of state of health of lithium-ion batteries

图4. 锂离子电池健康状态定义方式对比图

在近些年的研究中，部分研究人员开拓了新的应用方向，将工程领域中原本常用于材料寿命评估的雨流计数法应用到了锂离子电池等效循环次数的估算工作中。该方法的核心思路是把机械疲劳里的应力曲线和锂离子电池的SOC曲线进行等效替换。通过这样的替换，我们就能够对不同充放电深度下的循环次数进行折算，进而得出这些循环次数相对于满额充放电深度循环次数的等效数值[5]。

其核心步骤为：

1) 峰谷检测：识别SOC曲线中的局部极大值(峰)和极小值(谷)。

2) 循环配对：根据“雨流”规则匹配峰谷形成闭合循环。

3) 循环计数：统计等效循环次数并加权计算老化影响。

首先设电池的SOC时间序列为：

$\text{SOC}=\left\{s_1,s_2,\cdots ,s_n\right\}$

将其中的微小波动进行处理后，剩余较为明显的峰谷点。

第二步，定义峰点$P_i$ 和谷点$V_j$ ：

$P_i=\max \left(s_{i-1},s_i,s_{i+1}\right)$

且$s_i>s_{i-1},s_i>s_{i+1}$ ；

$V_j=\min \left(s_{j-1},s_j,s_{j+1}\right)$

且$s_j<s_{j-1},s_j<s_{j+1}$ ；

第三步，按照时间顺序挨个经过峰谷点，当某段路径形成闭合环路时，记为一次完整的循环。

最后对于每个闭合循环，计算其中的幅值和均值：

${\text{Range}}_k=\left|P_k-V_k\right|$

${\text{Mean}}_k=\frac{P_k+V_k}{2}$

还可以引入非线性加权系数来修正不同充放电深度(DOD)对电池老化带来的影响：

$N_{\text{eq},k}={\left(\frac{{\text{Range}}_k}{100}\right)}^{\alpha }\cdot {\left(\frac{T_k}{25}\right)}^{\beta }$

其中的$\alpha$ 为DOD指数；$\beta$ 为温度指数；$T_k$ 为循环期间的平均温度(℃)。

总的等效循环次数为：

$N_{\text{eq}}={\displaystyle \sum }_{k=1}^m{N}_{\text{eq},k}$

基于等效循环次数的SOH的衰减模型为：

$SOH=1-k\cdot N_{\text{eq}}$

其中k为衰减系数。

具体情况如下图5所示。

(a)

(b)

Figure 5. Schematic diagram of rainflow counting method. (a) SOC curve during battery charging and discharging; (b) Cyclic state of rainflow counting conversion

图5. 雨流计数法示意。(a) 电池充放电过程中SOC曲线；(b) 雨流计数法折算的循环状态

如图所示，在雨流计数法的计算下，可以将一段只有极值点的图像，分解为四个全循环，然而其等效精度仍有进一步提升的空间。因此，本文从锂离子电池全生命周期的放电电量角度出发，提出一种基于累计放电电量的SOH计算模式，同时注意到锂电池等效循环次数与最大循环次数可以用于表征锂电池的SOH [6]。因此，以锂电池的SOH作为纽带，能够建立起锂电池累计放电电量与等效循环次数之间的联系，进而构成如下的计算公式。

$SOH\left(k\right)=1-\frac{Q_s\left(k\right)}{Q_{\text{as}}\left(k\right)}=1-\frac{N_{\text{eq}}\left(k\right)}{N_{\text{cyc}}\left(k\right)}$

式中，Q_s和Q_as分别表示锂电池的累计放电电量以及全生命周期下的全部放电电量，而N_eq和N_cyc分别表示锂电池的等效循环次数以及最大循环次数，则$N_{\text{eq}}=Q_{\text{s}}÷Q_{\text{as}}\times N_{\text{cyc}}$ 。

结合上述用电池SOH来等效电池循环次数的计算方法，可以开展电池SOH的估计过程，其具体实现流程如下图6所示。

Figure 6. Estimation process of battery SOH based on rainflow counting method

图6. 基于雨流计数法的电池SOH估计流程