梯級水電站群高水頭巨型機組不規則多限制區自動規避方法

程春田，趙志鵬，靳曉雨，劉令軍，鐘儒鴻

（大連理工大學，遼寧 大連 116024）

1 研究背景

最近二十多年以來，我國在西南集中建成了以金沙江、瀾滄江、雅礱江、烏江、大渡河、紅水河等干流梯級水電站群為代表的巨型梯級水電站群，高水頭、大容量、巨型機組、遠距離、跨省跨區域輸送是其不同于以往中小流域的根本特征。由于高壓瞬變流，這些干流梯級巨型水電站普遍存在多個不規則機組限制區。當水電站或機組響應單個甚至多個受端電網負荷需求及水電電力市場化后時段出力波動頻繁時，由于水電站時段間和梯級上下游間緊密的水力聯系，容易引起水電站自身后續時段和梯級上下游水電站群出力、水頭和流量發生級聯變化，從而給水電站和梯級水電站群運行帶來巨大困難，嚴重影響水電站和電網安全經濟運行，成為制約西南水電安全經濟運行的技術瓶頸。

事實上，水電機組限制區問題一直是國內外水電機組組合（Hydro Unit Commitment，HUC）的經典問題，以往大多數研究主要集中在規則限制區建模與求解。Guan等[1]利用拉格朗日松弛法對規則限制區約束進行處理，將梯級水電調度問題轉化為一個由動態規劃和網絡流算法解決的雙層優化問題，取得了良好的結果。Bo等[2]、Borghetti等[3]和 Diniz等[4]均采用混合整數線性規劃（MILP）方法對規則限制區進行建模，對梯級水電機組組合問題進行了求解。隨著我國西南地區巨型梯級水電站群陸續建成，巨型機組高水頭不規則多限制區問題已經成為電網、發電企業面臨的棘手問題，基于工程中發現的問題，程春田等[5-6]提出了將機組不規則限制區通過組合數學的方式轉化成電站的限制區，通過廠站組合和啟發式策略規避機組限制區問題，但對于變化頻繁、要求更高的機組組合問題，上述方法還是略顯不足，需要更高效的多個高水頭不規則限制區快速自動規避方法。對于有復雜不規則機組限制區的HUC問題，國際上普遍采用MILP對問題進行建模求解[7-8]。MILP方法因為相對成熟的數學理論和良好的全局搜索能力、靈活的建模方式及有大量成熟的開源和商業求解器可以調用，是求解HUC問題最常用的數學規劃方法之一。Li等[7]采用基于人工剖分的方式對三峽電站存在的限制區進行MILP建模，實現了對32臺機組組合的高效求解，但該方法依賴特定電站的限制區特征和人工經驗，不適合自動化優化建模。Cheng等[8]提出一種基于分段線性化的不規則限制區MILP建模方法，實現了對烏江流域構皮灘單個電站的不規則限制區的建模與求解。但是其限制區邊界描述需要事先由人工確定，對限制區數據形式有一定要求，難以實現模型的自動化建模及應用。

為解決機組多不規則限制區快速規避和自動建模問題，本文以水電系統運行中常見的調峰需求為目標，提出基于Hertel-Mehlhorn凸剖分算法的多個不規則限制區約束自動解析技術、并應用凸優化理論及析取規劃理論方法構建復雜不規則限制區約束的MILP模型，然后應用商用求解器對問題求解，從而實現了復雜水電系統調度運行自動建模和問題求解。提出的方法在西南某干流流域梯級水電站群進行了方法驗證，下面將詳細介紹方法原理和具體的驗證情況。

2 高水頭巨型機組多個不規則限制區問題及其約束自動建模處理

據統計，世界上壩高200 m以上的水電站70%在中國，單機容量達到或者超過700 MW的水輪機組大多在中國，這樣就產生了非常獨特的西南地區高水頭巨型機組多個不規則限制區問題。不同于以往中低水頭和中下容量的機組規則限制區情況（圖1（a）），高水頭巨型機組多個不規則限制區邊界條件復雜（圖1（b）（c）），不能在工程應用通過簡單設定固定的出力條件就可快速解決此問題，這個問題隨著水電站或者機組對頻繁變化的負荷需求的響應而愈加困難，因此成為西南地區水電系統運行的瓶頸問題。為了自動識別和避開復雜的機組限制區，我們在此提出機組限制區約束的通用數學表達。

圖1 限制區示意圖

2.1 機組限制區約束的數學化定義首先，假設各限制區外邊界均為簡單多邊形，即任何不相鄰邊不相交。該假設符合目前已知的水電機組限制區特征。在實際工程中，可以直接獲取的數據包括：限制區邊界點，機組出力上下限，機組凈水頭上下限等。根據這些數據，進行如下定義：

式中：Rm表示機組第m個子限制區，R表示機組的限制區組合，Poly()表示由括號內點集依次相連，首尾相接組成的有界多邊形平面區域。(Hm，j，Pm，j)表示機組第m個限制區上的第j點。 M 表示機組含有的子限制區個數，Jm表示機組限制區m所包含的頂點數。為凈水頭和機組出力上下限形成的平面區域，其中“×”表示笛卡爾積乘法符號，和分別為機組凈水頭的下限和上限，分別為機組出力下限和上限。“”表示集合的減運算。 Asafe表示安全運行

區。在上述定義下，限制區約束可表示為：

式中：h、p分別為機組運行時的凈水頭和出力。

通過上述數學表達，我們將機組限制區約束統一描述為去除不規則限制區的相應安全區約束，從而為我們對不規則限制區識別及自動化建模打下基礎。

2.2 限制區約束自動建模方法通過限制區定義可以看出，去除限制區后的安全運行區域本質上是一個可能存在離散，有洞等復雜情形的極度不規則平面區域，如何實現對該區域的MILP自動化建模是本節及本文要解決的主要問題。為此，本文提出一種基于凸剖分算法及凸優化理論和析取規劃理論方法的不規則限制區自動化建模化方法。該方法首先需要對不規則安全區進行凸剖分，然后對剖分后的結果采用凸優化理論和析取規劃理論進行建模。

2.2.1 安全運行區Asafe的凸剖分 對安全區進行凸剖分是指將不規則安全運行區剖分成若干個互不重疊的凸多邊形區域的過程，而且要求這些剖分后的凸多邊形區域的并集與安全運行區相等。線性化建模過程實質上是對剖分結果的建模，剖分的結果會直接影響后續線性化建模質量。由于建模過程中引入的整數變量數目等于剖分后凸多邊形的個數（見2.2.2節），而整數變量數目又直接影響求解效率，因此為便于MILP求解，需要將不規則限制區剖分成盡量少的凸多邊形，這一問題在計算幾何中被稱為最優凸剖分（optimal convex decomposition，OCD）問題。OCD是典型的NP-hard問題，對該問題已有大量研究[9-10]，其中Hertel-Mehlhorn（HM）算法[10]因其計算復雜度低，且在大多數情形下可找到最優凸剖分結果，被廣泛應用到計算機游戲設計[11]、室內導航[12]等領域，本文在此采用HM算法實現對不規則限制區的凸剖分。使用HM算法對安全運行區進行凸剖分流程歸納如下：

圖2 安全區凸剖分示意圖

（1）預處理：由于限制區的復雜性，安全運行區可能是由多個多邊形組成，每個多邊形可能存在單個甚至多個洞（見圖2（a）），這些復雜的情形都不適用于HM凸剖分算法，因此需要對其進行預處理。預處理主要包括分離和去洞兩個操作。分離是指將包含有多個多邊形的情形分離成若干個多邊形，之后的所有操作均是對單個多邊形的操作。去洞是指將包含洞的分離后的多邊形轉化為不含洞的簡單多邊形的過程。該過程首先需要查找所有洞的最右側點，然后在距離該點較近的多邊形點之間進行分割，分割后多邊形洞的個數減少1。反復執行該過程，即可去除所有洞[13]。圖2（b）為去洞后安全區示意圖，從圖中可以看出，通過圖中右上角線段的分割，安全區不再包含洞。

（2）三角化：三角化是指將分離后的簡單多邊形劃分成若干互不重疊的三角形的過程。該過程采用耳切法（ear clipping，EC）[13]進行處理。對于簡單多邊形“耳”，指凸點與相鄰點圍成的三角形，且該三角形內部不可包含其他頂點。如圖3所示，圖中多邊形共包含4個耳，將三角形用其頂點構成的三元元組進行表示，則圖中多邊形的4個耳分別為（1，2，3）（2，3，4）（6，7，8）（7，8，9）。可以證明，任何超過3個頂點以上的簡單多邊形必然包括兩個以上的“耳”[13]，因此可以通過不斷切除多邊形耳的方式實現對簡單多邊形的三角剖分。圖2（c）為采用“耳切法”三角化后的安全區示意圖。

圖3 多邊形的“耳”

（3）去除非重要對角線：非重要對角線指去除后相鄰三角形的并集為凸多邊形的對角線，反之為重要對角線。

（4）重復步驟3，直到所有對角線均為重要對角線。在去除非重要對角線后，其他對角線的重要性可能會隨之改變，因此需要反復執行步驟2，直到不存在非重要對角線為止。圖2（d）為最終所得的安全區凸剖分結果，圖2（d）中比圖2（c）缺少的對角線，即為去除的非重要對角線。

2.2.2 線性化建模 本節進一步介紹如何根據前文所得的凸剖分結果實現對復雜限制區的MILP建模。為簡化表達，首先定義NN為不大于N的正整數集合，其中N為任意正整數。設經過凸剖分后，Asafe被剖分為凸多邊形集合根據 Acx，限制區約束（4）可進行如下轉化：

式中∨為邏輯“或”運算符號。

i的總邊數。根據凸優化理論，凸多邊形可以表示成以其各邊為界限的半平面的交集[14]，因此Acxi可轉化為下式：

式中：ai，j為 Acxi邊j的外法向量；bi，j為使上式成立的常數項。

根據式（6），式（5）可進一步轉化為：

上式中右側部分為典型的析取式結構，可有效被析取規劃方法處理[15-16]。因此進一步引入析取規劃建模方法對該析取式進行線性化建模。其中析取式是指由邏輯“或”（OR）運算符號∨連接的若干不等式或等式的關系結構。析取規劃方法是研究如何求解析取式約束的通用數學規劃方法。通過引入整數變量，將析取式結構轉化為求解器可以求解的合取式結構是析取規劃中最常用的處理方式。其中合取式是指由邏輯“與”連接的若干不等式或等式結構。析取規劃中對析取式轉化方法主要分為Big-M（BM）法[17-18]和凸包法兩種。其中BM法引入變量相對較少，計算效率通常較高[19]。因此本文引入BM法對式（7）進行線性化建模。

采用BM法進行構建，可得：

式中：yl為 Acxl的指標變量，如果 yl=1，式（8）中所有i≠l的約束將被BM常數松弛，僅保留i=l的約束項，此時為BM常數。式（8）現代求解器求解MILP問題通常采用分支定界或其變形方法，分支定界法首先求解原問題的線性松弛或部分線性松弛問題，因此其線性松弛問題的可行域與原問題越接近，越有利于問題的求解。顯然，較大的BM常數會導致其線性松弛問題可行域過大，從而降低分支定界算法的求解效率。因此，在滿足原問題結構的前提下，選擇盡量小的BM常數有利于進一步提高算法的求解效率。基于以上，本文BM常數取值方法如下：

該取值的可行性驗證如下：如果 yl′=1，則式（8）中滿足i=l′的約束組合等價于式（6），即對于式（8）中i≠l′的式子，將式（11）帶入其右手側常數項可得：

因此，該取值方式滿足式（8）結構。顯然，任何小于式（11）中的BM取值方法都有可能導致式（12）、式（8）的不成立，因此本文采用的BM值為該構建方法下的最優取值。至此，不規則限制區的線性化模型構建完畢。

3 具有多個不規則限制區的梯級水電站群優化調度模型

為驗證本文提出的處理機組復雜多不規則限制區方法的有效性，以調峰任務為目標，對梯級水電站群中多個水電站存在復雜限制區的調度任務進行建模。

3.1 目標函數調峰的目的是盡量使電網剩余負荷平緩，以減少系統整體峰谷差，降低火電等電源的爬坡壓力和燃料損耗。本文采用一階平均距作為目標函數以平緩剩余負荷。

式中：r、R為電站編號和電站總數，并規定編號較小者位于上游；t、TI為調度時段編號和總時段數，本文TI設置為24；u、Ur為機組編號和電站r的機組總數目。如無特殊說明，下文中r、t、u為下標時，均為電站，時段及機組編號；p′r，u，t為電站r機組u在時段t的出力；Dt、D′t分別為時段t的系統負荷和減去水庫調峰負荷的系統剩余負荷；----D′為系統剩余負荷平均值。上述形式無法直接采用求解器求解，本文引入輔助變量δt做如下等價轉換：

3.2 常規約束常規約束指水電機組組合問題中除了限制區約束以外的其他常見約束。

（1）水量平衡方程：

式中：qr，u，t表示機組開機時段平均發電流量;為時段末庫容下限及上限值；為時段平均出庫流量下限及上限值；為開機時段平均發電流量下限及上限值；為時段平均出力的下限和上限值。表示調度期初和末庫容是給定的。以上為機組開機時各變量的邊界約束，而機組實際出力和實際發電流量僅在開機時需要滿足上述約束，關機時需要設置為0。首先定義0-1機組狀態變量yr，i，t，如果yr，i，t=1表示相應機組處于開機狀態，否則 yr，i，t=0。則上述情形可以描述為下列各式：

式中：q′r，u，t為機組實際時段平均發電流量。式（26）、式（28）分別表示處于開機狀態時，機組實際發電流量（出力）等于機組開機發電流量（出力），若處于關機狀態，這兩個約束將被松弛。式（27）、式（29）表示處于關機狀態時，機組實際發電流量和出力必須為0，若處于關機狀態這兩個約束將被松弛。

（3）機組開停機持續時段約束：定義二元機組開機操作變量 gr，u，t，若 gr，u，t=1表示該機組在該時段進行開機操作，否則 gr，u，t=0。定義機組關機整數操作變量dr，u，t，若dr，u，t=1表示該機組在該時段進行關機操作，否則dr，u，t=0。則機組開停機約束可表示為：

式中TGr，u、TDr，u分別為相應機組的最小開機和關機持續時段數。

（4）電站出庫約束：

（5）凈水頭相關約束：本節介紹與凈水頭計算相關的約束。

（6）機組動力函數約束：

式中φr，u為機組動力性能函數。

3.3 限制區約束

3.4 常規約束構建方法常規約束中式（35）—式（37）為一維非凸非線性曲線約束，式（39）為二維非凸非線性曲面約束。這兩類約束無法直接利用求解器求解，需要進行線性化處理。而這兩類約束線性化已有大量成熟方法，本文對于一維非線性約束及二維非線性約束均采用文獻[8]中所述方法，在此不再贅述。

4 應用實例

為驗證本文所提方法的合理與有效性，選用中國西南地區W流域干流梯級水庫系統中包含復雜不規則限制區的高水頭巨型水庫A和水庫B作為重點研究對象。W流域梯級水庫系統總裝機達8 GW，最高水頭達到200 m，是我國的十三大水電基地之一。水庫A和水庫B是流域中總裝機最大的兩座高水頭巨型水庫。其中水庫A位于上游，調節性能為季調節，共包含5臺機組，總裝機達1250 MW，最高水頭為140 m。水庫B位于梯級系統的下游，為多年調節水庫，共包含5臺機組，總裝機達到3000 MW，最高水頭達200 m。其中水庫A所有機組，水庫B中除4#機組外其他機組均含有大范圍不規則限制區。選取某年7、8、10、11月份某典型日實際數據進行一日24點模擬計算，其中7、8月份方案作為汛期代表，10、11月份方案作為枯期代表。所有方案的最小開停機約束均為4 h。

表1 方案主要參數

各方案其他主要參數見表1。為進一步降低搜索空間提高求解精度和效率，表1中凈水頭范圍采用凈水頭動態設置方式確定，具體方法如下：

本文算法及模型構建均采用Python3.6語言編寫，程序運行的操作系統為Ubuntu16.4虛擬機，硬件配置為Inte（lR）Xeon（R）CPU E7-4850 v3@2.20GHz 96 logic CPU，32G RAM，并調用Gurobi8.1求解器分支定界算法進行求解，凸剖分涉及的計算幾何相關基礎算法采用Shapely[20]和PolyPartition[21]等開源庫中的相關算法。模型中對常規非線性約束進行線性化時，尾水位泄量約束、機組性能曲線約束各變量方向分段數設置3，水頭損失曲線分段數為4。對于巨型水庫，其一日內的水位變幅相對較小，因此本文僅在始末水位上下1m的范圍內進行離散，分段數設置為1。設置算法的停止準則為運行時間達到1800 s或gap值達到0.02。其中gap值指當前最優可行解和最優值下限的相對差值，gap值是描述當前解最優性的指標，gap值越小說明當前值與理論上的全局最優值越接近。

4.1 凸剖分結果分析圖4為水庫A及水庫B各機組的安全運行區及其凸剖分結果示意圖。值得注意的是，采用動態凈水頭區間設置方法后，其安全運行區會隨著凈水頭搜索區間的改變而改變。搜索區間越小，相應安全運行區越小，因此方案計算過程中采用的安全運行區要小于圖4所示的安全運行區。不失一般性，本節僅給出在設計最大凈水頭和最小凈水頭之間的安全區的剖分結果圖。可以看出，安全運行區呈現出了高度不規則甚至有洞等特性。其中水庫A中1#—3#機組，水庫B中1#—3#及5#機組存在不止一個不規則限制區，此時安全運行區的不規則性更加顯著。水庫B的1#—3#及5#機組在多限制區影響下，呈現出有洞的特性。水庫A中4#機組在水頭位于[112，129]時，存在一塊規則限制區，但在更大的區域內，安全運行區仍然具有不規則性，以往研究往往顯性或隱性假設機組運行不會超出該規則限制區的水頭范圍，這在大多數情況是成立的，但在某些時段機組仍有可能在超出規則限制區運行水頭范圍內運行，那么針對規則限制區的建模不再適用。水庫B的4#機組不包含限制區，因此其運行區為完整的一塊矩形區域。從剖分結果看，在離散、有洞、多不規則限制區、部分規則限制區影響下，算法均可以有效將安全運行區凸剖分為若干互不重疊的凸多邊形，這體現出凸剖分算法的通用性。

圖4 水庫安全運行區及其凸剖分結果

4.2 計算結果及調峰效果分析各月份方案計算結果見表2。從表中可以看出，汛期方案整體耗水量和發電量多于枯期方案。從變量及約束數量看，4個方案中連續變量均為3599個。不同方案之間離散變量個數及約束個數不同，這是由于采用了動態凈水頭搜索區間后，安全運行區也會動態的隨之改變，從而導致安全區剖分結果及線性化結果的不同。從結果gap值及計算時間看，4個典型方案均在給定時間內求出了近似最優解，其中汛期方案的gap值達到0.02，枯期方案中10、11月份的gap值也達到0.04和0.05。在實際應用過程中，對求解時間敏感的場景下，可以通過降低最大求解時間的方式或提高gap限制值的方式，減少總計算時間。相應的，對結果最優性要求較高的場景，則可以通過增加最大求解時間及進一步降低gap的方式，以期得到更優的解。

表2 各月份方案計算結果

從圖5中直觀看出，4個方案均達到較為顯著的削峰效果。水庫A和水庫B的發電過程也較好響應了調峰需求。進一步的，各方案結果的具體調峰指標值如表3所示。表中表示原負荷的平均爬坡，表示剩余負荷的平均爬坡。這兩個值能夠體現負荷過程的整體光滑程度，其越小則認為越平滑，也越有利于火電運行。表示平均爬坡的減少比例。從表3中可以看出，汛期兩方案的峰谷差減少比例分別達到0.65和0.70，平均爬坡減少比例也分別達到0.54和0.63。枯期兩方案由于整體發電量較少，調峰效果相較汛期方案相對較差。枯期兩方案峰谷差相對減少比例為0.24和0.35，平均爬坡減少比例也達到了0.11和0.35。由此可以看出，本文的調峰模型可以有效削減峰谷差，使調峰結果更加平緩。

圖5 各方案調峰效果

表3 各方案結果調峰指標值

4.3 限制區規避效果分析本節分析本文所提模型計算結果對限制區的規避效果，并與常規計算模型進行對比。其中常規計算模型指不考慮限制區約束的短期調峰模型。如圖6、圖7所示，本節分別選擇8月份和10月份方案作為汛枯期代表進行分析。圖中本文模型限制區指本文模型計算的運行過程中各時段機組平均凈水頭對應的出力限制區，常規模型限制區同理可得。從圖6中可以看出，常規模型由于沒有考慮限制區約束，在水庫A的3#機組，水庫B的2#和5#機組均出現落入限制區的情況，從而對電廠及電網的安全造成威脅。而本文模型在各個時段均避開限制區。在汛期電網用電高峰時刻，各機組也以接近裝機容量的狀態進行發電，以滿足電網的調峰需求。對于枯期10月份方案，從圖7可以看出，常規模型在水庫A的2#、3#機組，水庫B的1#、2#、5#機組出現了長時段運行在限制區的情況。本文模型計算結果均滿足限制區約束，保證了電站和電網的安全穩定運行。綜合對比圖6和圖7的出力限制區過程可以看出，在不同運行條件，甚至同一日內出力限制區都發生劇烈變化，傳統規則限制區考慮方式不再適用，對其的簡化也勢必造成誤差從而增加落入限制區的風險。綜上，本文所提模型可以有效考慮水電機組復雜限制區約束，在保證電廠電網安全運行的前提下，充分發揮水電的調峰能力。

圖6 8月份方案本文模型和常規模型計算結果對比

圖7 10月份方案本文模型和常規模型計算結果對比

5 結論

針對制約影響我國西南地區水電安全經濟運行的高水頭巨型機組復雜多不規則限制區瓶頸問題，本文提出了自動解析處理機組復雜限制區的建模優化方法，該方法結合了計算幾何凸剖分理論、凸優化理論、析取規劃方法及MILP數學規劃方法，并通過西南某干流梯級水電站群實際資料進行了驗證，結果展現了方法的有效性。由于所提方法不需要對復雜不規則限制區進行人工預處理，完全根據限制區的數學定義實現自動剖分，達到了自動建模與求解的效果。因此，解決了區域、省級電網大規模水電系統復雜不規則限制區自動優化建模難題，這對于將來響應電網差異化、自動化調度應用，響應市場化條件下的水電出力頻繁變化的新情況，意義特別重大。