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

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

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

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

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

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

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

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

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

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

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

2.1 反演模型

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

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

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

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

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

PSO的粒子速度更新方程為

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

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

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

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

2.3 尋源反演過(guò)程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

計(jì)算中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次.首先取ω為定值進(jìn)行計(jì)算，ω取值分別為0.0，0.1，0.2，…，1.0；然后再按式（7）取ω進(jìn)行對(duì)比計(jì)算.為了避免搜索結(jié)果的偶然性，每個(gè)工況計(jì)算10次，然后對(duì)計(jì)算結(jié)果取其平均值.采用數(shù)值模擬方法進(jìn)行計(jì)算，計(jì)算結(jié)果如表2所示.

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

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

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

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

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

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

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

4 結(jié) 論

采用粒子群優(yōu)化算法，建立了適合于尋源導(dǎo)熱反問(wèn)題的求解方法，并對(duì)粒子群算法中慣性系數(shù)ω的取值范圍進(jìn)行討論.ω設(shè)為定值時(shí)，其取值范圍在0.3～0.7能取得較好的求解結(jié)果；ω采用式（7）計(jì)算時(shí)，設(shè)ωmax=0.9、ωmin=0.3，可以在滿足全局最優(yōu)的同時(shí)，加快收斂的速度.數(shù)值計(jì)算結(jié)果表明，在參數(shù)設(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］ 楊海天，薛齊文.兩級(jí)敏度分析求解非線性穩(wěn)態(tài)多宗量熱傳導(dǎo)問(wèn)題［J］.工程熱物理學(xué)報(bào)，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］.計(jì)算機(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

國(guó)家自然科學(xué)基金資助項(xiàng)目（51176126）

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

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