Commit cc889944 authored by Andrea Angioni's avatar Andrea Angioni

ph.angle is also included as state in slack bus

parent a7194b80
......@@ -83,9 +83,21 @@ while max_delta > 1e-12 && iteration < Test_SetUp.limit2
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;
if isfield(GridData,'rm_column') == 1
if GridData.rm_column == 1
Volts(1,1)=Volts(1,1)+delta(n);
end
Volts(1,1)=Volts(1,1)+delta(n);
n=n+1;
Volts(1,1)=Volts(1,1)+1i*delta(n);
else
Volts(1,1)=Volts(1,1)+delta(n);
end
Vmagn_status(1,1) = abs(Volts(1,1));
Vph_status (1,1) = angle(Volts(1,1));
for i = 1:GridData.Lines_num
n=n+1;
Amps(i) = Amps(i) + delta(n,1);
......
......@@ -97,8 +97,13 @@ for n = 1 : GridData.MeasNum
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)=[];
if isfield(GridData,'rm_column') == 1
if GridData.rm_column == 1
H(:,dim_state_inv + 2)=[];
end
else
H(:,dim_state_inv + 2)=[];
end
elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridData.type_of_model,'three_phase_unbalance')==1
R1 = GridData.R1;
......
......@@ -138,7 +138,13 @@ if strcmp(GridData.type_of_model,'single_phase')==1
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)=[];
if isfield(GridData,'rm_column') == 1
if GridData.rm_column == 1
H(:,dim_state_inv + 2)=[];
end
else
H(:,dim_state_inv + 2)=[];
end
elseif strcmp(GridData.type_of_model,'three_phase_sequence')==1 || strcmp(GridData.type_of_model,'three_phase_unbalance')==1
......
......@@ -77,10 +77,20 @@ while max_delta > 1e-12 && iteration < Test_SetUp.limit2
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);
if isfield(GridData,'rm_column') == 1
if GridData.rm_column == 1
Volts(1,1)=Volts(1,1)+delta(n);
end
Volts(1,1)=Volts(1,1)+delta(n);
n=n+1;
Volts(1,1)=Volts(1,1)+1i*delta(n);
else
Volts(1,1)=Volts(1,1)+delta(n);
end
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);
......
......@@ -74,11 +74,18 @@ if strcmp(GridData.type_of_model,'single_phase') == 1
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);
rotPQ(1,2) = sin(PowerData.Vph(x))/(PowerData.Vmagn(x)^2);
rotPQ(2,1) = sin(PowerData.Vph(x))/(PowerData.Vmagn(x)^2);
rotPQ(2,2) =-cos(PowerData.Vph(x))/(PowerData.Vmagn(x)^2);
% rotPQ(1,1) = cos(PowerData.Vph(x))/(PowerData.Vmagn(x)^2);
% rotPQ(1,2) = sin(PowerData.Vph(x))/(PowerData.Vmagn(x)^2);
% rotPQ(2,1) = sin(PowerData.Vph(x))/(PowerData.Vmagn(x)^2);
% rotPQ(2,2) =-cos(PowerData.Vph(x))/(PowerData.Vmagn(x)^2);
rotPQ(1,1) = cos(PowerData.Vph(x))/(PowerData.Vmagn(x));
rotPQ(1,2) = sin(PowerData.Vph(x))/(PowerData.Vmagn(x));
rotPQ(2,1) = sin(PowerData.Vph(x))/(PowerData.Vmagn(x));
rotPQ(2,2) =-cos(PowerData.Vph(x))/(PowerData.Vmagn(x));
Rtemp = rotPQ*[(Accuracy.Accuracy_pseudo*PowerData.Pinj(x))^2, 0; 0 (Accuracy.Accuracy_pseudo*PowerData.Qinj(x))^2]*rotPQ';
if Rtemp(1,1) < LM; Rtemp(1,1) = LM; end
if Rtemp(2,2) < LM; Rtemp(2,2) = LM; end
R(n+1:n+2,n+1:n+2) = Rtemp;
......@@ -95,16 +102,51 @@ if strcmp(GridData.type_of_model,'single_phase') == 1
% 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*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);
if x == 1 % we remove the phase measurement from the first bus
if isfield(GridData,'rm_column') == 1
if GridData.rm_column == 0
rotV(1,1) = cos(PowerData.Vph(x));
rotV(1,2) = - sin(PowerData.Vph(x))*PowerData.Vmagn(x) ;
rotV(2,1) = sin(PowerData.Vph(x)) ;
rotV(2,2) = cos(PowerData.Vph(x))*PowerData.Vmagn(x);
Rtemp = rotV*[(Accuracy.Accuracy_Vmagn*PowerData.Vmagn(x))^2, 0; 0 (Accuracy.Accuracy_Vph)^2]*rotV';
if Rtemp(1,1) < LM; Rtemp(1,1) = LM; end
if Rtemp(2,2) < LM; Rtemp(2,2) = LM; end
R(n+1:n+2,n+1:n+2) = Rtemp;
W(n+1:n+2,n+1:n+2)= Rtemp^-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);
else
rotV(1,1) = cos(PowerData.Vph(x));
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);
end
else
rotV(1,1) = cos(PowerData.Vph(x));
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);
end
else
rotV(1,1) = cos(PowerData.Vph(x));
......
......@@ -99,8 +99,8 @@ if strcmp(GridData.type_of_model,'single_phase')==1
var_err_rectangular_mat_IRI = diag(Rx_rectangular_IRI);
std_err_rectangular_mat_IRI = sqrt(var_err_rectangular_mat_IRI);
rotV = calc_rotV( GridData,PowerData,use_pu);
rotI = calc_rotI( GridData,PowerData,use_pu);
[rotV,rm_column] = calc_rotV( GridData,PowerData,use_pu);
[rotI,rm_column] = calc_rotI( GridData,PowerData,use_pu);
%we apply the rotor
Rx_polar_VRI = rotV*Rx_rectangular_VRI*rotV';
Rx_polar_IRI = rotI*Rx_rectangular_IRI*rotI';
......@@ -109,20 +109,29 @@ if strcmp(GridData.type_of_model,'single_phase')==1
std_err_polar_mat_IRI = sqrt(var_err_polar_mat_IRI);
std_err_polar_mat_VRI = sqrt(var_err_polar_mat_VRI);
% standard deviation of the error in polar format
std_err_Vimag_mat(1) = 0;
std_err_Vreal_mat(1) = GridData.base_voltage*std_err_rectangular_mat_VRI(1);
std_err_Vreal_mat(2:GridData.Nodes_num,1) = GridData.base_voltage*std_err_rectangular_mat_VRI(2:2:2*GridData.Nodes_num-1);
std_err_Vimag_mat(2:GridData.Nodes_num,1) = GridData.base_voltage*std_err_rectangular_mat_VRI(3:2:2*GridData.Nodes_num-1);
if rm_column == 1
std_err_Vimag_mat(1) = 0;
else
std_err_Vimag_mat(1) = GridData.base_voltage*std_err_rectangular_mat_VRI(2);
end
std_err_Vreal_mat(2:GridData.Nodes_num,1) = GridData.base_voltage*std_err_rectangular_mat_VRI(2+(1-rm_column):2:2*GridData.Nodes_num-1+(1-rm_column));
std_err_Vimag_mat(2:GridData.Nodes_num,1) = GridData.base_voltage*std_err_rectangular_mat_VRI(3+(1-rm_column):2:2*GridData.Nodes_num-1+(1-rm_column));
std_err_Vmagnitude_mat(1) = GridData.base_voltage*std_err_polar_mat_VRI(1);
std_err_Vphase_mat(1) = 0;
std_err_Vmagnitude_mat(2:GridData.Nodes_num,1) = GridData.base_voltage*std_err_polar_mat_VRI(2:2:2*GridData.Nodes_num-1);
std_err_Vphase_mat(2:GridData.Nodes_num,1) = std_err_polar_mat_VRI(3:2:2*GridData.Nodes_num-1);
if rm_column == 1
std_err_Vphase_mat(1) = 0;
else
std_err_Vphase_mat(1) = std_err_rectangular_mat_VRI(2);
end
std_err_Vmagnitude_mat(2:GridData.Nodes_num,1) = GridData.base_voltage*std_err_polar_mat_VRI(2+(1-rm_column):2:2*GridData.Nodes_num-1+(1-rm_column));
std_err_Vphase_mat(2:GridData.Nodes_num,1) = std_err_polar_mat_VRI(3+(1-rm_column):2:2*GridData.Nodes_num-1+(1-rm_column));
std_err_Imagnitude_mat = GridData.base_current*std_err_polar_mat_IRI(2:2:2*GridData.Lines_num+1);
std_err_Iphase_mat = std_err_polar_mat_IRI(3:2:2*GridData.Lines_num+1);
std_err_Ireal_mat = GridData.base_current*std_err_rectangular_mat_IRI(2:2:2*GridData.Lines_num+1);
std_err_Iimag_mat = GridData.base_current*std_err_rectangular_mat_IRI(3:2:2*GridData.Lines_num+1);
std_err_Imagnitude_mat = GridData.base_current*std_err_polar_mat_IRI(2+(1-rm_column):2:2*GridData.Lines_num+1+(1-rm_column));
std_err_Iphase_mat = std_err_polar_mat_IRI(3+(1-rm_column):2:2*GridData.Lines_num+1+(1-rm_column));
std_err_Ireal_mat = GridData.base_current*std_err_rectangular_mat_IRI(2+(1-rm_column):2:2*GridData.Lines_num+1+(1-rm_column));
std_err_Iimag_mat = GridData.base_current*std_err_rectangular_mat_IRI(3+(1-rm_column):2:2*GridData.Lines_num+1+(1-rm_column));
std_err_Iri = GridData.base_current*std_err_rectangular_mat_IRI;
else
......@@ -183,9 +192,9 @@ if strcmp(GridData.type_of_model,'single_phase')==1
Rx_rectangular_IRI_rmsemp(2*(x-1)+1,2*(x-1)+1) = rms_err_Ireal_emp(x,1)^2;
Rx_rectangular_IRI_rmsemp(2*(x-1)+2,2*(x-1)+2) = rms_err_Iimag_emp(x,1)^2;
end
rotI = calc_rotI( GridData,PowerData,0);
Rx_polar_IRI_stdemp = rotI(2:end,2:end)*Rx_rectangular_IRI_stdemp*rotI(2:end,2:end)';
Rx_polar_IRI_rmsemp = rotI(2:end,2:end)*Rx_rectangular_IRI_rmsemp*rotI(2:end,2:end)';
[rotI,rm_column] = calc_rotI( GridData,PowerData,0);
Rx_polar_IRI_stdemp = rotI(2+(1-rm_column):end,2+(1-rm_column):end)*Rx_rectangular_IRI_stdemp*rotI(2+(1-rm_column):end,2+(1-rm_column):end)';
Rx_polar_IRI_rmsemp = rotI(2+(1-rm_column):end,2+(1-rm_column):end)*Rx_rectangular_IRI_rmsemp*rotI(2+(1-rm_column):end,2+(1-rm_column):end)';
for x = 1 : GridData.Lines_num
std_err_Imagnitude_empmat(x,1)= sqrt(Rx_polar_IRI_stdemp(2*(x-1)+1,2*(x-1)+1));
std_err_Iphase_empmat(x,1) = sqrt(Rx_polar_IRI_stdemp(2*(x-1)+2,2*(x-1)+2));
......
......@@ -20,10 +20,14 @@
% 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 [ rotI ] = calc_rotI( GridData,PowerData,use_pu)
function [rotI,rm_column] = calc_rotI( GridData,PowerData,use_pu)
% R = GridData.R2; %this are actually the inverses of the resistances
% X = GridData.X2; %this are actually the inverses of the reactances
rm_column = 1;
if isfield(GridData,'rm_column') == 1
rm_column = GridData.rm_column;
else
rm_column = 1;
end
if use_pu == 0
PowerDataVmagn = PowerData.Vmagn * GridData.base_voltage;
PowerDataImagn = PowerData.Imagn* GridData.base_current;
......
......@@ -20,8 +20,14 @@
% 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 [ rotV ] = calc_rotV( GridData,PowerData,use_pu)
rm_column = 1;
function [ rotV,rm_column] = calc_rotV( GridData,PowerData,use_pu)
if isfield(GridData,'rm_column') == 1
rm_column = GridData.rm_column;
else
rm_column = 1;
end
if rm_column == 1
rotV = zeros(2*GridData.Nodes_num-1,2*GridData.Nodes_num-1);
else
......@@ -43,7 +49,7 @@ for x = 1 : GridData.Nodes_num
rotV(2*(x-1)+1-1,2*(x-1)+2-1) = sin(PowerData.Vph(x,1));
rotV(2*(x-1)+2-1,2*(x-1)+1-1) = - sin(PowerData.Vph(x,1))/PowerDataVmagn(x,1);
rotV(2*(x-1)+2-1,2*(x-1)+2-1) = cos(PowerData.Vph(x,1))/PowerDataVmagn(x,1);
end
end
else
rotV(2*(x-1)+1,2*(x-1)+1) = cos(PowerData.Vph(x,1));
rotV(2*(x-1)+1,2*(x-1)+2) = sin(PowerData.Vph(x,1));
......
......@@ -28,11 +28,11 @@ close all
results_DSSE_WLS2=struct;results_DSSE_WLS2.err_Vmagnitude=[];results_DSSE_WLS2.err_Vphase=[];results_DSSE_WLS2.Vmagn_status=[];results_DSSE_WLS2.Vph_status=[];results_DSSE_WLS2.Imagn_status=[];results_DSSE_WLS2.Iph_status=[];results_DSSE_WLS2.Vreal_status=[];results_DSSE_WLS2.Vimag_status=[];results_DSSE_WLS2.Ireal_status=[];results_DSSE_WLS2.Iimag_status=[];results_DSSE_WLS2.Vmagn_true=[];results_DSSE_WLS2.Vph_true=[];results_DSSE_WLS2.Imagn_true=[];results_DSSE_WLS2.Iph_true=[];results_DSSE_WLS2.Vreal_true=[];results_DSSE_WLS2.Vimag_true=[];results_DSSE_WLS2.Ireal_true=[];results_DSSE_WLS2.Iimag_true=[];results_DSSE_WLS2.err_Imagnitude=[];results_DSSE_WLS2.err_Iphase=[];results_DSSE_WLS2.err_Ireal=[];results_DSSE_WLS2.err_Iimag=[];results_DSSE_WLS2.err_Vreal=[];results_DSSE_WLS2.err_Vimag=[];results_DSSE_WLS2.err_Vmagnitude_in =[];results_DSSE_WLS2.err_Vphase_in=[];results_DSSE_WLS2.err_Imagnitude_in =[];results_DSSE_WLS2.err_Iphase_in =[];results_DSSE_WLS2.err_Ireal_in =[];results_DSSE_WLS2.err_Iimag_in =[];results_DSSE_WLS2.Vmagn_measured =[];results_DSSE_WLS2.Vph_measured =[];
results_DSSE_WLS=struct;results_DSSE_WLS.err_Vmagnitude=[];results_DSSE_WLS.err_Vphase=[];results_DSSE_WLS.Vmagn_status=[];results_DSSE_WLS.Vph_status=[];results_DSSE_WLS.Imagn_status=[];results_DSSE_WLS.Iph_status=[];results_DSSE_WLS.Vreal_status=[];results_DSSE_WLS.Vimag_status=[];results_DSSE_WLS.Ireal_status=[];results_DSSE_WLS.Iimag_status=[];results_DSSE_WLS.Vmagn_true=[];results_DSSE_WLS.Vph_true=[];results_DSSE_WLS.Imagn_true=[];results_DSSE_WLS.Iph_true=[];results_DSSE_WLS.Vreal_true=[];results_DSSE_WLS.Vimag_true=[];results_DSSE_WLS.Ireal_true=[];results_DSSE_WLS.Iimag_true=[];results_DSSE_WLS.err_Imagnitude=[];results_DSSE_WLS.err_Iphase=[];results_DSSE_WLS.err_Ireal=[];results_DSSE_WLS.err_Iimag=[];results_DSSE_WLS.err_Vreal=[];results_DSSE_WLS.err_Vimag=[];results_DSSE_WLS.err_Vmagnitude_in =[];results_DSSE_WLS.err_Vphase_in=[];results_DSSE_WLS.err_Imagnitude_in =[];results_DSSE_WLS.err_Iphase_in =[];results_DSSE_WLS.err_Ireal_in =[];results_DSSE_WLS.err_Iimag_in =[];results_DSSE_WLS.Vmagn_measured =[];results_DSSE_WLS.Vph_measured =[];
%select type of model among 'single_phase' 'three_phase_sequence' 'three_phase_unbalance'
type_of_model = 'three_phase_unbalance';
type_of_model = 'single_phase';
[GridData] = generate_GridData(type_of_model);%in this function the static model is generated
[PowerData] = generate_PowerData(GridData); %in this function the power flow is run in order to obtain the reference values to test the state estimator
[Test_SetUp,Combination_devices,Accuracy] = generate_DSSEConfData(GridData);%here the test configuration data are set: measurement devices location and accuracy
GridData.rm_column = 0;%in this case the phase angle at 1st bus is also considered as state
[W,GridData,R] = Weight_m(GridData,PowerData,Combination_devices,Accuracy); %weight and covariance matrix of the state estimator
[Meas_true_vector] = calc_Mvector_external(GridData,PowerData); %vector with reference values of the measurements
......@@ -41,9 +41,9 @@ for z = 1 : Test_SetUp.N_MC %in this for loop the Monte Carlo tests are run
Meas_vector_external = mvnrnd(Meas_true_vector, R)'; %the measurements are corrupted following the covariance matrix
[Vmagn_status_WLS,Vph_status_WLS] = VRIDSSE(Meas_vector_external,W,GridData,Test_SetUp); %voltage state estimator
[Vmagn_status_WLS2,Vph_status_WLS2,Imagn_status_WLS2,Iph_status_WLS2] = IRIDSSE(Meas_vector_external,W,GridData,Test_SetUp);%current state estimator
[Vmagn_status_WLS2,Vph_status_WLS2,Imagn_status_WLS2,Iph_status_WLS2] = IRIDSSE(Meas_vector_external,W,GridData,Test_SetUp);%current state estimator
[results_DSSE_WLS] = Data_Output(Vmagn_status_WLS,Vph_status_WLS,results_DSSE_WLS,GridData,PowerData); %collects the metrics of the error
[results_DSSE_WLS2] = Data_Output(Vmagn_status_WLS2,Vph_status_WLS2,results_DSSE_WLS2,GridData,PowerData);
......
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