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

文章字数有点多，耐心不够的谨慎点击阅读。
下期继续讲解在线参数辨识方法。

背景

最近又在开始重新梳理锂电池建模仿真与SOC算法的知识框架，因此收集了大量文献资料，然后挑选几份相对靠谱的文献，挨个啃读解析。此篇推文选了其中一篇比较有代表性的论文《基于RLS方法的磷酸铁锂电池模型辨识及SOC估计策略研究》，作者为西南交通大学-郑卫同学。论文中实际用到的电池型号如下：

此外部分代码参考了下述链接的内容，在此一并列出。
下图为论文中识别出的一阶模型参数。文末附上使用论文中辨识后的参数做的simulink模型和相关代码。

• SOC-OCV
  注：论文中有一处笔误，k3系数应该是-1219.40。
  
• 阻容参数
  

模型简介

下图摘自论文内容，此为相对简单的等效模型。

往期的文章已经对上述模型有过解释与实际建模，内含simulink模型，可以看看：
【Matlab】如何使用Simulink进行电池建模和SOC预估
链接: link

数据收集

在开始这部分内容之前，引用论文中的几个基本概念：

首先对电池状态有如下定义：
（1）CC-CV充电模式：CC-CV（Constant Current-Constant
Voltage）表示恒流恒压充电，先以恒流充电至端电压到上截止电压Vmax，再以Vmax为充电电压进行恒压充电，当电流下降至截止电流时停止充电。
（2）满电状态：在一定温度下，电池恒流充电至上截止电压后，并转为恒压充电，至截止电流时停止充电后静置若干小时，电压趋于平稳，此时电池的状态被定义为满电状态，即SOC=1。
（3）放空状态：电池放电至下截止电压后，停止放电，静置若干小时，电压趋于平稳，此时电池的状态定义为放空状态，即SOC=0。
确定各状态定义后，进行的电池特性实验有：

1. 最大可用容量实验

参考《电动汽车用电池管理系统技术条件》对电池进行充放电，以标定其当前条件下的最大可用容量。电池容量会随着电池老化、工作温度等因素不断波动，因此有必要测试不同温度和不同循环次数情况下的最大可用容量。以标准电流使用恒流恒压
（CC-CV）方式将磷酸铁锂电池充满电，经充分静置后再恒流放电至下截止电压，连续测试三次，三次的结果与均值的偏差不得超过2%，否则此次测得的结果无效。此论文中实现选用1/3C（1C基于额定容量）为恒流段的充放电倍率，该试验的意义在于明确不同温度下的最大可用容量的分布情况。实验记录的电流电压数据如图2-2、2-3所示。

2. 开路电压实验

该实验的目的是获取不同温度、不同电流倍率情况下的电池OCV-SOC的函数关系曲线。每种电池体系都有自身特定的OCV曲线，与电池当前的SOC存在一定关系。
在同一温度条件下，电池的OCV也会受到老化因素的影响，也可以作为估算电池健康状态（SOH）的依据。OCV-SOC曲线对等效电路模型的参数辨识具有重要意义。OCV
分为充电OCV和放电OCV两组值。具体步骤为：
①作充电曲线时，以标准电流将电池放电至截止电压，静置5h后，测量端电压值，视为SOC=0时的开路电压值；作放电曲线时，先以标准电流充电至满电状态，静置5h，该值视为SOC=1时的开路电压值。
②作充电曲线时，对动力电池充电，充电量为5%的最大可用容量，直至截止电流，静置5h后测量电池端电压值；作放电曲线时，放电量为5%的最大可用容量，静置5h后测量电池端电压值，直至截止电压。
③ 跳到步骤②，循环进行步骤②和③直到磷酸铁锂电池达到截止条件。 实验记录的电流电压数据如图2-4至2-7所示。

小编备注：对于三元类电池，5%的DOD是够用的，但是对于磷酸铁锂电池而言，这个间隔大了一些，建议控制在≤2%，后面其他期的推文会讲到。


论文中采用一元八次方程进行SOC-OCV拟合，从工程应用角度不可取。建议使用线性插值，计算量小，而且只要SOC间隔控制在2%，绝对比高次方拟合效果好。

3. 混合动力脉冲特性实验

混合动力脉冲特性（HPPC）实验是使用连续的脉冲激励序列对电池进行充放电，以获得电池在不同SOC状态下的特性。结合动力电池建模理论，HPPC测试数据通常被用来进行模型参数的辨识，获得不同SOC点的模型参数表。具体测试步骤如下：
① 将环境温度设置到x℃，使用标准电流以CC-CV充电方式将待测试动力电池充至满电。
② 静置5h，使动力电池接近于平衡状态。
③ 加载混合脉冲电流，之后对电池进行一段时间的恒流放电，保证前后两次脉冲激励时电池的SOC相差5%，然后静置1h。
④ 重复步骤③，直到电池达到放电截止电压，可获得电池在SOC为100%、 95%、……、5%情况下的数据。

小编备注：同上，实际工程应用中，建议控制在≤2%，后面其他期的推文会讲到。另外，如果只是一阶模型，仅用开路电压的过程测试数据亦可以完成参数识别，HPPC的测试数据可以用于已识别参数的适应性和偏差验证。

离线辨识

作为一名专业的BMS软件工程师，我看到的大多数论文内容，缺乏工程化应用思维和更深层次的理解。基于这一点考量，所以由我来跟大家详细讲解一下。

本图摘自论文。

一阶的离线辨识相对比较简单，我们以开路电压过程数据即可完成识别。一起再更加深入地理解一遍，一阶模型的各项参数如何辨识：

1. a点的电压由于经过长时间静置，可以认为是此刻SOC对应的接近平衡态的电动势，即论文中提到的OCV；
2. 同样的e点也经过了长时间静置，也可以认为是另外一个SOC点的OCV；
3. (Ua - Ub) / I = R0为SOC-a对应的欧姆内阻，(Ud - Uc) / I = R0为SOC-c对应欧姆内阻。这里重点讲三点：

• 首先为什么可以这样算呢？我们知道电容的特性就是电容两端的电压在电流变化瞬间是不会突变的，但是纯电阻在电流流过瞬间就开始有压降，只要时间够短，我们可以认为从静置到有电流瞬间，端电压的变化完全由模型中的R0压降造成的，即ΔV=I*R0。
• 其次，在实际工程应用中，因为不同测试设备的差异（大部分电芯分容设备能做到100ms采样周期），a-b点采样时间不一定能做到毫秒级别，如果是桌面端的应用，能短则短，因为这样才能最为准确的识别出纯欧姆阻抗，时间长了会引入极化的影响；
• 如果是需要将离线辨识的参数，应用到嵌入式产品端，则最好与你嵌入式产品中，SOC算法运行的周期一致。例如，你的BMS产品中算法的运行频率是200ms，那你在离线测试电芯数据时也设置为200ms间隔去收集数据；此外BMS端的电芯电压和电流要尽可能同步采集，否则会影响你的电池模型算法递推的准确性；

4. d到e为极化R1C1的效果，那拟合这段曲线的底层逻辑又是什么？请往下看。电池充电结束后，由于极化效应，电芯电压会逐渐回落至接近平衡电动势，此时从曲线上看，就类似于电容的放电效应
   
   识别R1C1有两种方法：

• 方法1-手动找时间点：使用Ude段0.95压差处的时间作为3R1C1值。为什么是3倍的R1C1？我们知道电容的时间常数τ=RC。
  3个时间常数（3τ）后：充电时电压将接近95%（1 - e^(-3) ≈ 0.95），放电时则降至初始电压的5%（e^(-3) ≈ 0.05），从前人的经验来看，选用3τ时间拟合的效果比较好。
  那具体怎么选取呢？
  方法为以d点电压Ud为电压和时间基准，找到d-e线中，电压=Ud+（Ue-Ud）*0.95的相对时间T，那么3R1C1=T，R1C1即为T/3。再根据b-c时刻的电流I，算出极化电阻R1=(Ue-Ud)/I，再根据C1=T/3 / R1即可算出所有参数。
• 方法2-将数据丢给matlab，让它来帮忙拟合（直接用可视化cftool或者用现成的拟合库函数）：从d点开始，因为没有了电流，所以此刻:
• Ut = Uocv-Urc=Uocv-（Uocv-Ud）e^(-t/R1*C1)，其中Ut为实际设备采样到的电压，Uocv为放电结束时刻SOC值对应的OCV，t为设备采样电压数据的相对时间
• 将R1*C1作为一个整体时间常数进行辨识，那么R1和C1怎么单独算出来呢？答案藏在下述关系式中。在电流截止前的最后一刻,我们假想取一个微小时刻，那么dUc趋向于0，电容Uc产生的原因就是I(t)的影响，下述的第二个式子就变为了
• I = Uc / R1，其中Uc可以取静置最后时刻的电压，或者其他方式得到的该SOC态下对应的OCV，此处我们选择静置最后时刻的电压。R1有了，时间常数也有了，那么C1自然也可以算出。
  下图为论文中，作者通过辨识得到的参数，我规整了下。需要注意的是单位。我特意将毫欧和欧，千法拉和法拉列在一起。实际运算和建模都是用国际制标准单位，即V，Ω，F，时间单位为s。
  

对应模型

因为没有论文中的实际数据，所以使用论文中的参数，构建了一个电池模型，并通过设定电流进行模拟数据生成，然后再放入离线参数识别代码中，辨识出当前时刻的参数值。

下载链接：
https://download.csdn.net/download/wanff1204/89499217

对应代码

通过以下代码实际辨识出来的参数如下：
RC = 7154.10, R = 0.001661, C = 4306554.1
其中RC时间常数值与原始数据相当，R与C因为充电电流与时长的差异，会有差别。供读者自行优化。

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

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

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

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

% 步骤1: 运行模型，并提取所需数据用于其他步骤(模型的数据已经通过设置，输入到工作区)
modelname = 'Battery.slx'; 
sim(modelname);

Ut   = ScopeData4.signals(2).values; % Scope中的数据索引似乎从1开始
Ibat = ScopeData4.signals(3).values; 

cutOffIndex = 1;
% 截取电流截至后的电压
for i = 2 :  size(Ibat, 1)
    if Ibat(i - 1) > 0 &&  Ibat(i) <= 0
        cutOffIndex = i + 1;
        break;
    end
end

UtNoneCurrent = Ut(cutOffIndex : size(Ut, 1));
timeNoneCurrent = ScopeData4.time(cutOffIndex : size(Ut, 1)) - ScopeData4.time(cutOffIndex - 1); % 获取时间向量
startSoc = ScopeData4.signals(1).values(cutOffIndex);

% 获取电流，用于计算R1
dischgCurrent = ScopeData4.signals(3).values(cutOffIndex - 2);   

% 基于SOC获取OCV -- 见论文的拟合公式
ocv = -95.82*(startSoc^8)+549.26*(startSoc^7)-1219.4*startSoc^6+1387.01*startSoc^5-883.38*startSoc^4+320.4*startSoc^3-64.45*startSoc^2+6.89*startSoc+2.91;

% 步骤2: 定义指数模型函数-Uocv-（Uocv-Ud）e^(-t/R1*C1)
expModel = @(a,x) ocv - (ocv - UtNoneCurrent(1)) * exp(-x / a);


% 步骤3: 使用fitnlm进行拟合
fittedModel = fitnlm(timeNoneCurrent, UtNoneCurrent, expModel, 1800); 

% 步骤4: 获取拟合参数及置信区间
tau_rc = fittedModel.Coefficients.Estimate; % 拟合参数
r1 = (UtNoneCurrent(end) - UtNoneCurrent(1) ) / dischgCurrent;
c1 = tau_rc / r1;
fprintf('RC = %.2f, R = %f, C = %.1f\r\n', tau_rc, r1, c1);

ci = coefCI(fittedModel); % 置信区间

% 步骤5: 绘制原始数据与拟合曲线
figure;
plot(timeNoneCurrent, UtNoneCurrent, 'o', 'MarkerFaceColor', 'g', 'MarkerEdgeColor', 'k', 'LineWidth', 1.5, 'DisplayName', '原始静置曲线');
hold on;
xx = linspace(min(timeNoneCurrent), max(timeNoneCurrent)); % 生成x轴范围的等间距点用于绘图
yy = expModel(tau_rc, xx); % 计算拟合曲线上的y值
plot(xx, yy, 'r-', 'LineWidth', 2, 'DisplayName', '拟合后的静置曲线');
legend('show');
title('Exponential Fit of Data');
xlabel('X data');
ylabel('Y data');

% 步骤6: 绘制拟合曲线与原始数据的差值
figure;
residuals = UtNoneCurrent - expModel(tau_rc, timeNoneCurrent); % 残差（拟合值与实际值之差）
plot(timeNoneCurrent, residuals, 'bx', 'MarkerFaceColor', 'b', 'LineWidth', 1.5, 'DisplayName', 'Residuals');
xlabel('X data');
ylabel('Residuals');
title('Residual Plot');

总结

文章详细介绍了锂电池一阶模型离线辨识方法。此处离线辨识主要通过开路电压过程数据来完成，包括欧姆内阻R0和极化电阻R1、电容C1的辨识。文章提供了两种辨识方法：手动找时间点和使用MATLAB进行数据拟合。最后，通过模拟数据生成和离线参数识别代码，辨识出当前时刻的参数值，并附上MATLAB代码供读者参考。

下期预告

下期会继续基于该篇论文，详细展开说明如何在线识别上述参数。