0001 function t_get_losses(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 if nargin < 1
0013 quiet = 0;
0014 end
0015
0016 n_tests = 20;
0017
0018 t_begin(n_tests, quiet);
0019
0020
0021 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0022 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0023 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0024 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0025 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0026
0027 casefile = 't_case9_opf';
0028 if quiet
0029 verbose = 0;
0030 else
0031 verbose = 0;
0032 end
0033 if have_feature('octave')
0034 if have_feature('octave', 'vnum') >= 4
0035 file_in_path_warn_id = 'Octave:data-file-in-path';
0036 else
0037 file_in_path_warn_id = 'Octave:load-file-in-path';
0038 end
0039 s1 = warning('query', file_in_path_warn_id);
0040 warning('off', file_in_path_warn_id);
0041 end
0042
0043 mpopt = mpoption('opf.violation', 1e-6, 'mips.gradtol', 1e-8, ...
0044 'mips.comptol', 1e-8, 'mips.costtol', 1e-9);
0045 mpopt = mpoption(mpopt, 'out.all', 0, 'verbose', verbose, 'opf.ac.solver', 'MIPS');
0046 mpc = loadcase(casefile);
0047 mpc.branch(8, TAP) = 0.95;
0048 mpc.branch(4, SHIFT) = -2;
0049 r = runopf(mpc, mpopt);
0050 p_loss = [0; 0.250456; 0.944337; 0; 0.151545; 0.298009; 0; 1.600081; 0.216618];
0051 q_loss = [3.890347; 1.355409; 4.116342; 5.017914; 1.283672; 2.524309; 10.20216; 8.050407; 1.841254];
0052 the_fchg = [0; 9.476543; 20.751015; 0; 12.113509; 8.563392; 0; 20.056039; 10.431575];
0053 the_tchg = [0; 9.158269; 20.749455; 0; 12.011738; 8.813679; 0; 18.136716; 10.556150];
0054
0055
0056 t = 'loss = get_losses(results) : ';
0057 loss = get_losses(r);
0058 t_is(real(loss), p_loss, 6, [t 'P loss']);
0059 t_is(imag(loss), q_loss, 6, [t 'Q loss']);
0060
0061 t = 'loss = get_losses(baseMVA, bus, branch) : ';
0062 loss = get_losses(r.baseMVA, r.bus, r.branch);
0063 t_is(real(loss), p_loss, 6, [t 'P loss']);
0064 t_is(imag(loss), q_loss, 6, [t 'Q loss']);
0065
0066 t = '[loss, chg, chg] = get_losses(results) : ';
0067 [loss, chg] = get_losses(r);
0068 t_is(real(loss), p_loss, 6, [t 'P loss']);
0069 t_is(imag(loss), q_loss, 6, [t 'Q loss']);
0070 t_is(chg, the_fchg+the_tchg, 6, [t 'Q inj (total)']);
0071
0072 t = '[loss, fchg, tchg] = get_losses(results) : ';
0073 [loss, fchg, tchg] = 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 t_is(fchg, the_fchg, 6, [t 'Q inj (from)']);
0077 t_is(tchg, the_tchg, 6, [t 'Q inj (to)']);
0078 t_is(real(loss), r.branch(:, PF)+r.branch(:, PT), 12, [t 'P_loss = Pf+Pt']);
0079 t_is(imag(loss), r.branch(:, QF)+r.branch(:, QT)+tchg+fchg, 12, [t 'Q_loss = Qf+Qt+fchg+tchg']);
0080
0081
0082 t = '[loss, fchg, tchg, dloss_dV, dchg_dVm] = get_losses(r) : ';
0083 [loss, fchg, tchg, dloss_dV, dchg_dVm] = get_losses(r);
0084 epsilonVa = 1e-9;
0085 epsilonVm = 1e-8;
0086
0087 nb = size(r.bus, 1);
0088 nl = size(r.branch, 1);
0089 dloss_dVa = sparse(nl, nb);
0090 dloss_dVm = sparse(nl, nb);
0091 dfchg_dVm = sparse(nl, nb);
0092 dtchg_dVm = sparse(nl, nb);
0093 for j = 1:nb
0094 mpc = r;
0095 mpc.bus(j, VA) = mpc.bus(j, VA) + 180 / pi * epsilonVa;
0096 loss2 = get_losses(mpc);
0097 dloss_dVa(:, j) = (loss2 - loss) / epsilonVa;
0098 end
0099 for j = 1:nb
0100 mpc = r;
0101 mpc.bus(j, VM) = mpc.bus(j, VM) + epsilonVm;
0102 [loss2, fchg2, tchg2] = get_losses(mpc);
0103 dloss_dVm(:, j) = (loss2 - loss) / epsilonVm;
0104 dfchg_dVm(:, j) = (fchg2 - fchg) / epsilonVm;
0105 dtchg_dVm(:, j) = (tchg2 - tchg) / epsilonVm;
0106 end
0107 t_is(full(real(dloss_dV.a)), full(real(dloss_dVa)), 5, [t 'dPloss/dVa']);
0108 t_is(full(imag(dloss_dV.a)), full(imag(dloss_dVa)), 4, [t 'dQloss/dVa']);
0109 t_is(full(real(dloss_dV.m)), full(real(dloss_dVm)), 5, [t 'dPloss/dVm']);
0110 t_is(full(imag(dloss_dV.m)), full(imag(dloss_dVm)), 4, [t 'dQloss/dVm']);
0111 t_is(full(dchg_dVm.f), full(dfchg_dVm), 5, [t 'dfchg/dVm']);
0112 t_is(full(dchg_dVm.t), full(dtchg_dVm), 6, [t 'dtchg/dVm']);
0113
0114 t = 'Loss Sensitivity Factors (LSF)';
0115
0116 ri = ext2int(r);
0117 [loss, fchg, tchg, dloss_dV] = get_losses(ri);
0118 nb = size(ri.bus, 1);
0119 nl = size(ri.branch, 1);
0120 [ref, pv, pq] = bustypes(ri.bus, ri.gen);
0121 J = makeJac(ri);
0122 dL = real([dloss_dV.a(:, [pv;pq]) dloss_dV.m(:, pq)]) / ri.baseMVA;
0123 LSFi = zeros(nl, 2*nb);
0124 LSFi(:, [pv; pq; nb+pq]) = dL / J;
0125
0126
0127 nb = size(r.bus, 1);
0128 nl = size(r.branch, 1);
0129 LSF = zeros(nl, 2*nb);
0130 LSF(ri.order.branch.status.on, [ri.order.bus.status.on; nb+ri.order.bus.status.on]) = LSFi;
0131
0132 numLSF = zeros(nl, 2*nb);
0133 epsilon = 5e-6;
0134 mpopt = mpoption('out.all', 0, 'verbose', 0, 'pf.tol', 1e-12);
0135 for j = 1:nb
0136 mpc = r;
0137 mpc.bus(j, PD) = mpc.bus(j, PD) - epsilon;
0138 rr = runpf(mpc, mpopt);
0139 loss2 = get_losses(rr);
0140 numLSF(:, j) = real(loss2 - loss) / epsilon;
0141
0142 mpc = r;
0143 mpc.bus(j, QD) = mpc.bus(j, QD) - epsilon;
0144 rr = runpf(mpc, mpopt);
0145 loss2 = get_losses(rr);
0146 numLSF(:, nb+j) = real(loss2 - loss) / epsilon;
0147 end
0148 t_is(LSF, numLSF, 6, t);
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161 t_end;