席振翔， 于 帆

(北京科技大學(xué) 能源與環(huán)境工程學(xué)院，北京 100083)

1 引 言

根據(jù)測得的溫度場反演材料隨溫度變化的熱導(dǎo)率是一類典型的熱傳導(dǎo)反問題IHCP(inverse heat conduction problem)，在各類工業(yè)領(lǐng)域廣泛存在[1,2]。IHCP屬于不適定問題，較小的輸入誤差可能造成反演結(jié)果精度驟降，不恰當(dāng)?shù)某踔狄矔?huì)使得算法無法獲得有效結(jié)果。受限于測溫技術(shù)，測量誤差無法避免，這就要求解決方案能在一定的輸入誤差下穩(wěn)定工作，同時(shí)，若解決方案過多依賴于初值或是搜索范圍的選取，則會(huì)大大降低其有效性。

對于熱導(dǎo)率反演問題，國內(nèi)外學(xué)者提出的求解方法主要可分為梯度類方法與隨機(jī)類方法。梯度類方法主要包括Levenberg-Marquardt算法[3]、伴隨方程法[4]和共軛梯度法[5]等，此類方法收斂速度快，計(jì)算結(jié)果精度高，但會(huì)涉及到靈敏度矩陣的計(jì)算，且初值敏感度高，無法跳出局部極值點(diǎn)。隨機(jī)類方法中應(yīng)用最廣泛的為遺傳算法[6]，該方法具有較強(qiáng)的全局搜索能力，對搜索范圍的要求相對較低，但收斂速度緩慢。在先驗(yàn)信息不足且同時(shí)反演的熱導(dǎo)率參數(shù)較多的情況下，即便是全局搜索能力較強(qiáng)的隨機(jī)類方法也難以得出有效結(jié)果。

為解決以上問題，本文將一種改進(jìn)的量子行為粒子群優(yōu)化算法應(yīng)用于高溫隔熱材料熱導(dǎo)率函數(shù)估計(jì)問題中，并提出了一種多輪升維策略提升算法在該問題的搜索效率。

2 正問題

為確定材料熱導(dǎo)率而開展的實(shí)驗(yàn)往往會(huì)設(shè)置一個(gè)加熱面，通過大功率加熱使得試件升溫，而非加熱面既可以安裝隔熱材料防止熱量散失，也可以使其暴露于空氣中。經(jīng)過一段時(shí)間的加熱，獲取試件上各測點(diǎn)的溫度數(shù)據(jù)，以此來反演材料隨溫度變化的熱導(dǎo)率。本文考慮材料除上下表面外的四個(gè)端面為絕熱，在一維變物性、無源項(xiàng)和瞬態(tài)導(dǎo)熱問題中討論算法性能，其控制方程可表示為

(1)

初始條件為

t=0∶T(y,t)=T0

(2)

邊界條件為y=0∶T=f(t)

(3)

由外節(jié)點(diǎn)法將目標(biāo)區(qū)域離散為有限個(gè)網(wǎng)格，采用有限差分法推導(dǎo)時(shí)間隱式離散方程，并由三對角陣算法求解即可獲得各節(jié)點(diǎn)處的溫度值。

3 反問題

3.1 反問題的優(yōu)化目標(biāo)函數(shù)

在反問題中，優(yōu)化目標(biāo)函數(shù)設(shè)置為測點(diǎn)溫度計(jì)算值與測量值之間的殘差，即

(4)

x=(x1,x2,…,xD)=(λ1,λ2,…,λD)

(5)

將熱導(dǎo)率離散為常數(shù)從而進(jìn)行反演的方法稱為熱導(dǎo)率的函數(shù)估計(jì)方法(function estimation approach)，本文將待反演的溫度范圍均分成多個(gè)溫度區(qū)間，x的每一分量對應(yīng)某一溫度區(qū)間端點(diǎn)的熱導(dǎo)率值，溫度區(qū)間內(nèi)熱導(dǎo)率值由溫度區(qū)間端點(diǎn)處熱導(dǎo)率值線性插值獲得。

3.2 改進(jìn)的量子行為粒子群優(yōu)化算法

粒子群優(yōu)化PSO(particle swarm optimization)算法是一種基于種群的高性能優(yōu)化技術(shù)[7]，在IHCP領(lǐng)域已經(jīng)得到了廣泛應(yīng)用，并在一些傳統(tǒng)方法失效的情況下得以證明，PSO算法也能取得良好的數(shù)值結(jié)果[8]。在標(biāo)準(zhǔn)的PSO基礎(chǔ)上，Sun等[9，10]從量子力學(xué)的角度提出了量子行為粒子群優(yōu)化QPSO(quantum-behavior particle swarm optimization)算法，并通過實(shí)驗(yàn)結(jié)果證明了QPSO算法收斂性能相比于PSO算法以及遺傳算法有了很大的提升。QPSO算法的粒子進(jìn)化方程可表示為

xi,k + 1=ai,k±α|mk-xi,k|°ln(1/ui,k)

(6)

式中xi,k為第i個(gè)粒子在第k代時(shí)的位置,α為收縮擴(kuò)張系數(shù),ui,k為0～1上均勻分布的隨機(jī)變量組成的向量,ai,k為局部吸引子,mk為第k代的平均最好位置，°為兩向量對應(yīng)元素分別相乘。

ai,k=φi,k°pi,k+(e-φi,k)°gk

(7)

式中pi,k為第i個(gè)粒子在第k代時(shí)的個(gè)體最佳位置，gk為第k次迭代時(shí)的種群最佳位置，φi,k為 0～1上均勻分布的隨機(jī)變量組成的向量，e為與φi,k同維的單位向量。

(8)

式中N為種群規(guī)模。

對于最小化問題，當(dāng)k>0時(shí)，pi,k與gk按如下方法確定，

(9)

(10，11)

本文式(6)的參數(shù)α采用自適應(yīng)策略[11]，相比于普遍采用的線性遞減策略[9]，該策略能使QPSO的收斂速度明顯提高。其實(shí)現(xiàn)方法如下，

αk=C1-C2sk+C3jk

(12)

式中C1，C2和C3均為常數(shù)，分別取為1，0.4和0.1；sk與jk分別為進(jìn)化速度因子和聚集度因子。

在尋找極小值的問題中，當(dāng)k>0時(shí)，sk與jk的計(jì)算方法為

sk=F(gk)/F(gk - 1)

(13)

(14)

為了增加算法的全局搜索能力，合理利用多核處理器資源，本文基于多種群協(xié)同搜索策略[12]構(gòu)建一種粗粒度并行模型，將整個(gè)粒子群劃分為多個(gè)子種群，每個(gè)子種群由多核CPU的一個(gè)核進(jìn)行獨(dú)立計(jì)算，再按照一定的交流概率共享部分信息。其交流概率計(jì)算為

qk=(k/K)z

(15)

式中K為迭代總次數(shù)，z取0.8。

采用以上改進(jìn)策略的QPSO算法即為多種群并行自適應(yīng)量子行為粒子群優(yōu)化算法，即MAQPSO算法，其具有較強(qiáng)全局搜索能力，適合于處理高維問題。采用MAQPSO算法反演熱導(dǎo)率的計(jì)算流程如下。

(1) 置k=0，在搜索范圍(xmin,xmax)內(nèi)初始化各粒子群中粒子位置xi,0，置pi,0=xi,0，由式(8)計(jì)算m0，由式(4)計(jì)算F(xi,0)，由式(10，11)確定g0。

(2) 由式(15)計(jì)算qk，產(chǎn)生(0,1)內(nèi)均勻分布的隨機(jī)數(shù)μk，若μk

(3) 若k=0，置α0=1；否則，分別由式(13，14)計(jì)算進(jìn)化速度因子sk與聚集度因子jk，再由式(12)計(jì)算收縮擴(kuò)張因子αk。

(4) 由式(7)計(jì)算局部吸引子ai,k，由式(6)更新粒子位置xi,k + 1，由式(4)計(jì)算F(xi,k + 1)，由式(9～11)確定pi,k + 1與gk + 1，由式(8)計(jì)算mk + 1。

(5) 若k+1

3.3 基于熱導(dǎo)率函數(shù)估計(jì)的多輪升維策略

在反演較大溫度跨度下的材料熱導(dǎo)率時(shí)，為了更細(xì)致地捕捉熱導(dǎo)率隨溫度的變化規(guī)律，往往會(huì)劃分更小的溫度區(qū)間對熱導(dǎo)率進(jìn)行離散，這意味著更高的反演維度，這時(shí)即便采用上述MAQPSO算法，也難以在較少迭代次數(shù)下獲得有效結(jié)果。

熱導(dǎo)率反演也存在初值敏感度問題，梯度類算法雖然有更高的計(jì)算效率，但其計(jì)算表現(xiàn)卻嚴(yán)重依賴于初值的選取，且初值敏感度會(huì)隨著反演參數(shù)個(gè)數(shù)以及測量誤差水平的增加而增加。相比之下，基于種群的隨機(jī)類算法僅需輸入一個(gè)搜索范圍，這使得此類方法在材料先驗(yàn)信息不足的情況下有著天然的優(yōu)勢。然而在較高的反演維度下，過大的搜索范圍也會(huì)使得隨機(jī)類算法收斂速度大大降低。

為了解決以上問題，本文將MAQPSO的反演過程劃分為M輪(記m為反演輪次，dm與Km分別為第m輪反演的維度與最大迭代次數(shù))，dm從低到高逐輪提升，并在最后一輪開始時(shí)達(dá)到目標(biāo)維度。在每輪反演結(jié)束后將粒子的位置x通過線性插值轉(zhuǎn)化為下一輪反演的粒子初始位置，而個(gè)體最佳位置p與種群最佳位置g也做同樣的處理。

圖1 多輪升維策略實(shí)現(xiàn)方式Fig.1 Multi -round upgrading strategy implementation

即除了首輪反演中粒子位置需要在迭代前進(jìn)行隨機(jī)初始化以外，后續(xù)輪次反演開始時(shí)均繼承上一輪反演搜索到的有用信息。以上策略基于材料熱導(dǎo)率隨溫度變化相對連續(xù)的特點(diǎn)，利用算法在低維空間中獲取的信息指導(dǎo)高維搜索，從實(shí)現(xiàn)上來說則是在搜索過程中多次對粒子的位置信息進(jìn)行插值，使得搜索進(jìn)程由低維向高維推進(jìn)。

4 仿真實(shí)驗(yàn)與討論

4.1 仿真實(shí)驗(yàn)條件與算法設(shè)置

c(T)=1013+0.75×10-3T2

(16)

用離散數(shù)據(jù)點(diǎn)表示在300 K～1300 K的范圍內(nèi)熱導(dǎo)率隨溫度的變化，計(jì)算中將熱導(dǎo)率看作分段線性函數(shù)，每25 K取一個(gè)溫度節(jié)點(diǎn)，由41個(gè)離散數(shù)據(jù)點(diǎn)組成的向量即為待反演變量x的精確解。

一維模型為真實(shí)實(shí)驗(yàn)的簡化模型，在y=b面，取h=10 W/(m2·K)，Tf=300 K，而在y=0面進(jìn)行大功率加熱，其溫度隨時(shí)間變化為

(17)

(18)

在較大的搜索范圍下測試MAQPSO算法的搜索能力，取xmax=100，由于真實(shí)熱導(dǎo)率值總大于0，取xmin=0，采用4個(gè)種群并行計(jì)算，每個(gè)種群的粒子數(shù)為10?；诙噍喩S策略，將算法的反演過程分為4輪，d1～d4分別取2，6，21和41，K1～K4分別取30，30，30和50。

定義平均相對誤差，以描述反演結(jié)果的誤差水平為

(19)

本文數(shù)值實(shí)驗(yàn)環(huán)境為MATLAB R2019a，Intel(R) Core(TM) i7-6700HQ CPU@2.60GHz，為了便于對比，每次計(jì)算均固定ω，即在同一套測量噪音下進(jìn)行數(shù)值實(shí)驗(yàn)。

4.2 測量誤差及系統(tǒng)誤差下的反演結(jié)果

圖2 不同σ下熱導(dǎo)率的反演結(jié)果Fig.2 Inversion results of thermal conductivity with different σ

圖3 不同σ下利用反演結(jié)果計(jì)算出來的測點(diǎn)溫度變化歷程與其測量值比較Fig.3 Comparison of temperature history calculated with estimated thermal conductivity and experimental value under different σ

當(dāng)σ=0,2%,4%和6%，反演結(jié)果的平均相對誤差分別為2.172%,3.058%,2.773%和3.157%。可以看出，算法在不同誤差水平下都能得到有效結(jié)果，且對測量誤差的敏感度較低，隨著σ增加，計(jì)算結(jié)果偏離精確解的程度略有增加。

將熱導(dǎo)率真值設(shè)為分段函數(shù)，忽略了系統(tǒng)誤差的影響，理想情況下，若σ=0，式(4)的值可降低至0。為研究系統(tǒng)誤差帶來的影響，熱導(dǎo)率真值改為由以下隨溫度變化的連續(xù)函數(shù)獲得

λ(T)=0.5exp(T/1500)+0.2sin(T/150)

(20)

而算法反演時(shí)仍將熱導(dǎo)率看作分段線性函數(shù)，以此來產(chǎn)生系統(tǒng)誤差，在這種情況下，即便σ=0，式(4)的值也無法降低至0，精確解與反演結(jié)果如圖4所示。

圖4 不同σ下λ(T)的反演結(jié)果Fig.4 Inversion results of λ(T) with different σ

可以看出，在存在系統(tǒng)誤差的情況下，算法仍能較好地捕捉熱導(dǎo)率隨溫度的變化，反演結(jié)果在較低σ下與真值的函數(shù)曲線符合良好。隨著σ增加，剩余的目標(biāo)函數(shù)值增加，當(dāng)σ=0，2%，4%和6%，反演結(jié)果的目標(biāo)函數(shù)值分別為0.092，4.591，9.180和13.771，反演結(jié)果逐漸散布在了真值兩側(cè)。

4.3 多輪升維策略對算法的影響

將MAQPSO算法與PSO和QPSO算法進(jìn)行比較，分別取xmax為1，10，100和1000，σ=2%，總粒子數(shù)為40，迭代總次數(shù)為140，熱導(dǎo)率真值設(shè)為以分段線性函數(shù)表示的材料真實(shí)熱導(dǎo)率，算法將反演300 K～1300 K下41個(gè)離散熱導(dǎo)率參數(shù)值(圖2 中Exact)，按是否采用多輪升維策略分為兩組，若采用該策略，各輪次迭代的反演維度與次數(shù)設(shè)置同4.1節(jié)。

PSO算法為采用了慣性權(quán)重線性遞減策略的標(biāo)準(zhǔn)粒子群算法，其慣性權(quán)重值隨著迭代次數(shù)的增加從0.9線性減小至0.6，設(shè)置粒子的運(yùn)動(dòng)速度上限vmax=0.1xmax，學(xué)習(xí)因子c1=c2=1.5。

QPSO算法中α采用線性遞減策略，其值隨著迭代次數(shù)的增加從0.8線性減小至0.5；MAQPSO算法的設(shè)置與4.1節(jié)相同。

兩組數(shù)值實(shí)驗(yàn)中，同一搜索范圍下各算法均采用同一隨機(jī)數(shù)種子初始化種群，使得各算法在同一搜索范圍下有相同的起點(diǎn)。不同搜索范圍下的初始目標(biāo)函數(shù)值Finitial、最終目標(biāo)函數(shù)值Ffinal、反演結(jié)果的平均相對誤差Emr以及算法總運(yùn)行時(shí)間列入表1和表2。

表1 未采用多輪升維策略時(shí)各算法的反演結(jié)果Tab.1 Inversion results of each algorithm without using multi-round upgrading strategy

表2 采用多輪升維策略時(shí)各算法的反演結(jié)果Tab.2 Inversion results of each algorithm using multi-round upgrading strategy

由表1可知，在未采用多輪升維策略時(shí)，各算法均無法高效處理41個(gè)熱導(dǎo)率分量。隨著搜索范圍增大，剩余的目標(biāo)函數(shù)值也增大，即便在(0,1)的范圍中進(jìn)行搜索，各算法所得結(jié)果也均無法滿足工程應(yīng)用的要求。

由表2可知，采用了多輪升維策略后，各算法的搜索效率明顯提高，在(0,1)的搜索范圍下，PSO、QPSO和MAQPSO算法均在140次迭代后獲得了有效結(jié)果，如圖5所示。可以看出，即便是采用了多輪升維策略，隨著搜索范圍的加大，PSO算法的結(jié)果質(zhì)量也逐漸下滑，當(dāng)xmax=1000時(shí)，僅有采用了多輪升維策略的QPSO與MAQPSO算法能在140次迭代中獲得有效結(jié)果，且PSO算法在使用時(shí)需要設(shè)置粒子運(yùn)動(dòng)速度范圍，不恰當(dāng)?shù)膙max也會(huì)影響算法收斂，這限制了傳統(tǒng)粒子群優(yōu)化算法的工程應(yīng)用。采用多輪升維策略的QPSO與MAQPSO算法在不同搜索范圍下的計(jì)算結(jié)果如圖6和圖7所示，可以看出，其反演結(jié)果幾乎不受搜索范圍的影響。

圖5 不同算法在(0,1)搜索范圍下的反演結(jié)果Fig.5 Thermal conductivity inversion results obtained by different algorithms in the search range of (0,1)

圖6 基于多輪升維策略的QPSO算法在不同xmax下的反演結(jié)果Fig.6 Thermal conductivity inversion results of QPSO algorithm based on multi-round upgrading strategy at different xmax

圖7 基于多輪升維策略的MAQPSO算法在不同xmax下的反演結(jié)果Fig.7 Thermal conductivity inversion results of MAQPSO algorithm based on multi -round upgrading strategy at different xmax

由表2可知，當(dāng)采用了多輪升維策略后，MAQPSO算法的反演精度略低于QPSO算法，其部分原因在于多輪升維策略使得該問題對算法全局搜索能力的要求降低。但從耗時(shí)上考慮，MAQPSO算法的并行計(jì)算方式使其能夠利用多核處理器資源，相比于PSO和QPSO算法有著更短的計(jì)算時(shí)間。

5 結(jié) 論

本文將一種改進(jìn)的量子行為粒子群算法應(yīng)用于熱導(dǎo)率函數(shù)估計(jì)問題中，并提出了一種多輪升維策略提升粒子群優(yōu)化算法在該問題中的搜索性能。通過數(shù)值實(shí)驗(yàn)討論了測量誤差以及系統(tǒng)誤差對算法反演結(jié)果的影響，并分析了不同粒子群優(yōu)化算法在該問題下的表現(xiàn)，得到的主要結(jié)論如下。

(1) 基于多輪升維策略的MAQPSO算法能在單個(gè)測點(diǎn)下同時(shí)反演多個(gè)熱導(dǎo)率參數(shù)，對搜索范圍以及測量誤差的敏感度較低。

(2) 多輪升維策略對PSO、QPSO和MAQPSO算法在熱導(dǎo)率函數(shù)估計(jì)問題的計(jì)算表現(xiàn)均有較大提升，但使用該策略后的PSO算法對搜索范圍仍然較為敏感，相比之下，QPSO和MAQPSO算法有更佳的計(jì)算表現(xiàn)。