CN111597688A - 提取时变工频和次同步频率分量的方法及系统 - Google Patents

提取时变工频和次同步频率分量的方法及系统 Download PDF

Info

Publication number
CN111597688A
CN111597688A CN202010335054.9A CN202010335054A CN111597688A CN 111597688 A CN111597688 A CN 111597688A CN 202010335054 A CN202010335054 A CN 202010335054A CN 111597688 A CN111597688 A CN 111597688A
Authority
CN
China
Prior art keywords
subsynchronous
fundamental
state vector
frequency
error covariance
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
CN202010335054.9A
Other languages
English (en)
Other versions
CN111597688B (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.)
Yunnan Power Grid Co Ltd
Original Assignee
Tsinghua University
Yunnan Power Grid 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 Tsinghua University, Yunnan Power Grid Co Ltd filed Critical Tsinghua University
Priority to CN202010335054.9A priority Critical patent/CN111597688B/zh
Publication of CN111597688A publication Critical patent/CN111597688A/zh
Application granted granted Critical
Publication of CN111597688B publication Critical patent/CN111597688B/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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • 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
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/06—Wind turbines or wind farms

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Supply And Distribution Of Alternating Current (AREA)

Abstract

本发明涉及一种提取时变工频和次同步频率分量的方法及系统，属于电力系统动态监测技术领域。本发明包括以下内容：（1）构建量测模型；（2）状态向量初始化；（3）更新状态向量和误差协方差；（4）计算卡尔曼增益；（5）利用卡尔曼增益更新状态和误差协方差；（6）更新基波频率；（7）计算次同步分量。本发明方法能够适用于基波和次同步谐波具有时变特性的情况，可以准确地提取出次同步谐波，具有很好的有效性和实用性；同时，本发明在模拟测试信号以及存在次同步振荡现象的真实风电系统上得到了验证。

Description

提取时变工频和次同步频率分量的方法及系统
技术领域
本发明属于电力系统动态监测技术领域，具体涉及一种提取时变工频和次同步频率分量的方法及系统。
背景技术
风电机组控制器和串补电网相互作用会引起次同步控制相互作用(subsynchronous control interaction,SSCI)，与传统火电机组和串补作用产生的次同步振荡不同，其振荡频率随系统运行状态变化而变化(如风机并网数量、风速大小、风机控制结构及参数、串补度等)，且变化范围较广，具有很强的不确定性。SSCI在风电系统中主要表现为次/超同步电流、电压和谐波功率的持续增长(发散)或等幅振荡，严重危及大型风电系统可靠运行。
现有抑制次同步控制相互作用的方法往往需要从测量得到的信号中精确提取出次同步分量，但信号基波及次同步谐波频率具有时变特性，因此需要提出一种方法能够从监测数据中提取基波和次同步谐波频率的时变特性。
现有提取次同步谐波频率的方法，比如线性卡尔曼滤波器，对频率固定的基波和谐波分量具有良好的作用，但不能用于具有时变特性的基波和次同步谐波的分析。其原因在于现有方法不能追踪具有时变特征的基波和次同步分量的频率。因此如何克服现有技术的不足是目前电力系统动态监测技术领域亟需解决的问题。
发明内容
本发明的目的是为了解决现有技术的不足，针对风电系统发生的次同步控制相互作用，提供一种提取时变工频和次同步频率分量的方法及系统，该方法能够适用于基波和次同步谐波具有时变特性的情况，可以准确地提取出次同步谐波，具有很好的有效性和实用性。
为实现上述目的，本发明采用的技术方案如下：
本发明第一方面提供提取时变工频和次同步频率分量的方法，包含有基波和次同步频率分量的监测信号表示为：
z(k)＝X(k)+vss(k)；
式中，z(k)为测量信号，vss(k)表示次同步分量；X(k)为基波分量，k为采样时间点，k≥1的整数；
具体包括如下步骤：
步骤(1)，首先构建量测模型：
将基波分量X(k)表示为状态向量的形式，即：
X(k)＝[x1(k),x2(k)]*；
式中，x1为基波相量实部，x2为基波相量虚部；
x1(k)＝x1(k-1)cos(2πf0(k-1)T)+x2(k-1)sin(2πf0(k-1)T)；
x2(k)＝-x1(k-1)sin(2πf0(k-1)T)+x2(k-1)sin(2πf0(k-1)T)；
式中，T表示采样周期，上标*表示矩阵转置；f0表示基波频率；
基波分量X重新写为：
X(k)＝FX(k-1)
式中，F为转换矩阵，表示为：
Figure BDA0002466268560000021
因此，量测模型表示为
z(k)＝HX(k-1)+vss(k)；
式中，H表示测量矩阵H＝[1,0]，是将估计量转化为量测量；
步骤(2)，状态向量初始化：
状态向量的初始值：
Figure BDA0002466268560000022
误差协方差的初始值：P(0)＝10I2×2；
基波频率的初始值：f0(0)＝50Hz；
其中，I表示单位对角矩阵；
步骤(3)，更新状态向量和误差协方差
更新状态向量：
Figure BDA0002466268560000023
更新误差协方差：P-(k)＝FP(k-1)F*+Q；
式中，Q为噪声协方差矩阵，上标—表示在不利用该时刻的测量值的情况下计算状态向量和过程协方差矩阵；
步骤(4)，计算卡尔曼增益：
卡尔曼增益表示为：
Figure BDA0002466268560000031
式中，R＝0.1；
步骤(5)，计算状态向量和误差协方差：
利用卡尔曼增益和测量值更新状态和误差协方差：
状态向量计算：
Figure BDA0002466268560000032
误差协方差计算：P(k)＝P-(k)+K(k)HP-(k)；
其中，z(k)为测量信号；
步骤(6)：更新基波频率：
由步骤(5)中计算得到
Figure BDA0002466268560000033
后，将其表示为复数形式，即
Figure BDA0002466268560000034
因此，从相角变化能得到基波频率变化，其计算公式为：
相角公式为：
Figure BDA0002466268560000035
基波频率更新为：
Figure BDA0002466268560000036
步骤(7)：计算次同步分量：
根据步骤(6)中计算的
Figure BDA0002466268560000037
可计算残差得到次同步谐波分量为：
Figure BDA0002466268560000038
进一步，优选的是，噪声协方差矩阵Q为对角矩阵
Figure BDA0002466268560000039
Figure BDA0002466268560000041
通过最大化密度函数得到：
Figure BDA0002466268560000042
式中，函数f定义为
Figure BDA0002466268560000043
次同步谐波vss(k-1)是上一时刻通过残差计算得到的数值；
Figure BDA0002466268560000044
本发明第二方面提供提取时变工频和次同步频率分量的系统，包括：
量测模型构建模块，用于将基波分量X(k)表示为状态向量的形式后，将基波分量重新写为X(k)＝FX(k-1)，从而得到量测模型z(k)＝HX(k-1)+vss(k)；
其中，
Figure BDA0002466268560000045
H＝[1,0]，vss(k)表示次同步分量，k为采样时间点；
状态向量初始化模块，用于对状态向量、误差协方差、基波频率进行初始化；
第一处理模块，用于更新状态向量和误差协方差；
第二处理模块，用于计算卡尔曼增益；
第三处理模块，用于计算状态向量和误差协方差；
基波频率与次同步谐波分量提取模块，用于提取基波频率与次同步谐波分量并进行显示。
本发明第三方面提供一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如上述提取时变工频和次同步频率分量的方法的步骤。
本发明第四方面提供一种非暂态计算机可读存储介质，其上存储有计算机程序，其特征在于，该计算机程序被处理器执行时实现如上述提取时变工频和次同步频率分量的方法的步骤。
本发明既可以适用于连续信号，同时也可以适用于离散信号；
本发明方法具体实现过程中，可以参考卡尔曼滤波的硬件电路实现方式，依照本发明的方法加以改造。
本发明与现有技术相比，其有益效果为：
1、本发明方法及系统适用于次同步振荡频率时变的情况，也可以适用于一般情况，适用范围广，并且能够准确提取出基波及次同步谐波频率时变特性；
2、本发明方法及系统具有很好的实用性和有效性，已经在实际风电场中通过验证；
3、本发明方法及系统使用自适应协方差矩阵，具有很强的追踪性能。
附图说明
图1为本发明的算法流程图；
图2为某典型系统拓扑结构图；
图3为次同步振荡现象时域波形和FFT结果；其中，(a)为电流信号时域曲线；(b)为；FFT结果
图4为基波频率和次同步相量的估算结果；其中，(a)为基波频率变化曲线；(b)为次同步分量幅值；
图5是本发明提取时变工频和次同步频率分量的系统的结构示意图；
图6为本发明电子设备结构示意图。
具体实施方式
下面结合实施例对本发明作进一步的详细描述。
本领域技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限定本发明的范围。实施例中未注明具体技术或条件者，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。所用材料或设备未注明生产厂商者，均为可以通过购买获得的常规产品。
提取时变工频和次同步频率分量的方法，包括如下步骤：
包含有基波和次同步频率分量的监测信号表示为：
z(k)＝X(k)+vss(k)；
式中，z(k)为测量信号，vss(k)表示次同步分量；X(k)为基波分量，k为采样时间点，k≥1的整数；
步骤(1)，首先构建量测模型：
将基波分量X(k)表示为状态向量的形式，即
X(k)＝[x1(k),x2(k)]*；
式中，x1为基波相量实部，x2为基波相量虚部；
x1(k)＝x1(k-1)cos(2πf0(k-1)T)+x2(k-1)sin(2πf0(k-1)T)；
x2(k)＝-x1(k-1)sin(2πf0(k-1)T)+x2(k-1)sin(2πf0(k-1)T)；
式中，T表示采样周期，上标*表示矩阵转置；f0表示计算的基波频率。
基波分量X重新写为：
X(k)＝FX(k-1)
式中，F为转换矩阵，表示为：
Figure BDA0002466268560000061
因此，量测模型表示为
z(k)＝HX(k-1)+vss(k)；
式中，H表示测量矩阵H＝[1,0]，是将估计量转化为量测量；
步骤(2)，状态向量初始化：
状态向量的初始值：
Figure BDA0002466268560000062
误差协方差的初始值：P0＝10I2×2；
基波频率的初始值：f0(0)＝50Hz；
其中，I表示单位对角矩阵；
步骤(3)，更新状态向量和误差协方差
更新状态向量：
Figure BDA0002466268560000063
更新误差协方差：P-(k)＝FP(k-1)F*+Q；
式中，Q为噪声协方差矩阵，上标—表示在不利用该时刻的测量值的情况下计算状态向量和过程协方差矩阵。
步骤(4)，计算卡尔曼增益
卡尔曼增益表示为：
Figure BDA0002466268560000071
式中，R＝0.1；
步骤(5)，计算状态向量和误差协方差：
利用卡尔曼增益和测量值更新状态和误差协方差：
状态向量计算：
Figure BDA0002466268560000072
误差协方差计算：P(k)＝P-(k)+K(k)HP-(k)；
其中，z(k)为测量信号；
步骤(6)：更新基波频率
由步骤(5)中计算得到
Figure BDA0002466268560000073
后，将其表示为复数形式，即
Figure BDA0002466268560000074
因此，从相角变化能得到基波频率变化，其计算公式为：
相角公式为：
Figure BDA0002466268560000075
基波频率更新为：
Figure BDA0002466268560000076
步骤(7)：计算次同步分量
根据步骤(6)中计算的
Figure BDA0002466268560000077
可计算残差得到次同步谐波分量为
Figure BDA0002466268560000078
优选，噪声协方差矩阵Q为对角矩阵；
Figure BDA0002466268560000079
Figure BDA0002466268560000081
通过最大化密度函数得到：
Figure BDA0002466268560000082
式中，函数f定义为
Figure BDA0002466268560000083
次同步谐波vss(k-1)是上一时刻通过残差计算得到的数值；
Figure BDA0002466268560000084
重复进行式(3)～(7)可计算得到随时间的变化的基波频率和次同步频率分量。
本方法流程图见图1。
如图5所示，提取时变工频和次同步频率分量的系统，包括：
量测模型构建模块101，用于将基波分量X(k)表示为状态向量的形式后，将基波分量重新写为X(k)＝FX(k-1)，从而得到量测模型z(k)＝HX(k-1)+vss(k)；
其中，
Figure BDA0002466268560000085
H＝[1,0]，vss(k)表示次同步分量，k为采样时间点；
状态向量初始化模块102，用于对状态向量、误差协方差、基波频率进行初始化；
第一处理模块103，用于更新状态向量和误差协方差；
第二处理模块104，用于计算卡尔曼增益；
第三处理模块105，用于计算状态向量和误差协方差；
基波频率与次同步谐波分量提取模块106，用于提取基波频率与次同步谐波分量并进行显示。
本系统6个模块的具体计算和提取方法如提取时变工频和次同步频率分量的方法中所述。
本发明实施例提供的提取时变工频和次同步频率分量的系统，该系统能够适用于基波和次同步谐波具有时变特性的情况，可以准确地提取出次同步谐波，具有很好的有效性和实用性；同时，本发明在模拟测试信号以及存在次同步振荡现象的真实风电系统上得到了验证。
本发明实施例提供的系统是用于执行上述各方法实施例的，具体流程和详细内容请参照上述实施例，此处不再赘述。
图6为本发明实施例提供的电子设备结构示意图，参照图6，该电子设备可以包括：处理器(processor)201、通信接口(CommunicationsInterface)202、存储器(memory)203和通信总线204，其中，处理器201，通信接口202，存储器203通过通信总线204完成相互间的通信。处理器201可以调用存储器203中的逻辑指令，以执行如下方法：
z(k)＝X(k)+vss(k)；
式中，z(k)为测量信号，vss(k)表示次同步分量；X(k)为基波分量，k为采样时间点，k≥1的整数；
步骤(1)，首先构建量测模型：
将基波分量X(k)表示为状态向量的形式，即：
X(k)＝[x1(k),x2(k)]*；
式中，x1为基波相量实部，x2为基波相量虚部；
x1(k)＝x1(k-1)cos(2πf0(k-1)T)+x2(k-1)sin(2πf0(k-1)T)；
x2(k)＝-x1(k-1)sin(2πf0(k-1)T)+x2(k-1)sin(2πf0(k-1)T)；
式中，T表示采样周期，上标*表示矩阵转置；f0表示基波频率；
基波分量X重新写为：
X(k)＝FX(k-1)
式中，F为转换矩阵，表示为：
Figure BDA0002466268560000091
因此，量测模型表示为
z(k)＝HX(k-1)+vss(k)；
式中，H表示测量矩阵H＝[1,0]，是将估计量转化为量测量；
步骤(2)，状态向量初始化：
状态向量的初始值：
Figure BDA0002466268560000101
误差协方差的初始值：P(0)＝10I2×2；
基波频率的初始值：f0(0)＝50Hz；
其中，I表示单位对角矩阵；
步骤(3)，更新状态向量和误差协方差
更新状态向量：
Figure BDA0002466268560000102
更新误差协方差：P-(k)＝FP(k-1)F*+Q；
式中，Q为噪声协方差矩阵，上标—表示在不利用该时刻的测量值的情况下计算状态向量和过程协方差矩阵；
步骤(4)，计算卡尔曼增益：
卡尔曼增益表示为：
Figure BDA0002466268560000103
式中，R＝0.1；
步骤(5)，计算状态向量和误差协方差：
利用卡尔曼增益和测量值更新状态和误差协方差：
状态向量计算：
Figure BDA0002466268560000104
误差协方差计算：P(k)＝P-(k)+K(k)HP-(k)；
其中，z(k)为测量信号；
步骤(6)：更新基波频率：
由步骤(5)中计算得到
Figure BDA0002466268560000105
后，将其表示为复数形式，即
Figure BDA0002466268560000106
因此，从相角变化能得到基波频率变化，其计算公式为：
相角公式为：
Figure BDA0002466268560000107
基波频率更新为：
Figure BDA0002466268560000108
步骤(7)：计算次同步分量：
根据步骤(6)中计算的
Figure BDA0002466268560000111
可计算残差得到次同步谐波分量为：
Figure BDA0002466268560000112
此外，上述的存储器203中的逻辑指令可以通过软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
另一方面，本发明实施例还提供一种非暂态计算机可读存储介质，其上存储有计算机程序，该计算机程序被处理器执行时实现以执行上述各实施例提供的提取时变工频和次同步频率分量的方法，例如包括：z(k)＝X(k)+vss(k)；
式中，z(k)为测量信号，vss(k)表示次同步分量；X(k)为基波分量，k为采样时间点，k≥1的整数；
步骤(1)，首先构建量测模型：
将基波分量X(k)表示为状态向量的形式，即：
X(k)＝[x1(k),x2(k)]*；
式中，x1为基波相量实部，x2为基波相量虚部；
x1(k)＝x1(k-1)cos(2πf0(k-1)T)+x2(k-1)sin(2πf0(k-1)T)；
x2(k)＝-x1(k-1)sin(2πf0(k-1)T)+x2(k-1)sin(2πf0(k-1)T)；
式中，T表示采样周期，上标*表示矩阵转置；f0表示基波频率；
基波分量X重新写为：
X(k)＝FX(k-1)
式中，F为转换矩阵，表示为：
Figure BDA0002466268560000121
因此，量测模型表示为
z(k)＝HX(k-1)+vss(k)；
式中，H表示测量矩阵H＝[1,0]，是将估计量转化为量测量；
步骤(2)，状态向量初始化：
状态向量的初始值：
Figure BDA0002466268560000122
误差协方差的初始值：P(0)＝10I2×2；
基波频率的初始值：f0(0)＝50Hz；
其中，I表示单位对角矩阵；
步骤(3)，更新状态向量和误差协方差
更新状态向量：
Figure BDA0002466268560000123
更新误差协方差：P-(k)＝FP(k-1)F*+Q；
式中，Q为噪声协方差矩阵，上标—表示在不利用该时刻的测量值的情况下计算状态向量和过程协方差矩阵；
步骤(4)，计算卡尔曼增益：
卡尔曼增益表示为：
Figure BDA0002466268560000124
式中，R＝0.1；
步骤(5)，计算状态向量和误差协方差：
利用卡尔曼增益和测量值更新状态和误差协方差：
状态向量计算：
Figure BDA0002466268560000125
误差协方差计算：P(k)＝P-(k)+K(k)HP-(k)；
其中，z(k)为测量信号；
步骤(6)：更新基波频率：
由步骤(5)中计算得到
Figure BDA0002466268560000131
后，将其表示为复数形式，即
Figure BDA0002466268560000132
因此，从相角变化能得到基波频率变化，其计算公式为：
相角公式为：
Figure BDA0002466268560000133
基波频率更新为：
Figure BDA0002466268560000134
步骤(7)：计算次同步分量：
根据步骤(6)中计算的
Figure BDA0002466268560000135
可计算残差得到次同步谐波分量为：
Figure BDA0002466268560000136
以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性的劳动的情况下，即可以理解并实施。
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到各实施方式可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件。基于这样的理解，上述技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在计算机可读存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行各个实施例或者实施例的某些部分所述的方法。
应用实例
将本发明方法应用于某实际风电场。该风电系统考虑两个风电场，经过220kV变电站升压并通过串联补偿线路输出到大电网。仿真3s后投入C2产生次同步振荡。由于不同的系统参数，比如风速、投入风机数量，次同步振荡的频率和幅值不同。该系统拓扑结构图见图2，系统参数表如表1所示，时域仿真结果及频率分析结果如图3、4所示。
表1系统参数表
系统频率 50Hz
风机侧变压器TF<sub>w</sub> 0.69/220kV
变电站变压器TF<sub>g</sub> 220/500kV
风电场1电缆长度L<sub>w1</sub> 48km
风电场2电缆长度L<sub>w2</sub> 65km
并联传输线长度L<sub>1</sub>、L<sub>2</sub>: 220km
串联补偿电容C<sub>1</sub>、C<sub>2</sub> 90μF
风电场1风速 5.0m/s
风电场2风速 4.9m/s
从图3可以看出，3s后由幅值逐渐增大的次同步振荡出现，快速傅里叶变换(FFT)结果显示有50Hz的基波以及10Hz的次同步振荡，图4通过周期验证，表明本发明能够很好地提取次同步分量。
以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (5)

1.提取时变工频和次同步频率分量的方法，其特征在于，包含有基波和次同步频率分量的监测信号表示为：
z(k)＝X(k)+vss(k)；
式中，z(k)为测量信号，vss(k)表示次同步分量；X(k)为基波分量，k为采样时间点，k≥1的整数；
具体包括如下步骤：
步骤(1)，首先构建量测模型：
将基波分量X(k)表示为状态向量的形式，即：
X(k)＝[x1(k),x2(k)]*；
式中，x1为基波相量实部，x2为基波相量虚部；
x1(k)＝x1(k-1)cos(2πf0(k-1)T)+x2(k-1)sin(2πf0(k-1)T)；
x2(k)＝-x1(k-1)sin(2πf0(k-1)T)+x2(k-1)sin(2πf0(k-1)T)；
式中，T表示采样周期，上标*表示矩阵转置；f0表示基波频率；
基波分量X重新写为：
X(k)＝FX(k-1)
式中，F为转换矩阵，表示为：
Figure FDA0002466268550000011
因此，量测模型表示为
z(k)＝HX(k-1)+vss(k)；
式中，H表示测量矩阵H＝[1,0]，是将估计量转化为量测量；
步骤(2)，状态向量初始化：
状态向量的初始值：
Figure FDA0002466268550000012
误差协方差的初始值：P(0)＝10I2×2；
基波频率的初始值：f0(0)＝50Hz；
其中，I表示单位对角矩阵；
步骤(3)，更新状态向量和误差协方差
更新状态向量：
Figure FDA0002466268550000021
更新误差协方差：P-(k)＝FP(k-1)F*+Q；
式中，Q为噪声协方差矩阵，上标—表示在不利用该时刻的测量值的情况下计算状态向量和过程协方差矩阵；
步骤(4)，计算卡尔曼增益：
卡尔曼增益表示为：
Figure FDA0002466268550000022
式中，R＝0.1；
步骤(5)，计算状态向量和误差协方差：
利用卡尔曼增益和测量值更新状态和误差协方差：
状态向量计算：
Figure FDA0002466268550000023
误差协方差计算：P(k)＝P-(k)+K(k)HP-(k)；
其中，z(k)为测量信号；
步骤(6)：更新基波频率：
由步骤(5)中计算得到
Figure FDA0002466268550000024
后，将其表示为复数形式，即
Figure FDA0002466268550000025
因此，从相角变化能得到基波频率变化，其计算公式为：
相角公式为：
Figure FDA0002466268550000026
基波频率更新为：
Figure FDA0002466268550000027
步骤(7)：计算次同步分量：
根据步骤(6)中计算的
Figure FDA0002466268550000028
可计算残差得到次同步谐波分量为：
Figure FDA0002466268550000029
2.根据权利要求1所述的提取时变工频和次同步频率分量的方法，其特征在于，噪声协方差矩阵Q为对角矩阵；
Figure FDA0002466268550000031
Figure FDA0002466268550000032
通过最大化密度函数得到：
Figure FDA0002466268550000033
式中，函数f定义为
Figure FDA0002466268550000034
次同步谐波vss(k-1)是上一时刻通过残差计算得到的数值；
Figure FDA0002466268550000035
3.提取时变工频和次同步频率分量的系统，其特征在于，包括：
量测模型构建模块，用于将基波分量X(k)表示为状态向量的形式后，将基波分量重新写为X(k)＝FX(k-1)，从而得到量测模型z(k)＝HX(k-1)+vss(k)；
其中，
Figure FDA0002466268550000036
H＝[1,0]，vss(k)表示次同步分量，k为采样时间点；
状态向量初始化模块，用于对状态向量、误差协方差、基波频率进行初始化；
第一处理模块，用于更新状态向量和误差协方差；
第二处理模块，用于计算卡尔曼增益；
第三处理模块，用于计算状态向量和误差协方差；
基波频率与次同步谐波分量提取模块，用于提取基波频率与次同步谐波分量并进行显示。
4.一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如权利要求1至2任一项所述提取时变工频和次同步频率分量的方法的步骤。
5.一种非暂态计算机可读存储介质，其上存储有计算机程序，其特征在于，该计算机程序被处理器执行时实现如权利要求1至2任一项所述提取时变工频和次同步频率分量的方法的步骤。
CN202010335054.9A 2020-04-24 2020-04-24 提取时变工频和次同步频率分量的方法及系统 Active CN111597688B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010335054.9A CN111597688B (zh) 2020-04-24 2020-04-24 提取时变工频和次同步频率分量的方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010335054.9A CN111597688B (zh) 2020-04-24 2020-04-24 提取时变工频和次同步频率分量的方法及系统

Publications (2)

Publication Number Publication Date
CN111597688A true CN111597688A (zh) 2020-08-28
CN111597688B CN111597688B (zh) 2022-11-15

Family

ID=72190687

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010335054.9A Active CN111597688B (zh) 2020-04-24 2020-04-24 提取时变工频和次同步频率分量的方法及系统

Country Status (1)

Country Link
CN (1) CN111597688B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105223418A (zh) * 2015-09-22 2016-01-06 清华大学 次同步和超同步谐波相量的测量方法及测量装置
CN106786561A (zh) * 2017-02-20 2017-05-31 河海大学 一种基于自适应卡尔曼滤波的低频振荡模态参数辨识方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105223418A (zh) * 2015-09-22 2016-01-06 清华大学 次同步和超同步谐波相量的测量方法及测量装置
CN106786561A (zh) * 2017-02-20 2017-05-31 河海大学 一种基于自适应卡尔曼滤波的低频振荡模态参数辨识方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
李金等: "风电并网系统次/超同步振荡的动态监测方法研究", 《现代电力》 *
王茂海等: "电力系统次同步振荡分量的快速在线检测算法", 《电力系统自动化》 *

Also Published As

Publication number Publication date
CN111597688B (zh) 2022-11-15

Similar Documents

Publication Publication Date Title
Paternina et al. Identification of electromechanical oscillatory modes based on variational mode decomposition
CN106845010B (zh) 基于改进SVD降噪和Prony的低频振荡主导模式辨识方法
CN103487652B (zh) 一种频率自适应实时分次谐波检测方法
CN102707122A (zh) 基于箕舌线的变步长lms谐波电流检测方法
CN104833852A (zh) 一种基于人工神经网络的电力系统谐波信号估计测量方法
CN110867889B (zh) 风电场/机组接入交流电网的振荡稳定性判别方法和系统
CN112018784A (zh) 一种基于同步相量测量数据的次同步谐振溯源方法
CN111597688B (zh) 提取时变工频和次同步频率分量的方法及系统
CN109684937A (zh) 一种基于fft及数学形态法的信号去噪方法及装置
Gurusinghe et al. Implementation of smart DFT-based PMU model in the real-time digital simulator
CN113794198B (zh) 抑制宽频振荡的方法、装置、终端及存储介质
CN116316909A (zh) 基于armax模型的电力系统惯量在线辨识方法与系统
CN110647720A (zh) 一种嵌入式平台下非平稳信号电能计量的方法
CN112345826B (zh) 孤网失稳状态下频率和暂态谐波测量方法
CN112688325B (zh) 基于二阶段改进itd算法的风电场次同步振荡监测方法
Farrokhifard et al. Different approaches for estimation of dampings and frequencies of electromechanical modes from PMU ambient data
CN115219787A (zh) 基于改进矩阵束的电网相量移动测量方法、系统及介质
CN110890753B (zh) 一种基于vmd算法的发电机组扰动源定位方法
CN109390957B (zh) 一种风电功率波动诱发系统强迫振荡的检测方法
Siavashi et al. Detection of voltage sag using unscented Kalman smoother
Mondal et al. Modified Taylor-levenberg ADALINE using adaptive learning parameter for power system harmonic estimation
Peng et al. Detection of lightly damped inter-area power oscillations using extended complex kalman filter
Guo et al. Modeling and simulation of power grid voltage harmonic detection method based on adaptive Kalman filter
Chen et al. Online SSO stability analysis-based oscillation parameter estimation in converter-tied grids
Yi et al. An anti mode mixing EMD algorithm for detecting the characteristics of low frequency oscillations in power system

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
TA01 Transfer of patent application right

Effective date of registration: 20220718

Address after: No.73 Tuodong Road, Guandu District, Kunming, Yunnan 650011

Applicant after: YUNNAN POWER GRID Co.,Ltd.

Address before: No.73 Tuodong Road, Guandu District, Kunming, Yunnan 650011

Applicant before: YUNNAN POWER GRID Co.,Ltd.

Applicant before: Tsinghua University

TA01 Transfer of patent application right
GR01 Patent grant
GR01 Patent grant