0001 function t_most_spuc(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
0030 solvers = {'CPLEX', 'GLPK', 'GUROBI', 'MOSEK', 'OT'};
0031 fcn = {'cplex', 'glpk', 'gurobi', 'mosek', 'intlinprog'};
0032
0033
0034
0035
0036
0037
0038 ntests = 144;
0039 t_begin(ntests*length(solvers), quiet);
0040
0041 if quiet
0042 verbose = 0;
0043 else
0044 verbose = 0;
0045 end
0046
0047
0048 casefile = 'ex_case3b';
0049 mpopt = mpoption;
0050 mpopt = mpoption(mpopt, 'out.gen', 1);
0051 mpopt = mpoption(mpopt, 'verbose', verbose);
0052
0053
0054 mpopt = mpoption(mpopt, 'model', 'DC');
0055 mpopt = mpoption(mpopt, 'sopf.force_Pc_eq_P0', 0);
0056 mpopt = mpoption(mpopt, 'most.price_stage_warn_tol', 10);
0057
0058
0059 if have_feature('cplex')
0060
0061
0062 mpopt = mpoption(mpopt, 'cplex.lpmethod', 2);
0063
0064
0065 mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.mipgap', 0);
0066 mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.absmipgap', 0);
0067 mpopt = mpoption(mpopt, 'cplex.opts.threads', 2);
0068 end
0069 if have_feature('glpk')
0070 mpopt = mpoption(mpopt, 'glpk.opts.mipgap', 0);
0071 mpopt = mpoption(mpopt, 'glpk.opts.tolint', 1e-10);
0072 mpopt = mpoption(mpopt, 'glpk.opts.tolobj', 1e-10);
0073 end
0074 if have_feature('gurobi')
0075
0076
0077 mpopt = mpoption(mpopt, 'gurobi.method', 1);
0078
0079 mpopt = mpoption(mpopt, 'gurobi.threads', 2);
0080 mpopt = mpoption(mpopt, 'gurobi.opts.MIPGap', 0);
0081 mpopt = mpoption(mpopt, 'gurobi.opts.MIPGapAbs', 0);
0082 end
0083 if have_feature('mosek')
0084 sc = mosek_symbcon;
0085
0086
0087
0088 mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
0089
0090
0091 mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_NODE_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
0092 mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_ROOT_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
0093 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_RELAX_INT', 1e-9);
0094
0095 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_GAP', 0);
0096 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_GAP', 0);
0097 end
0098 if have_feature('intlinprog')
0099
0100
0101
0102 mpopt = mpoption(mpopt, 'linprog.Algorithm', 'dual-simplex');
0103
0104 mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'dual-simplex');
0105 mpopt = mpoption(mpopt, 'intlinprog.TolCon', 1e-9);
0106 mpopt = mpoption(mpopt, 'intlinprog.TolGapAbs', 0);
0107 mpopt = mpoption(mpopt, 'intlinprog.TolGapRel', 0);
0108 mpopt = mpoption(mpopt, 'intlinprog.TolInteger', 1e-6);
0109
0110
0111 mpopt = mpoption(mpopt, 'intlinprog.LPPreprocess', 'none');
0112 end
0113 if ~verbose
0114 mpopt = mpoption(mpopt, 'out.all', 0);
0115 end
0116
0117
0118
0119 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0120 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0121 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0122 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0123 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0124 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0125 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0126 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0127 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0128 [CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
0129 CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
0130 CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
0131 CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
0132 CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
0133 CT_MODCOST_X] = idx_ct;
0134
0135
0136 mpc = loadcase(casefile);
0137 mpc.gencost(:, STARTUP) = 0;
0138 mpc.gencost(:, SHUTDOWN) = 0;
0139
0140
0141 contab = loadgenericdata('ex_contab', 'array');
0142 pp = [1-sum(contab(:,2)); contab(1,2); contab(2,2)];
0143
0144 xgd = loadxgendata('ex_xgd_uc', mpc);
0145 [iwind, mpc, xgd] = addwind('ex_wind_uc', mpc, xgd);
0146 mpc.reserves.zones = [mpc.reserves.zones 0];
0147 mpc = scale_load(499, mpc, [], struct('scale', 'QUANTITY'));
0148 mpc0 = mpc;
0149 xgd0 = xgd;
0150
0151
0152 if create_plots
0153 j = 1;
0154 Pg = NaN(5, 7);
0155 Rp = NaN(5, 7);
0156 Rm = NaN(5, 7);
0157 lamP = NaN(3, 7);
0158 muF = zeros(3, 7);
0159 end
0160
0161 for s = 1:length(solvers)
0162 if ~have_feature(fcn{s})
0163 t_skip(ntests, sprintf('%s not installed', solvers{s}));
0164 else
0165 mpopt = mpoption(mpopt, 'opf.dc.solver', solvers{s});
0166 mpopt = mpoption(mpopt, 'most.solver', mpopt.opf.dc.solver);
0167
0168
0169 if verbose
0170 fprintf('\n\n');
0171 fprintf('--------------------------------------------\n');
0172 fprintf('----- economic dispatch (no network) -----\n');
0173 fprintf('--------------------------------------------\n');
0174 end
0175 t = sprintf('%s : economic dispatch (no network) : runopf ', solvers{s});
0176 mpc = mpc0;
0177 mpc.branch(:, RATE_A) = 0;
0178 r = runuopf(mpc, mpopt);
0179 t_ok(r.success, [t 'success']);
0180 t_is(r.f, -488030, 7, [t 'f']);
0181 t_is(r.gen(:, PG), [200; 199; 0; -499; 100], 7, [t 'Pg']);
0182 t_is(r.gen(:, GEN_STATUS), [1; 1; 0; 1; 1], 7, [t 'u']);
0183 t_is(r.bus(:, LAM_P), [30; 30; 30], 7, [t 'lam P']);
0184
0185
0186 t = sprintf('%s : economic dispatch (no network) : most ', solvers{s});
0187 mpc = mpc0;
0188 mpc.gen(1, GEN_STATUS) = 0;
0189 mpc.gen(2, GEN_STATUS) = 0;
0190
0191 mpopt = mpoption(mpopt, 'most.dc_model', 0);
0192 mdi = loadmd(mpc, [], xgd);
0193 mdo = most(mdi, mpopt);
0194 rr = mdo.flow(1,1,1).mpc;
0195 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0196 t_is(mdo.QP.f, -488030, 7, [t 'f']);
0197 t_is(rr.gen(:, PG), [200; 199; 0; -499; 100], 7, [t 'Pg']);
0198 t_is(rr.gen(:, GEN_STATUS), [1; 1; 0; 1; 1], 7, [t 'u']);
0199
0200 t_is(rr.bus(:, LAM_P), [30; 30; 30], 7, [t 'lam P']);
0201 if create_plots
0202 Pg(:, j) = mdo.results.ExpectedDispatch;
0203 Rp(:, j) = 0;
0204 Rm(:, j) = 0;
0205 lamP(:, j) = rr.bus(:, LAM_P);
0206 j = j + 1;
0207 end
0208
0209
0210
0211 if verbose
0212 fprintf('\n\n');
0213 fprintf('--------------------------------------------\n');
0214 fprintf('----- DC OPF -----\n');
0215 fprintf('--------------------------------------------\n');
0216 end
0217 t = sprintf('%s : DC OPF : runopf ', solvers{s});
0218 mpc = mpc0;
0219
0220 r = runuopf(mpc, mpopt);
0221 t_ok(r.success, [t 'success']);
0222 t_is(r.f, -486040, 7, [t 'f']);
0223 t_is(r.gen(:, PG), [200; 0; 199; -499; 100], 7, [t 'Pg']);
0224 t_is(r.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0225 t_is(r.bus(:, LAM_P), [40; 40; 40], 7, [t 'lam P']);
0226 t_is(r.branch(:, MU_SF) + r.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow']);
0227
0228
0229 t = sprintf('%s : DC OPF : most ', solvers{s});
0230 mpc = mpc0;
0231
0232 mpopt = mpoption(mpopt, 'most.dc_model', 1);
0233 mdi = loadmd(mpc, [], xgd);
0234 mdo = most(mdi, mpopt);
0235 rr = mdo.flow(1,1,1).mpc;
0236 if verbose
0237 printpf(rr, [], mpopt);
0238 end
0239 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0240 t_is(mdo.QP.f, -486040, 7, [t 'f']);
0241 t_is(rr.gen(:, PG), [200; 0; 199; -499; 100], 6, [t 'Pg']);
0242 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0243
0244 t_is(rr.bus(:, LAM_P), [40; 40; 40], 7, [t 'lam P']);
0245 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow']);
0246 if create_plots
0247 Pg(:, j) = mdo.results.ExpectedDispatch;
0248 Rp(:, j) = 0;
0249 Rm(:, j) = 0;
0250 lamP(:, j) = rr.bus(:, LAM_P);
0251 muF(:, j) = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0252 j = j + 1;
0253 end
0254
0255
0256
0257 if verbose
0258 fprintf('\n\n');
0259 fprintf('--------------------------------------------\n');
0260 fprintf('----- economic dispatch (w/reserves) -----\n');
0261 fprintf('--------------------------------------------\n');
0262 end
0263 t = sprintf('%s : economic dispatch (w/reserves) : runopf ', solvers{s});
0264 mpc = mpc0;
0265 mpc.branch(:, RATE_A) = 0;
0266 mpc = toggle_reserves(mpc, 'on');
0267 r = runuopf(mpc, mpopt);
0268 t_ok(r.success, [t 'success']);
0269 t_is(r.f, -486802, 5, [t 'f']);
0270 t_is(r.gen(:, PG), [200; 139; 60; -499; 100], 7, [t 'Pg']);
0271 t_is(r.gen(:, GEN_STATUS), [1; 1; 1; 1; 1], 7, [t 'u']);
0272 t_is(r.bus(:, LAM_P), [32; 32; 32], 7, [t 'lam P']);
0273 t_is(r.reserves.R, [0; 61; 89; 0; 0], 7, [t 'R']);
0274 t_is(r.reserves.prc, [5; 5; 5; 0; 0], 7, [t 'reserve prc']);
0275 t_is(r.reserves.mu.Pmax + r.gen(:, MU_PMAX), [7; 2; 0; 0; 32], 7, [t 'reserve muPmax']);
0276
0277
0278 t = sprintf('%s : economic dispatch (w/reserves) : most ', solvers{s});
0279 mpc = mpc0;
0280
0281 mpopt = mpoption(mpopt, 'most.dc_model', 0);
0282 mdi = loadmd(mpc, [], xgd);
0283 mdi.FixedReserves = mpc.reserves;
0284 mdo = most(mdi, mpopt);
0285 rr = mdo.flow(1,1,1).mpc;
0286 if verbose
0287 printpf(rr, [], mpopt);
0288 end
0289 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0290 t_is(mdo.QP.f, -486802, 5, [t 'f']);
0291 t_is(rr.gen(:, PG), [200; 139; 60; -499; 100], 7, [t 'Pg']);
0292 t_is(rr.gen(:, GEN_STATUS), [1; 1; 1; 1; 1], 7, [t 'u']);
0293
0294 t_is(rr.bus(:, LAM_P), [32; 32; 32], 7, [t 'lam P']);
0295 t_is(rr.reserves.R, [0; 61; 89; 0; 0], 7, [t 'R']);
0296 t_is(rr.reserves.prc, [5; 5; 5; 0; 0], 7, [t 'reserve prc']);
0297 t_is(rr.reserves.mu.Pmax + rr.gen(:, MU_PMAX), [7; 2; 0; 0; 32], 7, [t 'reserve muPmax']);
0298 if create_plots
0299 Pg(:, j) = mdo.results.ExpectedDispatch;
0300 Rp(:, j) = rr.reserves.R;
0301 Rm(:, j) = 0;
0302 lamP(:, j) = rr.bus(:, LAM_P);
0303 j = j + 1;
0304 end
0305
0306
0307
0308 if verbose
0309 fprintf('\n\n');
0310 fprintf('--------------------------------------------\n');
0311 fprintf('----- DC OPF (w/reserves) -----\n');
0312 fprintf('--------------------------------------------\n');
0313 end
0314 t = sprintf('%s : DC OPF (w/reserves) : runopf ', solvers{s});
0315 mpc = mpc0;
0316 r = runopf_w_res(mpc, mpopt);
0317 t_ok(r.success, [t 'success']);
0318 t_is(r.f, -485656, 5, [t 'f']);
0319 t_is(r.gen(:, PG), [156; 65; 178; -499; 100], 7, [t 'Pg']);
0320 t_is(r.gen(:, GEN_STATUS), [1; 1; 1; 1; 1], 7, [t 'u']);
0321 t_is(r.bus(:, LAM_P), [29; 40; 51], 7, [t 'lam P']);
0322 t_is(r.branch(:, MU_SF) + r.branch(:, MU_ST), [0; 33; 0], 7, [t 'mu flow']);
0323 t_is(r.reserves.R, [44; 100; 6; 0; 0], 7, [t 'R']);
0324 t_is(r.reserves.prc, [5; 5; 5; 0; 0], 7, [t 'reserve prc']);
0325 t_is(r.reserves.mu.Pmax + r.gen(:, MU_PMAX), [4; 0; 0; 0; 40], 7, [t 'reserve muPmax']);
0326
0327
0328 t = sprintf('%s : DC OPF (w/reserves) : most ', solvers{s});
0329 mpc = mpc0;
0330
0331 mpopt = mpoption(mpopt, 'most.dc_model', 1);
0332 mdi = loadmd(mpc, [], xgd);
0333 mdi.FixedReserves = mpc.reserves;
0334 mdo = most(mdi, mpopt);
0335 rr = mdo.flow(1,1,1).mpc;
0336 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0337 t_is(mdo.QP.f, -485656, 5, [t 'f']);
0338 t_is(rr.gen(:, PG), [156; 65; 178; -499; 100], 6, [t 'Pg']);
0339 t_is(rr.gen(:, GEN_STATUS), [1; 1; 1; 1; 1], 7, [t 'u']);
0340
0341 t_is(rr.bus(:, LAM_P), [29; 40; 51], 6, [t 'lam P']);
0342 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 33; 0], 7, [t 'mu flow']);
0343 t_is(rr.reserves.R, [44; 100; 6; 0; 0], 6, [t 'R']);
0344 t_is(rr.reserves.prc, [5; 5; 5; 0; 0], 7, [t 'reserve prc']);
0345 t_is(rr.reserves.mu.Pmax + rr.gen(:, MU_PMAX), [4; 0; 0; 0; 40], 7, [t 'reserve muPmax']);
0346 if create_plots
0347 Pg(:, j) = mdo.results.ExpectedDispatch;
0348 Rp(:, j) = rr.reserves.R;
0349 Rm(:, j) = 0;
0350 lamP(:, j) = rr.bus(:, LAM_P);
0351 muF(:, j) = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0352 j = j + 1;
0353 end
0354
0355 t = sprintf('%s : Secure DC OPF (w/cont,res,ramp) : c3sopf ', solvers{s});
0356 mpc = mpc0;
0357 mpc = scale_load(350, mpc, [], struct('scale', 'QUANTITY'));
0358 xgd_table.colnames = {
0359 'PositiveActiveReservePrice', ...
0360 'PositiveActiveReserveQuantity', ...
0361 'NegativeActiveReservePrice', ...
0362 'NegativeActiveReserveQuantity', ...
0363 'PositiveActiveDeltaPrice', ...
0364 'NegativeActiveDeltaPrice', ...
0365 };
0366 xgd_table.data = [
0367 5 250 10 0 1e-9 1e-9;
0368 1e-8 100 2e-8 0 1e-9 1e-9;
0369 1.5 600 3 0 1e-9 1e-9;
0370 1e-8 800 2e-8 0 1e-9 1e-9;
0371 1e-8 200 2e-8 0 1e-9 1e-9;
0372 ];
0373
0374 xgd.PositiveActiveReserveQuantity(2) = 100;
0375 xgd.PositiveActiveReserveQuantity(4) = 800;
0376 xgd.NegativeActiveReserveQuantity(:) = 0;
0377 xgd.PositiveActiveReservePrice([1;3]) = [5; 1.5];
0378 xgd.NegativeActiveReservePrice([1;3]) = [10; 3];
0379
0380 mpc1 = mpc;
0381 mpc1.gen(2, GEN_STATUS) = 0;
0382 r = c3sopf(mpc1, xgd_table.data, contab, mpopt);
0383
0384
0385
0386 t_ok(r.success, [t 'success']);
0387 t_is(r.opf_results.f, -340850, 4, [t 'f']);
0388 rr = r.base;
0389
0390
0391
0392 t_is(rr.gen(:, PG), [190; 0; 60; -350; 100], 6, [t 'Pg base']);
0393 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0394 t_is(rr.bus(:, LAM_P), [25; 25; 25], 7, [t 'lam P base']);
0395 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow base']);
0396 rr = r.cont(1);
0397
0398
0399
0400 t_ok(isempty(rr.gen), [t 'Pg 1']);
0401 t_ok(isempty(rr.bus), [t 'lam P 1']);
0402 t_ok(isempty(rr.branch), [t 'mu flow 1']);
0403 rr = r.cont(2);
0404
0405
0406
0407 t_is(rr.gen(:, PG), [190; 0; 60; -300; 50], 5, [t 'Pg 2']);
0408 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0409 t_is(rr.bus(:, LAM_P), [0; 0; 40], 7, [t 'lam P 2']);
0410 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 40], 7, [t 'mu flow 2']);
0411
0412 t = sprintf('%s : Secure DC OPF (w/cont,res,ramp) : most ', solvers{s});
0413 mdi = loadmd(mpc, [], xgd, [], contab);
0414 mdo = most(mdi, mpopt);
0415
0416
0417
0418 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0419 t_is(mdo.QP.f, -340850, 5, [t 'f']);
0420 rr = mdo.flow(1,1,1).mpc;
0421
0422
0423
0424
0425 t_is(rr.gen(:, PG), [190; 0; 60; -350; 100], 6, [t 'Pg base']);
0426 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0427
0428 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow base']);
0429 if create_plots
0430 lamP(:, j) = rr.bus(:, LAM_P);
0431 muF(:, j) = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0432 end
0433
0434 rr1 = rr;
0435 rr = mdo.flow(1,1,2).mpc;
0436
0437
0438
0439
0440 t_is(rr.gen(:, PG), [190; 0; 60; -350; 100], 5, [t 'Pg 1']);
0441 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0442 t_is(rr1.bus(:, LAM_P) + rr.bus(:, LAM_P), [25; 25; 25], 7, [t 'lam P base + lam P 1']);
0443 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 1']);
0444 if create_plots
0445 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0446 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0447 end
0448
0449 rr = mdo.flow(1,1,3).mpc;
0450
0451
0452
0453
0454 t_is(rr.gen(:, PG), [190; 0; 60; -300; 50], 5, [t 'Pg 2']);
0455 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0456 t_is(rr.bus(:, LAM_P), [0; 0; 40], 7, [t 'lam P 2']);
0457 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 40], 7, [t 'mu flow 2']);
0458 if create_plots
0459 Pg(:, j) = mdo.results.ExpectedDispatch;
0460 Rp(:, j) = mdo.results.Rpp + mdo.results.Pc - mdo.results.ExpectedDispatch;
0461 Rm(:, j) = mdo.results.Rpm - mdo.results.Pc + mdo.results.ExpectedDispatch;
0462 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0463 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0464 j = j + 1;
0465 end
0466
0467 t = sprintf('%s : Stochastic DC OPF (w/wind,res) : c3sopf ', solvers{s});
0468 nt = 1;
0469 nj = 3;
0470
0471 transmat = {[0.158655253931457; 0.682689492137086; 0.158655253931457]};
0472 profiles = getprofiles(uniformwindprofile(nt, nj), iwind);
0473
0474 pp = transmat{1};
0475
0476
0477 contab0 = [
0478 1 pp(1) profiles.table profiles.rows profiles.col profiles.chgtype profiles.values(1)/profiles.values(2);
0479 3 pp(3) profiles.table profiles.rows profiles.col profiles.chgtype profiles.values(3)/profiles.values(2);
0480 ];
0481 mpc1 = mpc;
0482 mpc1.gen(iwind, PMAX) = mpc1.gen(iwind, PMAX) * profiles.values(2);
0483 mpc1.gen(3, GEN_STATUS) = 0;
0484 wind = loadgenericdata('ex_wind_uc', 'struct', 'gen', 'wind');
0485 offers = [xgd_table.data; wind.xgd_table.data(:, [3:8])];
0486
0487
0488 r = c3sopf(mpc1, offers, contab0, mpopt);
0489 t_ok(r.success, [t 'success']);
0490 t_is(r.opf_results.f, -341928.605134, 5, [t 'f']);
0491 rr = r.cont(1);
0492
0493
0494
0495 t_is(rr.gen(:, PG), [200; 150; 0; -350; 0], 7, [t 'Pg 1']);
0496 t_is(rr.bus(:, LAM_P), [4.7596576; 4.7596576; 4.7596576], 7, [t 'lam P 1']);
0497 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 1']);
0498 rr = r.base;
0499 t_is(rr.gen(:, PG), [200; 100; 0; -350; 50], 5, [t 'Pg 2']);
0500 t_is(rr.bus(:, LAM_P), [20.4806848; 20.4806848; 20.4806848], 6, [t 'lam P 2']);
0501 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2']);
0502 rr = r.cont(2);
0503 t_is(rr.gen(:, PG), [200; 65; 0; -350; 85], 5, [t 'Pg 3']);
0504 t_is(rr.bus(:, LAM_P), [0; 0; 0], 6, [t 'lam P 3']);
0505 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 3']);
0506
0507 t = sprintf('%s : Stochastic DC OPF (w/wind,res) : most ', solvers{s});
0508 mdi = loadmd(mpc, transmat, xgd, [], [], profiles);
0509 mdo = most(mdi, mpopt);
0510
0511
0512
0513 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0514 t_is(mdo.QP.f, -341928.605134, 6, [t 'f']);
0515 rr = mdo.flow(1,1,1).mpc;
0516
0517
0518
0519
0520 t_is(rr.gen(:, PG), [200; 150; 0; -350; 0], 7, [t 'Pg 1']);
0521 t_is(rr.gen(:, GEN_STATUS), [1; 1; 0; 1; 1], 7, [t 'u']);
0522 t_is(rr.bus(:, LAM_P), [4.7596576; 4.7596576; 4.7596576], 7, [t 'lam P 1']);
0523 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 1']);
0524 if create_plots
0525 lamP(:, j) = rr.bus(:, LAM_P);
0526 muF(:, j) = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0527 end
0528
0529 rr = mdo.flow(1,2,1).mpc;
0530
0531
0532
0533
0534 t_is(rr.gen(:, PG), [200; 100; 0; -350; 50], 7, [t 'Pg 2']);
0535 t_is(rr.gen(:, GEN_STATUS), [1; 1; 0; 1; 1], 7, [t 'u']);
0536 t_is(rr.bus(:, LAM_P), [20.4806848; 20.4806848; 20.4806848], 6, [t 'lam P 2']);
0537 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2']);
0538 if create_plots
0539 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0540 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0541 end
0542
0543 rr = mdo.flow(1,3,1).mpc;
0544
0545
0546
0547
0548 t_is(rr.gen(:, PG), [200; 65; 0; -350; 85], 5, [t 'Pg 3']);
0549 t_is(rr.gen(:, GEN_STATUS), [1; 1; 0; 1; 1], 7, [t 'u']);
0550 t_is(rr.bus(:, LAM_P), [0; 0; 0], 6, [t 'lam P 3']);
0551 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 3']);
0552 if create_plots
0553 Pg(:, j) = mdo.results.ExpectedDispatch;
0554 Rp(:, j) = mdo.results.Rpp + mdo.results.Pc - mdo.results.ExpectedDispatch;
0555 Rm(:, j) = mdo.results.Rpm - mdo.results.Pc + mdo.results.ExpectedDispatch;
0556 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0557 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0558 j = j + 1;
0559 end
0560
0561
0562 t = sprintf('%s : Secure Stochastic DC OPF (w/wind,cont,res,ramp) : most ', solvers{s});
0563 mdi = loadmd(mpc, transmat, xgd, [], contab, profiles);
0564 mdo = most(mdi, mpopt);
0565
0566
0567
0568 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0569 t_is(mdo.QP.f, -338857.9224447, 4, [t 'f']);
0570
0571
0572 rr = mdo.flow(1,1,1).mpc;
0573
0574
0575
0576
0577 t_is(rr.gen(:, PG), [200; 0; 150; -350; 0], 7, [t 'Pg 1 base']);
0578 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0579
0580 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 2, [t 'mu flow 1 base']);
0581 if create_plots
0582 lamP(:, j) = rr.bus(:, LAM_P);
0583 muF(:, j) = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0584 end
0585
0586 rr1 = rr;
0587 rr = mdo.flow(1,1,2).mpc;
0588
0589
0590
0591
0592 t_is(rr.gen(:, PG), [200; 0; 150; -350; 0], 4, [t 'Pg 1 1']);
0593 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0594
0595 t_is(rr1.bus(:, LAM_P) + rr.bus(:, LAM_P), [7.5923617; 7.5923617; 7.5923617], 7, [t 'lam P 1 1']);
0596 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 1 1']);
0597 if create_plots
0598 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0599 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0600 end
0601
0602 rr = mdo.flow(1,1,3).mpc;
0603
0604
0605
0606
0607 t_is(rr.gen(:, PG), [200; 0; 100; -300; 0], 4, [t 'Pg 1 2']);
0608 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0609 t_is(rr.bus(:, LAM_P), [0.2538484; 0.2538484; 6.3462102], 6, [t 'lam P 1 2']);
0610 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 6.0923618], 6, [t 'mu flow 1 2']);
0611 if create_plots
0612 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0613 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0614 end
0615
0616 rr = mdo.flow(1,2,1).mpc;
0617
0618
0619
0620
0621 t_is(rr.gen(:, PG), [200; 0; 100; -350; 50], 6, [t 'Pg 2 base']);
0622 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0623 t_is(rr.bus(:, LAM_P), [24.5768217; 24.5768217; 24.5768217], 6, [t 'lam P 2 base']);
0624 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2 base']);
0625 if create_plots
0626 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0627 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0628 end
0629
0630 rr = mdo.flow(1,2,2).mpc;
0631
0632
0633
0634
0635 t_is(rr.gen(:, PG), [200; 0; 100; -350; 50], 6, [t 'Pg 2 1']);
0636 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0637 t_is(rr.bus(:, LAM_P), [1.6384548; 1.6384548; 1.6384548], 6, [t 'lam P 2 1']);
0638 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2 1']);
0639 if create_plots
0640 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0641 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0642 end
0643
0644 rr = mdo.flow(1,2,3).mpc;
0645
0646
0647
0648
0649 t_is(rr.gen(:, PG), [200; 0; 60; -300; 40], 6, [t 'Pg 2 2']);
0650 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0651 t_is(rr.bus(:, LAM_P), [0; 0; 27.3075797], 5, [t 'lam P 2 2']);
0652 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 27.3075797], 6, [t 'mu flow 2 2']);
0653 if create_plots
0654 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0655 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0656 end
0657
0658 rr = mdo.flow(1,3,1).mpc;
0659
0660
0661
0662
0663 t_is(rr.gen(:, PG), [200; 0; 60; -350; 90], 6, [t 'Pg 3 base']);
0664 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0665 t_is(rr.bus(:, LAM_P), [0; 0; 0], 6, [t 'lam P 3 base']);
0666 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 5, [t 'mu flow 3 base']);
0667 if create_plots
0668 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0669 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0670 end
0671
0672 rr = mdo.flow(1,3,2).mpc;
0673
0674
0675
0676
0677 t_is(rr.gen(:, PG), [200; 0; 60; -350; 90], 6, [t 'Pg 3 1']);
0678 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0679 t_is(rr.bus(:, LAM_P), [0; 0; 0], 6, [t 'lam P 3 1']);
0680 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 3 1']);
0681 if create_plots
0682 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0683 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0684 end
0685
0686 rr = mdo.flow(1,3,3).mpc;
0687
0688
0689
0690
0691 t_is(rr.gen(:, PG), [200; 0; 60; -300; 40], 4, [t 'Pg 3 2']);
0692 t_is(rr.gen(:, GEN_STATUS), [1; 0; 1; 1; 1], 7, [t 'u']);
0693 t_is(rr.bus(:, LAM_P), [0; 0; 6.3462102], 5, [t 'lam P 3 2']);
0694 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 6.3462102], 6, [t 'mu flow 3 2']);
0695
0696 end
0697
0698 if create_plots
0699 create_plots = 0;
0700
0701 Pg(:, j) = mdo.results.ExpectedDispatch;
0702 Rp(:, j) = mdo.results.Rpp + mdo.results.Pc - mdo.results.ExpectedDispatch;
0703 Rm(:, j) = mdo.results.Rpm - mdo.results.Pc + mdo.results.ExpectedDispatch;
0704 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0705 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0706
0707 R = Rp + Rm;
0708
0709 labels = {'Economic Dispatch'; 'DC OPF'; 'Economic Dispatch w/reserves'; 'DC OPF w/reserves'; 'secure DC OPF'; 'stochastic DC OPF'; 'secure stochastic DC OPF'};
0710
0711 figure(1);
0712
0713 subplot(4, 1, 1);
0714 bar(abs(Pg([1:3 5],:)'),'stacked');
0715 title('Generation');
0716 ylabel('MW');
0717 h = gca;
0718 h.XTickLabel = labels;
0719 h.XTickLabelRotation = 20;
0720
0721 subplot(4, 1, 2);
0722 bar(abs(R([1:3 5],:)'),'stacked');
0723 title('Reserves');
0724 ylabel('MW');
0725 h = gca;
0726 h.XTickLabel = labels;
0727 h.XTickLabelRotation = 20;
0728
0729 subplot(4, 1, 3);
0730 plot(lamP');
0731 title('Nodal Prices');
0732 ylabel('$/MW');
0733 h = gca;
0734 h.XTickLabel = {'', labels{:}, ''}';
0735 h.XTickLabelRotation = 20;
0736 h = [0 8 0 100];
0737 axis(h);
0738
0739 subplot(4, 1, 4);
0740 plot(muF');
0741 title('Flow Constraint Shadow Prices');
0742 ylabel('$/MW');
0743 h = gca;
0744 h.XTickLabel = {'', labels{:}, ''}';
0745 h.XTickLabelRotation = 20;
0746 h = [0 8 0 60];
0747 axis(h);
0748
0749 if create_pdfs
0750 fname = 'single-period-uc';
0751 h = gcf;
0752 set(h, 'PaperSize', [11 8.5]);
0753 set(h, 'PaperPosition', [0.25 0.25 10.5 8]);
0754 print('-dpdf', fullfile(savedir, fname));
0755 end
0756
0757 for j = 1:7;
0758 figure(j+1);
0759 if create_pdfs
0760 fname = sprintf('single-period-uc-%d', j);
0761 else
0762 fname = '';
0763 end
0764 plot_case(labels{j}, Pg([1:3 5], j), Rp([1:3 5], j), Rm([1:3 5], j), lamP(:, j), muF(:, j), 250, 100, savedir, fname);
0765 end
0766 end
0767 end
0768
0769 t_end;
0770
0771
0772 function h = plot_case(label, Pg, Rp, Rm, lamP, muF, maxq, maxp, mypath, fname)
0773
0774 subplot(1, 3, 1);
0775 h = bar([Pg-Rm Rm Rp], 'stacked');
0776 set(h(2), 'FaceColor', [0 0.35 0.33]);
0777 ah1 = gca;
0778 title('Generation & Reserves');
0779 ylabel('MW');
0780 xlabel('Gen');
0781 set(gca, 'FontSize', 14);
0782
0783 if nargin < 6
0784 maxq = ah1.YLim(2);
0785 end
0786 ah1.YLim(2) = maxq;
0787 ah1.YLim(1) = 0;
0788
0789 subplot(1, 3, 2);
0790 bar(lamP);
0791 ah3 = gca;
0792 title('Nodal Prices');
0793 ylabel('$/MW');
0794 xlabel('Bus');
0795 set(gca, 'FontSize', 14);
0796
0797 subplot(1, 3, 3);
0798 bar(muF);
0799 ah4 = gca;
0800 title('Flow Constraint Shadow Prices');
0801 ylabel('$/MW');
0802 xlabel('Line');
0803 set(gca, 'FontSize', 14);
0804
0805 if nargin < 7
0806 maxp = max(ah3.YLim(2), ah4.YLim(2));
0807 end
0808 ah3.YLim(1) = 0;
0809 ah4.YLim(1) = 0;
0810 ah3.YLim(2) = maxp;
0811 ah4.YLim(2) = maxp;
0812
0813 [ax,h] = suplabel(label, 't');
0814 set(h, 'FontSize', 18)
0815
0816 if nargin > 7 && ~isempty(fname)
0817 h = gcf;
0818 set(h, 'PaperSize', [11 8.5]);
0819 set(h, 'PaperPosition', [0.25 0.25 10.5 8]);
0820 print('-dpdf', fullfile(mypath, fname));
0821 end