基于多輪升維策略與改進(jìn)量子行為粒子群優(yōu)化算法的熱導(dǎo)率函數(shù)估計(jì)方法

席振翔， 于 帆

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

1 引 言

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

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

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

2 正問(wèn)題

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

(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í)間隱式離散方程，并由三對(duì)角陣算法求解即可獲得各節(jié)點(diǎn)處的溫度值。

3 反問(wèn)題

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

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

(4)

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

(5)

將熱導(dǎo)率離散為常數(shù)從而進(jìn)行反演的方法稱(chēng)為熱導(dǎo)率的函數(shù)估計(jì)方法(function estimation approach)，本文將待反演的溫度范圍均分成多個(gè)溫度區(qū)間，x的每一分量對(duì)應(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)算法，并通過(guò)實(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代的平均最好位置，°為兩向量對(duì)應(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ī)模。

對(duì)于最小化問(wèn)題，當(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)化速度因子和聚集度因子。

在尋找極小值的問(wè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)全局搜索能力，適合于處理高維問(wèn)題。采用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ū)間對(duì)熱導(dǎo)率進(jìn)行離散，這意味著更高的反演維度，這時(shí)即便采用上述MAQPSO算法，也難以在較少迭代次數(shù)下獲得有效結(jié)果。

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

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

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

即除了首輪反演中粒子位置需要在迭代前進(jìn)行隨機(jī)初始化以外，后續(xù)輪次反演開(kāi)始時(shí)均繼承上一輪反演搜索到的有用信息。以上策略基于材料熱導(dǎo)率隨溫度變化相對(duì)連續(xù)的特點(diǎn)，利用算法在低維空間中獲取的信息指導(dǎo)高維搜索，從實(shí)現(xiàn)上來(lái)說(shuō)則是在搜索過(guò)程中多次對(duì)粒子的位置信息進(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)的簡(jiǎn)化模型，在y=b面，取h=10 W/(m2·K)，Tf=300 K，而在y=0面進(jìn)行大功率加熱，其溫度隨時(shí)間變化為

(17)

(18)

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

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

(19)

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

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

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

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

當(dāng)σ=0,2%,4%和6%，反演結(jié)果的平均相對(duì)誤差分別為2.172%,3.058%,2.773%和3.157%?？梢钥闯觯惴ㄔ诓煌`差水平下都能得到有效結(jié)果，且對(duì)測(cè)量誤差的敏感度較低，隨著σ增加，計(jì)算結(jié)果偏離精確解的程度略有增加。

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

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

(20)

而算法反演時(shí)仍將熱導(dǎo)率看作分段線性函數(shù)，以此來(lái)產(chǎn)生系統(tǒng)誤差，在這種情況下，即便σ=0，式(4)的值也無(wú)法降低至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 多輪升維策略對(duì)算法的影響

將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é)果的平均相對(duì)誤差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í)，各算法均無(wú)法高效處理41個(gè)熱導(dǎo)率分量。隨著搜索范圍增大，剩余的目標(biāo)函數(shù)值也增大，即便在(0,1)的范圍中進(jìn)行搜索，各算法所得結(jié)果也均無(wú)法滿足工程應(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算法，其部分原因在于多輪升維策略使得該問(wèn)題對(duì)算法全局搜索能力的要求降低。但從耗時(shí)上考慮，MAQPSO算法的并行計(jì)算方式使其能夠利用多核處理器資源，相比于PSO和QPSO算法有著更短的計(jì)算時(shí)間。

5 結(jié) 論

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

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

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