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

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

0基于潮流计算的稀疏技术应用研究在计算电网的过程中，有大量的矩阵、矩阵和向量之间的操作过程。由于电网本身的连接线特点，这些矩阵和向量通常只有少量的非零元素。该矩阵和向量称为稀疏矩阵和稀疏向量。稀疏技术就是利用矩阵或者向量的非零元分布特点,从排零存储和排零运算2个方面来提高运算效率。稀疏技术在电力系统计算中已经得到广泛应用。1967年,美国学者W.F.Tinney等人将稀疏矩阵技术应用到潮流计算,大幅度提高了潮流计算速度。1985年,Tinney等人又在稀疏矩阵技术的基础上系统提出了稀疏向量法。当线性方程组右端项只存在少量非零元或者只需求解待求向量中少量分量时,应用稀疏向量法可以避免前代回代运算过程中无效的运算操作,从而进一步提高求解大规模电网应用问题的计算效率。在我国,对稀疏技术的研究也比较多,文献重点研究了稀疏矩阵的链表存储技术及其在潮流计算中的应用,文献研究了稀疏向量法在电网计算中的应用。目前稀疏技术已广泛应用于电力系统分析计算的各个领域,如潮流计算、静态安全分析、电网故障计算等。状态估计是能量管理系统(EMS)的核心功能,是EMS高级应用软件实用化的基础。与潮流计算相比,由于状态估计处理的量测数目众多,种类多样,因此雅可比矩阵规模更大,运算更为复杂,不仅要应用到稀疏矩阵的存储技术,还涉及矩阵的相乘、转置等操作。因此,在状态估计中提高稀疏技术应用技巧显得尤为重要。本文首先根据因子矩阵和消去树,讨论了稀疏向量法。然后基于给定的稀疏矩阵存储方法和符号因子化技术,提出了一种应用稀疏向量法进行LDLT分解的算法,最后结合量测残差方差计算给出了2种应用稀疏技术计算的方法。为验证算法有效性,本文使用IEEE118和IEEE300系统算例进行了计算,给出了测试结果。1节点状态估计给定量测模型在给定网络接线、支路参数的条件下,状态估计问题是求取使得如下目标函数达到最小的x的值式(1)(2)中:z是量测矢量,v是量测误差,m维;x是节点状态量矢量,n维,一般情况下m≥n;h(x)为非线性量测方程式矢量,包括支路潮流方程式、节点潮流注入方程式等;R-1是量测权重对角阵,是每个量测方差σ2倒数的对角阵。采用牛顿法求解上述方程,状态估计计算过程为反复迭代求解方程式(3)并修正状态量x,直至满足指定的收敛判据。其中:H为量测雅可比矩阵;G称为信息矩阵;Δz为量测残差矢量;Δx为状态量修正量。2稀疏技术介绍2.1基于链表的存储方法本文涉及的稀疏技术主要包含2方面的内容:一是稀疏矩阵存储技术,主要研究稀疏矩阵的存储策略、索引技巧等,力求节省内存的同时能够快速定位矩阵元素;二是稀疏向量技术,主要研究在前代回代运算中如何避免不必要的运算操作,提高计算效率。一般来讲,有2类稀疏矩阵的存储方法,一类是基于数组的传统存储方法。这是目前大多数计算程序特别是Fortran程序所采用的方法;另一类是基于链表的存储方法。文献比较了2种存储方法的效率,认为基于链表的存储方法虽然多占用了少量的内存,但计算效率高于传统存储方法。本文采用基于数组的存储方式,但在数组上根据计算需要加链,提高稀疏矩阵元素的查找效率。和基于链表存储方法相比,这种存储方法只需对传统程序作较少的修改,查找效率不变,并且更为节省内存。以导纳矩阵的存储为例,当系统中不含移相器支路时,导纳矩阵是对称矩阵,所以只需存储对角元和上三角阵非零元。另外,结合状态估计计算过程,雅可比矩阵H可按行形成,并且量测残差Δz计算只需按行取导纳矩阵的一行,所以在存储每行非零元起始位置、非零元列号及值的同时,再存储列方向上的行号纵向链和位置纵向链(简称列链),可以保证存取导纳矩阵的效率。其存储结构如下:DIAG为存储每行对角元位置的一维数组;COL为存储每个非零元列号;LINKN为列方向的行号纵向链;LINKP为列方向的位置纵向链;G、B为导纳值的一维数组。为了减少内存占用,状态估计过程中有些矩阵元素值可以不必存储,例如H的转置阵HT以及信息矩阵G。由于这些矩阵可以由H阵转置或者HTR-1H相乘获得,因此只需要在迭代前形成HT的索引并对HTH阵进行符号定位,等到真正需要这些元素参与计算时,再计算出它们的值。2.2快速前代和快速回代运算对于式(3)所示的方程式,记右端项为列向量b=HTR-1Δz,通过因子分解可将信息矩阵G表示为其中:U阵为对角元为1的上三角矩阵;D为n维对角阵。定义Uk为U阵中的保留第k行元素而其他对角元为1的n阶因子矩阵,即式中:mk是行向量,为U阵中除去对角元后的第k行;Ik是列向量,为单位阵的第k列。不难验证Uk-1=I-Ikmk,那么U阵及其逆阵可以表示成一序列因子矩阵的乘积,如式(6)因此,求解方程(3)可分为3步:其中:式(7)称为前代运算;式(8)称为规格化;式(9)称为回代运算。当右端项向量b中仅有少量非零元时,如果求解式(7)只在非零元的消去树路径决定的因子矩阵上进行时,可以避免零元参与运算,此时称为快速前代运算;同样,在求解式(9)的过程中如果我们仅对Δx中的少量分量感兴趣,就只需在这几个分量的消去树路径决定的因子矩阵上运算,这样可以避免不必要的运算操作,此时称为快速回代运算。例如,对于图1所示由信息矩阵G确定的网络,其消去树如图2所示(图1中虚线表示因子分解过程中增加的注入元)。图2中,边的方向是由小号指向大号,节点k代表U阵中k行确定的因子矩阵Uk,节点9是消去树的根节点,代表单位矩阵。除根节点外,每个节点都有一个父节点,定义为对应因子矩阵第一个非零非对角元的列号。每个节点的消去树路径是指包含该节点直到树根上的所有节点的集合,如节点5的消去树路径指的是{5,7,8,9}。如果右端项中只有非零元b1,根据式(7)在计算y向量时,我们只需在b1路径确定的因子矩阵U1、U4、U7、U8、U9上进行快速前代运算,即同样如果我们只对Δx向量中的Δx5感兴趣,根据式(9),我们只需在U9、U8、U7、U5上进行快速回代运算,即实际计算过程中,稀疏向量b中可能有多个非零元,或者对Δx向量中几个分量感兴趣,这时参与快速前代或快速回代运算的因子矩阵是这几个元素消去树路径的合集。如对b1和b5确定的消去树路径是U1、U4、U5、U7、U8、U9,其中节点7称为分叉点。快速前代运算时必须保证集合中比分叉点7小的因子矩阵按顺序运算完成后才能使U7参与运算;同样,快速回代运算时必须保证比分叉点大的因子矩阵按顺序完成后方可使U7参与运算。由此可见,提高稀疏向量法效率的措施可以从2方面着手:一方面尽量减少U阵中元素,使参与计算的元素尽可能的少,即减少因子分解过程中产生的注入元;另一方面可以尽可能减少参与计算的因子矩阵个数,也就是消去树路径的长度来减少总的参与计算的元素。计算前的节点优化编号正是从这2方面来提高稀疏向量法的效率。3稀疏技术在状态评估中的应用3.1方程的求解3.1.1基于ldlt因子表的求解策略采用直接法求解法方程式(3)时,可以采用LDLT因子分解法,也可以采用正交分解法或者其他方法。正交分解法数值稳定性好,但计算量很大,并且编程复杂,在需要多次分解G阵时计算效率较低;LDLT分解法当量测权重相差悬殊时数值稳定性较差,但其计算量比正交分解法小得多,是电力系统计算中对称矩阵分解常用的算法。本节研究了LDLT因子分解技术,包括存储策略分析、G阵和U阵稀疏结构的形成以及采用稀疏向量法求解LDLT因子表等。LDLT因子分解的计算公式如下其中:gii、gij为G阵中元素;uki、ukj是U阵的元素;dii是D阵对角元。以一次形成U阵的一个行向量为例,如第i行,由(10)可见,计算dii、uij时需要查询定位U阵的第i列元素以及前i-1行元素。如果U阵存储只保留每行非零非对角元起始位置、非零元列号及值,则查询定位U阵的i列元素具有O(n2)阶次的运算操作。对于大系统而言,这种存储策略查询定位效率很低。因此,本文存储U阵时,一并存储了U阵元素的位置纵向链和行号(也简称列链)。这样与U阵第i行元素有关信息如下:在D阵第i个对角元处不仅记录了U阵第i行第一个非零非对角元位置,而且记录了U阵中第i列行号最小非零元的位置,而U阵每个非零非对角元则记录了自己所在行号、列号以及同列下一个非零元位置;如果到了同列元素中最大行号元素处,则记为-1,表示链表末端。这样的存储策略可以保证从D阵任何对角元出发,从行和列2个方向顺序存取U阵元素。D阵和U阵示意图如图3,图中U方向第一个箭头表示D阵元素存储了指向U阵第一个元素的位置,其他箭头表示U阵元素的顺序存放关系。3.1.2生成g阵上三角阵的实际操作由式(3)可知,G阵稀疏结构可以采用模拟矩阵HTH相乘的方法获得,虽然因为G阵对称而只需作上三角阵的稀疏结构分析,但效率还是较低。文献对于每个状态量i,扫描与i有关的测点的其他状态量j(j>i),然后采用插入排号法形成G阵上三角阵的一行。这种算法在形成一行的过程中由于需要多次插入排号,影响了效率。文献介绍G阵结构和H阵间的关系,据此可以形成由量测决定的节点邻接矩阵。结合节点优化编号后新旧节点间的对应关系,这个邻接矩阵正好反映的是G阵的稀疏结构。因此可以直接在节点i的邻接节点中取出比i大的节点号,然后快速排序一次形成G阵一行的稀疏结构。待整个G阵行的非零分布确定后,扫描一遍G阵元素位置,形成G阵列链。3.1.3消去树对g阵的非零分布的认识U阵稀疏结构可以由一次模拟的LDLT分解形成,但其算法复杂度等同于数值分解,效率较低。文献介绍了基于消去树形成U阵稀疏结构的算法。令Struct(Gi)表示G阵第i行的非零分布,Struct(Ui)表示U阵第i行的非零分布,P(j)表示j节点在消去树中父节点,则U阵i行的非零分布可用下式表示:根据式(11)可形成U阵每行的非零分布信息,然后形成U阵的列链。由式(11)可见,U阵i行非零分布的形成依赖于G阵的消去树,因此必须首先形成G阵消去树。消去树的形成算法如下:算法中adj(xi)表示与节点xi相邻的节点集合。3.1.4在dlt的分解法本文提出一种应用稀疏向量法进行LDLT分解的算法。假设分解进行到了第k步,设Uk-1表示U阵的第k-1阶子式,U(1:k-1,k)表示U阵第k列不包括对角元的列向量,简记为v,相应的对于G阵有Gk-1子阵、G(1:k-1,k)(简记为w)以及gkk,对于D阵有Dk-1和dkk,此时,如下等式成立:即所以有如下等式(12)成立因此LDLT分解的算法可以描述如下:第一步:由UkT-1y=w采用稀疏向量法求y;第二步:v=D-1k-1y;第三步: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阵元素初始化工作行,由于在存储结构上没有形成U阵的列链和初始化的工作行没有为注入元预留位置,所以不得不在U阵列方向上搜索以及在U阵行方向上向工作行插入注入元,从而导致算法效率较低。算法2:仍然采用文献中按行LDLT分解方法,但根据3.1.2节和3.1.3节所述形成了U阵和G阵的列链。由于列链的作用,可以避免在列方向上搜索U阵元素;同时由于U阵结构已知,因而可以根据U阵结构分配工作行和根据G阵值初始化工作行,这样注入元的位置可以在工作行上预留出来。算法3:基于3.1.2和3.1.3节形成的G阵和U阵结构,采用稀疏向量法求解。由于每次稀疏向量求解元素的位置和顺序都是确定的,无需形成工作行,因而计算效率得以提高。表1是在IEEE118和IEEE300系统上测试上述各种算法效率的计算时间结果。两系统量测分别取1098个和2544个,节点优化编号方法采用最小度最少前趋节点数法(MD-MNP)。表中时间为计算10000次的平均时间,单位ms。测试结果表明,系统规模越大,相对而言算法3更具有优势。3.2确定残差方差使用标准化残差搜索法或者递归量测误差估计辨识法时,需要计算量测的方差。据文献,残差∆z的协方差矩阵为式(13)WR阵的对角元即对应量测残差方差为求解式(14)可以采用不同的稀疏技术。(1)方法1。将式(14)写成分量形式:由于hk是稀疏向量,因此可以采用快速前代技术先求出向量U-Thk,然后作向量的乘积即可求得对应的残差方差。(2)方法2。利用状态估计收敛后已有的因子表计算信息矩阵G的逆阵,然后做矩阵与向量的相乘,可以求得残差方差dk。据式(4)进而可以写成根据式(5)(6),Uk-1可以直接由Uk获得,所以G-1的求取只需要一序列因子矩阵乘积即可。并且由于G阵是对角占优的对称正定阵,可以利用G阵的“稀疏逆”近似替代G-1阵,即只计算G-1上三角阵中与U阵分布相同的非零元。令C=G-1,则C中元素可通过下