1 主要内容

该程序对应文章《Power System Dynamic State Estimation Using Extended and Unscented Kalman Filters》，电力系统状态的准确估计对于提高电力系统的可靠性、弹性、安全性和稳定性具有重要意义，虽然近年来测量设备和传输技术的发展大大降低了测量和传输误差，但这些测量仍然不能完全没有测量噪声。因此，需要对噪声测量值进行过滤，从而获得准确的电力系统运行动态。本程序采用两种方法，分别是扩展卡尔曼滤波(EKF)和无迹卡尔曼滤波(UKF)，对电力系统进行动态状态估计，以39节点系统为算例验证了方法的有效性。

2 部分代码

clear;
clc;

%% Power Flow calculation
% Y=Ybus_new(case9_new_Sauer); % 9 bus system data obtained from MATPOWER
% result=runpf(case9_new_Sauer); % run ac power flow, in this case default NR is used


%result= runpf(case5_Overbye); 
%Y=Ybus_new(case5_Overbye); 

% Y=Ybus_new(case14); 
% result=runpf(case14);
% 
Y=Ybus_new(case39);
result=runpf(case39);


Vmag=result.bus(:, 8); % Pu voltage magnitude of each buses 
Vph=result.bus(:, 9); % angle in degree
V=Vmag.*exp(1j*Vph*pi/180); 
P_jQ=conj(V).*(Y*V); % Net Power at each node
S=conj(P_jQ);
S=S/100; 
Sg=result.gen(:, 2)+1j*result.gen(:, 3); 
Sg=Sg/100;


%% machine data for 9 bus system
% Xd=[0.06080; 0.11980; 0.18130];
% R=[0;0;0];
% H=[23.64; 6.4; 3.010];
% M=H/(pi*60); 
%D=[0.0125;0.0034;0.0016];
% 
%% Data of 9 bus system from Peter Sauer.
% Xd=[0.06080; 0.11980; 0.18130];
% R=[0;0;0];
% H=[23.64; 6.4; 3.01];
% %H=[13.64; 6.4; 3.01]; 
% D=[0.0255; 0.00663; 0.00265]; 
% %D=[9.6; 2.5; 1]; % If we use this value need to devide the D term by 2*pi*60 
% f0=60; 
% w_syn=2*pi*f0; 
% M=2*H/w_syn; 
% gen_bus=result.gen(:, 1); 

%% machine data for 14 bus system 
% % Machine data 
% H=[5.1498; 6.54; 6.54; 5.06; 5.06];
% Xd=[0.2995; 0.185; 0.185; 0.232; 0.232];
% R=zeros(length(Xd), 1); 
%  
% f0=60; 
% w_syn=2*pi*f0; 
% D=[2; 2; 2; 2; 2]/w_syn;
% 
% M=2*H/w_syn; 
% 
% gen_bus=result.gen(:, 1); 

%% Overbye data for 5 bus system 
% Xd=[0.05; 0.025]; 
% R=[0; 0]; 
% H=[]; 
% D=[]; 
% f0=60; 
% w_syn=2*pi*f0; 
% M=2*H/w_syn; 
% gen_bus=result.gen(:, 1); 

%% Case 39 bus data 
Xd=[0.006; 0.0697; 0.0531; 0.0436; 0.132; 0.05; 0.049; 0.057; 0.057; 0.031]; 
H=[500; 30.3; 35.8;28.6; 26; 34.8; 26.4; 24.3; 34.5; 42]; 
R=zeros(length(Xd), 1); 
f0=60; 
w_syn=2*pi*f0; 

D=[0; 0;0 ;0; 0; 0; 0; 0; 0; 0]; 
D=D/w_syn;

M=2*H/w_syn; 

gen_bus=result.gen(:, 1); 

%% case 145



%% calculate Y22
Y22=diag(1./(1j*Xd)); 

%% Calculation of Y11
SL=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;