function [MVAbase, bus, gen, branch, success, et] = ...    			runpf(casename, mpopt, fname, solvedcase)%RUNPF  Runs a power flow.%%   [baseMVA, bus, gen, branch, success, et] = ...%           runpf(casename, mpopt, fname, solvedcase)%%   Runs a power flow (full AC Newton's method by default) and optionally%   returns the solved values in the data matrices, a flag which is true if%   the algorithm was successful in finding a solution, and the elapsed time%   in seconds. All input arguments are optional. If casename is provided it%   specifies the name of the input data file or struct (see also 'help%   caseformat' and 'help loadcase') containing the power flow data. The%   default value is 'case9'. If the mpopt is provided it overrides the%   default MATPOWER options vector and can be used to specify the solution%   algorithm and output options among other things (see 'help mpoption' for%   details). If the 3rd argument is given the pretty printed output will be%   appended to the file whose name is given in fname. If solvedcase is%   specified the solved case will be written to a case file in MATPOWER%   format with the specified name. If solvedcase ends with '.mat' it saves%   the case as a MAT-file otherwise it saves it as an M-file.%%   If the ENFORCE_Q_LIMS options 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.%   MATPOWER%   $Id: runpf.m,v 1.14 2006/09/29 19:23:07 ray Exp $%   by Ray Zimmerman, PSERC Cornell%   Enforcing of generator Q limits inspired by contributions%   from Mu Lin, Lincoln University, New Zealand (1/14/05).%   Copyright (c) 1996-2005 by Power System Engineering Research Center (PSERC)%   See http://www.pserc.cornell.edu/matpower/ for more info.%%-----  initialize  -----%% define named indices into bus, gen, branch matrices[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;%% default argumentsif nargin < 4    solvedcase = '';                %% don't save solved case    if nargin < 3        fname = '';                 %% don't print results to a file        if nargin < 2            mpopt = mpoption;       %% use default options            if nargin < 1                casename = 'case9'; %% default data file is 'case9.m'            end        end    endend%% optionsverbose = mpopt(31);qlim = mpopt(6);                    %% enforce Q limits on gens?dc = mpopt(10);                     %% use DC formulation?%% read data & convert to internal bus numbering[baseMVA, bus, gen, branch] = loadcase(casename);[i2e, bus, gen, branch] = ext2int(bus, gen, branch);%% get bus index lists of each type of bus[ref, pv, pq] = bustypes(bus, gen);%% generator infoon = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?gbus = gen(on, GEN_BUS);                %% what buses are they at?%%-----  run the power flow  -----t0 = clock;if dc                               %% DC formulation    %% initial state    Va0 = bus(:, VA) * (pi/180);        %% build B matrices and phase shift injections    [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);        %% compute complex bus power injections (generation - load)    %% adjusted for phase shifters and real shunts    Pbus = real(makeSbus(baseMVA, bus, gen)) - Pbusinj - bus(:, GS) / baseMVA;        %% "run" the power flow    Va = dcpf(B, Pbus, Va0, ref, pv, pq);        %% update data matrices with solution    branch(:, [QF, QT]) = zeros(size(branch, 1), 2);    branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;    branch(:, PT) = -branch(:, PF);    bus(:, VM) = ones(size(bus, 1), 1);    bus(:, VA) = Va * (180/pi);    %% update Pg for swing generator (note: other gens at ref bus are accounted for in Pbus)    %%      Pg = Pinj + Pload + Gs    %%      newPg = oldPg + newPinj - oldPinj    refgen = find(gbus == ref);             %% which is(are) the reference gen(s)?    gen(on(refgen(1)), PG) = gen(on(refgen(1)), PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA;        success = 1;else                                %% AC formulation    %% initial state    % V0    = ones(size(bus, 1), 1);            %% flat start    V0  = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA));    V0(gbus) = gen(on, VG) ./ abs(V0(gbus)).* V0(gbus);        if qlim        ref0 = ref;                         %% save index and angle of        Varef0 = bus(ref0, VA);             %%   original reference bus        limited = [];                       %% list of indices of gens @ Q lims        fixedQg = zeros(size(gen, 1), 1);   %% Qg of gens at Q limits    end    repeat = 1;    while (repeat)        %% build admittance matrices        [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);                %% compute complex bus power injections (generation - load)        Sbus = makeSbus(baseMVA, bus, gen);                %% run the power flow        alg = mpopt(1);        if alg == 1            [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt);        elseif alg == 2 | alg == 3            [Bp, Bpp] = makeB(baseMVA, bus, branch, alg);            [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt);        elseif alg == 4            [V, success, iterations] = gausspf(Ybus, Sbus, V0, ref, pv, pq, mpopt);        else            error('Only Newton''s method, fast-decoupled, and Gauss-Seidel power flow algorithms currently implemented.');        end                %% update data matrices with solution        [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq);                if qlim             %% enforce generator Q limits            %% find gens with violated Q constraints            mx = find( gen(:, GEN_STATUS) > 0 & gen(:, QG) > gen(:, QMAX) );            mn = find( gen(:, GEN_STATUS) > 0 & gen(:, QG) < gen(:, QMIN) );                        if ~isempty(mx) | ~isempty(mn)  %% we have some Q limit violations                if isempty(pv)                	if verbose                		if ~isempty(mx)			                	fprintf('Gen %d (only one left) exceeds upper Q limit : INFEASIBLE PROBLEM\n', mx);                		else		                	fprintf('Gen %d (only one left) exceeds lower Q limit : INFEASIBLE PROBLEM\n', mn);                		end	                end                	success = 0;                	break;                end                if verbose & ~isempty(mx)                    fprintf('Gen %d at upper Q limit, converting to PQ bus\n', mx);                end                if verbose & ~isempty(mn)                    fprintf('Gen %d at lower Q limit, converting to PQ bus\n', mn);                end                                %% save corresponding limit values                fixedQg(mx) = gen(mx, QMAX);                fixedQg(mn) = gen(mn, QMIN);                mx = [mx;mn];                                %% convert to PQ bus                gen(mx, QG) = fixedQg(mx);      %% set Qg to binding limit                gen(mx, GEN_STATUS) = 0;        %% temporarily turn off gen,                for i = 1:length(mx)            %% (one at a time, since                    bi = gen(mx(i), GEN_BUS);   %%  they may be at same bus)                    bus(bi, [PD,QD]) = ...      %% adjust load accordingly,                        bus(bi, [PD,QD]) - gen(mx(i), [PG,QG]);                end                bus(gen(mx, GEN_BUS), BUS_TYPE) = PQ;   %% & set bus type to PQ                                %% update bus index lists of each type of bus                ref_temp = ref;                [ref, pv, pq] = bustypes(bus, gen);                if verbose & ref ~= ref_temp                    fprintf('Bus %d is new slack bus\n', ref);                end                limited = [limited; mx];            else                repeat = 0; %% no more generator Q limits violated            end        else            repeat = 0;     %% don't enforce generator Q limits, once is enough        end    end    if qlim & ~isempty(limited)        %% restore injections from limited gens (those at Q limits)        gen(limited, QG) = fixedQg(limited);    %% restore Qg value,        for i = 1:length(limited)               %% (one at a time, since            bi = gen(limited(i), GEN_BUS);      %%  they may be at same bus)            bus(bi, [PD,QD]) = ...              %% re-adjust load,                bus(bi, [PD,QD]) + gen(limited(i), [PG,QG]);        end        gen(limited, GEN_STATUS) = 1;               %% and turn gen back on        if ref ~= ref0            %% adjust voltage angles to make original ref bus correct            bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0;        end    endendet = etime(clock, t0);%%-----  output results  -----%% convert back to original bus numbering & print results[bus, gen, branch] = int2ext(i2e, bus, gen, branch);if fname    [fd, msg] = fopen(fname, 'at');    if fd == -1        error(msg);    else        printpf(baseMVA, bus, gen, branch, [], success, et, fd, mpopt);        fclose(fd);    endendprintpf(baseMVA, bus, gen, branch, [], success, et, 1, mpopt);%% save solved caseif solvedcase    savecase(solvedcase, baseMVA, bus, gen, branch);end%% this is just to prevent it from printing baseMVA%% when called with no output argumentsif nargout, MVAbase = baseMVA; endreturn;