前言
电力系统状态估计是电力系统调度中心的能量管理系统（EMS）的核心功能之一，其功能是根据电力系统的各种量测信息，估计出电力系统当前的运行状态。现代电网的安全经济运行依赖于能量管理系统（EMS），而能量管理系统的众多功能又可分成针对电网实时变化进行分析的在线应用和针对典型潮流断面进行分析的离线应用两大部分。
电力系统状态估计可以说是大部分在线应用的高级软件的基础，如果电力系统状态估计结果不准确，后续的任何分析计算将不可能得到准确的结果。
分享一些基础电力系统同状态估计算法的优缺点
本程序主要利用扩展卡尔曼滤波法和无迹卡尔曼滤波法对电力系统进行动态状态估计，下面分别介绍这两种算法。
卡尔曼滤波简介
在控制系统和信号处理领域，卡尔曼滤波是一种强大的递归算法，用于估计动态系统的状态。由于其高效的状态估计能力，卡尔曼滤波在许多领域得到广泛应用，包括导航、目标跟踪、机器人技术等。
1.1 基本思想
卡尔曼滤波的基本思想是通过融合系统的动态模型和实际测量数据，得到对系统状态的最优估计。它在面对不完全和带有噪声的传感器数据时表现出色，通过迭代更新，逐步优化状态的估计值。
卡尔曼滤波通过两个主要步骤实现状态估计：
1）预测步骤：使用系统的动态模型进行状态预测。
2）更新步骤：使用测量值对预测的状态进行修正。
1.2 扩展卡尔曼滤波的引入
尽管卡尔曼滤波在处理线性系统方面表现出色，但对于非线性系统，特别是那些包含非线性状态转移或观测方程的系统，卡尔曼滤波不再适用。为了解决这个问题，扩展卡尔曼滤波应运而生。
1.3 扩展卡尔曼滤波与卡尔曼滤波的区别
卡尔曼滤波是基于线性动态模型和线性观测模型的，因此对于非线性系统不太适用。扩展卡尔曼滤波通过在线性化非线性系统的模型，将卡尔曼滤波的思想扩展到非线性问题中。这主要涉及到使用雅可比矩阵来近似非线性函数。
1.4 扩展卡尔曼滤波的数学原理
扩展卡尔曼滤波在数学上与卡尔曼滤波相似，但状态转移和观测方程的雅可比矩阵被引入，以适应非线性系统。预测和更新的公式仍然保持卡尔曼滤波的形式，但雅可比矩阵的引入使其能够处理更广泛的问题。
1.5 扩展卡尔曼滤波的算法思想
扩展卡尔曼滤波的算法思想可分为预测步骤和更新步骤：
1）预测步骤
状态预测：使用非线性状态转移函数预测系统的下一个状态。
协方差预测：利用状态方程的雅可比矩阵进行协方差预测。
2）更新步骤
卡尔曼增益计算：利用观测方程的雅可比矩阵计算卡尔曼增益。
状态更新：使用测量值更新状态估计。
协方差更新：利用卡尔曼增益更新协方差矩阵。
无迹卡尔曼滤波法简介
无迹卡尔曼滤波（Unscented Kalman Filter，UKF）是一种用于非线性系统状态估计的滤波器，它是对经典的卡尔曼滤波器的扩展，适用于非线性系统。UKF通过一种称为“无迹变换”的技术，将系统状态的不确定性进行传递和更新。
2.1 无迹卡尔曼滤波的关键特点和步骤：
1）非线性系统状态估计：传统的卡尔曼滤波器适用于线性系统，而UKF则可以处理非线性系统，如通过非线性方程描述的动力学系统。
2）无迹变换：UKF使用无迹变换来进行状态传递和更新。无迹变换通过一组称为“sigma点”的特殊采样点来近似概率分布的参数，并在这些点上执行状态传递和观测更新。这些sigma点的选择是通过在高斯分布上以状态和协方差为中心进行采样而确定的。
3）状态预测：在每个时间步骤中，UKF首先利用当前的状态估计和协方差估计生成一组sigma点。然后，这些sigma点通过系统的动态方程进行状态传递，从而预测下一个时间步骤的状态。
4）测量更新：一旦观测到新的测量数据，UKF使用预测的状态和测量模型来生成一组新的sigma点。然后，通过测量方程将这些sigma点映射到观测空间中，以进行状态更新。
5）估计更新：利用测量更新后的信息，UKF通过融合预测和观测信息来更新状态估计和协方差估计。
6）迭代：这些步骤随着时间的推移进行迭代，以提供对系统状态的连续估计。
程序介绍
电力系统的动态状态估计对于广域控制是必要的。在同步电机的状态中，关于转子角度和速度偏差的精确、准确和及时的信息有助于提高电力系统的可靠性和稳定性。程序采用两种常用的非线性估计方法：扩展卡尔曼滤波器（EKF）和无迹卡尔曼滤波器（UKF），利用相量测量单元（PMU）提供的高速率数据来估计简单电力系统的主要状态。程序算例丰富、注释清晰、干货满满，可扩展性和创新性很高！下面对文章和程序做简要介绍！
适用平台：Matlab+Matpower
程序结果
（1）EKF状态估计结果
（2）UKF状态估计结果
部分程序
%% 计算 Y11SL=result.bus(:, 3)+1j*result.bus(:, 4); SL=SL/100;YL=conj(SL)./(abs(V).^2); % Y11=Y+diag(YL);Y11(gen_bus, gen_bus)=Y11(gen_bus, gen_bus)+Y22;%% 计算 Y12 和 Y21Y12=zeros(length(result.bus(:,1)), length(result.gen(:,1)));%Y12(gen_bus, gen_bus)=Y12(gen_bus, gen_bus)-Y22;Y12(q,k)=-1/(R(k)+Xd(k)*1j);Y21=transpose(Y12);%% 计算故障前简化矩阵Ybf=Y22-Y21*inv(Y11)*Y12 ;% 重新构建母线矩阵RV(:, :, 1)=-inv(Y11)*Y12;%% 在此处输入故障以计算故障后和故障减少时间% 矩阵f11=4;F=[4 14];f1=F(1);f2=F(2);%% 故障中Y11df=Y11; Y11df(f11, :)=[];Y11df(:,f11)=[];Y12df=Y12; Y12df(f11, :)=[];Y21df=transpose(Y12df);% 在故障减少矩阵期间Ydf=Y22-Y21df*inv(Y11df)*Y12df;RV(:, :, 2)=zeros(size(RV(:, :, 1)));RV(1:end-1, :, 2)=RV(1:end-1, :, 2)-inv(Y11df)*Y12df;%% 故障后 Y11Y11after=Y11;Y11after(f1,f2)=0;Y11after(f2,f1)=0;if (f1==result.branch(i,1)&& f2==result.branch(i,2))||(f2==result.branch(i,1)&& f1==result.branch(i,2))Y11after(f1,f1)=Y11after(f1,f1)-result.branch(i,5)*1j/2-1/(result.branch(i,3)+result.branch(i,4)*1j);Y11after(f2,f2)=Y11after(f2,f2)-result.branch(i,5)*1j/2-1/(result.branch(i,3)+result.branch(i,4)*1j);% 故障后矩阵Yaf=Y22-Y21*inv(Y11after)*Y12 ;%RV_af=-inv(Y11after)*Y12 ; RV(:, :, 3)=-inv(Y11after)*Y12;%% 初始化deltt=0.0005; t_SW=1; t_FC=1.0333; t_max=10; Ig=conj(Sg./V(1:length(result.gen(:, 1))));E0=V(gen_bus)+Ig.*(R+1j*Xd); % 发电机端电压E_abs=abs(E0);I0=Ybf*E0;delta0=angle(E0)*180/pi;w0=zeros(length(Xd), 1);X_0=[angle(E0); w0];% 初始功率PG0=real(E0.*conj(I0)); PM=PG0;QG0=imag(E0.*conj(I0));YBUS(:, :, 1)=Ybf; YBUS(:, :, 2)=Ydf; YBUS(:, :, 3)=Yaf;n=length(gen_bus); s=length(result.bus(:, 1)); %% Estimated State:% 状态和测量的数量ns=2*n; nm=2*n+2*s; % 协方差矩阵P=sig^2*eye(ns); % 误差协方差矩阵Q=sig^2*eye(ns); % 系统噪声协方差矩阵R=sig^2*eye(nm); % 测量噪声协方差矩阵 X_hat=X_0;X_est=[]; X_mes=[]; % 初始状态l%扩展卡尔曼滤波（EKF）法for k=0:deltt:t_max% Ybus 根据需要重构
部分内容源自网络，侵权联系删除！
欢迎感兴趣的小伙伴点击文末阅读原文获取完整版代码，小编会不定期更新高质量的学习资料、文章和程序代码，为您的科研加油助力！
更多创新性预测及优化程序请点击公众号首页：《预测优化》|《综合能源》链接！
