電力系統(tǒng)狀態(tài)估計(jì)中稀疏技術(shù)的應(yīng)用

0基于潮流計(jì)算的稀疏技術(shù)應(yīng)用研究

在計(jì)算電網(wǎng)的過程中，有大量的矩陣、矩陣和向量之間的操作過程。由于電網(wǎng)本身的連接

線特點(diǎn)，這些矩陣和向量通常只有少量的非零元素。該矩陣和向量稱為稀疏矩陣和稀疏向

量。稀疏技術(shù)就是利用矩陣或者向量的非零元分布特點(diǎn),從排零存儲(chǔ)和排零運(yùn)算2個(gè)方面

來提高運(yùn)算效率。

稀疏技術(shù)在電力系統(tǒng)計(jì)算中已經(jīng)得到廣泛應(yīng)用。1967年,美國學(xué)者W.F.Tinney等人將稀疏

矩陣技術(shù)應(yīng)用到潮流計(jì)算,大幅度提高了潮流計(jì)算速度。1985年,Tinney等人又在稀疏矩

陣技術(shù)的基礎(chǔ)上系統(tǒng)提出了稀疏向量法。當(dāng)線性方程組右端項(xiàng)只存在少量非零元或者只需

求解待求向量中少量分量時(shí)，應(yīng)用稀疏向量法可以避免前代回代運(yùn)算過程中無效的運(yùn)算操

作,從而進(jìn)一步提高求解大規(guī)模電網(wǎng)應(yīng)用問題的計(jì)算效率。在我國，對稀疏技術(shù)的研究也比

較多，文獻(xiàn)重點(diǎn)研究了稀疏矩陣的鏈表存儲(chǔ)技術(shù)及其在潮流計(jì)算中的應(yīng)用，文獻(xiàn)研究了稀疏

向量法在電網(wǎng)計(jì)算中的應(yīng)用。目前稀疏技術(shù)已廣泛應(yīng)用于電力系統(tǒng)分析計(jì)算的各個(gè)領(lǐng)域，

如潮流計(jì)算、靜態(tài)安全分析、電網(wǎng)故障計(jì)算等。

狀態(tài)估計(jì)是能量管理系統(tǒng)(EMS)的核心功能,是EMS高級應(yīng)用軟件實(shí)用化的基礎(chǔ)。與潮流計(jì)

算相比，由于狀態(tài)估計(jì)處理的量測數(shù)目眾多，種類多樣,因此雅可比矩陣規(guī)模更大,運(yùn)算更為

復(fù)雜，不僅要應(yīng)用到稀疏矩陣的存儲(chǔ)技術(shù),還涉及矩陣的相乘、轉(zhuǎn)置等操作。因此,在狀態(tài)

估計(jì)中提高稀疏技術(shù)應(yīng)用技巧顯得尤為重要。本文首先根據(jù)因子矩陣和消去樹,討論了稀

疏向量法。然后基于給定的稀疏矩陣存儲(chǔ)方法和符號因子化技術(shù)，提出了一種應(yīng)用稀疏向

量法進(jìn)行LDLT分解的算法,最后結(jié)合量測殘差方差計(jì)算給出了2種應(yīng)用稀疏技術(shù)計(jì)算的方

法。為驗(yàn)證算法有效性，木文使用IEEE118和IEEE300系統(tǒng)算例進(jìn)行了計(jì)算，給出了測試

結(jié)果。

1節(jié)點(diǎn)狀態(tài)估計(jì)

給定量測模型

在給定網(wǎng)絡(luò)接線、支路參數(shù)的條件下，狀態(tài)估計(jì)問題是求取使得如下目標(biāo)函數(shù)達(dá)到最小的x

的值

式(1)(2)中：z是量測矢量,v是量測誤差,m維;x是節(jié)點(diǎn)狀態(tài)量矢量,n維,一般情況下

m2n;h(x)為非線性量測方程式矢量，包括支路潮流方程式、節(jié)點(diǎn)潮流注入方程式等;RT是

量測權(quán)重對角陣，是每個(gè)量測方差。2倒數(shù)的對角陣。

采用牛頓法求解上述方程,狀態(tài)估計(jì)計(jì)算過程為反復(fù)迭代求解方程式(3)并修正狀態(tài)量x,直

至滿足指定的收斂判據(jù)。

其中:H為量測雅可比矩陣;G稱為信息矩陣；Az為量測殘差矢量；Ax為狀態(tài)量修正量。

2稀疏技術(shù)介紹

2.1基于鏈表的存儲(chǔ)方法

本文涉及的稀疏技術(shù)主要包含2方面的內(nèi)容:一是稀疏矩陣存儲(chǔ)技術(shù),主要研究稀疏矩陣的

存儲(chǔ)策略、索引技巧等，力求節(jié)省內(nèi)存的同時(shí)能夠快速定位矩陣元素;二是稀疏向量技術(shù)，

主要研究在前代回代運(yùn)算中如何避免不必要的運(yùn)算操作，提高計(jì)算效率。

?般來講，有2類稀疏矩除的存儲(chǔ)方法，?類是基于數(shù)組的傳統(tǒng)存儲(chǔ)方法。這是目前大多數(shù)

計(jì)算程序特別是Fortran程序所采用的方法;另一類是基于鏈表的存儲(chǔ)方法。文獻(xiàn)比較了2

種存儲(chǔ)方法的效率，認(rèn)為基于鏈表的存儲(chǔ)方法雖然多占用了少量的內(nèi)存，但計(jì)算效率高于傳

統(tǒng)存儲(chǔ)方法。本文采用基于數(shù)組的存儲(chǔ)方式，但在數(shù)組上根據(jù)計(jì)算需要加鏈,提高稀疏矩陣

元素的杳找效率。和基于鏈表存儲(chǔ)方法相比,這種存儲(chǔ)方法只需對傳統(tǒng)程序作較少的修改，

查找效率不變,并且更為節(jié)省內(nèi)存。

以導(dǎo)納矩陣的存儲(chǔ)為例，當(dāng)系統(tǒng)中不含移相器支路時(shí),導(dǎo)納矩陣是對稱矩陣,所以只需存儲(chǔ)

對角元和上二角陣邪零元、另外,結(jié)合狀態(tài)估計(jì)計(jì)算過程，雅可比矩陣H可按行形成，并且

量測殘差A(yù)z計(jì)算只需按行取導(dǎo)納矩陣的一行,所以在存儲(chǔ)每行非零元起始位置、非零元

列號及值的同時(shí),再存儲(chǔ)列方向上的行號縱向鏈和位置縱向鏈(簡稱列鏈)，可以保證存取導(dǎo)

納矩陣的效率。其存儲(chǔ)結(jié)構(gòu)如下:DIAG為存儲(chǔ)每行對角元位置的一維數(shù)組;COL為存儲(chǔ)每個(gè)

非零元列號;LINKN為列方向的行號縱向鏈;LINKP為列方向的位置縱向鏈;G、B為導(dǎo)納值的

一維數(shù)組。

為了減少內(nèi)存占用，狀態(tài)估計(jì)過程中有些矩陣元素值可以不必存儲(chǔ)，例如H的轉(zhuǎn)置陣HT以

及信息矩陣Go由于這些矩陣可以由H陣轉(zhuǎn)置或者HTR-1H相乘獲得，因此只需要在迭弋前

形成HT的索引并對HTH必講行符號定位，等到真正需要這些元素參與計(jì)算時(shí)，再計(jì)算巴官

們的值。

2.2快速前代和快速回代運(yùn)算

對于式⑶所示的方程式,記右端項(xiàng)為列向量b=HTR-lAz,通過因子分解可將信息矩陣G表

示為

其中:U陣為對角元為1的上三角矩陣;D為n維對角陣。定義Lk為U陣中的保留第k行元

素而其他對角元為1的n階因子矩陣,即

式中:砒是行向量，為U陣中除去對角元后的第k行；Ik是列向量,為單位陣的第k歹人不

難驗(yàn)證UkT=ITkmk,那么U陣及其逆陣可以表示成一序列因子矩陣的乘積,如式(6)

因此，求解方程(3)可分為3步:

其中：式(7)稱為前代運(yùn)算;式(8)稱為規(guī)格化;式(9)稱為回代運(yùn)算。當(dāng)右端項(xiàng)向量b中僅有

少量非零元時(shí);如果求解式(7)只在非零元的消去樹路徑?jīng)Q定的因子矩陣上進(jìn)行時(shí),可以避

免零元參與運(yùn)算，此時(shí)稱為快速前代運(yùn)算；同樣，在求解式(9)的過程中如果我們僅對Ax中

的少量分量感興趣,就只需在這幾個(gè)分量的消去樹路徑?jīng)Q定的因子矩陣上運(yùn)算,這樣可以避

免不必要的運(yùn)算操作,此時(shí)稱為快速回代運(yùn)算。

例如,對于圖1所示由信息矩陣G確定的網(wǎng)絡(luò),其消去樹如圖2所示(圖1中虛線表示因子

分解過程中增加的注入元)。圖2中，邊的方向是由小號指向大號,節(jié)點(diǎn)k代表U陣中k行

確定的因子矩陣Uk,節(jié)點(diǎn)9是消去樹的根節(jié)點(diǎn)，代表單位矩陣。除根節(jié)點(diǎn)外,每個(gè)節(jié)點(diǎn)都有

一個(gè)父節(jié)點(diǎn),定義為對應(yīng)因子矩陣第一個(gè)非零非對角元的列號。每個(gè)節(jié)點(diǎn)的消去樹路徑是

指包含該節(jié)點(diǎn)直到樹根上的所有節(jié)點(diǎn)的集合，如節(jié)點(diǎn)5的消去樹路徑指的是｛5,7,8,9｝。如

果右端項(xiàng)中只有非零元bl,根據(jù)式(7)在計(jì)算y向量時(shí),我們只需在bl路徑確定的因子矩陣

UI、U4、U7、U8、U9上進(jìn)行快速前代運(yùn)算，即

同樣如果我們只對Ax向量中的Ax5感興趣，根據(jù)式(9),我們只需在U9、U8、U7、U5±

進(jìn)行快速回代運(yùn)算，即

實(shí)際計(jì)算過程中，稀疏向量b中可能有多個(gè)非零元，或者對'Ax向量中幾個(gè)分量感興趣，這

時(shí)參與快速前代或快速回代運(yùn)算的因了矩陣是這幾個(gè)元素消去樹路徑的合集。如對bl和

b5確定的消去樹路徑是UI、U4、U5、U7、U8、U9,其中節(jié)點(diǎn)7稱為分叉點(diǎn)。快速前代運(yùn)算

時(shí)必須保證集合中比分叉點(diǎn)7小的因子矩陣按順序運(yùn)算完成后才能使U7參與運(yùn)算；同樣，

快速回代運(yùn)算時(shí)必須保證比分叉點(diǎn)大的因子矩陣按順序完成后方可使U7參與運(yùn)算。

由此可見,提高稀疏向量法效率的措施可以從2方面著手:一方面盡量減少U陣中元素，使

參與計(jì)算的元素盡可能的少，即減少因子分解過程中產(chǎn)生的注入元；另一方面可以盡可能減

少參與計(jì)算的因子矩陣個(gè)數(shù),也就是消去樹路徑的長度來減少總的參與計(jì)算的元素。計(jì)算

前的節(jié)點(diǎn)優(yōu)化編號正是從這2方面來提高稀疏向量法的效率。

3稀疏技術(shù)在狀態(tài)評估中的應(yīng)用

3.1方程的求解

3.1.1基于Idlt因子表的求解策略

采用直接法求解法方程式(3)時(shí),可以采用LDLT因子分解法,也可以采用正交分解法或者其

他方法。正交分解法數(shù)值穩(wěn)定性好，但計(jì)算量很大,并且編程復(fù)雜,在需要多次分解G陣時(shí)

計(jì)算效率較低;LDLT分解法當(dāng)量測權(quán)重:相差懸殊時(shí)數(shù)值穩(wěn)定性較差，但其計(jì)算量比正交分解

法小得多，是電力系統(tǒng)計(jì)算中對稱矩陣分解常用的算法。本節(jié)研究了LDLT因子分解技術(shù)，

包括存儲(chǔ)策略分析、G陣和I：陣稀疏結(jié)構(gòu)的形成以及采用稀疏向量法求解LDLT因子表等。

LDLT因子分解的計(jì)算公式如下

本文提出一種應(yīng)用稀疏向量法進(jìn)行LDLT分解的算法。假設(shè)分解進(jìn)行到了第k步,設(shè)Uk-1

表示U陣的第k-1階子式:U(1:k-l,k)表示U陣第k列不包括對角元的列向量,簡記為v,相

應(yīng)的對于G陣有Gk-1子區(qū)、G(l:k-1,k)(簡記為w)以及gkk,對于D陣有Dk-1和dkk,此時(shí),

如下等式成立：

即

所以有如下等式(12)成立

因此LDLT分解的算法可以描述如下：

第一步：由UkT-ly=w采用稀疏向量法求y;

第二步:vDTk-ly;

第三步：dkk=gkk-vTy;

Endfor

其中第一步的計(jì)算量最大:其他步相對第一步而言很小。一般情況下,Uk-1.w和y是稀疏

矩陣或者稀疏向量,為求出y,必須先根據(jù)w的非零元決定y的非零分布。事實(shí)上,w和y的

非零分布分別是G陣和U陣的列鏈描述的非零分布,其形成過程在3.1.2節(jié)和3.1.3節(jié)中

己有敘述。所以,在每次用稀疏向量法求解y的過程中，無需形成Uk-1子矩陣的消去樹，而

且計(jì)算所得的結(jié)果向量v可按列鏈順序存入U(xiǎn)陣相應(yīng)位置。

3.1.5算法效率較低

為測試采用稀疏向量法進(jìn)行LDLT分解的效率，本文實(shí)現(xiàn)了幾種常用的LDLT分解算法，各種

算法簡單描述如下：

算法1：文獻(xiàn)中按行LDLT分解算法。該算法根據(jù)G陣元素初始化工作行，由于在存儲(chǔ)結(jié)構(gòu)h

沒有形成U陣的列鏈和初始化的工作行沒有為注入元預(yù)留位置,所以不得不在U陣列方向

上搜索以及在U陣行方向上向工作行插入注入元,從而導(dǎo)致算法效率較低。

算法2:仍然采用文獻(xiàn)中按行LDLT分解方法,但根據(jù)3.1,2節(jié)和3.1.3節(jié)所述形成了U陣和

G陣的列鏈。由于列鏈的作用，可以避免在列方向上搜索U陣元素；同時(shí)由于f陣結(jié)構(gòu)已知，

因而可以根據(jù)U陣結(jié)構(gòu)分配工作行和根據(jù)G陣值初始化工作行,這樣注入元的位置可以在

工作行上預(yù)留出來。

算法3:基于3.1.2和3.1,3節(jié)形成的G陣和U陣結(jié)構(gòu)，采用稀疏向量法求解。由于每次稀

疏向量求解元素的位置和順序都是確定的，無需形成工作行,因而計(jì)算效率得以提高。

表1是在IEEE118和IEEE300系統(tǒng)上測試上述各種算法效率的計(jì)算時(shí)間結(jié)果。兩系統(tǒng)量

測分別取1098個(gè)和2544個(gè),節(jié)點(diǎn)優(yōu)化編號方法采用最小度最少前趨節(jié)點(diǎn)數(shù)法(MD-MNP)。

表中時(shí)間為計(jì)算10000次的平均時(shí)間，單位ms。測試結(jié)果表明，系統(tǒng)規(guī)模越大，相對而言算

法3更具有優(yōu)勢。

3.2確定殘差方差

使用標(biāo)準(zhǔn)化殘差搜索法或者遞歸量測誤差估計(jì)辨識法時(shí)，需要計(jì)算量測的方差。據(jù)文獻(xiàn),殘

差E的協(xié)方差矩陣為

式(13)WR陣的對角元即對應(yīng)量測殘差方差為

求解式(14)可以采用不同的稀疏技術(shù)。

(1)方法io

將式(14)寫成分量形式：

由于hk是稀疏向量,因此可以采用快速前代技術(shù)先求出向量U-Thk,然后作向量的乘積即可

求得對應(yīng)的殘差方差。

⑵方法2。

利用狀態(tài)估計(jì)收斂后已有的因子表計(jì)算信息矩陣G的逆陣,然后做矩陣與向量的相乘，可以

求得殘差方差dk.據(jù)式(4)

進(jìn)而可以寫成

根據(jù)式(5)(6),UkT可以直接由Uk獲得,所以G-1的求取只需要一序列因子矩陣乘積即可。

并且由于G陣是對角占優(yōu)的對稱正定陣，可以利用G陣的“稀疏逆”近似替代GT陣，即只

計(jì)算GT上三角陣中與U陣分布相同的非零元。令OGT,則C中元素可通過下列公式求解:

實(shí)際計(jì)算效果表明，對于GT沒有參與計(jì)算的元素，其值非常小，忽略其影響對量測殘差方

差計(jì)算影響極小。

(3)兩種計(jì)算方法的比較。

以3.1.5中采用