卡爾曼濾波法（Kalman濾波）用于SOC估算卡爾曼濾波法（Kalman濾波）1卡爾曼濾波的由來卡爾曼濾波的由來卡爾曼，全名RudolfEmilKalman，匈牙利數(shù)學(xué)家，1930年出生于匈牙利首都布達(dá)佩斯。1953，1954年于麻省理工學(xué)院分別獲得電機(jī)工程學(xué)士及碩士學(xué)位。1957年于哥倫比亞大學(xué)獲得博士學(xué)位。我們在現(xiàn)代控制理論中要學(xué)習(xí)的卡爾曼濾波器，正是源于他的博士論文和1960年發(fā)表的論文《ANewApproachtoLinearFilteringandPredictionProblems》（線性濾波與預(yù)測問題的新方法）?？柭鼮V波的由來卡爾曼濾波的由來卡爾曼，全名RudolfE2Kalman濾波主要應(yīng)用（1）導(dǎo)航制導(dǎo)、目標(biāo)定位和跟蹤領(lǐng)域。（2）通信與信號處理、數(shù)字圖像處理、語音信號處理。（3）天氣預(yù)報(bào)、地震預(yù)報(bào)。（4）地質(zhì)勘探、礦物開采。（5）故障診斷、檢測。（6）證券股票市場預(yù)測。Kalman濾波主要應(yīng)用（1）導(dǎo)航制導(dǎo)、目標(biāo)定位和跟蹤領(lǐng)域。3卡爾曼濾波法：SOC在電流積分法基礎(chǔ)上進(jìn)行測量值（電流積分法）和估算值（初始值）利用前一時(shí)刻的估計(jì)值和現(xiàn)時(shí)刻的觀測值來更新對狀態(tài)變量的估計(jì)，求出現(xiàn)時(shí)刻的估計(jì)值。計(jì)算Covariance（協(xié)方差）Kalman增益：略Matlabmex語言DSPSystemToolbox模塊卡爾曼濾波法：SOC在電流積分法基礎(chǔ)上進(jìn)行測量值（電流積分法4對SOC初始值0.7不敏感對SOC初始值0.7不敏感5SOC初始值：0.95SOC初始值：0.956卡爾曼濾波法優(yōu)點(diǎn)：1.克服電流積分法對初始值依賴的嚴(yán)重缺點(diǎn)2.能夠消除采樣噪聲缺點(diǎn)：模型參數(shù)會(huì)隨時(shí)間變化，需修正計(jì)算量大，一個(gè)采樣周期難以完成計(jì)算計(jì)算機(jī)字長有限造成舍入、截?cái)嗾`差積累卡爾曼濾波法優(yōu)點(diǎn)：2.能夠消除采樣噪聲缺點(diǎn)：模型參數(shù)會(huì)隨7例子：房間溫度某一時(shí)刻的溫度：估算：估計(jì)偏差–高斯白噪聲用溫度計(jì)測量：測量偏差–高斯白噪聲==》用上面2個(gè)值估算下一時(shí)刻的溫度例子：房間溫度某一時(shí)刻的溫度：8KalmanFilter的實(shí)質(zhì)是一種數(shù)據(jù)處理算法利用測量數(shù)據(jù)來濾波數(shù)據(jù)濾波是去除噪聲還原真實(shí)數(shù)據(jù)的一種數(shù)據(jù)處理技術(shù)Kalman濾波在測量方差已知的情況下能夠從一系列存在測量噪聲的數(shù)據(jù)中，估計(jì)動(dòng)態(tài)系統(tǒng)的狀態(tài)“沒有時(shí)間把一件事情做好，卻有時(shí)間把一件事情反復(fù)做”KalmanFilter的實(shí)質(zhì)是一種數(shù)據(jù)處理算法9反饋控制法估計(jì)狀態(tài)反饋控制法估計(jì)狀態(tài)10符號慣例X：狀態(tài)變量U：輸入量（如電流）Z：測量值H：Z=H*X（H–系數(shù)）P：協(xié)方差K（Kg）：Kalman增益Q、R：估算與測量噪聲的方差符號慣例X：狀態(tài)變量11線性Kalman濾波：一般理論狀態(tài)方程 X(k)=A*X(k-1)+B*U(k)+W(k) Cov(W(k))=Q測量方程 Z(k)=H*X(k)+V(k) Cov(V(k))=RX：狀態(tài)變量U：輸入量Z：測量值A(chǔ)、B、H系數(shù)Q、R：噪聲線性Kalman濾波：一般理論狀態(tài)方程X：狀態(tài)變量12請記住從第k-1步到第k步估算值：下標(biāo)

k,k-1

最優(yōu)值：下標(biāo)

k,k或者k-1,k-1

請記住從第k-1步到第k步13Kalman方程（1-3）狀態(tài)預(yù)測

k|k-1表示“預(yù)測”

X(k|k-1)=A*X(k-1|k-1)+B*U(k)協(xié)方差預(yù)測

P(k|k-1)=A*P(k-1|k-1)*A’+Q狀態(tài)最優(yōu)

k|k（k-1|k-1）

表示“最優(yōu)”

X(k|k)=X(k|k-1)+Kg(k)*(Z(k)–H*X(k|k-1))Kalman增益協(xié)方差更新X：狀態(tài)變量U：輸入量P：協(xié)方差Kg：增益Z：測量值A(chǔ)、B、H系數(shù)Q：噪聲Kalman方程（1-3）狀態(tài)預(yù)測 k|k-1表示“預(yù)測14Kalman方程（4-5）狀態(tài)預(yù)測 X(k|k-1)=A*X(k-1|k-1)+B*U(k)協(xié)方差預(yù)測 P(k|k-1)=A*P(k-1|k-1)*A’+Q狀態(tài)最優(yōu)X(k|k)=X(k|k-1)+Kg(k)*(Z(k)–H*X(k|k-1))Kalman增益

Kg(k)=P(k|k-1)*H’/(H*P(k|k-1)*H’+R)協(xié)方差更新

P(k|k)=(I–Kg(k)*H)*P(k|k-1)X：狀態(tài)變量U：輸入量P：協(xié)方差Kg：增益Z：測量值A(chǔ)、B、H系數(shù)Q、R：噪聲Kalman方程（4-5）狀態(tài)預(yù)測 X(k|k-1)=15卡爾曼濾波法(Kalman濾波)用于SOC估算16Kalman方程：1維狀態(tài)預(yù)測 Xk,k-1=a*Xk-1+b*Uk協(xié)方差預(yù)測 Pk,k-1=a2*Pk-1+Q狀態(tài)最優(yōu) Xk=Xk,k-1+Kgk*(Zk–h*Xk,k-1)Kalman增益 Kgk=Pk,k-1*h/(h2*Pk,k-1+R)協(xié)方差更新 Pk=(I–Kgk*h)*Pk,k-1X：狀態(tài)變量U：輸入量P：協(xié)方差Kg：增益Z：測量值A(chǔ)、b、h系數(shù)Q、R：噪聲Kalman方程：1維狀態(tài)預(yù)測 Xk,k-1=a*Xk17X：估算溫度，Z：測量溫度(a=1,Uk=0,h=1)狀態(tài)預(yù)測 Xk,k-1=a*Xk-1+b*Uk=》Xk,k-1=Xk-1協(xié)方差預(yù)測Pk,k-1=a2*Pk-1+Q=》Pk,k-1=Pk-1+Q狀態(tài)最優(yōu) Xk=Xk,k-1+Kgk*(Zk–h*Xk,k-1) ==》 Xk=Xk,k-1+Kgk*(Zk–Xk,k-1)Kalman增益 Kgk=Pk,k-1*h/(h2*Pk,k-1+R) ==》 Kgk=Pk,k-1/(Pk,k-1+R)協(xié)方差更新 Pk=(I–Kgk*h)*Pk,k-1 ==》 Pk=(I–Kgk)*Pk,k-1例子：房間溫度X：狀態(tài)變量U：輸入量P：協(xié)方差Kg：增益Z：測量值Q、R：噪聲X：估算溫度，Z：測量溫度(a=1,Uk=0,18K,k-1表示“預(yù)測”，k或者k-1表示“最優(yōu)”狀態(tài)預(yù)測

Xk,k-1=Xk-1協(xié)方差預(yù)測

Pk,k-1=Pk-1+Q狀態(tài)最優(yōu)

Xk=Xk,k-1+Kgk*(Zk–Xk,k-1)Kalman增益

Kgk=Pk,k-1/(Pk,k-1+R)協(xié)方差更新

Pk=(I–Kgk)*Pk,k-1例子：房間溫度X：狀態(tài)變量P：協(xié)方差Kg：增益Z：測量值Q、R：噪聲K,k-1表示“預(yù)測”，k或者k-1表示“最優(yōu)”狀態(tài)預(yù)測19遞歸推演溫度方差初始值25測量值24.9R=0.25估算值25.1，Q=0.01，P(k-1)=0.1^2P(k,k-1)=P(k-1)+Q=0.02增益Kg(k-1)=P(k,k-1)/

(P(k,k-1)+R)=0.0741最優(yōu)值24.9+Kg(k-1)*

(25.1–24.9)=24.915P(k)=(1–Kg(k-1))*P(k,k-1)=0.0186K–1時(shí)刻溫度方差估算值24.915測量值25.5（更新值不需測量值）先計(jì)算P(k+1,k)及Kg(k)更新值P(k+1,k)=P(k)+QKg(k)=P(k+1,k)/

(P(k+1,k)+R)最優(yōu)值25.5+Kg(k)*(24.915–25.5)K時(shí)刻遞歸推演溫度方差初始值25測量值24.9R=0.25估算20卡爾曼濾波法(Kalman濾波)用于SOC估算21Kalman的初始估計(jì)誤差：收斂快Kalman的初始估計(jì)誤差：收斂快22清華專利介紹清華專利介紹23狀態(tài)方程與觀測方程狀態(tài)

SOCk+1

=SOCk+（?t

/

Q0）Ik

觀測

Ut=Uoc-UD-iLRi狀態(tài)方程與觀測方程狀態(tài)24狀態(tài)方程狀態(tài)

SOCk+1

=SOCk+（?t

/

Q0）Ik

X==》SOC

Xk+1= Xk+Ik*Dt/C （Ik<0） Xk+Ik*h*Dt/C （Ik>0）狀態(tài)方程狀態(tài)25電池等效電路Nernst模型電池等效電路Nernst模型26電池觀測方程Y(k)=K0–RIk–K1/xk–K2xk+K3ln(xk)+K4ln(1–xk)+VkY(k)是xk的非線性方程電池觀測方程Y(k)=K0–RIk–K1/x27Nernst模型系數(shù)K0=534.0017K1=2.6273K2=-131.7037K3=95.4526K4=-6.2601確定：放電實(shí)驗(yàn)+最小二乘法Nernst模型系數(shù)K0=534.001728對非線性的觀測方程做線性化Y(k)=f(Ik,xk)+Vkf(Ik,xk)對xk在某一時(shí)刻的xk0做泰勒展開其一次項(xiàng)系數(shù)為擴(kuò)展卡爾曼濾波（EKF）對非線性的觀測方程做線性化擴(kuò)展卡爾曼濾波（EKF）29Taylor展開f(x)=f(x0)+(x-x0)*f’(x0)+(x-x0)2*f’’(x0)+……略去高階小量(x-x0)2，得到 f(x)=F*X+GTaylor展開f(x)=f(x0)+(x-x0)30卡爾曼濾波方程Ck=K1/xk2–K2+K3/xk–K4/(1–xk)

取代了線性Kalman方程中的H卡爾曼增益Kk=Pk*Ck/(Ck2*Pk+Q)協(xié)方差Pk=（1–Kk

*Ck）*Pk-1遞推公式xk=xk-1+a*Kk*(yk-f(Ik,xk))

濾波放大倍數(shù)a

=5~20卡爾曼濾波方程Ck=K1/xk2–K2+K331卡爾曼濾波法(Kalman濾波)用于SOC估算32卡爾曼濾