clear;clc;
Cn=2228*3.6;%電池容量，單位mAh
delta_t=1;%采樣時間
ro=0.13/1.76;%電池內阻
rp1=0.0118;cp1=1.567;
rp2=0.0170;cp2=80.65;%RC環節的內阻和電容值，單位歐姆和法擬合數據1
%?rp1=0.000206;cp1=89.958;
%?rp2=0.0161;cp2=85.17;%RC環節的內阻和電容值，單位歐姆和法，擬合數據2
nita=2200/2228;%放電系數
i=zeros（13648）;
for?k=2:3648
????i（k）=2.2;
end
i（1）=0;%放電電流矩陣
fn=@（x）?78.86*x^7-297*x^6+459.5*x^5-381.6*x^4+187.3*x^3-54.08*x^2+...
????9.074*x+10.57;%擬合soc-ocv表達式
syms?t1;
g=fn（t1）;
g=matlabFunction（diff（g））;%fn對soc的偏導

A=[1?0?0;0?rp1*cp1/（delta_t+rp1*cp1）?0;0?0?rp2*cp2/（delta_t+rp2*cp2）];
B=[-nita*delta_t/Cn;delta_t*rp1/（delta_t+rp1）;delta_t*rp2/（delta_t+rp2）];%A，B矩陣
t=0:delta_t:3647;%仿真序列時間
N=length（t）;
size=[3N];%數組大小
Q=4e-9*eye（3）;%狀態方程中的誤差
R=1e-2;?%觀測方程中的誤差

X=zeros（size）;%狀態變量
Soc=zeros（1N）;%soc
Up1=zeros（1N）;
Up2=zeros（1N）;%RC回路的電壓
E=zeros（1N）;%電池電壓
Uoc=zeros（1N）;%電池端電壓
X（:1）=[1;0;0];%狀態變量初值

for?k=2:N
????X（:k）=A*X（:k-1）+B*i（k）+sqrtm（Q）*randn（31）;
end
for?k=1:N
????Soc（k）=X（1k）;
????Up1（k）=X（2k）;
????Up2（k）=X（3k）;
????E（k）=fn（Soc（k））;
????Uoc（k）=E（k）-i（k）*ro-Up1（k）-Up2（k）+R*randn（1）;
end%求取仿真所用的真值

Xekf=zeros（size）;%EKF下的狀態變量
P0=eye（3）;%協方差矩陣初值
Xekf（:1）=[1;0;0];%EKF狀態變量初值，可更改

for?k=2:N
????Xn=A*Xekf（:k-1）+B*i（k）;%計算先驗估計
????P1=A*P0*A‘+Q;%更新協方差矩陣
????dd=Xn（1）-Soc（k）;%計算先驗估計值與真值的差
????Zm=g（Soc（k））*dd+fn（Soc（k））-Xn（2）-Xn（3）-i（k）*ro;%計算觀測值
????H=[g（Xn（1））?0?0];%雅克比矩陣計算
????K=P1*H‘*inv（H*P1*H‘+R）;%計算卡爾曼增益
????Xekf（:k）=Xn+K*（Uoc（k）-Zm）;%更新狀態
????P0=（eye（3）-K*H）*P1;%更新協方差方程
end

Socekf=zeros（1N）;%Ekf下的soc
Uocekf=zeros（1N）;%Ekf下的端電壓
Eekf=zeros（1N）;%EKF下的電池電壓

for?k=1:N
????Socekf（k）=Xekf（1k）;
????Eekf（k）=fn（Socekf（k））;
????Uocekf（k）=fn（Socekf（k））-i（k）*ro-Xekf（2k）-Xekf（3k）;
end%計算三個值

%依次繪制先驗數據，后驗數據和實驗數據做比較
%這里繪制的是soc-ocv曲線
figure;
plot（X（1:）Uoc‘b‘）;
hold?on;
plot（Xekf（1:）Uocekf‘g‘）;
hold?on;
load（‘2.2socv.mat‘）
plot（socsocv‘r‘）;
grid?on;
legend（‘先驗數據‘‘后驗數據‘‘實驗數據‘）;
xlabel（‘soc‘）;ylabel（‘Uocv/V‘）;
title（‘拓展卡爾曼濾波仿真soc-ocv曲線對比‘）;

%依次繪制先驗數據，后驗數據和實驗數據做比較
%這里繪制的是t-ocv曲線
figure;
plot（tUoc‘r‘）;
hold?on;
load（‘2.2vol.mat‘）
plot（timevol‘b‘）;
hold?on;
plot（tUocekf‘g‘）;
grid?on;
xlabel（‘t/s‘）;ylabel（‘Uocv/V‘）;
legend（‘先驗數據‘‘實驗數據‘‘后驗數據‘）;
title（‘拓展卡爾曼濾波仿真t-Uov曲線對比‘）;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件??????16500??2018-12-27?11:28??2.2socv.mat

?????文件???????5739??2018-12-26?23:13??2.2vol.mat

?????文件???????2607??2018-12-27?15:20??EkfRcTwo.m

-----------?---------??----------?-----??----

????????????????24846????????????????????3