Untitled3_SOC的matlab程序_SOC_SOC估计_电池SOC_EKFSOC_

syms I V1 SOC T=100000;%length(simout.signals.values); Q=[10 0;0 0.1]; R=10; P=[100 0;0 1]; w=sqrt(Q)*randn(2,T);%产生过程噪声 v=sqrt(R)*randn(1,T);%产生观测噪声 OCV=zeros(1,T); R0=zeros(1,T); R1=zeros(1,T); C1=zeros(1,T); SOC=zeros(1,T); V1=zeros(1,T); I=zeros(1,T); y=zeros(1,T); for k=1:T; y(k)=simout.signals.values(k,4)+v(k);%加入观测噪声的输出值 end I(1)=simout.signals.values(1,1); SOC(1)=simout.signals.values(1,5); OCV(1)=485.76*SOC(1)^11-2044.70*SOC(1)^10+3000.00*SOC(1)^9-934.20*SOC(1)^8-2309.12*SOC(1)^7+3000.00*SOC(1)^6-1426.11*SOC(1)^5+148.51*SOC(1)^4+125.72*SOC(1)^3-53.31*SOC(1)^2+8.58*SOC(1)+3.04; R0(1)=0.065*SOC(1)^5-0.192*SOC(1)^4+0.215*SOC(1)^3-0.109*SOC(1)^2+0.022*SOC(1)+0.005; R1(1)=-0.396*SOC(1)^5+1.114*SOC(1)^4-1.180*SOC(1)^3+0.576*SOC(1)^2-0.128*SOC(1)+0.015; C1(1)=-40000.00*SOC(1)^5+1317.82*SOC(1)^4+40000.00*SOC(1)^3-40000.006*SOC(1)^2+19422.74*SOC(1)+12199.92; V1(1)=OCV(1)-I(1)*R0(1)-simout.signals.values(1,4); for k=1:T; I(k)=simout.signals.values(k,1); SOC(k)=simout.signals.values(k,5);%荷电状态真实值 OCV(k)=485.76*SOC(k)^11-2044.70*SOC(k)^10+3000.00*SOC(k)^9-934.20*SOC(k)^8-2309.12*SOC(k)^7+3000.00*SOC(k)^6-1426.11*SOC(k)^5+148.51*SOC(k)^4+125.72*SOC(k)^3-53.31*SOC(k)^2+8.58*SOC(k)+3.04; R0(k)=0.065*SOC(k)^5-0.192*SOC(k)^4+0.215*SOC(k)^3-0.109*SOC(k)^2+0.022*SOC(k)+0.005; V1(k)=OCV(k)-I(k)*R0(k)-simout.signals.values(k,4);%V1的真实值 end SOCekf=zeros(1,T); SOCekf(1)=SOC(1); V1ekf=zeros(1,T); V1ekf(1)=V1(1); Yekf=zeros(1,T); Yekf(1)=y(1); R1ekf=zeros(1,T); R1ekf(1)=R1(1); OCVekf=zeros(1,T); OCVekf(1)=OCV(1); R0ekf=zeros(1,T); R0ekf(1)=R0(1); C1ekf=zeros(1,T); C1ekf(1)=C1(1); t=zeros(1,T); e=zeros(1,T); f=zeros(1,T); g=zeros(1,T); n1=zeros(1,T); m1=zeros(1,T); t(1)=R1(1)*C1(1); f3=zeros(1,T); f4=zeros(1,T); h1=zeros(1,T); for k=2:T; t(k)=R1ekf(k-1)*C1ekf(k-1); e(k-1)=[-40000.00 1317.82 40000.00 -40000.006 +19422.74]*[5*SOCekf(k-1)^4;4*SOCekf(k-1)^3;3*SOCekf(k-1)^2;2*SOCekf(k-1);1]; %C1求导 f(k-1)=[-0.396 1.114 -1.180 0.576 -0.128]*[5*SOCekf(k-1)^4;4*SOCekf(k-1)^3;3*SOCekf(k-1)^2;2*SOCekf(k-1);1]; %R1求导 m1(k-1)=f(k-1)*C1(k-1)+R1(k-1)*e(k-1);%时间常数t求导 f3(k-1)=-(V1ekf(k-1)+I(k-1)*R1ekf(k-1))*m1(k-1)/t(k-1)+I(k-1)*f(k-1)*exp(-1/t(k-1)); f4(k-1)=exp(-1/t(k-1)); F=[1 0;f3(k-1) f4(k-1)]; P=F*P*F'+Q; SOCn(k)=SOCekf(k-1)-I(k-1)/(20*3600);%时间更新 V1n(k)=V1ekf(k-1)*exp(-1/t(k-1))+I(k-1)*R1(k-1)*(1-exp(-1/t(k-1))); OCVn(k)=485.76*SOCn(k)^11-2044.70*SOCn(k)^10+3000.00*SOCn(k)^9-934.20*SOCn(k)^8-2309.12*SOCn(k)^7+3000.00*SOCn(k)^6-1426.11*SOCn(k)^5+148.51*SOCn(k)^4+125.72*SOCn(k)^3-53.31*SOCn(k)^2+8.58*SOCn(k)+3.04; R0n(k)=0.065*SOCn(k)^5-0.192*SOCn(k)^4+0.215*SOCn(k)^3-0.109*SOCn(k)^2+0.022*SOCn(k)+0.005; g(k)=[485.76 -2044.70 3000.00 -934.2 -230912 300 -1426011 148.51 125.72 -3.31 8.58]*[11*SOCn(k)^10;10*SOCn(k)^9;9*SOCn(k)^8;8*SOCn(k)^7;7*SOCn(k)^6;6*SOCn(k)^5;5*SOCn(k)^4;4*SOCn(k)^3;3*SOCn(k)^2;2*SOCn(k);1]; %OCV求导 n1(k)=[0.065 -0.192 0.215 -0.109 0.022]*[5*SOCn(k)^4;4*SOCn(k)^3;3*SOCn(k)^2;2*SOCn(k);1]; %R0求导 h1(k)=g(k)-n1(k)*I(k); H=[h1(k) -1]; %求偏导 K=P*H'*inv(H*P*H'+R);%测量更新 SOCekf(k)=SOCn(k)+K(1,1)*(y(k)-OCVn(k)+V1n(k)+I(k)*R0n(k));%测量更新v(k);% V1ekf(k)=V1n(k)+K(2,1)*(y(k)-OCVn(k)+V1n(k)+I(k)*R0n(k)); P=(eye(1)-K*H)*P; OCVekf(k)=485.76*SOCekf(k)^11-2044.70*SOCekf(k)^10+3000.00*SOCekf(k)^9-934.20*SOCekf(k)^8-2309.12*SOCekf(k)^7+3000.00*SOCekf(k)^6-1426.11*SOCekf(k)^5+148.51*SOCekf(k)^4+125.72*SOCekf(k)^3-53.31*SOCekf(k)^2+8.58*SOCekf(k)+3.04; R0ekf(k)=0.065*SOCekf(k)^5-0.192*SOCekf(k)^4+0.215*SOCekf(k)^3-0.109*SOCekf(k)^2+0.022*SOCekf(k)+0.005; R1ekf(k)=-0.396*SOCekf(k)^5+1.114*SOCekf(k)^4-1.180*SOCekf(k)^3+0.576*SOCekf(k)^2-0.128*SOCekf(k)+0.015; C1ekf(k)=-40000.00*SOCekf(k)^5+1317.82*SOCekf(k)^4+40000.00*SOCekf(k)^3-40000.006*SOCekf(k)^2+19422.74*SOCekf(k)+12199.92; %Yekf(k)=OCVekf(k)-V1ekf(k)-I(k)*R0ekf(k); end SOCstd=zeros(1,T); V1std=zeros(1,T); for k=1:T; SOCstd(k)=abs(SOCekf(k)-SOC(k)); V1std(k)=abs(V1ekf(k)-V1(k)); end figure hold on;box on; plot(SOC,'MarkerFace','k'); plot(SOCekf,'MarkerFace','b'); legend('SOC真实值','EKF滤波估计值'); xlabel('时间/s'); ylabel('状态值x'); figure hold on;box on; plot(V1,'MarkerFace','k'); plot(V1ekf,'MarkerFace','b'); legend('V1真实值','EKF滤波估计值'); xlabel('时间/s'); ylabel('状态值x'); figure hold on;box on; plot(SOCstd,'MarkerFace','g'); xlabel('时间/s'); ylabel('SOC状态估计偏差'); figure hold on;box on; plot(V1std,'MarkerFace','g'); xlabel('时间/s'); ylabel('V1状态估计偏差')