張鵠志,郭建青

(長安大學(xué)環(huán)境科學(xué)與工程學(xué)院,陜西 西安 710054)

粒子群優(yōu)化(particle swarm optimization,PSO)算法[1-2]是由Eberhart R C和Kennedy J于1995年提出的一種基于迭代的智能隨機(jī)優(yōu)化算法,具有易于理解和實現(xiàn)、全局搜索能力強(qiáng)[3]和計算結(jié)果精度不易受人為因素影響等優(yōu)點。由于PSO算法具有這些優(yōu)點,其在水利工程規(guī)劃設(shè)計領(lǐng)域中已經(jīng)得到了一定的應(yīng)用。例如,國內(nèi)學(xué)者將PSO算法應(yīng)用于水庫(群)優(yōu)化調(diào)度中,為求解水庫調(diào)度問題提供了一種新的有效解決途徑[4-5];還有的學(xué)者將PSO算法引入大壩安全監(jiān)控領(lǐng)域,并結(jié)合多元回歸統(tǒng)計模型建立基于粒子群算法的混凝土壩變形預(yù)報模型,并通過優(yōu)化迭代計算確定壩體變形統(tǒng)計模型中各回歸系數(shù)[6]。在開發(fā)利用地下水資源過程中,利用不同的方法分析抽水試驗數(shù)據(jù),是確定含水層參數(shù)的主要途徑之一。標(biāo)準(zhǔn)曲線配線法是分析抽水試驗數(shù)據(jù)的常用方法之一,但是,其在應(yīng)用過程中存在著人為隨意性,對計算結(jié)果的精度影響較大。為了彌補這個缺陷,人們已將智能優(yōu)化算法應(yīng)用于分析抽水試驗數(shù)據(jù),確定滿足泰斯假設(shè)條件下的含水層參數(shù)[7]。但是,在利用PSO算法分析越流系統(tǒng)含水層中的抽水試驗數(shù)據(jù)并確定含水層參數(shù)方面,以及PSO算法中的3個控制參數(shù)對算法收斂性和計算結(jié)果精度的影響方面開展的研究工作較少。鑒于此,筆者將PSO算法應(yīng)用于確定第1類越流系統(tǒng)的含水層參數(shù),并通過數(shù)值試驗,初步討論算法控制參數(shù)c1,c2和 ω對算法收斂性的影響。

1 PSO算法

PSO算法假想一群粒子的每個粒子都在搜索空間中以一定速度飛行,并根據(jù)本身以及同伴的飛行經(jīng)驗動態(tài)調(diào)整自身速度,從而形成群體尋優(yōu)的機(jī)制,最終搜索到問題的最優(yōu)解[1-2]。人們用粒子的位置表示待優(yōu)化問題的解,將每個粒子位置的優(yōu)劣程度稱為它的適應(yīng)度,由待優(yōu)化問題的目標(biāo)函數(shù)值來定量地表征。每個粒子由1個速度矢量決定其飛行方向和速率大小。設(shè)在1個d維的目標(biāo)搜索空間中,由m個粒子組成1個群體,其中,在第 t次迭代時粒子i的位置表示為Xi(t)=(xi1(t),xi2(t),…,xid(t)),相應(yīng)的飛行速度表示為 V i(t)=(vi1(t),vi2(t),…,vid(t))。開始運算時,首先隨機(jī)初始化這m個粒子的位置和速度,然后通過迭代尋找最優(yōu)解。在每次迭代過程中,各粒子通過跟蹤2個極值來更新自己的速度和位置:1個是該粒子本身迄今搜索到的最優(yōu)解,稱為個體極值 pbest,表示為P i(t)=(pi1(t),pi2(t),…,pid(t));另1個是整個粒子群迄今找到的最優(yōu)解,稱為全局極值gbest,表示為Pg(t)=(pg1(t),pg2(t),…,pgd(t))。將每個粒子更新個體極值的過程稱為歷史尋優(yōu),將各不同粒子通過相互比較來更新全局極值的過程稱為全局尋優(yōu)。在第t+1次迭代計算時,粒子i根據(jù)下列規(guī)則來更新自己的速度和位置[3]:

式中:ω為慣性權(quán)重系數(shù);c1和c2為學(xué)習(xí)因子,文獻(xiàn)[3]中建議 c1=c2=2.0;rand1和 rand2分別為 2個在[0,1]區(qū)間內(nèi)變化的隨機(jī)數(shù);下標(biāo)k表示第k維。

式(1)右邊第1項為個體粒子飛行過程中的“慣性”部分,反映了粒子的運動“習(xí)慣”,代表粒子有維持自己先前速度的趨勢;第2項為“認(rèn)知”部分,反映了粒子對自身歷史經(jīng)驗的記憶,代表粒子有向自身個體極值逼近的趨勢;第3項為“社會”部分,反映了粒子間協(xié)同合作與知識共享的群體歷史經(jīng)驗,代表粒子有向全局極值逼近的趨勢[3]。另外,粒子每一維的速度vik都應(yīng)被一個最大速度所限制。而為粒子位置設(shè)定閾值,從而保證各粒子在限定的搜索空間內(nèi)飛行,也是有待進(jìn)一步研究的問題。

2 PSO算法在確定第1類越流系統(tǒng)含水層參數(shù)中的應(yīng)用

2.1 基本公式

本文討論的第1類越流系統(tǒng)的模型與文獻(xiàn)[8]中相同,抽水過程中的水頭降深公式為

式中:s為水頭降深,[L];Q為抽水流量,[L3?T-1];T為導(dǎo)水系數(shù),[L2?T-1];r為距抽水主井的距離,[L];1/B為越流補給因子,[L-1]。

在應(yīng)用PSO算法確定含水層參數(shù)時,要求式(4)表達(dá)的目標(biāo)函數(shù)值達(dá)到最小,即

式中:Xi(t)為第t次迭代時粒子i的位置,它以各含水層參數(shù)作為其各維上的坐標(biāo),在第1類越流系統(tǒng)中有3個待求的含水層參數(shù),分別是導(dǎo)水系數(shù)T、彈性釋水系數(shù)S和越流補給因子1/B,分別視為xi1(t),xi2(t),xi3(t),這樣粒子的位置坐標(biāo)就是X=(T,S,1/B);j為抽水試驗過程中水頭降深觀測數(shù)據(jù)的序列號,j=1,2,…,N;soj為實際觀測到的第j個降深值,[L];scj為利用式(3)計算的與 soj相同時刻的水頭降深值,[L]。

式(4)的意義為選取適當(dāng)?shù)牧Ｗ游恢?X=(T,S,1/B),使降深的計算值與其觀測值之間的殘差平方和的均值達(dá)到最小。此時的 T,S,1/B即為所求的最優(yōu)解。

2.2 算例

根據(jù)上述基本公式和PSO算法的流程,利用C語言編寫了相應(yīng)的計算程序,在運行開始時,輸入數(shù)據(jù)為相應(yīng)不同抽水時間的水頭降深值,輸出結(jié)果為含水層參數(shù)T,S,1/B的最優(yōu)值。

表1中給出了抽水流量Q=17.0m3/min時,距抽水主井距離r=12.19 m的觀察孔中觀測到的降深隨時間的變化數(shù)據(jù)[10]。

表1 觀測數(shù)據(jù)

根據(jù)多次試算的經(jīng)驗,采用的粒子個數(shù)為50個;收斂標(biāo)準(zhǔn)取目標(biāo)函數(shù)值小于0.00060842。權(quán)重系數(shù)按式(4)進(jìn)行線性衰減:式中:ωmax為最大慣性權(quán)重系數(shù),取 0.9;ωmin為最小慣性權(quán)重系數(shù),取0.1;t′為當(dāng)前迭代次數(shù);n為最大迭代次數(shù),取3000。

取學(xué)習(xí)因子c1=c2=2.0,經(jīng)過平均700次迭代得到計算結(jié)果,見表2。表2也給出了傳統(tǒng)配線法的計算結(jié)果。

表2 不同方法的計算結(jié)果

2.3 可靠性檢驗及與傳統(tǒng)配線法效果比較

分別利用由PSO算法和配線法求出的含水層參數(shù)來計算降深～時間曲線,并與實際的降深～時間曲線進(jìn)行比較,如圖1所示。由圖1可以看出,PSO算法所求結(jié)果對觀測數(shù)據(jù)的擬合程度與傳統(tǒng)配線法很接近,是可靠的。但配線法具有較大的主觀性,不同人員所得結(jié)果的精度可能不同,而PSO算法則可以穩(wěn)定地得到精度較高的結(jié)果。

圖1 PSO算法與配線法效果比較

3 PSO算法控制參數(shù)的作用和影響

由前文可以看出,有5個參數(shù)影響著PSO算法的運算過程和結(jié)果,分別是:粒子個數(shù)、待求參數(shù)初值的取值范圍、收斂標(biāo)準(zhǔn)、c1和c2以及 ω。擬通過數(shù)值試驗結(jié)果,分析 c1,c2和 ω對算法收斂性的影響。

3.1 學(xué)習(xí)因子

3.1.1 定性分析

以 c1為例,從式(1)和式(2)不難看出,c1rand 1[pik(t)-xik(t)]的作用是:在粒子更新位置的過程中,使粒子以“一定程度”向它的個體極值靠攏。而這個“一定程度”是由 c1rand1確定的。比如xik(t)+1.0rand1[pik(t)-xik(t)](即 c1=1.0)就表示粒子向其個體極值隨機(jī)地移動某一距離;而xik(t)+2.0rand1[pik(t)-xik(t)](即 c1=2.0)表示粒子不但可以飛向個體極值,而且也可以飛過個體極值位置到達(dá)另一側(cè)(但不會超出原來相差的距離),如圖2所示。因此,當(dāng)c1=1.0時,粒子直接以其個體極值為目標(biāo),靠攏的趨勢是最強(qiáng)的;當(dāng) c1=2.0時,新位置在個體極值兩側(cè)的分布幾率是最均衡的。

圖2 粒子飛行示意圖

同理,c2rand2確定了粒子向全局極值移動的程度:當(dāng)c2=1.0時,粒子向全局極值靠攏的趨勢最強(qiáng);當(dāng)c2=2.0時,粒子趨于在全局極值周圍幾率均衡地分布。

3.1.2 數(shù)值試驗結(jié)果分析

文中主要對最具代表性的 ci=1.0和 ci=2.0(i=1,2)情況進(jìn)行試算。另由于不同目標(biāo)函數(shù)值對應(yīng)著不同的粒子位置,因此用目標(biāo)函數(shù)值的變化來反映粒子位置的變化情況。

筆者在其他參數(shù)取值相同(c2=1.0,ωmax=0.8,ωmin=0.6,粒子個數(shù)取 5個,收斂標(biāo)準(zhǔn)取0.00060842)的條件下,含水層參數(shù)取值范圍為0.93 m2/min≤T≤3.72m2/min,0.002≤S≤0.005,0.0016 m-1≤1/B ≤0.0033 m-1時,對 c1=2.0和c1=1.0分別進(jìn)行3次試算,輸出每個粒子目標(biāo)函數(shù)值隨迭代次數(shù)增加的變化(篇幅所限未列出具體結(jié)果)。通過試算可以觀察到2種現(xiàn)象:①c1影響粒子的活躍程度,進(jìn)而影響粒子的搜索能力:當(dāng) c1=2.0時,目標(biāo)函數(shù)值在迭代次數(shù)達(dá)到2000次以上時仍然有所變化,即粒子仍然移動,保持著搜索能力;當(dāng)c1=1.0時,目標(biāo)函數(shù)值在迭代次數(shù)達(dá)到1000次以上時就已經(jīng)不再變化,即粒子基本失去搜索能力。②當(dāng)c1=2.0時,3次試算目標(biāo)函數(shù)值分別穩(wěn)定于0.00061086,0.00061082和0.00061098左右;c1=1.0時,3次試算目標(biāo)函數(shù)值分別穩(wěn)定于0.00086816,0.00071975和0.00101040,均大于c1=2.0時的結(jié)果。顯然 c1=2.0時的最終運算結(jié)果具有更高的精度。

產(chǎn)生上述現(xiàn)象的原因是,當(dāng) c1=1.0時,粒子在每次迭代中都強(qiáng)烈地趨向于它的某個個體極值,因此,在經(jīng)過較少次的迭代后粒子就會到達(dá)某個個體極值附近。這時,負(fù)責(zé)歷史尋優(yōu)的c2rand1[pik(t)-xik(t)]項就會趨于零(因為這時pik(t)-xik(t)趨于零),從而使歷史尋優(yōu)機(jī)制形同虛設(shè),粒子便停滯在這個個體極值而不再隨迭代而移動(實際上,試算中c2=1.0的取法加劇了這種作用,它使各粒子的個體極值迅速趨向于某一局部極值)。粒子這種過早的停滯使其不能再搜索到更優(yōu)的值,因此其最終結(jié)果的精度普遍較差。當(dāng)c1=2.0時,粒子在每次迭代中都能幾率均衡地在個體極值兩側(cè)進(jìn)行歷史尋優(yōu),所以不易停滯,從而保持了更強(qiáng)的歷史尋優(yōu)能力,可以找到更優(yōu)的值。

對于c2的取值,在其他參數(shù)取值相同(c1=1.0,ωmax=0.9,ωmin=0.7,粒子個數(shù)取 5 個,收斂標(biāo)準(zhǔn)取0.00060842)的條件下,含水層參數(shù)取值范圍為:0.93 m2/min≤T≤3.72m2/min,0.002≤S≤0.005,0.0016m-1≤1/B ≤0.0033m-1時 ,對 c2=2.0和c2=1.0分別進(jìn)行3次運算,其觀察到的現(xiàn)象以及原因分析都與 c1時非常相似,不同的是c2直接對全局尋優(yōu)而不是歷史尋優(yōu)起作用。

實際上,c1和c2的作用是相互交織的,歷史尋優(yōu)的效果和全局尋優(yōu)的效果是相互影響的。

3.2 慣性權(quán)重系數(shù)

由式(1)可以看出,ω決定了對前一速度值繼承的程度。它影響粒子的搜索能力并影響搜索的精細(xì)程度。

筆者在其他參數(shù)取值相同(c1=c2=1.0,粒子個數(shù)取5個,收斂標(biāo)準(zhǔn)取0.00060842)的條件下,含水層參數(shù)取值范圍為 0.93 m2/min≤T≤3.72m2/min,0.002≤S≤0.005,0.0016 m-1≤1/B≤0.0033 m-1時,對 0.5≤ω≤1.0和0≤ω≤0.5分別進(jìn)行3次試算,輸出每個粒子目標(biāo)函數(shù)值隨迭代次數(shù)增加的變化。通過試算觀察到這樣的現(xiàn)象:當(dāng) ω較大(0.5≤ω≤1.0)、迭代次數(shù)達(dá)到1000次以上時各粒子的目標(biāo)函數(shù)值仍有較大變化,即粒子仍保持搜索能力,3次運算的最終目標(biāo)函數(shù)值分別穩(wěn)定于0.00062041,0.00061140,0.00061102;當(dāng) ω較小(0≤ω≤0.5)、迭代次數(shù)不到1000次時各粒子已經(jīng)基本停滯,3次運算的最終目標(biāo)函數(shù)值分別穩(wěn)定于0.00110930,0.00066110,0.00064549。產(chǎn)生上述現(xiàn)象的原因是,ω值越大,前一速度對當(dāng)前的速度貢獻(xiàn)越大,當(dāng)粒子聚集到個體極值和全局極值附近時,式(1)中c1rand1[pik(t)-xik(t)]+c2rand2[pgk(t)-xik(t)]變得很小,這樣 ωvik(t)項就成為新速度的主要部分;如果 ω很小,那么新的速度就會迅速變小,導(dǎo)致粒子活動范圍越來越小,全局搜索能力降低。但有研究認(rèn)為,ω較小時可增強(qiáng)局部搜索能力[3]。

4 結(jié) 語

文中將PSO算法應(yīng)用于分析第1類越流系統(tǒng)含水層中的抽水試驗數(shù)據(jù),并確定含水層的參數(shù)。通過對數(shù)值試驗結(jié)果的初步分析討論,可以得出:①PSO算法原理簡單,易于編程運算,計算結(jié)果精度不受人為影響;②PSO算法的計算結(jié)果精度高,可靠性強(qiáng);③學(xué)習(xí)因子c1,c2和慣性權(quán)重系數(shù) ω綜合影響著粒子的活躍程度,進(jìn)而影響粒子群的搜索能力;④建議算法控制參數(shù)取值為 c1=c2=2.0,ω由0.9隨迭代次數(shù)線性遞減至0.4?？傊?PSO算法不失為一種分析抽水試驗資料、確定第1類越流系統(tǒng)含水層參數(shù)的有效方法。

[1]KENNEDY J,EBERHART R C.Particle swarm optimization[C]//Proceedings of IEEE International Conference on Neural Networks.Piscataway:IEBE Service Center,1995(4):1942-1948.

[2]DOROGO M,MANIEZZO V,COLORI A.Ant system:optimization by a colony of cooperating agents[J].IEEE Trans on System,Man,and Cybernetics,1996,26(1):28-41.

[3]范娜,云慶夏.粒子群優(yōu)化算法及其應(yīng)用[J].信息技術(shù),2006(1):53-56.

[4]鄧顯羽,彭勇,葉碎高,等.粒子群算法在水庫(群)優(yōu)化調(diào)度研究中的應(yīng)用綜述[J].水利水電科技進(jìn)展,2010,30(5):90-94.

[5]向波,紀(jì)昌明,羅慶松.免疫粒子群算法及其在水庫優(yōu)化調(diào)度中的應(yīng)用[J].河海大學(xué)學(xué)報:自然科學(xué)版,2008,36(2):198-202.

[6]王偉,王連慶.基于粒子群仿生算法的混凝土壩變形預(yù)報模型[J].水利水電科技進(jìn)展,2008,28(4):11-14.

[7]張娟娟.智能優(yōu)化算法在含水層參數(shù)反演中的應(yīng)用[D].西安:長安大學(xué),2006.

[8]陳崇希,林敏.地下水動力學(xué)[M].武漢:中國地質(zhì)大學(xué)出版社,1999:70-122.

[9]楊天行,傅澤周,劉金山,等.地下水向井的非穩(wěn)定流動的原理與計算方法[M].北京:地質(zhì)出版社,1980:195-199.

[10]U.S.Department of the Interior(USDI).Ground watermAnual[M].Washington,D.C.:Bureau of Reclamation,1981:121-127.