Home > matpower5.1 > 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 %   Copyright (c) 1996-2015 by Power System Engineering Research Center (PSERC)
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 %
0067 %   $Id: runpf.m 2644 2015-03-11 19:34:22Z ray $
0068 %
0069 %   This file is part of MATPOWER.
0070 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0071 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0072 
0073 %%-----  initialize  -----
0074 %% define named indices into bus, gen, branch matrices
0075 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0076     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0077 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0078     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0079     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0080 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0081     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0082     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0083 
0084 %% default arguments
0085 if nargin < 4
0086     solvedcase = '';                %% don't save solved case
0087     if nargin < 3
0088         fname = '';                 %% don't print results to a file
0089         if nargin < 2
0090             mpopt = mpoption;       %% use default options
0091             if nargin < 1
0092                 casedata = 'case9'; %% default data file is 'case9.m'
0093             end
0094         end
0095     end
0096 end
0097 
0098 %% options
0099 qlim = mpopt.pf.enforce_q_lims;         %% enforce Q limits on gens?
0100 dc = strcmp(upper(mpopt.model), 'DC');  %% use DC formulation?
0101 
0102 %% read data
0103 mpc = loadcase(casedata);
0104 
0105 %% add zero columns to branch for flows if needed
0106 if size(mpc.branch,2) < QT
0107   mpc.branch = [ mpc.branch zeros(size(mpc.branch, 1), QT-size(mpc.branch,2)) ];
0108 end
0109 
0110 %% convert to internal indexing
0111 mpc = ext2int(mpc);
0112 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0113 
0114 %% get bus index lists of each type of bus
0115 [ref, pv, pq] = bustypes(bus, gen);
0116 
0117 %% generator info
0118 on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
0119 gbus = gen(on, GEN_BUS);                %% what buses are they at?
0120 
0121 %%-----  run the power flow  -----
0122 t0 = clock;
0123 if mpopt.verbose > 0
0124     v = mpver('all');
0125     fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date);
0126 end
0127 if dc                               %% DC formulation
0128     if mpopt.verbose > 0
0129       fprintf(' -- DC Power Flow\n');
0130     end
0131     %% initial state
0132     Va0 = bus(:, VA) * (pi/180);
0133     
0134     %% build B matrices and phase shift injections
0135     [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
0136     
0137     %% compute complex bus power injections (generation - load)
0138     %% adjusted for phase shifters and real shunts
0139     Pbus = real(makeSbus(baseMVA, bus, gen)) - Pbusinj - bus(:, GS) / baseMVA;
0140     
0141     %% "run" the power flow
0142     Va = dcpf(B, Pbus, Va0, ref, pv, pq);
0143     
0144     %% update data matrices with solution
0145     branch(:, [QF, QT]) = zeros(size(branch, 1), 2);
0146     branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0147     branch(:, PT) = -branch(:, PF);
0148     bus(:, VM) = ones(size(bus, 1), 1);
0149     bus(:, VA) = Va * (180/pi);
0150     %% update Pg for slack generator (1st gen at ref bus)
0151     %% (note: other gens at ref bus are accounted for in Pbus)
0152     %%      Pg = Pinj + Pload + Gs
0153     %%      newPg = oldPg + newPinj - oldPinj
0154     refgen = zeros(size(ref));
0155     for k = 1:length(ref)
0156         temp = find(gbus == ref(k));
0157         refgen(k) = on(temp(1));
0158     end
0159     gen(refgen, PG) = gen(refgen, PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA;
0160     
0161     success = 1;
0162 else                                %% AC formulation
0163     alg = upper(mpopt.pf.alg);
0164     if mpopt.verbose > 0
0165         switch alg
0166             case 'NR'
0167                 solver = 'Newton';
0168             case 'FDXB'
0169                 solver = 'fast-decoupled, XB';
0170             case 'FDBX'
0171                 solver = 'fast-decoupled, BX';
0172             case 'GS'
0173                 solver = 'Gauss-Seidel';
0174             otherwise
0175                 solver = 'unknown';
0176         end
0177         fprintf(' -- AC Power Flow (%s)\n', solver);
0178     end
0179     %% initial state
0180     % V0    = ones(size(bus, 1), 1);            %% flat start
0181     V0  = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA));
0182     vcb = ones(size(V0));           %% create mask of voltage-controlled buses
0183     vcb(pq) = 0;                    %% exclude PQ buses
0184     k = find(vcb(gbus));            %% in-service gens at v-c buses
0185     V0(gbus(k)) = gen(on(k), VG) ./ abs(V0(gbus(k))).* V0(gbus(k));
0186     
0187     if qlim
0188         ref0 = ref;                         %% save index and angle of
0189         Varef0 = bus(ref0, VA);             %%   original reference bus(es)
0190         limited = [];                       %% list of indices of gens @ Q lims
0191         fixedQg = zeros(size(gen, 1), 1);   %% Qg of gens at Q limits
0192     end
0193 
0194     %% build admittance matrices
0195     [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0196     
0197     repeat = 1;
0198     while (repeat)
0199         %% compute complex bus power injections (generation - load)
0200         Sbus = makeSbus(baseMVA, bus, gen);
0201         
0202         %% run the power flow
0203         switch alg
0204             case 'NR'
0205                 [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt);
0206             case {'FDXB', 'FDBX'}
0207                 [Bp, Bpp] = makeB(baseMVA, bus, branch, alg);
0208                 [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt);
0209             case 'GS'
0210                 [V, success, iterations] = gausspf(Ybus, Sbus, V0, ref, pv, pq, mpopt);
0211             otherwise
0212                 error('Only Newton''s method, fast-decoupled, and Gauss-Seidel power flow algorithms currently implemented.');
0213         end
0214         
0215         %% update data matrices with solution
0216         [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq);
0217         
0218         if qlim             %% enforce generator Q limits
0219             %% find gens with violated Q constraints
0220             mx = find( gen(:, GEN_STATUS) > 0 ...
0221                     & gen(:, QG) > gen(:, QMAX) + mpopt.opf.violation );
0222             mn = find( gen(:, GEN_STATUS) > 0 ...
0223                     & gen(:, QG) < gen(:, QMIN) - mpopt.opf.violation );
0224             
0225             if ~isempty(mx) || ~isempty(mn)  %% we have some Q limit violations
0226                 %% first check for INFEASIBILITY (all remaining gens violating)
0227                 infeas = union(mx', mn')';  %% transposes handle fact that
0228                     %% union of scalars is a row vector
0229                 remaining = find( gen(:, GEN_STATUS) > 0 & ...
0230                                 ( bus(gen(:, GEN_BUS), BUS_TYPE) == PV | ...
0231                                   bus(gen(:, GEN_BUS), BUS_TYPE) == REF ));
0232                 if length(infeas) == length(remaining) && all(infeas == remaining)
0233                     if mpopt.verbose
0234                         fprintf('All %d remaining gens exceed their Q limits : INFEASIBLE PROBLEM\n', length(infeas));
0235                     end
0236                     success = 0;
0237                     break;
0238                 end
0239 
0240                 %% one at a time?
0241                 if qlim == 2    %% fix largest violation, ignore the rest
0242                     [junk, k] = max([gen(mx, QG) - gen(mx, QMAX);
0243                                      gen(mn, QMIN) - gen(mn, QG)]);
0244                     if k > length(mx)
0245                         mn = mn(k-length(mx));
0246                         mx = [];
0247                     else
0248                         mx = mx(k);
0249                         mn = [];
0250                     end
0251                 end
0252 
0253                 if mpopt.verbose && ~isempty(mx)
0254                     fprintf('Gen %d at upper Q limit, converting to PQ bus\n', mx);
0255                 end
0256                 if mpopt.verbose && ~isempty(mn)
0257                     fprintf('Gen %d at lower Q limit, converting to PQ bus\n', mn);
0258                 end
0259                 
0260                 %% save corresponding limit values
0261                 fixedQg(mx) = gen(mx, QMAX);
0262                 fixedQg(mn) = gen(mn, QMIN);
0263                 mx = [mx;mn];
0264                 
0265                 %% convert to PQ bus
0266                 gen(mx, QG) = fixedQg(mx);      %% set Qg to binding limit
0267                 gen(mx, GEN_STATUS) = 0;        %% temporarily turn off gen,
0268                 for i = 1:length(mx)            %% (one at a time, since
0269                     bi = gen(mx(i), GEN_BUS);   %%  they may be at same bus)
0270                     bus(bi, [PD,QD]) = ...      %% adjust load accordingly,
0271                         bus(bi, [PD,QD]) - gen(mx(i), [PG,QG]);
0272                 end
0273                 if length(ref) > 1 && any(bus(gen(mx, GEN_BUS), BUS_TYPE) == REF)
0274                     error('Sorry, MATPOWER cannot enforce Q limits for slack buses in systems with multiple slacks.');
0275                 end
0276                 bus(gen(mx, GEN_BUS), BUS_TYPE) = PQ;   %% & set bus type to PQ
0277                 
0278                 %% update bus index lists of each type of bus
0279                 ref_temp = ref;
0280                 [ref, pv, pq] = bustypes(bus, gen);
0281                 %% previous line can modify lists to select new REF bus
0282                 %% if there was none, so we should update bus with these
0283                 %% just to keep them consistent
0284                 if ref ~= ref_temp
0285                     bus(ref, BUS_TYPE) = REF;
0286                     bus( pv, BUS_TYPE) = PV;
0287                     if mpopt.verbose
0288                         fprintf('Bus %d is new slack bus\n', ref);
0289                     end
0290                 end
0291                 limited = [limited; mx];
0292             else
0293                 repeat = 0; %% no more generator Q limits violated
0294             end
0295         else
0296             repeat = 0;     %% don't enforce generator Q limits, once is enough
0297         end
0298     end
0299     if qlim && ~isempty(limited)
0300         %% restore injections from limited gens (those at Q limits)
0301         gen(limited, QG) = fixedQg(limited);    %% restore Qg value,
0302         for i = 1:length(limited)               %% (one at a time, since
0303             bi = gen(limited(i), GEN_BUS);      %%  they may be at same bus)
0304             bus(bi, [PD,QD]) = ...              %% re-adjust load,
0305                 bus(bi, [PD,QD]) + gen(limited(i), [PG,QG]);
0306         end
0307         gen(limited, GEN_STATUS) = 1;               %% and turn gen back on
0308         if ref ~= ref0
0309             %% adjust voltage angles to make original ref bus correct
0310             bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0;
0311         end
0312     end
0313 end
0314 mpc.et = etime(clock, t0);
0315 mpc.success = success;
0316 
0317 %%-----  output results  -----
0318 %% convert back to original bus numbering & print results
0319 [mpc.bus, mpc.gen, mpc.branch] = deal(bus, gen, branch);
0320 results = int2ext(mpc);
0321 
0322 %% zero out result fields of out-of-service gens & branches
0323 if ~isempty(results.order.gen.status.off)
0324   results.gen(results.order.gen.status.off, [PG QG]) = 0;
0325 end
0326 if ~isempty(results.order.branch.status.off)
0327   results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0;
0328 end
0329 
0330 if fname
0331     [fd, msg] = fopen(fname, 'at');
0332     if fd == -1
0333         error(msg);
0334     else
0335         if mpopt.out.all == 0
0336             printpf(results, fd, mpoption(mpopt, 'out.all', -1));
0337         else
0338             printpf(results, fd, mpopt);
0339         end
0340         fclose(fd);
0341     end
0342 end
0343 printpf(results, 1, mpopt);
0344 
0345 %% save solved case
0346 if solvedcase
0347     savecase(solvedcase, results);
0348 end
0349 
0350 if nargout == 1 || nargout == 2
0351     MVAbase = results;
0352     bus = success;
0353 elseif nargout > 2
0354     [MVAbase, bus, gen, branch, et] = ...
0355         deal(results.baseMVA, results.bus, results.gen, results.branch, results.et);
0356 % else  %% don't define MVAbase, so it doesn't print anything
0357 end

Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005