廖英祺，荊江平，葉婷

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

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

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

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

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

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

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

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

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

(1)

(2)

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

應(yīng)用偏微分方程形式的配氣網(wǎng)管道方程進(jìn)行狀態(tài)估計，會大大增加計算復(fù)雜度；因此，本研究對配氣網(wǎng)管道的暫態(tài)方程采用Lax-Wendroff有限元差分格式建模，差分后的配氣網(wǎ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é)點(包含真實節(jié)點和虛擬節(jié)點)編號；t為時刻(時間斷面)，用下標(biāo)t表示時刻t的相關(guān)變量，下同；Δx、Δt分別為差分空間步長、時間步長。如圖1所示，配氣網(wǎng)管道差分建模引入了虛擬節(jié)點、支路(將第j條真實支路差分成M段)，同時將連續(xù)的偏微分方程轉(zhuǎn)化為若干組離散的代數(shù)方程，以降低求解難度，便于精確獲取短時間尺度下配氣網(wǎ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)中的真實節(jié)點數(shù)為N，真實支路數(shù)量為b，虛擬支路數(shù)量為v(即所有真實支路差分段數(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分別為真實節(jié)點壓強行向量、虛擬節(jié)點壓強行向量；qt為虛擬支路質(zhì)量流量行向量；pBj，t、qBj，t分別為第j條真實支路的虛擬節(jié)點壓強行向量、虛擬支路質(zhì)量流量行向量；qBj，0，t、qBj，M，t可等價為第j條真實支路首、末端質(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ìn)一步考慮配氣網(wǎng)節(jié)點平衡關(guān)系，即任意時刻t下流入、流出任意節(jié)點的質(zhì)量流量之和等于節(jié)點注入質(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)移矩陣、控制變量矩陣。計算過程中：Fg為常系數(shù)矩陣，與原始管道參數(shù)和配氣網(wǎng)動態(tài)特性有關(guān)；Gg，t與下一時刻各真實節(jié)點注入質(zhì)量流量qn，t+1有關(guān)，由于配氣網(wǎng)短期負(fù)荷波動相對平緩，可通過Holt兩參數(shù)指數(shù)平滑法預(yù)測獲得。

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

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

(14)

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

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

(15)

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

2 RMES-SE建模

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

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

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

(16)

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

2.2 耦合元件建模

在RMES中，天然氣能與電能一般通過GT和P2G實現(xiàn)雙向傳遞。GT通過燃燒天然氣產(chǎn)生電能，具有啟停靈活、能量梯級利用的特點；而P2G通過電解水制備氫氣，并以氫氣為原料制造甲烷后，注入到現(xiàn)有配氣網(wǎng)進(jìn)行規(guī)?；鎯蛡鬏?，實現(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分別為等價后的GT、P2G能量轉(zhuǎn)換系數(shù)。

2.3 魯棒RMES-SE

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

a)預(yù)測步。

(18)

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

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

(19)

(20)

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

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為量測預(yù)測向量，其中ze，f，t+1=h(xe，f，t+1)，zg，f，t+1=Hgxg，f，t+1；Pk，f，z，t+1為量測預(yù)測協(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)獨立進(jìn)行魯棒估計的基礎(chǔ)上，將配氣網(wǎng)時序狀態(tài)約束和電-氣邊界耦合約束添加到RMES-SE中，在協(xié)調(diào)中心側(cè)采用帶線性等式約束的KF對全局狀態(tài)進(jìn)行融合估計[23]。根據(jù)估計投影法，融合估計可轉(zhuǎn)化為求解下列問題：

(24)

(25)

式(24)、(25)中：xc，t+1、xp，t+1和Pp，t+1分別為電-氣融合修正狀態(tài)向量、估計狀態(tài)向量和估計協(xié)方差矩陣；hc(·)為配電網(wǎng)耦合量測方程；η為RMES能量轉(zhuǎn)換系數(shù)矩陣；Hgc為配氣網(wǎng)耦合量測系數(shù)矩陣；04∈R(2m-1)×(N+2v)，05∈R(N+2v)×(2m-1)，其中m為配電網(wǎng)節(jié)點數(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)耦合量測對應(yīng)的雅可比矩陣。

3 算例分析

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

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

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

表1 量測配置信息Tab.1 Measurement configuration

3.1 濾波效果測試

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

(29)

(30)

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

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

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

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

圖4 RMES耦合量測估計精度對比Fig.4 Comparison of coupled measurement estimation accuracy in RMES

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

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

3.2 抗差性能測試

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

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

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

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

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

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

3.3 計算效率測試

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

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

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

4 結(jié)論

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

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

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

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