电池等效电路模型极化参数估计的收敛性分析

Convergence Analysis on Polarization Parameter Estimator of Battery Equivalent Circuit Models

摘要

为探究动力电池极化参数估计在不同工况下的误差收敛特性，采用李雅普诺夫第二法对估计器的稳定性和收敛性进行分析.基于等效电路模型和遗忘因子递归最小二乘法推导得到极化参数估计误差的状态方程.采用李雅普诺夫第二法对状态方程进行分析，得到估计器渐进收敛的必要条件，即持续变化的电流输入，并提出一种图解法用于分析误差的动态收敛特性和论证此必要条件的合理性.采用实验校准的电池模型产生的数据对理论分析过程和结果进行验证.结果表明，估计器在持续变化的电流输入下能逐渐收敛到真值附近，且在变化剧烈的正负交变工况下具备更好的收敛性.

Abstract

To investigate the error convergence characteristics of power battery polarization parameter estimation under different operating conditions， the stability and convergence of the estimators are rigorously analyzed in this paper by using Lyapunov’s second method. The state equations of polarized parameter estimation errors are firstly derived based on the equivalent circuit model and the forgetting factor recursive least squares （FFRLS）. Subsequently， the state equation is analyzed using Lyapunov’s second method to obtain the necessary condition for the asymptotic convergence of the estimator， i.e.， continuously varying current input. A graphical method is proposed for analyzing the dynamic convergence properties of the error and justifying this necessary condition. Finally， the theoretical analysis process and results are validated using data generated from an experimentally calibrated battery model. The results show that the estimator can gradually converge to the vicinity of the true value under continuously varying current inputs， and the drastically varying positive and negative alternating conditions have better convergence properties.

Graphical abstract

关键词

电池管理系统 / 等效电路模型 / 参数辨识 / 渐进收敛性 / 李雅普诺夫第二法

在电动汽车、智能电网等许多应用中，动力电池内部状态和参数估计是能源管理系统的关键. 准确的电池参数估计依赖于精确且计算高效的电池模型及算法^［1］. 近年来随着研究的深入，参数估计算法层出不穷，电池模型也朝着高精度、高效率的方向发展. 然而，模型参数的在线辨识仍是该领域的一大难题，研究者提出了多种在线估计等效电路模型（equivalent circuit model，ECM）参数的算法，包括最小二乘法及其变体（如递归最小二乘法和加权最小二乘滤波器）^［2］、双扩展卡尔曼滤波器 ^［3］、自适应无迹卡尔曼滤波器 ^［4］、粒子滤波器^［5］等. 但这些研究成果并未在工业界得到广泛应用，最重要的原因是参数在线辨识的鲁棒性欠佳，辨识结果在某些情况下携带较大的误差甚至发散. 鉴于此，有学者对在线辨识时参数的可辨识性及其与输入数据特性的关系展开了研究.

Sharma等^［6］采用Fisher信息矩阵分析了一阶ECM中内阻和容量参数的可辨识性，并通过克拉默-拉奥（Cramer-Rao， CR）边界量化了联合估计误差的下界与传感器误差、电流大小和电池开路电压（open circuit voltage， OCV）‐荷电状态 （state of charge， SoC） 曲线斜率的关系. Rothenberger等^［7］使用同样的方法分析了参数CR边界与电流轨迹的幅值和频率的关系，并据此设计了一条可提高辨识精度的电流轨迹. 在上述研究的基础上，Lin等^［8］分析了SoC、容量和内阻在典型电流输入下单独估计和联合估计的CR边界变化，并得到了联合估计相比单独估计时放大系数的表达式. Lin^［9］推导了极化参数CR边界关于电流的频率响应表达式. Song等^［10］在此基础上得到了OCV、欧姆内阻和极化参数的CR边界关于电流的频率响应表达式，并测试了不同电流输入以及不同参数联合估计时的估计精度. 研究者意识到输入数据特征对参数辨识精度的重要性，Song等^［11］提出了使用多频率滤波器分离电流以提高辨识精度的方法，并据此设计了一种顺序估计算法^［12］. 各种噪声和参数初值误差也会影响估计精度. Lin^［13］分析了存在测量偏差时SoC的稳态收敛情况. Shen等^［14］分析了容量和SoC的初值存在误差时SoC误差的稳态收敛值与修正增益之间的关系. Zhou等^［15］基于参数的置信区间评估了二阶ECM中各参数的可辨识性，并探究了数据量、采样频率和传感器噪声等对可辨识性的影响. Andersson等^［16］指出，并非所有模型参数都能被准确辨识，也并非所有参数都与模型性能相关，参数辨识不仅需要了解其所针对的物理系统，还应理解参数辨识背后的理论和概念.

上述研究虽量化了参数可辨识性与输入数据之间的关系，但并未解决参数估计器在不同工况下的收敛性这一重要问题. 为填补该项研究空白，本文基于李雅普诺夫第二法推导得到了估计器渐进收敛的必要条件，并创新性地提出了一种图解法，详细分析估计器的收敛机制. 本文的研究成果对电池模型参数在线估计算法的设计具有一定指导意义.

1 参数估计误差的状态方程

1.1 基于一阶ECM的参数估计

本节基于一阶ECM设计电池极化参数的在线估计器^［17-18］. 一阶ECM结构如图1所示.

根据基尔霍夫定律和电池SoC定义，可得一阶ECM的状态方程和输出方程如下：

V_{\mathrm{p},k}=V_{\mathrm{p},k-\mathrm{1}}{\mathrm{e}}^{-\mathrm{\Delta }t/R_{\mathrm{1}}{C}_{\mathrm{1}}}+I_{k-\mathrm{1}}{R}_{\mathrm{1}}\left(\mathrm{1}-{\mathrm{e}}^{-\mathrm{\Delta }t/R_{\mathrm{1}}{C}_{\mathrm{1}}}\right)（1）

{\xi }_k={\xi }_{k-\mathrm{1}}+{\eta }_{\mathrm{c}}{I}_{k-\mathrm{1}}\mathrm{\Delta }t{\left(\mathrm{3}\mathrm{600}C\right)}^{-\mathrm{1}}（2）

U_k=V_{\mathrm{o}\mathrm{c}}\left({\xi }_k\right)+V_{\mathrm{p},k}+I_k{R}_{\mathrm{0}}（3）

式中：\xi为电池的SoC； V_{\mathrm{o}\mathrm{c}}\left({\xi }_k\right)为电池的开路电压，是关于SoC的单调递增函数；V_{\mathrm{p},k}和U_k分别为极化电压和端电压；R_{\mathrm{1}}和C_{\mathrm{1}}表示极化电阻和极化电容；R_{\mathrm{0}}为欧姆内阻；C为电池容量；{\eta }_{\mathrm{c}}为库伦效率；\mathrm{\Delta }t为采样周期；I为电流，规定充电方向为电流正方向；k为采样时刻.

由式 （1） 可知，R_{\mathrm{1}}\mathrm{、}{C}_{\mathrm{1}}与极化电压之间的关系为非线性，直接辨识这两个参数十分困难^［19-20］. 因此，本文采用文献［21］中的方法，通过辨识两个与R_{\mathrm{1}}\mathrm{、}{C}_{\mathrm{1}}相关的参数来完成对R_{\mathrm{1}}\mathrm{、}{C}_{\mathrm{1}}的辨识. 令

\left\{\begin{array}{l}{p}_{\mathrm{1}}={\mathrm{e}}^{-\mathrm{\Delta }\mathrm{t}/{\mathrm{R}}_{\mathrm{1}}{\mathrm{C}}_{\mathrm{1}}}\\ p_{\mathrm{2}}=R_{\mathrm{1}}\left(\mathrm{1}-{\mathrm{e}}^{-\mathrm{\Delta }\mathrm{t}/{\mathrm{R}}_{\mathrm{1}}{\mathrm{C}}_{\mathrm{1}}}\right)\end{array}\right.（4）

令 y_k = V_{\mathrm{p},k}，将式（1）转化为如下辨识方程：

y_k=y_{k-\mathrm{1}}{p}_{\mathrm{1},k}+I_{k-\mathrm{1}}{p}_{\mathrm{2},k}（5）

很明显，式（5）是关于待辨识参数的线性方程. 为了便于阅读，将式（5）表示成矩阵形式：

y_k={\mathit{\phi }}_k{\mathit{\theta }}_k（6）

式中：{\mathit{\phi }}_k和{\mathit{\theta }}_k分别为辨识方程的输入向量和待辨识参数向量.

\left\{\begin{array}{c}{\mathit{\phi }}_k=\left[\begin{array}{cc}{y}_{k-\mathrm{1}}& I_{k-\mathrm{1}}\end{array}\right]\\ {\mathit{\theta }}_k={\left[\begin{array}{cc}{p}_{\mathrm{1},k}& p_{\mathrm{2},k}\end{array}\right]}^{\mathrm{T}}\end{array}\right.（7）

参数修正过程可写成：

{\widehat{\mathit{\theta }}}_k={\widehat{\mathit{\theta }}}_{k-\mathrm{1}}+{\mathit{K}}_k\left(y_k^m-{\widehat{y}}_k\right)（8）

式中：{\mathit{K}}_k={\left[K_{\mathrm{1},k},K_{\mathrm{2},k}\right]}^{\mathrm{T}}为估计器的更新增益；y_k^m为y_k的伪测量值，由电压、电流的测量值以及SoC、R₀的 估计值计算得到；右上标m表示变量的测量值；上 标 ^ 表示变量的估计值.

y_k^m=U_k-V_{\mathrm{o}\mathrm{c}}\left({\xi }_k\right)-I_k{R}_{\mathrm{0}}（9）

不失一般性，本文选择最基本的遗忘因子递归最小二乘法（forgetting factor recursive least squares，FFRLS）作为辨识算法. 算法的增益向量{\mathit{K}}_k和协方差矩阵{\mathit{P}}_{\mathrm{f},k}的更新方程如下：

\left\{\begin{array}{l}{\mathit{K}}_k={\mathit{P}}_{\mathrm{f},k-\mathrm{1}}{\mathit{\phi }}_k^{\mathrm{T}}{\left[{\lambda }_{\mathrm{f}}+{\mathit{\phi }}_k{\mathit{P}}_{\mathrm{f},k-\mathrm{1}}{\mathit{\phi }}_k^{\mathrm{T}}\right]}^{-\mathrm{1}}\\ {\mathit{P}}_{\mathrm{f},k}={\lambda }_{\mathrm{f}}^{-\mathrm{1}}\left[\mathit{I}-{\mathit{K}}_k{\mathit{\phi }}_k\right]{\mathit{P}}_{\mathrm{f},k-\mathrm{1}}\end{array}\right.（10）

式中：{\lambda }_{\mathrm{f}}为辨识过程的遗忘因子，是一个常数.

1.2 误差状态方程

为分析估计误差的动态收敛特性，需先得到其状态方程，再根据李雅普诺夫第二法^［22-23］对系统的稳定性进行分析. 首先，记各变量的误差如下：

\left\{\begin{array}{c}{\widehat{\xi }}_k={\xi }_k^{\mathrm{*}}-\mathrm{\Delta }{\widehat{\xi }}_k\\ {\widehat{V}}_{\mathrm{p},k}=V_{\mathrm{p},k}^{\mathrm{*}}-\mathrm{\Delta }{\widehat{V}}_{\mathrm{p},k}\\ {\widehat{\theta }}_{\mathrm{1},k}={\theta }_{\mathrm{1},k}^{\mathrm{*}}-\mathrm{\Delta }{\widehat{\theta }}_{\mathrm{1},k}\\ {\widehat{\theta }}_{\mathrm{2},k}={\theta }_{\mathrm{2},k}^{\mathrm{*}}-\mathrm{\Delta }{\widehat{\theta }}_{\mathrm{2},k}\\ {\widehat{R}}_{\mathrm{0},k}=R_{\mathrm{0},k}^{\mathrm{*}}-\mathrm{\Delta }{\widehat{R}}_{\mathrm{0},k}\\ I_k^m=I_k^{\mathrm{*}}-e_{I,k}\\ U_k^m=U_k^{\mathrm{*}}-e_{U,k}\end{array}\right.（11）

式中：右上标*表示变量的真值；\mathrm{\Delta }表示变量的误差；e_{I,k}和e_{U,k}分别为电流和电压的测量噪声.

联立式（1）、式（5）~式（7）、式 （11），经简化整理后可得：

y_k^m-{\widehat{y}}_k=\left(V_{\mathrm{p},k-\mathrm{1}}^{\mathrm{*}}+e_{\mathrm{v}\mathrm{p},k-\mathrm{1}}\right)\mathrm{\Delta }{\widehat{\theta }}_{\mathrm{1},k-\mathrm{1}}+

\left(I_{k-\mathrm{1}}^{\mathrm{*}}-e_{I,k-\mathrm{1}}\right)\mathrm{\Delta }{\widehat{\theta }}_{\mathrm{2},k-\mathrm{1}}+e_{\mathrm{v}\mathrm{p},k}+e_{I,k-\mathrm{1}}{\theta }_{\mathrm{2},k-\mathrm{1}}^{\mathrm{*}}-

e_{\mathrm{v}\mathrm{p},k-\mathrm{1}}{\theta }_{\mathrm{1},k-\mathrm{1}}^{\mathrm{*}}（12）

式中：

e_{\mathrm{v}\mathrm{p},k}=\chi \left(\mathrm{\Delta }{\widehat{\xi }}_{k-\mathrm{1}}+\frac{{\eta }_{\mathrm{c}}{I}_{k-\mathrm{1}}^{\mathrm{*}}\mathrm{\Delta }t}{\mathrm{3}\mathrm{600}{C}^{\mathrm{*}}}\right)+\mathrm{\Delta }{\widehat{R}}_{\mathrm{0},k}{I}_k^{\mathrm{*}}-e_{U,k}+

e_{I,k}\left(R_{\mathrm{0},k}^{\mathrm{*}}-\mathrm{\Delta }{\widehat{R}}_{\mathrm{0},k}\right)（13）

式中：\chi为\mathrm{O}\mathrm{C}\mathrm{V}-\mathrm{S}\mathrm{o}\mathrm{C}曲线在相应SoC区间内的斜率.

将式（12）代入式（8），可得参数估计误差的状态方程，并整理成如下矩阵形式：

\left[\begin{array}{c}\mathrm{\Delta }{\widehat{\theta }}_{\mathrm{1},k}\\ \mathrm{\Delta }{\widehat{\theta }}_{\mathrm{2},k}\end{array}\right]={\mathit{G}}_k\left[\begin{array}{c}\mathrm{\Delta }{\widehat{\theta }}_{\mathrm{1},k-\mathrm{1}}\\ \mathrm{\Delta }{\widehat{\theta }}_{\mathrm{2},k-\mathrm{1}}\end{array}\right]+\left[\begin{array}{c}{\theta }_{\mathrm{1},k}^{\mathrm{*}}-{\theta }_{\mathrm{1},k-\mathrm{1}}^{\mathrm{*}}\\ {\theta }_{\mathrm{2},k}^{\mathrm{*}}-{\theta }_{\mathrm{2},k-\mathrm{1}}^{\mathrm{*}}\end{array}\right]+

\left[\begin{array}{c}-K_{\mathrm{1},k}\\ -K_{\mathrm{2},k}\end{array}\right]\left(e_{\mathrm{v}\mathrm{p},k}+e_{I,k-\mathrm{1}}{\theta }_{\mathrm{2},k-\mathrm{1}}^{\mathrm{*}}-e_{\mathrm{v}\mathrm{p},k-\mathrm{1}}{\theta }_{\mathrm{1},k-\mathrm{1}}^{\mathrm{*}}\right)（14）

式中：

{\mathit{G}}_k=\left[\begin{array}{cc}\mathrm{1}-K_{\mathrm{1},k}{y}_{k-\mathrm{1}}^m& -K_{\mathrm{1},k}{I}_{k-\mathrm{1}}^m\\ -K_{\mathrm{2},k}{y}_{k-\mathrm{1}}^m& \mathrm{1}-K_{\mathrm{2},k}{I}_{k-\mathrm{1}}^m\end{array}\right]（15）

估计器的稳定性由系统矩阵决定，估计器稳定，误差可收敛. 然而，在误差收敛过程中，会受到各种外部误差的干扰，影响估计误差最终的收敛程度，也有可能影响估计器的稳定性. 从式（14）中可以看到，外部误差既会影响系统矩阵，又会引入额外估计误差，最终引起稳态误差. 虽然外部干扰会影响系统矩阵，但系统矩阵结构并未改变，系统稳定性是否受影响还需进一步分析.

2 误差动态特性的理论分析

2.1 估计器的渐近稳定性条件分析

得到估计器的误差状态方程的系统矩阵后，便可根据李雅普诺夫第二法分析系统的稳定性. 估计器本质上是一种线性时变的离散系统，该系统的能量函数可设为V\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)：

\begin{array}{l}V\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)=P_{\mathrm{11},k}\mathrm{\Delta }{\theta }_{\mathrm{1},\mathrm{k}}^{\mathrm{2}}+\mathrm{2}{P}_{\mathrm{12},k}\mathrm{\Delta }{\theta }_{\mathrm{1},k}\mathrm{\Delta }{\theta }_{\mathrm{2},k}+\\ P_{\mathrm{22},k}\mathrm{\Delta }{\theta }_{\mathrm{2},\mathrm{k}}^{\mathrm{2}}=\mathrm{\Delta }{\mathit{\theta }}_k^{\mathrm{T}}{\mathit{P}}_k\mathrm{\Delta }{\mathit{\theta }}_k\end{array}（16）

式中：{\mathit{P}}_k为待构造的正定对称矩阵.

结合式（14），可得误差在单步修正过程中的能量函数对时间的导数.

\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)=V\left(\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}\right)-V\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)=

\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}^{\mathrm{T}}{\mathit{P}}_{k+\mathrm{1}}\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}-\mathrm{\Delta }{\mathit{\theta }}_k^{\mathrm{T}}{\mathit{P}}_k\mathrm{\Delta }{\mathit{\theta }}_k=-\mathrm{\Delta }{\mathit{\theta }}_k^{\mathrm{T}}{\mathit{Q}}_k\mathrm{\Delta }{\mathit{\theta }}_{\mathit{k}}（17）

式中：

{\mathit{Q}}_k={\mathit{P}}_k-{\mathit{G}}_{k+\mathrm{1}}^{\mathrm{T}}{\mathit{P}}_{k+\mathrm{1}}{\mathit{G}}_{k+\mathrm{1}}（18）

通过\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)的表达式可知：若{\mathit{Q}}_k矩阵正定，\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)小于0，误差能量不断减小，则系统渐近稳定. 反之，系统不稳定.

将式（15）代入式（18），整理可得\mathit{Q}矩阵（为了便于阅读，此处省略了下标k）：

\mathit{Q}=

\left[\begin{array}{cc}\begin{array}{c}{K}_{P\mathrm{1}}{\phi }_{\mathrm{1}}+K_{P\mathrm{1}}{\phi }_{\mathrm{1}}-{\phi }_{\mathrm{1}}{\phi }_{\mathrm{1}}{K}_{PK}\\ K_{P\mathrm{2}}{\phi }_{\mathrm{1}}+K_{P\mathrm{1}}{\phi }_{\mathrm{2}}-{\phi }_{\mathrm{2}}{\phi }_{\mathrm{1}}{K}_{PK}\end{array}& \begin{array}{c}{K}_{P\mathrm{1}}{\phi }_{\mathrm{2}}+K_{P\mathrm{2}}{\phi }_{\mathrm{1}}-{\phi }_{\mathrm{1}}{\phi }_{\mathrm{2}}{K}_{PK}\\ K_{P\mathrm{2}}{\phi }_{\mathrm{2}}+K_{P\mathrm{2}}{\phi }_{\mathrm{2}}-{\phi }_{\mathrm{2}}{\phi }_{\mathrm{2}}{K}_{PK}\end{array}\end{array}\right]（19）

式中：

\left\{\begin{array}{c}{K}_{Pj}=K_{\mathrm{1}}{P}_{j\mathrm{1}}+K_{\mathrm{2}}{P}_{j\mathrm{2}},j=\mathrm{1,2}\\ K_{PK}=K_{\mathrm{1}}{K}_{P\mathrm{1}}+K_{\mathrm{2}}{K}_{P\mathrm{2}}\end{array}\right.

由此可得：

\left|\mathit{Q}\right|=\left(\mathrm{2}{K}_{P\mathrm{1}}{\phi }_{\mathrm{1}}-{\phi }_{\mathrm{1}}{\phi }_{\mathrm{1}}{K}_{PK}\right)\left(\mathrm{2}{K}_{P\mathrm{2}}{\phi }_{\mathrm{2}}-{\phi }_{\mathrm{2}}{\phi }_{\mathrm{2}}{K}_{PK}\right)-

{\left(K_{P\mathrm{1}}{\phi }_{\mathrm{2}}+K_{P\mathrm{2}}{\phi }_{\mathrm{1}}-{\phi }_{\mathrm{1}}{\phi }_{\mathrm{2}}{K}_{PK}\right)}^{\mathrm{2}}=

\mathrm{4}{\phi }_{\mathrm{1}}{\phi }_{\mathrm{2}}{K}_{P\mathrm{1}}{K}_{P\mathrm{2}}-{\left(K_{P\mathrm{1}}{\phi }_{\mathrm{2}}+K_{P\mathrm{2}}{\phi }_{\mathrm{1}}\right)}^{\mathrm{2}}=

-{\left(K_{P\mathrm{1}}{\phi }_{\mathrm{2}}-K_{P\mathrm{2}}{\phi }_{\mathrm{1}}\right)}^{\mathrm{2}}（20）

显然，无论\mathit{P}矩阵取何值，\left|\mathit{Q}\right|\le \mathrm{0}. 根据希尔维斯特判据^［24］，\mathit{Q}不是正定矩阵，系统在单步修正过程中的\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)不恒小于0，系统不稳定. 要使得系统渐近稳定，应找到一种\mathit{P}矩阵的构造方法，使得\mathit{Q}矩阵正定，\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)+\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}\right)在连续两步输入的修正过程中恒小于0.

要满足上述要求需要保证：系统在单步修正过程中的\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)小于或等于0；在连续两步不同输入下\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)和\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}\right)不会同时等于0.

2.2 估计器的渐近稳定性证明

首先，要使\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)在单步修正过程中小于或等 于0，对应的\mathit{Q}矩阵应保持半正定. 根据希尔维斯特判据，矩阵半正定应满足\left|\mathit{Q}\right|=0和顺序主子行列式Q_{\mathrm{\Delta }\mathrm{1}}>\mathrm{0}.

需要找到一种\mathit{P}矩阵的构造方法能使\left|\mathit{Q}\right|=\mathrm{0}. 经过分析，\mathit{P}矩阵可由协方差矩阵{\mathit{P}}_f通过式（21）变形后构造得到：

\mathit{P}=\left|{\mathit{P}}_f\right|{{\mathit{P}}_f}^{-\mathrm{1}}（21）

结合式（10），将其代入式（20），化简后可得：

\left|\mathit{Q}\right|=-\left[K_{\mathrm{1}}\left(P_{f\mathrm{12}}{\phi }_{\mathrm{1}}+P_{f\mathrm{22}}{\phi }_{\mathrm{2}}\right)-\right.

{\left.K_{\mathrm{2}}\left(P_{f\mathrm{12}}{\phi }_{\mathrm{2}}+P_{f\mathrm{11}}{\phi }_{\mathrm{1}}\right)\right]}^{\mathrm{2}}=

-{\left(K_{\mathrm{1}}{K}_{\mathrm{2}}-K_{\mathrm{2}}{K}_{\mathrm{1}}\right)}^{\mathrm{2}}{\left({\mathit{\phi }}_f{\mathit{P}}_f{\mathit{\phi }}_f^{\mathrm{T}}+{\lambda }_f\right)}^{\mathrm{2}}=\mathrm{0}（22）

同样地，结合式（10）、式（19），将构造的\mathit{P}矩阵代入Q_{\mathrm{\Delta }\mathrm{1}}，经化简后可得：

Q_{\mathrm{\Delta }\mathrm{1}}=K_{P\mathrm{1}}{\phi }_{\mathrm{1}}+K_{P\mathrm{1}}{\phi }_{\mathrm{1}}-{\phi }_{\mathrm{1}}{\phi }_{\mathrm{1}}{K}_{PK}=

K_{\mathrm{1}}{\phi }_{\mathrm{1}}\left(\mathrm{2}-K_{\mathrm{1}}{\phi }_{\mathrm{1}}\right)P_{f\mathrm{22}}-\mathrm{2}{K}_{\mathrm{2}}{\phi }_{\mathrm{1}}\left(\mathrm{1}-K_{\mathrm{1}}{\phi }_{\mathrm{1}}\right)P_{f\mathrm{12}}-

{\left(K_{\mathrm{2}}{\phi }_{\mathrm{1}}\right)}^{\mathrm{2}}{P}_{f\mathrm{11}}=\frac{{\phi }_{\mathrm{1}}^{\mathrm{2}}\left|{\mathit{P}}_f\right|\left(\mathit{\phi }{\mathit{P}}_f{\mathit{\phi }}^{\mathrm{T}}+\mathrm{2}{\lambda }_f\right)}{{\left(\mathit{\phi }{\mathit{P}}_f{\mathit{\phi }}^{\mathrm{T}}+{\lambda }_f\right)}^{\mathrm{2}}}>\mathrm{0}（23）

由此可知，若\mathit{P}矩阵由FFRLS算法中的协方差矩阵通过式（21）所示的方法构造得到，能使误差在单步修正过程中的\mathit{Q}矩阵半正定.

其次，需要保证在不同输入下\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)和\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}\right)不会同时等于0. 换言之，只有当输入相同时，\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)和\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}\right)才会同时等于0. 要证明上述结论，需先分析k时刻的误差修正方程，结合式（14）可得：

\mathrm{\Delta }{\widehat{\mathit{\theta }}}_k=\mathrm{\Delta }{\widehat{\mathit{\theta }}}_{k-\mathrm{1}}-{\mathit{K}}_k\left({\phi }_{\mathrm{1},k}\mathrm{\Delta }{\theta }_{\mathrm{1},k}+{\phi }_{\mathrm{2},k}\mathrm{\Delta }{\theta }_{\mathrm{2},k}\right)（24）

式中：

\mathrm{\Delta }{\mathit{\theta }}_k\in \{\mathrm{\Delta }{\mathit{\theta }}_k|{\mathit{\phi }}_k\mathrm{\Delta }{\mathit{\theta }}_k=\mathrm{0}\}

从式（24）中可以看到，参数误差在修正过程中表现为：估计值的误差\mathrm{\Delta }{\widehat{\mathit{\theta }}}_{k-\mathrm{1}}向目标值的误差\mathrm{\Delta }{\mathit{\theta }}_{\mathit{k}}所在的直线{\mathit{\phi }}_k\mathrm{\Delta }{\mathit{\theta }}_k=\mathrm{0}（后文简称目标直线）不断靠近（图2）. 经过k时刻的修正，能使\mathrm{\Delta }{\widehat{\mathit{\theta }}}_{\mathit{k}}到达目标直线附近. 若下一时刻的输入{\mathit{\phi }}_{k+\mathrm{1}}={\mathit{\phi }}_k，目标直线不变，那么目标值的误差也不变，即\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}=\mathrm{\Delta }{\mathit{\theta }}_k. 由式（17）可知，此时\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}\right)=\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)=\mathrm{0}，误差能量保持不变. 在电池参数辨识中，也有类似的结论^［10-11］——恒流输入时估计器无法进一步收敛；相反，若{\mathit{\phi }}_{k+\mathrm{1}}\ne {\mathit{\phi }}_k，目标直线变化，\mathrm{\Delta }{\widehat{\mathit{\theta }}}_k向目标值的误差\mathrm{\Delta }{\mathit{\theta }}_{\mathit{k}+\mathrm{1}}所在的直线{\mathit{\phi }}_{k+\mathrm{1}}\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}=\mathrm{0}不断靠近，此时\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}\right)\ne \dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right). 由此可知，在不同输入下，\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)和\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_{k+\mathrm{1}}\right)不会同时等于0.

从上述证明过程可知，在单步修正过程中会存在一种\mathit{P}矩阵使得\mathit{Q}矩阵半正定，即使输入\mathit{\phi }携带误差也不会对\mathit{Q}矩阵的正定性造成影响. 由此可知，虽然外部干扰影响系统矩阵，但系统稳定性不受其影响. 在连续两步修正过程中，不断变化的输入\mathit{\phi }是系统渐近稳定的必要条件.

2.3 基于图解法的渐近稳定性分析

为了进一步分析参数估计器的收敛过程，本节利用图解法来描述能量函数、误差修正方向和目标直线的几何特征. 由式（16）可知，能量函数的几何形状由{\mathit{P}}_k决定，为一系列同心椭圆，且误差\mathrm{\Delta }{\theta }_{\mathrm{1},k}、\mathrm{\Delta }{\theta }_{\mathrm{2},k}越大，能量椭圆面积越大，误差能量也就越大，反之亦然. 由式（14）可知，误差修正方向为两个时刻的误差点形成的向量，其斜率k_{\mathrm{m}}为：

k_{\mathrm{m}}=\frac{\mathrm{\Delta }{\widehat{\theta }}_{\mathrm{2},k}-\mathrm{\Delta }{\widehat{\theta }}_{\mathrm{2},k-\mathrm{1}}}{\mathrm{\Delta }{\widehat{\theta }}_{\mathrm{1},k}-\mathrm{\Delta }{\widehat{\theta }}_{\mathrm{1},k-\mathrm{1}}}=\frac{K_{\mathrm{2},k}}{K_{\mathrm{1},k}}（25）

当输入一定时，整个平面的误差修正方向都是相同的，且误差都沿着修正方向往目标直线{\mathit{\phi }}_k\mathrm{\Delta }{\mathit{\theta }}_k=\mathrm{0}靠近.

由式（20）可知，系统在单步修正过程中的\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)不恒小于0，即误差能量的变化是不确定的. 因此，本节将结合上述几何特征来定性分析误差能量的变化情况. 当{\mathit{\phi }}_k和{\mathit{P}}_k一定时，能量函数、误差修正方向和目标直线也随之被确定. 结合图3（a）进行分析，在k时刻，每一条误差修正方向总会与能量椭圆相切并形成一个切点，多个切点可以连成一条直线. 由切点连成的直线与目标直线将误差平面划分为3个区域（A、B、C）. 若初始误差\mathrm{\Delta }{\mathit{\theta }}_{k-\mathrm{1}}位于A、C两部分时，经过k时刻的修正，\mathrm{\Delta }{\mathit{\theta }}_{k-\mathrm{1}}将沿着误差的修正方向抵达\mathrm{\Delta }{\mathit{\theta }}_k. 显然，由\mathrm{\Delta }{\mathit{\theta }}_k形成的能量椭圆面积小于\mathrm{\Delta }{\mathit{\theta }}_{k-\mathrm{1}}，误差能量减小，误差渐近收敛；反之，误差渐近发散. 因此，根据构造的{\mathit{P}}_k和\mathrm{\Delta }{\mathit{\theta }}_{k-\mathrm{1}}所在位置的不同，误差有可能是渐近收敛的，也有可能是渐近发散的.

若改变{\mathit{P}}_k矩阵使得由切点连成的直线与目标直线所夹的B区域面积越小，误差就越有可能收敛，如图3（b）所示. 若在k时刻合理构造一种{\mathit{P}}_k\mathrm{矩阵}，使得由切点连成的直线与目标直线重合，如图3（c）所示，便可使除了目标直线上误差能量不变，其余位置的误差能量都减小，误差将最大可能地渐近收敛. 因此，通过构造{\mathit{P}}_k矩阵使单步修正过程中\dot{V}\left(\mathrm{\Delta }{\mathit{\theta }}_k\right)不大于0的方法是可行的，这与前文的分析是一致的.

结合误差能量的变化，误差收敛原理图解如图4所示. 当输入不断发生变化时，目标直线持续旋转，误差沿着修正方向朝目标直线靠近，误差能量逐渐减小，误差在趋近目标直线的过程中逐渐接近0.

3 动态工况下的仿真与实验验证

通过对估计器的误差收敛特性进行分析可知，变化的输入是估计器渐近收敛的必要条件，且误差能量在收敛的过程中会跟随修正方向不断减小. 本节选用2.6 Ah/3.6 V的三元锂电池设计实验与仿真，对上述理论分析进行验证. 实验采用新威CT-4004-5V100A型高性能电池检测系统作为电池监测设备. 在真实工况中无法实时获取极化参数的真值，为了更好地评价估计器的估计效果，需构建一个一阶ECM，用于生成仿真工况数据.

3.1 电池模型参数获取与实验验证

电池模型参数可通过连续充放电脉冲实验进行离线辨识获取初值，再由动态工况进行优化调整获取最终的阻抗参数，如图5所示. 其中，单个脉冲曲线包含10 s的1 C放电、40 s的静置以及10 s的1 C充电，通过对放电脉冲回弹的电压曲线进行拟合可得阻抗参数的初值. 间隔10% SoC重复上述流程，可得电池模型阻抗参数的初值. 随后采用US06工况进行动态测试，对阻抗参数进行优化调整获取最终电池阻抗参数及OCV-SoC曲线，分别如图6（a）和图6（b）所示. 调整的目的是使离线辨识的阻抗参数在动态工况中更为适用，进一步减小模型端电压误差.

调整的步骤为：

1）绘制US06工况下模型的端电压误差曲线.

2）设置电压误差阈值的绝对值为20 mV.

3）对于超过阈值的SoC区间，固定欧姆内阻、OCV-SoC曲线不变，通过调整对应SoC区间的极化电容与极化电阻，从而减小端电压误差.

经优化调整后，构建的一阶电池模型产生的仿真电压与实验电压的最大电压误差小于20 mV，如图6（c）所示. 为进一步验证该模型在其他工况下的适用性，采用动态压力测试工况（dynamic stress test， DST）、中国轻型车行驶工况（China light-duty vehicle test cycle，CLTC）对电池模型产生的端电压进行验证，如图7所示. 由图7可知，该电池模型在不同工况下产生的最大电压误差均小于30 mV，平均绝对误差小于6 mV，证明该模型对动态工况具有良好的适应性，可用于后续仿真验证.

3.2 估计器渐近收敛性的仿真验证

由电池模型生成动态工况下的电压数据作为电压真值. 为了模拟外部干扰，在电压真值中添加均值为2 mV、方差为3 mV的测量噪声^［25］作为电压测量值，并输入参数估计器中进行参数辨识. 同时，为了更好地对比不同工况下估计器的收敛程度，生成了两种典型动态工况，如图8所示.

图9展示了在不同工况中参数估计器参数误差变化轨迹. 由图9可知，在动态工况1中，经过一段时间的修正，估计器误差能收敛到0附近，表明该估计器在动态输入下具有良好的收敛性. 最后误差在0附近震荡，表明外部干扰的存在最终只会影响估计器的稳态特性. 由于工况2中的电流变化一开始较为平缓，其估计误差先是增加，而后才逐渐减小. 很明显，工况2中误差最终的收敛程度不如工况1，且p₁的误差显著小于0. 该结果验证了前文理论分析的正确性，即估计器在变化的输入下能渐近收敛. 同时，该结果进一步表明，工况的变化越剧烈，估计器的收敛性会越好.

以工况1为例对估计器的动态收敛过程进行详细分析，根据参数误差的动态变化轨迹，可绘制误差修正方向和误差能量椭圆图，如图10所示. 由图10可知，估计器的误差修正方向与理论分析的结果是一致的，整个误差平面的修正方向以目标直线为界同侧一致，异侧相反. 由于外部干扰的存在，修正方向无法完全平行于目标直线和能量椭圆等高线交点的切线，只能近似平行该切线，但外部干扰的存在并没有影响估计器的收敛性. 当输入不断变化时，误差沿着修正方向不断向目标直线靠近，直至误差减小到零附近. 工况2的收敛原理类似，故不再重复展示.

4 结 论

本文基于递归最小二乘法推导了锂电池一阶等效电路模型极化参数估计误差的状态方程，明确了系统矩阵是估计器稳定性的决定因素. 采用李雅普诺夫第二法对估计器的系统矩阵进行分析，得出持续变化的电流输入是估计器渐近收敛的必要条件. 提出了一种图解法对误差的修正方向和收敛轨迹进行了分析，分析表明误差在沿着修正方向朝目标直线靠近的过程中逐渐收敛，误差能量得到减小. 采用基于实验数据校准的电池模型对理论分析进行了仿真验证，结果表明估计器在持续变化的输入下能渐近收敛，且误差会跟随其修正方向不断收敛到零附近.

本文研究表明，估计器在持续变化的工况下，可展现出更好的收敛性. 后续的研究可基于这一结论开发自适应的工况筛选策略和算法，以提升估计器的收敛性能.

参考文献

[1] HAVANGI R.Adaptive robust unscented Kalman filter with recursive least square for state of charge estimation of batteries［J］.Electrical Engineering， 2022， 104（2）： 1001-1017．

[2] ZHENG Y J， CUI Y F， HAN X B，et al ．Lithium-ion battery capacity estimation based on open circuit voltage identification using the iteratively reweighted least squares at different aging levels［J］．Journal of Energy Storage，2021，44：103487．

[3] LIU Z， CHEN S H， JING B Q，et al ．Fractional variable-order calculus based state of charge estimation of Li-ion battery using dual fractional order Kalman filter［J］．Journal of Energy Storage，2022，52：104685．

[4] HOU J， LIU J W， CHEN F W，et al ．Robust lithium-ion state-of-charge and battery parameters joint estimation based on an enhanced adaptive unscented Kalman filter［J］．Energy，2023，271：126998．

[5] XING L K， LING L Y， GONG B，et al ．State-of-charge estimation for lithium-ion batteries using Kalman filters based on fractional-order models［J］．Connection Science，2022，34（1）：162-184．

[6] SHARMA A， FATHY H K. Fisher identifiability analysis for a periodically-excited equivalent-circuit lithium-ion battery model［C］//2014 American Control Conference， June 4-6，2014.Portland， OR， USA： IEEE， 2014： 274-280．

[7] ROTHENBERGER M J， ANSTROM J， BRENNAN S，et al ．Maximizing parameter identifiability of an equivalent-circuit battery model using optimal periodic input shaping［C］//ASME 2014 Dynamic Systems and Control Conference，October 22-24，2014. San Antonio，Texas，USA： American Society of Mechanical Engineers，2014.

[8] LIN X F， STEFANOPOULOU A G ．Analytic bound on accuracy of battery state and parameter estimation［J］．Journal of the Electrochemical Society，2015，162（9）：A1879-A1891．

[9] LIN X F ．Analytic analysis of the data-dependent estimation accuracy of battery equivalent circuit dynamics［J］．IEEE Control Systems Letters，2017，1（2）：304-309．

[10] SONG Z Y， HOFMANN H， LIN X F，et al ．Parameter identification of lithium-ion battery pack for different applications based on Cramer-Rao bound analysis and experimental study［J］．Applied Energy，2018，231：1307-1318．

[11] SONG Z Y， WANG H， HOU J，et al ．Combined state and parameter estimation of lithium-ion battery with active current injection［J］．IEEE Transactions on Power Electronics， 2020， 35（4）： 4439-4447．

[12] SONG Z Y， HOU J， LI X F，et al ．The sequential algorithm for combined state of charge and state of health estimation of lithium-ion battery based on active current injection［J］．Energy，2020，193：116732．

[13] LIN X F ．Theoretical analysis of battery SOC estimation errors under sensor bias and variance［J］．IEEE Transactions on Industrial Electronics，2018，65（9）：7138-7148．

[14] SHEN P， OUYANG M G， HAN X B，et al. Error analysis of the model-based state-of-charge observer for lithium-ion batteries［J］．IEEE Transactions on Vehicular Technology，2018，67（9）： 8055-8064．

[15] ZHOU W， HUANG R J， LIU K，et al ．A novel interval-based approach for quantifying practical parameter identifiability of a lithium-ion battery model［J］．International Journal of Energy Research， 2020， 44（5）： 3558-3573．

[16] ANDERSSON M， STREB M， KO J Y，et al ．Parametrization of physics-based battery models from input-output data：a review of methodology and current research［J］．Journal of Power Sources，2022，521：230859．

[17] NASERI F， SCHALTZ E， STROE D I，et al ．An enhanced equivalent circuit model with real-time parameter identification for battery state-of-charge estimation［J］．IEEE Transactions on Industrial Electronics， 2022， 69（4）： 3743-3751．

[18] HUA X， ZHANG C， OFFER G ．Finding a better fit for lithium ion batteries：a simple，novel，load dependent，modified equivalent circuit model and parameterization method［J］．Journal of Power Sources，2021，484：229117．

[19] SHU X， CHEN Z， SHEN J W，et al ．State of charge estimation for lithium-ion battery based on hybrid compensation modeling and adaptive H-infinity filter［J］. IEEE Transactions on Transportation Electrification， 2023， 9（1）： 945-957．

[20] WU J J， FANG C， JIN Z Y，et al ．A multi-scale fractional-order dual unscented Kalman filter based parameter and state of charge joint estimation method of lithium-ion battery［J］．Journal of Energy Storage，2022，50：104666．

[21] CHEN Z G， ZHOU J X， ZHOU F，et al ．State-of-charge estimation of lithium-ion batteries based on improved H infinity filter algorithm and its novel equalization method［J］．Journal of Cleaner Production，2021，290：125180．

[22] 王长江，姜涛，陈厚合，等 ．基于相位校正李雅普诺夫指数的电力系统暂态电压稳定评估［J］．电工技术学报，2021，36（15）：3221-3236．

[23] WANG C J， JIANG T， CHEN H H，et al ．Transient voltage stability assessment of power systems based on phase correction maximum Lyapunov exponent［J］．Transactions of China Electrotechnical Society，2021，36（15）：3221-3236．（in Chinese）

[24] ZHOU P， HU X K， ZHU Z G，et al ．What is the most suitable Lyapunov function？［J］．Chaos，Solitons & Fractals，2021，150：111154.

[25] WANG R N， WANG Q W， LIU L S ．Solving a system of sylvester-like quaternion matrix equations［J］．Symmetry，2022，14（5）：1056．

[26] CUI Z R， CUI N X， WANG C Y，et al ．A robust online parameter identification method for lithium-ion battery model under asynchronous sampling and noise interference［J］. IEEE Transactions on Industrial Electronics，2021，68（10）：9550-9560．

基金资助

湖南省自然科学基金资助项目（2022JJ30158）; 国家自然科学基金资助项目（51705139）