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.25 2010/04/26 19:45:25 ray 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 swing generator (note: other gens at ref bus are accounted for in Pbus) 0171 %% Pg = Pinj + Pload + Gs 0172 %% newPg = oldPg + newPinj - oldPinj 0173 refgen = find(gbus == ref); %% which is(are) the reference gen(s)? 0174 gen(on(refgen(1)), PG) = gen(on(refgen(1)), PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA; 0175 0176 success = 1; 0177 else %% AC formulation 0178 if verbose > 0 0179 fprintf(' -- AC Power Flow '); %% solver name and \n added later 0180 end 0181 %% initial state 0182 % V0 = ones(size(bus, 1), 1); %% flat start 0183 V0 = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA)); 0184 V0(gbus) = gen(on, VG) ./ abs(V0(gbus)).* V0(gbus); 0185 0186 if qlim 0187 ref0 = ref; %% save index and angle of 0188 Varef0 = bus(ref0, VA); %% original reference bus 0189 limited = []; %% list of indices of gens @ Q lims 0190 fixedQg = zeros(size(gen, 1), 1); %% Qg of gens at Q limits 0191 end 0192 repeat = 1; 0193 while (repeat) 0194 %% build admittance matrices 0195 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch); 0196 0197 %% compute complex bus power injections (generation - load) 0198 Sbus = makeSbus(baseMVA, bus, gen); 0199 0200 %% run the power flow 0201 alg = mpopt(1); 0202 if alg == 1 0203 [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt); 0204 elseif alg == 2 || alg == 3 0205 [Bp, Bpp] = makeB(baseMVA, bus, branch, alg); 0206 [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt); 0207 elseif alg == 4 0208 [V, success, iterations] = gausspf(Ybus, Sbus, V0, ref, pv, pq, mpopt); 0209 else 0210 error('Only Newton''s method, fast-decoupled, and Gauss-Seidel power flow algorithms currently implemented.'); 0211 end 0212 0213 %% update data matrices with solution 0214 [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq); 0215 0216 if qlim %% enforce generator Q limits 0217 %% find gens with violated Q constraints 0218 mx = find( gen(:, GEN_STATUS) > 0 & gen(:, QG) > gen(:, QMAX) ); 0219 mn = find( gen(:, GEN_STATUS) > 0 & gen(:, QG) < gen(:, QMIN) ); 0220 0221 if ~isempty(mx) || ~isempty(mn) %% we have some Q limit violations 0222 if isempty(pv) 0223 if verbose 0224 if ~isempty(mx) 0225 fprintf('Gen %d (only one left) exceeds upper Q limit : INFEASIBLE PROBLEM\n', mx); 0226 else 0227 fprintf('Gen %d (only one left) exceeds lower Q limit : INFEASIBLE PROBLEM\n', mn); 0228 end 0229 end 0230 success = 0; 0231 break; 0232 end 0233 0234 %% one at a time? 0235 if qlim == 2 %% fix largest violation, ignore the rest 0236 [junk, k] = max([gen(mx, QG) - gen(mx, QMAX); 0237 gen(mn, QMIN) - gen(mn, QG)]); 0238 if k > length(mx) 0239 mn = mn(k-length(mx)); 0240 mx = []; 0241 else 0242 mx = mx(k); 0243 mn = []; 0244 end 0245 end 0246 0247 if verbose && ~isempty(mx) 0248 fprintf('Gen %d at upper Q limit, converting to PQ bus\n', mx); 0249 end 0250 if verbose && ~isempty(mn) 0251 fprintf('Gen %d at lower Q limit, converting to PQ bus\n', mn); 0252 end 0253 0254 %% save corresponding limit values 0255 fixedQg(mx) = gen(mx, QMAX); 0256 fixedQg(mn) = gen(mn, QMIN); 0257 mx = [mx;mn]; 0258 0259 %% convert to PQ bus 0260 gen(mx, QG) = fixedQg(mx); %% set Qg to binding limit 0261 gen(mx, GEN_STATUS) = 0; %% temporarily turn off gen, 0262 for i = 1:length(mx) %% (one at a time, since 0263 bi = gen(mx(i), GEN_BUS); %% they may be at same bus) 0264 bus(bi, [PD,QD]) = ... %% adjust load accordingly, 0265 bus(bi, [PD,QD]) - gen(mx(i), [PG,QG]); 0266 end 0267 bus(gen(mx, GEN_BUS), BUS_TYPE) = PQ; %% & set bus type to PQ 0268 0269 %% update bus index lists of each type of bus 0270 ref_temp = ref; 0271 [ref, pv, pq] = bustypes(bus, gen); 0272 if verbose && ref ~= ref_temp 0273 fprintf('Bus %d is new slack bus\n', ref); 0274 end 0275 limited = [limited; mx]; 0276 else 0277 repeat = 0; %% no more generator Q limits violated 0278 end 0279 else 0280 repeat = 0; %% don't enforce generator Q limits, once is enough 0281 end 0282 end 0283 if qlim && ~isempty(limited) 0284 %% restore injections from limited gens (those at Q limits) 0285 gen(limited, QG) = fixedQg(limited); %% restore Qg value, 0286 for i = 1:length(limited) %% (one at a time, since 0287 bi = gen(limited(i), GEN_BUS); %% they may be at same bus) 0288 bus(bi, [PD,QD]) = ... %% re-adjust load, 0289 bus(bi, [PD,QD]) + gen(limited(i), [PG,QG]); 0290 end 0291 gen(limited, GEN_STATUS) = 1; %% and turn gen back on 0292 if ref ~= ref0 0293 %% adjust voltage angles to make original ref bus correct 0294 bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0; 0295 end 0296 end 0297 end 0298 mpc.et = etime(clock, t0); 0299 mpc.success = success; 0300 0301 %%----- output results ----- 0302 %% convert back to original bus numbering & print results 0303 [mpc.bus, mpc.gen, mpc.branch] = deal(bus, gen, branch); 0304 results = int2ext(mpc); 0305 0306 %% zero out result fields of out-of-service gens & branches 0307 if ~isempty(results.order.gen.status.off) 0308 results.gen(results.order.gen.status.off, [PG QG]) = 0; 0309 end 0310 if ~isempty(results.order.branch.status.off) 0311 results.branch(results.order.branch.status.off, [PF QF PT QT]) = 0; 0312 end 0313 0314 if fname 0315 [fd, msg] = fopen(fname, 'at'); 0316 if fd == -1 0317 error(msg); 0318 else 0319 printpf(results, fd, mpopt); 0320 fclose(fd); 0321 end 0322 end 0323 printpf(results, 1, mpopt); 0324 0325 %% save solved case 0326 if solvedcase 0327 savecase(solvedcase, results); 0328 end 0329 0330 if nargout == 1 || nargout == 2 0331 MVAbase = results; 0332 bus = success; 0333 elseif nargout > 2 0334 [MVAbase, bus, gen, branch, et] = ... 0335 deal(results.baseMVA, results.bus, results.gen, results.branch, results.et); 0336 % else %% don't define MVAbase, so it doesn't print anything 0337 end