基于Gauss-Markov模型的電力系統(tǒng)t型抗差狀態(tài)估計

顏全椿，衛(wèi)志農(nóng)，徐泰山，王勝明，孫國強

（1.河海大學 可再生能源發(fā)電技術教育部工程研究中心，江蘇 南京 210098；2.國網(wǎng)電力科學研究院/南京南瑞集團公司，江蘇 南京 210003）

0 引言

自1794年Gauss創(chuàng)立最小二乘LS（Least Square）法以來，特別是Gauss-Markov定理建立后，最小二乘估計一直被作為一個良好的估計方法而被廣泛采用。傳統(tǒng)的基于最小二乘及由其派生出的快速解耦法狀態(tài)估計在實際電網(wǎng)已有多年的運行經(jīng)驗，該方法是對服從高斯分布樣本的最優(yōu)估計，具有模型簡單、計算量小、收斂速度快等優(yōu)點[1-2]。然而，實際量測系統(tǒng)中往往存在誤差較大的不良數(shù)據(jù)，而這些不良數(shù)據(jù)并不服從高斯分布，若不能將這些不良數(shù)據(jù)剔除，將嚴重影響狀態(tài)估計結果的精度。因此，研究如何盡可能地降低不良數(shù)據(jù)對電力系統(tǒng)狀態(tài)估計的影響具有重要意義[3]。

針對該問題，國內(nèi)外學者的研究主要側重2個方面：一方面是在最小二乘估計的基礎上，嵌入一個含不良數(shù)據(jù)辨識的模塊[4-6]；另一方面則是以M估計為代表的抗差狀態(tài)估計[7-12]。對于前一方面研究，目前應用較廣泛的是估計辨識法，即首先找出殘差較大的量測作為不良數(shù)據(jù)，用估計值直接修正原來的狀態(tài)估計，修正后的狀態(tài)估計恰恰是不良數(shù)據(jù)的最優(yōu)估計。該方法無需重新進行狀態(tài)估計，實時性好，但存在僅可辨識多個弱相關的不良數(shù)據(jù)、較難處理含多個強相關的不良數(shù)據(jù)的不足。對于后一方面研究，M估計的抗差性和效率取決于狀態(tài)量初值的可靠性、等價權函數(shù)及其臨界值的合理性。電力系統(tǒng)狀態(tài)估計中使用比較廣泛的權函數(shù)有Huber權函數(shù)、Hampel權函數(shù)、Turkey權函數(shù)和含有相關量測的IGGШ權函數(shù)等。不同等價權函數(shù)的選取相當于建立了不同的誤差分布模型。統(tǒng)計界近年來提出并有較深入理論研究支撐的t型估計，以t分布對量測誤差建模，屬于帶有刻度參數(shù)的M估計[13-15]。當t分布取較大的自由度時，趨于高斯分布，對應最小二乘估計，在量測誤差為高斯分布時具有較高的效率；當取較小的自由度時，為柯西分布，相應的估計具有很好的抗差性。恰當?shù)剡x取自由度，可以使t型估計兼顧狀態(tài)估計的抗差性和效率。

本文在前人工作基礎上，將t型估計引入電力系統(tǒng)狀態(tài)估計中，針對t型估計求解較繁瑣的特點，提出近似模型，并以零注入節(jié)點方程為等式約束。隨后，利用牛頓法進行求解，結合一個簡單的電力系統(tǒng)算例分析該方法在實際應用中可能存在的局部最優(yōu)點問題，通過分步調整自由度保證取得全局最優(yōu)解。最后，借助IEEE標準系統(tǒng)和某實際電網(wǎng)測試該方法的有效性，結果表明，與含不良數(shù)據(jù)辨識功能加權最小二乘（WLS）估計（簡稱 WLS+BD）和二次-常數(shù)QC（Quadratic-Constant）估計進行比較，本文方法處理不良數(shù)據(jù)能力更強。

1 t型估計模型及其求解

1.1 電力系統(tǒng)t型估計模型

在給定網(wǎng)絡結線、支路參數(shù)和量測系統(tǒng)的條件下，電力系統(tǒng)的量測方程為[1]：

其中，x為n×1維狀態(tài)變量，包括除平衡節(jié)點外的所有節(jié)點電壓幅值和相角；zi、hi、ri分別為 m× 1 維量測向量、量測函數(shù)向量和殘差向量的第i個分量；m、n分別為量測量及狀態(tài)量的個數(shù)。

考慮量測誤差的Gauss-Markov模型為：

其中，σ0為未知的單位權方差因子；E（r）、D（r）分別為量測殘差的均值和方差；R為量測權重對角矩陣。

對于式（2），定義t型抗差狀態(tài)估計為如下的最優(yōu)化問題[13]：

由于在實際運行中，系統(tǒng)必須滿足零注入節(jié)點的等式約束，考慮這些約束則式（3）可寫為：

其中，g（x）為潮流約束及零注入節(jié)點功率量測的等式約束。

1.2 求解

利用內(nèi)點罰函數(shù)法對式（4）進行處理，則式（4）可轉化為以下無約束的優(yōu)化問題[9]：

其中，L（x，λ）為拉格朗日函數(shù)；λ為p維拉格朗日乘子向量；gj（x）表示第 j個等式約束。

式（5）中 L（x，λ）的一階最優(yōu)條件為：

其中，H為h（x）對x的雅可比矩陣；G為零注入節(jié)點等式約束g（x）對x的雅可比矩陣。

通過牛頓法求解式（6）的非線性方程，對式（6）求偏導可得增廣拉格朗日函數(shù)L（x，λ）的海森矩陣：

式中，?2L/?x2部分忽略了二階項，在形式上與傳統(tǒng)WLS類似，具有較好的兼容性。

令 x（k+1）=x（k）+Δ x（k+1），λ（k+1）=λ（k）+Δλ（k+1），即可對式（4）的t型抗差估計模型進行求解。

因此，可得如下修正公式：

2 局部最優(yōu)解情況分析

從數(shù)學上看，t型估計是一個典型的非線性優(yōu)化問題，可能存在局部最優(yōu)解。以2節(jié)點系統(tǒng)為例進行測試，該系統(tǒng)的潮流分布及網(wǎng)絡參數(shù)參見文獻[9]。由于2節(jié)點系統(tǒng)平衡節(jié)點B1的電壓為1.05∠0°p.u.，不參與迭代，因此優(yōu)化變量僅有節(jié)點B2的電壓幅值與相角。

分別取V=0.01和V=1，t型估計目標函數(shù)隨著節(jié)點B2電壓幅值和相角變化的曲面如圖1、2所示，圖中電壓幅值為標幺值，后同。為了便于觀察，將目標函數(shù)J（x）符號取反，圖3、4給出了相應的等值線。

圖1 V=0.01時t型估計目標函數(shù)J（x）的變化曲面Fig.1 Curved surface of objective function J（x）of t-type estimation when V is 0.01

圖2 V=1時t型估計目標函數(shù)J（x）的變化曲面Fig.2 Curved surface of objective function J（x） of t-type estimation when V is 1

圖3 V=0.01時2節(jié)點系統(tǒng)目標函數(shù)等值線圖Fig.3 Contour map of objective function of 2-bus system when V is 0.01

圖4 V=1時2節(jié)點系統(tǒng)目標函數(shù)等值線圖Fig.4 Contour map of objective function of 2-bus system when V is 1

由圖1及圖3可知，V=0.01時，全局最優(yōu)點為1.005∠-29.22°p.u.，與真實值 1.002∠-30°p.u.誤差較小，但此時目標函數(shù)出現(xiàn)多個局部最優(yōu)點。圖2中，V=1，此時局部最優(yōu)點消失了，但是全局最優(yōu)點為1.020∠-16.62°p.u.，偏離了真實值。進一步，本文對取不同自由度時全局最優(yōu)點的影響進行分析，結果如表1所示。由表可知，自由度V取值越小，估計結果的精度越高，全局最優(yōu)點與真實值越接近。

表1 自由度V對全局最優(yōu)解的影響Tab.1 Impact of V on globally optimal solution

實際上，V值的改變與文獻[9]中的Parzen窗寬度類似，在數(shù)學上其作用是等價的。由圖4可知，局部最優(yōu)點與全局最優(yōu)點距離較遠，當節(jié)點電壓在真實值的較大鄰近區(qū)域內(nèi)取初值時，不會陷入局部最優(yōu)點。因此，如何確定初值是本文算法的關鍵。對于工程應用，可將前一斷面的狀態(tài)估計結果作為后一斷面狀態(tài)估計的初值，根據(jù)電力系統(tǒng)運行狀態(tài)一般變化較平緩的特點，能有效減少迭代次數(shù)，從而提高程序收斂速度。同時，文獻[9]采用逐步調整Parzen窗寬的策略進行多次迭代求解，能有效克服局部最優(yōu)解問題，但存在收斂速度較慢的不足。本文對自由度進行動態(tài)設置[7]：

由式（9）可知，自由度的確定與量測系統(tǒng)的殘差大小有關。t型估計的計算流程如圖5所示。

3 方法特點

與以往的狀態(tài)估計方法相比，本文方法特點如下。

圖5 自適應t型狀態(tài)估計流程圖Fig.5 Flowchart of adaptive t-type estimation

a.傳統(tǒng)WLS狀態(tài)估計方法結果易受不良數(shù)據(jù)影響，本文方法以量測殘差對數(shù)函數(shù)值作為目標函數(shù)，能夠有效減小不良數(shù)據(jù)對狀態(tài)估計結果的影響，抗差性強，估計結果合理性好。

c.該方法在循環(huán)迭代求解過程中通過對自由度動態(tài)調整，可有效克服局部最優(yōu)解問題[9]。

4 算例分析

為了驗證本文方法的有效性，在Visual C++平臺上編制了自適應t型抗差估計。對IEEE標準系統(tǒng)和某實際輸電網(wǎng)的數(shù)據(jù)進行試驗，通過將潮流結果加2%高斯隨機誤差獲得量測生數(shù)據(jù)，通過對生數(shù)據(jù)置0的方式得到試驗用的不良數(shù)據(jù)。

4.1 抗差性測試

4.1.1 4節(jié)點系統(tǒng)

由殘差方程知，某量測上出現(xiàn)不良數(shù)據(jù)時，對自身及其他量測的殘差均產(chǎn)生影響。由于量測誤差的不確定性，反映在殘差上的情形也非常復雜，本文僅考慮2個量測間的相關性，r與量測誤差v存在以下關系：

其中，K為殘差靈敏度矩陣。

對于t型估計，若量測i的殘差較大，則對角元素 Fii（x）=2wi/（V+r2i）≈0，即對該量測降權處理，由此可以得出，t型估計模型中抑制不良數(shù)據(jù)的影響是可能的。

以下針對4節(jié)點系統(tǒng)[16]進行測試，分析含強相關不良數(shù)據(jù)時狀態(tài)估計結果的準確性。該系統(tǒng)的結構、量測分布和節(jié)點電壓真值如圖6所示。

圖6 4節(jié)點系統(tǒng)結構圖Fig.6 Structure of 4-bus system

在給定的量測系統(tǒng)和網(wǎng)絡結構下，找出具有強相關性的量測對（CM），并分析不良數(shù)據(jù)對該量測估計值的影響。相關性分析可通過統(tǒng)計K中的相關系數(shù)，量測i與量測j之間的相關度指標[12]為：

其中，Kij為K的第i行第j列元素值。

當相關度指標滿足0.752 1＜Pij＜0.9498時，為強相關量測對。該系統(tǒng)含4組強相關量測，各估計方法的強相關量測指標見表2。這4組量測對于WLS估計均為強相關量測，而采用t型估計后，相關度指標顯著減小。

表2 4節(jié)點系統(tǒng)強相關量測指標Tab.2 Strong correlation measurement index of 4-bus system

為進一步測試不同方法對強相關量測不良數(shù)據(jù)的估計精度，對上述4組強相關量測添加不良數(shù)據(jù)，不同方法的估計結果見表3，表中數(shù)據(jù)均為標幺值。

表3 4節(jié)點系統(tǒng)強相關量測不同方法的估計結果比較Tab.3 Comparison of estimation results among different strong correlation measurement methods of 4-bus system

以支路1-2有功功率P1-2為例說明不同方法的抗差性能，由表3可知，該量測的真值為0.3883 p.u.，WLS估計結果為 0.2465 p.u.，殘差為0.1418 p.u.，而 t型估計值為0.3876 p.u.，殘差僅為0.0007 p.u.，很好地辨識出了不良數(shù)據(jù)。此外，對于該支路P2-1功率量測，t型估計的誤差也僅為-0.0006 p.u.，表明t型估計成功辨識出了該組強相關不良數(shù)據(jù)。

4.1.2 IEEE 57節(jié)點系統(tǒng)

采用以下指標衡量估計的狀態(tài)變量（除平衡節(jié)點外其余節(jié)點的電壓幅值和相角）的準確性[10]：

其中，x*i為第i個狀態(tài)變量的估計值；xi為狀態(tài)變量的真實值。

表4給出了57節(jié)點系統(tǒng)的試驗結果，該系統(tǒng)含421個量測。由表4可知，本文方法隨著不良數(shù)據(jù)比例的增加，S1指標和S2指標變化較小，表明本文方法不易受不良數(shù)據(jù)的影響。

表4 57節(jié)點系統(tǒng)試驗結果Tab.4 Test results of 57-bus system

4.1.3 某實際省網(wǎng)

為了驗證本文方法在實際電網(wǎng)中的性能，對某輸電網(wǎng)進行測試，該系統(tǒng)含1 119條母線、1 992條線路（簡稱省網(wǎng)A）。由于實際系統(tǒng)中量測的真實值無法得到，本文對不同狀態(tài)估計結果進行合格率評估[10]。表5、6給出了各類型量測的考核基準值與合格門檻值，不同方法的狀態(tài)估計合格率如表7所示。

由表7可知，6個時間斷面中，WLS+BD估計的合格率最低為93.54%，而本文方法的合格率最低為96.24%，相比WLS+BD提高了2.70%。QC估計最低合格率為95.06%，較WLS+BD有一定程度的提高，但仍低于t估計。試驗結果表明本文方法具有較強的抗差性。

表5 狀態(tài)估計合格評估基準值Tab.5 Benchmark of state estimation conformity assessment

表6 狀態(tài)估計合格評估門檻值Tab.6 Threshold of state estimation conformity assessment

表7 某實際輸電網(wǎng)試驗結果Tab.7 Test results of a practical transmission network

4.2 自由度對估計結果的影響測試

第2節(jié)討論了本文方法的局部最優(yōu)解問題，提出了動態(tài)調整自由度的策略[7，9]，并對 2節(jié)點系統(tǒng)進行了分析，由試驗結果可知自由度的選擇及初值的選取方法將決定算法的最終結果。本文對IEEE 57節(jié)點系統(tǒng)進一步分析，采用定V方式（以WLS狀態(tài)估計結果作為系統(tǒng)狀態(tài)量初值）進行t型狀態(tài)估計計算。同時，為了驗證本文動態(tài)調整自由度V方式的有效性，表8給出了上述不同自由度下狀態(tài)估計的結果比較。

由表8可知，隨著自由度V的減小，S1及S2指標逐漸減小，即較小的自由度具有更強的抗差能力。當自由度V=1時，S1指標為0.27415，平均誤差為0.00243，最大誤差為0.00481。在試驗中發(fā)現(xiàn)，本文方法在第1次迭代后利用式（9）計算得到V=0.0129，繼續(xù)進行迭代計算，經(jīng)過5次迭代收斂后V=0.00268，滿足收斂條件。最終t型估計結果S1指標進一步減小至0.01638，僅為第1次迭代時的5.975%。同時，與定V方式比較可知，該方法并未陷入局部最優(yōu)解。在計算效率上，定V方式經(jīng)4次迭代總時間為26.014 ms，而本文僅為11.892 ms，表明本文混合策略在保證求解精度的前提下，顯著提高了收斂速率。

表8 不同V值的t型估計結果Tab.8 Results of t-type estimation for different V values

為了進一步驗證本文動態(tài)調整自由度策略在實際省網(wǎng)中的有效性，圖7給出省網(wǎng)A的120個量測斷面狀態(tài)估計時間。由圖7可知，采用文獻[9]的定V方式時，由于每次狀態(tài)估計需4次循環(huán)迭代求解，平均計算時間為0.4447 s，而采用本文方法后，平均計算時間為0.367 2 s，相比文獻[9]可減少17.43%。同時，本文方法估計結果具有與定V方式同樣的精度，這也進一步說明了本文動態(tài)調整自由度策略在保持狀態(tài)估計求解精度的前提下，可有效減少t型估計的計算時間。實際上，目前抗差估計的計算量往往較大，這也是制約其在實際中進一步應用的原因[9]。對于該省網(wǎng)，實際現(xiàn)場采用的WLS估計平均計算時間為0.1765 s，僅為本文方法的48.07%，因此，如何進一步提高計算速度是今后一個重要的研究內(nèi)容。

圖7 某實際省網(wǎng)t型抗差狀態(tài)估計計算時間比較Fig.7 Calculation time of t-type state estimation for a practical provincial network

5 結論

本文基于t型估計，建立了適用于電力系統(tǒng)狀態(tài)估計的簡化模型，利用其目標函數(shù)連續(xù)可微的特點給出類牛頓法的求解方法。該方法具有如下特點。

a.以t型分布對誤差建模，當取較小自由度時，具有良好的抗差性能，有效排除不良數(shù)據(jù)對狀態(tài)估計結果的影響。

b.提出自適應動態(tài)調整自由度的策略，兼顧了t型估計的效率和抗差性，滿足工程要求。

c.作為一種新的抗差估計，t型估計有很好的理論基礎，將t型估計引入電力系統(tǒng)狀態(tài)估計中，能夠豐富和發(fā)展抗差狀態(tài)估計的理論和方法。與傳統(tǒng)WLS估計相比，本文方法計算量仍然較大。