1. 电池SOC估算的背景与挑战
在电池管理系统(BMS)中,荷电状态(State of Charge, SOC)的准确估算一直是个核心难题。SOC就像是我们手机上的电量百分比,但实际要复杂得多——它直接关系到电池的安全使用、寿命管理和能量调度。我做过不少电池项目,最头疼的就是这个SOC估算不准的问题。
传统方法里,安时积分法(Ah-counting)是最直接的,就像用沙漏计时一样累计进出电池的电量。但这个方法有个致命缺陷:误差会随时间累积。就像用有误差的钟表,时间越长偏差越大。实际项目中,仅用安时积分法,一周后SOC误差可能超过10%,这对电动汽车等应用是完全不可接受的。
于是研究者们引入了卡尔曼滤波家族的方法。就像用GPS修正计步器的误差一样,EKF(扩展卡尔曼滤波)和UEKF(无迹扩展卡尔曼滤波)通过电压等观测值来修正安时积分的累积误差。但这里面的门道很多——电池的非线性特性、模型精度、噪声处理等都会显著影响结果。
2. 基础原理:从安时积分到卡尔曼滤波
2.1 安时积分法的本质与局限
安时积分的数学表达很简单:
SOC(t) = SOC₀ + (1/Qₙ) ∫ηi(τ)dτ
其中Qₙ是额定容量,η是库仑效率,i是电流。我在Matlab中实现这个公式时,关键要注意:
- 采样频率要足够高(至少1Hz)
- 电流传感器的零漂必须校准
- 初始SOC₀的准确性至关重要
matlab复制% 安时积分法简单实现示例
function soc = ah_integration(current, dt, soc_init, Qn)
persistent accumulated_charge;
if isempty(accumulated_charge)
accumulated_charge = soc_init * Qn;
end
accumulated_charge = accumulated_charge + current * dt;
soc = accumulated_charge / Qn;
end
但实际测试会发现,即使电流测量非常精确,容量衰减和温度影响也会导致误差。这就是为什么需要引入模型修正。
2.2 电池建模:二阶RC等效电路
在EKF/UEKF中,我们需要电池的数学模型。二阶RC模型是最常用的选择,它用电路元件模拟电池动态:
- 开路电压(OCV):理想电压源,与SOC有确定关系
- 欧姆内阻(R₀):瞬时电压降
- 两个RC环节(R₁C₁, R₂C₂):模拟极化效应
matlab复制% 二阶RC模型参数示例(锂离子电池)
params.R0 = 0.01; % Ohm
params.R1 = 0.005; % Ohm
params.C1 = 2000; % F
params.R2 = 0.01; % Ohm
params.C2 = 5000; % F
建立准确的OCV-SOC关系曲线是关键。我通常会在不同温度下进行充放电测试,用多项式拟合:
matlab复制% OCV-SOC关系拟合
soc_points = 0:0.1:1;
ocv_points = [3.0 3.3 3.45 3.5 3.55 3.6 3.65 3.7 3.8 3.9 4.1];
p = polyfit(soc_points, ocv_points, 5);
2.3 卡尔曼滤波框架概述
卡尔曼滤波的核心思想是"预测-修正"循环:
- 预测步:基于模型预测下一时刻状态
- 更新步:用实际测量值修正预测
对于非线性系统,EKF通过雅可比矩阵线性化,而UEKF采用sigma点采样,后者在处理强非线性时通常表现更好。
3. EKF实现细节与Matlab代码解析
3.1 状态空间方程建立
我们需要定义状态变量和观测变量。典型选择:
状态变量x = [SOC, V₁, V₂]ᵀ
观测变量z = 端电压
状态方程:
xₖ = f(xₖ₋₁, iₖ₋₁) + wₖ₋₁
zₖ = h(xₖ, iₖ) + vₖ
其中w和v是过程噪声和观测噪声。
matlab复制% 状态转移函数
function x_next = stateFcn(x, current, params, dt)
soc = x(1);
v1 = x(2);
v2 = x(3);
% SOC更新
soc_next = soc - (current * dt) / (params.Qn * 3600);
% RC环节电压更新
v1_next = exp(-dt/(params.R1*params.C1)) * v1 + ...
params.R1*(1-exp(-dt/(params.R1*params.C1))) * current;
v2_next = exp(-dt/(params.R2*params.C2)) * v2 + ...
params.R2*(1-exp(-dt/(params.R2*params.C2))) * current;
x_next = [soc_next; v1_next; v2_next];
end
% 观测函数
function voltage = measurementFcn(x, current, params)
soc = x(1);
v1 = x(2);
v2 = x(3);
ocv = polyval(params.ocv_poly, soc);
voltage = ocv - current*params.R0 - v1 - v2;
end
3.2 EKF核心算法实现
EKF的关键步骤包括:
- 状态预测
- 协方差预测
- 卡尔曼增益计算
- 状态更新
- 协方差更新
matlab复制% EKF主循环
for k = 2:length(t)
% 预测步
x_pred = stateFcn(x_est(:,k-1), current(k-1), params, dt);
F = computeJacobianF(x_est(:,k-1), current(k-1), params, dt); % 状态转移雅可比
P_pred = F * P_est(:,:,k-1) * F' + Q;
% 更新步
z_pred = measurementFcn(x_pred, current(k), params);
H = computeJacobianH(x_pred, current(k), params); % 观测雅可比
K = P_pred * H' / (H * P_pred * H' + R);
x_est(:,k) = x_pred + K * (voltage(k) - z_pred);
P_est(:,:,k) = (eye(3) - K*H) * P_pred;
end
雅可比矩阵的计算需要特别注意。以观测雅可比为例:
matlab复制function H = computeJacobianH(x, current, params)
soc = x(1);
% OCV对SOC的导数
docv_dsoc = polyval(polyder(params.ocv_poly), soc);
H = [docv_dsoc, -1, -1]; % 对SOC,V1,V2的偏导
end
3.3 参数调优经验
几个关键参数需要仔细调整:
-
过程噪声协方差Q:反映模型不确定性
- SOC噪声通常设为1e-6
- RC电压噪声设为1e-4
-
观测噪声协方差R:与电压测量精度相关
- 典型值1e-3
-
初始协方差P₀:反映初始状态的不确定性
- SOC不确定性可能设为0.01
- 电压不确定性设为0.1
实际项目中,我通常先用仿真数据调整这些参数,再在真实系统上微调。一个实用技巧是用不同驾驶循环数据测试,观察收敛速度与稳定性。
4. UEKF改进与实现
4.1 UEKF原理优势
EKF的线性化近似在强非线性区域(如SOC接近0%或100%)会引入显著误差。UEKF通过无迹变换(UT)解决了这个问题:
- 选择一组sigma点表征状态分布
- 直接通过非线性函数传播这些点
- 计算传播点的统计特性
这种方法避免了雅可比矩阵计算,且能更准确地捕获非线性特性。
4.2 Sigma点采样策略
常用的对称采样策略:
χ₀ = x̂
χᵢ = x̂ + (√(n+λ)P)ᵢ, i=1,...,n
χᵢ₊ₙ = x̂ - (√(n+λ)P)ᵢ, i=1,...,n
其中λ=α²(n+κ)-n,α和κ是调节参数。
matlab复制function [sigma_points, wm, wc] = generateSigmaPoints(x, P, alpha, beta, kappa)
n = length(x);
lambda = alpha^2 * (n + kappa) - n;
% 计算矩阵平方根
S = chol((n + lambda) * P)';
sigma_points = zeros(n, 2*n+1);
sigma_points(:,1) = x;
for i = 1:n
sigma_points(:,i+1) = x + S(:,i);
sigma_points(:,i+n+1) = x - S(:,i);
end
% 权重计算
wm = zeros(1, 2*n+1);
wc = zeros(1, 2*n+1);
wm(1) = lambda / (n + lambda);
wc(1) = wm(1) + (1 - alpha^2 + beta);
for i = 2:2*n+1
wm(i) = 1 / (2*(n + lambda));
wc(i) = wm(i);
end
end
4.3 UEKF完整实现
matlab复制% UEKF主循环
for k = 2:length(t)
% 生成sigma点
[sigma_points, wm, wc] = generateSigmaPoints(x_est(:,k-1), P_est(:,:,k-1), 1e-3, 2, 0);
% 预测步
sigma_points_pred = zeros(size(sigma_points));
for i = 1:size(sigma_points,2)
sigma_points_pred(:,i) = stateFcn(sigma_points(:,i), current(k-1), params, dt);
end
x_pred = sigma_points_pred * wm';
P_pred = zeros(size(P_est(:,:,k-1)));
for i = 1:size(sigma_points_pred,2)
P_pred = P_pred + wc(i) * (sigma_points_pred(:,i) - x_pred) * (sigma_points_pred(:,i) - x_pred)';
end
P_pred = P_pred + Q;
% 更新步
z_sigma = zeros(1, size(sigma_points_pred,2));
for i = 1:size(sigma_points_pred,2)
z_sigma(i) = measurementFcn(sigma_points_pred(:,i), current(k), params);
end
z_pred = z_sigma * wm';
Pzz = 0;
Pxz = zeros(size(x_pred,1), 1);
for i = 1:size(sigma_points_pred,2)
Pzz = Pzz + wc(i) * (z_sigma(i) - z_pred) * (z_sigma(i) - z_pred)';
Pxz = Pxz + wc(i) * (sigma_points_pred(:,i) - x_pred) * (z_sigma(i) - z_pred)';
end
Pzz = Pzz + R;
K = Pxz / Pzz;
x_est(:,k) = x_pred + K * (voltage(k) - z_pred);
P_est(:,:,k) = P_pred - K * Pzz * K';
end
4.4 UEKF性能对比
在实际测试中,UEKF相比EKF有几个明显优势:
- SOC在极端区域(如<10%或>90%)的估算精度提高约30%
- 对模型参数误差的鲁棒性更好
- 不需要计算复杂的雅可比矩阵
但代价是计算量增加约2-3倍。在资源有限的嵌入式系统中,需要权衡精度与计算开销。
5. 实际应用中的关键问题与解决方案
5.1 初始SOC不确定性问题
在实际应用中,初始SOC往往未知。我常用的解决方案:
- 开机时根据OCV估算初始SOC
- 需要电池静置足够时间(至少2小时)
- OCV-SOC曲线必须准确
matlab复制function soc_init = estimateInitialSOC(voltage, params)
% 查找OCV-SOC曲线上最接近的SOC值
soc_points = 0:0.01:1;
ocv_points = polyval(params.ocv_poly, soc_points);
[~, idx] = min(abs(ocv_points - voltage));
soc_init = soc_points(idx);
end
- 结合安时积分和EKF的混合方法
- 前几分钟主要依赖安时积分
- 当电压变化足够大时启用EKF
5.2 模型参数在线辨识
电池参数会随老化、温度变化而改变。实用的在线参数辨识方法:
- 递推最小二乘法(RLS)
- 带遗忘因子的参数更新
matlab复制% 简单的RLS参数辨识示例
function params = onlineParameterEstimation(params, voltage, current, soc, dt)
persistent P theta;
if isempty(P)
P = eye(3) * 1e6;
theta = [params.R0; params.R1; params.R2];
end
% 构建回归向量
ocv = polyval(params.ocv_poly, soc);
phi = [-current;
-current*(1-exp(-dt/(params.R1*params.C1)));
-current*(1-exp(-dt/(params.R2*params.C2)))];
% RLS更新
K = P * phi / (1 + phi' * P * phi);
error = (ocv - voltage) - phi' * theta;
theta = theta + K * error;
P = (eye(3) - K * phi') * P;
% 更新参数
params.R0 = theta(1);
params.R1 = theta(2);
params.R2 = theta(3);
end
5.3 温度补偿策略
温度对电池行为影响显著。必须实施温度补偿:
- OCV-SOC曲线随温度变化
- 内阻与温度呈指数关系
- 容量随温度降低而减小
matlab复制% 温度补偿示例
function params = temperatureCompensation(params, temp)
% 参考温度25°C
T_ref = 25;
% 内阻温度系数
params.R0 = params.R0_ref * exp(0.01*(T_ref - temp));
params.R1 = params.R1_ref * exp(0.008*(T_ref - temp));
params.R2 = params.R2_ref * exp(0.008*(T_ref - temp));
% 容量调整
params.Qn = params.Qn_ref * (1 - 0.005*(T_ref - temp));
end
5.4 代码优化与实时性考虑
在嵌入式平台实现时,需要优化计算:
- 预先计算常数项
- 采用定点数运算
- 简化矩阵操作
matlab复制% 优化的EKF预测步示例
function [x_pred, P_pred] = optimizedPredict(x_est, P_est, current, params)
persistent F Q;
if isempty(F)
% 预先计算常数部分
F = eye(3);
F(2,2) = exp(-dt/(params.R1*params.C1));
F(3,3) = exp(-dt/(params.R2*params.C2));
Q = diag([1e-6, 1e-4, 1e-4]);
end
% 状态预测
x_pred = F * x_est;
x_pred(1) = x_pred(1) - (current * dt) / (params.Qn * 3600);
% 协方差预测
P_pred = F * P_est * F' + Q;
end
6. 完整Matlab实现与测试案例
6.1 仿真数据生成
为了验证算法,首先需要生成仿真数据:
matlab复制% 生成动态应力测试(DST)电流曲线
function [current, voltage, soc_true] = generateSimulationData(params, duration)
t = 0:1:duration; % 1秒间隔
current = zeros(size(t));
% DST工况
for k = 1:length(t)
phase = mod(t(k), 60);
if phase < 10
current(k) = 1 * params.Qn;
elseif phase < 20
current(k) = -0.5 * params.Qn;
elseif phase < 30
current(k) = 0.75 * params.Qn;
elseif phase < 40
current(k) = -0.25 * params.Qn;
elseif phase < 50
current(k) = 0.5 * params.Qn;
else
current(k) = -1 * params.Qn;
end
end
% 生成真实SOC和电压
soc_true = zeros(size(t));
soc_true(1) = 0.5; % 初始SOC 50%
v1 = 0;
v2 = 0;
voltage = zeros(size(t));
for k = 2:length(t)
dt = t(k) - t(k-1);
soc_true(k) = soc_true(k-1) - (current(k-1) * dt) / (params.Qn * 3600);
v1 = exp(-dt/(params.R1*params.C1)) * v1 + ...
params.R1*(1-exp(-dt/(params.R1*params.C1))) * current(k-1);
v2 = exp(-dt/(params.R2*params.C2)) * v2 + ...
params.R2*(1-exp(-dt/(params.R2*params.C2))) * current(k-1);
ocv = polyval(params.ocv_poly, soc_true(k));
voltage(k) = ocv - current(k)*params.R0 - v1 - v2;
end
% 添加噪声
voltage = voltage + 0.01 * randn(size(voltage));
end
6.2 完整SOC估算流程
matlab复制% 主测试脚本
clear all; close all;
% 电池参数初始化
params.Qn = 2.3 * 3600; % Ah转As
params.R0 = 0.01;
params.R1 = 0.005;
params.C1 = 2000;
params.R2 = 0.01;
params.C2 = 5000;
soc_points = 0:0.1:1;
ocv_points = [3.0 3.3 3.45 3.5 3.55 3.6 3.65 3.7 3.8 3.9 4.1];
params.ocv_poly = polyfit(soc_points, ocv_points, 5);
% 生成仿真数据
[current, voltage, soc_true] = generateSimulationData(params, 600); % 10分钟数据
% EKF初始化
x_est = zeros(3, length(current));
x_est(:,1) = [0.6; 0; 0]; % 初始SOC设为60%(有误差)
P_est = zeros(3, 3, length(current));
P_est(:,:,1) = diag([0.1, 0.1, 0.1]); % 初始协方差
Q = diag([1e-6, 1e-4, 1e-4]); % 过程噪声
R = 1e-3; % 观测噪声
% UEKF初始化
x_est_ukf = x_est;
P_est_ukf = P_est;
% 主循环
for k = 2:length(current)
dt = 1; % 1秒间隔
% EKF估算
[x_est(:,k), P_est(:,:,k)] = ekfStep(x_est(:,k-1), P_est(:,:,k-1), ...
current(k-1), voltage(k), params, dt);
% UEKF估算
[x_est_ukf(:,k), P_est_ukf(:,:,k)] = ukfStep(x_est_ukf(:,k-1), P_est_ukf(:,:,k-1), ...
current(k-1), voltage(k), params, dt);
end
% 结果可视化
figure;
subplot(2,1,1);
plot(1:length(current), soc_true, 'k', ...
1:length(current), x_est(1,:), 'b', ...
1:length(current), x_est_ukf(1,:), 'r');
legend('真实值', 'EKF估算', 'UEKF估算');
xlabel('时间(s)'); ylabel('SOC');
title('SOC估算结果对比');
subplot(2,1,2);
plot(1:length(current), abs(soc_true - x_est(1,:)), 'b', ...
1:length(current), abs(soc_true - x_est_ukf(1,:)), 'r');
legend('EKF误差', 'UEKF误差');
xlabel('时间(s)'); ylabel('SOC误差');
6.3 性能评估指标
为了量化算法性能,我通常计算以下指标:
- 最大绝对误差(MAE)
- 均方根误差(RMSE)
- 收敛时间(从初始误差到<2%的时间)
matlab复制% 计算评估指标
mae_ekf = max(abs(soc_true - x_est(1,:)));
mae_ukf = max(abs(soc_true - x_est_ukf(1,:)));
rmse_ekf = sqrt(mean((soc_true - x_est(1,:)).^2));
rmse_ukf = sqrt(mean((soc_true - x_est_ukf(1,:)).^2));
% 收敛时间计算
threshold = 0.02;
ekf_converge_idx = find(abs(soc_true - x_est(1,:)) < threshold, 1);
ukf_converge_idx = find(abs(soc_true - x_est_ukf(1,:)) < threshold, 1);
fprintf('EKF性能: MAE=%.4f, RMSE=%.4f, 收敛时间=%ds\n', mae_ekf, rmse_ekf, ekf_converge_idx);
fprintf('UEKF性能: MAE=%.4f, RMSE=%.4f, 收敛时间=%ds\n', mae_ukf, rmse_ukf, ukf_converge_idx);
在典型测试中,UEKF相比EKF能有20-30%的精度提升,特别是在SOC极端区域。但计算时间会增加约2倍,需要根据应用场景权衡。
