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 vector 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 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('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 vector 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 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('ENFORCE_Q_LIMS', 1)); 0058 % 0059 % See also RUNDCPF. 0060 0061 % MATPOWER 0062 % $Id: runpf.m,v 1.27 2011/12/14 17:05:18 cvs Exp $ 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 verbose = mpopt(31); 0119 qlim = mpopt(6); %% enforce Q limits on gens? 0120 dc = mpopt(10); %% use DC formulation? 0121 0122 %% read data 0123 mpc = loadcase(casedata); 0124 0125 %% add zero columns to branch for flows if needed 0126 if size(mpc.branch,2) < QT 0127 mpc.branch = [ mpc.branch zeros(size(mpc.branch, 1), QT-size(mpc.branch,2)) ]; 0128 end 0129 0130 %% convert to internal indexing 0131 mpc = ext2int(mpc); 0132 [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch); 0133 0134 %% get bus index lists of each type of bus 0135 [ref, pv, pq] = bustypes(bus, gen); 0136 0137 %% generator info 0138 on = find(gen(:, GEN_STATUS) > 0); %% which generators are on? 0139 gbus = gen(on, GEN_BUS); %% what buses are they at? 0140 0141 %%----- run the power flow ----- 0142 t0 = clock; 0143 if verbose > 0 0144 v = mpver('all'); 0145 fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date); 0146 end 0147 if dc %% DC formulation 0148 if verbose > 0 0149 fprintf(' -- DC Power Flow\n'); 0150 end 0151 %% initial state 0152 Va0 = bus(:, VA) * (pi/180); 0153 0154 %% build B matrices and phase shift injections 0155 [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch); 0156 0157 %% compute complex bus power injections (generation - load) 0158 %% adjusted for phase shifters and real shunts 0159 Pbus = real(makeSbus(baseMVA, bus, gen)) - Pbusinj - bus(:, GS) / baseMVA; 0160 0161 %% "run" the power flow 0162 Va = dcpf(B, Pbus, Va0, ref, pv, pq); 0163 0164 %% update data matrices with solution 0165 branch(:, [QF, QT]) = zeros(size(branch, 1), 2); 0166 branch(:, PF) = (Bf * Va + Pfinj) * baseMVA; 0167 branch(:, PT) = -branch(:, PF); 0168 bus(:, VM) = ones(size(bus, 1), 1); 0169 bus(:, VA) = Va * (180/pi); 0170 %% update Pg for slack generator (1st gen at ref bus) 0171 %% (note: other gens at ref bus are accounted for in Pbus) 0172 %% Pg = Pinj + Pload + Gs 0173 %% newPg = oldPg + newPinj - oldPinj 0174 refgen = zeros(size(ref)); 0175 for k = 1:length(ref) 0176 temp = find(gbus == ref(k)); 0177 refgen(k) = on(temp(1)); 0178 end 0179 gen(refgen, PG) = gen(refgen, PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA; 0180 0181 success = 1; 0182 else %% AC formulation 0183 alg = mpopt(1); 0184 if verbose > 0 0185 if alg == 1 0186 solver = 'Newton'; 0187 elseif alg == 2 0188 solver = 'fast-decoupled, XB'; 0189 elseif alg == 3 0190 solver = 'fast-decoupled, BX'; 0191 elseif alg == 4 0192 solver = 'Gauss-Seidel'; 0193 else 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 V0(gbus) = gen(on, VG) ./ abs(V0(gbus)).* V0(gbus); 0202 0203 if qlim 0204 ref0 = ref; %% save index and angle of 0205 Varef0 = bus(ref0, VA); %% original reference bus(es) 0206 limited = []; %% list of indices of gens @ Q lims 0207 fixedQg = zeros(size(gen, 1), 1); %% Qg of gens at Q limits 0208 end 0209 repeat = 1; 0210 while (repeat) 0211 %% build admittance matrices 0212 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch); 0213 0214 %% compute complex bus power injections (generation - load) 0215 Sbus = makeSbus(baseMVA, bus, gen); 0216 0217 %% run the power flow 0218 alg = mpopt(1); 0219 if alg == 1 0220 [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt); 0221 elseif alg == 2 || alg == 3 0222 [Bp, Bpp] = makeB(baseMVA, bus, branch, alg); 0223 [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt); 0224 elseif alg == 4 0225 [V, success, iterations] = gausspf(Ybus, Sbus, V0, ref, pv, pq, mpopt); 0226 else 0227 error('Only Newton''s method, fast-decoupled, and Gauss-Seidel power flow algorithms currently implemented.'); 0228 end 0229 0230 %% update data matrices with solution 0231 [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq); 0232 0233 if qlim %% enforce generator Q limits 0234 %% find gens with violated Q constraints 0235 mx = find( gen(:, GEN_STATUS) > 0 & gen(:, QG) > gen(:, QMAX) ); 0236 mn = find( gen(:, GEN_STATUS) > 0 & gen(:, QG) < gen(:, QMIN) ); 0237 0238 if ~isempty(mx) || ~isempty(mn) %% we have some Q limit violations 0239 if isempty(pv) 0240 if verbose 0241 if ~isempty(mx) 0242 fprintf('Gen %d (only one left) exceeds upper Q limit : INFEASIBLE PROBLEM\n', mx); 0243 else 0244 fprintf('Gen %d (only one left) exceeds lower Q limit : INFEASIBLE PROBLEM\n', mn); 0245 end 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 verbose && ~isempty(mx) 0265 fprintf('Gen %d at upper Q limit, converting to PQ bus\n', mx); 0266 end 0267 if 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('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 if verbose && ref ~= ref_temp 0293 fprintf('Bus %d is new slack bus\n', ref); 0294 end 0295 limited = [limited; mx]; 0296 else 0297 repeat = 0; %% no more generator Q limits violated 0298 end 0299 else 0300 repeat = 0; %% don't enforce generator Q limits, once is enough 0301 end 0302 end 0303 if qlim && ~isempty(limited) 0304 %% restore injections from limited gens (those at Q limits) 0305 gen(limited, QG) = fixedQg(limited); %% restore Qg value, 0306 for i = 1:length(limited) %% (one at a time, since 0307 bi = gen(limited(i), GEN_BUS); %% they may be at same bus) 0308 bus(bi, [PD,QD]) = ... %% re-adjust load, 0309 bus(bi, [PD,QD]) + gen(limited(i), [PG,QG]); 0310 end 0311 gen(limited, GEN_STATUS) = 1; %% and turn gen back on 0312 if ref ~= ref0 0313 %% adjust voltage angles to make original ref bus correct 0314 bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0; 0315 end 0316 end 0317 end 0318 mpc.et = etime(clock, t0); 0319 mpc.success = success; 0320 0321 %%----- output results ----- 0322 %% convert back to original bus numbering & print results 0323 [mpc.bus, mpc.gen, mpc.branch] = deal(bus, gen, branch); 0324 results = int2ext(mpc); 0325 0326 %% zero out result fields of out-of-service gens & branches 0327 if ~isempty(results.order.gen.status.off) 0328 results.gen(results.order.gen.status.off, [PG QG]) = 0; 0329 end 0330 if ~isempty(results.order.branch.status.off) 0331 results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0; 0332 end 0333 0334 if fname 0335 [fd, msg] = fopen(fname, 'at'); 0336 if fd == -1 0337 error(msg); 0338 else 0339 printpf(results, fd, mpopt); 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