使用EKF扩展卡尔曼滤波进行二阶RC模型的参数辨识

原创已于 2025-03-14 17:03:46 修改

· 853 阅读

12 ·

版权

文章标签：

#算法 #matlab #能源

于 2025-03-11 18:02:17 首次发布

锂离子电池二阶RC模型同时被 3 个专栏收录

13 篇文章

订阅专栏

EKF参数辨识

1 篇文章

订阅专栏

扩展卡尔曼滤波辨识电池

1 篇文章

订阅专栏

这次来写一篇使用EKF扩展卡尔曼滤波进行二阶RC模型的参数辨识的程序及讲解，目前使用的比较多的参数辨识方法是离线参数辨识（指数拟合HPPC工况端电压数据）、遗忘因子最小二乘法这两种算法，而EKF用的多的地方是在参数辨识完成后的后续SOC估计中，下面就来介绍一下EKF算法进行参数辨识与最小二乘法的一些核心区别：

EKF是一种针对非线性系统的状态估计算法，通过泰勒展开对非线性模型进行局部线性化，结合卡尔曼滤波框架实现参数和状态的联合估计。在二阶RC模型中，EKF将电路参数（如R0、R1、R2、C1、C2）和状态变量（如SOC、极化电压）共同作为状态向量，利用系统动态方程和观测方程进行递推更新其核心步骤包括状态预测、协方差预测、卡尔曼增益计算和状态更新，涉及雅可比矩阵的推导；

最小二乘法是一种基于误差平方和最小化的优化方法，通过递推公式逐步更新参数估计值。递推最小二乘法（RLS）通过引入遗忘因子（如FRLS）降低旧数据的影响，解决数据饱和问题。在二阶RC模型中，RLS通常仅用于参数辨识，如辨识R0、R1等，而不涉及状态估计；

EKF在线参数辨识的原理及推导过程

EKF的基本思想是将非线性系统线性化后，对线性化后的系统进行卡尔曼滤波，实现状态变量的估计。非线性系统的模型状态空间表达式进行离散化得到状态方程和输出方程：

为适应EKF算法建立在线辨识锂离子电池参数的状态空间方程，引入离散状态方程和输出方程如下：

EKF在线参数辨识方法的具体步骤如表下：

步骤一：初始化参数状态变量和参数误差协方差

步骤二：参数状态估计和参数误差协方差时间更新

步骤三：计算卡尔曼增益矩阵

步骤四：参数状态估计和参数误差协方差测量更新

基于卡尔曼滤波算法的在线辨识方法主要由两个阶段构成，第一阶段为时间更新，第二阶段为测量更新。时间更新过程在卡尔曼滤波中被称为预测过程，因为它对当前状态变量进行推算，以提供下一个时刻的先验估计。而测量更新过程在卡尔曼滤波中被称为校正过程，因为它通过反馈观测值来进行修正。基于 EKF 的模型在线参数辨识方法的工作原理流程图如图下图所示。

EKF参数辨识的实现代码分析

1.数据导入与仿真参数设置

load curvol
para_true.Cn = 2.3*3600; % 额定容量(A·s)
OCV_SOC = OCV_data;
Ts = 1;                  % 采样时间(s)
sim_time = 1179;         % 总仿真时间(s)
t = 0:Ts:sim_time;       % 时间序列
N = length(t);           % 数据点数
I = cur;
V_meas = vol;

2.EKF参数初始化：待辨识参数初始猜测

para_est.R0 = 0.015;     
para_est.R1 = 0.003;    
para_est.R2 = 0.01;     
para_est.C1 = 5000;     
para_est.C2 = 10000;    
para_est.Cn = 2.3*3600;

3.状态向量: [SOC; V1; V2; R0; R1; R2; C1; C2]

x_est = [0.5; 0; 0; para_est.R0; para_est.R1; para_est.R2;...
        para_est.C1; para_est.C2];
n_states = length(x_est);

4.协方差矩阵与噪声矩阵初始化

P = diag([1e-4, 1e-4, 1e-4, 1e-6, 1e-6, 1e-6, 1e-3, 1e-3]);
Q = diag([1e-6, 1e-6, 1e-6, 1e-8, 1e-8, 1e-8, 1e-4, 1e-4]);
R = 1e-4;

5.主循环 - EKF参数辨识

for k = 1:N
    % 获取当前输入
    current = I(k);
    
    % EKF预测步骤
    [x_pred, F] = ekfPredict(x_est, current, Ts);
    
    % 协方差预测
    P_pred = F*P*F' + Q;
    
    % EKF更新步骤
    [V_est, H] = getObsModel(x_pred, current, OCV_SOC);
    y_res = V_meas(k) - V_est;
    
    % 卡尔曼增益
    S = H*P_pred*H' + R;
    K = P_pred*H'/S;
    
    % 状态更新
    x_est = x_pred + K*y_res;
    
    % 协方差更新
    P = (eye(n_states) - K*H)*P_pred;
    
    % 参数限幅（保证物理合理性）
    x_est(4:8) = max(x_est(4:8), [0.001; 0.001; 0.001; 100; 100]);
    
    % 保存参数和电压估计
    para_history(k,:) = x_est(4:8)';
    voltage_est(k) = V_est;
end

6.结果可视化，绘图。

figure('Position', [100,100,1000,800])
% 电压对比
V_meas=V_meas(1:1180);
subplot(3,1,1)
plot(t, V_meas, 'b', t, voltage_est, 'r')
title(sprintf('端电压对比 (RMSE=%.4f V)', sqrt(mean((V_meas-voltage_est).^2))))
ylabel('Voltage(V)'), legend('Measured', 'Estimated'), grid on
% 电阻参数收敛
subplot(3,1,2)
plot(t, para_history(:,1:3), 'LineWidth',1.5)
hold on
plot(t,para_history(:,1), 'k--');
plot(t,para_history(:,2), 'k--');
plot(t,para_history(:,3), 'k--');
% plot([t(1) t(end)], [para_true.R0 para_true.R0], 'k--')
% plot([t(1) t(end)], [para_true.R1 para_true.R1], 'k--')
% plot([t(1) t(end)], [para_true.R2 para_true.R2], 'k--')
title('电阻参数收敛过程'), ylabel('R(Ω)')
legend('R0','R1','R2','True Values'), grid on
% 电容参数收敛
subplot(3,1,3)
plot(t, para_history(:,4:5), 'LineWidth',1.5)
hold on
plot(t,para_history(:,4), 'k--');
plot(t,para_history(:,5), 'k--');
title('电容参数收敛过程'), xlabel('Time(s)'), ylabel('C(F)')
legend('C1','C2','True Values'), grid on

7.function函数——电池模型仿真函数

function [V, SOC] = simBatteryModel(para, I, Ts, OCV_SOC)
    n = length(I);
    SOC = zeros(n,1);
    V1 = zeros(n,1);
    V2 = zeros(n,1);
    
    for k = 2:n
        % SOC更新
        SOC(k) = SOC(k-1) - I(k-1)*Ts/para.Cn;
        SOC(k) = max(min(SOC(k),1),0);
        
        % RC动态
        V1(k) = exp(-Ts/(para.R1*para.C1))*V1(k-1) + ...
                para.R1*(1-exp(-Ts/(para.R1*para.C1)))*I(k-1);
        V2(k) = exp(-Ts/(para.R2*para.C2))*V2(k-1) + ...
                para.R2*(1-exp(-Ts/(para.R2*para.C2)))*I(k-1);
    end
    % 端电压计算
    OCV = interp1(OCV_SOC(:,1), OCV_SOC(:,2), SOC);
    V = OCV - I*para.R0 - V1 - V2;
end

8.EKF预测模型

function [x_pred, F] = ekfPredict(x, I, Ts)
    % 状态分解
    SOC = x(1);
    V1 = x(2);
    V2 = x(3);
    R0 = x(4);
    R1 = x(5);
    R2 = x(6);
    C1 = x(7);
    C2 = x(8);
    
    % SOC预测
    SOC_pred = SOC - I*Ts/C2;
    
    % RC动态预测
    tau1 = R1*C1;
    alpha1 = exp(-Ts/tau1);
    beta1 = R1*(1-alpha1);
    
    tau2 = R2*C2;
    alpha2 = exp(-Ts/tau2);
    beta2 = R2*(1-alpha2);
    
    V1_pred = alpha1*V1 + beta1*I;
    V2_pred = alpha2*V2 + beta2*I;
    
    % 参数动态（假设缓慢变化）
    R0_pred = R0;
    R1_pred = R1;
    R2_pred = R2;
    C1_pred = C1;
    C2_pred = C2;
    
    x_pred = [SOC_pred; V1_pred; V2_pred; R0_pred; R1_pred; R2_pred; C1_pred; C2_pred];
    
    % 雅可比矩阵计算
    F = zeros(8,8);
    
    % SOC相关
    F(1,1) = 1;          % ∂SOC/∂SOC
    F(1,8) = I*Ts/C2^2;  % ∂SOC/∂C2
    
    % V1相关
    F(2,2) = alpha1;     % ∂V1/∂V1
    F(2,5) = (Ts*V1/(R1^2*C1)) * alpha1 + (I*Ts/C1) * alpha1;  % ∂V1/∂R1
    F(2,7) = (Ts*V1/(R1*C1^2)) * alpha1 + (I*Ts/(R1*C1^2)) * (1-alpha1);  % ∂V1/∂C1
    F(2,5) = F(2,5) - beta1*I/R1;  % 修正项
    
    % V2相关
    F(3,3) = alpha2;     % ∂V2/∂V2
    F(3,6) = (Ts*V2/(R2^2*C2)) * alpha2 + (I*Ts/C2) * alpha2;  % ∂V2/∂R2
    F(3,8) = (Ts*V2/(R2*C2^2)) * alpha2 + (I*Ts/(R2*C2^2)) * (1-alpha2);  % ∂V2/∂C2
    F(3,6) = F(3,6) - beta2*I/R2;  % 修正项
    
    % 参数自相关
    F(4:8,4:8) = eye(5);
end

辨识结果及分析

最后来演示下最终的参数辨识结果，由于MATLAB版本的问题，导致显示不出标题，下图中从上到下的标题分别为：仿真的端电压与真实端电压的对比、电阻参数收敛过程、电容参数收敛过程。从图中可以看出，所提出的EKF算法可以迅速收敛到真实值附近，参数辨识的精度也能够满足后续SOC估计的需求。通过锂离子电池动态工况实验中的FUDS动态工况对锂离子电池二阶RC等效电路模型进行EKF在线参数辨识。锂离子电池二阶RC等效电路模型在线辨识的各电阻、电容参数辨识结果如下图所示。通过EKF辨识结果可以看出欧姆内阻和极化内阻在辨识初期参数变化比较剧烈，这是因为模型在初始化阶段选取初值不当。R0在整个工况下辨识基本在40毫欧左右；R1和 R2保持在0.1Ω左右；C1和C2电容基本在5000法拉和10000法拉左右。参数的变化说明锂离子电池是一个时变系统。