廖英祺，荊江平，葉婷

(1.國(guó)網(wǎng)江蘇省電力有限公司南京供電分公司，江蘇 南京 210019；2.國(guó)網(wǎng)江蘇省電力有限公司，江蘇 南京 210024)

區(qū)域多能源系統(tǒng)(regional multi-energy system，RMES)是集電能、天然氣等異質(zhì)能源的產(chǎn)生、傳輸、分配與消費(fèi)于一體的新型能源系統(tǒng)，能夠有效提升用戶側(cè)終端綜合能效和風(fēng)光發(fā)電消納水平，是清潔低碳、安全高效的能源體系的重要組成部分[1-3]。為了保障RMES的安全可靠運(yùn)行、提高用戶的用能質(zhì)量，以及適應(yīng)主動(dòng)管理的需求，調(diào)度管理人員需實(shí)時(shí)感知RMES的運(yùn)行狀態(tài)。作為態(tài)勢(shì)感知的基礎(chǔ)，狀態(tài)估計(jì)(state estimation，SE)能夠?qū)α繙y(cè)數(shù)據(jù)進(jìn)行分析和計(jì)算，為RMES提供更加精確的狀態(tài)信息。

當(dāng)下配電網(wǎng)估計(jì)領(lǐng)域的學(xué)術(shù)研究已較為廣泛[4-6]，眾多學(xué)者在此基礎(chǔ)上展開(kāi)了一系列面向電-氣多能源系統(tǒng)(electric-gas multi-energy system，EGMES)的研究：文獻(xiàn)[7]考慮復(fù)雜天然氣管網(wǎng)的穩(wěn)態(tài)建模及壞數(shù)據(jù)辨識(shí)，提出了穩(wěn)態(tài)EGMES-SE方法；針對(duì)EGMES-SE中存在的配氣網(wǎng)初值和求解效率問(wèn)題，文獻(xiàn)[8-9]提出了抗差及雙線性抗差EGMES-SE方法；文獻(xiàn)[10]計(jì)及實(shí)時(shí)量測(cè)誤差的不確定性，提出一種區(qū)間EGMES-SE方法；文獻(xiàn)[11]考慮EGMES的多管理主體特征和通信隱私問(wèn)題，提出基于交替方向乘子法的分布式估計(jì)求解方法。上述文獻(xiàn)中配氣網(wǎng)均根據(jù)Weymouth方程進(jìn)行穩(wěn)態(tài)建模，忽略了管道內(nèi)天然氣的流量和壓強(qiáng)隨時(shí)間和空間發(fā)生變化的特性，會(huì)產(chǎn)生較大的估計(jì)誤差。而配氣管道相比于輸氣管道儲(chǔ)氣能力較弱，在分鐘級(jí)至小時(shí)級(jí)時(shí)間尺度下，配氣管道的動(dòng)態(tài)特征更為明顯，因此亟需構(gòu)建精確的天然氣動(dòng)態(tài)模型[12]。

天然氣動(dòng)態(tài)模型特性遵循流體力學(xué)定律，通常由質(zhì)量守恒與動(dòng)量守恒偏微分方程組描述，難以用解析法直接求解，目前較為主流的處理方法是通過(guò)有限差分，將連續(xù)的偏微分方程轉(zhuǎn)化為若干組離散的代數(shù)方程來(lái)降低求解難度[13-14]。文獻(xiàn)[15]考慮電、氣、熱子系統(tǒng)間的耦合約束，提出了滿足多場(chǎng)景需求的順序估計(jì)算法；文獻(xiàn)[16]將擴(kuò)展卡爾曼濾波(Kalman filter，KF)算法引入EGMES-SE中，并通過(guò)線性插值法生成偽量測(cè)來(lái)解決電、氣多時(shí)間尺度問(wèn)題；文獻(xiàn)[17]考慮多源量測(cè)設(shè)備采樣周期的延遲問(wèn)題，提出了計(jì)及時(shí)序運(yùn)行特性的EGMES-SE算法；文獻(xiàn)[18]提出異步分布式EGMES-SE策略，根據(jù)配氣網(wǎng)靜態(tài)安全域自主調(diào)整估計(jì)周期。以上研究大大促進(jìn)了動(dòng)態(tài)EGMES-SE的發(fā)展，但仍存在一定的不足：①現(xiàn)有的動(dòng)態(tài)估計(jì)方法均不具備抗差特性；②配氣網(wǎng)有限差分模型會(huì)引入大量的時(shí)、空微元，顯著增加了狀態(tài)變量的數(shù)量，在有限的量測(cè)配置下，配氣網(wǎng)量測(cè)冗余度過(guò)低，導(dǎo)致估計(jì)精度下降；③研究對(duì)象多為跨區(qū)級(jí)系統(tǒng)，以大型輸電、氣系統(tǒng)作為骨干網(wǎng)架，缺乏對(duì)接入間歇性可再生能源的RMES場(chǎng)景的考慮。

針對(duì)上述問(wèn)題，本研究提出基于約束卡爾曼的魯棒RMES-SE方法。首先，采用Lax-Wendroff方法對(duì)配氣網(wǎng)動(dòng)態(tài)過(guò)程進(jìn)行數(shù)值求解，推導(dǎo)出配氣網(wǎng)狀態(tài)空間模型；其次，構(gòu)建配電網(wǎng)狀態(tài)空間模型，通過(guò)耦合約束搭建電、氣子系統(tǒng)之間能量轉(zhuǎn)換的橋梁，為RMES-SE提供模型支撐；最后，以KF算法為基礎(chǔ)，添加量測(cè)噪聲自適應(yīng)算法，增強(qiáng)算法的抗差性能，并計(jì)及配氣網(wǎng)時(shí)序狀態(tài)約束及電-氣邊界耦合約束修正估計(jì)結(jié)果，大幅提升配氣網(wǎng)量測(cè)冗余度和估計(jì)精度。仿真算例根據(jù)改進(jìn)的IEEE 33節(jié)點(diǎn)配電網(wǎng)與10節(jié)點(diǎn)配氣網(wǎng)耦合而成，考慮分布式電源接入的場(chǎng)景，由此驗(yàn)證所提方法的有效性。

1 配氣網(wǎng)狀態(tài)空間建模

相較于配電網(wǎng)而言，配氣網(wǎng)具有更長(zhǎng)的暫態(tài)過(guò)程，主要是因?yàn)樘烊粴獾膫鬏斔俣认鄬?duì)較慢，導(dǎo)致配氣網(wǎng)的響應(yīng)速度較慢。天然氣的流量受到管道首末端壓力差的驅(qū)動(dòng)，并且受到管道長(zhǎng)度、內(nèi)徑、粗糙度、傳輸路徑高度和邊界條件等多方面因素影響；因此，正確地描述配氣網(wǎng)管道內(nèi)氣體流量和壓強(qiáng)隨時(shí)間和空間的變化關(guān)系，是配氣網(wǎng)狀態(tài)空間建模的關(guān)鍵[12]。

1.1 配氣網(wǎng)管道差分建模

在管道水平條件下，不考慮重力加速度和傾角的影響，管道內(nèi)部的暫態(tài)變化過(guò)程，可用沿天然氣管道軸線的一維等溫傳輸動(dòng)態(tài)方程來(lái)描述，即根據(jù)質(zhì)量及動(dòng)量守恒定律得到的一組偏微分方程，具體為[16]：

(1)

(2)

式(1)、(2)中：ta、x分別為時(shí)間、管道距離；p、q分別為天然氣的氣體壓強(qiáng)、質(zhì)量流量；va為氣體聲速；S、D和λ分別為管道的橫截面積、內(nèi)徑和摩擦系數(shù)；w為管道平均流速。

應(yīng)用偏微分方程形式的配氣網(wǎng)管道方程進(jìn)行狀態(tài)估計(jì)，會(huì)大大增加計(jì)算復(fù)雜度；因此，本研究對(duì)配氣網(wǎng)管道的暫態(tài)方程采用Lax-Wendroff有限元差分格式建模，差分后的配氣網(wǎng)管道方程能夠表征相鄰時(shí)刻間流量和壓強(qiáng)的變化關(guān)系，具體可將式(1)、(2)轉(zhuǎn)化為[16]：

pi+1，t+1+pi，t+1-pi+1，t-pi，t+

(3)

(4)

式(3)、(4)中：i為管道差分節(jié)點(diǎn)(包含真實(shí)節(jié)點(diǎn)和虛擬節(jié)點(diǎn))編號(hào)；t為時(shí)刻(時(shí)間斷面)，用下標(biāo)t表示時(shí)刻t的相關(guān)變量，下同；Δx、Δt分別為差分空間步長(zhǎng)、時(shí)間步長(zhǎng)。如圖1所示，配氣網(wǎng)管道差分建模引入了虛擬節(jié)點(diǎn)、支路(將第j條真實(shí)支路差分成M段)，同時(shí)將連續(xù)的偏微分方程轉(zhuǎn)化為若干組離散的代數(shù)方程，以降低求解難度，便于精確獲取短時(shí)間尺度下配氣網(wǎng)波動(dòng)的狀態(tài)參數(shù)。

圖1 配氣網(wǎng)管道有限元差分建模Fig.1 Finite difference model of gas network pipeline

進(jìn)一步將式(3)、(4)轉(zhuǎn)化為統(tǒng)一的矩陣形式，即

(5)

(6)

式(5)、(6)中A1—A5均為常系數(shù)矩陣的元素。

1.2 配氣網(wǎng)狀態(tài)方程推導(dǎo)

假設(shè)配氣網(wǎng)中的真實(shí)節(jié)點(diǎn)數(shù)為N，真實(shí)支路數(shù)量為b，虛擬支路數(shù)量為v(即所有真實(shí)支路差分段數(shù)之和)，則時(shí)刻t的配氣網(wǎng)狀態(tài)列向量

xg，t=(pt，qt)T=(pn，t，pB，t，qt)T∈RN+2v.

(7)

其中：

(8)

(9)

式(7)—(9)中：pt=(pn，t，pB，t)，pn，t、pB，t分別為真實(shí)節(jié)點(diǎn)壓強(qiáng)行向量、虛擬節(jié)點(diǎn)壓強(qiáng)行向量；qt為虛擬支路質(zhì)量流量行向量；pBj，t、qBj，t分別為第j條真實(shí)支路的虛擬節(jié)點(diǎn)壓強(qiáng)行向量、虛擬支路質(zhì)量流量行向量；qBj，0，t、qBj，M，t可等價(jià)為第j條真實(shí)支路首、末端質(zhì)量流量。

每段虛擬支路的配氣網(wǎng)暫態(tài)特性皆可由式(5)描述，因此xg，t+1與xg，t間的關(guān)系式可表示為

J1xg，t+1=J2xg，t，

(10)

式中J1，J2∈R2v×(N+2v)為常系數(shù)矩陣，具體的矩陣元素可根據(jù)式(5)、(6)計(jì)算得到。

進(jìn)一步考慮配氣網(wǎng)節(jié)點(diǎn)平衡關(guān)系，即任意時(shí)刻t下流入、流出任意節(jié)點(diǎn)的質(zhì)量流量之和等于節(jié)點(diǎn)注入質(zhì)量流量，根據(jù)拓?fù)浼s束關(guān)系可構(gòu)建方程

(11)

聯(lián)合式(10)、(11)，可得

(12)

式中零矩陣02∈RN×(N+2v)，03∈R2v。

由式(12)可轉(zhuǎn)化得到配氣網(wǎng)的暫態(tài)狀態(tài)方程

xg，t+1=Fgxg，t+Gg，t=

(13)

式中：Fg∈R(N+2v)×(N+2v)、Gg，t∈RN+2v分別為配氣網(wǎng)的狀態(tài)轉(zhuǎn)移矩陣、控制變量矩陣。計(jì)算過(guò)程中：Fg為常系數(shù)矩陣，與原始管道參數(shù)和配氣網(wǎng)動(dòng)態(tài)特性有關(guān)；Gg，t與下一時(shí)刻各真實(shí)節(jié)點(diǎn)注入質(zhì)量流量qn，t+1有關(guān)，由于配氣網(wǎng)短期負(fù)荷波動(dòng)相對(duì)平緩，可通過(guò)Holt兩參數(shù)指數(shù)平滑法預(yù)測(cè)獲得。

1.3 配氣網(wǎng)狀態(tài)空間建模

配氣網(wǎng)量測(cè)向量zg，t+1主要包括真實(shí)節(jié)點(diǎn)壓強(qiáng)pn，t+1，真實(shí)支路首、末端質(zhì)量流量q0，t+1、qM，t+1，以及真實(shí)節(jié)點(diǎn)注入質(zhì)量流量qn，t+1，可構(gòu)建量測(cè)方程

(14)

式中：zg，t+1∈R2N+2b；Hg∈R(2N+2b)×(N+2v)、J4∈R(N+2b)×(N+2v)為常系數(shù)矩陣，其中J4各行向量中與量測(cè)量相對(duì)應(yīng)的狀態(tài)量位置元素為1，其余元素為0。

由式(13)、(14)可構(gòu)建配氣網(wǎng)狀態(tài)空間模型

(15)

式中wg，t、vg，t+1分別為配氣網(wǎng)過(guò)程向量、量測(cè)噪聲向量，一般服從互不相關(guān)的零均值高斯分布。

2 RMES-SE建模

在配氣網(wǎng)狀態(tài)空間建模的基礎(chǔ)上，本節(jié)進(jìn)一步構(gòu)建配電網(wǎng)狀態(tài)空間模型，并以燃?xì)廨啓C(jī)(gas turbine，GT)和電轉(zhuǎn)氣(power to gas，P2G)為接口，建立電-氣雙向耦合環(huán)節(jié)，為RMES-SE提供模型支撐。本研究采用的RMES-SE方法以KF算法為基礎(chǔ)，通過(guò)自適應(yīng)更新量測(cè)噪聲統(tǒng)計(jì)特性，保證估計(jì)的魯棒性，并計(jì)及配氣網(wǎng)時(shí)序狀態(tài)約束及電-氣邊界耦合約束修正估計(jì)結(jié)果，在提高量測(cè)冗余度的同時(shí)，確保電-氣邊界條件的一致性，實(shí)現(xiàn)對(duì)RMES運(yùn)行狀態(tài)的實(shí)時(shí)感知。

2.1 配電網(wǎng)狀態(tài)空間建模

本研究建立的配電網(wǎng)狀態(tài)空間模型為準(zhǔn)穩(wěn)態(tài)模型[19]，主要面向配電網(wǎng)負(fù)荷波動(dòng)較小、穩(wěn)態(tài)運(yùn)行的場(chǎng)景。配電網(wǎng)的狀態(tài)方程通過(guò)時(shí)間序列預(yù)測(cè)方法建立，即采用Holt兩參數(shù)指數(shù)平滑法；而量測(cè)方程則與常規(guī)配電網(wǎng)穩(wěn)態(tài)估計(jì)相同。配電網(wǎng)狀態(tài)量選取節(jié)點(diǎn)α的電壓幅值Uα和電壓相角θα，量測(cè)量包括節(jié)點(diǎn)電壓幅值Uα，節(jié)點(diǎn)注入有功功率Pα、無(wú)功功率Qα，以及節(jié)點(diǎn)α、β之間支路有功功率Pαβ、無(wú)功功率Qαβ。構(gòu)建配電網(wǎng)狀態(tài)空間模型為

(16)

式中：xe，t+1、ze，t+1分別為配電網(wǎng)狀態(tài)向量、量測(cè)向量；Fe，t、Ge，t為配電網(wǎng)的狀態(tài)轉(zhuǎn)移矩陣、控制變量矩陣，可通過(guò)Holt兩參數(shù)指數(shù)平滑法求解獲得；h(·)為配電網(wǎng)量測(cè)方程；we，t、ve，t+1分別為配電網(wǎng)過(guò)程向量、量測(cè)噪聲向量，誤差分布原理與上述配氣網(wǎng)噪聲相同。

2.2 耦合元件建模

在RMES中，天然氣能與電能一般通過(guò)GT和P2G實(shí)現(xiàn)雙向傳遞。GT通過(guò)燃燒天然氣產(chǎn)生電能，具有啟停靈活、能量梯級(jí)利用的特點(diǎn)；而P2G通過(guò)電解水制備氫氣，并以氫氣為原料制造甲烷后，注入到現(xiàn)有配氣網(wǎng)進(jìn)行規(guī)?；鎯?chǔ)和傳輸，實(shí)現(xiàn)節(jié)能減排[20-21]。本研究中GT和P2G輸入/生成氣流量和輸出/消耗電功率的函數(shù)關(guān)系為[20-21]

(17)

式中：PGT、PP2G分別為GT、P2G輸出/消耗的電功率；qGT、qP2G分別為GT、P2G燃燒/生成的氣流量；GT、P2G分別為GT、P2G的能量轉(zhuǎn)換效率，非熱電聯(lián)產(chǎn)形式下GT效率可達(dá)40%左右，P2G完整化學(xué)反應(yīng)綜合能量轉(zhuǎn)換效率約為45%～60%；LHV為天然氣低熱值[21]，標(biāo)準(zhǔn)狀況下其值為35.40～39.12 MJ/m3，采用其平均值37.26 MJ/m3，標(biāo)準(zhǔn)狀況下天然氣密度0.717 4 kg/m3，經(jīng)轉(zhuǎn)換后LHV=51.94 MJ/kg；ηGT、ηP2G分別為等價(jià)后的GT、P2G能量轉(zhuǎn)換系數(shù)。

2.3 魯棒RMES-SE

受現(xiàn)場(chǎng)工況、通信干擾等因素影響，RMES量測(cè)量中不可避免地出現(xiàn)非高斯噪聲和不良數(shù)據(jù)，若量測(cè)噪聲統(tǒng)計(jì)特性與真實(shí)噪聲情況不匹配，將導(dǎo)致常規(guī)KF算法無(wú)法準(zhǔn)確跟蹤系統(tǒng)狀態(tài)。本文所提方法在KF算法的基礎(chǔ)上，運(yùn)用Huber的魯棒M估計(jì)理論[22]，動(dòng)態(tài)修正量測(cè)噪聲協(xié)方差陣，保證量測(cè)噪聲統(tǒng)計(jì)特性符合實(shí)際工況。在分別獲得電、氣子系統(tǒng)本地估計(jì)值后，考慮電-氣邊界耦合約束及配氣網(wǎng)時(shí)序狀態(tài)約束，采用帶線性等式約束的KF[23]進(jìn)行融合估計(jì)，修正得到全局狀態(tài)一致解。本文所提方法主要分為以下4個(gè)步驟。

a)預(yù)測(cè)步。

(18)

式中：下標(biāo)k∈{e，g}，表示電、氣子系統(tǒng)之一，當(dāng)k=“g”時(shí)，F(xiàn)g，t=Fg，t+1=Fg；xk，f，t+1、Pk，f，t+1分別為時(shí)刻t+1配電/氣網(wǎng)狀態(tài)預(yù)測(cè)向量、預(yù)測(cè)協(xié)方差矩陣；xk，p，t、Pk，p，t分別為時(shí)刻t配電/氣網(wǎng)狀態(tài)估計(jì)向量、估計(jì)協(xié)方差矩陣；Qk，t為過(guò)程噪聲協(xié)方差矩陣。

b)量測(cè)噪聲更新[22]。

(19)

(20)

式(19)、(20)中：Rk，t+1、Rk，c，t+1分別為修正前、后的量測(cè)噪聲方差矩陣；φk，t+1為量測(cè)噪聲尺度因子，φk，s，t+1為φk，t+1的第s個(gè)對(duì)角元素；c為殘差閾值，通常取值為1.3～2.0；rk，s，t+1為歸一化量測(cè)預(yù)測(cè)殘差向量rk，t+1的第s個(gè)元素。rk，t+1計(jì)算公式為：

rk，t+1=diag(zk，t+1-zk，f，t+1)[diag(Pk，f，z，t+1)]-1/2，

(21)

(22)

式(21)、(22)中：zk，f，t+1為量測(cè)預(yù)測(cè)向量，其中ze，f，t+1=h(xe，f，t+1)，zg，f，t+1=Hgxg，f，t+1；Pk，f，z，t+1為量測(cè)預(yù)測(cè)協(xié)方差矩陣；He，t+1=?h(xe，f，t+1)/?xe，f，t+1為配電網(wǎng)雅可比矩陣。

c)濾波步。

(23)

式中Kk，t+1為卡爾曼增益矩陣。

d)約束KF狀態(tài)修正。

在配電、氣網(wǎng)獨(dú)立進(jìn)行魯棒估計(jì)的基礎(chǔ)上，將配氣網(wǎng)時(shí)序狀態(tài)約束和電-氣邊界耦合約束添加到RMES-SE中，在協(xié)調(diào)中心側(cè)采用帶線性等式約束的KF對(duì)全局狀態(tài)進(jìn)行融合估計(jì)[23]。根據(jù)估計(jì)投影法，融合估計(jì)可轉(zhuǎn)化為求解下列問(wèn)題：

(24)

(25)

式(24)、(25)中：xc，t+1、xp，t+1和Pp，t+1分別為電-氣融合修正狀態(tài)向量、估計(jì)狀態(tài)向量和估計(jì)協(xié)方差矩陣；hc(·)為配電網(wǎng)耦合量測(cè)方程；η為RMES能量轉(zhuǎn)換系數(shù)矩陣；Hgc為配氣網(wǎng)耦合量測(cè)系數(shù)矩陣；04∈R(2m-1)×(N+2v)，05∈R(N+2v)×(2m-1)，其中m為配電網(wǎng)節(jié)點(diǎn)數(shù)量。式(24)的解為

xc，t+1=xp，t+1-

(26)

(27)

(28)

式(26)—(28)中：Ht+1、D(·)分別為電-氣融合系數(shù)矩陣、等式約束向量；06∈R2v×(2k-1)；Hec，t+1為配電網(wǎng)耦合量測(cè)對(duì)應(yīng)的雅可比矩陣。

3 算例分析

為驗(yàn)證本文所提方法的性能效果，采用改進(jìn)的IEEE 33節(jié)點(diǎn)配電網(wǎng)和10節(jié)點(diǎn)配氣網(wǎng)[24]構(gòu)建RMES算例作為本研究的測(cè)試系統(tǒng)，共包含32條配電線路、11條配氣管道，并通過(guò)1臺(tái)GT和1臺(tái)P2G實(shí)現(xiàn)電氣耦合，系統(tǒng)結(jié)構(gòu)如圖2所示，其中：配電網(wǎng)中2臺(tái)光伏發(fā)電裝置(PV)經(jīng)電壓源換流器并網(wǎng)，分別接于節(jié)點(diǎn)18、31，均采用有功功率-無(wú)功功率控制策略[24]，即控制交流側(cè)注入總有功、無(wú)功功率為定值；2臺(tái)風(fēng)電機(jī)組(WG)采用直接并網(wǎng)方式，分別接于節(jié)點(diǎn)9、21[25]。

圖2 RMES測(cè)試系統(tǒng)Fig.2 Test system of RMES

仿真實(shí)驗(yàn)中：配電網(wǎng)的基準(zhǔn)電壓為12.66 kV，基準(zhǔn)容量為1 MV·A；配氣網(wǎng)的基準(zhǔn)氣壓為1 bar，基準(zhǔn)流量為1 kg/s，配氣網(wǎng)管道摩擦系數(shù)設(shè)為0.004，氣體聲速設(shè)為400 m/s；配氣網(wǎng)差分空間步長(zhǎng)、時(shí)間步長(zhǎng)分別設(shè)置為500 m、15 min。在潮流真值的基礎(chǔ)上添加服從高斯分布的隨機(jī)噪聲獲得量測(cè)值，相較潮流真值，各量測(cè)值最大誤差分別為：1.5%(電壓幅值實(shí)時(shí)量測(cè))、3%(支路功率、負(fù)荷功率實(shí)時(shí)量測(cè))、30%(根據(jù)歷史數(shù)據(jù)、實(shí)時(shí)天氣預(yù)測(cè)的負(fù)荷功率和直接并網(wǎng)分布式電源出力偽量測(cè))、3%(節(jié)點(diǎn)壓強(qiáng)實(shí)時(shí)量測(cè))、4.5%(支路首末端質(zhì)量流量、節(jié)點(diǎn)注入質(zhì)量流量實(shí)時(shí)量測(cè))。量測(cè)配置信息見(jiàn)表1。仿真環(huán)境的建立采用Intel Core i7-10700 CPU和16 GiB RAM的計(jì)算機(jī)，通過(guò)MATLAB 2016b進(jìn)行求解。

表1 量測(cè)配置信息Tab.1 Measurement configuration

3.1 濾波效果測(cè)試

根據(jù)RMES負(fù)荷曲線和分布式電源出力曲線，以15 min為采樣間隔時(shí)長(zhǎng)，進(jìn)行連續(xù)96個(gè)時(shí)刻的測(cè)試，在無(wú)壞數(shù)據(jù)干擾情況下使用KF(無(wú)量測(cè)噪聲更新和狀態(tài)修正步驟)和本文所提的魯棒約束卡爾曼濾波(robust constrained Kalman filter，RCKF)算法進(jìn)行濾波效果測(cè)試。采用狀態(tài)量的最大估計(jì)誤差Emax和平均估計(jì)誤差Eavg衡量估計(jì)精度：

(29)

(30)

式(29)、(30)中：T為時(shí)刻總數(shù)；n為狀態(tài)變量數(shù)；xt為狀態(tài)向量真實(shí)值，xh，t為其第h個(gè)元素；xc，h，t為狀態(tài)向量融合修正值的第h個(gè)元素。表2展示了不同方法的狀態(tài)量估計(jì)精度對(duì)比結(jié)果。

在某一時(shí)刻分別采用KF與RCKF測(cè)試，使用蒙特卡洛模擬法進(jìn)行1 000次抽樣，計(jì)算電-氣狀態(tài)量的平均估計(jì)值，算法精度對(duì)比結(jié)果如圖3所示。結(jié)合表2和圖3可知：相較KF，RCKF對(duì)配氣網(wǎng)估計(jì)精度有明顯的提升作用；而配電網(wǎng)估計(jì)精度雖有所下降，但基本處于同一數(shù)量級(jí)，對(duì)整體濾波效果影響不大。

圖3 RMES各節(jié)點(diǎn)狀態(tài)估計(jì)精度對(duì)比Fig.3 Comparison of state estimation accuracy for each node in RMES

為進(jìn)一步比較KF、RCKF對(duì)RMES局部濾波效果的影響，選取算例中與GT相關(guān)聯(lián)的電氣耦合量測(cè)，即配電網(wǎng)節(jié)點(diǎn)7注入有功功率和配氣網(wǎng)節(jié)點(diǎn)7注入質(zhì)量流量，并采用相對(duì)估計(jì)誤差來(lái)衡量量測(cè)估計(jì)精度，各時(shí)刻測(cè)試結(jié)果如圖4所示。

圖4 RMES耦合量測(cè)估計(jì)精度對(duì)比Fig.4 Comparison of coupled measurement estimation accuracy in RMES

配氣網(wǎng)估計(jì)精度提升明顯的原因主要是：差分模型下配氣網(wǎng)狀態(tài)量大幅增加，量測(cè)數(shù)目有限條件下，量測(cè)冗余度較低；而本文所提的RCKF通過(guò)添加配氣網(wǎng)時(shí)序狀態(tài)約束和電-氣邊界耦合約束變相提升量測(cè)冗余，能夠在KF基礎(chǔ)上修正得到更為精確的濾波結(jié)果。量測(cè)冗余度對(duì)比見(jiàn)表3。

配電網(wǎng)估計(jì)精度下降的原因是：一方面配電網(wǎng)量測(cè)冗余改變不足以明顯提升估計(jì)精度；另一方面RCKF添加了量測(cè)噪聲更新算法，當(dāng)部分正常量測(cè)數(shù)據(jù)時(shí)序波動(dòng)較大時(shí)，為保障狀態(tài)估計(jì)的魯棒性，會(huì)主動(dòng)修正數(shù)據(jù)統(tǒng)計(jì)特性，降低其對(duì)估計(jì)結(jié)果的影響，導(dǎo)致估計(jì)誤差略微增大。

3.2 抗差性能測(cè)試

為測(cè)試RCKF在壞數(shù)據(jù)干擾下的抗差性能，假設(shè)在第32—第40個(gè)采樣時(shí)刻，RMES中配電網(wǎng)節(jié)點(diǎn)5-6支路首端有功功率和配氣網(wǎng)節(jié)點(diǎn)2-4支路首端質(zhì)量流量量測(cè)通信丟失(即量測(cè)值為0)，在此條件下進(jìn)行動(dòng)態(tài)估計(jì)。圖5展示了RMES量測(cè)通信丟失下的估計(jì)結(jié)果。

圖5 RMES量測(cè)丟失下的估計(jì)結(jié)果Fig.5 Estimation results under measurement loss in RMES

常規(guī)KF算法缺乏對(duì)量測(cè)噪聲特性的自適應(yīng)處理能力，當(dāng)量測(cè)量統(tǒng)計(jì)特性偏離先驗(yàn)統(tǒng)計(jì)特性時(shí)，丟失量測(cè)對(duì)應(yīng)的估計(jì)值大幅偏離真實(shí)值；而本文所提方法在此期間仍能較好跟蹤系統(tǒng)狀態(tài)變化。

除了較為嚴(yán)重的量測(cè)數(shù)據(jù)丟失情況外，當(dāng)量測(cè)受到通信干擾或儀表故障影響，其誤差也會(huì)遠(yuǎn)大于正常的量測(cè)噪聲。為此在各時(shí)刻下，設(shè)壞數(shù)據(jù)由疊加10倍標(biāo)準(zhǔn)差的高斯白噪聲隨機(jī)生成，占子系統(tǒng)量測(cè)的2%～10%，針對(duì)不同比例的壞數(shù)據(jù)場(chǎng)景，分別進(jìn)行1 000次蒙特卡洛仿真實(shí)驗(yàn)，并使用狀態(tài)量的Eavg和Emax來(lái)評(píng)估本文所提方法抗差性能的普適性。表4展示了相應(yīng)的仿真結(jié)果。

表4 不同壞數(shù)據(jù)比例下的估計(jì)結(jié)果Tab.4 Estimation results under different bad data ratios

由表4可知，隨著壞數(shù)據(jù)占比的增加，RMES狀態(tài)量的Eavg和Emax也會(huì)相應(yīng)上升，但仍能夠維持較為穩(wěn)定的數(shù)據(jù)水平。相較于正常量測(cè)下的估計(jì)結(jié)果，數(shù)據(jù)差異不大，這主要是因?yàn)楸疚乃岱椒ú捎没贖uber的魯棒M估計(jì)理論，能夠根據(jù)量測(cè)新息計(jì)算時(shí)變多維的量測(cè)噪聲尺度因子。當(dāng)發(fā)現(xiàn)可疑的量測(cè)數(shù)據(jù)時(shí)，本文所提方法會(huì)修正噪聲統(tǒng)計(jì)特性，進(jìn)而調(diào)整濾波增益，降低可疑數(shù)據(jù)對(duì)狀態(tài)估計(jì)的影響，即使在面對(duì)壞數(shù)據(jù)干擾時(shí)也能夠保持魯棒性能。

3.3 計(jì)算效率測(cè)試

本文所提的RMES-SE方法分為在線估計(jì)和離線計(jì)算2個(gè)部分：在線估計(jì)指接收到實(shí)時(shí)量測(cè)信息后進(jìn)行配電網(wǎng)和配氣網(wǎng)魯棒估計(jì)，并通過(guò)約束修正達(dá)成全局一致解和精度提升；離線計(jì)算指RMES-SE可根據(jù)歷史狀態(tài)信息或負(fù)荷預(yù)測(cè)信息對(duì)下一時(shí)刻的狀態(tài)量進(jìn)行預(yù)測(cè)，預(yù)測(cè)步驟無(wú)需占用系統(tǒng)實(shí)時(shí)算力資源。表5列出了本文所提RCKF方法各步驟計(jì)算耗時(shí)。

表5 RMES-SE計(jì)算效率統(tǒng)計(jì)數(shù)據(jù)Tab.5 Statistical data of calculation efficiency for RMES-SE

由表5可知，本文所提方法整體在線計(jì)算時(shí)間僅需約3×10-3s，能夠滿足實(shí)時(shí)跟蹤RMES動(dòng)態(tài)變化的需求。相較常規(guī)KF方法，RCKF方法增加了量測(cè)噪聲更新和約束修正步驟，能增強(qiáng)算法魯棒性，提高電氣量測(cè)冗余度，達(dá)成邊界條件約束一致，且不會(huì)增加大時(shí)間成本，具有良好的工程應(yīng)用前景。

4 結(jié)論

伴隨著以天然氣分布式能源與可再生能源為主的新型供能技術(shù)的飛速發(fā)展，RMES需求側(cè)的能源多樣化促使子系統(tǒng)間耦合性增強(qiáng)，連鎖故障引發(fā)的風(fēng)險(xiǎn)系數(shù)也在不斷提高；因此，需要實(shí)時(shí)感知RMES的運(yùn)行狀態(tài)，預(yù)知未來(lái)變化趨勢(shì)，保障RMES安全穩(wěn)定運(yùn)行。鑒于此，本文提出了基于約束KF的魯棒RMES-SE方法，得到以下結(jié)論：

a)本文所提方法采用KF算法進(jìn)行各子系統(tǒng)獨(dú)立估計(jì)，并基于Huber的魯棒M估計(jì)理論實(shí)時(shí)修正量測(cè)噪聲信息，原理簡(jiǎn)單，易于實(shí)現(xiàn)，有效保障了估計(jì)的魯棒性。

b)本文所提方法利用耦合元件邊界約束和配氣網(wǎng)時(shí)序狀態(tài)約束，協(xié)同修正各子系統(tǒng)估計(jì)值，獲得全局一致解，可提升量測(cè)冗余度，具有較好的估計(jì)精度與計(jì)算效率。

本文采用管道平均流速對(duì)配氣網(wǎng)偏微分方程進(jìn)行線性化處理，適用于各管道流速在設(shè)定流速附近波動(dòng)的場(chǎng)景；但面對(duì)負(fù)荷突變等場(chǎng)景，可能會(huì)產(chǎn)生較大估計(jì)誤差。后續(xù)將在估計(jì)過(guò)程中添加兼顧計(jì)算效率與結(jié)果精度的流速自適應(yīng)校正步驟。此外，電、氣子系統(tǒng)存在的量測(cè)傳輸非同步性及耦合元件非線性建模等因素對(duì)RMES-SE的影響，也值得加入到今后的探索研究中。