張 瀚，黃 琨，劉鳳全，王成福，周萬鵬，肖少華

（1.國家電網(wǎng)全球能源互聯(lián)網(wǎng)集團(tuán)有限公司，北京 100031；2.國網(wǎng)浙江省電力有限公司嘉興供電公司，浙江 嘉興 314033；3.山東大學(xué) 電氣工程學(xué)院，山東 濟(jì)南 250061）

0 引言

當(dāng)前各國將能源轉(zhuǎn)型、合作作為能源戰(zhàn)略的焦點(diǎn)，旨在促進(jìn)能源資源優(yōu)化配置和提高可再生能源消納水平[1]～[3]。以綜合能源系統(tǒng)（Integrated Energy System，IES）[4]概 念 為 依 托，針 對 電、氣 互 聯(lián)系統(tǒng)，國內(nèi)外學(xué)者在可靠性評估[5]、協(xié)同規(guī)劃和運(yùn)行[6]、能 量 流 分 析[7]，[8]、集 成 建 模 和 隨 機(jī) 優(yōu) 化[9]，[10]等方面取得眾多成果。

現(xiàn)有成果大多集中在區(qū)域型IES中氣、電、冷、熱等多品類能源系統(tǒng)的轉(zhuǎn)換及耦合機(jī)理分析，系統(tǒng)性的跨區(qū)域IES研究還很缺乏[11]，[12]。此外，現(xiàn)有成果較少考慮風(fēng)電隨機(jī)性對IES運(yùn)行調(diào)控的影響，而利用確定性的單一風(fēng)電預(yù)測場景數(shù)據(jù)則容易產(chǎn)生較大的計算誤差[13]。因此，迫切需要改進(jìn)調(diào)控方法，以應(yīng)對跨區(qū)型IES調(diào)控信息不完全互通而影響集中式調(diào)控。由于分布式優(yōu)化理論分解系統(tǒng)整體的優(yōu)化目標(biāo)，多個自主體彼此交換部分必要信息，同時進(jìn)行優(yōu)化調(diào)控以實(shí)現(xiàn)有效的協(xié)調(diào)合作[14]，[15]，因此，該理論與跨 區(qū)型IES協(xié)調(diào) 調(diào)控 需求不謀而合?，F(xiàn)有文獻(xiàn)中分布式優(yōu)化理論大多用于IES中不同能源類別間的分布式協(xié)調(diào)[16]，[17]，對于不同能源區(qū)域間以及二者同時考慮的研究仍然較少。在算法層面，對偶上升法的可分解性和乘子法的上界收斂屬性，在ADMM算法中實(shí)現(xiàn)了有機(jī)結(jié)合，使得該算法在算法收斂性以及解決大規(guī)模復(fù)雜問題方面具有顯著的優(yōu)勢[18]，[19]，因而成為IES分布式協(xié)調(diào)調(diào)控的重要算法。

在相互間信息不完全已知條件下，本文結(jié)合ADMM算法提出一種考慮風(fēng)電出力隨機(jī)性的CRIEGS分散協(xié)調(diào)調(diào)控策略。該策略可通過分散協(xié)調(diào)方式應(yīng)對信息不完全互通所帶來的挑戰(zhàn)。首先，在聯(lián)絡(luò)管線公共節(jié)點(diǎn)和耦合元件端口處進(jìn)行分解撕裂，構(gòu)建基于ADMM算法的分布式求解框架，并基于拉丁超立方分層采樣方法和后向場景削減技術(shù)生成典型風(fēng)電場景。進(jìn)而，以各個具有耦合關(guān)系的能源主體運(yùn)行成本的期望值最小為目標(biāo)，利用潘漢德爾修正穩(wěn)態(tài)潮流模型計算天然氣系統(tǒng)潮流分布，并考慮網(wǎng)絡(luò)中各個節(jié)點(diǎn)電壓、氣壓安全約束條件等。最后，通過模型線性化方法和變懲罰參數(shù)的交替迭代方式，以加快整體求解效率，并通過所搭建的典型算例系統(tǒng)對所提調(diào)控策略的可行性與有效性進(jìn)行驗(yàn)證。

1 CRIEGS協(xié)調(diào)調(diào)控策略

1.1 CRIEGS組成結(jié)構(gòu)及協(xié)調(diào)特點(diǎn)

本文的CRIEGS由耦合元件、天然氣和電力系統(tǒng)組成，典型結(jié)構(gòu)如圖1所示。各區(qū)域內(nèi)部的天然氣系統(tǒng)與電力系統(tǒng)的雙向耦合關(guān)系由電轉(zhuǎn)氣（Power to Gas，P2G）和 燃 氣 機(jī) 組（Gas-Fired Unit，GFU）設(shè)備構(gòu)成，區(qū)域間的天然氣、電力傳輸以管道、電力線路為紐帶。由此，跨區(qū)域能源傳輸與區(qū)域內(nèi)能源交互的內(nèi)外兩層結(jié)構(gòu)已建立成功。

圖1 CRIEGS結(jié)構(gòu)示意圖Fig.1 Schematic diagram of CRIEGS

跨區(qū)協(xié)調(diào)可分為有上級和無上級協(xié)調(diào)機(jī)構(gòu)兩種類型，分別對應(yīng)跨省及以下和跨境范圍。在我國，跨省電力交易的順利進(jìn)行基于統(tǒng)一調(diào)度、分級管理模式，在此情況下，上級協(xié)調(diào)機(jī)構(gòu)可以減少互聯(lián)區(qū)域間信息交互。而在南方電網(wǎng)與東南亞國家的電力交易中，上級協(xié)調(diào)機(jī)構(gòu)卻沒有出現(xiàn)，雙方通過簽訂合同，明確電價、電量和送電周期等信息以及聯(lián)絡(luò)線技術(shù)參數(shù)信息，在滿足輸電協(xié)議以及系統(tǒng)安全要求的前提下進(jìn)行電力輸送。此外，天然氣系統(tǒng)和電力系統(tǒng)之間也尚未形成上級協(xié)調(diào)機(jī)構(gòu)，兩者長期的獨(dú)立運(yùn)營模式屬于無上級協(xié)調(diào)類型。

1.2 分布式求解框架

ADMM算法求解問題的基本形式以及對應(yīng)的縮放型迭代過程為[21]

式中：x和z分別為優(yōu)化問題F和G的變量；A，B和c為線性等式約束系數(shù)矩陣和向量；λ，μ分別為懲罰參數(shù)和對偶向量；k為迭代次數(shù)。

在跨區(qū)層面分解互聯(lián)的2個區(qū)域，基于天然氣、電力分別通過管線和電路線路進(jìn)行傳輸?shù)膶?shí)際情況，利用兩者的節(jié)點(diǎn)運(yùn)行狀態(tài)信息反映能源傳輸情況，既可方便信息互通，又不影響私密性。因此本文利用節(jié)點(diǎn)撕裂方法對聯(lián)絡(luò)管線上的公共節(jié)點(diǎn)進(jìn)行分解，如圖2所示。

圖2 公共節(jié)點(diǎn)分解Fig.2 Decomposition of public nodes

當(dāng)聯(lián)絡(luò)管線上沒有公共節(jié)點(diǎn)時，可以以管線中點(diǎn)作為虛擬節(jié)點(diǎn)進(jìn)行分解。分解后節(jié)點(diǎn)及管線兩端運(yùn)行狀態(tài)變量應(yīng)滿足：

式中：U，θ分別為公共節(jié)點(diǎn)的電壓幅值和相角；P，Q分別為聯(lián)絡(luò)線路傳輸?shù)挠泄蜔o功功率；f為聯(lián)絡(luò)管道氣流量。

在區(qū)域內(nèi)部，以相互耦合的2個能源類別為分解對象，考慮電力、天然氣系統(tǒng)通過耦合元件實(shí)現(xiàn)相互轉(zhuǎn)換，將GFU和P2G設(shè)備分別看作電力系統(tǒng)和天然氣系統(tǒng)的組成元件，進(jìn)而在GFU和P2G設(shè)備的輸入端口處進(jìn)行分解，如圖3所示。

圖3 耦合端口分解Fig.3 Decomposition of coupling ports

分解后耦合元件輸入端口兩端的運(yùn)行狀態(tài)變量應(yīng)滿足：

式中：fELE，fGAS分別為GFU耗氣量和天然氣系統(tǒng)供氣量；PGAS，PELE分別為P2G設(shè)備耗電量和電力系統(tǒng)供電量。

式（3），（4）為 多 能 系 統(tǒng) 中 耦 合 變 量 的 線 性 等式約束，將其帶入式（2）后進(jìn)行迭代，對系統(tǒng)進(jìn)行分布式求解。迭代收斂條件滿足：

式中：μDUAL，μPRI分別為相對容忍度和絕對容忍度；s，r分別為對偶?xì)埐詈驮細(xì)埐睢?/p>

1.3 場景法

利用多個離散場景表征風(fēng)電的隨機(jī)性是一種簡單便捷的方法，主要包括場景生成和場景削減2個過程。本文基于拉丁超立方采樣（Latin Hypercube Sampling，LHS）[22]方 法 對 日 前 各 時 段 風(fēng)電預(yù)測誤差概率分布進(jìn)行分層采樣。假設(shè)在任意t時段風(fēng)電預(yù)測出力Pwt的誤差服從正態(tài)分布N（u，σ2），其 中u為 出 力 期 望，σ為 標(biāo) 準(zhǔn) 差，則 對 應(yīng)的累積概率分布函數(shù)為

將Ft（Pwt）曲線的縱軸均勻分成N等份，其中N為采樣數(shù)。構(gòu)造[0，1]上T×N階隨機(jī)數(shù)矩陣，T為時段數(shù)。任意元素 εt，n對應(yīng)t時段、第n個采樣區(qū)間的采樣點(diǎn)。風(fēng)電預(yù)測出力的采樣值由累積概率分布函數(shù)的反函數(shù)計算得到，即：

由式（7）即可得到風(fēng)電預(yù)測出力初始場景集，然后結(jié)合Gram-Schmidt序列正交化方法[24]對初始場景集進(jìn)行排列，并利用向后場景削減技術(shù)[24]對場景集進(jìn)行削減，以降低求解規(guī)模。

2 日前優(yōu)化調(diào)控

2.1 目標(biāo)函數(shù)

跨區(qū)能源合作通常是為了促進(jìn)送端系統(tǒng)可再生能源消納、保障受端系統(tǒng)能源供應(yīng)安全，最終降低系統(tǒng)整體運(yùn)行成本。因此本文以機(jī)組的燃料成本、可控設(shè)備的運(yùn)行成本以及棄風(fēng)懲罰成本作為目標(biāo)函數(shù)，在保證系統(tǒng)經(jīng)濟(jì)運(yùn)行的前提下，確保風(fēng)電的最大程度消納，各個能源區(qū)域內(nèi)部電力主體和天然氣主體目標(biāo)函數(shù)的期望E分別為

式 中：S為 場 景 數(shù)；ψs為 場 景s的 概 率；PTUi，t，PGFUi，t分別為t時段火電機(jī)組和GFU的出力；Pci，t為風(fēng)電機(jī)組棄風(fēng)量，是風(fēng)電預(yù)測出力與實(shí)際出力PWU的 差 值 ；fWEi，t，fP2Gi，t分 別 為 氣 源 和P2G設(shè) 備 供 氣量；下標(biāo)C，N分別對應(yīng)不同設(shè)備的成本系數(shù)和數(shù)量。

2.2 約束條件

本文通過添加約束條件來確保系統(tǒng)的功率平衡和安全運(yùn)行。對于CRIEGS任意節(jié)點(diǎn)k的電功率或氣流量平衡約束條件為

式 中：Pk，Qk為 節(jié) 點(diǎn)k注 入 的 有 功 和 無 功 功 率；fTRi，fPRi各自指代區(qū)域間天然氣交換量和壓縮機(jī)耗氣量；fLi，PLi為 氣 負(fù) 荷 和 有 功 電 負(fù) 荷 ；fGFUi，PP2Gi分 別為GFU和P2G設(shè)備的耗氣量和耗電量，并且分別 滿 足fGFUi=ηGFUiPGFUi，PP2Gi=ηP2GifP2Gi；wi為 標(biāo) 識 符 ，天然氣供應(yīng)端取1，接收端取-1；Q其他下標(biāo)所對應(yīng)的設(shè)備或負(fù)荷與P相同。

線性化交流潮流算法在文獻(xiàn)[25]，[26]中被應(yīng)用，本文將其應(yīng)用于計算節(jié)點(diǎn)注入功率，對節(jié)點(diǎn)注入功率方程和支路潮流方程進(jìn)行線性化處理后，可得：

式 中：bki，gki為 各 自 互 電 納 和 互 電 導(dǎo)；Bki，Gki為 支路ki上的電納和電導(dǎo)；v為有載調(diào)壓變壓器的變比，此處視為常量并取值為1。

本文使用潘漢德爾修正穩(wěn)態(tài)潮流模型表示氣壓與遠(yuǎn)距離輸氣管道的關(guān)系：

式中：Kp為管道的傳輸系數(shù)；η為管道的輸氣效率，取值一般為0.9～0.96；T為管道的輸氣溫度；ρ0，C為天然氣的相對密度和壓縮因子；L，D分別為管道的長度和內(nèi)徑。

先對式（9）左右兩邊同時進(jìn)行平方近似處理，得到：

采 用 分 段 線 性 化 方 法[26]，[27]對 式（14）等 號 兩 邊的二次項(xiàng)進(jìn)行線性化處理，以y=f2為例，分段線性函數(shù)表示為

式中：?為自變量f的分段節(jié)點(diǎn)數(shù)；fi和yi為第i個分段節(jié)點(diǎn)及其對應(yīng)的函數(shù)值；ζi為0～1的連續(xù)變量；εi為二進(jìn)制變量，確保fi左右兩端 ζi的值分別 取1和0。

式（14）等號右邊氣壓二次項(xiàng)的線性化形式與式（15）～（18）一 致，只 需 定 義 不 同 的 輔 助 變 量 即可，這里不再贅述。

機(jī)組和設(shè)備的出力上、下限約束，爬坡約束，管網(wǎng)電壓、氣壓和支路潮流安全約束等安全運(yùn)行約束條件參考文獻(xiàn)[26]。

2.3 迭代求解過程

為進(jìn)一步加快迭代的收斂速度，本文采用變懲罰參數(shù)的迭代方式[29]：

式 中 ：γDEC，γINC為 大 于1常 數(shù) ；φDEC，φINC為 大 于0常數(shù)。

本文所提算法迭代過程如圖4所示。

圖4 迭代過程Fig.4 Iterative process

3 算例分析

本文構(gòu)建如圖5所示的12bus—12node算例系統(tǒng)進(jìn)行仿真分析。

圖5 算例結(jié)構(gòu)圖Fig.5 Structure of the case

由圖5可知，2條聯(lián)絡(luò)線路、聯(lián)絡(luò)管道將2個具有相同網(wǎng)架結(jié)構(gòu)的電、氣互聯(lián)系統(tǒng)連接起來。各區(qū)域中的GFU，P2G設(shè)備將6節(jié)點(diǎn)天然氣系統(tǒng)和6節(jié)點(diǎn)電力系統(tǒng)進(jìn)行耦合。其中，天然氣系統(tǒng)為比利 時20節(jié) 點(diǎn) 系 統(tǒng)[30]的 第9～12，17和18管 段，電 力系統(tǒng)由標(biāo)準(zhǔn)IEEE6節(jié)點(diǎn)系統(tǒng)修改得到。

將風(fēng)電預(yù)測出力的15%作為標(biāo)準(zhǔn)差，取N=1000，利用LHS方法生成24×1000維的風(fēng)電預(yù)測出力場景集，然后進(jìn)行后向削減為24×10維的場景集，詳細(xì)數(shù)據(jù)及場景概率見表1，2。

表1 風(fēng)電出力情況Table1 Parameters of wind turbines kW

續(xù)表1

表2 場景概率分布Table2 Parameters of wind turbines

循環(huán)迭代關(guān)鍵參數(shù)設(shè)置情況見表3。

表3 迭代關(guān)鍵參數(shù)Table3 Key iteration parameters

3.1 分散協(xié)調(diào)調(diào)控策略可行性分析

暫不考慮風(fēng)電隨機(jī)性，利用表1數(shù)據(jù)參與系統(tǒng)分散迭代求解。仿真程序在Matlab編寫，每次迭代過程調(diào)用CPLEX進(jìn)行求解。經(jīng)過26次交替迭代過程，程序最終收斂。耦合變量殘差以及目標(biāo)函數(shù)值，如圖6，7所示。

圖6 殘差變化情況Fig.6 Changes in residuals

由圖6可知，在前8次迭代過程中，原始?xì)埐畹臄?shù)值明顯小于對偶?xì)埐?，說明按照當(dāng)前運(yùn)行狀態(tài)，部分能源主體可以達(dá)到個體運(yùn)行經(jīng)濟(jì)性最優(yōu)，但是與全局最優(yōu)狀態(tài)偏離較大，無法實(shí)現(xiàn)全局最優(yōu)，因此必須改變運(yùn)行狀態(tài)。為了不犧牲個體利益，運(yùn)行狀態(tài)并沒有發(fā)生實(shí)質(zhì)性的改變，該結(jié)論可從圖7的目標(biāo)函數(shù)變化情況得出。由圖7可知，前8次迭代完成以后目標(biāo)函數(shù)值緩慢變化，說明各個能源主體運(yùn)行狀態(tài)改變幅度較小。從第9～14次迭代過程開始，各個能源主體開始犧牲個體利益以尋找全局最優(yōu)運(yùn)行狀態(tài)。該結(jié)論也是由圖7目標(biāo)函數(shù)變化情況得出，即目標(biāo)函數(shù)值在第9～14次迭代過程中開始攀升，并逐漸趨于穩(wěn)定。在這個過程中，原始?xì)埐顢?shù)值先增后減，對偶?xì)埐罴眲∠陆?，可以表征決策沖突從出現(xiàn)到緩和的變化過程。迭代過程進(jìn)行到第15次以后，原始?xì)埐詈蛯ε細(xì)埐钜鸦揪S持在0.1以下的數(shù)值，迭代逐漸逼近最優(yōu)解。

圖7 目標(biāo)函數(shù)變化情況Fig.7 Changes in the objective function

3.2 變懲罰參數(shù)迭代方式收斂性分析

為驗(yàn)證本文所采用的變懲罰參數(shù)的迭代方法在提高交替迭代收斂性方面的作用，將迭代過程中懲罰參數(shù)不變（設(shè)置為常數(shù)1）的方法作為對比方法，將其結(jié)果與本文所提變懲罰參數(shù)方法的結(jié)果進(jìn)行對比，殘差變化情況和其他關(guān)鍵數(shù)據(jù)對比情況如圖8和表4所示。

表4 收斂情況對比Table4 Convergence comparison

圖8 殘差變化情況對比Fig.8 Comparison of residual changes

由圖8可見，在前6次迭代過程中，兩種迭代方式原始?xì)埐詈蛯ε細(xì)埐畹臄?shù)值近似相同，且殘差數(shù)值相對較大，原因在于邊界耦合變量的初始值均設(shè)置為1，與最終優(yōu)化結(jié)果的數(shù)值相差較大。在第6～12次的迭代過程中，采用懲罰參數(shù)不變的更新方式時，原始?xì)埐詈蛯ε細(xì)埐罹霈F(xiàn)了明顯的波動，而采用變懲罰參數(shù)的更新方式時殘差的波動較小。造成殘差波動的原因與增廣拉格朗日松弛算法本身特性有關(guān)。結(jié)合表4以迭代過程中殘差逆向（上升）變化為依據(jù)，統(tǒng)計的振蕩次數(shù)可見，變懲罰參數(shù)的更新方式可以有效減少殘差的波動幅度和次數(shù)，有利于殘差向著收斂條件變化，從而加快收斂速度。在第12次迭代以后，兩種方式下殘差已經(jīng)降低到較低水平，盡管仍然存在振蕩，但波動幅度較小。最終，兩種迭代方式分別在第33次和第26次迭代以后達(dá)到收斂。由于是否采用變懲罰參數(shù)的方式，都不會影響每次迭代區(qū)域內(nèi)部各自的求解過程，因此在總的求解時間方面，變懲罰參數(shù)方式所減少的求解時間與減少的迭代次數(shù)緊密相關(guān)。

3.3 場景法及分散協(xié)調(diào)調(diào)控策略有效性分析

通過上述算例結(jié)果和分析，已經(jīng)可以驗(yàn)證本文所提分散協(xié)調(diào)調(diào)控策略的可行性，以及變懲罰參數(shù)的迭代更新方式在加快收斂速度等方面的優(yōu)勢。為了進(jìn)一步分析本文采用的場景法在表征風(fēng)電出力隨機(jī)性方面的作用，以及整個分散協(xié)調(diào)調(diào)控策略計算結(jié)果的正確性，在前文的基礎(chǔ)上，利用表1，2的風(fēng)電場景數(shù)據(jù)進(jìn)行目標(biāo)函數(shù)期望的最優(yōu)值求解，并假設(shè)風(fēng)電預(yù)測出力概率最大的場景3為實(shí)際出力，分別計算實(shí)際出力和預(yù)測出力情況下的目標(biāo)函數(shù)值，并以實(shí)際出力的目標(biāo)函數(shù)值為基準(zhǔn)計算其他目標(biāo)函數(shù)及其期望的誤差，如表5所示。然后，將系統(tǒng)各個能源主體看作一個整體，即將式（8）和式（9）的目標(biāo)函數(shù)整合，將各類約束條件集中羅列，進(jìn)行統(tǒng)一集中式求解，與分散協(xié)調(diào)求解結(jié)果對比如表6所示。

表5 目標(biāo)函數(shù)對比Table5 Comparison of objective functions

表6 與集中式調(diào)控結(jié)果對比Table6 Comparison with the results of centralized regulation

由表5可見，在假定風(fēng)電預(yù)測出力概率最高的場景為實(shí)際出力時，采用單一確定性預(yù)測出力的目標(biāo)函數(shù)計算誤差較場景法的誤差要大，即給優(yōu)化結(jié)果造成的負(fù)面影響，同時考慮到場景法較其他解析方法原理和實(shí)現(xiàn)方法更加簡潔，而且需要的數(shù)據(jù)較其他抽樣方法更少，因此在跨區(qū)協(xié)調(diào)調(diào)控中具有很好的適用性。由表6可見，在求解時間方面，分散式調(diào)控的求解時間要遠(yuǎn)遠(yuǎn)大于集中式調(diào)控，這是由于分散式框架中各個能源自主體之間需要時刻交換需要的信息，以同時進(jìn)行優(yōu)化，最終迭代求得全局最優(yōu)解，即本文提出的基于ADMM的分散求解框架本質(zhì)屬于一個大型循環(huán)框架，因此，運(yùn)算的冗余程度提高。這也從側(cè)面體現(xiàn)了本文中調(diào)度模型線性化處理的必要性。

4 結(jié)論

本文提出一種考慮風(fēng)電隨機(jī)性的跨區(qū)多能系統(tǒng)分散協(xié)調(diào)調(diào)控策略，基于ADMM分布式優(yōu)化算法和耦合系統(tǒng)分解方式，實(shí)現(xiàn)了CRIEGS的分散協(xié)調(diào)調(diào)控，并以場景法應(yīng)對單一風(fēng)電出力預(yù)測場景存在較大計算誤差的問題。算例結(jié)果表明，所提分散協(xié)調(diào)調(diào)控策略能夠?qū)崿F(xiàn)多區(qū)、多能源類別間調(diào)控信息不能完全互通條件下的協(xié)調(diào)調(diào)控，并通過交替迭代以達(dá)到最優(yōu)運(yùn)行狀態(tài)。最終，通過與傳統(tǒng)的集中式調(diào)控方式進(jìn)行比對，證明了所提策略的可行性與有效性。