CN109754013B - 一种基于无迹卡尔曼滤波的电力系统混合量测融合方法 - Google Patents

一种基于无迹卡尔曼滤波的电力系统混合量测融合方法 Download PDF

Info

Publication number
CN109754013B
CN109754013B CN201811653402.6A CN201811653402A CN109754013B CN 109754013 B CN109754013 B CN 109754013B CN 201811653402 A CN201811653402 A CN 201811653402A CN 109754013 B CN109754013 B CN 109754013B
Authority
CN
China
Prior art keywords
measurement
value
scada
pmu
formula
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.)
Expired - Fee Related
Application number
CN201811653402.6A
Other languages
English (en)
Other versions
CN109754013A (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.)
Zhejiang University ZJU
State Grid Shanghai Electric Power Co Ltd
Original Assignee
Zhejiang University ZJU
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 Zhejiang University ZJU, State Grid Shanghai Electric Power Co Ltd filed Critical Zhejiang University ZJU
Priority to CN201811653402.6A priority Critical patent/CN109754013B/zh
Publication of CN109754013A publication Critical patent/CN109754013A/zh
Application granted granted Critical
Publication of CN109754013B publication Critical patent/CN109754013B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

本发明公开了一种基于无迹卡尔曼滤波的电力系统混合量测融合方法，针对目前μPMU量测与SCADA量测长期共存但难以同时使用的现状，引入无迹卡尔曼滤波算法实现混合量测的融合。分析两种量测采样频率及量测类型存在的差异，实现混合量测的同步化。提出了一种自适应无迹卡尔曼滤波器AUKF，该滤波器通过采用比例修正最小偏度单形采样策略解决传统UKF计算量大且易产生采样的非局部效应等问题，并通过实现自适应选取比例因子来提高滤波精度。本发明方法能够实现混合量测融合，且融合后数据有较高精度。

Description

一种基于无迹卡尔曼滤波的电力系统混合量测融合方法
技术领域
本发明属于电力系统领域，尤其涉及一种基于无迹卡尔曼滤波的电力系统混合量测融合方法。
背景技术
目前，基于GPS的微型同步相量测量装置(micro Phasor Measurement Unit，μPMU)逐步应用于配电网中，但短期内μPMU无法取代监控及数据采集系统(SupervisoryControl And Data Acquisition，SCADA)实现全网可观测要求，因此由μPMU和SCADA组成的混合量测是改善以单一SCADA量测进行状态估计及其他高级应用的一种有效手段。由于μPMU和SCADA系统是采用不同的方法为其技术支撑，两种量测数据存在较多差异，若不经任何处理就混合，必将出现一系列兼容性问题，不但不能利用μPMU量测改善数据质量，反而会降低应用的性能，因此，需对混合量测进行融合处理后将其提供给后续应用。
在数据融合领域中，广泛应用的融合技术有经典推理法、卡尔曼滤波法、贝叶斯估计法、物理模型法、支持向量机等，而应用最广泛的是卡尔曼滤波法。对非线性系统，多为扩展卡尔曼滤波(extend Kalman filter，EKF)，虽然EKF有许多改进方法，但EKF有其先天局限性，即在计算Jacobian矩阵时不可避免地引入线性化误差。
由于EKF存在的线性化误差问题，Julier等提出了无迹卡尔曼滤波(unscentedKalman filter，UKF)方法。UKF摒弃了对非线性函数线性化的做法，对非线性函数的概率密度分布进行近似，用一系列确定样本来逼近状态的后验概率密度，不依靠Jacobian矩阵进行线性化。UKF没有忽略高阶项，因此对于非线性分布的统计量有较高的计算精度，有效克服了EKF的滤波精度低、稳定性差的缺陷。
发明内容
本发明的目的在于针对现有技术的不足，提供一种基于无迹卡尔曼滤波的电力系统混合量测融合方法。
本发明的目的是通过以下技术方案来实现的：一种基于无迹卡尔曼滤波的电力系统混合量测融合方法，该方法包括以下步骤：
(1)通过μPMU及SCADA采集电力系统节点的实时量测，μPMU量测中包含相量量测，SCADA量测中只包含幅值量测。
(2)通过计算两类量测的相关度及量测时标对齐实现异步量测同步化，具体包括以下子步骤：
(2.1)计算μPMU量测及SCADA量测之间相同部分的等级差数d，由量测变换方法计算混合量测异同部分的等效等级差数D，通过两等级差数计算混合量测的相关度ρ。
(2.2)将相关度最大的μPMU量测时标赋予对应的SCADA量测，实现量测时标对齐。
(3)将同步化后的μPMU量测与SCADA量测采用自适应无迹卡尔曼滤波算法进行融合，具体包括以下子步骤；
(3.1)以同步化后的μPMU量测与SCADA量测共有部分为状态量，采用自适应最小偏度单形采样策略对混合量测进行UT变换，变换公式如下：
选取权值初值
Figure GDA0002729408270000021
设置比例修正因子初值α、先验系数β，则第i个状态量对应的采样点权值为
Figure GDA0002729408270000022
Figure GDA0002729408270000023
计算初始向量
Figure GDA0002729408270000024
则状态量为n时的向量由j＝2,3,…,n迭代计算得到，向量递推公式如下
Figure GDA0002729408270000025
Figure GDA0002729408270000026
表示经j轮递推后第i个状态值的采样点偏差范围。
由式(4)生成的Sigma点集为
Figure GDA0002729408270000027
式中
Figure GDA0002729408270000028
为协方差矩阵P的Cholesky因子，x为UT变换前的值。
(3.2)通过Holt指数平滑法计得到系统状态转移函数f(·)，计算k时刻第i个状态量的sigma采样点集χi,k|k的一步预测值χi,k+1|k及系统一步预测协方差Pk+1|k，公式如下：
χi,k+1|k＝f(χi,k|k)+wk (6)
Figure GDA0002729408270000029
Figure GDA0002729408270000031
式中，wk为系统过程噪声，xk+1|k为系统状态的一步预测值，Qk为过程噪声协方差矩阵。
(3.3)对一步预测值xk+1|k再次进行UT变换，产生新的Sigma点集χ'i,k+1,k，将χ'i,k+1,k点集带入观测方程h(·)得到Sigma点集的观测预测值zi,k+1,k，并加权求得观测预测均值
Figure GDA0002729408270000032
自协方差
Figure GDA0002729408270000033
及互协方差
Figure GDA0002729408270000034
通过修正因子α对自协方差进行修正
zi,k+1,k＝h(χ'i,k+1,k)+vk (9)
Figure GDA0002729408270000035
Figure GDA0002729408270000036
Figure GDA0002729408270000037
(3.4)计算Kalman增益矩阵Kk+1
Figure GDA0002729408270000038
(3.5)更新状态及协方差
Figure GDA0002729408270000039
Figure GDA00027294082700000310
(3.6)自适应比例修正因子选取方法
由式(16)估算出xk|k与x之间的距离dk
Figure GDA00027294082700000311
设置变化因子μ，则αk+1的值式(17)、(18)确定
Figure GDA00027294082700000312
αk+1＝αk(1-μn)1/n (18)
由式(16)、(17)、(18)计算出αk+1，将αk+1作为下一时刻UT过程的比例因子，这样可把Sigma采样点限制在合理范围，从而避免产生采样的非局部效应；
本发明以SCADA量测为滤波状态向量结构，利用一系列确定样本点来逼近状态量的后验概率分布，因此滤波结果xk+1|k+1即为混合量测融合后的输出结果。
进一步地，所述步骤(2.1)中等级差数d及等效等级差数D公式如下：
Figure GDA0002729408270000041
Figure GDA0002729408270000042
式中ZS(i)为节点i的SCADA量测，
Figure GDA0002729408270000043
为节点i的PMU量测中与SCADA量测结构相同的部分，
Figure GDA0002729408270000044
为节点i的PMU量测中与SCADA量测异同部分经量测变换后的等效量测，Z'S(i)为SCADA量测与变换后PMU量测结构相同的部分，变换公式如下：
Figure GDA0002729408270000045
Figure GDA0002729408270000046
式中
Figure GDA0002729408270000047
为节点iPMU量测中的电压相量量测，
Figure GDA0002729408270000048
为节点iPMU量测中的电流相量量测的共轭值，Yij为节点i、j之间的支路导纳，Yi0为节点i对地导纳，
Figure GDA0002729408270000049
为变换后节点i到j的等效有功和无功，Uj为变换后节点j的等效电压。
进一步地，所述步骤(2.1)中μPMU及SCADA采集的量测的相关度ρ公式如下：
Figure GDA00027294082700000410
式中n为两列量测中相同元素的等级个数，m为两列量测中异同元素的等级个数，ω1、ω2分别为相同部分及异同部分的权值，di为量测相同部分对应元素的等级差数，Dj为量测异同部分对应元素的等效等级差数。
进一步地，所述步骤(3.2)中的Holt指数平滑法具体如下：
假设k时刻系统状态的一步预测值为xk|k-1，系统滤波值为xk|k，则k+1时刻的预测值为：
Figure GDA00027294082700000411
由Kalman滤波可知一步预测表达式为
xk+1|k＝Fkxk|k+Gk (25)
由式(24)、(25)可得
Figure GDA00027294082700000412
式中ak为水平分量，bk为倾斜分量，αH、βH为平滑参数。
由Sigma点集χi，k|k、式(6)～(8)和式(24)～(26)计算一步预测值xk+1|k及一步预测协方差Pk+1|k。
进一步地，所述步骤(3.3)中的自协方差修正方法如下：
Figure GDA0002729408270000051
λ＝1+β-α2 (28)
式中λ为协方差修正系数，z0,k+1,k为sigma点集观测预测值zi,k+1,k中第一个点的值。
本发明的有益效果是：为解决μPMU量测无法直接同SCADA量测同时用于电力系统，本发明分析两类量测的非同步性误差，计算μPMU量测与SCADA量测的相关度，将相关度最大μPMU量测时标赋予SCADA量测，实现混合量测的同步化，引入UKF对两种量测进行融合处理，采用最小偏度单行采样策略以及改进比例参数选取，提出AUKF算法，通过自适应选取比例因子，有效提高了滤波精度，同时也降低了UT变换的非局部效应。
附图说明
图1为同步化与非同步化量测误差比较图；
图2为AUKF滤波误差及量测误差对比图；
图3为IEEE 118节点系统中三种算法的RMSE对比图；
图4为融合值用于IEEE 118节点系统状态估计后的估计值RMSE比较，(a)为状态估计电压幅值RMSE，(b)为状态估计电压相角RMSE。
具体实施方式
下面结合附图和具体实施例对本发明作进一步详细说明。
本发明提供的一种基于无迹卡尔曼滤波的电力系统混合量测融合方法，引入UKF滤波方法进行电力系统混合量测融合，由于SCADA与μPMU量测的更新周期及结构均有所不同，因此在量测融合之前，需对混合量测进行同步化；同时改进UT变换过程中的比例因子选取，提出一种自适应无迹卡尔曼滤波器(adaptive unscented Kalman filter，AUKF)；最后在IEEE-14节点系统与IEEE-118节点系统中进行仿真试验，验证了AUKF算法的有效性。
1.无迹卡尔曼滤波原理具体如下：
考虑如下离散的非线性系统：
xk+1＝f(xk)+wk (1)
zk+1＝h(xk+1)+vk+1 (2)
式中k为离散时间，xk为系统在k时刻的状态，zk为对应状态的观测信号，f(xk)为描述系统状态变化的状态转移方程，h(xk+1)为描述系统状态与量测量关系的量测方程，wk为k时刻系统过程噪声，vk为k时刻系统观测噪声。
设wk和vk是互不相关的非均值高斯白噪声，且具有时变统计特性，即
Figure GDA0002729408270000061
式中δkj为Kronecker函数，qk、rk分别为噪声wk、vk的期望值，Qk为过程噪声协方差矩阵，Rk为观测噪声协方差矩阵。
1.1UT变换过程
UT变换是在原状态分布中选择适合的采样策略来选取一定采样点，使选取的采样点均值和协方差等于原状态分布的均值和协方差，将采样点进行非线性变换，使得变换后的均值和协方差具有三阶甚至更高阶精度。采样点的选择基于先验均值和先验协方差矩阵的平方根相关列实现。
对于一个n维分布的状态空间，对称采样需要2n+1个采样点，而最小偏度单形采样只需要n+1个采样点(若考虑中心点为n+2个采样点)即可，虽然单形采样精度略低于对称采样，但二者皆可达到三阶精度，且单形采样实时性高于对称采样。对于非线性变换y＝f(x)，状态向量x为n维随机变量，且已知均值
Figure GDA00027294082700000610
和方差P，具体Sigma点采样策略如下：
1)计算采样点相应的权值
选取权值
Figure GDA0002729408270000062
剩余权值为
Figure GDA0002729408270000063
式中
Figure GDA0002729408270000064
为计算系统状态预测值及观测预测值的采样点权值，
Figure GDA0002729408270000065
为计算误差协方差的采样点权值。
2)计算n+2个Sigma点
初始向量(状态维数j＝1时情况)为
Figure GDA0002729408270000066
输入维数j＝2,3,…,n时，向量递推公式如下
Figure GDA0002729408270000067
Figure GDA0002729408270000068
表示经j轮递推后第i个状态值的采样点偏差范围。
由式(6)生成的Sigma点集为
Figure GDA0002729408270000069
式中
Figure GDA0002729408270000071
为协方差矩阵P的Cholesky因子。
1.2UKF滤波过程
1)预测
xk为k时刻状态量，Pk为k时刻状态量的协方差，由式(4)-(7)得Sigma点集{χi,k}及其相应权值
Figure GDA0002729408270000072
由非线性变换y＝f(x)及当前系统状态sigma点χi,k|k计算k时刻第i个状态量的sigma采样点集χi,k|k的一步预测值χi,k+1|k及系统一步预测协方差Pk+1|k
χi,k+1|k＝f(χi,k|k)+wk (8)
Figure GDA0002729408270000073
Figure GDA0002729408270000074
2)更新
依据一步预测值，再次进行UT变换，产生新的Sigma点集χ'i,k+1,k，将χ'i,k+1,k点集带入观测方程得到Sigma点集的观测预测值zi,k+1,k，并加权求得观测预测均值
Figure GDA0002729408270000075
自协方差
Figure GDA0002729408270000076
及互协方差
Figure GDA0002729408270000077
zi,k+1,k＝h(χ'i,k+1,k)+vk (11)
Figure GDA0002729408270000078
Figure GDA0002729408270000079
Figure GDA00027294082700000710
计算Kalman增益矩阵Kk+1
Figure GDA00027294082700000711
状态及协方差更新
Figure GDA00027294082700000712
Figure GDA00027294082700000713
1.3混合量测时延误差分析
目前μPMU量测的采样频率一般为100测量点/s，而SCADA量测采样频率达到0.1-5测量点/s，假设μPMU及SCADA系统采集量测的时刻均为t，μPMU量测时延为ΔtP，SCADA量测时延为ΔtS，且有ΔtP＜ΔtS。由于SCADA量测时延远大于μPMU量测，那么当t时刻的SCADA量测ZS(t)到达时，μPMU量测的到达值为Zp(t')，又因为目前系统中只有少部分节点配置μPMU装置，而SCADA量测不具备时标，无法将两类量测在时间上直接同步。因此，若直接将两类量测同时使用即同时取用最新到达的两类量测值，μPMU量测反而会增大融合值的误差，具体分析如下：
考虑混合量测的非同步性误差，那么式(3)中的观测噪声矩阵可表示为
Figure GDA0002729408270000081
式中
Figure GDA0002729408270000082
为测量噪声矩阵，
Figure GDA0002729408270000083
为量测不同步形成的噪声矩阵，其中
Figure GDA0002729408270000084
Figure GDA0002729408270000085
式中Z'P为与SCADA量测同时到达的PMU量测，Zp为与SCADA量测同时采集的PMU量测，
Figure GDA0002729408270000086
为采集时刻的量测真实值，由(18)～(20)计算并整理得
Figure GDA0002729408270000087
若单独考虑非同步性误差
Figure GDA0002729408270000088
那么有
Figure GDA0002729408270000089
对比式(21)、(22)不难发现
Figure GDA00027294082700000810
因此，量测不同步不仅会引入非同步性误差，同时还会放大系统原有误差，从而降低滤波精度。
2.混合量测融合算法
2.1异步量测同步化
由1.3分析可知，要实现混合量测融合，须对两种量测进行同步化。由于SCADA量测不具备时标，且μPMU的采样频率远高于SCADA，那么可以确定在SCADA采集时刻的一个较小时间差内一定存在一个μPMU量测。因此运用Spearman相关系数法分析混合量测中的同类型量测，将相关系数最大的μPMU量测时标赋予SCADA量测，实现异步量测的同步化。具体函数如下：
Figure GDA00027294082700000811
式中n为两列量测中相同元素的等级个数，m为两列量测中异同元素的等级个数，di为量测相同部分对应元素的等级差数，Dj为量测异同部分对应元素的等效等级差数，ω1、ω2分别为相同部分及异同部分的权值，其中等级差数d及等效等级差数D公式如下：
Figure GDA00027294082700000812
Figure GDA0002729408270000091
Figure GDA0002729408270000092
式中ZS(i)为节点i的SCADA量测，
Figure GDA0002729408270000093
为节点i的PMU量测中与SCADA量测结构相同的部分，
Figure GDA0002729408270000094
为节点i的PMU量测中与SCADA量测异同部分经量测变换后的等效量测，Z'S(i)为SCADA量测与变换后PMU量测结构相同的部分，变换公式如下：
Figure GDA0002729408270000095
Figure GDA0002729408270000096
式中
Figure GDA0002729408270000097
为节点i PMU量测中的电压相量量测，
Figure GDA0002729408270000098
为节点i PMU量测中的电流相量量测的共轭值，Yij为节点i、j之间的支路导纳，Yi0为节点i对地导纳，
Figure GDA0002729408270000099
为变换后节点i到j的等效有功和无功，Uj为变换后节点j的等效电压。
为降低同步化过程的计算量，提高算法的实时性，取初始时刻至第一个SCADA量测到达时刻之间所有μPMU量测进行同步化。对后续同步化过程，取前一次同步后SCADA量测时标加上SCADA系统的采集周期为下一次同步时刻预测值，计算预测时刻前后各十个μPMU量测与SCADA量测的相关系数，实现同步化。
对配置μPMU的节点，将与SCADA量测相关度最大的μPMU量测的时标赋予SCADA量测，实现两类量测的同步化；对未配置μPMU的节点，可求取某时刻所有配置μPMU的节点同步化后SCADA量测时标的期望，将期望时标赋予其他SCADA量测。
2.2自适应比例修正的最小偏度单形采样
经同步化过程后，采用UKF算法实现混合量测融合，而UKF算法的核心是采样策略的选取，因此选取单行采样策略并进行适当改进。在单形采样策略中，随着状态维数增加，Sigma点到中心点
Figure GDA00027294082700000914
的距离也相应增加，因此随着迭代次数的增加，便会产生采样的非局部效应，因此采用比例修正采样策略。单形采样修正部分为
Figure GDA00027294082700000910
Figure GDA00027294082700000911
Figure GDA00027294082700000912
Figure GDA00027294082700000913
其中α∈(0,1]为比例因子，通过选取合适的α值可降低UT过程的非局部效应；β用于描述状态的先验分布信息，在高斯分布时，通常取β＝2；
Figure GDA0002729408270000101
为修正前采样点权值，
Figure GDA0002729408270000102
和
Figure GDA0002729408270000103
为修正后采样点均值和协方差权值。
传统比例修正为选取固定α值，较未修正算法性能有一定的提高，但随着迭代次数增加，固定α值对系统性能的提升作用也随之减小，因此对修正策略进行改进，提出一种自适应选取α值的修正方法。对不同类型的数据，在取得最优滤波效果时α的取值也有较大差异，因此选取合适的α初值后，依据运行状态修正的α值，实现滤波效果的提升。
由UKF滤波过程可知，k时刻的滤波值xk|k为k+1时刻一步预测采样中心点，而k时刻滤波更新后的误差协方差矩阵Pk|k描述了系统滤波值xk|k与真实值x之间的关系，即
Pk|k＝E[(xk|k-x)(xk|k-x)T] (33)
由式(33)可估算出xk|k与x之间的距离
Figure GDA0002729408270000104
UKF算法过程是将xk|k取代x进行k+1时刻一步预测时Sigma采样。加入变化因子μ，则αk+1的值由下式确定
αk+1＝αk(1-μn)1/n (34)
Figure GDA0002729408270000105
由式(33)～(35)可计算出αk+1，将αk+1作为该时刻UT过程的比例因子，这样可把Sigma采样点限制在合理范围，从而避免产生采样的非局部效应。
2.3AUKF算法过程
将改进的单形采样策略运用于UKF算法，取为AUKF算法，AUKF算法过程与传统UKF算法相似，主要由预测和更新两部分组成，在每次Sigma点采集之前进行修正因子αk的计算。
1)预测
取前两次误差估计协方差Pk|k及Pk-1|k-1，由式(33)～(35)计算变化因子μ和比例修正因子αk。
以前一次滤波值xk|k为中心点，由(29)～(31)及比例修正因子αk构造Sigma点集χi，k|k及相应权值
Figure GDA0002729408270000106
系统动态模型采用Holt指数平滑法。假设k时刻系统状态的一步预测值为xk|k-1，系统滤波值为xk|k，则k+1时刻的预测值为：
Figure GDA0002729408270000107
由Kalman滤波可知一步预测表达式为
xk+1|k＝Fkxk|k+Gk (37)
由式(36)、(37)可得
Figure GDA0002729408270000111
式中ak为水平分量，bk为倾斜分量，αH、βH为平滑参数。
由Sigma点集χi，k|k、式(8)～(10)和式(36)～(38)计算一步预测值xk+1|k及一步预测协方差Pk+1|k。
2)更新
由一步预测值依据式(29)～(32)及修正因子αk构造新的Sigma点集χ'i,k+1,k，依据式(11)将预测状态Sigma值χ'i,k+1,k转换为量测估计Sigma值zi,k+1,k，依据式(12)～(14)计算预测量测均值及协方差，由式(15)～(17)计算Kalman增益并对状态量和协方差估计值进行更新。
本发明在IEEE-14节点系统及IEEE-118节点系统对同步化过程及AUKF算法进行仿真验证，通过电力系统分析软件获得两测试系统的潮流计算数据(标幺值)并作为真值，在真值的基础上，添加服从高斯分布的随机扰动形成量测值。
为了验证本发明方法的有效性，选取最小偏度单行采样UKF、对称采样UKF与AUKF算法进行比较，并定义如下几种性能指标。
滤波误差绝对值、量测误差绝对值：
Figure GDA0002729408270000112
Figure GDA0002729408270000113
最大绝对滤波误差：
Figure GDA0002729408270000114
及均方根误差：
Figure GDA0002729408270000115
式中：
Figure GDA0002729408270000116
xk,i、
Figure GDA0002729408270000117
zk,i分别为k时刻第i个滤波值、真实值、量测值、量测真实值，n为滤波值维数。其中滤波误差绝对值ε1及量测误差绝对值ε2用于对同步化过程及AUKF的滤波效果进行评价；最大绝对误差
Figure GDA0002729408270000118
和均方根误差RMSE用于同步化后三种算法的性能对比。
3.2IEEE 14节点系统测试结果
使用IEEE 14节点系统分析同步化过程效果，并验证AUKF算法的滤波效果，测试过程参数选取如下：SCADA量测中电压量测的误差标准差取0.01，均值为0；功率量测的误差标准差取0.02，均值为0。μPMU电压量测及电流幅值量测的误差标准差取0.002，均值为0，相角量测的误差标准差取0.001，均值为0；μPMU功率量测的误差标准差取0.004，均值为0。假设在节点4配置μPMU。Holt两参数法的取值为αH＝0.775，βH＝0.7。采样修正系数初值取值为α＝0.9、β＝2，权值为w0＝0.5。
将未同步化的混合量测与同步化后混合量测分别用传统UKF算法进行量测融合，融合结果如图1所示，可以看出，同步化过程能有效降低量测的非同步性误差。
图2给出了AUKF滤波后的各类滤波值误差及量测误差对比，电压量测误差标准差及功率量测标准差分别为给定值0.002和0.004，经计算滤波值均方根误差分别为0.000868和0.001418。
表1均方根误差比较
算法 最大电压误差 电压误差均方根 最大功率误差 功率误差均方根
对称采样UKF 0.003062 0.001239 0.014191 0.002805
单形采样UKF 0.004216 0.001541 0.012384 0.002936
AUKF 0.002376 0.000868 0.004114 0.001418
表1给出了对称采样UKF、单行采样UKF和AUKF的均方根误差及最大误差，对比三种滤波算法结果，其误差均值都在一定范围之内，但AUKF算法对误差的滤波程度相较于其他两种算法都有显著提升，说明AUKF可以随着时间变化选取最优比例因子从而使滤波值更接近真值。
3.3 IEEE-118节点系统测试结果
为了验证AUKF在大规模系统的适用性，在IEEE 118节点系统分别用对称采样UKF、单行采样UKF和AUKF对同步化后量测进行同条件测试，测试参数同IEEE 14节点系统一致，同时将测试结果分别用于状态估计，通过状态估计结果，对算法的性能进行进一步验证。配置μPMU的节点为：4，8，17，19，25，28，32，38，40，44，48，55，60，65，72，75，79，80，85，90，101，106，111。
图3为三种算法在IEEE 118节点系统中的性能比较，图中性能指标RMSE也进一步验证了理论分析，即AUKF滤波效果优于其他两种算法，说明AUKF能有效减小量测噪声。
量测融合后，将三种算法的融合结果同未配置PMU节点的量测值组成IEEE 118节点系统量测，用状态估计进行测试，以性能指标RMSE对运算结果进行评价，评价结果示于图4。
从图4所示结果可以看出，经AUKF算法滤波后的融合值用于状态估计均能取得较优效果，但从电压幅值及相角的RMSE结果比较可知，经AUKF算法滤波后的融合值用于状态估计取得的效果明显优于其他两种算法，与理论分析一致。
上述实施例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。

Claims (5)

1.一种基于无迹卡尔曼滤波的电力系统混合量测融合方法，其特征在于，该方法包括以下步骤：
(1)通过μPMU及SCADA采集电力系统节点的实时量测，μPMU量测中包含相量量测，SCADA量测中只包含幅值量测；
(2)通过计算两类量测的相关度及量测时标对齐实现异步量测同步化，具体包括以下子步骤：
(2.1)计算μPMU量测及SCADA量测之间相同部分的等级差数dg，由量测变换方法计算混合量测异同部分的等效等级差数Dh，通过两等级差数计算混合量测的相关度ρ；
(2.2)将相关度最大的μPMU量测时标赋予对应的SCADA量测，实现量测时标对齐；
(3)将同步化后的μPMU量测与SCADA量测采用自适应无迹卡尔曼滤波算法进行融合，具体包括以下子步骤；
(3.1)以同步化后的μPMU量测与SCADA量测共有部分为状态量，采用自适应最小偏度单形采样策略对混合量测进行UT变换，变换公式如下：
选取权值初值
Figure FDA0003071698140000011
设置比例修正因子初值α、先验系数β，则第i个采样点对应的修正后的状态权值
Figure FDA0003071698140000012
与协方差权值
Figure FDA0003071698140000013
分别为
Figure FDA0003071698140000014
Figure FDA0003071698140000015
计算初始向量
Figure FDA0003071698140000016
则状态量为n时的向量由j＝2,3,…,n迭代计算得到，向量递推公式如下
Figure FDA0003071698140000017
Figure FDA0003071698140000018
表示经j轮递推后第i个状态值的采样点偏差范围；
由式(4)生成的Sigma点集为
Figure FDA0003071698140000021
式中
Figure FDA0003071698140000022
为协方差矩阵P的Cholesky因子，x为UT变换前的值；
(3.2)通过Holt指数平滑法计算得到系统状态转移函数f(·)，计算k时刻第i个状态量的sigma采样点集χi,k|k的k+1时刻一步预测值χi,k+1|k及系统一步预测协方差Pk+1|k，公式如下：
χi,k+1|k＝f(χi,k|k)+wk (6)
Figure FDA0003071698140000023
Figure FDA0003071698140000024
式中，wk为系统过程噪声，xk+1|k为k+1时刻系统状态的一步预测值，Qk为过程噪声协方差矩阵；
(3.3)对一步预测值xk+1|k再次进行UT变换，产生新的Sigma点集χ′i,k+1,k，将χ′i,k+1,k点集带入观测方程h(·)得到Sigma点集的观测预测值zi,k+1,k，并加权求得观测预测均值
Figure FDA0003071698140000025
自协方差
Figure FDA0003071698140000026
及互协方差
Figure FDA0003071698140000027
通过修正因子α对自协方差进行修正
zi,k+1,k＝h(χ′i,k+1,k)+vk (9)
式中，vk为k时刻系统观测噪声；
Figure FDA0003071698140000028
Figure FDA0003071698140000029
Figure FDA00030716981400000210
(3.4)计算Kalman增益矩阵Kk+1
Figure FDA00030716981400000211
(3.5)更新状态及协方差
Figure FDA00030716981400000212
Figure FDA00030716981400000213
(3.6)自适应比例修正因子选取方法
由式(16)估算出xk|k与x之间的距离dk
Figure FDA00030716981400000214
设置变化因子μ，则αk+1的值式(17)、(18)确定
Figure FDA0003071698140000031
αk+1＝αk(1-μn)1/n (18)
由式(16)、(17)、(18)计算出αk+1，将αk+1作为下一时刻UT过程的比例因子；
以SCADA量测为滤波状态向量结构，利用一系列确定样本点来逼近状态量的后验概率分布，因此滤波结果xk+1|k+1即为混合量测融合后的输出结果。
2.根据权利要求1所述的一种基于无迹卡尔曼滤波的电力系统混合量测融合方法，其特征在于，所述步骤(2.1)中等级差数dg及等效等级差数Dh公式如下：
Figure FDA0003071698140000032
Figure FDA0003071698140000033
式中ZS(g)为节点g的SCADA量测，
Figure FDA0003071698140000034
为节点g的PMU量测中与SCADA量测结构相同的部分，
Figure FDA0003071698140000035
为节点g的PMU量测中与SCADA量测异同部分经量测变换后的等效量测，Z′S(g)为SCADA量测与变换后PMU量测结构相同的部分，变换公式如下：
Figure FDA0003071698140000036
Figure FDA0003071698140000037
式中
Figure FDA0003071698140000038
为节点gPMU量测中的电压相量量测，
Figure FDA0003071698140000039
为节点gPMU量测中的电流相量量测的共轭值，Ygh为节点g、h之间的支路导纳，Yh0为节点g对地导纳，
Figure FDA00030716981400000310
为变换后节点g到h的等效有功和无功，Uh为变换后节点h的等效电压。
3.根据权利要求1所述的一种基于无迹卡尔曼滤波的电力系统混合量测融合方法，其特征在于，所述步骤(2.1)中μPMU及SCADA采集的量测的相关度ρ公式如下：
Figure FDA00030716981400000311
式中p为两列量测中相同元素的等级个数，q为两列量测中异同元素的等级个数，ω1、ω2分别为相同部分及异同部分的权值，dg为量测相同部分对应元素的等级差数，Dh为量测异同部分对应元素的等效等级差数。
4.根据权利要求1所述的一种基于无迹卡尔曼滤波的电力系统混合量测融合方法，其特征在于，所述步骤(3.2)中的Holt指数平滑法具体如下：
k时刻系统状态的一步预测值为xk|k-1，系统滤波值为xk|k，则k+1时刻的预测值为：
Figure FDA0003071698140000041
由Kalman滤波可知一步预测表达式为
xk+1|k＝Fkxk|k+Gk (25)
由式(24)、(25)可得
Figure FDA0003071698140000042
式中ak为水平分量，bk为倾斜分量，αH、βH为平滑参数；
由Sigma点集χi，k|k、式(6)～(8)和式(24)～(26)计算一步预测值xk+1|k及一步预测协方差Pk+1|k。
5.根据权利要求1所述的一种基于无迹卡尔曼滤波的电力系统混合量测融合方法，其特征在于，所述步骤(3.3)中的自协方差修正方法如下：
Figure FDA0003071698140000043
λ＝1+β-α2 (28)
式中λ为协方差修正系数，z0,k+1,k为sigma点集观测预测值zi,k+1,k中第一个点的值。
CN201811653402.6A 2018-12-31 2018-12-31 一种基于无迹卡尔曼滤波的电力系统混合量测融合方法 Expired - Fee Related CN109754013B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811653402.6A CN109754013B (zh) 2018-12-31 2018-12-31 一种基于无迹卡尔曼滤波的电力系统混合量测融合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811653402.6A CN109754013B (zh) 2018-12-31 2018-12-31 一种基于无迹卡尔曼滤波的电力系统混合量测融合方法

Publications (2)

Publication Number Publication Date
CN109754013A CN109754013A (zh) 2019-05-14
CN109754013B true CN109754013B (zh) 2021-09-10

Family

ID=66404995

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811653402.6A Expired - Fee Related CN109754013B (zh) 2018-12-31 2018-12-31 一种基于无迹卡尔曼滤波的电力系统混合量测融合方法

Country Status (1)

Country Link
CN (1) CN109754013B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110676940B (zh) * 2019-10-21 2021-02-02 国网上海市电力公司 一种提高参数辨识精度的配电网pmu配置方法及系统
CN110752622B (zh) * 2019-12-12 2023-12-05 燕山大学 一种配电网仿射状态估计方法
CN111181529B (zh) * 2020-01-17 2022-02-08 中山大学 一种应用于非线性高斯模型的平滑约束扩展卡尔曼滤波方法
CN116303786B (zh) * 2023-03-18 2023-10-27 上海圈讯科技股份有限公司 一种基于多维数据融合算法的区块链金融大数据管理系统
CN117713750A (zh) * 2023-12-14 2024-03-15 河海大学 一种基于分数次幂的一致性卡尔曼滤波状态估计方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107679133A (zh) * 2017-09-22 2018-02-09 电子科技大学 一种实用于海量实时pmu数据的挖掘方法
CN108574291A (zh) * 2018-04-23 2018-09-25 河海大学 一种基于集合卡尔曼滤波发电机动态状态估计方法
CN108759838A (zh) * 2018-05-23 2018-11-06 安徽科技学院 基于秩卡尔曼滤波器的移动机器人多传感器信息融合方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103633739B (zh) * 2013-11-28 2015-05-20 中国科学院广州能源研究所 一种微电网能量管理系统和方法
US9746511B2 (en) * 2015-11-25 2017-08-29 Hitachi, Ltd. Estimating the locations of power system events using PMU measurements
CN106707061A (zh) * 2016-12-16 2017-05-24 湖南大学 基于混合量测的配电网动态状态估计方法
CN106844952A (zh) * 2017-01-20 2017-06-13 河海大学 基于无迹粒子滤波理论的发电机动态状态估计方法
CN107607977B (zh) * 2017-08-22 2020-12-08 哈尔滨工程大学 一种基于最小偏度单形采样的自适应ukf组合导航方法
CN107565553A (zh) * 2017-09-19 2018-01-09 贵州大学 一种基于ukf的配电网抗差动态状态估计方法
CN108448568B (zh) * 2018-03-08 2020-03-27 国网山东省电力公司潍坊供电公司 基于多种时间周期测量数据的配电网混合状态估计方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107679133A (zh) * 2017-09-22 2018-02-09 电子科技大学 一种实用于海量实时pmu数据的挖掘方法
CN108574291A (zh) * 2018-04-23 2018-09-25 河海大学 一种基于集合卡尔曼滤波发电机动态状态估计方法
CN108759838A (zh) * 2018-05-23 2018-11-06 安徽科技学院 基于秩卡尔曼滤波器的移动机器人多传感器信息融合方法

Also Published As

Publication number Publication date
CN109754013A (zh) 2019-05-14

Similar Documents

Publication Publication Date Title
CN109754013B (zh) 一种基于无迹卡尔曼滤波的电力系统混合量测融合方法
CN110659722B (zh) 基于AdaBoost-CBP神经网络的电动汽车锂离子电池健康状态估算方法
CN107563550B (zh) 一种基于pmu的配电网实时分布式状态估计及pmu的优化配置方法
CN108155648A (zh) 基于自适应h无穷扩展卡尔曼滤波的状态估计方法
CN107658881A (zh) 基于戴维南等值方法的电压稳定临界点判断方法
CN111426957B (zh) 一种模拟车辆工况动力电池荷电状态soc估算优化方法
CN108985373B (zh) 一种多传感器数据加权融合方法
CN110907702A (zh) 一种改进动态谐波估计方法和系统
CN107025609A (zh) 基于奇异值分解cdkf的电力系统动态状态估计方法
CN114567288B (zh) 基于变分贝叶斯的分布协同非线性系统状态估计方法
CN113536623B (zh) 一种材料不确定性结构稳健性拓扑优化设计方法
CN107204616B (zh) 基于自适应稀疏伪谱法的电力系统随机状态估计方法
CN110490378A (zh) 基于云scada大数据的电网状态估计精度的计算方法
CN114239796A (zh) 一种基于扩展卡尔曼滤波的电力系统状态估计方法
CN110632521B (zh) 一种锂离子电池容量的融合估计方法
CN109782188B (zh) 航天器电池组模拟器的soc特性校准方法
CN109447512B (zh) 基于均匀设计的大电网可靠性评估方法
CN108228959A (zh) 利用删失数据估计系统实际状态的方法及应用其的滤波器
CN109638811B (zh) 基于模型等值的配电网电压功率灵敏度鲁棒估计方法
CN113553538B (zh) 一种递推修正混合线性状态估计方法
CN106100609B (zh) 单状态变量和两级Kalman滤波器时间尺度算法
CN107967395A (zh) 一种基于beta小波基函数展开的时变非线性系统快速辨识方法
Gokpinar et al. Generalization of inclusion probabilities in ranked set sampling
Mercieca et al. Unscented Rauch-Tung-Striebel smoothing for nonlinear descriptor systems
CN111597766A (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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210910

CF01 Termination of patent right due to non-payment of annual fee