基于无迹粒子滤波理论的发电机动态状态估计方法与流程

发明涉及一种基于无迹粒子滤波理论的发电机动态状态估计方法，属于电力系统监测、分析和控制技术领域。

背景技术：

随着电网规模的增大和电力系统复杂度增加，同时由于现有的测量设备存在测量误差，直接测量手段难以获取电力系统真实的状态信息。电力系统真实状态是电力系统分析、控制和决策的重要参考，而状态估计可以滤除电力系统量测数据的误差，得到尽可能接近系统真实状态的近似值。

传统的状态估计常采用以最小二乘算法为代表的静态状态估计方法迭代求解得到电力系统某一时间断面的状态近似值。由于电力系统规模大、复杂度高，电力系统暂态故障难以避免，因此，急需寻找精度和计算效率更高的状态估计方法从而更快得到更加准确的状态信息，缩短故障发生到采取控制保护的时间。相比于静态状态估计，动态状态估计不仅可以滤除量测噪声，且具有良好的预测能力，能够为电力系统的安全评估、状态预测、预防控制等在线功能提供支撑。而广域测量系统中的PMU设备能够提供高精度、高频刷新的量测数据，也为动态状态估计的实时性提供了有力保障。

建立合理的发电机动态模型和选择性能优良的滤波器是电力系统发电机动态状态估计的首要任务。针对发电机动态模型，不同学者根据研究需求选择了不同阶次的发电机动态方程建立发电机动态模型，主要分为二阶、四阶和六阶动态方程。而考虑到发电机动态模型的非线性，多位学者分别提出了基于卡尔曼估计和贝叶斯估计等框架下的滤波算法，主要代表有扩展卡尔曼滤波(EKF)、无迹卡尔曼滤波(UKF)、容积卡尔曼滤波(CKF)、粒子滤波(PF)等。EKF存在线性化截断误差大的问题，UKF及CKF采用不同的采样方式，利用采样点传播非线性方程的均值和方差，无需进行线性化。然而，UKF对初值敏感，参数无确定选取原则，CKF虽无需选择参数，但滤波提升有限。PF滤波精度高，但需要大量粒子进行运算，计算效率低，且当预测先验与似然函数重叠较少或量测模型精度较高时，可能导致滤波失效。

技术实现要素：
：

本发明所要解决的技术问题是针对现有技术存在的不足而提供一种基于无迹粒子滤波理论的发电机动态状态估计方法。

本发明为实现上述目的，采用如下技术方案：

一种基于无迹粒子滤波理论的发电机动态状态估计方法，所述方法是在计算机中依次按以下步骤实现的：

1)获取所需动态状态估计的系统中发电机组的参数信息；

2)利用电力系统分析软件模拟PMU设备获取状态估计所需的量测数据；

3)状态估计器初始化；

4)建立发电机动态状态估计模型；采用发电机的四阶动态方程建立发电机的状态方程，根据获取的PMU数据建立发电机的量测方程；

5)在状态初值附近生成初始粒子，启动滤波算法；

6)通过比例修正采样的UKF生成重要性密度函数及新的采样粒子；

7)更新粒子权重并归一化；

8)判断是否需要重采样；若是，则转到重采样步骤，若否，则继续下一步。

9)输出当前时刻的状态估计结果。

10)判断滤波运算是否结束，若是，则输出最后全部状态估计结果；若否，则转到步骤6)继续下一步。

步骤1)中，发电机的参数信息包括发电机的机械转矩、有功和无功额定功率、同步电角速度、阻尼系数、惯性时间常数、定子励磁电压额定值以及发电机的机组总数。

步骤2)中的量测数据包括：发电机的绝对功角、角速度变化值、端口电压相量的幅值和相角、端口电磁功率和机械转矩。

步骤3)中的初始化包括输入参数和量测数据，设置状态初值，设置过程噪声和量测噪声协方差阵，设置预测协方差初值，设置UPF粒子数和滤波参数以及采样间隔和采样周期；

步骤4)中发电机的状态方程形式为：

式中，δ为dq0坐标系下发电机转子的绝对功角(弧度)，ω₀为额定同步转速(电弧度/秒)，Δω为发电机角速度变化量(标幺值)，TJ为发电机时间惯性常数,D为发电机阻尼系数，E_f为定子励磁电压，T_m为发电机的机械转矩(标幺值)，T_e(T_e＝P_e/ω)为发电机的电磁转矩(标幺值)，忽略定子阻抗并假设ω≈1，得T_e≈P_e(P_e为发电机的电磁功率)，E′_d和E′_q分别为发电机的d轴和q轴暂态电动势，U和为发电机端口电压的幅值和相角，X_d和X_q分别为发电机d轴和q轴的同步电抗，X′_d和X′_q分别为发电机d轴和q轴的暂态电抗，T′_d0和T′_q0分别为发电机d轴和q轴的暂态开路时间常数；

发电机量测量分别为发电机绝对功角、角速度、电磁功率，由PMU设备直接获得，即：y＝[δ_z,ω_z,P_ez]^T＝[y₁,y₂,y₃]^T，发电机量测方程形式为：

步骤6)中利用比例采样的UKF生成重要性密度函数包括:

比例修正采样的权值满足：

式中：λ＝α²(L+κ)-L是微调参数，用来控制点到均值的距离，其中L是状态变量个数，κ是次级采样因子，取为3-L；α是比例修正因子，取为0.12；β是待选参数，调节β提高方差精度，取β＝2。

由UKF生成重要性密度函数的步骤包括：

由上一时刻每个粒子的状态滤波值和滤波协方差矩阵生成新的采样点，代入发电机的状态方程，获取每个粒子在UKF中的第一步状态预测均值和第一步预测协方差矩阵；

由每个粒子的第一步预测值生成新的采样点，利用量测方程和当前时刻的量测信息修正每个粒子在UKF中的第二步状态预测均值和第二步预测协方差矩阵；

根据第二步状态预测均值和第二步预测协方差矩阵生成PF中的重要性密度函数并生成新的采样粒子；

步骤8)中的判断重采样的有效粒子阈值选为总粒子数的三分之一，重采样方法选为计算量较小，性能良好的残差重采样。

有益效果：本发明与现有技术相比：能够对机电暂态过程中的发电机状态量提供更准确的估计，提高了对发电机状态量滤波的准确度。本发明提出的方法减小了对粒子的需求量，计算效率优于传统的PF估计器，同时增加了粒子的分散性，使得本发明提出的方法的鲁棒性也优于传统的PF估计器，且本发明提出的方法可根据粒子数调整滤波精度，滤波精度和灵活性均优于UKF估计器。

附图说明：

图1为本发明的方法流程图；

图2为WSCC三机九节点系统示意图；

图3为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机1的功角估计曲线与BPA仿真软件真实值曲线的对比图；

图4为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机1的角速度估计曲线与BPA仿真软件真实值曲线的对比图；

图5为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机1的q轴暂态电动势估计曲线与BPA仿真软件真实值曲线的对比图；

图6为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机1的d轴暂态电动势估计曲线与BPA仿真软件真实值曲线的对比图；

图7为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机2的功角估计曲线与BPA仿真软件真实值曲线的对比图；

图8为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机2的角速度估计曲线与BPA仿真软件真实值曲线的对比图；

图9为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机2的q轴暂态电动势估计曲线与BPA仿真软件真实值曲线的对比图；

图10为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机2的d轴暂态电动势估计曲线与BPA仿真软件真实值曲线的对比图；

图11为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机3的功角估计曲线与BPA仿真软件真实值曲线的对比图；

图12为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机3的角速度估计曲线与BPA仿真软件真实值曲线的对比图；

图13为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机3的q轴暂态电动势估计曲线与BPA仿真软件真实值曲线的对比图；

图14为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机3的d轴暂态电动势估计曲线与BPA仿真软件真实值曲线的对比图；

图15为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机1的功角估计曲线与UKF、PF算法的对比示意图；

图16为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机1的角速度估计曲线与UKF、PF算法的对比示意图；

图17为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机1的q轴暂态电动势估计曲线与UKF、PF算法的对比示意图；

图18为以WSCC三机九节点系统为基础测试系统，本发明方法下发电机1的d轴暂态电动势估计曲线与UKF、PF算法的对比示意图；

图19为国内某实际大电网系统为有偏噪声测试系统，本发明方法在四种不同的噪声情况下对发电机的功角估计曲线与BPA仿真软件真实值曲线的对比示意图；

图20为国内某实际大电网系统为有偏噪声测试系统，本发明方法在四种不同的噪声情况下对发电机的角速度估计曲线与BPA仿真软件真实值曲线的对比示意图；

图21为国内某实际大电网系统为有偏噪声测试系统，本发明方法在四种不同的噪声情况下对发电机的q轴暂态电动势估计曲线与BPA仿真软件真实值曲线的对比示意图；

图22为国内某实际大电网系统为有偏噪声测试系统，本发明方法在四种不同的噪声情况下对发电机的d轴暂态电动势估计曲线与BPA仿真软件真实值曲线的对比示意图。

具体实施方式：

下面结合附图对发明的技术流程进行详细说明：

动态状态估计

动态状态估计首先采用电力系统分析软件BPA模拟PMU设备进行数据采集，并将BPA运行结果经过转换与计算作为真实值，将真实值叠加随机误差作为量测值送入状态估计器。

近年来，PMU设备的广泛运用给电力系统动态状态估计带来变革。PMU设备可提供高精度，高频刷新的量测数据，为动态状态估计的实时性提供了保障。由于发电机转子存在惯性，发电机的动态过程中功角等状态量不会发生突变。已知发电机前一时刻状态量的滤波值，利用状态方程得到当前时刻预测值，结合PMU量测值，可得当前时刻状态量的滤波值，以此实现时间轴上的递推，最终得到发电机的动态状态估计结果。发电机的状态方程和量测方程的形式为：

式中，x、u、y分别定义为状态变量、控制变量和量测变量，f是发电机状态方程，h是发电机量测方程，w是过程噪声，v是量测噪声，一般假设w～N(0,Q)，v～N(0,R)，其中Q和R分别为w和v的误差方差阵，w和v相互独立且独立于状态变量。

发电机状态方程和量测方程的离散化形式为：

式中，F(x_k,u_k,w_k)＝x_k+f(x_k,u_k,w_k)Δt，k为采样时刻，Δt为采样步长。

发电机的模型

发电机的动态状态估计模型由状态方程和量测方程组成，发电机状态方程的个数即发电机的阶数。由于发电机次暂态过程短暂，现有PMU设备难以准确获取次暂态过程量测量，忽略与次暂态过程对应的D、Q绕组及定子动态过程，将发电机模型降维，得到同步发电机四阶非线性模型。利用PMU设备测得发电机出口的电压、电流等相量将发电机与其连接的外部电网进行解耦，使发电机状态能够独立于电网进行动态状态估计。发电机动态状态估计的状态方程如下：

式中，δ为dq0坐标系下发电机转子的绝对功角(弧度)，ω₀为额定同步转速(电弧度/秒)，Δω为发电机角速度变化量(标幺值)，TJ为发电机时间惯性常数,D为发电机阻尼系数，E_f为定子励磁电压，T_m为发电机的机械转矩(标幺值)，T_e(T_e＝P_e/ω)为发电机的电磁转矩(标幺值)，忽略定子阻抗并假设ω≈1，得T_e≈P_e(P_e为发电机的电磁功率)，E′_d和E′_q分别为发电机的d轴和q轴暂态电动势，U和为发电机端口电压的幅值和相角，X_d和X_q分别为发电机d轴和q轴的同步电抗，X′_d和X′_q分别为发电机d轴和q轴的暂态电抗，T′_d0和T′_q0分别为发电机d轴和q轴的暂态开路时间常数；

量测量y＝[δ_z,ω_z,P_ez]^T＝[y₁,y₂,y₃]^T，动态状态估计的量测方程如下：

本发明方法的步骤

利用UPF算法实施动态状态估计。UPF基于PF框架，利用UKF生成重要性密度函数代替PF中的重要性密度函数，算法具体步骤如下：

步骤一：初始化：

输入发电机参数及两侧数据等信息，程序初始化，启动状态估计器，令k＝0，在初始状态变量x₀附近生成原始粒子集，相应的权值系数(M为粒子个数)。

当k>0时，令k＝k+1，将上步迭代中得到的一组新的粒子集代入下一步。步骤二：利用UKF生成重要性密度函数及采样粒子：

设和是每个粒子k时刻的状态量的均值和协方差，采用sigma点比例修正采样方式，构造权值和sigma采样点：

式中：λ＝α²(L+κ)-L是微调参数，用来控制点到均值的距离，其中L是状态变量个数，κ是次级采样因子，取为3-L；α是比例修正因子，取为0.12；β是待选参数，调节β提高方差精度，β＝2。

式中，是矩阵的第i行(P＝A^TA,A是P通过柯林斯基分解得到的低阶三角矩阵)或i列(P＝AA^T)。

对每个粒子进行预测，计算状态量一步预测值及其预测协方差

根据和再次构造每个粒子sigma采样点：

对每个粒子状态和方差更新：

至此，对每个粒子利用UKF算法得到当前时刻该粒子的估计值和协方差从而得到重要性密度函数为：

其中N(·)是服从高斯分布的概率密度函数。

则重要密度函数中的采样粒子为：

步骤三：计算权值并归一化：

计算权值：

归一化：

步骤四：判断是否需要进行重采样，若是，则进入重采样步骤，若否进入下一步：

计算若(N_t为设定的阈值，一般取为M/3)则说明粒子的权值已经退化严重，需要进行重采样，否则直接进入下一步。重采样时，按权值高低复制与重排粒子，复制的个数与其权值成正比，剔除权值小的粒子，并将处理后的粒子映射为等权重的M个粒子，即

步骤五：输出状态估计值和协方差矩阵P_k+1：

步骤六：判断循环是否结束，若是，则输出全部结果，滤波结束，若否，转入步骤二。本发明方法的优点

由于EKF截断误差大，而UKF利用采样点处理非线性函数，虽然保留了函数的非线性，滤波精度高，但UKF初值选取不准确可能导致收敛缓慢甚至不收敛，且UKF参数的选取缺乏严格的理论证明，其滤波效果受参数制约，无法根据滤波场合的不同要求调整滤波精度，缺乏灵活性。基于贝叶斯理论的PF可改变粒子数调整滤波精度，且随着粒子数的增加，估计值能无限接近真实值。但PF使用大量的粒子对状态空间进行搜索，加大了计算负担，且PF以先验概率为建议分布，当似然函数恰好在先验概率分布的尾部或者似然函数太狭窄(由于量测误差很小导致的尖峰)时，滤波就容易发散或失效。

因此，本发明方法做了以下改进：

首先，本发明方法采用PMU设备获取发电机量测数据，实现了发电机与外部电网解耦，减少传输误差影响。其次，UPF致力于解决卡尔曼框架下的滤波方法对强非线性非高斯系统的滤波精度有限，而传统PF在采样过程中由于重要性密度函数不准确容易出现粒子退化的问题。本发明方法利用UKF生成PF中的重要性密度函数，它具备了最优重要性密度函数产生预测粒子的两大特点：一是已有当前时刻的粒子信息，二是利用最新的量测数据。因此，本发明方法利用UKF产生的重要性密度函数与真实后验概率密度函数支集重叠更多，由于在生成重要性密度函数的过程中，本发明方法使用了最新的量测信息，可将粒子转移至高似然区域，减少了描述后验概率密度的粒子需求量，提高了滤波精度和滤波效率。此外，本发明方法易于采样实现，又引入重采样，增加了有效粒子的权重，估计精度高，能准确跟踪状态量的变化，且对噪声具有鲁棒性。

实施例

本发明测试的算例为WSCC三机九节点系统和某224节点海南电网系统。WSCC三机九节点系统节点量测数据为BPA仿真软件真值叠加随机噪声模拟获取，仿真时采用发电机四阶模型，并考虑发电机2水轮调速器作用，获得发电机2机械转矩的真实值曲线。仿真故障设置为第40周波至第45周波(即0.8～0.9s)，Bus5-Bus7支路线路出口处发生三相金属性短路，而后故障清除。发电机经过一段时间后过渡到新的稳态。仿真时长为6s。

为了验证本发明方法的性能，选择跟踪速度较快的UKF和滤波精度较高的PF作为本发明的方法的对比方法，将三种算法的性能进行了对比。

为了更加直观的比较三种算法的滤波效果，选取平均相对误差和最大绝对误差x_m作为算法对比的性能评价指标。即：

式中，N为采样的次数，取为300，xi为第i次采样的真实值，为第i次滤波的估计值。

图3～图5为WSCC三机九节点系统中三台发电机的状态量估计值和BPA仿真结果对比图。可见在0～40周波(系统稳态时)，UPF算法能够快速收敛，实时跟踪，当系统突然发生暂态故障，UPF仍能精确跟踪故障时和故障后的发电机的各个状态变量，估计精度满足要求。

图7～图10为发电机1的四个状态变量在本发明方法和UKF以及PF算法下的滤波曲线。其中UKF算法中，其滤波协方差初值设为单位阵，以发电机稳态时的初值作为状态变量初值。UKF中κ＝3-L，α＝0.12，β＝2。UPF中参数与协方差设置与UKF相同，本次测试中UPF与PF均采用残差重采样法，粒子数均为100。可见，UPF算法在滤波启动后能够快速收敛，在故障发生前后均能准确跟踪发电机状态量，滤波精度高。UKF算法在滤波启动后收敛较慢，故障发生后不能准确跟踪，滤波精度低于UPF。PF算法在滤波启动后能较快收敛，在故障时刻也能实时跟踪，但滤波精度低于UPF。

比较表1中三种算法的指标数据可知，UPF在平均相对误差和最大绝对误差指标上滤波效果均明显优于UKF和PF。UPF滤波精度高，能准确跟踪发电机状态变量，在同样的参数设置下滤波精度优于UKF，在同样的粒子数下滤波精度也明显优于PF。

为研究不同粒子数对本发明方法滤波效果的影响，在实际电网系统中选取编号为HAID4G的发电机状态进行跟踪和滤波，研究在不同粒子数下的本发明方法滤波性能。由表2可见，本发明方法的滤波精度随着粒子数的增加呈变优趋势，但其在提高滤波精度的同时存在着计算量不断增大的问题。本发明方法可以根据不同的滤波精度要求选择不同的粒子数，即在滤波精度相对较低的情况下选取少量粒子来减少计算负担，在滤波精度较高而对计算时间要求不大的场合可以适当增加粒子数保证精度，灵活性强。PF在达到本发明方法相同滤波精度时需要用到更多的粒子，计算量也相应的增加，计算负担变重。测试表明，相比于PF，在达到同样精度要求的情况下，本发明方法的计算量约为PF的一半，比PF效率高。

为了验证本发明方法在方差增大的有偏噪声下的滤波效果，在海南某实际电网中选取编号为HAID4G的发电机状态进行跟踪和滤波，并将UKF与PF在相同噪声条件下的滤波性能与本发明方法进行了对比。本次测试中UPF与PF粒子数均为100。在50周波时设置某发电机出口线路三相短路故障。量测噪声设置为逐步增大的高斯白噪声和均值非零的高斯噪声，为以下4种情况：Case1，均值为0，标准差为δ₀；Case2，均值为0，标准差为3δ₀；Case3，均值为E₀，标准差为δ₀；Case4，均值为3E₀，标准差为3δ₀。

由表3可见，随着噪声标准差增大或噪声均值偏离零值尺度增加，三种滤波方法下发电机四个状态变量的平均相对估计误差逐步增大，其中，在前三种情形时呈现UPF滤波效果优于PF，PF滤波效果优于UKF的情形，而当噪声均值偏离零值较远，误差标准差也增大到原标准差的3倍时，100个粒子下的PF容易出现滤波不能准确跟踪，滤波部分曲线成常数曲线，滤波精度急剧变差，此时滤波失效。这是由于发电机发生故障期间的过程噪声尤其是故障初始时段的过程噪声远大于稳定期间及故障消失后的过程噪声。PF在故障初始时段低估了预测步的过程噪声方差，导致先验状态量过程噪声的协方差小于真实过程噪声协方差，粒子的分散性弱，在噪声有偏和过程噪声偏差的双重影响下，100个粒子下的PF滤波效果急剧变坏。通过增大PF的过程噪声协方差并增加PF的粒子数可以相应改善PF的滤波效果。UPF通过选择UKF生成重要性密度函数利用了最新的量测信息，滤波精度更高，且在UPF滤波过程中，对每个粒子均进行sigma点采样，增加了粒子的分散性，因此在故障过程中受过程噪声的影响较小。通过三种方法的对比发现，UPF的滤波精度在标准差增大或噪声均值非零等噪声情况下仍能满足滤波精度和收敛性要求，对噪声具有较强的鲁棒性。

表1不同算法下的发电机动态状态估计评价指标

表2不同粒子数下发电机动态状态估计评价指标

表3不同噪声下3种算法动态状态估计评价指标