t_most_uc

T_MOST_UC  Tests of deteministic unit commitment optimizations

   T_MOST_UC(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_uc(0, 1, 1, '~/Downloads/uc_plots')

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