使用粒子群优化算法（PSO）辨识锂电池二阶RC模型参数（附MATLAB代码）

使用粒子群优化算法来辨识锂离子电池二阶RC模型的参数。将粒子群优化算法寻找到的最优参数代入二阶RC模型进行仿真，经过验证，端电压的估计误差小于0.1%，说明粒子群优化算法辨识得到的参数具有较高的精度，为锂离子电池SOC的估计做铺垫。

新能源佬大

4487人浏览 · 2023-06-07 16:24:55

新能源佬大 · 2023-06-07 16:24:55 发布

文章最下方的长方形小方框内有博主的QQ名片获取本文同款程序

或者关注博主成为粉丝后聊天框会自动发送小企鹅哦！！！！或可直接＋企鹅

硕博期间所有的程序代码，一共2个多g，可以给你指导，赠送半个小时的语音电话答疑。电池数据+辨识程序+各种卡尔曼滤波算法都在里面了，后续还会有新模型的更新。快速入门BMS软件。

一、原理部分

二、代码详解部分

三、结果及分析

一、原理部分

PSO算法由美国学者于 1995 年提出，因其算法简单、效果良好，而在很多领域得到了广泛应用。该算法的起源是模拟鸟群的觅食过程，形成一种群体智能搜索算法。

其核心是，通过迭代的方式来寻求最优解的过程。在算法运行过程中，在每一次迭代中会筛选出微粒最优解和群体最优解。然后，通过向最优解学习的方式，每个粒子都能不断地更新自身的速度和位置，从而逐渐使得整个群体中的所有微粒，都趋近目标函数的最优解，最终达到优化要求。

“粒子”在不同的情景下，代表不同的群体，但是其没有质量和体积，算法中重点强调的是其代表的位置和速度。空间中的所有粒子,都有相对应的表达式，迭代到第t代的第i个粒子，它在空间中的位置和速度如下式所示:

个体极值用来描述， t代的第个微粒运行到当前时刻，所经历的历史最优位置，用下式表示：

全局极值用来描述，整个粒子群算法运行到当下时刻，所搜索到的最优位置，用下式表示：

找出全局极值后，需要对下一代粒子的速度和位置进行更新，以往最优目标方向引导整个粒子群。第i个粒子的速度及位置更新的公式如下：

下标表示搜索空间的维度。c 和c为粒子的学习因子，一般按经验值来取，为非负常数。 c和c分别代表粒子向个体极值和全局值靠拢的一次步长;w为惯性权重其可以代表第t代粒子的速度对当前第 +1 代粒子速度的影响，如果w 偏大，则粒子全局搜索能力强、局部搜索能力太弱:如果w偏小，则粒子全局能力弱，局部搜索能力太强。改变 w 值可以相应的调整粒子的全局搜索能力和局部搜索能力;和为随机数均匀分布在[0，1]区间上。

本文运用粒子群算法进行参数辨识的流程为:设定目标函数，目标函数是关于4个电池参数的函数;对目标函数寻优，有同学可能会好奇为什么是4个参数而不是5个参数，这里解释一下，R0参数的辨识是通过充放电瞬间的压降来识别的，即

，因此本套代码只适用于HPPC工况下的参数辨识，当目标函数收敛至最优值时，相应的参数也收敛至稳定值。

在本项目中，目标函数设置为:

U(k+1)为 k+1时刻，电池的实际输出电压。Y(k+1)为经过粒子群算法输出的跟随电压，它是关于电池5个参数的函数，如下式所示：

当目标函数取得稳定收敛值时，参数也取得了相应合理值。本质为通过目标函数寻优，进行参数辨识。

二、代码详解部分

1.清空环境变量并读取HPPC的电流、电压，时间单位为1s。

clc;clear;close all

HPPC_25_05C = xlsread('013005-2021-07-01 07-30-1209fc.xls','记录列表1');
ZI_SOC_05 = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1];
ZI_OCV_05 = [3.24899 3.36034 3.48125 3.56082 3.60988 3.66923 3.77507 3.86842 3.94964 4.05713 4.16831];

2.拟合OCV-SOC曲线，利用8次多项式拟合

p = polyfit(ZI_SOC_05,ZI_OCV_05,8);   %  多项式拟合参数调用的格式，8次多项式拟合公式

3.数据的处理，如将数据中为NAN的数据删除

V_05C = HPPC_25_05C(:,4)/1000;
I_05C = HPPC_25_05C(:,5);
V_05C(isnan(V_05C)==1) = [];   %  将数据中为NAN的数据删除
I_05C(isnan(I_05C)==1) = [];

4.辨识R0的参数，利用两个电压值来辨识

U = V_05C(11583:end,:)';
I = I_05C(11583:end,:)';

A=[];B=[];C=[];D=[];
for j=1:length(U)-1
    if I(j)==0&&I(j+1)<0
        A=[A j];%即将有电流的点
        B=[B j+1];%刚刚有电流的点
    end
end
for j=1:length(U)-1
    if I(j)<0&&I(j+1)==0
        C=[C j];%电流即将变为0的点
        D=[D j+1];%电流刚变为0的点
    end
end
M1=U(1,A);%即将有电流的点
M2=U(1,B);%刚刚有电流的点
N1=U(1,C);%电流即将变为0的点
N2=U(1,D);%电流刚变为0的点
%   R=(M1-M2+N2-N1)/2/1.375;%计算R0
R=(N2-N1)/1.375;%计算R0
R0_Pingjun=mean(R(2:8)); % 取SOC>20%的R0，平均值认为是欧姆内阻的值。

5.各种参数设置

timer1=cputime;
Ts=1;%系统采样时间
Cn=2.75;
soc_chushizhi=0.5; %SOC的初始值（必须知道脉冲放电或者充电前的SOC，而且是准确的，可根据开路电压曲线知道）   
model_RC=2;     %  模型的选择 
if model_RC==1  %R  int模型
    %  参数上下值界限     Canshu = [R1,C1];
    lb = [1e-03,1e-02];
    ub = [5e-02,1e+06];     
    %   ub = [2e-01,1e+06];     %  0的时候的参数值
    dim = 2;
elseif model_RC==2
    lb = [1e-03,1e-03,1e+01,1e+01];
    ub = [5e-02,3e-02,1e+06,1e+06];
    %  ub = [2e-01,2e-01,1e+06,1e+06];    %  0的时候的参数值
    dim = 4;
elseif model_RC==3
    lb = [1e-03,1e-03,1e-03,1e-02,1e-02,1e-02];
    ub = [5e-02,3e-02,3e-02,1e+06,1e+06,1e+06];
    %  ub = [2e-01,2e-01,2e-01,1e+06,1e+06,1e+06];    %  0的时候的参数值
    dim = 6;
end

6.导入数据，计算ocv、soc真实值、开路电压多项拟合参数选择；实验数据是两列，第一列是电压、第二列是电流（mA）其实，只需要其中辨识脉冲那部分的电池端电压、电流数据即可。

匹配SOC初始值。

if soc_chushizhi==1
    R0 = R(1);
    U = V_05C_1_L;
    I = I_05C_1_L;    
elseif soc_chushizhi==0.9
    R0 = R(2);
    U = V_05C_09_L;
    I = I_05C_09_L;
elseif soc_chushizhi==0.8
    R0 = R(3);
    U = V_05C_08_L;
    I = I_05C_08_L;
elseif soc_chushizhi==0.7
    R0 = R(4);
    U = V_05C_07_L;
    I = I_05C_07_L;
elseif soc_chushizhi==0.6
    R0 = R(5);
    U = V_05C_06_L;
    I = I_05C_06_L;
elseif soc_chushizhi==0.5
    R0 = R(6);
    U = V_05C_05_L;
    I = I_05C_05_L;
elseif soc_chushizhi==0.4
    R0 = R(7);
    U = V_05C_04_L;
    I = I_05C_04_L;
elseif soc_chushizhi==0.3
    R0 = R(8);
    U = V_05C_03_L;
    I = I_05C_03_L;
elseif soc_chushizhi==0.2
    R0 = R(9);
    U = V_05C_02_L;
    I = I_05C_02_L;
elseif soc_chushizhi==0.1
    R0 = R(10);
    U = V_05C_01_L;
    I = I_05C_01_L;
end

7.安时积分法计算soc

N = length(I);
soc(1)=soc_chushizhi; 
ocv(1)=U(1);   %   你自己初始ocv-soc数据 
for i=2:N
    soc(i)=soc(i-1)+I(i)*1/(Cn*3600);   %  额定容量为2Ah
    ocv(i)=polyval(p,soc(i)); 
end

8.PSO部分参数的初始化。

test_number=2;%测试次数，测试多次，去最好那个
SearchAgents_no=100;%种群，最好选择100
Max_iter=500;%迭代次数
Alpha_score=zeros(test_number,1);
Alpha_pos=zeros(test_number,dim);
Convergence_curve=zeros(test_number,Max_iter);
R0_R1_C1_R2_C2_curve=zeros(dim,Max_iter/SearchAgents_no,test_number);%三维矩阵

9.非常重要的function，输入为：种群数量、迭代次数、电流初始值、电压初始值、辨识参数的维度（4维，即R1、R2、C1、C2）、电压序列、电流序列、OCV值、R0、模型阶数；输出为最优粒子群适应度函数值、适应度迭代曲线、端电压误差、端电压值。

for i=1:test_number
        [Alpha_score,Alpha_pos,Convergence_curve,X_error,U_duan]=PSO_R0(SearchAgents_no,Max_iter,lb,ub,dim,U,I,ocv,R0,model_RC);
end