电力系统状态估计中稀疏技术的应用

电力系统状态估计中稀疏技术的应用

0基于潮流计算的稀疏技术应用研究

在计算电网的过程中，有大量的矩阵、矩阵和向量之间的操作过程。由于电网本身的连接

线特点，这些矩阵和向量通常只有少量的非零元素。该矩阵和向量称为稀疏矩阵和稀疏向

量。稀疏技术就是利用矩阵或者向量的非零元分布特点,从排零存储和排零运算2个方面

来提高运算效率。

稀疏技术在电力系统计算中已经得到广泛应用。1967年,美国学者W.F.Tinney等人将稀疏

矩阵技术应用到潮流计算,大幅度提高了潮流计算速度。1985年,Tinney等人又在稀疏矩

阵技术的基础上系统提出了稀疏向量法。当线性方程组右端项只存在少量非零元或者只需

求解待求向量中少量分量时，应用稀疏向量法可以避免前代回代运算过程中无效的运算操

作,从而进一步提高求解大规模电网应用问题的计算效率。在我国，对稀疏技术的研究也比

较多，文献重点研究了稀疏矩阵的链表存储技术及其在潮流计算中的应用，文献研究了稀疏

向量法在电网计算中的应用。目前稀疏技术已广泛应用于电力系统分析计算的各个领域，

如潮流计算、静态安全分析、电网故障计算等。

状态估计是能量管理系统(EMS)的核心功能,是EMS高级应用软件实用化的基础。与潮流计

算相比，由于状态估计处理的量测数目众多，种类多样,因此雅可比矩阵规模更大,运算更为

复杂，不仅要应用到稀疏矩阵的存储技术,还涉及矩阵的相乘、转置等操作。因此,在状态

估计中提高稀疏技术应用技巧显得尤为重要。本文首先根据因子矩阵和消去树,讨论了稀

疏向量法。然后基于给定的稀疏矩阵存储方法和符号因子化技术，提出了一种应用稀疏向

量法进行LDLT分解的算法,最后结合量测残差方差计算给出了2种应用稀疏技术计算的方

法。为验证算法有效性，木文使用IEEE118和IEEE300系统算例进行了计算，给出了测试

结果。

1节点状态估计

给定量测模型

在给定网络接线、支路参数的条件下，状态估计问题是求取使得如下目标函数达到最小的x

的值

式(1)(2)中：z是量测矢量,v是量测误差,m维;x是节点状态量矢量,n维,一般情况下

m2n;h(x)为非线性量测方程式矢量，包括支路潮流方程式、节点潮流注入方程式等;RT是

量测权重对角阵，是每个量测方差。2倒数的对角阵。

采用牛顿法求解上述方程,状态估计计算过程为反复迭代求解方程式(3)并修正状态量x,直

至满足指定的收敛判据。

其中:H为量测雅可比矩阵;G称为信息矩阵；Az为量测残差矢量；Ax为状态量修正量。

2稀疏技术介绍

2.1基于链表的存储方法

本文涉及的稀疏技术主要包含2方面的内容:一是稀疏矩阵存储技术,主要研究稀疏矩阵的

存储策略、索引技巧等，力求节省内存的同时能够快速定位矩阵元素;二是稀疏向量技术，

主要研究在前代回代运算中如何避免不必要的运算操作，提高计算效率。

•般来讲，有2类稀疏矩除的存储方法，•类是基于数组的传统存储方法。这是目前大多数

计算程序特别是Fortran程序所采用的方法;另一类是基于链表的存储方法。文献比较了2

种存储方法的效率，认为基于链表的存储方法虽然多占用了少量的内存，但计算效率高于传

统存储方法。本文采用基于数组的存储方式，但在数组上根据计算需要加链,提高稀疏矩阵

元素的杳找效率。和基于链表存储方法相比,这种存储方法只需对传统程序作较少的修改，

查找效率不变,并且更为节省内存。

以导纳矩阵的存储为例，当系统中不含移相器支路时,导纳矩阵是对称矩阵,所以只需存储

对角元和上二角阵邪零元、另外,结合状态估计计算过程，雅可比矩阵H可按行形成，并且

量测残差Az计算只需按行取导纳矩阵的一行,所以在存储每行非零元起始位置、非零元

列号及值的同时,再存储列方向上的行号纵向链和位置纵向链(简称列链)，可以保证存取导

纳矩阵的效率。其存储结构如下:DIAG为存储每行对角元位置的一维数组;COL为存储每个

非零元列号;LINKN为列方向的行号纵向链;LINKP为列方向的位置纵向链;G、B为导纳值的

一维数组。

为了减少内存占用，状态估计过程中有些矩阵元素值可以不必存储，例如H的转置阵HT以

及信息矩阵Go由于这些矩阵可以由H阵转置或者HTR-1H相乘获得，因此只需要在迭弋前

形成HT的索引并对HTH必讲行符号定位，等到真正需要这些元素参与计算时，再计算巴官

们的值。

2.2快速前代和快速回代运算

对于式⑶所示的方程式,记右端项为列向量b=HTR-lAz,通过因子分解可将信息矩阵G表

示为

其中:U阵为对角元为1的上三角矩阵;D为n维对角阵。定义Lk为U阵中的保留第k行元

素而其他对角元为1的n阶因子矩阵,即

式中:砒是行向量，为U阵中除去对角元后的第k行；Ik是列向量,为单位阵的第k歹人不

难验证UkT=ITkmk,那么U阵及其逆阵可以表示成一序列因子矩阵的乘积,如式(6)

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

其中：式(7)称为前代运算;式(8)称为规格化;式(9)称为回代运算。当右端项向量b中仅有

少量非零元时;如果求解式(7)只在非零元的消去树路径决定的因子矩阵上进行时,可以避

免零元参与运算，此时称为快速前代运算；同样，在求解式(9)的过程中如果我们仅对Ax中

的少量分量感兴趣,就只需在这几个分量的消去树路径决定的因子矩阵上运算,这样可以避

免不必要的运算操作,此时称为快速回代运算。

例如,对于图1所示由信息矩阵G确定的网络,其消去树如图2所示(图1中虚线表示因子

分解过程中增加的注入元)。图2中，边的方向是由小号指向大号,节点k代表U阵中k行

确定的因子矩阵Uk,节点9是消去树的根节点，代表单位矩阵。除根节点外,每个节点都有

一个父节点,定义为对应因子矩阵第一个非零非对角元的列号。每个节点的消去树路径是

指包含该节点直到树根上的所有节点的集合，如节点5的消去树路径指的是｛5,7,8,9｝。如

果右端项中只有非零元bl,根据式(7)在计算y向量时,我们只需在bl路径确定的因子矩阵

UI、U4、U7、U8、U9上进行快速前代运算，即

同样如果我们只对Ax向量中的Ax5感兴趣，根据式(9),我们只需在U9、U8、U7、U5±

进行快速回代运算，即

实际计算过程中，稀疏向量b中可能有多个非零元，或者对'Ax向量中几个分量感兴趣，这

时参与快速前代或快速回代运算的因了矩阵是这几个元素消去树路径的合集。如对bl和

b5确定的消去树路径是UI、U4、U5、U7、U8、U9,其中节点7称为分叉点。快速前代运算

时必须保证集合中比分叉点7小的因子矩阵按顺序运算完成后才能使U7参与运算；同样，

快速回代运算时必须保证比分叉点大的因子矩阵按顺序完成后方可使U7参与运算。

由此可见,提高稀疏向量法效率的措施可以从2方面着手:一方面尽量减少U阵中元素，使

参与计算的元素尽可能的少，即减少因子分解过程中产生的注入元；另一方面可以尽可能减

少参与计算的因子矩阵个数,也就是消去树路径的长度来减少总的参与计算的元素。计算

前的节点优化编号正是从这2方面来提高稀疏向量法的效率。

3稀疏技术在状态评估中的应用

3.1方程的求解

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

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

他方法。正交分解法数值稳定性好，但计算量很大,并且编程复杂,在需要多次分解G阵时

计算效率较低;LDLT分解法当量测权重:相差悬殊时数值稳定性较差，但其计算量比正交分解

法小得多，是电力系统计算中对称矩阵分解常用的算法。本节研究了LDLT因子分解技术，

包括存储策略分析、G阵和I：阵稀疏结构的形成以及采用稀疏向量法求解LDLT因子表等。

LDLT因子分解的计算公式如下

本文提出一种应用稀疏向量法进行LDLT分解的算法。假设分解进行到了第k步,设Uk-1

表示U阵的第k-1阶子式:U(1:k-l,k)表示U阵第k列不包括对角元的列向量,简记为v,相

应的对于G阵有Gk-1子区、G(l:k-1,k)(简记为w)以及gkk,对于D阵有Dk-1和dkk,此时,

如下等式成立：

即

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

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

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

第二步:vDTk-ly;

第三步：dkk=gkk-vTy;

Endfor

其中第一步的计算量最大:其他步相对第一步而言很小。一般情况下,Uk-1.w和y是稀疏

矩阵或者稀疏向量,为求出y,必须先根据w的非零元决定y的非零分布。事实上,w和y的

非零分布分别是G阵和U阵的列链描述的非零分布,其形成过程在3.1.2节和3.1.3节中

己有叙述。所以,在每次用稀疏向量法求解y的过程中，无需形成Uk-1子矩阵的消去树，而

且计算所得的结果向量v可按列链顺序存入U阵相应位置。

3.1.5算法效率较低

为测试采用稀疏向量法进行LDLT分解的效率，本文实现了几种常用的LDLT分解算法，各种

算法简单描述如下：

算法1：文献中按行LDLT分解算法。该算法根据G阵元素初始化工作行，由于在存储结构h

没有形成U阵的列链和初始化的工作行没有为注入元预留位置,所以不得不在U阵列方向

上搜索以及在U阵行方向上向工作行插入注入元,从而导致算法效率较低。

算法2:仍然采用文献中按行LDLT分解方法,但根据3.1,2节和3.1.3节所述形成了U阵和

G阵的列链。由于列链的作用，可以避免在列方向上搜索U阵元素；同时由于f阵结构已知，

因而可以根据U阵结构分配工作行和根据G阵值初始化工作行,这样注入元的位置可以在

工作行上预留出来。

算法3:基于3.1.2和3.1,3节形成的G阵和U阵结构，采用稀疏向量法求解。由于每次稀

疏向量求解元素的位置和顺序都是确定的，无需形成工作行,因而计算效率得以提高。

表1是在IEEE118和IEEE300系统上测试上述各种算法效率的计算时间结果。两系统量

测分别取1098个和2544个,节点优化编号方法采用最小度最少前趋节点数法(MD-MNP)。

表中时间为计算10000次的平均时间，单位ms。测试结果表明，系统规模越大，相对而言算

法3更具有优势。

3.2确定残差方差

使用标准化残差搜索法或者递归量测误差估计辨识法时，需要计算量测的方差。据文献,残

差E的协方差矩阵为

式(13)WR阵的对角元即对应量测残差方差为

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

(1)方法io

将式(14)写成分量形式：

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

求得对应的残差方差。

⑵方法2。

利用状态估计收敛后已有的因子表计算信息矩阵G的逆阵,然后做矩阵与向量的相乘，可以

求得残差方差dk.据式(4)

进而可以写成

根据式(5)(6),UkT可以直接由Uk获得,所以G-1的求取只需要一序列因子矩阵乘积即可。

并且由于G阵是对角占优的对称正定阵，可以利用G阵的“稀疏逆”近似替代GT阵，即只

计算GT上三角阵中与U阵分布相同的非零元。令OGT,则C中元素可通过下列公式求解:

实际计算效果表明，对于GT没有参与计算的元素，其值非常小，忽略其影响对量测残差方

差计算影响极小。

(3)两种计算方法的比较。

以3.1.5中采用