0001 function t_most_sp(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 t_begin(156, quiet);
0031
0032 if quiet
0033 verbose = 0;
0034 else
0035 verbose = 0;
0036 end
0037
0038
0039 casefile = 'ex_case3a';
0040 mpopt = mpoption;
0041 mpopt = mpoption(mpopt, 'out.gen', 1);
0042 mpopt = mpoption(mpopt, 'verbose', verbose);
0043
0044
0045 mpopt = mpoption(mpopt, 'model', 'DC');
0046 mpopt = mpoption(mpopt, 'sopf.force_Pc_eq_P0', 0);
0047 mpopt = mpoption(mpopt, 'opf.dc.solver', 'MIPS');
0048 mpopt = mpoption(mpopt, 'most.solver', mpopt.opf.dc.solver);
0049 if ~verbose
0050 mpopt = mpoption(mpopt, 'out.all', 0);
0051 end
0052
0053
0054
0055 if have_feature('octave')
0056 s = warning('query', 'Octave:nearly-singular-matrix');
0057 warning('off', 'Octave:nearly-singular-matrix');
0058 end
0059
0060
0061 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0062 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0063 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0064 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0065 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0066 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0067 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0068 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0069 [CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
0070 CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
0071 CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
0072 CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
0073 CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
0074 CT_MODCOST_X] = idx_ct;
0075
0076
0077 mpc = loadcase(casefile);
0078
0079
0080 xgd_table = loadgenericdata('ex_xgd', 'struct');
0081
0082
0083 wind = loadgenericdata('ex_wind', 'struct');
0084
0085
0086 contab = loadgenericdata('ex_contab', 'array');
0087 pp = [1-sum(contab(:,2)); contab(1,2); contab(2,2)];
0088
0089 xgd = loadxgendata(xgd_table, mpc);
0090
0091 mpc0 = mpc;
0092 xgd0 = xgd;
0093
0094
0095 if create_plots
0096 j = 1;
0097 Pg = NaN(5, 7);
0098 Rp = NaN(5, 7);
0099 Rm = NaN(5, 7);
0100 lamP = NaN(3, 7);
0101 muF = zeros(3, 7);
0102 end
0103
0104
0105
0106 if verbose
0107 fprintf('\n\n');
0108 fprintf('--------------------------------------------\n');
0109 fprintf('----- economic dispatch (no network) -----\n');
0110 fprintf('--------------------------------------------\n');
0111 end
0112 t = 'economic dispatch (no network) : runopf ';
0113 mpc.branch(:, RATE_A) = 0;
0114 r = runopf(mpc, mpopt);
0115 t_ok(r.success, [t 'success']);
0116 t_is(r.f, -443250, 5, [t 'f']);
0117 t_is(r.gen(:, PG), [150; 150; 150; -450], 7, [t 'Pg']);
0118 t_is(r.bus(:, LAM_P), [30; 30; 30], 7, [t 'lam P']);
0119
0120
0121 t = 'economic dispatch (no network) : most ';
0122 mpc = mpc0;
0123 xgd = xgd0;
0124 mpopt = mpoption(mpopt, 'most.dc_model', 0);
0125 mdi = loadmd(mpc);
0126 mdo = most(mdi, mpopt);
0127 rr = mdo.flow(1,1,1).mpc;
0128 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0129 t_is(mdo.QP.f, -443250, 5, [t 'f']);
0130 t_is(rr.gen(:, PG), [150; 150; 150; -450], 7, [t 'Pg']);
0131 t_is(rr.bus(:, LAM_P), [30; 30; 30], 7, [t 'lam P']);
0132 if create_plots
0133 Pg(1:4, j) = mdo.results.ExpectedDispatch;
0134 Rp(1:4, j) = 0;
0135 Rm(1:4, j) = 0;
0136 lamP(:, j) = rr.bus(:, LAM_P);
0137 j = j + 1;
0138 end
0139
0140
0141
0142 if verbose
0143 fprintf('\n\n');
0144 fprintf('--------------------------------------------\n');
0145 fprintf('----- DC OPF -----\n');
0146 fprintf('--------------------------------------------\n');
0147 end
0148 t = 'DC OPF : runopf ';
0149 mpc = mpc0;
0150 r = runopf(mpc, mpopt);
0151 t_ok(r.success, [t 'success']);
0152 t_is(r.f, -443115, 5, [t 'f']);
0153 t_is(r.gen(:, PG), [135; 135; 180; -450], 7, [t 'Pg']);
0154 t_is(r.bus(:, LAM_P), [27; 36; 45], 7, [t 'lam P']);
0155 t_is(r.branch(:, MU_SF) + r.branch(:, MU_ST), [0; 27; 0], 7, [t 'mu flow']);
0156
0157
0158 t = 'DC OPF : most ';
0159 mpc = mpc0;
0160 mpopt = mpoption(mpopt, 'most.dc_model', 1);
0161 mdi = loadmd(mpc);
0162 mdo = most(mdi, mpopt);
0163 rr = mdo.flow(1,1,1).mpc;
0164 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0165 t_is(mdo.QP.f, -443115, 5, [t 'f']);
0166 t_is(rr.gen(:, PG), [135; 135; 180; -450], 6, [t 'Pg']);
0167 t_is(rr.bus(:, LAM_P), [27; 36; 45], 7, [t 'lam P']);
0168 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 27; 0], 6, [t 'mu flow']);
0169 if create_plots
0170 Pg(1:4, j) = mdo.results.ExpectedDispatch;
0171 Rp(1:4, j) = 0;
0172 Rm(1:4, j) = 0;
0173 lamP(:, j) = rr.bus(:, LAM_P);
0174 muF(:, j) = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0175 j = j + 1;
0176 end
0177
0178
0179
0180 if verbose
0181 fprintf('\n\n');
0182 fprintf('--------------------------------------------\n');
0183 fprintf('----- economic dispatch (w/reserves) -----\n');
0184 fprintf('--------------------------------------------\n');
0185 end
0186 t = 'economic dispatch (w/reserves) : runopf ';
0187 mpc = mpc0;
0188 mpc.branch(:, RATE_A) = 0;
0189 r = runopf_w_res(mpc, mpopt);
0190 t_ok(r.success, [t 'success']);
0191 t_is(r.f, -442820, 5, [t 'f']);
0192 t_is(r.gen(:, PG), [140; 150; 160; -450], 7, [t 'Pg']);
0193 t_is(r.bus(:, LAM_P), [32; 32; 32], 7, [t 'lam P']);
0194 t_is(r.reserves.R, [60; 50; 40; 0], 7, [t 'R']);
0195 t_is(r.reserves.prc, [5; 5; 5; 0], 7, [t 'reserve prc']);
0196 t_is(r.reserves.mu.Pmax + r.gen(:, MU_PMAX), [4; 2; 0; 0], 7, [t 'reserve muPmax']);
0197
0198
0199 t = 'economic dispatch (w/reserves) : most ';
0200 mpc = mpc0;
0201 mpopt = mpoption(mpopt, 'most.dc_model', 0);
0202 mdi = loadmd(mpc);
0203 mdi.FixedReserves = mpc.reserves;
0204 mdo = most(mdi, mpopt);
0205 rr = mdo.flow(1,1,1).mpc;
0206 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0207 t_is(mdo.QP.f, -442820, 5, [t 'f']);
0208 t_is(rr.gen(:, PG), [140; 150; 160; -450], 7, [t 'Pg']);
0209 t_is(rr.bus(:, LAM_P), [32; 32; 32], 7, [t 'lam P']);
0210 t_is(rr.reserves.R, [60; 50; 40; 0], 7, [t 'R']);
0211 t_is(rr.reserves.prc, [5; 5; 5; 0], 7, [t 'reserve prc']);
0212 t_is(rr.reserves.mu.Pmax + rr.gen(:, MU_PMAX), [4; 2; 0; 0], 7, [t 'reserve muPmax']);
0213 if create_plots
0214 Pg(1:4, j) = mdo.results.ExpectedDispatch;
0215 Rp(1:4, j) = rr.reserves.R;
0216 Rm(1:4, j) = 0;
0217 lamP(:, j) = rr.bus(:, LAM_P);
0218 j = j + 1;
0219 end
0220
0221
0222
0223 if verbose
0224 fprintf('\n\n');
0225 fprintf('--------------------------------------------\n');
0226 fprintf('----- DC OPF (w/reserves) -----\n');
0227 fprintf('--------------------------------------------\n');
0228 end
0229 t = 'DC OPF (w/reserves) : runopf ';
0230 mpc = mpc0;
0231 r = runopf_w_res(mpc, mpopt);
0232 t_ok(r.success, [t 'success']);
0233 t_is(r.f, -442760, 5, [t 'f']);
0234 t_is(r.gen(:, PG), [130; 140; 180; -450], 7, [t 'Pg']);
0235 t_is(r.bus(:, LAM_P), [30; 36; 42], 7, [t 'lam P']);
0236 t_is(r.branch(:, MU_SF) + r.branch(:, MU_ST), [0; 18; 0], 7, [t 'mu flow']);
0237 t_is(r.reserves.R, [70; 60; 20; 0], 7, [t 'R']);
0238 t_is(r.reserves.prc, [5; 5; 5; 0], 7, [t 'reserve prc']);
0239 t_is(r.reserves.mu.Pmax + r.gen(:, MU_PMAX), [4; 2; 0; 0], 7, [t 'reserve muPmax']);
0240
0241
0242 t = 'DC OPF (w/reserves) : most ';
0243 mpc = mpc0;
0244 mpopt = mpoption(mpopt, 'most.dc_model', 1);
0245 mdi = loadmd(mpc);
0246 mdi.FixedReserves = mpc.reserves;
0247 mdo = most(mdi, mpopt);
0248 rr = mdo.flow(1,1,1).mpc;
0249 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0250 t_is(mdo.QP.f, -442760, 5, [t 'f']);
0251 t_is(rr.gen(:, PG), [130; 140; 180; -450], 6, [t 'Pg']);
0252 t_is(rr.bus(:, LAM_P), [30; 36; 42], 6, [t 'lam P']);
0253 t_is(rr.reserves.R, [70; 60; 20; 0], 6, [t 'R']);
0254 t_is(rr.reserves.prc, [5; 5; 5; 0], 7, [t 'reserve prc']);
0255 t_is(rr.reserves.mu.Pmax + rr.gen(:, MU_PMAX), [4; 2; 0; 0], 7, [t 'reserve muPmax']);
0256 if create_plots
0257 Pg(1:4, j) = mdo.results.ExpectedDispatch;
0258 Rp(1:4, j) = rr.reserves.R;
0259 Rm(1:4, j) = 0;
0260 lamP(:, j) = rr.bus(:, LAM_P);
0261 muF(:, j) = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0262 j = j + 1;
0263 end
0264
0265
0266
0267 if verbose
0268 fprintf('\n\n');
0269 fprintf('--------------------------------------------\n');
0270 fprintf('----- DC OPF (w/contingencies) -----\n');
0271 fprintf('--------------------------------------------\n');
0272 end
0273 t = 'DC OPF (w/contingencies) : runopf ';
0274 ff = [-443115; -439750; -297000];
0275 mpc = mpc0;
0276 r = runopf(mpc, mpopt);
0277 t_is(r.f, ff(1), 5, [t 'f']);
0278 t_is(r.gen(:, PG), [135; 135; 180; -450], 7, [t 'Pg']);
0279 t_is(r.bus(:, LAM_P), [27; 36; 45], 7, [t 'lam P']);
0280 t_is(r.branch(:, MU_SF) + r.branch(:, MU_ST), [0; 27; 0], 7, [t 'mu flow']);
0281
0282 mpc = mpc0;
0283 mpc.gen(2, GEN_STATUS) = 0;
0284 r = runopf(mpc, mpopt);
0285 t_is(r.f, ff(2), 5, [t 'f']);
0286 t_ok(r.success, [t 'success']);
0287 t_is(r.gen(:, PG), [200; 0; 250; -450], 7, [t 'Pg']);
0288 t_is(r.bus(:, LAM_P), [50; 50; 50], 7, [t 'lam P']);
0289 t_is(r.branch(:, MU_SF) + r.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow']);
0290
0291 mpc = mpc0;
0292 mpc.branch(2, BR_STATUS) = 0;
0293 r = runopf(mpc, mpopt);
0294 t_ok(r.success, [t 'success']);
0295 t_is(r.f, ff(3), 5, [t 'f']);
0296 t_is(r.gen(:, PG), [100; 100; 100; -300], 7, [t 'Pg']);
0297 t_is(r.bus(:, LAM_P), [20; 20; 1000], 7, [t 'lam P']);
0298 t_is(r.branch(:, MU_SF) + r.branch(:, MU_ST), [0; 0; 980], 7, [t 'mu flow']);
0299
0300
0301
0302 t = 'DC OPF (w/contingencies) : c3sopf ';
0303 mpc = mpc0;
0304 r = c3sopf(mpc, xgd_table.data, contab, mpopt);
0305
0306
0307
0308 t_ok(r.success, [t 'success']);
0309 t_is(r.opf_results.f, sum(pp.*ff), 5, [t 'f']);
0310 rr = r.base;
0311 t_is(rr.gen(:, PG), [135; 135; 180; -450], 7, [t 'Pg base']);
0312 t_is(rr.bus(:, LAM_P), [27; 36; 45]*pp(1), 7, [t 'lam P base']);
0313 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 27; 0]*pp(1), 7, [t 'mu flow base']);
0314 rr = r.cont(1);
0315 t_is(rr.gen(:, PG), [200; 0; 250; -450], 7, [t 'Pg 1']);
0316 t_is(rr.bus(:, LAM_P), [50; 50; 50]*pp(2), 7, [t 'lam P 1']);
0317 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0]*pp(2), 7, [t 'mu flow 1']);
0318 rr = r.cont(2);
0319 t_is(rr.gen(:, PG), [100; 100; 100; -300], 6, [t 'Pg 2']);
0320 t_is(rr.bus(:, LAM_P), [20; 20; 1000]*pp(3), 7, [t 'lam P 2']);
0321 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 980]*pp(3), 7, [t 'mu flow 2']);
0322
0323 t = 'DC OPF (w/contingencies) : most ';
0324 mdi = loadmd(mpc, [], xgd, [], contab);
0325 mdo = most(mdi, mpopt);
0326
0327
0328
0329 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0330 t_is(mdo.QP.f, sum(pp.*ff), 3, [t 'f']);
0331 rr = mdo.flow(1,1,1).mpc;
0332 t_is(rr.gen(:, PG), [135; 135; 180; -450], 5, [t 'Pg base']);
0333 t_is(rr.bus(:, LAM_P), [27; 36; 45]*pp(1), 5, [t 'lam P base']);
0334 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 27; 0]*pp(1), 5, [t 'mu flow base']);
0335 rr = mdo.flow(1,1,2).mpc;
0336 t_is(rr.gen(:, PG), [200; 0; 250; -450], 4, [t 'Pg 1']);
0337 t_is(rr.bus(:, LAM_P), [50; 50; 50]*pp(2), 5, [t 'lam P 1']);
0338 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0]*pp(2), 5, [t 'mu flow 1']);
0339 rr = mdo.flow(1,1,3).mpc;
0340 t_is(rr.gen(:, PG), [100; 100; 100; -300], 4, [t 'Pg 2']);
0341 t_is(rr.bus(:, LAM_P), [20; 20; 1000]*pp(3), 6, [t 'lam P 2']);
0342 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 980]*pp(3), 6, [t 'mu flow 2']);
0343
0344
0345
0346
0347 t = 'Secure DC OPF (w/cont,res,ramp) : c3sopf ';
0348 xgd_table = loadgenericdata('ex_xgd_res', 'struct');
0349 xgd = loadxgendata(xgd_table, mpc);
0350
0351
0352
0353
0354
0355 r = c3sopf(mpc, xgd_table.data, contab, mpopt);
0356
0357
0358
0359 t_ok(r.success, [t 'success']);
0360 t_is(r.opf_results.f, -436686.1245, 5, [t 'f']);
0361 rr = r.base;
0362
0363
0364
0365 t_is(rr.gen(:, PG), [142.15; 127.85; 180; -450], 7, [t 'Pg base']);
0366 t_is(rr.bus(:, LAM_P), [23.6958; 32.4; 41.1042], 7, [t 'lam P base']);
0367 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 26.1126; 0], 7, [t 'mu flow base']);
0368 rr = r.cont(1);
0369 t_is(rr.gen(:, PG), [142.15; 0; 307.85; -450], 6, [t 'Pg 1']);
0370 t_is(rr.bus(:, LAM_P), [5.1942; 5.1942; 5.1942], 7, [t 'lam P 1']);
0371 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 1']);
0372 rr = r.cont(2);
0373 t_is(rr.gen(:, PG), [142.15; 27.85; 130; -300], 5, [t 'Pg 2']);
0374 t_is(rr.bus(:, LAM_P), [-0.46; -0.46; 40], 7, [t 'lam P 2']);
0375 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 40.46], 7, [t 'mu flow 2']);
0376
0377 t = 'Secure DC OPF (w/cont,res,ramp) : most ';
0378 mdi = loadmd(mpc, [], xgd, [], contab);
0379 mdo = most(mdi, mpopt);
0380
0381
0382
0383 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0384 t_is(mdo.QP.f, -436686.1245, 2, [t 'f']);
0385 rr = mdo.flow(1,1,1).mpc;
0386 t_is(rr.gen(:, PG), [142.15; 127.85; 180; -450], 4, [t 'Pg base']);
0387 t_is(rr.bus(:, LAM_P), [23.6958; 32.4; 41.1042], 5, [t 'lam P base']);
0388 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 26.1126; 0], 5, [t 'mu flow base']);
0389 if create_plots
0390 lamP(:, j) = rr.bus(:, LAM_P);
0391 muF(:, j) = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0392 end
0393
0394 rr = mdo.flow(1,1,2).mpc;
0395 t_is(rr.gen(:, PG), [142.15; 0; 307.85; -450], 4, [t 'Pg 1']);
0396 t_is(rr.bus(:, LAM_P), [5.1942; 5.1942; 5.1942], 6, [t 'lam P 1']);
0397 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 1']);
0398 if create_plots
0399 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0400 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0401 end
0402
0403 rr = mdo.flow(1,1,3).mpc;
0404 t_is(rr.gen(:, PG), [142.15; 27.85; 130; -300], 3, [t 'Pg 2']);
0405 t_is(rr.bus(:, LAM_P), [-0.46; -0.46; 40], 6, [t 'lam P 2']);
0406 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 40.46], 6, [t 'mu flow 2']);
0407 if create_plots
0408 Pg(1:4, j) = mdo.results.ExpectedDispatch;
0409 Rp(1:4, j) = mdo.results.Rpp + mdo.results.Pc - mdo.results.ExpectedDispatch;
0410 Rm(1:4, j) = mdo.results.Rpm - mdo.results.Pc + mdo.results.ExpectedDispatch;
0411 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0412 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0413 j = j + 1;
0414 end
0415
0416 t = 'Stochastic DC OPF (w/wind,res) : c3sopf ';
0417 [iwind, mpc, xgd] = addwind(wind, mpc, xgd);
0418 nt = 1;
0419 nj = 3;
0420
0421 transmat = {[0.158655253931457; 0.682689492137086; 0.158655253931457]};
0422 profiles = getprofiles(uniformwindprofile(nt, nj), iwind);
0423
0424 pp = transmat{1};
0425
0426
0427 contab0 = [
0428 1 pp(1) profiles.table profiles.rows profiles.col profiles.chgtype profiles.values(1)/profiles.values(2);
0429 3 pp(3) profiles.table profiles.rows profiles.col profiles.chgtype profiles.values(3)/profiles.values(2);
0430 ];
0431 mpc0 = mpc;
0432 mpc0.gen(iwind, PMAX) = mpc0.gen(iwind, PMAX) * profiles.values(2);
0433 offers = [xgd_table.data; wind.xgd_table.data(:, [3:8])];
0434
0435
0436 r = c3sopf(mpc0, offers, contab0, mpopt);
0437 t_ok(r.success, [t 'success']);
0438 t_is(r.opf_results.f, -444525.605052, 5, [t 'f']);
0439 rr = r.cont(1);
0440
0441
0442
0443 t_is(rr.gen(:, PG), [128.7823267; 141.2176733; 180; -450; 0], 7, [t 'Pg 1']);
0444 t_is(rr.bus(:, LAM_P), [4.4809852; 7.2115891; 9.9421931], 7, [t 'lam P 1']);
0445 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 8.1918119; 0], 7, [t 'mu flow 1']);
0446 rr = r.base;
0447 t_is(rr.gen(:, PG), [128.7823267; 135.6088362; 135.6088370; -450; 50], 5, [t 'Pg 2']);
0448 t_is(rr.bus(:, LAM_P), [18.5157455; 18.5157455; 18.5157455], 6, [t 'lam P 2']);
0449 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2']);
0450 rr = r.cont(2);
0451 t_is(rr.gen(:, PG), [128.7823267; 86.9726845; 134.2449888; -450; 100], 5, [t 'Pg 3']);
0452 t_is(rr.bus(:, LAM_P), [2.7597346; 2.7597346; 2.7597346], 6, [t 'lam P 3']);
0453 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 3']);
0454
0455 t = 'Stochastic DC OPF (w/wind,res) : most ';
0456 mdi = loadmd(mpc, transmat, xgd, [], [], profiles);
0457 mdo = most(mdi, mpopt);
0458
0459
0460
0461 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0462 t_is(mdo.QP.f, -444525.605052, 5, [t 'f']);
0463 rr = mdo.flow(1,1,1).mpc;
0464
0465
0466
0467 t_is(rr.gen(:, PG), [128.7823267; 141.2176733; 180; -450; 0], 7, [t 'Pg 1']);
0468 t_is(rr.bus(:, LAM_P), [4.4809852; 7.2115891; 9.9421931], 7, [t 'lam P 1']);
0469 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 8.1918119; 0], 7, [t 'mu flow 1']);
0470 if create_plots
0471 lamP(:, j) = rr.bus(:, LAM_P);
0472 muF(:, j) = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0473 end
0474
0475 rr = mdo.flow(1,2,1).mpc;
0476
0477
0478
0479 t_is(rr.gen(:, PG), [128.7823267; 135.6088362; 135.6088370; -450; 50], 5, [t 'Pg 2']);
0480 t_is(rr.bus(:, LAM_P), [18.5157455; 18.5157455; 18.5157455], 6, [t 'lam P 2']);
0481 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2']);
0482 if create_plots
0483 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0484 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0485 end
0486
0487 rr = mdo.flow(1,3,1).mpc;
0488
0489
0490
0491 t_is(rr.gen(:, PG), [128.7823267; 86.9726845; 134.2449888; -450; 100], 5, [t 'Pg 3']);
0492 t_is(rr.bus(:, LAM_P), [2.7597346; 2.7597346; 2.7597346], 6, [t 'lam P 3']);
0493 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 7, [t 'mu flow 3']);
0494 if create_plots
0495 Pg(:, j) = mdo.results.ExpectedDispatch;
0496 Rp(:, j) = mdo.results.Rpp + mdo.results.Pc - mdo.results.ExpectedDispatch;
0497 Rm(:, j) = mdo.results.Rpm - mdo.results.Pc + mdo.results.ExpectedDispatch;
0498 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0499 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0500 j = j + 1;
0501 end
0502
0503
0504
0505 t = 'Secure Stochastic DC OPF (w/wind,cont,res,ramp) : most ';
0506 mdi = loadmd(mpc, transmat, xgd, [], contab, profiles);
0507 mdo = most(mdi, mpopt);
0508
0509
0510
0511 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0512 t_is(mdo.QP.f, -438182.429, 5, [t 'f']);
0513
0514 rr = mdo.flow(1,1,1).mpc;
0515
0516
0517
0518 t_is(rr.gen(:, PG), [137.6930222; 132.3069578; 180; -450; 0], 4, [t 'Pg 1 base']);
0519 t_is(rr.bus(:, LAM_P), [3.9958603; 5.1404304; 6.2850005], 3, [t 'lam P 1 base']);
0520 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 3.4337103; 0], 2, [t 'mu flow 1 base']);
0521 if create_plots
0522 lamP(:, j) = rr.bus(:, LAM_P);
0523 muF(:, j) = rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0524 end
0525
0526 rr = mdo.flow(1,1,2).mpc;
0527
0528
0529
0530 t_is(rr.gen(:, PG), [137.6930928; 0; 312.3069041; -450; 0], 5, [t 'Pg 1 1']);
0531 t_is(rr.bus(:, LAM_P), [2.0945890; 2.0945890; 2.0945892], 6, [t 'lam P 1 1']);
0532 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 1 1']);
0533 if create_plots
0534 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0535 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0536 end
0537
0538 rr = mdo.flow(1,1,3).mpc;
0539
0540
0541
0542 t_is(rr.gen(:, PG), [137.6930139; 45.2766586; 117.0303189; -300; 0], 4, [t 'Pg 1 2']);
0543 t_is(rr.bus(:, LAM_P), [0.0574664; 0.0574664; 6.3462102], 6, [t 'lam P 1 2']);
0544 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 6.2887437], 6, [t 'mu flow 1 2']);
0545 if create_plots
0546 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0547 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0548 end
0549
0550 rr = mdo.flow(1,2,1).mpc;
0551
0552
0553
0554 t_is(rr.gen(:, PG), [137.6929448; 131.1533869; 131.1535795; -450; 50], 4, [t 'Pg 2 base']);
0555 t_is(rr.bus(:, LAM_P), [16.1166803; 16.1166890; 16.1166978], 6, [t 'lam P 2 base']);
0556 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 4, [t 'mu flow 2 base']);
0557 if create_plots
0558 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0559 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0560 end
0561
0562 rr = mdo.flow(1,2,2).mpc;
0563
0564
0565
0566 t_is(rr.gen(:, PG), [137.6930711; 0; 262.3068485; -450; 50], 4, [t 'Pg 2 1']);
0567 t_is(rr.bus(:, LAM_P), [2.1488905; 2.1488905; 2.1488907], 6, [t 'lam P 2 1']);
0568 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 2 1']);
0569 if create_plots
0570 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0571 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0572 end
0573
0574 rr = mdo.flow(1,2,3).mpc;
0575
0576
0577
0578 t_is(rr.gen(:, PG), [137.6929611; 32.3066174; 117.0299763; -300; 12.9704443], 6, [t 'Pg 2 2']);
0579 t_is(rr.bus(:, LAM_P), [0; 0; 27.3075797], 6, [t 'lam P 2 2']);
0580 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 27.3075797], 6, [t 'mu flow 2 2']);
0581 if create_plots
0582 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0583 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0584 end
0585
0586 rr = mdo.flow(1,3,1).mpc;
0587
0588
0589
0590 t_is(rr.gen(:, PG), [137.6929341; 95.2771176; 117.0299561; -450; 100], 4, [t 'Pg 3 base']);
0591 t_is(rr.bus(:, LAM_P), [2.7209142; 2.7209144; 2.7209147], 6, [t 'lam P 3 base']);
0592 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 5, [t 'mu flow 3 base']);
0593 if create_plots
0594 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0595 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0596 end
0597
0598 rr = mdo.flow(1,3,2).mpc;
0599
0600
0601
0602 t_is(rr.gen(:, PG), [137.6930070; 0.0000000; 212.3070512; -450; 100], 4, [t 'Pg 3 1']);
0603 t_is(rr.bus(:, LAM_P), [0.4042034; 0.4042034; 0.4042036], 6, [t 'lam P 3 1']);
0604 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 0], 6, [t 'mu flow 3 1']);
0605 if create_plots
0606 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0607 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0608 end
0609
0610 rr = mdo.flow(1,3,3).mpc;
0611
0612
0613
0614 t_is(rr.gen(:, PG), [137.6930131; 32.3071796; 117.0302086; -300; 12.9695914], 5, [t 'Pg 3 2']);
0615 t_is(rr.bus(:, LAM_P), [0; 0; 6.3462102], 6, [t 'lam P 3 2']);
0616 t_is(rr.branch(:, MU_SF) + rr.branch(:, MU_ST), [0; 0; 6.3462102], 6, [t 'mu flow 3 2']);
0617
0618
0619 if create_plots
0620 Pg(:, j) = mdo.results.ExpectedDispatch;
0621 Rp(:, j) = mdo.results.Rpp + mdo.results.Pc - mdo.results.ExpectedDispatch;
0622 Rm(:, j) = mdo.results.Rpm - mdo.results.Pc + mdo.results.ExpectedDispatch;
0623 lamP(:, j) = lamP(:, j) + rr.bus(:, LAM_P);
0624 muF(:, j) = muF(:, j) + rr.branch(:, MU_SF) + rr.branch(:, MU_ST);
0625
0626 R = Rp + Rm;
0627
0628 labels = {'Economic Dispatch'; 'DC OPF'; 'Economic Dispatch w/reserves'; 'DC OPF w/reserves'; 'secure DC OPF'; 'stochastic DC OPF'; 'secure stochastic DC OPF'};
0629
0630 figure(1);
0631
0632 subplot(4, 1, 1);
0633 bar(abs(Pg([1:3 5],:)'),'stacked');
0634 title('Generation');
0635 ylabel('MW');
0636 h = gca;
0637 h.XTickLabel = labels;
0638 h.XTickLabelRotation = 20;
0639
0640 subplot(4, 1, 2);
0641 bar(abs(R([1:3 5],:)'),'stacked');
0642 title('Reserves');
0643 ylabel('MW');
0644 h = gca;
0645 h.XTickLabel = labels;
0646 h.XTickLabelRotation = 20;
0647
0648 subplot(4, 1, 3);
0649 plot(lamP');
0650 title('Nodal Prices');
0651 ylabel('$/MW');
0652 h = gca;
0653 h.XTickLabel = {'', labels{:}, ''}';
0654 h.XTickLabelRotation = 20;
0655 h = [0 8 0 100];
0656 axis(h);
0657
0658 subplot(4, 1, 4);
0659 plot(muF');
0660 title('Flow Constraint Shadow Prices');
0661 ylabel('$/MW');
0662 h = gca;
0663 h.XTickLabel = {'', labels{:}, ''}';
0664 h.XTickLabelRotation = 20;
0665 h = [0 8 0 60];
0666 axis(h);
0667
0668 if create_pdfs
0669 fname = 'single-period-continuous';
0670 h = gcf;
0671 set(h, 'PaperSize', [11 8.5]);
0672 set(h, 'PaperPosition', [0.25 0.25 10.5 8]);
0673 print('-dpdf', fullfile(savedir, fname));
0674 end
0675
0676 for j = 1:7;
0677 figure(j+1);
0678 if create_pdfs
0679 fname = sprintf('single-period-cont-%d', j);
0680 else
0681 fname = '';
0682 end
0683 plot_case(labels{j}, Pg([1:3 5], j), Rp([1:3 5], j), Rm([1:3 5], j), lamP(:, j), muF(:, j), 320, 100, savedir, fname);
0684 end
0685 end
0686
0687
0688 if have_feature('octave')
0689 warning(s.state, 'Octave:nearly-singular-matrix');
0690 end
0691
0692 t_end;
0693
0694
0695 function h = plot_case(label, Pg, Rp, Rm, lamP, muF, maxq, maxp, mypath, fname)
0696
0697 subplot(1, 3, 1);
0698 h = bar([Pg-Rm Rm Rp], 'stacked');
0699 set(h(2), 'FaceColor', [0 0.35 0.33]);
0700 ah1 = gca;
0701 title('Generation & Reserves');
0702 ylabel('MW');
0703 xlabel('Gen');
0704 set(gca, 'FontSize', 14);
0705
0706 if nargin < 6
0707 maxq = ah1.YLim(2);
0708 end
0709 ah1.YLim(2) = maxq;
0710 ah1.YLim(1) = 0;
0711
0712 subplot(1, 3, 2);
0713 bar(lamP);
0714 ah3 = gca;
0715 title('Nodal Prices');
0716 ylabel('$/MW');
0717 xlabel('Bus');
0718 set(gca, 'FontSize', 14);
0719
0720 subplot(1, 3, 3);
0721 bar(muF);
0722 ah4 = gca;
0723 title('Flow Constraint Shadow Prices');
0724 ylabel('$/MW');
0725 xlabel('Line');
0726 set(gca, 'FontSize', 14);
0727
0728 if nargin < 7
0729 maxp = max(ah3.YLim(2), ah4.YLim(2));
0730 end
0731 ah3.YLim(1) = 0;
0732 ah4.YLim(1) = 0;
0733 ah3.YLim(2) = maxp;
0734 ah4.YLim(2) = maxp;
0735
0736 [ax,h] = suplabel(label, 't');
0737 set(h, 'FontSize', 18)
0738
0739 if nargin > 7 && ~isempty(fname)
0740 h = gcf;
0741 set(h, 'PaperSize', [11 8.5]);
0742 set(h, 'PaperPosition', [0.25 0.25 10.5 8]);
0743 print('-dpdf', fullfile(mypath, fname));
0744 end