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 %   $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;