基于SA-FOA的水庫防洪優(yōu)化調(diào)度研究

王成民，王 昊，刁艷芳

（1.日照水庫管理運(yùn)行中心，山東 日照 276800； 2.山東農(nóng)業(yè)大學(xué)水利土木工程學(xué)院，山東 泰安 271018）

0 引 言

洪水災(zāi)害是目前發(fā)生頻率最高的自然災(zāi)害，據(jù)統(tǒng)計，僅2018年全球洪水災(zāi)害事件共發(fā)生109起，死亡人口和受災(zāi)人口分別為1 995 人和1 262 萬人，直接經(jīng)濟(jì)損失45億美元［1］，其中，中國、印度、印度尼西亞、美國等國家是洪水災(zāi)害發(fā)生次數(shù)和遭受損失最多的國家。雖然在近30年間，全球洪水災(zāi)害死亡人口及受災(zāi)規(guī)模都呈現(xiàn)持續(xù)下降趨勢，但經(jīng)濟(jì)損失卻呈上升趨勢。水庫是在河流上筑壩進(jìn)而起到防洪效果的一種重要水利工程，且隨著城市供水規(guī)模的不斷擴(kuò)大和農(nóng)業(yè)用水灌溉方式的多樣化，水庫在防洪、供水、發(fā)電及生態(tài)環(huán)境保護(hù)等方面發(fā)揮了巨大的社會和經(jīng)濟(jì)效益。科學(xué)的水庫防洪調(diào)度設(shè)計是安全經(jīng)濟(jì)建設(shè)水庫與實(shí)時調(diào)度發(fā)揮設(shè)計效益的關(guān)鍵。目前，常見的兩類水庫調(diào)度方法包括水庫常規(guī)調(diào)度和水庫優(yōu)化調(diào)度，后者由于可以通過設(shè)定優(yōu)化目標(biāo)來制定優(yōu)化作業(yè)方案，因此在水庫防洪調(diào)度中廣泛應(yīng)用。起初，水庫優(yōu)化調(diào)度所采用的方法主要有線性規(guī)劃［2］、非線性規(guī)劃、動態(tài)規(guī)劃［3］、逐次優(yōu)化算法［4］及大系統(tǒng)分解協(xié)調(diào)原理［5］等，然而，面對由水庫群、蓄滯洪區(qū)、湖泊及其他防洪工程組成的復(fù)雜防洪系統(tǒng)，上述優(yōu)化算法表現(xiàn)出收斂效率低、維度低等明顯局限性。

現(xiàn)代計算技術(shù)的進(jìn)步使得應(yīng)用生物學(xué)、物理學(xué)和人工智能原理的啟發(fā)式智能算法得到了廣泛發(fā)展，這些算法可以有效解決上述問題，主要包括遺傳算法［6］、蟻群算法［7］、人工魚群算法［8］及粒子群算法［9］等。由于智能啟發(fā)式算法相比傳統(tǒng)的優(yōu)化算法能夠解決非線性、間斷、不可微和多維的復(fù)雜問題，因此在水庫優(yōu)化調(diào)度中得到了廣泛的應(yīng)用。然而，啟發(fā)式智能算法也存在不足，例如遺傳算法易早熟且計算量大；蟻群算法、人工魚群算法過于復(fù)雜且計算量大；粒子群算法易陷入局部最優(yōu)值。上述智能算法因其各自尋優(yōu)時存在的弊端，使得應(yīng)用受到限制。果蠅優(yōu)化算法（FOA）［10］是近年來提出的一種新型智能優(yōu)化算法，與其他算法相比，它具有原理簡單、參數(shù)少、計算量小及易于實(shí)現(xiàn)等優(yōu)點(diǎn)，故受到了國內(nèi)外眾多學(xué)科領(lǐng)域?qū)W者的廣泛關(guān)注。通過對該算法的大量研究發(fā)現(xiàn)，該算法在遇到高維優(yōu)化問題時，存在過早收斂和易于陷入局部最優(yōu)等問題，影響了算法的穩(wěn)定性［11］。為彌補(bǔ)FOA的不足，采用設(shè)置變動的搜索距離L和模擬退火算法（Simulated Annealing algorithm，SA）的Metropolis準(zhǔn)則［12］對FOA的迭代過程進(jìn)行擾動，有效避免上述缺陷。

基于優(yōu)先考慮水庫及下游防護(hù)區(qū)防洪安全的前提，建立水庫雙目標(biāo)防洪優(yōu)化調(diào)度模型；然后，提出耦合模擬退火算法的果蠅優(yōu)化算法（SA-FOA）對模型進(jìn)行求解；最后，以日照水庫2020年的兩場洪水調(diào)度過程為例，驗證了該算法的可行性和合理性。

1 水庫雙目標(biāo)防洪優(yōu)化調(diào)度模型的建立和求解

1.1 水庫雙目標(biāo)防洪優(yōu)化調(diào)度模型的建立

水庫自身、庫區(qū)和下游防護(hù)區(qū)是水庫防洪調(diào)度的3個主要防護(hù)對象，其中，為保證水庫及庫區(qū)安全，希望水庫汛期盡可能貼近汛限水位運(yùn)行，故常以水庫調(diào)洪最高水位最低或最大防洪庫容最小作為目標(biāo)函數(shù)；為保障下游防護(hù)區(qū)安全，希望水庫能夠最大限度地削減洪峰流量，故常以最大削峰準(zhǔn)則作為目標(biāo)函數(shù)［13］。本文綜合考慮3個主要防護(hù)對象，建立水庫雙目標(biāo)防洪優(yōu)化調(diào)度模型，以實(shí)現(xiàn)調(diào)洪最高水位最低和最大限度削減洪峰流量。

1.1.1 目標(biāo)函數(shù)

（1）水庫及庫區(qū)防洪安全目標(biāo)。調(diào)洪最高水位越小，水庫及庫區(qū)越安全，因此以壩前最高水位作為目標(biāo)函數(shù)，即：

式中：T為調(diào)洪歷時；Z(t)為t時刻水庫壩上水位，m。

（2）下游防護(hù)區(qū)安全目標(biāo)。經(jīng)證明，水庫最大削峰準(zhǔn)則的目標(biāo)是水庫出庫流量過程平方和最小，即下泄流量最均勻［14］。利用該準(zhǔn)則所構(gòu)建的目標(biāo)函數(shù)的表達(dá)式為［14-16］：

式中：q(t)為t時刻水庫的泄流量，m3/s；R（t）為t時刻水庫到下游防護(hù)區(qū)之間的區(qū)間來水流量，如無區(qū)間來水，其取值為0，m3/s。

1.1.2 約束條件

（1）水量平衡方程。

（2）水庫下泄流量約束。

（3）水庫泄流能力約束。

（4）庫容約束。

（5）非負(fù)約束：以上變量均為非負(fù)。

式中：Q(t)為第t時段水庫的入庫流量，m3/s；V（t）為第t時段水庫的庫容，m3；Δt為時段長，h；Φ[V(t)]為水庫蓄水量為V(t)時的最大下泄能力，m3/s；q(t)min、q(t)max分別為第t時段水庫的泄流量下限、上限，m3/s；V(t)min、V(t)max分別為第t時段水庫的庫容下限、上限，m3。

1.2 模型的求解

1.2.1 FOA計算步驟

FOA由潘文超于2011年提出［17］，是一種通過模擬果蠅種群覓食行為尋求最優(yōu)解的方法，其尋優(yōu)過程如下：

（1）初始參數(shù)的設(shè)置，包括果蠅群里初始位置X0和Y0，最大覓食代數(shù)Gmax，搜索距離L，種群規(guī)模size。

（2）果蠅個體通過嗅覺隨機(jī)搜尋食物，其位置為：

式中：i=1，2，…，size。

（3）計算果蠅個體與原點(diǎn)之間的距離Disti，再計算味道濃度判斷值Si，該值可取為距離的倒數(shù)，即：

（4）將Si代入味道濃度判斷函數(shù)中，求解果蠅個體的味道濃度Smelli：

（5）找出果蠅種群中味道濃度最佳值Smellb，如果目標(biāo)函數(shù)是求解最小值，則味道濃度最小值的個體為最佳，反之亦然。

（6）保留最佳味道濃度Smellb及其位置坐標(biāo)，此時果蠅種群利用視覺向該位置飛行。

（7）迭代尋優(yōu)，重復(fù)步驟（2）～（4），判斷當(dāng)前味道濃度是否優(yōu)于前一次迭代的最佳味道濃度，并判斷當(dāng)前迭代次數(shù)是否小于Gmax，若是執(zhí)行步驟（6）；否則，結(jié)束迭代輸出最優(yōu)值。

1.2.2 SA算法

SA算法是Kirkpatrick等人［18］受到固體或熱力學(xué)系統(tǒng)的物理退火法的啟發(fā)而提出的一種概率優(yōu)化方法。在該方法中，熱力系統(tǒng)的當(dāng)前狀態(tài)、能量方程和基態(tài)類似于優(yōu)化問題的當(dāng)前調(diào)度解、目標(biāo)函數(shù)和全局最優(yōu)解。SA算法使用基于概率的Me‐tropolis準(zhǔn)則搜索全局最優(yōu)并跳出局部最優(yōu)，該準(zhǔn)則設(shè)置接受劣解的概率［18，19］，其定義如下：

式中：G為迭代次數(shù)；pG是第G次迭代中的接受概率；EG是第G次迭代中的目標(biāo)函數(shù)值；Eg是歷史最優(yōu)目標(biāo)函數(shù)值；KbO是Boltzmann系數(shù)；TG是第G次迭代中的退火溫度，該溫度最初設(shè)置為較大的值，然后通過以下溫度控制函數(shù)減小到較小的值：

式中：α為退火系數(shù)，其取值區(qū)間為［0.8，1.0］。因此，隨著TG的逐漸減小，劣解的接受概率將接近0。

1.2.3 SA-FOA

FOA中果蠅移動的搜索距離L會影響其尋優(yōu)的結(jié)果，L大有利于跳出局部極值進(jìn)而尋求全局最優(yōu)，并且提高收斂效率，但其局部搜索能力將下降；L過小會導(dǎo)致陷入局部極值，且降低收斂效率。在算法迭代的過程中，對L大小的需求不同。在算法迭代初期，需要較大的L以提高收斂速度；在迭代后期，需要較小的L以避免錯過最優(yōu)值。因此，本文采用隨著迭代次數(shù)逐漸減小的搜索距離，公式如下：

再者引入SA算法，將目標(biāo)函數(shù)增量?s=SmellbG+1-SmellbG作為衡量標(biāo)準(zhǔn)，其中，SmellbG為第G次迭代中果蠅種群中味道濃度最佳值，根據(jù)Metropolis準(zhǔn)則依概率對產(chǎn)生的新個體進(jìn)行選擇，選取進(jìn)入下一代的群體。具體方法見1.2.4節(jié)的第（6）步。該方法可以很好地確定搜索方向，又能依概率舍棄，跳出局部最優(yōu)，使算法能夠較快的搜索到最優(yōu)位置。

1.2.4 SA-FOA確定水庫防洪優(yōu)化調(diào)度方案的步驟

SA-FOA確定水庫防洪優(yōu)化調(diào)度方案的步驟如下：

（1）初始參數(shù)的設(shè)置，包括最大覓食代數(shù)Gmax，初始搜索距離L0，當(dāng)前覓食代數(shù)G=1，種群規(guī)模size，空間維度dim=T，即洪水時段數(shù)，初始溫度TG，退火系數(shù)α。

（2）隨機(jī)生成初始種群，即水庫初始泄流過程線，即：

（3）將種群代入味道濃度判斷函數(shù)中，求出果蠅個體的味道濃度值，即：

式中：為第i個種群在第G次覓食時的最高水位，m；Zmax和Zmin分別為調(diào)度期水庫允許的最高水位和最低水位，m；為第i個種群在第G次覓食時泄流量的平方和分別為第G次覓食時，所有種群泄流量平方和的最大值和最小值。

（4）尋求果蠅種群中味道濃度最佳的個體。本文是求最小值問題，則果蠅種群中味道濃度最小值的個體為最佳，即SmellbG。

（5）采用式（13）計算LG，果蠅群體的位置變?yōu)椋?/p>

再由第（3）和第（4）步驟求出SmellbG+1。

（6）計算味道濃度差值?s=SmellbG+1-SmellbG。依據(jù)Metropolis準(zhǔn)則，當(dāng)Δs≤0時，SmellbG+1被接受，保存該果蠅個體；當(dāng)Δs>0時，依據(jù)式（11）的概率pG=exp(-Δs/TG)收斂。若此時的pG>rand，則新的最小值被接受，將(t)替代作為下次尋優(yōu)的初始位置；否則，保持不變，開始進(jìn)行下一次尋優(yōu)。

（7）令G=G+1，如果G≤Gmax，轉(zhuǎn)入第（3）步繼續(xù)進(jìn)行迭代，否則停止迭代，輸出最優(yōu)解。

2 實(shí)例研究

2.1 日照水庫簡介

日照水庫位于山東省日照市東港區(qū)，控制流域面積548 km2，總庫容3.18億m3，興利庫容1.82 億m3，是一座集防洪、供水、發(fā)電、灌溉等多種功能于一體的大（II）型水庫，同時也是日照市區(qū)重要飲用水水源地。日照水庫保護(hù)著下游后村鎮(zhèn)、高興鎮(zhèn)、南湖鎮(zhèn)等7個鄉(xiāng)鎮(zhèn)、街道辦事處90個自然村，9.06萬人民群眾的生命財產(chǎn)安全，及204國道、兗石鐵路、日蘭和沈海高速公路和青連高鐵等交通干線以及臨港大中型企業(yè)等防洪重點(diǎn)保護(hù)對象。日照水庫防洪標(biāo)準(zhǔn)為100年一遇設(shè)計、5 000年一遇校核，汛限水位42.50 m、興利水位43.00 m、防洪高水位43.50 m、設(shè)計洪水位43.97 m、校核洪水位46.52 m。汛期洪水調(diào)度原則為，當(dāng)雨前或預(yù)計雨后水位低于42.50 m時，不泄洪；當(dāng)雨前或預(yù)計雨后水位42.50～43.50 m時，控泄流量1 000 m3/s；當(dāng)雨前或預(yù)計雨后水位高于43.50 m時，溢洪閘自由敞泄；當(dāng)水位高于43.97 m時，北放水洞參與泄洪。

2.2 結(jié)果對比

收集了日照水庫20200813和20200826兩場洪水資料，包括水庫的水位、庫容、入庫流量和泄流量，時段長為2 h，同時收集了水庫的水位～庫容～泄量關(guān)系等數(shù)據(jù)，滿足建立水庫防洪優(yōu)化調(diào)度模型的要求。參考已有文獻(xiàn)［11，17-19］及對兩場洪水的計算，確定SA-FOA的初始參數(shù)為Gmax=2 000，L0=2，size=100，T1=106，α=0.9。采用SA-FOA對日照水庫雙目標(biāo)防洪優(yōu)化調(diào)度模型進(jìn)行求解，得到非劣解集中各調(diào)度方案的目標(biāo)函數(shù)見表1和表2，同時采用FOA進(jìn)行求解，兩個算法求得的調(diào)度方案非劣解前端見圖1和圖2，兩場洪水部分方案的泄流過程見圖3和圖4。

圖1 20200813洪水調(diào)度非劣方案集Fig.1 Non-inferior scheme set of 20200813 flood

圖2 20200826洪水調(diào)度非劣方案集Fig.2 Non-inferior scheme set of 20200826 flood

圖3 20200813洪水部分調(diào)度方案下泄流量過程Fig.3 Discharges of partial operation schemes of 20200813 flood

圖4 20200826洪水部分調(diào)度方案下泄流量過程Fig.4 Discharges of partial operation schemes of 20200826 flood

表1 20200813洪水調(diào)度非劣方案集Tab.1 Non-inferior scheme set of 20200813 flood

表2 20200826洪水調(diào)度非劣方案集Tab.2 Non-inferior scheme set of 20200826 flood

（1）由表1、2和圖1、2可見，由SA-FOA和FOA求得的調(diào)度方案集的最高水位和泄流量平方和近似成反比關(guān)系，這是因為兩者代表的目標(biāo)存在沖突導(dǎo)致的。由圖1和圖2可見，在同一最高水位時，SA-FOA求得的泄流量平方和小于FOA，在同一泄流量平方和時，SA-FOA求得的最高水位小于FOA，由此可見，SA-FOA能夠比FOA求得更優(yōu)的目標(biāo)值。

（2）20200821洪水的實(shí)測最高水位和泄流量平方和分別為42.74 m和3 477 174，由表1可見，均大于SA-FOA求得方案集的目標(biāo)值，比方案集中最高水位的最大值和最小值分別大0.17%和0.25%，比方案集中的泄流量平方和的最大值和最小值分別大14.87%和22.35%。20200826洪水的實(shí)測最高水位和泄流量平方和分別為42.90 m和1 133 554，由表2可見，同樣均大于SA-FOA求得方案集的目標(biāo)值，比方案集中最高水位的最大值和最小值分別大0.27%和0.35%，比方案集中的泄流量平方和的最大值和最小值分別大5.63%和71.54%。綜上，SAFOA求得的泄流方案優(yōu)于實(shí)測值。

（3）由圖3可見，20200821洪水的4個優(yōu)化方案的最大泄流量均小于實(shí)測值，且比實(shí)測泄流過程提前加大了泄流，這是導(dǎo)致優(yōu)化方案最高水位低于實(shí)測值的原因。由圖4可見，20200826洪水的4個優(yōu)化方案的最大泄流量均大于實(shí)測值，然而其不超過日照水庫的控泄流量1 000 m3/s，故對下游防護(hù)區(qū)也是安全的。

3 結(jié) 論

（1）鑒于FOA存在過早收斂、易于陷入局部最優(yōu)的不足，本文提出了SA和FOA耦合算法SA-FOA，該算法采用SA的Metropolis準(zhǔn)則引入FOA，使其在迭代過程中依概率接受最優(yōu)解，同時采用逐漸遞減的搜索距離L，并將該算法求解水庫防洪優(yōu)化調(diào)度模型。

（2）由日照水庫兩場洪水調(diào)度的實(shí)例可以看出，SA-FOA求解的調(diào)度方案集的最高水位和泄流量平方和均小于實(shí)測值和FOA求得的調(diào)度方案集，由此表明SA-FOA求得的優(yōu)化調(diào)度方案集既能夠保證水庫和庫區(qū)安全，又能提高下游防護(hù)區(qū)安全。綜上所述，本文建立的水庫防洪優(yōu)化調(diào)度模型及SA-FOA能夠制定出滿意的洪水調(diào)度決策方案，為決策者提供技術(shù)支持。