量調節方式下區域熱電系統的聯合最優潮流

韓 赫，張沛超，杜 煒，劉學智，徐博強，姚 程

（1.電力傳輸與功率變換控制教育部重點實驗室（上海交通大學）,上海市200240；2.國電南瑞科技股份有限公司,江蘇省南京市211106；3.國網天津市電力公司電力科學研究院,天津市300384）

0 引言

區域供熱網（district heating network,DHN）在現代社會中具有重要作用。DHN 可以與其他能源形式構成綜合能源系統［1-2］,例如通過熱電聯供（combined heat and power,CHP）機組構成區域熱電系統（district heat-electric system,DHES）,從而實現多能耦合,有效提高能源利用效率［3］。

DHN 有質調節和量調節這2 種基本調節方式。前者采用“變溫恒流”,而后者采用“恒溫變流”［4］。由于質調節的水力工況穩定、易于控制,現有研究多集中于質調節方式。例如,文獻［5］假設管道熱水流量恒定,通過改變CHP 機組供熱蒸汽流量來調節管網供熱溫度；文獻［6］在質調節模式下提出考慮機組調峰主動性的電熱協調調度方法；文獻［7-8］基于廣義電路分析理論,建立多能源系統的支路及網絡模型,描述熱網的動態特性。然而,質調節存在熱傳輸時延且單純使用時易造成能源浪費［9］。對此,文獻［10］提出應在供暖期不同階段采用不同的調節方式。例如,在供暖初、末期負荷較小時,系統以最小設計流量按質調節運行；在供暖中期,當供熱溫度達到設計值后按量調節運行?？梢?實際運行中量調節與質調節方式應構成有效互補。但是,針對質調節的研究以恒流作為基本假設,難以應用于量調節。為此,本文專門研究量調節方式。

量調節方式下由于熱力工況穩定,無須考慮熱傳輸時延問題。但因流量可變,其水力模型難以線性化處理。同時,其熱力模型也為非線性,并與水力模型存在耦合。故采用量調節的文獻均使用非線性、非凸模型［11-12］。但一般的非凸優化問題缺乏穩定、高效的解法,且難以得到全局最優解。為此,在保證足夠精度的前提下,本文研究如何將上述一般的非凸優化問題轉換為易于求解的問題。本文的主要貢獻如下。

1）提出水力支路壓降的二階錐松弛方法。水路管道壓降與流量呈平方關系,文獻［13］假設水力工況為層流,從而實現壓降方程線性化。但DHN 為避免流量過小引發水力失調,實際工況常處于紊流阻力平方區［14］,線性近似誤差較大。本文提出基于凸包絡的壓降方程松弛方法,并證明該方法為精確松弛。

2）提出熱力支路的等效熱損方法。針對DHN水力、熱力耦合的非線性因素,文獻［13,15］先采用定流量求解熱力支路模型,再迭代更新流量直至算法收斂。本文根據量調節特點,提出熱網等效熱損模型,具有較高精度且無須交替求解水、熱模型。

3）在上述基礎上,建立DHES 聯合最優潮流的凸凹規劃（convex-concave programming,CCP）模型,并將其轉換為二階錐規劃（second-order cone programming,SOCP）問題進行序貫求解,顯著提高了求解效率。

1 DHN 模型及其凸化方法

DHN 由空間同構的供、回水網組成［11］,如圖1所示。本文DHN 采用量調節方式,熱源溫度Tsrc和熱負荷出口溫度To恒定。Ts1,Ts2,Ts3分別表示節點1,2,3 的供熱溫度,熱水流經熱負荷后,溫度從供熱溫度降為To。Tr1,Tr2,Tr3分別表示回水網中節點1,2,3 的回熱溫度,由于回水網溫度低、熱損小,可假設各節點回熱溫度近似一致且等于To［11,13］,這樣研究中僅需關注供水網。本章研究如何對DHN 中的非線性因素做松弛和近似變換。

圖1 DHN 結構Fig.1 Structure of DHN

1.1 支路模型

1.1.1 水力支路模型松弛方法

管道壓降可用Darcy-Weisbach 公式描述［13］,即

式中：h和mb分別為管道壓降和流量（即質量流率）；L和D分別為管道長度和內徑；ρ為熱水密度；g為重力加速度；λ為Darcy 摩擦系數；Λ為阻力系數。

由式（1）可知,h與mb符號一致。在紊流阻力平方區下,λ由Prandtl-Karman 公式確定,為只取決于管道自身特性、不隨流量變化的常數［14］,即

式中：ε為管道絕對粗糙度。

式（1）為非仿射等式約束。為保證水力工況的相對穩定,本文假設管道流向在一個調度周期內保持不變,且可以在日前獲知［16］。設該流向為正方向,然后將式（1）松弛為二階錐不等式約束,即

上述松弛擴大了問題的可行域,對此將在1.2 節解決。

1.1.2 熱力支路模型近似方法

管道溫度下降方程描述了熱水與外界環境的熱交換過程,表現為熱水流經管道后溫度下降,即

式中：Tstart,Tend,Ta分別為管道始端溫度、末端溫度及環境溫度；β為管道熱傳導系數；σ為熱水比熱容。

進一步可得管道熱損Φloss為：

流量和溫度的乘積項表現了DHN 中水、熱網的非線性耦合關系。本文提出熱網等效熱損模型,方法如下。

首先,通常情況下有βL?σmb,因此可將式（4）中指數項做一階泰勒展開［16-18］,即

將其代入式（5）,可得：

其次,由于量調節方式下熱網熱力工況穩定,且隨著DHN 向第4 代發展［19］,熱損較低,節點供熱溫度Ts一般僅比熱源溫度Tsrc低1～2 ℃［4,15］,遠小于供給熱負荷的溫降。因此,本文將各節點供熱溫度均近似為Tsrc,從而式（7）可改寫為式（8）,式（8）實際上代表管道對環境的穩定熱傳導。

最后,將上述熱損等效疊加至管道末端的熱負荷上,如圖2 所示。經過以上處理,在量調節方式下可實現DHN 中水、熱的解耦分析。需要指出的是,上述熱力支路近似模型也會給管道流量帶來一定的計算誤差,詳細分析見附錄A。

圖2 熱力支路等效熱損模型Fig.2 Equivalent heat loss model of thermal branch

1.2 拓撲約束模型

1.2.1 水力拓撲約束

類似基爾霍夫第一定律,各節點應滿足如下流量平衡方程。

式 中：mni為 節 點i的 凈 注 出 流 量；mbk為 管 道k的 流量；aik為節點關聯矩陣A中的元素,定義為

類似基爾霍夫第二定律,環路壓降平衡方程要求任意環路的管道壓降代數和應等于0,即

式中：hk為管道k的壓降；bk為環路關聯矩陣B中的元素,定義為

為收緊式（3）造成的松弛間隙,本文在目標函數中添加懲罰項fh,通過最小化fh可實現零間隙,即

式中：Ω為懲罰系數；Λk為支路k的阻力系數。

收緊管道壓降方程松弛如圖3 所示,由圖3 可見,式（3）的實質是用凸包絡（圖中藍色區域）代替原可行域以達到凸化目的；而懲罰項又重新收緊了可行域,實現了精確松弛。詳細證明見附錄B。

圖3 收緊管道壓降方程松弛Fig.3 Tightening relaxation of pipe pressuredrop equation

1.2.2 熱力拓撲約束

在管道的連結節點處,根據焓守恒原理可以導出如下方程。

式中：min和mout分別為流入、流出混合節點的管道流量；Tin為流入節點管道的末端溫度；Tout為流出節點管道的始端溫度。

但在量調節方式下,由于各節點供熱溫度近似一致,式（14）可以用式（9）的流量平衡方程代替。

2 DHES 聯合潮流模型

2.1 配電網模型

交流潮流模型為典型的非線性方程組。節點i的功率注入模型為：

式中：Pni和Qni分別為節點i的凈注入有功和無功功率；Vi,Vj,θij分別為節點i和j的節點電壓和相角差；Gij和Bij分別為節點導納陣第i行、第j列元素的實部和虛部。

對于輻射型配電網,可將節點注入模型式（15）—式（16）改寫成DistFlow 支路潮流模型［20］。同時,可仿照熱網,對線損做相同等效處理,可得到如下潮流方程。

式中：Pbk和Qbk分別為支路k的有功和無功功率,其中支路k由節點i流向節點j；J為節點j的流出支路k'的集合；Rk和Xk分別為支路k的電阻和電抗；lk為支路k的電流平方值；ui為節點i的電壓平方值。

將式（20）松弛為具有二階錐形式的不等式約束式（21）,從而建立電網潮流方程的凸優化模型［21］。

文獻［21］證明了對于輻射型配電網,上述支路潮流松弛為精確松弛。

2.2 DHES 聯合潮流模型的矩陣形式

至此,可將DHES 聯合潮流模型統一表達為矩陣形式。定義DHES 的聯合節點關聯矩陣為：

式中：AE和AH分別為配電網和DHN 的節點關聯矩陣?？梢?矩陣A仍保留了強稀疏性。

定義Apos為矩陣A的非負部分,表征節點與流入節點的支路之間的對應關系,其元素apos,ik為：

DHES 聯合潮流模型統一將支路的功率損失等效到支路末端的負荷上,這樣,DHES 的系統功率平衡約束可以統一寫成：

式中：Pb,Pn,Ploss,Psrc,Pload分別為有功功率列向量P的支路潮流、節點潮流、支路損耗、節點電源、節點負荷列向量；無功功率列向量Q及熱功率列向量Φ以此類推。

各類支路損耗分別為：

式中：R和X分別為支路電阻和電抗的列向量；l為支路電流的平方列向量；L和β分別為支路管道長度和管道熱傳導系數列向量。向量中無該類參數及變量的元素取為0；“°”運算符表示2 個向量按元素相乘。

DHN 環路壓降平衡約束為：

式中：h,Λ,mb分別為支路管道壓降、支路阻力系數、支路管道流量的向量形式。

其余電、熱潮流約束也寫為矩陣形式,即

式中：u為節點電壓平方值列向量。

最后,各變量還需滿足如下運行約束。

式中：Pb,max,Qb,max,Imax分別為支路有功潮流、無功潮流、電 流 上 限 列 向 量；Vmax,Vmin和mb,max,mb,min分 別為節點電壓和支路流量的上、下限列向量。

3 DHES 聯合最優潮流

3.1 CCP 模型

DHES 優化的目標函數為：

式中：fP,t,fΦ,t,fchp,t分別為t時段的購電成本、購熱成本和CHP 成本；fh,t為式（13）中t時 段的壓降懲罰成本。

購電與購熱成本皆為線性函數,即

式中：CP,t和Pbuy,t分別為t時段的電價和購電量；CΦ,t和Φbuy,t分別為t時段的熱價和購熱量。

CHP 運行成本包括發電、產熱和耦合成本項［16］,可表示為：

式中：α0至α5為成本系數；Pchp,t和Φchp,t分別為t時段CHP 機組電、熱出力。

考慮如下約束條件。

1）熱電聯合潮流約束

熱電聯合潮流約束見式（24）—式（35）。

2）設備運行約束

CHP 設備運行約束見文獻［22］。儲能（包括電儲能和熱儲能）運行約束見文獻［23］,針對儲能的非線性充放互斥約束,該文獻提出了一種線性松弛方法。

3）購能約束

式中：Pbuy,max和Φbuy,max分別為購電和購熱上限。

分析上述優化問題（記為問題P1）可知,其約束條件皆為仿射約束或二階錐約束；在目標函數式（36）中,懲罰項fh,t為凹函數,其余為凸函數。這樣,基于第1 章和第2 章的工作,本文將DHES 聯合最優潮流轉換為一類特殊優化問題——CCP 問題［24］。雖其仍為非凸問題,但卻有可靠的求解方法,并已在電力系統中獲得應用［20］。

3.2 序貫SOCP 求解方法

本文將上述CCP 問題轉換為可序貫求解的凸優化子問題,具體步驟如下。

步驟1：首先忽略式（36）中的懲罰項fh,t,得到凸優化子問題（記為問題P2）,求解問題P2并獲得各支路流量的初值。問題P2保留了問題P1的全部約束條件,則所求得初值必為問題P1的可行解,這能確保問題P1在后續迭代過程中收斂［24］。

步驟2：考慮懲罰項fh,t,但將其在當前流量值Mb處做一階泰勒展開,得到凸優化子問題并求解。該凸優化子問題的目標函數（記為問題P3）如式（42）所示。

式中：hk,t和mbk,t分別為t時段管道k的壓降和流量；Mbk,t為Mb中 元 素。

在這一步中,僅對懲罰項做了線性處理,盡量保留了問題P1的信息,這有助于快速收斂。

步驟3：根據優化結果更新流量值,重復步驟2,直至目標函數值改變量小于容差（本文設為10-3）。至此得到問題P1的最優解,該最優解已重新收緊水力支路壓降的可行域。

但在凸優化子問題P2和P3中,目標函數中仍包含CHP 的二次成本函數項。本文引入輔助變量δt,并增加二階錐不等式約束,將最小化fchp,t轉變為如下等價問題。

這樣,子問題P2和P3被進一步變換為SOCP 問題,變換過程詳見附錄C。最終,問題P1被轉換為序貫SOCP 問題。相比一般凸優化問題,SOCP 問題具有非常高效的求解方法。

4 算例驗證

4.1 算例設置

DHES 算例的網絡拓撲如附錄D 圖D1 所示,包含IEEE 33 節 點 配 電 網 和 巴 厘 島32 節 點DHN［11］。熱電耦合設備為背壓式CHP 機組和抽凝式CHP 機組,其熱電可行域參數見文獻［22］,電熱出力范圍均為0.2～1 MW,成本系數見附錄D 表D1。

電網參數為：節點1 購電功率上限為3 MW；節點25 連接風電場；節點6 配置電儲能,參數見附錄D表D2；節點3,8,13,28 均配置無功補償裝置,最大容量均為0.6 Mvar；線路潮流有功和無功功率的上限分別為3 MW 和2 Mvar；設置節點1 為平衡節點,基準電壓為12.66 kV,其余節點電壓標幺值范圍為0.95～1.05；負荷及主網電價見圖D2,風電預測曲線見圖D3。

熱網參數為：節點1 向一級熱網購熱功率的上限為1 MW；節點5 連接太陽能熱廠；節點14 配置熱儲能,參數見附錄D 表D2；管道流量下限取為2.2 倍內徑,以保證雷諾系數不低于104（紊流區）,上限取為6 kg/s,其余參數見文獻［11］；熱源溫度為設計溫度70 ℃,熱負荷出口溫度為30 ℃；熱負荷及主網熱價見圖D2,環境溫度及光熱預測曲線見圖D3。

總調度周期為24 h,單個滾動優化周期包含6 個時段,每時段為1 h。

4.2 模型精度分析

本節驗證本文第3 章所建立的CCP 模型的求解精度。對比模型為未做松弛與近似變換的非凸精確模型,上述2 類模型的區別見附錄C 表C1。本節設置如下3 個算例：算例a 采用CCP 模型；算例b 采用非凸精確模型,取零初值；算例c 同樣采用非凸精確模型,但以算例a 的結果作為初值。由于初值選取恰當,認為算例c 求得了全局最優解。在求解方面,算例a 采用本文3.2 節的序貫SOCP 方法,SOCP子 問 題 的 求 解 器 為SeDuMi［25］；而 算 例b 和c 均 使 用內點法求解。

1）水路壓降松弛方法驗證

CCP 模型中環路管道計算壓降的代數和（取正值）如附錄E 圖E1 所示。計算壓降特指由流量mb優化結果所計算的壓降值Λm2b。由圖E1 可見,各時段環路壓降代數和均在10-8～10-6量級之間,近似為0,證明本文對水路壓降的松弛方法能有效滿足環路壓降平衡方程。表E1 中還給出了CCP 模型與非凸精確模型關于流量與壓降的詳細對比結果。

2）熱網等效熱損模型驗證

算例a 與c 的管網總熱損對比如圖4 所示。由圖4 可知,由于CCP 模型將節點供熱溫度近似為熱源溫度,導致熱損始終偏高。進一步對照附錄D 圖D2 中熱負荷曲線及圖D3 中環境溫度趨勢可以發現,外界環境溫度越高,DHN 熱損越??；熱負荷越小,CCP 模型的熱損相對于精確熱損的誤差越大,但最大誤差不超過2%?？紤]到總熱損僅占總負荷的3%～5%,上述絕對誤差很小。而附錄E 表E1 已表明水力計算結果誤差很小。綜上可見,本文雖然對水力、熱力做解耦建模并獨立分析,但等效保留了管道熱水向外散熱的損耗影響,計算結果足以滿足工程需要。

圖4 各時段總熱損對比Fig.4 Comparison of total heat loss in each hour

為了驗證熱網在量調節方式下熱力工況穩定的特點,本文基于非凸精確模型,給出了調度期間各節點溫度的實際變化與分布情況,見附錄E 圖E2。

3）優化結果驗證

各時段總成本對比如圖5 所示。比較算例a 和c可見,兩者求得的系統總成本幾乎一致,絕對誤差均保持在10-2量級,最大和最小相對誤差分別出現在時段5（0.046%）和時段20（0.004%）。這說明本文方法保留了很高的精度。設備出力與成本詳細對比結果見附錄E 表E2。

圖5 各時段總成本對比Fig.5 Comparison of total cost in each hour

進一步對比算例b 和c,發現部分時段算例b 的總成本明顯偏高。這說明雖然算例b 和c 均采用精確模型,但由于所建立的優化問題非凸,其結果受初值影響較大,容易陷入局部最優。

4.3 模型計算效率分析

本節基于算例a 和c 進一步驗證模型的計算效率。統計各算例在每個調度周期的平均求解時間發現,算例c（非凸精確模型+內點法）的平均求解時間 為49.96 s,而 算 例a（CCP 模 型+序 貫SOCP 方法）的平均求解時間為6.21 s。可見,本文方法在保證精度的同時,還顯著提高了計算效率,這非常有利于大規模DHES 的優化計算。

附錄E 圖E3 為部分時段序貫SOCP 的收斂過程。由收斂曲線可見,求解初值時（問題P2）因不考慮壓降懲罰項,目標函數值誤差很大,說明經松弛后壓降求解結果遠遠偏離與流量的二次關系曲線；當進入序貫SOCP 迭代后（問題P3）,懲罰項迫使壓降松弛間隙逐漸逼近0,目標函數值經不超過10 次的迭代后快速收斂。

5 結語

量調節是DHN 的一種基本調節方式,本文針對其特點提出了DHES 的聯合最優潮流模型及其高效解法,主要結論如下。

1）量調節方式下難以對水力模型做線性化處理。可利用凸包絡對熱網的水力壓降方程進行松弛,并通過設計懲罰函數來收緊松弛間隙。

2）量調節方式具有熱力工況穩定的特點,據此可實現穩態水力、熱力模型的解耦,結果表明該近似具有很高精度,并能等效計及管網熱損。

3）DHES 的聯合最優潮流模型可被表征為特殊的非凸優化問題,CCP 問題進而轉換為序貫SCOP問題。結果表明該求解方法能顯著提高計算效率,有利于大規模DHES 的運行優化。

本文和已有文獻均假定在一個調度周期內僅有一種調節方式。下一步將研究量調節與質調節的最優協調與方式切換問題。

附錄見本刊網絡版（http：//www.aeps-info.com/aeps/ch/index.aspx），掃英文摘要后二維碼可以閱讀網絡全文。