CN113553538A - 一种递推修正混合线性状态估计方法 - Google Patents

一种递推修正混合线性状态估计方法 Download PDF

Info

Publication number
CN113553538A
CN113553538A CN202110532099.XA CN202110532099A CN113553538A CN 113553538 A CN113553538 A CN 113553538A CN 202110532099 A CN202110532099 A CN 202110532099A CN 113553538 A CN113553538 A CN 113553538A
Authority
CN
China
Prior art keywords
measurement
scada
system state
measurement data
pmu
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110532099.XA
Other languages
English (en)
Other versions
CN113553538B (zh
Inventor
孙康
卫志农
黄蔓云
孙国强
臧海祥
陈�胜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202110532099.XA priority Critical patent/CN113553538B/zh
Publication of CN113553538A publication Critical patent/CN113553538A/zh
Application granted granted Critical
Publication of CN113553538B publication Critical patent/CN113553538B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/06—Energy or water supply
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J3/00—Circuit arrangements for ac mains or ac distribution networks
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J2203/00—Indexing scheme relating to details of circuit arrangements for AC mains or AC distribution networks
    • H02J2203/20—Simulating, e g planning, reliability check, modelling or computer assisted design [CAD]
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E40/00—Technologies for an efficient electrical power generation, transmission or distribution
    • Y02E40/70—Smart grids as climate change mitigation technology in the energy generation sector
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
    • Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
    • Y04S10/00—Systems supporting electrical power generation, transmission or distribution
    • Y04S10/22—Flexible AC transmission systems [FACTS] or power factor or reactive power compensating or correcting units
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
    • Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
    • Y04S10/00—Systems supporting electrical power generation, transmission or distribution
    • Y04S10/50—Systems or methods supporting the power network operation or management, involving a certain degree of interaction with the load-side end user applications

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Human Resources & Organizations (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Economics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Strategic Management (AREA)
  • Development Economics (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Educational Administration (AREA)
  • Algebra (AREA)
  • Health & Medical Sciences (AREA)
  • Entrepreneurship & Innovation (AREA)
  • General Business, Economics & Management (AREA)
  • Tourism & Hospitality (AREA)
  • General Engineering & Computer Science (AREA)
  • Quality & Reliability (AREA)
  • Power Engineering (AREA)
  • Game Theory and Decision Science (AREA)
  • Computing Systems (AREA)
  • Public Health (AREA)
  • Water Supply & Treatment (AREA)
  • General Health & Medical Sciences (AREA)
  • Primary Health Care (AREA)
  • Measurement Of Resistance Or Impedance (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种递推修正混合线性状态估计方法，包括：获得电力系统的拓扑结构、网络参数以及量测数据；根据SCADA量测数据计算相对于参考节点的SCADA复数伪量测，估计粗略系统状态；根据PMU量测数据，对粗略系统状态进行递推修正，获取精确系统状态；根据获取的精确系统状态，进行权重递推校正；利用修正区域调整算法进行修正区域调整；建立基于流计算的多线程处理框架，将粗略系统状态和精确系统状态的估计过程拆分为若干个并行执行的进程，获取最终的精确系统状态输出。本发明提高了状态估计的估计精度和计算效率，提高了状态量更新的频率和状态估计的实时性，并且可以并行执行，保证了算法的可行性，使得估计性能更加稳定。

Description

一种递推修正混合线性状态估计方法
技术领域
本发明涉及一种递推修正混合线性状态估计方法，属于电力系统运行与调度的技术领域。
背景技术
状态估计是电力系统发展的基石，它通过对量测数据的滤波计算，确定当前时间断面下电力系统运行状态的最优估计值。近年来，随着电力系统规模的不断扩大和可再生能源的迅速渗透，电力系统的运行方式愈加复杂，状态变化也更为频繁。为了应对这些变化，电力系统状态估计算法也必须有更高的估计精度、计算效率和执行频率。确保状态数据的准确性和实时性是保证电力系统安全经济运行的首要条件。
状态估计器采用的原始量测数据一般来自于数据采集与监视控制系统(SCADA)或相量量测装置(PMU)。SCADA量测覆盖率高，但是其采样率低，无法保证时效性，且SCADA量测和系统状态变量之间是非线性的关系，因此往往需要利用牛顿法迭代求解。相比较之下，PMU采样率较高，且和系统状态变量之间是线性关系，无需迭代求解。但是PMU在系统中覆盖率低，往往只有关键节点安装了PMU。实际应用中，为了保证全系统的可观性，往往只采用SCADA量测作为原始数据，从而此类状态估计器相应的最多只能以SCADA量测相同频率提供估计值，在新能源大量接入导致系统状态急剧变化的今天，传统的状态估计算法已经无法满足实时跟踪系统状态的需要，状态估计作为能量管理系统中的关键应用，不及时更新的状态估计值势必也会对后续的高级应用造成影响，使得监测和控制能力无法满足现代电力系统的需要。
传统的电力系统状态估计方法有加权最小二乘法(Weighted Least Square,WLS)，加权最小绝对值法(Weight Least Absolut Value,WLAV)等。这些方法已在电力系统中成熟运用多年，但随着电网技术的不断发展，这些方法的准确度和实时性已经得不到很好的保证。半正定规划法(Semidefinite Programming)通过秩松弛重构了状态估计问题，从理论上改善了状态估计问题的求解复杂度，但是其依然是一个迭代算法，计算效率仍然偏低，其估计精度较传统的WLS等算法没有本质的提高。双线性状态估计通过将一步非线性迭代转化为两步线性迭代和两步非线性变换，提高了状态估计的效率，但是其两步线性迭代过程仍然是非常耗时的，且其速度的提高是以牺牲一定的估计精度来实现的。也有算法使用了SCADA和PMU混合量测，进行状态估计计算，提高了信息冗余度，增加了状态估计的精度，但是其计算效率反而下降了。基于数据驱动的状态估计也被提出，其利用历史量测和估计数据离线训练机器学习模型，并将模型在线应用于状态估计值的获取，无需模型计算，在线应用实时性好，然而其难以应对拓扑变化，状态估计精度也受突发实践在线模型更新频率等因素的影响。且由于SCADA量测采样率的限制，以上状态估计器的执行频率很低，无法很好的跟踪系统状态。有算法仅利用PMU的量测实现了线性状态估计，实现了更高的估计精度和计算效率，但是由于目前PMU装配率不足以实行全系统的状态估计，故仅使用PMU的线性状态估计算法仅能停留在理论层面，无法真正实施。
发明内容
为了解决当前电力系统状态估计结果更新周期较长，当前状态估计难以保证实时跟踪效果的问题，本发明提供一种递推修正混合线性状态估计方法，在使用混合SCADA/PMU量测的同时，将非线性的状态估计问题转化为非迭代的线性状态估计问题，实现了状态量的高精度、高效率求解；同时，设计了基于流数据的多线程运行模式，可以利用所有时间断面的SCADA和PMU数据，进一步提高了算法的计算效率和跟踪性能。
本发明具体采用以下技术方案解决上述技术问题：
一种递推修正混合线性状态估计方法，包括以下步骤：
步骤1：获得电力系统的拓扑结构、网络参数以及量测数据,其中量测数据包括SCADA量测数据和PMU量测数据；
步骤2：根据获取的SCADA量测数据计算相对于参考节点的SCADA复数伪量测；
步骤3：根据计算的SCADA复数伪量测，估计粗略系统状态；
步骤4：根据获取的PMU量测数据，利用系统状态递推修正算法对粗略系统状态进行递推修正，获取精确系统状态；
步骤5：根据获取的精确系统状态，进行权重递推校正；
步骤6：根据获取的精确系统状态，利用修正区域调整算法进行修正区域调整；
步骤7：建立基于流计算的多线程处理框架，将粗略系统状态和精确系统状态的估计过程拆分为若干个并行执行的进程，通过并行的进程，获取最终的精确系统状态输出。
进一步地，作为本发明的一种优选技术方案，所述步骤1中获得的网络参数包括输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、变压器变比和阻抗；以及，获得的量测数据中SCADA量测数据包括节点电压幅值、节点注入有功功率、节点注入无功功率、节点注入电流幅值、线路首端有功功率、线路首端无功功率、线路末端有功功率、线路末端无功功率、线路首端电流幅值以及线路末端电流幅值；以及，获得的量测数据中PMU量测数据包括节点复电压、线路首端复电流、线路末端复电流。
进一步地，作为本发明的一种优选技术方案，所述步骤2中计算相对于参考节点的SCADA复数伪量测，具体为：
1)建立SCADA复电压
Figure BDA0003066047150000031
量测方程：
Figure BDA0003066047150000032
其中，Uf,S是节点f的SCADA电压幅值量测值；
Figure BDA0003066047150000033
和
Figure BDA0003066047150000034
为节点f的复电压和相角指数值；
2)建立SCADA支路复电流
Figure BDA0003066047150000035
量测方程：
Figure BDA0003066047150000036
其中，
Figure BDA0003066047150000037
和
Figure BDA0003066047150000038
分别为线路f-t的并联导纳和支路导纳；
Figure BDA0003066047150000039
和
Figure BDA00030660471500000310
表示节点f和节点t的复电压，
Figure BDA00030660471500000311
为局部支路复电流伪量测；
3)建立SCADA节点复电流
Figure BDA00030660471500000312
量测方程：
Figure BDA00030660471500000313
其中，Lf表示表示和节点f通过线路直接相连的所有其他节点；
定义下标k,t表示第k次SCADA量测数据上传以来的第t次PMU量测数据上传对应的时刻，简称为(k,t)时刻或k周期的t时刻；zk,t和mk,t分别为SCADA复数伪量测和PMU量测数据，T和K分别为SCADA复数伪量测和PMU量测数据对应的量测方程系数矩阵，Wk,t和Vk,t分别为SCADA复数伪量测和PMU量测数据对应的量测权重矩阵，则相应的，SCADA复数伪量测的增益矩阵Gk,t和PMU量测数据的增益矩阵Fk,t分别定义为Gk,t＝TTWk,tT和Fk,t＝KTVk,tK；
根据以上SCADA复数伪量测方程，建立矩阵形式的量测方程：
Figure BDA0003066047150000041
其中，Lk,1是(k,1)时刻量测方程的系数矩阵，
Figure BDA0003066047150000042
是(k,1)时刻临时状态变量，εS是SCADA量测的误差向量，
Figure BDA0003066047150000043
是(k,1)时刻平衡节点的量测向量；
并且，若用
Figure BDA0003066047150000044
和δk,1表示(k,1)时刻的复电压和相角，则上述Lk,1和
Figure BDA0003066047150000045
可表示为：
Figure BDA0003066047150000046
Figure BDA0003066047150000047
其中，n表示节点数量，ES为常数阵，
Figure BDA0003066047150000048
为
Figure BDA0003066047150000049
和
Figure BDA00030660471500000410
组成的常数矩阵，
Figure BDA00030660471500000411
为导纳矩阵；使用
Figure BDA00030660471500000412
和
Figure BDA00030660471500000413
分别对应代表ES、
Figure BDA00030660471500000414
-Uf,S、
Figure BDA00030660471500000415
和
Figure BDA00030660471500000416
根据上述矩阵形式的量测方程建立线性最小二乘模型，求解可得临时状态变量为：
Figure BDA00030660471500000417
计算完成后，根据如下公式计算(k,1)时刻相对于参考节点的SCADA复数伪量测：
Figure BDA00030660471500000418
其中，
Figure BDA00030660471500000419
是Lk,1的后n-1列，
Figure BDA00030660471500000420
为相角指数量，即
Figure BDA00030660471500000421
的后n-1列。
进一步地，作为本发明的一种优选技术方案，所述步骤3中估计粗略系统状态，具体为：
利用下述线性WLS状态估计方程获取粗略系统状态：
Figure BDA00030660471500000422
其中，
Figure BDA00030660471500000423
为相对于参考节点的SCADA复数伪量测的量测矩阵，J为一个常数矩阵；
Figure BDA00030660471500000424
为粗略系统状态；Wk,t为SCADA复数伪量测的量测权重矩阵；zk,1为(k,1)时刻相对于参考节点的SCADA复数伪量测；
定义如下中间系数矩阵Ak,t＝Gk,t，Bk,t＝JTWk,t，则上述线性WLS状态估计方程表述为：
Figure BDA0003066047150000051
进一步地，作为本发明的一种优选技术方案，所述步骤4中利用系统状态递推修正算法对粗略系统状态进行递推修正，具体为：
当PMU量测数据mk,t到来后，通过如下的系统状态递推修正算法将粗略系统状态
Figure BDA0003066047150000052
修正为精确系统状态
Figure BDA0003066047150000053
Figure BDA0003066047150000054
其中，Wk,t和Vk,t分别为SCADA复数伪量测和PMU量测数据对应的量测权重矩阵，则SCADA复数伪量测的增益矩阵Gk,t和PMU量测的增益矩阵Fk,t分别定义为Gk,t＝JTWk,tJ和Fk,t＝KTVk,tK，公式中J和K分别为SCADA复数伪量测和PMU量测对应的量测方程系数矩阵；
定义如下中间系数矩阵Ck,t＝Gk,t+Fk,t，Dk,t＝KTVk,t，则上述系统状态递推修正算法表述为：
Figure BDA0003066047150000055
Figure BDA0003066047150000056
进一步地，作为本发明的一种优选技术方案，所述步骤5中进行权重递推校正，具体为：
获取k周期t时刻的SCADA复数伪量测和PMU的量测残差rk,t和ek,t如下：
Figure BDA0003066047150000057
Figure BDA0003066047150000058
其中，J和K分别为SCADA复数伪量测和PMU量测数据对应的量测方程系数矩阵，
Figure BDA0003066047150000059
为精确系统状态，zk,1为(k,1)时刻相对于参考节点的SCADA复数伪量测，mk,t为(k,1)时刻相对于参考节点的PMU量测数据；
在量测误差服从高斯分布时，在当前周期SCADA和PMU协方差
Figure BDA0003066047150000061
和
Figure BDA0003066047150000062
的基础上，通过如下的公式对rk,t和ek,t下一周期SCADA和PMU更新对应的协方差
Figure BDA0003066047150000063
和
Figure BDA0003066047150000064
进行递推估计：
Figure BDA0003066047150000065
Figure BDA0003066047150000066
其中，med{(·)}表示取中位数操作，Le是估计窗口长度，c＝1.483(1+5/(Le-1))是采样校正系数，α是修正遗忘因子；rk,t,rk-1,t,…,rk-Le+1,t分别表示k至k-Le+1周期下t时刻的SCADA复数伪量测的残差；ek,t,ek-1,t,…,ek-Le+1,t分别表示k至k-Le+1周期下t时刻的PMU量测数据的残差；
对量测权重Wk,t和Vk,t递推修正，得到量测权重的计算式为：
Figure BDA0003066047150000067
Figure BDA0003066047150000068
其中，i表示第i个元素，即第i个量测。
进一步地，作为本发明的一种优选技术方案，所述步骤6中利用修正区域调整算法进行修正区域调整，具体为：
当(k,1)时刻的SCADA复数伪量测和PMU量测数据到来时，通过如下的修正效果指标
Figure BDA0003066047150000069
来评估所有节点在k-1周期的修正效果：
Figure BDA00030660471500000610
其中，
Figure BDA00030660471500000611
表示(k,t)时刻的精确系统状态向量的第i个元素，即第i个节点，下标(·)k-1,τ表示k-1周期的最后一个时刻。此时，通过如下的调整指标
Figure BDA00030660471500000612
来选取k周期修正区域：
Figure BDA00030660471500000613
其中，Lm表示监控窗口的长度；
并且，将调整指标
Figure BDA0003066047150000071
的节点划分至修正区域内，只有对修正区域内的节点进行系统状态递推校正。
进一步地，作为本发明的一种优选技术方案，所述步骤7中将粗略系统状态和精确系统状态的估计过程拆分为若干个并行执行的进程，具体为：
进程1：当SCADA量测数据到来时，计算SCADA复数伪量测zk,t；
进程2：根据SCADA复数伪量测zk,t，计算粗略系统状态
Figure BDA0003066047150000072
若zk,t没有计算完成，则用k-1周期t时刻的SCADA复数伪量测zk-1,t代替；
进程3：当PMU量测数据到来时，修正粗略系统状态
Figure BDA0003066047150000073
为精确系统状态
Figure BDA0003066047150000074
若
Figure BDA0003066047150000075
没有计算完成，则用k-1周期t时刻的粗略系统状态
Figure BDA0003066047150000076
代替，以获取最终的精确系统状态输出；
进程4：当精确系统状态
Figure BDA0003066047150000077
计算完成时，进行k+1周期t时刻的SCADA复数伪量测的量测权重矩阵Wk+1,t和PMU量测数据的量测权重矩阵Vk+1,t的递推校正，基于校正后的量测权重矩阵，计算下一个周期的中间系数矩阵Ak+1,t，Bk+1,t，Ck+1,t，Dk+1,t，并对这四个矩阵进行LU分解用于下一个周期其他进程的计算；
进程5：当(k,1)时刻的SCADA复数伪量测和PMU量测数据到来时，进行修正区域调整计算，获取第k周期采用的修正区域。
本发明采用上述技术方案，能产生如下技术效果：
本发明的方法首先仅根据数据采集与监视控制系统(SCADA)的实际量测计算SCADA复数形式伪量测；接着，根据SCADA复数形式伪量测和相量量测装置(PMU)的实际量测对系统状态变量进行实时递推估计；最后，对量测权重递推校正，并对修正区域做适应性调整。此外，本发明所提算法的计算可在基于流计算的多线程处理框架下执行。IEEE标准系统的测试结果表明，由于流计算和递推修正理论的引入，本发明提出的方法估计精度和计算效率均高于已有的状态估计算法，而对PMU量测的有效利用很好的保证了状态估计器的实时跟踪能力，使得本发明提出的方法可以以PMU相同的刷新率提供状态估计结果。
因此，本发明方法提高了状态估计的估计精度和计算效率，提高了状态量更新的频率，提高了状态估计的实时性。所提的基于流计算的多线程处理框架，使得算法可以并行执行，保证了算法的可行性，使得估计性能更加稳定。故本发明在提高了状态估计精度的同时，大大提高了状态估计的计算效率，保证了状态估计的实时性。
附图说明
图1为本发明的递推修正混合线性状态估计方法的流程图。
图2为本发明中多线程并行处理示意图。
图3为现有技术中传统的传统流计算框架示意图。
图4为本发明中状态估计的流计算框架示意图。
图5为本发明中三节点系统算例图。
图6为本发明中负荷波动曲线示意图。
图7为本发明中不同状态估计算法下IEEE14节点状态估计的估计误差比较示意图。
图8为本发明中不同状态估计算法下IEEE14节点9状态估计结果比较示意图。
图9为本发明中不同PMU安装率下状态估计的估计误差比较示意图。
图10为本发明中不同PMU安装率下状态估计的计算时间比较示意图。
具体实施方式
下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
如图1所示，本发明涉及一种递推修正混合线性状态估计方法，该方法具体包括以下步骤：
步骤1：获得电力系统的拓扑结构、网络参数以及量测数据，其中量测数据包括SCADA量测数据和PMU量测数据，具体如下：
①获取拓扑结构；
②获取网络参数，包括输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、变压器变比和阻抗；
③获取量测数据，包括SCADA量测数据和PMU量测数据；其中SCADA量测数据包括：节点电压幅值、节点注入有功功率、节点注入无功功率、节点注入电流幅值、线路首端有功功率、线路首端无功功率、线路末端有功功率、线路末端无功功率、线路首端电流幅值以及线路末端电流幅值；其中PMU量测数据包括：节点复电压、线路首端复电流、线路末端复电流。
步骤2：根据获取的SCADA量测数据计算相对于参考节点的SCADA复数伪量测，具体如下：
定义局部支路复电流伪量测
Figure BDA0003066047150000091
如下：
Figure BDA0003066047150000092
其中，上标local表示该量测为局部量测，即其角度是相对于节点f而非参考节点，下标f和t分别代表节点f和节点t，下标ft表示连接节点f和节点t的线路，下标S代表SCADA量测数据。Ift,S表示支路f-t的SCADA支路电流量测。P和Q分别为有功量测值和无功量测值，并有线路f-t的功率因数角θft＝tg-1(-Qft/Pft)和节点f的功率因数角θf＝tg-1(-Qf/Pf)。此时相对于参考节点的复电流可以表示为：
Figure BDA0003066047150000093
其中，δf是节点f相对于参考节点的相角，它的值是未知的，可通过后续的算法获得。
建立以复电压和相角指数值为状态变量的SCADA量测方程如下：
1)建立SCADA复电压
Figure BDA0003066047150000094
量测方程：
Figure BDA0003066047150000095
其中，Uf,S是节点f的SCADA电压幅值量测值。
Figure BDA0003066047150000096
和
Figure BDA0003066047150000097
为节点f的复电压和相角指数值，共同构成量测方程的状态变量。
2)建立SCADA支路复电流
Figure BDA0003066047150000098
量测方程：
Figure BDA0003066047150000099
其中，
Figure BDA00030660471500000910
和
Figure BDA00030660471500000911
分别为线路f-t的并联导纳和支路导纳。
Figure BDA00030660471500000912
和
Figure BDA00030660471500000913
表示节点f和节点t的复电压，
Figure BDA00030660471500000914
为相角指数值，共同构成量测方程的状态变量。
3)建立SCADA节点复电流
Figure BDA00030660471500000915
量测方程：
Figure BDA00030660471500000916
其中，Lf表示表示和节点f通过线路直接相连的所有其他节点。
定义下标k,t表示第k次SCADA数据上传以来的第t次PMU数据上传对应的时刻，可简称为(k,t)时刻或k周期的t时刻，zk,t和mk,t分别为SCADA复数伪量测和PMU量测，J和K分别为SCADA复数伪量测和PMU量测对应的量测方程系数矩阵，Wk,t和Vk,t分别为SCADA复数伪量测和PMU量测对应的量测权重矩阵，则相应的，SCADA复数伪量测的增益矩阵Gk,t和PMU量测的增益矩阵Fk,t分别定义为Gk,t＝JTWk,tJ和Fk,t＝KTVk,tK。
根据以上SCADA复数伪量测方程，建立矩阵形式的量测方程：
Figure BDA0003066047150000101
其中，Lk,1是(k,1)时刻量测方程的系数矩阵，
Figure BDA0003066047150000102
是(k,t)时刻临时状态变量，εS是SCADA量测的误差向量，
Figure BDA0003066047150000103
是(k,1)时刻平衡节点的量测向量。
并且，若用
Figure BDA0003066047150000104
和δk,1表示(k,1)时刻的复电压和相角，上述Lk,1和
Figure BDA0003066047150000105
可以表示为：
Figure BDA0003066047150000106
Figure BDA0003066047150000107
其中，n表示节点数量，ES为常数阵，有相对应的SCADA电压幅值量测的位置为1，其余位置为0，若所有节点均有SCADA电压幅值量测，则ES为单位矩阵，
Figure BDA0003066047150000108
为
Figure BDA0003066047150000109
和
Figure BDA00030660471500001010
组成的常数矩阵，
Figure BDA00030660471500001011
为导纳矩阵。为方便表述，使用
Figure BDA00030660471500001012
和
Figure BDA00030660471500001013
分别对应代表ES、
Figure BDA00030660471500001014
-Uf,S、
Figure BDA00030660471500001015
和
Figure BDA00030660471500001016
为解释上述矩阵形式的量测方程，以图5所示的三节点系统为例，图中1、2、3分别对应于节点1、2、3，则有：
Figure BDA0003066047150000111
Figure BDA0003066047150000112
Figure BDA0003066047150000113
其中，(k,1)时刻的状态变量
Figure BDA0003066047150000114
中有
Figure BDA0003066047150000115
而
Figure BDA0003066047150000116
Figure BDA0003066047150000117
分别是节点1，2，3的电压幅值状态变量，
Figure BDA0003066047150000118
和
Figure BDA0003066047150000119
分别是节点2和3的相角指数值状态变量。
根据上述量测方程建立线性最小二乘模型，求解可得临时状态变量为：
Figure BDA00030660471500001110
计算完成后，根据如下公式计算(k,1)时刻相对于参考节点的SCADA复数伪量测：
Figure BDA00030660471500001111
其中，
Figure BDA00030660471500001112
是Lk,1的后n-1列，
Figure BDA00030660471500001113
为相角指数量，即
Figure BDA00030660471500001114
的后n-1列。
步骤3：根据计算的SCADA复数伪量测，估计粗略系统状态，具体如下：
利用下述线性WLS状态估计方程获取粗略系统状态：
Figure BDA00030660471500001115
其中，
Figure BDA00030660471500001116
为相对于参考节点的SCADA复数伪量测的量测矩阵，
Figure BDA00030660471500001117
为粗略系统状态。进一步定义如下中间系数矩阵Ak,t＝Gk,t，Bk,t＝JTWk,t，则上述线性WLS状态估计方程可表述为：
Figure BDA00030660471500001118
其中，J是一个常数矩阵，当网络拓扑的参数不变时，不需要对其重新计算，Wk,t可通过后续算法获得。
步骤4：根据获取的PMU量测数据，利用系统状态递推修正算法对粗略系统状态进行递推修正，获取精确系统状态，具体如下：
当新的PMU量测数据mk,t到来后，可以通过如下的系统状态递推修正算法将粗略系统状态
Figure BDA0003066047150000121
修正为精确系统状态
Figure BDA0003066047150000122
Figure BDA0003066047150000123
进一步定义中间系数矩阵Ck,t＝Gk,t+Fk,t，Dk,t＝KTVk,t，则上述系统状态递推修正算法可表述为：
Figure BDA0003066047150000124
以及，中间系数矩阵：
Figure BDA0003066047150000125
其中，对公式
Figure BDA0003066047150000126
的推导过程如下：
对于SCADA和PMU混合量测的线性状态估计问题，其(k,t)时刻的量测矩阵
Figure BDA0003066047150000127
权重矩阵
Figure BDA0003066047150000128
和量测向量ζk,t为：
Figure BDA0003066047150000129
Figure BDA00030660471500001210
Figure BDA00030660471500001211
此时，状态估计结果
Figure BDA00030660471500001212
和相应的增益矩阵
Figure BDA00030660471500001213
为：
Figure BDA00030660471500001214
Figure BDA0003066047150000131
根据
Figure BDA0003066047150000132
则有：
Figure BDA0003066047150000133
Figure BDA0003066047150000134
此时，则有：
Figure BDA0003066047150000135
将此式带入式
Figure BDA0003066047150000136
可以得到最终的系统状态递推修正算法：
Figure BDA0003066047150000137
步骤5：根据获取的精确系统状态，进行权重递推校正，具体如下：
在执行了状态估计之后，可获取k周期t时刻下SCADA复数伪量测和PMU量测数据的残差rk,t和ek,t如下：
Figure BDA0003066047150000138
Figure BDA0003066047150000139
在量测误差服从高斯分布时，在当前周期SCADA复数伪量测和PMU量测数据的协方差
Figure BDA00030660471500001310
和
Figure BDA00030660471500001311
的基础上，可以通过如下的公式对rk,t和ek,t下一周期SCADA复数伪量测和PMU量测数据更新对应的协方差
Figure BDA00030660471500001312
和
Figure BDA00030660471500001313
进行递推估计：
Figure BDA00030660471500001314
Figure BDA00030660471500001315
其中，med{(·)}表示取中位数操作，Le是估计窗口长度，c＝1.483(1+5/(Le-1))是采样校正系数，α是修正遗忘因子。作为一种优选方案，Le可选为10，α可选为0.95。
rk,t,rk-1,t,…,rk-Le+1,t分别表示k至k-Le+1周期下t时刻的SCADA复数伪量测的残差；ek,t,ek-1,t,…,ek-Le+1,t分别表示k至k-Le+1周期下t时刻的PMU量测数据的残差。
PMU采样率远远大于SCADA采样率，如果不加以处理，时间偏移将会造成估计精度的明显下降。理论上，时间偏移造成的误差会时间偏移的增大而增大，幸而，SCADA和PMU量测之间的时间偏移大致是呈周期变化的，对于不同的k，只要t相同，他们的时间偏移就是相同的。故时间偏移造成的负面影响可以通过对k周期t时刻的量测权重矩阵Wk,t和Vk,t的递推修正来消除，量测权重矩阵的计算式为：
Figure BDA0003066047150000141
Figure BDA0003066047150000142
其中，i表示第i个元素，即第i个量测。
步骤6：根据获取的精确系统状态，利用修正区域调整算法进行修正区域调整，具体如下：
当SCADA和PMU采样率差异过大，且PMU不可观区域的状态变化过大时，仍然会有部分节点的时间偏移误差无法完全消除，此时，可以通过修正区域调整算法来选定需要修正系统状态的区域。
当(k,1)时刻的SCADA复数伪量测和PMU量测数据到来时，可以通过如下的修正效果指标
Figure BDA0003066047150000143
来评估所有节点在k-1周期的修正效果：
Figure BDA0003066047150000144
其中，
Figure BDA0003066047150000145
表示(k,t)时刻的精确系统状态向量的第i个元素，即第i个节点，下标(·)k-1,τ表示k-1周期的最后一个时刻。此时，可以通过如下的调整指标
Figure BDA0003066047150000146
来选取k周期修正区域：
Figure BDA0003066047150000147
其中，Lm表示监控窗口的长度，Lm的值越大，随机误差对修正区域调整造成的影响越小，但同时
Figure BDA0003066047150000151
调整的灵敏度会降低，调整指标
Figure BDA0003066047150000152
的节点将会被划分至修正区域内，只有修正区域内的节点会进行系统状态递推校正。作为一种优选方案，可选取Lm为10。
步骤7：建立基于流计算的多线程处理框架，将粗略系统状态和精确系统状态的估计过程拆分为若干个并行执行的进程，获取最终的精确系统状态输出，具体如下：
为了进一步提高算法的计算效率，将粗略系统状态和精确系统状态的估计过程拆分为如图2所示的5个可并行执行的进程：
进程1：当SCADA量测数据到来时，计算SCADA复数伪量测zk,t；
进程2：根据SCADA复数伪量测zk,t，计算粗略系统状态
Figure BDA0003066047150000153
若zk,t没有计算完成，则用k-1周期t时刻的SCADA复数伪量测zk-1,t代替；
进程3：当PMU量测数据到来时，修正粗略系统状态
Figure BDA0003066047150000154
为精确系统状态
Figure BDA0003066047150000155
若
Figure BDA0003066047150000156
没有计算完成，则用k-1周期t时刻的粗略系统状态
Figure BDA0003066047150000157
代替，以获取最终的精确系统状态输出；
进程4：当精确系统状态
Figure BDA0003066047150000158
计算完成时，进行k+1周期t时刻的SCADA复数伪量测的量测权重矩阵Wk+1,t和PMU量测数据的量测权重矩阵Vk+1,t的递推校正，基于校正后的量测权重矩阵，计算下一个周期的中间系数矩阵Ak+1,t，Bk+1,t，Ck+1,t，Dk+1,t，并对这四个矩阵进行LU分解用于下一个周期其他进程的计算；
进程5：当(k,1)时刻的SCADA复数伪量测和PMU量测数据到来时，进行修正区域调整计算，获取第k周期采用的修正区域。
可以看出上述进程1、2、4、5不要求实时完成，最终系统状态更新的时间仅取决于进程3。
此外，状态估计器需要容纳并处理大规模的数据流，这些数据是实时产生，实时到达，且为了保证状态估计结果的实时性，必须要实时处理并实时反馈状态估计结果。如图3所示，传统的分层数据处理框架先将量测数据存储在硬盘的数据库中，然后才能被内存中的状态估计器等应用程序调用。然而，在硬盘上存储所有的量测数据是不现实且没有必要的，况且，硬盘和内存之间频繁的数据交互是十分耗时的，为此，本发明提出的状态估计的流计算框架如图4所示，其可分为如下几个步骤：
①从量测装置获取量测数据：SCADA和PMU量测数据；
②将SCADA和PMU量测数据通过流数据管道传递给状态估计器进行状态估计计算；
③状态估计器将状态估计结果通过状态估计器传递给后续的高级应用，也可以将重要时间断面的估计结果保存在硬盘中的数据库中。
其中，所述流数据管道作用是提供缓冲，暂存无法及时处理的数据。可以看出，主要的数据交互和处理均在内存中进行，避免了硬盘和内存之间频繁的数据交互，大大提高了程序的执行速度。
为了验证本发明方法能够快速跟踪系统状态，下面以IEEE14、118、2383节点标准系统，进行验证说明。
IEEE标准节点量测数据由潮流真值添加随机噪声得到，其中RTU功率量测及电流幅值量测噪声的标准差为0.01pu，RTU电压幅值量测噪声的标准差为0.005pu，PMU复电压量测和复电流量测的标准差为0.001pu。作为对比，本发明选取WLS算法及文献《Robusthybrid linear state estimator utilizing SCADA and PMU measurements》提出的HLSE算法与本发明提出的方法进行性能比较。估计精度和计算效率均经过100次蒙特卡洛实验得到验证，每一组实验中，PMU的安装率均为20％，有5％的负荷存在负荷波动，负荷波动曲线如图6所示，每个节点的负荷波动从5组负荷波动曲线中随机选取一组，每一组实验都有575个时间断面。
为更好地比较状态量的估计精度，采用如下的平均绝对值误差(MAE)作为评价指标:
Figure BDA0003066047150000161
其中，xesti和xtrue分别表示状态量估计值和真值，状态量为电压幅值和电压相角。
对于IEEE 14节点系统，各算法根据时间和节点号的电压幅值及电压相角估计误差如图7所示，节点9的电压幅值和相角的变化情况如图8所示，可以看出，本专利所提算法的估计精度、跟踪性能高于WLS算法和HLSE算法。图9和图10分别展示了IEEE 118节点系统中PMU安装率对估计精度和计算效率的影响，可以看出，本发明所提算法随着PMU安装率的上升，估计误差也相应的降低，但对计算效率几乎没有影响。
表1不同估计器的估计误差
Figure BDA0003066047150000162
Figure BDA0003066047150000171
表2不同估计器的计算效率
Figure BDA0003066047150000172
表1展示了三种算法在分别在IEEE14，118，2383节点系统中的估计精度，表2展示了三种算法在分别在IEEE14，118，2383节点系统中的计算时间，可以看出，本发明所提算法在提高了状态估计精度的同时，大大提高了状态估计的计算效率，保证了状态估计的实时性。
以上仿真结果验证了本发明方法的有效性和实用性。因此，本发明方法提高了状态估计的估计精度和计算效率，提高了状态量更新的频率，提高了状态估计的实时性。所提的基于流计算的多线程处理框架，使得算法可以并行执行，保证了方法的可行性，使得估计性能更加稳定。
上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (8)

1.一种递推修正混合线性状态估计方法，其特征在于，包括以下步骤：
步骤1：获得电力系统的拓扑结构、网络参数以及量测数据,其中量测数据包括SCADA量测数据和PMU量测数据；
步骤2：根据获取的SCADA量测数据计算相对于参考节点的SCADA复数伪量测；
步骤3：根据计算的SCADA复数伪量测，估计粗略系统状态；
步骤4：根据获取的PMU量测数据，利用系统状态递推修正算法对粗略系统状态进行递推修正，获取精确系统状态；
步骤5：根据获取的精确系统状态，进行权重递推校正；
步骤6：根据获取的精确系统状态，利用修正区域调整算法进行修正区域调整；
步骤7：建立基于流计算的多线程处理框架，将粗略系统状态和精确系统状态的估计过程拆分为若干个并行执行的进程，通过并行的进程，获取最终的精确系统状态。
2.根据权利要求1所述递推修正混合线性状态估计方法，其特征在于，所述步骤1中获得的网络参数包括输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、变压器变比和阻抗；以及，获得的量测数据中SCADA量测数据包括节点电压幅值、节点注入有功功率、节点注入无功功率、节点注入电流幅值、线路首端有功功率、线路首端无功功率、线路末端有功功率、线路末端无功功率、线路首端电流幅值以及线路末端电流幅值；以及，获得的量测数据中PMU量测数据包括节点复电压、线路首端复电流、线路末端复电流。
3.根据权利要求1所述递推修正混合线性状态估计方法，其特征在于，所述步骤2中计算相对于参考节点的SCADA复数伪量测，具体为：
1)建立SCADA复电压
Figure FDA0003066047140000011
量测方程：
Figure FDA0003066047140000012
其中，Uf,S是节点f的SCADA电压幅值量测值；
Figure FDA0003066047140000013
和
Figure FDA0003066047140000014
为节点f的复电压和相角指数值；
2)建立SCADA支路复电流
Figure FDA0003066047140000015
量测方程：
Figure FDA0003066047140000016
其中，
Figure FDA0003066047140000021
和
Figure FDA0003066047140000022
分别为线路f-t的并联导纳和支路导纳；
Figure FDA0003066047140000023
和
Figure FDA0003066047140000024
表示节点f和节点t的复电压，
Figure FDA0003066047140000025
为局部支路复电流伪量测；
3)建立SCADA节点复电流
Figure FDA0003066047140000026
量测方程：
Figure FDA0003066047140000027
其中，Lf表示表示和节点f通过线路直接相连的所有其他节点；
定义下标k,t表示第k次SCADA量测数据上传以来的第t次PMU量测数据上传对应的时刻，简称为(k,t)时刻或k周期的t时刻；zk,t和mk,t分别为SCADA复数伪量测和PMU量测数据，T和K分别为SCADA复数伪量测和PMU量测数据对应的量测方程系数矩阵，Wk,t和Vk,t分别为SCADA复数伪量测和PMU量测数据对应的量测权重矩阵，则相应的，SCADA复数伪量测的增益矩阵Gk,t和PMU量测数据的增益矩阵Fk,t分别定义为Gk,t＝TTWk,tT和Fk,t＝KTVk,tK；
根据以上SCADA复数伪量测方程，建立矩阵形式的量测方程：
Figure FDA0003066047140000028
其中，Lk,1是(k,1)时刻量测方程的系数矩阵，
Figure FDA0003066047140000029
是(k,1)时刻临时状态变量，εS是SCADA量测的误差向量，
Figure FDA00030660471400000222
是(k,1)时刻平衡节点的量测向量；
并且，若用
Figure FDA00030660471400000210
和δk,1表示(k,1)时刻的复电压和相角，则上述Lk,1和
Figure FDA00030660471400000211
可表示为：
Figure FDA00030660471400000212
Figure FDA00030660471400000213
其中，n表示节点数量，ES为常数阵，
Figure FDA00030660471400000214
为
Figure FDA00030660471400000215
和
Figure FDA00030660471400000216
组成的常数矩阵，
Figure FDA00030660471400000217
为导纳矩阵；使用
Figure FDA00030660471400000218
和
Figure FDA00030660471400000219
分别对应代表ES、
Figure FDA00030660471400000220
-Uf,S、
Figure FDA00030660471400000221
和
Figure FDA0003066047140000031
根据上述矩阵形式的量测方程建立线性最小二乘模型，求解可得临时状态变量为：
Figure FDA0003066047140000032
计算完成后，根据如下公式计算(k,1)时刻相对于参考节点的SCADA复数伪量测：
Figure FDA0003066047140000033
其中，
Figure FDA0003066047140000034
是Lk,1的后n-1列，
Figure FDA0003066047140000035
为相角指数量，即
Figure FDA0003066047140000036
的后n-1列。
4.根据权利要求1所述递推修正混合线性状态估计方法，其特征在于，所述步骤3中估计粗略系统状态，具体为：
利用下述线性WLS状态估计方程获取粗略系统状态：
Figure FDA0003066047140000037
其中，
Figure FDA0003066047140000038
为相对于参考节点的SCADA复数伪量测的量测矩阵，J为一个常数矩阵；
Figure FDA0003066047140000039
为粗略系统状态；Wk,t为SCADA复数伪量测的量测权重矩阵；zk,1为(k,1)时刻相对于参考节点的SCADA复数伪量测；
定义如下中间系数矩阵Ak,t＝Gk,t，Bk,t＝JTWk,t，则上述线性WLS状态估计方程表述为：
Figure FDA00030660471400000310
5.根据权利要求1所述递推修正混合线性状态估计方法，其特征在于，所述步骤4中利用系统状态递推修正算法对粗略系统状态进行递推修正，具体为：
当PMU量测数据mk,t到来后，通过如下的系统状态递推修正算法将粗略系统状态
Figure FDA00030660471400000311
修正为精确系统状态
Figure FDA00030660471400000312
Figure FDA00030660471400000313
其中，Wk,t和Vk,t分别为SCADA复数伪量测和PMU量测数据对应的量测权重矩阵，则SCADA复数伪量测的增益矩阵Gk,t和PMU量测的增益矩阵Fk,t分别定义为Gk,t＝JTWk,tJ和Fk,t＝KTVk, tK，公式中J和K分别为SCADA复数伪量测和PMU量测数据对应的量测方程系数矩阵；
定义中间系数矩阵Ck,t＝Gk,t+Fk,t，Dk,t＝KTVk,t，和
Figure FDA0003066047140000041
则上述系统状态递推修正算法表述为：
Figure FDA0003066047140000042
6.根据权利要求1所述递推修正混合线性状态估计方法，其特征在于，所述步骤5中进行权重递推校正，具体为：
获取k周期t时刻的SCADA复数伪量测和PMU量测数据的残差rk,t和ek,t如下：
Figure FDA0003066047140000043
Figure FDA0003066047140000044
其中，J和K分别为SCADA复数伪量测和PMU量测数据对应的量测方程系数矩阵，
Figure FDA0003066047140000045
为精确系统状态，zk,1为(k,1)时刻相对于参考节点的SCADA复数伪量测，mk,t为(k,1)时刻相对于参考节点的PMU量测数据；
在量测误差服从高斯分布时，在当前周期SCADA和PMU协方差
Figure FDA0003066047140000046
和
Figure FDA0003066047140000047
的基础上，通过如下的公式对rk,t和ek,t下一周期SCADA和PMU更新对应的协方差
Figure FDA0003066047140000048
和
Figure FDA0003066047140000049
进行递推估计：
Figure FDA00030660471400000410
Figure FDA00030660471400000411
其中，med{(·)}表示取中位数操作，Le是估计窗口长度，c＝1.483(1+5/(Le-1))是采样校正系数，α是修正遗忘因子；rk,t,rk-1,t,…,rk-Le+1,t分别表示k至k-Le+1周期下t时刻的SCADA复数伪量测的残差；ek,t,ek-1,t,…,ek-Le+1,t分别表示k至k-Le+1周期下t时刻的PMU量测数据的残差；
对量测权重矩阵Wk,t和Vk,t递推修正，得到量测权重矩阵的计算式为：
Figure FDA00030660471400000412
Figure FDA00030660471400000413
其中，i表示第i个元素，即第i个量测。
7.根据权利要求1所述递推修正混合线性状态估计方法，其特征在于，所述步骤6中利用修正区域调整算法进行修正区域调整，具体为：
当(k,1)时刻的SCADA复数伪量测和PMU量测数据到来时，通过如下的修正效果指标
Figure FDA0003066047140000051
来评估所有节点在k-1周期的修正效果：
Figure FDA0003066047140000052
其中，
Figure FDA0003066047140000053
表示(k,t)时刻的精确系统状态向量的第i个元素，即第i个节点，下标(·)k-1,τ表示k-1周期的最后一个时刻；此时，通过如下的调整指标
Figure FDA0003066047140000054
来选取k周期修正区域：
Figure FDA0003066047140000055
其中，Lm表示监控窗口的长度；
并且，将调整指标
Figure FDA0003066047140000056
的节点划分至修正区域内，只对修正区域内的节点进行系统状态递推校正。
8.根据权利要求1所述递推修正混合线性状态估计方法，其特征在于，所述步骤7中将粗略系统状态和精确系统状态的估计过程拆分为若干个并行执行的进程，具体为：
进程1：当SCADA量测数据到来时，计算SCADA复数伪量测zk,t；
进程2：根据SCADA复数伪量测zk,t，计算粗略系统状态
Figure FDA0003066047140000057
若zk,t没有计算完成，则用k-1周期t时刻的SCADA复数伪量测zk-1,t代替；
进程3：当PMU量测数据到来时，修正粗略系统状态
Figure FDA0003066047140000058
为精确系统状态
Figure FDA0003066047140000059
若
Figure FDA00030660471400000510
没有计算完成，则用k-1周期t时刻的粗略系统状态
Figure FDA00030660471400000511
代替，以获取最终的精确系统状态输出；
进程4：当精确系统状态
Figure FDA00030660471400000512
计算完成时，进行k+1周期t时刻的SCADA复数伪量测的量测权重矩阵Wk+1,t和PMU量测数据的量测权重矩阵Vk+1,t的递推校正，基于校正后的量测权重矩阵，计算下一个周期的中间系数矩阵Ak+1,t，Bk+1,t，Ck+1,t，Dk+1,t，并对这四个矩阵进行LU分解用于下一个周期其他进程的计算；
进程5：当(k,1)时刻的SCADA复数伪量测和PMU量测数据到来时，进行修正区域调整计算，获取第k周期采用的修正区域。
CN202110532099.XA 2021-05-14 2021-05-14 一种递推修正混合线性状态估计方法 Active CN113553538B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110532099.XA CN113553538B (zh) 2021-05-14 2021-05-14 一种递推修正混合线性状态估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110532099.XA CN113553538B (zh) 2021-05-14 2021-05-14 一种递推修正混合线性状态估计方法

Publications (2)

Publication Number Publication Date
CN113553538A true CN113553538A (zh) 2021-10-26
CN113553538B CN113553538B (zh) 2023-12-01

Family

ID=78101837

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110532099.XA Active CN113553538B (zh) 2021-05-14 2021-05-14 一种递推修正混合线性状态估计方法

Country Status (1)

Country Link
CN (1) CN113553538B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115146538A (zh) * 2022-07-11 2022-10-04 河海大学 基于消息传递图神经网络的电力系统状态估计方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103326358A (zh) * 2013-06-17 2013-09-25 西南交通大学 基于同步相角测量装置的电力系统动态状态估计方法
WO2016078477A1 (zh) * 2014-11-18 2016-05-26 国电南瑞科技股份有限公司 一种变电站三相线性广义状态估计方法
CN110299762A (zh) * 2019-06-21 2019-10-01 三峡大学 基于pmu准实时数据的主动配电网抗差估计方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103326358A (zh) * 2013-06-17 2013-09-25 西南交通大学 基于同步相角测量装置的电力系统动态状态估计方法
WO2016078477A1 (zh) * 2014-11-18 2016-05-26 国电南瑞科技股份有限公司 一种变电站三相线性广义状态估计方法
CN110299762A (zh) * 2019-06-21 2019-10-01 三峡大学 基于pmu准实时数据的主动配电网抗差估计方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115146538A (zh) * 2022-07-11 2022-10-04 河海大学 基于消息传递图神经网络的电力系统状态估计方法

Also Published As

Publication number Publication date
CN113553538B (zh) 2023-12-01

Similar Documents

Publication Publication Date Title
CN108199375B (zh) 基于同步相量量测的智能配电网拓扑辨识方法
CN107590317B (zh) 一种计及模型参数不确定性的发电机动态估计方法
CN108155648B (zh) 基于自适应h无穷扩展卡尔曼滤波的状态估计方法
CN107133406B (zh) 一种电力系统静态电压稳定域边界的快速搜索方法
CN106505557B (zh) 一种遥测错误辨识方法及装置
CN102377180B (zh) 基于电能质量监测系统的电力系统负荷建模方法
Wang et al. On multi-event co-calibration of dynamic model parameters using soft actor-critic
CN109088407B (zh) 基于深度信念网络伪量测建模的配电网状态估计方法
CN110333389A (zh) 基于插值dft的正弦信号频率估计方法
CN110543720B (zh) 基于sdae-elm伪量测模型的状态估计方法
CN107658881A (zh) 基于戴维南等值方法的电压稳定临界点判断方法
US5627760A (en) Method and apparatus for real time recursive parameter energy management system
CN109494724B (zh) 基于lu分解的大电网戴维南等值参数在线辨识方法
CN109711662B (zh) 一种基于多源数据融合的电网抗差状态估计方法
CN112149879A (zh) 一种计及宏观波动性分类的新能源中长期电量预测方法
CN115453193B (zh) 基于pqm、ttu和sm量测数据协同的配电网谐波状态估计方法
CN111046327A (zh) 适用于低频振荡与次同步振荡辨识的Prony分析方法
CN104036356A (zh) 一种利用分形算法对电网未来运行状态进行预测的方法
CN104537233B (zh) 一种基于核密度估计的配电网伪量测生成方法
CN106372440B (zh) 一种并行计算的配电网自适应抗差状态估计方法及装置
CN113553538B (zh) 一种递推修正混合线性状态估计方法
Sun et al. High-refresh-rate robust state estimation based on recursive correction for large-scale power systems
CN110021931B (zh) 一种计及模型不确定性的电力系统辅助预测状态估计方法
CN109638811B (zh) 基于模型等值的配电网电压功率灵敏度鲁棒估计方法
Liu et al. A data-driven method for online constructing linear power flow model

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant