SOP1st_order_sop_电池SOP_一阶SOP_源码

%---SOP方法对比 1.在线辨识部分 2.估计部分 clc clear % load @25-FUDS-1C.mat % load model1.mat % load T25.mat % global P Q R Ts s=[95,90,85,80,75,70,65,60,55,50,45,40,35,30,25,20,15,10,5]; Ro=[0.000725714597082389,6.9808e-04,6.8306e-04,6.7193e-04,6.6453e-04,6.5868e-04,6.5454e-04,6.5223e-04,6.5030e-04,6.4885e-04,6.4750e-04,6.4645e-04,6.4559e-04,6.4452e-04,6.4363e-04,6.4332e-04,6.4460e-04,6.4714e-04,6.5438e-04]; R1=[7.1758e-04,4.7916e-04,4.5125e-04,4.5191e-04,4.6056e-04,4.6156e-04,4.5810e-04,4.4810e-04,4.3796e-04,4.3107e-04,4.3786e-04,4.5775e-04,4.8103e-04,4.9498e-04,5.0863e-04,5.9945e-04,0.0015,0.0032,0.0059]; C1=[5.6508e+04,4.6925e+04,4.1741e+04,3.8877e+04,3.7354e+04,3.7045e+04,3.7160e+04,3.7190e+04,3.7327e+04,3.8149e+04,4.0256e+04,4.3169e+04,4.6186e+04,4.8294e+04,5.0398e+04,5.5481e+04,6.5701e+04,6.7782e+04,6.5320e+04]; % s=[5,10,20,30,40,50,60,70,80,90,100]; % Ro=[0.001049,0.001049,0.0009991,0.0009905,0.000988,0.0009921,0.000999,0.001007,0.001028,0.001082,0.001386]; % R1=[0.01315,0.01315,0.007177,0.004355,0.003025,0.002305,0.001567,0.001572,0.001372,0.00126,0.001641]; % C1=[9.827e4,9.827e4,1.201e5,1.234e5,1.088e5,9.424e4,7.156e4,5.794e4,4.717e4,3.857e4,1.2e4]; c=[1.49887792449226e-15,-6.94926134223021e-13,1.36202271065627e-10,-1.46927185588170e-08,9.52265175450693e-07,-3.81018144485913e-05,0.000939879077959291,-0.0142045617688577,0.141404575917457,2.59254268623833]; % ss=[95;90;85;80;75;70;65;60;55;50;45;40;35;30;25;20;15;10;5]; load('@25-FUDS-1C.mat'); xdata=a(2,:); % ydata=xlsread('清洗后的电压电流数据','B2:B83938'); %votage % soc1=[100,90,80,70,60]; % soc2=[60,50,40,30]; % soc3=[30,20,10]; % ocv1=[4.182357788 4.080653667 3.99410057067871 3.90577459335327 3.79891276359558 ]; % ocv2=[3.79891276359558 3.70397853851318 3.6343491077423 3.54199361801147]; % ocv3=[3.54199361801147 3.4551179409027 3.22253680229187]; % p1=polyfit(soc1,ocv1,2); % p2=polyfit(soc2,ocv2,2); % p3=polyfit(soc3,ocv3,2); % t=a(1,:);xdata=a(2,:);ydata=a(3,:); len=length(xdata);Ts=1; z(1)=1;Cn=112; for k=2:len z(k)=z(k-1)+xdata(k-1)*Ts/Cn/3600; end z=z*100; % P=diag([1 0.2 0.0001 0.0001 100]); % R=1e0; % Q=diag([100 100 0.0001 0.0001 1e1]); % 离线数据进行初始模型---误差太大 Imin(1)=-3*Cn;Imax(1)=3*Cn;Utmin=2.5;Utmax=4.2;smin=5;smax=95;L=1;yita=1; Up(1)=0; % cycle=[]; for k=1:len % [uu,y,Ip1,OCV1,Ro1,Rp1,Cp1,P1]=fcn1(ydata(k),xdata(k),xdata(k-1),Ip,Uoc(k),R0(k),Rp(k),Cp(k)); R0(k)=interp1(s,Ro,z(k)); Rp(k)=interp1(s,R1,z(k)); Cp(k)=interp1(s,C1,z(k)); % if z(k)>60 % ptemp=p1; % elseif z(k)>30 % ptemp=p2; % else % ptemp=p3; % end % Uoc(k)=polyval(ptemp,z(k)); f=@(c,x)c(1)*x.^9+c(2)*x.^8+c(3)*x.^7+c(4)*x.^6+c(5)*x.^5+c(6)*x.^4++c(7)*x.^3+c(8)*x.^2++c(9)*x+c(10); Uoc(k)=f(c,z(k)); if k~=1%不等于 Up(k)=xdata(k-1)*Rp(k)*(1-exp(-Ts/(Rp(k)*Cp(k)))); Imax(k)=Imax(k-1); Imin(k)=Imin(k-1); end % Ut(k)=Uoc(k)-xdata(k)*R0(k)-Up(k); % error(k)=(Ut(k)-ydata(k))*1000;%离线数据误差大 % Current-based % Voltage-based sum1=0; for j=0:1:L-1 sum1=sum1+(exp(-Ts/(Rp(k)*Cp(k))))^(L-1-j); end % if z(k)>60 % dptemp=[2*p1(1),p1(2)]; % elseif z(k)>30 % dptemp=[2*p2(1),p2(2)]; % else % dptemp=[2*p3(1),p3(2)]; % end % dUoc(k)=polyval(dptemp,z(k)); f1=@(c,x)9*c(1)*x.^8+8*c(2)*x.^7+7*c(3)*x.^6+6*c(4)*x.^5+5*c(5)*x.^4+4*c(6)*x.^3+3*c(7)*x.^2+2*c(8)*x+c(9); dUoc(k)=f1(c,z(k)); Idu(k)=(Uoc(k)-Up(k)*exp(-Ts/(Rp(k)*Cp(k)))^L-Utmin)/(L*yita*Ts/(3600*Cn)*dUoc(k)+R0(k)+Rp(k)*(1-exp(-Ts/(Rp(k)*Cp(k))))*sum1); Icu(k)=(Uoc(k)-Up(k)*exp(-Ts/(Rp(k)*Cp(k)))^L-Utmax)/(L*yita*Ts/(3600*Cn)*dUoc(k)+R0(k)+Rp(k)*(1-exp(-Ts/(Rp(k)*Cp(k))))*sum1); %(5.19) % SOC-based if z(k)<=smin Ids(k)=0; else Ids(k)=(z(k)-smin)/(yita*L*Ts/3600/Cn); end if z(k)>smax Ics(k)=0; else Ics(k)=(z(k)-smax)/(yita*L*Ts/3600/Cn);%(5.20) end Idmax(k)=min([Ids(k),Idu(k),Imax(k)]); Icmax(k)=max([Ics(k),Icu(k),Imin(k)]);%(5.21) Udt(k)=Uoc(k)-Up(k)*exp(-Ts/(Rp(k)*Cp(k)))^L-Idmax(k)*(L*yita*Ts/Cn*dUoc(k)+R0(k)+Rp(k)*(1-exp(-Ts/(Rp(k)*Cp(k))))*sum1); Uct(k)=Uoc(k)-Up(k)*exp(-Ts/(Rp(k)*Cp(k)))^L-Idmax(k)*(L*yita*Ts/Cn*dUoc(k)+R0(k)+Rp(k)*(1-exp(-Ts/(Rp(k)*Cp(k))))*sum1);%(5.23) Pd(k)=Idmax(k)*Udt(k); Pc(k)=Icmax(k)*Uct(k); % for j=1:L % ss(j)=z(k)+yita*Ts*Idmax(k)*L/3600/Cn; % R01(j)=interp1(s,Ro,ss(j)); % Rp1(j)=interp1(s,R1,ss(j)); % Cp1(j)=interp1(s,C1,ss(j)); % f=@(c,x)c(1)*x.^9+c(2)*x.^8+c(3)*x.^7+c(4)*x.^6+c(5)*x.^5+c(6)*x.^4++c(7)*x.^3+c(8)*x.^2++c(9)*x+c(10); % Uoc1(j)=f(c,ss(j)); % if j==1 % Up1(j)=Up(k); % else % Up1(j)=Idmax(k)*Rp1(j)*(1-exp(-Ts/(Rp1(j)*Cp1(j)))); % end % Ut1(j)=Uoc1(j)-Idmax(k)*R01(j)-Up1(j); % if Ut1(j)<=2.5 % cycle=[cycle,k]; % % break % end % end %最优化方法 end % figure(1) % plot(z,Icu,'DisplayName','Voltage-based','color','blue','linewidth',2); % hold on; % plot(z,Ics,'DisplayName','SOC-based','color','g','linewidth',2); % plot(z,Icmax,'DisplayName','Current-based','color','r','linewidth',2,'linestyle','--');hold off; % axis([0,100,-500,0]) % legend % % figure(2) % plot(z,Idu,'DisplayName','Voltage-based','color','blue','linewidth',2); % hold on; % plot(z,Ids,'DisplayName','SOC-based','color','g','linewidth',2); % plot(z,Idmax,'DisplayName','Current-based','color','r','linewidth',2,'linestyle','--');hold off; % axis([0,100,0,1100]) % legend figure plot(z,Icu,'DisplayName','Model-based','color','blue','linewidth',2); hold on; plot(z,Ics,'DisplayName','SOC-based','color','g','linewidth',2); plot(z,Imin,'DisplayName','3C limitation','color','c','linewidth',2); plot(z,Icmax,'DisplayName','Multiple constraints','color','r','linewidth',2,'linestyle','--');hold off; % axis([0,len ,-500,0]) xlabel('SOC'); ylabel('Peak charge current/A'); legend figure plot(z,Idu,'DisplayName','Model-based','color','blue','linewidth',2); hold on; plot(z,Ids,'DisplayName','SOC-based','color','g','linewidth',2); plot(z,Imax,'DisplayName','3C limitation','color','c','linewidth',2); plot(z,Idmax,'DisplayName','Multiple constraints','color','r','linewidth',2,'linestyle','--');hold off; % axis([0,len,0,1100]) xlabel('SOC'); ylabel('Peak discharge current/A'); legend figure plot(z,Pd,'color','blue','linewidth',2); xlabel('SOC'); ylabel('Peak discharge power/W'); figure plot(z,Pc,'color','blue','linewidth',2); xlabel('SOC'); ylabel('Peak charge power/W'); % % Ts=1;yita=1;Imax=2.5*Ca;L=60; % % 60s的预测电压 % B0 = 3.04086666982140 ;B1 = 3.55348775822911 ;B2 = -1.49223282237204 ;B3 = -80.6122634694122 ;B4 = 451.723565232495 ;B5 = -1186.32880527585 ; % B6 = 1758.40347731021 ;B7 = -1491.43004268353 ;B8 = 668.249290201307 ;B9 = -120.929433849773 ; % dUoc=9*B9*s^8+8*B8*s^7+7*B7*s^6+6*B6*s^5+5*B5*s^4+4*B4*s^3+3*B3*s^2+2*B2*s+B1; % % Up=-Ik*Rp*(exp(-Ts/(Rp*Cp))-1); % % sum1=0; % for j=0:1:L-1 % sum1=sum1+exp(-Ts/(Rp*Cp))^(L-1-j); % end % % Ut=Uoc-Imax*yita*Ts/Ca*dUoc-Up*exp(-Ts/(Rp*Cp))^L-Imax*(R0+Rp*(1-exp(-Ts/(Rp*Cp)))*sum1); % SOPdis=Imax*Ut; % % % function [SOPdis,Imax,Ut] = peakpower1(Uoc,R0,Rp,Cp,s,Ca,Up) % Ts=1;yita=1;L=1;Utmin=2.5;smin=0.05; % % Imax=2.5*Ca; % % 60s的预测电压 % % B0 = 3.04086666982140 ; % B1 = 3.55348775822911 ;B2 = -1.49223282237204 ;B3 = -80.6122634694122 ;B4 = 451.723565232495 ;B5 = -1186.32880527585 ; % B6 = 1758.40347731021 ;B7 = -1491.43004268353 ;B8 = 668.249290201307 ;B9 = -120.929433849773 ; % dUoc=9*