Commit 006d1ce9 authored by Angioni's avatar Angioni

added recursive current DSSE and small updates to weight matrix

parent bbb210ca
function [Vmagn_status,Vph_status,Imagn_status,Iph_status,inj_status,P_post] = IRIDSSE_rec(Meas_vector_external,W,GridData,Test_SetUp,Imagn_status_p,Iph_status_p,inj_status_p,P_pre)
%1 the Jacobian and the gain matrix is calculated
[H] = Jacobian_m_IRIDSSE(GridData);
HW = H'*W;
G1 = HW*H;
%2 the input covariance, from the previous step, is scaled down, to convert
% them into pu values
for x = 1 : 11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 2*GridData.Lines_num
P_pre(:,x) = P_pre(:,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
P_pre(x,:) = P_pre(x,:) .* ([GridData.base_status; GridData.base_current*ones(2*GridData.Lines_num,1)] .^-1)';
end
%P_pre = diag(diag(P_pre));
% 3 the Jacobian of the previous state is an identity matrix
H_pre = eye(11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 1 + 2*GridData.Lines_num);
H_pre(11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 1,11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 1) = 0;
H_pre(11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 1,:) = [];
% 4 the previous matrix is combined with the one of the new measurements
R = inv(W);
R = blkdiag(R,P_pre);
W = inv(R);
H = [H;H_pre];
HW = H'*W;
G1 = HW*H;
%4 the output covariance is converted from pu to normal unit
P_post = inv(G1);
if GridData.inj_status == 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
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
P_post(x,:) = P_post(x,:) .* ([GridData.base_status; GridData.base_current*ones(2*GridData.Lines_num,1)] )';
end
end
%5 the input states are the currents, they are converted in real and
%imaginary parts
i_lineDQ = zeros(2*GridData.Lines_num,1);
for x = 1 : GridData.Lines_num
I_phasor = Imagn_status_p(x,1)*(cos(Iph_status_p(x,1)) + 1i * sin(Iph_status_p(x,1)));
i_lineDQ(2*x-1,1) = real(I_phasor);
i_lineDQ(2*x,1) = imag(I_phasor);
end
Meas_integrate = [inj_status_p;i_lineDQ];
old_states = zeros(11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 2*GridData.Lines_num,1);
if GridData.inj_status == 1 %in case we consider injection statuses
inj_status = zeros(11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload,1);
else
inj_status = 0;
end
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);
max_delta=10; %dummy inizialitazion
iteration=1;
% min_it = min(Ibase,GridData.base_voltage);
% Newton Rapson calculation of the state
while max_delta(iteration) > 1e-12 && iteration < Test_SetUp.limit2
iteration=iteration+1;
% 1) Calculate vector of measurements
%[Meas_vector] = build_Mvector(Vmagn_status,Vph_status,GridData,PowerData);
[Meas_vector] = calc_Mvector(Vmagn_status,Vph_status,GridData,Meas_vector_external);
Meas_vector = [Meas_vector;Meas_integrate];
% 2) Calculate h(x) vector (measurements calculated from the estimated
% states)
[hx] = calc_hx_IRIDSSE(Vmagn_status,Vph_status,Imagn_status,Iph_status,GridData,inj_status);
hx = [hx;old_states];
% build the residual vector
res = Meas_vector-hx;
% based on Jacobian, weight matrix and residual we can calculate the G
% matrix and therefore the delta to update the state
HWres = HW*res;
%dependent_columns = [[1:size(null(full(G1)),1)]',null(full(G1))]
delta = G1\HWres;
% if rcond(G1) < 1e-15
% conditioning_number = rcond(G1)
% end
% We update the state
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) ;
old_states(1:11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload ,1) = inj_status;
n = 11*GridData.DM.NGF+10*GridData.DM.NGS + 2*GridData.DM.Nload + 1;
end
Volts = Vmagn_status(:,1).*exp(sqrt(-1)*Vph_status(:,1));
Amps = Imagn_status(:,1).*exp(sqrt(-1)*Iph_status(:,1));
Volts(1,1) = Volts(1,1)+delta(n);
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);
n=n+1;
Amps(i,1)=Amps(i,1)+1i*delta(n);
Imagn_status(i,1)=abs(Amps(i,1));
Iph_status(i,1)=phase(Amps(i,1));
i_lineDQ(2*i-1,1) = Imagn_status(i,1)*cos(Iph_status(i,1));
i_lineDQ(2*i,1) = Imagn_status(i,1)*sin(Iph_status(i,1));
in1 = GridData.topology(2,i);
fin1 = GridData.topology(3,i);
Volts(fin1,1) = Volts(in1,1) - ( GridData.R1(1,i)+1i*GridData.X1(1,i))*Amps(i,1);
Vmagn_status(fin1,1) = abs(Volts(fin1,1));
Vph_status (fin1,1) = (phase(Volts(fin1,1)));
end %the state is updated
if GridData.inj_status == 1
old_states(11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 1 : 11*GridData.DM.NGF + 10*GridData.DM.NGS + 2*GridData.DM.Nload + 2*GridData.Lines_num,1) = i_lineDQ;
else
old_states(1 : 2*GridData.Lines_num,1) = i_lineDQ;
end
% We calculate the maximum delta so that we verify if we need to stop the
% Newton Rapson iterations
max_delta(iteration,1)=max(abs(delta));
end
Vph_status = wrapToPi(Vph_status);
if iteration == Test_SetUp.limit2
max_delta(end)
end
end
\ No newline at end of file
......@@ -23,7 +23,7 @@
function [H] = Jacobian_m_IRIDSSE(GridData)
if strcmp(GridData.type_of_model,'single_phase')==1
R1 = GridData.R1;
R1 = GridData.R1;
X1 = GridData.X1;
if GridData.inj_status == 0
H = zeros(GridData.MeasNum, 2 + 2*GridData.Lines_num);
......
......@@ -47,6 +47,7 @@ if strcmp(GridData.type_of_model,'single_phase') == 1
n = n +1;
LocationMeas(n,1) = x;
TypeMeas(n,1) = 100;
DelayMeas(n,1) = Combination_devices.xstatus_measure(2,x);
end
end
end
......@@ -66,9 +67,11 @@ if strcmp(GridData.type_of_model,'single_phase') == 1
n=n+1;
LocationMeas(n,1) = x;
TypeMeas(n,1) = 1;% 4 represent measurements of type voltage magnitude
DelayMeas(n,1) = Combination_devices.P_measure(2,x);
n=n+1;
LocationMeas(n,1) = x;
TypeMeas(n,1) = 2;% 4 represent measurements of type voltage phase angle
DelayMeas(n,1) = Combination_devices.Q_measure(2,x);
elseif Combination_devices.Pseudo_measure(1,x) == 1
rotPQ(1,1) = cos(PowerData.Vph(x))/(PowerData.Vmagn(x)^2);
......@@ -83,22 +86,25 @@ if strcmp(GridData.type_of_model,'single_phase') == 1
n=n+1;
LocationMeas(n,1) = x;
TypeMeas(n,1) = 1;% 4 represent measurements of type voltage magnitude
DelayMeas(n,1) = Combination_devices.Pseudo_measure(2,x);
n=n+1;
LocationMeas(n,1) = x;
TypeMeas(n,1) = 2;% 4 represent measurements of type voltage phase angle
DelayMeas(n,1) = Combination_devices.Pseudo_measure(2,x);
end
% voltage real part - here it is converted from the covariance of magnitude and phase angle
if Combination_devices.Vmagn_measure(1,x)==1 && Combination_devices.Vph_measure(1,x)==1
if x == 1 % we remove the phase measurement from the first bus
rotV(1,1) = cos(PowerData.Vph(x));
Rtemp = rotV*(Accuracy.Accuracy_Vmagn*PowerData.Vmagn(x))^2;
Rtemp = rotV*(Accuracy.Accuracy_Vmagn*PowerData.Vmagn(x))^2*rotV';
if Rtemp(1,1) < LM; Rtemp(1,1) = LM; end
R(n+1,n+1) = Rtemp;
W(n+1,n+1) = Rtemp^-1;
n=n+1;
LocationMeas(n,1) = x;
TypeMeas(n,1) = 3;% 3 represent measurements of type voltage magnitude
DelayMeas(n,1) = Combination_devices.Vmagn_measure(2,x);
else
rotV(1,1) = cos(PowerData.Vph(x));
......@@ -113,9 +119,11 @@ if strcmp(GridData.type_of_model,'single_phase') == 1
n=n+1;
LocationMeas(n,1) = x;
TypeMeas(n,1) = 3;% 4 represent measurements of type voltage magnitude
DelayMeas(n,1) = Combination_devices.Vmagn_measure(2,x);
n=n+1;
LocationMeas(n,1) = x;
TypeMeas(n,1) = 4;% 4 represent measurements of type voltage phase angle
DelayMeas(n,1) = Combination_devices.Vph_measure(2,x);
end
end
end
......@@ -138,29 +146,33 @@ if strcmp(GridData.type_of_model,'single_phase') == 1
n=n+1;
LocationMeas(n,1) = b;%in the 1st column we assign the line
TypeMeas(n,1) = 5;% 5 represent measurements of type current magnitude
DelayMeas(n,1) = Combination_devices.Imagn_measure(2,b);
n=n+1;
LocationMeas(n,1)=b;%in the 1st column we assign the line
TypeMeas(n,1)=6;% 6 represent measurements of type current ph angle
DelayMeas(n,1) = Combination_devices.Iph_measure(2,b);
end
if Combination_devices.Pflow_measure(1,b)~=0 %active power flow measurement
Rtemp = (Accuracy.Accuracy_Pflow * PowerData.Pflow(b)^2 )^2;
Rtemp = (Accuracy.Accuracy_Pflow * PowerData.Pflow(b))^2;
if Rtemp(1,1) < LM; Rtemp(1,1) = LM; end
R(n+1,n+1) = Rtemp;
W(n+1,n+1) = Rtemp^-1;
n=n+1;
LocationMeas(n,1) = b;%in the 1st column we assign the line
TypeMeas(n,1) = 7; % 7 represent measurements of type power active flow
DelayMeas(n,1) = Combination_devices.Pflow_measure(2,b);
end
if Combination_devices.Qflow_measure(1,b)~=0 %reactive power flow measurement
Rtemp = (Accuracy.Accuracy_Qflow * PowerData.Qflow(b)^2)^2;
Rtemp = (Accuracy.Accuracy_Qflow * PowerData.Qflow(b))^2;
if Rtemp(1,1) < LM; Rtemp(1,1) = LM; end
R(n+1,n+1) = Rtemp;
W(n+1,n+1) = Rtemp^-1;
n=n+1;
LocationMeas(n,1) = b;%in the 1st column we assign the line
TypeMeas(n,1) = 8; % 8 represent measurements of type power reactive flow
DelayMeas(n,1) = Combination_devices.Qflow_measure(2,b);
end
end
......@@ -183,6 +195,7 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
LocationMeas(n,1) = x;
TypeMeas(n,1) = 1;% 1 represent measurements of type active power
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.P_measure(2,x);
elseif Combination_devices.Pseudo_measure(1,x) == 1 %pseudo measurements, but we do not consider pseudo measurements at the slack bus
Rtemp = (Accuracy.Accuracy_pseudo*PowerData.Pinj(f,x))^2;
......@@ -193,6 +206,7 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
LocationMeas(n,1) = x;
TypeMeas(n,1) = 1;% 1 represent measurements of type active power
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.Pseudo_measure(2,x);
end
if Combination_devices.Q_measure(1,x)==1%reactive power measurement
......@@ -204,7 +218,7 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
LocationMeas(n,1) = x;
TypeMeas(n,1) = 2;% 2 represent measurements of type reactive power
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.Q_measure(2,x);
elseif Combination_devices.Pseudo_measure(1,x) == 1 %pseudo measurements, but we do not consider pseudo measurements at the slack bus
Rtemp = (Accuracy.Accuracy_pseudo*PowerData.Qinj(f,x))^2;
......@@ -215,6 +229,7 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
LocationMeas(n,1) = x;
TypeMeas(n,1) = 2;% 2 represent measurements of type reactive power
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.Pseudo_measure(2,x);
end
% voltage real part - here it is converted from the covariance of magnitude and phase angle
......@@ -229,6 +244,7 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
LocationMeas(n,1) = x;
TypeMeas(n,1) = 3;% 3 represent measurements of type voltage magnitude
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.Vmagn_measure(2,x);
else
rotV(1,1) = cos(PowerData.Vph(f,x));
......@@ -246,10 +262,12 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
LocationMeas(n,1) = x;
TypeMeas(n,1) = 3;% 4 represent measurements of type voltage magnitude
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.Vmagn_measure(2,x);
n=n+1;
LocationMeas(n,1) = x;
TypeMeas(n,1) = 4;% 4 represent measurements of type voltage phase angle
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.Vph_measure(2,x);
end
end
end
......@@ -275,11 +293,12 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
LocationMeas(n,1) = b;
TypeMeas(n,1) = 5;% 5 represent measurements of type current magnitude
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.Imagn_measure(2,b);
n=n+1;
LocationMeas(n,1)=b;
TypeMeas(n,1)=6;% 6 represent measurements of type current ph angle
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.Iph_measure(2,b);
end
if Combination_devices.Pflow_measure(1,b)~=0
Rtemp = (Accuracy.Accuracy_Pflow*PowerData.Pflow(f,b))^2;
......@@ -290,6 +309,7 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
LocationMeas(n,1) = b;
TypeMeas(n,1) = 7;% 1 represent measurements of type active power
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.Pflow_measure(2,b);
end
if Combination_devices.Qflow_measure(1,b)~=0
......@@ -301,6 +321,7 @@ elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridD
LocationMeas(n,1) = b;
TypeMeas(n,1) = 8;% 1 represent measurements of type active power
PhaseMeas(n,1) = f;
DelayMeas(n,1) = Combination_devices.Qflow_measure(2,b);
end
end
end
......@@ -310,4 +331,5 @@ end
GridData.MeasNum = n;
GridData.TypeMeas = TypeMeas;
GridData.LocationMeas = LocationMeas;
GridData.DelayMeas = DelayMeas;
end
\ No newline at end of file
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