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-2019, 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-SC'
0165                 mpopt = mpoption(mpopt, 'pf.current_balance', 0, 'pf.v_cartesian', 1);
0166             case 'NR-IP'
0167                 mpopt = mpoption(mpopt, 'pf.current_balance', 1, 'pf.v_cartesian', 0);
0168             case 'NR-IC'
0169                 mpopt = mpoption(mpopt, 'pf.current_balance', 1, 'pf.v_cartesian', 1);
0170         end
0171         if mpopt.verbose > 0
0172             switch alg
0173                 case 'NR'
0174                     solver = 'Newton';
0175                 case 'NR-SC'
0176                     solver = 'Newton-SC';
0177                 case 'NR-IP'
0178                     solver = 'Newton-IP';
0179                 case 'NR-IC'
0180                     solver = 'Newton-IC';
0181                 case 'FDXB'
0182                     solver = 'fast-decoupled, XB';
0183                 case 'FDBX'
0184                     solver = 'fast-decoupled, BX';
0185                 case 'GS'
0186                     solver = 'Gauss-Seidel';
0187                 case 'PQSUM'
0188                     solver = 'Power Summation';
0189                 case 'ISUM'
0190                     solver = 'Current Summation';
0191                 case 'YSUM'
0192                     solver = 'Admittance Summation';
0193                 otherwise
0194                     solver = 'unknown';
0195             end
0196             fprintf(' -- AC Power Flow (%s)\n', solver);
0197         end
0198         switch alg
0199             case {'NR', 'NR-SC', 'NR-IP', 'NR-IC'}  %% all 4 variants supported
0200             otherwise                   %% only power balance, polar is valid
0201                 if mpopt.pf.current_balance || mpopt.pf.v_cartesian
0202                     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.');
0203                 end
0204         end
0205         if have_zip_loads(mpopt)
0206             if mpopt.pf.current_balance || mpopt.pf.v_cartesian
0207                 warnstr = 'Newton algorithm (current or cartesian versions) do';
0208             elseif strcmp(alg, 'GS')
0209                 warnstr = 'Gauss-Seidel algorithm does';
0210             else
0211                 warnstr = '';
0212             end
0213             if warnstr
0214                 warning('runpf: %s not support ZIP load model. Converting to constant power loads.', warnstr);
0215                 mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads', ...
0216                                 struct('pw', [], 'qw', []));
0217             end
0218         end
0219 
0220         %% initial state
0221         % V0    = ones(size(bus, 1), 1);            %% flat start
0222         V0  = bus(:, VM) .* exp(1j * pi/180 * bus(:, VA));
0223         vcb = ones(size(V0));           %% create mask of voltage-controlled buses
0224         vcb(pq) = 0;                    %% exclude PQ buses
0225         k = find(vcb(gbus));            %% in-service gens at v-c buses
0226         V0(gbus(k)) = gen(on(k), VG) ./ abs(V0(gbus(k))).* V0(gbus(k));
0227 
0228         if qlim
0229             ref0 = ref;                         %% save index and angle of
0230             Varef0 = bus(ref0, VA);             %%   original reference bus(es)
0231             limited = [];                       %% list of indices of gens @ Q lims
0232             fixedQg = zeros(size(gen, 1), 1);   %% Qg of gens at Q limits
0233         end
0234 
0235         %% build admittance matrices
0236         [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0237 
0238         repeat = 1;
0239         while (repeat)
0240             %% function for computing V dependent complex bus power injections
0241             %% (generation - load)
0242             Sbus = @(Vm)makeSbus(baseMVA, bus, gen, mpopt, Vm);
0243 
0244             %% run the power flow
0245             switch alg
0246                 case {'NR', 'NR-SC', 'NR-IP', 'NR-IC'}
0247                     if mpopt.pf.current_balance
0248                         if mpopt.pf.v_cartesian     %% current, cartesian
0249                             newtonpf_fcn = @newtonpf_I_cart;
0250                         else                        %% current, polar
0251                             newtonpf_fcn = @newtonpf_I_polar;
0252                         end
0253                     else
0254                         if mpopt.pf.v_cartesian     %% power, cartesian
0255                             newtonpf_fcn = @newtonpf_S_cart;
0256                         else                        %% default - power, polar
0257                             newtonpf_fcn = @newtonpf;
0258                         end
0259                     end
0260                     [V, success, iterations] = newtonpf_fcn(Ybus, Sbus, V0, ref, pv, pq, mpopt);
0261                 case {'FDXB', 'FDBX'}
0262                     [Bp, Bpp] = makeB(baseMVA, bus, branch, alg);
0263                     [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt);
0264                 case 'GS'
0265                     [V, success, iterations] = gausspf(Ybus, Sbus([]), V0, ref, pv, pq, mpopt);
0266                 case {'PQSUM', 'ISUM', 'YSUM'}
0267                     [mpc, success, iterations] = radial_pf(mpc, mpopt);
0268                 otherwise
0269                     error('runpf: ''%s'' is not a valid power flow algorithm. See ''pf.alg'' details in MPOPTION help.', alg);
0270             end
0271             its = its + iterations;
0272 
0273             %% update data matrices with solution
0274             switch alg
0275                 case {'NR', 'NR-SC', 'NR-IP', 'NR-IC', 'FDXB', 'FDBX', 'GS'}
0276                     [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt);
0277                 case {'PQSUM', 'ISUM', 'YSUM'}
0278                     bus = mpc.bus;
0279                     gen = mpc.gen;
0280                     branch = mpc.branch;
0281             end
0282 
0283             if success && qlim      %% enforce generator Q limits
0284                 %% find gens with violated Q constraints
0285                 mx = find( gen(:, GEN_STATUS) > 0 ...
0286                         & gen(:, QG) > gen(:, QMAX) + mpopt.opf.violation );
0287                 mn = find( gen(:, GEN_STATUS) > 0 ...
0288                         & gen(:, QG) < gen(:, QMIN) - mpopt.opf.violation );
0289 
0290                 if ~isempty(mx) || ~isempty(mn)  %% we have some Q limit violations
0291                     %% first check for INFEASIBILITY
0292                     infeas = union(mx', mn')';  %% transposes handle fact that
0293                         %% union of scalars is a row vector
0294                     remaining = find( gen(:, GEN_STATUS) > 0 & ...
0295                                     ( bus(gen(:, GEN_BUS), BUS_TYPE) == PV | ...
0296                                       bus(gen(:, GEN_BUS), BUS_TYPE) == REF ));
0297                     if length(infeas) == length(remaining) && all(infeas == remaining) && ...
0298                             (isempty(mx) || isempty(mn))
0299                         %% all remaining PV/REF gens are violating AND all are
0300                         %% violating same limit (all violating Qmin or all Qmax)
0301                         if mpopt.verbose
0302                             fprintf('All %d remaining gens exceed their Q limits : INFEASIBLE PROBLEM\n', length(infeas));
0303                         end
0304                         success = 0;
0305                         break;
0306                     end
0307 
0308                     %% one at a time?
0309                     if qlim == 2    %% fix largest violation, ignore the rest
0310                         [junk, k] = max([gen(mx, QG) - gen(mx, QMAX);
0311                                          gen(mn, QMIN) - gen(mn, QG)]);
0312                         if k > length(mx)
0313                             mn = mn(k-length(mx));
0314                             mx = [];
0315                         else
0316                             mx = mx(k);
0317                             mn = [];
0318                         end
0319                     end
0320 
0321                     if mpopt.verbose && ~isempty(mx)
0322                         fprintf('Gen %d at upper Q limit, converting to PQ bus\n', mx);
0323                     end
0324                     if mpopt.verbose && ~isempty(mn)
0325                         fprintf('Gen %d at lower Q limit, converting to PQ bus\n', mn);
0326                     end
0327 
0328                     %% save corresponding limit values
0329                     fixedQg(mx) = gen(mx, QMAX);
0330                     fixedQg(mn) = gen(mn, QMIN);
0331                     mx = [mx;mn];
0332 
0333                     %% convert to PQ bus
0334                     gen(mx, QG) = fixedQg(mx);      %% set Qg to binding limit
0335                     gen(mx, GEN_STATUS) = 0;        %% temporarily turn off gen,
0336                     for i = 1:length(mx)            %% (one at a time, since
0337                         bi = gen(mx(i), GEN_BUS);   %%  they may be at same bus)
0338                         bus(bi, [PD,QD]) = ...      %% adjust load accordingly,
0339                             bus(bi, [PD,QD]) - gen(mx(i), [PG,QG]);
0340                     end
0341                     if length(ref) > 1 && any(bus(gen(mx, GEN_BUS), BUS_TYPE) == REF)
0342                         error('runpf: Sorry, MATPOWER cannot enforce Q limits for slack buses in systems with multiple slacks.');
0343                     end
0344                     bus(gen(mx, GEN_BUS), BUS_TYPE) = PQ;   %% & set bus type to PQ
0345 
0346                     %% update bus index lists of each type of bus
0347                     ref_temp = ref;
0348                     [ref, pv, pq] = bustypes(bus, gen);
0349                     %% previous line can modify lists to select new REF bus
0350                     %% if there was none, so we should update bus with these
0351                     %% just to keep them consistent
0352                     if ref ~= ref_temp
0353                         bus(ref, BUS_TYPE) = REF;
0354                         bus( pv, BUS_TYPE) = PV;
0355                         if mpopt.verbose
0356                             fprintf('Bus %d is new slack bus\n', ref);
0357                         end
0358                     end
0359                     limited = [limited; mx];
0360                 else
0361                     repeat = 0; %% no more generator Q limits violated
0362                 end
0363             else
0364                 repeat = 0;     %% don't enforce generator Q limits, once is enough
0365             end
0366         end
0367         if qlim && ~isempty(limited)
0368             %% restore injections from limited gens (those at Q limits)
0369             gen(limited, QG) = fixedQg(limited);    %% restore Qg value,
0370             for i = 1:length(limited)               %% (one at a time, since
0371                 bi = gen(limited(i), GEN_BUS);      %%  they may be at same bus)
0372                 bus(bi, [PD,QD]) = ...              %% re-adjust load,
0373                     bus(bi, [PD,QD]) + gen(limited(i), [PG,QG]);
0374             end
0375             gen(limited, GEN_STATUS) = 1;               %% and turn gen back on
0376             if ref ~= ref0
0377                 %% adjust voltage angles to make original ref bus correct
0378                 bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0;
0379             end
0380         end
0381     end
0382 else
0383     t0 = tic;
0384     success = 0;
0385     its = 0;
0386     if mpopt.verbose
0387         fprintf('Power flow not valid : MATPOWER case contains no connected buses');
0388     end
0389 end
0390 mpc.et = toc(t0);
0391 mpc.success = success;
0392 mpc.iterations = its;
0393 
0394 %%-----  output results  -----
0395 %% convert back to original bus numbering & print results
0396 [mpc.bus, mpc.gen, mpc.branch] = deal(bus, gen, branch);
0397 results = int2ext(mpc);
0398 
0399 %% zero out result fields of out-of-service gens & branches
0400 if ~isempty(results.order.gen.status.off)
0401   results.gen(results.order.gen.status.off, [PG QG]) = 0;
0402 end
0403 if ~isempty(results.order.branch.status.off)
0404   results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0;
0405 end
0406 
0407 if fname
0408     [fd, msg] = fopen(fname, 'at');
0409     if fd == -1
0410         error(msg);
0411     else
0412         if mpopt.out.all == 0
0413             printpf(results, fd, mpoption(mpopt, 'out.all', -1));
0414         else
0415             printpf(results, fd, mpopt);
0416         end
0417         fclose(fd);
0418     end
0419 end
0420 printpf(results, 1, mpopt);
0421 
0422 %% save solved case
0423 if solvedcase
0424     savecase(solvedcase, results);
0425 end
0426 
0427 if nargout == 1 || nargout == 2
0428     MVAbase = results;
0429     bus = success;
0430 elseif nargout > 2
0431     [MVAbase, bus, gen, branch, et] = ...
0432         deal(results.baseMVA, results.bus, results.gen, results.branch, results.et);
0433 % else  %% don't define MVAbase, so it doesn't print anything
0434 end
0435 
0436 function TorF = have_zip_loads(mpopt)
0437 if (~isempty(mpopt.exp.sys_wide_zip_loads.pw) && ...
0438         any(mpopt.exp.sys_wide_zip_loads.pw(2:3))) || ...
0439         (~isempty(mpopt.exp.sys_wide_zip_loads.qw) && ...
0440         any(mpopt.exp.sys_wide_zip_loads.qw(2:3)))
0441     TorF = 1;
0442 else
0443     TorF = 0;
0444 end
runpf

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
