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 % $Id: runpf.m 2414 2014-11-11 18:30:25Z ray $ 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 % Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC) 0067 % 0068 % This file is part of MATPOWER. 0069 % See http://www.pserc.cornell.edu/matpower/ for more info. 0070 % 0071 % MATPOWER is free software: you can redistribute it and/or modify 0072 % it under the terms of the GNU General Public License as published 0073 % by the Free Software Foundation, either version 3 of the License, 0074 % or (at your option) any later version. 0075 % 0076 % MATPOWER is distributed in the hope that it will be useful, 0077 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0078 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0079 % GNU General Public License for more details. 0080 % 0081 % You should have received a copy of the GNU General Public License 0082 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0083 % 0084 % Additional permission under GNU GPL version 3 section 7 0085 % 0086 % If you modify MATPOWER, or any covered work, to interface with 0087 % other modules (such as MATLAB code and MEX-files) available in a 0088 % MATLAB(R) or comparable environment containing parts covered 0089 % under other licensing terms, the licensors of MATPOWER grant 0090 % you additional permission to convey the resulting work. 0091 0092 %%----- initialize ----- 0093 %% define named indices into bus, gen, branch matrices 0094 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0095 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0096 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0097 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0098 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0099 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0100 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0101 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0102 0103 %% default arguments 0104 if nargin < 4 0105 solvedcase = ''; %% don't save solved case 0106 if nargin < 3 0107 fname = ''; %% don't print results to a file 0108 if nargin < 2 0109 mpopt = mpoption; %% use default options 0110 if nargin < 1 0111 casedata = 'case9'; %% default data file is 'case9.m' 0112 end 0113 end 0114 end 0115 end 0116 0117 %% options 0118 qlim = mpopt.pf.enforce_q_lims; %% enforce Q limits on gens? 0119 dc = strcmp(upper(mpopt.model), 'DC'); %% use DC formulation? 0120 0121 %% read data 0122 mpc = loadcase(casedata); 0123 0124 %% add zero columns to branch for flows if needed 0125 if size(mpc.branch,2) < QT 0126 mpc.branch = [ mpc.branch zeros(size(mpc.branch, 1), QT-size(mpc.branch,2)) ]; 0127 end 0128 0129 %% convert to internal indexing 0130 mpc = ext2int(mpc); 0131 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch); 0132 0133 %% get bus index lists of each type of bus 0134 [ref, pv, pq] = bustypes(bus, gen); 0135 0136 %% generator info 0137 on = find(gen(:, GEN_STATUS) > 0); %% which generators are on? 0138 gbus = gen(on, GEN_BUS); %% what buses are they at? 0139 0140 %%----- run the power flow ----- 0141 t0 = clock; 0142 if mpopt.verbose > 0 0143 v = mpver('all'); 0144 fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date); 0145 end 0146 if dc %% DC formulation 0147 if mpopt.verbose > 0 0148 fprintf(' -- DC Power Flow\n'); 0149 end 0150 %% initial state 0151 Va0 = bus(:, VA) * (pi/180); 0152 0153 %% build B matrices and phase shift injections 0154 [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch); 0155 0156 %% compute complex bus power injections (generation - load) 0157 %% adjusted for phase shifters and real shunts 0158 Pbus = real(makeSbus(baseMVA, bus, gen)) - Pbusinj - bus(:, GS) / baseMVA; 0159 0160 %% "run" the power flow 0161 Va = dcpf(B, Pbus, Va0, ref, pv, pq); 0162 0163 %% update data matrices with solution 0164 branch(:, [QF, QT]) = zeros(size(branch, 1), 2); 0165 branch(:, PF) = (Bf * Va + Pfinj) * baseMVA; 0166 branch(:, PT) = -branch(:, PF); 0167 bus(:, VM) = ones(size(bus, 1), 1); 0168 bus(:, VA) = Va * (180/pi); 0169 %% update Pg for slack generator (1st gen at ref bus) 0170 %% (note: other gens at ref bus are accounted for in Pbus) 0171 %% Pg = Pinj + Pload + Gs 0172 %% newPg = oldPg + newPinj - oldPinj 0173 refgen = zeros(size(ref)); 0174 for k = 1:length(ref) 0175 temp = find(gbus == ref(k)); 0176 refgen(k) = on(temp(1)); 0177 end 0178 gen(refgen, PG) = gen(refgen, PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA; 0179 0180 success = 1; 0181 else %% AC formulation 0182 alg = upper(mpopt.pf.alg); 0183 if mpopt.verbose > 0 0184 switch alg 0185 case 'NR' 0186 solver = 'Newton'; 0187 case 'FDXB' 0188 solver = 'fast-decoupled, XB'; 0189 case 'FDBX' 0190 solver = 'fast-decoupled, BX'; 0191 case 'GS' 0192 solver = 'Gauss-Seidel'; 0193 otherwise 0194 solver = 'unknown'; 0195 end 0196 fprintf(' -- AC Power Flow (%s)\n', solver); 0197 end 0198 %% initial state 0199 % V0 = ones(size(bus, 1), 1); %% flat start 0200 V0 = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA)); 0201 vcb = ones(size(V0)); %% create mask of voltage-controlled buses 0202 vcb(pq) = 0; %% exclude PQ buses 0203 k = find(vcb(gbus)); %% in-service gens at v-c buses 0204 V0(gbus(k)) = gen(on(k), VG) ./ abs(V0(gbus(k))).* V0(gbus(k)); 0205 0206 if qlim 0207 ref0 = ref; %% save index and angle of 0208 Varef0 = bus(ref0, VA); %% original reference bus(es) 0209 limited = []; %% list of indices of gens @ Q lims 0210 fixedQg = zeros(size(gen, 1), 1); %% Qg of gens at Q limits 0211 end 0212 0213 %% build admittance matrices 0214 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch); 0215 0216 repeat = 1; 0217 while (repeat) 0218 %% compute complex bus power injections (generation - load) 0219 Sbus = makeSbus(baseMVA, bus, gen); 0220 0221 %% run the power flow 0222 switch alg 0223 case 'NR' 0224 [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt); 0225 case {'FDXB', 'FDBX'} 0226 [Bp, Bpp] = makeB(baseMVA, bus, branch, alg); 0227 [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt); 0228 case 'GS' 0229 [V, success, iterations] = gausspf(Ybus, Sbus, V0, ref, pv, pq, mpopt); 0230 otherwise 0231 error('Only Newton''s method, fast-decoupled, and Gauss-Seidel power flow algorithms currently implemented.'); 0232 end 0233 0234 %% update data matrices with solution 0235 [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq); 0236 0237 if qlim %% enforce generator Q limits 0238 %% find gens with violated Q constraints 0239 mx = find( gen(:, GEN_STATUS) > 0 ... 0240 & gen(:, QG) > gen(:, QMAX) + mpopt.opf.violation ); 0241 mn = find( gen(:, GEN_STATUS) > 0 ... 0242 & gen(:, QG) < gen(:, QMIN) - mpopt.opf.violation ); 0243 0244 if ~isempty(mx) || ~isempty(mn) %% we have some Q limit violations 0245 %% first check for INFEASIBILITY (all remaining gens violating) 0246 infeas = union(mx', mn')'; %% transposes handle fact that 0247 %% union of scalars is a row vector 0248 remaining = find( gen(:, GEN_STATUS) > 0 & ... 0249 ( bus(gen(:, GEN_BUS), BUS_TYPE) == PV | ... 0250 bus(gen(:, GEN_BUS), BUS_TYPE) == REF )); 0251 if length(infeas) == length(remaining) && all(infeas == remaining) 0252 if mpopt.verbose 0253 fprintf('All %d remaining gens exceed their Q limits : INFEASIBLE PROBLEM\n', length(infeas)); 0254 end 0255 success = 0; 0256 break; 0257 end 0258 0259 %% one at a time? 0260 if qlim == 2 %% fix largest violation, ignore the rest 0261 [junk, k] = max([gen(mx, QG) - gen(mx, QMAX); 0262 gen(mn, QMIN) - gen(mn, QG)]); 0263 if k > length(mx) 0264 mn = mn(k-length(mx)); 0265 mx = []; 0266 else 0267 mx = mx(k); 0268 mn = []; 0269 end 0270 end 0271 0272 if mpopt.verbose && ~isempty(mx) 0273 fprintf('Gen %d at upper Q limit, converting to PQ bus\n', mx); 0274 end 0275 if mpopt.verbose && ~isempty(mn) 0276 fprintf('Gen %d at lower Q limit, converting to PQ bus\n', mn); 0277 end 0278 0279 %% save corresponding limit values 0280 fixedQg(mx) = gen(mx, QMAX); 0281 fixedQg(mn) = gen(mn, QMIN); 0282 mx = [mx;mn]; 0283 0284 %% convert to PQ bus 0285 gen(mx, QG) = fixedQg(mx); %% set Qg to binding limit 0286 gen(mx, GEN_STATUS) = 0; %% temporarily turn off gen, 0287 for i = 1:length(mx) %% (one at a time, since 0288 bi = gen(mx(i), GEN_BUS); %% they may be at same bus) 0289 bus(bi, [PD,QD]) = ... %% adjust load accordingly, 0290 bus(bi, [PD,QD]) - gen(mx(i), [PG,QG]); 0291 end 0292 if length(ref) > 1 && any(bus(gen(mx, GEN_BUS), BUS_TYPE) == REF) 0293 error('Sorry, MATPOWER cannot enforce Q limits for slack buses in systems with multiple slacks.'); 0294 end 0295 bus(gen(mx, GEN_BUS), BUS_TYPE) = PQ; %% & set bus type to PQ 0296 0297 %% update bus index lists of each type of bus 0298 ref_temp = ref; 0299 [ref, pv, pq] = bustypes(bus, gen); 0300 %% previous line can modify lists to select new REF bus 0301 %% if there was none, so we should update bus with these 0302 %% just to keep them consistent 0303 if ref ~= ref_temp 0304 bus(ref, BUS_TYPE) = REF; 0305 bus( pv, BUS_TYPE) = PV; 0306 if mpopt.verbose 0307 fprintf('Bus %d is new slack bus\n', ref); 0308 end 0309 end 0310 limited = [limited; mx]; 0311 else 0312 repeat = 0; %% no more generator Q limits violated 0313 end 0314 else 0315 repeat = 0; %% don't enforce generator Q limits, once is enough 0316 end 0317 end 0318 if qlim && ~isempty(limited) 0319 %% restore injections from limited gens (those at Q limits) 0320 gen(limited, QG) = fixedQg(limited); %% restore Qg value, 0321 for i = 1:length(limited) %% (one at a time, since 0322 bi = gen(limited(i), GEN_BUS); %% they may be at same bus) 0323 bus(bi, [PD,QD]) = ... %% re-adjust load, 0324 bus(bi, [PD,QD]) + gen(limited(i), [PG,QG]); 0325 end 0326 gen(limited, GEN_STATUS) = 1; %% and turn gen back on 0327 if ref ~= ref0 0328 %% adjust voltage angles to make original ref bus correct 0329 bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0; 0330 end 0331 end 0332 end 0333 mpc.et = etime(clock, t0); 0334 mpc.success = success; 0335 0336 %%----- output results ----- 0337 %% convert back to original bus numbering & print results 0338 [mpc.bus, mpc.gen, mpc.branch] = deal(bus, gen, branch); 0339 results = int2ext(mpc); 0340 0341 %% zero out result fields of out-of-service gens & branches 0342 if ~isempty(results.order.gen.status.off) 0343 results.gen(results.order.gen.status.off, [PG QG]) = 0; 0344 end 0345 if ~isempty(results.order.branch.status.off) 0346 results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0; 0347 end 0348 0349 if fname 0350 [fd, msg] = fopen(fname, 'at'); 0351 if fd == -1 0352 error(msg); 0353 else 0354 if mpopt.out.all == 0 0355 printpf(results, fd, mpoption(mpopt, 'out.all', -1)); 0356 else 0357 printpf(results, fd, mpopt); 0358 end 0359 fclose(fd); 0360 end 0361 end 0362 printpf(results, 1, mpopt); 0363 0364 %% save solved case 0365 if solvedcase 0366 savecase(solvedcase, results); 0367 end 0368 0369 if nargout == 1 || nargout == 2 0370 MVAbase = results; 0371 bus = success; 0372 elseif nargout > 2 0373 [MVAbase, bus, gen, branch, et] = ... 0374 deal(results.baseMVA, results.bus, results.gen, results.branch, results.et); 0375 % else %% don't define MVAbase, so it doesn't print anything 0376 end