趙志鵬，廖勝利，程春田，鐘儒鴻，王昱倩

（1.大連理工大學(xué)，遼寧 大連 116024；2.華能瀾滄江水電股份有限公司，云南 昆明 650206）

1 研究背景

截止到2016年底，我國水電裝機容量和年發(fā)電量已經(jīng)突破3億kW和1萬億kW·h，分別占全部能源的20.9%和19.4%［1］，水電已經(jīng)成為中國占比最大的可再生能源［2］。充分利用水電資源，進(jìn)行水庫群優(yōu)化調(diào)度是構(gòu)建清潔低碳，安全高效的現(xiàn)代能源體系的重要措施［1］。梯級水電站群優(yōu)化調(diào)度是一個典型的高維度、多階段、非線性、非凸優(yōu)化問題［3］，其求解難度隨著電站數(shù)量的增加，相互制約約束條件的引入，復(fù)雜電力和水力聯(lián)系的耦合而急劇增加。目前常用的求解方法包括線性規(guī)劃［4］、非線性規(guī)劃［5］、網(wǎng)絡(luò)流［6］、大系統(tǒng)分解協(xié)調(diào)方法［7］、動態(tài)規(guī)劃及其改進(jìn)算法［8-10］，以及以遺傳算法、粒子群算法、差分進(jìn)化算法等為代表的新興智能算法［11-14］等。線性規(guī)劃，大系統(tǒng)分解協(xié)調(diào)，網(wǎng)絡(luò)流等方法都需要對系統(tǒng)進(jìn)行一定的近似處理，容易產(chǎn)生偏差；非線性規(guī)劃對約束條件和目標(biāo)函數(shù)的函數(shù)形式有一定的要求，難以應(yīng)用于實際；智能算法，由于隨機因素的存在，其解穩(wěn)定性較差。動態(tài)規(guī)劃因其對約束條件和目標(biāo)函數(shù)的函數(shù)形式限制較少，結(jié)果穩(wěn)定等特性，是水庫調(diào)度領(lǐng)域應(yīng)用最廣泛的算法之一，但是隨著計算規(guī)模的增加面臨著嚴(yán)重的“維數(shù)災(zāi)”問題。

逐步優(yōu)化算法［15］（Progressive Optimality Algorithm，POA）作為動態(tài)規(guī)劃的改進(jìn)算法，將多階段問題轉(zhuǎn)化為多個兩階段問題，降低了階段間的狀態(tài)組合數(shù)目，被廣泛應(yīng)用于水電優(yōu)化調(diào)度［16-19］。在求解兩階段子問題時，POA采用給定離散步長的網(wǎng)格搜索法，確保了在求解子問題上的全局收斂性，但隨著規(guī)模的增加，仍然存在嚴(yán)重的“維數(shù)災(zāi)”問題。針對此問題，已有文獻(xiàn)將POA算法與正交抽樣［18］、逐次逼近算法（DPSA）［20-21］、并行計算［22］等方法進(jìn)行組合優(yōu)化，使“維數(shù)災(zāi)”問題得到緩解。POA算法中的兩階段子優(yōu)化問題實質(zhì)上是一個多變量單值函數(shù)的最優(yōu)化問題。本文提出離散梯度（Discrete Gradient）的概念，并利用離散梯度下降法對子優(yōu)化問題求解，從而提出離散梯度逐步優(yōu)化算法（DGPOA）。該方法可在不直接求導(dǎo)的情況下，使用離散梯度的信息，確定給定離散步長下的最優(yōu)搜索方向，并利用線搜索迭代方程進(jìn)行迭代求解，可以顯著減少計算規(guī)模，快速獲得優(yōu)化結(jié)果。

2 調(diào)度模型

2.1 目標(biāo)函數(shù)以發(fā)電為主的中長期優(yōu)化調(diào)度一般采用發(fā)電量最大模型。本文計算步長選為1個月，調(diào)度周期為1年。

式中：i、j分別為水庫和時段序號；N為參與計算的水庫總數(shù)；T為總時段數(shù)，本文設(shè)置為12；E為總發(fā)電量（kW·h）； Ai為水庫i的出力系數(shù)；Hi，j為水庫i時段j的平均凈水頭（m）；Δtj為時段j的小時數(shù)（h）；Qi，j為水庫i在時段j的發(fā)電流量（m3/s）。

2.2 約束條件

（1）始末水位約束：

（2）水位上下限：

（3）出力上下限：

（4）水量平衡方程：

其中

式中：Vi，j為水庫i時段j的末庫容（m3）； Ii，j、 qi，j、 Si，j分別為水庫i 時段j的入庫流量、區(qū)間流量以及棄水流量（m3/s），其中qi，j為已知；Ωi為水庫i的直接上游水庫集合，對于龍頭水庫Ωi為空集。

（5）發(fā)電流量限制：

（6）出庫流量限制：

3 逐步優(yōu)化算法

POA算法將多階段問題轉(zhuǎn)化為多個兩階段問題，求解時固定其他時段，僅對求解的兩時段問題進(jìn)行優(yōu)化計算。計算結(jié)果作為其他時段的初始解，如此逐個時段求解，直到收斂。計算步驟如下：（1）確定各水庫在調(diào)度期內(nèi)的初始水位過程，調(diào)度階段總數(shù)設(shè)為T，令E等于初始水位過程發(fā)電量。（2）設(shè)t=T-1。（3）固定其他調(diào)度時段，優(yōu)化各水庫水位Zi，t，使t階段和t+1階段梯級發(fā)電總量E(Z1，t，Z2，t，…，Zn，t) 達(dá)到最優(yōu)。（4）令t=t-1，如果 t≤ 0，則轉(zhuǎn)至步驟5，否則轉(zhuǎn)至步驟3。（5）計算本輪優(yōu)化所得的總發(fā)電量E′，判斷| E-E′|＜ε是否成立，如果成立，認(rèn)為算法收斂，終止計算，輸出結(jié)果。否則令E=E′，轉(zhuǎn)至步驟（2）。

步驟（3）實質(zhì)上是一個有約束的多變量單值函數(shù)的最優(yōu)化問題。對于此問題，POA算法采用給定離散步長的網(wǎng)格搜索法求解子優(yōu)化問題。當(dāng)水庫數(shù)量較少時可以在子優(yōu)化問題中獲得較為精確的全局最優(yōu)解，然而隨著水庫數(shù)目的增加，計算規(guī)模呈現(xiàn)指數(shù)增長［23］，面臨嚴(yán)重的“維數(shù)災(zāi)”問題。為解決該問題，本文使用離散梯度下降法解決步驟（3）的子優(yōu)化問題，并提出離散梯度逐步優(yōu)化方法（DGPOA）。

4 離散梯度逐步優(yōu)化算法

4.1 梯度下降法梯度下降方法是無約束優(yōu)化問題中最常用的方法之一。設(shè)函數(shù) f：Rn→R，如果在每個變量xj均具有偏導(dǎo)數(shù)，則對于無約束最優(yōu)化問題：

負(fù)梯度方向-?f()xn為其最速下降方向，序列x0，x1，x2，…，滿足：

可得：

對于凸函數(shù)，序列（xn）會收斂到全局最優(yōu)解，非凸函數(shù)會收斂到局部最優(yōu)解。

4.2 離散梯度下降法使用梯度下降法的一個重要前提是函數(shù) f(x)在求解域中是可導(dǎo)的。然而，對于水庫梯級優(yōu)化，大部分特征曲線（如水位庫容曲線，水頭耗水率曲線，尾水位泄量曲線等）為實測點集形成的分段函數(shù)，而約束條件和目標(biāo)函數(shù)則為這些特征曲線的非線性組合，難以直接求導(dǎo)。為了充分利用局部區(qū)域的梯度信息，優(yōu)化搜索方向，加快搜索速度，本文結(jié)合水庫調(diào)度的特點，提出給定離散步長下的離散梯度的定義及計算方法，以及其對應(yīng)的離散梯度下降法。

4.2.1 離散梯度計算 對于多變量單值函數(shù) f(x1，x2，…，xn)，給定離散步長h＞0的情況下，令：

定義1。離散步長h下的近似偏導(dǎo)數(shù)：

該近似偏導(dǎo)數(shù)有如下性質(zhì)：

在定義1的基礎(chǔ)上，可以進(jìn)一步定義離散梯度。

定義2。給定離散步長h下離散梯度：

如果 f在x=a點存在梯度，則有：

4.2.2 搜索步長的確定 對于水庫調(diào)度，優(yōu)化的狀態(tài)變量（如水位）和目標(biāo)（如發(fā)電量）常常具有不同的量綱，導(dǎo)致步長容易出現(xiàn)過大或過小的情況，繼而影響算法收斂和搜索速度。本文搜索步長根據(jù)離散步長h以及離散梯度值進(jìn)行動態(tài)設(shè)置。

設(shè)：

則設(shè)搜索步長為

線搜索迭代方程為

4.2.3 停止準(zhǔn)則

或者

對于離散情況，在某個精度下，有可能出現(xiàn)沿離散梯度負(fù)方向搜索函數(shù)值變大的情況，此時最優(yōu)值一般位于an，an+1之間，認(rèn)為已滿足該精度要求，可以停止搜索。

4.3 庫群中長期優(yōu)化調(diào)度總體求解流程

（1）以水庫水位作為狀態(tài)變量，設(shè)置水位離散步長h，調(diào)度階段總數(shù)T，發(fā)電量收斂精度ε。

（2）根據(jù)常規(guī)方法確定初始水位過程，令E為初始條件下的發(fā)電量。

（3）設(shè)置t=T-1。

（5）時段t的兩階段問題為有約束的最優(yōu)化問題：

式中：多變量標(biāo)量值函數(shù)P∶Rn→R，為階段t的發(fā)電量函數(shù)，C為約束集。

使用罰函數(shù)為常數(shù)值函數(shù)的內(nèi)罰函數(shù)法對約束進(jìn)行處理。即當(dāng)水位及由水位計算的關(guān)鍵量不滿足約束時，發(fā)電量設(shè)置為-M，其中M為足夠大的常數(shù)。

因此可定義函數(shù) f∶Rn→R：

并將（25）轉(zhuǎn)化為無約束最優(yōu)化問題：

（6）令n=0，設(shè)置離散梯度下降法初始解a0=。

（8）判斷t≤ 0，如果滿足跳到步驟（9），否則返回步驟（4）。

（9）根據(jù)優(yōu)化的水位計算發(fā)電量E′，判斷|E′-E|≤ε，如果滿足則輸出結(jié)果，結(jié)束計算，否則令E=E′，返回步驟（3）。

圖1DGPOA計算流程圖

流程見圖1。

5 應(yīng)用算例及分析

5.1 工程背景及參數(shù)選擇以云南省瀾滄江干流5座主要電站為例對本文模型算法進(jìn)行驗證。瀾滄江流域水能資源豐富，是我國重點開發(fā)的13個水電基地之一。包含的5座電站從上游到下游依次為小灣、漫灣、大朝山、糯扎渡和景洪，總裝機容量達(dá)到14700 MW，是云南重要的調(diào)峰調(diào)頻電源及西電東送電源點。除小灣、糯扎渡為多年調(diào)節(jié)外，其余均為季調(diào)節(jié)，各水庫基本情況和始末水位設(shè)置值見表1。模型算法采用Python語言編寫，如果沒有特殊說明，初始解均由等流量法生成，設(shè)置發(fā)電量收斂精度ε=0.001億kW·h。

5.2 不同精度下的優(yōu)化結(jié)果分析利用平水年數(shù)據(jù)分別使用POA、POA-DPSA［20-21］和DGPOA對瀾滄江流域優(yōu)化調(diào)度。方案1將瀾滄江流域簡化為小灣、糯扎渡兩庫系統(tǒng)，計算結(jié)果見表2。方案2對小灣、漫灣、大朝山、糯扎渡和景洪等五水庫系統(tǒng)進(jìn)行計算，計算結(jié)果見表3。兩方案三種算法計算發(fā)電量隨離散精度變化見圖2。定義POA計算發(fā)電量e0，POA-DPSA計算發(fā)電量e1，DGPOA計算發(fā)電量e。POA-DPSA與POA計算發(fā)電量的相對偏差，DGPOA與POA計算發(fā)電量的相對偏差DGPOA與POA-DPSA的相對偏差

表1 瀾滄江梯級水庫基本情況與始末水位設(shè)置值

表2 方案1不同精度下POA、POA-DPSA、DGPOA發(fā)電量對比

表3 方案2不同精度下POA，POA-DPSA，DGPOA發(fā)電量對比

從計算結(jié)果上可以看出，方案1兩庫系統(tǒng)中，POA發(fā)電量e0在各離散步長下均處于最優(yōu)。POADPSA計算發(fā)電量e1在不同離散步長下依次為0.19%、0.24%、0.15%和0.21%，與POA相對差距總體在0.20%左右。DGPOA計算發(fā)電量e2與POA的相對差值則從0.53%逐漸減少到了0.17%。如圖2（a）所示，DGPOA與POA-DPSA計算結(jié)果差距隨著精度的提高逐漸減少，并且在精度為0.05 m時DGPOA計算結(jié)果要優(yōu)于POA-DPSA。值得注意的是，并非離散步長越小，優(yōu)化結(jié)果越好，POA在離散步長從1 m減少到0.5 m后優(yōu)化結(jié)果反而出現(xiàn)了下降，而POA-DPSA最優(yōu)值出現(xiàn)在離散步長為0.5 m的時候。這是由梯級水庫優(yōu)化調(diào)度問題的非凸特性導(dǎo)致的［24］。

方案2五水庫系統(tǒng)中，POA在精度較高時，已經(jīng)表現(xiàn)出了明顯的“維數(shù)災(zāi)”問題，僅在h=2 m時可以在有限時間（1d）內(nèi)計算出結(jié)果e0=795.01億kW·h。如圖2（b）所示，DGPOA計算結(jié)果在較高精度下要優(yōu)于POA-DPSA計算結(jié)果。當(dāng)離散精度h≤0.25m時，DGPOA計算結(jié)果要優(yōu)于POA在離散精度h=2 m時的值。因此，DGPOA可以通過提高精度獲得與POA算法相近甚至更好解。

5.3 不同精度下的求解效率分析方案1和方案2各精度下不同算法計算效率分別如表4和表5所示。將方案1和方案2離散精度和相應(yīng)算法計算時間利用冪函數(shù)進(jìn)行擬合，可得圖3。定義POA計算時間為t0，POA-DPSA計算時間為t1，DGPOA計算時間為t2。

從計算時間上可以看出，方案1中各算法計算時間均隨著離散步長h的減少而增加。在不同精度下，DGPOA算法計算時間分別是POA算法計算時間的2.01%、1.03%、0.40%、0.27%、0.09%和0.04%，相對POA算法的加速比t0/t2達(dá)到了50～250倍。DGPOA與POA-DPSA算法計算時間t2，t1的比值相對穩(wěn)定，依次達(dá)到了13.70%、14.07%、11.33%、15.53%、12.99%和11.84%，加速比t1/t2可達(dá)7到9倍。方案2中，POA算法已經(jīng)表現(xiàn)出了明顯的“維數(shù)災(zāi)”問題，在精度h=2 m時的計算時間已經(jīng)達(dá)到了4523.2s。在精度進(jìn)一步提高后，其計算時間已經(jīng)無法滿足常規(guī)的使用要求。DGPOA計算時間相對于POA-DPSA分別占用17.89%、16.65%、13.01%、10.45%、5.22%和6.24%，計算加速比隨著精度的提高有明顯增加的趨勢。

圖2 各算法計算發(fā)電量隨離散精度變化

表4 方案1各精度下不同算法計算時間對比

表5 方案2各離散步長下不同算法計算時間對比

如圖3（a）所示，在方案1中，POA算法的計算時間大致與離散步長的平方成反比，POA-DPSA、DGPOA大致與離散步長的一次方成反比。這說明隨著離散步長的減少，POA算法計算時間的增長速度將遠(yuǎn)超另外兩種算法。這也體現(xiàn)出了線搜索類算法的優(yōu)勢。

如圖3（b）所示，方案2中DGPOA計算時間大致與離散步長的0.9次冪成反比，DGPOA計算時間大致與離散步長的1.25次冪成反比。由此可得，隨著精度的提高DGPOA相對POA-DPSA的計算速度優(yōu)勢將更加明顯。

3種算法搜索過程有所不同，POA-DPSA算法采用了坐標(biāo)方向的線搜索，不需要對網(wǎng)格內(nèi)所有解進(jìn)行遍歷，“維數(shù)災(zāi)”問題相較于POA有所緩解。DGPOA相對POA-DPSA和POA，利用了局部離散梯度信息，在優(yōu)化方向上進(jìn)行優(yōu)化，搜索路徑更短，極大提高了搜索速度。

5.4 不同初始解下的結(jié)果分析采用不同配比的出庫流量生成小灣，糯扎渡雙庫系統(tǒng)的初始水位解，并以此為基礎(chǔ)分別用POA、POA-DPSA和DGPOA算法進(jìn)行水庫優(yōu)化調(diào)度。其中離散步長h均選為0.05 m。瀾滄江流域根據(jù)徑流特征可以分為汛前期（1—5月）、汛期（6—10月）和汛后期（11—12月）。初始解根據(jù)各時期各月給定配比的出庫流量求出。各時期各月出庫流量配比及計算結(jié)果見表6。其中初始解e的流量配比等價于等流量法，即同一水庫各月出庫流量均相同。

圖3 計算時間與離散精度關(guān)系

表6 初始解生成方式及計算結(jié)果

不同流量配比會生成不同的初始解，不同初始解下的計算結(jié)果可以反映算法的全局搜索能力。從計算結(jié)果上看，POA算法在初始解a—h下計算結(jié)果均為524.7億kW·h。POA-DPSA算法在初始解a—f下計算結(jié)果均為523.59億kW·h，在初始解g和h下，計算結(jié)果有少量提高，分別為523.72億kW·h和523.86億kW·h。DGPOA算法在初始解a—h下隨著汛期各月份流量配比提高，計算結(jié)果從522.70億kW·h提高到524.04億kW·h，總體優(yōu)于POA-DPSA。在初始解i和j下，3種算法計算結(jié)果均出現(xiàn)了明顯下降，下降幅度DGPOA最大，POA-DPSA其次，POA最小。由此可得，POA算法求解結(jié)果總體更優(yōu)，對初始解依賴相對較少結(jié)果相對穩(wěn)定，全局搜索能力最強。POA-DPSA算法和DGPOA各有優(yōu)劣。如果初始解為隨機生成，在與最優(yōu)解差距較大時，POA-DPSA結(jié)果相對穩(wěn)定。而在POA算法常用的初始解生成方法（等流量法）及其相鄰搜索域（初始解d—h），DGPOA搜索效果相對更好。

5.5 不同來水條件下的結(jié)果分析分別采用枯水年（75%）、平水年（50%）和豐水年（25%）徑流資料和POA-DPSA和DGPOA算法對瀾滄江流域5個主要水庫進(jìn)行優(yōu)化調(diào)度。由于POA在5個水庫系統(tǒng)中已經(jīng)表現(xiàn)出了嚴(yán)重的“維數(shù)災(zāi)”問題，此處不再對其進(jìn)行計算。離散步長h都選為0.1 m。優(yōu)化結(jié)果見表7。從計算結(jié)果上看，DGPOA算法和POA-DPSA算法在枯水年分別為710.24億和710.31億kW·h，豐水年分別是877.93億和878.07億kW·h，POA-DPSA整體較優(yōu)但差距不大。在平水年，DGPOA計算結(jié)果則要優(yōu)于POA-DPSA，發(fā)電量分別為794.85億和793.19億kW·h。

表7 不同來水條件下各算法計算結(jié)果對比

從計算效率上看，在枯水年、平水年和豐水年水平下，DGPOA算法的計算時間t2分別為28.49、36.07和27.73 s，POA-DPSA算法的計算時間t1分別為320.57、444.25和264.40 s，加速比t1/t2達(dá)到了9～12倍。年平均水平DGPOA算法優(yōu)化調(diào)度結(jié)果見圖4。

圖4 平水年下DGPOA算法的優(yōu)化調(diào)度結(jié)果

從圖4可以看出，兩個多年調(diào)節(jié)水庫小灣和糯扎渡在汛前均加大水位削落，充分發(fā)揮了調(diào)節(jié)庫容，對其他水庫進(jìn)行補償調(diào)度。進(jìn)入汛期后，小灣和糯扎渡開始蓄水以保證汛后期的高水位運行。漫灣、大朝山和景洪3個季調(diào)節(jié)水庫，在小灣和糯扎渡的調(diào)蓄下，始終處于高水位運行狀態(tài)，確保了梯級的整體發(fā)電效益。

6 結(jié)論

為解決水電站群中長期優(yōu)化調(diào)度的“維數(shù)災(zāi)”問題，在提出離散梯度的概念基礎(chǔ)上，使用離散梯度下降法對POA兩階段子優(yōu)化問題進(jìn)行求解，提出了離散梯度逐步優(yōu)化算法。該方法對待求問題的函數(shù)形式?jīng)]有特別要求，可利用局部離散梯度信息，最優(yōu)化搜索方向，快速獲得優(yōu)化結(jié)果。在瀾滄江流域五水庫系統(tǒng)多方案的計算結(jié)果表明，在工程常規(guī)精度下，DGPOA在不顯著降低全局搜索能力的情況下，計算速度可達(dá)POA-DPSA算法的8～12倍，POA算法的50～250倍。并且隨著水庫數(shù)目的增加和計算精度的提高，其相對POA-DPSA和POA速度提升會更加顯著。因此，DGPOA是解決水電站群中長期優(yōu)化調(diào)度“維數(shù)災(zāi)”問題的有效算法。