简化单粒子SPM电化学模型,ESP,SP,包含测试数据,参数辨识代码以及验证的简化电化学模型P2D,锂离子电池,降阶电化学模型附Matlab代码

简化单粒子SPM电化学模型,ESP,SP,包含测试数据,参数辨识代码以及验证的简化电化学模型P2D,锂离子电池,降阶电化学模型附Matlab代码 ✅作者简介热爱科研的Matlab仿真开发者擅长毕业设计辅导、数学建模、数据处理、程序设计科研仿真。往期回顾关注个人主页Matlab科研工作室 关注我领取海量matlab电子书和数学建模资料个人信条做科研博学之、审问之、慎思之、明辨之、笃行之是为博学慎思明辨笃行。 内容介绍一、引言锂离子电池作为现代能源存储的关键设备其性能的精确预测和优化对于众多应用领域至关重要。电化学模型是理解和预测电池行为的有效工具其中简化单粒子模型SPM以其计算效率和相对准确性在电池研究中备受关注。本文聚焦于基于进化算法EA的 SPM 参数化涵盖简化单粒子 SPM 电化学模型包括 ESP、SP 等变体并包含测试数据、参数辨识代码以及经过验证的简化电化学模型 P2D伪二维模型旨在为锂离子电池的研究和应用提供全面且实用的参考。二、简化单粒子 SPM 电化学模型一基本原理单粒子模型基础单粒子模型SPM将锂离子电池的电极简化为单个球形粒子忽略电极内部的空间分布重点关注锂离子在固相中的扩散以及电极 - 电解液界面的电化学反应。这一简化假设使得模型计算复杂度大幅降低同时仍能捕捉电池行为的主要特征。二模型优势与局限性优势计算效率高适用于实时应用和系统级模拟。由于简化了电极结构减少了空间维度的计算能够快速预测电池的端电压、容量等关键性能指标为电池管理系统BMS的实时控制提供支持。局限性忽略了电极内部的详细结构和电解液中离子浓度的空间分布对于一些需要高精度描述电池内部过程的场景如研究电极材料的微观失效机制模型精度可能不足。三、测试数据一数据采集实验设置选择具有代表性的锂离子电池搭建电池测试平台。该平台包括高精度的恒流充放电设备、温度控制箱以及数据采集系统。在不同的温度条件如 25∘C、40∘C、−10∘C和充放电倍率如 0.5C、1C、2C下对电池进行充放电实验。测量参数记录电池的端电压、充放电电流、温度以及时间等参数。同时为了获取更全面的数据以支持模型参数辨识还可以通过电化学阻抗谱EIS测量电池的内阻随频率的变化以及通过循环伏安法CV获取电极反应的电化学动力学信息。二数据预处理⛳️ 运行结果 部分代码%% Balancing and alignment of half-cells to a full-cell% Method is similar to the method published by Honkura and Horiba in https://doi.org/10.1016/j.jpowsour.2014.04.036% Published in November 2020 by Nikolaos Wassiliadis and Matthias Wanzeladdpath(genpath(pwd));%% Selection of different measurement techniques% Mode 1 - Halfcell GITT Fullcell GITT with DVA objective function% Mode 2 - Halfcell GITT fit Fullcell GITT with DVA objective function% Mode 3 - Halfcell GITT Fullcell GITT with OCV objective function% Mode 4 - Halfcell GITT fit Fullcell GITT with OCV objective function deployedmode 4;%% Main routine with balancing and alignment of half cell potential curves to the full cell potential% Load OCV regressions into workspacetransfer load(OCV_anode); % Anode GITT regressionp_a transfer.p;transfer load(OCV_cathode); % Cathode GITT regressionp_c transfer.p;% Load OCV measurementsHalfcell_GITT load(Halfcell_GITT); % Half cell GITT measurementFullcell_GITT load(Fullcell_GITT); % Full cell GITT measurement% Configure balancing and alignmentreso 10000; % Fitting resolutionDVAstart 0.05; % Considered DVA range starting point as SOCDVAend 0.80; % Considered DVA range ending point as SOC% Mode 1 - Halfcell GITT - fullcell GITT with DVA objective functionif mode 1% Prepare full cell data to chosen resolution step sizeQ_fullcell (0:max(Fullcell_GITT.Results_Fullcell_20.q_ch)/reso:max(Fullcell_GITT.Results_Fullcell_20.q_ch))*3600; % Scale full cell chargeOCV_fullcell interp1(Fullcell_GITT.Results_Fullcell_20.OCV_q_ch*3600, Fullcell_GITT.Results_Fullcell_20.OCV_U_ch, Q_fullcell,spline); % Scale full cell OCV(by spline interpolation)q_dot_fullcell Q_fullcell(2:length(Q_fullcell)); % Calculate first derivative of full cell charge for DVAu_dot_fullcell diff(OCV_fullcell)./diff(Q_fullcell); % Calculate first derivative of full cell OCV for DVAu_rel_dot_fullcell diff(OCV_fullcell)./diff(Q_fullcell/max(Q_fullcell)); % Calculate normalized first derivative of full cell OCV for DVAq_rel_dot_fullcell (Q_fullcell(2:length(Q_fullcell))/max(Q_fullcell)); % Calculate normalized first derivative of full cell charge for DVA% Prepare half cell dataOCV_anode Halfcell_GITT.Results_Anode_20.OCV_U_ch;OCV_cathode Halfcell_GITT.Results_Cathode_20.OCV_U_ch;Q_anode (0:max(abs(Halfcell_GITT.Results_Anode_20.OCV_q_ch))/reso:max(abs(Halfcell_GITT.Results_Anode_20.OCV_q_ch)))*3600; % Scale anode chargeQ_cathode (0:max(Halfcell_GITT.Results_Cathode_20.OCV_q_ch)/reso:max(Halfcell_GITT.Results_Cathode_20.OCV_q_ch))*3600; % Scale cathode charge[x, index] unique(abs(Halfcell_GITT.Results_Anode_20.OCV_q_ch*3600));OCV_anode interp1(x, OCV_anode(index), Q_anode,spline);[x, index] unique(Halfcell_GITT.Results_Cathode_20.OCV_q_ch*3600);OCV_cathode interp1(x, OCV_cathode(index), Q_cathode,spline);% Get relevant DVA rangerg find(DVAstart*max(Q_fullcell)Q_fullcell Q_fullcellDVAend*max(Q_fullcell));% Balancing and alignment (objective function)% Optimization routine minimizes error between measured OCV and calculated OCV from both electrode potentials% x(1) and x(2) offset electrode capacities% x(3) and x(4) scale electrode capacitiesfun (x) ...(diff(OCV_fullcell(rg))./diff(Q_fullcell(rg))) ... % Full-cell DVA-(diff(interp1(Q_cathode*x(4)-x(2), OCV_cathode, Q_fullcell(rg),linear) ... % Cathode DVA-interp1(Q_anode*x(3)-x(1), OCV_anode, Q_fullcell(rg),linear))./diff(Q_fullcell(rg))); % Anode DVA% Mode 2 - Halfcell GITT fitted - fullcell GITT with DVA objective functionelseif mode 2% Prepare full cell data to chosen resolution step sizeQ_fullcell (0:max(Fullcell_GITT.Results_Fullcell_20.q_ch)/reso:max(Fullcell_GITT.Results_Fullcell_20.q_ch))*3600; % Scale full cell chargeOCV_fullcell interp1(Fullcell_GITT.Results_Fullcell_20.OCV_q_ch*3600, Fullcell_GITT.Results_Fullcell_20.OCV_U_ch, Q_fullcell,spline); % Scale full cell OCV(by spline interpolation)q_dot_fullcell Q_fullcell(2:length(Q_fullcell)); % Calculate first derivative of full cell charge for DVAu_dot_fullcell diff(OCV_fullcell)./diff(Q_fullcell); % Calculate first derivative of full cell OCV for DVAu_rel_dot_fullcell diff(OCV_fullcell)./diff(Q_fullcell/max(Q_fullcell)); % Calculate normalized first derivative of full cell OCV for DVAq_rel_dot_fullcell (Q_fullcell(2:length(Q_fullcell))/max(Q_fullcell)); % Calculate normalized first derivative of full cell charge for DVA% Prepare half cell data to chosen resolution step sizetheta 0:1/reso:1; % Calculate normalized resolutionOCV_anode p_a(1)p_a(2)*exp(p_a(3)*theta)p_a(4)*tanh((thetap_a(5))/p_a(6))p_a(7)*tanh((thetap_a(8))/p_a(9))p_a(10)*tanh((thetap_a(11))/p_a(12))p_a(13)*tanh((thetap_a(14))/p_a(15)); % Scale anode OCV (hyperbolic fit)OCV_cathode polyval(p_c, theta); % Scale cathode OCV (poly-fit)Q_anode (0:max(abs(Halfcell_GITT.Results_Anode_20.q_ch))/reso:max(abs(Halfcell_GITT.Results_Anode_20.q_ch)))*3600; % Scale anode chargeQ_cathode (0:max(Halfcell_GITT.Results_Cathode_20.q_ch)/reso:max(Halfcell_GITT.Results_Cathode_20.q_ch))*3600; % Scale cathode charge% Get relevant DVA rangerg find(DVAstart*max(Q_fullcell)Q_fullcell Q_fullcellDVAend*max(Q_fullcell));% Balancing and alignment (objective function)% Optimization routine minimizes error between measured OCV and calculated OCV from both electrode potentials% x(1) and x(2) offset electrode capacities% x(3) and x(4) scale electrode capacitiesfun (x) ...(diff(OCV_fullcell(rg))./diff(Q_fullcell(rg))) ... % Full-cell DVA-(diff(interp1(Q_cathode*x(4)-x(2), OCV_cathode, Q_fullcell(rg),linear) ... % Cathode DVA-interp1(Q_anode*x(3)-x(1), OCV_anode, Q_fullcell(rg),linear))./diff(Q_fullcell(rg))); % Anode DVA% Mode 3 - Halfcell GITT - fullcell GITT with DVA objective functionelseif mode 3% Prepare full cell data to chosen resolution step sizeQ_fullcell (0:max(Fullcell_GITT.Results_Fullcell_20.q_ch)/reso:max(Fullcell_GITT.Results_Fullcell_20.q_ch))*3600; % Scale full cell chargeOCV_fullcell interp1(Fullcell_GITT.Results_Fullcell_20.OCV_q_ch*3600, Fullcell_GITT.Results_Fullcell_20.OCV_U_ch, Q_fullcell,spline); % Scale full cell OCV(by spline interpolation)q_dot_fullcell Q_fullcell(2:length(Q_fullcell)); % Calculate first derivative of full cell charge for DVAu_dot_fullcell diff(OCV_fullcell)./diff(Q_fullcell); % Calculate first derivative of full cell OCV for DVAu_rel_dot_fullcell diff(OCV_fullcell)./diff(Q_fullcell/max(Q_fullcell)); % Calculate normalized first derivative of full cell OCV for DVAq_rel_dot_fullcell (Q_fullcell(2:length(Q_fullcell))/max(Q_fullcell)); % Calculate normalized first derivative of full cell charge for DVA% Prepare half cell dataOCV_anode Halfcell_GITT.Results_Anode_20.OCV_U_ch;OCV_cathode Halfcell_GITT.Results_Cathode_20.OCV_U_ch;Q_anode (0:max(abs(Halfcell_GITT.Results_Anode_20.OCV_q_ch))/reso:max(abs(Halfcell_GITT.Results_Anode_20.OCV_q_ch)))*3600; % Scale anode chargeQ_cathode (0:max(Halfcell_GITT.Results_Cathode_20.OCV_q_ch)/reso:max(Halfcell_GITT.Results_Cathode_20.OCV_q_ch))*3600; % Scale cathode charge[x, index] unique(abs(Halfcell_GITT.Results_Anode_20.OCV_q_ch*3600));OCV_anode interp1(x, OCV_anode(index), Q_anode,spline);[x, index] unique(Halfcell_GITT.Results_Cathode_20.OCV_q_ch*3600);OCV_cathode interp1(x, OCV_cathode(index), Q_cathode,spline);%Balancing and alignment (objective function)% Optimization routine minimizes error between measured OCV and calculated OCV from both electrode potentials% x(1) and x(2) offset electrode capacities% x(3) and x(4) scale electrode capacitiesfun (x) ...OCV_fullcell ... % Full-cell OCV-(interp1(Q_cathode*x(4)-x(2), OCV_cathode, Q_fullcell,linear) ... % Cathode OCV-interp1(Q_anode*x(3)-x(1), OCV_anode, Q_fullcell,linear)); % Anode OCV% Mode 4 - Halfcell GITT fitted - fullcell GITT with OCV objective functionelseif mode 4% Prepare full cell data to chosen resolution step sizeQ_fullcell (0:max(Fullcell_GITT.Results_Fullcell_20.q_ch)/reso:max(Fullcell_GITT.Results_Fullcell_20.q_ch))*3600; % Scale full cell chargeOCV_fullcell interp1(Fullcell_GITT.Results_Fullcell_20.OCV_q_ch*3600, Fullcell_GITT.Results_Fullcell_20.OCV_U_ch, Q_fullcell,spline); % Scale full cell OCV(by spline interpolation)q_dot_fullcell Q_fullcell(2:length(Q_fullcell)); % Calculate first derivative of full cell charge for DVAu_dot_fullcell diff(OCV_fullcell)./diff(Q_fullcell); % Calculate first derivative of full cell OCV for DVAu_rel_dot_fullcell diff(OCV_fullcell)./diff(Q_fullcell/max(Q_fullcell)); % Calculate normalized first derivative of full cell OCV for DVAq_rel_dot_fullcell (Q_fullcell(2:length(Q_fullcell))/max(Q_fullcell)); % Calculate normalized first derivative of full cell charge for DVA% Prepare half cell data to chosen resolution step sizetheta 0:1/reso:1; % Calculate normalized resolutionOCV_anode p_a(1)p_a(2)*exp(p_a(3)*theta)p_a(4)*tanh((thetap_a(5))/p_a(6))p_a(7)*tanh((thetap_a(8))/p_a(9))p_a(10)*tanh((thetap_a(11))/p_a(12))p_a(13)*tanh((thetap_a(14))/p_a(15)); % Scale anode OCV (hyperbolic fit)OCV_cathode polyval(p_c, theta); % Scale cathode OCV (poly-fit)Q_anode (0:max(abs(Halfcell_GITT.Results_Anode_20.q_ch))/reso:max(abs(Halfcell_GITT.Results_Anode_20.q_ch)))*3600; % Scale anode chargeQ_cathode (0:max(Halfcell_GITT.Results_Cathode_20.q_ch)/reso:max(Halfcell_GITT.Results_Cathode_20.q_ch))*3600; % Scale cathode charge% Balancing and alignment (objective function)% Optimization routine minimizes error between measured OCV and calculated OCV from both electrode potentials% x(1) and x(2) offset electrode capacities% x(3) and x(4) scale electrode capacitiesfun (x) ...OCV_fullcell ... % Full-cell OCV-(interp1(Q_cathode*x(4)-x(2), OCV_cathode, Q_fullcell,linear) ... % Cathode OCV-interp1(Q_anode*x(3)-x(1), OCV_anode, Q_fullcell,linear)); % Anode OCVend% Set starting parameter for optimizationsigma_neg 2; % Negative shiftsigma_pos 1250; % Positive shifts_neg 530; % Negative scales_pos 505; % Positive scalex0 [sigma_neg, sigma_pos, s_neg, s_pos]; % Assign to optimization% Optimization configurationoptions optimoptions(lsqnonlin, ...Algorithm,levenberg-marquardt, ...Display,iter-detailed, ...FiniteDifferenceType,central, ...MaxFunctionEvaluations,50000, ...FunctionTolerance, 1e-18, ...StepTolerance, 1e-24);% Objective function hand-overx lsqnonlin(fun, x0,[],[] , options);% Fitting parametersResults.sigma_neg x(1);Results.sigma_pos x(2);Results.s_neg x(3);Results.s_pos x(4);% Capacity and OCV of fullcellResults.Q_fullcell Q_fullcell;Results.OCV_fullcell OCV_fullcell;% Capacity and OCV of fitted half cellsResults.Q_anode Q_anode*x(3)-x(1);Results.Q_cathode Q_cathode*x(4)-x(2);Results.OCV_anode OCV_anode;Results.OCV_cathode OCV_cathode;% Used half cell OCV range (q Q_fullcell)Results.OCV_anode_fit interp1(Q_anode*x(3)-x(1), OCV_anode, Q_fullcell,linear);Results.OCV_cathode_fit interp1(Q_cathode*x(4)-x(2), OCV_cathode, Q_fullcell,linear);% Superposition of the used half cell OCV range (q Q_fullcell)Results.OCV_fullcell_fit interp1(Q_cathode*x(4)-x(2), OCV_cathode, Q_fullcell,linear)-interp1(Q_anode*x(3)-x(1), OCV_anode, Q_fullcell,linear);% Relative start- and endpoint of the used capacity ranges of the fitted half cell capacitiesResults.Q_anode_SP (min(Q_fullcell)-min(Q_anode*x(3)-x(1)))/(max(Q_anode*x(3)-x(1))-min(Q_anode*x(3)-x(1)));Results.Q_anode_EP (max(Q_fullcell)-min(Q_anode*x(3)-x(1)))/(max(Q_anode*x(3)-x(1))-min(Q_anode*x(3)-x(1)));Results.Q_cathode_SP (min(Q_fullcell)-min(Q_cathode*x(4)-x(2)))/(max(Q_cathode*x(4)-x(2))-min(Q_cathode*x(4)-x(2)));Results.Q_cathode_EP (max(Q_fullcell)-min(Q_cathode*x(4)-x(2)))/(max(Q_cathode*x(4)-x(2))-min(Q_cathode*x(4)-x(2)));% Goodness of FitResults.NRMSE goodnessOfFit(Results.OCV_fullcell_fit,Results.OCV_fullcell,NRMSE);Results.r2 1 - (norm(Results.OCV_fullcell_fit-Results.OCV_fullcell)/norm(Results.OCV_fullcell-mean(Results.OCV_fullcell))).^2;% Recalculate fullcell OCV from balanced and aligned electrode OCVsu_fullcell_fit interp1(Q_cathode*x(4)-x(2), OCV_cathode, Q_fullcell,linear)-interp1(Q_anode*x(3)-x(1), OCV_anode, Q_fullcell,linear); % Calculate fullcell OCVu_dot_fullcell_fit diff(u_fullcell_fit)./diff(Q_fullcell); % Calculate first derivative of fullcell OCVu_rel_dot_fullcell_fit diff(u_fullcell_fit)./diff(Q_fullcell/max(Q_fullcell)); % Calculate normalized first derivative of fullcell OCVq_rel_dot_fullcell_fit (Q_fullcell(2:length(Q_fullcell))/max(Q_fullcell)); % Calculate state-of-charge of fullcell OCV% Recalculate halfcell OCVs from balanced and aligned electrode OCVsOCV_anode_fit interp1(Q_anode*x(3)-x(1), OCV_anode, Q_fullcell,linear); % Shift and scale anode OCVu_dot_anode_fit diff(OCV_anode_fit)./diff(Q_fullcell); % Calculate first derivative of anode OCVu_rel_dot_anode_fit diff(OCV_anode_fit)./diff(Q_fullcell/max(Q_fullcell)); % Calculate state-of-charge of anode OCVOCV_cathode_fit interp1(Q_cathode*x(4)-x(2), OCV_cathode, Q_fullcell,linear); % Shift and scale cathode OCVu_dot_cathode_fit diff(OCV_cathode_fit)./diff(Q_fullcell); % Calculate first derivative of cathode OCVu_rel_dot_cathode_fit diff(OCV_cathode_fit)./diff(Q_fullcell/max(Q_fullcell)); % Calculate state-of-charge of cathode OCV% Fitting resultsResults.q_dot_fullcell q_dot_fullcell;Results.u_dot_fullcell u_dot_fullcell;Results.u_dot_fullcell_fit u_dot_fullcell_fit;Results.u_dot_anode_fit u_dot_anode_fit;Results.u_dot_cathode_fit u_dot_cathode_fit;Results.q_rel_dot_fullcell q_rel_dot_fullcell;Results.u_rel_dot_fullcell u_rel_dot_fullcell;Results.u_rel_dot_fullcell_fit u_rel_dot_fullcell_fit;Results.u_rel_dot_anode_fit u_rel_dot_anode_fit;Results.u_rel_dot_cathode_fit u_rel_dot_cathode_fit;% Calculate OCV goodness of fit for 0-100% SOC rangeOCV_NRMSE goodnessOfFit(u_fullcell_fit,OCV_fullcell,NRMSE); % Calculate OCV RMSEOCV_r2 1 - (norm(u_fullcell_fit-OCV_fullcell)/norm(OCV_fullcell-mean(OCV_fullcell))).^2; % Calculate OCV R square% Calculate DVA goodness of fit for 10-80% fast charging rangeindbeg find(q_dot_fullcellDVAstart*max(q_dot_fullcell), 1,first); % 10% SOC indexindend find(q_dot_fullcellDVAend*max(q_dot_fullcell), 1,first); % 80% SOC indexDVA_NRMSE goodnessOfFit(u_dot_fullcell_fit(indbeg:indend),u_dot_fullcell(indbeg:indend),NRMSE); % Calculate DVA RMSEDVA_r2 1 - (norm(u_dot_fullcell_fit(indbeg:indend)-u_dot_fullcell(indbeg:indend))/norm(u_dot_fullcell(indbeg:indend)-mean(u_dot_fullcell(indbeg:indend)))).^2; % Calculate DVA R square% Visualize datafigure(units,normalized,outerposition,[0 0 1 1]);set(gcf,color,w);% Plot 1: OCV over charge throughputsubplot(2,1,1);hold on;title(Halfcell GITT aligned and balanced to fullcell GITT - DVA objective function)xlabel(Q in As)ylabel(U in V)plot(Results.Q_fullcell,Results.OCV_fullcell,k);plot(Results.Q_fullcell,Results.OCV_fullcell_fit,g);plot(Results.Q_anode,Results.OCV_anode,k--);plot(Results.Q_cathode,Results.OCV_cathode,k--);plot(Results.Q_fullcell,Results.OCV_anode_fit,b);plot(Results.Q_fullcell,Results.OCV_cathode_fit,r);legend({Fullcell,Fullcell fit,Unused anode capacity,Unused cathode capacity,Anode,Cathode},interpreter,latex,Location,East);xlim([-1500 11500]);NE [max(xlim) max(ylim)]-[diff(xlim) diff(ylim)]*0.05;text(NE(1),NE(2),{[NRMSE num2str(OCV_NRMSE,4)], [R2 num2str(OCV_r2,4)]}, VerticalAlignment,top, HorizontalAlignment,right);% Plot 2: DVA over charge throughputsubplot(2,1,2);plot(q_dot_fullcell, u_dot_fullcell, k, q_dot_fullcell, u_dot_fullcell_fit,g, q_dot_fullcell, u_dot_anode_fit, b, q_dot_fullcell, u_dot_cathode_fit, r);hold on;ylabel(dU/dQ in V/As);xlabel(Q in As);title(DVA);ylim([-2.5e-4 3.5e-4]);xlim([-1500 11500]);legend({Fullcell,Fullcell Fit,Anode,Cathode},interpreter,latex,Location,East);NE [max(xlim) max(ylim)]-[diff(xlim) diff(ylim)]*0.05;text(NE(1),NE(2),{[NRMSE num2str(DVA_NRMSE,4)], [R2 num2str(DVA_r2,4)]}, VerticalAlignment,top, HorizontalAlignment,right);% Save resultsclearvars -except Fullcell Anode Cathode Results DVA_NRMSE DVA_r2 OCV_NRMSE OCV_r2 C100_DVA_r2save B_A Results DVA_NRMSE DVA_r2 OCV_NRMSE OCV_r2 参考文献更多免费数学建模和仿真教程关注领取