风电机组并网潮流计算的程序

function [pfitness] = ... runpf1(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. % % 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-2020, Power Systems 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 https://matpower.org 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 = 'case014'; %% default data file is 'case14.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); %mpc.bus(1,8) =1.065; %mpc.gen(1,6) =1.065; %mpc.bus(2,8) =1.0468; %mpc.gen(2,6) =1.0468; %mpc.bus(3,8) =1.009; %mpc.gen(3,6) =1.009; %mpc.bus(6,8) =1.0377; %mpc.gen(4,6) =1.0377; %mpc.bus(8,8) =1.0413; %mpc.gen(5,6) =1.0413; %mpc.bus(9,6) =0.1; %mpc.branch(8,9) =1; %mpc.branch(9,9) =1.02; %mpc.branch(10,9) =0.96; %% 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, mpopt); [baseMVA, bus, gen, branch] = deal(mpc.baseMVA, mpc.bus, mpc.gen, mpc.branch); %在这里添加风电场进来 NO_NUM=14; v=11;%风速 n=10;%风机台数 [Pe_WF,Ps,Pr,Qe,W]=QE_F(mpc.bus(NO_NUM,8),v,n); B_gen =[NO_NUM+1,(Pe_WF)/1000,Qe/1000,40,0,1.036,100,1,100,0,0,0,0,0,0,0,0,0,0,0,0]; gen(3,2)= gen(3,2)+(Pe_WF)/1000; gen(3,3)= gen(3,3)+Qe/1000; if ~isempty(mpc.bus) %% 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 = tic; its = 0; %% total iterations 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); its = 1; %% 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); switch alg case 'NR-SP' mpopt = mpoption(mpopt, 'pf.current_balance', 0, 'pf.v_cartesian', 0); case 'NR-SC' mpopt = mpoption(mpopt, 'pf.current_balance', 0, 'pf.v_cartesian', 1); case 'NR-SH' mpopt = mpoption(mpopt, 'pf.current_balance', 0, 'pf.v_cartesian', 2); case 'NR-IP' mpopt = mpoption(mpopt, 'pf.current_balance', 1, 'pf.v_cartesian', 0); case 'NR-IC' mpopt = m