作者简介:江 辉(1968—),女,深圳大学教授.研究方向:智能电网、电能质量分析与控制.E-mail:huijiang@szu.edu.cn
中文责编:英 子; 英文责编:雨 辰
1)深圳大学光电工程学院,广东深圳 518060; 2)深圳大学机电与控制工程学院,广东深圳 518060
1)College of Optoelectronic Engineering, Shenzhen University, Shenzhen 518060, Guangdong Province, P.R.China2)College of Mechatronics and Control Engineering, Shenzhen University, Shenzhen 518060, Guangdong Province, P.R.China
power system; power quality; dynamic harmonic estimation; particle filter; unscented Kalman filter; unscented particle filter
DOI: 10.3724/SP.J.1249.2016.01080
提出一种基于无迹粒子滤波(unscented particle filter, UPF)算法的电网动态谐波估计方法.通过无迹卡尔曼滤波算法得到电网动态谐波状态量的估计值和协方差,运用这些结果改进传统粒子滤波算法的重要密度函数,采用粒子滤波算法得到电网动态谐波的最优估计值.该方法克服了无迹卡尔曼滤波算法(unscented Kalman filter, UKF)对噪声要求为高斯分布的限制和传统粒子滤波(particle filter,PF)算法易退化的缺点,保留了UKF对非线性问题的较好处理和PF强抗干扰性能力.仿真结果表明,在高斯噪声和非高斯噪声情况下,UPF算法得到的电网动态谐波幅值、相位的估计值都更接近真实值.
This paper proposes a new method for estimating dynamic harmonics in power systems based on the unscented particle filter(UPF)algorithm. The UPF algorithm is used to estimate the values and covariance of the state variables of the dynamic harmonics in power systems. These estimated values are used to generate the importance density function of particle filter(PF)algorithm. The optimal estimation of dynamic harmonics in the power system is achieved by using the derived PF algorithm. The proposed method overcomes not only the restriction of Gaussian distribution noise required in unscented Kalman filters(UKF)but also the drawback of easy degeneration for conventional PF. In addition, it retains the high performance of UKF in processing nonlinearity and the strong anti-interference capability of PF.Experiments show that the estimates of dynamic harmonic amplitude and phase by the proposed method are closer to the true values under both Gaussian noise and non-Gaussian noise situations.
电网动态谐波是衡量电网电能质量的一个重要指标[1],对动态谐波的有效治理有助于电网电压的实时监控,准确掌握电压的动态变化,从而采取有效的技术手段对电网电压进行调节.为此,肖湘宁等[2]应用峰值电压法和基波分量法等进行电压凹陷的检测.江辉等[3]提出一种基于小波变换和改进S变换的电压扰动分类方法.丁宁等[4]采用有效值法(root mean square,RMS)分析不同情形下电压凹陷的情况[4].王效孟等[5]提出一种改进的有效值法——半周期有效值(half a cycle of RMS,HRMS)法对电压扰动进行检测,该方法相比传统的有效值法,检测误差更小,延时更短.秦英林等[6]给出了基于奇异值分解技术的电压扰动检测方法,对信号进行线性分解,得到扰动发生的起点和终点时刻.侯世英等[7]采用αβ求导检测法来检测单相电压骤降特征量.Singh等[8]利用双线性递归最小二乘法(bilinear recursive least square, BRLS)对静态谐波和动态谐波的幅值、相位进行估计.任文林等[9]提出一种基于卡尔曼滤波(Kalman filter,KF)的电压扰动检测方法.闵伟等[10]给出一种改进型Sage-Husa卡尔曼滤波电压暂降检测方法,实时更新系统噪声和量测噪声方差,提高了滤波的稳定性和计算速度.卡尔曼滤波是基于线性系统提出的,但实际中广泛存在的是非线性系统,从而卡尔曼滤波在电能质量分析中存在困难.扩展卡尔曼(extended Kalman filter,EKF)[11]将非线性系统线性化,文献[12-13]运用EKF算法对电网频率进行估计.但EKF算法是对非线性函数的泰勒展开式截取实现线性化,且需要计算Jacobian 矩阵的导数,对于高维的数据处理非常复杂.因而,Julier等[14]提出无迹卡尔曼滤波(unscented Kalman filter,UKF),通过一组精确选择的Sigma点匹配随机量的统计特性,UKF没有Jacobian矩阵计算问题,使算法的实现比EKF 更为容易,估计精度更高.
UKF模型要求系统的状态噪声和量测噪声为高斯分布,才能得到最优的状态估计,当条件不符合时其滤波和预测精度就很难得到保证.基于粒子群优化的无迹卡尔曼滤波算法(particle swarm optimized unscented Kalman filter, PSOUKF)[15]采用粒子群算法优化高斯噪声下的状态噪声和量测噪声协方差,提高了估计精度,但是仍对非高斯噪声情形下预测精度难以保证.粒子滤波(particle filter,PF)[16]突破了卡尔曼滤波的理论框架,系统的状态噪声和量测噪声不再受高斯分布限制,被广泛用于机动目标跟踪[17]、人导航与定位[18]、工业过程与监视[19]等非线性、非高斯的情形中.无迹粒子滤波(unscented particle filter,UPF),结合UKF和PF算法的优点,利用UKF对非线性的处理能力以及PF的强抗干扰能力,因而对强干扰性的非线性系统有良好的滤波效果.本研究运用UPF对电网动态谐波进行估计.
一个含噪声和衰减的N次谐波信号,离散化后可表示为[20]
y(k)=∑Nn=1Ansin(ωnkTs+φn)+
Adcexp(-αdckTs)+μ(k)(1)
其中, N为谐波的次数; An为幅值; ωn=n2πf0, f0为基波频率; φn为谐波的相位; Ts为采样周期; k=t/Ts; Adcexp(-αdckTs)为衰减部分, Adc和αdc是常量; μ(k)是加入的随机噪声.将Adcexp(-αdckTs)用二阶泰勒公式近似展开,可得
Adcexp(-αdckTs)=Adc-AdcαdckTs(2)
将式(2)代入式(1),则
y(k)=∑Nn=1Ansin(ωnkTs+φn)+
Adc-AdcαdckTs+μ(k)(3)
式(3)可表示为式(4)的量测方程形式
y(k)=H(k)X(k)+μ(k)(4)
其中,H(k)=[sin(ω1kTs)cos(ω1kTs)sin(ω2kTs) cos(ω2kTs)… sin(ωNkTs)cos(ωNkTs)1 -kTs]为测量转移矩阵; 状态量可表示为
X(k)=[X1(k)X2(k)… X2N-1(k)
X2N(k)Adc Adcαdc]T(5)
该处X2N-1(k)=ANcos φN, X2N(k)=ANsin φN. 根据式(5)可将谐波信号状态方程表示为
X(k+1)=F(k)X(k)(6)
其中, F(k)=[1 0 … 0
0 1 … 0
0 0 … 1]为状态转移矩阵.
由状态量X(k), 可得各次谐波的幅值及相位分别为
An=(X22n+X22n-1)1/2(7)
φn=arctan((X2n)/(X2n-1))(8)
UKF由EKF发展而来,是对卡尔曼滤波的一种改进算法,不需要像EKF那样要求把非线性系统线性化,也不需要计算Jacobian矩阵.其主要思想是,在确保采样均值和协方差前提下,选择一组样本点(sigma点集),每个样本点都有自己的状态值和加权值,将非线性变换应用于采样的每个sigma点,得到非线性转换后的sigma点集,从而得到新的均值和协方差[14].
PF[21]是一种基于蒙特卡罗的近似贝叶斯滤波算法,其主要思想是,根据系统状态量的经验条件分布,寻找一组在状态空间里传播的随机样本集合(这样的样本集合称为粒子),然后根据量测不断地改变粒子的权重和位置来修正最初的经验条件分布,实质是用粒子及其权重组成的离散随机测度来近似系统随机变量的概率密度函数,从而获得最小方差估计.
在常规算法中,常选择先验概率密度作为重要密度函数,这样虽然便于实现,但也易丢失k时刻的测量值,使状态量过分依赖模型,若模型不准确将不能有效估计真实分布.随着计算的推进,易引起粒子算法退化.为避免发生退化,在引入重采样的同时,还可选择一个适合的重要密度函数.
UPF利用UKF得到一个建议分布函数,然后利用该函数来代替PF算法中的重要密度函数, UPF算法主要包括选取sigma采样点、状态以及量测的更新、粒子集和权值的计算[22]4个方面.
非线性系统状态方程和量测方程分别为
xk=F(xk-1, νk)(9)
yk=H(xk-1, wk)(10)
其中, xk, νk, yk和wk分别是系统在k时刻的状态量、状态噪声、量测量和量测噪声; F和 H分别是系统状态转移矩阵和量测转移矩阵.
设 xik-1和 Pik-1是k-1时刻状态量 xk-1的均值和误差协方差, L是 xik-1的维数,粒子滤波中取M个粒子数, i=1,2,…,M, j=0,1,…,L. UPF算法计算步骤为
1)选取每个粒子的sigma采样点.
χik-1,j=[xik-1, xik-1+(((L+λ)Pik-1)1/2)j,
xik-1-(((L+λ)Pik-1)1/2)j](11)
其中,(((L+λ)Pik-1)1/2)j为矩阵(L+λ)Pik-1平方根的第j行, λ=α2(L+κ)-L是调节参数,控制sigma点与均值的距离.按照经验κ一般取0或3-L,10-4≤α≤1.
2)对每个粒子进行时间更新.利用式(9)传递sigma点χik|k-1,j,计算状态量预测值(-overx)ik|k-1和其协方差pik|k-1.
χik|k-1,j=F(χik|k-1,j)(12)
(-overx)ik|k-1=∑2Lj=0ωmjχik|k-1,j(13)
pik|k-1=∑2Lj=0[ωcj(χik|k-1-(-overx)ik|k-1)×
(χik|k-1-(-overx)ik|k-1)T]+Q(14)
通过式(10)传递量测量的采样点Y ik|k-1,j,计算出量测量预测值(-overy)ik|k-1及其协方差Pykyk
Y ik|k-1,j=H(χik|k-1,j)(15)
(-overy)ik|k-1=∑2Lj=0ωmjY ik|k-1,j(16)
Pykyk=∑2Lj=0[ωcj(Y ik|k-1,j-(-overy)ik|k-1)×
(Y ik|k-1,j-(-overy)ik|k-1)T]+R(17)
Q和R分别为两互不相干的系统状态噪声协方差和量测噪声协方差,其均值为0. ωmj和ωcj分别为分配到sigma点的均值和协方差的权重系数,按下式计算
{ωm0=λ/(L+λ),
ωc0=λ/(L+λ)+(1-α2+β)
ωmj=ωcj=λ/(2(L+λ))(18)
3)量测更新.计算状态量测协方差Pxkyk和UPF算法的增益KUPF,得到估计值 xik及误差协方差 Pik分别为
Pxkyk=∑2Lj=0[ωcj(χik|k-1,j-(-overx) ik|k-1)×
(Yik|k-1,j-(-overy)ik|k-1)T](19)
KUPF=PxkykP-1ykyk(20)
xik=(-overx)ik|k-1+KUPF(yk-(-overy)ik|k-1)(21)
Pik=Pik|k-1-KUPFPykykKUPFT(22)
至此,利用UKF算法得到每个粒子在k时刻的估计值和误差协方差 xik和 Pik, 从而得到重要密度函数为
q(xik|xi0:k-1,y1:k)=N(xik,Pik)(23)
其中, N(·)表示高斯函数,由此重要密度函数中采样粒子,计算权值并归一化.采样粒子为
x^ik~q(xik|xi0:k-1,y1:k)=N(xik,Pik)(24)
令x^i0:k~(xi0:k-1,x^ik), 计算权值
w*ik=w*ik-1(p(yk|x^ik)p(x^ik|x^ik-1))/(q(x^ik|x^ik-1,yk))(25)
归一化权值 wik=w*ik/∑Mi=1w*ik, p(·)为状态转移概率密度,取高斯分布[23].
4)重采样按照下式近似计算有效粒子数.
Neff=1/(∑Mi=1(wik)2)(26)
若Neff>Nth(Nth为设定的有效样本数,一般取2M/3), 则进行下一步; 否则,重新采样[ 21],并将原带权值的粒子映射为等权的M个粒子,即wik=1/M.
5)状态估计.
状态估计为 x^k=∑Mi=1wikxik(27)
方差估计为
Pk=∑Mi=1wik(xik-x^ik)(xik-x^ik)T(28)
返回步骤2)进行递推计算,直到结束.
1)根据前述谐波模型,对式(5)中的状态量设置初始值,从每个状态量中抽取M个粒子.设置状态噪声协方差Q、 量测噪声协方差R和状态误差协方差初始值P0.
2)对每个粒子,根据式(11)至式(22)计算k刻的状态量的估计值 xik和协方差 Pik.
3)由式(23)得出粒子滤波的重要性密度函数,根据该函数抽取采样粒子,并通过式(25)至式(28)得到最终的状态量最优估计值 x^k和协方差Pk.
4)最优估计值x^k中的分量按式(7)和式(8)计算出该时刻谐波的幅值和相位.令k值加1, 并转到本节步骤2)进行下一时刻的参数估计.
为比较分别采用KF、UKF和UPF的谐波估计方法的性能,本研究采用均方误差(mean square error,MSE)来衡量估计的效果.第k个时刻的均方误差定义为
MSE=1/N∑Mi=1(yi'-yi)2(29)
其中, N为状态量的维数; yi'和yi分别是该时刻预测值和理论值的第i个分量.显然,均方误差越小说明估计效果越好.
本研究是在Matlab R2007版本下进行仿真,采样频率fs=2.5 kHz.
含有噪声的谐波电压信号表达式描述为
y(t)=1.2sin(ωt+75°)+0.8sin(3ωt+55°)+
0.2sin(5ωt+45°)+μ(t)(30)
其中, μ(t)是随机高斯噪声,其均值为0,方差为1; t=0.2 s,采样间隔ts=1/2 500 s.根据式(4)和式(6), L=6, M取100,3种方法得到的5次谐波幅值及相位的结果如图1和图2.
在图1和图2中,可明显看出,基于UPF的谐波估计方法较其他两种方法能更迅速地估计出真实值,且更接近真实值.表1为在0.05 s时刻时,各方法静态谐波幅值与相位估计结果的对比.由表1可见,在0.05 s时,5次谐波的幅值的实际值为0.2 p.u..基于KF、 UKF 和UPF估计的幅值分别为0.191 3、0.201 5和0.201 0 p.u.,其相对误差分别为4.35%、0.75%和0.51%,可见基于UPF方法的估计效果要优于其他两种方法.对基波、3次谐波幅值和相位估计与分析有相同结果.
动态谐波电压信号(幅值随时间变化的信号),可描述为
y(t)=[1.2+a1(t)]sin(ωt+75°)+
[0.7+a2(t)]sin(3ωt+55°)+
[0.3+a3(t)]sin(5ωt+45°)+μ(t)(31)
其中,ω=2πf; f=50 Hz; μ(t)是随机高斯噪声,其均值为0,方差为1.
{a1(t)=0.12sin(2πf1t)+0.07sin(2πf5t)
a2(t)=0.07sin(2πf3t)+0.03sin(2πf5t)
a3(t)=0.03sin(2πf1t)+0.05sin(2πf5t)(32)
其中, f1=2.0 Hz; f3=5.0 Hz; f5=7.0 Hz.
采用KF、UKF和UPF对该含有噪声的动态谐波基波、3次谐波和5次谐波幅值的估计和其分别对应的MSE值如图3至图5,本研究仅给出了系统中含3次和5次动态谐波时,采用本方法的谐波幅值估计,对于系统含更多种谐波分量时,本方法同样适用,仿真都能在2~3个周波时间(0.04~0.06 s)内较好估计出动态谐波信号幅值.限于篇幅从略.
表2为0.05 s时采用KF、UKF和UPF的动态幅值估计结果对比.由表2可见,基于UPF的方法估计的幅值更接近实际值,其MSE值也均小于其他两种方法.
分别加入20和40 dB信噪比的高斯噪声后,在0.05 s时基于KF、UKF、UPF的谐波估计方法对基波幅值估计的MSE值分别为1.7×10-3、1.4×10-3、0.4×10-3和1.5×10-3、1.11×10-3、0.22×10-3,此时理论幅值应为1.327 2 p.u.,基于UPF估计的MSE值小于另外两种方法.为体现本研究所提算法的优越性,对于非高斯噪声情况下,将本方法与基KF、基于UKF和文献[15]提出的基PSOUKF方法进行比较,加入非高斯噪声u(t)=Γ(0.3, 0.2), 4种方法对基波在0.05 s时刻的幅值估计对比如表3.
加入的非高斯噪声为伽马分布时,与加入高斯噪声相比,各方法估计的幅值(真实值为1.327 2 p.u.)精度降低,但从各种方法估计幅值的MSE值可以看出,本研究提出的基于UPF的谐波估计方法相比基于PSOUKF的谐波估计方法在非高斯情形下,对非高斯(伽马分布)的情形具有更强的处理能力,是各种方法中效果最好、精度最高的一种算法.对于其他非高斯噪声情况,我们也进行了仿真,结论相同,限于篇幅,不再列举.
动态谐波信号(相位随时间变化而变化)的信号可描述为
y(t)=1.2sin[ωt+a1(t)+75°]+
0.5sin[3ωt+a2(t)+55°]+
μ(t)(33)
其中, ω=2πf, f=50 Hz; μ(t)是随机高斯噪声,其均值为0,方差为1;
{a1(t)=0.12sin(2πf1t)+0.07sin(2πf5t)
a2(t)=0.07sin(2πf3t)+0.03sin(2πf5t)(34)
这里, f1=2.0 Hz; f3=5.0 Hz; f5=7.0 Hz.
表4为t=0.05 s时基波和3次谐波的动态相位估计效果对比,图6和图7为基波和3次谐波相位的估计波形图.
算法 相位/(°)基波 3次谐波 MES基波 3次谐波KF 75.010 3 55.012 2 0.006 3 0.111 2UKF 75.008 9 55.016 7 0.003 4 0.067 3UPF 75.006 4 55.009 2 0.000 5 0.020 7
由表4、图6和图7可见,对于谐波的相位估计,本算法也优于KF和UKF方法.加入非高斯噪声u(t)=Γ(0.3, 0.2)后, t=0.05 s时基波相位实际值应为75.127 2°,本研究所提出的基于UPF的谐波估计方法与文献[15]提出的基于PSOUKF的谐波估计方法对该时间基波相位的估计值分别为75.128 4°和75.095 1°,MSE值分别为0.001 1和0.001 4,相对误差分别为0.002%和0.043%,各项数据表明,在非高斯噪声下,基于UPF的估计效果要优于基于PSOUF的方法.
仿真结果说明,在同样的高斯噪声下,基于UPF的谐波估计方法对谐波各时刻幅值以及相位的估计所得到的MSE值均要小于其他两种方法,其效果要好于基于KF和UKF的方法,在非高斯噪声下,也表现出更强的估计能力.
本研究提出一种基于UPF的动态谐波幅值及相位的估计方法,仿真结果表明,本研究所提出的方法其性能总体优于基于KF和UKF的方法,对于静态和动态谐波以及在不同噪声下都具有良好的估计精度,对于非高斯噪声,所提方法较文中所提其他方法表现出更强的估计能力.目前,如何用UPF算法对电网频率以及间谐波的进行估计,如何选择最合适的重要密度函数和算法的收敛性研究,尚有待深入研究.
深圳大学学报理工版
JOURNAL OF SHENZHEN UNIVERSITY SCIENCE AND ENGINEERING
(1984年创刊 双月刊)
主 管 深圳大学
主 办 深圳大学
编辑出版 深圳大学学报理工版编辑部
主 编 阮双琛
国内发行 深圳市邮电局
国外发行 中国国际图书贸易集团有限公司(北京399信箱)
地 址 北京东黄城根北街16号
邮 编 100717
电 话 0755-26732266
0755-26538306
Email journal@szu.edu.cn
标准刊号 ISSN 1000-2618
CN 44-1401/N