CN110162871A - 一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法 - Google Patents

一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法 Download PDF

Info

Publication number
CN110162871A
CN110162871A CN201910409884.9A CN201910409884A CN110162871A CN 110162871 A CN110162871 A CN 110162871A CN 201910409884 A CN201910409884 A CN 201910409884A CN 110162871 A CN110162871 A CN 110162871A
Authority
CN
China
Prior art keywords
moment
particle
weight
value
population
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
CN201910409884.9A
Other languages
English (en)
Other versions
CN110162871B (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.)
Northeastern University China
Original Assignee
Northeastern University China
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 Northeastern University China filed Critical Northeastern University China
Priority to CN201910409884.9A priority Critical patent/CN110162871B/zh
Publication of CN110162871A publication Critical patent/CN110162871A/zh
Application granted granted Critical
Publication of CN110162871B publication Critical patent/CN110162871B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/06—Energy or water supply
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Theoretical Computer Science (AREA)
  • Economics (AREA)
  • Physics & Mathematics (AREA)
  • Strategic Management (AREA)
  • General Physics & Mathematics (AREA)
  • Tourism & Hospitality (AREA)
  • Health & Medical Sciences (AREA)
  • Marketing (AREA)
  • General Business, Economics & Management (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Water Supply & Treatment (AREA)
  • Game Theory and Decision Science (AREA)
  • Educational Administration (AREA)
  • Development Economics (AREA)
  • Operations Research (AREA)
  • Public Health (AREA)
  • Quality & Reliability (AREA)
  • General Health & Medical Sciences (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Primary Health Care (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明的一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法，涉及电力系统监测、分析和控制技术领域。步骤为：步骤1：将当前电力系统状态变量进行初始化，计算初始状态均值、方差；步骤2：令k＝1；步骤3：重要性采样得到k时刻的粒子群和权重；步骤4：对粒子群进行粒子群分裂和权值调整；步骤5:判断Neff<Nth是否成立，是转到步骤6，否转到步骤8；步骤6：对粒子群进行复制和淘汰得到新的粒子群和权重；步骤7：重新进行权重归一化；步骤8：计算电力系统的状态估计值；步骤9：判断k≥Ω是否成立，不成立令k＝k+1转至步骤3，成立结束电力系统的动态估计。本方法改善了建议密度分布问题，可有效提高估计精度，有效解决粒子匮乏问题。

Description

一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法
技术领域
本发明涉及电力系统监测、分析和控制技术领域，尤其涉及一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法。
背景技术
随着社会的不断发展，人们对电力系统稳定性的要求越来越高。为了提高电力系统的稳定性，我们必须提高电力系统的调度、控制、安全评估等方面。电力系统动态估计是电力系统的调度、控制、安全评估等方面的基础，动态状态估计通过状态预测能对电力系统进行安全评估,实现经济调度、预防控制等在线功能,重要性不言而喻,鉴于此，我们必须快速、准确的进行电力系统的动态估计。
上世纪70年代初，Debs等人提出了卡尔曼滤波算法，并用最简单的系统模型建立动态状态估计，卡尔曼滤波处理的是线性问题，而电力系统的动态估计属于非线性问题，此方法的状态估计精度低。然后引出了扩展的卡尔曼滤波算法(EKF)，扩展的卡尔曼滤波算法处理非线性系统的思想是将非线性函数在估计点附近进行泰勒级数展开，并用一个等价于常规卡尔曼滤波方程的近似矩阵来代替非线性函数。此方法的估计精度较低，然后又引入了无迹卡尔曼滤波算法(UKF))来进行状态估计，无迹卡尔曼滤波也是地柜最小方差估计器，核心思想是无迹变换。如将非线性方程采用泰勒级数展开式表达，可以看出无迹卡尔曼滤波算法能够精确到与三阶泰勒级数展开式相当的均值和方差。
但EKF和UKF都是针对非线性系统的线性卡尔曼滤波方法的变形和改进。因此受线性卡尔曼滤波算法的条件制约。粒子滤波(PF)能够很好的进行非线性系统的状态估计，在一定程度上提高了状态估计精度，但是PF选择转移概率密度作为重要性密度函数，在计算重要性概率密度函数时没有考虑到最新的量测信息，当预测先验与似然函数重叠少或者量测模型密度较高时，可能偏离真实的后验分布，导致粒子滤波存在着退化现象和粒子匮乏问题。
发明内容
本发明要解决的技术问题是针对上述现有技术的不足，提供一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法，本方法改善了建议密度分布问题，可以有效的提高估计精度，有效的解决了粒子匮乏问题。
为解决上述技术问题，本发明所采取的技术方案是：
本发明提供一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法，包括如下步骤：
步骤1：将当前电力系统状态变量进行初始化，获取k＝0时刻的粒子群Φ,其中k代表时间；粒子个数为N,Fori＝1：N,从先验概率密度函数p(X0)中抽取初始化状态作为电力系统的初始状态，计算初始状态均值方差
其中X0＝[v0θ0]，X0为0时刻电力系统的初始状态，v0表示0时刻的节点电压的幅值，θ0表示0时刻的节点电压的相角，为0时刻第i个粒子的初始值，P0为初始状态估计误差的协方差，Q为过程噪声的协方差，R为测量噪声的协方差，E(·)代表求期望值的函数，(·)T代表(·)的转置；
步骤2：令k＝1；
步骤3：计算k-1时刻下粒子群中所有粒子的Sigma采样点集合，根据该集合预测出k时刻的Sigma采样点集合，并通过k时刻的Sigma采样点集合得出k时刻的状态粒子群Φk和粒子权重wi(k)；
步骤3.1：计算k-1时刻下粒子i的Sigma采样点集合，遍历k-1时刻的粒子群，得到所有粒子的Sigma采样点集合；
其中为k-1时刻第i个粒子的Sigma采样点，维数为n，为k-1时刻第i个粒子的方差，为k-1时刻第i个粒子的状态均值,λ为采样因子，λ＝αk-1 2(n+b)-n，b为自由参数，用来捕捉给定的高阶矩阵的信息，为矩阵的第j列，α为比例修正因子，取值为[0.0001,1]，设定初始值为α0＝0.01，其中α1到αt的计算步骤如下：
步骤3.1.1：令则上式可以转化为：其中ck-1表示k-1时刻的均值到sigma采样点的距离；
步骤3.1.2：计算
步骤3.1.3：令其中I为sigma采样点集合，为ck-1的最大值；
步骤3.1.4：计算k刻的比例修正因子其中tr(·)代表求矩阵的迹；
步骤3.2：根据步骤3.1中得到k-1时刻下粒子群中所有粒子的Sigma采样点集合预测出k时刻中所有粒子的Sigma采样点集合；
k时刻粒子i Sigma采样点集合的公式如下：
其中为第i个粒子k-1时刻对k时刻的预测值，f(·)为预测函数，为第i个粒子k-1时刻对k时刻的预测值的均值，为第i个粒子k-1时刻对k时刻的第a个预测值，为第i个粒子k-1时刻对k时刻预测值的方差，为第i个粒子k-1时刻对k时刻预测的量测值，为第i个粒子k-1时刻对k时刻预测的量测值的均值，h(·)为量测函数，为第i个粒子k-1时刻对k时刻预测的第a个量测值，代表sigma采样点求期望的权值、代表sigma采样点求方差的权值；
其中计算如下：
其中，a＝1,…,2n，为sigma采样点求方差的初始权值；为sigma采样点求期望的初始权值；α0为比例修正因子设定的初始值，β为高阶误差采样因子；
步骤3.3：将步骤3.2中得到的和融入到k时刻已知的实际量测值中并对k时刻粒子i的预测值的均值和k时刻粒子i的预测协方差进行更新；重复本步骤遍历k-1时刻下粒子群中所有粒子，得到k时刻所有粒子的预测值的均值和k时刻所有粒子的预测协方差
其中为已知量测zk时的方差，为已知量测zk时Xk的方差，Xk为电力系统k时刻的状态值，Kk为计算增益，Zk为k时刻的量测值；
步骤3.4：根据步骤3.3求出k时刻的采样更新粒子最终得到k时刻的粒子群Φk；
其中代表以为均值、以为方差的正态分布函数；
步骤3.5：通过PMU的量测和电力系统的静态估计形成混合量测，计算粒子的权重；其中PMU的量测和电力系统的静态估计形成混合量测具体步骤如下：
步骤3.5.1：PMU的量测是通过安装PMU的节点获取该节点的电压幅值和相角和电流；通过计算得到与其相连节点的电压的伪量测值，根据该值最粒子的权重进行更新；如图2所示,具体计算如下：
其中Y10为节点对地导纳，Y12为节点1、2直接互导纳，为流向节点1的电流，为节点1的电压，为节点2的电压；
步骤3.5.2：获取已有的电力系统的静态估计结果，选取静态估计精度高的结果用来对粒子的权重进行更新；
其中,d＝1,2,…M,代表静态估计值，V代表电压幅值，θ代表电压相角；代表安装PMU节点的量测值，通过向量节点计算得到的伪量测值,为k时刻第i个粒子第d个量测的初始权重，Zd(k)代表k时刻的第d个量测值，为k时刻第i个粒子第d个量测的初始权重的权重系数，σd为第d个量测的量测噪声，M代表量测的总数；
步骤3.5.3：最终权重为：
步骤3.6：将权重归一化，得到归一化权重wi(k)；归一化方法如下：
步骤4：对k时刻的粒子群进行粒子群分裂和权值调整；将粒子的权值按照高权值矩阵和低权值矩阵排列，将高权值矩阵中的后1/2的粒子分裂成两个权重减少1/2的粒子，把低权值的后1/2的粒子去除,形成新的粒子群Φ′k，重新进行权重归一化，得到权重w′i(k)；
步骤5:计算有效粒子Neff，定义阈值判断Neff<Nth是否成立，如果成立转到步骤6，如果不成立转到步骤8；其中Neff的计算如下：
步骤6：利用随机重采样方法，根据归一化权值w′i(k)大小，对粒子群Φ′k进行复制和淘汰得到粒子群For i＝1:N,重新设置权值权重
步骤7：利用Markov链蒙特卡罗方法保持重采样后粒子的多样性，对粒子群Φ″k进行更新,得到新的粒子群Φ″′k将权重wi″(k)重新进行权重归一化，得到权重w″′i(k)；
步骤8：计算得到电力系统的状态估计值，所述状态估计值包括当前电力系统电压的幅值和相角；判断是否为步骤4跳转至本步骤，若是，则电力系统的状态估计值的公式为：若否，则电力系统的状态估计值的公式为
步骤9：判断k≥Ω是否成立，如果不成立，令k＝k+1转至步骤2，如果成立，结束电力系统的动态估计，其中Ω为预测截止时间。
所述步骤6中随机重采样方法的具体步骤如下：
步骤6.1:产生[0,1]均匀分布的随机数组其中N为粒子数，g代表代表生成的随机数，并按升序排序；
步骤6.2:产生粒子权累积函数cdf,满足
步骤6.3:复制权值大的粒子，具体计算如下：
步骤6.3.1：令g＝1，i＝1；
步骤6.3.2：判断cdf(i)≤ug是否成立，成立则转至下一步，否则转至6.3.4；
步骤6.3.3：i＝i+1,返回步骤6.3.2；
步骤6.3.4：第i个粒子复制在第g个位置；
步骤6.3.4：g＝g+1；
步骤6.3.5：判断g≤N是否满足，满足转至步骤6.3.2，不满足则终止。
所述步骤7的具体步骤如下：
步骤7.1：获取已经得到的
步骤7.2：令i＝0；
步骤7.3:从[0,1]中抽样得到随机数e；
步骤7.4:根据建议分布和当前变量值采样得到候选值计算接收概率
其中为的目标分布，为的建议分布；
步骤7.5:判断是否成立，如果成立，则输出否则，输出
步骤7.6:i＝i+1，判断i<N是否成立，成立转至步骤7.3；否则，将步骤7.5中输出的粒子形成新的粒子群Φ″′k，转至7.7；
步骤7.7:对粒子的权重w″i(k)进行归一化，得到w″′i(k)。
采用上述技术方案所产生的有益效果在于：本发明提供的一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法，本方法在粒子滤波中引入了无迹卡尔曼滤波，改善了建议密度分布问题，减少了粒子滤波的退化，在粒子滤波中加入了粒子的分裂和权值更新，可以有效的提高估计精度，在粒子滤波中引入了重采样，有效的解决了粒子匮乏问题。由于无迹卡尔曼滤波的引入会增加估计时间，在保障精度的情况下，适当减少粒子数可以缩短状态估计的时间。
附图说明
图1为本发明实施例提供的电力系统动态状态估计的方法流程图；
图2为本发明实施例提供的PMU计算伪量测的图，其中1代表节点1,2代表节点2；
图3为本发明实施例提供的IEEE-3节点的拓扑图，其中，1代表发电机1，2代表代表发电机2，3代表母线3，4代表母线2，5代表连接线1-3，6代表连接线2-3，7代表连接线1-2，8代表母线1；
图4为本发明实施例提供的中节点1电压幅值动态估计仿真图；
图5为本发明实施例提供的电力系统动态电压幅值均方根误差对比图；
图6为本发明实施例提供的节点1电压相角动态估计仿真图；
图7为本发明实施例提供的电力系统动态电压相角均方根误差对比图；
图8为本发明实施例提供的电力系统动态估计时间对比图；
具体实施方式
下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
如图1所示，本实施例的方法如下所述。
本发明提供一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法，包括如下步骤：
步骤1：将当前电力系统状态变量进行初始化，获取k＝0时刻的粒子群Φ,其中k代表时间；粒子个数为N,Fori＝1：N,从先验概率密度函数p(X0)中抽取初始化状态作为电力系统的初始状态，计算初始状态均值方差
其中X0＝[v0θ0]，X0为0时刻电力系统的初始状态，v0表示0时刻的节点电压的幅值，θ0表示0时刻的节点电压的相角，为0时刻第i个粒子的初始值，P0为初始状态估计误差的协方差，Q为过程噪声的协方差，R为测量噪声的协方差，E(·)代表求期望值的函数，(·)T代表(·)的转置。
步骤2：令k＝1；
步骤3：计算k-1时刻下粒子群中所有粒子的Sigma采样点集合，根据该集合预测出k时刻的Sigma采样点集合，并通过k时刻的Sigma采样点集合得出k时刻的状态粒子群Φk和粒子权重wi(k)；
步骤3.1：计算k-1时刻下粒子i的Sigma采样点集合，遍历k-1时刻的粒子群，得到所有粒子的Sigma采样点集合；
其中为k-1时刻第i个粒子的Sigma采样点，维数为n，为k-1时刻第i个粒子的方差，为k-1时刻第i个粒子的状态均值,λ为采样因子，λ＝αk-1 2(n+b)-n，b为自由参数，用来捕捉给定的高阶矩阵的信息，在这里b＝0；为矩阵的第j列，α为比例修正因子，取值为[0.0001,1]，设定初始值为α0＝0.01，其中α1到αt的计算步骤如下：
步骤3.1.1：令则上式可以转化为：其中ck-1表示k-1时刻的均值到sigma采样点的距离；
步骤3.1.2：计算
步骤3.1.3：令其中I为sigma采样点集合，为ck-1的最大值；
步骤3.1.4：计算k刻的比例修正因子其中tr(·)代表求矩阵的迹；
步骤3.2：根据步骤3.1中得到k-1时刻下粒子群中所有粒子的Sigma采样点集合预测出k时刻中所有粒子的Sigma采样点集合；
k时刻粒子i Sigma采样点集合的公式如下：
其中为第i个粒子k-1时刻对k时刻的预测值，f(·)为预测函数，为第i个粒子k-1时刻对k时刻的预测值的均值，为第i个粒子k-1时刻对k时刻的第a个预测值，(Sigma点集合的个数是2n+1,a在这里是Sigma点集的第a个)，为第i个粒子k-1时刻对k时刻预测值的方差，为第i个粒子k-1时刻对k时刻预测的量测值，为第i个粒子k-1时刻对k时刻预测的量测值的均值，h(·)为量测函数，为第i个粒子k-1时刻对k时刻预测的第a个量测值，代表sigma采样点求期望的权值、代表sigma采样点求方差的权值；
其中计算如下：
其中，a＝1,…,2n，为sigma采样点求方差的初始权值；为sigma采样点求期望的初始权值；α0为比例修正因子设定的初始值，β为高阶误差采样因子；
步骤3.3：将步骤3.2中得到的和融入到k时刻已知的实际量测值中并对k时刻粒子i的预测值的均值和k时刻粒子i的预测协方差进行更新；重复本步骤遍历k-1时刻下粒子群中所有粒子，得到k时刻所有粒子的预测值的均值和k时刻所有粒子的预测协方差
其中为已知量测zk时的方差，为已知量测zk时Xk的方差，Xk为电力系统k时刻的状态值，Kk为计算增益，Zk为k时刻的量测值；
步骤3.4：根据步骤3.3求出k时刻的采样更新粒子最终得到k时刻的粒子群Φk；
其中代表以为均值、以为方差的正态分布函数；
步骤3.5：通过PMU的量测和电力系统的静态估计形成混合量测，计算粒子的权重；其中PMU的量测和电力系统的静态估计形成混合量测具体步骤如下：
步骤3.5.1：PMU的量测是通过安装PMU的节点获取该节点的电压幅值和相角和电流；通过计算得到与其相连节点的电压的伪量测值，根据该值最粒子的权重进行更新；具体计算如下：
其中Y10为节点对地导纳，Y12为节点1、2直接互导纳，为流向节点1的电流，为节点1的电压，为节点2的电压；
步骤3.5.2：获取已有的电力系统的静态估计结果，选取静态估计精度高的结果用来对粒子的权重进行更新；
其中,d＝1,2,…M,代表静态估计值，V代表电压幅值，θ代表电压相角。代表安装PMU节点的量测值，通过向量节点计算得到的伪量测值,为k时刻第i个粒子第d个量测的初始权重，Zd(k)代表k时刻的第d个量测值，为k时刻第i个粒子第d个量测的初始权重的权重系数，σd为第d个量测的量测噪声，M代表量测的总数；
步骤3.5.3：最终权重为：
步骤3.6：将权重归一化，得到归一化权重wi(k)；归一化方法如下：
步骤4：对k时刻的粒子群进行粒子群分裂和权值调整；将粒子的权值按照较高权值矩阵和较低权值矩阵排列，将较高权值矩阵中的后1/2的粒子分裂成两个权重减少1/2的粒子，把较低权值的后1/2的粒子去除，形成新的粒子群Φ′k，N的大小没有改变，按照步骤2.6的方法，重新进行权重归一化，得到权重w′i(k)；
步骤5:计算有效粒子Neff，定义阈值判断Neff<Nth是否成立，如果成立转到步骤6，如果不成立转到步骤8；其中Neff的计算如下：
步骤6：利用随机重采样方法，根据归一化权值w′i(k)大小，对粒子群Φ′k进行复制和淘汰得到粒子群For i＝1:N,重新设置权值权重所述随机重采样方法的具体步骤如下：
步骤6.1:产生[0,1]均匀分布的随机数组其中N为粒子数，g代表代表生成的随机数，并按升序排序；
步骤6.2:产生粒子权累积函数cdf,满足
步骤6.3:复制权值大的粒子，具体计算如下：
步骤6.3.1：令g＝1，i＝1；
步骤6.3.2：判断cdf(i)≤ug是否成立，成立则转至下一步，否则转至6.3.4；
步骤6.3.3：i＝i+1,返回步骤6.3.2；
步骤6.3.4：第i个粒子复制在第g个位置；
步骤6.3.4：g＝g+1；
步骤6.3.5：判断g≤N是否满足，满足转至步骤6.3.2，不满足则终止。
步骤7：利用Markov链蒙特卡罗方法(MCMC)保持重采样后粒子的多样性，对粒子群Φ″k进行更新,得到新的粒子群Φ″′k将权重w″i(k)重新进行权重归一化，得到权重w″′i(k)；具体步骤如下：
步骤7.1：获取已经得到的
步骤7.2：令i＝0；
步骤7.3:从[0,1]中抽样得到随机数e；
步骤7.4:根据建议分布和当前变量值采样得到候选值计算接收概率
其中为的目标分布，为的建议分布；
步骤7.5:判断是否成立，如果成立，则输出否则，输出
步骤7.6:i＝i+1，判断i<N是否成立，成立转至步骤7.3；否则，将步骤7.5中输出的粒子形成新的粒子群Φ″′k，转至7.7；
步骤7.7:对粒子的权重w″i(k)进行归一化，得到w″′i(k)。
步骤8：计算得到电力系统的状态估计值，所述状态估计值包括当前电力系统电压的幅值和相角；判断是否为步骤4跳转至本步骤，若是，则电力系统的状态估计值的公式为：若否，则电力系统的状态估计值的公式为
步骤9：判断k≥Ω是否成立，如果不成立，令k＝k+1转至步骤2，如果成立，结束电力系统的动态估计，其中Ω为预测截止时间。电力系统动态估计能够帮助调度中心分析和预测系统电力系统的运行趋势，对运行中发生的各种问题提出对策。
本实施例的系统描述：电力系统包括2个发电机和一个负荷，如图3所示，系统参数：功率基值：100MW，电压基值：230kV
其他参数见表1-表5：
表1线路阻抗参数(标么值)
表2线路量测(流出母线为正)
线路 首端量测 首端真值 末端量测 末端真值
连接线1-2 -61.3+j1.2 -63.32+j1.24 60+j2.4 63.99+j2.35
连接线1-3 -46.7-j14.8 -47.69-j14.79 45.9+j16.5 47.9+j16.54
连接线2-3 24-j6.6 24.93-j6.6 -24+j7.2 -24.86+j7.24
表3母线电压量测
母线 电压量测 真值(幅值/角度)
1 232 233/0
2 234.56 235.56/3.23
3 236.46 236.46/1.82
表4负荷量测(流出母线为正)
负荷 量测 真值
Load1 111+j13.50 111.02+j13.55
表4发电量测(流入母线为正)
发电机 量测 真值
发电机2 88-j4.24 88.92-j4.24
发电机3 23+j24 23.04+j23.78
本实施例通过以上系统进行了仿真验证，通过图4至图8，可以看出在稍微牺牲计算时间的情况下，无论是针对系统的某一个节点还是针对整个系统，在电压的幅值、相角的估计上，与扩展的卡尔曼滤波(EKF)、无极卡尔曼滤波(UKF)、粒子滤波(PF)、扩展的卡尔曼粒子滤波(EPF)相比，所提的无迹卡尔曼粒子滤波算法(UPF)在估计精度上有明显优势，说明了所提方法的有效性。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (3)

1.一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法，其特征在于：包括如下步骤：
步骤1：将当前电力系统状态变量进行初始化，获取k＝0时刻的粒子群Φ,其中k代表时间；粒子个数为N,For i＝1：N,从先验概率密度函数p(X0)中抽取初始化状态作为电力系统的初始状态，计算初始状态均值方差
其中X0＝[v0θ0]，X0为0时刻电力系统的初始状态，v0表示0时刻的节点电压的幅值，θ0表示0时刻的节点电压的相角，为0时刻第i个粒子的初始值，P0为初始状态估计误差的协方差，Q为过程噪声的协方差，R为测量噪声的协方差，E(·)代表求期望值的函数，(·)T代表(·)的转置；
步骤2：令k＝1；
步骤3：计算k-1时刻下粒子群中所有粒子的Sigma采样点集合，根据该集合预测出k时刻的Sigma采样点集合，并通过k时刻的Sigma采样点集合得出k时刻的状态粒子群Φk和粒子权重wi(k)；
步骤3.1：计算k-1时刻下粒子i的Sigma采样点集合，遍历k-1时刻的粒子群，得到所有粒子的Sigma采样点集合；
其中为k-1时刻第i个粒子的Sigma采样点，维数为n，为k-1时刻第i个粒子的方差，为k-1时刻第i个粒子的状态均值,λ为采样因子，b为自由参数，用来捕捉给定的高阶矩阵的信息，为矩阵的第j列，α为比例修正因子，取值为[0.0001,1]，设定初始值为α0＝0.01，其中α1到αt的计算步骤如下：
步骤3.1.1：令则上式可以转化为：其中ck-1表示k-1时刻的均值到sigma采样点的距离；
步骤3.1.2：计算
步骤3.1.3：令其中I为sigma采样点集合，为ck-1的最大值；
步骤3.1.4：计算k刻的比例修正因子其中tr(·)代表求矩阵的迹；
步骤3.2：根据步骤3.1中得到k-1时刻下粒子群中所有粒子的Sigma采样点集合预测出k时刻中所有粒子的Sigma采样点集合；
k时刻粒子i Sigma采样点集合的公式如下：
其中为第i个粒子k-1时刻对k时刻的预测值，f(·)为预测函数，为第i个粒子k-1时刻对k时刻的预测值的均值，为第i个粒子k-1时刻对k时刻的第a个预测值，为第i个粒子k-1时刻对k时刻预测值的方差，为第i个粒子k-1时刻对k时刻预测的量测值，为第i个粒子k-1时刻对k时刻预测的量测值的均值，h(·)为量测函数，为第i个粒子k-1时刻对k时刻预测的第a个量测值，代表sigma采样点求期望的权值、代表sigma采样点求方差的权值；
其中计算如下：
其中，a＝1,…,2n，为sigma采样点求方差的初始权值；为sigma采样点求期望的初始权值；α0为比例修正因子设定的初始值，β为高阶误差采样因子；
步骤3.3：将步骤3.2中得到的和融入到k时刻已知的实际量测值中并对k时刻粒子i的预测值的均值和k时刻粒子i的预测协方差进行更新；重复本步骤遍历k-1时刻下粒子群中所有粒子，得到k时刻所有粒子的预测值的均值和k时刻所有粒子的预测协方差
其中为已知量测zk时的方差，为已知量测zk时Xk的方差，Xk为电力系统k时刻的状态值，Kk为计算增益，Zk为k时刻的量测值；
步骤3.4：根据步骤3.3求出k时刻的采样更新粒子最终得到k时刻的粒子群Φk；
其中代表以为均值、以为方差的正态分布函数；
步骤3.5：通过PMU的量测和电力系统的静态估计形成混合量测，计算粒子的权重；其中PMU的量测和电力系统的静态估计形成混合量测具体步骤如下：
步骤3.5.1：PMU的量测是通过安装PMU的节点获取该节点的电压幅值和相角和电流；通过计算得到与其相连节点的电压的伪量测值，根据该值最粒子的权重进行更新；具体计算如下：
其中Y10为节点对地导纳，Y12为节点1、2直接互导纳，为流向节点1的电流，为节点1的电压，为节点2的电压；
步骤3.5.2：获取已有的电力系统的静态估计结果，选取静态估计精度高的结果用来对粒子的权重进行更新；
其中,d＝1,2,…M,代表静态估计值，V代表电压幅值，θ代表电压相角；代表安装PMU节点的量测值，通过向量节点计算得到的伪量测值,为k时刻第i个粒子第d个量测的初始权重，Zd(k)代表k时刻的第d个量测值，为k时刻第i个粒子第d个量测的初始权重的权重系数，σd为第d个量测的量测噪声，M代表量测的总数；
步骤3.5.3：最终权重为：
步骤3.6：将权重归一化，得到归一化权重wi(k)；归一化方法如下：
步骤4：对k时刻的粒子群进行粒子群分裂和权值调整；将粒子的权值按照高权值矩阵和低权值矩阵排列，将高权值矩阵中的后1/2的粒子分裂成两个权重减少1/2的粒子，把低权值的后1/2的粒子去除,形成新的粒子群Φ′k，重新进行权重归一化，得到权重w′i(k)；
步骤5:计算有效粒子Neff，定义阈值判断Neff<Nth是否成立，如果成立转到步骤6，如果不成立转到步骤8；其中Neff的计算如下：
步骤6：利用随机重采样方法，根据归一化权值wi′(k)大小，对粒子群Φ′k进行复制和淘汰得到粒子群For i＝1:N,重新设置权值权重
步骤7：利用Markov链蒙特卡罗方法保持重采样后粒子的多样性，对粒子群Φ″k进行更新,得到新的粒子群Φ″′k将权重w″i(k)重新进行权重归一化，得到权重w″′i(k)；
步骤8：计算得到电力系统的状态估计值，所述状态估计值包括当前电力系统电压的幅值和相角；判断是否为步骤4跳转至本步骤，若是，则电力系统的状态估计值的公式为：若否，则电力系统的状态估计值的公式为
步骤9：判断k≥Ω是否成立，如果不成立，令k＝k+1转至步骤2，如果成立，结束电力系统的动态估计，其中Ω为预测截止时间。
2.根据权利要求1所述的一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法，其特征在于：所述步骤6中随机重采样方法的具体步骤如下：
步骤6.1:产生[0,1]均匀分布的随机数组其中N为粒子数，g代表代表生成的随机数，并按升序排序；
步骤6.2:产生粒子权累积函数cdf,满足
步骤6.3:复制权值大的粒子，具体计算如下：
步骤6.3.1：令g＝1，i＝1；
步骤6.3.2：判断cdf(i)≤ug是否成立，成立则转至下一步，否则转至6.3.4；
步骤6.3.3：i＝i+1,返回步骤6.3.2；
步骤6.3.4：第i个粒子复制在第g个位置；
步骤6.3.4：g＝g+1；
步骤6.3.5：判断g≤N是否满足，满足转至步骤6.3.2，不满足则终止。
3.根据权利要求1所述的一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法，其特征在于：所述步骤7的具体步骤如下：
步骤7.1：获取已经得到的
步骤7.2：令i＝0；
步骤7.3:从[0,1]中抽样得到随机数e；
步骤7.4:根据建议分布和当前变量值采样得到候选值计算接收概率
其中为的目标分布，为的建议分布；
步骤7.5:判断是否成立，如果成立，则输出否则，输出
步骤7.6:i＝i+1，判断i<N是否成立，成立转至步骤7.3；否则，将步骤7.5中输出的粒子形成新的粒子群Φ″′k，转至7.7；
步骤7.7:对粒子的权重w″i(k)进行归一化，得到w″′i(k)。
CN201910409884.9A 2019-05-16 2019-05-16 一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法 Active CN110162871B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910409884.9A CN110162871B (zh) 2019-05-16 2019-05-16 一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910409884.9A CN110162871B (zh) 2019-05-16 2019-05-16 一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法

Publications (2)

Publication Number Publication Date
CN110162871A true CN110162871A (zh) 2019-08-23
CN110162871B CN110162871B (zh) 2023-01-31

Family

ID=67631090

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910409884.9A Active CN110162871B (zh) 2019-05-16 2019-05-16 一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法

Country Status (1)

Country Link
CN (1) CN110162871B (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110501686A (zh) * 2019-09-19 2019-11-26 哈尔滨工程大学 基于一种新型自适应高阶无迹卡尔曼滤波的状态估计方法
CN110806733A (zh) * 2019-10-30 2020-02-18 中国神华能源股份有限公司国华电力分公司 火电厂设备监测方法、装置及电子设备
CN111061151A (zh) * 2019-11-21 2020-04-24 东北大学 一种基于多元卷积神经网络的分布式能源状态监测方法
CN111353398A (zh) * 2020-02-24 2020-06-30 华北电力大学（保定） 一种电力信号同步相量参数估计方法
CN111797974A (zh) * 2020-06-01 2020-10-20 武汉大学 一种联合粒子滤波和卷积神经网络的电力系统状态估计方法
CN112613222A (zh) * 2021-01-04 2021-04-06 重庆邮电大学 基于改进粒子滤波的倾斜探测电离层muf短期预报方法
CN112688596A (zh) * 2020-12-04 2021-04-20 西安理工大学 基于upf永磁直线同步电动机的无位置传感器的控制方法
CN115070765A (zh) * 2022-06-27 2022-09-20 江南大学 一种基于变分推断的机器人状态估计方法及系统
CN118171050A (zh) * 2024-05-14 2024-06-11 河海大学 用于组合导航系统的低计算成本非线性卡尔曼滤波方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017035964A1 (zh) * 2015-08-31 2017-03-09 中车大连电力牵引研发中心有限公司 一种电力系统的负荷特性确定方法及系统
CN106844952A (zh) * 2017-01-20 2017-06-13 河海大学 基于无迹粒子滤波理论的发电机动态状态估计方法
CN107765179A (zh) * 2017-06-26 2018-03-06 河海大学 一种适用于量测丢失的发电机动态状态估计方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017035964A1 (zh) * 2015-08-31 2017-03-09 中车大连电力牵引研发中心有限公司 一种电力系统的负荷特性确定方法及系统
CN106844952A (zh) * 2017-01-20 2017-06-13 河海大学 基于无迹粒子滤波理论的发电机动态状态估计方法
CN107765179A (zh) * 2017-06-26 2018-03-06 河海大学 一种适用于量测丢失的发电机动态状态估计方法

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110501686A (zh) * 2019-09-19 2019-11-26 哈尔滨工程大学 基于一种新型自适应高阶无迹卡尔曼滤波的状态估计方法
CN110806733A (zh) * 2019-10-30 2020-02-18 中国神华能源股份有限公司国华电力分公司 火电厂设备监测方法、装置及电子设备
CN111061151A (zh) * 2019-11-21 2020-04-24 东北大学 一种基于多元卷积神经网络的分布式能源状态监测方法
CN111353398A (zh) * 2020-02-24 2020-06-30 华北电力大学（保定） 一种电力信号同步相量参数估计方法
CN111797974A (zh) * 2020-06-01 2020-10-20 武汉大学 一种联合粒子滤波和卷积神经网络的电力系统状态估计方法
CN112688596A (zh) * 2020-12-04 2021-04-20 西安理工大学 基于upf永磁直线同步电动机的无位置传感器的控制方法
CN112613222A (zh) * 2021-01-04 2021-04-06 重庆邮电大学 基于改进粒子滤波的倾斜探测电离层muf短期预报方法
CN112613222B (zh) * 2021-01-04 2023-09-15 重庆邮电大学 基于改进粒子滤波的倾斜探测电离层muf短期预报方法
CN115070765A (zh) * 2022-06-27 2022-09-20 江南大学 一种基于变分推断的机器人状态估计方法及系统
CN115070765B (zh) * 2022-06-27 2023-06-13 江南大学 一种基于变分推断的机器人状态估计方法及系统
CN118171050A (zh) * 2024-05-14 2024-06-11 河海大学 用于组合导航系统的低计算成本非线性卡尔曼滤波方法

Also Published As

Publication number Publication date
CN110162871B (zh) 2023-01-31

Similar Documents

Publication Publication Date Title
CN110162871A (zh) 一种基于无迹卡尔曼粒子滤波的电力系统动态估计方法
Yang et al. A reliability assessment approach for electric power systems considering wind power uncertainty
CN111353652B (zh) 一种风电出力短期区间预测方法
CN103326358A (zh) 基于同步相角测量装置的电力系统动态状态估计方法
CN103472398A (zh) 基于扩展卡尔曼粒子滤波算法的动力电池soc估计方法
CN107145707B (zh) 一种计及光伏出力不确定性和全寿命周期成本的配电网变压器规划方法
CN110601174B (zh) 一种基于深度学习的负荷建模与在线修正方法
CN107831448A (zh) 一种并联型电池系统的荷电状态估计方法
CN111064180A (zh) 基于ami潮流匹配的中压配电网拓扑检测与辨识方法
CN109754013A (zh) 一种基于无迹卡尔曼滤波的电力系统混合量测融合方法
CN109921426A (zh) 基于cv-kde的风电并网系统概率潮流计算方法
CN110222309A (zh) 一种基于鲁棒容积卡尔曼滤波的发电机动态估计方法
CN115081299A (zh) 一种基于upf的电力系统鲁棒辅助预测状态估计方法
CN115603320A (zh) 一种基于广域量测数据的实时感知与在线评估方法及系统
CN104537233B (zh) 一种基于核密度估计的配电网伪量测生成方法
CN110210690B (zh) 一种配电系统微型同步相量测量单元优化配置方法
CN116542021A (zh) 一种水文-水动力学耦合的河道型水库调洪演算方法
CN116865343B (zh) 分布式光伏配电网的无模型自适应控制方法、装置及介质
CN107958120B (zh) 一种基于幂级数展开的系统戴维南等值参数计算方法
CN105071388A (zh) 一种基于极大似然估计的配电网状态估计方法
CN110474323B (zh) 一种电力系统惯性时间常数测量方法
CN106599541B (zh) 一种动态电力负荷模型的结构和参数在线辨识方法
CN115293090A (zh) 基于重构数据处理的多谐波源责任量化方法
CN115409335A (zh) 基于深度学习考虑未知扰动类型的电力系统扰动识别方法
CN109256773A (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