基于擴(kuò)展Kalman濾波的熱耗散系數(shù)識別

王玉嬋,梅凌霜

南京信息工程大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院,江蘇 南京 210044

由于計算機(jī)的CPU高速運行時不斷產(chǎn)生熱量,因此需要在裝置中安裝電扇以保證CPU在正常工作的溫度范圍內(nèi),其工作原理見圖1[1]。數(shù)學(xué)上,上述熱傳導(dǎo)過程可以簡化為下面的一維熱傳導(dǎo)模型：

圖1 CPU熱對流模型[1]

其中a2、φ(t)分別表征了介質(zhì)內(nèi)的熱擴(kuò)散率及邊界x=0上的熱流；σ(t)是邊界上的阻尼系數(shù),描述了邊界x=d與外界的熱耗散能力。

CPU不斷產(chǎn)生的熱量,通過電扇等散熱裝置將熱量散發(fā)出去,確保高速運轉(zhuǎn)的CPU正常工作。一旦散熱裝置突發(fā)故障時,導(dǎo)致CPU無法及時散熱,這對高速運轉(zhuǎn)的CPU的是非常危險的,因此實時檢測CPU的工作溫度是非常重要且有意義的。但是實際中受客觀條件的限制,邊界x=d的溫度無法直接測量。事實上,CPU正常工作時,耗散系數(shù)在數(shù)學(xué)上是關(guān)于時間的光滑函數(shù),但是突發(fā)故障可能導(dǎo)致邊界x=d上的耗散系數(shù)發(fā)生突變,在數(shù)學(xué)上意味著耗散系數(shù)σ(t)為非光滑函數(shù)。因此待識別的內(nèi)部邊界耗散系數(shù)在數(shù)學(xué)上具有不同的光滑性質(zhì)。即需要借助于其他間接可測量數(shù)據(jù),例如外邊界的測量信息,來識別內(nèi)部具有不同光滑性的邊界耗散系數(shù)的不確定信息,進(jìn)而用戶可以根據(jù)檢測信息來控制電子器件工作。

基于上述實際背景,本文考慮利用邊界x=0上的溫度數(shù)據(jù)來識別邊界x=d時熱耗散系數(shù)σ(t)的不確定信息,即概率密度函數(shù)的統(tǒng)計信息。該問題是一個非線性的反邊值問題,屬于拋物方程反問題。WANG等[2]利用帶確定噪音模型的邊界溫度數(shù)據(jù)來反演非光滑的熱耗散系數(shù)的單值估計,證明了反問題的唯一性,并基于全變差罰項(TV)構(gòu)造了正則化重建方案。不同于確定誤差模型關(guān)注未知量的單值估計,隨機(jī)誤差模型側(cè)重評估未知變量的不確定性。近年來,有學(xué)者開始關(guān)注隨機(jī)誤差模型下相關(guān)反問題的研究[3]。而(擴(kuò)展)Kalman濾波法是研究一類非穩(wěn)態(tài)問題的廣泛方法[4]；IGLESIAS等[5]建立了Kalman濾波法處理線性反問題的收斂性,進(jìn)一步系統(tǒng)研究了Kalman濾波法與經(jīng)典Tikhonov正則化方法之間的聯(lián)系。而對于非線性的問題,無法直接使用Kalman濾波法處理,當(dāng)非線性的問題在局部內(nèi)可以線性化時,可以用擴(kuò)展Kalman濾波法來處理該非線性問題。因此,基于上述兩方面原因,本文考慮在隨機(jī)誤差模型下,基于擴(kuò)展Kalman濾波反演熱耗散系數(shù)。

此外,逆時問題[6]、源項識別問題[7]、期權(quán)隱含波動率識別問題[8]等都屬于拋物方程反問題,目前已經(jīng)了比較廣泛的研究成果,更多相關(guān)理論分析和重構(gòu)方案可參考文獻(xiàn)[9]。

本文首先基于問題(1)的伴隨問題,將反問題線性化,建立近似線性的動力系統(tǒng),再根據(jù)待識別耗散系數(shù)的不同光滑性,結(jié)合不同的空間先驗信息設(shè)計兩個反演算法,分別反演光滑的和分片光滑的耗散系數(shù)。

1 反問題描述與離散動力系統(tǒng)

1.1 反問題描述

在熱傳導(dǎo)模型(1)中,設(shè)附加數(shù)據(jù)為邊界x=0上帶隨機(jī)誤差的噪音數(shù)據(jù)

hε(t):=uε(0,t)=u(0,t)+ζ(t),

(2)

其中ζ(t)是服從零均值、方差ε2的高斯正態(tài)分布的隨機(jī)變量。

設(shè)熱傳導(dǎo)模型(1)中(a>0,f,φ,u0)已知,本文研究的反問題為,利用隨機(jī)噪音數(shù)據(jù)hε(t),識別邊界x=d上熱耗散系數(shù)的概率分布信息

σγ(t):=σ(t)+η(t),

(3)

其中η(t)是服從零均值、方差γ2的高斯正態(tài)分布的隨機(jī)變量。

1.2 離散動力系統(tǒng)

反問題(1)～(3)為非線性的發(fā)展模型,本文采用擴(kuò)展Kalman濾波方法求解,因此本小節(jié)首先將反問題(1)～(3)轉(zhuǎn)換為非線性動力系統(tǒng)

其中:映射Fk+1:σk+1→u[σk+1](0,t)是可微的,狀態(tài)噪聲Wk和觀測噪聲Vk+1滿足Gauss型分布,零均值且相互獨立,假設(shè)隨機(jī)變量σ0的初始條件密度函數(shù)已知。

Fk+1(σk+1)≈Fk+1(σk+1|k)+δw[σk+1|k,nk](0,t),

(5)

接下來將離散的公式(5)改寫為矩陣形式,定義方陣

則動力系統(tǒng)(4)中的第二個式子在σk+1|k局部,沿nk方向近似為線性方程：

hk+1=Fk+1(σk+1)+Vk+1≈Vk+1+Fk+1(σk+1|k)+DFk+1(σk+1|k)δ

≈Vk+1+Fk+1(σk+1|k)+DFk+1(σk+1|k)(σk+1-σk+1|k),Vk+1～N(0,ε2I),k=0,1,2,…。

(7)

此外,根據(jù)待識別熱耗散系數(shù)σ的光滑性改進(jìn)動力系統(tǒng)(4)。

參考文獻(xiàn)[4]中4.4節(jié)構(gòu)造L2型空間先驗函數(shù)的思想,將動力系統(tǒng)(4)中的第一個式子修正為

σk+1=G-1σk+Uk,k=0,1,2,…,

(8)

其中G=I+α2LTL,Uk～N(0,α2G-1),σk表示對應(yīng)向量t的阻尼系數(shù)的向量

至此,本文已建立了熱耗散系數(shù)σ的局部近似線性離散動力系統(tǒng)(7)～(8)。

又熱耗散系數(shù)σ為非光滑函數(shù),故結(jié)合TV正則化和交替迭代的思想,在離散系統(tǒng)(7)～(8)的基礎(chǔ)上,極小化TV罰項的近似泛函：

這是下一節(jié)擴(kuò)展Kalman濾波法算法的基礎(chǔ)。

2 擴(kuò)展Kalman濾波法的離散格式

基于Kalman濾波法可以求解線性動力系統(tǒng)(7)～(8)。參考文獻(xiàn)[4]的4.2節(jié),并利用定理3.8[4]、定理4.3[4],簡單計算直接可得下面定理。

定理1設(shè)動力系統(tǒng)系統(tǒng)(7)～(8)構(gòu)成一個發(fā)展-觀測模型。所有測量值集合記為Dk={h1,h2,…,hk},狀態(tài)向量σk的條件概率記π(σk|Dk),并記σk|1=E(σk|D1),Γk|1=cov(σk|D1),并定義π(σ0)=π(σ0|D0),則近似離散動力系統(tǒng)(8)和(7)的Kalman濾波更新公式具有如下形式：

(1)公式(8)更新：假設(shè)已知Gauss密度函數(shù)

π(σk|Dk)～N(σk|k,Γk|k),

那么π(σk+1|Dk)～N(σk+1|k,Γk+1|k),其中

矩陣L同(8)式中L。

(2)公式(7)更新：假設(shè)已知Gauss密度函數(shù)

π(σk+1|Dk)～N(σk+1|k,Γk+1|k),

那么π(σk+1|Dk+1)～N(σk+1|k+1,Γk+1|k+1),其中

σk+1|k=E(σk+1|Dk)=E(σk+1|σk),

Γk+1|k=cov(σk+1|Dk)=cov(σk+1|σk)。

基于定理1中動力系統(tǒng)的更新公式,給出反問題(1)～(3)基于Kalman濾波法識別非光滑σγ的算法1：具體計算步驟如下：

Step 1：給定α2、ε2、τ,初始迭代σ0|0、Γ0|0,最大迭代次數(shù)N和NTV,令k=0。

Step 3：σ=σk+1|k時求解定解問題(1),可得Fk+1(σk+1|k)=u[σk+1|k](d,t),再由定解問題(6)求解w[σk+1|k,nk](0,t),再計算方陣

Step 4：由公式(11)計算σk+1|k+1,Γk+1|k+1。

Step 5:由公式：

更新NTV次,并令k=k+1。

3 數(shù)值實現(xiàn)

本小節(jié)基于第二節(jié)中的算法進(jìn)行數(shù)值實驗,算法中的正問題和伴隨問題采用有限差分法計算。

設(shè)u(x,t)滿足下面的定解問題：

其精確解為：u(x,t)=sin(x)exp(-a2t)+x2+2a2t。

例1 設(shè)精確的阻尼系數(shù)與測量數(shù)據(jù)滿足：

hε(t):=u[σ](0,t)+N(0,ε2)。

本算例利用算法1來識別噪音水平ε2=1下阻尼系數(shù)的均值σk+1|k+1和方差Γk+1|k+1。

不同參數(shù)t、a2反演的σN|N和不同迭代步數(shù)反演的σk|k(0.5)分別見圖2和圖3。從圖2的紅線和黑線來看,改進(jìn)的(8)式確實可以降低擾動,但是不能較好的反演不光滑的角點；再與綠線對比發(fā)現(xiàn),(9)式的引入一方面降低擾動另一方面較好的保留了反演系數(shù)的不光滑性。從圖3來看,開始迭代時反演的σk|k(0.5)擾動著接近精確解,并且在較少步數(shù)時已得到較好的數(shù)值解。數(shù)值算例表明本文的算法是有效的。

圖2 不同τ、α2反演的σN|N

圖3 不同τ、α2反演的σk|k(0.5)

4 結(jié) 論

從數(shù)值算例來看,沒有引入TV罰項(9)式,反演的參數(shù)有較大擾動；引入TV罰項后反演參數(shù)(綠色)有了顯著提升。同時與該算法對噪音數(shù)據(jù)在較大的噪音方差下依然可以較好的反演。說明上述算法是有用的。