Commit 1b1423aa authored by Andrea Angioni's avatar Andrea Angioni

minor updates

parent 388998ae
......@@ -27,20 +27,21 @@ function [Vmagn_status,Vph_status,Imagn_status,Iph_status,inj_status,P_post]=I
HW = H'*W;
G1 = HW*H;%gain matrix
G2 = G1;
P_post = inv(G2); %estimated state error covariance
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,:)=[];
P_post(:,11*GridData.DM.NGF+10*GridData.DM.NGS+2*GridData.DM.Nload+1)=[];
P_post(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);
P_post(:,x) = P_post(:,x) .* ([GridData.base_status; GridData.base_current*ones(2*GridData.Lines_num,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)';
P_post(x,:) = P_post(x,:) .* ([GridData.base_status; GridData.base_current*ones(2*GridData.Lines_num,1)] )';
end
end
P_post = inv(G2); %estimated state error covariance
%we initialize the state
%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);
......@@ -54,8 +55,8 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
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);
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
......@@ -67,11 +68,11 @@ iteration=1;
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
[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) ;
......@@ -93,22 +94,22 @@ while max_delta > 1e-12 && iteration < Test_SetUp.limit2
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
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);
......@@ -150,7 +151,7 @@ while max_delta > 1e-12 && iteration < Test_SetUp.limit2
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));
max_delta=max(abs(delta));
end
Vph_status = wrapToPi(Vph_status);
......
......@@ -51,7 +51,7 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
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);
inj_status = zeros(11*GridData.DM.NGF+10*GridData.DM.NGS + 2*GridData.DM.Nload,1);
else
inj_status = 0;
end
......@@ -70,7 +70,7 @@ while max_delta > 1e-12 && iteration < Test_SetUp.limit2
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) ;
inj_status = inj_status + delta(1:11*GridData.DM.NGF+10*GridData.DM.NGS + 2*GridData.DM.Nload) ;
n = 11*GridData.DM.NGF+10*GridData.DM.NGS + 2*GridData.DM.Nload + 1;
end
......
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