most_ex6_uc

MOST_EX6_UC  Examples of deterministic unit commitment problems.

0001 function most_ex6_uc(quiet)
0002 %MOST_EX6_UC  Examples of deterministic unit commitment problems.
0003 
0004 %   MOST
0005 %   Copyright (c) 2015-2016, Power Systems Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   This file is part of MOST.
0009 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0010 %   See https://github.com/MATPOWER/most for more info.
0011 
0012 %% set up options
0013 define_constants;
0014 verbose = 1;
0015 mpopt = mpoption('verbose', verbose);
0016 mpopt = mpoption(mpopt, 'out.gen', 1);
0017 mpopt = mpoption(mpopt, 'model', 'DC');
0018 mpopt = mpoption(mpopt, 'most.solver', 'DEFAULT');
0019 if ~verbose
0020     mpopt = mpoption(mpopt, 'out.all', 0);
0021 end
0022 mpopt = mpoption(mpopt, 'most.price_stage_warn_tol', 1e-5);
0023 
0024 %% solver options
0025 if have_fcn('cplex')
0026     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 0);       %% automatic
0027     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 1);       %% primal simplex
0028     mpopt = mpoption(mpopt, 'cplex.lpmethod', 2);       %% dual simplex
0029     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 3);       %% network simplex
0030     %mpopt = mpoption(mpopt, 'cplex.lpmethod', 4);       %% barrier
0031     mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.mipgap', 0);
0032     mpopt = mpoption(mpopt, 'cplex.opts.mip.tolerances.absmipgap', 0);
0033     mpopt = mpoption(mpopt, 'cplex.opts.threads', 2);
0034 end
0035 if have_fcn('glpk')
0036     mpopt = mpoption(mpopt, 'glpk.opts.mipgap', 0);
0037     mpopt = mpoption(mpopt, 'glpk.opts.tolint', 1e-10);
0038     mpopt = mpoption(mpopt, 'glpk.opts.tolobj', 1e-10);
0039 end
0040 if have_fcn('gurobi')
0041     %mpopt = mpoption(mpopt, 'gurobi.method', -1);       %% automatic
0042     %mpopt = mpoption(mpopt, 'gurobi.method', 0);        %% primal simplex
0043     mpopt = mpoption(mpopt, 'gurobi.method', 1);        %% dual simplex
0044     %mpopt = mpoption(mpopt, 'gurobi.method', 2);        %% barrier
0045     mpopt = mpoption(mpopt, 'gurobi.threads', 2);
0046     mpopt = mpoption(mpopt, 'gurobi.opts.MIPGap', 0);
0047     mpopt = mpoption(mpopt, 'gurobi.opts.MIPGapAbs', 0);
0048 end
0049 if have_fcn('mosek')
0050     sc = mosek_symbcon;
0051     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_FREE);            %% default
0052     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_INTPNT);          %% interior point
0053     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_PRIMAL_SIMPLEX);  %% primal simplex
0054     mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);     %% dual simplex
0055     %mpopt = mpoption(mpopt, 'mosek.lp_alg', sc.MSK_OPTIMIZER_FREE_SIMPLEX);    %% automatic simplex
0056     %mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_X', 0);
0057     mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_NODE_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
0058     mpopt = mpoption(mpopt, 'mosek.opts.MSK_IPAR_MIO_ROOT_OPTIMIZER', sc.MSK_OPTIMIZER_DUAL_SIMPLEX);
0059     mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_RELAX_INT', 1e-9);
0060     %mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_RELAX_INT', 0);
0061     mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_REL_GAP', 0);
0062     mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_MIO_TOL_ABS_GAP', 0);
0063 end
0064 if have_fcn('intlinprog')
0065     %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'interior-point');
0066     %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'active-set');
0067     %mpopt = mpoption(mpopt, 'linprog.Algorithm', 'simplex');
0068     mpopt = mpoption(mpopt, 'linprog.Algorithm', 'dual-simplex');
0069     %mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'primal-simplex');
0070     mpopt = mpoption(mpopt, 'intlinprog.RootLPAlgorithm', 'dual-simplex');
0071     mpopt = mpoption(mpopt, 'intlinprog.TolCon', 1e-9);
0072     mpopt = mpoption(mpopt, 'intlinprog.TolGapAbs', 0);
0073     mpopt = mpoption(mpopt, 'intlinprog.TolGapRel', 0);
0074     mpopt = mpoption(mpopt, 'intlinprog.TolInteger', 1e-6);
0075     %% next line is to work around a bug in intlinprog
0076     % (Technical Support Case #01841662)
0077     % (except actually in this case it triggers it rather than working
0078     %  around it, so we comment it out)
0079     %mpopt = mpoption(mpopt, 'intlinprog.LPPreprocess', 'none');
0080 end
0081 
0082 casefile = 'ex_case3b';
0083 mpc = loadcase(casefile);
0084 xgd = loadxgendata('ex_xgd_uc', mpc);
0085 [iwind, mpc, xgd] = addwind('ex_wind_uc', mpc, xgd);
0086 profiles = getprofiles('ex_wind_profile_d', iwind);
0087 profiles = getprofiles('ex_load_profile', profiles);
0088 nt = size(profiles(1).values, 1);       % number of periods
0089 mpc_full = mpc;
0090 xgd_full = xgd;
0091 mpc.gencost(:, [STARTUP SHUTDOWN]) = 0; % remove startup/shutdown costs
0092 xgd.MinUp(2) = 1;                       % remove min up-time constraint
0093 xgd.PositiveLoadFollowReserveQuantity(3) = 250; % remove ramp reserve
0094 xgd.PositiveLoadFollowReservePrice(3) = 1e-6;   % constraint and costs
0095 xgd.NegativeLoadFollowReservePrice(3) = 1e-6;
0096 
0097 
0098 %%-----  Base : No Network  -----
0099 mpopt = mpoption(mpopt, 'most.dc_model', 0);    % use model with no network
0100 mdi = loadmd(mpc, nt, xgd, [], [], profiles);
0101 mdo = most(mdi, mpopt);
0102 if verbose
0103     ms = most_summary(mdo);
0104 end
0105 
0106 %%-----  + DC Network  -----
0107 mpopt = mpoption(mpopt, 'most.dc_model', 1);    % use DC network model (default)
0108 mdo = most(mdi, mpopt);
0109 if verbose
0110     ms = most_summary(mdo);
0111 end
0112 
0113 %%-----  + startup/shutdown costs  -----
0114 mpc.gencost(2, [STARTUP SHUTDOWN]) = [ 200 200];
0115 mpc.gencost(3, [STARTUP SHUTDOWN]) = [3000 600];
0116 % equivalent to doing: mpc = mpc_full;
0117 mdi = loadmd(mpc, nt, xgd, [], [], profiles);
0118 mdo = most(mdi, mpopt);
0119 if verbose
0120     ms = most_summary(mdo);
0121 end
0122 
0123 %%-----  + min up/down time constraints  -----
0124 xgd.MinUp(2) = 3;
0125 mdi = loadmd(mpc, nt, xgd, [], [], profiles);
0126 mdo = most(mdi, mpopt);
0127 if verbose
0128     ms = most_summary(mdo);
0129 end
0130 
0131 %%-----  + ramp constraint/ramp res cost  -----
0132 xgd.PositiveLoadFollowReserveQuantity(3) = 100; % restore ramp reserve
0133 xgd.PositiveLoadFollowReservePrice(3) = 10;     % constraint and costs
0134 xgd.NegativeLoadFollowReservePrice(3) = 10;
0135 % equivalent to doing: xgd = xgd_full;
0136 mdi = loadmd(mpc, nt, xgd, [], [], profiles);
0137 mdo = most(mdi, mpopt);
0138 if verbose
0139     ms = most_summary(mdo);
0140 end
0141 
0142 %%-----  + storage  -----
0143 mpopt = mpoption(mpopt, 'most.storage.cyclic', 1);
0144 [iess, mpc, xgd, sd] = addstorage('ex_storage', mpc, xgd);
0145 mdi = loadmd(mpc, nt, xgd, sd, [], profiles);
0146 mdo = most(mdi, mpopt);
0147 if verbose
0148     ms = most_summary(mdo);
0149 end