Commit c69304ab authored by Angioni's avatar Angioni

populated the code from private repository

parents
This diff is collapsed.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Current(I) Real Imaginary Distribution System State Estimator (VRIDSSE)
%
% @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 [Vmagn_status,Vph_status,Imagn_status,Iph_status,inj_status,P_post]=IRIDSSE...
(Meas_vector_external,W,GridData,Test_SetUp)
% Jacobian matrix
[H] = Jacobian_m_IRIDSSE(GridData);
HW = H'*W;
G1 = HW*H;%gain matrix
G2 = G1;
if GridData.inj_status == 1
G2(:,11*GridData.DM.NGF+10*GridData.DM.NGS+2*GridData.DM.Nload+1)=[];
G2(11*GridData.DM.NGF+10*GridData.DM.NGS+2*GridData.DM.Nload+1,:)=[];
for x = 1 : 11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 2*GridData.Lines_num
G2(:,x) = G2(:,x) .* ([GridData.base_status; GridData.base_current*ones(2*GridData.Lines_num,1)].^-1);
end
for x = 1 : 11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 2*GridData.Lines_num
G2(x,:) = G2(x,:) .* ([GridData.base_status; GridData.base_current*ones(2*GridData.Lines_num,1)].^-1)';
end
end
P_post = inv(G2); %estimated state error covariance
%we initialize the state
if strcmp(GridData.type_of_model,'single_phase')==1
Vmagn_status = ones(GridData.Nodes_num,1);
Vph_status = zeros(GridData.Nodes_num,1);
Imagn_status = zeros(GridData.Lines_num,1);
Iph_status = zeros(GridData.Lines_num,1);
elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridData.type_of_model,'three_phase_unbalance')==1
Vmagn_status=ones(3,GridData.Nodes_num);
Vph_status=zeros(3,GridData.Nodes_num);
Vph_status(2,:)=(4/3)*pi;
Vph_status(3,:)=(2/3)*pi;
Imagn_status = zeros(3,GridData.Lines_num);
Iph_status = zeros(3,GridData.Lines_num);
end
if GridData.inj_status == 1 %in case we consider dynamic measurements
inj_status = zeros(11*GridData.DM.NGF+10*GridData.DM.NGS + 2*GridData.DM.Nload,1);
else
inj_status = 0;
end
max_delta=10; %dummy inizialitazion
iteration=1;
% Newton Rapson calculation of the state
while max_delta > 1e-12 && iteration < Test_SetUp.limit2
iteration=iteration+1;
[Meas_vector] = calc_Mvector(Vmagn_status,Vph_status,GridData,Meas_vector_external);% 1) Update vector of measurements
[hx] = calc_hx_IRIDSSE(Vmagn_status,Vph_status,Imagn_status,Iph_status,GridData,inj_status);% 2) Calculate h(x) vector
res = Meas_vector - hx;% build the residual vector
HWres = HW*res;
delta = G1\HWres;
n = 1;
if GridData.inj_status == 1 %in case we consider dynamic measuremetns
inj_status = inj_status + delta(1:11*GridData.DM.NGF+10*GridData.DM.NGS + 2*GridData.DM.Nload,1) ;
n = 11*GridData.DM.NGF+10*GridData.DM.NGS + 2*GridData.DM.Nload + 1;
end
Volts = Vmagn_status.*exp(sqrt(-1)*Vph_status);
Amps = Imagn_status.*exp(sqrt(-1)*Iph_status);
%we update the status using delta
if strcmp(GridData.type_of_model,'single_phase')==1
Volts(1,1) = Volts(1,1)+delta(n,1);
Vmagn_status(1,1) = Volts(1,1);
Vph_status (1,1) = 0;
for i = 1:GridData.Lines_num
n=n+1;
Amps(i) = Amps(i) + delta(n,1);
n=n+1;
Amps(i,1) = Amps(i,1) + 1i * delta(n,1);
Imagn_status(i,1) = abs(Amps(i,1));
Iph_status (i,1) = phase(Amps(i,1));
end
Vnode = zeros(GridData.Nodes_num,1);
Vnode(1,1) = 1;
while min(Vnode) < 1
for b = 1 : GridData.Lines_num
if Vnode(GridData.topology(2,b),1) == 1
if Vnode(GridData.topology(3,b),1) == 0
x = GridData.topology(3,b);
Vnode(x,1) = Vnode(x,1) + 1;
Volts(x,1) = Volts(GridData.topology(2,b),1) - ( GridData.R1(1,b)+1i*GridData.X1(1,b))*Amps(b,1);
Vmagn_status(x,1) = abs(Volts(x,1));
Vph_status (x,1) = (phase(Volts(x,1)));
end
end
end
end
elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridData.type_of_model,'three_phase_unbalance')==1
Vreal(1,1)=real(Volts(1,1))+delta(1);
Volts(1,1)=Vreal(1,1);
Vreal(2,1)=real(Volts(2,1))+delta(2);
Volts(2,1)=Vreal(2,1)+1i*Vreal(2,1)*3/sqrt(3);
Vreal(3,1) = real(Volts(3,1))+delta(3);
Volts(3,1) = Vreal(3,1)+1i*Vreal(3,1)*-3/sqrt(3);
Vmagn_status (:,1) = abs(Volts(:,1));
Vph_status (:,1) = phase(Volts(:,1));
Amps = Imagn_status.*exp(sqrt(-1)*Iph_status);
n = 3;
for x = 1 : GridData.Lines_num
for f = 1 : 3
if GridData.present_line(f,x) > 0
n = n + 1;
Amps(f,x) = Amps(f,x) + delta(n);
n = n + 1;
Amps(f,x) = Amps(f,x) + 1i * delta(n);
Imagn_status(f,x) = abs (Amps(f,x));
Iph_status (f,x) = phase(Amps(f,x));
end
end
end
Vnode = zeros(GridData.Nodes_num,1);
Vnode(1,1) = 1;
while min(Vnode) < 1
for b = 1 : GridData.Lines_num
if Vnode(GridData.topology(2,b),1) == 1
if Vnode(GridData.topology(3,b),1) == 0
x = GridData.topology(3,b);
Vnode(x,1) = Vnode(x,1) + 1;
Volts(:,x) = Volts(:,GridData.topology(2,b)) - ( GridData.R1(:,3*b-2:3*b)+1i*GridData.X1(:,3*b-2:3*b))*Amps(:,b);
Vmagn_status(:,x) = abs(Volts(:,x));
Vph_status (:,x) = (phase(Volts(:,x)));
end
end
end
end
end
% We calculate the maximum delta so that we verify if we need to stop the Newton Rapson iterations
max_delta=max(abs(delta));
end
Vph_status = wrapToPi(Vph_status);
if iteration == Test_SetUp.limit2 %in case we reach the maximum number of iterations
max_delta
end
end
\ No newline at end of file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the Jacobian matrix for the Current Real Imaginary Distribution System State Estimator (IRIDSSE)
%
% @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(GridData)
if strcmp(GridData.type_of_model,'single_phase')==1
R1 = GridData.R1;
X1 = GridData.X1;
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;
X1 = GridData.X1;
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Voltage Real Imaginary Distribution System State Estimator (VRIDSSE)
%
% @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 [Vmagn_status,Vph_status,Imagn_status,Iph_status,inj_status,P_post] = VRIDSSE ...
(Meas_vector_external,W,GridData,Test_SetUp)
% Jacobian matrix
[H] = Jacobian_m_VRIDSSE(GridData);
HW = H'*W;
G1 = HW*H; %gain matrix
G2 = G1;
if GridData.inj_status
G2(:,11*GridData.DM.NGF+10*GridData.DM.NGS+2*GridData.DM.Nload+1)=[];
G2(11*GridData.DM.NGF+10*GridData.DM.NGS+2*GridData.DM.Nload+1,:)=[];
for x = 1 : 11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 2*GridData.Lines_num
G2(:,x) = G2(:,x) .* ([GridData.base_status; GridData.base_current*ones(2*GridData.Lines_num,1)].^-1);
end
for x = 1 : 11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 2*GridData.Lines_num
G2(x,:) = G2(x,:) .* ([GridData.base_status; GridData.base_current*ones(2*GridData.Lines_num,1)].^-1)';
end
end
P_post = inv(G2); %estimated state error covariance
%we initialize the state
if strcmp(GridData.type_of_model,'single_phase')==1
Vmagn_status=ones(GridData.Nodes_num,1);
Vph_status=zeros(GridData.Nodes_num,1);
elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridData.type_of_model,'three_phase_unbalance')==1
Vmagn_status=ones(3,GridData.Nodes_num);
Vph_status=zeros(3,GridData.Nodes_num);
Vph_status(2,:)=(4/3)*pi;
Vph_status(3,:)=(2/3)*pi;
end
if GridData.inj_status == 1 %in case we consider dynamic measuremetns
inj_status = zeros(size(PowerData.x_status,1),1);
else
inj_status = 0;
end
max_delta=10; %dummy inizialitazion
iteration=0;
% Newton Rapson calculation of the state
while max_delta > 1e-12 && iteration < Test_SetUp.limit2
iteration=iteration+1;
[Meas_vector] = calc_Mvector(Vmagn_status,Vph_status,GridData,Meas_vector_external);% 1) Update vector of measurements based on the voltage estimated
[hx] = calc_hx_VRIDSSE(Vmagn_status,Vph_status,GridData,inj_status); % 2) Calculate h(x) vector
res = Meas_vector-hx; % build the residual vector
HWres = HW*res;
delta = G1\HWres;
n = 1;
if GridData.inj_status == 1 %in case we consider dynamic status of the injection
inj_status = inj_status + delta(1:11*GridData.GridData.DM.NGF+10*GridData.GridData.DM.NGS + 2*GridData.GridData.DM.Nload) ;
n = 11*GridData.DM.NGF+10*GridData.DM.NGS + 2*GridData.DM.Nload + 1;
end
Volts = Vmagn_status.*exp(sqrt(-1)*Vph_status);
%we update the status using delta
if strcmp(GridData.type_of_model,'single_phase')==1
Vreal(1,1) = real(Volts(1,1))+delta(n); %at the slack bus we update only the magnitude because the phase angle is always zero;
Volts(1,1) = Vreal(1,1);
Vmagn_status(1,1)=abs(Volts(1,1));
Vph_status (1,1)=angle(Volts(1,1));
for x = 2:GridData.Nodes_num
n=n+1;
Volts(x)=Volts(x)+delta(n);
n=n+1;
Volts(x,1)=Volts(x,1)+1i*delta(n);
Vmagn_status(x,1)=abs(Volts(x,1));
Vph_status (x,1)=angle(Volts(x,1));
end %the state is updated
elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridData.type_of_model,'three_phase_unbalance')==1
Vreal(1,1)=real(Volts(1,1))+delta(1); %at the slack bus we update only the magnitude because the phase angle is always zero;
Volts(1,1)=Vreal(1,1);
Vmagn_status(1,1)=abs(Volts(1,1));
Vph_status (1,1)=angle(Volts(1,1));
Vreal(2,1)=real(Volts(2,1))+delta(2);
Volts(2,1)=Vreal(2,1)+1i*Vreal(2,1)*3/sqrt(3);
Vmagn_status(2,1)=abs(Volts(2,1));
Vph_status (2,1)=angle(Volts(2,1));
Vreal(3,1) = real(Volts(3,1))+delta(3);
Volts(3,1) = Vreal(3,1)+1i*Vreal(3,1)*-3/sqrt(3);
Vmagn_status (3,1) = abs(Volts(3,1));
Vph_status (3,1) = angle(Volts(3,1));
n = 3;
for x = 2 : GridData.Nodes_num
for f = 1 : 3
if GridData.present_node(f,x) > 0
n = n + 1;
Volts(f,x) = Volts(f,x) + delta(n);
n = n + 1;
Volts(f,x) = Volts(f,x) + 1i*delta(n);
end
Vmagn_status(f,x) = abs(Volts(f,x));
Vph_status (f,x) = angle(Volts(f,x));
end
end
end
% We calculate the maximum delta so that we verify if we need to stop the Newton Rapson iterations
max_delta=max(abs(delta));
end
%based on the estimated voltages we calculate also the currents
if strcmp(GridData.type_of_model,'single_phase')==1
I_branch = zeros(GridData.Lines_num,1);
Volts = Vmagn_status.*exp(sqrt(-1)*Vph_status);
for m=1:GridData.Lines_num
n_i=GridData.topology(2,m);
n_f=GridData.topology(3,m);
V_i=Volts(n_i);
V_f=Volts(n_f);
R=GridData.R2(m);
X=GridData.X2(m);
B=GridData.B1(m);
G=GridData.G1(m);
I_branch(m) = 0.5*(G+1i*B)*(V_i+V_f) + (R+1i*X)*(V_i-V_f);
end %branch currents %some problems using this type of calculation
Imagn_status = abs(I_branch);
Iph_status = phase(I_branch);
Vph_status = wrapToPi(Vph_status);
Iph_status = wrapToPi(Iph_status);
elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridData.type_of_model,'three_phase_unbalance')==1
Imagn_status=zeros(3,GridData.Lines_num);
Iph_status=zeros(3,GridData.Lines_num);
I_branch=zeros(3,GridData.Lines_num);
for m = 1 : GridData.Lines_num
n_i=GridData.topology(2,m);
n_f=GridData.topology(3,m);
V_i = Volts(:,n_i);
V_f = Volts(:,n_f);
R=GridData.R2(:,3*m-2:3*m);
X=GridData.X2(:,3*m-2:3*m);
B=GridData.B1(:,3*m-2:3*m);
G=GridData.G1(:,3*m-2:3*m);
I_branch(:,m) = 0.5*(G+1i*B)*(V_i+V_f) +(R+1i*X)*(V_i-V_f); %there is the current going to the admittances at beginnend and end of the line
Imagn_status(:,m)=abs(I_branch(:,m)); %amplitudes and phases of the branch currents are considered at the end of the line
Iph_status(:,m)=phase(I_branch(:,m));
end %branch currents %some problems using this type of calculation
Vph_status = wrapToPi(Vph_status);
Iph_status = wrapToPi(Iph_status);
end
if iteration == Test_SetUp.limit2 %in case we reach the maximum number of iterations
max_delta
end
end
\ No newline at end of file
This diff is collapsed.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Update the measurement vector based on the estimated voltages
%
% @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 [Meas_vector] = calc_Mvector(Vmagn_status,Vph_status,GridData,Meas_vector_external)
Meas_vector = Meas_vector_external;
if strcmp(GridData.type_of_model,'single_phase')==1
for n=1:GridData.MeasNum
if GridData.TypeMeas(n,1)==1 %real current
for m = 1 : GridData.MeasNum
if GridData.LocationMeas(m,1) == GridData.LocationMeas(n,1) && GridData.TypeMeas(m,1) == 2
S1 = (Meas_vector_external(n,1)+1i*Meas_vector_external(m,1));
end
end
V1 = (Vmagn_status(GridData.LocationMeas(n,1))*exp(sqrt(-1)*Vph_status(GridData.LocationMeas(n,1))));
I_load1 = conj(S1/V1);
Meas_vector(n,1)= + real(I_load1);
end
if GridData.TypeMeas(n,1)==2 %immaginary current
for m = 1 : GridData.MeasNum
if GridData.LocationMeas(m,1) == GridData.LocationMeas(n,1) && GridData.TypeMeas(m,1) == 1
S1 = (Meas_vector_external(m,1)+1i*Meas_vector_external(n,1));
end
end
Meas_vector(n,1)= - imag((1i*Meas_vector_external(n,1) + Meas_vector_external(n-1,1))/...
(Vmagn_status(GridData.LocationMeas(n,1))*exp(sqrt(-1)*Vph_status(GridData.LocationMeas(n,1)))));
end
end
elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridData.type_of_model,'three_phase_unbalance')==1
for n = 1 : GridData.MeasNum
if GridData.TypeMeas(n,1)==1 %real current
for m = 1 : GridData.MeasNum
if GridData.LocationMeas(m,1) == GridData.LocationMeas(n,1) && GridData.TypeMeas(m,1) == 2 && GridData.PhaseMeas(n,1) == GridData.PhaseMeas(m,1)
S1 = (Meas_vector_external(n,1)+1i*Meas_vector_external(m,1));
end
end
V1 = Vmagn_status(GridData.PhaseMeas(n,1),GridData.LocationMeas(n,1))*exp(sqrt(-1)*Vph_status(GridData.PhaseMeas(n,1),GridData.LocationMeas(n,1)));
I_load1 = conj(S1/V1);
Meas_vector(n,1)= + real( I_load1 );
end
if GridData.TypeMeas(n,1)==2 %immaginary current
for m = 1 : GridData.MeasNum
if GridData.LocationMeas(m,1) == GridData.LocationMeas(n,1) && GridData.TypeMeas(m,1) == 1 && GridData.PhaseMeas(n,1) == GridData.PhaseMeas(m,1)
S1 = (Meas_vector_external(m,1)+1i*Meas_vector_external(n,1));
end
end
V1 = Vmagn_status(GridData.PhaseMeas(n,1),GridData.LocationMeas(n,1))*exp(sqrt(-1)*Vph_status(GridData.PhaseMeas(n,1),GridData.LocationMeas(n,1)));
I_load1 = conj(S1/V1);
Meas_vector(n,1)= imag( I_load1 );
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to generate the measurement vector from PowerData structure
%
% @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 [Meas_vector] = calc_Mvector_external(GridData,PowerData)
Meas_vector=zeros(GridData.MeasNum,1);
inj_status = 0;
if strcmp(GridData.type_of_model,'single_phase')==1
for n=1:GridData.MeasNum
if GridData.TypeMeas(n,1)==1 %active power
Meas_vector(n,1)= PowerData.Pinj(GridData.LocationMeas(n,1));
end
if GridData.TypeMeas(n,1)==2 %reactive power
Meas_vector(n,1)= PowerData.Qinj(GridData.LocationMeas(n,1));
end
if GridData.TypeMeas(n,1)==3 %voltage magnitude - translated to Vreal
Meas_vector(n,1)= + real(PowerData.Vmagn(GridData.LocationMeas(n,1))*exp(sqrt(-1)*PowerData.Vph(GridData.LocationMeas(n,1))));
end
if GridData.TypeMeas(n,1)==4 %voltage phase angle - translated to Vimag
Meas_vector(n,1)= + imag(PowerData.Vmagn(GridData.LocationMeas(n,1))*exp(sqrt(-1)*PowerData.Vph(GridData.LocationMeas(n,1))));
end
if GridData.TypeMeas(n,1)==5 %current magnitude - translated to Ireal
Meas_vector(n,1)= + real((PowerData.Imagn(GridData.LocationMeas(n,1)))*exp(sqrt(-1)*PowerData.Iph(GridData.LocationMeas(n,1))));
end
if GridData.TypeMeas(n,1)==6 %current phase - translated to Iimag
Meas_vector(n,1)= + imag((PowerData.Imagn(GridData.LocationMeas(n,1)))*exp(sqrt(-1)*PowerData.Iph(GridData.LocationMeas(n,1))));
end
if GridData.TypeMeas(n,1)==7 %pflow - translated to Ireal
Meas_vector(n,1)= PowerData.Pflow(GridData.LocationMeas(n,1));
end
if GridData.TypeMeas(n,1)==8 %qflow - - translated to Ireal
Meas_vector(n,1)=PowerData.Qflow(GridData.LocationMeas(n,1));
end
if GridData.TypeMeas(n,1)==100 %injection states
inj_status = inj_status +1;
Meas_vector(n,1) = PowerData.x_status(inj_status,1) ;
end
end
elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridData.type_of_model,'three_phase_unbalance')==1
for n = 1 : GridData.MeasNum
if GridData.TypeMeas(n,1)==1 %active power
Meas_vector(n,1)= PowerData.Pinj(GridData.PhaseMeas(n,1),GridData.LocationMeas(n,1));