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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
