梯級水庫短期優化調度模型的精細化與GPU并行實現

紀昌明，馬皓宇，吳嘉杰，俞洪杰，彭 楊

（華北電力大學可再生能源學院，北京 102206）

1 研究背景

梯級水庫優化調度通過對入庫徑流進行調節以實現水電站的效益增長，一般分為中長期優化調度與短期優化調度，受限于當前水文與氣象預報技術的發展水平，中長期優化調度僅具有規劃意義，能指導水電站實際運營的只有短期優化調度［1］。而近幾十年的研究中，短期優化調度的理論與實際間的差距一直存在［2］，其主要原因在于為滿足發電計劃制定的時效性要求，以往的研究中人為引入了一系列策略以簡化計算模型，從而將計算時長降至可接受范圍內，但這也導致所得計劃在實際運行中出現偏差，而水電站出力的近似計算是造成偏差的主要原因［3］。為簡化廠內經濟運行，電站的出力可通過恒定出力系數的出力方程或忽略機組振動區的等微增率法計算獲得，而這種簡化計算方式將導致可用的水資源在時空分配上的不合理，通常不為水庫操作人員所接受。故本文構建同時考慮水庫優化調度與水電站優化運行的精細模型，通過經典算法—動態規劃（DP）［4-6］的二重嵌套形式，求解給定規模與離散精度下的優化調度方案。

DP 在求解精細模型時將面臨更為嚴重的“維數災”問題，如何處理該問題成為計算的關鍵。“維數災”本質上是程序計算量與內存占用量的指數型增長，水庫調度領域中許多學者僅針對前者提供了解決方案，主要分為兩類：（1）算法改進，通常需犧牲解的全局最優性以提高計算效率：馮仲愷等［7］為緩解“維數災”問題，結合均勻試驗設計提出均勻動態規劃；史亞軍等［8］提出一種基于灰色系統預測方法的離散微分動態規劃，以求解梯級水庫群優化調度問題；紀昌明等［9］為解決傳統動態規劃在處理水庫群聯合優化調度時面臨的約束處理機制選擇和計算時間長的問題，引入映射思想，基于映射和集合論知識構建可行域搜索映射模型；（2）CPU 并行加速，通常需花費高昂的成本以構筑計算集群：萬新宇等［10］在分析傳統串行動態規劃算法計算特點的基礎上，建立基于主從模式的并行動態規劃模型；Li等［11］利用高達350核的高性能CPU計算系統，實現了對等并行模式下的并行動態規劃；蔣志強等［12］基于.NET4 并行拓展庫，構建了時段間、時段內離散組合間以及兩者混合模式下的并行多維動態規劃算法。

區別于上述的研究成果，本文為全面且妥善地處理“維數災”，首先針對極少提及的內存占用量增長問題，采用數據壓縮與數據庫技術以減輕程序運行對計算機的內存負擔；其次針對重點研究的計算量與計算時間問題，采用在水動力模擬與分布式水文模型中已廣泛應用［13-15］，而在水庫調度中未引進的圖形處理器（GPU）并行加速［16］，以大幅降低程序的運行時間。最后以潘口、小漩兩庫聯合調度為研究實例，驗證所提模型、算法及優化策略的科學與實用性。

2 精細優化調度模型

本文選擇“以水定電”作為精細優化調度模型中各電站的運行模式，即未來一個調度期內，在各庫調度期初末庫容、各電站調度期初機組狀態、各階段區間徑流已知的前提下，同時考慮水庫調度的庫容、下泄流量等約束與水電站運行的機組振動區、啟停機耗水等約束，尋求最大化梯級系統總發電量的運行方式。

2.1 目標函數

式中：E*為調度期內梯級水電站的最大發電量；t、i分別為時段編號和水庫（水電站）編號；T、n分別為調度期階段數和水庫總數；Ni,t為t階段i庫對應的電站最優出力，其計算過程如圖1所示，水庫調度決策階段平均泄流Ri,t，由此推得階段平均水頭Hi,t，將兩者作為水電站運行的輸入，而電站通過決策本階段投入運行的機組與機組間的流量分配，以最大化電站出力，該過程表示為Ni,t=φi(Hi,t,Ri,t)；Δt為單位階段長度。

2.2 約束條件

（1）水庫調度約束。水量平衡約束：

式中：Vi,t、Vi,t+1分別為t階段和（t+1）階段i庫的階段初庫容；Qi,t、Ri,t分別為t 階段i 庫的階段平均入流和泄流值。

庫容約束：

下泄流量約束：

邊界約束：

式中：Vi,start、Vi,end為i庫給定的調度期初末庫容。

圖1 出力計算簡圖

（2）水電站運行約束。發電流量約束：

出力約束：

機組振動區約束：

機組啟停機耗水約束：

邊界約束：

2.3 模型求解該模型為二重嵌套優化結構，在調度領域也屬于相當復雜的模型，其外部結構（水庫調度）逐階段決策本階段的下泄流量，將泄流值與推得的階段平均水頭作為內部結構（水電站運行）的輸入，內部結構決策該階段電站的機組組合和流量分配，在可用流量與工作水頭確定下優化電站輸出，并將最大出力及相應的機組組合與流量分配返回至外部結構。為保證計算結果的質量，選擇動態規劃作為求解算法，針對該模型結構使用二重嵌套DP，同時進行水庫優化調度與廠內經濟運行。

對于外部結構，通常按1 h或15 min劃分調度期（1 d），將水庫優化調度轉化為多階段決策問題。假設各階段各庫的庫容均劃分為M級，則梯級水庫任意階段的狀態組合數為Mn，記t階段初組合號為pt（pt∈［1，Mn］）的梯級水庫蓄水狀態的n維向量為外部結構的DP計算流程圖如圖2所示，其具體的求解步驟如下：

（1）順序遞推計算。由調度期初逐階段向后遞推至調度期末，對于某階段t的階段末狀態組合號pt+1，其對應的蓄水狀態向量為Vt+1（pt+1），遍歷所有可能的階段初狀態Vt（pt），由式（11）求得對應于pt+1的最大累積效益及相應的最優轉換機組組合與流量分配并保存計算結果。當階段末所有組合號對應的計算均完成后，令t=t+1，進行下階段計算直至最后一個階段的計算完成。

圖2 外部結構的DP計算流程

圖3 內部結構的DP計算流程

（2）逆序遞推計算。由調度期末逐階段向前遞推至調度期初，根據所保存的推得梯級系統的調度期總發電量E及相應的各階段泄流蓄水狀態機組組合流量分配最后將計算結果輸出。

對于內部結構，本文在機組流量優化分配的過程中尋找恰當的機組組合方案，而不采用先確定投運機組再進行固定機組間流量分配的傳統分層求解方式。按機組編號將每臺可運行機組作為階段變量，從而將廠內經濟運行轉化為多階段決策問題。假設各階段的狀態（機組累積引流）均劃分為N級，記j階段初序號為qj（qj∈［1，N］）的累積引流狀態為Qj（qj）。內部結構的DP計算流程圖如圖3所示，具體的求解步驟如下：

（1）順序遞推計算。由階段1（1機組）逐階段向后遞推至階段m（m機組），對于某階段j的階段末狀態序號qj+1，其對應的累積引流為Qj+1(qj+1)，遍歷所有可能的階段初狀態Qj（qj），由式（12）求得對應于qj+1的最大累積出力及相應的最優轉換并保存計算結果。當階段末所有狀態序號對應的計算均完成后，令j=j+1，進行下階段計算直至最后一個階段的計算完成。

（2）逆序遞推計算。由階段m逐階段向前遞推至階段1，根據所保存的與推得電站優化出力Ni,t及相應的流量分配與機組組合最后將計算結果返回至外部結構。

3 求解算法性能優化

3.1 嵌套DP 算法性能分析本節從計算時間和內存占用量這兩個最為關鍵的性能參數入手，對求解精細優化模型的嵌套DP算法進行分析。首先對計算時長進行估計，算法推求精細模型的計算量集中在外部結構順序遞推的多重嵌套循環處，假設該循環任意（t，pt+1，pt）下的迭代步的運行時間均為Δτ，則整個模型的計算時間τ為：

接著對算法實現所需的內存大小進行估計。算法在計算過程中主要保存4種變量，其存儲形式如圖4所示（梯級系統調度期初末的庫容狀態固定，故相應的狀態組合號p1，pT+1∈{1}）：

本文所涉及的代碼均通過C語言進行編寫，浮點型變量通常定義為雙精度類型（double），占據8個字節，整型變量通常定義為普通類型（int），占據4個字節。故算法運行的內存占用量RAM為：

由上述分析可得，嵌套DP的時間復雜度為O（TM2n），空間復雜度為O（TMn），時空復雜度均呈指數型增長，這即是“維數災”出現的原因，同時也說明對“維數災”的處理包括減少計算時間與縮減內存占用兩方面，下面將針對這兩方面進行算法優化。

3.2 基于數據壓縮與數據庫技術的內存占用縮減以往的水庫調度研究通常僅著眼于算法求解時間方面的優化，通過算法改進減少計算量或通過CPU多核并行提升計算速度，使得模型計算時間滿足相應要求，而針對算法內存空間方面卻幾乎沒有提出相應的優化方案，且時間方面優化策略的實現往往增加了算法的空間復雜度，當問題規模較大或模擬精度較高時，計算機很容易出現內存瓶頸問題。由于實際問題的計算規模與離散精度等條件是未知的，故時空兩方面的問題均需考慮，本節將通過數據壓縮與數據庫技術，對3.1節所提的算法主要變量的存儲規模進行有效縮減，在不影響計算結果的前提下減少程序運行中的內存占用，優化后的變量存儲形式如圖5所示：

圖4 主要變量存儲形式

圖5 優化后的變量存儲形式

上述算法內存占用量優化的實現，一方面基于數據壓縮技術，在不損失有效信息的前提下完成存放變量的數組降維，另一方面基于數據庫技術，將逆序遞推計算才使用到的數據先存儲至計算機外存，以實現內存負擔向外存的轉移。同時將浮點型變量均定義為單精度類型（float），其僅占據4個字節，能在保證結果精確性的前提下減少內存占用并提升計算速度（CPU與GPU的單精度浮點計算能力均優于雙精度）。優化后程序運行的內存占用量RAM′為：

3.3 基于OpenACC的GPU并行加速近年來GPU并行加速技術發展迅速且已相對成熟，在硬件方面GPU的單精度浮點計算能力達到同期CPU的10倍以上，在軟件方面NVIDIA等公司已發布較為完善的并行編程標準CUDA和OpenACC。在水利計算上，GPU并行通常用于加速水動力學模型或分布式水文模型的求解，取得了較好效果［16］，而在梯級水庫優化調度中尚無文章明確實現GPU并行，故本節提出基于OpenACC標準的GPU并行嵌套DP，以期較以往的CPU并行進一步提升算法的計算效率。

OpenACC 是一種以編譯器導語形式實現的GPU 并行標準［17-19］，其可對復雜的C或Fortran 代碼加速。該標準中CPU稱為主機（host），GPU稱為設備（device），并行的目的是將原本由主機負責的計算稠密區域卸載至設備上多線程執行以提升運行速度，其基本步驟為：分配設備內存→將計算所需數據傳至設備內存→設備進行并行計算→將計算結果傳至主機內存→釋放設備內存。本節利用kernels構件創建并行區域，在區域內利用loop 構件實現外部順推計算的嵌套循環并行化。算法的并行流程如圖6所示，其主要步驟為：

圖6 算法并行化流程

（1）調用函數的處理。算法中水位與庫容間的值轉換、由下泄流量推求下游水位等功能已封裝為函數形式，通過調用特定的函數滿足相應的計算要求，但在并行計算區域內設備無法直接調用主機端存儲的函數，需添加routine導語以支持區域內自定義函數的調用。

具體的實現方式為：在被調函數定義與聲明的緊鄰上方添加#pragma acc routine seq，定義處的導語告知編譯器創建該函數的設備例程，聲明處的導語告知編譯器該函數可被設備端單線程直接調用。

（2）數據管理。主機內存與設備內存通常是完全分離的，這就導致主機與設備間必須進行數據傳輸，將主機內存中存儲的計算所需數據傳至設備內存以供并行計算使用，將設備內存中存儲的計算結果傳至主機內存以供逆序遞推計算使用。然而頻繁的大規模數據傳輸將消耗大量時間，為減少耗時需盡可能降低數據傳輸的量與次數，故本節利用enter data/exit data創建非結構化數據區域以減少數據傳輸耗時。在并行計算開始前，添加enter data導語為設備內存中存儲的變量分配空間，并將尾水位流量曲線、機組動力特性曲線等基礎數據傳入；在并行計算結束后，添加exit data導語將并行計算的最終結果傳回主機，之后釋放設備內存。

具體的實現方式為：在外部順推計算的嵌套循環前，添加#pragma acc enter data create（變量列表）copyin（變量列表），一次性完成主要變量的內存分配與數據傳入；在嵌套循環后，添加#pragma acc exit data copyout（變量列表）delete（變量列表），一次性完成計算結果的傳回與內存釋放。

（3）循環并行化。嵌套DP的計算量集中在外部順推計算的嵌套循環處，故應將該循環納入并行計算區域以多線程執行。由于嵌套循環最外層的階段循環具有順序依賴性，因此選擇對階段末狀態循環實現并行。通過kernels導語創建包含階段末循環的并行計算區域，并在該導語中添加present子語以告知編譯器所列變量已存在于設備內存，以此防止重復的數據操作；通過loop 導語指示編譯器在階段末循環處并行，并在該導語中添加private子語與independent子語，前者用于創建循環計算所使用的中間變量在設備各線程中的私有副本，以此防止競態錯誤，后者用于告知編譯器該層循環的迭代步間不存在數據依賴，滿足并行化要求。

上述操作具體的實現方式為：在嵌套循環的階段循環內，添加#pragma acc kernels present（變量列表），并通過{}將階段末循環及其內層循環包含在內；在階段末循環的緊鄰上方，添加#pragma acc loop private（變量列表）independent，以實現階段末循環的并行化。

（4）主機與設備的計算重疊。3.2節中為減輕計算機的內存負擔，將部分數據傳輸至外存的數據庫中存儲，當內存占用優化后的算法實現GPU并行時，則要求外部順推計算在每個階段的并行計算結束時，將該階段需存至數據庫中的計算成果先由設備內存傳輸至主機內存，再轉至主機外存中，這些數據操作將增加程序的非計算耗時。故本節先是通過update 導語將階段計算結果由設備傳至主機；接著為減少耗時，利用async子語創建異步隊列，以實現主機數據傳輸與設備并行計算的同時執行，避免計算資源的浪費，同時將數據傳輸的時間開銷隱藏于并行計算之下；最后利用wait 導語使主機線程等待設備操作的完成，防止程序出錯。

上述操作具體的實現方式為：在外部順推計算的階段循環內，添加#pragma acc update host（變量列表），將設備計算得出的上階段計算成果與傳至主機內存；在kernels導語中添加async 子語以將并行計算操作壓入異步隊列中執行；在并行區域后，將主機內存中的傳至數據庫，之后添加#pragma acc wait以等待設備操作的完成。

4 實例分析

堵河作為漢江上游南岸的一級支流，位于湖北省西北部，全長354 km，總落差500 余米。潘口與小漩水電站位于堵河干流的上游河段，開發任務以發電為主，梯級水電站的基本參數如表1所示。本文以該梯級系統的日優化調度為研究對象，選取2016年12月某日為調度期，以15 min為階段時長將該日劃分為96段。梯級系統的運行模式為“以水定電”，故系統調度期初末的庫容狀態、調度期初的機組啟停狀態及調度期各階段的入庫流量均為該日實際值，算法的模擬精度設置為：潘口各階段庫容的離散點數M1=600；小漩各階段庫容的離散點數M2=100。

實例研究所使用的軟硬件配置為：軟件方面，利用C 語言進行代碼編寫；選擇同時支持OpenMP、MPI 與OpenACC 并行標準的PGI 編譯器（PGI Professional Edition18.10）實現編譯。硬件方面，選擇Windows 10 操作系統；使用的CPU 型號為Intel Xeon E5-4640 v4，主頻2.1 GHz，物理內核數12；使用的GPU型號為NVIDIA Tesla K80，CUDA核心數4992。

表1 梯級水電站基本參數

4.1 精細模型與傳統模型對比本節分別構建梯級水電站的傳統優化模型與精細優化模型，兩者的區別在于水電站出力的計算方式不同：傳統模型分別采用恒定出力系數的出力方程和忽略機組振動區的等微增率法計算電站出力，前者標記為傳統模型1，后者標記為傳統模型2，兩者的外層水庫優化調度均通過DP求解；精細模型則將出力計算作為一個多階段優化決策問題，通過本文所提的嵌套DP求解。

表2為3個模型在計算條件相同的情況下求解所得的計劃發電量與模擬發電量，其中計劃發電量為模型求解所得的日發電效益的理論值，模擬發電量為模擬實際中調度人員按流量執行調度方案所得的日發電效益的仿真值。由表2 可以看出：（1）傳統模型1 和2 的計劃發電量分別為175.93萬kW·h與174.97萬kW·h，均高于模擬發電量的158.10萬kW·h與157.35萬kW·h，其原因在于傳統模型的調度計劃在實際執行中，原本未計量的機組啟停機耗水將消耗電站的發電水量，且為防止機組在振動區工作，須對未考慮振動區約束的機組運行方式進行調整，因此傳統模型發電量的模擬值均低于理論值；（2）精細模型的計劃發電量與模擬發電量相同，均為170.53萬kW·h，其原因在于精細模型的出力計算充分考慮了機組的啟停機耗水與振動區約束，故所制定的調度計劃貼合實際；（3）傳統模型的計劃發電量雖高于精細模型的計劃發電量，但其模擬發電量均低于精細模型的對應值，其原因在于本文所提出的精細模型將實際生產中需考慮的調度規則作為模型的約束條件，通過二重嵌套優化同時進行廠間與廠內的水量優化分配以盡可能發掘出梯級系統的發電效益，并利用嵌套DP保證求解結果的質量。

4.2 GPU 并行與CPU 并行對比4.1 節中的對比實驗證明了精細優化模型及其求解算法的有效性，之后的計算均針對精細模型并使用嵌套DP進行求解。本節首先測試3.2節中提出的基于數據壓縮與數據庫技術的內存占用縮減策略的效果，通過診斷工具檢測該策略引入前與引入后的程序運行的內存占用情況，檢測結果顯示引入后程序的內存占用量較引入前減少了99.5%，縮減幅度顯著，表明了所提空間優化策略的科學與有效性。

接下來通過不同的并行方式實現內存占用優化后算法的并行加速，并對加速效果進行比較。CPU 并行分別通過OpenMP 與MPI 實現，GPU 并行通過OpenACC 實現，并行化均是針對外部順推計算嵌套循環的階段末狀態循環層。需要注意的是階段末循環的不同迭代步的計算量可能相差很大，這意味著采用靜態調度的任務分配方式可能導致嚴重的負載不均進而影響加速效果，故本節中的3種并行方法均采用動態調度的方式以平衡負載從而保證并行效率。

基于OpenMP的CPU并行通過parallel導語創建多線程并行計算區域，通過for導語實現循環的并行化，在for導語中添加schedule（dynamic）子語以選擇動態調度的任務分配方式；基于MPI的CPU并行選擇master/worker 模式以實現動態調度下的多進程并行計算，其中master 進程負責計算任務向worker 進程的分配與管理，worker 進程負責處理計算任務并將計算結果返回至master 進程；基于OpenACC的GPU并行的流程詳見3.3節，其默認為動態調度方式。對于同一研究實例，內存占用縮減后的算法的串行版本與3種并行版本的結果如表3所示，為保證對比實驗的公平性，已盡可能挖掘出上文所列配置下的3種不同并行方法的最優結果，表3中所列的計算時長為相應并行手段下的最短時長，加速比為最優加速比。

表2 3種模型的計算結果 （單位：萬kW ·

表3 實測最優結果對比

由表3的計算結果可以看出，OpenACC 實現的GPU 并行的加速效果最優，可將算法的求解時間由串行的167.37 min壓縮至4.27 min，取得39.20的最大加速比，而OpenMP與MPI實現的CPU并行僅將計算時間縮短為19.31 min 與23.18 min，得到8.67 與7.22 的加速比。因并行計算所使用的CPU 與GPU處于同一價位，故可得在硬件成本相近時，GPU的并行計算能力遠高于CPU。通過上述分析可以得出，GPU并行較CPU并行更適合作為梯級水庫優化調度模型的加速策略。

5 結論

本文為縮減水庫優化調度的理論與實際間的差距，提出了精細至電站機組負載分配的優化調度模型，同時考慮水庫優化調度與廠內經濟運行并通過二重嵌套DP進行求解，針對“維數災”問題，提出了基于數據壓縮與數據庫技術的內存占用縮減策略以降低算法的空間復雜度，提出了基于OpenACC的GPU并行策略以減少算法的計算時間。最后以潘口小漩梯級水庫為研究實例，進行精細模型與傳統模型、優化算法內存與未優化算法內存以及GPU并行與CPU并行的對比分析，計算結果表明：本文所構建的精細模型較傳統模型更貼近電站的實際運行情況；所提出的嵌套DP能保證精細模型求解結果的質量；所提出的內存占用縮減策略能有效減輕算法對計算機的內存負擔；首次引入水庫調度領域的GPU并行較傳統的CPU并行更適合作為優化調度大規模計算的并行加速策略。