Hi Community,
I have developed the *DC State Estimation Code* for Approximating the power
flow.
I took the simplest case 'case 4gs' and run ac pf. Then i tried to
approximate the state estimation using DC (WLS function) . The Estimated
Measurements are correct but the phase angles are not.
I would request you to guide me to rectify my fault.
1) Is the construction of H Matrix is right, please
2) What is my mistake
m file is attached herewith, please.
Regards
Saeed Ahmed
%% Select a case and run AC Power Flow.
mpc = loadcase('case4gs');% Load case data
pfresult = rundcpf(mpc);% Runs an AC Power Flow based on Newton-Raphsons Method
%-----------------------------------------------------------------------------------------------------------%
% Start DC State Estimation
%-----------------------------------------------------------------------------------------------------------%
%% Step-1: Make the Measurement Matrix( Real Bus Power Injections & Real
Line/Branch Power Flows)
%-----------------------------------------------------------------------------------------------------------%
n_lines = size(pfresult.branch,1); % Number of Branches
% (I). Real Power Injection Matrix
Z_Pbus = real(makeSbus(mpc.baseMVA, mpc.bus, mpc.gen)) * mpc.baseMVA; % Real
Power Injections at the bus nodes
%-----------------------------------------------------------------------------------------------------------%
% (II). Active/Real Power Flow of Sending Nodes(From Bus)
PF_S = [ones(n_lines,1), 4*ones(n_lines,1), pfresult.branch(:,14)/mpc.baseMVA,
pfresult.branch(:,1:2), 0.005*ones(n_lines,1)];
% (III). Active/Real Power Flows Recieving Nodes(To Bus)
PF_R = [ones(n_lines,1),4*ones(n_lines,1), pfresult.branch(:,16)/mpc.baseMVA,
pfresult.branch(:,1:2), 0.005*ones(n_lines,1)];
Z_PF= [PF_S;PF_R] %(Combining I & II: Line/Branch Power Flows)
%-----------------------------------------------------------------------------------------------------------%
%(III). Introducing random meaurement errors in power flow according to
variances
for i = 1:size(Z_PF, 1)
Z_PF(i,3) = normrnd( Z_PF(i,3), Z_PF(i,6))
end
Z_pf=Z_PF(:,3)*mpc.baseMVA
%-----------------------------------------------------------------------------------------------------------%
% (III) The Measurement Data Matrix obtained by combining Active Power
Injections and Active Power Flows)
Z_Meter = [Z_Pbus; Z_pf]; % Only Real Power Injections and Power Flows due to
DC State Estimation
% Z_M(1)=136.81;
Z_M(1,:)=[]
%-----------------------------------------------------------------------------------------------------------%
%% Step-2: Make H Matrix (H Matrix has three parts)
%-----------------------------------------------------------------------------------------------------------%
% Method: Using B-Matrix (B(DC) )
[Bbus, Bf, Pbusinj, Pfinj] = makeBdc(mpc.baseMVA, mpc.bus, mpc.branch);
H_P1 = full(Bbus)%------ (1st Part of H-Matrix)[Ybus, YF, YT] = makeYbus(mpc);
H_P2 = full(Bf)
H_P3 = -full(Bf)
H = [H_P1; H_P2; H_P3]
H(:,1)=[]
H(1,:)=[]
%-----------------------------------------------------------------------------------------------------------%
%% Step-3 Make the Co-Var Matrix
%-----------------------------------------------------------------------------------------------------------%
R = 0.01^2*diag([1 1 1 1 1 1 1 1 1 1 1]); % Co-Variance Matrix (0.01)
W = R^-1; % Reciprocal of Diagonal Matrix Elements
%-----------------------------------------------------------------------------------------------------------%
%% Step-4 Estimate the State(Voltage Angles) WLS Criteria
%-----------------------------------------------------------------------------------------------------------%
% Estimated Angle
Ang_estd=((H'*W *H)^-1* H' * W * Z_M)
Ang_true=[-1.059;-2.032;1.783] % From IEEE (case 4gs; Just to compare with
estimated angles)
%-----------------------------------------------------------------------------------------------------------%
% Estimated Measurements
Z_est=H*Ang_estd
%%-----------------------------------------------------------------------------------------------------------%
% Testing/ Checking
K = H'*R^-1*H;
L =(H'*R^-1*H)^-1;
K*L;
%% Residual
r = Z_M - Z_est
