Commit 388998ae authored by Angioni's avatar Angioni

added code to calculate jacobian and covariance in no-per units

parent 006d1ce9
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the Jacobian matrix for the Current Real Imaginary Distribution
% System State Estimator (IRIDSSE) in normal units (no per units)
%
% @author Andrea Angioni <aangioni@eonerc.rwth-aachen.de>
% @copyright 2018, Institute for Automation of Complex Power Systems, EONERC
% @license GNU General Public License (version 3)
%
% dsse
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [H] = Jacobian_m_IRIDSSE_nopu(GridData)
if strcmp(GridData.type_of_model,'single_phase')==1
R1 = GridData.R1*GridData.base_impedance;
X1 = GridData.X1*GridData.base_impedance;
if GridData.inj_status == 0
H = zeros(GridData.MeasNum, 2 + 2*GridData.Lines_num);
dim_state_inv = 0;
else
H = zeros(GridData.MeasNum, 11*GridData.DM.NGF+10*GridData.DM.NGS+2*GridData.DM.Nload + 2 + 2*GridData.Lines_num);
dim_state_inv = 11*GridData.DM.NGF+10*GridData.DM.NGS+2*GridData.DM.Nload;
end
for n = 1 : GridData.MeasNum
if GridData.TypeMeas(n,1)==1 %derivative of injection real and imaginary current
for x = 1:GridData.Lines_num
if GridData.topology(3,x) == GridData.LocationMeas(n,1)
H(n,dim_state_inv + 2+2*x-1) = 1;
end%derivative with respect to entering current
if GridData.topology(2,x) == GridData.LocationMeas(n,1)
H(n,dim_state_inv + 2+2*x-1) = - 1;
end%derivative with respect to exiting current
end
end
if GridData.TypeMeas(n,1)==2 %reactive power measurements, both SM and pseudo
for x = 1:GridData.Lines_num
if GridData.topology(3,x) == GridData.LocationMeas(n,1)
H(n,dim_state_inv + 2+2*x) = 1;
end%derivative with respect to entering current
if GridData.topology(2,x) == GridData.LocationMeas(n,1)
H(n,dim_state_inv + 2+2*x) = - 1;
end%derivative with respect to exiting current
end
end
if GridData.TypeMeas(n,1)==3 %voltage magnitude measurements - converted to real voltage
for x = 1 : GridData.Lines_num
if GridData.A(x,GridData.LocationMeas(n,1)) == 1
H(n,dim_state_inv + 2+ 2*x-1) = - R1(1,x);
H(n,dim_state_inv + 2+ 2*x) = + X1(1,x);
%Vj = Vi - Zline*Iline: derivarive with respect to Iline ->
%it concern only the second bus of the line
end
end
%if GridData.LocationMeas(n,1) == 1
H(n,dim_state_inv + 1) = 1;
%end %derivative with respect to slack bus
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if GridData.TypeMeas(n,1)==4 %voltage phase angle measurements - converted to imaginary voltage
for x = 1 : GridData.Lines_num
if GridData.A(x,GridData.LocationMeas(n,1)) == 1
H(n,dim_state_inv + 2+ 2*x-1) = - X1(1,x);
H(n,dim_state_inv + 2+ 2*x) = - R1(1,x);
end
end
%if GridData.LocationMeas(n,1) == 1
H(n,dim_state_inv + 2) = 1;
%end %derivative with respect to slack bus
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if GridData.TypeMeas(n,1)== 5 || GridData.TypeMeas(n,1)== 7%current magnitude measurements - converted to real current and active power flow measurements
H(n,dim_state_inv + 2+ 2*GridData.LocationMeas(n,1) -1 ) = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if GridData.TypeMeas(n,1)== 6 || GridData.TypeMeas(n,1)== 8%current phase measurements - converted to imaginary current and reactive power flow measurements
H(n,dim_state_inv + 2+ 2*GridData.LocationMeas(n,1)) = 1;
end
if GridData.inj_status == 1
if GridData.TypeMeas(n,1)==100
H(n, GridData.LocationMeas(n,1)) = 1; %dynamic state
end
end
end
%the imaginary part of the slack bus is not a state (phase angle = 0),
%therefore it can be deleted
H(:,dim_state_inv + 2)=[];
elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridData.type_of_model,'three_phase_unbalance')==1
R1 = GridData.R1*GridData.base_impedance;
X1 = GridData.X1*GridData.base_impedance;
H = zeros(GridData.MeasNum, 6 + 6*GridData.Lines_num);
for n = 1 : GridData.MeasNum
if GridData.TypeMeas(n,1)==1 %derivative of injection real and imaginary current
for x = 1:GridData.Lines_num
if GridData.topology(3,x) == GridData.LocationMeas(n,1)
H(n, 6+6*(x-1)+ 2*(GridData.PhaseMeas(n,1)-1) + 1 ) = 1;
end%derivative with respect to entering current
if GridData.topology(2,x) == GridData.LocationMeas(n,1)
H(n, 6+6*(x-1)+ 2*(GridData.PhaseMeas(n,1)-1) + 1 ) = - 1;
end%derivative with respect to exiting current
end
end
if GridData.TypeMeas(n,1)==2 %reactive power measurements, both SM and pseudo
for x = 1:GridData.Lines_num
if GridData.topology(3,x) == GridData.LocationMeas(n,1)
H(n, 6+6*(x-1)+ 2*(GridData.PhaseMeas(n,1)) ) = 1;
end%derivative with respect to entering current
if GridData.topology(2,x) == GridData.LocationMeas(n,1)
H(n, 6+6*(x-1)+ 2*(GridData.PhaseMeas(n,1)) ) = - 1;
end%derivative with respect to exiting current
end
end
if GridData.TypeMeas(n,1)==3 %voltage magnitude measurements - converted to real voltage
for x = 1 : GridData.Lines_num
if GridData.A(x,GridData.LocationMeas(n,1)) == 1 %if the line x supplies the node m
if GridData.PhaseMeas(n,1)==1; a=2; b=3; end
if GridData.PhaseMeas(n,1)==2; a=1; b=3; end
if GridData.PhaseMeas(n,1)==3; a=1; b=2; end
H(n, 6+6*(x-1)+ 2*(GridData.PhaseMeas(n,1)-1) + 1) = - R1(GridData.PhaseMeas(n,1),3*(x-1) + GridData.PhaseMeas(n,1));
H(n, 6+6*(x-1)+ 2*(a-1) + 1) = - R1(GridData.PhaseMeas(n,1),3*(x-1) + a);
H(n, 6+6*(x-1)+ 2*(b-1) + 1) = - R1(GridData.PhaseMeas(n,1),3*(x-1) + b);
H(n, 6+6*(x-1)+ 2*(GridData.PhaseMeas(n,1)) ) = + X1(GridData.PhaseMeas(n,1),3*(x-1) + GridData.PhaseMeas(n,1));
H(n, 6+6*(x-1)+ 2*(a) ) = + X1(GridData.PhaseMeas(n,1),3*(x-1) + a);
H(n, 6+6*(x-1)+ 2*(b) ) = + X1(GridData.PhaseMeas(n,1),3*(x-1) + b);
%it concerns only the second bus of the line
end
end
H(n, 2*GridData.PhaseMeas(n,1)-1) = 1;
%end %derivative with respect to slack bus
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if GridData.TypeMeas(n,1)==4 %voltage phase angle measurements - converted to imaginary voltage
for x = 1 : GridData.Lines_num
if GridData.A(x,GridData.LocationMeas(n,1)) == 1
if GridData.PhaseMeas(n,1)==1; a=2; b=3; end
if GridData.PhaseMeas(n,1)==2; a=1; b=3; end
if GridData.PhaseMeas(n,1)==3; a=1; b=2; end
H(n, 6+6*(x-1)+ 2*(GridData.PhaseMeas(n,1)-1) + 1) = - X1(GridData.PhaseMeas(n,1),3*(x-1) + GridData.PhaseMeas(n,1));
H(n, 6+6*(x-1)+ 2*(a-1) + 1) = - X1(GridData.PhaseMeas(n,1),3*(x-1) + a);
H(n, 6+6*(x-1)+ 2*(b-1) + 1) = - X1(GridData.PhaseMeas(n,1),3*(x-1) + b);
H(n, 6+6*(x-1)+ 2*(GridData.PhaseMeas(n,1)) ) = - R1(GridData.PhaseMeas(n,1),3*(x-1) + GridData.PhaseMeas(n,1));
H(n, 6+6*(x-1)+ 2*(a) ) = - R1(GridData.PhaseMeas(n,1),3*(x-1) + a);
H(n, 6+6*(x-1)+ 2*(b) ) = - R1(GridData.PhaseMeas(n,1),3*(x-1) + b);
end
end
H(n, 2*GridData.PhaseMeas(n,1)) = 1;
%end %derivative with respect to slack bus
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if GridData.TypeMeas(n,1)== 5 || GridData.TypeMeas(n,1)== 7%current magnitude measurements - converted to real current and active power flow measurements
H(n, 3*(2+ 2*(GridData.LocationMeas(n,1) -1) ) + 2*(GridData.PhaseMeas(n,1)-1) + 1 ) = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if GridData.TypeMeas(n,1)== 6 || GridData.TypeMeas(n,1)== 8%current phase measurements - converted to imaginary current and reactive power flow measurements
H(n, 3*(2+ 2*(GridData.LocationMeas(n,1) -1) ) + 2*GridData.PhaseMeas(n,1) ) = 1;
end
end
%the imaginary part of the slack bus is not a state (phase angle = 0),
%therefore it can be deleted
H(:,3)=H(:,3)+3/sqrt(3)*H(:,4);
H(:,5)=H(:,5)-3/sqrt(3)*H(:,6);
% also the status of the missing phases can be deleted
delete_column = [[2,4,6],GridData.missing_line_phase];
H(:,delete_column)=[];
end
end
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment