锂离子电池等效电路模型在线参数辨识

% 电模型参数在线参数辨识
% 模型采用的是一阶RC等效电路模型
% 参数辨识算法：带遗忘因子的递推最小二乘法
% 电池采用的是8Ah的软包型三元锂离子电池
% 参数辨识工况采用的是25℃下NEDC工况

% 创建时间：2023/10/18

close all
clear
clc

%% 实验数据载入及参数初始化
% 实验数据载入
NEDC = importdata('NEDC_25deg.mat');        % 电池NEDC工况实验数据
I = NEDC(:,1);      % 电池电流数据（这里规定充电为正，放电为负）
V = NEDC(:,2);      % 电池端电压

L_data = length(I);     %实验数据长度
T = 1;      % 采样时间步长

% 模型未知参数初始化
Uoc = zeros(L_data,1);  % 用于存储不同时刻的开路电压
Ro = zeros(L_data,1);   % 用于存储不同时刻的欧姆内阻
Rp = zeros(L_data,1);   % 用于存储不同时刻的极化内阻
Cp = zeros(L_data,1);   % 用于存储不同时刻的极化电容

% 辨识算法初始化
u = 0.98;               % 遗忘因子
Phi = zeros(1,4);       % 数据向量
thita = zeros(4,1);     % 参数向量
P = 1e6*eye(4);         % 协方差矩阵
K = zeros(4,1);         % 增益矩阵

%% 模型参数在线辨识
for k = 2:L_data
    % FFRLS递推过程
    Phi = [1 V(k-1) I(k) I(k-1)];
    K = P*Phi'/(Phi*P*Phi'+u);
    thita = thita+K*(V(k)-Phi*thita);
    P = (eye(4)-K*Phi)*P/u;
    % 参数解析
    Uoc(k,:) = thita(1)/(1-thita(2));   % 开路电压
    Ro(k,:) = (thita(3)-thita(4))/(1+thita(2));     % 欧姆内阻
    Rp(k,:) = (thita(3)+thita(4))/(1-thita(2))-Ro(k,:);  % 极化内阻
    Cp(k,:) = (1+thita(2))/(2-2*thita(2))/Rp(k,:);       % 极化电容
end

% 这里假设初始时刻的未知参数等于后一时刻的
Uoc(1,:) = Uoc(2,:);
Ro(1,:) = Ro(2,:);
Rp(1,:) = Rp(2,:);
Cp(1,:) = Cp(2,:);

%% 模型精度验证
% 计算工况开始多少个时间步的电流为0
count = 1;
for step = 1:L_data
    if I(step+1,:)-I(step,:) == 0
        count = count+1;
    else
        count = count+1;
        break
    end
end

% 电压存储初始化
Vp = zeros(L_data,1);      % 用于存储不同时刻的极化电压
V_model = zeros(L_data,1);      % 用于存储不同时刻的模型端电压
V_model(1:count,:) = V(1:count,:);
V_err = zeros(L_data,1);        % 用于存储不同时刻的模型端电压误差

for k1 = count:L_data
    % 计算极化端电压
    Vp(k1,:) = Vp(k1-1,:)*exp(-T/(Rp(k1,:)*Cp(k1,:)))+I(k1,:)*Rp(k1,:)*(1-exp(-T/(Rp(k1,:)*Cp(k1,:))));
    % 计算模型端电压
    V_model(k1,:) = Uoc(k1,:)+Vp(k1,:)+I(k1,:)*Ro(k1,:);
    % 模型端电压误差
    V_err(k1,:) = (V_model(k1,:)-V(k1,:))*1000;
end

%% 绘图
% 模型端电压和实际端电压
figure
hold on;
box on;
plot(V_model,'LineWidth',2);
plot(V,'LineWidth',2);
legend('模型端电压','实际端电压');
xlabel('时间/s');
ylabel('电压/V');
% 误差分析
figure
hold on;
box on;
plot(V_err,'LineWidth',2);
xlabel('时间/s');
ylabel('电压误差/mV');

更多内容请关注【Future Car】