尹偉石, 劉曉奇, 徐 軒

(長(zhǎng)春理工大學(xué) 理學(xué)院, 長(zhǎng)春 130022)

熱傳導(dǎo)反問(wèn)題在航空航天、 核能工程、 醫(yī)療檢測(cè)等領(lǐng)域應(yīng)用廣泛. 在熱傳導(dǎo)正問(wèn)題中, 溫度場(chǎng)是在物體熱物性系數(shù)、 邊界條件和初始條件已知的情況下, 由一些數(shù)值算法求解得到, 常用的方法有差分法、 邊界元法[1]等; 在熱傳導(dǎo)反問(wèn)題中, 先已知溫度觀測(cè)值, 然后用已知觀測(cè)值對(duì)未知的熱參數(shù)等條件進(jìn)行求解, 常見(jiàn)的方法主要有正則化方法[2-3]、 有限元[4]、 遺傳算法[5]等. 熱源參數(shù)識(shí)別問(wèn)題作為熱傳導(dǎo)反問(wèn)題的一個(gè)重要研究方向, 目前已被廣泛關(guān)注. 例如： Yang等[6]利用Tikhonov正則化方法對(duì)空間熱源問(wèn)題進(jìn)行了求解； Shah等[7]利用共軛梯度法結(jié)合伴隨問(wèn)題對(duì)一維熱傳導(dǎo)方程的瞬時(shí)熱源進(jìn)行了反演； Lee[8]將排斥粒子群算法應(yīng)用于熱傳導(dǎo)反問(wèn)題的熱源參數(shù)估計(jì). 近年來(lái), 隨著計(jì)算機(jī)技術(shù)的發(fā)展, 以Bayes推理為基礎(chǔ)的熱傳導(dǎo)反問(wèn)題求解方法備受關(guān)注. 例如： Zeng等[9]提出了一種非參數(shù)總體Monte Carlo的自適應(yīng)Bayes近似計(jì)算方法, 用于線性和非線性熱傳導(dǎo)反問(wèn)題的計(jì)算； Silva等[10]利用Bayes濾波器和Markov鏈Monte Carlo方法對(duì)時(shí)變邊界熱流進(jìn)行了估計(jì)； Jakkareddy等[11]利用Bayes方法結(jié)合液晶熱像儀對(duì)自然對(duì)流中的局部換熱系數(shù)進(jìn)行反演. 本文利用Bayes遺傳算法對(duì)熱源位置和產(chǎn)生時(shí)間參數(shù)進(jìn)行求解.

1 正問(wèn)題

二維瞬態(tài)熱傳導(dǎo)問(wèn)題的控制方程可表示為

(1)

其中:u表示溫度；x,y表示空間坐標(biāo);Ω表示求解域;k表示導(dǎo)熱系數(shù);a=k/(ρc)表示熱擴(kuò)散系數(shù),ρ和c為常數(shù), 分別表示材料密度和比熱容；t表示時(shí)間;f(x,y,t)為Ω內(nèi)的瞬時(shí)熱源, 大小隨時(shí)間變化. 該問(wèn)題的邊界條件和初始條件分別為

其中:Γ1和Γ2表示邊界,Γ1∪Γ2=Γ,Γ1∩Γ2=?;n表示邊界Γ的外法向量.

2 Bayes遺傳算法

2.1 Bayes理論

Bayes公式為

式中：p(α|d)為參數(shù)的后驗(yàn)概率密度函數(shù),α為模型參數(shù),d為觀測(cè)數(shù)據(jù)；p(α)為參數(shù)的先驗(yàn)密度函數(shù)；P(d|α)為條件概率密度；p(d)為積分常數(shù).

根據(jù)Bayes公式可得參數(shù)的后驗(yàn)概率密度函數(shù). 本文研究的未知參數(shù)后驗(yàn)概率密度函數(shù)[12]為

(5)

式中：α=(α1,α2,…,αm)為模型未知參數(shù), 共有m個(gè)；d=(d1，d2,…,dn)為觀測(cè)值, 共有n個(gè)；ui(α)表示與參數(shù)α相對(duì)應(yīng)的第i個(gè)預(yù)測(cè)值,εi=di-ui(α)表示測(cè)量誤差,i=1,2,…,n, 并且假設(shè)觀測(cè)誤差服從均值為零、 標(biāo)準(zhǔn)差為σ的正態(tài)分布N(0,σ2), 且相互獨(dú)立；λ為常數(shù), 與參數(shù)α無(wú)關(guān).

(6)

(7)

由于式(7)中的數(shù)據(jù)通常都是高維的, 因此其形式非常復(fù)雜, 直接對(duì)其進(jìn)行計(jì)算很困難, 本文利用遺傳算法對(duì)式(7)進(jìn)行刻畫(huà), 進(jìn)而得到式(6)的最大后驗(yàn)概率密度.

2.2 遺傳算法

遺傳算法是一種具有全局尋優(yōu)能力的自適應(yīng)算法, 其用種群作為函數(shù)解, 算法具有高效、 實(shí)用、 穩(wěn)定性強(qiáng)的特點(diǎn), 可以有效地用于解決非線性、 多峰值、 多目標(biāo)函數(shù)等問(wèn)題. 遺傳算法包括適應(yīng)度計(jì)算、 選擇、 交叉、 變異4個(gè)步驟, 其中, 交叉是群體即解的信息隨機(jī)交換的過(guò)程, 變異是對(duì)解添加隨機(jī)擾動(dòng)的過(guò)程.

1) 適應(yīng)度計(jì)算： 把目標(biāo)函數(shù)作為適應(yīng)度函數(shù), 計(jì)算每個(gè)個(gè)體的適應(yīng)度值.

2) 選擇： 選擇是根據(jù)個(gè)體適應(yīng)度值, 對(duì)有效解個(gè)體進(jìn)行一定比例的選擇保留. 本文求解的問(wèn)題是使式(6)達(dá)到最大, 因此, 在計(jì)算過(guò)程中需按一定比例保留使適應(yīng)度值達(dá)到最大的個(gè)體, 被保留的這些個(gè)體稱為有效解, 復(fù)制這些解使其傳遞到子代.

在選擇過(guò)程中, 有多種選擇策略, 常用的是輪盤(pán)賭法. 輪盤(pán)賭法是依概率按照Monte Carlo法對(duì)每個(gè)個(gè)體進(jìn)行選取, 在選取過(guò)程中, 可以對(duì)個(gè)體的適應(yīng)度值進(jìn)行線性縮放, 即種群中最差個(gè)體的目標(biāo)函數(shù)值減去種群中每個(gè)個(gè)體的目標(biāo)函數(shù)值：

fj=max{uj|j=1,2,…,m}-uj.

(8)

在選擇過(guò)程中, 每個(gè)個(gè)體的選擇概率可根據(jù)適應(yīng)度值確定：

(9)

為有效減輕選擇比例的壓力, 減小算法早熟的問(wèn)題, 可采用排序選擇策略, 其基本思想是： 將M個(gè)個(gè)體的適應(yīng)度值按從小到大的升序進(jìn)行排列, 記為1,2,…,m, 排序后的個(gè)體被選擇概率為

(10)

3) 交叉： 交叉是對(duì)選中的一對(duì)候選解, 將其所包含的信息按一定的概率進(jìn)行交流互換的過(guò)程, 該過(guò)程是隨機(jī)的, 用公式表示為

(11)

4) 變異： 變異是對(duì)交叉后的種群進(jìn)行隨機(jī)擾動(dòng), 用公式表示為

xj,D+1=xr1,D+F(xr2,D-xr3,D),

(12)

其中:xj,D+1,xr1,D,xr2,D,xr3,D為解個(gè)體;r1,r2,r3∈{1,2,…,Np},Np為種群數(shù)目;D為進(jìn)化代數(shù);F∈[0,1]為變異因子. 變異是為提高遺傳算法對(duì)解空間局部搜索能力, 豐富種群多樣性, 增大其達(dá)到全局最優(yōu)可能而進(jìn)行的操作, 變異通常是以非常小的預(yù)定概率進(jìn)行的.

通過(guò)上述分析, 可得Bayes推理與遺傳算法相結(jié)合用于熱源參數(shù)識(shí)別問(wèn)題的求解步驟如下:

1) 根據(jù)Bayes原理和熱傳導(dǎo)方程確定目標(biāo)函數(shù)；

2) 根據(jù)待估參數(shù)的先驗(yàn)信息確定參數(shù)的取值范圍；

3) 初始化產(chǎn)生一定規(guī)模的種群, 即初始解；

4) 計(jì)算每個(gè)個(gè)體的適應(yīng)度值；

5) 根據(jù)合適的選擇策略對(duì)種群進(jìn)行選擇保留；

6) 選擇適合的交叉策略, 對(duì)保留的種群進(jìn)行交叉操作；

7) 根據(jù)相應(yīng)的變異策略對(duì)種群進(jìn)行變異操作；

8) 執(zhí)行終止判斷準(zhǔn)則, 如果滿足終止準(zhǔn)則, 則選擇適應(yīng)度值高的個(gè)體作為參數(shù)識(shí)別問(wèn)題的最終解, 停止計(jì)算; 否則, 返回步驟4), 繼續(xù)計(jì)算, 直到滿足終止判別條件為止.

3 數(shù)值算例

下面以包含瞬時(shí)熱源的二維熱傳導(dǎo)方程為例, 分析Bayes遺傳算法的反演效果.

3.1 實(shí) 例

對(duì)于二維熱傳導(dǎo)方程：

(13)

假設(shè)熱源為瞬時(shí)熱源, 建立熱源的二維熱傳導(dǎo)模型. 設(shè)熱源在坐標(biāo)中心處, 熱源產(chǎn)生的時(shí)刻記為t=0, 觀測(cè)點(diǎn)B的位置為(x0,y0), 首次測(cè)量時(shí)刻為熱源產(chǎn)生后的t0時(shí)刻, 則待反演參數(shù)的真實(shí)解為α=(x0,y0,t0). 取熱源參數(shù)真實(shí)值橫坐標(biāo)x0=0.35 m, 縱坐標(biāo)y0=0.55 m, 第一次檢測(cè)距熱源開(kāi)始時(shí)間t0=0.3 s. 用x,y,t作為參數(shù)真實(shí)值x0,y0和時(shí)刻t0的估計(jì)值, 記為α1=(x,y,t).

設(shè)從首次測(cè)量開(kāi)始后每隔0.02 s測(cè)量一次, 直到0.68 s, 利用熱源參數(shù)真實(shí)值和式(13)計(jì)算得觀測(cè)點(diǎn)B的溫度值, 共20個(gè)計(jì)算數(shù)據(jù), 結(jié)果列于表1.

表1 觀測(cè)點(diǎn)不同時(shí)刻測(cè)量的溫度值

3.2 Bayes遺傳算法的參數(shù)反演結(jié)果

根據(jù)上述數(shù)據(jù), 利用Bayes遺傳算法對(duì)熱源參數(shù)進(jìn)行估計(jì). 根據(jù)正問(wèn)題的先驗(yàn)信息, 得到α1的分布范圍, 其中每個(gè)參數(shù)都服從均勻分布, 且相互獨(dú)立, 然后根據(jù)Bayes定理, 由式(5)得模型后驗(yàn)概率密度函數(shù)為

(14)

設(shè)種群規(guī)模Np=100, 交叉率為Pc=0.9, 變異率Pm=0.1, 按照遺傳算法迭代50次達(dá)到穩(wěn)定, 參數(shù)的估計(jì)值為：x=0.351,y=0.551,t=0.301.

Bayes遺傳算法迭代300次得到參數(shù)x,y,t的后驗(yàn)概率直方圖和迭代曲線分別如圖1和圖2所示, 其中, 參數(shù)x,y,t的后驗(yàn)概率大小分別用p(x),p(y),p(t)表示.

圖1 基于Bayes遺傳算法的模型參數(shù)后驗(yàn)概率直方圖

由圖1可見(jiàn), 參數(shù)x,y,t的后驗(yàn)概率分布呈近似正態(tài)分布, 并且分別在0.35,0.55,0.3附近頻率達(dá)到最大. 由圖2可見(jiàn), 當(dāng)?shù)螖?shù)超過(guò)50次時(shí), 參數(shù)x,y,t穩(wěn)定地收斂到參數(shù)真實(shí)值附近.

3.3 反演結(jié)果對(duì)比

作為對(duì)比實(shí)驗(yàn), 采用Bayes微分進(jìn)化算法對(duì)模型參數(shù)進(jìn)行求解, 得到參數(shù)x,y,t的反演結(jié)果如圖3所示. 由圖3可見(jiàn), 當(dāng)?shù)芜_(dá)到130次時(shí), 參數(shù)x,y,t穩(wěn)定地收斂到參數(shù)真實(shí)解附近, 結(jié)合上述對(duì)Bayes遺傳算法收斂速度的分析, 可知Bayes遺傳算法的收斂速度更快.

圖2 基于Bayes遺傳算法的模型參數(shù)迭代曲線

圖3 基于Bayes微分進(jìn)化算法的模型參數(shù)迭代曲線

為進(jìn)一步分析Bayes遺傳算法的反演效果, 從Bayes遺傳算法的反演結(jié)果中剔除前50次不穩(wěn)定的結(jié)果, 只對(duì)后250次反演結(jié)果進(jìn)行統(tǒng)計(jì)分析； 同理, 剔除Bayes微分進(jìn)化算法前130次不穩(wěn)定的反演結(jié)果, 只對(duì)后170次反演結(jié)果進(jìn)行統(tǒng)計(jì)分析, 統(tǒng)計(jì)結(jié)果列于表2. 由表2可見(jiàn), Bayes遺傳算法得到參數(shù)x,y,t的均值誤差都在1%以下, Bayes微分進(jìn)化算法得到參數(shù)x,y,t的均值誤差在4%以下, 說(shuō)明與Bayes微分進(jìn)化算法相比, Bayes遺傳算法的反演精度更高.

表2 Bayes遺傳算法與Bayes微分進(jìn)化算法的模型參數(shù)后驗(yàn)統(tǒng)計(jì)結(jié)果比較

為驗(yàn)證Bayes遺傳算法的穩(wěn)定性, 先對(duì)測(cè)量數(shù)據(jù)分別添加5%和10%的噪聲, 然后對(duì)反演結(jié)果每隔10個(gè)點(diǎn)取一個(gè)值, 得到參數(shù)x,y,t序列的相對(duì)誤差和統(tǒng)計(jì)分析結(jié)果, 分別如圖4和表3所示.

圖4 Bayes遺傳算法在不同噪聲下的參數(shù)反演相對(duì)誤差序列

表3 Bayes遺傳算法在不同噪聲下參數(shù)反演的統(tǒng)計(jì)分析結(jié)果

由圖4可見(jiàn), 對(duì)觀測(cè)數(shù)據(jù)分別添加5%和10%的噪聲, 參數(shù)的相對(duì)誤差均控制在9%以下. 由表3可見(jiàn), 當(dāng)對(duì)觀測(cè)數(shù)據(jù)添加5%的噪聲時(shí), 參數(shù)的均值誤差都控制在2%～5%, 當(dāng)對(duì)觀測(cè)數(shù)據(jù)添加10%的噪聲時(shí), 參數(shù)的均值誤差都控制在4%～9%, 說(shuō)明Bayes遺傳算法具有較好的穩(wěn)定性.

綜上可見(jiàn), Bayes遺傳算法用概率語(yǔ)言對(duì)參數(shù)信息進(jìn)行描述, 算法具有較好的穩(wěn)定性, 能有效地對(duì)熱源參數(shù)進(jìn)行估計(jì), 為基于Bayes方法的熱傳導(dǎo)反問(wèn)題提供了一種新的計(jì)算方法.