[BL2, BL1, BL0] = MAKEBLOSS(MPC) MPC must have consecutive bus numbering (i.e. internal bus numbering) LOSS = (0.5 * Pg' * BL2 * Pg + BL1' * Pg + BL0 ) Example: mpc = loadcase('case30'); results = rundcopf(mpc); Pg = results.gen(:, PG); [BL2, BL1, BL0] = makeBloss2(mpc); loss = (0.5 * Pg' * BL2 * Pg + BL1' * Pg + BL0 ) Pf = results.branch(:, PF) / mpc.baseMVA; r = mpc.branch(:, BR_R); loss2 = sum( Pf .^ 2 .* r ) * mpc.baseMVA
0001 function [BL2, BL1, BL0] = makeBloss(mpc) 0002 % [BL2, BL1, BL0] = MAKEBLOSS(MPC) 0003 % 0004 % MPC must have consecutive bus numbering (i.e. internal bus numbering) 0005 % 0006 % LOSS = (0.5 * Pg' * BL2 * Pg + BL1' * Pg + BL0 ) 0007 % 0008 % Example: 0009 % mpc = loadcase('case30'); 0010 % results = rundcopf(mpc); 0011 % Pg = results.gen(:, PG); 0012 % [BL2, BL1, BL0] = makeBloss2(mpc); 0013 % loss = (0.5 * Pg' * BL2 * Pg + BL1' * Pg + BL0 ) 0014 % 0015 % Pf = results.branch(:, PF) / mpc.baseMVA; 0016 % r = mpc.branch(:, BR_R); 0017 % loss2 = sum( Pf .^ 2 .* r ) * mpc.baseMVA 0018 0019 % MATPOWER 0020 % $Id: makeBloss.m 2471 2014-12-16 15:32:35Z ray $ 0021 % by Ray Zimmerman, PSERC Cornell 0022 % Copyright (c) 2014 by Power System Engineering Research Center (PSERC) 0023 % 0024 % This file is part of MATPOWER. 0025 % See http://www.pserc.cornell.edu/matpower/ for more info. 0026 % 0027 % MATPOWER is free software: you can redistribute it and/or modify 0028 % it under the terms of the GNU General Public License as published 0029 % by the Free Software Foundation, either version 3 of the License, 0030 % or (at your option) any later version. 0031 % 0032 % MATPOWER is distributed in the hope that it will be useful, 0033 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0034 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0035 % GNU General Public License for more details. 0036 % 0037 % You should have received a copy of the GNU General Public License 0038 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0039 % 0040 % Additional permission under GNU GPL version 3 section 7 0041 % 0042 % If you modify MATPOWER, or any covered work, to interface with 0043 % other modules (such as MATLAB code and MEX-files) available in a 0044 % MATLAB(R) or comparable environment containing parts covered 0045 % under other licensing terms, the licensors of MATPOWER grant 0046 % you additional permission to convey the resulting work. 0047 0048 %% define named indices into data matrices 0049 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0050 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0051 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ... 0052 RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch; 0053 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ... 0054 GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen; 0055 0056 if ischar(mpc) 0057 mpc = loadcase(mpc); 0058 end 0059 [Bbus, Bf, Pbusinj, Pfinj] = makeBdc(mpc.baseMVA, mpc.bus, mpc.branch); 0060 0061 %% form generator connection matrix: element i,j is 1 if gen j at bus i is ON 0062 nb = size(mpc.bus, 1); 0063 ng = size(mpc.gen, 1); 0064 nl = size(mpc.branch, 1); 0065 Cg = sparse(mpc.gen(:, GEN_BUS), (1:ng)', ones(ng, 1), nb, ng); 0066 0067 %% get indices k of non-reference buses 0068 ref = find(mpc.bus(:, BUS_TYPE) == REF); 0069 ref = ref(1); 0070 k = (1:nb)'; 0071 k(ref) = []; 0072 0073 %% Pbus = Cg * Pg - Pd - Ps 0074 %% Pbus(k) = Cg(k,:) * Pg - Pd(k) - Ps(k) 0075 %% Pbus(k) = Bbus(k,k) * Va(k) + Pbusinj(k) 0076 %% Va(k) = inv( Bbus(k,k) ) * (Pbus(k) - Pbusinj(k)) 0077 %% Pf = Bf * Va + Pfinj 0078 %% Pf = Bf(:,k) * Va(k) + Bf(:,ref) * Va(ref) + Pfinj 0079 %% = Bf(:,k) * inv( Bbus(k,k) ) * (Pbus(k) - Pbusinj(k)) + Bf(:,ref) * Va(ref) + Pfinj 0080 %% = W * (Pbus(k) - Pbusinj(k)) + Bf(:,ref) * Va(ref) + Pfinj 0081 %% = W * (Cg(k,:) * Pg - Pd(k) - Ps(k) - Pbusinj(k)) + Bf(:,ref) * Va(ref) + Pfinj 0082 %% = Wg * Pg + Pz; 0083 %% W = Bf(:,k) * inv( Bbus(k,k) ); 0084 %% Wg = W * Cg(k,:); 0085 %% Pz = Bf(:,ref) * Va(ref) + Pfinj - W * (Pd(k) + Ps(k) + Pbusinj(k)); 0086 %% loss = Pf' * R * Pf 0087 %% = (Wg * Pg + Pz)' * R * (Wg * Pg + Pz) 0088 %% = Pg' * Wg' * R * Wg * Pg 0089 %% + 2 * Pz' * R * Wg * Pg 0090 %% + Pz' * R * Pz 0091 0092 r = mpc.branch(:, BR_R); 0093 Varef = mpc.bus(ref, VA) * pi/180; 0094 Pd = mpc.bus(:, PD) / mpc.baseMVA; 0095 Ps = mpc.bus(:, GS) / mpc.baseMVA; 0096 W = Bf(:, k) / Bbus(k, k); %% = Bf(:, k) * inv(Bbus(k, k)) 0097 Wg = W * Cg(k,:); 0098 Pz = Bf(:,ref) * Varef + Pfinj - W * (Pd(k) + Ps(k) + Pbusinj(k)); 0099 R = sparse(1:nl, 1:nl, r, nl, nl); %% nl x nl matrix w/ line resistences on diag 0100 0101 %% form B coefficients 0102 BL2 = 2 * Wg' * R * Wg / mpc.baseMVA; 0103 BL1 = 2 * Wg' * R * Pz; 0104 BL0 = Pz' * R * Pz * mpc.baseMVA;