Description of runpf

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
runpf

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
