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,v 1.11 2010/04/26 19:45:26 ray Exp $
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(2);
0048 max_it  = mpopt(3);
0049 verbose = mpopt(31);
0050 
0051 %% initialize
0052 converged = 0;
0053 i = 0;
0054 nb = length(V0);
0055 nbr = size(Yf, 1);
0056 nref = [pv;pq];             %% indices of all non-reference buses
0057 f = branch(:, F_BUS);       %% list of "from" buses
0058 t = branch(:, T_BUS);       %% list of "to" buses
0059 
0060 %%-----  evaluate Hessian  -----
0061 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V0);
0062 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V0);
0063 H = [
0064     real(dSf_dVa)   real(dSf_dVm);
0065     real(dSt_dVa)   real(dSt_dVm);
0066     real(dSbus_dVa) real(dSbus_dVm);
0067     speye(nb)       sparse(nb,nb);
0068     imag(dSf_dVa)   imag(dSf_dVm);
0069     imag(dSt_dVa)   imag(dSt_dVm);
0070     imag(dSbus_dVa) imag(dSbus_dVm);
0071     sparse(nb,nb)   speye(nb);
0072 ];
0073 
0074 %% true measurement
0075 z = [
0076     real(Sf);
0077     real(St);
0078     real(Sbus);
0079     angle(V0);
0080     imag(Sf);
0081     imag(St);
0082     imag(Sbus);
0083     abs(V0);
0084 ];
0085 
0086 %% create inverse of covariance matrix with all measurements
0087 fullscale = 30;
0088 sigma = [
0089     0.02 * abs(Sf)      + 0.0052 * fullscale * ones(nbr,1);
0090     0.02 * abs(St)      + 0.0052 * fullscale * ones(nbr,1);
0091     0.02 * abs(Sbus)    + 0.0052 * fullscale * ones(nb,1);
0092     0.2 * pi/180 * 3*ones(nb,1);
0093     0.02 * abs(Sf)      + 0.0052 * fullscale * ones(nbr,1);
0094     0.02 * abs(St)      + 0.0052 * fullscale * ones(nbr,1);
0095     0.02 * abs(Sbus)    + 0.0052 * fullscale * ones(nb,1);
0096     0.02 * abs(V0)      + 0.0052 * 1.1 * ones(nb,1);
0097 ] ./ 3;
0098 ns = length(sigma);
0099 W = sparse(1:ns, 1:ns ,  sigma .^ 2, ns, ns );
0100 WInv = sparse(1:ns, 1:ns ,  1 ./ sigma .^ 2, ns, ns );
0101 
0102 %% covariance of measurement residual
0103 %R = H * inv( H' * WInv * H ) * H';
0104 
0105 %% measurement with error
0106 err = normrnd( zeros(size(sigma)), sigma );
0107 % err = zeros(size(z));
0108 % save err err
0109 % load err
0110 % err(10) = 900 * W(10,10);     %% create a bad measurement
0111 z = z + err;
0112     
0113 %% use flat start for intial estimate
0114 V = ones(nb,1);
0115 
0116 %% compute estimated measurement
0117 Sfe = V(f) .* conj(Yf * V);
0118 Ste = V(t) .* conj(Yt * V);
0119 Sbuse = V .* conj(Ybus * V);
0120 z_est = [
0121     real(Sfe);
0122     real(Ste);
0123     real(Sbuse);
0124     angle(V);
0125     imag(Sfe);
0126     imag(Ste);
0127     imag(Sbuse);
0128     abs(V);
0129 ];
0130 
0131 %% measurement residual
0132 delz = z - z_est;
0133 normF = delz' * WInv * delz;  
0134 chusqu = err' * WInv * err;  
0135      
0136 %% check tolerance
0137 if verbose > 1
0138     fprintf('\n it     norm( F )       step size');
0139     fprintf('\n----  --------------  --------------');
0140     fprintf('\n%3d    %10.3e      %10.3e', i, normF, 0);
0141 end
0142 if normF < tol
0143     converged = 1;
0144     if verbose > 1
0145         fprintf('\nConverged!\n');
0146     end
0147 end
0148 
0149 %% index vector for measurements that are to be used
0150 %%%%%% NOTE: This is specific to the 30-bus system   %%%%%%
0151 %%%%%%       where bus 1 is the reference bus which  %%%%%%
0152 %%%%%%       is connected to branches 1 and 2        %%%%%%
0153 vv=[(3:nbr), ...                    %% all but 1st two Pf
0154     (nbr+1:2*nbr), ...              %% all Pt
0155     (2*nbr+2:2*nbr+nb), ...         %% all but 1st Pbus
0156     (2*nbr+nb+2:2*nbr+2*nb), ...    %% all but 1st Va
0157     (2*nbr+2*nb+3:3*nbr+2*nb), ...  %% all but 1st two Qf
0158     (3*nbr+2*nb+1:4*nbr+2*nb), ...  %% all Qt
0159     (4*nbr+2*nb+2:4*nbr+3*nb), ...  %% all but 1st Qbus
0160     (4*nbr+3*nb+2:4*nbr+4*nb)]';    %% all but 1st Vm
0161 %% index vector for state variables to be updated
0162 ww = [ nref; nb+nref ];
0163 
0164 %% bad data loop
0165 one_at_a_time = 1; max_it_bad_data = 50;
0166 % one_at_a_time = 0; max_it_bad_data = 5;
0167 ibd = 1;
0168 while (~converged && ibd <= max_it_bad_data)
0169     nm = length(vv);
0170     baddata = 0;
0171 
0172     %% find reduced Hessian, covariance matrix, measurements
0173     HH = H(vv,ww);
0174     WWInv = WInv(vv,vv);
0175     ddelz = delz(vv);
0176     VVa = angle(V(nref));
0177     VVm = abs(V(nref));
0178     
0179 %   B0 = WWInv * (err(vv) .^ 2);
0180 %   B00 = WWInv * (ddelz .^ 2);
0181 %   [maxB0,i_maxB0] = max(B0)
0182 %   [maxB00,i_maxB00] = max(B00)
0183     
0184     %%-----  do Newton iterations  -----
0185     max_it = 100;
0186     i = 0;
0187     while (~converged && i < max_it)
0188         %% update iteration counter
0189         i = i + 1;
0190         
0191         %% compute update step
0192         F = HH' * WWInv * ddelz;
0193         J = HH' * WWInv * HH;
0194         dx = (J \ F);
0195         
0196         %% update voltage
0197         VVa = VVa + dx(1:nb-1);
0198         VVm = VVm + dx(nb:2*nb-2);
0199         V(nref) = VVm .* exp(1j * VVa);
0200 
0201         %% compute estimated measurement
0202         Sfe = V(f) .* conj(Yf * V);
0203         Ste = V(t) .* conj(Yt * V);
0204         Sbuse = V .* conj(Ybus * V);
0205         z_est = [
0206             real(Sfe);
0207             real(Ste);
0208             real(Sbuse);
0209             angle(V);
0210             imag(Sfe);
0211             imag(Ste);
0212             imag(Sbuse);
0213             abs(V);
0214         ];
0215         
0216         %% measurement residual
0217         delz = z - z_est;
0218         ddelz = delz(vv);
0219         normF = ddelz' * WWInv * ddelz;  
0220 
0221         %% check for convergence
0222         step = dx' * dx;
0223         if verbose > 1
0224             fprintf('\n%3d    %10.3e      %10.3e', i, normF, step);
0225         end
0226         if (step < tol) 
0227             converged = 1;
0228             if verbose
0229                 fprintf('\nState estimator converged in %d iterations.\n', i);
0230             end
0231         end
0232     end
0233     if verbose
0234         if ~converged
0235             fprintf('\nState estimator did not converge in %d iterations.\n', i);
0236         end
0237     end
0238     
0239     %%-----  Chi squared test for bad data and bad data rejection  -----
0240     B = zeros(nm,1);
0241     bad_threshold = 6.25;       %% the threshold for bad data = sigma squared
0242     RR = inv(WWInv) - 0.95 * HH * inv(HH' * WWInv * HH) * HH';
0243 %   RI = inv( inv(WWInv) - HH * inv(HH' * WWInv * HH) * HH' );
0244 %   find(eig(full(inv(WWInv) - HH * inv(HH' * WWInv * HH) * HH')) < 0)
0245 %   chi = ddelz' * RR * ddelz
0246     rr = diag(RR);
0247 
0248     B = ddelz .^ 2 ./ rr;
0249     [maxB,i_maxB] = max(B);
0250     if one_at_a_time
0251         if maxB >= bad_threshold
0252             rejected = i_maxB;
0253         else
0254             rejected = [];
0255         end
0256     else
0257         rejected = find( B >= bad_threshold );
0258     end
0259     if length(rejected)
0260         baddata = 1;
0261         converged = 0;
0262         if verbose
0263             fprintf('\nRejecting %d measurement(s) as bad data:\n', length(rejected));
0264             fprintf('\tindex\t      B\n');
0265             fprintf('\t-----\t-------------\n');
0266             fprintf('\t%4d\t%10.2f\n', [ vv(rejected), B(rejected) ]' );
0267         end
0268         
0269         %% update measurement index vector
0270         k = find( B < bad_threshold );
0271         vv = vv(k);
0272         nm = length(vv);
0273     end
0274 
0275     if (baddata == 0) 
0276         converged = 1;
0277         if verbose
0278             fprintf('\nNo remaining bad data, after discarding data %d time(s).\n', ibd-1);
0279             fprintf('Largest value of B = %.2f\n', maxB);
0280         end
0281     end
0282     ibd = ibd + 1;
0283 end
state_est

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
