makeBloss

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