Home > matpower7.1 > extras > state_estimator > state_est.m

state_est

PURPOSE ^

STATE_EST Solves a state estimation problem.

SYNOPSIS ^

function [V, converged, i] = state_est(branch, Ybus, Yf, Yt, Sbus, V0, ref, pv, pq, mpopt)

DESCRIPTION ^

STATE_EST  Solves a state estimation problem.
   [V, CONVERGED, I] = STATE_EST(BRANCH, YBUS, YF, YT, SBUS, ...
                                   V0, REF, PV, PQ, MPOPT)
   State estimator (under construction) based on code from James S. Thorp.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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