半不变量法概率潮流_基于半不变量法的概率潮流计算_概率潮流

function [MVAbase, bus, gen, branch, success, et, dSbus_dVm, dSbus_dVa] = ... runpf(casedata, mpopt, fname, solvedcase) %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. % 运行潮流（默认情况下为完全AC牛顿方法），可以选择返回RESULTS结构和SUCCESS标志。 % 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. % MATPOWER % Copyright (c) 1996-2016 by Power System Engineering Research Center (PSERC) % by Ray Zimmerman, PSERC Cornell % Enforcing of generator Q limits inspired by contributions % from Mu Lin, Lincoln University, New Zealand (1/14/05). % % This file is part of MATPOWER. % Covered by the 3-clause BSD License (see LICENSE file for details). % 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 arguments if 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 casedata = 'case9'; %% default data file is 'case9.m' end end end end %% options qlim = mpopt.pf.enforce_q_lims; %% enforce Q limits on gens? dc = strcmp(upper(mpopt.model), 'DC'); %% use DC formulation? %% read data mpc = loadcase(casedata); %% add zero columns to branch for flows if needed if size(mpc.branch,2) < QT mpc.branch = [ mpc.branch zeros(size(mpc.branch, 1), QT-size(mpc.branch,2)) ]; end %% convert to internal indexing mpc = ext2int(mpc); [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch); %% get bus index lists of each type of bus [ref, pv, pq] = bustypes(bus, gen); %% generator info on = 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 mpopt.verbose > 0 v = mpver('all'); fprintf('\nMATPOWER Version %s, %s', v.Version, v.Date); end if dc %% DC formulation if mpopt.verbose > 0 fprintf(' -- DC Power Flow\n'); end %% 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, success] = 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 slack generator (1st gen at ref bus) %% (note: other gens at ref bus are accounted for in Pbus) %% Pg = Pinj + Pload + Gs %% newPg = oldPg + newPinj - oldPinj refgen = zeros(size(ref)); for k = 1:length(ref) temp = find(gbus == ref(k)); refgen(k) = on(temp(1)); end gen(refgen, PG) = gen(refgen, PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA; else %% AC formulation alg = upper(mpopt.pf.alg); if mpopt.verbose > 0 switch alg case 'NR' solver = 'Newton'; case 'FDXB' solver = 'fast-decoupled, XB'; case 'FDBX' solver = 'fast-decoupled, BX'; case 'GS' solver = 'Gauss-Seidel'; otherwise solver = 'unknown'; end fprintf(' -- AC Power Flow (%s)\n', solver); end %% initial state % V0 = ones(size(bus, 1), 1); %% flat start V0 = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA)); vcb = ones(size(V0)); %% create mask of voltage-controlled buses vcb(pq) = 0; %% exclude PQ buses k = find(vcb(gbus)); %% in-service gens at v-c buses V0(gbus(k)) = gen(on(k), VG) ./ abs(V0(gbus(k))).* V0(gbus(k)); if qlim ref0 = ref; %% save index and angle of Varef0 = bus(ref0, VA); %% original reference bus(es) limited = []; %% list of indices of gens @ Q lims fixedQg = zeros(size(gen, 1), 1); %% Qg of gens at Q limits end %% build admittance matrices [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch); repeat = 1; while (repeat) %% function for computing V dependent complex bus power inj