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] = makeBloss(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] = makeBloss(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 Extras.
0024 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0025 %   See https://github.com/MATPOWER/matpower-extras 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;