基于異質(zhì)能源統(tǒng)一建模方法的天然氣網(wǎng)安全域穩(wěn)態(tài)模型

宋晨輝，肖 峻，秋澤楷，焦 衡，屈玉清

（智能電網(wǎng)教育部重點(diǎn)實(shí)驗(yàn)室（天津大學(xué)），天津市 300072）

0 引言

綜合能源系統(tǒng)（integrated energy system，IES）是實(shí)現(xiàn)多能互補(bǔ)與碳中和的關(guān)鍵載體［1］，中國(guó)能源局已將綜合能源服務(wù)納入發(fā)展規(guī)劃［2］。但多能互聯(lián)也帶來(lái)了新的安全問(wèn)題，即擾動(dòng)易在異質(zhì)能源系統(tǒng)傳遞，引發(fā)大范圍連鎖故障［3-4］。IES 安全性已成為國(guó)內(nèi)外研究焦點(diǎn)。

IES 安全性研究主要基于“逐點(diǎn)法”［3，5-10］，包括穩(wěn)態(tài)［3，5-7］與動(dòng)態(tài)［8-9］兩方面。針對(duì)穩(wěn)態(tài)，文獻(xiàn)［3］提出了多能流安全分析的概念；文獻(xiàn)［5］提出一種靈敏度分析方法，可分析注入電功率對(duì)燃?xì)鈮毫Φ挠绊懀晃墨I(xiàn)［6］提出一種考慮管道N-1 的多能流安全分析方法；文獻(xiàn)［7］研究了不同工況下耦合元件對(duì)IES 安全性的影響。文獻(xiàn)［8-9］進(jìn)一步考慮了動(dòng)態(tài)，文獻(xiàn)［8］建立了考慮動(dòng)態(tài)的多能流模型；文獻(xiàn)［9］進(jìn)行IES 狀態(tài)估計(jì)時(shí)，考慮了天然氣暫態(tài)前后時(shí)間斷面的相關(guān)性。

相較逐點(diǎn)法，電力系統(tǒng)安全域方法具有分析效率高、安全信息豐富的優(yōu)勢(shì)［10］，逐漸被應(yīng)用于IES 及天然氣網(wǎng)的穩(wěn)態(tài)分析［11-16］。文獻(xiàn)［11］提出了IES 安全域的概念，并提出一種高精度的邊界擬合方法。文獻(xiàn)［12］提出了基于凸包的IES 魯棒安全域，并分析了忽略管存的穩(wěn)態(tài)域的可用性。文獻(xiàn)［13］在區(qū)域IES 安全域建模時(shí)考慮了N-1。對(duì)于天然氣網(wǎng)，文獻(xiàn)［14］提出了天然氣網(wǎng)安全域的概念、模型與安全控制方法；文獻(xiàn)［15］提出了天然氣網(wǎng)在給定注入流量下的壓力安全域?，F(xiàn)有研究表明安全域?qū)τ谔烊粴饩W(wǎng)的在線安全評(píng)估具有重要作用，能夠判斷供需關(guān)系下運(yùn)行方案的可行性，為調(diào)度提供決策［14-15］。

但無(wú)論是逐點(diǎn)法，還是域方法，均以異質(zhì)能源網(wǎng)各自固有模型為基礎(chǔ):如電力網(wǎng)采用潮流方程、天然氣網(wǎng)采用管道壓降方程［14］、熱力網(wǎng)采用水力與熱力方程［7］。這些模型的數(shù)學(xué)差異很大，缺乏統(tǒng)一框架，為異質(zhì)能源的綜合分析帶來(lái)了壁壘［16-17］。

異質(zhì)能源網(wǎng)本質(zhì)上都是能量網(wǎng)絡(luò)，可以從統(tǒng)一的視角進(jìn)行建模分析［18］，文獻(xiàn)［16-17，19-23］對(duì)此展開了研究。文獻(xiàn)［19-20］最早將電路思想用于異質(zhì)能源網(wǎng)分析，建立了多能源網(wǎng)絡(luò)的廣義電路理論［16，21］。文獻(xiàn)［17，22-23］提出統(tǒng)一能路理論，建立了異質(zhì)能源的能路模型，研究已涉及天然氣網(wǎng)的氣路［17］與狀態(tài)估計(jì)［23］、潮流計(jì)算［22］等。

廣義電路理論與統(tǒng)一能路理論的思想都是將異質(zhì)能源系統(tǒng)的穩(wěn)態(tài)與動(dòng)態(tài)分析納入統(tǒng)一框架，本文首次將該思想引入天然氣網(wǎng)安全分析，提出了天然氣網(wǎng)的安全域模型。較已有模型［14］，本文模型在數(shù)學(xué)形式上實(shí)現(xiàn)了天然氣網(wǎng)和電力網(wǎng)在安全域上的統(tǒng)一，為后續(xù)建立IES 安全域通用化模型奠定了基礎(chǔ)。同時(shí)，顯著提升了天然氣網(wǎng)安全域的計(jì)算效率，發(fā)掘了異質(zhì)能源統(tǒng)一建模方法在天然氣網(wǎng)安全分析中的應(yīng)用價(jià)值。

1 天然氣網(wǎng)安全域的基礎(chǔ)概念

1.1 工作點(diǎn)

電力系統(tǒng)采用工作點(diǎn)［10］描述安全性。參考電力系統(tǒng)定義天然氣管網(wǎng)工作點(diǎn):能描述管網(wǎng)運(yùn)行狀態(tài)安全性的最小變量集合［14］。工作點(diǎn)選取運(yùn)營(yíng)商關(guān)心的流量作為變量［14］:

式中:W為工作點(diǎn)向量；Gl，i和Gs，j分別為負(fù)荷i、氣源j的流量；m和n分別為負(fù)荷和氣源的節(jié)點(diǎn)個(gè)數(shù)；Ga為節(jié)點(diǎn)a的流量，其中a=1，2，…，m+n。

1.2 安全性

天然氣管網(wǎng)的運(yùn)行安全性［14］定義為:對(duì)于某工作點(diǎn)，其所有狀態(tài)量是否滿足運(yùn)行約束。若滿足，則運(yùn)行安全，該點(diǎn)是安全工作點(diǎn)，記為Ws；若不滿足，則存在隱患，工作點(diǎn)不安全。

天然氣管網(wǎng)的臨界安全性［14］定義為:對(duì)于一個(gè)安全工作點(diǎn)，是否至少存在一個(gè)負(fù)荷節(jié)點(diǎn)，在其流量增加或減少后，形成的新工作點(diǎn)將不安全。若存在，則運(yùn)行臨界安全，原工作點(diǎn)是臨界安全工作點(diǎn)，簡(jiǎn)稱臨界點(diǎn)，記為Wb。

臨界安全性定義中，流量增加導(dǎo)致的臨界安全性為正臨界性，流量減少導(dǎo)致的臨界安全性為負(fù)臨界性。

1.3 安全域

天然氣管網(wǎng)系統(tǒng)安全域［14］（security region for natural gas pipeline network system，NGS-SR）定義為:管網(wǎng)運(yùn)行時(shí)，所有安全工作點(diǎn)構(gòu)成的集合，記為ΩSR。NGS-SR 在狀態(tài)空間中為封閉區(qū)域，域內(nèi)點(diǎn)安全，域外點(diǎn)存在隱患。

安全邊界［14］定義為:NGS-SR 中所有臨界點(diǎn)構(gòu)成的集合，記為?SR。安全邊界分為上邊界和下邊界:具有正臨界性的全部臨界點(diǎn)構(gòu)成上邊界；具有負(fù)臨界性的全部臨界點(diǎn)構(gòu)成下邊界。

輸氣能力（gas transmission capability，GTC）［14］定義為:天然氣管網(wǎng)安全運(yùn)行時(shí)的最大輸氣量。GTC 點(diǎn)是管網(wǎng)輸氣量最大時(shí)的運(yùn)行狀態(tài)，代表了NGS-SR 中最高效的工作點(diǎn)。

2 基于氣路的天然氣網(wǎng)安全域建模

2.1 穩(wěn)態(tài)氣路模型

異質(zhì)能源統(tǒng)一建模時(shí)，基于氣路建模天然氣管網(wǎng)［17］，其氣壓與流量對(duì)應(yīng)電路的電壓與電流。本節(jié)在動(dòng)態(tài)氣路模型［17］基礎(chǔ)上，建立穩(wěn)態(tài)氣路模型，包括管道穩(wěn)態(tài)氣路與網(wǎng)絡(luò)穩(wěn)態(tài)方程兩部分。

2.1.1 管道的穩(wěn)態(tài)氣路

線性化天然氣動(dòng)量守恒方程，再推導(dǎo)得到管道穩(wěn)態(tài)氣路，如圖1 所示，具體過(guò)程見(jiàn)附錄A。

圖1 天然氣管道的穩(wěn)態(tài)氣路Fig.1 Steady gas circuit of natural gas pipeline

圖1（a）為分布參數(shù)氣路，用于刻畫dx長(zhǎng)度管道微元的氣壓降和流量差。圖1（b）為圖1（a）等效的π形集中參數(shù)氣路［17］。

圖1（a）中，Rg、Lg、kg、Cg、p分別為分布參數(shù)氣阻、氣感、受控氣壓源與氣容、氣壓源氣壓。穩(wěn)態(tài)時(shí)，“氣感短路，氣容斷路”［17］，各元件公式如下:

式中:λ、vb、s、d、θ分別為管道的摩擦系數(shù)、天然氣流速基值、橫截面積、內(nèi)徑和傾角；g 為重力加速度；R和T分別為天然氣氣體常數(shù)和溫度。需注意:vb為氣路建模的主要誤差來(lái)源，其設(shè)定可根據(jù)經(jīng)驗(yàn)或參考文獻(xiàn)［22］。

穩(wěn)態(tài)時(shí)，圖1（b）中的集中參數(shù)計(jì)算公式如下:

式中:Z、kb、Y1、Y2分別為π 形氣路集中參數(shù)的支路阻抗、受控氣壓源和兩端對(duì)地導(dǎo)納；l為管長(zhǎng)。

2.1.2 網(wǎng)絡(luò)的穩(wěn)態(tài)方程

在得到管道π 形集中參數(shù)氣路后，進(jìn)一步考慮管網(wǎng)支路特性［17］與網(wǎng)絡(luò)拓?fù)?，采用文獻(xiàn)［17］方法可推導(dǎo)得到天然氣網(wǎng)絡(luò)穩(wěn)態(tài)方程。文獻(xiàn)［17］推導(dǎo)時(shí)將壓縮機(jī)建模為氣壓源，對(duì)實(shí)際運(yùn)行描述偏理想。在天然氣領(lǐng)域，用壓縮機(jī)壓比模擬含壓縮機(jī)管道的壓力變化更為常見(jiàn)［23］，因此本文借鑒文獻(xiàn)［21］，將壓縮機(jī)建模為氣路變壓器:

式中:pc′、pc、K分別為壓縮機(jī)首、末端氣壓和壓比。

此外，本文規(guī)定節(jié)點(diǎn)流出為正方向。原因如下:管網(wǎng)分析時(shí)，若更關(guān)注負(fù)荷，通常規(guī)定流出為正［13-14］，而本文算例以負(fù)荷為主。

天然氣網(wǎng)絡(luò)穩(wěn)態(tài)方程的具體推導(dǎo)過(guò)程見(jiàn)附錄A，推導(dǎo)結(jié)果如下:

式中:Gn和pn分別為節(jié)點(diǎn)流量、氣壓構(gòu)成的列向量；Ag、Ag+、Ag-分別為節(jié)點(diǎn)-支路關(guān)聯(lián)矩陣、節(jié)點(diǎn)-流出支路關(guān)聯(lián)矩陣和節(jié)點(diǎn)-流入支路關(guān)聯(lián)矩陣，其含義見(jiàn)附錄A；K為壓縮機(jī)壓比方陣，對(duì)于第i行、第i列元素(K)i，i，用(K)i，i=K表示支路i接入了壓比為K的壓縮機(jī)，否則取1；yb和kb分別為支路導(dǎo)納和受控氣壓源構(gòu)成的對(duì)角矩陣，矩陣元素采用π 形氣路集中參數(shù)的穩(wěn)態(tài)值。

引入廣義節(jié)點(diǎn)導(dǎo)納矩陣［17］Y'g:

則天然氣網(wǎng)絡(luò)穩(wěn)態(tài)方程將表示成與電力網(wǎng)絡(luò)穩(wěn)態(tài)方程相統(tǒng)一的數(shù)學(xué)形式:

需注意，上述氣路模型中，天然氣流量均采用質(zhì)量流量；若采用體積流量，需在得到天然氣質(zhì)量流量結(jié)果后，除以標(biāo)準(zhǔn)狀態(tài)下的天然氣密度。

2.2 安全域模型

2.2.1 安全域

得到網(wǎng)絡(luò)穩(wěn)態(tài)方程后，可建模得到與電力網(wǎng)安全域數(shù)學(xué)形式相統(tǒng)一的天然氣網(wǎng)安全域ΩSR。

式中:Gx、Gminn、pax、pin分別為節(jié)點(diǎn)流量和氣壓上、下限構(gòu)成的列向量；A為Ag的左逆矩陣；Cb為管道容量構(gòu)成的列向量；Kmax和Kmin分別為將矩陣K中的K替換為其上、下限值構(gòu)成的矩陣；Ws為安全工作點(diǎn)，需滿足平衡約束h(Ws)和安全約束g(Ws)，依次為:網(wǎng)絡(luò)方程、節(jié)點(diǎn)流量和壓力的上下限、管道流量上限、壓縮機(jī)壓比上下限。

上述安全域模型適用于輸配氣管網(wǎng)，原因如下:對(duì)于不同級(jí)別的管網(wǎng)，其模型都可采用本文的天然氣網(wǎng)絡(luò)方程，只是摩擦系數(shù)的計(jì)算公式不同，并不影響建模過(guò)程與結(jié)果。

2.2.2 安全邊界

模型如式（15）所示，含義如下:Wb∈ΩSR表示工作點(diǎn)Wb位于域內(nèi)，[Gl，1，…，Gl，i，…，Gl，j，…，Gl，m]是Wb中由負(fù)荷構(gòu)成的向量。負(fù)荷i增加后形成新工作點(diǎn)W*，其中G*l，i為其元素。若?ε*≠0，?i=1，2，…，m，使得W*?ΩSR，則Wb是一個(gè)邊界點(diǎn)，全部Wb構(gòu)成安全邊界。若ε*＞0，則Wb具有正臨界性，位于上邊界；若ε*＜0，則Wb位于下邊界。

2.2.3 輸氣能力

輸氣能力的模型如式（16）所示:

式中:GGTC為GTC 點(diǎn)流量；W∈ΩSR表示W(wǎng)是位于域中的安全工作點(diǎn)。

2.3 模型的統(tǒng)一性

附錄B 表B1 對(duì)比了電力網(wǎng)安全域模型［10］、天然氣網(wǎng)安全域傳統(tǒng)模型［14］和本文模型?？梢钥闯觯疚哪Ｐ驮诋愘|(zhì)能源安全域建模中，體現(xiàn)出很好的統(tǒng)一性:傳統(tǒng)模型的平衡約束采用管道壓降方程，數(shù)學(xué)形式與電力網(wǎng)不同；而本文模型平衡約束采用網(wǎng)絡(luò)方程，與電力網(wǎng)模型的數(shù)學(xué)形式一致。本文模型平衡約束能采用網(wǎng)絡(luò)方程的原因是:在氣路中，天然氣流量與氣壓分別對(duì)應(yīng)功率（電流）與電壓，天然氣網(wǎng)可表示為氣阻等元件構(gòu)成的氣路圖。網(wǎng)絡(luò)方程只與氣路圖元件參數(shù)與拓?fù)湎嚓P(guān)，不受限于能量形式。

2.4 模型的保守性

本文的天然氣網(wǎng)穩(wěn)態(tài)安全域未考慮動(dòng)態(tài)過(guò)程。天然氣網(wǎng)傳輸較慢，其管存效應(yīng)可緩解天然氣生產(chǎn)與使用的不平衡［24］，使得域外點(diǎn)可能是安全的。雖然穩(wěn)態(tài)域在判斷安全性上具有一定保守性，但其在安全分析與調(diào)度中仍具備可用性。原因如下:

首先，電力系統(tǒng)存在類似先例。經(jīng)典的基于穩(wěn)定域邊界的主導(dǎo)不穩(wěn)定平衡點(diǎn)法［25］也具有保守性:通過(guò)該方法判斷為穩(wěn)定的情況，一定穩(wěn)定；判斷為不穩(wěn)定的，再通過(guò)時(shí)域仿真判斷。這將大大減少時(shí)域仿真次數(shù)，提高安全評(píng)估效率［25］。

其次，管網(wǎng)運(yùn)行的實(shí)際流量通常小于設(shè)計(jì)流量，使得穩(wěn)態(tài)域可滿足多數(shù)情況下的安全分析與調(diào)度。

2.5 模型的擴(kuò)展性

本文的天然氣網(wǎng)穩(wěn)態(tài)安全域，為擴(kuò)展到動(dòng)態(tài)安全域奠定了基礎(chǔ)。后續(xù)可考慮如下路徑:

路徑1 是直接建立動(dòng)態(tài)安全域模型，再研究求解方法。動(dòng)態(tài)偏微分方程存在求解困難，可以在建模中，借鑒“統(tǒng)一能路理論”［17，22］和“廣義電路分析理論”［16，21］將偏微分方程簡(jiǎn)化為代數(shù)方程。

路徑2 是對(duì)穩(wěn)態(tài)域結(jié)果進(jìn)行修正。動(dòng)態(tài)情況下，安全域是依賴于時(shí)間和分析節(jié)點(diǎn)的函數(shù)，但對(duì)于某時(shí)間斷面的安全域，仍可采用本文提出的穩(wěn)態(tài)建模方法。但穩(wěn)態(tài)建模因未考慮管存而具有保守性，可采用動(dòng)態(tài)仿真工具校驗(yàn)域外點(diǎn)的安全性，也可設(shè)計(jì)修正因子修正域。

在擴(kuò)展到動(dòng)態(tài)安全域時(shí)，還需考慮IES 的多時(shí)間尺度特性［26］。首先，電力系統(tǒng)也具有多時(shí)間尺度特性，可借鑒現(xiàn)有電力系統(tǒng)動(dòng)態(tài)安全域的處理方法［10］。其次，還可借鑒廣義電路理論和統(tǒng)一能路理論，將偏微分方程轉(zhuǎn)化為代數(shù)方程，使得在分鐘級(jí)與小時(shí)級(jí)的時(shí)間尺度內(nèi)，能用統(tǒng)一模型來(lái)表征不同時(shí)間尺度能源網(wǎng)絡(luò)的動(dòng)態(tài)特性［16，22］。

本文模型能擴(kuò)展適應(yīng)IES 的多主體性［26］。當(dāng)簡(jiǎn)化分析市場(chǎng)多主體機(jī)制時(shí)，安全域的統(tǒng)一分析可通過(guò)交換異質(zhì)能源系統(tǒng)的邊界信息，實(shí)現(xiàn)更高層面的多能源全局運(yùn)行的安全和高效；當(dāng)充分考慮市場(chǎng)多主體機(jī)制時(shí)，可通過(guò)對(duì)運(yùn)行工作點(diǎn)和異質(zhì)能源耦合關(guān)系進(jìn)行約束，建立兼顧多主體利益的域模型。

3 安全域求解與觀測(cè)

由于本文安全域模型不是顯式的解析式，難以直接求解，故需進(jìn)行仿真求解:先求解域的邊界點(diǎn)，再對(duì)邊界點(diǎn)擬合得到擬合表達(dá)式。求解到域后，還能進(jìn)一步觀測(cè)其二維/三維視圖。上述方法的具體流程見(jiàn)附錄C 圖C1（a），包括3 步。

3.1 基于氣路的安全域建模

1）輸入管網(wǎng)參數(shù):氣體常數(shù)R，溫度T，管道參數(shù)（摩擦系數(shù)λ等），節(jié)點(diǎn)流量與壓力上下限Gmaxn、Gminn、pmaxn、pminn，壓比上下限Kmax、Kmin。

2）將天然氣管網(wǎng)等效為氣路。

3）跟據(jù)式（2）至式（9）計(jì)算氣路參數(shù)。

4）考慮管網(wǎng)支路特性和拓?fù)浼s束，根據(jù)式（11）至式（13）表示廣義節(jié)點(diǎn)導(dǎo)納矩陣和網(wǎng)絡(luò)方程。

5）根據(jù)管網(wǎng)元件范圍表示安全約束。

6）得到形如式（14）所示的安全域模型。

3.2 安全邊界點(diǎn)計(jì)算

以上邊界點(diǎn)為例說(shuō)明，下邊界點(diǎn)求解類似。

循環(huán)生成初始點(diǎn)W=[Gi，Gj]，對(duì)每個(gè)初始點(diǎn)迭代修正求邊界點(diǎn)；修正時(shí)，修正Gi，保持Gj不變。對(duì)于M（M＞2）維安全域情況，同理選定一個(gè)負(fù)荷流量進(jìn)行修正，其余M-1 個(gè)負(fù)荷流量保持不變。

1）參數(shù)初始化和生成首個(gè)初始點(diǎn)。設(shè)定初始點(diǎn)采樣間距α，對(duì)Gj等間距采樣。計(jì)數(shù)k=1，生成首個(gè)初始點(diǎn)，令Gi，1=Gmaxi或(Gmini+Gmaxi)/2，Gj，1=Gmaxj，其中Gmaxi和Gmini分別為節(jié)點(diǎn)i流量的上、下限。

2）由初始點(diǎn)求邊界點(diǎn)。單個(gè)初始點(diǎn)求邊界點(diǎn)子流程見(jiàn)附錄C 圖C1（b），這是一個(gè)基于二分法的修正過(guò)程，步驟如下。

步驟1:輸入第k個(gè)初始點(diǎn)。

步驟2:設(shè)定修正步長(zhǎng)初值β0、收斂精度ε、修正點(diǎn)安全性標(biāo)志位Fs。ε是初始點(diǎn)迭代修正到邊界點(diǎn)時(shí)修正步長(zhǎng)需滿足的精度，是迭代結(jié)束的一個(gè)條件。Fs是記錄修正點(diǎn)是否安全的變量:若修正點(diǎn)安全，則Fs=1，否則Fs=-1。

步驟3:修正步長(zhǎng)βi=β0，迭代次數(shù)ki=0；校驗(yàn)初始點(diǎn)安全性，根據(jù)安全性結(jié)果設(shè)置Fs初值。

步驟4:將初始點(diǎn)設(shè)為修正點(diǎn)迭代的初值。

步驟5:開始迭代，對(duì)修正點(diǎn)進(jìn)行安全性校驗(yàn)。計(jì)算修正點(diǎn)的天然氣流速基值vb。將修正點(diǎn)的負(fù)荷流量、氣源壓力、壓縮機(jī)壓比、天然氣流速基值代入網(wǎng)絡(luò)方程，計(jì)算修正點(diǎn)狀態(tài)量。

若狀態(tài)量滿足安全約束，則為收斂性校驗(yàn)做準(zhǔn)備，即ki=ki+1，F(xiàn)s(ki)=1；若Fs(ki)Fs(ki-1)＜0，則令βi=βi/2，執(zhí)行步驟6。

若狀態(tài)量不滿足安全約束，則該點(diǎn)位于域外，在邊界上方，需對(duì)Gi進(jìn)行負(fù)向修正，使該點(diǎn)向域內(nèi)移動(dòng)，即ki=ki+1，F(xiàn)s(ki)=-1；若Fs(ki)Fs(ki-1)＜0，則令βi=βi/2，Gi，k(ki)=Gi，k(ki-1)-βi，生成下一個(gè)修正點(diǎn)，執(zhí)行步驟7。

步驟6:收斂性校驗(yàn)。若βi≤ε，滿足收斂性，則輸出Gi，k、Gj，k，得到一個(gè)邊界點(diǎn)，迭代結(jié)束；否則該點(diǎn)在域內(nèi)，位于邊界下方，需對(duì)Gi正向修正，使之向邊界移動(dòng)，即ki=ki+1，Gi，k(ki)=Gi，k(ki-1)-βi，生成下一個(gè)修正點(diǎn)，并令βi=βi/2，返回步驟4。

步驟7:若Gi，k(ki)=Gi，k(ki-1)-βi，初始點(diǎn)在Gi方向上無(wú)對(duì)應(yīng)邊界點(diǎn)，迭代結(jié)束；否則，返回步驟4。

3）若子流程2）修正到邊界點(diǎn)，則記錄該邊界點(diǎn)Wb，k=[Gi，k，Gj，k]，生成下一個(gè)初始點(diǎn):k=k+1，Gi，k=Gi，k-1，Gj，k=Gj，k+α，執(zhí)行子流程4）；否則，已得到全部邊界點(diǎn)，執(zhí)行3.3 節(jié)步驟。

4）若Gj，k＞Gmaxj，則初始點(diǎn)在Gi方向上超出狀態(tài)空間邊界，已得到全部邊界點(diǎn)，執(zhí)行3.3 節(jié)步驟；否則返回子流程2）的步驟2，修正該初始點(diǎn)。

3.3 安全域的表達(dá)式擬合與觀測(cè)

1）利用邊界點(diǎn)對(duì)Gi、Gj分段線性擬合［11］，結(jié)合節(jié)點(diǎn)流量約束，得到安全域擬合表達(dá)式。

2）將求解結(jié)果繪于Gi-O-Gj平面，得到安全域視圖。對(duì)于三維安全域，通過(guò)三維視圖進(jìn)行觀測(cè)；對(duì)于M（M＞3）維安全域，利用文獻(xiàn)［14］方法降維觀測(cè)。

4 算例分析

4.1 算例原始數(shù)據(jù)

先采用5 節(jié)點(diǎn)天然氣管網(wǎng)系統(tǒng)［14］作為測(cè)試算例進(jìn)行驗(yàn)證。系統(tǒng)結(jié)構(gòu)如圖2 所示，其中N1～N5表示節(jié)點(diǎn)，C1表示壓縮機(jī)，參數(shù)見(jiàn)附錄D。

圖2 5 節(jié)點(diǎn)天然氣管網(wǎng)系統(tǒng)結(jié)構(gòu)Fig.2 Structure of 5-bus natural gas network system

4.2 安全域建模求解

4.2.1 安全域建模

基于網(wǎng)絡(luò)方程進(jìn)行NGS-SR 建模，包括2 步。

1）將管網(wǎng)系統(tǒng)表示為如圖3 所示的等值氣路。系統(tǒng)將等效為3 段π 形等值氣路串聯(lián)形式。其中，Zbj、kbj、Ybj，1、Ybj，2（j=1，2，3）分別表示管道1、2、3 的支路阻抗、受控氣壓源和對(duì)地導(dǎo)納。穩(wěn)態(tài)分析時(shí)，氣容和氣感分別作斷路和短路處理。

圖3 5 節(jié)點(diǎn)天然氣管網(wǎng)系統(tǒng)的等值氣路Fig.3 Equivalent gas circuit of 5-bus natural gas network system

對(duì)應(yīng)的形如式（13）所示的天然氣網(wǎng)絡(luò)方程如下:

式中:G1、G2、G3、G5分別為節(jié)點(diǎn)N1、N2、N3、N5的流量；p1、p2、p3、p5分別為節(jié)點(diǎn)N1、N2、N3、N5的壓力。Y'g的詳細(xì)推導(dǎo)過(guò)程與表達(dá)式見(jiàn)附錄E。

2）根據(jù)3.1 節(jié)，得到基于網(wǎng)絡(luò)方程的NGS-SR模型如式（18）所示:

4.2.2 模型求解過(guò)程

首先，對(duì)工作點(diǎn)的安全校驗(yàn)過(guò)程進(jìn)行說(shuō)明。將待校驗(yàn)工作點(diǎn)對(duì)應(yīng)的p1、G2、G3、G5、G3、G5、K與vb代入式（18），得到含4 個(gè)變量（G1、p2、p3、p5）的線性方程組，其系數(shù)矩陣與增廣矩陣的秩相等，存在唯一解，可求得工作點(diǎn)全部狀態(tài)量:若滿足安全約束，則工作點(diǎn)安全；否則不安全。

然后，根據(jù)3.2 節(jié)，求解安全邊界點(diǎn):利用上述安全分析結(jié)果，修正初始點(diǎn)到邊界點(diǎn)。圖4 以上邊界點(diǎn)為例對(duì)此過(guò)程進(jìn)行展示。

以圖4 中Wb，1=[80，297，-377]為例，示意了單個(gè)上邊界點(diǎn)的迭代修正過(guò)程:

圖4 上邊界點(diǎn)求解過(guò)程Fig.4 Solution process of upper boundary points

1）初始點(diǎn)的G2取其上限值，即300 m3/s，修正步長(zhǎng)β的初值取5 m3/s，收斂精度ε取1 m3/s；

2）校驗(yàn)修正點(diǎn)的安全性，根據(jù)校驗(yàn)結(jié)果對(duì)G2進(jìn)行正/負(fù)向修正；

3）經(jīng)過(guò)4 次迭代修正，修正步長(zhǎng)依次為5，2.5，1.25，0.625 m3/s（+/-表示修正方向），修正點(diǎn)安全且β=0.625 m3/s＜1 m3/s，求解得到邊界點(diǎn)Wb，1。

需說(shuō)明，邊界點(diǎn)求解平均迭代7 次，收斂很快。這是由于修正第k個(gè)初始點(diǎn)時(shí)，G2初值取第k-1 個(gè)邊界點(diǎn)G2值，以保證其與邊界點(diǎn)Wb，k很接近。附錄F 表F1 展示了上邊界點(diǎn)詳細(xì)迭代次數(shù)。

最后，根據(jù)3.3 節(jié)，得到NGS-SR 求解結(jié)果，包括:1）上、下邊界點(diǎn)；2）安全域擬合表達(dá)式。同時(shí)，還進(jìn)一步給出了可視化圖像。

4.2.3 求解結(jié)果

1）邊界點(diǎn)

算例共求得上邊界點(diǎn)22 個(gè)，下邊界點(diǎn)7 個(gè)，詳見(jiàn)附錄F。上邊界點(diǎn)包括15 個(gè)GTC 點(diǎn)和7 個(gè)輸氣量非GTC 的臨界點(diǎn)（簡(jiǎn)稱非GTC 點(diǎn)）。

2）NGS-SR 表達(dá)式與可視化圖像

利用邊界點(diǎn)對(duì)安全邊界進(jìn)行分段線性擬合，得到安全域的擬合表達(dá)式如表1 所示，包括4 段。

表1 5 節(jié)點(diǎn)天然氣管網(wǎng)系統(tǒng)NGS-SR 的擬合表達(dá)式Table 1 Fitting expression of NGS-SR of 5-bus natural gas network system

本文擬合區(qū)域與精確域結(jié)果并無(wú)嚴(yán)格包含關(guān)系，但擬合結(jié)果偏保守，有助于保證區(qū)域內(nèi)工作點(diǎn)的安全性。此外，利用文獻(xiàn)［11］方法可計(jì)算得到擬合結(jié)果最大誤差為0.41 m3/s，在0.1% 以內(nèi)，精度較高。

圖5 給出了NGS-SR 的觀測(cè)結(jié)果，與表1 所示4 段表達(dá)式對(duì)應(yīng)。

圖5 5 節(jié)點(diǎn)天然氣管網(wǎng)系統(tǒng)安全域Fig.5 NGS-SR of 5-bus natural gas network system

4.3 與傳統(tǒng)方法對(duì)比

傳統(tǒng)方法［14］基于管道壓降方程對(duì)NGS-SR 建模，本文采用網(wǎng)絡(luò)方程對(duì)NGS-SR 建模。本節(jié)對(duì)比2 種方法的求解結(jié)果和時(shí)間，以驗(yàn)證本文方法的正確性與高效性。首先對(duì)比狀態(tài)量求解及安全分析，因?yàn)槠涫前踩蛴?jì)算的基礎(chǔ)；然后對(duì)比安全域求解。本文的實(shí)驗(yàn)環(huán)境如下:2.50 GHz Intel Core i5-2450 處理器，內(nèi)存4 GB，仿真平臺(tái)MATLAB R2014a。

4.3.1 工作點(diǎn)狀態(tài)量求解及安全分析的對(duì)比

針對(duì)某工作點(diǎn)，通過(guò)2 種方法計(jì)算其狀態(tài)量，再進(jìn)行安全性和臨界性校驗(yàn)。以邊界點(diǎn)Wb，1=[80，297，-377]、隨機(jī)生成的安全工作點(diǎn)［58，155，-213］與不安全工作點(diǎn)［160，223，-383］為例對(duì)比，如表2 所示。表中，Gb1、Gb2、Gb3分別表示管道1～3的天然氣流量。

由表2 可以得出以下結(jié)論。

表2 典型工作點(diǎn)的狀態(tài)變量計(jì)算結(jié)果、安全性與臨界性校驗(yàn)結(jié)果對(duì)比Table 2 Result comparison of state variable calculation, security and criticality verification for typical operating points

1）本文計(jì)算的狀態(tài)量有很小誤差:僅氣壓存在誤差，最大為0.03%。誤差主要由于氣路引入了天然氣流速基值:若基值等于實(shí)際流速，則網(wǎng)絡(luò)方程與管道壓降方程等價(jià)，理論上將不存在誤差；通過(guò)合理取值或修正基值可有效減小誤差［22］。

2）由于狀態(tài)量計(jì)算誤差很小，不會(huì)影響到工作點(diǎn)安全性和臨界性結(jié)果，具體驗(yàn)證見(jiàn)附錄G。這為本文方法求解安全域邊界點(diǎn)奠定了基礎(chǔ)。

3）本文方法狀態(tài)量計(jì)算效率提高約10 倍，原因是建模時(shí)線性化了天然氣流動(dòng)過(guò)程的動(dòng)量守恒方程［17］，避免了傳統(tǒng)方法采用非線性管道壓降方程求解時(shí)的迭代過(guò)程。

4.3.2 安全域求解的對(duì)比

首先，對(duì)比建模結(jié)果。傳統(tǒng)方法［14］結(jié)果如式（19）所示。相較本文模型式（18），式（19）含管道壓降方程這一非線性約束（約束第一式），更為復(fù)雜。

其次，對(duì)比模型求解結(jié)果。以N2、N5為觀測(cè)節(jié)點(diǎn)，以10 m3/s 為采樣間隔，采用傳統(tǒng)方法［14］求解邊界點(diǎn)，與本文結(jié)果對(duì)比，如表3 所示。

由表3 看出，2 種方法求解的邊界點(diǎn)相同，故由此得到的擬合表達(dá)式及NGS-SR 圖像也一致，說(shuō)明了本文方法的正確性。求得邊界點(diǎn)相同的原因如下:網(wǎng)絡(luò)方程引入的狀態(tài)量求解誤差很小，未影響到工作點(diǎn)的安全性與臨界性校驗(yàn)結(jié)果。

表3 傳統(tǒng)方法和本文方法求得的邊界點(diǎn)對(duì)比Table 3 Comparison of boundary points obtained by traditional method and proposed method

最后，比較NGS-SR 求解時(shí)間。對(duì)于測(cè)試算例，傳統(tǒng)方法耗時(shí)582.58 s，本文方法僅耗時(shí)22.07 s，求解效率提高26.40 倍。原因是:1）本文基于網(wǎng)絡(luò)方程的模型求解不存在非線性方程求解的迭代過(guò)程；2）引入二分法，提升了邊界點(diǎn)的修正效率。

此外，本文NGS-SR 計(jì)算效率的提升是由模型改進(jìn)帶來(lái)的，而非由算法改進(jìn)帶來(lái)的。換而言之，其他算法可與本文模型結(jié)合，進(jìn)一步提升計(jì)算效率。

4.4 實(shí)際算例驗(yàn)證

4.4.1 算例原始數(shù)據(jù)

以圖6 的比利時(shí)東南地區(qū)管網(wǎng)［14］為例進(jìn)行驗(yàn)證，包含9 個(gè)節(jié)點(diǎn)N1～N9、11 條管道與3 臺(tái)壓縮機(jī)C1～C3，參數(shù)見(jiàn)附錄H。

圖6 比利時(shí)東南地區(qū)天然氣管網(wǎng)系統(tǒng)結(jié)構(gòu)Fig.6 Structure of southeastern Belgium gas network system

4.4.2 安全域結(jié)果

首先，邊界點(diǎn)求解參數(shù)如下:以5 m3/s 為采樣間隔，對(duì)節(jié)點(diǎn)N3、N5、N8的負(fù)荷在狀態(tài)空間中采樣；修正節(jié)點(diǎn)N9負(fù)荷。對(duì)于首個(gè)初始點(diǎn)，N9節(jié)點(diǎn)處負(fù)荷初值取30 m3/s；修正步長(zhǎng)初值取4 m3/s，其收斂精度取1 m3/s。

上述參數(shù)條件下，共求解到上邊界點(diǎn)338 個(gè)，包括155 個(gè)GTC 點(diǎn)和183 個(gè)非GTC 點(diǎn)，詳見(jiàn)附錄I。下邊界僅為一個(gè)點(diǎn)，即［0，0，125，69，-194］，該點(diǎn)是狀態(tài)空間邊界下限的交點(diǎn)。

接著，利用邊界點(diǎn)進(jìn)行擬合，結(jié)果如表4 所示。

表4 比利時(shí)東南地區(qū)天然氣管網(wǎng)系統(tǒng)的NGS-SR 擬合表達(dá)式Table 4 Fitting expression of NGS-SR of southeastern Belgium gas network system

最后，由于本算例安全域維度大于3，可通過(guò)降維觀測(cè)可視化，附錄J 進(jìn)行了展示。

4.4.3 與傳統(tǒng)方法對(duì)比

關(guān)于求解結(jié)果，傳統(tǒng)方法［14］求得的邊界點(diǎn)與本文相同，進(jìn)而擬合表達(dá)式也與表4 相同。本文方法與傳統(tǒng)方法的差異僅在求解過(guò)程:由于引入流速基值，工作點(diǎn)的狀態(tài)量計(jì)算存在微小誤差。但該誤差并未影響到工作點(diǎn)的安全性和臨界性校驗(yàn)結(jié)果。

關(guān)于求解時(shí)間，傳統(tǒng)方法耗時(shí)164.06 h，本文耗時(shí)3.23 h，效率提高約50.79 倍。可見(jiàn)，隨著算例規(guī)模增加，本文方法的計(jì)算效率優(yōu)勢(shì)更加明顯。

5 結(jié)語(yǔ)

本文借鑒異質(zhì)能源統(tǒng)一建模的思想，對(duì)天然氣網(wǎng)安全域的建模與求解進(jìn)行了研究主要貢獻(xiàn)如下:

1）提出了基于異質(zhì)能源統(tǒng)一建模方法的天然氣網(wǎng)安全域穩(wěn)態(tài)模型，在安全分析上實(shí)現(xiàn)了天然氣網(wǎng)和電力網(wǎng)的數(shù)學(xué)形式統(tǒng)一；

2）提出了天然氣網(wǎng)安全域的全維求解方法，能得到域邊界點(diǎn)與擬合表達(dá)式，并實(shí)現(xiàn)可視化觀測(cè)；

3）通過(guò)經(jīng)典5 節(jié)點(diǎn)算例和比利時(shí)管網(wǎng)實(shí)際算例驗(yàn)證了本文模型及方法的正確性；本文方法效率相對(duì)現(xiàn)有方法提升了1 個(gè)以上數(shù)量級(jí)。

本文工作表明了異質(zhì)能源統(tǒng)一建模方法在天然氣網(wǎng)安全分析中的應(yīng)用前景，為電、氣、熱異質(zhì)能源的綜合安全分析奠定了基礎(chǔ)。本文方法可擴(kuò)展到天然氣網(wǎng)動(dòng)態(tài)安全域，擴(kuò)展后將更具備實(shí)際價(jià)值。下一步研究主要包括考慮動(dòng)態(tài)特性和多主體特性的IES 安全域、IES 安全域的解析表達(dá)式等。

附錄見(jiàn)本刊網(wǎng)絡(luò)版（http：//www.aeps-info.com/aeps/ch/index.aspx），掃英文摘要后二維碼可以閱讀網(wǎng)絡(luò)全文。