目录
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;
3 程序结果
原文结果：
