基于温度采样数据的锂电池执行器热故障诊断

Thermal Fault Diagnosis of Li-Ion Battery Actuator Based on Temperature Sampling Data

摘要: 针对锂电池执行器热故障问题，提出了一种基于集员滤波的故障诊断方法。通过建立锂电池电热耦合模型，基于电池温度采样数据，采用全对称多胞体卡尔曼滤波算法对采样数据进行分析，得到故障估计的区间，从而实现锂电池执行器热故障诊断。

Abstract: A fault diagnosis method based on set membership filtering is proposed for thermal fault of Li-ion battery. Through constructing electrothermal coupled model, based on the battery temperature sampling data, the Zonotopic Kalman filter algorithm is used to analyze the sampling data to obtain the interval of fault estimation, so as to realize the thermal fault diagnosis of Li-ion battery actuator.

1. 引言

锂电池是一种常用的动力电池，因其寿命长、自放率低、能量密度高等优势成为最具有发展前景的储能电源，被广泛应用于储能领域的各个方面 [1]。锂电池应用领域的不断拓展对电池的性能提出了较高的要求，安全性和可靠性问题是电池储能系统最关心的问题。锂电池热故障是电池故障中最重要的故障之一，它具有潜在的灾难性，会致使鼓包、漏液等问题，进而导致整个储能系统故障严重 [2]。因此，为保证电池储能系统的安全运行，对锂电池进行快速而准确的热故障诊断有着重要的意义。

现有的电池故障诊断方法可以分为两类，一种是基于电池模型的故障诊断方法，另一种是无电池模型的故障诊断方法 [3]。基于电池模型的故障诊断方法是通过比较电池模型的预测值与测量值之间的差异获取残差信号，进一步对信号处理以获取故障信息。该诊断方法可以实时进行电池故障检测和诊断，在电池的故障诊断技术领域中占据重要的地位。然而在实际情况中，系统中噪声的存在会导致实际系统建模不准确，基于电池模型的故障诊断方法的性能也因此受到影响。

本文主要针对上述基于电池模型的故障诊断方法中存在的不足，提出一种基于集员滤波思想的电池故障诊断方法，通过建立锂电池电热耦合模型，获取锂电池核心温度和表面温度的采样数据，将其作为全对称多胞体卡尔曼滤波算法的输入，从而得到锂电池执行器热故障的估计区间，根据故障估计区间可判断系统的故障状况。

2. 模型建立

锂电池热故障指导致电池温度异常的故障，包括电池滥用、风机故障和电池过充电故障等。热故障的普遍表现为电池核心温度和表面温度异常。因此需要建立电池热模型来模拟核心温度和表面温度的动态行为。考虑到电池的电学参数如电池电流、端电压等会影响到电池的产热，结合锂电池二阶Thevenin等效电路模型，建立电热耦合模型。

Figure 1. Schematic diagram of second-order Thevenin model

图1. 二阶Thevenin模型示意图

图1所示为二阶Thevenin模型， $R_1$ 和 $C_{p1}$ 分别为电化学极化电阻和电容， $R_2$ 和 $C_{p2}$ 分别为浓差化电阻和电容。 $U_1,U_2$ 为两个RC电路的压降， $U_o$ 为电池两端电压， $U_{oc}$ 为电池内部电压。R为电池内阻，I为电池电流。模型的数学表达式如下 [4]：

$\begin{array}{l}U=U_{oc}-R_{o}I-U_1-U_2\hfill \\ \{\begin{array}{l}{\stackrel{\dot{}}{U}}_1=-\frac{1}{R_1{C}_{p1}}{U}_1+\frac{1}{C_{p1}}I\\ {\stackrel{\dot{}}{U}}_2=-\frac{1}{R_2{C}_{p2}}{U}_2+\frac{1}{C_{p2}}I\end{array}\hfill \end{array}$ (1)

式中， ${\stackrel{\dot{}}{U}}_1$ 和 ${\stackrel{\dot{}}{U}}_2$ 分别表示电压 $U_1$ 和 $U_2$ 的变化率。

根据电池内外的能量守恒方程、傅里叶热定律和牛顿散热定律，将产热模型定义为 [5]：

$\{\begin{array}{l}{C}_c{\stackrel{\dot{}}{T}}_c=Q_{gen}+\frac{T_s-T_c}{R_c}\hfill \\ C_s{\stackrel{\dot{}}{T}}_{\text{s}}=\frac{T_e-T_s}{R_u}-\frac{T_s-T_c}{R_c}\hfill \end{array}$ (2)

式中， $C_c$ 和 $C_s$ 分别表示电池内部材料的热容系数和电池表面的热容系数， $R_c$ 表示电池核心与表面之间的等效热传导电阻， $R_u$ 表示电池表面与冷却空气之间的等效对流电阻。 $T_e$ 表示环境温度， $Q_{gen}$ 表示电池核心的发热功率。

根据Bernardi公式，电池的发热功率 $Q_{gen}$ 为：

$Q_{gen}=I\left(U_{oc}-U\right)=I\left(R_{o}I+U_1+U_2\right)$ (3)

以产热率 $Q_{gen}$ 和环境温度 $T_e$ 作为电热耦合模型的输入，以表面温度的测量值 ${\bar{T}}_s$ 和核心温度的测量值 ${\bar{T}}_c$ 作为电热耦合模型的输出。

离散化后，电热耦合模型的状态空间表达式为：

$\{\begin{array}{l}{x}_{k+1}=A{x}_k+B{u}_k+D_1{w}_k\hfill \\ y_k=C{x}_k+D_2{v}_k\hfill \end{array}$ (4)

式中， $x_k={\left[T_{c,k},T_{s,k}\right]}^{\text{T}}\in {ℝ}^n$ 为系统的状态向量， $u_k={\left[Q_{gen,k},T_{e,k}\right]}^{\text{T}}\in {ℝ}^r$ 为系统的输入向量， $y_k={\left[{\bar{T}}_{c,k},{\bar{T}}_{s,k}\right]}^{\text{T}}\in {ℝ}^p$ 为系统的输出向量， $w_k\in {ℝ}^w$ 是系统中的扰动噪声， $v_k\in {ℝ}^v$ 是系统中的测量噪声。 $\Delta t$ 是采样时间，并且

$A=\left[\begin{array}{cc}1-\frac{\text{Δ}t}{R_c{C}_c}& \frac{\text{Δ}t}{R_c{C}_c}\\ \frac{\text{Δ}t}{R_c{C}_s}& 1-\frac{\text{Δ}t}{R_c{C}_s}-\frac{\text{Δ}t}{R_u{C}_s}\end{array}\right],B=\left[\begin{array}{cc}\frac{\text{Δ}t}{C_c}& 0\\ 0& \frac{\text{Δ}t}{R_u{C}_s}\end{array}\right],C=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],D_1=D_2=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$

根据构建的电热耦合模型，设置相应的系统参数，可得到相对应的电池核心温度和表面温度的采样数据。

3. 集员滤波方法

集员滤波是集员估计的一种。集员估计是一类利用几何空间求解状态估计的方法，相比于传统的状态估计算法，集员估计方法不需要知道系统噪声和干扰的统计特性，只要求系统噪声有界 [6]。集员估计方法的关键之处在于找到一个包含状态或参数的可行集的简单几何体来近似可行集。目前常用的几何空间有区间，椭球和全对称多胞体。由于区间分析本身的缺陷，基于区间的集员算法收敛速度较慢。基于椭球的集员算法迭代简单，但包裹性差，因此保守性有限。基于全对称多胞体的集员算法可以转换为简单的矩阵运算，相比于其他几何空间，所需的计算量较少，保守性低 [7]。

全对称多胞体卡尔曼滤波(Zonotopic Kalman filter, ZKF)算法是将卡尔曼滤波与全对称多胞体集员估计结合。与卡尔曼滤波相似，ZKF可分为状态预测和测量更新。状态预测使用前一时刻的状态多胞体，干扰和噪声多胞体，以及系统的输入和测量输出来生成下一时刻状态的预测多胞体。测量更新步骤中，将预测多胞体与来自测量输出信息的带状空间相交，以获得当前时刻的状态多胞体。基于ZKF的故障诊断算法将系统故障向量扩展到系统状态向量，构造增广系统，通过ZKF算法获得状态的估计区间。当故障发生时，故障的估计区间也可以得到，由此可以判断系统是否发生故障以及故障类型。

4. 基于ZKF的锂电池执行器热故障诊断

基于电热耦合模型，根据状态空间表达式(4)，锂电池执行器热故障模型可写为：

$\{\begin{array}{l}{x}_{k+1}=A{x}_k+B{u}_k+f_k+D_1{w}_k\hfill \\ y_k=C{x}_k+D_2{v}_k\hfill \end{array}$ (5)

式中， $f_k={\left[f_{1,k},f_{2,k}\right]}^{\text{T}}\in {ℝ}^p$ 表示热故障向量，当锂电池发生产热故障时，例如过充\过放以及电池短路时， $f_1\ne 0,f_2=0$。当锂电池发生冷却系统故障时， $f_1=0,f_2\ne 0$。当锂电池内部导热电阻发生故障时， $f_1\ne 0,f_2\ne 0$。 $w_k$ 和 $v_k$ 分别表示系统未知但有界的扰动和测量噪声，即 $\left|w_k\right|\le \tilde{w}$， $\left|v_k\right|\le \tilde{v}$。

首先将故障向量扩展到电热耦合模型的状态向量，得到增广的状态向量如下式：

${\bar{x}}_k=\left[\begin{array}{c}{x}_k\\ f_{k-1}\end{array}\right]$ (6)

因此锂电池执行器热故障模型可以重新写成：

$\{\begin{array}{l}E{\bar{x}}_{k+1}=\bar{A}{\bar{x}}_k+\bar{B}{u}_k+{\bar{D}}_1{w}_k\hfill \\ y_k=\bar{C}{\bar{x}}_k+D_2{v}_k\hfill \end{array}$ (7)

其中

$\begin{array}{l}E=\left[\begin{array}{cccc}1& 0& -1& 0\\ 0& 1& 0& -1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],\bar{A}=\left[\begin{array}{cccc}1-\frac{\text{Δ}t}{R_c{C}_c}& \frac{\text{Δ}t}{R_c{C}_c}& 0& 0\\ \frac{\text{Δ}t}{R_c{C}_s}& 1-\frac{\text{Δ}t}{R_c{C}_s}-\frac{\text{Δ}t}{R_u{C}_s}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],\bar{B}=\left[\begin{array}{cc}\frac{\text{Δ}t}{C_c}& 0\\ 0& \frac{\text{Δ}t}{R_u{C}_s}\\ 0& 0\\ 0& 0\end{array}\right],\\ {\bar{D}}_1=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\\ 0& 0\end{array}\right],\bar{C}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]\end{array}$

不失一般性，假设系统(5)的初始状态 ${\bar{x}}_0$ 对应的全对称多胞体 ${\mathcal{Z}}_0={\bar{p}}_0\oplus {\bar{G}}_0{B}^s$， ${\bar{p}}_0$ 代表初始时刻对应全对称多胞体的中心点， ${\bar{G}}_0$ 代表初始时刻对应全对称多胞体的形状矩阵， ${{\rm B}}^s$ 为s个单位区间 $\left[-1,1\right]$ 构成的单位盒子。假设扰动噪声 $w_k$ 和测量噪声 $v_k$ 都可以被一个全对称多胞体包裹：

$\begin{array}{l}{w}_k\in \mathcal{W}=\langle 0,W\rangle \\ v_k\in \mathcal{V}=\langle 0,V\rangle \end{array}$ (8)

其中， $\mathcal{W},\mathcal{V}$ 分别为扰动和测量噪声所在的全对称多胞体集合。 $W=\text{diag}\left(\tilde{w}\right),V=\text{diag}\left(\tilde{v}\right)$。

假设k时刻的状态变量 ${\bar{x}}_k$ 对应的全对称多胞体 ${\mathcal{Z}}_k={\bar{p}}_k\oplus {\bar{G}}_k{B}^s=\langle {\bar{p}}_k,{\bar{G}}_k\rangle$，根据全对称多胞体的性质及系统状态方程可知，第 $k+1$ 时刻状态变量的预测集为

${\widehat{\stackrel{¯}{\chi }}}_{k+1}=\langle {\widehat{\stackrel{¯}{p}}}_{k+1},{\widehat{\stackrel{¯}{G}}}_{k+1}\rangle$ (9)

${\widehat{\stackrel{¯}{p}}}_{k+1}=T\bar{A}{\bar{p}}_k+T\bar{B}{u}_k+N{y}_{k+1}$ (10)

${\widehat{\stackrel{¯}{G}}}_{k+1}=\left[T\bar{A}{\downarrow }_{re}{\bar{G}}_{k}T{\bar{D}}_{1}W-N{D}_{2}V\right]$ (11)

其中，T和N可以看作加权矩阵，T和N的含义分别是从模型和测量输出中选择信息。T和N的可行解为：

$T={\text{Θ}}^{†}{\alpha }_1+S\text{Ψ}{\alpha }_1,N={\text{Θ}}^{†}{\alpha }_2+S\text{Ψ}{\alpha }_2$ (12)

其中，S为任意选取的矩阵， $\text{Θ}=\left[\begin{array}{c}E\\ \bar{C}\end{array}\right]$， ${\text{Θ}}^{†}$ 是 $\text{Θ}$ 的伪逆， $\text{Ψ}=I_{n+2p}-{\text{ΘΘ}}^{†}$， ${\alpha }_1=\left[\begin{array}{c}{I}_{n+p}\\ 0\end{array}\right]$， ${\alpha }_2=\left[\begin{array}{c}0\\ I_p\end{array}\right]$ ； ${\downarrow }_{re}{\bar{G}}_k$ 表示将k时刻状态所在全对称多胞体降阶后的生成矩阵，通过式(13)~(16)求得

$\mathcal{Z}=\langle \bar{p},\bar{G}\rangle \subseteq \bar{p}\oplus rs\left(\bar{G}\right)B^n$ (13)

$rs\left(\bar{G}\right)=\left[\begin{array}{cccc}{\displaystyle \underset{j=1}{\overset{s}{\sum }}\left|{\bar{G}}_{1,j}\right|}& 0& \cdots & 0\\ 0& {\displaystyle \underset{j=1}{\overset{s}{\sum }}\left|{\bar{G}}_{2,j}\right|}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\displaystyle \underset{j=1}{\overset{s}{\sum }}\left|{\bar{G}}_{n,j}\right|}\end{array}\right]$ (14)

$\mathcal{Z}=\langle \bar{p},\bar{G}\rangle \subseteq \langle \bar{p},{\downarrow }_{re}\bar{G}\rangle$ (15)

${\downarrow }_{re}\bar{G}=\{\begin{array}{ll}\bar{G},\hfill & s\le q\hfill \\ \left[\begin{array}{cc}{\bar{G}}_{>}& rs\left({\bar{G}}_{<}\right)\end{array}\right],\hfill & s>q\hfill \end{array}$ (16)

求得状态预测集 ${\widehat{\stackrel{¯}{\chi }}}_{k+1}$ 后，寻求 $k+1$ 时刻的带状空间 ${\mathcal{S}}_{k+1}$ 与状态预测集的交集，进而得到 $k+1$ 时刻的状态估计范围。

带状空间 ${\mathcal{S}}_{k+1}$ 是根据 $k+1$ 时刻系统输出变量构造得到，如下式所示：

${\mathcal{S}}_{k+1}=\left\{{\bar{x}}_{k+1}\in {ℝ}^{n+s}:\left|\bar{C}{\bar{x}}_{k+1}-y_{k+1}\right|\le D_2\tilde{v}\right\}$ (17)

根据 $k+1$ 时刻的带状空间 ${\mathcal{S}}_{k+1}$ 和状态变量预测集 ${\widehat{\stackrel{¯}{\chi }}}_{k+1}$ 的全对称多胞体集合，结合式(18)~(21)计算得到 $k+1$ 时刻的系统状态变量 ${\bar{x}}_{k+1}$ 的全对称多胞体集合：

${\bar{x}}_{k+1}\in {\bar{\chi }}_{k+1}=\langle {\bar{p}}_{k+1},{\bar{G}}_{k+1}\rangle$ (18)

${\bar{p}}_{k+1}={\widehat{\stackrel{¯}{p}}}_{k+1}+L_{k+1}\left(y_{k+1}-\bar{C}{\widehat{\stackrel{¯}{p}}}_{k+1}\right)$ (19)

${\bar{G}}_{k+1}=\left[\left(I_{n+p}-L_{k+1}\bar{C}\right){\widehat{\stackrel{¯}{G}}}_{k+1}{L}_{k+1}{D}_v\right]$ (20)

$L_{k+1}={\widehat{\stackrel{¯}{G}}}_{k+1}{\widehat{\stackrel{¯}{G}}}_{k+1}^{\text{T}}{\bar{C}}^{\text{T}}{\left(\bar{C}{\widehat{\stackrel{¯}{G}}}_{k+1}{\widehat{\stackrel{¯}{G}}}_{k+1}^{\text{T}}{\bar{C}}^{\text{T}}+D_v{D}_v^{\text{T}}\right)}^{-1}$ (21)

其中， $D_v=\text{diag}\left(D_2\tilde{v}\right)$。

对于一个全对称多胞体 $\mathcal{Z}=\langle p,G\rangle \subset {ℝ}^n$，能够被一个区间盒子 $\text{Box}\left(\mathcal{Z}\right)=\left[z^{-},z^{+}\right]$ 限定。 $z^{-}$ 和 $z^{+}$ 可以通过如下公式计算：

$\{\begin{array}{l}{z}^{-}\left(i\right)=p\left(i\right)-{\displaystyle \underset{j=1}{\overset{s}{\sum }}\left|G\left(i,j\right)\right|},i=1,\cdots ,n\\ z^{+}\left(i\right)=p\left(i\right)-{\displaystyle \underset{j=1}{\overset{s}{\sum }}\left|G\left(i,j\right)\right|},i=1,\cdots ,n\end{array}$ (22)

因此，结合全对称多胞体卡尔曼滤波算法，可以得到 ${\bar{x}}_{k+1}$ 的区间估计，计算公式如下：

$\{\begin{array}{l}{\bar{x}}_{k+1}^{+}\left(i\right)={\bar{p}}_{k+1}\left(i\right)+{\displaystyle \underset{j=1}{\overset{q}{\sum }}\left|{\bar{G}}_{k+1}\left(i,j\right)\right|},i=1,\cdots ,n+p\hfill \\ {\bar{x}}_{k+1}^{-}\left(i\right)={\bar{p}}_{k+1}\left(i\right)-{\displaystyle \underset{j=1}{\overset{q}{\sum }}\left|{\bar{G}}_{k+1}\left(i,j\right)\right|},i=1,\cdots ,n+p\hfill \end{array}$ (23)

式中， ${\bar{x}}_{k+1}^{+}\left(i\right),{\bar{x}}_{k+1}^{-}\left(i\right)$ 分别表示 $k+1$ 时刻状态估计的最小上界和最大下界。

根据式(24)可以确定时刻的故障估计区间

$\begin{array}{l}{f}_k^{+}=\left[\begin{array}{cc}0& I_p\end{array}\right]{\bar{x}}_{k+1}^{+}\\ f_k^{-}=\left[\begin{array}{cc}0& I_p\end{array}\right]{\bar{x}}_{k+1}^{-}\end{array}$ (24)

式中， $f_k^{+},f_k^{-}$ 分别表示k时刻故障估计的最小上界和最大下界。

当锂电池执行器不存在热故障时，则故障估计的上下界位于0的两侧。而当锂电池执行器发生热故障时，则故障估计的上下界同时位于0的一侧，故障的真实值则处于故障估计的区间范围内。故障类型不同，故障估计的区间不同，因此可以根据故障估计区间判断执行器热故障类型。

依据上述算法，基于温度采样数据，能够准确判断锂电池在工作过程中是否产生热故障，以及判断热故障类型，有效地保证了锂电池储能系统的安全性、可靠性和持久高效性。

5. 仿真示例

基于锂电池电热耦合模型，设置仿真模型参数如表1所示：

Table 1. Main parameters of the electrothermal coupling model of lithium battery

表1. 锂电池电热耦合模型主要参数

环境温度 $T_e$ 设为25℃，核心温度 $T_c$ 和表面温度 $T_s$ 的初始值都设为25℃。扰动和测量噪声限制为 $\left|w_k\right|\le 0.01,\left|v_k\right|\le 0.01$。锂电池执行器热故障 $f_k={\left[\begin{array}{cc}{f}_{1,k}& f_{2,k}\end{array}\right]}^{\text{T}}$ 具体形式如式(25)所示：

$f_k={\left[f_{1,k}{f}_{2,k}\right]}^{\text{T}}=\{\begin{array}{l}{\left[\begin{array}{cc}0& 0\end{array}\right]}^{\text{T}},0\le k<100\hfill \\ {\left[\begin{array}{cc}0.1& 0\end{array}\right]}^{\text{T}},100\le k<300\hfill \\ {\left[\begin{array}{cc}0.2& 0.2\end{array}\right]}^{\text{T}},300\le k<500\hfill \\ {\left[\begin{array}{cc}0.2& 0.4\end{array}\right]}^{\text{T}},k\ge 500\hfill \end{array}$ (25)

按照设置的参数运行，采样到的锂电池核心温度值和表面温度值如图2所示。

Figure 2. Core and surface temperature samples

图2. 核心温度和表面温度采样值

基于温度采样数据，采用全对称多胞体卡尔曼滤波算法对采样数据进行分析。设置 $p_0={\left[\begin{array}{cc}25& 25\end{array}\right]}^{\text{T}}$， $G_0=0.1I_4$，设置s矩阵：

$S=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\end{array}\right]$ (26)

根据式(12)，可以得到矩阵T和N如下式所示：

$T=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ -1& 0& 1& 0\\ 0& -1& 0& 1\end{array}\right],N=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 1& 0\\ 0& 1\end{array}\right]$ (27)

锂电池执行器热故障和它们的区间估计结果如图3所示。实际的故障用红线表示。从图3中可以看出，全对称多胞体卡尔曼滤波算法可以快速跟踪实际的执行器故障，并给出准确的区间估计结果。对于锂电池电热耦合模型而言，当系统出现故障时，ZKF可以检测到故障，并且可以迅速、准确地判断故障的类型，同时，根据故障估计的最大下界和最小上界，在一定的误差范围内，可以得到锂电池热故障的估计值。由此可见，ZKF算法准确地实现了锂电池执行器热故障诊断。

(a) (b)

Figure 3. Lithium battery thermal fault and its estimation intervals

图3. 锂电池热故障和故障估计区间

6. 结论

本文针对锂电池执行器热故障诊断问题，提出了一种基于集员滤波的电池故障诊断方法，通过获得锂电池的温度采样数据，采用全对称多胞体卡尔曼滤波算法得到系统故障的估计区间，若锂电池运行过程中某一时刻有某类热故障发生，则该故障估计的上下界同时位于0的一侧，从而实现锂电池执行器热故障诊断。

基于集员滤波的故障诊断方法有着准确性好、鲁棒性强、应用范围广等优点，因此本文将基于集员滤波的故障诊断方法运用到锂电池执行器的热故障诊断中，能够有效提高锂电池热故障诊断的质量和性能。

NOTES

^*通讯作者。