[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 % Copyright (c) 2014-2016, Power Systems Engineering Research Center (PSERC) 0021 % by Ray Zimmerman, PSERC Cornell 0022 % 0023 % This file is part of MATPOWER. 0024 % Covered by the 3-clause BSD License (see LICENSE file for details). 0025 % See http://www.pserc.cornell.edu/matpower/ for more info. 0026 0027 %% define named indices into data matrices 0028 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... 0029 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; 0030 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ... 0031 RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch; 0032 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ... 0033 GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen; 0034 0035 if ischar(mpc) 0036 mpc = loadcase(mpc); 0037 end 0038 [Bbus, Bf, Pbusinj, Pfinj] = makeBdc(mpc.baseMVA, mpc.bus, mpc.branch); 0039 0040 %% form generator connection matrix: element i,j is 1 if gen j at bus i is ON 0041 nb = size(mpc.bus, 1); 0042 ng = size(mpc.gen, 1); 0043 nl = size(mpc.branch, 1); 0044 Cg = sparse(mpc.gen(:, GEN_BUS), (1:ng)', ones(ng, 1), nb, ng); 0045 0046 %% get indices k of non-reference buses 0047 ref = find(mpc.bus(:, BUS_TYPE) == REF); 0048 ref = ref(1); 0049 k = (1:nb)'; 0050 k(ref) = []; 0051 0052 %% Pbus = Cg * Pg - Pd - Ps 0053 %% Pbus(k) = Cg(k,:) * Pg - Pd(k) - Ps(k) 0054 %% Pbus(k) = Bbus(k,k) * Va(k) + Pbusinj(k) 0055 %% Va(k) = inv( Bbus(k,k) ) * (Pbus(k) - Pbusinj(k)) 0056 %% Pf = Bf * Va + Pfinj 0057 %% Pf = Bf(:,k) * Va(k) + Bf(:,ref) * Va(ref) + Pfinj 0058 %% = Bf(:,k) * inv( Bbus(k,k) ) * (Pbus(k) - Pbusinj(k)) + Bf(:,ref) * Va(ref) + Pfinj 0059 %% = W * (Pbus(k) - Pbusinj(k)) + Bf(:,ref) * Va(ref) + Pfinj 0060 %% = W * (Cg(k,:) * Pg - Pd(k) - Ps(k) - Pbusinj(k)) + Bf(:,ref) * Va(ref) + Pfinj 0061 %% = Wg * Pg + Pz; 0062 %% W = Bf(:,k) * inv( Bbus(k,k) ); 0063 %% Wg = W * Cg(k,:); 0064 %% Pz = Bf(:,ref) * Va(ref) + Pfinj - W * (Pd(k) + Ps(k) + Pbusinj(k)); 0065 %% loss = Pf' * R * Pf 0066 %% = (Wg * Pg + Pz)' * R * (Wg * Pg + Pz) 0067 %% = Pg' * Wg' * R * Wg * Pg 0068 %% + 2 * Pz' * R * Wg * Pg 0069 %% + Pz' * R * Pz 0070 0071 r = mpc.branch(:, BR_R); 0072 Varef = mpc.bus(ref, VA) * pi/180; 0073 Pd = mpc.bus(:, PD) / mpc.baseMVA; 0074 Ps = mpc.bus(:, GS) / mpc.baseMVA; 0075 W = Bf(:, k) / Bbus(k, k); %% = Bf(:, k) * inv(Bbus(k, k)) 0076 Wg = W * Cg(k,:); 0077 Pz = Bf(:,ref) * Varef + Pfinj - W * (Pd(k) + Ps(k) + Pbusinj(k)); 0078 R = sparse(1:nl, 1:nl, r, nl, nl); %% nl x nl matrix w/ line resistences on diag 0079 0080 %% form B coefficients 0081 BL2 = 2 * Wg' * R * Wg / mpc.baseMVA; 0082 BL1 = 2 * Wg' * R * Pz; 0083 BL0 = Pz' * R * Pz * mpc.baseMVA;