Description of t_most_30b_3_1

0001 function t_most_30b_3_1_17(quiet)
0002 %T_MOST_30B_3_1_17  Tests for MOST.
0003 
0004 %   MOST
0005 %   Copyright (c) 2009-2020, 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 https://github.com/MATPOWER/most for more info.
0011 
0012 if nargin < 1
0013     quiet = 0;
0014 end
0015 
0016 n_tests = 96;
0017 
0018 t_begin(n_tests, quiet);
0019 
0020 casename = 't_case30_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     = {'DEFAULT'};  %% opf.dc.solver sequence to try for c3sopf (DC run)
0029 algs.ac     = {'DEFAULT'};  %% 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 if have_feature('gurobi')
0034     mpopt = mpoption(mpopt, 'gurobi.method', 1);    %% dual-simplex
0035 end
0036 if have_feature('mosek')
0037     sc = mosek_symbcon;
0038     mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);     %% dual simplex
0039 end
0040 if have_feature('linprog')
0041     if have_feature('linprog_ds')
0042         mpopt = mpoption(mpopt, 'linprog.Algorithm', 'dual-simplex');
0043     else
0044         mpopt = mpoption(mpopt, 'linprog.Algorithm', 'simplex');
0045     end
0046 end
0047 mpoptac = mpoption(mpopt, 'model', 'AC');
0048 mpoptdc = mpoption(mpopt, 'model', 'DC');
0049 mpopt = mpoption(mpopt, 'most.solver', algs.dc{1});
0050 
0051 %% turn off warnings
0052 s7 = warning('query', 'MATLAB:nearlySingularMatrix');
0053 s6 = warning('query', 'MATLAB:nearlySingularMatrixUMFPACK');
0054 warning('off', 'MATLAB:nearlySingularMatrix');
0055 warning('off', 'MATLAB:nearlySingularMatrixUMFPACK');
0056 
0057 %% define named indices into data matrices
0058 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0059     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0060     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0061 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0062     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0063     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0064 [CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
0065     CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
0066     CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
0067     CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
0068     CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
0069     CT_MODCOST_X] = idx_ct;
0070 
0071 %% reserve and delta offers
0072 xgd_table.colnames = {
0073     'PositiveActiveReservePrice', ...
0074             'PositiveActiveReserveQuantity', ...
0075                     'NegativeActiveReservePrice', ...
0076                             'NegativeActiveReserveQuantity', ...
0077                                     'PositiveActiveDeltaPrice', ...
0078                                             'NegativeActiveDeltaPrice', ...
0079 };
0080 xgd_table.data = [
0081     10.1    15      10.0    15      0.1     0.0;
0082     10.3    30      10.2    30      0.3     0.2;  
0083     10.5    20      10.4    20      0.5     0.4;  
0084     10.7    25      10.6    25      0.7     0.6;
0085     20.1    25      20.0    25      60.1    60.0; 
0086     20.3    15      20.2    15      15.1    15.0;
0087     20.5    30      20.4    30      60.3    60.2;
0088     20.7    15      20.6    15      15.3    15.2;
0089     30.1    15      30.0    15      60.3    60.4;
0090     30.3    30      30.2    30      30.1    30.0;
0091     30.5    25      30.4    25      60.7    60.6;
0092     30.7    30      30.6    30      30.3    30.2;
0093     0.001   50      0.002   50      0       0;
0094     0.001   50      0.002   50      0       0;
0095     0.001   50      0.002   50      0       0;
0096     0.001   50      0.002   50      0       0;
0097     0.001   50      0.002   50      0       0;
0098     0.001   50      0.002   50      0       0;
0099     0.001   50      0.002   50      0       0;
0100     0.001   50      0.002   50      0       0;
0101     0.001   50      0.002   50      0       0;
0102     0.001   50      0.002   50      0       0;
0103     0.001   50      0.002   50      0       0;
0104     0.001   50      0.002   50      0       0;
0105     0.001   50      0.002   50      0       0;
0106     0.001   50      0.002   50      0       0;
0107     0.001   50      0.002   50      0       0;
0108     0.001   50      0.002   50      0       0;
0109     0.001   50      0.002   50      0       0;
0110     0.001   50      0.002   50      0       0;
0111     0.001   50      0.002   50      0       0;
0112     0.001   50      0.002   50      0       0;
0113 ];
0114 
0115 %% contingency table
0116 % label probty  type        row column      chgtype newvalue
0117 contab = [
0118     1   0.002   CT_TBRCH    1   BR_STATUS   CT_REP  0;      %% line 1-2
0119     2   0.002   CT_TBRCH    2   BR_STATUS   CT_REP  0;      %% line 1-3, all power from gen1 flows via gen2
0120     3   0.002   CT_TBRCH    3   BR_STATUS   CT_REP  0;      %% line 2-4, a path to loads @ buses 7 & 8
0121     4   0.002   CT_TBRCH    5   BR_STATUS   CT_REP  0;      %% line 2-5, a path to loads @ buses 7 & 8
0122     5   0.002   CT_TBRCH    6   BR_STATUS   CT_REP  0;      %% line 2-6, a path to loads @ buses 7 & 8
0123     6   0.002   CT_TBRCH    36  BR_STATUS   CT_REP  0;      %% line 28-27, tie line between areas 1 & 3
0124     7   0.002   CT_TBRCH    15  BR_STATUS   CT_REP  0;      %% line 4-12, tie line between areas 1 & 2
0125     8   0.002   CT_TBRCH    12  BR_STATUS   CT_REP  0;      %% line 6-10, tie line between areas 1 & 3
0126     9   0.002   CT_TBRCH    14  BR_STATUS   CT_REP  0;      %% line 9-10, tie line between areas 1 & 3
0127     10  0.002   CT_TGEN     1   GEN_STATUS  CT_REP  0;      %% gen 1 at bus 1
0128     11  0.002   CT_TGEN     2   GEN_STATUS  CT_REP  0;      %% gen 2 at bus 2
0129     12  0.002   CT_TGEN     3   GEN_STATUS  CT_REP  0;      %% gen 3 at bus 22
0130     13  0.002   CT_TGEN     4   GEN_STATUS  CT_REP  0;      %% gen 4 at bus 27
0131     14  0.002   CT_TGEN     5   GEN_STATUS  CT_REP  0;      %% gen 5 at bus 23
0132     15  0.002   CT_TGEN     6   GEN_STATUS  CT_REP  0;      %% gen 6 at bus 13
0133     20  0.010   CT_TGEN     0   PMIN        CT_REL  1.1;    %% 10% load increase
0134     20  0.010   CT_TGEN     0   QMIN        CT_REL  1.1;
0135     21  0.010   CT_TGEN     0   PMIN        CT_REL  0.9;    %% 10% load decrease
0136     21  0.010   CT_TGEN     0   QMIN        CT_REL  0.9;
0137 ];
0138 clist = unique(contab(:, CT_LABEL));
0139 nc = length(clist);
0140 
0141 %% load the case
0142 mpc = loadcase(casename);
0143 gbus = mpc.gen(:, GEN_BUS);
0144 
0145 %%-----  get c3sopf results  -----
0146 rdc = c3sopf_retry(algs.dc, mpc, xgd_table.data, contab, mpoptdc);
0147 % rac = c3sopf_retry(algs.ac, mpc, xgd_table.data, contab, mpoptac);
0148 % save t_most2_soln rdc rac -v6
0149 % s = load('t_most2_soln');
0150 s.rdc = rdc;
0151 % s.rac = rac;
0152 
0153 %%-----  set up data for DC run (most)  -----
0154 ng = size(mpc.gen, 1);      %% number of gens
0155 nt = 3;
0156 xgd = loadxgendata(xgd_table, mpc);
0157 md = loadmd(mpc, nt, xgd, [], contab);
0158 
0159 %%-----  do DC run (most)  -----
0160 r = most(md, mpopt);
0161 
0162 %%-----  test the results  -----
0163 t = 'success1';
0164 t_ok(s.rdc.opf_results.success, t);
0165 t = 'success2';
0166 t_ok(r.QP.exitflag, t);
0167 
0168 t = 'f';
0169 t_is(r.results.f/sum(r.StepProb), s.rdc.opf_results.f, 4, t);
0170 
0171 for tt = 1:nt
0172 
0173 t = sprintf('(t=%d) Pg : base', tt);
0174 t_is(r.flow(tt,1,1).mpc.gen(:, PG), s.rdc.base.gen(:, PG), 5, t);
0175 t = sprintf('(t=%d) Pg : cont ', tt);
0176 for k = 1:nc
0177     t_is(r.flow(tt,1,k+1).mpc.gen(:, PG), s.rdc.cont(k).gen(:, PG), 5, sprintf('%s %d', t, k));
0178 end
0179 
0180 % t = sprintf('(t=%d) gen : base', tt);
0181 % t_is(r.flow(tt,1,1).mpc.gen(:,1:MU_PMIN), s.rdc.base.gen(:,1:MU_PMIN), 3, t);
0182 % t = sprintf('(t=%d) gen : cont ', tt);
0183 % for k = 1:nc
0184 %     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));
0185 % end
0186 
0187 t = sprintf('(t=%d) energy prices', tt);
0188 t_is(r.results.GenPrices(:,tt)/r.StepProb(tt), s.rdc.energy.prc.sum_bus_lam_p(gbus), 6, t);
0189 
0190 t = sprintf('(t=%d) Pc', tt);
0191 t_is(r.results.Pc(:,tt), s.rdc.energy.Pc, 4, t);
0192 
0193 t = sprintf('(t=%d) Gmin', tt);
0194 t_is(r.results.Pc(:,tt) - r.results.Rpm(:,tt), s.rdc.energy.Gmin, 4, t);
0195 
0196 t = sprintf('(t=%d) Gmax', tt);
0197 t_is(r.results.Pc(:,tt) + r.results.Rpp(:,tt), s.rdc.energy.Gmax, 4, t);
0198 
0199 t = sprintf('(t=%d) upward contingency reserve quantities', tt);
0200 t_is(r.results.Rpp(:,tt), s.rdc.reserve.qty.Rp_pos, 4, t);
0201 
0202 t = sprintf('(t=%d) downward contingency reserve quantities', tt);
0203 t_is(r.results.Rpm(:,tt), s.rdc.reserve.qty.Rp_neg, 4, t);
0204 
0205 t = sprintf('(t=%d) upward contingency reserve prices', tt);
0206 t_is(r.results.RppPrices(:,tt)/r.StepProb(tt), s.rdc.reserve.prc.Rp_pos, 6, t);
0207 
0208 t = sprintf('(t=%d) downward contingency reserve prices', tt);
0209 t_is(r.results.RpmPrices(:,tt)/r.StepProb(tt), s.rdc.reserve.prc.Rp_neg, 6, t);
0210 
0211 t = sprintf('(t=%d) contingency physical ramp price', tt);
0212 [vv, ll] = r.om.get_idx();
0213 Ramp_P_max = zeros(ng, nc);
0214 sum_muPmax = zeros(ng, 1);
0215 sum_muPmin = zeros(ng, 1);
0216 for k = 1:nc+1
0217     ii = find(r.flow(tt,1,k).mpc.gen(:, GEN_STATUS) > 0);
0218     if k > 1
0219         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;
0220     end
0221     sum_muPmax(ii) = sum_muPmax(ii) + r.flow(tt,1,k).mpc.gen(ii, MU_PMAX);
0222     sum_muPmin(ii) = sum_muPmin(ii) + r.flow(tt,1,k).mpc.gen(ii, MU_PMIN);
0223 end
0224 t_is(sum(Ramp_P_max/r.StepProb(tt), 2), sum(s.rdc.energy.mu.Ramp_P_max, 2), 6, t);
0225 % t_is(Ramp_P_max/r.StepProb(tt), s.rdc.energy.mu.Ramp_P_max, 6, t);
0226 
0227 t = sprintf('(t=%d) sum_muPmax', tt);
0228 t_is(sum_muPmax/r.StepProb(tt), s.rdc.energy.sum_muPmax, 1, t);
0229 
0230 t = sprintf('(t=%d) sum_muPmin', tt);
0231 t_is(sum_muPmin/r.StepProb(tt), s.rdc.energy.sum_muPmin, 0, t);
0232 
0233 t = sprintf('(t=%d) Rpmax_pos', tt);
0234 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;
0235 t_is(Rpmax_pos, s.rdc.reserve.mu.Rpmax_pos, 6, t);
0236 
0237 t = sprintf('(t=%d) Rpmax_neg', tt);
0238 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;
0239 t_is(Rpmax_neg, s.rdc.reserve.mu.Rpmax_neg, 6, t);
0240 
0241 end
0242 
0243 % g1 = s.rdc.base.gen(:, PG);
0244 % g2 = r.flow(1,1,1).mpc.gen(:, PG);
0245 % for k = 1:nc
0246 %     g1 = [ g1 s.rdc.cont(k).gen(:, PG) ];
0247 %     g2 = [ g2 r.flow(1,1,k+1).mpc.gen(:, PG) ];
0248 % end
0249 % [m,n] = size(g1);
0250 % for j = 1:n
0251 %     fprintf('\n');
0252 %     for i = 1:m
0253 %         fprintf('%9.2f  %9.2f\n', g1(i,j), g2(i,j));
0254 %     end
0255 % end
0256 
0257 %%-----  do AC run (most)  -----
0258 %mostac;
0259 
0260 
0261 %% turn warnings back on
0262 warning(s7.state, 'MATLAB:nearlySingularMatrix');
0263 warning(s6.state, 'MATLAB:nearlySingularMatrixUMFPACK');
0264 
0265 t_end;
t_most_30b_3_1_17

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
