一种基于信息矩阵稀疏求解的电力系统状态估计方法和系统与流程

本发明涉及电力系统调度自动化技术领域，特别是一种对状态估计信息矩阵的求解过程进行优化的电力系统状态估计方法和系统。

背景技术：

目前调度控制系统中实际生产运行系统的状态估计程序大部分都是基于最小二乘原理的快速分解状态估计算法(fastdecoupledstateestimation，fdse)，其迭代计算公式如下式所示。

式中：z为量测矢量(m维)，为状态量(n维)，r为量测误差方差阵，h为非线性量测矢量函数，为雅可比矩阵，g＝h^tr^-1h为信息矩阵，z一般包括支路功率、节点注入功率和节点电压幅值等量测，为状态量修正量。

计算信息矩阵g＝h^tr^-1h可描述为稀疏矩阵相乘问题(sparsematrixmultiplication，spmm)；式(1-a)为线性方程组求解问题，数学上普遍基于三角分解法进行求解，由于g为对称实矩阵，一般将其分解为ldu，其中l＝ut,为单位下三角阵，d为对角阵。

随着调控云实时数据平台的建设，电力系统状态估计计算规模日益增大，对于国分和省地实时数据平台，状态估计计算节点数目达到1万～2万节点，量测数目达到5万～7万个，使得已有状态计算程序在效率上难以满足实时计算分析的需求。在fdse中信息矩阵的计算以及迭代公式线性方程组系数矩阵的lu因子分解占据了状态估计大部分时间，提高状态估计信息矩阵和lu分解的计算速度对整个电力系统状态的加速起着举足轻重的作用。

传统状态估计程序中对信息矩阵的求解虽然采用稀疏存储，但在矩阵乘法、因子分解过程中仍按稠密向量进行运算，从而导致计算效率的降低。在状态估计中，信息矩阵的计算可以用式c＝ab表示，其中，a、b是实数稀疏矩阵。设稀疏矩阵a、b的维数分别为m1*n1和m2*n2，且n1＝m2，a、b均以csr(compressedsparserow)格式进行存储，ai_nz，bi_nz分别为a、b第i行非零元数目。设c(i)为c第i行向量，c(i,j)为c(i)第j列非零元值，采用逐行乘积累积求和的方式求解c＝ab，有矩阵c一般也为稀疏矩阵。由于不能预先知道稀疏向量c(i)/矩阵c采用csr格式存储占用的空间大小，传统的实现方法一般采用稠密向量c(i)，维数定义为n2，在计算c(i,j)前对c(i)进行初始化或者预先虚拟计算得到c(i)需要的存储空间，但这样会消耗很长时间，影响计算效率。例如，采用稠密向量c(i)，由于预先未知c(i)的非零元分布，采用逐行乘积累积求和前需要将c(i)中所有元素初始化为零，c(i)运算结束后再采用稀疏压缩存储从而实现c的csr格式存储，这样操作存在的主要问题是稀疏矩阵c中原本不需要存储的零元在稠密向量c(i)中需要初始化为零。

在状态估计问题求解中，当计算节点数目n较大时，需要做大量的内存初始化(或者内存分配释放)工作，从而会降低计算效率，对3万节点规模的电网在计算信息矩阵时，采用单精度浮点float类型、非零索引采用4字节整型时，求解c(i)需要初始化的内存达到0.228mb，求解c需要初始化的内存达到68.66mb，当节点规模更大或者采用双精度时需要初始化的内存更大，会大大增加高速缓存与内存之间的数据传递量，导致高速缓存命中率较低，严重影响了信息矩阵的计算效率。

在稀疏矩阵因子分解过程中会产生非零注入元，在按行/列消去运算前无法确定因子表非零元的分布，常规的按行/列分解方法是将系数矩阵行稀疏向量展开为稠密向量，并基于稠密向量进行消去运算，每行/列运算结束后再采用稀疏压缩存储从而实现因子表的csr/csc格式存储。类似于spmm问题，这种实现方法存在的主要问题是基于稀疏向量展开为稠密向量时原本不需要存储的零元在稠密向量中进行了存储，当计算节点数目n较大时，即使程序能够利用空间局部性提高命中率，但大量的数据仍需要在缓存和主存之间进行传送，影响了三角分解的效率。

技术实现要素：

本发明的目的是，提供一种电力系统状态估计方法和系统，可减少求解过程中的内存操作，减小计算负荷，提升电网状态估计的计算效率。

本发明采取的技术方案如下：

一方面，本发明提出一种电力系统状态估计方法，包括：

实时获取电力系统一次设备遥测遥信数据；

基于获取到的数据构建电力系统节点支路模型；

根据节点支路模型确定有功雅可比矩阵的非零元分布；

根据有功雅可比矩阵的非零元分布，确定有功信息矩阵的非零元分布；

基于有功信息矩阵非零元分布结构进行节点优化编号；

根据新的节点顺序计算有功雅可比矩阵和无功雅可比矩阵，进而得到新的有功信息矩阵和无功信息矩阵；

对新的有功信息矩阵和无功信息矩阵分别进行因子分解，得到有功信息矩阵因子表和无功信息矩阵因子表；

基于有功信息矩阵因子表和无功信息矩阵因子表进行电力系统状态估计。

本发明通过基于有功信息矩阵非零元分布进行节点优化编号，而非基于节点出线度和节点量测数进行节点优化编号，可减少后续信息矩阵因子分解过程中非零元产生的数目；在此基础上，本发明对新的有功信息矩阵和无功信息矩阵分别进行因子分解，避免了将稀疏向量展开为稠密向量的工作，以及分解过程中需要频繁初始化数组内存并运算引起的效率低下，可有效提高高维度稀疏矩阵三角分解的计算效率，也即减少了求解过程中的内存操作，减小了计算负荷，提升电网状态估计的计算效率。

可选的，所述有功雅可比矩阵非零元分布和有功信息矩阵非零元分布，以及节点优化编号后的有功信息矩阵和无功信息矩阵，分别采用c++标准模板库关联容器定义结构进行稀疏存储；有功信息矩阵因子表和无功信息矩阵因子表存入采用c++标准模板库关联容器定义的数据结构，或者转换为按行/按列稀疏储存。基于stl关联容器实现稀疏矩阵的存储和运算，降低了spmm、节点优化编号和三角分解等三个过程中不必要的内存操作，可有效降低内存和高速缓存之间数据传递量，显著提高了状态估计的计算效率。

可选的，各电力系统一次设备遥测遥信数据，并行从调度自动化系统中分别针对各一次设备类型设置的实时库单表中获取。

可选的，定义m为有功量测数，n为电网计算节点数，则有功雅可比矩阵hp为m*(n-1)维稀疏矩阵，通过模拟乘法运算形成有功信息矩阵为：其中为有功量测误差方差阵；

稀疏矩阵hp、gp非零元分布采用c++stl提供的关联容器map和set定义结构进行存储，包括：定义一结构：typedefstd::set＜int＞sparserownzset用于存储某一行非零元列号集合，列号由set中int型(整型变量)的索引值表示；定义一结构typedefstd::map＜int,sparserownzset＞sparsematrixnzmap用于存储稀疏矩阵hp、gp非零元的分布，行号由int型的键索引表示。

在模拟乘法运算时，不进行实际的数值运算，只进行gp非零元结构的确定，对于hp、gp中零元素不进行存储。

由于hp中各元素直接完全相互独立，因此可通过并行方式形成有功雅可比矩阵hp非零元分布。

可选的，定义有功信息矩阵gp的非零元分布为g_nz_map，其存储于sparsematrixnzmap结构中，基于有功信息矩阵非零元分布结构进行节点优化编号包括：

基于g_nz_map获得初始系数矩阵各出线度的节点组合，定义为d_bus_map，存储于基于c++stl关联容器预先定义的typedefstd::map＜int,std::set＜int＞＞degreebusmap结构中；

采用模拟高斯消元法和半动态最小度算法对有功信息矩阵gp的节点进行排序，将排序过程中各出线度的节点集合分别存储至degreebusmap结构中；

degreebusmap结构中的各元素以出线度值作为键值，实际值为出线度等于键值的节点set集合。节点优化编号过程中只进行网络结构变化的处理，不进行实际的数值运算。

可选的，对有功信息矩阵gp的节点进行排序的过程中，使用c++stl提供的insert、erase操作实现对set集合中新增注入元的插入和已消去节点的删除。可避免自行设计复杂数据结构及其排序、增加、删除等算法操作，提高开发编程效率。

可选的，采用模拟高斯消元法和半动态最小度算法对有功信息矩阵gp的节点进行排序包括：

s1，确定当前g_nz_map和d_bus_map中的节点元素；

s2，从当前d_bus_map中选取节点出线度最小的节点作为当前次消去节点，记为k；

s3，根据g_nz_map修正与节点k及其关联节点的关联节点集合set，记为link_set；

s4，根据节点k的键值，从g_nz_map中删除节点k，并对当前删除的节点k进行编号；

s5，根据节点删除后的g_nz_map修正d_bus_map，作为更新后的当前d_bus_map；

s6，重复步骤s1至s6，直至初始d_bus_map中的所有节点排序完成，结束排序操作。

半动态最小度排序算法的核心思想是在分解的过程中将非零元最少的节点优先编号，在消去节点时修改剩余节点所在行的非零元数目，已编号的节点及其出发的边不再参与后面的模拟消去运算。

可选的，节点优化编号后，根据新的节点顺序并行计算有功雅可比矩阵hp和无功雅可比矩阵hq，则有功信息矩阵和无功信息矩阵其中，为无功测误差方差阵，gq为n×n维对称稀疏矩阵；

将hp、hq、gp和gq采用c++stl提供的关联容器定义结构进行稀疏矩阵的存储，矩阵中的零元素不进行存储。

可选的，对于节点优化编号后采用sparserowmap结构存储的gp和gq，并行采用高斯消元法进行因子分解，新产生的非零注入元调用关联容器的insert操作，或者将map作为关联式数组进行非零元的插入。在基于sparserowmap进行消去运算后，因子表非零元的列号索引将自动进行排序。

本发明还提供一种电力系统状态估计系统，包括：

数据获取模块，用于实时获取电力系统一次设备遥测遥信数据；

节点支路模型构建模块，用于基于获取到的数据构建电力系统节点支路模型；

有功雅克比矩阵确定模块，用于根据节点支路模型确定有功雅可比矩阵的非零元分布；

有功信息矩阵确定模块，根据有功雅可比矩阵的非零元分布，确定有功信息矩阵的非零元分布；

节点排序模块，用于基于有功信息矩阵非零元分布结构进行节点优化编号；

信息矩阵确定模块，用于根据新的节点顺序计算有功雅可比矩阵和无功雅可比矩阵，进而得到新的有功信息矩阵和无功信息矩阵；

因子分解模块，用于对新的有功信息矩阵和无功信息矩阵分别进行因子分解，得到有功信息矩阵因子表和无功信息矩阵因子表；

以及后续处理模块，用于基于有功信息矩阵因子表和无功信息矩阵因子表进行电力系统状态估计。

有益效果

与现有技术相比，本发明具有以下进步：

(1)通过基于有功信息矩阵非零元分布进行节点优化编号，并对节点排序后的有功信息矩阵和无功信息矩阵采用高斯消元法进行因子分解，可减少信息矩阵因子分解过程中非零元产生的数目，避免了将稀疏向量展开为稠密向量的工作，以及分解过程中需要频繁初始化数组内存并运算引起的效率低下，提高高维度稀疏矩阵三角分解的计算效率；

(2)通过基于stl关联容器实现稀疏矩阵的存储和运算，降低了spmm、节点优化编号和三角分解等三个过程中不必要的内存操作，可有效降低内存和高速缓存之间数据传递量，显著提高了状态估计的计算效率；

(3)基于stl关联容器编程可有效利用stl先进高效的算法实现内存数据的管理，避免了自行编写复杂的数据结构及算法，降低了编程的复杂度，有效提升了编程效率和程序品质；

(4)本发明方法能够有效兼容目前调度控制系统中已有的软硬件平台，不改变状态估计应用的部署架构和总体流程，可以复用已有的程序功能。

附图说明

图1所示为本发明状态估计方法中信息矩阵求解流程示意图；

图2所示为雅可比矩阵非零元分布示意图；

图3所示为信息矩阵求解示意图；

图4所示为基于stl关联容器的快速节点排序流程示意图。

具体实施方式

以下结合附图和具体实施例进一步描述。

传统状态估计程序中对信息矩阵的求解虽然采用稀疏存储，但在矩阵乘法、因子分解过程中仍按稠密向量进行运算，从而导致计算效率的降低。而本发明的发明构思为：通过在对信息矩阵的求解和矩阵乘法、因子分解过程中皆采用稀疏算法，提高状态估计的计算效率。

实施例1

本实施例为一种电力系统状态估计方法，参考图1所示，方法包括：

实时获取电力系统一次设备遥测遥信数据；

基于获取到的数据构建电力系统节点支路模型；

根据节点支路模型确定有功雅可比矩阵的非零元分布；

根据有功雅可比矩阵的非零元分布，确定有功信息矩阵的非零元分布；

基于有功信息矩阵非零元分布结构进行节点优化编号；

根据新的节点顺序计算有功雅可比矩阵和无功雅可比矩阵，进而得到新的有功信息矩阵和无功信息矩阵；

对新的有功信息矩阵和无功信息矩阵分别进行因子分解，得到有功信息矩阵因子表和无功信息矩阵因子表；

基于有功信息矩阵因子表和无功信息矩阵因子表进行电力系统状态估计。

通过基于有功信息矩阵非零元分布进行节点优化编号，而非基于节点出线度和节点量测数进行节点优化编号，可减少后续信息矩阵因子分解过程中非零元产生的数目；在此基础上，本发明对新的有功信息矩阵和无功信息矩阵分别进行因子分解，可采用高斯消元法，能够避免将稀疏向量展开为稠密向量的工作，以及分解过程中需要频繁初始化数组内存并运算引起的效率低下，可有效提高高维度稀疏矩阵三角分解的计算效率，也即减少了求解过程中的内存操作，减小了计算负荷，提升电网状态估计的计算效率。

在上述基础上，本实施例中，对于有功雅可比矩阵非零元分布和有功信息矩阵非零元分布，以及节点优化编号后的有功信息矩阵和无功信息矩阵，分别采用c++标准模板库关联容器定义结构进行稀疏存储；有功信息矩阵因子表和无功信息矩阵因子表存入采用c++标准模板库关联容器定义的数据结构，或者转换为按行/按列稀疏储存。基于stl关联容器实现稀疏矩阵的存储和运算，降低了spmm、节点优化编号和三角分解等三个过程中不必要的内存操作，可有效降低内存和高速缓存之间数据传递量，显著提高了状态估计的计算效率。

实施例1-1

基于实施例1的基本流程，本实施例对电力系统状态估计方法的流程进行具体描述。

一、电力系统一次设备遥信遥测数据获取

各电力系统一次设备遥测遥信数据，并行从调度自动化系统中分别针对各一次设备类型设置的实时库单表中读取。电力系统一次设备包括发电机、变压器、线路、负荷、容抗器、断路器、隔离开关等。

二、节点支路模型构建

根据调度自动化系统中的关系型电网模型，包含一次设备的连接关系、设备参数等数据，进行拓扑分析，进而形成节点支路模型，具体可参考现有技术。

三、雅可比矩阵形成

根据节点支路模型并行形成有功雅可比矩阵hp非零元分布，如图2所示，通过模拟乘法计算形成有功信息矩阵非零元分布，如图3所示。

有功雅可比矩阵hp为m*(n-1)维稀疏矩阵，m为有功量测数，n为电网计算节点数，有功信息矩阵中，为有功量测误差方差阵；

稀疏矩阵hp、gp非零元分布采用c++stl提供的关联容器map和set定义结构进行存储，包括：定义一结构：typedefstd::set＜int＞sparserownzset用于存储某一行非零元列号集合，列号由set中int型的索引值表示；定义一结构：typedefstd::map＜int,sparserownzset＞sparsematrixnzmap用于存储稀疏矩阵hp、gp非零元的分布，行号由int型的键索引表示；

在模拟乘法运算时，不进行实际的数值运算，只进行gp非零元结构的确定，对于hp、gp中零元素不进行存储。

由于hp中各元素直接完全相互独立，因此可通过并行方式形成有功雅可比矩阵hp非零元分布。

三、快速节点最小度节点优化编号

在快速分解状态估计中，传统方法是基于节点出线度(基于节点导纳矩阵y)、节点量测数(雅可比矩阵h)进行节点优化编号，这种编号方法虽然简单，但y、h非零元分布和fdse迭代方程中线性方程组系数g的非零元分布并不一致。为了减少信息矩阵因子分解过程中非零元产生的数目，本发明基于有功信息矩阵非零元结构进行节点排序。

在排序过程中，基于c++stl关联容器定义结构typedefstd::map＜int,std::set＜int＞＞degreebusmap用于记录排序过程中各出线度的节点集合，键值为出线度值，实际值为出线度等于键值的节点set集合。半动态最小度排序算法的核心思想是在分解的过程中将非零元最少的节点优先编号，在消去节点时修改剩余节点所在行的非零元数目，已编号的节点及其出发的边不再参与后面的模拟消去运算。在节点排序过程中使用c++标准模板库提供的insert、erase操作实现对set集合中新增注入元的插入和已消去节点的删除，避免了自行设计复杂数据结构及其排序、增加、删除等算法操作。在节点优化编号过程中只进行网络结构变化的处理，不进行实际的数值运算。

具体的，节点优化排序流程参考图4所示，定义有功信息矩阵gp的非零元分布为g_nz_map，其存储于sparsematrixnzmap结构中，基于有功信息矩阵非零元分布结构进行节点优化编号包括：

基于g_nz_map获得初始系数矩阵各出线度的节点组合，定义为d_bus_map，存储于基于c++stl关联容器预先定义的typedefstd::map＜int,std::set＜int＞＞degreebusmap结构中；

采用模拟高斯消元法和半动态最小度算法对有功信息矩阵gp的节点进行排序，将排序过程中各出线度的节点集合分别存储至degreebusmap结构中；

degreebusmap结构中的各元素以出线度值作为键值，实际值为出线度等于键值的节点set集合。

其中，采用模拟高斯消元法和半动态最小度算法对有功信息矩阵gp的节点进行排序包括：

s1，确定当前g_nz_map和d_bus_map中的节点元素；

s2，从当前d_bus_map中选取节点出线度最小的节点作为当前次消去节点，记为k；

s3，根据g_nz_map修正与节点k及其关联节点的关联节点集合set，记为link_set；

s4，根据节点k的键值，从g_nz_map中删除节点k，并对当前删除的节点k进行编号；

s5，根据节点删除后的g_nz_map修正d_bus_map，作为更新后的当前d_bus_map；

s6，重复步骤s1至s6，直至初始d_bus_map中的所有节点排序完成，结束排序操作。

四、有功和无功雅可比矩阵以及信息矩阵形成

节点优化编号后，根据新的节点顺序并行计算有功雅可比矩阵hp和无功雅可比矩阵hq，则有功信息矩阵和无功信息矩阵其中，为无功测误差方差阵，gq为n×n维对称稀疏矩阵；

将hp、hq、gp和gq采用c++stl提供的关联容器定义结构进行稀疏矩阵的存储，矩阵中的零元素不进行存储。gp和gq采用c++stl提供的关联容器map和set定义结构进行稀疏存储。定义结构：typedefstd::map＜int,float＞sparserowmap用于存储稀疏矩阵行非零元，列号由int型的键索引表示，稀疏矩阵非零元由float型的关联值表示；定义结构typedefstd::map＜int,sparserowmap＞sparsematrixmap用于存储稀疏矩阵，行号由int型的键索引表示。

在快速解耦状态估计中有功雅可比矩阵hp和无功雅可比矩阵hq相互独立，采用并行计算方法实现hp和hq的并行计算。

五、信息矩阵因子分解

由于有功信息矩阵gp和gq的因子表相互独立，因此对于节点优化编号后采用sparserowmap结构存储的gp和gq，可并行采用高斯消元法进行因子分解，得到相应的因子表。

由于消去运算前无法确定因子表非零元的分布，传统的方法是按行将系数矩阵行稀疏向量展开为稠密向量，并基于稠密向量进行消去运算，每行运算结束后再采用稀疏压缩存储从而实现因子表的稀疏格式存储。

本实施例对采用sparserowmap结构的信息矩阵gp和gq采用高斯消元法按行进行直接分解，新产生的非零注入元调用关联容器的insert操作或者将map作为关联式数组实现非零元的高效插入，并且基于sparserowmap进行消去运算后，因子表非零元的列号索引自动进行了排序。这样避免了将稀疏向量展开为稠密向量的工作以及按行分解需要频繁初始化数组内存及其运算引起的效率低下，有效提高了高维度稀疏矩阵三角分解的计算效率。

六、信息矩阵因子表的存储

本实施例将有功信息矩阵因子表和无功信息矩阵因子表存入采用c++stl关联容器定义的数据结构或者转换为按行/按列稀疏储存以备后续程序使用。

七、后续状态估计计算

后续状态估计计算基于步骤六存储的信息矩阵因子表进行，包括迭代计算节点电压等，可采用现有技术，不作为本发明的研究内容，不予赘述。

实施例2

与实施例1基于相同的发明构思，本实施例为一种电力系统状态估计系统，包括：

数据获取模块，用于实时获取电力系统一次设备遥测遥信数据；

节点支路模型构建模块，用于基于获取到的数据构建电力系统节点支路模型；

有功雅克比矩阵确定模块，用于根据节点支路模型确定有功雅可比矩阵的非零元分布；

有功信息矩阵确定模块，根据有功雅可比矩阵的非零元分布，确定有功信息矩阵的非零元分布；

节点排序模块，用于基于有功信息矩阵非零元分布结构进行节点优化编号；

信息矩阵确定模块，用于根据新的节点顺序计算有功雅可比矩阵和无功雅可比矩阵，进而得到新的有功信息矩阵和无功信息矩阵；

因子分解模块，用于对新的有功信息矩阵和无功信息矩阵分别采用高斯消元法进行因子分解，并存储分解得到的有功信息矩阵因子表和无功信息矩阵因子表；

以及后续处理模块，用于基于有功信息矩阵因子表和无功信息矩阵因子表进行电力系统状态估计。

上述各模块的具体功能实现参考实施例1-1中的具体实施方式。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的系统。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令系统的制造品，该指令系统实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。