Home > matpower4.1 > extras > state_estimator > runse.m

runse

PURPOSE ^

RUNSE Runs a state estimator.

SYNOPSIS ^

function [MVAbase, bus, gen, branch, success, et] = runse(casedata, mpopt, fname, solvedcase)

DESCRIPTION ^

RUNSE  Runs a state estimator.
   [BASEMVA, BUS, GEN, BRANCH, SUCCESS, ET] = ...
           RUNSE(CASEDATA, MPOPT, FNAME, SOLVEDCASE)

   Runs a state estimator (after a Newton power flow). Under construction with
   parts based on code from James S. Thorp.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [MVAbase, bus, gen, branch, success, et] = runse(casedata, mpopt, fname, solvedcase)
0002 %RUNSE  Runs a state estimator.
0003 %   [BASEMVA, BUS, GEN, BRANCH, SUCCESS, ET] = ...
0004 %           RUNSE(CASEDATA, MPOPT, FNAME, SOLVEDCASE)
0005 %
0006 %   Runs a state estimator (after a Newton power flow). Under construction with
0007 %   parts based on code from James S. Thorp.
0008 
0009 %   MATPOWER
0010 %   $Id: runse.m,v 1.14 2010/04/26 19:45:26 ray Exp $
0011 %   by Ray Zimmerman, PSERC Cornell
0012 %   parts based on code by James S. Thorp, June 2004
0013 %   Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC)
0014 %
0015 %   This file is part of MATPOWER.
0016 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0017 %
0018 %   MATPOWER is free software: you can redistribute it and/or modify
0019 %   it under the terms of the GNU General Public License as published
0020 %   by the Free Software Foundation, either version 3 of the License,
0021 %   or (at your option) any later version.
0022 %
0023 %   MATPOWER is distributed in the hope that it will be useful,
0024 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0025 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0026 %   GNU General Public License for more details.
0027 %
0028 %   You should have received a copy of the GNU General Public License
0029 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0030 %
0031 %   Additional permission under GNU GPL version 3 section 7
0032 %
0033 %   If you modify MATPOWER, or any covered work, to interface with
0034 %   other modules (such as MATLAB code and MEX-files) available in a
0035 %   MATLAB(R) or comparable environment containing parts covered
0036 %   under other licensing terms, the licensors of MATPOWER grant
0037 %   you additional permission to convey the resulting work.
0038 
0039 %%-----  initialize  -----
0040 %% define named indices into bus, gen, branch matrices
0041 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0042     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0043 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0044     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0045     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0046 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0047     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0048     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0049 
0050 %% default arguments
0051 if nargin < 4
0052     solvedcase = '';                %% don't save solved case
0053     if nargin < 3
0054         fname = '';                 %% don't print results to a file
0055         if nargin < 2
0056             mpopt = mpoption;       %% use default options
0057             if nargin < 1
0058                 casedata = 'case9'; %% default data file is 'case9.m'
0059             end
0060         end
0061     end
0062 end
0063 
0064 %% options
0065 dc = mpopt(10);                     %% use DC formulation?
0066 
0067 %% read data & convert to internal bus numbering
0068 [baseMVA, bus, gen, branch] = loadcase(casedata);
0069 [i2e, bus, gen, branch] = ext2int(bus, gen, branch);
0070 
0071 %% get bus index lists of each type of bus
0072 [ref, pv, pq] = bustypes(bus, gen);
0073 
0074 %% generator info
0075 on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
0076 gbus = gen(on, GEN_BUS);                %% what buses are they at?
0077 
0078 %%-----  run the power flow  -----
0079 t0 = clock;
0080 if dc                               %% DC formulation
0081     %% initial state
0082     Va0 = bus(:, VA) * (pi/180);
0083     
0084     %% build B matrices and phase shift injections
0085     [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
0086     
0087     %% compute complex bus power injections (generation - load)
0088     %% adjusted for phase shifters and real shunts
0089     Pbus = real(makeSbus(baseMVA, bus, gen)) - Pbusinj - bus(:, GS) / baseMVA;
0090     
0091     %% "run" the power flow
0092     Va = dcpf(B, Pbus, Va0, ref, pv, pq);
0093     
0094     %% update data matrices with solution
0095     branch(:, [QF, QT]) = zeros(size(branch, 1), 2);
0096     branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0097     branch(:, PT) = -branch(:, PF);
0098     bus(:, VM) = ones(size(bus, 1), 1);
0099     bus(:, VA) = Va * (180/pi);
0100     %% update Pg for swing generator (note: other gens at ref bus are accounted for in Pbus)
0101     %%      Pg = Pinj + Pload + Gs
0102     %%      newPg = oldPg + newPinj - oldPinj
0103     refgen = find(gbus == ref);             %% which is(are) the reference gen(s)?
0104     gen(on(refgen(1)), PG) = gen(on(refgen(1)), PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA;
0105     
0106     success = 1;
0107 else                                %% AC formulation
0108     %% initial state
0109     % V0    = ones(size(bus, 1), 1);            %% flat start
0110     V0  = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA));
0111     V0(gbus) = gen(on, VG) ./ abs(V0(gbus)).* V0(gbus);
0112     
0113     %% build admittance matrices
0114     [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0115     
0116     %% compute complex bus power injections (generation - load)
0117     Sbus = makeSbus(baseMVA, bus, gen);
0118     
0119     %% run the power flow
0120         alg = mpopt(1);
0121     if alg == 1
0122         [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt);
0123     elseif alg == 2 || alg == 3
0124         [Bp, Bpp] = makeB(baseMVA, bus, branch, alg);
0125         [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt);
0126     elseif alg == 4
0127         [V, success, iterations] = gausspf(Ybus, Sbus, V0, ref, pv, pq, mpopt);
0128     else
0129         error('Only Newton''s method, fast-decoupled, and Gauss-Seidel power flow algorithms currently implemented.');
0130     end
0131     
0132     %% update data matrices with solution
0133     [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq);
0134 end
0135 et = etime(clock, t0);
0136 
0137 %%--------------------  begin state estimator code  --------------------
0138 %% save some values from load flow solution
0139 Pflf=branch(:,PF);
0140 Qflf=branch(:,QF);
0141 Ptlf=branch(:,PT);
0142 Qtlf=branch(:,QT);
0143 Sbuslf = V .* conj(Ybus * V);
0144 Vlf=V;
0145 
0146 %% run state estimator
0147 [V, converged, i] = state_est(branch, Ybus, Yf, Yt, Sbuslf, Vlf, ref, pv, pq, mpopt);
0148 
0149 %% update data matrices to match estimator solution ...
0150 %% ... bus injections at PQ buses
0151 Sbus = V .* conj(Ybus * V);
0152 bus(pq, PD) = -real(Sbus(pq)) * baseMVA;
0153 bus(pq, QD) = -imag(Sbus(pq)) * baseMVA;
0154 %% ... gen outputs at PV buses
0155 on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
0156 gbus = gen(on, GEN_BUS);                %% what buses are they at?
0157 gen(on, PG) = real(Sbus(gbus)) * baseMVA + bus(gbus, PD);   %% inj P + local Pd
0158 %% ... line flows, reference bus injections, etc.
0159 [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq);
0160 
0161 %% plot differences from load flow solution
0162 Pfe=branch(:,PF);
0163 Qfe=branch(:,QF);
0164 Pte=branch(:,PT);
0165 Qte=branch(:,QT);
0166 nbr = length(Pfe);
0167 subplot(3,2,1), plot(180/pi*(angle(Vlf)-angle(V)),'.'), title('Voltage Angle (deg)');
0168 subplot(3,2,2), plot(abs(Vlf)-abs(V),'.'), title('Voltage Magnitude (p.u.)');
0169 subplot(3,2,3), plot((1:nbr),(Pfe-Pflf),'r.',(1:nbr),(Pte-Ptlf),'b.'), title('Real Flow (MW)');
0170 subplot(3,2,4), plot((1:nbr),(Qfe-Qflf),'r.',(1:nbr),(Qte-Qtlf),'b.'), title('Reactive Flow (MVAr)');
0171 subplot(3,2,5), plot(baseMVA*real(Sbuslf-Sbus), '.'), title('Real Injection (MW)');
0172 subplot(3,2,6), plot(baseMVA*imag(Sbuslf-Sbus), '.'), title('Reactive Injection (MVAr)');
0173 %%--------------------  end state estimator code  --------------------
0174 
0175 %%-----  output results  -----
0176 %% convert back to original bus numbering & print results
0177 [bus, gen, branch] = int2ext(i2e, bus, gen, branch);
0178 if fname
0179     [fd, msg] = fopen(fname, 'at');
0180     if fd == -1
0181         error(msg);
0182     else
0183         printpf(baseMVA, bus, gen, branch, [], success, et, fd, mpopt);
0184         fclose(fd);
0185     end
0186 end
0187 printpf(baseMVA, bus, gen, branch, [], success, et, 1, mpopt);
0188 
0189 %% save solved case
0190 if solvedcase
0191     savecase(solvedcase, baseMVA, bus, gen, branch);
0192 end
0193 
0194 %% this is just to prevent it from printing baseMVA
0195 %% when called with no output arguments
0196 if nargout, MVAbase = baseMVA; end
Generated on Mon 26-Jan-2015 15:00:13 by m2html © 2005