Home > matpower7.1 > most > lib > t > t_most_spuc.m

t_most_spuc

PURPOSE ^

T_MOST_SPUC Tests of single-period unit commitment optimizations

SYNOPSIS ^

function t_most_spuc(quiet, create_plots, create_pdfs, savedir)

DESCRIPTION ^

T_MOST_SPUC  Tests of single-period unit commitment optimizations

   T_MOST_SPUC(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_spuc(0, 1, 1, '~/Downloads/spuc_plots')

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function t_most_spuc(quiet, create_plots, create_pdfs, savedir)
0002 %T_MOST_SPUC  Tests of single-period unit commitment optimizations
0003 %
0004 %   T_MOST_SPUC(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_spuc(0, 1, 1, '~/Downloads/spuc_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 
0030 solvers = {'CPLEX', 'GLPK', 'GUROBI', 'MOSEK', 'OT'};
0031 fcn = {'cplex', 'glpk', 'gurobi', 'mosek', 'intlinprog'};
0032 % solvers = {'OT'};
0033 % fcn = {'intlinprog'};
0034 % solvers = {'GUROBI'};
0035 % fcn = {'gurobi'};
0036 % solvers = {'GLPK'};
0037 % fcn = {'glpk'};
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 % verbose = 2;
0047 
0048 casefile = 'ex_case3b';
0049 mpopt = mpoption;
0050 mpopt = mpoption(mpopt, 'out.gen', 1);
0051 mpopt = mpoption(mpopt, 'verbose', verbose);
0052 % mpopt = mpoption(mpopt, 'opf.violation', 1e-6, 'mips.gradtol', 1e-8, ...
0053 %         'mips.comptol', 1e-8, 'mips.costtol', 1e-8);
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 %% solver options
0059 if have_feature('cplex')
0060     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 0);       %% automatic
0061     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 1);       %% primal simplex
0062     mpopt = mpoption(mpopt, 'cplex.lpmethod', 2);       %% dual simplex
0063     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 3);       %% network simplex
0064     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 4);       %% barrier
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     %mpopt = mpoption(mpopt, 'gurobi.method', -1);       %% automatic
0076     %mpopt = mpoption(mpopt, 'gurobi.method', 0);        %% primal simplex
0077     mpopt = mpoption(mpopt, 'gurobi.method', 1);        %% dual simplex
0078     %mpopt = mpoption(mpopt, 'gurobi.method', 2);        %% barrier
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     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_FREE);            %% default
0086     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_INTPNT);          %% interior point
0087     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_PRIMAL_SIMPLEX);  %% primal simplex
0088     mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);     %% dual simplex
0089     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_FREE_SIMPLEX);    %% automatic simplex
0090     %mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_X', 0);
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     %mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_RELAX_INT', 0);
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     %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'interior-point');
0100     %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'active-set');
0101     %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'simplex');
0102     mpopt = mpoption(mpopt, 'linprog.Algorithm', 'dual-simplex');
0103     %mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'primal-simplex');
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     %% next line is to work around a bug in intlinprog
0110     % (Technical Support Case #01841662)
0111     mpopt = mpoption(mpopt, 'intlinprog.LPPreprocess', 'none');
0112 end
0113 if ~verbose
0114     mpopt = mpoption(mpopt, 'out.all', 0);
0115 end
0116 % mpopt = mpoption(mpopt, 'out.all', -1);
0117 
0118 %% define named indices into data matrices
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 %% load base case file
0136 mpc = loadcase(casefile);
0137 mpc.gencost(:, STARTUP) = 0;
0138 mpc.gencost(:, SHUTDOWN) = 0;
0139 
0140 %%-----  contingencies  -----
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 %% data structures for results for plotting
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})    %% check if we have the solver
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 %%-----  economic dispatch (no network)  -----
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 %% most
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 % xgd = xg
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 % rr.gen(:, GEN_STATUS)
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 %%-----  DC OPF  -----
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 % mpc = scale_load(500, mpc, [], struct('scale', 'QUANTITY'));
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 %% most
0229 t = sprintf('%s : DC OPF : most   ', solvers{s});
0230 mpc = mpc0;
0231 % mpc = scale_load(500, mpc, [], struct('scale', 'QUANTITY'));
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 % rr.gen(:, GEN_STATUS)
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 %%-----  economic dispatch (w/reserves)  -----
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 %% most
0278 t = sprintf('%s : economic dispatch (w/reserves) : most   ', solvers{s});
0279 mpc = mpc0;
0280 % mpc = scale_load(350.8, mpc, [], struct('scale', 'QUANTITY'));
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 % rr.gen(:, GEN_STATUS)
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 %%-----  DC OPF  -----
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 %% most
0328 t = sprintf('%s : DC OPF (w/reserves) : most   ', solvers{s});
0329 mpc = mpc0;
0330 % mpc = scale_load(350.8, mpc, [], struct('scale', 'QUANTITY'));
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 % rr.gen(:, GEN_STATUS)
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 % xgd = loadxgendata(xgd_table, mpc);
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 % r
0384 % r.opf_results
0385 % r.opf_results.f
0386 t_ok(r.success, [t 'success']);
0387 t_is(r.opf_results.f, -340850, 4, [t 'f']);
0388 rr = r.base;
0389 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0390 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0391 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0398 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0399 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0405 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0406 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % mdo
0416 % mdo.QP
0417 % mdo.QP.f
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0422 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0423 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0424 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % t_is(rr.bus(:, LAM_P), [20; 20; 20], 7, [t 'lam P base']);
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0437 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0438 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0439 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0451 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0452 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0453 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % transmat = transmatnormalpersistent(nt, nj);
0471 transmat = {[0.158655253931457; 0.682689492137086; 0.158655253931457]};
0472 profiles = getprofiles(uniformwindprofile(nt, nj), iwind);
0473 
0474 pp = transmat{1};
0475 %% contingency table
0476 % label probty  type            row             column          chgtype newvalue
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 % mpopt.verbose = 2;
0487 % mpopt.out.all = -1;
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 % rr.gen(:, PG)
0493 % rr.bus(:, LAM_P)
0494 % rr.branch(:, MU_SF) + rr.branch(:, MU_ST)
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 % mdo
0511 % mdo.QP
0512 % mdo.QP.f
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0517 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0518 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0519 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0531 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0532 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0533 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0545 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0546 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0547 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % keyboard;
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 % mdo
0566 % mdo.QP
0567 % mdo.QP.f
0568 t_ok(mdo.QP.exitflag > 0, [t 'success']);
0569 t_is(mdo.QP.f, -338857.9224447, 4, [t 'f']);
0570 % fprintf('%.7f\n', mdo.QP.f);
0571 
0572 rr = mdo.flow(1,1,1).mpc;
0573 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0574 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0575 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0576 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % t_is(rr.bus(:, LAM_P), [7.2115891; 7.2115891; 7.2115891], 6, [t 'lam P 1 base']);
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0589 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0590 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0591 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % t_is(rr.bus(:, LAM_P), [0.3807726; 0.3807726; 0.3807726], 7, [t 'lam P 1 1']);
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0604 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0605 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0606 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0618 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0619 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0620 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0632 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0633 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0634 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0646 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0647 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0648 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0660 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0661 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0662 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0674 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0675 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0676 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0688 % fprintf('[%d; %d; %d; %d; %d]\n', rr.gen(:, GEN_STATUS));
0689 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0690 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % keyboard;
0696 end
0697 
0698 if create_plots
0699     create_plots = 0;   %% don't do them again
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 %    R(4:5, :) = NaN;
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

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005