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-2016, 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 http://www.pserc.cornell.edu/matpower/ 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); 0110 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch); 0111 0112 %% get bus index lists of each type of bus 0113 [ref, pv, pq] = bustypes(bus, gen); 0114 0115 %% generator info 0116 on = find(gen(:, GEN_STATUS) > 0); %% which generators are on? 0117 gbus = gen(on, GEN_BUS); %% what buses are they at? 0118 0119 %%----- run the power flow ----- 0120 t0 = clock; 0121 its = 0; %% total iterations 0122 if mpopt.verbose > 0 0123 v = mpver('all'); 0124 fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date); 0125 end 0126 if dc %% DC formulation 0127 if mpopt.verbose > 0 0128 fprintf(' -- DC Power Flow\n'); 0129 end 0130 %% initial state 0131 Va0 = bus(:, VA) * (pi/180); 0132 0133 %% build B matrices and phase shift injections 0134 [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch); 0135 0136 %% compute complex bus power injections (generation - load) 0137 %% adjusted for phase shifters and real shunts 0138 Pbus = real(makeSbus(baseMVA, bus, gen)) - Pbusinj - bus(:, GS) / baseMVA; 0139 0140 %% "run" the power flow 0141 [Va, success] = dcpf(B, Pbus, Va0, ref, pv, pq); 0142 its = 1; 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 else %% AC formulation 0161 alg = upper(mpopt.pf.alg); 0162 if mpopt.verbose > 0 0163 switch alg 0164 case 'NR' 0165 solver = 'Newton'; 0166 case 'FDXB' 0167 solver = 'fast-decoupled, XB'; 0168 case 'FDBX' 0169 solver = 'fast-decoupled, BX'; 0170 case 'GS' 0171 solver = 'Gauss-Seidel'; 0172 otherwise 0173 solver = 'unknown'; 0174 end 0175 fprintf(' -- AC Power Flow (%s)\n', solver); 0176 end 0177 %% initial state 0178 % V0 = ones(size(bus, 1), 1); %% flat start 0179 V0 = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA)); 0180 vcb = ones(size(V0)); %% create mask of voltage-controlled buses 0181 vcb(pq) = 0; %% exclude PQ buses 0182 k = find(vcb(gbus)); %% in-service gens at v-c buses 0183 V0(gbus(k)) = gen(on(k), VG) ./ abs(V0(gbus(k))).* V0(gbus(k)); 0184 0185 if qlim 0186 ref0 = ref; %% save index and angle of 0187 Varef0 = bus(ref0, VA); %% original reference bus(es) 0188 limited = []; %% list of indices of gens @ Q lims 0189 fixedQg = zeros(size(gen, 1), 1); %% Qg of gens at Q limits 0190 end 0191 0192 %% build admittance matrices 0193 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch); 0194 0195 repeat = 1; 0196 while (repeat) 0197 %% function for computing V dependent complex bus power injections 0198 %% (generation - load) 0199 Sbus = @(Vm)makeSbus(baseMVA, bus, gen, mpopt, Vm); 0200 0201 %% run the power flow 0202 switch alg 0203 case 'NR' 0204 [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt); 0205 case {'FDXB', 'FDBX'} 0206 [Bp, Bpp] = makeB(baseMVA, bus, branch, alg); 0207 [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt); 0208 case 'GS' 0209 if (~isempty(mpopt.exp.sys_wide_zip_loads.pw) && ... 0210 any(mpopt.exp.sys_wide_zip_loads.pw(2:3))) || ... 0211 (~isempty(mpopt.exp.sys_wide_zip_loads.qw) && ... 0212 any(mpopt.exp.sys_wide_zip_loads.qw(2:3))) 0213 warning('runpf: Gauss-Seidel algorithm does not support ZIP load model. Converting to constant power loads.') 0214 mpopt = mpoption(mpopt, 'exp.sys_wide_zip_loads', ... 0215 struct('pw', [], 'qw', [])); 0216 end 0217 [V, success, iterations] = gausspf(Ybus, Sbus([]), V0, ref, pv, pq, mpopt); 0218 otherwise 0219 error('runpf: Only Newton''s method, fast-decoupled, and Gauss-Seidel power flow algorithms currently implemented.'); 0220 end 0221 its = its + iterations; 0222 0223 %% update data matrices with solution 0224 [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt); 0225 0226 if qlim %% enforce generator Q limits 0227 %% find gens with violated Q constraints 0228 mx = find( gen(:, GEN_STATUS) > 0 ... 0229 & gen(:, QG) > gen(:, QMAX) + mpopt.opf.violation ); 0230 mn = find( gen(:, GEN_STATUS) > 0 ... 0231 & gen(:, QG) < gen(:, QMIN) - mpopt.opf.violation ); 0232 0233 if ~isempty(mx) || ~isempty(mn) %% we have some Q limit violations 0234 %% first check for INFEASIBILITY 0235 infeas = union(mx', mn')'; %% transposes handle fact that 0236 %% union of scalars is a row vector 0237 remaining = find( gen(:, GEN_STATUS) > 0 & ... 0238 ( bus(gen(:, GEN_BUS), BUS_TYPE) == PV | ... 0239 bus(gen(:, GEN_BUS), BUS_TYPE) == REF )); 0240 if length(infeas) == length(remaining) && all(infeas == remaining) && ... 0241 (isempty(mx) || isempty(mn)) 0242 %% all remaining PV/REF gens are violating AND all are 0243 %% violating same limit (all violating Qmin or all Qmax) 0244 if mpopt.verbose 0245 fprintf('All %d remaining gens exceed their Q limits : INFEASIBLE PROBLEM\n', length(infeas)); 0246 end 0247 success = 0; 0248 break; 0249 end 0250 0251 %% one at a time? 0252 if qlim == 2 %% fix largest violation, ignore the rest 0253 [junk, k] = max([gen(mx, QG) - gen(mx, QMAX); 0254 gen(mn, QMIN) - gen(mn, QG)]); 0255 if k > length(mx) 0256 mn = mn(k-length(mx)); 0257 mx = []; 0258 else 0259 mx = mx(k); 0260 mn = []; 0261 end 0262 end 0263 0264 if mpopt.verbose && ~isempty(mx) 0265 fprintf('Gen %d at upper Q limit, converting to PQ bus\n', mx); 0266 end 0267 if mpopt.verbose && ~isempty(mn) 0268 fprintf('Gen %d at lower Q limit, converting to PQ bus\n', mn); 0269 end 0270 0271 %% save corresponding limit values 0272 fixedQg(mx) = gen(mx, QMAX); 0273 fixedQg(mn) = gen(mn, QMIN); 0274 mx = [mx;mn]; 0275 0276 %% convert to PQ bus 0277 gen(mx, QG) = fixedQg(mx); %% set Qg to binding limit 0278 gen(mx, GEN_STATUS) = 0; %% temporarily turn off gen, 0279 for i = 1:length(mx) %% (one at a time, since 0280 bi = gen(mx(i), GEN_BUS); %% they may be at same bus) 0281 bus(bi, [PD,QD]) = ... %% adjust load accordingly, 0282 bus(bi, [PD,QD]) - gen(mx(i), [PG,QG]); 0283 end 0284 if length(ref) > 1 && any(bus(gen(mx, GEN_BUS), BUS_TYPE) == REF) 0285 error('runpf: Sorry, MATPOWER cannot enforce Q limits for slack buses in systems with multiple slacks.'); 0286 end 0287 bus(gen(mx, GEN_BUS), BUS_TYPE) = PQ; %% & set bus type to PQ 0288 0289 %% update bus index lists of each type of bus 0290 ref_temp = ref; 0291 [ref, pv, pq] = bustypes(bus, gen); 0292 %% previous line can modify lists to select new REF bus 0293 %% if there was none, so we should update bus with these 0294 %% just to keep them consistent 0295 if ref ~= ref_temp 0296 bus(ref, BUS_TYPE) = REF; 0297 bus( pv, BUS_TYPE) = PV; 0298 if mpopt.verbose 0299 fprintf('Bus %d is new slack bus\n', ref); 0300 end 0301 end 0302 limited = [limited; mx]; 0303 else 0304 repeat = 0; %% no more generator Q limits violated 0305 end 0306 else 0307 repeat = 0; %% don't enforce generator Q limits, once is enough 0308 end 0309 end 0310 if qlim && ~isempty(limited) 0311 %% restore injections from limited gens (those at Q limits) 0312 gen(limited, QG) = fixedQg(limited); %% restore Qg value, 0313 for i = 1:length(limited) %% (one at a time, since 0314 bi = gen(limited(i), GEN_BUS); %% they may be at same bus) 0315 bus(bi, [PD,QD]) = ... %% re-adjust load, 0316 bus(bi, [PD,QD]) + gen(limited(i), [PG,QG]); 0317 end 0318 gen(limited, GEN_STATUS) = 1; %% and turn gen back on 0319 if ref ~= ref0 0320 %% adjust voltage angles to make original ref bus correct 0321 bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0; 0322 end 0323 end 0324 end 0325 mpc.et = etime(clock, t0); 0326 mpc.success = success; 0327 mpc.iterations = its; 0328 0329 %%----- output results ----- 0330 %% convert back to original bus numbering & print results 0331 [mpc.bus, mpc.gen, mpc.branch] = deal(bus, gen, branch); 0332 results = int2ext(mpc); 0333 0334 %% zero out result fields of out-of-service gens & branches 0335 if ~isempty(results.order.gen.status.off) 0336 results.gen(results.order.gen.status.off, [PG QG]) = 0; 0337 end 0338 if ~isempty(results.order.branch.status.off) 0339 results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0; 0340 end 0341 0342 if fname 0343 [fd, msg] = fopen(fname, 'at'); 0344 if fd == -1 0345 error(msg); 0346 else 0347 if mpopt.out.all == 0 0348 printpf(results, fd, mpoption(mpopt, 'out.all', -1)); 0349 else 0350 printpf(results, fd, mpopt); 0351 end 0352 fclose(fd); 0353 end 0354 end 0355 printpf(results, 1, mpopt); 0356 0357 %% save solved case 0358 if solvedcase 0359 savecase(solvedcase, results); 0360 end 0361 0362 if nargout == 1 || nargout == 2 0363 MVAbase = results; 0364 bus = success; 0365 elseif nargout > 2 0366 [MVAbase, bus, gen, branch, et] = ... 0367 deal(results.baseMVA, results.bus, results.gen, results.branch, results.et); 0368 % else %% don't define MVAbase, so it doesn't print anything 0369 end