CN107271768A - 一种最小二乘拟合动态频率测量方法 - Google Patents
一种最小二乘拟合动态频率测量方法 Download PDFInfo
- Publication number
- CN107271768A CN107271768A CN201710383645.1A CN201710383645A CN107271768A CN 107271768 A CN107271768 A CN 107271768A CN 201710383645 A CN201710383645 A CN 201710383645A CN 107271768 A CN107271768 A CN 107271768A
- Authority
- CN
- China
- Prior art keywords
- mrow
- frequency
- msub
- msup
- delta
- 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
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
-
- 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Theoretical Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Algebra (AREA)
- Operations Research (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Probability & Statistics with Applications (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Life Sciences & Earth Sciences (AREA)
- Evolutionary Biology (AREA)
- Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
- Measuring Frequencies, Analyzing Spectra (AREA)
Abstract
本发明公开一种最小二乘拟合动态频率测量方法，依次包括以下步骤：将电网连续信号通过A/D转换采样，形成等时间间隔离散信号；采用FIR数字滤波器提取较为纯净的基波离散信号并形成采样矩阵；应用二元函数泰勒展开式，提取电网信号中的频率偏差量和频率变化率偏差量，离线计算常系数矩阵；应用最小二乘法，求出含有频率偏差量和频率变化率偏差量的未知参数矩阵；求出该数据窗内每个采样点的动态频率。本发明以多元泰勒级数为基础，建立电网信号测量矩阵方程，并且采用最小二乘法解决方程线性拟合问题，在跟踪频率的同时，还能对频率变化率进行监测，具有精确度高，实时性好和抗干扰性强等特点。
Description
技术领域
本发明属于电力系统动态频率测量技术，具体涉及一种最小二乘拟合动态频率测量方法。
背景技术
频率是电能质量的重要参数之一，能够反映电力系统瞬时运行状态的变化。从电力供应网络的全局来评价，频率是电网有功功率平衡状态的指示量。测量频率有助于监测电力系统电能发、输、配、用的全过程。从电能质量的角度来分析，频率是衡量电能质量的重要标准之一。频率偏高或者偏低都会影响设备的运行效率，甚至会对设备造成不可逆转的损害。以电力系统安全稳定运行为出发点和归宿点，频率是电力系统继电保护设备的重要控制参数之一。当电力系统处于非正常运行状态时，频率监测有助于辨识电力系统运行状态，并为继电保护装置、低频减载装置的动作提供判断依据。总的来说，实时、精确地频率测量有助于电网安全、稳定、经济运行。
目前，国内外大量学者已经在数字化测量技术的基础上提出了很多软件测频理论和算法。与硬件测量法相比，软件测频法不依赖于硬件电路，研究的重点主要在数学方法上。因而软件测频法灵活多变，适用场景广泛，能够匹配多样化的电网运行状态。一般而言，频率测量包括三个步骤：信号预处理、频率测量、结果再处理。其中，信号预处理主要是滤波环节和筛选数据，可以利用硬件电路或者滤波辅助算法达到预处理的目标。结果再处理同样为辅助算法，能够检验异常值或者通过辅助算法优化测量结果，使测量结果达到工程应用的要求。频率测量算法是频率测量步骤的核心部分。经典的测频算法主要有过零点检测法、函数解析法、函数逼近法、离散傅里叶变换算法及其改进算法、小波分析法等等。
上述算法都是在基本观测模型的基础上，不断改进电网信号的数学模型，目的在于逼近真实的物理信号。同时，上述算法在适用范围、测量精度、响应时间、算法复杂度、谐波抑制和噪声抗干扰等方面各有特点。但是，上述算法只对稳定频率测量有较好的效果，而对动态频率的测量缺乏支撑。
发明内容
发明目的：本发明的目的在于解决现有技术中存在的不足，提供一种最小二乘拟合动态频率测量方法，改进传统的最小二乘测频法，不仅可用于稳定频率的测量，而且在频率动态变化特别是频率线性变化的状态下，能够实现频率的实时测量；还能直接测量出频率变化率，增强对频率变化速率的监测。
技术方案：本发明所述的一种最小二乘拟合动态频率测量方法，依次包括以下步骤：
一种最小二乘拟合动态频率测量方法，其特征在于：依次包括以下步骤：
(1)设置移动数据窗的时间长度T0和采样时间间隔△t，将单个移动数据窗的电网连续信号通过A/D转换采样，并将第一个采样点的时间参数置零，从而形成该数据窗内的等时间间隔离散信号；
(2)利用FIR数字滤波器对待测离散信号进行工频基波提取；
(3)从工频基波提取后的离散信号中提取频率偏差和频率变化量偏差，对该数据窗内的采样点进行二元函数泰勒展开，构建量测状态模型；该模型的矩阵方程为：
[A][Y]＝[X]
[A]为由采样时间间隔△t、采样点序数n和预估频率f0组成的常系数矩阵，[Y]为含有待求频率参数组成的10×1的矩阵，[X]为N个采样点组成的N×1的矩阵，
(4)利用预估的频率值，离线算出常系数矩阵[A]；
(5)应用最小二乘法，按最小平方误差原理进行电网信号曲线拟合，求解量测状态矩阵方程，得到[Y]；
(6)应用公式求取频率偏差量△f的绝对值，频率偏差量的正负可以由公式进行判定；
应用公式求取频率变化率偏差量△k；
其中，y1，y2，y3，y4，y5，y6依次分别为[Y]中的第1至6列；
(7)该数据窗的最后一个采样点的动态频率可以根据该点的采样时间、频率偏差量和频率变化率偏差量求得，即公式f＝f0+△f+N△k△t；f0为预估频率，f为实际频率；
(8)当n<M，时间向前推进一个采样时间间隔，数据窗也随之剔除第1个采样点和加入新的采样点，返回步骤(3)，形成新的量测状态矩阵，进行最新采样点动态频率测量流程；n是指采样点序数，M是指采样点的数量；
(9)当n≥M，结束测量流程。
进一步的，实现实时性与精确度的统一，所述步骤(1)中的移动数据窗的时间长度T0取值范围是0.01-0.04s，T0不宜过长，会加大算法的复杂度。
进一步的，所述步骤(2)中的FIR数字滤波器采用通带频率范围在40～60Hz的8阶FIR数字滤波器，该滤波器的差分方程表达式为
z(n)＝0.02712x(n)+0.09165x(n-1)+0.17275x(n-2)+
0.23402x(n-3)+0.23402x(n-4)+0.17275x(n-5)
+0.09165x(n-6)+0.02712x(n-7)
式中，x(n)为第n个采样点的数值。
进一步的，所述步骤(3)中，基波状态下的量测状态模型为：
(3.1)在电力系统出现大量有功缺额的情况下，电网的频率变化是一个非常复杂的动态过程，若将电力系统视为等值的单机系统，频率特性可以用数学公式表达为：
Tf为系统频率变化过程中的时间常数，一般在4s～6s间变化；△f为频率偏差，满足公式△f＝f-f0，f为实际频率，f0为基波频率；△P为该单机系统的有功缺额；KL为负荷的频率调节系数；又因为对该式两边进行积分可得：
在非常短的时间内，将上式进行线性化处理，令频率变化率为k，则有f＝f'+kt，从而短时间内的频率动态变化问题可以转化为频率线性变化问题；电网信号函数可以表示为
x(t)为电网波形的单相电压或者电流，A1为基波峰值，f‘为待测量的基波频率，k为频率变化率，为信号基波的初相角；
(3.2)电网信号函数是关于频率变化率和频率的二元函数，由于频率关于时间的函数为正弦函数的隐函数，自变量f‘和k很难直接通过正弦函数解出，然后通过该二元函数的泰勒展开公式提取自变量△f和△k，从而便于建立自变量f‘和k的线性方程，单个采样点的线性方程可以表示为：
n为采样点的序号，△t为采样时间间隔，将数据窗内所有采样点的方程联立起来，分离未知量和常系数量，就构成了量测状态矩阵方程[A][Y]＝[X]，[X]为单个数据窗所有采样点数值组成的矩阵；
常系数矩阵[A]中的第n行可以表示为
未知矩阵[Y]中的参数为
(3.3)常系数矩阵[A]是由采样时间间隔△t、采样点序数n和预估频率f0构成，第一个数据窗的初始预估频率采用三点测频法进行估计，频率公式为：
其中，x(k-1)、x(k)和x(k+1)为相邻三个采样点的数值，第二个数据窗的初始预估频率采用第一个数据窗所测量的频率值，以此类推，每个数据窗的常系数矩阵可以离线求出，从而减少测量方法的复杂度。
进一步的，所述步骤(5)中，利用最小二乘法进行曲线拟合，通过最小化误差的平方和寻找矩阵方程的最优解，即使用公式[Y]＝{[A]T[A]}-1[A]T[X]确定含有频率和频率变化率的未知参数矩阵[Y]。
有益效果：发明不仅可以用于稳定频率的测量，而且在频率动态变化，特别是频率线性变化的状态下，能够实现频率的实时测量。同时，该方法还能直接测量出频率变化率，增强对频率变化速率的监测。
附图说明
图1为本发明的流程示意图。
具体实施方式
下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。
本发明利用二元函数泰勒展开提取电网信号模型的频率和频率变化率参数，在获取等时间间隔的离散信号后建立量测矩阵方程。然后，通过最小二乘法解决线性方程组的拟合问题，求取矩阵方程的参数。为了避免采样数据的递增对计算量带来的影响，设定固定时间的数据窗进行参数计算，在保证精度的同时，能够提高算法的实时性和对突变信号的敏感性。
实施例1：
如图1所示，本实施例步骤具体如下：
1、设置数据窗时间T＝0.02s和采样频率fs＝1600Hz，对电网信号进行采样。每个数据窗的采样点为24个，起始采样时间置零。
2、利用FIR数字滤波器对数据窗的信号进行处理，得到较为纯净的基波信号并形成信号矩阵[X]。
3、采用二元函数泰勒级数从待测离散信号中提取频率偏差和频率变化量偏差，对该数据窗内的采样点进行二元函数泰勒展开。单个采样点展开公式为：
式中，n为采样点的序列号，f0为预估频率。联立该数据窗所有采样点的展开公式，构建量测状态矩阵方程[A][Y]＝[X]。[A]为由时间和采样间隔函数组成的常系数矩阵，[Y]为含有待求频率参数组成的10×1的矩阵，[X]为24个采样点组成的24×1的矩阵。第一个数据窗的初始频率利用三点测频法进行估计，频率公式为式中，x(k-1)、x(k)和x(k+1)为相邻三个采样点的数值。计算结果f0＝50.0001Hz。第二个数据窗的初始预估频率采用第一个数据窗所测量的频率值，以此类推。每个数据窗的常系数矩阵可以离线求出，从而减少测量方法的复杂度。
4、应用最小二乘法拟合量测状态矩阵方程，求解未知参数矩阵[Y]。求解公式为:[Y]＝{[A]T[A]}-1[A]T[X]。
5、由于未知参数矩阵[Y]包含了频率偏差量和频率变化率偏差量，利用公式分别求取频率偏差量和频率变化率偏差量：
式中，y1，y2，y3，y4，y5，y6分别为[Y]中的第1、2、3、4、5、6列。再根据式f＝f0+△f+0.02*△k就可以得到该数据窗内第24个采样点的频率。
6、当n<N，时间向前推进一个采样时间间隔，数据窗也随之剔除第n-23个采样点和加入第n+1个采样点，返回步骤(2)，进行该数据窗的频率测量。N为所需测量频率的最后一个采样点的序号。
7、当n≥N，结束流程。
至此，完成了电网信号动态频率的实时测量和电网信号频率变化率的跟踪。
从上述实施例可以看出，本发明不仅可以适用于稳定频率的测量，而且在频率动态变化，特别是频率线性变化的状态下，能够实现频率的实时测量；本发明以多元泰勒级数为基础，建立电网信号测量矩阵方程，并且采用最小二乘法解决方程线性拟合问题，在跟踪频率的同时，还能对频率变化率进行监测，具有精确度高，实时性好和抗干扰性强等特点。经验证，本发明的误差在10-3Hz以内，测量速度略大于0.02s。此外，该算法还具有一定的抗干扰性和鲁棒性。
Claims (5)
1.一种最小二乘拟合动态频率测量方法，其特征在于：依次包括以下步骤：
(1)设置移动数据窗的时间长度T0和采样时间间隔Δt，将单个移动数据窗的电网连续信号通过A/D转换采样，并将第一个采样点的时间参数置零，从而形成该数据窗内的等时间间隔离散信号；
(2)利用FIR数字滤波器对待测离散信号进行工频基波提取；
(3)从工频基波提取后的离散信号中提取频率偏差和频率变化量偏差，对该数据窗内的采样点进行二元函数泰勒展开，构建量测状态模型；该模型的矩阵方程为：
[A][Y]＝[X]
[A]为由采样时间间隔Δt、采样点序数n和预估频率f0组成的常系数矩阵，[Y]为含有待求频率参数组成的10×1的矩阵，[X]为该数据窗口内N个采样点组成的N×1的矩阵，
(4)利用预估的频率值，离线算出常系数矩阵[A]；
(5)应用最小二乘法，按最小平方误差原理进行电网信号曲线拟合，求解量测状态矩阵方程，得到[Y]；
(6)应用公式求取频率偏差量Δf的绝对值，频率偏差量的正负可以由公式进行判定；
应用公式求取频率变化率偏差量Δk；
其中，y1，y2，y3，y4，y5，y6依次分别为[Y]中的第1至6列；
(7)该数据窗的最后一个采样点的动态频率可以根据该点的采样时间、频率偏差量和频率变化率偏差量求得，即公式f＝f0+Δf+NΔkΔt；f0为预估频率，f为实际频率；
(8)当n<M，时间向前推进一个采样时间间隔，数据窗也随之剔除第1个采样点和加入新的采样点，返回步骤(3)，形成新的量测状态矩阵，进行最新采样点动态频率测量流程；n是指采样点序数，M是指采样点的数量；
(9)当n≥M，结束测量流程。
2.根据权利要求1所述的最小二乘法拟合动态频率测量方法，其特征在于：所述步骤(1)中的移动数据窗的时间长度T0取值范围是0.01-0.04s。
3.根据权利要求1所述的最小二乘法拟合动态频率测量方法，其特征在于：所述步骤(2)中的FIR数字滤波器采用通带频率范围在40～60Hz的8阶FIR数字滤波器，该滤波器的差分方程表达式为
z(n)＝0.02712x(n)+0.09165x(n-1)+0.17275x(n-2)+
0.23402x(n-3)+0.23402x(n-4)+0.17275x(n-5)
+0.09165x(n-6)+0.02712x(n-7)
式中，x(n)为第n个采样点的数值。
4.根据权利要求1所述的最小二乘法拟合动态频率测量方法，其特征在于：所述步骤(3)中，基波状态下的量测状态模型为：
(3.1)在电力系统出现大量有功缺额的情况下，电网的频率变化是一个非常复杂的动态过程，若将电力系统视为等值的单机系统，频率特性可以用数学公式表达为：
<mrow>
<msub>
<mi>T</mi>
<mi>f</mi>
</msub>
<mfrac>
<mrow>
<mi>d</mi>
<mi>&Delta;</mi>
<mi>f</mi>
</mrow>
<mrow>
<mi>d</mi>
<mi>t</mi>
</mrow>
</mfrac>
<mo>+</mo>
<mi>&Delta;</mi>
<mi>f</mi>
<mo>=</mo>
<mfrac>
<mrow>
<mi>&Delta;</mi>
<mi>P</mi>
</mrow>
<msub>
<mi>K</mi>
<mi>L</mi>
</msub>
</mfrac>
</mrow>
Tf为系统频率变化过程中的时间常数，一般在4s～6s间变化；Δf为频率偏差，满足公式Δf＝f-f0，f为实际频率，f0为基波频率；ΔP为该单机系统的有功缺额；KL为负荷的频率调节系数；又因为对该式两边进行积分可得：
<mrow>
<mi>f</mi>
<mo>=</mo>
<msub>
<mi>f</mi>
<mn>0</mn>
</msub>
<mo>+</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>f</mi>
<mo>&Proportional;</mo>
</msub>
<mo>-</mo>
<msub>
<mi>f</mi>
<mn>0</mn>
</msub>
<mo>)</mo>
</mrow>
<mo>&lsqb;</mo>
<mn>1</mn>
<mo>-</mo>
<mi>exp</mi>
<mrow>
<mo>(</mo>
<mo>-</mo>
<mfrac>
<mi>t</mi>
<msub>
<mi>T</mi>
<mi>f</mi>
</msub>
</mfrac>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
在非常短的时间内，将上式进行线性化处理，令频率变化率为k，则有f＝f'+kt，从而短时间内的频率动态变化问题可以转化为频率线性变化问题；电网信号函数可以表示为
x(t)为电网波形的单相电压或者电流，A1为基波峰值，f‘为待测量的基波频率，k为频率变化率，为信号基波的初相角；
(3.2)电网信号函数是关于频率变化率和频率的二元函数，由于频率关于时间的函数为正弦函数的隐函数，自变量f‘和k很难直接通过正弦函数解出，然后通过该二元函数的泰勒展开公式提取自变量Δf和Δk，从而便于建立自变量f‘和k的线性方程，单个采样点的线性方程可以表示为：
n为采样点的序号，Δt为采样时间间隔，将数据窗内所有采样点的方程联立起来，分离未知量和常系数量，就构成了量测状态矩阵方程[A][Y]＝[X]，[X]为单个数据窗所有采样点数值组成的矩阵；
常系数矩阵[A]中的第n行可以表示为
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>n</mi>
<mn>1</mn>
</mrow>
</msub>
<mo>=</mo>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msub>
<mi>&pi;f</mi>
<mn>0</mn>
</msub>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>n</mi>
<mn>2</mn>
</mrow>
</msub>
<mo>=</mo>
<mi>cos</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msub>
<mi>&pi;f</mi>
<mn>0</mn>
</msub>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>n</mi>
<mn>3</mn>
</mrow>
</msub>
<mo>=</mo>
<mn>2</mn>
<mi>&pi;</mi>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
<mi> </mi>
<mi>cos</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msub>
<mi>&pi;f</mi>
<mn>0</mn>
</msub>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>n</mi>
<mn>4</mn>
</mrow>
</msub>
<mo>=</mo>
<mo>-</mo>
<mn>2</mn>
<mi>&pi;</mi>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
<mi> </mi>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msub>
<mi>&pi;f</mi>
<mn>0</mn>
</msub>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>n</mi>
<mn>5</mn>
</mrow>
</msub>
<mo>=</mo>
<mn>2</mn>
<msup>
<mi>&pi;n</mi>
<mn>2</mn>
</msup>
<msup>
<mi>&Delta;t</mi>
<mn>2</mn>
</msup>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msub>
<mi>&pi;f</mi>
<mn>0</mn>
</msub>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>n</mi>
<mn>6</mn>
</mrow>
</msub>
<mo>=</mo>
<mn>2</mn>
<msup>
<mi>&pi;n</mi>
<mn>2</mn>
</msup>
<msup>
<mi>&Delta;t</mi>
<mn>2</mn>
</msup>
<mi>cos</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msub>
<mi>&pi;f</mi>
<mn>0</mn>
</msub>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>n</mi>
<mn>7</mn>
</mrow>
</msub>
<mo>=</mo>
<mo>-</mo>
<mn>4</mn>
<msup>
<mi>&pi;</mi>
<mn>2</mn>
</msup>
<msup>
<mi>n</mi>
<mn>3</mn>
</msup>
<msup>
<mi>&Delta;t</mi>
<mn>3</mn>
</msup>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msub>
<mi>&pi;f</mi>
<mn>0</mn>
</msub>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>n</mi>
<mn>8</mn>
</mrow>
</msub>
<mo>=</mo>
<mo>-</mo>
<mn>4</mn>
<msup>
<mi>&pi;</mi>
<mn>2</mn>
</msup>
<msup>
<mi>n</mi>
<mn>3</mn>
</msup>
<msup>
<mi>&Delta;t</mi>
<mn>3</mn>
</msup>
<mi>cos</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msub>
<mi>&pi;f</mi>
<mn>0</mn>
</msub>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>n</mi>
<mn>9</mn>
</mrow>
</msub>
<mo>=</mo>
<mo>-</mo>
<mn>2</mn>
<msup>
<mi>&pi;</mi>
<mn>2</mn>
</msup>
<msup>
<mi>n</mi>
<mn>4</mn>
</msup>
<msup>
<mi>&Delta;t</mi>
<mn>4</mn>
</msup>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msub>
<mi>&pi;f</mi>
<mn>0</mn>
</msub>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>a</mi>
<mrow>
<mi>n</mi>
<mn>10</mn>
</mrow>
</msub>
<mo>=</mo>
<mo>-</mo>
<mn>2</mn>
<msup>
<mi>&pi;</mi>
<mn>2</mn>
</msup>
<msup>
<mi>n</mi>
<mn>4</mn>
</msup>
<msup>
<mi>t</mi>
<mn>4</mn>
</msup>
<mi>cos</mi>
<mrow>
<mo>(</mo>
<mrow>
<mn>2</mn>
<msub>
<mi>&pi;f</mi>
<mn>0</mn>
</msub>
<mi>n</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
未知矩阵[Y]中的参数为
(3.3)常系数矩阵[A]是由采样时间间隔Δt、采样点序数n和预估频率f0构成，第一个数据窗的初始预估频率采用三点测频法进行估计，频率公式为：
<mrow>
<msub>
<mi>f</mi>
<mn>0</mn>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
</mfrac>
<mi>arccos</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mn>2</mn>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
其中，x(k-1)、x(k)和x(k+1)为相邻三个采样点的数值，第二个数据窗的初始预估频率采用第一个数据窗所测量的频率值，以此类推，每个数据窗的常系数矩阵可以离线求出，从而减少测量方法的复杂度。
5.根据权利要求1所述的最小二乘法拟合动态频率测量方法，其特征在于：所述步骤(5)中，利用最小二乘法进行曲线拟合，通过最小化误差的平方和寻找矩阵方程的最优解，即使用公式[Y]＝{[A]T[A]}-1[A]T[X]确定含有频率和频率变化率的未知参数矩阵[Y]。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710383645.1A CN107271768B (zh) | 2017-05-26 | 2017-05-26 | 一种最小二乘拟合动态频率测量方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710383645.1A CN107271768B (zh) | 2017-05-26 | 2017-05-26 | 一种最小二乘拟合动态频率测量方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107271768A true CN107271768A (zh) | 2017-10-20 |
| CN107271768B CN107271768B (zh) | 2019-06-21 |
Family
ID=60064296
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710383645.1A Active CN107271768B (zh) | 2017-05-26 | 2017-05-26 | 一种最小二乘拟合动态频率测量方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107271768B (zh) |
Cited By (16)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108627694A (zh) * | 2018-04-17 | 2018-10-09 | 中国农业大学 | 一种电网电压频率的检测方法、装置和设备 |
| CN108776262A (zh) * | 2018-06-04 | 2018-11-09 | 西南交通大学 | 一种考虑带外干扰的电力系统频率测量方法 |
| CN109374966A (zh) * | 2018-10-23 | 2019-02-22 | 国网重庆市电力公司电力科学研究院 | 一种电网频率估计方法 |
| CN109490630A (zh) * | 2018-11-22 | 2019-03-19 | 华北电力大学 | 一种基于矩阵束的动态相量测量方法 |
| CN109541303A (zh) * | 2018-12-10 | 2019-03-29 | 华北电力大学 | 一种相角正弦调制信号频率和频率变化率的补偿方法 |
| CN110470903A (zh) * | 2019-07-31 | 2019-11-19 | 山东建筑大学 | 一种电压频率软测量装置和方法 |
| CN112109583A (zh) * | 2020-09-11 | 2020-12-22 | 国网山东省电力公司临沂供电公司 | 一种电网谐波畸变处理方法 |
| CN112269054A (zh) * | 2020-09-16 | 2021-01-26 | 国网安徽省电力有限公司六安供电公司 | 基于改进Prony的功率自适应算法 |
| CN112414526A (zh) * | 2020-11-13 | 2021-02-26 | 迪比（重庆）智能科技研究院有限公司 | 散袋小包装中药的快速称量方法 |
| CN112557751A (zh) * | 2020-12-03 | 2021-03-26 | 东南大学 | 一种基于dft迭代法的谐波参数估计方法 |
| CN112948755A (zh) * | 2021-04-01 | 2021-06-11 | 中国空空导弹研究院 | 一种遥测正弦参数判读方法 |
| CN112965365A (zh) * | 2021-02-23 | 2021-06-15 | 浙江中智达科技有限公司 | Pid控制回路的模型辨识方法、装置、系统及存储介质 |
| CN114034927A (zh) * | 2021-11-01 | 2022-02-11 | 南京国电南自电网自动化有限公司 | 一种跟频插值采样的信号测量方法及系统 |
| CN114217126A (zh) * | 2021-12-09 | 2022-03-22 | 广西电网有限责任公司电力科学研究院 | 一种电力系统暂态信号瞬时频率识别算法 |
| CN114236232A (zh) * | 2021-12-16 | 2022-03-25 | 广州城市理工学院 | 一种考虑频率变化趋势的小水电频率预测方法 |
| CN114355360A (zh) * | 2021-12-14 | 2022-04-15 | 中船航海科技有限责任公司 | 一种多普勒计程仪回波频率估计方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101661069A (zh) * | 2009-09-25 | 2010-03-03 | 北京四方继保自动化股份有限公司 | 不依赖状态矩阵的弱可观非pmu测点动态过程实时估计方法 |
| CN102305891A (zh) * | 2011-07-04 | 2012-01-04 | 武汉大学 | 一种电力系统低频振荡在线监测方法 |
| CN102495281A (zh) * | 2011-12-14 | 2012-06-13 | 北京易事特电源有限公司 | 一种电力系统相量频率测量方法 |
| CN102565531A (zh) * | 2012-02-07 | 2012-07-11 | 天津大学 | 压电换能器动态参数测量仪及测量方法 |
| CN103323855A (zh) * | 2012-03-22 | 2013-09-25 | 中国科学院电子学研究所 | 一种基线动态测量系统的精度获取方法 |
-
2017
- 2017-05-26 CN CN201710383645.1A patent/CN107271768B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101661069A (zh) * | 2009-09-25 | 2010-03-03 | 北京四方继保自动化股份有限公司 | 不依赖状态矩阵的弱可观非pmu测点动态过程实时估计方法 |
| CN102305891A (zh) * | 2011-07-04 | 2012-01-04 | 武汉大学 | 一种电力系统低频振荡在线监测方法 |
| CN102495281A (zh) * | 2011-12-14 | 2012-06-13 | 北京易事特电源有限公司 | 一种电力系统相量频率测量方法 |
| CN102565531A (zh) * | 2012-02-07 | 2012-07-11 | 天津大学 | 压电换能器动态参数测量仪及测量方法 |
| CN103323855A (zh) * | 2012-03-22 | 2013-09-25 | 中国科学院电子学研究所 | 一种基线动态测量系统的精度获取方法 |
Non-Patent Citations (2)
| Title |
|---|
| 孔庆鹏等: ""基于分段最小二乘拟合的瞬时频率估计方法"", 《农业机械学报》 * |
| 王仲根等: ""基于最小二乘拟合和特征基函数法的目标宽带RCS快速计算"", 《安徽师范大学学报》 * |
Cited By (24)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108627694B (zh) * | 2018-04-17 | 2020-02-11 | 中国农业大学 | 一种电网电压频率的检测方法、装置和设备 |
| CN108627694A (zh) * | 2018-04-17 | 2018-10-09 | 中国农业大学 | 一种电网电压频率的检测方法、装置和设备 |
| CN108776262A (zh) * | 2018-06-04 | 2018-11-09 | 西南交通大学 | 一种考虑带外干扰的电力系统频率测量方法 |
| CN109374966A (zh) * | 2018-10-23 | 2019-02-22 | 国网重庆市电力公司电力科学研究院 | 一种电网频率估计方法 |
| CN109490630A (zh) * | 2018-11-22 | 2019-03-19 | 华北电力大学 | 一种基于矩阵束的动态相量测量方法 |
| CN109490630B (zh) * | 2018-11-22 | 2020-11-10 | 华北电力大学 | 一种基于矩阵束的动态相量测量方法 |
| CN109541303A (zh) * | 2018-12-10 | 2019-03-29 | 华北电力大学 | 一种相角正弦调制信号频率和频率变化率的补偿方法 |
| CN110470903A (zh) * | 2019-07-31 | 2019-11-19 | 山东建筑大学 | 一种电压频率软测量装置和方法 |
| CN110470903B (zh) * | 2019-07-31 | 2021-08-17 | 山东建筑大学 | 一种电压频率软测量装置和方法 |
| CN112109583A (zh) * | 2020-09-11 | 2020-12-22 | 国网山东省电力公司临沂供电公司 | 一种电网谐波畸变处理方法 |
| CN112109583B (zh) * | 2020-09-11 | 2021-09-10 | 国网山东省电力公司临沂供电公司 | 一种电网谐波畸变处理方法 |
| CN112269054A (zh) * | 2020-09-16 | 2021-01-26 | 国网安徽省电力有限公司六安供电公司 | 基于改进Prony的功率自适应算法 |
| CN112414526A (zh) * | 2020-11-13 | 2021-02-26 | 迪比（重庆）智能科技研究院有限公司 | 散袋小包装中药的快速称量方法 |
| CN112414526B (zh) * | 2020-11-13 | 2022-07-01 | 迪比（重庆）智能科技研究院有限公司 | 散袋小包装中药的快速称量方法 |
| CN112557751A (zh) * | 2020-12-03 | 2021-03-26 | 东南大学 | 一种基于dft迭代法的谐波参数估计方法 |
| CN112965365A (zh) * | 2021-02-23 | 2021-06-15 | 浙江中智达科技有限公司 | Pid控制回路的模型辨识方法、装置、系统及存储介质 |
| CN112948755A (zh) * | 2021-04-01 | 2021-06-11 | 中国空空导弹研究院 | 一种遥测正弦参数判读方法 |
| CN112948755B (zh) * | 2021-04-01 | 2023-06-23 | 中国空空导弹研究院 | 一种遥测正弦参数判读方法 |
| CN114034927A (zh) * | 2021-11-01 | 2022-02-11 | 南京国电南自电网自动化有限公司 | 一种跟频插值采样的信号测量方法及系统 |
| CN114217126A (zh) * | 2021-12-09 | 2022-03-22 | 广西电网有限责任公司电力科学研究院 | 一种电力系统暂态信号瞬时频率识别算法 |
| CN114217126B (zh) * | 2021-12-09 | 2023-10-24 | 广西电网有限责任公司电力科学研究院 | 一种电力系统暂态信号瞬时频率识别算法 |
| CN114355360A (zh) * | 2021-12-14 | 2022-04-15 | 中船航海科技有限责任公司 | 一种多普勒计程仪回波频率估计方法 |
| CN114355360B (zh) * | 2021-12-14 | 2024-05-03 | 中船航海科技有限责任公司 | 一种多普勒计程仪回波频率估计方法 |
| CN114236232A (zh) * | 2021-12-16 | 2022-03-25 | 广州城市理工学院 | 一种考虑频率变化趋势的小水电频率预测方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107271768B (zh) | 2019-06-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107271768B (zh) | 一种最小二乘拟合动态频率测量方法 | |
| CN107478990B (zh) | 一种发电机机电暂态过程动态估计方法 | |
| CN106199183B (zh) | 一种实现次同步振荡在线辨识告警的pmu和方法 | |
| Khoshkhoo et al. | On-line dynamic voltage instability prediction based on decision tree supported by a wide-area measurement system | |
| CN102163844B (zh) | 基于相量测量装置的电力系统状态检测方法 | |
| CN103401238B (zh) | 一种基于总体测辨法的电力负荷建模方法 | |
| CN103941088A (zh) | 一种基于三相信号的电力系统频率快速测量方法 | |
| CN106405230B (zh) | 频率测量方法和装置 | |
| CN108574291A (zh) | 一种基于集合卡尔曼滤波发电机动态状态估计方法 | |
| CN104865474A (zh) | 一种基于pmu数据实时监测低频振荡源的方法 | |
| Khodaparast | A review of dynamic phasor estimation by non-linear Kalman filters | |
| CN107565553A (zh) | 一种基于ukf的配电网抗差动态状态估计方法 | |
| CN103983847B (zh) | 一种同步相量测量中基于rls的自适应频率跟踪测量方法 | |
| CN112564157B (zh) | 一种抑制连续换相失败的定熄弧角控制改进方法 | |
| CN108333426A (zh) | 基于傅氏算法的电力系统频率测量方法 | |
| CN102967779A (zh) | 一种输电线路分布参数的辨识方法 | |
| CN104749488A (zh) | 一种基于连续数据窗的直流线路时域故障测距方法 | |
| CN106443253A (zh) | 一种基于pmu数据的输电线路参数辨识方法 | |
| Ren et al. | Use of recursive wavelet transform for estimating power system frequency and phasors | |
| CN102522742B (zh) | 基于单点量测信息的外网戴维南等值参数的估计方法 | |
| CN103972889B (zh) | 一种配电线路阻抗在线辨识方法 | |
| CN109218073A (zh) | 一种计及网络攻击和参数不确定性的动态状态估计方法 | |
| CN109039208A (zh) | 一种开关磁阻电机增量电感特性在线检测方法 | |
| CN111965409B (zh) | 基于分段差分波形有效值的电压暂态扰动检测方法 | |
| CN103308759A (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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |