0001 function t_most_uc(quiet, create_plots, create_pdfs, savedir)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 if nargin < 4
0018 savedir = '.';
0019 if nargin < 3
0020 create_pdfs = 0;
0021 if nargin < 2
0022 create_plots = 0;
0023 if nargin < 1
0024 quiet = 0;
0025 end
0026 end
0027 end
0028 end
0029 if create_plots
0030 if create_pdfs
0031 fname = 'uc-ex';
0032 else
0033 fname = '';
0034 end
0035 pp = 0;
0036 end
0037
0038 solvers = {'CPLEX', 'GLPK', 'GUROBI', 'MOSEK', 'OT'};
0039 fcn = {'cplex', 'glpk', 'gurobi', 'mosek', 'intlinprog'};
0040
0041
0042
0043
0044 ntests = 66;
0045 t_begin(ntests*length(solvers), quiet);
0046
0047 if quiet
0048 verbose = 0;
0049 else
0050 verbose = 0;
0051 end
0052
0053
0054 if have_feature('octave')
0055 if have_feature('octave', 'vnum') >= 4
0056 file_in_path_warn_id = 'Octave:data-file-in-path';
0057 else
0058 file_in_path_warn_id = 'Octave:load-file-in-path';
0059 end
0060 s1 = warning('query', file_in_path_warn_id);
0061 warning('off', file_in_path_warn_id);
0062 end
0063
0064 casefile = 'ex_case3b';
0065 solnfile = 't_most_uc_soln';
0066 soln = load(solnfile);
0067 mpopt = mpoption;
0068 mpopt = mpoption(mpopt, 'out.gen', 1);
0069 mpopt = mpoption(mpopt, 'verbose', verbose);
0070
0071
0072 mpopt = mpoption(mpopt, 'model', 'DC');
0073 mpopt = mpoption(mpopt, 'most.price_stage_warn_tol', 1e1);
0074
0075
0076 if have_feature('cplex')
0077
0078
0079 mpopt = mpoption(mpopt, 'cplex.lpmethod', 2);
0080
0081
0082 mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.mipgap', 0);
0083 mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.absmipgap', 0);
0084 mpopt = mpoption(mpopt, 'cplex.opts.simplex.tolerances.optimality', 1e-9);
0085 mpopt = mpoption(mpopt, 'cplex.opts.simplex.tolerances.feasibility', 1e-9);
0086 mpopt = mpoption(mpopt, 'cplex.opts.emphasis.numerical', 1);
0087 mpopt = mpoption(mpopt, 'cplex.opts.threads', 2);
0088 end
0089 if have_feature('glpk')
0090 mpopt = mpoption(mpopt, 'glpk.opts.mipgap', 0);
0091 mpopt = mpoption(mpopt, 'glpk.opts.tolint', 1e-10);
0092 mpopt = mpoption(mpopt, 'glpk.opts.tolobj', 1e-10);
0093 end
0094 if have_feature('gurobi')
0095
0096
0097 mpopt = mpoption(mpopt, 'gurobi.method', 1);
0098
0099 mpopt = mpoption(mpopt, 'gurobi.threads', 2);
0100 mpopt = mpoption(mpopt, 'gurobi.opts.MIPGap', 0);
0101 mpopt = mpoption(mpopt, 'gurobi.opts.MIPGapAbs', 0);
0102 end
0103 if have_feature('mosek')
0104 sc = mosek_symbcon;
0105
0106
0107
0108 mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
0109
0110
0111 mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_NODE_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
0112 mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_ROOT_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
0113 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_RELAX_INT', 1e-9);
0114
0115 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_GAP', 0);
0116 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_GAP', 0);
0117 end
0118 if have_feature('intlinprog')
0119
0120
0121
0122 mpopt = mpoption(mpopt, 'linprog.Algorithm', 'dual-simplex');
0123
0124 mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'dual-simplex');
0125 mpopt = mpoption(mpopt, 'intlinprog.TolCon', 1e-9);
0126 mpopt = mpoption(mpopt, 'intlinprog.TolGapAbs', 0);
0127 mpopt = mpoption(mpopt, 'intlinprog.TolGapRel', 0);
0128 mpopt = mpoption(mpopt, 'intlinprog.TolInteger', 1e-6);
0129
0130
0131
0132
0133
0134 end
0135 if ~verbose
0136 mpopt = mpoption(mpopt, 'out.all', 0);
0137 end
0138
0139
0140
0141 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0142 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0143 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0144 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0145 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0146 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0147 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0148 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0149 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0150 [CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
0151 CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
0152 CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
0153 CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
0154 CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
0155 CT_MODCOST_X] = idx_ct;
0156
0157
0158 mpc = loadcase(casefile);
0159
0160 nb = size(mpc.bus, 1);
0161 nl = size(mpc.branch, 1);
0162 ng = size(mpc.gen, 1);
0163
0164 xgd = loadxgendata('ex_xgd_uc', mpc);
0165 [iwind, mpc, xgd] = addwind('ex_wind_uc', mpc, xgd);
0166 profiles = getprofiles('ex_wind_profile_d', iwind);
0167 profiles = getprofiles('ex_load_profile', profiles);
0168 nt = size(profiles(1).values, 1);
0169
0170 mpc_full = mpc;
0171 xgd_full = xgd;
0172
0173 mpc.gencost(:, [STARTUP SHUTDOWN]) = 0;
0174 xgd.MinUp(2) = 1;
0175 xgd.PositiveLoadFollowReserveQuantity(3) = 250;
0176 xgd.PositiveLoadFollowReservePrice(3) = 1e-6;
0177 xgd.NegativeLoadFollowReservePrice(3) = 1e-6;
0178 mpc0 = mpc;
0179 xgd0 = xgd;
0180
0181 for s = 1:length(solvers)
0182 if ~have_feature(fcn{s})
0183 t_skip(ntests, sprintf('%s not installed', solvers{s}));
0184 else
0185 mpopt = mpoption(mpopt, 'opf.dc.solver', solvers{s});
0186 mpopt = mpoption(mpopt, 'most.solver', mpopt.opf.dc.solver);
0187 mpopt = mpoption(mpopt, 'most.storage.cyclic', 1);
0188
0189 t = sprintf('%s : base (econ disp, no network) : ', solvers{s});
0190 mpc = mpc0;
0191 xgd = xgd0;
0192 mpopt = mpoption(mpopt, 'most.dc_model', 0);
0193 mdi = loadmd(mpc, nt, xgd, [], [], profiles);
0194 mdo = most(mdi, mpopt);
0195 ms = most_summary(mdo);
0196 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0197 ex = soln.ed;
0198 t_is(ms.f, ex.f, 8, [t 'f']);
0199 t_is(ms.Pg, ex.Pg, 8, [t 'Pg']);
0200 t_is(ms.Rup, ex.Rup, 8, [t 'Rup']);
0201 t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
0202 t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
0203 t_is(ms.u, ex.u, 8, [t 'u']);
0204 t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
0205 t_is(ms.muF, ex.muF, 8, [t 'muF']);
0206
0207 if create_plots
0208 pp = pp + 1;
0209 plot_case('Base : No Network', mdo, ms, 500, 100, savedir, pp, fname);
0210 end
0211
0212
0213 t = sprintf('%s : + DC OPF constraints : ', solvers{s});
0214 mpc = mpc0;
0215
0216 mpopt = mpoption(mpopt, 'most.dc_model', 1);
0217 mdi = loadmd(mpc, nt, xgd, [], [], profiles);
0218 mdo = most(mdi, mpopt);
0219 ms = most_summary(mdo);
0220 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0221 ex = soln.dcopf;
0222 t_is(ms.f, ex.f, 8, [t 'f']);
0223 t_is(ms.Pg, ex.Pg, 8, [t 'Pg']);
0224 t_is(ms.Rup, ex.Rup, 8, [t 'Rup']);
0225 t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
0226 t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
0227 t_is(ms.u, ex.u, 8, [t 'u']);
0228 t_is(ms.lamP, ex.lamP, 5, [t 'lamP']);
0229 t_is(ms.muF, ex.muF, 8, [t 'muF']);
0230
0231 if create_plots
0232 pp = pp + 1;
0233 plot_case('+ DC Network', mdo, ms, 500, 100, savedir, pp, fname);
0234 end
0235
0236
0237 t = sprintf('%s : + startup/shutdown costs : ', solvers{s});
0238 if mpopt.out.all
0239 fprintf('Add STARTUP and SHUTDOWN costs\n');
0240 end
0241 mpc = mpc_full;
0242
0243
0244
0245
0246
0247
0248 mdi = loadmd(mpc, nt, xgd, [], [], profiles);
0249 mdo = most(mdi, mpopt);
0250 ms = most_summary(mdo);
0251 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0252 ex = soln.wstart;
0253 t_is(ms.f, ex.f, 8, [t 'f']);
0254 t_is(ms.Pg, ex.Pg, 8, [t 'Pg']);
0255 t_is(ms.Rup, ex.Rup, 8, [t 'Rup']);
0256 t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
0257 t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
0258 t_is(ms.u, ex.u, 8, [t 'u']);
0259 t_is(ms.lamP, ex.lamP, 8, [t 'lamP']);
0260 t_is(ms.muF, ex.muF, 8, [t 'muF']);
0261
0262 if create_plots
0263 pp = pp + 1;
0264 plot_case('+ Startup/Shutdown Costs', mdo, ms, 500, 100, savedir, pp, fname);
0265 end
0266
0267
0268 t = sprintf('%s : + min up/down time constraints : ', solvers{s});
0269 if mpopt.out.all
0270 fprintf('Add MinUp time\n');
0271 end
0272 xgd.MinUp(2) = 3;
0273 mdi = loadmd(mpc, nt, xgd, [], [], profiles);
0274 mdo = most(mdi, mpopt);
0275 ms = most_summary(mdo);
0276 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0277 ex = soln.wminup;
0278 t_is(ms.f, ex.f, 6, [t 'f']);
0279 t_is(ms.Pg, ex.Pg, 7, [t 'Pg']);
0280 t_is(ms.Rup, ex.Rup, 7, [t 'Rup']);
0281 t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
0282 t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
0283 t_is(ms.u, ex.u, 8, [t 'u']);
0284 t_is(ms.lamP, ex.lamP, 8, [t 'lamP']);
0285 t_is(ms.muF, ex.muF, 8, [t 'muF']);
0286
0287 if create_plots
0288 pp = pp + 1;
0289 plot_case('+ Min Up/Down Time Constraints', mdo, ms, 500, 100, savedir, pp, fname);
0290 end
0291
0292
0293 t = sprintf('%s : + ramp constraint/ramp res cost : ', solvers{s});
0294 if mpopt.out.all
0295 fprintf('Restrict ramping and add ramp reserve costs\n');
0296 end
0297 xgd = xgd_full;
0298 mdi = loadmd(mpc, nt, xgd, [], [], profiles);
0299 mdo = most(mdi, mpopt);
0300 ms = most_summary(mdo);
0301 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0302 ex = soln.wramp;
0303 t_is(ms.f, ex.f, 8, [t 'f']);
0304 t_is(ms.Pg, ex.Pg, 8, [t 'Pg']);
0305 t_is(ms.Rup, ex.Rup, 8, [t 'Rup']);
0306 t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
0307 t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
0308 t_is(ms.u, ex.u, 8, [t 'u']);
0309 t_is(ms.lamP, ex.lamP, 8, [t 'lamP']);
0310 t_is(ms.muF, ex.muF, 8, [t 'muF']);
0311
0312 if create_plots
0313 pp = pp + 1;
0314 plot_case('+ Ramping Constraints/Ramp Reserve Costs', mdo, ms, 500, 100, savedir, pp, fname);
0315 end
0316
0317
0318 t = sprintf('%s : + storage : ', solvers{s});
0319 if mpopt.out.all
0320 fprintf('Add storage\n');
0321 end
0322 [iess, mpc, xgd, sd] = addstorage('ex_storage', mpc, xgd);
0323 mdi = loadmd(mpc, nt, xgd, sd, [], profiles);
0324 mdo = most(mdi, mpopt);
0325 ms = most_summary(mdo);
0326 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0327 ex = soln.wstorage;
0328 t_is(ms.f, ex.f, 8, [t 'f']);
0329 t_is(ms.Pg, ex.Pg, 8, [t 'Pg']);
0330 t_is(ms.Rup, ex.Rup, 8, [t 'Rup']);
0331 t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
0332 t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
0333 t_is(ms.u, ex.u, 8, [t 'u']);
0334
0335
0336
0337 if create_plots
0338 pp = pp + 1;
0339 plot_case('+ Storage', mdo, ms, 500, 100, savedir, pp, fname);
0340 create_plots = 0;
0341 end
0342
0343
0344 t = sprintf('%s : + storage2 : ', solvers{s});
0345 if mpopt.out.all
0346 fprintf('Add storage\n');
0347 end
0348 mpopt = mpoption(mpopt, 'most.storage.cyclic', 0);
0349 mdi.Storage.rho = 1;
0350 mdo = most(mdi, mpopt);
0351 ms = most_summary(mdo);
0352 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0353 ex = soln.wstorage2;
0354 t_is(ms.f, ex.f, 8, [t 'f']);
0355 t_is(ms.Pg, ex.Pg, 8, [t 'Pg']);
0356 t_is(ms.Rup, ex.Rup, 8, [t 'Rup']);
0357 t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
0358 t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
0359 t_is(ms.u, ex.u, 8, [t 'u']);
0360
0361
0362
0363 if create_plots
0364 pp = pp + 1;
0365 plot_case('+ Storage2', mdo, ms, 500, 100, savedir, pp, fname);
0366 create_plots = 0;
0367 end
0368
0369
0370 t = sprintf('%s : + storage3 : ', solvers{s});
0371 if mpopt.out.all
0372 fprintf('Add storage\n');
0373 end
0374 mdi.Storage.rho = 0;
0375 mdo = most(mdi, mpopt);
0376 ms = most_summary(mdo);
0377 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0378 ex = soln.wstorage3;
0379 t_is(ms.f, ex.f, 8, [t 'f']);
0380 t_is(ms.Pg, ex.Pg, 8, [t 'Pg']);
0381 t_is(ms.Rup, ex.Rup, 8, [t 'Rup']);
0382 t_is(ms.Rdn, ex.Rdn, 8, [t 'Rdn']);
0383 t_is(ms.Pf, ex.Pf, 8, [t 'Pf']);
0384 t_is(ms.u, ex.u, 8, [t 'u']);
0385
0386
0387
0388 if create_plots
0389 pp = pp + 1;
0390 plot_case('+ Storage3', mdo, ms, 500, 100, savedir, pp, fname);
0391 create_plots = 0;
0392 end
0393
0394 end
0395 end
0396
0397 if have_feature('octave')
0398 warning(s1.state, file_in_path_warn_id);
0399 end
0400
0401 t_end;
0402
0403
0404
0405 function h = plot_case(label, md, ms, maxq, maxp, mypath, pp, fname)
0406
0407 if nargin < 8
0408 fname = '';
0409 end
0410
0411
0412 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]};
0413
0414 ig = (1:3)';
0415 id = 4;
0416 iw = 5;
0417 is = 6;
0418
0419 subplot(3, 1, 1);
0420 md.mpc = rmfield(md.mpc, 'genfuel');
0421 plot_uc(md, [], 'title', label);
0422 ylabel('Unit Commitment', 'FontSize', 16);
0423 ah = gca;
0424 ah.YAxisLocation = 'left';
0425
0426 subplot(3, 1, 2);
0427 x = (1:ms.nt)';
0428 y1 = ms.Pg(ig, :)';
0429 if ms.ng == 6
0430 y1 = [y1 max(-ms.Pg(is, :), 0)' max(ms.Pg(is, :), 0)'];
0431 end
0432 y2 = -sum(ms.Pg([id; iw], :), 1)';
0433 [ah1, h1, h2] = plotyy(x, y1, x, y2);
0434 axis(ah1(1), [0.5 12.5 0 maxq]);
0435 axis(ah1(2), [0.5 12.5 0 maxq]);
0436
0437
0438
0439
0440 ah1(1).YTickMode = 'auto';
0441 ah1(2).YTickMode = 'auto';
0442 ah1(1).XTick = 1:12;
0443 nn = 3;
0444 for j = 1:3
0445 h1(j).LineWidth = 2;
0446 h1(j).Color = cc{j};
0447 end
0448 if ms.ng == 6
0449 h1(4).LineWidth = 2;
0450 h1(4).Color = cc{5};
0451 h1(4).LineStyle = ':';
0452 h1(5).LineWidth = 2;
0453 h1(5).Color = cc{5};
0454 end
0455 h2.LineWidth = 2;
0456 h2.Color = cc{4};
0457 h2.LineStyle = ':';
0458 ah1(2).YColor = cc{4};
0459
0460 ylabel(ah1(1), 'Generation, MW', 'FontSize', 16);
0461 ylabel(ah1(2), 'Net Load, MW', 'FontSize', 16);
0462 xlabel('Period', 'FontSize', 16);
0463 set(ah1(1), 'FontSize', 14);
0464 set(ah1(2), 'FontSize', 14);
0465 if ms.ng == 6
0466 legend('Gen 1', 'Gen 2', 'Gen 3', 'Storage Charge', 'Storage Discharge', 'Location', [0.7 0.6 0 0]);
0467 else
0468 legend('Gen 1', 'Gen 2', 'Gen 3', 'Location', [0.7 0.58 0 0]);
0469 end
0470
0471 subplot(3, 1, 3);
0472 y1 = ms.lamP';
0473 plot(x, y1, 'LineWidth', 2);
0474
0475 ylabel('Nodal Price, $/MWh', 'FontSize', 16);
0476 xlabel('Period', 'FontSize', 16);
0477 axis([0.5 12.5 0 maxp]);
0478 ah = gca;
0479 set(ah, 'FontSize', 14);
0480 ah.XTick = 1:12;
0481 legend('Bus 1', 'Bus 2', 'Bus 3', 'Location', [0.7 0.28 0 0]);
0482
0483 if nargin > 7 && ~isempty(fname)
0484 h = gcf;
0485 set(h, 'PaperSize', [11 8.5]);
0486 set(h, 'PaperPosition', [0.25 0.25 10.5 8]);
0487 print('-dpdf', fullfile(mypath, sprintf('%s-%d', fname, pp)));
0488 end