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

t_most_sp

PURPOSE ^

T_MOST_SP Tests of single-period continuous optimizations

SYNOPSIS ^

function t_most_sp(quiet, create_plots, create_pdfs, savedir)

DESCRIPTION ^

T_MOST_SP  Tests of single-period continuous optimizations

   T_MOST_SP(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_sp(0, 1, 1, '~/Downloads/sp_plots')

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function t_most_sp(quiet, create_plots, create_pdfs, savedir)
0002 %T_MOST_SP  Tests of single-period continuous optimizations
0003 %
0004 %   T_MOST_SP(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_sp(0, 1, 1, '~/Downloads/sp_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 t_begin(156, quiet);
0031 
0032 if quiet
0033     verbose = 0;
0034 else
0035     verbose = 0;
0036 end
0037 % verbose = 2;
0038 
0039 casefile = 'ex_case3a';
0040 mpopt = mpoption;
0041 mpopt = mpoption(mpopt, 'out.gen', 1);
0042 mpopt = mpoption(mpopt, 'verbose', verbose);
0043 % mpopt = mpoption(mpopt, 'opf.violation', 1e-6, 'mips.gradtol', 1e-8, ...
0044 %         'mips.comptol', 1e-8, 'mips.costtol', 1e-8);
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 % mpopt = mpoption(mpopt, 'out.all', -1);
0053 
0054 %% turn off warnings
0055 if have_feature('octave')
0056     s = warning('query', 'Octave:nearly-singular-matrix');
0057     warning('off', 'Octave:nearly-singular-matrix');
0058 end
0059 
0060 %% define named indices into data matrices
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 %% load base case file
0077 mpc = loadcase(casefile);
0078 
0079 %% initial xGenData
0080 xgd_table = loadgenericdata('ex_xgd', 'struct');
0081 
0082 %%-----  wind  -----
0083 wind = loadgenericdata('ex_wind', 'struct');
0084 
0085 %%-----  contingencies  -----
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 %% data structures for results for plotting
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 %%-----  economic dispatch (no network)  -----
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 %% most
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 %%-----  DC OPF  -----
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 %% most
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 %%-----  economic dispatch (w/reserves)  -----
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 %% most
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 %%-----  DC OPF  -----
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 %% most
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 %%-----  DC OPF w/contingencies  -----
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 % mpopt = mpoption(mpopt, 'mips.comptol', 1e-8);
0301 
0302 t = 'DC OPF (w/contingencies) : c3sopf ';
0303 mpc = mpc0;
0304 r = c3sopf(mpc, xgd_table.data, contab, mpopt);
0305 % r
0306 % r.opf_results
0307 % r.opf_results.f
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 % mdo
0327 % mdo.QP
0328 % mdo.QP.f
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 % mpopt.verbose = 2;
0345 % mpopt.out.all = -1;
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 % xgd.PositiveActiveReserveQuantity(2) = 100;
0351 % xgd.PositiveActiveReserveQuantity(4) = 500;
0352 % xgd.NegativeActiveReserveQuantity(:) = 0;
0353 % xgd.PositiveActiveReservePrice([1;3]) = [5; 1.5];
0354 % xgd.NegativeActiveReservePrice([1;3]) = [10; 3];
0355 r = c3sopf(mpc, xgd_table.data, contab, mpopt);
0356 % r
0357 % r.opf_results
0358 % r.opf_results.f
0359 t_ok(r.success, [t 'success']);
0360 t_is(r.opf_results.f, -436686.1245, 5, [t 'f']);
0361 rr = r.base;
0362 % rr.gen(:, PG)
0363 % rr.bus(:, LAM_P)
0364 % rr.branch(:, MU_SF) + rr.branch(:, MU_ST)
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 % mdo
0381 % mdo.QP
0382 % mdo.QP.f
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 % transmat = transmatnormalpersistent(nt, nj);
0421 transmat = {[0.158655253931457; 0.682689492137086; 0.158655253931457]};
0422 profiles = getprofiles(uniformwindprofile(nt, nj), iwind);
0423 
0424 pp = transmat{1};
0425 %% contingency table
0426 % label probty  type            row             column          chgtype newvalue
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 % mpopt.verbose = 2;
0435 % mpopt.out.all = -1;
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0441 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0442 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % mdo
0459 % mdo.QP
0460 % mdo.QP.f
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0465 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0466 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0477 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0478 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0489 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0490 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % keyboard;
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 % mdo
0509 % mdo.QP
0510 % mdo.QP.f
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0516 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0517 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0528 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0529 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0540 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0541 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0552 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0553 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0564 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0565 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0576 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0577 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0588 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0589 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0600 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0601 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % fprintf('[%.7f; %.7f; %.7f; %.7f; %.7f]\n', rr.gen(:, PG));
0612 % fprintf('[%.7f; %.7f; %.7f]\n', rr.bus(:, LAM_P));
0613 % fprintf('[%.7f; %.7f; %.7f]\n', rr.branch(:, MU_SF) + rr.branch(:, MU_ST));
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 % keyboard;
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 %    R(4:5, :) = NaN;
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 %% turn warnings back on
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

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