Home > matpower5.0 > runpf.m

runpf

PURPOSE ^

RUNPF Runs a power flow.

SYNOPSIS ^

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

DESCRIPTION ^

RUNPF  Runs a power flow.
   [RESULTS, SUCCESS] = RUNPF(CASEDATA, MPOPT, FNAME, SOLVEDCASE)

   Runs a power flow (full AC Newton's method by default), optionally
   returning a RESULTS struct and SUCCESS flag.

   Inputs (all are optional):
       CASEDATA : either a MATPOWER case struct or a string containing
           the name of the file with the case data (default is 'case9')
           (see also CASEFORMAT and LOADCASE)
       MPOPT : MATPOWER options struct to override default options
           can be used to specify the solution algorithm, output options
           termination tolerances, and more (see also MPOPTION).
       FNAME : name of a file to which the pretty-printed output will
           be appended
       SOLVEDCASE : name of file to which the solved case will be saved
           in MATPOWER case format (M-file will be assumed unless the
           specified name ends with '.mat')

   Outputs (all are optional):
       RESULTS : results struct, with the following fields:
           (all fields from the input MATPOWER case, i.e. bus, branch,
               gen, etc., but with solved voltages, power flows, etc.)
           order - info used in external <-> internal data conversion
           et - elapsed time in seconds
           success - success flag, 1 = succeeded, 0 = failed
       SUCCESS : the success flag can additionally be returned as
           a second output argument

   Calling syntax options:
       results = runpf;
       results = runpf(casedata);
       results = runpf(casedata, mpopt);
       results = runpf(casedata, mpopt, fname);
       results = runpf(casedata, mpopt, fname, solvedcase);
       [results, success] = runpf(...);

       Alternatively, for compatibility with previous versions of MATPOWER,
       some of the results can be returned as individual output arguments:

       [baseMVA, bus, gen, branch, success, et] = runpf(...);

   If the pf.enforce_q_lims option is set to true (default is false) then, if
   any generator reactive power limit is violated after running the AC power
   flow, the corresponding bus is converted to a PQ bus, with Qg at the
   limit, and the case is re-run. The voltage magnitude at the bus will
   deviate from the specified value in order to satisfy the reactive power
   limit. If the reference bus is converted to PQ, the first remaining PV
   bus will be used as the slack bus for the next iteration. This may
   result in the real power output at this generator being slightly off
   from the specified values.

   Examples:
       results = runpf('case30');
       results = runpf('case30', mpoption('pf.enforce_q_lims', 1));

   See also RUNDCPF.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [MVAbase, bus, gen, branch, success, et] = ...
0002                 runpf(casedata, mpopt, fname, solvedcase)
0003 %RUNPF  Runs a power flow.
0004 %   [RESULTS, SUCCESS] = RUNPF(CASEDATA, MPOPT, FNAME, SOLVEDCASE)
0005 %
0006 %   Runs a power flow (full AC Newton's method by default), optionally
0007 %   returning a RESULTS struct and SUCCESS flag.
0008 %
0009 %   Inputs (all are optional):
0010 %       CASEDATA : either a MATPOWER case struct or a string containing
0011 %           the name of the file with the case data (default is 'case9')
0012 %           (see also CASEFORMAT and LOADCASE)
0013 %       MPOPT : MATPOWER options struct to override default options
0014 %           can be used to specify the solution algorithm, output options
0015 %           termination tolerances, and more (see also MPOPTION).
0016 %       FNAME : name of a file to which the pretty-printed output will
0017 %           be appended
0018 %       SOLVEDCASE : name of file to which the solved case will be saved
0019 %           in MATPOWER case format (M-file will be assumed unless the
0020 %           specified name ends with '.mat')
0021 %
0022 %   Outputs (all are optional):
0023 %       RESULTS : results struct, with the following fields:
0024 %           (all fields from the input MATPOWER case, i.e. bus, branch,
0025 %               gen, etc., but with solved voltages, power flows, etc.)
0026 %           order - info used in external <-> internal data conversion
0027 %           et - elapsed time in seconds
0028 %           success - success flag, 1 = succeeded, 0 = failed
0029 %       SUCCESS : the success flag can additionally be returned as
0030 %           a second output argument
0031 %
0032 %   Calling syntax options:
0033 %       results = runpf;
0034 %       results = runpf(casedata);
0035 %       results = runpf(casedata, mpopt);
0036 %       results = runpf(casedata, mpopt, fname);
0037 %       results = runpf(casedata, mpopt, fname, solvedcase);
0038 %       [results, success] = runpf(...);
0039 %
0040 %       Alternatively, for compatibility with previous versions of MATPOWER,
0041 %       some of the results can be returned as individual output arguments:
0042 %
0043 %       [baseMVA, bus, gen, branch, success, et] = runpf(...);
0044 %
0045 %   If the pf.enforce_q_lims option is set to true (default is false) then, if
0046 %   any generator reactive power limit is violated after running the AC power
0047 %   flow, the corresponding bus is converted to a PQ bus, with Qg at the
0048 %   limit, and the case is re-run. The voltage magnitude at the bus will
0049 %   deviate from the specified value in order to satisfy the reactive power
0050 %   limit. If the reference bus is converted to PQ, the first remaining PV
0051 %   bus will be used as the slack bus for the next iteration. This may
0052 %   result in the real power output at this generator being slightly off
0053 %   from the specified values.
0054 %
0055 %   Examples:
0056 %       results = runpf('case30');
0057 %       results = runpf('case30', mpoption('pf.enforce_q_lims', 1));
0058 %
0059 %   See also RUNDCPF.
0060 
0061 %   MATPOWER
0062 %   $Id: runpf.m 2414 2014-11-11 18:30:25Z ray $
0063 %   by Ray Zimmerman, PSERC Cornell
0064 %   Enforcing of generator Q limits inspired by contributions
0065 %   from Mu Lin, Lincoln University, New Zealand (1/14/05).
0066 %   Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC)
0067 %
0068 %   This file is part of MATPOWER.
0069 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0070 %
0071 %   MATPOWER is free software: you can redistribute it and/or modify
0072 %   it under the terms of the GNU General Public License as published
0073 %   by the Free Software Foundation, either version 3 of the License,
0074 %   or (at your option) any later version.
0075 %
0076 %   MATPOWER is distributed in the hope that it will be useful,
0077 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0078 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0079 %   GNU General Public License for more details.
0080 %
0081 %   You should have received a copy of the GNU General Public License
0082 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0083 %
0084 %   Additional permission under GNU GPL version 3 section 7
0085 %
0086 %   If you modify MATPOWER, or any covered work, to interface with
0087 %   other modules (such as MATLAB code and MEX-files) available in a
0088 %   MATLAB(R) or comparable environment containing parts covered
0089 %   under other licensing terms, the licensors of MATPOWER grant
0090 %   you additional permission to convey the resulting work.
0091 
0092 %%-----  initialize  -----
0093 %% define named indices into bus, gen, branch matrices
0094 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0095     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0096 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0097     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0098     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0099 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0100     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0101     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0102 
0103 %% default arguments
0104 if nargin < 4
0105     solvedcase = '';                %% don't save solved case
0106     if nargin < 3
0107         fname = '';                 %% don't print results to a file
0108         if nargin < 2
0109             mpopt = mpoption;       %% use default options
0110             if nargin < 1
0111                 casedata = 'case9'; %% default data file is 'case9.m'
0112             end
0113         end
0114     end
0115 end
0116 
0117 %% options
0118 qlim = mpopt.pf.enforce_q_lims;         %% enforce Q limits on gens?
0119 dc = strcmp(upper(mpopt.model), 'DC');  %% use DC formulation?
0120 
0121 %% read data
0122 mpc = loadcase(casedata);
0123 
0124 %% add zero columns to branch for flows if needed
0125 if size(mpc.branch,2) < QT
0126   mpc.branch = [ mpc.branch zeros(size(mpc.branch, 1), QT-size(mpc.branch,2)) ];
0127 end
0128 
0129 %% convert to internal indexing
0130 mpc = ext2int(mpc);
0131 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0132 
0133 %% get bus index lists of each type of bus
0134 [ref, pv, pq] = bustypes(bus, gen);
0135 
0136 %% generator info
0137 on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
0138 gbus = gen(on, GEN_BUS);                %% what buses are they at?
0139 
0140 %%-----  run the power flow  -----
0141 t0 = clock;
0142 if mpopt.verbose > 0
0143     v = mpver('all');
0144     fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date);
0145 end
0146 if dc                               %% DC formulation
0147     if mpopt.verbose > 0
0148       fprintf(' -- DC Power Flow\n');
0149     end
0150     %% initial state
0151     Va0 = bus(:, VA) * (pi/180);
0152     
0153     %% build B matrices and phase shift injections
0154     [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
0155     
0156     %% compute complex bus power injections (generation - load)
0157     %% adjusted for phase shifters and real shunts
0158     Pbus = real(makeSbus(baseMVA, bus, gen)) - Pbusinj - bus(:, GS) / baseMVA;
0159     
0160     %% "run" the power flow
0161     Va = dcpf(B, Pbus, Va0, ref, pv, pq);
0162     
0163     %% update data matrices with solution
0164     branch(:, [QF, QT]) = zeros(size(branch, 1), 2);
0165     branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0166     branch(:, PT) = -branch(:, PF);
0167     bus(:, VM) = ones(size(bus, 1), 1);
0168     bus(:, VA) = Va * (180/pi);
0169     %% update Pg for slack generator (1st gen at ref bus)
0170     %% (note: other gens at ref bus are accounted for in Pbus)
0171     %%      Pg = Pinj + Pload + Gs
0172     %%      newPg = oldPg + newPinj - oldPinj
0173     refgen = zeros(size(ref));
0174     for k = 1:length(ref)
0175         temp = find(gbus == ref(k));
0176         refgen(k) = on(temp(1));
0177     end
0178     gen(refgen, PG) = gen(refgen, PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA;
0179     
0180     success = 1;
0181 else                                %% AC formulation
0182     alg = upper(mpopt.pf.alg);
0183     if mpopt.verbose > 0
0184         switch alg
0185             case 'NR'
0186                 solver = 'Newton';
0187             case 'FDXB'
0188                 solver = 'fast-decoupled, XB';
0189             case 'FDBX'
0190                 solver = 'fast-decoupled, BX';
0191             case 'GS'
0192                 solver = 'Gauss-Seidel';
0193             otherwise
0194                 solver = 'unknown';
0195         end
0196         fprintf(' -- AC Power Flow (%s)\n', solver);
0197     end
0198     %% initial state
0199     % V0    = ones(size(bus, 1), 1);            %% flat start
0200     V0  = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA));
0201     vcb = ones(size(V0));           %% create mask of voltage-controlled buses
0202     vcb(pq) = 0;                    %% exclude PQ buses
0203     k = find(vcb(gbus));            %% in-service gens at v-c buses
0204     V0(gbus(k)) = gen(on(k), VG) ./ abs(V0(gbus(k))).* V0(gbus(k));
0205     
0206     if qlim
0207         ref0 = ref;                         %% save index and angle of
0208         Varef0 = bus(ref0, VA);             %%   original reference bus(es)
0209         limited = [];                       %% list of indices of gens @ Q lims
0210         fixedQg = zeros(size(gen, 1), 1);   %% Qg of gens at Q limits
0211     end
0212 
0213     %% build admittance matrices
0214     [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0215     
0216     repeat = 1;
0217     while (repeat)
0218         %% compute complex bus power injections (generation - load)
0219         Sbus = makeSbus(baseMVA, bus, gen);
0220         
0221         %% run the power flow
0222         switch alg
0223             case 'NR'
0224                 [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt);
0225             case {'FDXB', 'FDBX'}
0226                 [Bp, Bpp] = makeB(baseMVA, bus, branch, alg);
0227                 [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt);
0228             case 'GS'
0229                 [V, success, iterations] = gausspf(Ybus, Sbus, V0, ref, pv, pq, mpopt);
0230             otherwise
0231                 error('Only Newton''s method, fast-decoupled, and Gauss-Seidel power flow algorithms currently implemented.');
0232         end
0233         
0234         %% update data matrices with solution
0235         [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq);
0236         
0237         if qlim             %% enforce generator Q limits
0238             %% find gens with violated Q constraints
0239             mx = find( gen(:, GEN_STATUS) > 0 ...
0240                     & gen(:, QG) > gen(:, QMAX) + mpopt.opf.violation );
0241             mn = find( gen(:, GEN_STATUS) > 0 ...
0242                     & gen(:, QG) < gen(:, QMIN) - mpopt.opf.violation );
0243             
0244             if ~isempty(mx) || ~isempty(mn)  %% we have some Q limit violations
0245                 %% first check for INFEASIBILITY (all remaining gens violating)
0246                 infeas = union(mx', mn')';  %% transposes handle fact that
0247                     %% union of scalars is a row vector
0248                 remaining = find( gen(:, GEN_STATUS) > 0 & ...
0249                                 ( bus(gen(:, GEN_BUS), BUS_TYPE) == PV | ...
0250                                   bus(gen(:, GEN_BUS), BUS_TYPE) == REF ));
0251                 if length(infeas) == length(remaining) && all(infeas == remaining)
0252                     if mpopt.verbose
0253                         fprintf('All %d remaining gens exceed their Q limits : INFEASIBLE PROBLEM\n', length(infeas));
0254                     end
0255                     success = 0;
0256                     break;
0257                 end
0258 
0259                 %% one at a time?
0260                 if qlim == 2    %% fix largest violation, ignore the rest
0261                     [junk, k] = max([gen(mx, QG) - gen(mx, QMAX);
0262                                      gen(mn, QMIN) - gen(mn, QG)]);
0263                     if k > length(mx)
0264                         mn = mn(k-length(mx));
0265                         mx = [];
0266                     else
0267                         mx = mx(k);
0268                         mn = [];
0269                     end
0270                 end
0271 
0272                 if mpopt.verbose && ~isempty(mx)
0273                     fprintf('Gen %d at upper Q limit, converting to PQ bus\n', mx);
0274                 end
0275                 if mpopt.verbose && ~isempty(mn)
0276                     fprintf('Gen %d at lower Q limit, converting to PQ bus\n', mn);
0277                 end
0278                 
0279                 %% save corresponding limit values
0280                 fixedQg(mx) = gen(mx, QMAX);
0281                 fixedQg(mn) = gen(mn, QMIN);
0282                 mx = [mx;mn];
0283                 
0284                 %% convert to PQ bus
0285                 gen(mx, QG) = fixedQg(mx);      %% set Qg to binding limit
0286                 gen(mx, GEN_STATUS) = 0;        %% temporarily turn off gen,
0287                 for i = 1:length(mx)            %% (one at a time, since
0288                     bi = gen(mx(i), GEN_BUS);   %%  they may be at same bus)
0289                     bus(bi, [PD,QD]) = ...      %% adjust load accordingly,
0290                         bus(bi, [PD,QD]) - gen(mx(i), [PG,QG]);
0291                 end
0292                 if length(ref) > 1 && any(bus(gen(mx, GEN_BUS), BUS_TYPE) == REF)
0293                     error('Sorry, MATPOWER cannot enforce Q limits for slack buses in systems with multiple slacks.');
0294                 end
0295                 bus(gen(mx, GEN_BUS), BUS_TYPE) = PQ;   %% & set bus type to PQ
0296                 
0297                 %% update bus index lists of each type of bus
0298                 ref_temp = ref;
0299                 [ref, pv, pq] = bustypes(bus, gen);
0300                 %% previous line can modify lists to select new REF bus
0301                 %% if there was none, so we should update bus with these
0302                 %% just to keep them consistent
0303                 if ref ~= ref_temp
0304                     bus(ref, BUS_TYPE) = REF;
0305                     bus( pv, BUS_TYPE) = PV;
0306                     if mpopt.verbose
0307                         fprintf('Bus %d is new slack bus\n', ref);
0308                     end
0309                 end
0310                 limited = [limited; mx];
0311             else
0312                 repeat = 0; %% no more generator Q limits violated
0313             end
0314         else
0315             repeat = 0;     %% don't enforce generator Q limits, once is enough
0316         end
0317     end
0318     if qlim && ~isempty(limited)
0319         %% restore injections from limited gens (those at Q limits)
0320         gen(limited, QG) = fixedQg(limited);    %% restore Qg value,
0321         for i = 1:length(limited)               %% (one at a time, since
0322             bi = gen(limited(i), GEN_BUS);      %%  they may be at same bus)
0323             bus(bi, [PD,QD]) = ...              %% re-adjust load,
0324                 bus(bi, [PD,QD]) + gen(limited(i), [PG,QG]);
0325         end
0326         gen(limited, GEN_STATUS) = 1;               %% and turn gen back on
0327         if ref ~= ref0
0328             %% adjust voltage angles to make original ref bus correct
0329             bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0;
0330         end
0331     end
0332 end
0333 mpc.et = etime(clock, t0);
0334 mpc.success = success;
0335 
0336 %%-----  output results  -----
0337 %% convert back to original bus numbering & print results
0338 [mpc.bus, mpc.gen, mpc.branch] = deal(bus, gen, branch);
0339 results = int2ext(mpc);
0340 
0341 %% zero out result fields of out-of-service gens & branches
0342 if ~isempty(results.order.gen.status.off)
0343   results.gen(results.order.gen.status.off, [PG QG]) = 0;
0344 end
0345 if ~isempty(results.order.branch.status.off)
0346   results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0;
0347 end
0348 
0349 if fname
0350     [fd, msg] = fopen(fname, 'at');
0351     if fd == -1
0352         error(msg);
0353     else
0354         if mpopt.out.all == 0
0355             printpf(results, fd, mpoption(mpopt, 'out.all', -1));
0356         else
0357             printpf(results, fd, mpopt);
0358         end
0359         fclose(fd);
0360     end
0361 end
0362 printpf(results, 1, mpopt);
0363 
0364 %% save solved case
0365 if solvedcase
0366     savecase(solvedcase, results);
0367 end
0368 
0369 if nargout == 1 || nargout == 2
0370     MVAbase = results;
0371     bus = success;
0372 elseif nargout > 2
0373     [MVAbase, bus, gen, branch, et] = ...
0374         deal(results.baseMVA, results.bus, results.gen, results.branch, results.et);
0375 % else  %% don't define MVAbase, so it doesn't print anything
0376 end

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005