MATLAB|【鲁棒电力系统状态估计】基于投影统计的电力系统状态估计

部分代码：
zdata = zconv(nbus); % Get Conventional Measurement data..
[bsh g b] = line_mat_func(nbus); % Get conductance and susceptance matrix 
type = zdata(:,2); 
% Type of measurement,
% type =1 voltage magnitude p.u
% type =2 Voltage phase angle in degree
% type =3 Real power injections
% type =4 Reactive power injection
% type =5 Real power flow 
% type =6 Reactive power flow 
z = zdata(:,3); % Measurement values
Z=z;% for ploting figures
fbus = zdata(:,4); % From bus
tbus = zdata(:,5); % To bus
Ri = diag(zdata(:,6)); % Measurement Error Covariance matrix
e = ones(nbus,1); % Initialize the real part of bus voltages
f = zeros(nbus,1);% Initialize the imaginary part of bus voltages
E = [f;e];  % State Vector comprising of imaginary and real part of voltage
G = real(ybus);
B = imag(ybus);
ei = find(type == 1); % Index of voltage magnitude measurements..
fi = find(type == 2); % Index of voltage angle measurements..
ppi = find(type == 3); % Index of real power injection measurements..
qi = find(type == 4); % Index of reactive power injection measurements..
pf = find(type == 5); % Index of real power flow measurements..
qf = find(type == 6); % Index of reactive power flow measurements..
Vm=z(ei);
Thm=z(fi);
z(ei)=Vm.*cosd(Thm); % converting voltage from polar to Cartesian
z(fi)=Vm.*sind(Thm);
nei = length(ei); % Number of Voltage measurements(real)
nfi = length(fi); % Number of Voltage measurements(imaginary)
npi = length(ppi); % Number of Real Power Injection measurements..
nqi = length(qi); % Number of Reactive Power Injection measurements..
npf = length(pf); % Number of Real Power Flow measurements..
nqf = length(qf); % Number of Reactive Power Flow measurements..
nm=nei+nfi+npi+nqi+npf+nqf; % total number of measurements
% robust parameters
tol=1;
maxiter=30;% maximal iteration for iteratively reweighted least squares (IRLS) algorithm
c=1.5; % for Huber-estimator
bm=mad_factor(nm); % correction factor to achieve unbiasness under Gaussian measurement noise
%%%%%%% GM-estimator%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% flat initialization
    iter=1;
    s=1;
%% For the GM-estimator to be able to handle two conforming outliers located on the same bus
%% the local redundancy must be large enough
%% add outliers %%