Givens正交變換法狀態(tài)估計中基于VPAIR消元的改進(jìn)策略

劉 穎 ，童偉林 ，鄒 宇 ，王建全，夏彥輝 ，肖譚南

（1.南京國電南自電網(wǎng)自動化有限公司，江蘇 南京 211153；2.浙江大學(xué) 電氣工程學(xué)院，浙江 杭州 310027）

0 引言

電力系統(tǒng)狀態(tài)估計是能量管理系統(tǒng)（EMS）的重要組成部分，也是穩(wěn)定控制、經(jīng)濟(jì)調(diào)度等其他電力系統(tǒng)應(yīng)用的重要基礎(chǔ)，在電網(wǎng)安全運(yùn)行中承擔(dān)著重要作用［1-2］。它利用實(shí)時測量系統(tǒng)的冗余度來提高數(shù)據(jù)精度，給出系統(tǒng)運(yùn)行狀態(tài)的估計值。

狀態(tài)估計一般用非線性加權(quán)最小二乘法問題加以描述。正則方程法是早期狀態(tài)估計的主要方法，包括加權(quán)最小二乘法、快速分解法等。但是，正則方程系數(shù)矩陣的條件數(shù)是量測雅可比矩陣條件數(shù)的平方，方程的病態(tài)性導(dǎo)致計算精度下降。在處理虛擬量測的情況下，方程病態(tài)問題更加突出。

對此，一類研究較多的方法是利用拉格朗日乘子法來處理等式約束［3］，但是該方法存在信息矩陣不正定和計算量大的不足。Cholesky分解法［4］能夠提升算法的數(shù)值穩(wěn)定性、克服信息矩陣不正定的情況，但計算效率仍然有待提高。使用修正牛頓法［5-7］處理等式約束計算簡單，可以有效規(guī)避上述問題。

另一類研究方法則是利用正交變換法。正交變換法［8-9］是一種可以有效提高數(shù)值穩(wěn)定性的計算方法，在電力系統(tǒng)狀態(tài)估計中應(yīng)用廣泛。Givens旋轉(zhuǎn)法及Household鏡像法是正交變換的2種主要方法。Household鏡像法適用于滿矩陣的情況，而狀態(tài)估計的加權(quán)雅可比矩陣一般比較稀疏，使用Givens正交變換法具有更高的效率。

在進(jìn)行Givens正交變換過程中會產(chǎn)生大量的中間注入元素，增加了Givens旋轉(zhuǎn)的次數(shù)和計算量，影響計算效率。選擇合適的列編號順序和消元策略，可以減少Givens正交變換的計算量。列編號順序一般采用最小度原則，可以保證分解后矩陣的稀疏性。

而Givens正交變換的消元策略可分為逐行消元和逐列消元。早期的消元策略一般為逐行消元策略和定轉(zhuǎn)軸逐列消元策略。而后續(xù)研究發(fā)現(xiàn)，變轉(zhuǎn)軸逐列消元策略可以在很大程度上減少中間注入元素的數(shù)量，有效地改善算法的性能。文獻(xiàn)［10］提出將最稀疏的一行作為固定轉(zhuǎn)軸，其余行與轉(zhuǎn)軸比較，按照注入元素的個數(shù)依次進(jìn)行Givens旋轉(zhuǎn)。文獻(xiàn)［11］提出基于矩陣稀疏性的變轉(zhuǎn)軸策略，選擇最稀疏的2行進(jìn)行Givens旋轉(zhuǎn)。文獻(xiàn)［12］提出基于最少注入元的VPAIR變轉(zhuǎn)軸策略，選擇產(chǎn)生中間注入元素最少的2行進(jìn)行Givens旋轉(zhuǎn)。文獻(xiàn)［13］對幾種策略進(jìn)行了比較，證明采用最小度原則的列排序和VPAIR逐列變轉(zhuǎn)軸策略是一種效率更高的方法。變轉(zhuǎn)軸消元策略的正交變換法在狀態(tài)估計中有廣泛的應(yīng)用［14-16］。

本文在VPAIR逐列變轉(zhuǎn)軸策略的基礎(chǔ)上，提出在其消元過程中動態(tài)調(diào)整量測排序的改進(jìn)策略。與傳統(tǒng)的策略相比，所提策略對量測排序進(jìn)行了在線動態(tài)優(yōu)化，可以進(jìn)一步有效地減少注入元數(shù)量，提高計算效率。

1 正交變換法狀態(tài)估計

狀態(tài)估計用非線性加權(quán)最小二乘法問題描述如下：

其中，x為n維狀態(tài)向量；z為m維量測向量；h（x）為m維量測方程向量；W為m×m階量測權(quán)重矩陣（對角矩陣）。

對量測方程h（x）進(jìn)行泰勒級數(shù)展開，取線性項(xiàng)，有：

其中，H為量測雅可比矩陣。將式（2）代入式（1）可得到：

其中，r=z-h（x0）為 m 維量測殘差向量；rw=W1/2r為m維加權(quán)量測殘差向量；Hw=W1/2H為m×n階加權(quán)量測雅可比矩陣；‖·‖2為歐幾里得范數(shù)。

假設(shè)存在 m×m 階正交矩陣 Q，使得式（4）、（5）成立，那么可得式（6）。

其中，R為n×n階上三角矩陣；Z為（m-n）×n階全零矩陣；b1為n維向量；b2為m-n維向量。

由于范數(shù)大于或等于 0，當(dāng)式（7）成立時，J（Δx）達(dá)到最小值‖b2‖2。

因此，只需要構(gòu)造矩陣Q就可以變換出R和b1，求解式（7）可以得到狀態(tài)向量修正量Δx。矩陣Q可以通過對Hw的Givens變換或者Household變換獲得。

2 Givens旋轉(zhuǎn)

Givens正交變換過程中每次Givens旋轉(zhuǎn)僅消去1個非零元素，直到對角線左邊沒有非零元素為止。以第i行和第j行2行元素為例說明。

當(dāng)對第j行第k列以第i行第k列為轉(zhuǎn)軸進(jìn)行Givens旋轉(zhuǎn)時，第i、j行第k列左邊的元素均為0，以使第j行第k列的元素變?yōu)?。初始的第i、j行元素如下：

經(jīng)過單次Givens變換，得到如下結(jié)果：

Givens變換推導(dǎo)后得到計算公式如下：

3 VPAIR消元策略

VPAIR消元策略是Robey提出的基于最少注入元的變轉(zhuǎn)軸逐列消元策略，該策略經(jīng)過大量算例驗(yàn)證，計算效率高，是目前被廣泛采用的方法。其原則如下：

a.檢查是否存在結(jié)構(gòu)完全一致（無中間注入元）的2行，若有則直接進(jìn)行Givens旋轉(zhuǎn)；

b.計算所有旋轉(zhuǎn)對產(chǎn)生的中間注入元個數(shù)，并對各旋轉(zhuǎn)對產(chǎn)生注入元的個數(shù)進(jìn)行比較，選擇產(chǎn)生注入元個數(shù)最少的2行作為確定的旋轉(zhuǎn)對進(jìn)行Givens旋轉(zhuǎn)；

c.在中間注入元個數(shù)一樣的情況下，選擇最稀疏的2行作為旋轉(zhuǎn)對進(jìn)行Givens旋轉(zhuǎn)。

下面通過簡單的4×3階矩陣A來演示用VPAIR策略進(jìn)行逐列變轉(zhuǎn)軸Givens旋轉(zhuǎn)的過程。原始矩陣A如式（11）所示。

其中，×表示非零元素。

步驟1：第1列共計有個旋轉(zhuǎn)對可供旋轉(zhuǎn)。旋轉(zhuǎn)對［R3，R4］無注入元，故選之進(jìn)行變換，元素（4，1）被消去，所得結(jié)果如式（12）所示。

其中，［Ri，Rj］表示第i行、第j行組成的旋轉(zhuǎn)對；（i，j）表示矩陣第i行第j列元素。

步驟2：第1列有3個旋轉(zhuǎn)對可供旋轉(zhuǎn)，旋轉(zhuǎn)對［R1，R2］和［R2，R3］均只有 1 個非零注入元，而旋轉(zhuǎn)對［R1，R2］更稀疏，故選之。此時元素（2，1）被消去，同時產(chǎn)生注入元（1，1）（注入元用?表示），所得結(jié)果如式（13）所示。

步驟 3：第1 列只剩下旋轉(zhuǎn)對［R1，R3］，故選之進(jìn)行變換。此時元素（3，1）被消去，同時產(chǎn)生注入元（1，2）。至此第1列旋轉(zhuǎn)結(jié)束。所得結(jié)果如式（14）所示。

步驟 4：進(jìn)行第2 列旋轉(zhuǎn)。旋轉(zhuǎn)對［R3，R4］無注入元，故直接選之進(jìn)行變換，元素（4，3）被消去。所得結(jié)果如式（15）所示。

步驟 5：第2列只剩下旋轉(zhuǎn)對［R2，R3］，故選之進(jìn)行變換。元素（3，2）被消去，產(chǎn)生注入元（2，2）。至此第2列旋轉(zhuǎn)結(jié)束。所得結(jié)果如式（16）所示。

步驟 6：第3列只剩下旋轉(zhuǎn)對［R3，R4］，故選之進(jìn)行變換。上三角矩陣R形成，旋轉(zhuǎn)過程結(jié)束。所得結(jié)果如式（17）所示。

由上述可知，VPAIR策略為了實(shí)現(xiàn)嚴(yán)格的矩陣上三角結(jié)構(gòu)，進(jìn)行了6次旋轉(zhuǎn)操作。

4 基于VPAIR消元的改進(jìn)策略

4.1 改進(jìn)策略的基本原理

VPAIR消元策略采用尋找產(chǎn)生最少注入元的旋轉(zhuǎn)對進(jìn)行變換，很大程度上減少了注入元的數(shù)量。但是，由4×3階矩陣A不難發(fā)現(xiàn)，當(dāng)對第i列進(jìn)行旋轉(zhuǎn)時，元素（i，i）出現(xiàn)了為0的情況。為了滿足上三角矩陣R的結(jié)構(gòu)特點(diǎn)，仍然要將該行納入旋轉(zhuǎn)對進(jìn)行變換，這在一定程度上增加了注入元素的數(shù)量，降低了計算效率。由于電力系統(tǒng)狀態(tài)估計的加權(quán)雅可比矩陣的大規(guī)模稀疏性，矩陣對角元素為0的情況更為普遍。對此，本文提出了基于VPAIR消元的改進(jìn)策略。仍然采用逐列消元，具體步驟如下。

a.對第i列進(jìn)行變換時，統(tǒng)計該列第i行以下（包括第i行）所有非零元素個數(shù)和非零元素所在行的信息。由非零元素所在行組成子矩陣。

b.對子矩陣進(jìn)行VPAIR策略的消元，即尋找非零注入元個數(shù)最少的旋轉(zhuǎn)對進(jìn)行旋轉(zhuǎn)。記錄該列最后一次旋轉(zhuǎn)的旋轉(zhuǎn)軸所在行信息。

c.將該旋轉(zhuǎn)軸所在行與第i行元素位置互換。d.進(jìn)行第i+1列的變換。

上述策略的關(guān)鍵在于旋轉(zhuǎn)過程中動態(tài)調(diào)整行排序（即量測排序），實(shí)現(xiàn)排序優(yōu)化。由式（4）可知：

由式（18）可知，加權(quán)雅可比矩陣 Hw經(jīng)Givens正交變換后得到的最終上三角矩陣R和信息矩陣平方根分解后得到的上三角矩陣相同。調(diào)整行排序的過程中雖然Hw發(fā)生變化，但是不會發(fā)生變化，所以最后的上三角矩陣R也是一致的。因此，改進(jìn)策略不改變計算結(jié)果，從理論上而言是可行的。

4.2 改進(jìn)策略的具體實(shí)現(xiàn)

為了說明改進(jìn)策略的變換過程，同樣采用上述4×3階矩陣A進(jìn)行演示。

步驟1：第1列有3個非零元素，共計個旋轉(zhuǎn)對可供旋轉(zhuǎn)。旋轉(zhuǎn)對［R3，R4］無注入元，故選之進(jìn)行變換，元素（4，1）被消去，所得結(jié)果如式（19）所示。

步驟2：A1第1列還有2個非零元素，組成唯一可供選擇的旋轉(zhuǎn)對［R2，R3］，故選之進(jìn)行變換，元素（3，1）被消去。此時第1列只剩下第2行一個非零元素，該列旋轉(zhuǎn)過程結(jié)束。將第2行元素與第1行元素位置互換，所得結(jié)果如式（20）所示。

步驟3：第2列還有2個非零元素，組成唯一可供選擇的旋轉(zhuǎn)對［R3，R4］，選之進(jìn)行變換，元素（4，1）被消去。此時第2列只剩下第3行一個非零元素，該列旋轉(zhuǎn)過程結(jié)束。將第3行元素與第2行元素（如備注所示為原第1行）位置互換，所得結(jié)果如式（21）所示。

步驟 4：對第3 列唯一旋轉(zhuǎn)對［R3，R4］（如備注所示為原［R1，R4］）進(jìn)行變換，形成上三角矩陣 R（A4）。

通過演示可以發(fā)現(xiàn)對于4×3階矩陣A而言，VPAIR策略需要進(jìn)行6次旋轉(zhuǎn)操作；改進(jìn)策略進(jìn)行了4次旋轉(zhuǎn)操作和2次矩陣行位置互換操作。由第3節(jié)Givens旋轉(zhuǎn)公式可知，每一次的旋轉(zhuǎn)操作是復(fù)雜的，運(yùn)行時間也遠(yuǎn)遠(yuǎn)大于位置互換操作的時間。因此，利用改進(jìn)策略的Givens變換可以在很大程度上減少計算時間，提高計算效率。

由式（22）的結(jié)果不難發(fā)現(xiàn)，改進(jìn)策略在得到最后上三角矩陣的同時，也給出了雅可比矩陣所在行（即量測量）的新排序。這種排序是在利用VPAIR策略進(jìn)行Givens變換的過程中得到的更加優(yōu)化的排序。該排序可以有效地規(guī)避旋轉(zhuǎn)過程中處理對角元元素為0的情況，很大程度上提高了計算效率。

4.3 策略比較與分析

與一般加權(quán)最小二乘法處理對稱矩陣不同，在正交變換法狀態(tài)估計中，Givens旋轉(zhuǎn)由于需要直接對加權(quán)雅可比矩陣進(jìn)行處理，所以存在量測排序的優(yōu)化問題。量測排序的合理與否，對旋轉(zhuǎn)過程中產(chǎn)生注入元的數(shù)量會產(chǎn)生重要的影響。

VPAIR策略在選擇最少注入元所在的行對進(jìn)行變換的過程中，本質(zhì)上也對量測排序進(jìn)行了一定程度的優(yōu)化。但是在逐列消元的過程中，該列最后一次Givens旋轉(zhuǎn)的旋轉(zhuǎn)軸所在行號與列號一致，是已經(jīng)固定的。所以這種量測排序的優(yōu)化只是半動態(tài)的優(yōu)化策略，部分量測量順序沒有得到優(yōu)化。

而改進(jìn)策略在逐列消元過程中，避免了固定最后旋轉(zhuǎn)軸所在行號的弊端，同時將最后的旋轉(zhuǎn)軸所在行與和該列列號一致的行號進(jìn)行在線動態(tài)的互換，調(diào)整了量測排序。這種量測排序的優(yōu)化手段是一種動態(tài)的優(yōu)化策略，因此具有更高的計算效率。

5 算例驗(yàn)證

為了驗(yàn)證所提改進(jìn)策略的計算效率，并與VPAIR策略進(jìn)行比較，2種策略均采用按照最小度原則的列排序，保證矩陣R稀疏性的同時，便于策略之間的比較。在IEEE 39、57、118節(jié)點(diǎn)系統(tǒng)和實(shí)際某區(qū)域23節(jié)點(diǎn)、區(qū)域433節(jié)點(diǎn)上分別進(jìn)行了測試。選用主頻為3.4 GHz、內(nèi)存為8 G的臺式機(jī)運(yùn)行程序。測試算例的配置信息如表1所示。

表1 測試系統(tǒng)的配置信息Table 1 Configuration information of test systems

將原策略和改進(jìn)策略分別在上述系統(tǒng)中進(jìn)行運(yùn)行和驗(yàn)證，計算結(jié)果如表2所示。其中旋轉(zhuǎn)次數(shù)為狀態(tài)估計中對加權(quán)雅可比矩陣Hw完成一次完整的正交變換，形成最終上三角矩陣R所進(jìn)行的Givens旋轉(zhuǎn)總次數(shù)，而CPU時間為完成上述正交變換的計算時間。

由表2可知，改進(jìn)策略在Givens旋轉(zhuǎn)總次數(shù)和計算時間上都有50%以上大幅度的減少。因此，改進(jìn)策略的計算效率得到很大的提高。

表2 VPAIR策略與改進(jìn)策略測試結(jié)果比較Table 2 Comparison of testing result between VPAIR strategy and improved strategy

6 結(jié)論

在電力系統(tǒng)狀態(tài)估計中，針對正交變換法在變換過程中出現(xiàn)大量非零注入元的情況，本文在VPAIR變轉(zhuǎn)軸逐列消元策略的基礎(chǔ)上，提出消元過程中動態(tài)調(diào)整量測排序的改進(jìn)策略。將所提改進(jìn)策略與VPAIR策略在不同規(guī)模仿真系統(tǒng)運(yùn)行比較，表明改進(jìn)策略具有更高的計算效率，滿足在線狀態(tài)估計的要求。

本文提出的改進(jìn)策略結(jié)合了量測排序的動態(tài)優(yōu)化，簡單有效，切實(shí)可行，能夠有效提高計算效率。

參考文獻(xiàn)：

［1］林偉芳，湯涌，孫華東，等.巴西“2·4”大停電事故及對電網(wǎng)安全穩(wěn)定運(yùn)行的啟示［J］.電力系統(tǒng)自動化，2011，35（9）：1-5.LIN Weifang，TANG Yong，SUN Huadong，etal.Blackoutin Brazilpowergrid on February 4，2011 and inspirationsfor stable operation of power grid［J］.Automation of Electric Power Systems，2011，35（9）：1-5.

［2］鐘慶，張哲，許中，等.廣州配電網(wǎng)故障停電事故的自組織臨界特征［J］.電力自動化設(shè)備，2017，37（4）：109-113.ZHONG Qing，ZHANG Zhe，XU Zhong，et al.SOC characteristics of power-supply failure of Guangzhou distribution network［J］.Electric Power Automation Equipment，2017，37（4）：109-113.

［3］倪小平，張步涵.一種帶有等式約束的狀態(tài)估計新算法［J］.電力系統(tǒng)自動化，2001，25（21）：42-44.NIXiaoping，ZHANG Buhan.New stateestimation algorithm with equality constraints［J］.Automation of Electric Power Systems，2001，25（21）：42-44.

［4］KORRES G N.A robustalgorithm forpowersystem state estimation with equality constraints［J］.IEEE Transactions on Power Systems，2010，25（3）：1531-1541.

［5］郭燁，張伯明，吳文傳，等.直角坐標(biāo)下含零注入約束的電力系統(tǒng)狀態(tài)估計修正牛頓法［J］.中國電機(jī)工程學(xué)報，2012，32（19）：96-100.GUO Ye，ZHANG Boming，WU Wenchuan，et al.Power system state estimation solution with zero injection constraints using modified Newton method［J］.Proceedings of the CSEE，2012，32（19）：96-100.

［6］郭燁，張伯明，吳文傳，等.極坐標(biāo)下含零注入約束的電力系統(tǒng)狀態(tài)估計的修正牛頓法與快速解耦估計［J］.中國電機(jī)工程學(xué)報，2012，32（19）：113-117.GUO Ye，ZHANG Boming，WU Wenchuan，et al.Power system state estimation solution with zero injection constraints using modified Newton method and fast decoupled method in polar coordinate［J］.Proceedings of the CSEE，2012，32（19）：113-117.

［7］常鮮戎，樊瑞.計及零注入節(jié)點(diǎn)約束的混合量測分區(qū)狀態(tài)估計方法［J］.電網(wǎng)技術(shù)，2015，39（8）：2254-2257.CHANG Xianrong，F(xiàn)AN Rui.A mixedmeasurementpartition state estimation method taking zero injection node constraint into account［J］.Power System Technology，2015，39 （8）：2254-2257.

［8］劉廣一，胡錫龍，于爾鏗，等.快速正交變換阻尼最小二乘法在電力系統(tǒng)狀態(tài)估計中的應(yīng)用［J］.中國電機(jī)工程學(xué)報，1991，11（6）：34-40.LIU Guangyi，HU Xilong，YU Erkeng，et al.Application of fast orthogonal transformation with damping factor to power system state estimation［J］.Proceedings of the CSEE，1991，11（6）：34-40.

［9］VEMPATI N，SLUTSKER I W，TINNEY W F.Orthogonal sparse vector methods［J］.IEEE Transactions on Power Systems，1992，7（2）：926-932.

［10］DUFF I S.Pivot selection and row ordering in Givens reduction on sparse matrices［J］.Computing，1974，13（3-4）：239-248.

［11］GENTLEMAN W M.Row elimination for solving sparse linear systems and least squares problems［J］.Springer Berlin Heidelberg，1976，506：122-133.

［12］ROBEY T H，SULSKY D L.Row ordering for a sparse QR decomposition［J］.SIAM Journal of Matrix Analysis and Applications，1994，15（4）：1208-1225.

［13］PANDIAN A，PARTHASARATHY K，SOMAN S A.Towards fast givens rotations based power system state estimator［J］.IEEE Transactions on Power Systems，1999，14（3）：837-843.

［14］杜正春，牛振勇，方萬良.基于分塊QR分解的一種狀態(tài)估計算法［J］.中國電機(jī)工程學(xué)報，2003，23（8）：50-55.DU Zhengchun，NIU Zhenyong，F(xiàn)ANG Wanliang.A block QR based power system state estimation algorithm［J］.Proceedings of the CSEE，2003，23（8）：50-55.

［15］郭瑞鵬，邵學(xué)儉，韓禎祥.基于分塊吉文斯旋轉(zhuǎn)的電力系統(tǒng)狀態(tài)估計［J］.中國電機(jī)工程學(xué)報，2006，26（12）：26-31.GUO Ruipeng，SHAO Xuejian，HAN Zhenxiang. A blocked givens rotations based algorithm for power system state estimation［J］.Proceedings of the CSEE，2006，26（12）：26-31.

［16］WANG G，ALI K，XU J.An efficient engine for orthogonal SE［C］∥Proceedingsofthe Third InternationalConference on Control Automation&Systems Engineering.［S.l.］：CASE，2013：13-17.