doSE

DOSE  Do state estimation.
   created by Rui Bo on 2007/11/12

0001 function [V, converged, iterNum, z, z_est, error_sqrsum] = doSE(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V0, ref, pv, pq, measure, idx, sigma)
0002 %DOSE  Do state estimation.
0003 %   created by Rui Bo on 2007/11/12
0004 
0005 %   MATPOWER
0006 %   $Id: doSE.m,v 1.5 2010/04/26 19:45:26 ray Exp $
0007 %   by Rui Bo
0008 %   and Ray Zimmerman, PSERC Cornell
0009 %   Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC)
0010 %   Copyright (c) 2009-2010 by Rui Bo
0011 %
0012 %   This file is part of MATPOWER.
0013 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0014 %
0015 %   MATPOWER is free software: you can redistribute it and/or modify
0016 %   it under the terms of the GNU General Public License as published
0017 %   by the Free Software Foundation, either version 3 of the License,
0018 %   or (at your option) any later version.
0019 %
0020 %   MATPOWER is distributed in the hope that it will be useful,
0021 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0022 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0023 %   GNU General Public License for more details.
0024 %
0025 %   You should have received a copy of the GNU General Public License
0026 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0027 %
0028 %   Additional permission under GNU GPL version 3 section 7
0029 %
0030 %   If you modify MATPOWER, or any covered work, to interface with
0031 %   other modules (such as MATLAB code and MEX-files) available in a
0032 %   MATLAB(R) or comparable environment containing parts covered
0033 %   under other licensing terms, the licensors of MATPOWER grant
0034 %   you additional permission to convey the resulting work.
0035 
0036 %% define named indices into bus, gen, branch matrices
0037 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0038     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0039 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ...
0040     RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch;
0041 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
0042     GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;
0043 
0044 %% options
0045 tol     = 1e-5; % mpopt(2);
0046 max_it  = 100;  % mpopt(3);
0047 verbose = 0;    % mpopt(31);
0048 
0049 %% initialize
0050 j = sqrt(-1);
0051 converged = 0;
0052 i = 0;
0053 V = V0;
0054 Va = angle(V);
0055 Vm = abs(V);
0056 
0057 nb = size(Ybus, 1);
0058 f = branch(:, F_BUS);       %% list of "from" buses
0059 t = branch(:, T_BUS);       %% list of "to" buses
0060 
0061 %% get non reference buses
0062 nonref = [pv;pq];
0063 
0064 %% form measurement vector 'z'. NOTE: all are p.u. values
0065 z = [
0066     measure.PF
0067     measure.PT
0068     measure.PG
0069     measure.Va
0070     measure.QF
0071     measure.QT
0072     measure.QG
0073     measure.Vm
0074     ];
0075 
0076 %% form measurement index vectors
0077 idx_zPF = idx.idx_zPF;
0078 idx_zPT = idx.idx_zPT;
0079 idx_zPG = idx.idx_zPG;
0080 idx_zVa = idx.idx_zVa;
0081 idx_zQF = idx.idx_zQF;
0082 idx_zQT = idx.idx_zQT;
0083 idx_zQG = idx.idx_zQG;
0084 idx_zVm = idx.idx_zVm;
0085 
0086 %% get R inverse matrix
0087 sigma_vector = [
0088     sigma.sigma_PF*ones(size(idx_zPF, 1), 1)
0089     sigma.sigma_PT*ones(size(idx_zPT, 1), 1)
0090     sigma.sigma_PG*ones(size(idx_zPG, 1), 1)
0091     sigma.sigma_Va*ones(size(idx_zVa, 1), 1)
0092     sigma.sigma_QF*ones(size(idx_zQF, 1), 1)
0093     sigma.sigma_QT*ones(size(idx_zQT, 1), 1)
0094     sigma.sigma_QG*ones(size(idx_zQG, 1), 1)
0095     sigma.sigma_Vm*ones(size(idx_zVm, 1), 1)
0096     ]; % NOTE: zero-valued elements of simga are skipped
0097 sigma_square = sigma_vector.^2;
0098 R_inv = diag(1./sigma_square);
0099 
0100 %% do Newton iterations
0101 while (~converged & i < max_it)
0102     %% update iteration counter
0103     i = i + 1;
0104     
0105     %% --- compute estimated measurement ---
0106     Sfe = V(f) .* conj(Yf * V);
0107     Ste = V(t) .* conj(Yt * V);
0108     %% compute net injection at generator buses
0109     gbus = gen(:, GEN_BUS);
0110     Sgbus = V(gbus) .* conj(Ybus(gbus, :) * V);
0111     Sgen = Sgbus * baseMVA + (bus(gbus, PD) + j*bus(gbus, QD));   %% inj S + local Sd
0112     Sgen = Sgen/baseMVA;
0113     z_est = [ % NOTE: all are p.u. values
0114         real(Sfe(idx_zPF));
0115         real(Ste(idx_zPT));
0116         real(Sgen(idx_zPG));
0117         angle(V(idx_zVa));
0118         imag(Sfe(idx_zQF));
0119         imag(Ste(idx_zQT));
0120         imag(Sgen(idx_zQG));
0121         abs(V(idx_zVm));
0122     ];
0123 
0124     %% --- get H matrix ---
0125     [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);
0126     [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
0127 %     genbus_row = findBusRowByIdx(bus, gbus);
0128     genbus_row = gbus;  %% rdz, this should be fine if using internal bus numbering
0129 
0130     %% get sub-matrix of H relating to line flow
0131     dPF_dVa = real(dSf_dVa); % from end
0132     dQF_dVa = imag(dSf_dVa);   
0133     dPF_dVm = real(dSf_dVm);
0134     dQF_dVm = imag(dSf_dVm);
0135     dPT_dVa = real(dSt_dVa);% to end
0136     dQT_dVa = imag(dSt_dVa);   
0137     dPT_dVm = real(dSt_dVm);
0138     dQT_dVm = imag(dSt_dVm);   
0139     %% get sub-matrix of H relating to generator output
0140     dPG_dVa = real(dSbus_dVa(genbus_row, :));
0141     dQG_dVa = imag(dSbus_dVa(genbus_row, :));
0142     dPG_dVm = real(dSbus_dVm(genbus_row, :));
0143     dQG_dVm = imag(dSbus_dVm(genbus_row, :));
0144     %% get sub-matrix of H relating to voltage angle
0145     dVa_dVa = eye(nb);
0146     dVa_dVm = zeros(nb, nb);
0147     %% get sub-matrix of H relating to voltage magnitude
0148     dVm_dVa = zeros(nb, nb);
0149     dVm_dVm = eye(nb);
0150     H = [
0151         dPF_dVa(idx_zPF, nonref)   dPF_dVm(idx_zPF, nonref);
0152         dPT_dVa(idx_zPT, nonref)   dPT_dVm(idx_zPT, nonref);
0153         dPG_dVa(idx_zPG, nonref)   dPG_dVm(idx_zPG, nonref);
0154         dVa_dVa(idx_zVa, nonref)   dVa_dVm(idx_zVa, nonref);
0155         dQF_dVa(idx_zQF, nonref)   dQF_dVm(idx_zQF, nonref);
0156         dQT_dVa(idx_zQT, nonref)   dQT_dVm(idx_zQT, nonref);
0157         dQG_dVa(idx_zQG, nonref)   dQG_dVm(idx_zQG, nonref);
0158         dVm_dVa(idx_zVm, nonref)   dVm_dVm(idx_zVm, nonref);
0159         ];
0160     
0161     %% compute update step
0162     J = H'*R_inv*H;
0163     F = H'*R_inv*(z-z_est); % evalute F(x)
0164     if ~isobservable(H, pv, pq)
0165         error('doSE: system is not observable');
0166     end
0167     dx = (J \ F);
0168 
0169     %% check for convergence
0170     normF = norm(F, inf);
0171     if verbose > 1
0172         fprintf('\niteration [%3d]\t\tnorm of mismatch: %10.3e', i, normF);
0173     end
0174     if normF < tol
0175         converged = 1;
0176     end
0177     
0178     %% update voltage
0179     Va(nonref) = Va(nonref) + dx(1:size(nonref, 1));
0180     Vm(nonref) = Vm(nonref) + dx(size(nonref, 1)+1:2*size(nonref, 1));
0181     V = Vm .* exp(j * Va); % NOTE: angle is in radians in pf solver, but in degree in case data
0182     Vm = abs(V);            %% update Vm and Va again in case
0183     Va = angle(V);          %% we wrapped around with a negative Vm
0184 end
0185 
0186 iterNum = i;
0187 
0188 %% get weighted sum of squared errors
0189 error_sqrsum = sum((z - z_est).^2./sigma_square);