一篇就够了，为你答疑解惑：锂电池一阶模型-在线参数辨识（附代码）

背景

锂电池一阶戴维南等效模型的基础知识和离线辨识方法，已经在上一期非常详细地讲解了一轮（上期文章请戳此处），本期继续讲解一下如何进行在线辨识。
此篇推文继续使用论文《基于RLS方法的磷酸铁锂电池模型辨识及SOC估计策略研究》中的方法，作者为西南交通大学-郑卫同学；用到的模型在上一期可以获取。

在线 VS 离线 参数辨识

• 首先，我们来详细说说在线辨识与离线辨识的基础概念。

两者本质的区别为，辨识的过程是否在系统正常运行过程中实时或周期性地估计参数。

离线辨识是在系统不运行或实验条件下进行的参数估计过程。它通常涉及收集系统的静态或历史数据，然后通过复杂的计算和数据分析来确定模型参数。离线辨识的优势在于能够使用更复杂的算法和模型，因为它不受实时性约束，可以花更多的时间和计算资源来提高辨识精度。这种方法常用于研发阶段，帮助工程师理解系统的工作原理，优化设计，或建立精确的模型用于仿真。

在线辨识是在系统正常运行过程中实时或周期性地估计参数的技术。它通过实时监测系统的输入输出信号（如电流、电压、速度等），使用快速算法实时更新模型参数，以适应系统随时间和环境变化的动态特性。在线辨识强调实时性和鲁棒性，适合需要动态调整控制策略的应用场景。

• 其次我们想一想，为什么需要在线辨识？
  锂电池特性的影响因素很多，如焊接工艺，温度，电流倍率和电池本身的劣化等。上述在线辨识的基本概念，已经体现了它的优势-可以适应系统随时间和环境变化的动态特性。因此，纯粹通过离线辨识多个参数组，不仅耗费大量时间和金钱，而且在实际使用中会因为单个或者多个因素叠加，导致离线参数实用性差，从而导致模型算法难以精准表征电池状态。
• 那是不是离线辨识就没有意义了？
  并不是。通过离线辨识，我们可以提前知道该电池参数的大致范围，可以用于在线辨识开始前的参数初始化以及在线辨识结果的上下限判断。

递推最小二乘法

论文中在线辨识的方法为带遗忘因子的递推最小二乘法，以下内容摘自论文

对电池系统而言，根据系统辨识的理论建立数学模型，利用实验的输入输出数据 辨识所需要的参数，就是电池模型的参数辨识。本文选用递推最小二乘法（Recursive Least Squares Method，RLS）进行模型的参数辨识。递推最小二乘法是一种线性无偏估算方法，使观测误差的平方和达到极小，主要特点是适用性非常广，线性系统和非线性系统均适用；动态系统和静态系统也均适用；离线估算和在线估算均可用该方法；在随机环境下，观测得到的数据即使没有相关的概率统计信息，使用递推最小二乘法 的估算结果仍具有良好的统计特性。
递推最小二乘法（RLS）是在最小二乘法（LS）的基础上发展来的一种估算方法，是对一次最小二乘法的估算值修正后的一种在线辨识方法。当辨识系统启动后，每取得一次新的观测数据，就在前一次估算结果的基础上，引入此次新的观测数据对前一次估算的结果作修正，进而递推得出新的参数估值。如此算法进行下去，随着新的观测数据不断被引入进来，就实现了一次接一次的参数辨识，参数估算值直至收敛到合格精度为止。
递推最小二乘法的基本思想：新估算值=旧估算值+修正项 基本计算公式为：

RLS方法具有无限记忆长度，随着递推运算过程的不断进行，旧的数据越来越多导致出现数据饱和的现象，使得递推结果逐渐失效。因此为避免这种情况，引入遗忘因子，当时，就退变为普通的递推最小二乘法，越小表示遗忘速度越快，跟
踪能力越强，但波动也越大，一般取[0.95 1]。

这里需要说明一点，递推最小二乘法本身没有在线或者离线属性，取决于我们实际如何用它。你把它放在的嵌入式代码中跑，每新采集一组数据时算一次，就成了在线辨识；你已经有了数据，再去按步骤递推就是离线辨识。

一阶戴维南Z域离散表达式

在使用递推前，首先要推导出参数矩阵与数据矩阵的离散关系式。
基础表达式还是下面这个，只不过咱们需要对它进行离散化。


注：这里面涉及到一个电流哪个方向为正方向问题，论文此处两个笔误，以免误导：
- 上图的uc(s)与uc(t)是两个概念，不是同一个，但为了和论文统一，我接下来的uc统一指的是内部压降，而非电容的电压
- 若以放电方向为正方向，则uc = E(k) - V(k) ,其中V(k)为端电压，即实际采集到的电池电压，所以论文中的式子3-2是错误的

但是这个式子怎么来的呢？对于没有学过信号处理或者控制理论的人来说，在看类似论文时，会一头雾水。如果只想应用结论的，则可以跳开此部分内容，直接看本章最后的内容。
在s域分析中，电阻和电容并联的等效阻抗可以用复频域的运算来表示。在时域中，电阻的阻抗是纯实数，而电容的阻抗是纯虚数并且与频率有关。将它们转换到s域（拉普拉斯变换的领域），电阻R的阻抗仍然是R（因为其即时响应特性），而电容C的阻抗变为1/cs.
因此，电阻R与电容C并联的总阻抗Z可以用以下公式表示：
$$Z=\frac{1}{\frac{1}{R}+\frac{1}{\frac{1}{sC}}}$$
简化后得到：
$$Z=\frac{R\cdot \frac{1}{sC}}{R+\frac{1}{sC}}$$
$$Z=\frac{R}{1+RsC}=\frac{R}{1+\tau s}$$
传递函数一般都会写成 输出/输入 的形式，因此再加上R0与R1C1串联，则可以得到上述s域的表达式。
但是s域主要应用于连续时间信号和系统的分析，要放在递推中，还需要转为z域形式。z域是通过Z变换（Z Transform）从离散时间域转换而来，其中z也是一个复数变量，主要关注的是离散时间序列的频域特性。z变换适用于数字信号处理、数字控制以及任何涉及采样数据系统的分析，比如数字滤波器设计、数字信号处理器(DSP)算法开发等。
从s域转到z域，通常使用双线性变换，即通过下述方法，将s替换。
$$s=\frac{2}{T}\ast \frac{1-z^{-}1}{1+z^{-}1}或\frac{2}{T}\ast \frac{z-1}{z+1}$$
两者是等效的，前者是后者分子分母同时除了z而已，其中T为数据采样周期。则最终可以得到形如这样的标准的z域传递函数形式：

具体推导过程可以自行再推一次，下图仅供参考。需要注明一点的是Z^-1乘上I(k)表示I(k)前一个采样周期的值，即为递推过程的I(k-1)值。
最终得到的：
$$a0=T+2\ast \tau =T+2\ast R1\ast C1$$
$$a1=T-2\ast \tau =T-2\ast R1\ast C1$$
$$b0=T\ast (R0+R1)+2\ast R0\ast \tau$$
$$b1=T\ast (R0+R1)-2\ast R0\ast \tau$$
再令
$$c1=-\frac{a1}{a0}$$ $$c2=\frac{b0}{a0}$$$$c3=\frac{b1}{a0}$$
可得
$$Uc(k)=c1\ast Uc(k-1)+c2\ast I(k)+c3\ast I(k-1)$$

注：上图推导过程将c1的符号位忘了移到Uc(k-1)侧了，懒得重新写了
由c1，c2，c3，我们可以逆推R0，R1，C1的表达式

回到刚才的递推最小二乘法表达式，采用矩阵的方式进行表达

论文部分符号的含义没有解释清楚，在此一并补充说明：

• K-增益矩阵，影响了新数据对参数估计的影响程度。可以不用给初值，每次迭代过程中产生
• P-协方差矩阵，需要根据实际情况给定初值。其对角线上的元素表示各个参数估计值的不确定性的平方，这其实就是协方差的定义。初值的选择会影响算法的收敛速度和参数估计的准确性。一般来说，初值越大，说明预设的参数与真实参数的偏差比较大，则算法在初始阶段对观测数据的响应越灵敏，但也可能导致算法在收敛过程中出现较大的波动。
• z(t) - 有些文献写的是y(k)，离散化表达式的左值，对应于我们的模型就是Uc(k)观测值=电池端电压-当前SOC态下Uoc
  则根据前面的推导，套入到RLS递推，有如下表达式
  $$\theta =[c1c2c3]$$$$\varphi =[-Uc(k-1)I(k)I(k-1){]}^T$$
  接下来就会用到我们上一篇离线参数辨识的结果了。而且我们需要给协方差矩阵P一个初值，我们认为经过了离线辨识，参数已经比较准了，可以给个对角元素均为1的3*3矩阵，同时根据表达式，可以给出c1-c3的初值。代码如下：

% 执行该脚本前，请确认simulink模型与该脚本文件是否在同一路径

clear all;  % 清除工作区中的所有变量
close all; % 关闭所有已打开的图形窗口
clc;           % 清空命令窗口的内容

% 打印脚本开始的信息（可选）
fprintf('Script started.\n');

% 这里开始编写你的MATLAB脚本内容...

% 步骤1: 运行模型，并提取所需数据用于其他步骤(模型的数据已经通过设置，输入到工作区)
% 具体方法为用Scope记录模型仿真过程数据，然后配置为输出到工作区
modelname = 'Battery.slx';
sim(modelname);

% 步骤2: 根据之前离线辨识的结果，确定参数给定初值，注意：递推用的是c1-c3
% 获取模型的运行步长，用于递推的输入，因为我们的模型是定步长.，所以可以用如下接口获取
T = 0.1;%str2double(get_param('Battery.slx', 'FixedStep'));
R1C1_init = 7154.10;
R0_init = 0.001661;
C1_init = 4306554.1;
R1_init = R1C1_init / C1_init;

c1_init = (T - 2 * R1C1_init) / (T + 2 * R1C1_init);
c2_init = ((R0_init + R1_init) * T + 2 * R1C1_init) / (T + 2 * R1C1_init);
c3_init = ((R0_init + R1_init) * T - 2 * R1C1_init) / (T + 2 * R1C1_init);

theta = [c1_init c2_init c3_init]'; %参数向量初值
phi = zeros(1, 3);                 %数据向量初值
P = 1*eye(3);                      %协方差矩阵
K = zeros(3,1);                    %增益矩阵
lambda = 0.99;

% 步骤3: 递推
% 数据准备
Ut   = ScopeData4.signals(2).values; % Scope中的数据索引似乎从1开始
Ibat = ScopeData4.signals(3).values; 
Soc = ScopeData4.signals(1).values;
time_val = ScopeData4.time;

% 见博客说明，先处理z(k) = Ut(k) - Uocv(k)
% 基于SOC获取OCV -- 见论文的拟合公式
ocv = -95.82*Soc.^8  + 549.26*Soc.^7 ...
      -1219.4*Soc.^6 + 1387.01*Soc.^5 ...
	  -883.38*Soc.^4 + 320.4*Soc.^3 ...
	  -64.45*Soc.^2  + 6.89*Soc + 2.91;
Uc = Ut - ocv;
R0 = zeros(1, size(Uc, 1) - 1);
R1 = zeros(1, size(Uc, 1) - 1);
C1 = zeros(1, size(Uc, 1) - 1);

% 递推主体
for i = 2 :  size(Uc, 1)
	Phi   = [-Uc(i - 1) Ibat(i) Ibat(i - 1)];
	K     = P * Phi' / (Phi * P * Phi' + lambda);
	theta = theta + K * (Uc(i) - Phi * theta);
	P     = (eye(3) - K*Phi) * P / lambda;
	
	% 基于递推结果，反推R0, R1, C1
	c1_k = theta(1);
	c2_k = theta(2);
	c3_k = theta(3);

	R0(i) = (c2_k - c3_k)/(1 - c1_k);
	R1(i) = (c2_k + c3_k) / (1 + c1_k) - R0(i);
	C1(i) = T * (1 - c1_k) / (2 * R1(i) * (1 + c1_k));
end


% 步骤4: 画图查看递推过程情况图
figure;
subplot(2,2,1); % 创建一个2行2列的子图，并激活第1个子图
plot(time_val, Ut, 'o', 'MarkerFaceColor', 'g', 'MarkerEdgeColor', 'k', 'LineWidth', 1, 'DisplayName', '电压曲线');
title('原始曲线');
hold on;
plot(time_val, Ibat, 'r-', 'LineWidth', 1, 'DisplayName', '电流曲线');
legend('show');
xlabel('Time (s)');
ylabel('Voltage/Current');
hold off; % 关闭保持状态，避免下一个子图影响当前子图

subplot(2,2,2); % 激活第2个子图,显示R0
plot(time_val, R0 * 1000, 'bx', 'MarkerFaceColor', 'b', 'LineWidth', 1, 'DisplayName', 'R0');
legend('show');
xlabel('Time (s)');
ylabel('R0(mO)');
title('R0-递推');
hold off; % 关闭保持状态，避免下一个子图影响当前子图

subplot(2,2,3); % 激活第3个子图,显示R1
plot(time_val, R1 * 1000, 'bx', 'MarkerFaceColor', 'b', 'LineWidth', 1, 'DisplayName', 'R1');
legend('show');
xlabel('Time (s)');
ylabel('R1(mO)');
title('R1-递推');
hold off; % 关闭保持状态，避免下一个子图影响当前子图

subplot(2,2,4); % 激活第4个子图,显示R0
plot(time_val, C1, 'bx', 'MarkerFaceColor', 'b', 'LineWidth', 1, 'DisplayName', 'C0');
legend('show');
xlabel('Time (s)');
ylabel('C1');
title('C1-递推');

R0= 0.0049， R1 =0.0997， C1 =7.1732e+04
最终辨识结果如上