Home > matpower5.0 > t > t_get_losses.m

t_get_losses

PURPOSE ^

T_GET_LOSSES Tests for code in GET_LOSSES.

SYNOPSIS ^

function t_get_losses(quiet)

DESCRIPTION ^

T_GET_LOSSES  Tests for code in GET_LOSSES.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function t_get_losses(quiet)
0002 %T_GET_LOSSES  Tests for code in GET_LOSSES.
0003 
0004 %   MATPOWER
0005 %   $Id: t_get_losses.m 2449 2014-12-05 20:39:29Z ray $
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %   Copyright (c) 2014 by Power System Engineering Research Center (PSERC)
0008 %
0009 %   This file is part of MATPOWER.
0010 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0011 %
0012 %   MATPOWER is free software: you can redistribute it and/or modify
0013 %   it under the terms of the GNU General Public License as published
0014 %   by the Free Software Foundation, either version 3 of the License,
0015 %   or (at your option) any later version.
0016 %
0017 %   MATPOWER is distributed in the hope that it will be useful,
0018 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0019 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0020 %   GNU General Public License for more details.
0021 %
0022 %   You should have received a copy of the GNU General Public License
0023 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0024 %
0025 %   Additional permission under GNU GPL version 3 section 7
0026 %
0027 %   If you modify MATPOWER, or any covered work, to interface with
0028 %   other modules (such as MATLAB code and MEX-files) available in a
0029 %   MATLAB(R) or comparable environment containing parts covered
0030 %   under other licensing terms, the licensors of MATPOWER grant
0031 %   you additional permission to convey the resulting work.
0032 
0033 if nargin < 1
0034     quiet = 0;
0035 end
0036 
0037 n_tests = 20;
0038 
0039 t_begin(n_tests, quiet);
0040 
0041 %% define named indices into data matrices
0042 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0043     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0044 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0045     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0046     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0047 
0048 casefile = 't_case9_opf';
0049 if quiet
0050     verbose = 0;
0051 else
0052     verbose = 0;
0053 end
0054 if have_fcn('octave')
0055     s1 = warning('query', 'Octave:load-file-in-path');
0056     warning('off', 'Octave:load-file-in-path');
0057 end
0058 
0059 mpopt = mpoption('opf.violation', 1e-6, 'mips.gradtol', 1e-8, ...
0060         'mips.comptol', 1e-8, 'mips.costtol', 1e-9);
0061 mpopt = mpoption(mpopt, 'out.all', 0, 'verbose', verbose, 'opf.ac.solver', 'MIPS');
0062 mpc = loadcase(casefile);
0063 mpc.branch(8, TAP) = 0.95;
0064 mpc.branch(4, SHIFT) = -2;
0065 r = runopf(mpc, mpopt);
0066 p_loss = [0; 0.250456; 0.944337; 0; 0.151545; 0.298009; 0; 1.600081; 0.216618];
0067 q_loss = [3.890347; 1.355409; 4.116342; 5.017914; 1.283672; 2.524309; 10.20216; 8.050407; 1.841254];
0068 the_fchg = [0; 9.476543; 20.751015; 0; 12.113509; 8.563392; 0; 20.056039; 10.431575];
0069 the_tchg = [0; 9.158269; 20.749455; 0; 12.011738; 8.813679; 0; 18.136716; 10.556150];
0070 
0071 %%-----  all load  -----
0072 t = 'loss = get_losses(results) : ';
0073 loss = get_losses(r);
0074 t_is(real(loss), p_loss, 6, [t 'P loss']);
0075 t_is(imag(loss), q_loss, 6, [t 'Q loss']);
0076 
0077 t = 'loss = get_losses(baseMVA, bus, branch) : ';
0078 loss = get_losses(r.baseMVA, r.bus, r.branch);
0079 t_is(real(loss), p_loss, 6, [t 'P loss']);
0080 t_is(imag(loss), q_loss, 6, [t 'Q loss']);
0081 
0082 t = '[loss, chg, chg] = get_losses(results) : ';
0083 [loss, chg] = get_losses(r);
0084 t_is(real(loss), p_loss, 6, [t 'P loss']);
0085 t_is(imag(loss), q_loss, 6, [t 'Q loss']);
0086 t_is(chg, the_fchg+the_tchg, 6, [t 'Q inj (total)']);
0087 
0088 t = '[loss, fchg, tchg] = get_losses(results) : ';
0089 [loss, fchg, tchg] = get_losses(r);
0090 t_is(real(loss), p_loss, 6, [t 'P loss']);
0091 t_is(imag(loss), q_loss, 6, [t 'Q loss']);
0092 t_is(fchg, the_fchg, 6, [t 'Q inj (from)']);
0093 t_is(tchg, the_tchg, 6, [t 'Q inj (to)']);
0094 t_is(real(loss), r.branch(:, PF)+r.branch(:, PT), 12, [t 'P_loss = Pf+Pt']);
0095 t_is(imag(loss), r.branch(:, QF)+r.branch(:, QT)+tchg+fchg, 12, [t 'Q_loss = Qf+Qt+fchg+tchg']);
0096 
0097 
0098 t = '[loss, fchg, tchg, dloss_dV, dchg_dVm] = get_losses(r) : ';
0099 [loss, fchg, tchg, dloss_dV, dchg_dVm] = get_losses(r);
0100 epsilon = 1e-8;
0101 %% build numerical versions of dloss_dVa and dloss_dVm
0102 nb = size(r.bus, 1);    %% number of buses
0103 nl = size(r.branch, 1); %% number of branches
0104 dloss_dVa = sparse(nl, nb);
0105 dloss_dVm = sparse(nl, nb);
0106 dfchg_dVm = sparse(nl, nb);
0107 dtchg_dVm = sparse(nl, nb);
0108 for j = 1:nb
0109     mpc = r;
0110     mpc.bus(j, VA) = mpc.bus(j, VA) + 180 / pi * epsilon;
0111     loss2 = get_losses(mpc);
0112     dloss_dVa(:, j) = (loss2 - loss) / epsilon;
0113 end
0114 for j = 1:nb
0115     mpc = r;
0116     mpc.bus(j, VM) = mpc.bus(j, VM) + epsilon;
0117     [loss2, fchg2, tchg2] = get_losses(mpc);
0118     dloss_dVm(:, j) = (loss2 - loss) / epsilon;
0119     dfchg_dVm(:, j) = (fchg2 - fchg) / epsilon;
0120     dtchg_dVm(:, j) = (tchg2 - tchg) / epsilon;
0121 end
0122 t_is(full(real(dloss_dV.a)), full(real(dloss_dVa)), 5, [t 'dPloss/dVa']);
0123 t_is(full(imag(dloss_dV.a)), full(imag(dloss_dVa)), 4, [t 'dQloss/dVa']);
0124 t_is(full(real(dloss_dV.m)), full(real(dloss_dVm)), 5, [t 'dPloss/dVm']);
0125 t_is(full(imag(dloss_dV.m)), full(imag(dloss_dVm)), 4, [t 'dQloss/dVm']);
0126 t_is(full(dchg_dVm.f), full(dfchg_dVm), 5, [t 'dfchg/dVm']);
0127 t_is(full(dchg_dVm.t), full(dtchg_dVm), 5, [t 'dtchg/dVm']);
0128 
0129 t = 'Loss Sensitivity Factors (LSF)';
0130 %% convert to internal indexing to use makeJac()
0131 ri = ext2int(r);            %% results with internal indexing
0132 [loss, fchg, tchg, dloss_dV] = get_losses(ri);
0133 nb = size(ri.bus, 1);
0134 nl = size(ri.branch, 1);
0135 [ref, pv, pq] = bustypes(ri.bus, ri.gen);
0136 J = makeJac(ri);
0137 dL = real([dloss_dV.a(:, [pv;pq]) dloss_dV.m(:, pq)]) / ri.baseMVA;
0138 LSFi = zeros(nl, 2*nb);
0139 LSFi(:, [pv; pq; nb+pq]) = dL / J;  %% loss sensitivity factors in internal indexing
0140 
0141 %% convert to external indexing
0142 nb = size(r.bus, 1);
0143 nl = size(r.branch, 1);
0144 LSF = zeros(nl, 2*nb);
0145 LSF(ri.order.branch.status.on, [ri.order.bus.status.on; nb+ri.order.bus.status.on]) = LSFi;
0146 
0147 numLSF = zeros(nl, 2*nb);
0148 epsilon = 1e-4;
0149 mpopt = mpoption('out.all', 0, 'verbose', 0, 'pf.tol', 1e-12);
0150 for j = 1:nb
0151     mpc = r;
0152     mpc.bus(j, PD) = mpc.bus(j, PD) - epsilon;
0153     rr = runpf(mpc, mpopt);
0154     loss2 = get_losses(rr);
0155     numLSF(:, j) = real(loss2 - loss) / epsilon;
0156     
0157     mpc = r;
0158     mpc.bus(j, QD) = mpc.bus(j, QD) - epsilon;
0159     rr = runpf(mpc, mpopt);
0160     loss2 = get_losses(rr);
0161     numLSF(:, nb+j) = real(loss2 - loss) / epsilon;
0162 end
0163 t_is(LSF, numLSF, 1e-6, t);
0164 
0165 % norm(numLSF - LSF, Inf)
0166 % norm(numLSF(:, 1:nb) - LSF(:, 1:nb), Inf)
0167 
0168 % LSF
0169 % numLSF
0170 
0171 % sum(LSF)
0172 % sum(numLSF)
0173 
0174 % norm(sum(numLSF) - sum(LSF), Inf)
0175 
0176 t_end;
Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005