王成民，王 昊，刁艷芳

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

1.1.2 約束條件

（1）水量平衡方程。

（2）水庫(kù)下泄流量約束。

（3）水庫(kù)泄流能力約束。

（4）庫(kù)容約束。

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

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

1.2 模型的求解

1.2.1 FOA計(jì)算步驟

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

（6）計(jì)算味道濃度差值?s=SmellbG+1-SmellbG。依據(jù)Metropolis準(zhǔn)則，當(dāng)Δs≤0時(shí)，SmellbG+1被接受，保存該果蠅個(gè)體；當(dāng)Δs>0時(shí)，依據(jù)式（11）的概率pG=exp(-Δs/TG)收斂。若此時(shí)的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 日照水庫(kù)簡(jiǎn)介

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

2.2 結(jié)果對(duì)比

收集了日照水庫(kù)20200813和20200826兩場(chǎng)洪水資料，包括水庫(kù)的水位、庫(kù)容、入庫(kù)流量和泄流量，時(shí)段長(zhǎng)為2 h，同時(shí)收集了水庫(kù)的水位～庫(kù)容～泄量關(guān)系等數(shù)據(jù)，滿足建立水庫(kù)防洪優(yōu)化調(diào)度模型的要求。參考已有文獻(xiàn)［11，17-19］及對(duì)兩場(chǎng)洪水的計(jì)算，確定SA-FOA的初始參數(shù)為Gmax=2 000，L0=2，size=100，T1=106，α=0.9。采用SA-FOA對(duì)日照水庫(kù)雙目標(biāo)防洪優(yōu)化調(diào)度模型進(jìn)行求解，得到非劣解集中各調(diào)度方案的目標(biāo)函數(shù)見表1和表2，同時(shí)采用FOA進(jìn)行求解，兩個(gè)算法求得的調(diào)度方案非劣解前端見圖1和圖2，兩場(chǎng)洪水部分方案的泄流過程見圖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)系，這是因?yàn)閮烧叽淼哪繕?biāo)存在沖突導(dǎo)致的。由圖1和圖2可見，在同一最高水位時(shí)，SA-FOA求得的泄流量平方和小于FOA，在同一泄流量平方和時(shí)，SA-FOA求得的最高水位小于FOA，由此可見，SA-FOA能夠比FOA求得更優(yōu)的目標(biāo)值。

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

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

3 結(jié) 論

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

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