CN110907702A - 一种改进动态谐波估计方法和系统 - Google Patents

一种改进动态谐波估计方法和系统 Download PDF

Info

Publication number
CN110907702A
CN110907702A CN201911045986.3A CN201911045986A CN110907702A CN 110907702 A CN110907702 A CN 110907702A CN 201911045986 A CN201911045986 A CN 201911045986A CN 110907702 A CN110907702 A CN 110907702A
Authority
CN
China
Prior art keywords
harmonic
measurement
covariance
representing
current
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.)
Pending
Application number
CN201911045986.3A
Other languages
English (en)
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.)
State Grid Corp of China SGCC
China Electric Power Research Institute Co Ltd CEPRI
State Grid Shanghai Electric Power Co Ltd
Original Assignee
State Grid Corp of China SGCC
China Electric Power Research Institute Co Ltd CEPRI
State Grid Shanghai Electric Power Co Ltd
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 State Grid Corp of China SGCC, China Electric Power Research Institute Co Ltd CEPRI, State Grid Shanghai Electric Power Co Ltd filed Critical State Grid Corp of China SGCC
Priority to CN201911045986.3A priority Critical patent/CN110907702A/zh
Publication of CN110907702A publication Critical patent/CN110907702A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16—Spectrum analysis; Fourier analysis
    • G01R23/165—Spectrum analysis; Fourier analysis using filters

Landscapes

  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Remote Monitoring And Control Of Power-Distribution Networks (AREA)

Abstract

本发明提供了一种改进动态谐波估计方法和系统，包括：采集电网各支路或母线的同步相量测量装置的谐波电压；将谐波电压带入非线性状态方程与量测方程模型；采用改进Sage‑Husa平方根无迹卡尔曼滤波算法，根据非线性状态方程与量测方程模型，对待估测支路进行谐波电流估算，得到待估测支路谐波电流的估算值；其中，改进Sage‑Husa平方根无迹卡尔曼滤波算法包括：基于量测更新和遗忘因子，对非线性状态方程与量测方程模型进行状态更新。本发明针对平方根无迹卡尔曼滤波进行电力系统谐波状态估计中无法从未知噪声中得到准确估计结果问题，引入Sage‑Husa滤波算法中基于极大后验估计的思想，利用Sage‑Husa滤波的遗忘因子，可以实现未知噪声实时动态估测，得到高精度谐波估计结果。

Description

一种改进动态谐波估计方法和系统
技术领域
本发明属于电力系统状态估计技术领域，具体涉及一种改进动态谐波估计方法和系统。
背景技术
谐波状态估计是近年来逐步发展起来的一项有效谐波监控与分析技术。谐波状态估计技术利用安装在系统部分母线及线路上的同步相量测量装置(Phasor MeasurementUnit,PMU)设备所提供的数据，依据相应的估计准则，推断出整个电网的支路谐波电流与节点谐波电压状态。其中，动态谐波状态估计算法是根据电力系统谐波的运动方程以及某一时刻的测量数据为初值进行下一个时刻状态量的估计算法。相比较静态谐波状态估计而言，动态谐波状态估计既可以滤波也可以预测。
卡尔曼滤波(Kalman Filter，KF)是一种高效的自回归滤波器，它能在存在诸多不确定性情况的组合信息中估计动态系统的状态。传统KF算仅适用于线性系统状态估计，针对非线性系统状态估计中的过程噪声与观测噪声估计问题，诸如扩展卡尔曼(ExtendedKalman Filter,EKF)与无迹卡尔曼(Unscented Kalman Filter,UKF)等基于KF算法等衍生算法不断被提出。但EKF算法是通过泰勒级数展开取低阶项的方式，将非线性模型近似为线性模型，该处理过程将引入近似误差。UKF相对于EKF的估计精度有一定的提高，但存在难以确定状态噪声协方差矩阵正定性的缺陷。
考虑到电力系统是一个时变的非线性系统，谐波亦是动态变化的，对谐波进行估计时系统的过程噪声方差与观测噪声方差是未知的，具有强不确定性，通常只能通过经验进行估计，而错误的参数估计又常导致滤波的发散，继而无法得到准确的估计结果。平方根无迹卡尔曼滤波(Square-Root Unscented Kalman Filter,SRUKF)利用协方差平方根代替协方差参加递推运算，可一定程度的解决由于协方差矩阵负定而导致的滤波结果发散问题，然而其本质是采用正态分布来逼近系统状态的后验概率密度，并未有效解决电力系统噪声的强不确定性将影响估计结果精度的问题。为了提高电力系统谐波状态估计鲁棒性与精确性，应当科学的改进现有状态估计算法的缺点，使得其能够滤除未知噪声得到准确的估测谐波估测。
发明内容
为克服上述现有技术的不足，本发明提出一种改进动态谐波估计方法，其改进之处在于，包括：
采集电网各支路或母线的同步相量测量装置的谐波电压；
将所述谐波电压带入预先建立的非线性状态方程与量测方程模型；
采用改进Sage-Husa平方根无迹卡尔曼滤波算法，根据所述非线性状态方程与量测方程模型，对待估测支路进行谐波电流估计计算，得到待估测支路谐波电流的估算值；
其中，所述改进Sage-Husa平方根无迹卡尔曼滤波算法包括：基于量测更新和遗忘因子，对所述非线性状态方程与量测方程模型进行状态更新。
本发明提供的第一优选技术方案，其改进之处在于，所述非线性状态方程与量测方程模型如下式所示：
Figure BDA0002254155310000021
其中，下标k表示下一采样时刻，下标k-1表示当前采样时刻；Xk表示下一采样时刻的谐波电流的估算值，Xk-1表示表示当前采样时刻的谐波电流的估算值，Uk表示下一采样时刻的谐波电压，F为系统状态函数，H为系统观测函数，Ak表示下一采样时刻的系统矩阵，Bk表示下一采样时刻的控制矩阵，Ck表示下一采样时刻的观测矩阵，Dk表示下一采样时刻的直接转递矩阵，Wk表示下一采样时刻的系统过程噪声，Vk表示下一采样时刻的系统观测噪声，Zk表示下一采样时刻的谐波电压观测值。
本发明提供的第二优选技术方案，其改进之处在于，所述用改进Sage-Husa平方根无迹卡尔曼滤波算法，根据所述非线性状态方程与量测方程模型，对预设待估测支路进行谐波电流估计计算，得到待估测支路谐波电流的估算值，包括：
基于所述非线性状态方程与量测方程模型，采用无迹卡尔曼滤波算法进行时间更新；
基于时间更新的结果和Sage-Husa方法对所述非线性状态方程与量测方程模型进行量测更新；所述非线性状态方程与量测方程模型中的量测量序列的初始值由初次采样时刻采集的谐波电压初始化得到；
基于量测更新的结果和预设遗忘因子，对所述非线性状态方程与量测方程模型进行状态更新，得到谐波电流的估算值。
本发明提供的第三优选技术方案，其改进之处在于，初始化非线性状态方程与量测方程模型中量测量序列的初始值，如下式所示：
Figure BDA0002254155310000031
其中，下标0表示的初始采样时刻，X0表示谐波电流的初值，Q0表示初始采样时刻系统过程噪声的协方差向量，R0表示初始采样时刻系统观测噪声的协方差向量，
Figure BDA0002254155310000032
为Q0的平方根，
Figure BDA0002254155310000033
为R0的平方根，
Figure BDA0002254155310000034
是X0的数学期望；E表示求取数学期望，Chol函数为矩阵Cholesky因子分解函数，S0是
Figure BDA0002254155310000035
和X0的协方差；
所述谐波电流的初值X0与所述谐波电压的关系式如下：
Z0＝H[X0,U0]＝C0X0+D0U0+V0
其中，Z0表示初始采样时刻经过量测变换后的谐波电压，U0表示初始采样时刻的谐波电压，V0表示初始采样时刻的系统观测噪声，H为观测函数，C为观测矩阵，D为直接转递矩阵。
本发明提供的第四优选技术方案，其改进之处在于，所述基于所述非线性状态方程与量测方程模型，采用无迹卡尔曼滤波算法进行时间更新，包括：
对初始化后的所述非线性状态方程与量测方程模型中的谐波电流进行无迹变换，得到包含Sigma点集的Sigma矩阵；
根据所述非线性状态方程对所述Sigma矩阵进行非线性变换，并利用系统过程噪声的协方差的平方根参与无迹卡尔曼滤波算法进行递推运算，对所述谐波电流与谐波电流协方差进行预测。
本发明提供的第五优选技术方案，其改进之处在于，所述递推运算的计算式如下：
Figure BDA0002254155310000036
式中，
Figure BDA0002254155310000041
表示下一采样时刻的谐波电流中间值，χi,k-1表示当前采样时刻的谐波电流对应的第i个Sigma点的值，i的取值范围为0到2n，n为谐波电流的维度数，Uk-1表示当前采样时刻的谐波电压，Ak-1表示当前采样时刻的系统矩阵，Bk-1表示当前采样时刻的控制矩阵，F为系统状态函数，
Figure BDA0002254155310000042
表示第i个Sigma点在计算均值时的权重，
Figure BDA0002254155310000043
表示第i个Sigma点在计算协方差时的权重，χ* i,k-1表示对应χi,k-1的中间过程变量，Sxx,k表示下一采样时刻的谐波电流的协方差，Qk表示下一采样时刻的系统过程噪声协方差向量，
Figure BDA0002254155310000045
表示第0个Sigma点在计算协方差时的权重，χ* 0,k-1表示当前采样时刻的谐波电流对应的第0个Sigma点对应的中间过程变量，
Figure BDA0002254155310000046
表示对应Sxx,k的中间过程变量，qr为正交三角矩阵分解运算，Cholupdata为Cholesky分解运算；
当前采样时刻的谐波电流对应的第i个Sigma点的值χi,k-1如下式所示：
Figure BDA0002254155310000047
式中，σk-1表示当前采样时刻的无迹变化参数，
Figure BDA0002254155310000048
表示当前采样时刻的谐波电流中间值；
所述当前采样时刻的无迹变化参数的计算式如下：
Figure BDA0002254155310000049
式中，Sxx,k-1表示表示当前采样时刻的谐波电流的协方差，λ为缩放比例系数。
本发明提供的第六优选技术方案，其改进之处在于，所述基于时间更新的结果和预设遗忘因子，基于Sage-Husa方法对所述非线性状态方程与量测方程模型进行量测更新，包括：
基于时间更新的结果，对Sigma点进行重采样，得到Sigma点；
通过测量方程对Sigma点进行非线性变换并计算测量残差；
根据所述测量残差并考虑预设遗忘因子，基于Sage-Husa采用平方根无迹卡尔曼滤波方法对系统观测噪声协方差的平方根进行更新；
基于系统观测噪声协方差的平方根，对谐波电压和谐波电流间的协方差以及谐波电压的协方差进行更新；
根据谐波电压和谐波电流间的协方差以及谐波电压的协方差，计算卡尔曼滤波增益。
本发明提供的第七优选技术方案，其改进之处在于，所述通过测量方程对重采样的Sigma点集进行非线性变换并计算测量残差，如下式所示：
Figure BDA0002254155310000051
式中，H为观测函数，χi,k表示重采样的下一采样时刻的谐波电流对应的第i个Sigma点的值，Uk-1表示当前采样时刻的谐波电压，
Figure BDA0002254155310000052
第i个Sigma点在计算均值时的权重，
Figure BDA0002254155310000053
表示下一采样时刻的添加权重后修正后且经量测变换后的谐波电压，Zk表示对应下一采样时刻采集的谐波电压观测值，ek表示下一采样时刻谐波电压的量测残差，Z* i,k表示对应
Figure BDA0002254155310000054
的第i个Sigma点的中间过程变量。
本发明提供的第八优选技术方案，其改进之处在于，所述根据所述测量残差并考虑预设遗忘因子，基于Sage-Husa采用平方根无迹卡尔曼滤波方法对系统观测噪声协方差的平方根进行更新，如下式所示：
Figure BDA0002254155310000055
式中，
Figure BDA0002254155310000056
表示下一采样时刻系统观测噪声协方差的平方根，R** k表示对应
Figure BDA0002254155310000057
的第一中间过程变量，R* k表示对应
Figure BDA0002254155310000058
的第二中间过程变量，
Figure BDA0002254155310000059
表示当前采样时刻系统观测噪声协方差的平方根，ek表示下一采样时刻谐波电压的量测残差，dk为根据遗忘因子计算的下一采样时刻的遗忘修正因子，
Figure BDA00022541553100000510
表示第i个Sigma点在计算协方差时的权重，
Figure BDA00022541553100000511
表示下一采样时刻的添加权重后修正后第i个Sigma点且经量测变换后的谐波电压，
Figure BDA00022541553100000512
表示当前采样时刻的添加权重后修正后且经量测变换后第i个Sigma点的谐波电压，
Figure BDA00022541553100000513
表示R* k的转置，Cholupdata为Cholesky分解运算，diag为构造对角矩阵运算；
下一采样时刻的遗忘修正因子dk的计算式如下：
dk＝(1-b)/(1-bk+1)
式中，b为遗忘因子。
本发明提供的第九优选技术方案，其改进之处在于，所述基于系统观测噪声协方差的平方根，对谐波电压和谐波电流间的协方差以及谐波电压的协方差进行更新，如下式所示：
Figure BDA0002254155310000061
式中，Pxz,k表示下一采样时刻谐波电压和谐波电流间的协方差，Szz,k表示下一采样时刻谐波电压的协方差，S* zz,k表示对应Szz,k的中间过程变量，
Figure BDA0002254155310000062
表示下一采样时刻的添加权重后修正后且经量测变换后的谐波电压，Z* i,k表示对应
Figure BDA0002254155310000063
的第i个Sigma点的中间过程变量，
Figure BDA0002254155310000064
表示下一采样时刻得到的谐波电流中间值，χ* i,k表示对应χi,k的中间过程变量，χi,k表示下一采样时刻的谐波电流对应的第i个Sigma点的值，
Figure BDA00022541553100000611
表示第i个Sigma点在计算协方差时的权重，R* k表示对应
Figure BDA0002254155310000065
的第二中间过程变量，
Figure BDA0002254155310000066
表示第0个Sigma点在计算协方差时的权重，Z* 0,k-1表示对应
Figure BDA0002254155310000067
的第0个Sigma点的中间过程变量，
Figure BDA0002254155310000068
表示当前采样时刻的添加权重后修正后且经量测变换后的谐波电压，qr为正交三角矩阵分解运算，Cholupdata为Cholesky分解运算。
本发明提供的第十优选技术方案，其改进之处在于，所述卡尔曼滤波增益的计算式如下：
Figure BDA0002254155310000069
式中，Kk表示下一采样时刻的卡尔曼滤波增益，ST zz,k表示Szz,k的转置。
本发明提供的第十一优选技术方案，其改进之处在于，所述基于量测更新的结果，对所述非线性状态方程与量测方程模型进行状态更新，得到谐波电流的估算值，包括：
根据所述非线性状态方程与量测方程模型中的谐波电流中间值、测量残差和卡尔曼滤波增益进行计算，得到谐波电流的估算值。
本发明提供的第十二优选技术方案，其改进之处在于，所述谐波电流的估算值的计算式如下：
Figure BDA00022541553100000610
其中，Xk表示下一采样时刻的谐波电流的估算值，
Figure BDA0002254155310000071
表示下一采样时刻得到的谐波电流中间值，Kk表示下一采样时刻的卡尔曼滤波增益，ek表示下一采样时刻谐波电压的量测残差。
本发明提供的第十三优选技术方案，其改进之处在于，所述得到谐波电流的估算值之后，还包括：
对所述谐波电流的协方差估计值进行更新，以及对系统过程噪声协方差向量的平方根进行更新；
将当前采样时刻非线性状态方程与量测方程模型的参数传递至下一采样时刻的非线性状态方程与量测方程模型。
本发明提供的第十四优选技术方案，其改进之处在于，所述谐波电流的协方差估计值的计算式如下：
Sk＝cholupdate{Sk-1,KkSzz,k,-1}
其中，Sk表示下一采样时刻谐波电流的协方差估计值，Sk-1表示当前采样时刻的谐波电流的协方差估计值，Kk表示下一采样时刻的卡尔曼滤波增益，Szz,k表示下一采样时刻谐波电压的协方差，Cholupdata为Cholesky分解运算；
所述系统过程噪声协方差向量的平方根的计算式如下：
Figure BDA0002254155310000072
其中，
Figure BDA0002254155310000073
表示下一采样时刻的系统过程噪声协方差向量的平方根，Qk-1表示第k-1次递推后的系统过程噪声协方差向量，dk为下一采样时刻的遗忘修正因子，Kk表示下一采样时刻的卡尔曼滤波增益，ek表示下一采样时刻谐波电压的量测残差，Szz,k表示第k次递推后谐波电压的协方差，Q** k表示对应
Figure BDA0002254155310000074
的第一中间过程变量，Q* k表示对应
Figure BDA0002254155310000075
的第二中间过程变量，Cholupdata为Cholesky分解运算，diag为构造对角矩阵运算。
本发明提供的第十五优选技术方案，其改进之处在于，所述将当前采样时刻非线性状态方程与量测方程模型的参数传递至下一采样时刻的非线性状态方程与量测方程模型，包括：
将当前采样时刻的谐波电流中间值作为下一采样时刻非线性状态方程与量测方程模型的谐波电流中间值，将当前采样时刻的谐波电流的协方差估计值作为下一采样时刻非线性状态方程与量测方程模型的谐波电流的协方差估计值，将当前采样时刻的系统过程噪声协方差向量的平方根作为下一采样时刻非线性状态方程与量测方程模型的谐波电流的协方差估计值，将当前采样时刻的系统观测噪声协方差的平方根作为下一采样时刻非线性状态方程与量测方程模型的系统观测噪声协方差的平方根。
基于同一发明构思，本发明还提供了一种改进动态谐波估计系统，其改进之处在于，包括：数据采集模块、数据输入模块和谐波估计模块；
所述数据采集模块，用于采集电网各支路或母线的同步相量测量装置的谐波电压；
所述数据输入模块，用于将所述谐波电压带入预先建立的非线性状态方程与量测方程模型；
所述谐波估计模块，用于采用改进Sage-Husa平方根无迹卡尔曼滤波算法，根据所述非线性状态方程与量测方程模型，对待估测支路进行谐波电流估计计算，得到待估测支路谐波电流的估算值；
其中，所述改进Sage-Husa平方根无迹卡尔曼滤波算法包括：基于量测更新和遗忘因子，对所述非线性状态方程与量测方程模型进行状态更新。
本发明提供的第十六优选技术方案，其改进之处在于，所述谐波估计模块包括：时间更新模块、量测更新模块和状态更新模块；
所述时间更新模块，用于基于所述非线性状态方程与量测方程模型，采用无迹卡尔曼滤波算法进行时间更新；所述非线性状态方程与量测方程模型中的量测量序列的初始值由初次采样时刻采集的谐波电压初始化得到；
所述量测更新模块，用于基于时间更新的结果和Sage-Husa方法对所述非线性状态方程与量测方程模型进行量测更新；
所述状态更新模块，用于基于量测更新的结果和预设遗忘因子，对所述非线性状态方程与量测方程模型进行状态更新，得到谐波电流的估算值。
与最接近的现有技术相比，本发明具有的有益效果如下：
本发明提供了一种改进动态谐波估计方法和系统，包括：采集电网各支路或母线的同步相量测量装置的谐波电压；将谐波电压带入预先建立的非线性状态方程与量测方程模型；采用改进Sage-Husa平方根无迹卡尔曼滤波算法，根据非线性状态方程与量测方程模型，对待估测支路进行谐波电流估计计算，得到待估测支路谐波电流的估算值；其中，改进Sage-Husa平方根无迹卡尔曼滤波算法包括：基于量测更新和遗忘因子，对非线性状态方程与量测方程模型进行状态更新。本发明针对采用平方根无迹卡尔曼滤波进行电力系统谐波状态估计中无法从未知噪声中得到准确估计结果问题，引入Sage-Husa滤波算法中基于极大后验估计的思想，提出的改进Sage-Husa平方根无迹卡尔曼滤波谐波动态估计算法，利用Sage-Husa滤波中的遗忘因子，可以实现未知噪声实时动态估测，得到高精度谐波估计结果。
本发明提供的改进动态谐波估计方法和系统，可进一步通过递推滤波不断修正过程噪声与观测噪声的统计特性，从而得到更高精确谐波估计结果，并实现谐波状态的动态跟踪。
附图说明
图1为本发明提供的一种改进动态谐波估计方法流程示意图；
图2为本发明提供的一种改进动态谐波估计方法简化流程示意图；
图3为本发明提供的一种改进Sage-Husa平方根无迹卡尔曼滤波实现状态估计的通用流程图；
图4为本发明提供的一种改进动态谐波估计方法详细流程示意图；
图5为本发明提供的一种改进动态谐波估计系统基本结构示意图；
图6为本发明提供的一种改进动态谐波估计系统详细本结构示意图。
具体实施方式
下面结合附图对本发明的具体实施方式做进一步的详细说明。
实施例1：
本发明提供的一种改进动态谐波估计方法流程示意图如图1所示，包括：
步骤1：采集电网各支路或母线的同步相量测量装置的谐波电压；
步骤2：将谐波电压带入预先建立的非线性状态方程与量测方程模型；
步骤3：采用改进Sage-Husa平方根无迹卡尔曼滤波算法，根据非线性状态方程与量测方程模型，对待估测支路进行谐波电流估计计算，得到待估测支路谐波电流的估算值；
其中，改进Sage-Husa平方根无迹卡尔曼滤波算法包括：基于量测更新和遗忘因子，对非线性状态方程与量测方程模型进行状态更新。
具体的，本实施例中提供的一种改进动态谐波估计方为一种改进Sage-Husa平方根无迹卡尔曼滤波的动态谐波估计方法，其实现步骤如图2所示，包括：
步骤101：依据电网拓扑结构建立非线性状态方程与量测方程模型；
步骤102：输入当前时刻量测数据；
步骤103：采用改进Sage-Husa平方根无迹卡尔曼滤波算法进行动态谐波估计计算；
步骤104：参数传递，重复步骤102至步骤104。
本发明中，一种改进Sage-Husa平方根无迹卡尔曼滤波实现状态估计的通用流程图如图3所示，将该方法运用到具体的动态谐波估计上提供的一种改进动态谐波估计方法详细流程如图4所示。
所述的一种改进Sage-Husa平方根无迹卡尔曼滤波的动态谐波估计方法，步骤101中的量测方程为建立非线性谐波量测方程，依据电力系统网络架构，设置量测量为母线谐波电压，状态估计量为待估计的节点注入各次谐波电流，状态方程与量测方程的一般形式为：
Figure BDA0002254155310000101
式中，I为谐波电流状态估计序列；
Y为导纳矩阵函数；
Ia为额外谐波电流注入向量，通常为变压器励磁电流；
Z为量测变换后的谐波电压量测向量；
U为谐波电压量测向量；
M为谐波电压量测矩阵函数；
ε、η分别为过程噪声向量与观测噪声向量；
下标k与k-1分布代表第k和第k-1采样时刻，或者说第k和第k-1次递推。
所述的一种改进Sage-Husa平方根无迹卡尔曼滤波的动态谐波估计方法，在应用改进Sage-Husa平方根无迹卡尔曼滤波进行动态谐波估计前，需将量测方程(1)改下为如式(15)所示的离散系统标准模型：
Figure BDA0002254155310000102
其中，X和Z分别为n维状态变量和m维量测变量。F和H分别为系统状态函数和观测函数。电力系统为非线性系统，因此F和H为非线性函数。W和V分别为n维系统过程噪声与m维系统观测噪声序列。A为系统矩阵，B为控制矩阵，C为观测矩阵，D为直接转递矩阵。
本实施例提供的一种改进Sage-Husa平方根无迹卡尔曼滤波的动态谐波估计方法，当电力系统网络架构、支路传输线路参数与负荷确定后，需选取待估测的谐波电流构成式(15)中的状态变量序列Xk；选取合适的母线电压构成式(15)中的量测量序列Zk；函数Y为Zk映射到Xk的导纳矩阵。配置完成的导纳矩阵需满足其秩为满秩，以及rank(Y)＝n，若不满足，则需调整量测量序列，直到满足满秩要求为止。其中，n为状态向量的维度数目。为即本实施例中，状态变量序列Xk即为波电流状态估计序列Ik。
一种改进Sage-Husa平方根无迹卡尔曼滤波的动态谐波估计方法，量测数据来源为布置在电力系统网络不同支路或母线的PMU量测装置。通常电力系统网络不会在所有支路或母线均布置PMU量测装置，因此需要通过建立系统量测方程，来估计未布置PMU量测装置的支路或母线谐波状态。
本实施例的一种改进Sage-Husa平方根无迹卡尔曼滤波的动态谐波估计方法，步骤102中改进Sage-Husa平方根无迹卡尔曼滤波的实现过程为：
a.初始化。
Figure BDA0002254155310000111
式中，下标0表示的初始采样时刻，X0表示状态向量即谐波电流的初值，Q0表示初始采样时刻系统过程噪声W0的协方差向量，R0表示初始采样时刻系统观测噪声V0的协方差向量，
Figure BDA0002254155310000112
为Q0的平方根，
Figure BDA0002254155310000113
为R0的平方根，
Figure BDA0002254155310000114
是X0的数学期望；E表示求取数学期望，Chol函数为矩阵Cholesky因子分解函数，S0是
Figure BDA0002254155310000115
和X0的协方差。
初始化仅在k＝0时的初始采样时刻进行，当k>0时，不进行初始化。
b.时间更新：对于k＝1,2，……进行迭代。
首先进行无迹变换，无迹变换包含Sigma计算与其权重的计算两部分。计算Sigma采样点即Sigma点的值:设采样点维数为状态量的维度n,需计算2n+1个Sigma采样点，计算方法如式(3)所示。
Figure BDA0002254155310000121
式中，χi表示第i个Sigma点的值，
Figure BDA0002254155310000122
和S分别为n维状态向量X的均值和协方差，λ为缩放比例系数，取λ＝α2(n+ρ)-n。α控制Sigma点的分布状态，取值范围为(10-4,1)。ρ尺度参数，取值的前提是保证(n+λ)S为半正定矩阵。
计算Sigma点的权值ω：
Figure BDA0002254155310000123
式中，上标m与c分别表示均值计算和协方差计算时用到的不同权值，下标i表示维数。参数β≥0，取值条件视是否需考虑高阶项影响而定。
经过无迹变换后，得到包含Sigma点集的矩阵：
Figure BDA0002254155310000124
其中，σk-1表示第k-1次递推后的无迹变化参数，
Figure BDA0002254155310000125
式中Sxx,k-1表示表示第k-1次递推后的谐波电流的协方差。
根据系统公式(15)对Sigma矩阵进行非线性变换，并对状态量与协方差进行一步预测，取系统过程噪声协方差的平方根代替普通UKF算法中的协方差参加递推运算。
Figure BDA0002254155310000126
式中，
Figure BDA0002254155310000131
表示第k次递推后得到的谐波电流中间值，χi,k-1表示第k-1次递推后的谐波电流对应的第i个Sigma点的值，i的取值范围为0到2n，n为谐波电流的维度数，Uk-1表示第k-1次递推后的谐波电压，Ak-1表示第k-1次递推后的系统矩阵，Bk-1表示第k-1次递推后的控制矩阵，F为系统状态函数，
Figure BDA0002254155310000132
表示第i个Sigma点在计算均值时的权重，
Figure BDA0002254155310000133
表示第i个Sigma点在计算协方差时的权重，χ* i,k-1表示对应χi,k-1的中间过程变量，Sxx,k表示第k次递推后的谐波电流的协方差，Qk表示第k次递推后的系统过程噪声协方差向量，
Figure BDA0002254155310000134
表示第0个Sigma点在计算协方差时的权重，χ* 0，k-1表示第k-1次递推后的谐波电流对应的第0个Sigma点对应的中间过程变量，
Figure BDA0002254155310000136
表示对应Sxx,k的中间过程变量，qr为正交三角矩阵分解运算，Cholupdata为Cholesky分解运算。
c.测量更新。
首先对第k个采样周期也即第k次递推的Sigma点进行重采样。
Figure BDA0002254155310000137
其中，
Figure BDA0002254155310000138
通过测量方程对Sigma点集进行非线性变换并计算测量残差ek：
Figure BDA0002254155310000139
式中，H为观测函数，χi,k表示第k次递推后的谐波电流对应的第i个Sigma点的值，Uk-1表示第k-1次递推后的谐波电压，
Figure BDA00022541553100001310
第i个Sigma点在计算均值时的权重，
Figure BDA00022541553100001311
表示第k次递推后的添加权重后修正后且经量测变换后的谐波电压，Zk表示对应第k次递推后采集的谐波电压观测值，ek表示第k次递推后谐波电压的量测残差，Z* i,k表示对应
Figure BDA00022541553100001312
的第i个Sigma点的中间过程变量。
引入Sage-Husa滤波思想，在更新估计测量噪声统计特性时加入遗忘因子，对SRUKF进行改进，系统观测噪声协方差平方根
Figure BDA00022541553100001313
更新为
Figure BDA0002254155310000141
式中，
Figure BDA0002254155310000142
表示第k次递推后系统观测噪声协方差的平方根，R** k表示对应
Figure BDA0002254155310000143
的第一中间过程变量，R* k表示对应
Figure BDA0002254155310000144
的第二中间过程变量，
Figure BDA0002254155310000145
表示第k-1次递推后系统观测噪声协方差的平方根，ek表示第k次递推后谐波电压的量测残差，dk为第k次递推后的遗忘修正因子，
Figure BDA0002254155310000146
表示第i个Sigma点在计算协方差时的权重，
Figure BDA0002254155310000147
表示第k次递推后的添加权重后修正后第i个Sigma点且经量测变换后的谐波电压，
Figure BDA0002254155310000148
表示第k-1次递推后的添加权重后修正后且经量测变换后第i个Sigma点的谐波电压，
Figure BDA00022541553100001418
表示R* k的转置，Cholupdata为Cholesky分解运算，diag为构造对角矩阵运算；dk＝(1-b)/(1-bk+1)，b为遗忘因子，取值范围为0.9到1。
更新量测量与状态量的协方差Pxz,k与量测量方差矩阵Szz,k：
Figure BDA0002254155310000149
式中，Pxz,k表示第k次递推后谐波电压和谐波电流间的协方差，Szz,k表示第k次递推后谐波电压的协方差，S* zz,k表示对应Szz,k的中间过程变量，
Figure BDA00022541553100001410
表示第k次递推后的添加权重后修正后且经量测变换后的谐波电压，Z* i,k表示对应
Figure BDA00022541553100001411
的第i个Sigma点的中间过程变量，
Figure BDA00022541553100001412
表示第k次递推后得到的谐波电流中间值，χ* i,k表示对应χi,k的中间过程变量，χi,k表示第k次递推后的谐波电流对应的第i个Sigma点的值，
Figure BDA00022541553100001413
表示第i个Sigma点在计算协方差时的权重，R* k表示对应
Figure BDA00022541553100001414
的第二中间过程变量，
Figure BDA00022541553100001415
表示第0个Sigma点在计算协方差时的权重，Z* 0,k-1表示对应
Figure BDA00022541553100001416
的第0个Sigma点的中间过程变量，
Figure BDA00022541553100001417
表示第k-1次递推后的添加权重后修正后且经量测变换后的谐波电压，qr为正交三角矩阵分解运算，Cholupdata为Cholesky分解运算。
计算卡尔曼滤波增益Kk：
Figure BDA0002254155310000151
式中，Kk表示第k次递推后的卡尔曼滤波增益，ST zz,k表示Szz,k的转置。
d.状态更新。
估算校正后的状态值Xk和状态值的协方差估计值Sk：
Figure BDA0002254155310000152
其中，Xk表示第k次递推后的谐波电流的估算值，
Figure BDA0002254155310000153
表示第k次递推后得到的谐波电流中间值，Sk表示第k次递推后谐波电流的协方差估计值，Sk-1表示第k-1次递推后的谐波电流的协方差估计值。
新估测过程噪声统计特性
Figure BDA0002254155310000154
Figure BDA0002254155310000155
其中，
Figure BDA0002254155310000156
表示第k次递推后的系统过程噪声协方差向量的平方根，Qk-1表示第k-1次递推后的系统过程噪声协方差向量，dk为第k次递推后的遗忘修正因子，Kk表示第k次递推后的卡尔曼滤波增益，ek表示第k次递推后谐波电压的量测残差，Szz,k表示第k次递推后谐波电压的协方差，Q** k表示对应
Figure BDA0002254155310000157
的第一中间过程变量，Q* k表示对应
Figure BDA0002254155310000158
的第二中间过程变量。
步骤104中，参数传递如式(14)所示：
Figure BDA0002254155310000159
即新的时刻k+1到来后，将上一时刻k的谐波电流中间值作为当前时刻k+1非线性状态方程与量测方程模型的谐波电流中间值，将上一时刻k的谐波电流的协方差估计值作为当前时刻k+1非线性状态方程与量测方程模型的谐波电流的协方差估计值，将上一时刻k的系统过程噪声协方差向量的平方根作为当前时刻k+1非线性状态方程与量测方程模型的谐波电流的协方差估计值，将上一时刻k的系统观测噪声协方差的平方根作为当前时刻k+1非线性状态方程与量测方程模型的系统观测噪声协方差的平方根。
实施例2：
基于同一发明构思，本发明还提供了一种改进动态谐波估计系统，由于这些设备解决技术问题的原理与改进动态谐波估计方法相似，重复之处不再赘述。
该系统基本结构如图5所示，包括：数据采集模块、数据输入模块和谐波估计模块；
数据采集模块，用于采集电网各支路或母线的同步相量测量装置的谐波电压；
数据输入模块，用于将谐波电压带入预先建立的非线性状态方程与量测方程模型；
谐波估计模块，用于采用改进Sage-Husa平方根无迹卡尔曼滤波算法，根据非线性状态方程与量测方程模型，对待估测支路进行谐波电流估计计算，得到待估测支路谐波电流的估算值；
其中，改进Sage-Husa平方根无迹卡尔曼滤波算法包括：基于量测更新和遗忘因子，对非线性状态方程与量测方程模型进行状态更新。
一种改进动态谐波估计系统详细结构如图6所示。
其中，谐波估计模块包括：时间更新模块、量测更新模块和状态更新模块；
时间更新模块，用于基于非线性状态方程与量测方程模型，采用无迹卡尔曼滤波算法进行时间更新；非线性状态方程与量测方程模型中的量测量序列的初始值由初次采样时刻采集的谐波电压初始化得到；
量测更新模块，用于基于时间更新的结果和Sage-Husa方法对非线性状态方程与量测方程模型进行量测更新；
状态更新模块，用于基于量测更新的结果和预设遗忘因子，对非线性状态方程与量测方程模型进行状态更新，得到谐波电流的估算值。
其中，时间更新模块包括：无迹变换单元和时间更新单元；
无迹变换单元，用于对初始化后的非线性状态方程与量测方程模型中的谐波电流进行无迹变换，得到包含Sigma点集的Sigma矩阵；
时间更新单元，用于根据非线性状态方程对Sigma矩阵进行非线性变换，并利用系统过程噪声的协方差的平方根参与无迹卡尔曼滤波算法进行递推运算，对谐波电流与谐波电流协方差进行预测。
其中，量测更新模块包括：重采样单元、测量残差单元、观测噪声更新单元、量测更新单元和滤波增益单元；
重采样单元，用于基于时间更新的结果，对Sigma点进行重采样，得到Sigma点；
测量残差单元，用于通过测量方程对Sigma点进行非线性变换并计算测量残差；
观测噪声更新单元，用于根据测量残差并考虑预设遗忘因子，基于Sage-Husa采用平方根无迹卡尔曼滤波方法对系统观测噪声协方差的平方根进行更新；
量测更新单元，用于基于系统观测噪声协方差的平方根，对谐波电压和谐波电流间的协方差以及谐波电压的协方差进行更新；
滤波增益单元，用于根据谐波电压和谐波电流间的协方差以及谐波电压的协方差，计算卡尔曼滤波增益。
其中，该系统还包括电流及过程噪声更新模块和参数传递模块；
电流及过程噪声更新模块，用于对谐波电流的协方差估计值进行更新，以及对系统过程噪声协方差向量的平方根进行更新；
参数传递模块，用于将当前采样时刻非线性状态方程与量测方程模型的参数传递至下一采样时刻的非线性状态方程与量测方程模型。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
最后应当说明的是:以上实施例仅用于说明本申请的技术方案而非对其保护范围的限制,尽管参照上述实施例对本申请进行了详细的说明,所属领域的普通技术人员应当理解:本领域技术人员阅读本申请后依然可对申请的具体实施方式进行种种变更、修改或者等同替换，但这些变更、修改或者等同替换，均在申请待批的权利要求保护范围之内。

Claims (18)

1.一种改进动态谐波估计方法，其特征在于，包括：
采集电网各支路或母线的同步相量测量装置的谐波电压；
将所述谐波电压带入预先建立的非线性状态方程与量测方程模型；
采用改进Sage-Husa平方根无迹卡尔曼滤波算法，根据所述非线性状态方程与量测方程模型，对待估测支路进行谐波电流估计计算，得到待估测支路谐波电流的估算值；
其中，所述改进Sage-Husa平方根无迹卡尔曼滤波算法包括：基于量测更新和遗忘因子，对所述非线性状态方程与量测方程模型进行状态更新。
2.如权利要求1所述的方法，其特征在于，所述非线性状态方程与量测方程模型如下式所示：
Figure FDA0002254155300000011
其中，下标k表示下一采样时刻，下标k-1表示当前采样时刻；Xk表示下一采样时刻的谐波电流的估算值，Xk-1表示表示当前采样时刻的谐波电流的估算值，Uk表示下一采样时刻的谐波电压，F为系统状态函数，H为系统观测函数，Ak表示下一采样时刻的系统矩阵，Bk表示下一采样时刻的控制矩阵，Ck表示下一采样时刻的观测矩阵，Dk表示下一采样时刻的直接转递矩阵，Wk表示下一采样时刻的系统过程噪声，Vk表示下一采样时刻的系统观测噪声，Zk表示下一采样时刻的谐波电压观测值。
3.如权利要求1所述的方法，其特征在于，所述采用改进Sage-Husa平方根无迹卡尔曼滤波算法，根据所述非线性状态方程与量测方程模型，对预设待估测支路进行谐波电流估计计算，得到待估测支路谐波电流的估算值，包括：
基于所述非线性状态方程与量测方程模型，采用无迹卡尔曼滤波算法进行时间更新；所述非线性状态方程与量测方程模型中的量测量序列的初始值由初次采样时刻采集的谐波电压初始化得到；
基于时间更新的结果和Sage-Husa方法对所述非线性状态方程与量测方程模型进行量测更新；
基于量测更新的结果和预设遗忘因子，对所述非线性状态方程与量测方程模型进行状态更新，得到谐波电流的估算值。
4.如权利要求3所述的方法，其特征在于，初始化非线性状态方程与量测方程模型中量测量序列的初始值，如下式所示：
Figure FDA0002254155300000021
其中，下标0表示的初始采样时刻，X0表示谐波电流的初值，Q0表示初始采样时刻系统过程噪声的协方差向量，R0表示初始采样时刻系统观测噪声的协方差向量，
Figure FDA0002254155300000022
为Q0的平方根，
Figure FDA0002254155300000023
为R0的平方根，
Figure FDA0002254155300000024
是X0的数学期望；E表示求取数学期望，Chol函数为矩阵Cholesky因子分解函数，S0是
Figure FDA0002254155300000025
和X0的协方差；
所述谐波电流的初值X0与所述谐波电压的关系式如下：
Z0＝H[X0,U0]＝C0X0+D0U0+V0
其中，Z0表示初始采样时刻经过量测变换后的谐波电压，U0表示初始采样时刻的谐波电压，V0表示初始采样时刻的系统观测噪声，H为观测函数，C为观测矩阵，D为直接转递矩阵。
5.如权利要求3所述的方法，其特征在于，所述基于所述非线性状态方程与量测方程模型，采用无迹卡尔曼滤波算法进行时间更新，包括：
对初始化后的所述非线性状态方程与量测方程模型中的谐波电流进行无迹变换，得到包含Sigma点集的Sigma矩阵；
根据所述非线性状态方程对所述Sigma矩阵进行非线性变换，并利用系统过程噪声的协方差的平方根参与无迹卡尔曼滤波算法进行递推运算，对所述谐波电流与谐波电流协方差进行预测。
6.如权利要求5所述的方法，其特征在于，所述递推运算的计算式如下：
Figure FDA0002254155300000026
式中，
Figure FDA0002254155300000031
表示下一采样时刻的谐波电流中间值，χi,k-1表示当前采样时刻的谐波电流对应的第i个Sigma点的值，i的取值范围为0到2n，n为谐波电流的维度数，Uk-1表示当前采样时刻的谐波电压，Ak-1表示当前采样时刻的系统矩阵，Bk-1表示当前采样时刻的控制矩阵，F为系统状态函数，
Figure FDA0002254155300000032
表示第i个Sigma点在计算均值时的权重，
Figure FDA0002254155300000033
表示第i个Sigma点在计算协方差时的权重，χ* i,k-1表示对应χi,k-1的中间过程变量，Sxx,k表示下一采样时刻的谐波电流的协方差，Qk表示下一采样时刻的系统过程噪声协方差向量，
Figure FDA0002254155300000034
表示第0个Sigma点在计算协方差时的权重，χ* 0,k-1表示当前采样时刻的谐波电流对应的第0个Sigma点对应的中间过程变量，
Figure FDA0002254155300000035
表示对应Sxx,k的中间过程变量，qr为正交三角矩阵分解运算，Cholupdata为Cholesky分解运算；
当前采样时刻的谐波电流对应的第i个Sigma点的值χi,k-1如下式所示：
Figure FDA0002254155300000036
式中，σk-1表示当前采样时刻的无迹变化参数，
Figure FDA0002254155300000037
表示当前采样时刻的谐波电流中间值；
所述当前采样时刻的无迹变化参数的计算式如下：
Figure FDA0002254155300000038
式中，Sxx,k-1表示表示当前采样时刻的谐波电流的协方差，λ为缩放比例系数。
7.如权利要求6所述的方法，其特征在于，所述基于时间更新的结果和预设遗忘因子，基于Sage-Husa方法对所述非线性状态方程与量测方程模型进行量测更新，包括：
基于时间更新的结果，对Sigma点进行重采样，得到Sigma点；
通过测量方程对Sigma点进行非线性变换并计算测量残差；
根据所述测量残差并考虑预设遗忘因子，基于Sage-Husa采用平方根无迹卡尔曼滤波方法对系统观测噪声协方差的平方根进行更新；
基于系统观测噪声协方差的平方根，对谐波电压和谐波电流间的协方差以及谐波电压的协方差进行更新；
根据谐波电压和谐波电流间的协方差以及谐波电压的协方差，计算卡尔曼滤波增益。
8.如权利要求7所述的方法，其特征在于，所述通过测量方程对重采样的Sigma点集进行非线性变换并计算测量残差，如下式所示：
Figure FDA0002254155300000041
式中，H为观测函数，χi,k表示重采样的下一采样时刻的谐波电流对应的第i个Sigma点的值，Uk-1表示当前采样时刻的谐波电压，
Figure FDA0002254155300000042
第i个Sigma点在计算均值时的权重，
Figure FDA0002254155300000043
表示下一采样时刻的添加权重后修正后且经量测变换后的谐波电压，Zk表示对应下一采样时刻采集的谐波电压观测值，ek表示下一采样时刻谐波电压的量测残差，Z* i,k表示对应
Figure FDA0002254155300000044
的第i个Sigma点的中间过程变量。
9.如权利要求8所述的方法，其特征在于，所述根据所述测量残差并考虑预设遗忘因子，基于Sage-Husa采用平方根无迹卡尔曼滤波方法对系统观测噪声协方差的平方根进行更新，如下式所示：
Figure FDA0002254155300000045
式中，
Figure FDA0002254155300000046
表示下一采样时刻系统观测噪声协方差的平方根，R** k表示对应
Figure FDA0002254155300000047
的第一中间过程变量，R* k表示对应
Figure FDA0002254155300000048
的第二中间过程变量，
Figure FDA0002254155300000049
表示当前采样时刻系统观测噪声协方差的平方根，ek表示下一采样时刻谐波电压的量测残差，dk为根据遗忘因子计算的下一采样时刻的遗忘修正因子，
Figure FDA00022541553000000410
表示第i个Sigma点在计算协方差时的权重，
Figure FDA00022541553000000411
表示下一采样时刻的添加权重后修正后第i个Sigma点且经量测变换后的谐波电压，
Figure FDA00022541553000000412
表示当前采样时刻的添加权重后修正后且经量测变换后第i个Sigma点的谐波电压，R* k T表示R* k的转置，Cholupdata为Cholesky分解运算，diag为构造对角矩阵运算；
下一采样时刻的遗忘修正因子dk的计算式如下：
dk＝(1-b)/(1-bk+1)
式中，b为遗忘因子。
10.如权利要求9所述的方法，其特征在于，所述基于系统观测噪声协方差的平方根，对谐波电压和谐波电流间的协方差以及谐波电压的协方差进行更新，如下式所示：
Figure FDA0002254155300000051
式中，Pxz,k表示下一采样时刻谐波电压和谐波电流间的协方差，Szz,k表示下一采样时刻谐波电压的协方差，S* zz,k表示对应Szz,k的中间过程变量，
Figure FDA0002254155300000052
表示下一采样时刻的添加权重后修正后且经量测变换后的谐波电压，Z* i,k表示对应
Figure FDA0002254155300000053
的第i个Sigma点的中间过程变量，
Figure FDA0002254155300000054
表示下一采样时刻得到的谐波电流中间值，χ* i,k表示对应χi,k的中间过程变量，χi,k表示下一采样时刻的谐波电流对应的第i个Sigma点的值，
Figure FDA0002254155300000055
表示第i个Sigma点在计算协方差时的权重，R* k表示对应
Figure FDA0002254155300000056
的第二中间过程变量，
Figure FDA0002254155300000057
表示第0个Sigma点在计算协方差时的权重，Z* 0,k-1表示对应
Figure FDA0002254155300000058
的第0个Sigma点的中间过程变量，
Figure FDA0002254155300000059
表示当前采样时刻的添加权重后修正后且经量测变换后的谐波电压，qr为正交三角矩阵分解运算，Cholupdata为Cholesky分解运算。
11.如权利要求10所述的方法，其特征在于，所述卡尔曼滤波增益的计算式如下：
Figure FDA00022541553000000510
式中，Kk表示下一采样时刻的卡尔曼滤波增益，ST zz,k表示Szz,k的转置。
12.如权利要求7所述的方法，其特征在于，所述基于量测更新的结果，对所述非线性状态方程与量测方程模型进行状态更新，得到谐波电流的估算值，包括：
根据所述非线性状态方程与量测方程模型中的谐波电流中间值、测量残差和卡尔曼滤波增益进行计算，得到谐波电流的估算值。
13.如权利要求12所述的方法，其特征在于，所述谐波电流的估算值的计算式如下：
Figure FDA00022541553000000511
其中，Xk表示下一采样时刻的谐波电流的估算值，
Figure FDA00022541553000000512
表示下一采样时刻得到的谐波电流中间值，Kk表示下一采样时刻的卡尔曼滤波增益，ek表示下一采样时刻谐波电压的量测残差。
14.如权利要求12所述的方法，其特征在于，所述得到谐波电流的估算值之后，还包括：
对所述谐波电流的协方差估计值进行更新，以及对系统过程噪声协方差向量的平方根进行更新；
将当前采样时刻非线性状态方程与量测方程模型的参数传递至下一采样时刻的非线性状态方程与量测方程模型。
15.如权利要求14所述的方法，其特征在于，所述谐波电流的协方差估计值的计算式如下：
Sk＝cholupdate{Sk-1,KkSzz,k,-1}
其中，Sk表示下一采样时刻谐波电流的协方差估计值，Sk-1表示当前采样时刻的谐波电流的协方差估计值，Kk表示下一采样时刻的卡尔曼滤波增益，Szz,k表示下一采样时刻谐波电压的协方差，Cholupdata为Cholesky分解运算；
所述系统过程噪声协方差向量的平方根的计算式如下：
Figure FDA0002254155300000061
其中，
Figure FDA0002254155300000062
表示下一采样时刻的系统过程噪声协方差向量的平方根，Qk-1表示第k-1次递推后的系统过程噪声协方差向量，dk为下一采样时刻的遗忘修正因子，Kk表示下一采样时刻的卡尔曼滤波增益，ek表示下一采样时刻谐波电压的量测残差，Szz,k表示第k次递推后谐波电压的协方差，Q** k表示对应
Figure FDA0002254155300000063
的第一中间过程变量，Q* k表示对应
Figure FDA0002254155300000064
的第二中间过程变量，Cholupdata为Cholesky分解运算，diag为构造对角矩阵运算。
16.如权利要求14所述的方法，其特征在于，所述将当前采样时刻非线性状态方程与量测方程模型的参数传递至下一采样时刻的非线性状态方程与量测方程模型，包括：
将当前采样时刻的谐波电流中间值作为下一采样时刻非线性状态方程与量测方程模型的谐波电流中间值，将当前采样时刻的谐波电流的协方差估计值作为下一采样时刻非线性状态方程与量测方程模型的谐波电流的协方差估计值，将当前采样时刻的系统过程噪声协方差向量的平方根作为下一采样时刻非线性状态方程与量测方程模型的谐波电流的协方差估计值，将当前采样时刻的系统观测噪声协方差的平方根作为下一采样时刻非线性状态方程与量测方程模型的系统观测噪声协方差的平方根。
17.一种改进动态谐波估计系统，其特征在于，包括：数据采集模块、数据输入模块和谐波估计模块；
所述数据采集模块，用于采集电网各支路或母线的同步相量测量装置的谐波电压；
所述数据输入模块，用于将所述谐波电压带入预先建立的非线性状态方程与量测方程模型；
所述谐波估计模块，用于采用改进Sage-Husa平方根无迹卡尔曼滤波算法，根据所述非线性状态方程与量测方程模型，对待估测支路进行谐波电流估计计算，得到待估测支路谐波电流的估算值；
其中，所述改进Sage-Husa平方根无迹卡尔曼滤波算法包括：基于量测更新和遗忘因子，对所述非线性状态方程与量测方程模型进行状态更新。
18.如权利要求17所述的系统，其特征在于，所述谐波估计模块包括：时间更新模块、量测更新模块和状态更新模块；
所述时间更新模块，用于基于所述非线性状态方程与量测方程模型，采用无迹卡尔曼滤波算法进行时间更新；所述非线性状态方程与量测方程模型中的量测量序列的初始值由初次采样时刻采集的谐波电压初始化得到；
所述量测更新模块，用于基于时间更新的结果和Sage-Husa方法对所述非线性状态方程与量测方程模型进行量测更新；
所述状态更新模块，用于基于量测更新的结果和预设遗忘因子，对所述非线性状态方程与量测方程模型进行状态更新，得到谐波电流的估算值。
CN201911045986.3A 2019-10-30 2019-10-30 一种改进动态谐波估计方法和系统 Pending CN110907702A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911045986.3A CN110907702A (zh) 2019-10-30 2019-10-30 一种改进动态谐波估计方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911045986.3A CN110907702A (zh) 2019-10-30 2019-10-30 一种改进动态谐波估计方法和系统

Publications (1)

Publication Number Publication Date
CN110907702A true CN110907702A (zh) 2020-03-24

Family

ID=69816007

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911045986.3A Pending CN110907702A (zh) 2019-10-30 2019-10-30 一种改进动态谐波估计方法和系统

Country Status (1)

Country Link
CN (1) CN110907702A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111551785A (zh) * 2020-04-29 2020-08-18 南京理工大学 基于无迹卡尔曼滤波的频率与谐波检测方法
CN114236236A (zh) * 2021-12-17 2022-03-25 福州大学 一种基于区间动态状态估计的谐波源定位方法
CN114545162A (zh) * 2022-02-15 2022-05-27 国网湖南省电力有限公司 一种分布式输电线路绝缘性能在线监测方法
CN115032458A (zh) * 2022-06-02 2022-09-09 北京妙微科技有限公司 一种基于电网系统的谐波阻抗估计方法及计算机存储介质

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111551785A (zh) * 2020-04-29 2020-08-18 南京理工大学 基于无迹卡尔曼滤波的频率与谐波检测方法
CN111551785B (zh) * 2020-04-29 2022-08-12 南京理工大学 基于无迹卡尔曼滤波的频率与谐波检测方法
CN114236236A (zh) * 2021-12-17 2022-03-25 福州大学 一种基于区间动态状态估计的谐波源定位方法
CN114236236B (zh) * 2021-12-17 2024-02-06 福州大学 一种基于区间动态状态估计的谐波源定位方法
CN114545162A (zh) * 2022-02-15 2022-05-27 国网湖南省电力有限公司 一种分布式输电线路绝缘性能在线监测方法
CN115032458A (zh) * 2022-06-02 2022-09-09 北京妙微科技有限公司 一种基于电网系统的谐波阻抗估计方法及计算机存储介质
CN115032458B (zh) * 2022-06-02 2023-08-11 北京妙微科技有限公司 一种基于电网系统的谐波阻抗估计方法及计算机存储介质

Similar Documents

Publication Publication Date Title
CN110907702A (zh) 一种改进动态谐波估计方法和系统
Debs et al. A dynamic estimator for tracking the state of a power system
CN107590317B (zh) 一种计及模型参数不确定性的发电机动态估计方法
Schweppe et al. Power system static-state estimation, Part I: Exact model
CN108155648B (zh) 基于自适应h无穷扩展卡尔曼滤波的状态估计方法
CN114124033A (zh) 卡尔曼滤波器的实现方法、装置、存储介质和设备
Nakamori et al. Recursive estimators of signals from measurements with stochastic delays using covariance information
CN112713587A (zh) 一种基于平方根容积卡尔曼滤波器的配电网动态状态估计方法及系统
CN109754013B (zh) 一种基于无迹卡尔曼滤波的电力系统混合量测融合方法
CN110490366A (zh) 基于变分模态分解和迭代决策树的径流量预测方法
CN111445108A (zh) 数据驱动的配电网线变关系诊断方法、装置及系统
CN113608121A (zh) 基于模糊分数阶无迹卡尔曼滤波的锂电池soc估计方法
CN114567288B (zh) 基于变分贝叶斯的分布协同非线性系统状态估计方法
Kulikov et al. Accurate state estimation of stiff continuous-time stochastic models in chemical and other engineering
CN110929835B (zh) 一种新型碳化硅基航空功率变换器故障诊断方法及系统
CN115130770A (zh) 基于张量特征重构的工业废水排放水质预测方法
Sadhu et al. Central difference formulation of risk-sensitive filter
CN112381279B (zh) 一种基于vmd和bls组合模型的风电功率预测方法
CN114239796A (zh) 一种基于扩展卡尔曼滤波的电力系统状态估计方法
CN107204616B (zh) 基于自适应稀疏伪谱法的电力系统随机状态估计方法
CN106100609B (zh) 单状态变量和两级Kalman滤波器时间尺度算法
CN113484818B (zh) 基于滑动窗口的抗高频采集异常电能表精准定位方法
Mercieca et al. Unscented Rauch-Tung-Striebel smoothing for nonlinear descriptor systems
Nerger Parallel filter algorithms for data assimilation in oceanography
CN113884914A (zh) 一种渐消并行卡尔曼滤波动力电池荷电状态估计方法

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