runpf

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.

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-2020, Power Systems 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 %   This file is part of MATPOWER.
0068 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0069 %   See https://matpower.org for more info.
0070 
0071 %%-----  initialize  -----
0072 %% define named indices into bus, gen, branch matrices
0073 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0074     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0075 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0076     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0077     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0078 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0079     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0080     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0081 
0082 %% default arguments
0083 if nargin < 4
0084     solvedcase = '';                %% don't save solved case
0085     if nargin < 3
0086         fname = '';                 %% don't print results to a file
0087         if nargin < 2
0088             mpopt = mpoption;       %% use default options
0089             if nargin < 1
0090                 casedata = 'case9'; %% default data file is 'case9.m'
0091             end
0092         end
0093     end
0094 end
0095 
0096 %% options
0097 qlim = mpopt.pf.enforce_q_lims;         %% enforce Q limits on gens?
0098 dc = strcmp(upper(mpopt.model), 'DC');  %% use DC formulation?
0099 
0100 %% read data
0101 mpc = loadcase(casedata);
0102 
0103 %% add zero columns to branch for flows if needed
0104 if size(mpc.branch,2) < QT
0105   mpc.branch = [ mpc.branch zeros(size(mpc.branch, 1), QT-size(mpc.branch,2)) ];
0106 end
0107 
0108 %% convert to internal indexing
0109 mpc = ext2int(mpc, mpopt);
0110 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch);
0111 
0112 if ~isempty(mpc.bus)
0113     %% get bus index lists of each type of bus
0114     [ref, pv, pq] = bustypes(bus, gen);
0115 
0116     %% generator info
0117     on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
0118     gbus = gen(on, GEN_BUS);                %% what buses are they at?
0119 
0120     %%-----  run the power flow  -----
0121     t0 = tic;
0122     its = 0;            %% total iterations
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, success] = dcpf(B, Pbus, Va0, ref, pv, pq);
0143         its = 1;
0144 
0145         %% update data matrices with solution
0146         branch(:, [QF, QT]) = zeros(size(branch, 1), 2);
0147         branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0148         branch(:, PT) = -branch(:, PF);
0149         bus(:, VM) = ones(size(bus, 1), 1);
0150         bus(:, VA) = Va * (180/pi);
0151         %% update Pg for slack generator (1st gen at ref bus)
0152         %% (note: other gens at ref bus are accounted for in Pbus)
0153         %%      Pg = Pinj + Pload + Gs
0154         %%      newPg = oldPg + newPinj - oldPinj
0155         refgen = zeros(size(ref));
0156         for k = 1:length(ref)
0157             temp = find(gbus == ref(k));
0158             refgen(k) = on(temp(1));
0159         end
0160         gen(refgen, PG) = gen(refgen, PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA;
0161     else                                %% AC formulation
0162         alg = upper(mpopt.pf.alg);
0163         switch alg
0164             case 'NR-SP'
0165                 mpopt = mpoption(mpopt, 'pf.current_balance', 0, 'pf.v_cartesian', 0);
0166             case 'NR-SC'
0167                 mpopt = mpoption(mpopt, 'pf.current_balance', 0, 'pf.v_cartesian', 1);
0168             case 'NR-SH'
0169                 mpopt = mpoption(mpopt, 'pf.current_balance', 0, 'pf.v_cartesian', 2);
0170             case 'NR-IP'
0171                 mpopt = mpoption(mpopt, 'pf.current_balance', 1, 'pf.v_cartesian', 0);
0172             case 'NR-IC'
0173                 mpopt = mpoption(mpopt, 'pf.current_balance', 1, 'pf.v_cartesian', 1);
0174             case 'NR-IH'
0175                 mpopt = mpoption(mpopt, 'pf.current_balance', 1, 'pf.v_cartesian', 2);
0176         end
0177         if mpopt.verbose > 0
0178             switch alg
0179                 case {'NR', 'NR-SP'}
0180                     solver = 'Newton';
0181                 case 'NR-SC'
0182                     solver = 'Newton-SC';
0183                 case 'NR-SH'
0184                     solver = 'Newton-SH';
0185                 case 'NR-IP'
0186                     solver = 'Newton-IP';
0187                 case 'NR-IC'
0188                     solver = 'Newton-IC';
0189                 case 'NR-IH'
0190                     solver = 'Newton-IH';
0191                 case 'FDXB'
0192                     solver = 'fast-decoupled, XB';
0193                 case 'FDBX'
0194                     solver = 'fast-decoupled, BX';
0195                 case 'GS'
0196                     solver = 'Gauss-Seidel';
0197                 case 'PQSUM'
0198                     solver = 'Power Summation';
0199                 case 'ISUM'
0200                     solver = 'Current Summation';
0201                 case 'YSUM'
0202                     solver = 'Admittance Summation';
0203                 otherwise
0204                     solver = 'unknown';
0205             end
0206             fprintf(' -- AC Power Flow (%s)\n', solver);
0207         end
0208         switch alg
0209             case {'NR', 'NR-SP', 'NR-SC', 'NR-SH', 'NR-IP', 'NR-IC', 'NR-IH'}  %% all 6 variants supported
0210             otherwise                   %% only power balance, polar is valid
0211                 if mpopt.pf.current_balance || mpopt.pf.v_cartesian
0212                     error('runpf: power flow algorithm ''%s'' only supports power balance, polar version\nI.e. both ''pf.current_balance'' and ''pf.v_cartesian'' must be set to 0.');
0213                 end
0214         end
0215         if have_zip_loads(mpopt)
0216             if mpopt.pf.current_balance || mpopt.pf.v_cartesian
0217                 warnstr = 'Newton algorithm (current or cartesian/hybrid versions) do';
0218             elseif strcmp(alg, 'GS')
0219                 warnstr = 'Gauss-Seidel algorithm does';
0220             else
0221                 warnstr = '';
0222             end
0223             if warnstr
0224                 warning('runpf: %s not support ZIP load model. Converting to constant power loads.', warnstr);
0225                 mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads', ...
0226                                 struct('pw', [], 'qw', []));
0227             end
0228         end
0229 
0230         %% initial state
0231         % V0    = ones(size(bus, 1), 1);            %% flat start
0232         V0  = bus(:, VM) .* exp(1j * pi/180 * bus(:, VA));
0233         vcb = ones(size(V0));           %% create mask of voltage-controlled buses
0234         vcb(pq) = 0;                    %% exclude PQ buses
0235         k = find(vcb(gbus));            %% in-service gens at v-c buses
0236         V0(gbus(k)) = gen(on(k), VG) ./ abs(V0(gbus(k))).* V0(gbus(k));
0237 
0238         if qlim
0239             ref0 = ref;                         %% save index and angle of
0240             Varef0 = bus(ref0, VA);             %%   original reference bus(es)
0241             limited = [];                       %% list of indices of gens @ Q lims
0242             fixedQg = zeros(size(gen, 1), 1);   %% Qg of gens at Q limits
0243         end
0244 
0245         %% build admittance matrices
0246         [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0247 
0248         repeat = 1;
0249         while (repeat)
0250             %% function for computing V dependent complex bus power injections
0251             %% (generation - load)
0252             Sbus = @(Vm)makeSbus(baseMVA, bus, gen, mpopt, Vm);
0253 
0254             %% run the power flow
0255             switch alg
0256                 case {'NR', 'NR-SP', 'NR-SC', 'NR-SH', 'NR-IP', 'NR-IC', 'NR-IH'}
0257                     if mpopt.pf.current_balance
0258                         switch mpopt.pf.v_cartesian
0259                             case 0                  %% current, polar
0260                                 newtonpf_fcn = @newtonpf_I_polar;
0261                             case 1                  %% current, cartesian
0262                                 newtonpf_fcn = @newtonpf_I_cart;
0263                             case 2                  %% current, hybrid
0264                                 newtonpf_fcn = @newtonpf_I_hybrid;
0265                         end
0266                     else
0267                         switch mpopt.pf.v_cartesian
0268                             case 0                  %% default - power, polar
0269                                 newtonpf_fcn = @newtonpf;
0270                             case 1                  %% power, cartesian
0271                                 newtonpf_fcn = @newtonpf_S_cart;
0272                             case 2                  %% power, hybrid
0273                                 newtonpf_fcn = @newtonpf_S_hybrid;
0274                         end
0275                     end
0276                     [V, success, iterations] = newtonpf_fcn(Ybus, Sbus, V0, ref, pv, pq, mpopt);
0277                 case {'FDXB', 'FDBX'}
0278                     [Bp, Bpp] = makeB(baseMVA, bus, branch, alg);
0279                     [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt);
0280                 case 'GS'
0281                     [V, success, iterations] = gausspf(Ybus, Sbus([]), V0, ref, pv, pq, mpopt);
0282                 case {'PQSUM', 'ISUM', 'YSUM'}
0283                     [mpc, success, iterations] = radial_pf(mpc, mpopt);
0284                 otherwise
0285                     error('runpf: ''%s'' is not a valid power flow algorithm. See ''pf.alg'' details in MPOPTION help.', alg);
0286             end
0287             its = its + iterations;
0288 
0289             %% update data matrices with solution
0290             switch alg
0291                 case {'NR', 'NR-SP', 'NR-SC', 'NR-SH', 'NR-IP', 'NR-IC', 'NR-IH', 'FDXB', 'FDBX', 'GS'}
0292                     [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt);
0293                 case {'PQSUM', 'ISUM', 'YSUM'}
0294                     bus = mpc.bus;
0295                     gen = mpc.gen;
0296                     branch = mpc.branch;
0297             end
0298 
0299             if success && qlim      %% enforce generator Q limits
0300                 %% find gens with violated Q constraints
0301                 mx = find( gen(:, GEN_STATUS) > 0 ...
0302                         & gen(:, QG) > gen(:, QMAX) + mpopt.opf.violation );
0303                 mn = find( gen(:, GEN_STATUS) > 0 ...
0304                         & gen(:, QG) < gen(:, QMIN) - mpopt.opf.violation );
0305 
0306                 if ~isempty(mx) || ~isempty(mn)  %% we have some Q limit violations
0307                     %% first check for INFEASIBILITY
0308                     infeas = union(mx', mn')';  %% transposes handle fact that
0309                         %% union of scalars is a row vector
0310                     remaining = find( gen(:, GEN_STATUS) > 0 & ...
0311                                     ( bus(gen(:, GEN_BUS), BUS_TYPE) == PV | ...
0312                                       bus(gen(:, GEN_BUS), BUS_TYPE) == REF ));
0313                     if length(infeas) == length(remaining) && all(infeas == remaining) && ...
0314                             (isempty(mx) || isempty(mn))
0315                         %% all remaining PV/REF gens are violating AND all are
0316                         %% violating same limit (all violating Qmin or all Qmax)
0317                         if mpopt.verbose
0318                             fprintf('All %d remaining gens exceed their Q limits : INFEASIBLE PROBLEM\n', length(infeas));
0319                         end
0320                         success = 0;
0321                         break;
0322                     end
0323 
0324                     %% one at a time?
0325                     if qlim == 2    %% fix largest violation, ignore the rest
0326                         [junk, k] = max([gen(mx, QG) - gen(mx, QMAX);
0327                                          gen(mn, QMIN) - gen(mn, QG)]);
0328                         if k > length(mx)
0329                             mn = mn(k-length(mx));
0330                             mx = [];
0331                         else
0332                             mx = mx(k);
0333                             mn = [];
0334                         end
0335                     end
0336 
0337                     if mpopt.verbose && ~isempty(mx)
0338                         fprintf('Gen %d at upper Q limit, converting to PQ bus\n', mx);
0339                     end
0340                     if mpopt.verbose && ~isempty(mn)
0341                         fprintf('Gen %d at lower Q limit, converting to PQ bus\n', mn);
0342                     end
0343 
0344                     %% save corresponding limit values
0345                     fixedQg(mx) = gen(mx, QMAX);
0346                     fixedQg(mn) = gen(mn, QMIN);
0347                     mx = [mx;mn];
0348 
0349                     %% convert to PQ bus
0350                     gen(mx, QG) = fixedQg(mx);      %% set Qg to binding limit
0351                     if length(ref) > 1 && any(bus(gen(mx, GEN_BUS), BUS_TYPE) == REF)
0352                         error('runpf: Sorry, MATPOWER cannot enforce Q limits for slack buses in systems with multiple slacks.');
0353                     end
0354                     bus(gen(mx, GEN_BUS), BUS_TYPE) = PQ;   %% & set bus type to PQ
0355 
0356                     %% update bus index lists of each type of bus
0357                     ref_temp = ref;
0358                     [ref, pv, pq] = bustypes(bus, gen);
0359                     %% previous line can modify lists to select new REF bus
0360                     %% if there was none, so we should update bus with these
0361                     %% just to keep them consistent
0362                     if ref ~= ref_temp
0363                         bus(ref, BUS_TYPE) = REF;
0364                         bus( pv, BUS_TYPE) = PV;
0365                         if mpopt.verbose
0366                             fprintf('Bus %d is new slack bus\n', ref);
0367                         end
0368                     end
0369                     limited = [limited; mx];
0370                 else
0371                     repeat = 0; %% no more generator Q limits violated
0372                 end
0373             else
0374                 repeat = 0;     %% don't enforce generator Q limits, once is enough
0375             end
0376         end
0377         if qlim && ~isempty(limited)
0378             if ref ~= ref0
0379                 %% adjust voltage angles to make original ref bus correct
0380                 bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0;
0381             end
0382         end
0383     end
0384 else
0385     t0 = tic;
0386     success = 0;
0387     its = 0;
0388     if mpopt.verbose
0389         fprintf('Power flow not valid : MATPOWER case contains no connected buses\n');
0390     end
0391 end
0392 mpc.et = toc(t0);
0393 mpc.success = success;
0394 mpc.iterations = its;
0395 
0396 %%-----  output results  -----
0397 %% convert back to original bus numbering & print results
0398 [mpc.bus, mpc.gen, mpc.branch] = deal(bus, gen, branch);
0399 results = int2ext(mpc);
0400 
0401 %% zero out result fields of out-of-service gens & branches
0402 if ~isempty(results.order.gen.status.off)
0403   results.gen(results.order.gen.status.off, [PG QG]) = 0;
0404 end
0405 if ~isempty(results.order.branch.status.off)
0406   results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0;
0407 end
0408 
0409 if fname
0410     [fd, msg] = fopen(fname, 'at');
0411     if fd == -1
0412         error(msg);
0413     else
0414         if mpopt.out.all == 0
0415             printpf(results, fd, mpoption(mpopt, 'out.all', -1));
0416         else
0417             printpf(results, fd, mpopt);
0418         end
0419         fclose(fd);
0420     end
0421 end
0422 printpf(results, 1, mpopt);
0423 
0424 %% save solved case
0425 if solvedcase
0426     savecase(solvedcase, results);
0427 end
0428 
0429 if nargout == 1 || nargout == 2
0430     MVAbase = results;
0431     bus = success;
0432 elseif nargout > 2
0433     [MVAbase, bus, gen, branch, et] = ...
0434         deal(results.baseMVA, results.bus, results.gen, results.branch, results.et);
0435 % else  %% don't define MVAbase, so it doesn't print anything
0436 end
0437 
0438 function TorF = have_zip_loads(mpopt)
0439 if (~isempty(mpopt.exp.sys_wide_zip_loads.pw) && ...
0440         any(mpopt.exp.sys_wide_zip_loads.pw(2:3))) || ...
0441         (~isempty(mpopt.exp.sys_wide_zip_loads.qw) && ...
0442         any(mpopt.exp.sys_wide_zip_loads.qw(2:3)))
0443     TorF = 1;
0444 else
0445     TorF = 0;
0446 end