浦陽江流域梯級水庫防洪調度兩階段啟發式搜索算法研究

陳 佳 吳國強

(諸暨市水利局，浙江 諸暨 311800)

動態規劃(DP)類算法(包括DP[1]、逐次漸進DP[2]、離散微分DP[3]、折疊DP[4]等)是解決多階段序貫決策問題的經典方法，在梯級水庫優化調度中獲得了廣泛應用。應用DP類算法求解梯級水庫優化調度問題須滿足目標函數階段可分和狀態過程無后效性兩個基本條件。對于梯級水庫調度而言，前一條件通常容易滿足，而后一條件受河道水流傳播延遲(水庫的下泄流至其緊鄰下游水庫需要花費一定的時間)影響，嚴格來說不能滿足。當調度決策時段為旬、月，甚至更長時，對于大多數梯級水庫系統而言，河道水流傳播時間遠遠小于調度決策時段長度，使河道水流傳播延遲對調度的影響基本可以忽略不計，此時，DP類算法仍具有適用性。隨著調度決策時段變短，這種影響將逐漸增強，當調度決策時段長度變至與河道水流傳播時間基本相當或者更短時，這種影響將十分顯著，此時，DP類算法已失去適用性。

在梯級水庫防洪調度中，河道水流傳播延遲對調度的影響通常不能忽略不計，梯級系統最優控制方案不僅與系統當前時刻的蓄水狀態有關，而且與當前時刻以前若干時刻的蓄水狀態有關(即“有后效性”問題)[5]。對于梯級水庫防洪調度問題，通常采用智能算法(如遺傳算法[6]、粒子群優化法[7]、水循環算法[8]等)或者逐步優化算法(POA)[9-12]進行求解。智能算法能夠處理非線性、不可微、水流傳播延遲等各種復雜因素，通用性強，且此類算法大多采用隨機搜索機制，理論上可以產生全局解，但在實際應用中往往出現“早熟”收斂和“不穩定解”等問題[13]，應用效果受到影響。POA算法將多階段問題分解為若干兩階段問題進行逐次求解和迭代逼近，雖然能夠處理“有后效性”問題，但受限于“盲目搜索”問題和“維數障礙”[14]，對于梯級水庫防洪調度問題的求解效率不高。

因此，本文以浦陽江流域防洪控制性水庫聯合調度為研究對象，在POA算法的基礎上，設計一種兩階段啟發式搜索算法，以期為解決梯級水庫防洪調度“有后效性”問題和“維數障礙”問題提供一條新的途徑。

1 問題描述

1.1 流域簡介

浦陽江又名“浣江”，系錢塘江支流，發源于浙江省金華市浦江縣，干流自南向北流經浦江縣、諸暨市、蕭山區，沿途接納大陳江、開化江、五泄江、楓橋江、凰桐江、永興河等主要支流，在蕭山區境內注入錢塘江[15]。浦陽江主流長150km，流域面積3455km2，流域地處亞熱帶季風性氣候區，受地形、天文、氣象等綜合因素影響，流域飽受洪、澇、臺等災害困擾，是浙江省防汛防臺重點流域之一。浦陽江干、支流大部分河段已修筑堤防，但汛期發生洪水時，水面高程往往大于兩岸地面高程，易發生水災，素有“浙江小黃河”之稱，夏季的梅雨和臺風雨是造成流域洪澇災害的主要因素。

1949年后，按照“上蓄、中分、下泄”的治理方針，浦陽江上游興建了通濟橋、安華、陳蔡、石壁等控制性水庫，中游興建了高湖蓄滯洪區，下游拓浚河道、裁彎取直和修建堤防，形成了以通濟橋、安華、陳蔡、石壁等水庫和高湖蓄滯洪區為骨干控制性工程，輔以中小型水庫、干支流堤防、電排泵站等的流域防洪工程體系。流域防洪控制性水庫拓撲關系和特征參數見圖1和表1。

圖1 浦陽江流域防洪控制性水庫拓撲關系

表1 浦陽江流域防洪控制性水庫特征參數 單位：億m3

1.2 目標函數

本文以通濟橋、安華、陳蔡、石壁4座水庫為對象，建立浦陽江流域梯級水庫聯合防洪調度模型。浦陽江流域防洪調度主要目的是在保證各水庫運行安全的基礎上，利用水庫防洪庫容進行攔洪錯峰，盡可能地削減下游重要防洪控制斷面的洪峰流量。為了便于公式化描述，將通濟橋、安華、陳蔡、石壁4座水庫依次編號為1、2、3、4。目標函數表示為

(1)

(2)

式中：Rt為t時段通過諸暨太平橋控制斷面的流量，m3/s；T為調度期總時段數；max(·)為取括號中集合的最大元素；min{·}為最小化括號中的項；ni為i水庫與下游諸暨太平橋控制斷面之間河道洪水演進水流流達時間，s；Qi,t-ni為i水庫t-ni時段的下泄流量，m3/s。

1.3 約束條件

浦陽江流域梯級水庫防洪調度應滿足以下約束條件：

a.水庫水量平衡方程：

Vi,t=Vi,t-1+(Ii,t-Qi,t)Δt

(3)

b.梯級水庫水力聯系：

I2,t=At+Q1,t-n

(4)

c.水庫蓄水量上、下限：

(5)

d.水庫泄流能力限制：

(6)

e.下游河道安全泄量限制：

(7)

f.水庫下泄流量變幅限制：

(8)

g.邊界條件：

(9)

h.流域洪水調度原則要求：各水庫控制運用應遵循浦陽江干流洪水調度原則，詳見文獻[16]。

2 問題求解

2.1 啟發式搜索策略

POA采用兩階段優化策略在問題的解空間內進行迭代搜索，因無須求解貝爾曼方程，且能夠適應河道水流傳播延遲引起的“有后效性”問題，因此在梯級水庫短期調度中獲得了廣泛應用。然而，采用POA求解目標函數式(1)所示問題時發現，當且僅當所求解的兩階段問題涵蓋洪峰流量所在時段時，才能達到削減洪峰的目的。但是，這樣的兩階段問題在POA所有兩階段問題中的占比很小，僅為2/(T-1)(考慮洪峰流量為單值的情況)。換句話說，POA大部分兩階段問題的求解達不到削減洪峰的目的。因此，POA求解目標函數式(1)所示問題的效率并不高。

為了提高POA的求解效率，對POA兩階段優化策略進行改進。研究發現，目標函數式(1)的預期效果是洪水總量在不同時段之間分配盡可能均勻、洪峰流量盡可能小。因此，在峰、谷時段之間進行水量優化分配，可以達到快速改善目標函數式(1)的效果。因此，固定其他T-2個時段的流量值，僅對峰、谷兩時段的流量值進行優化。不妨假定峰、谷流量對應時段號分別為p和q，由水量平衡原理可知，p、q兩時段的流量之和應當為定值。因此，兩階段啟發式搜索算法(以下簡稱“改進算法”)的優化問題可以表示為

s.t.Rp+Rq=C

(10)

式中：Rp、Rq分別為p、q時段(峰、谷流量時段)諸暨太平橋控制斷面的流量，m3/s；C為常數。

本文的改進算法正是通過反復構建并求解式(10)所示兩階段問題，而達到實現洪峰流量最小化的目的的。

2.2 兩階段問題求解

對于不同形式的POA兩階段問題，因無通用的方法可供使用，通常采用窮舉法進行求解。即先將連續的水庫狀態區間離散為由若干狀態值構成的有限集合，然后對所有可能的離散狀態組合方案進行比選，獲得最優狀態組合方案。采用窮舉法求解POA兩階段問題，最大的障礙為所需計算量和存儲量會隨水庫數量呈指數增長，即“維數障礙”。因此，采用增量法對式(10)所示問題進行求解，主要步驟如下：

a.確定各水庫待優化的變量。將式(2)代入式(10)可得，安華、陳蔡、石壁3座水庫待優化的變量為Qi,p-ni和Qi,q-ni(i=2,3,4)；通濟橋水庫待優化的變量為Q1,p-n2-n和Q1,q-n2-n。

b.固定其他變量的值，對各水庫待優化的變量的值進行優化。根據水量平衡原理可知，i(i=1,2,3,4)水庫2個待優化的變量之和為定值，因此，只需對1個量進行尋優，不妨假定為Qi,m,在Qi,m上、下各變動一個增量ΔQi(ΔQi取較小的正數，且對于不同的水庫，ΔQi取值可以不相同)，形成3個流量值，即Qi,m+ΔQi、Qi,m和Qi,m-ΔQi。對4座水庫81種流量組合方案進行比選，獲得max(Rp，Rq)值最小的方案，并按此方案更新4座水庫待優化的變量的值和Rp、Rq的值。

c.反復執行b步驟，直到|Rp-Rq|滿足指定的精度要求為止。

2.3 主要步驟和流程圖

本文的改進算法的主要步驟如下：

a.產生一條可行的水庫系統初始下泄軌線，并計算諸暨太平橋控制斷面的流量過程{Rt|t=1,2,…,T}。水庫系統初始下泄軌線通常可以自上而下依次對各水庫采用等流量調節方法產生。

b.構建式(10)所示兩階段優化問題。值得注意的是，若式(10)中的Rq已無法上調(到達約束邊界上)，則應當用{Rt|t=1,2,…,T}序列中第2小流量對應時段替換q，若替換后依然無法上調，再用序列中第3小流量對應時段進行替換，依次類推，直到某個時段的流量可以上調為止。

c.采用增量法對b步驟構建的兩階段優化問題進行求解。

d.反復執行b～c步驟，直到達到最大迭代次數或者相鄰兩代最優下泄軌線的偏差滿足指定的精度要求為止。

將上述步驟在計算機中進行編程實現。程序運行之前，需要對相關參數和變量進行初始化。參數初始化是指給搜索步長ΔQi、最大迭代次數Num和收斂精度Lim賦初值；變量初始化是指為水庫系統指定一條可行的初始下泄軌線。初始化完成之后，算法將按圖2所示流程運行。

圖2 兩階段啟發式搜索算法流程

3 案例研究

將本文的改進算法應用于浦陽江流域2021年“煙花”臺風洪水優化調度中。“煙花”臺風期間，浦陽江流域自7月22日開始出現降雨，至7月29日中午降雨基本結束，降雨主要集中在23—25日，流域最大3日降雨量為258.0mm，最大5日降雨量為294.6mm。該案例調度期取7月23日9時至28日7時，以1h為時段長度，總計118個時段，洪水總量為2.40億m3。在不考慮水庫調蓄的情況下，諸暨太平橋控制斷面天然洪水洪峰流量為1330.50m3/s，峰現時間為26日6時，峰谷差為1302.88m3/s。為了更好地展示本文的改進算法的應用效果，將其優化結果與POA的優化結果進行對比分析。各水庫調度期初蓄水量取實測值，調度期末蓄水量取各自汛限水位對應的蓄水值。POA兩階段問題采用窮舉法進行求解。為簡單起見，各水庫搜索步長ΔQi統一取10m3/s，且與POA的離散精度相同。兩種算法優化后，諸暨太平橋控制斷面洪水過程對比見圖3。

圖3 不同算法優化后諸暨太平橋控制斷面洪水過程對比

從峰值看，本文的改進算法與POA優化后，諸暨太平橋控制斷面洪峰流量分別為593.47m3/s和923.47m3/s(峰現時間均為25日23時)，與天然洪水相比，削峰率分別為55.39%和30.59%，“削峰”效果顯著，但前者效果更好，更大程度地保障了水庫系統下游城區的防洪安全。

從洪水過程看，本文的改進算法與POA優化后，諸暨太平橋控制斷面洪水過程峰谷差分別為369.58m3/s和889.37m3/s，與天然洪水過程相比，分別下降了71.63%和31.74%，“削峰填谷”效果明顯，但前者效果更好，洪水過程更平穩、波動更小，方案實用性更強。

從計算耗時看，本文的改進算法耗時為3.4s，遠遠小于POA(POA耗時為14808.8s)，說明與后者相比，本文的改進算法效率更高(兩種算法在同一臺PC上進行測試，PC配置Intel Core i5-7500 3.40GHz CPU和8GB RAM)。研究發現，POA耗時長主要有兩個方面的原因：?采用窮舉法求解兩階段問題，計算工作量隨水庫數量呈指數增長，存在“維數障礙”；?大部分兩階段問題不涵蓋洪峰流量所在時段，對其求解達不到削減洪峰的目的。

4 結 語

為了更好地解決梯級水庫防洪調度“有后效性”問題和“維數障礙”問題，本文對POA進行改進，設計了一種POA改良版本。本文的改進算法首先采用逐次逼近策略在解空間內進行迭代搜索，解決水庫系統狀態過程有后效性導致難以使用DP類算法進行求解的難題；其次，對POA兩階段優化策略進行改進，采用啟發式搜索策略對峰、谷流量時段進行兩階段優化，提高POA求解效率和解的質量；最后，采用增量法對兩階段問題進行求解，緩解POA的“維數障礙”。根據浦陽江流域防洪控制性水庫聯合調度應用結果，本文得出以下結論：

a.本文的改進算法通過迭代搜索，逼近問題的解，因無須求解貝爾曼方程，能夠適應梯級系統狀態過程有后效性的特點，且通過應用啟發式搜索策略，與POA相比，削峰率提高了35.73%，峰谷差降低了58.44%，說明本文改進算法能夠更好地解決梯級水庫防洪調度“有后效性”問題。

b.本文的改進算法采用啟發式搜索策略構建兩階段問題，使用增量法對兩階段問題進行求解，與POA相比，在不犧牲解的質量的情況下，算法耗時縮短了99.98%，說明本文改進算法能夠更好地解決梯級水庫防洪調度“維數障礙”問題。

此外，本文的改進算法概念簡單，易于實現，其適用范圍不局限于梯級水庫防洪調度，亦可推廣應用于求解其他具有min{max(·)}或者max{min(·)}形式目標函數的水資源系統優化問題(如水電調峰問題)。因本文的關注點在于改進POA的性能，模型中河道洪水演進采用簡易平移法，雖然會引進一定的偏差，但不影響本文的主要結論。下一步研究將致力于完善河道洪水演進模型，進一步提高調度方案的可靠性和實用性。