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-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