張 濤， 盧 玫， 陶 亮， 李博漢

（上海理工大學(xué)能源與動力工程學(xué)院，上海 200093）

基于粒子群優(yōu)化算法的尋源導(dǎo)熱反問題研究

張 濤， 盧 玫， 陶 亮， 李博漢

（上海理工大學(xué)能源與動力工程學(xué)院，上海 200093）

給出了描述導(dǎo)熱問題的數(shù)學(xué)模型，根據(jù)最小二乘法原理建立導(dǎo)熱反問題的目標(biāo)函數(shù)，并采用從鳥群捕食行為演化而來的粒子群優(yōu)化算法對含有熱源項的導(dǎo)熱反問題進行熱源位置的反演求解，同時對粒子群優(yōu)化算法中慣性系數(shù)的取值范圍進行了討論.結(jié)果表明：采用粒子群算法反演熱源的位置可以取得較好的結(jié)果，使用隨迭代次數(shù)變化的慣性系數(shù)可以加快算法的收斂速度.

熱傳導(dǎo)；反問題；熱源；粒子群優(yōu)化算法

尋源導(dǎo)熱反問題是指利用已知研究對象邊界上或內(nèi)部若干測點的溫度，反演求解熱源位置或熱源強度.在生物醫(yī)學(xué)、航空航天、化工制藥、熱工測量、無損探傷、電子散熱等領(lǐng)域都會涉及到此類反問題，如探尋電腦芯片中大量熱產(chǎn)生的熱源問題、微波爐的傳熱過程、化學(xué)反應(yīng)過程如何確定反應(yīng)產(chǎn)熱的熱源問題等［1］.反問題一般具有不適定性，較小的測量誤差就可能導(dǎo)致問題多解或者無解，因此導(dǎo)熱反問題成為許多學(xué)者的研究重點.導(dǎo)熱反問題的研究主要有熱物性參數(shù)的識別［2］、邊界形狀的識別［3］、邊界條件的識別［4］以及對各種反演算法的研究等.關(guān)于熱源項的反演研究相對比較少.Huang等［5］利用共軛梯度法反演熱源強度大小，Silva等［6］用同樣方法研究了兩個熱源的強度反演問題.早期求解導(dǎo)熱反問題的算法主要有正則化法、共軛梯度法、高斯牛頓法等.Huang［7］采用共軛梯度法識別二維一階非定常導(dǎo)熱問題的不規(guī)則幾何形狀，指出了共軛梯度法抗噪性能較好，但收斂速度較慢；楊海天等［8］采用高斯牛頓法對非線性導(dǎo)熱反問題進行了求解.近年來隨著計算機技術(shù)和仿生學(xué)的發(fā)展，一些學(xué)者將遺傳算法、神經(jīng)網(wǎng)絡(luò)法、蟻群算法等引入導(dǎo)熱反問題領(lǐng)域，取得了很好的結(jié)果，豐富了導(dǎo)熱反問題的求解方法.粒子群優(yōu)化算法（PSO）是由美國的Kenney和Eberhart［9］在1995年提出的，該算法是以粒子對群體中最優(yōu)粒子的追隨進行解空間的搜索.PSO的優(yōu)點在于程序流程簡單易實現(xiàn)、算法參數(shù)少易調(diào)整、收斂速度快、且有很多措施可以避免陷入局部最優(yōu).因此PSO從出現(xiàn)后，很快應(yīng)用到函數(shù)優(yōu)化、模糊控制等領(lǐng)域.

1 含內(nèi)熱源的二維穩(wěn)態(tài)導(dǎo)熱問題數(shù)學(xué)模型

式中，λ為導(dǎo)熱系數(shù)，W／（m·K）；T為溫度，K；Φ為點熱源強度，W，其對應(yīng)的熱源位置為（xΦ，yΦ）；Г為求解物體的邊界，Г=Г1+Г2+Г3為第一類邊界條件溫度分布，K；q為第二類邊界條件的熱流密度，（W／m2）；h為第三類邊界條件的表面對流換熱系數(shù)，W／（m2·K）；Tw，Tf分別為壁面溫度和外界流體的溫度，K.

2 反演模型及粒子群優(yōu)化算法

2.1 反演模型

對于反演熱源位置的導(dǎo)熱反問題，其熱源位置的確定必須附加物體表面或者內(nèi)部測點的測量信息來確定.其求解可以歸結(jié)為優(yōu)化問題

式中，N為溫度測點數(shù)目；T0（x，y）為邊界測點的實測溫度；Tc（x，y）為反演的熱源位置代入方程式（1）、邊界條件式（2）求得的對應(yīng)測點的計算值.目標(biāo)函數(shù)J隨反演過程中搜尋的熱源位置（xi，yj）而變化，當(dāng)（xi，yj）趨于真實熱源位置（xΦ，yΦ）時，式（3）所示目標(biāo)函數(shù)就趨于最小值，此時熱源位置（xi，yj）即為反演所求結(jié)果.

2.2 粒子群優(yōu)化算法

PSO來源于對鳥群搜尋食物的行為研究.鳥群在一個區(qū)域內(nèi)隨機搜索食物源時，所有鳥都不知道食物源在哪里，但是它們知道當(dāng)前離食物源最近的鳥在哪里，所以找到食物最好的策略就是搜尋目前離食物源最近鳥的周圍區(qū)域.PSO即是從這種模型中得到啟示并用于解決優(yōu)化問題.在PSO中，搜尋食物的鳥被稱為“粒子”.粒子飛行到一個位置后都有一個描述該位置距離食物源遠(yuǎn)近的適應(yīng)值，并且會根據(jù)最優(yōu)粒子的位置調(diào)整飛行到下一位置的速度矢量，以此追隨當(dāng)前的最優(yōu)粒子在解的空間中搜索.如在一個D維區(qū)域內(nèi)，由m個粒子構(gòu)成一個群體.設(shè)zi=（zi1，zi2，…，ziD）為第i個粒子（i=1，2，…，m）飛行到某一位置的空間坐標(biāo)，根據(jù)事先設(shè)定的目標(biāo)函數(shù)計算zi當(dāng)前的適應(yīng)值，以衡量粒子所在位置的優(yōu)劣；vi=（vi1，vi2，…，vid，…，viD）為粒子i的下一次的飛行速度，即粒子下一次的飛行距離和方向；每個粒子迄今為止找到的最優(yōu)位置為pbi=（pi1，pi2，…，pid，…，piD）；整個群體迄今為止找到的最優(yōu)解位置為gb=（g1，g2，…，gd，…，gD）.

對于所要研究的尋源導(dǎo)熱反問題，熱源位置等同于鳥群捕食的食物源；導(dǎo)熱區(qū)域為搜索熱源的空間范圍；式（3）所示目標(biāo)函數(shù)值的大小表示鳥離食物源的距離.某個粒子當(dāng)前已飛過的所有位置中，對應(yīng)目標(biāo)函數(shù)最小的那個位置為該粒子的個體極值，即pbi；整個粒子群所找到的對應(yīng)目標(biāo)函數(shù)最小的位置為全局極值，即gb.粒子通過更新方程改變每次搜尋的速度，進行下一次迭代（即飛行）.

PSO的粒子速度更新方程為

式中，k為當(dāng)前的迭代代數(shù)；c1，c2為學(xué)習(xí)因子；r1，r2為［0，1］范圍內(nèi)的隨機數(shù).迭代終止條件設(shè)為最大迭代次數(shù)或者搜索到滿足精度要求的最優(yōu)解.所示PSO速度更新方程中，等號右邊由三項組成：第一項表示粒子當(dāng)前的速度，第二項是根據(jù)粒子本身的個體極值對速度的調(diào)整，第三項是根據(jù)粒子群的全局極值對速度的調(diào)整.

在基本PSO的基礎(chǔ)上，通過添加慣性系數(shù)ω來調(diào)整粒子當(dāng)前速度對下一次飛行速度的影響程度.粒子速度更新方程為

其中，ω的取值可以為常數(shù)，亦可根據(jù)迭代次數(shù)的增加而減小，即計算式為

式中，ωmax，ωmin為最大和最小慣性系數(shù)；k，kmax為當(dāng)前迭代次數(shù)和最大迭代次數(shù).

2.3 尋源反演過程

圖1所示，為所編制的采用PSO算法反演熱源位置的導(dǎo)熱反問題流程框圖.

圖1 PSO算法求解導(dǎo)熱反問題流程圖Fig.1 Flow chart of solving IHCP by PSO

具體實現(xiàn)步驟如下：

Step1 初始化粒子群，在搜索空間內(nèi)隨機給定粒子的初始位置（xi，yi）和初始速度vi；

Step2 根據(jù)粒子的位置代入導(dǎo)熱正問題，求出測點位置的計算溫度值，代入目標(biāo)函數(shù)求其適應(yīng)值；

Step3 根據(jù)各粒子的適應(yīng)值，更新粒子的個體極值pbi、和全局極值gb；

Step4 根據(jù)式（5）-（7）更新各粒子的速度和位置；

Step5 判斷g是否滿足設(shè)定的要求，滿足則計算結(jié)束；不滿足進行下一步；

Step6 判斷迭代次數(shù)是否達到設(shè)定的最大值，若是，則迭代結(jié)束，否則，轉(zhuǎn)到step2.

3 計算結(jié)果及分析

采用上述方法，對含有內(nèi)熱源二維穩(wěn)態(tài)導(dǎo)熱反問題進行了反演計算.反問題物理模型為一個截面形狀為正方形的二維導(dǎo)熱物體，邊長為0.1 m，截面上有一個強度為Φ的點熱源，熱源位置未知，如圖2所示.

圖2 二維帶點熱源導(dǎo)熱反問題示意圖Fig.2 Sketch map of 2-D IHCP with point heat source

熱源位置反演模型的數(shù)學(xué)描述為式（1），邊界條件為

其中，外界流體溫度tf=24℃，對流換熱系數(shù)h= 13 W／（m2·K），導(dǎo)熱系數(shù)λ=1.5 W／（m·K）.邊界上選取6個測點溫度，并根據(jù)導(dǎo)熱正問題模型計算求得.計算中取點熱源強度Φ=72 W，熱源位置為（0.039 5，0.029 5）.表1給出了對應(yīng)測點的溫度值.

計算中PSO算法參數(shù)設(shè)置為粒子規(guī)模m= 20，速度v=（v1，v2）中v1、v2∈［-0.03，0.03］，zmin=0，zmax=0.1，取c1=1.5、c2=2.5［10］，設(shè)置最大迭代次數(shù)為20次.首先取ω為定值進行計算，ω取值分別為0.0，0.1，0.2，…，1.0；然后再按式（7）取ω進行對比計算.為了避免搜索結(jié)果的偶然性，每個工況計算10次，然后對計算結(jié)果取其平均值.采用數(shù)值模擬方法進行計算，計算結(jié)果如表2所示.

表2 計算結(jié)果Tab.2 Results of calculation

從表2所列數(shù)據(jù)可以看出，在0～1范圍內(nèi)，不論慣性系數(shù)ω取何值，粒子群算法最終幾乎都能找到熱源位置，但對搜索速度有影響.當(dāng)ω=0時，式（6）第一項為零，則速度自身無記憶性，粒子飛行速度僅由自身當(dāng)前位置、個體最優(yōu)位置和群體最好位置決定，計算結(jié)果表明，其平均迭代次數(shù)基本達到最大；當(dāng)ω過大時，粒子原有的速度起較大的作用，也就不斷向新的區(qū)域探索，這時粒子需要更多的迭代才能搜尋到熱源位置.

從表2數(shù)據(jù)可以看出，ω取為定值，其取值范圍在0.3～0.7時搜索速度比較快.表2最后一行數(shù)據(jù)為ω按式（7）計算結(jié)果，式（7）中，取ωmax=0.9，ωmin=0.3，ωmax選擇大于0.7，目的是為了使粒子在搜索前期偏向于全局搜索.計算結(jié)果表明，取隨迭代次數(shù)變化的ω可以加快粒子群算法的收斂速度.

反演過程中，每個粒子飛行到一個位置，就要進行一次導(dǎo)熱正問題的計算，以計算目標(biāo)函數(shù).圖3所示為反演過程中目標(biāo)函數(shù)隨正問題計算次數(shù)tm變化曲線.可以看出，導(dǎo)熱計算區(qū)域網(wǎng)格劃分為100× 100時，正問題計算次數(shù)不足200次就能找到導(dǎo)熱區(qū)域中的熱源位置.圖4所示為初始時刻粒子群在求解區(qū)域內(nèi)的位置分布，圖5所示分別為第二次、第四次、第六次和第八次迭代時刻粒子群的位置分布.可以看出，粒子群中的每個粒子每次迭代后根據(jù)式（6）調(diào)整飛行速度的方向和大小，從而使粒子群從初始的隨機分布逐漸趨于集中，最終逼近熱源位置，目標(biāo)函數(shù)達到最小.

圖3 目標(biāo)函數(shù)隨正問題計算次數(shù)的變化Fig.3 Objective function variation with time

圖4 初始粒子群的分布Fig.4 Distribution of particles at the beginning

圖5 迭代過程中粒子的分布Fig.5 Distribution of particles during the process of iterations

4 結(jié) 論

采用粒子群優(yōu)化算法，建立了適合于尋源導(dǎo)熱反問題的求解方法，并對粒子群算法中慣性系數(shù)ω的取值范圍進行討論.ω設(shè)為定值時，其取值范圍在0.3～0.7能取得較好的求解結(jié)果；ω采用式（7）計算時，設(shè)ωmax=0.9、ωmin=0.3，可以在滿足全局最優(yōu)的同時，加快收斂的速度.數(shù)值計算結(jié)果表明，在參數(shù)設(shè)置合適時，利用粒子群優(yōu)化算法反演導(dǎo)熱區(qū)域的熱源位置具有較高的精度和較好的穩(wěn)定性.

［1］ Yang CY.The determination of two moving heat sources in two-dimensional inverse heat problem［J］.Applied Mathematical Modelling，2005，30（3）：278-292.

［2］ Beck J V，Blackwell B，Haji-Sheikh A.Comparison of inverse heat conduction methods using experimental data［J］.International Journal of Heat Mass Transfer，1996，39（17）：3649-3657.

［3］ Huang C H，Chao B H.An inverse geometry problem in identifying irregular boundary configurations［J］. International Journal of Heat and Mass Transfer， 1997，40：2045-2053.

［4］ Refahi Abou khachfe，Yvon Jarny.Determination of heat sources and heat transfer coefficient for twodimensional heat floe-numerical and experimental study［J］.International Journal of Heat and Mass Transfer，2001，4：1309-1322.

［5］ Huang C H，Ozisik M N.Inverse problem of determining the unknown strength of an internal plate heat source［J］.Franklin Inst，1992，329：751-764.

［6］ Silva Neto A J，Ozisik M N.Inverse problem of simultaneously estimating the time-varying strengths of two plane heat sources［J］.Appl Phys，1993，73（5）：2132-2137.

［7］ Huang C H.Transient inverse two-dimensional geometry problem in estimating time-dependent irregular boundary configurations［J］.International Journal of Heat and Mass Transfer，1998，41（12）：1707 -1718.

［8］ 楊海天，薛齊文.兩級敏度分析求解非線性穩(wěn)態(tài)多宗量熱傳導(dǎo)問題［J］.工程熱物理學(xué)報，2003，24（3）：463 -465.

［9］ Kennedy J，Eberhart R.Particle Swarm Optimization［C］／／Proceedings of IEEE International Conference on Neural Networks.Perth：IEEE，1995.

［10］ 范培蕾，張曉今，楊濤.克服早熟收斂現(xiàn)象的粒子群優(yōu)化算法［J］.計算機應(yīng)用，2009，29（6）：122-124.

（編輯：金 虹）

Seeking Heat Source in Inverse Heat Conduction Problem by Using Particle Swarm Optimization

ZHANGTao， LUMei， TAOLiang， LIBo-han
（School of Energy and Power Engineering，University of Shanghai for Science and Technology，Shanghai 200093，China）

A numerical model was built for the solution of heat conduction problems.In solving inverse problems，a least-square based objective function was minimized by PSO.A new robust solving method to determine the heat source in IHCP.was developed.With a number of numerical examples，the evaluation of inertia coefficient in PSO was studied.The results show that the proposed method is an accurate and efficient method to determine the location of heat source in IHCP.

heat conduction；inverse problem；heat source；particle awarm optimization（PSO）

TK 124

A

1007-6735（2013）04-0377-05

2012-12-17

國家自然科學(xué)基金資助項目（51176126）

張 濤（1987-），男，碩士研究生.研究方向：導(dǎo)熱反問題.E-mail：zhangtaobeyond＠126.com

盧 玫（1958-），女，教授.研究方向：強化傳熱和能量有效利用.E-mail：rose_luu＠usst.edu.cn