t_most_3b_3_1_2

T_MOST_3B_1_1_2  Tests for MOST.

0001 function result = t_most_3b_3_1_2(quiet)
0002 %T_MOST_3B_1_1_2  Tests for MOST.
0003 
0004 %   MOST
0005 %   Copyright (c) 2009-2016, Power Systems Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   This file is part of MOST.
0009 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0010 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0011 
0012 if nargin < 1
0013     quiet = 0;
0014 end
0015 
0016 n_tests = 51;
0017 
0018 t_begin(n_tests, quiet);
0019 
0020 casename = 't_case3_most';
0021 fudging = struct( ...       %% paramters for fudging reserve contract for sopf2
0022     'fudge',    0.05, ...   %% initial value (MW)
0023     'step',     0.01, ...   %% if necessary, increase by this amount and retry (MW)
0024     'lim',      0.1);       %% upper limit (MW), give up if no convergence
0025                             %% with fudge equal to this limit
0026 
0027 %% options
0028 algs.dc     = {'MIPS'}; %% opf.dc.solver sequence to try for c3sopf (DC run)
0029 algs.ac     = {'MIPS'}; %% opf.ac.solver sequence to try for c3sopf (AC run)
0030 mpopt = mpoption('verbose', 0, 'out.all', 0);
0031 mpopt = mpoption(mpopt, 'opf.violation', 5e-7, 'mips.comptol', 5e-8);
0032 mpopt = mpoption(mpopt, 'sopf.force_Pc_eq_P0', 0);  %% don't constrain contracted == base case dispatch
0033 mpoptac = mpoption(mpopt, 'model', 'AC');
0034 mpoptdc = mpoption(mpopt, 'model', 'DC');
0035 mpopt = mpoption(mpopt, 'most.solver', algs.dc{1});
0036 
0037 %% turn off warnings
0038 if have_fcn('octave')
0039     s = warning('query', 'Octave:nearly-singular-matrix');
0040     warning('off', 'Octave:nearly-singular-matrix');
0041 else
0042     s7 = warning('query', 'MATLAB:nearlySingularMatrix');
0043     s6 = warning('query', 'MATLAB:nearlySingularMatrixUMFPACK');
0044     warning('off', 'MATLAB:nearlySingularMatrix');
0045     warning('off', 'MATLAB:nearlySingularMatrixUMFPACK');
0046 end
0047 
0048 %% define named indices into data matrices
0049 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0050     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0051     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0052 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0053     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0054     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0055 [CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
0056     CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
0057     CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
0058     CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
0059     CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
0060     CT_MODCOST_X] = idx_ct;
0061 
0062 %% reserve and delta offers
0063 xgd_table.colnames = {
0064     'PositiveActiveReservePrice', ...
0065             'PositiveActiveReserveQuantity', ...
0066                     'NegativeActiveReservePrice', ...
0067                             'NegativeActiveReserveQuantity', ...
0068                                     'PositiveActiveDeltaPrice', ...
0069                                             'NegativeActiveDeltaPrice', ...
0070 };
0071 xgd_table.data = [
0072     1       400     2       400     0.01    0.01;
0073     3       300     4       300     0.01    0.01;
0074     0.001   450     0.002   450     0       0;
0075 ];
0076 
0077 %% contingency table
0078 % label probty  type        row column      chgtype newvalue
0079 contab = [
0080     1   0.01    CT_TBRCH    1   BR_STATUS   CT_REP  0;      %% line 1-2
0081     2   0.01    CT_TGEN     1   GEN_STATUS  CT_REP  0;      %% gen 1 at bus 1
0082 %   2   0.00049621965   CT_TGEN     1   GEN_STATUS  CT_REP  0;      %% gen 1 at bus 1
0083 ];
0084 clist = unique(contab(:, CT_LABEL));
0085 nc = length(clist);
0086 
0087 %% load the case
0088 mpc = loadcase(casename);
0089 gbus = mpc.gen(:, GEN_BUS);
0090 
0091 %%-----  get c3sopf results  -----
0092 rdc = c3sopf_retry(algs.dc, mpc, xgd_table.data, contab, mpoptdc);
0093 % rac = c3sopf_retry(algs.ac, mpc, xgd_table.data, contab, mpoptac);
0094 % save t_most3_soln rdc rac -v6
0095 % s = load('t_most3_soln');
0096 s.rdc = rdc;
0097 % s.rac = rac;
0098 
0099 %%-----  set up data for DC run (most)  -----
0100 ng = size(mpc.gen, 1);      %% number of gens
0101 nt = 3;
0102 xgd = loadxgendata(xgd_table, mpc);
0103 md = loadmd(mpc, nt, xgd, [], contab);
0104 
0105 %%-----  do DC run (most)  -----
0106 r = most(md, mpopt);
0107 
0108 %%-----  test the results  -----
0109 t = 'success1';
0110 t_ok(s.rdc.opf_results.success, t);
0111 t = 'success2';
0112 t_ok(r.QP.exitflag, t);
0113 
0114 t = 'f';
0115 t_is(r.results.f/sum(r.StepProb), s.rdc.opf_results.f, 4, t);
0116 
0117 for tt = 1:nt
0118 
0119 t = sprintf('(t=%d) Pg : base', tt);
0120 t_is(r.flow(tt,1,1).mpc.gen(:, PG), s.rdc.base.gen(:, PG), 5, t);
0121 t = sprintf('(t=%d) Pg : cont ', tt);
0122 for k = 1:nc
0123     t_is(r.flow(tt,1,k+1).mpc.gen(:, PG), s.rdc.cont(k).gen(:, PG), 5, sprintf('%s %d', t, k));
0124 end
0125 
0126 % t = sprintf('(t=%d) gen : base', tt);
0127 % t_is(r.flow(tt,1,1).mpc.gen(:,1:MU_PMIN), s.rdc.base.gen(:,1:MU_PMIN), 3, t);
0128 % t = sprintf('(t=%d) gen : cont ', tt);
0129 % for k = 1:nc
0130 %     t_is(r.flow(tt,1,k+1).mpc.gen(:,1:MU_PMIN), s.rdc.cont(k).gen(:,1:MU_PMIN), 3, sprintf('%s %d', t, k));
0131 % end
0132 
0133 t = sprintf('(t=%d) energy prices', tt);
0134 t_is(r.results.GenPrices(:,tt)/r.StepProb(tt), s.rdc.energy.prc.sum_bus_lam_p(gbus), 6, t);
0135 
0136 t = sprintf('(t=%d) Pc', tt);
0137 t_is(r.results.Pc(:,tt), s.rdc.energy.Pc, 4, t);
0138 
0139 t = sprintf('(t=%d) Gmin', tt);
0140 t_is(r.results.Pc(:,tt) - r.results.Rpm(:,tt), s.rdc.energy.Gmin, 4, t);
0141 
0142 t = sprintf('(t=%d) Gmax', tt);
0143 t_is(r.results.Pc(:,tt) + r.results.Rpp(:,tt), s.rdc.energy.Gmax, 4, t);
0144 
0145 t = sprintf('(t=%d) upward contingency reserve quantities', tt);
0146 t_is(r.results.Rpp(:,tt), s.rdc.reserve.qty.Rp_pos, 4, t);
0147 
0148 t = sprintf('(t=%d) downward contingency reserve quantities', tt);
0149 t_is(r.results.Rpm(:,tt), s.rdc.reserve.qty.Rp_neg, 4, t);
0150 
0151 t = sprintf('(t=%d) upward contingency reserve prices', tt);
0152 t_is(r.results.RppPrices(:,tt)/r.StepProb(tt), s.rdc.reserve.prc.Rp_pos, 6, t);
0153 
0154 t = sprintf('(t=%d) downward contingency reserve prices', tt);
0155 t_is(r.results.RpmPrices(:,tt)/r.StepProb(tt), s.rdc.reserve.prc.Rp_neg, 6, t);
0156 
0157 t = sprintf('(t=%d) contingency physical ramp price', tt);
0158 [vv, ll] = get_idx(r.om);
0159 Ramp_P_max = zeros(ng, nc);
0160 sum_muPmax = zeros(ng, 1);
0161 sum_muPmin = zeros(ng, 1);
0162 for k = 1:nc+1
0163     ii = find(r.flow(tt,1,k).mpc.gen(:, GEN_STATUS) > 0);
0164     if k > 1
0165         Ramp_P_max(ii,k-1) = (r.QP.lambda.mu_u(ll.i1.rampcont(tt,1,k):ll.iN.rampcont(tt,1,k)) - r.QP.lambda.mu_l(ll.i1.rampcont(tt,1,k):ll.iN.rampcont(tt,1,k))) / mpc.baseMVA;
0166     end
0167     sum_muPmax(ii) = sum_muPmax(ii) + r.flow(tt,1,k).mpc.gen(ii, MU_PMAX);
0168     sum_muPmin(ii) = sum_muPmin(ii) + r.flow(tt,1,k).mpc.gen(ii, MU_PMIN);
0169 end
0170 t_is(Ramp_P_max/r.StepProb(tt), s.rdc.energy.mu.Ramp_P_max, 6, t);
0171 
0172 t = sprintf('(t=%d) sum_muPmax', tt);
0173 t_is(sum_muPmax/r.StepProb(tt), s.rdc.energy.sum_muPmax, 2, t);
0174 
0175 t = sprintf('(t=%d) sum_muPmin', tt);
0176 t_is(sum_muPmin/r.StepProb(tt), s.rdc.energy.sum_muPmin, 2, t);
0177 
0178 t = sprintf('(t=%d) Rpmax_pos', tt);
0179 Rpmax_pos = (r.QP.lambda.upper(vv.i1.Rpp(1):vv.iN.Rpp(1)) - r.QP.lambda.lower(vv.i1.Rpp(1):vv.iN.Rpp(1))) / mpc.baseMVA;
0180 t_is(Rpmax_pos, s.rdc.reserve.mu.Rpmax_pos, 6, t);
0181 
0182 t = sprintf('(t=%d) Rpmax_neg', tt);
0183 Rpmax_neg = (r.QP.lambda.upper(vv.i1.Rpm(1):vv.iN.Rpm(1)) - r.QP.lambda.lower(vv.i1.Rpm(1):vv.iN.Rpm(1))) / mpc.baseMVA;
0184 t_is(Rpmax_neg, s.rdc.reserve.mu.Rpmax_neg, 6, t);
0185 
0186 end
0187 
0188 % g1 = s.rdc.base.gen(:, PG);
0189 % g2 = r.flow(1,1,1).mpc.gen(:, PG);
0190 % for k = 1:nc
0191 %     g1 = [ g1 s.rdc.cont(k).gen(:, PG) ];
0192 %     g2 = [ g2 r.flow(1,1,k+1).mpc.gen(:, PG) ];
0193 % end
0194 % [m,n] = size(g1);
0195 % for j = 1:n
0196 %     fprintf('\n');
0197 %     for i = 1:m
0198 %         fprintf('%9.2f  %9.2f\n', g1(i,j), g2(i,j));
0199 %     end
0200 % end
0201 
0202 %%-----  do AC run (most)  -----
0203 %mostac;
0204 
0205 
0206 %% turn warnings back on
0207 if have_fcn('octave')
0208     warning(s.state, 'Octave:nearly-singular-matrix');
0209 else
0210     warning(s7.state, 'MATLAB:nearlySingularMatrix');
0211     warning(s6.state, 'MATLAB:nearlySingularMatrixUMFPACK');
0212 end
0213 
0214 t_end;
0215 
0216 if nargout
0217     result = r;
0218 end