Description of state

0001 function [V, converged, i] = state_est(branch, Ybus, Yf, Yt, Sbus, V0, ref, pv, pq, mpopt)
0002 %STATE_EST  Solves a state estimation problem.
0003 %   [V, CONVERGED, I] = STATE_EST(BRANCH, YBUS, YF, YT, SBUS, ...
0004 %                                   V0, REF, PV, PQ, MPOPT)
0005 %   State estimator (under construction) based on code from James S. Thorp.
0006 
0007 %   MATPOWER
0008 %   Copyright (c) 1996-2015 by Power System Engineering Research Center (PSERC)
0009 %   by Ray Zimmerman, PSERC Cornell
0010 %   based on code by James S. Thorp, June 2004
0011 %
0012 %   $Id: state_est.m 2644 2015-03-11 19:34:22Z ray $
0013 %
0014 %   This file is part of MATPOWER.
0015 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0016 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0017 
0018 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0019     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0020     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0021 
0022 %% default arguments
0023 if nargin < 10
0024     mpopt = mpoption;
0025 end
0026 
0027 %% options
0028 tol     = mpopt.pf.tol;
0029 max_it  = mpopt.pf.nr.max_it;
0030 
0031 %% initialize
0032 converged = 0;
0033 i = 0;
0034 nb = length(V0);
0035 nbr = size(Yf, 1);
0036 nref = [pv;pq];             %% indices of all non-reference buses
0037 f = branch(:, F_BUS);       %% list of "from" buses
0038 t = branch(:, T_BUS);       %% list of "to" buses
0039 
0040 %%-----  evaluate Hessian  -----
0041 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V0);
0042 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V0);
0043 H = [
0044     real(dSf_dVa)   real(dSf_dVm);
0045     real(dSt_dVa)   real(dSt_dVm);
0046     real(dSbus_dVa) real(dSbus_dVm);
0047     speye(nb)       sparse(nb,nb);
0048     imag(dSf_dVa)   imag(dSf_dVm);
0049     imag(dSt_dVa)   imag(dSt_dVm);
0050     imag(dSbus_dVa) imag(dSbus_dVm);
0051     sparse(nb,nb)   speye(nb);
0052 ];
0053 
0054 %% true measurement
0055 z = [
0056     real(Sf);
0057     real(St);
0058     real(Sbus);
0059     angle(V0);
0060     imag(Sf);
0061     imag(St);
0062     imag(Sbus);
0063     abs(V0);
0064 ];
0065 
0066 %% create inverse of covariance matrix with all measurements
0067 fullscale = 30;
0068 sigma = [
0069     0.02 * abs(Sf)      + 0.0052 * fullscale * ones(nbr,1);
0070     0.02 * abs(St)      + 0.0052 * fullscale * ones(nbr,1);
0071     0.02 * abs(Sbus)    + 0.0052 * fullscale * ones(nb,1);
0072     0.2 * pi/180 * 3*ones(nb,1);
0073     0.02 * abs(Sf)      + 0.0052 * fullscale * ones(nbr,1);
0074     0.02 * abs(St)      + 0.0052 * fullscale * ones(nbr,1);
0075     0.02 * abs(Sbus)    + 0.0052 * fullscale * ones(nb,1);
0076     0.02 * abs(V0)      + 0.0052 * 1.1 * ones(nb,1);
0077 ] ./ 3;
0078 ns = length(sigma);
0079 W = sparse(1:ns, 1:ns ,  sigma .^ 2, ns, ns );
0080 WInv = sparse(1:ns, 1:ns ,  1 ./ sigma .^ 2, ns, ns );
0081 
0082 %% covariance of measurement residual
0083 %R = H * inv( H' * WInv * H ) * H';
0084 
0085 %% measurement with error
0086 err = normrnd( zeros(size(sigma)), sigma );
0087 % err = zeros(size(z));
0088 % save err err
0089 % load err
0090 % err(10) = 900 * W(10,10);     %% create a bad measurement
0091 z = z + err;
0092     
0093 %% use flat start for intial estimate
0094 V = ones(nb,1);
0095 
0096 %% compute estimated measurement
0097 Sfe = V(f) .* conj(Yf * V);
0098 Ste = V(t) .* conj(Yt * V);
0099 Sbuse = V .* conj(Ybus * V);
0100 z_est = [
0101     real(Sfe);
0102     real(Ste);
0103     real(Sbuse);
0104     angle(V);
0105     imag(Sfe);
0106     imag(Ste);
0107     imag(Sbuse);
0108     abs(V);
0109 ];
0110 
0111 %% measurement residual
0112 delz = z - z_est;
0113 normF = delz' * WInv * delz;  
0114 chusqu = err' * WInv * err;  
0115      
0116 %% check tolerance
0117 if mpopt.verbose > 1
0118     fprintf('\n it     norm( F )       step size');
0119     fprintf('\n----  --------------  --------------');
0120     fprintf('\n%3d    %10.3e      %10.3e', i, normF, 0);
0121 end
0122 if normF < tol
0123     converged = 1;
0124     if mpopt.verbose > 1
0125         fprintf('\nConverged!\n');
0126     end
0127 end
0128 
0129 %% index vector for measurements that are to be used
0130 %%%%%% NOTE: This is specific to the 30-bus system   %%%%%%
0131 %%%%%%       where bus 1 is the reference bus which  %%%%%%
0132 %%%%%%       is connected to branches 1 and 2        %%%%%%
0133 vv=[(3:nbr), ...                    %% all but 1st two Pf
0134     (nbr+1:2*nbr), ...              %% all Pt
0135     (2*nbr+2:2*nbr+nb), ...         %% all but 1st Pbus
0136     (2*nbr+nb+2:2*nbr+2*nb), ...    %% all but 1st Va
0137     (2*nbr+2*nb+3:3*nbr+2*nb), ...  %% all but 1st two Qf
0138     (3*nbr+2*nb+1:4*nbr+2*nb), ...  %% all Qt
0139     (4*nbr+2*nb+2:4*nbr+3*nb), ...  %% all but 1st Qbus
0140     (4*nbr+3*nb+2:4*nbr+4*nb)]';    %% all but 1st Vm
0141 %% index vector for state variables to be updated
0142 ww = [ nref; nb+nref ];
0143 
0144 %% bad data loop
0145 one_at_a_time = 1; max_it_bad_data = 50;
0146 % one_at_a_time = 0; max_it_bad_data = 5;
0147 ibd = 1;
0148 while (~converged && ibd <= max_it_bad_data)
0149     nm = length(vv);
0150     baddata = 0;
0151 
0152     %% find reduced Hessian, covariance matrix, measurements
0153     HH = H(vv,ww);
0154     WWInv = WInv(vv,vv);
0155     ddelz = delz(vv);
0156     VVa = angle(V(nref));
0157     VVm = abs(V(nref));
0158     
0159 %   B0 = WWInv * (err(vv) .^ 2);
0160 %   B00 = WWInv * (ddelz .^ 2);
0161 %   [maxB0,i_maxB0] = max(B0)
0162 %   [maxB00,i_maxB00] = max(B00)
0163     
0164     %%-----  do Newton iterations  -----
0165     max_it = 100;
0166     i = 0;
0167     while (~converged && i < max_it)
0168         %% update iteration counter
0169         i = i + 1;
0170         
0171         %% compute update step
0172         F = HH' * WWInv * ddelz;
0173         J = HH' * WWInv * HH;
0174         dx = (J \ F);
0175         
0176         %% update voltage
0177         VVa = VVa + dx(1:nb-1);
0178         VVm = VVm + dx(nb:2*nb-2);
0179         V(nref) = VVm .* exp(1j * VVa);
0180 
0181         %% compute estimated measurement
0182         Sfe = V(f) .* conj(Yf * V);
0183         Ste = V(t) .* conj(Yt * V);
0184         Sbuse = V .* conj(Ybus * V);
0185         z_est = [
0186             real(Sfe);
0187             real(Ste);
0188             real(Sbuse);
0189             angle(V);
0190             imag(Sfe);
0191             imag(Ste);
0192             imag(Sbuse);
0193             abs(V);
0194         ];
0195         
0196         %% measurement residual
0197         delz = z - z_est;
0198         ddelz = delz(vv);
0199         normF = ddelz' * WWInv * ddelz;  
0200 
0201         %% check for convergence
0202         step = dx' * dx;
0203         if mpopt.verbose > 1
0204             fprintf('\n%3d    %10.3e      %10.3e', i, normF, step);
0205         end
0206         if (step < tol) 
0207             converged = 1;
0208             if mpopt.verbose
0209                 fprintf('\nState estimator converged in %d iterations.\n', i);
0210             end
0211         end
0212     end
0213     if mpopt.verbose
0214         if ~converged
0215             fprintf('\nState estimator did not converge in %d iterations.\n', i);
0216         end
0217     end
0218     
0219     %%-----  Chi squared test for bad data and bad data rejection  -----
0220     B = zeros(nm,1);
0221     bad_threshold = 6.25;       %% the threshold for bad data = sigma squared
0222     RR = inv(WWInv) - 0.95 * HH * inv(HH' * WWInv * HH) * HH';
0223 %   RI = inv( inv(WWInv) - HH * inv(HH' * WWInv * HH) * HH' );
0224 %   find(eig(full(inv(WWInv) - HH * inv(HH' * WWInv * HH) * HH')) < 0)
0225 %   chi = ddelz' * RR * ddelz
0226     rr = diag(RR);
0227 
0228     B = ddelz .^ 2 ./ rr;
0229     [maxB,i_maxB] = max(B);
0230     if one_at_a_time
0231         if maxB >= bad_threshold
0232             rejected = i_maxB;
0233         else
0234             rejected = [];
0235         end
0236     else
0237         rejected = find( B >= bad_threshold );
0238     end
0239     if length(rejected)
0240         baddata = 1;
0241         converged = 0;
0242         if mpopt.verbose
0243             fprintf('\nRejecting %d measurement(s) as bad data:\n', length(rejected));
0244             fprintf('\tindex\t      B\n');
0245             fprintf('\t-----\t-------------\n');
0246             fprintf('\t%4d\t%10.2f\n', [ vv(rejected), B(rejected) ]' );
0247         end
0248         
0249         %% update measurement index vector
0250         k = find( B < bad_threshold );
0251         vv = vv(k);
0252         nm = length(vv);
0253     end
0254 
0255     if (baddata == 0) 
0256         converged = 1;
0257         if mpopt.verbose
0258             fprintf('\nNo remaining bad data, after discarding data %d time(s).\n', ibd-1);
0259             fprintf('Largest value of B = %.2f\n', maxB);
0260         end
0261     end
0262     ibd = ibd + 1;
0263 end
state_est

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
