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 2397 2014-10-16 20:37:40Z ray $ 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 gens at PV/slack buses and Pg for slack bus(es) ----- 0055 %% generator info 0056 on = find(gen(:, GEN_STATUS) > 0 & ... %% which generators are on? 0057 bus(gen(:, GEN_BUS), BUS_TYPE) ~= PQ); %% ... and not at PQ buses 0058 off = find(gen(:, GEN_STATUS) <= 0); %% which generators are off? 0059 gbus = gen(on, GEN_BUS); %% what buses are they at? 0060 0061 %% compute total injected bus powers 0062 Sbus = V(gbus) .* conj(Ybus(gbus, :) * V); 0063 0064 %% update Qg for generators at PV/slack buses 0065 gen(off, QG) = zeros(length(off), 1); %% zero out off-line Qg 0066 %% don't touch the ones at PQ buses 0067 gen(on, QG) = imag(Sbus) * baseMVA + bus(gbus, QD); %% inj Q + local Qd 0068 %% ... at this point any buses with more than one generator will have 0069 %% the total Q dispatch for the bus assigned to each generator. This 0070 %% must be split between them. We do it first equally, then in proportion 0071 %% to the reactive range of the generator. 0072 0073 if length(on) > 1 0074 %% build connection matrix, element i, j is 1 if gen on(i) at bus j is ON 0075 nb = size(bus, 1); 0076 ngon = size(on, 1); 0077 Cg = sparse((1:ngon)', gbus, ones(ngon, 1), ngon, nb); 0078 0079 %% divide Qg by number of generators at the bus to distribute equally 0080 ngg = Cg * sum(Cg)'; %% ngon x 1, number of gens at this gen's bus 0081 gen(on, QG) = gen(on, QG) ./ ngg; 0082 0083 %% divide proportionally 0084 Cmin = sparse((1:ngon)', gbus, gen(on, QMIN), ngon, nb); 0085 Cmax = sparse((1:ngon)', gbus, gen(on, QMAX), ngon, nb); 0086 Qg_tot = Cg' * gen(on, QG); %% nb x 1 vector of total Qg at each bus 0087 Qg_min = sum(Cmin)'; %% nb x 1 vector of min total Qg at each bus 0088 Qg_max = sum(Cmax)'; %% nb x 1 vector of max total Qg at each bus 0089 ig = find(Cg * Qg_min == Cg * Qg_max); %% gens at buses with Qg range = 0 0090 Qg_save = gen(on(ig), QG); 0091 gen(on, QG) = gen(on, QMIN) + ... 0092 (Cg * ((Qg_tot - Qg_min)./(Qg_max - Qg_min + eps))) .* ... 0093 (gen(on, QMAX) - gen(on, QMIN)); %% ^ avoid div by 0 0094 gen(on(ig), QG) = Qg_save; 0095 end %% (terms are mult by 0 anyway) 0096 0097 %% update Pg for slack gen(s) 0098 for k = 1:length(ref) 0099 refgen = find(gbus == ref(k)); %% which is(are) the reference gen(s)? 0100 gen(on(refgen(1)), PG) = real(Sbus(refgen(1))) * baseMVA ... 0101 + bus(ref(k), PD); %% inj P + local Pd 0102 if length(refgen) > 1 %% more than one generator at this ref bus 0103 %% subtract off what is generated by other gens at this bus 0104 gen(on(refgen(1)), PG) = gen(on(refgen(1)), PG) ... 0105 - sum(gen(on(refgen(2:length(refgen))), PG)); 0106 end 0107 end 0108 0109 %%----- update/compute branch power flows ----- 0110 out = find(branch(:, BR_STATUS) == 0); %% out-of-service branches 0111 br = find(branch(:, BR_STATUS)); %% in-service branches 0112 Sf = V(branch(br, F_BUS)) .* conj(Yf(br, :) * V) * baseMVA; %% complex power at "from" bus 0113 St = V(branch(br, T_BUS)) .* conj(Yt(br, :) * V) * baseMVA; %% complex power injected at "to" bus 0114 branch(br, [PF, QF, PT, QT]) = [real(Sf) imag(Sf) real(St) imag(St)]; 0115 branch(out, [PF, QF, PT, QT]) = zeros(length(out), 4);