趙志鵬，劉 杰，程春田，廖勝利，靳曉雨

（大連理工大學(xué)，遼寧 大連 116024）

1 研究背景

我國正處于經(jīng)濟結(jié)構(gòu)改革的過渡期，加之清潔能源的大規(guī)模投產(chǎn)，電力系統(tǒng)面臨巨大的調(diào)峰壓力，水電作為一種可調(diào)節(jié)的清潔能源需要承擔(dān)更大的調(diào)峰任務(wù)［1-4］。水調(diào)系統(tǒng)信息化程度的提高對水電短期調(diào)度精細化的要求越來越嚴格［5］，尤其對于有靈敏水力聯(lián)系的梯級庫群，不僅水庫運行狀況（如出力、振動區(qū)、最大發(fā)電流量、水頭損失等）對水頭極為敏感，而且梯級之間還存在敏感的水力耦合，以往將一日內(nèi)水頭設(shè)為定值或采用平均水頭進行日調(diào)度的粗放模式難以滿足精細化調(diào)度的要求，因此需要充分考慮日內(nèi)的水頭變化。其中回水頂托是存在于耦合梯級庫群間的一種復(fù)雜的水力聯(lián)系，正常情況下，水庫的尾水位和出庫流量存在一條穩(wěn)定的關(guān)系曲線，但當(dāng)梯級水庫上下游壩址較近時，較高的下游庫水位產(chǎn)生回水，上游水庫的尾水位于回水區(qū)間內(nèi)時，上述穩(wěn)定的水位流量關(guān)系曲線將被破壞，這就是所謂的回水頂托影響［6-8］，回水頂托在我國西南水電庫群間普遍存在，如三峽-葛洲壩，溪洛渡-向家壩，天生橋一級-二級等?；厮斖械闹苯佑绊懯翘Ц呱嫌嗡畮煳菜?，如果在日前調(diào)度計劃中未將其考慮在內(nèi)，將導(dǎo)致上游水頭計算偏高與實際運行工況不符，對于上下游電站有利益沖突的梯級，甚至?xí)?dǎo)致計劃發(fā)電量與實際發(fā)電量存在沖突而產(chǎn)生經(jīng)濟糾紛［9］。因而，針對有靈敏水力聯(lián)系的梯級庫群，如何將下游回水頂托的影響考慮到模型中，并實現(xiàn)模型的高效求解，是當(dāng)前水庫精細化調(diào)度的重點及難點之一。

水電日前調(diào)峰優(yōu)化調(diào)度需要同時考慮電網(wǎng)、水庫及機組多重約束，本質(zhì)上屬于復(fù)雜的多維、動態(tài)、非凸、非線性規(guī)劃、強耦合問題［10］。將回水頂托影響考慮到模型之后，該問題的耦合性、非凸、非線性等特征進一步凸顯，極大地增加了求解難度。目前短期調(diào)峰模型的研究較多，通常采用余荷均方差最小或余荷最大值最小為目標(biāo)建立模型，由于目標(biāo)函數(shù)的特殊性多采用啟發(fā)式搜索方法［11-13］，解的穩(wěn)定性難以保障，且大多研究簡化水頭計算，將回水頂托對發(fā)電水頭的影響考慮在內(nèi)的研究少之又少。文獻［14］在解決三峽機組組合時考慮葛洲壩回水頂托影響，在葛洲壩水位已知的情況下擬合二維水位流量曲線并采用迭代法求解尾水位；文獻［15］雖然考慮了電站出力的水頭敏感性，但其尾水位的計算或采用定值或采用下游電站壩上水位；文獻［16］提出了回水頂托的三維關(guān)系曲線，根據(jù)下游壩上水位的范圍采用二進制變量選取某個水位流量曲線，進而將三維關(guān)系轉(zhuǎn)換成二維關(guān)系進行求解，但其求解精度存在一定的偏差。目前調(diào)度機構(gòu)在進行日前調(diào)峰計劃編制時，大多將一日劃分成96 個時段進行調(diào)節(jié)計算。而以上研究雖在一定程度上考慮了回水頂托對調(diào)度的影響，但由于調(diào)峰問題及復(fù)雜水力耦合本身存在的復(fù)雜性或所采用算法本身的局限性，大多將一日簡化為24 個時段進行優(yōu)化計算，計算結(jié)果的實用性相對較差。

混合整數(shù)線性規(guī)劃（MILP）是水庫發(fā)電調(diào)度領(lǐng)域最為常用的數(shù)學(xué)規(guī)劃算法之一［17-19］。相較于傳統(tǒng)的啟發(fā)式算法，MILP 具有理論成熟，結(jié)果穩(wěn)定性及全局收斂性好，以及有大量先進的開源及商業(yè)求解器可以直接調(diào)用等特點。因此，本文建立考慮回水頂托影響的梯級庫群日前調(diào)度的96 點混合整數(shù)線性規(guī)劃（MILP）模型。模型以余荷平均距平絕對值最小為目標(biāo)，通過函數(shù)聚合方式減少常規(guī)及有回水頂托水庫的非線性因素。聚合后，常規(guī)水庫和有回水頂托水庫的非線性因素分別轉(zhuǎn)化為三維曲面約束和四維超曲面約束，并引入SOS2（Special Order Set of type 2）約束建模方法對其進行線性化處理。最后以西南流域某梯級為例進行驗證所提算法及模型的有效性。

2 數(shù)學(xué)模型

2.1 目標(biāo)函數(shù)采用余荷平均距平絕對值最小為目標(biāo)函數(shù)以平緩負荷峰谷差：

式中： r 為電站自上游至下游編號； R 為總電站數(shù)； t 為時段編號；T 為總時段數(shù)，本文將一天劃分為96 個時段； Dt為時段t 系統(tǒng)面臨負荷； D′t為時段t余荷；為余荷平均值； pr,t為電站r 在時段t 出力，如無特殊說明，下文r 均代表電站編號，t 均代表時段編號。

目標(biāo)函數(shù)中包含絕對值符號難以直接利用MILP 求解，引入非負輔助變量δt后，目標(biāo)函數(shù)可等價轉(zhuǎn)化為：

2.2 約束條件

（1）梯級日總電量約束：

式中： Δt 為時段間隔； E 為梯級日總電量。

（2）水量平衡約束：

式中：vr，t為庫容；Ir，t為區(qū)間入庫流量；ur，t為出庫流量，包括發(fā)電流量qr，t及棄水流量

（3）出力曲線約束：

（4）水位庫容約束：

（5）尾水位曲線約束：

（6）凈水頭計算公式及水頭損失約束：

式中： Δhr,t為電站r 在時段t 的水頭損失； ar為水頭損失系數(shù)。

（7）電站出力約束及爬坡約束：

（9）棄水流量約束：

（10）出庫流量約束：

（11）庫容約束：

（12）日調(diào)節(jié)水庫末庫容約束：

3 求解方法

3.1 多約束曲線聚合模型中存在諸多非線性因素，難于直接利用MILP 進行求解，分段線性化是MILP 求解非線性約束時常用的處理方法。模型中非線性約束包括式（10）—（12）和式（14）等。以往研究大多對這些約束直接進行線性化，每個約束的線性化過程均會引入大量0-1 整數(shù)變量，從而極大影響求解效率，甚至無法求解。為了避免引入過多的0-1 整數(shù)變量，分別針對常規(guī)無回水頂托水庫和有回水頂托水庫提出如下約束聚合策略：

3.1.1 無回水頂推水庫約束聚合 對于常規(guī)無回水頂托水庫，將式（10）—（14）進行聚合可得：

實際調(diào)度中要避免發(fā)生棄水，因而本文設(shè)定棄水流量為0，則ur,t=qr,t，對式（22）中相應(yīng)部分進行替換可得：

3.1.2 有回水頂推水庫約束聚合 對于有回水頂托水庫，將式（10）—（14）進行聚合可得：

式中，φ′i( )將( vr,t,qr,t,vr+1,t,pr,t)約束在四維空間的超曲面上，并包含了有回水頂托水庫模型中所有的非線性項，本文稱之有回水頂托水庫的等效動力函數(shù)，簡稱頂托動力函數(shù)。本文利用約束式（26）對有回水頂托水庫模型約束式（10）—（14）進行等效替換。

3.2 等效動力函數(shù)線性化常規(guī)動力函數(shù)為三維曲面函數(shù)，此類函數(shù)的線性化方法以往已有大量文獻進行研究［20］。對于頂托動力函數(shù)，本質(zhì)為四維超曲面函數(shù)，極大的增加了建模難度。本文引入SOS2 約束建模方法統(tǒng)一對這兩種動力函數(shù)進行線性化。為簡化表達，首先定義NI為不大于I的正整數(shù)集合，?0I為不大于I的非負整數(shù)集合，

3.2.1 SOS2 約束 SOS2 是一種特殊的非負實數(shù)集合，該集合只允許出現(xiàn)至多兩個非零實數(shù)，且這兩個非負實數(shù)必須是在順序上是相鄰的。定義有序?qū)崝?shù)集合對ηS施加SOS2 約束（表示為SOS2（ηS）），ηS即成為集合。當(dāng)S ＞1 時，SOS2（ηS）可以如下方式建模：

式中： bs為指示變量，如果bs=1，則表示除ηs-1≥0 ，ηs≥0 外，ηS中其他元素必須為零。如果S=1，則SOS2（ηS）退化為下式：

通過SOS2 約束建模方式，可以進一步構(gòu)建常規(guī)及頂托動力函數(shù)的線性化方法。

3.2.2 常規(guī)等效動力函數(shù)線性化 常規(guī)等效動力函數(shù)是一個典型的非凸三維曲面，常用的線性化方法主要包括兩類：一類是將曲面在其中一個維度上劃分成若干區(qū)間，然后在相應(yīng)區(qū)間上用一個一維曲線近似表示，該方法建模比較直觀，引入的二進制數(shù)量相對較少，但精度不足［18，21-22］。另一類是利用三角柵格化技術(shù)將曲面劃分成若干三角網(wǎng)格，利用三角網(wǎng)格頂點的凸組合近似表示相應(yīng)區(qū)域的曲面，該方法求解精度高，但是會引入大量0-1 整數(shù)變量，從而影響求解效率［14，23］。本文采用基于四邊形網(wǎng)格柵格化技術(shù)與SOS2 約束建模方法對常規(guī)等效動力函數(shù)進行線性化處理。

圖1 常規(guī)動力函數(shù)圖像在庫容-發(fā)電流量平面投影

圖1為常規(guī)動力函數(shù)圖像在庫容-發(fā)電流量平面的投影示意圖。為便于表達，變量、函數(shù)及常量均省略電站編號r 及時段編號t 下標(biāo)（下同）。如圖所示，柵格化時設(shè)置為庫容及發(fā)電流量坐標(biāo)方向上柵格化的格點集合，且滿足其中I 和J 分別為庫容及發(fā)電流量方向柵格化區(qū)間數(shù)。圖中表示對應(yīng)格點的出力值。引入λi,j非負輔助變量為格點的權(quán)重系數(shù)。則具體構(gòu)建方法如下：

圖2 頂托動力函數(shù)在庫容-發(fā)電流量-下游庫容三維空間投影

3.2.3 頂托等效動力函數(shù)線性化 頂托動力函數(shù)本質(zhì)是一個四維非凸超曲面函數(shù)，這進一步增加了建模難度。本文利用直六面體柵格化技術(shù)和SOS2 約束建模法對該函數(shù)進行線性化處理。

圖2為頂托動力函數(shù)四維曲面在庫容-發(fā)電流量-下游庫容三維空間的投影示意圖，其中v′為下游水庫庫容。與常規(guī)動力函數(shù)類似，柵格化時除設(shè)置為庫容及發(fā)電流量坐標(biāo)方向上柵格化格點集合外，進一步設(shè)置為下游庫容坐標(biāo)方向的柵格化格點集合，且滿足其中K 為下游庫容方向的柵格化間隔數(shù)，分別為下游水庫庫容下限及上限值。引入λi,j,k非負輔助變量為格點的權(quán)重系數(shù)。具體構(gòu)建方法如下：

其中式（40）—（43）將( v ,q,v′,p )表示成四維點集的凸組合。αi，βj，γk為引入的表達式變量，分別表示圖2中所示的庫容方向法平面i ，發(fā)電流量方向法平面 j ，下游庫容方向法平面k 的上格點對應(yīng)權(quán)重之和。式（47）依次表示對庫容，發(fā)電流量，下游庫容方向上各法平面權(quán)重施加SOS2 約束。在SOS2 約束下，至多允許庫容方向相鄰兩個法平面，發(fā)電流量方向相鄰兩個法平面以及下游庫容方向相鄰兩個法平面上格點對應(yīng)權(quán)重非零。取交集后，最終僅允許相鄰8 個格點對應(yīng)權(quán)重非零。此時點( v,q,v′,p )被限制在這8 個格點最小凸包內(nèi)部。當(dāng)柵格化格點足夠多時，點( v ,q,v′,p )將無限趨近于頂托等效動力函數(shù)四維超曲面。

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

4.1 場景選擇與模型參數(shù)設(shè)置本文選擇西南某流域電站A、電站B 組成的梯級庫群為研究對象，電站A、B 是南方電網(wǎng)重要的調(diào)峰調(diào)頻電源，壩址距離僅為7 km，電站A 位于電站B 的回水區(qū)間內(nèi)。電站A 調(diào)節(jié)性能為不完全多年調(diào)節(jié)，總裝機1200 MW；電站B 調(diào)節(jié)性能為日調(diào)節(jié)，總裝機1320 MW。模型均采用Python3.6 語言進行編寫，并調(diào)用Gurobi8.1 商業(yè)求解器MILP 分支定界算法進行求解。程序運行環(huán)境為Ubuntu16.4 虛擬機，硬件配置為Intel（R）Xeon（R）CPU E7-4850 v3 @2.20GHz 96 logic cpu，32G RAM。所有模型設(shè)置最大計算時間為3600 s，gap 參數(shù)設(shè)置為0.01，gap計算公式見下式：

采用如下累積誤差公式衡量模型線性化精度：

本文所提模型命名為本文模型，設(shè)置常規(guī)模型為對照組，常規(guī)模型直接采用傳統(tǒng)方法進行建模。傳統(tǒng)方法是指不經(jīng)過聚合，直接對2.2 節(jié)中非線性約束進行分段線性化建模的方法［19，24］。其中單變量非線性約束采用文獻［24］中的方法進行建模，雙變量非線性約束采用3.2.1 節(jié)中介紹的方法進行建模。模型中非線性約束相應(yīng)自變量分段線性的分段數(shù)設(shè)置見表1。由于水庫A 調(diào)節(jié)性能為不完全年調(diào)節(jié)，在一天內(nèi)的庫容變化較小，因此本文僅在其初始水位上下1 m 的庫容范圍內(nèi)進行離散，分段數(shù)設(shè)為1。其余變量均在其最大與最小值之間進行均勻分段。經(jīng)過聚合后，本文模型各水庫不再包含水位、凈水頭等變量，表格中相應(yīng)內(nèi)容用“-”表示。

表1 各變量分段數(shù)目

為體現(xiàn)算法通用性，選取該梯級歷史2月、3月、7月、8月典型日實際數(shù)據(jù)進行一天96 點模擬計算，其中2、3月為枯期代表，7、8月為汛期代表。電站B 為日調(diào)節(jié)電站，末水位等于始水位，實現(xiàn)日內(nèi)周期調(diào)節(jié)，各月份方案參數(shù)見表2。

4.2 結(jié)果分析對4 種典型日場景下均采用本文模型和常規(guī)模型進行建模計算，表3為本文模型及常規(guī)模型中連續(xù)變量，約束及離散變量總個數(shù)。從表3可以看出，通過聚合方式，本文模型的變量和約束數(shù)量均顯著小于常規(guī)模型，這為本文模型的高效求解提供了基本保障。

表2 各方案參數(shù)設(shè)置

表3 本文模型和常規(guī)模型中約束及變量總個數(shù)

4.2.1 調(diào)峰效果分析 圖3給出了本文模型在4 種場景的削峰結(jié)果圖，均達到不同程度調(diào)峰效果。進一步分析結(jié)果指標(biāo)，如表4所示，由計算結(jié)果可知，本文模型在滿足電量約束且不發(fā)生棄水的條件下均能找到可行解?？萜谔菁壙傠娏恳筝^小，均達到較好的調(diào)峰效果，2月、3月峰谷差減小率達31.41%和35.05%，達到較好的調(diào)峰效果；而汛期在不發(fā)生棄水的環(huán)境下，調(diào)峰能力有限，7月、8月峰谷差減小率僅為17.00%和1.00%。因而，在梯級日總電量較大的情況下，不棄水與調(diào)峰存在矛盾，而本文模型均能在不棄水條件下找到可行解。

表4 典型日場景調(diào)峰結(jié)果

圖3 本文模型各月份調(diào)峰效果

4.2.2 回水頂托效果分析 選取3月的計算結(jié)果對回水頂托影響進行分析。圖4（a）以出庫流量為參照，考慮回水頂托和不考慮回水頂托尾水位對比圖，圖4（b）以下游壩上水位為參照，考慮回水頂推和不考慮回水頂推尾水位對比圖。考慮回水頂托的結(jié)果由本文模型進行計算，不考慮回水頂托的結(jié)果是根據(jù)本文模型計算結(jié)果的下泄流量過程進行反推計算出的結(jié)果。從圖4可以看出，不考慮頂托的尾水位趨勢與出庫流量一致，不受下游壩上水位影響?？紤]頂托之后，尾水位同時受下泄流量及下游壩上水位影響，整體趨勢與下游壩上水位趨勢一致，因而電站A 的尾水位受下游電站壩上水位影響極大，在日常調(diào)度模型構(gòu)建中不可忽略。表5給出了詳細的指標(biāo)計算結(jié)果，可以看出，不考慮頂托與考慮頂托，水庫A 的水頭差距達4.67 m，累計多發(fā)電量433.14 MW·h，累計誤差達2.84%，實際操作若不考慮回水頂托，必然造成計劃與實際工況不符，不能按預(yù)期完成調(diào)度任務(wù)。

圖4 考慮回水頂托與不考慮回水頂托尾水位對比

表5 水庫A 考慮頂托與不考慮回水頂托計算差異對比

4.2.3 計算效率與精度分析 表6列出了本文模型，常規(guī)模型在4 種典型日場景中計算結(jié)果。設(shè)置最長計算時間為3600 s，表中g(shù)ap 為“-”時，表示沒有在設(shè)定時間內(nèi)算出可行解。從兩個模型的求解時間可以看出，本文模型在各個月份情形下均能夠算出有較高精度的可行解，且大部分月份（2、3、7月份）可以在較短時間內(nèi)算出滿足gap 要求的最優(yōu)解；常規(guī)模型僅2月份在3006 s 算出滿足gap 的可行解，3月份在設(shè)定時間內(nèi)甚至沒有得到可行解，由此可見，本文模型的計算效率要顯著優(yōu)于常規(guī)模型。對比目標(biāo)值可以看出，本文模型在大部分月份求解的目標(biāo)值要優(yōu)于傳統(tǒng)模型。對比累積誤差可以看出，各場景下常規(guī)模型的累積誤差均大于本文模型的累積誤差。因此，本文所建模型無論在精度還是計算時間及算法穩(wěn)定性上都優(yōu)于常規(guī)線性化模型，對MILP 在實際工程中的應(yīng)用具有重要意義。在實際應(yīng)用時，對于求解時間敏感的應(yīng)用場景，可以通過適當(dāng)減少最長運行時間限制或者提高gap 值的方式進一步降低總體運行時間。

表6 典型日場景計算結(jié)果

5 結(jié)論

為處理梯級庫群間存在的強耦合水力聯(lián)系，采用考慮下游壩上水位的尾水位計算公式將回水頂托影響考慮到模型當(dāng)中，把梯級庫群調(diào)度中多個非線性約束聚合成單個等效動力特性曲線，借助SOS2 約束對模型進行線性化，與傳統(tǒng)線性化方法相比，顯著減少了模型中的連續(xù)變量及離散變量數(shù)目，提高了模型的求解效率及精度。以西南某存在回水頂托的梯級庫群為例對所提模型及算法進行驗證，無論在求解精度還是求解效率上都優(yōu)于常規(guī)MILP 模型，對比考慮與不考慮回水頂托模型計算結(jié)果，不考慮回水頂托影響所制定的計劃造成較大累積誤差，因此，針對存在回水頂托的梯級庫群制定調(diào)峰日計劃時必須將其考慮在內(nèi)。