0001 function t_get_losses(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033 if nargin < 1
0034 quiet = 0;
0035 end
0036
0037 n_tests = 20;
0038
0039 t_begin(n_tests, quiet);
0040
0041
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
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
0102 nb = size(r.bus, 1);
0103 nl = size(r.branch, 1);
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
0131 ri = ext2int(r);
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;
0140
0141
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
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176 t_end;