基于解析法的電-熱互聯(lián)綜合能源系統(tǒng)概率潮流計算

李紅，王文學(xué)，伏祥運(yùn)，何維祥，王煒，袁斌，王舒捷

(國網(wǎng)江蘇省電力有限公司連云港供電分公司，江蘇 連云港 222000)

0 引言

區(qū)域綜合能源系統(tǒng)以電力系統(tǒng)為核心，綜合考慮電、氣、熱等多種用能需求，形成具有多能協(xié)同特征的綜合系統(tǒng)[1—4]。但目前電力系統(tǒng)、熱力系統(tǒng)等普遍存在獨(dú)立規(guī)劃、獨(dú)立運(yùn)行的特點(diǎn)，各種能源網(wǎng)絡(luò)間缺乏有效交互，難以發(fā)揮多能源系統(tǒng)的耦合互補(bǔ)的優(yōu)勢[5—7]。電-熱系統(tǒng)穩(wěn)態(tài)潮流屬于確定性分析，是對綜合能源系統(tǒng)深入研究的基礎(chǔ)。綜合能源系統(tǒng)穩(wěn)態(tài)潮流主要從模型和求解方法兩方面進(jìn)行研究：在模型方面，電網(wǎng)穩(wěn)態(tài)潮流方程是由節(jié)點(diǎn)功率方程組成的非線性方程組[8]。氣網(wǎng)穩(wěn)態(tài)模型包括管道流量方程、加壓站方程和節(jié)點(diǎn)流量平衡方程[9]。熱網(wǎng)穩(wěn)態(tài)模型主要包括水力模型和熱力模型[10]。目前主要利用牛拉法求解非線性方程組，一般分為統(tǒng)一求解和順序求解[11]。

隨著新能源的大量并網(wǎng)，負(fù)荷的不確定性增大，同時能源轉(zhuǎn)化設(shè)備應(yīng)用日益廣泛，增強(qiáng)了不同能源系統(tǒng)間的耦合程度，因此綜合能源系統(tǒng)包含大量不確定性因素[12]，如各類負(fù)荷或新能源出力的波動、故障等，且不同能源網(wǎng)絡(luò)之間存在相互影響，某網(wǎng)絡(luò)的波動可能造成綜合能源系統(tǒng)潮流分布產(chǎn)生根本性的變化[13]。僅通過確定性分析難以把握綜合能源系統(tǒng)在不確定環(huán)境下的運(yùn)行狀態(tài)，概率潮流能夠充分考慮這些不確定性因素，從而更全面地反映系統(tǒng)運(yùn)行狀況。

目前，概率潮流的求解方法一般包括模擬法、解析法和近似法。模擬法中常用蒙特卡洛法，通過生成大量場景，尋找狀態(tài)量的概率信息[14]，文獻(xiàn)[15]提出一種電-熱互聯(lián)系統(tǒng)概率潮流計算方法，定量評估電轉(zhuǎn)熱技術(shù)對系統(tǒng)概率潮流影響。解析法中常用的是半不變量法[16—17]，文獻(xiàn)[18]提出一種基于多點(diǎn)線性半不變量法的電-熱互聯(lián)系統(tǒng)概率潮流計算方法，定量得到相關(guān)狀態(tài)量的概率信息。近似法中常用點(diǎn)估計法，優(yōu)點(diǎn)為無需求解輸入輸出的具體關(guān)系，但相比其他2種方法，計算精度一般較低[19]。文獻(xiàn)[20]提出一種計及相關(guān)性的電-氣-熱綜合能源系統(tǒng)概率最優(yōu)能量流計算方法，并采用基于Nataf變換的三點(diǎn)估計法對該模型進(jìn)行求解。

熱網(wǎng)方程具有強(qiáng)非線性，以往半不變量法將狀態(tài)量描述為關(guān)于電-熱負(fù)荷的線性方程，一般無法同時滿足計算復(fù)雜度、精度和計算速度的要求。文中基于解析法和近似法，提出一種輻射狀供熱網(wǎng)絡(luò)概率潮流快速計算方法，將狀態(tài)量描述為關(guān)于電-熱負(fù)荷的非線性方程，可直接獲得管道流量和節(jié)點(diǎn)溫度的概率信息，顯著提高了計算速度和求解精度。

1 熱網(wǎng)模型

1.1 熱網(wǎng)非線性模型

熱網(wǎng)模型包括水力模型和熱力模型[21]，水力模型包括流量平衡方程和管道壓降方程：

Am=mq

(1)

hf=Km|m|

(2)

式中：A為系數(shù)矩陣；m，mq分別為管道流量向量和熱負(fù)荷流量向量，kg/s；hf為環(huán)路管道壓降組成的向量，m；K為環(huán)路管道阻尼系數(shù)組成的向量；|m|為管道流量的絕對值。

熱力模型包括熱負(fù)荷方程、管道溫降方程和熱功率守恒方程。

Φ=Cpmq(Ts-To)

(3)

Tend=(Tstart-Ta)e-λL/(Cpm)+Ta

(4)

Tout∑mout=∑(minTin)

(5)

式中：Φ為熱負(fù)荷向量，MW；Cp為水比熱容，且Cp=4 182 J/(kg·℃)；Ts，To分別為節(jié)點(diǎn)供水溫度向量和回水溫度向量，℃；Tstart，Tend，Ta分別為管道首端溫度、管道末端溫度和環(huán)境溫度，℃；λ為傳熱系數(shù)，W/(m·K)；L為管道長度，m；min，mout分別為注入節(jié)點(diǎn)流量向量和流出節(jié)點(diǎn)流量向量，kg/s；Tin，Tout分別為流入節(jié)點(diǎn)溫度和節(jié)點(diǎn)混合溫度，℃。

熱網(wǎng)方程式(1)—式(5)中含有指數(shù)方程，非線性程度較高，且管道流量與節(jié)點(diǎn)溫度為乘積關(guān)系，耦合程度較高。目前常采用牛頓法求解穩(wěn)態(tài)潮流，計算較為復(fù)雜，甚至可能存在不收斂等問題。

1.2 電熱耦合元件模型

熱電聯(lián)產(chǎn)機(jī)組分為定熱電比機(jī)組和變熱電比機(jī)組2種，分別用式(6)和式(7)表示。

Cm=ΦCHP/PCHP

(6)

Cz=ΦCHP/(ηeFin-PCHP)

(7)

式中：Cm，Cz分別為機(jī)組的定熱電比和變熱電比；ΦCHP，PCHP分別為機(jī)組的熱出力和電出力；ηe，F(xiàn)in分別為機(jī)組的冷凝效率和燃料輸入速率。

1.3 輻射狀熱網(wǎng)概率潮流模型

兩節(jié)點(diǎn)供熱網(wǎng)絡(luò)如圖1所示，圖中TH，TH1，TH2分別為熱源溫度、節(jié)點(diǎn)1溫度和節(jié)點(diǎn)2溫度；mH為熱負(fù)荷流量。利用管道溫降方程式(4)與熱負(fù)荷功率方程式(3)推導(dǎo)[22]可得式(8)。

圖1 兩節(jié)點(diǎn)供熱網(wǎng)絡(luò)Fig.1 District heating network with two nodes

(8)

式中：TR為熱負(fù)荷回水溫度；Φ2為節(jié)點(diǎn)2的熱負(fù)荷。將式(8)推廣至多根管道，如式(9)所示。

(9)

多節(jié)點(diǎn)供熱網(wǎng)絡(luò)如圖2所示，對于含有分支的節(jié)點(diǎn)，此時m1≠m2≠…≠mi，其中m1，m2，…，mi-1為管道流量，mi為熱負(fù)荷流量。

圖2 多節(jié)點(diǎn)供熱網(wǎng)絡(luò)Fig.2 District heating network with multiple nodes

圖2中，引入變量n1，n2，…，ni，令n1=m1/mi，n2=m2/mi，…，ni-1=mi-1/mi，ni=mi/mi=1，代入式(10)可得mi為：

(10)

由于潮流計算前管道流量未知，令：

(11)

式中：k為熱負(fù)荷流量mi經(jīng)過的干路管道編號；Φi為節(jié)點(diǎn)i的熱負(fù)荷；Qk為管道k所流入熱能的期望值，因此n1，n2,…,nk將為定值[22]。

因此，由隨機(jī)變量的數(shù)字特征，mi的均值和方差可表示為：

(12)

(13)

式中：E(mi)為管道i中熱負(fù)荷流量mi的均值；D(mi)為管道i中熱負(fù)荷流量mi的方差。由于比熱容Cp=4 182 J/(kg·℃)，則當(dāng)式(12)中ni為定值時誤差較小，因此由熱負(fù)荷的均值和方差可得管道流量的均值和方差。將熱負(fù)荷功率方程式(3)代入式(10)，移項(xiàng)變形可得式(14)。

(14)

再由隨機(jī)變量的數(shù)字特征可得：

(15)

(16)

由式(15)、式(16)可知，得到流量倒數(shù)的均值和流量倒數(shù)的方差后，可獲得節(jié)點(diǎn)溫度THi的均值與方差。設(shè)正態(tài)分布隨機(jī)變量x=mi的均值和標(biāo)準(zhǔn)差分別為μi，σi，其概率密度函數(shù)為：

(17)

由連續(xù)型隨機(jī)變量的概率密度函數(shù)性質(zhì)可知，y=1/x=1/mi，概率密度函數(shù)為：

(18)

綜上所述，由式(12)和式(13)可得供熱網(wǎng)絡(luò)管道流量的均值與方差，由式(15)—式(18)可得節(jié)點(diǎn)溫度的均值與方差。求得熱電聯(lián)產(chǎn)(combined heating and power，CHP)熱出力概率分布后，由機(jī)組熱出力與電出力關(guān)系可得機(jī)組電出力概率分布，進(jìn)而求得電網(wǎng)中各狀態(tài)變量的概率分布。

2 輻射狀熱網(wǎng)概率潮流算例測試

23節(jié)點(diǎn)輻射狀熱網(wǎng)算例如圖3所示，CHP源溫度恒定為80 ℃，負(fù)荷節(jié)點(diǎn)回水溫度恒定為45 ℃，環(huán)境溫度Ta為10 ℃，測試采用Matlab編程，仿真環(huán)境為1.9 GHz CPU，4.0 GB RAM的PC機(jī)。

圖3 電-熱互聯(lián)綜合能源系統(tǒng)Fig.3 The diagram of combing electricity and heat networks

2.1 典型狀態(tài)量的均值與方差

若每根管道長度設(shè)置為500 m，所有負(fù)荷節(jié)點(diǎn)熱負(fù)荷均設(shè)置為0.5 MW，熱負(fù)荷波動均在±10%以內(nèi)。文中方法和蒙特卡洛法所得流量均值、流量標(biāo)準(zhǔn)差分別為μm，σm，μm,mcs，σm,mcs；流量均值誤差百分?jǐn)?shù)、流量標(biāo)準(zhǔn)差誤差分別為δμ,m，δσ,m。因此有:

δμ,m=(|μm-μm,mcs|/μm,mcs)×100%

(19)

δσ,m=|σm-σm,mcs|

(20)

溫度均值、溫度標(biāo)準(zhǔn)差分別為μT，σT，μT,mcs，σT,mcs；溫度均值誤差百分?jǐn)?shù)、溫度標(biāo)準(zhǔn)差誤差分別為δμ,T，δσ,T。因此有：

δμ,T=(|μT-μT,mcs|/μT,mcs)×100%

(21)

δσ,T=|σT-σT,mcs|

(22)

其中蒙特卡洛結(jié)果是采用熱網(wǎng)非線性模型計算所得，采樣20 000次。此時計算出的流量和溫度的均值與標(biāo)準(zhǔn)差如表1和表2所示，文中方法耗時0.019 6 s，蒙特卡洛法耗時1 720.105 4 s。

表1 所有管道流量均值和標(biāo)準(zhǔn)差的比較Table 1 Comparison of mean and standard deviation of all pipe flow

表2 所有節(jié)點(diǎn)溫度均值和標(biāo)準(zhǔn)差的比較Table 2 Comparison of mean and standard deviation of all node temperatures

由表1可知，管道流量均值誤差百分?jǐn)?shù)最大值為0.024 4%，平均值為0.008 0%，管道流量標(biāo)準(zhǔn)差誤差最大值為0.003 1 kg/s，平均值為0.001 3 kg/s。由表2可知節(jié)點(diǎn)溫度均值誤差百分?jǐn)?shù)最大值為0.001 5%，平均值為0.000 5%，節(jié)點(diǎn)溫度標(biāo)準(zhǔn)差誤差最大值為0.001 1 ℃，平均值為0.000 2 ℃。可見文中方法具有較高的計算精度和計算速度。

由表1可知，文中方法和蒙特卡洛法所得管道1的流量均值分別為42.266 9 kg/s和42.264 5 kg/s，2種方法求得的熱源節(jié)點(diǎn)回水溫度均值均為44.712 2 ℃，由式(3)計算可得熱源的熱出力相差354.176 9 W，若式(6)中定熱電比取1.3，則2種方法求得的熱出力誤差引起的電壓和相角誤差分別如圖4和圖5所示。

圖4 熱出力誤差引起的電網(wǎng)電壓誤差Fig.4 Error of voltage amplitude caused by thermal output

圖5 熱出力誤差引起的電網(wǎng)相角誤差Fig.5 Error of voltage phase angle caused by thermal output

由圖4和圖5可知，文中方法和蒙特卡洛法引起的電網(wǎng)電壓誤差和相角誤差極小，表明文中方法產(chǎn)生的誤差對電網(wǎng)影響極小，驗(yàn)證了該方法的有效性與合理性。

2.2 控制變量法

(1)設(shè)置每個負(fù)荷節(jié)點(diǎn)熱負(fù)荷均為1 MW，所有熱負(fù)荷波動量為±10%且保持不變，逐漸同時增大每根管道長度，節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差見圖6(節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差誤差最大)。

圖6 節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差隨管道長度變化情況Fig.6 The temperature standard deviation of node 19 varying with the length of the pipeline

由圖6可知，節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差與管道長度近似呈線性正相關(guān)，且δσ，T小于0.05%。

(2)設(shè)置每根管道長度為500 m，熱負(fù)荷波動量為±10%且保持不變，逐漸同時增大每個負(fù)荷節(jié)點(diǎn)熱負(fù)荷值，節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差如圖7所示。

圖7 節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差隨熱負(fù)荷值變化情況Fig.7 The temperature standard deviation of node 19 varying with the heat load value

由圖7可知，節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差與熱負(fù)荷值呈負(fù)相關(guān)，且曲線的切線逐漸平緩，熱負(fù)荷值對節(jié)點(diǎn)溫度標(biāo)準(zhǔn)差的影響相對較大。同時，隨著熱負(fù)荷值的增大，文中方法的誤差逐漸縮小，當(dāng)熱負(fù)荷值為1 MW時，節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差誤差小于0.02%。

(3)設(shè)置管道長度為500 m，熱負(fù)荷為1 MW且保持不變，逐漸同時增大熱負(fù)荷波動量，節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差如圖8所示。

圖8 節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差隨熱負(fù)荷波動量變化情況Fig.8 The temperature standard deviation of node 19 varying with the fluctuation of heat load

由圖8可知，節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差與熱負(fù)荷波動量近似呈線性正相關(guān)，且δσ，T小于0.1%。

2.3 典型狀態(tài)量的三維關(guān)系圖

設(shè)置熱負(fù)荷波動量為±10%，當(dāng)改變管道長度和熱負(fù)荷值時，采用蒙特卡洛法所得管道1流量標(biāo)準(zhǔn)差和節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差分別見圖9和圖10。

圖9 管道1流量標(biāo)準(zhǔn)差與管道長度和熱負(fù)荷的關(guān)系Fig.9 The relationship among the flow stand-ard deviation of pipe 1,the pipe length and heat load value

圖10 節(jié)點(diǎn)19溫度標(biāo)準(zhǔn)差與管道長度和熱負(fù)荷關(guān)系Fig.10 The relationship among the temperature standard deviation of node 19,the pipe length and heat load value

由圖9可知，管道長度與流量標(biāo)準(zhǔn)差沒有呈現(xiàn)相關(guān)性，熱負(fù)荷值與流量標(biāo)準(zhǔn)差近似呈線性正相關(guān)。由圖10可知，管道長度與節(jié)點(diǎn)溫度標(biāo)準(zhǔn)差呈正相關(guān)，而熱負(fù)荷值與節(jié)點(diǎn)溫度標(biāo)準(zhǔn)差呈負(fù)相關(guān)。

3 結(jié)論

文中提出了一種輻射狀供熱網(wǎng)絡(luò)概率潮流快速計算方法，可同時獲得熱力網(wǎng)的穩(wěn)態(tài)潮流與概率潮流，無需迭代，計算準(zhǔn)確度高且不存在收斂性問題，結(jié)論如下：

(1)管道長度與流量標(biāo)準(zhǔn)差沒有呈現(xiàn)相關(guān)性，管道流量標(biāo)準(zhǔn)差與熱負(fù)荷值近似呈線性正相關(guān)。

(2)節(jié)點(diǎn)溫度標(biāo)準(zhǔn)差與熱負(fù)荷波動量近似呈線性正相關(guān)，節(jié)點(diǎn)溫度標(biāo)準(zhǔn)差與管道長度近似呈線性正相關(guān)，節(jié)點(diǎn)溫度標(biāo)準(zhǔn)差與熱負(fù)荷值呈負(fù)相關(guān)。

(3)文中方法引起的熱源熱出力誤差對電網(wǎng)電壓和相角的影響極小。

文中方法在保證計算精度的同時，極大地提高了計算速度，可為綜合能源系統(tǒng)不確定性分析提供參考。下一步將研究大規(guī)模可再生能源并網(wǎng)對電-熱互聯(lián)綜合能源系統(tǒng)概率潮流的影響。

本文得到國網(wǎng)江蘇省電力有限公司連云港供電分公司伏祥運(yùn)勞模創(chuàng)新工作室2021年研究項(xiàng)目“電-熱互聯(lián)綜合能源系統(tǒng)快速潮流計算方法與優(yōu)化調(diào)度方法技術(shù)的研究與應(yīng)用”(B110D02188D7)資助，謹(jǐn)此致謝！