DC State Estimation

[Saeed Ahmed]

Hi Community,
I need you to help and guide me to improve my results .... I have writtent
the code but the estimated measurement matrix is not correct.

First i run an AC power flow case '4gs' and then make the DC State
estimation code to compare the state(angles) and estimated measurements
%%              Select a case and run AC Power Flow.
mpc = loadcase('case4gs');% Load case data
pfres = runpf(mpc);% Runs an AC Power Flow based on Newton-Raphsons Method
casesystem = size(mpc.bus, 1);
%-----------------------------------------------------------------------------------------------------------%
                                           % Start DC State Estimation
%-----------------------------------------------------------------------------------------------------------%
%% Step-1: Make the Measurement Matrix( Real Power Injections and Real
Power Transmitted)
%-----------------------------------------------------------------------------------------------------------%
% (I). Real Power Injection Matrix
p_bus = real(makeSbus(mpc.baseMVA, mpc.bus, mpc.gen)) * mpc.baseMVA; % Real
Power Injections at the bus nodes
% p_bus(1,:)=[]; % Eleminate the Referece Bus
%-----------------------------------------------------------------------------------------------------------%
% (II). Real Power Flows: Power Transmitted from bus i to bus j)
n_lines = size(pfres.branch,1);

% Active Power Flow of Sending Nodes(From Bus)
PF_S = [ones(n_lines,1), 4*ones(n_lines,1), pfres.branch(:,14)/mpc.baseMVA,
pfres.branch(:,1:2), 0.005*ones(n_lines,1)];
pf_s=PF_S(:,3)*mpc.baseMVA;

% Real Power Flows(To Bus) of Recieving Nodes(To Bus)
PF_R = [ones(n_lines,1),4*ones(n_lines,1), pfres.branch(:,16)/mpc.baseMVA,
pfres.branch(:,1:2), 0.005*ones(n_lines,1)];
pf_r=PF_R(:,3)*mpc.baseMVA;

%-----------------------------------------------------------------------------------------------------------%
% (III) The Measurement Data Matrix obtained by combining Active Power
Injections and Active Power Flows)
Z_M = [p_bus; pf_s; pf_r]; % Only Real Power Injections and Power Flows due
to DC State Estimation

% (IV) Introducing random meaurement errors according to variances
% for i = 1:size(z, 1)
%     Z(i,3) = normrnd( zdata(i,3), zdata(i,6))
% end

%% Step-2:   Make H Matrix (H Matrix has three parts)
%( H=Part (1) B Matrix;Part (2) Connection Mat(Real Tx Pwr From bus i>j
);Part (3) Connection Mat(Real Tx Pwr From bus i>j)
%-----------------------------------------------------------------------------------------------------------%
% Stage-1:   Make Ybus(Bus Admittence Matrix)
[Ybus, YF, YT] = makeYbus(mpc);
Ybus=full(imag(Ybus))
% Stage-2 Make the B Matrix and phase shift injection vectors needed for a
DC power flow.
[Bbus, Bf, Pbusinj, Pfinj] = makeBdc(mpc.baseMVA, mpc.bus, mpc.branch);
H_P1=full(Bbus)%------ (1st Part of H-Matrix)
% Stage-3 Make the Connection Matrix or Incidence Matrix
[Ainc] = makeIncidence(mpc.bus, mpc.branch);
Ainc = full(Ainc)
% Stage-3 Make the 2nd Part & 3rd of H Matrix
[row col]= size (Ainc);
for i=1:row
    for j=1:col
        H_P2(i,j)=Ainc(i,j)*Ybus(i,j);%------ (2nd Part of H-Matrix)
        H_P3(i,j)=(-1)*Ainc(i,j)*Ybus(i,j);%------ (3rd Part of H-Matrix)
    end
end
% Stage-4 Make the H-Matrix (Combining 1st,2nd and Thrird Part)

H=[H_P1;H_P2;H_P3];

%% Step-3 Make the Co-Var Matrix
R=  0.1^2*diag([1 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
Ang_est=(H'*W *H)^-1* H' * W * Z_M
%% Step-4 Make the Estimated Measurments Matrix
Z_est=H*ang_est; % Estimated Measurements

Regards
Saeed Ahmed