t_most_suc

T_MOST_SUC  Tests of stochastic unit commitment optimizations.

   T_MOST_SUC(QUIET, CREATE_PLOTS, CREATE_PDFS, SAVEDIR)
   Can generate summary plots and save them as PDFs in a directory of
   your choice.
   E.g. t_most_suc(0, 1, 1, '~/Downloads/suc_plots')

0001 function t_most_suc(quiet, create_plots, create_pdfs, savedir)
0002 %T_MOST_SUC  Tests of stochastic unit commitment optimizations.
0003 %
0004 %   T_MOST_SUC(QUIET, CREATE_PLOTS, CREATE_PDFS, SAVEDIR)
0005 %   Can generate summary plots and save them as PDFs in a directory of
0006 %   your choice.
0007 %   E.g. t_most_suc(0, 1, 1, '~/Downloads/suc_plots')
0008 
0009 %   MOST
0010 %   Copyright (c) 2015-2020, Power Systems Engineering Research Center (PSERC)
0011 %   by Ray Zimmerman, PSERC Cornell
0012 %
0013 %   This file is part of MOST.
0014 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0015 %   See https://github.com/MATPOWER/most for more info.
0016 
0017 if nargin < 4
0018     savedir = '.';              %% save in current working directory by default
0019     if nargin < 3
0020         create_pdfs = 0;        %% do NOT save plots to PDF files
0021         if nargin < 2
0022             create_plots = 0;   %% do NOT create summary plots of results
0023             if nargin < 1
0024                 quiet = 0;      %% verbose by default
0025             end
0026         end
0027     end
0028 end
0029 if create_plots
0030     if create_pdfs
0031         fname = 'suc-ex';
0032     else
0033         fname = '';
0034     end
0035     pp = 0;     %% plot counter
0036 end
0037 
0038 solvers = {'CPLEX', 'GLPK', 'GUROBI', 'MOSEK', 'OT'};
0039 fcn = {'cplex', 'glpk', 'gurobi', 'mosek', 'intlinprog'};
0040 % solvers = {'CPLEX'};
0041 % fcn = {'cplex'};
0042 % solvers = {'OT'};
0043 % fcn = {'intlinprog'};
0044 % solvers = {'GUROBI'};
0045 % fcn = {'gurobi'};
0046 % solvers = {'DEFAULT'};
0047 % fcn = {'most'};
0048 ntests = 37;
0049 t_begin(ntests*length(solvers), quiet);
0050 
0051 if quiet
0052     verbose = 0;
0053 else
0054     verbose = 0;
0055 end
0056 % verbose = 2;
0057 
0058 if have_feature('octave')
0059     if have_feature('octave', 'vnum') >= 4
0060         file_in_path_warn_id = 'Octave:data-file-in-path';
0061     else
0062         file_in_path_warn_id = 'Octave:load-file-in-path';
0063     end
0064     s1 = warning('query', file_in_path_warn_id);
0065     warning('off', file_in_path_warn_id);
0066 end
0067 
0068 casefile = 'ex_case3b';
0069 solnfile =  't_most_suc_soln';
0070 soln = load(solnfile);
0071 mpopt = mpoption;
0072 mpopt = mpoption(mpopt, 'out.gen', 1);
0073 mpopt = mpoption(mpopt, 'verbose', verbose);
0074 % mpopt = mpoption(mpopt, 'opf.violation', 1e-6, 'mips.gradtol', 1e-8, ...
0075 %         'mips.comptol', 1e-8, 'mips.costtol', 1e-8);
0076 mpopt = mpoption(mpopt, 'model', 'DC');
0077 mpopt = mpoption(mpopt, 'most.price_stage_warn_tol', 10);
0078 
0079 %% solver options
0080 if have_feature('cplex')
0081     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 0);       %% automatic
0082     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 1);       %% primal simplex
0083     mpopt = mpoption(mpopt, 'cplex.lpmethod', 2);       %% dual simplex
0084     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 3);       %% network simplex
0085     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 4);       %% barrier
0086     mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.mipgap', 0);
0087     mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.absmipgap', 0);
0088     mpopt = mpoption(mpopt, 'cplex.opts.threads', 2);
0089 end
0090 if have_feature('glpk')
0091     mpopt = mpoption(mpopt, 'glpk.opts.mipgap', 0);
0092     mpopt = mpoption(mpopt, 'glpk.opts.tolint', 1e-10);
0093     mpopt = mpoption(mpopt, 'glpk.opts.tolobj', 1e-10);
0094 end
0095 if have_feature('gurobi')
0096     %mpopt = mpoption(mpopt, 'gurobi.method', -1);       %% automatic
0097     %mpopt = mpoption(mpopt, 'gurobi.method', 0);        %% primal simplex
0098     mpopt = mpoption(mpopt, 'gurobi.method', 1);        %% dual simplex
0099     %mpopt = mpoption(mpopt, 'gurobi.method', 2);        %% barrier
0100     mpopt = mpoption(mpopt, 'gurobi.threads', 2);
0101     mpopt = mpoption(mpopt, 'gurobi.opts.MIPGap', 0);
0102     mpopt = mpoption(mpopt, 'gurobi.opts.MIPGapAbs', 0);
0103 end
0104 if have_feature('mosek')
0105     sc = mosek_symbcon;
0106     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_FREE);            %% default
0107     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_INTPNT);          %% interior point
0108     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_PRIMAL_SIMPLEX);  %% primal simplex
0109     mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);     %% dual simplex
0110     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_FREE_SIMPLEX);    %% automatic simplex
0111     %mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_X', 0);
0112     mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_NODE_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
0113     mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_ROOT_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
0114     mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_RELAX_INT', 1e-9);
0115     %mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_RELAX_INT', 0);
0116     mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_GAP', 0);
0117     mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_GAP', 0);
0118 end
0119 if have_feature('intlinprog')
0120     %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'interior-point');
0121     %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'active-set');
0122     %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'simplex');
0123     mpopt = mpoption(mpopt, 'linprog.Algorithm', 'dual-simplex');
0124     %mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'primal-simplex');
0125     mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'dual-simplex');
0126     mpopt = mpoption(mpopt, 'intlinprog.TolCon', 1e-9);
0127     mpopt = mpoption(mpopt, 'intlinprog.TolGapAbs', 0);
0128     mpopt = mpoption(mpopt, 'intlinprog.TolGapRel', 0);
0129     mpopt = mpoption(mpopt, 'intlinprog.TolInteger', 1e-6);
0130     %% next line is to work around a bug in intlinprog
0131     % (Technical Support Case #01841662)
0132     mpopt = mpoption(mpopt, 'intlinprog.LPPreprocess', 'none');
0133 end
0134 if ~verbose
0135     mpopt = mpoption(mpopt, 'out.all', 0);
0136 end
0137 % mpopt = mpoption(mpopt, 'out.all', -1);
0138 
0139 %% define named indices into data matrices
0140 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0141     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0142 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0143     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0144     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0145 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0146     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0147     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0148 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0149 [CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
0150     CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
0151     CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
0152     CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
0153     CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
0154     CT_MODCOST_X] = idx_ct;
0155 
0156 %% load base case file
0157 mpc = loadcase(casefile);
0158 
0159 nb = size(mpc.bus, 1);
0160 nl = size(mpc.branch, 1);
0161 ng = size(mpc.gen, 1);
0162 
0163 xgd = loadxgendata('ex_xgd_uc', mpc);
0164 [iwind, mpc, xgd] = addwind('ex_wind_uc', mpc, xgd);
0165 profiles_d = getprofiles('ex_wind_profile_d', iwind);
0166 profiles_d = getprofiles('ex_load_profile', profiles_d);
0167 profiles_s = getprofiles('ex_wind_profile', iwind);
0168 profiles_s = getprofiles('ex_load_profile', profiles_s);
0169 nt = size(profiles_d(1).values, 1);
0170 
0171 mpc0 = mpc;
0172 xgd0 = xgd;
0173 
0174 for s = 1:length(solvers)
0175     if ~have_feature(fcn{s})    %% check if we have the solver
0176         t_skip(ntests, sprintf('%s not installed', solvers{s}));
0177     else
0178         mpopt = mpoption(mpopt, 'opf.dc.solver', solvers{s});
0179         mpopt = mpoption(mpopt, 'most.solver', mpopt.opf.dc.solver);
0180         mpopt = mpoption(mpopt, 'most.storage.cyclic', 1);
0181 
0182         mpc = mpc0;
0183         xgd = xgd0;
0184 
0185         t = sprintf('%s : deterministic : ', solvers{s});
0186         mdi = loadmd(mpc, nt, xgd, [], [], profiles_d);
0187         mdo = most(mdi, mpopt);
0188         ms = most_summary(mdo);
0189         t_ok(mdo.QP.exitflag > 0, [t 'success']);
0190         ex = soln.determ;
0191         t_is(ms.f, ex.f, 8, [t 'f']);
0192         t_is(ms.Pg, ex.Pg, 8, [t 'Pg']);
0193         t_is(ms.Rup, ex.Rup, 8, [t 'Rup']);
0194         t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
0195         t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
0196         t_is(ms.u, ex.u, 8, [t 'u']);
0197         t_is(ms.lamP, ex.lamP, 8, [t 'lamP']);
0198         t_is(ms.muF, ex.muF, 8, [t 'muF']);
0199         % determ = most_summary(mdo);
0200         if create_plots
0201             pp = pp + 1;
0202             plot_case('Base : Deterministic', mdo, ms, 500, 150, savedir, pp, fname);
0203         end
0204         % keyboard;
0205 
0206         t = sprintf('%s : individual trajectories : ', solvers{s});
0207         transmat_s = cell(1, nt);
0208         I = speye(3);
0209         [transmat_s{:}] = deal(I);
0210         transmat_s{1} = [ 0.158655253931457; 0.682689492137086; 0.158655253931457 ];
0211         mdi = loadmd(mpc, transmat_s, xgd, [], [], profiles_s);
0212         mdi = filter_ramp_transitions(mdi, 0.1);
0213         mdo = most(mdi, mpopt);
0214         ms = most_summary(mdo);
0215         t_ok(mdo.QP.exitflag > 0, [t 'success']);
0216         ex = soln.transprob1;
0217         t_is(ms.f, ex.f, 5, [t 'f']);
0218         t_is(ms.Pg, ex.Pg, 6, [t 'Pg']);
0219         t_is(ms.Rup, ex.Rup, 6, [t 'Rup']);
0220         t_is(ms.Rdn, ex.Rdn, 6, [t 'Rdn']);
0221         t_is(ms.Pf, ex.Pf, 6, [t 'Pf']);
0222         t_is(ms.u, ex.u, 8, [t 'u']);
0223         % t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
0224         % t_is(ms.muF, ex.muF, 5, [t 'muF']);
0225         % transprob1 = most_summary(mdo);
0226         if create_plots
0227             pp = pp + 1;
0228             plot_case('Individual Trajectories', mdo, ms, 500, 150, savedir, pp, fname);
0229         end
0230         % keyboard;
0231 
0232         t = sprintf('%s : full transition probabilities : ', solvers{s});
0233 %        transmat_sf = transmat_s;
0234         transmat_sf = ex_transmat(nt);
0235 %         transmat_sf = cell(1, nt);
0236 %         T = [ 0.158655253931457; 0.682689492137086; 0.158655253931457 ];
0237 %         [transmat_sf{:}] = deal(T * ones(1,3));
0238 %         transmat_sf{1} = T;
0239         mdi = loadmd(mpc, transmat_sf, xgd, [], [], profiles_s);
0240 %        mdi = filter_ramp_transitions(mdi, 0.9);
0241         mdo = most(mdi, mpopt);
0242         ms = most_summary(mdo);
0243         t_ok(mdo.QP.exitflag > 0, [t 'success']);
0244         ex = soln.transprobfull;
0245         t_is(ms.f, ex.f, 3, [t 'f']);
0246         t_is(ms.Pg, ex.Pg, 3, [t 'Pg']);
0247         t_is(ms.Rup, ex.Rup, 3, [t 'Rup']);
0248         t_is(ms.Rdn, ex.Rdn, 3, [t 'Rdn']);
0249         t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
0250         t_is(ms.u, ex.u, 8, [t 'u']);
0251         % t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
0252         % t_is(ms.muF, ex.muF, 5, [t 'muF']);
0253         % transprobfull = most_summary(mdo);
0254         if create_plots
0255             pp = pp + 1;
0256             plot_case('Full Transition Probabilities', mdo, ms, 500, 150, savedir, pp, fname);
0257         end
0258         % keyboard;
0259 
0260         t = sprintf('%s : full transition probabilities + cont : ', solvers{s});
0261         mdi = loadmd(mpc, transmat_sf, xgd, [], 'ex_contab', profiles_s);
0262 %        mdi = filter_ramp_transitions(mdi, 0.9);
0263         mdo = most(mdi, mpopt);
0264         ms = most_summary(mdo);
0265         t_ok(mdo.QP.exitflag > 0, [t 'success']);
0266         ex = soln.transprobcont;
0267         t_is(ms.f, ex.f, 4, [t 'f']);
0268         t_is(ms.Pg, ex.Pg, 6, [t 'Pg']);
0269         t_is(ms.Rup, ex.Rup, 6, [t 'Rup']);
0270         t_is(ms.Rdn, ex.Rdn, 6, [t 'Rdn']);
0271         t_is(ms.Pf, ex.Pf, 6, [t 'Pf']);
0272         t_is(ms.u, ex.u, 8, [t 'u']);
0273         % t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
0274         % t_is(ms.muF, ex.muF, 5, [t 'muF']);
0275         % transprobcont = most_summary(mdo);
0276         if create_plots
0277             pp = pp + 1;
0278             plot_case('+ Contingencies', mdo, ms, 500, 150, savedir, pp, fname);
0279         end
0280         % keyboard;
0281 
0282         t = sprintf('%s : + storage : ', solvers{s});
0283         if mpopt.out.all
0284             fprintf('Add storage\n');
0285         end
0286         [iess, mpc, xgd, sd] = addstorage('ex_storage', mpc, xgd);
0287         mdi = loadmd(mpc, transmat_sf, xgd, sd, 'ex_contab', profiles_s);
0288         mdo = most(mdi, mpopt);
0289         ms = most_summary(mdo);
0290         t_ok(mdo.QP.exitflag > 0, [t 'success']);
0291         ex = soln.wstorage;
0292         t_is(ms.f, ex.f, 3, [t 'f']);
0293         t_is(ms.Pg, ex.Pg, 3, [t 'Pg']);
0294         t_is(ms.Rup, ex.Rup, 3, [t 'Rup']);
0295         t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
0296         t_is(ms.Pf, ex.Pf, 3, [t 'Pf']);
0297         t_is(ms.u, ex.u, 8, [t 'u']);
0298         % t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
0299         % t_is(ms.muF, ex.muF, 5, [t 'muF']);
0300         % wstorage = most_summary(mdo);
0301         if create_plots
0302             pp = pp + 1;
0303             plot_case('+ Storage', mdo, ms, 500, 150, savedir, pp, fname);
0304             create_plots = 0;   %% don't do them again
0305         end
0306         % keyboard;
0307     end
0308 end
0309 
0310 if have_feature('octave')
0311     warning(s1.state, file_in_path_warn_id);
0312 end
0313 
0314 t_end;
0315 
0316 % save t_most_suc_soln determ transprob1 transprobfull transprobcont wstorage
0317 % determ.u
0318 % transprob1.u
0319 % transprobcont.u
0320 % transprobfull.u
0321 % wstorage.u
0322 
0323 
0324 function h = plot_case(label, md, ms, maxq, maxp, mypath, pp, fname)
0325 
0326 if nargin < 8
0327     fname = '';
0328 end
0329 
0330 %% colors:  blue     red               yellow           purple            green
0331 cc = {[0 0.45 0.74], [0.85 0.33 0.1], [0.93 0.69 0.13], [0.49 0.18 0.56], [0.47 0.67 0.19]};
0332 
0333 ig = (1:3)';
0334 id = 4;
0335 iw = 5;
0336 is = 6;
0337 
0338 subplot(3, 1, 1);
0339 md.mpc = rmfield(md.mpc, 'genfuel');
0340 plot_uc(md, [], 'title', label);
0341 ylabel('Unit Commitment', 'FontSize', 16);
0342 ah = gca;
0343 ah.YAxisLocation = 'left';
0344 
0345 subplot(3, 1, 2);
0346 x = (1:ms.nt)';
0347 Pg = md.results.ExpectedDispatch;
0348 y1 = Pg(ig, :)';
0349 if ms.ng == 6
0350     y1 = [y1 max(-Pg(is, :), 0)' max(Pg(is, :), 0)'];
0351 end
0352 y2 = -sum(Pg([id; iw], :), 1)';
0353 [ah1, h1, h2] = plotyy(x, y1, x, y2);
0354 axis(ah1(1), [0.5 12.5 0 maxq]);
0355 axis(ah1(2), [0.5 12.5 0 maxq]);
0356 % ah1(1).XLim = [0.5 12.5];
0357 % ah1(2).XLim = [0.5 12.5];
0358 % ah1(1).YLim = [0 300];
0359 % ah1(2).YLim = [0 450];
0360 ah1(1).YTickMode = 'auto';
0361 ah1(2).YTickMode = 'auto';
0362 ah1(1).XTick = 1:12;
0363 nn = 3;
0364 for j = 1:3
0365     h1(j).LineWidth = 2;
0366     h1(j).Color = cc{j};
0367 end
0368 if ms.ng == 6
0369     h1(4).LineWidth = 2;
0370     h1(4).Color = cc{5};
0371     h1(4).LineStyle = ':';
0372     h1(5).LineWidth = 2;
0373     h1(5).Color = cc{5};
0374 end
0375 h2.LineWidth = 2;
0376 h2.Color = cc{4};
0377 h2.LineStyle = ':';
0378 ah1(2).YColor = cc{4};
0379 %title('Generation & Net Load', 'FontSize', 16);
0380 ylabel(ah1(1), 'Generation, MW', 'FontSize', 16);
0381 ylabel(ah1(2), 'Net Load, MW', 'FontSize', 16);
0382 xlabel('Period', 'FontSize', 16);
0383 set(ah1(1), 'FontSize', 14);
0384 set(ah1(2), 'FontSize', 14);
0385 if ms.ng == 6
0386     legend('Gen 1', 'Gen 2', 'Gen 3', 'Storage Charge', 'Storage Discharge', 'Location', [0.7 0.6 0 0]);
0387 else
0388     legend('Gen 1', 'Gen 2', 'Gen 3', 'Location', [0.7 0.58 0 0]);
0389 end
0390 
0391 subplot(3, 1, 3);
0392 if length(size(ms.lamP)) == 4
0393     elamP = sum(sum(ms.lamP, 4), 3) ./ (ones(ms.nb,1) * md.StepProb);
0394 %     elamP = sum(sum(ms.lamP, 4), 3);
0395 elseif length(size(ms.lamP)) == 3
0396     elamP = sum(ms.lamP, 3);
0397 else
0398     elamP = ms.lamP;
0399 end
0400 
0401 y1 = elamP';
0402 plot(x, y1, 'LineWidth', 2);
0403 % title('Nodal Price', 'FontSize', 16);
0404 ylabel('Nodal Price, $/MWh', 'FontSize', 16);
0405 xlabel('Period', 'FontSize', 16);
0406 axis([0.5 12.5 0 maxp]);
0407 ah = gca;
0408 set(ah, 'FontSize', 14);
0409 ah.XTick = 1:12;
0410 legend('Bus 1', 'Bus 2', 'Bus 3', 'Location', [0.7 0.28 0 0]);
0411 
0412 if nargin > 7 && ~isempty(fname)
0413     h = gcf;
0414     set(h, 'PaperSize', [11 8.5]);
0415     set(h, 'PaperPosition', [0.25 0.25 10.5 8]);
0416     print('-dpdf', fullfile(mypath, sprintf('%s-%d', fname, pp)));
0417 end