PFSOLN Updates bus, gen, branch data structures to match power flow soln.
[BUS, GEN, BRANCH] = PFSOLN(BASEMVA, BUS0, GEN0, BRANCH0, ...
YBUS, YF, YT, V, REF, PV, PQ)0001 function [bus, gen, branch] = pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq) 0002 %PFSOLN Updates bus, gen, branch data structures to match power flow soln. 0003 % [BUS, GEN, BRANCH] = PFSOLN(BASEMVA, BUS0, GEN0, BRANCH0, ... 0004 % YBUS, YF, YT, V, REF, PV, PQ) 0005 0006 % MATPOWER 0007 % $Id: pfsoln.m,v 1.20 2011/05/17 15:16:05 cvs Exp $ 0008 % by Ray Zimmerman, PSERC Cornell 0009 % Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC) 0010 % 0011 % This file is part of MATPOWER. 0012 % See http://www.pserc.cornell.edu/matpower/ for more info. 0013 % 0014 % MATPOWER is free software: you can redistribute it and/or modify 0015 % it under the terms of the GNU General Public License as published 0016 % by the Free Software Foundation, either version 3 of the License, 0017 % or (at your option) any later version. 0018 % 0019 % MATPOWER is distributed in the hope that it will be useful, 0020 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0021 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0022 % GNU General Public License for more details. 0023 % 0024 % You should have received a copy of the GNU General Public License 0025 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0026 % 0027 % Additional permission under GNU GPL version 3 section 7 0028 % 0029 % If you modify MATPOWER, or any covered work, to interface with 0030 % other modules (such as MATLAB code and MEX-files) available in a 0031 % MATLAB(R) or comparable environment containing parts covered 0032 % under other licensing terms, the licensors of MATPOWER grant 0033 % you additional permission to convey the resulting work. 0034 0035 %% define named indices into bus, gen, branch matrices 0036 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0037 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0038 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... 0039 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... 0040 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; 0041 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... 0042 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... 0043 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; 0044 0045 %% initialize return values 0046 bus = bus0; 0047 gen = gen0; 0048 branch = branch0; 0049 0050 %%----- update bus voltages ----- 0051 bus(:, VM) = abs(V); 0052 bus(:, VA) = angle(V) * 180 / pi; 0053 0054 %%----- update Qg for all gens and Pg for slack bus(es) ----- 0055 %% generator info 0056 on = find(gen(:, GEN_STATUS) > 0); %% which generators are on? 0057 gbus = gen(on, GEN_BUS); %% what buses are they at? 0058 0059 %% compute total injected bus powers 0060 Sbus = V(gbus) .* conj(Ybus(gbus, :) * V); 0061 0062 %% update Qg for all generators 0063 gen(:, QG) = zeros(size(gen, 1), 1); %% zero out all Qg 0064 gen(on, QG) = imag(Sbus) * baseMVA + bus(gbus, QD); %% inj Q + local Qd 0065 %% ... at this point any buses with more than one generator will have 0066 %% the total Q dispatch for the bus assigned to each generator. This 0067 %% must be split between them. We do it first equally, then in proportion 0068 %% to the reactive range of the generator. 0069 0070 if length(on) > 1 0071 %% build connection matrix, element i, j is 1 if gen on(i) at bus j is ON 0072 nb = size(bus, 1); 0073 ngon = size(on, 1); 0074 Cg = sparse((1:ngon)', gbus, ones(ngon, 1), ngon, nb); 0075 0076 %% divide Qg by number of generators at the bus to distribute equally 0077 ngg = Cg * sum(Cg)'; %% ngon x 1, number of gens at this gen's bus 0078 gen(on, QG) = gen(on, QG) ./ ngg; 0079 0080 %% divide proportionally 0081 Cmin = sparse((1:ngon)', gbus, gen(on, QMIN), ngon, nb); 0082 Cmax = sparse((1:ngon)', gbus, gen(on, QMAX), ngon, nb); 0083 Qg_tot = Cg' * gen(on, QG); %% nb x 1 vector of total Qg at each bus 0084 Qg_min = sum(Cmin)'; %% nb x 1 vector of min total Qg at each bus 0085 Qg_max = sum(Cmax)'; %% nb x 1 vector of max total Qg at each bus 0086 ig = find(Cg * Qg_min == Cg * Qg_max); %% gens at buses with Qg range = 0 0087 Qg_save = gen(on(ig), QG); 0088 gen(on, QG) = gen(on, QMIN) + ... 0089 (Cg * ((Qg_tot - Qg_min)./(Qg_max - Qg_min + eps))) .* ... 0090 (gen(on, QMAX) - gen(on, QMIN)); %% ^ avoid div by 0 0091 gen(on(ig), QG) = Qg_save; 0092 end %% (terms are mult by 0 anyway) 0093 0094 %% update Pg for slack gen(s) 0095 for k = 1:length(ref) 0096 refgen = find(gbus == ref(k)); %% which is(are) the reference gen(s)? 0097 gen(on(refgen(1)), PG) = real(Sbus(refgen(1))) * baseMVA ... 0098 + bus(ref(k), PD); %% inj P + local Pd 0099 if length(refgen) > 1 %% more than one generator at this ref bus 0100 %% subtract off what is generated by other gens at this bus 0101 gen(on(refgen(1)), PG) = gen(on(refgen(1)), PG) ... 0102 - sum(gen(on(refgen(2:length(refgen))), PG)); 0103 end 0104 end 0105 0106 %%----- update/compute branch power flows ----- 0107 out = find(branch(:, BR_STATUS) == 0); %% out-of-service branches 0108 br = find(branch(:, BR_STATUS)); %% in-service branches 0109 Sf = V(branch(br, F_BUS)) .* conj(Yf(br, :) * V) * baseMVA; %% complex power at "from" bus 0110 St = V(branch(br, T_BUS)) .* conj(Yt(br, :) * V) * baseMVA; %% complex power injected at "to" bus 0111 branch(br, [PF, QF, PT, QT]) = [real(Sf) imag(Sf) real(St) imag(St)]; 0112 branch(out, [PF, QF, PT, QT]) = zeros(length(out), 4);