Home > matpower5.1 > extras > se > doSE.m

doSE

PURPOSE ^

DOSE Do state estimation.

SYNOPSIS ^

function [V, converged, iterNum, z, z_est, error_sqrsum] = doSE(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V0, ref, pv, pq, measure, idx, sigma)

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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 %   Copyright (c) 1996-2015 by Power System Engineering Research Center (PSERC)
0007 %   by Rui Bo
0008 %   and Ray Zimmerman, PSERC Cornell
0009 %
0010 %   $Id: doSE.m 2644 2015-03-11 19:34:22Z ray $
0011 %
0012 %   This file is part of MATPOWER.
0013 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0014 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0015 
0016 %% define named indices into bus, gen, branch matrices
0017 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0018     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0019 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ...
0020     RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch;
0021 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
0022     GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;
0023 
0024 %% options
0025 tol     = 1e-5; % mpopt.pf.tol;
0026 max_it  = 100;  % mpopt.pf.nr.max_it;
0027 verbose = 0;
0028 
0029 %% initialize
0030 j = sqrt(-1);
0031 converged = 0;
0032 i = 0;
0033 V = V0;
0034 Va = angle(V);
0035 Vm = abs(V);
0036 
0037 nb = size(Ybus, 1);
0038 f = branch(:, F_BUS);       %% list of "from" buses
0039 t = branch(:, T_BUS);       %% list of "to" buses
0040 
0041 %% get non reference buses
0042 nonref = [pv;pq];
0043 
0044 %% form measurement vector 'z'. NOTE: all are p.u. values
0045 z = [
0046     measure.PF
0047     measure.PT
0048     measure.PG
0049     measure.Va
0050     measure.QF
0051     measure.QT
0052     measure.QG
0053     measure.Vm
0054     ];
0055 
0056 %% form measurement index vectors
0057 idx_zPF = idx.idx_zPF;
0058 idx_zPT = idx.idx_zPT;
0059 idx_zPG = idx.idx_zPG;
0060 idx_zVa = idx.idx_zVa;
0061 idx_zQF = idx.idx_zQF;
0062 idx_zQT = idx.idx_zQT;
0063 idx_zQG = idx.idx_zQG;
0064 idx_zVm = idx.idx_zVm;
0065 
0066 %% get R inverse matrix
0067 sigma_vector = [
0068     sigma.sigma_PF*ones(size(idx_zPF, 1), 1)
0069     sigma.sigma_PT*ones(size(idx_zPT, 1), 1)
0070     sigma.sigma_PG*ones(size(idx_zPG, 1), 1)
0071     sigma.sigma_Va*ones(size(idx_zVa, 1), 1)
0072     sigma.sigma_QF*ones(size(idx_zQF, 1), 1)
0073     sigma.sigma_QT*ones(size(idx_zQT, 1), 1)
0074     sigma.sigma_QG*ones(size(idx_zQG, 1), 1)
0075     sigma.sigma_Vm*ones(size(idx_zVm, 1), 1)
0076     ]; % NOTE: zero-valued elements of simga are skipped
0077 sigma_square = sigma_vector.^2;
0078 R_inv = diag(1./sigma_square);
0079 
0080 %% do Newton iterations
0081 while (~converged & i < max_it)
0082     %% update iteration counter
0083     i = i + 1;
0084     
0085     %% --- compute estimated measurement ---
0086     Sfe = V(f) .* conj(Yf * V);
0087     Ste = V(t) .* conj(Yt * V);
0088     %% compute net injection at generator buses
0089     gbus = gen(:, GEN_BUS);
0090     Sgbus = V(gbus) .* conj(Ybus(gbus, :) * V);
0091     Sgen = Sgbus * baseMVA + (bus(gbus, PD) + j*bus(gbus, QD));   %% inj S + local Sd
0092     Sgen = Sgen/baseMVA;
0093     z_est = [ % NOTE: all are p.u. values
0094         real(Sfe(idx_zPF));
0095         real(Ste(idx_zPT));
0096         real(Sgen(idx_zPG));
0097         angle(V(idx_zVa));
0098         imag(Sfe(idx_zQF));
0099         imag(Ste(idx_zQT));
0100         imag(Sgen(idx_zQG));
0101         abs(V(idx_zVm));
0102     ];
0103 
0104     %% --- get H matrix ---
0105     [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);
0106     [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
0107 %     genbus_row = findBusRowByIdx(bus, gbus);
0108     genbus_row = gbus;  %% rdz, this should be fine if using internal bus numbering
0109 
0110     %% get sub-matrix of H relating to line flow
0111     dPF_dVa = real(dSf_dVa); % from end
0112     dQF_dVa = imag(dSf_dVa);   
0113     dPF_dVm = real(dSf_dVm);
0114     dQF_dVm = imag(dSf_dVm);
0115     dPT_dVa = real(dSt_dVa);% to end
0116     dQT_dVa = imag(dSt_dVa);   
0117     dPT_dVm = real(dSt_dVm);
0118     dQT_dVm = imag(dSt_dVm);   
0119     %% get sub-matrix of H relating to generator output
0120     dPG_dVa = real(dSbus_dVa(genbus_row, :));
0121     dQG_dVa = imag(dSbus_dVa(genbus_row, :));
0122     dPG_dVm = real(dSbus_dVm(genbus_row, :));
0123     dQG_dVm = imag(dSbus_dVm(genbus_row, :));
0124     %% get sub-matrix of H relating to voltage angle
0125     dVa_dVa = eye(nb);
0126     dVa_dVm = zeros(nb, nb);
0127     %% get sub-matrix of H relating to voltage magnitude
0128     dVm_dVa = zeros(nb, nb);
0129     dVm_dVm = eye(nb);
0130     H = [
0131         dPF_dVa(idx_zPF, nonref)   dPF_dVm(idx_zPF, nonref);
0132         dPT_dVa(idx_zPT, nonref)   dPT_dVm(idx_zPT, nonref);
0133         dPG_dVa(idx_zPG, nonref)   dPG_dVm(idx_zPG, nonref);
0134         dVa_dVa(idx_zVa, nonref)   dVa_dVm(idx_zVa, nonref);
0135         dQF_dVa(idx_zQF, nonref)   dQF_dVm(idx_zQF, nonref);
0136         dQT_dVa(idx_zQT, nonref)   dQT_dVm(idx_zQT, nonref);
0137         dQG_dVa(idx_zQG, nonref)   dQG_dVm(idx_zQG, nonref);
0138         dVm_dVa(idx_zVm, nonref)   dVm_dVm(idx_zVm, nonref);
0139         ];
0140     
0141     %% compute update step
0142     J = H'*R_inv*H;
0143     F = H'*R_inv*(z-z_est); % evalute F(x)
0144     if ~isobservable(H, pv, pq)
0145         error('doSE: system is not observable');
0146     end
0147     dx = (J \ F);
0148 
0149     %% check for convergence
0150     normF = norm(F, inf);
0151     if verbose > 1
0152         fprintf('\niteration [%3d]\t\tnorm of mismatch: %10.3e', i, normF);
0153     end
0154     if normF < tol
0155         converged = 1;
0156     end
0157     
0158     %% update voltage
0159     Va(nonref) = Va(nonref) + dx(1:size(nonref, 1));
0160     Vm(nonref) = Vm(nonref) + dx(size(nonref, 1)+1:2*size(nonref, 1));
0161     V = Vm .* exp(j * Va); % NOTE: angle is in radians in pf solver, but in degree in case data
0162     Vm = abs(V);            %% update Vm and Va again in case
0163     Va = angle(V);          %% we wrapped around with a negative Vm
0164 end
0165 
0166 iterNum = i;
0167 
0168 %% get weighted sum of squared errors
0169 error_sqrsum = sum((z - z_est).^2./sigma_square);
Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005